Preconditioners: IC(0)

IC(0) is the SPD incomplete-factor preconditioner for PCG. It preserves the lower sparse pattern without fill and applies L y = r, then L.T z = y using native triangular solves. Default checks reject nonsymmetric input and non-positive pivots; shifts must be explicit.

import mlx.core as mx
import numpy as np
import scipy.sparse
import mlx_sparse as ms
from mlx_sparse import linalg
from mlx_sparse.linalg import preconditioners

# Use CPU execution throughout.
ms.use_cpu(require_available=False)
np.set_printoptions(precision=4, suppress=True)

2-D Poisson PCG

The five-point Poisson operator is a standard SPD benchmark. IC(0) is stronger than Jacobi on many PDE-like problems, but it also costs CPU setup plus two triangular solves per application.

grid = 12
T = scipy.sparse.diags(
    [-np.ones(grid - 1), 4.0 * np.ones(grid), -np.ones(grid - 1)],
    [-1, 0, 1],
    format="csr",
    dtype=np.float32,
)
I = scipy.sparse.eye(grid, format="csr", dtype=np.float32)
Y = scipy.sparse.diags([-np.ones(grid - 1), -np.ones(grid - 1)], [-1, 1], shape=(grid, grid), format="csr")
A_sp = (scipy.sparse.kron(I, T) + scipy.sparse.kron(Y, I)).astype(np.float32)
A = ms.from_scipy(A_sp)
b = mx.ones((A.shape[0],), dtype=mx.float32)

x0, info0 = linalg.cg(A, b, rtol=1e-4, maxiter=512, return_info=True)
Mj = preconditioners.jacobi(A, check=True)
xj, infoj = linalg.cg(A, b, M=Mj, rtol=1e-4, maxiter=512, return_info=True)
Mic = preconditioners.ichol0(A)
xic, infoic = linalg.cg(A, b, M=Mic, rtol=1e-4, maxiter=512, return_info=True)

print("IC0 nnz_L/fill:", Mic.nnz_L, round(Mic.nnz / A.nnz, 3))
print("none  :", info0.status, info0.iterations, f"{info0.residual_norm:.3e}")
print("jacobi:", infoj.status, infoj.iterations, f"{infoj.residual_norm:.3e}")
print("ic0   :", infoic.status, infoic.iterations, f"{infoic.residual_norm:.3e}")

IC0 nnz_L/fill: 408 0.607
none  : 0 16 2.630e-04
jacobi: 0 16 2.630e-04
ic0   : 0 8 6.199e-04

Shift policy

If a weakly diagonally dominant SPD-like matrix has a non-positive incomplete pivot, pass an explicit positive shift. mlx-sparse does not silently perturb pivots or alter missing diagonal structure.