Ste*_*eve 4 python math scipy matrix-inverse numerical-stability
在我在C中工作的数值解算器中,我需要反转2x2矩阵,然后在右侧乘以另一个矩阵:
C = B . inv(A)
我一直在使用倒置2x2矩阵的以下定义:
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);
在我解算器的前几次迭代中,这似乎给出了正确的答案,然而,经过几个步骤后,事情开始增长并最终爆炸.
现在,与使用SciPy的实现相比,我发现相同的数学不会爆炸.我能找到的唯一区别是SciPy代码使用scipy.linalg.inv(),它在内部使用LAPACK来执行反转.
当我用inv()上面的计算替换调用时,Python版本确实会爆炸,所以我很确定这是问题所在.计算中的微小差异正在蔓延,这使我相信它是一个数值问题 - 对于反演操作并不完全令人惊讶.
我正在使用双精度浮点数(64位),希望数值问题不会成为问题,但显然情况并非如此.
但是:我想在我的C代码中解决这个问题,而不需要调用像LAPACK这样的库,因为将它移植到纯C的全部原因是让它在目标系统上运行.此外,我想了解这个问题,而不仅仅是呼唤黑匣子.最终,如果可能的话,我也希望它以单精度运行.
所以,我的问题是,对于这样一个小矩阵,有一个数值上更稳定的方法来计算A的倒数吗？
谢谢.
编辑:目前试图弄清楚我是否可以通过解决来避免反演C.
计算行列式并不稳定。更好的方法是使用带有部分旋转的 Gauss-Jordan，您可以在此处轻松明确地计算出来。
让我们求解该系统（使用 c, f = 1, 0，然后使用 c, f = 0, 1 求逆）
a * x + b * y = c
d * x + e * y = f
在伪代码中，这读作
if a == 0 and d == 0 then "singular"
if abs(a) >= abs(d):
alpha <- d / a
beta <- e - b * alpha
if beta == 0 then "singular"
gamma <- f - c * alpha
y <- gamma / beta
x <- (c - b * y) / a
else
swap((a, b, c), (d, e, f))
restart
这比行列式 + comatrix 更稳定（beta行列式 * 是某个常数，以稳定的方式计算）。您可以计算出完整的旋转等效值（即可能交换 x 和 y，以便第一个除法a是aa、b、d、e 中幅度最大的数字），并且在某些情况下这可能更稳定，但上述方法对我来说效果很好。
这相当于执行 LU 分解（如果要存储此 LU 分解，请存储 gamma、beta、a、b、c）。
计算 QR 分解也可以显式完成（并且如果正确执行，也非常稳定），但速度较慢（并且涉及求平方根）。这是你的选择。
如果您需要更好的精度（上述方法是稳定的，但存在一些舍入误差，与特征值的比率成正比），您可以“求解校正”。
事实上，假设您用上述方法解决A * x = b了x。现在你计算A * x，你发现它并不完全等于b，有一个小错误：
A * x - b = db
现在，如果你求解dxin A * dx = db，你有
A * (x - dx) = b + db - db - ddb = b - ddb
其中ddb是 的数值求解引起的误差A * dx = db，通常远小于db（因为db远小于b）。
您可以重复上述过程，但通常需要一步才能恢复完整的机器精度。
不要颠倒矩阵.几乎总是,在不反转矩阵的情况下,您可以更快,更准确地完成使用逆转完成的事情.矩阵求逆本质上是不稳定的,将其与浮点数混合会产生麻烦.
说法C = B . inv(A)是一样的话,你想解决AC = B的C.您可以通过拆分每个做到这一点B,并C为两列.解决A C1 = B1和A C2 = B2会产生C.
你的代码很好; 但是,它可能会导致四次减法中任何一次的精度损失.
考虑使用更高级的技术,例如matfunc.py中使用的技术.该代码使用通过Householder反射实现的QR分解来执行反演.通过迭代细化进一步改善了结果.