Numpy矩陣行列式精度問題

Question

我想在python上編寫一個腳本來確定使用高斯方法的矩陣行列式。它工作正常，但精度對我來說還不夠。 我的代碼是：

import scipy.linalg as sla 
import numpy as np 
def my_det(X): 
    n = len(X) 
    s = 0 
    if n != len(X[0]): 
     return ValueError 
    for i in range(0, n): 
     maxElement = abs(X[i][i]) 
     maxRow = i 
     for k in range(i+1, n): 
      if abs(X[k][i]) > maxElement: 
       maxElement = abs(X[k][i]) 
       maxRow = k 
     if maxRow != i: 
      s += 1 
     for k in range(i, n): 
      X[i][k], X[maxRow][k] = X[maxRow][k], X[i][k] 
     for k in range(i+1, n): 
      c = -X[k][i]/X[i][i] 
      for j in range(i, n): 
       if i == j: 
        X[k][j] = 0 
       else: 
        X[k][j] += c * X[i][j] 
    det = (-1)**s 
    for i in range(n): 
     det *= X[i][i] 
    return det 

而且我對這個代碼測試：

for x in range(10): 
X = np.random.rand(3,3) 
if np.abs(my_det(X) - sla.det(X)) > 1e-6: 
    print('FAILED') 

我的函數失敗的所有測試。我嘗試了小數，但沒有幫助。 有什麼問題？

Answer 1

爲什麼不使用numpy.linalg.det或scipy.linalg.det函數？ 這些函數使用LU分解和LAPACK計算行列式。它會比任何'手動'功能都快。

Answer 2

爲什麼代碼失敗的測試條件，abs(my_det(X) - sla.det(X)) < 1e-6究其原因，是不是因爲缺乏精確度，而在於符號的變化 帶來的my_det意想不到的副作用變異X：

X[i][k], X[maxRow][k] = X[maxRow][k], X[i][k] 

這行交換改變了行列式的符號。 代碼使用s來調整符號的變化，但X本身以改變行列式符號的方式更改爲 。

因此，傳遞給my_det的X與隨後傳遞給sla.det的X不一樣。這裏是哪裏的X的變化改變了行列式的符號的例子：

In [55]: X = np.random.rand(3, 3); X 
Out[55]: 
array([[ 0.38062719, 0.41892961, 0.88277747], 
     [ 0.39881724, 0.00188804, 0.79258322], 
     [ 0.40195279, 0.3950311 , 0.32771527]]) 

In [56]: my_det(X) 
Out[56]: 0.098180005266934267 

In [57]: X 
Out[57]: 
array([[ 0.40195279, 0.3950311 , 0.32771527], 
     [ 0.  , -0.39006151, 0.46742438], 
     [ 0.  , 0.  , 0.62620267]]) 

In [58]: sla.det(X) 
Out[58]: -0.09818000526693427 

---

您可以通過的X副本修復內部my_det問題：

def my_det(X): 
    X = np.array(X, copy=True) # copy=True is the default; shown here for emphasis 
    ... 

因此，後續在my_det之內更改爲X不再影響X之外的 my_det。

---

import scipy.linalg as sla 
import numpy as np 


def my_det(X): 
    X = np.array(X, dtype='float64', copy=True) 
    n = len(X) 
    s = 0 
    if n != len(X[0]): 
     return ValueError 
    for i in range(0, n): 
     maxElement = abs(X[i, i]) 
     maxRow = i 
     for k in range(i + 1, n): 
      if abs(X[k, i]) > maxElement: 
       maxElement = abs(X[k, i]) 
       maxRow = k 
     if maxRow != i: 
      s += 1 
     for k in range(i, n): 
      X[i, k], X[maxRow, k] = X[maxRow, k], X[i, k] 
     for k in range(i + 1, n): 
      c = -X[k, i]/X[i, i] 
      for j in range(i, n): 
       if i == j: 
        X[k, j] = 0 
       else: 
        X[k, j] += c * X[i, j] 
    det = (-1)**s 
    for i in range(n): 
     det *= X[i, i] 
    return det 


for i in range(10): 
    X = np.random.rand(3, 3) 
    diff = abs(my_det(X) - sla.det(X)) 
    if diff > 1e-6: 
     print('{} FAILED: {:0.8f}'.format(i, diff)) 

---

還要注意的是D型細胞事項：

In [88]: my_det(np.arange(9).reshape(3,3)) 
Out[88]: 6 

而正確答案是

In [89]: my_det(np.arange(9).reshape(3,3).astype(float)) 
Out[89]: 0.0 

由於my_det使用部門（在c = -X[k, i]/X[i, i]），我們需要X具有浮點數dtype，因此/執行浮點除法，而不是整數除法。 因此，要解決這個問題，使用X = np.asarray(X, dtype='float64')確保X具有D型float64：

def my_det(X): 
    X = np.array(X, dtype='float64', copy=True) 
    ... 

隨着這一變化，

In [91]: my_det(np.arange(9).reshape(3,3)) 
Out[91]: 0.0 

現在給出正確的答案。