lean-eval Aristotle (Harmonic) submission for oppenheim_inequality

Challenge.lean

import Mathlib

open scoped MatrixOrder Matrix

theorem oppenheim_inequality {n : Type*} [Fintype n] [DecidableEq n]
    {A B : Matrix n n ℝ} (hA : A.PosSemidef) (hB : B.PosSemidef) :
    A.det * ∏ i, B i i ≤ (A ⊙ B).det := by
  sorry

config.json

{
  "challenge_module": "Challenge",
  "solution_module": "Solution",
  "theorem_names": [
    "oppenheim_inequality"
  ],
  "permitted_axioms": [
    "propext",
    "Quot.sound",
    "Classical.choice"
  ],
  "enable_nanoda": false
}

Helpers.lean

namespace Submission.Helpers

end Submission.Helpers

lakefile.toml

name = "oppenheim_inequality"
testDriver = "workspace_test"
defaultTargets = ["Challenge", "Solution", "Submission"]

[leanOptions]
autoImplicit = false

[[require]]
name = "mathlib"
git = "https://github.com/leanprover-community/mathlib4.git"
rev = "5450b53e5ddc"

[[lean_lib]]
name = "Challenge"

[[lean_lib]]
name = "Solution"

[[lean_lib]]
name = "Submission"

[[lean_exe]]
name = "workspace_test"
root = "WorkspaceTest"

lean-toolchain

leanprover/lean4:v4.30.0-rc2

Solution.lean

import Mathlib
import Submission

open scoped MatrixOrder Matrix

theorem oppenheim_inequality {n : Type*} [Fintype n] [DecidableEq n]
    {A B : Matrix n n ℝ} (hA : A.PosSemidef) (hB : B.PosSemidef) :
    A.det * ∏ i, B i i ≤ (A ⊙ B).det := by
  exact Submission.oppenheim_inequality hA hB

Submission.lean

import Mathlib

open scoped MatrixOrder Matrix
open Matrix Finset

namespace Submission

/-! ## Schur complement -/

/-- Schur complement of a matrix with respect to its (0,0) entry -/
noncomputable def schurComp {m : ℕ} (M : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ) :
    Matrix (Fin m) (Fin m) ℝ :=
  fun i j => M i.succ j.succ - M i.succ 0 * M 0 j.succ / M 0 0

/-! ## Helper lemmas -/

lemma hadamard_submatrix {α : Type*} [Mul α] {n m : Type*}
    (A B : Matrix n n α) (e : m → n) :
    (A ⊙ B).submatrix e e = (A.submatrix e e) ⊙ (B.submatrix e e) := by
  ext i j; simp [Matrix.hadamard_apply, Matrix.submatrix_apply]

lemma PosSemidef.entry_eq_zero_of_diag_eq_zero
    {k : ℕ} {M : Matrix (Fin k) (Fin k) ℝ} (hM : M.PosSemidef)
    {i : Fin k} (hi : M i i = 0) (j : Fin k) :
    M i j = 0 := by
  by_cases hij : i = j
  · rw [hij]; exact hi
  · exact hM.eq_zero_of_ne_diag_of_diag_eq_zero hi hij

lemma PosSemidef.submatrix_succ {m : ℕ}
    {M : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ} (hM : M.PosSemidef) :
    (M.submatrix Fin.succ Fin.succ).PosSemidef :=
  hM.submatrix _

lemma prod_fin_succ' {m : ℕ} (f : Fin (m + 1) → ℝ) :
    ∏ i, f i = f 0 * ∏ i : Fin m, f i.succ :=
  Fin.prod_univ_succ f

lemma hadamard_zero_zero {m : ℕ}
    (A B : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ) :
    (A ⊙ B) 0 0 = A 0 0 * B 0 0 :=
  Matrix.hadamard_apply A B 0 0

/-! ## Schur complement properties -/

/-
Determinant formula via Schur complement
-/
lemma det_eq_mul_det_schurComp {m : ℕ}
    (M : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ)
    (h : M 0 0 ≠ 0) :
    M.det = M 0 0 * (schurComp M).det := by
  -- Let's create the forward substitution matrix L.
  let L : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ := Matrix.of (fun i j => if i = 0 then if j = 0 then 1 else 0 else if j = 0 then -(M i 0 / M 0 0) else if i = j then 1 else 0);
  -- Define the matrix $N$ such that its entries are the same as $M$ except for the first column.
  let N : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ := L * M;
  -- Now let's compute the determinant of $N$.
  have h_det_N : Matrix.det N = M 0 0 * Matrix.det (Matrix.submatrix N (Fin.succ) (Fin.succ)) := by
    rw [ Matrix.det_succ_column_zero ];
    simp +zetaDelta at *;
    simp_all +decide [ Fin.sum_univ_succ, Matrix.mul_apply ];
  convert h_det_N using 1;
  · rw [ Matrix.det_mul, show Matrix.det L = 1 from ?_ ];
    · ring;
    · rw [ ← Matrix.det_transpose, Matrix.det_of_upperTriangular ];
      · aesop;
      · intro i j hij; aesop;
  · congr! 2;
    ext i j; simp +decide [ L, N, Matrix.mul_apply ] ; ring;
    simp +decide [ Finset.sum_ite, Finset.filter_eq', Finset.filter_ne', h, schurComp ];
    ring

/-
Schur complement of PSD matrix is PSD
-/
lemma schurComp_posSemidef {m : ℕ}
    {M : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ}
    (hM : M.PosSemidef) (h0 : 0 < M 0 0) :
    (schurComp M).PosSemidef := by
  refine' ⟨ _, fun x => _ ⟩;
  · ext i j; simp +decide [ *, schurComp ] ; ring;
    have := congr_fun ( congr_fun hM.1 i.succ ) j.succ; simp_all +decide [ Matrix.IsHermitian, mul_assoc, mul_comm, mul_left_comm ] ; ring;
    have := congr_fun ( congr_fun hM.1 0 ) i.succ; have := congr_fun ( congr_fun hM.1 0 ) j.succ; simp_all +decide [ Matrix.IsHermitian, mul_assoc, mul_comm, mul_left_comm ] ;
  · -- Let $y$ be the vector defined by $y_0 = -\frac{1}{M_{00}} \sum_{i=1}^m M_{0i} x_i$ and $y_i = x_i$ for $i > 0$.
    set y : Fin (m + 1) → ℝ := Fin.cons (- (∑ i, M 0 (Fin.succ i) * x i) / M 0 0) x;
    -- Then $y^T M y \geq 0$ (since $M$ is PSD).
    have hyM : 0 ≤ ∑ i, y i * ∑ j, M i j * y j := by
      have := hM.2;
      convert this ( Finsupp.equivFunOnFinite.symm y ) using 1 ; simp +decide [ Finsupp.sum_fintype, Finset.mul_sum _ _ _, mul_assoc, mul_comm, mul_left_comm ];
    simp +zetaDelta at *;
    simp_all +decide [ Fin.sum_univ_succ, Finset.sum_add_distrib, mul_add, add_mul, mul_assoc, mul_comm, mul_left_comm, div_eq_inv_mul, ne_of_gt ];
    convert hyM using 1 ; simp +decide [ Finset.mul_sum _ _ _, mul_assoc, mul_comm, mul_left_comm, Finset.sum_mul, schurComp ] ; ring!;
    simp +decide [ Finsupp.sum_fintype, mul_assoc, mul_comm, mul_left_comm, Finset.mul_sum _ _ _, Finset.sum_add_distrib ] ; ring!;

/-! ## Key identity for Oppenheim -/

/-
The key Schur complement identity for the Oppenheim inequality
-/
lemma oppenheim_schurComp_psd {m : ℕ}
    {A B : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ}
    (hA : A.PosSemidef) (hB : B.PosSemidef)
    (hA0 : 0 < A 0 0) (hB0 : 0 < B 0 0) :
    (schurComp (A ⊙ B) -
      (schurComp A) ⊙ (B.submatrix Fin.succ Fin.succ)).PosSemidef := by
  refine' ⟨ _, _ ⟩;
  · ext i j; simp +decide [ schurComp ] ; ring;
    have := hA.1; have := hB.1; simp_all +decide [ Matrix.IsHermitian, mul_assoc, mul_comm, mul_left_comm ] ;
    rw [ ← Matrix.ext_iff ] at *;
    simp_all +decide [ Matrix.transpose_apply, mul_assoc, mul_comm, mul_left_comm ] ; ring;
  · intro x
    have h_sum : ∑ i : Fin m, ∑ j : Fin m, x i * ((schurComp (A ⊙ B) - schurComp A ⊙ B.submatrix Fin.succ Fin.succ) i j) * x j = (1 / A 0 0) * ∑ i : Fin m, ∑ j : Fin m, (A (Fin.succ i) 0 * x i) * (schurComp B) i j * (A (Fin.succ j) 0 * x j) := by
      simp +decide [ schurComp, Finset.mul_sum _ _ _, Finset.sum_mul, mul_assoc, mul_comm, mul_left_comm, div_eq_inv_mul ];
      have h_symm : ∀ i j, A i j = A j i ∧ B i j = B j i := by
        have := hA.1; have := hB.1; simp_all +decide [ Matrix.IsHermitian, ← Matrix.ext_iff ] ;
      grind;
    -- Since $schurComp B$ is positive semidefinite, we have $\sum_{i,j} (A (Fin.succ i) 0 * x i) * (schurComp B) i j * (A (Fin.succ j) 0 * x j) \geq 0$.
    have h_schurComp_B_posSemidef : ∀ (y : Fin m → ℝ), 0 ≤ ∑ i : Fin m, ∑ j : Fin m, y i * (schurComp B) i j * y j := by
      have h_schurComp_B_posSemidef : (schurComp B).PosSemidef := by
        apply_rules [ schurComp_posSemidef ];
      have := h_schurComp_B_posSemidef.2;
      intro y; specialize this ( Finsupp.equivFunOnFinite.symm y ) ; simp_all +decide [ Finsupp.sum_fintype ] ;
    simp_all +decide [ Finsupp.sum_fintype ]

/-! ## Determinant monotonicity -/

/-
schurComp(1+Q) entry simplification
-/
lemma schurComp_one_add_entry {m : ℕ}
    (Q : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ)
    (h : (1 : ℝ) + Q 0 0 ≠ 0) (i j : Fin m) :
    schurComp (1 + Q) i j =
      (if i = j then 1 else 0) +
      (Q i.succ j.succ - Q i.succ 0 * Q 0 j.succ / (1 + Q 0 0)) := by
  -- By definition of Schur complement, we have:
  simp [schurComp];
  simp +decide [ Fin.ext_iff, Matrix.one_apply ] ; ring

/-
The correction matrix Q' - qq^T/(1+Q₀₀) is PSD
-/
lemma correction_psd {m : ℕ}
    {Q : Matrix (Fin (m + 1)) (Fin (m + 1)) ℝ}
    (hQ : Q.PosSemidef) :
    (Matrix.of fun i j : Fin m =>
      Q i.succ j.succ - Q i.succ 0 * Q 0 j.succ / (1 + Q 0 0)).PosSemidef := by
  simp +decide [ Matrix.PosSemidef ] at hQ ⊢;
  by_cases h : 1 + Q 0 0 = 0 <;> simp_all +decide [ Matrix.IsHermitian, Fin.sum_univ_succ ];
  · exact absurd ( hQ.2 ( Finsupp.single 0 1 ) ) ( by norm_num; nlinarith );
  · constructor;
    · rw [ ← Matrix.ext_iff ] at *;
      simp_all +decide [ Matrix.transpose_apply, mul_comm ];
    · intro x
      have h_sum : x.sum (fun i xi => x.sum (fun j xj => xi * Q i.succ j.succ * xj)) - (x.sum (fun i xi => xi * Q i.succ 0))^2 / (1 + Q 0 0) ≥ 0 := by
        have h_sum : ∀ (y : ℝ), (x.sum (fun i xi => x.sum (fun j xj => xi * Q i.succ j.succ * xj))) + 2 * y * (x.sum (fun i xi => xi * Q i.succ 0)) + y^2 * (1 + Q 0 0) ≥ 0 := by
          intro y
          have h_sum : (x.sum (fun i xi => x.sum (fun j xj => xi * Q i.succ j.succ * xj))) + 2 * y * (x.sum (fun i xi => xi * Q i.succ 0)) + y^2 * (1 + Q 0 0) = (Finsupp.equivFunOnFinite.symm (Fin.cons y x)).sum (fun i xi => (Finsupp.equivFunOnFinite.symm (Fin.cons y x)).sum (fun j xj => xi * Q i j * xj)) + y^2 := by
            simp +decide [ Finsupp.sum_fintype, Fin.sum_univ_succ ];
            simp +decide [ Finset.sum_add_distrib, mul_assoc, mul_comm, mul_left_comm, Finset.mul_sum _ _ _, Finset.sum_mul _ _ _, hQ.1 ] ; ring;
            simp +decide [ ← Matrix.ext_iff ] at *;
            simp +decide [ ← hQ.1, mul_two, Finset.sum_add_distrib ];
            ring;
          exact h_sum.symm ▸ add_nonneg ( hQ.2 _ ) ( sq_nonneg _ );
        contrapose! h_sum;
        use - (x.sum (fun i xi => xi * Q i.succ 0)) / (1 + Q 0 0);
        convert h_sum using 1 ; ring;
        grind +splitImp;
      convert h_sum.le using 1 ; simp +decide [ Finset.sum_add_distrib, mul_sub, sub_mul, mul_assoc, mul_comm, mul_left_comm, div_eq_mul_inv, h ] ; ring;
      simp +decide [ sq, mul_assoc, mul_comm, mul_left_comm, Finset.mul_sum _ _ _, Finset.sum_mul, Finsupp.sum ];
      exact Finset.sum_congr rfl fun i hi => Finset.sum_congr rfl fun j hj => by rw [ ← Matrix.ext_iff ] at hQ; aesop;

/-- det(I + Q) ≥ 1 for PSD Q -/
lemma det_one_add_posSemidef_ge_one :
    ∀ (m : ℕ) (Q : Matrix (Fin m) (Fin m) ℝ),
      Q.PosSemidef → 1 ≤ (1 + Q).det := by
  intro m
  induction m with
  | zero => intro Q _; simp [Matrix.det_isEmpty]
  | succ m ih =>
    intro Q hQ
    have hQ00 : 0 ≤ Q 0 0 := hQ.diagonal_nonneg 0
    have h1Q00 : (1 : ℝ) + Q 0 0 ≠ 0 := by linarith
    have h1Q00_pos : 1 ≤ 1 + Q 0 0 := by linarith
    -- det(1+Q) = (1 + Q₀₀) * det(schurComp(1+Q))
    have hdet := det_eq_mul_det_schurComp (1 + Q) (by simp; linarith)
    rw [show (1 + Q) 0 0 = 1 + Q 0 0 from by simp] at hdet
    -- schurComp(1+Q) = 1 + R with R PSD
    have hR_psd := correction_psd hQ
    have hSC_eq : schurComp (1 + Q) = 1 + Matrix.of (fun i j : Fin m =>
        Q i.succ j.succ - Q i.succ 0 * Q 0 j.succ / (1 + Q 0 0)) := by
      ext i j
      rw [schurComp_one_add_entry Q h1Q00]
      simp [Matrix.add_apply, Matrix.one_apply, Matrix.of_apply]
    -- By induction: det(1 + R) ≥ 1
    have hR_det := ih _ hR_psd
    rw [hSC_eq] at hdet
    -- Combine
    linarith [mul_le_mul_of_nonneg_left hR_det (by linarith : (0 : ℝ) ≤ 1 + Q 0 0)]

/-
det(M + P) ≥ det(M) for PSD M, P
-/
lemma det_add_posSemidef_le {n : Type*} [Fintype n] [DecidableEq n]
    {M P : Matrix n n ℝ} (hM : M.PosSemidef) (hP : P.PosSemidef) :
    M.det ≤ (M + P).det := by
  by_cases h_det_M_zero : M.det = 0;
  · exact h_det_M_zero.symm ▸ Matrix.PosSemidef.det_nonneg ( hM.add hP );
  · -- Since $M$ is positive definite, there exists an invertible matrix $B$ such that $M = B^T B$.
    obtain ⟨B, hB⟩ : ∃ B : Matrix n n ℝ, M = B.transpose * B ∧ IsUnit B.det := by
      have := Matrix.posSemidef_iff_eq_conjTranspose_mul_self.mp hM;
      obtain ⟨ B, rfl ⟩ := this; use B; simp_all +decide [ Matrix.det_mul ] ;
    have h_det_one_add_posSemidef_ge_one : 1 ≤ (1 + (B⁻¹).transpose * P * B⁻¹).det := by
      have h_det_one_add_posSemidef_ge_one : (B⁻¹).transpose * P * B⁻¹ |> Matrix.PosSemidef := by
        convert hP.conjTranspose_mul_mul_same ( B⁻¹ ) using 1;
      convert det_one_add_posSemidef_ge_one ( Fintype.card n ) ( ( B⁻¹ᵀ * P * B⁻¹ ).submatrix ( fun i => Fintype.equivFin n |>.symm i ) ( fun i => Fintype.equivFin n |>.symm i ) ) _ using 1;
      · convert Matrix.det_reindex_self ( Fintype.equivFin n |> Equiv.symm ) _ using 2;
        ext i j; simp +decide [ Matrix.one_apply ] ;
      · convert h_det_one_add_posSemidef_ge_one.submatrix _ using 1;
    have h_det_M_plus_P : (M + P).det = B.det ^ 2 * (1 + (B⁻¹).transpose * P * B⁻¹).det := by
      have h_det_M_plus_P : (M + P).det = (B.transpose * (1 + (B⁻¹).transpose * P * B⁻¹) * B).det := by
        simp +decide [ hB.1, mul_add, add_mul, mul_assoc, hB.2 ];
        simp +decide [ Matrix.transpose_nonsing_inv, hB.2 ];
        cases hB.2.nonempty_invertible ; aesop;
      simp_all +decide [ sq, mul_assoc, mul_comm, mul_left_comm ];
    simp_all +decide [ Matrix.det_mul ];
    nlinarith

/-! ## Oppenheim inequality for Fin m -/

lemma oppenheim_fin :
    ∀ (m : ℕ) (A B : Matrix (Fin m) (Fin m) ℝ),
      A.PosSemidef → B.PosSemidef →
      A.det * ∏ i, B i i ≤ (A ⊙ B).det := by
  intro m
  induction m with
  | zero =>
    intro A B _ _
    simp [Matrix.det_isEmpty, Finset.univ_eq_empty]
  | succ m ih =>
    intro A B hA hB
    by_cases hA0 : A 0 0 = 0
    · -- A₀₀ = 0 ⟹ det(A) = 0
      have hrow : ∀ j, A 0 j = 0 := fun j =>
        PosSemidef.entry_eq_zero_of_diag_eq_zero hA hA0 j
      simp [Matrix.det_eq_zero_of_row_eq_zero 0 hrow]
      exact (hA.hadamard hB).det_nonneg
    · by_cases hB0 : B 0 0 = 0
      · -- B₀₀ = 0 ⟹ ∏ B_{ii} = 0
        simp [Finset.prod_eq_zero (Finset.mem_univ (0 : Fin (m+1))) hB0]
        exact (hA.hadamard hB).det_nonneg
      · -- Main case: A₀₀ > 0, B₀₀ > 0
        have hA0' : 0 < A 0 0 :=
          lt_of_le_of_ne (hA.diagonal_nonneg 0) (Ne.symm hA0)
        have hB0' : 0 < B 0 0 :=
          lt_of_le_of_ne (hB.diagonal_nonneg 0) (Ne.symm hB0)
        have hAB00 : (A ⊙ B) 0 0 = A 0 0 * B 0 0 := hadamard_zero_zero A B
        have hAB00_pos : 0 < (A ⊙ B) 0 0 := by rw [hAB00]; positivity
        -- Schur complement facts
        have hSA := schurComp_posSemidef hA hA0'
        have hB' := hB.submatrix_succ
        have hkey := oppenheim_schurComp_psd hA hB hA0' hB0'
        -- Inductive hypothesis on Schur complement
        have hind := ih (schurComp A) (B.submatrix Fin.succ Fin.succ) hSA hB'
        -- Det monotonicity
        have hdet_mono := det_add_posSemidef_le (hSA.hadamard hB') hkey
        rw [sub_add_cancel] at hdet_mono
        -- Det formulas
        have hdetA := det_eq_mul_det_schurComp A hA0
        have hdetAB := det_eq_mul_det_schurComp (A ⊙ B) (ne_of_gt hAB00_pos)
        have hprod := prod_fin_succ' (fun i => B i i)
        -- Combine everything
        rw [hdetA, hdetAB, hAB00, hprod]
        have h1 : (schurComp (A ⊙ B)).det ≥ (schurComp A).det * ∏ i : Fin m, B i.succ i.succ :=
          le_trans hind hdet_mono
        nlinarith [hA0'.le, hB0'.le, hSA.det_nonneg,
          Finset.prod_nonneg (fun i _ => hB.diagonal_nonneg i.succ)]

/-! ## Main theorem -/

theorem oppenheim_inequality {n : Type*} [Fintype n] [DecidableEq n]
    {A B : Matrix n n ℝ} (hA : A.PosSemidef) (hB : B.PosSemidef) :
    A.det * ∏ i, B i i ≤ (A ⊙ B).det := by
  let e := Fintype.equivFin n
  have hAe : (A.submatrix e.symm e.symm).PosSemidef := by
    rwa [Matrix.PosSemidef.submatrix_equiv]
  have hBe : (B.submatrix e.symm e.symm).PosSemidef := by
    rwa [Matrix.PosSemidef.submatrix_equiv]
  rw [← Matrix.det_submatrix_equiv_self e,
      ← Matrix.det_submatrix_equiv_self e,
      ← hadamard_submatrix,
      ← Finset.prod_equiv e (by simp) (by simp)]
  exact oppenheim_fin _ _ _ hAe hBe

end Submission

WorkspaceTest.lean

import Lean

open Lean

def comparatorExists (comparatorBin : String) : IO Bool := do
  if comparatorBin.contains '/' then
    return (← System.FilePath.pathExists comparatorBin)
  try
    let child ← IO.Process.spawn {
      cmd := "sh"
      args := #["-c", "command -v \"$1\" >/dev/null 2>&1", "sh", comparatorBin]
    }
    let exitCode ← child.wait
    return exitCode == 0
  catch _ =>
    return false

def main : IO UInt32 := do
  let comparatorBin := (← IO.getEnv "COMPARATOR_BIN").getD "comparator"
  if !(← comparatorExists comparatorBin) then
    IO.eprintln s!"Failed to run comparator via `{comparatorBin}`."
    IO.eprintln "Make sure `comparator` is installed and on your `PATH`, or set `COMPARATOR_BIN=/path/to/comparator`."
    IO.eprintln "See the root repository README for comparator setup details, including landrun and lean4export."
    pure 1
  else
    try
      let child ← IO.Process.spawn {
        cmd := "lake"
        args := #["env", comparatorBin, "config.json"]
      }
      let exitCode ← child.wait
      pure exitCode
    catch err =>
      IO.eprintln s!"Failed to run comparator via `{comparatorBin}`."
      IO.eprintln "Make sure `comparator` is installed and on your `PATH`, or set `COMPARATOR_BIN=/path/to/comparator`."
      IO.eprintln "See the root repository README for comparator setup details, including landrun and lean4export."
      IO.eprintln s!"Original error: {err}"
      pure 1