Compute the determinant
Challenge
A simple challenge: Given a two-dimensional matrix (an array of arrays) of real numbers, compute the determinant.
The determinant of a matrix is a mathematical construct used in many applications, such as solving polynomial equations, identifying the invertibility of the matrix, and finding the scaling factor under a matrix transformation. For more information about it, see this Wikipedia entry.
There are a couple of different ways to compute the determinant, and it is up to you how you implement it.
For instance, you may compute it using the Laplace expansion, a recursive algorithm which goes
- Pick a row or column.
- Start with a sum of zero.
- For each entry of the row/column:
- Create a new matrix with the row and column of the entry removed. This new matrix is a square matrix of size one less that the original.
- Compute the determinant of that smaller matrix.
- Multiply that determinant by the entry.
- If the row index plus the column index is even[1], add it to the sum, otherwise, subtract it.
- The final sum is the determinant.
As an example, here is an ungolfed implementation along the first column.
function laplaceDet(matrix) {
if (matrix.length === 1) return matrix[0][0];
let sum = 0;
for (let rowIndex = 0; rowIndex < matrix.length; ++rowIndex) {
let minorMatrix = matrix.filter((_, index) => index !== rowIndex)
.map(row => row.slice(1));
sum += ((-1) ** rowIndex) * matrix[rowIndex][0] * laplaceDet(minorMatrix);
}
return sum;
}
This is code-golf, so the program with the lowest byte-count wins!
Test cases
$$ \begin{aligned} \det\begin{bmatrix}1&0\\0&1\end{bmatrix}&=1 \\ \det\begin{bmatrix}1.5&2\\-3&4.5\end{bmatrix}&=12.75 \\ \det\begin{bmatrix}3&7\\1&-4\end{bmatrix}&=-19 \\ \det\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}&=1 \\ \det\begin{bmatrix}1&2&3\\4&5&6\\7&8&9\end{bmatrix}&=0 \end{aligned} $$
In text form,
[[1,0],[0,1]] -> 1
[[1.5,2],[-3,4.5]] -> 12.75
[[3,7],[1,-4]] -> -19
[[1,0,0],[0,1,0],[0,0,1]] -> 1
[[1,2,3],[4,5,6],[7,8,9]] -> 0
-
Note that it doesn't matter if it is 1-indexed or 0-indexed, as 1+1 and 0+0 are both the same parity. ↩︎
CJam, 45 bytes {:A, …
3y ago
[Python 3.8 (pre-release)], 10 …
3y ago
[Haskell], 75 73 bytes -2 b …
3y ago
C, 147 bytes ```c float d( …
3y ago
Scala, 130 125 bytes Saved …
3y ago
J, 5 bytes ``` -/ . ``` …
3y ago
[Python 3], 29 bytes …
3y ago
[Wolfram Language (Mathematica …
2y ago
Ruby, 35 bytes ->m{requi …
3y ago
9 answers
You are accessing this answer with a direct link, so it's being shown above all other answers regardless of its score. You can return to the normal view.
Python 3, 29 bytes
import numpy
numpy.linalg.det
Yeah, this is really boring, but I'm curious if this can be beat in Python.
3 comment threads
CJam, 45 bytes
{:A_,{1$_,,.=:+\)/:CAff*A@zf{\f.*::+}..-}/;C}
This implementation is an anonymous block (~function). Online test suite.
Dissection
This implements the Faddeev-LeVerrier algorithm.
The objective is to calculate the coefficients $c_k$ of the characteristic polynomial of the $n\times n$ matrix $A$, $$\begin{eqnarray}p(\lambda)\equiv \det(\lambda I_{n}-A) = \sum_{k=0}^{n}c_{k}\lambda^{k}\end{eqnarray}$$ where, evidently, $c_n = 1$ and $c_0 = (-1)^n \det A$.
The coefficients are determined recursively from the top down, by dint of the auxiliary matrices $M$, $$\begin{aligned}M_{0}&\equiv 0&c_{n}&=1\qquad &(k=0)\\M_{k}&\equiv AM_{k-1}+c_{n-k+1}I\qquad \qquad &c_{n-k}&=-{\frac {1}{k}}\mathrm {tr} (AM_{k})\qquad &k=1,\ldots ,n~.\end{aligned}$$
The code never works directly with $c_{n-k}$ and $M_k$, but always with $(-1)^k c_{n-k}$ and $(-1)^{k+1}AM_k$, so the recurrence is $$\begin{eqnarray*}(-1)^k c_{n-k} &=& \frac1k \mathrm{tr} ((-1)^{k+1} AM_{k}) \\ (-1)^{k+2} AM_{k+1} &=& (-1)^k c_{n-k}A - A((-1)^{k+1}A M_k)\end{eqnarray*}$$
{ e# Define a block
:A e# Store the input matrix in A
_, e# Take the length of a copy
{ e# for i = 0 to n-1
e# Stack: (-1)^{i+2}AM_{i+1} i
1$_,,.=:+ e# Calculate tr((-1)^{i+2}AM_{i+1})
\)/:C e# Divide by (i+1) and store in C
Aff* e# Multiply by A
A@ e# Push a copy of A, bring (-1)^{i+2}AM_{i+1} to the top
zf{\f.*::+} e# Matrix multiplication
..- e# Matrix subtraction
}/
; e# Pop (-1)^{n+2}AM_{n+1} (which incidentally is 0)
C e# Fetch the last stored value of C
}
0 comment threads
Python 3.8 (pre-release), 106 95 bytes
Saved 11 bytes thanks to Peter Taylor!
f=lambda m:sum((-1)**i*x*f([r[:i]+r[i+1:]for r in m[1:]])for i,x in enumerate(m[0]))if m else 1
Here's an attempt at beating Quintec's answer. Only 78 67 more bytes to go!
Haskell, 75 73 bytes
-2 bytes thanks to @user
f[]=1
f(x:y)=foldr(-)0$zipWith(\i->(*f[take i r++drop(i+1)r|r<-y]))[0..]x
C, 147 bytes
float d(float**m,int r){if(r<2)return**m;int i=r,j;float s=0,*n[--r];for(;i--;s+=(i%2?-1:1)**m[i]*d(n,r))for(j=r;j--;)n[j]=m[j+(j>=i)]+1;return s;}
Basically a C version of my example code. Takes an array of pointers and the size.
Prettified and commented version:
float d(float **m, int r) {
if(r < 2) return **m; // Base case: 1x1 matrix
int i = r, j; // Variable initialization, i and j are counters
float s = 0, *n[--r]; // s is the sum, n is the minor matrix
for(; i--; s += (i % 2 ? -1 : 1) * *m[i] * d(n, r)) // Outer loop: loop over the rows of M.
// After each iteration, add the term to the sum
for(j = r; j--; ) // Inner loop to fill in the minor matrix
n[j] = m[j + (j >= i)] + 1; // (j >= i) term skips over the ith row of M.
// Addition of 1 effectively removes the first element
return s;
}
0 comment threads
Scala, 130 125 bytes
Saved 5 bytes by returning 1 like Hakerh400's great answer
def f(m:Seq[Seq[Double]]):Double=if(m.size<1)1 else(m.indices:\.0){(i,a)=>m(0)(i)*f(m.tail.map(r=>r.take(i)++r.drop(i+1)))-a}
A recursive function that uses the first row for the Laplace expansion. Explanation coming soon.
def f(m: Seq[Seq[Double]]): Double =
if (m.size < 1) 1 //If it's empty, return 1
else
//Otherwise, fold right over [0..m.size-1]
//The initial accumulator is 0.0
(m.indices :\ .0) {
(i, a) =>
//Multiply the entry by
m(0)(i) *
//The determinant of the smaller matrix
f(
//Drop the first row
m.tail.map(
//And for each row r, drop column i (0-indexed)
r => r.take(i) ++ r.drop(i + 1)
)
)
//Subtract the accumulator from that to invert
//its sign each time
- a
}
0 comment threads
1 comment thread