LU factorization

factorizing by $A=LU$ where $L$ is lower triangular and $U$ is upper triangular matrix.

Basic factorization method

Make all diagonal elements of $L$ to one and manually calculate it with $U$.

3x3 example

$$L=\left[\begin{array}{ccc}1& 0& 0\\ x& 1& 0\\ y& z& 1\end{array}\right]$$ $$U=\left[\begin{array}{ccc}a& b& c\\ 0& d& e\\ 0& 0& f\end{array}\right]$$

So we can say

$$A=\left[\begin{array}{ccc}a& b& c\\ xa& xb+d& xc+e\\ ya& yb+zd& yc+ze+f\end{array}\right]$$

$\square$

Multiplying factorization method (like Gaussian elimination?)

Started with $A=IA$ where $I$ is identical matrix.

3x3 example

\begin{eqs} A &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ \frac{d}{a} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a & b & c \\ 0 & e - \frac{bd}{a} & f - \frac{cd}{a} \\ g & h & i \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ \frac{d}{a} & 1 & 0 \\ \frac{g}{a} & 0 & 1 \end{bmatrix} \begin{bmatrix} a & b & c \\ 0 & e - \frac{bd}{a} & f - \frac{cd}{a} \\ 0 & h - \frac{bg}{a} & i - \frac{cg}{a} \end{bmatrix} \end{eqs}

then eliminate $h-\frac{bg}{a}$ term by first or second row etc. $\square$

Partial Pivot

Exchange rows to let the biggest element in first row.
Given by the formula $PA=LU$.
But there are types that no need to pivot at all:

• Diagonally dominant
• SPD, where we can use Cholesky factorization instead

Applications

After factorization, we can solve $Ax=b$ by Forward-backward substitution.
Basically $Ax=b$ is equivalent to $x=A^{-1}b$.

Sherman-Morrison

$x=(A-u{v}^T{)}^{-1}b=A^{-1}b+A^{-1}u(1-v^T{A}^{-1}u{)}^{-1}{v}^T{A}^{-1}b$
We can solve those $A^{-1}b$ and $A^{-1}u$ by Forward-backward substitution.

Costs

generally total cost in $O(n^3)\approx \frac{2n^3}{3}$ operations
But after enable partial pivoting (only row exchange), total cost is $O(n^2)<<\frac{2n^3}{3}$
However for full pivoting, it is still $O(n^3)$ for additional searching cost

Code

credit: UTexas.

import numpy as np
 
def LU(A):
	n = A.shape[0]
	U = A.copy()
	L = np.eye(n, dtype=np.double)
	for i in range(n):
		factor = U[i + 1:, i] / U[i, i]
		L[i + 1:, i] = factor
		U[i + 1:] -= factor[:, np.newaxis] * U[i]
	return L, U