Numerical Linear Algebra: LU, QR and Iterative Methods

Almost every scientific computation, from structural analysis to machine learning, eventually becomes "solve \(Ax = b\)" or "find a least-squares fit". Numerical linear algebra is the study of doing this quickly and accurately on a computer. This guide explains the three workhorses taught in BS and MS courses, LU, QR and iterative methods, with worked examples and practical advice.

Why Not Just Compute \(A^{-1}\)?

In theory \(x = A^{-1}b\). In practice, forming the inverse costs more operations than necessary and amplifies rounding errors. Numerical analysts instead factorise \(A\) into simpler pieces, typically triangular or orthogonal matrices, and solve with those. A triangular system is solved by substitution in about \(n^2\) operations.

LU Decomposition

LU decomposition writes \(A = LU\), where \(L\) is unit lower triangular (1s on the diagonal) and \(U\) is upper triangular. It is simply Gaussian elimination with the multipliers stored in \(L\). Once you have it, solving \(Ax = b\) is two easy steps:

1. Forward substitution: solve \(Ly = b\).
2. Back substitution: solve \(Ux = y\).

Factorising costs about \(\tfrac{2}{3}n^3\) floating-point operations; each extra right-hand side then costs only about \(2n^2\). That is why LU is ideal when you solve with the same matrix many times.

Worked Example 1: LU Factorisation and Solving

Factorise and solve \(Ax = b\) with

$$ A = \begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{pmatrix}, \qquad b = \begin{pmatrix} 4 \\ 10 \\ 24 \end{pmatrix}. $$

Step 1: Eliminate column 1. The pivot is \(2\). Multipliers: \(l_{21} = 4/2 = 2\), \(l_{31} = 8/2 = 4\).

• Row 2 minus \(2 \times\) Row 1: \((4, 3, 3) - (4, 2, 2) = (0, 1, 1)\).
• Row 3 minus \(4 \times\) Row 1: \((8, 7, 9) - (8, 4, 4) = (0, 3, 5)\).

Step 2: Eliminate column 2. The pivot is \(1\). Multiplier \(l_{32} = 3/1 = 3\).

• Row 3 minus \(3 \times\) Row 2: \((0, 3, 5) - (0, 3, 3) = (0, 0, 2)\).

Step 3: Record the factors.

$$ L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix}, \qquad U = \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{pmatrix}. $$

Check: the third row of \(LU\) is \(4(2,1,1) + 3(0,1,1) + (0,0,2) = (8, 7, 9)\). Correct.

Step 4: Forward substitution for \(Ly = b\): \(y_1 = 4\); \(y_2 = 10 - 2(4) = 2\); \(y_3 = 24 - 4(4) - 3(2) = 2\).

Step 5: Back substitution for \(Ux = y\): \(2x_3 = 2\) gives \(x_3 = 1\); \(x_2 + x_3 = 2\) gives \(x_2 = 1\); \(2x_1 + x_2 + x_3 = 4\) gives \(x_1 = 1\).

So \(x = (1, 1, 1)^T\). Sanity check: \(A(1,1,1)^T = (4, 10, 24)^T = b\).

Pivoting

If a pivot is zero or very small, elimination fails or loses accuracy. Partial pivoting swaps rows so the largest available entry in the column becomes the pivot, giving \(PA = LU\) with a permutation matrix \(P\). For symmetric positive definite matrices, the Cholesky factorisation \(A = LL^T\) is about twice as fast and needs no pivoting.

QR Factorisation

QR factorisation writes \(A = QR\), where \(Q\) has orthonormal columns (\(Q^TQ = I\)) and \(R\) is upper triangular. Orthogonal matrices preserve lengths, so they do not amplify errors. QR is the standard method for least-squares problems: minimising \(\|Ax - b\|_2\) reduces to solving

$$ Rx = Q^T b. $$

QR can be computed by Gram–Schmidt (easy by hand), modified Gram–Schmidt (more stable), or Householder reflections (the professional standard).

Worked Example 2: QR by Gram–Schmidt

Find the QR factorisation of \(A = \begin{pmatrix} 3 & 1 \\ 4 & 2 \end{pmatrix}\).

Step 1: First column. \(a_1 = (3, 4)^T\) has length \(r_{11} = \sqrt{9 + 16} = 5\), so \(q_1 = (0.6, 0.8)^T\).

Step 2: Project the second column. \(a_2 = (1, 2)^T\), and \(r_{12} = q_1^T a_2 = 0.6 + 1.6 = 2.2\).

Step 3: Remove the component.

$$ v = a_2 - r_{12}q_1 = (1 - 1.32,\; 2 - 1.76)^T = (-0.32,\; 0.24)^T. $$

Step 4: Normalise. \(r_{22} = \sqrt{0.1024 + 0.0576} = \sqrt{0.16} = 0.4\), so \(q_2 = (-0.8, 0.6)^T\).

Step 5: Assemble.

$$ Q = \begin{pmatrix} 0.6 & -0.8 \\ 0.8 & 0.6 \end{pmatrix}, \qquad R = \begin{pmatrix} 5 & 2.2 \\ 0 & 0.4 \end{pmatrix}. $$

Check: \(q_1 \cdot q_2 = -0.48 + 0.48 = 0\), and the second column of \(QR\) is \(2.2(0.6, 0.8) + 0.4(-0.8, 0.6) = (1.32 - 0.32,\; 1.76 + 0.24) = (1, 2)\). Correct.

Iterative Methods

For very large sparse systems, such as those from finite difference or finite element models with millions of unknowns, direct factorisation is too expensive. Iterative methods start from a guess and improve it.

Split \(A = D + L + U\) (diagonal, strictly lower, strictly upper parts):

• Jacobi: \(x^{(k+1)} = D^{-1}\big(b - (L + U)x^{(k)}\big)\), using only old values.
• Gauss–Seidel: uses each new component as soon as it is computed.

Both converge if \(A\) is strictly diagonally dominant. More generally, convergence holds if the spectral radius of the iteration matrix is less than 1.

Worked Example 3: Jacobi versus Gauss–Seidel

Solve \(4x + y = 6\), \(x + 3y = 7\) (exact solution \(x = 1, y = 2\)), starting from \((0, 0)\). The system is diagonally dominant since \(4 > 1\) and \(3 > 1\).

The update formulas are \(x = (6 - y)/4\) and \(y = (7 - x)/3\).

For example, Gauss–Seidel step 1 computes \(x = 6/4 = 1.5\), then immediately uses it: \(y = (7 - 1.5)/3 = 1.8333\). Gauss–Seidel reaches four-decimal accuracy roughly twice as fast here. For large symmetric positive definite systems, the conjugate gradient method and Krylov methods such as GMRES are the modern choices, often combined with preconditioning.

Conditioning: How Trustworthy Is the Answer?

The condition number \(\kappa(A) = \|A\|\,\|A^{-1}\|\) measures how much relative errors in \(b\) can be magnified in \(x\). A useful rule of thumb: you may lose about \(\log_{10}\kappa(A)\) significant digits. The matrix \(\begin{pmatrix}1 & 1\\ 1 & 1.0001\end{pmatrix}\) has \(\kappa_2 \approx 4 \times 10^4\), so roughly four digits of accuracy can vanish, regardless of the algorithm. A stable algorithm cannot fix an ill-conditioned problem.

Did You Know?

• Alan Turing introduced the term "condition number" and analysed rounding errors in LU decomposition in a 1948 paper.
• The QR algorithm for eigenvalues, found independently around 1961 by John Francis and Vera Kublanovskaya, was named one of the top ten algorithms of the twentieth century in a well-known list published in 2000.
• The TOP500 supercomputer ranking is based on the LINPACK benchmark, which times the solution of a huge dense system by LU decomposition.

Key Takeaways

• Solve linear systems by factorising, not by computing \(A^{-1}\).
• LU (with partial pivoting) is the default for square systems; Cholesky for symmetric positive definite ones.
• QR uses orthogonal matrices and is the stable choice for least squares.
• Iterative methods suit large sparse systems; Gauss–Seidel usually beats Jacobi, and conjugate gradient beats both for SPD matrices.
• The condition number tells you how many digits you can trust.

For step-by-step solutions to numerical analysis problems or help with MATLAB and Python assignments, contact Amathify's experts and PhD scholars on WhatsApp or through the Ask a Question form.