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)