Fast Privacy Loss Accounting for Random Allocation via Multiscale FFT

Introduction

We consider privacy accounting for random allocation, a sampling scheme in which the data of each user is used in $k$ steps chosen randomly and uniformly from a sequence of $t$ steps (independently of the other users). This ensures that each user contributes a fixed number of times while still benefiting from the privacy amplification due to sampling. For $k=1$, this is the balls-and-bins sampling scheme used to define the batches in DP-SGD [CGHLKKMSZ24, CCGHST25]. More generally, the scheme can be used to enforce a communication constraint in distributed aggregation by limiting the number of coordinate blocks that each user sends [AFKRT26]. In both of these applications random allocation has been shown to have utility advantages over Poisson sampling [CGHLKKMSZ24, CCGHST25, AFKRT26]. These applications have also motivated increasingly accurate analyses of its privacy guarantees [FS25ra, DCO25, FS26pld]. For the Gaussian mechanism, the most accurate of these analyses shows that the privacy guarantees of random allocation are slightly better than those of Poisson sampling with the same noise scale and expected batch size in nearly all parameter regimes [FS26pld]. Thus in these applications random allocation achieves a better privacy-utility trade-off than Poisson sampling.

As shown in [FS26pld] (FS26 below), the privacy of random allocation can be analyzed using privacy loss distributions (PLDs). Given a PLD of a dominating pair of distributions of the per-step algorithm, their accountant computes a PLD that dominates random allocation applied to that algorithm (an allocation PLD). They reduce the allocation PLD to the distribution of a sum of $t$ independent likelihood ratios (that is, exponentiated privacy losses), and compute this distribution by repeated squaring. Each pairwise sum is computed by a direct convolution on a geometric grid and is rounded to the grid in the direction that preserves stochastic domination. We refer to this algorithm as direct convolution. Its running time is quadratic in the inverse of $\alpha$ the accuracy parameter in $\varepsilon$ that bounds the error in the privacy loss. While direct convolution is efficient for a moderate level of accuracy, this dependence leads to a high running time in the high-accuracy regime. A highly accurate PLD is needed for the approximation error to remain small after a large number of compositions. For example, the analysis of DP-SGD with PREAMBLE in FS26 includes thousands of inner compositions and more than ten thousand training updates. In this regime, computing a PLD for which the final $\varepsilon$ is accurate to within $0.05$ typically takes more than an hour on a single Apple M4 core. Finding the noise scale that achieves a target $\varepsilon$ requires multiple such computations, and a careful choice of the accuracy parameters of direct convolution.

The standard FFT-based approach to the composition of PLDs [KJH20, GLW21] does not resolve this problem. Composition adds privacy losses, whereas random allocation adds their exponentials, that is, likelihood ratios. A uniform grid over likelihood ratios needs to resolve small ratios over a range that can be very large. As a result, ensuring accuracy $\alpha$ using FFT requires an impractically large grid, and smaller grids have been demonstrated to have a poor trade-off between accuracy and running time [FS26pld Appendix D]. Thus, given the practical benefits of random allocation, we address the question of whether there is a more efficient approach to privacy accounting for random allocation.

Our Contribution

Our main contribution is a multiscale FFT algorithm with connect-the-dots projections. This algorithm computes a dominating allocation PLD in time that is nearly linear in the inverse of the accuracy parameter. Our algorithm approximates the PLD of the same dominating pair of distributions as in FS26, accepts a general PLD as its input and outputs a PLD. In particular, the output PLD supports additional accounting operations such as subsampling and composition.

Multiscale FFT.

We replace the geometric grid of FS26 by a grid with $m$ uniformly spaced points in each dyadic interval $[2^j,2^{j+1})$ for $j \in \mathbb{Z}$, and compute every pairwise sum by multiscale FFT (§ 3). Its output is identical to rounding the exact sum onto the grid. For the same accuracy, the dyadic grid has a comparable number of points. Its main advantage is for computation of each merge: direct convolution’s runtime is quadratic in grid size since it sums all pairs of grid points, whereas multiscale FFT takes time nearly linear in the number of grid points. As in FS26, the input PLD is only approximated between its $\beta/2t$ and $1-\beta/2t$ quantiles, where $\beta$ bounds the probability of the tails (translating to an increase in $\delta$). Let $\Delta$ be the length of the smallest interval that contains these quantiles of the log-likelihood ratio $\ln(P(\omega)/Q(\omega))$ (where $(P,Q)$ is the dominating pair that defines the input PLD) both for $\omega\sim P$ and for $\omega\sim Q$, that is, of the input PLD and of the negation of its dual. Then for $t$ steps the running time of our algorithm (up to reading the input; Corollary 4.5 and Remark 4.6) is

For example, for the Gaussian mechanism with noise scale $\sigma$, $\Delta=O(1/\sigma^2+\sqrt{\log(t/\beta)}/\sigma)$. Thus, for a fixed number of steps, the running time is linear in $\Delta/\alpha$ up to a logarithmic factor, whereas the running time of direct convolution is quadratic in the analogous ratio [FS26pld Theorem 4.6].

Connect-the-dots projections.

Conservative rounding in FS26 moves every value in the same direction, and thus increases the expected privacy loss and the $\chi^2$ divergence of the computed pair by amounts of the order of $\alpha$. We refer to the increase in the expected privacy loss as the bias. We replace it by the connect-the-dots (CtD) projection of [DGKKM22], which splits the mass of each value between its two neighboring grid points so as to preserve its mean. We demonstrate that the bias of our algorithm and its error measured as the increase in the $\chi^2$ divergence of the computed pair both drop to $O(\alpha^2/\log(t))$, with the same running time and the same $(\alpha,\beta)$ guarantee (§ 4). Our focus on measuring the increase in $\chi^2$ relies on recent work showing that $\chi^2$ of the dominating pair tightly controls the privacy loss for (subsampled) composition and for random allocation [FS26clt]. Our algorithm applies CtD projection at multiple levels and, as we argue briefly in Remark 4.4, the results in [FS26clt] apply to all but a small number of the highest-level projections.

CtD projection gives valid upper bounds only when it is applied to likelihood ratios under $Q$. Therefore we compute only the add direction and obtain the remove direction as its dual.

Numerical results.

Numerical evaluation using an implementation that can be found at [PLD_acc_msfft] shows that our algorithm is faster than direct convolution of FS26 by several orders of magnitude in the regimes that require high accuracy. As a quick numerical example, in the setting of DP-SGD with PREAMBLE above, multiscale FFT computes the allocation PLD in $0.012$ seconds and, for the same noise scale, improves the upper bound on $\varepsilon$ at $\delta=10^{-6}$ from $0.985$ to $0.92$ (§ 5.0.0.3). Naturally, the speedups are smaller when the input PLD is concentrated, since direct convolution is already fast on the resulting small grid. We report the full range of these comparisons in § 5, together with improved privacy bounds for DP-SGD with PREAMBLE after Poisson subsampling and composition.

We also examine the effect of the CtD projection on the accuracy of the final bound and compare it with the error resulting from the rounding approach in FS26. Our evaluation suggests that the CtD improves the dependence of the running time on the accuracy from about $1/\alpha$ to about $\sqrt{1/\alpha}$ (§ 5). This suggests that the accuracy of the final bound behaves essentially like the bias and $\chi^2$ error.

Technical Overview

We represent distributions of likelihood ratios on a grid that is uniform within each dyadic interval $[2^j,2^{j+1})$ (or band) and whose spacing doubles between consecutive bands (Figure 1). Such a grid gives a relative approximation at every scale while allowing the use of FFT-based convolutions within a band. The main step of our algorithm computes the distribution of the sum of two independent random variables supported on this grid and projects it back onto the same grid using the CtD projection (we refer to this step as a merge).

Computing a merge by direct convolution requires going over all pairs of grid points, whereas FFT requires both values to lie on the same uniform grid. The obvious issue with using FFT on the dyadic grid is that two values from different bands lie on grids with different spacings.

However, the sum of such values does not require the finer resolution of the lower band. Specifically, for a value $x$ in band $j$ and a value $y$ in a lower band, the sum $x+y$ lies in bands $j$ or $j+1$, and the spacing of the grid above $2^j$ is at least the spacing of band $j$. We therefore project $y$ onto the uniform grid with the spacing of band $j$ (extended down to zero), on which $x$ also lies. A crucial question is whether this can be done without increasing the approximation error. We show that the answer is yes as long as the uniform grids of different bands are nested, which holds for the dyadic grid (every grid point in band $j$ or above is a multiple of the spacing of band $j$). In fact, in this case the result is identical to the projection of $x+y$ itself onto the grid.

Thus for every band $j$ we apply this projection to all the mass of each distribution below band $j$. The resulting projected prefix occupies only $m+1$ grid points, regardless of the number of lower bands. All the pairs of values in which the larger value is in band $j$ then lie on the same uniform grid. Their contribution to the sum is given by three convolutions on this grid (band $j$ of both distributions with each other, and band $j$ of each distribution with the projected prefix of the other), which we compute using FFTs of length $4m$. Pairs in which both values lie below band $j$ are accounted for, at a finer resolution, in a lower band. Overall, every pair of values is accounted for exactly once. Therefore a merge computes the projection of the exact sum in time that is nearly linear in the size of the grid (Theorem 3.1).

[x=0.8cm,y=0.72cm,
  knot/.style={circle,fill,inner sep=0pt,minimum size=3.4pt},
  latt/.style={circle,draw,fill=white,inner sep=0pt,minimum size=4.2pt},
  pref/.style={circle,draw,fill=black!40,inner sep=0pt,minimum size=4.2pt},
  lab/.style={font=\scriptsize},
  side/.style={font=\scriptsize,anchor=east}]
\draw[->] (-0.3,1.7) -- (16.7,1.7);
\foreach \x in {0,0.5,0.625,0.75,0.875,1,1.25,1.5,1.75,2,2.5,3,3.5,4,5,6,7,8,10,12,14,16}
  \node[knot] at (\x,1.7) {};
\foreach \x in {0,1,2,4,8,16} \node[lab,below=2pt] at (\x,1.7) {$\x$};
\foreach \a/\b/\j in {0.5/1/-1,1/2/0,2/4/1,4/8/2,8/16/3}{
  \draw[decorate,decoration={brace,amplitude=2.5pt}] (\a,1.88) -- (\b,1.88);
  \node[lab] at ({(\a+\b)/2},2.3) {$\j$};}
\node[side] at (-0.45,2.3) {band};
\node[side] at (-0.45,1.7) {grid $\gridOf{4}$};
\draw[->,thin,shorten <=2pt,shorten >=3pt] (2.5,1.7) -- (2,0);
\draw[->,thin,shorten <=2pt,shorten >=3pt] (2.5,1.7) -- (3,0);
\node[lab] at (2.95,1.0) {$y$};
\foreach \x in {4,5,6,7,8,10,12,14} \draw[densely dashed,gray] (\x,1.12) -- (\x,0.14);
\draw[->] (-0.3,0) -- (16.7,0);
\foreach \x in {0,...,15} \node[latt] at (\x,0) {};
\foreach \x in {0,...,4} \node[pref] at (\x,0) {};
\node[side] at (-0.45,0) {grid $h_2\mathbb Z_{\ge0}$};
\draw[decorate,decoration={brace,mirror,amplitude=3pt,raise=5pt}] (0,0) -- (4,0)
  node[lab,midway,below=8pt,align=center] {mass below band $2$\\ (projected prefixes)};
\draw[decorate,decoration={brace,mirror,amplitude=3pt,raise=5pt}] (4.1,0) -- (14,0)
  node[lab,midway,below=8pt,align=center] {support of the three convolutions at band $2$};

As in FS26, we compute the distribution of the sum of $t$ identically distributed likelihood ratios by repeated squaring, which requires only logarithmically many merges. Each projection changes a value by a factor of at most $1+1/m$, resulting in an error of at most $1/m$ in the log space. Adding nonnegative values preserves a common bound on their multiplicative error. Thus the error of the computed sum grows with the number of merges that a summand passes through (at most $\lceil\log_2(t)\rceil$) and not with the number of summands. In particular, $m=O(\log(t)/\alpha)$ points per band suffice for accuracy $\alpha$. Every distribution in the computation has $O(\Delta+\log(t))$ bands (a merge adds at most one band at the top), and each band is merged in $O(m\log(m))$ time, which gives the running time in (1).

Viewed as a random map, the CtD projection preserves the conditional mean of the likelihood ratio, $\underset{}{\mathbb E}\left[\widetilde Y\mid Y\right]=Y$, where $\widetilde Y$ is the projection of a random variable $Y$. This makes the increase of the $\chi^2$ divergence caused by a projection quadratic in the spacing of the grid, due to the cancellation of the first-order errors (§ 4.0.0.4), as for the mean and the variance of the privacy loss [FS26clt Section 6]. With a spacing of the order of $\alpha/\log(t)$ and $O(\log(t))$ levels of projections, this gives the bound $O(\alpha^2/\log(t))$ on the $\chi^2$ error and on the bias stated above. This is in contrast to the rounding of FS26, which moves every value in the same direction and shifts the expected privacy loss by $O(\alpha)$.

At the same time, the CtD projection moves some of the mass down and need not give stochastic domination. Instead, preservation of the conditional mean increases $Y$ in the convex order, which is preserved under sums of independent random variables. Under $Q$, both privacy profiles of random allocation are expectations of convex functions of the average likelihood ratio, which gives the validity of both directions (Lemma 4.3).

As discussed, we obtain the remove direction as the dual of the add direction. For tightness we show that in both directions the exact and the computed losses can be coupled so that their finite values differ by at most $\gamma=(\lceil\log_2(t)\rceil+1)\ln(1+1/m)$. To construct the coupling, we implement each projection as a random map that moves a value to one of its two neighboring grid points with the probabilities that preserve its mean (§ 4.0.0.3). For an input with an unbounded range of likelihood ratios (such as the Gaussian mechanism) we first truncate its tails, at a negligible cost in the tail parameter, and then project it onto the grid by one more CtD projection (Remark 4.6).

Related Work

The dominating pair of distributions for random allocation with $k=1$ and Gaussian noise was given in the analyses of balls-and-bins sampling for DP-SGD [CGHLKKMSZ24] and for the matrix mechanism [CCGHST25]. It was extended in [FS25ra] to algorithms whose steps are dominated by an arbitrary pair of distributions, and [DCO25] derived it independently for Gaussian noise. [CGHLKKMSZ24] and [CCGHST25] estimate the privacy parameters of this pair using Monte Carlo simulations, and [DG26] base Monte Carlo accounting for $b$-min-sep subsampling (a generalization of balls-and-bins sampling) in banded matrix factorization via dynamic programming. These simulations provide strong evidence that the privacy guarantees of $1$-out-of-$t$ allocation for Gaussian noise are similar to those of $t$ compositions of Poisson sampling with rate $1/t$. The primary issue with Monte Carlo estimation is its computational cost when $\delta$ is small and a high confidence level is required. In addition, it estimates the privacy parameters of the pair rather than its full privacy profile or PLD, and thus does not support subsequent accounting operations such as subsampling and composition.

The analysis in [FS25ra] gives provable bounds on the privacy parameters of $k$-out-of-$t$ allocation via a reduction to a $k$-fold composition of $1$-out-of-$t/k$ allocation (Lemma A.1), including bounds on the Rényi divergence of the dominating pair. [DCO25] and [SK26] give bounds on the Rényi divergence for Gaussian noise and for matrix mechanisms, respectively, the latter computed by dynamic programming. The primary disadvantage of bounds on Rényi divergence is that their conversion to $(\varepsilon,\delta)$ bounds is known to be somewhat lossy. Additionally their bounds apply only to Rényi orders of at least $2$ and therefore provide a loose bound on the expected privacy loss.

Our work builds on the PLD accountant of FS26, which computes the PLD of the dominating pair numerically and thus supports subsequent accounting operations. Specifically, we use its representation of the allocation PLD via sums of likelihood ratios, its repeated squaring schedule (which gives logarithmic dependence on the number of steps), and its accounting operations on PLDs. We improve the cost of each merge at the resolution needed for an accurate PLD, reduce the $\chi^2$ error and the bias by replacing the rounding with the CtD projection, and obtain the remove direction as the dual of the add direction. In the remove direction, FS26 also merge a likelihood ratio under $P$ and thus discretize both the input PLD and its dual. Since stochastic domination is not preserved by the dual, they discretize the dual separately and compute it from the exact input [FS26pld Section 4 and Remark C.2].

Numerical composition of PLDs was introduced by [MM18] and [SMM19], who quantize the PLD into privacy buckets and compose them by direct convolution, using exponentiation by squaring for a logarithmic dependence of the running time on the number of compositions. [KJH20, KJPH21, KH21] compose discretized PLDs using FFT, and [GLW21] give a tight analysis of the dependence of the running time on the accuracy. This approach is now the standard for the analysis of DP-SGD and is supported by several libraries [Google20, Microsoft21, Meta21]. [GKKM22] give a faster self-composition algorithm in which the grid spacing and the truncation interval grow with the number of composed copies, so that each of its stages uses only polylogarithmically many grid points. These algorithms compute the distribution of a sum of privacy losses. In contrast, random allocation requires adding likelihood ratios, and the dyadic grid is designed to address their large dynamic range.

The mean-preserving projection we use is the pessimistic CtD approximation of [DGKKM22], expressed in terms of likelihood ratios (Claim C.2). Viewed as a random map, it is the one-dimensional case of the dual quantization of [PW12], which preserves the conditional mean and therefore gives upper bounds on the expectations of convex functions and preserves the convex order. [PW09] apply dual quantization recursively to the partial sums of a random walk and bound the resulting error for smooth functions.

The closest technique to our multiscale FFT merge is the fast projected convolution of [Hac07], which computes the exact $L^2$-orthogonal projection of the convolution of two piecewise-polynomial functions onto a locally refined grid (for example, one refined towards zero) in nearly linear time, under certain conditions on the grids. As in our algorithm, the grids are decomposed into levels of uniform grids, and the projection is incorporated before the discrete convolutions on each level, which are computed by FFT. Our merge applies this principle to the CtD projection of atomic measures, for which the nesting of the uniform grids of the bands makes the projection of the lower bands exact (Theorem 3.1).

Preliminaries

We follow the notation of FS26. The statements that we use, with the proof of the new identity (4), appear in Appendix § 7. For a random variable $Y$ we denote its distribution by $f_Y$ and its complementary CDF by $\overline F_Y(x)\coloneqq\underset{}{\mathbb P}\left(Y>x\right)$. Random variables take values in $[-\infty,\infty]$. We denote by $\ln$ and $\log_2$ the natural and binary logarithms, and write $\log$ inside $O(\cdot)$, where the base is immaterial.

Random allocation and adjacency.

We consider $t$-step algorithms, which apply an algorithm $M$ to $t$ subsets of the dataset in sequence, where each step also receives the outputs of the previous steps (its view). In the $k$-out-of-$t$ random allocation scheme, each element of the dataset is added to $k$ of the $t$ subsets, chosen uniformly at random and independently of the other elements. We denote the resulting algorithm by $\mathcal A_{t,k}\left(M\right)$, and by $\mathcal A_{t}\left(M\right)$ when $k=1$. To define adjacency for sampling schemes that operate on a fixed number of elements, the domain $\mathcal X$ of the elements is augmented with a null element $\bot$, which $M$ treats as absent, that is, $M((\boldsymbol s,\bot),\boldsymbol v)=M(\boldsymbol s,\boldsymbol v)$ for every dataset $\boldsymbol s$ and view $\boldsymbol v$. Two datasets are adjacent if one of them is obtained from the other by replacing a single element by $\bot$. Thus replacing an element by $\bot$ has the same effect on the output as removing it, while the number of elements, and thus the sampling scheme, remain unchanged.

Privacy loss distributions and privacy profiles.

For distributions $P$ and $Q$ over a domain $\Omega$, the privacy loss random variable $L_{P,Q}$ is $\ln(P(\omega)/Q(\omega))$ for $\omega\sim P$, and its distribution is the privacy loss distribution (PLD) [DR16]. For $\kappa\ge0$, the $\kappa$-hockey-stick divergence between $P$ and $Q$ is $\boldsymbol H_{\kappa}\left(P\,\Vert\,Q\right)\coloneqq\underset{}{\mathbb E}\left[(1-\kappa e^{-L_{P,Q}})_+\right]$, where $(x)_+\coloneqq\max\{0,x\}$ [BKOZB12]. For a random variable $L$ we refer to $\boldsymbol H_{\kappa}\left(L\right)\coloneqq\underset{}{\mathbb E}\left[(1-\kappa e^{-L})_+\right]$ as its hockey-stick functional. The privacy profile $\delta_{M}(\varepsilon)$ of an algorithm $M$ is the largest $e^\varepsilon$-hockey-stick divergence between its output distributions on two adjacent datasets with the same view [BBG18]. Since the hockey-stick divergence is asymmetric, we denote by $\vec{\delta}_{M}$ the privacy profile in the remove direction, in which the second dataset is obtained from the first by replacing an element by $\bot$, and by \reflectbox{$\vec{\reflectbox{$\delta$}}$}_{M} the one in the add direction, in which the first dataset is obtained from the second in this way. The algorithm is $(\varepsilon,\delta)$-differentially private if \delta_{M}(\varepsilon)=\max\{\vec{\delta}_{M}(\varepsilon),\reflectbox{$\vec{\reflectbox{$\delta$}}$}_{M}(\varepsilon)\}\le\delta [DKMMN06].

Dominating pairs.