Implement the Cholesky decomposition

[astrozot]

Doing the Cholesky decomposition with Tensors.jl is significantly slower than using StaticArrays.jl, presumably because there is no specific code and everything is converted into a standard Matrix:

julia> using Tensors, StaticArrays, LinearAlgebra, BenchmarkTools

julia> A = SymmetricTensor{2,2,Float64}((1.0, 2.0, 5.0))
2×2 SymmetricTensor{2, 2, Float64, 3}:
 1.0  2.0
 2.0  5.0

julia> B = SMatrix{2,2}(1.0, 2.0, 2.0, 5.0)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 2.0  5.0

julia> @benchmark cholesky($A)
BenchmarkTools.Trial: 10000 samples with 992 evaluations per sample.
 Range (min … max):  38.306 ns …  3.228 μs  ┊ GC (min … max): 0.00% … 98.18%
 Time  (median):     43.557 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   46.372 ns ± 71.473 ns  ┊ GC (mean ± σ):  6.35% ±  4.48%

                   ▇█▆▂▂
  ▁▅█▇▇▇▅▄▃▃▃▄▃▃▃▄▅█████▆▄▃▃▃▃▃▃▂▂▂▂▃▃▃▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  38.3 ns         Histogram: frequency by time        55.4 ns <

 Memory estimate: 112 bytes, allocs estimate: 2.

julia> @benchmark cholesky($B)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  1.708 ns … 13.750 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.792 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.864 ns ±  0.199 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

       ▅     █     ▁     ▁     ▂    ▄     ▇                  ▁
  ▅▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁██▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁█ █
  1.71 ns      Histogram: log(frequency) by time     2.12 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[KnutAM]

Out of curiosity - what is the use case for this instead of using inv directly?

[astrozot]

@KnutAM I am not sure I understand your question: the matrix inverse and the Cholesky decomposition are quite different operations. The latter is used, for example, to generate multivariate normal (Gaussian) random vectors with a given covariance matrix.

[KnutAM]

Thanks for the clarification! I've only used it for solving equation systems.

AFAIU, we would need a custom overload for all factorizations (as we do with eigen already), and the list is quite long. Or maybe we can do something clever- that would be nice.

I believe we have the internal tosmatrix in src/eigen.jl, perhaps a viable way could be to make this public.