递归计算矩阵 (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屋！