Description
nonzeros
is not defined for several structured matrix types (Symmetric
, Tridiagonal
, etc. as shown below).
It also has behavior that is accurately described by the docs, but counter to what would be expected, for Triangular
matrices. Specifically, the underlying matrix .data
is called directly by nonzeros
, with no attempt to filter out non-zero values that appear in the wrong half of the matrix.
The docs state:
nonzeros(A)
Return a vector of the structural nonzero values in sparse array A. This includes zeros that are explicitly stored in the sparse array. The returned vector points
directly to the internal nonzero storage of A, and any modifications to the returned vector will mutate A as well. See rowvals and nzrange.
which is what the function does, but not what I think users would expect from the name (as they would expect it to only return non-structural-zeros that fall in the correct half).
See the following examples:
julia> x = sprand(10, 10, .05)
10×10 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
⋅ ⋅ ⋅ ⋅ 0.711598 ⋅ 0.200245 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.884361 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.831008 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
julia> u = UpperTriangular(x)
10×10 UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}:
0.0 0.0 0.0 0.0 0.711598 0.0 0.200245 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 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> nonzeros(x)
4-element Vector{Float64}:
0.8843607078139244
0.8310084043926247
0.7115978155585387
0.20024470975196884
julia> nonzeros(u)
4-element Vector{Float64}:
0.8843607078139244
0.8310084043926247
0.7115978155585387
0.20024470975196884
julia> nonzeros(sparse(Matrix(u)))
2-element Vector{Float64}:
0.7115978155585387
0.20024470975196884
julia> s = Symmetric(x)
10×10 Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}:
0.0 0.0 0.0 0.0 0.711598 0.0 0.200245 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 0.0 0.0 0.0 0.0 0.0
0.711598 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.200245 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 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> nonzeros(s)
ERROR: MethodError: no method matching nonzeros(::Symmetric{Float64, SparseMatrixCSC{Float64, Int64}})
Closest candidates are:
nonzeros(::SparseMatrixCSC) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:153
nonzeros(::SubArray{Tv, 2, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {Tv, Ti, I<:AbstractUnitRange, var"#s814"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:154
nonzeros(::UpperTriangular{var"#s814", var"#s813"} where {var"#s814", var"#s813"<:(Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {I<:AbstractUnitRange, var"#s814"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}} where {Tv, Ti})}) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:155
...
Stacktrace:
[1] top-level scope
@ REPL[24]:1
julia> t = Tridiagonal(sprand(10, .1), sprand(11, .1), sprand(10, .1))
11×11 Tridiagonal{Float64, SparseVector{Float64, Int64}}:
0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.0 0.0138954 0.0127756 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.285739 0.0 0.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.487423 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0
julia> nonzeros(t)
ERROR: MethodError: no method matching nonzeros(::Tridiagonal{Float64, SparseVector{Float64, Int64}})
Closest candidates are:
nonzeros(::SparseMatrixCSC) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:153
nonzeros(::SubArray{Tv, 2, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {Tv, Ti, I<:AbstractUnitRange, var"#s814"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:154
nonzeros(::UpperTriangular{var"#s814", var"#s813"} where {var"#s814", var"#s813"<:(Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {I<:AbstractUnitRange, var"#s814"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}} where {Tv, Ti})}) at /home/mc/github/julia/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:155
...
Stacktrace:
[1] top-level scope
@ REPL[27]:1
Adding support for the structured matrix types isn't too hard. I'll do it, unless someone else is interested.
I am curious about if the decision regarding the behavior for Triangular
should be revisited. I think it should be, but of course I am happy to hear why the current behavior was decided on in the first place.