CubicPathwidth.T2.FcubicBounds

Source file CubicPathwidth/T2/FcubicBounds.lean of the Lean project (how it is built).

import CubicPathwidth.T2.Bridge
import CubicPathwidth.C3.MinBound
import CubicPathwidth.C3.IdealStar
import CubicPathwidth.C3.Constants

/-!
# Two-sided bounds on the cubic plain-sweep constant `F`

`1/50 ≤ Fcubic ≤ 103/400`.  The ideal 3-star is realized by explicit vectors in `ℝ³`;
by the bridge its frontier probability is a standard Gaussian probability.  Lower bound:
`Fcubic ≥ P(Y₀ > 0 ≥ Y₁) = arccos(r)/(2π)`.  Upper bound: a union bound over the three
leaves, and for a pair of unit vectors with inner product `r` the sum/difference
decomposition `Y₀ = αM + βD`, `Y₁ = αM − βD` with independent standard normals `M, D`.
-/

noncomputable section
open scoped RealInnerProductSpace
open MeasureTheory ProbabilityTheory Filter Topology Set

namespace CubicPathwidth
namespace T2

open GaussianField (Sample)

/-! ## The ideal 3-star as explicit vectors -/

/-- The ideal 3-star vectors: center `e₀`, leaves `r e₀ + √(1-r²) u_j`. -/
def star3Vec (r : ℝ) : Fin 4 → Sample 3 :=
  Fin.cons (C3.Ideal.e 0) (C3.Ideal.bstar r)

theorem r33_pos : 0 < r33 := by
  unfold r33; positivity

theorem r33_sq : r33 ^ 2 = 8 / 9 := by
  unfold r33
  have : Real.sqrt 2 ^ 2 = 2 := Real.sq_sqrt (by norm_num)
  nlinarith

theorem r33_lt_one : r33 < 1 := by
  have := r33_sq
  nlinarith [r33_pos]

theorem gram_star3Vec : gram (star3Vec r33) = idealStarCov 3 r33 := by
  have h1 : r33 ≤ 1 := r33_lt_one.le
  have h2 : -1 ≤ r33 := by linarith [r33_pos]
  have hee : ∀ i j : Fin 3, i ≠ j → inner ℝ (C3.Ideal.e i) (C3.Ideal.e j) = 0 := by
    intro i j h; simpa [h] using C3.Ideal.inner_e i j
  have h00 : inner ℝ (C3.Ideal.e 0) (C3.Ideal.e 0) = 1 := by simpa using C3.Ideal.inner_e 0 0
  have hsq := r33_sq
  ext i j
  simp only [gram, idealStarCov, star3Vec]
  rcases Fin.eq_zero_or_eq_succ i with rfl | ⟨a, rfl⟩ <;>
    rcases Fin.eq_zero_or_eq_succ j with rfl | ⟨b, rfl⟩
  · simp only [Fin.cons_zero, h00, ite_true]
  · simp only [Fin.cons_zero, Fin.cons_succ, C3.Ideal.inner_e0_bstar]
    simp [(Fin.succ_ne_zero b).symm]
  · simp only [Fin.cons_zero, Fin.cons_succ]
    rw [real_inner_comm, C3.Ideal.inner_e0_bstar]
    simp [Fin.succ_ne_zero]
  · have := C3.Ideal.inner_bstar h1 h2 a b
    simp only [Fin.cons_succ]
    rw [this]
    by_cases hab : a = b
    · simp [hab]
    · have h3 : (a.succ : Fin 4) ≠ b.succ := fun h => hab (Fin.succ_injective _ h)
      simp [hab, h3, C3.pairCorr, Fin.succ_ne_zero]
      nlinarith

/-! ## A Gaussian integral and the sum/difference pair bound -/

theorem integral_posPart_gaussian :
    ∫ d, max d 0 ∂gaussianReal 0 1 = (Real.sqrt (2 * Real.pi))⁻¹ := by
  rw [integral_gaussianReal_eq_integral_smul one_ne_zero]
  set c : ℝ := (Real.sqrt (2 * Real.pi))⁻¹ with hc
  have hform : ∀ x : ℝ, gaussianPDFReal 0 1 x = c * Real.exp (-(x ^ 2 / 2)) := by
    intro x
    rw [gaussianPDFReal_def]
    simp only [NNReal.coe_one, mul_one, sub_zero, hc]
    congr 2
    ring
  have hfun : (fun x : ℝ => gaussianPDFReal 0 1 x • max x 0) =
      Set.indicator (Ioi 0) (fun x => c * (x * Real.exp (-(x ^ 2 / 2)))) := by
    funext x
    rw [Set.indicator_apply, hform]
    by_cases hx : x ∈ Ioi (0:ℝ)
    · simp only [hx, ite_true, smul_eq_mul]
      have : max x 0 = x := max_eq_left (le_of_lt hx)
      rw [this]; ring
    · simp only [hx, ite_false, smul_eq_mul]
      have : max x 0 = 0 := max_eq_right (not_lt.mp hx)
      rw [this]; ring
  rw [hfun, integral_indicator measurableSet_Ioi, integral_const_mul]
  have hder : ∀ x ∈ Ici (0:ℝ), HasDerivAt (fun x => -Real.exp (-(x ^ 2 / 2)))
      (x * Real.exp (-(x ^ 2 / 2))) x := by
    intro x _
    have h1 : HasDerivAt (fun x : ℝ => -(x ^ 2 / 2)) (-x) x := by
      have := ((hasDerivAt_pow 2 x).div_const 2).neg
      convert this using 1
      simp
    have := (h1.exp).neg
    convert this using 1
    ring
  have hint : IntegrableOn (fun x : ℝ => x * Real.exp (-(x ^ 2 / 2))) (Ioi 0) := by
    have := (integrable_mul_exp_neg_mul_sq (b := 1 / 2) (by norm_num)).integrableOn (s := Ioi 0)
    refine this.congr_fun (fun x _ => ?_) measurableSet_Ioi
    simp only
    congr 2
    ring
  have hlim : Tendsto (fun x : ℝ => -Real.exp (-(x ^ 2 / 2))) atTop (𝓝 0) := by
    have h : Tendsto (fun x : ℝ => -(x ^ 2 / 2)) atTop atBot := by
      apply tendsto_neg_atTop_atBot.comp
      exact (tendsto_pow_atTop (by norm_num)).atTop_div_const (by norm_num)
    have := (Real.tendsto_exp_atBot.comp h).neg
    simpa using this
  rw [integral_Ioi_of_hasDerivAt_of_tendsto' hder hint hlim]
  simp [hc]

theorem integrable_posPart_gaussian : Integrable (fun d : ℝ => max d 0) (gaussianReal 0 1) := by
  have h : Integrable (fun d : ℝ => d) (gaussianReal 0 1) :=
    (memLp_id_gaussianReal 1).integrable (by simp)
  refine (h.norm).mono' (by fun_prop) (ae_of_all _ fun d => ?_)
  rw [Real.norm_eq_abs, abs_of_nonneg (le_max_right _ _)]
  simp only [Real.norm_eq_abs]
  exact max_le (le_abs_self d) (abs_nonneg d)

variable {E : Type*} [NormedAddCommGroup E] [InnerProductSpace ℝ E]
  [FiniteDimensional ℝ E] [MeasurableSpace E] [BorelSpace E]

/-- For orthonormal `a, b`, the vectors `p = αa + βb` and `q = αa - βb` give
`P(τ < ⟨p,x⟩, ⟨q,x⟩ ≤ τ) ≤ (β/α)/√(2π)`: with the independent standard normals
`D = ⟨b,x⟩`, `M = ⟨a,x⟩` the event forces `D > 0` and `M` in an interval of length `2βD/α`. -/
theorem pair_core {a b : E} (ha : ⟪a, a⟫ = 1) (hb : ⟪b, b⟫ = 1) (hab : ⟪a, b⟫ = 0)
    {α β : ℝ} (hα : 0 < α) (hβ : 0 < β) (τ : ℝ) :
    (stdGaussian E).real
        {x | τ < ⟪α • a + β • b, x⟫ ∧ ⟪α • a - β • b, x⟫ ≤ τ} ≤
      β / α * (Real.sqrt (2 * Real.pi))⁻¹ := by
  set μ := stdGaussian E
  have hba : ⟪b, a⟫ = 0 := by rw [real_inner_comm]; exact hab
  set T : Set (ℝ × ℝ) := {q | τ < α * q.2 + β * q.1 ∧ α * q.2 - β * q.1 ≤ τ} with hT
  have hTm : MeasurableSet T := by
    refine MeasurableSet.inter ?_ ?_
    · exact measurableSet_lt (f := fun _ : ℝ × ℝ => τ) (g := fun q : ℝ × ℝ => α * q.2 + β * q.1)
        measurable_const (by fun_prop)
    · exact measurableSet_le (f := fun q : ℝ × ℝ => α * q.2 - β * q.1) (g := fun _ : ℝ × ℝ => τ)
        (by fun_prop) measurable_const
  have hset : {x : E | τ < ⟪α • a + β • b, x⟫ ∧ ⟪α • a - β • b, x⟫ ≤ τ} =
      (fun x => (⟪b, x⟫, ⟪a, x⟫)) ⁻¹' T := by
    ext x
    simp only [mem_ofPred_eq, mem_preimage, hT, inner_add_left, inner_sub_left,
      real_inner_smul_left]
  have hmap := GaussianAngle.orthogonal_pair_map hb ha hba
  rw [hset, ← map_measureReal_apply (by fun_prop) hTm, hmap]
  change ((gaussianReal 0 1).prod (gaussianReal 0 1)).real T ≤ _
  rw [C3.real_prod_eq_integral _ _ hTm]
  have hle : ∀ d : ℝ, (gaussianReal 0 1).real (Prod.mk d ⁻¹' T) ≤ β / α * max d 0 := by
    intro d
    by_cases hd : 0 < d
    · have hsub : Prod.mk d ⁻¹' T ⊆ Ioc ((τ - β * d) / α) ((τ + β * d) / α) := by
        intro m hm
        simp only [mem_preimage, hT, mem_ofPred_eq] at hm
        simp only [mem_Ioc]
        constructor
        · rw [div_lt_iff₀ hα]; linarith [hm.1]
        · rw [le_div_iff₀ hα]; linarith [hm.2]
      have hle1 : (τ - β * d) / α ≤ (τ + β * d) / α :=
        div_le_div_of_nonneg_right (by nlinarith) hα.le
      calc _ ≤ (gaussianReal 0 1).real (Ioc ((τ - β * d) / α) ((τ + β * d) / α)) :=
            measureReal_mono hsub
        _ ≤ ((τ + β * d) / α - (τ - β * d) / α) / 2 := C3.gaussianReal_Ioc_le hle1
        _ = β / α * max d 0 := by
            rw [max_eq_left hd.le]; field_simp; ring
    · have hempty : Prod.mk d ⁻¹' T = ∅ := by
        ext m
        simp only [mem_preimage, hT, mem_ofPred_eq, mem_empty_iff_false, iff_false, not_and, not_le]
        intro h1
        have : d ≤ 0 := not_lt.mp hd
        nlinarith
      rw [hempty, max_eq_right (not_lt.mp hd)]
      simp
  calc _ ≤ ∫ d, β / α * max d 0 ∂gaussianReal 0 1 := by
        apply integral_mono_of_nonneg (ae_of_all _ fun d => measureReal_nonneg)
        · exact integrable_posPart_gaussian.const_mul _
        · exact ae_of_all _ hle
    _ = _ := by rw [integral_const_mul, integral_posPart_gaussian]

/-- Two unit vectors with inner product `r ∈ (-1,1)`: sum/difference decomposition. -/
theorem pair_le {p q : E} (hp : ⟪p, p⟫ = 1) (hq : ⟪q, q⟫ = 1) {r : ℝ} (hpq : ⟪p, q⟫ = r)
    (hr0 : -1 < r) (hr1 : r < 1) (τ : ℝ) :
    (stdGaussian E).real {x | τ < ⟪p, x⟫ ∧ ⟪q, x⟫ ≤ τ} ≤
      Real.sqrt ((1 - r) / 2) / Real.sqrt ((1 + r) / 2) * (Real.sqrt (2 * Real.pi))⁻¹ := by
  have hqp : ⟪q, p⟫ = r := by rw [real_inner_comm]; exact hpq
  have h1 : 0 < (1 + r) / 2 := by linarith
  have h2 : 0 < (1 - r) / 2 := by linarith
  set α := Real.sqrt ((1 + r) / 2) with hαdef
  set β := Real.sqrt ((1 - r) / 2) with hβdef
  have hα : 0 < α := Real.sqrt_pos.mpr h1
  have hβ : 0 < β := Real.sqrt_pos.mpr h2
  have hα2 : α * α = (1 + r) / 2 := Real.mul_self_sqrt h1.le
  have hβ2 : β * β = (1 - r) / 2 := Real.mul_self_sqrt h2.le
  set a : E := (2 * α)⁻¹ • (p + q) with hadef
  set b : E := (2 * β)⁻¹ • (p - q) with hbdef
  have ha : ⟪a, a⟫ = 1 := by
    simp only [hadef, real_inner_smul_left, real_inner_smul_right, inner_add_left,
      inner_add_right, hp, hq, hpq, hqp]
    field_simp
    nlinarith
  have hb : ⟪b, b⟫ = 1 := by
    simp only [hbdef, real_inner_smul_left, real_inner_smul_right, inner_sub_left,
      inner_sub_right, hp, hq, hpq, hqp]
    field_simp
    nlinarith
  have hab : ⟪a, b⟫ = 0 := by
    simp only [hadef, hbdef, real_inner_smul_left, real_inner_smul_right, inner_sub_right,
      inner_add_left, hp, hq, hpq, hqp]
    ring
  have hpe : α • a + β • b = p := by
    simp only [hadef, hbdef, smul_smul]
    have e1 : α * (2 * α)⁻¹ = 1 / 2 := by field_simp
    have e2 : β * (2 * β)⁻¹ = 1 / 2 := by field_simp
    rw [e1, e2]
    module
  have hqe : α • a - β • b = q := by
    simp only [hadef, hbdef, smul_smul]
    have e1 : α * (2 * α)⁻¹ = 1 / 2 := by field_simp
    have e2 : β * (2 * β)⁻¹ = 1 / 2 := by field_simp
    rw [e1, e2]
    module
  have := pair_core ha hb hab hα hβ τ
  rwa [hpe, hqe] at this

/-! ## Identification of the ideal 3-star -/

theorem idealStarFn_eq_star3 (τ : ℝ) :
    idealStarFn 3 r33 τ =
      (stdGaussian (Sample 3)).real
        {x | τ < ⟪star3Vec r33 0, x⟫ ∧ ∃ i : Fin 4, i ≠ 0 ∧ ⟪star3Vec r33 i, x⟫ ≤ τ} := by
  rw [idealStarFn, ← gram_star3Vec, bridge 3 4 (star3Vec r33) _ (measurableSet_frontierEvent 3 τ)]
  rfl

theorem idealStarFn_le_one (τ : ℝ) : idealStarFn 3 r33 τ ≤ 1 := by
  rw [idealStarFn_eq_star3]
  exact measureReal_le_one

theorem idealStarFn_nonneg (τ : ℝ) : 0 ≤ idealStarFn 3 r33 τ := by
  rw [idealStarFn_eq_star3]
  exact measureReal_nonneg

theorem inner_e0_e0 : ⟪(C3.Ideal.e 0 : Sample 3), C3.Ideal.e 0⟫ = 1 := by
  simpa using C3.Ideal.inner_e 0 0

theorem inner_bstar_self (j : Fin 3) : ⟪(C3.Ideal.bstar r33 j : Sample 3), C3.Ideal.bstar r33 j⟫ = 1 := by
  have := C3.Ideal.inner_bstar r33_lt_one.le (by linarith [r33_pos]) j j
  simpa using this

theorem sqrt_two_pi_inv_le : (Real.sqrt (2 * Real.pi))⁻¹ ≤ 2 / 5 := by
  have hpi : (3.14 : ℝ) < Real.pi := Real.pi_gt_d2
  have h : (5 / 2 : ℝ) ≤ Real.sqrt (2 * Real.pi) := by
    apply Real.le_sqrt_of_sq_le
    nlinarith
  calc (Real.sqrt (2 * Real.pi))⁻¹ ≤ (5 / 2 : ℝ)⁻¹ := inv_anti₀ (by norm_num) h
    _ = 2 / 5 := by norm_num

/-- Each leaf term of the union bound. -/
theorem leaf_le (j : Fin 3) (τ : ℝ) :
    (stdGaussian (Sample 3)).real
        {x | τ < ⟪(C3.Ideal.e 0 : Sample 3), x⟫ ∧ ⟪(C3.Ideal.bstar r33 j : Sample 3), x⟫ ≤ τ} ≤
      2 / 25 := by
  have hr := r33_sq
  have h := pair_le (E := Sample 3) inner_e0_e0 (inner_bstar_self j)
    (C3.Ideal.inner_e0_bstar r33 j) (by linarith [r33_pos]) r33_lt_one τ
  refine h.trans ?_
  have hc := sqrt_two_pi_inv_le
  have hratio : Real.sqrt ((1 - r33) / 2) / Real.sqrt ((1 + r33) / 2) ≤ 1 / 5 := by
    have h1 : 0 < (1 + r33) / 2 := by linarith [r33_pos]
    rw [div_le_iff₀ (Real.sqrt_pos.mpr h1)]
    have h12 : 12 / 13 ≤ r33 := by nlinarith [r33_pos]
    apply Real.sqrt_le_iff.mpr
    refine ⟨by positivity, ?_⟩
    have := Real.sq_sqrt h1.le
    nlinarith
  have hnn : 0 ≤ Real.sqrt ((1 - r33) / 2) / Real.sqrt ((1 + r33) / 2) := by positivity
  calc _ ≤ 1 / 5 * (2 / 5) := mul_le_mul hratio hc (by positivity) (by norm_num)
    _ = 2 / 25 := by norm_num

/-! ## The upper bound -/

theorem idealStarFn_le (τ : ℝ) : idealStarFn 3 r33 τ ≤ 6 / 25 := by
  rw [idealStarFn_eq_star3]
  have hsub : {x : Sample 3 | τ < ⟪star3Vec r33 0, x⟫ ∧
        ∃ i : Fin 4, i ≠ 0 ∧ ⟪star3Vec r33 i, x⟫ ≤ τ} ⊆
      ⋃ j : Fin 3, {x | τ < ⟪(C3.Ideal.e 0 : Sample 3), x⟫ ∧
        ⟪(C3.Ideal.bstar r33 j : Sample 3), x⟫ ≤ τ} := by
    rintro x ⟨h0, i, hi, hle⟩
    obtain ⟨j, rfl⟩ := Fin.exists_succ_eq.mpr hi
    simp only [star3Vec, Fin.cons_zero, Fin.cons_succ] at h0 hle
    exact mem_iUnion.mpr ⟨j, h0, hle⟩
  calc _ ≤ (stdGaussian (Sample 3)).real (⋃ j : Fin 3, {x : Sample 3 |
          τ < ⟪(C3.Ideal.e 0 : Sample 3), x⟫ ∧ ⟪(C3.Ideal.bstar r33 j : Sample 3), x⟫ ≤ τ}) :=
        measureReal_mono hsub
    _ ≤ ∑ j : Fin 3, (stdGaussian (Sample 3)).real {x : Sample 3 |
          τ < ⟪(C3.Ideal.e 0 : Sample 3), x⟫ ∧ ⟪(C3.Ideal.bstar r33 j : Sample 3), x⟫ ≤ τ} :=
        measureReal_iUnion_fintype_le _
    _ ≤ ∑ _j : Fin 3, (2 / 25 : ℝ) := Finset.sum_le_sum fun j _ => leaf_le j τ
    _ = 6 / 25 := by simp; norm_num

theorem Fcubic_le : Fcubic ≤ 103 / 400 := by
  unfold Fcubic
  refine csSup_le (Set.range_nonempty _) ?_
  rintro _ ⟨τ, rfl⟩
  exact (idealStarFn_le τ).trans (by norm_num)

/-! ## The lower bound -/

theorem lower_zero : Real.arccos r33 / (2 * Real.pi) ≤ idealStarFn 3 r33 0 := by
  rw [idealStarFn_eq_star3]
  set p : Sample 3 := C3.Ideal.e 0 with hp
  set q : Sample 3 := C3.Ideal.bstar r33 0 with hq
  have hpp : ⟪p, p⟫ = 1 := inner_e0_e0
  have hqq : ⟪q, q⟫ = 1 := inner_bstar_self 0
  have hpq : ⟪p, q⟫ = r33 := C3.Ideal.inner_e0_bstar r33 0
  have hm : ∀ (u v : Sample 3), MeasurableSet {x : Sample 3 | 0 < ⟪u, x⟫ ∧ 0 < ⟪v, x⟫} := by
    intro u v
    exact (measurableSet_lt (f := fun _ : Sample 3 => (0 : ℝ)) (g := fun x => ⟪u, x⟫)
        measurable_const (by fun_prop)).inter
      (measurableSet_lt (f := fun _ : Sample 3 => (0 : ℝ)) (g := fun x => ⟪v, x⟫)
        measurable_const (by fun_prop))
  have h1 := C3.bivariate_orthant hpp hpp
  have h2 := C3.bivariate_orthant hpp hqq
  rw [hpp, Real.arccos_one] at h1
  rw [hpq] at h2
  have hsub : {x : Sample 3 | 0 < ⟪p, x⟫ ∧ 0 < ⟪p, x⟫} \ {x | 0 < ⟪p, x⟫ ∧ 0 < ⟪q, x⟫} ⊆
      {x : Sample 3 | 0 < ⟪star3Vec r33 0, x⟫ ∧
        ∃ i : Fin 4, i ≠ 0 ∧ ⟪star3Vec r33 i, x⟫ ≤ 0} := by
    rintro x ⟨⟨h0, -⟩, hn⟩
    refine ⟨by simpa [star3Vec] using h0, (1 : Fin 4), by decide, ?_⟩
    simp only [star3Vec]
    by_contra hc
    exact hn ⟨h0, not_le.mp hc⟩
  have hinc : {x : Sample 3 | 0 < ⟪p, x⟫ ∧ 0 < ⟪q, x⟫} ⊆ {x | 0 < ⟪p, x⟫ ∧ 0 < ⟪p, x⟫} :=
    fun x hx => ⟨hx.1, hx.1⟩
  have hd := measureReal_sdiff (μ := stdGaussian (Sample 3)) hinc (hm p q)
  calc Real.arccos r33 / (2 * Real.pi)
      = (stdGaussian (Sample 3)).real {x : Sample 3 | 0 < ⟪p, x⟫ ∧ 0 < ⟪p, x⟫} -
          (stdGaussian (Sample 3)).real {x | 0 < ⟪p, x⟫ ∧ 0 < ⟪q, x⟫} := by
        rw [h1, h2]; ring
    _ = (stdGaussian (Sample 3)).real
          ({x : Sample 3 | 0 < ⟪p, x⟫ ∧ 0 < ⟪p, x⟫} \ {x | 0 < ⟪p, x⟫ ∧ 0 < ⟪q, x⟫}) := hd.symm
    _ ≤ _ := measureReal_mono hsub

theorem Fcubic_ge : 1 / 50 ≤ Fcubic := by
  have hb : BddAbove (Set.range fun τ : ℝ => idealStarFn 3 r33 τ) :=
    ⟨1, by rintro _ ⟨τ, rfl⟩; exact idealStarFn_le_one τ⟩
  refine le_trans ?_ (le_csSup hb ⟨0, rfl⟩)
  refine le_trans ?_ lower_zero
  have hpi : Real.pi < 3.15 := Real.pi_lt_d2
  have hpi0 : 0 < Real.pi := Real.pi_pos
  have hac := C3.lt_arccos_correlation
  have : r33 = 2 * Real.sqrt 2 / 3 := rfl
  rw [this, le_div_iff₀ (by positivity)]
  nlinarith


end T2
end CubicPathwidth

#print axioms CubicPathwidth.T2.Fcubic_ge
#print axioms CubicPathwidth.T2.Fcubic_le