Exclude zero terms from the sum

Stream: new members

[jsodd]

I've been trying to solve the problem suggested recently here:
Eric Wieser said:

Here's the statement, in case anyone feels like being sniped:

import Mathlib

variable {m n R} [CommRing R] [StarRing R] [PartialOrder R] [StarOrderedRing R]


/-- https://en.wikipedia.org/wiki/Sylvester%27s_criterion -/
theorem posDef_fin_iff_det {A : Matrix (Fin n) (Fin n) R} :
    A.PosDef ↔ ∀ i (h : i ≤ n), 0 < (A.submatrix (Fin.castLE h) (Fin.castLE h)).det :=
  -- this is the interesting mathematical result
  sorry

It turned out that even the "easy" direction of Sylvester's criterion is hard for me to formalise. I decided to start with a preliminary lemma showing that A.PosDef implies that the minors of A are posdef

lemma posDef_minor {A : Matrix (Fin n) (Fin n) R} :
  A.PosDef → ∀ i (h : i ≤ n), (A.submatrix (Fin.castLE h) (Fin.castLE h)).PosDef := by
    intro h₀ i h₁
    constructor
    · exact Matrix.IsHermitian.submatrix (h₀.1) (Fin.castLE h₁)
    · intro x hx
      let y : Fin n → R := λ j ↦ if h : j < i then x ⟨j, h⟩ else 0
      -- prove that y A y = x A' x and use y A y > 0
      sorry

I have no idea how to approach this. I've tried using Mathlib lemmas, such as this one

theorem submatrix_mulVec_equiv [Fintype n] [Fintype o] [NonUnitalNonAssocSemiring α]
    (M : Matrix m n α) (v : o → α) (e₁ : l → m) (e₂ : o ≃ n) :
    M.submatrix e₁ e₂ *ᵥ v = (M *ᵥ (v ∘ e₂.symm)) ∘ e₁ :=
  funext fun _ => Eq.symm (dotProduct_comp_equiv_symm _ _ _)

but unfortunately the condition (e₂ : o ≃ n) is too restrictive to handle the "extension by zero" idea. I've also tried expanding the sum, but then I don't know how to throw away all the zeroes from it...

[jsodd]

Here's the code up to a point which I cannot prove:

lemma posDef_minor {A : Matrix (Fin n) (Fin n) R} :
  A.PosDef → ∀ i (h : i ≤ n), (A.submatrix (Fin.castLE h) (Fin.castLE h)).PosDef := by
    intro h₀ i h₁
    constructor
    · exact Matrix.IsHermitian.submatrix (h₀.1) (Fin.castLE h₁)
    · intro x hx
      let y : Fin n → R := λ j ↦ if h : j < i then x ⟨j, h⟩ else 0
      have h₂ : y ≠ 0 := by
        intro hy
        apply hx
        ext j
        have : y (Fin.castLE h₁ j) = x j :=  dif_pos j.isLt
        rw [← this, hy]
        simp only [Pi.zero_apply]
      calc
        0 < Matrix.dotProduct (star y) (A.mulVec y) := h₀.2 y h₂
        _ = Matrix.dotProduct (star x)
          ((A.submatrix (Fin.castLE h₁) (Fin.castLE h₁)).mulVec x) := by
          -- ??
          sorry

[Eric Wieser]

The key lemma is docs#Finset.sum_filter_add_sum_filter_not, I think. This works, but is a mess:

import Mathlib

variable {m n R} [Fintype m] [Fintype n] [CommRing R] [PartialOrder R] [StarRing R]

lemma Matrix.PosDef.submatrix {A : Matrix m m R} (hA : A.PosDef) (f : n ↪ m):
    (A.submatrix f f).PosDef := by
  constructor
  · exact hA.isHermitian.submatrix f
  · intro x hx
    let y : m → R := Function.extend f x 0
    have hyf (i) : y (f i) = x i := by simp [y, (f.injective.factorsThrough _).extend_apply]
    have h₂ : y ≠ 0 := by
      intro hy
      apply hx
      ext j
      replace h := congr_fun hy (f j)
      simpa [y, (f.injective.factorsThrough _).extend_apply] using h
    calc
      0 < Matrix.dotProduct (star y) (A.mulVec y) := hA.2 y h₂
      _ = Matrix.dotProduct (star x) ((A.submatrix f f).mulVec x) := by
        -- ??
        simp only [Matrix.dotProduct, Matrix.mulVec]
        classical
        simp only [← Finset.sum_filter_add_sum_filter_not Finset.univ (· ∈ Set.range f), mul_add, Finset.mul_sum,
          Finset.sum_add_distrib]
        simp only [Set.mem_range, Finset.univ_filter_exists, Pi.star_apply, Finset.mem_univ,
          EmbeddingLike.apply_eq_iff_eq, imp_self, implies_true, Finset.sum_image, hyf,
          submatrix_apply]
        simp (config := {contextual := true}) [Finset.sum_filter, y, -not_exists]

[jsodd]

Wow, thanks a lot! How do you even find all these lemmas! I looked through Loogle and Moogle, but couldn't come up with it.

[jsodd]

So, Finset.sum_filter_add_sum_filter_not is what I was looking for. If I understand correctly, it allows to split a sum in two using some decidable proposition-valued function?

[Eric Wieser]

Yep, and classical makes any such function decidable

[Eric Wieser]

@jsodd, are you interested in PRing the above to mathlib?

[Eric Wieser]

No need to get the full result you were aiming for straight away, that lemma by itself would be a good one to have

[jsodd]

Eric Wieser said:

jsodd, are you interested in PRing the above to mathlib?

Thanks, I'll think about it! In the meantime, I've managed to prove one more lemma for Sylvester's criterion. Here it is:

lemma Matrix.PosDef.submatrix.det_pos {A : Matrix n n 𝕜}
  (hA : A.PosDef) (e : m ↪ n) :
    0 < (A.submatrix e e).det :=
      Matrix.PosDef.det_pos (Matrix.PosDef.submatrix hA e)

def Matrix.toBilin_block₁₁
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :=
    Matrix.dotProduct (x ∘ Sum.inl) (A.toBlocks₁₁.mulVec (y ∘ Sum.inl))

def Matrix.toBilin_block₁₂
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :=
    Matrix.dotProduct (x ∘ Sum.inl) (A.toBlocks₁₂.mulVec (y ∘ Sum.inr))

def Matrix.toBilin_block₂₁
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :=
    Matrix.dotProduct (x ∘ Sum.inr) (A.toBlocks₂₁.mulVec (y ∘ Sum.inl))

def Matrix.toBilin_block₂₂
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :=
    Matrix.dotProduct (x ∘ Sum.inr) (A.toBlocks₂₂.mulVec (y ∘ Sum.inr))

def Matrix.toBilin_block
  (A : Matrix (n ⊕ o) (l ⊕ m) R) :=
    A.toBilin_block₁₁ + A.toBilin_block₁₂ + A.toBilin_block₂₁ + A.toBilin_block₂₂

def Matrix.toBilin_block'
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :=
    A.toBilin_block x y

lemma Matrix.toBilin_block_apply
  (A : Matrix (n ⊕ o) (l ⊕ m) R) (x : (n ⊕ o) → R) (y : (l ⊕ m) → R) :
    Matrix.dotProduct x (A.mulVec y) = A.toBilin_block x y := by
      dsimp only [dotProduct, mulVec, toBilin_block, Pi.add_apply, toBilin_block₁₁,
        Function.comp_apply, toBlocks₁₁, of_apply, toBilin_block₁₂, toBlocks₁₂, toBilin_block₂₁,
        toBlocks₂₁, toBilin_block₂₂, toBlocks₂₂]
      simp only [Fintype.sum_sum_type]
      simp only [mul_add, Finset.sum_add_distrib, Finset.mul_sum]
      ring

It may also be useful...

[jsodd]

Here's my attempt at formalising Sylvester's criterion. A bit messy and I still struggle to close some of the goals, but here it is:

import Mathlib

universe u
variable {m n l o R 𝕜 : Type u}
variable [Fintype m] [Fintype n] [DecidableEq m] [DecidableEq n] [Nonempty n]
variable [Fintype l] [Fintype o] [DecidableEq l] [DecidableEq o]
variable [CommRing R] [PartialOrder R] [StarRing R]
variable [RCLike 𝕜]

open Classical ComplexOrder

lemma Matrix.PosDef.submatrix {A : Matrix n n R} (hA : A.PosDef) (f : m ↪ n):
    (A.submatrix f f).PosDef := by
  constructor
  · exact hA.isHermitian.submatrix f
  · intro x hx
    let y : n → R := Function.extend f x 0
    have hyf : ∀ i, y (f i) = x i := by simp [y, (f.injective.factorsThrough _).extend_apply]
    have h₂ : y ≠ 0 := by
      intro hy
      apply hx
      ext j
      replace h := congr_fun hy (f j)
      simpa [y, (f.injective.factorsThrough _).extend_apply] using h
    calc
      0 < Matrix.dotProduct (star y) (A.mulVec y) := hA.2 y h₂
      _ = Matrix.dotProduct (star x) ((A.submatrix f f).mulVec x) := by
        simp only [Matrix.dotProduct, Matrix.mulVec]
        simp only [← Finset.sum_filter_add_sum_filter_not Finset.univ (· ∈ Set.range f), mul_add, Finset.mul_sum,
          Finset.sum_add_distrib]
        simp only [Set.mem_range, Finset.univ_filter_exists, Pi.star_apply, Finset.mem_univ,
          EmbeddingLike.apply_eq_iff_eq, imp_self, implies_true, Finset.sum_image, hyf,
          submatrix_apply]
        simp (config := {contextual := true}) [Finset.sum_filter, y, -not_exists]


lemma Matrix.PosDef.submatrix.det_pos {A : Matrix n n 𝕜}
  (hA : A.PosDef) (e : m ↪ n) :
    0 < (A.submatrix e e).det :=
      Matrix.PosDef.det_pos (Matrix.PosDef.submatrix hA e)

variable [MulPosReflectLT 𝕜] [PosMulReflectLT 𝕜]

-- /-- https://en.wikipedia.org/wiki/Sylvester%27s_criterion -/
theorem posDef_fin_iff_det {A : Matrix n n 𝕜} (hA : A.IsHermitian) :
    A.PosDef ↔ ∀ (m : Type u) [Fintype m] [DecidableEq m] (e : m ↪ n),
      0 < (A.submatrix e e).det := by
      constructor
      · intro hA m hm e
        exact Matrix.PosDef.submatrix.det_pos hA
      · intro h_sub
        constructor
        exact hA
        generalize hN : Fintype.card n = N
        cases N with
          | zero =>
            exact absurd hN Fintype.card_ne_zero
          | succ N =>
            induction N generalizing n with
            | zero =>
              rw [zero_add] at hN
              have h_det : 0 < A.det := by
                apply h_sub n ⟨id, Function.injective_id⟩
              intro x hx
              dsimp only [Matrix.dotProduct, Matrix.mulVec]
              obtain ⟨k, hk⟩ := Fintype.card_eq_one_iff.mp hN
              rw [Matrix.det_eq_elem_of_card_eq_one hN k] at h_det
              simp only [hN, hk, Finset.sum_const, Finset.card_univ, nsmul_eq_mul]
              simp only [Nat.cast_one, one_mul, mul_comm, mul_assoc]
              apply (mul_pos_iff_of_pos_left h_det).2
              apply mul_star_self_pos
              apply isRegular_of_ne_zero
              intro hk_neg
              apply hx
              ext j
              simp [hk]
              exact hk_neg
            | succ N ih =>
              intro x hx
              rename_i _ _ _ hn _ _
              let k := Classical.arbitrary n
              let s := hn.1 \ {k}
              have hs : Disjoint s {k} := Finset.sdiff_disjoint
              -- hs₁ follows from card n ≥ 2
              have hs₁ : Nonempty s := sorry
              let e := fun (a : s) ↦ (a : n)
              have he : Function.Injective e := Subtype.val_injective
              -- hN₁ follows from s = n \ {k}, where N = card n
              have hN₁ : Fintype.card s = N + 1 := sorry
              have hn_u : Finset.univ = s ∪ {k} := by
                refine Eq.symm (Finset.sdiff_union_of_subset ?h)
                simp only [Finset.singleton_subset_iff]
                exact Fintype.complete k
              let v : s → 𝕜 := fun j ↦ A j k
              let x' := x ∘ e
              let A' := A.submatrix e e
              let y := (x k) • (A'⁻¹.mulVec v) + x'
              let C := (star x k) * x k * (A k k - Matrix.dotProduct (star v) (A'⁻¹.mulVec v))
              -- Proof by computation
              have H : Matrix.dotProduct (star x) (A.mulVec x)
                = Matrix.dotProduct (star y) (A'.mulVec y) + C := sorry
              rw [H]
              have H₁ : 0 ≤ C := by
                dsimp only [C]
                have H₁₁ : 0 < (A k k - Matrix.dotProduct (star v) (A'⁻¹.mulVec v)) := sorry
                have H₁₂ : 0 ≤ star x k * x k := sorry
                sorry
              have H₂ : 0 < Matrix.dotProduct (star y) (A'.mulVec y) := sorry

              exact Right.add_pos_of_pos_of_nonneg H₂ H₁

[Eric Wieser]

Do you need to split into 0, 1, and n + 2, or would 0 and n + 1 work as a base case?

[Eric Wieser]

Note that this is not quite Sylvester's criterion; you've shown it for all possible principal minors, not the stronger case of just the _leading_ principal ones

[jsodd]

Eric Wieser said:

Do you need to split into 0, 1, and n + 2, or would 0 and n + 1 work as a base case?

Well, if n = 0 it's just not true...

[jsodd]

Eric Wieser said:

Note that this is not quite Sylvester's criterion; you've shown it for all possible principal minors, not the stronger case of just the _leading_ principal ones

Yes, I know. I am still shocked by how tedious and hard it turns out to formalize something that obvious as the Sylvester's criterion. I hope it gets better with experience...

[Eric Wieser]

@jsodd, do you want to PR your Matrix.PosDef.submatrix above to mathlib?