Computational determinant of a matrix (nxn) recursively

Question

I am going to write code that computes the determinant of a square matrix (nxn) using the Laplace algorithm (value of a recursive algorithm), as written by Wikipedia The Laplace extension .

I already have a Matrix class that includes init , setitem , getitem , repr and everything I need to calculate the determinant (including minor(i,j) ).

So, I tried the code below:

 def determinant(self,i=0) # i can be any of the matrix 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 

As expected, with each recursive call, self turns into a suitable minor. But, returning from a recursive call, it does not return to the original matrix. This causes problems when in a for loop (when a function arrives at (n,m)==(1,1) , this one matrix value is returned, but in a for self loop it is still a 1x1 matrix - why?)

Answer 1

Are you sure your minor returns a new object, 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 great for me.

Below is a short-term implementation of your matrix class, since I do not have your implementation. For brevity, I decided to implement it only for square matrices, which in this case does not matter, since we are dealing with determinants. Pay attention to the det method, which matches yours, and the minor method (other methods are designed to facilitate 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__ 

Now for testing

 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 

The remainder in the case of numpy is that it calculates the determinant by the (Gaussian) method, not the Laplace extension. You can also compare the results on random matrices to see that the difference between your determinant function and numpy does not exceed the precision of the 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 

Answer 2

Here is the function in python 3.

Note. I used a one-dimensional list to place the matrix, and the size array is the number of rows or columns in the square array. It uses a recursive algorithm to find the determinant.

 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

Hi, I wrote code in MATLAB using a recursive function. This may be useful for you.

 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 4

I posted this code because I could not distinguish it on the Internet, how to solve the n * n qualifier using only the standard library. the goal is to share it with those who find it useful. I started by calculating the submatrix Ai associated with a (0, i). and I used a recursive qualifier to make it short.

  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 

sorry to not comment in front of the guys :) if you need any further explanation, feel free to ask.