Preconditioners: ILU(0)

ILU(0) is the first stronger native preconditioner for nonsymmetric systems in v0.0.5b0. Setup is natural-order and no-fill: no reordering, no pivoting, and no new off-diagonal entries. Application is two native CSR triangular solves, so it has real setup/apply cost and should be compared by time-to-solution, not only iteration count.

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)

Nonsymmetric convection-diffusion GMRES

This tridiagonal operator has directional transport, so CG is not appropriate. GMRES with ILU(0) uses left preconditioning and checks the true residual.

n = 96
diffusion = np.float32(0.18)
convection = np.float32(0.55)
h = np.float32(1.0 / (n + 1))
lower = (-diffusion / h**2 - convection / h) * np.ones(n - 1, dtype=np.float32)
diag = (2.0 * diffusion / h**2 + convection / h + 1.0) * np.ones(n, dtype=np.float32)
upper = (-diffusion / h**2) * np.ones(n - 1, dtype=np.float32)
A_sp = scipy.sparse.diags([lower, diag, upper], [-1, 0, 1], format="csr", dtype=np.float32)
A = ms.from_scipy(A_sp)
b = mx.cos(mx.linspace(0.0, 2.0, n))

x0, info0 = linalg.gmres(A, b, rtol=5e-4, restart=8, maxiter=128, return_info=True)
M = preconditioners.ilu0(A)
xi, infoi = linalg.gmres(A, b, M=M, rtol=5e-4, restart=8, maxiter=128, return_info=True)

print("ILU0 nnz_L/nnz_U/fill:", M.nnz_L, M.nnz_U, round(M.nnz / A.nnz, 3))
print("none:", info0.status, info0.iterations, f"{info0.residual_norm:.3e}")
print("ilu0:", infoi.status, infoi.iterations, f"{infoi.residual_norm:.3e}")

ILU0 nnz_L/nnz_U/fill: 191 191 1.336
none: 128 128 2.144e+00
ilu0: 0 8 4.268e-04

When to use it

Use ILU(0) when a nonsymmetric problem needs a stronger preconditioner than Jacobi and the natural-order no-fill pattern is acceptable. The setup is CPU native; repeated application can reuse native triangular-solve infrastructure. reuse_analysis=True is available for repeated apply workloads but should be benchmarked for the target matrix and device.