Sparse Cholesky of Gram Matrix (SuiteSparse?)

Specific Domains Numerics
,
17

I did a little digging, and you can access this functionality from Julia without too much difficulty. In Julia 1.4, do:

using SparseArrays, LinearAlgebra
using SuiteSparse.CHOLMOD

# factorize A*A'
function cholgram(A::Sparse)
    cm = CHOLMOD.defaults(CHOLMOD.common_struct[Threads.threadid()])
    CHOLMOD.set_print_level(cm, 0)

    # Compute the symbolic factorization
    F = CHOLMOD.analyze(A, cm)

    # Compute the numerical factorization
    cholesky!(F, A)

    return F
end
cholgram(A::SparseMatrixCSC) = cholgram(Sparse(A, 0)) # stype == 0

and then cholgram(A) \ b is equivalent to (but much faster than) (A*A') \ b. (Note that it’s AA^T, not A^TA.)

It would be nice to support this as a keyword argument directly in Julia’s standard library. Essentially, all I did was to disable these two lines and be sure to pass 0 for stype in the Sparse constructor, so it would be an easy PR to add a keyword argument for that if someone wants to take it up.

5 Likes