🧮Linear Systems: Direct Methods

Sparse Matrices

CSR storage, assembly from elements, fill-in, reordering and sparse direct solvers

Seattle, 1956

Boeing's engineers are designing jet aircraft with thin, swept, delta-shaped wings. Beam theory, which treats a wing as a single long beam, cannot predict how such a wing bends and twists. In the Journal of the Aeronautical Sciences that September, M. J. Turner, R. W. Clough, H. C. Martin and L. J. Topp describe another way. Cut the wing skin into small triangular panels. For one triangle, the relation between the forces and displacements at its three corners is a small stiffness matrix that can be written down by hand. Then assemble: at every corner shared by several triangles, add their contributions. The result is one large system

K\mathbf u = \mathbf f

for the whole wing. Four years later Clough named it the finite element method. Today it designs bridges, dams, engine blades, foundations and spacecraft. Each corner is shared by only a few elements, so each row of

K

has only a handful of nonzeros — perhaps 20 out of a million. A dense million-by-million matrix would need 8 terabytes. Stored as nonzeros only, the same matrix fits in a few hundred megabytes.

A matrix is sparse when most of its entries are zero and — the practical definition — when it pays to exploit that. Sparse matrices come from local physics: a finite-element node couples only to nodes of the elements it touches, a bus in a power grid only to the buses it has lines to, a cell in a CFD mesh only to its faces' neighbours. The number of nonzeros grows like

n

, not

n^2

💡Compressed sparse row (CSR) storage

Store three arrays:

data — the nonzero values, row by row;
indices — the column of each value;
indptr — where each row starts in data (length $n+1$ ).

Row

i

is data[indptr[i]:indptr[i+1]] in columns indices[indptr[i]:indptr[i+1]]. A matrix–vector product touches each nonzero once:

O(\text{nnz})

work instead of

O(n^2)

. CSC is the same by columns (what sparse factorisations prefer); COO is three lists of (row, column, value) — the natural format for assembly, because duplicates are simply summed when converting to CSR. All three are in scipy.sparse.

Derivation

The five-point Laplacian, from two second differences

On a square grid with spacing

h

, write

u_{i,j}

for the value at

(ih, jh)

. Lesson 3-3's second difference works in each direction separately:

-u_{xx} \approx \frac{-u_{i-1,j} + 2u_{i,j} - u_{i+1,j}}{h^2}, \qquad -u_{yy} \approx \frac{-u_{i,j-1} + 2u_{i,j} - u_{i,j+1}}{h^2}.

Add them:

-\nabla^2u \approx \frac{4u_{i,j} - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1}}{h^2}.

Each unknown is coupled to itself and four neighbours: at most five nonzeros per row, however large the grid. Numbering the

m \times m

unknowns row by row, the 1-D matrix

T = \operatorname{tridiag}(-1, 2, -1)

acts within each grid row (giving

I \otimes T

) and between neighbouring rows (giving

T \otimes I

); their sum is the kron(I, T) + kron(T, I) used in the lab.

📐Storage for a plate

The 5-point Laplacian on an

m \times m

grid has about 5 nonzeros per row:

$m = 100$ : $10^4$ unknowns, 49,600 nonzeros, 0.6 MB in CSR, 0.8 GB dense;
$m = 300$ : $9 \times 10^4$ unknowns, 0.45 million nonzeros, 5.7 MB versus 65 GB;
$m = 1000$ : $10^6$ unknowns, 5 million nonzeros, 64 MB versus 8,000 GB.

Solving the uniformly heated square plate,

-\nabla^2 u = 1

with

u = 0

on the edges, with scipy.sparse.linalg.spsolve gives centre temperatures 0.073446, 0.073657 and 0.073670 for

m = 15, 63, 255

— converging to the series solution 0.073671 — and the largest of those systems, 65,025 unknowns, takes well under a second.

Essentially, all models are wrong, but some are useful.

– George E. P. Box, Empirical Model-Building and Response Surfaces (1987)

Fill-in and ordering. Factorising a sparse matrix creates nonzeros where

A

had zeros. In graph terms: the unknowns are nodes, nonzeros are edges, and eliminating a node connects all of its remaining neighbours to each other. Fill-in is the new edges.

The amount of fill depends entirely on the order in which unknowns are eliminated — that is, on how they are numbered. The matrix, the physics and the answer are identical; the cost is not. On a

30 \times 30

grid the nonzeros in

L + U

are:

row-by-row numbering: 54,058;
random numbering (as a mesh generator might produce): 123,208;
reverse Cuthill–McKee applied to the random numbering: 38,630;
minimum degree: 21,572 — out of 810,000 entries in a dense factor.

Derivation

Why eliminating a node joins its neighbours

Eliminating unknown

i

subtracts multiples of row

i

from every row

j

with

a_{ji} \ne 0

. The entry in row

j

, column

k

becomes

a_{jk} \;\leftarrow\; a_{jk} - \frac{a_{ji}\,a_{ik}}{a_{ii}}.

The correction is nonzero exactly when

a_{ji} \ne 0

and

a_{ik} \ne 0

— when

j

and

k

are both neighbours of

i

in the graph. So if

a_{jk}

was zero it generally becomes nonzero:

j

and

k

, which were not connected, now are. After the elimination all of

i

's remaining neighbours are connected to one another, and every new connection is a fill-in entry. A node with

d

neighbours can create up to

d(d-1)/2

of them, which is why minimum degree eliminates the node with the fewest neighbours first.

💡Reordering algorithms

Reverse Cuthill–McKee (Cuthill and McKee 1969; reversed by Alan George, 1971): a breadth-first search from a peripheral node, numbering neighbours by increasing degree. It produces a narrow band — ideal for banded and 'skyline' solvers.
Minimum degree (Tinney and Walker, 1967 — the same William Tinney whose Newton power flow closed Unit 2): at every step eliminate the node with the fewest current neighbours, so the clique it creates is small. Modern versions (AMD, COLAMD) are the default in SuperLU, UMFPACK and CHOLMOD.
Nested dissection (George, 1973): cut the mesh in half with a small separator, number the separator last, recurse. For 2-D grids it is provably near-optimal: about $n\log n$ fill and $n^{3/2}$ work.

A sparse direct solver runs in three phases: analyse (choose the ordering, predict fill symbolically), factorise (numerically), solve. The first is reused whenever the pattern repeats — every Newton step, every time step.

🕰Power grids and sparsity

Sparse matrix technology grew up in two places at once: structural analysis and electrical power systems. Tinney and Walker's 1967 ordering made Newton power flow practical for networks of thousands of buses, where the Jacobian has a handful of nonzeros per row. Harry Markowitz — later a Nobel laureate in economics, for portfolio theory — had proposed choosing pivots to limit fill for large linear programming models in 1957.

⚠When direct methods stop scaling

Even with the best ordering, fill in 3-D grows fast: nested dissection gives about

n^{4/3}

nonzeros and

n^2

work for an

n

-node 3-D mesh. At a million unknowns a sparse direct solve is routine; at a hundred million, memory for the factors runs out. Beyond that, codes stop factorising and use iterative methods that only ever multiply by

A

— the subject of the next unit.

💡Where you will use this

Structures, Geotechnics and Mechanical design: every finite-element model is a sparse SPD system, assembled in COO, solved by sparse Cholesky with a fill-reducing ordering. Power Systems: sparse LU inside Newton load flow. CFD and Hydraulics: sparse pressure equations at every time step. GNSS Navigation: large network adjustments and precise orbit determination are sparse least-squares problems.

📚Further reading

Timothy Davis, Direct Methods for Sparse Linear Systems (SIAM, 2006) — sparse storage, elimination graphs, orderings and factorisation, with complete C code.
Alan George and Joseph Liu, Computer Solution of Large Sparse Positive Definite Systems (Prentice-Hall, 1981) — the classic account of fill-in, reverse Cuthill–McKee, minimum degree and nested dissection.
M. J. Turner, R. W. Clough, H. C. Martin and L. J. Topp, 'Stiffness and deflection analysis of complex structures', Journal of the Aeronautical Sciences 23(9), 805–823 (1956).

↓But there is a problem

A sparse direct solver is exact (to rounding), robust and — with a good ordering — fast up to a few million unknowns. Then the factors no longer fit in memory. There is another way to use a matrix: never factorise it at all. Start with a guess, use the equations to improve it, and repeat. Each step costs one sweep through the nonzeros. Carl Friedrich Gauss described such a method in a letter in 1823, and remarked that it could be carried out 'half asleep, or while thinking about other things'. Next unit: how to make that idea converge — and converge fast.