LUP分解求线性方程组

L、U、P是三个矩阵，满足PA = LU。其中L矩阵是一个单位下三角矩阵，U是一个上三角矩阵，P是一个置换矩阵。

我们要求解的是A x = b

$\begin{aligned} & Ax = b\\ & PAx = Pb\\ & LUx = Pb(根据上面的式子替换） \\ & 令Ux = y 得 \\ & Ly = Pb \\ \end{aligned}$

正向替换

正向替换的目标是通过Ly = Pb来求y

我们可以使用一个$\pi [1…n]$来替换P。对i = 1, 2, …, n, 元素$\pi [i]$表示$P_{i,\pi [i]} = 1$(第i行第$\pi [i]$列）。由于L是单位下三角矩阵，所以可以将式子写成

$\begin{matrix} y_1 & & & &= b_{\pi [1]}\\ y_2 l_{21} &+ y_2 & & &= b_{\pi [2]}\\ y_3 l_{31} &+ y_3 l_{32} &+ y_3 & &= b_{\pi [3]}\\ y_4 l_{41} &+ y_4 l_{42} &+ y_4 l_{43} &+ y_4 &= b_{\pi [4]}\\ ... \end{matrix}$

可以求出 $y1 = b{\pi [1]}$，之后将$y_1$代入第二个式子可以求出$y_2$，然后依次代入可以求出所有的y

然后可以求出$y1, y_2…$ 得 $y_i = b{\pi [i]} - \sum{j=1}^{i-1} l{ij}y_j$

反向替换

反向替换目标是通过 Ux = y求x

与之对应，过程不再细说。结果是$xi = ( y_i - \sum{j=i+1}^{n} u{ij}x_j ) / u{ii}$

LU分解

我们计算出 A = LU，则称这两个矩阵是A的一个LU分解

$\begin{aligned} \begin{bmatrix} a_{11}& a_{12} & ... & a_{1n}\\ a_{21}& a_{22} & ... & a_{2n}\\ \vdots& \vdots& \ddots & \vdots \\ a_{n1}& a_{n2} & ... & a_{nn} \end{bmatrix} &= \begin{bmatrix} a_{11} & \omega^T \\ v & A^{'} \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0\\ v/a_{11} & E_{n-1} \end{bmatrix} \begin{bmatrix} a_{11} & \omega^T\\ 0 & A^{'} - v \omega^T / a_{11} \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0\\ v/a_{11} & L^{'} \end{bmatrix} \begin{bmatrix} a_{11} & \omega^T\\ 0 & U^{'} \end{bmatrix}\\ \end{aligned}$

其中v是n-1的列向量， $\omega^T$是n-1的行向量.

最后一步是因为$A^{‘} - v \omega^T / a_{11}$也可以进行分解，并且分解的结果是$L^{‘}U^{‘}$

for(int i=0; i<n; i++)
{
    u[i][i] = a[i][i];
    for(int k=i+1; k<n; k++)
    {
        l[i][k] = a[i][k] / u[i][i];
        u[k][i] = a[k][i];
        for(int k=i+1; k<n; k++)
        {
            for(int j=k+1; j<n; j++)
            {
                a[i][j] = a[i][j] - l[i][k] * u[k][j];
            }
        }
    }
}

他有一个限制条件u[i][i]不能是0，否则会出现除0错误

LUP分解

为了保证a[i][i] != 0, 我们可以使用该列中最大元素a[k][1]替换a[1][1],为了使方程组还成立，我们可以将第1行和第k行互换，等价于用一个置换矩阵Q乘在A的左边

$QA = \begin{bmatrix} 1 & 0\\ v/a_{k1} & E_{n-1} \end{bmatrix} \begin{bmatrix} a_{k1} & \omega^T\\ 0 & A^{'} - v \omega ^{T} / a_{k1} \end{bmatrix}\\$

和上面类似，就是$a{11}$变成了$a{k1}$

矩阵求逆

矩阵求逆也可以使用LUP分解我们可以令 AX = $E_n$ 其中 $E_n$是n阶单位矩阵.X是n个列向量的和

我们可以把它看成n个$AX_i = E_i$来进行求解，然后最后将$X_i$合并即可