CubicPathwidth.T2.Slepian

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

import CubicPathwidth.T2.SlepianSmooth
import Mathlib.Analysis.SpecialFunctions.SmoothTransition

/-!
# Slepian's inequality for Gaussian upper orthants (O3)

`slepian : SlepianClaim`.

Route: Gaussian interpolation with smooth product test functions
(`CubicPathwidth.T2.SlepianSmooth`), followed by

* approximation of the indicator of `(0, ∞)` by `s ↦ smoothTransition ((n + 1) s)`, which is
  eventually constant at every point, so dominated convergence applies;
* reduction of the inner-product model on `Sample m` to the coordinate model on
  `Fin (m + m) → ℝ`, with the `a`-vectors on the first block and the `b`-vectors on the second
  block, so that the two families are orthogonal.
-/

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

namespace CubicPathwidth
namespace T2
namespace Slepian

/-! ## The smooth step and its derivatives -/

local notation "ST" => Real.smoothTransition

theorem st_smooth : Differentiable ℝ ST ∧ Differentiable ℝ (deriv ST) ∧
    Continuous (deriv (deriv ST)) := by
  have h0 : ContDiff ℝ ∞ ST := Real.smoothTransition.contDiff (n := ⊤)
  have h1 := contDiff_infty_iff_deriv.1 h0
  have h2 := contDiff_infty_iff_deriv.1 h1.2
  exact ⟨h1.1, h2.1, h2.2.continuous⟩

theorem st1_zero : ∀ x : ℝ, x < 0 ∨ 1 < x → deriv ST x = 0 := by
  rintro x (hx | hx)
  · have : ST =ᶠ[𝓝 x] fun _ => 0 := by
      filter_upwards [Iio_mem_nhds hx] with y hy
      exact Real.smoothTransition.zero_of_nonpos hy.le
    rw [this.deriv_eq]; simp
  · have : ST =ᶠ[𝓝 x] fun _ => 1 := by
      filter_upwards [Ioi_mem_nhds hx] with y hy
      exact Real.smoothTransition.one_of_one_le hy.le
    rw [this.deriv_eq]; simp

theorem st2_zero : ∀ x : ℝ, x < 0 ∨ 1 < x → deriv (deriv ST) x = 0 := by
  rintro x (hx | hx)
  · have : deriv ST =ᶠ[𝓝 x] fun _ => 0 := by
      filter_upwards [Iio_mem_nhds hx] with y hy
      exact st1_zero y (Or.inl hy)
    rw [this.deriv_eq]; simp
  · have : deriv ST =ᶠ[𝓝 x] fun _ => 0 := by
      filter_upwards [Ioi_mem_nhds hx] with y hy
      exact st1_zero y (Or.inr hy)
    rw [this.deriv_eq]; simp

theorem bdd_of_vanish {h : ℝ → ℝ} (hc : Continuous h) (hz : ∀ x, x < 0 ∨ 1 < x → h x = 0) :
    ∃ C, ∀ x, |h x| ≤ C := by
  obtain ⟨C, hC⟩ := (isCompact_Icc (a := (0 : ℝ)) (b := 1)).exists_bound_of_continuousOn
    hc.continuousOn
  refine ⟨max C 0, fun x => ?_⟩
  by_cases hx : x ∈ Icc (0 : ℝ) 1
  · exact (by simpa [Real.norm_eq_abs] using hC x hx : |h x| ≤ C).trans (le_max_left _ _)
  · have : x < 0 ∨ 1 < x := by
      simp only [Set.mem_Icc, not_and_or, not_le] at hx; exact hx
    rw [hz x this]; simp

/-- Scaled smooth step. -/
def gs (c s : ℝ) : ℝ := ST (c * s)

/-- Its derivative. -/
def gs1 (c s : ℝ) : ℝ := c * deriv ST (c * s)

/-- Its second derivative. -/
def gs2 (c s : ℝ) : ℝ := c * (c * deriv (deriv ST) (c * s))

theorem hasDerivAt_gs (c s : ℝ) : HasDerivAt (gs c) (gs1 c s) s := by
  have := (st_smooth.1 (c * s)).hasDerivAt.comp s ((hasDerivAt_id s).const_mul c)
  unfold gs gs1
  convert this using 1 <;> first | rfl | ring

theorem hasDerivAt_gs1 (c s : ℝ) : HasDerivAt (gs1 c) (gs2 c s) s := by
  have := ((st_smooth.2.1 (c * s)).hasDerivAt.comp s ((hasDerivAt_id s).const_mul c)).const_mul c
  unfold gs1 gs2
  convert this using 1 <;> first | rfl | ring

theorem continuous_gs2 (c : ℝ) : Continuous (gs2 c) := by
  unfold gs2
  exact continuous_const.mul (continuous_const.mul
    (st_smooth.2.2.comp (continuous_const.mul continuous_id)))

/-! ## The coordinate model -/

theorem measurableSet_orthant {N k : ℕ} (A : Fin k → Fin N → ℝ) (t : Fin k → ℝ) :
    MeasurableSet {x | ∀ i, t i < ip (A i) x} := by
  simp only [Set.ofPred_forall]
  exact MeasurableSet.iInter fun i =>
    measurableSet_lt measurable_const (continuous_ip _).measurable

theorem tendsto_smooth_orthant {N k : ℕ} (A : Fin k → Fin N → ℝ) (t : Fin k → ℝ) :
    Tendsto (fun n : ℕ => ∫ x, ∏ i, gs ((n : ℝ) + 1) (ip (A i) x - t i) ∂gpi N) atTop
      (𝓝 ((gpi N).real {x | ∀ i, t i < ip (A i) x})) := by
  rw [← integral_indicator_one (measurableSet_orthant A t)]
  refine tendsto_integral_of_dominated_convergence (fun _ => 1) (fun n => ?_)
    (integrable_const 1) (fun n => ae_of_all _ fun x => ?_) (ae_of_all _ fun x => ?_)
  · unfold gs
    exact (continuous_finsetProd _ fun i _ => Real.smoothTransition.continuous.comp
      (continuous_const.mul ((continuous_ip _).sub continuous_const))).aestronglyMeasurable
  · simp only [gs]
    rw [Real.norm_eq_abs, abs_of_nonneg (Finset.prod_nonneg fun i _ =>
      Real.smoothTransition.nonneg _)]
    exact Finset.prod_le_one₀ (fun i _ => Real.smoothTransition.nonneg _)
      fun i _ => Real.smoothTransition.le_one _
  · have hev : ∀ i, ∀ᶠ n : ℕ in atTop,
        gs ((n : ℝ) + 1) (ip (A i) x - t i) = if t i < ip (A i) x then 1 else 0 := by
      intro i
      split_ifs with h
      · have hp : 0 < ip (A i) x - t i := sub_pos.2 h
        obtain ⟨n0, hn0⟩ := exists_nat_ge (1 / (ip (A i) x - t i))
        filter_upwards [eventually_ge_atTop n0] with n hn
        apply Real.smoothTransition.one_of_one_le
        have h1 : (1 / (ip (A i) x - t i)) ≤ (n : ℝ) + 1 := by
          have : (n0 : ℝ) ≤ n := by exact_mod_cast hn
          linarith
        rw [div_le_iff₀ hp] at h1
        exact h1
      · refine Eventually.of_forall fun n => Real.smoothTransition.zero_of_nonpos ?_
        exact mul_nonpos_of_nonneg_of_nonpos (by positivity) (by linarith [not_lt.1 h])
    refine tendsto_const_nhds.congr' ?_
    filter_upwards [eventually_all.2 hev] with n hn
    rw [Finset.prod_congr rfl fun i _ => hn i, Fintype.prod_boole]
    simp only [Set.indicator_apply, Set.mem_ofPred_eq, Pi.one_apply]
    convert rfl

/-- Slepian's inequality in the coordinate model, with `A` orthogonal to `B`. -/
theorem pi_slepian {N k : ℕ} (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)) (t : Fin k → ℝ) :
    (gpi N).real {x | ∀ i, t i < ip (A i) x} ≤ (gpi N).real {x | ∀ i, t i < ip (B i) x} := by
  obtain ⟨C1, hC1⟩ := bdd_of_vanish (st_smooth.2.1.continuous) st1_zero
  obtain ⟨C2, hC2⟩ := bdd_of_vanish st_smooth.2.2 st2_zero
  refine le_of_tendsto_of_tendsto' (tendsto_smooth_orthant A t) (tendsto_smooth_orthant B t)
    fun n => ?_
  set c : ℝ := (n : ℝ) + 1
  have hc : 0 ≤ c := by positivity
  refine smooth_slepian A B horth hnorm hoff (hasDerivAt_gs c) (hasDerivAt_gs1 c)
    (continuous_gs2 c) (fun s => Real.smoothTransition.nonneg _)
    (fun s => Real.smoothTransition.le_one _)
    (fun s => mul_nonneg hc (Real.smoothTransition.monotone.deriv_nonneg))
    (C1 := c * C1) (C2 := c * (c * C2)) (fun s => ?_) (fun s => ?_) t
  · exact mul_le_mul_of_nonneg_left ((le_abs_self _).trans (hC1 _)) hc
  · unfold gs2
    rw [abs_mul, abs_mul, abs_of_nonneg hc]
    exact mul_le_mul_of_nonneg_left (mul_le_mul_of_nonneg_left (hC2 _) hc) hc

/-! ## Reduction from the inner-product model -/

theorem map_castAdd (m m' : ℕ) :
    (gpi (m + m')).map (fun x (j : Fin m) => x (Fin.castAdd m' j)) = gpi m := by
  refine (Measure.pi_eq fun s hs => ?_).symm
  rw [Measure.map_apply (by fun_prop) (MeasurableSet.univ_pi hs)]
  have : (fun x (j : Fin m) => x (Fin.castAdd m' j)) ⁻¹' Set.pi univ s =
      Set.pi univ (Fin.append s fun _ => univ) := by
    ext x
    simp [Fin.forall_fin_add]
  rw [this, Measure.pi_pi, Fin.prod_univ_add]
  simp

theorem map_natAdd (m m' : ℕ) :
    (gpi (m + m')).map (fun x (j : Fin m') => x (Fin.natAdd m j)) = gpi m' := by
  refine (Measure.pi_eq fun s hs => ?_).symm
  rw [Measure.map_apply (by fun_prop) (MeasurableSet.univ_pi hs)]
  have : (fun x (j : Fin m') => x (Fin.natAdd m j)) ⁻¹' Set.pi univ s =
      Set.pi univ (Fin.append (fun _ => univ) s) := by
    ext x
    simp [Fin.forall_fin_add]
  rw [this, Measure.pi_pi, Fin.prod_univ_add]
  simp

theorem inner_eq_ip {m : ℕ} (u v : GaussianField.Sample m) :
    ⟪u, v⟫ = ip (fun j => u j) (fun j => v j) := by
  simp [PiLp.inner_apply, ip, mul_comm]

/-- Transfer of an orthant probability along a block embedding. -/
theorem orthant_transfer {m M k : ℕ} (a : Fin k → GaussianField.Sample m)
    (e : (Fin M → ℝ) → (Fin m → ℝ)) (he : Measurable e) (hmap : (gpi M).map e = gpi m)
    (A : Fin k → Fin M → ℝ) (hA : ∀ i x, ip (A i) x = ip (fun j => a i j) (e x))
    (t : Fin k → ℝ) :
    (stdGaussian (GaussianField.Sample m)).real {x | ∀ i, t i < ⟪a i, x⟫} =
      (gpi M).real {x | ∀ i, t i < ip (A i) x} := by
  have hS : MeasurableSet {x : GaussianField.Sample m | ∀ i, t i < ⟪a i, x⟫} := by
    simp only [Set.ofPred_forall]
    exact MeasurableSet.iInter fun i =>
      measurableSet_lt measurable_const (continuous_const.inner continuous_id).measurable
  have hmeas : Measurable fun x => WithLp.toLp 2 (e x) := by fun_prop
  have hstd : (gpi M).map (fun x => WithLp.toLp 2 (e x)) =
      stdGaussian (GaussianField.Sample m) := by
    rw [← map_pi_eq_stdGaussian]
    change _ = Measure.map (WithLp.toLp 2) (gpi m)
    rw [← hmap, Measure.map_map (by fun_prop) he]
    rfl
  rw [← hstd, map_measureReal_apply hmeas hS]
  congr 1
  ext x
  simp only [Set.mem_preimage, Set.mem_ofPred_eq, inner_eq_ip, hA]

theorem ip_append {m m' : ℕ} (u v : Fin m → ℝ) (u' v' : Fin m' → ℝ) :
    ip (Fin.append u u') (Fin.append v v') = ip u v + ip u' v' := by
  simp [ip, Fin.sum_univ_add]

end Slepian

open Slepian in

/-- **Slepian's inequality** for Gaussian upper orthants (O3). -/
theorem slepian : SlepianClaim := by
  intro m k a b hn ho t
  set A : Fin k → Fin (m + m) → ℝ := fun i => Fin.append (fun j => a i j) (0 : Fin m → ℝ)
    with hAdef
  set B : Fin k → Fin (m + m) → ℝ := fun i => Fin.append (0 : Fin m → ℝ) (fun j => b i j)
    with hBdef
  have split : ∀ x : Fin (m + m) → ℝ,
      x = Fin.append (fun j => x (Fin.castAdd m j)) (fun j => x (Fin.natAdd m j)) := by
    intro x; ext j
    refine Fin.addCases (fun j => ?_) (fun j => ?_) j
    · rw [Fin.append_left]
    · rw [Fin.append_right]
  have hA : ∀ i x, ip (A i) x = ip (fun j => a i j) (fun j => x (Fin.castAdd m j)) := by
    intro i x
    conv_lhs => rw [split x]
    rw [hAdef, ip_append]; simp [ip]
  have hB : ∀ i x, ip (B i) x = ip (fun j => b i j) (fun j => x (Fin.natAdd m j)) := by
    intro i x
    conv_lhs => rw [split x]
    rw [hBdef, ip_append]; simp [ip]
  rw [orthant_transfer a _ (by fun_prop) (map_castAdd m m) A hA t,
    orthant_transfer b _ (by fun_prop) (map_natAdd m m) B hB t]
  refine pi_slepian A B (fun i l => ?_) (fun i => ?_) (fun i l hil => ?_) t
  · rw [hAdef, hBdef, ip_append]; simp [ip]
  · rw [hAdef, hBdef, ip_append, ip_append, ← inner_eq_ip, ← inner_eq_ip,
      real_inner_self_eq_norm_sq, real_inner_self_eq_norm_sq, hn i]
    simp [ip]
  · rw [hAdef, hBdef, ip_append, ip_append, ← inner_eq_ip, ← inner_eq_ip]
    simpa [ip] using ho i l hil

end T2
end CubicPathwidth

#print axioms CubicPathwidth.T2.slepian