extract L from CholeskySolverF

[schenkch]

I am trying to get the actual lower triangular matrix L out from CholeskySolverF but I have not succeeded so far. Is it possible?

[bathal1]

Hi,
At the moment it is not possible, because the sparse representation of L lives on the GPU.
It should be easy to do it by adding a function to copy the arrays to the CPU and returning a copy of them.

That said, if you care only about L and not solving triangular systems, then you might be better off using some other tools doing Cholesky decomposition on the CPU, e.g. SciPy

[schenkch]

Thank you. I see but I am looking for a sparse solver where you can extract L. Do you know of any other sparse Cholesky python solvers apart from scikit-sparse?

[99991]

I know of three Python implementations that compute sparse Cholesky decompositions where you can extract L:

1. sksparse.cholmod.cholesky
2. scipy.sparse.linalg.spilu (technically, this is an LU decomposition, but it can also be used to compute a Cholesky decomposition if the input matrix is symmetric)
3. pymatting.ichol with discard_threshold=0 (PyMatting uses an incomplete Cholesky factorization in conjunction with conjugate gradient descent to solve large linear systems, which is more efficient, but ichol can also be used to compute the full Cholesky factorization)

Here are examples for all three:

import numpy as np
import scipy.sparse

def factorize_with_sksparse(A):
    from sksparse.cholmod import cholesky

    decomposition = cholesky(A)

    L = decomposition.L()

    assert isinstance(L, scipy.sparse.csc_matrix)

    check_close(A, L @ L.T)

def check_close(A, B, threshold=1e-10):
    difference = A - B

    error = difference.multiply(difference).sum()

    assert error < threshold, f"error {error} exceeds threshold {threshold}"

def factorize_with_scipy(A):
    from scipy.sparse.linalg import spilu

    decomposition = spilu(A)

    L = decomposition.L
    U = decomposition.U

    # Can also compute L from U if A is symmetric
    d = U.diagonal()
    L_from_U = U.T / d

    # For numerical reasons, this is not very precise
    check_close(L, L_from_U, threshold=1e-3)

    # LDL' decomposition
    check_close(A, L @ scipy.sparse.diags(d) @ L.T)

    L = L @ scipy.sparse.diags(np.sqrt(d))

    # LL' decomposition
    check_close(A, L @ L.T)

def factorize_with_pymatting(A):
    from pymatting import ichol

    decomposition = ichol(A, discard_threshold=0)

    L = decomposition.L

    check_close(A, L @ L.T)

def main():
    A = np.array([
        [1, 1, 1],
        [1, 2, 2],
        [1, 2, 3],
    ], np.float64)

    A = scipy.sparse.csc_matrix(A)

    factorize_with_sksparse(A)
    factorize_with_scipy(A)
    factorize_with_pymatting(A)

if __name__ == "__main__":
    main()