diagm and spdiagm

[fredrikekre]

Current state of affairs

There is currently a lot of inconsistencies between diagm and spdiagm. Some methods of diagm and its behavior:

• diagm(v::AbstractVector, k): place vector v on the kth diagonal
• diagm(v::SparseMatrixCSC): Note SparseMatrixCSC, not SparseVector, and it has to be of size (n,1) or (1,n). This places places the entries on the 0th diagonal (no choice).
• diagm(v::BitMatrix): Same thing here, matrix of size (n,1) or (1,n) and no choice of placement.

spdiagm is a bit more flexible: we can supply multiple vectors for multiple diagonals, e.g. spdiagm((v1, v2), (k1,k2)) places vector v1 and v2 on diagonal k1 and k2 respectively. We can also specify the resulting matrix size, e.g. spdiagm((v1, v2), (k1,k2), m, n) which will create a m×n matrix. Is this actually useful? Looking at the code, it seems like its there to kind of mirror the sparse(I,J,V) / sparse(I,J,V,m,n) methods.

Action plan

The most obvious things are:

• Deprecate diagm(A::SparseMatrixCSC) in favor of spdiagm(sparsevec(A)) (#23341)
• Deprecate diagm(A::BitMatrix) in favor of diagm(vec(A)) (#23373)

Change to nicer Pair API instead of the diagm(tuple_of_vecs, tuple_of_diags) methods.

• Deprecate spdiagm(tuple_of_vecs, tuple_of_diags) in favor of spdiagm(p::Pair{<:Integer,<:AbstractVector}...) (#23757)
• Implement diagm(p::Pair{<:Integer,<:AbstractVector}...) to match the spdiagm method. (#24047)

Then some other things:

• Remove the option to supply resulting matrix size for spdiagm and always return a square matrix like diagm. (#23757)

So when we are here, we have diagm and spdiagm that match each other. Where could we go from here? We delete spdiagm! 🎉

• Deprecate spdiagm(vecs, diags) in favor of diagm(sparsevec.(vecs), diags)

So now we only have diagm. Lets delete that too?

• Deprecate diagm(v::AbstractVector) in favor of Matrix(Diagonal(v))
• Deprecate diagm(v::SparseVector) in favor of sparse(Diagonal(v))
• Deprecate diagm(v::BitVector) in favor of BitMatrix(Diagonal(v))

The only problem with this is that we can't really specify diagonals easily anymore, and lose the option to supply more than one vector. So perhaps this last step was just too much.

[stevengj]

Requiring diagm(sparsevec.(vecs), diags) seems odd. In most uses of spdiagm, I suspect that the diagonals themselves are not sparse, so it is weird to require a SparseVector copy.

[fredrikekre]

Yea, I know and agree that it is suboptimal, but it would be nice to get rid of the sp* functions... Perhaps diagm(SparseMatrixCSC, vecs, diags) would be an alternative spelling?

[tkoolen]

If this is being redesigned, could JuliaArrays/OffsetArrays.jl#25 be taken into account?

[StefanKarpinski]

Would diagm(k₁ => d₁, k₂ => d₂, k₃ => d₃) be a nicer API?

[tkoolen]

I'm in favor of those last three deprecations in the issue description, as that would allow us to not have to think about what the indices for diagm(::OffsetVector, ::Int) should be, which is unclear (see OffsetArrays issue).

OffsetArrays.jl could then choose to have diagm(::OffsetVector, ::Int) (or diagm(::Pair{Int, <:AbstractVector}...)) error, while still having the most common case (constructing a truly diagonal matrix) work by overloading Diagonal(::OffsetVector) and indices(::Diagonal{<:OffsetVector}) in OffsetArrays. Or perhaps better, change the Diagonal code in Base so that Diagonal{<:OffsetVector} automagically has the right behavior.

cc: @timholy

[stevengj]

@fredrikekre, spdiagm is a lot easier to type, we have a lot of sp* functions, and some variant of spdiag is fairly common (Matlab, Numpy, R). diagm(SparseMatrixCSC, vecs, diags) doesn't seem like an improvement.

[ViralBShah]

It is unsatisfying to have the sp* functions, but passing SparseMatrixCSC as an input argument is certainly not an improvement.

[fredrikekre]

Alternative suggestions welcome.

[stevengj]

My suggestion is to keep spdiagm

[fredrikekre]

Here is another inconsistency:

julia> diagm(rand(3), 4)
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.519002  0.0       0.0    
 0.0  0.0  0.0  0.0  0.0       0.812178  0.0    
 0.0  0.0  0.0  0.0  0.0       0.0       0.64331
 0.0  0.0  0.0  0.0  0.0       0.0       0.0    
 0.0  0.0  0.0  0.0  0.0       0.0       0.0    
 0.0  0.0  0.0  0.0  0.0       0.0       0.0    
 0.0  0.0  0.0  0.0  0.0       0.0       0.0    

julia> full(spdiagm(rand(3), 4))
3×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.696544  0.0       0.0     
 0.0  0.0  0.0  0.0  0.0       0.100363  0.0     
 0.0  0.0  0.0  0.0  0.0       0.0       0.979094

Which behavior do we want here?

[KristofferC]

Imo the first one.

[fredrikekre]

From searching a bit it seems that the option to supply the resulting size for spdiagm (point four in the top post) is mostly used to make sure a square matrix is returned. So if we always return a square matrix, then we can probably remove the option to supply the resulting size.

[fredrikekre]

I think we are done here -- now we have diagm and spdiagm that behaves identically. Leaving the rest to JuliaLang/julia#11557, JuliaLang/julia#16740.

[KristofferC]

Good job!