高效的4x4矩阵求逆（仿射变换）

Question

我希望有人能指出一种高效的4x4仿射矩阵变换公式。目前我的代码使用余子式展开，并为每个余子式分配一个临时数组。虽然易于阅读，但速度比应该慢。

请注意，这不是作业，我知道如何使用4x4余子式展开手动计算，但这很麻烦，对我来说并不是一个有趣的问题。此外，我已经谷歌过了，找到了一些已经给出公式的网站 (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm)。然而，这个公式可能可以通过预先计算某些乘积进行进一步优化。我相信在某个时间点上，有人会找到适用于此的“最佳”公式吧？

- Budric

---

为什么不使用一些现有的库呢？很可能这些库已经被优化过了。 - kennytm

没错。不幸的是，那段矩阵代码是用Java编写的，然后由GWT编译。大多数库都无法使用。而且它是一个相当狭窄的应用程序。我只需要处理4x4矩阵。我不想为了获得inverse()和multiply()功能而链接一个庞大的线性代数库。 - Budric

@kennytm 如果我们想学习怎么办？ - Justin Meiners

Answer 1

你应该能够利用矩阵是仿射的这一事实，以加速全反演。也就是说，如果你的矩阵看起来像这样

A = [ M   b  ]
    [ 0   1  ]

假设A为4x4的矩阵，M为3x3的矩阵，b为3x1的列向量，并且A的底部一行为(0,0,0,1)，那么

inv(A) = [ inv(M)   -inv(M) * b ]
         [   0            1     ]

根据你的情况，计算 inv(A)*x 的结果可能比实际形成 inv(A) 更快。在这种情况下，事情会变得简单：

inv(A) * [x] = [ inv(M) * (x - b) ]
         [1] = [        1         ] 

其中x是一个3x1的向量（通常是一个三维点）。

最后，如果M表示旋转（即其列是正交的），则可以利用inv(M) = transpose(M)这个事实。然后计算A的逆只需要减去平移分量，并乘以3x3部分的转置。

请注意，矩阵是否正交是您应该从问题分析中知道的事情。在运行时检查它将是相当昂贵的，尽管您可能希望在调试版本中执行此操作以检查您的假设是否成立。

希望所有这些都很清楚...

- celion

---

你从“分块求逆”（http://en.wikipedia.org/wiki/Invertible_matrix）得到的第一个公式是什么？还是有其他的心理技巧？我不太确定inv(A) * x = inv(M) * (x - b)。首先，它们的大小不同——你是在左边删除一行A还是在右边添加一行？其次，我不确定这个方程是怎么来的。第三，我不知道你在这个方程中要解决什么问题。Oliver一直提到不要符号计算逆矩阵，但我不知道这是什么意思——我需要逆矩阵才能进行逆变换。如果你有时间，我想知道。 - Budric

我修改了inv(A) * x公式，以使维度更清晰。第一个公式来自http://en.wikipedia.org/wiki/Affine_transformation。暂且不谈仿射变换，在一般情况下，当您解决A*x = b时，您需要解出inv(A)*b。但是通常情况下，您不需要实际形成inv(A)，只需计算出乘积即可。回到仿射变换，对于3D应用程序，您可能实际上并不需要矩阵的逆，而只是想要逆变换作用于(相乘)向量。如果是这种情况，则使用该公式可能会更快。 - celion

即使您需要存储矩阵的逆，也可以利用它是仿射的事实来减少计算逆的工作量，因为您只需要反转一个3x3矩阵而不是4x4矩阵。如果您知道它是旋转，则计算转置比计算逆快得多，在这种情况下，它们是等效的。 - celion

这里有一个更好的解释，关于我所说的计算inv(A) * x的含义：http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/（作者是SO上的常客）。 - celion

现在清除。非常有趣。 - Budric

Answer 2

在跟进pkhaler和Robin Hilliard上面的精彩回答后，这是Robin的ActionScript 3代码转换成C#方法。希望这能为其他C#开发人员以及需要4x4矩阵求逆函数的C/C++和Java开发人员节省一些打字时间：

public static double[,] GetInverse(double[,] a)
{
    var s0 = a[0, 0] * a[1, 1] - a[1, 0] * a[0, 1];
    var s1 = a[0, 0] * a[1, 2] - a[1, 0] * a[0, 2];
    var s2 = a[0, 0] * a[1, 3] - a[1, 0] * a[0, 3];
    var s3 = a[0, 1] * a[1, 2] - a[1, 1] * a[0, 2];
    var s4 = a[0, 1] * a[1, 3] - a[1, 1] * a[0, 3];
    var s5 = a[0, 2] * a[1, 3] - a[1, 2] * a[0, 3];

    var c5 = a[2, 2] * a[3, 3] - a[3, 2] * a[2, 3];
    var c4 = a[2, 1] * a[3, 3] - a[3, 1] * a[2, 3];
    var c3 = a[2, 1] * a[3, 2] - a[3, 1] * a[2, 2];
    var c2 = a[2, 0] * a[3, 3] - a[3, 0] * a[2, 3];
    var c1 = a[2, 0] * a[3, 2] - a[3, 0] * a[2, 2];
    var c0 = a[2, 0] * a[3, 1] - a[3, 0] * a[2, 1];

    // Should check for 0 determinant
    var invdet = 1.0 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    var b = new double[4, 4];

    b[0, 0] = ( a[1, 1] * c5 - a[1, 2] * c4 + a[1, 3] * c3) * invdet;
    b[0, 1] = (-a[0, 1] * c5 + a[0, 2] * c4 - a[0, 3] * c3) * invdet;
    b[0, 2] = ( a[3, 1] * s5 - a[3, 2] * s4 + a[3, 3] * s3) * invdet;
    b[0, 3] = (-a[2, 1] * s5 + a[2, 2] * s4 - a[2, 3] * s3) * invdet;

    b[1, 0] = (-a[1, 0] * c5 + a[1, 2] * c2 - a[1, 3] * c1) * invdet;
    b[1, 1] = ( a[0, 0] * c5 - a[0, 2] * c2 + a[0, 3] * c1) * invdet;
    b[1, 2] = (-a[3, 0] * s5 + a[3, 2] * s2 - a[3, 3] * s1) * invdet;
    b[1, 3] = ( a[2, 0] * s5 - a[2, 2] * s2 + a[2, 3] * s1) * invdet;

    b[2, 0] = ( a[1, 0] * c4 - a[1, 1] * c2 + a[1, 3] * c0) * invdet;
    b[2, 1] = (-a[0, 0] * c4 + a[0, 1] * c2 - a[0, 3] * c0) * invdet;
    b[2, 2] = ( a[3, 0] * s4 - a[3, 1] * s2 + a[3, 3] * s0) * invdet;
    b[2, 3] = (-a[2, 0] * s4 + a[2, 1] * s2 - a[2, 3] * s0) * invdet;

    b[3, 0] = (-a[1, 0] * c3 + a[1, 1] * c1 - a[1, 2] * c0) * invdet;
    b[3, 1] = ( a[0, 0] * c3 - a[0, 1] * c1 + a[0, 2] * c0) * invdet;
    b[3, 2] = (-a[3, 0] * s3 + a[3, 1] * s1 - a[3, 2] * s0) * invdet;
    b[3, 3] = ( a[2, 0] * s3 - a[2, 1] * s1 + a[2, 2] * s0) * invdet;

    return b;
}

- Anders Gustafsson

Answer 3

万一有人想要省点打字，这是我写的一个AS3版本，基于phkahler在上面发布的链接的第9页（Laplace Expansion Theorem的更有效率的版本）：

public function invert() : Matrix4 {
    var m : Matrix4 = new Matrix4();

    var s0 : Number = i00 * i11 - i10 * i01;
    var s1 : Number = i00 * i12 - i10 * i02;
    var s2 : Number = i00 * i13 - i10 * i03;
    var s3 : Number = i01 * i12 - i11 * i02;
    var s4 : Number = i01 * i13 - i11 * i03;
    var s5 : Number = i02 * i13 - i12 * i03;

    var c5 : Number = i22 * i33 - i32 * i23;
    var c4 : Number = i21 * i33 - i31 * i23;
    var c3 : Number = i21 * i32 - i31 * i22;
    var c2 : Number = i20 * i33 - i30 * i23;
    var c1 : Number = i20 * i32 - i30 * i22;
    var c0 : Number = i20 * i31 - i30 * i21;

    // Should check for 0 determinant

    var invdet : Number = 1 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    m.i00 = (i11 * c5 - i12 * c4 + i13 * c3) * invdet;
    m.i01 = (-i01 * c5 + i02 * c4 - i03 * c3) * invdet;
    m.i02 = (i31 * s5 - i32 * s4 + i33 * s3) * invdet;
    m.i03 = (-i21 * s5 + i22 * s4 - i23 * s3) * invdet;

    m.i10 = (-i10 * c5 + i12 * c2 - i13 * c1) * invdet;
    m.i11 = (i00 * c5 - i02 * c2 + i03 * c1) * invdet;
    m.i12 = (-i30 * s5 + i32 * s2 - i33 * s1) * invdet;
    m.i13 = (i20 * s5 - i22 * s2 + i23 * s1) * invdet;

    m.i20 = (i10 * c4 - i11 * c2 + i13 * c0) * invdet;
    m.i21 = (-i00 * c4 + i01 * c2 - i03 * c0) * invdet;
    m.i22 = (i30 * s4 - i31 * s2 + i33 * s0) * invdet;
    m.i23 = (-i20 * s4 + i21 * s2 - i23 * s0) * invdet;

    m.i30 = (-i10 * c3 + i11 * c1 - i12 * c0) * invdet;
    m.i31 = (i00 * c3 - i01 * c1 + i02 * c0) * invdet;
    m.i32 = (-i30 * s3 + i31 * s1 - i32 * s0) * invdet;
    m.i33 = (i20 * s3 - i21 * s1 + i22 * s0) * invdet;

    return m;
}

通过使用从该方法返回的逆矩阵，我成功地将各种3D变换矩阵乘以单位矩阵。我相信您可以搜索/替换以得到您所需语言的结果。

- Robin Hilliard

---

非常感谢发布，@Robin，这对我的C#项目帮助很大。我在上面的代码中发现了一个小拼写错误：在定义c5时应该更正为i31 * i23。在修复这个问题之后，矩阵反转对我来说运行得非常好。 - Anders Gustafsson

嗨@AndersGustafsson，我想你是指c4的定义-感谢更正-Robin将修复原始内容。 - Johnus

@Johnus，你说得完全正确，当我在评论一个打字错误时犯了一个打字错误，真是太傻了 :-) 谢谢你指出这一点。 - Anders Gustafsson

Answer 4

我记得你可以通过预计算一堆（12个？）2x2行列式来大大缩小代码和时间。将矩阵在垂直方向上分成两半，在上半部分和下半部分中计算每个2x2。这些较小的行列式中的一个将用于您需要进行的更大计算中的每个术语，并且它们每个都会被重复使用。

此外，不要使用单独的行列式函数-重复使用您为伴随矩阵计算的子行列式以获取行列式。

哦，刚刚找到this.。

还有一些改进措施，因为它是某种类型的变换。

- phkahler

---

谢谢，这让我节省了很多时间！ - Budric

非常快，解释得很好。看起来可以工作（还没有对完整的回归测试运行）。再次感谢。 - Budric

+1 对于这个链接；然而，我认为以符号方式计算这些逆是一个错误... 你必须意识到你正在执行多少不必要的乘法/加法。只要这部分代码不是瓶颈，可能还可以接受。 - Olivier Verdier

我在很多事情上都使用4x4，所以我更喜欢广义逆。就像我说的，你可以用特定类型的变换做得更好。链接的论文仍然有用，可以用于执行提问者似乎正在使用的3x3逆运算。如果你知道3x3是一个纯旋转，那么你甚至可以做得更好——我记得它的逆矩阵是转置矩阵。 - phkahler

Answer 5

我认为计算逆矩阵的唯一方法是解决n次方程：Ax = y，其中y跨越单位向量，即第一个是（1,0,0,0），第二个是（0,1,0,0），以此类推。

（使用余子式（克莱姆法则）是不好的选择，除非您想要逆矩阵的符号公式。）

大多数线性代数库都可以让您解决这些线性系统，甚至计算逆矩阵。在Python中的示例（使用numpy）：

from numpy.linalg import inv
inv(A) # here you go

- Olivier Verdier