Sparse matrices in depth

Build COO matrices, convert to CSR or CSC for arithmetic, solve systems with direct and iterative methods, and avoid the memory blow-up of fill-in.

Formats and conversion

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve, cg

rows = [0, 0, 1, 2, 2, 2]
cols = [0, 3, 1, 0, 1, 2]
vals = [4.0, 1.0, 3.0, 1.0, 2.0, 5.0]

A = sparse.coo_matrix((vals, (rows, cols)), shape=(3, 4))
A_csr = A.tocsr()            # row slicing and matrix-vector products
A_csc = A.tocsc()            # column slicing and most factorisations

print(A_csr.nnz, A_csr.shape)          # stored non-zeros, not total entries
A_csr.data *= 2                        # scale in place, no dense copy

Format	Build by	Good for
`coo_matrix`	Triplets, fast to assemble	Construction only; convert before arithmetic
`csr_matrix`	Rows compressed	Matrix products, row slicing, most solvers
`csc_matrix`	Columns compressed	Column slicing, LU and Cholesky factorisation
`lil_matrix`	Lists of lists	Slow incremental assignment during setup
`dia_matrix`	Diagonals	Regular stencil matrices from finite differences

Solving sparse systems

# direct solver: exact, but can fill in badly
x = spsolve(A_csc, np.ones(3))

# reuse a factorisation across right-hand sides
lu = sparse.linalg.splu(A_csc.tocsc())
x1 = lu.solve(np.ones(3))
x2 = lu.solve(np.arange(3.0))

# iterative solver: little memory, needs a well-behaved matrix
sym = (A_csr @ A_csr.T).tocsr()
x_it, info = cg(sym, np.ones(3), rtol=1e-10, atol=0.0)
print("converged" if info == 0 else "failed")

Iterative solvers need symmetry and positive definiteness for CG; use gmres or bicgstab for general matrices.
Preconditioning is what makes iterative methods practical: an incomplete LU (spilu) can cut iterations by orders of magnitude.
info != 0 means the tolerance was not reached. Treat it as an error rather than using the returned vector.
Reordering the matrix with reverse_cuthill_mckee often reduces fill-in for direct factorisation.

Sparsity-preserving practice

# WRONG: dense intermediate, memory explodes
bad = np.linalg.inv(A_csr.toarray())

# RIGHT: stay sparse end to end
good = sparse.linalg.spsolve(A_csc, b)

# check before you commit: how much would dense cost?
dense_gb = A_csr.shape[0] * A_csr.shape[1] * 8 / 1e9
print(f"dense would need {dense_gb:.1f} GB, stored nnz={A_csr.nnz}")

💡

A single dense operation such as toarray(), multiplication with a dense array, or an element-wise power can convert your matrix to a full array and exhaust memory. Check nnz and shape before any operation whose output could plausibly be dense.

FAQ

Why is my sparse solver slower than the dense one?

Direct sparse factorisation suffers fill-in: zeros become non-zeros during elimination. For a small or fairly dense matrix the overhead wins out, so use a dense solver below a few hundred rows and try reordering on larger ones.

How do I build a matrix when I do not know the size in advance?

Collect triplets in plain Python lists and construct a coo_matrix once at the end. Appending to a lil_matrix in a loop is simpler but far slower for large inputs.

Linear algebra with scipy.linalg SciPy arrays and NumPy interop

Last refreshed 2026-09-18.

Sparse matrices in depth

Formats and conversion

Solving sparse systems

Sparsity-preserving practice

FAQ

Related