【课堂笔记】LU分解，Cholesky分解
文章目录
- 回顾
- L U LU LU分解
- 高斯消元法(Gaussian Elimination)
- 转动LU分解(pivoting LU)
- 稳定性分析
- Cholesky 分解
回顾
我们已经有了三种矩阵分解的方法：
- Q R QR QR分解， A = Q R A=QR A=QR，其中 Q Q Q是正交矩阵， R R R是上三角矩阵
- 特征值分解， A = V D V − 1 A = VDV^{-1} A=VDV−1，其中 A A A是可对角化的， V V V是奇异的， D D D是对角的
- S V D SVD SVD奇异值分解， A = U ∑ V ∗ A=U\sum V^* A=U∑V∗，其中 U , V U, V U,V是正交的， ∑ \sum ∑是对角的
而这篇将介绍另外两种分解方法。
L U LU LU分解
令 A ∈ C m × m A \in \mathbb{C}^{m \times m} A∈Cm×m，最简单的形式为 A = L U A = LU A=LU，其中 L L L是单位下三角矩阵， U U U是上三角矩阵。更一般的可以写成 P A = L U PA = LU PA=LU，其中 P P P 是置换矩阵(permutation)
高斯消元法(Gaussian Elimination)
在解线性方程组 A x = b Ax = b Ax=b时，我们会先利用高斯消元将 A A A转化为上三角矩阵。
- 例如 { 2 x + 4 y − 2 z = 2 4 x + 9 y − 3 z = 8 − 2 x − 3 y + 7 z = 10 \left\{\begin{matrix} 2x+4y-2z=2 \\ 4x+9y-3z=8 \\ -2x-3y+7z=10 \end{matrix}\right. ⎩ ⎨ ⎧2x+4y−2z=24x+9y−3z=8−2x−3y+7z=10
[ 2 4 − 2 2 4 9 − 3 8 − 2 − 3 7 10 ] → L 1 [ 2 4 − 2 2 0 1 1 4 0 1 5 12 ] → L 2 [ 2 4 − 2 2 0 1 1 4 0 0 4 8 ] \left[ \begin{array}{ccc|c} 2 & 4 & -2 & 2 \\ 4 & 9 & -3 & 8 \\ -2 & -3 & 7 & 10 \\ \end{array} \right] \xrightarrow{L_1} \left[ \begin{array}{ccc|c} 2 & 4 & -2 & 2 \\ 0 & 1 & 1 & 4 \\ 0 & 1 & 5 & 12 \\ \end{array} \right] \xrightarrow{L_2} \left[ \begin{array}{ccc|c} 2 & 4 & -2 & 2 \\ 0 & 1 & 1 & 4 \\ 0 & 0 & 4 & 8 \\ \end{array} \right] 24−249−3−2−372810 L1 200411−2152412 L2 200410−214248
我们记 L k = [ 1 1 1 − l k + 1 , k ⋮ ⋱ − l m , k 1 ] L_k = \begin{bmatrix} 1 & & & & & \\ & 1 & & & & \\ & & 1 & & & \\ & & -l_{k+1, k} & & & \\ & & \vdots & \ddots & & \\ & & -l_{m, k} & & & 1 \end{bmatrix} Lk= 111−lk+1,k⋮−lm,k⋱1 ，其中 l j k = x j k / x k k l_{jk} = x_{jk}/x_{kk} ljk=xjk/xkk，我们有 L k x k = [ x 1 k ⋮ x k k 0 ⋮ 0 ] L_kx_k = \begin{bmatrix} x_{1k} \\ \vdots \\ x_{kk} \\ 0 \\ \vdots \\ 0 \end{bmatrix} Lkxk= x1k⋮xkk0⋮0
令 l k ⃗ \vec{l_k} lk为 L k L_k Lk的第 k k k列，则 L k = I − l k ⋅ e k ∗ L_k = I - l_k \cdot e_k^* Lk=I−lk⋅ek∗
若 j > k j > k j>k，则 L k L j = ( I − l k e k ∗ ) ( I − l j e j ∗ ) = I − l k e k ∗ − l j e j ∗ L_kL_j = (I - l_ke_k^*)(I-l_je^*_j)=I-l_ke_k^*-l_je_j^* LkLj=(I−lkek∗)(I−ljej∗)=I−lkek∗−ljej∗
我们由上面的变化可知， L m − 1 ⋯ L 2 L 1 A = U L_{m-1}\cdots L_2L_1A=U Lm−1⋯L2L1A=U， A = ( L 1 − 1 L 2 − 1 ⋯ L m − 1 − 1 ) ⋅ U A = (L_1^{-1}L_2^{-1}\cdots L_{m-1}^{-1})\cdot U A=(L1−1L2−1⋯Lm−1−1)⋅U。
令 L = L 1 − 1 L 2 − 1 ⋯ L m − 1 − 1 L = L_1^{-1}L_2^{-1}\cdots L_{m-1}^{-1} L=L1−1L2−1⋯Lm−1−1，则 A = L U A = LU A=LU。可验证 L = [ 1 l 21 1 l 31 l 32 1 ⋮ ⋮ ⋱ l m 1 l m 2 ⋯ ⋯ 1 ] L = \begin{bmatrix} 1 & & & & \\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & & \ddots & \\ l_{m1} & l_{m2} & \cdots & \cdots & 1 \end{bmatrix} L= 1l21l31⋮lm11l32⋮lm21⋯⋱⋯1
下面是实现算法的伪代码，精度为 2 3 m 3 \frac{2}{3}m^3 32m3 flops
'''
Algorithm for A=LU
Input A
'''
U, L = A, I
for k=1: m-1for j=k+1 : mL[j, k] = U[j, k] / U[k, k]U[j, k:m] = U[j, k:m] - L[j,k]U[k, k:m]
转动LU分解(pivoting LU)
上述伪代码中有除以 U [ k , k ] U[k, k] U[k,k]的操作，当它很小或者为 0 0 0时就会出一些问题。
- 对矩阵 A A A有性质上的要求，例如 A = [ 0 1 1 0 ] A = \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} A=[0110]就没有 L U LU LU分解。
- 会产生数值不稳定问题，例如 A = [ ϵ 1 1 0 ] = [ 1 0 1 / ϵ 1 ] [ ϵ 1 0 − 1 / ϵ ] A = \begin{bmatrix} \epsilon & 1\\ 1 & 0 \end{bmatrix}=\begin{bmatrix} 1 & 0\\ 1/\epsilon & 1 \end{bmatrix}\begin{bmatrix} \epsilon & 1\\ 0 & -1/\epsilon \end{bmatrix} A=[ϵ110]=[11/ϵ01][ϵ01−1/ϵ]，数值溢出达到 O ( ϵ − 2 ) O(\epsilon^{-2}) O(ϵ−2)
我们引入置换矩阵 P P P，通过换主元的方式让所有非奇异的矩阵都能写为 P A = L U PA = LU PA=LU。
在第 k k k步中，我们找出第 k k k列从第 k k k行的最大元素，将它与第 k k k行交换，然后再进行消元。
将交换步骤写成矩阵 P i P_i Pi，则 L m − 1 P m − 1 . . . L 2 P 2 L 1 P 1 A = U L_{m-1}P_{m-1}...L_2P_2L_1P_1A = U Lm−1Pm−1...L2P2L1P1A=U。
记 L m − 1 ′ = L m − 1 , L j ′ = P m − 1 P m − 2 . . . P j + 1 L j P j + 1 − 1 . . . P m − 2 − 1 P m − 1 − 1 L_{m-1}' = L_{m-1}, L_{j}' = P_{m-1}P_{m-2}...P_{j+1}L_jP_{j+1}^{-1}...P_{m-2}^{-1}P_{m-1}^{-1} Lm−1′=Lm−1,Lj′=Pm−1Pm−2...Pj+1LjPj+1−1...Pm−2−1Pm−1−1，就可以写成：
L m − 1 ′ L m − 2 ′ . . . L 1 ′ ⋅ ( P m − 1 P m − 2 . . . P 1 ) A = U L_{m-1}'L_{m-2}'...L_1' \cdot (P_{m-1}P_{m-2}...P_1)A = U Lm−1′Lm−2′...L1′⋅(Pm−1Pm−2...P1)A=U
稳定性分析
设 A ∈ R m × m A \in \mathbb{R}^{m \times m} A∈Rm×m，计算的 L ~ , U ~ \tilde{L}, \tilde{U} L~,U~满足 L ~ U ~ = A + E \tilde{L}\tilde{U} = A + E L~U~=A+E，则有：
∣ E ∣ ≤ 2 ( m − 1 ) ϵ math ( ∣ A ∣ + ∣ L ∣ ∣ U ∣ ) |E| \le 2(m-1)\epsilon_{\text{math}}(|A|+|L||U|) ∣E∣≤2(m−1)ϵmath(∣A∣+∣L∣∣U∣)
设 P A = L U ∈ R m × m PA = LU \in \mathbb{R}^{m \times m} PA=LU∈Rm×m，计算的 P ~ , L ~ , U ~ \tilde{P}, \tilde{L}, \tilde{U} P~,L~,U~满足
L ~ U ~ = P ~ A + δ A ∥ δ A ∥ ∥ A ∥ = O ( ρ ϵ math ) ρ = max i j ∣ u i j ∣ max i j ∣ a i j ∣ \tilde{L}\tilde{U} = \tilde{P}A + \delta A \\ \frac{\|\delta A\|}{\|A\|} = O(\rho \epsilon_{\text{math}}) \\ \rho = \frac{\max_{ij}|u_{ij}|}{\max_{ij}|a_{ij}|} L~U~=P~A+δA∥A∥∥δA∥=O(ρϵmath)ρ=maxij∣aij∣maxij∣uij∣
其中 ρ \rho ρ是增长因子(growth factor)
Cholesky 分解
若 A ∈ C m × m A \in \mathbb{C}^{m\times m} A∈Cm×m是Hermitian的，即 a i j = a j i ‾ a_{ij} = \overline{a_{ji}} aij=aji，且 A A A是正定的，则 A A A可以分解为 A = R ∗ R A=R^*R A=R∗R
A = [ a 11 w ∗ w k ] = [ α w / α I ] [ 1 k − w w ∗ / a 11 ] [ α w ∗ / α I ] A = \left[ \begin{array}{c|c} a_{11} & w^* \\ \hline w & k \end{array} \right] = \begin{bmatrix} \alpha & \\ w/\alpha & I \end{bmatrix} \begin{bmatrix} 1 & \\ & k-ww^*/a_{11} \end{bmatrix} \begin{bmatrix} \alpha & w^*/\alpha\\ & I \end{bmatrix} A=[a11ww∗k]=[αw/αI][1k−ww∗/a11][αw∗/αI]
其中 α = a 11 \alpha = \sqrt{a_{11}} α=a11， k − w w ∗ / a 11 k - ww^*/a_{11} k−ww∗/a11是半正定的，
x ∗ [ 1 k − w w ∗ / a 11 ] x = x ∗ [ α w / α I ] − 1 ⏟ y ∗ A [ α w ∗ / α I ] − 1 x ⏟ y x^*\begin{bmatrix} 1 & \\ & k-ww^*/a_{11} \end{bmatrix}x = \underbrace{ x^* \begin{bmatrix} \alpha & \\ w/\alpha & I \end{bmatrix}^{-1} }_{y^*} A \underbrace{ \begin{bmatrix} \alpha & w^*/\alpha \\ & I \end{bmatrix}^{-1} x }_{y} x∗[1k−ww∗/a11]x=y∗ x∗[αw/αI]−1Ay [αw∗/αI]−1x
重复上述过程直到中间变成 I I I
'''
Algorithm for Cholesky factorization
Input A
'''
R=A
for k=1:mfor j=k+1:mR[j, j:m] = R[j, j:m] - R[k, j:m] * conj(R[k, j]) / R[k, k]R[k, k:m] = R[k, k:m] / sqrt(R[k, k])
若计算结果为 R ~ \tilde{R} R~，有
R ~ ∗ R ~ = A + δ A ∥ δ A ∥ ∥ A ∥ = O ( ϵ math ) \tilde{R}^*\tilde{R} = A + \delta A \\ \frac{\|\delta A\|}{\|A\|} = O(\epsilon_{\text{math}}) R~∗R~=A+δA∥A∥∥δA∥=O(ϵmath)