# Explicit Inversion of a random matrix performs the same as Symmetric & SPD

**URL:** <https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732>\
**Category:** Performance\
**Tags:** question, linearalgebra\
**Created:** [October 24, 2018, 12:28pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732 "2018-10-24T12:28:53Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![cpsarras](https://avatars.discourse-cdn.com/v4/letter/c/bcef8e/32.png) [@cpsarras](https://discourse.julialang.org/u/cpsarras)\
**Post date:** [October 24, 2018, 12:28pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/1 "2018-10-24T12:28:53Z")

</div>

Hi all,

so I was testing how Julia (1.1.0-DEV.452) handles Matrices with different kinds of properties (Symmetric, SPD, UpperTriangular etc.) and how much faster certain operations (i.e. inv(A)) perform based on those properties. I noticed that for Random, Symmetric & SPD (of size 1000x1000 Float64) performance (in terms of execution time, using the BenchmarkTools module) was pretty much the same, which is odd, considering that, for example, using the Cholesky factorization on an SPD and then computing the inverse is much faster.

Then I decided to use the Profile module to see exactly what functions run in each case. As you can see in the screenshot, in both cases LU factorization is called, even though that is not optimal for SPD.

 ![screen](https://global.discourse-cdn.com/julialang/original/3X/6/3/63163bb3577249d6921ae26b31abe1361e9dd291.png)

Is this expected behavior?

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [October 24, 2018, 7:57pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/2 "2018-10-24T19:57:45Z")

</div>

Digging into the [source](https://github.com/JuliaLang/julia/blob/44f2563d11f3a6fa56e9a059a8349df53cdc8eb8/stdlib/LinearAlgebra/src/symmetric.jl#L465-L466),

> Bunch-Kaufman solves can not utilize BLAS-3 for multiple right hand sides, so using LU is faster for AbstractMatrix right hand side.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [October 24, 2018, 7:59pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/3 "2018-10-24T19:59:41Z")

</div>

Please post the code that you used to generate the matrices. In general, allowing people to reproduce what you are seeing is crucial if you want to get any real help.

---

<div class="post-metadata">

**Author:** ![cpsarras](https://avatars.discourse-cdn.com/v4/letter/c/bcef8e/32.png) [@cpsarras](https://discourse.julialang.org/u/cpsarras)\
**Post date:** [October 24, 2018, 9:46pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/4 "2018-10-24T21:46:08Z")

</div>

I use the following code:

```julia
using BenchmarkTools
using LinearAlgebra

A = rand(1000, 1000)
B = rand(1000, 1000)
C = rand(1000, 1000)

B = B*B'
B[1,1] = -1
B = Symmetric(B)
println(typeof(B))
println(issymmetric(B))
println(isposdef(B))

C = C*C'
C = C+1000*Array{Float64}(I, (1000, 1000))
C = Symmetric(C)
println(typeof(C))
println(issymmetric(C))
println(isposdef(C))

println(@benchmark inv($A))
println(@benchmark inv($B))
println(@benchmark inv($C))

```

---

<div class="post-metadata">

**Author:** ![simonbyrne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonbyrne/32/19_2.png) [@simonbyrne](https://discourse.julialang.org/u/simonbyrne)\
**Post date:** [October 24, 2018, 10:06pm UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/5 "2018-10-24T22:06:27Z")

</div>

> [@cpsarras](#):
>
> using the Cholesky factorization on an SPD and then computing the inverse is much faster.

You’re correct, but the problem is how do you determine if a symmetric the matrix is positive-definite? You need to try to compute the Cholesky factor.

On your 1000x1000 matrices I get roughly

- 10ms to compute Cholesky, 20ms to invert
- 20ms to compute LU, 40ms to invert
- 17ms to compute Bunch-Kaufman, 55 ms to invert

(it also takes ~1.4ms to check symmetry if you don’t use a `Symmetric` wrapper)

So the current choice gives the best “worst-case” behaviour (i.e. if your matrix is not positive definite). Another option would be to try to compute the Cholesky, then fallback on the LU if that fails. That would slightly increase the worst-case time (by 10ms), but give a large decrease to the positive-definite case (by 30ms), so that might be worth trying.

Of course, if you know that your matrix is positive definite, you can just call `inv(cholesky(C))` which will give you the optimal performance.

---

<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:** [October 25, 2018, 12:58am UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/6 "2018-10-25T00:58:21Z")

</div>

> [@simonbyrne](#):
>
> You’re correct, but the problem is how do you determine if a symmetric the matrix is positive-definite? You need to try to compute the Cholesky factor.

In real applications I’ve never used Cholesky unless I knew _a priori_ that the matrix was SPD. Usually SPD matrices are SPD for a reason — there is something about your problem that ensures this property, and it is a good idea to know about it if you want to understand what is going on.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [October 25, 2018, 1:43am UTC](https://discourse.julialang.org/t/explicit-inversion-of-a-random-matrix-performs-the-same-as-symmetric-spd/16732/7 "2018-10-25T01:43:18Z")

</div>

I would also never just use `\` or `inv` with an SPD matrix (actually, probably `inv` in general). _Now_. But I think it is important to consider the use cases of `\` and `inv`. They are not for seasoned users, who will just use the best factorization directly given the structure of the matrix. Rather, I think the main use cases for `\` and `inv` are prototyping and language promotion. In the latter case, they can be a litmus test for the language.

So why do MathWorks choose to have a deeper decision tree (see [Solve systems of linear equations Ax = B for x - MATLAB mldivide \](https://www.mathworks.com/help/matlab/ref/mldivide.html), algorithms section)? What is the right thing to optimize for, if not worst-case performance?

Edit: note that Matlab’s `inv` for dense matrices uses LDL for hermitian, LU otherwise, [Matrix inverse - MATLAB inv](https://www.mathworks.com/help/matlab/ref/inv.html?searchHighlight=inv&s_tid=doc_srchtitle).