Privacy-Preserving Experiment Design Questions

import torch
from opacus import PrivacyEngine

model = MyModel().to(device)
optimizer = torch.optim.SGD(model.parameters(), lr=0.1)
privacy_engine = PrivacyEngine(
    model,
    sample_rate=batch_size / len(train_dataset),  # if using Poisson sampling: batch_size/len
    alphas=[10, 100],
    noise_multiplier=noise_multiplier,
    max_grad_norm=clip_norm,
)
privacy_engine.attach(optimizer)

eps, best_delta = privacy_engine.accountant.get_privacy_spent(delta=1e-5)
print(f"After epoch {epoch}: ε = {eps:.3f}, δ = {best_delta}")

import torch
from opacus import PrivacyEngine
# model, dataloader, optimizer_d defined
privacy_engine = PrivacyEngine(
    module=disc,
    sample_rate=batch_size / len(dataset),
    alphas=[1 + x/10. for x in range(1, 100)],  # RDP orders
    noise_multiplier=noise_multiplier,
    max_grad_norm=clip_norm,
)
privacy_engine.attach(optimizer_d)

for epoch in range(epochs):
    for real_batch in loader:
        # standard GAN discriminator step
        optimizer_d.zero_grad()
        loss_d = discriminator_loss(real_batch, gen)
        loss_d.backward()
        optimizer_d.step()
    eps, best_alpha = privacy_engine.get_privacy_spent(delta=target_delta)
    print(f"Epoch {epoch} ε = {eps:.3f}, δ = {target_delta}")

from functorch import make_functional_with_buffers, vmap, grad

# make functional for generator
f_gen, params_gen, buffers = make_functional_with_buffers(gen)
def loss_fn(params, buffers, z, real_labels):
    fake = f_gen(params, buffers, z)
    return generator_loss_from_disc(fake, disc)

# compute per-sample grads
per_sample_grads = vmap(grad(loss_fn), in_dims=(None, None, 0, 0))(params_gen, buffers, z_batch, real_labels_batch)
# clip per-sample grads (L2 across parameters)
def clip_and_aggregate(per_sample_grads, max_norm):
    norms = torch.stack([torch.sqrt(sum((g**2).sum() for g in sample))) for sample in zip(*per_sample_grads)])
    scales = (max_norm / (norms + 1e-6)).clamp(max=1.0)
    # scale each sample's grads then sum and add noise
    aggregated = [ (g * scales.view(-1, *([1]* (g.dim()-1)))).sum(0) for g in per_sample_grads ]
    return aggregated

agg_grads = clip_and_aggregate(per_sample_grads, clip_norm)
# add Gaussian noise
for g in agg_grads:
    g += torch.normal(0, noise_multiplier*clip_norm, size=g.shape, device=g.device) / batch_size
# load into optimizer & step (manual param update or set gradients then optimizer.step)

import tensorflow as tf
from tensorflow_privacy.privacy.optimizers.dp_optimizer_keras import DPKerasAdamOptimizer
from tensorflow_privacy.privacy.analysis.rdp_accountant import compute_rdp, get_privacy_spent

optimizer_d = DPKerasAdamOptimizer(
    l2_norm_clip=clip_norm,
    noise_multiplier=noise_multiplier,
    num_microbatches=batch_size,  # or set microbatches <= batch_size
    learning_rate=lr
)
# compile discriminator with optimizer_d and train via gradient tape
# compute RDP epsilon each epoch:
orders = [1 + x/10. for x in range(1, 100)]
sampling_prob = batch_size / dataset_size
rdp = compute_rdp(q=sampling_prob, noise_multiplier=noise_multiplier, steps=steps, orders=orders)
eps, opt_order = get_privacy_spent(orders, rdp, target_delta)
print(eps)

import math
import numpy as np

def log_comb(n, k):
    # log binomial using log gamma
    return math.lgamma(n+1) - math.lgamma(k+1) - math.lgamma(n-k+1)

def logsumexp(arr):
    a = max(arr)
    return a + math.log(sum(math.exp(x - a) for x in arr))

def rdp_poisson_subsampled_gaussian(alpha, q, sigma):
    """
    Compute RDP at integer alpha >= 2 for Poisson subsampling + Gaussian noise.
    Uses log-space sums for stability.
    Formula: RDP = 1/(alpha-1) * log( sum_{j=0}^alpha C(alpha,j) q^j (1-q)^{alpha-j} exp((j^2 - j)/(2 sigma^2)) )
    """
    if alpha < 2 or int(alpha) != alpha:
        raise ValueError("This implementation supports integer alpha >= 2")
    alpha = int(alpha)
    log_terms = []
    log_q = math.log(q) if q > 0 else float("-inf")
    log_1_q = math.log(1 - q) if q < 1 else float("-inf")
    coeff = 1.0  # placeholder
    for j in range(0, alpha+1):
        # log C(alpha, j) + j*log q + (alpha-j)*log(1-q) + ((j^2 - j)/(2 sigma^2))
        log_p = log_comb(alpha, j)
        log_p += j * log_q if log_q != float("-inf") else (0.0 if j == 0 else float("-inf"))
        log_p += (alpha - j) * log_1_q if log_1_q != float("-inf") else (0.0 if alpha - j == 0 else float("-inf"))
        log_p += (j*j - j) / (2.0 * sigma * sigma)
        log_terms.append(log_p)
    log_A = logsumexp(log_terms)
    return (log_A) / (alpha - 1.0)

def get_epsilon(q, sigma, T, delta, max_alpha=200):
    """
    Compose RDP over T steps (additive) and convert to (epsilon, delta).
    Search over integer alphas from 2..max_alpha and a small set of fractional orders near minima.
    """
    best_eps = float("inf")
    best_alpha = None
    # precompute RDP per step
    for alpha in range(2, max_alpha+1):
        rdp = rdp_poisson_subsampled_gaussian(alpha, q, sigma)
        rdp_total = rdp * T
        eps = (rdp_total - math.log(delta)) / (alpha - 1.0)
        if eps < best_eps:
            best_eps = eps
            best_alpha = alpha
    return best_eps, best_alpha

# Small brute-force test to validate log-space calculation for small alphas
def brute_rdp_integer(alpha, q, sigma):
    # direct arithmetic (no log-sum-exp) to compare
    s = 0.0
    for j in range(0, alpha+1):
        comb = math.comb(alpha, j)
        p = comb * (q**j) * ((1 - q)**(alpha - j)) * math.exp((j*j - j) / (2.0 * sigma * sigma))
        s += p
    return math.log(s) / (alpha - 1.0)

def test():
    for alpha in [2,3,4,5,10]:
        q = 0.1
        sigma = 1.0
        v1 = rdp_poisson_subsampled_gaussian(alpha, q, sigma)
        v2 = brute_rdp_integer(alpha, q, sigma)
        assert abs(v1 - v2) < 1e-10, (alpha, v1, v2)
    # small composition test vs naive: one step T=1
    eps, a = get_epsilon(q=0.1, sigma=1.0, T=1, delta=1e-5, max_alpha=50)
    print("Test passed. Example epsilon:", eps, "best alpha:", a)

if __name__ == "__main__":
    test()