LU factorization and computational linear algebra

Engineering software solves linear systems constantly. A structural analysis program assembles a stiffness matrix with hundreds of thousands of unknowns and solves it for each load case; a circuit simulator solves a linear system at every time step; a weather model solves enormous systems many times an hour. Lesson 3 showed how Gaussian elimination solves a system by hand, and Lesson 4 introduced the matrix algebra that describes it. This lesson looks at how computers actually solve linear systems efficiently and reliably. It shows that elimination is equivalent to factoring the matrix into lower and upper triangular parts, explains why that factorization is so useful, counts the cost of the operations, introduces special methods for the symmetric, banded and sparse matrices that dominate engineering, and explains how the condition number warns when a computed solution cannot be trusted.

Elimination as a factorization

Gaussian elimination reduces a matrix \(A\) to an upper triangular matrix \(U\) by subtracting multiples of earlier rows from later rows. Each such step is equivalent to multiplying by a simple elementary matrix, and the multipliers used can be stored in a lower triangular matrix \(L\) with ones on its diagonal. The result is the LU factorization:

where \(L\) is lower triangular, holding the multipliers \(\ell_{ij}\) (the amount of row \(j\) subtracted from row \(i\)), and \(U\) is the upper triangular matrix that elimination produces. Nothing new is computed: the factorization is simply a record of the elimination.

Solving with the factorization

Once \(A = LU\) is known, the system \(A\mathbf{x} = \mathbf{b}\) becomes \(L(U\mathbf{x}) = \mathbf{b}\). Write \(\mathbf{y} = U\mathbf{x}\) and solve in two easy stages:

• Forward substitution: solve \(L\mathbf{y} = \mathbf{b}\) from the top row down.
• Back substitution: solve \(U\mathbf{x} = \mathbf{y}\) from the bottom row up.

Triangular systems are cheap to solve, because each step involves only one new unknown.

In practice, rows are interchanged during elimination to put the largest available entry in the pivot position (partial pivoting, Lesson 3), which keeps rounding errors under control. The interchanges are recorded in a permutation matrix \(P\), and the factorization becomes \(PA = LU\). Every numerical library computes this factorization when asked to solve a general linear system.

The cost of solving linear systems

Counting arithmetic operations shows why factorization matters. For an \(n \times n\) matrix:

• LU factorization takes about \(\tfrac{2}{3}n^3\) floating-point operations.
• Each forward and back substitution pair takes about \(2n^2\) operations.
• Computing the inverse \(A^{-1}\) takes about \(2n^3\) operations, three times as many as the factorization, and multiplying \(A^{-1}\mathbf{b}\) costs another \(2n^2\).

The \(n^3\) growth is severe: doubling the size of a problem multiplies the factorization cost by eight. For \(n = 10{,}000\), the factorization requires about \(6.7 \times 10^{11}\) operations, a fraction of a second to a few seconds on a modern processor, while each additional right-hand side needs only \(2 \times 10^8\). This is why engineers almost never compute matrix inverses: solving with an LU factorization is cheaper and more accurate. The notation \(A^{-1}\mathbf{b}\) means "the solution of \(A\mathbf{x} = \mathbf{b}\)", not an instruction to form the inverse.

Exploiting structure

The matrices that arise in engineering are rarely arbitrary. Recognizing their structure can reduce the work by orders of magnitude.

Symmetric positive definite matrices and Cholesky factorization

Stiffness matrices, conductance matrices and many others are symmetric (\(A^T = A\)) and positive definite: \(\mathbf{x}^T A\mathbf{x} > 0\) for every non-zero \(\mathbf{x}\), which physically means that any non-zero displacement stores positive energy. Such matrices have a Cholesky factorization \(A = R^T R\), with \(R\) upper triangular, which needs no pivoting and costs about \(\tfrac{1}{3}n^3\) operations, half the cost of LU. Failure of the Cholesky factorization is itself informative: in a structural model it usually means the structure is not properly supported, so that some displacement stores no energy (a mechanism).

Banded and tridiagonal matrices

In many problems each unknown interacts only with its neighbours: the temperatures at points along a rod, the displacements of masses connected in a chain, the nodes of a long transmission line. The resulting matrices are banded, with non-zero entries only near the diagonal. For a tridiagonal matrix, elimination involves only a few operations per row, and the Thomas algorithm solves the system in about \(8n\) operations instead of \(\tfrac{2}{3}n^3\).

Sparse matrices

Large engineering models produce sparse matrices, in which almost all entries are zero: a finite element model with a million unknowns might have only about 30 non-zero entries per row. Storing only the non-zeros saves enormous memory, and sparse direct solvers reorder the unknowns to limit fill-in, the creation of new non-zeros during elimination. Alternatively, iterative methods avoid factorization entirely.

Iterative methods

Iterative methods start from a guess and improve it repeatedly, using only multiplications by the matrix, which are cheap for sparse matrices. The simplest, the Jacobi method, solves each equation for its own diagonal unknown using the previous values of the others; the Gauss-Seidel method uses each new value as soon as it is computed. Both converge if the matrix is strictly diagonally dominant, meaning each diagonal entry is larger in magnitude than the sum of the other entries in its row, a property many physical models have.

Conditioning: when can a solution be trusted?

Lesson 3 showed that a nearly singular system can magnify small changes in its data into large changes in its solution. The condition number \(\kappa(A)\) measures this sensitivity. Roughly, if the data (the entries of \(A\) or \(\mathbf{b}\)) have relative errors of size \(\epsilon\), the solution can have relative errors up to about \(\kappa(A)\,\epsilon\):

where \(\|\cdot\|\) denotes a vector or matrix norm (a measure of size). The condition number is at least 1; a value of \(10^k\) means that up to \(k\) significant digits can be lost. Since double-precision arithmetic carries about 16 digits, a matrix with \(\kappa = 10^{12}\) leaves only about four trustworthy digits, and if the input data are known to only three digits, the answer may be meaningless.

A small residual \(\mathbf{r} = \mathbf{b} - A\hat{\mathbf{x}}\) does not guarantee an accurate solution \(\hat{\mathbf{x}}\) when the matrix is ill-conditioned: many very different vectors nearly satisfy the equations. Numerical libraries can estimate the condition number cheaply alongside the solution, and responsible engineering software warns when it is large.