LU Decomposition

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:

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:

with $i = 1, \ldots, r$, row exchange matrix $P_i$ or $P_i = I$ and elimination matrix:

The inverse of $L_i$ is:

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$.

Structure of Decomposition Results

From Gaussian elimination with steps $j_1, \ldots, j_r$, we get three matrices with special structure:

Matrix $P$ is obtained from applying permutation vector $p$ 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$

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.

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

   $$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.

Final result:

From this elimination result, we obtain:

The proof shows that $L \cdot U = P \cdot A$ for the initially given matrix $A$.

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}$.

Algorithm Implementation in Python

Here is the implementation of the LU decomposition algorithm in Python:

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

Solving Linear Systems

After obtaining LU decomposition, solving the system $A \cdot x = b$ 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:

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.

For each row $i = 1, \ldots, m$, we calculate:

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.

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$.

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.

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$.

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.

1. If column j is a pivot step, meaning $j = j_i$ for step $i$, then we calculate the solution for variable $x_j$:

   $$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:

   $$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)

This algorithm produces solution $x$ and matrix $K$ that satisfy $U \cdot x = d$ and $U \cdot 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.

The complete solution set is:

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$.

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:

   $$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:

   $$y_3 = c_3 - l_{31} \cdot y_1 - l_{32} \cdot y_2$$

   $$= 3 - \frac{1}{3} \cdot 3 - \frac{1}{2} \cdot 4$$

   $$= 3 - 1 - 2 = 0$$

   So we obtain:

   $$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, and 4. Column 3 has no pivot so it becomes a free variable.

   From the backward substitution process, we obtain the particular solution:

   $$x = \begin{pmatrix} -1 \\ 2 \\ 0 \\ 0 \end{pmatrix}$$

   and kernel matrix:

   $$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:

   $$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$$

   $$= \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.

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.