Solution of Linear Algebraic Equations

Circuit analysis, structural frames, the Google PageRank vector and the normal equations of Unit 2 all lead to systems of linear equations Ax = b, sometimes with thousands of unknowns.

This unit covers direct methods (Gauss elimination, Gauss–Jordan, LU), iterative methods (Jacobi, Gauss–Seidel), the dangers of ill-conditioning, and the power method for eigenvalues.

Syllabus note. The mirror of the TU CACS252 outline used for this book (see bca/SOURCES.md) lists Units 1, 2, 3, 5, 6 but mentions "Linear Algebraic Equations" in the course description and lab list.

The unit is included here as Unit 4, as in other TU outlines; its 10 hours are the remainder of 45.

4.1 Gauss elimination

Forward elimination converts Ax = b to an upper-triangular system by subtracting multiples m_{ik} = \frac{a_{ik}}{a_{kk}} of the pivot row from the rows below. Back substitution then solves from the last equation upwards.

The cost is about \frac{n^{3}}{3} multiplications.

Partial pivoting. Before eliminating column k, swap the row with the largest |a_{ik}| (i \geq k) into the pivot position. This avoids dividing by a zero or tiny pivot and keeps the multipliers |m| \leq 1, which controls round-off. Always use it on a computer.

Singular and inconsistent systems. If a zero pivot remains in a column after trying all rows below, the matrix is singular: either no solution (a row 0 = nonzero) or infinitely many (a row 0 = 0).

The determinant equals the product of the pivots times (−1) to the number of swaps.

4.2 Gauss–Jordan method and matrix inversion

Gauss–Jordan continues eliminating above the pivots too, reducing A to the identity; the right-hand side becomes the solution. Applied to the augmented matrix [A | I] it produces A^{-1}.

Gauss–Jordan costs about \frac{n^{3}}{2}, more than Gauss elimination, so use it for inverses, not for single systems.

4.3 LU factorisation

If no row swaps are needed, A = L U with L unit-lower-triangular (1s on the diagonal, the multipliers below) and U upper-triangular (the result of elimination).

Then Ax = b becomes two cheap triangular solves: Ly = b (forward substitution), Ux = y (back substitution).

Once A is factorised, new right-hand sides cost only O (n^{2}), which is why LU is the workhorse of numerical libraries (with pivoting: P A = L U). The determinant is the product of U's diagonal.

Variants: Crout (U has the unit diagonal), Cholesky A = L L^{T} for symmetric positive-definite matrices (half the work, no pivoting needed).

4.4 Iterative methods: Jacobi and Gauss–Seidel

For large sparse systems direct methods fill in zeros and cost too much. Iterative methods start from a guess and improve it. Solve equation \ast i \ast for x_{i}:

• Jacobi: use only the old iterate x^{(k)} on the right (all components updated "simultaneously").
• Gauss–Seidel: use the newest values as soon as they are available (components updated one after another), so it normally converges about twice as fast.

Convergence. Both converge for any start if A is strictly diagonally dominant (|a_{ii}| > \sum_{j \neq i} |a_{ij}| in every row).

That is sufficient, not necessary; the exact condition is that the iteration matrix has spectral radius below 1. Rearranging the equations so that large coefficients lie on the diagonal often helps.

4.5 Ill-conditioning

A system is ill-conditioned if tiny changes in A or b cause large changes in x. The culprit is a nearly singular matrix, and no algorithm can repair it; pivoting only avoids additional damage.

The condition number κ(A) = ‖A‖·‖A⁻¹‖ measures the amplification: relative error in x \leq \kappa (A) × relative error in the data.

For the matrix above \kappa = 40002, so a relative data error of 10^{-4} may produce a relative solution error of order one.

The 3 \times 3 Hilbert matrix, with entries \frac{1}{i + j - 1}, has \kappa = 524; it is moderately conditioned for n = 3 but becomes terrible as n grows. Rules of thumb: \kappa \approx 10^{k} loses about k significant digits.

Remedies: scaling, pivoting, higher precision, iterative refinement, or reformulating the problem.

4.6 Eigenvalues: the power method

The power method finds the dominant eigenvalue (largest in magnitude) \lambda_{1} and its eigenvector. Start with any v_{0} (with a component along the dominant eigenvector).

Repeat: w = A v_{k-1} ;\ \lambda_{k} = the entry of w with the largest magnitude; v_{k} = \frac{w}{\lambda_{k}}. Then \lambda_{k} \to \lambda_{1} and v_{k} → eigenvector, at a rate |\frac{\lambda_{2}}{\lambda_{1}}|.

4.7 Seeing the ideas: elimination, iteration speed and the geometry of ill-conditioning

Carl Friedrich Gauss gave elimination its name, and his portrait opens Plate 4.1. The plate also shows why some systems are dangerous: when two lines are nearly parallel, a tiny change moves their intersection a long way.

(a) Gauss

(b) Well and badly conditioned

4.7.1 Eight uses of linear systems and their numerical solution

1. Electrical networks. Kirchhoff's laws turn a circuit with many loops into a set of linear equations. A hydropower grid has thousands of nodes and is solved by computer. Elimination or iteration gives all currents at once.
2. Structures and bridges. Engineers divide a truss or bridge into joints and write force-balance equations for each. The unknown forces in every member come from one big linear system. The same method checks buildings against earthquake loads.
3. Economics and planning. Input-output tables describe how each industry uses the products of others. Solving a linear system shows how much each sector must produce to meet final demand. Governments use it in planning.
4. Curve fitting. The normal equations of least squares are linear equations in the unknown coefficients. Every trend line and calibration uses them. A spreadsheet solves them in a moment.
5. Computer graphics. Moving, rotating and projecting a 3D model involves solving linear systems for the transformation. Games solve millions of such equations per second. Graphics cards are built for it.
6. Search and ranking. The ranking of web pages is an eigenvector problem solved by repeating a matrix multiplication, a method close to the power method. Recommendation systems use the same idea. Linear algebra runs behind the screen.
7. Weather and flow simulation. Dividing the air or a river into small cells gives millions of linked equations, solved by fast iterative methods. Direct elimination would be far too slow. Choosing the right method decides whether a forecast arrives in time.
8. Network and traffic analysis. Flow of traffic, water in a pipe network and data in a computer network follow balance equations. Solving them locates bottlenecks. Planners in Kathmandu can test a road change on paper first.

(a) A slide rule: three digits of precision

(b) Digits lost = log10 of the condition number

Formula and fact sheet

• Multiplier m_{ik} = \frac{a_{ik}}{a_{kk}}; partial pivoting: swap to the largest |a_{ik}| ;\ \det = \pm \Pi pivots.
• Gauss–Jordan on [A | I] gives A^{-1}. LU: Ly = b then Ux = y.
• Jacobi: old values; Gauss–Seidel: newest values; converge if diagonally dominant.
• κ(A) = ‖A‖‖A⁻¹‖; digits lost ≈ \log_{10} \kappa.
• Power method: w = Av, λ = largest entry, v = \frac{w}{\lambda}; rate |\frac{\lambda_{2}}{\lambda_{1}}|.

Practice

1. Solve 2 x + y - z = 8 ,\ - 3 x - y + 2 z = - 11 ,\ - 2 x + y + 2 z = - 3 by Gauss elimination with partial pivoting.
2. Find the inverse of \begin{bmatrix} 2 & 1 \\ 5 & 3 \end{bmatrix} by Gauss–Jordan.
3. Factorise A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix} as LU and solve Ax = (4, 10, 24).
4. Do two Jacobi and two Gauss–Seidel iterations for the 4 \times 4 system of Section 4.4 from zero; which gives the better x_{2}?
5. State the sufficient condition for convergence of Jacobi and Gauss–Seidel. Is it necessary?
6. Why does Gauss–Seidel diverge for x + 3 y = 5 ,\ 4 x + y = 6 but converge after swapping the equations?
7. What is an ill-conditioned system? Illustrate with an example and name the condition number.
8. Apply the power method to A = \begin{bmatrix} 2 & - 1 & 0 \\ - 1 & 2 & - 1 \\ 0 & - 1 & 2 \end{bmatrix} and report the eigenvalue after four iterations.
9. Compare Gauss elimination and LU factorisation when there are many right-hand sides.

Model answers: medium (5 marks) and long (10 marks) questions

Q1 (5 marks). Solve 10 x + y + z = 12 ,\ 2 x + 10 y + z = 13 ,\ 2 x + 2 y + 10 z = 14 by Gauss elimination.

Augmented matrix. [10 1 1 | 12; 2 10 1 | 13; 2 2 10 | 14].

*Step 1 (pivot 10 ) . \ast m21 = \frac{2}{10} = 0.2 and m31 = 0.2: R2 = R2 - 0.2R1 gives [0 9.8 0.8 | 10.6]; R3 = R3 - 0.2R1 gives [0 1.8 9.8 | 11.6].

*Step 2 (pivot 9.8 ) . \ast m32 = 1.8 / 9.8 = 0.18367: R3 = R3 - 0.18367R2 gives [0 0 9.6531 | 9.6531] (Figure 4.1).

Back substitution. z = 9.6531 / 9.6531 = 1; y = \frac{10.6 - 0.8 x 1}{9.8} = 1; x = \frac{12 - 1 - 1}{10} = 1.

Check 10 + 1 + 1 = 12 ,\ 2 + 10 + 1 = 13 ,\ 2 + 2 + 10 = 14.

Note with a zero or very small pivot, rows must be swapped (partial pivoting) to keep rounding errors small. The work is about \frac{n^{3}}{3} multiplications.

Conclusion the solution is x = y = z = 1.

Q2 (5 marks). Solve the same system by the Jacobi and Gauss-Seidel methods (four iterations from zero). State the condition for convergence.

Formulas. x = \frac{12 - y - z}{10} ,\ y = \frac{13 - 2 x - z}{10} ,\ z = \frac{14 - 2 x - 2 y}{10}. Jacobi uses only old values of an iteration; Gauss-Seidel uses each new value immediately.

Gauss-Seidel is already correct to four decimals after four sweeps, while Jacobi still has an error of 0.008 (Figure 4.2).

Convergence guaranteed when the matrix is strictly diagonally dominant, that is |a (i i)| > sum of |a (i j)| for j not equal i in every row; here 10 > 2 ,\ 10 > 3 and 10 > 4.

Conclusion iterative methods suit large sparse systems and Gauss-Seidel is usually about twice as fast.

Q3 (10 marks). Compare direct and iterative methods for linear systems. Explain LU factorisation and ill-conditioning with an example.

Introduction. To solve Ax = b, direct methods (Gauss elimination, Gauss-Jordan, LU) give the answer after a fixed number of operations, while iterative methods (Jacobi, Gauss-Seidel) improve a guess until it is accurate enough.

Direct methods. Gauss elimination reduces A to upper-triangular form and back-substitutes, with cost about \frac{n^{3}}{3}.

LU factorisation writes A = L U with L lower triangular (holding the multipliers, here 0.2, 0.2, 0.18367) and U upper triangular (the reduced matrix).

To solve for any b, first solve Ly = b by forward substitution, then Ux = y by back substitution, each costing only n^{2}. LU is therefore ideal when the same matrix is used with many right-hand sides.

Gauss-Jordan continues to a diagonal form and can find an inverse but costs more.

Iterative methods. They store only the nonzero entries and need diagonal dominance (Figure 4.2). For a system with a million unknowns, such as a weather grid, they are the only practical choice.

Ill-conditioning. A system is ill-conditioned when small changes in the data cause large changes in the solution (Plate 4.1b).

Example x + y = 2 and x + 1.0001y = 2.0001 has the solution x = 1 ,\ y = 1. Changing the right side to 2.0002 gives y = 0.0002 / 0.0001 = 2 and x = 0: a change of 0.0001 in one datum changed the answer from (1, 1) to (0, 2).

The condition number is \mid |A| \mid x \mid |A^{- 1}| \mid = 2.0001 x 20,001, about 40,000, so relative errors in the data can be magnified up to 40,000 times.

Conclusion. Use LU with partial pivoting for ordinary systems, iterative methods for huge sparse ones, and always check the condition number before trusting an answer.