Cannot solve sylvester equation with sparse matrix

General Usage
1

Hello, I have a peculiar situation where I am trying to solve AX + XB = C but A is sparse, for example think of A as a differential operator. But, this isn’t allowed, I have a small working example here:

begin 
    A = sparse(diagm(0=>rand(10),1=>rand(9),-1=>rand(9)))
    B = rand(2,2)
    C = rand(10,2)
    sylvester(A,B,C)
end

to which I get the following error

ERROR: MethodError: no method matching schur!(::SparseMatrixCSC{Float64, Int64})
The function `schur!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  schur!(::StridedMatrix{T}, ::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/schur.jl:348
  schur!(::UpperHessenberg{T, S} where S<:AbstractMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/schur.jl:105
  schur!(::StridedMatrix{var"#s5030"} where var"#s5030"<:Union{Float32, Float64, ComplexF64, ComplexF32})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/schur.jl:103

Stacktrace:
 [1] schur(A::SparseMatrixCSC{Float64, Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/schur.jl:157
 [2] sylvester(A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64}, C::Matrix{Float64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:1692
 [3] top-level scope
   @ REPL[81]:5
 [4] top-level scope
   @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/REPL/src/REPL.jl:1660

of course, if I omit the sparse function call everything works smoothly. From what I understand sparse matrices are just not implemented for the schur function, but is this expected or should this be a github issue, so that it is included?

1 Like
2

The algorithm in Julia’s standard library only works for dense matrices, and it’s not practical to make Schur factorization work for general sparse matrices since the result is dense. However, there are other algorithms you could try depending on your problem.

Since this is just a linear system of equations, you can pass it to any of the iterative algorithms in Julia by wrapping it in a linear operator acting on x = vec(X). In particular, if X is m \times n, you just need the linear operator L = x -> let X = reshape(x, m, n); vec(A*X + X*B); end to solve the equation Lx = \text{vec}(C). e.g. use LinearMaps.jl or LinearOperators.jl to wrap the closure L in a matrix-like object and pass it to LinearSolve.jl, Krylov.jl, KrylovKit.jl, or similar.

Alternatively, if n is sufficiently small, you can convert it into an ordinary sparse system of linear equations by Kronecker products. If x = vec(X), then (I_n \otimes A + B^T \otimes I_m) x = \text{vec}(C), which you can solve by:

x = (kron(spdiagm(ones(n)), A) + kron(B', spdiagm(ones(m)))) \ vec(C)
X = reshape(x, m, n)

With your example above, m, n = size(C) and this works fine — the Kronecker-product matrix is still quite sparse, only about 4 entries per row.

PS. sparse(diagm(0=>rand(10),1=>rand(9),-1=>rand(9))) is a terrible way to construct a sparse matrix, because it creates a dense matrix first and then converts it to a sparse one. Use spdiagm instead.

2 Likes
3

Thank you for the detailed reply. I guessed that it wasn’t an oversight but I just wanted to be sure about it.
I guess in my case since n is in fact small, constructing a “larger” sparse system makes more sense.
Also, thanks for pointing out the spdiagm!

4

Yes. In general, my advice is always to try a direct solver first before trying iterative methods. Direct methods are more straightforward — they either work or they run out of memory. Iterative methods often require more trial-and-error to get rapid convergence.

1 Like