🧮Linear Systems: Direct Methods

Cholesky, Banded and Tridiagonal Systems

Symmetric positive definite matrices and energy, A = LLᵀ, band storage and the Thomas algorithm

Crete, 1907

André-Louis Cholesky, a graduate of the École Polytechnique, is an officer of the French artillery assigned to the army's geographic service. Between 1907 and 1908 he surveys Crete; later, North Africa. A survey network is a web of triangles, with angles and distances measured more times than strictly necessary, and every measurement slightly wrong. Fitting the best coordinates to all of them by least squares produces a system of linear equations whose matrix is symmetric — and has one more property that Cholesky exploits. Around 1910 he writes up a method that solves such systems in about half the work of elimination, by taking square roots. He never publishes it. Cholesky is mortally wounded in battle in the last summer of the First World War and dies on 31 August 1918. In 1924 a fellow officer, Commandant Benoît, publishes the method in the Bulletin géodésique under Cholesky's name. A century later, it is how almost every finite-element code factorises a stiffness matrix.

A matrix is symmetric positive definite (SPD) if

A = A^{\mathsf T}

and

\mathbf x^{\mathsf T}A\mathbf x > 0

for every nonzero

\mathbf x

. That second condition sounds abstract. In engineering it is usually energy:

for a structure, $\tfrac12\mathbf u^{\mathsf T}K\mathbf u$ is the strain energy stored by displacements $\mathbf u$ — positive for any deformation of a properly supported structure;
for a resistor network, $\mathbf v^{\mathsf T}G\mathbf v$ is the power dissipated;
for heat conduction, it is the integral of $|\nabla T|^2$ .

Remove the supports and the structure can move as a rigid body with zero energy:

K

becomes only semi-definite and singular. Positive definiteness is the matrix form of 'the structure is stable'.

💡Cholesky factorisation

Every SPD matrix can be written uniquely as

A = LL^{\mathsf T}

with

L

lower triangular and positive diagonal. Column by column:

l_{jj} = \sqrt{a_{jj} - \textstyle\sum_{k<j} l_{jk}^2}, \qquad l_{ij} = \dfrac{1}{l_{jj}}\Big(a_{ij} - \textstyle\sum_{k<j} l_{ik}l_{jk}\Big),\ i > j.

Symmetry means only half the matrix is stored and computed: $\tfrac13 n^3$ operations, half of LU. Positive definiteness guarantees the quantity under every square root is positive and that no pivoting is needed — Cholesky is backward stable as it stands. And if a square root of a negative number appears, the matrix was not positive definite: Cholesky doubles as the cheapest test for it.

Derivation

Cholesky's formulas, by matching entries

The $2 \times 2$ case. Write out the product and match it with

A

entry by entry:

\begin{bmatrix} a & b \\ b & c \end{bmatrix} = \begin{bmatrix} l_{11} & 0 \\ l_{21} & l_{22} \end{bmatrix}\begin{bmatrix} l_{11} & l_{21} \\ 0 & l_{22} \end{bmatrix} = \begin{bmatrix} l_{11}^2 & l_{11}l_{21} \\ l_{11}l_{21} & l_{21}^2 + l_{22}^2 \end{bmatrix}.

Top left:

l_{11} = \sqrt a

. Off diagonal:

l_{21} = b/\sqrt a

. Bottom right:

l_{22} = \sqrt{c - b^2/a}

. The last square root is real exactly when

ac - b^2 > 0

— when the determinant is positive, which positive definiteness guarantees. When a quantity under a root comes out exactly zero — as it will for the springs below once their support is removed — the matrix is singular.

The general case. Entry

(i, j)

LL^{\mathsf T}

is row

i

L

times row

j

L

\sum_k l_{ik}l_{jk}

, and since both rows end at their diagonal the sum stops at

k = j

when

i \ge j

. Split off the last term:

a_{ij} = \sum_{k<j} l_{ik}l_{jk} + l_{ij}\,l_{jj}.

For

i = j

this is

a_{jj} = \sum_{k<j} l_{jk}^2 + l_{jj}^2

, which gives

l_{jj}

; for

i > j

it gives

l_{ij}

— the two formulas above. Working column by column, everything on the right is already known. The quantity under each square root equals the pivot that Gaussian elimination would produce, and for an SPD matrix every pivot is positive (each leading block has positive determinant, the product of the pivots so far).

📐Three springs

Springs of 150, 150 and 100 kN/m in series, fixed at one end, 5 kN at the free end. The stiffness matrix and its Cholesky factor are

K = \begin{bmatrix} 300 & -150 & 0 \\ -150 & 250 & -100 \\ 0 & -100 & 100 \end{bmatrix}, \qquad L = \begin{bmatrix} 17.32 & 0 & 0 \\ -8.660 & 13.23 & 0 \\ 0 & -7.559 & 6.547 \end{bmatrix}.

Forward and back substitution give displacements 0.0333, 0.0667 and 0.1167 m — just

5/150

, then

+5/150

, then

+5/100

, as springs in series must. Remove the support (set

K_{11} = 150

) and the third pivot becomes exactly zero: the structure is a mechanism.

It is not knowledge, but the act of learning, not possession but the act of getting there, which grants the greatest enjoyment.

– Carl Friedrich Gauss, letter to Farkas Bolyai (1808)

Bandwidth. Number the unknowns along a rod and each one is coupled only to its neighbours:

a_{ij} = 0

whenever

|i - j| > w

, the half-bandwidth. The key fact:

In Gaussian elimination or Cholesky without pivoting, the factors stay inside the band. Zeros outside the band are never touched.

So only the band needs to be stored —

n(2w+1)

numbers, or

n(w+1)

for a symmetric matrix — and the work drops from

\tfrac13 n^3

to about

n w^2

. For the beam below,

w = 2

: the work is proportional to

n

, not

n^3

. scipy.linalg.solve_banded and cholesky_banded take matrices in exactly this compressed form.

Derivation

Why elimination never leaves the band

At stage

k

, row

r

is changed only if it has a nonzero in column

k

below the pivot, which in a band means

k < r \le k + w

. The change subtracts a multiple of row

k

, whose nonzeros lie in columns

k

k + w

. So row

r

can only gain nonzeros in columns

\le k + w

, and

k + w < r + w

: inside the band on the right. On the left, columns before

k

have already been cleared in both rows. Nothing can appear outside

|i - j| \le w

The work follows: at each of the

n

stages, at most

w

rows are updated in at most

w + 1

columns, about

2w^2

operations, so a banded LU costs about

2nw^2

and banded Cholesky, using symmetry, about

nw^2

— instead of

\tfrac23 n^3

and

\tfrac13 n^3

💡The Thomas algorithm

For a tridiagonal system (

w = 1

) with sub-diagonal

a_i

, diagonal

b_i

, super-diagonal

c_i

, elimination reduces to two short sweeps: forward:

\;c'_i = \dfrac{c_i}{b_i - a_i c'_{i-1}}, \quad d'_i = \dfrac{d_i - a_i d'_{i-1}}{b_i - a_i c'_{i-1}}

backward:

\;x_n = d'_n, \quad x_i = d'_i - c'_i\,x_{i+1}

About

8n

operations and

4n

numbers of storage. It is stable without pivoting when the matrix is diagonally dominant or SPD — which the matrices of heat conduction, diffusion and cubic splines always are. Llewellyn Thomas described it in 1949.

Derivation

The Thomas algorithm is elimination, written for three diagonals

Row

i

of the system reads

a_ix_{i-1} + b_ix_i + c_ix_{i+1} = d_i

(with

a_1 = 0

and

c_n = 0

Row 1. Divide by

b_1

x_1 + c'_1x_2 = d'_1

, with

c'_1 = c_1/b_1

and

d'_1 = d_1/b_1

Row $i$ . Suppose the row above has already been reduced to

x_{i-1} + c'_{i-1}x_i = d'_{i-1}

. Use it to eliminate

x_{i-1}

: substitute

x_{i-1} = d'_{i-1} - c'_{i-1}x_i

into row

i

a_i\big(d'_{i-1} - c'_{i-1}x_i\big) + b_ix_i + c_ix_{i+1} = d_i \quad\Longrightarrow\quad (b_i - a_ic'_{i-1})\,x_i + c_ix_{i+1} = d_i - a_id'_{i-1},

and divide by

b_i - a_ic'_{i-1}

. The result has the same two-term form,

x_i + c'_ix_{i+1} = d'_i

, with exactly the forward-sweep formulas above.

Back substitution. The last row has no

x_{n+1}

, so

x_n = d'_n

; each earlier row then gives

x_i = d'_i - c'_ix_{i+1}

. Every row costs a handful of operations:

O(n)

in total. Diagonal dominance keeps each divisor

b_i - a_ic'_{i-1}

away from zero, which is why no pivoting is needed.

Derivation

The second-difference stencil, from Taylor's theorem

Lesson 1-3 built first-derivative formulas by subtracting the expansions of

u(x + h)

and

u(x - h)

. Adding them instead cancels every odd derivative:

u(x + h) + u(x - h) = 2u(x) + h^2u''(x) + \frac{h^4}{12}u''''(x) + \dots

Solve for

u''

u''(x) = \frac{u(x - h) - 2u(x) + u(x + h)}{h^2} - \frac{h^2}{12}u''''(x) - \dots

-u'' \approx (-u_{i-1} + 2u_i - u_{i+1})/h^2

— the rows

(-1,\ 2,\ -1)/h^2

of the matrix below — with a truncation error proportional to

h^2

. Each equation involves a point and its two neighbours, and that is the whole reason the matrix is tridiagonal.

Applying it twice gives the fourth derivative of the beam equation. The second difference of the second differences has weights

(1, -2, 1)

combined with

(1, -2, 1)

u''''_i \approx \frac{u_{i-2} - 4u_{i-1} + 6u_i - 4u_{i+1} + u_{i+2}}{h^4},

five diagonals:

1 \cdot 1

, then

1 \cdot (-2) + (-2) \cdot 1 = -4

, then

1 + 4 + 1 = 6

, and symmetrically.

📐Heat in a rod, and a beam

Steady conduction with internal heating,

-u'' = f

[0, 1]

, discretised with central differences on

n

interior points, gives the tridiagonal system with rows

(-1,\ 2,\ -1)/h^2

. With

f = \pi^2\sin\pi x

the exact answer is

\sin\pi x

, and the Thomas algorithm gives maximum errors

8.1 \times 10^{-5}

8.2 \times 10^{-7}

and

8.2 \times 10^{-9}

for

n = 100, 1000, 10\,000

: the

O(h^2)

of lesson 1-3, at a cost proportional to

n

. Stored densely,

n = 100\,000

would need 80 GB; as three diagonals, 2.4 MB. A simply supported beam,

EI\,w'''' = q

, gives five diagonals

(1, -4, 6, -4, 1)/h^4

. For

L = 10

EI = 2 \times 10^7

N m²,

q = 10

kN/m, the banded solver's midspan deflection converges to 65.1042 mm — exactly the textbook

5qL^4/(384EI)

⚠When the band is not enough

Number a 2-D plate's grid row by row and each node is coupled to the node one full row away: the half-bandwidth is the number of nodes per row,

m

, and the band is mostly holes. Cholesky fills those holes — every zero inside the band becomes a nonzero in

L

, called fill-in. For

m^2

unknowns that is about

m^3

entries in

L

and

m^4

operations: far better than dense, far worse than the five nonzeros per row the matrix actually has.

💡Where you will use this

Structures and Geotechnics: stiffness matrices are SPD and banded; skyline and banded Cholesky solvers ran the first generations of finite-element software. Hydraulic and Thermal engineering: implicit time steps for 1-D diffusion and open-channel flow are tridiagonal, solved by Thomas at every step. Attitude Control and Navigation: Kalman filters keep covariance matrices as Cholesky factors ('square-root filters') so they stay positive definite in finite precision. Interpolation: cubic splines (Unit 7) solve a tridiagonal system.

📚Further reading

Gene Golub and Charles Van Loan, Matrix Computations (4th ed., Johns Hopkins, 2013), chapter 4, 'Special Linear Systems' — Cholesky, banded and tridiagonal solvers.
Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lecture 23, 'Cholesky Factorization'.
Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 2 — the second-difference matrix of a two-point boundary value problem, its accuracy and its stability.

↓But there is a problem

Banded solvers rescue 1-D problems completely. In 2-D and 3-D the band is wide and full of holes, and fill-in eats most of the savings. In 1956 four engineers at Boeing — M. J. Turner, R. W. Clough, H. C. Martin and L. J. Topp — needed the stiffness of a swept delta wing. They broke it into small triangular panels, wrote a stiffness matrix for each, and assembled them into one huge system. The method became the finite element method. Its matrices have millions of rows and only a handful of nonzeros in each — and the order in which you number the unknowns decides whether it can be solved at all.