Fixing rdp2approx delta conversion

autodp/autodp_core.py

@@ -200,7 +200,6 @@ def approx_dp_func(delta1):
         elif type_of_update == 'RDP':
             # function output RDP eps as a function of alpha
             self.RenyiDP = converter.pointwise_minimum(self.RenyiDP, func)
-            self.approx_delta = converter.pointwise_minimum(self.approx_delta, converter.rdp_to_delta(self.RenyiDP))
             if fDP_based_conversion:
 
                 fdp_log, fdp_grad_log = converter.rdp_to_fdp_and_fdp_grad_log(func)
@@ -234,12 +233,16 @@ def approx_dp_func(delta1):
                                                             converter.fdp_fdp_grad_to_approxdp(
                                                                 fdp_log, fdp_grad_log,
                                                                 log_flag=True))
+                self.approx_delta = converter.pointwise_minimum(self.approx_delta,
+                          converter.rdp_to_delta(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
 
                 # self.approxDP = converter.pointwise_minimum(self.approxDP,
                 #                                            converter.fdp_to_approxdp(self.fDP))
             else:
                 self.approxDP = converter.pointwise_minimum(self.approxDP,
                          converter.rdp_to_approxdp(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
+                self.approx_delta = converter.pointwise_minimum(self.approx_delta,
+                          converter.rdp_to_delta(self.RenyiDP, BBGHS_conversion=BBGHS_conversion))
                 self.fDP = converter.pointwise_maximum(self.fDP,
                                                        converter.approxdp_func_to_fdp(
                                                            self.approxDP))

autodp/converter.py

@@ -60,7 +60,7 @@ def approxdp(delta):
     return approxdp
 
 
-def rdp_to_delta(rdp):
+def rdp_to_delta(rdp, BBGHS_conversion=True):
     """
     from RDP to delta with a fixed epsilon
     """
@@ -70,48 +70,64 @@ def approx_delta(eps, naive=False):
         """
         approximate delta as a function of epsilon
         """
-
-        def get_eps(delta):
-            if delta == 0:
-                return rdp(np.inf)
+        def get_delta(alpha):
+            if BBGHS_conversion:
+                # Also by Canonne et al., 2020
+                return np.minimum(
+                    (np.exp((alpha-1)*(rdp(alpha)-eps))/alpha)*((1-1/alpha)**(alpha-1)), 1.0)
             else:
-                def fun(x):  # the input the RDP's \alpha
-                    if x <= 1:
-                        return np.inf
-                    else:
-
-                        if naive:
-                            return np.log(1 / delta) / (x - 1) + rdp(x)
-                        bbghs = np.maximum(rdp(x) + np.log((x-1)/x)
-                                          - (np.log(delta) + np.log(x))/(x-1), 0)
-                        """
-                        The following is for optimal conversion
-                        1/(alpha -1 )log(e^{(alpha-1)*rdp -1}/(alpha*delta) +1 )
-                        """
-                        sign, term_1= utils.stable_log_diff_exp((x-1)*rdp(x),0)
-                        result = utils.stable_logsumexp_two(term_1 - np.log(x)- np.log(delta),0)
-                        return min(result*1.0/(x - 1), bbghs)
-
-                results = minimize_scalar(fun, method='Brent', bracket=(1, 2), bounds=[1, 100000])
-                if results.success:
-                   # print('delta', delta,'eps under rdp', results.fun)
-                    return results.fun
-                else:
-                    return np.inf
+                # Naive conversion from Mironov
+                return np.minimum(np.exp((alpha-1)*(rdp(alpha)-eps)),1.0)
 
-
-        def err(delta):
-            current_eps = get_eps(delta)
-            #print('current delta', delta, 'eps', current_eps)
-            return abs(eps - current_eps)
-
-        results = minimize_scalar(err, method='bounded', bounds=[0, 0.1],options={'xatol':1e-14})
+        results = minimize_scalar(get_delta, method = 'Brent', bracket=(1,2))
         if results.success:
-            #print('results', results.x)
-            return results.x
+            # print('delta', delta,'eps under rdp', results.fun)
+            return results.fun
         else:
             print('not found')
-            return 1
+            return 1.0
+
+        # def get_eps(delta):
+        #     if delta == 0:
+        #         return rdp(np.inf)
+        #     else:
+        #         def fun(x):  # the input the RDP's \alpha
+        #             if x <= 1:
+        #                 return np.inf
+        #             else:
+        #
+        #                 if naive:
+        #                     return np.log(1 / delta) / (x - 1) + rdp(x)
+        #                 bbghs = np.maximum(rdp(x) + np.log((x-1)/x)
+        #                                   - (np.log(delta) + np.log(x))/(x-1), 0)
+        #                 """
+        #                 The following is for optimal conversion
+        #                 1/(alpha -1 )log(e^{(alpha-1)*rdp -1}/(alpha*delta) +1 )
+        #                 """
+        #                 sign, term_1= utils.stable_log_diff_exp((x-1)*rdp(x),0)
+        #                 result = utils.stable_logsumexp_two(term_1 - np.log(x)- np.log(delta),0)
+        #                 return min(result*1.0/(x - 1), bbghs)
+        #
+        #         results = minimize_scalar(fun, method='Brent', bracket=(1, 2))#, bounds=[1, 100000])
+        #         if results.success:
+        #            # print('delta', delta,'eps under rdp', results.fun)
+        #             return results.fun
+        #         else:
+        #             return np.inf
+        #
+        #
+        # def err(delta):
+        #     current_eps = get_eps(delta)
+        #     #print('current delta', delta, 'eps', current_eps)
+        #     return abs(eps - current_eps)
+        #
+        # results = minimize_scalar(err, method='bounded', bounds=[0, 0.1],options={'xatol':1e-14})
+        # if results.success:
+        #     #print('results', results.x)
+        #     return results.x
+        # else:
+        #     print('not found')
+        #     return 1
 
 
     return approx_delta
@@ -146,7 +162,7 @@ def fun(x):  # the input the RDP's \alpha
                     else:
                         return np.log(1 / delta) / (x - 1) + rdp(x)
 
-            results = minimize_scalar(fun, method='Brent', bracket=(1,2), bounds=[1, alpha_max])
+            results = minimize_scalar(fun, method='Brent', bracket=(1,2))#, bounds=[1, alpha_max])
             if results.success:
                 return results.fun
             else:
@@ -244,13 +260,13 @@ def fdp(x):
 
         def fun(alpha):
             if alpha < 0.5:
-                return np.inf
+                return 0
             else:
                 single_fdp = single_rdp_to_fdp(alpha, rdp(alpha))
                 return -single_fdp(x)
 
         # This will use brent to start with 1,2.
-        results = minimize_scalar(fun, bracket=(0.5, 2), bounds=(0.5, alpha_max))
+        results = minimize_scalar(fun, bracket=(0.5, 2))#, bounds=(0.5, alpha_max))
         if results.success:
             return -results.fun
         else:
@@ -527,7 +543,7 @@ def fun(alpha):
                 return log_one_minus_fdp_alpha(logx)
 
         # This will use brent to start with 1,2.
-        results = minimize_scalar(fun, bracket=(0.5, 2), bounds=(0.5, alpha_max))
+        results = minimize_scalar(fun, bracket=(0.5, 2))#, bounds=(0.5, alpha_max))
         if results.success:
             return [results.fun, results.x]
         else:
@@ -720,7 +736,7 @@ def normal_equation_loglogx(loglogx):
             bound1 = np.log(-tmp - tmp**2 / 2 - tmp**3 / 6)
         else:
             bound1 = np.log(1-np.exp(fun1(np.log(1-delta))))
-        #results = minimize_scalar(normal_equation, bounds=[-np.inf,0], bracket=[-1,-2])
+        #results = minimize_scalar(normal_equation, bracket=[-1,-2])
         results = minimize_scalar(normal_equation, method="Bounded", bounds=[bound1,0],
                                   options={'xatol': 1e-10, 'maxiter': 500, 'disp': 0})
         if results.success:

autodp/dp_bank.py

@@ -85,7 +85,7 @@ def fun(x):
     else:
         raise RuntimeError(f"Failed to find epsilon: {results.flag}")
 
-def eps_generalized_gaussian(x, sigma, delta,k, c, c_tilde):
+def get_eps_gaussian_svt(x, sigma, delta, k, c, c_tilde):
     """
     submodule for generalized SVT with Gaussian noise
     we want to partition c into [c/c'] parts, each part using (k choose c')

autodp/mechanism_zoo.py

@@ -431,7 +431,7 @@ def __init__(self, params=None,approxDP_off=False, name='StageWiseMechanism'):
         self.delta0 = 0
 
         if not approxDP_off:  # Direct implementation of approxDP
-            new_approxdp = lambda x: dp_bank.get_generalized_gaussian(params, x)
+            new_approxdp = lambda x: dp_bank.get_eps_gaussian_svt(params, x)
             self.propagate_updates(new_approxdp, 'approxDP_func')
 
 

example/example_amplification_by_sampling.py

@@ -38,9 +38,9 @@
 
 
 # Now let's do subsampling. First we need to use replace-one version of the base mechanisms.
-gm1.replace_one = True
-gm2.replace_one = True
-SVT.replace_one = True
+gm1.neighboring = "replace_one"
+gm2.neighboring = "replace_one"
+SVT.neighboring = "replace_one"
 
 composed_subsampled_mech = compose([subsample(gm1,prob),
                                     subsample(gm2,prob),

setup.py

@@ -18,14 +18,14 @@ def _parse_requirements(path):
 
 setup(
    name='autodp',
-   version='0.2.1b',
+   version='0.2.3',
    description='Automating Differential Privacy Computation',
    license="Apache",
    long_description="The package helps researchers and developers to correctly use advanced methods in differential privacy and obtain provable DP guarantees. The core of the package is an analytical moments accountant that keeps track of Renyi Differential Privacy in analytical forms.",
    author='Yu-Xiang Wang',
    author_email='[email protected]',
    url='https://github.com/yuxiangw/autodp',
-   download_url = 'https://github.com/yuxiangw/autodp/archive/refs/tags/v0.2.1b.tar.gz',   
+   download_url = 'https://github.com/yuxiangw/autodp/archive/refs/tags/v0.2.3.tar.gz',   
    keywords = ['Differential Privacy','Moments Accountant','Renyi Differential Privacy'],
    packages=['autodp'],  #same as name
    install_requires=[_parse_requirements('requirements.txt')], #external packages as dependencies

tutorials/tutorial_PATE_with_autoDP.ipynb

@@ -54,15 +54,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "2.2540846502197414 1e-06\n",
+      "2.2540846502197396 1e-06\n",
       "4.886554117462211 1e-06\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -111,8 +111,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "PATE_eps=2 {'sigma': 22.30476293897311} 1.999999977711552\n",
-      "PATE_eps=0.5 {'sigma': 80.57618519751232} 0.49999999740790046\n"
+      "PATE_eps=2 {'sigma': 22.304762938973138} 1.9999999777115467\n",
+      "PATE_eps=0.5 {'sigma': 80.57618519751239} 0.49999999740789647\n"
      ]
     }
    ],
@@ -177,7 +177,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,

tutorials/tutorial_calibrator.ipynb

@@ -33,8 +33,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "GM {'sigma': 36.304688756269286} 0.10000000492560832\n",
-      "Ana_GM {'sigma': 36.304691899114694} 0.09999999565548537\n"
+      "GM {'sigma': 36.30468875626822} 0.10000000492561108\n",
+      "Ana_GM {'sigma': 36.304691899114694} 0.09999999565548348\n"
      ]
     }
    ],
@@ -80,7 +80,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2555: RuntimeWarning: invalid value encountered in double_scalars\n",
+      "/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py:2884: RuntimeWarning: invalid value encountered in scalar divide\n",
       "  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom\n"
      ]
     }
@@ -307,7 +307,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,

tutorials/tutorial_fdp_of_basic_mechanisms.ipynb

@@ -321,7 +321,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.16"
   }
  },
  "nbformat": 4,