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

General Usage Performance
10 
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.

1 Like