diff --git a/LMLExplore/MarkovDecisionProcess/MarkovPolicy.lean b/LMLExplore/MarkovDecisionProcess/MarkovPolicy.lean index c3546e76..8626ba79 100644 --- a/LMLExplore/MarkovDecisionProcess/MarkovPolicy.lean +++ b/LMLExplore/MarkovDecisionProcess/MarkovPolicy.lean @@ -32,8 +32,9 @@ policy `π`, these marginals are given by powers of two kernels: `(M.saKernel π ^ n) ∘ₘ (μ ⊗ₘ π)`. Integrating a function against `M.stateKernel π s` averages over the actions `a` of `π s` its -integral against `M.P (s, a)` (`integral_stateKernel`); integrating a function against -`M.saKernel π p` integrates over `M.P p` its average over the actions of `π` (`integral_saKernel`). +integral against `M.P (s, a)` (`integral_stateKernel`, and `integral_stateKernel_eq_sum` for +finitely many actions); integrating a function against `M.saKernel π p` integrates over `M.P p` +its average over the actions of `π` (`integral_saKernel`). The deterministic stationary policies `policyAlg π hπ` are Markov policies (`policyAlg_eq_markov`). @@ -139,6 +140,16 @@ lemma integral_stateKernel {s : 𝓢} {v : 𝓢 → ℝ} (hv : Integrable v (M.s rw [Kernel.integral_comp hv, Kernel.prod_apply, Kernel.id_apply, Measure.dirac_prod, integral_map measurable_prodMk_left.aemeasurable h] +omit [IsSFiniteKernel π] in +variable {M π} in +/-- For finitely many actions, integrating against the state kernel of `π` from `s` is the average +over the actions `a`, with the weights `π(a | s)`, of the integrals against `M.P (s, a)`. -/ +lemma integral_stateKernel_eq_sum [Fintype 𝓐] [MeasurableSingletonClass 𝓐] [IsFiniteKernel π] + {s : 𝓢} {v : 𝓢 → ℝ} (hv : Integrable v (M.stateKernel π s)) : + ∫ s', v s' ∂(M.stateKernel π s) = ∑ a, (π s).real {a} * ∫ s', v s' ∂(M.P (s, a)) := by + rw [integral_stateKernel hv, integral_fintype .of_finite] + simp only [smul_eq_mul] + @[simp] lemma saKernel_apply (p : 𝓢 × 𝓐) : M.saKernel π p = M.P p ⊗ₘ π := by rw [saKernel_def, Kernel.comp_apply, Measure.compProd_eq_comp_prod] diff --git a/LMLExplore/MarkovDecisionProcess/Simulation.lean b/LMLExplore/MarkovDecisionProcess/Simulation.lean index 78f8266d..bc39986a 100644 --- a/LMLExplore/MarkovDecisionProcess/Simulation.lean +++ b/LMLExplore/MarkovDecisionProcess/Simulation.lean @@ -366,8 +366,8 @@ private lemma abs_integral_resolvent_le_of_abs_le_sqrt_variance (hγ : γ < 1) set σ : 𝓢 × 𝓐 → ℝ := fun q ↦ √(Var[M.value γ π; M.P q]) have hσ : Measurable σ := (hv.variance_kernel (η := M.P)).sqrt have hσB (q : 𝓢 × 𝓐) : |σ q| ≤ C / (1 - γ) := by - rw [abs_of_nonneg (Real.sqrt_nonneg _), Real.sqrt_le_left hB] - exact variance_le_sq_of_abs_le hv.aemeasurable (.of_forall hvB) + rw [abs_of_nonneg (Real.sqrt_nonneg _)] + exact sqrt_variance_le_of_abs_le hv.aemeasurable hB (.of_forall hvB) have hσi := (integrable_of_abs_le hσ hσB ρ).const_mul a have h_int : ∫ q, σ q ∂ρ ≤ C * √(2 / (1 - γ) ^ 3) := integral_sqrt_variance_value_resolvent_le hγ hR p @@ -537,8 +537,7 @@ theorem abs_optQValue_sub_optQValue_le_of_le_sqrt_variance (hγ : γ < 1) a * (C / (1 - γ)) + b := by refine (hε q).trans ?_ gcongr - rw [Real.sqrt_le_left hB] - exact variance_le_sq_of_abs_le hv.aemeasurable (.of_forall hvB) + exact sqrt_variance_le_of_abs_le hv.aemeasurable hB (.of_forall hvB) -- the greedy policies for the optimal values are optimal stationary policies have hg := measurable_greedy (M := M) (γ := γ) hv have hg' := measurable_greedy (M := M') (γ := γ) (measurable_optValue hγ hR') diff --git a/LMLExplore/Probability/Kernel/Empirical.lean b/LMLExplore/Probability/Kernel/Empirical.lean index 2ecf8f6b..8ba826b4 100644 --- a/LMLExplore/Probability/Kernel/Empirical.lean +++ b/LMLExplore/Probability/Kernel/Empirical.lean @@ -82,6 +82,8 @@ This is the decoupling device of [AKY20, Lemma 6]. `LMLExplore.le_streamMeasure_forall_abs_integral_sub_integral_empiricalKernel_lt_sqrt_add`: the empirical Bernstein inequality for the empirical means of a bounded function, for one key and uniformly over finitely many keys. +* `LMLExplore.measureReal_sqrt_add_le_abs_sub_inv_mul_sum_le`: Bernstein's inequality for the + empirical means `N⁻¹ ∑_{i < N} Z (ω i)` of a bounded function `Z` of the samples. * `LMLExplore.measureReal_sqrt_variance_add_le_abs_integral_sub_integral_empiricalKernel_le`: Bernstein's inequality for the empirical means of a bounded function, with its variance under `ν a`. @@ -296,6 +298,58 @@ theorem measureReal_sqrt_add_le_abs_integral_sub_integral_empiricalKernel_le (ν integral_empiricalKernel N ω a hf.stronglyMeasurable, smul_eq_mul, inv_mul_eq_div] rw [abs_sub_comm] +omit [MeasurableSingletonClass α] [Countable α] in +/-- **Bernstein's inequality for the empirical means of a function of the samples.** Under the +stream model `streamMeasure ν` of a Markov kernel `ν`, let `Z` be a measurable function of the +samples `ρ : α → β` with mean `m`, `|Z - m| ≤ b` almost everywhere for `b > 0`, and second moment +`∫ (Z - m) ^ 2 ≤ σ` under the law `Measure.infinitePi ν` of the samples. For `N ≠ 0` and every +`x`, the mean `N⁻¹ ∑_{i < N} Z (ω i)` deviates from `m` by at least `√(2 σ x / N) + b x / (3 N)` +with probability at most `2 exp (-x)`. -/ +theorem measureReal_sqrt_add_le_abs_sub_inv_mul_sum_le (ν : Kernel α β) [IsMarkovKernel ν] + {Z : (α → β) → ℝ} (hZ : Measurable Z) {m b σ : ℝ} (hb : 0 < b) + (hZm : ∫ ρ, Z ρ ∂(Measure.infinitePi ν) = m) + (hZb : ∀ᵐ ρ ∂(Measure.infinitePi ν), |Z ρ - m| ≤ b) + (hZσ : ∫ ρ, (Z ρ - m) ^ 2 ∂(Measure.infinitePi ν) ≤ σ) (N : ℕ) [NeZero N] (x : ℝ) : + (streamMeasure ν).real {ω | √(2 * σ * x / N) + b * x / (3 * N) ≤ + |m - (N : ℝ)⁻¹ * ∑ i ∈ range N, Z (ω i)|} ≤ 2 * exp (-x) := by + have hN : (0 : ℝ) < N := by exact_mod_cast Nat.pos_of_ne_zero (NeZero.ne N) + -- the centred samples `Z (ω i) - m` are independent, with mean zero and second moment `≤ σ` + have h_law (i : ℕ) := hasLaw_eval_streamMeasure ν i + have h_indep : iIndepFun (fun i (ω : ℕ → α → β) ↦ Z (ω i) - m) (streamMeasure ν) := + (iIndepFun_eval_streamMeasure' ν).comp (fun _ ρ ↦ Z ρ - m) fun _ ↦ hZ.sub_const m + have hZi : Integrable Z (Measure.infinitePi ν) := + .of_bound hZ.aestronglyMeasurable (b + |m|) (hZb.mono fun ρ hρ ↦ by + rw [Real.norm_eq_abs] + linarith [abs_sub_abs_le_abs_sub (Z ρ) m]) + have h_mean (i : ℕ) : ∫ ω, Z (ω i) - m ∂(streamMeasure ν) = 0 := by + refine ((h_law i).integral_comp (f := fun ρ ↦ Z ρ - m) + (hZ.sub_const m).aestronglyMeasurable).trans ?_ + rw [integral_sub hZi (integrable_const m), hZm] + simp + have h_sq (i : ℕ) : ∫ ω, (Z (ω i) - m) ^ 2 ∂(streamMeasure ν) ≤ σ := + ((h_law i).integral_comp (f := fun ρ ↦ (Z ρ - m) ^ 2) + ((hZ.sub_const m).pow_const 2).aestronglyMeasurable).trans_le hZσ + have h := h_indep.measureReal_sqrt_add_le_abs_sum_le_of_abs_le (s := range N) (b := b) + (v := N * σ) (L := x) (fun i _ ↦ (hZ.sub_const m).comp_aemeasurable (h_law i).aemeasurable) + (fun i _ ↦ h_mean i) hb + (fun i _ ↦ ae_of_ae_map (p := fun ρ ↦ |Z ρ - m| ≤ b) (h_law i).aemeasurable + ((h_law i).map_eq ▸ hZb)) + ((sum_le_sum fun i _ ↦ h_sq i).trans_eq (by simp)) + refine le_trans (measureReal_mono fun ω (hω : _ ≤ _) ↦ ?_) h + -- the centred sum is `-N` times the deviation of the mean + have h_sum : ∑ i ∈ range N, (Z (ω i) - m) = + -(N * (m - (N : ℝ)⁻¹ * ∑ i ∈ range N, Z (ω i))) := by + simp only [sum_sub_distrib, sum_const, card_range, nsmul_eq_mul] + field_simp + ring + have h_sqrt : √(2 * (N * σ) * x) = N * √(2 * σ * x / N) := by + rw [show 2 * (N * σ) * x = N ^ 2 * (2 * σ * x / N) by field_simp, + Real.sqrt_mul (sq_nonneg _), Real.sqrt_sq hN.le] + simp only [Set.mem_ofPred_eq, h_sum, abs_neg, abs_mul, abs_of_pos hN, h_sqrt] + calc N * √(2 * σ * x / N) + b * x / 3 = N * (√(2 * σ * x / N) + b * x / (3 * N)) := by + field_simp + _ ≤ N * |m - (N : ℝ)⁻¹ * ∑ i ∈ range N, Z (ω i)| := by gcongr + /-- **Bernstein's inequality for the empirical kernel.** Under the stream model `streamMeasure ν` of a Markov kernel `ν`, let `f` be a measurable function with values in `[l, u]`, `l < u`, `ν a`-almost everywhere. For `N ≠ 0` and every `x`, the mean of `f` over the first `N` samples of @@ -306,49 +360,24 @@ theorem measureReal_sqrt_variance_add_le_abs_integral_sub_integral_empiricalKern {a : α} (hfb : ∀ᵐ y ∂(ν a), f y ∈ Set.Icc l u) (N : ℕ) [NeZero N] (x : ℝ) : (streamMeasure ν).real {ω | √(2 * Var[f; ν a] * x / N) + (u - l) * x / (3 * N) ≤ |∫ y, f y ∂(ν a) - ∫ y, f y ∂(empiricalKernel N ω a)|} ≤ 2 * exp (-x) := by - have hN : (0 : ℝ) < N := by exact_mod_cast Nat.pos_of_ne_zero (NeZero.ne N) - set m := ∫ y, f y ∂(ν a) + -- the samples `f (ρ a)` have mean `∫ f d(ν a)` and variance `Var[f; ν a]` + have h_law := hasLaw_eval_infinitePi ν a have hf_int : Integrable f (ν a) := Integrable.of_mem_Icc l u hf.aemeasurable hfb - have hm : m ∈ Set.Icc l u := by + have hm : ∫ y, f y ∂(ν a) ∈ Set.Icc l u := by constructor · simpa using integral_mono_ae (integrable_const l) hf_int (hfb.mono fun y hy ↦ hy.1) · simpa using integral_mono_ae hf_int (integrable_const u) (hfb.mono fun y hy ↦ hy.2) - -- the centred samples `f (ω i a) - m` are independent, with mean zero and variance `Var[f; ν a]` - have h_law (i : ℕ) := hasLaw_eval_eval_streamMeasure ν i a - have h_indep : iIndepFun (fun i (ω : ℕ → α → β) ↦ f (ω i a) - m) (streamMeasure ν) := - (iIndepFun_eval_streamMeasure'' ν a).comp (fun _ y ↦ f y - m) fun _ ↦ hf.sub_const m - have h_mean (i : ℕ) : ∫ ω, f (ω i a) - m ∂(streamMeasure ν) = 0 := by - refine ((h_law i).integral_comp (f := fun y ↦ f y - m) - (hf.sub_const m).aestronglyMeasurable).trans ?_ - rw [integral_sub hf_int (integrable_const m)] - simp [m] - have h_sq (i : ℕ) : ∫ ω, (f (ω i a) - m) ^ 2 ∂(streamMeasure ν) = Var[f; ν a] := by - refine ((h_law i).integral_comp (f := fun y ↦ (f y - m) ^ 2) - ((hf.sub_const m).pow_const 2).aestronglyMeasurable).trans ?_ - rw [variance_eq_integral hf.aemeasurable] - have h := h_indep.measureReal_sqrt_add_le_abs_sum_le_of_abs_le (s := range N) (b := u - l) - (v := N * Var[f; ν a]) (L := x) - (fun i _ ↦ (hf.sub_const m).comp_aemeasurable (h_law i).aemeasurable) (fun i _ ↦ h_mean i) - (sub_pos.2 hlu) - (fun i _ ↦ (ae_of_ae_map (h_law i).aemeasurable ((h_law i).map_eq ▸ hfb)).mono fun ω hω ↦ - abs_le.2 ⟨by linarith [hω.1, hm.2], by linarith [hω.2, hm.1]⟩) - (by simp [h_sq]) - refine le_trans (measureReal_mono fun ω (hω : _ ≤ _) ↦ ?_) h - -- the centred sum is `-N` times the deviation of the empirical mean - have h_sum : ∑ i ∈ range N, (f (ω i a) - m) = - -(N * (m - ∫ y, f y ∂(empiricalKernel N ω a))) := by - rw [integral_empiricalKernel N ω a hf.stronglyMeasurable] - simp [sum_sub_distrib] - field_simp - ring - have h_sqrt : √(2 * (N * Var[f; ν a]) * x) = N * √(2 * Var[f; ν a] * x / N) := by - rw [show 2 * (N * Var[f; ν a]) * x = N ^ 2 * (2 * Var[f; ν a] * x / N) by field_simp, - Real.sqrt_mul (sq_nonneg _), Real.sqrt_sq hN.le] - simp only [Set.mem_ofPred_eq, h_sum, abs_neg, abs_mul, abs_of_pos hN, h_sqrt] - calc N * √(2 * Var[f; ν a] * x / N) + (u - l) * x / 3 = - N * (√(2 * Var[f; ν a] * x / N) + (u - l) * x / (3 * N)) := by - field_simp - _ ≤ N * |m - ∫ y, f y ∂(empiricalKernel N ω a)| := by gcongr + have h := measureReal_sqrt_add_le_abs_sub_inv_mul_sum_le ν (Z := fun ρ ↦ f (ρ a)) + (hf.comp (measurable_pi_apply a)) (sub_pos.2 hlu) + (h_law.integral_comp hf.aestronglyMeasurable) + ((ae_of_ae_map h_law.aemeasurable (h_law.map_eq ▸ hfb)).mono fun ρ hρ ↦ + abs_le.2 ⟨by linarith [hρ.1, hm.2], by linarith [hρ.2, hm.1]⟩) + ((h_law.integral_comp (f := fun y ↦ (f y - ∫ y, f y ∂(ν a)) ^ 2) + ((hf.sub_const _).pow_const 2).aestronglyMeasurable).trans_le + (variance_eq_integral hf.aemeasurable).ge) N x + refine le_trans (measureReal_mono fun ω hω ↦ ?_) h + rwa [Set.mem_ofPred_eq, integral_empiricalKernel N ω a hf.stronglyMeasurable, smul_eq_mul] + at hω /-- **Bernstein's inequality for the empirical means of a random function independent of the samples.** Under the stream model `streamMeasure ν` with countably many values, let `W ω` be a diff --git a/LMLExplore/Probability/Kernel/Resolvent.lean b/LMLExplore/Probability/Kernel/Resolvent.lean index 6a241e3d..13c9c0a7 100644 --- a/LMLExplore/Probability/Kernel/Resolvent.lean +++ b/LMLExplore/Probability/Kernel/Resolvent.lean @@ -15,6 +15,7 @@ public import Mathlib.Probability.Kernel.Composition.MeasureComp public import Mathlib.Probability.Kernel.WithDensity public import Mathlib.Probability.Moments.Variance import LMLExplore.Analysis.Convex.Integral +import Mathlib.Analysis.Complex.ExponentialBounds import LMLExplore.Probability.Moments.Variance /-! @@ -55,7 +56,11 @@ function `2 B V - V ^ 2` is a supersolution for `γ ^ 2 Var_κ(V)`, so that `R_γ √Var_κ(V) ≤ B / (γ √(1 - γ))` [LWCGC24, Lemma 11]. Bounds of this kind control the error `γ R_γ(κ') (κ' - κ) V` of the simulation identity when the deviation `(κ' - κ) V` is bounded by a multiple of the standard deviation `√Var_κ(V)`, as Bernstein's inequality bounds it for an -empirical kernel `κ'`. +empirical kernel `κ'`. Iterating this control along the functions `V 0 = R_γ f` and +`V (l + 1) = R_γ √Var_κ(V l)`, each bounded through Lemma 11 by the previous one, compares +`R_γ(κ') f` with `R_γ(κ) f` up to `3 C √(ε / (1 - γ) ^ 3)` for `0 ≤ f ≤ C`: for `0 < γ < 1` and +`9 ε ≤ 1 - γ`, as soon as `κ'` deviates from `κ` by `√(ε Var_κ(V l))` plus a lower-order term on +the levels `l ≤ m`, for some `m ≥ log (1 - γ)⁻¹` [LWCGC24, Lemma 2]. In a Markov decision process, the value of a stationary policy is the resolvent of the state kernel of the policy applied to its reward, and the discounted occupancy measure of a stationary @@ -101,10 +106,19 @@ equation of these objects. supersolutions, `f + γ κ u ≤ u → R_γ f ≤ u` and `u ≤ f + γ κ u → u ≤ R_γ f`. * `ProbabilityTheory.Kernel.integral_resolvent_sub_integral_resolvent`: the simulation identity with two rewards `R_γ(κ) f - R_γ(κ') f' = R_γ(κ') ((f - f') + γ (κ - κ') R_γ(κ) f)`. +* `ProbabilityTheory.Kernel.abs_integral_resolvent_sub_integral_resolvent_le`: + the simulation bound `|R_γ(κ') f - R_γ(κ) f| ≤ γ R_γ(κ') e` when `|(κ' - κ) R_γ(κ) f| ≤ e`. * `ProbabilityTheory.Kernel.mul_integral_variance_integral_resolvent_le`, `ProbabilityTheory.Kernel.integral_sqrt_variance_integral_resolvent_le`: for a nonnegative bounded measurable `f` with `V = R_γ f ≤ B`, `γ ^ 2 R_γ Var_κ(V) ≤ 2 B V - V ^ 2` and `R_γ √Var_κ(V) ≤ B / (γ √(1 - γ))` [LWCGC24, Lemma 11]. +* `ProbabilityTheory.Kernel.measurable_of_eq_integral_resolvent_sqrt_variance`, + `ProbabilityTheory.Kernel.mem_Icc_of_eq_integral_resolvent_sqrt_variance`: the iterated + resolvents `V (l + 1) = R_γ √Var_κ(V l)` of the one-step standard deviations of `V 0 = R_γ f` are + measurable, with `0 ≤ V l ≤ C / ((1 - γ) (γ √(1 - γ)) ^ l)` for `0 ≤ f ≤ C`. +* `ProbabilityTheory.Kernel.abs_integral_resolvent_sub_integral_resolvent_le_of_abs_sub_le`: the + comparison of the resolvents of two kernels with Bernstein-type deviations on these functions + [LWCGC24, Lemma 2]. * `MeasureTheory.Measure.resolvent_comp_eq_add`: for every kernel and every `γ`, `resolvent γ κ ∘ₘ μ` solves the flow equation `ν = μ + γ • (κ ∘ₘ ν)`; `MeasureTheory.Measure.eq_resolvent_comp_of_eq_add`: for a Markov kernel and `γ < 1`, it is the @@ -120,8 +134,8 @@ equation of these objects. * [A. Agarwal, N. Jiang, S. M. Kakade and W. Sun, *Reinforcement Learning: Theory and Algorithms*, 2022][AJKS], Lemma 2.2 * [G. Li, Y. Wei, Y. Chi, Y. Gu and Y. Chen, *Breaking the sample size barrier in model-based - reinforcement learning with a generative model*, Operations Research, 2024][LWCGC24], Lemma 11 - (arXiv:2005.12900) + reinforcement learning with a generative model*, Operations Research, 2024][LWCGC24], Lemmas 2 + and 11 (arXiv:2005.12900) -/ public section @@ -466,6 +480,24 @@ theorem integral_resolvent_sub_integral_resolvent (κ' : Kernel α α) [IsMarkov congr with z rw [integral_resolvent_eq_add (κ := κ) hγ hf hfB z] ring +/-- **The simulation bound**: for Markov kernels `κ` and `κ'`, `γ < 1` and a bounded measurable +`f`, if the one-step deviation `|κ' V - κ V|` of `V = R_γ(κ) f` is bounded by a function `e` +integrable against `R_γ(κ') x`, then `|R_γ(κ') f - R_γ(κ) f| ≤ γ R_γ(κ') e` at `x`. -/ +theorem abs_integral_resolvent_sub_integral_resolvent_le + (κ' : Kernel α α) [IsMarkovKernel κ'] (hγ : γ < 1) {f e : α → ℝ} {B : ℝ} (hf : Measurable f) + (hfB : ∀ z, |f z| ≤ B) {x : α} (he : Integrable e (resolvent γ κ' x)) + (h : ∀ z, |∫ w, (∫ y, f y ∂(resolvent γ κ w)) ∂(κ' z) - + ∫ w, (∫ y, f y ∂(resolvent γ κ w)) ∂(κ z)| ≤ e z) : + |∫ z, f z ∂(resolvent γ κ' x) - ∫ z, f z ∂(resolvent γ κ x)| ≤ + γ * ∫ z, e z ∂(resolvent γ κ' x) := by + rw [abs_sub_comm, integral_resolvent_sub_integral_resolvent κ' hγ hf hfB hf hfB x, + ← integral_const_mul] + simp_rw [sub_self, zero_add] + refine abs_integral_le_integral_abs.trans (integral_mono_of_nonneg + (.of_forall fun _ ↦ abs_nonneg _) (he.const_mul _) (.of_forall fun z ↦ ?_)) + dsimp only + rw [abs_mul, NNReal.abs_eq, abs_sub_comm] + exact mul_le_mul_of_nonneg_left (h z) γ.coe_nonneg /-! ### The resolvent of the one-step variances -/ @@ -548,6 +580,263 @@ theorem integral_sqrt_variance_integral_resolvent_le (hγ₀ : 0 < γ) (hγ : γ ring] gcongr +/-! ### The resolvent of a kernel with Bernstein-type deviations -/ + +omit [IsMarkovKernel κ] in +/-- The iterated resolvents of the one-step standard deviations are measurable: if +`V 0 = R_γ f` for a measurable `f` and `V (l + 1) = R_γ √Var_κ(V l)` for every `l`, then every +`V l` is measurable. -/ +theorem measurable_of_eq_integral_resolvent_sqrt_variance [IsSFiniteKernel κ] {f : α → ℝ} + (hf : Measurable f) {V : ℕ → α → ℝ} (hV₀ : ∀ x, V 0 x = ∫ y, f y ∂(resolvent γ κ x)) + (hV : ∀ l x, V (l + 1) x = ∫ y, √Var[V l; κ y] ∂(resolvent γ κ x)) (l : ℕ) : + Measurable (V l) := by + induction l with + | zero => rw [funext hV₀]; exact measurable_integral_resolvent hf + | succ l ih => rw [funext (hV l)]; exact measurable_integral_resolvent ih.variance_kernel.sqrt + +/-- **A priori bounds on the iterated resolvents of the one-step standard deviations** +[LWCGC24, (73)]: for `0 < γ < 1` and a measurable `f` with `0 ≤ f ≤ C`, if `V 0 = R_γ f` and +`V (l + 1) = R_γ √Var_κ(V l)` for every `l`, then `0 ≤ V l ≤ C / ((1 - γ) (γ √(1 - γ)) ^ l)`. +The bound of `V 0` is `C / (1 - γ)`, and each step multiplies it by `1 / (γ √(1 - γ))` +(`ProbabilityTheory.Kernel.integral_sqrt_variance_integral_resolvent_le`). -/ +theorem mem_Icc_of_eq_integral_resolvent_sqrt_variance (hγ₀ : 0 < γ) (hγ : γ < 1) {f : α → ℝ} + {C : ℝ} (hf : Measurable f) (hf₀ : ∀ x, 0 ≤ f x) (hfC : ∀ x, f x ≤ C) {V : ℕ → α → ℝ} + (hV₀ : ∀ x, V 0 x = ∫ y, f y ∂(resolvent γ κ x)) + (hV : ∀ l x, V (l + 1) x = ∫ y, √Var[V l; κ y] ∂(resolvent γ κ x)) (l : ℕ) (x : α) : + V l x ∈ Icc 0 (C / ((1 - γ) * (γ * √(1 - γ)) ^ l)) := by + have hVm := measurable_of_eq_integral_resolvent_sqrt_variance hf hV₀ hV + -- `V l` is the resolvent of a nonnegative bounded measurable function `g` + suffices h : ∃ g : α → ℝ, Measurable g ∧ (∀ y, 0 ≤ g y) ∧ (∃ B, ∀ y, |g y| ≤ B) ∧ + (∀ x, V l x = ∫ y, g y ∂(resolvent γ κ x)) ∧ + ∀ x, V l x ∈ Icc 0 (C / ((1 - γ) * (γ * √(1 - γ)) ^ l)) from h.choose_spec.2.2.2.2 x + induction l with + | zero => + refine ⟨f, hf, hf₀, ⟨C, fun y ↦ abs_le.2 ⟨by linarith [hf₀ y, hfC y], hfC y⟩⟩, hV₀, + fun x ↦ ⟨?_, ?_⟩⟩ + · rw [hV₀] + exact integral_nonneg hf₀ + · rw [hV₀, pow_zero, mul_one] + exact (le_abs_self _).trans (abs_integral_resolvent_le hγ + (fun y ↦ abs_le.2 ⟨by linarith [hf₀ y, hfC y], hfC y⟩) x) + | succ l ih => + obtain ⟨g, hg, hg₀, ⟨B, hgB⟩, hVg, hVl⟩ := ih + have hb : 0 ≤ C / ((1 - γ) * (γ * √(1 - γ)) ^ l) := (hVl x).1.trans (hVl x).2 + have hσ (y : α) : √Var[V l; κ y] ≤ C / ((1 - γ) * (γ * √(1 - γ)) ^ l) := + sqrt_variance_le_of_abs_le (hVm l).aemeasurable hb + (.of_forall fun z ↦ abs_le.2 ⟨by linarith [(hVl z).1], (hVl z).2⟩) + refine ⟨fun y ↦ √Var[V l; κ y], (hVm l).variance_kernel.sqrt, fun y ↦ Real.sqrt_nonneg _, + ⟨_, fun y ↦ (abs_of_nonneg (Real.sqrt_nonneg _)).trans_le (hσ y)⟩, hV l, + fun x ↦ ⟨?_, ?_⟩⟩ + · rw [hV] + exact integral_nonneg fun y ↦ Real.sqrt_nonneg _ + · -- Lemma 11 for `V l = R_γ g ≤ C / ((1 - γ) (γ √(1 - γ)) ^ l)` + have h := integral_sqrt_variance_integral_resolvent_le (κ := κ) hγ₀ hγ hg hgB hg₀ + (fun y ↦ (hVg y).symm.trans_le (hVl y).2) x + rw [← funext hVg] at h + rw [hV, pow_succ, ← mul_assoc, ← div_div] + exact h + +/-- The scalar recursion behind [LWCGC24, Lemma 2, (76)–(78)]: for `0 < γ < 1`, `9 ε ≤ 1 - γ`, +`log (1 - γ)⁻¹ ≤ m` and `b l = C / ((1 - γ) (γ √(1 - γ)) ^ l)`, if the errors `d l` and the +next-level terms `e l` satisfy `d l ≤ γ √ε e l + γ ε b l / (1 - γ)` for `l ≤ m`, +`e l ≤ b (l + 1) + d (l + 1)` and `e m ≤ b m / (1 - γ)`, then `d 0 ≤ 3 C √(ε / (1 - γ) ^ 3)`. -/ +private lemma le_three_mul_sqrt_of_le_mul_add (hγ₀ : 0 < γ) (hγ : γ < 1) {C ε : ℝ} (hC : 0 ≤ C) + {m : ℕ} (hm : Real.log (1 - γ)⁻¹ ≤ m) (hε₀ : 0 ≤ ε) (hε : 9 * ε ≤ 1 - γ) {d e : ℕ → ℝ} + (hd : ∀ l ≤ m, d l ≤ + γ * √ε * e l + γ * ε * (C / ((1 - γ) * (γ * √(1 - γ)) ^ l)) / (1 - γ)) + (he : ∀ l, e l ≤ C / ((1 - γ) * (γ * √(1 - γ)) ^ (l + 1)) + d (l + 1)) + (hem : e m ≤ C / ((1 - γ) * (γ * √(1 - γ)) ^ m) / (1 - γ)) : + d 0 ≤ 3 * C * √(ε / (1 - γ) ^ 3) := by + have hγ₁ : (γ : ℝ) < 1 := hγ + have hu : 0 < 1 - (γ : ℝ) := sub_pos.2 hγ₁ + set b : ℕ → ℝ := fun l ↦ C / ((1 - γ) * (γ * √(1 - γ)) ^ l) with hb_def + -- rescaling: with `a = γ √ε` and `q = √(ε / (1 - γ))`, `a ^ l b l = (C / (1 - γ)) q ^ l` + set a : ℝ := γ * √ε with ha_def + set q : ℝ := √(ε / (1 - γ)) with hq_def + set B : ℝ := C / (1 - γ) with hB_def + have ha₀ : 0 ≤ a := by positivity + have hq₀ : 0 ≤ q := Real.sqrt_nonneg _ + have hB₀ : 0 ≤ B := by positivity + have hq : q = √ε / √(1 - γ) := Real.sqrt_div' _ hu.le + have hsqrt : 0 < √(1 - (γ : ℝ)) := Real.sqrt_pos.2 hu + have hab (l : ℕ) : a ^ l * b l = B * q ^ l := by + rw [hB_def, hq, ha_def] + simp only [hb_def] + rw [mul_pow, mul_pow, div_pow] + field_simp + have hq2 : γ * ε / (1 - γ) ≤ q ^ 2 := by + rw [hq_def, Real.sq_sqrt (div_nonneg hε₀ hu.le), mul_div_assoc] + exact mul_le_of_le_one_left (by positivity) hγ₁.le + have haq : a / (1 - γ) ≤ q / √(1 - γ) := by + rw [hq, div_div, Real.mul_self_sqrt hu.le, ha_def] + exact div_le_div_of_nonneg_right (mul_le_of_le_one_left (Real.sqrt_nonneg _) hγ₁.le) hu.le + have hq13 : q ≤ 1 / 3 := by + rw [hq_def, show (1 / 3 : ℝ) = √(1 / 9) by + rw [show (1 / 9 : ℝ) = (1 / 3) ^ 2 by norm_num, Real.sqrt_sq (by norm_num)]] + exact Real.sqrt_le_sqrt (by rw [div_le_iff₀ hu]; linarith) + set θ : ℝ := q * (1 + q) / (1 - q) with hθ_def + have hq1 : 0 < 1 - q := by linarith + have hθ : q + θ * q + q ^ 2 = θ := by + rw [hθ_def] + field_simp + ring + have hθq : q ^ 2 ≤ θ := by + rw [hθ_def, le_div_iff₀ hq1] + nlinarith + -- backward induction over the levels: `a ^ l d l ≤ B (θ q ^ l + q ^ m q / √(1 - γ))` + have hind (k : ℕ) (hk : k ≤ m) : + a ^ (m - k) * d (m - k) ≤ B * (θ * q ^ (m - k) + q ^ m * (q / √(1 - γ))) := by + induction k with + | zero => + rw [Nat.sub_zero] + have h₁ := hd m le_rfl + calc a ^ m * d m ≤ a ^ m * (a * (b m / (1 - γ)) + γ * ε * b m / (1 - γ)) := by + gcongr + exact h₁.trans (by gcongr) + _ = (a ^ m * b m) * (a / (1 - γ) + γ * ε / (1 - γ)) := by ring + _ ≤ (B * q ^ m) * (q / √(1 - γ) + θ) := by + rw [hab] + gcongr + exact hq2.trans hθq + _ = B * (θ * q ^ m + q ^ m * (q / √(1 - γ))) := by ring + | succ k ih => + have hl : m - k = m - (k + 1) + 1 := by omega + have h₁ := hd (m - (k + 1)) (Nat.sub_le _ _) + have h₂ := he (m - (k + 1)) + have h₃ := ih (by omega) + rw [hl] at h₃ + set l := m - (k + 1) + calc a ^ l * d l ≤ a ^ l * (a * (b (l + 1) + d (l + 1)) + γ * ε * b l / (1 - γ)) := by + gcongr + exact h₁.trans (by gcongr) + _ = a ^ (l + 1) * b (l + 1) + a ^ (l + 1) * d (l + 1) + + (a ^ l * b l) * (γ * ε / (1 - γ)) := by ring + _ ≤ B * q ^ (l + 1) + B * (θ * q ^ (l + 1) + q ^ m * (q / √(1 - γ))) + + (B * q ^ l) * q ^ 2 := by + rw [hab, hab] + gcongr + _ = B * (q ^ l * (q + θ * q + q ^ 2) + q ^ m * (q / √(1 - γ))) := by ring + _ = B * (θ * q ^ l + q ^ m * (q / √(1 - γ))) := by rw [hθ]; ring + -- after `m ≥ log (1 - γ)⁻¹` levels, `q ^ m ≤ exp (-m) ≤ 1 - γ` + have hqm : q ^ m ≤ 1 - γ := by + have hqe : q ≤ Real.exp (-1) := + hq13.trans (le_of_lt (by + rw [Real.exp_neg, one_div] + exact inv_strictAnti₀ (Real.exp_pos 1) Real.exp_one_lt_three)) + calc q ^ m ≤ Real.exp (-1) ^ m := pow_le_pow_left₀ hq₀ hqe m + _ = Real.exp (-m) := by rw [← Real.exp_nat_mul]; ring_nf + _ ≤ Real.exp (Real.log (1 - γ)) := Real.exp_le_exp.2 (by rw [Real.log_inv] at hm; linarith) + _ = 1 - γ := Real.exp_log hu + have h₀ := hind m le_rfl + rw [Nat.sub_self, pow_zero, one_mul, pow_zero, mul_one] at h₀ + have hfin : q ^ m * (q / √(1 - γ)) ≤ q := by + rw [mul_div_assoc', div_le_iff₀ hsqrt] + calc q ^ m * q ≤ (1 - γ) * q := by gcongr + _ = q * (√(1 - γ) * √(1 - γ)) := by rw [Real.mul_self_sqrt hu.le]; ring + _ ≤ q * √(1 - γ) := by + gcongr + rw [mul_le_iff_le_one_left hsqrt, Real.sqrt_le_one] + linarith [γ.coe_nonneg] + have hθ2 : θ ≤ 2 * q := by + rw [hθ_def, div_le_iff₀ hq1] + nlinarith + have hBq : B * q = C * √(ε / (1 - γ) ^ 3) := by + rw [hB_def, hq_def, show ε / (1 - (γ : ℝ)) ^ 3 = ε / (1 - γ) / (1 - γ) ^ 2 by + rw [div_div]; ring, Real.sqrt_div' _ (sq_nonneg _), Real.sqrt_sq hu.le] + ring + calc d 0 ≤ B * (θ + q ^ m * (q / √(1 - γ))) := h₀ + _ ≤ B * (2 * q + q) := by gcongr + _ = 3 * (B * q) := by ring + _ = 3 * C * √(ε / (1 - γ) ^ 3) := by rw [hBq]; ring + +/-- **The resolvent of a kernel with Bernstein-type deviations** [LWCGC24, Lemma 2]. Let `κ` and +`κ'` be Markov kernels, `0 < γ < 1`, and `f` a measurable function with `0 ≤ f ≤ C`. Let +`V 0 = R_γ(κ) f` and `V (l + 1) = R_γ(κ) √Var_κ(V l)` be the iterated resolvents of the one-step +standard deviations, so that `0 ≤ V l ≤ b l := C / ((1 - γ) (γ √(1 - γ)) ^ l)` +(`ProbabilityTheory.Kernel.mem_Icc_of_eq_integral_resolvent_sqrt_variance`). Suppose that `κ'` +deviates from `κ` on the first levels by a Bernstein-type bound +`|κ' V l - κ V l| ≤ √(ε Var_κ(V l)) + ε b l` for every `l ≤ m`, with `log (1 - γ)⁻¹ ≤ m` and +`9 ε ≤ 1 - γ`. Then the resolvents of `f` satisfy +`|R_γ(κ') f - R_γ(κ) f| ≤ 3 C √(ε / (1 - γ) ^ 3)`. + +With `ε = β / N`, [LWCGC24] assumes the deviation bound with the sup norm of `V l`, +which is at most `b l`, in place of `b l`, and concludes +`|R_γ(κ') f - R_γ(κ) f| ≤ 6 C √(ε / (1 - γ) ^ 3)` under the stronger condition +`16 e ^ 2 ε < 1 - γ`; the bound `γ R_γ √Var_κ(V) ≤ sup V / √(1 - γ)` of Lemma 11 with constant +`1` instead of `4` gives the constant `3`. -/ +theorem abs_integral_resolvent_sub_integral_resolvent_le_of_abs_sub_le (κ' : Kernel α α) + [IsMarkovKernel κ'] (hγ₀ : 0 < γ) (hγ : γ < 1) {f : α → ℝ} {C ε : ℝ} (hf : Measurable f) + (hf₀ : ∀ x, 0 ≤ f x) (hfC : ∀ x, f x ≤ C) {V : ℕ → α → ℝ} + (hV₀ : ∀ x, V 0 x = ∫ y, f y ∂(resolvent γ κ x)) + (hV : ∀ l x, V (l + 1) x = ∫ y, √Var[V l; κ y] ∂(resolvent γ κ x)) {m : ℕ} + (hm : Real.log (1 - γ)⁻¹ ≤ m) (hε₀ : 0 ≤ ε) (hε : 9 * ε ≤ 1 - γ) + (h : ∀ l ≤ m, ∀ x, |∫ y, V l y ∂(κ' x) - ∫ y, V l y ∂(κ x)| ≤ + √(ε * Var[V l; κ x]) + ε * (C / ((1 - γ) * (γ * √(1 - γ)) ^ l))) (x : α) : + |∫ y, f y ∂(resolvent γ κ' x) - ∫ y, f y ∂(resolvent γ κ x)| ≤ + 3 * C * √(ε / (1 - γ) ^ 3) := by + -- The error `D l = |R_γ(κ') (r l) - V l|` at level `l`, where `V l = R_γ(κ) (r l)`, is bounded + -- by the simulation bound by `γ √ε` times `R_γ(κ') √Var_κ(V l)`, which is `V (l + 1)` up to + -- the error at level `l + 1` [LWCGC24, (74)–(76)]. The scalar recursion + -- `le_three_mul_sqrt_of_le_mul_add` unrolls these bounds over the levels `l ≤ m`. + have := isFiniteKernel_resolvent (κ := κ') hγ + have hC : 0 ≤ C := (hf₀ x).trans (hfC x) + set b : ℕ → ℝ := fun l ↦ C / ((1 - γ) * (γ * √(1 - γ)) ^ l) with hb_def + have hVm := measurable_of_eq_integral_resolvent_sqrt_variance hf hV₀ hV + have hVb : ∀ l y, V l y ∈ Icc 0 (b l) := + mem_Icc_of_eq_integral_resolvent_sqrt_variance hγ₀ hγ hf hf₀ hfC hV₀ hV + have hb₀ (l : ℕ) : 0 ≤ b l := (hVb l x).1.trans (hVb l x).2 + have hσ (l : ℕ) (y : α) : √Var[V l; κ y] ≤ b l := + sqrt_variance_le_of_abs_le (hVm l).aemeasurable (hb₀ l) + (.of_forall fun z ↦ abs_le.2 ⟨by linarith [(hVb l z).1, hb₀ l], (hVb l z).2⟩) + -- the rewards `r l` of the levels: `V l = R_γ(κ) (r l)` + set r : ℕ → α → ℝ := fun l ↦ Nat.casesOn l f fun k y ↦ √Var[V k; κ y] with hr_def + have hr₀ : r 0 = f := rfl + have hr (l : ℕ) : r (l + 1) = fun y ↦ √Var[V l; κ y] := rfl + have hr_meas (l : ℕ) : Measurable (r l) := by + cases l + exacts [hf, (hVm _).variance_kernel.sqrt] + have hrB (l : ℕ) : ∃ B, ∀ y, |r l y| ≤ B := by + cases l with + | zero => rw [hr₀]; exact ⟨C, fun y ↦ abs_le.2 ⟨by linarith [hf₀ y, hfC y], hfC y⟩⟩ + | succ k => + rw [hr] + exact ⟨b k, fun y ↦ (abs_of_nonneg (Real.sqrt_nonneg _)).trans_le (hσ k y)⟩ + have hVr (l : ℕ) (x : α) : V l x = ∫ y, r l y ∂(resolvent γ κ x) := by + cases l + exacts [hV₀ x, hV _ x] + -- the error `D l = |R_γ(κ') (r l) - V l|` at level `l` + set D : ℕ → ℝ := fun l ↦ |∫ y, r l y ∂(resolvent γ κ' x) - V l x| with hD_def + have hstep (l : ℕ) (hl : l ≤ m) : + D l ≤ γ * √ε * ∫ y, r (l + 1) y ∂(resolvent γ κ' x) + γ * ε * b l / (1 - γ) := by + obtain ⟨B, hB⟩ := hrB l + obtain ⟨B', hB'⟩ := hrB (l + 1) + have hi : Integrable (fun z ↦ √ε * r (l + 1) z) (resolvent γ κ' x) := + (integrable_of_abs_le (hr_meas (l + 1)) hB' _).const_mul _ + have hdev (z : α) : |∫ w, (∫ y, r l y ∂(resolvent γ κ w)) ∂(κ' z) - + ∫ w, (∫ y, r l y ∂(resolvent γ κ w)) ∂(κ z)| ≤ √ε * r (l + 1) z + ε * b l := by + simp_rw [← hVr l] + rw [hr, ← Real.sqrt_mul hε₀] + exact h l hl z + simp only [hD_def] + rw [hVr l x] + refine (abs_integral_resolvent_sub_integral_resolvent_le κ' hγ + (hr_meas l) hB (hi.add (integrable_const (ε * b l))) hdev).trans_eq ?_ + rw [Pi.add_def, integral_add hi (integrable_const _), integral_const_mul, integral_const, + smul_eq_mul, measureReal_resolvent_univ hγ] + field_simp + -- the next level: `R_γ(κ') (r (l + 1)) ≤ V (l + 1) + D (l + 1) ≤ b (l + 1) + D (l + 1)` + have hnext (l : ℕ) : ∫ y, r (l + 1) y ∂(resolvent γ κ' x) ≤ b (l + 1) + D (l + 1) := by + have h₁ := le_abs_self (∫ y, r (l + 1) y ∂(resolvent γ κ' x) - V (l + 1) x) + linarith [(hVb (l + 1) x).2] + -- the last level: `R_γ(κ') (r (m + 1)) ≤ b m / (1 - γ)` + have hlast : ∫ y, r (m + 1) y ∂(resolvent γ κ' x) ≤ b m / (1 - γ) := + (le_abs_self _).trans (abs_integral_resolvent_le hγ + (fun y ↦ (abs_of_nonneg (Real.sqrt_nonneg _)).trans_le (hσ m y)) x) + have hD₀ : D 0 = |∫ y, f y ∂(resolvent γ κ' x) - ∫ y, f y ∂(resolvent γ κ x)| := by + simp only [hD_def, hr₀, hV₀] + rw [← hD₀] + exact le_three_mul_sqrt_of_le_mul_add hγ₀ hγ hC hm hε₀ hε hstep hnext hlast + end Markov end ProbabilityTheory.Kernel diff --git a/LMLExplore/Probability/Moments/Variance.lean b/LMLExplore/Probability/Moments/Variance.lean index 03db15fe..c5f7540b 100644 --- a/LMLExplore/Probability/Moments/Variance.lean +++ b/LMLExplore/Probability/Moments/Variance.lean @@ -20,7 +20,8 @@ This is what makes one-step variances such as `κ (V ^ 2) - (κ V) ^ 2` of a Mar themselves integrable against kernels, in variance recursions of Markov chains. A real random variable with `|X| ≤ B` almost surely has variance at most `B ^ 2`, the symmetric -case of Mathlib's `ProbabilityTheory.variance_le_sq_of_bounded`. +case of Mathlib's `ProbabilityTheory.variance_le_sq_of_bounded`, and standard deviation at most +`B`. The standard deviation `√Var[X; μ]` is a seminorm on square-integrable random variables: the Cauchy–Schwarz inequality `cov[X, Y] ^ 2 ≤ Var[X] Var[Y]` gives the triangle inequality @@ -42,6 +43,7 @@ Apache 2.0). * `Measurable.variance_kernel`: `b ↦ Var[X; η b]` is measurable. * `ProbabilityTheory.variance_le_sq_of_abs_le`: `|X| ≤ B` almost surely gives `Var[X; μ] ≤ B ^ 2`. +* `ProbabilityTheory.sqrt_variance_le_of_abs_le`: `|X| ≤ B` almost surely gives `√Var[X; μ] ≤ B`. * `ProbabilityTheory.sq_covariance_le_variance_mul_variance`: the Cauchy–Schwarz inequality `cov[X, Y; μ] ^ 2 ≤ Var[X; μ] * Var[Y; μ]`. * `ProbabilityTheory.sqrt_variance_add_le`: `√Var[X + Y; μ] ≤ √Var[X; μ] + √Var[Y; μ]`. @@ -123,6 +125,12 @@ theorem variance_le_sq_of_abs_le (hX : AEMeasurable X μ) (hXB : ∀ᵐ b ∂μ, (variance_le_sq_of_bounded (a := -B) (b := B) (hXB.mono fun _ ↦ abs_le.1) hX).trans_eq (by ring) +/-- A real random variable bounded by `B` in absolute value almost surely has standard deviation +at most `B`. -/ +theorem sqrt_variance_le_of_abs_le (hX : AEMeasurable X μ) (hB : 0 ≤ B) + (hXB : ∀ᵐ b ∂μ, |X b| ≤ B) : √Var[X; μ] ≤ B := + (Real.sqrt_le_left hB).2 (variance_le_sq_of_abs_le hX hXB) + /-- Two square-integrable real random variables within `D` of each other almost surely have standard deviations within `D`: `√Var[X; μ] ≤ √Var[Y; μ] + D`. -/ theorem sqrt_variance_le_sqrt_variance_add_of_abs_sub_le (hX : MemLp X 2 μ) (hY : MemLp Y 2 μ) @@ -135,8 +143,7 @@ theorem sqrt_variance_le_sqrt_variance_add_of_abs_sub_le (hX : MemLp X 2 μ) (hY rw [add_sub_cancel] at h refine h.trans ?_ gcongr - rw [Real.sqrt_le_left hD₀] - exact variance_le_sq_of_abs_le (hX.sub hY).aemeasurable hD + exact sqrt_variance_le_of_abs_le (hX.sub hY).aemeasurable hD₀ hD /-- **Perturbation of a Bernstein-type bound.** If `Y` satisfies the Bernstein-type bound `|∫ Y dμ - ∫ Y dν| ≤ a √Var[Y; μ] + b`, with `a ≥ 0`, between two probability measures, and `X` is diff --git a/LMLExplore/ReinforcementLearning/Generative/PolicyEvaluation.lean b/LMLExplore/ReinforcementLearning/Generative/PolicyEvaluation.lean new file mode 100644 index 00000000..907647a6 --- /dev/null +++ b/LMLExplore/ReinforcementLearning/Generative/PolicyEvaluation.lean @@ -0,0 +1,320 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import LMLExplore.ReinforcementLearning.Generative.Basic +import LMLExplore.Algebra.Order.Ring.Abs +import LMLExplore.Analysis.Convex.Integral +import LMLExplore.Analysis.Convex.Simplex + +/-! +# Planning with a generative model: policy evaluation over the full range of sample sizes + +Let `M` be an MDP with `S` states, finitely many actions and rewards in `[0, C]`, `γ < 1`, and `π` +a stationary policy fixed independently of the data. Under the generative model +`Bandits.streamMeasure M.P`, the value `v̂^π` of `π` in the empirical MDP `M.empirical N ω` of `N` +samples per state-action pair estimates the value `v^π` of `π` in `M` with accuracy +[LWCGC24, Lemma 1]: for `δ > 0` and `L = log (4 S log (e / (1 - γ)) / δ)`, as soon as +`N ≥ 18 L / (1 - γ)`, with probability at least `1 - δ`, + +`|v̂^π - v^π| ≤ 3 C √(2 L / (N (1 - γ) ^ 3))`. + +The accuracy is of order `√(L / (N (1 - γ) ^ 3))` over the whole range of sample sizes +`N ≥ 18 L / (1 - γ)`: at `N = 18 L / (1 - γ)`, the bound is the trivial bound `C / (1 - γ)`. + +The bound has a deterministic part and a probabilistic part. +* The **deterministic part** [LWCGC24, Lemma 2], + `ProbabilityTheory.Kernel.abs_integral_resolvent_sub_integral_resolvent_le_of_abs_sub_le` in + `Probability/Kernel/Resolvent.lean`, compares the resolvents `v^π = R_γ(P_π) r_π` and + `v̂^π = R_γ(P̂_π) r_π` of two state kernels when `P̂_π` deviates from `P_π` by a Bernstein-type + bound on the first `m ≥ log (1 - γ)⁻¹` functions `V 0 = v^π` and + `V (l + 1) = R_γ(P_π) √Var_{P_π}(V l)`. +* The **probabilistic part** is Bernstein's inequality for the empirical state kernel of `π` + (`measureReal_sqrt_variance_add_le_abs_integral_sub_integral_stateKernel_empirical_le`). + The functions `V l` depend on `M` and `π` only, so that Bernstein's inequality applies to each of + them, and a union bound over the `m + 1 ≤ 2 log (e / (1 - γ))` levels and the `S` states gives + the Bernstein-type deviations with `ε = 2 L / N`. + +For a stochastic policy, the empirical state kernel at `s` averages over the actions `a` the means +of the samples of `(s, a)`: its `i`-th summand `∑ a, π(a | s) V (ω i (s, a))` is a mixture of +independent samples, whose variance is at most the variance `Var_{P_π s}(V)` of the mixture law +by Jensen's inequality. [LWCGC24] states the lemma for deterministic policies. + +## Main statements + +* `measureReal_sqrt_variance_add_le_abs_integral_sub_integral_stateKernel_empirical_le`: Bernstein's + inequality for the empirical state kernel of a stationary policy. +* `Learning.MDP.le_streamMeasure_forall_abs_value_empirical_sub_value_le`: model-based policy + evaluation with a generative model [LWCGC24, Lemma 1]. + +## References + +* [G. Li, Y. Wei, Y. Chi, Y. Gu and Y. Chen, *Breaking the sample size barrier in model-based + reinforcement learning with a generative model*, Operations Research, 2024][LWCGC24], Lemmas 1 + and 2 (arXiv:2005.12900) +-/ + +public section + +open MeasureTheory ProbabilityTheory LMLExplore Learning Bandits Finset +open Real (log exp) +open scoped NNReal + +namespace Learning.MDP + +variable {𝓢 𝓐 𝓡 : Type*} {m𝓢 : MeasurableSpace 𝓢} {m𝓐 : MeasurableSpace 𝓐} + {m𝓡 : MeasurableSpace 𝓡} [MeasurableSingletonClass 𝓢] [MeasurableSingletonClass 𝓐] + +/-- **Bernstein's inequality for the empirical state kernel of a stationary policy.** Under the +generative model `streamMeasure M.P` of an MDP `M` with finitely many states and actions, let `π` +be a stationary policy and `f` a function with values in `[l, u]`, `l < u`, `M.P (s, a)`-almost +everywhere for every action `a`. For `N ≠ 0` and every `x`, the integral of `f` against the state +kernel of `π` in the empirical MDP `M.empirical N ω` deviates from its integral against the state +kernel of `π` in `M` by at least `√(2 Var[f; P_π s] x / N) + (u - l) x / (3 N)` with probability at +most `2 exp (-x)`. -/ +theorem measureReal_sqrt_variance_add_le_abs_integral_sub_integral_stateKernel_empirical_le + [Finite 𝓢] [Finite 𝓐] (M : MDP 𝓢 𝓐 𝓡) (π : Kernel 𝓢 𝓐) [IsMarkovKernel π] + {f : 𝓢 → ℝ} {l u : ℝ} (hlu : l < u) {s : 𝓢} + (hfb : ∀ a, ∀ᵐ y ∂(M.P (s, a)), f y ∈ Set.Icc l u) (N : ℕ) [NeZero N] (x : ℝ) : + (streamMeasure M.P).real {ω | √(2 * Var[f; M.stateKernel π s] * x / N) + + (u - l) * x / (3 * N) ≤ + |∫ y, f y ∂(M.stateKernel π s) - ∫ y, f y ∂((M.empirical N ω).stateKernel π s)|} ≤ + 2 * exp (-x) := by + -- The empirical integral is the mean of the `N` independent variables + -- `Z (ω i) = ∑ a, π(a | s) f (ω i (s, a))`, `i < N`, with values in `[l, u]` and mean + -- `∫ f d(P_π s)`. By Jensen's inequality, their variance is at most the variance `Var[f; P_π s]` + -- of the mixture `P_π s = ∑ a, π(a | s) M.P (s, a)`, and the bound is Bernstein's inequality for + -- their mean. + have := Fintype.ofFinite 𝓐 + have hf : Measurable f := .of_discrete + -- the weights `w a = π(a | s)` of the actions + set w : 𝓐 → ℝ := fun a ↦ (π s).real {a} with hw_def + have hw₀ (a : 𝓐) : 0 ≤ w a := measureReal_nonneg + have hw₁ : ∑ a, w a = 1 := by simp [hw_def] + have h_comb {y : 𝓐 → ℝ} (hy : ∀ a, y a ∈ Set.Icc l u) : ∑ a, w a * y a ∈ Set.Icc l u := by + simpa only [smul_eq_mul] using (convex_Icc l u).sum_mem (fun a _ ↦ hw₀ a) hw₁ fun a _ ↦ hy a + -- the mean `m = ∫ f d(P_π s)` and the variable `Z` + set m := ∫ y, f y ∂(M.stateKernel π s) with hm_def + have hm : m = ∑ a, w a * ∫ y, f y ∂(M.P (s, a)) := integral_stateKernel_eq_sum .of_finite + have hm_mem : m ∈ Set.Icc l u := + hm ▸ h_comb fun a ↦ integral_mem_Icc_of_ae_mem_Icc hf.aemeasurable (hfb a) + set Z : (𝓢 × 𝓐 → 𝓢) → ℝ := fun ρ ↦ ∑ a, w a * f (ρ (s, a)) with hZ_def + have hZ : Measurable Z := + Finset.measurable_sum _ fun a _ ↦ (hf.comp (measurable_pi_apply (s, a))).const_mul _ + -- the law of the samples `ρ (s, a)` under `Measure.infinitePi M.P` is `M.P (s, a)` + set μ := Measure.infinitePi M.P with hμ_def + have h_law (a : 𝓐) := hasLaw_eval_infinitePi M.P (s, a) + have h_ae : ∀ᵐ ρ ∂μ, ∀ a, f (ρ (s, a)) ∈ Set.Icc l u := + ae_all_iff.2 fun a ↦ ae_of_ae_map (h_law a).aemeasurable ((h_law a).map_eq ▸ hfb a) + have h_int (a : 𝓐) {g : ℝ → ℝ} (hg : Measurable g) {c d : ℝ} + (hgb : ∀ y ∈ Set.Icc l u, g y ∈ Set.Icc c d) : + Integrable (fun ρ : 𝓢 × 𝓐 → 𝓢 ↦ g (f (ρ (s, a)))) μ := + Integrable.of_mem_Icc c d ((hg.comp hf).comp_aemeasurable (h_law a).aemeasurable) + (h_ae.mono fun ρ hρ ↦ hgb _ (hρ a)) + have hZm : ∫ ρ, Z ρ ∂μ = m := by + have h₁ (a : 𝓐) : ∫ ρ, f (ρ (s, a)) ∂μ = ∫ y, f y ∂(M.P (s, a)) := + (h_law a).integral_comp hf.aestronglyMeasurable + have hi (a : 𝓐) : Integrable (fun ρ : 𝓢 × 𝓐 → 𝓢 ↦ f (ρ (s, a))) μ := + h_int a (g := fun y ↦ y) measurable_id' fun y hy ↦ hy + simp only [hZ_def] + rw [integral_finsetSum _ fun a _ ↦ (hi a).const_mul _] + simp_rw [integral_const_mul, h₁, ← hm] + -- the centred variable is bounded by `u - l` + have hZb : ∀ᵐ ρ ∂μ, |Z ρ - m| ≤ u - l := by + filter_upwards [h_ae] with ρ hρ + have hZρ := h_comb hρ + rw [abs_le] + constructor <;> simp only [hZ_def] <;> linarith [hZρ.1, hZρ.2, hm_mem.1, hm_mem.2] + -- by Jensen's inequality, the variance of `Z` is at most `Var[f; P_π s]` + have hZσ : ∫ ρ, (Z ρ - m) ^ 2 ∂μ ≤ Var[f; M.stateKernel π s] := by + have h_pt (ρ : 𝓢 × 𝓐 → 𝓢) : (Z ρ - m) ^ 2 ≤ ∑ a, w a * (f (ρ (s, a)) - m) ^ 2 := by + have hZρ : Z ρ - m = ∑ a, w a * (f (ρ (s, a)) - m) := by + simp only [hZ_def, mul_sub, sum_sub_distrib, ← sum_mul, hw₁, one_mul] + rw [hZρ] + exact sq_dotProduct_le_dotProduct_sq_of_mem_simplex (mem_simplex_iff.2 ⟨hw₀, hw₁⟩) + fun a ↦ f (ρ (s, a)) - m + have h_int_sq (a : 𝓐) := h_int a (g := fun y ↦ (y - m) ^ 2) (by fun_prop) + fun y hy ↦ sq_sub_mem_Icc_of_mem_Icc hy hm_mem + have h_int_Z : Integrable (fun ρ ↦ (Z ρ - m) ^ 2) μ := + Integrable.of_mem_Icc 0 ((u - l) ^ 2) ((hZ.sub_const m).pow_const 2).aemeasurable + (hZb.mono fun ρ hρ ↦ + ⟨sq_nonneg _, by rw [← sq_abs]; exact pow_le_pow_left₀ (abs_nonneg _) hρ 2⟩) + calc ∫ ρ, (Z ρ - m) ^ 2 ∂μ + ≤ ∫ ρ, ∑ a, w a * (f (ρ (s, a)) - m) ^ 2 ∂μ := + integral_mono h_int_Z (integrable_finsetSum _ fun a _ ↦ (h_int_sq a).const_mul _) h_pt + _ = ∑ a, w a * ∫ y, (f y - m) ^ 2 ∂(M.P (s, a)) := by + rw [integral_finsetSum _ fun a _ ↦ (h_int_sq a).const_mul _] + refine sum_congr rfl fun a _ ↦ ?_ + rw [integral_const_mul] + congr 1 + exact (h_law a).integral_comp (f := fun y ↦ (f y - m) ^ 2) + ((hf.sub_const m).pow_const 2).aestronglyMeasurable + _ = Var[f; M.stateKernel π s] := by + rw [← integral_stateKernel_eq_sum (v := fun y ↦ (f y - m) ^ 2) .of_finite, + variance_eq_integral hf.aemeasurable] + refine le_trans (measureReal_mono fun ω (hω : _ ≤ _) ↦ ?_) + (measureReal_sqrt_add_le_abs_sub_inv_mul_sum_le M.P hZ (sub_pos.2 hlu) hZm hZb hZσ N x) + -- the empirical integral is the mean of the `Z (ω i)`, `i < N` + have h_emp : ∫ y, f y ∂((M.empirical N ω).stateKernel π s) = + (N : ℝ)⁻¹ * ∑ i ∈ range N, Z (ω i) := by + rw [integral_stateKernel_eq_sum .of_finite, mul_sum] + simp_rw [empirical_P, integral_empiricalKernel N ω _ hf.stronglyMeasurable, smul_eq_mul, + hZ_def, mul_sum] + rw [sum_comm] + exact sum_congr rfl fun i _ ↦ sum_congr rfl fun a _ ↦ by ring + rw [Set.mem_ofPred_eq, ← h_emp] + exact hω + +/-- **Model-based policy evaluation with a generative model** [LWCGC24, Lemma 1]. Under the +generative model of an MDP `M` with `S` states, finitely many actions and rewards in `[0, C]`, let +`γ < 1`, `π` a stationary policy, `δ > 0` and `L = log (4 S log (e / (1 - γ)) / δ)`. If +`18 L ≤ (1 - γ) N`, then with probability at least `1 - δ`, the value `v̂^π` of `π` in the empirical +MDP `M.empirical N ω` satisfies `|v̂^π - v^π| ≤ 3 C √(2 L / (N (1 - γ) ^ 3))` at every state. + +[LWCGC24] assumes `N ≥ 32 e ^ 2 L / (1 - γ)`, which implies `18 L ≤ (1 - γ) N`, and proves the +bound `6 √(2 L / (N (1 - γ) ^ 3))` for rewards in `[0, 1]` and deterministic policies; the constant +`3` comes from the constant `1` of +`ProbabilityTheory.Kernel.integral_sqrt_variance_integral_resolvent_le`, which improves the +constant `4` of [LWCGC24, Lemma 11]. -/ +theorem le_streamMeasure_forall_abs_value_empirical_sub_value_le [Fintype 𝓢] [Finite 𝓐] + {M : MDP 𝓢 𝓐 ℝ} {γ : ℝ≥0} + (hγ : γ < 1) {C : ℝ} (hR : ∀ p, ∀ᵐ x ∂(M.R p), x ∈ Set.Icc 0 C) (π : Kernel 𝓢 𝓐) + [IsMarkovKernel π] (N : ℕ) [NeZero N] {δ : ℝ} (hδ : 0 < δ) + (hN : 18 * log (4 * Fintype.card 𝓢 * log (exp 1 / (1 - γ)) / δ) ≤ (1 - γ) * N) : + 1 - δ ≤ (streamMeasure M.P).real {ω | ∀ s, + |(M.empirical N ω).value γ π s - M.value γ π s| ≤ + 3 * C * √(2 * log (4 * Fintype.card 𝓢 * log (exp 1 / (1 - γ)) / δ) / + (N * (1 - γ) ^ 3))} := by + -- We apply the deterministic comparison of resolvents + -- `Kernel.abs_integral_resolvent_sub_integral_resolvent_le_of_abs_sub_le` with `ε = 2 L / N` + -- and `m = ⌈log (1 - γ)⁻¹⌉₊` levels. The functions `V l` of the levels only depend on `M` and + -- `π`, so that Bernstein's inequality for the empirical state kernel of `π`, with a union bound + -- over the `(m + 1) S ≤ 2 S log (e / (1 - γ))` levels and states, gives the deviations it needs. + rcases le_or_gt 1 δ with hδ1 | hδ1 + · exact (sub_nonpos.2 hδ1).trans measureReal_nonneg + have hN₀ : (0 : ℝ) < N := by exact_mod_cast Nat.pos_of_ne_zero (NeZero.ne N) + have hγ₁ : (γ : ℝ) < 1 := hγ + have hu : 0 < 1 - (γ : ℝ) := sub_pos.2 hγ₁ + have hRabs (p : 𝓢 × 𝓐) : ∀ᵐ x ∂(M.R p), |x| ≤ C := + (hR p).mono fun x hx ↦ abs_le.2 ⟨by linarith [hx.1, hx.2], hx.2⟩ + set f := M.policyReward π with hf_def + have hf : Measurable f := .of_discrete + have hf₀ (s : 𝓢) : 0 ≤ f s := by + rw [hf_def, policyReward_def] + exact integral_nonneg fun a ↦ by + rw [meanReward_def] + exact integral_nonneg_of_ae ((hR _).mono fun x hx ↦ hx.1) + have hfC (s : 𝓢) : f s ≤ C := (le_abs_self _).trans (abs_policyReward_le hRabs s) + -- both values are resolvents of the same reward `f = r_π` + set P := M.stateKernel π with hP_def + have hv (ω : ℕ → 𝓢 × 𝓐 → 𝓢) (s : 𝓢) : + (M.empirical N ω).value γ π s - M.value γ π s = + ∫ y, f y ∂(Kernel.resolvent γ ((M.empirical N ω).stateKernel π) s) - + ∫ y, f y ∂(Kernel.resolvent γ P s) := by + have hfe : (M.empirical N ω).policyReward π = f := + funext fun s ↦ by rw [policyReward_def, empirical_R, hf_def, policyReward_def] + rw [value_eq_integral_resolvent hγ hRabs, value_eq_integral_resolvent hγ (by simpa using hRabs), + hfe] + -- for `γ = 0`, both values are the reward + rcases eq_or_ne γ 0 with rfl | hγ₀ + · calc 1 - δ ≤ (streamMeasure M.P).real Set.univ := by simp [hδ.le] + _ ≤ _ := measureReal_mono fun ω _ ↦ by + simp only [Set.mem_ofPred_eq] + intro s + rw [hv, Kernel.resolvent_zero, Kernel.resolvent_zero, Kernel.id_apply, sub_self, + abs_zero] + have hC : 0 ≤ C := (hf₀ s).trans (hfC s) + positivity + have hγ₀' : 0 < γ := pos_iff_ne_zero.2 hγ₀ + -- the confidence level `L`, the accuracy `ε = 2 L / N` and the number `m` of levels + set S : ℝ := (Fintype.card 𝓢 : ℝ) with hS_def + set lg : ℝ := log (exp 1 / (1 - γ)) with hlg_def + have hlg : lg = 1 + log (1 - γ)⁻¹ := by + rw [hlg_def, Real.log_div (Real.exp_pos 1).ne' hu.ne', Real.log_exp, Real.log_inv] + ring + have hy : 0 ≤ log (1 - (γ : ℝ))⁻¹ := + Real.log_nonneg (one_le_inv₀ hu |>.2 (by linarith [γ.coe_nonneg])) + set L := log (4 * S * lg / δ) with hL_def + set ε := 2 * L / N with hε_def + set m := ⌈log (1 - (γ : ℝ))⁻¹⌉₊ with hm_def + have hm : log (1 - (γ : ℝ))⁻¹ ≤ m := Nat.le_ceil _ + have hm' : (m : ℝ) + 1 ≤ 2 * lg := by + have := Nat.ceil_lt_add_one hy + rw [hlg] + linarith + have hS (s : 𝓢) : 1 ≤ S := Nat.one_le_cast.2 (Fintype.card_pos_iff.2 ⟨s⟩) + have hlg1 : 1 ≤ lg := by rw [hlg]; linarith + have hL (s : 𝓢) : 0 ≤ L := Real.log_nonneg (by rw [le_div_iff₀ hδ]; nlinarith [hS s]) + -- the levels `V l` and their bounds `b l` + set V : ℕ → 𝓢 → ℝ := fun l ↦ Nat.rec (fun x ↦ ∫ y, f y ∂(Kernel.resolvent γ P x)) + (fun _ W x ↦ ∫ y, √Var[W; P y] ∂(Kernel.resolvent γ P x)) l with hV_def + have hV₀ (x : 𝓢) : V 0 x = ∫ y, f y ∂(Kernel.resolvent γ P x) := rfl + have hV (l : ℕ) (x : 𝓢) : V (l + 1) x = ∫ y, √Var[V l; P y] ∂(Kernel.resolvent γ P x) := rfl + set b : ℕ → ℝ := fun l ↦ C / ((1 - γ) * (γ * √(1 - γ)) ^ l) with hb_def + have hVb : ∀ l y, V l y ∈ Set.Icc 0 (b l) := + Kernel.mem_Icc_of_eq_integral_resolvent_sqrt_variance hγ₀' hγ hf hf₀ hfC hV₀ hV + -- the bad events: a Bernstein deviation fails at some level `l ≤ m` and state `s` + set bad : Fin (m + 1) × 𝓢 → Set (ℕ → 𝓢 × 𝓐 → 𝓢) := fun i ↦ {ω | + √(ε * Var[V i.1; P i.2]) + ε * b i.1 < + |∫ y, V i.1 y ∂((M.empirical N ω).stateKernel π i.2) - ∫ y, V i.1 y ∂(P i.2)|} + with hbad_def + have h_bad (i : Fin (m + 1) × 𝓢) : + (streamMeasure M.P).real (bad i) ≤ δ / Fintype.card (Fin (m + 1) × 𝓢) := by + have hS := hS i.2 + have hL := hL i.2 + have hpos : 0 < 4 * S * lg / δ := by positivity + have hcard : (Fintype.card (Fin (m + 1) × 𝓢) : ℝ) = (m + 1) * S := by + simp [hS_def, Fintype.card_prod] + have hb₀ : 0 ≤ b i.1 := (hVb i.1 i.2).1.trans (hVb i.1 i.2).2 + rcases hb₀.eq_or_lt with hb0 | hbpos + · -- if `b l = 0`, then `V l = 0` and the bad event is empty + have hV0 : V i.1 = 0 := + funext fun y ↦ le_antisymm ((hVb i.1 y).2.trans_eq hb0.symm) (hVb i.1 y).1 + have h_empty : bad i = ∅ := by + ext ω + simp only [hbad_def, Set.mem_ofPred_eq, hV0, Pi.zero_apply, integral_zero, sub_self, + abs_zero, ← hb0, mul_zero, add_zero, Set.mem_empty_iff_false, iff_false, not_lt] + exact Real.sqrt_nonneg _ + rw [h_empty, measureReal_empty] + positivity + -- Bernstein's inequality at `x = L`, with `2 exp (-L) ≤ δ / ((m + 1) S)` + refine le_trans (measureReal_mono (fun ω (hω : _ < _) ↦ ?_) (measure_ne_top _ _)) + ((measureReal_sqrt_variance_add_le_abs_integral_sub_integral_stateKernel_empirical_le M π + (l := 0) (u := b i.1) (s := i.2) hbpos (fun a ↦ .of_forall (hVb i.1)) N L).trans ?_) + · rw [Set.mem_ofPred_eq, abs_sub_comm] + refine le_trans (add_le_add (le_of_eq ?_) ?_) hω.le + · rw [hε_def] + congr 1 + ring + · have hbL : 0 ≤ b i.1 * L / N := by positivity + calc (b i.1 - 0) * L / (3 * N) = b i.1 * L / N / 3 := by + rw [sub_zero, div_div, mul_comm (3 : ℝ)] + _ ≤ b i.1 * L / N * 2 := by linarith + _ = ε * b i.1 := by rw [hε_def]; ring + · rw [hcard, Real.exp_neg, hL_def, Real.exp_log hpos, inv_div] + calc 2 * (δ / (4 * S * lg)) = δ / (2 * S * lg) := by field_simp; ring + _ ≤ δ / ((m + 1) * S) := + div_le_div_of_nonneg_left hδ.le (by positivity) (by nlinarith) + -- outside the bad events, the deterministic comparison of resolvents applies + refine (one_sub_le_measureReal_compl_iUnion hδ.le h_bad).trans + (measureReal_mono fun ω hω ↦ ?_) + simp only [Set.mem_compl_iff, Set.mem_iUnion, not_exists] at hω + simp only [Set.mem_ofPred_eq] + intro s + have hε₀ : 0 ≤ ε := by have := hL s; positivity + have hε9 : 9 * ε ≤ 1 - γ := by + rw [hε_def, ← mul_div_assoc, div_le_iff₀ hN₀] + linarith + rw [hv] + refine (Kernel.abs_integral_resolvent_sub_integral_resolvent_le_of_abs_sub_le + ((M.empirical N ω).stateKernel π) hγ₀' hγ hf hf₀ hfC hV₀ hV hm hε₀ hε9 + (fun l hl x ↦ ?_) s).trans_eq ?_ + · have h := hω (⟨l, Nat.lt_succ_of_le hl⟩, x) + simp only [hbad_def, Set.mem_ofPred_eq, not_lt] at h + exact h + · rw [hε_def, div_div] + +end Learning.MDP