AutoSamplingTheory.ExampleCases.SmoothedPicardHMC.ReferenceCarryingKernel

No declaration docstring.

private theorem state_iterates {g : E → E} (hg : Measurable g)
    {b step : X → ℝ} {u initial : X → E}
    (hb : Measurable b) (hs : Measurable step) (hu : Measurable u)
    (hi : Measurable initial) (n : ℕ) :
    Measurable (fun s => (fun x => x-step s • (g x+b s • (x-u s)))^[n] (initial s)) := by
  induction n with
  | zero => simpa only [Function.iterate_zero,id_eq] using hi
  | succ n ih =>
    simp only [Function.iterate_succ_apply']
    exact ih.sub (hs.smul ((hg.comp ih).add (hb.smul (ih.sub hu))))

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:44published source at 0e31a3cda412

No declaration docstring.

private def firstIndex (q : ℕ → X → ℝ) (b : X → ℝ) (s : X) : ℕ := by
  classical
  exact if h : ∃ n, q n s ≤ b s then Nat.find h else 0

omit [NormedAddCommGroup E] [InnerProductSpace ℝ E] [FiniteDimensional ℝ E] [BorelSpace E] in

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:55published source at 0e31a3cda412

No declaration docstring.

private theorem stopped_family (x : ℕ → X → E) (q : ℕ → X → ℝ) (b : X → ℝ)
    (hx : ∀ n, Measurable (x n))
    (hq : ∀ n, MeasurableSet {s | q n s ≤ b s})
    (ht : ∀ s, ∃ n, q n s ≤ b s) :
    let N := firstIndex q b
    Measurable N ∧ Measurable (fun s => x (N s) s) ∧
      ∀ s, q (N s) s ≤ b s ∧ ∀ j < N s, b s < q j s := by
  classical
  let N := firstIndex q b
  have hn (s : X) : N s=Nat.find (ht s) := by
    simp only [N,firstIndex,dif_pos (ht s)]
  have hm : Measurable N := by
    rw [show N=(fun s => Nat.find (ht s)) from funext hn]
    exact measurable_find ht hq
  have ho : Measurable (fun s => x (N s) s) := by
    simp_rw [hn]
    exact Measurable.find hx hq ht
  refine ⟨hm,ho,fun s => ?_⟩
  change q (N s) s ≤ b s ∧ ∀ j < N s, b s < q j s
  rw [hn]
  exact ⟨Nat.find_spec (ht s),fun j hj => lt_of_not_ge (Nat.find_min (ht s) hj)⟩

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:60published source at 0e31a3cda412

No declaration docstring.

private theorem variable_reference {V : E → ℝ} {α β : ℝ≥0}
    (hV : ContDiff ℝ 2 V)
    (hH : ∀ x v : E, (α:ℝ)*‖v‖^2 ≤ fderiv ℝ (fderiv ℝ V) x v v ∧
      fderiv ℝ (fderiv ℝ V) x v v ≤ (β:ℝ)*‖v‖^2)
    (hαβ : α ≤ β) (hd : 0 < (Module.finrank ℝ E : ℝ))
    {b : X → ℝ} {u initial : X → E} (hb : Measurable b)
    (hb0 : ∀ s, 0 < b s) (hu : Measurable u) (hi : Measurable initial) :
    let d : ℝ := Module.finrank ℝ E
    let F := fun s x => V x+b s/2*‖x-u s‖^2
    let T := fun s x => x-((β:ℝ)+b s)⁻¹ • gradient (F s) x
    let q := fun n s => ‖gradient (F s) ((T s)^[n] (initial s))‖^2
    let N := firstIndex q (fun s => ((α:ℝ)+b s)*d)
    Measurable N ∧ Measurable (fun s => (T s)^[N s] (initial s)) ∧
      (∀ s x, gradient (F s) x=gradient V x+b s • (x-u s)) ∧
      ∀ s, q (N s) s ≤ ((α:ℝ)+b s)*d ∧
        ∀ j < N s, ((α:ℝ)+b s)*d < q j s := by
  let d : ℝ := Module.finrank ℝ E
  let F := fun s x => V x+b s/2*‖x-u s‖^2
  let T := fun s x => x-((β:ℝ)+b s)⁻¹ • gradient (F s) x
  let q := fun n s => ‖gradient (F s) ((T s)^[n] (initial s))‖^2
  have hgrad (s : X) (x : E) : gradient (F s) x=gradient V x+b s • (x-u s) := by
    have h := TerminalReferenceGradientDescent.terminal_reference_gradient_descent
      hV hH hαβ (inv_pos.mpr (hb0 s)) hd rfl (u s)
    simpa only [inv_inv] using h.2.2.2.1 x
  have hgV : Measurable (gradient V) :=
    (TechnicalLemmas.Analysis.Calculus.Gradient.continuous_gradient_of_contDiff_one
      (hV.of_le (by norm_num))).measurable
  have hTeq : T=(fun s x => x-((β:ℝ)+b s)⁻¹ • (gradient V x+b s • (x-u s))) := by
    funext s x
    dsimp [T]
    rw [hgrad]
  have hm (n : ℕ) : Measurable (fun s => (T s)^[n] (initial s)) := by
    rw [hTeq]
    exact state_iterates hgV hb ((measurable_const.add hb).inv) hu hi n
  have hq (n : ℕ) : Measurable (q n) := by
    change Measurable (fun s => ‖gradient (F s) ((T s)^[n] (initial s))‖^2)
    simp_rw [hgrad]
    exact ((hgV.comp (hm n)).add (hb.smul ((hm n).sub hu))).norm.pow_const 2
  have hevent (n : ℕ) : MeasurableSet {s | q n s ≤ ((α:ℝ)+b s)*d} :=
    measurableSet_le (hq n) ((measurable_const.add hb).mul_const d)
  have ht (s : X) : ∃ n, q n s ≤ ((α:ℝ)+b s)*d := by
    have h := TerminalReferenceGradientDescent.terminal_reference_gradient_descent
      hV hH hαβ (inv_pos.mpr (hb0 s)) hd rfl (u s)
    simp only [inv_inv] at h
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:82published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.

No declaration docstring.

private theorem quadratic_source_alignment (V : E → ℝ) {r t : ℝ}
    (hr : 0 ≤ r) (ht : 0 < t) (u y : E) :
    let w := (r+t)⁻¹ • (r • u+t • y)
    let F := fun x => V x+r/2*‖x-u‖^2+t/2*‖x-y‖^2
    let G := fun x => V x+(r+t)/2*‖x-w‖^2
    (∀ x, F x=G x+r*t/(2*(r+t))*‖u-y‖^2) ∧
      (∀ x, gradient F x=gradient G x) ∧
      ∀ (h : ℝ) (n : ℕ) (x : E),
        (fun z => z-h • gradient F z)^[n] x =
          (fun z => z-h • gradient G z)^[n] x := by
  let w := (r+t)⁻¹ • (r • u+t • y)
  let F := fun x => V x+r/2*‖x-u‖^2+t/2*‖x-y‖^2
  let G := fun x => V x+(r+t)/2*‖x-w‖^2
  have hrt : 0 < r+t := add_pos_of_nonneg_of_pos hr ht
  have he (x : E) : F x=G x+r*t/(2*(r+t))*‖u-y‖^2 := by
    dsimp [F,G,w]
    simp only [← real_inner_self_eq_norm_sq,inner_sub_left,inner_sub_right,
      inner_add_left,inner_add_right,real_inner_smul_left,real_inner_smul_right]
    rw [real_inner_comm y u,real_inner_comm u x,real_inner_comm y x]
    field_simp
    ring
  have hg (x : E) : gradient F x=gradient G x := by
    rw [show F=(fun x => G x+r*t/(2*(r+t))*‖u-y‖^2) from funext he]
    simp only [gradient,fderiv_add_const]
  refine ⟨he,hg,fun h n x => ?_⟩
  have hh : (fun z => z-h • gradient F z)=(fun z => z-h • gradient G z) := by
    funext z
    rw [hg]
  rw [hh]

omit [MeasurableSpace E] [BorelSpace E] [MeasurableSpace X] in

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:137published source at 0e31a3cda412

No declaration docstring.

private theorem scaled_reference_bound {F : E → ℝ} (hF : Differentiable ℝ F)
    {c a d : ℝ} (hc : 0 < c) {x : E} (hx : ‖gradient F x‖^2 ≤ a*d) :
    ‖gradient (fun z => F ((Real.sqrt c)⁻¹ • z)) (Real.sqrt c • x)‖^2 ≤ d/(c/a) := by
  let s := (Real.sqrt c)⁻¹
  let L : E →L[ℝ] E := s • ContinuousLinearMap.id ℝ E
  have hL (z : E) : HasFDerivAt (fun w : E => s • w) L z :=
    (hasFDerivAt_id z).const_smul s
  have hs : Real.sqrt c ≠ 0 := (Real.sqrt_pos.mpr hc).ne'
  have hscale : s • (Real.sqrt c • x)=x := by
    simp [s,smul_smul,inv_mul_cancel₀ hs]
  have hg : gradient (fun z => F (s • z)) (Real.sqrt c • x)=s • gradient F x := by
    apply HasGradientAt.gradient
    rw [hasGradientAt_iff_hasFDerivAt]
    have he : (InnerProductSpace.toDual ℝ E (gradient F x)).comp L =
        InnerProductSpace.toDual ℝ E (s • gradient F x) := by
      ext v
      simp [L]
    rw [← he]
    have h : HasFDerivAt (fun z => F (s • z))
        ((InnerProductSpace.toDual ℝ E (gradient F (s • (Real.sqrt c • x)))).comp L)
        (Real.sqrt c • x) :=
      (hF (s • (Real.sqrt c • x))).hasGradientAt.hasFDerivAt.comp _ (hL _)
    simpa only [hscale] using h
  rw [hg,norm_smul,Real.norm_eq_abs,mul_pow,sq_abs]
  have he : s^2=c⁻¹ := by dsimp [s]; rw [inv_pow,Real.sq_sqrt hc.le]
  rw [he]
  calc
    c⁻¹*‖gradient F x‖^2 ≤ c⁻¹*(a*d) := mul_le_mul_of_nonneg_left hx (inv_pos.mpr hc).le
    _ = d/(c/a) := by field_simp

local notation "RefState" E:arg => ℝ≥0 × E × E × ℕ × (ℕ → E)
local notation "OldState" E:arg => ℝ≥0 × E × ℕ × (ℕ → E)

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:168published source at 0e31a3cda412

No declaration docstring.

theorem reference_carrying_kernel {V : E → ℝ} {α β : ℝ≥0}
    (hα : 0 < α) (hαβ : α ≤ β) (hV : ContDiff ℝ 2 V)
    (hH : ∀ x w : E, (α:ℝ)*‖w‖^2 ≤ fderiv ℝ (fderiv ℝ V) x w w ∧
      fderiv ℝ (fderiv ℝ V) x w w ≤ (β:ℝ)*‖w‖^2)
    (hd : 0 < (Module.finrank ℝ E : ℝ))
    (η τ : RefState E → ℝ) (hη : Measurable η) (hτ : Measurable τ)
    (hη0 : ∀ s, 0 < η s) (hτ0 : ∀ s, 0 < τ s)
    (M : Kernel (RefState E) E) [IsMarkovKernel M] (threshold : ℝ≥0) :
    let d : ℝ := Module.finrank ℝ E
    let v := fun s : RefState E => (η s+τ s)/((β:ℝ)+s.1)
    let bp := fun s : RefState E => s.1+Real.toNNReal (v s)⁻¹
    let obs := fun p : RefState E × (E × E) => p.2.1+Real.sqrt (τ p.1/((β:ℝ)+p.1.1)) • p.2.2
    let center := fun p : RefState E × (E × E) => (bp p.1:ℝ)⁻¹ •
      ((p.1.1:ℝ) • p.1.2.1+(v p.1)⁻¹ • obs p)
    let F := fun p x => V x+(bp p.1:ℝ)/2*‖x-center p‖^2
    let W := fun p x => V x+(p.1.1:ℝ)/2*‖x-p.1.2.1‖^2+(v p.1)⁻¹/2*‖x-obs p‖^2
    let T := fun p x => x-((β:ℝ)+bp p.1)⁻¹ • gradient (F p) x
    let q := fun n p => ‖gradient (F p) ((T p)^[n] p.2.1)‖^2
    let N := firstIndex q (fun p => ((α:ℝ)+bp p.1)*d)
    let out := fun p => (T p)^[N p] p.2.1
    let update := fun p : RefState E × (E × E) =>
      (bp p.1,center p,out p,p.1.2.2.2.1+1,fun n => Nat.casesOn n (obs p) p.1.2.2.2.2)
    let project := fun s : RefState E => (s.1,s.2.1,s.2.2.2.1,s.2.2.2.2)
    let oldUpdate := fun p : RefState E × (E × E) =>
      (bp p.1,center p,p.1.2.2.2.1+1,fun n => Nat.casesOn n (obs p) p.1.2.2.2.2)
    (∀ s, 0 < v s ∧ (bp s:ℝ)=(s.1:ℝ)+(v s)⁻¹) ∧
    Measurable N ∧ Measurable out ∧ Measurable update ∧
    (∀ p, q (N p) p ≤ ((α:ℝ)+bp p.1)*d ∧
      (∀ j < N p, ((α:ℝ)+bp p.1)*d < q j p) ∧
      (∀ x, gradient (F p) x=gradient V x+(bp p.1:ℝ) • (x-center p)) ∧
      (∀ x, W p x=F p x+
        (p.1.1:ℝ)*(v p.1)⁻¹/(2*((p.1.1:ℝ)+(v p.1)⁻¹))*‖p.1.2.1-obs p‖^2) ∧
      (∀ x, gradient (W p) x=gradient (F p) x) ∧
      (∀ n, (fun x => x-((β:ℝ)+bp p.1)⁻¹ • gradient (W p) x)^[n] p.2.1 =
        (T p)^[n] p.2.1) ∧
      ‖gradient (fun x => F p ((Real.sqrt ((β:ℝ)+bp p.1))⁻¹ • x))
          (Real.sqrt ((β:ℝ)+bp p.1) • out p)‖^2 ≤
        d/(((β:ℝ)+bp p.1)/((α:ℝ)+bp p.1))) ∧
    ∃ P : Kernel (RefState E) (RefState E), IsMarkovKernel P ∧
      (∀ s, P s=if threshold ≤ s.1 then Measure.dirac s else
        ((M s).prod (stdGaussian E)).map (fun z => update (s,z))) ∧
      (∀ s, (P s).map project=if threshold ≤ s.1 then Measure.dirac (project s) else
        ((M s).prod (stdGaussian E)).map (fun z => oldUpdate (s,z))) := by
  classical
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/ReferenceCarryingKernel.lean:201published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.