arXiv CSopen access

A Sieve-Accelerated Quadrature Method for Exact Privacy Accounting in the 2020 U.S. Decennial Census

A Sieve-Accelerated Quadrature Method for Exact Privacy Accounting in the 2020 U.S. Decennial Census Buxin Su∗

Weijie Su†

Chendi Wang‡

arXiv:2606.29835v1 [cs.CR] 29 Jun 2026

June 30, 2026

Abstract In 2020, the U.S. Census Bureau adopted differential privacy for the Decennial Census by injecting integer-valued Gaussian noise into published census tabulations. Exactly evaluating the privacy guarantees of these data releases would enable the Bureau to determine the absolute minimum noise required to satisfy a given privacy budget, preventing the injection of unnecessary excess noise and thereby substantially enhancing the statistical utility of the data for downstream applications such as federal funding allocation and political redistricting. In this paper, we introduce a computationally efficient and mathematically rigorous quadrature method to evaluate the exact privacy profile of practical, large-scale census releases under the composition of heterogeneous discrete Gaussian mechanisms. Mathematically, this problem reduces to evaluating the tail probabilities of high-dimensional convolutions of integer-valued random variables sampled from heterogeneous discrete Gaussian distributions under exceptionally stringent numerical error tolerances (e.g., 10−35 ). By recasting the exact privacy accounting as a numerical integration problem via the discrete Fourier transform, we explicitly exploit the exponential convergence of the trapezoidal rule for complex analytic, periodic characteristic functions. Furthermore, to overcome the computational bottleneck of evaluating highly oscillatory integrands in high dimensions, we develop a sieve algorithm that identifies and prunes negligible quadrature nodes, accelerating the computation by three orders of magnitude. Taken together, these numerical innovations enable the first exact, assumption-free privacy accounting for the 2020 Census Demographic and Housing Characteristics File, achieving a 1,824-fold speedup over prior methods while maintaining census-mandated error tolerances.

1

Introduction

The U.S. Decennial Census is among the most influential statistical programs in the world, shaping federal funding allocation [31, 59], political redistricting [14, 34], congressional apportionment [25, 54], and a wide range of labor and economic research [7, 55]. Decennial Census products, such as the 2020 Census Demographic and Housing Characteristics (DHC) File, contain highly ∗

University of Pennsylvania; Email: [email protected]. University of Pennsylvania; Email: [email protected]. ‡ Xiamen University; Email: [email protected]. § Authors are listed in alphabetical order. †

1

sensitive demographic information [29]. Consequently, the direct release of exact census counts creates severe privacy risks through reconstruction and re-identification attacks [1, 18, 30, 39]. To rigorously mitigate these concerns, the Census Bureau adopted differential privacy [23, 24] as the official confidentiality protection framework for the 2020 Census. Differential privacy provides a mathematical guarantee that changing any single individual’s record has only a strictly bounded effect on the released output. For implementation, the Census Bureau utilized the TopDown algorithm [3, 4, 15] within its Disclosure Avoidance System (DAS). The DAS injects integer-valued noise, sampled from discrete Gaussian distributions, into a massive collection of counting queries across hierarchical geographic levels, followed by post-processing to ensure data consistency and non-negativity. The Census Bureau currently provides formal privacy guarantees for the 2020 Census DHC File using zero-concentrated differential privacy (zCDP) [12, 44]. However, the privacy guarantees derived from zCDP are upper bounds on the (unknown) exact privacy loss. The gap between the bounds and the exact privacy loss leads to a suboptimal privacy–utility trade-off, resulting in the injection of unnecessarily high levels of noise into the census counts for a given privacy budget. This degrades the statistical accuracy of the data and, in turn, impacts downstream applications such as redistricting and social science research [51, 33, 35, 5, 9]. In addition, zCDP parameters lack a straightforward statistical interpretation, often requiring technical expertise to inform public policy discussions [36]. A crucial fact that is frequently overlooked in the literature is that evaluating the exact privacy loss for the Census DHC File can be cast as a well-defined mathematical problem. Indeed, recent work by [48] demonstrates that the total privacy loss can be expressed exactly as a composition of trade-off curves associated with the census queries across hierarchical geographic levels under the f -differential privacy (f -DP) framework [19]. This equivalently reduces the exact accounting problem to computing the tail probabilities of high-dimensional convolutions of heterogeneous discrete Gaussian distributions, that is, discretized counterparts of the continuous Gaussian distribution supported solely on the integers. While the exact privacy level is conceptually straightforward to express from a mathematical perspective, numerically evaluating it via quadrature presents significant computational challenges. To achieve numerical tractability, prior methods [48] relied on theoretical simplifications of the dataset’s sensitivity by adopting the add-remove model instead of the replacement model for neighboring datasets [16] and assuming homogeneity of the injected integer noise under composition. While this effort resolved an open question posed by the Census Bureau [36], demonstrating that noise variances could be reduced while maintaining (nearly) the same level of privacy protection, the computational cost proved prohibitive. Specifically, evaluating a single privacy budget allocation required approximately 1,000 CPU hours.1 Even worse, the 2020 Census DHC File involves 946 distinct pairs of privacy budget allocations. Requiring 1,000 CPU hours for each of the 946 evaluations renders comprehensive, full-scale accounting computationally impractical. From a mathematical viewpoint, the difficulty of this numerical evaluation stems from several intertwined issues. Unlike continuous Gaussian mechanisms, where convolutions remain Gaussian and admit closed-form expressions, convolutions of discrete Gaussian random variables yield proba1

All CPU times in this paper are measured on an AWS EC2 c5.metal instance.

2

bility mass functions that oscillate rapidly around their continuous approximations. This oscillatory behavior is compounded by high-dimensional, multi-fold compositions (e.g., 80-fold compositions along geographic paths in the DHC File).2 This renders standard continuous analytic composition approaches [32, 68] inapplicable. Moreover, the discrete Gaussian mechanism used in the 2020 Census is heterogeneous, because the injected noise variances are not identical across queries and geographic levels. This heterogeneity makes accounting substantially more difficult, as it destroys the algebraic structure that would otherwise permit closed-form or low-complexity approximations. Furthermore, the design of the 2020 Census DHC File demands exceptionally stringent numerical error tolerances (e.g., bounded by 10−25 or 10−35 ). Such extreme precision requirements render standard numerical approaches, including fast Fourier transform (FFT)-based quadrature methods [37, 27], infeasible due to accumulated truncation errors.

1.1

Contributions

In this paper, we develop a highly efficient, mathematically rigorous quadrature-based numerical integration method for the exact privacy accounting of the 2020 Census DHC File that explicitly accommodates fully practical settings. Under identical problem settings, including the same numerical precision requirements, our proposed method requires only 0.5 CPU hours to evaluate a single privacy budget allocation, compared to the roughly 1,000 CPU hours required by the prior state-of-the-art [48]. This acceleration by over three orders of magnitude enables the first exact, assumption-free privacy accounting for all 946 pairs of privacy budget allocations in the 2020 Census DHC File. Consequently, we provide a comprehensively improved privacy–utility trade-off compared to the Bureau’s zCDP-based approach [3], demonstrating that exact accounting permits noise variance reductions of 15.08% to 24.82% without sacrificing privacy guarantees. This addresses the open problem posed by the Census Bureau in a substantially more comprehensive manner [36]. Our methodological and computational advances are enabled by three key numerical innovations, each focused on resolving the complexities inherent in composing heterogeneous mechanisms: Exponentially convergent quadrature via the trapezoidal rule. By reformulating the discrete privacy accounting problem as a numerical integration task via the discrete Fourier transform, we explicitly leverage the exponential convergence of the trapezoidal rule for complex analytic, periodic integrands. This approach achieves the required precision accuracy using exponentially fewer quadrature nodes. Notably, this exponential convergence rate is not automatic; we establish it only after a careful, rigorous complex-analytic analysis of the discrete Gaussian characteristic function, explicitly bounding its growth within a strip in the complex plane. Sub-Gaussian tail bounds via lattice mapping for heterogeneous weights. To exploit the exponential convergence of the trapezoidal rule, we must first truncate the infinite sum involved in the discrete Fourier transform. An obstacle in doing so arises from the heterogeneity of discrete Gaussian mechanisms under composition, where the weighted sum of the integer-valued noise no 2 When the compositional dimension is extremely large, continuous approximations suffice; when the dimension is very small, naive numerical approximations remain viable. An intermediate dimension such as 80 falls into a difficult regime where neither approach succeeds without substantial error.

3

longer naturally resides on a simple lattice support. We overcome this by recognizing that the privacy budgets allocated by the Census Bureau are rational numbers. By deriving a unified scaling factor that maps the heterogeneous weighted sum onto a common integer lattice via Bézout’s identity, we preserve the discrete probability mass function without resorting to lossy binning or continuous interpolation. This enables us to establish precise sub-Gaussian concentration inequalities for heterogeneous discrete Gaussian mechanisms, with bounds that guarantee the tail mass falls below the prescribed error margins. Accelerated integration via a sieve algorithm. Even with exponential convergence, achieving a 10−35 tolerance in high-dimensional settings demands the evaluation of millions of quadrature nodes. Inspired by sieve methods in number theory, we develop a novel algorithmic procedure that exploits the periodic peak structure of the characteristic functions to rigorously identify and discard nodes where the highly oscillatory integrand is guaranteed to be negligible. This sieve drastically reduces the number of required integrand evaluations from over ∼ 106 nodes down to just 203. As an aside, we also provide a theoretical explanation, grounded in the central limit theorem, for why exact privacy accounting for compositions of the discrete Gaussian mechanism strictly diverges from that of the continuous Gaussian mechanism. Specifically, using an Edgeworth expansion, we prove that a central-limit-theorem-type approximation for the convolution of discrete Gaussian distributions involves a non-negligible, higher-order oscillatory discrepancy term. We show that this term is explicitly bounded by O(m−3/2 ), where m denotes the number of folds in the composition. Our code is available at https://github.com/BuxinSu/Exact-Privacy-Accounting-for-2020-U. S.-Census.git.

1.2

Organization of the Paper

The remainder of the paper is organized as follows. Section 2 reviews preliminaries on differential privacy and the 2020 U.S. Census, with particular emphasis on the 2020 Census DHC File. Section 3 reformulates the privacy accounting problem as a numerical computation problem. Section 4 presents our new privacy accounting method, which enables exact computation of the (ε, δ)-DP and f -DP curves for the 2020 Census DHC File; the corresponding empirical results are reported in Section 5. Finally, Section 6 provides theoretical insight into why the discrete Gaussian mechanism differs from the continuous one from the perspective of the central limit theorem.

2

Preliminaries

In this section, we review the differential privacy frameworks used in our analysis and in the 2020 Census DHC File. We then outline how differential privacy is implemented in the Census Bureau’s DAS via the discrete Gaussian mechanism and hierarchical privacy budget allocation.

2.1

Differential Privacy

We first introduce three different variants of differential privacy. The most well-known and standard definition is (ε, δ)-DP [23, 24]. For ease of computing compositions, the method currently adopted by the Census Bureau uses ρ-zCDP [12, 22, 44, 13, 3] to determine the privacy level. Our improved 4

approach is largely based on f -DP [19], which is defined from a statistical perspective and provides a natural interpretation of the trade-off between privacy level and the power of re-identification inferences [19, 36]. (ε, δ)-DP. The (ε, δ)-DP notion, introduced by [23, 24], is the first formal definition of differential f is said to satisfy (ε, δ)-DP for ε ≥ 0 and 0 ≤ δ ≤ 1 if, for any privacy. A randomized mechanism M pair of neighboring datasets D and D′ —which differ in a single individual record—and any event S, the following holds: f(D) ∈ S) ≤ eε · P(M f(D′ ) ∈ S) + δ. P(M (2.1) f is regarded as more private when ε and δ are small. In particular, when ε = δ = 0, the M f(D) and M f(D′ ) are identical, meaning that the mechanism achieves perfect distributions of M privacy. ρ-zCDP. The privacy budget for the 2020 Census DHC File is measured using zero-Concentrated Differential Privacy [12], which is based on Rényi divergence. For two distributions P and Q with probability density functions p and q, respectively, the Rényi divergence of order α > 1 is defined R 1 log p(x)α q(x)1−α dx. The quantities R1 (P ∥Q) and R∞ (P ∥Q) are defined as as Rα (P ∥Q) = α−1 the limits of Rα (P ∥Q) as α → 1 and α → ∞, respectively. Based on Rényi divergence, one obtains the definition of zCDP, where the divergence between two random variables is understood as the divergence between their corresponding distributions. f is said to satisfy ρ-zCDP if Definition 2.1 (zCDP, [12]). A randomized mechanism M f(D)∥M f(D′ )) ≤ ρα, Rα (M

for all α > 1,

for any neighboring datasets D and D′ . f(D) and M f(D′ ) are closer to each other. The privacy When ρ is small, the distributions of M budget allocated by the Census Bureau [57] is ρ = 4.9622 for the 2020 Census DHC File. f -DP. In this paper, we employ the recently developed f -DP framework [19], which has been shown to be well suited for privacy analysis under composition [10, 63, 49, 41]. To define f -DP, consider formulating the problem of distinguishing between a pair of neighboring datasets D and D′ as a hypothesis testing problem: H0 : the true dataset is D

H1 : the true dataset is D′ .

versus

Let 0 ≤ ϕ ≤ 1 be any (possibly randomized) rejection rule, and denote by αϕ = EH0 [ϕ] and f(D), M f(D′ )) : βϕ = 1−EH1 [ϕ] the type I and type II errors, respectively. The trade-off function T (M ′ [0, 1] → [0, 1] between D and D is defined as f(D), M f(D′ ))(α) = inf {βϕ : αϕ ≤ α} T (M ϕ

5

for any 0 ≤ α ≤ 1 [19].3 That is, it characterizes the minimal achievable type II error for a given f satisfies f -DP if level of type I error. We say that a mechanism M f(D), M f(D′ ))(α) ≥ f (α) T (M for any neighboring D and D′ and any α ∈ [0, 1]. Typically, we require f : [0, 1] → [0, 1] to be a valid trade-off function for some pair of distributions. This holds if and only if f is continuous, convex, non-increasing, and satisfies f (α) ≤ 1 − α. A larger trade-off function f indicates that it is more difficult to distinguish between H0 and H1 , and hence the mechanism provides stronger privacy. It is worth noting that f -DP not only enables exact privacy analysis, but also offers a more meaningful semantic interpretation of differential privacy [36, 26] compared with (ε, δ)-DP. Mathematically, an f -DP guarantee is equivalent to an infinite collection of guarantees specified by the (ε, δ) privacy curve for all ε > 0 [19]. However, in practice, computing both the exact trade-off function f and the corresponding (ε, δ)-DP curve is f may be complex, especially under composition. In generally challenging, as the distribution of M this paper, we present an efficient and practical method for computing the exact f -DP curve and the corresponding (ε, δ)-DP curve. Composition. Composition describes the cumulative privacy loss incurred when multiple mechfi }k satisfy (εi , δi )-DP anisms are applied to the same dataset. For (ε, δ)-DP, if mechanisms {M i=1 Pk Pk individually, then their sequential composition satisfies ( i=1 εi , i=1 δi )-DP. For zCDP, Rényi divergences are additive under composition, so the total privacy parameter ρ is the sum of the indiP vidual RDP parameters ρi of each mechanism, i.e., ρ = ki=1 ρi . For f -DP, composition corresponds to the tensor product of privacy trade-off functions. Precisely, we have fi (D)}ki=1 , {M fi (D′ )}ki=1 ) = T T ({M

k Y i=1

Pi ,

k Y

! Qi

,

i=1

f(D) and M f(D′ ), correspondingly. This yields where Pi and Qi are the underlying distributions of M an exact characterization of cumulative privacy through hypothesis-testing trade-offs. Parallel composition captures the privacy loss when multiple mechanisms are applied to disjoint subsets of the data. Suppose the dataset D is partitioned into k disjoint blocks D(1) , . . . , D(k) , and each fi accesses only D(i) . The composition of the sequence M f1 (D(1) ), . . . , M fk (D(k) ) is mechanism M characterized in Proposition 3.1. Sensitivity. as

For a deterministic query M taking values in Rd , the ℓ2 -sensitivity of M is defined

∆2 = sup

D,D′

n

o M (D) − M (D′ ) ℓ2 ,

f(D) and M f(D′ ), correspondingly. Formally, the trade-off Let P and Q denote the probability distributions of M ′ f(D), M f(D )) should be defined through P and Q, thereby being expressed as T (P, Q). function T (M 3

6

where ∥ · ∥ℓ2 denotes the ℓ2 norm of a vector, and the supremum is taken over all neighboring datasets D and D′ that differ in at most one data record. In the 2020 Census DHC File, a bounded (replacement) sensitivity is used for each marginal query. This reflects the difference between two neighboring datasets under a single-record replacement. The DAS [58, 3] considers the sensitivity of coarsened counting queries with binary categories (e.g., “18 and older” vs. “17 and younger”), as detailed in Section 2.2. Replacing one individual with another whose attributes fall into different categories can change each of the counts (“18 and older” and “17 and younger”) by 1.

2.2

The 2020 Census DHC File

Data structure. The 2020 Census DHC File is a major official Census product released by the U.S. Census Bureau [53]. It is a collection of integer-valued counting statistics, with each corresponding to an item in {geographic units} × {marginal queries} [16]. Geographically, the 2020 Census DHC File partitions the United States into eight geographic levels: US, State, County, PRIM, Tract Subset, Tract Subset Group, Block Group, Block.

(2.2)

Each node, or unit, within a geographic level is referred to as a geographic unit. For example, Pennsylvania is a geographic unit at the State level. The smallest geographic unit in the 2020 Census DHC File is the Block, which typically corresponds to a very small area, such as a few street blocks in a city. All other geographic units are defined as aggregations of blocks. Therefore, for any fixed block, there is at most one geographic unit at each geographic level (2.2) that contains this block. Using this property, we define the directed path Pb of a block b as the collection of all unique geographic units containing this block across the geographic levels (2.2). For example, the directed path of a block near the University of Pennsylvania would begin with US, Pennsylvania, Philadelphia, · · ·. For each geographic unit, the 2020 Census DHC File applies the same set of ten marginal queries, including topics on age, sex, race, Hispanic or Latino origin, household type, family type, relationship to householder, group quarters population, housing occupancy, and housing tenure. The collection of integer-valued counting statistics for all items in {geographic units} × {marginal queries} constitutes the unprivatized data. An example of the counting statistics for items in {geographic units} × {marginal queries} is the number of Hispanic people aged 18 or older in Pennsylvania. The final step in producing the 2020 Census DHC File is to add noise to each entry in the unprivatized data to protect privacy using the DAS [3], which we introduce later in this section. Discrete Gaussian Mechanism. The noise distribution used in the DAS is the discrete Gaussian distribution [43, 13], a distribution on the integers Z whose probability mass function is proportional to the Gaussian density. The discrete Gaussian distribution with noise parameter σ 2 ,

7

denoted by NZ (0, σ 2 ), has probability mass function 2

2

e−x /2σ pσ (x) = P −i2 /2σ 2 i∈Z e

(2.3)

 T for any x ∈ Z. For any parameter vector σ 2 = σ12 , · · · , σn2 , we define the multivariate discrete  T Gaussian distribution with parameters σ 2 as N(0, σ 2 ) = NZ (0, σ12 ), · · · , NZ (0, σn2 ) . Let X denote the sample space, and let D ⊂ X be a dataset. Consider a d-dimensional deterministic query M : X → Zd that takes integer values. The discrete Gaussian mechanism adds noise to M(D) as f M(D) = M(D) + NZ (0, σ 2 ).

(2.4)

where NZ (0, σ 2 ) denotes the d-dimensional multivariate discrete Gaussian distribution. Disclosure Avoidance System. Recall that the unprivatized data form a collection of integervalued counting statistics for each item in {geographic units} × {marginal queries}. Mathematically, for any fixed directed path Pb , we denote Mb = {M(g,q)1 , M(g,q)2 , · · · , M(g,q)n }

(2.5)

as the collection of integer-valued counting statistics for each (g, q), where g denotes the unique geographic unit containing block b at each geographic level, q denotes the marginal query, and n = |geographic units containing block b| × |number of marginal query| = 80.

(2.6)

A typical pair (g, q) could be (the state containing b, “Hispanic people aged 18 and older”), in which case M(g,q) denotes the number of people aged 18 and older in the tract containing block b. The DAS protects data privacy by adding independent discrete Gaussian noise with noise parameter σi2 to each M(g,q)i in (2.5). The noise parameter σi2 may vary across different (g, q)i and is determined as follows. Along each directed path Pb , the 2020 Census DHC File allocates a privacy budget ρi to each item in (2.5), such that the sum of all ρi equals a fixed overall budget ρ. Mathematically, along each directed path Pb , the 2020 Census DHC File assigns a sequence of privacy budgets corresponding one-to-one with (2.5): ρb = {ρ1 , ρ2 , · · · , ρn },

(2.7)

P where ρ = i ρi = 4.9622, and the sequence ρb is referred to as the privacy budget allocation along Pb . Note that the allocation ρb may depend on b and differ across directed paths. One example of ρb is shown in Table 1. Among millions of directed paths [3], a total of 43 distinct privacy budget allocations are used in the 2020 Census DHC File. Given the privacy budget ρi for each item in (2.5), the 2020 Census DHC File injects independent

8

discrete Gaussian noise with noise parameter σi2 into each M(g,q)i as f(g,q) = M(g,q) + NZ (0, σi2 ), M i i

σi2 =

1 . ρi

(2.8)

This produces the privatized counting statistics f b = {M f(g,q) , M f(g,q) , · · · , M f(g,q) } M n 1 2

with noise parameters σb2 = [σ12 , σ22 , · · · , σn2 ].

(2.9)

Under the above mechanism, the 2020 Census DHC File satisfies ρ-zCDP, as established in Theorem 2 of [3], due to the composition properties of the discrete Gaussian mechanism [13].4 In the f 0 to denote the privatized counting statistics with σ 2 determined remainder of the paper, we use M i by Table 1 and (2.8). Table 1: One of the ρi allocations used in the 2020 Census DHC File is shown here. Each row represents a marginal query, and each column corresponds to a geographic level. The remaining allocations in (2.7) are provided in Section E of the Supplementary Material. Block 11/10000 11/10000 11/10000 11/10000 11/10000 11/10000 11/10000 11/10000 11/10000 11/10000

2.3

Block Group 43/1000 43/1000 43/1000 43/1000 43/1000 43/1000 43/1000 43/1000 43/1000 43/1000

County 31/1000 31/1000 31/1000 31/1000 31/1000 31/1000 31/1000 31/1000 31/250 31/250

PRIM 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/1250

State 999/10000 999/10000 999/10000 999/10000 999/10000 999/10000 999/10000 999/10000 999/10000 999/2500

Tract Subset 217/2500 217/2500 217/2500 217/2500 217/2500 217/2500 217/2500 217/2500 217/625 217/625

Tract Subset Group 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/5000 239/1250

US 73/10000 73/10000 73/10000 73/10000 73/10000 73/10000 73/10000 73/10000 73/10000 73/2500

Notation

We summarize here the notation used throughout the paper. We denote by P the probability measure, and by Z, N, Q+ , and R the sets of integers, natural numbers, positive rational numbers, and real numbers, respectively. Bold symbols such as X or x denote matrices or vectors, while plain symbols such as x denote scalar quantities. NZ (0, σ 2 ) denotes the univariate discrete Gaussian distribution (2.3) with noise parameter σ 2 . NZ (0, σ 2 ) =  T NZ (0, σ12 ), · · · , NZ (0, σn2 ) denotes the multivariate discrete Gaussian distribution with parameter  T vector σ 2 = σ12 , · · · , σn2 . The 2020 Census DHC File. We consider the 2020 Census DHC File [53]. Let b denote an individual census block, the smallest geographic unit in the 2020 Census DHC File. The 2020 Census DHC File includes eight geographic levels, as listed in (2.2), where at each level there is at most one geographic unit containing block b. We define the directed path Pb as the collection of 4

A post-processing step is applied before releasing the 2020 Census DHC File (see [3]). However, since privacy guarantees are invariant under post-processing [21], we omit the details to avoid confusion.

9

all geographic units that contain block b. For any fixed directed path Pb , the same set of marginal queries is applied to every geographic unit, producing the integer-valued counting statistics Mb in (2.5). The 2020 Census DHC File assigns a privacy budget allocation ρb as in (2.7) along Pb , with the total budget equal to the overall privacy level ρ = 4.9622. Finally, the 2020 Census DHC File injects independent discrete Gaussian noise with noise parameter σi2 into each M(g,q)i as in (2.8), f b in (2.9). producing the privatized counting statistics M f for its privatized (noise-added) version Privacy accounting. For a mechanism M, we write M and M for the corresponding unprivatized mechanism. We use ρ to denote the overall privacy budget under the zCDP framework. The pair (ε, δ) denotes the privacy parameters in (ε, δ)-DP. The pair (α, β) denotes the type I and type II errors on the f -DP (trade-off) curve.

3

Privacy Accounting Formulation

In this section, we formulate privacy accounting for the 2020 Census DHC File as a statistical hypothesis testing and numerical computation problem. We first express the overall privacy guarantee through pairwise compositions of neighboring datasets, showing that exact privacy levels reduce to evaluating hypothesis-testing trade-off functions under composed discrete Gaussian mechanisms. We then transform both the (ε, δ)-DP and f -DP privacy profiles into a unified numerical task—accurately computing tail probabilities of weighted convolutions of discrete Gaussian distributions—which serves as the foundation for our methodology.

3.1

From preliminaries to privacy accounting

We now reformulate the preliminaries as a privacy accounting problem. Consider two neighboring datasets D and D′ obtained by replacing a record x on directed path Pk with another record x′ on directed path Pl . For example, this corresponds to replacing an individual living in (US, Pennsylvania, Philadelphia, · · · , the University of Pennsylvania) with another individual living in f k and M f l denote the privatized (US, Massachusetts, Cambridge, · · · , Harvard University). Let M counting statistics along directed path Pk and Pl , respectively. Let m = 2n = 160 denote the fk , M f l ], with n as defined in (2.6). number of dimensions of the composed mechanism [M   NZ (0, σ12 ) " #     f k (D) M Mk (D)  NZ (0, σ22 )  = + ,  f l (D) Ml (D) ···   M 2 ) NZ (0, σm 2 ] = [σ 2 , σ 2 ] is the concatenation of the two sequences of noise parameters where [σ12 , σ22 , · · · , σm k l f k and M f l , as specified in (2.8) and (2.9). used in M Let µ = 1 denote the sensitivity of each coordinate of Mk and Ml . Then the composed mechfk , M f l ] satisfies fkl -DP, where fkl is the trade-off function associated with the following anism [M

10

hypothesis testing problem: 

 NZ (0, σ12 )    N (0, σ22 )  H0 :  Z ,  ··· ,  2 ) NZ (0, σm



 NZ (µ, σ12 )    N (µ, σ22 )  H1 :  Z .  ··· ,  2 ) NZ (µ, σm

(3.1)

The following proposition characterizes the overall privacy level of the 2020 Census DHC File. Its proof follows directly from the fact that the overall privacy guarantee of the 2020 Census DHC File is given by the parallel composition (cf., [42] and [47]) over all pairs of composed mechanisms fk , M f l ]. Note that the parallel composition requires taking minimum over all fkl while taking [M maximum over εkl . Proposition 3.1. The overall privacy level of the 2020 Census DHC File is (mink,l fkl )∗∗ where (mink,l fkl )∗∗ is the lower convex envelope [62] of mink,l fkl . Moreover, suppose the composition fk , M f l ] satisfies (εkl , δ)-DP, then overall privacy budget in the 2020 Census of two mechanisms [M DHC File is (maxk,l εkl , δ)-DP. Proposition 3.1 indicates that, in order to obtain the exact privacy guarantees for the 2020 Census DHC File, it suffices to compute the exact privacy accounting quantities fkl and the (εkl , δ) curve for all pairs of directed paths k and l. To this end, we leverage the hypothesis testing formulation (3.1) of differential privacy [19, 64].

3.2

From privacy accounting to numerical computing

For any given k and l, we now cast the privacy accounting of εkl and fkl as a numerical computation problem, as stated in Propositions 3.3 and 3.4. This formulation serves as the cornerstone of our overall privacy accounting analysis for the 2020 Census DHC File. Section 4, which contains the main contribution of this paper, introduces a novel method for computing the (εkl , δ) and fkl curves; the corresponding results and practical recommendations are presented in Section 5. We first leverage Proposition 3.2 in [64] to characterize the (ε, δ)-DP curve in terms of the distribution functions of the privacy-loss log-likelihood ratios. Lemma 3.2 (Proposition 3.2 in [64]). Let Xi ∼ Pi = NZ (0, σi2 ) and Yi ∼ Qi = NZ (µ, σi2 ). The fk , M f l ] satisfy (εkl , δ)-DP curve satisfying composition of two mechanisms [M # "m # "m X X dQi (Xi ) dQi (Yi ) εkl > εkl − e · P log > εkl . δ(εkl ) = P log dPi (Yi ) dPi (Xi ) i=1

i=1

A straightforward calculation yields the following proposition, which reduces the privacy accounting problem to a numerical computation problem involving the accurate evaluation of tail probabilities of weighted convolutions of discrete Gaussian distributions. Section 4 focuses on addressing this numerical challenge. Proposition 3.3 (Formulation of the exact (ε, δ)-DP curve). Recall that ρ = 4.9622 is the total privacy budget in the 2020 Census DHC File. Let ai = ρi /ρ. Then, the (εkl , δ)-DP curve is exactly 11

characterized by the following equation. δ(εkl ) = PXi ∼NZ (0,σ2 )

m X

i

i=1

εkl −1 ai Xi > ρ

! − eεkl · PXi ∼NZ (0,σ2 )

m X

i

i=1

! εkl +1 . ai Xi > ρ

By the Neyman–Pearson lemma (see [40] or Lemma A.1 in [19]), we can reduce the privacy accounting of fkl to a similar numerical computation problem without loss of optimality. Proposition 3.4 (Numerical formulation of the exact f -DP curve). The trade-off function fkl is uniquely determined by the type I error αkl (ζ) and the type II error βkl (ζ). ! ! m m X X αkl (ζ) = PXi ∼NZ (0,σ2 ) ai Xi > ζ + 1 + c · PXi ∼NZ (0,σ2 ) ai Xi = ζ + 1 i

βkl (ζ) = PXi ∼NZ (0,σ2 )

i

i=1 m X

i

i=1

! ai Xi ≤ ζ − 1

− c · PXi ∼NZ (0,σ2 )

i=1 m X

i

i=1

! ai Xi = ζ − 1

for some constant c. The proofs of Propositions 3.3 and 3.4 are provided in Section B. A direct consequence of Propositions 3.3 and 3.4 is that privacy accounting for the composition of discrete or continuous Gaussian mechanisms is mathematically equivalent to computing the tail probability of a convolution of discrete Gaussian distributions, namely ! m X PXi ∼NZ (0,σ2 ) ai Xi > t0 (3.2) i

i=1

for any threshold t0 . In the next section, we propose an efficient method that can compute the (εkl , δ)-DP and fkl -DP curves under a prescribed error tolerance. Before introducing our method, we emphasize the difficulty of numerically computing (3.2). The first challenge arises from the stringent accuracy requirement. In the 2020 Census DHC File [3], an illustrative privacy level is given by (ε, δ) = (26.34, 10−10 ). According to Proposition 3.3, an effective numerical error tolerance should be much smaller than δ/eε . In Section 5.2, we therefore set the tolerance level to ∆ = 10−35 to ensure that numerical errors do not affect the accuracy of the resulting (ε, δ)-DP curve for the 2020 Census DHC File. Although the f -DP curve imposes less stringent accuracy requirements, we set ∆ = 10−25 in Section 5.1 to maintain high precision, particularly near the extreme points of the curve. The second difficulty arises from the distinction between discrete and continuous Gaussian distributions. Empirically, Figure 1 illustrates the probability mass function of the convolution of two discrete Gaussian distributions and compares it with the corresponding continuous Gaussian density, both in the form of (3.2). Notably, the probability mass function of the discrete Gaussian convolution oscillates around the continuous Gaussian density, and a visible discrepancy can be observed between the two curves. Adopting ρ-zCDP for privacy accounting, as in the 2020 Census DHC File, is intuitively equivalent to approximating the discrete Gaussian mechanism by a continuous Gaussian one, which generally yields only an upper bound on the true privacy level [36]. 12

Section 6 provides a more detailed discussion from a theoretical perspective.

Figure 1: Probability mass function of the convolution of two discrete Gaussian distributions and its continuous Gaussian approximation.

4

Method

We now propose a new method for computing the (εkl , δ)-DP curve and the fkl -DP curve to any prescribed error tolerance for fixed k and l. By Propositions 3.3 and 3.4, this is equivalent to developing a method that computes the upper tail probability of a weighted sum of independent, heterogeneous discrete Gaussian random variables with a uniform error tolerance ∆, as in (3.2), for any threshold t0 . Our method consists of four steps, each with a provable guarantee that the overall error does not exceed ∆. At the end of this section, we summarize the complete method in Theorem 4.7. P Step 1: Tail bound. The goal of this step is to truncate m i=1 ai LXi so that the tail probability beyond the truncation is at most ∆/4. To this end, note that all ai ∈ Q+ in the 2020 Census DHC File. Hence, there exists a unique L ∈ Q+ such that ai L ∈ N for all i and gcd(a1 L, · · · , am L) = 1. By Bézout’s identity, it follows that ! m X support ai LXi = Z. i=1

We now state a practical choice of the truncation level U depending on L. The proof, given in Section C.1, is based on a sub-Gaussian concentration inequality. Proposition 4.1 (Practical choice of truncation U ). Let U be the smallest integer larger than v ! u m X u 2 2 t L · −2 (4.1) ai σi · log (∆/8), i=1

13

we have PXi ∼NZ (0,σ2 )

m X

i

! ai Xi > t0

i=1

− PXi ∼NZ (0,σ2 )

U≥

i

m X i=1

! ai LXi > t0 · L

≤ ∆/4,

P With this truncation, instead of computing the probability mass function of m over an i=1 ai LX Pmi infinite support beyond t0 , it suffices to consider only the finite set of values for which i=1 ai LXi lies between t0 L and U . Step 2: Discrete Fourier transform. The goal of this step is to express the truncated cumulaP tive probability of m i=1 ai LXi as a single integral, so that it can be evaluated numerically in later steps. We begin with the following decomposition: ! ! m U m X X X PXi ∼NZ (0,σ2 ) U ≥ ai LXi > t0 · L = PXi ∼NZ (0,σ2 ) ai LXi = t i

i

i=1

i=1

t=⌈t0 ·L⌉

1 2π

(a)

=

Z π U X

e−iζt

t=⌈t0 ·L⌉ −π

m Y

fai LXi (ζ)dζ,

i=1

where fai LXi (ζ) = Eeiζai LXi denotes the characteristic function of ai LXi , and step (a) follows from the inverse discrete Fourier transform (see Exercise 3.3.2(iii) in [20]). The characteristic function has the closed-form expression 2

fai LXi (ζ) = Eeiζ·ai LXi = fXi (ai Lζ) =

P∞

iu·ai Lζ · e u=−∞ e

P∞

u=−∞ e

Let X =

Pm

iζX = i=1 ai LXi . Since Ee

1 2π

Z π U X t=⌈t0 ·L⌉ −π

−iζt

e

m Y

2 − u2 2σ i

− u2 2σ

i

.

Qm

i=1 fai LXi (ζ) is an even function of ζ, we obtain U X

Z π

1 −π 2π

fai LXi (ζ)dζ =

i=1

1 −π 2π

U X t=⌈t0 ·L⌉

m Y

fai LXi (ζ)dζ

i=1

t=⌈t0 ·L⌉

Z π =

e−iζt

cos(ζt) · EeiζX .

Define  F (ζ) =

1  · 2π

U X

t=⌈t0 ·L⌉

We conclude Step 2 with the following result.

14

 cos(ζt) · EeiζX .

(4.2)

Proposition 4.2. With F (ζ) defined in (4.2), we have U≥

PXi ∼NZ (0,σ2 ) i

m X i=1

! ai LXi > t0 · L

Z π F (ζ)dζ.

= −π

where U is chosen according to Proposition 4.1. Step 3: Trapezoidal rule. With the integral representation in Proposition 4.2, we next compute Rπ −π F (ζ)dζ numerically. Since F (ζ) is 2π-periodic, we may equivalently write Z 2π

Z π

F (ζ)dζ.

F (ζ)dζ = −π

0

We adopt the trapezoidal rule, together with the error bound in Lemma 4.3. For a number of nodes N ∈ N, the trapezoidal approximation takes the form N

IN =

2π X F (ζk ), N

(4.3)

k=1

where ζk = 2πk/N . The following lemma, due to [52], shows that for periodic complex analytic functions the trapezoidal rule achieves exponential convergence in N . Lemma 4.3 (Theorem 3.2 in [52]). Consider a 2π-periodic complex analytic function F (ζ) satisfying |F (ζ)| ≤ M on the strip −a < Im ζ < a for some a > 0. Then, for any N ≥ 1, N

2π X F (ζk ) − N k=1

Z 2π 0

F (ζ)dζ ≤

4πM eaN − 1

,

and the constant 4π is optimal. The trapezoidal rule is simple to implement and converges dramatically faster. Empirically, the convergence is often even faster than exponential (see, e.g., [45, 52]). To our knowledge, this phenomenon has not been explicitly exploited in differential privacy, and it suggests a general acceleration strategy for previous numerical accountants [37, 27]. Based on Lemma 4.3, we obtain the following practical choice of N . Proposition 4.4 (Practical Choice of N ). For any N ∈ N satisfying s P N ≥ 2 · (U + 1) + dL log(2 · dL) we have

N

2π X F (ζk ) − N k=1

Z 2π 0

with

F (ζ)dζ ≤ ∆/4,

15

m 2 2 i=1 ai σi

d=

−

with

ζk =

2 log(∆/8)

2πk . N

,

(4.4)

Step 4: Sieve-type truncation. The remaining challenge is that evaluating F (ζ) at all N nodes may be computationally expensive. With the stringent tolerance level ∆ = 10−35 used in our experiments in Section 5.2, the required N is often ≫ 106 , depending on the choice of k and l. To improve efficiency, we introduce an additional truncation of the trapezoidal sum by identifying a subset C ⊆ [N ] that satisfies the following desiderata: First, F (ζk ) is negligible for all k ∈ [N ] \ C. Specifically, we want 2π X F (ζk ) ≤ ∆/4, N k∈[N ]\C

P PN so that 2π approximation to 2π k∈C F (ζk ) serves as an accurateP k=1 F (ζk ). Second, the set C N N is as small as possible, so that evaluating 2π F (ζ ) is substantially faster than evaluating the k k∈C N full sum over [N ]. Finally, the procedure for constructing C is itself computationally efficient. To this end, we propose Algorithm 1 to construct C ⊆ [N ], which is theoretically motivated by sieve methods in number theory. In our implementation, standard algorithmic tools (e.g., dynamic 2 ) programming and binary search) are used to accelerate the selection of C. Let X(i) ∼ NZ (0, σ(i) denote the random variable among {Xi ∼ NZ (0, σi2 )}m i=1 with the ith largest noise parameter σi2 . Let a(i) be the corresponding weight in {ai }m , and write fX(i) (a(i) Lζ) for the characteristic i=1 function of a(i) LX(i) . Algorithm 1 proceeds as follows. Starting from X(1) with noise parameter 2 , we discard all nodes k ∈ [N ] for which the characteristic function is sufficiently small, namely σ(1) those satisfying fX(1) (a(1) Lζ) < ∆/(8 · U ). By the argument in Section C, for any ζ satisfying fX(1) (a(1) Lζ) < ∆/(8 · U ), we have 2π · |F (ζ)| ≤ ∆/4. The algorithm then repeats the same filtering step sequentially for X(2) , . . . , X(m) , and the surviving indices constitute the final set C. Mathematically, C is defined to be C=

m \



[

i=1 0≤j≤a(i) L

2π 2π · j − ζ (i) , · j + ζ (i) a(i) L a(i) L



The following proposition concludes Algorithm 1 by showing that restricting the trapezoidal sum to the indices in C incurs an error of at most ∆/4. The proof of Proposition 4.5 is deferred to Appendix C.3. Proposition 4.5. With C output by Algorithm 1, we have N

2π X 2π X F (ζk ) − F (ζk ) ≤ ∆/4, N N k=1

k∈C

with

ζk =

2πk . N

We now provide further intuition behind Algorithm 1. We begin with several basic properties of the characteristic function fXi (ζ) for all i. Lemma 4.6. The characteristic function fXi (ζ) satisfies the following properties: 1. |fXi (ζ)| ≤ 1, and fXi (ζ) attains its maximum at ζ = 0 with fXi (0) = 1. 2. fXi (ζ) is 2π-periodic; consequently, fXi (ai Lζ) is periodic with period 2π/(ai L). 16

Algorithm 1: Truncation: computationally efficient to find C. Input: Number of nodes N ; collection of all nodes {2πk/N : k ∈ [N ]}; discrete Gaussian 2 ; f distribution X(i) with i-th largest parameter σ(i) X(i) (a(i) Lζ) to be characteristic function of a(i) LX(i) . Candidate Nodes = [N]

// initialization

for i = 1 to m do ∆ Search for the smallest a such that ζ (i) = a · a π L satisfies fX(i) (a(i) Lζ (i) ) < 8U ; (i)

// If X(i) has r i.i.d. copies among {Xj }m j=1 , then replace the test by fX(i) (a(i) Lζ (i) )

r

∆ < 8U .

Update Candidate Nodes by

// Keep ONLY the nodes that are close to peak

Candidate Nodes = Candidate Nodes  [ \ 2πk ∈ k ∈ [N ] :  N

0≤j≤a(i) L



 

2π 2π . · j − ζ (i) , · j + ζ (i)  a(i) L a(i) L

(4.5)

end Output: C = Candidate Nodes 3. fXi (ζ) is strictly increasing on (−π, 0) and strictly decreasing on (0, π). Lemma 4.6 shows that |fXi (ai Lζ)| exhibits periodic peaks at a 2πL · k for all integers k. The (i) key idea of Algorithm 1 is therefore to sieve (i.e., discard) those nodes at which at least one characteristic function fX(i) (a(i) Lζ) is sufficiently small, exploiting this periodic structure. In each iteration, we first determine a radius ζ (i) such that |fXi (ai Lζ)| < ∆/(8 · U ) whenever ζ ∈ / [0, ζ (i) ]. 2π By periodicity, any ζk whose distance to one of the points a L · j is greater than ζ (i) satisfies (i) |fXi (ai Lζ)| < ∆/(8 · U ), which leads to the update rule in (4.5). In Section C, we show that removing all nodes ζk outside the union described in (4.5) incurs an error of at most ∆/4. Hence, Algorithm 1 enables the efficient yet rigorous elimination of nodes whose values of F are negligible. We now provide further intuition for starting the iteration with X(1) , which has the largest noise 2 . When N , or the number of Candidate Nodes, is large, it can be computationally parameter σ(1) expensive to identify all nodes whose distance to the set { a 2πL · j : 0 ≤ j ≤ a(i) L} is less than (i)

ζ (i) .

2 . Observe that a larger noise We therefore begin with X(1) , corresponding to the largest σ(1) 2 corresponds to a smaller weight a . Consequently, the update rule in (4.5) involves parameter σ(1) (1) the smallest number of periodic unions, each with the largest interval length. Empirically, this ordering provides the fastest reduction in the set of Candidate Nodes. Algorithm 1 substantially accelerates the computation. Empirically, the dominant computational cost lies in evaluating the function F (ζ). Under the allocation in Table 1, computing Rπ 6 −π F (ζ)dζ without Algorithm 1 would require evaluating F (ζk ) at all k ∈ [N ] with N = 10 . 7 Consequently, each subplot in Figure 4 would involve more than 4 × 10 evaluations of F (ζ) at

17

high precision, which is infeasible on a single CPU. With Algorithm 1, our method evaluates F (ζk ) only for k ∈ C, where |C| = 203. As a result, each subplot in Figure 4 requires only about 8,120 evaluations of F (ζ), and can be completed within approximately 0.5 hours. We now summarize the entire procedure in the following theorem. Theorem 4.7. Equation (3.2) can be numerically evaluated with error tolerance ∆ using ! m X 2πk 2π X F (ζk ) ≤ ∆ with ζk = , PXi ∼NZ (0,σ2 ) ai Xi > t0 − i N N i=1

(4.6)

k∈C

where  F (ζ) =

1  2π

U X

t=⌈t0 ·L⌉

 1 cos(ζt) · EeiζX = · weight(ζ) · char(ζ), 2π

and weight(ζ) and char(ζ) have an explicit evaluable expression as the following when ζ ̸= 0 or 2π. U X

  cos(ζ/2) 1 cos(⌈t0 L⌉ · ζ) + cos(U · ζ) + (sin(U · ζ) − sin(⌈t0 L⌉ · ζ)) , weight(ζ) = cos(ζt) = 2 sin(ζ/2) t=⌈t0 ·L⌉ ) ( ∞ )−1 ( ∞ m Y X X 2 2 2 2 . char(ζ) = EeiζX = e−xi /2σi fXi (ai Lζ), with fXi (ai Lζ) = eiai Lζxi · e−xi /2σi · xi =−∞

i=1

xi =−∞

Here U , N , and C are chosen according to Propositions 4.1, 4.4, and 4.5, respectively. Theorem 4.7 has value beyond the specific application to the 2020 Census DHC File. It introduces a novel numerical method for privacy accounting that provably achieves exponential convergence, representing a substantial improvement over first-order methods [38, 27] and existing higher-order approaches [48], making it broadly applicable to a wide range of privacy accounting problems. Moreover, for privacy accounting tasks that require a large number of nodes, strategies analogous to Algorithm 1 can be employed to dramatically reduce computational cost.

5

Results

We extract the privacy budget allocations for all directed paths from the 2020 Census DHC File. Although there are millions of directed paths [4], the 2020 Census DHC File employs only 43 distinct privacy budget allocations across all paths. Section 4 provides an efficient and practical procedure for computing fkl and (εkl , δ) for any pair of directed paths k and l. Recall from Proposition 3.1 that the overall privacy guarantee is given by (mink,l fkl )∗∗ -DP and (maxk,l εkl , δ)-DP. To compute these quantities, we leverage parallel computation on a CPU-based infrastructure and launch 43 × 42/2 + 43 = 946 jobs simultaneously, where each job evaluates the privacy level fkl or (εkl , δ) for a given pair k, l ∈ [43]. For the f -DP guarantee, we compute the lower convex envelope of mink,l fkl using the monotone chain algorithm, with results presented in Section 5.1. The overall (ε, δ)-DP curve of the 2020 Census DHC File is then obtained by taking the pointwise maximum over (εkl , δ) across 18

all pairs (k, l) and the results are reported in Section 5.2. In practice, the U.S. Census Bureau tests multiple privacy budgets during the tuning stage and evaluates their performance on downstream tasks [56, 2, 60, 61]. Motivated by our empirical findings, Section 5.3 provides practical guidance for efficiently applying our method during privacy budget tuning, thereby avoiding repeated full-scale parallel computation of the overall privacy accounting. 1.0

0

Bureau’s Ours

−2

0.8

0.6 log10(δ)

type II error

−4

0.4

−6 −8

−10

0.2

0.0 0.0

−12 −14

0.2

0.4 0.6 type I error

0.8

1.0

10

Bureau’s Ours 15

20 

25

30

Figure 2: The f -DP curves and (ε, δ)-DP curve for the 2020 Census DHC File computed using our method (blue) and the Census Bureau’s accounting method (red).

5.1

f -DP curves

The left panel of Figure 2 characterizes the overall privacy guarantee of the 2020 Census DHC File in terms of f -DP curves. The blue curve represents the exact f -DP guarantee computed using our method, while the red curve is obtained by converting the ρ-zCDP guarantee to f -DP using AutoDP [65, 66, 67]. Under the same mechanism and injected noise, our method yields a strictly better privacy–utility trade-off. Figure 3 presents 43 subplots illustrating the privacy guarantees of the f0, M f l ], where M f 0 corresponds to the privacy budget allocation in Table composed mechanisms [M fk , M f l ] are available at https://github. 1. The privacy guarantees of all composed mechanisms [M com/BuxinSu/Exact-Privacy-Accounting-for-2020-U.S.-Census.git. Although there are 946 distinct composed mechanisms, their privacy profiles are different yet highly similar, as illustrated in Figure 3.

5.2

(ε, δ)-DP curves

Since the (ε, δ)-DP curve is equivalent to the f -DP curve and has been suggested as an alternative way to report privacy guarantees [26], we present the exact (ε, δ)-DP curve for the 2020 Census DHC File in the right panel of Figure 2. Similarly, Figure 4 illustrates the privacy guarantees of f0, M f l ] for all l, expressed in terms of their (ε, δ)-DP curves. the composed mechanisms [M

19

0.0 0.0

0.0 0.0

0.0 0.0

0.4 0.2

0.5 type I error

0.0 0.0

1.0

0.5 type I error

0.6 0.4

0.4

0.0 0.0

1.0

Bureau’s Ours

0.8

0.6 0.4

Bureau’s Ours

0.8

0.6 0.4 0.2

0.5 type I error

0.0 0.0

1.0

1.0 type II error

type II error

0.8

0.6

0.2

1.0 Bureau’s Ours

0.8

1.0 Bureau’s Ours

0.8

0.6 0.4

0.8

0.6 0.4 0.2

0.2

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

Bureau’s Ours

0.8

0.6 0.4

0.5 type I error

1.0

1.0 Bureau’s Ours

0.8

0.6 0.4

0.5 type I error

1.0

1.0 Bureau’s Ours

0.8

0.6 0.4

Bureau’s Ours

0.8

0.6 0.4

0.2

0.2

0.2

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

0.2

type II error

0.8

0.2 0.5 type I error

Bureau’s Ours

0.4

0.5 type I error

0.6 0.4

0.2 0.0 0.0

0.5 type I error

1.0

type II error

0.8

0.6 0.4

0.2 0.0 0.0

0.5 type I error

1.0

1.0

0.5 type I error

1.0

1.0

0.5 type I error

0.4

type II error

0.0 0.0

1.0

0.5 type I error

Bureau’s Ours

Bureau’s Ours

0.8 0.6 0.4 0.2

0.5 type I error

0.0 0.0

1.0

0.5 type I error

0.5 type I error

Bureau’s Ours

Bureau’s Ours

0.8 0.6 0.4 0.2

0.5 type I error

1.0

0.0 0.0

0.5 type I error

1.0

1.0

Bureau’s Ours

0.6

Bureau’s Ours

0.6 0.4 0.2

0.4

0.0 0.0

0.5 type I error

1.0

0.0 0.00

0.25

1.0

1.0

0.4

0.8

1.0

1.0

0.6

0.0 0.0

0.6

0.2

0.2

Bureau’s Ours

0.2 0.5 type I error

type II error

0.4

1.0

Bureau’s Ours

0.8

0.4

0.8

0.6

0.0 0.0

Bureau’s Ours

1.0 Bureau’s Ours

0.5 type I error

1.0

0.6

0.0 0.0

1.0

0.2 0.5 type I error

0.0 0.0

1.0

0.2 0.5 type I error

0.4

0.4

0.8

0.4

0.0 0.0

0.8

1.0 Bureau’s Ours

0.6

0.2

0.4

0.6

0.0 0.0

1.0

1.0 Bureau’s Ours

Bureau’s Ours

0.2 0.5 type I error

0.8

1.0

0.4

0.8

type II error

1.0

Bureau’s Ours

0.4

0.0 0.0

Bureau’s Ours

0.6

0.0 0.0

1.0

0.5 type I error

1.0

1.0

0.6

0.2

0.5 type I error

0.6

0.0 0.0

1.0

1.0

0.2

0.2

type II error

type II error

0.8

0.4

0.8

1.0 Bureau’s Ours

0.8

1.0

0.6

0.0 0.0

1.0

Bureau’s Ours

0.6

0.0 0.0

1.0

0.2

1.0

0.8

0.5 type I error

0.5 type I error

1.0

0.2

type II error

Bureau’s Ours

0.4

0.8

0.4

1.0

0.6

0.0 0.0

0.8

0.6

0.0 0.0

1.0

1.0 0.8

Bureau’s Ours

0.2 0.5 type I error

1.0

type II error

0.4

0.5 type I error

1.0

type II error

0.8

0.6

0.0 0.0

1.0

type II error

Bureau’s Ours

type II error

0.8

0.5 type I error

1.0

type II error

1.0

type II error

0.5 type I error

0.6

0.2

0.4

0.2

1.0

Bureau’s Ours

0.8

0.6

0.2

1.0

Bureau’s Ours

1.0 type II error

1.0 type II error

0.8

1.0

type II error

1.0

0.5 type I error

type II error

1.0

0.5 type I error

1.0

0.4

0.2

0.4

0.0 0.0

1.0

0.6

0.2

0.6

0.2 0.5 type I error

1.0 Bureau’s Ours

1.0

Bureau’s Ours

0.8

0.4

0.2 0.5 type I error

Bureau’s Ours

0.6

0.0 0.0

1.0

0.5 type I error

1.0

0.2 0.5 type I error

0.4

0.0 0.0

1.0

type II error

0.6

Bureau’s Ours

0.5 type I error

1.0

type II error

0.8

1.0

0.6

0.2

type II error

0.4

Bureau’s Ours

0.5 type I error

1.0 type II error

0.6

1.0

type II error

0.8

0.5 type I error

1.0

type II error

1.0

type II error

Bureau’s Ours

type II error

type II error

0.5 type I error

1.0

type II error

1.0

Bureau’s Ours

0.8

0.4

0.0 0.0

0.5 type I error

1.0 Bureau’s Ours

0.6

0.0 0.0

0.8 type II error

0.4

0.2

1.0

type II error

0.6

0.2

0.0 0.0

type II error

0.4

0.8

0.2

0.2

type II error

0.6

1.0 Bureau’s Ours

0.2

0.8

type II error

0.4

0.8

0.2

1.0

type II error

0.6

1.0 Bureau’s Ours

type II error

0.4

0.8

type II error

0.6

1.0 Bureau’s Ours

type II error

0.8

type II error

1.0 Bureau’s Ours

type II error

type II error

0.8

type II error

1.0

0.50 0.75 type I error

1.00

f0, M f l ] for all l, where M f 0 corresponds Figure 3: The f -DP curves for the composed mechanisms [M to the privacy budget allocation in Table 1. Although the curves appear visually similar, they differ by amounts significantly larger than the numerical error introduced by floating-point arithmetic or the mpmath package, as quantified in Figure 5. The enlarged panel in the lower-right corner displays f0, M f 0 ]. the privacy profile of the composed mechanism [M 20

1.0

−5

−5

−5

−5

−5

−5

−10

−15 10

20 

30

−15 10

20 

30

−15 10

20 

−5

−10

20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

log10(δ)

−5

log10(δ)

−5

log10(δ)

−5

log10(δ)

−5

log10(δ)

−5

Bureau’s Ours

−10

Bureau’s Ours 20 

30

20 

30

−5

−5

−5

−5

−5

−5

−10

20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

−10

Bureau’s Ours 20 

30

−15 10

20 

30

−5

−5

−5

−5

−5

−5

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

−10

Bureau’s Ours 20 

30

−15 10

20 

−5

−5

−5

−5

−5

−5

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

−10

Bureau’s Ours 20 

30

−15 10

20 

0

−5

−5

−5

−5

−5

−5

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

0

0.0

−5

−5

−5

−2.5

−15 10

−10

Bureau’s Ours 20 

30

−15 10

−10

Bureau’s Ours 20 

30

−15 10

0

0

−5

−5

−5

−10 −15 10

−10

Bureau’s Ours 20 

30

−15 10

log10(δ)

0

−10

Bureau’s Ours 20 

30

−15 10

Bureau’s Ours 20 

30

log10(δ)

0

log10(δ)

0

−10

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

−15 10

20 

20 

−10

Bureau’s Ours

20 

−5.0 −7.5

−10.0

Bureau’s Ours 20 

30

−12.5

Bureau’s Ours

−15.0 10

15

20 

25

30

Bureau’s Ours

−15 10

30

30

Bureau’s Ours

−15 10

30

0

−10

20 

−10

Bureau’s Ours

30

Bureau’s Ours

−15 10

30

0

−10

20 

−10

Bureau’s Ours

30

Bureau’s Ours

−15 10

0

−10

20 

−10

Bureau’s Ours

30

Bureau’s Ours

−15 10

0

Bureau’s Ours

20 

−10

Bureau’s Ours

−15 10

Bureau’s Ours

−15 10

30

0

log10(δ)

log10(δ) log10(δ)

30

0

−15 10

log10(δ)

20 

−10

Bureau’s Ours

0

−10

log10(δ)

−15 10

−10

Bureau’s Ours

0

−15 10

log10(δ)

30

−10

Bureau’s Ours

0

−10

log10(δ)

20 

−10

Bureau’s Ours

0

log10(δ)

log10(δ)

−15 10

−10

Bureau’s Ours

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

0

log10(δ)

log10(δ)

0

30

f0, M f l ] for all l, where M f 0 corFigure 4: The (ε, δ)-DP curves for the composed mechanisms [M responds to the privacy budget allocation in Table 1. Although the curves appear highly similar, they differ by amounts significantly larger than the numerical error introduced by floating-point arithmetic or the mpmath package, as detailed in Figure 5. The enlarged panel in the lower-right f0, M f 0 ]. corner displays the privacy guarantee of the composed mechanism [M 21

30

5.3

Practical advice on privacy budget tuning

In practice, the U.S. Census Bureau considers multiple privacy budget allocations prior to deployment and evaluates their performance on downstream tasks [56, 2, 60, 61]. However, determining the privacy budget based on the overall privacy guarantee can be computationally expensive, as it requires large-scale parallel computation. We therefore provide practical guidance to avoid repeatedly recomputing the full privacy accounting. As shown in Figures 3 and 4, the privacy guarantees fk , M f l ] are not identical but are highly similar. Motivated of different composed mechanisms [M f0, M f 0 ], with the privacy by this observation, we recommend using the composed mechanism [M budget allocation given in Table 1, as a proxy during the privacy budget tuning stage.5 After the tuning process is completed, the exact overall privacy guarantee can then be computed to ensure full mathematical rigor. To make this recommendation more precise, we first quantify how the privacy guarantee of the f0, M f 0 ] differs from the overall privacy guarantee. We measure the relative composed mechanism [M difference between the two as f00 − (mink,l fkl )∗∗ , f00

and

ε00 − maxk,l εkl . ε00

Figure 5 shows that, although the overall privacy level differs slightly from the value computed using the allocation in Table 1, the discrepancy is empirically confined to a very small range. Specifically, the overall privacy guarantee is observed to uniformly lie within max εkl ∈ [ε00 , (1 + 0.005)ε00 ] k,l   (min fkl )∗∗ ∈ (1 − 1.0 × 10−4 )f00 , f00 .

(5.1)

k,l

Although (5.1) is not a formal theoretical guarantee, it provides strong empirical evidence that the overall privacy level is relatively insensitive to the choice among different privacy budget allocations. This observation motivates the following practical strategy for privacy budget tuning. In practice, the Census Bureau may use the privacy budget allocation in which no geographic level is bypassed, f0, M f 0 ], as a proxy during the tuning phase. After tuning is complete, the composed mechanism [M a final evaluation of the full overall privacy guarantee can be performed, which is expected to yield very similar performance.

6

Approximation Gap from the Continuous Gaussian

In this section, we discuss why exact privacy accounting for compositions of the discrete Gaussian mechanism, both in this paper and in prior work [48], differs substantially from that of the (continuous) Gaussian mechanism [21, 8, 19]. Propositions 3.3, 3.4 and Theorem 8 in [8] show that privacy accounting for compositions of discrete and continuous Gaussian mechanisms is mathematically equivalent to evaluating the tail probability of weighted convolutions of discrete and continuous Gaussian distributions, respectively. Empirically, Figure 1 illustrates the probability mass function 5

More precisely, we recommend using a privacy budget allocation in which no geographic level is bypassed.

22

8

×10−5

0.000 −0.001

6

Percentage difference

percentage difference

7

5 4 3 2 1 0 0.0

−0.002 −0.003 −0.004 −0.005

0.2

0.4 0.6 type I error

0.8

1.0

10

15

20 

25

30

f0, M f0] Figure 5: Percentage difference between the privacy level of the composed mechanism [M and the overall privacy level in terms of f -DP (left) and (ε, δ)-DP (right). of the convolution of two discrete Gaussian distributions and its continuous Gaussian counterpart, revealing a pronounced discrepancy between the discrete convolution and the smooth Gaussian density curve. To gain further intuition, we examine the difference from the perspective of the central limit theorem. Even under i.i.d. compositions and equal weight, the convolution of discrete Gaussian variables produces a distribution that differs from its continuous Gaussian counterpart. For i.i.d. continuous Gaussian variables, the normalized sum has density ! 2 m X 1 e−y /2 d PXi ∼N (0,σ2 ) p Xi = y = √ . Pm 2π var ( i=1 Xi ) i=1 In contrast, a second-order CLT refinement—known as the Edgeworth expansion [28, 11, 17]—yields the following approximation for the m-fold convolution of discrete Gaussian variables: !   2 m X 1 1 e−y /2 λ4 2 4 PXi ∼NZ (0,σ2 ) p · √ 1+ (3 − 6y + y ) , Xi = y ≈ p P P 24m 2π var ( m var ( m i=1 Xi ) i=1 Xi ) i=1

λ4 (3 − 6y 2 + y 4 ) captures the discrepancy between the discrete and where the additional term 24m continuous Gaussian mechanisms. For the noise parameters used in the 2020 Census DHC File, the magnitude of λ4 /(24m) is not negligible. The illustrative result in Fact A.2 shows that λ4 /(24m) ≫ ∆ under some of the adopted noise parameters. Theorem 6.1 formalizes this observation by establishing a rigorous uniform error bound of order m−3/2 . p P P i.i.d. Let Xi ∼ NZ (0, σ 2 ) and define Sm = B1m m var( m i=1 Xi , where Bm = i=1 Xi ). As illustrated in Fact A.3 of [48], var(Xi ) is often close to, but not identical to, the noise parameter σ 2 . Let

23

κr = κr (Xi ) denote the r-th cumulant of Xi , and define r/2

λr = κr /κ2 ,

r/2

Kr = E|Xi |r /κ2 .

Using these quantities, we can derive an explicit uniform error bound for the Edgeworth expansion. Theorem 6.1 (Uniform error bound for the Edgeworth expansion). For any y in the support of Sm , define the local error   2 λ4 1 e−y /2 2 4 · √ 1+ (3 − 6y + y ) . Rm (y) := P (Sm = y) − Bm 24m 2π Then, there exists a constant C > 0 such that sup Rm (y) ≤

y∈Z/Bm

C . m3/2

(6.1)

More precisely, let U m (c) and Um (ζ) be defined in (D.6) and (D.10). For any c satisfying U m (c) < 1, we have,   Z ∞ 1 λ4 ζ 4 −1/2 −ζ 2 /2 · sup Rm (y) ≤ m dζ e 1+ √ 1/2 24m c· log(m) y∈Z/Bm πκ2  2   Z c·√log(m) 1 ζ −1/2 ′ ′ ′ +m · exp − × Λ5 + Λ6 + Λ7 dζ (6.2) 1/2 2 0 πκ2 1/2 !  3 c2 cos(θ0 ) − 1 + θ02 /2 log (m) + exp − log(m) + · 8c3 K3 · , 2 m 2θ03 where Λ′5 , Λ′6 , Λ′7 are given by     ζ4 1 K6 ζ 6 1 K6 ζ 6 K8 ζ 8 K10 ζ 10 K12 ζ 12 · +1 + + · + + + 8m (1 − U m (c))2 6! m2 (1 − U m (c))2 48m2 640m3 17280m4 1036800m5  2 1 K4 ζ 4 K6 ζ 6 Um (ζ) ′ + + Λ6 = 2 4! m 6! m2 2(1 − U m (c))2 !  p 3 1 c4 K4 log(m)2 c6 K6 log(m)3 Um (c log(m)) K4 ζ 4 K6 ζ 6 Um (ζ) ′ Λ7 = exp + + × + + , 6 4! m 6! m2 2(1 − U m (c))2 4! m 6! m2 2(1 − U m (c))2

Λ′5 =

and θ0 is the unique root of the equation θ2 + 2θ sin(θ) + 6(cos(θ) − 1) = 0,

(6.3)

within the interval (0, 2π). The proof of Theorem 6.1 is postponed to Section D. We now briefly explain why the right-hand side of Equation (6.2) is of order O(m−3/2 ), as stated in Equation (6.1). Recall that Equation (6.2) 24

holds for any c satisfying U m (c) = c2 ·

log(m) K4 + c4 · 2m 4!



log(m) m

2

+ c6 ·

K6 6!



log(m) m

3 <1

(6.4)

When m is sufficiently large, we may choose c = 2. With this choice, we have   Z ∞ λ4 ζ 4 1 −ζ 2 /2 1+ e dζ = o(m−3/2 ) √ 1/2 24m c· log(m) πκ2  3 1/2 ! c2 log (m) cos(θ0 ) − 1 + θ02 /2 3 exp − log(m) + · 8c K3 · = o(m−3/2 ) 2 m 2θ03       1 1 1 ′ ′ ′ , Λ6 = O , Λ7 = O . Λ5 = O m m2 m3 Combining these bounds, we conclude that Rm (y) is bounded above by O(m−3/2 ).

7

Discussion

In this paper, we provide, for the first time, an exact characterization of the overall privacy guarantee of the 2020 Census DHC File in terms of f -DP and its corresponding (ε, δ)-DP profile. Compared with the U.S. Census Bureau’s current upper-bounding approach, these curves offer a more direct and policy-relevant interpretation of re-identification risk and are easier to communicate in empirical research and policymaking across the social sciences [50], political science [6, 14], and economics [7, 55], while achieving the optimal privacy–utility trade-off with the exact computation of the composition bounds. Achieving exact privacy accounting for the discrete Gaussian mechanism presents two key challenges: analytic composition methods developed for the continuous Gaussian setting [32, 68] do not extend to the discrete case, and exact (ε, δ)-DP and f -DP curves require extremely high numerical precision together with large-scale parallel composition (e.g., 946 mechanisms for the 2020 Census DHC File), making existing numerical approaches [37, 27, 48] computationally impractical. To address this, we develop a new framework that integrates numerical analysis and numbertheoretic ideas. Our acceleration rests on three numerical innovations tailored to the heterogeneity and extreme-precision regime of the discrete Gaussian mechanism. We recast exact privacy accounting into a discrete Fourier–based numerical integration problem, and then explicitly exploit the exponential convergence of the trapezoidal rule for complex-analytic periodic integrands. To make this exponential regime usable in the fully heterogeneous setting, we map rational weights onto a common integer lattice, enabling sharp sub-Gaussian truncation bounds. Finally, inspired by sieve methods in number theory, we prune quadrature nodes by rigorously discarding regions where the oscillatory integrand must be negligible, reducing the required evaluations from ∼ 106 nodes down to as few as 203 in representative DHC settings; combined with practical tuning guidance, this makes exact privacy characterization feasible at scale. More broadly, our contribution is a fast and numerically reliable framework for exact privacy accounting in large-composition settings. This makes accurate privacy characterization practically 25

feasible and enables policymakers to treat privacy cost as a tunable resource when allocating budget across releases. Beyond the discrete Gaussian mechanism and the Census setting, several directions merit further study. For instance, developing a direct analytical construction of the sieve (e.g., via lattice reduction) could eliminate the current sequential search and further accelerate the algorithm. In addition, although the Census Bureau’s current allocation weights are rational, extending the method to handle irrational weights and almost-periodic integrands could enable more flexible privacy budget allocations. Relatedly, integrating the framework with other privacy accounting paradigms may provide accurate privacy budgets with tighter privacy–utility trade-offs at scale. Finally, it may be fruitful to explore whether these techniques can inform the design of improved integer-valued mechanisms for official statistics.

Acknowledgments This work was supported in part by a Meta Faculty Research Award, and Wharton AI for Business. Chendi Wang was supported by the Fundamental Research Funds for the Central Universities (No. 20720251050) and NSFC Basic Science Center Project for Econometric Modeling and Economic Policy Studies (No. 71988101).

References [1] J. M. Abowd. Staring down the database reconstruction theorem, 2019. URL https://www2. census.gov/programs-surveys/decennial/2020/resources/presentations-publications/ 2019-02-16-abowd-db-reconstruction.pdf. [2] J. M. Abowd. census.gov. https://www.census.gov/content/dam/Census/newsroom/press-kits/ 2020/jsm/trying-to-be-a-good-data-steward-in-the-21st-century.pdf, 2020. [3] J. M. Abowd, R. Ashmead, R. Cumings-Menon, S. Garfinkel, M. Heineck, C. Heiss, R. Johns, D. Kifer, P. Leclerc, A. Machanavajjhala, B. Moran, W. Sexton, M. Spence, and P. Zhuravlev. The 2020 census disclosure avoidance system TopDown algorithm. Harvard Data Science Review, (Special Issue 2), 2022. [4] J. M. Abowd, R. Ashmead, R. Cumings-Menon, S. L. Garfinkel, M. Heineck, C. Heiss, R. Johns, D. Kifer, P. Leclerc, A. Machanavajjhala, B. Moran, W. Sexton, M. Spence, and P. Zhuravlev. Invited lecture: The u.s. census bureau adopts differential privacy. In KDD ’18: Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, 2022. [5] M. J. Anderson. The American census: A social history. Yale University Press, 2015. [6] S. Ansolabehere and J. Snyder. The End of Inequality: One Person, One Vote and the Transformation of American Politics. Issues in American democracy. Norton, 2008. ISBN 9780393931037. [7] D. H. Autor and M. G. Duggan. The rise in the disability rolls and the decline in unemployment. The Quarterly Journal of Economics, 118(1):157–206, 2003. [8] B. Balle and Y.-X. Wang. Improving the gaussian mechanism for differential privacy: Analytical calibration and optimal denoising. In International conference on machine learning, pages 394–403. PMLR, 2018. [9] D. Boyd and J. Sarathy. Differential Perspectives: Epistemic Disconnects Surrounding the U.S. Census Bureau’s Use of Differential Privacy. Harvard Data Science Review, (Special Issue 2), 2022.

26

[10] Z. Bu, J. Dong, Q. Long, and W. J. Su. Deep learning with Gaussian differential privacy. Harvard Data Science Review, 2020(23):10–1162, 2020. [11] J. P. Buhler, A. C. Gamst, R. Graham, and A. W. Hales. Explicit error bounds for lattice edgeworth expansions. Connections in Discrete Mathematics: A Celebration of the Work of Ron Graham, pages 321–352, 2018. [12] M. Bun and T. Steinke. Concentrated differential privacy: Simplifications, extensions, and lower bounds. In Theory of Cryptography Conference, pages 635–658. Springer, 2016. [13] C. L. Canonne, G. Kamath, and T. Steinke. The discrete Gaussian for differential privacy. In Advances in Neural Information Processing Systems, volume 33, pages 15676–15688. Curran Associates, Inc., 2020. [14] A. Cohen, M. Duchin, J. Matthews, and B. Suwal. Census TopDown: The impacts of differential privacy on redistricting. In 2nd Symposium on Foundations of Responsible Computing, FORC 2021, June 9-11, 2021, Virtual Conference, volume 192 of LIPIcs, pages 5:1–5:22. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021. [15] R. Cumings-Menon, R. Ashmead, D. Kifer, P. Leclerc, M. Spence, P. Zhuravlev, and J. M. Abowd. Disclosure avoidance for the 2020 census demographic and housing characteristics file. arXiv preprint arXiv:2312.10863, 2023. [16] R. Cumings-Menon, R. Ashmead, D. Kifer, P. Leclerc, J. Ocker, M. Ratcliffe, P. Zhuravlev, and J. Abowd. Geographic spines in the 2020 census disclosure avoidance system. Journal of Privacy and Confidentiality, 14(3), Aug. 2024. [17] A. Derumigny, L. Girard, and Y. Guyonvarch. Explicit non-asymptotic bounds for the distance to the first-order edgeworth expansion. Sankhya A, 86(1):261–336, 2024. [18] T. Dick, C. Dwork, M. Kearns, T. Liu, A. Roth, G. Vietri, and Z. S. Wu. Confidence-ranked reconstruction of census microdata from published statistics. Proceedings of the National Academy of Sciences, 120(8):e2218605120, 2023. [19] J. Dong, A. Roth, and W. J. Su. Gaussian differential privacy. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 84(1):3–37, 2022. [20] R. Durrett. Probability: theory and examples, volume 49. Cambridge university press, 2019. [21] C. Dwork and A. Roth. The algorithmic foundations of differential privacy. Foundations and Trends® in Theoretical Computer Science, 9(3–4):211–407, 2014. [22] C. Dwork and G. N. Rothblum. Concentrated differential privacy. arXiv preprint arXiv:1603.01887, 2016. [23] C. Dwork, K. Kenthapadi, F. McSherry, I. Mironov, and M. Naor. Our data, ourselves: Privacy via distributed noise generation. In Advances in Cryptology-EUROCRYPT 2006: 24th Annual International Conference on the Theory and Applications of Cryptographic Techniques, St. Petersburg, Russia, May 28-June 1, 2006. Proceedings 25, pages 486–503. Springer, 2006. [24] C. Dwork, F. McSherry, K. Nissim, and A. Smith. Calibrating noise to sensitivity in private data analysis. Theory Of Cryptography, Proceedings, 3876:265–284, 2006.

27

[25] S. J. Eckman. Apportionment and redistricting following the 2020 census. https://sgp.fas.org/crs/ misc/IN11360.pdf, 2021. [26] J. F. Gomez, B. Kulynych, G. Kaissis, J. Hayes, B. Balle, and A. Honkela. Gaussian dp for reporting differential privacy guarantees in machine learning. arXiv preprint arXiv:2503.10945, 2025. [27] S. Gopi, Y. T. Lee, and L. Wutschitz. Numerical composition of differential privacy. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, pages 11631–11642, 2021. [28] P. Hall. The bootstrap and Edgeworth expansion. Springer Science & Business Media, 2013. [29] S. Haney, S. Berghel, B. Carlson, R. Cumings-Menon, L. Hartman, M. Hay, A. Machanavajjhala, G. Miklau, A. Pai, S. Rajpal, D. Pujol, W. Sexton, R. Shrestha, and D. Simmons-Marengo. Safetabp: Disclosure avoidance for the 2020 census detailed demographic and housing characteristics file a (detailed dhc-a). arXiv preprint arXiv:2505.01472, 2025. [30] M. Hawes. Reconstruction and re-identification of the Demographic and Housing Characteristics File (DHC), 2022. URL https://www2.census.gov/about/partners/cac/sac/meetings/2022-03/ presentation-reconstruction-and-reidentification-of-the-dhc.pdf. [31] M. Hotchkiss and J. Phelan. Uses of Census Bureau data in federal funds distribution: A new design for the 21st century. United States Census Bureau, 2017. [32] P. Kairouz, Z. Liu, and T. Steinke. The distributed discrete gaussian mechanism for federated learning with secure aggregation. In Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learning Research, pages 5201–5212. PMLR, 2021. [33] C. T. Kenny, S. Kuriwaki, C. McCartan, E. T. R. Rosenman, T. Simko, and K. Imai. The use of differential privacy for census data and its impact on redistricting: The case of the 2020 U.S. census. Science Advances, 7(41):eabk3283, 2021. [34] C. T. Kenny, S. Kuriwaki, C. McCartan, E. T. R. Rosenman, T. Simko, and K. Imai. Comment: The Essential Role of Policy Evaluation for the 2020 Census Disclosure Avoidance System. Harvard Data Science Review, (Special Issue 2), 2023. [35] C. T. Kenny, C. McCartan, T. Simko, and K. Imai. Census officials must constructively engage with independent evaluations. Proceedings of the National Academy of Sciences, 121(11):e2321196121, 2024. [36] D. Kifer, J. M. Abowd, R. Ashmead, R. Cumings-Menon, P. Leclerc, A. Machanavajjhala, W. Sexton, and P. Zhuravlev. Bayesian and frequentist semantics for common variations of differential privacy: Applications to the 2020 census. arXiv preprint arXiv:2209.03310, 2022. [37] A. Koskela, J. Jälkö, and A. Honkela. Computing tight differential privacy guarantees using FFT. In The 23rd International Conference on Artificial Intelligence and Statistics, AISTATS 2020, 26-28 August 2020, Online [Palermo, Sicily, Italy], volume 108 of Proceedings of Machine Learning Research, pages 2560–2569. PMLR, 2020. [38] A. Koskela, J. Jälkö, L. Prediger, and A. Honkela. Tight differential privacy for discrete-valued mechanisms and for the subsampled Gaussian mechanism using FFT. In The 24th International Conference on Artificial Intelligence and Statistics, AISTATS 2021, April 13-15, 2021, Virtual Event, volume 130 of Proceedings of Machine Learning Research, pages 3358–3366. PMLR, 2021.

28

[39] B. Kulynych, J. F. Gomez, G. Kaissis, J. Hayes, B. Balle, F. Calmon, and J. L. Raisaro. Unifying re-identification, attribute inference, and data reconstruction risks in differential privacy. In The Thirtyninth Annual Conference on Neural Information Processing Systems, 2025. URL https://openreview. net/forum?id=rem4dgVrFg. [40] E. L. Lehmann and J. P. Romano. Testing statistical hypotheses. Springer Texts in Statistics. Springer, New York, third edition, 2005. ISBN 0-387-98864-5. [41] X. Li, B. Su, C. Wang, Q. Long, and W. J. Su. Mitigating privacy-utility trade-off in decentralized federated learning via f -differential privacy. arXiv preprint arXiv:2510.19934, 2025. [42] F. McSherry. Privacy integrated queries: an extensible platform for privacy-preserving data analysis. Commun. ACM, 53(9):89–97, 2010. [43] D. Micciancio and O. Regev. Worst-case to average-case reductions based on Gaussian measures. SIAM Journal on Computing, 37(1):267–302, 2007. [44] I. Mironov. Rényi differential privacy. In 2017 IEEE 30th computer security foundations symposium (CSF), pages 263–275. IEEE, 2017. [45] S. Sachdeva and N. K. Vishnoi. Faster algorithms via approximation theory. Foundations and Trends® in Theoretical Computer Science, 9(2):125–210, 2014. ISSN 1551-305X. doi: 10.1561/0400000065. URL http://dx.doi.org/10.1561/0400000065. [46] I. Shevtsova. Moment-type estimates with asymptotically optimal structure for the accuracy of the normal approximation. Ann. Math. Inform, 39:241–307, 2012. [47] J. Smith, H. J. Asghar, G. Gioiosa, S. Mrabet, S. Gaspers, and P. Tyler. Making the most of parallel composition in differential privacy. Proc. Priv. Enhancing Technol., 2022(1):253–273, 2022. [48] B. Su, W. J. Su, and C. Wang. The 2020 US Decennial Census is more private than you (might) think. Proceedings of the National Academy of Sciences, 122(45):e2500337122, 2025. doi: 10.1073/pnas. 2500337122. URL https://www.pnas.org/doi/abs/10.1073/pnas.2500337122. [49] W. J. Su. A Statistical Viewpoint on Differential Privacy: Hypothesis Testing, Representation, and Blackwell’s theorem. Annual Review of Statistics and Its Application, 12(1):157–175, 2025. [50] T. A. Sullivan. Coming to Our Census: How Social Statistics Underpin Our Democracy (and Republic). Harvard Data Science Review, 2(1), 2020. [51] The American Redistricting Project. Phillips v. U.S. Census Bureau — thearp.org. https://thearp. org/litigation/phillips-v-us-census-bureau/. [52] L. N. Trefethen and J. A. C. Weideman. The exponentially convergent trapezoidal rule. SIAM Review, 56(3):385–458, 2014. doi: 10.1137/130932132. URL https://doi.org/10.1137/130932132. [53] US Census Bureau. 2020 Census Demographic and Housing Characteristics File (DHC) — census.gov. https://www.census.gov/data/tables/2023/dec/2020-census-dhc.html, 2020. [54] US Census Bureau. 2020 census apportionment results. https://www.census.gov/data/tables/2020/ dec/2020-apportionment-data.html, 2021. [55] US Census Bureau. Guidance for labor force statistics data users. https://www.census.gov/topics/ employment/labor-force/guidance.html, 2021.

29

[56] US Census Bureau. New Demonstration Data Will Feature Higher Privacy-loss Budget, Satisfies Redistricting Accuracy Targets — census.gov. https://www.census.gov/programs-surveys/ decennial-census/decade/2020/planning-management/process/disclosure-avoidance/ 2020-das-updates/2021-04-19.html, 2021. [57] US Census Bureau. Privacy-loss budget allocation 2022-08-25. https://www2.census. gov/programs-surveys/decennial/2020/program-management/data-product-planning/ 2010-demonstration-data-products/02-Demographic_and_Housing_Characteristics/ 2022-08-25_Summary_File/2022-08-25_Privacy-Loss_Budget_Allocations.pdf, 2022. [58] US Census Bureau. DAS-implementation-details. U.S. Census Bureau, 2023. https://github.com/uscensusbureau/DAS_2020_DHC_Production_Code/blob/main/wiki/ DAS-Implementation-Details.md.

URL

[59] US Census Bureau. Census bureau data guide more than $2.8 trillion in federal funding in fiscal year 2021, 2023. URL https://www.census.gov/newsroom/press-releases/2023/ decennial-census-federal-funds-distribution.html. [60] US Census Bureau. Announcing the 2030 Census Disclosure Avoidance Research Program — census.gov. https://www.census.gov/newsroom/blogs/research-matters/2025/ 2030-census-disclosure-avoidance.html, 2025. [61] US Census Bureau. Decennial Census of Population and Housing Disclosure Avoidance — census.gov. https://www.census.gov/programs-surveys/decennial-census/disclosure-avoidance. html, 2025. [62] D. T. S. Vu. Numerical resolution of algebraic systems with complementarity conditions. Application to the thermodynamics of compositional multiphase mixtures. Theses, Université Paris-Saclay, Oct. 2020. URL https://ifp.hal.science/tel-02965421. [63] C. Wang, B. Su, J. Ye, R. Shokri, and W. J. Su. Unified enhancement of privacy bounds for mixture mechanisms via f -differential privacy. Advances in Neural Information Processing Systems, 36, 2024. [64] H. Wang, S. Gao, H. Zhang, M. Shen, and W. J. Su. Analytical composition of differential privacy via the Edgeworth accountant. arXiv preprint arXiv:2206.04236, 2022. [65] Y.-X. Wang, B. Balle, and S. P. Kasiviswanathan. Subsampled rényi differential privacy and analytical moments accountant. In The 22nd international conference on artificial intelligence and statistics, pages 1226–1235. PMLR, 2019. [66] Y. Zhu and Y.-X. Wang. Poission subsampled rényi differential privacy. In International Conference on Machine Learning, pages 7634–7642. PMLR, 2019. [67] Y. Zhu and Y.-X. Wang. Improving sparse vector technique with renyi differential privacy. Advances in neural information processing systems, 33:20249–20258, 2020. [68] Y. Zhu, J. Dong, and Y.-X. Wang. Optimal accounting of differential privacy via characteristic function. In International Conference on Artificial Intelligence and Statistics, pages 4782–4817. PMLR, 2022.

30

Supplementary Materials for A Sieve-Accelerated Quadrature Method for Exact Privacy Accounting in the 2020 U.S. Decennial Census A

Preliminaries

A.1

Useful facts for discrete Gaussian distributions

We introduce some useful facts for discrete Gaussian distributions in this section, including the sub-Gaussian tail bound for discrete Gaussian distributions and properties of the characteristic functions. The proof of this section can be found in [13, 48]. Properties of discrete Gaussian distributions. The following proposition on the moment generating function and tail bound shows that the discrete Gaussian is sub-Gaussian. Proposition A.1 (Lemma 16 & Corollary 17 in [13]). A discrete Gaussian random variable X ∼ NZ (0, σ 2 ) is sub-Gaussian with variance proxy σ 2 . For any t ∈ R, we have 2 2

EetX ≤ et σ /2 . Moreover, we have the tail bound 2

2

P(X > t) ≤ e−t /2σ . The following proposition records the numerical values of λ4 defined in Section 6 for some examples of noise parameters σ 2 . Fact A.2. The λ4 defined in Section 6 takes the following values. For σ 2 = 0.1, λ4 ≈ 72.21. For σ 2 = 0.2, λ4 ≈ 4.12. For σ 2 = 0.4, λ4 ≈ 0.19. For σ 2 = 1.0, λ4 ≈ 8.34 × 10−6 . For σ 2 = 2.0, λ4 ≈ 8.92 × 10−14 . For σ 2 = 4.0, λ4 ≈ 2.55 × 10−30 . For σ 2 = 10.0, λ4 < 10−50 . Characteristic functions of discrete Gaussian distribution. Sm defined in Section 6 can be represented as follows. fSm (t) = EeitSm = (a)

=

The characteristic function of

! −u2 /2σ 2 ei·t/Bm ·u m e u=−∞ P∞ −u2 /2σ 2 u=−∞ e ! P∞ −σ 2 (t/Bm −2πu)2 /2 m e u=−∞ P∞ , −2σ 2 π 2 u2 u=−∞ e

P∞

(A.1)

where Equality (a) holds due to Poisson summation formula. For X ∼ NZ (0, σ 2 ), characteristic functions fX (ζ) have the following properties. Lemma A.3 (Proposition A.6 in [48]). We have the following observation on fX (ζ): 31

1. |fX (ζ)| ≤ 1, fX (ζ) achieves its maximum at ζ = 0 with a maximum value of fX (ζ) = 1. 2. fX (ζ) is periodic with period 2π, therefore fX (ai Lζ) is periodic with period 2π/ai L. 3. fX (ζ) is strictly increasing on (−π, 0) and is strictly decreasing on (0, π).

B

Proof of Section 3

B.1

Proof of Proposition 3.3

The probability mass function of Pi and Qi is given by the following. m Y i=1

dPi (x1 , · · · , xm ) =

m Y i=1

2

m Y

2 2 m Y e−(xi −µ) /2σi . dQi (x1 , · · · , xm ) = P −x2 /2σi2 x∈Z e i=1 i=1

2

e−xi /2σi

−x2 /2σi2 x∈Z e

P

,

Given that µ = 1 is an integer, we have m X

2 2 2  m 2 m Y Y dQi (Yi ) e−(Yi −µ) /2σi e−Yi /2σi log = log P P −x2 /2σi2 −x2 /2σi2 dPi (Yi ) x∈Z e x∈Z e i=1 i=1 i=1 ! 2 2  2 2 m X e−(Yi −µ) /2σi e−Yi /2σi = log P P −x2 /2σi2 −x2 /2σi2 x∈Z e x∈Z e i=1

m X Yi

!

µ2

2 − 2σ 2 σ i i i=1 m X Xi µ2 d = + , σ2 2σi2 i=1 i

=

and m X i=1

m

X Xi µ2 dQi (Xi ) − . = log dPi (Xi ) σi2 2σi2 i=1

By Lemma 3.2, we have δ(εkl ) = P

m X Xi

σ2 i=1 i

> εkl −

m X µ2

2σi2 i=1

! − eεkl · P

32

m X Xi i=1

σi2

> εkl +

m X µ2 i=1

2σi2

! .

Let ai = ρi /ρ. Recall that µ = 1, and Equation (2.8), we have ! ! m m m m 2 X X X εkl X µ2 ε µ kl δ(εkl ) = P ai Xi > − − eεkl · P + ai Xi > 2 ρ ρ 2ρσ 2ρσi2 i i=1 i=1 i=1 i=1 ! ! m m m m X X εkl X ai εkl X ai εkl −e ·P − ai Xi > + =P ai Xi > ρ 2 ρ 2 i=1 i=1 i=1 i=1 ! ! m m X X εkl εkl εkl =P ai Xi > −1 −e ·P ai Xi > +1 . ρ ρ i=1

i=1

fk , M fl] The last equation follows from the fact that we are considering the composed composition [M 2 ) denotes parameters of noises injected in M f k and M f l determined by Equation and (σ12 , σ22 , · · · , σm P P (2.8). Therefore, i ρi = 2ρ and i ai = 2.

B.2

Proof of Proposition 3.4

By Neyman-Pearson Lemma, we have αkl (ζ) = PXi ∼NZ (0,σ2 ) i

βkl (ζ) = PYi ∼NZ (µ,σ2 ) i

! ! m X dQi (Xi ) dQi (Xi ) > ζ + c · PXi ∼NZ (0,σ2 ) log =ζ , log i dPi (Xi ) dPi (Xi ) i=1 i=1 ! ! m m X X dQi (Yi ) dQi (Yi ) log ≤ ζ − c · PYi ∼NZ (µ,σ2 ) =ζ log i dPi (Yi ) dPi (Yi ) m X

i=1

i=1

for some c. Equivalently, after change of variable, as µ = 1 we have the parametric equation ! ! m m m m X X X X ai ai αkl (ζ) = PXi ∼NZ (0,σ2 ) ai Xi > ζ + + c · PXi ∼NZ (0,σ2 ) ai Xi = ζ + , i i 2 2 i=1 i=1 i=1 i=1 ! ! m m m m X X X X ai ai βkl (ζ) = PXi ∼NZ (0,σ2 ) ai Xi ≤ ζ − − c · PXi ∼NZ (0,σ2 ) ai Xi = ζ − . i i 2 2 i=1

As

i=1

i=1

i=1

P

i ai = 2, the trade-off function is determined by the following parametric equation:

αkl (ζ) = PXi ∼NZ (0,σ2 )

m X

i

βkl (ζ) = PXi ∼NZ (0,σ2 )

i=1 m X

i

i=1

! ai Xi > ζ + 1

+ c · PXi ∼NZ (0,σ2 ) i

! ai Xi ≤ ζ − 1

33

m X

− c · PXi ∼NZ (0,σ2 )

i=1 m X

i

i=1

! ai Xi = ζ + 1 , ! ai Xi = ζ − 1 .

C

Proof of Section 4

C.1

Proof of Proposition 4.1

Proof. By Lemma 16 in [13], we have   2 2 EXi ∼NZ (0,σ2 ) etXi ≤ et σi /2 . i

Therefore, m h Pm i Y Pm 2 2 2   2 EXi ∼NZ (0,σ2 ) et i=1 ai LXi = EXi ∼NZ (0,σ2 ) etai LXi ≤ et · i=1 ai L σi /2 , i

i

i=1

which further implies that PXi ∼NZ (0,σ2 )

m X

i

! 2

≤ 2 · e−U /2

ai LXi > U

i=1

Pm

2 2 2 i=1 ai L σi

.

For U satisfying Equation (4.1), by Corollary 17 in [13], we have ! m Pm 2 2 2 X 2 PXi ∼NZ (0,σ2 ) ai LXi > U ≤ 2 · e−U /2 i=1 ai L σi ≤ ∆/4. i

i=1

This completes the proof of Proposition 4.1 that "m # X PXi ∼NZ (0,σ2 ) ai Xi > t0 − PXi ∼NZ (0,σ2 ) i

i

i=1

U≥

m X i=1

! ai LXi > t0 · L

≤ ∆/4

for any U satisfying Equation (4.1).

C.2

Proof of Proposition 4.4

Proof. Consider Lemma 4.3 with the complex strip boundary a = 1/(dL) for an optimized parameter d > 0. Let ζ = x + iy for x, y ∈ R restricted to the complex strip −1/(dL) ≤ y ≤ 1/(dL). We aim to rigorously establish an upper bound M such that |F (ζ)| ≤ M across this strip. Recall the definition of F (ζ) from Equation (4.2):   U 1  X F (ζ) = cos(ζt) · E[eiζX ]. 2π t=⌈t0 ·L⌉

By the triangle inequality, we can bound the magnitude of the product by the product of the magnitudes:   U X 1  |F (ζ)| ≤ | cos(ζt)| · |E[eiζX ]| (C.1) 2π t=⌈t0 ·L⌉

34

Step 1: Bounding the Expectation. We must bound the magnitude of the complex expectation without improperly factoring dependent variables. Using Jensen’s inequality, we have |E[ei(x+iy)X ]| = |E[eixX e−yX ]| ≤ E[|eixX e−yX |] = E[|eixX |e−yX ] = E[e−yX ]. P Because the discrete Gaussian convolution X = m i=1 ai LXi is symmetrically distributed around −yX yX 0, we have E[e ] = E[e ]. The moment generating function is an even, convex function of y, so its maximum on the interval y ∈ [−1/(dL), 1/(dL)] is achieved at the boundary y = 1/(dL). Thus, uniformly across the strip: |E[eiζX ]| ≤ E[eX/(dL) ]. P By definition of X = m i=1 ai LXi , we have EeX/dL =

m Y

e(ai /d)·Xi =

i=1

m Y

(

)

xi =−∞

i=1 m Y

∞ X

a2i σi2 /2d2

e(ai /d)·xi · e   ∞  X

−x2i /2σi2

−

√xi 2 2σ

·

m Y

(

∞ X

)−1 e

−x2i /2σi2

i=1 xi =−∞ !2  (  ai σi m  Y √ − d 2

∞ X

)−1 −x2i /2σi2

i e · e     xi =−∞ i=1 xi =−∞ i=1  !2  ( ∞ )−1   x m ∞ m  X  Y − √i X (c) Y 2 2 2 2 2 2 2σ i < eai σi /2d · e · e−xi /2σi   xi =−∞  i=1 xi =−∞ i=1

=

<

m Y

e

2 2

·

2

Pm

eai σi /2d = e

2 2 2 i=1 ai σi /2d

,

i=1

where Equation (c) follows from Fact A.4 in [48]. Step 2: Bounding the Finite Geometric Sum. Next, we bound the sum of the cosine terms. For ζ = x + iy: | cos(ζt)| =

e−yt + eyt ei(x+iy)t + e−i(x+iy)t ≤ = cosh(yt) ≤ e|y||t| ≤ e|t|/(dL) . 2 2

This forms a finite geometric series with common ratio r = e1/(dL) > 1. Its exact value is: U X t=⌈t0 ·L⌉

| cos(ζt)| ≤ 2

U  X t=0

e1/(dL)

t

=2·

e(U +1)/(dL) − 1 e(U +1)/(dL) < 2 · . e1/(dL) − 1 e1/(dL) − 1

Using the fundamental exponential inequality ez − 1 > z for all z > 0, we substitute z = 1/(dL) to 1 < dL. Thus, the summation is strictly bounded by: get e1/(dL) −1 U X t=⌈t0 ·L⌉

| cos(ζt)| < 2 · dL · e(U +1)/(dL) .

35

Step 3: Establishing M and the Trapezoidal Error Bound. Substituting the rigorous bounds from Step 1 and 2 back into Equation (C.1), the absolute upper bound M for |F (ζ)| on the complex strip is: Pm 2 2   Pm 2 2    a σ dL 1  U +1 (U +1)/(dL) i=1 ai σi = M= 2 · dL · e exp exp + i=1 2 i i . 2 2π 2d π dL 2d 4πM By Lemma 4.3, the approximation error of the trapezoidal rule is bounded by eN/(dL) . To guar−1 antee this error is at most ∆/4, we require:

4πM eN/(dL) − 1

≤

∆ 16πM ⇐⇒ eN/(dL) ≥ 1 + , 4 ∆

yielding the sufficient condition: e

N/(dL)

16πM 16dL ≥1+ =1+ exp ∆ ∆



U +1 + dL

Pm

2 2 i=1 ai σi . 2d2

Taking the natural logarithm of both sides and multiplying by dL, we obtain the explicit requirement for N : P   2 2 L m 16dL i=1 ai σi N ≥U +2+ + dL log . 2d ∆ We expand the logarithm: N ≥U +2+

L

Pm

2 2 i=1 ai σi

2d

 − dL log

∆ 8

 + dL log(2 · dL).

(C.2)

Step 4: Optimizing d via AM-GM. To minimize the required number of integration nodes N , we find the optimal parameter d > 0. The threshold in Equation (C.2) is dominated by the two terms containing d. By the AM-GM inequality, their sum is globally minimized when the terms are perfectly equal: qP m Pm 2 2   2 2 i=1 ai σi L i=1 ai σi ∆ = −dL log ⇐⇒ d = p . 2d 8 −2 log(∆/8) Substituting this optimal value of d back into the dominant terms evaluates to exactly twice their balanced value: v ! u      m X u ∆ ∆ 2 2 t 2 −dL log = L −2 ai σi log ≤ U. 8 8 i=1

Substituting this combined minimum back into Equation (C.2) yields the final, formally rigorous condition for N : v ! u   m X u ∆ 2 2 t N ≥ 2 · (U + 1) + dL log(2 · dL) ≥ U + 2 + L −2 ai σi log + dL log(2 · dL). 8 i=1

36

This explicitly corrects the required parameters to guarantee a numerical tolerance of ∆/4, completing the proof.

C.3

Proof of Proposition 4.5

Proof. By triangle inequality, it suffices to prove for each k ∈ [N ] − C, we have 2π · |F (ζk )| ≤ ∆/4. Recall that Equation (4.2) that   U m X Y 1   cos(ζt) · F (ζ) = · fai LXi (ζ). 2π i=1

t=⌈t0 ·L⌉

According to Equation (4.5), we have for each k ∈ [N ] − C, we have |fai LXi (ζk )| < ∆/(8 · U ) for some i. This implies that for each k ∈ [N ] − C, we have   U m X Y   2π · |F (ζk )| = cos(ζt) · fai LXi (ζ) i=1

t=⌈t0 ·L⌉

 ≤ 

U X

t=⌈t0 ·L⌉

 ≤ 

U X

t=⌈t0 ·L⌉

 cos(ζt) · min {|fX1 (ai Lζ)|, · · · , |fXm (ai Lζ)|}  cos(ζt) · ∆/(8 · U ) ≤ ∆/4.

By triangle inequality, this completes the proof of Proposition 4.5.

D

Proof of Theorem 6.1

p P P I.I.D. Recall the notation Xi ∼ NZ (0, σ 2 ) and Sm = B1m m var( m i=1 Xi , where Bm = i=1 Xi ). th κr = κr (Xi ) is the abbreviation for the r cumulant of Xi and λr and Kr as: λr = κr /var(Xi )r/2 ,

Kr = E|Xi |r /var(Xi )r/2 .

By Exercise 3.3.2 (iii) in [20], the Fourier inversion for discrete random variables is given by 1 P(Sm = y) = 2πBm

Z πBm

e−iζy fSm (ζ)dζ.

(D.1)

−πBm

By Fourier Inversion, the edgeworth expansion can be written as     Z ∞ 2 λ4 1 λ4 ζ 4 1 e−y /2 2 4 −iζy −ζ 2 /2 · √ 1+ (3 − 6y + y ) = e e 1+ dζ Bm 24m 2πBm −∞ 24m 2π 37

(D.2)

Combining (D.1) and (D.2), we have Z ∞

  Z πBm λ4 ζ 4 2 e−iζy fSm (ζ)dζ e−iζy e−ζ /2 1 + dζ − 24m −πBm −∞   Z πBm Z πBm 4 λ4 ζ 1 2 e−iζy e−ζ /2 1 + e−iζy fSm (ζ)dζ ≤ dζ − 2πBm −πBm 24m −πBm   Z ∞ 4 1 λ ζ 2 4 + e−iζy e−ζ /2 1 + dζ πBm πBm 24m

Rm ≤

1 2πBm

(D.3)

=: Λ1 + Λ2 Lemma D.1 (Upper bound on Λ2 ). By triangle inequality, we have 1 Λ2 ≤ πBm

Z ∞ e

−ζ 2 /2

πBm

  λ4 ζ 4 1+ dζ 24m

It remains to provide an upper bound on Λ1 .   Z πBm λ4 ζ 4 dζ − e−iζy fSm (ζ)dζ e e 1+ 24m −πBm −πBm   Z πBm 4 λ4 ζ 1 2 e−ζ /2 1 + − fSm (ζ) dζ ≤ 2πBm −πBm 24m   Z πBm λ4 ζ 4 1 −ζ 2 /2 e 1+ ≤ − fSm (ζ) dζ πBm 0 24m

1 Λ1 ≤ 2πBm

Z πBm

−iζy −ζ 2 /2

where the last line follows from the fact that fSm is even function. This further implies that   Z c·√log(m) 1 λ4 ζ 4 −ζ 2 /2 − fSm (ζ) dζ Λ1 ≤ e 1+ πBm 0 24m   Z πBm 1 λ4 ζ 4 (D.4) −ζ 2 /2 + − fSm (ζ) dζ e 1+ √ πBm c· log(m) 24m = Λ 3 + Λ4 for c > 0 satisfying the condition in Lemma D.2. Obtaining a sharp upper bound of Λ3 and Λ4 is technically involved. We summarize them into the following two lemmas and the proof is deferred to Section D.1 and D.2. Lemma D.2 (Upper bound on Λ3 ). For any c that U m (c) < 1, we have  2   Z c·√log(m) ζ 1 ′ ′ ′ Λ3 ≤ exp − × Λ5 + Λ6 + Λ7 dζ πBm 0 2

38

where     ζ4 1 K6 ζ 6 1 K6 ζ 6 K8 ζ 8 K10 ζ 10 K12 ζ 12 · +1 + + · + + + 8m (1 − U m (c))2 6! m2 (1 − U m (c))2 48m2 640m3 17280m4 1036800m5  2 1 K4 ζ 4 K6 ζ 6 Um (ζ) ′ Λ6 = + + 2 4! m 6! m2 2(1 − U m (c))2 !  p 3 4 K log(m)2 6 K log(m)3 U (c log(m)) K4 ζ 4 K6 ζ 6 Um (ζ) 1 c c m 4 6 ′ + × + Λ7 = exp + + 6 4! m 6! m2 2(1 − U m (c))2 4! m 6! m2 2(1 − U m (c))2 Λ′5 =

and Um (ζ), U m (c) are given in equation (D.10) and (D.6). Now, we remain to give an upper bound on Λ4 : Z πBm

 λ4 ζ 4 1+ e − fSm (ζ) dζ √ 24m c· log(m)   Z πBm Z πBm 1 λ4 ζ 4 1 −ζ 2 /2 ≤ e 1+ + |fS (ζ)|dζ πBm c·√log(m) 24m πBm c·√log(m) m

1 Λ4 = πBm

−ζ 2 /2



Lemma D.3 (Upper bound on Λ4 ). Let θ0 be unique root of the Equation (6.3) within the interval (0, 2π). We have 1 Λ4 ≤ πBm

Z πBm c·

√

−ζ 2 /2

e log(m)

  λ4 ζ 4 1+ 24m

cos(θ0 ) − 1 + θ02 /2 c2 + exp − log(m) + · 8c3 K3 · 2 2θ03



log3 (m) m

1/2 !

Combining Lemma D.1, D.2, D.3 completes the proof of Theorem 6.1.

D.1

Proof of Lemma D.2

Proof. The characteristic function fXj is given by EXj2 ζ 2 EXj4 ζ 4 EXj6 ζ 6 fXj (ζ/Bm ) = 1 − + − + ··· 2 4 6 2! Bm 4! Bm 6! Bm Applying a Taylor-Lagrange expansion, there exists a complex number θ1,j,m satisfying |θ1,j,m | < 1 such that Uj,m (ζ) := fXj (ζ/Bm ) − 1 = −

EXj2 ζ 2 EXj4 ζ 4 θ1,j,m EXj6 ζ 6 + − 2 4 6 2! Bm 4! Bm 6! Bm

39

The triangle inequality gives the following upper bound on |Uj,m (ζ)| as follows: EXj4 ζ 4 EXj6 ζ 6 EXj2 ζ 2 + + 2 4 6 2! Bm 4! Bm 6! Bm h p i Denote U m (c) to be the upper bound of |Uj,m (ζ)| for all ζ ∈ 0, c log(m) .

(D.5)

|Uj,m (ζ)| ≤

log(m) K4 U (c) := c · + c4 · 2m 4! m

2



log(m) m

2

K6 +c · 6! 6



log(m) m

3 > |Uj,m (ζ)|

(D.6)

We choose c so that the right hand side of Equation (D.5) is smaller than 1. Specifically, we define c satisfying U m (c) < 1. This ensures existence of a complex number θ2,j,m such that |θ2,j,m | < 1 and log(fXj (ζ/Bm )) = log(1 + Uj,m (ζ)) = Uj,m (ζ) − Summing of all j, we obtain  m X  fSm (ζ) = exp Uj,m (ζ) − j=1



Uj,m (ζ)2 2(1 + θ2,j,m (ζ)Uj,m (ζ))2



(ζ)2

Uj,m  2(1 + θ2,j,m (ζ)Uj,m (ζ))2

 2 EXj4 ζ 4 θ1,j,m EXj6 ζ 6 EXj2 ζ 2 Uj,m (ζ)  = exp  + − − − 2 4 6 2! Bm 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1    2 m 6 6 4 4 2 X θ EX EX Uj,m (ζ) ζ 1,j,m j ζ j ζ  = exp − × exp  − − 4 6 2 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 m X

j=1

Applying Taylor expansion to exp(x), there exists a complex number θ3,m with   m 2 X EXj4 ζ 4 EXj6 ζ 6 |Uj,m (ζ)|  |θ3,m | ≤ sup exp  + + 4 6 2 √ 4! B 6! B 2|1 + θ (ζ)U (ζ)| 2,j,m j,m m m ζ∈[0,c log(m)] j=1   m 2 4 6 X |U (ζ)| K ζ K ζ j,m 4 6  = sup exp  + + 2 2 √ 4! m 6! m 2|1 + θ 2,j,m (ζ)Uj,m (ζ)| ζ∈[0,c log(m)] j=1

40

(D.7)

such that  2 ( m X EXj4 ζ 4 θ1,j,m EXj6 ζ 6 Uj,m (ζ)2 ζ × 1+ fSm (ζ) = exp − − − 4 6 2 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1 2  m θ1,j,m EXj6 ζ 6 Uj,m (ζ)2 1 X EXj4 ζ 4  + − − 4 6 2 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1  3 ) m 4 4 6 6 2 X EX θ EX θ3,m  Uj,m (ζ) 1,j,m j ζ j ζ  + − − 4 6 6 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1  2 ( m θ1,j,m EXj6 ζ 6 Uj,m (ζ)2 ζ λ4 4 X 3κ2 2 ζ 4 = exp − × 1+ ζ + − − 4 6 2 4!m 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1  2 m 6 6 4 4 2 X θ EX EX Uj,m (ζ) 1 1,j,m j ζ j ζ  +  − − 4 6 2 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1  3 ) m 4 6 2 4 6 X EXj ζ θ1,j,m EXj ζ Uj,m (ζ) θ3,m   − − + 4 6 6 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1

Using the triangle inequality, we have   λ4 ζ 4 −ζ 2 /2 1+ fSm − e 24m  2 ( X m θ1,j,m EXj6 ζ 6 Uj,m (ζ)2 ζ 3κ2 2 ζ 4 ≤ exp − × − − 4 6 2 4! Bm 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 j=1  2 m EXj6 ζ 6 |Uj,m (ζ)|2 1 X EXj4 ζ 4  + + + 4 6 2 4! Bm 6! Bm 2|1 + θ2,j,m (ζ)Uj,m (ζ)|2 j=1  3 ) m 4 4 6 6 2 X EX EX |θ3,m |  |Uj,m (ζ)| j ζ j ζ  + + + 4 6 6 4! Bm 6! Bm 2|1 + θ2,j,m (ζ)Uj,m (ζ)|2 j=1  2   ζ × Λ5 + Λ 6 + Λ 7 =: exp − 2

41

(D.8)

In order to give explicit upper bound on Λ5 , Λ6 and Λ7 , (D.5) implies that m X j=1

2

|Uj,m (ζ)| ≤ ≤

m X j=1

EXj4 ζ 4 EXj6 ζ 6 EXj2 ζ 2 + + 2 4 6 2! Bm 4! Bm 6! Bm

m X (E|Xj |2 )2 ζ 4 j=1

+

+

4 4Bm

!2

(E|Xj |4 )2 ζ 8 (E|Xj |6 )2 ζ 12 (E|Xj |2 )(E|Xj |4 )ζ 6 + + 8 12 6 576Bm 518400Bm 24Bm

(E|Xj |2 )(E|Xj |6 )ζ 8 (E|Xj |4 )(E|Xj |6 )ζ 10 + 8 10 720Bm 8640Bm

By Jensen’s inequality, for any k, l ≤ 6, m X

1

k

B k+l j=1 m

l

(E|Xj | )(E|Xj | ) ≤

m X

1

B k+l j=1 m

(E|Xj |k+l )k/k+l (E|Xj |k+l )l/k+l

(D.9)

Kk+l,m = (k+l−2)/2 m Therefore, Equation (D.9) allows us to conclude m X j=1

|Uj,m (ζ)|2 ≤

K6 ζ 6 K8 ζ 8 K8 ζ 8 K10 ζ 10 K12 ζ 12 K4 ζ 4 + + + + + 4m 24m2 720m3 576m3 8640m4 518400m5

K4 ζ 4 K6 ζ 6 K8 ζ 8 K10 ζ 10 K12 ζ 12 + + + + 4m 24m2 320m3 8640m4 518400m5 := Um (ζ) =

42

(D.10)

Equation (D.10) allows us to bound Λ5 , Λ6 and Λ7 as follows: m X 3κ2 2 ζ 4

Λ5 =

j=1

− 4

4! Bm

θ1,j,m EXj6 ζ 6 Uj,m (ζ)2 − 6 6! Bm 2(1 + θ2,j,m (ζ)Uj,m (ζ))2

m

X Uj,m (ζ)2 ζ4 K6 ζ 6 + + 2 8m 2(1 + θ2,j,m (ζ)Uj,m (ζ)) 6! m2 j=1  m X (E|Xj |2 )2 ζ 4 (E|Xj |4 )2 ζ 8 (E|Xj |6 )2 ζ 12 ζ4 1 + × + + ≤ 4 8 12 8m 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 4Bm 576Bm 518400Bm j=1  (E|Xj |2 )(E|Xj |4 )ζ 6 (E|Xj |2 )(E|Xj |6 )ζ 8 (E|Xj |4 )(E|Xj |6 )ζ 10 K6 ζ 6 + + + + 6 8 10 24Bm 720Bm 8640Bm 6! m2  m X (E|Xj |2 )2 ζ 4 (E|Xj |4 )2 ζ 8 (E|Xj |6 )2 ζ 12 ζ4 1 ≤ + + + × 4 8 12 8m 2(1 + θ2,j,m (ζ)Uj,m (ζ))2 4Bm 576Bm 518400Bm j=1  (E|Xj |2 )(E|Xj |4 )ζ 6 (E|Xj |2 )(E|Xj |6 )ζ 8 (E|Xj |4 )(E|Xj |6 )ζ 10 K6 ζ 6 + + + + 6 8 10 24Bm 720Bm 8640Bm 6! m2   m X 1 1  ζ4 ≤ 1 + (1 + θ2,j,m (ζ)Uj,m (ζ))2 m 8m ≤

j=1

+

m X j=1

1 · 2(1 + θ2,j,m (ζ)Uj,m (ζ))2

(E|Xj |4 )2 ζ 8 (E|Xj |6 )2 ζ 12 + 8 12 576Bm 518400Bm

(E|Xj |2 )(E|Xj |4 )ζ 6 (E|Xj |2 )(E|Xj |6 )ζ 8 (E|Xj |4 )(E|Xj |6 )ζ 10 + + + 6 8 10 24Bm 720Bm 8640Bm

! +

K6 ζ 6 6! m2

+

K6 ζ 6 (D.11) 6! m2

Recall Equation (D.6) and the fact that |θ2,j,m | < 1, we conclude  Λ5 (ζ) ≤ +

1 +1 (1 − U m (c))2 m 1X 1 2

j=1



ζ4 8m

· (1 − U m (c))2

(E|Xj |4 )2 ζ 8 (E|Xj |6 )2 ζ 12 + 8 12 576Bm 518400Bm

(E|Xj |2 )(E|Xj |4 )ζ 6 (E|Xj |2 )(E|Xj |6 )ζ 8 (E|Xj |4 )(E|Xj |6 )ζ 10 + + + 6 8 10 24Bm 720Bm 8640Bm   ζ4 1 K6 ζ 6 ≤ · + 1 + 8m (1 − U m (c))2 6! m2   1 K6 ζ 6 K8 ζ 8 K10 ζ 10 K12 ζ 12 + · + + + (1 − U m (c))2 48m2 640m3 17280m4 1036800m5

43

!

Similarly, we also give explicit upper bound of Λ6 and Λ7 .  2 m EXj6 ζ 6 |Uj,m (ζ)|2 1 X EXj4 ζ 4  Λ6 = + + 4 6 2 4! Bm 6! Bm 2|1 + θ2,j,m (ζ)Uj,m (ζ)|2 j=1  2 m 4 6 X 1 K4 ζ K6 ζ 1 ≤  + + |Uj,m (ζ)|2  2 4! m 6! m2 2(1 − U m (c))2

(D.12)

j=1

≤

1 2



ζ4

ζ6

K4 K6 Um (ζ) + + 2 4! m 6! m 2(1 − U m (c))2

2

with Um (ζ) defined in (D.10). To have an upper bound for Λ7 , we first need an upper bound for |θ3,m |. Recall Equation (D.7) that   m 2 4 6 X |U (ζ)| ζ K ζ K j,m 6 4  |θ3,m | ≤ sup exp  + + 2 √ 4! m 6! m 2|1 + θ (ζ)Uj,m (ζ)|2 2,j,m ζ∈[0,c log(m)] j=1 ! √ c4 K4 log(m)2 c6 K6 log(m)3 supζ∈[0,c log(m)] Um (ζ) + + ≤ exp 4! m 6! m2 2(1 − U m (c))2 ! p c4 K4 log(m)2 c6 K6 log(m)3 Um (c log(m)) + + = exp 4! m 6! m2 2(1 − U m (c))2 Having the upper bound on |θ3,m |, Λ7 yields the following upper bound. 3  m EXj6 ζ 6 |θ3,m | X EXj4 ζ 4 |Uj,m (ζ)|2  Λ7 = + + 4 6 6 4! Bm 6! Bm 2|1 + θ2,j,m (ζ)Uj,m (ζ)|2 j=1 ! p 1 c4 K4 log(m)2 c6 K6 log(m)3 Um (c log(m)) ≤ exp + + 6 4! m 6! m2 2(1 − U m (c))2  3 m 4 6 2 4 6 X EXj ζ EXj ζ |Uj,m (ζ)|  × + + 4 6 4! Bm 6! Bm 2|1 + θ2,j,m (ζ)Uj,m (ζ)|2 j=1 ! p 1 c4 K4 log(m)2 c6 K6 log(m)3 Um (c log(m)) ≤ exp + + 6 4! m 6! m2 2(1 − U m (c))2  3 K4 ζ 4 K6 ζ 6 Um (ζ) × + + 4! m 6! m2 2(1 − U m (c))2 Combining (D.11), (D.12) and (D.13), we conclude that fSm − e

−ζ 2 /2

   2   ζ λ4 ζ 4 ′ ′ ′ ≤ exp − 1+ × Λ5 + Λ 6 + Λ 7 24m 2

44

(D.13)

where     ζ4 1 K6 ζ 6 1 K6 ζ 6 K8 ζ 8 K10 ζ 10 K12 ζ 12 · +1 + + · + + + 8m (1 − U m (c))2 6! m2 (1 − U m (c))2 48m2 640m3 17280m4 1036800m5  2 1 K4 ζ 4 K6 ζ 6 Um (ζ) ′ Λ6 = + + 2 4! m 6! m2 2(1 − U m (c))2 !  p 3 4 K log(m)2 6 K log(m)3 U (c log(m)) K4 ζ 4 K6 ζ 6 Um (ζ) 1 c c m 4 6 ′ + × + Λ7 = exp + + 6 4! m 6! m2 2(1 − U m (c))2 4! m 6! m2 2(1 − U m (c))2 Λ′5 =

This completes the proof of this lemma.

D.2

Proof of Lemma D.3

Proof. We would like to give an upper bound on |fSm (t)|. The idea is inspired by Theorem 2.2 in [46]. Lemma D.4 ([46]). Let θ0 be unique root of the Equation (6.3) within the interval (0, 2π). For any x ∈ R, and θ0 < θ ≤ 2π, then cos(x) ≤ 1 − a(θ)x2 + b(θ)|x|3 where 1 − cos(θ) sin(θ) − θ2 θ 1 − cos(θ) sin(θ) b(θ) = 2 · − θ3 θ2

a(θ) = 3 ·

ei be an independent copy of Xi , then we have: Now, let X 2

|fSm (ζ)| =

m Y i=1

 fXi

ζ Bm

2 =

m Y

E cos

i=1

ei ) ζ(Xi − X Bm

ei )2 = 2κ2 , we obtain By Lemma D.4 and relation E(Xi − X m Y

e i )2 e i |3 ζ 2 E(Xi − X |ζ|3 E|Xi − X |fSm (ζ)| ≤ 1 − a(θ) + b(θ) 2 3 Bm Bm i=1 ! m Y ei |3 2ζ 2 κ2 |ζ|3 E|Xi − X ≤ 1 − a(θ) 2 + b(θ) 3 Bm Bm 2

i=1

45

!

By AM-GM inequality, we have !m m e i |3 1 X 2ζ 2 κ2 |ζ|3 E|Xi − X 1− a(θ) 2 + b(θ) 3 m Bm Bm i=1 !!m Pm e 3 2 2 3 i=1 E|Xi − Xi | 1+ −a(θ)ζ + b(θ)|ζ| 3 m 2Bm

|fSm (ζ)|2 ≤ ≤ This implies that

|fSm (ζ)|2 ≤ exp −2a(θ)ζ 2 + b(θ)|ζ|3

Pm

e 3 i=1 E|Xi − Xi | 3 Bm

!

For convenience, we choose a specific θ = θ0 . We have upper bound cos(θ0 ) − 1 + θ02 /2 · |fSm (ζ)|2 ≤ exp −ζ 2 + θ03

Pm

e 3 i=1 E|Xi − Xi | · |ζ|3 3 Bm

By Lemma 4.6, the characteristic function is strictly decreasing between c p Therefore, for any ζ ∈ [c log(m), πBm ], we have p c2 cos(θ0 ) − 1 + θ02 /2 |fSm (ζ)| ≤ |fSm (c log(m))| ≤ exp − log(m) + · 2 2θ03

!

p log(m) and πBm .

Pm

e 3 i=1 E|Xi − Xi | · c3 (log(m))3/2 3 Bm

This further implies that 1 Λ4 ≤ πBm

Z πBm c·

√

e

−ζ 2 /2

log(m)

  λ4 ζ 4 1+ 24m

c2 cos(θ0 ) − 1 + θ02 /2 + exp − log(m) + · 2 2θ03   Z πBm λ4 ζ 4 1 −ζ 2 /2 = e 1+ πBm c·√log(m) 24m

Pm

e 3 i=1 E|Xi − Xi | · c3 (log(m))3/2 3 Bm

c2 cos(θ0 ) − 1 + θ02 /2 8mE|Xi |3 3 + exp − log(m) + · 3/2 · c (log(m))3/2 2 2θ03 κ2 m3/2   Z πBm 1 λ4 ζ 4 −ζ 2 /2 e 1+ = πBm c·√log(m) 24m  3 1/2 ! c2 cos(θ0 ) − 1 + θ02 /2 log (m) + exp − log(m) + · 8c3 K3 · 2 m 2θ03 This gives an upper bound on Λ4 .

46

!

!

!

E

Supplementary Figures and Tables

The supplemental figures include 946 figures with corresponding raw data tables, which is the prifk , M f l ]. All figures follow the same pattern as vacy accounting of all the composed mechanisms [M in Figures 3 and 4. For brevity, we provide the Figures and Tables in the folder “results” in GitHub repository (see https://github.com/BuxinSu/Exact-Privacy-Accounting-for-2020-U.S.-Census. git). We also include all privacy budget allocations ρi used in the 2020 Census DHC File in the folder “privacy budget allocation”.

47