python - 递归计算矩阵 (nxn) 的行列式

Question

我将要编写一些代码来计算方阵 (nxn) 的行列式，使用拉普拉斯算法(意思是递归算法)，如 Wikipedia's Laplace Expansion 所写.

我已经有了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 矩阵 - 为什么？)

Answer 1

您确定您的 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