figurec

# Numerical Composition of Differential Privacy111Author ordering is
alphabetical. Code is available at
https://github.com/microsoft/prv_accountant.

Sivakanth Gopi  Microsoft
<EMAIL_ADDRESS>Yin Tat Lee  University of Washington
<EMAIL_ADDRESS>Lukas Wutschitz  Microsoft
<EMAIL_ADDRESS>
###### Abstract

We give a fast algorithm to optimally compose privacy guarantees of
differentially private (DP) algorithms to arbitrary accuracy. Our method is
based on the notion of _privacy loss random variables_ to quantify the privacy
loss of DP algorithms. The running time and memory needed for our algorithm to
approximate the privacy curve of a DP algorithm composed with itself $k$ times
is $\tilde{O}(\sqrt{k})$. This improves over the best prior method by Koskela
et al. [KH21] which requires $\tilde{\Omega}(k^{1.5})$ running time. We
demonstrate the utility of our algorithm by accurately computing the privacy
loss of DP-SGD algorithm of Abadi et al. [ACG+16] and showing that our
algorithm speeds up the privacy computations by a few orders of magnitude
compared to prior work, while maintaining similar accuracy.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Contributions
  2. 2 DP Preliminaries
  3. 3 Privacy Loss Random Variables (PRVs)
  4. 4 Numerical composition of privacy curves
  5. 5 Error analysis
  6. 6 Experiments
    1. 6.1 Comparison with [KJPH21]
  7. A Effect of floating point arithmetic
  8. B Privacy Loss Random Variables
    1. B.1 Examples of privacy loss random variables
    2. B.2 Subsampling
  9. C Missing Proofs in Error Analysis
    1. C.1 Facts about Coupling Approximation
    2. C.2 Bounding the error using tail bounds of PRVs
    3. C.3 Tail Bound for PRVs
    4. C.4 Proof of Theorem 5.5

##  1 Introduction

Differential privacy (DP) introduced by [DMNS06] provides a provable and
quantifiable guarantee of privacy when the results of an algorithm run on
private data are made public. Formally, we can define an
$(\varepsilon,\delta)$-differentially private algorithm as follows.

######  Definition 1.1 ($(\varepsilon,\delta)$-DP [DMNS06, DKM+06]).

An algorithm $\mathcal{M}$ is $(\varepsilon,\delta)$-DP if for any two
neighboring databases $D,D^{\prime}$ differing in exactly one user and any
subset $S$ of outputs, we have $\Pr[\mathcal{M}(D)\in S]\leq
e^{\varepsilon}\Pr[\mathcal{M}(D^{\prime})\in S]+\delta.$

Intuitively, it says that looking at the outcome of $\mathcal{M}$, we cannot
tell whether it was run on $D$ or $D^{\prime}$. Hence an adversary cannot
infer the existence of any particular user in the input database, and
therefore cannot learn any personal data of any particular user.

DP algorithms have an important property called _composition_. Suppose $M_{1}$
and $M_{2}$ are DP algorithms and say $M(D)=(M_{1}(D),M_{2}(D))$, i.e., $M$
runs both the algorithms on $D$ and outputs their results. Then $M$ is also a
DP algorithm.

######  Proposition 1.2 (Simple composition [DKM+06, DL09]).

If $M_{1}$ is $(\varepsilon_{1},\delta_{1})$-DP and $M_{2}$ is
$(\varepsilon_{2},\delta_{2})$-DP, then $M(D)=(M_{1}(D),M_{2}(D))$ is
$(\varepsilon_{1}+\varepsilon_{2},\delta_{1}+\delta_{2})$-DP.

This also holds under _adaptive composition_ (denoted by $M=M_{2}\circ
M_{1}$), where $M_{2}$ can look at both the database and the output of
$M_{1}$.222Here $M(D)=(M_{1}(D),M_{2}(D,M_{1}(D)))$. It turns out that both
compositions enjoy much better DP guarantees than this simple composition
rule. Let $M$ be an $(\varepsilon,\delta)$-DP algorithm and let $M^{\circ k}$
denote the (adaptive) composition of $M$ with itself $k$ times. The naive
composition rule shows that $M^{\circ k}$ is $(k\varepsilon,k\delta)$-DP. This
was significantly improved in [DRV10].

######  Proposition 1.3 (Advanced composition [DRV10, DR+14]).

If $M$ is $(\varepsilon,\delta)$-DP, then $M^{\circ k}$ is
$(\varepsilon^{\prime},k\delta+\delta^{\prime})$-DP where

$\varepsilon^{\prime}=\varepsilon\sqrt{2k\log\left(\frac{1}{\delta^{\prime}}\right)}+k\varepsilon(e^{\varepsilon}-1).$

Note that if $\varepsilon=O\left(\frac{1}{\sqrt{k}}\right)$ and
$\delta=o\left(\frac{1}{k}\right)$, then $M^{\circ k}$ satisfies
$(O_{\delta^{\prime}}(1),\delta^{\prime})$-DP. Using simple composition
(Proposition 1.2), we can only claim that $M^{\circ k}$ is
$(O(\sqrt{k}),o(1))$-DP. Thus advanced composition often results in
$\sqrt{k}$-factor savings in privacy which is significant in practice. The
optimal DP guarantees for $k$-fold composition of an $(\varepsilon,\delta)$-DP
algorithm were finally obtained by [KOV15]. For composing different
algorithms, the situation is more complicated. If $M_{1},M_{2},\dots,M_{k}$
are DP algorithms such that $M_{i}$ is $(\varepsilon_{i},\delta_{i})$-DP, then
it is shown by [MV16] that computing the _exact_ DP guarantees for
$M=M_{1}\circ M_{2}\circ\dots\circ M_{k}$ is #P-complete. They also give an
algorithm to approximate the DP guarantees of $M$ to desired accuracy $\eta$
which runs in

$\tilde{O}\left(\frac{k^{3}\bar{\varepsilon}(1+\bar{\varepsilon})}{\eta}\right)$
(1)

time where
$\bar{\varepsilon}=(\sum_{i=1}^{k}\varepsilon_{i})/k.$333$\varepsilon$ has an
additive error of $\eta$ and $\delta$ has a multiplicative error of $\eta$. If
each $\varepsilon_{i}\approx\frac{1}{\sqrt{k}}$ (so that $M$ will satisfy
reasonable privacy guarantees by advanced composition), then the running time
is $\tilde{O}(k^{2.5}/\eta).$

In most situations, DP algorithms come with a collection of
$(\varepsilon,\delta)$-DP guarantees, i.e., for each value of $\varepsilon$,
there exists $\delta$ such that the algorithm is $(\varepsilon,\delta)$-DP.

######  Definition 1.4 (Privacy curve).

A DP algorithm $M$ is said to have privacy curve $\delta:\mathbb{R}\to[0,1]$,
if for every $\varepsilon\in\mathbb{R}$, $M$ is
$(\varepsilon,\delta(\varepsilon))$-DP.

For example the privacy curve of a Gaussian mechanism (with sensitivity $1$
and noise scale $\sigma$) is given by
$\delta(\varepsilon)=\Phi\left(-\varepsilon\sigma+1/2\sigma\right)-e^{\varepsilon}\Phi\left(-\varepsilon\sigma-1/2\sigma\right)$
where $\Phi(\cdot)$ is the Gaussian CDF [BW18]. Suppose we want to compose
several Gaussian mechanisms, which $(\varepsilon,\delta)$-DP guarantee should
we choose for each mechanism? Any choice will lead to suboptimal DP guarantees
for the final composition. Instead, we need a way to compose the privacy
curves directly. This was suggested through the use of privacy region in
[KOV15] and explicitly studied in the $f$-DP framework of [DRS19]. $f$-DP is a
dual way (and equivalent) to look at the privacy curve $\delta(\varepsilon).$

Independently, an algorithm called _Privacy Buckets_ for approximately
composing privacy curves using the notion of was initiated in [MM18]. This
algorithm depends on the notion of _Privacy Loss Random Variable_ (PRV)
[DR16], whose distribution is called _Privacy Loss Distribution_ (PLD). For
any DP-algorithm, one can associate a PRV and the privacy curve of that
algorithm can be easily obtained from the PRV. The really useful property of
PRVs is that under adaptive composition, they just add up; the PRV $Y$of the
composition $M=M_{1}\circ M_{2}\circ\dots\circ M_{k}$ is given by
$Y=\sum_{i=1}^{k}Y_{i}$ where $Y_{i}$ is the PRV of $M_{i}.$444[KJH+20] only
state this for non-adaptive composition. In this paper we show how to extend
this to adaptive composition as well. Therefore, one can find the distribution
of $Y$ by the convolution of the distributions of $Y_{1},Y_{2},\dots,Y_{k}$.
In an important paper, [KJH+20] proposed that one can speed up the
convolutions using Fast Fourier Transform (FFT). Explicit error bounds were
obtained for the approximation obtained by their algorithm in [KJH+20, KJPH21,
KH21]. The running time of this algorithm was analyzed in [KH21] where it was
shown that the privacy curve $\delta_{M}(\varepsilon)$ of $M=M_{1}\circ
M_{2}\circ\dots\circ M_{k}$ can be computed up to an additive error of
$\delta_{\mathrm{error}}$ in time

$\tilde{O}\left(\frac{k^{3}\bar{\varepsilon}}{\delta_{\mathrm{error}}}\right),$
(2)

if each algorithm $M_{i}$ is satisfies $(\varepsilon_{i},0)$-DP and
$\bar{\varepsilon}=\frac{1}{k}\sum_{i=1}^{k}\varepsilon_{i}$. Assuming that
each $\varepsilon_{i}\approx\frac{1}{\sqrt{k}}$, we get
$\tilde{O}(k^{2.5}/\delta_{\mathrm{error}})$ running time. Note that this is
slightly worse than (1), where the denominator $\eta$ is the multiplicative
error in $\delta_{M}$. When composing the same algorithm with itself for $k$
times, the running time can be improved to
$\tilde{O}\left(\frac{k^{2}\bar{\varepsilon}}{\delta_{\mathrm{error}}}\right)$,
which is $\tilde{O}\left(\frac{k^{1.5}}{\delta_{\mathrm{error}}}\right)$ when
$\bar{\varepsilon}=\frac{1}{\sqrt{k}}$.

#### Moments Accountant and Renyi DP

In an influential paper where they introduce Differentially Private Deep
Learning, [ACG+16] proposed a method called the Moments Accountant (MA) for
giving an upper bound the privacy curve of a composition of DP algorithms.
They applied their method to bound the privacy loss of differentially private
Stochastic Gradient Descent (DP-SGD) algorithm which they introduced.
Analyzing the privacy loss of DP-SGD involves composing the privacy curve of
each iteration of training with itself $k$ times, where $k$ is the total of
number of training iterations. Typical values of $k$ range from $1000$ to
$300000$ (such as when training large models like GPT3). The Moments
Accountant was subsumed into the framework of Renyi Differential Privacy (RDP)
introduced by [Mir17]. The running time of these accountants are independent
of $k$, but they only give an upper bound and cannot approximate the privacy
curve to arbitrary accuracy.

DP-SGD is one of the most important DP algorithms in practice, because one can
use it to train neural networks to achieve good privacy-vs-utility tradeoffs.
Therefore obtaining accurate and tight privacy guarantees for DP-SGD is
important. For example reducing $\varepsilon$ from $2$ to $1$, can mean that
one can train the network for 4 times more epochs while staying within the
same privacy budget. Therefore DP-SGD is one of the main motivations for this
work.

There are also situations when the PRVs do not have bounded moments and so
Moments Accountant or Renyi DP cannot be applied for analyzing privacy. An
example of such an algorithm is the DP-SGD-JL algorithm of [BGK+21] which uses
numerical composition of PRVs to analyze privacy.

#### GDP Accountant

[DRS19, BDLS19] introduced the notion of Gaussian Differential Privacy (GDP)
and used it to develop an accountant for DP-SGD. The accountant is based on
central limit theorem and only gives an approximation to the true privacy
curve, where the approximation gets better with $k$. But as we show in Figure
1, GDP accountant can significantly underreport the true epsilon value.

Several different notions of privacy were introduced for obtaining good upper
bounds on the privacy curve of composition of DP algorithms such as
Concentrated DP (CDP) [DR16, BS16], Truncated CDP [BDRS18] etc. None of these
methods can approximate the privacy curve of compositions to arbitrary
accuracy. The notion of $f$-DP introduced by [DRS19], allows for a lossless
composition theorem, but computing the privacy curve of composition seems
computationally hard and they do not give any algorithms for doing it.

###  1.1 Our Contributions

The main contribution of this work is a new algorithm with an improved
analysis for computing the privacy curve of the composition of a large number
of DP algorithms.

######  Theorem 1.5 (Informal version of Theorem 5.5).

Suppose $M_{1},M_{2},\dots,M_{k}$ are DP algorithms. Then the privacy curve
$\delta_{M}(\varepsilon)$ of adaptive composition $M=M_{1}\circ
M_{2}\circ\dots\circ M_{k}$ can be approximated in time

$O\left(\frac{\varepsilon_{\mathrm{upper}}\
k^{1.5}\log{k}\sqrt{\log\frac{1}{\delta_{\mathrm{error}}}}}{\varepsilon_{\mathrm{error}}}\right),$
(3)

where $\varepsilon_{\mathrm{error}}$ is the additive error in $\varepsilon$,
$\delta_{\mathrm{error}}$ is the additive error in $\delta$ and
$\varepsilon_{\mathrm{upper}}$ is an upper bound on
$\max\left\\{\varepsilon_{M}(\delta_{\mathrm{error}}),\max_{i}\varepsilon_{M_{i}}\left(\frac{\delta_{\mathrm{error}}}{k}\right)\right\\}.$555$\varepsilon_{M}(\delta)$
is the inverse of $\delta_{M}(\varepsilon)$.

If each $M_{i}$ satisfies $\left(\frac{1}{\sqrt{k}},\frac{o(1)}{k}\right)$-DP,
then by advanced composition (Proposition 1.3), we can set
$\varepsilon_{\mathrm{upper}}=O(1)$. Therefore the running time of our
algorithm in this case is
$\tilde{O}\left(\frac{k^{1.5}\sqrt{\log\frac{1}{\delta_{\mathrm{error}}}}}{\varepsilon_{\mathrm{error}}}\right).$
We can save a factor of $k$, when we compose the same algorithm with itself
$k$ times.

######  Theorem 1.6.

Suppose $M$ is a DP algorithm. Then the privacy curve $\delta_{M^{\circ
k}}(\varepsilon)$ of $M$ (adaptively) composed with itself $k$ times can be
approximated in time

$O\left(\frac{\varepsilon_{\mathrm{upper}}\
k^{\frac{1}{2}}\log{k}\sqrt{\log\frac{1}{\delta_{\mathrm{error}}}}}{\varepsilon_{\mathrm{error}}}\right),$
(4)

where $\varepsilon_{\mathrm{error}}$ is the additive error in $\varepsilon$,
$\delta_{\mathrm{error}}$ is the additive error in $\delta$ and
$\varepsilon_{\mathrm{upper}}$ is an upper bound on
$\max\left\\{\varepsilon_{M^{\circ
k}}(\delta_{\mathrm{error}}),\varepsilon_{M}\left(\frac{\delta_{\mathrm{error}}}{k}\right)\right\\}.$

Thus we improve the state-of-the-art by at least a factor of $k$ in running
time. We also note that our algorithm improves the memory required by a factor
of $k.$ See Figure 1 for a comparison of our algorithm with that of [KJPH21].
Also note that RDP Accountant (equivalent to the Moments Accountant)
significantly overestimates the true $\varepsilon$, while the GDP Accountant
significantly underestimates the true $\varepsilon.$ In contrast, the upper
and lower bounds provided by our algorithm lie very close to each other.

$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$$10$$15$$20$$25$$30$$35$$40$$45$$50$

$\varepsilon$

DPSGD
steps$\varepsilon_{\textnormal{up}}$Ours$\varepsilon_{\textnormal{low}}$Ours$\varepsilon_{\textnormal{up}}$[KJPH21]$\varepsilon_{\textnormal{low}}$[KJPH21]

(a) Our algorithm gives much closer upper and lower bounds on the true privacy
curve compared to [KJPH21], under the same mesh size of $4\times 10^{-5}$. Our
upper and lower bounds are nearly coinciding.

$0$$0.5$$1$$1.5$$2$$2.5$$0$$500$$1000$$1500$$2000$

$\varepsilon$

DPSGD stepsOurs UpperOurs LowerGDPRDP

(b) Our algorithm can improve significantly over the RDP Accountant. We also
see that GDP Accountant can significantly underreport the true $\varepsilon$.
We have set $\varepsilon_{\mathrm{error}}=0.1$,
$\delta_{\mathrm{error}}=\delta/1000$ here.

Figure 1: Case study on DP-SGD. Sampling probability $p=10^{-3}$, noise scale
$\sigma=0.8$, $\delta=10^{-7}$.

#### Our Techniques

Our algorithm (also the prior work of [KJH+20]) proceeds by approximating the
privacy loss random variables (PRVs) by truncating and discretizing them. We
then use Fast Fourier Transform (FFT) to convolve the distributions
efficiently. The main difference is in the approximation procedure and the
error analysis. In the approximation procedure, we correct the approximation
so that the expected value of the discretization matches with the expected
value of the PRV.

To analyze the approximation error, we introduce the concept of _coupling
approximation_ (Definition 5.1), which is a variant of Wasserstein (optimal
transport) distance specifically tailored to this application. We first show
that the approximation output by our algorithm to each privacy random variable
is a good coupling approximation. We then show that when independent coupling
approximations are added, cancellation happens between the errors due to
Hoeffding bound, producing a much better coupling approximation than one
naively expects from the triangle inequality. This allows us to choose the
mesh size in our discretization to be $\approx\frac{1}{\sqrt{k}}$, whereas
[KH21] choose a mesh size of $\approx\frac{1}{k}$. The other improvement is
the truncation procedure. We give a tight tail bound of the PRVs (Lemma 5.4).
This allows us to choose the domain size for in truncation to be
$\approx\tilde{O}(1)$, whereas [KH21] choose $\approx\tilde{O}(\sqrt{k})$.
Both ideas together saves a factor of $k$ in the run time and memory.

For the analysis, the previous paper analyzes the discretization error by
studying the stability of convolution. This leads to complicated calculations
with the runtime linear in $1/{\delta_{\mathrm{error}}}$ (see (2)). Since
$\delta_{\mathrm{error}}\ll\delta\ll 1/N$ is required to give meaningful
privacy guarantee ($N$ is the number of users), this term
$1/{\delta_{\mathrm{error}}}$ is huge. In this paper, we show various facts
about how coupling approximation accumulates and use them to give a runtime
depending only on $\sqrt{\log(1/{\delta_{\mathrm{error}}})}$.

##  2 DP Preliminaries

Given a DP algorithm $\mathcal{M}$, for each value of $\varepsilon\geq 0$,
there exists some $\delta\in[0,1]$ such that $\mathcal{M}$ is
$(\varepsilon,\delta)$-DP. We can represent all these privacy guarantees by a
function $\delta_{\mathcal{M}}(\varepsilon):\mathbb{R}^{\geq 0}\to[0,1]$ and
say that $\delta_{\mathcal{M}}(\cdot)$ is the _privacy curve_ of
$\mathcal{M}.$ This inspires the following definition of a privacy curve
between two random variables.

######  Definition 2.1 (Privacy curve).

Given two random variables $X,Y$ supported on some set $\Omega$, define
$\delta(X||Y):\mathbb{R}\to[0,1]$ as:

$\delta(X||Y)(\varepsilon)=\sup_{S\subset\Omega}\Pr[Y\in
S]-e^{\varepsilon}\Pr[X\in S].$

Therefore an algorithm $\mathcal{M}$ is $(\varepsilon,\delta)$-DP iff
$\delta\left(\mathcal{M}(D)||\mathcal{M}(D^{\prime})\right)(\varepsilon)\leq\delta$
for all neighboring databases $D,D^{\prime}.$

######  Remark 2.2.

Note that not all functions $\delta:\mathbb{R}\to[0,1]$ are privacy curves. A
characterization of privacy curves can be obtained using the $f$-DP framework
of [DRS19]. The notion of privacy curve $\delta(X||Y)$ and tradeoff function
$T(X||Y)$ are dual to each other via convex duality [DRS19]. This implies a
characterization of privacy curves as shown in [ZDW21].

######  Definition 2.3 (Composition of privacy curves [DRS19]).

Let $\delta_{1}\equiv\delta(X_{1}||Y_{1})$ and
$\delta_{2}\equiv\delta(X_{2}||Y_{2})$ be any two privacy curves. The
composition of the privacy curves, denoted by $\delta_{1}\otimes\delta_{2}$,
is defined as

$\delta_{1}\otimes\delta_{2}\equiv\delta\left((X_{1},X_{2})||(Y_{1},Y_{2})\right)$

where $X_{1},X_{2}$ are independently sampled and $Y_{1},Y_{2}$ are
independently sampled.

Note that there can be many pairs of random variables which have the same
privacy curve, but the above operation is well-defined. If
$\delta(X_{1}||Y_{1})\equiv\delta(X_{1}^{\prime}||Y_{1}^{\prime})$ and
$\delta(X_{2}||Y_{2})\equiv\delta(X_{2}^{\prime}||Y_{2}^{\prime})$, then it
was shown by [DRS19] that

$\delta\left((X_{1},X_{2})||(Y_{1},Y_{2})\right)=\delta\left((X_{1}^{\prime},X_{2}^{\prime})||(Y_{1}^{\prime},Y_{2}^{\prime})\right).$

[DRS19] also show that $\otimes$ is a commutative and associative operation.

Given two DP algorithms $M_{1}$ and $M_{2}$, the adaptive composition
$(M_{2}\circ M_{1})(D)$ is an algorithm which outputs
$(M_{1}(D),M_{2}(D,M_{1}(D))$, i.e., $M_{2}$ can look at the database $D$ and
also the output of the previous algorithm $M_{1}(D)$. Adaptive composition of
more than two algorithms is similarly defined. Suppose $M_{1}$ has privacy
curve $\delta_{1}$ and $M_{2}$ has privacy curve $\delta_{2}$ (i.e.,
$M_{2}(\cdot,y)$ is a DP algorithm with privacy curve $\delta_{2}$ for any
fixed $y.$). The following composition theorem shows how to get the privacy
curve of $M_{2}\circ M_{1}.$

######  Theorem 2.4 (Composition theorem [DRS19]).

Let $M_{1},M_{2},\dots,M_{k}$ be DP algorithms with privacy curves given by
$\delta_{1},\delta_{2},\dots,\delta_{k}$ respectively. The privacy curve of
the adaptive composition $M_{k}\circ M_{k-1}\circ\dots\circ M_{1}$ is given by
$\delta_{1}\otimes\delta_{2}\otimes\cdots\otimes\delta_{k}.$

##  3 Privacy Loss Random Variables (PRVs)

The notion of _privacy loss random variables_ (PRVs) is a unique way to assign
a pair $(X,Y)$ for any privacy curve $\delta$ such that
$\delta\equiv\delta(X||Y).$ PRVs allow us to compute composition of two
algorithms via summing random variables (Theorem 3.5) (equivalently,
convolving their distributions). Thus PRVs can be thought of as a
_reparametrization of privacy curves_ where composition becomes convolution.
In this paper, we differ from the usual definition of PRVs given in [DR16,
KJH+20], which are tied to a specific algorithm. Instead we think of them as a
reparametrization of privacy curves and study them directly. This allows us to
succinctly prove many useful properties of PRVs.

Let $\overline{\mathbb{R}}=\mathbb{R}\cup\\{-\infty,\infty\\}$ be the extended
real line where we define $\infty+x=\infty$ and $-\infty+x=-\infty$ for
$x\in\mathbb{R}.$

######  Definition 3.1 (Privacy loss random variables (PRVs)).

Given a privacy curve $\delta:\mathbb{R}\to[0,1]$, we say that $(X,Y)$ are
_privacy loss random variables_ for $\delta$, if they satisfy the following
conditions:

  * •

$X,Y$ are supported on $\overline{\mathbb{R}}$,

  * •

$\delta(X||Y)\equiv\delta$,

  * •

$Y(t)=e^{t}X(t)$ for every $t\in\mathbb{R}$ and

  * •

$Y(-\infty)=0$ and $X(\infty)=0$

where $X(t),Y(t)$ are probability density functions of $X,Y$ respectively.

Mathematically, the correct way to write the condition $Y(t)=e^{t}X(t)$ is to
say that
$\operatorname{\mathbb{E}}_{Y}[\phi(Y)]=\operatorname{\mathbb{E}}_{X}[\phi(X)e^{X}]$
for all test functions $\phi:\overline{\mathbb{R}}\to[0,1]$ with
$\phi(\infty)=\phi(-\infty)=0.$ This will generalize to all situations where
$X,Y$ are continuous or discrete or both. For ease of exposition, we ignore
this complication and assume that $X(t),Y(t)$ represent the PDFs if $X,Y$ are
continuous at $t$, or the probability masses if they have point masses at $t$.

The following theorem shows that the PRVs for a privacy curve
$\delta=\delta(P||Q)$ are given by the log-likelihood random variables of
$P,Q.$

######  Theorem 3.2.

Let $\delta:\mathbb{R}\to[0,1]$ be a privacy curve given by
$\delta\equiv\delta(P||Q)$ where $P,Q$ are two random variables supported on
$\Omega$. The PRVs $(X,Y)$ for the privacy curve $\delta$ are given by666Here
$Q(\omega)$ and $P(\omega)$ are the probability density functions of $Q,P$
respectively. Note that the mathematically precise way is to replace the ratio
$\frac{Q(\omega)}{P(\omega)}$ by the Radon-Nikodym derivative
$\frac{dQ}{dP}(\omega).$:

$X=\log\left(\frac{Q(\omega)}{P(\omega)}\right)\text{ where }\omega\sim P,$
$Y=\log\left(\frac{Q(\omega)}{P(\omega)}\right)\text{ where }\omega\sim Q.$

The following theorem provides a formula for computing the privacy curve
$\delta$ in terms of the PRVs and conversely a formula for PRVs in terms of
the privacy curve. A similar statement appears in [SMM19, KJH+20].

######  Theorem 3.3.

The privacy curve $\delta$ can be expressed in terms of PRVs $(X,Y)$ as:

$\delta(\varepsilon)=\Pr[Y>\varepsilon]-e^{\varepsilon}\Pr[X>\varepsilon]=\operatorname{\mathbb{E}}_{Y}[(1-e^{\varepsilon-Y})_{+}]=\Pr[Y\geq\varepsilon+Z].$
(5)

where $Z$ is an exponential random variable.777For
$x\in\mathbb{R},x_{+}=\max\\{x,0\\}$. Conversely, given a privacy curve
$\delta:\mathbb{R}\to[0,1]$, we can compute the PDFs of its PRVs $(X,Y)$ as:

$Y(t)=\delta^{\prime\prime}(t)-\delta^{\prime}(t)\text{ and
}X(t)=e^{t}(\delta^{\prime\prime}(t)-\delta^{\prime}(t)).$ (6)

######  Remark 3.4.

Theorem 3.3 shows that the PRVs $X,Y$ do not depend on the particular $P,Q$
used to represent the privacy curve $\delta$ in Theorem 3.2. So we should
think of the PDF of of the PRV $Y$ (or $X$) as an equivalent reparametrization
of the privacy curve $\delta:\mathbb{R}\to[0,1]$, just as the notion of $f$-DP
[DRS19] is a reparametrization of the privacy curve $\delta$.

PRVs are useful in computing privacy curves because the composition of two
privacy curves can be computed by adding the corresponding pairs of PRVs. A
similar statement appears in [DR16].

######  Theorem 3.5.

Let $\delta_{1},\delta_{2}$ be two privacy curves with PRVs $(X_{1},Y_{1})$
and $(X_{2},Y_{2})$ respectively. Then the PRVs for
$\delta_{1}\otimes\delta_{2}=\delta(X_{1},X_{2}||Y_{1},Y_{2})$ are given by
$(X_{1}+X_{2},Y_{1}+Y_{2})$. In particular,

$\delta_{1}\otimes\delta_{2}=\delta(X_{1}+X_{2}||Y_{1}+Y_{2}).$

###### Proof.

Let $(X,Y)$ be the privacy random variables for
$\delta(X_{1},X_{2}||Y_{1},Y_{2})$. By Theorem 3.2,

$\displaystyle X$
$\displaystyle=\log\left(\frac{(Y_{1},Y_{2})(t_{1},t_{2})}{(X_{1},X_{2})(t_{1},t_{2})}\right)\text{
where }(t_{1},t_{2})\sim(X_{1},X_{2})$
$\displaystyle=\log\left(\frac{Y_{1}(t_{1})Y_{2}(t_{2})}{X_{1}(t_{1})X_{2}(t_{2})}\right)\text{
where }t_{1}\sim X_{1},t_{2}\sim X_{2}$ (By independence of $X_{1},X_{2}$ and
indpendence of $Y_{1},Y_{2}$) $\displaystyle=\log\left(e^{t_{1}}\cdot
e^{t_{2}}\right)\text{ where }t_{1}\sim X_{1},t_{2}\sim X_{2}$
$\displaystyle=t_{1}+t_{2}\text{ where }t_{1}\sim X_{1},t_{2}\sim X_{2}$
$\displaystyle=X_{1}+X_{2}.$

Similarly,

$\displaystyle Y$
$\displaystyle=\log\left(\frac{(Y_{1},Y_{2})(t_{1},t_{2})}{(X_{1},X_{2})(t_{1},t_{2})}\right)\text{
where }(t_{1},t_{2})\sim(Y_{1},Y_{2})$ $\displaystyle=t_{1}+t_{2}\text{ where
}t_{1}\sim Y_{1},t_{2}\sim Y_{2}$ $\displaystyle=Y_{1}+Y_{2}.$

∎

In Appendix B, we provide a proof of Theorems 3.2 and 3.3. We also discuss how
to compute the PRVs for a subsampled mechanism given the PRVs for the original
mechanism and give examples of PRVs for few standard mechanisms. These are
used in our experiments to calculate the PRVs for DP-SGD.

##  4 Numerical composition of privacy curves

In this section, we present an efficient and numerically accurate method,
ComposePRV (Algorithm 1), for composing privacy guarantees by utilizing the
notion of PRVs.

Input: $\mathrm{CDF}$s of PRVs $Y_{1},Y_{2},\dots,Y_{k}$, mesh size $h$,
Truncation parameter $L\in\frac{h}{2}+h\mathbb{Z}^{>0}$

Output: PDF of an approximation $\widetilde{Y}$ for $Y=\sum_{i=1}^{k}Y_{i}$.
$\widetilde{Y}$ will be supported on $\mu+(h\mathbb{Z}\cap[-L,L])$ for some
$\mu\in[0,\frac{h}{2}]$.

for  _$\ell=1$ to $k$_ do

$\widetilde{Y}_{i}\leftarrow\mathsf{DiscretizePRV}(Y_{i},L,h)$;

end for

Compute PDF Of
$\widetilde{Y}=\widetilde{Y}_{1}\oplus_{L}\widetilde{Y}_{2}\oplus_{L}\cdots\oplus_{L}\widetilde{Y}_{k}$
by convolving PDFs of
$\widetilde{Y}_{1},\widetilde{Y}_{2},\dots,\widetilde{Y}_{k}$ using FFT;

Compute
$\delta_{\widetilde{Y}}(\varepsilon)=\operatorname{\mathbb{E}}_{\widetilde{Y}}\left[\left(1-e^{\varepsilon-\widetilde{Y}}\right)_{+}\right]$
for all $\varepsilon\in[0,L]$;

Return $\widetilde{Y},\delta_{\widetilde{Y}}(\cdot)$

Algorithm 1 ComposePRV: Composing privacy curves using PRVs

In the algorithm ComposePRV, we compute the circular convolution $\oplus_{L}$
using Fast Fourier Transform (FFT). Fix some $L>0$. For $x\in\mathbb{R}$, we
define $x\pmod{2L}=x-2Ln$ where $n\in\mathbb{Z}$ is chosen such that
$x-2Ln\in(-L,L].$ Given $x,y$, we define the circular addition

$x\oplus_{L}y=x+y\pmod{2L}.$

When we use FFT to compute the convolution of two discrete distributions
$Y_{1},Y_{2}$ supported on $h\mathbb{Z}\cap[-L,L]$, we are implicitly
calculating the the distribution of $Y_{1}\oplus_{L}Y_{2}$. In the appendix,
we show that
$\widetilde{Y}_{1}\oplus_{L}\widetilde{Y}_{2}\oplus_{L}\dots\oplus_{L}\widetilde{Y}_{k}$
is a good approximation of $Y_{1}+Y_{2}+\dots+Y_{k}$.

The subroutine DiscretizePRV (Algorithm 2) is used to truncate and discretize
PRVs. In this subroutine, we shift the discretized random variables such that
it has the same mean as the original variables. This is one of main
differences between our algorithm and the algorithm in [KJPH21, KH21]. We show
that this significantly decreases the discretization error and allow us to use
much coarser mesh $h\approx 1/\sqrt{k}$ instead of $h\approx 1/k$.

Input: $\mathrm{CDF}_{Y}(\cdot)$ of a PRV $Y$, mesh size $h$, Truncation
parameter $L\in\frac{h}{2}+h\mathbb{Z}^{>0}$

Output: PDF of an approximation $\widetilde{Y}$ supported on
$\mu+(h\mathbb{Z}\cap[-L,L])$ for some $\mu\in[0,\frac{h}{2}].$

$n\leftarrow\frac{L-\frac{h}{2}}{h}$;

for  _$i=-n$ to $n$_ do

$q_{i}\leftarrow\mathrm{CDF}_{Y}(ih+h/2)-\mathrm{CDF}_{Y}(ih-h/2)$;

end for

$q\leftarrow q/\left(\sum_{i=-n}^{n}q_{i}\right)$ ;

// Normalize $q$ to make it a probability distribution

$Y^{L}\leftarrow Y\big{|}_{|Y|\leq L}$ (i.e., $Y$ conditioned on $|Y|\leq L$);

$\mu\leftarrow\operatorname{\mathbb{E}}[Y^{L}]-\sum_{i=-n}^{n}ih\cdot q_{i}$;

$\widetilde{Y}\leftarrow\begin{cases}ih+\mu&\text{ w.p. }q_{i}\text{ for
}-n\leq i\leq n\end{cases}$;

Return $\widetilde{Y}$;

Algorithm 2 DiscretizePRV: Discretize and truncate a PRV

For simplicity, throughout this paper, we will assume that the PRVs
$Y_{1},Y_{2},\dots,Y_{k}$ do not have any mass at $\infty$. This is with out
loss of generality. Suppose $\Pr[Y_{i}=\infty]=\delta_{i}$ for each $i.$ Let
$Y_{i}^{\prime}=Y_{i}|_{Y_{i}\neq\infty}$. Then

$\displaystyle
Y_{1}+Y_{2}+\dots+Y_{k}=\begin{cases}Y_{1}^{\prime}+Y_{2}^{\prime}+\dots+Y_{k}^{\prime}&\text{
w.p. }(1-\delta_{1})(1-\delta_{2})\cdots(1-\delta_{k})\\\ \infty&\text{ w.p.
}1-(1-\delta_{1})(1-\delta_{2})\cdots(1-\delta_{k}).\end{cases}$

Therefore we can use Algorithm 1 to approximate the distribution of
$Y_{1}^{\prime}+Y_{2}^{\prime}+\dots+Y_{k}^{\prime}$, and use it to
approximate the distribution of $Y_{1}+Y_{2}+\dots+Y_{k}$.

##  5 Error analysis

To analyze the discretization error, we introduce the notion of _coupling
approximation_ , a variant of Wasserstein distance. Intuitively, a good
coupling approximation is a coupling where the two random variables are close
to each other with high probability.

######  Definition 5.1 (coupling approximation).

Given two random variables $Y_{1},Y_{2}$, we write $|Y_{1}-Y_{2}|\leq_{\eta}h$
if there exists a coupling between $Y_{1},Y_{2}$ such that
$\Pr[|Y_{1}-Y_{2}|>h]\leq\eta.$

The following lemma shows that if we have a good coupling approximation
$\widetilde{Y}$ to a PRV $Y$, then the privacy curves
$\delta_{Y}(\varepsilon)$ and $\delta_{\widetilde{Y}}(\varepsilon)$ should be
close.

######  Lemma 5.2.

If $Y$ and $\widetilde{Y}$ are two random variables such that
$|Y-\widetilde{Y}|\leq_{\eta}h$, then for every $\varepsilon\in\mathbb{R}$,

$\delta_{\widetilde{Y}}(\varepsilon+h)-\eta\leq\delta_{Y}(\varepsilon)\leq\delta_{\widetilde{Y}}(\varepsilon-h)+\eta.$

###### Proof.

By Theorem 3.2, $\delta_{Y}(\varepsilon)=\Pr[Y\geq\varepsilon+Z]$ and hence

$\displaystyle\delta_{Y}(\varepsilon)$
$\displaystyle=\Pr[Y-\widetilde{Y}+\widetilde{Y}\geq\varepsilon+Z]$
$\displaystyle\leq\Pr[Y-\widetilde{Y}\geq
h]+\Pr[\widetilde{Y}\geq\varepsilon-h+Z]$
$\displaystyle\leq\eta+\delta_{\widetilde{Y}}(\varepsilon-h).$

Similarly, we have
$\delta_{\widetilde{Y}}(\varepsilon)\leq\eta+\delta_{Y}(\varepsilon-h)$ for
all $\varepsilon\in\mathbb{R}.$ ∎

Therefore the goal of our analysis is to show that the ComposePRV algorithm
finds a good coupling approximation $\widetilde{Y}$ to
$Y=\sum_{i=1}^{k}Y_{i}.$ We first show that the DiscretizePRV algorithm
computes a good coupling approximation to the PRVs and crucially, it preserves
the expected value after truncation. Lemma C.5 shows that
$|\widetilde{Y}-Y^{L}|\leq_{0}h$ where $\widetilde{Y}$ is the approximation of
a PRV $Y$ output by Algorithm 2 and $Y^{L}$ is the truncation of $Y$ to
$[-L,L].$

We then use the following key lemma which shows that when we add independent
coupling approximations (where expected values match), we get a much better
coupling approximation than what the triangle inequality predicts.

######  Lemma 5.3.

Suppose $Y_{1},Y_{2},\dots,Y_{k}$ and
$\widetilde{Y}_{1},\widetilde{Y}_{2},\dots,\widetilde{Y}_{k}$ are two
collections of independent random variables such that
$|Y_{i}-\widetilde{Y}_{i}|\leq_{0}h$ and
$\operatorname{\mathbb{E}}[Y_{i}]=\operatorname{\mathbb{E}}[\widetilde{Y}_{i}]$
for all $i$, then

$\left|\sum_{i=1}^{k}Y_{i}-\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\leq_{\eta}h\sqrt{2k\log{\frac{2}{\eta}}}.$

###### Proof.

Let $X_{i}=Y_{i}-\widetilde{Y}_{i}$ where $(Y_{i},\widetilde{Y}_{i})$ are
coupled such that $|Y_{i}-\widetilde{Y}_{i}|\leq h$ w.p. $1$. Then
$X_{i}\in[-h,h]$ w.p. $1$. Note that $X_{1},X_{2},\dots,X_{k}$ are independent
of each other. By Hoeffding’s inequality,

$\Pr\left[\left|\sum_{i}X_{i}\right|\geq t\right]\leq
2\exp\left(-\frac{2t^{2}}{k(2h)^{2}}\right)=\eta$

if we set $t=h\sqrt{2k\log{\frac{2}{\eta}}}$. ∎

This lemma shows that the error of $k$ times composition is around
$\sqrt{k}\cdot h$ and hence setting $h\approx 1/\sqrt{k}$ gives small enough
error. Next, we bound the domain size $L$. Naively, the domain size $L$ should
be of the order of $\sqrt{k}$ because $Y$ is the sum of $k$ independent random
variables with each bounded by a constant. In the appendix, we give a tighter
tail bound of $Y$.

######  Lemma 5.4.

Let $(X,Y)$ be the privacy random variables for a $(\varepsilon,\delta)$-DP
algorithm, then for any $t\geq 0$, we have

$\Pr[|Y|\geq\varepsilon+t]\leq\frac{\delta\left(1+e^{-\varepsilon-t}\right)}{1-e^{-t}}.$

This shows that $\Pr[|Y|\geq\varepsilon+2]\leq\frac{4}{3}\delta$ and hence
truncating the domain with $L=2+\varepsilon$ only introduces an additive
$\delta$ error in the privacy curve. Therefore, if the composition satisfies a
good privacy guarantee (namely $\varepsilon=O(1)$ for small enough $\delta$),
we can truncate the domain at $L=\Theta(1)$. Together with the fact that mesh
size is $1/\sqrt{k}$, this gives a $O(\sqrt{k})$-time algorithm for computing
the privacy curve when we compose the same mechanism with itself $k$ times.
The following theorem gives a formal statement of the error bounds of our
algorithm, it is proved in Appendix 5.

######  Theorem 5.5.

Let $\varepsilon_{\mathrm{error}},\delta_{\mathrm{error}}>0$ be some fixed
error terms. Let $\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{k}$ be DP
algorithms with privacy curves $\delta_{\mathcal{M}_{i}}(\varepsilon)$. Let
$Y_{i}$ be the PRV corresponding to $\mathcal{M}_{i}$ such that
$\delta_{\mathcal{M}_{i}}(\varepsilon)=\delta_{Y_{i}}(\varepsilon)$ for
$\varepsilon\geq 0$. Let $\mathcal{M}$ be the (adaptive) composition of
$\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{k}$ and let
$\delta_{\mathcal{M}}(\varepsilon)$ be its privacy curve. Set $L\geq
2+\varepsilon_{\mathrm{error}}$ sufficiently large such that

$\sum_{i=1}^{k}\delta_{\mathcal{M}_{i}}(L-2)\leq\frac{\delta_{\mathrm{error}}}{8}\text{
and
}\delta_{\mathcal{M}}(L-2-\varepsilon_{\mathrm{error}})\leq\frac{\delta_{\mathrm{error}}}{4}.$
(7)

Let $\widetilde{Y}$ be the approximation of $Y=\sum_{i=1}^{k}Y_{i}$ produced
by ComposePRV algorithm with mesh size

$h=\frac{\varepsilon_{\mathrm{error}}}{\sqrt{\frac{k}{2}\log\frac{12}{\delta_{\mathrm{error}}}}}.$

Then

$\delta_{\widetilde{Y}}(\varepsilon+\varepsilon_{\mathrm{error}})-\delta_{\mathrm{error}}\leq\delta_{Y}(\varepsilon)=\delta_{\mathcal{M}}(\varepsilon)\leq\delta_{\widetilde{Y}}(\varepsilon-\varepsilon_{\mathrm{error}})+\delta_{\mathrm{error}}.$
(8)

Furthermore, our algorithm takes
$O\left(b\frac{L}{h}\log\left(\frac{L}{h}\right)\right)$ time where $b$ is the
number of distinct algorithms among
$\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{k}$.

######  Remark 5.6.

A simple way to set $L$ such that the condition (7) holds is by choosing an
$L$ such that:

$L\geq
2+\max\left\\{\varepsilon_{\mathrm{error}}+\varepsilon_{\mathcal{M}}\left(\frac{\delta_{\mathrm{error}}}{4}\right),\max_{i\in[k]}\
\varepsilon_{\mathcal{M}_{i}}\left(\frac{\delta_{\mathrm{error}}}{8k}\right)\right\\}$
(9)

where $\varepsilon_{\mathcal{A}}(\delta)$ is the inverse of
$\delta_{\mathcal{A}}(\varepsilon)$. To set the value of $L$, we do not need
the exact value of $\varepsilon_{\mathcal{M}}$ (or
$\varepsilon_{\mathcal{M}_{i}}$). We only need an upper bound on
$\varepsilon_{\mathcal{M}}$, which can often by obtained by using the RDP
Accountant or any other method to derive upper bounds on privacy.

##  6 Experiments

In this section, we demonstrate the utility of our composition method by
computing the privacy curves for the DP-SGD algorithm which is one of the most
important algorithms in differential privacy.

The DP-SGD algorithm [ACG+16] is a variant of stochastic gradient descent with
$k$ steps. In each step, the algorithm selects a $p$ fraction of training
examples uniformly at random. The algorithm adds a Gaussian vector with
variance $\propto\sigma^{2}$ to the clipped gradient of the selected batch.
Then it performs a gradient step (or any other iterative methods) using the
noisy gradient computed. The privacy loss of DP-SGD involves composing the
privacy curve of each iteration with itself $k$ times. The PRVs for each
iteration have a closed form and depend only $p,\sigma$ (see Appendix). Our
algorithms use this closed form of PRVs.

See Figure 1(b) for the comparison between our algorithm and the GDP and RDP
Accountant. Our method provides a lower and upper bound of the privacy curve
according to (8). In Figure 1(a), we compare our algorithm with [KJPH21]
(implemented in [KP21]). Under the same mesh size, our algorithm computes a
much closer upper and lower bound.

We validate our program for the case $p=1$. When $p=1$, we have an exact
formula for

$\delta(\varepsilon)=\Phi\left(-\frac{\varepsilon}{\mu}+\frac{\mu}{2}\right)-e^{\varepsilon}\Phi\left(-\frac{\varepsilon}{\mu}-\frac{\mu}{2}\right)$
(10)

where $\mu=\frac{\sqrt{k}}{\sigma}$. In Figure 2, we show that the true
privacy curve is indeed sandwiched between the bounds we compute and that the
vertical distance between our bounds is indeed $2\varepsilon_{\mathrm{error}}$
with a neglible $\delta_{\mathrm{error}}$ of $10^{-10}$.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$

$\varepsilon$

$\delta$$\varepsilon_{\textnormal{exact}}$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=0.1$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=0.1$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=0.5$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=0.5$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=1.0$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=1.0$

Figure 2: Setting $p=1$ and comparing to the analytical solution (10).

#### Floating point errors

Note that our error analysis in Section 5 ignores floating point errors. This
is because they are negligible compared to the discretization and truncation
errors we analyzed in Section 5 for the range of $\delta$ we are interested
in. Our implementation uses `ong } {\verb doub`e floating point format which
is platform dependent, however, it guarantees a precision at least as good as
double precision which has a resolution of $10^{-15}$. Computations involving
$\delta$ of these orders of magnitude suffer from floating point inaccuracies.
Our implementation therefore only allows $\delta$ values which are greater
than $10^{-10}$ which sufficies for practical use cases. See Appendix A for
more details.

###  6.1 Comparison with [KJPH21]

$10^{-10}$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{3}$$10^{4}$$10^{5}$$10^{6}$$10^{7}$$10^{8}$$10^{9}$

$\delta$

Discretisation
points$\delta_{\textnormal{up}}$[KJPH21]$\delta_{\textnormal{low}}$[KJPH21]$\delta_{\textnormal{up}}$Ours$\delta_{\textnormal{low}}$Ours

Figure 3: Comparison of error bounds of $\delta$ with varying number of
discretisation points for $p=4\times
10^{-3},\sigma=0.8,\varepsilon=1.5,k=1000$.

$0$$0.5$$1$$1.5$$2$$10^{-9}$$10^{-6}$$10^{-3}$1$\delta$vs
$\varepsilon$$\varepsilon_{\textnormal{up}}$[KJPH21]$\varepsilon_{\textnormal{low}}$[KJPH21]$\varepsilon_{\textnormal{up}}$Ours$\varepsilon_{\textnormal{low}}$Ours

(a) $k=10$

$0$$1$$2$$3$$4$$10^{-9}$$10^{-6}$$10^{-3}$1$\delta$vs $\varepsilon$

(b) $k=100$

$0$$1$$2$$3$$4$$5$$10^{-9}$$10^{-6}$$10^{-3}$1$\delta$vs $\varepsilon$

(c) $k=1000$

Figure 4:  Comparing different error bounds using the same mesh size $8\times
10^{-4}$ under different number of steps $k=10,100,1000$. (With $p=10^{-2}$,
$\sigma=0.8$.)

$0$$2$$4$$6$$8$$10$$12$$100$$1000$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

Runtime [s]

$\times$Speed up

DPSGD stepsOursKJPH21Speed up

$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$0$$500$$1000$$1500$$2000$

$\delta$

DPSGD stepsOurs, UpperOurs, LowerKJPH21, UpperKJPH21, Lower

Figure 5: (a) Comparing runtimes for our algorithm with that of [KJPH21] when
aligned on accuracy for $\sigma=0.8$, $p=4\times 10^{-3}$. We can see a
significant reduction in runtime in particular for large number of DPSGD
steps. We were not able to run the algorithm of [KJPH21] beyond 2,000 steps,
since it becomes unstable beyond that point.999We are using the implementation
of [KJPH21] from [KP21].. We also plot the speed up directly on the secondary
$y$-axis. (b) Verification of the alignment of the error bounds of both
algorithms at $\varepsilon=1.5$.

In this section, we provide more results demonstrating the practical use of
our algorithm. We compare runtimes of our algorithm with [KJPH21], which is
the state-of-the-art, for typical values of privacy parameters ($\sigma=0.8$,
$p=4\times 10^{-3}$, $\varepsilon=1.5$).

See Figure 3 for the effect of the number of discretisation points $n$ on the
accuracy of $\delta$. Our algorithm requires about a few orders of magnitude
smaller number of discretization points to converge compared to the algorithm
of [KJPH21]. A similar picture can be seen in Figure 4. While for a small
number of compositions, the algorithm of [KJPH21] gives reasonable estimates,
for a large number of compositions, their error bounds worsen quickly.

We note that runtimes are directly proportional to the memory required by the
algorithms and so a separate memory analysis is not required; the runtime and
memory are dominated by the number of points in the discretization of PRV. All
experiments are performed on a Intel Xeon W-2155 CPU with 3.30GHz with 128GB
of memory.

In order to compare runtimes, we align the accuracy of both FFT algorithms. We
find sets of numerical parameters (number of discretization bins and domain
length) such that both algorithms give similarly accurate bounds and verify it
visually (see Figure 9 (b)). Figure 9 illustrates the runtimes for varying
numbers of DPSGD steps. We observe a significant reduction in the runtime
using our algorithms.

### Acknowledgements

We would like to thank Janardhan Kulkarni and Sergey Yekhanin for several
useful discussions and encouraging us to work on this problem. L.W. would like
to thank Daniel Jones and Victor Rühle for fruitful discussions and helpful
guidance.

## References

  * [ACG+16] Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang.  Deep learning with differential privacy.  In Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, pages 308–318, 2016.
  * [BDLS19] Zhiqi Bu, Jinshuo Dong, Qi Long, and Weijie J. Su.  Deep learning with gaussian differential privacy, 2019.
  * [BDRS18] Mark Bun, Cynthia Dwork, Guy N Rothblum, and Thomas Steinke.  Composable and versatile privacy via truncated cdp.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 74–86, 2018.
  * [BGK+21] Zhiqi Bu, Sivakanth Gopi, Janardhan Kulkarni, Yin Tat Lee, Judy Hanwen Shen, and Uthaipon Tantipongpipat.  Fast and memory efficient differentially private-sgd via jl projections.  arXiv preprint arXiv:2102.03013, 2021.
  * [BS16] Mark Bun and Thomas Steinke.  Concentrated differential privacy: Simplifications, extensions, and lower bounds.  In Theory of Cryptography Conference, pages 635–658. Springer, 2016\.
  * [BW18] Borja Balle and Yu-Xiang Wang.  Improving the gaussian mechanism for differential privacy: Analytical calibration and optimal denoising.  In International Conference on Machine Learning, pages 403–412, 2018.
  * [DKM+06] Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor.  Our data, ourselves: Privacy via distributed noise generation.  In Annual International Conference on the Theory and Applications of Cryptographic Techniques, pages 486–503. Springer, 2006.
  * [DL09] Cynthia Dwork and Jing Lei.  Differential privacy and robust statistics.  In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 371–380, 2009.
  * [DMNS06] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith.  Calibrating noise to sensitivity in private data analysis.  In Theory of Cryptography Conference, pages 265–284. Springer, 2006\.
  * [DR+14] Cynthia Dwork, Aaron Roth, et al.  The algorithmic foundations of differential privacy.  Foundations and Trends® in Theoretical Computer Science, 9(3–4):211–407, 2014.
  * [DR16] Cynthia Dwork and Guy N Rothblum.  Concentrated differential privacy.  arXiv preprint arXiv:1603.01887, 2016.
  * [DRS19] Jinshuo Dong, Aaron Roth, and Weijie J. Su.  Gaussian differential privacy, 2019.
  * [DRV10] Cynthia Dwork, Guy N Rothblum, and Salil Vadhan.  Boosting and differential privacy.  In 2010 IEEE 51st Annual Symposium on Foundations of Computer Science, pages 51–60. IEEE, 2010.
  * [KH21] Antti Koskela and Antti Honkela.  Computing differential privacy guarantees for heterogeneous compositions using fft.  arXiv preprint arXiv:2102.12412, 2021.
  * [KJH+20] Antti Koskela, Joonas Jälkö, Antti Honkela, et al.  Computing tight differential privacy guarantees using fft.  In Proceedings of the Twenty Third 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 Arindam Banerjee and Kenji Fukumizu, editors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 3358–3366. PMLR, 13–15 Apr 2021.
  * [KOV15] Peter Kairouz, Sewoong Oh, and Pramod Viswanath.  The composition theorem for differential privacy.  In International conference on machine learning, pages 1376–1385. PMLR, 2015.
  * [KP21] Antti Koskela and Lukas Prediger.  Github repository for fourier accountant.  https://github.com/DPBayes/PLD-Accountant, 2021.
  * [Mir17] Ilya Mironov.  Rényi differential privacy.  In 2017 IEEE 30th Computer Security Foundations Symposium (CSF), pages 263–275. IEEE, 2017.
  * [MM18] Sebastian Meiser and Esfandiar Mohammadi.  Tight on budget? tight bounds for r-fold approximate differential privacy.  In Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security, pages 247–264, 2018.
  * [MV16] Jack Murtagh and Salil Vadhan.  The complexity of computing the optimal composition of differential privacy.  In Theory of Cryptography Conference, pages 157–175. Springer, 2016\.
  * [SMM19] David M Sommer, Sebastian Meiser, and Esfandiar Mohammadi.  Privacy loss classes: The central limit theorem in differential privacy.  Proceedings on privacy enhancing technologies, 2019(2):245–269, 2019.
  * [ZDW21] Yuqing Zhu, Jinshuo Dong, and Yu-Xiang Wang.  Optimal accounting of differential privacy via characteristic function.  arXiv preprint arXiv:2106.08567, 2021.

##  Appendix A Effect of floating point arithmetic

In this section, we demonstrate the effect of floating point inaccuracies on
the computed privacy parameters. Figure 6 compares lower and upper bounds of
the privacy curve with the analytical solution for small values of $\delta$.
As mentioned in section 6, we use a floating point representation with a
resolution of at least $10^{-15}$. The number of discretization points in this
examples are on the order of $10^{4}$. Consequently, we expect floating point
inaccuracies to become dominant for values on the order of $10^{-11}$. This
can be also seen in the illustration, where the lower and upper bound fail to
produce meaningful results for $\delta<2\times 10^{-11}$.

$5.5$$6$$6.5$$7$$7.5$$8$$10^{-13}$$10^{-11}$$10^{-9}$$10^{-7}$$10^{-5}$

$\varepsilon$

$\delta$$\varepsilon_{\textnormal{exact}}$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=0.1$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=0.1$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=0.5$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=0.5$$\varepsilon_{\textnormal{up}}$,
$\varepsilon_{\textnormal{shift}}=1.0$$\varepsilon_{\textnormal{low}}$,
$\varepsilon_{\textnormal{shift}}=1.0$

Figure 6: Setting $p=1$ and comparing to the analytical solution (10) for
values of $\delta$ beyond expected floating point accuracy.

##  Appendix B Privacy Loss Random Variables

In this section, we continue the discussion on privacy random variables in
Section 3. First, we give the proof of the formula for PRVs of $\delta(P||Q)$
and the formula for a privacy curve given its PRVs (Theorem 3.2).

See 3.2

###### Proof.

We will first verify that $Y(t)=e^{t}X(t).$ This is equivalent to proving that
$\operatorname{\mathbb{E}}_{Y}[\phi(Y)]=\operatorname{\mathbb{E}}_{X}[\phi(X)e^{X}]$
for any test function $\phi:\overline{\mathbb{R}}\to[0,1]$. This is true since

$\displaystyle\operatorname{\mathbb{E}}_{Y}[\phi(Y)]$
$\displaystyle=\operatorname{\mathbb{E}}_{\omega\sim
Q}\left[\phi\left(\log\left(\frac{Q(\omega)}{P(\omega)}\right)\right)\right]$
$\displaystyle=\operatorname{\mathbb{E}}_{\omega\sim
P}\left[\phi\left(\log\left(\frac{Q(\omega)}{P(\omega)}\right)\right)\frac{Q(\omega)}{P(\omega)}\right]$
$\displaystyle=\operatorname{\mathbb{E}}_{X}[\phi(X)e^{X}].$

We will now prove that $\delta(X||Y)=\delta(P||Q).$ We have

$\displaystyle\delta(P||Q)(\varepsilon)$
$\displaystyle=\sup_{S\subset\Omega}\Pr[Q\in S]-e^{\varepsilon}\Pr[P\in S]$
$\displaystyle=\Pr[Q\in S_{\varepsilon}]-e^{\varepsilon}\Pr[P\in
S_{\varepsilon}]$

where

$S_{\varepsilon}=\left\\{\omega\in\Omega:\frac{Q(\omega)}{P(\omega)}>e^{\varepsilon}\right\\}=\left\\{\omega\in\Omega:\log\left(\frac{Q(\omega)}{P(\omega)}\right)>\varepsilon\right\\}.$

Therefore $\Pr[Q\in S_{\varepsilon}]=\Pr[Y>\varepsilon]$ and $\Pr[P\in
S_{\varepsilon}]=\Pr[X>\varepsilon]$. To complete the proof, note that

$\displaystyle\delta(X||Y)(\varepsilon)$
$\displaystyle=\sup_{T\subset\overline{\mathbb{R}}}\Pr[Y\in
T]-e^{\varepsilon}\Pr[X\in T]$ $\displaystyle=\Pr[Y\in
T_{\varepsilon}]-e^{\varepsilon}\Pr[X\in T_{\varepsilon}]$

where

$T_{\varepsilon}=\left\\{t\in\overline{\mathbb{R}}:\frac{Y(t)}{X(t)}>e^{\varepsilon}\right\\}=\left\\{t\in\overline{\mathbb{R}}:e^{t}>e^{\varepsilon}\right\\}=(\varepsilon,\infty].$

Putting it all together, we have:

$\displaystyle\delta(P||Q)(\varepsilon)$
$\displaystyle=\Pr[Y>\varepsilon]-e^{\varepsilon}\Pr[X>\varepsilon]=\delta(X||Y).$

∎

See 3.3

###### Proof.

Since the PDFs of PRVs $(X,Y)$ satisfy the relation $Y(t)=e^{t}X(t)$, we can
rewrite the equation 5 in terms of just $Y$ or just $X$.

$\displaystyle\delta(\varepsilon)$
$\displaystyle=\Pr[Y\geq\varepsilon]-e^{\varepsilon}\Pr[X\geq\varepsilon]$
$\displaystyle=\int_{\varepsilon}^{\infty}Y(t)dt-\int_{\varepsilon}^{\infty}e^{\varepsilon}X(t)dt$
$\displaystyle=\int_{\varepsilon}^{\infty}Y(t)dt-\int_{\varepsilon}^{\infty}e^{\varepsilon-t}Y(t)dt$
(Since $Y(t)=e^{t}X(t)$)
$\displaystyle=\int_{\varepsilon}^{\infty}Y(t)(1-e^{\varepsilon-t})dt$
$\displaystyle=\int_{-\infty}^{\infty}Y(t)(1-e^{\varepsilon-t})_{+}dt$
$\displaystyle=\operatorname{\mathbb{E}}_{Y}[(1-e^{\varepsilon-Y})_{+}]$

To get the other form for $\delta(\varepsilon)$, we use the integration by
parts formula.

$\displaystyle\delta(\varepsilon)$
$\displaystyle=\int_{\varepsilon}^{\infty}Y(t)(1-e^{\varepsilon-t})dt$
$\displaystyle=\int_{\varepsilon}^{\infty}Y(t)dt+\int_{\varepsilon}^{\infty}\left(-Y(t)\right)e^{\varepsilon-t}dt$
$\displaystyle=\Pr[Y\geq\varepsilon]+\left(\Pr[Y\geq
t]e^{\varepsilon-t}\Big{|}_{\varepsilon}^{\infty}-\int_{\varepsilon}^{\infty}\Pr[Y\geq
t]\left(-e^{\varepsilon-t}\right)dt\right)$
$\displaystyle=\Pr[Y\geq\varepsilon]-\Pr[Y\geq\varepsilon]+\int_{\varepsilon}^{\infty}\Pr[Y\geq
t]e^{\varepsilon-t}dt$
$\displaystyle=\int_{\varepsilon}^{\infty}e^{\varepsilon-t}\Pr[Y\geq t]dt$
$\displaystyle=\int_{0}^{\infty}e^{-z}\Pr[Y\geq\varepsilon+z]dz$ (Substituting
$z=t-\varepsilon$) $\displaystyle=\Pr[Y\geq\varepsilon+Z].$ (where $Z$ is an
exponential random variable)

We now prove the converse relation by differentiating the expression for
$\delta(\varepsilon)$ twice. We have:

$\displaystyle\delta(\varepsilon)=\int_{\varepsilon}^{\infty}Y(t)dt-e^{\varepsilon}\int_{\varepsilon}^{\infty}e^{-t}Y(t)dt$
$\displaystyle\implies$
$\displaystyle\delta^{\prime}(\varepsilon)=-Y(\varepsilon)+e^{\varepsilon}\cdot
e^{-\varepsilon}Y(\varepsilon)-e^{\varepsilon}\cdot\int_{\varepsilon}^{\infty}e^{-t}Y(t)dt=-e^{\varepsilon}\cdot\int_{\varepsilon}^{\infty}e^{-t}Y(t)dt$
$\displaystyle\implies$ $\displaystyle
e^{-\varepsilon}\delta^{\prime}(\varepsilon)=-\int_{\varepsilon}^{\infty}e^{-t}Y(t)dt$
$\displaystyle\implies$ $\displaystyle
e^{-\varepsilon}\delta^{\prime\prime}(\varepsilon)-e^{-\varepsilon}\delta^{\prime}(\varepsilon)=e^{-\varepsilon}Y(\varepsilon)$
$\displaystyle\implies$ $\displaystyle
Y(\varepsilon)=\delta^{\prime\prime}(\varepsilon)-\delta^{\prime}(\varepsilon).$

∎

###  B.1 Examples of privacy loss random variables

In this section, we state the PRVs for a few standard mechanisms.

######  Proposition B.1 (Gaussian Mechanism).

The PRVs for $\delta(\mathcal{N}(\mu,1)||\mathcal{N}(0,1))$ are:

$X=\mathcal{N}(-\mu^{2}/2,\mu^{2})\text{ and
}Y=\mathcal{N}(\mu^{2}/2,\mu^{2}).$

###### Proof.

Let $P=\mathcal{N}(\mu,1)$ and $Q=\mathcal{N}(0,1)$. By Theorem 3.2,

$\displaystyle Y$ $\displaystyle\sim\log\left(\frac{Q(t)}{P(t)}\right)\text{
where }t\sim Q$
$\displaystyle\sim\log\left(\frac{\exp(-t^{2}/2)}{\exp(-(t-\mu)^{2}/2)}\right)\text{
where }t\sim\mathcal{N}(0,1)$
$\displaystyle\sim\frac{(t-\mu)^{2}}{2}-\frac{t^{2}}{2}\text{ where
}t\sim\mathcal{N}(0,1)$ $\displaystyle\sim\frac{\mu^{2}}{2}-\mu t\text{ where
}t\sim\mathcal{N}(0,1)$
$\displaystyle=\mathcal{N}\left(\frac{\mu^{2}}{2},\mu^{2}\right).$

A similar calculation shows that
$X=\mathcal{N}\left(-\frac{\mu^{2}}{2},\mu^{2}\right)$ ∎

######  Proposition B.2 (Laplace Mechanism).

The PRVs for the privacy curve
$\delta\left(\mathsf{Lap}\left(\mu,1\right)||\mathsf{Lap}\left(0,1\right)\right)$
are:

$X=|Z|-|Z-\mu|\text{ and }Y=|Z-\mu|-|Z|$

where $Z\sim\mathsf{Lap}(0,1).$

###### Proof.

Let $P=\mathsf{Lap}(\mu,1)$ and $Q=\mathsf{Lap}(0,1)$. By Theorem 3.2,

$\displaystyle Y$ $\displaystyle\sim\log\left(\frac{Q(t)}{P(t)}\right)\text{
where }t\sim Q$
$\displaystyle\sim\log\left(\frac{\exp(-|t|)}{\exp(-|t-\mu|}\right)\text{
where }t\sim\mathsf{Lap}(0,1)$ $\displaystyle\sim|t-\mu|-|t|\text{ where
}t\sim\mathsf{Lap}(0,1)$ $\displaystyle=|Z-\mu|-|Z|\text{ where
}Z\sim\mathsf{Lap}(0,1).$

A similar calculation shows that $X=|Z|-|Z-\mu|$ where
$Z\sim\mathsf{Lap}(0,1).$ ∎

######  Proposition B.3 ($(\varepsilon,\delta)$-DP).

The PRVs for the privacy curve of a $(\varepsilon,\delta)$-DP algorithm are

$X=\begin{cases}-\infty&\text{w.p. }\delta\\\ -\varepsilon&\text{w.p.
}\frac{(1-\delta)e^{\varepsilon}}{e^{\varepsilon}+1}\\\ \varepsilon&\text{w.p.
}\frac{1-\delta}{e^{\varepsilon}+1},\end{cases}$
$Y=\begin{cases}-\varepsilon&\text{w.p. }\frac{1-\delta}{e^{\varepsilon}+1}\\\
\varepsilon&\text{w.p. }\frac{(1-\delta)e^{\varepsilon}}{e^{\varepsilon}+1}\\\
\infty&\text{w.p. }\delta.\end{cases}$

###### Proof.

It is easy to verify that $Y(t)=e^{t}X(t)$ for all $t\in\mathbb{R}$. We can
also verify that

$\delta(\varepsilon)=\Pr[Y>\varepsilon]-e^{\varepsilon}\Pr[X>\varepsilon]=\delta.$

Morever $X=-Y$, therefore the privacy curve $\delta(X||Y)$ is symmetric by
Proposition C.9, i.e., $\delta(X||Y)=\delta(Y||X)$. These conditions together
imply that $X,Y$ are PRVs for the $(\varepsilon,\delta)$-DP curve. ∎

Note that in the all the above examples, we have $X=-Y$ as the privacy curves
are symmetric.

###  B.2 Subsampling

In this section, we calculate the PRVs for a subsampled mechanism given the
PRVs for the original mechanism. Given two random variables $P,Q$ and a
sampling probability $p\in[0,1]$, $p\cdot P+(1-p)\cdot Q$ denotes the mixture
where we sample $P$ w.p. $p$ and $Q$ w.p. $1-p.$

######  Proposition B.4.

Let $(X,Y)$ be the PRVs for a privacy curve $\delta(P||Q)$. Let
$(X_{p},Y_{p})$ be the PRVs for $\delta_{p}=\delta(P||\ p\cdot P+(1-p)\cdot
Q)$. Then

$\displaystyle X_{p}$ $\displaystyle=\log(1+p(e^{X}-1)),$ $\displaystyle
Y_{p}$ $\displaystyle=\begin{cases}\log(1+p(e^{Y}-1))\text{ w.p. }p\\\
\log(1+p(e^{X}-1))\text{ w.p. }1-p.\end{cases}$

The CDFs of $X_{p}$ and $Y_{p}$ are given by:

$\displaystyle\mathrm{CDF}_{X_{p}}(t)$
$\displaystyle=\begin{cases}\mathrm{CDF}_{X}\left(\log\left(\frac{e^{t}-(1-p)}{p}\right)\right)&\text{
if }t\geq\log(1-p)\\\ 0&\text{ if }t<\log(1-p)\end{cases}$
$\displaystyle\mathrm{CDF}_{Y_{p}}(t)$
$\displaystyle=\begin{cases}p\cdot\mathrm{CDF}_{Y}\left(\log\left(\frac{e^{t}-(1-p)}{p}\right)\right)+(1-p)\cdot\mathrm{CDF}_{X}\left(\log\left(\frac{e^{t}-(1-p)}{p}\right)\right)&\text{
if }t\geq\log(1-p)\\\ 0&\text{ if }t<\log(1-p).\end{cases}$

###### Proof.

By Theorem 3.2,

$\displaystyle X_{p}$
$\displaystyle=\log\left(\frac{pY(t)+(1-p)X(t)}{X(t)}\right)\text{ where
}t\sim X$ $\displaystyle=\log\left(pe^{t}+1-p\right)\text{ where }t\sim X$
$\displaystyle=\log\left(pe^{X}+1-p\right).$

Similarly,

$\displaystyle Y_{p}$
$\displaystyle=\log\left(\frac{pY(t)+(1-p)X(t)}{X(t)}\right)\text{ where
}t\sim pY+(1-p)X$ $\displaystyle=\log\left(pe^{t}+1-p\right)\text{ where
}t\sim pY+(1-p)X$ $\displaystyle=\begin{cases}\log(1+p(e^{Y}-1))\text{ w.p.
}p\\\ \log(1+p(e^{X}-1))\text{ w.p. }1-p.\end{cases}$

The CDF of $X_{p}$ is given by:

$\displaystyle\Pr[X_{p}\leq t]$
$\displaystyle=\Pr\left[\log\left(pe^{X}+1-p\right)\leq t\right]$
$\displaystyle=\Pr\left[X\leq\log\left(\frac{e^{t}-(1-p)}{p}\right)\right]$

The CDF of $Y_{p}$ is given by:

$\displaystyle\Pr[Y_{p}\leq t]$
$\displaystyle=p\Pr\left[\log\left(pe^{Y}+1-p\right)\leq
t\right]+(1-p)\Pr\left[\log\left(pe^{X}+1-p\right)\leq t\right]$
$\displaystyle=p\Pr\left[Y\leq\log\left(\frac{e^{t}-(1-p)}{p}\right)\right]+(1-p)\Pr\left[X\leq\log\left(\frac{e^{t}-(1-p)}{p}\right)\right].$

∎

##  Appendix C Missing Proofs in Error Analysis

###  C.1 Facts about Coupling Approximation

Here we collect some useful properties of coupling approximations. The
following lemma shows that the coupling approximations satisfy a triangle
inequality.

######  Lemma C.1 (Triangle inequality for couplings).

Suppose $X,Y,Z$ are random variables such that $|X-Y|\leq_{\eta_{1}}h_{1}$ and
$|Y-Z|\leq_{\eta_{2}}h_{2}$. Then $|X-Z|\leq_{\eta_{1}+\eta_{2}}h_{1}+h_{2}.$

###### Proof.

There exists couplings $(X,Y)$ and $(Y,Z)$ such that

$\Pr[|X-Y|\geq h_{1}]\leq\eta_{1}\text{ and }\Pr[|Y-Z|\geq
h_{2}]\leq\eta_{2}.$

From these two couplings, we can construct a coupling between $(X,Z)$: sample
$X$, sample $Y$ from $Y|X$ (given by coupling $(X,Y)$) and finally sample $Z$
from $Z|Y$ (given by coupling $(Y,Z)$). Therefore for this coupling, we have:

$\displaystyle\Pr[|X-Z|\geq h_{1}+h_{2}]$
$\displaystyle\leq\Pr[|(X-Y)+(Y-Z)|\geq h_{1}+h_{2}]$
$\displaystyle\leq\Pr[|X-Y|+|Y-Z|\geq h_{1}+h_{2}]$
$\displaystyle\leq\Pr[|X-Y|\geq h_{1}]+\Pr[|Y-Z|\geq h_{2}]$
$\displaystyle\leq\eta_{1}+\eta_{2}.$

∎

The following lemma shows that small total variation distance implies good
coupling approximation.

######  Lemma C.2.

If the total variation distance $d_{TV}(X,Y)\leq\eta$, then
$|X-Y|\leq_{\eta}0$.

###### Proof.

It is well known that for any two random variables $X,Y$, there exists a
coupling such that $d_{TV}(X,Y)=\Pr[X\neq Y]$. This immediately implies what
we want. ∎

###  C.2 Bounding the error using tail bounds of PRVs

The goal of this section is to bound the error of ComposePRV in terms of the
tail bounds of the underlying PRVs.

######  Theorem C.3.

Let $Y_{1},Y_{2},\dots,Y_{k}$ be PRVs and let $\widetilde{Y}$ be the
approximation produced by the ComposePRV algorithm (Algorithm 1) for
$Y=\sum_{i=1}^{k}Y_{i}$ with truncation parameter $L$ and mesh size

$h=\frac{\varepsilon_{\mathrm{error}}}{\sqrt{\frac{k}{2}\log\frac{2}{\eta_{0}}}}.$

Then

$\delta_{\widetilde{Y}}(\varepsilon+\varepsilon_{\mathrm{error}})-\delta_{\mathrm{error}}\leq\delta_{Y}(\varepsilon)\leq\delta_{\widetilde{Y}}(\varepsilon-\varepsilon_{\mathrm{error}})+\delta_{\mathrm{error}}$

where

$\displaystyle\delta_{\mathrm{error}}$
$\displaystyle=\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq
L\right]+\sum_{i=1}^{k}\Pr[|Y_{i}|\geq L]+\eta_{0}$
$\displaystyle\leq\Pr\left[\left|\sum_{i=1}^{k}Y_{i}\right|\geq
L-\varepsilon_{\mathrm{error}}\right]+2\sum_{i=1}^{k}\Pr[|Y_{i}|\geq
L]+2\eta_{0}.$

######  Remark C.4.

We can directly bound
$\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq L\right]$ using
moment generating functions as

$\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq
L\right]\leq\inf_{\lambda>0}\frac{\prod_{i=1}^{k}\operatorname{\mathbb{E}}[\exp(\lambda\widetilde{Y}_{i})]}{e^{\lambda
L}}+\inf_{\lambda>0}\frac{\prod_{i=1}^{k}\operatorname{\mathbb{E}}[\exp(-\lambda\widetilde{Y}_{i})]}{e^{\lambda
L}}.$

Sometimes, if we already have good upper bound for
$\Pr\left[\left|\sum_{i}Y_{i}\right|\geq L\right]$, then the second bound on
$\delta_{\mathrm{error}}$ is useful.

The following key lemma shows that the DiscretizePRV algorithm (Algorithm 2)
produces a good coupling approximation to the PRV and preserves the mean.

######  Lemma C.5.

Given a PRV $Y$, let $Y^{L}=Y\big{|}_{|Y|\leq L}$ be its truncation. The
approximation $\widetilde{Y}$ returned by DiscretizePRV satisfies
$\operatorname{\mathbb{E}}[\widetilde{Y}]=\operatorname{\mathbb{E}}[Y^{L}]$
and $|Y^{L}-(\widetilde{Y}-\mu)|\leq_{0}\frac{h}{2}$ for some $\mu$ where $h$
is the mesh size. We also have $|Y^{L}-Y|\leq_{\eta}0$ where $\eta=\Pr[|Y|\geq
L].$

###### Proof.

Since $d_{TV}(Y,Y^{L})\leq\Pr[|Y|\geq L]=\eta$, by Lemma C.2,
$|Y-Y^{L}|\leq_{\eta}0$. It is clear that by construction
$\operatorname{\mathbb{E}}[\widetilde{Y}]=\operatorname{\mathbb{E}}[Y^{L}],$

$\operatorname{\mathbb{E}}[\widetilde{Y}]=\mu+\sum_{i=-n}^{n}ih\cdot
q_{i}=\left(\operatorname{\mathbb{E}}[Y^{L}]-\sum_{i=-n}^{n}ih\cdot
q_{i}\right)+\sum_{i=-n}^{n}ih\cdot q_{i}=\operatorname{\mathbb{E}}[Y^{L}].$

We will now construct the coupling between $(Y^{L},\widetilde{Y})$ such that
$|Y^{L}-(\widetilde{Y}-\mu)|\leq\frac{h}{2}$. The coupling is as follows:
First sample $y\sim Y^{L}$. Suppose $y\in(ih-\frac{h}{2},ih+\frac{h}{2}]$ for
some integer $i$ such that $-n\leq i\leq n$, then return
$\widetilde{y}=\mu+ih$. Clearly, the distribution of $\widetilde{y}$ matches
with $\widetilde{Y}$ and $|y-(\widetilde{y}-\mu)|=|y-ih|\leq\frac{h}{2}$.

∎

Since our error bound on $\widetilde{Y}$ is slightly different from the
assumption in Lemma 5.3, we need the following generalization using the same
proof.

######  Lemma C.6.

Suppose $Y_{1},Y_{2},\dots,Y_{k}$ and
$\widetilde{Y}_{1},\widetilde{Y}_{2},\dots,\widetilde{Y}_{k}$ are two
collections of independent random variables such that
$|Y_{i}-(\widetilde{Y}_{i}-\mu_{i})|\leq_{0}h$ for some $\mu_{i}$ and
$\operatorname{\mathbb{E}}[Y_{i}]=\operatorname{\mathbb{E}}[\widetilde{Y}_{i}]$
for all $i$, then

$\left|\sum_{i=1}^{k}Y_{i}-\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\leq_{\eta}h\sqrt{2k\log{\frac{2}{\eta}}}.$

In the algorithm, we only calculate the distribution of $Y_{1}\oplus
Y_{2}\oplus\dots\oplus Y_{k}$ instead of $Y_{1}+Y_{2}+\dots+Y_{k}$. The
following simple lemma shows that this is still a good approximation as long
as $\sum_{i}Y_{i}$ stays within $[-L,L]$ with high probability.

######  Lemma C.7.

Let $Y_{1},Y_{2},\dots,Y_{k}$ be random variables supported on $(-L,L]$. Then

$\left|\sum_{i=1}^{k}Y_{i}-\left(Y_{1}\oplus_{L}Y_{2}\oplus_{L}\cdots\oplus_{L}Y_{k}\right)\right|\leq_{\eta}0$

where

$\eta=\Pr\left[\left|\sum_{i=1}^{k}Y_{i}\right|\geq L\right].$

###### Proof.

$\Pr\left[\sum_{i=1}^{k}Y_{i}\neq\left(Y_{1}\oplus_{L}Y_{2}\oplus_{L}\cdots\oplus_{L}Y_{k}\right)\right]\leq\Pr\left[\left|\sum_{i=1}^{k}Y_{i}\right|\geq
L\right].$

This clearly implies what we want. ∎

Combining all the above lemmas, we get the following corollary.

######  Corollary C.8.

Let $Y_{1},Y_{2},\dots,Y_{k}$ be random variables supported on and let
$\widetilde{Y}_{i}$ be the discretization of $Y_{i}$ produced by DiscretizePRV
algorithm with mesh size
$h=\frac{h_{0}}{\sqrt{\frac{k}{2}\log\frac{2}{\eta_{0}}}}$ and truncation
parameter $L$. Then

$\left|(Y_{1}+Y_{2}+\dots+Y_{k})-(\widetilde{Y}_{1}\oplus\widetilde{Y}_{2}\oplus\dots\oplus\widetilde{Y}_{k})\right|\leq_{\eta}h_{0}$

where

$\eta=\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq
L\right]+\sum_{i=1}^{k}\Pr[|Y_{i}|\geq L]+\eta_{0}.$

Furthermore, we can bound

$\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq
L\right]\leq\Pr\left[\left|\sum_{i=1}^{k}Y_{i}\right|\geq
L-h_{0}\right]+\sum_{i=1}^{k}\Pr[|Y_{i}|\geq L]+\eta_{0}.$

###### Proof.

Let $Y^{L}\equiv Y_{i}\big{|}_{|Y_{i}|\leq L}$ be the truncation of $Y_{i}.$
By Lemma C.5, $|Y_{i}^{L}-(\widetilde{Y}_{i}-\mu_{i})|\leq_{0}\frac{h}{2}$ for
some $\mu_{i}$ and $|Y_{i}^{L}-Y_{i}|_{\xi_{i}}\leq 0$ where
$\xi_{i}=\Pr[|Y_{i}|\geq L].$ Now applying the triangle inequality for
coupling approximations (Lemma C.1), we have

$\left|\sum_{i}Y_{i}-\sum_{i}Y_{i}^{L}\right|\leq_{\eta_{1}}0$

where $\eta_{1}=\sum_{i}\xi_{i}=\sum_{i}\Pr[|Y_{i}|\geq L].$ By Lemma C.6, we
have

$\left|\sum_{i}Y_{i}^{L}-\sum_{i}\widetilde{Y}_{i}\right|\leq_{\eta_{0}}\frac{h}{2}\cdot\sqrt{2k\log\frac{2}{\eta_{0}}}=h\sqrt{\frac{k}{2}\log\frac{2}{\eta_{0}}}=h_{0}.$

By Lemma C.7,

$\left|\sum_{i=1}^{k}\widetilde{Y}_{i}-\left(\widetilde{Y}_{1}\oplus_{L}\widetilde{Y}_{2}\oplus_{L}\cdots\oplus_{L}\widetilde{Y}_{k}\right)\right|\leq_{\eta_{2}}0$

where $\eta_{2}=\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq
L\right].$ Finally applying triangle inequality (Lemma C.1) once again, we
get:

$\left|(Y_{1}+Y_{2}+\dots+Y_{k})-(\widetilde{Y}_{1}\oplus_{L}\widetilde{Y}_{2}\oplus_{L}\dots\oplus_{L}\widetilde{Y}_{k})\right|\leq_{\eta}h_{0}$

where $\eta=\eta_{0}+\eta_{1}+\eta_{2}$. We can bound
$\Pr\left[\left|\sum_{i=1}^{k}\widetilde{Y}_{i}\right|\geq L\right]$ as:

$\displaystyle\Pr\left[\left|\sum_{i}\widetilde{Y}_{i}\right|\geq L\right]$
$\displaystyle=\Pr\left[\left|\sum_{i}(\widetilde{Y}_{i}-Y_{i}^{L})+\sum_{i}(Y_{i}^{L}-Y_{i})+\sum_{i}Y_{i}\right|\geq
L\right]$
$\displaystyle\leq\Pr\left[\left|\sum_{i}(\widetilde{Y}_{i}-Y_{i}^{L})\right|+\left|\sum_{i}(Y_{i}^{L}-Y_{i})\right|+\left|\sum_{i}Y_{i}\right|\geq
h_{0}+0+L-h_{0}\right]$
$\displaystyle\leq\Pr\left[\left|\sum_{i}(\widetilde{Y}_{i}-Y_{i}^{L})\right|>h_{0}\right]+\Pr\left[\left|\sum_{i}(Y_{i}^{L}-Y_{i})\right|>0\right]+\Pr\left[\left|\sum_{i}Y_{i}\right|\geq
L-h_{0}\right]$
$\displaystyle\leq\eta_{0}+\eta_{1}+\Pr\left[\left|\sum_{i=1}^{k}Y_{i}\right|\geq
L-h_{0}\right].$

∎

###### Proof of Theorem C.3.

Combining Corollary C.8 (with $h_{0}=\varepsilon_{\mathrm{error}}$) and Lemma
5.2, we have Theorem C.3. ∎

###  C.3 Tail Bound for PRVs

To finish the proof of our main theorem (Theorem 5.5, we need a tail bound on
PRVs in terms of their privacy curves. First, we need a lemma relating the
PRVs of a privacy curve $\delta(P||Q)$ with the PRVs of $\delta(Q||P)$.

######  Proposition C.9.

Let $(X,Y)$ be the PRVs for a privacy curve $\delta(P||Q)$. Then the PRVs for
the privacy curve $\delta(Q||P)$ are $(-Y,-X).$

###### Proof.

Let $(\widetilde{X},\widetilde{Y})$ be the PRVs for $\delta(Q||P)$. We know
that $\delta(P||Q)=\delta(X||Y)$. So $\delta(Q||P)=\delta(Y||X).$ Then by
Theorem 3.2,

$\displaystyle\widetilde{X}$
$\displaystyle=\log\left(\frac{X(t)}{Y(t)}\right)\text{ where }t\sim Y$
$\displaystyle=\log\left(e^{-t}\right)\text{ where }t\sim Y$
$\displaystyle=-Y.$ $\displaystyle\widetilde{Y}$
$\displaystyle=\log\left(\frac{X(t)}{Y(t)}\right)\text{ where }t\sim X$
$\displaystyle=\log\left(e^{-t}\right)\text{ where }t\sim X$
$\displaystyle=-X.$

∎

Now, we show our tail bound, which shows the PRVs $(X,Y)$ for a
$(\varepsilon,\delta)$-DP algorithm satisfies roughly that
$\Pr(|Y|\geq\varepsilon+2)\leq 2\delta$.

See 5.4

###### Proof.

We have $\delta(X||Y)\leq f_{\varepsilon,\delta}$ and $\delta(Y||X)\leq
f_{\varepsilon,\delta}$ where $f_{\varepsilon,\delta}$ is the privacy curve of
a $(\varepsilon,\delta)$-DP algorithm. By Theorem 3.2, we have

$\displaystyle\delta$
$\displaystyle\geq\int_{0}^{\infty}\Pr[Y\geq\varepsilon+s]e^{-s}ds$
$\displaystyle\geq\int_{0}^{t}\Pr[Y\geq\varepsilon+s]e^{-s}ds$
$\displaystyle\geq\Pr[Y\geq\varepsilon+t]\int_{0}^{t}e^{-s}ds$
$\displaystyle\geq\Pr[Y\geq\varepsilon+t](1-e^{-t}).$

By Proposition C.9, the PRVs for $\delta(Y||X)$ are $(-Y,-X)$. Therefore by a
similar argument, we have

$\Pr[X\leq-\varepsilon-t]=\Pr[-X\geq\varepsilon+t]\leq\frac{\delta}{1-e^{-t}}.$

Finally, note that $Y(s)=e^{s}X(s)$ for all $s\in\mathbb{R}$ and
$Y(-\infty)=0$ by the definition of PRVs. Therefore

$\Pr[Y\leq-\varepsilon-t]\leq e^{-\varepsilon-t}\Pr[X\leq-\varepsilon-t].$

Therefore we have:

$\displaystyle\Pr[|Y|\geq\varepsilon+t]$
$\displaystyle=\Pr[Y\geq\varepsilon+t]+\Pr[Y\leq-\varepsilon-t]$
$\displaystyle\leq\Pr[Y\geq\varepsilon+t]+e^{-\varepsilon-t}\Pr[X\leq-\varepsilon-t]$
$\displaystyle\leq\left(1+e^{-\varepsilon-t}\right)\frac{\delta}{1-e^{-t}}.$

∎

###  C.4 Proof of Theorem 5.5

Now, we can prove our main theorem. See 5.5

###### Proof.

By Lemma 5.4,

$\displaystyle\Pr[|Y_{i}|\geq L]$ $\displaystyle=\Pr[|Y_{i}|\geq L-2+2]$
$\displaystyle\leq\delta_{Y_{i}}(L-2)\cdot\frac{1+e^{-L}}{1-e^{-2}}$
$\displaystyle\leq\delta_{Y_{i}}(L-2)\cdot\frac{1+e^{-2}}{1-e^{-2}}$
$\displaystyle\leq\delta_{Y_{i}}(L-2)\cdot\frac{4}{3}.$

Therefore we have

$\sum_{i=1}^{k}\Pr[|Y_{i}|\geq
L]\leq\frac{4}{3}\sum_{i=1}^{k}\delta_{Y_{i}}(L-2)\leq\frac{4}{3}\cdot\frac{\delta_{\mathrm{error}}}{8}=\frac{\delta_{\mathrm{error}}}{6}.$

Similarly

$\Pr[|Y|\geq
L-\varepsilon_{\mathrm{error}}]\leq\frac{4}{3}\delta_{Y}(L-2-\varepsilon_{\mathrm{error}})\leq\frac{\delta_{\mathrm{error}}}{3}.$

Therefore by Theorem C.3, setting
$\eta_{0}=\frac{\delta_{\mathrm{error}}}{6}$, we get the desired result.

For the runtime, we note that the bottleneck of our algorithm is to compute
the convolution, which can be done using FFT. In total, we need to compute
$b+1$ many FFT for $b$ distinct algorithms, one for each for computing the
Fourier transform and one of computing the inverse Fourier transform. Since
the length of the array for the FFT is bounded by $O(L/h)$, this costs
$O(bL/h\log(L/h))$ in total.

The step
$\delta_{\widetilde{Y}}(\varepsilon)=\operatorname{\mathbb{E}}_{\widetilde{Y}}\left[\left(1-e^{\varepsilon-\widetilde{Y}}\right)_{+}\right]$
can be computed in linear time by first computing the CDF of $\widetilde{Y}$
and the prefix sum
$\operatorname{\mathbb{E}}_{\widetilde{Y}\leq\alpha}\left[e^{-\widetilde{Y}}\right]$
for all $\alpha$. ∎