在我用C工作的数值解算器的数值稳定逆，我需要反转一个2x2矩阵，它随后被另一矩阵相乘右侧：一个2x2矩阵

C = B . inv(A)

我有在使用倒置×2矩阵的定义如下：

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的逆？

谢谢。

编辑：目前试图找出我是否可以只avoid the inversion通过解决C。

来源

2011-12-31 Steve

它在什么意义上爆炸？数字溢出？矩阵元素是什么类型，有哪些值？ – fge

回答fge的问题会有所帮助。这里也有可能被零除，这可能会导致你的“爆炸”。 – greg

对不起，如果我不清楚，爆炸不是直接来自这个操作，而是从这个操作引入的错误到一个反馈函数，我没有在这里描述。 – Steve

请勿颠倒矩阵。几乎总是，你使用逆来完成的事情可以更快，更准确地完成，而不需要反转矩阵。矩阵求逆本质上是不稳定的，并且将其与浮点数混合在一起会造成麻烦。

说C = B . inv(A)是等于说你要解决AC = B为C. 您可以通过每个B和C分裂成两列做到这一点。解决A C1 = B1和A C2 = B2将产生C.

来源

2011-12-31 22:25:29 redxaxder

不错的想法，谢谢我会尝试。 – Steve

似乎已经做到了，谢谢！虽然我会注意到，我很担心，因为我发现解决方案仍然包含一个减法除法，就像在逆的公式中一样，如Raymond H.提到的那样可能是不好的。但是，它似乎不会导致相同的分歧，所以我认为它不太可能引入错误。这是有道理的，因为我不需要乘以一个非常小的倒数，而是直接转到最终的解决方案'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

你的代码没问题;但是，它会从四个减法中的任何一个中减去loss of precision。

考虑使用更高级的技术，如matfunc.py中使用的技术。该代码使用使用Householder reflections实现的QR decomposition执行反演。使用iterative refinement可以进一步提高结果。

来源

2011-12-31 20:50:05

是的，我想这是由于一小部分的划分，但减法可能是更糟的元素在这里！我会检查这些代码，非常感谢你。这可能需要一点时间才能转化为C. – Steve

使用Jacobi方法，这是一种迭代方法，它涉及到只反演A的主对角线，这非常简单并且不易反转整个矩阵的数值不稳定性。

来源

2012-01-02 20:15:49

计算行列式不稳定。更好的方法是使用高斯 - 乔丹的部分透视，这样你就可以在这里明确地计算出来。

求解一个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的第一个除法是a是a，b，d，e中最大的数量级），并且这可能会更稳定有些情况下，但上述方法对我来说工作得很好。

这相当于执行LU分解（如果要存储此LU分解，请存储gamma，beta，a，b，c）。

计算QR分解也可以明确地进行（并且如果你正确地做它也是非常稳定的），但是它更慢（并涉及取平方根）。这是你的选择。

提高精度

如果需要更好的精确度（在上述方法是稳定的，但有一些舍入误差，正比于特征值的比），可以在“解决对于校正”。

事实上，假设你用上述方法解决了A * x = b的x。现在计算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）。

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

来源

2012-04-17 09:09:45

我同意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-13 01:07:40 moldovean

一个2x2矩阵

回答

求解一个2x2系统

提高精度

相关问题