python实现LU分解与LUP分解

计算机求解线性方程组Ax=bAx=bAx=b过程中，更多的是采用数值计算方法求解而取代数学意义上效率更高的求逆运算，其中一个重要的问题是数值的稳定性。

上述线性方程组中AAA为nnn阶方阵，其中实际求解问题中只针对非奇异矩阵的情况下，这里首先介绍一种较为常见的LUPLUPLUP分解方式求解方法。

LUPLUPLUP方法求解

原理为找出满足条件的三个nnn阶方阵LUPLUPLUP使得
PA=LUPA=LU PA=LU
其中LLL为下三角矩阵，UUU为上三角矩阵，PPP为置换矩阵，在原方程Ax=bAx=bAx=b中会得到
PAx=PbLUx=Pb\begin{aligned} PAx&=Pb\\ LUx&=Pb\\ \end{aligned} PAxLUx=Pb=Pb
其中定义Ux=yUx=yUx=y得到Ly=PbLy=PbLy=Pb这时该位置向量yyy会被更容易的求得，之后将Ux=yUx=yUx=y以类似方式求解，其中做回推：
A=P−1LUAx=P−1LUxAx=P−1LyAx=P−1PbAx=b\begin{aligned} A&=P^{-1}LU\\ Ax&=P^{-1}LUx\\ Ax&=P^{-1}Ly\\ Ax&=P^{-1}Pb\\ Ax&=b\\ \end{aligned} AAxAxAxAx=P−1LU=P−1LUx=P−1Ly=P−1Pb=b
上述两步Ly=PbLy=PbLy=Pb和Ux=yUx=yUx=y为三角矩阵求解未知数，用前者表示说明：
y1=bπ[1],l21y1+y2=bπ[2],l31y1+l32y2+y3=bπ[3],ln1y1+ln2y2+ln3y3+⋯+yn=bπ[n].\begin{aligned} &y_{1} \quad&=b_{\pi[1]},\\ &l_{21} y_{1}+y_{2}&=b_{\pi[2]},\\ &l_{31} y_{1}+l_{32} y_{2}+y_{3}&=b_{\pi[3]},\\ &l_{n 1} y_{1}+l_{n 2} y_{2}+l_{n 3} y_{3}+\cdots+y_{n}&=b_{\pi[n]} . \end{aligned} y1l21y1+y2l31y1+l32y2+y3ln1y1+ln2y2+ln3y3+⋯+yn=bπ[1],=bπ[2],=bπ[3],=bπ[n].
由于PPP为置换矩阵所以PbPbPb只是被替换顺序的原先列向量，而LLL是对角线元素为111的下三角矩阵，那么求解过程可以很容易的表示为：
yi=bπ[i]−∑j=1i−1lijyjy_i=b_{\pi[i]}-\sum^{i-1}_{j=1}l_{ij}y_j yi=bπ[i]−j=1∑i−1lijyj
可以看出方式为朴素的从上到下迭代求解，时间复杂度为O(n2)O(n^2)O(n2)，第二步的上三角矩阵反向推得可知：
u11x1+u12x2+⋯+u1,n−2xn−2+u1,n−1xn−1+u1nxn=y1u22x2+⋯+u2,n−2xn−2+u2,n−1xn−1+u2nxn=y2⋮un−2,n−2xn−2+un−2,n−1xn−1+un−2,nxn=yn−2un−1,n−1xn−1+un−1,nxn=yn−1un,nxn=yn\begin{aligned} u_{11} x_{1}+u_{12} x_{2}+\cdots+u_{1, n-2} x_{n-2}+u_{1, n-1} x_{n-1}+u_{1 n} x_{n} &=y_{1} \\ u_{22} x_{2}+\cdots+u_{2, n-2} x_{n-2}+u_{2, n-1} x_{n-1}+u_{2 n} x_{n} &=y_{2} \\ & \vdots \\ u_{n-2, n-2} x_{n-2}+u_{n-2, n-1} x_{n-1}+u_{n-2, n} x_{n} &=y_{n-2} \\ u_{n-1, n-1} x_{n-1}+u_{n-1, n} x_{n} &=y_{n-1} \\ u_{n, n} x_{n} &=y_{n} \end{aligned} u11x1+u12x2+⋯+u1,n−2xn−2+u1,n−1xn−1+u1nxnu22x2+⋯+u2,n−2xn−2+u2,n−1xn−1+u2nxnun−2,n−2xn−2+un−2,n−1xn−1+un−2,nxnun−1,n−1xn−1+un−1,nxnun,nxn=y1=y2⋮=yn−2=yn−1=yn
这里UUU为对角线不为1的上三角矩阵，那么求解过程不仅仅要自下而上并且需要有除法进行，其中计算机内进行过程中会涉及到数值精度问题。
xi=(yi−∑j=i+1nuijxj)/uiix_i=(y_i-\sum_{j=i+1}^nu_{ij}x_j)/u_{ii} xi=(yi−j=i+1∑nuijxj)/uii
这里的置换矩阵PPP提高了解决问题的数值稳定性，主要解决的是分解矩阵时存在uiiu_{ii}uii趋近于零的情况，本文暂时讨论LUPLUPLUP分解计算方法，这里给出PythonPythonPython上述对应的的求解方法：

import numpy as np
def LUPSolve(L, U, pi, b):n = L.shape[0]y = np.zeros((n, 1))x = np.zeros((n, 1))for i in range(n):y[i] = b[pi[i]] - np.dot(L[i], y)for i in range(n-1, -1, -1):x[i] = (y[i]-np.dot(U[i], x))/U[i, i]return x

其中各项数据及其输出为：

L = np.array([[1, 0, 0], [0.2, 1, 0], [0.6, 0.5, 1]])
U = np.array([[5, 6, 3], [0, 0.8, -0.6], [0, 0, 2.5]])
b = np.array([3, 7, 8])
pi = np.array([2, 0, 1])
print(LUPSolve(L, U, pi, b))
->[[-1.4] [ 2.2] [ 0.6]]

前三个变量以知，pipipi的角色扮演为索引数组取代了运算中置换矩阵的作用

计算LULULU分解

首先考虑基础的不含PPP矩阵的LULULU分解，那么可以由高斯消元法推得：
A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)=(a11wTvA′),\begin{aligned} A &=\left(\begin{array}{c|ccc}a_{11} & a_{12} & \cdots & a_{1 n} \\ \hline a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n}\end{array}\right) \\ &=\left(\begin{array}{cc}a_{11} & w^{\mathrm{T}} \\ v & A^{\prime}\end{array}\right), \end{aligned} A=⎝⎜⎜⎜⎛a11a21⋮an1a12a22⋮an2⋯⋯⋱⋯a1na2n⋮ann⎠⎟⎟⎟⎞=(a11vwTA′),
此处在方阵中采用的递推式消元，上述公式指明了消元的一个状态，而每个状态都能成功分解一层，由矩阵乘法验证可知：
A=(10v/a11In−1)(a11wT0A′−vwT/a11)=(10v/a11In−1)(a11wT0L′U′)=(10v/a11L′)(a11wT0U′)=LU\begin{aligned} A &=\left(\begin{array}{cc}1 & 0 \\ v / a_{11} & I_{n-1}\end{array}\right)\left(\begin{array}{cc}a_{11} & w^{\mathrm{T}} \\ 0 & A^{\prime}-v w^{\mathrm{T}} / a_{11}\end{array}\right) \\ &=\left(\begin{array}{cc}1 & 0 \\ v / a_{11} & I_{n-1}\end{array}\right)\left(\begin{array}{cc}a_{11} & w^{\mathrm{T}} \\ 0 & L^{\prime} U^{\prime}\end{array}\right) \\ &=\left(\begin{array}{cc}1 & 0 \\ v / a_{11} & L^{\prime}\end{array}\right)\left(\begin{array}{cc}a_{11} & w^{\mathrm{T}} \\ 0 & U^{\prime}\end{array}\right) \\ &=L U \end{aligned} A=(1v/a110In−1)(a110wTA′−vwT/a11)=(1v/a110In−1)(a110wTL′U′)=(1v/a110L′)(a110wTU′)=LU
其中上述式子推导类似高等代数中分块矩阵乘法所应用的情景，最终通过迭代nnn次求解LULULU矩阵，可以看出并不用每一层处理后都进行分解，因为在递推的的状态下只需要将每层的数据保留继续计算下一层的AAA直至结束即可，最终上三角部分即为矩阵UUU，下三角部分使用111填补对角线即为矩阵LLL。

def LUDecomposition(A):n = A.shape[0]L, U = np.eye(n), np.zeros((n, n))for k in range(n):U[k, k] = A[k, k]for i in range(k+1, n):L[i, k] = A[i, k]/U[k, k]U[k, i] = A[k, i]for i in range(k + 1, n):for j in range(k+1, n):A[i, j] -= L[i, k]*U[k, j]return L, U

由于LLL对角线为111所以提前设置为单位矩阵，之后每行遍历对应第kkk行：

计算AkkA_{kk}Akk到AnnA_{nn}Ann的子矩阵；
保留AkkA_{kk}Akk到矩阵UUU；
维护子矩阵第kkk列数据与第kkk行数据；
计算下一阶段子矩阵:
A′=A′−vwTAkk∵Lik=AikUkk∴Aij′=Aij′−LikUkjA'=A'-\frac{vw^T}{A_{kk}}\\ \because L_{ik}=\frac{A_{ik}}{U_{kk}}\\ \therefore A_{ij}'=A_{ij}'-L_{ik}U_{kj}\\ A′=A′−AkkvwT∵Lik=UkkAik∴Aij′=Aij′−LikUkj

计算LUPLUPLUP分解

为了解决上述uiiu_{ii}uii有可能为零的情况，LUPLUPLUP分解再分解过程中加入了选择步骤：

选择当前状态当前列的绝对值最大的行索引交换到当前行，并将过程中的交换内容保存在置换矩阵中，这里我们在数据实现部分将交换保存在了数组中，主要作用是减少占用空间且作为索引数组使用，记作π[index]\pi[index]π[index]。

其中推导过程下节结合舒尔补做进一步说明，这里给出实现程序：

def LUPDecomposition(A):n = A.shape[0]pi = np.zeros((n, 1))for i in range(n):pi[i] = ifor k in range(n):p = 0for i in range(k, n):if np.abs(A[i, k]) > p:p = np.abs(A[i, k])k_ = iif p == 0:error("singular matrix")pi[k], pi[k_] = pi[k_], pi[k]A[[k, k_], :] = A[[k_, k], :]for i in range(k+1, n):A[i, k] /= A[k, k]for j in range(k+1, n):A[i, j] -= A[i, k]*A[k, j]return A

全部程序点击下方阅读原文获取