Iterative solvers and preconditioning

When direct factorisation runs out of memory, iterative methods are the only way to solve the system.

Direct solvers like spsolve factorise the matrix, and the factors can fill in far more non-zeros than the matrix had. For very large problems (3-D grids, big graphs) Krylov iterative methods in scipy.sparse.linalg are used instead. They only need matrix-vector products:

cg: symmetric positive definite matrices,
minres: symmetric indefinite,
gmres, bicgstab, lgmres: general non-symmetric.

They return (x, info): info == 0 means converged, info > 0 means the iteration limit was hit. Control accuracy with rtol and atol (the old tol keyword was removed) and cap work with maxiter.

Convergence speed depends on the matrix's condition number, so preconditioning is the key skill: pass M, an approximation of the inverse of A that is cheap to apply. Options range from the Jacobi preconditioner (inverse diagonal), to an incomplete LU from spilu wrapped as a LinearOperator, to algebraic multigrid from the pyamg package.

A `LinearOperator(shape, matvec=...)` represents a matrix you never store: anything you can multiply by a vector, such as a convolution, an FFT-based operator or a product of factors. Iterative solvers, eigsh and svds all accept it.

For eigenvalues, eigsh(A, k, which="LM") finds the largest; for the smallest of a large sparse matrix, use shift-invert mode, sigma=0, which converges far faster than which="SM".

Example

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import LinearOperator, cg, eigsh

n = 10_000
d = np.linspace(1, 1000, n)                         # widely varying diagonal
A = sparse.diags_array([d, -0.5 * np.ones(n - 1), -0.5 * np.ones(n - 1)],
                       offsets=[0, -1, 1]).tocsr()
b = np.ones(n)

def run(M=None):
    its = []
    x, info = cg(A, b, rtol=1e-10, M=M, callback=lambda xk: its.append(1))
    return info, len(its)

jacobi = LinearOperator((n, n), matvec=lambda v: v / d)   # matrix-free preconditioner
print(run())            # (0, 414) (converged, iterations)
print(run(M=jacobi))    # (0, 18) (far fewer with preconditioning)

vals = eigsh(A, k=2, sigma=0, which="LM", return_eigenvectors=False)
print(np.sort(vals).round(3))                       # [0.294 0.582] (two smallest eigenvalues)

Remember

Krylov solvers need only matrix-vector products; pick by matrix symmetry
Preconditioning (M=) is what makes iterative solvers converge fast
Use LinearOperator for matrix-free problems and sigma=0 for small eigenvalues

Watch out for

Ignoring the returned info and using an unconverged solution.

Check yourself

cg returns info = 1000 for a large system. What does it mean and what helps most?

Official docs: docs.scipy.org/doc/scipy/reference/sparse.linalg.html

Example

Remember

Watch out for

Check yourself

Builds on

Unlocks