diag specializations for structured matrices yield dense vectors

[Sacha0]

In #411 and related triage, discussion settled on returning SparseVector from diag(A::SparseMatrixCSC) and diag(A::SparseMatrixCSC, k::Integer). JuliaLang/julia#23261 implemented that behavior. In contrast, the diag specializations for structured matrices still return Vector. Should this discrepancy be remedied, or should diag continue to return Vector for structured matrices?r Best!

[fredrikekre]

I don't think we should return SparseVector or Vector unconditionally, but rather always return a vector of the wrapped type. I think this is just something I forgot to update in the JuliaLang/julia#22718, JuliaLang/julia#22925, JuliaLang/julia#23035, JuliaLang/julia#23154 PR series and should be updated. Consider the following:

x, y, z = sparse(rand(5)), sparse(rand(4)), sparse(rand(4))
D = Diagonal(x)
B = Bidiagonal(x, y, :U)
T = Tridiagonal(y, x, z)
S = SymTridiagonal(x, y)

TT = typeof(x)

# Diagonal
typeof(diag(D))    == TT # true
typeof(diag(D, 1)) == TT # false

# Bidiagonal
typeof(diag(B))    == TT # true
typeof(diag(B, 1)) == TT # true
typeof(diag(B, 2)) == TT # false

# Tridiagonal
typeof(diag(T))    == TT # true
typeof(diag(T, 1)) == TT # true
typeof(diag(T, 2)) == TT # false

# SymTridiagonal
typeof(diag(S))    == TT # true
typeof(diag(S, 1)) == TT # true
typeof(diag(S, 2)) == TT # false

I think all of these should be true. Because I did not update that, we also have:

julia> @inferred diag(S, 1)
ERROR: return type SparseVector{Float64,Int64} does not match inferred return type Union{Array{Float64,1}, SparseVector{Float64,Int64}}
Stacktrace:
 [1] error(::String) at ./error.jl:33

[Sacha0]

That suggestion sounds great! Two questions:

(1) Vector should be the most common wrapped type by far. diag(..., k) for all k outside the structured matrix's bandwidth will then return a Vector of all zeros. That result seems subideal, but perhaps fine (having been the status quo so far). Thoughts?

(2) How would you handle the wrapped type being a <:AbstractRange? Specifically, what would you return for diagonals outside the structured matrix's bandwidth?

Best!

[fredrikekre]

(1) Yes, that is true, and could be solved by always returning a SparseVector, but then we should do that even for diagonal 0. This would be less performant than simply returning D.diag. On the other hand, I suspect that it is most common to extract the diagonals which are actually stored, and it that case it does not really matter what you get for calls like diag(D::Diagonal, 4) since you never will call it.

(2) AbstractRanges has other limitations too when wrapped in structured matrices, so if you want everything to work properly you need to do something like Diagonal(Vector(1:3)). For instance convert(Tridiagonal, ::Diagonal) hits a call to similar with UnitRange, which yields a Vector, see https://github.com/JuliaLang/julia/blob/4253e5622e7172ea1d95aaff448176a7dc8e37b5/base/linalg/special.jl#L13-L15. Such a call to similar is also what is missing in diag(::[Bi][Tri][SymTri]diagonal, k::Integer). Instead we call zeros there https://github.com/JuliaLang/julia/blob/4253e5622e7172ea1d95aaff448176a7dc8e37b5/base/linalg/tridiag.jl#L141

[Sacha0]

(1) Yes, that is true, and could be solved by always returning a SparseVector, but then we should do that even for diagonal 0. This would be less performant than simply returning D.diag. On the other hand, I suspect that it is most common to extract the diagonals which are actually stored, and it that case it does not really matter what you get for calls like diag(D::Diagonal, 4) since you never will call it.

👍

A tangential question comes to mind. What is diag's contract with respect to aliasing? Best!

[fredrikekre]

A tangential question comes to mind. What is diag's contract with respect to aliasing? Best!

diag(A, k::Integer=0) should probably have the same behaviour as getindex? In fact, for Matrix, that is already the case (of course). We can still use D.diag in places where that is OK, but perhaps diag should follow the getindex semantics? diag is, after all, just a special case of getindex. Thoughts? If we change getindex to return a view for Array some day we should of course do the same for structured matrices.

[Sacha0]

To make certain I follow, you propose that diag(A[, k]) not alias A? If so, that sounds reasonable on this end :).

[fredrikekre]

To make certain I follow, you propose that diag(A[, k]) not alias A? If so, that sounds reasonable on this end :).

Yes, basically I think diag should behave the same as getindex. If we do that, then we could also reconsider the original issue here: return SparseVector or the vector type which is wrapped. I still think that for structured matrices the most common operation is to extract the stored diagonals, so I'm in favor of returning a VectorType for diag(D::Diagonal{T,VectorType}[, k]) even if this is a waste of memory for non-stored k's, but as I said, I don't think this is such a common operation.

[Sacha0]

Sounds good on this end! :)