Solving sparse linear system returns a Matrix

[A-C1]

Solving the sparse system is returning a Matrix while converting to dense and then solving returns a Vector as expected. Discourse discussion here.

using SparseArrays

# Constructing sparse array 
function hess_sparse(x)
    return [-sin(x[1] + x[2]) + 1, -sin(x[1] + x[2]), -sin(x[1] + x[2]), -sin(x[1] + x[2]) + 1, 1.0, 1.0, 12*x[5]^2 + 1.0, 1.0]
end
rowval = [1, 1, 2, 2, 3, 4, 5, 6]
colval = [1, 2, 1, 2, 3, 4, 5, 6]

# Constructing sparse vec
function grad_sparse(x::Vector{T}) where T<: Number
    return [cos(x[1] + x[2]), cos(x[1] + x[2]), 2*x[3], 1/2, 4*x[5]^3, 1/2]
end
gradinds = [1, 2, 3, 4, 5, 6]


x0 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
n = length(x0)
hess_mat = sparse(rowval, colval, hess_sparse(x0), n, n) 
grad_vec = sparsevec(gradinds, grad_sparse(x0), n)

# Dense linear solve returns vector 
@show x = Array(hess_mat) \ Array(grad_vec)

# Sparse linear solve returns  matrix
@show x = hess_mat \ grad_vec

[inkydragon]

julia V1.6.7

julia> hess_mat \ Vector(grad_vec)
6-element Vector{Float64}:
 0.5083672767925531
 0.5083672767925533
 2.0
 0.5
 0.3076923076923077
 0.5

julia> hess_mat \ grad_vec
ERROR: MethodError: no method matching lu!(::SparseMatrixCSC{Float64, Int64}, ::Val{true}; check=true)
Closest candidates are:
  lu!(::StridedMatrix{T}, ::Union{Val{true}, Val{false}}; check) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:79
  lu!(::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S}, ::Union{Val{true}, Val{false}}; check) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:88
  lu!(::StridedMatrix{T} where T, ::Union{Val{true}, Val{false}}; check) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:130
  ...
Stacktrace:

[LilithHafner]

This may or may not be relevant, but the docstring for \ includes "When A is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices." and A is indefinite in this case

julia> x = [1,1,0,0,0,0];

julia> transpose(x) * hess_mat * x
-1.6371897073027268

julia> y = [1,0,0,0,0,0];

julia> transpose(y) * hess_mat * y
0.09070257317431829

[dkarrasch]

First of all, I think it's great that this is already kind of working now. This error also occurs on v1.8.4, and I think this was "fixed" in the StridedMatrix relaxation PR. Second, I have a fix that makes the solve return a sparse vector.

[dkarrasch]

This may or may not be relevant, but the docstring for \ includes "When A is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices." and A is indefinite in this case

That's correct, that's why it uses an LDLt factorization which is successful in this case.

[A-C1]

In case it was missed, in some cases \ already returns a sparse vector. See example below.

using SparseArrays
using LinearSolve
using LinearAlgebra

function grad_sparse(x::Vector{T}) where T<: Number
    return [2x[1], 2x[2]]
end
gradinds = [1, 2]

function hess_sparse(x)
    vals = [2, 2]
    return vals
end
rowval = [1, 2]
colval = [1, 2]

x0 = [1.0, 1.0] 
hess_mat = sparse(rowval, colval, hess_sparse(x0))
grad_vec = sparsevec(gradinds, grad_sparse(x0))
x = hess_mat \ grad_vec # Here x is a sparse vector