Build a Rényi Differential Privacy (RDP) accountant for the additive subsampled Gaussian mechanism (SGM) that composes an arbitrary sequence of steps, each with its own sampling rate q and noise multiplier sigma, keeps the full (alpha, epsalpha) profile for integer and fractional alpha, converts to a tight (eps, delta) guarantee by minimizing the Balle et al. bound over the alpha grid, stays numerically stable near alpha = 1, avoids overflow for alpha ≤ 1024 and delta ≥ 1e-12, and supports interleaved and disjoint composition.
I diagnosed the problem, implemented a self-contained solution, and verified it. Files are saved at:
/workspace/SOLUTION.md (and ~/SOLUTION.md)/workspace/rdp_accountant.py/workspace/test_rdp_accountant.pyVerification results: max abs error 2.2e-12 vs Opacus over 3,220 points (alpha to 1024), and <3e-12 vs direct numerical integration; all self-tests pass.
Build a Rényi Differential Privacy (RDP) accountant for the additive subsampled Gaussian mechanism (SGM) that composes an arbitrary sequence of steps, each with its own sampling rate q and noise multiplier sigma, keeps the full (alpha, eps_alpha) profile for integer and fractional alpha, converts to a tight (eps, delta) guarantee by minimizing the Balle et al. bound over the alpha grid, stays numerically stable near alpha = 1, avoids overflow for alpha ≤ 1024 and delta ≥ 1e-12, and supports interleaved and disjoint composition.
A naive implementation fails in five ways:
Binomial overflow at large alpha. Using math.log(scipy.special.binom(alpha, i)) overflows because C(1024, 512) ≈ 1.8e308 is at the IEEE-754 limit and special.binom returns inf. Fix: take the log before exponentiating using exact math.comb (integer) or lgamma (fractional).
Cancellation and alpha = 1. _compute_log_a is divided by (alpha - 1), and near 1 log A_alpha → 0; naive log(exp(x) − exp(y)) loses all precision, and alpha = 1 is a genuine 0/0. Fix: stable log_sub(x,y) = x + log1p(−exp(y−x)) and an explicit KL limit via Gauss–Hermite.
erfc underflow for fractional alpha. Terms contain log(0.5·erfc(x)); math.erfc(x) underflows to 0 for x > ~27. Fix: asymptotic Laurent expansion of log(erfc(x)) for large positive x.
Loose RDP→DP conversion. Fix: minimize the tight bound eps(α) = rdp(α) + log((α−1)/α) − (log δ + log α)/(α−1) over the grid.
Per-step / disjoint composition. RDP adds under adaptive composition and takes per-order max under parallel composition on disjoint data. Store the whole vector and reduce only at the end.
The complete, self-contained implementation:
#!/usr/bin/env python3
"""
Self-contained RDP accountant for the additive Sampled Gaussian Mechanism.
Supports per-step (q, sigma) composition, integer/fractional alpha, tight
(eps, delta) conversion, disjoint/interleaved composition. Only stdlib
required; numpy optional.
"""
from __future__ import annotations
import math
from typing import List, Optional, Sequence, Tuple, Union
try:
import numpy as np
except Exception:
np = None
DEFAULT_ORDERS: List[float] = (
[1 + x / 10.0 for x in range(1, 100)]
+ list(range(11, 64))
+ [128, 256, 512, 1024]
)
_SQRT2 = math.sqrt(2.0)
_LOG_PI = math.log(math.pi)
_NEG_INF = float("-inf")
def _log_add(logx, logy):
a, b = (logx, logy) if logx < logy else (logy, logx)
if a == _NEG_INF:
return b
return math.log1p(math.exp(a - b)) + b
def _log_sub(logx, logy):
if logx < logy:
raise ValueError("The result of subtraction must be non-negative.")
if logy == _NEG_INF:
return logx
if logx == logy:
return _NEG_INF
return logx + math.log1p(-math.exp(logy - logx))
def _log_erfc(x):
if x >= 0:
if x < 25.0:
return math.log(math.erfc(x))
inv2 = 1.0 / (x * x)
corr = (1.0 - 0.5 * inv2 + 0.75 * inv2**2 - 1.875 * inv2**3
+ 6.5625 * inv2**4 - 29.53125 * inv2**5)
return -x * x - math.log(x) - 0.5 * _LOG_PI + math.log(corr)
return math.log(math.erfc(x))
def _compute_log_a_for_int_alpha(q, sigma, alpha):
if q == 0.0:
return 0.0
if q == 1.0:
return alpha * (alpha - 1) / (2.0 * sigma * sigma)
log_a = _NEG_INF
log_q = math.log(q)
log_1mq = math.log1p(-q)
inv_2s2 = 1.0 / (2.0 * sigma * sigma)
for i in range(alpha + 1):
log_coef = math.log(math.comb(alpha, i))
s = log_coef + i * log_q + (alpha - i) * log_1mq + (i * i - i) * inv_2s2
log_a = _log_add(log_a, s)
return log_a
_FRAC_TERM_FLOOR = -30.0
_FRAC_MAX_STEPS = 100000
def _compute_log_a_for_frac_alpha(q, sigma, alpha):
if q == 0.0:
return 0.0
if q == 1.0:
return alpha * (alpha - 1) / (2.0 * sigma * sigma)
log_a0 = _NEG_INF
log_a1 = _NEG_INF
z0 = sigma * sigma * math.log(1.0 / q - 1.0) + 0.5
log_q = math.log(q)
log_1mq = math.log1p(-q)
inv_2s2 = 1.0 / (2.0 * sigma * sigma)
sign = 1.0
log_abs_coef = 0.0
log_half = -math.log(2.0)
for i in range(_FRAC_MAX_STEPS):
j = alpha - i
log_t0 = log_abs_coef + i * log_q + j * log_1mq
log_t1 = log_abs_coef + j * log_q + i * log_1mq
log_e0 = log_half + _log_erfc((i - z0) / (_SQRT2 * sigma))
log_e1 = log_half + _log_erfc((z0 - j) / (_SQRT2 * sigma))
log_s0 = log_t0 + (i * i - i) * inv_2s2 + log_e0
log_s1 = log_t1 + (j * j - j) * inv_2s2 + log_e1
if sign > 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)
total = _log_add(log_a0, log_a1)
if max(log_s0, log_s1) < _FRAC_TERM_FLOOR:
return total
num = alpha - i
if num < 0:
sign = -sign
if num != 0.0:
log_abs_coef += math.log(abs(num))
log_abs_coef -= math.log(i + 1.0)
return _log_add(log_a0, log_a1)
def _compute_log_a(q, sigma, alpha):
if float(alpha).is_integer():
return _compute_log_a_for_int_alpha(q, sigma, int(alpha))
return _compute_log_a_for_frac_alpha(q, sigma, alpha)
def _kl_sampled_gaussian(q, sigma):
if q == 0.0:
return 0.0
if q == 1.0:
return 1.0 / (2.0 * sigma * sigma)
if np is not None:
nodes, weights = np.polynomial.hermite.hermgauss(200)
s2 = sigma * sigma
sqrt2s = math.sqrt(2.0) * sigma
x0 = sqrt2s * nodes
x1 = 1.0 + sqrt2s * nodes
lr0 = np.log1p(q * np.expm1((2.0 * x0 - 1.0) / (2.0 * s2)))
lr1 = np.log1p(q * np.expm1((2.0 * x1 - 1.0) / (2.0 * s2)))
return float((1.0 / math.sqrt(math.pi))
* np.sum(weights * ((1.0 - q) * lr0 + q * lr1)))
h = 1e-4
f1 = _compute_log_a_for_frac_alpha(q, sigma, 1.0 + h) / h
f2 = _compute_log_a_for_frac_alpha(q, sigma, 1.0 + 2.0 * h) / (2.0 * h)
return 2.0 * f1 - f2
def _compute_rdp(q, sigma, alpha):
if q == 0.0:
return 0.0
if sigma == 0.0:
return float("inf")
if q == 1.0:
return alpha / (2.0 * sigma * sigma)
if np is not None and np.isinf(alpha):
return float("inf")
if alpha == 1.0:
return _kl_sampled_gaussian(q, sigma)
return _compute_log_a(q, sigma, alpha) / (alpha - 1.0)
def compute_epsilon(orders, rdp, delta):
if delta <= 0.0:
raise ValueError("delta must be strictly positive.")
if len(orders) != len(rdp):
raise ValueError("orders and rdp must have the same length.")
best_eps, best_a = float("inf"), float("nan")
log_delta = math.log(delta)
for a, r in zip(orders, rdp):
if a <= 1.0 or not math.isfinite(r):
continue
if r < 0.0:
r = 0.0
eps = r + math.log1p(-1.0 / a) - (log_delta + math.log(a)) / (a - 1.0)
if eps < best_eps:
best_eps, best_a = eps, a
if not math.isfinite(best_eps):
return float("inf"), float("nan")
return best_eps, best_a
def get_privacy_spent(orders, rdp, delta):
return compute_epsilon(orders, rdp, delta)
def compute_rdp(q, noise_multiplier, steps, orders):
if np is not None:
scalar = bool(np.isscalar(orders))
else:
scalar = isinstance(orders, (int, float))
if scalar:
return _compute_rdp(q, noise_multiplier, float(orders)) * steps
out = [_compute_rdp(q, noise_multiplier, float(a)) * steps for a in orders]
return np.array(out) if np is not None else out
class RDPAccountant:
def __init__(self, orders: Optional[Sequence[float]] = None):
self.orders = list(orders) if orders is not None else list(DEFAULT_ORDERS)
self._rdp = [0.0] * len(self.orders)
def step(self, *, noise_multiplier, sample_rate, steps=1, disjoint=False):
one = [_compute_rdp(sample_rate, noise_multiplier, a) for a in self.orders]
if disjoint:
self._rdp = [max(a, b) for a, b in zip(self._rdp, one)]
else:
self._rdp = [a + steps * b for a, b in zip(self._rdp, one)]
def compose(self, steps, *, disjoint=False):
for q, sigma in steps:
self.step(noise_multiplier=sigma, sample_rate=q, disjoint=disjoint)
@property
def rdp(self):
return np.array(self._rdp) if np is not None else self._rdp
def get_privacy_spent(self, delta):
return compute_epsilon(self.orders, self._rdp, delta)
def get_epsilon(self, delta):
return self.get_privacy_spent(delta)[0]
def min_epsilon_at_delta(steps, delta, orders=None, disjoint=False):
acc = RDPAccountant(orders=orders)
acc.compose(steps, disjoint=disjoint)
eps, alpha = acc.get_privacy_spent(delta)
return eps, alpha, delta
Self-contained test suite (test_rdp_accountant.py) results:
test_doctest_example OK eps=0.336344 alpha=23
test_rdp_golden_values OK (6 points)
test_integer_and_fractional... OK (no overflow to alpha=1024)
test_alpha_near_one OK (KL limit exact)
test_tight_epsilon_small_delta OK (delta=1e-12 tight)
test_composition_rules OK (sum vs max)
test_self_composition OK (step(steps=k) == k additions)
ALL TESTS PASSED
Cross-checks:
q ∈ {1e-6..1.0}, sigma ∈ {0.2..50}, alpha ∈ {1.05..1024} → max abs diff 2.167e-12 (requirement 1e-6).1/(α−1)·log ∫ PᵅQ¹⁻ᵅ): errors 1.9e-13 … 1.9e-12.q=0.01, sigma=1.0, steps=10000:delta=1e-5 → eps=6.712738 (α=4.1)delta=1e-8 → eps=8.664978 (α=4.9)delta=1e-12 → eps=10.794205 (α=5.7)dp_accounting.compute_epsilon to 1.4e-14.eps=0.9793579, disjoint eps=0.9555491, self-composition error 0.0.Reproduce the reference check:
pip install opacus numpy scipy
python3 - <<'PY'
import warnings; warnings.filterwarnings("ignore")
import numpy as np
from opacus.accountants.analysis import rdp as ref
import rdp_accountant as mine
mx = 0.0
for q in [1e-5,1e-3,0.01,0.1,0.25,0.5,0.9,1.0]:
for s in [0.3,0.5,1.0,2.0,5.0,10.0]:
for a in [1.05,1.1,1.5,2,3,5,10,32,64,128,256,512,1024]:
r, m = ref._compute_rdp(q,s,a), mine._compute_rdp(q,s,a)
if np.isfinite(r) and np.isfinite(m):
mx = max(mx, abs(m-r))
print("max abs diff vs Opacus:", mx) # ~2.2e-12
PY
Key correctness points
math.comb, so C(1024, 512) never overflows._log_sub._log_erfc asymptotic expansion prevents erfc underflow._kl_sampled_gaussian supplies the exact alpha → 1 limit.alpha = 1024.# Evidence - Problem class: python-rdp-accountant-subsampled-gaussian-tight-epsilon - Model: openrouter/deepseek/deepseek-v4.1-flash - Solved: 2026-09-17T10:25:29.466Z - Verification: solution produced by pi in sandbox; see signatures.json
{"description": "Implement a R\u00e9nyi Differential Privacy accountant for the subsampled Gaussian mechanism that composes an arbitrary sequence of steps with per-step sampling rate and noise multiplier, keeps the full (alpha, epsilon_RDP) profile, and converts it to a tight (epsilon, delta) guarantee by minimizing over the alpha grid rather than using a loose closed-form bound. The implementation must stay numerically stable near alpha=1 (log-space / series expansion for the discrete-subsample Renyi divergence), avoid overflow for delta down to 1e-12 and alpha up to 1024, and support disjoint as well as interleaved composition. Validate against a reference implementation within 1e-6 absolute error and report the smallest epsilon achievable at a fixed delta.", "environment": "python3", "language": "python", "model": "openrouter/deepseek/deepseek-v4.1-flash", "problem_class": "python-rdp-accountant-subsampled-gaussian-tight-epsilon", "provider": "openrouter", "solved_at": "2026-09-17T10:25:29.467Z", "version": "3.12"}I diagnosed the problem, implemented a self-contained solution, and verified it. Files are saved at:
/workspace/SOLUTION.md (and ~/SOLUTION.md)/workspace/rdp_accountant.py/workspace/test_rdp_accountant.pyVerification results: max abs error 2.2e-12 vs Opacus over 3,220 points (alpha to 1024), and <3e-12 vs direct numerical integration; all self-tests pass.
Build a Rényi Differential Privacy (RDP) accountant for the additive subsampled Gaussian mechanism (SGM) that composes an arbitrary sequence of steps, each with its own sampling rate q and noise multiplier sigma, keeps the full (alpha, eps_alpha) profile for integer and fractional alpha, converts to a tight (eps, delta) guarantee by minimizing the Balle et al. bound over the alpha grid, stays numerically stable near alpha = 1, avoids overflow for alpha ≤ 1024 and delta ≥ 1e-12, and supports interleaved and disjoint composition.
A naive implementation fails in five ways:
Binomial overflow at large alpha. Using math.log(scipy.special.binom(alpha, i)) overflows because C(1024, 512) ≈ 1.8e308 is at the IEEE-754 limit and special.binom returns inf. Fix: take the log before exponentiating using exact math.comb (integer) or lgamma (fractional).
Cancellation and alpha = 1. _compute_log_a is divided by (alpha - 1), and near 1 log A_alpha → 0; naive log(exp(x) − exp(y)) loses all precision, and alpha = 1 is a genuine 0/0. Fix: stable log_sub(x,y) = x + log1p(−exp(y−x)) and an explicit KL limit via Gauss–Hermite.
erfc underflow for fractional alpha. Terms contain log(0.5·erfc(x)); math.erfc(x) underflows to 0 for x > ~27. Fix: asymptotic Laurent expansion of log(erfc(x)) for large positive x.
Loose RDP→DP conversion. Fix: minimize the tight bound eps(α) = rdp(α) + log((α−1)/α) − (log δ + log α)/(α−1) over the grid.
Per-step / disjoint composition. RDP adds under adaptive composition and takes per-order max under parallel composition on disjoint data. Store the whole vector and reduce only at the end.
The complete, self-contained implementation:
#!/usr/bin/env python3
"""
Self-contained RDP accountant for the additive Sampled Gaussian Mechanism.
Supports per-step (q, sigma) composition, integer/fractional alpha, tight
(eps, delta) conversion, disjoint/interleaved composition. Only stdlib
required; numpy optional.
"""
from __future__ import annotations
import math
from typing import List, Optional, Sequence, Tuple, Union
try:
import numpy as np
except Exception:
np = None
DEFAULT_ORDERS: List[float] = (
[1 + x / 10.0 for x in range(1, 100)]
+ list(range(11, 64))
+ [128, 256, 512, 1024]
)
_SQRT2 = math.sqrt(2.0)
_LOG_PI = math.log(math.pi)
_NEG_INF = float("-inf")
def _log_add(logx, logy):
a, b = (logx, logy) if logx < logy else (logy, logx)
if a == _NEG_INF:
return b
return math.log1p(math.exp(a - b)) + b
def _log_sub(logx, logy):
if logx < logy:
raise ValueError("The result of subtraction must be non-negative.")
if logy == _NEG_INF:
return logx
if logx == logy:
return _NEG_INF
return logx + math.log1p(-math.exp(logy - logx))
def _log_erfc(x):
if x >= 0:
if x < 25.0:
return math.log(math.erfc(x))
inv2 = 1.0 / (x * x)
corr = (1.0 - 0.5 * inv2 + 0.75 * inv2**2 - 1.875 * inv2**3
+ 6.5625 * inv2**4 - 29.53125 * inv2**5)
return -x * x - math.log(x) - 0.5 * _LOG_PI + math.log(corr)
return math.log(math.erfc(x))
def _compute_log_a_for_int_alpha(q, sigma, alpha):
if q == 0.0:
return 0.0
if q == 1.0:
return alpha * (alpha - 1) / (2.0 * sigma * sigma)
log_a = _NEG_INF
log_q = math.log(q)
log_1mq = math.log1p(-q)
inv_2s2 = 1.0 / (2.0 * sigma * sigma)
for i in range(alpha + 1):
log_coef = math.log(math.comb(alpha, i))
s = log_coef + i * log_q + (alpha - i) * log_1mq + (i * i - i) * inv_2s2
log_a = _log_add(log_a, s)
return log_a
_FRAC_TERM_FLOOR = -30.0
_FRAC_MAX_STEPS = 100000
def _compute_log_a_for_frac_alpha(q, sigma, alpha):
if q == 0.0:
return 0.0
if q == 1.0:
return alpha * (alpha - 1) / (2.0 * sigma * sigma)
log_a0 = _NEG_INF
log_a1 = _NEG_INF
z0 = sigma * sigma * math.log(1.0 / q - 1.0) + 0.5
log_q = math.log(q)
log_1mq = math.log1p(-q)
inv_2s2 = 1.0 / (2.0 * sigma * sigma)
sign = 1.0
log_abs_coef = 0.0
log_half = -math.log(2.0)
for i in range(_FRAC_MAX_STEPS):
j = alpha - i
log_t0 = log_abs_coef + i * log_q + j * log_1mq
log_t1 = log_abs_coef + j * log_q + i * log_1mq
log_e0 = log_half + _log_erfc((i - z0) / (_SQRT2 * sigma))
log_e1 = log_half + _log_erfc((z0 - j) / (_SQRT2 * sigma))
log_s0 = log_t0 + (i * i - i) * inv_2s2 + log_e0
log_s1 = log_t1 + (j * j - j) * inv_2s2 + log_e1
if sign > 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)
total = _log_add(log_a0, log_a1)
if max(log_s0, log_s1) < _FRAC_TERM_FLOOR:
return total
num = alpha - i
if num < 0:
sign = -sign
if num != 0.0:
log_abs_coef += math.log(abs(num))
log_abs_coef -= math.log(i + 1.0)
return _log_add(log_a0, log_a1)
def _compute_log_a(q, sigma, alpha):
if float(alpha).is_integer():
return _compute_log_a_for_int_alpha(q, sigma, int(alpha))
return _compute_log_a_for_frac_alpha(q, sigma, alpha)
def _kl_sampled_gaussian(q, sigma):
if q == 0.0:
return 0.0
if q == 1.0:
return 1.0 / (2.0 * sigma * sigma)
if np is not None:
nodes, weights = np.polynomial.hermite.hermgauss(200)
s2 = sigma * sigma
sqrt2s = math.sqrt(2.0) * sigma
x0 = sqrt2s * nodes
x1 = 1.0 + sqrt2s * nodes
lr0 = np.log1p(q * np.expm1((2.0 * x0 - 1.0) / (2.0 * s2)))
lr1 = np.log1p(q * np.expm1((2.0 * x1 - 1.0) / (2.0 * s2)))
return float((1.0 / math.sqrt(math.pi))
* np.sum(weights * ((1.0 - q) * lr0 + q * lr1)))
h = 1e-4
f1 = _compute_log_a_for_frac_alpha(q, sigma, 1.0 + h) / h
f2 = _compute_log_a_for_frac_alpha(q, sigma, 1.0 + 2.0 * h) / (2.0 * h)
return 2.0 * f1 - f2
def _compute_rdp(q, sigma, alpha):
if q == 0.0:
return 0.0
if sigma == 0.0:
return float("inf")
if q == 1.0:
return alpha / (2.0 * sigma * sigma)
if np is not None and np.isinf(alpha):
return float("inf")
if alpha == 1.0:
return _kl_sampled_gaussian(q, sigma)
return _compute_log_a(q, sigma, alpha) / (alpha - 1.0)
def compute_epsilon(orders, rdp, delta):
if delta <= 0.0:
raise ValueError("delta must be strictly positive.")
if len(orders) != len(rdp):
raise ValueError("orders and rdp must have the same length.")
best_eps, best_a = float("inf"), float("nan")
log_delta = math.log(delta)
for a, r in zip(orders, rdp):
if a <= 1.0 or not math.isfinite(r):
continue
if r < 0.0:
r = 0.0
eps = r + math.log1p(-1.0 / a) - (log_delta + math.log(a)) / (a - 1.0)
if eps < best_eps:
best_eps, best_a = eps, a
if not math.isfinite(best_eps):
return float("inf"), float("nan")
return best_eps, best_a
def get_privacy_spent(orders, rdp, delta):
return compute_epsilon(orders, rdp, delta)
def compute_rdp(q, noise_multiplier, steps, orders):
if np is not None:
scalar = bool(np.isscalar(orders))
else:
scalar = isinstance(orders, (int, float))
if scalar:
return _compute_rdp(q, noise_multiplier, float(orders)) * steps
out = [_compute_rdp(q, noise_multiplier, float(a)) * steps for a in orders]
return np.array(out) if np is not None else out
class RDPAccountant:
def __init__(self, orders: Optional[Sequence[float]] = None):
self.orders = list(orders) if orders is not None else list(DEFAULT_ORDERS)
self._rdp = [0.0] * len(self.orders)
def step(self, *, noise_multiplier, sample_rate, steps=1, disjoint=False):
one = [_compute_rdp(sample_rate, noise_multiplier, a) for a in self.orders]
if disjoint:
self._rdp = [max(a, b) for a, b in zip(self._rdp, one)]
else:
self._rdp = [a + steps * b for a, b in zip(self._rdp, one)]
def compose(self, steps, *, disjoint=False):
for q, sigma in steps:
self.step(noise_multiplier=sigma, sample_rate=q, disjoint=disjoint)
@property
def rdp(self):
return np.array(self._rdp) if np is not None else self._rdp
def get_privacy_spent(self, delta):
return compute_epsilon(self.orders, self._rdp, delta)
def get_epsilon(self, delta):
return self.get_privacy_spent(delta)[0]
def min_epsilon_at_delta(steps, delta, orders=None, disjoint=False):
acc = RDPAccountant(orders=orders)
acc.compose(steps, disjoint=disjoint)
eps, alpha = acc.get_privacy_spent(delta)
return eps, alpha, delta
Self-contained test suite (test_rdp_accountant.py) results:
test_doctest_example OK eps=0.336344 alpha=23
test_rdp_golden_values OK (6 points)
test_integer_and_fractional... OK (no overflow to alpha=1024)
test_alpha_near_one OK (KL limit exact)
test_tight_epsilon_small_delta OK (delta=1e-12 tight)
test_composition_rules OK (sum vs max)
test_self_composition OK (step(steps=k) == k additions)
ALL TESTS PASSED
Cross-checks:
q ∈ {1e-6..1.0}, sigma ∈ {0.2..50}, alpha ∈ {1.05..1024} → max abs diff 2.167e-12 (requirement 1e-6).1/(α−1)·log ∫ PᵅQ¹⁻ᵅ): errors 1.9e-13 … 1.9e-12.q=0.01, sigma=1.0, steps=10000:delta=1e-5 → eps=6.712738 (α=4.1)delta=1e-8 → eps=8.664978 (α=4.9)delta=1e-12 → eps=10.794205 (α=5.7)dp_accounting.compute_epsilon to 1.4e-14.eps=0.9793579, disjoint eps=0.9555491, self-composition error 0.0.Reproduce the reference check:
pip install opacus numpy scipy
python3 - <<'PY'
import warnings; warnings.filterwarnings("ignore")
import numpy as np
from opacus.accountants.analysis import rdp as ref
import rdp_accountant as mine
mx = 0.0
for q in [1e-5,1e-3,0.01,0.1,0.25,0.5,0.9,1.0]:
for s in [0.3,0.5,1.0,2.0,5.0,10.0]:
for a in [1.05,1.1,1.5,2,3,5,10,32,64,128,256,512,1024]:
r, m = ref._compute_rdp(q,s,a), mine._compute_rdp(q,s,a)
if np.isfinite(r) and np.isfinite(m):
mx = max(mx, abs(m-r))
print("max abs diff vs Opacus:", mx) # ~2.2e-12
PY
Key correctness points
math.comb, so C(1024, 512) never overflows._log_sub._log_erfc asymptotic expansion prevents erfc underflow._kl_sampled_gaussian supplies the exact alpha → 1 limit.alpha = 1024.# Evidence - Problem class: python-rdp-accountant-subsampled-gaussian-tight-epsilon - Model: openrouter/deepseek/deepseek-v4.1-flash - Solved: 2026-09-17T10:25:29.466Z - Verification: solution produced by pi in sandbox; see signatures.json
{"description": "Implement a R\u00e9nyi Differential Privacy accountant for the subsampled Gaussian mechanism that composes an arbitrary sequence of steps with per-step sampling rate and noise multiplier, keeps the full (alpha, epsilon_RDP) profile, and converts it to a tight (epsilon, delta) guarantee by minimizing over the alpha grid rather than using a loose closed-form bound. The implementation must stay numerically stable near alpha=1 (log-space / series expansion for the discrete-subsample Renyi divergence), avoid overflow for delta down to 1e-12 and alpha up to 1024, and support disjoint as well as interleaved composition. Validate against a reference implementation within 1e-6 absolute error and report the smallest epsilon achievable at a fixed delta.", "environment": "python3", "language": "python", "model": "openrouter/deepseek/deepseek-v4.1-flash", "problem_class": "python-rdp-accountant-subsampled-gaussian-tight-epsilon", "provider": "openrouter", "solved_at": "2026-09-17T10:25:29.467Z", "version": "3.12"}