num ← ##.det nmat                           ⍝ Determinant of square matrix.

From Roger Hui,  [det] uses Gaussian elimination & Laplace expansion to find the
determinant of square matrix [nmat].

[det] finds the maximal entry in the entire matrix, not just the  maximal  entry
in column 0 (i.e. it does "full pivoting".).  It doesn't do row or column inter-
changes at all.  Instead,  after  Gaussian  elimination to zero out the i-th row
(other than column j) or zero out the j-th column (other than row i),  it does a
Laplace expansion on the i-th row (or the j-th column).  Post-Gaussian eliminat-
ion, either expansion has just one non-zero  co-factor, that is,  just one minor
∇ ⍵[k~i;k~j] with a non-zero coefficient ⍵[i;j]ׯ1*i+j,  and that coefficient is
multiplied by ⍺ on the tail call.

Alternative codings
-------------------
Also from Roger, an earlier draft of [det]:

    det_4a←{⎕io←0
        ⍺←1                     ⍝ product of diagonal entries so far
        0=n←⊃⍴⍵:⍺               ⍝ answer for 0-by-0
        p←{⍵⍳⌈/⍵}|⍵[;0]         ⍝ index of pivot row
        p=n:0                   ⍝ ⍵ is singular if column is all 0
        k←⍳n                    ⍝ row indices
        k[⌽i]←k[i←(×p)/p,0]     ⍝ interchange p and 0, if different
        (⍺×dׯ1*×p)∇ ⍵[1↓k;1↓⍳n]-(⍵[1↓k;0]÷d←⍵[⊃k;0])∘.×⍵[⊃k;1↓⍳n]
    }

The  following  rather  unpleasant  one-line  alternative,  is included just for
amusement. It is interesting only in that it:

- Is both ⎕IO and ⎕ML independent.

- Has  no  guards  or  local variables. It is a single pure APL expression that,
  without testing or looping, denotes the determinant of its argument matrix.

- Codes  the mathematical definition of determinant, which is an O(!n) algorithm
  and so is viable for only the smallest of argument matrices.

- Is the sort of thing that gets APL a bad name.

    detriment←{∇{⍵+(¯1*+/∧\⍺)×⍵⍵[⎕io;⍺⍳0]×⍺⍺ ⍺/1 0↓⍵⍵}⍵/(↓≠/¨⍳⍴⍵),0 0≡⍴⍵}

VMJ suggests this version using the Gaussian method:

    detG←{                                  ⍝ Gaussian determinant.
        (⎕io ⎕ml)←1 3 ⋄ r←⍴⍵
        1≥×/r:↑⍵ ⋄ 2≠⍴r:0 ⋄ ≠/r:0           ⍝ check for trivial cases
        1{                                  ⍝ inner loop
            2>↑⍴⍵:↑⍺×⍵                      ⍝ end? -> result
            (m a c)←{                       ⍝ calculate matrix, anchor & coeff
                0≠1⍴⍵:⍵(1⍴⍵)1               ⍝ default: original matrix
                z←⍵ ⋄ j←{⍵⍳⌈/⍵}|z[;1]       ⍝ look for 1st non-zero
                z[1,j;]←z[j,1;]             ⍝ reorder matrix
                z(1⍴z)¯1                    ⍝ NB. coeff=¯1
            }⍵
            a=0:0                           ⍝ row of zeroes..
            (⍺×a×c)∇ 1 1↓m-m[;1]∘.×m[1;]÷a  ⍝ re-curses!
        }⍵
    }

and ...

    det_mv←{        ⍝ Matrix determinant, based on Mahajan-Vinay algorithm:

        ⍝|  Meena Mahajan and V Vinay.
        ⍝|   Determinant: Combinatorics, Algorithms, and Complexity
        ⍝|    https://www.imsc.res.in/~meena/publ.html

        (⎕io ⎕ml)←0 3 ⋄ mat←⍵ ⋄ n←↑⍴mat
        V←(2,3/n)⍴0 ⋄ V[2|n;;;0]←1

        ↑{i←⍵
            ↑{v←⍵
                ↑{u←⍵ ⋄ n=⍵+1:0
                    ↑{
                        V[0;u;⍵;i+1]+←V[0;u;v;i]×mat[v;⍵]
                        V[1;⍵;⍵;i+1]+←V[0;u;v;i]×mat[v;u]
                        V[1;u;⍵;i+1]+←V[1;u;v;i]×mat[v;⍵]
                        V[0;⍵;⍵;i+1]+←V[1;u;v;i]×mat[v;u]
                        0
                    }¨(⍵+1)↓⍳n
                }¨⍳⍵+1
            }¨⍳n
        }¨⍳n-1:

        +/{-+/mat[⍵;⍳⍵+1]×-⌿V[;⍳⍵+1;⍵;n-1]}¨⍳n
    }

Examples:

    det 2 2⍴2 3 4 5
¯2
    det 2 2⍴0 1 1 0
¯1
    det 1 1⍴3
3
    det 0 0⍴7
1
    hil←{÷1+∘.+⍨(⍳⍵)-⎕io}   ⍝ Order ⍵ Hilbert matrix.

    det hil 10
2.1644805E¯53

See also: gauss_jordan Cholesky

Back to: contents

Back to: Workspaces