CLRS-Lean
Chapters 1–35 · CLRS Fourth Edition
Chapter 28. Matrix Operations

Chapter 28 — Matrix Operations

CLRS, fourth edition · Lean 4 formalization

The proofs below use the models and assumptions described in the scope and implementation notes.

Imports

import Mathlib
import CLRSLean.Chapter_03.Section_03_1_Asymptotic_Notation

28.1. Solving Systems of Linear Equations

This section formalizes the LUP decomposition (CLRS Theorem 28.1): every nonsingular matrix A over a field admits a factorization P·A = L·U where P is a permutation matrix, L is unit lower-triangular, and U is upper-triangular. The decomposition underlies Gaussian elimination, the determinant, and matrix inversion.

Main results:

Theorem exists_lup_decomposition: every nonsingular matrix admits an LUP decomposition σ.permMatrix · A = L · U. The proof is by induction on n: the pivot row is swapped into row 0, one Gaussian-elimination step (elimination) zeroes the subdiagonal, the induction hypothesis is applied to the nonsingular Schur complement (nonsingularity via det_block_schur), and the block factors are assembled with the finOneSumFin reindexing. Supporting lemmas cover pivot selection (exists_col_zero_ne_zero), the elimination step (elimination_unitLowerTriangular / elimination_mul_zero_zero / elimination_mul_col_zero), the unit-triangular determinant (det_unitLowerTriangular), and the permutation-matrix bookkeeping (fromBlocks_one_zero_zero_permMatrix / conjPermMatrix).
Theorem lup_solve_correct (CLRS §28.1, Algorithm LUP-SOLVE): if σ.permMatrix · A = L · U is an LUP decomposition and the substitution equations L·y = σ.permMatrix·b and U·x = y hold, then A·x = b — forward then backward substitution through the factors solves the system.
Theorem forwardSubst_spec (CLRS Lemma 28.1): forwardSubst L b, forward substitution through a unit lower-triangular L, satisfies L·(forwardSubst L b) = b.
Theorem backSubst_spec (CLRS Lemma 28.2): backSubst U y, backward substitution through an upper-triangular U with nonzero diagonal, satisfies U·(backSubst U y) = y.
Theorem lupSolve_correct: the constructive solver lupSolve σ L U b (forward-then-back substitution through the factors) solves A·x = b given an LUP decomposition of A.
Theorem exists_solution_of_nonsingular: a nonsingular matrix over a field solves every linear system (∃ x, A·x = b).
Theorems unique_solution_of_nonsingular, unique_solution_unitLowerTriangular (Lemma 28.1), and unique_solution_upperTriangular (Lemma 28.2): nonsingular, unit-lower-triangular, and upper-triangular systems with nonzero diagonal have at most one solution.
Theorem det_eq_sign_mul_det_of_lup (Corollary to Theorem 28.1): from an LUP decomposition σ.permMatrix · A = L · U with L unit lower-triangular, det U = sign σ · det A — determinants agree up to sign; det_ne_zero_of_lup gives U nonsingular when A is.
The legacy substitutionCost_isBigO, lupDecompositionCost_isBigO, matrixInversionCost_isBigO, and choleskyCost_isBigO prove upper bounds for abstract numerical budgets, not two-sided execution bounds. The ExecutableLUP companion separately proves actual decomposition and solve counters with exact-field correctness. Inversion and Cholesky budgets here are not measured execution counters.

Notation conventions:

A : an n × n matrix over a field F.
σ.permMatrix F : the permutation matrix of σ.
IsUpperTriangular / IsLowerTriangular: zero above/below the main diagonal.
IsUnitLowerTriangular: lower-triangular with diagonal ones.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F]

Upper-triangular: every entry strictly below the main diagonal is zero.

def IsUpperTriangular {n : ℕ} (M : Matrix (Fin n) (Fin n) F) : Prop :=
  ∀ ⦃i j : Fin n⦄, j < i → M i j = 0

Lower-triangular: every entry strictly above the main diagonal is zero.

def IsLowerTriangular {n : ℕ} (M : Matrix (Fin n) (Fin n) F) : Prop :=
  ∀ ⦃i j : Fin n⦄, i < j → M i j = 0

Unit lower-triangular: lower-triangular with every diagonal entry equal to one (the shape of the L factor in an LUP decomposition).

def IsUnitLowerTriangular {n : ℕ} (M : Matrix (Fin n) (Fin n) F) : Prop :=
  IsLowerTriangular M ∧ ∀ i : Fin n, M i i = 1

The determinant of a matrix with an all-zero column is zero.


lemma det_eq_zero_of_col_zero {n : ℕ} {A : Matrix (Fin (n + 1)) (Fin (n + 1)) F}
    (h : ∀ i : Fin (n + 1), A i 0 = 0) : A.det = 0 := by
  classical
  -- det A = Σ σ, sign σ · ∏ i, A (σ i) i; column 0 zero kills every term.
  rw [Matrix.det_apply]
  apply Finset.sum_eq_zero
  intro σ hσ
  have hfactor : A (σ 0) 0 = 0 := h (σ 0)
  have hprod : (∏ i : Fin (n + 1), A (σ i) i) = 0 := by
    exact Finset.prod_eq_zero (Finset.mem_univ _) hfactor
  simp [hprod]

A nonsingular matrix has a nonzero entry in its first column (a pivot).


lemma exists_col_zero_ne_zero {n : ℕ} {A : Matrix (Fin (n + 1)) (Fin (n + 1)) F}
    (hA : A.det ≠ 0) : ∃ p : Fin (n + 1), A p 0 ≠ 0 := by
  classical
  by_contra h
  have hzero : ∀ i : Fin (n + 1), A i 0 = 0 := by
    intro i
    exact not_ne_iff.mp (fun hne => h ⟨i, hne⟩)
  exact hA (det_eq_zero_of_col_zero hzero)

Multiplying by a permutation matrix on the left permutes the rows: (σ.permMatrix * A) i j = A (σ i) j.


lemma permMatrix_mul_apply {n : ℕ} (σ : Equiv.Perm (Fin n))
    (A : Matrix (Fin n) (Fin n) F) (i j : Fin n) :
    (σ.permMatrix F * A) i j = A (σ i) j := by
  rw [Matrix.mul_apply]
  rw [Finset.sum_eq_single (σ i)]
  · simp
  · intro b hb hbne
    simp [Equiv.Perm.permMatrix, PEquiv.toMatrix_apply, hbne.symm]
  · simp

Swapping the pivot row p into row 0 makes the leading entry nonzero: (swap 0 p).permMatrix * A has a nonzero (0,0) entry.


lemma perm_mul_zero_zero_ne_zero {n : ℕ} {A : Matrix (Fin (n + 1)) (Fin (n + 1)) F}
    {p : Fin (n + 1)} (hp : A p 0 ≠ 0) :
    ((Equiv.swap 0 p).permMatrix F * A) 0 0 ≠ 0 := by
  have h : ((Equiv.swap 0 p).permMatrix F * A) 0 0 = A p 0 := by
    rw [permMatrix_mul_apply]
    simp
  rw [h]
  exact hp

The Gaussian-elimination step matrix: unit lower-triangular, with the subdiagonal entries of column 0 chosen to zero them out in E * B.

def elimination {n : ℕ} (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (_h : B 0 0 ≠ 0) :
    Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
  fun i j => if i = j then 1 else if j = 0 then -B i 0 / B 0 0 else 0

The elimination matrix is unit lower-triangular.


lemma elimination_unitLowerTriangular {n : ℕ} (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (h : B 0 0 ≠ 0) : IsUnitLowerTriangular (elimination B h) := by
  constructor
  · intro i j hij
    unfold elimination
    by_cases hEq : i = j
    · exfalso; exact (ne_of_lt hij) hEq
    · have hj0 : j ≠ 0 := by
        intro hj
        exact (not_lt_of_ge (by simp [hj]) (hij : i < j))
      simp [hEq, hj0]
  · intro i
    unfold elimination
    simp

The pivot entry (0,0) is unchanged by the elimination step.


lemma elimination_mul_zero_zero {n : ℕ} (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (h : B 0 0 ≠ 0) : (elimination B h * B) 0 0 = B 0 0 := by
  rw [Matrix.mul_apply]
  rw [Finset.sum_eq_single (0 : Fin (n + 1))]
  · unfold elimination
    simp
  · intro b hb hb0
    unfold elimination
    by_cases hEq : (0 : Fin (n + 1)) = b
    · exfalso; exact hb0 hEq.symm
    · have hb0' : b ≠ (0 : Fin (n + 1)) := by
        intro hb'
        exact hEq hb'.symm
      simp [hEq, hb0']
  · simp

The elimination step zeroes out column 0 below the pivot.


lemma elimination_mul_col_zero {n : ℕ} (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (h : B 0 0 ≠ 0) (i : Fin n) : (elimination B h * B) (Fin.succ i) 0 = 0 := by
  rw [Matrix.mul_apply]
  have hterm : ∀ k : Fin (n + 1), k ≠ Fin.succ i → k ≠ 0 →
      elimination B h (Fin.succ i) k * B k 0 = 0 := by
    intro k hk1 hk0
    unfold elimination
    have hne : Fin.succ i ≠ k := by exact Ne.symm hk1
    simp [hne, hk0]
  have hself : elimination B h (Fin.succ i) (Fin.succ i) * B (Fin.succ i) 0 =
      B (Fin.succ i) 0 := by
    unfold elimination
    simp
  have hzero : elimination B h (Fin.succ i) 0 * B 0 0 =
      (-B (Fin.succ i) 0 / B 0 0) * B 0 0 := by
    unfold elimination
    simp
  let S : Finset (Fin (n + 1)) := ({Fin.succ i, 0} : Finset (Fin (n + 1)))
  have hsubset : S ⊆ Finset.univ := by intro k hk; simp
  have hsubsum : (∑ k ∈ S, elimination B h (Fin.succ i) k * B k 0) =
      (∑ k : Fin (n + 1), elimination B h (Fin.succ i) k * B k 0) := by
    refine Finset.sum_subset hsubset ?_
    intro k hk hnot
    have hneq : ¬ k = Fin.succ i ∧ ¬ k = 0 := not_or.mp (by simpa [S] using hnot)
    exact hterm k hneq.1 hneq.2
  have hS : (∑ k ∈ S, elimination B h (Fin.succ i) k * B k 0) =
      B (Fin.succ i) 0 + (-B (Fin.succ i) 0 / B 0 0) * B 0 0 := by
    have hpair : ({Fin.succ i, 0} : Finset (Fin (n + 1))) =
        insert (Fin.succ i) ({0} : Finset (Fin (n + 1))) := by
      ext k
      simp [Finset.mem_insert]
    dsimp [S]
    rw [hpair]
    rw [Finset.sum_insert (by simp)]
    rw [Finset.sum_singleton]
    rw [hself, hzero]
  have hsum : (∑ k : Fin (n + 1), elimination B h (Fin.succ i) k * B k 0) =
      B (Fin.succ i) 0 + (-B (Fin.succ i) 0 / B 0 0) * B 0 0 := by
    rw [← hsubsum]
    exact hS
  rw [hsum]
  field_simp [h]
  ring

A permutation of Fin n other than the identity sends some index strictly below itself (so the identity is the only ≤-monotone bijection).


lemma exists_lt_of_perm_ne_id {n : ℕ} {σ : Equiv.Perm (Fin n)} (hσ : σ ≠ 1) :
    ∃ i : Fin n, σ i < i := by
  classical
  by_contra h
  have hle : ∀ i : Fin n, i ≤ σ i := by
    intro i
    exact le_of_not_gt (fun hgt => h ⟨i, hgt⟩)
  have hsum : (∑ i : Fin n, (σ i).val) = ∑ i : Fin n, i.val := by
    simpa using (Equiv.sum_comp σ (fun i : Fin n => i.val))
  by_cases hall : ∀ i : Fin n, σ i = i
  · exact hσ (Equiv.ext (fun i => hall i))
  · rcases not_forall.mp hall with ⟨i, hne⟩
    have hgt : i < σ i := lt_of_le_of_ne (hle i) (Ne.symm hne)
    have hlt : (∑ i : Fin n, i.val) < (∑ i : Fin n, (σ i).val) := by
      exact Finset.sum_lt_sum (fun j hj => hle j) ⟨i, by simp, hgt⟩
    exact (not_lt_of_ge (le_of_eq hsum)) hlt

A unit lower-triangular matrix has determinant one: in the expansion, every non-identity permutation σ picks a factor M (σ i) i with σ i < i, which is zero.


lemma det_unitLowerTriangular {n : ℕ} {M : Matrix (Fin n) (Fin n) F}
    (hM : IsUnitLowerTriangular M) : M.det = 1 := by
  classical
  rw [Matrix.det_apply]
  have hId : Equiv.Perm.sign (1 : Equiv.Perm (Fin n)) • ∏ i : Fin n, M i i = 1 := by
    simp [hM.2]
  have hone : ∀ σ : Equiv.Perm (Fin n), σ ≠ 1 →
      Equiv.Perm.sign σ • ∏ i : Fin n, M (σ i) i = 0 := by
    intro σ hσ
    rcases exists_lt_of_perm_ne_id hσ with ⟨i, hi⟩
    have hzero : M (σ i) i = 0 := hM.1 hi
    have hprod : (∏ i : Fin n, M (σ i) i) = 0 := by
      exact Finset.prod_eq_zero (Finset.mem_univ i) hzero
    simp [hprod]
  rw [Finset.sum_eq_single (1 : Equiv.Perm (Fin n))]
  · simpa using hId
  · intro σ hσ hne
    exact hone σ hne
  · intro hσ1
    exfalso
    exact hσ1 (by simp)

Reindex Fin (n + 1) as Fin 1 ⊕ Fin n: row 0 ↦ inl (), row succ i ↦ inr i. This is the bookkeeping that lets an (n + 1) × (n + 1) matrix be viewed as a 2 × 2 block matrix with a 1 × 1 top-left block.


noncomputable def finOneSumFin (n : ℕ) : Fin (n + 1) ≃ Fin 1 ⊕ Fin n where
  toFun i := if h : i = 0 then Sum.inl ⟨0, by omega⟩ else Sum.inr (i.pred h)
  invFun x := Sum.elim (fun _ => (0 : Fin (n + 1))) Fin.succ x
  left_inv := by
    intro i
    by_cases h : i = 0
    · simp [h]
    · simp [h]
  right_inv := by
    intro x
    cases x with
    | inl j =>
        have hj : (⟨0, by omega⟩ : Fin 1) = j := by
          ext
          omega
        simp [hj]
    | inr i =>
        simp

@[simp]
lemma finOneSumFin_symm_inl (n : ℕ) (j : Fin 1) :
    (finOneSumFin n).symm (Sum.inl j) = (0 : Fin (n + 1)) := by
  rfl

@[simp]
lemma finOneSumFin_symm_inr (n : ℕ) (i : Fin n) :
    (finOneSumFin n).symm (Sum.inr i) = Fin.succ i := by
  rfl

The determinant of a 1×1 matrix is its single entry.

lemma det_const_fin_one (x : F) : Matrix.det (fun _ _ : Fin 1 => x) = x := by
  simp

The determinant of a block-lower-triangular matrix: if the first column of C below the pivot is all zero, then det C = C 0 0 · det M where M is the Schur complement (the succ/succ submatrix). This is the n + 1 case of CLRS Lemma 28.1's determinant step.


lemma det_block_schur {n : ℕ} (C : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (h : ∀ i : Fin n, C (Fin.succ i) 0 = 0) :
    C.det = C 0 0 * Matrix.det (fun i j : Fin n => C (Fin.succ i) (Fin.succ j)) := by
  classical
  let e : Fin (n + 1) ≃ Fin 1 ⊕ Fin n := finOneSumFin n
  let C' : Matrix (Fin 1 ⊕ Fin n) (Fin 1 ⊕ Fin n) F := C.reindex e e
  have hdet : C'.det = C.det := by
    simp [C']
  have h21 : C'.toBlocks₂₁ = 0 := by
    ext i j
    simp [C', Matrix.reindex_apply, e, Matrix.toBlocks₂₁]
    exact h i
  have h11 : C'.toBlocks₁₁ = (fun _ _ : Fin 1 => C 0 0) := by
    ext i j
    simp [C', Matrix.reindex_apply, e, Matrix.toBlocks₁₁]
  have h22 : C'.toBlocks₂₂ = (fun i j : Fin n => C (Fin.succ i) (Fin.succ j)) := by
    ext i j
    simp [C', Matrix.reindex_apply, e, Matrix.toBlocks₂₂]
  calc
    C.det = C'.det := by rw [hdet]
    _ = (Matrix.fromBlocks C'.toBlocks₁₁ C'.toBlocks₁₂ C'.toBlocks₂₁ C'.toBlocks₂₂).det := by
      rw [Matrix.fromBlocks_toBlocks]
    _ = C'.toBlocks₁₁.det * C'.toBlocks₂₂.det := by
      rw [h21]
      rw [Matrix.det_fromBlocks_zero₂₁]
    _ = C 0 0 * Matrix.det (fun i j : Fin n => C (Fin.succ i) (Fin.succ j)) := by
      rw [h11, h22]
      rw [det_const_fin_one]

The block-diagonal matrix diag(1, P₁) with a 1 × 1 top-left block is the permutation matrix of the sum-of-permutations sumCongr 1 σ₁.

lemma fromBlocks_one_zero_zero_permMatrix {n : ℕ} (σ₁ : Equiv.Perm (Fin n)) :
    Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 0 (σ₁.permMatrix F) =
      (Equiv.Perm.sumCongr (1 : Equiv.Perm (Fin 1)) σ₁).permMatrix F := by
  ext i j
  cases i <;> cases j <;>
    simp [Matrix.fromBlocks, Matrix.one_apply, Equiv.Perm.permMatrix, PEquiv.toMatrix_apply,
      Equiv.toPEquiv_apply, Pi.single_apply, eq_comm]

Conjugating a permutation by an equivalence reindexes its permutation matrix: (e.trans σ.trans e.symm).permMatrix = σ.permMatrix.reindex e.symm e.symm.

lemma conjPermMatrix {m n : Type*} [DecidableEq m] [DecidableEq n] (e : m ≃ n)
    (σ : Equiv.Perm n) :
    Equiv.Perm.permMatrix F ((e.trans σ).trans e.symm) = (σ.permMatrix F).reindex e.symm e.symm := by
  ext i j
  simp [Equiv.Perm.permMatrix, PEquiv.toMatrix_apply, Equiv.toPEquiv_apply, Matrix.reindex_apply,
    Matrix.submatrix_apply, Equiv.trans_apply, Equiv.symm_symm, Equiv.apply_eq_iff_eq_symm_apply]

If i < j in Fin (n + 1) then j is not the zero row.

set_option linter.unnecessarySimpa false in


lemma ne_zero_of_lt_succ {n : ℕ} {i j : Fin (n + 1)} (hij : i < j) : j ≠ 0 := by
  intro hj0
  have hi0 : i < (0 : Fin (n + 1)) := by simpa [hj0] using hij
  have hi0v : i.val < 0 := by simpa using hi0
  exact (Nat.not_lt_zero i.val) hi0v

If j < i in Fin (n + 1) then i is not the zero row.

set_option linter.unnecessarySimpa false in


lemma ne_zero_of_lt {n : ℕ} {i j : Fin (n + 1)} (hji : j < i) : i ≠ 0 := by
  intro hi0
  have hj0 : j < (0 : Fin (n + 1)) := by simpa [hi0] using hji
  have hj0v : j.val < 0 := by simpa using hj0
  exact (Nat.not_lt_zero j.val) hj0v

Theorem 28.1 (LUP decomposition). Every nonsingular n × n matrix A over a field admits an LUP decomposition: a permutation σ, a unit lower-triangular L, and an upper-triangular U such that σ.permMatrix · A = L · U.

The proof is by induction on n. In the inductive step we pivot the first column so the leading entry is nonzero, perform one Gaussian-elimination step to zero the subdiagonal, apply the induction hypothesis to the Schur complement, and assemble the resulting block matrices.


theorem exists_lup_decomposition {n : ℕ} (A : Matrix (Fin n) (Fin n) F) (hA : A.det ≠ 0) :
    ∃ (σ : Equiv.Perm (Fin n)) (L U : Matrix (Fin n) (Fin n) F),
      IsUnitLowerTriangular L ∧ IsUpperTriangular U ∧ σ.permMatrix F * A = L * U := by
  classical
  induction n with
  | zero =>
      refine ⟨1, 1, 1, ?_, ?_, ?_⟩
      · simp [IsUnitLowerTriangular, IsLowerTriangular]
      · simp [IsUpperTriangular]
      · simp
        ext i j
        exact Fin.elim0 i
  | succ n ih =>
      obtain ⟨p, hp⟩ := exists_col_zero_ne_zero hA
      let σ : Equiv.Perm (Fin (n + 1)) := Equiv.swap 0 p
      let B : Matrix (Fin (n + 1)) (Fin (n + 1)) F := σ.permMatrix F * A
      have hB : B 0 0 ≠ 0 := by
        dsimp [B, σ]
        exact perm_mul_zero_zero_ne_zero hp
      let E : Matrix (Fin (n + 1)) (Fin (n + 1)) F := elimination B hB
      let D : Matrix (Fin (n + 1)) (Fin (n + 1)) F := E * B
      let M : Matrix (Fin n) (Fin n) F := fun i j => D (Fin.succ i) (Fin.succ j)
      have hDcol : ∀ i : Fin n, D (Fin.succ i) 0 = 0 := by
        intro i
        dsimp [D, E]
        exact elimination_mul_col_zero B hB i
      have hMdet : M.det ≠ 0 := by
        have hBdet : B.det ≠ 0 := by
          dsimp [B, σ]
          rw [Matrix.det_mul]
          have hrinv : (Equiv.swap 0 p)⁻¹.permMatrix F * (Equiv.swap 0 p).permMatrix F = 1 := by
            rw [← Matrix.permMatrix_mul]
            simp
          have hdetσ : ((Equiv.swap 0 p).permMatrix F).det ≠ 0 := by
            exact Matrix.det_ne_zero_of_left_inverse hrinv
          exact mul_ne_zero hdetσ hA
        have hDdet : D.det ≠ 0 := by
          dsimp [D, E]
          rw [Matrix.det_mul]
          rw [det_unitLowerTriangular (elimination_unitLowerTriangular B hB)]
          simpa using hBdet
        have hblock : D.det = B 0 0 * M.det := by
          dsimp [M, D, E]
          rw [det_block_schur (elimination B hB * B) hDcol]
          rw [elimination_mul_zero_zero B hB]
        have hprod : B 0 0 * M.det ≠ 0 := by
          exact hblock ▸ hDdet
        exact (mul_ne_zero_iff.mp hprod).2
      obtain ⟨σ₁, L₁, U₁, hL₁, hU₁, hMfac⟩ := ih M hMdet
      let re : Fin (n + 1) ≃ Fin 1 ⊕ Fin n := finOneSumFin n
      let α : Matrix (Fin 1) (Fin 1) F := fun _ _ => D 0 0
      let v : Matrix (Fin 1) (Fin n) F := fun _ j => D 0 (Fin.succ j)
      let mult : Matrix (Fin n) (Fin 1) F := fun i _ => -E (Fin.succ i) 0
      let P₁ : Matrix (Fin n) (Fin n) F := σ₁.permMatrix F
      let diagP : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
        (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 0 P₁).reindex re.symm re.symm
      let E' : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
        (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 mult
          (1 : Matrix (Fin n) (Fin n) F)).reindex re.symm re.symm
      let L : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
        (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 (P₁ * mult) L₁).reindex re.symm re.symm
      let U : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
        (Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) U₁).reindex re.symm re.symm
      let σP : Equiv.Perm (Fin (n + 1)) :=
        (re.trans (Equiv.Perm.sumCongr (1 : Equiv.Perm (Fin 1)) σ₁)).trans re.symm
      let σ₀ : Equiv.Perm (Fin (n + 1)) := σ * σP
      have hDblocks : D.reindex re re =
          Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) M := by
        rw [Matrix.ext_iff_blocks]
        constructor
        · ext i j
          simp [re, α, Matrix.reindex_apply, Matrix.toBlocks₁₁]
        · constructor
          · ext i j
            simp [re, v, Matrix.reindex_apply, Matrix.toBlocks₁₂]
          · constructor
            · ext i j
              simp [re, Matrix.reindex_apply, Matrix.toBlocks₂₁]
              exact hDcol i
            · ext i j
              simp [re, M, Matrix.reindex_apply, Matrix.toBlocks₂₂]
      have hEblocks : E.reindex re re =
          Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F)
            (0 : Matrix (Fin 1) (Fin n) F) (-mult) (1 : Matrix (Fin n) (Fin n) F) := by
        rw [Matrix.ext_iff_blocks]
        constructor
        · ext i j
          simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₁₁, Matrix.one_apply]
          exact Subsingleton.elim i j
        · constructor
          · ext i j
            simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₁₂, eq_comm]
          · constructor
            · ext i j
              simp [re, E, elimination, mult, Matrix.reindex_apply, Matrix.toBlocks₂₁]
            · ext i j
              simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₂₂, Matrix.one_apply]
      have hEinv : E' * E = 1 := by
        apply (reindexRingEquiv F re).injective
        rw [map_mul, map_one]
        rw [show (reindexRingEquiv F re) E' = Matrix.fromBlocks 1 0 mult
            (1 : Matrix (Fin n) (Fin n) F) by
          simp [E', Matrix.reindex_apply, Matrix.submatrix_submatrix]]
        rw [show (reindexRingEquiv F re) E = Matrix.fromBlocks 1 0 (-mult)
            (1 : Matrix (Fin n) (Fin n) F) by
          simpa using hEblocks]
        rw [Matrix.fromBlocks_multiply]
        simp
      have hσP : Equiv.Perm.permMatrix F σP = diagP := by
        dsimp [σP, diagP]
        rw [conjPermMatrix re (Equiv.Perm.sumCongr (1 : Equiv.Perm (Fin 1)) σ₁)]
        rw [← fromBlocks_one_zero_zero_permMatrix σ₁]
        dsimp [P₁]
      have hDiagE'D : diagP * E' * D = L * U := by
        apply (reindexRingEquiv F re).injective
        rw [map_mul, map_mul]
        rw [show (reindexRingEquiv F re) D = Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) M by
          simpa using hDblocks]
        rw [map_mul]
        rw [show (reindexRingEquiv F re) diagP = Matrix.fromBlocks 1 0 0 P₁ by
          simp [diagP, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
        rw [show (reindexRingEquiv F re) E' = Matrix.fromBlocks 1 0 mult
            (1 : Matrix (Fin n) (Fin n) F) by
          simp [E', Matrix.reindex_apply, Matrix.submatrix_submatrix]]
        rw [show (reindexRingEquiv F re) L = Matrix.fromBlocks 1 0 (P₁ * mult) L₁ by
          simp [L, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
        rw [show (reindexRingEquiv F re) U = Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) U₁ by
          simp [U, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
        simp [Matrix.fromBlocks_multiply, P₁, hMfac]
      have hBD : B = E' * D := by
        calc
          B = 1 * B := by rw [one_mul]
          _ = (E' * E) * B := by rw [hEinv]
          _ = E' * (E * B) := by rw [Matrix.mul_assoc]
          _ = E' * D := by rfl
      have hDiagB : diagP * B = L * U := by
        calc
          diagP * B = diagP * (E' * D) := by rw [hBD]
          _ = (diagP * E') * D := by rw [Matrix.mul_assoc]
          _ = L * U := hDiagE'D
      have hEq : σ₀.permMatrix F * A = L * U := by
        calc
          σ₀.permMatrix F * A = (σ * σP).permMatrix F * A := by rfl
          _ = (σP.permMatrix F * σ.permMatrix F) * A := by rw [Matrix.permMatrix_mul]
          _ = σP.permMatrix F * (σ.permMatrix F * A) := by rw [Matrix.mul_assoc]
          _ = σP.permMatrix F * B := by rfl
          _ = diagP * B := by rw [hσP]
          _ = L * U := hDiagB
      have hL : IsUnitLowerTriangular L := by
        constructor
        · intro i j hij
          by_cases hi : i = 0
          · have hj : j ≠ 0 := ne_zero_of_lt_succ hij
            dsimp only [L]
            rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
            have hri : re i = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
              dsimp [re]
              simp [hi, finOneSumFin]
            have hrj : re j = Sum.inr (j.pred hj) := by
              dsimp [re]
              simp [hj, finOneSumFin]
            rw [hri, hrj]
            simp
          · have hj : j ≠ 0 := ne_zero_of_lt_succ hij
            dsimp only [L]
            rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
            have hri : re i = Sum.inr (i.pred hi) := by
              dsimp [re]
              simp [hi, finOneSumFin]
            have hrj : re j = Sum.inr (j.pred hj) := by
              dsimp [re]
              simp [hj, finOneSumFin]
            rw [hri, hrj]
            simp
            apply hL₁.1
            have hsi : i = Fin.succ (i.pred hi) := (Fin.succ_pred i hi).symm
            have hsj : j = Fin.succ (j.pred hj) := (Fin.succ_pred j hj).symm
            have hlt : Fin.succ (i.pred hi) < Fin.succ (j.pred hj) := by
              rw [hsi, hsj] at hij
              exact hij
            exact Fin.succ_lt_succ_iff.mp hlt
        · intro i
          by_cases hi : i = 0
          · dsimp only [L]
            rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
            have hri : re i = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
              dsimp [re]
              simp [hi, finOneSumFin]
            rw [hri]
            simp
          · dsimp only [L]
            rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
            have hri : re i = Sum.inr (i.pred hi) := by
              dsimp [re]
              simp [hi, finOneSumFin]
            rw [hri]
            simp
            exact hL₁.2 (i.pred hi)
      have hU : IsUpperTriangular U := by
        intro i j hji
        by_cases hj : j = 0
        · have hi : i ≠ 0 := ne_zero_of_lt hji
          dsimp only [U]
          rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
          have hrj : re j = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
            dsimp [re]
            simp [hj, finOneSumFin]
          have hri : re i = Sum.inr (i.pred hi) := by
            dsimp [re]
            simp [hi, finOneSumFin]
          rw [hrj, hri]
          simp
        · have hi : i ≠ 0 := ne_zero_of_lt hji
          dsimp only [U]
          rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
          have hrj : re j = Sum.inr (j.pred hj) := by
            dsimp [re]
            simp [hj, finOneSumFin]
          have hri : re i = Sum.inr (i.pred hi) := by
            dsimp [re]
            simp [hi, finOneSumFin]
          rw [hrj, hri]
          simp
          apply hU₁
          have hsj : j = Fin.succ (j.pred hj) := (Fin.succ_pred j hj).symm
          have hsi : i = Fin.succ (i.pred hi) := (Fin.succ_pred i hi).symm
          have hlt : Fin.succ (j.pred hj) < Fin.succ (i.pred hi) := by
            rw [hsj, hsi] at hji
            exact hji
          exact Fin.succ_lt_succ_iff.mp hlt
      exact ⟨σ₀, L, U, hL, hU, hEq⟩

LUP-SOLVE correctness (CLRS §28.1, Algorithm LUP-SOLVE). If σ.permMatrix · A = L · U is an LUP decomposition and y, x are obtained by forward and backward substitution (L·y = σ.permMatrix·b, U·x = y), then x solves the linear system A·x = b.

The proof composes the two substitution equations through the factorization: σ.permMatrix·(A·x) = (σ.permMatrix·A)·x = (L·U)·x = L·(U·x) = L·y = σ.permMatrix·b, then cancels the permutation matrix σ.permMatrix (whose mulVec is the bijection v ↦ v∘σ).


theorem lup_solve_correct {n : ℕ} {A L U : Matrix (Fin n) (Fin n) F} {σ : Equiv.Perm (Fin n)}
    (hLUP : σ.permMatrix F * A = L * U) (b x y : Fin n → F)
    (hLy : L *ᵥ y = σ.permMatrix F *ᵥ b) (hUx : U *ᵥ x = y) :
    A *ᵥ x = b := by
  have hPAx : σ.permMatrix F *ᵥ (A *ᵥ x) = σ.permMatrix F *ᵥ b := by
    rw [Matrix.mulVec_mulVec]
    rw [hLUP]
    rw [← Matrix.mulVec_mulVec]
    rw [hUx]
    exact hLy
  have hcomp : (A *ᵥ x) ∘ σ = b ∘ σ := by
    rw [Matrix.permMatrix_mulVec, Matrix.permMatrix_mulVec] at hPAx
    exact hPAx
  have h := congrArg (fun f : Fin n → F => f ∘ σ.symm) hcomp
  simpa [Function.comp_assoc] using h

A nonsingular matrix over a field solves every linear system: if A.det ≠ 0 then for every right-hand side b there is x with A·x = b.


theorem exists_solution_of_nonsingular {n : ℕ} (A : Matrix (Fin n) (Fin n) F) (hA : A.det ≠ 0)
    (b : Fin n → F) : ∃ x : Fin n → F, A *ᵥ x = b := by
  have hunit : IsUnit A := by
    rw [Matrix.isUnit_iff_isUnit_det]
    exact isUnit_iff_ne_zero.mpr hA
  exact Matrix.mulVec_surjective_iff_isUnit.2 hunit b

Forward substitution (CLRS Lemma 28.1). forwardSubst L b is the vector y computed by forward substitution through the unit lower-triangular matrix L; the recursion splits off the first component y₀ = b₀ and substitutes the tail through the trailing block L₂.

The constructor forwardSubst is total (any L); its correctness (forwardSubst_spec) requires L to be unit lower-triangular.

noncomputable def forwardSubst : ∀ {n : ℕ}, Matrix (Fin n) (Fin n) F → (Fin n → F) → (Fin n → F)
  | 0, _, _ => fun i => Fin.elim0 i
  | n + 1, L, b =>
      let L2 : Matrix (Fin n) (Fin n) F := fun i j => L (Fin.succ i) (Fin.succ j)
      Fin.cons (b 0) (forwardSubst L2 (fun i => b (Fin.succ i) - L (Fin.succ i) 0 * b 0))

The trailing block of a unit lower-triangular matrix is unit lower-triangular.

lemma unitLowerTriangular_tail {n : ℕ} (L : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (hL : IsUnitLowerTriangular L) :
    IsUnitLowerTriangular (fun i j : Fin n => L (Fin.succ i) (Fin.succ j)) := by
  constructor
  · intro i j hij
    exact hL.1 (Fin.succ_lt_succ_iff.mpr hij)
  · intro i
    exact hL.2 (Fin.succ i)

Forward substitution solves L·y = b (CLRS Lemma 28.1). For a unit lower-triangular L, forwardSubst L b satisfies L·(forwardSubst L b) = b.


theorem forwardSubst_spec : ∀ {n : ℕ} (L : Matrix (Fin n) (Fin n) F) (b : Fin n → F),
    IsUnitLowerTriangular L → L *ᵥ forwardSubst L b = b
  | 0, _, _, _ => by
      ext i
      exact Fin.elim0 i
  | n + 1, L, b, hL =>
      by
      let L2 : Matrix (Fin n) (Fin n) F := fun i j => L (Fin.succ i) (Fin.succ j)
      let b2 : Fin n → F := fun i => b (Fin.succ i) - L (Fin.succ i) 0 * b 0
      have hL2 : IsUnitLowerTriangular L2 := by
        simpa [L2] using (unitLowerTriangular_tail L hL)
      have hspec : L2 *ᵥ forwardSubst L2 b2 = b2 := forwardSubst_spec L2 b2 hL2
      change L *ᵥ Fin.cons (b 0) (forwardSubst L2 b2) = b
      ext i
      rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨i', rfl⟩
      · unfold Matrix.mulVec dotProduct
        rw [Fin.sum_univ_succ]
        have h00 : L 0 0 = 1 := hL.2 0
        have h0k : ∀ k : Fin n, L 0 (Fin.succ k) = 0 := fun k => hL.1 (Fin.succ_pos k)
        simp [Fin.cons_zero, Fin.cons_succ, h00, h0k]
      · unfold Matrix.mulVec dotProduct
        rw [Fin.sum_univ_succ]
        simp [Fin.cons_zero, Fin.cons_succ]
        change L (Fin.succ i') 0 * b 0 + (L2 *ᵥ forwardSubst L2 b2) i' = b (Fin.succ i')
        rw [hspec]
        simp [b2]

Backward substitution (CLRS Lemma 28.2). backSubst U y is the vector x computed by backward substitution through the upper-triangular matrix U; the recursion splits off the last component xₙ = yₙ/Uₙₙ and substitutes the tail through the leading block U₁.

The constructor backSubst is total (any U); its correctness (backSubst_spec) requires U to be upper-triangular with nonzero diagonal.

noncomputable def backSubst : ∀ {n : ℕ}, Matrix (Fin n) (Fin n) F → (Fin n → F) → (Fin n → F)
  | 0, _, _ => fun i => Fin.elim0 i
  | n + 1, U, y =>
      let last : Fin (n + 1) := Fin.last n
      let U1 : Matrix (Fin n) (Fin n) F := fun i j => U (Fin.castSucc i) (Fin.castSucc j)
      let xl : F := y last / U last last
      Fin.snoc (backSubst U1 (fun i => y (Fin.castSucc i) - U (Fin.castSucc i) last * xl)) xl

The leading block of an upper-triangular matrix is upper-triangular.

lemma upperTriangular_tail {n : ℕ} (U : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (hU : IsUpperTriangular U) :
    IsUpperTriangular (fun i j : Fin n => U (Fin.castSucc i) (Fin.castSucc j)) := by
  intro i j hji
  exact hU (Fin.castSucc_lt_castSucc_iff.mpr hji)

Backward substitution solves U·x = y (CLRS Lemma 28.2). For an upper-triangular U with nonzero diagonal, backSubst U y satisfies U·(backSubst U y) = y.


theorem backSubst_spec : ∀ {n : ℕ} (U : Matrix (Fin n) (Fin n) F) (y : Fin n → F),
    IsUpperTriangular U → (∀ i : Fin n, U i i ≠ 0) → U *ᵥ backSubst U y = y
  | 0, _, _, _, _ => by ext i; exact Fin.elim0 i
  | n + 1, U, y, hU, hdiag =>
      by
      let last : Fin (n + 1) := Fin.last n
      let U1 : Matrix (Fin n) (Fin n) F := fun i j => U (Fin.castSucc i) (Fin.castSucc j)
      let xl : F := y last / U last last
      let y1 : Fin n → F := fun i => y (Fin.castSucc i) - U (Fin.castSucc i) last * xl
      have hU1 : IsUpperTriangular U1 := by
        simpa [U1] using (upperTriangular_tail U hU)
      have hdiag1 : ∀ i : Fin n, U1 i i ≠ 0 := by
        intro i
        exact hdiag (Fin.castSucc i)
      have hspec : U1 *ᵥ backSubst U1 y1 = y1 := backSubst_spec U1 y1 hU1 hdiag1
      change U *ᵥ Fin.snoc (backSubst U1 y1) xl = y
      ext i
      rcases Fin.eq_castSucc_or_eq_last i with ⟨i', rfl⟩ | rfl
      · unfold Matrix.mulVec dotProduct
        rw [Fin.sum_univ_castSucc]
        simp [Fin.snoc_castSucc, Fin.snoc_last]
        change (U1 *ᵥ backSubst U1 y1) i' + U (Fin.castSucc i') last * xl = y (Fin.castSucc i')
        rw [hspec]
        simp [y1]
      · unfold Matrix.mulVec dotProduct
        rw [Fin.sum_univ_castSucc]
        have hUlast : ∀ j' : Fin n, U last (Fin.castSucc j') = 0 := by
          intro j'
          exact hU (Fin.castSucc_lt_last j')
        simp [Fin.snoc_castSucc, Fin.snoc_last]
        have hz : (∑ x : Fin n, U (Fin.last n) (Fin.castSucc x) * backSubst U1 y1 x) = 0 := by
          refine Finset.sum_eq_zero ?_
          intro x hx
          have h0 : U (Fin.last n) (Fin.castSucc x) = 0 := hU (Fin.castSucc_lt_last x)
          simp [h0]
        rw [hz]
        simp [xl]
        have hln : U last last ≠ 0 := hdiag last
        field_simp [hln]
        ring

LUP-SOLVE. lupSolve σ L U b solves A·x = b through an LUP decomposition σ.permMatrix · A = L · U by forward-substituting L·y = σ.permMatrix·b and then back-substituting U·x = y. It is the constructive counterpart of the compositional theorem lup_solve_correct.

noncomputable def lupSolve {n : ℕ} (σ : Equiv.Perm (Fin n)) (L U : Matrix (Fin n) (Fin n) F)
    (b : Fin n → F) : Fin n → F :=
  backSubst U (forwardSubst L (σ.permMatrix F *ᵥ b))

LUP-SOLVE correctness. If σ.permMatrix · A = L · U is an LUP decomposition with L unit lower-triangular and U upper-triangular with nonzero diagonal, then lupSolve σ L U b solves A·x = b.


theorem lupSolve_correct {n : ℕ} {A L U : Matrix (Fin n) (Fin n) F} {σ : Equiv.Perm (Fin n)}
    (hLUP : σ.permMatrix F * A = L * U) (hL : IsUnitLowerTriangular L)
    (hU : IsUpperTriangular U) (hUdiag : ∀ i : Fin n, U i i ≠ 0) (b : Fin n → F) :
    A *ᵥ lupSolve σ L U b = b := by
  unfold lupSolve
  let Pb : Fin n → F := σ.permMatrix F *ᵥ b
  let y : Fin n → F := forwardSubst L Pb
  have hLy : L *ᵥ y = σ.permMatrix F *ᵥ b := by
    simpa [y, Pb] using (forwardSubst_spec L Pb hL)
  have hUx : U *ᵥ (backSubst U y) = y := by
    exact backSubst_spec U y hU hUdiag
  exact lup_solve_correct hLUP b (backSubst U y) y hLy hUx

A nonsingular matrix over a field has a unique solution for every linear system: if A·x₁ = b and A·x₂ = b then x₁ = x₂.


theorem unique_solution_of_nonsingular {n : ℕ} (A : Matrix (Fin n) (Fin n) F) (hA : A.det ≠ 0)
    (b x1 x2 : Fin n → F) (h1 : A *ᵥ x1 = b) (h2 : A *ᵥ x2 = b) : x1 = x2 := by
  have hunit : IsUnit A := by
    rw [Matrix.isUnit_iff_isUnit_det]
    exact isUnit_iff_ne_zero.mpr hA
  apply (Matrix.mulVec_injective_iff_isUnit.2 hunit)
  calc
    A *ᵥ x1 = b := h1
    _ = A *ᵥ x2 := h2.symm

A unit lower-triangular system L·x = b has at most one solution (CLRS Lemma 28.1: forward substitution computes the unique one).


theorem unique_solution_unitLowerTriangular {n : ℕ} (L : Matrix (Fin n) (Fin n) F)
    (hL : IsUnitLowerTriangular L) (b x1 x2 : Fin n → F)
    (h1 : L *ᵥ x1 = b) (h2 : L *ᵥ x2 = b) : x1 = x2 := by
  have hunit : IsUnit L := by
    rw [Matrix.isUnit_iff_isUnit_det]
    rw [det_unitLowerTriangular hL]
    exact isUnit_one
  apply (Matrix.mulVec_injective_iff_isUnit.2 hunit)
  calc
    L *ᵥ x1 = b := h1
    _ = L *ᵥ x2 := h2.symm

An upper-triangular system U·x = y with nonzero diagonal has at most one solution (CLRS Lemma 28.2: backward substitution computes the unique one).


theorem unique_solution_upperTriangular {n : ℕ} (U : Matrix (Fin n) (Fin n) F)
    (hU : IsUpperTriangular U) (hdiag : ∀ i : Fin n, U i i ≠ 0)
    (y x1 x2 : Fin n → F) (h1 : U *ᵥ x1 = y) (h2 : U *ᵥ x2 = y) : x1 = x2 := by
  have hdet : U.det ≠ 0 := by
    rw [Matrix.det_of_upperTriangular]
    · exact Finset.prod_ne_zero_iff.2 (fun i hi => hdiag i)
    · intro i j hij
      exact hU hij
  have hunit : IsUnit U := by
    rw [Matrix.isUnit_iff_isUnit_det]
    exact isUnit_iff_ne_zero.mpr hdet
  apply (Matrix.mulVec_injective_iff_isUnit.2 hunit)
  calc
    U *ᵥ x1 = y := h1
    _ = U *ᵥ x2 := h2.symm

Determinant from the LUP decomposition (CLRS Corollary to Theorem 28.1). If σ.permMatrix · A = L · U with L unit lower-triangular, then det U = sign σ · det A: the determinants of A and U agree up to sign, so a determinant is computable from an LUP decomposition in cubic time.


theorem det_eq_sign_mul_det_of_lup {n : ℕ} {A L U : Matrix (Fin n) (Fin n) F}
    {σ : Equiv.Perm (Fin n)} (hLUP : σ.permMatrix F * A = L * U) (hL : IsUnitLowerTriangular L) :
    (Equiv.Perm.sign σ : F) * A.det = U.det := by
  calc
    (Equiv.Perm.sign σ : F) * A.det = (σ.permMatrix F).det * A.det := by rw [Matrix.det_permutation]
    _ = (σ.permMatrix F * A).det := by rw [Matrix.det_mul]
    _ = (L * U).det := by rw [hLUP]
    _ = L.det * U.det := by rw [Matrix.det_mul]
    _ = 1 * U.det := by rw [det_unitLowerTriangular hL]
    _ = U.det := by ring

In an LUP decomposition of a nonsingular matrix, the upper factor U is nonsingular.


theorem det_ne_zero_of_lup {n : ℕ} {A L U : Matrix (Fin n) (Fin n) F}
    {σ : Equiv.Perm (Fin n)} (hLUP : σ.permMatrix F * A = L * U) (hL : IsUnitLowerTriangular L)
    (hA : A.det ≠ 0) : U.det ≠ 0 := by
  intro hU0
  have h := det_eq_sign_mul_det_of_lup hLUP hL
  have hz : (Equiv.Perm.sign σ : F) * A.det = 0 := by
    rw [h, hU0]
  have hc_ne : (Equiv.Perm.sign σ : F) ≠ 0 := by
    rcases Int.units_eq_one_or (Equiv.Perm.sign σ) with h | h
    · rw [h]
      norm_num
    · rw [h]
      norm_num
  have hA0 : A.det = 0 := by
    calc
      A.det = (Equiv.Perm.sign σ : F)⁻¹ * ((Equiv.Perm.sign σ : F) * A.det) := by
        field_simp [hc_ne]
      _ = 0 := by rw [hz]; ring
  exact hA hA0

An upper-triangular matrix with nonzero determinant has every diagonal entry nonzero.


lemma upperTriangular_diag_ne_zero_of_det_ne_zero {n : ℕ} (U : Matrix (Fin n) (Fin n) F)
    (hU : IsUpperTriangular U) (hdet : U.det ≠ 0) : ∀ i : Fin n, U i i ≠ 0 := by
  intro i
  have hdiag : U.det = ∏ i, U i i := Matrix.det_of_upperTriangular (by
    intro a b hab
    exact hU hab)
  intro hzero
  apply hdet
  rw [hdiag]
  exact Finset.prod_eq_zero (Finset.mem_univ i) hzero

section Cost

Abstract forward+backward substitution cost for LUP-SOLVE: the numerical envelope n²/2. This is not the exact two-loop operation count.

noncomputable def substitutionCost (n : ℕ) : ℝ := (n : ℝ) * (n : ℝ) / 2

Abstract cubic LUP-decomposition envelope n³/3, separate from the execution counter in the ExecutableLUP companion.

noncomputable def lupDecompositionCost (n : ℕ) : ℝ := (n : ℝ) ^ 3 / 3

Abstract cubic matrix-inversion budget; no executed inversion counter is attached here.

noncomputable def matrixInversionCost (n : ℕ) : ℝ := (n : ℝ) ^ 3

Abstract cubic Cholesky budget; no executed Cholesky counter is attached here.

noncomputable def choleskyCost (n : ℕ) : ℝ := (n : ℝ) ^ 3

The substitution envelope is at most n²: n²/2 ≤ n².

theorem substitutionCost_quadratic_bound (n : ℕ) : substitutionCost n ≤ (n : ℝ) ^ 2 := by
  unfold substitutionCost
  nlinarith [sq_nonneg (n : ℝ)]

The numerical substitution envelope has an O(n²) upper bound.


theorem substitutionCost_isBigO :
    CLRS.Chapter03.isBigO substitutionCost (fun n => (n : ℝ) ^ 2) := by
  rw [CLRS.Chapter03.isBigO_iff]
  refine ⟨1, by norm_num, 0, fun n hn => ?_⟩
  have hnonneg : 0 ≤ substitutionCost n := by
    unfold substitutionCost
    positivity
  have hnonneg2 : 0 ≤ (n : ℝ) ^ 2 := sq_nonneg (n : ℝ)
  rw [abs_of_nonneg hnonneg, abs_of_nonneg hnonneg2]
  simpa using (substitutionCost_quadratic_bound n)

The LUP-decomposition cost is at most n³: n³/3 ≤ n³.


theorem lupDecompositionCost_cubic_bound (n : ℕ) : lupDecompositionCost n ≤ (n : ℝ) ^ 3 := by
  unfold lupDecompositionCost
  have h3 : (0 : ℝ) ≤ (n : ℝ) ^ 3 := by positivity
  nlinarith

The numerical LUP-decomposition envelope has an O(n³) upper bound.


theorem lupDecompositionCost_isBigO :
    CLRS.Chapter03.isBigO lupDecompositionCost (fun n => (n : ℝ) ^ 3) := by
  rw [CLRS.Chapter03.isBigO_iff]
  refine ⟨1, by norm_num, 0, fun n hn => ?_⟩
  have hnonneg : 0 ≤ lupDecompositionCost n := by
    unfold lupDecompositionCost
    positivity
  have hnonneg2 : 0 ≤ (n : ℝ) ^ 3 := by positivity
  rw [abs_of_nonneg hnonneg, abs_of_nonneg hnonneg2]
  simpa using (lupDecompositionCost_cubic_bound n)

The abstract matrix-inversion budget has an O(n³) upper bound.

theorem matrixInversionCost_isBigO :
    CLRS.Chapter03.isBigO matrixInversionCost (fun n => (n : ℝ) ^ 3) := by
  exact CLRS.Chapter03.isBigO_refl matrixInversionCost

The abstract Cholesky budget has an O(n³) upper bound.

theorem choleskyCost_isBigO :
    CLRS.Chapter03.isBigO choleskyCost (fun n => (n : ℝ) ^ 3) := by
  exact CLRS.Chapter03.isBigO_refl choleskyCost

              end Costend Chapter28end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Basic

CLRS Section 28.1 - Executable LUP foundations

Result records and the concrete pivot/elimination executions used by the dimension-recursive LUP decomposition.

              namespace CLRSnamespace Chapter28open Matrix

The three factors returned by a successful LUP decomposition.

structure LUPFactors (n : Nat) (F : Type) where
  perm : Equiv.Perm (Fin n)
  lower : Matrix (Fin n) (Fin n) F
  upper : Matrix (Fin n) (Fin n) F

Result and exact-algebra work of a total LUP execution: pivot comparisons plus the field operations performed by elimination and factor assembly.

structure LUPExecution (n : Nat) (F : Type) where
  result : Option (LUPFactors n F)
  work : Nat

Result and comparison count of scanning a matrix's first column.

structure PivotExecution {n : Nat} {F : Type} [Zero F]
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) where
  pivot : Option {p : Fin (n + 1) // A p 0 ≠ 0}
  comparisons : Nat

Result and work of one field-valued arithmetic expression.

structure FieldExecution (F : Type) where
  value : F
  work : Nat

Result and field-operation work of a square-matrix construction.

structure MatrixExecution (n : Nat) (F : Type) where
  value : Matrix (Fin n) (Fin n) F
  work : Nat

section Pivotvariable {F : Type} [Zero F] [DecidableEq F]

Scan candidate rows from left to right, charging one comparison for each tested first-column entry.

def findPivotListWithCost {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) :
    List (Fin (n + 1)) → PivotExecution A
  | [] => ⟨none, 0⟩
  | p :: ps =>
      if h : A p 0 = 0 then
        let rest := findPivotListWithCost A ps
        ⟨rest.pivot, rest.comparisons + 1⟩
      else
        ⟨some ⟨p, h⟩, 1⟩

Scan every row for the first nonzero pivot in column zero.

def findPivotWithCost {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) : PivotExecution A :=
  findPivotListWithCost A (List.finRange (n + 1))

              end Pivotsection Eliminationvariable {F : Type} [Field F]

One entry of the direct Gaussian-elimination update. A copied pivot-row entry is free in the field-operation model; every other entry performs one division, one multiplication, and one subtraction.

def eliminateEntryWithCost {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (_h : B 0 0 ≠ 0)
    (i j : Fin (n + 1)) : FieldExecution F :=
  if i = 0 then
    ⟨B i j, 0⟩
  else
    ⟨B i j - (B i 0 / B 0 0) * B 0 j, 3⟩

Pointwise Gaussian elimination together with the sum of the entry-level field-operation counters.

def eliminateWithCost {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0) :
    MatrixExecution (n + 1) F :=
  ⟨fun i j => (eliminateEntryWithCost B h i j).value,
    ∑ i : Fin (n + 1), ∑ j : Fin (n + 1), (eliminateEntryWithCost B h i j).work⟩

              end Eliminationend Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Correctness

CLRS Section 28.1 - Executable LUP correctness

Successful recursive results satisfy the exact triangular and factorization contract. Nonsingular inputs succeed, so failure characterizes singularity.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F] [DecidableEq F]

omit [DecidableEq F] in
private theorem assembleLUPFactors_correct {n : Nat}
    (A B D : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (p : Fin (n + 1))
    (hB0 : B 0 0 ≠ 0)
    (hB : B = (Equiv.swap 0 p).permMatrix F * A)
    (hD : D = elimination B hB0 * B)
    (child : LUPFactors n F)
    (hchildL : IsUnitLowerTriangular child.lower)
    (hchildU : IsUpperTriangular child.upper)
    (hchildFac :
      let M : Matrix (Fin n) (Fin n) F := fun i j => D (Fin.succ i) (Fin.succ j)
      child.perm.permMatrix F * M = child.lower * child.upper) :
    let factors := assembleLUPFactors B D (Equiv.swap 0 p) child
    IsUnitLowerTriangular factors.lower ∧
      IsUpperTriangular factors.upper ∧
      factors.perm.permMatrix F * A = factors.lower * factors.upper := by
  let re : Fin (n + 1) ≃ Fin 1 ⊕ Fin n := execFinOneSumFin n
  let E : Matrix (Fin (n + 1)) (Fin (n + 1)) F := elimination B hB0
  let M : Matrix (Fin n) (Fin n) F := fun i j => D (Fin.succ i) (Fin.succ j)
  let α : Matrix (Fin 1) (Fin 1) F := fun _ _ => D 0 0
  let v : Matrix (Fin 1) (Fin n) F := fun _ j => D 0 (Fin.succ j)
  let mult : Matrix (Fin n) (Fin 1) F := fun i _ => B (Fin.succ i) 0 / B 0 0
  let P₁ : Matrix (Fin n) (Fin n) F := child.perm.permMatrix F
  let permutedMult : Matrix (Fin n) (Fin 1) F :=
    fun i j => mult (child.perm i) j
  have hPermutedMult : permutedMult = P₁ * mult := by
    ext i j
    dsimp only [permutedMult, P₁]
    simp [Matrix.mul_apply, Equiv.Perm.permMatrix, PEquiv.toMatrix_apply]
  have hMfac : P₁ * M = child.lower * child.upper := by
    simpa [P₁, M] using hchildFac
  let diagP : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 0 P₁).reindex re.symm re.symm
  let E' : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 mult
      (1 : Matrix (Fin n) (Fin n) F)).reindex re.symm re.symm
  let L : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 permutedMult child.lower).reindex
      re.symm re.symm
  let U : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) child.upper).reindex
      re.symm re.symm
  let σP : Equiv.Perm (Fin (n + 1)) := liftTrailingPerm child.perm
  let σ₀ : Equiv.Perm (Fin (n + 1)) := Equiv.swap 0 p * σP
  have hDcol : ∀ i : Fin n, D (Fin.succ i) 0 = 0 := by
    intro i
    rw [hD]
    exact elimination_mul_col_zero B hB0 i
  have hDblocks : D.reindex re re =
      Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) M := by
    rw [Matrix.ext_iff_blocks]
    constructor
    · ext i j
      simp [re, α, Matrix.reindex_apply, Matrix.toBlocks₁₁]
    · constructor
      · ext i j
        simp [re, v, Matrix.reindex_apply, Matrix.toBlocks₁₂]
      · constructor
        · ext i j
          simp [re, Matrix.reindex_apply, Matrix.toBlocks₂₁]
          exact hDcol i
        · ext i j
          simp [re, M, Matrix.reindex_apply, Matrix.toBlocks₂₂]
  have hEblocks : E.reindex re re =
      Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F)
        (0 : Matrix (Fin 1) (Fin n) F) (-mult) (1 : Matrix (Fin n) (Fin n) F) := by
    rw [Matrix.ext_iff_blocks]
    constructor
    · ext i j
      simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₁₁, Matrix.one_apply]
      exact Subsingleton.elim i j
    · constructor
      · ext i j
        simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₁₂, eq_comm]
      · constructor
        · ext i j
          simp [re, E, elimination, mult, Matrix.reindex_apply, Matrix.toBlocks₂₁]
          ring
        · ext i j
          simp [re, E, elimination, Matrix.reindex_apply, Matrix.toBlocks₂₂, Matrix.one_apply]
  have hEinv : E' * E = 1 := by
    apply (reindexRingEquiv F re).injective
    rw [map_mul, map_one]
    rw [show (reindexRingEquiv F re) E' = Matrix.fromBlocks 1 0 mult
        (1 : Matrix (Fin n) (Fin n) F) by
      simp [E', Matrix.reindex_apply, Matrix.submatrix_submatrix]]
    rw [show (reindexRingEquiv F re) E = Matrix.fromBlocks 1 0 (-mult)
        (1 : Matrix (Fin n) (Fin n) F) by
      simpa using hEblocks]
    rw [Matrix.fromBlocks_multiply]
    simp
  have hσP : Equiv.Perm.permMatrix F σP = diagP := by
    dsimp [σP, liftTrailingPerm, diagP]
    rw [conjPermMatrix re (Equiv.Perm.sumCongr (1 : Equiv.Perm (Fin 1)) child.perm)]
    rw [← fromBlocks_one_zero_zero_permMatrix child.perm]
    rfl
  have hDiagE'D : diagP * E' * D = L * U := by
    apply (reindexRingEquiv F re).injective
    rw [map_mul, map_mul]
    rw [show (reindexRingEquiv F re) D = Matrix.fromBlocks α v
        (0 : Matrix (Fin n) (Fin 1) F) M by simpa using hDblocks]
    rw [map_mul]
    rw [show (reindexRingEquiv F re) diagP = Matrix.fromBlocks 1 0 0 P₁ by
      simp [diagP, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
    rw [show (reindexRingEquiv F re) E' = Matrix.fromBlocks 1 0 mult
        (1 : Matrix (Fin n) (Fin n) F) by
      simp [E', Matrix.reindex_apply, Matrix.submatrix_submatrix]]
    rw [show (reindexRingEquiv F re) L = Matrix.fromBlocks 1 0 (P₁ * mult) child.lower by
      rw [← hPermutedMult]
      simp [L, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
    rw [show (reindexRingEquiv F re) U = Matrix.fromBlocks α v
        (0 : Matrix (Fin n) (Fin 1) F) child.upper by
      simp [U, Matrix.reindex_apply, Matrix.submatrix_submatrix]]
    simp [Matrix.fromBlocks_multiply, hMfac]
  have hBD : B = E' * D := by
    calc
      B = 1 * B := by rw [one_mul]
      _ = (E' * E) * B := by rw [hEinv]
      _ = E' * (E * B) := by rw [Matrix.mul_assoc]
      _ = E' * D := by rw [hD]
  have hDiagB : diagP * B = L * U := by
    calc
      diagP * B = diagP * (E' * D) := by rw [hBD]
      _ = (diagP * E') * D := by rw [Matrix.mul_assoc]
      _ = L * U := hDiagE'D
  have hEq : σ₀.permMatrix F * A = L * U := by
    calc
      σ₀.permMatrix F * A = ((Equiv.swap 0 p) * σP).permMatrix F * A := by rfl
      _ = (σP.permMatrix F * (Equiv.swap 0 p).permMatrix F) * A := by
        rw [Matrix.permMatrix_mul]
      _ = σP.permMatrix F * ((Equiv.swap 0 p).permMatrix F * A) := by
        rw [Matrix.mul_assoc]
      _ = σP.permMatrix F * B := by rw [hB]
      _ = diagP * B := by rw [hσP]
      _ = L * U := hDiagB
  have hL : IsUnitLowerTriangular L := by
    constructor
    · intro i j hij
      by_cases hi : i = 0
      · have hj : j ≠ 0 := ne_zero_of_lt_succ hij
        dsimp only [L]
        rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
        have hri : re i = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
          dsimp [re, execFinOneSumFin]
          simp [hi]
        have hrj : re j = Sum.inr (j.pred hj) := by
          dsimp [re, execFinOneSumFin]
          simp [hj]
        rw [hri, hrj]
        simp
      · have hj : j ≠ 0 := ne_zero_of_lt_succ hij
        dsimp only [L]
        rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
        have hri : re i = Sum.inr (i.pred hi) := by
          dsimp [re, execFinOneSumFin]
          simp [hi]
        have hrj : re j = Sum.inr (j.pred hj) := by
          dsimp [re, execFinOneSumFin]
          simp [hj]
        rw [hri, hrj]
        simp
        apply hchildL.1
        have hsi : i = Fin.succ (i.pred hi) := (Fin.succ_pred i hi).symm
        have hsj : j = Fin.succ (j.pred hj) := (Fin.succ_pred j hj).symm
        rw [hsi, hsj] at hij
        exact Fin.succ_lt_succ_iff.mp hij
    · intro i
      by_cases hi : i = 0
      · dsimp only [L]
        rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
        have hri : re i = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
          dsimp [re, execFinOneSumFin]
          simp [hi]
        rw [hri]
        simp
      · dsimp only [L]
        rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
        have hri : re i = Sum.inr (i.pred hi) := by
          dsimp [re, execFinOneSumFin]
          simp [hi]
        rw [hri]
        simp
        exact hchildL.2 (i.pred hi)
  have hU : IsUpperTriangular U := by
    intro i j hji
    by_cases hj : j = 0
    · have hi : i ≠ 0 := ne_zero_of_lt hji
      dsimp only [U]
      rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
      have hrj : re j = Sum.inl (⟨0, by omega⟩ : Fin 1) := by
        dsimp [re, execFinOneSumFin]
        simp [hj]
      have hri : re i = Sum.inr (i.pred hi) := by
        dsimp [re, execFinOneSumFin]
        simp [hi]
      rw [hrj, hri]
      simp
    · have hi : i ≠ 0 := ne_zero_of_lt hji
      dsimp only [U]
      rw [Matrix.reindex_apply, Matrix.submatrix_apply, Equiv.symm_symm]
      have hrj : re j = Sum.inr (j.pred hj) := by
        dsimp [re, execFinOneSumFin]
        simp [hj]
      have hri : re i = Sum.inr (i.pred hi) := by
        dsimp [re, execFinOneSumFin]
        simp [hi]
      rw [hrj, hri]
      simp
      apply hchildU
      have hsj : j = Fin.succ (j.pred hj) := (Fin.succ_pred j hj).symm
      have hsi : i = Fin.succ (i.pred hi) := (Fin.succ_pred i hi).symm
      rw [hsj, hsi] at hji
      exact Fin.succ_lt_succ_iff.mp hji
  change IsUnitLowerTriangular
      (assembleLUPFactors B D (Equiv.swap 0 p) child).lower ∧
    IsUpperTriangular (assembleLUPFactors B D (Equiv.swap 0 p) child).upper ∧
    (assembleLUPFactors B D (Equiv.swap 0 p) child).perm.permMatrix F * A =
      (assembleLUPFactors B D (Equiv.swap 0 p) child).lower *
        (assembleLUPFactors B D (Equiv.swap 0 p) child).upper
  simpa [assembleLUPFactors, re, α, v, mult, P₁, permutedMult, L, U, σP, σ₀,
    liftTrailingPerm] using And.intro hL (And.intro hU hEq)

Direct row lookup agrees with multiplying by the selected swap matrix.

omit [DecidableEq F] in


theorem pivotedMatrix_eq_permMatrix_mul {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (p : Fin (n + 1)) :
    pivotedMatrix A p = (Equiv.swap 0 p).permMatrix F * A := by
  ext i j
  rw [permMatrix_mul_apply]
  rfl

Every factor triple returned by the recursive execution is a valid LUP decomposition.


theorem lupDecomposeWithCost_sound : ∀ {n : Nat}
    (A : Matrix (Fin n) (Fin n) F) (factors : LUPFactors n F),
    (lupDecomposeWithCost n A).result = some factors →
      IsUnitLowerTriangular factors.lower ∧
        IsUpperTriangular factors.upper ∧
        factors.perm.permMatrix F * A = factors.lower * factors.upper
  | 0, A, factors, hresult => by
      have hf : (⟨1, 1, 1⟩ : LUPFactors 0 F) = factors := by
        exact Option.some.inj (by simpa [lupDecomposeWithCost] using hresult)
      subst factors
      constructor
      · simp [IsUnitLowerTriangular, IsLowerTriangular]
      · constructor
        · simp [IsUpperTriangular]
        · ext i j
          exact Fin.elim0 i
  | n + 1, A, factors, hresult => by
      cases hpivot : (findPivotWithCost A).pivot with
      | none =>
          simp [lupDecomposeWithCost, hpivot] at hresult
      | some p =>
          let pivotPerm : Equiv.Perm (Fin (n + 1)) := Equiv.swap 0 p.1
          let B := pivotedMatrix A p.1
          have hp0 : A p.1 0 ≠ 0 := p.2
          have hB0 : B 0 0 ≠ 0 := by
            simpa [B, pivotedMatrix] using hp0
          let eliminated := eliminateWithCost B hB0
          let M : Matrix (Fin n) (Fin n) F :=
            fun i j => eliminated.value (Fin.succ i) (Fin.succ j)
          let child := lupDecomposeWithCost n M
          cases hchild : child.result with
          | none =>
              have hchild' :
                  (lupDecomposeWithCost n M).result = none := by
                simpa [child] using hchild
              simp [lupDecomposeWithCost, hpivot, B, eliminated, M, hchild'] at hresult
          | some childFactors =>
              have hchild' :
                  (lupDecomposeWithCost n M).result = some childFactors := by
                simpa [child] using hchild
              have hfactor :
                  assembleLUPFactors B eliminated.value pivotPerm childFactors = factors := by
                apply Option.some.inj
                simpa [lupDecomposeWithCost, hpivot, pivotPerm, B, hp0, hB0,
                  eliminated, M, child, hchild'] using hresult
              have hchildResult :
                  (lupDecomposeWithCost n M).result = some childFactors := by
                exact hchild'
              have hchildCorrect := lupDecomposeWithCost_sound M childFactors hchildResult
              subst factors
              apply assembleLUPFactors_correct A B eliminated.value p.1 hB0
              · exact pivotedMatrix_eq_permMatrix_mul A p.1
              · exact eliminateWithCost_value B hB0
              · exact hchildCorrect.1
              · exact hchildCorrect.2.1
              · simpa [M] using hchildCorrect.2.2

Every nonsingular matrix makes the explicit pivot/elimination recursion return a factor triple.


theorem lupDecomposeWithCost_nonsingular : ∀ {n : Nat}
    (A : Matrix (Fin n) (Fin n) F), A.det ≠ 0 →
      ∃ factors, (lupDecomposeWithCost n A).result = some factors
  | 0, _A, _hA => by
      exact ⟨⟨1, 1, 1⟩, rfl⟩
  | n + 1, A, hA => by
      cases hpivot : (findPivotWithCost A).pivot with
      | none =>
          have hcol : ∀ p : Fin (n + 1), A p 0 = 0 :=
            (findPivotWithCost_eq_none A).1 hpivot
          exact False.elim (hA (det_eq_zero_of_col_zero hcol))
      | some p =>
          let B := pivotedMatrix A p.1
          have hB0 : B 0 0 ≠ 0 := by
            simpa [B, pivotedMatrix] using p.2
          let eliminated := eliminateWithCost B hB0
          let M : Matrix (Fin n) (Fin n) F :=
            fun i j => eliminated.value (Fin.succ i) (Fin.succ j)
          have hBdet : B.det ≠ 0 := by
            dsimp [B]
            rw [pivotedMatrix_eq_permMatrix_mul, Matrix.det_mul]
            have hrinv :
                (Equiv.swap 0 p.1)⁻¹.permMatrix F * (Equiv.swap 0 p.1).permMatrix F = 1 := by
              rw [← Matrix.permMatrix_mul]
              simp
            have hdetPerm : ((Equiv.swap 0 p.1).permMatrix F).det ≠ 0 :=
              Matrix.det_ne_zero_of_left_inverse hrinv
            exact mul_ne_zero hdetPerm hA
          have hDdet : eliminated.value.det ≠ 0 := by
            rw [eliminateWithCost_value, Matrix.det_mul]
            rw [det_unitLowerTriangular (elimination_unitLowerTriangular B hB0)]
            simpa using hBdet
          have hDcol : ∀ i : Fin n, eliminated.value (Fin.succ i) 0 = 0 := by
            exact eliminateWithCost_col_zero B hB0
          have hblock : eliminated.value.det = B 0 0 * M.det := by
            rw [det_block_schur eliminated.value hDcol]
            have h00 : eliminated.value 0 0 = B 0 0 := by
              rw [eliminateWithCost_value]
              exact elimination_mul_zero_zero B hB0
            rw [h00]
          have hMdet : M.det ≠ 0 := by
            have hprod : B 0 0 * M.det ≠ 0 := hblock ▸ hDdet
            exact (mul_ne_zero_iff.mp hprod).2
          obtain ⟨childFactors, hchild⟩ := lupDecomposeWithCost_nonsingular M hMdet
          refine ⟨assembleLUPFactors B eliminated.value (Equiv.swap 0 p.1) childFactors, ?_⟩
          simp [lupDecomposeWithCost, hpivot, B, eliminated, M, hchild]

A successful execution certifies nonsingularity; a singular matrix cannot pass every nonzero-pivot stage.


theorem lupDecomposeWithCost_success_det_ne_zero : ∀ {n : Nat}
    (A : Matrix (Fin n) (Fin n) F) (factors : LUPFactors n F),
    (lupDecomposeWithCost n A).result = some factors → A.det ≠ 0
  | 0, _A, _factors, _hresult => by simp
  | n + 1, A, factors, hresult => by
      cases hpivot : (findPivotWithCost A).pivot with
      | none =>
          simp [lupDecomposeWithCost, hpivot] at hresult
      | some p =>
          let B := pivotedMatrix A p.1
          have hB0 : B 0 0 ≠ 0 := by
            simpa [B, pivotedMatrix] using p.2
          let eliminated := eliminateWithCost B hB0
          let M : Matrix (Fin n) (Fin n) F :=
            fun i j => eliminated.value (Fin.succ i) (Fin.succ j)
          let child := lupDecomposeWithCost n M
          cases hchild : child.result with
          | none =>
              have hchild' : (lupDecomposeWithCost n M).result = none := by
                simpa [child] using hchild
              simp [lupDecomposeWithCost, hpivot, B, eliminated, M, hchild'] at hresult
          | some childFactors =>
              have hchild' :
                  (lupDecomposeWithCost n M).result = some childFactors := by
                simpa [child] using hchild
              have hMdet : M.det ≠ 0 :=
                lupDecomposeWithCost_success_det_ne_zero M childFactors hchild'
              have hDcol : ∀ i : Fin n, eliminated.value (Fin.succ i) 0 = 0 :=
                eliminateWithCost_col_zero B hB0
              have hblock : eliminated.value.det = B 0 0 * M.det := by
                rw [det_block_schur eliminated.value hDcol]
                have h00 : eliminated.value 0 0 = B 0 0 := by
                  rw [eliminateWithCost_value]
                  exact elimination_mul_zero_zero B hB0
                rw [h00]
              have hDdet : eliminated.value.det ≠ 0 := by
                rw [hblock]
                exact mul_ne_zero hB0 hMdet
              have hBdet : B.det ≠ 0 := by
                rw [eliminateWithCost_value, Matrix.det_mul] at hDdet
                rw [det_unitLowerTriangular (elimination_unitLowerTriangular B hB0)] at hDdet
                simpa using hDdet
              have hdetEq : B.det =
                  ((Equiv.swap 0 p.1).permMatrix F).det * A.det := by
                dsimp [B]
                rw [pivotedMatrix_eq_permMatrix_mul, Matrix.det_mul]
              have hprod : ((Equiv.swap 0 p.1).permMatrix F).det * A.det ≠ 0 :=
                hdetEq ▸ hBdet
              exact (mul_ne_zero_iff.mp hprod).2

Every successful execution has nonzero upper-triangular diagonal entries, as required by backward substitution.

theorem lupDecomposeWithCost_upper_diag_ne_zero {n : Nat}
    (A : Matrix (Fin n) (Fin n) F) (factors : LUPFactors n F)
    (hresult : (lupDecomposeWithCost n A).result = some factors) :
    ∀ i : Fin n, factors.upper i i ≠ 0 := by
  have hcorrect := lupDecomposeWithCost_sound A factors hresult
  have hA := lupDecomposeWithCost_success_det_ne_zero A factors hresult
  have hUdet := det_ne_zero_of_lup hcorrect.2.2 hcorrect.1 hA
  exact upperTriangular_diag_ne_zero_of_det_ne_zero factors.upper hcorrect.2.1 hUdet

The total execution returns failure exactly for singular matrices.


theorem lupDecomposeWithCost_eq_none_iff {n : Nat}
    (A : Matrix (Fin n) (Fin n) F) :
    (lupDecomposeWithCost n A).result = none ↔ A.det = 0 := by
  constructor
  · intro hnone
    by_contra hdet
    obtain ⟨factors, hsome⟩ := lupDecomposeWithCost_nonsingular A hdet
    rw [hnone] at hsome
    simp at hsome
  · intro hdet
    cases hresult : (lupDecomposeWithCost n A).result with
    | none => rfl
    | some factors =>
        exact False.elim
          ((lupDecomposeWithCost_success_det_ne_zero A factors hresult) hdet)

end Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Cost

CLRS Section 28.1 - Executable LUP work

The cubic bound is derived from the counter stored by the recursive execution, including pivot comparisons, pointwise elimination, successful factor assembly, and both successful and early-failure paths.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F] [DecidableEq F]private theorem lup_step_arithmetic (n : Nat) :
    (n + 1) + 3 * (n + 1) ^ 2 + 4 * n ^ 3 + n ≤ 4 * (n + 1) ^ 3 := by
  nlinarith [sq_nonneg (n : Int)]

Every execution path uses at most 4n³ pivot comparisons and field operations in the declared exact-algebra unit-cost model.


theorem lupDecomposeWithCost_work_le : ∀ {n : Nat}
    (A : Matrix (Fin n) (Fin n) F),
    (lupDecomposeWithCost n A).work ≤ 4 * n ^ 3
  | 0, _A => by simp [lupDecomposeWithCost]
  | n + 1, A => by
      have hpivotWork := findPivotWithCost_comparisons_le A
      cases hpivot : (findPivotWithCost A).pivot with
      | none =>
          have hwork : (lupDecomposeWithCost (n + 1) A).work =
              (findPivotWithCost A).comparisons := by
            simp [lupDecomposeWithCost, hpivot]
          rw [hwork]
          have hone : n + 1 ≤ 4 * (n + 1) ^ 3 := by
            nlinarith [sq_nonneg (n : Int)]
          exact le_trans hpivotWork hone
      | some p =>
          let B := pivotedMatrix A p.1
          have hB0 : B 0 0 ≠ 0 := by
            simpa [B, pivotedMatrix] using p.2
          let eliminated := eliminateWithCost B hB0
          let M : Matrix (Fin n) (Fin n) F :=
            fun i j => eliminated.value (Fin.succ i) (Fin.succ j)
          let child := lupDecomposeWithCost n M
          have helimWork : eliminated.work ≤ 3 * (n + 1) ^ 2 := by
            simpa [eliminated] using eliminateWithCost_work_le B hB0
          have hchildWork : child.work ≤ 4 * n ^ 3 := by
            simpa [child] using lupDecomposeWithCost_work_le M
          have hsum :
              (findPivotWithCost A).comparisons + eliminated.work + child.work + n ≤
                (n + 1) + 3 * (n + 1) ^ 2 + 4 * n ^ 3 + n := by
            omega
          have hwork : (lupDecomposeWithCost (n + 1) A).work ≤
              (findPivotWithCost A).comparisons + eliminated.work + child.work + n := by
            cases hchild : child.result <;>
              simp [lupDecomposeWithCost, hpivot, B, eliminated, M, child, hchild]
          exact hwork.trans (hsum.trans (lup_step_arithmetic n))

Public executable Theorem 28.1 bundle: one run returns certified LUP factors and satisfies the cubic work bound.

theorem lupDecomposeWithCost_correct {n : Nat}
    (A : Matrix (Fin n) (Fin n) F) (hA : A.det ≠ 0) :
    ∃ factors,
      (lupDecomposeWithCost n A).result = some factors ∧
      IsUnitLowerTriangular factors.lower ∧
      IsUpperTriangular factors.upper ∧
      (∀ i : Fin n, factors.upper i i ≠ 0) ∧
      factors.perm.permMatrix F * A = factors.lower * factors.upper ∧
      (lupDecomposeWithCost n A).work ≤ 4 * n ^ 3 := by
  obtain ⟨factors, hresult⟩ := lupDecomposeWithCost_nonsingular A hA
  have hcorrect := lupDecomposeWithCost_sound A factors hresult
  exact ⟨factors, hresult, hcorrect.1, hcorrect.2.1,
    lupDecomposeWithCost_upper_diag_ne_zero A factors hresult, hcorrect.2.2,
    lupDecomposeWithCost_work_le A⟩

end Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Elimination

CLRS Section 28.1 - Costed direct elimination

The pointwise execution is identified with multiplication by the existing Gaussian-elimination matrix and its field-operation sum is bounded.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F]

private theorem elimination_mul_row_zero {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0)
    (j : Fin (n + 1)) :
    (elimination B h * B) 0 j = B 0 j := by
  rw [Matrix.mul_apply]
  rw [Finset.sum_eq_single (0 : Fin (n + 1))]
  · simp [elimination]
  · intro k _hk hk0
    have h0k : (0 : Fin (n + 1)) ≠ k := by exact Ne.symm hk0
    simp [elimination, h0k, hk0]
  · simp

private theorem elimination_mul_row_succ {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0)
    (i j : Fin (n + 1)) (hi : i ≠ 0) :
    (elimination B h * B) i j = B i j - (B i 0 / B 0 0) * B 0 j := by
  rw [Matrix.mul_apply]
  let S : Finset (Fin (n + 1)) := {i, 0}
  have hsubset : S ⊆ Finset.univ := by intro k _hk; simp
  have hsubsum :
      (∑ k ∈ S, elimination B h i k * B k j) =
        ∑ k : Fin (n + 1), elimination B h i k * B k j := by
    refine Finset.sum_subset hsubset ?_
    intro k _hk hnot
    have hne : k ≠ i ∧ k ≠ 0 := by
      simpa [S, eq_comm] using hnot
    have hik : i ≠ k := Ne.symm hne.1
    simp [elimination, hik, hne.2]
  rw [← hsubsum]
  have hpair : ({i, 0} : Finset (Fin (n + 1))) = insert i {0} := by
    ext k
    simp [eq_comm]
  dsimp [S]
  rw [hpair, Finset.sum_insert (by simpa using hi), Finset.sum_singleton]
  simp [elimination, hi]
  ring

Erasing the entry counters gives multiplication by the existing elimination matrix.

theorem eliminateWithCost_value {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0) :
    (eliminateWithCost B h).value = elimination B h * B := by
  ext i j
  by_cases hi : i = 0
  · subst i
    simp only [eliminateWithCost, eliminateEntryWithCost, ↓reduceIte]
    exact (elimination_mul_row_zero B h j).symm
  · simp only [eliminateWithCost, eliminateEntryWithCost, hi, ↓reduceIte]
    exact (elimination_mul_row_succ B h i j hi).symm

The direct elimination execution zeroes column zero below the pivot.


theorem eliminateWithCost_col_zero {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0)
    (i : Fin n) : (eliminateWithCost B h).value (Fin.succ i) 0 = 0 := by
  rw [eliminateWithCost_value]
  exact elimination_mul_col_zero B h i

At most three field operations are charged for each matrix entry.

theorem eliminateWithCost_work_le {n : Nat}
    (B : Matrix (Fin (n + 1)) (Fin (n + 1)) F) (h : B 0 0 ≠ 0) :
    (eliminateWithCost B h).work ≤ 3 * (n + 1) ^ 2 := by
  unfold eliminateWithCost
  simp only
  calc
    (∑ i : Fin (n + 1), ∑ j : Fin (n + 1),
        (eliminateEntryWithCost B h i j).work)
        ≤ ∑ _i : Fin (n + 1), ∑ _j : Fin (n + 1), 3 := by
          apply Finset.sum_le_sum
          intro i _hi
          apply Finset.sum_le_sum
          intro j _hj
          by_cases hi0 : i = 0
          · simp [eliminateEntryWithCost, hi0]
          · simp [eliminateEntryWithCost, hi0]
    _ = 3 * (n + 1) ^ 2 := by simp [pow_two]; ring

end Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Execution

CLRS Section 28.1 - Recursive executable LUP decomposition

This file contains only data-producing definitions. The recursive execution uses the concrete pivot scan and direct elimination before assembling child factors by block reindexing. Permutations are direct index lookups, and the counter includes the multiplier divisions performed during factor assembly.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F] [DecidableEq F]

Computable reindexing of Fin (n + 1) as the pivot coordinate followed by the trailing Fin n coordinates.


def execFinOneSumFin (n : Nat) : Fin (n + 1) ≃ Fin 1 ⊕ Fin n where
  toFun i := if h : i = 0 then Sum.inl ⟨0, by omega⟩ else Sum.inr (i.pred h)
  invFun x := Sum.elim (fun _ => (0 : Fin (n + 1))) Fin.succ x
  left_inv := by
    intro i
    by_cases h : i = 0
    · simp [h]
    · simp [h]
  right_inv := by
    intro x
    cases x with
    | inl j =>
        have hj : (⟨0, by omega⟩ : Fin 1) = j := Subsingleton.elim _ _
        simp [hj]
    | inr i => simp

@[simp]
theorem execFinOneSumFin_symm_inl (n : Nat) (j : Fin 1) :
    (execFinOneSumFin n).symm (Sum.inl j) = (0 : Fin (n + 1)) := by
  rfl

@[simp]
theorem execFinOneSumFin_symm_inr (n : Nat) (i : Fin n) :
    (execFinOneSumFin n).symm (Sum.inr i) = Fin.succ i := by
  rfl

Move pivot row p to row zero by direct row lookup.

def pivotedMatrix {n : Nat} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (p : Fin (n + 1)) : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
  fun i j => A (Equiv.swap 0 p i) j

Lift a child permutation so it fixes the pivot coordinate.

def liftTrailingPerm {n : Nat} (σ : Equiv.Perm (Fin n)) :
    Equiv.Perm (Fin (n + 1)) :=
  let re := execFinOneSumFin n
  (re.trans (Equiv.Perm.sumCongr (1 : Equiv.Perm (Fin 1)) σ)).trans re.symm

Assemble the parent factors from a pivoted matrix, its eliminated value, and a successful decomposition of the trailing block.

def assembleLUPFactors {n : Nat}
    (B D : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (pivotPerm : Equiv.Perm (Fin (n + 1))) (child : LUPFactors n F) :
    LUPFactors (n + 1) F :=
  let re := execFinOneSumFin n
  let α : Matrix (Fin 1) (Fin 1) F := fun _ _ => D 0 0
  let v : Matrix (Fin 1) (Fin n) F := fun _ j => D 0 (Fin.succ j)
  let mult : Matrix (Fin n) (Fin 1) F := fun i _ => B (Fin.succ i) 0 / B 0 0
  -- Multiplication by a permutation matrix is implemented as direct row
  -- lookup, so factor assembly does not hide an `n × n` matrix product.
  let permutedMult : Matrix (Fin n) (Fin 1) F :=
    fun i j => mult (child.perm i) j
  let L : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks (1 : Matrix (Fin 1) (Fin 1) F) 0 permutedMult child.lower).reindex
      re.symm re.symm
  let U : Matrix (Fin (n + 1)) (Fin (n + 1)) F :=
    (Matrix.fromBlocks α v (0 : Matrix (Fin n) (Fin 1) F) child.upper).reindex
      re.symm re.symm
  ⟨pivotPerm * liftTrailingPerm child.perm, L, U⟩

Total LUP execution. Failure is data: no factor triple is returned, but the comparisons and field operations already performed remain recorded. A successful size-n+1 assembly charges the n multiplier divisions used in the lower-left block.

def lupDecomposeWithCost : ∀ (n : Nat),
    Matrix (Fin n) (Fin n) F → LUPExecution n F
  | 0, _A =>
      ⟨some ⟨1, 1, 1⟩, 0⟩
  | n + 1, A =>
      let pivotRun := findPivotWithCost A
      match pivotRun.pivot with
      | none => ⟨none, pivotRun.comparisons⟩
      | some p =>
          let pivotPerm : Equiv.Perm (Fin (n + 1)) := Equiv.swap 0 p.1
          let B := pivotedMatrix A p.1
          have hp0 : A p.1 0 ≠ 0 := p.2
          have hB : B 0 0 ≠ 0 := by
            simpa [B, pivotedMatrix] using hp0
          let eliminated := eliminateWithCost B hB
          let M : Matrix (Fin n) (Fin n) F :=
            fun i j => eliminated.value (Fin.succ i) (Fin.succ j)
          let child := lupDecomposeWithCost n M
          match child.result with
          | none =>
              ⟨none, pivotRun.comparisons + eliminated.work + child.work⟩
          | some factors =>
              ⟨some (assembleLUPFactors B eliminated.value pivotPerm factors),
                pivotRun.comparisons + eliminated.work + child.work + n⟩

end Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.Pivot

CLRS Section 28.1 - Executable pivot scan

Soundness, failure characterization, and comparison work for the concrete first-column scan.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Zero F] [DecidableEq F]theorem findPivotListWithCost_eq_none {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (ps : List (Fin (n + 1))) :
    (findPivotListWithCost A ps).pivot = none ↔
      ∀ p ∈ ps, A p 0 = 0 := by
  induction ps with
  | nil => simp [findPivotListWithCost]
  | cons q qs ih =>
      by_cases hq : A q 0 = 0
      · simp [findPivotListWithCost, hq, ih]
      · simp [findPivotListWithCost, hq]theorem findPivotListWithCost_comparisons_le {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    (ps : List (Fin (n + 1))) :
    (findPivotListWithCost A ps).comparisons ≤ ps.length := by
  induction ps with
  | nil => simp [findPivotListWithCost]
  | cons q qs ih =>
      by_cases hq : A q 0 = 0
      · simp [findPivotListWithCost, hq]
        omega
      · simp [findPivotListWithCost, hq]

A returned pivot is a genuinely nonzero first-column entry.

theorem findPivotWithCost_found {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F)
    {p : {p : Fin (n + 1) // A p 0 ≠ 0}}
    (_hp : (findPivotWithCost A).pivot = some p) : A p.1 0 ≠ 0 :=
  p.2

The pivot scan fails exactly when the complete first column is zero.


theorem findPivotWithCost_eq_none {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) :
    (findPivotWithCost A).pivot = none ↔ ∀ p : Fin (n + 1), A p 0 = 0 := by
  rw [findPivotWithCost, findPivotListWithCost_eq_none]
  constructor
  · intro h p
    exact h p (List.mem_finRange p)
  · intro h p _hp
    exact h p

At most one pivot comparison is charged per row.

theorem findPivotWithCost_comparisons_le {n : Nat}
    (A : Matrix (Fin (n + 1)) (Fin (n + 1)) F) :
    (findPivotWithCost A).comparisons ≤ n + 1 := by
  simpa [findPivotWithCost] using
    findPivotListWithCost_comparisons_le A (List.finRange (n + 1))

end Chapter28end CLRS

CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP.SolveCost

CLRS Section 28.1 - Costed LUP-SOLVE

Forward and backward substitution carry counters produced by the same dimension recursion as their returned vectors.

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F]

Vector result with exact-algebra field-operation work.

structure VectorExecution (n : Nat) (F : Type) where
  value : Fin n → F
  work : Nat

Costed forward substitution. Each trailing right-hand-side entry charges one multiplication and one subtraction.

def forwardSubstWithCost : ∀ {n : Nat},
    Matrix (Fin n) (Fin n) F → (Fin n → F) → VectorExecution n F
  | 0, _L, _b => ⟨fun i => Fin.elim0 i, 0⟩
  | n + 1, L, b =>
      let L2 : Matrix (Fin n) (Fin n) F := fun i j => L (Fin.succ i) (Fin.succ j)
      let b2 : Fin n → F := fun i => b (Fin.succ i) - L (Fin.succ i) 0 * b 0
      let tail := forwardSubstWithCost L2 b2
      ⟨Fin.cons (b 0) tail.value, tail.work + 2 * n⟩

Costed backward substitution. A level charges one pivot division plus the two field operations used to update each leading right-hand-side entry.

def backSubstWithCost : ∀ {n : Nat},
    Matrix (Fin n) (Fin n) F → (Fin n → F) → VectorExecution n F
  | 0, _U, _y => ⟨fun i => Fin.elim0 i, 0⟩
  | n + 1, U, y =>
      let last : Fin (n + 1) := Fin.last n
      let U1 : Matrix (Fin n) (Fin n) F := fun i j => U (Fin.castSucc i) (Fin.castSucc j)
      let xl : F := y last / U last last
      let y1 : Fin n → F := fun i => y (Fin.castSucc i) - U (Fin.castSucc i) last * xl
      let head := backSubstWithCost U1 y1
      ⟨Fin.snoc head.value xl, head.work + 2 * n + 1⟩

              theorem forwardSubstWithCost_value : ∀ {n : Nat}
    (L : Matrix (Fin n) (Fin n) F) (b : Fin n → F),
    (forwardSubstWithCost L b).value = forwardSubst L b
  | 0, _L, _b => by
      funext i
      exact Fin.elim0 i
  | n + 1, L, b => by
      simp [forwardSubstWithCost, forwardSubst, forwardSubstWithCost_value]theorem backSubstWithCost_value : ∀ {n : Nat}
    (U : Matrix (Fin n) (Fin n) F) (y : Fin n → F),
    (backSubstWithCost U y).value = backSubst U y
  | 0, _U, _y => by
      funext i
      exact Fin.elim0 i
  | n + 1, U, y => by
      simp [backSubstWithCost, backSubst, backSubstWithCost_value]theorem forwardSubstWithCost_work_le : ∀ {n : Nat}
    (L : Matrix (Fin n) (Fin n) F) (b : Fin n → F),
    (forwardSubstWithCost L b).work ≤ n ^ 2
  | 0, _L, _b => by simp [forwardSubstWithCost]
  | n + 1, L, b => by
      let L2 : Matrix (Fin n) (Fin n) F := fun i j => L (Fin.succ i) (Fin.succ j)
      let b2 : Fin n → F := fun i => b (Fin.succ i) - L (Fin.succ i) 0 * b 0
      have ih := forwardSubstWithCost_work_le L2 b2
      simp only [forwardSubstWithCost]
      change (forwardSubstWithCost L2 b2).work + 2 * n ≤ (n + 1) ^ 2
      simp [pow_two]
      nlinariththeorem backSubstWithCost_work_le : ∀ {n : Nat}
    (U : Matrix (Fin n) (Fin n) F) (y : Fin n → F),
    (backSubstWithCost U y).work ≤ n ^ 2
  | 0, _U, _y => by simp [backSubstWithCost]
  | n + 1, U, y => by
      let last : Fin (n + 1) := Fin.last n
      let U1 : Matrix (Fin n) (Fin n) F := fun i j => U (Fin.castSucc i) (Fin.castSucc j)
      let xl : F := y last / U last last
      let y1 : Fin n → F := fun i => y (Fin.castSucc i) - U (Fin.castSucc i) last * xl
      have ih := backSubstWithCost_work_le U1 y1
      simp only [backSubstWithCost]
      change (backSubstWithCost U1 y1).work + 2 * n + 1 ≤ (n + 1) ^ 2
      simp [pow_two]
      nlinarith

Apply a permutation to a vector by direct index lookup. This is the zero-field-operation implementation of multiplication by a permutation matrix.

def permuteVector {n : Nat} (σ : Equiv.Perm (Fin n)) (b : Fin n → F) : Fin n → F :=
  b ∘ σ

theorem permuteVector_eq_permMatrix_mulVec {n : Nat} (σ : Equiv.Perm (Fin n))
    (b : Fin n → F) :
    permuteVector σ b = σ.permMatrix F *ᵥ b := by
  simpa [permuteVector] using (Matrix.permMatrix_mulVec (σ := σ) (v := b)).symm

Costed LUP-SOLVE, including the direct permutation and both triangular substitutions.

def lupSolveWithCost {n : Nat} (σ : Equiv.Perm (Fin n))
    (L U : Matrix (Fin n) (Fin n) F) (b : Fin n → F) : VectorExecution n F :=
  let forward := forwardSubstWithCost L (permuteVector σ b)
  let backward := backSubstWithCost U forward.value
  ⟨backward.value, forward.work + backward.work⟩

Erasing the solver counter gives the existing lupSolve.


theorem lupSolveWithCost_value {n : Nat} (σ : Equiv.Perm (Fin n))
    (L U : Matrix (Fin n) (Fin n) F) (b : Fin n → F) :
    (lupSolveWithCost σ L U b).value = lupSolve σ L U b := by
  rw [lupSolveWithCost]
  simp only
  rw [backSubstWithCost_value, forwardSubstWithCost_value]
  rw [permuteVector_eq_permMatrix_mulVec]
  rfl

The actual two-substitution execution uses at most 2n² field operations.


theorem lupSolveWithCost_work_le {n : Nat} (σ : Equiv.Perm (Fin n))
    (L U : Matrix (Fin n) (Fin n) F) (b : Fin n → F) :
    (lupSolveWithCost σ L U b).work ≤ 2 * n ^ 2 := by
  rw [lupSolveWithCost]
  simp only
  have hf := forwardSubstWithCost_work_le L (permuteVector σ b)
  have hb := backSubstWithCost_work_le U
    (forwardSubstWithCost L (permuteVector σ b)).value
  omega

The costed solver returns the same certified solution as lupSolve.


theorem lupSolveWithCost_correct {n : Nat}
    {A L U : Matrix (Fin n) (Fin n) F} {σ : Equiv.Perm (Fin n)}
    (hLUP : σ.permMatrix F * A = L * U) (hL : IsUnitLowerTriangular L)
    (hU : IsUpperTriangular U) (hdiag : ∀ i : Fin n, U i i ≠ 0)
    (b : Fin n → F) :
    A *ᵥ (lupSolveWithCost σ L U b).value = b := by
  rw [lupSolveWithCost_value]
  exact lupSolve_correct hLUP hL hU hdiag b

end Chapter28end CLRS

Imports

import CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations

28.2. Inverting Matrices

This section formalizes the matrix-inversion consequence of the LUP decomposition (CLRS Theorem 28.2). Once A has an LUP factorization σ.permMatrix · A = L · U (Theorem 28.1), its inverse is obtained by inverting the triangular factors and undoing the row permutation: A⁻¹ = U⁻¹ · L⁻¹ · σ.permMatrix.

Main results:

Theorem inv_eq_lup: from σ.permMatrix · A = L · U (with σ.permMatrix the permutation matrix), A⁻¹ = U⁻¹ · L⁻¹ · σ.permMatrix.
Lemma permMatrix_inv: the inverse of a permutation matrix is the permutation matrix of the inverse permutation.
Lemma permMatrix_mul_inv: a permutation matrix left-multiplied by its inverse is the identity.

Notation conventions:

A : an n × n matrix over a field F.
σ.permMatrix F : the permutation matrix of σ.
M⁻¹ : the matrix inverse of M (zero for singular matrices).

              namespace CLRSnamespace Chapter28open Matrixvariable {F : Type} [Field F]

The inverse of a permutation matrix is the permutation matrix of the inverse permutation: (σ.permMatrix)⁻¹ = σ⁻¹.permMatrix.


lemma permMatrix_inv {n : ℕ} (σ : Equiv.Perm (Fin n)) :
    (σ.permMatrix F)⁻¹ = σ⁻¹.permMatrix F := by
  refine Matrix.inv_eq_right_inv ?_
  rw [← Matrix.permMatrix_mul]
  simp

A permutation matrix left-multiplied by its inverse is the identity.


lemma permMatrix_mul_inv {n : ℕ} (σ : Equiv.Perm (Fin n)) :
    (σ.permMatrix F)⁻¹ * σ.permMatrix F = 1 := by
  rw [permMatrix_inv]
  rw [← Matrix.permMatrix_mul]
  simp

Theorem 28.2 (inversion from the LUP decomposition). If σ.permMatrix · A = L · U is an LUP factorization, then A⁻¹ = U⁻¹ · L⁻¹ · σ.permMatrix: the inverse is obtained by inverting the upper factor, then the lower factor, then undoing the row permutation.

Proof: apply the inverse to both sides of the factorization and use Matrix.mul_inv_rev and the fact that a permutation matrix is invertible with inverse σ⁻¹.permMatrix.


theorem inv_eq_lup {n : ℕ} {A L U : Matrix (Fin n) (Fin n) F} {σ : Equiv.Perm (Fin n)}
    (hLUP : σ.permMatrix F * A = L * U) :
    A⁻¹ = U⁻¹ * L⁻¹ * σ.permMatrix F := by
  have hσ : (σ.permMatrix F)⁻¹ * σ.permMatrix F = 1 := permMatrix_mul_inv σ
  calc
    A⁻¹ = A⁻¹ * 1 := by rw [Matrix.mul_one]
    _ = A⁻¹ * ((σ.permMatrix F)⁻¹ * σ.permMatrix F) := by rw [hσ]
    _ = A⁻¹ * (σ.permMatrix F)⁻¹ * σ.permMatrix F := by rw [← Matrix.mul_assoc]
    _ = (σ.permMatrix F * A)⁻¹ * σ.permMatrix F := by rw [← Matrix.mul_inv_rev]
    _ = (L * U)⁻¹ * σ.permMatrix F := by rw [hLUP]
    _ = U⁻¹ * L⁻¹ * σ.permMatrix F := by rw [Matrix.mul_inv_rev]

end Chapter28end CLRS

Imports

import CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations

28.3. Symmetric Positive-Definite Matrices and Least Squares

This section formalizes CLRS §28.3: symmetric positive-definite (SPD) matrices, the Cholesky decomposition (Theorem 28.3), and least-squares approximation via the normal equations (Theorem 28.4).

Main results:

Definition IsSymPosDef: a real matrix is SPD if it is symmetric and xᵀAx > 0 for every nonzero vector x.
Theorem isSymPosDef_iff_posDef: SPD coincides with Mathlib's Matrix.PosDef, giving nonsingularity (IsSymPosDef.det_pos), positive diagonal (IsSymPosDef.diag_pos), and injectivity of mulVec (IsSymPosDef.mulVec_injective).
Theorem posDef_mul_transpose: if A has full column rank, then AᵀA is SPD (the Gram matrix is nonsingular).
Theorem cholesky_decomposition (Theorem 28.3): every SPD matrix factors as A = L·Lᵀ with L lower-triangular with positive diagonal (IsLowerTriangularPosDiag). The recursive construction uses the block factor L = [[√a, 0],[v/√a, L₂]] (choleskyFactor) with the positive-definite Schur complement (cholesky_schur_complement).
Theorem cholesky_unique: the Cholesky factor is unique — if L₁ and L₂ are lower-triangular with positive diagonal and L₁·L₁ᵀ = L₂·L₂ᵀ, then L₁ = L₂ (CLRS §28.3: the decomposition is unique).
Theorem normal_equations_minimizes (Theorem 28.4): if xh satisfies the normal equations Aᵀ(A·xh - b) = 0, then xh minimizes the squared residual (A·x - b)⬝ᵥ(A·x - b).
Theorem normal_equations_unique: when A has full column rank the minimizer is unique.
Theorem least_squares_closed_form: xh = (AᵀA)⁻¹·(Aᵀ·b) satisfies the normal equations — the closed-form least-squares solution.

This section now covers the whole of CLRS §28.3: the SPD foundations, the Cholesky decomposition (Theorem 28.3), and least-squares approximation (Theorem 28.4).

Notation conventions:

A : an m × n matrix over the reals; in the SPD sections, an n × n matrix.
v ⬝ᵥ w : the Euclidean dot product ∑ᵢ vᵢwᵢ; over the reals v ⬝ᵥ v is the squared 2-norm ‖v‖₂².
A *ᵥ x : matrix-vector product.

              namespace CLRSnamespace Chapter28open Matrixsection SymPosDef

A real matrix is symmetric positive-definite (SPD) when it is symmetric and xᵀAx > 0 for every nonzero vector x (CLRS §28.3). This is the CLRS definition; isSymPosDef_iff_posDef shows it coincides with Mathlib's Matrix.PosDef.

def IsSymPosDef {n : ℕ} (A : Matrix (Fin n) (Fin n) ℝ) : Prop :=
  A.IsSymm ∧ ∀ x : Fin n → ℝ, x ≠ 0 → 0 < x ⬝ᵥ (A *ᵥ x)

The CLRS notion of symmetric positive-definite coincides with Mathlib's Matrix.PosDef (Hermitian with xᴴAx > 0); over the reals xᴴ = x and Hermitian means symmetric.

theorem isSymPosDef_iff_posDef {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} :
    IsSymPosDef A ↔ A.PosDef := by
  constructor
  · intro hA
    apply Matrix.PosDef.of_dotProduct_mulVec_pos
    · simpa using hA.1
    · intro x hx
      exact hA.2 x hx
  · intro hA
    constructor
    · exact Matrix.isHermitian_iff_isSymm.mp hA.isHermitian
    · intro x hx
      simpa using hA.dotProduct_mulVec_pos hx

namespace IsSymPosDef

An SPD matrix is symmetric.

theorem isSymm {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) : A.IsSymm :=
  hA.1

An SPD matrix satisfies xᵀAx > 0 for every nonzero x.

theorem dotProduct_pos {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A)
    {x : Fin n → ℝ} (hx : x ≠ 0) : 0 < x ⬝ᵥ (A *ᵥ x) :=
  hA.2 x hx

An SPD matrix is positive definite in the Mathlib sense.

theorem posDef {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) : A.PosDef :=
  isSymPosDef_iff_posDef.mp hA

SPD matrices are nonsingular: the determinant is positive.

theorem det_pos {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) : 0 < A.det := by
  exact hA.posDef.det_pos

An SPD matrix has nonzero determinant.

theorem det_ne_zero {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) : A.det ≠ 0 :=
  hA.det_pos.ne'

Every diagonal entry of an SPD matrix is positive.

theorem diag_pos {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) (i : Fin n) :
    0 < A i i := by
  exact hA.posDef.diag_pos

An SPD matrix is invertible.

theorem isUnit {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) : IsUnit A := by
  exact hA.posDef.isUnit

An SPD matrix is injective on vectors: A·x = 0 implies x = 0.

theorem mulVec_injective {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : IsSymPosDef A) :
    Function.Injective A.mulVec := by
  exact A.mulVec_injective_iff_isUnit.2 hA.isUnit

end IsSymPosDef

AᵀA is symmetric for every A.


lemma mul_transpose_isSymm {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) : (Aᵀ * A).IsSymm := by
  rw [← Matrix.isHermitian_iff_isSymm]
  simpa using (Matrix.isHermitian_conjTranspose_mul_self A)

The adjoint identity over the reals: (A·v) ⬝ᵥ w = v ⬝ᵥ (Aᵀ·w).


lemma dot_mulVec_dotProduct {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (v : Fin n → ℝ)
    (w : Fin m → ℝ) : (A *ᵥ v) ⬝ᵥ w = v ⬝ᵥ (Aᵀ *ᵥ w) := by
  rw [← Matrix.vecMul_transpose]
  rw [Matrix.dotProduct_mulVec]

If A has full column rank (A.mulVec injective), then AᵀA is symmetric positive-definite. This is the nonsingularity of the Gram matrix that makes the normal equations uniquely solvable.


theorem posDef_mul_transpose {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (hA : Function.Injective A.mulVec) : IsSymPosDef (Aᵀ * A) := by
  constructor
  · exact mul_transpose_isSymm A
  · intro x hx
    rw [← Matrix.mulVec_mulVec]
    rw [← dot_mulVec_dotProduct A x (A *ᵥ x)]
    have hAx : A *ᵥ x ≠ 0 := by
      intro h
      exact hx (hA (by simpa using h))
    simpa using (Matrix.dotProduct_self_star_pos_iff).2 hAx

end SymPosDefsection LeastSquares

The squared Euclidean 2-norm of the residual A·x - b: the objective that least-squares approximation minimizes.

def residualSq {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ) (x : Fin n → ℝ) : ℝ :=
  (A *ᵥ x - b) ⬝ᵥ (A *ᵥ x - b)

If xh satisfies the normal equations, the residual A·xh - b is orthogonal to the column space of A.


lemma residual_orthogonal {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ)
    (xh : Fin n → ℝ) (h : Aᵀ *ᵥ (A *ᵥ xh - b) = 0) (x : Fin n → ℝ) :
    (A *ᵥ (x - xh)) ⬝ᵥ (A *ᵥ xh - b) = 0 := by
  rw [dot_mulVec_dotProduct A (x - xh) (A *ᵥ xh - b)]
  rw [h]
  simp

Pythagorean decomposition. When xh satisfies the normal equations, the squared residual at any x splits as the squared residual at xh plus the squared norm of A·(x - xh): ‖A·x - b‖₂² = ‖A·xh - b‖₂² + ‖A·(x - xh)‖₂².


lemma residual_sq_decomposition {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ)
    (xh : Fin n → ℝ) (h : Aᵀ *ᵥ (A *ᵥ xh - b) = 0) (x : Fin n → ℝ) :
    residualSq A b x = residualSq A b xh + (A *ᵥ (x - xh)) ⬝ᵥ (A *ᵥ (x - xh)) := by
  have hdec : ∀ r d : Fin m → ℝ,
      (r + d) ⬝ᵥ (r + d) = r ⬝ᵥ r + 2 * (r ⬝ᵥ d) + d ⬝ᵥ d := by
    intro r d
    rw [add_dotProduct, dotProduct_add, dotProduct_add, dotProduct_comm d r]
    ring
  have hx : A *ᵥ x - b = (A *ᵥ xh - b) + A *ᵥ (x - xh) := by
    rw [Matrix.mulVec_sub]
    abel
  rw [residualSq, hx]
  rw [hdec (A *ᵥ xh - b) (A *ᵥ (x - xh))]
  have hcross : (A *ᵥ xh - b) ⬝ᵥ (A *ᵥ (x - xh)) = 0 := by
    rw [dotProduct_comm]
    exact residual_orthogonal A b xh h x
  rw [hcross]
  rw [residualSq]
  ring

Theorem 28.4 (least-squares approximation). If xh satisfies the normal equations Aᵀ·(A·xh - b) = 0, then xh minimizes the squared residual ‖A·x - b‖₂² over all x.


theorem normal_equations_minimizes {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ)
    (xh : Fin n → ℝ) (h : Aᵀ *ᵥ (A *ᵥ xh - b) = 0) (x : Fin n → ℝ) :
    residualSq A b xh ≤ residualSq A b x := by
  rw [residual_sq_decomposition A b xh h x]
  have hnonneg : 0 ≤ (A *ᵥ (x - xh)) ⬝ᵥ (A *ᵥ (x - xh)) := by
    simpa using (dotProduct_self_star_nonneg (A *ᵥ (x - xh)))
  linarith

When A has full column rank, the minimizer of the squared residual is unique.


theorem normal_equations_unique {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ)
    (xh : Fin n → ℝ) (hA : Function.Injective A.mulVec)
    (h : Aᵀ *ᵥ (A *ᵥ xh - b) = 0) (x : Fin n → ℝ) (heq : residualSq A b x = residualSq A b xh) :
    x = xh := by
  have hdec := residual_sq_decomposition A b xh h x
  have hd : (A *ᵥ (x - xh)) ⬝ᵥ (A *ᵥ (x - xh)) = 0 := by
    linarith
  have hAx : A *ᵥ (x - xh) = 0 := by
    exact (dotProduct_self_star_eq_zero.mp (by simpa using hd))
  have hsub : x - xh = 0 := hA (by simpa using hAx)
  ext i
  have hi : x i - xh i = 0 := congr_fun hsub i
  linarith

Closed-form least-squares solution. When A has full column rank, xh = (AᵀA)⁻¹·(Aᵀ·b) satisfies the normal equations.

set_option linter.unnecessarySimpa false in


theorem least_squares_closed_form {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (b : Fin m → ℝ)
    (hA : Function.Injective A.mulVec) :
    Aᵀ *ᵥ (A *ᵥ ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b)) - b) = 0 := by
  have hPD : (Aᵀ * A).PosDef := by
    simpa using (Matrix.PosDef.conjTranspose_mul_self A hA)
  have hcancel : (Aᵀ * A) * (Aᵀ * A)⁻¹ = 1 := by
    letI := hPD.isUnit.invertible
    simpa using (Matrix.mul_inv_cancel_right_of_invertible (1 : Matrix (Fin n) (Fin n) ℝ))
  have hx : (Aᵀ * A) *ᵥ ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b)) = Aᵀ *ᵥ b := by
    rw [Matrix.mulVec_mulVec]
    rw [hcancel]
    simp
  calc
    Aᵀ *ᵥ (A *ᵥ ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b)) - b)
        = (Aᵀ * A) *ᵥ ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b)) - Aᵀ *ᵥ b := by
          rw [Matrix.mulVec_sub, Matrix.mulVec_mulVec]
    _ = Aᵀ *ᵥ b - Aᵀ *ᵥ b := by rw [hx]
    _ = 0 := by abel

The closed-form least-squares solution (AᵀA)⁻¹·Aᵀ·b minimizes the squared residual.

theorem least_squares_closed_form_minimizes {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (b : Fin m → ℝ) (hA : Function.Injective A.mulVec) (x : Fin n → ℝ) :
    residualSq A b ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b)) ≤ residualSq A b x := by
  exact normal_equations_minimizes A b ((Aᵀ * A)⁻¹ *ᵥ (Aᵀ *ᵥ b))
    (least_squares_closed_form A b hA) x

end LeastSquaressection Cholesky

The Schur complement of the leading 1×1 block of A: with a = A 0 0 and the first-column vector v i = A (Fin.succ i) 0, it is S = A₂₂ - (v·vᵀ)/a, the trailing block after one step of Gaussian elimination.

noncomputable def choleskySchur {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ) :
    Matrix (Fin n) (Fin n) ℝ :=
  fun i j => A (Fin.succ i) (Fin.succ j) - (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0

The dot product over Fin 1 is a product.


@[simp]
lemma dotProduct_const_fin_one (t w : ℝ) :
    (fun _ : Fin 1 => t) ⬝ᵥ (fun _ : Fin 1 => w) = t * w := by
  unfold dotProduct
  rw [Fin.sum_univ_succ]
  simp

Dot product distributes over a pointwise scalar multiple.


@[simp]
lemma dotProduct_add_scalar {n : ℕ} (y v : Fin n → ℝ) (t : ℝ) (w : Fin n → ℝ) :
    y ⬝ᵥ (fun i => v i * t + w i) = t * (v ⬝ᵥ y) + y ⬝ᵥ w := by
  unfold dotProduct
  rw [show (∑ i, y i * (v i * t + w i)) = ∑ i, (t * (y i * v i) + y i * w i) by
    apply Finset.sum_congr rfl
    intro i hi
    ring]
  rw [Finset.sum_add_distrib]
  rw [← Finset.mul_sum]
  congr 1
  rw [show (∑ i, y i * v i) = ∑ i, v i * y i by
    apply Finset.sum_congr rfl
    intro i hi
    ring]

Extracting a scalar factor from a dot product.


@[simp]
lemma dotProduct_smul_div {n : ℕ} (y v : Fin n → ℝ) (c : ℝ) :
    y ⬝ᵥ (fun i => v i * (v ⬝ᵥ y) / c) = (v ⬝ᵥ y) ^ 2 / c := by
  change (∑ i : Fin n, y i * (v i * (v ⬝ᵥ y) / c)) = (v ⬝ᵥ y) ^ 2 / c
  rw [show (∑ i, y i * (v i * (v ⬝ᵥ y) / c)) = (v ⬝ᵥ y) / c * (∑ i, v i * y i) by
    rw [Finset.mul_sum]
    apply Finset.sum_congr rfl
    intro i hi
    ring]
  rw [show (∑ i, v i * y i) = v ⬝ᵥ y by rfl]
  ring

The quadratic form of the reindexed vector z = u ∘ re equals the quadratic form of u under the reindexed matrix.


lemma reindex_quadratic {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (re : Fin (n + 1) ≃ Fin 1 ⊕ Fin n) (u : Fin 1 ⊕ Fin n → ℝ) :
    (u ∘ re) ⬝ᵥ (A *ᵥ (u ∘ re)) = u ⬝ᵥ ((A.reindex re re) *ᵥ u) := by
  have hmul : A *ᵥ (u ∘ re) = (A.reindex re re *ᵥ u) ∘ re := by
    funext i
    change (∑ j : Fin (n + 1), A i j * (u ∘ re) j) =
      ∑ s : Fin 1 ⊕ Fin n, (A.reindex re re) (re i) s * u s
    rw [← Equiv.sum_comp re.symm]
    apply Finset.sum_congr rfl
    intro s hs
    simp [Matrix.reindex_apply, Equiv.symm_apply_apply, Equiv.apply_symm_apply]
  rw [hmul]
  unfold dotProduct
  rw [← Equiv.sum_comp re.symm]
  apply Finset.sum_congr rfl
  intro s hs
  simp [Equiv.apply_symm_apply]

Reindexing A by finOneSumFin views it as the block matrix with entries [[A 0 0, A 0 (succ j)], [A (succ i) 0, A (succ i) (succ j)]].


lemma reindex_fromBlocks {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ) :
    A.reindex (finOneSumFin n) (finOneSumFin n) = Matrix.fromBlocks
      (fun _ _ : Fin 1 => A 0 0)
      (fun (_ : Fin 1) j => A 0 (Fin.succ j))
      (fun i (_ : Fin 1) => A (Fin.succ i) 0)
      (fun i j : Fin n => A (Fin.succ i) (Fin.succ j)) := by
  rw [Matrix.ext_iff_blocks]
  constructor
  · ext i j
    simp [Matrix.reindex_apply, finOneSumFin, Matrix.toBlocks₁₁]
  · constructor
    · ext i j
      simp [Matrix.reindex_apply, finOneSumFin, Matrix.toBlocks₁₂]
    · constructor
      · ext i j
        simp [Matrix.reindex_apply, finOneSumFin, Matrix.toBlocks₂₁]
      · ext i j
        simp [Matrix.reindex_apply, finOneSumFin, Matrix.toBlocks₂₂]

Block quadratic form. For z = (t, y) (the vector 0 ↦ t, succ i ↦ y i), zᵀAz expands as a·t² + 2·t·(v ⬝ᵥ y) + yᵀA₂₂y where a = A 0 0, v i = A (Fin.succ i) 0, and A₂₂ is the trailing block. This is the algebraic core of the Cholesky recursion.


lemma schur_quadratic_form {n : ℕ} {A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ}
    (hA : A.IsSymm) (t : ℝ) (y : Fin n → ℝ) :
    let a : ℝ := A 0 0
    let v : Fin n → ℝ := fun i => A (Fin.succ i) 0
    let A22 : Matrix (Fin n) (Fin n) ℝ := fun i j => A (Fin.succ i) (Fin.succ j)
    let z : Fin (n + 1) → ℝ := (Sum.elim (fun _ : Fin 1 => t) y) ∘ (finOneSumFin n)
    z ⬝ᵥ (A *ᵥ z) = a * t ^ 2 + 2 * t * (v ⬝ᵥ y) + y ⬝ᵥ (A22 *ᵥ y) := by
  intro a v A22 z
  dsimp [z]
  rw [reindex_quadratic A (finOneSumFin n) (Sum.elim (fun _ : Fin 1 => t) y)]
  rw [reindex_fromBlocks A]
  let A11 : Matrix (Fin 1) (Fin 1) ℝ := fun _ _ => A 0 0
  let A12 : Matrix (Fin 1) (Fin n) ℝ := fun _ j => A 0 (Fin.succ j)
  let A21 : Matrix (Fin n) (Fin 1) ℝ := fun i _ => A (Fin.succ i) 0
  let A22b : Matrix (Fin n) (Fin n) ℝ := fun i j => A (Fin.succ i) (Fin.succ j)
  rw [Matrix.fromBlocks_mulVec]
  rw [sumElim_dotProduct_sumElim]
  have hu_inl : (Sum.elim (fun _ : Fin 1 => t) y) ∘ Sum.inl = (fun _ : Fin 1 => t) := by
    funext x
    simp
  have hu_inr : (Sum.elim (fun _ : Fin 1 => t) y) ∘ Sum.inr = y := by
    funext x
    simp
  rw [hu_inl, hu_inr]
  have h11 : A11 *ᵥ (fun _ : Fin 1 => t) = fun _ : Fin 1 => a * t := by
    ext i
    simp [A11, Matrix.mulVec, a]
  have h12 : A12 *ᵥ y = fun _ : Fin 1 => v ⬝ᵥ y := by
    ext i
    simp [A12, Matrix.mulVec, v]
    congr 1
    funext j
    exact (Matrix.IsSymm.ext_iff.mp hA) (Fin.succ j) 0
  have h21 : A21 *ᵥ (fun _ : Fin 1 => t) = fun i : Fin n => v i * t := by
    ext i
    simp [A21, Matrix.mulVec, v]
  have h22 : A22b *ᵥ y = A22 *ᵥ y := by
    rfl
  rw [h11, h12, h21, h22]
  change (fun _ : Fin 1 => t) ⬝ᵥ (fun _ : Fin 1 => a * t + (v ⬝ᵥ y)) +
    y ⬝ᵥ (fun i : Fin n => v i * t + (A22 *ᵥ y) i) =
    a * t ^ 2 + 2 * t * (v ⬝ᵥ y) + y ⬝ᵥ (A22 *ᵥ y)
  rw [dotProduct_const_fin_one t (a * t + (v ⬝ᵥ y))]
  rw [dotProduct_add_scalar y v t (A22 *ᵥ y)]
  ring

Entrywise expansion of the Schur complement acting on y: (S·y) i = (A₂₂·y) i - v i·(v ⬝ᵥ y)/a.


lemma schur_mulVec {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ) (hA : A.IsSymm)
    (y : Fin n → ℝ) :
    (choleskySchur A *ᵥ y) = fun i =>
      ((fun j : Fin n => A (Fin.succ i) (Fin.succ j)) ⬝ᵥ y) -
        (A (Fin.succ i) 0 * ((fun j : Fin n => A (Fin.succ j) 0) ⬝ᵥ y)) / A 0 0 := by
  funext i
  unfold choleskySchur Matrix.mulVec dotProduct
  rw [show (∑ j, (A (Fin.succ i) (Fin.succ j) - (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0) * y j) =
      (∑ j, A (Fin.succ i) (Fin.succ j) * y j) -
        A (Fin.succ i) 0 * (∑ j, A (Fin.succ j) 0 * y j) / A 0 0 by
    rw [show (∑ j, (A (Fin.succ i) (Fin.succ j) - (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0) * y j) =
        ∑ j, (A (Fin.succ i) (Fin.succ j) * y j -
          (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0 * y j) by
      apply Finset.sum_congr rfl
      intro j hj
      ring]
    rw [Finset.sum_sub_distrib]
    congr 1
    rw [show (∑ j, (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0 * y j) =
          A (Fin.succ i) 0 * (∑ j, A (Fin.succ j) 0 * y j) / A 0 0 by
        rw [show (∑ j, (A (Fin.succ i) 0 * A 0 (Fin.succ j)) / A 0 0 * y j) =
            (A (Fin.succ i) 0 / A 0 0) * (∑ j, A 0 (Fin.succ j) * y j) by
          rw [Finset.mul_sum]
          apply Finset.sum_congr rfl
          intro j hj
          ring]
        rw [show (∑ j, A 0 (Fin.succ j) * y j) = ∑ j, A (Fin.succ j) 0 * y j by
          apply Finset.sum_congr rfl
          intro j hj
          exact congrArg (fun x => x * y j) (Matrix.IsSymm.ext_iff.mp hA (Fin.succ j) 0)]
        ring]]

The quadratic form of the Schur complement: yᵀSy = yᵀA₂₂y - (v ⬝ᵥ y)²/a.


lemma schur_residual {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ) (hA : A.IsSymm)
    (y : Fin n → ℝ) :
    y ⬝ᵥ (choleskySchur A *ᵥ y) =
      y ⬝ᵥ ((fun i j : Fin n => A (Fin.succ i) (Fin.succ j)) *ᵥ y) -
        ((fun i : Fin n => A (Fin.succ i) 0) ⬝ᵥ y) ^ 2 / A 0 0 := by
  rw [schur_mulVec A hA y]
  let W : Fin n → ℝ := fun i => (fun j : Fin n => A (Fin.succ i) (Fin.succ j)) ⬝ᵥ y
  let Z : Fin n → ℝ := fun i =>
    A (Fin.succ i) 0 * ((fun j : Fin n => A (Fin.succ j) 0) ⬝ᵥ y) / A 0 0
  change y ⬝ᵥ (W - Z) =
    y ⬝ᵥ ((fun i j : Fin n => A (Fin.succ i) (Fin.succ j)) *ᵥ y) -
      ((fun i : Fin n => A (Fin.succ i) 0) ⬝ᵥ y) ^ 2 / A 0 0
  rw [dotProduct_sub]
  congr 1
  rw [dotProduct_smul_div y (fun i => A (Fin.succ i) 0) (A 0 0)]

The Schur complement of an SPD matrix is SPD. This is the strict version of the PSD Matrix.PosSemidef.fromBlocks₁₁ fact, proved directly: for y ≠ 0, the vector z = (t, y) with t = -(v ⬝ᵥ y)/a satisfies zᵀAz = yᵀSy, so positivity of A transfers to S. This is the induction step of the Cholesky recursion (CLRS §28.3).


theorem cholesky_schur_complement {n : ℕ} {A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ}
    (hA : IsSymPosDef A) : IsSymPosDef (choleskySchur A) := by
  constructor
  · rw [Matrix.IsSymm.ext_iff]
    intro i j
    unfold choleskySchur
    have hs := Matrix.IsSymm.ext_iff.mp hA.isSymm
    rw [hs (Fin.succ j) (Fin.succ i), hs 0 (Fin.succ i), hs (Fin.succ j) 0]
    ring
  · intro y hy
    let a : ℝ := A 0 0
    let v : Fin n → ℝ := fun i => A (Fin.succ i) 0
    let A22 : Matrix (Fin n) (Fin n) ℝ := fun i j => A (Fin.succ i) (Fin.succ j)
    let S := choleskySchur A
    have ha : a ≠ 0 := ne_of_gt (hA.diag_pos 0)
    let t : ℝ := -(v ⬝ᵥ y) / a
    let z : Fin (n + 1) → ℝ := (Sum.elim (fun _ : Fin 1 => t) y) ∘ (finOneSumFin n)
    have hz : z ≠ 0 := by
      intro hz0
      apply hy
      funext i
      have hzi : z (Fin.succ i) = 0 := congr_fun hz0 (Fin.succ i)
      simpa [z, finOneSumFin] using hzi
    have ht2 : a * t ^ 2 + 2 * t * (v ⬝ᵥ y) = -(v ⬝ᵥ y) ^ 2 / a := by
      dsimp [t]
      field_simp [ha]
      ring
    have hres : y ⬝ᵥ (S *ᵥ y) = y ⬝ᵥ (A22 *ᵥ y) - (v ⬝ᵥ y) ^ 2 / a := by
      dsimp [S, A22, v, a]
      rw [schur_residual A hA.isSymm y]
    have hquad : z ⬝ᵥ (A *ᵥ z) = a * t ^ 2 + 2 * t * (v ⬝ᵥ y) + y ⬝ᵥ (A22 *ᵥ y) := by
      dsimp [a, v, A22, t, z]
      exact schur_quadratic_form hA.isSymm t y
    have hexp : z ⬝ᵥ (A *ᵥ z) = y ⬝ᵥ (S *ᵥ y) := by
      rw [hquad, ht2]
      rw [hres]
      ring
    have hpos : 0 < z ⬝ᵥ (A *ᵥ z) := hA.dotProduct_pos hz
    rw [hexp] at hpos
    exact hpos

The Cholesky factor for the recursion step: with a = A 0 0 and the first-column vector v i = A (Fin.succ i) 0, choleskyFactor A L2 is the block lower-triangular matrix [[√a, 0], [v/√a, L2]] over Fin (n+1) built from the scalar √a, the column v/√a, and the trailing factor L2.

noncomputable def choleskyFactor {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ :=
  fun i j =>
    Fin.cases
      (fun j : Fin (n + 1) => if j = 0 then Real.sqrt (A 0 0) else 0)
      (fun i' : Fin n => fun j : Fin (n + 1) =>
        Fin.cases (A (Fin.succ i') 0 / Real.sqrt (A 0 0)) (fun j' : Fin n => L2 i' j') j)
      i j

The top-left entry of the Cholesky factor is √a.

@[simp]
lemma choleskyFactor_00 {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) : choleskyFactor A L2 0 0 = Real.sqrt (A 0 0) := by
  simp [choleskyFactor]

The first row of the Cholesky factor (right of the diagonal) is zero.

@[simp]
lemma choleskyFactor_0_succ {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (j : Fin n) : choleskyFactor A L2 0 (Fin.succ j) = 0 := by
  simp [choleskyFactor]

The first column of the Cholesky factor (below the diagonal) is v/√a.

@[simp]
lemma choleskyFactor_succ_0 {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (i : Fin n) :
    choleskyFactor A L2 (Fin.succ i) 0 = A (Fin.succ i) 0 / Real.sqrt (A 0 0) := by
  simp [choleskyFactor]

The trailing block of the Cholesky factor is L2.

@[simp]
lemma choleskyFactor_succ_succ {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (i j : Fin n) :
    choleskyFactor A L2 (Fin.succ i) (Fin.succ j) = L2 i j := by
  simp [choleskyFactor]

L is lower-triangular with positive diagonal: every entry above the main diagonal is zero and every diagonal entry is strictly positive.

def IsLowerTriangularPosDiag {n : ℕ} (L : Matrix (Fin n) (Fin n) ℝ) : Prop :=
  (∀ ⦃i j : Fin n⦄, i < j → L i j = 0) ∧ (∀ i : Fin n, 0 < L i i)

The (0,0) entry of L·Lᵀ for the Cholesky factor is a = A 0 0.


lemma choleskyFactor_mul_00 {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (hpos : 0 < A 0 0) :
    (choleskyFactor A L2 * (choleskyFactor A L2)ᵀ) 0 0 = A 0 0 := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp
  simpa [pow_two] using (Real.sq_sqrt hpos.le)

The (0, j) entries of L·Lᵀ for the Cholesky factor match A.


lemma choleskyFactor_mul_0_succ {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (hA : A.IsSymm) (hpos : 0 < A 0 0) (j : Fin n) :
    (choleskyFactor A L2 * (choleskyFactor A L2)ᵀ) 0 (Fin.succ j) = A 0 (Fin.succ j) := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp
  have hs_ne : Real.sqrt (A 0 0) ≠ 0 := (Real.sqrt_pos.2 hpos).ne'
  have hv : A 0 (Fin.succ j) = A (Fin.succ j) 0 := by
    exact (Matrix.IsSymm.ext_iff.mp hA) (Fin.succ j) 0
  rw [hv]
  field_simp [hs_ne]

The (i, 0) entries of L·Lᵀ for the Cholesky factor match A.


lemma choleskyFactor_mul_succ_0 {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (hpos : 0 < A 0 0) (i : Fin n) :
    (choleskyFactor A L2 * (choleskyFactor A L2)ᵀ) (Fin.succ i) 0 = A (Fin.succ i) 0 := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp
  have hs_ne : Real.sqrt (A 0 0) ≠ 0 := (Real.sqrt_pos.2 hpos).ne'
  field_simp [hs_ne]

Entrywise expansion of the trailing block of L·Lᵀ: (L·Lᵀ)ᵢⱼ = (vᵢ/√a)·(vⱼ/√a) + (L2·L2ᵀ)ᵢⱼ.


lemma choleskyFactor_mul_succ_succ {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (i j : Fin n) :
    (choleskyFactor A L2 * (choleskyFactor A L2)ᵀ) (Fin.succ i) (Fin.succ j) =
      (A (Fin.succ i) 0 / Real.sqrt (A 0 0)) * (A (Fin.succ j) 0 / Real.sqrt (A 0 0)) +
        (L2 * L2ᵀ) i j := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp [Matrix.mul_apply, Matrix.transpose_apply]

When L2·L2ᵀ is the Schur complement, the trailing block of L·Lᵀ matches the trailing block of A.


lemma choleskyFactor_mul_succ_succ_eq {n : ℕ} (A : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (L2 : Matrix (Fin n) (Fin n) ℝ) (hA : A.IsSymm) (hpos : 0 < A 0 0) (i j : Fin n)
    (hL2 : choleskySchur A = L2 * L2ᵀ) :
    (choleskyFactor A L2 * (choleskyFactor A L2)ᵀ) (Fin.succ i) (Fin.succ j) = A (Fin.succ i) (Fin.succ j) := by
  rw [choleskyFactor_mul_succ_succ A L2 i j]
  rw [← hL2]
  unfold choleskySchur
  have hsym : A 0 (Fin.succ j) = A (Fin.succ j) 0 := by
    exact (Matrix.IsSymm.ext_iff.mp hA) (Fin.succ j) 0
  rw [hsym]
  have hs_ne : Real.sqrt (A 0 0) ≠ 0 := (Real.sqrt_pos.2 hpos).ne'
  have hs_sq : (Real.sqrt (A 0 0)) ^ 2 = A 0 0 := Real.sq_sqrt hpos.le
  field_simp [hs_ne, hpos.ne']
  ring_nf
  rw [hs_sq]
  ring

Cholesky decomposition (CLRS Theorem 28.3). Every symmetric positive-definite matrix A factors as A = L·Lᵀ with L lower-triangular with positive diagonal.

The proof is the block recursion: writing a = A 0 0, v i = A (Fin.succ i) 0, and S for the Schur complement A₂₂ - v·vᵀ/a, the strict positivity of A makes S SPD (cholesky_schur_complement), so by induction S = L₂·L₂ᵀ. The factor L = [[√a, 0],[v/√a, L₂]] (choleskyFactor) is lower-triangular with positive diagonal and satisfies L·Lᵀ = A.


theorem cholesky_decomposition {n : ℕ} (A : Matrix (Fin n) (Fin n) ℝ) (hA : IsSymPosDef A) :
    ∃ L : Matrix (Fin n) (Fin n) ℝ, IsLowerTriangularPosDiag L ∧ A = L * Lᵀ := by
  induction n with
  | zero =>
      refine ⟨1, ?tri, ?eq⟩
      · constructor
        · intro i j hij
          exact Fin.elim0 i
        · intro i
          exact Fin.elim0 i
      · ext i j
        exact Fin.elim0 i
  | succ n ih =>
      let a : ℝ := A 0 0
      let S : Matrix (Fin n) (Fin n) ℝ := choleskySchur A
      have ha : 0 < a := by simpa [a] using hA.diag_pos 0
      have hS : IsSymPosDef S := by simpa [S] using cholesky_schur_complement hA
      rcases ih S hS with ⟨L2, ⟨hL2lt, hL2pd⟩, hL2eq⟩
      let L : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ := choleskyFactor A L2
      refine ⟨L, ?tri', ?eq'⟩
      · constructor
        · intro i j hij
          rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨i', rfl⟩
          · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
            · simp at hij
            · simp [L]
          · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
            · simp at hij
            · have hij' : i' < j' := by simpa using hij
              simpa [L] using (hL2lt hij')
        · intro i
          rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨i', rfl⟩
          · simpa [L, a] using (Real.sqrt_pos.2 ha)
          · simpa [L] using (hL2pd i')
      · ext i j
        rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨i', rfl⟩
        · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
          · simpa [L, a] using (choleskyFactor_mul_00 A L2 ha).symm
          · simpa [L] using (choleskyFactor_mul_0_succ A L2 hA.isSymm ha j').symm
        · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
          · simpa [L] using (choleskyFactor_mul_succ_0 A L2 ha i').symm
          · simpa [L] using (choleskyFactor_mul_succ_succ_eq A L2 hA.isSymm ha i' j' hL2eq).symm

The trailing block L₂ i j = L (Fin.succ i) (Fin.succ j) of a matrix L (its restriction to rows and columns 1..n).

def trailingBlock {n : ℕ} (L : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ) : Matrix (Fin n) (Fin n) ℝ :=
  fun i j => L (Fin.succ i) (Fin.succ j)

The (0,0) entry of L·Lᵀ for a lower-triangular L is L 0 0 · L 0 0 (the 1×1 block).


lemma lowerTri_mul_transpose_00 {n : ℕ} (L : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (hL : IsLowerTriangular L) : (L * Lᵀ) 0 0 = L 0 0 * L 0 0 := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp [Matrix.transpose_apply]
  have h0k : ∀ k : Fin n, L 0 (Fin.succ k) = 0 := fun k => hL (Fin.succ_pos k)
  simp [h0k]

The (i,0) entry of L·Lᵀ for a lower-triangular L is L (succ i) 0 · L 0 0 (the first column of L against the 1×1 block).


lemma lowerTri_mul_transpose_succ_0 {n : ℕ} (L : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (hL : IsLowerTriangular L) (i : Fin n) : (L * Lᵀ) (Fin.succ i) 0 = L (Fin.succ i) 0 * L 0 0 := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp [Matrix.transpose_apply]
  have h0k : ∀ k : Fin n, L 0 (Fin.succ k) = 0 := fun k => hL (Fin.succ_pos k)
  simp [h0k]

The trailing (i,j) entry of L·Lᵀ for a lower-triangular L splits as L (succ i) 0 · L (succ j) 0 + (L₂·L₂ᵀ) i j with L₂ the trailing block.


lemma lowerTri_mul_transpose_succ_succ {n : ℕ} (L : Matrix (Fin (n + 1)) (Fin (n + 1)) ℝ)
    (Variable name `hL` is not explicitly referenced.

The binding can be removed (if unused) or named `_` (if used implicitly).

Note: This linter can be disabled with `set_option linter.unusedVariables false`hL : IsLowerTriangular L) (i j : Fin n) :
    (L * Lᵀ) (Fin.succ i) (Fin.succ j) =
      L (Fin.succ i) 0 * L (Fin.succ j) 0 + (trailingBlock L * (trailingBlock L)ᵀ) i j := by
  rw [Matrix.mul_apply]
  rw [Fin.sum_univ_succ]
  simp [Matrix.transpose_apply]
  unfold trailingBlock
  simp [Matrix.mul_apply, Matrix.transpose_apply]

Uniqueness of the Cholesky decomposition. If L₁ and L₂ are both lower-triangular with positive diagonal and L₁·L₁ᵀ = L₂·L₂ᵀ, then L₁ = L₂.

The proof is the block recursion of cholesky_decomposition: the (0,0) entries give L₁ 0 0 = L₂ 0 0 (equal squares, both positive), the first column then matches, and the trailing blocks satisfy L₁₂·L₁₂ᵀ = L₂₂·L₂₂ᵀ, so the induction hypothesis applies.


theorem cholesky_unique {n : ℕ} (L1 L2 : Matrix (Fin n) (Fin n) ℝ)
    (h1 : IsLowerTriangularPosDiag L1) (h2 : IsLowerTriangularPosDiag L2)
    (h : L1 * L1ᵀ = L2 * L2ᵀ) : L1 = L2 := by
  induction n with
  | zero =>
      ext i
      exact Fin.elim0 i
  | succ n ih =>
      have h00 : L1 0 0 = L2 0 0 := by
        have hsq : L1 0 0 * L1 0 0 = L2 0 0 * L2 0 0 := by
          calc
            L1 0 0 * L1 0 0 = (L1 * L1ᵀ) 0 0 := by rw [lowerTri_mul_transpose_00 L1 h1.1]
            _ = (L2 * L2ᵀ) 0 0 := by rw [h]
            _ = L2 0 0 * L2 0 0 := by rw [lowerTri_mul_transpose_00 L2 h2.1]
        have hsq2 : L1 0 0 ^ 2 = L2 0 0 ^ 2 := by simpa [pow_two] using hsq
        have hp1 : 0 ≤ L1 0 0 := le_of_lt (h1.2 0)
        have hp2 : 0 ≤ L2 0 0 := le_of_lt (h2.2 0)
        rcases (sq_eq_sq_iff_eq_or_eq_neg.mp hsq2) with e | e
        · exact e
        · nlinarith
      have hcol : ∀ i : Fin n, L1 (Fin.succ i) 0 = L2 (Fin.succ i) 0 := by
        intro i
        have hc : L1 (Fin.succ i) 0 * L1 0 0 = L2 (Fin.succ i) 0 * L2 0 0 := by
          calc
            L1 (Fin.succ i) 0 * L1 0 0 = (L1 * L1ᵀ) (Fin.succ i) 0 := by rw [lowerTri_mul_transpose_succ_0 L1 h1.1 i]
            _ = (L2 * L2ᵀ) (Fin.succ i) 0 := by rw [h]
            _ = L2 (Fin.succ i) 0 * L2 0 0 := by rw [lowerTri_mul_transpose_succ_0 L2 h2.1 i]
        rw [h00] at hc
        exact mul_right_cancel₀ (h2.2 0).ne' hc
      have htrail : trailingBlock L1 * (trailingBlock L1)ᵀ = trailingBlock L2 * (trailingBlock L2)ᵀ := by
        ext i j
        have hfull : (L1 * L1ᵀ) (Fin.succ i) (Fin.succ j) = (L2 * L2ᵀ) (Fin.succ i) (Fin.succ j) := by
          simpa using congr_fun (congr_fun h (Fin.succ i)) (Fin.succ j)
        rw [lowerTri_mul_transpose_succ_succ L1 h1.1 i j] at hfull
        rw [lowerTri_mul_transpose_succ_succ L2 h2.1 i j] at hfull
        rw [hcol i, hcol j] at hfull
        linarith
      have htrail1 : IsLowerTriangularPosDiag (trailingBlock L1) := by
        constructor
        · intro i j hij
          exact h1.1 (Fin.succ_lt_succ_iff.mpr hij)
        · intro i
          exact h1.2 (Fin.succ i)
      have htrail2 : IsLowerTriangularPosDiag (trailingBlock L2) := by
        constructor
        · intro i j hij
          exact h2.1 (Fin.succ_lt_succ_iff.mpr hij)
        · intro i
          exact h2.2 (Fin.succ i)
      have heq_trail : trailingBlock L1 = trailingBlock L2 :=
        ih (trailingBlock L1) (trailingBlock L2) htrail1 htrail2 htrail
      ext i j
      rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨i', rfl⟩
      · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
        · exact h00
        · exact (h1.1 (Fin.succ_pos j')).trans (h2.1 (Fin.succ_pos j')).symm
      · rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨j', rfl⟩
        · exact hcol i'
        · exact congr_fun (congr_fun heq_trail i') j'

              end Choleskyend Chapter28end CLRS

Scope and implementation notes

Imports

import CLRSLean.Chapter_28
import CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations
import CLRSLean.FourthEdition.Chapter_28.Section_28_1_Linear_Equations.ExecutableLUP
import CLRSLean.FourthEdition.Chapter_28.Section_28_2_Inverting_Matrices
import CLRSLean.FourthEdition.Chapter_28.Section_28_3_Symmetric_Positive_Definite

Current source

Sections 28.1--28.3 are native fourth-edition sections (solving systems of linear equations via LUP decomposition, inverting matrices, and symmetric positive-definite matrices with least-squares approximation), imported directly from Section 28.1, Section 28.2, and Section 28.3. Declarations keep their current namespaces; the third-edition-numbered imports CLRSLean.Chapter_28 and CLRSLean.Chapter_28.Section_28_* forward to these sources.

Coverage boundary

The native sections supply the represented fourth-edition matrix-operations sections (CLRS Theorem 28.1 and Lemmas 28.1--28.2). The constructive Theorem 28.1 layer is exposed by lupDecomposeWithCost: it explicitly scans for a nonzero pivot, swaps rows, performs pointwise Gaussian elimination, recurses on the Schur block, directly reindexes permutation rows, counts the factor-assembly multiplier divisions, and returns either factors or failure. lupDecomposeWithCost_correct proves that every nonsingular input returns a unit-lower-triangular factor, an upper-triangular factor with nonzero diagonal, and the exact equation P·A=L·U, while the same execution's counter is at most 4n³. Failure is equivalent to a zero determinant.

The costed solver lupSolveWithCost erases to the existing lupSolve, inherits its solution theorem, and records at most 2n² field operations; its permutation is implemented by direct vector indexing rather than a hidden matrix-vector product. These are exact-field unit-cost results with decidable zero testing; floating-point stability, mutable storage, allocation, and bit/RAM costs remain outside this boundary. The older numerical substitution, decomposition, inversion and Cholesky budgets have separate isBigO proofs; those statements are upper bounds. Only the costed LUP decomposition/solve results described above are attached to those executions.

See docs/clrs-fourth-edition-map.csv for the section-level mapping and docs/migrations/clrs4.md for compatibility and deprecation policy.

CLRS, fourth edition · Chapter 28 of 35