我要写一些代码来计算方阵(nxn)的行列式，使用维基百科拉普拉斯展开所写的拉普拉斯算法(意思是递归算法)。
我已经有了Matrix类，其中包括init， setitem， getitem， repr和所有我需要计算行列式的东西(包括minor(i,j))。
所以我尝试了下面的代码:
def determinant(self,i=0) # i can be any of the matrix's rows
assert isinstance(self,Matrix)
n,m = self.dim() # Q.dim() returns the size of the matrix Q
assert n == m
if (n,m) == (1,1):
return self[0,0]
det = 0
for j in range(n):
det += ((-1)**(i+j))*(self[i,j])*((self.minor(i,j)).determinant())
return det
正如预期的那样，在每次递归调用中，self都变成了适当的次要变量。但是当从递归调用返回时，它不会变回原来的矩阵。这在for循环中会引起麻烦(当函数到达(n,m)==(1,1)时，返回矩阵的这一个值，但在for循环中，self仍然是1x1矩阵-为什么?)
您确定您的minor返回一个新对象，而不是对原始矩阵对象的引用吗?我使用了您的确切行列式方法，并为您的类实现了minor方法，它对我来说很好。
下面是你的矩阵类的快速/肮脏的实现，因为我没有你的实现。为简洁起见，我选择只在方阵中实现它，在这种情况下，我们处理行列式时，这应该无关紧要。注意det方法，它与您的方法相同，以及minor方法(其余的方法是为了方便实现和测试):
class matrix:
def __init__(self, n):
self.data = [0.0 for i in range(n*n)]
self.dim = n
@classmethod
def rand(self, n):
import random
a = matrix(n)
for i in range(n):
for j in range(n):
a[i,j] = random.random()
return a
@classmethod
def eye(self, n):
a = matrix(n)
for i in range(n):
a[i,i] = 1.0
return a
def __repr__(self):
n = self.dim
for i in range(n):
print str(self.data[i*n: i*n+n])
return ''
def __getitem__(self,(i,j)):
assert i < self.dim and j < self.dim
return self.data[self.dim*i + j]
def __setitem__(self, (i, j), val):
assert i < self.dim and j < self.dim
self.data[self.dim*i + j] = float(val)
#
def minor(self, i,j):
n = self.dim
assert i < n and j < n
a = matrix(self.dim-1)
for k in range(n):
for l in range(n):
if k == i or l == j: continue
if k < i:
K = k
else:
K = k-1
if l < j:
L = l
else:
L = l-1
a[K,L] = self[k,l]
return a
def det(self, i=0):
n = self.dim
if n == 1:
return self[0,0]
d = 0
for j in range(n):
d += ((-1)**(i+j))*(self[i,j])*((self.minor(i,j)).det())
return d
def __mul__(self, v):
n = self.dim
a = matrix(n)
for i in range(n):
for j in range(n):
a[i,j] = v * self[i,j]
return a
__rmul__ = __mul__
现在开始测试
import numpy as np
a = matrix(3)
# same matrix from the Wikipedia page
a[0,0] = 1
a[0,1] = 2
a[0,2] = 3
a[1,0] = 4
a[1,1] = 5
a[1,2] = 6
a[2,0] = 7
a[2,1] = 8
a[2,2] = 9
a.det() # returns 0.0
# trying with numpy the same matrix
A = np.array(a.data).reshape([3,3])
print np.linalg.det(A) # returns -9.51619735393e-16
numpy的残差是因为它通过(高斯)消去法而不是拉普拉斯展开来计算行列式。您还可以比较随机矩阵的结果，以查看您的行列式函数和numpy的行列式函数之间的差异不会超过float精度:
import numpy as np
a = 10*matrix.rand(4)
A = np.array( a.data ).reshape([4,4])
print (np.linalg.det(A) - a.det())/a.det() # varies between zero and 1e-14
使用Sarrus规则(非递归方法)下面链接的示例是用Javascript编写的，但可以很容易地用python编写https://github.com/apanasara/Faster_nxn_Determinant
import numpy as np
def smaller_matrix(original_matrix,row, column):
for ii in range(len(original_matrix)):
new_matrix=np.delete(original_matrix,ii,0)
new_matrix=np.delete(new_matrix,column,1)
return new_matrix
def determinant(matrix):
"""Returns a determinant of a matrix by recursive method."""
(r,c) = matrix.shape
if r != c:
print("Error!Not a square matrix!")
return None
elif r==2:
simple_determinant = matrix[0][0]*matrix[1][1]-matrix[0][1]*matrix[1][0]
return simple_determinant
else:
answer=0
for j in range(r):
cofactor = (-1)**(0+j) * matrix[0][j] * determinant(smaller_matrix(matrix, 0, j))
answer+= cofactor
return answer
#test the function
#Only works for numpy.array input
np.random.seed(1)
matrix=np.random.rand(5,5)
determinant(matrix)
下面是python 3中的函数。
注意:我使用一维列表来存放矩阵，数组的大小是正方形数组中的行或列的数量。它使用递归算法来求行列式。
def solve(matrix,size):
c = []
d = 0
print_matrix(matrix,size)
if size == 0:
for i in range(len(matrix)):
d = d + matrix[i]
return d
elif len(matrix) == 4:
c = (matrix[0] * matrix[3]) - (matrix[1] * matrix[2])
print(c)
return c
else:
for j in range(size):
new_matrix = []
for i in range(size*size):
if i % size != j and i > = size:
new_matrix.append(matrix[i])
c.append(solve(new_matrix,size-1) * matrix[j] * ((-1)**(j+2)))
d = solve(c,0)
return d
我张贴了这个代码，因为我不能在互联网上罚款，如何解决n*n行列式只使用标准库。目的是与那些觉得有用的人分享。我从计算a(0,i)的子矩阵Ai开始。我用递归行列式来简化。
def submatrix(M, c):
B = [[1] * len(M) for i in range(len(M))]
for l in range(len(M)):
for k in range(len(M)):
B[l][k] = M[l][k]
B.pop(0)
for i in range(len(B)):
B[i].pop(c)
return B
def det(M):
X = 0
if len(M) != len(M[0]):
print('matrice non carrée')
else:
if len(M) <= 2:
return M[0][0] * M[1][1] - M[0][1] * M[1][0]
else:
for i in range(len(M)):
X = X + ((-1) ** (i)) * M[0][i] * det(submatrix(M, i))
return X
很抱歉没有在大家之前评论:)如果你需要进一步的解释，请尽管问。