4

我将编写一些代码来计算方阵 (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 矩阵 - 为什么？）

4

5 回答 5

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
于 2016-12-22T01:01:34.813 回答
0 
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)
于 2020-09-16T22:59:36.880 回答
0

这是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
于 2016-01-20T02:14:30.300 回答
0

在下面的链接中使用 Sarrus 的规则（非递归方法）示例是用 Javascript 编写的，但可以很容易地用 python https://github.com/apanasara/Faster_nxn_Determinant编写

于 2019-12-14T19:04:55.417 回答
-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

很抱歉没有在伙计们之前发表评论：）如果您需要任何进一步的解释，请随时询问。

于 2016-12-21T13:35:11.057 回答