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 , we get:
A
P⋅A=L⋅U
Where P is a permutation matrix (row order arranger), L is a lower triangular matrix, and U is an upper triangular matrix.
Gaussian elimination computes this LU decomposition through a series of systematic row operations. This process can be written as:
U=Lr⋅Pr⋅…⋅L3⋅P3⋅L2⋅P2⋅L1⋅P1⋅A
with i=1,…,r, row exchange matrix Pi or Pi=I and elimination matrix:
Li=1⋱λi+1i⋮λmi⋱1
The inverse of Li is:
Li−1=1⋱−λi+1i⋮−λmi⋱1
Since Pi−1=Pi, we can rewrite this sequence of operations into a simpler form until we finally get L−1⋅P⋅A=U.
The Gaussian elimination algorithm works like systematic cleaning. For each column, we perform:
Find the best pivot element
Swap rows if necessary
Eliminate elements below the pivot
Store multiplier information
The time complexity for computing LU decomposition is 31n3+O(n2) arithmetic operations. The proof shows that multiplication and addition occur simultaneously in the elimination row transformation aij:=aij+λir⋅arl.
Here is the implementation of the LU decomposition algorithm in Python:
lu_decomposition.py
import numpy as npdef 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, Pdef 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 ydef 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 xdef 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
The forward substitution algorithm solves L⋅y=P⋅b with regular lower triangular matrix L and vector c=P⋅b. This process is like filling stairs from bottom to top, where each step depends on the previous step.
For each row i=1,…,m, we calculate:
yi:=lii1⋅(ci−j=1∑i−1lij⋅yj)
This algorithm is efficient because we only need to calculate one value at each step. The diagonal element lii is always 1 for matrix L from LU decomposition, so division becomes simple.
The backward substitution algorithm solves U⋅x=y with matrix U in row echelon form with r steps on columns j1,…,jr and vector d=y.
First, we check whether the system can be solved. If r<m and there exists di=0 for i∈{r+1,…,m}, then the system has no solution.
If the system can be solved, we initialize kernel matrix K with size n×(n−r), and start with k=0 and i=r.
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.
If column j is a pivot step, meaning j=ji for step i, then we calculate the solution for variable xj:
xj:=uij1⋅di−l=j+1∑nuil⋅xl
And we also calculate the contribution of this variable to the kernel matrix:
Kjq:=uij1⋅−l=j+1∑nuil⋅Klq
for q=1,…,n−r, then we decrease i by 1.
If column j is not a pivot step, meaning there is no pivot in that column, then variable xj is a free variable. We set:
k:=k+1
xj:=0 (particular solution)
This algorithm produces solution x and matrix K that satisfy U⋅x=d and U⋅K=0.
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.
Let's see how LU decomposition is used to solve linear systems concretely. Suppose we have vector b=433 and want to solve system A⋅x=b.
The first step is to calculate c=P⋅b. Since permutation matrix P swaps the first row with the second row, we get:
c=P⋅b=010100001433=343
The second step is forward substitution to solve L⋅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 l11=1, we get .
The third step is backward substitution to solve U⋅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.
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⋅x.
Mathematical verification shows that for any parameter value t1, the general solution x+K⋅t1 still satisfies the original equation.
Let's prove by calculating A⋅(x+K⋅t1). The result is:
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.
From the backward substitution process, we obtain the particular solution:
x=−1200
and kernel matrix:
K=1−210
A⋅(x+K⋅t1)=A⋅−1200+A⋅1−210⋅t1
=433+000⋅t1=b
Note that A⋅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.
