递归计算矩阵 (nxn) 的行列式 [英] computing determinant of a matrix (nxn) recursively
问题描述
我将要编写一些代码来计算方阵 (nxn) 的行列式,使用 Laplace 算法(含义递归算法)作为编写 维基百科的拉普拉斯扩展.
I'm about to write some code that computes the determinant of a square matrix (nxn), using the Laplace algorithm (Meaning recursive algorithm) as written Wikipedia's Laplace Expansion.
我已经有了 Matrix
类,其中包括 init、setitem、getitem、repr 以及计算行列式所需的所有内容(包括 minor(i,j)
).
I already have the class Matrix
, which includes init, setitem, getitem, repr and all the things I need to compute the determinant (including minor(i,j)
).
所以我尝试了下面的代码:
So I've tried the code below:
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 矩阵 - 为什么?)
As expected, in every recursive call, self
turns into an appropriate minor. But when coming back from the recursive call, it doesn't change back to it's original matrix.
This causes trouble when in the for
loop (when the function arrives at (n,m)==(1,1)
, this one value of the matrix is returned, but in the for
loop, self
is still a 1x1 matrix - why?)
推荐答案
您确定您的 minor
返回的是新对象而不是对原始矩阵对象的引用吗?我使用了你的确切行列式方法并为你的类实现了一个 minor
方法,它对我来说很好用.
Are you sure that your minor
returns the a new object and not a reference to your original matrix object? I used your exact determinant method and implemented a minor
method for your class, and it works fine for me.
以下是您的矩阵类的快速/脏实现,因为我没有您的实现.为简洁起见,我选择仅对方阵实施它,在这种情况下,这无关紧要,因为我们正在处理行列式.注意det
方法,和你的一样,还有minor
方法(剩下的方法都是为了方便实现和测试):
Below is a quick/dirty implementation of your matrix class, since I don't have your implementation. For brevity I have chosen to implement it for square matrices only, which in this case shouldn't matter as we are dealing with determinants. Pay attention to det
method, that is the same as yours, and minor
method (the rest of the methods are there to facilitate the implementation and testing):
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
精度:
The residual in case of numpy is because it calculates the determinant through (Gaussian) elimination method rather than the Laplace expansion. You can also compare the results on random matrices to see that the difference between your determinant function and numpy's doesn't grow beyond float
precision:
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
这篇关于递归计算矩阵 (nxn) 的行列式的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!