在我用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系統

提高精度

相關問題