引擎-关于添加矩阵求逆函数的一些事宜

具体实现

自研引擎今天添加了矩阵求逆函数，首先先来讲一下我的矩阵数据结构。
首先我的矩阵是由一个二层嵌套的 StaticArray 构成的。
StaticArray 也是一个自建类型，类似于 std::array，在矩阵中表示为 StaticArray<StaticArray<T, N>, M> Data;
由于这个矩阵利用模板去规定维度，所以必须在编译期对矩阵的维度进行设置
也就是 Math::Matrix<double, 4, 4> 其中4可以自由修改维度。

了解了我的矩阵数据结构的基本原理后，就可以开始着手做矩阵的求逆了。
首先了解一下初等变换，初等行变换是利用一种规则变换矩阵的方式，一般需要用到：

• 任意两行交换
• 任意一行乘以一个任意数
• 任意一行加上另外一行

当初等变换结束后，矩阵的秩与矩阵的行列式的值不会发生变化。
我们需要利用初等行变换将一个矩阵变为 RREF（最简行阶梯形式），如果一个矩阵可逆，那么他的秩必须等于阶数，可以通俗的理解为行数与列数需要相等。
利用初等行变换，先将矩阵化为 RREF。

    inline constexpr Matrix GetRREF() const requires (CArithmetic<T>)
    {
        Matrix Result = *this;
        int32 Lead = 0;
        for (int32 Row = 0; Row < M; Row++)
        {
            if (N <= Lead)
            {
                break;
            }
            int32 i = Row;
            while (MathAlgo::NearlyEqual(Result[i][Lead], T(0)))
            {
                i++;
                if (i == M)
                {
                    i = Row;
                    Lead++;
                    if (N == Lead)
                    {
                        Lead--;
                        break;
                    }
                }
            }

            if (i != Row)
            {
                RangedSwap(&Result[i][0], &Result[i][N], &Result[Row][0]);
            }

            double Lv = Result[Row][Lead];
            for (int32 j = 0; j < N; j++)
            {
                Result[Row][j] /= Lv;
            }

            for (int32 k = 0; k < M; k++)
            {
                if (k != Row)
                {
                    double TempLV = Result[k][Lead];
                    for (int32 j = 0; j < N; j++)
                    {
                        Result[k][j] -= TempLV * Result[Row][j];
                    }
                }
            }
            Lead++;
        }

        return Result;
    }

当矩阵变为 RREF 形式后，我们需要用到求逆的方法。

对秩等于阶数的矩阵（原本的矩阵，非RREF）的右侧新增一个维度相等的矩阵，使这个矩阵变为增广矩阵。
维度将会变化为：Matrix<T, M, N> -> Matrix<T, M, N * 2>;
新增加在右侧的矩阵，必须是一个单位矩阵，也就是

   \begin{bmatrix}
   1 & 0 & 0 \\
   0 & 1 & 0 \\
   0 & 0 & 1 \\
   \end{bmatrix}

的形式，假设我们的矩阵一开始是:

   \begin{bmatrix}
   10 & 0 & 9 \\
   3 & 1 & 2 \\
   0 & 7 & 1 \\
   \end{bmatrix}

那么在增广之后的结果应该为：

   \begin{bmatrix}
   10 & 0 & 9 & 1 & 0 & 0 \\
   3 & 1 & 2 & 0 & 1 & 0\\
   0 & 7 & 1 & 0 & 0 & 1\\
   \end{bmatrix}

接下来，对增广矩阵进行 RREF 运算，使其原矩阵部分变成最简矩阵。
此时增广的额外部分，也就是添加在右侧的单位矩阵，便是矩阵的逆矩阵。

   \begin{bmatrix}
   1 & 0 & 0 & -0.22 & 1.06 & -0.15 \\
   0 & 1 & 0 & -0.05 & 0.16 & 0.11\\
   0 & 0 & 1 & 0.35 & -1.18 & 0.17\\
   \end{bmatrix}

剥离出增广部分，并创建成新的矩阵，与原有矩阵进行乘法运算。
原矩阵 * 逆矩阵

    Matrix<double, 3, 3> A = {
        {10, 0, 9},
        {3, 1, 2},
        {0, 7, 1}
    };
    const Matrix<double, 3, 3> B = A;
    A.Inverse();
    Matrix<double, 3, 3> C = B * A;

如代码所示，C为单位矩阵，验证成功。

弄个实例

众所周知，笛卡尔坐标系由一个原点和三个坐标轴。
假设现在有一个向量：

\begin{bmatrix}
x \\
y \\
z \\
\end{bmatrix}

则同一坐标系下，另外一个向量可以表示为原向量坐标的线性函数。

\begin{bmatrix}
x` \\
y` \\
z` \\
\end{bmatrix}

也就是说在实现原坐标系到新坐标系的线性变换的时候，可以表示为如下形式。

\begin{bmatrix}
x` \\
y` \\
z` \\
\end{bmatrix}
=
\begin{bmatrix}
U_{1} & V_{1} & W_{1} \\
U_{2} & V_{2} & W_{2} \\
U_{3} & V_{3} & W_{3} \\
\end{bmatrix}
\begin{bmatrix}
x \\
y \\
z \\
\end{bmatrix}
+
\begin{bmatrix}
T_{1} \\
T_{2} \\
T_{3} \\
\end{bmatrix}

其中向量T表示原向量到新坐标系原点的平移向量，而 U V W 矩阵则表示了坐标系进行线性变换的时候坐标系的坐标轴方向的变换方式。
若此线性变换可逆，则从新坐标变换回原始坐标系的公式为：

\begin{bmatrix}
x \\
y \\
z \\
\end{bmatrix}
=
\begin{bmatrix}
U_{1} & V_{1} & W_{1} \\
U_{2} & V_{2} & W_{2} \\
U_{3} & V_{3} & W_{3} \\
\end{bmatrix}^{-1}
\begin{pmatrix}
\begin{bmatrix}
x` \\
y` \\
z` \\
\end{bmatrix}
-
\begin{bmatrix}
T_{1} \\
T_{2} \\
T_{3} \\
\end{bmatrix}
\end{pmatrix}

此时利用矩阵的逆可逆转矩阵，使其恢复成原矩阵（假设矩阵可逆）