Prove using Cayley-Hamilton

Stream: new members

[Nick_adfor]

I want to prove the following theorem in lean 4, leaving some sorry wondering how to fill in:

Theorem:

We are given that B is an invertible n × n matrix over a field F, and its characteristic polynomial is:

p_B(X) = Xⁿ + cₙ₋₁ Xⁿ⁻¹ + ... + c₁ X + c₀.

We are to prove that B⁻¹ can be written as a linear combination of Iₙ, B, ..., Bⁿ⁻¹, specifically:

B⁻¹ = b₀ Iₙ + b₁ B + b₂ B² + ... + bₙ₋₁ Bⁿ⁻¹,

where bᵢ = -c_{i+1} / c₀ for i = 0 to n−1.

import Mathlib

variable {F : Type*} [Field F] {n : Type*} [DecidableEq n] [Fintype n]

theorem inverse_eq_poly (B : Matrix n n F) (hB : IsUnit B) :
    ∃ (b : Fin (Fintype.card n) → F), B⁻¹ = Finset.sum Finset.univ (fun i : Fin (Fintype.card n) => b i • (B ^ (i.val))) := by
  let p := Matrix.charpoly B
  have hdeg : Polynomial.natDegree p = Fintype.card n := Matrix.charpoly_natDegree_eq_dim B
  have hpoly : p = Finset.sum (Finset.range (Fintype.card n + 1)) (fun i => Polynomial.monomial i (Polynomial.coeff p i)) := by
    refine Polynomial.as_sum_range' p (Fintype.card n + 1) ?_
    rw [hdeg, @Nat.lt_add_right_iff_pos, @Nat.lt_one_iff]
  have hCH : Polynomial.aeval B p = 0 := Matrix.aeval_self_charpoly B
  have h_expanded :
      B^(Fintype.card n) +
      (∑ i ∈ Finset.range (Fintype.card n), (Polynomial.coeff p i) • B ^ i) = 0 := by
    rw [← Polynomial.sum_monomial_eq p]
    have : Polynomial.aeval B p =
        Polynomial.aeval B (Finset.sum (Finset.range (Fintype.card n + 1))
                              (fun i => Polynomial.monomial i (Polynomial.coeff p i))) := by
      rw [hpoly]
      sorry
    sorry
  have hc0 : Polynomial.coeff p 0 ≠ 0 := by
    have hdet := Matrix.det_eq_sign_charpoly_coeff B
    have hcoeff : (-1)^(Fintype.card n) * p.coeff 0 = B.det := by
      rw [hdet, mul_comm]
    have hdet_unit : IsUnit B.det := by
      exact (Matrix.isUnit_iff_isUnit_det B).mp hB
    by_contra h
    rw [h, mul_zero] at hcoeff
    have : B.det = 0 := by
      rw [← hcoeff]
    exact hdet_unit.ne_zero this
  have hI : (1 : Matrix n n F) = (- (Polynomial.coeff p 0)⁻¹) • (Finset.sum (Finset.range (Fintype.card n)) (fun i => Polynomial.coeff p (i + 1) • (B ^ (i + 1)))) := by
    sorry
  have hBinv : B⁻¹ = (- (Polynomial.coeff p 0)⁻¹) • (Finset.sum (Finset.range (Fintype.card n)) (fun i => Polynomial.coeff p (i + 1) • (B ^ i))) := by
    rw [← Matrix.mul_one B⁻¹, hI, Matrix.mul_smul, Matrix.mul_sum]
    sorry
  let b (i : Fin (Fintype.card n)) := - (Polynomial.coeff p (i.val + 1)) * (Polynomial.coeff p 0)⁻¹
  use b
  rw [hBinv]
  simp only [Fin.sum_univ_eq_sum_range, smul_smul, mul_comm _ (Polynomial.coeff p 0)⁻¹, b]
  ext i
  sorry

[Aaron Liu]

For the first sorry, you can do

    have : Polynomial.aeval B p =
        Polynomial.aeval B (Finset.sum (Finset.range (Fintype.card n + 1))
                              (fun i => Polynomial.monomial i (Polynomial.coeff p i))) := by
      rw [← hpoly]

[Nick_adfor]

Aaron Liu said:

← hpoly

Thanks for your correction! I even make a mistake using the wrong side of hpoly : (

[Nick_adfor]

Now this

import Mathlib

variable {F : Type*} [Field F] {n : Type*} [DecidableEq n] [Fintype n][Nonempty n]

lemma coeff_sum_monomial_eq (p : Polynomial F) (i : ℕ) :
  (p.sum (fun n a => Polynomial.monomial n a)).coeff i = p.coeff i := by simp

theorem inverse_eq_poly (B : Matrix n n F) (hB : IsUnit B) :
    ∃ (b : Fin (Fintype.card n) → F), B⁻¹ = Finset.sum Finset.univ (fun i : Fin (Fintype.card n) => b i • (B ^ (i.val))) := by
  let p := Matrix.charpoly B
  have hdeg : Polynomial.natDegree p = Fintype.card n := Matrix.charpoly_natDegree_eq_dim B
  have hpoly : p = Finset.sum (Finset.range (Fintype.card n + 1)) (fun i => Polynomial.monomial i (Polynomial.coeff p i)) := by
    refine Polynomial.as_sum_range' p (Fintype.card n + 1) ?_
    rw [hdeg, @Nat.lt_add_right_iff_pos, @Nat.lt_one_iff]
  have hCH : Polynomial.aeval B p = 0 := Matrix.aeval_self_charpoly B
  have h_expanded :
      B^(Fintype.card n) +
      (∑ i ∈ Finset.range (Fintype.card n), (Polynomial.coeff p i) • B ^ i) = 0 := by
    rw [← Polynomial.sum_monomial_eq p]
    rw [hpoly] at hCH
    simp only [Polynomial.aeval_def] at hCH
    have h_coeff_eq := fun i => coeff_sum_monomial_eq p i
    have eq_sum : ∑ i ∈ Finset.range (Fintype.card n), (p.sum fun n a => Polynomial.monomial n a).coeff i • B ^ i
               = ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i := by
      rw [@Polynomial.sum_monomial_eq]
    rw [eq_sum]
    have sum_decomp :
      ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i =
      p.coeff 0 • B^0 + ∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B^(i + 1) := by
        have h_card : Fintype.card n = (Fintype.card n - 1).succ := by
          rw [Nat.succ_eq_add_one]
          have hpos : 0 < Fintype.card n := by
            exact Fintype.card_pos
          have h := Nat.succ_pred_eq_of_pos hpos
          exact id (Eq.symm h)
        rw [h_card, Finset.sum_range_succ']
        exact
          Eq.symm
            (AddCommMagma.add_comm (p.coeff 0 • B ^ 0)
              (∑ i ∈ Finset.range ((Fintype.card n - 1).succ - 1), p.coeff (i + 1) • B ^ (i + 1)))
    rw [sum_decomp]
    simp only [pow_zero]
    rw [← add_zero (B^(Fintype.card n)), add_assoc, ← neg_eq_iff_add_eq_zero, zero_add]
    rw [neg_eq_iff_eq_neg, neg_add]
    rw [add_comm] at hCH
    have h_eval_expanded : Polynomial.eval₂ (algebraMap F (Matrix n n F)) B
        (∑ i ∈ Finset.range (1 + Fintype.card n), Polynomial.monomial i (p.coeff i)) =
        ∑ i ∈ Finset.range (1 + Fintype.card n),
          (Polynomial.eval₂ (algebraMap F (Matrix n n F)) B (Polynomial.monomial i (p.coeff i))) :=
      by rw [@Polynomial.eval₂_finset_sum]
    have h_mono_eval : ∀ i, Polynomial.eval₂ (algebraMap F (Matrix n n F)) B
        (Polynomial.monomial i (p.coeff i)) = (algebraMap F (Matrix n n F)) (p.coeff i) * B^i := by
      intro i
      rw [Polynomial.eval₂_monomial]
    rw [h_eval_expanded] at hCH
    simp only [h_mono_eval, Algebra.algebraMap_eq_smul_one] at hCH
    have h_range_split : Finset.range (1 + Fintype.card n) =
        Finset.range (Fintype.card n + 1) := by rw [add_comm]
    rw [h_range_split, Finset.sum_range_succ] at hCH
    simp only [Algebra.smul_mul_assoc, one_mul] at hCH
    have h_leading_coeff : p.coeff (Fintype.card n) = 1 := by
      rw [←Matrix.charpoly_natDegree_eq_dim B]
      exact (Matrix.charpoly_monic B).coeff_natDegree
    rw [h_leading_coeff, one_smul] at hCH
    have h_eq : B^(Fintype.card n) = - ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i := by
      exact Eq.symm (neg_eq_of_add_eq_zero_right hCH)
    have h_split :
      ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i =
        p.coeff 0 • B^0 + ∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1) := by
      rw [sum_decomp]
    rw [h_split] at h_eq
    rw [pow_zero] at h_eq
    rw [h_eq]
    exact
      neg_add (p.coeff 0 • 1)
        (∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1))
  have hc0 : Polynomial.coeff p 0 ≠ 0 := by
      have hdet := Matrix.det_eq_sign_charpoly_coeff B
      have hcoeff : (-1)^(Fintype.card n) * p.coeff 0 = B.det := by
        rw [hdet, mul_comm]
      have hdet_unit : IsUnit B.det := by
        exact (Matrix.isUnit_iff_isUnit_det B).mp hB
      by_contra h
      rw [h, mul_zero] at hcoeff
      have : B.det = 0 := by
        rw [← hcoeff]
      exact hdet_unit.ne_zero this
  have hI :
    (1 : Matrix n n F) =
      (- (Polynomial.coeff p 0)⁻¹) •
        (Finset.sum (Finset.range (Fintype.card n)) fun i =>
          Polynomial.coeff p (i + 1) • B ^ (i + 1)) := by
              sorry
  have hBinv : B⁻¹ = (- (Polynomial.coeff p 0)⁻¹) • (Finset.sum (Finset.range (Fintype.card n)) (fun i => Polynomial.coeff p (i + 1) • (B ^ i))) := by
    rw [← Matrix.mul_one B⁻¹, hI, Matrix.mul_smul, Matrix.mul_sum]
    congr
    ext i
    rw [Matrix.mul_smul, pow_succ']
    sorry
  let b (i : Fin (Fintype.card n)) := - (Polynomial.coeff p (i.val + 1)) * (Polynomial.coeff p 0)⁻¹
  use b
  rw [hBinv]
  simp only [mul_comm _ (Polynomial.coeff p 0)⁻¹, b]
  ext i j
  simp only [Matrix.smul_apply, Finset.sum_apply, smul_eq_mul]
  let ι := Fin (Fintype.card n)
  have h_equiv : Finset.range (Fintype.card n) = Finset.univ.image (fun (x : ι) => x.val) := by
    apply Finset.ext
    intro x
    simp only [Finset.mem_range, Finset.mem_image, Finset.mem_univ, true_and]
    constructor
    · intro hx
      exact CanLift.prf x hx
    · rintro ⟨y, _, rfl⟩
      exact y.is_lt
  rw [h_equiv, Finset.sum_image]
  · sorry
  · intro x _ y _ h
    exact Fin.ext h

[Nick_adfor]

Finished

import Mathlib.Algebra.Lie.OfAssociative
import Mathlib.LinearAlgebra.Matrix.Charpoly.Coeff
import Mathlib.RingTheory.Henselian

-- Declare variables: F is a field, n is a finite nonempty type with decidable equality
variable {F : Type*} [Field F] {n : Type*} [DecidableEq n] [Fintype n] [Nonempty n]

/-- `coeff_sum_monomial_eq`: For any polynomial `p`, summing its monomials explicitly and extracting the i-th coefficient yields the same result as `p.coeff i`. -/
lemma coeff_sum_monomial_eq (p : Polynomial F) (i : ℕ) :
  (p.sum (fun n a => Polynomial.monomial n a)).coeff i = p.coeff i := by simp

/-- `inverse_eq_poly`: If `B` is an invertible matrix over a field `F`, then its inverse can be expressed as a linear combination of powers of `B` with coefficients in `F`. -/
theorem inverse_eq_poly (B : Matrix n n F) [Invertible B] (hB : IsUnit B) :
    ∃ (b : Fin (Fintype.card n) → F), B⁻¹ = Finset.sum Finset.univ (fun i : Fin (Fintype.card n) => b i • (B ^ (i.val))) := by
  let p := Matrix.charpoly B  -- Let p be the characteristic polynomial of B

  -- The degree of the characteristic polynomial equals the dimension of the matrix
  have hdeg : Polynomial.natDegree p = Fintype.card n := Matrix.charpoly_natDegree_eq_dim B

  -- Express the polynomial `p` as a sum of monomials explicitly
  have hpoly : p = Finset.sum (Finset.range (Fintype.card n + 1)) (fun i => Polynomial.monomial i (Polynomial.coeff p i)) := by
    refine Polynomial.as_sum_range' p (Fintype.card n + 1) ?_
    rw [hdeg, Nat.lt_add_right_iff_pos, Nat.lt_one_iff]

  -- Cayley-Hamilton: `p` evaluated at `B` gives zero
  have hCH : Polynomial.aeval B p = 0 := Matrix.aeval_self_charpoly B

  -- Convert card(n) to a form that makes index splitting easier
  have h_card : Fintype.card n = (Fintype.card n - 1).succ := by
    rw [Nat.succ_eq_add_one]
    have hpos : 0 < Fintype.card n := by exact Fintype.card_pos
    have h := Nat.succ_pred_eq_of_pos hpos
    exact id (Eq.symm h)

  -- Break the sum over the polynomial into constant and non-constant parts
  have sum_decomp :
      ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i =
      p.coeff 0 • B^0 + ∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B^(i + 1) := by
    rw [h_card, Finset.sum_range_succ']
    exact Eq.symm (AddCommMagma.add_comm (p.coeff 0 • B ^ 0)
              (∑ i ∈ Finset.range ((Fintype.card n - 1).succ - 1), p.coeff (i + 1) • B ^ (i + 1)))

  -- Use the Cayley-Hamilton theorem and expand p(B) = 0 into matrix powers
  have h_expanded :
      B^(Fintype.card n) +
      (∑ i ∈ Finset.range (Fintype.card n), (Polynomial.coeff p i) • B ^ i) = 0 := by
    rw [← Polynomial.sum_monomial_eq p]
    rw [hpoly] at hCH
    simp only [Polynomial.aeval_def] at hCH
    have eq_sum : ∑ i ∈ Finset.range (Fintype.card n), (p.sum fun n a => Polynomial.monomial n a).coeff i • B ^ i = ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i := by
      rw [Polynomial.sum_monomial_eq]
    rw [eq_sum]
    rw [sum_decomp]
    simp only [pow_zero]
    rw [← add_zero (B^(Fintype.card n)), add_assoc, ← neg_eq_iff_add_eq_zero, zero_add]
    rw [neg_eq_iff_eq_neg, neg_add]
    rw [add_comm] at hCH

    -- Expand the evaluation of the sum of monomials at B
    have h_eval_expanded : Polynomial.eval₂ (algebraMap F (Matrix n n F)) B
        (∑ i ∈ Finset.range (1 + Fintype.card n), Polynomial.monomial i (p.coeff i)) =
        ∑ i ∈ Finset.range (1 + Fintype.card n),
          (Polynomial.eval₂ (algebraMap F (Matrix n n F)) B (Polynomial.monomial i (p.coeff i))) :=
      by rw [Polynomial.eval₂_finset_sum]

    -- Evaluate each monomial under eval₂
    have h_mono_eval : ∀ i, Polynomial.eval₂ (algebraMap F (Matrix n n F)) B
        (Polynomial.monomial i (p.coeff i)) = (algebraMap F (Matrix n n F)) (p.coeff i) * B^i := by
      intro i
      rw [Polynomial.eval₂_monomial]

    rw [h_eval_expanded] at hCH
    simp only [h_mono_eval, Algebra.algebraMap_eq_smul_one] at hCH

    -- Normalize index range
    have h_range_split : Finset.range (1 + Fintype.card n) = Finset.range (Fintype.card n + 1) := by rw [add_comm]
    rw [h_range_split, Finset.sum_range_succ] at hCH
    simp only [Algebra.smul_mul_assoc, one_mul] at hCH

    -- Leading coefficient of characteristic polynomial is 1
    have h_leading_coeff : p.coeff (Fintype.card n) = 1 := by
      rw [←Matrix.charpoly_natDegree_eq_dim B]
      exact (Matrix.charpoly_monic B).coeff_natDegree

    rw [h_leading_coeff, one_smul] at hCH

    -- Solve for B^n in terms of lower powers
    have h_eq : B^(Fintype.card n) = - ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i := by
      exact Eq.symm (neg_eq_of_add_eq_zero_right hCH)

    -- Break the sum again
    have h_split :
      ∑ i ∈ Finset.range (Fintype.card n), p.coeff i • B ^ i =
        p.coeff 0 • B^0 + ∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1) := by
      rw [sum_decomp]

    rw [h_split] at h_eq
    rw [pow_zero] at h_eq
    rw [h_eq]
    exact neg_add (p.coeff 0 • 1)
        (∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1))

  -- The constant term of the characteristic polynomial is nonzero if B is invertible
  have hc0 : Polynomial.coeff p 0 ≠ 0 := by
      have hdet := Matrix.det_eq_sign_charpoly_coeff B
      have hcoeff : (-1)^(Fintype.card n) * p.coeff 0 = B.det := by
        rw [hdet, mul_comm]
      have hdet_unit : IsUnit B.det := by
        exact (Matrix.isUnit_iff_isUnit_det B).mp hB
      by_contra h
      rw [h, mul_zero] at hcoeff
      have : B.det = 0 := by
        rw [← hcoeff]
      exact hdet_unit.ne_zero this

  -- Construct the identity matrix as a linear combination of B^i
  have hI :
    (1 : Matrix n n F) = (- (Polynomial.coeff p 0)⁻¹) • (Finset.sum (Finset.range (Fintype.card n)) fun i => Polynomial.coeff p (i + 1) • B ^ (i + 1)) := by
    have h_sum_zero := h_expanded
    have h_split := sum_decomp
    rw [h_split, pow_zero] at h_sum_zero
    have h_leading : p.coeff (Fintype.card n) = 1 := by
      rw [← Matrix.charpoly_natDegree_eq_dim B]
      exact (Matrix.charpoly_monic B).coeff_natDegree
    have h_Bpow : B ^ (Fintype.card n) = p.coeff (Fintype.card n) • B ^ (Fintype.card n) := by rw [h_leading, one_smul]
    rw [h_Bpow] at h_sum_zero

    -- Collect constant and nonconstant terms from Cayley-Hamilton
    have h_combined :
      0 = p.coeff 0 • 1 + (∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1)
          + p.coeff (Fintype.card n) • B ^ (Fintype.card n)) := by
      rw [← h_sum_zero]; abel

    have h_sumeq : ∑ i ∈ Finset.range (Fintype.card n), p.coeff (i + 1) • B ^ (i + 1) =
        ∑ i ∈ Finset.range (Fintype.card n - 1), p.coeff (i + 1) • B ^ (i + 1)
          + p.coeff (Fintype.card n) • B ^ (Fintype.card n) := by
      rw [h_card]
      exact Finset.sum_range_succ (fun x => p.coeff (x + 1) • B ^ (x + 1)) (Fintype.card n - 1)

    have h_combined' : 0 = p.coeff 0 • 1 + ∑ i ∈ Finset.range (Fintype.card n), p.coeff (i + 1) • B ^ (i + 1) :=
      by rw [← h_sumeq] at h_combined; exact h_combined

    -- Rearranged to express 1 as a linear combination
    have h_neg : (p.coeff 0 • 1) = - ∑ i ∈ Finset.range (Fintype.card n), p.coeff (i + 1) • B ^ (i + 1) :=
      eq_neg_of_add_eq_zero_left (eq_comm.mp h_combined')
    have h_neg' : - (p.coeff 0 • 1) = ∑ i ∈ Finset.range (Fintype.card n), p.coeff (i + 1) • B ^ (i + 1) := by rw [h_neg, neg_neg]
    rw [← h_neg', ← smul_neg]; simp
    rw [propext (eq_inv_smul_iff₀ hc0)]

  -- Multiply both sides by B⁻¹ to get B⁻¹ as linear combination of B^i
  have hBinv : B⁻¹ = (- (Polynomial.coeff p 0)⁻¹) • (Finset.sum (Finset.range (Fintype.card n)) (fun i => Polynomial.coeff p (i + 1) • (B ^ i))) := by
    rw [← Matrix.mul_one B⁻¹, hI, Matrix.mul_smul, Matrix.mul_sum]
    congr; ext i
    rw [Matrix.mul_smul, pow_succ']; simp

  -- Define coefficients b : Fin n → F such that B⁻¹ = ∑ b_i • B^i
  let b (i : Fin (Fintype.card n)) := - (Polynomial.coeff p (i.val + 1)) * (Polynomial.coeff p 0)⁻¹
  use b

  -- Match the shape with sum over Fin n
  rw [hBinv]; simp only [mul_comm _ (Polynomial.coeff p 0)⁻¹, b]
  ext i j
  simp only [Matrix.smul_apply, smul_eq_mul]

  -- Justify changing Finset.range n to Fin.univ mapped to indices
  let ι := Fin (Fintype.card n)
  have h_equiv : Finset.range (Fintype.card n) = Finset.univ.image (fun (x : ι) => x.val) := by
    apply Finset.ext
    intro x
    simp only [Finset.mem_range, Finset.mem_image, Finset.mem_univ, true_and]
    constructor
    · intro hx; exact CanLift.prf x hx
    · rintro ⟨y, _, rfl⟩; exact y.is_lt

  -- Rearrange sum over Fin to match expected format
  rw [h_equiv, Finset.sum_image]
  · have : ∀ i j, (-(p.coeff 0)⁻¹ * (∑ x : ι, p.coeff (↑x + 1) • B ^ (↑x : ℕ)) i j)
    = (∑ x : Fin (Fintype.card n), ((p.coeff 0)⁻¹ * -p.coeff (↑x + 1)) • B ^ (↑x : ℕ)) i j := by
      have h_matrix_eq :
        -(p.coeff 0)⁻¹ • ∑ x : ι, p.coeff (↑x + 1) • B ^ (↑x : ℕ)
        = ∑ x : ι, ((p.coeff 0)⁻¹ * -p.coeff (↑x + 1)) • B ^ (↑x : ℕ) := by
        rw [Finset.smul_sum]; simp; congr with x; rw [smul_smul]
      rw [<-h_matrix_eq]; exact fun i j => rfl
    rw [this]
  · intro x _ y _ h
    exact Fin.ext h