LUP分解求线性方程组

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

我们要求解的是A x = b

\begin{aligned} A x = b \\ P A x = P b \\ L U x = P b (根 据 上 面 的 式 子 替 换 ） \\ 令 U x = y 得 \\ L y = P b \end{aligned}

正向替换

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

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

\begin{matrix} y_{1} & = b_{π [1]} \\ y_{2} l_{21} & + y_{2} & = b_{π [2]} \\ y_{3} l_{31} & + y_{3} l_{32} & + y_{3} & = b_{π [3]} \\ y_{4} l_{41} & + y_{4} l_{42} & + y_{4} l_{43} & + y_{4} & = b_{π [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{array}{c} a_{11} & a_{12} & . . . & a_{1 n} \\ a_{21} & a_{22} & . . . & a_{2 n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n 1} & a_{n 2} & . . . & a_{n n} \end{array}] & = [\begin{array}{c} a_{11} & ω^{T} \\ v & A^{^{'}} \end{array}] \\ = [\begin{array}{c} 1 & 0 \\ v / a_{11} & E_{n - 1} \end{array}] [\begin{array}{c} a_{11} & ω^{T} \\ 0 & A^{^{'}} - v ω^{T} / a_{11} \end{array}] \\ = [\begin{array}{c} 1 & 0 \\ v / a_{11} & L^{^{'}} \end{array}] [\begin{array}{c} a_{11} & ω^{T} \\ 0 & U^{^{'}} \end{array}] \end{aligned}

其中v是n-1的列向量， $ω^{T}$ 是n-1的行向量.

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

inform7

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的左边

Q A = [\begin{matrix} 1 & 0 \\ v / a_{k 1} & E_{n - 1} \end{matrix}] [\begin{matrix} a_{k 1} & ω^{T} \\ 0 & A^{^{'}} - v ω^{T} / a_{k 1} \end{matrix}]

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

矩阵求逆

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

我们可以把它看成n个 $A X_{i} = E_{i}$ 来进行求解，然后最后将 $X_{i}$ 合并即可