Bird's Algorithm for Computing Determinants

Question

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]]

Answer 1

K (ngn/k), ⁴³ 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

Answer 2

R, 113 109 103 96 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!

Answer 3

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.

Answer 4

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.

Answer 5

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
        0                   //         start with q = 0
      )                     //       end of reduce()
    )                       //     end of map()
  ),                        //   end of map()
  m = A                     //   start with m[] = A[]
)                           // end of map()

Answer 6

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

Answer 7

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

Answer 8

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.)

Answer 9

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!