# Large matrix operations involving inversion

**URL:** <https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993>\
**Category:** Performance\
**Tags:** memory, matrices, inverse\
**Created:** [May 13, 2022, 4:10am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993 "2022-05-13T04:10:38Z")\
**Posts on this page:** 6\
**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:** [May 13, 2022, 4:10am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/1 "2022-05-13T04:10:38Z")

</div>

In a particular problem, I have the following set of three evaluations involving matrices `A`, `B`, `C` of the `SparseMatrixCSC` type as

`B' * A^{-} * B`  
`B' * A^{-} * C`  
`C' * A^{-} * C`

Here,  
`shape(A) = (n, n)`  
`shape(B) = (n, p)`  
`shape(C) = (n, 1)`

I have to deal with `n > 1_000_000` and `p` up to 15.

I am currently using the following piece of code,

`B' * (A \ B)`  
`B' * (A \ C)`  
`C' * (A \ C)`

In another case, with `n=700_000` and `p = 9` evaluation took about 20 minutes on a system with 64GB RAM, 32 logical processors and 2.4GHz specifications.

On the same machine, `n = 500_000` along with `p = 12` is giving `OutOfMemoryError()`. Then, I tried on a node of the HPC cluster with 128GB RAM, where it didn’t either give `outofMemoryError()` or didn’t finish even after 2 hours.  
This is despite the fact that the number of elements in the first and second cases is in the same range (about 6\_000\_000).

Please let me know  
1. Is the procedure I am using computationally efficient?  
2. Any tricks/workarounds for dealing with large matrices where `n>1_000_000` and `p>10`.

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:** [May 13, 2022, 4:25am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/2 "2022-05-13T04:25:24Z")

</div>

How many nonzeros?

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 13, 2022, 4:29am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/3 "2022-05-13T04:29:53Z")

</div>

Is `A` structured at all? Like, is it symmetric, antisymmetric, something? If it is, there might be some tricks you could do.

Otherwise, one improvement would be to form `A_lu = lu(A)` and use that instead of `A` in the operations. Basically it does the slow and hard part of the equation solving upfront, so you don’t have to do it three times.

---

<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:** [May 13, 2022, 4:37am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/4 "2022-05-13T04:37:12Z")

</div>

Thanks, Gustaphe.

Yes, `A` is the symmetric and positive definite matrix. Also, non-zero terms of `A` are expected to form a banded pattern about the main diagonal.

---

<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 13, 2022, 4:38am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/5 "2022-05-13T04:38:25Z")

</div>

then use `BandedMatrices.jl` it will be way faster (at least 10x).

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 13, 2022, 5:57am UTC](https://discourse.julialang.org/t/large-matrix-operations-involving-inversion/80993/6 "2022-05-13T05:57:35Z")

</div>

For general SPD matrices there’s the cholesky decomposition (replace `lu` by `chol` above). I don’t know how that compares to/interacts with bandedness. There is an implementation in BandedMatrices.jl for LU composition though, so I would probably use that.