Sparse Cholesky of Gram Matrix (SuiteSparse?)

Specific Domains Numerics
,
21

I’m curious as well. I discarded this “expanded system” solution some time ago since the matrices are quite large, but I can definitely give it a try to see how quickly the problem is solved!

I actually haven’t tried HSL at all for this (I actually forgot it existed…), so I will do that at some point in the near future. Thanks for the heads up!

22

You’re right. And if the matrix is very non-square, then forming AA^T first is much faster:

A = sprand(1000, 10^6, 1e-6);
A += SparseMatrixCSC(I, size(A)...)
cholgram2(A) = cholesky(A*A')
@btime cholgram($A);
@btime cholgram2($A);

gives

21.030 ms (38 allocations: 23.15 MiB)
1.880 ms (47 allocations: 339.10 KiB)

whereas a square matrix

A = sprand(10^6, 10^6, 1e-9);
A += SparseMatrixCSC(I, size(A)...)
@btime cholgram($A);
@btime cholgram2($A);

gives

207.157 ms (38 allocations: 213.74 MiB)
200.582 ms (52 allocations: 240.73 MiB)

I wonder what the point of this feature of CHOLMOD is, if it isn’t faster? Maybe it is more accurate if A is badly conditioned, since it avoids squaring the condition number (I guess)?

23

no need for an interpolation mark ($) with @btime ?

24

Interpolating global arguments to functions is always a good idea with @btime, so that there is no dynamic-dispatch overhead.

25

That is why I wonder about the impact of missing $ marks in @guille measurements.

26

These matrices are big enough that the overhead of dynamic dispatch is probably negligible, so $ interpolation is just a good habit.

2 Likes