What & How & Why

您在这里: start » math » linear_algebra » laff » week_7

More Gaussian Elimination

−目录

LAFF Week 7 Notes

本章的重要内容：

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

When Gaussian Elimination Works

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

那么反过来想，如果主对角线所有元素都不为 00，那么该算法是有解的。鉴于在 Lz=bLz = b 的算法中不包含除法，因此 Ux=zUx=z 有解也就意味着整个 Ax=bAx=b 有解。那么这里有一个问题：Ax=bAx=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
Plain text
来看一下 Lz=0Lz=0 意味着什么：

根据上图，我们发现如果要 Lz=0Lz = 0，那么：

  1. ζ0=0\zeta_0 = 0
  2. λ10ζ0+1∗ζ1=0\lambda_{10}\zeta_0+1*\zeta_1 = 0

以此类推，我们可以得出 zz 必须为 00 的结论。

再来看看通过 z=0z=0，也就是 Uw=0Uw=0 意味着什么：

根据之前的假设，我们假设 Ux=zUx=z 有解，也就是图中所有的 vv 都不为 0，那么有：

  1. vn−1,n−1wn−1=0v_{n-1,n-1}w_{n-1} = 0，因为所有 vv 不为 0，那么 wn−1=0w_{n-1} = 0
  2. vn−2,n−2wn−2+vn−1,n−1wn−1=0v_{n-2,n-2}w_{n-2}+v_{n-1,n-1}w_{n-1} = 0，也就是 wn−2=0w_{n-2} = 0

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

主对角线有 0 时的问题

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

(0110)(x0x1)=(β0β1) \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 的方阵；比如下图的矩阵就是一个置换矩阵：

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

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

(010001100)(−212321−10−3)=(321−10−3−212) \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)

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

<html>

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

</html>

置换矩阵的一般形式

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

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

现在使用 PP 来表示对应 pp 的置换矩阵，那么有：

P(p)=(e~k0Te~k1T⋮e~kn−1T) 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)

假设我们的目的是用 PP 对 AA 中的行为单位排序，那么我们需要将 PP 和 AA 都按照行来划分，因此有： PA=P(p)A=(e~k0Te~k1T⋮e~kn−1T)(a~0Ta~1T⋮a~n−1T)=(a~k0Ta~k1T⋮a~kn−1T) 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)

注意这里 eiTAe^T_iA 代表的就是 AA 的第 ii 行。



而如果需要对 AA 按列为单位排序，为了得到与 AA 维度相同的矩阵，我们需要使用 PTP^T 对 AA 进行右乘：

APT=A(e~k0Te~k1T⋮e~kn−1T)T=A(ek0ek1…ekn−1)=(Aek0Aek1…Aekn−1)=(ak0ak1…akn−1) 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 的情况，我们的处理思路就是将这一行与其他的行进行交换；这样的交换有可能使高斯变换避免计算中除零的情况。

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



可以注意到的是：

  • P~(π)\tilde P( \pi ) 的转置矩阵与其相等
  • P~(π)\tilde P( \pi ) 左乘目标矩阵是换行
  • P~(π)\tilde P( \pi ) 右乘目标矩阵是换列（因为 P~(π)=P~(π)T\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 分解，并求出 pp。该过程手动计算的过程大致如下：

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

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

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

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



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

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

Inverse Matrix

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

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

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

逆矩阵的基本性质

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

特殊矩阵的逆矩阵

详情如下图：



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

参考资料

  • 最后更改: 2021/11/11 08:07
  • ©2022 YGwiki.com