Skip to content

nonzeros is missing for structured matrix types and has unusual behavior for (Upper|Lower)Triangular #64

Open
@mcognetta

Description

@mcognetta

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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions