递归计算一个矩阵（nxn）的行列式

Question

我即将编写一些代码，用拉普拉斯算法（即递归算法），计算一个方阵（nxn）的行列式，具体实现可以参考维基百科的 拉普拉斯展开。

我已经有了类 Matrix，并包括初始化、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的矩阵 - 为什么？）

- user2375340

---

self.minor(i,j) 实际上是否会改变 self？ - Radio-

是的，你说得对。我会更改并再次检查。谢谢。 - user2375340

我建议不要像这样使用“asserts”。第一个（assert isinstance(self, Matrix)）没有任何作用 - 它是Matrix上的一个方法，所以您确定self是一个Matrix。第二个会更好地使用类似于``if n！= m：raise ValueError（"Matrix is not square"）`的东西，因为异常更具体地说明了问题（并且可以更容易地被捕获）。 - Benjamin Hodgson

另外，您真的希望用户能够将“i”作为参数传递吗？行列式是矩阵的一个属性，与您选择哪一行或哪一列无关。因此，让用户选择并没有什么意义（因为无论他们选择什么，结果都将是相同的）。 - Benjamin Hodgson

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

- s.yadegari

Answer 2

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)

- Tanveer Mahtab

Answer 3

这是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

- Advanderar

Answer 4

使用Sarrus法则（非递归方法）。 下面链接中的示例是用Javascript编写的，但很容易改写成Python。 https://github.com/apanasara/Faster_nxn_Determinant

- Amit Panasara

Answer 5

我发布了这段代码，因为在互联网上找不到如何仅使用标准库解决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

抱歉之前没有评论，伙计们 :) 如果您需要进一步的解释，请不要犹豫问。

- yasser zairi

---

在SO上，仅发布代码答案而不注释或解释其如何解决问题是一种不鼓励的做法。 - Akshat Mahajan