Composition Accounting with RDP for Spatial Queries

A spatial pipeline does not make one release; it makes a heatmap every hour, an origin–destination matrix every night, and a handful of ad-hoc answers in between. Basic sequential composition adds the ε\varepsilon of each of those, and after a week of ordinary operation the sum is a number no one can defend. Rényi differential privacy (RDP) is the accounting that makes the same workload affordable without weakening the guarantee — it tracks a curve rather than a scalar, and converts to (ε,δ)(\varepsilon, \delta) only at the moment of reporting. This guide implements that accountant for spatial workloads, under privacy budget management in Core Fundamentals & Architecture for Spatial Privacy.

Parameter Configuration and Calibration

  • orders — the Rényi orders tracked. The accountant maintains ρ(α)\rho(\alpha) at a fixed grid of orders and reports the minimum converted ε\varepsilon across them. Too few orders and the reported budget is inflated; too many and every debit costs more arithmetic than the release. A grid of roughly 30 orders from 1.25 to 256, denser at the low end, is the well-trodden choice and is what most libraries ship.
  • delta — fixed once, per privacy unit. The conversion is only meaningful at a stated δ\delta, and δ\delta must be well below the inverse of the population: 10610^{-6} for a metro of a million devices, 10910^{-9} if the same ledger covers a national dataset. Changing δ\delta later invalidates every previously reported number.
  • sensitivity per query family. The accountant is agnostic to what a release is; it only knows the mechanism and its sensitivity. A per-cell count has Δ=1\Delta = 1 under a device-day unit with clipping; a dwell-minute sum has Δ=C\Delta = C, the clipping bound. Register the sensitivity with the query family so a debit cannot be recorded without it.
  • sampling_rate for subsampled releases. Subsampling amplifies privacy, and RDP is where that amplification is cheapest to express. A release computed over a 1% sample of devices costs roughly q2q^2 of the full-sample RDP at small qq, which is the single largest saving available to a high-cadence spatial pipeline.
Mechanism RDP at order α Notes
Gaussian, noise multiplier σ α/(2σ2)\alpha / (2\sigma^2) the canonical case; composes by addition
Subsampled Gaussian, rate q 2q2α/σ2\approx 2 q^2 \alpha / \sigma^2 valid for small q; use a library bound otherwise
Laplace, scale b numerically evaluated no clean closed form; tabulate per α
Post-processing (clamp, smooth) 0 free — never debit it
Why a high-cadence pipeline needs Rényi accounting Two curves of total reported epsilon against the number of releases, on a logarithmic axis. Sequential composition grows linearly and passes 1,600 at 3,200 releases. Renyi accounting over subsampled Gaussian releases grows far more slowly, staying in single digits across the whole range. The gap is the difference between a workload that can be defended and one that cannot. Why a high-cadence pipeline needs Rényi accounting δ = 1e-6, noise multiplier 1.1 10 100 1000 0 500 1000 1500 releases total ε reported sequential (ε = 0.5 each) RDP, subsampled q = 0.01
The same workload, two accountants: one reports an ε nobody can defend, the other reports single digits. Nothing about the mechanism changed.

Reference Implementation

python
from __future__ import annotations

import math
from dataclasses import dataclass, field
from typing import Iterable, Sequence

# A dense-at-the-bottom order grid: small alphas dominate for tight (eps, delta),
# large alphas matter only for very small delta.
DEFAULT_ORDERS: tuple[float, ...] = tuple(
    [1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 16.0,
     20.0, 24.0, 32.0, 48.0, 64.0, 96.0, 128.0, 192.0, 256.0]
)


@dataclass
class RDPAccountant:
    """Tracks the Renyi curve for one privacy unit and converts on demand."""

    orders: Sequence[float] = DEFAULT_ORDERS
    rho: list[float] = field(default_factory=list)

    def __post_init__(self) -> None:
        if not self.rho:
            self.rho = [0.0] * len(self.orders)

    def spend_gaussian(self, sigma: float, *, q: float = 1.0, steps: int = 1) -> None:
        """Debit `steps` Gaussian releases at noise multiplier sigma, sampling rate q.

        The subsampled bound below is the small-q approximation. It is accurate for
        the q < 0.01 regime a device-sampled spatial release runs in, and it is NOT
        safe for q approaching 1 — fall back to a library's exact bound there.
        """
        if not 0 < q <= 1.0:
            raise ValueError("sampling rate must be in (0, 1]")
        if q > 0.1:
            # Full-sample bound: exact, and always an upper bound on the subsampled one.
            per_step = [a / (2.0 * sigma * sigma) for a in self.orders]
        else:
            per_step = [2.0 * q * q * a / (sigma * sigma) for a in self.orders]
        self.rho = [r + steps * p for r, p in zip(self.rho, per_step)]

    def spend_laplace(self, scale_ratio: float, *, steps: int = 1) -> None:
        """Debit Laplace releases. `scale_ratio` is Delta/b, i.e. the usual epsilon.

        The Laplace mechanism's RDP has no compact closed form; this is the standard
        numerically stable expression, evaluated per order.
        """
        e = scale_ratio
        per_step = []
        for a in self.orders:
            if abs(a - 1.0) < 1e-9:
                per_step.append(e + math.expm1(-e))
                continue
            term = (
                a * math.exp((a - 1.0) * e) + (a - 1.0) * math.exp(-a * e)
            ) / (2.0 * a - 1.0)
            per_step.append(max(0.0, math.log(term) / (a - 1.0)))
        self.rho = [r + steps * p for r, p in zip(self.rho, per_step)]

    def epsilon(self, delta: float) -> tuple[float, float]:
        """Convert the curve to (epsilon, delta), returning the minimising order too."""
        if not 0 < delta < 1:
            raise ValueError("delta must be in (0, 1)")
        best, best_alpha = float("inf"), self.orders[0]
        for rho_a, a in zip(self.rho, self.orders):
            if a <= 1.0:
                continue
            eps = rho_a + math.log1p(-1.0 / a) - (math.log(delta) + math.log(a)) / (a - 1.0)
            if eps < best:
                best, best_alpha = eps, a
        return best, best_alpha

    def merge(self, other: "RDPAccountant") -> "RDPAccountant":
        """Compose two curves — e.g. a heatmap ledger and a trajectory ledger."""
        if tuple(self.orders) != tuple(other.orders):
            raise ValueError("cannot merge accountants with different order grids")
        return RDPAccountant(self.orders, [a + b for a, b in zip(self.rho, other.rho)])

Validation Checkpoint

python
def _validate() -> None:
    delta = 1e-6

    # 1. RDP composition must never exceed naive sequential composition.
    naive = RDPAccountant()
    naive.spend_laplace(0.5, steps=20)
    eps_rdp, _ = naive.epsilon(delta)
    assert eps_rdp <= 20 * 0.5 + 1e-9, "RDP must not be worse than sequential"

    # 2. Composition is additive on the curve, not on epsilon.
    a1 = RDPAccountant(); a1.spend_gaussian(1.1, q=0.01, steps=100)
    a2 = RDPAccountant(); a2.spend_gaussian(1.1, q=0.01, steps=100)
    merged = a1.merge(a2)
    both = RDPAccountant(); both.spend_gaussian(1.1, q=0.01, steps=200)
    assert abs(merged.epsilon(delta)[0] - both.epsilon(delta)[0]) < 1e-9

    # 3. Doubling the steps must cost strictly less than doubling epsilon.
    one = RDPAccountant(); one.spend_gaussian(1.1, q=0.01, steps=100)
    two = RDPAccountant(); two.spend_gaussian(1.1, q=0.01, steps=200)
    assert two.epsilon(delta)[0] < 2 * one.epsilon(delta)[0]

    # 4. Subsampling must be cheaper than the full-sample release.
    full = RDPAccountant(); full.spend_gaussian(1.1, q=1.0, steps=100)
    sub = RDPAccountant(); sub.spend_gaussian(1.1, q=0.01, steps=100)
    assert sub.epsilon(delta)[0] < full.epsilon(delta)[0] / 10

    # 5. The minimising order must be interior, not at the grid edge — otherwise
    #    the order grid is too narrow and the reported epsilon is inflated.
    _, alpha = sub.epsilon(delta)
    assert DEFAULT_ORDERS[0] < alpha < DEFAULT_ORDERS[-1], f"alpha at grid edge: {alpha}"

    # 6. A tighter delta costs more, but sublinearly in log(1/delta).
    loose, _ = sub.epsilon(1e-5)
    tight, _ = sub.epsilon(1e-9)
    assert tight > loose and tight < 3 * loose

    print("RDP accountant: all assertions passed")


_validate()
The minimising order moves with the sampling rate Three U-shaped curves of converted epsilon against Renyi order for three sampling rates. Each has a distinct minimum, and the minimum moves to a lower order as the sampling rate rises. A fixed order grid that stops too early would report the value at its edge rather than at the true minimum, inflating the budget. The minimising order moves with the sampling rate 500 releases, σ = 1.1, δ = 1e-6 10 100 0 200 400 Rényi order α converted ε q = 0.001 q = 0.01 q = 0.05
Assert that the minimising α is interior to your grid. When it lands on the edge, the accountant is reporting the best number it can see, not the best available.

Assertion 5 deserves to run on every debit in production, not just in tests. When the minimising order lands on the edge of the grid, the accountant is reporting the best number it can see, which may be far from the best number available — and the symptom is a budget that mysteriously exhausts early.

Incident Response and Edge Cases

  • The reported ε jumps after a routine deploy. Check whether the order grid changed or a mechanism was re-registered with a different sensitivity. RDP curves from different order grids cannot be merged, and a silent grid change makes historical entries incomparable — which is why merge refuses rather than interpolating.
  • A subsampled release is debited at q = 1 by mistake. The budget drains 100× faster than planned. Make the sampling rate a required argument at the call site with no default, so that forgetting it is a TypeError rather than a silent over-charge.
  • Laplace and Gaussian releases in the same ledger. Legitimate and common; the curves simply add. What is not legitimate is converting each to (ε,δ)(\varepsilon, \delta) separately and adding the epsilons — that discards the whole benefit of RDP and typically doubles the reported budget.
  • Someone debits post-processing. Clamping negatives, smoothing, rounding and re-projection are all free. A ledger that charges for them will exhaust early and, worse, will teach the team that the accountant is arbitrary. Make the free operations explicit in the API so nobody has to guess.
  • The curve is stored but the unit is not. A ledger entry without its privacy unit cannot be audited later. Store the unit, the mechanism, the sensitivity and the sampling rate with every debit — the curve alone is not evidence.
Tightening δ is cheap; tightening ε is not A bar chart of reported epsilon across five delta values spanning five orders of magnitude. Epsilon rises only modestly as delta tightens, from about 1.0 at one in ten thousand to about 1.6 at one in a billion. The practical reading is that delta should be set conservatively from the population size and then left alone. Tightening δ is cheap; tightening ε is not 400 subsampled releases at q = 0.008, σ = 1.1 0 1 2 ε at the minimising order 1.03 1e-4 1.20 1e-5 1.35 1e-6 1.50 1e-7 1.74 1e-9
Five orders of magnitude of δ cost about 60% more ε. Set δ well below 1/n, stop negotiating it, and spend the argument on ε.

Frequently Asked Questions

Why track a curve instead of a single epsilon?

Because composition is additive on the Rényi curve and only sub-additive on epsilon. Keeping the curve lets the accountant convert once, at the end, at the order that minimises the reported epsilon — which for a realistic release cadence is often two to five times better than adding epsilons as you go.

Can I convert to (ε, δ) at every release for a dashboard?

You can display it, but do not store it as the ledger value. Conversion is lossy: once you have collapsed the curve to a scalar you cannot compose further without giving up the tightness. Store the curve, convert for display.

Does RDP change the guarantee I can claim?

No. It is an accounting technique, not a different definition. After conversion you claim ordinary (ε, δ)-differential privacy at your stated δ. What changes is that the ε you can honestly claim is smaller for the same workload.

What order grid should I use?

Around 20–30 orders spanning roughly 1.25 to 256, denser below 10. Verify with an assertion that the minimising order is interior; if it repeatedly lands at the top of the grid, extend the grid rather than accepting the inflated number.

Up one level: Privacy Budget Management · Section: Core Fundamentals & Architecture for Spatial Privacy.