Note on Relation to Previous Math Courses
- Some parts in this topic are also covered in your previous math courses
- But math courses focus on solving linear systems by hands (i.e., analytical methods)
- But here we focus on solving them using computers (i.e., numerical methods)
- So a few steps of some methods may be different from what you might have learned
- Because what is best for solving linear systems by hands may not be the same as what is best for solving them with computers
- There are new materials, too.
Analytical vs Numerical Methods
- No concrete definitions for them
- Roughly, suppose all data (input) are arbitrary, analytical methods can still be used to get closed-form solutions exactly
- Numerical methods work only when all data (input) are concrete numbers as the values of the data do affect the sequence of operations of the methods
- Many problems cannot be solved analytically while many other problems can be solved analytically but too slowly to be practical
Block Matrix Notation
Partitioning matrices into smaller submatrices
- Notation for partitioning a matrix along rows and columns into smaller matrices
- E.g., the following 3-by-5 matrix
1 & 2 & 3 & 4 & 5 \\
6 & 7 & 8 & 9 & 10 \\
11 & 12 & 13 & 14 & 15
\end{bmatrix} = \begin{bmatrix}
A & B \\
C & D
\end{bmatrix}$$
where:
- $A = \begin{bmatrix} 1 & 2 \\ 6 & 7 \end{bmatrix}$
- $B = \begin{bmatrix} 3 & 4 & 5 \\ 8 & 9 & 10 \end{bmatrix}$
- $C = \begin{bmatrix} 11 & 12 \end{bmatrix}$
- $D = \begin{bmatrix} 13 & 14 & 15 \end{bmatrix}$
### Important Notes about Block Matrices
- Note that the dimensions of $A$, $B$, $C$, $D$ are compatible
- I.e., $A$ and $B$ have the same number of rows, $B$ and $D$ have the same number of columns, and so on
- Also, $A$, $B$, $C$, $D$ are called **submatrices**
- Can partition to as many rows/columns as you want. E.g.
$$\begin{bmatrix}
A & B \\
C & D \\
E & F
\end{bmatrix}, \quad \begin{bmatrix}
A & B & C & D \\
E & F & G & H
\end{bmatrix}$$
## Operations on Block Matrices
- See [[02 Operations with Matrices]]
- Adding/multiplying two block matrices are just like adding/multiplying two normal matrices *(assuming the submatrices have compatible dimensions)*
### Addition:
$$\begin{bmatrix}
A & B \\
C & D
\end{bmatrix} + \begin{bmatrix}
E & F \\
G & H
\end{bmatrix} = \begin{bmatrix}
A + E & B + F \\
C + G & D + H
\end{bmatrix}$$
### Multiplication:
$$\begin{bmatrix}
A & B & C \\
D & E & F
\end{bmatrix} \cdot \begin{bmatrix}
G \\
H \\
I
\end{bmatrix} = \begin{bmatrix}
AG + BH + CI \\
DG + EH + FI
\end{bmatrix}$$
$$\begin{bmatrix}
A & B \\
C & D
\end{bmatrix} \cdot \begin{bmatrix}
E & F \\
G & H
\end{bmatrix} = \begin{bmatrix}
AE + BG & AF + BH \\
CE + DG & CF + DH
\end{bmatrix}$$
- Note that, say, $AE$ is computed by matrix multiplication!
## Block Matrices in [[1.5 - Getting Started with MATLAB|MATLAB]]
You can define a matrix using block matrix notation in MATLAB:
```matlab
>> A = [4];
>> B = [-1 3];
>> C = [2; 4];
>> D = [1 2; 3 4];
>> S = [A B; C D];
```
## Systems of Linear Equations
### The Problem
- Want to find the values of $x_1, x_2, \ldots, x_n$ such that:
$$\begin{align}
a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &= b_1 \\
a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &= b_2 \\
&\vdots \\
a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n &= b_n
\end{align}$$
- All $a_{ij}$ and $b_i$ are given constants
- This is $n$ equations with $n$ unknowns
### Matrix Form
Let:
$$A = \begin{bmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{bmatrix}, \quad x = \begin{bmatrix}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{bmatrix}, \quad b = \begin{bmatrix}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{bmatrix}$$
- Restate the problem as: Find $x \in \mathbb{R}^n$ satisfying
$$\boxed{Ax = b}$$where $A \in \mathbb{R}^{n \times n}$, $b \in \mathbb{R}^n$ are given
- Recall that $\mathbb{R}^n$ is the set of all real vectors of length $n$
- In scientific computing, we always use column vectors
### Solution Existence
- If $A$ is nonsingular, then $A^{-1}$ exists. The system $Ax = b$ **always** has a unique solution $$\boxed{x = A^{-1}b}$$
- If $A$ is singular, the number of solutions to $Ax = b$ depends on $b$ (either no solutions or infinitely many solutions)
## Lower Triangular Matrices
- A matrix $L$ is called a **lower triangular matrix** if all of its entries above the main diagonal are zero, i.e. if $l_{ij} = 0$ for $i < j$
### Examples of Lower Triangular Matrices:
$$\begin{bmatrix}
1 & 0 \\
4 & -2
\end{bmatrix}, \quad \begin{bmatrix}
0 & 0 & 0 \\
4 & 1 & 0 \\
0 & 6 & -2
\end{bmatrix}, \quad \begin{bmatrix}
1 & 0 & 0 & 0 \\
1 & -1 & 0 & 0 \\
0 & 4 & 2 & 0 \\
2 & 2 & -10 & -6
\end{bmatrix}$$
## Exercise 1: Solving Lower-Triangular Linear Systems
Solve:
$$\begin{bmatrix}
-1 & 0 & 0 \\
-1 & 2 & 0 \\
2 & 0 & 3
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3
\end{bmatrix} = \begin{bmatrix}
1 \\
0 \\
-2
\end{bmatrix}$$
### Solution:
- $-x_1 = 1 \Rightarrow x_1 = -1$
- $-x_1 + 2x_2 = 0 \Rightarrow 2x_2 = x_1 \Rightarrow 2x_2 = -1 \Rightarrow x_2 = -1/2$
- $2x_1 + 3x_3 = -2 \Rightarrow 3x_3 = -2 - 2x_1 \Rightarrow 3x_3 = -2 - 2(-1) = 0 \Rightarrow x_3 = 0$
Therefore, $x_1 = -1$, $x_2 = -1/2$, and $x_3 = 0$
>โอเค อันนี้คือง่ายมาก มองออกเลย คือเขาให้มาเป็น Lower-Triangular อยู่แล้ว มันก็สามารถหา $x_1$ ได้เลย แล้วก็เอาไปหา $x_2$ ต่อ จนได้คำตอบทั้งหมด ง่ายใช่มั้ยล่ะ! — เราจะเรียกสิ่งนี้ว่า Forward Substitution
## Solving Lower-Triangular Linear Systems in General
- Consider solving $Lx = b$, where $L$ is lower triangular
- The system is:
$$\begin{bmatrix}
l_{11} & 0 & 0 & \cdots & 0 \\
l_{21} & l_{22} & 0 & \cdots & 0 \\
l_{31} & l_{32} & l_{33} & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn}
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\vdots \\
x_n
\end{bmatrix} = \begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\vdots \\
b_n
\end{bmatrix}$$
### Solution Process:
- Multiply the first row of $L$ with $x$ yields: $l_{11}x_1 = b_1$
- So: $x_1 = \frac{b_1}{l_{11}}$
- Multiplying the second row of $L$ with $x$ yields: $l_{21}x_1 + l_{22}x_2 = b_2$
- Since now $x_1$ is known, we have:
- $l_{22}x_2 = b_2 - l_{21}x_1$
- $x_2 = \frac{b_2 - l_{21}x_1}{l_{22}}$
- Similarly, doing the same thing for the next row and so on yields:
- $x_3 = \frac{b_3 - l_{31}x_1 - l_{32}x_2}{l_{33}}$
$\vdots$
- $x_n = \frac{b_n - l_{n1}x_1 - l_{n2}x_2 - \cdots - l_{n,n-1}x_{n-1}}{l_{nn}}$
### Summary Formula:
$$\boxed{x_i = \left(b_i - \sum_{j=1}^{i-1} l_{ij}x_j\right) / l_{ii}, \quad i = 1, \ldots, n}$$
Note that $\sum_{j=1}^{0} l_{ij}x_j = 0$ by convention.
> [!NOTE] Theorem
> A square upper or lower triangular matrix is nonsingular if and only if all of its diagonal entries are nonzeros.
>
> ง่าย ๆ ถ้ามี diagonal element = 0 แค่ตัวเดียว ผลคูณทั้งหมดจะกลายเป็น 0 → determinant = 0 → matrix singular → ไม่มี inverse → แก้สมการเชิงเส้นบางแบบจะหาคำตอบไม่ได้
## Forward Substitution Algorithm Ver. 1
[[Forward Substitution Algorithms in Swift]]
```
for i = 1 to n
if l_ii = 0 then stop. (The matrix is singular in this case)
for j = 1 to i-1
b_i = b_i - l_ij * x_j
end
x_i = b_i / l_ii
end
```
- What is the running time of this algorithm? $O(n^2)$, where $n$ is the number of rows/columns of $L$
## Forward Substitution Algorithm Ver. 2
```
for j = 1 to n
if l_jj = 0 then stop. (The matrix is singular in this case)
x_j = b_j / l_jj
for i = j+1 to n
b_i = b_i - l_ij * x_j
end
end
```
- Note that the ordering of the loop indices is different from ver. 1 (may not be intuitive)
- Same running time as ver. 1
## How Version 2 works on Exercise 1
$$\begin{bmatrix}
-1 & 0 & 0 \\
-1 & 2 & 0 \\
2 & 0 & 3
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3
\end{bmatrix} = \begin{bmatrix}
1 \\
0 \\
-2
\end{bmatrix}$$
- $x_1 = b_1/l_{11} = 1/(-1) = -1$
- $b_2 := b_2 - l_{21}x_1 = 0 - (-1)(-1) = -1$
- $b_3 := b_3 - l_{31}x_1 = -2 - 2(-1) = 0$
- $x_2 = b_2/l_{22} = -1/2$
- $b_3 := b_3 - l_{32}x_2 = 0 - 0(-1/2) = 0$
- $x_3 = b_3/l_{33} = 0/3 = 0$
## Is Version 2 Better?
### The inner for-loop of Version 1 is:
• $b_i = b_i - l_{i1}x_1$
• $b_i = b_i - l_{i2}x_2$
• $\vdots$
• $b_i = b_i - l_{i,i-1}x_{i-1}$
### The inner for-loop of Version 2 is:
• $b_{j+1} = b_{j+1} - l_{j+1,j}x_j$
• $b_{j+2} = b_{j+2} - l_{j+2,j}x_j$
• $\vdots$
• $b_n = b_n - l_{nj}x_j$
### Key Difference:
- In Version 1, the inner for-loop must be performed **sequentially**
- In Version 2, the inner for-loop can be done **concurrently**
- (Optional) BLAS (Basic Linear Algebra Subprograms) levels
## Upper Triangular Matrices
- A matrix $U$ is called an **upper triangular matrix** if all of its entries below the main diagonal are zero, i.e. if $u_{ij} = 0$ for $i > j$
### Examples of Upper Triangular Matrices:
$$\begin{bmatrix}
2 & 4 \\
0 & -2
\end{bmatrix}, \quad \begin{bmatrix}
0 & 4 & 1 \\
0 & 1 & 5 \\
0 & 0 & 0
\end{bmatrix}, \quad \begin{bmatrix}
-1 & 0 & 0 & 4 \\
0 & -2 & 1 & 2 \\
0 & 0 & -3 & 9 \\
0 & 0 & 0 & -4
\end{bmatrix}$$
## Exercise 2: Solving Upper-Triangular Linear Systems
Solve:
$$\begin{bmatrix}
3 & 4 & 0 \\
0 & 1 & -1 \\
0 & 0 & -2
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3
\end{bmatrix} = \begin{bmatrix}
0 \\
2 \\
2
\end{bmatrix}$$
### Solution:
- $-2x_3 = 2 \Rightarrow x_3 = 2/(-2) = -1$
- $x_2 - x_3 = 2 \Rightarrow x_2 = 2 + x_3 = 2 + (-1) = 1$
- $3x_1 + 4x_2 = 0 \Rightarrow 3x_1 = -4x_2 \Rightarrow 3x_1 = -4(1) = -4 \Rightarrow x_1 = -4/3$
Therefore, $x_1 = -4/3$, $x_2 = 1$, and $x_3 = -1$
>แบบนี้ก็จะคล้ายข้างบนเรยร่ะ แต่ว่ามันเป็น Upper Triangular แทน ก็จะเริ่มหา $x$ ได้จากล่าง ๆ ก่อนนะ ก็จะเรียกว่า Backward Substitution แทน
## Solving Upper-Triangular Linear Systems in General
- Consider solving $Ux = b$, where $U$ is upper triangular
- The system is:
$$\begin{bmatrix}
u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\
0 & u_{22} & u_{23} & \cdots & u_{2n} \\
0 & 0 & u_{33} & \cdots & u_{3n} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & u_{nn}
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\vdots \\
x_n
\end{bmatrix} = \begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\vdots \\
b_n
\end{bmatrix}$$
### Solution Process:
- Multiplying the last row of $U$ with $x$ yields: $u_{nn}x_n = b_n$
- So: $x_n = \frac{b_n}{u_{nn}}$
- Multiplying the second-to-last row of $U$ with $x$ yields: $u_{n-1,n-1}x_{n-1} + u_{n-1,n}x_n = b_{n-1}$
- Since now $x_n$ is known, we have:
- $u_{n-1,n-1}x_{n-1} = b_{n-1} - u_{n-1,n}x_n$
- $x_{n-1} = \frac{b_{n-1} - u_{n-1,n}x_n}{u_{n-1,n-1}}$
- Similarly, doing the same thing for the row above and so on yields:
- $x_{n-2} = \frac{b_{n-2} - u_{n-2,n-1}x_{n-1} - u_{n-2,n}x_n}{u_{n-2,n-2}}$
$\vdots$
- $x_1 = \frac{b_1 - u_{12}x_2 - u_{13}x_3 - \cdots - u_{1n}x_n}{u_{11}}$
### Summary Formula:
$$\boxed{x_i = \left(b_i - \sum_{j=i+1}^{n} u_{ij}x_j\right) / u_{ii}, \quad i = n, \ldots, 1}$$
Note that $\sum_{j=n+1}^{n} u_{ij}x_j = 0$ by convention.
> [!quote] END OF WEEK 1
### Back Substitution Algorithm Version 1
```
for i = n to 1
if u_ii = 0 then stop. (The matrix is singular)
for j = i+1 to n
b_i = b_i - u_ij * x_j
end
x_i = b_i / u_ii
end
```
- **Running time**: $O(n^2)$
### Back Substitution Algorithm Version 2
```
for j = n to 1
if u_jj = 0 then stop. (The matrix is singular)
x_j = b_j / u_jj
for i = 1 to j-1
b_i = b_i - u_ij * x_j
end
end
```
- **Running time**: $O(n^2)$
- **Note**: Version 2 has better concurrency than Version 1
- **Important**: Pay attention to the ordering of loop indices
### Example: Version 2 on Exercise 2
Given system:
$$\begin{bmatrix} 3 & 4 & 0 \\ 0 & 1 & -1 \\ 0 & 0 & -2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 2 \\ 2 \end{bmatrix}$$
**Steps:**
- $x_3 = b_3/u_{33} = 2/(-2) = -1$
- $b_1 := b_1 - u_{13}x_3 = 0 - 0(-1) = 0$
- $b_2 := b_2 - u_{23}x_3 = 2 - (-1)(-1) = 1$
- $x_2 = b_2/u_{22} = 1/1 = 1$
- $b_1 := b_1 - u_{12}x_2 = 0 - 4(1) = -4$
- $x_1 = b_1/u_{11} = -4/3$
## Gaussian Elimination (GE)
### Main Idea
- Transform matrix $A$ into upper triangular form, which we know how to solve
- The transformation is known as ==**Gaussian Elimination (GE)**==
### Gaussian Elimination Concepts
**Basic Process:**
- Eliminate (zero out) entries below the main diagonal one column at a time from left to right
**Algorithm Structure:**
```
For k = 1 to n (column to eliminate)
For i = k+1 to n (row to eliminate)
Eliminate a_ik
End
End
```
**Elimination Formula:**
$$\boxed{\text{(Row } i\text{)} := \text{(Row } i\text{)} - \left(\frac{a_{ik}}{a_{kk}}\right) \cdot \text{(Row } k\text{)}}$$
- $a_{kk}$ is called a **pivot**
- $\frac{a_{ik}}{a_{kk}}$ is called a **multiplier**
### Example 3: Solving a System with GE
**System:**
$$\begin{align}
x_1 + 2x_2 + x_3 &= -1 \\
-3x_1 + x_2 + x_3 &= 0 \\
x_1 + 3x_3 &= 1
\end{align}$$
**Matrix form:**
$$\begin{bmatrix} 1 & 2 & 1 \\ -3 & 1 & 1 \\ 1 & 0 & 3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix}$$
**Step 1 - Eliminate first column:**
- $\text{(Row 2)} := \text{(Row 2)} - \frac{-3}{1} \cdot \text{(Row 1)}$
- $\text{(Row 3)} := \text{(Row 3)} - \frac{1}{1} \cdot \text{(Row 1)}$
Result:
$$\begin{bmatrix} 1 & 2 & 1 \\ 0 & 7 & 4 \\ 0 & -2 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} -1 \\ -3 \\ 2 \end{bmatrix}$$
**Step 2 - Eliminate second column:**
- $\text{(Row 3)} := \text{(Row 3)} - \frac{-2}{7} \cdot \text{(Row 2)}$
Final upper triangular form:
$$\begin{bmatrix} 1 & 2 & 1 \\ 0 & 7 & 4 \\ 0 & 0 & \frac{22}{7} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} -1 \\ -3 \\ \frac{8}{7} \end{bmatrix}$$
**Solution by back substitution:**
- $x_1 = -\frac{1}{11}$
- $x_2 = -\frac{7}{11}$
- $x_3 = \frac{4}{11}$
___
## LU Factorization
### The L Matrix Construction
Create lower-triangular matrix $L$ with same dimension as $A$:
- **0** for entries above diagonal
- **1** for entries on diagonal
- **Multipliers** for entries below diagonal (put multiplier used to zero out $a_{ij}$ in position $l_{ij}$)
>ก็คือการจะ Construct L matrix ได้ก็ต้องทำ GE มาก่อน แล้วเอา Multiplier จากด้านบนมาใส่เลยนะ!
For Example 3:
$$L = \begin{bmatrix} 1 & 0 & 0 \\ -3 & 1 & 0 \\ 1 & -\frac{2}{7} & 1 \end{bmatrix}$$
### LU Factorization Theorem
$$\boxed{A = LU}$$
Where:
- $L$ is **unit lower triangular** (lower triangular with ones on main diagonal)
- $U$ is the upper triangular matrix from GE
>บางคนอาจจะถามแล้ว $U$ ไปเอามาจากไหนล่ะ?? ก็ที่ได้หลังจาก GE เสร็จแล้วนั่นเอง เริ่ดเรยหร่ะ
**Why this works:** (ไม่ต้องรู้ก็ได้น้า ทำเป็นพอ)
The elimination steps can be written as:
- $u_{1\bullet} = a_{1\bullet}$
- $u_{2\bullet} = a_{2\bullet} - (-3)u_{1\bullet}$
- $u_{3\bullet} = a_{3\bullet} - u_{1\bullet} - (-\frac{2}{7})u_{2\bullet}$
![[Slides_CSS322_2_Linear_Systems_handout-41-45.pdf.pdf]]
### Procedure for Solving $Ax = b$ using LU Factorization
$$\boxed{
\begin{align}
&\text{1. Factor } A = LU \\
&\text{2. Solve } Lw = b \text{ for } w \text{ (forward substitution)} \\
&\text{3. Solve } Ux = w \text{ for } x \text{ (back substitution)}
\end{align}
}$$
**Verification:**
$$Ax = (LU)x = L(Ux) = Lw = b$$
![[Slides_CSS322_2_Linear_Systems_handout-Example3Steps.pdf.pdf]]
### LU-Factorization Algorithm
```
For k = 1 to n (column to eliminate)
piv = a_kk
For i = k+1 to n (row to eliminate)
mult = a_ik / piv
// (Row i) := (Row i) - mult * (Row k)
For j = k to n
a_ij = a_ij - mult * a_kj
End
l_ik = mult
End
l_kk = 1
End
```
- **Running time**: $O(n^3)$
- Matrix $A$ is replaced by $U$ on exit
### Advantages of LU Factorization
- **Multiple right-hand sides**: If solving $Ax_1 = b_1$, $Ax_2 = b_2$, $Ax_3 = b_3$
- Factor $A = LU$ once
- Repeat steps 2 and 3 for each right-hand side vector
- More efficient than performing GE three separate times
### Time Complexity Analysis
**Total time for solving $Ax = b$:**
- $O(n^3)$ to factor $A = LU$
- $O(n^2)$ for forward substitution
- $O(n^2)$ for back substitution
- **Total**: $O(n^3)$ operations
**Note**: Same operation count as plain GE on $[A \, b]$ followed by back substitution
## Gaussian Elimination with Partial Pivoting (GEPP)
### Problem with Plain GE
- A pivot may be 0 even if the matrix is nonsingular
- **Example**: $\begin{bmatrix} 0 & 2 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 2 \\ 4 \end{bmatrix}$
### Pivoting Strategy
Before each elimination step:
- Exchange rows so that element with ==**maximum absolute value**== in pivot column (among uneliminated entries) is brought to pivot position
### Example 5: GEPP in Action
**System:**
$$\begin{bmatrix} 2 & 4 & 5 \\ 1 & -1 & -1 \\ 3 & -2 & 1 \end{bmatrix} x = \begin{bmatrix} 3 \\ -1 \\ -5 \end{bmatrix}$$
**First column**: $|3|$ is maximum among $\{|2|, |1|, |3|\}$, so swap Rows 1 and 3:
$$\begin{bmatrix} 3 & -2 & 1 \\ 1 & -1 & -1 \\ 2 & 4 & 5 \end{bmatrix} x = \begin{bmatrix} -5 \\ -1 \\ 3 \end{bmatrix}$$
**Eliminate first column:**
$$\begin{bmatrix} 3 & -2 & 1 \\ 0 & -\frac{1}{3} & -\frac{4}{3} \\ 0 & \frac{16}{3} & \frac{13}{3} \end{bmatrix} x = \begin{bmatrix} -5 \\ \frac{2}{3} \\ \frac{19}{3} \end{bmatrix}$$
**Second column**: $|\frac{16}{3}| > |-\frac{1}{3}|$, so swap Rows 2 and 3:
$$\begin{bmatrix} 3 & -2 & 1 \\ 0 & \frac{16}{3} & \frac{13}{3} \\ 0 & -\frac{1}{3} & -\frac{4}{3} \end{bmatrix} x = \begin{bmatrix} -5 \\ \frac{19}{3} \\ \frac{2}{3} \end{bmatrix}$$
**Eliminate second column:**
$$\begin{bmatrix} 3 & -2 & 1 \\ 0 & \frac{16}{3} & \frac{13}{3} \\ 0 & 0 & -\frac{17}{16} \end{bmatrix} x = \begin{bmatrix} -5 \\ \frac{19}{3} \\ \frac{17}{16} \end{bmatrix}$$
- The we finish finding $x$ by backward substitution.
### Permutation Matrices
- **Definition**: Square matrix with entries 0 or 1, exactly one 1 per row and column
- **Example**:
$$\begin{bmatrix} 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$
- **Properties:**
1. **Left multiplication** by permutation matrix permutes rows:
$$\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix} = \begin{bmatrix} 3 & 4 \\ 5 & 6 \\ 1 & 2 \end{bmatrix}$$
2. **Right multiplication** by permutation matrix permutes columns:
$$\begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \end{bmatrix} \cdot \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix} = \begin{bmatrix} 2 & 1 & 4 & 3 \\ 6 & 5 & 8 & 7 \end{bmatrix}$$
3. $PP^T = I$, i.e., $P^{-1} = P^T$
### GEPP and Matrix Factorization
- **Key Theorem**: GEPP computes factorization: $\boxed{PA = LU}$
- Or equivalently: $\boxed{A = P^T LU}$
Where $P$ is the permutation matrix representing all row swaps performed.
![[Slides_CSS322_2_Linear_Systems_Example6.pdf.pdf]]
### GEPP Algorithm
```
P = I
For k = 1 to n (column to eliminate)
Find p = argmax{|a_pk| : p = k, k+1, ..., n}
Swap rows k and p of A, L, and P
piv = a_kk
If piv = 0 then A is singular. Stop.
For i = k+1 to n (row to eliminate)
mult = a_ik / piv
For j = k to n (update the whole row)
a_ij = a_ij - mult * a_kj
End
l_ik = mult
End
l_kk = 1
End
```
### Finding Permutation Matrix P
**Easiest method:**
- Initially set $P = I$
- As you swap rows of $A$ and $L$ during pivoting, swap rows of $P$ accordingly
### Procedure for Solving $Ax = b$ with GEPP
$$\boxed{
\begin{align}
&\text{1. Factor } A = P^T LU \text{ by GEPP} \\
&\text{2. Set } \hat{b} = Pb \\
&\text{3. Solve } Lw = \hat{b} \text{ for } w \text{ (forward substitution)} \\
&\text{4. Solve } Ux = w \text{ for } x \text{ (back substitution)}
\end{align}
}$$
**Verification:**
$Ax = (P^T LU)x = P^T L(Ux) = P^T(Lw) = P^T\hat{b} = b$
>ถ้าถามว่า เอ๊ะ แล้วอยู่ ๆ $P$ หรือ $P^T$ โผล่มาได้ยังไงเนี่ย ก็ต้องตอบว่า จริง ๆ แล้วมันไม่ได้ GEPP ไม่ได้ compute $A$ ตรง ๆ เนอะ มันคือตรงสลับ Row ก่อน จึงเป็น $PA=LU$
> [!NOTE] Theorem
> GEPP encounters a zero pivot if and only if $A$ is singular
### Running Time and Comments
- **GEPP running time**: Still $O(n^3)$ (pivoting takes $O(n^2)$ total)
- **Why factor $A = P^T LU$ rather than GEPP on $[A \, b]$?**
- Same operation count for single system
- More efficient for multiple right-hand sides
![[Slides_CSS322_2_Linear_Systems_handout-ex7.pdf.pdf]]
![[Slides_CSS322_2_Linear_Systems_handout-ex8-full.pdf.pdf]]
### Partial vs Complete Pivoting
- **Partial pivoting**: Only swap rows (what we've discussed)
- **Complete pivoting (GECP)**: Exchange both rows and columns to bring largest uneliminated entry to pivot position
- Needs $O(n^3)$ operatio ns for pivoting
- Tiny improvement in accuracy
- Generally not needed
### Notes on Implementation Differences
**Mathematical approach** (often taught in math courses):
- "Multiply rows by constants and add/subtract rows to eliminate entries"
- Useful for small systems by hand and proving theorems
- **Problem**: Won't give correct matrix $L$ for factorization
- No advantages for computer implementation
**Computational approach** (what we use):
- Specific elimination formula with multipliers
- Gives correct $L$ matrix
- Optimal for computer implementation
### Comments on Pivoting Strategy
- **Can we swap rows arbitrarily** (as long as we avoid zero pivots)?
- Yes, in exact arithmetic
- May be helpful for proving properties
- **In floating-point arithmetic**: Using maximum absolute value entry as pivot yields lowest error
- **If goal is only to solve**: Additional row exchanges are unnecessary and increase running time
## Computing Matrix Inverses
### Why Not Gauss-Jordan Elimination?
**Problems with GJE:**
- Less accurate than GEPP-based method in floating-point arithmetic
- About 50% more operations than GEPP method
- **However**: GJE is mathematically elegant and useful for:
- Proving theorems
- Parallel computing
### Method Using GEPP
To compute $A^{-1}$ where $A \in \mathbb{R}^{n \times n}$:
**Setup:**
- Let $x_i$ = $i$-th column of $A^{-1}$
- Let $e_i$ = $i$-th column of $I_n$
**Examples of standard basis vectors:**
- For $A \in \mathbb{R}^{3 \times 3}$: $e_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$, $e_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}$, $e_3 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$
**Key insight**: Since $AA^{-1} = I$:
$A \cdot [x_1 \, x_2 \, \cdots \, x_n] = [e_1 \, e_2 \, \cdots \, e_n]$
This means: $Ax_i = e_i$ for all $i$
**Algorithm:**
1. Factor $A = P^T LU$ once using GEPP
2. Solve $Ax_1 = e_1, Ax_2 = e_2, \ldots, Ax_n = e_n$ using the factorization
3. Combine columns $[x_1 \, x_2 \, \cdots \, x_n]$ to get $A^{-1}$
### Example 9: Computing Inverse
Find inverse of:
$$A = \begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix}$$
**Step 1**: GEPP gives:
$L = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ \frac{1}{4} & \frac{1}{2} & 1 \end{bmatrix}, \quad U = \begin{bmatrix} 4 & 4 & 2 \\ 0 & 2 & 2 \\ 0 & 0 & \frac{1}{2} \end{bmatrix}, \quad P = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}$
**For finding $x_1$ (first column of $A^{-1}$):**
- Set $\hat{e_1} = Pe_1 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$
- Solve $Lw_1 = \hat{e_1}\xrightarrow{}w_1 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$
- Solve $Ux_1 = w_1\xrightarrow{}x_1 = \begin{bmatrix} 1 \\ -2 \\ 2 \end{bmatrix}$
**For finding $x_2$ (second column of $A^{-1}$):**
- Set $\hat{e_2} = Pe_2 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$
- Solve $Lw_2 = \hat{e_2} \xrightarrow{} w_2 = \begin{bmatrix} 1 \\ -1 \\ \frac{1}{4} \end{bmatrix}$
- Solve $Ux_2 = w_2 \xrightarrow{} x_2 = \begin{bmatrix} 1 \\ -1 \\ \frac{1}{2} \end{bmatrix}$
**For finding $x_3$ (third column of $A^{-1}$):**
- Following similar steps → $x_3 = \begin{bmatrix} -1 \\ 1.5 \\ -1 \end{bmatrix}$
**Result:**
$$A^{-1} = \begin{bmatrix} 1 & 1 & -1 \\ -2 & -1 & 1.5 \\ 2 & 0.5 & -1 \end{bmatrix}$$
## Important Computational Advice
### Computing $A^{-1}v$ Efficiently
**Question**: Given known vector $v$ and matrix $A$, what's the best way to compute $A^{-1}v$?
$$\color{}\large\boxed{\text{Best approach: Solve } Ax = v \text{ for } x \text{ using GEPP}}$$
**Why?**
- If $Ax = v$, then $x = A^{-1}v$
- **DO NOT** explicitly compute $A^{-1}$ and multiply by $v$
- Solving $Ax = v$ is faster and more accurate than computing $A^{-1}$ first
**Reason**: GEPP for solving $Ax = v$ is:
- Faster than computing $A^{-1}$ then multiplying
- More numerically accurate in floating-point arithmetic
- Uses fewer operations overall
### Example 10: Computing $A^{-1}v$ Correctly
Find:
$$\begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix}^{-1} \cdot \begin{bmatrix} -2 \\ -1 \\ -2 \end{bmatrix}$$
**Correct approach**: Solve:
$\begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix} x = \begin{bmatrix} -2 \\ -1 \\ -2 \end{bmatrix}$
Using the 4-step GEPP procedure:
1. **Factor** $A = P^T LU$ (same as Example 9)
2. **Set** $\hat{b} = Pb = \begin{bmatrix} -1 \\ -2 \\ -2 \end{bmatrix}$
3. **Solve** $Lw = \hat{b}$ → $w = \begin{bmatrix} -1 \\ -1 \\ -1.25 \end{bmatrix}$
4. **Solve** $Ux = w$ → $x = \begin{bmatrix} -1 \\ 2 \\ -2.5 \end{bmatrix}$
___
## MATLAB Implementation
### Solving Linear Systems
```matlab
>> x = A \ b; % MATLAB chooses appropriate algorithm automatically
```
MATLAB automatically:
- Detects matrix structure (triangular, symmetric, etc.)
- Chooses optimal algorithm
- For general matrices, uses GEPP
### LU Factorization
```matlab
>> [L,U,P] = lu(A); % Returns P^T LU factorization
>> [PtL,U] = lu(A); % Returns P^T L as single matrix
>> [L,U,pvec] = lu(A, 'vector'); % P in vector form where A(pvec,:) = L*U
```
### Multiple Systems with Same Matrix
```matlab
>> [L,U,P] = lu(A);
>> b1hat = P*b1; w1 = L \ b1hat; x1 = U \ w1;
>> b2hat = P*b2; w2 = L \ b2hat; x2 = U \ w2;
>> b3hat = P*b3; w3 = L \ b3hat; x3 = U \ w3;
```
### Matrix Inverse
```matlab
>> Ai = inv(A); % Uses GEPP-based algorithm we discussed
```
### Computing $A^{-1}v$ in MATLAB
```matlab
>> x = A \ v; % CORRECT - solve directly
>> x = inv(A) * v; % AVOID - less efficient and accurate
```
**For documentation:**
```matlab
>> doc mldivide; % or help mldivide
>> doc lu; % or help lu
```
## Key Takeaways
### Algorithm Complexity Summary
- **Back substitution**: $O(n^2)$
- **Forward substitution**: $O(n^2)$
- **LU factorization**: $O(n^3)$
- **GEPP**: $O(n^3)$
- **Solving $Ax = b$**: $O(n^3)$ total
### Best Practices
1. **For single system**: Use GEPP (4-step procedure) or MATLAB's `A \ b`
2. **For multiple systems with same $A$**: Factor once, reuse for each $b$
3. **For computing $A^{-1}v$**: Solve $Ax = v$ directly, don't compute $A^{-1}$
4. **For numerical stability**: Always use partial pivoting
5. **For implementation**: Use computational approach with multipliers, not mathematical row operations
### Important Distinctions
- **Plain GE** vs **GEPP**: GEPP handles singular matrices and improves numerical stability
- **Partial** vs **Complete pivoting**: Partial pivoting (row swaps only) is sufficient for most applications
- **Mathematical** vs **Computational approach**: Computational approach gives correct $L$ matrix and better implementation
> [!quote] END OF WEEK 2
# Solving Modified Linear Systems
>อันนี้คือแล้วถ้า Matrix มันไม่ใช่ $A$ ธรรมดาแต่เป้น $A-uv^T$ ล่ะ จะแก้ยังไงดีเอ่ย
>จริง ๆ แล้วจะ compute ให้กลายเป็น Matrix ก้อนเดียวก่อนก็ได้ แล้วหาคำตอบ
>แต่มันสามารถใช้ $P^TLU$ factor to use to find $x$ in order of $n^2$ population
>- **You already have** the PTLUPTLU factorization of AA
>- **Don't throw it away** by computing (A−uvT)(A−uvT) explicitly
>- **Reuse the factorization** through Sherman-Morrison
>- **Solve in O(n2)O(n2)** operations instead of O(n3)O(n3)
## Rank-One Modification Problem
- **Problem Statement**: Given the $P^T LU$ factorization of matrix $A$, how to solve $(A + uv^T)x = b$ efficiently?
- **Key Concepts**:
- This is known as a **rank-one modification**
- The matrix $uv^T$ has rank one
- Any rank-one matrix can be expressed as $uv^T$ for some vectors $u$ and $v$
- **Example**: If a single entry of $A$ changes from $a_{jk}$ to $\tilde{a}_{jk}$, then:
- New matrix: $A + \alpha e_j e_k^T$
- Where $\alpha = a_{jk} - \tilde{a}_{jk}$
## Sherman-Morrison Formula
- **Formula**:
$$\boxed{(A + uv^T)^{-1} = A^{-1} + A^{-1}u(1 - v^T A^{-1}u)^{-1}v^T A^{-1}}$$
• **Solution to $(A + uv^T)x = b$**:
$$x = (A + uv^T)^{-1}b = A^{-1}b + A^{-1}u(1 - v^T A^{-1}u)^{-1}v^T A^{-1}b$$
- **Implementation Notes**:
- Involves computing $A^{-1}b$ and $A^{-1}u$
- Should solve linear systems $Ay = b$ and $Az = u$ instead of computing inverses
- $y$ may already be known from previous computations
## Rank-One Updating Algorithm
- **Step-by-step process**:
1. Solve $Az = u$ for $z$
2. Solve $Ay = b$ for $y$
3. Compute: $\boxed{x = y + \frac{v^T y}{1 - v^T z}z}$
- **Efficiency**: If factors of $A$ are already known, this algorithm requires $O(n^2)$ operations
## Woodbury Formula (Generalization)
- **Formula**:
$$\boxed{(A + UV^T)^{-1} = A^{-1} + A^{-1}U(I - V^T A^{-1}U)^{-1}V^T A^{-1}}$$
- Where $U$ and $V$ are $n \times k$ matrices
- Generalizes Sherman-Morrison formula to rank-$k$ modifications
## Example 11: Numerical Application
>ถ้าเห็นจากตัวอย่างอะนะ คือมันไม่ใช่ Matrix $A$ นะที่เอาไป Compute อะ แต่มันต่างกันแค่นิดเดียว จาก 6 เป็น 4 แบบนี้จะทำใหม่ทั้งหมดเลยรึ? No!!
- **Given**: $P^T LU$ factors of matrix
$$A = \begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix}$$
- **Solve**:
$$\begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 4 & 4 \end{bmatrix} x = \begin{bmatrix} 3 \\ 6 \\ 10 \end{bmatrix}$$
- **Solution Process**:
- New matrix = $A - \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \times [0, 2, 0]$
- So $u = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$ and $v = \begin{bmatrix} 0 \\ 2 \\ 0 \end{bmatrix}$
- Solve $Az = u$: $z = \begin{bmatrix} -1 \\ 1.5 \\ -1 \end{bmatrix}$
- Solve $Ay = b$: $y = \begin{bmatrix} -1 \\ 3 \\ -1 \end{bmatrix}$
- Final solution: $x = y + \frac{v^T y}{1 - v^T z}z = \begin{bmatrix} -1 \\ 3 \\ -1 \end{bmatrix} + \frac{6}{1-3)}\begin{bmatrix} -1 \\ 1.5 \\ -1 \end{bmatrix} = \begin{bmatrix} 2 \\ -1.5 \\ 2 \end{bmatrix}$
## Example 12: Non-existence of LU Factorization
- **Problem**: Prove that $A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$ has no LU factorization
- **Proof by Contradiction**:
- Assume LU factorization exists: $\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} l_{11} & 0 \\ l_{21} & l_{22} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} \\ 0 & u_{22} \end{bmatrix}$
>อันนี้มัน Assume ง่าย ๆ ได้เลยล่ะ $A=LU$
- From first row × first column: $0 = l_{11}u_{11}$ → either $l_{11} = 0$ or $u_{11} = 0$
- From first row × second column: $1 = l_{11}u_{12}$ → $l_{11} \neq 0$
- Therefore $u_{11} = 0$
- From second row × first column: $1 = l_{21}u_{11}$ → contradiction since $u_{11} = 0$
- **Conclusion**: No LU factorization exists
## Example 13: General 2×2 LU Factorization
- **Matrix**: $A = \begin{bmatrix} 1 & a \\ c & b \end{bmatrix}$
- **LU Factorization**:
- $L = \begin{bmatrix} 1 & 0 \\ c & 1 \end{bmatrix}$
- $U = \begin{bmatrix} 1 & a \\ 0 & b-ac \end{bmatrix}$
- **Singularity Condition**: Matrix $A$ is singular when $b - ac = 0$ (จริง ๆ แค่นี้พอแล้ว แต่ใช้ Theorem อื่นได้ ข้างล่าง)
- Since $L$ is nonsingular (diagonal entries are nonzero)
- $A$ is singular ⟺ $U$ is singular ⟺ $b - ac = 0$
## Example 14: Efficient Formula Implementation
- **Problem**: Best way to implement $x = A^{-1}(L^{-1} + B)v$ for efficiency and accuracy
- **Solution Steps**:
1. Use forward substitution to solve $Lu = v$ for $u$
2. Set $w := u + Bv$
3. Use GEPP to solve $Ax = w$ for $x$
- **Reasoning**:
- $x = A^{-1}(L^{-1} + B)v$
- $Ax = (L^{-1} + B)v = L^{-1}v + Bv$
- Let $u = L^{-1}v$ (solve $Lu = v$)
- Let $w = u + Bv$
- Solve $Ax = w$ → Solve by GEPP
- **MATLAB Implementation**:
```matlab
function x = solve_ex14(A, B, L, v)
u = L \ v; % MATLAB will use forward substitution
w = u + B*v;
x = A \ w; % MATLAB will use GEPP
end
```
## Example 15: Unconventional Factorization
- **Given**: $A = UPL$ (where $U$ is upper triangular, $P$ is permutation, $L$ is lower triangular)
- **Solve**: $Ax = b$ in $O(n^2)$ operations
- **Solution Steps**:
1. Solve $Uv = b$ for $v$ by back substitution
2. Set $w = P^T v$ (permute entries of $v$)
3. Solve $Lx = w$ for $x$ by forward substitution
- **Reasoning**:
- $Ax = UPLx = b$
- Let $v = PLx$, then $Uv = b$
- From $PLx = v$, we get $Lx = P^T v$
- Set $w = P^T v$ and solve $Lx = w$
- **Algorithm Summary**:
- Back substitution: $Uv = b$
- Permutation: $w = P^T v$
- Forward substitution: $Lx = w$
> [!quote] END OF WEEK 4