# What is the most efficient way to compute the inverse of a sparse matrix?

**URL:** <https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167>\
**Category:** Performance\
**Tags:** question, linearalgebra, sparse\
**Created:** [May 1, 2023, 4:56pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167 "2023-05-01T16:56:13Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 1, 2023, 4:56pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/1 "2023-05-01T16:56:14Z")

</div>

I have a sparse matrix Q that is a covariance matrix for a Gaussian Process. I need to invert Q for computations in Bayesian full-conditional distributions. For example, I need to compute terms like \boldsymbol{\beta}^T\mathbf{X}^TQ^{-1}, where \mathbf{X}\boldsymbol{\beta} has the typical linear regression definition. I have found that `inv(sparse(Q))` returns an error stating `"The inverse of a sparse matrix can often be dense and can cause the computer to run out of memory..."` The fastest solution I have to compute the inverse of Q is using `inv(cholesky(Q))`. However, this is still taking much longer than I’d like it to.

Throughout your Julia journeys, how have you handled inverting sparse matrices?

---

<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:** [May 1, 2023, 5:05pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/2 "2023-05-01T17:05:40Z")

</div>

Rather than inverting, you almost always should use a linear solve. So, for example `Q\X`. If you want to do a bunch of solves, you can use `c=cholesky(Q)` and `c\X`.

---

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 1, 2023, 5:08pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/3 "2023-05-01T17:08:18Z")

</div>

So, if I want to compute \boldsymbol{\beta}^T\mathbf{X}^TQ^{-1}, would I be able to use `beta'*X'*Q\(beta'*X')`? Or am I completely off?

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 1, 2023, 5:08pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/4 "2023-05-01T17:08:37Z")

</div>

As the warning indicates, sparse matrices seldom have sparse inverses. The best way to perform these sorts of calculations is with a factorization of your `Q` matrix, rather than an explicit inverse. Factorizations of sparse matrices often can be sparse, which is one reason why this is preferable. Since your `Q` is positive definite, you’ve already stumbled across the `cholesky` function. You can use this factorization object directly in calculations. For example, your \beta X^T Q^{-1} calculation can be done as `cholQ = cholesky(Q)` followed by `β * X' / cholQ`.

---

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 1, 2023, 5:11pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/5 "2023-05-01T17:11:15Z")

</div>

I see. What if I wanted to do a calculation such as (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})'\*Q^{-1}\*(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})? Would I be able to use `(y-X*beta)'*(y-X*beta)/cholQ`, or am I violating some rules by bringing `cholQ` to the end?

---

<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:** [May 1, 2023, 5:14pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/6 "2023-05-01T17:14:05Z")

</div>

Matrix multiplication isn’t commutative so you have to apply it in the right order, but other than that, this all should work.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 1, 2023, 5:20pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/7 "2023-05-01T17:20:53Z")

</div>

Matrix multiplication is not commutative, so `A * B != B * A` except under specific circumstances that you should not usually count on. By extension, you cannot move the Q^{-1} to the end.

The correct way to write (y - X\beta)^\top Q^{-1} (y - X\beta) in code is `((y - X * beta)' / cholQ) * (y - X * beta)` or `(y - X * beta)' * (cholQ \ (y - X * beta))`. Note that `/` and `\` are used to invert the matrix on the right or left of a multiplication, respectively. In this case I would write this as

```julia
cholQ = cholesky(Q)
r = y - X * beta
z = r' * (cholQ \ r) # or (r' / cholQ) * r

```

to avoided computing `y - X * beta` twice, since the compiler will probably not be smart enough to see that redundancy on its own.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 4, 2023, 5:56am UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/8 "2023-05-04T05:56:03Z")

</div>

If you need to calculate a lot of times the value for different y-X \beta but fixed and already Cholesky factorize Q, then you could save some time by exploiting the symmetry of the factorization (one triangular backward substitution instead of two)

```julia
r = y - X * beta
w = cholQ.factors.PtL \ r
z=w'*w

```

and of course you can try to save some allocations using in place calculations.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 4, 2023, 2:13pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/9 "2023-05-04T14:13:16Z")

</div>

The `.factors.PtL` property must be new in `v1.9`, or is otherwise unavailable in my version of `v1.8` at least (`.factors` is simply returning the raw underlying storage and has no special properties). In the meantime, if one is using a non-pivoted Cholesky decomposition one can use `w = cholQ.L \ r`. If pivoted, one can use `w = cholQ.L \ r[cholQ.p]`.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 4, 2023, 2:55pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/10 "2023-05-04T14:55:42Z")

</div>

Sorry, I was copying a piece of code on a version on which we added a layer.

```julia
w = cholQ.PtL \ r

```

should work

---

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 5, 2023, 3:32pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/11 "2023-05-05T15:32:12Z")

</div>

Is there any way that I can use Cholesky factorizations to reduce the amount of computation (RAM) needed to compute rQr'? Right now, my `r` is `576x165600`, so, `w` is `576x165600`, and my 64GB of RAM is still not enough to compute `z=w'*w`.

Maybe @mikmoore would have an idea, too.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 5, 2023, 3:44pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/12 "2023-05-05T15:44:26Z")

</div>

I hadn’t realized that `r` (and `y` and `beta`) were matrices and not vectors.

The issue is that you’re trying to compute a `165600x165600` matrix, which has 27 billion elements. For a dense matrix at `Float64` precision this is 220GB. If `w` is sparse, you can call `sparse(w)` to make it sparse again and then `w'*w` will be sparse and might be much smaller.

But I’m not sure I would expect `w` to actually be sparse. Further, you’re computing a big matrix which is `165600x165600` but is only rank `576` (or less), which is almost always the wrong thing to do even if it is sparse.

I would try to refactor your following calculations to exploit this. For example, instead of computing `(w' * w) * v`, which will fail on `w' * w`, compute `w' * (w * v)` which will succeed unless `v` has a large number of columns.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 5, 2023, 3:45pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/13 "2023-05-05T15:45:57Z")

</div>

Can you clarify ? Is it r' Q^{1} r or r Q^{-1} r' you calculate, or in other words what is the size of Q.

Also, what do you want to use z for ? If you multiply it with a vector for example, z\* v you could do (w'\*(w\*v)) and save memory

---

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 5, 2023, 4:01pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/14 "2023-05-05T16:01:06Z")

</div>

Surprisingly, `w` is sparse, but I agree that computing the big matrix seems wrong.

Are you suggesting that `v` is a variable that I already mentioned, or was that just for an example?

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 5, 2023, 4:03pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/15 "2023-05-05T16:03:15Z")

</div>

Just an example. It is rare that one is interested in a huge matrix z on itself. Generally you want to do something with it and the suggestion is to look how to organize the following calculations, exemplified here with a multiplication by a vector.

---

<div class="post-metadata">

**Author:** ![schwob](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schwob/32/49423_2.png) [@schwob](https://discourse.julialang.org/u/schwob)\
**Post date:** [May 5, 2023, 4:06pm UTC](https://discourse.julialang.org/t/what-is-the-most-efficient-way-to-compute-the-inverse-of-a-sparse-matrix/98167/16 "2023-05-05T16:06:57Z")

</div>

I am trying to compute r'Q^{-1}r, so Q is `576x576`.

Ultimately, I am going to use `z` for a symmetric Woodbury inverse. So [for reference](https://en.wikipedia.org/wiki/Woodbury_matrix_identity), r=(A^{-1}U)', Q=C^{-1} + U'A^{-1}U, and V=U'.