解线性方程组

约 4247 个字预计阅读时间 14 分钟

本章默认 \(n=m\)，先从方阵的消元入手，之后扩展到所有矩阵．

1. Gauss-Jordan 消元

如果我们有上三角矩阵 \(U\)，我们可以从后往前回代得到答案．因为此时第 \(n\) 行退化为一元一次方程 \(a_{nn}x_{n}=b_{n}\)，解出 \(x_{n}\)，再带回到 \(n-1\) 行继续该操作．

Gauss 就是 \(A\to U\) 的过程，而 Jordan 就是 \(U\to R\) 的过程．如果矩阵可逆，这里的 \(R\) 即为对角阵 \(D\)．

Row Echelon Form 行阶梯形矩阵：简称 \(\text{ref}\) 或 \(R\)，表示主元所在列上下方元素均为 \(0\) 的矩阵，其可以通过 \(U\) 经过消元得到．有时 \(R\) 也指 Reduced Row Echelon Form 简化行阶梯形矩阵，其在 \(\text{ref}\) 基础上满足了主元均为 \(1\)．显然，一个矩阵的 \(\text{rref}\) 是唯一的，而 \(\text{ref}\) 不是，因为其的每一行可以有任意的系数．

Gauss的过程：用 \(a_{ii}\) 的元素将其下方所有值消去．例如我们想在三阶矩阵消去 \(a_{21}\)，首先得到 \(l_{21}=\dfrac{a_{21}}{a_{11}}\)，然后将第二行整体减去第一行的 \(l_{21}\) 倍，从而将 \(a_{21}\) 消为 \(0\)．接着按照从上到下、从左到右的顺序对 \(a_{31}\) 与 \(a_{32}\) 进行消元，得到 \(U\)．

在对角线上用于消去其下方元素的元素被称为主元．若出现主元为 \(0\) 时，我们需要与下方该位置非 \(0\) 的行进行行交换，继续消元过程．如果后面已经没有行可以交换，那么 Gauss 过程结束．

得到上三角矩阵 \(U\) 后，使用上述 Jordan 消元过程得到 \(R\) 矩阵．

Gauss:
for i in [1, n]:
    if a[i][i] == 0:
        find k, s.t. k > i and a[k][i] != 0
        swap(row[i], row[k])
    for j in [i + 1, n]:
        l[j][i] = a[j][i] / a[i][i]
        row[j] -= l[j][i] * row[i]

Jordan:
for i in [n, 1]:
    if a[i][i] != 0:
        for j in [i - 1, 1]:
            row[j] -= a[j][i] / a[i][i] * row[i]

我们发现消元过程中只涉及行交换、行倍加操作，而行变换不改变方程组的解且可以用左乘变换矩阵来表示，这为 \(LU\) 分解打下了基础．

2. 矩阵乘法

2.1 视角

见矩阵乘积．

2.2 增广矩阵

在线性代数中我们常做行变换而少做列变换，是因为对于方程组 \(A\boldsymbol{x=b}\)，行变换的过程相当于方程相加减的过程，其不改变方程组的解．

但是在方程相加减的过程中，右侧的 \(\boldsymbol{b}\) 也进行了加减，为了让其与 \(A\) 的行变换同步，我们将 \(A\) 与 \(\boldsymbol{b}\) 拼在一起得到增广矩阵 \([A, \boldsymbol{b}]\)．这样行变换的过程就与解方程完全统一了．

2.3 行变换矩阵

我们知道消元的过程就是行变换的过程，而行变换可以用矩阵做到．回顾上述“行向量乘行向量”视角，想要让新矩阵为原矩阵行向量的线性组合，那么行变换矩阵必须左乘原矩阵；反之，右乘为列变换．

在 Gauss-Jordan 消元中我们需要用到两种矩阵：行倍增矩阵和行交换矩阵．下面来介绍这两种矩阵．

2.3.1 消元矩阵

Elimination Matrix 或 Elementary Matrix，用符号 \(E\) 表示．我们用 \(E_{ij}\) 表示我们想要对矩阵 \(A\) 消去 \(a_{ij}\) 元素．

\(E_{ij}\) 需要让第 \(i\) 行减去第 \(j\) 行的 \(l_{ij}\) 倍而不改变其他值，我们可以用行向量乘行向量的视角来看待：除了第 \(i\) 行外，其他行不做变换，因此 \(E_{ij}\) 需要基于行单位矩阵 \(I\)；在 \(i\) 行处，要让第 \(i\) 行减去第 \(j\) 行的 \(l_{ij}\) 倍，那么 \(e_{ij}\) 处的值就要成为 \(-l_{ij}\)．

例

想让第三行减去 \(3\) 倍的第一行，有

\[ E_{31}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -3 & 0 & 1 \end{bmatrix} \]

\(E_{31}A\) 矩阵的最后一行就是 \(A\) 第一行的 \(-3\) 倍与第三行的 \(1\) 倍的和，达到了消元的目的．

2.3.2 置换矩阵

Permutation Matrix，用符号 \(P\) 表示．\(P_{ij}\) 表示对矩阵第 \(i\) 行与第 \(j\) 行作交换．

仍然是用行向量乘行向量的视角：要让 \(P_{ij}A\) 矩阵的第 \(i\) 行变为 \(A\) 的第 \(j\) 行，只需在 \(I\) 的基础上，将元素 \(p_{ii}\) 设为 \(0\) 而 \(p_{ij}\) 设为 \(1\)，同理 \(p_{jj}\) 为 \(0\) 而 \(p_{ji}\) 为 \(1\)．

例

\[ \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} =\begin{bmatrix} 3 & 4 \\ 1 & 2 \end{bmatrix} \]

2.3.3 对角矩阵

Diagonal Matrix，用符号 \(D\) 表示．其可用于分别将矩阵各行添加不同的系数，一般用于将 \(U\) 矩阵的对角线元素化为 \(1\)．

仍然是用行向量乘行向量的视角：要当 \(DA\) 的矩阵第 \(i\) 行各元素附上系数 \(k\)，只需在单位矩阵的基础上令 \(d_{ii}=k\) 即可．

例

\[ \begin{bmatrix} 2 & 4 \\ 0 & 5 \end{bmatrix}= \begin{bmatrix} 2 & 0 \\ 0 & 5 \end{bmatrix}\cdot \begin{bmatrix} 1 & 2 \\ 0 & 1 \end{bmatrix} \]

2.4 运算法则

对于矩阵乘法而言，其不满足交换律，满足结合律、分配律．

2.5 分块矩阵

对矩阵用横线、竖线将矩阵分块，可以得到分块矩阵．对于矩阵 \(A,B\)，如果： + \(A,B\) 本身满足矩阵乘法条件 + \(A,B\) 的分块矩阵的切分结构满足矩阵数量的条件（如 \(A\) 按宽度分成 \(3\) 块，而 \(B\) 按高度分成 \(3\) 块） + \(A,B\) 的各个小块满足矩阵乘法的条件

那么 \(A,B\) 就可以进行分块矩阵乘法，其得到的结果与原矩阵乘法相同．

舒尔补（Schur Complement）：对分块矩阵进行 Gauss 消元操作．对于矩阵

\[ M= \begin{bmatrix} A & B \\ C & D \end{bmatrix} \]

要想消去 \(C\)，第二行需要减去 \(CA^{-1}Row_{1}\)，即

\[ \begin{bmatrix} I & 0 \\ -CA^{-1} & I \end{bmatrix} \cdot \begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} A & B \\ 0 & D-CA^{-1}B \end{bmatrix} \]

其中 \(S=D-CA^{-1}B\) 就称为 \(A\) 在 \(M\) 中的舒尔补．

3. 逆矩阵

我们希望得到一个变换 \(A^{-1}\)，使得其与 \(A\) 的变换完全相反，即 \(A^{-1}A=I\)，这样就可以通过给 \(A\boldsymbol{x=b}\) 左乘 \(A^{-1}\) 得到结果 \(\boldsymbol{x=}A^{-1}\boldsymbol{b}\)．

3.1 性质

我们先来探索一下逆矩阵性质．

首先是左右逆等效．若 \(BA=I\) 且 \(AC=I\)，则 \((BA)C=B(AC) \implies C=B\)，即左逆等于右逆，因此 \(A^{-1}A\) 与 \(AA^{-1}\) 等价．同时我们也证明了：一个矩阵至多有一个逆矩阵．显然，矩阵可逆的必要条件是方阵，不然 \(A^{-1}A\) 与 \(AA^{-1}\) 的阶数不同．
矩阵乘积的逆．可以将矩阵乘法看作栈，对于矩阵 \(AB\)，想要将其变换为 \(I\)，首先需要消除栈顶的 \(A\)，因此先对其左乘 \(A^{-1}\)；此时栈中剩下 \(B\)，再次左乘 \(B^{-1}\) 即可，得到 \((AB)^{-1} =B^{-1}A^{-1}\)．
对于上/下三角矩阵而言，其逆矩阵依然是上/下三角矩阵．以上三角矩阵为例，考虑逆矩阵对其右乘：逆矩阵的第一列为原矩阵列向量的线性组合系数，该线性组合得到单位向量的第一列；显然逆矩阵的元素 \(A^{-1}_{n1}\) 必须为 \(0\)，否则 \(A\) 的第 \(n\) 列加入后，得到的列向量最后一个分量不为 \(0\)（其他向量最后一个分量均为 \(0\)，无法将其消去）．然后递归得到 \(A^{{-1}}_{(n-1)1}\) 为 \(0\)，直到 \(A^{-1}_{21}\) 都为 \(0\)，第一列只有 \(A^{-1}_{11}\) 不为 \(0\)．其他列同理．

3.2 计算

同时对矩阵 \(A\) 与 \(I\) 进行 Gauss-Jordan，试图将 \(A\) 的部分化为 \(I\)．若可以化为 \(I\)，说明行变换矩阵等价于 \(A^{-1}\)，此时对右侧 \(I\) 而言其相当于 \(A^{-1}I=A^{-1}\)，即为其逆矩阵．

显然，如果 \(A\) 的部分无法化为 \(I\)，则不存在逆矩阵．而 \(A\) 的部分能化为 \(I\) 的充要条件为 \(A\) 在进行 Gauss过程中每一行都存在主元．因此若 \(A\) 的主元数小于 \(n\)，其就是不可逆矩阵．

4. A的LU分解

注：在讨论 \(LU\) 分解时，默认 \(A\) 可逆，并且在进行 Gauss 操作时无需进行行置换，也就是只有 \(E\) 的参与；若需要行置换，我们可以用置换矩阵 \(P\) 将所需要的行置换提前换好，得到矩阵 \(PA\)，保证 \(PA\) 在 Gauss 操作时无需进行行置换．

为什么我们需要LU分解

在实际应用中，我们很可能需要对同一个 \(A\)、不同的 \(\boldsymbol{b}\) 解许多次 \(A\boldsymbol{x=b}\)．我们希望不用每次都对 \([A,\boldsymbol{b}]\) 这个增广矩阵重新做 Gauss 操作，而是对 \(A\) 进行一次 Gauss 操作后得到 \(U\)，然后将所做的行变换记录下来．当所给的 \(\boldsymbol{b}\) 确定时，直接对其做我们记录下的行变换得到 \(U\boldsymbol{x=c}\)，然后用 Jordan解出 \(\boldsymbol{x}\) 即可．

那么我们应该去哪里寻找我们做过操作的记忆矩阵呢？

考虑如下两个消元矩阵

\[ E_{21}=\begin{bmatrix} 1 & 0 & 0 \\ -5 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\quad E_{32}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -4 & 0 \end{bmatrix} \]

由于矩阵的消元是从右往左操作，因此与这两个消元矩阵等效的操作矩阵为

\[ E=E_{32}E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -4 & 0 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 \\ -5 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ -5 & 1 & 0 \\ -20 & -4 & 0 \end{bmatrix} \]

我们发现出现了复杂的数字 \(-20\)，这不是我们想要的．

为什么会出现这种结果？因为我们是按照“从上到下”的顺序进行消元的，而后面步骤的消元会用到之前消元得到的结果．类似于动态规划中，我们使用了同一轮更新后的数据而不是上一轮的数据．

也就是说想要解决这个问题，只需要按照“从下到上”的顺序进行消元即可，即将我们所有的消元矩阵顺序全部反转．而什么操作可以将操作顺序反转而不影响结果？没错，就是逆矩阵．

在此，先证明一个引理：

消元矩阵的逆

对于消元矩阵 \(E_{ij}\)，若想找到其逆矩阵，可以从实际意义出发：当 \(E_{ij}\) 一个矩阵后，它会使该矩阵第 \(i\) 行加上自身第 \(j\) 行的 \(e_{ij}\) 倍；而想要抵消该操作，当然是把加上的这行减回去，即逆矩阵满足元素 \(e^{-1}_{ij}=-e_{ij}\)．

回代证明后也正确

\[ E_{32}\cdot E_{32}^{-1}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -l_{32} & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & l_{32} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} =I \]

由于假设消元过程中不存在行交换，以 \(3\) 阶矩阵为例，我们有

\[ (E_{32}E_{31}E_{21})A=U \]

等式两边乘以 \((E_{32}E_{31}E_{21})^{-1}\) 得到

\[ A=(E_{21}^{-1}E_{31}^{-1}E_{32}^{-1})U \]

而 \(E_{21}^{-1}E_{31}^{-1}E_{32}^{-1}\) 就是我们要找的记忆矩阵 \(L\)．看看它有什么特点：

\[ \begin{aligned} L&=E_{21}^{-1}E_{31}^{-1}E_{32}^{-1}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ l_{31} & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & l_{32} & 1 \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ l_{31} & l_{32} & 1 \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{bmatrix} \end{aligned} \]

由于逆矩阵按照相反的顺序进行行变换，先做的行变换不会影响到后做的行变换结果，顺序确保了每次用以倍加的行从未被修改过；因此我们得到了一个漂亮的下三角矩阵 \(L\)：它的主对角线元素都是 \(1\)，而下方对应位置元素都是 Gauss 消元时用到的倍数．

进一步化简

我们可以用对角阵对 \(A\) 进行进一步化简．\(L\) 的主对角线元素均为 \(1\)，而 \(U\) 的主对角线元素很乱．根据之前提到的对角矩阵化简，将 \(U\) 拆分成 \(DU'\)，我们可以得到 \(A=LDU'\)，其中 \(L\) 与 \(U'\) 的主对角线元素均为 \(1\)，即单位三角矩阵．

\[ \text{Split U into} \begin{bmatrix} d_{1} & & & \\ & d_{2} & & \\ & & \ddots & \\ & & & d_{n} \end{bmatrix} \cdot \begin{bmatrix} 1 & u_{12}/d_{1} & u_{13}/d_{1} & \\ & 1 & u_{23}/d_{2} & \\ & & \ddots & \vdots \\ & & & 1 \end{bmatrix} \]

LDU矩阵的好处是：对于可逆矩阵 \(A\)，其分解为 \(LDU\) 的方式是唯一的．假设存在不用的两组分解 \(A=L_{1}D_{1}U_{1}=L_{2}D_{2}U_{2}\)，由于 \(A\) 可逆则子矩阵均可逆，可得到

\[ (L_{2}^{-1}L_{1})D_{1}=D_{2}(U_{2}U_{1}^{-1}) \]

由于 \(L_{2}\) 为单位下三角矩阵，其逆 \(L_{2}^{-1}\) 也为单位下三角矩阵，乘积 \(L_{2}^{-1}L_{1}\) 也为单位下三角矩阵，在乘以对角矩阵 \(D_{1}\)，得到的是下三角矩阵．同理，右边的结果是上三角矩阵．

下三角矩阵 \(=\) 上三角矩阵，说明为对角矩阵．而单位下三角矩阵与对角矩阵的乘积为对角矩阵，说明这个单位下三角矩阵只有对角线上有元素，即其为对角矩阵．即为单位下三角矩阵又为对角矩阵，说明 \(L_{2}^{-1}L_{1}\) 为单位矩阵．因此 \(L_{1}=L_{2}\)，同理 \(D_{1}=D_{2}\)，即只有一个分解方式．

当我们得到 \(A=LU\) 后，原来的方形方程组就化简为了两个三角形方程组．而三角形方程组是很容易解的．

解 \(A\boldsymbol{x=b}\)，等价于解

\[ LU\boldsymbol{x=b} \]

根据结合律，我们可以改变括号的位置

\[ L(U\boldsymbol{x})=\boldsymbol{b} \]

引入中间变量 \(\boldsymbol{c}=U\boldsymbol{x}\)，等价于解两个三角形方程组

\[ \begin{cases} L\boldsymbol{c=b}\\ U\boldsymbol{x=c} \end{cases} \]

两个三角形方程组的解法，此处就不赘述了．

时间复杂度计算

由于乘法的计算时间远大于加法的计算时间，因此我们以乘法的计算次数为准．

首先考虑对 \(A\) 的 Gauss 操作．首先我们要消去第一列除 \(a_{11}\) 外的所有元素，因此我们要消去 \(n-1\) 行；每次消去操作都是对一整行的元素进行的，因此第一轮消去需要进行 \(n(n-1)\) 步，渐近看作 \(n^{2}\) 步；接下去的过程为递归消元 \(n-1\) 阶矩阵，需要 \((n-1)^{2}\) 步，总步数为

\[ \sum_{1}^{n}n^{2}=\frac{n(n+1)(2n+1)}{6} \]

复杂度为 \(O(n^{3})\)．

接着考虑对 \(\boldsymbol{b}\) 的操作．此时 \(A=LU\) 已知，计算 \(\boldsymbol{b}\) 的复杂度相当于两个三角形方程组的和，即计算单个三角形方程组的复杂度．以 \(U\boldsymbol{x=c}\) 为例，首先消去最后一列除 \(u_{nn}\) 的元素，但此时只需要对单元素操作而不需要对整行操作，时间复杂度 \(O(n)\)，总复杂度 \(O(n^2)\)．

因此，主要的复杂度就在 \(LU\) 分解上，后续对于每个 \(\boldsymbol{b}\) 的计算是很快速的．

5. 转置矩阵

将矩阵 \(A\) 转置（Transpose）得到 \(A^{T}\)，满足 \(a_{ij}=a^{T}_{ji}\)．

5.1 性质

和的转置：\((A+B)^{T}=A^{T}+B^{T}\)，因为只是元素位置换了而已．
乘积的转置：\((AB)^{T}=B^{T}A^{T}\)．将 \(B\) 拆分成列向量，对于单个列向量 \(\boldsymbol{x}\) 而言，\(A\boldsymbol{x}\) 为 \(A\) 列向量的线性组合，其转置结果应该为 \(A^{T}\) 行向量的线性组合，即 \(\boldsymbol{x}^{T}A^{T}\)，因此 \(A\boldsymbol{x}=\boldsymbol{x}^{T}A^{T}\)，推广得到 \((AB)^{T}=B^{T}A^{T}\)．
逆矩阵的转置：\((A^{-1})^{T} = (A^{T})^{-1}\)．对 \(AA^{-1}=I\) 两边同时取转置得到 \((A^{-1})^{T}A^{T}=I\)，由逆矩阵的唯一性得到 \((A^{-1})^{T} = (A^{T})^{-1}\)．

5.2 应用

一般而言，我们主要是用转置将列向量转为行向量．如 \(\boldsymbol{x^{T}y}\) 表示行向量与列向量的乘积，即内积；而 \(\boldsymbol{xy^{T}}\) 表示列向量与行向量的乘积，即外积．

6. 对称矩阵

转置矩阵与自身相等的矩阵为对称矩阵 Symmetric Matrix，用 \(S\) 表示．

\((S^{-1})^{T}=(S^{T})^{-1}=S^{-1}\)，即 \(S\) 的逆矩阵也为对称矩阵（如果存在）．

对于任意矩阵 \(A\)，我们可以用其生成对称矩阵 \(S\)：\(AA^{T}\) 与 \(A^{T}A\)．不难证明这两个矩阵均为对称矩阵．

当我们对 \(S\) 进行 \(LDU\) 分解时，我们得到

\[ \begin{cases} S=LDU \\ S^{T}=(LDU)^{T}=U^{T}D^{T}L^{T}=U^{T}DL^{T} \end{cases} \]

即 \(LDU=U^{T}DL^{T}\)，由于 \(LDU\) 分解是唯一的，得到 \(L=U^{T}\)，因此 \(S=LDL^{T}\)．