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

内容目录

具体实现

自研引擎今天添加了矩阵求逆函数,首先先来讲一下我的矩阵数据结构。
首先我的矩阵是由一个二层嵌套的 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}

此时利用矩阵的逆可逆转矩阵,使其恢复成原矩阵(假设矩阵可逆)

同步Sch

Recent Posts

UE-Rider 无法解析 SDK 解决方案

无法解析 SDK“Micros…

5 天 ago

杂项 – 关于Ascii字符画生成相关网站

网站 如果想实现这样的注释效果…

2 周 ago

C++-分数泛型模板类实现

前言 在引擎库中,有时候可能会…

3 周 ago

C++-HashMap 数据结构的实现

前要说明 HashMap 是一…

1 月 ago

杂项 – Markdown

跳转到 Markdown – …

1 月 ago

杂谈-我都玩过什么游戏

在面试的过程中,我发现游戏公司…

3 月 ago