CubicPathwidth.T2.SlepianSmooth

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

import CubicPathwidth.T2.SlepianStein

/-!
# Slepian's inequality for smooth product test functions

On `Fin N → ℝ` with the product standard normal law, let `A i` and `B i` be coefficient vectors
with every `A i` orthogonal to every `B l`, equal norms `‖A i‖ = ‖B i‖`, and
`⟨A i, A l⟩ ≤ ⟨B i, B l⟩` for `i ≠ l`. For a bounded, nondecreasing, twice differentiable `g`
with values in `[0, 1]`,
`E ∏ g(⟨A i, x⟩ - t i) ≤ E ∏ g(⟨B i, x⟩ - t i)`.

The proof interpolates along `c_l(θ) = cos θ • A l + sin θ • B l`, `θ ∈ [0, π/2]`, and
computes the derivative with the coordinatewise Stein identity.
-/

noncomputable section
open scoped Real
open MeasureTheory ProbabilityTheory Set Filter Topology Finset

namespace CubicPathwidth
namespace T2
namespace Slepian

/-- Coordinate inner product on `Fin N → ℝ`. -/
def ip {N : ℕ} (a x : Fin N → ℝ) : ℝ := ∑ j, a j * x j

theorem continuous_ip {N : ℕ} (a : Fin N → ℝ) : Continuous (ip a) := by
  unfold ip; fun_prop

theorem integrable_ip {N : ℕ} (a : Fin N → ℝ) : Integrable (ip a) (gpi N) := by
  unfold ip
  exact integrable_finsetSum _ fun j _ => (integrable_coord N j).const_mul (a j)

theorem ip_update {N : ℕ} (a x : Fin N → ℝ) (j : Fin N) (s : ℝ) :
    ip a (Function.update x j s) = ip a x + a j * (s - x j) := by
  unfold ip
  have : ∀ r, a r * Function.update x j s r = a r * x r + a r * (Pi.single j (s - x j) : Fin N → ℝ) r := by
    intro r
    by_cases h : r = j
    · subst h; simp; ring
    · simp [h]
  simp only [this, Finset.sum_add_distrib]
  congr 1
  rw [Finset.sum_eq_single j]
  · simp
  · intro r _ hr; simp [hr]
  · simp

variable {k : ℕ} (g g1 g2 : ℝ → ℝ)

/-- First partial derivative of `∏ g(z l)` in the `i`-th coordinate. -/
def P (i : Fin k) (z : Fin k → ℝ) : ℝ := (∏ l ∈ univ.erase i, g (z l)) * g1 (z i)

/-- Mixed second partial derivative (for `l ≠ i`). -/
def Q (i l : Fin k) (z : Fin k → ℝ) : ℝ :=
  (∏ r ∈ (univ.erase i).erase l, g (z r)) * g1 (z l) * g1 (z i)

/-- Diagonal second partial derivative. -/
def R (i : Fin k) (z : Fin k → ℝ) : ℝ := (∏ l ∈ univ.erase i, g (z l)) * g2 (z i)

variable {g g1 g2}

theorem hasDerivAt_P_line (hg : ∀ s, HasDerivAt g (g1 s) s) (hg1 : ∀ s, HasDerivAt g1 (g2 s) s)
    (i : Fin k) (u w : Fin k → ℝ) (s0 : ℝ) :
    HasDerivAt (fun s => P g g1 i (fun l => u l + s * w l))
      (∑ l ∈ univ.erase i, w l * Q g g1 i l (fun l => u l + s0 * w l) +
        w i * R g g2 i (fun l => u l + s0 * w l)) s0 := by
  have hlin : ∀ l, HasDerivAt (fun s : ℝ => u l + s * w l) (w l) s0 := by
    intro l
    simpa using ((hasDerivAt_id s0).mul_const (w l)).const_add (u l)
  have hl : ∀ l ∈ univ.erase i, HasDerivAt (fun s => g (u l + s * w l))
      (g1 (u l + s0 * w l) * w l) s0 := fun l _ => (hg _).comp s0 (hlin l)
  have hprod := HasDerivAt.fun_finsetProd hl
  have hi : HasDerivAt (fun s => g1 (u i + s * w i)) (g2 (u i + s0 * w i) * w i) s0 :=
    (hg1 _).comp s0 (hlin i)
  have := hprod.mul hi
  unfold P
  convert this using 1
  unfold Q R
  simp only [smul_eq_mul, Finset.sum_mul]
  congr 1
  · refine Finset.sum_congr rfl fun l _ => ?_
    ring
  · ring

section bounds

variable (h0 : ∀ s, 0 ≤ g s) (h1 : ∀ s, g s ≤ 1) (hg1n : ∀ s, 0 ≤ g1 s)
include h0 h1

theorem prod_mem01 (s : Finset (Fin k)) (z : Fin k → ℝ) :
    0 ≤ ∏ l ∈ s, g (z l) ∧ ∏ l ∈ s, g (z l) ≤ 1 :=
  ⟨Finset.prod_nonneg fun _ _ => h0 _, Finset.prod_le_one₀ (fun _ _ => h0 _) fun _ _ => h1 _⟩

include hg1n

theorem P_bound {C1 : ℝ} (hC1 : ∀ s, g1 s ≤ C1) (i : Fin k) (z : Fin k → ℝ) :
    |P g g1 i z| ≤ C1 := by
  obtain ⟨a, b⟩ := prod_mem01 h0 h1 (univ.erase i) z
  unfold P
  rw [abs_of_nonneg (mul_nonneg a (hg1n _))]
  calc _ ≤ 1 * g1 (z i) := mul_le_mul_of_nonneg_right b (hg1n _)
    _ ≤ C1 := by rw [one_mul]; exact hC1 _

theorem Q_nonneg (i l : Fin k) (z : Fin k → ℝ) : 0 ≤ Q g g1 i l z := by
  unfold Q
  exact mul_nonneg (mul_nonneg (prod_mem01 h0 h1 _ z).1 (hg1n _)) (hg1n _)

theorem Q_bound {C1 : ℝ} (hC1 : ∀ s, g1 s ≤ C1) (i l : Fin k) (z : Fin k → ℝ) :
    |Q g g1 i l z| ≤ C1 * C1 := by
  rw [abs_of_nonneg (Q_nonneg h0 h1 hg1n i l z)]
  obtain ⟨a, b⟩ := prod_mem01 h0 h1 ((univ.erase i).erase l) z
  unfold Q
  have c0 : 0 ≤ C1 := (hg1n 0).trans (hC1 0)
  calc _ ≤ 1 * C1 * C1 := by
        refine mul_le_mul (mul_le_mul b (hC1 _) (hg1n _) zero_le_one) (hC1 _) (hg1n _) ?_
        exact mul_nonneg zero_le_one c0
    _ = C1 * C1 := by ring

omit hg1n in

theorem R_bound {C2 : ℝ} (hC2 : ∀ s, |g2 s| ≤ C2) (i : Fin k) (z : Fin k → ℝ) :
    |R g g2 i z| ≤ C2 := by
  obtain ⟨a, b⟩ := prod_mem01 h0 h1 (univ.erase i) z
  unfold R
  rw [abs_mul, abs_of_nonneg a]
  calc _ ≤ 1 * C2 := mul_le_mul b (hC2 _) (abs_nonneg _) zero_le_one
    _ = C2 := one_mul _

end bounds

theorem integrable_of_bdd {N : ℕ} {f : (Fin N → ℝ) → ℝ} (hf : Continuous f) {C : ℝ}
    (hb : ∀ x, |f x| ≤ C) : Integrable f (gpi N) :=
  Integrable.of_bound hf.aestronglyMeasurable C (ae_of_all _ fun x => by simpa using hb x)

/-- Stein identity for `P i` composed with affine maps. -/
theorem stein_P {n : ℕ} (hg : ∀ s, HasDerivAt g (g1 s) s) (hg1 : ∀ s, HasDerivAt g1 (g2 s) s)
    (hg2 : Continuous g2) (h0 : ∀ s, 0 ≤ g s) (h1 : ∀ s, g s ≤ 1) (hg1n : ∀ s, 0 ≤ g1 s)
    {C1 C2 : ℝ} (hC1 : ∀ s, g1 s ≤ C1) (hC2 : ∀ s, |g2 s| ≤ C2)
    (c : Fin k → Fin (n + 1) → ℝ) (t : Fin k → ℝ) (d : Fin (n + 1) → ℝ) (i : Fin k) :
    ∫ x, P g g1 i (fun l => ip (c l) x - t l) * ip d x ∂gpi (n + 1) =
      ∑ l ∈ univ.erase i, ip d (c l) * ∫ x, Q g g1 i l (fun l => ip (c l) x - t l) ∂gpi (n + 1) +
        ip d (c i) * ∫ x, R g g2 i (fun l => ip (c l) x - t l) ∂gpi (n + 1) := by
  have hgc : Continuous g := continuous_iff_continuousAt.2 fun s => (hg s).continuousAt
  have hg1c : Continuous g1 := continuous_iff_continuousAt.2 fun s => (hg1 s).continuousAt
  set Z : (Fin (n + 1) → ℝ) → Fin k → ℝ := fun x l => ip (c l) x - t l with hZ
  have hZc : Continuous Z := by
    rw [hZ]; exact continuous_pi fun l => (continuous_ip _).sub continuous_const
  have hPc : Continuous fun x => P g g1 i (Z x) := by
    unfold P
    exact (continuous_finsetProd _ fun l _ => hgc.comp ((continuous_apply l).comp hZc)).mul
      (hg1c.comp ((continuous_apply i).comp hZc))
  have hQc : ∀ l, Continuous fun x => Q g g1 i l (Z x) := by
    intro l; unfold Q
    exact ((continuous_finsetProd _ fun r _ => hgc.comp ((continuous_apply r).comp hZc)).mul
      (hg1c.comp ((continuous_apply l).comp hZc))).mul (hg1c.comp ((continuous_apply i).comp hZc))
  have hRc : Continuous fun x => R g g2 i (Z x) := by
    unfold R
    exact (continuous_finsetProd _ fun l _ => hgc.comp ((continuous_apply l).comp hZc)).mul
      (hg2.comp ((continuous_apply i).comp hZc))
  have hQi : ∀ l, Integrable (fun x => Q g g1 i l (Z x)) (gpi (n + 1)) := fun l =>
    integrable_of_bdd (hQc l) fun x => Q_bound h0 h1 hg1n hC1 i l (Z x)
  have hRi : Integrable (fun x => R g g2 i (Z x)) (gpi (n + 1)) :=
    integrable_of_bdd hRc fun x => R_bound h0 h1 hC2 i (Z x)
  -- the coordinate derivative
  set H' : Fin (n + 1) → (Fin (n + 1) → ℝ) → ℝ := fun j x =>
    ∑ l ∈ univ.erase i, c l j * Q g g1 i l (Z x) + c i j * R g g2 i (Z x) with hH'
  have hstein : ∀ j, ∫ x, x j * P g g1 i (Z x) ∂gpi (n + 1) = ∫ x, H' j x ∂gpi (n + 1) := by
    intro j
    refine stein_coord j hPc ?_ (fun x => P_bound h0 h1 hg1n hC1 i (Z x))
      (C' := ∑ l ∈ univ.erase i, |c l j| * (C1 * C1) + |c i j| * C2) ?_ ?_
    · rw [hH']
      exact ((continuous_finsetSum _ fun l _ => continuous_const.mul (hQc l)).add
        (continuous_const.mul hRc)).measurable
    · intro x
      rw [hH']
      refine (abs_add_le _ _).trans (add_le_add ((abs_sum_le_sum_abs _ _).trans
        (Finset.sum_le_sum fun l _ => ?_)) ?_)
      · rw [abs_mul]; exact mul_le_mul_of_nonneg_left (Q_bound h0 h1 hg1n hC1 i l _) (abs_nonneg _)
      · rw [abs_mul]; exact mul_le_mul_of_nonneg_left (R_bound h0 h1 hC2 i _) (abs_nonneg _)
    · intro x
      have := hasDerivAt_P_line hg hg1 i (fun l => ip (c l) x - t l - c l j * x j)
        (fun l => c l j) (x j)
      have e1 : (fun l => ip (c l) x - t l - c l j * x j + x j * c l j) = Z x := by
        funext l; rw [hZ]; ring
      have e2 : (fun s => P g g1 i (fun l => ip (c l) x - t l - c l j * x j + s * c l j)) =
          fun s => P g g1 i (Z (Function.update x j s)) := by
        funext s; congr 1; funext l; rw [hZ]; simp only [ip_update]; ring
      rw [e1, e2] at this
      exact this
  -- assemble
  change ∫ x, P g g1 i (Z x) * ip d x ∂gpi (n + 1) =
    ∑ l ∈ univ.erase i, ip d (c l) * ∫ x, Q g g1 i l (Z x) ∂gpi (n + 1) +
      ip d (c i) * ∫ x, R g g2 i (Z x) ∂gpi (n + 1)
  have hd : ∀ x, P g g1 i (Z x) * ip d x = ∑ j, d j * (x j * P g g1 i (Z x)) := by
    intro x; unfold ip; rw [Finset.mul_sum]
    exact Finset.sum_congr rfl fun j _ => by ring
  simp_rw [hd]
  rw [integral_finsetSum _ fun j _ =>
    (((integrable_coord _ j).mul_bdd hPc.aestronglyMeasurable
      (ae_of_all _ fun x => by simpa using P_bound h0 h1 hg1n hC1 i (Z x))).const_mul (d j))]
  simp_rw [integral_const_mul, hstein]
  have hH'int : ∀ j, ∫ x, H' j x ∂gpi (n + 1) =
      ∑ l ∈ univ.erase i, c l j * ∫ x, Q g g1 i l (Z x) ∂gpi (n + 1) +
        c i j * ∫ x, R g g2 i (Z x) ∂gpi (n + 1) := by
    intro j
    rw [hH', integral_add (integrable_finsetSum _ fun l _ => (hQi l).const_mul _)
      (hRi.const_mul _), integral_finsetSum _ fun l _ => (hQi l).const_mul _]
    simp_rw [integral_const_mul]
  simp_rw [hH'int, mul_add, Finset.sum_add_distrib, Finset.mul_sum]
  congr 1
  · rw [Finset.sum_comm]
    refine Finset.sum_congr rfl fun l _ => ?_
    unfold ip
    rw [Finset.sum_mul]
    exact Finset.sum_congr rfl fun j _ => by ring
  · unfold ip
    rw [Finset.sum_mul]
    exact Finset.sum_congr rfl fun j _ => by ring


section interp

variable {N : ℕ} (A B : Fin k → Fin N → ℝ)

/-- Interpolating coefficient vectors `cos θ • A l + sin θ • B l`. -/
def cvec (θ : ℝ) (l : Fin k) : Fin N → ℝ := fun j => Real.cos θ * A l j + Real.sin θ * B l j

/-- Their `θ`-derivatives. -/
def dvec (θ : ℝ) (l : Fin k) : Fin N → ℝ := fun j => -Real.sin θ * A l j + Real.cos θ * B l j

theorem ip_comm (a x : Fin N → ℝ) : ip a x = ip x a := by
  unfold ip; exact Finset.sum_congr rfl fun j _ => mul_comm _ _

theorem ip_cvec (θ : ℝ) (l : Fin k) (x : Fin N → ℝ) :
    ip (cvec A B θ l) x = Real.cos θ * ip (A l) x + Real.sin θ * ip (B l) x := by
  unfold ip cvec
  rw [Finset.mul_sum, Finset.mul_sum, ← Finset.sum_add_distrib]
  exact Finset.sum_congr rfl fun j _ => by ring

theorem ip_dvec (θ : ℝ) (l : Fin k) (x : Fin N → ℝ) :
    ip (dvec A B θ l) x = -Real.sin θ * ip (A l) x + Real.cos θ * ip (B l) x := by
  unfold ip dvec
  rw [Finset.mul_sum, Finset.mul_sum, ← Finset.sum_add_distrib]
  exact Finset.sum_congr rfl fun j _ => by ring

theorem ip_dvec_cvec (horth : ∀ i l, ip (A i) (B l) = 0) (θ : ℝ) (i l : Fin k) :
    ip (dvec A B θ i) (cvec A B θ l) =
      Real.sin θ * Real.cos θ * (ip (B i) (B l) - ip (A i) (A l)) := by
  rw [ip_comm, ip_cvec, ip_comm (A l), ip_comm (B l), ip_dvec, ip_dvec, horth i l,
    ip_comm (B i) (A l), horth l i]
  ring

end interp

/-- Slepian's inequality for smooth product test functions, with `A` orthogonal to `B`. -/
theorem smooth_slepian {N : ℕ} (A B : Fin k → Fin N → ℝ) (horth : ∀ i l, ip (A i) (B l) = 0)
    (hnorm : ∀ i, ip (A i) (A i) = ip (B i) (B i))
    (hoff : ∀ i l, i ≠ l → ip (A i) (A l) ≤ ip (B i) (B l))
    (hg : ∀ s, HasDerivAt g (g1 s) s) (hg1 : ∀ s, HasDerivAt g1 (g2 s) s)
    (hg2 : Continuous g2) (h0 : ∀ s, 0 ≤ g s) (h1 : ∀ s, g s ≤ 1) (hg1n : ∀ s, 0 ≤ g1 s)
    {C1 C2 : ℝ} (hC1 : ∀ s, g1 s ≤ C1) (hC2 : ∀ s, |g2 s| ≤ C2) (t : Fin k → ℝ) :
    ∫ x, ∏ i, g (ip (A i) x - t i) ∂gpi N ≤ ∫ x, ∏ i, g (ip (B i) x - t i) ∂gpi N := by
  cases N with
  | zero => simp [ip]
  | succ n =>
  have hgc : Continuous g := continuous_iff_continuousAt.2 fun s => (hg s).continuousAt
  have hg1c : Continuous g1 := continuous_iff_continuousAt.2 fun s => (hg1 s).continuousAt
  set Z : ℝ → (Fin (n + 1) → ℝ) → Fin k → ℝ := fun θ x l => ip (cvec A B θ l) x - t l with hZ
  have hZc : ∀ θ, Continuous (Z θ) := fun θ =>
    continuous_pi fun l => (continuous_ip _).sub continuous_const
  set F : ℝ → ℝ := fun θ => ∫ x, ∏ l, g (Z θ x l) ∂gpi (n + 1) with hF
  set F' : ℝ → ℝ := fun θ =>
    ∫ x, ∑ i, P g g1 i (Z θ x) * ip (dvec A B θ i) x ∂gpi (n + 1) with hF'
  have hPc : ∀ θ i, Continuous fun x => P g g1 i (Z θ x) := by
    intro θ i
    unfold P
    exact (continuous_finsetProd _ fun l _ => hgc.comp ((continuous_apply l).comp (hZc θ))).mul
      (hg1c.comp ((continuous_apply i).comp (hZc θ)))
  have hPdi : ∀ θ i, Integrable (fun x => P g g1 i (Z θ x) * ip (dvec A B θ i) x)
      (gpi (n + 1)) := fun θ i =>
    (integrable_ip _).bdd_mul (hPc θ i).aestronglyMeasurable
      (ae_of_all _ fun x => by simpa using P_bound h0 h1 hg1n hC1 i (Z θ x))
  have hderiv : ∀ θ0, HasDerivAt F (F' θ0) θ0 := by
    intro θ0
    have key := hasDerivAt_integral_of_dominated_loc_of_deriv_le (μ := gpi (n + 1))
      (F := fun θ x => ∏ l, g (Z θ x l))
      (F' := fun θ x => ∑ i, P g g1 i (Z θ x) * ip (dvec A B θ i) x) (x₀ := θ0)
      (bound := fun x => ∑ i, C1 * (|ip (A i) x| + |ip (B i) x|)) (s := univ) univ_mem
      (Eventually.of_forall fun θ =>
        (continuous_finsetProd _ fun l _ =>
          hgc.comp ((continuous_apply l).comp (hZc θ))).aestronglyMeasurable)
      (integrable_of_bdd (continuous_finsetProd _ fun l _ =>
          hgc.comp ((continuous_apply l).comp (hZc θ0))) (C := 1) fun x => by
        obtain ⟨a, b⟩ := prod_mem01 h0 h1 univ (Z θ0 x)
        rw [abs_of_nonneg a]; exact b)
      (integrable_finsetSum _ fun i _ => hPdi θ0 i).aestronglyMeasurable
      (ae_of_all _ fun x θ _ => by
        refine (norm_sum_le _ _).trans (Finset.sum_le_sum fun i _ => ?_)
        rw [Real.norm_eq_abs, abs_mul, ip_dvec]
        have hd : |-Real.sin θ * ip (A i) x + Real.cos θ * ip (B i) x| ≤
            |ip (A i) x| + |ip (B i) x| := by
          refine (abs_add_le _ _).trans (add_le_add ?_ ?_)
          · rw [abs_mul, abs_neg]
            exact mul_le_of_le_one_left (abs_nonneg _) (Real.abs_sin_le_one θ)
          · rw [abs_mul]
            exact mul_le_of_le_one_left (abs_nonneg _) (Real.abs_cos_le_one θ)
        exact mul_le_mul (P_bound h0 h1 hg1n hC1 i _) hd (abs_nonneg _)
          ((abs_nonneg _).trans (P_bound h0 h1 hg1n hC1 i (Z θ x))))
      (integrable_finsetSum _ fun i _ =>
        (((integrable_ip _).abs.add (integrable_ip _).abs)).const_mul C1)
      (ae_of_all _ fun x θ _ => by
        have hl : ∀ l ∈ (univ : Finset (Fin k)), HasDerivAt (fun θ => g (Z θ x l))
            (g1 (Z θ x l) * ip (dvec A B θ l) x) θ := by
          intro l _
          have e : (fun θ => Z θ x l) =
              fun θ => Real.cos θ * ip (A l) x + Real.sin θ * ip (B l) x - t l := by
            funext θ; rw [hZ]; simp only [ip_cvec]
          have hlin : HasDerivAt (fun θ => Z θ x l) (ip (dvec A B θ l) x) θ := by
            rw [e, ip_dvec]
            have := (((Real.hasDerivAt_cos θ).mul_const (ip (A l) x)).add
              ((Real.hasDerivAt_sin θ).mul_const (ip (B l) x))).sub_const (t l)
            convert this using 1
          exact (hg _).comp θ hlin
        have := HasDerivAt.fun_finsetProd hl
        convert this using 1
        refine Finset.sum_congr rfl fun i _ => ?_
        unfold P; rw [smul_eq_mul]; ring)
    exact key.2
  have hnn : ∀ θ ∈ Icc 0 (π / 2), 0 ≤ F' θ := by
    intro θ hθ
    have hsc : 0 ≤ Real.sin θ * Real.cos θ :=
      mul_nonneg (Real.sin_nonneg_of_nonneg_of_le_pi hθ.1 (by linarith [hθ.2, Real.pi_pos]))
        (Real.cos_nonneg_of_mem_Icc ⟨by linarith [hθ.1, Real.pi_pos], hθ.2⟩)
    simp only [hF']
    rw [integral_finsetSum _ fun i _ => hPdi θ i]
    refine Finset.sum_nonneg fun i _ => ?_
    have := stein_P hg hg1 hg2 h0 h1 hg1n hC1 hC2 (cvec A B θ) t (dvec A B θ i) i
    rw [this, ip_dvec_cvec A B horth, hnorm i, sub_self, mul_zero, zero_mul, add_zero]
    refine Finset.sum_nonneg fun l hl => ?_
    rw [ip_dvec_cvec A B horth]
    refine mul_nonneg (mul_nonneg hsc (sub_nonneg.2 (hoff i l (Finset.ne_of_mem_erase hl).symm)))
      (integral_nonneg fun x => Q_nonneg h0 h1 hg1n i l _)
  have mono : MonotoneOn F (Icc 0 (π / 2)) :=
    monotoneOn_of_hasDerivWithinAt_nonneg (convex_Icc _ _)
      (HasDerivAt.continuousOn fun θ _ => hderiv θ)
      (fun θ _ => (hderiv θ).hasDerivWithinAt) fun θ hθ => hnn θ (interior_subset hθ)
  have hpi : (0 : ℝ) ≤ π / 2 := by linarith [Real.pi_pos]
  have := mono (left_mem_Icc.2 hpi) (right_mem_Icc.2 hpi) hpi
  simp only [hF, hZ, ip_cvec, Real.cos_zero, Real.sin_zero, Real.cos_pi_div_two,
    Real.sin_pi_div_two, one_mul, zero_mul, add_zero, zero_add] at this
  exact this

end Slepian
end T2
end CubicPathwidth