ナンシー行列の行列式の精度の問題

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') 

My機能は、すべてのテストに失敗します。私はDecimalsを試しましたが、それは助けになりませんでした。 どうしたの？

Answer 1

なぜnumpy.linalg.detまたはscipy.linalg.det機能を使用しないのですか？ これらの関数は、LU分解とLAPACKを使用して行列式を計算します。それはどんな「手動」機能よりも速くなります。

Answer 2

コードは、テスト条件を失敗した理由、abs(my_det(X) - sla.det(X)) < 1e-6は、精度の不足によるものではなく、符号の変化 はXを変異my_detの意図しない副作用をもたらし：

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 

---

あなたはmy_det内Xのコピーを作成することにより、問題を解決することができます

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

したがって、それに続きますX内のmy_detへの変更は、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)) 

---

また

注意してDTYPE事項：

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

my_detは（c = -X[k, i]/X[i, i]で）除算を使用しているので、正しい答えは

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

ている間は、私たちはを必要とします0は浮動小数点dtypeを持つため、/は整数除算ではなく浮動小数点除算を実行します。

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

は今、正しい答えを与え、この変更により

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

：これ、XはDTYPE float64を持っていることを確認するためにX = np.asarray(X, dtype='float64')を使用し修正する 。