5

在我正在使用 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。

4

5 回答 5

6

不要反转矩阵。几乎总是，您使用逆来完成的事情可以更快，更准确地完成，而无需反转矩阵。矩阵求逆本质上是不稳定的，将其与浮点数混合是自找麻烦。

说C = B . inv(A)与说你想解决AC = BC 是一样的。你可以通过将每一列分成两列来完成此B操作C。求解A C1 = B1并将A C2 = B2产生 C.

于 2011-12-31T22:25:29.720 回答
5

你的代码很好；但是，它可能会从四个减法中的任何一个中损失精度。

考虑使用更高级的技术，例如matfunc.py中使用的技术。该代码使用通过Householder 反射实现的QR 分解来执行反转。结果通过迭代细化得到进一步改进。

于 2011-12-31T20:50:05.020 回答
4

计算行列式是不稳定的。更好的方法是使用带有部分旋转的 Gauss-Jordan，您可以在此处轻松地明确计算。

求解 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

这比行列式+ 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）。

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

于 2012-04-17T09:09:45.933 回答
1

使用 Jacobi 方法，这是一种迭代方法，只涉及“反转” A 的主对角线，与反转整个矩阵相比，它非常简单且不易出现数值不稳定。

于 2012-01-02T20:15:49.027 回答
1

我同意 Jean-Vicotr 的观点，您可能应该使用 Jacobbian 方法。这是我的例子：

#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
于 2014-09-13T01:07:40.240 回答