A naive accountant fails for five distinct reasons:
I've implemented, verified, and documented the solution. Files written to /workspace/:
shuffled_ldp_accountant.py — the accountanttest_accountant.py — verification suite (all 8 tests pass)SOLUTION.md — the full write-upBelow is the solution document.
A naive accountant fails for five distinct reasons:
| # | Naive approach | Why it breaks |
|---|---|---|
| 1 | Basic composition eps = Σ eps_i |
Valid but massively loose; discards shuffle/subsampling amplification. |
| 2 | Advanced composition of (eps,δ) |
Components are not pure-DP; their privacy-loss distributions matter. |
| 3 | Homogeneous shuffle formula with max eps0 |
Not a dominating pair for heterogeneous channels; a mixed assignment of other users can be strictly worse than all-clones. |
| 4 | Use only the count S under subsampling |
The shuffled output is (N,S) = (messages, ones). Conditioning on N increases divergence, so S alone is not certified. |
| 5 | Linear-space arithmetic | At eps0=25, e^25 ≈ 7e10; probabilities underflow to 0 → 0/0. Must work in log-space. |
The right object is the privacy-loss distribution L = log(dP_A/dP_B). Independent releases compose by convolution; the certified epsilon for target delta is the smallest eps with
delta(eps) = E_A[(1 - exp(eps - L))_+] ≤ delta.
The hard part is a valid per-release PLD over all other-user assignments. The "all-clones" pair is not an upper bound for heterogeneous/subsampled channels. The certified fix is a max-survival envelope over candidate worlds:
S~(t) = max over candidates c of P_c(L_c > t).
S~ is a valid survival function, stochastically dominates every candidate, and stochastic dominance is preserved by convolution — so composed envelopes certify all assignments at once.
shuffled_ldp_accountant.py"""
shuffled_ldp_accountant.py
==========================
Tight, numerically stable privacy accountant for a pipeline of *shuffled*
local-DP releases with heterogeneous per-user binary randomized-response
(RR) channels, per-user subsampling, and adaptive repeated releases.
Every release contributes a privacy-loss distribution (PLD): the law of
L = log(dP_A/dP_B) under world A. Releases compose by convolving PLDs.
The certified epsilon is the smallest value whose hockey-stick divergence
E_A[(1 - exp(eps - L))_+] is at most the target delta.
For the worst-case over the other users' inputs we take the *max-survival
envelope* over a certified candidate family, which stochastically dominates
every candidate and composes correctly. Only numpy + stdlib are used.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from typing import List, Optional, Sequence, Tuple
import numpy as np
__all__ = ["BinaryRRChannel", "ShuffleRound", "PRVAccountant",
"CertificationResult", "poisson_binomial_logpmf"]
def poisson_binomial_logpmf(probs: Sequence[float]) -> np.ndarray:
"""log P(S = s) for a sum of independent Bernoulli variables (log-space)."""
probs = np.asarray(probs, dtype=np.float64)
logpmf = np.array([0.0])
with np.errstate(divide="ignore", invalid="ignore"):
for p in probs:
lp = math.log(p) if p > 0.0 else -np.inf
lq = math.log1p(-p) if p < 1.0 else -np.inf
nxt = np.full(logpmf.size + 1, -np.inf)
nxt[: logpmf.size] = np.logaddexp(nxt[: logpmf.size], logpmf + lq)
nxt[1:] = np.logaddexp(nxt[1:], logpmf + lp)
logpmf = nxt
return logpmf
def _joint_participation_logpmf(
user_probs: Sequence[Tuple[float, float]]
) -> np.ndarray:
"""log P(N=n, S=s) for independent (subsample, report) users.
user_probs[j] = (gamma_j, r_j): inclusion probability and P(report 1).
Tracking both N and S is essential for a certified subsampling bound.
"""
m = len(user_probs)
dp = np.full((m + 1, m + 1), -np.inf)
dp[0, 0] = 0.0
for gamma, r in user_probs:
ndp = np.full((m + 1, m + 1), -np.inf)
pz = math.log1p(-gamma) if gamma < 1.0 else -np.inf
pr = math.log(gamma * r) if gamma > 0.0 and r > 0.0 else -np.inf
pnr = (math.log(gamma * (1.0 - r))
if gamma > 0.0 and r < 1.0 else -np.inf)
for n in range(m + 1):
for s in range(n + 1):
v = dp[n, s]
if v == -np.inf:
continue
if pz > -np.inf:
ndp[n, s] = np.logaddexp(ndp[n, s], v + pz)
if n + 1 <= m:
if pnr > -np.inf:
ndp[n + 1, s] = np.logaddexp(ndp[n + 1, s], v + pnr)
if pr > -np.inf:
ndp[n + 1, s + 1] = np.logaddexp(
ndp[n + 1, s + 1], v + pr)
dp = ndp
return dp
@dataclass(frozen=True)
class BinaryRRChannel:
"""Binary randomized response that is eps-locally-DP."""
eps: float
@property
def report_truth(self) -> float:
if self.eps >= 0.0:
return 1.0 / (1.0 + math.exp(-self.eps))
e = math.exp(self.eps)
return e / (1.0 + e)
@property
def report_other(self) -> float:
return 1.0 - self.report_truth
@dataclass
class ShuffleRound:
"""A single shuffled release.
channels : one BinaryRRChannel per user
sample_rate : per-user independent inclusion probability in [0, 1]
exact_limit : if len(channels)-1 <= exact_limit enumerate all
assignments exactly; for large homogeneous groups use
O(m) candidates; for large heterogeneous groups use the
safe local-DP fallback.
"""
channels: Sequence[BinaryRRChannel]
sample_rate: float = 1.0
exact_limit: int = 11
def __post_init__(self) -> None:
self.channels = list(self.channels)
if not 0.0 <= self.sample_rate <= 1.0:
raise ValueError("sample_rate must lie in [0, 1]")
def _world_logpmfs(
self,
designated: int,
others_bits: Sequence[int],
eps_override: Optional[Sequence[float]] = None,
) -> Tuple[np.ndarray, np.ndarray]:
eps = ([c.eps for c in self.channels]
if eps_override is None else list(eps_override))
gamma = float(self.sample_rate)
m = len(eps)
def build(world_bit: int) -> np.ndarray:
rs: List[float] = []
j_other = 0
for i in range(m):
ch = BinaryRRChannel(eps[i])
if i == designated:
r = (ch.report_other if world_bit == 0
else ch.report_truth)
else:
r = (ch.report_truth if others_bits[j_other] == 1
else ch.report_other)
j_other += 1
rs.append(r)
if gamma >= 1.0:
return poisson_binomial_logpmf(rs)
return _joint_participation_logpmf([(gamma, r) for r in rs])
return build(0), build(1)
def candidate_family(self) -> List[Tuple[np.ndarray, np.ndarray]]:
"""[(losses, probs), ...] for every candidate neighbouring pair."""
m = len(self.channels)
if m == 0:
return [(np.zeros(1), np.ones(1))]
eps = [c.eps for c in self.channels]
homogeneous = (max(eps) - min(eps)) <= 1e-12
family: List[Tuple[np.ndarray, np.ndarray]] = []
if (m - 1) <= self.exact_limit:
import itertools
for designated in range(m):
for bits in itertools.product([0, 1], repeat=m - 1):
family.append(self._pair(designated, bits))
elif homogeneous:
for k in range(m): # O(m) candidates
family.append(self._pair(0, [1] * k + [0] * (m - 1 - k)))
else:
# safe no-amplification fallback (shuffle is post-processing)
d = int(np.argmax(eps))
ch = self.channels[d]
a, b = ch.report_truth, ch.report_other
if a <= 0.0 or b <= 0.0:
family.append((np.array([0.0]), np.array([1.0])))
else:
L = math.log(a / b)
family.append((np.array([L, -L]), np.array([a, b])))
return family
def _pair(self, designated, others_bits, eps_override=None):
la, lb = self._world_logpmfs(designated, others_bits, eps_override)
finite = np.isfinite(la) & np.isfinite(lb)
losses = la[finite] - lb[finite]
w = np.exp(la[finite])
tot = w.sum()
return losses, (w / tot if tot > 0 else w)
def max_survival_envelope(family):
"""Distribution whose survival function is pointwise max over family.
It stochastically dominates every member, so it is a valid dominating
PLD for the whole family.
"""
vals = np.unique(np.concatenate([L for L, _ in family]))
if vals.size == 0:
return np.array([0.0]), np.array([1.0])
prepared = []
for L, w in family:
order = np.argsort(L, kind="stable")
Ls, ws = L[order], w[order]
prepared.append((Ls, np.concatenate([[0.0], np.cumsum(ws)])))
ge = np.zeros(vals.size)
gt = np.zeros(vals.size)
for Ls, cum in prepared:
ge = np.maximum(ge, 1.0 - cum[np.searchsorted(Ls, vals, "left")])
gt = np.maximum(gt, 1.0 - cum[np.searchsorted(Ls, vals, "right")])
mass = np.maximum(ge - gt, 0.0)
s = mass.sum()
return (vals, mass / s) if s > 0 else (np.array([0.0]), np.array([1.0]))
class _GridPLD:
"""PLD discretised on a uniform grid; losses are rounded UP (conservative)."""
__slots__ = ("lo", "h", "probs")
def __init__(self, lo, h, probs):
self.lo, self.h, self.probs = float(lo), float(h), probs
@property
def values(self):
return self.lo + self.h * np.arange(self.probs.size)
@classmethod
def from_distribution(cls, values, probs, h):
if values.size == 0:
return cls(0.0, h, np.array([1.0]))
lo = math.floor(float(values.min()) / h + 1e-12) * h
# tolerance prevents float noise at a grid point from adding a bin
idx = np.maximum(np.ceil((values - lo) / h - 1e-9).astype(np.int64), 0)
arr = np.zeros(int(idx.max()) + 1)
np.add.at(arr, idx, probs)
tot = arr.sum()
return cls(lo, h, arr / tot if tot > 0 else arr)
def compose(self, other):
n = self.probs.size + other.probs.size - 1
nfft = 1 << (n - 1).bit_length()
conv = np.fft.irfft(
np.fft.rfft(self.probs, nfft) * np.fft.rfft(other.probs, nfft),
nfft)[:n]
np.clip(conv, 0.0, None, out=conv)
tot = conv.sum()
return _GridPLD(self.lo + other.lo, self.h,
conv / tot if tot > 0 else conv)
def delta(self, eps):
vals = self.values
mask = vals > eps
if not np.any(mask):
return 0.0
return float(np.sum(self.probs[mask] * (-np.expm1(eps - vals[mask]))))
def epsilon(self, delta):
vals = self.values
hi = float(vals[-1])
lo = min(float(vals[0]), 0.0) - 1.0
while self.delta(lo) <= delta:
lo -= 1.0
for _ in range(200):
mid = 0.5 * (lo + hi)
if self.delta(mid) <= delta:
hi = mid
else:
lo = mid
return hi
@dataclass
class CertificationResult:
epsilon: float
delta: float
per_step_epsilon: List[float]
leave_one_out_gain: List[float]
dominant_step: Optional[int]
grid_size: float
def summary(self):
lines = [f"epsilon = {self.epsilon:.6f} (delta = {self.delta:g},"
f" grid = {self.grid_size:g})"]
for i, (e, g) in enumerate(zip(self.per_step_epsilon,
self.leave_one_out_gain)):
star = " <-- dominant" if i == self.dominant_step else ""
lines.append(f" step {i:>3}: standalone eps = {e:.6f}, "
f"leave-one-out gain = {g:.6f}{star}")
return "\n".join(lines)
class PRVAccountant:
"""PLD accountant over an ordered transcript of shuffled releases."""
def __init__(self, grid_size: float = 1e-3):
self.h = float(grid_size)
def certify(self, transcript, delta, return_dominant=True):
if not 0.0 < delta < 1.0:
raise ValueError("delta must lie in (0, 1)")
rounds = list(transcript)
base = self._compose(rounds)
eps = base.epsilon(delta)
per_step = [self._compose([r]).epsilon(delta) for r in rounds]
dominant, gains = None, []
if return_dominant and rounds:
for i in range(len(rounds)):
without = self._compose(rounds[:i] + rounds[i + 1:])
gains.append(eps - without.epsilon(delta))
dominant = int(np.argmax(gains))
return CertificationResult(eps, delta, per_step, gains, dominant,
self.h)
def round_grid(self, rnd):
vals, probs = max_survival_envelope(rnd.candidate_family())
return _GridPLD.from_distribution(vals, probs, self.h)
def _compose(self, rounds):
if not rounds:
return _GridPLD(0.0, self.h, np.array([1.0]))
acc = self.round_grid(rounds[0])
for rnd in rounds[1:]:
acc = acc.compose(self.round_grid(rnd))
return acc
def homogeneous_round(n, eps, sample_rate=1.0, **kw):
return ShuffleRound([BinaryRRChannel(eps) for _ in range(n)],
sample_rate, **kw)
def heterogeneous_round(epsilons, sample_rate=1.0, **kw):
return ShuffleRound([BinaryRRChannel(e) for e in epsilons],
sample_rate, **kw)
(N,S) DP — no underflow at extreme ratios.(N,S) under subsampling — using S alone is not certified.+grid_size bug seen at exact grid points.eps - eps_without_i, argmax = dominating step.Run:
python3 test_accountant.py
Recorded result (all 8 tests pass):
[1] max-survival envelope dominates every candidate PLD: OK
[2] no-amplification fallback is a valid local-DP certificate: OK
[3] delta=0.1: grid=2.041573 exact=2.041510 gap=6.26e-05
delta=0.01: grid=3.226450 exact=3.226450 gap=2.66e-15
delta=0.001: grid=3.661549 exact=3.661549 gap=1.78e-15
[4] n=3: eps=4.0000 | n=5: eps=4.0000 | n=20: eps=4.0000 | n=100: eps=3.9999
[5] gamma=1.0 -> eps=1.9795; gamma=0.2 -> eps=0.4966
[6] epsilon = 7.999442 (delta=1e-05, grid=0.001)
step 0: standalone 1.999969, gain 1.999620
step 1: standalone 4.999989, gain 4.999967 <-- dominant
step 2: standalone 0.999832, gain 0.999475
[7] local eps=25, 2 releases -> eps=50.0000 (finite, no NaN)
[8] epsilon monotone non-increasing in delta: OK
ALL TESTS PASSED
End-to-end demo:
from shuffled_ldp_accountant import PRVAccountant, homogeneous_round, heterogeneous_round
acc = PRVAccountant(grid_size=1e-3)
transcript = [
homogeneous_round(40, 1.2, sample_rate=0.3),
heterogeneous_round([0.5]*7 + [2.8] + [0.6]*32, sample_rate=1.0),
homogeneous_round(40, 0.8, sample_rate=0.5),
]
print(acc.certify(transcript, delta=1e-5).summary())
epsilon = 3.309248 (delta = 1e-05, grid = 0.001)
step 0: standalone eps = 0.372516, leave-one-out gain = 0.187800
step 1: standalone eps = 2.799989, leave-one-out gain = 2.797851 <-- dominant
step 2: standalone eps = 0.322807, leave-one-out gain = 0.138181
Interpretation. Binary RR with large eps0 amplifies only weakly by shuffling (the local blanket is far from uniform), and the accountant correctly reports this rather than an optimistic value. Subsampling yields strong amplification (1.98 → 0.50). Composition is additive in privacy loss, and the eps=5 release is correctly identified as dominant. Grid losses are rounded up, so the returned epsilon is always a certified upper bound; reduce grid_size for extra tightness.
# Evidence - Problem class: python-shuffled-ldp-privacy-amplification-accountant - Model: openrouter/deepseek/deepseek-v4.1-flash - Solved: 2026-10-10T10:27:58.692Z - Verification: solution produced by pi in sandbox; see signatures.json
{"description": "Implement a tight privacy accountant for a shuffled local-DP pipeline with heterogeneous per-user randomized-response channels, subsampling, and adaptive repeated releases. Given a target delta and an ordered transcript of mechanisms, compute a certified epsilon bound and identify the composition step that dominates while preserving numerical stability at extreme likelihood ratios.", "environment": "python3", "language": "python", "model": "openrouter/deepseek/deepseek-v4.1-flash", "problem_class": "python-shuffled-ldp-privacy-amplification-accountant", "provider": "openrouter", "solved_at": "2026-10-10T10:27:58.692Z", "version": "3.11"}