我试图在python上编写一个脚本来确定使用Gauss方法的矩阵行列式。它工作正常,但对我来说精度还不够。我的代码是:
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')
我的功能未通过所有测试。我尝试过Decimals,但它没有帮助。怎么了?
代码未通过测试条件的原因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
本身的改变方式改变了行列式的符号。
因此传递给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))
另请注意dtype很重要:
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
有dtype 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
现在给出了正确的答案。
为什么不使用numpy.linalg.det
或scipy.linalg.det
函数?这些函数使用LU分解和LAPACK计算行列式。它会比任何“手动”功能更快。
即使对于3x3矩阵或更多,使用numpy也会更简单:
import numpy as np
a = np.array([[2,1,4],[4,2,1],[5,1,3]])
print(np.linalg.det(a))