import math
import numpy as np
from scipy import special
from scipy.special import logsumexp
from .utils import _log_add, _log_coef_i, _log_e, _log_s, _log_sub, _log_t
[docs]def eps_gaussian(alpha, params):
if alpha == math.inf:
return min(
4 * (np.exp(eps_gaussian(2, params) - 1)),
2 * np.exp(eps_gaussian(2, params)),
)
return alpha / (2 * (params["sigma"] ** 2))
[docs]def eps_laplace(alpha, params):
if alpha <= 1:
return 1 / params["b"] + np.exp(-1 / params["b"]) - 1
elif alpha == math.inf:
return 1 / params["b"]
return (1 / (alpha - 1)) * logsumexp(
[
np.log(alpha / (2 * alpha - 1)) + (alpha - 1) / params["b"],
np.log((alpha - 1) / (2 * alpha - 1)) + (-alpha) / params["b"],
],
b=[1, 1],
)
[docs]def eps_randresp(alpha, params):
if params["p"] == 1 or params["p"] == 0:
return math.inf
if alpha <= 1:
return (2 * params["p"] - 1) * np.log(params["p"] / (1 - params["p"]))
elif alpha == math.inf:
return np.abs(np.log((1.0 * params["p"] / (1 - params["p"]))))
return (1 / (alpha - 1)) * logsumexp(
[
alpha * np.log(params["p"]) + (1 - alpha) * np.log(1 - params["p"]),
alpha * np.log(1 - params["p"]) + (1 - alpha) * params["p"],
],
b=[1, 1],
)
[docs]def calc_upperbound_of_rdp_with_Sampled_Gaussian_Mechanism(
alpha, params, sampling_rate, _eps
):
"""Compute log(A_alpha) for any positive finite alpha."""
if float(alpha).is_integer():
return calc_upperbound_of_rdp_with_Sampled_Gaussian_Mechanism_int(
int(alpha), params, sampling_rate
)
else:
return calc_upperbound_of_rdp_with_Sampled_Gaussian_Mechanism_float(
alpha, params, sampling_rate
)
[docs]def calc_upperbound_of_rdp_with_Sampled_Gaussian_Mechanism_int(
alpha, params, sampling_rate
):
"""Renyi Differential Privacy of the Sampled Gaussian Mechanism
3.3 Numerically Stable Computatio
"""
log_a = -np.inf
for i in range(alpha + 1):
log_coef_i = (
np.log(special.binom(alpha, i))
+ i * np.log(sampling_rate)
+ (alpha - i) * np.log(1 - sampling_rate)
)
s = log_coef_i + (i * i - i) / (2 * (params["sigma"] ** 2))
log_a = _log_add(log_a, s)
return float(log_a) / (alpha - 1)
[docs]def calc_upperbound_of_rdp_with_Sampled_Gaussian_Mechanism_float(
alpha, params, sampling_rate
):
"""Compute log(A_alpha) for fractional alpha. 0 < q < 1."""
# The two parts of A_alpha, integrals over (-inf,z0] and [z0, +inf), are
# initialized to 0 in the log space:
log_a0, log_a1 = -np.inf, -np.inf
i = 0
z0 = params["sigma"] ** 2 * np.log(1 / sampling_rate - 1) + 0.5
while True: # do ... until loop
coef = special.binom(alpha, i)
log_coef = np.log(abs(coef))
j = alpha - i
log_t0 = _log_t(log_coef, i, j, sampling_rate)
log_t1 = _log_t(log_coef, j, i, sampling_rate)
log_e0 = _log_e(i, z0, params["sigma"])
log_e1 = _log_e(z0, j, params["sigma"])
log_s0 = _log_s(i, log_t0, log_e0, params["sigma"])
log_s1 = _log_s(j, log_t1, log_e1, params["sigma"])
if coef > 0:
log_a0 = _log_add(log_a0, log_s0)
log_a1 = _log_add(log_a1, log_s1)
else:
log_a0 = _log_sub(log_a0, log_s0)
log_a1 = _log_sub(log_a1, log_s1)
i += 1
if max(log_s0, log_s1) < -30:
break
return _log_add(log_a0, log_a1) / (alpha - 1)
[docs]def calc_upperbound_of_rdp_with_theorem27_of_wang_2019(
alpha, params, sampling_rate, _eps
):
def B(el):
res = 0
for i in range(el + 1):
res += (
((-1) ** (i)) * special.binom(el, i) * np.exp((i - 1) * _eps(i, params))
)
return abs(res)
"""
def logAbsB(el):
ts = []
ss = []
for i in range(el + 1):
ts.append(np.log(special.binom(el, i)) + (i - 1) * _eps(i, params))
ss.append((-1) ** (i + el % 2))
return logsumexp(ts, b=ss)
"""
terms = []
signs = []
terms.append(np.log(1))
signs.append(1)
second = (
2 * np.log(sampling_rate)
+ np.log(special.binom(alpha, 2))
+ min(
np.log(4) + logsumexp([_eps(2, params), np.log(1)], b=[1, -1]),
_eps(2, params)
+ min(
np.log(2),
2 * logsumexp([_eps(math.inf, params), np.log(1)], b=[1, -1]),
),
)
)
terms.append(second)
signs.append(1)
for j in range(3, alpha + 1):
terms.append(
np.log(4)
+ j * np.log(sampling_rate)
+ np.log(special.binom(alpha, j))
+ (1 / 2)
* (np.log(B(2 * math.floor(j / 2))) + np.log(B(2 * math.ceil(j / 2))))
)
signs.append(1)
return (1 / (alpha - 1)) * logsumexp(terms, b=signs)
[docs]def calc_first_term_of_general_upper_bound_of_rdp(alpha, sampling_rate):
return ((alpha - 1) * np.log(1 - sampling_rate)) + np.log(
alpha * sampling_rate - sampling_rate + 1
)
[docs]def calc_general_upperbound_of_rdp_with_theorem5_of_zhu_2019(
alpha, params, sampling_rate, _eps
):
terms = []
signs = []
first = calc_first_term_of_general_upper_bound_of_rdp(alpha, sampling_rate)
terms.append(first)
signs.append(1)
if alpha >= 2:
second = (
np.log(special.binom(alpha, 2))
+ (2 * np.log(sampling_rate))
+ ((alpha - 2) * np.log(1 - sampling_rate))
+ _eps(2, params)
)
terms.append(second)
signs.append(1)
if alpha >= 3:
for el in range(3, alpha + 1):
third = np.log(3) + (
_log_coef_i(alpha, el, sampling_rate) + ((el - 1) * _eps(el, params))
)
terms.append(third)
signs.append(1)
return (1 / (alpha - 1)) * logsumexp(terms, b=signs)
[docs]def calc_tightupperbound_lowerbound_of_rdp_with_theorem6and8_of_zhu_2019(
alpha, params, sampling_rate, _eps
):
terms = []
signs = []
first = calc_first_term_of_general_upper_bound_of_rdp(alpha, sampling_rate)
terms.append(first)
signs.append(1)
if alpha >= 2:
for el in range(2, alpha + 1):
second = _log_coef_i(alpha, el, sampling_rate) + (
(el - 1) * _eps(el, params)
)
terms.append(second)
signs.append(1)
return (1 / (alpha - 1)) * logsumexp(terms, b=signs)
[docs]def calc_tightupperbound_lowerbound_of_rdp_with_theorem6and8_of_zhu_2019_with_tau_estimation(
alpha, params, sampling_rate, _eps, tau=10
):
terms = []
signs = []
eps_alpha_minus_tau = _eps(alpha - tau, params)
first = alpha * np.log(1 - sampling_rate) + logsumexp(
[np.log(1), -eps_alpha_minus_tau], b=[1, -1]
)
terms.append(first)
signs.append(1)
second = -eps_alpha_minus_tau + alpha * logsumexp(
[
np.log(1),
np.log(sampling_rate),
np.log(sampling_rate) + eps_alpha_minus_tau,
],
b=[1, -1, 1],
)
terms.append(second)
signs.append(1)
for el in range(2, tau):
third = _log_coef_i(alpha, el, sampling_rate) + logsumexp(
[(el - 1) * eps_alpha_minus_tau, (el - 1) * _eps(el, params)], b=[1, -1]
)
terms.append(third)
signs.append(-1)
for el in range(alpha - tau + 1, alpha + 1):
fourth = _log_coef_i(alpha, el, sampling_rate) + logsumexp(
[(el - 1) * _eps(el, params), (el - 1) * eps_alpha_minus_tau],
b=[1, -1],
)
terms.append(fourth)
signs.append(1)
return (1 / (alpha - 1)) * logsumexp(terms, b=signs)
