A division-free determinant algorithm

This file defines birdDetand Spec.birdDet, implementations of an division-free algorithm for computing determinants. The algorithm runs in O(n^4) for an n-by-n matrix.

This determinant algorithm comes from [Richard S. Bird, A simple division-free algorithm for computing determinants][bird2011].

Main definitions

• BirdDet.birdDet: The entrypoint for the determinant calculation.
• BirdDet.get: matrix entry lookup.
• BirdDet.sumFrom: The sum f lo + ... + f (n - 1).
• BirdDet.stepEntry: One scalar recurrence step.
• BirdDet.Spec.birdDet: An implementation of Bird's algorithm using Matrix.

Main lemmas

The lemmas in this file are unfolding equations.

def BirdDet.get {R : Type u_1} [CommRing R] (n : ℕ) (A : Array R) (i j : ℕ) : R

get n A i j returns the (i, j)th entry of the n × n matrix whose entries are stored in A in row-major order.

The function does not check the matrix index bounds.

@[irreducible]

def BirdDet.sumFrom {R : Type u_1} [CommRing R] (n lo : ℕ) (f : ℕ → R) : R

Sum f lo + ... + f (n - 1). Returns zero when n <= lo.

def BirdDet.stepEntry {R : Type u_1} [CommRing R] (n : ℕ) (A : Array R) (F : ℕ → ℕ → R) (i j : ℕ) : R

One entry of one scalar Bird recurrence step.

Bird's paper defines a matrix recursion for an n × n matrix A:

F_0 = A
F_{t+1} = μ(F_t) * A

where μ(F_t) is obtained from F_t by replacing each diagonal entry F_t k k with the negative sum of the diagonal entries below it, setting the entries in the lower triangular part to 0, and leaving all other entries unchanged:

μ(F_t) =
  0                                   if i >= j
  - ∑ k from i+1 to n-1, F_t k k      if i = j
  F_t i j                             if i < j

If we write out the entry-wise matrix multiplication F_{t+1} i j = (μ(F_t) * A) i j we obtain:

F_{t+1} i j =
  - (∑ k from i+1 to n-1, F_t k k) * (A i j)
  + ∑ k from i+1 to n-1, (F_t i k) * (A k j)

def BirdDet.birdDet {R : Type u_1} [CommRing R] (n : ℕ) (A : Array R) : R

birdDet n A computes the determinant of the n × n matrix whose entries are stored in A in row-major order.

theorem BirdDet.get_eq {R : Type u_1} [CommRing R] (n : ℕ) (A : Array R) (i j : ℕ) : BirdDet.get n A i j = A.getD (n * i + j) 0

Unfold a row-major matrix entry lookup.

theorem BirdDet.sumFrom_step {R : Type u_1} [CommRing R] (n lo : ℕ) (f : ℕ → R) (h : lo < n) : BirdDet.sumFrom n lo f = f lo + BirdDet.sumFrom n (lo + 1) f

theorem BirdDet.sumFrom_stop {R : Type u_1} [CommRing R] (n lo : ℕ) (f : ℕ → R) (h : ¬lo < n) : BirdDet.sumFrom n lo f = 0

theorem BirdDet.sumFrom_induct (n : ℕ) (motive : ℕ → Prop) (step : ∀ lo < n, motive (lo + 1) → motive lo) (stop : ∀ (lo : ℕ), ¬lo < n → motive lo) (lo : ℕ) : motive lo

Induction following the recursive structure of sumFrom.

theorem BirdDet.stepEntry_eq {R : Type u_1} [CommRing R] (n : ℕ) (A : Array R) (F : ℕ → ℕ → R) (i j : ℕ) : stepEntry n A F i j = (-BirdDet.sumFrom n (i + 1) fun (k : ℕ) => F k k) * BirdDet.get n A i j + BirdDet.sumFrom n (i + 1) fun (k : ℕ) => F i k * BirdDet.get n A k j

Unfold one scalar Bird recurrence step to the entry-wise formula.

theorem BirdDet.birdDet_zero {R : Type u_1} [CommRing R] (A : Array R) : birdDet 0 A = 1

theorem BirdDet.birdDet_succ {R : Type u_1} [CommRing R] (k : ℕ) (A : Array R) : birdDet (k + 1) A = (-1) ^ k * (stepEntry (k + 1) A)^[k] (BirdDet.get (k + 1) A) 0 0

Unfold birdDet at a successor dimension.

theorem BirdDet.birdDet_eq {R : Type u_1} [CommRing R] (n k : ℕ) (A : Array R) (hn : n = k + 1) : birdDet n A = (-1) ^ k * (stepEntry n A)^[k] (BirdDet.get n A) 0 0

def BirdDet.Spec.stepEntry {R : Type u_1} [CommRing R] {n : ℕ} (A F : Matrix (Fin n) (Fin n) R) : Matrix (Fin n) (Fin n) R

One entry of one Matrix/Fin Bird recurrence step.

def BirdDet.Spec.birdDet {R : Type u_1} [CommRing R] {n : ℕ} (A : Matrix (Fin n) (Fin n) R) : R

A version of the Bird determinant algorithm that is stated in terms of Matrix.

theorem BirdDet.Spec.stepEntry_eq {R : Type u_1} [CommRing R] {n : ℕ} (A F : Matrix (Fin n) (Fin n) R) : stepEntry A F = Matrix.of fun (i j : Fin n) => (-∑ k > i, F k k) * A i j + ∑ k > i, F i k * A k j

@[simp]

theorem BirdDet.Spec.birdDetSpec_zero {R : Type u_1} [CommRing R] (A : Matrix (Fin 0) (Fin 0) R) : birdDet A = 1

theorem BirdDet.Spec.birdDetSpec_succ {R : Type u_1} [CommRing R] {k : ℕ} (A : Matrix (Fin (k + 1)) (Fin (k + 1)) R) : birdDet A = (-1) ^ k * (stepEntry A)^[k] A 0 0