# Inversion and multiplication of sparse matrices

**URL:** <https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819>\
**Category:** Performance\
**Tags:** sparse, matrices\
**Created:** [December 30, 2021, 3:07pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819 "2021-12-30T15:07:12Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Rajesh\_Nakka](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rajesh_nakka/32/28445_2.png) [@Rajesh\_Nakka](https://discourse.julialang.org/u/Rajesh_Nakka)\
**Post date:** [December 30, 2021, 3:07pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/1 "2021-12-30T15:07:12Z")

</div>

I have encountered a problem of the following kind.

> C = A^-1 \* B

where `A` is a very large `nxn` sparse banded symmetric matrix (of `SparseMatrixCSC` type) and `B` is `nx5` matrix.

I have read on this forum and somewhere else that inverting a large sparse matrix gives a dense matrix and is computationally expensive.  
So, I am thinking of taking each column of `B` at a time and solving individually and then assembling solutions into a matrix `C`.

```julia
A * x1 = b1
A * x2 = b2
A * x3 = b3
A * x4 = b4
A * x5 = b5
C = [x1 x2 x3 x4 x5]

```

Is this the right approach?  
Is there any efficient way in Julia or in general to get matrix C?

Thanks in advance.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 30, 2021, 3:19pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/3 "2021-12-30T15:19:38Z")

</div>

Just use `C = B\A` That will automatically use a sparse solve.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [December 30, 2021, 3:36pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/4 "2021-12-30T15:36:48Z")

</div>

Do you mean `A \ B`?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 30, 2021, 3:39pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/5 "2021-12-30T15:39:26Z")

</div>

Oops

---

<div class="post-metadata">

**Author:** ![Rajesh\_Nakka](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rajesh_nakka/32/28445_2.png) [@Rajesh\_Nakka](https://discourse.julialang.org/u/Rajesh_Nakka)\
**Post date:** [December 30, 2021, 5:03pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/6 "2021-12-30T17:03:35Z")

</div>

Note: As I am interested in checking if the procedure works, I have used the fully populated matrices in sparse matrix form, ignoring the computational cost.

> julia\> using SparseArrays  
> julia\> A = sparse(rand(4,4))  
> 4×4 SparseMatrixCSC{Float64, Int64} with 16 stored entries:  
> 0.320401 0.966415 0.164792 0.222872  
> 0.53768 0.54209 0.197565 0.0362941  
> 0.897137 0.242278 0.3575 0.882516  
> 0.250789 0.872864 0.955418 0.0661719
> 
> julia\> b = sparse(rand(4,3))  
> 4×3 SparseMatrixCSC{Float64, Int64} with 12 stored entries:  
> 0.0658833 0.46049 0.835914  
> 0.331482 0.703552 0.915843  
> 0.936023 0.367614 0.172021  
> 0.783369 0.502601 0.214066

> julia\> A \ b  
> ERROR: MethodError: no method matching ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}, ::SparseMatrixCSC{Float64, Int64})  
> Closest candidates are:  
> ldiv!(::LinearAlgebra.LU{T, LinearAlgebra.Tridiagonal{T, V}}, ::AbstractVecOrMat) where {T, V} at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:595  
> ldiv!(::Number, ::AbstractArray) at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\LinearAlgebra\src\generic.jl:271  
> ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{T}, ::StridedVecOrMat{T}) where T\<:Union{Float64, ComplexF64} at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\SuiteSparse\src\umfpack.jl:751

I am getting a similar error with `b \ A`

> julia\> b \ A  
> ERROR: MethodError: no method matching ldiv!(::SuiteSparse.SPQR.QRSparse{Float64, Int64}, ::SparseMatrixCSC{Float64, Int64})  
> Closest candidates are:  
> ldiv!(::LinearAlgebra.LU{T, LinearAlgebra.Tridiagonal{T, V}}, ::AbstractVecOrMat) where {T, V} at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:595  
> ldiv!(::Number, ::AbstractArray) at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\LinearAlgebra\src\generic.jl:271  
> ldiv!(::LinearAlgebra.LDLt{\<:Any, \<:LinearAlgebra.SymTridiagonal}, ::AbstractVecOrMat) at C:\Users\rajeshnakka\AppData\Local\Programs\Julia170\share\julia\stdlib\v1.7\LinearAlgebra\src\ldlt.jl:172

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [December 30, 2021, 8:11pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/7 "2021-12-30T20:11:51Z")

</div>

Coming from Matlab, I was expecting the same behavior too.

As far as I can tell, Julia’s `\` accepts as right-hand side `b` only vectors like `rand(4)` or single-column matrices like `rand(4,1)`. The result is of the same type as `b`.  
The vector form seems to be preferred, not all linear solvers accept a single-column matrix.

See also [Avoiding memory allocation for solves with multiple right-hand sides - General Usage / Performance - JuliaLang](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879)

```julia
using LinearAlgebra, SparseArrays

function solve_multiRHS!(A, b)
    lu!(A)
    for i = 1:length(b)
        b[i] = A\b[i]
    end
    return b
end

A = rand(4,4)
b = Vector{Float64}[]
for i = 1:3
    push!(b, rand(4))
end 

solve_multiRHS!(A, b)

```

The `sparse` operator, when needed, is then applied to `A` and the _individual_ vectors of `b`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 30, 2021, 9:43pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/8 "2021-12-30T21:43:23Z")

</div>

> [@Bardo](#):
>
> As far as I can tell, Julia’s `\` accepts as right-hand side `b` only vectors like `rand(4)` or single-column matrices like `rand(4,1)` .

Nope. `sparse(rand(4,4)) \ rand(4,3)` (i.e. a matrix of right-hand sides) works just fine.

The key point here is that while the system matrix `A` can be sparse, the right-hand-side vector (or matrix!) must be dense — Julia doesn’t have any method to solve with a sparse right-hand-side.

The basic reason for this is not specific to Julia: you _can’t_ generally gain anything from sparse right-hand sides in linear algebre, since once you perform the solve the solution is no longer sparse. In fact, this is exactly the same reason why you don’t want to compute `inv(A)` or `A^-1` for a sparse `A` — matrix inversion is equivalent to solving with a right-hand side equal to the identity matrix (which is very sparse!) but the resulting inverse matrix is almost always dense.

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [December 30, 2021, 9:52pm UTC](https://discourse.julialang.org/t/inversion-and-multiplication-of-sparse-matrices/73819/9 "2021-12-30T21:52:21Z")

</div>

> [@stevengj](#):
>
> The key point here is that while the system matrix `A` can be sparse, the right-hand-side vector (or matrix!) must be dense — Julia doesn’t have any method to solve with a sparse right-hand-side.

Thanks again, Steven!