|
| 1 | +/- |
| 2 | +Copyright (c) 2026 Paul Cadman. All rights reserved. |
| 3 | +Released under Apache 2.0 license as described in the file LICENSE. |
| 4 | +Authors: Paul Cadman |
| 5 | +-/ |
| 6 | +module |
| 7 | + |
| 8 | +public import Mathlib.LinearAlgebra.Matrix.Determinant.Basic |
| 9 | +public import Mathlib.LinearAlgebra.Matrix.Determinant.Bird.Defs |
| 10 | +import Mathlib.Algebra.Order.BigOperators.Group.LocallyFinite |
| 11 | +import Mathlib.Data.Fintype.Order |
| 12 | +import Mathlib.Order.Preorder.Finite |
| 13 | + |
| 14 | +/-! |
| 15 | +# Correctness of Bird's determinant algorithm |
| 16 | +
|
| 17 | +This file contains a proof that Bird's division-free algorithm computes |
| 18 | +`Matrix.det`, in both its matrix form `BirdDet.Spec.birdDet` |
| 19 | +(`birdDetSpec_eq_det`) and its flat-array form `BirdDet.birdDet` |
| 20 | +(`det_eq_birdDet`), formalizing the combinatorial argument of |
| 21 | +[Richard S. Bird, *A simple division-free algorithm for computing determinants*][bird2011]. |
| 22 | +
|
| 23 | +## Correspondence with the paper |
| 24 | +
|
| 25 | +* A word of length `p` is a tuple `Fin p → Fin n`, using the indexing convention |
| 26 | + above (NB: Indices in Bird's paper start at 1). |
| 27 | +* `f[α, β]`, the minor on rows `α` and columns `β`, is `(A.submatrix α β).det`. |
| 28 | +* `f[iα, jα]`, a bordered minor, is `bminor A i j α`, with the word `iα` spelled |
| 29 | + `Fin.cons i α`. |
| 30 | +* `f[α, α]`, a principal minor, is `pminor A α`. |
| 31 | +* If `i : Fin n` represents Bird's symbol `r = i.val + 1`, then Bird's |
| 32 | + `βᵣ = [r + 1, ..., n]` is represented by `Finset.Ioi i`. |
| 33 | +* Bird's `Sₚ(βᵣ)`, the length `p` subsequences of `βᵣ`, is represented by `S p i`. |
| 34 | +
|
| 35 | +The theorem names `paper_eq1`, ..., `paper_eq5` follow Bird's numbering. |
| 36 | +
|
| 37 | +## Main results |
| 38 | +
|
| 39 | +- `BirdDet.birdDetSpec_eq_det`: `Matrix.det` computes the same determinant as `BirdDet.Spec.birdDet` |
| 40 | +- `BirdDet.det_eq_birdDet`: `Matrix.det` computes the same determinant as `BirdDet.birdDet` |
| 41 | +-/ |
| 42 | + |
| 43 | +namespace BirdDet |
| 44 | + |
| 45 | +open Function |
| 46 | +variable {R : Type*} [CommRing R] {m n : ℕ} |
| 47 | + |
| 48 | +/-- `sumFrom n lo f` is the sum of `f` over the half-open interval `[lo, n)`. -/ |
| 49 | +theorem sumFrom_eq_sum_Ico {lo : ℕ} (f : ℕ → R) : |
| 50 | + BirdDet.sumFrom n lo f = ∑ k ∈ Finset.Ico lo n, f k := by |
| 51 | + induction lo using BirdDet.sumFrom_induct n with |
| 52 | + | step lo hlo ih => rw [sumFrom_step n lo f hlo, ih, ← Finset.sum_eq_sum_Ico_succ_bot hlo f] |
| 53 | + | stop lo hlo => rw [sumFrom_stop n lo f hlo, Finset.Ico_eq_empty hlo, Finset.sum_empty] |
| 54 | + |
| 55 | +theorem sumFrom_fin_tail (i : Fin n) (f : ℕ → R) : |
| 56 | + BirdDet.sumFrom n (i.val + 1) f = ∑ k ∈ Finset.Ioi i, f k.val := calc |
| 57 | + _ = ∑ k ∈ Finset.Ico (i.val + 1) n, f k := by rw [sumFrom_eq_sum_Ico] |
| 58 | + _ = ∑ k ∈ (Finset.range n).filter (i.val < ·), f k := by congr; ext; aesop |
| 59 | + _ = ∑ k ∈ Finset.range n, if i.val < k then f k else 0 := by rw [Finset.sum_filter] |
| 60 | + _ = ∑ k : Fin n, if i.val < k.val then f k.val else 0 := by rw [← Fin.sum_univ_eq_sum_range] |
| 61 | + _ = ∑ k ∈ Finset.Ioi i, f k.val := by simp [← Finset.sum_filter, Finset.filter_lt_eq_Ioi] |
| 62 | + |
| 63 | +/-- The scalar recurrence initialized by array lookup agrees pointwise |
| 64 | + with the matrix recurrence. -/ |
| 65 | +theorem iterate_stepEntry_get_eq_spec (A : Array R) (hA : A.size = n * n) (t : ℕ) (i j : Fin n) : |
| 66 | + ((stepEntry n A)^[t] (BirdDet.get n A)) i.val j.val = |
| 67 | + (Spec.stepEntry (.ofArray A hA))^[t] (.ofArray A hA) i j := by |
| 68 | + rw [Matrix.ofArray_eq_of_getD] |
| 69 | + induction t generalizing i j with |
| 70 | + | zero => simp [get_eq] |
| 71 | + | succ t ih => |
| 72 | + simp_rw [iterate_succ_apply', stepEntry_eq, Spec.stepEntry_eq, sumFrom_fin_tail, ih, |
| 73 | + Matrix.of_apply, get_eq] |
| 74 | + |
| 75 | +/-- The flat-array algorithm `BirdDet.birdDet` computes the same determinant as |
| 76 | + `BirdDet.Spec.birdDet`. -/ |
| 77 | +theorem birdDet_eq_birdDetSpec (A : Array R) (hA : A.size = n * n) : |
| 78 | + birdDet n A = Spec.birdDet (.ofArray A hA) := by |
| 79 | + cases n with |
| 80 | + | zero => rw [birdDet_zero, Spec.birdDetSpec_zero] |
| 81 | + | succ k => simp [birdDet_succ, Spec.birdDetSpec_succ, ← iterate_stepEntry_get_eq_spec A hA k] |
| 82 | + |
| 83 | +variable (A : Matrix (Fin n) (Fin n) R) {p : ℕ} |
| 84 | + |
| 85 | +/-- Bird's bordered minor `f[iα, jα]`. -/ |
| 86 | +abbrev bminor (i j : Fin n) (α : Fin p → Fin n) : R := |
| 87 | + (A.submatrix (Fin.cons i α) (Fin.cons j α)).det |
| 88 | + |
| 89 | +/-- Bird's principal minor `f[α, α]`. -/ |
| 90 | +abbrev pminor (α : Fin p → Fin n) : R := |
| 91 | + (A.submatrix α α).det |
| 92 | + |
| 93 | +lemma det_submatrix_removeNth_eq_sign_mul_bminor |
| 94 | + (α : Fin (p + 1) → Fin n) (i : Fin n) (s : Fin (p + 1)) : |
| 95 | + (A.submatrix (Fin.cons i (s.removeNth α)) α).det = |
| 96 | + (-1 : R) ^ s.val * bminor A i (α s) (s.removeNth α) := |
| 97 | + calc |
| 98 | + _ = (-1 : R) ^ s.val * ((A.submatrix (Fin.cons i (s.removeNth α)) α) |
| 99 | + |>.submatrix id (Fin.cycleRange s).symm).det := by |
| 100 | + rw [Matrix.det_permute'] |
| 101 | + simp [← mul_assoc, ← pow_add] |
| 102 | + _ = (-1 : R) ^ s.val * bminor A i (α s) (s.removeNth α) := by |
| 103 | + congrm _ * Matrix.det ?_ |
| 104 | + simp [Fin.cons_removeNth_eq_comp_cycleRange_symm] |
| 105 | + |
| 106 | +/-- First-column Laplace expansion of a bordered minor -/ |
| 107 | +theorem det_bordered_expand (α : Fin (p + 1) → Fin n) (i j : Fin n) : |
| 108 | + bminor A i j α = |
| 109 | + pminor A α * A i j - ∑ s : Fin (p + 1), bminor A i (α s) (s.removeNth α) * A (α s) j := calc |
| 110 | + _ = A i j * pminor A α + |
| 111 | + ∑ s : Fin (p + 1), (-1 : R) ^ (s.val + 1) * A (α s) j * |
| 112 | + (A.submatrix (Fin.cons i (s.removeNth α)) α).det := by |
| 113 | + rw [bminor, Matrix.det_succ_column_zero, Fin.sum_univ_succ]; simp |
| 114 | + _ = pminor A α * A i j + |
| 115 | + ∑ s : Fin (p + 1), |
| 116 | + ((-1 : R) ^ (s.val + 1) * |
| 117 | + (A.submatrix (Fin.cons i (s.removeNth α)) α).det) * A (α s) j := by |
| 118 | + simp only [mul_comm (A i j), mul_right_comm] |
| 119 | + _ = pminor A α * A i j + |
| 120 | + ∑ s : Fin (p + 1), -(bminor A i (α s) (s.removeNth α) * A (α s) j) := by |
| 121 | + simp only [det_submatrix_removeNth_eq_sign_mul_bminor, ← mul_assoc, ← pow_add]; aesop |
| 122 | + _ = pminor A α * A i j - |
| 123 | + ∑ s : Fin (p + 1), bminor A i (α s) (s.removeNth α) * A (α s) j := by |
| 124 | + simp only [Finset.sum_neg_distrib, sub_eq_add_neg] |
| 125 | + |
| 126 | +/-- A bordered minor is zero when its border column already occurs in the word. -/ |
| 127 | +theorem bminor_eq_zero_of_mem_range |
| 128 | + {k : Fin n} (α : Fin p → Fin n) (i : Fin n) (hk : k ∈ Set.range α) : |
| 129 | + bminor A i k α = 0 := by |
| 130 | + obtain ⟨q, rfl⟩ := hk |
| 131 | + -- The repeated columns in the submatrix used in bminor are `0` and `q + 1`. |
| 132 | + exact Matrix.det_zero_of_column_eq q.succ_ne_zero <| by simp |
| 133 | + |
| 134 | +/-- `S p i` is Bird's `Sₚ(βᵢ)`: words `α` of length `p` over the alphabet `βᵢ`. -/ |
| 135 | +def S (p : ℕ) (i : Fin n) : Finset (Fin p → Fin n) := |
| 136 | + {α : Fin p → Fin n | StrictMono (Fin.cons i α)} |
| 137 | + |
| 138 | +/-- Membership in `S p i` is strict monotonicity of the bordered word. -/ |
| 139 | +theorem mem_S_iff {p : ℕ} {i : Fin n} {α : Fin p → Fin n} : |
| 140 | + α ∈ S p i ↔ StrictMono (Fin.cons i α) := |
| 141 | + Finset.mem_filter_univ α |
| 142 | + |
| 143 | +/-- The base case of equation (1): `S₀(α) = {ε}`, the singleton of the empty |
| 144 | +word `ε`. -/ |
| 145 | +theorem S_zero (i : Fin n) : S 0 i = {![]} := by |
| 146 | + ext; simp [mem_S_iff, Fin.strictMono_iff_lt_succ, eq_iff_true_of_subsingleton] |
| 147 | + |
| 148 | +/-- The unique maximum-length word over the symbols above `0` is `Fin.succ`. -/ |
| 149 | +@[simp] lemma S_zero_eq_singleton {p : ℕ} : S p 0 = {Fin.succ} := by |
| 150 | + ext; simp [mem_S_iff] |
| 151 | + |
| 152 | +/-! ## Decomposition `S_{p+1}(βᵢ) = { kα | k ∈ βᵢ, α ∈ S_p(β_k) }` -/ |
| 153 | + |
| 154 | +/-- `S (p + 1) i` can be written as the image of a `biUnion` -/ |
| 155 | +theorem S_succ_eq_biUnion {p : ℕ} (i : Fin n) : |
| 156 | + S (p + 1) i = (Finset.Ioi i).biUnion fun k ↦ (S p k).image (Fin.cons k) := by |
| 157 | + ext α |
| 158 | + simp only [Finset.mem_biUnion, Finset.mem_image, Finset.mem_Ioi, mem_S_iff] |
| 159 | + refine ⟨fun hα ↦ ⟨α 0, (Fin.strictMono_cons.mp hα).1 0, Fin.tail α, ?_, Fin.cons_self_tail α⟩, ?_⟩ |
| 160 | + · simp only [Fin.cons_self_tail] |
| 161 | + exact hα.comp Fin.strictMono_succ |
| 162 | + · rintro ⟨k, hk, u, hu, rfl⟩ |
| 163 | + exact StrictMono.vecCons hu hk |
| 164 | + |
| 165 | +/-- A symbol above `i` that does not occur in a word in `S p i` can be inserted while |
| 166 | +preserving strict monotonicity. -/ |
| 167 | +lemma exists_insertNth_mem_S {p : ℕ} {i : Fin n} {α : Fin p → Fin n} {k : Fin n} |
| 168 | + (hα : α ∈ S p i) (hik : i < k) (hk : k ∉ Set.range α) : |
| 169 | + ∃ t : Fin (p + 1), t.insertNth k α ∈ S (p + 1) i := by |
| 170 | + set t := ⨅ j ∈ {j | k < α j}, j.castSucc with t_eq |
| 171 | + use t |
| 172 | + simp only [mem_S_iff, Fin.strictMono_cons] at ⊢ hα |
| 173 | + refine ⟨fun j ↦ Fin.succAboveCases t ?_ ?_ j, ?_⟩ |
| 174 | + · simp only [t_eq, Set.mem_ofPred_eq, Fin.strictMono_insertNth_iff, hα.2, lt_iInf_iff, |
| 175 | + le_iInf_iff, forall_exists_index, and_imp, iInf_le_iff_forall_lt, iInf_lt_iff, |
| 176 | + exists_prop, true_and] |
| 177 | + refine ⟨fun j x hjx h ↦ ?_, fun j h ↦ ?_⟩ |
| 178 | + · contrapose! h |
| 179 | + have k_ne (j : Fin p) : k ≠ α j := fun hj ↦ hk ⟨j, hj.symm⟩ |
| 180 | + exact ⟨j, h.lt_of_ne (k_ne _), hjx⟩ |
| 181 | + · obtain ⟨q, hkq, hqj⟩ := h j.succ j.castSucc_lt_succ |
| 182 | + exact hkq.trans_le <| hα.2.monotone <| Fin.castSucc_lt_succ_iff.mp hqj |
| 183 | + · simpa |
| 184 | + · simpa using hα.1 |
| 185 | + |
| 186 | +variable (p) in |
| 187 | +/-- Bird's equation (1) : `x^(p)_ij = (-1)^p ∑ { f[iα, jα] | α ∈ S_p(βᵢ) }`. -/ |
| 188 | +abbrev Eq1 : Prop := |
| 189 | + (Spec.stepEntry A)^[p] A = .of fun i j ↦ (-1) ^ p * ∑ α ∈ S p i, bminor A i j α |
| 190 | + |
| 191 | +/-! ## Equations (2) and (3): substituting the induction hypothesis -/ |
| 192 | + |
| 193 | +/-- Bird's equation (2), assuming equation (1) at `p` as the induction hypothesis. -/ |
| 194 | +theorem paper_eq2 (i : Fin n) (hEq1 : Eq1 A p) : |
| 195 | + (-∑ k ∈ Finset.Ioi i, (Spec.stepEntry A)^[p] A k k) = |
| 196 | + (-1) ^ (p + 1) * ∑ α ∈ S (p + 1) i, pminor A α := by |
| 197 | + calc |
| 198 | + (-∑ k ∈ Finset.Ioi i, (Spec.stepEntry A)^[p] A k k) = |
| 199 | + (-1) ^ (p + 1) * ∑ k ∈ Finset.Ioi i, ∑ α ∈ S p k, bminor A k k α := by |
| 200 | + simp only [hEq1, Matrix.of_apply, ← Finset.mul_sum] |
| 201 | + ring |
| 202 | + _ = (-1) ^ (p + 1) * ∑ α ∈ S (p + 1) i, pminor A α := by |
| 203 | + rw [S_succ_eq_biUnion, Finset.sum_biUnion] |
| 204 | + · congrm (((-1) ^ (p + 1) * ∑ k ∈ Finset.Ioi i, ?_)) |
| 205 | + symm |
| 206 | + exact Finset.sum_image fun _ _ _ _ hαβ => (Fin.cons_inj.mp hαβ).2 |
| 207 | + · grind [Set.PairwiseDisjoint, Set.Pairwise, Finset.disjoint_left, Fin.cons_inj] |
| 208 | + |
| 209 | +/-- Bird's equation (3), assuming equation (1) at `p` as the induction hypothesis. -/ |
| 210 | +theorem paper_eq3 (i j : Fin n) (hEq1 : Eq1 A p) : |
| 211 | + ((Spec.stepEntry A)^[p + 1] A) i j = |
| 212 | + (-1) ^ (p + 1) * (∑ α ∈ S (p + 1) i, pminor A α * A i j - |
| 213 | + ∑ k ∈ Finset.Ioi i, ∑ α ∈ S p i, bminor A i k α * A k j) := by |
| 214 | + simp_rw [iterate_succ_apply', Spec.stepEntry_eq, Matrix.of_apply, paper_eq2 _ _ hEq1, hEq1, |
| 215 | + Matrix.of_apply, mul_assoc, Finset.sum_mul, ← Finset.mul_sum] |
| 216 | + ring |
| 217 | + |
| 218 | +/-! ## Equation (5): first-column Laplace expansion -/ |
| 219 | + |
| 220 | +/-- Bird's equation (5) -/ |
| 221 | +theorem paper_eq5 (i j : Fin n) : |
| 222 | + ∑ α ∈ S (p + 1) i, bminor A i j α = |
| 223 | + ∑ α ∈ S (p + 1) i, pminor A α * A i j - |
| 224 | + ∑ α ∈ S (p + 1) i, ∑ t : Fin (p + 1), bminor A i (α t) (t.removeNth α) * A (α t) j := calc |
| 225 | + _ = ∑ α ∈ S (p + 1) i, (pminor A α * A i j - |
| 226 | + ∑ t : Fin (p + 1), bminor A i (α t) (t.removeNth α) * A (α t) j) := by |
| 227 | + exact Finset.sum_congr rfl <| by simp [det_bordered_expand] |
| 228 | + _ = ∑ α ∈ S (p + 1) i, pminor A α * A i j - |
| 229 | + ∑ α ∈ S (p + 1) i, ∑ t : Fin (p + 1), bminor A i (α t) (t.removeNth α) * A (α t) j := by |
| 230 | + rw [Finset.sum_sub_distrib] |
| 231 | + |
| 232 | +/-! ## Comparing equations (3) and (5): reindex by sorted insert/delete -/ |
| 233 | + |
| 234 | +/-- The off-diagonal sums in Bird's equations (3) and (5) agree. -/ |
| 235 | +theorem paper_eq3_eq5_off_diag (i j : Fin n) : |
| 236 | + ∑ k ∈ Finset.Ioi i, ∑ α ∈ S p i, bminor A i k α * A k j = |
| 237 | + ∑ α' ∈ S (p + 1) i, ∑ t : Fin (p + 1), bminor A i (α' t) (t.removeNth α') * A (α' t) j := by |
| 238 | + rw [Finset.sum_comm, ← Finset.sum_product', ← Finset.sum_product'] |
| 239 | + -- The right-hand summand is the left-hand summand composed with the deletion map |
| 240 | + -- |
| 241 | + -- d (α, t) := (t.removeNth α, α t). |
| 242 | + -- |
| 243 | + -- This map is injective, and every left-hand summand outside its image is zero, |
| 244 | + -- so `sum_of_injOn` applies. |
| 245 | + symm |
| 246 | + refine Finset.sum_of_injOn (fun ⟨α, k⟩ ↦ ⟨k.removeNth α, α k⟩) ?_ ?_ ?_ ?_ |
| 247 | + · simp only [Set.InjOn, Finset.coe_product, Finset.coe_univ, Set.mem_prod, Set.mem_univ, |
| 248 | + and_true, Finset.mem_coe, Prod.mk.injEq, and_imp, Prod.forall, mem_S_iff, Fin.strictMono_cons] |
| 249 | + intros α k hi hiα α' k' hj hiα' hremove hvalue |
| 250 | + suffices hrange : Set.range α = Set.range α' by |
| 251 | + rw [hiα.range_inj hiα'] at hrange |
| 252 | + subst α' |
| 253 | + exact ⟨rfl, hiα.injective hvalue⟩ |
| 254 | + calc |
| 255 | + _ = Set.insert (α k) (Set.range (k.removeNth α)) := by |
| 256 | + rw [← Fin.range_insertNth, Fin.insertNth_self_removeNth] |
| 257 | + _ = Set.insert (α' k') (Set.range (k'.removeNth α')) := by |
| 258 | + rw [hvalue, hremove] |
| 259 | + _ = Set.range α' := by |
| 260 | + rw [← Fin.range_insertNth, Fin.insertNth_self_removeNth] |
| 261 | + · rintro ⟨α, t⟩ hα |
| 262 | + simp only [Finset.coe_product, Finset.coe_univ, Set.mem_prod, Set.mem_univ, and_true, |
| 263 | + Finset.mem_coe, Finset.coe_Ioi, Set.mem_Ioi, mem_S_iff, Fin.strictMono_cons] at hα ⊢ |
| 264 | + obtain ⟨hbound, hmono⟩ := hα |
| 265 | + exact ⟨⟨fun q => hbound (t.succAbove q), hmono.removeNth t⟩, hbound t⟩ |
| 266 | + · rintro ⟨α, k⟩ htarget hnotmem |
| 267 | + simp only [Finset.mem_product, Finset.mem_Ioi] at htarget |
| 268 | + obtain ⟨hα, hk⟩ := htarget |
| 269 | + by_cases hoccurs : k ∈ Set.range α |
| 270 | + · -- The border column `k` is repeated among the columns indexed by `α` and |
| 271 | + -- so the bordered minor is 0. |
| 272 | + rw [bminor_eq_zero_of_mem_range A α i hoccurs, zero_mul] |
| 273 | + · contrapose hnotmem |
| 274 | + obtain ⟨t, ht⟩ := exists_insertNth_mem_S hα hk hoccurs |
| 275 | + exact ⟨(t.insertNth k α, t), by simpa, by simp⟩ |
| 276 | + · simp |
| 277 | + |
| 278 | +/-! ## Bird's Equation (1) -/ |
| 279 | +theorem paper_eq1 : Eq1 A p := by |
| 280 | + induction p with |
| 281 | + | zero => |
| 282 | + ext i j |
| 283 | + simp [iterate_zero_apply, S_zero, bminor] |
| 284 | + | succ p ih => |
| 285 | + ext i j |
| 286 | + rw [Matrix.of_apply, paper_eq3 A i j ih, paper_eq5 A, paper_eq3_eq5_off_diag A] |
| 287 | + |
| 288 | +/-! ## instantiating equation (1) to prove Theorem 1 -/ |
| 289 | + |
| 290 | +/-- Bird's Theorem 1 -/ |
| 291 | +theorem birdDetSpec_eq_det (A : Matrix (Fin n) (Fin n) R) : |
| 292 | + Matrix.det A = Spec.birdDet A := by |
| 293 | + cases n with |
| 294 | + | zero => simp |
| 295 | + | succ k => |
| 296 | + have : ∑ α ∈ S k 0, bminor A 0 0 α = A.det := by simp [bminor] |
| 297 | + rw [Spec.birdDetSpec_succ, paper_eq1, Matrix.of_apply, ← mul_assoc, ← pow_add]; aesop |
| 298 | + |
| 299 | +/-- `BirdDet.birdDet n A` computes the determinant of the `n × n` matrix whose |
| 300 | + entries are stored in row-major order in `A`. -/ |
| 301 | +public theorem det_eq_birdDet (A : Array R) (hA : A.size = n * n) : |
| 302 | + Matrix.det (.ofArray A hA) = birdDet n A := by |
| 303 | + rw [birdDet_eq_birdDetSpec, birdDetSpec_eq_det] |
| 304 | + |
| 305 | +end BirdDet |
0 commit comments