Matrix Determinant And Inverse
May 23, 2017
Today’s exercise is preliminary to the exercise we will have later this week. You are to write programs that calculate the determinant and inverse of a matrix. I won’t go into the math involved behind the matrix arithmetic, as there are many sources on the internet that know far more about the process than I. Google for “matrix determinant” or “matrix inverse”; I used YouTube for my instruction.
Your task is to write programs that calculate the determinant and inverse of a matrix. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
No worse than addend, summand, subtrahend, minuend, multiplicand.
In Python, the determinant part.
def determinant(matrix): if len(matrix) == 2: return (matrix[0][0]*matrix[1][1]-matrix[0][1]*matrix[1][0]) else: negator = 1 result = 0 for col in range(len(matrix)): sub_matrix = [[] for i in range(len(matrix)-1)] for sub_row in range(1,len(matrix)): for sub_col in range(len(matrix)): if sub_col != col: sub_matrix[sub_row-1].append(matrix[sub_row][sub_col]) result += negator * matrix[0][col] * determinant(sub_matrix) negator *= -1 return result print(-240 == determinant([[1,-4,9],[-6,7,3],[1,2,3]])) print(441 == determinant([[6,24,1],[13,16,10],[20,17,15]]))I’ve found these two links with explanations and code. It’s quite a hassle for larger dimensions….
Inverse of a matrix: http://www.geeksforgeeks.org/adjoint-inverse-matrix/
Determinant of a matrix: http://www.geeksforgeeks.org/determinant-of-a-matrix/
Though, I haven’t tried their code yet.
Here in Python an implementation of an LU decomposition. Once the matrix is decomposed, the determinant and inverse are straightforward.
from copy import deepcopy from functools import reduce from operator import mul def trans(A): return [list(a) for a in zip(*A)] def dot(a, b): return sum(ai*bi for ai, bi in zip(a, b)) def matvec(A, b): return [dot(a, b) for a in A] def matmul(A, B): return trans(matvec(A, b) for b in zip(*B)) def zeros(m, n): return [[0]*n for i in range(m)] def lu_decompose(A, tol): """ returns matrix LU with upper and lower matrix vector P: P[:n] contains the ordering of rows P[n] = n + numbeer of swaps in P[:n] """ LU = deepcopy(A) # the original matrix is not changed n = len(LU) P = list(range(n+1)) for i in range(n): imax = max(range(i, n), key=lambda m: abs(LU[m][i])) if abs(LU[i][imax]) < tol: raise ValueError("Matrix is degenerate") if imax != i: P[i], P[imax] = P[imax], P[i] LU[i], LU[imax] = LU[imax], LU[i] P[n] += 1 for j in range(i+1, n): LU[j][i] /= LU[i][i] for k in range(i+1, n): LU[j][k] -= LU[j][i] * LU[i][k]; return LU, P def lu_solve(LU, P, b): """solve Ax=b A is te original matrix""" n = len(LU) x = [b[p] for p in P[:-1]] for i in range(n): x[i] = b[P[i]] for k in range(i): x[i] -= LU[i][k] * x[k] for i in range(n-1, -1, -1): for k in range(i + 1, n): x[i] -= LU[i][k] * x[k] x[i] /= LU[i][i] return x def lu_invert(LU, P): """returns the inverse of the original matrix""" n = len(LU) IA = zeros(n, n) for j in range(n): for i in range(n): IA[i][j] = 1.0 if P[i] == j else 0.0 for k in range(i): IA[i][j] -= LU[i][k] * IA[k][j] for i in range(n-1, -1, -1): for k in range(i+1, n): IA[i][j] -= LU[i][k] * IA[k][j] IA[i][j] /= LU[i][i] return IA def lu_determinant(LU, P): """returns the determinant of the original matrix""" n = len(LU) det = reduce(mul, (LU[i][i] for i in range(n))) return det if (P[n]-n) % 2 == 0 else -det def print_matrix(A): print() for row in A: for ai in row: print("{:8.4f}".format(ai), end=" ") print() b = [[6, 24, 1], [13, 16, 10], [20, 17, 15]] A, P = lu_decompose(b, 1e-8) ib = lu_invert(A, P) print_matrix(ib) print_matrix(matmul(b, ib)) print(lu_determinant(A, P))Klong (for determinant)
steve@steve-Satellite-L555D:~$ rlwrap kg Klong 20161212 :" Beginning of code" det::{:[2=#x; cal1(x); cal2(x)]} cal1::{((x:@[0 0])*(x:@[1 1]))-((x:@[0 1])*(x:@[1 0]))} mults::{[a b]; a::1; b::[]; {b::b,(x*a); a::-a}'*x; b} makerow::{[a pos rn row t]; a::[]; t::x; rn::y; pos::z; row::t@rn; ((1+&#row):=0,pos){:[x=1; a::a,y; a]}'row; a} maketbl::{[pos tbl]; tbl::x; pos::y; {makerow(tbl; x; pos)}'(1+!(#tbl)-1)} maketbls::{[t]; t::x; {maketbl(t; x)}'!#t} cal2::{[a t]; a::0; t::x; mults(t){a::a+(x*det(y))}'maketbls(t); a} :" End of code" t2::[[1 2] [3 4]] t3::[[1 2 3] [2 3 4] [3 4 5]] t4::[[1 2 3 4] [2 3 4 5] [3 4 55 666] [7 88 999 1234]] :" Should be -2" det(t2) -2 :" Should be 0" det(t3) 0 :" Should be 498600" det(t4) 498600 :" Should be 18 per Wikipedia" det([[-2 2 -3] [-1 1 3] [2 0 -1]]) 18 :" Per @Rutger should be -240" det([[1 -4 9] [-6 7 3] [1 2 3]]) -240 :" Per @Rutger shuld be 441" det([[6 24 1] [13 16 10] [20 17 15]]) 441