Bird's Algorithm for Computing Determinants

There is a division-free algorithm for computing determinants published by R.S.Bird in 2011 that uses only matrix multiplications. Given a $$\n×n\$$ matrix $$\X\$$, the matrix $$\Y=μ(X)\$$ is another $$\n×n\$$ matrix which entries are given by

$$Y_{i,j} = \begin{cases} 0 & \text{ if } j < i \\ X_{i,j} & \text{ if } j > i \\ -(X_{i+1,i+1} + X_{i+2,i+2} + \cdots + X_{n,n}) & \text{ if } j = i \end{cases}$$

Thus entries of $$\X\$$ below the diagonal are made zero, those above the diagonal are left unchanged, and each diagonal entry is replaced by the negated sum of the elements in the diagonal below it. Note that $$\Y_{n,n} = 0\$$.

Define the operation $$\A⊗X\$$ between two $$\n×n\$$ matrices $$\A\$$ and $$\X\$$ as $$A⊗X=-μ(X) \cdot A,$$ where $$\\cdot\$$ denotes matrix multiplication.

Now, to obtain $$\\det A\$$, the determinant of $$\A\$$, compute

$$A⊗(A⊗(A⊗(...⊗A)))$$

where the operation $$\⊗\$$ is applied $$\n-1\$$ times. The resulting matrix is everywhere zero except for its leading entry, which equals $$\\det A\$$.

Example

$$A=\left[\begin{array}{rrr}-2 & -4 & -1 \\0 & 1 & 3 \\-2 & 3 & 1 \\\end{array}\right]$$ $$\mu(A)=\left[\begin{array}{rrr}-2 & -4 & -1 \\0 & -1 & 3 \\0 & 0 & 0 \\\end{array}\right]$$ $$B=A⊗A=-\mu(A)\cdot A= \left[\begin{array}{rrr}-6 & -1 & 11 \\6 & -8 & 0 \\0 & 0 & 0 \\\end{array}\right]$$ $$\mu(B)=\left[\begin{array}{rrr}8 & -1 & 11 \\0 & 0 & 0 \\0 & 0 & 0 \\\end{array}\right]$$ $$A⊗(A⊗A)=A⊗B=-\mu(B)\cdot A=\left[\begin{array}{rrr}38 & 0 & 0 \\0 & 0 & 0 \\0 & 0 & 0 \\\end{array}\right]$$ $$\det A=38$$

Challenge

Implement Bird's algorithm for computing determinants. For a given input $$\A\$$ of size $$\n > 1\$$, your program must output the intermediate steps $$\A⊗A\$$, $$\A⊗(A⊗A)\$$, and so on, with $$\n-1\$$ applications of $$\⊗\$$. Standard code-golf rules apply. The shortest code in bytes wins.

Test Cases

A:  [[0,3,1],[0,2,1],[0,2,1]]
A⊗A:    [[0,1,-1],[0,0,0],[0,0,0]]
A⊗A⊗A:  [[0,0,0],[0,0,0],[0,0,0]]

A:  [[0,3,-1],[3,0,-2],[-2,-3,0]]
A⊗A:    [[-11,-3,6],[-4,-6,0],[0,0,0]]
A⊗A⊗A:  [[21,0,0],[0,0,0],[0,0,0]]

A:  [[3,1,0],[-2,1,-2],[3,-1,-2]]
A⊗A:    [[-1,-2,2],[10,-4,0],[0,0,0]]
A⊗A⊗A:  [[-22,0,0],[0,0,0],[0,0,0]]

A:  [[2,3,3,2],[-2,-2,2,0],[-3,-3,1,2],[-3,3,-3,0]]
A⊗A:    [[19,6,-6,-8],[4,4,0,-4],[6,-6,6,0],[0,0,0,0]]
A⊗A⊗A:  [[-10,48,0,32],[-24,0,0,0],[0,0,0,0],[0,0,0,0]]
A⊗A⊗A⊗A:        [[192,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]

A:  [[3,-2,2,0],[0,2,-2,2],[0,2,-3,0],[2,0,3,-2]]
A⊗A:    [[-9,6,-4,4],[-4,-6,-2,-6],[0,-4,6,0],[0,0,0,0]]
A⊗A⊗A:  [[-8,-4,-12,-4],[12,16,0,0],[0,0,0,0],[0,0,0,0]]
A⊗A⊗A⊗A:        [[56,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]

A:  [[-3,-1,0,0],[-1,2,3,2],[3,0,-3,0],[2,0,2,-3]]
A⊗A:    [[11,6,3,2],[-7,-12,-13,-6],[-9,0,9,0],[0,0,0,0]]
A⊗A⊗A:  [[2,-9,-13,-6],[42,18,0,0],[0,0,0,0],[0,0,0,0]]
A⊗A⊗A⊗A:        [[-12,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]

• While an interesting challenge, it would seem hard to verify that a particular answer uses Bird's exact algorithm and not a variation thereof or some other algorithm altogether. Think twice if you're writing a challenge around an algorithm you found on the Internet because there's a consensus to disallow non-observable requirements for code golf challenges.
Nov 3, 2021 at 12:26
• Please consider posting challenges to Sandbox first, so issues like above can be resolved. Nov 3, 2021 at 12:43
• I would suggest to output all the steps, not just the final determinant. Then it is easy to see if it uses the right algorithm or not, as it is part of the challenge output itself. Nov 3, 2021 at 13:35
• May we assume n>=2? Nov 4, 2021 at 9:36
• Yes @pajonk, I think we can agree that a 1x1 matrix is not really a matrix. Nov 4, 2021 at 12:27

K (ngn/k), 43 41 bytes

{{+/'y*/:(u*|+\|+/x*u)-x*++\u:=#x}\x:\:x}


Try it online!

Brutally golfed by ngn.

-2 bytes because it is allowed to output starting with A itself.

input: matrix A with size n * n

{{...}\x:\:x}    scan over the list of n copies of A
{...}            inner fn: x is running matrix, y is A
u:=#x            u: identity matrix of size n
x*++\            x * transpose of row-wise cumsum of u
(keep upper-triangular part of x)
(...)-           subtract the above from...
+/x*u            diagonal entries of x as a flat list
u*|+\|           cumsum in reverse, placed back on the diagonal
the subtraction gives -mu(x)
+/'y*/:          backwards matrix product with A


K (ngn/k), 62 bytes

{1_(#1_x){+/'@[;;:;]'[-y*|\'=#y;!#y;|0,+\|1_y@'!#y]*\:x}[x]\x}


Try it online!

How it works

The first time actually using ternary each (given a ternary function and three input lists, map over the triples) in code golf. It works so well that I can't think of any other way to do the job.

input: matrix A

{1_(#1_x){f}[x]\x}    repeat applying curried function f[A] to A
(size of A - 1) times, and then remove the first item

f:{+/'bigThingOnY*\:x}    matrix multiply -mu(y) with x, where
y is the current running matrix and x is original A

bigThingOnY: @[;;:;]'[-y*|\'=#y;!#y;|0,+\|1_y@'!#y]
@[;;:;]'[u;v;w]    ternary amend each; for each row in u and each number in
v and w, replace v-th value in u with w
u:-y*|\'=#y    negation of y, with all entries below diagonal set to 0
v:!#y          indices of y (when zipped, overwrites diagonal entries)
w:|0,+\|1_y@'!#y    diagonal of y, remove 1st, reverse, cumsum,
prepend zero, and reverse back

• I don't speak K, but wouldn't my approach of cumsum-sum on the diagonal be shorter here as well? Nov 4, 2021 at 13:18
• @pajonk reversed_cumsum-self turns out to be much shorter, probably because reversal is just 1 byte. Nov 4, 2021 at 23:03

R, 11310910396 95 bytes

Or R>=4.1, 88 bytes by replacing the word function with \.

-4 bytes and another -7 thanks to @Dominic van Essen.

function(A,B=A)for(i in 2:nrow(A))B=print((diag(sum(b<-diag(B))-cumsum(b))-B*upper.tri(B))%*%A)


Try it online!

Straightforward approach.

The terms on the diagonal are calculated as cumsum-sum.

Solution shorter for R>=4.1, by @Dominic van Essen again.

R, 97 bytes

Or R>=4.1, 83 bytes by replacing two function appearances with \s.

function(A,B=A)Map(function(i)B<<-(diag(sum(b<-diag(B))-cumsum(b))-B*upper.tri(B))%*%A,2:nrow(A))


Try it online!

• Assuming you can't have a 1x1 matrix, I think you can swap n=nrow(A) and while(n<-n-1) for for(i in 2:nrow(A)) to save a bit... try it Nov 4, 2021 at 9:24
• 96 bytes... Nov 4, 2021 at 14:04
• ...or down to 83 bytes using R>=4.1... Nov 4, 2021 at 16:16

Python 3 + numpy, 117, 99, 97, 91, 89(@Dingus), 88 bytes

from numpy import*
def f(M):
N=M;T=triu(M==M)
for _ in M:N=T@N[T&T.T,print(N)]*M-T*N@M


Try it online!

Old version

Old version

Old version

Old version

Old version

Detailed explanation

(refers to the latest and greatest versison (88 bytes))

We start with a sizeable upfront investment to create auxiliary matrix T which is just an upper triangular mask.

T is then used

• to extract the upper triangle of other matrices (obvious),
• to form another mask for extraction of the diagonal (Still rather obvious by and-ing with its transpose T&T.T) and
• to do the cumsum: matrix-vector multiplication Tv with a triangular matrix of ones yields the cumsum of v. Conveniently, we can choose between left-to-right and right-to-left summmation by using upper or lower triangles.

To avoid having to reinsert the cumsum into the diagonal we utilise the fact that matrix-matrix multiplication DM with a diagonal matrix D is the same as python style broadcasted pointwise multiplication with the diagonal reshaped into a column.

Extracting and reshaping the diagonal is done in one indexing operation N[T&T.T,print(N)]. The purpose of placing the print statement in this particular spot is that we can make good use of its return value, None. None as an index inserts a new axis which here has the effect of creating a column instead of a row.

Avoiding explicit import, 84 bytes

def f(M):
N=M;T=(M==M).cumsum(1);T=T>=T.T
for _ in M:N=T@N[T&T.T,print(N)]*M-T*N@M


Try it online!

Bypassing the call to triu the last reference to the toplevel numpy namespace has been removed; the import statement is therefore no longer required. Feels a bit ugly to me but seems to be legal.

• I think you can remove a space to save another byte: from numpy import*. Nov 5, 2021 at 5:35

J, 67 bytes

1}.(+/ .*~0-]((*0+/\"1@,~1&}.)+]*$@[$[:(-+/\.)*&,)=@i.@#)^:(<@#@])~


Try it online!

1} …^:(<@#@]) repeat n-1 times, returning the intermediate results. =@i.@# identity matrix. ]*$@[$[:(-+/\.)*&, running sum of the diagonal. (*0+/\"1@,~1&}.) everything on the right of the diagonal stays the same. +/ .*~0- matrix multiplication of the negative with the original matrix.

• 63 Nov 4, 2021 at 22:45

JavaScript (ES6), 121 bytes

Essentially the same as below, but returns an array of $$\n\$$ matrices, the first one being the input matrix.

A=>A.map((_,i)=>m=i?m.map((r,y)=>r.map((_,x)=>r.reduce((q,v,k)=>y>k?q:m.map((r,z)=>v-=y<k|z<k?0:r[z])|q-A[k][x]*v,0))):A)


Try it online!

JavaScript (ES6),  132 128  126 bytes

A naive implementation of Bird's algorithm. Returns an array of $$\n-1\$$ matrices.

A=>A.slice(1).map(_=>m=m.map((r,y)=>r.map((_,x)=>r.reduce((q,v,k)=>y>k?q:m.map((r,z)=>v-=y<k|z<k?0:r[z])|q-A[k][x]*v,0))),m=A)


Try it online!

Commented

A =>                        // A[] = input matrix of size n x n
A.slice(1).map(_ =>         // repeat n - 1 times:
m =                       //   update m[]
m.map((r, y) =>           //   for each row r[] at position y in m[]:
r.map((_, x) =>         //     for each value at position x in r[]:
r.reduce((q, v, k) => //       for each value v at position k in r[],
//       using q as an accumulator:
y > k ?             //         if y is greater than k:
q                 //           leave q unchanged
:                   //         else:
m.map((r, z) =>   //           for each row r[] at position z in m[]:
v -=            //             subtract from v:
y < k |       //               if y is less than k
z < k ?       //               or z is less than k:
0           //                 leave v unchanged
:             //               else:
r[z]        //                 subtract r[z] = m[z][z]
) |               //           end of map()
q - A[k][x] * v,  //           subtract A[k][x] * v from q
)                     //       end of reduce()
)                       //     end of map()
),                        //   end of map()
)                           // end of map()

• Interestingly, your previous answer produced different intermediate results but arrived at the same determinant! Nov 4, 2021 at 12:31
• @ViníciusMello Yes, I think I messed up the indices in the 1st version, but I'm almost certain it was eventually always giving the correct answer. (Although the new version is pretty close, it's a complete rewrite based on cleaner code.) Nov 4, 2021 at 13:32

Lua (LuaJIT), 451 392 bytes

function det(A)
local n,Y,X,y,yl,x=#A,{},{}
for i=1,n do x={} for j=1,n do x[#x+1]=A[i][j] end Y[#Y+1],X[#X+1]={},x end
for l=1,n-1 do
o(l,X)
for i=1,n do for j=1,n do Y[i][j]=0 end end
for i=n,1,-1 do for j=n,i,-1 do
y = j>i and -X[i][j] or (i==n and 0 or yl+X[i+1][i+1])
yl = i==j and y or yl
for k=1,n do Y[i][k]=Y[i][k]+y*A[j][k] end
end end
Y,X=X,Y
end
o(n,X)
end


Try it online!

Old version

• Welcome to Code Golf, and nice first answer! Make sure to check out our tips for golfing in Lua to see if there are any ways for you to golf your code. Nov 4, 2021 at 13:42
• For reference only. I did a lot of golfing in order to build matrix $\mu(X)$ bottom up. Nov 4, 2021 at 13:44
• Oh, whoops, I didn't see that you were the OP. Nov 4, 2021 at 13:48

APL (Dyalog Extended), 62 38 bytes

{-⍺+.×⍨(+\-+/)@(∘.=⍨⍳∘≢)⍵×∘.≤⍨⍳≢⍵}\≢⍴⊂


Try it online!

Breaking from right to left:

1. ≢⍴⊂ : $$\[A,A,\cdots,A]\$$ (produces $$\n\$$ copies of input);
2. {...}\ : $$\[A,A⊗A,A⊗(A⊗A),\cdots]\$$ ({...} implements $$\⊗\$$);
3. ⍵×∘.≤⍨⍳≢⍵: upper triangular matrix;
4. (+\-+/)@(∘.=⍨⍳∘≢) : replace diagonal by the negated sum. Thus, (+\-+/)@(∘.=⍨⍳∘≢)⍵×∘.≤⍨⍳≢⍵ is $$\\mu(X)\$$. (+\-+/) is a very nice trick to compute the new diagonal, using a train: (+\diag)-(+/diag);
5. -⍺+.×⍨...: $$\-\mu(X)\cdot A=A⊗X\$$.

Many thanks to @rak1507, @Adám and @Marshall for this nice improvement!

{((⌽0,+\¯1↓⌽1 1⍉⍵)(×⍤0 1)⍺)-((⍵×(∘.<)⍨(⍳≢⍵))+.×⍺)}\((≢⊢)⍴(⊂⊢))


Old version

• I started this challenge after watching the wonderful APL videos from the code_report channel on youyube, so I was hoping to see some solution in Dyalog APL. As the answer hasn't appeared yet, my attempt follows. Nov 7, 2021 at 23:59
• Quick golf on your solution gives 43 bytes. And there was some golfing in the APL orchard, giving a 38-byte solution. Nov 8, 2021 at 8:18
• Thanks, @Bubbler! I have so much to learn! Nov 8, 2021 at 10:38

Charcoal, 47 bytes

≔θηＦ⊖Ｌθ«≔ＥηＥκΣＥθ×⎇⁼πλΣＥη∧›ςλ§ρς∧›πλ±§κπ§ξνη⟦⭆¹η


Try it online! Link is to verbose version of code. Explanation:

≔θη


Start with the output X equal to the input.

Ｆ⊖Ｌθ«


Repeat n-1 times.

≔ＥηＥκΣＥθ×...§ξνη


Matrix multiply -μ(X) by A.

⎇⁼πλΣＥη∧›ςλ§ρς∧›πλ±§κπ


Dynamically generate -μ(X) as part of the multiplication; for elements on the main diagonal, sum the elements further down on the diagonal, otherwise negate or zero the element depending on whether it is above or below the diagonal.

⟦⭆¹η


Pretty print the latest result. (Normal printing would save 1 byte but be almost impossible to read.)

Pari/GP, 80 bytes

a->[x=matrix(n,,i,j,if(i>j,0,i-j,-x[i,j],sum(k=i+1,n,x[k,k])))*a|l<-[2..n=#x=a]]


Try it online!

• With this little change of the minus sign, from the definition of $⊗$ to the definition of $\mu(X)$, the program becomes clearer. Your answer is succinct, precise and literal at the same time! Nov 5, 2021 at 11:10