Systems of linear equations · Lecture 13

LU, Cholesky and Tridiagonal Systems

Elimination done once and reused: LU factorisation (Doolittle and Crout) for many right-hand sides, Cholesky for symmetric positive definite matrices at half the cost, and the Thomas algorithm that solves tridiagonal systems in linear time.

Md Mahbub E Noor · Assistant Professor
CSE-3107 Numerical Methods · B.Sc. 3rd Year 1st Semester
Department of CSE, University of Barishal

Why this lecture

The problem: same plate, a new temperature every minute

The top edge of the plate from Lecture 12 is heated and cooled during the day. The matrix $A$ never changes; only the right-hand side $\mathbf{b}$ does.

One reading per minute for a day: 1440 right-hand sides, one matrix
Gaussian elimination repeats the $\tfrac23 n^3$ elimination every time, although $A$ is the same
A realistic model with $n = 1000$ unknowns: $1440 \times 0.67$ s = 16 minutes
Can we do the expensive part once and reuse it?

top edge	$T_1$	$T_2$	$T_3$	$T_4$
100°	43.75	56.25	18.75	31.25
80°	?	?	?	?
60°	?	?	?	?

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

2 / 25

Today's goals

Explain why Gaussian elimination is a factorisation $A = LU$
Compute $L$ and $U$ by Doolittle and Crout, and solve with forward and back substitution
Recognise symmetric positive definite matrices and factor them by Cholesky, $A = LL^T$
Solve tridiagonal systems with the Thomas algorithm in $O(n)$ operations
Choose the cheapest safe solver for a given matrix

Last time

Gaussian elimination with partial pivoting: $\tfrac23 n^3$ for each system.

Today

Factor once, solve many times; exploit symmetry and zeros.

Systems of linear equations · you are here

L11 · Linear Systems, Cramer's Rule and the Inverse→L12 · Gaussian and Gauss–Jordan Elimination→L13 · LU, Cholesky and Tridiagonal Systems→L14 · Iterative Methods

Builds on: L12

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

3 / 25

Start here

The idea in one picture

A full matrix is hard to solve with.Split it once into a lower triangle L times an upper triangle U.Now take a right-hand side b.Triangle 1: solve L y = b from the top down (easy).Triangle 2: solve U x = y from the bottom up (easy).Another b? Skip the splitting: L and U are reused.

In plain words

Elimination does the same work every time, even if only b changes
Do the hard part once: split A into two triangles L and U
Every new b then needs only two quick triangular solves

New words today

factorisation: writing a matrix as a product of simpler ones
L, U: lower and upper triangular factors
forward substitution: solving from the top down
positive definite (SPD): symmetric, and xᵀAx > 0 for every x ≠ 0
tridiagonal: non-zeros only on three diagonals

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

4 / 25

The idea: elimination is a factorisation

Keep the multipliers in L; the eliminated matrix is U

\underbrace{\begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9\end{pmatrix}}_{A} = \underbrace{\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1\end{pmatrix}}_{L \text{ (multipliers)}}\underbrace{\begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2\end{pmatrix}}_{U}

Then Ax = b splits into two triangular solves

L\mathbf{y} = \mathbf{b} \;\;(\text{forward}), \qquad U\mathbf{x} = \mathbf{y} \;\;(\text{backward})

$l_{ik} = m_{ik}$ : exactly the multipliers of Lecture 12
Factor once: $\tfrac23 n^3$ operations
Each new $\mathbf{b}$ : two triangular solves, $2n^2$ operations
$n = 1000$ : 2 ms per right-hand side instead of 0.67 s

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

5 / 25

The method, step by step

Doolittle's LU method, step by step

Doolittle: $L$ has ones on its diagonal. Compute row $i$ of $U$ , then column $i$ of $L$ , for $i = 1, 2, \dots, n$ .

Row $i$ of $U$ : $\;u_{ij} = a_{ij} - \sum_{k<i} l_{ik}u_{kj}$ for $j \ge i$
Column $i$ of $L$ : $\;l_{ji} = \big(a_{ji} - \sum_{k<i} l_{jk}u_{ki}\big)/u_{ii}$ for $j > i$ (needs $u_{ii} \ne 0$ )
Forward: $y_i = b_i - \sum_{k<i} l_{ik}y_k$
Backward: $x_i = \big(y_i - \sum_{k>i} u_{ik}x_k\big)/u_{ii}$
For another $\mathbf{b}$ , repeat only steps 3 and 4

Pseudocode

for i = 1 … n:              # factor once
    for j = i … n:
        u[i][j] ← a[i][j] − Σ_{k<i} l[i][k]·u[k][j]
    l[i][i] ← 1
    for j = i+1 … n:
        l[j][i] ← (a[j][i] − Σ_{k<i} l[j][k]·u[k][i]) / u[i][i]
for each right-hand side b:   # solve often
    y ← forward(L, b)
    x ← backward(U, y)

Same operations as elimination, organised so that each entry of $L$ and $U$ is computed once.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

6 / 25

Problem 1: Doolittle LU, then solve

Factor

A = \begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9\end{pmatrix}

as

LU

and solve

A\mathbf{x} = (3, 7, 19)^T

.

1

Row 1 of

U

:

(2, 1, 1)

. Column 1 of

L

:

l_{21} = 4/2 = 2

,

l_{31} = 8/2 = 4

2

u_{22} = 3 - 2(1) = 1

,

u_{23} = 3 - 2(1) = 1

;

\;l_{32} = (7 - 4(1))/1 = 3

;

\;u_{33} = 9 - 4(1) - 3(1) = 2

3

Forward

L\mathbf{y} = \mathbf{b}

:

y_1 = 3

,

y_2 = 7 - 2(3) = 1

,

y_3 = 19 - 4(3) - 3(1) = 4

4

Backward

U\mathbf{x} = \mathbf{y}

:

x_3 = 4/2 = 2

,

x_2 = 1 - 2 = -1

,

x_1 = (3 - (-1) - 2)/2 = 1

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

,

U = \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2\end{pmatrix}

,

\mathbf{x} = (1, -1, 2)

.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

7 / 25

Doolittle or Crout? Where the ones go

Ones on the diagonal of L

Doolittle: $l_{ii} = 1$

$U$ = the matrix left by Gaussian elimination
Rows of $U$ first, then columns of $L$

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

Ones on the diagonal of U

Crout: $u_{ii} = 1$

$l_{ij} = a_{ij} - \sum_{k<j} l_{ik}u_{kj}$ (column of $L$ first)
$u_{ij} = (a_{ij} - \sum_{k<i} l_{ik}u_{kj})/l_{ii}$ (then row of $U$ )

\begin{pmatrix} 2 & 0 & 0 \\ 4 & 1 & 0 \\ 8 & 3 & 2\end{pmatrix}\begin{pmatrix} 1 & 0.5 & 0.5 \\ 0 & 1 & 1 \\ 0 & 0 & 1\end{pmatrix}

Both cost the same and solve the same way. Crout's $L$ is Doolittle's $L$ with each column multiplied by the pivot; Crout's $U$ is Doolittle's $U$ with each row divided by it. A problem that names a method expects that placement of the ones.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

8 / 25

Symmetric positive definite matrices and Cholesky

Symmetric positive definite (SPD)

A = A^T \quad\text{and}\quad \mathbf{x}^T A\mathbf{x} > 0 \;\text{ for every } \mathbf{x} \ne \mathbf{0}

Then A has a "square root" factorisation

A = LL^T: \quad l_{jj} = \sqrt{a_{jj} - \sum_{k<j} l_{jk}^2}, \quad l_{ij} = \frac{a_{ij} - \sum_{k<j} l_{ik}l_{jk}}{l_{jj}}

Where SPD matrices come from: energies, covariance matrices, least squares $A^TA$ (Lecture 23), the plate and circuit matrices
Quick test: symmetric, positive diagonal, all pivots (or all $l_{jj}^2$ ) positive
Only one triangle is computed: $\tfrac13 n^3$ operations, half of LU
No pivoting needed: the square roots are of positive numbers and stay well-behaved

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

9 / 25

The method, step by step

Cholesky, step by step

Check $A = A^T$
For column $j = 1, \dots, n$ : $\;l_{jj} = \sqrt{a_{jj} - \sum_{k<j} l_{jk}^2}$ ; if the number under the root is $\le 0$ , stop: not SPD
Below it: $l_{ij} = \big(a_{ij} - \sum_{k<j} l_{ik}l_{jk}\big)/l_{jj}$ for $i > j$
Forward: $L\mathbf{y} = \mathbf{b}$
Backward: $L^T\mathbf{x} = \mathbf{y}$

Pseudocode

for j = 1 … n:
    s ← a[j][j] − Σ_{k<j} l[j][k]²
    if s ≤ 0: stop (not SPD)
    l[j][j] ← √s
    for i = j+1 … n:
        l[i][j] ← (a[i][j] − Σ_{k<j} l[i][k]·l[j][k]) / l[j][j]
y ← forward(L, b)
x ← backward(Lᵀ, y)

$n$ square roots, about $\tfrac13 n^3$ other operations.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

10 / 25

Problem 2: Cholesky factorisation

Factor

A = \begin{pmatrix} 4 & 2 & -2 \\ 2 & 10 & 2 \\ -2 & 2 & 6\end{pmatrix}

as

LL^T

and solve

A\mathbf{x} = (-2, -4, 8)^T

.

1

A

is symmetric. Column 1:

l_{11} = \sqrt4 = 2

,

l_{21} = 2/2 = 1

,

l_{31} = -2/2 = -1

2

Column 2:

l_{22} = \sqrt{10 - 1^2} = 3

,

l_{32} = (2 - (-1)(1))/3 = 1

;

\;

column 3:

l_{33} = \sqrt{6 - 1 - 1} = 2

3

Forward

L\mathbf{y} = \mathbf{b}

:

y_1 = -2/2 = -1

,

y_2 = (-4 - 1(-1))/3 = -1

,

y_3 = (8 - (-1)(-1) - 1(-1))/2 = 4

4

Backward

L^T\mathbf{x} = \mathbf{y}

:

x_3 = 4/2 = 2

,

x_2 = (-1 - 1(2))/3 = -1

,

x_1 = (-1 - 1(-1) - (-1)(2))/2 = 1

\mathbf{x} = (1, -1, 2)

. Every square root was of a positive number, so

A

is SPD.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

11 / 25

Tridiagonal systems

Only three diagonals are nonzero

\begin{pmatrix} b_1 & c_1 & & \\ a_2 & b_2 & c_2 & \\ & \ddots & \ddots & \ddots \\ & & a_n & b_n\end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n\end{pmatrix} = \begin{pmatrix} d_1 \\ d_2 \\ \vdots \\ d_n\end{pmatrix}

Elimination touches only one entry per row (Thomas algorithm)

\begin{aligned} c'_1 &= \frac{c_1}{b_1}, & d'_1 &= \frac{d_1}{b_1} \\ c'_i &= \frac{c_i}{b_i - a_i c'_{i-1}}, & d'_i &= \frac{d_i - a_i d'_{i-1}}{b_i - a_i c'_{i-1}} \\ x_n &= d'_n, & x_i &= d'_i - c'_i\,x_{i+1} \end{aligned}

Where they come from: heat in a rod, cubic splines (Lecture 18), implicit ODE methods
Store three vectors, not $n^2$ numbers
Cost $8n - 7$ operations instead of $\tfrac23 n^3$
Safe without pivoting when $|b_i| > |a_i| + |c_i|$ (diagonally dominant)

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

12 / 25

The method, step by step

The Thomas algorithm, step by step

Store the sub-, main and super-diagonals $a_i$ , $b_i$ , $c_i$ and the right side $d_i$
Start: $c'_1 = c_1/b_1$ , $d'_1 = d_1/b_1$
Sweep down $i = 2, \dots, n$ : $\;w = b_i - a_i c'_{i-1}$ , $\;c'_i = c_i/w$ , $\;d'_i = (d_i - a_i d'_{i-1})/w$
Sweep up: $x_n = d'_n$ , then $x_i = d'_i - c'_i x_{i+1}$ for $i = n-1, \dots, 1$
Check one equation

Pseudocode

input a, b, c, d (length n)
c′[1] ← c[1]/b[1];  d′[1] ← d[1]/b[1]
for i = 2 … n:            # down
    w ← b[i] − a[i]·c′[i−1]
    if w = 0: stop
    c′[i] ← c[i]/w
    d′[i] ← (d[i] − a[i]·d′[i−1])/w
x[n] ← d′[n]
for i = n−1 … 1:          # up
    x[i] ← d′[i] − c′[i]·x[i+1]

Two loops, no nested loop: the work grows like $n$ .

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

13 / 25

Problem 3: heat in a rod by the Thomas algorithm

A rod has ends at 100° and 0°; four interior points satisfy

-T_{i-1} + 2T_i - T_{i+1} = 0

. So

2T_1 - T_2 = 100

,

-T_1 + 2T_2 - T_3 = 0

,

-T_2 + 2T_3 - T_4 = 0

,

-T_3 + 2T_4 = 0

.

$i$	$a_i$	$b_i$	$c_i$	$d_i$	$w = b_i - a_ic'_{i-1}$	$c'_i$	$d'_i$	$T_i$
1	—	2	−1	100	2	−0.5	50	80
2	−1	2	−1	0	1.5	−0.666667	33.333333	60
3	−1	2	−1	0	1.333333	−0.75	25	40
4	−1	2	—	0	1.25	—	20	20

Down sweep fills $c'$ and $d'$ ; up sweep: $T_4 = 20$ , $T_3 = 25 + 0.75(20) = 40$ , $T_2 = 33.33 + 0.6667(40) = 60$ , $T_1 = 50 + 0.5(60) = 80$ . A straight line from 100° to 0°, as physics predicts.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

14 / 25

Classical vs numerical

Same problems, two approaches

Formula

Classical (inverse matrix)

Compute $A^{-1}$ once ( $2n^3$ ), then $\mathbf{x} = A^{-1}\mathbf{b}$ ( $2n^2$ ) per reading
Plate, 3 readings: same temperatures
Tridiagonal $A$ : $A^{-1}$ is full, so the zeros are lost

vs

Algorithm

Numerical (LU, Thomas)

Factor once ( $\tfrac23 n^3$ ), then $2n^2$ per reading: three times cheaper, more accurate
Plate, 3 readings: one factorisation, three pairs of substitutions
Tridiagonal $A$ : Thomas keeps the zeros, $8n$ operations

When the classical approach cannot work

A rod with a million grid points

A full matrix of $10^6 \times 10^6$ doubles needs 8 TB of memory, and dense elimination $6.7\times10^{17}$ operations (21 years at $10^9$ per second). The Thomas algorithm: $8\times10^6$ operations (8 ms) and four vectors (32 MB).

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

15 / 25

LU once, many right-hand sides

Doolittle factorisation of the plate matrix ( $4\times4$ )
Then three right-hand sides: top edge 100°, 80°, 60°
Each solve is a forward and a backward substitution
The factorisation loop runs once

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

16 / 25

The Thomas algorithm in Python and C

Three diagonals and $d$ stored as vectors
Down sweep: $c'$ and $d'$ ; up sweep: $x$
Change $n$ to 9 or 99: the code is unchanged
Lab: add a heat source (nonzero $d_i$ )

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

17 / 25

When the factorisations break

Limitations shown failing

LU without pivoting: $\begin{pmatrix} 0 & 1 \\ 1 & 1\end{pmatrix}$ has no $LU$ at all ( $u_{11} = 0$ ); with row swaps $PA = LU$ works
Cholesky on $\begin{pmatrix} 1 & 2 \\ 2 & 1\end{pmatrix}$ : $l_{22} = \sqrt{1 - 4}$ , not SPD
Thomas without diagonal dominance can meet $w = 0$ , just like a zero pivot
None of them helps when $A$ is ill-conditioned: the factors inherit $\kappa(A)$

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

18 / 25

Practice

Practice problem

Problem 4"Factorise

A = \begin{pmatrix} 1 & 1 & 1 \\ 4 & 3 & -1 \\ 3 & 5 & 3\end{pmatrix}

as

LU

by Doolittle's method and solve

A\mathbf{x} = (1, 6, 4)^T

."

How to write a complete answer

State the form: $L$ unit lower triangular, $U$ upper triangular
Row 1 of $U$ = $(1, 1, 1)$ ; $l_{21} = 4$ , $l_{31} = 3$ ; then $u_{22} = -1$ , $u_{23} = -5$ , $l_{32} = -2$ , $u_{33} = -10$
Forward: $\mathbf{y} = (1, 2, 5)$ ; backward: $x_3 = -\tfrac12$ , $x_2 = \tfrac12$ , $x_1 = 1$
Verify $LU = A$ (one row is enough) and substitute $\mathbf{x}$ into one equation

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

19 / 25

When to use it

LU, Cholesky and Thomas: method profile

Best case

LU: one matrix, many right-hand sides (time steps, load cases, Newton with a frozen Jacobian)
Cholesky: SPD matrices, half the work and no pivoting
Thomas: tridiagonal systems in $O(n)$ , millions of unknowns in milliseconds

Worst case

LU without pivoting: zero or tiny $u_{ii}$
Cholesky on a matrix that is not SPD: a negative number under the root
Sparse matrices that are not tridiagonal: LU fills the zeros in ("fill-in")

Best suited for

Implicit time stepping (the same $A$ every step)
Statistics and ML: covariance matrices, least squares, Gaussian processes (Cholesky)
Splines and 1-D heat or diffusion problems (Thomas)

Avoid when

Just one system (plain elimination is equally good)
Huge 2-D or 3-D grids: fill-in makes LU too big; use iterative methods (Lecture 14)
Forming $A^{-1}$ "to save time": LU does the same job better

Rule of thumb: tridiagonal → Thomas; SPD → Cholesky; many right-hand sides → LU once; otherwise Gaussian elimination with pivoting.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

20 / 25

Comparative analysis: direct solvers

Method	Needs	Factor	Each new $\mathbf{b}$	Memory	$n = 1000$ , 1440 b's
Gauss each time	any nonsingular $A$	—	$\tfrac23 n^3$	$n^2$	16 min
Inverse, then multiply	any nonsingular $A$	$2n^3$	$2n^2$	$n^2$	4.9 s
LU (with pivoting)	any nonsingular $A$	$\tfrac23 n^3$	$2n^2$	$n^2$	3.5 s
Cholesky	SPD	$\tfrac13 n^3$	$2n^2$	$\tfrac12 n^2$	3.2 s
Thomas	tridiagonal, diag. dominant	$\approx 3n$	$\approx 5n$	$4n$	7 ms (if tridiagonal)

Times at $10^9$ operations per second. The inverse row is slower and less accurate than LU for the same job.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

21 / 25

Roadmap

Limitations, and where we fix them

Direct factorisations are fast for dense, banded and tridiagonal matrices. Some problems are still out of reach:

What it does well (keep this)

Factor once, solve many times
Cholesky for SPD: half the work, no pivoting
Thomas: $O(n)$ for tridiagonal systems
Exact up to rounding, in a known number of steps

LimitationHow it is overcomeWhere

2-D and 3-D grids: fill-in, memory grows like

n^{1.5}

or worse

→

Iterative methods: only the nonzeros are stored

Lecture 14 →

Splines need a tridiagonal solve

→

Thomas, used directly

Lecture 18 →

Least squares gives SPD normal equations

A^TA\mathbf{x} = A^T\mathbf{y}

→

Cholesky (and QR for better conditioning)

Lecture 23 →

Implicit ODE methods solve a system every step

→

Factor once, reuse each step

Lecture 22 →

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

22 / 25

Applications

Where this method runs in the real world

Click a card to see where the factorisation happens in a real computation.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

23 / 25

Checkpoint

Checkpoint · think, then discuss with your neighbourDoolittle LU of

A = \begin{pmatrix} 2 & 1 \\ 6 & 8\end{pmatrix}

. What are

l_{21}

and

u_{22}

?

A

l_{21} = 3

,

u_{22} = 5

B

l_{21} = 3

,

u_{22} = 8

C

l_{21} = \tfrac13

,

u_{22} = 5

D

l_{21} = 6

,

u_{22} = 2

A.

l_{21} = a_{21}/u_{11} = 6/2 = 3

;

u_{22} = a_{22} - l_{21}u_{12} = 8 - 3(1) = 5

. Check:

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

.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

24 / 25

Summary

Key points

Gaussian elimination is $A = LU$ : $L$ holds the multipliers, $U$ the eliminated matrix
Factor once ( $\tfrac23 n^3$ ), then $L\mathbf{y} = \mathbf{b}$ , $U\mathbf{x} = \mathbf{y}$ ( $2n^2$ ) per right-hand side
SPD matrices: Cholesky $A = LL^T$ , half the cost, no pivoting, and a test for positive definiteness
Tridiagonal: the Thomas algorithm, $8n - 7$ operations

Opening problem: solved

One LU of the plate matrix, then three solves: top 100° gives $(43.75, 56.25, 18.75, 31.25)$ ; 80° gives $(36.25, 48.75, 16.25, 28.75)$ ; 60° gives $(28.75, 41.25, 13.75, 26.25)$ . For $n = 1000$ and 1440 readings: 3.5 s instead of 16 minutes.

Homework

Factor $\begin{pmatrix} 3 & 2 & 1 \\ 6 & 6 & 3 \\ 9 & 10 & 6\end{pmatrix}$ by Doolittle's method and solve with $\mathbf{b} = (6, 15, 25)$ . Then repeat with Crout's method.
Find the Cholesky factor of $\begin{pmatrix} 4 & -2 & 0 \\ -2 & 5 & -2 \\ 0 & -2 & 5\end{pmatrix}$ .
Solve $2x_1 - x_2 = 1$ , $-x_{i-1} + 2x_i - x_{i+1} = 0$ ( $i = 2, 3$ ), $-x_3 + 2x_4 = 1$ by the Thomas algorithm.
A model has $n = 2000$ and 500 right-hand sides. Compare the time of Gaussian elimination for each with one LU factorisation.

Next lecture

Iterative methods: Jacobi, Gauss–Seidel and SOR. When the matrix is huge and sparse, stop factoring and start iterating, just like fixed-point iteration in Lecture 8.

Md Mahbub E Noor

Lecture 13: LU, Cholesky and Tridiagonal Systems

University of Barishal

25 / 25

LU, Cholesky and Tridiagonal Systems

The problem: same plate, a new temperature every minute

Today's goals

The idea in one picture

The idea: elimination is a factorisation

Doolittle's LU method, step by step

Problem 1: Doolittle LU, then solve

Doolittle or Crout? Where the ones go

Doolittle: $l_{ii} = 1$

Crout: $u_{ii} = 1$

Symmetric positive definite matrices and Cholesky

Cholesky, step by step

Problem 2: Cholesky factorisation

Tridiagonal systems

The Thomas algorithm, step by step

Problem 3: heat in a rod by the Thomas algorithm

Same problems, two approaches

Classical (inverse matrix)

Numerical (LU, Thomas)

A rod with a million grid points

LU once, many right-hand sides

The Thomas algorithm in Python and C

Limitations shown failing

Practice problem

LU, Cholesky and Thomas: method profile

Comparative analysis: direct solvers

Limitations, and where we fix them

Where this method runs in the real world

🔥 Heat flow in time (implicit method)

🎲 Correlated random numbers (Monte Carlo, ML)

〰️ Smooth curves through points (splines)

🏗️ Bridges and buildings: many load cases

Checkpoint

Summary

LU, Cholesky and Tridiagonal Systems

The problem: same plate, a new temperature every minute

Today's goals

The idea in one picture

The idea: elimination is a factorisation

Doolittle's LU method, step by step

Problem 1: Doolittle LU, then solve

Doolittle or Crout? Where the ones go

Doolittle: lii=1l_{ii} = 1lii​=1

Crout: uii=1u_{ii} = 1uii​=1

Symmetric positive definite matrices and Cholesky

Cholesky, step by step

Problem 2: Cholesky factorisation

Tridiagonal systems

The Thomas algorithm, step by step

Problem 3: heat in a rod by the Thomas algorithm

Same problems, two approaches

Classical (inverse matrix)

Numerical (LU, Thomas)

A rod with a million grid points

LU once, many right-hand sides

The Thomas algorithm in Python and C

Limitations shown failing

Practice problem

LU, Cholesky and Thomas: method profile

Comparative analysis: direct solvers

Limitations, and where we fix them

Where this method runs in the real world

🔥 Heat flow in time (implicit method)

🎲 Correlated random numbers (Monte Carlo, ML)

〰️ Smooth curves through points (splines)

🏗️ Bridges and buildings: many load cases

Checkpoint

Summary

Doolittle: $l_{ii} = 1$

Crout: $u_{ii} = 1$