Fast Privacy Loss Accounting for Random Allocation via Multiscale FFT

Fast Privacy Loss Accounting for Random Allocation via Multiscale FFT Moshe Shenfeld∗

arXiv:2610.42117v1 [cs.LG] 30 Oct 2026

Vitaly Feldman

Abstract In the random allocation sampling scheme each user’s data is used in k steps chosen randomly and uniformly from a sequence of t steps. This scheme is known to improve on Poisson sampling [Chu+24; CGHSG25] in private optimization and in communication-efficient aggregation [AFKRT26]. The existing privacy loss distribution (PLD) accountant for this scheme runs in time quadratic in the inverse of the accuracy parameter α [FS26]. We give an algorithm, multiscale FFT, that represents likelihood ratios on a grid that is uniform within each dyadic interval and computes the exact rounding of every pairwise sum using FFTs within each interval. It keeps the guarantees of the accountant and improves the dependence of its running time on α from O(α−2 ) to O(α−1 log(1/α)). In addition, we demonstrate that rounding can be performed via the mean-preserving connect-the-dots projection of Doroshenko et al. [DGKKM22]. This reduces the errors that rounding introduces in the mean of the privacy loss and in the χ2 divergence of the computed pair from the order of α to O(α2 / log(t)). In our experiments, our algorithm achieves speedups of several orders of magnitude at high accuracy and gives improved privacy bounds for DP-SGD with PREAMBLE.

1

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 [Chu+24; CGHSG25]. 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 [Chu+24; CGHSG25; AFKRT26]. These applications have also motivated increasingly accurate analyses of its privacy guarantees [FS25; DCO25; FS26]. 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 [FS26]. Thus in these applications random allocation achieves a better privacy-utility trade-off than Poisson sampling. As shown in [FS26] (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 ∗

Work partially done while author was an intern at Apple.

1

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 α – the accuracy parameter in ε 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 ε 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 ε 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 α using FFT requires an impractically large grid, and smaller grids have been demonstrated to have a poor trade-off between accuracy and running time [FS26, 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.

1.1

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 [2j , 2j+1 ) for j ∈ Z, and compute every pairwise sum by multiscale FFT (Section 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 β/2t and 1 − β/2t quantiles, where β bounds the probability of the tails (translating to an increase in δ). Let ∆ be the length of the smallest interval that contains these quantiles of the log-likelihood ratio ln(P (ω)/Q(ω)) (where (P, Q) is the dominating pair that defines the input PLD) both for ω ∼ P and for ω ∼ 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 (∆ + log(t)) log2 (t) log(t) O log α α 

!

.

(1)

For example, for the Gaussian mechanism with noise scale σ, ∆ = O(1/σ 2 + log(t/β)/σ). Thus, for a fixed number of steps, the running time is linear in ∆/α up to a logarithmic factor, whereas the running time of direct convolution is quadratic in the analogous ratio [FS26, Theorem 4.6]. p

2

Connect-the-dots projections. Conservative rounding in FS26 moves every value in the same direction, and thus increases the expected privacy loss and the χ2 divergence of the computed pair by amounts of the order of α. We refer to the increase in the expected privacy loss as the bias. We replace it by the connect-the-dots (CtD) projection of Doroshenko et al. [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 χ2 divergence of the computed pair both drop to O(α2 / log(t)), with the same running time and the same (α, β) guarantee (Section 4). Our focus on measuring the increase in χ2 relies on recent work showing that χ2 of the dominating pair tightly controls the privacy loss for (subsampled) composition and for random allocation [FS26b]. Our algorithm applies CtD projection at multiple levels and, as we argue briefly in Remark 4.4, the results in [FS26b] 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] 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 ε at δ = 10−6 from 0.985 to 0.92 (Section 5). 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 Section 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 p the CtD improves the dependence of the running time on the accuracy from about 1/α to about 1/α (Section 5). This suggests that the accuracy of the final bound behaves essentially like the bias and χ2 error.

1.2

Technical Overview

We represent distributions of likelihood ratios on a grid that is uniform within each dyadic interval [2j , 2j+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 2j 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

3

−1

band grid G4 0

1

0

1 2

y

2

3

4

8

16

grid h2 Z≥0 mass below band 2 (projected prefixes)

support of the three convolutions at band 2

Figure 1: The dyadic grid and the merge step for m = 4. Top: the grid G4 and its bands Bj = [2j , 2j+1 ) (bands −2, −3, . . . accumulate at zero), where band 2 is [4, 8) with spacing h2 = 1. Bottom: the uniform grid of band 2 extended down to zero. A value y below band 2 is projected onto this grid (arrows), after which the mass below band 2 occupies m + 1 points. The knots of G4 in [4, 16) (dashed, in bands 2 and 3) are a subset of this grid, and [4, 16) contains the support of the three convolutions at band 2 (band 2 of both inputs with each other, and band 2 of each input with the projected prefix of the other). Thus projecting the sums onto G4 gives the same measure as projecting the exact sums (Theorem 3.1). 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). 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 ⌈log2 (t)⌉) and not with the number of summands. In particular, m = O(log(t)/α) points per band suffice for accuracy α. Every distribution in the computation has O(∆ + 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 i random map, the CtD projection preserves the conditional mean of the likelihood h as a ratio, E Ye | Y = Y , where Ye is the projection of a random variable Y . This makes the increase of the χ2 divergence caused by a projection quadratic in the spacing of the grid, due to the cancellation of the first-order errors (Section 4), as for the mean and the variance of the privacy loss [FS26b, Section 6]. With a spacing of the order of α/ log(t) and O(log(t)) levels of projections, this gives the bound O(α2 / log(t)) on the χ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(α). 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

4

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 γ = (⌈log2 (t)⌉ + 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 (Section 4). 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).

1.3

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 [Chu+24] and for the matrix mechanism [CGHSG25]. It was extended in [FS25] to algorithms whose steps are dominated by an arbitrary pair of distributions, and Dong, Chen, and Ozgur [DCO25] derived it independently for Gaussian noise. Chua et al. [Chu+24] and Choquette-Choo et al. [CGHSG25] estimate the privacy parameters of this pair using Monte Carlo simulations, and Dong and Ganesh [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 δ 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 [FS25] 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. Dong, Chen, and Ozgur [DCO25] and Schuchardt and Kalinin [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 (ε, δ) 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 χ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 [FS26, Section 4 and Remark C.2]. Numerical composition of PLDs was introduced by Meiser and Mohammadi [MM18] and Sommer, Meiser, and Mohammadi [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. Koskela, Jälkö, and Honkela [KJH20], Koskela et al. [KJPH21], 5

and Koskela and Honkela [KH21] compose discretized PLDs using FFT, and Gopi, Lee, and Wutschitz [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 [Goo22; Mic21; Met21]. Ghazi et al. [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 Doroshenko et al. [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 Pagès and Wilbertz [PW12b], which preserves the conditional mean and therefore gives upper bounds on the expectations of convex functions and preserves the convex order. Pagès and Wilbertz [PW12a] 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 Hackbusch [Hac07], which computes the exact L2 -orthogonal projection of the convolution of two piecewisepolynomial 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).

2

Preliminaries

We follow the notation of FS26. The statements that we use, with the proof of the new identity (4), appear in Appendix A. For a random variable Y we denote its distribution by fY and its complementary CDF by F Y (x) := P (Y > x). Random variables take values in [−∞, ∞]. We denote by ln and log2 the natural and binary logarithms, and write log inside O(·), 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 At,k (M ), and by At (M ) when k = 1. To define adjacency for sampling schemes that operate on a fixed number of elements, the domain X of the elements is augmented with a null element ⊥, which M treats as absent, that is, M ((s, ⊥), v) = M (s, v) for every dataset s and view v. Two datasets are adjacent if one of them is obtained from the other by replacing a single element by ⊥. Thus replacing an element by ⊥ 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 Ω, the privacy loss random variable LP,Q is ln(P (ω)/Q(ω)) for ω ∼ P , and its distribution is the privacy loss distribution (PLD) [DR16]. For κ ≥ 0, the κ-hockey-stick divergence between P and Q h i is H κ (P ∥ Q) := E (1 − κe−LP,Q )+ , where (x)+ := max{0, x} [BKOZ12]. For a random variable L h

i

we refer to H κ (L) := E (1 − κe−L )+ as its hockey-stick functional. The privacy profile δM (ε) of

6

an algorithm M is the largest eε -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 ⃗δM the privacy profile in the remove direction, in which the second dataset is obtained from the first by replacing an element by ⊥, and by δ⃗M the one in the add direction, in which the first dataset is obtained from the second in this way. The algorithm is (ε, δ)-differentially private if δM (ε) = max{⃗δM (ε), δ⃗M (ε)} ≤ δ [DKMMN06]. Dominating pairs. We say that a PLD LP1 ,Q1 dominates a PLD LP2 ,Q2 if H κ (P2 ∥ Q2 ) ≤ H κ (P1 ∥ Q1 ) for all κ ≥ 0. A pair (P, Q) is a dominating pair for M in the remove (add) direction if ⃗δM (ε) ≤ H eε (P ∥ Q) (respectively, δ⃗M (ε) ≤ H eε (P ∥ Q)) for all ε [ZDW22]. Domination is transitive, and it is preserved under composition, whose PLD is the sum of the PLDs of the components [ZDW22, Theorem 10]. If every step of a t-step algorithm M is dominated by (P, Q) in the remove direction, then At (M ) is dominated by (P̄t , Qt ) in the remove direction and by (Qt , P̄t ) in the add direction, where t 1X Qi−1 × P × Qt−i , P̄t := t i=1 and the k-out-of-t allocation reduces to compositions of such allocations (Lemma A.1). We refer to the PLDs LP̄t ,Qt and LQt ,P̄t as the allocation PLDs. h

i

PLD realizations. A random variable L is a PLD realization if E e−L ≤ 1 and P (L = −∞) = 0 [FS26, Definition 3.1]. Every PLD is a PLD realization, and conversely, every PLD realization is the PLD of some pair of distributions. The dual D (L) of a PLD realization L is the PLD of the reversed pair, so that D (LP,Q )h= LQ,P i (A.2). Its distribution is that of −L with probabilities weighted by −L −L e , with mass 1 − E e at +∞ [FS26, Definition 3.1]. A PLD realization U dominates a PLD realization L if H κ (L) ≤ H κ (U ) for all κ ≥ 0, and in this case D (U ) dominates D (L) (A.3). The privacy profiles of L and of D (L) are both expectations of functions of e−L : h

i

h

i

h

i

H eε (D (L)) = 1 − E e−L + E (e−L − eε )+ .

H eε (L) = E (1 − eε−L )+ ,

(2)

The first formulah is the i definition, and the second one h holds since the i mass h of D (L) iat +∞ −L −L ε L contributes 1 − E e and its finite part contributes E e (1 − e e )+ = E (e−L − eε )+ . Both functions are convex in e−L . For the PLD L = LQ,P in the add direction, e−L = P (ω)/Q(ω) for ω ∼ Q is the likelihood ratio under Q. The allocation PLDs. ⃗t (L) := ln ψ

For a PLD realization L let   1

t

L(0)

e

+

t−1 X

e

−D(L)(i)



,

i=1

t  X

1 (i) ψ⃗t (L) := − ln e−L , t i=1 

(3)

⃗t (LP,Q ) and L t = where the superscripts denote independent copies. FS26 show that LP̄t ,Qt = ψ Q ,P̄t ψ⃗t (LQ,P ), and that these transformations preserve domination (Theorem A.5 and [FS26, Section 4]). ⃗t , the factor eL(0) is the likelihood ratio under P and the other factors are likelihood ratios In ψ under Q. Thus FS26 compute the remove direction from both the input PLD and its dual. However, the two allocation PLDs are the PLDs of the same pair in the two directions. Therefore, for every PLD realization L (A.7),   ⃗t (D (L)) = D ψ⃗ (L) . ψ (4) t 7

(i)

Thus both allocation PLDs are determined by the average Rt := 1t ti=1 e−L of t likelihood ratios under Q, since ψ⃗t (L) = − ln(Rt ). By (2), both privacy profiles of random allocation are expectations of convex functions of Rt . P

Stochastic domination and tightness. For any two random variables V, U , we write V ⪯ U if F V ≤ F U , and V ⪯(α,β) U if F V (x) ≤ F U (x − α) + β for all x. Following FS26 we say that U (α, β)-tightly stochastically dominates V if V ⪯ U ⪯(α,β) V . Stochastic domination implies domination, and V ⪯(α,β) U implies H eε (V ) ≤ H eε−α (U ) + β (A.4). The output of our accountant dominates the exact PLD but need not stochastically dominate it. Therefore we state its guarantees directly in terms of privacy profiles. Definition 2.1 (Tight domination). Given α ≥ 0 and β ∈ [0, 1], a PLD realization U (α, β)-tightly dominates a PLD realization L if for every ε ∈ R, H eε (L) ≤ H eε (U ) ≤ H eε−α (L) + β. That is, U dominates L, and the privacy profile of U is at most that of L shifted by α, plus β. By the implication above, if U dominates L and U ⪯(α,β) L then U (α, β)-tightly dominates L, and in particular tight stochastic domination implies tight domination with the same parameters (A.4). h

i

Definition 2.2 (χ2 error and bias). For a PLD realization L we denote χ2 (L) := E eL − 1 ∈ [0, ∞]. For PLD realizations U and L with χ2 (L) < h∞, the χ2ierrorh of U relative to L is χ2 (U ) − χ2 (L). i The bias of U relative to L is the increase E U 1{U <∞} − E L1{L<∞} in the expectation of the finite part of the loss. h

i

h

i

Note that χ2 (L) ≥ 0, since E eL ≥ 1/E e−L ≥ 1 by Jensen’s inequality, and that χ2 (L) = ∞ if L has mass at +∞. For a pair (P, Q) with P ≪ Q, χ2 (LP,Q ) = E

ω∼Q

(P (ω)/Q(ω))2 − 1 is the χ2





divergence of P from Q. The χ2 divergence of a single step controls the Gaussian approximation of the privacy losses of random allocation [FS26b, Section 5].

3

Multiscale FFT

In this section we describe our accountant (Algorithm 1) and the multiscale FFT merge that it relies on. The accountant computes the add direction from the distribution of the sum of t independent likelihood ratios under Q, and it obtains the remove direction as the dual of the result (by (4)). After projecting the likelihood ratio of its input onto a grid that is uniform within each dyadic interval, it then computes the distribution of the sum by repeated squaring, in which every convolution is followed by the projection onto the grid (a merge). We define the grid and the projection in Section 3.1 and describe the accountant in Section 3.2. In Section 3.3 we show that the multiscale FFT merge computes every merge exactly in time nearly linear in the size of the grid (Theorem 3.1). We analyze the accountant in Section 4.

3.1

The Grid and the Projection

For an integer m ≥ 2, let 

Gm := {0} ∪ 2

j



i−1 1+ m 8





: j ∈ Z, i ∈ [m] .

We refer to its points as knots, to Bj := [2j , 2j+1 ) as band j, to hj := 2j /m as the spacing of band j, and to hj Z≥0 as the lattice of band j. The knots in band j are m consecutive points of its lattice (Figure 1). The analysis uses two properties of the grid (Claim B.2): (G1) Resolution: Every x > 0 satisfies u ≤ x ≤ v for two consecutive knots 0 < u < v ≤ (1 + 1/m)u. (G2) Dyadic scaling: For every j ∈ Z, every knot in band j or above is an integer multiple of hj . The multiscale FFT merge relies on the scaling, and the analysis of the accuracy in Section 4 on the resolution. Our algorithm relies on CtD projections [DGKKM22] onto Gm , which split the mass of each value between its two neighboring knots so as to preserve its mean. Formally, for a set of knots S ⊆ [0, ∞) that contains 0 and a value x between adjacent knots u < v of S, let ΠS (δx ) :=

v−x x−u δu + δv , v−u v−u

and let ΠS (δx ) := δx for x ∈ S ∪ {∞}. We extend ΠS linearly to finite measures and abbreviate Π := ΠGm and Πh := ΠhZ≥0 for h > 0 (Definition B.1). Properties of the projections. Our analysis uses only the following properties of the projections. We state them, and Theorem 3.1, for general projections, since it is possible to use the multiscale FFT merge with rounding to get lower bounds (Remark 3.2). A projection onto a set of knots S is a map from finite measures on [0, ∞] to finite measures on S ∪ {∞} that fixes every measure supported on S ∪ {∞}. Let ξ and ζ be finite measures on [0, ∞], h > 0 and j ∈ Z. (P1) Locality: ΠS (δx ) is supported on {u, v} for every x between adjacent knots u < v of S. (P2) Linearity: ΠS (ξ + ζ) = ΠS (ξ) + ΠS (ζ). (P3) Commutativity: Πh (ξ ∗ ζ) = Πh (ξ) ∗ ζ if ζ is supported on hZ≥0 . (P4) Projection hierarchy: Π2h (Πh (ξ)) = Π2h (ξ), and Π(Πhj (ξ)) = Π(ξ) if ξ is supported on [2j , ∞). (P5) Mean preservation: ΠS (δx ) has mean x for every x ≥ 0. The multiscale FFT merge relies on (P1)–(P4), and the analysis of the accountant in Section 4 on (P1) and (P5). The CtD projections satisfy all five properties (Lemma B.3). The second part of the projection hierarchy is specific to the dyadic grid and relies on its scaling and the fact 2j ∈ Gm .

3.2

The Accountant

⃗ in the Algorithm 1 receives a PLD realization L in the add direction. (For a PLD realization L   ⃗ , which only reweights its atoms.) It chooses the number m of remove direction we use L = D L points per band from the accuracy parameter (line 1; see Corollary 4.5) and projects e−L onto Gm , ′ which gives a PLD realization L′ for which e−L is supported on Gm (lines 2–3). It then computes b which approximates the sum of t by repeated squaring the distribution fSb of the computed sum S,   b ⃗t := − ln(S/t) ⃗ t := D L ⃗t independent copies of e−L (line 4). Finally, it returns L and its dual L ⃗t has mass f (z) at − ln(z/t), and L ⃗ t has mass z f (z) at ln(z/t) for every z > 0 (line 5). That is, L b b t S S h

i

b and the remaining mass 1 − E S/t at +∞.

9

Algorithm 1 CtD accountant rand-alloc-ctd(L, t, α) Input: L a PLD realization (add direction) with bounded finite values; t ∈ N, t ≥ 2; α ∈ (0, 1] ⃗t (D (L)) ⃗t and L ⃗ t that (α, 0)-tightly dominate ψ⃗t (L) and ψ Output: L 1: m ← ⌈(⌈log2 (t)⌉ + 1)/α⌉ 2: j− ← ⌊ℓ− / ln(2)⌋; j+ ← ⌊ℓ+ / ln(2)⌋ ′ 3: L′ ← dyadic-projection(L; j− , j+ , m) {supp(e−L ) ⊆ Gm } 4: f b ← self-conv(fe−L′ , t) S   b ⃗t ← − ln(S/t); ⃗t ⃗t ← D L 5: L L 6: return

⃗t , L ⃗t L

Discretizing the input. In line 2 of Algorithm 1, ℓ− and ℓ+ are the infimum and the supremum of the finite values of −L (an input with unbounded losses is first truncated, as described in Remark 4.6). Since band j contains ex for j = ⌊x/ ln(2)⌋, the indices j− and j+ are the bands that contain eℓ− and eℓ+ , the infimum and the supremum of the positive values of the likelihood ratio e−L , and j+ − j− ≤ (ℓ+ − ℓ− )/ ln(2) + 1. The range [2j− , 2j+ +1 ], from the start of band j− to the end of band j+ , contains every positive value of e−L , and the discretization dyadicprojection(L; j− , j+ , m) returns the PLD realization L′ for which fe−L′ = Π(fe−L ), that is, the CtD projection of the likelihood ratio onto Gm . It applies to discrete and continuous inputs alike and needs only the probability and the first moment of e−L on each cell between adjacent knots. Repeated squaring. As in FS26, we compute the n-fold self-convolution of a distribution µ by repeated squaring, and replace every convolution by a merge, which we compute by the P multiscale FFT merge ms-fft-merge (Algorithm 2). Let n = sk=1 2bk with b1 < · · · < bs , let µ(0) := µ and µ(b) := ms-fft-merge(µ(b−1) , µ(b−1) ) for 1 ≤ b ≤ bs , and let ξ1 := µ(b1 ) and ξk := ms-fft-merge(ξk−1 , µ(bk ) ) for 2 ≤ k ≤ s. We denote the output ξs by self-conv(µ, n). Repeated squaring uses bs + s − 1 ≤ 2⌈log2 (n)⌉ merges. Since a convolution is the distribution of a sum of independent random variables, the computation corresponds to a binary tree with n leaves, which we refer to as the unfolded tree. Its leaves are independent random variables with distribution µ, and every internal node is a merge applied to the values computed at its two children (a merge of µ(b−1) with itself is a node whose children are two independent copies of the subtree of µ(b−1) ). The depth of this tree is ⌈log2 (n)⌉ (Claim B.4). Unlike FS26, we do not truncate the tails between merges, since a merge adds at most one band (Theorem 3.1).

3.3

The Multiscale FFT Merge

We now describe the multiscale FFT merge (Algorithm 2) and show that it is exact, that is, that it computes the projection onto Gm of the exact convolution of its inputs. It accounts for every pair of values in the band of the larger value, and it replaces the mass of each input below that band by its projection onto the lattice of the band. This projection does not change the result (Theorem 3.1), and it reduces the contribution of each band to three convolutions on its lattice, which we compute by FFTs. For a probability distribution µ on Gm ∪ {∞}, let µj be its restriction to Bj , µ<j its restriction to (0, 2j ), µ+ its positive finite part, and µ(0) and µ(∞) its masses at 0 and at ∞. We refer to µ̂<j := Πhj (µ<j ) as its projected prefix. By locality (P1), µ̂<j is supported on {0, hj , . . . , mhj }, regardless of the number of bands below j, since it moves every value in (0, 2j ) = (0, mhj ) to the neighboring knots of hj Z≥0 . Algorithm 2 computes the projection of the convolution µ ∗ ν of two

10

such distributions band by band, and it maintains the projected prefixes of the current band in µ̂ and ν̂ (line 5). Algorithm 2 Multiscale FFT merge ms-fft-merge(µ, ν) max Input: µ, ν probability distributions on Gm ∪ {∞}, with µ+ , ν+ supported on jj=j Bj min Output: Π(µ ∗ ν) 1: ξ ← 0; µ̂ ← 0; ν̂ ← 0 2: for j = jmin , . . . , jmax do 3: C ← FFTconv(µj , νj ) + FFTconv(µj , ν̂) + FFTconv(µ̂, νj ) {FFTs of length 4m} 4: ξ ← ξ + Π(C) 5: µ̂ ← Πhj+1 (µ̂ + µj ); ν̂ ← Πhj+1 (ν̂ + νj ) 6: end for
7: return ξ + µ(0) ν+ + ν(0) µ+ + µ(0)ν(0) δ0 + 1 − (1 − µ(∞))(1 − ν(∞)) δ∞

S

Theorem 3.1 (Multiscale FFT merge). Let Π and Πh for h > 0 be projections onto Gm and onto hZ≥0 that satisfy properties (P1)–(P4), let µ and ν be probability distributions on Gm ∪ {∞} whose positive finite parts are supported on the bands jmin , . . . , jmax , and let J := jmax − jmin + 1. Then Algorithm 2 returns Π(µ ∗ ν), uses O(Jm log(m)) arithmetic operations and O(Jm) space, and the positive finite part of its output is supported on the bands jmin , . . . , jmax + 1. Proof outline. The full proof appears in Appendix B.2. Every pair of positive finite values is counted in exactly one term of µ+ ∗ ν+ =

jX max

µj ∗ νj + µj ∗ ν<j + µ<j ∗ νj ,

j=jmin

namely, in the term of the larger of their bands. Since µj ∗ ν<j is supported on [2j , ∞) and µj on hj Z≥0 , by the dyadic scaling, (P4)

Π(µj ∗ ν<j ) = Π Πhj (µj ∗ ν<j )

(P3)

= Π(µj ∗ ν̂<j ),

and similarly for µ<j ∗ νj . By the linearity of the projection, this gives Π(µ ∗ ν) =

jX max

Π µj ∗ νj + µj ∗ ν̂<j + µ̂<j ∗ νj

(5)

j=jmin

+ µ(0) ν+ + ν(0) µ+ + µ(0)ν(0) δ0 + 1 − (1 − µ(∞))(1 − ν(∞)) δ∞ . Its second line collects the terms with the atoms at 0 and at ∞. By linearity and the projection hierarchy, µ̂<j+1 = Πhj+1 (µ̂<j + µj ), which line 5 computes. Finally, since the projected prefixes are supported on {0, hj , . . . , mhj }, in units of hj the convolutions in line 3 are supported on the indices m, . . . , 4m − 2, and thus FFTs of length 4m compute them exactly. Together with locality for the projections in lines 4 and 5, this gives the bounds on the running time, the space and the support of the output. Remark 3.2 (Multiscale FFT merge with rounding). The projections that round every value up, or every value down, to a neighboring knot (Definition B.1) also satisfy (P1)–(P4) (Lemma B.3). Thus the multiscale FFT merge also computes the exact rounding of the sum, with the same bounds on 11

its running time, its space and the support of its output (Theorem 3.1). Rounding moves every value in one direction instead of preserving its mean, which preserves stochastic domination in that direction. For upper bounds, the sums are rounded up in the remove direction, whose loss is b b ln(S/t), and down in the add direction, whose loss is − ln(S/t), and the input is discretized in the matching directions. For lower bounds, which the CtD projection does not give, all these directions are reversed. Thus replacing the direct convolutions in the accountant of FS26 by multiscale FFT merges with the same rounding keeps the guarantees of FS26 (Remark 4.2) and makes the running time nearly linear in 1/α. We use this accountant for comparisons in our experiments (Section 5).

4

Analysis

We now state the guarantees of Algorithm 1. ⃗t and L ⃗ t be the outputs of Algorithm 1 on an input L, let m, Theorem 4.1 (CtD accountant). Let L j− and j+ be the resolution and the lowest and the highest bands that it computes, let J := j+ −j− +1, ⃗t and L ⃗ t are PLD realizations such that: and let γ := (⌈log2 (t)⌉ + 1) ln(1 + 1/m). Then L ⃗t (D (L)). ⃗t dominates ψ⃗t (L), and L ⃗ t dominates ψ 1. Validity: L 2. Tightness: These dominations are (γ, 0)-tight. 3. χ2 error and bias: If m ≥ ⌈log2 (t)⌉ + 1, then in each direction in which the allocation PLD Ψ that the output dominates satisfies χ2 (Ψ) < ∞, the χ2 error of the output relative to Ψ is between 0 and (⌈log2 (t)⌉ + 1)(1 + χ2 (Ψ))/(4m2 ). In each direction, the output and the allocation PLD that it dominates have the same mass at infinity, and the bias of the output relative to this PLD is between 0 and (⌈log2 (t)⌉ + 1) ln2 (1 + 1/m)/8. 4. Computational complexity: The running time is O((J + log(t)) m log(m) log(t)), apart from lines 2 and 3, whose cost depends on the representation of the input, and the space is O((J + log(t)) m). We outline the proofs of the first three parts below and give the full proof in Appendix C.2. The bounds on the running time and the space follow from Theorem 3.1, since the positive finite part of every distribution in the computation is supported on at most J + ⌈log2 (t)⌉ + 1 bands (Claim C.5). The proofs do not use the choice of m in line 1. Thus the theorem holds for every m ≥ 2. Remark 4.2 (Accountant with upper and lower rounding). The guarantees of Theorem 4.1, except for the bounds on the χ2 error and the bias, also hold for the accountant with rounding of Remark 3.2, which computes each direction separately, as in FS26. Validity and tightness then hold with γ replaced by (⌈log2 (t − 1)⌉ + 2) ln(1 + 1/m) in the remove direction. With the directions of the rounding reversed, the same holds for lower bounds, with each output dominated by the allocation PLD instead of dominating it, except that the outputs need not be PLD realizations. The outputs can still be used to provide valid lower bounds on the privacy profile as in FS26. Thus we also obtain a method to compute lower bounds in time that is nearly quadratically better than that of FS26. The proofs are minor variants of the proof of Theorem 4.1 and we omit them for brevity. Validity uses the monotonicity of rounding instead of the convex order, and tightness results from the grid resolution and the locality of the projection.

12

The projection as a random map. In the analysis we view the CtD projection Π onto Gm as a random map that is applied independently of everything else and maps a value x to a random knot x̃ with the distribution Π(δx ). The image Ye of a random variable Y then has the distribution Π(fY ). The proofs of validity, tightness, and the bounds on the χ2 error and the bias use only the fact that the merges are exact (Theorem 3.1), the grid resolution, and the locality (P1) and the mean preservation (P5) of the projection. For the random map they give x̃ ∈ {u, v} for every value x between adjacent knots u < v, and E [x̃] = x for every x ≥ 0. Validity. Validity follows from the convex order, where Y ≤cx Y ′ if E [φ(Y )] ≤ E [φ(Y ′ )] for every convex function φ [SS07]. By Jensen’s inequality, Y ≤cx Y ′ whenever E [Y ′ | Y ] = Y , and two random variables that are ordered in the convex order have the same mean, since x 7→ x and x 7→ −x are convex. (i)

Lemma 4.3 (Validity). Let L be a PLD realization for which e−L is bounded, let Rt := 1t ti=1 e−L , b be a bounded nonnegative random variable such that Rt ≤cx R. b Then − ln(R) b is a PLD and let R ⃗t (D (L)). realization that dominates ψ⃗t (L), and its dual dominates ψ P

h

i

b they have the same mean E e−L ≤ 1, and since R b is bounded, − ln(R) b is a Proof. Since Rt ≤cx R, PLD realization. By (2), both privacy profiles are expectations of convex functions of the likelihood b gives domination in the add direction. Domination ratio, and thus the convex order between Rt and R in the remove direction follows by (4), since the dual preserves domination (Claim A.3).

We apply Lemma 4.3 to the input L. We view the projection of the input in line 3 as projections at the leaves of the unfolded tree, and every projection as an independent random map (by Claim B.5, this gives the distributions that repeated squaring computes). By the mean preservation of every b := S/t b satisfies projection, that the computed average R h i induction over the unfolded tree shows Pt (i) 1 −L b b E R | Rt = Rt for the exact average Rt := t i=1 e (Lemma C.3). Thus Rt ≤cx R, and since 



b and L ⃗t = − ln(R) ⃗t = D L ⃗t , Lemma 4.3 gives the first part of Theorem 4.1. L This argument does not apply to the likelihood ratio under P . In the remove direction of FS26, the likelihood ratio x under P enters the privacy profile through the function x 7→ (1 − eε t/(x + s))+ , where s is the sum of the other factors. This function is strictly concave where it is positive, since its second derivative there is −2eε t/(x + s)3 . Thus a mean-preserving projection of a value x between two adjacent knots in this region strictly decreases the privacy profile.

Tightness. Tightness bounds the privacy profiles of the output from above, and we obtain it from the same coupling. By locality of the projection and the grid resolution, every projection, including b ≤ eγ Rt , those at the leaves, changes a value by a factor of at most 1 + 1/m, and thus e−γ Rt ≤ R since the depth of the tree is ⌈log2 (t)⌉ (Claim B.4 and Lemma C.3). In the add direction, the b are either both infinite or both finite, and in ⃗t = − ln(R) coupled losses ψ⃗t (L) = − ln(Rt ) and L the latter case they differ by at most γ (Lemma C.4). The tightness of the add direction does not b necessarily imply the tightness of its dual. In the remove direction, we h i weight this coupling by R b | Rt = Rt , this gives a coupling and add the remaining mass 1 − E [Rt ] at (+∞, +∞). Since E R of the remove losses, whose finite values again differ by at most γ (Lemma C.4). These couplings ⃗t (D (L)). Together with validity, Claim A.4 gives the second ⃗t ⪯(γ,0) ψ⃗t (L) and L ⃗ t ⪯(γ,0) ψ give L part of Theorem 4.1.

13

χ2 error and bias. Rounding moves every value in the same direction, and thus increases the χ2 divergence (Claim C.6) and the expected privacy loss by amounts of the first order in the spacing of the grid. In contrast, by locality of the projection and the grid resolution, a CtD projection changes a positive value by a factor of at most h1 + 1/m, i and together with mean preservation, for e every nonnegative random variable Y we get E Y | Y = Y and (1 + 1/m)−1 Y ≤ Ye ≤ (1 + 1/m)Y (Claim C.1). Thus the first-order effects of a projection cancel, and its χ2 error and bias are of the second order in the spacing, as are its effects on the mean and the variance of the privacy loss [FS26b, Section 6]. In Appendices C.4 and C.5 we show that these errors add up over the at most ⌈log2 (t)⌉ + 1 projections on a path from a leaf to the root of the unfolded tree, which gives the bounds in the third part of Theorem 4.1. Remark 4.4 (Effect of the lower levels on ε). The loss parameter γ of Theorem 4.1 charges ln(1+1/m) to every level of projections. For all but the several highest levels this appears to be overly conservative. At height j, the CtD projections are applied to about t/2j independent sums of 2j likelihood ratios, and these projected sums are then added up. If the higher levels added them exactly, the result would be random allocation over t/2j steps, with the distribution of the projected sums in place of that of the exact ones. The privacy losses of random allocation over many steps are approximately Gaussian, with variance equal to the χ2 divergence of the allocation PLD and mean equal to half of it [FS26b, Theorem 5.2]. Thus the projections at such a level affect ε essentially only through the relative increase in this χ2 divergence, which is of the order of α2 /χ2 (Ψ) for the allocation PLD Ψ. Only the several highest levels, at which few sums are added up, are not covered by this argument. Guarantees for a given accuracy. Substituting the resolution chosen in line 1 of Algorithm 1 into Theorem 4.1 gives the following guarantees in terms of α (Appendix C.6). Corollary 4.5 (CtD accountant for a given accuracy). Let α ∈ (0, 1] and let t ≥ 2 be an integer. ⃗ and L ⃗ be PLD realizations with a total of S atoms, and let ∆ be the width of the smallest Let L ⃗ and of −L. ⃗ There exists an algorithm that receives α, t, interval that contains the finite values of L ⃗t (L) ⃗ ⃗ ⃗ ⃗ ⃗ and ψ⃗t (L), ⃗ in L and L, and returns PLD realizations Lt and Lt that (α, 0)-tightly dominate ψ time  ! (∆ + log(t)) log2 (t) log(t) O S+ log , α α and with space O(S + (∆ + log(t)) log(t)/α). Moreover, in each direction the bias of the output relative to the allocation PLD Ψ that it dominates is at most α2 /(8(⌈log2 (t)⌉+1)), and if χ2 (Ψ) < ∞, then the χ2 error of the output relative to Ψ is at most α2 (1 + χ2 (Ψ))/(4(⌈log2 (t)⌉ + 1)). Remark 4.6 (Truncation of unbounded inputs). An input whose loss is unbounded, such as the PLD of the Gaussian mechanism, is first truncated, separately from the accountant. The truncation chooses a range of likelihood ratios that leaves out at most β/(2t) of them under Q and under P on each side, projects the likelihood ratios below this range onto zero and its lower end by CtD, and moves those above it to its upper end, so that the truncated loss and its dual still dominate those of the input. Each factor of (3) changes only on an event of probability at most β/t, and thus the truncation affects the final bounds only through the tail parameter β, which is negligible when β is much smaller than the values of δ of interest. For the Gaussian mechanism with noise scale σ, the corresponding range of losses has width 1/σ 2 + 2Φ−1 (1 − β/(2t))/σ, since −L and D (L) are distributed as N (∓1/(2σ 2 ), 1/σ 2 ), and the projection in line 3 needs only O(Jm) evaluations of the standard normal CDF Φ, since the probability and the first moment of the likelihood ratio on 14

=1, N=1

Running time (s)

103 102 101 100 10 1 10 2

103 102 101 100 10 1 10 2 10 3

=1, N=100 103 102 101 100 10 1 10 2

10 1

10 3

10 5

10 7

=4, N=1

10 9

10 11

100

10 2

103

10 4

10 6

10 8

10 7

10 9

=4, N=100

102 101 100 10 1 10 2 10 1

10 3

10 5

10 7

10 9

_ref

Direct convolution, rounding (t=10) Multiscale FFT, CtD (t=10)

10 11

101

10 13

10 1

10 3

Direct convolution, rounding (t=1000) Multiscale FFT, CtD (t=1000)

10 5

_ref

Reference resolution (t=10) Reference resolution (t=1000)

Figure 2: Running time against the discrepancy ε − εref of the computed upper bound ε at δ = 10−6 from a reference εref , for direct convolution with rounding and multiscale FFT with CtD projections, for t ∈ {10, 1000}, a single allocation (left) and its 100-fold composition (right), with σ = 1 (top) and σ = 4 (bottom). Here εref is the bound of multiscale FFT with at least twice the largest m used in the same setting. The vertical lines mark the resolution of εref for each t, its difference from the bound with half its m. Each curve ends once its bound comes within twice this resolution of εref . A bound within this resolution of εref , or below it, is drawn on the line. each cell are differences of two values of Φ. Thus the accountant computes PLD realizations that (α, β)-tightly dominate the allocation PLDs of the Gaussian mechanism using log2 (t) 1 O log(t) + 2 + α σ

p

!

log(t/β) log(t) log σ α 

!

arithmetic operations and evaluations of Φ.

5

Numerical Results

We evaluate our accountant (Algorithm 1), referred to as multiscale FFT, using an implementation in the PLD accounting library [PLD], which also contains the implementation of FS26. All the running times are of single-threaded runs on one core of a MacBook Pro M4 Max, and include the composition, when applicable. Running time and accuracy. The two algorithms use their accuracy parameters differently, and CtD projections give only upper bounds. Therefore we compare the algorithms by the discrepancy of their upper bounds from a reference computed by multiscale FFT on a finer grid (Figure 2). The 15

=1, N=1

103

=1, N=100 103 102 101

102

Running time (s)

101 100 10 1 10 2

100 10 1 10 2 10 1

10 3

103 102 101 100 10 1 10 2

10 5

10 7

=4, N=1

10 9

100

10 2

103

10 4

10 6

10 8

10 7

10 9

=4, N=100

102 101 100 10 1 10 2 10 1

10 3

10 5

10 7

_ref

Direct convolution, rounding Direct convolution, CtD

101

10 9

Multiscale FFT, CtD Multiscale FFT, rounding

10 1

10 3

10 5

_ref

Reference resolution (t=1000)

Figure 3: The comparison of Figure 2 at t = 1000 for the four combinations of the convolution method (direct convolution on the geometric grid or multiscale FFT on the dyadic grid) and the projection (rounding or CtD). The discrepancy ε − εref is as in Figure 2, and the dotted line marks the resolution of εref , its difference from the bound with half its m. regimes in the figure cover a wide range of resulting ε values at δ = 10−6 , from ≈ 0.027 (when σ = 4 and t = 1000, for a single allocation) to ≈ 26 (when σ = 1 and t = 10, for the 100-fold composition). At the finest accuracy of direct convolution, where a run takes about 1,000 seconds, multiscale FFT reaches the same discrepancy in at most 0.03 seconds, that is, it is faster by a factor of 104 to 105 . The running time grows roughly as the square of the inverse discrepancy for direct convolution and as its square root for multiscale FFT. Fitting T = c + a xk , where T is the running time and x the inverse of the discrepancy, by least squares on ln(T ) gives exponents k between 1.96 and 2.02 for direct convolution and between 0.44 and 0.55 for multiscale FFT, over the points of the eight settings in Figure 2. Grid and projection. These speedups correspond to the complete accountant, which changes both the grid and the projection. Figure 3 separates the two changes at t = 1000 by combining each convolution method (direct on a geometric grid or multiscale-FFT on the dyadic) with each projection (stochastic domination rounding or hockey-stick domination via connect-the-dots). For each projection, the dyadic grid with multiscale FFT merges reaches the same or better accuracy faster than the geometric grid with direct convolution: at the finest accuracy of the geometric grid that the dyadic grid reaches, it is faster by two to three orders of magnitude. For each grid, CtD projections reduce the discrepancy by orders of magnitude compared with rounding with the same accuracy parameters, and the difference grows with the number of compositions, which is consistent

16

Batch Size: 512

2.2

2.2

2.0

2.0

1.8 1.6 1.4 1.2

1.8 1.6 1.4 1.2

1.0

1.0 101

102

B

101

103

Batch Size: 4096

2.4 2.2

2.2

2.0

2.0

1.8 1.6 1.4 1.2

102

B

103

Batch Size: 600000

2.4

sigma ratio

sigma ratio

Batch Size: 1024

2.4

sigma ratio

sigma ratio

2.4

1.8 1.6 1.4 1.2

1.0

1.0 101

102

B

101

103

PREAMBLE RDP (FS25) PREAMBLE PLD (direct convolution, FS26)

102

B

103

PREAMBLE PLD (multiscale FFT) Gaussian

Figure 4: The ratio between the noise scales required for (1, 10−6 )-DP by DP-SGD with PREAMBLE and by DP-SGD with Gaussian noise addition, as a function of the block size B, in the four settings of FS26. The results for RDP [FS25] and for direct convolution are from FS26. with the χ2 error of rounding accumulating over the compositions. Thus the grid and the merge determine the running time, and the projection determines the χ2 error. These routes also differ in other implementation choices, so this comparison indicates, but does not isolate, the contribution of each change. DP-SGD with PREAMBLE. We follow the setting of FS26 for DP-SGD with Poisson subsampling of the users and PREAMBLE aggregation [AFKRT26], with n = 6 · 105 users, dimension d = 220 , communication budget C = 215 and 10 epochs. Each selected user sends k = C/B out of t = d/B blocks of size B. The analysis reduces to 1-out-of-32 allocation composed k times (up to 8192 times for B = 4), followed by Poisson subsampling and composition over the updates. Since the allocation PLD is composed thousands of times, it has to be computed to a high accuracy. For B = 4 and expected batch size 1024, the complete analysis at a fixed noise scale takes 0.06 seconds with multiscale FFT and about an hour with direct convolution. Accordingly, finding the noise scale that ensures (1, 10−6 )-DP takes at most a few minutes with multiscale FFT, instead of several hours with direct convolution. Figure 4 shows these noise scales, relative to that of DP-SGD in which each user sends the entire gradient, for the expected batch sizes 512, 1024, 4096 and 600000 and each B. The resulting bounds are also slightly tighter: multiscale FFT never requires a larger noise scale than direct convolution, and requires a smaller one in 35 of the 40 settings, mostly for small B, where the number of inner compositions is large. The example in the introduction uses B = 4, expected batch size 1024, the noise scale that direct convolution finds for (1, 10−6 )-DP and 17

multiscale FFT with m = 1024. Multiscale FFT then computes the allocation PLD in 0.012 seconds and gives the upper bound 0.92094 on ε at δ = 10−6 , compared with 0.985 for direct convolution. The direct-convolution column is the run of FS26, with accuracy parameter 0.005, tail parameter 10−8 and at most 106 grid points. FS26 used the expected batch size 1028 instead of 1024, so for batch size 1024 we ran their implementation with the same configuration. Multiscale FFT was run with accuracy parameter 0.01, tail parameter 10−8 and at most 108 points in the output grid. After the division of the accuracy parameter among the compositions, the cap m ≤ 218 binds for B ≤ 512, that is, in 32 of the 40 settings. In our runs, the allocation PLD is composed k times by FFT, subsampled by the PLD subsampling of FS26, and composed over the updates by dp_accounting. The noise scale is found by a bisection that accepts any noise scale for which the bound on ε lies in [0.98, 1]. The running times at a fixed noise scale above were measured in the same way for both algorithms, whereas the search times are indicative only: the searches of multiscale FFT ran four at a time, and those of direct convolution are from its batch-1024 run.

6

Discussion

Multiscale FFT makes the running time of computing the allocation PLD nearly linear in the inverse accuracy, with the same (α, β) bound on the error of the privacy profiles as FS26. Empirically, the accuracy is even better: the error of the final bound on ε behaves essentially like the bias and the χ2 error, which are of the second order in the spacing of the grid (Section 5). As we argue in Remark 4.4, the projections at all but the several highest levels of the unfolded tree affect ε essentially only through the relative increase in the χ2 divergence. A natural question is whether this also holds for the highest levels, at which few independent sums are projected.

AI Usage The authors formulated the problem and proposed the two key elements of the algorithm. The development of these ideas and writing were performed through an iterative process in which AI (Claude Opus 5.5 and ChatGPT Astra 6) generated the implementation, details of proofs and wrote parts of the prose based on author-provided guidance and structure. The resulting text was also extensively edited and proofread by the authors. The implementation of the algorithm and evaluation code was written by AI under close supervision and auditing by the authors. The authors take full responsibility for the content of this paper.

Acknowledgments Shenfeld’s work was supported in part by ERC grant 101125913, Simons Foundation Collaboration 733792, Israel Science Foundation (ISF) grant 2861/20, and Apple. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

18

A

Preliminary Statements

In this section we give the statements from FS26 and ZDW22 that we use, following their notation, together with the proof of (4).

A.1

Domination

The pair (N (1, σ 2 ), N (0, σ 2 )) dominates the Gaussian mechanism with sensitivity one and noise scale σ in both directions, and its PLD is N (1/(2σ 2 ), 1/σ 2 ) [ZDW22]. We say that a pair (P, Q) dominates a t-step algorithm M if it is a dominating pair in the remove direction for every step of M and every view (in the terms of [FS26, Definition 2.6], M is dominated by a randomizer with output distributions P and Q). Lemma A.1 ([FS26, Lemmas 2.7 and 2.8 and Claim 2.9]). Let (P, Q) dominate a t-step algorithm M . Then At (M ) is dominated by (P̄t , Qt ) in the remove direction and by (Qt , P̄t ) in the add mf mc × P̄⌈t/k⌉ , Qt ) in the remove direction and by direction. Moreover, At,k (M ) is dominated by (P̄⌊t/k⌋ the reversed pair in the add direction, where mf + mc = k, mf ⌊t/k⌋ + mc ⌈t/k⌉ = t, and powers denote product distributions. Claim A.2 ([FS26, Section 3 and Claims A.1 and A.2]). Every PLD LP,Q is a PLD realization, and its dual is LQ,P . Conversely, a PLD realization L is the PLD of the pair consisting of the distribution of L and the distribution of −D (L), and D (L) is a PLD realization and the PLD of the reversed pair. Claim A.3 ([ZDW22, Lemma 28]). If a PLD realization U dominates a PLD realization L, then D (U ) dominates D (L). Claim A.4 ([FS26, Claim 4.2]). If V ⪯(α,β) U then H eε (V ) ≤ H eε−α (U ) + β for every ε. In particular, stochastic domination implies domination.

A.2

Random Allocation and Duality

⃗t (L) and ψ⃗ (L) be as in (3). Theorem A.5 ([FS26, Theorem 4.4]). For a PLD realization L, let ψ t ⃗ Then for any two distributions P and Q, LP̄t ,Qt = ψt (LP,Q ) and LQt ,P̄t = ψ⃗t (LQ,P ). Lemma A.6 ([FS26, Lemma 4.5]). Let V1 , . . . , Vn be independent, let U1 , . . . , Un be independent, P P and let Vi ⪯(α,βi ) Ui for every i. Then ln( i eVi ) ⪯(α,P βi ) ln( i eUi ). i

FS26 state this lemma for two summands, and the general case follows by induction. The same proof applies when α = 0 or βi = 0. For α = βi = 0 it states that stochastic domination is preserved P by ln( i eVi ). The normalization by t in (3) is a common shift. 



⃗t (D (L)) = D ψ⃗ (L) . Claim A.7. For every PLD realization L and integer t ≥ 1, ψ t Proof. By Claim A.2, L is the PLD of a pair of distributions, that is, L = LQ,P for some P and ⃗t (D (L)) = L Q, and then D (L) = LP,Q . By Theorem A.5, ψ⃗t (L) = LQt ,P̄t and ψ P̄t ,Qt , which is the dual of LQt ,P̄t .

19

B

Missing Details from Section 3

B.1

Projections and Their Lattice Identities

Definition B.1 (Projections). Let S ⊆ [0, ∞) be a set of knots that contains 0 and has no accumulation point in (0, ∞), such as Gm or the lattice hZ≥0 for h > 0. For x ∈ [0, sup S], let u ≤ x ≤ v be the largest knot that is at most x and the smallest knot that is at least x (so that u = v = x if x ∈ S). For u < v we define Π↑S (δx ) := δv ,

Π↓S (δx ) := δu ,

Πctd S (δx ) :=

v−x x−u δu + δv , v−u v−u

and for x ∈ S all three are equal to δx . We extend each projection linearly to finite measures on [0, sup S], and to measures with an atom at +∞ by leaving that atom unchanged. Claim B.2. For every integer m ≥ 2, the grid Gm satisfies (G1) and (G2). Proof. Let x > 0 and let j ∈ Z be such that x ∈ Bj . Then u ≤ x ≤ v for two consecutive knots u < v with u ∈ Bj and v = u + hj , also when v = 2j+1 . Since u ≥ 2j = mhj , we have v/u = 1 + hj /u ≤ 1 + 1/m, which proves (G1). For (G2), every knot that is at least 2j is of the ′ ′ form 2j (1 + i/m) with j ′ ≥ j and 0 ≤ i ≤ m − 1, and it is equal to 2j −j (m + i)hj . Lemma B.3. Given ρ ∈ {↑, ↓, ctd}, the maps Π := ΠρGm and Πh := ΠρhZ≥0 for h > 0 are projections in the sense of Section 3.1, and for every h > 0, j ∈ Z and all finite measures ξ and ζ, each of these projections ΠS satisfies: ΠS (δx ) is supported on {u, v} for every x between adjacent knots u < v of S (P1), ΠS (ξ + ζ) = ΠS (ξ) + ΠS (ζ) (P2), Πh (ξ ∗ ζ) = Πh (ξ) ∗ ζ if ζ is supported on hZ≥0 (P3), and Π2h (Πh (ξ)) = Π2h (ξ) and Π(Πhj (ξ)) = Π(ξ) if ξ is supported on [2j , ∞) (P4). Moreover, for ρ = ctd, ΠS (δx ) has mean x for every x ≥ 0 (P5). The second part of (P4) does not hold for measures that are not supported on [2j , ∞), since below 2j the knots of Gm are finer than the lattice hj Z. Proof. Each map ΠρS sends a finite measure to a measure on S ∪ {∞} and fixes the knots and +∞. Since it is extended linearly, it satisfies (P2) and fixes every measure supported on S ∪ {∞}, that is, it is a projection. By definition, it satisfies (P1). For ρ = ctd and x between adjacent knots u < v of S, the mean of ΠS (δx ) is u(v − x)/(v − u) + v(x − u)/(v − u) = x, which is (P5). For a set of knots S, a knot r ∈ S and a projection ρ, let kSρ (x, r) := ΠρS (δx )({r}) denote the mass that the projection of x assigns to r. As a function of x, kSctd (·, r) is affine between every two adjacent knots of S, kS↑ (·, r) is constant on every interval (s, s′ ] between adjacent knots, and kS↓ (·, r) is constant on every interval [s, s′ ) between adjacent knots. We first show that if T ⊆ S are two sets of knots then ΠT (ΠS (ξ)) = ΠT (ξ) for every finite measure ξ on which both projections are defined. By (P2) it suffices to consider ξ = δx . If x ∈ S the claim is trivial. Otherwise, let x lie between adjacent knots u < v of S. Since T ⊆ S, the interval [u, v] is contained in an interval between adjacent knots of T . For ρ = ctd and r ∈ T , ΠT (ΠS (δx ))({r}) =

v−x x−u kT (u, r) + kT (v, r) = kT (x, r), v−u v−u

since kT (·, r) is affine on [u, v]. For ρ = ↑ we have ΠS (δx ) = δv with x ∈ (u, v], and kT (v, r) = kT (x, r) since kT (·, r) is constant on an interval that contains (u, v]. The case ρ = ↓ is symmetric. The first part of (P4) follows by taking T = 2hZ≥0 ⊆ S = hZ≥0 . For its second part, let := T Gm ∩ [2j , ∞) and S := hj Z ∩ [2j , ∞), so that T ⊆ S by (G2). Since 2j ∈ hj Z, the projection 20

onto hj Z≥0 of a value that is at least 2j is supported on [2j , ∞) and is equal to its projection onto S. Similarly, since 2j ∈ Gm , the projection onto Gm of a measure supported on [2j , ∞) is equal to its projection onto T . Thus Π(Πhj (ξ)) = ΠT (ΠS (ξ)) = ΠT (ξ) = Π(ξ). For (P3), for y ∈ hZ≥0 the adjacent knots of hZ≥0 around x + y are the adjacent knots around x translated by y, and the weights (for CtD) or the choice of the knot (for rounding) are the same. Thus Πh (δx+y ) = Πh (δx )∗δy , and integrating over ξ(dx) ζ(dy) gives the claim. Atoms at +∞ are not affected by any of the operations.

B.2

The Multiscale FFT Merge

Proof of Theorem 3.1. The identity. Writing µ = µ(0) δ0 + µ+ + µ(∞) δ∞ and similarly for ν, we get µ ∗ ν = µ+ ∗ ν+ + ρ, where ρ denotes the second line of (5). Since µ+ and ν+ are supported on Gm , the measure ρ is supported on Gm ∪ {∞}, and thus Π(ρ) = ρ. Every pair of values from the positive bands i and k is counted in exactly one term of jX max

µ+ ∗ ν+ =

µj ∗ νj + µj ∗ ν<j + µ<j ∗ νj ,

(6)

j=jmin

namely, in the term with j = max{i, k}. For every j, the measure µj ∗ ν<j is supported on [2j , ∞), and µj is supported on hj Z≥0 by (G2). Thus (P4) and (P3) with ζ = µj give (P4)

Π(µj ∗ ν<j ) = Π Πhj (µj ∗ ν<j )

(P3)

= Π µj ∗ Πhj (ν<j ) = Π(µj ∗ ν̂<j ),

and similarly Π(µ<j ∗ νj ) = Π(µ̂<j ∗ νj ). Therefore (P2)

Π(µ ∗ ν) = Π(µ+ ∗ ν+ ) + Π(ρ) jX max

(a)

=

Π(µj ∗ νj ) + Π(µj ∗ ν<j ) + Π(µ<j ∗ νj ) + ρ

j=jmin jX max

(b)

=

Π(µj ∗ νj ) + Π(µj ∗ ν̂<j ) + Π(µ̂<j ∗ νj ) + ρ

j=jmin (P2)

=

jX max

Π µj ∗ νj + µj ∗ ν̂<j + µ̂<j ∗ νj + ρ,

j=jmin

where (a) follows from (6), (P2) and Π(ρ) = ρ, and (b) from the two identities above. This proves (5). The algorithm. We show by induction on j that µ̂ = µ̂<j and ν̂ = ν̂<j at the beginning of iteration j of Algorithm 2. For j = jmin we have µ<j = 0, and thus µ̂<j = Πhj (0) = 0 = µ̂, since a projection fixes the zero measure. If µ̂ = µ̂<j at the beginning of iteration j, then line 5 sets µ̂ to (P2)

Πhj+1 (µ̂<j + µj ) = Πhj+1 Πhj (µ<j ) + Πhj+1 (µj ) (P4)

(P2)

= Πhj+1 (µ<j ) + Πhj+1 (µj ) = Πhj+1 (µ<j+1 ) = µ̂<j+1 ,

where (P4) applies since hj+1 = 2hj , and the last two steps use µ<j+1 = µ<j + µj . The same holds for ν̂. Thus in iteration j, line 3 computes the argument of the jth projection in (5), and line 4 adds this projection to ξ. Therefore Algorithm 2 returns the right-hand side of (5), which is equal to Π(µ ∗ ν) by the identity. 21

Running time, space and support. By (P1), each projection moves a value between adjacent knots to one or both of these knots. In iteration j we represent measures supported on hj Z ∩ [0, 4mhj ) by the arrays of their masses at the indices 0, . . . , 4m − 1 in units of hj . By the definition of Gm , the knots in band j are the indices m, . . . , 2m − 1, and those in band j + 1 are the even indices 2m, . . . , 4m − 2, since hj+1 = 2hj . Since mhj = 2j is a knot of hj Z≥0 , the projection µ̂<j of the values of µ in (0, 2j ) is supported on the indices 0, . . . , m, and so is ν̂<j . By the previous part, these are the values of µ̂ and ν̂ in iteration j. Thus the measure C in line 3 is supported on the indices m, . . . , 4m − 2: the term µj ∗ νj on the indices 2m, . . . , 4m − 2, and the other two terms on the indices m, . . . , 3m − 1. Since 4m − 2 < 4m, the cyclic convolutions of length 4m of the zero-padded arrays are equal to the linear ones, and thus C is computed by a constant number of FFTs of length 4m, that is, in O(m log(m)) operations. The masses of C at the indices m, . . . , 2m − 1 and at the even indices 2m, . . . , 4m − 2 lie on knots of Gm and are fixed by Π, and every odd index i with 2m < i < 4m − 2 lies between the adjacent knots i − 1 and i + 1. Thus line 4 obtains Π(C) by moving (for rounding) or splitting (for CtD) the mass at each such index to these two indices, in O(m) operations, and Π(C) is supported on the bands j and j + 1. Similarly, µ̂ + µj is supported on the indices 0, . . . , 2m − 1, and line 5 obtains its projection onto hj+1 Z≥0 , that is, onto the even indices, by moving or splitting the mass at each odd index to its two neighbors, in O(m) operations. The result is supported on the indices 0, . . . , m in units of hj+1 , as needed for the next iteration. Line 7 takes O(Jm) operations, since µ+ and ν+ are supported on the Jm knots of J bands. Summing over the J iterations, Algorithm 2 uses O(Jm log(m)) operations. The positive finite part of its output is the sum of the measures Π(C) over the J iterations and of µ(0) ν+ + ν(0) µ+ , and thus it is supported on the bands jmin , . . . , jmax + 1. Finally, the inputs, ξ and the output are supported on at most (J + 1)m + 2 points, and the prefixes and the buffers of the FFT take O(m) space, which gives the space bound O(Jm).

B.3

The Unfolded Tree of Repeated Squaring

By Lemma A.6, the loss parameter of a sum is the larger of the loss parameters of its summands. Thus the errors of the merges accumulate along the depth of the unfolded tree (Section 3.2) rather than over all of its merges. The same observation applies to the accountant of FS26, whose analysis charges its spacing α′ for the discretization of the input and for every merge. Its remove direction merges the sum of t − 1 copies with one more factor, so its unfolded tree has depth ⌈log2 (t − 1)⌉ + 1 ≤ ⌈log2 (t)⌉ + 1, and its guarantees also hold with α′ = α/(⌈log2 (t)⌉ + 2) in place of α/(2⌈log2 (t)⌉ + 1). Claim B.4. For n ≥ 1, repeated squaring for the sum of n copies uses at most 2⌈log2 (n)⌉ merges, and its unfolded tree has depth ⌈log2 (n)⌉. Proof. The case n = 1 is trivial, so let n ≥ 2 and let n = sk=1 2bk with b1 < · · · < bs , so that bs = ⌊log2 (n)⌋. Repeated squaring performs bs squarings and s − 1 further merges. Since b1 , . . . , bs are distinct integers in {0, . . . , bs }, we have s ≤ bs + 1, and thus bs + s − 1 ≤ 2⌊log2 (n)⌋ ≤ 2⌈log2 (n)⌉. Repeated squaring computes the sum of 2b copies for b ≤ bs by successive merges of a sum with itself. Thus the unfolded tree of the sum of 2b copies is a complete binary tree of depth b. The accumulated sum is initialized to the sum of 2b1 copies (of depth b1 ) and, for k = 2, . . . , s, it is merged with the sum of 2bk copies. After this merge its depth is max{Dk−1 , bk } + 1 = bk + 1, where Dk−1 ≤ bk−1 + 1 ≤ bk is its depth before the merge. Thus the depth of the result is bs = log2 (n) if s = 1 and bs + 1 = ⌊log2 (n)⌋ + 1 otherwise, which is ⌈log2 (n)⌉ in both cases. P

The next claim justifies viewing every projection in the unfolded tree as an independent random map, as the analysis in Section 4 does. 22

Claim B.5. Let µ be a probability distribution on Gm whose positive values lie in finitely many bands, and let n ≥ 1. Consider the unfolded tree of repeated squaring for the sum of n copies, whose leaves are independent random variables with distribution µ, and in which every internal node applies the CtD projection onto Gm , as an independent random map, to the sum of the values of its children. Then the value at every node has the distribution that repeated squaring computes at the corresponding merge, and in particular the value at the root has the distribution self-conv(µ, n). Proof. We use induction on the height of the node. A leaf has the distribution µ = µ(0) . Consider an internal node whose children have, by the induction hypothesis, the distributions ξ1 and ξ2 computed by repeated squaring. The subtrees of the two children have disjoint sets of leaves and of projections. Thus the values at the children are independent, and their sum has the distribution ξ1 ∗ ξ2 . The image of an independent random variable with distribution f under the random map has the distribution Π(f ) (Claim C.1). Therefore the value at the node has the distribution Π(ξ1 ∗ ξ2 ), which is equal to ms-fft-merge(ξ1 , ξ2 ) by Theorem 3.1 and Lemma B.3, that is, to the distribution computed at the corresponding merge.

C

Missing Details from Section 4

C.1

The CtD Projection and the Pessimistic CtD Approximation

We first prove the properties of the CtD projection as a random map that Section 4 states. Claim C.1. Let x̃ denote the image of a value x ≥ 0 under the CtD projection onto Gm , viewed as a random map. Then x̃ has the distribution Π(δx ), and the image of a random variable with distribution f under an independent random map has the distribution Π(f ). Moreover, for every nonnegative random variable Y and its image Ye under an independent random map, h

i

E Ye | Y = Y, 





1+

 1 1 −1 Y ≤ Ye ≤ 1 + Y, m m

(7)

h i

and in particular P Ye = 0 = P (Y = 0) and E Ye = E [Y ]. Proof. If x is a knot then x̃ = x. Otherwise, let u < x < v be the adjacent knots around x. The map moves x to v with probability (x − u)/(v − u) and to u otherwise, which gives the distribution R Π(δx ). The distribution of the image of a random variable with distribution f is Π(δx ) f (dx), which is Π(f ) by the linear extension of Π. Since x̃ ∈ {u, v} by (P1) and v ≤ (1 + 1/m)u by (G1) (Claim B.2), we have x/(1 + 1/m) ≤ x̃ ≤ (1 + 1/m)x for every x ≥ 0. Conditioning on Y and e applying this bound and (P5) with x = Y gives h(7). i In particular, Y = 0 if and only if Y = 0, and taking expectations in the first identity gives E Ye = E [Y ]. Next we show that the CtD projection of likelihood ratios gives the pessimistic CtD approximation of Doroshenko et al. [DGKKM22]. Given a PLD and a grid of losses, this approximation is the PLD supported on the grid and on +∞ whose privacy profile, as a function of eε , is equal to that of the given PLD at the points of the grid and is affine between them. Claim C.2. Let S ⊆ [0, ∞) be a set of knots that contains 0, let L be a PLD realization such e e be the PLD realization for which e−L that e−L ≤ smax for a knot smax of S, and let L is the CtD −L ε projection of e onto S. Then, as functions of e :

23

 

e is equal to that of D (L) at every point eε ∈ S, is affine between 1. The privacy profile of D L every two consecutive points of S (including 0, where both profiles tend to 1), and is constant for eε ≥ smax . e is equal to that of L at every point eε = 1/s with s ∈ S \ {0}, is affine 2. The privacy profile of L between every two such consecutive points and for eε ≤ 1/smax , and is constant for eε ≥ 1/smin if S has a smallest positive knot smin .  

e are the pessimistic CtD approximations of D (L) and of L on the grids of losses e and L Thus D L ln(S \ {0}) and − ln(S \ {0}), respectively. h i

e e := e−L e = E [X]. Let C(s) := 1 − E [X] + E [(X − s)+ ], Proof. Let X := e−L and X , so that E X  

e By (2), the privacy profiles of D (L) and of D L e and let Ce be defined in the same way for X.

e ε ). For a knot s ∈ S, the function x 7→ (x − s)+ is affine between every at eε are C(eε ) and C(e e two adjacent knots, and the CtD projection preserves the mean of every value. Thus C(s) = C(s). e Between two consecutive knots u < v, every value of X is either at most u or at least v. Thus Ce is e e e ≤ smax . affine on [u, v]. In addition, C(0) = C(0) = 1, and C(s) = 1 − E [X] for s ≥ smax since X This proves the first part.
For the second part, (1 − eε x)+ = eε (x − e−ε )+ − x + e−ε . Thus by (2) the privacy profile of L e and C. e Thus the two privacy profiles are equal at eε is 1 − eε + eε C(e−ε ), and the same holds for L ε e at every point e = 1/s with s ∈ S \ {0}. If C(s) = A + Bs for s = e−ε between two consecutive e is (1 + B) + (A − 1)eε , which is affine in eε . For eε ≤ 1/smax it knots, then the privacy profile of L e is 1 − E [X] eε . If S has a smallest positive knot smin , then Ce is affine on [0, smin ] with C(0) = 1, ε that is, A = 1 there, and the privacy profile is constant for e ≥ 1/smin .   e is supported on − ln(S \ {0}) ∪ {∞} and, by the definition of the dual, D L e is Since L supported on ln(S \ {0}) ∪ {∞}, the two parts of the claim identify them with the pessimistic CtD approximations defined above (see also [FS26, Remark C.2]).

In particular, the CtD approximations of L and of D (L) obtained in this way are duals of each other. FS26 require this property when they use CtD to discretize the input of the remove direction [FS26, Remark C.2].

C.2

Couplings and the Proof of Theorem 4.1

Lemma C.3 (Coupling along the tree). Consider an unfolded tree of depth d whose leaves are independent nonnegative random variables, in which every leaf applies the CtD projection onto Gm , as an independent random map, to its value, and every internal node applies it to the sum of the values of its children. Let L denote the values of the leaves, let Y be their sum and let Yb be the value at the root. Then h

i

E Yb | L = Y,



1+

 1 d+1 1 −(d+1) Y ≤ Yb ≤ 1 + Y, m m

(8)

and Yb = 0 if and only if Y = 0. Proof. We prove the claims for every node n by induction on its height h(n), where Yn and Ybn denoteh the exact and the computed sums at n. For a leaf, Ybn is the projection of Yn , and i by (7), E Ybn | L = Yn and Ybn is within a factor of 1 + 1/m of Yn . For an internal node with 24

children c1 and c2 , the input of the projection is yn = Ybc1 + Ybc2 . By the induction hypothesis, h i h(n) b E [yn | L] = Yc1 +Yc2 = Yn and yn is within a factor of (1+1/m) of Yn . By (7), E Yn | yn , L = yn and Ybn is within a factor of 1 + 1/m of yn , and thus within a factor of (1 + 1/m)h(n)+1 of Yn . This gives (8) at the root. Finally, the projection maps zero to itself and positive values to positive values (Claim C.1). Next we show how to convert a coupling of the likelihood ratios under Q into couplings of the PLDs in both directions. Note that this conversion does not introduce a multiplicative factor in the error in the remove direction. Lemma C.4 (Two directions). Let Y and Yb be bounded nonnegative random variables on a common probability space such that E [Y ] ≤ 1 and h

i

E Yb | Y = Y,

e−γ Y ≤ Yb ≤ eγ Y.

b = (− ln(Y ), − ln(Yb )) and for Then − ln(Y ) and − ln(Yb ) are PLD realizations, and for (A, A)  b = (D (− ln(Y )) , D − ln(Yb ) ) there is a coupling of A and A b under which A = ∞ if and only (A, A) b ≤ γ otherwise. In particular, A ⪯ b b if Ab = ∞, and |A − A| (γ,0) A and A ⪯(γ,0) A. h i

Proof. Since E Yb = E [Y ] ≤ 1, both are PLD realizations. Under a coupling as in the statement, Ab ≤ A + γ and A ≤ Ab + γ, which gives the last statement. We now construct the couplings. Let π be the joint distribution of (Y, Yb ). For − ln(Y ) and − ln(Yb ), the coupling is the joint distribution of − ln(Y ) and − ln(Yb ) under π. Indeed, Y = 0 if and only hif Yb =i0, and the finite values differ by at most γ. For the duals, we weight π(dy, dŷ) by ŷ. Since E Yb | Y = Y , the marginal of the resulting measure over y is the distribution of Y weighted by y, and its marginal over ŷ is the distribution of Yb weighted by ŷ. This measure has mass E [Y ], and under it both coordinates are positive. Taking the logarithm of both coordinates and adding mass 1 − E [Y ] at (∞, ∞), we obtain a coupling of the two duals (by the definition of the dual) in which the finite values differ by at most γ. Claim C.5. In the setting of Theorem 4.1, the positive finite part of the distribution computed at a node of height h of the unfolded tree of repeated squaring is supported on the bands j− , . . . , j+ + h + 1. In particular, every distribution in the computation is supported on at most J + ⌈log2 (t)⌉ + 1 bands. ′

Proof. We use induction on h. For a leaf (h = 0) this holds since the positive values of e−L are knots of Gm in [2j− , 2j+ +1 ] (Section 3.2), and 2j+ +1 is the smallest knot of band j+ + 1. The children of a node of height h ≥ 1 have heights at most h − 1, and thus by the induction hypothesis their positive finite parts are supported on the bands j− , . . . , j+ + h. By Theorem 3.1, the positive finite part of their merge is supported on the bands j− , . . . , j+ + h + 1. The last statement holds since h ≤ ⌈log2 (t)⌉ (Claim B.4) and J = j+ − j− + 1. Proof of Theorem 4.1. Setup. By line 2, the positive values of e−L lie in [2j− , 2j+ +1 ], and thus e−L is a bounded nonnegative random variable. Consider the unfolded tree of repeated squaring for the sum of t copies, with independent copies of e−L at its leaves, in which every leaf projects its value onto Gm (line 3) and every projection is applied as an independent random map. The projected ′ leaves are independent copies of e−L , whose positive values are knots of Gm in [2j− , 2j+ +1 ]. The depth of the tree is d := ⌈log2 (t)⌉ (Claim B.4), and the value Sb at its root has the distribution P (i) fSb computed in line 4 of Algorithm 1 (Claim B.5). Let Rt := 1t ti=1 e−L be the average of 25





b := S/t. b b and L b , and by (3) and (4), ⃗t = − ln(R) ⃗ t = D − ln(R) the leaves and let R By line 5, L ⃗t (D (L)) = D (− ln(Rt )). ψ⃗ (L) = − ln(Rt ) and ψ t

h

(1)

(t)

Validity. Let γ = (d + 1) ln(1 + 1/m). By Lemma C.3, E Sb | e−L , . . . , e−L

i

= tRt and

e−γ tR

γ b t ≤iS ≤ e tRt . Dividing by t, and since Rt is a function of the leaves, the tower property gives h b | Rt = Rt and e−γ Rt ≤ R b ≤ eγ Rt . In particular, Rt ≤cx R b by Jensen’s inequality. Since e−L E R

b are bounded and nonnegative, Lemma 4.3 shows that L ⃗t is a PLD realization that dominates and R ⃗ ⃗ ⃗ ψ t (L), and that Lt dominates ψt (D (L)). Since the dual of a PLD realization is a PLD realization ⃗ t is also a PLD realization. (Claim A.2), L h i b are bounded and nonnegative, and E [Rt ] = E e−L ≤ 1. Tightness. The variables Rt and R b gives L ⃗t (D (L)). By Claim A.4, ⃗t ⪯(γ,0) ψ⃗ (L) and L ⃗ t ⪯(γ,0) ψ Thus Lemma C.4 with (Y, Yb ) = (Rt , R) 









 t





⃗t (D (L)) for every ε ∈ R. Together with ⃗t ≤ H eε−γ ψ⃗t (L) and H eε L ⃗ t ≤ H eε−γ ψ H eε L validity, these are the two inequalities of Definition 2.1 with the parameters (γ, 0). χ2 error and bias. By the setup, Lemma C.8 applies to this tree, and its exact and computed sums b are (ψ b Thus its two pairs (A, A) ⃗ (L), L ⃗t (D (L)), L ⃗t ) and (ψ ⃗ t ), at the root are Y = tRt and Yb = S. t 2 2 d+1 2 and for each of them with χ (A) < ∞ the χ error is at most ((1 + τ ) − 1)(1 + χ (A)) for := τ 1/(4m(m + 1)). For m ≥ d + 1 we have (d + 1)τ ≤ 1, and thus (1 + τ )d+1 − 1 ≤ e(d+1)τ − 1 ≤ (d + 1)τ + (d + 1)2 τ 2 ≤

d+1 , 4m2

where the last inequality holds since d + 1 ≤ 4(m + 1). This gives the bound on the χ2 error, and Appendix C.5 gives the bound on the bias. Computational complexity. The positive finite part of every distribution in the computation is supported on at most J + d + 1 bands (Claim C.5). Repeated squaring performs at most 2d merges (Claim B.4), and by Theorem 3.1, each of them uses O((J + d)m log(m)) operations and O((J + d)m) space. Besides the input, repeated squaring keeps at most two intermediate distributions, each supported on O((J + d)m) points. Line 5 computes both outputs from fSb in time linear in the size of its support, that is, O((J + d)m). Since d = O(log(t)), the running time is O((J + log(t)) m log(m) log(t)) and the space is O((J + log(t)) m), apart from the computation of L′ in lines 2 and 3, which depends on the representation of the input.

C.3

χ2 Error of a CtD Projection

We first bound the χ2 error of rounding, which Section 4 uses for the comparison with direct convolution. Claim C.6. Let L and U be PLD realizations with χ2 (L) < ∞. 1. If L ⪯ U ⪯(α,0) L for some α ≥ 0, then the χ2 error of U relative to L is between 0 and (eα − 1)(1 + χ2 (L)). 2. If U is obtained by rounding L up to the grid α′ (Z + θ), where α′ > 0 and θ is uniform on [0, 1) ′ and independent of L, then the χ2 error of U relative to L is ((eα − 1)/α′ − 1)(1 + χ2 (L)) ≥ (α′ /2)(1 + χ2 (L)). Proof. For the h first i part, h Fi L (x) h≤ F Ui(x) ≤ F L h (xi− α) for all x. Since the exponential function is L U L+α α increasing, E e ≤ E e ≤ E e = e E eL .

26

For the second part, L is finite since χ2 (L) < ∞. Let Z := L/α′ − θ, so that U − L = α′ (⌈Z⌉ − Z). Conditioned on L, the fractional part of Z is uniform thus ⌈Z⌉ − Z is uniform on h on [0, i 1), and R ′ ′ [0, 1) up to an event of probability zero. Therefore E eU | L = eL 01 eα v dv = eL (eα − 1)/α′ , and ′

(eα − 1)/α′ − 1 =

P

k≥1 α

′k /(k + 1)! ≥ α′ /2.

In contrast, the CtD projection changes the χ2 divergence only to the second order. Lemma C.7 (χ2 error of a CtD projection). Let α′ > 0 and let S ⊆ [0, ∞) be a set of knots that ′ contains 0 and in which every two consecutive positive knots u < v satisfy v ≤ eα u. Let X be a e be its nonnegative random variable whose positive values lie between positive knots of S, and let X CtD projection onto S (as a random map). Then h

i

h

i

h

e 2 ≤ cosh2 (α′ /2) E X 2 E X2 ≤ E X

and

h

i

h

i

i

h

i

e −1 1 E X −1 1{X>0} ≤ E X ≤ cosh2 (α′ /2) E X −1 1{X>0} . e {X>0} e be the PLD realization for which f In particular, let L be a PLD realization, and let L = ΠS (fe−L ). −e L e

 

e = (L, L) e and for (A, A) e = (D (L) , D L e ) with χ2 (A) < ∞, the χ2 error of A e Then for (A, A)

relative to A is between 0 and sinh2 (α′ /2)(1 + χ2 (A)). The factor cosh2 (α′ /2) cannot be improved. ′ If v = eα u are adjacent knots of S and X ≡ (u + v)/2, then the second upper bound holds with equality. Proof. The projection fixes 0 and the knots, and by (P1) it maps positive values to positive values. ′ Thus it suffices to consider a value x of X between adjacent positive knots u < x < v, where v ≤ eα u. The projection maps x to u and to v with the probabilities (v − x)/(v − u) and (x − u)/(v − u), and thus with mean x and variance (x − u)(v − x). Direct calculation gives h

i

e2 | X = x E X

x2

h

=1+

(x − u)(v − x) x2

and

i

e −1 | X = x E X

x−1

=1+

(x − u)(v − x) . uv

Since (x − u)(v − x) ≤ (v − u)2 /4, the second fraction after 1 is at most (v − u)2 /(4uv), and so is the first, since (v − u)2 x2 − 4uv(x − u)(v − x) = ((u + v)x − 2uv)2 ≥ 0. In addition, 1 + (v − u)2 /(4uv) = cosh2 (ln(v/u)/2) ≤ cosh2 (α′ /2). Integrating over the values of X gives the two ′ chains of inequalities. For x = (u + v)/2 and v = eα u, the second fraction is equal to (v − u)2 /(4uv), which gives the equality. 2 For the second statement, let X := e−L   h . iSince χ (A) < ∞, A has no mass at infinity. By (P1) e = 0 = P (X = 0) and E X e = E [X], and thus A e has no mass at infinity either. and (P5), P X h

i

h

i

Therefore the definition of the dual gives E eL = E X −1 for A = L and E eD(L) = E X 2 







e Thus the two chains of inequalities give the claim, since for A = D (L), and similarly for A. 2 ′ 2 ′ cosh (α /2) − 1 = sinh (α /2).

For the grid Gm we can take α′ = ln(1 + 1/m), for which sinh2 (α′ /2) = 1/(4m(m + 1)). Thus for a single projection onto Gm , the χ2 error of Ae relative to A is at most (1 + χ2 (A))/(4m(m + 1)).

27

C.4

Proof of the χ2 Error Bound

Lemma C.8 (χ2 error along the tree). Consider the unfolded tree of Lemma C.3 of depth d with t leaves, which are independent copies of e−L for a PLD realization L such that e−L takes values in {0} ∪ [2j− , 2j+ +1 ] for integers j− ≤ j+ . Let Y and Yb denote the exact and the computed sums at the root, and let τ := 1/(4m(m + 1)). Then h

i

h

i

h

E Y 2 ≤ E Yb 2 ≤ (1 + τ )d+1 E Y 2 and

h

i

h

i

i

h

i

E Y −1 1{Y >0} ≤ E Yb −1 1{Yb >0} ≤ (1 + τ )d+1 E Y −1 1{Y >0} . 



b = (− ln(Y /t), − ln(Yb /t)) and for (A, A) b = (D (− ln(Y /t)) , D − ln(Yb /t) ), In particular, for (A, A)

if χ2 (A) < ∞, then the χ2 error of Ab relative to A is between 0 and ((1 + τ )d+1 − 1)(1 + χ2 (A)). We prove Lemma C.8 by induction on the subtrees of the unfolded tree. All the random variables below are bounded, and on the event that they are positive they are bounded away from zero. Thus all the expectations are finite. Let φ(w) = w2 or φ(w) = 1/w, which are convex on (0, ∞) and satisfy φ(ab) = φ(a)φ(b). For a node n, let Yn and Ybn denote the exact and the computed sums at n, and let L denote the values of the leaves, as in Lemma C.3. Claim C.9 (Moments along the subtrees). For every node n of height h, on the event Yn > 0, h

i

1 ≤ E φ(Ybn /Yn ) | L ≤ (1 + τ )h+1 . h

i

Proof. The lower bound follows from Jensen’s inequality, since E Ybn | L = Yn by Lemma C.3 applied to the subtree of n. By Lemma C.7 applied to a constant y > 0, with α′ = ln(1 + 1/m) by (G1), the projection ỹ of y satisfies E [φ(ỹ/y)] ≤ cosh2 (α′ /2) = 1 + τ . At a leaf, this is the upper bound for h = 0. At an internal node n of height h, conditioned on L, the sums yi at its children ni are fixed, and the computed value at a child with yi = 0 is zero (Lemma C.3). Thus the input s of the projection at n satisfies X yi Ybn s = · i, Yn i: y >0 Yn yi i

which is a convex combination. Since φ is convex and the heights of the children are at most h − 1, the induction hypothesis gives E [φ(s/Yn ) | L] ≤ (1 + τ )h . Conditioned also on s, the projection at h

i

n gives E φ(Ybn /s) | s, L ≤ 1 + τ . Since Ybn /Yn = (Ybn /s)(s/Yn ) and φ is multiplicative, the tower property gives the upper bound. Proof of Lemma C.8. The root has height at most d. Multiplying Claim C.9 at the root by Y 2 , or by Y −1 on the event Y > 0, and taking expectations gives the two chains of inequalities, where the reciprocal moments are restricted to the event Y > 0, which is the event Yb > 0. For the second statement, by the definition of the dual, the masses of − ln(Y /t) and of its dual at infinity are P (Y = 0) and 1 − E [Y ] /t, and the expectations of e− ln(Y /t) and of the exponential of its dual h

i

over their finite parts are t E Y −1 1{Y >0} and E Y 2 /t2 , and similarly for Yb . By Lemma C.3, 







h i

P Yb = 0 = P (Y = 0) and E Yb = E [Y ], and thus A and Ab have the same mass at infinity. If χ2 (A) < ∞, this mass is zero, and the two chains of inequalities give the bounds on the χ2 errors. Lemma C.8 gives the bound on the χ2 error in the third part of Theorem 4.1, since − ln(Y /t) = ⃗t (D (L)) by (4). ψ⃗t (L) and D (− ln(Y /t)) = ψ 28

C.5

Proof of the Bias Bound

Proof of the bound on the bias in Theorem 4.1. We use the unfolded tree of the proof of Theorem 4.1, whose depth is d := ⌈log2 (t)⌉, and the notation Yn , Ybn and L of Appendix C.4, with Y and Yb the sums at the root. Let α′ := ln(1 + 1/m). By line 5 and (4), the allocation PLD and the output are A = − ln(Y /t) and Ab = − ln(Yb /t) in the add direction, and their duals in the remove direction. As in the proof of Lemma C.8, A and Ab have the same mass at infinity. We first bound the bias of a single projection. For the projection ỹ of a fixed value y > 0, let Z := ỹ/y. By (P1), (P5) and (G1), − ln(ỹ) is the CtD discretization of the loss − ln(y) on a grid with spacing at most α′ [FS26b, Definition 6.1], which changes this loss by − ln(Z). In the dual, the loss ln(y) changes by ln(Z), with the distribution of Z weighted by Z. Thus the bound on the bias of a CtD discretization [FS26b, Lemma 6.2 and its proof], applied in both directions, gives 0 ≤ E [− ln(Z)] ≤

α′2 8

0 ≤ E [Z ln(Z)] ≤

and

α′2 . 8

Let φ(w) = − ln(w) or φ(w) = w ln(w), which are convex on (0, ∞) and satisfy φ(1) = 0. We show that for every node n of height h, on the event Yn > 0, h

i

0 ≤ E φ(Ybn /Yn ) | L ≤

(h + 1)α′2 . 8 h

(9)

i

The lower bound follows from Jensen’s inequality, since E Ybn | L = Yn by Lemma C.3 applied to the subtree of n. We prove the upper bound by induction on h, as in the proof of Claim C.9. At a leaf, it is the bound for a single projection. At an internal node n, conditioned on L, the ratio V := s/Yn for the input s of the projection at n is a convex combination of the ratios Ybni /Yni at the children ni of n with Yni > 0. Since φ is convex and the heights of the children are at most h − 1, the induction hypothesis gives E [φ(V ) | L] ≤ hα′2 /8. Conditioned also on V , the ratio Ybn /s is the ratio Z for y = s, and Ybn /Yn = V Z. Therefore E [− ln(V Z) | V, L] ≤ − ln(V ) +

α′2 , 8

E [V Z ln(V Z) | V, L] ≤ V ln(V ) +

V α′2 . 8

Since E [V | L] = 1, taking expectations conditioned on L gives the upper bound. The root has height at most d. In the add direction, both losses are finite exactly when Y > 0, h i b b b and then A − A = − ln(Y /Y ). Thus the bias is the expectation of E φ(Y /Y ) | L 1{Y >0} for φ(w) = − ln(w). In the remove direction, the finite parts of A and Ab are as ln(Y /t) and i h distributed b b b ln(Y /t), weighted by Y /t and Y /t, respectively (Section 2). Since E Y | L = Y , the bias is then h

i

the expectation of (Y /t) E φ(Yb /Y ) | L 1{Y >0} for φ(w) = w ln(w). Therefore (9) and E [Y /t] ≤ 1 give the bounds on the bias.

C.6

Guarantees for a Given Accuracy

⃗ and the remove Proof of Corollary 4.5. Thealgorithm returns add output of Algorithm 1 on L   the   ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ output of Algorithm 1 on D L , since ψt (D D L ) = ψt (L). (If L and L are duals of each other, then a single call gives both outputs.) Consider one of these calls, and denote its input by L and ⃗t and L ⃗ t . The bands chosen in line 2 satisfy j+ − j− ≤ (ℓ+ − ℓ− )/ ln(2) + 1 ≤ its outputs by L   ⃗ are those of L. ⃗ Let α′ = ln(1 + 1/m), so that ∆/ ln(2) + 1, since the finite values of −D L 29

⃗t and L ⃗t (⌈log2 (t)⌉ + 1)α′ ≤ (⌈log2 (t)⌉ + 1)/m ≤ α. Thus by the first two parts of Theorem 4.1, L ⃗ ⃗ (α, 0)-tightly dominate ψ t (L) and ψt (D (L)). χ2 error and bias. Since m ≥ (⌈log2 (t)⌉ + 1)/α ≥ ⌈log2 (t)⌉ + 1, the third part of Theorem 4.1 bounds the χ2 error of each output relative to the corresponding allocation PLD Ψ of L with χ2 (Ψ) < ∞ by (⌈log2 (t)⌉ + 1)(1 + χ2 (Ψ))/(4m2 ) ≤ α2 (1 + χ2 (Ψ))/(4(⌈log2 (t)⌉ + 1)), and the bias of each output by (⌈log2 (t)⌉ + 1)α′2 /8 ≤ α2 /(8(⌈log2 (t)⌉ + 1)). −L′ lie in the J + 1 = O(1 + ∆) Complexity. Since j+ − j− ≤ ∆/ ln(2)  + 1, the positive values of e ⃗ and L′ takes O(S + Jm) operations, since every atom bands j− , . . . , j+ + 1. Computing D L contributes to at most two knots. Substituting m = O(log(t)/α) into Theorem 4.1 gives the bounds on the running time and the space.

References [AFKRT26]

Hilal Asi, Vitaly Feldman, Hannah Keller, Guy Rothblum, and Kunal Talwar. “PREAMBLE: Private and Efficient Aggregation via Block Sparse Vectors”. In: NeurIPS. 2026.

[BBG18]

Borja Balle, Gilles Barthe, and Marco Gaboardi. “Privacy amplification by subsampling: Tight analyses via couplings and divergences”. In: NeurIPS. 2018.

[BKOZ12]

Gilles Barthe, Boris Köpf, Federico Olmedo, and Santiago Zanella Beguelin. “Probabilistic relational reasoning for differential privacy”. In: ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages (POPL). 2012.

[CGHSG25]

Christopher Choquette-Choo, Arun Ganesh, Saminul Haque, Thomas Steinke, and Abhradeep Guha Thakurta. “Near-Exact Privacy Amplification for Matrix Mechanisms”. In: International Conference on Learning Representations (ICLR). 2025.

[Chu+24]

Lynn Chua, Badih Ghazi, Charlie Harrison, Pritish Kamath, Ravi Kumar, Ethan Jacob Leeman, Pasin Manurangsi, Amer Sinha, and Chiyuan Zhang. “Balls-and-Bins Sampling for DP-SGD”. In: International Conference on Artificial Intelligence and Statistics (AISTATS). 2024.

[DCO25]

Andy Dong, Wei-Ning Chen, and Ayfer Ozgur. “Leveraging Randomness in Model and Data Partitioning for Privacy Amplification”. In: ICML. 2025.

[DG26]

Andy Dong and Arun Ganesh. “Privacy Amplification for BandMF via b-Min-Sep Subsampling”. 2026. arXiv: 2602.09338.

[DGKKM22]

Vadym Doroshenko, Badih Ghazi, Pritish Kamath, Ravi Kumar, and Pasin Manurangsi. “Connect the Dots: Tighter Discrete Approximations of Privacy Loss Distributions”. In: Proceedings on Privacy Enhancing Technologies (2022).

[DKMMN06] Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor. “Our data, ourselves: Privacy via distributed noise generation”. In: Advances in Cryptology (EUROCRYPT). 2006. [DR16]

Cynthia Dwork and Guy N Rothblum. “Concentrated differential privacy”. 2016. arXiv: 1603.01887.

[FS25]

Vitaly Feldman and Moshe Shenfeld. “Privacy amplification by random allocation”. In: NeurIPS. 2025.

[FS26]

Vitaly Feldman and Moshe Shenfeld. “Efficient privacy loss accounting for subsampling and random allocation”. 2026. arXiv: 2602.17284. 30

[FS26b]

Vitaly Feldman and Moshe Shenfeld. “Gaussian Approximation for Composition and Subsampling of Differentially Private Algorithms”. 2026. arXiv: 2610.33333.

[GKKM22]

Badih Ghazi, Pritish Kamath, Ravi Kumar, and Pasin Manurangsi. “Faster privacy accounting via evolving discretization”. In: ICML. 2022.

[GLW21]

Sivakanth Gopi, Yin Tat Lee, and Lukas Wutschitz. “Numerical composition of differential privacy”. In: NeurIPS. 2021.

[Goo22]

Google. dp-accounting. https://github.com/google/differential- privacy/ tree/main/python/dp_accounting. 2022.

[Hac07]

Wolfgang Hackbusch. “Fast and exact projected convolution for non-equidistant grids”. In: Computing (2007).

[KH21]

Antti Koskela and Antti Honkela. “Computing differential privacy guarantees for heterogeneous compositions using FFT”. 2021. arXiv: 2102.12412.

[KJH20]

Antti Koskela, Joonas Jälkö, and Antti Honkela. “Computing tight differential privacy guarantees using FFT”. In: International Conference on Artificial Intelligence and Statistics (AISTATS). 2020.

[KJPH21]

Antti Koskela, Joonas Jälkö, Lukas Prediger, and Antti Honkela. “Tight differential privacy for discrete-valued mechanisms and for the subsampled gaussian mechanism using FFT”. In: International Conference on Artificial Intelligence and Statistics (AISTATS). 2021.

[Met21]

Meta. Opacus: Differential Privacy Library for PyTorch. https://github.com/ pytorch/opacus. 2021.

[Mic21]

Microsoft. PRV Accountant. https://github.com/microsoft/prv_accountant. 2021.

[MM18]

Sebastian Meiser and Esfandiar Mohammadi. “Tight on budget? tight bounds for r-fold approximate differential privacy”. In: ACM SIGSAC Conference on Computer and Communications Security (CCS). 2018.

[PLD]

PLD accounting for subsampling and random allocation. Code: https://github. com/moshenfeld/PLD_accounting.git, branch MS-FFT-paper (the version used in this paper), package: https://pypi.org/project/PLD-accounting/. 2026.

[PW12a]

Gilles Pagès and Benedikt Wilbertz. “Dual Quantization for random walks with application to credit derivatives”. In: The Journal of Computational Finance (2012).

[PW12b]

Gilles Pagès and Benedikt Wilbertz. “Intrinsic stationarity for vector quantization: Foundation of dual quantization”. In: SIAM Journal on Numerical Analysis (2012).

[SK26]

Jan Schuchardt and Nikita Kalinin. “Sampling-Free Privacy Accounting for Matrix Mechanisms under Random Allocation”. 2026. arXiv: 2601.21636.

[SMM19]

David M Sommer, Sebastian Meiser, and Esfandiar Mohammadi. “Privacy Loss Classes: The Central Limit Theorem in Differential Privacy”. In: Proceedings on Privacy Enhancing Technologies (2019).

[SS07]

Moshe Shaked and J George Shanthikumar. Stochastic orders. Springer, 2007.

[ZDW22]

Yuqing Zhu, Jinshuo Dong, and Yu-Xiang Wang. “Optimal accounting of differential privacy via characteristic function”. In: International Conference on Artificial Intelligence and Statistics (AISTATS). 2022.

31

Related documents

• Limiting Access to Certain Anonymous Information: From the Group Right to Privacy to the Principle of Protecting the Vulnerable#87484
• A privacy policies dataset in Greek in the GDPR era#964851
• Privacy Without the Right to Privacy#116781
• Oblivious Pseudorandom Functions (OPRFs) Using Prime-Order Groups#132392
• Near-Optimal Generalized Private Testing#216744