======Gaussian Elimination====== LAFF Week 6 Notes ---- 本章的主要内容: * 高斯消元、高斯变换的推导及其算法实现 * LU 分解的推导及其算法实现 * 相关算法复杂度的掌握 ====Linear Systems & Gaussian Elimination==== 在中学的时候我们都学习过如何解三元一次方程组。来看看下面三元一次方程组是怎么解的: $$ \left\{\begin{matrix} 2x_0& &+&4x_1 &-&2x_2 &=& -10&\\ 4x_0& &-&2x_1 &+&6x_2 &=& 20& \\ 6x_0& &-&4x_1 &+&2x_2 &=&18& \end{matrix}\right. $$ 解的方法很简单:**消元**;也就是用方程组内的几个线性关系通过**加减**、**与标量相乘** 和**等式的顺序变换**,可以达到减少未知数的效果。上面的例子,我们可以观察到: * 第二行加上 ''-2'' 乘以第一行的结果,可以消去第二行中的 $x_0$ * 第三行加上 ''-3'' 乘以第一行的结果,可以消去第三行中的 $x_0$ 完成上面两步以后,我们的方程组变成了这个样子: $$ \left\{\begin{matrix} 2x_0& &+&4x_1 &-&2x_2 &= &-10\\ & &-&10x_1 &+&10x_2 &= &40 \\ & &-&16x_1 &+&8x_2 &= &48 \end{matrix}\right. $$ 我们发现第二行和第三行的方程只存在 $x_1$、$x_2$ 两个未知数了。按照之前的方法,我们在将第三行加上第二行乘以 ''-1.6'' 的结果,发现可以消除第三行中的 $x_1$。得到的结果如下: $$ \left\{\begin{matrix} 2x_0& &+&4x_1 &-&2x_2 &= &-10\\ & &-&10x_1 &+&10x_2 &= &40 \\ & && &-&8x_2 &= &-16 \end{matrix}\right. $$ ===Appended Matrices=== 我们可以将上面的方程组写成如下不带未知变量的形式: \\ $$ \left( \begin{array}{r r r | r} 2 & 4 & -2 & -10 \\ 4 & -2 & 6 & 20 \\ 6 & -4 & 2 & -18 \end{array} \right) $$ \\ 这样的形式被称为 //Appended Matrices//,这种形式可以更方便的将方程组与矩阵联系起来。让我们把上面这个流程写成该矩阵的形式: \\ $$ \left( \begin{array}{r r r | r} 2 & 4 & -2 & -10 \\ & -2 & 6 & 40 \\ & -4 & 2 & -48 \end{array} \right)\tag{step.1} $$ $$ \left( \begin{array}{r r r | r} 2 & 4 & -2 & -10 \\ & -2 & 6 & 40 \\ & & 2 & -16 \end{array} \right)\tag{step.2} $$ \\ 如果我们把每一步消去的未知变量的位置都用 ''0'' 填充的话,我们会发现整个消元的过程实际上是**目标矩阵变换为一个三角矩阵的过程(此处是上三角矩阵): ** $$ \left( \begin{array}{r r r | r} 2 & 4 & -2 & -10 \\ 0 & -2 & 6 & 40 \\ 0& 0 & 2 & -16 \end{array} \right) $$ \\ ===Gauss Transforms=== 我们再来看看当初将例子中的原矩阵变换为三角矩阵的(以 $L$ 表示方程组的行): - $L_2 + (-2) * L_1$ - $L_3 + (-3) * L_1$ - $L_3 + (-1.6) * L_2$ 现在我们将上面例子中的 Appended Matrices 的左边部分表示为一个单位矩阵与其的乘积: \\ \\ $$ \left( \begin{array}{r r r} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \left( \begin{array}{r r r | r} 2 & 4 & -2 & -10 \\ 4 & -2 & 6 & 20 \\ 6 & -4 & 2 & -18 \end{array} \right) $$ \\ \\ 如果我们注意观察的话,我们第一步和第二步是需要消除第二行和第三行中的 $x_0$。先来看看第二行,这里发生的点积如下图: \\ \\
a21 = a21 / a11 //it means l21=a21/a11, then a21 = l21;
A22 = A22 - a21 * a12t // means A22=A22-l21*a12t , rank 1 update
下面是完整的算法:
\\
\\
{{ math:linear_algebra:laff:6_2_5_1_ge.png?400 |}}
\\
\\
===Forward Substitution 的算法设计===
Forward Substitution 就是将表示消元步骤的单位矩阵与目标矩阵的右边表示结果的向量相乘,得到与变换后的三角矩阵匹配的结果向量。利用之前的表示方法,我们也可以将该迭代过程表示成如下的矩阵乘法:
\\
\\
$$
\left( \begin{array}{r | r r}
\color{blue}{I} & 0 & 0 \\
\hline
0 & 1 & 0 \\
0 & \color{red}{-l_{21}} & \color{blue}{I}
\end{array} \right)
\left( \begin{array}
y_0 \\
\hline
\color{blue}{\beta_1} \\
\hline
b_2
\end{array} \right)
=
\left( \begin{array}
y_0 \\
\hline
\color{blue}{\chi_1} \\
\hline
b_2^{new}
\end{array} \right)
$$
\\
化简左边以后可以得到:
\\
\\
$$
\left( \begin{array}
y_0 \\
\hline
\color{blue}{\beta_1} \\
\hline
\color{#0F0}-l_{21}\beta{11}+b_2
\end{array} \right)
=
\left( \begin{array}
y_0 \\
\hline
\color{blue}{\chi_1} \\
\hline
b_2^{new}
\end{array} \right)
$$
\\
\\
可以看到的是,每次更新中,只有 $b_2$ 的值发生了更新,也就是每次更新只对结果变量的余下部分造成影响。因此每次迭代的内容实际上就是更新 $b_2$。又因为 $l_{21}$ 实际上就是更新后的 $a_{21}$,因此整个迭代过程可以记做:
b_2 = b_2 - a21b11
整个算法演示如下:
\\
\\
{{ math:linear_algebra:laff:6_2_5_1_fs.png?400 |}}
\\
\\
====LU factorization====
===什么是LU分解===
先来了解一下 LU 分解的概念:
* $A$ 是一个 ''n*n'' 的矩阵
* $L$ 是一个 ''n*n'' 的**下**三角**单位**矩阵
* $U$ 是一个 ''n*n'' 的**上**三角矩阵
对于以上条件,如果存在 $A = LU$,那么我们就称 $LU$ 是 $A$ 的 $LU$ 分解。
\\
\\
===LU 分解的算法===
\\
a21 = a21/a11
A22 = A22 - a21a12t
整个算法如下:
\\
\\
{{ math:linear_algebra:laff:6_3_1_1_lu.png?400 |}}
\\
\\
是不是跟高斯消元很像?$L$ 是不是很像我们用于表示消元的高斯变换,$U$ 像不像消元后的目标矩阵?如果还不明白的话,来从头看一遍为什么要进行 $LU$ 分解:
- 我们想要算出 $Ax = b$
- 我们用 $LU$ 替换 $A$,也就是 $A = LU$
- 那么就有 $(LU)x = b, L(Ux) = b$
- 令 $Ux = z$,因此有 $L(z) = b$
- 通过 $L(z) = b$ 算出 $z$
- 通过 $Ux = z$ 算出 $x$
- 因为是三角矩阵相关运算,因此比直接算容易很多
===LU分解的 forward substitution:Lz = b===
这一个分解是解出 $z$ 是多少,并提供给后面的 $Ux = z$。
\\
\\
回想一下之前高斯变换的 forward substitution,我们用一个代表消元的矩阵乘以
按照老办法,我们将 $L$ 也按四象分割。再自上往下分割 向量 $z$ 与 目标向量 $b$。根据 $Lz = b$ 我们可以得到以下等式:
\\
\\
$$
\left( \begin{array}{c | c}
1 & 0 \\
\hline
l_{21} & L_{22} \\
\end{array}
\right)
\left( \begin{array} {r}
\zeta_1 \\
\hline
z_2
\end{array}
\right)
=
\left(\begin{array}{r} \beta_1 \\ \hline b_2 \end{array}\right)
$$
\\
\\
化简后:
\\
\\
$$
\left( \begin{array} {c}
\zeta_1 \\
\hline
l_{21}\zeta_1 + L_{22}z{2}
\end{array}
\right)
=
\left(\begin{array}{c} \beta_1 \\ \hline b_2 \end{array}\right)
$$
\\
\\
因此可以得到关系:
* $\zeta_1 = \beta_1$
* $L_{22} z_2 = b_2 - \zeta_1 l_{21}$
注意 $L_{22}z_{2}$,是一个下三角单元矩阵与向量的乘积。这里与 $LU$ 分解类似,该乘积代表了更新后 $b$ 的值。而对其更新的是 $-l_{21}\zeta_1$。因此每次迭代中的更新的内容如下:
b2 = b2 -beta1 * l21
迭代开头的 ''b2'' 就是更新过的 $b$,也就是 $L_{22}z_{2}$。下面是整个算法:
\\
\\
{{ math:linear_algebra:laff:6.3.2ltrsvpicture.png?400 |}}
\\
\\
===LU分解的 back substitution:Ux = z===
之前解决了 $Lz = b$ 这个方程组,我们得到了现在我们需要求解 $Ux = z$ 。这一步是在求解 $x$ 是多少。按照之前的手法,将参与运算的单位都进行对应的分解,然后根据 $Ux = z$ 列出等式,有:
\\
\\
$$
\left( \begin{array}{c | c}
\nu_{11} & u_{12}^T \\
\hline
0 & U_{22} \\
\end{array}
\right)
\left( \begin{array} {r}
\chi_1 \\
\hline
x_2
\end{array}
\right)
=
\left(\begin{array}{r} \zeta_1 \\ \hline z_2 \end{array}\right)
$$
\\
\\
化简以后有:
\\
\\
$$
\left( \begin{array} {c}
\nu_{11}\chi_1 + u_{12}^Tx_2\\
\hline
U_{22}x_{2}
\end{array}
\right)
=
\left(\begin{array}{c} \zeta_1 \\ \hline z_2 \end{array}\right)
$$
\\
\\
这里跟之前有一些不同。想象一下之前高斯消元中的 back substitution,求解的过程是**自下而上**的。也就是说,我们会先解出 $x_2$,再代入 $x_2$ 解出 $x_1$,以此类推。那么我们再来看看上面的式子,其中第二行的关系是:
\\
$$
U_{22}x_2 = z_2
$$
\\
假设当前的循环开始于这么一个状态:
* 通过上一轮的的循环我们解出了 $x_2$ 部分
* 其他的运算还没有发生
那么作为结果的 $z$ 向量中,对应 $x_2$ 部分的分量 $z_2$,存放的就是上一轮计算出来的 $x_2$。因为此时并没有通过 $U_{22}x_2$ 来更新 $z_2$ 中的内容。因为 $z_2$ 是已知的,那么可以确定的是:
$$z_2 = x_2$$
\\
再将这个结果带入到
$$\nu_{11}\chi_1 + u_{12}^Tx_2 = \zeta_1$$
可以得到:
$$ \chi_1 = (\zeta_1 - u_{12}^Tz_2) / \nu_{11} $$
\\
\\
那么可以看出来,整个循环的过程是这样的:
- 通过上一轮中得到的 $x_2$,我们更新了当前的 $\zeta_1$ 的内容
- 接着再用 $\zeta_1$ 的内容与 $\nu_{11}$ 相除,我们得到了新的 $chi_{1}$
接下来就是算法部分了。我们注意到,在逐步求出 $x$ 的过程中,每一次迭代中我们只会使用对应大小的 $z$ 的分量来求 $x$。这就意味着之前使用的 $z$ 中对应分量不会再使用到了,一个例子就是我们之前所提到的 $x_2 = z_2$。因此我们可以利用 $z$ 的空间来存储 $x$,也就是将每次求出的 $x$ 直接存储到 $z$ 中。那么循环部分的算法实际上就可以写成:
zeta1 = zeta1 - u12t * zeta2 // x1 = zeta1 -u12t *zeta2, zeta1 = x1;
zeta1 = zeta1 / u11 // update chi1
完整的算法如下。注意该算法是自左下往右上的,是一个逐步用求得未知数代入上面的方程,求得新的未知数的过程。
\\
\\
{{ math:linear_algebra:laff:3.2.6.3_symmetrizepicture.png?400 |}}
\\
\\
===LU分解算法的消耗===
* $LU$ 分解的消耗:$(2/3)n^3$
* 计算 $Lz=b$:$n^2$
* 计算 $Ux=z$:$n^2$
====参考资料====
* [[https://www.youtube.com/watch?v=m3EojSAgIao|Solve a System of Linear Equations Using LU Decomposition]]
* 本章所有非 svg 图片来源于 LAFF 课件