# Inverting a symmetric matrix is not faster than inverting a random one

**URL:** <https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [March 27, 2023, 11:27am UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654 "2023-03-27T11:27:39Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [March 27, 2023, 11:27am UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/1 "2023-03-27T11:27:39Z")

</div>

I thought that inverting a symmetric matrix would be faster, since there are optimised algorithms and I also saw that there is a method explicity for `Symmetric`, but maybe I am overlooking some trivial things, so please help 😉

This is the method I found with `methods(inv)`

```julia
[9] inv(A::LinearAlgebra.Symmetric{<:Any, <:StridedMatrix{T} where T}) in LinearAlgebra at /Users/tamasgal/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/symmetric.jl:647

```

and in the example below `M` is a subtype of that `Symmetric` type:

```julia
julia> N = 1000; A = rand(N, N); M = Symmetric(A'A);

julia> typeof(M) <: LinearAlgebra.Symmetric{<:Any, <:StridedMatrix{T} where T}
true

```

So why is there no difference between a random matrix and a symmetric one?

```julia
julia> using BenchmarkTools, LinearAlgebra

julia> N = 100; A = rand(N, N); M = A'A;

julia> @benchmark inv(M)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 415.816 μs … 5.373 ms ┊ GC (min … max): 0.00% … 90.70%
 Time (median): 473.493 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 481.312 μs ± 105.862 μs ┊ GC (mean ± σ): 0.54% ± 2.66%

              ▂▂▃▅▅▆▇▇██▅▆▄▅▄▄▃▂▁
  ▂▁▂▂▂▃▃▄▅▆▇█████████████████████▇▇▇▅▅▅▅▄▄▄▄▃▃▃▃▃▃▂▂▂▂▃▂▂▂▂▂▂▂ ▅
  416 μs Histogram: frequency by time 575 μs <

 Memory estimate: 129.17 KiB, allocs estimate: 5.

julia> N = 100; A = rand(N, N); M = Symmetric(A'A);

julia> @benchmark inv(M)
BenchmarkTools.Trial: 9563 samples with 1 evaluation.
 Range (min … max): 437.397 μs … 5.713 ms ┊ GC (min … max): 0.00% … 90.19%
 Time (median): 509.356 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 518.800 μs ± 133.612 μs ┊ GC (mean ± σ): 0.91% ± 3.63%

                   ▂▃▄▆▅▆▇▇█▇▇▅▅▅▄▂▃▁▁
  ▁▁▁▁▁▁▁▁▂▂▂▃▄▅▆█████████████████████▇▇▇▇▅▅▄▄▄▃▃▃▃▃▂▂▂▂▂▂▂▂▁▁▁ ▄
  437 μs Histogram: frequency by time 600 μs <

 Memory estimate: 129.20 KiB, allocs estimate: 6.

```

also I see O(N^2), so I guess it’s using the Coppersmith–Winograd algorithm and not something specialised like LDU:

```julia
julia> N = 10; A = rand(N, N); M = Symmetric(A'A);

julia> @btime inv(M);
  1.442 μs (5 allocations: 6.17 KiB)

julia> N = 100; A = rand(N, N); M = Symmetric(A'A);

julia> @btime inv(M);
  146.709 μs (6 allocations: 129.20 KiB)

julia> N = 1000; A = rand(N, N); M = Symmetric(A'A);

julia> @btime inv(M);
  16.408 ms (6 allocations: 8.13 MiB)

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 27, 2023, 12:22pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/2 "2023-03-27T12:22:51Z")

</div>

Not an expert here, but when you use BenchmarkTools, don’t forget to interpolate all global variables with a dollar sign to avoid skewed results. I’m assuming it doesn’t change much in this case?

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [March 27, 2023, 12:29pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/3 "2023-03-27T12:29:46Z")

</div>

As far as I can tell from the source code, both code paths lead to an `lu` call and the Symmetric method mostly makes sure that the output also becomes Symmetric.

I suspect that you are drawing too strong complexity conclusions from your samples.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 27, 2023, 12:42pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/4 "2023-03-27T12:42:39Z")

</div>

If you follow the chain of calls (tip: use `@edit inv(Symmetric(randn(2,2)))`) you’ll see that in the end the same method (lu) is called. You can access LAPACK’s specialized symmetric routines with something like

```julia
function inv_sym(A)
    A, ipiv, info = LAPACK.sytrf!('U', A)
    LAPACK.sytri!('U', A, ipiv)
end

```

I tested on 1k x 1k, 100x100 and 20x20

With OpenBLAS, for 1k and 100, both this and `inv` took about the same time; for 20, the second (`inv`, so just lu) was about 3x faster.

With MKL, for 1k it took the same time, for 100 inv was about twice as slow, and for 20 inv was about 3x as slow.

In my experience OpenBLAS regularly displayed very strange performance characteristics. I’d start by using MKL, then maybe investigating direct LAPACK calls if you find the julia codepaths are suboptimal. (the other thing is that there are very few cases where you need explicit inverses, so maybe start with seeing if you can’t get rid of them)

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 27, 2023, 12:49pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/5 "2023-03-27T12:49:55Z")

</div>

> [@tamasgal](#):
>
> also I see O(N^2), so I guess it’s using the Coppersmith–Winograd algorithm and not something specialised like LDU:

This has nothing to do with the underlying algorithm (which is not coppersmith-winograd or any other fancy algorithm, and is just plain old N^3), and is probably just a sign the problem is bandwidth-bound. You should see N^3 scaling for much larger matrices

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [March 27, 2023, 12:49pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/6 "2023-03-27T12:49:57Z")

</div>

No worries, I checked that before of course. No problems here 😉

---

<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:** [March 27, 2023, 12:51pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/7 "2023-03-27T12:51:27Z")

</div>

> [@tamasgal](#):
>
> `M = Symmetric(A'A)`

This isn’t just symmetric, it’s symmetric positive-definite, in which case you could call the `cholesky` factorization, which really is faster than `lu`.

(Of course, the usual caveats about matrix inversion apply — most of the time, you don’t need to compute an inverse explicitly, because you can just work with the factorization instead, e.g. to repeatedly solve different right-hand-sides. When you do `A \ b` it never computes a matrix inverse.)

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [March 27, 2023, 12:55pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/8 "2023-03-27T12:55:59Z")

</div>

Thanks, I am a bit confused since the `inv_sym()` is twice as slow compared to `inv()`, at least on my M1. I don’t have access to MKL on this machine since it’s not Intel.

Anyways, I think it would make sense to use the power of (multiple) dispatch to call the optimised algorithm, otherwise I don’t really see the point of a special method for `int()`. I mean, the first thing which comes to my mind is that it’s probably choosing a better algorithm…

```julia
julia> N = 100; A = rand(N, N); M = A'A;

julia> @benchmark inv($M)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 385.735 μs … 14.834 ms ┊ GC (min … max): 0.00% … 96.37%
 Time (median): 467.223 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 477.210 μs ± 208.396 μs ┊ GC (mean ± σ): 1.18% ± 3.48%

                    ▂▄▄▄▄▅▇▇█▇█▅▄▅▄▃▂▂
  ▁▁▁▁▁▂▂▂▂▂▃▃▄▃▅▅▇▇█████████████████████▇▆▆▄▅▄▄▃▃▃▂▂▂▂▂▁▂▂▁▁▁▁ ▄
  386 μs Histogram: frequency by time 564 μs <

 Memory estimate: 129.17 KiB, allocs estimate: 5.

julia> @benchmark inv_sym($M)
BenchmarkTools.Trial: 6344 samples with 1 evaluation.
 Range (min … max): 676.431 μs … 10.389 ms ┊ GC (min … max): 0.00% … 92.72%
 Time (median): 768.999 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 784.270 μs ± 195.641 μs ┊ GC (mean ± σ): 0.64% ± 2.66%

               ▁▁▄▄▆▆▇▇▆█▆▅▅▃▂▂
  ▁▁▁▁▁▂▂▂▃▄▅▆███████████████████▆▆▆▅▄▄▃▄▃▃▂▂▃▂▂▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁ ▄
  676 μs Histogram: frequency by time 928 μs <

 Memory estimate: 51.88 KiB, allocs estimate: 4.

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 27, 2023, 1:00pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/9 "2023-03-27T13:00:06Z")

</div>

> Anyways, I think it would make sense to use the power of (multiple) dispatch to call the optimised algorithm, otherwise I don’t really see the point of a special method for `int()` .

The problem is that which algorithm is faster very much depends on your matrix size, BLAS implementation and hardware, so it’s tricky for julia to know the right thing to do. In this case, I suspect it’s just that OpenBLAS’s LDLT is just not very good. There’s the open issue [BLAS support for M1 ARM64 via Apple's Accelerate · Issue #42312 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/42312) for a faster BLAS on apple

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 27, 2023, 1:00pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/10 "2023-03-27T13:00:22Z")

</div>

```julia
julia> N = 4000; A = rand(N, N); M = Symmetric(A'A);

julia> @benchmark inv($M)
BenchmarkTools.Trial: 10 samples with 1 evaluation.
 Range (min … max): 537.033 ms … 539.106 ms ┊ GC (min … max): 0.14% … 0.08%
 Time (median): 538.214 ms ┊ GC (median): 0.14%
 Time (mean ± σ): 538.154 ms ± 673.172 μs ┊ GC (mean ± σ): 0.12% ± 0.03%

  █ █ █ █ █ █ █ █ █ █  
  █▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁█▁█▁▁▁▁▁▁▁▁█▁▁▁▁█ ▁
  537 ms Histogram: frequency by time 539 ms <

 Memory estimate: 139.68 MiB, allocs estimate: 6.

julia> N = 8000; A = rand(N, N); M = Symmetric(A'A);

julia> @benchmark inv($M)
BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max): 4.096 s … 4.107 s ┊ GC (min … max): 0.01% … 0.01%
 Time (median): 4.101 s ┊ GC (median): 0.01%
 Time (mean ± σ): 4.101 s ± 7.605 ms ┊ GC (mean ± σ): 0.01% ± 0.00%

  █ █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  4.1 s Histogram: frequency by time 4.11 s <

 Memory estimate: 529.36 MiB, allocs estimate: 6.

```

As others have noted, this uses LU which calls LAPACK `dgetrf`. LAPACK does blocking to improve performance where appropriate, but the routines opt out of this for small matrices. Somewhere in your increasing n, LAPACK would have switched to using blocking, effectively changing the details of the algorithm. And in general LAPACK is not especially optimized for smaller matrices. You need to go with some larger matrices if you want to observe the O(n^3) behavior. I’ve never heard of Coppersmith-Winograd actually being used for real matrix computations.

---

<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:** [March 27, 2023, 2:37pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/11 "2023-03-27T14:37:22Z")

</div>

WRT the quadratic time you’re seeing, that’s just because LU factorizations have high quadratic overhead and look like they have quadratic looking runtime for a shockingly long time. Looking at the plot below, you can see that for sizes between 8x8 and 512x512, the time for `inv` is almost completely quadratic, but for 2048x2048 and above, the time is almost completely cubic.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/a/fa60d4c9603535feea96764a88912707324eaa78.png)

---

<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:** [March 28, 2023, 12:28pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/12 "2023-03-28T12:28:50Z")

</div>

7 posts were split to a new topic: [Asymptotically faster matrix-multiplication algorithms](https://discourse.julialang.org/t/asymptotically-faster-matrix-multiplication-algorithms/96710)

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [March 27, 2023, 3:11pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/14 "2023-03-27T15:11:29Z")

</div>

> [@tamasgal](#):
>
> Anyways, I think it would make sense to use the power of (multiple) dispatch to call the optimised algorithm

As you observed yourself it’s not clear that the specialized algorithms necessarily are faster, so that might explain why it hasn’t been implemented.

> [@tamasgal](#):
>
> otherwise I don’t really see the point of a special method for `int()`

It allows propagation of the Symmetric type, which seems worthwhile in itself since later computations might be able make good use of it.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 27, 2023, 3:15pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/15 "2023-03-27T15:15:36Z")

</div>

> [@Oscar\_Smith](#):
>
> WRT the quadratic time you’re seeing, that’s just because LU factorizations have high quadratic overhead and look like they have quadratic looking runtime for a shockingly long time

That’s a very interesting plot. I tried to predict the crossing point with a back of the envelope computation, using data from [https://freecontent.manning.com/performance-limits-and-profiling/](https://freecontent.manning.com/performance-limits-and-profiling/) which quotes 260Gflops/s and a bandwidth of 22GB/s. Given we transfer 8 N^2 bytes and do N^3 flops on them, that suggests the crossing point is around N=100, whereas we observe a much higher crossing point. Any idea what is wrong in the above computation? I ignored caches, but cache effects if anything should decrease the crossing point, not increase it. Or is memory latency rather than bandwidth the limiting factor?

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [March 27, 2023, 3:22pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/16 "2023-03-27T15:22:26Z")

</div>

Absolutely, propagating the type (symmetric) makes of course sense, I did not think about that!

---

<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:** [March 27, 2023, 3:23pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/17 "2023-03-27T15:23:24Z")

</div>

There’s a few things you’re missing. The first is multi-threading which lets you gain extra efficiency in the 100-400 range. Also, pivoting and computing the `inv` on top of the `LU` add a decent amount of n^2 overhead.

---

<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:** [March 27, 2023, 3:25pm UTC](https://discourse.julialang.org/t/inverting-a-symmetric-matrix-is-not-faster-than-inverting-a-random-one/96654/18 "2023-03-27T15:25:13Z")

</div>

> [@Oscar\_Smith](#):
>
> computing the `inv` on top of the `LU` add a decent amount of n^2 overhead.

Computing the `inv` on top of `LU` is O(n^3).

In general, multithreading makes this a lot more difficult to make sense of, because as the problem gets bigger it starts to benefit from more threads, somewhat counteracting the O(n^3) growth. To separate out this effect, I would benchmark with a single thread:

```julia
julia> using LinearAlgebra

julia> BLAS.set_num_threads(1)

```

With a single thread, I’ve typically gotten pretty clear \Theta(n^3) LU behavior for n \gtrsim 100 in the past, e.g. in [this notebook](https://nbviewer.org/github/mitmath/1806/blob/fall22/notes/Complexity.ipynb).