封面：Epsilons, no. 3: The LU decomposition - by Tivadar Danka

Definition

考虑A∈Rn×nA\in\mathbb{R}^{n\times n}，满足所有顺序主子式（leading principal minors）不为0，则存在以下分解：

A=LUA = LU

其中，L∈Rn×nL\in\mathbb R^{n\times n}是一个单位下三角矩阵，U∈Rn×nU\in\mathbb{R}^{n\times n}是一个上三角阵

L=[10l1,21⋮⋮⋱ln,1ln,2⋯1]U=[u1,1⋯u1,n−1u1,n⋱⋮⋮un−1,n−1un−1,n0un,n]L = \begin{bmatrix}1&&&0\\l_{1,2}&1&&\\\vdots&\vdots&\ddots\\ l_{n,1}&l_{n,2}&\cdots&1\end{bmatrix}\quad U = \begin{bmatrix}u_{1,1}&\cdots&u_{1,n-1}&u_{1,n}\\&\ddots&\vdots&\vdots\\&&u_{n-1,n-1}&u_{n-1,n}\\ 0&&&u_{n,n}\end{bmatrix}

Application

1. 解线性方程

对于线性方程组Ax=bA\mathbf x=b, A∈Rn×n,x∈Rn,b∈RnA\in\mathbb R^{n\times n},\mathbf x\in\mathbb R^n,b\in\mathbb R^{n}. 如果有A=LUA = LU，则

Ax=b⇔LUx=b⇔{Ux=yLy=bA\mathbf x=b\Leftrightarrow LU\mathbf x=b\Leftrightarrow \left\{\begin{aligned}U\mathbf x &= \mathbf{y}\\L\mathbf y&=b\end{aligned}\right.

由于LL是一个下三角阵，可以轻松的得到Ly=bL\mathbf y=b的解

Ly=b⇔[1l211l31l321⋮⋮⋱⋱ln1ln2⋯lnn−11][y1y2y3⋮yn]=[b1b2b3⋮bn]L\mathbf y=b\Leftrightarrow\begin{bmatrix}1&&&&\\l_{21}&1&&&\\l_{31}&l_{32}&1&&\\\vdots&\vdots&\ddots&\ddots\\l_{n1}&l_{n2}&\cdots&l_{nn-1}&1\end{bmatrix}\begin{bmatrix}y_{1}\\y_{2}\\y_{3}\\\vdots\\y_{n}\end{bmatrix}=\begin{bmatrix}b_{1}\\b_{2}\\b_{3}\\\vdots\\b_{n}\end{bmatrix}

按照y1→yny_1\rightarrow y_n的顺序：

y1=b1l21y1+y2=b2⋮ln1y1+⋯+lnn−1yn−1+yn=bn\begin{gather*} y_1 = b_1 \\ l_{21}y_1 + y_2 = b_2 \\ \vdots \\ l_{n1}y_1 + \cdots + l_{nn-1}y_{n-1} + y_n = b_n \end{gather*}

对于yky_k，y1:k−1y_{1:k-1}都是已知的，每个方程只有yky_k一个未知数。

由于UU是一个上三角矩阵，同样的，可以很轻松的得到Ux=yU\mathbf x=\mathbf y的解

ux=y⇔[u1,1⋯⋯u1,n⋱⋱⋮un−1,n−1un−1,nun,n][x1⋮xn−1xn]=[y1⋮yn−1yn]u_{x}=y\Leftrightarrow\begin{bmatrix}u_{1,1}&\cdots&\cdots&u_{1,n}\\&\ddots&\ddots& \vdots\\&& u_{n-1,n-1}&u_{n-1,n}\\&&& u_{n,n}\end{bmatrix}\begin{bmatrix}x_{1}\\\vdots\\x_{n-1}\\x_{n}\end{bmatrix}=\begin{bmatrix}y_{1}\\\vdots\\y_{n-1}\\y_{n}\end{bmatrix}

从xn→x1x_n\rightarrow x_1的顺序。

2. 计算行列式

det⁡(A)=det⁡(LU)=det⁡(L)det⁡(U)=∏i=1nuii\det(A) = \det(LU) = \det(L)\det(U)=\prod_{i=1}^n u_{ii}

Algorithm

思路：通过一系列简单的矩阵变换，逐步将AA转化为上三角矩阵UU。在每一步中，使用一个单位下三角矩阵LjL_j来消去AA中的部分元素，直到得到最终的UU矩阵。

A=LU⇔L−1A=U⇔Ln−1⋯L2L1A=UA = LU\Leftrightarrow L^{-1}A = U\Leftrightarrow L_{n-1}\cdots L_2L_1A = U

其中，L1:n−1L_{1:n-1}都是简单的单位下三角阵

令

Lj=I−liej⊤L_j = I - l_ie_j^\top

其中，ljl_j满足lij=0,i≤jl_{ij}=0, i\leq j

lj=[0⋮0lj+1,j⋮ln,j],ej⊤=[0⋯01j↓ 0⋯0]l_j = \begin{bmatrix}0\\\vdots\\0\\ l_{j+1,j}\\\vdots\\l_{n,j}\end{bmatrix},\quad e_j^\top = \begin{bmatrix}0&\cdots&0& \overset{\smash{\begin{aligned}j\\\downarrow\\~\end{aligned}}}{1}&0&\cdots0\end{bmatrix}

Lj=I−ljej⊤=[1⋱11−lj+1,j⋱⋮⋱−ln,j1]L_j = I - l_j e_j^\top = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & 1 & & & \\ & & & 1 & & \\ & & & -l_{j+1,j} & \ddots & \\ & & & \vdots & & \ddots \\ & & & -l_{n,j} & & & 1 \\ \end{bmatrix}

令aia_i和a~i\tilde a_i分别表示AA的第ii列和第ii行

A=[a1⋯an]=[a~1⋮a~n]A = \begin{bmatrix}a_1&\cdots&a_n\end{bmatrix}= \begin{bmatrix}\tilde a_1\\\vdots\\\tilde a_n\end{bmatrix}

希望经过L1L_1的变换之后，a1a_1这一列可以变成一个上三角阵的形状（a1,2:n=0a_{1,2:n}=0）

L1a1=c1e1, for some c1∈RL_1a_1=c_1e_1,\quad \text{ for some } c_1\in\mathbb R

(I−l1e1⊤)a1=c1e1⇒a1−l1(e1⊤a1)=c1e1⇔a1−a11l1=c1e1⇔[a11a21⋮an1]−[0a11l21⋮a11ln1]=[c10⋮0]\begin{gather*} (I - l_1e_1^\top)a_1=c_1e_1\Rightarrow a_1-l_1(e_1^\top a_1)=c_1e_1\\ \Leftrightarrow a_1-a_{11}l_1=c_1e_1\\ \Leftrightarrow \begin{bmatrix}a_{11}\\a_{21}\\\vdots\\ a_{n1}\end{bmatrix}-\begin{bmatrix}0\\a_{11}l_{21}\\\vdots\\ a_{11}l_{n1}\end{bmatrix} = \begin{bmatrix}c_1\\0\\\vdots\\ 0\end{bmatrix}\\ \end{gather*}

选择

l1=[0a21/a11⋮an1/a11]l_1 = \begin{bmatrix}0\\a_{21}/a_{11}\\\vdots\\ a_{n1}/a_{11}\end{bmatrix}

A(1)=L1A=(I−l1e1⊤)A=A−l1a~1=A−[0l21⋮ln1][a11a12⋯a1n]=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann]−[00⋯0l21a11⋮li1⋅a1jln1a11]=[a11a12⋯a1n0⋮aij−li1a1j0]\begin{align} A^{(1)}&=L_1A = (I-l_1e_1^\top)A = A-l_1\tilde a_1\\ &=A - \begin{bmatrix}0\\l_{21}\\\vdots\\ l_{n1}\end{bmatrix}\begin{bmatrix}a_{11}&a_{12}&\cdots& a_{1n}\end{bmatrix}\\ &=\begin{bmatrix} \begin{array}{c|ccc} a_{11} & a_{12} & \cdots & a_{1n} \\ \hline a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots&\vdots&\ddots&\vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{array} \end{bmatrix}- \begin{bmatrix} \begin{array}{c|ccc} 0 & 0 & \cdots & 0 \\ \hline l_{21}a_{11} & & & \\ \vdots&&l_{i1}\cdot a_{1j}&\\ l_{n1}a_{11} & & & \\ \end{array} \end{bmatrix} \\ &=\begin{bmatrix} \begin{array}{c|ccc} a_{11} & a_{12} & \cdots & a_{1n} \\ \hline 0 & & & \\ \vdots&&a_{ij}-l_{i1}a_{1j}&\\ 0 & & & \\ \end{array} \end{bmatrix} \end{align}

对于L2L_2，类似的，只需要解一个n−1n-1维的子问题

对于LkL_k和Ak−1A_{k-1}

lik=aik(k−1)/akk(k)aij(k)=aij(k−1)−likakj(k−1)\begin{align} l_{ik}&=a^{(k-1)}_{ik}/a_{kk}^{(k)}\\ a_{ij}^{(k)} &= a_{ij}^{(k-1)}-l_{ik}a_{kj}^{(k-1)} \end{align}

image-20250216233101038

重复以上步骤直到k=n−1k=n-1，最终，有：

Ln−1Ln−1⋯L2L1A=A(n−1)≡UL_{n-1}L_{n-1}\cdots L_2L_1A = A^{(n-1)}\equiv U

L=(Ln−1Ln−1⋯L2L1)−1L =(L_{n-1}L_{n-1}\cdots L_2L_1)^{-1}

Theorem:

L=I+l1e1⊤+⋯+ln−1en−1⊤L = I+l_1e_1^\top+\cdots+l_{n-1}e_{n-1}^\top

proof:

L=L1−1⋯Ln−1−1L = L_1^{-1}\cdots L_{n-1}^{-1}

首先证明Lj−1=I+ljej⊤L_{j}^{-1} = I+l_je_j^\top，recall that Lj=I−ljej⊤L_j = I-l_je_j^\top

(I+ljej⊤)(I−ljej⊤)=I−ljej⊤ljej⊤=I−lj(ej⊤lj)ej=I−0=I\begin{align} (I+l_je_j^\top)(I-l_je_j^\top) &= I-l_je_j^\top l_je_j^\top\\ &=I-l_j(e_j^\top l_j)e_j\\ &=I-0=I \end{align}

根据定义 ljl_j的第 jj位是 0，而 eje_{j}只有第 jj位是1，因此内积为0

然后使用归纳法：

L1−1L2−1=(I+l1e1⊤)(I+l2e2⊤)=I+l1e1⊤+l2e2⊤+l1e1⊤l2e2⊤=I+l1e1⊤+l2e2⊤+l1(e1⊤l2)e2⊤=I+l1e1⊤+l2e2⊤\begin{align} L_1^{-1}L_2^{-1}&=(I+l_1e_1^\top)(I+l_2e_2^\top)\\ &=I+l_1e_1^\top+l_2e_2^\top+l_1e_1^\top l_2e_2^\top\\ &=I+l_1e_1^\top+l_2e_2^\top+l_1(e_1^\top l_2)e_2^\top\\ &=I+l_1e_1^\top+l_2e_2^\top \end{align}

和前面同样的道理

如果L1−1⋯Lk−1=I+l1e1⊤+⋯+lkek⊤L_{1}^{-1}\cdots L_{k}^{-1} = I+l_1e_1^\top+\cdots+l_ke_k^\top，则

L1−1⋯Lk+1−1=(I+l1e1⊤+⋯+lkek⊤)(I+lk+1ek+1⊤)=I+l1e1⊤+⋯+lkek⊤+lk+1ek+1⊤+l1(e1⊤lk+1)ek+1⊤+⋯+lk(ek⊤lk+1)ek+1⊤=I+l1e1⊤+⋯+lkek⊤+lk+1ek+1⊤\begin{align} L_{1}^{-1}\cdots L_{k+1}^{-1} &=(I+l_1e_1^\top+\cdots+l_ke_k^\top)(I+l_{k+1}e_{k+1}^\top)\\ &=I+l_1e_1^\top+\cdots+l_ke_k^\top+l_{k+1}e_{k+1}^\top+l_1(e_1^\top l_{k+1})e_{k+1}^\top+\cdots+l_k(e_k^\top l_{k+1})e_{k+1}^\top\\ &=I+l_1e_1^\top+\cdots+l_ke_k^\top+l_{k+1}e_{k+1}^\top \end{align}

检查LL是不是单位上三角阵：

[L]ij=ei⊤Lej=ei⊤(I+l1e1⊤+⋯+ln−1en−1⊤)ej=Iij+ei⊤lj(ei⊤ej=0 if i≠j)={1+0=1if i=j0+0=0if i<j0+lij=lijif i>j\begin{align} [L]_{ij} &= e_i^\top Le_j=e_i^\top (I+l_1e_1^\top+\cdots+l_{n-1}e_{n-1}^\top)e_j\\ &=I_{ij}+e_i^\top l_j\quad(e_i^\top e_j=0\text{ if }i\neq j)\\ &=\left\{\begin{aligned} &1+0=1&&\text{if }i=j\\ &0+0=0&&\text{if }i<j\\ &0+l_{ij}=l_{ij}&&\text{if }i>j \end{aligned}\right. \end{align}

因此

L=[1l211l31l321⋮⋮⋱⋱ln1ln2⋯lnn−11]L =\begin{bmatrix}1&&&&\\l_{21}&1&&&\\l_{31}&l_{32}&1&&\\\vdots&\vdots&\ddots&\ddots\\l_{n1}&l_{n2}&\cdots&l_{nn-1}&1\end{bmatrix}

最终算法

image-20250217001329032

image-20250217000955585

计算复杂度

∑k=1n−1((n−k)+2(n−k)2)=∑k′=1n−1(k′+2(k′)2)=O(n3)\sum_{k=1}^{n-1}\big((n-k)+2(n-k)^2\big)=\sum_{k^\prime=1}^{n-1}\big(k^\prime+2(k^\prime)^2\big) = \mathcal O(n^3)

Example

A=[21−145−3−25−2]A=\begin{bmatrix}2&1&-1\\4&5&-3\\-2&5&-2\end{bmatrix}

LU example

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np


def LU(A):
    assert A.shape[0] == A.shape[1]
    for k in range(A.shape[0] - 1):
        l = np.zeros((A.shape[1], 1))
        # l_{ik} = a_{ik}/a{kk}, i\in[k+1,n]
        l[k + 1:, 0] = A[k + 1:, k] / A[k, k]
        # A_{ij} = a_{ij} - l_{ik}a_{kj}
        A[k + 1:, k + 1:] = A[k + 1:, k + 1:] - l[k + 1:, :] @ A[k:k + 1, k + 1:]
        # use the lower triangular part of A to save l
        A[k + 1:, k:k + 1] = l[k + 1:, :]
    L = np.tril(A, k=-1) + np.eye(A.shape[0])
    U = np.triu(A, k=0)
    return L, U


A = np.array([[2, 1, -1], [4, 5, -3], [-2, 5, -2]])
L, U = LU(A)
print(L @ U)

变体

Pivoted LU decomposition

在LU分解中，有可能遇到akk=0a_{kk}=0的情况，此时算法无法继续。Pivoted LU decomposition 引入主元，即n−k+1n-k+1矩阵左上角的元素。每次通过行交换选择最大的元素作为主元。使用矩阵PP记录行交换过程，即

PA=LUPA = LU

最终再交换回去。

同理可以进行列交换

PAQ=LUPAQ = LU

Cholesky分解

Theorem:

如果A∈Rn×nA\in\mathbb R^{n\times n}是一个半正定矩阵，则存在一个分解：

A=LL⊤A = LL^\top

其中LL是一个下三角矩阵（不一定是单位下三角）

proof:

因为AA时一个半正定矩阵，因此所有的顺序主子式都可逆。所以存在一个LULU使得A=L0UA=L_0U

令D=diag(U)D={\rm diag}(U)，DD可逆（我们假设主元不为0，这样的话UU的对角线元素中也不会有0）

A=L0U=LD(D−1U)=L0DU0A = L_0U=LD(D^{-1}U)=L_0DU_0

其中U0U_0是一个单位上三角阵。且U0=L0⊤U_0=L_0^\top，因为：

A⊤=A⇒U0⊤DL0=L0DU0L0−1U0⊤D=DU0L0−⊤\begin{gather} A^\top = A \Rightarrow U_0^\top DL_0 = L_0DU_0\\ L_0^{-1}U_0^\top D = DU_0L_0^{-\top} \end{gather}

其中，L0−1U0⊤L_0^{-1}U_0^\top是一个下三角阵，U0L0−⊤U_0L_0^{-\top}是一个上三角阵

上三角阵的逆是上三角阵，下三角阵的逆是下三角阵

因此L0−1U0⊤L_0^{-1}U_0^\top和U0L0−⊤U_0L_0^{-\top}是对角阵（一个上三角阵等于一个下三角阵）。因此L0−1U0⊤=U0L0−⊤=I⇒U0=L0⊤L_0^{-1}U_0^\top=U_0L_0^{-\top}=I\Rightarrow U_0=L_0^\top

因此

A=L0DL0⊤A = L_0DL_0^\top

接下来证明dii>0d_{ii}>0:

令x=L0−⊤ei≠0x=L_0^{-\top}e_i\neq 0

x⊤Ax=ei⊤L0−1(L0DL0⊤)L−⊤ei=ei⊤Dei=dii>0x^\top Ax =e_i^\top L_0^{-1}(L_0DL_0^\top)L^{-\top}e_i=e^\top_iDe_i = d_{ii}>0

令

D1/2=[d11⋯0⋮⋱⋮0⋯dnn]D^{1/2} = \begin{bmatrix} \sqrt{d_{11}}&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&\sqrt{d_{nn}} \end{bmatrix}

L=L0D1/2L = L_0D^{1/2}

A=L0D1/2D1/2L0⊤=LL⊤A = L_0D^{1/2}D^{1/2}L_0^\top = LL^\top

Algorithom

image-20250305184229912

由于L=L0D12L = L_0D^{\frac12}，因此在每一步都把ll乘一个uiiu_{ii}就可以了。
由于L0=U0L_0 = U_0因此在第二个for循环中只需要计算一半。
之所以第二个循环j=i to nj=i~{\rm to}~n部分要用ajka_{jk}或ljkl_{jk}是因为A(1)=L1A=(I−l1e1⊤)A=A−l1a~1A^{(1)}=L_1A = (I-l_1e_1^\top)A = A-l_1\tilde a_1中L1L_1发生了变化，在Cholesky分解中，L大概长这样，所以之前的一些性质失效了。

Li:=(Ii−1000ai,i001ai,ibiIn−i),\begin{gathered}\mathbf{L}_i:=\begin{pmatrix}\mathbf{I}_{i-1}&0&0\\0&\sqrt{a_{i,i}}&0\\0&\frac{1}{\sqrt{a_{i,i}}}\mathbf{b}_i&\mathbf{I}_{n-i}\end{pmatrix},\end{gathered}

现在的算法基于以下事实：

A(i)=(Ii−1000ai,ibi∗0biB(i)),A(i+1)=(Ii−10001000B(i)−1ai,ibibi∗),A(i)=LiA(i+1)Li∗\begin{gather} \mathbf{A}^{(i)}=\begin{pmatrix}\mathbf{I}_{i-1}&0&0\\0&a_{i,i}&\mathbf{b}_i^*\\0&\mathbf{b}_i&\mathbf{B}^{(i)}\end{pmatrix},\\ \mathbf{A}^{(i+1)}=\begin{pmatrix}\mathbf{I}_{i-1}&0&0\\0&1&0\\0&0&\mathbf{B}^{(i)}-\frac{1}{a_{i,i}}\mathbf{b}_i\mathbf{b}_i^*\end{pmatrix},\\ \mathbf{A}^{(i)}=\mathbf{L}_i\mathbf{A}^{(i+1)}\mathbf{L}_i^* \end{gather}

Example

image-20250305211604957

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def Cholesky_Decomposition(A):
    assert A.shape[0] == A.shape[1]
    for k in range(A.shape[0] - 1):
        l = np.zeros((A.shape[1], 1))
        A[k,k] = np.sqrt(A[k,k])
        l[k + 1:, 0] = A[k + 1:, k] / A[k, k]
        # A_{ij} = a_{ij} - l_{ik}a_{kj}
        A[k + 1:, k + 1:] = A[k + 1:, k + 1:] - l[k + 1:, :] @ l[k+1:,:].T #A[k:k + 1, k + 1:]
        # use the lower triangular part of A to save l
        A[k + 1:, k:k + 1] = l[k + 1:, :]
    A[-1,-1] = np.sqrt(A[-1,-1])
    L = np.tril(A, k=0)
    U = np.triu(A, k=0)
    return L, U


A = np.array([[5, 2, 5], [2, 4, 3], [5, 3, 10]],dtype=np.float64)
L, U = Cholesky_Decomposition(A)