Part 1. Direct Methods for Linear Systems

Last Update: 2026-10-17

This part introduces the numerical solution of linear systems through Gaussian elimination and LU factorization. The purpose is not to review linear algebra or to practice hand calculations. Instead, the goal is to understand how a mathematically simple problem \(Au=b,\) becomes a computational problem when $A$ is large, the data are represented in finite precision, and the system must be solved repeatedly.

The main topics are:

• Gaussian elimination as an algorithm;
• LU factorization;
• computational cost and storage;
• triangular solves;
• partial pivoting;
• residuals, errors, and numerical stability;
• factorization reuse for multiple right-hand sides.

Throughout this part, use dense matrices unless otherwise stated. Sparse matrices and iterative methods will be introduced later.

Suggested References

The following books offer further explanations and examples for this part. Please consult only the sections on linear systems; both texts cover many additional topics that are not required here. Neither book is required reading.

Timothy Sauer. Numerical Analysis. Recommended as the primary introductory reference. Its exposition, worked examples, and computational exercises are well suited to students encountering numerical methods for the first time.

Michael T. Heath. Scientific Computing: An Introductory Survey. Recommended as a complementary reference for the conceptual underpinnings of numerical computation, especially conditioning, stability, and computational cost. Selected sections help clarify why algorithms are designed and assessed this way; sequential reading of the whole book is unnecessary.

1.1 From a Mathematical System to a Numerical Problem

Consider the linear system \(Au=b,\qquad A\in\mathbb R^{n\times n},\quad u,b\in\mathbb R^n.\)

In elementary linear algebra, one may write \(u=A^{-1}b.\)

This expression is useful for theoretical purposes, but it is generally not an appropriate numerical algorithm.

In practical computation:

1. The inverse $A^{-1}$ is usually not formed explicitly.
2. The entries of $A$ and $b$ are represented with finite precision.
3. The cost of the method depends strongly on $n$.
4. The storage requirement may become significant.
5. Small perturbations in the input may lead to large changes in the solution.
6. A computed vector may satisfy the equations approximately without being close to the exact solution.

The central computational task is therefore:

Given $A$ and $b$, compute a reliable approximation $u_h$ to the solution of $Au=b$, while controlling computational cost and assessing the quality of the result.

A numerical method should not be evaluated only by whether it returns a vector. We also need to ask:

• How much work does it require?
• How much memory does it use?
• Is it numerically stable?
• How accurate is the computed solution?
• Can the factorization be reused if $b$ changes?
• What happens when the matrix is nearly singular?

1.2 Gaussian Elimination

Gaussian elimination is the method one has studied in the linear algebra course, sometimes refers as row reduction algorithm. It transforms the original system into an upper triangular system.

Let \(A^{(1)}=A,\qquad b^{(1)}=b.\)

At elimination step $k$, assuming $a_{kk}^{(k)}\neq 0$, define \(\ell_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}, \qquad i=k+1,\ldots,n.\)

Then update the remaining rows: \(a_{ij}^{(k+1)} = a_{ij}^{(k)} - \ell_{ik}a_{kj}^{(k)}, \qquad j=k+1,\ldots,n,\)

and \(b_i^{(k+1)} = b_i^{(k)} - \ell_{ik}b_k^{(k)}.\)

After $n-1$ elimination steps, the system has the form \(Uu=\widetilde b,\)

where $U$ is upper triangular. The solution is then computed by backward substitution: \(u_n=\frac{\widetilde b_n}{u_{nn}},\) and, for $i=n-1,\ldots,1$, \(u_i = \frac{1}{u_{ii}}\left(\widetilde b_i - \sum_{j=i+1}^n u_{ij}u_j\right).\)

Important observation

Gaussian elimination consists of two logically distinct stages:

1. Factorization or elimination: transform $A$ into triangular form;
2. Triangular solve: solve the resulting triangular system.

This distinction becomes important when the same matrix $A$ is used with several right-hand sides.

1.3 Implementation Exercise: Gaussian Elimination

Implement Gaussian elimination for a dense matrix.

Your implementation should:

1. Accept $A$ and $b$ as input.
2. Avoid explicitly forming $A^{-1}$.
3. Perform the elimination in a numerically reasonable way.
4. Use backward substitution after elimination.
5. Return the computed solution and, if possible, the triangular factor.
6. Preserve the original input when requested, so that the residual can be evaluated afterward.

Test your implementation on systems with known exact solutions. For example, generate a random matrix $A$, choose a random vector $u_{\mathrm{exact}}$, and define

Then compare the computed solution $u_h$ with $u_{\mathrm{exact}}$.

Report: \(\|u_h-u_{\mathrm{exact}}\|_2,\) and the relative error \(\frac{\|u_h-u_{\mathrm{exact}}\|_2}{\|u_{\mathrm{exact}}\|_2}.\)

Do not use only randomly generated test problems. Also construct examples that illustrate the effect of small pivots and ill-conditioning.

1.4 LU Factorization

Gaussian elimination can be interpreted as an LU factorization.

For a suitable nonsingular matrix $A$, \(A=LU,\) where

• $L$ is lower triangular with unit diagonal;
• $U$ is upper triangular.

The entries eliminated during Gaussian elimination become the subdiagonal entries of $L$.

Instead of solving $Au=b$ directly, we solve two triangular systems: \(Ly=b,\) followed by \(Uu=y.\)

The first system is solved by forward substitution; the second is solved by backward substitution.

Multiple right-hand sides

If $A$ is fixed and several systems must be solved, \(Au^{(1)}=b^{(1)},\qquad Au^{(2)}=b^{(2)},\qquad \ldots,\) then the factorization of $A$ needs to be computed only once. Each additional right-hand side requires only: \(Ly^{(j)}=b^{(j)},\qquad Uu^{(j)}=y^{(j)}.\)

This gives the practical decomposition \(\text{total cost} = \text{factorization cost} + \text{number of right-hand sides} \times \text{triangular solve cost}.\)

This distinction is important in many applications, where the system matrix remains fixed while the data vector changes.

1.5 Computational Cost

For a dense $n\times n$ matrix:

• Gaussian elimination requires approximately $O(n^3)$ arithmetic operations;
• forward or backward substitution requires $O(n^2)$ operations;
• storing a dense matrix requires $O(n^2)$ memory.

Thus, for a single right-hand side, the factorization usually dominates the cost. For many right-hand sides, the factorization cost can be amortized.

Exercise 1: Scaling with problem size

For a sequence of dense random systems, for example \(n=50,100,200,400,800,\) measure separately:

1. the factorization time;
2. the triangular solve time;
3. the total time;
4. the memory usage, if available.

Plot the measured factorization time against $n$ on a log-log scale. Discuss whether the observed behavior is consistent with $O(n^3)$ complexity.

Repeat the experiment with several right-hand sides: \(AX=B,\qquad B\in\mathbb R^{n\times m}.\)

Compare:

• solving each right-hand side independently;
• factoring $A$ once and reusing the factors.

The report should distinguish the cost of factorization from the cost of solving.

1.6 Why Pivoting Is Necessary

The basic Gaussian elimination procedure assumes that the pivot $a_{kk}^{(k)}$ is suitable for division. A nonzero pivot is mathematically sufficient, but it may be numerically unsafe if the pivot is very small.

Consider a step in which \(|a_{kk}^{(k)}|\ll |a_{ik}^{(k)}|\)

for some $i>k$. The elimination multiplier \(\ell_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}\) may then be very large. In finite-precision arithmetic, this can lead to substantial round-off error and loss of accuracy.

This motivates pivoting.

1.7 Partial Pivoting

At step $k$, partial pivoting selects the entry of largest magnitude in the current column among the remaining rows: \(p = \mathop{\arg\max}_{i=k,\ldots,n}|a_{ik}^{(k)}|.\)

If $p\neq k$, interchange rows $p$ and $k$, including the corresponding entries of the right-hand side.

Then use the new $a_{kk}^{(k)}$ as the pivot.

The resulting factorization can be written as \(PA=LU,\) where $P$ is a permutation matrix.

The solution process is therefore: \(PAu=LUu=Pb.\)

First solve \(Ly=Pb,\) then solve \(Uu=y.\)

Terminology

Partial pivoting is sometimes described as column-based pivot selection, because the pivot is selected from the current column. The actual operation is a row interchange.

It should be distinguished from complete pivoting, in which both rows and columns may be interchanged. Complete pivoting is not required in this introductory training.

1.8 Pivoting Exercise

Implement Gaussian elimination with partial pivoting.

Compare the following methods:

1. Gaussian elimination without pivoting;
2. Gaussian elimination with partial pivoting;
3. a reliable library solver, such as a standard solve routine.

For each method, report:

• relative solution error;
• relative residual;
• execution time;
• whether the method encounters a zero or very small pivot.

Use several types of test problems:

Test A: A matrix with a small initial pivot

Construct a nonsingular matrix whose first diagonal entry is very small while another entry in the first column is of moderate size.

Test B: A nearly singular matrix

Use a matrix with nearly linearly dependent columns or rows.

Test C: A random matrix with controlled conditioning

Generate matrices with different approximate condition numbers, if possible.

The purpose is not merely to show that pivoting “works.” Explain:

• what problem pivoting addresses;
• why pivoting does not remove ill-conditioning;
• why a well-designed algorithm can still produce an inaccurate solution for an ill-conditioned problem.

1.9 Residual, Error, and Accuracy

Let $u$ denote the exact solution and $u_h$ the computed solution. Two different quantities should be distinguished.

Solution error If the exact solution is available, \(e=u-u_h.\) The relative solution error is \(\frac{\|u-u_h\|}{\|u\|}.\)

Residual The residual of the computed solution is \(r=b-Au_h.\) The relative residual is commonly measured as \(\frac{\|b-Au_h\|}{\|A\|\,\|u_h\|+\|b\|},\) or, in a simpler form, \(\frac{\|b-Au_h\|}{\|b\|}.\) The residual measures how well $u_h$ satisfies the computed linear equations. It is not, by itself, the same as the error in the solution.

A small residual does not necessarily imply a small solution error. The relationship between the two depends on the conditioning of $A$.

1.10 Conditioning and Numerical Stability

The condition number of a nonsingular matrix in a chosen norm is \(\kappa(A)=\|A\|\,\|A^{-1}\|.\)

A large condition number means that the problem is sensitive to perturbations in the data. This is a property of the problem itself, not necessarily a defect of the algorithm.

Three concepts should be distinguished:

Conditioning How sensitive is the exact solution to perturbations in $A$ or $b$?

Stability How much additional error does the numerical algorithm introduce?

Accuracy How close is the computed solution to the exact solution?

A stable algorithm cannot generally produce a highly accurate answer when the underlying problem is severely ill-conditioned and the input data are already inaccurate.

For a computed solution $u_h$, a typical residual-based estimate has the form \(\frac{\|u-u_h\|}{\|u\|} \lesssim \kappa(A) \frac{\|b-Au_h\|}{\|b\|},\) under suitable assumptions and up to norm-dependent constants.

This estimate should be used as an interpretation tool rather than treated as an exact equality.

1.11 Conditioning Exercise

Construct a family of matrices $A_\varepsilon$ whose condition number increases as $\varepsilon$ decreases. For example, \(A_\varepsilon = \begin{pmatrix} 1&1\\ 1&1+\varepsilon \end{pmatrix}, \qquad \varepsilon>0.\)

For a chosen exact solution $u$, define \(b_\varepsilon=A_\varepsilon u.\)

Then perturb the right-hand side: \(\widetilde b_\varepsilon = b_\varepsilon+\delta b.\)

Solve both systems and compare:

• $\kappa(A_\varepsilon)$;
• $|b_\varepsilon-\widetilde b_\varepsilon|$;
• $|u-\widetilde u|$;
• the residual of $\widetilde u$;
• the predicted and observed error amplification.

Repeat the experiment in single and double precision if possible.

The report should explain whether the observed loss of accuracy is caused primarily by:

• the conditioning of the problem;
• the numerical implementation;
• the finite precision;
• or a combination of these effects.

1.12 A Complete Direct-Solver Experiment

Implement a complete direct solver based on LU factorization with partial pivoting.

For a family of test systems, report:

1. matrix dimension $n$;
2. number of nonzero entries, if relevant;
3. estimated or computed condition number;
4. factorization time;
5. solve time;
6. relative residual;
7. relative solution error when the exact solution is known;
8. memory usage, if available;
9. whether pivoting was required;
10. whether the factorization was reused.

Use both randomly generated systems and at least one structured family of matrices.

The final discussion should answer:

• When is a dense direct method appropriate?
• Why is explicitly forming $A^{-1}$ usually undesirable?
• When does factorization reuse become advantageous?
• What does partial pivoting improve?
• What does partial pivoting not improve?
• Why can a small residual coexist with a large solution error?
• What features of a problem may motivate the use of iterative methods instead?

1.13 Completion Criteria

Part 1 is complete when you can:

• explain Gaussian elimination as an algorithm rather than a hand-calculation procedure;
• implement or correctly use LU factorization;
• solve triangular systems;
• explain the distinction between factorization cost and solve cost;
• use partial pivoting and explain why it is needed;
• compute and interpret residuals;
• distinguish conditioning, stability, and accuracy;
• compare a computed solution with a known exact solution;
• explain why direct methods may become unsuitable for large or sparse systems.

The next part will use the residual \(r_k=b-Au_k\) as the basis for iterative correction methods.