python - 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

代码未通过测试条件的原因，abs(my_det(X) -
sla.det(X)) < 1e-6 ，不是由于缺乏精度，而是由于变化in 符号带来了 my_det 的意外副作用突变X :

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

此行交换会更改行列式的符号。代码使用s调整符号的变化，但是 X本身被改变以改变行列式符号的方式。

所以X传递至my_det与 X 不一样随后传递给 sla.det 。这是一个示例，其中更改 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
    ...

因此，后续更改为 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))

<小时/>

另请注意，数据类型很重要:

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具有 float 据类型，以便 /执行浮点除法，而不是整数除法。因此，要解决此问题，请使用 X = np.asarray(X, dtype='float64')确保X有数据类型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

现在给出正确答案。