Switched Renyi div to ENNReal

SampCert/DiffPrivacy/ConcentratedBound.lean

@@ -75,12 +75,14 @@ theorem SG_Renyi_simplify {σ : ℝ} (h : σ ≠ 0) (μ ν : ℤ) (α : ℝ) :
     left
     rw [← mul_assoc]
 
+noncomputable def RenyiDivergence' (p q : T → ℝ) (α : ℝ) : ℝ :=
+  (1 / (α - 1)) * Real.log (∑' x : T, (p x)^α  * (q x)^(1 - α))
 
-theorem RenyiDivergenceBound {σ : ℝ} (h : σ ≠ 0) (μ : ℤ) (α : ℝ) (h' : α > 1) :
-  RenyiDivergence (fun (x : ℤ) => (gauss_term_ℝ σ μ) x / ∑' x : ℤ, (gauss_term_ℝ σ μ) x)
+theorem RenyiDivergenceBound {σ : ℝ} (h : σ ≠ 0) (μ : ℤ) (α : ℝ) (h' : 1 < α) :
+  RenyiDivergence' (fun (x : ℤ) => (gauss_term_ℝ σ μ) x / ∑' x : ℤ, (gauss_term_ℝ σ μ) x)
                   (fun (x : ℤ) => (gauss_term_ℝ σ (0 : ℤ)) x / ∑' x : ℤ, (gauss_term_ℝ σ (0 : ℤ)) x)
                   α ≤ α * (μ^2 / (2 * σ^2)) := by
-  unfold RenyiDivergence
+  unfold RenyiDivergence'
   have A : 0 < 1 / (α - 1) := by
     simp [h']
   rw [← le_div_iff' A]
@@ -300,9 +302,9 @@ theorem sg_mul_simplify (ss : ℝ) (x μ ν : ℤ) :
   rw [← neg_add]
 
 theorem SG_Renyi_shift {σ : ℝ} (h : σ ≠ 0) (α : ℝ) (μ ν τ : ℤ) :
-  RenyiDivergence (fun (x : ℤ) => (gauss_term_ℝ σ μ) x / ∑' x : ℤ, (gauss_term_ℝ σ μ) x) (fun (x : ℤ) => (gauss_term_ℝ σ ν) x / ∑' x : ℤ, (gauss_term_ℝ σ ν) x) α
-    = RenyiDivergence (fun (x : ℤ) => (gauss_term_ℝ σ ((μ + τ) : ℤ)) x / ∑' x : ℤ, (gauss_term_ℝ σ ((μ + τ) : ℤ)) x) (fun (x : ℤ) => (gauss_term_ℝ σ ((ν + τ) : ℤ)) x / ∑' x : ℤ, (gauss_term_ℝ σ ((ν + τ) : ℤ)) x) α := by
-  unfold RenyiDivergence
+  RenyiDivergence' (fun (x : ℤ) => (gauss_term_ℝ σ μ) x / ∑' x : ℤ, (gauss_term_ℝ σ μ) x) (fun (x : ℤ) => (gauss_term_ℝ σ ν) x / ∑' x : ℤ, (gauss_term_ℝ σ ν) x) α
+    = RenyiDivergence' (fun (x : ℤ) => (gauss_term_ℝ σ ((μ + τ) : ℤ)) x / ∑' x : ℤ, (gauss_term_ℝ σ ((μ + τ) : ℤ)) x) (fun (x : ℤ) => (gauss_term_ℝ σ ((ν + τ) : ℤ)) x / ∑' x : ℤ, (gauss_term_ℝ σ ((ν + τ) : ℤ)) x) α := by
+  unfold RenyiDivergence'
   congr 2
   conv =>
     left
@@ -410,11 +412,144 @@ theorem SG_Renyi_shift {σ : ℝ} (h : σ ≠ 0) (α : ℝ) (μ ν τ : ℤ) :
     apply Summable.mul_right
     apply summable_gauss_term' h
 
-theorem RenyiDivergenceBound' {σ α : ℝ} (h : σ ≠ 0) (h' : 1 < α) (μ ν : ℤ)   :
-  RenyiDivergence (fun (x : ℤ) => discrete_gaussian σ μ x)
+theorem RenyiDivergenceBound_pre {σ α : ℝ} (h : σ ≠ 0) (h' : 1 < α) (μ ν : ℤ)   :
+  RenyiDivergence' (fun (x : ℤ) => discrete_gaussian σ μ x)
                   (fun (x : ℤ) => discrete_gaussian σ ν x)
                   α ≤ α * (((μ - ν) : ℤ)^2 / (2 * σ^2)) := by
   unfold discrete_gaussian
   rw [SG_Renyi_shift h α μ ν (-ν)]
   rw [add_right_neg]
   apply  RenyiDivergenceBound h (μ + -ν) α h'
+
+theorem RenyiSumSG_nonneg {σ α : ℝ} (h : σ ≠ 0) (μ ν n : ℤ) :
+  0 ≤ discrete_gaussian σ μ n ^ α * discrete_gaussian σ ν n ^ (1 - α) := by
+  have A := discrete_gaussian_nonneg h μ n
+  have B := discrete_gaussian_nonneg h ν n
+  rw [mul_nonneg_iff]
+  left
+  constructor
+  . apply Real.rpow_nonneg A
+  . apply Real.rpow_nonneg B
+
+theorem SummableRenyiGauss {σ : ℝ} (h : σ ≠ 0) (μ ν : ℤ) (α : ℝ) :
+  Summable fun (x : ℤ) => discrete_gaussian σ μ x ^ α * discrete_gaussian σ ν x ^ (1 - α) := by
+  simp [discrete_gaussian]
+  have B : ∀ μ : ℤ, ∀ x : ℝ, 0 ≤ (gauss_term_ℝ σ μ) x := by
+    intro μ x
+    unfold gauss_term_ℝ
+    apply exp_nonneg
+  have C : ∀ μ : ℤ, 0 ≤ (∑' (x : ℤ), (gauss_term_ℝ σ μ) x)⁻¹ := by
+    intro μ
+    rw [inv_nonneg]
+    apply le_of_lt
+    apply sum_gauss_term_pos h
+  conv =>
+    right
+    intro x
+    rw [division_def]
+    rw [division_def]
+    rw [mul_rpow (B μ x) (C μ)]
+    rw [mul_rpow (B ν x) (C ν)]
+  conv =>
+    right
+    intro x
+    rw [mul_assoc]
+    right
+    rw [← mul_assoc]
+    left
+    rw [mul_comm]
+  conv =>
+    right
+    intro x
+    ring_nf
+  apply Summable.mul_right
+  apply Summable.mul_right
+  unfold gauss_term_ℝ
+  conv =>
+    right
+    intro x
+    rw [← Real.exp_mul]
+    rw [← Real.exp_mul]
+    rw [← exp_add]
+    rw [← mul_div_right_comm]
+    rw [← mul_div_right_comm]
+    rw [div_add_div_same]
+    rw [mul_sub_left_distrib]
+    rw [mul_one]
+    right
+    left
+    ring_nf
+
+  have X : ∀ x : ℤ,
+    ↑x * ↑μ * α * 2 - ↑x * α * ↑ν * 2 + ↑x * ↑ν * 2 + (-↑x ^ 2 - ↑μ ^ 2 * α) + (α * ↑ν ^ 2 - ↑ν ^ 2)
+           = - ((↑x ^ 2 - 2 * x * (↑μ * α - α * ↑ν + ↑ν)) + (↑μ ^ 2 * α - α * ↑ν ^ 2 + ↑ν ^ 2)) := by
+        intro x
+        ring_nf
+  conv =>
+    right
+    intro x
+    rw [X]
+  clear X
+
+  have X : (↑μ ^ 2 * α - α * ↑ν ^ 2 + ↑ν ^ 2)
+    = (↑μ * α - α * ↑ν + ↑ν)^2 + (- ↑μ * α * ↑ν * 2 + ↑μ * α ^ 2 * ↑ν * 2 -
+        ↑μ ^ 2 * α ^ 2 + α * ↑ν ^ 2 - α ^ 2 * ↑ν ^ 2 + α * ↑μ ^ 2) := by
+    ring_nf
+  conv =>
+    right
+    intro x
+    rw [X]
+    rw [← add_assoc]
+  clear X
+
+  have X : ∀ x : ℤ, (x - (↑μ * α - α * ↑ν + ↑ν))^2 = ↑x ^ 2 - 2 * ↑x * (↑μ * α - α * ↑ν + ↑ν) + (↑μ * α - α * ↑ν + ↑ν) ^ 2 := by
+    intro x
+    ring_nf
+  conv =>
+    right
+    intro x
+    rw [← X]
+    rw [neg_add]
+    rw [← div_add_div_same]
+    rw [exp_add]
+  clear X
+  apply Summable.mul_right
+  apply summable_gauss_term' h
+
+theorem RenyiDivergenceBound' {σ α : ℝ} (h : σ ≠ 0) (h' : 1 < α) (μ ν : ℤ)   :
+  RenyiDivergence (fun (x : ℤ) => ENNReal.ofReal (discrete_gaussian σ μ x))
+                  (fun (x : ℤ) => ENNReal.ofReal (discrete_gaussian σ ν x))
+                  α ≤ α * (((μ - ν) : ℤ)^2 / (2 * σ^2)) := by
+  have A : RenyiDivergence (fun (x : ℤ) => ENNReal.ofReal (discrete_gaussian σ μ x))
+                  (fun (x : ℤ) => ENNReal.ofReal (discrete_gaussian σ ν x))
+                  α  = RenyiDivergence' (fun (x : ℤ) => discrete_gaussian σ μ x)
+                  (fun (x : ℤ) => discrete_gaussian σ ν x)
+                  α  := by
+    unfold RenyiDivergence
+    unfold RenyiDivergence'
+    congr
+    simp
+    have A₁ : ∀ x : ℤ, 0 ≤ discrete_gaussian σ μ x ^ α := by
+      intro x
+      apply Real.rpow_nonneg
+      apply discrete_gaussian_nonneg h μ x
+    conv =>
+      left
+      right
+      right
+      intro x
+      rw [ENNReal.ofReal_rpow_of_pos (discrete_gaussian_pos h μ x)]
+      rw [ENNReal.ofReal_rpow_of_pos (discrete_gaussian_pos h ν x)]
+      rw [← ENNReal.ofReal_mul (A₁ x)]
+    rw [← ENNReal.ofReal_tsum_of_nonneg]
+    . simp
+      apply tsum_nonneg
+      intro i
+      apply RenyiSumSG_nonneg h
+    . apply RenyiSumSG_nonneg h
+    . apply SummableRenyiGauss h
+
+
+
+  rw [A]
+  apply RenyiDivergenceBound_pre h h'

SampCert/DiffPrivacy/DP.lean

@@ -14,14 +14,11 @@ import SampCert.DiffPrivacy.Sensitivity
 
 noncomputable section
 
-open Classical Nat Int Real
+open Classical Nat Int Real ENNReal
 
 def DP (q : List T → SLang U) (ε : ℝ) : Prop :=
   ∀ α : ℝ, 1 < α → ∀ l₁ l₂ : List T, Neighbour l₁ l₂ →
-  RenyiDivergence
-    (fun x : U => (q l₁ x).toReal)
-    (fun x : U => (q l₂ x).toReal) α
-    ≤ (1/2) * ε ^ 2 * α
+  RenyiDivergence (q l₁) (q l₂) α ≤ (1/2) * ε ^ 2 * α
 
 namespace SLang
 
@@ -76,7 +73,7 @@ theorem NoisedQueryDP (query : List T → ℤ) (Δ ε₁ ε₂ : ℕ+) (bounded_
       trivial
     . apply sq_nonneg
   . rw [pow_two]
-    rw [mul_pos_iff]
+    rw [_root_.mul_pos_iff]
     left
     simp
 
@@ -93,7 +90,7 @@ theorem DPCompose {nq1 nq2 : List T → SLang ℤ} {ε₁ ε₂ ε₃ ε₄ : 
   replace h1 := h1 α h3 l₁ l₂ h4
   replace h2 := h2 α h3 l₁ l₂ h4
   simp [RenyiDivergence] at h1 h2
-  rw [tsum_prod']
+  rw [tsum_prod' ENNReal.summable (fun b : ℤ => ENNReal.summable)]
   . simp
     have A : ∀ b c : ℤ, ∀ l : List T, (∑' (a : ℤ), nq1 l a * ∑' (a_1 : ℤ), if b = a ∧ c = a_1 then nq2 l a_1 else 0) = nq1 l b * nq2 l c := by
       sorry
@@ -102,15 +99,14 @@ theorem DPCompose {nq1 nq2 : List T → SLang ℤ} {ε₁ ε₂ ε₃ ε₄ : 
       right
       right
       right
+      right
       intro b
       right
       intro c
       rw [A]
       rw [A]
-      rw [ENNReal.toReal_mul]
-      rw [ENNReal.toReal_mul]
-      rw [Real.mul_rpow ENNReal.toReal_nonneg ENNReal.toReal_nonneg]
-      rw [Real.mul_rpow ENNReal.toReal_nonneg ENNReal.toReal_nonneg]
+      rw [ENNReal.mul_rpow_of_ne_zero sorry sorry]
+      rw [ENNReal.mul_rpow_of_ne_zero sorry sorry]
       rw [mul_assoc]
       right
       rw [mul_comm]
@@ -122,38 +118,35 @@ theorem DPCompose {nq1 nq2 : List T → SLang ℤ} {ε₁ ε₂ ε₃ ε₄ : 
       right
       right
       right
+      right
       intro b
       right
       intro c
       rw [← mul_assoc]
-    have B : Summable fun c => (nq2 l₁ c).toReal ^ α * (nq2 l₂ c).toReal ^ (1 - α) := sorry
-    have C : Summable fun b => (nq1 l₁ b).toReal ^ α * (nq1 l₂ b).toReal ^ (1 - α) := sorry
     conv =>
       left
       right
       right
       right
+      right
       intro b
-      rw [@Summable.tsum_mul_left ℤ ℝ _ _ _ (fun c => (nq2 l₁ c).toReal ^ α * (nq2 l₂ c).toReal ^ (1 - α)) _ ((nq1 l₁ b).toReal ^ α * (nq1 l₂ b).toReal ^ (1 - α)) B]
-    rw [Summable.tsum_mul_right ]
-    . rw [Real.log_mul]
-      . rw [mul_add]
-        have D := _root_.add_le_add h1 h2
-        apply le_trans D
-        clear h1 h2 A B C D
-        rw [← add_mul]
-        rw [mul_le_mul_iff_of_pos_right]
-        . rw [← mul_add]
-          rw [mul_le_mul_iff_of_pos_left]
-          . ring_nf
-            simp
-          . sorry -- easy
+      rw [ENNReal.tsum_mul_left]
+    rw [ENNReal.tsum_mul_right]
+    rw [ENNReal.toReal_mul]
+    rw [Real.log_mul]
+    . rw [mul_add]
+      have D := _root_.add_le_add h1 h2
+      apply le_trans D
+      rw [← add_mul]
+      rw [mul_le_mul_iff_of_pos_right]
+      . rw [← mul_add]
+        rw [mul_le_mul_iff_of_pos_left]
+        . ring_nf
+          simp
         . sorry -- easy
-      . sorry -- ≠ 0 , not obvious
-      . sorry -- ≠ 0
-    . apply C
-  . sorry
-  . sorry
+      . sorry -- easy
+    . sorry -- ≠ 0 , not obvious
+    . sorry -- ≠ 0
 
 def PostProcess (nq : List T → SLang U) (pp : U → ℤ) (l : List T) : SLang ℤ := do
   let A ← nq l

SampCert/DiffPrivacy/RenyiDivergence.lean

@@ -6,11 +6,11 @@ Authors: Jean-Baptiste Tristan
 
 import Mathlib.Topology.Algebra.InfiniteSum.Ring
 import Mathlib.NumberTheory.ModularForms.JacobiTheta.TwoVariable
---import Mathlib.Analysis.SpecialFunctions.Pow.NNReal
+import Mathlib.Analysis.SpecialFunctions.Pow.NNReal
 
 open Real ENNReal
 
 variable {T : Type}
 
-noncomputable def RenyiDivergence (p q : T → ℝ) (α : ℝ) : ℝ :=
-  (1 / (α - 1)) * Real.log (∑' x : T, (p x)^α  * (q x)^(1 - α))
+noncomputable def RenyiDivergence (p q : T → ENNReal) (α : ℝ) : ℝ :=
+  (1 / (α - 1)) * Real.log (∑' x : T, (p x)^α  * (q x)^(1 - α)).toReal

SampCert/Samplers/GaussianGen/Properties.lean

@@ -48,7 +48,7 @@ theorem DiscreteGaussianGenSample_apply (num : PNat) (den : PNat) (μ x : ℤ) :
   . rw [SG_periodic' A]
 
 theorem DiscreteGaussianGenSampleZeroConcentrated {α : ℝ} (h : 1 < α) (num : PNat) (den : PNat) (μ ν : ℤ) :
-  RenyiDivergence (fun x : ℤ => ((DiscreteGaussianGenSample num den μ) x).toReal) (fun x : ℤ => ((DiscreteGaussianGenSample num den ν) x).toReal) α ≤
+  RenyiDivergence ((DiscreteGaussianGenSample num den μ)) (DiscreteGaussianGenSample num den ν) α ≤
   α * (((μ - ν) : ℤ)^2 / (2 * ((num : ℝ) / (den : ℝ))^2)) := by
   have A : (num : ℝ) / (den : ℝ) ≠ 0 := by
     simp only [ne_eq, div_eq_zero_iff, cast_eq_zero, PNat.ne_zero, or_self, not_false_eq_true]
@@ -57,10 +57,8 @@ theorem DiscreteGaussianGenSampleZeroConcentrated {α : ℝ} (h : 1 < α) (num :
     congr
     . intro x
       rw [DiscreteGaussianGenSample_apply]
-      rw [ENNReal.toReal_ofReal (discrete_gaussian_nonneg A μ x)]
     . intro x
       rw [DiscreteGaussianGenSample_apply]
-      rw [ENNReal.toReal_ofReal (discrete_gaussian_nonneg A ν x)]
     . skip
   apply RenyiDivergenceBound' A h