SC-06

Course 6

Direct Methods for Linear Systems

Learn to solve linear systems $Ax = b$ using direct methods. From Gaussian elimination to LU decomposition, understand the algorithms, their stability, and when to use them.

Learning Objectives

Implement Gaussian elimination with back substitution
Understand and apply pivoting strategies
Compute LU decomposition and use it for multiple right-hand sides
Apply Cholesky factorization for symmetric positive definite matrices
Analyze condition numbers and their effect on solution accuracy
Exploit special matrix structures for efficiency

1. Gaussian Elimination

Gaussian elimination transforms a linear system into upper triangular form, which is then solved by back substitution.

Algorithm Gaussian Elimination

Input: Matrix $A$ , vector $b$

Forward Elimination: For $k = 1, \ldots, n-1$ :

For $i = k+1, \ldots, n$ : compute multiplier $m_{ik} = a_{ik}/a_{kk}$
Update row $i$ : $a_{ij} \leftarrow a_{ij} - m_{ik} a_{kj}$ for $j = k, \ldots, n$
Update $b_i \leftarrow b_i - m_{ik} b_k$

Back Substitution: For $k = n, n-1, \ldots, 1$ :

x_k = \frac{1}{a_{kk}} \left( b_k - \sum_{j=k+1}^{n} a_{kj} x_j \right)

Theorem 6.1: Operation Count

Gaussian elimination requires:

Forward elimination: $\frac{2n^3}{3} + O(n^2)$ operations
Back substitution: $n^2$ operations

Example: Gaussian Elimination

Solve the system:

\begin{pmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 8 \\ -11 \\ -3 \end{pmatrix}

After forward elimination: $x_3 = -1, x_2 = 3, x_1 = 2$

2. Pivoting Strategies

When the pivot $a_{kk}$ is zero or very small, we need to reorder rows (and possibly columns) to maintain numerical stability.

Definition 6.1: Partial Pivoting

At step $k$ , find the row $p \geq k$ with $\max_{i \geq k} |a_{ik}|$ , then swap rows $k$ and $p$ .

This ensures $|m_{ik}| \leq 1$ for all multipliers.

Definition 6.2: Complete Pivoting

Search both rows and columns: find $\max_{i,j \geq k} |a_{ij}|$ and swap both rows and columns. More stable but rarely needed in practice.

Remark:

Scaled partial pivoting: Considers relative element sizes by comparing $|a_{ik}|/s_i$ where $s_i = \max_j |a_{ij}|$ is the row scale factor.

Note:

Partial pivoting is almost always sufficient and is the default in production software. The growth factor (maximum element during elimination divided by maximum in original matrix) is typically small with partial pivoting.

3. LU Decomposition

Definition 6.3: LU Factorization

Express $A$ as $A = LU$ where:

$L$ is lower triangular (with 1s on diagonal for Doolittle)
$U$ is upper triangular

Theorem 6.2: Existence of LU

If all leading principal minors of $A$ are nonzero, then $A$ has a unique LU factorization with unit diagonal on $L$ .

Algorithm Solving with LU

To solve $Ax = b$ given $A = LU$ :

Forward substitution: Solve $Ly = b$ for $y$
Back substitution: Solve $Ux = y$ for $x$

Each solve is $O(n^2)$ , efficient for multiple right-hand sides.

Definition 6.4: LU with Pivoting

With partial pivoting, we get $PA = LU$ where $P$ is a permutation matrix.

To solve: $Ax = b \Rightarrow PAx = Pb \Rightarrow LUx = Pb$

Example: LU Decomposition

A = \begin{pmatrix} 2 & 1 \\ 4 & 3 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 2 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 \\ 0 & 1 \end{pmatrix} = LU

4. Special Matrix Structures

Definition 6.5: Cholesky Factorization

For symmetric positive definite $A$ : $A = LL^T$ where $L$ is lower triangular with positive diagonal.

l_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2}, \quad l_{ij} = \frac{1}{l_{jj}}\left(a_{ij} - \sum_{k=1}^{j-1} l_{ik}l_{jk}\right)

Remark:

Cholesky factorization requires only $\frac{n^3}{6}$ operations (half of LU) and is numerically stable without pivoting for SPD matrices.

Definition 6.6: Tridiagonal Systems

For tridiagonal matrix $A$ with diagonals $a_i, b_i, c_i$ , theThomas algorithm solves in $O(n)$ :

\begin{pmatrix} b_1 & c_1 \\ a_2 & b_2 & c_2 \\ & \ddots & \ddots & \ddots \\ & & a_n & b_n \end{pmatrix}

Algorithm Thomas Algorithm

Forward sweep: For $i = 2, \ldots, n$ :

w = a_i / b_{i-1}, \quad b_i \leftarrow b_i - w c_{i-1}, \quad d_i \leftarrow d_i - w d_{i-1}

Back substitution: $x_n = d_n/b_n$ , then for $i = n-1, \ldots, 1$ :

x_i = (d_i - c_i x_{i+1}) / b_i

5. Error Analysis and Condition Numbers

Definition 6.7: Matrix Norms

Common matrix norms:

$\|A\|_1 = \max_j \sum_i |a_{ij}|$ (max column sum)
$\|A\|_\infty = \max_i \sum_j |a_{ij}|$ (max row sum)
$\|A\|_2 = \sqrt{\rho(A^T A)}$ (spectral norm)
$\|A\|_F = \sqrt{\sum_{i,j} |a_{ij}|^2}$ (Frobenius norm)

Definition 6.8: Condition Number

\kappa(A) = \|A\| \cdot \|A^{-1}\|

For the 2-norm: $\kappa_2(A) = \sigma_{\max} / \sigma_{\min}$ (ratio of largest to smallest singular value).

Theorem 6.3: Error Bound

If $Ax = b$ and $(A + \delta A)(x + \delta x) = b + \delta b$ , then:

\frac{\|\delta x\|}{\|x\|} \leq \kappa(A) \left( \frac{\|\delta A\|}{\|A\|} + \frac{\|\delta b\|}{\|b\|} \right)

Example: Ill-Conditioned System

The Hilbert matrix $H_{ij} = 1/(i+j-1)$ is extremely ill-conditioned:

$\kappa(H_5) \approx 4.8 \times 10^5$

$\kappa(H_{10}) \approx 1.6 \times 10^{13}$

Solving systems with such matrices loses many digits of accuracy.

Remark:

Rule of thumb: If $\kappa(A) \approx 10^k$ and you work in $d$ -digit precision, expect about $d - k$ correct digits in the solution.

Practice Quiz

Direct Linear Systems Quiz

Questions

Correct

Accuracy

Gaussian elimination without pivoting transforms

Ax = b

into:

Easy

Not attempted

The computational complexity of Gaussian elimination for an

n \times n

system is:

Easy

Not attempted

LU decomposition factors

A

as:

Easy

Not attempted

Partial pivoting in Gaussian elimination selects the pivot by:

Medium

Not attempted

For a symmetric positive definite matrix, the Cholesky factorization gives:

Medium

Not attempted

The condition number

\kappa(A)

measures:

Medium

Not attempted

A tridiagonal system can be solved in:

Medium

Not attempted

\kappa(A) = 10^{10}

and we use double precision, we can expect to lose approximately:

Hard

Not attempted

The Doolittle LU decomposition sets:

Hard

Not attempted

For the

p

-norm, which inequality relates the condition number to relative errors?

Hard

Not attempted

Frequently Asked Questions

When should I use LU vs. Gaussian elimination?

Single right-hand side: Both are equivalent in cost.

Multiple right-hand sides: Use LU. Factor once ( $O(n^3)$ ), then solve each system in $O(n^2)$ .

Why is pivoting necessary?

Without pivoting, a zero or near-zero pivot causes division by zero or severe error amplification.

Partial pivoting keeps multipliers bounded ( $|m_{ik}| \leq 1$ ), limiting error growth. It's essential for numerical stability.

When should I use Cholesky vs. LU?

Cholesky: When $A$ is symmetric positive definite. It's twice as fast, more stable, and needs no pivoting.

LU: For general (non-symmetric or indefinite) matrices.

How do I know if my matrix is ill-conditioned?

Compute $\kappa(A)$ using your linear algebra library. If $\kappa \approx 1/\epsilon_{\text{mach}}$ (about $10^{16}$ for double precision), the matrix is effectively singular. For $\kappa > 10^{10}$ , be cautious about the accuracy of results.

What if the matrix is sparse?

Use sparse matrix data structures and algorithms (sparse LU with fill-reducing orderings). Libraries like UMFPACK, SuperLU, or SuiteSparse are designed for this. Direct methods may still produce significant fill-in; iterative methods (next chapter) may be better.