# Nakafa Learning Content
> For AI agents: use [llms.txt](https://nakafa.com/llms.txt) for the site index. Markdown versions are available by appending `.md` to content URLs or sending `Accept: text/markdown`.
URL: https://nakafa.com/en/subjects/ai-ds/linear-methods/lu-decomposition
Source: https://raw.githubusercontent.com/nakafaai/nakafa.com/refs/heads/main/packages/contents/material/lesson/ai-ds/linear-methods/lu-decomposition/en.mdx
Learn LU decomposition: factorize matrices into lower/upper triangular forms using Gaussian elimination, with Python implementation and system solving.
---
## What is LU Decomposition?
LU decomposition is an abbreviation for **Lower-Upper decomposition**. This name comes from the two matrices resulting from decomposition:
- **L (Lower)** = Lower triangular matrix that has zero elements above the main diagonal
- **U (Upper)** = Upper triangular matrix that has zero elements below the main diagonal
Why do we need this decomposition? Because solving linear equation systems with triangular matrices is much easier than with regular matrices.
## Basic Concepts of LU Decomposition
LU decomposition is a way to break down a matrix into the multiplication of two simpler matrices. Imagine we have a complex puzzle, then we disassemble it into pieces that are easier to reassemble.
For a given matrix $$A$$, we get:
Visible text: For a given matrix , we get:
```math
P \cdot A = L \cdot U
```
Where $$P$$ is a permutation matrix (row order arranger), $$L$$ is a lower triangular matrix, and $$U$$ is an upper triangular matrix.
Visible text: Where is a permutation matrix (row order arranger), is a lower triangular matrix, and is an upper triangular matrix.
Gaussian elimination computes this LU decomposition through a series of systematic row operations. This process can be written as:
```math
U = L_r \cdot P_r \cdot \ldots \cdot L_3 \cdot P_3 \cdot L_2 \cdot P_2 \cdot L_1 \cdot P_1 \cdot A
```
with $$i = 1, \ldots, r$$, row exchange matrix $$P_i$$ or $$P_i = I$$ and elimination matrix:
Visible text: with , row exchange matrix or and elimination matrix:
```math
L_i = \begin{pmatrix} 1 & & & \\ & \ddots & & \\ & & \lambda_{i+1i} & \\ & & \vdots & \ddots \\ & & \lambda_{mi} & & 1 \end{pmatrix}
```
The inverse of $$L_i$$ is:
Visible text: The inverse of is:
```math
L_i^{-1} = \begin{pmatrix} 1 & & & \\ & \ddots & & \\ & & -\lambda_{i+1i} & \\ & & \vdots & \ddots \\ & & -\lambda_{mi} & & 1 \end{pmatrix}
```
Since $$P_i^{-1} = P_i$$, we can rewrite this sequence of operations into a simpler form until we finally get $$L^{-1} \cdot P \cdot A = U$$.
Visible text: Since , we can rewrite this sequence of operations into a simpler form until we finally get .
## Structure of Decomposition Results
From Gaussian elimination with steps $$j_1, \ldots, j_r$$, we get three matrices with special structure:
Visible text: From Gaussian elimination with steps , we get three matrices with special structure:
```math
\begin{pmatrix} 0 & u_{1j_1} & \cdots & & \cdots & u_{1n} \\ l_{21} & 0 & 0 & u_{2j_2} & \cdots & u_{2n} \\ l_{31} & l_{32} & 0 & 0 & \ddots & \\ \vdots & \vdots & \ddots & 0 & 0 & u_{rj_r} \cdots u_{rm} \\ & & & l_{(r+1)r} & 0 & 0 & 0 & 0 \\ \vdots & \vdots & & \vdots & 0 & 0 & 0 & 0 \\ l_{m1} & l_{m2} & \cdots & l_{mr} & 0 & 0 & 0 & 0 \end{pmatrix}
```
Matrix $$P$$ is obtained from applying permutation vector $$p$$ to the rows of the identity matrix. The three matrices have properties:
Visible text: Matrix is obtained from applying permutation vector to the rows of the identity matrix. The three matrices have properties:
1. $$U$$ is an upper triangular matrix in row echelon form with rank $$r$$
2. $$L$$ is a regular lower triangular matrix with $$1$$ on the diagonal
3. $$P$$ is a permutation matrix with $$P^{-1} = P^T$$
Visible text: 1. is an upper triangular matrix in row echelon form with rank
2. is a regular lower triangular matrix with on the diagonal
3. is a permutation matrix with
## Gaussian Elimination Steps
Let's look at a concrete example. Imagine we are cleaning house floor by floor, starting from the top.
1. **First column** with pivot element $$a_{21} = 3$$. We swap the first and second rows like rearranging furniture positions.
```math
A = \begin{pmatrix} 3 & 3 & 3 & 3 \\ 0 & 2 & 4 & 2 \\ 1 & 2 & 3 & 1 \end{pmatrix}
```
```math
r = 1, \text{ Steps } = \{1\}, p = (2,1,3), v = 1
```
Eliminate other entries in that column: $$\lambda_{21} = 0$$, $$\lambda_{31} = -\frac{1}{3}$$
2. **Second column** with pivot element $$a_{22} = 2$$. No exchange needed.
Elimination: $$\lambda_{32} = -\frac{1}{2}$$
3. **Third column** has no pivot element that can be used.
4. **Fourth column** with pivot element $$a_{34} = -1$$. No more operations needed.
Visible text: 1. **First column** with pivot element . We swap the first and second rows like rearranging furniture positions.
Eliminate other entries in that column: ,
2. **Second column** with pivot element . No exchange needed.
Elimination:
3. **Third column** has no pivot element that can be used.
4. **Fourth column** with pivot element . No more operations needed.
Final result:
```math
r = 3, \text{Steps} = \{1, 2, 4\}, p = (2,1,3), v = 1
```
From this elimination result, we obtain:
Component: MathContainer
Children:
```math
U = \begin{pmatrix} 3 & 3 & 3 & 3 \\ 0 & 2 & 4 & 2 \\ 0 & 0 & 0 & -1 \end{pmatrix}
```
```math
L = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1/3 & 1/2 & 1 \end{pmatrix}
```
```math
P = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}
```
The proof shows that $$L \cdot U = P \cdot A$$ for the initially given matrix $$A$$.
Visible text: The proof shows that for the initially given matrix .
## Algorithm and Complexity
The Gaussian elimination algorithm works like systematic cleaning. For each column, we perform:
1. Find the best pivot element
2. Swap rows if necessary
3. Eliminate elements below the pivot
4. Store multiplier information
The time complexity for computing LU decomposition is $$\frac{1}{3}n^3 + O(n^2)$$ arithmetic operations. The proof shows that multiplication and addition occur simultaneously in the elimination row transformation $$a_{ij} := a_{ij} + \lambda_{ir} \cdot a_{rl}$$.
Visible text: The time complexity for computing LU decomposition is arithmetic operations. The proof shows that multiplication and addition occur simultaneously in the elimination row transformation .
Component: MathContainer
Children:
```math
\sum_{j=1}^n \sum_{i=j+1}^n \sum_{l=j+1}^n 1 = \sum_{j=1}^n (n-j)^2
```
```math
= \sum_{j=1}^n (n^2 - 2nj + j^2)
```
```math
= n^3 - 2n \sum_{j=1}^n j + \sum_{j=1}^n j^2
```
```math
= n^3 - 2n \frac{n(n+1)}{2} + \frac{n(n+1)(2n+1)}{6}
```
```math
= n^3 - n^3 + \frac{2}{6}n^3 + O(n^2)
```
```math
= \frac{1}{3}n^3 + O(n^2)
```
## Algorithm Implementation in Python
Here is the implementation of the LU decomposition algorithm in Python:
File: lu_decomposition.py
```python
import numpy as np
def lu_decomposition(A):
"""
Perform LU decomposition with partial pivoting
Input: Matrix A (n x n)
Output: L, U, P (Lower, Upper, and Permutation matrices)
"""
n = A.shape[0]
U = A.copy().astype(float)
L = np.zeros((n, n))
P = np.eye(n)
for i in range(n):
max_row = i
for k in range(i + 1, n):
if abs(U[k, i]) > abs(U[max_row, i]):
max_row = k
if max_row != i:
U[[i, max_row]] = U[[max_row, i]]
P[[i, max_row]] = P[[max_row, i]]
if i > 0:
L[[i, max_row], :i] = L[[max_row, i], :i]
for k in range(i + 1, n):
if U[i, i] != 0:
factor = U[k, i] / U[i, i]
L[k, i] = factor
U[k, i:] = U[k, i:] - factor * U[i, i:]
np.fill_diagonal(L, 1)
return L, U, P
def forward_substitution(L, b):
"""
Forward substitution to solve L * y = b
"""
n = len(b)
y = np.zeros(n)
for i in range(n):
y[i] = b[i] - np.dot(L[i, :i], y[:i])
return y
def backward_substitution(U, y):
"""
Backward substitution to solve U * x = y
"""
n = len(y)
x = np.zeros(n)
for i in range(n - 1, -1, -1):
if U[i, i] != 0:
x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
return x
def solve_linear_system(A, b):
"""
Solve system Ax = b using LU decomposition
"""
L, U, P = lu_decomposition(A)
pb = np.dot(P, b)
y = forward_substitution(L, pb)
x = backward_substitution(U, y)
return x, L, U, P
```
## Solving Linear Systems
After obtaining LU decomposition, solving the system $$A \cdot x = b$$ becomes like solving a puzzle that has been sorted.
Visible text: After obtaining LU decomposition, solving the system becomes like solving a puzzle that has been sorted.
The system $$A \cdot x = b$$ transforms into $$U \cdot x = L^{-1} \cdot P \cdot b$$ through two stages:
Visible text: The system transforms into through two stages:
### Forward Substitution
The forward substitution algorithm solves $$L \cdot y = P \cdot b$$ with regular lower triangular matrix $$L$$ and vector $$c = P \cdot b$$. This process is like filling stairs from bottom to top, where each step depends on the previous step.
Visible text: The forward substitution algorithm solves with regular lower triangular matrix and vector . This process is like filling stairs from bottom to top, where each step depends on the previous step.
For each row $$i = 1, \ldots, m$$, we calculate:
Visible text: For each row , we calculate:
```math
y_i := \frac{1}{l_{ii}} \cdot \left( c_i - \sum_{j=1}^{i-1} l_{ij} \cdot y_j \right)
```
This algorithm is efficient because we only need to calculate one value at each step. The diagonal element $$l_{ii}$$ is always $$1$$ for matrix $$L$$ from LU decomposition, so division becomes simple.
Visible text: This algorithm is efficient because we only need to calculate one value at each step. The diagonal element is always for matrix from LU decomposition, so division becomes simple.
### Backward Substitution
The backward substitution algorithm solves $$U \cdot x = y$$ with matrix $$U$$ in row echelon form with $$r$$ steps on columns $$j_1, \ldots, j_r$$ and vector $$d = y$$.
Visible text: The backward substitution algorithm solves with matrix in row echelon form with steps on columns and vector .
First, we check whether the system can be solved. If $$r < m$$ and there exists $$d_i \neq 0$$ for $$i \in \{r + 1, \ldots, m\}$$, then the system has no solution.
Visible text: First, we check whether the system can be solved. If and there exists for , then the system has no solution.
If the system can be solved, we initialize kernel matrix $$K$$ with size $$n \times (n-r)$$, and start with $$k = 0$$ and $$i = r$$.
Visible text: If the system can be solved, we initialize kernel matrix with size , and start with and .
The algorithm works backward from column $$j = n$$ to $$j = 1$$. For each column, we check whether that column is a pivot step or not.
Visible text: The algorithm works backward from column to . For each column, we check whether that column is a pivot step or not.
1. **If column j is a pivot step**, meaning $$j = j_i$$ for step $$i$$, then we calculate the solution for variable $$x_j$$:
```math
x_j := \frac{1}{u_{ij}} \cdot \left( d_i - \sum_{l=j+1}^n u_{il} \cdot x_l \right)
```
And we also calculate the contribution of this variable to the kernel matrix:
```math
K_{jq} := \frac{1}{u_{ij}} \cdot \left( - \sum_{l=j+1}^n u_{il} \cdot K_{lq} \right)
```
for $$q = 1, \ldots, n - r$$, then we decrease $$i$$ by 1.
2. **If column j is not a pivot step**, meaning there is no pivot in that column, then variable $$x_j$$ is a free variable. We set:
- $$k := k + 1$$
- $$x_j := 0$$ (particular solution)
- $$K_{jk} := 1$$ (kernel basis)
Visible text: 1. **If column j is a pivot step**, meaning for step , then we calculate the solution for variable :
And we also calculate the contribution of this variable to the kernel matrix:
for , then we decrease by 1.
2. **If column j is not a pivot step**, meaning there is no pivot in that column, then variable is a free variable. We set:
-
- (particular solution)
- (kernel basis)
This algorithm produces solution $$x$$ and matrix $$K$$ that satisfy $$U \cdot x = d$$ and $$U \cdot K = 0$$.
Visible text: This algorithm produces solution and matrix that satisfy and .
The columns of matrix $$K$$ form a basis for the kernel (null space) of matrix $$U$$. This means that each column $$K$$ is a vector that when multiplied by $$U$$ produces a zero vector.
Visible text: The columns of matrix form a basis for the kernel (null space) of matrix . This means that each column is a vector that when multiplied by produces a zero vector.
The complete solution set is:
```math
\left\{ x + K \cdot \begin{pmatrix} t_1 \\ \vdots \\ t_{n-r} \end{pmatrix} : t_1, \ldots, t_{n-r} \in \mathbb{R} \right\}
```
## Application Example
Let's see how LU decomposition is used to solve linear systems concretely. Suppose we have vector $$b = \begin{pmatrix} 4 \\ 3 \\ 3 \end{pmatrix}$$ and want to solve system $$A \cdot x = b$$.
Visible text: Let's see how LU decomposition is used to solve linear systems concretely. Suppose we have vector and want to solve system .
1. **The first step** is to calculate $$c = P \cdot b$$. Since permutation matrix $$P$$ swaps the first row with the second row, we get:
```math
c = P \cdot b = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 4 \\ 3 \\ 3 \end{pmatrix} = \begin{pmatrix} 3 \\ 4 \\ 3 \end{pmatrix}
```
2. **The second step** is forward substitution to solve $$L \cdot y = c$$. We use the lower triangular matrix $$L$$ that we have obtained. This process is done from top to bottom like filling stairs one by one.
For the first row with $$l_{11} = 1$$, we get $$y_1 = c_1 = 3$$.
For the second row with $$l_{22} = 1$$, we get $$y_2 = c_2 = 4$$.
For the third row with $$l_{33} = 1$$, we calculate:
```math
y_3 = c_3 - l_{31} \cdot y_1 - l_{32} \cdot y_2
```
```math
= 3 - \frac{1}{3} \cdot 3 - \frac{1}{2} \cdot 4
```
```math
= 3 - 1 - 2 = 0
```
So we obtain:
```math
y = L^{-1} \cdot c = \begin{pmatrix} 3 \\ 4 \\ 0 \end{pmatrix}
```
3. **The third step** is backward substitution to solve $$U \cdot x = y$$. Matrix $$U$$ in row echelon form has pivots in columns $$1, 2, 4$$. Column $$3$$ has no pivot so it becomes a free variable.
From the backward substitution process, we obtain the particular solution:
```math
x = \begin{pmatrix} -1 \\ 2 \\ 0 \\ 0 \end{pmatrix}
```
and kernel matrix:
```math
K = \begin{pmatrix} 1 \\ -2 \\ 1 \\ 0 \end{pmatrix}
```
4. **Result interpretation** is very important to understand. Solution $$x$$ is one particular solution of the equation system. Matrix $$K$$ shows the direction in which we can move in the solution space without changing the result of multiplication $$A \cdot x$$.
5. **Mathematical verification** shows that for any parameter value $$t_1$$, the general solution $$x + K \cdot t_1$$ still satisfies the original equation.
Let's prove by calculating $$A \cdot (x + K \cdot t_1)$$. The result is:
```math
A \cdot (x + K \cdot t_1) = A \cdot \begin{pmatrix} -1 \\ 2 \\ 0 \\ 0 \end{pmatrix} + A \cdot \begin{pmatrix} 1 \\ -2 \\ 1 \\ 0 \end{pmatrix} \cdot t_1
```
```math
= \begin{pmatrix} 4 \\ 3 \\ 3 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \cdot t_1 = b
```
Note that $$A \cdot K = 0$$, which means vector $$K$$ is in the null space of matrix $$A$$. This explains why adding multiples of $$K$$ to the solution does not change the result.
Visible text: 1. **The first step** is to calculate . Since permutation matrix swaps the first row with the second row, we get:
2. **The second step** is forward substitution to solve . We use the lower triangular matrix that we have obtained. This process is done from top to bottom like filling stairs one by one.
For the first row with , we get .
For the second row with , we get .
For the third row with , we calculate:
So we obtain:
3. **The third step** is backward substitution to solve . Matrix in row echelon form has pivots in columns . Column has no pivot so it becomes a free variable.
From the backward substitution process, we obtain the particular solution:
and kernel matrix:
4. **Result interpretation** is very important to understand. Solution is one particular solution of the equation system. Matrix shows the direction in which we can move in the solution space without changing the result of multiplication .
5. **Mathematical verification** shows that for any parameter value , the general solution still satisfies the original equation.
Let's prove by calculating . The result is:
Note that , which means vector is in the null space of matrix . This explains why adding multiples of to the solution does not change the result.
Compared to computing LU decomposition completely, this forward and backward substitution process is much more efficient for solving linear systems with different right-hand side vectors.