What & How & Why

More Gaussian Elimination

LAFF Week 7 Notes


本章的重要内容:

  • 置换矩阵的概念
  • 使用置换矩阵处理 $LU$ 分解中可能存在的对角线元素为 0 的问题
  • 非奇异矩阵与逆矩阵
  • 逆矩阵的概念及其性质

When Gaussian Elimination Works

先来回顾一下 $LU$ 分解中 关于 $Ux = z$ 求解的循环部分,我们发现其中有这么一个过程:
$$\beta_{21} = \beta_{21} / v_{11}$$
这一部分代表了对 multiplier 的计算。如果仔细观察,我们注意到如果 $v_{11}=0$,也就是主对角线元素为 $0$ 的时候,该算法就无法进行下去了。

那么反过来想,如果主对角线所有元素都不为 $0$,那么该算法是有解的。鉴于在 $Lz = b$ 的算法中不包含除法,因此 $Ux=z$ 有解也就意味着整个 $Ax=b$ 有解。那么这里有一个问题:$Ax=b$ 的解是唯一的吗?

我们试着通过反证法来证明一下:

/* part 1*/
假设:Ax = b 有两个解 u 和 v
那么根据 Ax = b,有
Au = b
Ax = b
因为 A 是 矩阵(线性变换)
因此 A (u - v) = b => Au - Av = b - b = 0
令 u - v = w,即可得 Aw = 0
通过 LU 分解, A = LU
那么 Aw = (LU)w = L(Uw) = 0
来看一下 $Lz=0$ 意味着什么:

根据上图,我们发现如果要 $Lz = 0$,那么:

  1. $\zeta_0 = 0$
  2. $\lambda_{10}\zeta_0+1*\zeta_1 = 0$

以此类推,我们可以得出 $z$ 必须为 $0$ 的结论。

再来看看通过 $z=0$,也就是 $Uw=0$ 意味着什么:

根据之前的假设,我们假设 $Ux=z$ 有解,也就是图中所有的 $v$ 都不为 0,那么有:

  1. $v_{n-1,n-1}w_{n-1} = 0$,因为所有 $v$ 不为 0,那么 $w_{n-1} = 0$
  2. $v_{n-2,n-2}w_{n-2}+v_{n-1,n-1}w_{n-1} = 0$,也就是 $w_{n-2} = 0$

依次类推,可以得出 $w=0$ 的结论,也就是 $u=v$;因此可以推断出只要通 $LU$ 分解求解中,如果主对角线均不为 0 , $Ax=b$ 必定有且只有一个唯一解

主对角线有 0 时的问题

当 $U$ 的主对角线上有 0 的时候,麻烦就来了。根据 $Ux=z$ 的算法,我们不可避免的遭遇到除 0 的情况。该情况会导致算法无法继续进行;但问题在于,有时算法的无法进行并不代表我们无法得出解。下面是一个很明显的例子:

$$ \displaystyle \left( \begin{array}{c c} 0 & 1 \\ 1 & 0 \end{array} \right) \left( \begin{array}{c} x_0\\ x_1 \end{array} \right) = \left( \begin{array}{c} \beta_0\\ \beta_1 \end{array} \right) $$
上式很显然有解,但我们并不能通过之前学习的算法来求解。因此,我们需要引进置换矩阵的概念来解决这个问题。

置换矩阵 Permutations Matrix

置换矩阵(Permutations Matrix)可以简要的概述为每行每列有且只有一个 1,其他元素都为 0 的方阵;比如下图的矩阵就是一个置换矩阵:

$$\left(\begin{array}{c}0&1&0\\0&0&1\\1&0&0\end{array}\right)$$

置换这个名字非常形象。以上面的置换矩阵为例,来看一个矩阵的乘法:

$$ \left(\begin{array}{c}0&1&0\\0&0&1\\1&0&0\end{array}\right) \left(\begin{array}{c}-2&1&2\\3&2&1\\-1&0&-3\end{array}\right) = \left(\begin{array}{c}3&2&1\\-1&0&-3\\-2&1&2\end{array}\right) $$

这个例子有趣的地方在于,如果将置换矩阵改为单位矩阵 $I$,再对两个乘法进行对比,我们会发现其实置换矩阵改变的是目标矩阵行的排序,而该排序正好对应置换矩阵中行之余单位矩阵中的行的位置:

<html>

<img src=“/_media/math/linear_algebra/laff/premult_mat.svg” width=“450”>

</html>

置换矩阵的一般形式

大致知道置换矩阵的功能之后,我们来考虑一下置换矩阵的一般形式。我们注意到,无论是从行还是列来看待置换矩阵,置换矩阵都可以看作是一个由 Unit basis vector 组成的向量。如果我们使用向量 $p$ 代表置换矩阵中 Unit basis vector 的顺序,那么 $p$ 可以表示为:

$$ p = \left( \begin{array}{c} k_0, \cdots, k_{n-1} \end{array} \right)^T $$

现在使用 $P$ 来表示对应 $p$ 的置换矩阵,那么有:

$$ P(p) = \left( \begin{array}{c} \widetilde e_{k_0}^T \\ \hline \widetilde e_{k_1}^T \\ \hline \vdots \\ \hline \widetilde e_{k_{n-1}}^T \end{array} \right) $$

假设我们的目的是用 $P$ 对 $A$ 中的为单位排序,那么我们需要将 $P$ 和 $A$ 都按照行来划分,因此有: $$ PA = P(p)A = \left( \begin{array}{c} \widetilde e_{k_0}^T \\ \hline \widetilde e_{k_1}^T \\ \hline \vdots \\ \hline \widetilde e_{k_{n-1}}^T \end{array} \right) \left( \begin{array}{c} \widetilde a_0^T \\ \hline \widetilde a_1^T \\ \hline \vdots \\ \hline \widetilde a_{n-1}^T \end{array} \right) = \left( \begin{array}{c} \widetilde a_{k_0}^T \\ \hline \widetilde a_{k_1}^T \\ \hline \vdots \\ \hline \widetilde a_{k_{n-1}}^T \end{array} \right) $$

注意这里 $e^T_iA$ 代表的就是 $A$ 的第 $i$ 行。



而如果需要对 $A$ 按为单位排序,为了得到与 $A$ 维度相同的矩阵,我们需要使用 $P^T$ 对 $A$ 进行右乘

$$ AP^T = A \left( \begin{array}{c} \widetilde e_{k_0}^T \\ \hline \widetilde e_{k_1}^T \\ \hline \vdots \\ \hline \widetilde e_{k_{n-1}}^T \end{array} \right)^T =A\left( \begin{array}{c|c|c|c} e_{k_0} & e_{k_1} & \dots & e_{k_{n-1}} \end{array} \right)\\ =\left( \begin{array}{c|c|c|c} Ae_{k_0} &A e_{k_1} & \dots & Ae_{k_{n-1}} \end{array} \right) =\left( \begin{array}{c|c|c|c} a_{k_0} & a_{k_1} & \dots & a_{k_{n-1}} \end{array} \right) $$

使用转置矩阵解决除零问题

通过上面对转置矩阵的了解,我们知道可以通过转置矩阵控制目标矩阵中行(列)的位置。这一点对于解决高斯变换(LU 分解)中除零的问题是非常有效的。对于对角线中元素是 0 的情况,我们的处理思路就是将这一行与其他的行进行交换;这样的交换有可能使高斯变换避免计算中除零的情况。

理想中对目标矩阵行的交换只牵涉到两行,因此我们需要引进一种更加具体的置换矩阵 $\tilde P( \pi )$,我们称之为 Pivot Matrix(中文: 枢轴矩阵,翻译可能不准确)。在 $\tilde P( \pi )$ 中,只有第 $\pi$ 行与第一行发生了交换,其他行都保持不变。具体的演示如下图:



可以注意到的是:

  • $\tilde P( \pi )$ 的转置矩阵与其相等
  • $\tilde P( \pi )$ 左乘目标矩阵是换行
  • $\tilde P( \pi )$ 右乘目标矩阵是换列(因为 $\tilde P( \pi ) = \tilde P( \pi )^T$ )

算法:LU Factorization with Partial Pivoting

该算法基本上高斯变换的算法相同,只需要添加一个条件分支即可:

<html>

<img src=“/_media/math/linear_algebra/laff/pivot_al.svg” width=“250”>

</html>

总的说来,本算法分两个大步骤。第一个步骤是对目标矩阵进行 LU 分解,并求出 $p$。该过程手动计算的过程大致如下:

  1. 每一个循环都需要应用和更新代表高斯变换的矩阵($L$)。
  2. 查看目标对角线上的当前循环的元素,如果不为 0,无需交换行。也就是说使用 $I$ 作为置换矩阵;并将 $p$ 中对应的分量设置为 0($p$ 对应分量负责存储 $\pi$ 的值,0 代表当前没有置换的行)
  3. 如果当前循环对角线上元素为 0,查看当前行下方,选取当前行正下方元素不为 0 的第一行进行交换(应用对应的置换矩阵),并记录下用于交换行的行数到 $p$ 对应的分量中。
  4. 如此反复,直到完成 $LU$ 分解。

当完成这部分计算以后,第二个步骤是使用 $p$ 来对 $b$ 进行更新。简单的说来,因为换过行,因此不但 $A$ 中的行要随之变化,对应的 $b$ 中的分量也要随之变化。(这个很好理解:某一个方程等式左边和右边应该视作一个整体一起交换)。也就是说,

因为 Ax = b
所以 PAx = Pb //保持相等

整个完整的算法如下图所示:



跟之前的 LU 分解相比,上图唯一的变化就是先求出 $p$,然后使用 $p$ 更新 $b$。接着就是我们熟悉的两步: $Lz=b$ 与 $Ux = z$ 了。

当然,我们还会遇到一种情况:矩阵中根本就找不到合适的一行与当前含 0 行进行替换。这象征着 $Ax=b$ 无解或者有无限多解,该内容将在第八章中详述。

Inverse Matrix

在高中我们都学过反函数。反函数可以理解为一个映射的逆映射,也就是还原当前映射。令原函数为 $f$,反函数为 $f^{-1}$,则下方数学表达式表达了如下关系:
$$ f^{-1}(f(x)) = x $$
很显然,原函数中的自变量和因变量必须要一一对应,才能完成这样的“还原”。

下面再来看看矩阵。根据之前的知识,我们知道矩阵也是映射的一种。与反函数相似,如果存在这么一个矩阵,可以还原某个矩阵表示的映射,那么我们就称这个矩阵是原矩阵的逆矩阵。而原矩阵,我们也称为非奇异矩阵Nonsingular Matrix)。

当然,相较于 1D 函数,矩阵是向量的函数,映射前与映射后维度可能发生变化 $f: \mathbb{R}^n \to \mathbb{R}^m)$。但因为逆矩阵必须要求映射的定义域与值域有着一一对应的关系,那么只有当 $m=n$,也就是原矩阵是方阵的时候,才有可能存在逆矩阵。

逆矩阵的基本性质

  1. 逆矩阵也是线性变换
  2. $AA^{-1} = I$
  3. $(aB)^{-1} = \frac{1}{a}B^{-1} $
  4. $(AB)^{-1} = B^{-1}A^{-1}$
  5. $(ABC)^{-1} = C^{-1}B^{-1}A^{-1}$
  6. $(A^T)^{-1} = (A^{-1})^T$

特殊矩阵的逆矩阵

详情如下图:



当然,不管是什么矩阵,求其逆矩阵的主要方法都是通过 $AA^{-1} = I$ 这个关系来求。进一步的说,我们可以将 $I$ 按列划分,按列得出关系,一列一列的求逆矩阵。(例子待添加)FIXME

参考资料