Fix rdp

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

@@ -21,8 +21,8 @@ def get_eps_rdp(func, delta):
     :param delta:
     :return: The corresponding epsilon
     """
-    assert(delta >= 0)
-    acct = rdp_acct.anaRDPacct(m=10,m_max=10)
+    assert (delta >= 0)
+    acct = rdp_acct.anaRDPacct(m=10, m_max=10)
     acct.compose_mechanism(func)
     return acct.get_eps(delta)
 
@@ -34,38 +34,38 @@ def get_eps_rdp_subsampled(func, delta, prob):
     :param delta:
     :return: The corresponding epsilon
     """
-    assert(delta >= 0)
-    assert(prob >=0)
-    if prob==0:
+    assert (delta >= 0)
+    assert (prob >= 0)
+    if prob == 0:
         return 0
     elif prob == 1:
-        return get_eps_rdp(func,delta)
+        return get_eps_rdp(func, delta)
     else:
         acct = rdp_acct.anaRDPacct()
-        acct.compose_subsampled_mechanism(func,prob)
+        acct.compose_subsampled_mechanism(func, prob)
         return acct.get_eps(delta)
 
 
 # Get the eps and delta for a single Gaussian mechanism
 def get_eps_gaussian(sigma, delta):
     """ This function calculates the eps for Gaussian Mech given sigma and delta"""
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_gaussian({'sigma':sigma},x)
-    return get_eps_rdp(func,delta)
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_gaussian({'sigma': sigma}, x)
+    return get_eps_rdp(func, delta)
 
 
-def get_logdelta_ana_gaussian(sigma,eps):
+def get_logdelta_ana_gaussian(sigma, eps):
     """ This function calculates the delta parameter for analytical gaussian mechanism given eps"""
-    assert(eps>=0)
+    assert (eps >= 0)
     s, mag = utils.stable_log_diff_exp(norm.logcdf(0.5 / sigma - eps * sigma),
-                                       eps + norm.logcdf(-0.5/sigma - eps * sigma))
+                                       eps + norm.logcdf(-0.5 / sigma - eps * sigma))
     return mag
 
 
 def get_eps_ana_gaussian(sigma, delta):
     """ This function calculates the gaussian mechanism given sigma and delta using analytical GM"""
     # Basically inverting the above function by solving a nonlinear equation
-    assert(delta >=0 and delta <=1)
+    assert (delta >= 0 and delta <= 1)
 
     if delta == 0:
         return np.inf
@@ -78,71 +78,94 @@ def fun(x):
         else:
             return get_logdelta_ana_gaussian(sigma, x) - np.log(delta)
 
-    eps_upperbound = 1/2/sigma**2+1/sigma*np.sqrt(2*np.log(1/delta))
-    results = root_scalar(fun,bracket=[0, eps_upperbound])
+    eps_upperbound = 1 / 2 / sigma ** 2 + 1 / sigma * np.sqrt(2 * np.log(1 / delta))
+    results = root_scalar(fun, bracket=[0, eps_upperbound])
     if results.converged:
         return results.root
     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):
+def get_eps_gaussian_svt(params, delta):
     """
-    submodule for generalized SVT with Gaussian noise
+    Implements Theorem 12 (stage-wise length-capped Gaussian SVT) in Zhu et.al., NeurIPS-20.
+
     we want to partition c into [c/c'] parts, each part using (k choose c')
     need to check whether (k choose c') > log(1/delta')
-    k is the maximam number of queries to answer for each chunk
-    x is log delta for each chunk, it needs to be negative
+    k is the maximum number of queries to answer for each chunk
+    x is delta for each chunk, it needs to be negative
     :param x:
     :param sigma:
     :param delta:
     :return:
     """
-    acct = dp_acct.DP_acct()
-    per_delta = np.exp(x) # per_delta for each c' chunk
-    coeff = comb(k,c_tilde)
-    assert per_delta < 1.0/(coeff)
-    #compute the eps per step with 1/(sigma_1**2) + sqrt(2/simga_1**2 *(log k + log(1/epr_delta)))
-    # compose eps for each chunk
-    while c:
-        if c>= c_tilde:
-            c = c - c_tilde
-        else:
-            # the remaining part c // c_tilde
-            c = 0
-            c_tilde = c
-            coeff = comb(k, c_tilde)
-        per_eps = (1.0+c_tilde) / (2*sigma ** 2) + np.sqrt((1.0 +c_tilde) / (2*sigma ** 2) * (np.log(coeff) - x))
-        acct.update_DPlosses(per_eps, per_delta)
-
-    compose_eps = acct.get_eps(delta)
-    return compose_eps
-
-
-def get_eps_laplace(b,delta):
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_laplace({'b':b},x)
-    return get_eps_rdp(func,delta)
+    # {'sigma': params['sigma'], 'k': params['k'], 'c': params['c']}
+    sigma = params['sigma']
+    k = params.get('k', 1)
+    c = params.get('c', 1)
+    c_tilde = params.get('c_tilde', int(np.sqrt(c)))  # default c_tilde is sqrt(c)
+    coeff = comb(k, c_tilde)
+
+    # the function below returns the composed epsilon if each chunk is assigned a e^x budget of delta (the delta'
+    # used in Theorem 12).
+    def search_per_delta_each_chunk(x, sigma, c, c_tilde, k, coeff):
+        acct = dp_acct.DP_acct()
+        per_delta = np.exp(x)  # per_delta for each c' chunk
+        assert per_delta < 1.0 / (coeff)
+        # compute the eps per step with 1/(sigma_1**2) + sqrt(2/simga_1**2 *(log k + log(1/epr_delta)))
+        # compose eps for each chunk
+        while c:
+            if c >= c_tilde:
+                c = c - c_tilde
+            else:
+                # the remaining part c // c_tilde
+                c = 0
+                c_tilde = c
+                coeff = comb(k, c_tilde)
+            per_eps = (1.0 + c_tilde) / (2 * sigma ** 2) + np.sqrt(
+                (1.0 + c_tilde) / (2 * sigma ** 2) * (np.log(coeff) - x))
+            acct.update_DPlosses(per_eps, per_delta)
+        compose_eps = acct.get_eps(delta)
+        return compose_eps  # composed epsilon over c/c' chunks.
+
+    fun = lambda x: search_per_delta_each_chunk(x, sigma, c, c_tilde, k, coeff)
+    # the bound of per chunk delta (delta').
+    const = 10
+    right_bound = min(delta / (c_tilde + 1), 1.0 / coeff)
+    left_bound =delta / (c_tilde * 10)
+    # ensures the left bound is smaller than the right bound
+    while 10 * left_bound > right_bound:
+        const = const * 10
+        left_bound = left_bound * 0.1
+    results = minimize_scalar(fun, method='Bounded', bounds=[np.log(left_bound), np.log(right_bound)], options={'disp': False})
+    return results.fun
+
+
+def get_eps_laplace(b, delta):
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_laplace({'b': b}, x)
+    return get_eps_rdp(func, delta)
 
 
-def get_eps_randresp(p,delta):
-    assert(delta >= 0)
-    func = lambda x: rdp_bank.RDP_randresponse({'p':p},x)
+def get_eps_randresp(p, delta):
+    assert (delta >= 0)
+    func = lambda x: rdp_bank.RDP_randresponse({'p': p}, x)
     return get_eps_rdp(func, delta)
 
 
+def get_eps_randresp_optimal(p, delta):
+    assert (delta >= 0)
+    if p < 0.5:
+        p = 1 - p
 
-def get_eps_randresp_optimal(p,delta):
-    assert(delta >= 0)
-    if p<0.5:
-        p = 1-p
-
-    if p==0 or p==1:
+    if p == 0 or p == 1:
         return np.inf
 
-    eps_max = np.log(p) - np.log(1-p)
+    eps_max = np.log(p) - np.log(1 - p)
     if delta == 0:
         return eps_max
-    elif delta >= 2*p - 1:
+    elif delta >= 2 * p - 1:
         return 0.0
     else:
-        return np.log(p-delta) - np.log(1 - p)
+        return np.log(p - delta) - np.log(1 - p)