遞歸計算矩陣（nxn）的行列式

Question

我正在編寫一些代碼，使用拉普拉斯算法（含義遞歸算法）計算方陣矩陣（nxn）的行列式，如Wikipedia's Laplace Expansion所述。

我已經有了類Matrix，其中包括初始化，setitem，的GetItem，再版和所有我需要計算行列式的東西（包括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

嘿，我已經寫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 

Answer 2

下面是在python的功能3.

注：I中使用的一維列表以容納基質和尺寸數組是行或列的正方形陣列中的量。它使用遞歸算法來查找行列式。

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 

Answer 3

我發佈此代碼，因爲我無法在互聯網上罰款，如何解決使用標準庫n * n行列式。 的目的是與那些會發現它有用的人分享。 我通過計算與（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 

對不起球員之前不評論:) 如果您需要任何進一步的解釋不要猶豫，問。

Answer 4

你確定你的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