LU 分解和矩阵实现

2024-5-1

计算机图形学

这一节我们实现矩阵运算。矩阵的乘法等操作都是很好实现的，唯独矩阵的逆比较复杂。之前大学线性代数的课程中，矩阵求逆都是用解线性方程组的方式进行手算的。但是觉得“模拟”的方式不太适合计算机运行，特别涉及行交换的时候。

查资料选了 LU 分解的方式求逆矩阵。选择它的原因是，LU 分解本质和手算使用的高斯消去法是一样的，让我觉得应该不会太难，先够用就行。

当然，空间变换时维度低的矩阵进行求逆，可以直接使用伴随矩阵的方法。不过为了通用性，决定采用 LU 分解实现求逆。

1. 高斯消去法：矩阵变换

我们先以矩阵变换的角度来表示高斯消去法。比如，矩阵使用高斯消去法，首先消去第一列，即：第一行乘以 -2 加到第二行，第一行乘以 -4 加到第三行，第一行乘以 -3 加到第四行。写成矩阵变换的形式为：

\( L_1A= \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ -2 & 1 & \phantom{0} & \phantom{0} \\ -4 & \phantom{0} & 1 & \phantom{0}\\ -3 & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 1 & 0 \\ 4 & 3 & 3 & 1 \\ 8 & 7 & 9 & 5\\ 6 & 7 & 9 & 8 \end{pmatrix} = \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 3 & 5 & 5\\ 0 & 4 & 6 & 8 \end{pmatrix} \)

同理，消去第二列：

\( L_2(L_1A)= \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ \phantom{0} & -3 & 1 & \phantom{0}\\ \phantom{0} & -4 & \phantom{0} & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 3 & 5 & 5\\ 0 & 4 & 6 & 8 \end{pmatrix} = \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 2 & 2\\ 0 & 0 & 2 & 4 \end{pmatrix} \)

同理，消去第三列：

\( L_3(L_2L_1A)= \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ \phantom{0} & \phantom{0} & 1 & \phantom{0}\\ \phantom{0} & \phantom{0} & -1 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 2 & 2\\ 0 & 0 & 2 & 4 \end{pmatrix} = \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 2 & 2\\ 0 & 0 & 0 & 2 \end{pmatrix} \)

这边最终的消去结果就是 \(U\)。就是上三角矩阵的缩写，Upper triangular matrix。

现在我们看一下 \(L_1\)、\(L_2\)、\(L_3\) 的性质。\(L_1\) 的逆为：

\( L_1^{-1}= \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ 2 & 1 & \phantom{0} & \phantom{0} \\ 4 & \phantom{0} & 1 & \phantom{0}\\ 3 & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \)

可以看到是“系数是相反数”的形式。同样 \(L_2\)、\(L_3\) 也是。

直观上也容易理解：第二行已经加上了第一行乘以 -2 的值，现在求逆，即需要把它恢复回去，只需要第二行再加上了第一行乘以 2 的值。

我们再看一个很好的性质：

\(L_1^{-1}L_2^{-1}= \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ 2 & 1 & \phantom{0} & \phantom{0} \\ 4 & 3 & 1 & \phantom{0}\\ 3 & 4 & \phantom{0} & 1 \end{pmatrix} \)

可以发现就是“内容的叠加”。

直观上也容易理解：\(L_1^{-1}L_2^{-1}\) 就是 \(L_2^{-1}\) 作用上 \(L_1^{-1}\) 的行变化。

\(L_1^{-1}\) 只有第一行中，只有第一列不为0，所以变化只会影响到 \(L_2^{-1}\) 的第一列。

而 \(L_2^{-1}\) 的第一列全为 0，所以反应在结果上就是“内容的叠加”。

同样 \(L_1^{-1}L_2^{-1}L_3^{-1}\) 也是这样。现在我们已经有

\(L_3L_2L_1A=U\)

用逆做变换

\(L_1^{-1}L_2^{-1}L_3^{-1}L_3L_2L_1A=L_1^{-1}L_2^{-1}L_3^{-1}U\)

\(A=L_1^{-1}L_2^{-1}L_3^{-1}U\)

这边的 \(L_1^{-1}L_2^{-1}L_3^{-1}\)，我们记作 \(L\)，即下三角矩阵的缩写，Lower triangular matrix。这就是 LU 分解名称的由来，它把一个矩阵分解成下三角和上三角矩阵相乘的形式。

1.1 LU 分解的代码实现

因为 LU 分解的性质非常好，代码实现也很简洁。如代码清单 1 所示，U 基于原矩阵，做高斯消去；L 基于单位矩阵，记录各个高斯消去系数。

代码清单 1 LU 分解
  1. void LUDecomposition(Matrix* L, Matrix* U)
  2. {
  3.     assert(ROWS == COLS);
  4.  
  5.     *U = *this;
  6.     L->ToIdentity();
  7.  
  8.     for (size_t k = 0; k < ROWS; k++)
  9.     {
  10.         for (size_t i = k + 1; i < ROWS; i++)
  11.         {
  12.             assert(U->at(k, k) != 0);
  13.             L->at(i, k) = U->at(i, k) / U->at(k, k);
  14.  
  15.             for (size_t j = k; j < ROWS; j++)
  16.                 U->at(i, j) -= L->at(i, k) * U->at(k, j);
  17.         }
  18.     }
  19. }

2. PLU 分解

之前的 LU 分解还有一些缺陷。比如，对角线上有元素为 0，0 乘以任何数都为 0，消去不了任何项；对角线上的元素很小，消去的时候大数除以小数，影响精度。

针对以上问题，我们可以改进选择消去法时主元的策略：我们可以选取绝对值最大的数当作主元。

我们看到例子。首先交换第一行和第三行，有置换矩阵 \(P_1\)

\( \begin{pmatrix} \phantom{0} & \phantom{0} & 1 & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 1 & 0 \\ 4 & 3 & 3 & 1 \\ 8 & 7 & 9 & 5\\ 6 & 7 & 9 & 8 \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 4 & 3 & 3 & 1 \\ 2 & 1 & 1 & 0 \\ 6 & 7 & 9 & 8 \end{pmatrix} \)

置换矩阵形式上就是单位矩阵交换行。

接着进行消元，有变换矩阵 \(L_1\)

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ -\frac{1}{2} & 1 & \phantom{0} & \phantom{0} \\ -\frac{1}{4} & \phantom{0} & 1 & \phantom{0}\\ -\frac{3}{4} & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \begin{pmatrix} 8 & 7 & 9 & 5 \\ 4 & 3 & 3 & 1 \\ 2 & 1 & 1 & 0 \\ 6 & 7 & 9 & 8 \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & -\frac{1}{2} & -\frac{3}{2} & -\frac{3}{2} \\ 0 & -\frac{3}{4} & -\frac{5}{4} & -\frac{5}{4} \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \end{pmatrix} \)

接着交换第二行和第四行交换，有置换矩阵 \(P_2\)

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & \phantom{0} & \phantom{0} & 1 \\ \phantom{0} & \phantom{0} & 1 & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \end{pmatrix} \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & -\frac{1}{2} & -\frac{3}{2} & -\frac{3}{2} \\ 0 & -\frac{3}{4} & -\frac{5}{4} & -\frac{5}{4} \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & -\frac{3}{4} & -\frac{5}{4} & -\frac{5}{4} \\ 0 & -\frac{1}{2} & -\frac{3}{2} & -\frac{3}{2} \end{pmatrix} \)

接着进行消元，有变换矩阵 \(L_2\)

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ \phantom{0} & \frac{3}{7} & 1 & \phantom{0}\\ \phantom{0} & \frac{2}{7} & \phantom{0} & 1 \end{pmatrix} \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & -\frac{3}{4} & -\frac{5}{4} & -\frac{5}{4} \\ 0 & -\frac{1}{2} & -\frac{3}{2} & -\frac{3}{2} \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & 0 & -\frac{2}{7} & \frac{4}{7} \\ 0 & 0 & -\frac{6}{7} & -\frac{2}{7} \end{pmatrix} \)

接着交换第三行和第四行交换，有置换矩阵 \(P_3\)

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ \phantom{0} & \phantom{0} & \phantom{0} & 1\\ \phantom{0} & \phantom{0} & 1 & \phantom{0} \end{pmatrix} \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & 0 & -\frac{2}{7} & \frac{4}{7} \\ 0 & 0 & -\frac{6}{7} & -\frac{2}{7} \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & 0 & -\frac{6}{7} & -\frac{2}{7} \\ 0 & 0 & -\frac{2}{7} & \frac{4}{7} \end{pmatrix} \)

接着进行消元，有变换矩阵 \(L_3\)

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0}\\ \phantom{0} & 1 & \phantom{0} & \phantom{0} \\ \phantom{0} & \phantom{0} & 1 & \phantom{0}\\ \phantom{0} & \phantom{0} & -\frac{1}{3} & 1 \end{pmatrix} \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & 0 & -\frac{6}{7} & -\frac{2}{7} \\ 0 & 0 & -\frac{2}{7} & \frac{4}{7} \end{pmatrix} = \begin{pmatrix} 8 & 7 & 9 & 5 \\ 0 & \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\ 0 & 0 & -\frac{6}{7} & -\frac{2}{7} \\ 0 & 0 & 0 & \frac{2}{3} \end{pmatrix} \)

最终的形式如下

\(L_3P_3L_2P_2L_1P_1A=U\)

但是以上这形式，既不确定 \(L_3P_3L_2P_2L_1P_1\) 是否是下三角矩阵，也不好求逆。因为 \(P^{-1}=P\)，我们重新构造一下：

\(L_3P_3L_2(P_3P_3)P_2L_1(P_2P_3P_3P_2)P_1A=U\)

重新组合一下

\(L_3(P_3L_2P_3)(P_3P_2L_1P_2P_3)(P_3P_2P_1)A=U\)

我们令

\(P=P_3P_2P_1\)

\(\widetilde{L}_2=P_3L_2P_3\)

\(\widetilde{L}_1=P_3P_2L_1P_2P_3\)

这边 \(\widetilde{L}_1\) 、\(\widetilde{L}_2\) 的形式和 \(L_1\) 、\(L_2\) 的形式是一样的。我们以 \(\widetilde{L}_1\) 为例子感受一下。

左乘是行交换，右乘是列交换。\(P_2\) 作用于第二行和第四行，首先行交换，得

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0} \\ -\frac{3}{4} & \phantom{0} & \phantom{0} & 1 \\ -\frac{1}{4} & \phantom{0} & 1 & \phantom{0}\\ -\frac{1}{2} & 1 & \phantom{0} & \phantom{0} \end{pmatrix} \)

再交换列，得

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0} \\ -\frac{3}{4} & 1 & \phantom{0} & \phantom{0} \\ -\frac{1}{4} & \phantom{0} & 1 & \phantom{0}\\ -\frac{1}{2} & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \)

可以发现又重新变成了行变化的形式，只是第一列中局部的行进行了交换。

\(P_3\) 作用于第三行和第四行。同理，交换行列，得

\( \begin{pmatrix} 1 & \phantom{0} & \phantom{0} & \phantom{0} \\ -\frac{3}{4} & 1 & \phantom{0} & \phantom{0} \\ -\frac{1}{2} & \phantom{0} & 1 & \phantom{0}\\ -\frac{1}{4} & \phantom{0} & \phantom{0} & 1 \end{pmatrix} \)

同样也只是第一列中局部的行进行了交换。一般的有

\(\widetilde{L}_k=P_{n-1}\cdots P_{k+1}L_kP_{k+1}\cdots P_{n-1}\)

从下标中我们可以感受到，因为是方阵，\(L_k\) 对于比它下标大的操作，比如 \(P_{k+1}\)，对于行索引，作用范围肯定是大于 k。这就表示，作用不到单位矩阵对角线上的最小索引的 1 位置；因为行变换，破坏了对角线上 1 的位置，但是列变化也会换过来，对于列索引，同样作用范围肯定是大于 k，这时除了对角线上的 1 以外，其余都是 0，并不影响最重要的 k 列结果。

我不知道这个证明的书面表示是怎么样的。但是感觉严谨的公式表达也会非常繁琐。

有时候有的问题要“站得高处”思考更容易；有时候有的问题，像此处，反而以特例试验一下，更容易理解。

不过也只是方便理解。“第一个人”是怎么从无到有想到这种组合方式的，还是非常好奇。

至此，PLU 分解就说明完毕了，我们记

\(L=\widetilde{L}_1^{-1}\widetilde{L}_2^{-1}L_3^{-1}\)

得

\(PA=LU\)

2.1 PLU 分解的代码实现

已经知道了具体的步骤，我们来看如何代码实现。如代码清单 2 所示，第 15 至 26 行算出绝对值最大主元。

我们不用“严格按照”上述的公式一次性得到所求的值，可以逐步完成。如 \(P\)，我们可以依次交换得到。同样 \(L\)，我们也可以后续逐步“修正”它的位置。

代码清单 2 PLU 分解
  1. void PLUDecomposition(Matrix* P, Matrix* L, Matrix* U)
  2. {
  3.     assert(ROWS == COLS);
  4.  
  5.     T maxEle, ele;
  6.     size_t i, j, k;
  7.     size_t maxRow;
  8.  
  9.     *U = *this;
  10.     L->ToIdentity();
  11.     P->ToIdentity();
  12.  
  13.     for (k = 0; k < ROWS; k++)
  14.     {
  15.         maxEle = std::abs(U->at(k, k));
  16.         maxRow = k;
  17.  
  18.         for (i = k + 1; i < ROWS; i++)
  19.         {
  20.             ele = std::abs(U->at(i, k));
  21.             if (ele > maxEle)
  22.             {
  23.                 maxEle = ele;
  24.                 maxRow = i;
  25.             }
  26.         }
  27.  
  28.         if (k != maxRow)
  29.         {
  30.             U->SwapRows(k, maxRow);
  31.             P->SwapRows(k, maxRow);
  32.  
  33.             // SwapRowsPartial
  34.             for (j = 0; j < k; j++)
  35.                 std::swap(L->at(k, j), L->at(maxRow, j));
  36.         }
  37.  
  38.         for (i = k + 1; i < ROWS; i++)
  39.         {
  40.             assert(U->at(k, k) != 0);
  41.             L->at(i, k) = U->at(i, k) / U->at(k, k);
  42.             for (size_t j = k; j < ROWS; j++)
  43.                 U->at(i, j) -= L->at(i, k) * U->at(k, j);
  44.         }
  45.     }
  46. }

3. 基于 PLU 分解求逆

现在我们已经可以得到 PLU 分解的结果了，我们来看如何通过它算矩阵的逆。和求解方程组一样，我们求 \(x\)，使得

\(Ax=b\)

两边同乘以 \(P\)，得

\(PAx=Pb\)

代入 LU 分解，得

\(LUx=Pb\)

上下三角求解是方便的，所以我们 \(Y=Ux\)，即求上三角形式

\(L(Ux)=Pb \Leftrightarrow LY=Pb\)

得到 \(Y\) 后，最后再求下三角形式

\(Ux=Y\)

针对求逆的情况，我们只需要令 \(b=E\) 即可。

3.1 代码实现

首先，如代码清单 3.1 所示，我们实现上、下三角矩阵形式的方程求解。

代码清单 3.1 求解方程
  1. // LY=B
  2. template<size_t N>
  3. static void ForwardSubstitute(const Matrix& L, const Matrix<T, ROWS, N>& B, Matrix<T, COLS, N>* Y)
  4. {
  5.     assert(ROWS == COLS);
  6.     T sum;
  7.     for (size_t i = 0; i < COLS; i++)
  8.     {
  9.         for (size_t j = 0; j < N; j++)
  10.         {
  11.             sum = B.at(i, j);
  12.  
  13.             for (size_t k = 0; k < i; k++)
  14.                 sum -= L.at(i, k) * Y->at(k, j);
  15.  
  16.             Y->at(i, j) = sum / L.at(i, i);
  17.         }
  18.     }
  19. }
  20.  
  21. // UX=Y
  22. template<size_t N>
  23. static void BackwardSubstitute(const Matrix& U, const Matrix<T, ROWS, N>& Y, Matrix<T, COLS, N>* X)
  24. {
  25.     assert(ROWS == COLS);
  26.     T sum;
  27.     for (int i = COLS - 1; i >= 0; i--)
  28.     {
  29.         for (size_t j = 0; j < N; j++)
  30.         {
  31.             sum = Y.at(i, j);
  32.             for (size_t k = i + 1; k < COLS; k++)
  33.                 sum -= U.at(i, k) * X->at(k, j);
  34.  
  35.             X->at(i, j) = sum / U.at(i, i);
  36.         }
  37.     }
  38. }

矩阵求逆的过程如代码清单 3.2 所示，和我们上述介绍的过程一样，首先进行 PLU 分解，接着依次求解上、下三角矩阵形式的方程组。

代码清单 3.2 求逆
  1. Matrix Inverse()
  2. {
  3.     assert(ROWS == COLS);
  4.     TMatrix P, L, U;
  5.     PLUDecomposition(&P, &L, &U);
  6.  
  7.     //LY=P
  8.     TMatrix X, Y;
  9.     ForwardSubstitute(L, P, &Y);
  10.     //UX=Y
  11.     BackwardSubstitute(U, Y, &X);
  12.  
  13.     return X;
  14. }

最后我们写一个测试用例，如代码清单 3.3 所示，我们计算矩阵的逆，并进行打印显示。同时，我们也测试一下奇异矩阵的情况，目前我们的处理流程中会发生断言。

本章完整代码在 tag/matrix_ops。

代码清单 3.3 测试用例
  1. TMatrixOpsTask::TMatrixOpsTask()
  2. {
  3.     tmath::Mat4f m = { 8, 7, 9, 5,
  4.                        4, 3, 3, 1,
  5.                        2, 1, 1, 0,
  6.                        6, 7, 9, 8 };
  7.     m.Print("m");
  8.  
  9.     tmath::Mat4f mi = m.Inverse();
  10.     mi.Print("inverse of m");
  11.  
  12. #if 0
  13.     tmath::Mat4f s = { 1, 2, 3, 4,
  14.                        2, 4, 6, 8,
  15.                        1, 1, 1, 1,
  16.                        0, 1, 2, 3 };
  17.     s.Print("singular");
  18.  
  19.     tmath::Mat4f si = s.Inverse();
  20. #endif
  21. }