Mathlib.Probability.Gaussian­Mean­Embedding

lemma gaussianFeature_memℓp (r : ℝ) : Memℓp (gaussianFeatureCoefficient r) 2 := by apply memℓp_gen norm_num simp only [gaussianFeatureCoefficient, gaussianFeatureNormalization] have hs : Summable (fun m : ℕ => (2 * r ^ 2) ^ m / m.factorial) := Real.summable_pow_div_factorial _ apply (hs.mul_left (exp (-r ^ 2) ^ 2)).congr intro m have hnonneg : 0 ≤ (2 : ℝ) ^ m / m.factorial := by positivity calc _ = exp (-r ^ 2) ^ 2 * ((2 : ℝ) ^ m / m.factorial) * r ^ (2 * m) := by ring _ = (sqrt ((2 : ℝ) ^ m / m.factorial)) ^ 2 * exp (-r ^ 2) ^ 2 * r ^ (2 * m) := by rw [Real.sq_sqrt hnonneg] ac_rfl _ = _ := by ring

lemma gaussianFeature_inner (a b : ℝ) : ⟪gaussianFeature a, gaussianFeature b⟫_ℝ = gaussianKernel a b := by rw [lp.inner_eq_tsum] simp_rw [RCLike.inner_apply, starRingEnd_apply, star_trivial] change (∑' m : ℕ, gaussianFeatureCoefficient b m * gaussianFeatureCoefficient a m) = _ have hterm : (fun m : ℕ => gaussianFeatureCoefficient b m * gaussianFeatureCoefficient a m) = fun m => (exp (-a ^ 2) * exp (-b ^ 2)) * ((2 * a * b) ^ m / m.factorial) := by funext m unfold gaussianFeatureCoefficient gaussianFeatureNormalization have hnonneg : 0 ≤ (2 : ℝ) ^ m / m.factorial := by positivity have hsqrt := Real.mul_self_sqrt hnonneg calc _ = (sqrt ((2 : ℝ) ^ m / m.factorial) * sqrt ((2 : ℝ) ^ m / m.factorial)) * (exp (-a ^ 2) * exp (-b ^ 2)) * (a ^ m * b ^ m) := by ring _ = _ := by rw [hsqrt]; ring rw [hterm, tsum_mul_left, (NormedSpace.expSeries_div_hasSum_exp (2 * a * b)).tsum_eq, ← Real.exp_eq_exp_ℝ, ← Real.exp_add, ← Real.exp_add] unfold gaussianKernel congr 1 ring

lemma gaussianFeature_norm (r : ℝ) : ‖gaussianFeature r‖ = 1 := by have hsq : ‖gaussianFeature r‖ ^ 2 = 1 := by rw [← real_inner_self_eq_norm_sq, gaussianFeature_inner] simp [gaussianKernel] nlinarith [norm_nonneg (gaussianFeature r)]

lemma UnitNormGaussianFeatureMap.continuous {H : Type*} [NormedAddCommGroup H] [InnerProductSpace ℝ H] (U : UnitNormGaussianFeatureMap H) : Continuous U.feature := by rw [continuous_iff_continuousAt] intro a change Tendsto U.feature (nhds a) (nhds (U.feature a)) rw [Metric.tendsto_nhds] intro ε hε have hcont : Tendsto (fun b : ℝ => ‖U.feature b - U.feature a‖ ^ 2) (nhds a) (nhds 0) := by have h : ContinuousAt (fun b : ℝ => 2 - 2 * gaussianKernel b a) a := by unfold gaussianKernel fun_prop have heq : (fun b : ℝ => ‖U.feature b - U.feature a‖ ^ 2) = fun b => 2 - 2 * gaussianKernel b a := by funext b exact U.dist_sq b a rw [heq] change Tendsto (fun b : ℝ => 2 - 2 * gaussianKernel b a) (nhds a) (nhds ((fun b : ℝ => 2 - 2 * gaussianKernel b a) a)) at h simpa [gaussianKernel] using h have hev : ∀ᶠ b in nhds a, ‖U.feature b - U.feature a‖ ^ 2 < ε ^ 2 := (tendsto_order.1 hcont).2 _ (sq_pos_of_pos hε) filter_upwards [hev] with b hb rw [dist_eq_norm] have hnonneg := norm_nonneg (U.feature b - U.feature a) nlinarith

lemma UnitNormGaussianFeatureMap.integrable {H : Type*} [NormedAddCommGroup H] [InnerProductSpace ℝ H] [CompleteSpace H] (U : UnitNormGaussianFeatureMap H) (μ : Measure ℝ) [IsFiniteMeasure μ] : Integrable U.feature μ := by exact Integrable.of_bound U.continuous.stronglyMeasurable.aestronglyMeasurable 1 (Filter.Eventually.of_forall fun r => by rw [U.norm_eq_one])

lemma gaussianCoordinate_meanEmbedding (μ : Measure ℝ) [IsFiniteMeasure μ] (m : ℕ) : gaussianCoordinate m (meanEmbedding gaussianFeatureMap μ) = ∫ r, gaussianFeatureCoefficient r m ∂μ := by unfold meanEmbedding rw [← (gaussianCoordinate m).integral_comp_comm (gaussianFeatureMap.integrable μ)] rfl

lemma integrable_gaussianWeightedMonomial (μ : Measure ℝ) [IsFiniteMeasure μ] (m : ℕ) : Integrable (gaussianWeightedMonomial m) μ := by let f : ZeroAtInftyContinuousMap ℝ ℝ := { toFun := gaussianWeightedMonomial m continuous_toFun := by unfold gaussianWeightedMonomial gaussianWeight fun_prop zero_at_infty' := by rw [tendsto_zero_iff_norm_tendsto_zero] have h := tendsto_rpow_abs_mul_exp_neg_mul_sq_cocompact (a := (1 : ℝ)) (by positivity) (m : ℝ) have h' : Tendsto (fun x : ℝ => |x| ^ m * exp (-1 * x ^ 2)) (cocompact ℝ) (nhds 0) := by simpa only [Real.rpow_natCast] using h convert h' using 1 simp [gaussianWeightedMonomial, gaussianWeight, Real.norm_eq_abs, abs_of_pos (Real.exp_pos _), mul_comm] } exact f.toBCF.integrable μ

lemma gaussianWeightedPolynomial_integral_eq (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (hmom : ∀ m : ℕ, ∫ r, gaussianWeightedMonomial m r ∂μ = ∫ r, gaussianWeightedMonomial m r ∂ν) (p : ℝ[X]) : ∫ r, gaussianWeight r * p.eval r ∂μ = ∫ r, gaussianWeight r * p.eval r ∂ν := by induction p using Polynomial.induction_on' with | add p q hp hq => have hpμ := integrable_gaussianWeightedPolynomial μ p have hqμ := integrable_gaussianWeightedPolynomial μ q have hpν := integrable_gaussianWeightedPolynomial ν p have hqν := integrable_gaussianWeightedPolynomial ν q simp only [Polynomial.eval_add, mul_add] rw [integral_add hpμ hqμ, integral_add hpν hqν, hp, hq] | monomial n a => have hfun : (fun r : ℝ => gaussianWeight r * ((Polynomial.monomial n a).eval r)) = fun r => a * gaussianWeightedMonomial n r := by funext r simp only [Polynomial.eval_monomial, gaussianWeightedMonomial] ring rw [hfun, integral_const_mul, integral_const_mul, hmom n]

theorem continuousIntegral_eq_of_gaussianWeightedMoments_eq (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (B : ℝ) (hB : 0 ≤ B) (hμ : μ (Icc (-B) B)ᶜ = 0) (hν : ν (Icc (-B) B)ᶜ = 0) (hmom : ∀ m : ℕ, ∫ r, gaussianWeightedMonomial m r ∂μ = ∫ r, gaussianWeightedMonomial m r ∂ν) (f : ℝ → ℝ) (hf : Continuous f) : ∫ r, f r ∂μ = ∫ r, f r ∂ν := by have _hB : -B ≤ B := by linarith apply eq_of_forall_dist_le intro ε hε let M : ℝ := μ.real univ + ν.real univ + 1 have hMpos : 0 < M := by dsimp [M] have hμnonneg : 0 ≤ μ.real univ := ENNReal.toReal_nonneg have hνnonneg : 0 ≤ ν.real univ := ENNReal.toReal_nonneg linarith let δ : ℝ := ε / M have hδ : 0 < δ := div_pos hε hMpos have hquot : Continuous (fun r : ℝ => f r / gaussianWeight r) := by apply hf.div · unfold gaussianWeight fun_prop · intro r unfold gaussianWeight exact (Real.exp_pos _).ne' obtain ⟨p, hp⟩ := exists_polynomial_near_of_continuousOn (-B) B (fun r : ℝ => f r / gaussianWeight r) hquot.continuousOn δ hδ let g : ℝ → ℝ := fun r => gaussianWeight r * p.eval r have hg : Continuous g := by dsimp [g] unfold gaussianWeight fun_prop have hμmem : ∀ᵐ r ∂μ, r ∈ Icc (-B) B := by rw [ae_iff] exact hμ have hνmem : ∀ᵐ r ∂ν, r ∈ Icc (-B) B := by rw [ae_iff] exact hν have hμres : μ.restrict (Icc (-B) B) = μ := Measure.restrict_eq_self_of_ae_mem hμmem have hνres : ν.restrict (Icc (-B) B) = ν := Measure.restrict_eq_self_of_ae_mem hνmem have hfμ : Integrable f μ := by rw [← hμres] exact hf.integrableOn_Icc have hfν : Integrable f ν := by rw [← hνres] exact hf.integrableOn_Icc have hgμ : Integrable g μ := by rw [← hμres] exact hg.integrableOn_Icc have hgν : Integrable g ν := by rw [← hνres] exact hg.integrableOn_Icc have herror : ∀ r ∈ Icc (-B) B, ‖f r - g r‖ ≤ δ := by intro r hr have hwpos : 0 < gaussianWeight r := by unfold gaussianWeight exact Real.exp_pos _ have hwle : gaussianWeight r ≤ 1 := by unfold gaussianWeight rw [Real.exp_le_one_iff] exact neg_nonpos.mpr (sq_nonneg r) have happ := hp r hr rw [Real.norm_eq_abs] dsimp [g] have hid : f r - gaussianWeight r * p.eval r = -gaussianWeight r * (p.eval r - f r / gaussianWeight r) := by field_simp ring rw [hid, abs_mul, abs_neg, abs_of_pos hwpos] have habs := abs_nonneg (p.eval r - f r / gaussianWeight r) nlinarith have herrμ : ‖∫ r, (f r - g r) ∂μ‖ ≤ δ * μ.real univ := by apply norm_integral_le_of_norm_le_const filter_upwards [hμmem] with r hr exact herror r hr have herrν : ‖∫ r, (f r - g r) ∂ν‖ ≤ δ * ν.real univ := by apply norm_integral_le_of_norm_le_const filter_upwards [hνmem] with r hr exact herror r hr rw [integral_sub hfμ hgμ] at herrμ rw [integral_sub hfν hgν] at herrν have hpoly : (∫ r, g r ∂μ) = ∫ r, g r ∂ν := by dsimp [g] exact gaussianWeightedPolynomial_integral_eq μ ν hmom p have hmassμ : 0 ≤ μ.real univ := ENNReal.toReal_nonneg have hmassν : 0 ≤ ν.real univ := ENNReal.toReal_nonneg calc dist (∫ r, f r ∂μ) (∫ r, f r ∂ν) = ‖((∫ r, f r ∂μ) - ∫ r, g r ∂μ) + ((∫ r, g r ∂ν) - ∫ r, f r ∂ν)‖ := by rw [dist_eq_norm, hpoly] congr 1 ring _ ≤ ‖(∫ r, f r ∂μ) - ∫ r, g r ∂μ‖ + ‖(∫ r, g r ∂ν) - ∫ r, f r ∂ν‖ := norm_add_le _ _ _ = ‖(∫ r, f r ∂μ) - ∫ r, g r ∂μ‖ + ‖(∫ r, f r ∂ν) - ∫ r, g r ∂ν‖ := by congr 1 exact norm_sub_rev _ _ _ ≤ δ * μ.real univ + δ * ν.real univ := add_le_add herrμ herrν _ = δ * (μ.real univ + ν.real univ) := by ring _ ≤ δ * (μ.real univ + ν.real univ + 1) := by nlinarith _ = ε := by dsimp [δ, M] field_simp

theorem measure_eq_of_continuousIntegral_eq (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (h : ∀ (f : ℝ → ℝ), Continuous f → ∫ r, f r ∂μ = ∫ r, f r ∂ν) : μ = ν := by apply Measure.ext_of_integral_eq_on_compactlySupported intro f exact h f f.continuous

theorem measure_eq_of_gaussianWeightedMoments_eq (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (B : ℝ) (hB : 0 ≤ B) (hμ : μ (Icc (-B) B)ᶜ = 0) (hν : ν (Icc (-B) B)ᶜ = 0) (hmom : ∀ m : ℕ, ∫ r, gaussianWeightedMonomial m r ∂μ = ∫ r, gaussianWeightedMonomial m r ∂ν) : μ = ν := by apply measure_eq_of_continuousIntegral_eq μ ν intro f hf exact continuousIntegral_eq_of_gaussianWeightedMoments_eq μ ν B hB hμ hν hmom f hf

lemma gaussianWeightedMoments_eq_of_meanEmbedding_eq (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (hemb : meanEmbedding gaussianFeatureMap μ = meanEmbedding gaussianFeatureMap ν) (m : ℕ) : ∫ r, gaussianWeightedMonomial m r ∂μ = ∫ r, gaussianWeightedMonomial m r ∂ν := by have hcoord := congrArg (gaussianCoordinate m) hemb rw [gaussianCoordinate_meanEmbedding, gaussianCoordinate_meanEmbedding] at hcoord have hscaled : gaussianFeatureNormalization m * ∫ r, gaussianWeightedMonomial m r ∂μ = gaussianFeatureNormalization m * ∫ r, gaussianWeightedMonomial m r ∂ν := by simpa only [gaussianFeatureCoefficient, gaussianWeightedMonomial, gaussianWeight, mul_assoc, integral_const_mul] using hcoord exact mul_left_cancel₀ (ne_of_gt (gaussianFeatureNormalization_pos m)) hscaled

theorem secondMoment_eq_of_meanEmbedding_eq_of_boundedSupport (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (B : ℝ) (hB : 0 ≤ B) (hμ : μ (Icc (-B) B)ᶜ = 0) (hν : ν (Icc (-B) B)ᶜ = 0) (hemb : meanEmbedding gaussianFeatureMap μ = meanEmbedding gaussianFeatureMap ν) : (∫ r, r ^ 2 ∂μ) = ∫ r, r ^ 2 ∂ν := by rw [measure_eq_of_meanEmbedding_eq_of_boundedSupport μ ν B hB hμ hν hemb]

theorem norm_meanEmbedding_sub_pos_of_secondMoment_ne_of_boundedSupport (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (B : ℝ) (hB : 0 ≤ B) (hμ : μ (Icc (-B) B)ᶜ = 0) (hν : ν (Icc (-B) B)ᶜ = 0) (hmom : (∫ r, r ^ 2 ∂μ) ≠ ∫ r, r ^ 2 ∂ν) : 0 < ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ := by rw [norm_pos_iff] intro hzero apply hmom apply secondMoment_eq_of_meanEmbedding_eq_of_boundedSupport μ ν B hB hμ hν exact sub_eq_zero.mp hzero

theorem abs_gaussianWeightedMoment_sub_le (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (m : ℕ) : |(∫ r, gaussianWeightedMonomial m r ∂μ) - ∫ r, gaussianWeightedMonomial m r ∂ν| ≤ ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ / gaussianFeatureNormalization m := by have hcoord : gaussianCoordinate m (meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν) = gaussianFeatureNormalization m * ((∫ r, gaussianWeightedMonomial m r ∂μ) - ∫ r, gaussianWeightedMonomial m r ∂ν) := by rw [map_sub, gaussianCoordinate_meanEmbedding, gaussianCoordinate_meanEmbedding] simp only [gaussianFeatureCoefficient, gaussianWeightedMonomial, gaussianWeight, mul_assoc, integral_const_mul] ring have hbound : |gaussianFeatureNormalization m * ((∫ r, gaussianWeightedMonomial m r ∂μ) - ∫ r, gaussianWeightedMonomial m r ∂ν)| ≤ ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ := by rw [← hcoord, ← Real.norm_eq_abs] exact lp.norm_apply_le_norm (by norm_num) (meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν) m apply (le_div_iff₀ (gaussianFeatureNormalization_pos m)).2 simpa only [abs_mul, abs_of_pos (gaussianFeatureNormalization_pos m), mul_comm] using hbound

theorem abs_secondMomentWeightedSum_sub_le (μ ν : Measure ℝ) [IsFiniteMeasure μ] [IsFiniteMeasure ν] (N : ℕ) : |∑ k ∈ range (N + 1), (1 / (k.factorial : ℝ)) * ((∫ r, gaussianWeightedMonomial (2 * k + 2) r ∂μ) - ∫ r, gaussianWeightedMonomial (2 * k + 2) r ∂ν)| ≤ secondMomentRecoveryCoefficient N * ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ := by simpa [secondMomentRecoveryCoefficient] using (abs_sum_gaussianWeightedMoment_sub_le μ ν (range (N + 1)) (fun k => 2 * k + 2) (fun k => 1 / (k.factorial : ℝ)))

theorem secondMomentRecoveryCoefficient_term_le (k : ℕ) : (1 / (k.factorial : ℝ)) / gaussianFeatureNormalization (2 * k + 2) ≤ (k + 1 : ℕ) := by have hchoose : (2 * (k + 1)).choose (k + 1) ≤ 2 ^ (2 * (k + 1)) := Nat.choose_le_two_pow _ _ have hfac_nat : (2 * (k + 1)).factorial ≤ 2 ^ (2 * (k + 1)) * (k + 1).factorial * (k + 1).factorial := by calc (2 * (k + 1)).factorial = (2 * (k + 1)).choose (k + 1) * (k + 1).factorial * (k + 1).factorial := by rw [show 2 * (k + 1) = (k + 1) + (k + 1) by omega] exact (Nat.add_choose_mul_factorial_mul_factorial (k + 1) (k + 1)).symm _ ≤ 2 ^ (2 * (k + 1)) * (k + 1).factorial * (k + 1).factorial := by gcongr have hfac : ((2 * (k + 1)).factorial : ℝ) ≤ (2 : ℝ) ^ (2 * (k + 1)) * ((k + 1).factorial : ℝ) * ((k + 1).factorial : ℝ) := by exact_mod_cast hfac_nat have hnorm_pos : 0 < gaussianFeatureNormalization (2 * k + 2) := gaussianFeatureNormalization_pos _ have hkfac_pos : 0 < (k.factorial : ℝ) := by positivity have hnorm_sq : gaussianFeatureNormalization (2 * k + 2) ^ 2 = (2 : ℝ) ^ (2 * k + 2) / ((2 * k + 2).factorial : ℝ) := by unfold gaussianFeatureNormalization rw [Real.sq_sqrt] positivity have hsquare : ((1 / (k.factorial : ℝ)) / gaussianFeatureNormalization (2 * k + 2)) ^ 2 ≤ ((k + 1 : ℕ) : ℝ) ^ 2 := by rw [Nat.factorial_succ] at hfac rw [div_pow, div_pow, one_pow, hnorm_sq] field_simp norm_num [Nat.cast_add, Nat.cast_one] at hfac ⊢ convert hfac using 1 <;> ring_nf have hlhs_nonneg : 0 ≤ (1 / (k.factorial : ℝ)) / gaussianFeatureNormalization (2 * k + 2) := by positivity simp only [Nat.cast_add, Nat.cast_one] at hsquare ⊢ nlinarith

theorem secondMomentRecoveryCoefficient_oneHundred_le : secondMomentRecoveryCoefficient 100 ≤ 5151 := by unfold secondMomentRecoveryCoefficient calc ∑ k ∈ range (100 + 1), (1 / (k.factorial : ℝ)) / gaussianFeatureNormalization (2 * k + 2) ≤ ∑ k ∈ range (100 + 1), ((k + 1 : ℕ) : ℝ) := by apply Finset.sum_le_sum intro k hk exact secondMomentRecoveryCoefficient_term_le k _ = 5151 := by norm_num

theorem secondMomentTaylorPolynomial_eval (N : ℕ) (r : ℝ) : (secondMomentTaylorPolynomial N).eval r = ∑ k ∈ range (N + 1), r ^ (2 * k + 2) / (k.factorial : ℝ) := by rw [secondMomentTaylorPolynomial, Polynomial.eval_finsetSum] apply Finset.sum_congr rfl intro k _ rw [Polynomial.eval_monomial] ring

theorem secondMomentTaylorPolynomial_natDegree_oneHundred : (secondMomentTaylorPolynomial 100).natDegree = 202 := by apply Polynomial.natDegree_eq_of_le_of_coeff_ne_zero · apply Polynomial.natDegree_sum_le_of_forall_le intro k hk rw [Polynomial.natDegree_monomial] split_ifs · omega · have hk' : k < 101 := Finset.mem_range.mp hk omega · set_option maxRecDepth 10000 in rw [secondMomentTaylorPolynomial] rw [Polynomial.finsetSum_coeff] simp only [Polynomial.coeff_monomial] rw [Finset.sum_eq_single 100] · norm_num · intro b hb hne split_ifs with h · omega · rfl · simp

theorem exp_series_tail_le_geometric (N : ℕ) (x : ℝ) (hx : 0 ≤ x) (hden : x < N + 2) : exp x - ∑ k ∈ range (N + 1), x ^ k / (k.factorial : ℝ) ≤ (x ^ (N + 1) / ((N + 1).factorial : ℝ)) * ((N + 2 : ℝ) / ((N + 2 : ℝ) - x)) := by let f : ℕ → ℝ := fun k => x ^ k / (k.factorial : ℝ) have hf : Summable f := by simpa [f] using NormedSpace.expSeries_div_summable x have hsplit := hf.sum_add_tsum_nat_add (N + 1) have hexp : exp x = ∑' k, f k := by rw [Real.exp_eq_exp_ℝ] simpa [f] using (NormedSpace.expSeries_div_hasSum_exp x).tsum_eq.symm have htail : exp x - ∑ k ∈ range (N + 1), x ^ k / (k.factorial : ℝ) = ∑' j : ℕ, f (j + (N + 1)) := by rw [hexp] change (∑' k, f k) - ∑ k ∈ range (N + 1), f k = _ linarith rw [htail] let a : ℝ := x ^ (N + 1) / ((N + 1).factorial : ℝ) let q : ℝ := x / (N + 2 : ℝ) have hq : |q| < 1 := by rw [abs_of_nonneg (div_nonneg hx (by positivity))] exact (div_lt_one (by positivity)).2 hden have hfshift : Summable (fun j : ℕ => f (j + (N + 1))) := by exact hf.comp_injective (fun _ _ h => Nat.add_right_cancel h) have hg : Summable (fun j : ℕ => a * q ^ j) := (summable_geometric_of_norm_lt_one (by simpa [Real.norm_eq_abs] using hq)).mul_left a calc (∑' j : ℕ, f (j + (N + 1))) ≤ ∑' j : ℕ, a * q ^ j := by apply Summable.tsum_le_tsum _ hfshift hg intro j dsimp [f, a, q] have hfac : ((N + 1).factorial : ℝ) * (N + 2 : ℝ) ^ j ≤ ((N + 1 + j).factorial : ℝ) := by exact_mod_cast factorial_mul_pow_le_factorial_add N j have hpos : 0 < ((N + 1).factorial : ℝ) * (N + 2 : ℝ) ^ j := by positivity have hnum : 0 ≤ x ^ (N + 1 + j) := by positivity rw [show j + (N + 1) = N + 1 + j by omega] calc x ^ (N + 1 + j) / ((N + 1 + j).factorial : ℝ) ≤ x ^ (N + 1 + j) / (((N + 1).factorial : ℝ) * (N + 2 : ℝ) ^ j) := div_le_div_of_nonneg_left hnum hpos hfac _ = (x ^ (N + 1) / ((N + 1).factorial : ℝ)) * (x / (N + 2 : ℝ)) ^ j := by rw [div_pow, pow_add] field_simp _ = a * (1 - q)⁻¹ := by rw [tsum_mul_left, tsum_geometric_of_abs_lt_one hq] _ = (x ^ (N + 1) / ((N + 1).factorial : ℝ)) * ((N + 2 : ℝ) / ((N + 2 : ℝ) - x)) := by dsimp [a, q] have hn : (N + 2 : ℝ) ≠ 0 := by positivity have hd : (N + 2 : ℝ) - x ≠ 0 := ne_of_gt (sub_pos.mpr hden) field_simp

theorem abs_sq_sub_gaussianWeightedTaylor_oneHundred_le {r : ℝ} (hr : r ∈ Icc (0 : ℝ) 5) : |r ^ 2 - exp (-r ^ 2) * (secondMomentTaylorPolynomial 100).eval r| ≤ 25 * (25 ^ 101 / ((101 : ℕ).factorial : ℝ)) * (102 / 77) := by have hx0 : 0 ≤ r ^ 2 := sq_nonneg r have hprev := abs_secondMoment_exp_sub_taylor_oneHundred_le hr have hexp : exp (-r ^ 2) * exp (r ^ 2) = 1 := by rw [← Real.exp_add] simp have hid : r ^ 2 - exp (-r ^ 2) * (secondMomentTaylorPolynomial 100).eval r = exp (-r ^ 2) * (r ^ 2 * exp (r ^ 2) - (secondMomentTaylorPolynomial 100).eval r) := by calc r ^ 2 - exp (-r ^ 2) * (secondMomentTaylorPolynomial 100).eval r = (exp (-r ^ 2) * exp (r ^ 2)) * r ^ 2 - exp (-r ^ 2) * (secondMomentTaylorPolynomial 100).eval r := by rw [hexp]; ring _ = _ := by ring rw [hid, abs_mul, abs_of_pos (Real.exp_pos _)] calc exp (-r ^ 2) * |r ^ 2 * exp (r ^ 2) - (secondMomentTaylorPolynomial 100).eval r| ≤ 1 * |r ^ 2 * exp (r ^ 2) - (secondMomentTaylorPolynomial 100).eval r| := mul_le_mul_of_nonneg_right (Real.exp_le_one_iff.mpr (by linarith)) (abs_nonneg _) _ ≤ 25 * (25 ^ 101 / ((101 : ℕ).factorial : ℝ)) * (102 / 77) := by simpa using hprev

theorem secondMoment_sub_le_of_taylor_error (μ ν : Measure ℝ) [IsProbabilityMeasure μ] [IsProbabilityMeasure ν] (R ε : ℝ) (N : ℕ) (hμ : μ (Icc (0 : ℝ) R)ᶜ = 0) (hν : ν (Icc (0 : ℝ) R)ᶜ = 0) (happrox : ∀ r ∈ Icc (0 : ℝ) R, |r ^ 2 - gaussianWeight r * (secondMomentTaylorPolynomial N).eval r| ≤ ε) : |(∫ r, r ^ 2 ∂μ) - ∫ r, r ^ 2 ∂ν| ≤ secondMomentRecoveryCoefficient N * ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ + 2 * ε := by let S : Set ℝ := Icc 0 R let q : ℝ → ℝ := fun r => gaussianWeight r * (secondMomentTaylorPolynomial N).eval r have hμS : ∀ᵐ r ∂μ, r ∈ S := by rw [ae_iff] exact hμ have hνS : ∀ᵐ r ∂ν, r ∈ S := by rw [ae_iff] exact hν have hμrestrict : μ.restrict S = μ := Measure.restrict_eq_self_of_ae_mem hμS have hνrestrict : ν.restrict S = ν := Measure.restrict_eq_self_of_ae_mem hνS have hsquare_cont : Continuous (fun r : ℝ => r ^ 2) := by fun_prop have hq_cont : Continuous q := by dsimp [q] unfold gaussianWeight fun_prop have hsquareμ : Integrable (fun r : ℝ => r ^ 2) μ := by have h := hsquare_cont.integrableOn_Icc (a := (0 : ℝ)) (b := R) (μ := μ) change Integrable (fun r : ℝ => r ^ 2) (μ.restrict S) at h rwa [hμrestrict] at h have hsquareν : Integrable (fun r : ℝ => r ^ 2) ν := by have h := hsquare_cont.integrableOn_Icc (a := (0 : ℝ)) (b := R) (μ := ν) change Integrable (fun r : ℝ => r ^ 2) (ν.restrict S) at h rwa [hνrestrict] at h have hqμ : Integrable q μ := by have h := hq_cont.integrableOn_Icc (a := (0 : ℝ)) (b := R) (μ := μ) change Integrable q (μ.restrict S) at h rwa [hμrestrict] at h have hqν : Integrable q ν := by have h := hq_cont.integrableOn_Icc (a := (0 : ℝ)) (b := R) (μ := ν) change Integrable q (ν.restrict S) at h rwa [hνrestrict] at h have herrorμ : |(∫ r, r ^ 2 ∂μ) - ∫ r, q r ∂μ| ≤ ε := by have hbound : ∀ᵐ r ∂μ, ‖r ^ 2 - q r‖ ≤ ε := hμS.mono fun r hr => by simpa only [S, q, Real.norm_eq_abs] using happrox r hr rw [← integral_sub hsquareμ hqμ, ← Real.norm_eq_abs] simpa using norm_integral_le_of_norm_le_const hbound have herrorν : |(∫ r, r ^ 2 ∂ν) - ∫ r, q r ∂ν| ≤ ε := by have hbound : ∀ᵐ r ∂ν, ‖r ^ 2 - q r‖ ≤ ε := hνS.mono fun r hr => by simpa only [S, q, Real.norm_eq_abs] using happrox r hr rw [← integral_sub hsquareν hqν, ← Real.norm_eq_abs] simpa using norm_integral_le_of_norm_le_const hbound have hq_integral (m : Measure ℝ) [IsFiniteMeasure m] : (∫ r, q r ∂m) = ∑ k ∈ range (N + 1), (1 / (k.factorial : ℝ)) * ∫ r, gaussianWeightedMonomial (2 * k + 2) r ∂m := by simp only [q, secondMomentTaylorPolynomial_eval, Finset.mul_sum] rw [integral_finset_sum] · apply Finset.sum_congr rfl intro k hk have hfun : (fun r : ℝ => gaussianWeight r * (r ^ (2 * k + 2) / (k.factorial : ℝ))) = fun r => (1 / (k.factorial : ℝ)) * gaussianWeightedMonomial (2 * k + 2) r := by funext r simp only [gaussianWeightedMonomial] ring rw [hfun, integral_const_mul] · intro k hk rw [show (fun r : ℝ => gaussianWeight r * (r ^ (2 * k + 2) / (k.factorial : ℝ))) = fun r => (1 / (k.factorial : ℝ)) * gaussianWeightedMonomial (2 * k + 2) r by funext r simp only [gaussianWeightedMonomial] ring] exact (integrable_gaussianWeightedMonomial m (2 * k + 2)).const_mul _ have hmiddle : |(∫ r, q r ∂μ) - ∫ r, q r ∂ν| ≤ secondMomentRecoveryCoefficient N * ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ := by rw [hq_integral μ, hq_integral ν, ← Finset.sum_sub_distrib] simpa only [mul_sub] using abs_secondMomentWeightedSum_sub_le μ ν N calc |(∫ r, r ^ 2 ∂μ) - ∫ r, r ^ 2 ∂ν| ≤ |(∫ r, r ^ 2 ∂μ) - ∫ r, q r ∂μ| + |(∫ r, q r ∂μ) - ∫ r, q r ∂ν| + |(∫ r, q r ∂ν) - ∫ r, r ^ 2 ∂ν| := by calc _ ≤ |(∫ r, r ^ 2 ∂μ) - ∫ r, q r ∂μ| + |(∫ r, q r ∂μ) - ∫ r, r ^ 2 ∂ν| := abs_sub_le _ _ _ _ ≤ _ := by have htri : |(∫ r, q r ∂μ) - ∫ r, r ^ 2 ∂ν| ≤ |(∫ r, q r ∂μ) - ∫ r, q r ∂ν| + |(∫ r, q r ∂ν) - ∫ r, r ^ 2 ∂ν| := abs_sub_le _ _ _ linarith _ ≤ ε + secondMomentRecoveryCoefficient N * ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ + ε := by gcongr simpa only [abs_sub_comm] using herrorν _ = secondMomentRecoveryCoefficient N * ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ + 2 * ε := by ring

theorem gaussian_meanEmbedding_norm_lowerBound_of_secondMoment_gap (μ ν : Measure ℝ) [IsProbabilityMeasure μ] [IsProbabilityMeasure ν] (hμ : μ (Icc (0 : ℝ) 5)ᶜ = 0) (hν : ν (Icc (0 : ℝ) 5)ᶜ = 0) (hgap : (1 / 10 ^ 15 : ℝ) < |(∫ r, r ^ 2 ∂μ) - ∫ r, r ^ 2 ∂ν|) : 0 < (|(∫ r, r ^ 2 ∂μ) - ∫ r, r ^ 2 ∂ν| - (1 / 10 ^ 15 : ℝ)) / 5151 ∧ (|(∫ r, r ^ 2 ∂μ) - ∫ r, r ^ 2 ∂ν| - (1 / 10 ^ 15 : ℝ)) / 5151 ≤ ‖meanEmbedding gaussianFeatureMap μ - meanEmbedding gaussianFeatureMap ν‖ := by have hbound := gaussian_meanEmbedding_secondMoment_stability_Icc_zero_five μ ν hμ hν constructor · positivity · nlinarith