一个稳定的二维矩阵的数值逆矩阵

Question

在我正在编写的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来避免求逆。

- Steve

---

它以什么意义上爆炸？数值溢出吗？矩阵元素是什么类型的，其中有什么值？ - fge

回答fge的问题会有所帮助。此外，这里还存在除以零的潜在风险，可能会导致您的“爆炸”。 - greg

对不起如果我没有表达清楚，这次操作并不直接导致崩溃，而是由此操作引入到一个我没有在这里描述的反馈函数中的错误所导致的。 - Steve

Answer 1

计算行列式是不稳定的。更好的方法是使用带有部分主元选取的高斯-约旦法，您可以在此处轻松地明确计算出来。

解决2x2系统

让我们解决这个系统（使用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

这个方法比行列式+余子式更加稳定（beta是行列式*某个常数，在计算中以一种稳定的方式计算）。您可以计算全主元的等价形式（即可能交换x和y，使得第一个被a除的数在a、b、d、e中是数量级最大的），在某些情况下，这可能更加稳定，但上述方法对我来说已经很有效了。

这相当于进行LU分解（如果您想存储此LU分解，则存储gamma、beta、a、b、c）。

计算QR分解也可以明确地进行（并且只要您正确地执行就非常稳定），但速度较慢（并涉及取平方根）。选择权在您手中。

提高精度

如果您需要更好的精度（上述方法稳定，但存在一些舍入误差，与特征值的比率成比例），您可以“求解修正”。

实际上，假设您使用上述方法为x解决了A * x = b。现在，您计算A * x，并发现它并不完全等于b，有一些轻微的误差：

A * x - b = db

现在，如果你解决 A * dx = db 中的 dx，你就会得到

A * (x - dx) = b + db - db - ddb = b - ddb

其中ddb是通过对A * dx = db进行数值求解产生的误差，通常比db要小得多（因为db比b小得多）。

您可以迭代上述过程，但通常需要一步来恢复完整的机器精度。

- Alexandre C.

---

谢谢您，这太棒了！快点提醒一下：对我的应用程序来说，“单数”情况很重要，而检查beta == 0一旦值超过 ~30位精度就可能得到错误的答案。对我来说，检查a * e == b * d的结果更好。 - Nemo

Answer 2

您的代码很好，但是它存在风险精度损失，这可能来自其中任意一个减法。

考虑使用更先进的技术，例如matfunc.py中使用的技术。该代码使用QR分解实现反转，使用Householder反射进行实现。接着，使用迭代细化对结果进行进一步改善。

- Raymond Hettinger

---

是的，我想这是由于除以一个小数导致的，但是这里的减法可能是更糟糕的因素！我会查看那个代码，非常感谢。将其转换为C语言可能需要一点时间。 - Steve

Answer 3

不要求逆矩阵。几乎总是，你使用逆矩阵所要达到的目的可以更快且更准确地完成，而不需要求逆矩阵。矩阵求逆本身就不稳定，并将其与浮点数混合使用会带来问题。

说 C = B . inv(A) 和说你想解决 AC = B 得到 C 是一样的。 你可以通过将每个 B 和 C 分成两列来实现这一点。解决 A C1 = B1 和 A C2 = B2 将产生 C。

- redxaxder

---

似乎解决了问题，谢谢！系统看起来更加稳定了一些，尽管我注意到解决方案仍然包含一个除法减法的操作，就像求逆公式中所提到的那样，这可能是不好的。然而，它似乎没有引起相同的发散，所以我想它不太容易引入错误。这是有道理的，因为我不需要乘以一个非常小的倒数，而是直接得到最终解C。 - Steve

我将在此为后人发布解决方案： c00 = -(a11*b00-a10*b01) / (a01*a10-a00*a11); c01 = (a01*b00-a00*b01) / (a01*a10-a00*a11); c10 = (a10*b11-a11*b10) / (a01*a10-a00*a11); c11 = -(a00*b11-a01*b10) / (a01*a10-a00*a11); c20 = (a10*b21-a11*b20) / (a01*a10-a00*a11); c21 = -(a00*b21-a01*b20) / (a01*a10-a00*a11); - Steve

哎呀，发布解决方案有更好的方法吗？评论不允许代码块？也许我应该将其添加到redxaxder的答案中？ - Steve

Answer 4

使用雅可比方法，这是一种迭代方法，涉及到仅反转A的主对角线，这非常简单，并且比反转整个矩阵更不容易出现数值不稳定性。

- Jean-Victor Côté

Answer 5

我同意Jean-Victor的看法，你应该使用雅可比方法。这是我的示例：

#Helper functions:
def check_zeros(A,I,row, col=0):
"""
returns recursively the next non zero matrix row A[i]
"""
if A[row, col] != 0:
    return row
else:
    if row+1 == len(A):
        return "The Determinant is Zero"
    return check_zeros(A,I,row+1, col)

def swap_rows(M,I,row,index):
"""
swaps two rows in a matrix
"""
swap = M[row].copy()
M[row], M[index] = M[index], swap
swap = I[row].copy()
I[row], I[index] = I[index], swap

# Your Matrix M
M = np.array([[0,1,5,2],[0,4,9,23],[5,4,3,5],[2,3,1,5]], dtype=float)
I = np.identity(len(M))

M_copy = M.copy()
rows = len(M)

for i in range(rows):
index =check_zeros(M,I,i,i)
while index>i:
    swap_rows(M, I, i, index)
    print "swaped"
    index =check_zeros(M,I,i,i) 

I[i]=I[i]/M[i,i]
M[i]=M[i]/M[i,i]   

for j in range(rows):
    if j !=i:
        I[j] = I[j] - I[i]*M[j,i]
        M[j] = M[j] - M[i]*M[j,i]
print M
print I  #The Inverse Matrix

- moldovean