[Question❓] about pymatting

[99991]

Sounds like a bad idea, but if you know what you are doing, you can just use the ichol(A).L.todense(). Note that this is a colossal waste of memory and if all you want is to solve a linear system like $A x = L L^T x = b$, then you probably want to do decomposition = ichol(A); x = decomposition(b) instead, which executes backwards substitution on the sparse $L$ and $L^T$ efficiently:

pymatting/pymatting/preconditioner/ichol.py

Line 177 in b5f03a2

def __call__(self, b):

Disclaimer: Bad example, don't do this!

import numpy as np
import scipy.sparse
from pymatting import ichol

# Create a symmetric, positive definite matrix
n = 100
np.random.seed(0)
A = np.random.rand(n, n)
# Remove a few values of A to make it more sparse
A[np.random.rand(*A.shape) < 0.98] = 0
A += A.T + n * np.eye(n)

A_sparse = scipy.sparse.csr_matrix(A)

decomposition = ichol(A_sparse, discard_threshold=0.0)

# This constructs a scipy.sparse.csr_matrix
L_sparse = decomposition.L

# The following is not recommended because it constructs
# a dense matrix from a sparse matrix, which is a colossal waste of memory.
# If you just want to solve a linear system like
# A x = L L^T x = b,
# call the decomposition object like `x = decomposition(b)`
L_dense = L_sparse.toarray()

# Compute element-wise mean squared error to dense Cholesky decomposition
mse = np.mean(np.square(L_dense - np.linalg.cholesky(A)))

# Should be very small
print(mse)