4

Description(this is a hwk question):

I am not sure where to start here. I plan to use Laplace's Expansion but I am not sure how to implement it for nxn matrices. Any help would be appreciated.

Note: I already have a function to generate random matrices for a nxn matrix. Also the timing the calculation isn't a problem. The only thing I have an issue is how to calculate the determinant.

Had to delete the question description b/c of my class policy.

theAngryPhysicist
  • 129
  • 1
  • 5
  • 12

5 Answers5

8

Here is recursive python code for the adjucate method for finding a matrix's determinant. 

def getMatrixMinor(m,i,j):
    return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])]

def getMatrixDeternminant(m):
    #base case for 2x2 matrix
    if len(m) == 2:
        return m[0][0]*m[1][1]-m[0][1]*m[1][0]

    determinant = 0
    for c in range(len(m)):
        determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c))
    return determinant

Note that the input is an array of arrays representing the nxn matrix

stackPusher
  • 6,076
  • 3
  • 35
  • 36
3

Ok, here's a hint.

  1. write a function to calculate the minor matrices. (hint, use slices)
  2. write a function to calculate the cofactors (this should call the first function, and the determinate function)
  3. the determinate function calls the function in step two and adds the results together. (hint: use sum)

viola, you have a determinant.

Also, don't forget that because of the way we write lists in python, the indices get reversed. That is if 

M = [[1, 2, 3],
     [4, 5, 6],
     [7, 8, 9]]

then m0,1 is 2 and not 4 as it would be in normal notation. you can think of it as a transpose or use zip

aaronasterling
  • 68,820
  • 20
  • 127
  • 125
  • So for the indexes the columns are counted first and then the rows? – theAngryPhysicist Sep 29 '10 at 07:46
  • @shagster no, first the rows then the columns. Look at the code and think about it. (for me at least, the normal way to indice is column, row. "First you go down the hall then you go up the elevator") it's really a moot point here though because D(A^t) = D(A) – aaronasterling Sep 29 '10 at 07:52
1
def minor(array,i,j):
    c = array
    c = c[:i] + c[i+1:]
    for k in range(0,len(c)):
        c[k] = c[k][:j]+c[k][j+1:]
    return c
def det(array,n):
    if n == 1 :return array[0][0]
    if n == 2 :return array[0][0]*array[1][1] - array[0][1]*array[1][0]
    sum = 0
    for i in range(0,n):
        m = minor(array,0,i)
        sum =sum + ((-1)**i)*array[0][i] * det(m,n-1)
    return sum
Peyman Naseri
  • 11
  • 1
1

For large matrices, it is not recommended to use Laplace's method for determinant calculation, since it is computationally expensive (recursive functions). Instead, a better approach is to use the Gauss Elimination method to convert the original matrix into an upper triangular matrix. The determinant of a lower or an upper triangular matrix is simply the product of the diagonal elements.

Here we show an example.

import numpy as np
from scipy import linalg

def determinant(a):
    assert len(a.shape) == 2 # check if a is a two diamentional matrix
    assert a.shape[0] == a.shape[1] # check if matrix is square
    n = a.shape[0]
   
    for k in range(0, n-1):
       
        for i in range(k+1, n):
            if a[i,k] != 0.0:
                lam = a [i,k]/a[k,k]
                a[i,k:n] = a[i,k:n] - lam*a[k,k:n]

               
    # the matrix (a) is now in the upper triangular form
                
    return np.prod(np.diag(a))

matrix = np.random.rand(50, 50)


print(linalg.det(matrix))
print(determinant(matrix))

The result is similar to the one obtained from the Scipy determinant!

-3032.573716363944

-3032.573716915967

However, the function can still be formulated to include pivoting. Additionally, it can be brought to a similar efficiency as the one from scipy by using numba library.

Khalil Al Hooti
  • 4,207
  • 5
  • 23
  • 40
0

I've used as base idea of Peyman Naseri and want to share with my implementation,

import numpy as np
import time

def minor(arr,i,j):
    c = arr[:]
    c = np.delete(c, (i),axis=0)
    return [np.delete(row, (j),axis=0) for row in (c)]
    

def det(arr):
    n = len(arr)
    if n == 1 :return arr[0][0]
    if n == 2 :return arr[0][0]*arr[1][1] - arr[0][1]*arr[1][0]
    sum = 0
    for i in range(0,n):
        m = minor(arr,0,i)
        sum =sum + ((-1)**i)*arr[0][i] * det(m)
    return sum

matrix = np.random.randint(-5, 5, size=(10, 10)) # martix nxn with integer values in interval [-5, 5)

print(matrix)
start_time = time.time()
print("started:", start_time)
det(matrix)
end_time = time.time()
print("finished:", end_time)
print("--- %s seconds ---" % (end_time - start_time))

as result determinant calculation for matrix 10*10 on my laptop takes about 1 minute, I understand that my code not optimal but main reason for this implementation (maybe I lost something) I just need to get working code base on Peyman Naseri solution which seems to me very pretty.

[[ 2 -1 -1  0 -2  0  4  4  3  4]
 [-3  1 -1  3  0 -3 -2  0  3 -1]
 [ 2 -1 -4  3  0 -2 -2 -5  3 -5]
 [-2 -1  2 -2  4 -3 -5 -1 -5  3]
 [ 1 -4  1 -5 -5  4 -3 -5  3  1]
 [ 2  4  0 -1 -1 -5 -2 -2 -3 -5]
 [ 1  4 -3 -4 -5  0  0  0 -5 -1]
 [ 0 -5 -5  4 -3 -2  2 -4  2 -5]
 [-3  1 -1 -4  4 -5  3 -3 -4  0]
 [ 0 -2  2 -3  1  3  2  0 -1  4]]
started: 1636896388.0213237
finished: 1636896442.846928
--- 54.82560420036316 seconds ---
GoAl
  • 1
  • 1