矩阵 LUP 分解解线性方程组求行列式值矩阵求逆算法说解

算法：矩阵 LUP 分解

本文着笔于矩阵 LUP 分解算法，以及利用矩阵的 LUP 分解来解线性方程组、求矩阵对应行列式的值、求逆矩阵。

对于矩阵的定义代码如下：

struct Matrix
{double dat[MAX_N][MAX_N],det,LUdat[MAX_N][MAX_N],rev[MAX_N][MAX_N];int dgr,pi[MAX_N];bool LUP_Decomposition();VectorColumn LUP_Solve(VectorColumn);void GetDeterminant();void GetReverse();
};

矩阵的 LUP 分解

矩阵 A 的 LUP 分解即为找到下三角矩阵 L、上三角矩阵 U、行置换矩阵 P 使得 PA = LU ，实现方法为高斯消元。

E=I=(10⋯001⋯0⋮⋮⋱⋮00⋯1),eij={1i=j0i≠j\begin{matrix} E \end{matrix} = \begin{matrix} I \end{matrix} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix}, e_{ij} = \begin{cases} 1 & i = j \\0 & i \neq j \end{cases} E=I=⎝⎜⎜⎜⎛10⋮001⋮0⋯⋯⋱⋯00⋮1⎠⎟⎟⎟⎞,eij={10i=ji̸=j

初始矩阵 P = E。对于矩阵 A，我们可以将它分解成这样两个矩阵的乘积（这里直接借用了 Introduction to Algorithms 里的描述）
A=(a11ωTνA′)=(10ν/a11In−1)(a11ωT0A′−νωT/a11)\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & \omega^T \\ \nu & A' \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ \nu / a_{11} & I_{n-1} \end{pmatrix} \begin{pmatrix} a_{11} & \omega^T \\ 0 & A' - \nu \omega^T / a_{11} \end{pmatrix} A=(a11νωTA′)=(1ν/a110In−1)(a110ωTA′−νωT/a11)

这样，将这个规模为 n 的问题分解为一个规模为 n - 1 的子问题，迭代之即可求得矩阵 PA 的 LU 分解。最后得到的矩阵 U 相当于高斯消元后的结果，因为求矩阵 U 的实质是用当前规模矩阵的第一行的相应倍数减其它行，将其它行首列元素消为 0。

为了避免除数为零或过小引起精度问题，在迭代到矩阵第 i 行时，我们可以找到 j 使得 ∣Aj,i∣\vert A_{j, i} \vert∣Aj,i∣ 尽可能大，其中 i≤j≤ni \leq j \leq ni≤j≤n。当 i 不等于 j 时，交换矩阵 PA 的第 i 行与第 j 行，相当于矩阵 A 不变，矩阵 P 的第 i 行与第 j 行交换。

对于矩阵的 LUP 分解代码如下：

#define ABS(x) \
({ \typeof(x) __tmp_x=(x); \__tmp_x<0?-__tmp_x:__tmp_x; \
})#define EXCHANGE(a,b) \
{ \typeof(a) __tmp_a=(a); \a=b,b=__tmp_a; \
}bool Matrix::LUP_Decomposition()
{int p;for(int i=1;i<=dgr;i++){pi[i]=i;for(int j=1;j<=dgr;j++)LUdat[i][j]=dat[i][j];}det=1;for(int i=1;i<=dgr;i++){p=i;for(int j=i+1;j<=dgr;j++)if(ABS(LUdat[j][i])>ABS(LUdat[p][i]))p=j;if(LUdat[p][i]==0)return false;if(p!=i){det=-det;EXCHANGE(pi[i],pi[p]);}for(int j=1;j<=dgr;j++)EXCHANGE(LUdat[i][j],LUdat[p][j]);for(int j=i+1;j<=dgr;j++){LUdat[j][i]/=LUdat[i][i];for(int k=i+1;k<=dgr;k++)LUdat[j][k]-=LUdat[j][i]*LUdat[i][k];}}return true;
}

时间复杂度 O(n3)O(n^3)O(n3)。

矩阵求解线性方程组

对于线性方程组
{a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋮⋮⋮⋮an1x1+an2x2+⋯+annxn=bn\begin{cases} a_{11} x_1 + & a_{12} x_2 + \cdots + & a_{1n} x_n = & b_1 \\ a_{21} x_1 + & a_{22} x_2 + \cdots + & a_{2n} x_n = & b_2 \\ \vdots & \vdots & \vdots & \vdots \\ a_{n1} x_1 + & a_{n2} x_2 + \cdots + & a_{nn} x_n = & b_n \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧a11x1+a21x1+⋮an1x1+a12x2+⋯+a22x2+⋯+⋮an2x2+⋯+a1nxn=a2nxn=⋮annxn=b1b2⋮bn

构建矩阵 A 和列向量 B，使得
A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann),B=(b1b2⋮bn)\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix}, \begin{matrix} B \end{matrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} A=⎝⎜⎜⎜⎛a11a21⋮an1a12a22⋮an2⋯⋯⋱⋯a1na2n⋮ann⎠⎟⎟⎟⎞,B=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞

矩阵求解线性方程组即对于矩阵 A 和列向量 B，找到列向量 X，使得
AX=B\begin{matrix} A \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} B \end{matrix} AX=B

两边同乘 P，得
LUX=PAX=PB\begin{matrix} L \end{matrix} \begin{matrix} U \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} A \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} B \end{matrix} LUX=PAX=PB

设列向量 Y = UX，得
LY=PB\begin{matrix} L \end{matrix} \begin{matrix} Y \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} B \end{matrix} LY=PB

设数组 π1..n\pi_{1 .. n}π1..n，其中 πi\pi_iπi 表示矩阵 P 中第 i 行第 πi\pi_iπi 列元素为 1，第 i 行其余元素均为 0，则 LY = PB 对应的方程组为
{l11y1=bπ1l21y1+l22y2=bπ2⋮⋮⋮ln1y1+ln2y2+⋯+lnnyn=bπn\begin{cases} l_{11} y_1 & & & = & b_{\pi_1} \\ l_{21} y_1 + & l_{22} y_2 & & = & b_{\pi_2} \\ \vdots & \vdots & & & \vdots \\ l_{n1} y_1 + & l_{n2} y_2 + \cdots + & l_{nn} y_n & = & b_{\pi_n} \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧l11y1l21y1+⋮ln1y1+l22y2⋮ln2y2+⋯+lnnyn===bπ1bπ2⋮bπn

由定义，矩阵 L 中对角线元素均为 1。从上往下依次解出 y1..ny_{1 .. n}y1..n，这样解出 Y 后，利用方程 UX = Y，解出 X 即可。对应方程组为
{u11x1+u12x2+⋯+u1nxn=y1u22x2+⋯+u2nxn=y2⋮⋮unnxn=yn\begin{cases} u_{11} x_1 + & u_{12} x_2 + \cdots + & u_{1n} x_n & = & y_1 \\ & u_{22} x_2 + \cdots + & u_{2n} x_n & = & y_2 \\ & & \vdots & & \vdots \\ & & u_{nn} x_n & = & y_n \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧u11x1+u12x2+⋯+u22x2+⋯+u1nxnu2nxn⋮unnxn===y1y2⋮yn

从下往上依次解出 xn..1x_{n .. 1}xn..1。在具体实现过程中，由于结果中 Y 无需保留，列向量 X 和 Y 可共用存储空间。

对于求解线性方程组的代码如下：

struct VectorColumn
{double dat[MAX_N];int raw;
};VectorColumn Matrix::LUP_Solve(VectorColumn b)
{VectorColumn x;x.raw=dgr;for(int i=1;i<=dgr;i++){x.dat[i]=b.dat[pi[i]];for(int j=1;j<i;j++)x.dat[i]-=LUdat[i][j]*x.dat[j];}for(int i=dgr;i>=1;i--){for(int j=dgr;j>i;j--)x.dat[i]-=LUdat[i][j]*x.dat[j];x.dat[i]/=LUdat[i][i];}return x;
}

时间复杂度 O(n2)O(n^2)O(n2)。（若算上矩阵 LUP 分解，总时间复杂度 O(n3)O(n^3)O(n3)）

矩阵的行列式值求解

由行列式性质可知
∣P∣∣A∣=∣L∣∣U∣\begin{vmatrix} P \end{vmatrix} \begin{vmatrix} A \end{vmatrix} = \begin{vmatrix} L \end{vmatrix} \begin{vmatrix} U \end{vmatrix} ∣∣P∣∣∣∣A∣∣=∣∣L∣∣∣∣U∣∣

其中，L 为对角线元素均为 1 的下三角矩阵，故其行列式值为 1；U 为上三角矩阵，其行列式值为对角线元素之积；P 初始行列式值为 1，每对 P 进行一次行交换，相当于行列式值取其相反数。

在具体实现过程中，行列式值的初值为 1 或 -1 在 LUP 分解过程中决定，此部分代码见上文。

对于矩阵的行列式值求解代码如下：

void Matrix::GetDeterminant()
{for(int i=1;i<=dgr;i++)det*=LUdat[i][i];
}

时间复杂度 O(n)O(n)O(n)。（若算上矩阵 LUP 分解，总时间复杂度 O(n3)O(n^3)O(n3)）

矩阵求逆

求矩阵 AAA 的逆矩阵，即对于矩阵 AAA 和单位矩阵 EEE，找到矩阵 A−1A^{-1}A−1，使得 AA−1=EA A^{-1} = EAA−1=E。

对于矩阵 A−1A^{-1}A−1 和矩阵 EEE 的任一列，设列数为 iii，构建列向量
Acol(i)−1=(a1ia2i⋮ani),Ecol(i)=(e1ie2i⋮eni)\begin{matrix} A_{col(i)}^{-1} \end{matrix} = \begin{pmatrix} a_{1i} \\ a_{2i} \\ \vdots \\ a_{ni} \end{pmatrix}, \begin{matrix} E_{col(i)} \end{matrix} = \begin{pmatrix} e_{1i} \\ e_{2i} \\ \vdots \\ e_{ni} \end{pmatrix} Acol(i)−1=⎝⎜⎜⎜⎛a1ia2i⋮ani⎠⎟⎟⎟⎞,Ecol(i)=⎝⎜⎜⎜⎛e1ie2i⋮eni⎠⎟⎟⎟⎞

则有
AAcol(i)−1=Ecol(i)\begin{matrix} A \end{matrix} \begin{matrix} A_{col(i)}^{-1} \end{matrix} = \begin{matrix} E_{col(i)} \end{matrix}AAcol(i)−1=Ecol(i)

套用矩阵求解线性方程组的方法即可。

对于矩阵求逆的代码如下：

void Matrix::GetReverse()
{VectorColumn curcol,revcol;curcol.raw=dgr;for(int i=1;i<=dgr;i++)curcol.dat[i]=0;for(int i=1;i<=dgr;i++){curcol.dat[i-1]=0,curcol.dat[i]=1;revcol=LUP_Solve(curcol);for(int j=1;j<=dgr;j++)rev[j][i]=revcol.dat[j];}
}

时间复杂度 O(n3)O(n^3)O(n3)。