computing determinant of a matrix (nxn) recursively

Question

I'm about to write some code that computes the determinant of a square matrix (nxn), using the Laplace algorithm (Meaning recursive algorithm) as written Wikipedia's Laplace Expansion.

I already have the class Matrix, which includes init, setitem, getitem, repr and all the things I need to compute the determinant (including minor(i,j)).

So I've tried the code below:

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

As expected, in every recursive call, self turns into an appropriate minor. But when coming back from the recursive call, it doesn't change back to it's original matrix. This causes trouble when in the for loop (when the function arrives at (n,m)==(1,1), this one value of the matrix is returned, but in the for loop, self is still a 1x1 matrix - why?)

Answer 1

Are you sure that your minor returns the a new object and 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 fine for me.

Below is a quick/dirty implementation of your matrix class, since I don't have your implementation. For brevity I have chosen to implement it for square matrices only, which in this case shouldn't matter as we are dealing with determinants. Pay attention to det method, that is the same as yours, and minor method (the rest of the methods are there to facilitate the 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 residual in case of numpy is because it calculates the determinant through (Gaussian) elimination method rather than the Laplace expansion. You can also compare the results on random matrices to see that the difference between your determinant function and numpy's doesn't grow beyond float precision:

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

use Sarrus' Rule (non recursive method) example on below link is in Javascript, but easily can be written in python https://github.com/apanasara/Faster_nxn_Determinant

Answer 3

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)

Answer 4

Here's the function in python 3.

Note: I used a one-dimensional list to house the matrix and the size array is the amount 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 5

i posted this code because i couldn't fine it on the internet, how to solve n*n determinant using only standard library. the purpose is to share it with those who will find it useful. i started by calculating the submatrix Ai related to a(0,i). and i used recursive determinant 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 for not commenting before guys :) if you need any further explanation don't hesitate to ask .