2x2矩阵的数值稳定逆 [英] numerically stable inverse of a 2x2 matrix
问题描述
在我使用C进行的数值求解器中，我需要求一个2x2矩阵的倒数，然后在右边将其乘以另一个矩阵:
In a numerical solver I am working on in C, I need to invert a 2x2 matrix and it then gets multiplied on the right side by another matrix:
C = B . inv(A)
我一直在使用以下2x2逆矩阵的定义:
I have been using the following definition of an inverted 2x2 matrix:
a = A[0][0];
b = A[0][1];
c = A[1][0];
d = A[1][1];
invA[0][0] = d/(a*d-b*c);
invA[0][1] = -b/(a*d-b*c);
invA[1][0] = -c/(a*d-b*c);
invA[1][1] = a/(a*d-b*c);
在我的求解器的前几次迭代中，这似乎给出了正确的答案，但是，经过几步，事情开始发展并最终爆炸.
In the first few iterations of my solver this seems to give the correct answers, however, after a few steps things start to grow and eventually explode.
现在，与使用SciPy的实现相比，我发现相同的数学不会爆炸.我能找到的唯一区别是SciPy代码使用scipy.linalg.inv()，内部使用LAPACK来执行反转.
Now, comparing to an implementation using SciPy, I found that the same math does not explode. The only difference I can find is that the SciPy code uses scipy.linalg.inv(), which internally uses LAPACK internally to perform the inversion.
当我用上述计算替换对inv()的调用时，Python版本确实爆炸了，因此，我很确定这是问题所在.计算中的细微差异正在悄悄蔓延，这使我相信这是一个数值问题，对于反演操作来说并不完全令人惊讶.
When I replace the call to inv() with the above calculations the Python version does explode, so I'm pretty sure this is the problem. Small differences in the calculations are creeping in, which leads me to believe it is a numerical problem--not entirely surprising for an inversion operation.
我正在使用双精度浮点数(64位)，希望数值问题不会成为问题，但显然不是这种情况.
I am using double-precision floats (64-bit), hoping that numerical issues would not be a problem, but apparently that is not the case.
但是:我想在我的C代码中解决此问题而无需调用LAPACK之类的库，因为将其移植到纯C的全部原因是使它在目标系统上运行.而且，我想了解问题，而不仅仅是打电话给黑匣子.最终，如果可能的话，我也希望它也能以单精度运行.
But: I would like to solve this in my C code without needing to call out to a library like LAPACK, because the whole reason for porting it to pure C is to get it running on a target system. Moreover, I'd like to understand the problem, not just call out to a black box. Eventually I'd like to it run with single-precision too, if possible.
所以，我的问题是，对于这么小的矩阵，是否存在一种数值上更稳定的方法来计算A的逆?
So, my question is, for such a small matrix, is there a numerically more stable way to calculate the inverse of A?
谢谢.
Currently trying to figure out if I can just avoid the inversion by solving for C.
推荐答案
不要反转矩阵.几乎总是，您使用逆函数完成的事情可以更快，更准确地完成，而无需反转矩阵.矩阵求逆本质上是不稳定的，将其与浮点数混合会带来麻烦.
Don't invert the matrix. Almost always, the thing you're using the inverse to accomplish can be done faster and more accurately without inverting the matrix. Matrix inversion is inherently unstable, and mixing that with floating point numbers is asking for trouble.
说C = B . inv(A)与说要解决C的AC = B相同.
您可以通过将每个B和C分成两列来完成此操作.解决A C1 = B1和A C2 = B2会产生C.
Saying C = B . inv(A) is the same as saying you want to solve AC = B for C.
You can accomplish this by splitting up each B and C into two columns. Solving A C1 = B1 and A C2 = B2 will produce C.
这篇关于2x2矩阵的数值稳定逆的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
