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 

现在给出正确的答案。