首页 文章

以递归方式计算矩阵(nxn)的行列式

提问于
浏览 1702
4

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

我已经有了类 Matrix ,其中包括 initsetitemgetitemrepr 以及计算行列式所需的所有东西(包括 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矩阵 - 为什么?)

4 回答

  • 1

    这是python 3中的函数 .

    注意:我使用一维列表来容纳矩阵,而size数组是方形数组中的行数或列数 . 它使用递归算法来查找行列式 .

    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
    
  • 0

    嘿,我已经使用递归函数在MATLAB中编写了代码 . 这可能对您有所帮助 .

    function value = determinant(A)
    % Calculates determinant of a square matrix A.
    % This is a recursive function. Not suitable for large matrices.
    
    [rows, columns] = size(A);
    if rows ~= columns
        error('input matrix is not a square matrix.')
    end
    
    value = 0;
    if rows == 2
        for i = 1:rows
            value = A(1,1)*A(2,2) - A(1,2)*A(2,1);
        end
    else
        for i = 1:rows
            columnIndices = [1:i-1 i+1:rows];
            value = value + (-1)^(i+1)*A(1,i)*... 
                determinant(A(2:rows, columnIndices));
        end
    end
    end
    
  • -2

    你确定 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 's doesn' t之间的差异超出 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
    
  • -2

    我发布了这段代码,因为我无法在互联网上使用它,如何使用标准库来解决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
    

    对不起之前没有评论:)如果你需要任何进一步的解释,请不要犹豫 .

相关问题