LU Decomposition

LU decomposition is a common method to solve polynomial \(Ax=b\), and this method is consist of three steps:

1. LU factorization: Gaussian elimination on matrix \(A\), factorize \(A\) to \(L,U\) such that \(PA=LU\)

2. Forward substitution: Solve lower triangular system \(Ld = Pb\)

3. Backward substitution: Solve upper triangular system \(Ux = d\)

In the following part, I would use the example of \(A\) and \(b\) to explain the calculation of LU decomposition. \[ A=\left(\begin{array}{ccc} 1 & 2 & 3 \\ 2 & 5 & 10 \\ 3 & 10 & 10 \end{array}\right), \quad b=\left(\begin{array}{c} 3 \\ 7 \\ 13 \end{array}\right) \]

(1) LU factorization

In LU factorization, we need to factorize \(A\) into \(L\) and \(U\), and let \(LU=PA\) \[ L=\left(\begin{array}{ccc} a_1 & 0 & 0 \\ a_2 & a_3 & 0 \\ a_4 & a_5 & a_6 \end{array}\right), \quad U=\left(\begin{array}{ccc} b_1 & b_2 & b_3 \\ 0 & b_4 & b_5 \\ 0 & 0 & b_6 \end{array}\right) \]

Sometimes, we can get the right result through the decomposition of \(LU=A\), but this decomposition isn't stable that might cause a large error in the calculation. Therefore, we need to do the "pivoting" procedure (\(P\)) during the factorization step. "pivoting" means that we should keep the value of each non-zero element of each row to be larger the other elements of the same column below this one, and we can achieve this goal by row exchange.

\(U\) is obtained from Gaussian elimination of \(A\), and \(L\) is computed from \(U\). \(P\) is the identity matrix at first, and if we exchange row during the Gaussian elimination, the same exchange need to be applied on \(P\). And \(A, P, L, U\) all have the same size.

The calculation below is LU factorization example:

\[ \begin{aligned} &A=\left(\begin{array}{ccc} 1 & 2 & 3 \\ 2 & 5 & 10 \\ 3 & 10 & 16 \end{array}\right), \quad P=\left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right) \\ &\stackrel{R_{1} \leftrightarrow R_{3}}{\longrightarrow}\left(\begin{array}{ccc} 3 & 10 & 16 \\ 2 & 5 & 10 \\ 1 & 2 & 3 \end{array}\right), \quad P=\left(\begin{array}{lll} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{array}\right) \\ &\stackrel{R_{2}-\frac{2}{3} R_{1}}{\longrightarrow}\left(\begin{array}{ccc} 3 & 10 & 10 \\ 0 & -\frac{5}{3} & -\frac{2}{3} \\ 0 & -\frac{4}{3} & -\frac{1}{3} \end{array}\right), \quad P=\left(\begin{array}{lll} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{array}\right) \\ &\stackrel{R_{3}-\frac{1}{5} R_{1} R_{2}}{\longrightarrow}\left(\begin{array}{ccc} 3 & 10 & 16 \\ 0 & -\frac{5}{3} & -\frac{2}{3} \\ 0 & 0 & -\frac{5}{5} \end{array}\right)=U, \quad P=\left(\begin{array}{lll} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{array}\right) \end{aligned} \]

\[ \begin{aligned} L&=\left(\left(\begin{array}{l} 3 \\ 2 \\ 1 \end{array}\right) / 3,\left(\begin{array}{c} 0 \\ -\frac{5}{3} \\ -\frac{4}{3} \end{array}\right) / -\frac{5}{3},\left(\begin{array}{c} 0 \\ 0 \\ -\frac{9}{5} \end{array}\right) /-\frac{9}{5}\right) \\ &=\left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{2}{3} & 1 & 0 \\ \frac{1}{3} & \frac{4}{5} & 1 \end{array}\right) \end{aligned} \]

(2) Forward and Backward substitution

Forward and Backward substitution is used to calculate the \(x\) and \(y\) in the following functions. \[ \begin{aligned} &\left\{\begin{array} { l } { A x = b } \\ { L U = P A } \\ { U x = y } \end{array} \rightarrow \left\{\begin{array}{l} L y=P b \\ U x=y \end{array}\right.\right.\\ \end{aligned} \] Forward substitution: \[ \begin{aligned} &P b=\left(\begin{array}{lll} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{array}\right) \cdot\left(\begin{array}{l} 3 \\ 7 \\ 13 \end{array}\right)=\left(\begin{array}{c} 13 \\ 7 \\ 3 \end{array}\right)\\ &L y=\left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{2}{3} & 1 & 0 \\ \frac{1}{3} & \frac{4}{5} & 1 \end{array}\right) \cdot\left(\begin{array}{l} y_{1} \\ y_{2} \\ y_{3} \end{array}\right)=\left(\begin{array}{c} 13 \\ 7 \\ 3 \end{array}\right)\\ &\rightarrow\left\{\begin{array} { l } { y _ { 1 } = 1 3 } \\ { \frac { 2 } { 3 } y _ { 1 } + y _ { 2 } = 1 } \\ { \frac { 1 } { 3 } y _ { 1 } + \frac { 4 } { 5 } y _ { 2 } + y _ { 3 } = 3 } \end{array} \rightarrow \left\{\begin{array}{l} y_{1}=13 \\ y_{2}=-\frac{5}{3} \\ y_{3}=0 \end{array} \quad \rightarrow y=\left(\begin{array}{c} 13 \\ -\frac{5}{3} \\ 0 \end{array}\right)\right.\right.\\ \end{aligned} \] Backward substitution: \[ \begin{aligned} &Ux=\left(\begin{array}{ccc} 3 & 10 & 16 \\ 0 & -\frac{5}{3} & -\frac{2}{3} \\ 0 & 0 & -\frac{9}{5} \end{array}\right) \cdot\left(\begin{array}{l} x_{1} \\ x_{2} \\ x_{3} \end{array}\right)=\left(\begin{array}{c} 13 \\ -\frac{5}{3} \\ 0 \end{array}\right)\\ &\rightarrow\left\{\begin{array} { c } { 3 x _ { 1 } + 1 0 x _ { 2 } + 1 6 x _ { 3 } = 1 3 } \\ { - \frac { 5 } { 3 } x _ { 2 } - \frac { 2 } { 3 } x _ { 3 } = - \frac { 5 } { 3 } } \\ { - \frac { 9 } { 5 } x _ { 3 } = 0 } \end{array} \rightarrow \left\{\begin{array}{l} x_{1}=1 \\ x_{2}=1 \\ x_{3}=0 \end{array} \rightarrow x=\left(\begin{array}{l} 1 \\ 1 \\ 0 \end{array}\right)\right.\right. \end{aligned} \]

(3) LU decomposition in Python

import numpy as np
from scipy.linalg import lu
from scipy.linalg import solve_triangular
A=np.array([[3 -1 2],[1 0 -1],[4 2 -3]])
P, L, U =lu(A)
b=np.array([[8],[-1],[-4]])
b = P.T @ b # Row changes in b
# Forward and backward substitution
d=solve_triangular(L, b, lower=True)
x=solve_triangular(U, y)
# Check if the result from LU-decomposition is the same as direct computation
x2 = np.linalg.solve(A,b)