# Integer Subspace Differential Privacy

Prathamesh Dharangutte Rutgers University<EMAIL_ADDRESS>Jie Gao Rutgers
University<EMAIL_ADDRESS>Ruobin Gong Rutgers University,
<EMAIL_ADDRESS>Fang-Yi Yu George Mason University<EMAIL_ADDRESS>
###### Abstract

We propose new differential privacy solutions for when external _invariants_
and _integer_ constraints are simultaneously enforced on the data product.
These requirements arise in real world applications of private data curation,
including the public release of the 2020 U.S. Decennial Census. They pose a
great challenge to the production of provably private data products with
adequate statistical usability. We propose _integer subspace differential
privacy_ to rigorously articulate the privacy guarantee when data products
maintain both the invariants and integer characteristics, and demonstrate the
composition and post-processing properties of our proposal. To address the
challenge of sampling from a potentially highly restricted discrete space, we
devise a pair of unbiased additive mechanisms, the generalized Laplace and the
generalized Gaussian mechanisms, by solving the Diophantine equations as
defined by the constraints. The proposed mechanisms have good accuracy, with
errors exhibiting sub-exponential and sub-Gaussian tail probabilities
respectively. To implement our proposal, we design an MCMC algorithm and
supply empirical convergence assessment using estimated upper bounds on the
total variation distance via $L$-lag coupling. We demonstrate the efficacy of
our proposal with applications to a synthetic problem with intersecting
invariants, a sensitive contingency table with known margins, and the 2010
Census county-level demonstration data with mandated fixed state population
totals.

##  1 Introduction

Motivation. Differential privacy (DP) is a formal mathematical framework that
quantifies the extent to which an adversary can learn about an individual from
sanitized data products. However, data curators may be mandated, by law or by
policy, to release sanitized data products that obey certain externally
determined constraints, in a manner that existing differential privacy
solutions are not designed to address. Two challenging types of constraints
are _invariants_ (Ashmead et al.,, 2019), which are exact statistics
calculated from the confidential data (such as sum, margins of a contingency
table, etc), and _integer_ characteristics, pertaining to data of count
nature. The challenge is exacerbated when the sanitized data product is
required to respect both.

Data privatization that simultaneously satisfies mandated invariants and
preserves integral characteristics is of paramount importance to a myriad of
data products by official statistical agencies. One of the prominent example
is the Disclosure Avoidance System (DAS) of the 2020 Decennial Census of the
United States. The Census Bureau is mandated to observe a set of invariants
for its decennial census data products, such as its constitutional obligation
to enumerate state populations for the apportionment house seats. As the
decennial census data products are count data, the bureau must ensure that the
sanitized release are also integer-valued. In addition, the bureau strives to
maintain the consistency of the data products with common knowledge to various
degrees, such as the the non-negativity of count data and possible logical
relationships between different tabulated quantities. The challenge thus
remains: how to design sanitized data products that respect the pre-specified
constraints, while maintaining rigorous privacy and good statistical
properties? This is the central problem we address in this paper.

Our contribution. In this paper, we develop the _integer subspace differential
privacy_ scheme, through which we formulate and implement mathematically
rigorous privacy mechanisms that meet the following _utility objectives_ : 1)
respecting mandated invariants, 2) maintaining integral characteristics, and
3) achieving good statistical properties including _unbiasedness_ , _accuracy_
, as well as _probabilistic transparency_ of the mechanism specification.
Invariants not only restrict the universe of databases of interest, but also
partially destroys the relevance of neighbors, since databases differing by
precisely one entry may or may not meet the same invariants. To this end, we
adapt the classical differential privacy definition to incorporate invariant
constraints with two modifications (i) we limit the comparison of databases
only to those that satisfy the same invariants, in a manner similar to
pufferfish privacy (Kifer and Machanavajjhala,, 2014); and (ii) for databases
that meet the same invariants, we use their metric distance as a
multiplicative factor to the privacy loss parameter under the more general
smooth differential privacy framework (Chatzikokolakis et al.,, 2013;
Desfontaines and Pejó,, 2019). We show that integer subspace differential
privacy preserves _composition_ and _post-processing_ properties.

To design a differential privacy mechanism that simultaneously respects
invariants and integer requirement poses a number of new challenges. For
additive mechanisms, the perturbation noise vector needs to be limited to the
“null space” of the invariants so as not to violate the constraints. When the
invariants are linear, this requires solving a linear system with integer
coefficients, i.e. Diophantine equations, limiting the privacy noise to a
(transformed) discrete lattice space. To this end, we propose generalized
discrete Laplace and Gaussian mechanisms for the discrete lattice space. Both
mechanisms are transparently specified additive mechanisms that enjoy
demonstrated unbiasedness and accuracy via bounds on tail probabilities. The
tail bounds call for extra technicality as the lattice spaces we deal with are
generally not spherically symmetric. To implement the proposed mechanisms, we
design a Gibbs-within-Metropolis sampler, with a transition kernel that is
always limited to the desired subspace to achieve sampling efficiency. To
provide empirical guidance on the Markov chain’s convergence to target, we
supply an assessment scheme using the $L$-lag coupling method proposed by
(Biswas et al.,, 2019). We demonstrate the efficacy of our proposal with
applications to a synthetic problem with intersecting invariants, a
contingency table with known margins concerning delinquent children (Federal
Committee on Statistical Methodology,, 2005), and the 2010 Census county-level
demonstration data with mandated fixed state population totals.

Related work. There has been an extensive line of work on differential privacy
and its alternative definitions; see a recent survey and the references
therein (Desfontaines and Pejó,, 2019). On the other hand, invariants pose a
relatively recent challenge to formal privacy, and we focus our review below
on prior work concerning invariants. It has been recognized that invariants
which are non-trivial functions of the confidential data compromise the
differential privacy guarantee in its classic sense, and that the invariants
must be incorporated into the privacy guarantee (Gong and Meng,, 2020; Seeman
et al.,, 2022). Towards designing formally private mechanisms that meet
prescribed invariant constraints, the only prior work in this direction is
named _subspace differential privacy_ (Gao et al.,, 2022) where the inputs
satisfy a set of linear constraints and the sanitized outcome shall meet the
same invariants. However, subspace differential privacy only considers data
products that are real-valued. To impose the additional integral requirement
on the data outputs is not trivial: as the TopDown algorithm (Abowd et al.,,
2022) already demonstrates, simple rounding or projection will either violate
the invariants or introduce biases. He et al., (2014) consider known
constraints of the dataset (e.g. population counts for each state) which
restricts the possible set of datasets, but the output of their mechanism need
not satisfy the known constraints. We focus on mechanisms whose output also
satisfy the counting constraints. Another line of work (Cormode et al.,, 2017)
consider additional structural properties for the mechanism’s output (e.g.
monotone, symmetric, fair). However, those properties are independent of the
private data and can be seen as internal distributional consistency. Lastly we
note a related, but different line of previous work concerning _internal
consistency_ of the data outputs (Barak et al.,, 2007; Hay et al.,, 2009),
such as maintaining the sum of counts over some partitioning of a set to be
identical as the total sum. This type of consistency requirement is
independent of the value of the private data, and can be satisfied by certain
types of DP mechanisms, such as local DP schemes when each data item is
independently perturbed and normal algebraic operations are applied on the
perturbed output. In contrast, invariant constraints considered in this work
are in general non-trivial functions of the confidential data.

To impose invariants and preserve integrality for Census data, the TopDown
algorithm (Abowd et al.,, 2022) employs a discrete Gaussian mechanism (Canonne
et al.,, 2020) to inject integer-valued noise during the so-called
_measurement_ phase, to form noise-infused counts to be tabulated into multi-
way contingency tables at various levels of geographic resolution. This is
followed by the so-called _estimation_ phase, a constrained optimization to
impose invariants (such as state population totals, etc) and controlled
rounding to integer solutions (see Table 1 of Abowd et al.,, 2022). In
general, this approach first imposes unconstrained noise and performs post-
processing to meet the additional constraints. One limitation of this solution
is that it is not easy to characterize the probabilistic distribution of the
resulting privacy noise after post-processing, as this distribution crucially
depends on the particular values of the input data. Thus, the post-processing
approach may destroy the probabilistic transparency of the privatized data
product (Gong,, 2022). Our solution maintains the transparency of the privacy
mechanism as well as the statistical independence between the confidential
data and the privacy errors, which are of paramount importance to downstream
applications derived from these noisy data products. Furthermore, since
invariants and integer value constraints are our first priority besides
privacy, we supply a rigorously updated definition of differential privacy
under this setting by limiting the discussion only among databases that
satisfy these constraints.

##  2 Integer Subspace Differential Privacy

In this work, we model private data as a histogram
${\vec{x}}=\left(x_{1},\dots,x_{d}\right)^{\top}\in\mathbb{N}^{d}$, where
$x_{i}\in\mathbb{N}$ is the number of individuals with feature
$i\in[d]:=\\{1,\dots,d\\}$. A trusted curator holds the database ${\vec{x}}$,
and provides an interface to the database through a randomized mechanism
$M:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ where we require the output to always take
integer values. As motivated in introduction, our goal is to design good
mechanisms whose output is close to the original histogram that satisfies not
only certain privacy notions, but also invariant constraints and integral
requirements.

The notion of differential privacy ensures that no individual’s data has much
effect on the probabilistic characteristic of the output from the mechanism
$M$ (Dwork et al.,, 2006; Barber and Duchi,, 2014). Formally, we say that a
randomized mechanism $M:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ is
_$\left(\epsilon,\delta\right)$ -differentially private_ if for some
$\epsilon,\delta\geq 0$, all neighboring databases ${\vec{x}}$ and
${\vec{x}}^{\prime}$ with $\|{\vec{x}}-{\vec{x}}^{\prime}\|=1$, and all event
$S\subseteq\mathbb{Z}^{d}$,

$\Pr\left[M\left({\vec{x}}\right)\in S\right]\leq
e^{\epsilon}\cdot\Pr\left[M\left({\vec{x}}^{\prime}\right)\in
S\right]+\delta,$ (1)

where $\|{\vec{x}}-{\vec{x}}^{\prime}\|$ is a distance norm of two databases
$\left({\vec{x}},{\vec{x}}^{\prime}\right)$. Depending on the application, the
norm may be chosen either as $\ell_{1}$ norm or $\ell_{2}$ norm.

The notion of differential privacy in Eq. 1 satisfies a more general notion of
_smooth differential privacy_ (Chatzikokolakis et al.,, 2013). $M$ is
$(\epsilon,\delta)$-_smooth differentially private_ if for all datasets
${\vec{x}},{\vec{x}}^{\prime}$ and all event $S$ ,

$\Pr[M({\vec{x}})\in S]\leq
e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M({\vec{x}}^{\prime})\in
S]+\frac{e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon}-1}\delta.$
(2)

Previous works mostly consider the pure smooth differential privacy with
$\delta=0$, and we generalize it to the approximated case.111When $\delta=0$,
Eq. 2 reduces to $\Pr[M({\vec{x}})\in E]\leq
e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}\cdot\Pr[M({\vec{x}}^{\prime})\in
E]$. However, for the approximated smooth differential privacy, the error term
grows as the distance $\|{\vec{x}}-{\vec{x}}^{\prime}\|$ increases. For
instance, as the distance increases, we can see the error term of the Gaussian
mechanism also grows. Finally, when $\|\cdot\|$ is a norm, a mechanism on the
histogram is $(\epsilon,\delta)$-differentially private if and only if the
mechanism is $(\epsilon,\delta)$-smooth differentially private. Compared to
Eq. 1, the definition in Eq. 2 replaces the neighboring relationship of
$\left({\vec{x}},{\vec{x}}^{\prime}\right)$ by a norm function.

From the data curator’s perspective, in addition to privacy concerns, there
often exists external constraints that the privatized output $M$ must meet.
These constraints can often be represented as counting functions defined
below.

######  Definition 1 (name = counting invariant constraints, label =
def:count_inv).

Given a collection of subsets $\mathcal{A}=\\{A_{1},\dots,A_{k}\\}$ on $[d]$
with $k\leq d$, we say a function $f:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ is
$\mathcal{A}$-invariant if

$\sum_{i\in A}x_{i}=\sum_{i\in A}f({\vec{x}})_{i}\text{, for all
}{\vec{x}}\in\mathbb{N}^{n}\text{ and }A\in\mathcal{A}.$

Alternatively, given a collection of counting invariant constraints
$\mathcal{A}$, we define an equivalence relationship $\equiv_{\mathcal{A}}$ so
that ${\vec{x}}\equiv_{\mathcal{A}}{\vec{x}}^{\prime}$ if $\sum_{i\in
A}x_{i}=\sum_{i\in A}x_{i}^{\prime}$ for all $A\in\mathcal{A}$. Then an
$\mathcal{A}$-invariant $f$ satisfies
${\vec{x}}\equiv_{\mathcal{A}}f({\vec{x}})$ for all ${\vec{x}}$.

Note that when a privatized output is subject to counting invariant
constraints as above, we may not use the standard neighboring relation because
two neighboring datasets differing by one element may not satisfy the same
counting constraint. In particular, the feasible outputs of an
$\mathcal{A}$-invariant mechanism lie within an integer subspace. Therefore,
we employ smooth differential privacy, and only control the privacy loss on
datasets that satisfy the same counting constraints, i.e., only “secret pairs”
of databases within the same equivalent classes. This formulation was also
used in previous work such as pufferfish privacy (Kifer and Machanavajjhala,,
2014).

######  Definition 2.

We say a mechanism $M$ is $(\epsilon,\delta)$-differentially private on an
equivalence relation $\equiv$ if for all equivalent pair
${\vec{x}}\equiv{\vec{x}}^{\prime}$ and event $E$ on the output space

$\Pr[M({\vec{x}})\in E]\leq
e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M({\vec{x}}^{\prime})\in
E]+\frac{e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon}-1}\delta.$

Moreover, given a collection of counting constraint $\mathcal{A}$, we say $M$
is _$\mathcal{A}$ -induced integer subspace differentially private_ with
_privacy loss budget_ $(\epsilon,\delta)$, if $M$ is $\mathcal{A}$-invariant
and $(\epsilon,\delta)$-differentially private on an equivalence relation
$\equiv_{\mathcal{A}}$.

Many invariant constraints can be formulated using counting constraints. Below
are some examples.

######  Example 3.

Suppose $\mathcal{A}=\\{[d]\\}$ is a singleton. All datasets with the same
number of agents are equivalent under $\equiv_{\\{[d]\\}}$, so the
differential privacy on $\equiv_{\\{[d]\\}}$ reduces to the original (bounded)
differential privacy.

######  Example 4.

Privatized census data products must ensure that the total population of each
state is reported exactly as enumerated. Given $k$ states, and
$A_{l}\subset[d]$ be the collection of features that belong to state $l$. we
can define counting invariants $\mathcal{A}=\\{A_{1},\dots,A_{k}\\}$ where
$A_{1},\dots,A_{k}$ forms a partition on the universe $[d]$.

######  Example 5.

We can encode the invariant condition on a $\sqrt{d}\times\sqrt{d}$ two-
dimensional contingent table where the sum of each row and each column sums
are fixed as a collection of counting constraints. Formally, given
$\sqrt{d}\in\mathbb{N}$, the set $\mathcal{A}$ contains the following subsets
on $[d]$, for all $i$ and $j$
$R_{i}=\\{i\sqrt{d}+\ell:0\leq\ell<\sqrt{d}\\}\text{ and
}C_{j}=\\{j+\ell\sqrt{d}:0\leq\ell<\sqrt{d}\\}.$ That is
$\mathcal{A}=\\{R_{i},C_{j}:0\leq i,j<\sqrt{d}\\}$.

Composition and post-processing properties for Definition 2 also hold due to
similar arguments as the differential privacy, but with the additional
requirement of equivalence over databases. Let
$\equiv_{(\mathcal{A}_{1},\mathcal{A}_{2})}$ be the equivalence relation such
that $x\equiv_{(\mathcal{A}_{1},\mathcal{A}_{2})}x^{\prime}$ implies
$x\equiv_{\mathcal{A}_{1}}x^{\prime}$ and
$x\equiv_{\mathcal{A}_{2}}x^{\prime}$.

######  Proposition 6 (Composition).

Let $M_{1}$ be $(\epsilon_{1},\delta_{1})$\- differentially private on
equivalence relation $\equiv_{\mathcal{A}_{1}}$ and $M_{2}$ be
$(\epsilon_{2},\delta_{2})$\- differentially private on equivalence relation
$\equiv_{\mathcal{A}_{2}}$. For any pair of databases $x,x^{\prime}$ having
$x\equiv_{(\mathcal{A}_{1},\mathcal{A}_{2})}x^{\prime}$, the composed
mechanism $M_{1,2}(x)=(M_{1}(x),M_{2}(x))$ is
$(\epsilon_{1}+\epsilon_{2},\delta_{1}+\delta_{2})$-differentially private on
$\equiv_{\mathcal{A}_{1},\mathcal{A}_{2}}$.

######  Proposition 7 (Post-processing).

If $M$ is $(\epsilon,\delta)$-differentially private over $\equiv$, then for a
arbitrary randomized mapping $F:\mathcal{Y}\to\mathcal{Z}$, $F\circ M$ is
$(\epsilon,\delta)$-differentially private over $\equiv$.

Details of Proposition 6 and 7 can be found in Appendix A. Note that when
privacy mechanisms that obey different invariant constraints are imposed, it
might be possible to infer aspects of the confidential data following logical
consequences from both. It is not clear how to define post-processing for
invariant constraints in those situations.

##  3 Generalized Laplace and Gaussian Mechanisms

Now we study $\mathcal{A}$-induced integer differentially private mechanisms
that satisfy counting invariant constraints. The main challenge is that
counting invariant constraints significantly restrict the feasible output
space, especially when the output is required to take only integer values. To
resolve this issue, we first show that counting invariant constraints on
integer-valued output can be written as a lattice space. Then we propose
revised Laplace and Gaussian mechanisms on a lattice space. We discuss how to
implement these mechanisms via sampling and MCMC in the Implementation
section.

###  3.1 Counting invariants and lattice spaces

In this section, we show the output space of integer datasets satisfying
counting invariant constraints are lattice spaces. A _lattice_ $\Lambda$ is a
discrete additive subgroup of $\mathbb{R}^{d}$. Given a matrix
${\bm{B}}\in\mathbb{R}^{d\times m}$ consisting of $m$ basis
${\vec{b}}_{1},\dots,{\vec{b}}_{m}$ with $1\leq m\leq d$, a lattice generated
by the basis is
$\Lambda({\bm{B}})=\left\\{{\bm{B}}{\vec{v}}=\sum_{i=1}^{m}v_{i}{\vec{b}}_{i}:{\vec{v}}\in\mathbb{Z}^{m}\right\\}.$

######  Proposition 8.

Given a collection of counting invariants
$\mathcal{A}=\\{A_{1},\dots,A_{k}\\}$, there exists a lattice
$\Lambda_{\mathcal{A}}=\Lambda({\bm{C}}_{\mathcal{A}})$with basis
${\bm{C}}_{\mathcal{A}}\in\mathbb{Z}^{d\times(d-k)}$ so that a function
$f:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ is $\mathcal{A}$-invariant if and only if
for all ${\vec{x}}\in\mathbb{N}^{d}$

$f({\vec{x}})-{\vec{x}}\in\Lambda_{\mathcal{A}}.$

Note that as ${\bm{C}}_{\mathcal{A}}$ is an integer-valued matrix,
$\Lambda_{\mathcal{A}}$ is a subset of the integer grid $\mathbb{Z}^{d}$. As
we impose more counting invariants, the lattice $\Lambda_{\mathcal{A}}$
becomes sparser, making it harder to find a feasible solution for an
$\mathcal{A}$-invariant mechanism.

###### Proof.

Given a collection of counting invariant constraints
$\mathcal{A}=\\{A_{1},\dots,A_{k}\\}$ and a dataset ${\vec{x}}$, finding a
feasible output ${\vec{y}}$ with ${\vec{y}}\equiv_{\mathcal{A}}{\vec{x}}$ is
equivalent to solving the following linear equations in Eq. 3, where
${\bm{A}}\in\\{0,1\\}^{k\times d}$ is an _incidence matrix_ , for all
$i\in[d]$ and $l\leq k$ $A_{l,i}=1$ if $i\in A_{l}$ and zero otherwise. Then
${\vec{y}}\equiv_{\mathcal{A}}{\vec{x}}$ if and only if

${\bm{A}}{\vec{y}}={\bm{A}}{\vec{x}}.$ (3)

Since we require ${\vec{y}}\in\mathbb{Z}^{d}$, the problem of solving
${\vec{y}}$ is known as solving the linear Diophantine equation (Gilbert and
Pathria,, 1990; Schrijver,, 1998; Greenberg,, 1971), and can be done by
computing the _Smith normal form_ of ${\bm{A}}$. Specifically, given an
integer-valued matrix ${\bm{A}}\in\mathbb{Z}^{k\times d}$ of rank $k$, there
exists ${\bm{U}}\in\mathbb{Z}^{k\times k},{\bm{V}}\in\mathbb{Z}^{d\times d}$
and ${\bm{D}}\in\mathbb{Z}^{k\times d}$ with

${\bm{U}}{\bm{A}}{\bm{V}}={\bm{D}}\quad$ (4)

so that ${\bm{U}}$ and ${\bm{V}}$ are unimodular matrices that are invertible
over the integers, and ${\bm{D}}$ is a diagonal matrix (i.e. the Smith normal
form of ${\bm{A}}$) with $D_{l,l}\neq 0$ for all $l\in[k]$.

With the above decomposition, we can characterize the integer solutions of Eq.
3. ${\vec{y}}$ is a solution of Eq. 3, if and only if
${\bm{A}}({\vec{y}}-{\vec{x}})={\mathbf{0}}$ and it is equivalent to
${\bm{D}}{\bm{V}}^{-1}({\vec{y}}-{\vec{x}})={\mathbf{0}}$, since ${\bm{U}}$ is
unimodular. Because ${\bm{D}}$ is diagonal and has rank $k$, ${\vec{y}}$ is a
solution if and only if there exists ${\vec{w}}\in\mathbb{Z}^{d}$ so that
${\vec{y}}={\vec{x}}+{\bm{V}}{\vec{w}}$ and $w_{l}=0$ for all $l\leq k$.
Additionally, if we define
${\bm{C}}_{\mathcal{A}}\in\mathbb{Z}^{d\times(d-k)}$ that consists the bottom
$d-k$ rows of ${\bm{V}}$, then ${\vec{y}}$ is a solution if and only if
${\vec{y}}-{\vec{x}}\in\\{{\bm{C}}_{\mathcal{A}}{\vec{v}}:{\vec{v}}\in\mathbb{Z}^{d-k}\\}=\Lambda({\bm{C}}_{\mathcal{A}})$.
We call $\mathcal{A}$ _full rank_ if the rank of
${\bm{A}}\in\\{0,1\\}^{k\times d}$ is $k$, the associated
${\bm{C}}_{\mathcal{A}}$ as the _basis of $\mathcal{A}$_, and
$\Lambda_{\mathcal{A}}:=\Lambda({\bm{C}}_{\mathcal{A}})$. The integer
solutions of Eq. 3 is the lattice $\Lambda_{\mathcal{A}}$ shifted by
${\vec{x}}$. Therefore, if $M:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ is
$\mathcal{A}$-invariants the output given input ${\vec{x}}$ is contained in
${\vec{x}}+\Lambda_{\mathcal{A}}:=\\{{\vec{x}}+{\vec{z}}:{\vec{z}}\in\Lambda_{\mathcal{A}}\\}$,
$M({\vec{x}})\in{\vec{x}}+\Lambda_{\mathcal{A}}.$ ∎

###  3.2 Generalized Laplace mechanism

Now we can define a generalized Laplace mechanism whose output space satisfies
Proposition 8. We adopt the classical exponential mechanism to this restricted
output space and show the resulting mechanisms are as required. Because the
output space Proposition 8 is unbounded, we show the utility guarantee by
controlling the dimension of the lattice spaces which may be of interest by
itself. First, we formally define the generalized Laplace mechanism.

######  Definition 9 (label = def:lap).

Given a collection of counting invariant condition $\mathcal{A}$ of full rank
and $\epsilon>0$, the _generalized Laplace mechanism_ is
$M_{Lap,\mathcal{A},\epsilon}({\vec{x}})={\vec{x}}+{\vec{z}}$ with ${\vec{z}}$
sampled from

$q_{\epsilon}({\vec{z}})\propto\exp\left(-\epsilon\|{\vec{z}}\|\right)\mathbf{1}[{\vec{z}}\in\Lambda_{\mathcal{A}}]$
(5)

where $\Lambda_{\mathcal{A}}$ is defined in Proposition 8.

Note that if the counting constraint is _vacuous_ , i.e.
$\mathcal{A}=\emptyset$ with $d=1$, and the output can be real number, the
generalized Laplace mechanism reduces to the original Laplace mechanism on
histograms. Additionally, if we require integer-valued output with
$\mathcal{A}=\emptyset$, the generalized Laplace mechanism becomes the double
geometric mechanism.

######  Theorem 10 (label = thm:lap).

Given full rank $\mathcal{A}$ and $\epsilon>0$, mechanism
$M_{Lap,\mathcal{A},\epsilon}:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ in
LABEL:def:lap is $\mathcal{A}$-induced integer subspace with $(\epsilon,0)$.

The proof is similar to the privacy guarantee of Laplace mechanisms, and is
presented in Appendix B. The implementation of the generalized Laplace
mechanisms in LABEL:def:lap is non-trivial and is elaborated in Implementation
section.

After the privacy guarantee, we turn to discuss the utility of the generalized
Laplace mechanism. First, Eq. 6 establishes the unbiasedness of the
generalized Laplace mechanism. The tail bound needs additional work, because
the output spaces of the conventional exponential mechanism are often finite,
which is not the case here. We resolve this issue using the bounded-dimension
property of lattices spaces (LABEL:lem:doubling), and prove the error bound is
sub-exponential in Eq. 7. There is a long line of work (Landau,, 1915; Tsuji,,
1953; Bentkus and Gotze,, 1999) that approximate the number of lattice points
in sphere by the volume. Using these ideas, we prove the following lemma (with
proof in Appendix B) for LABEL:thm:lap_acc.

######  Lemma 11 (label = lem:doubling).

Let ${\bm{B}}\in\mathbb{R}^{d\times m}$ be a rank $m\leq d$ matrix. For all
$r>0$ large enough, $|\\{{\vec{z}}\in\Lambda({\bm{B}}):\|{\vec{z}}\|_{2}\leq
r\\}|\leq\frac{2V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}r^{m}$ where
$V_{m}$ is the volume of the $m$ dimensional unit sphere.

######  Theorem 12 (name = unbiasedness and accuracy,label = thm:lap_acc).

Given $\mathcal{A}$ of $k$ counting constraints and $\epsilon>0$,
$M_{Lap,\mathcal{A},\epsilon}$ in LABEL:def:lap is unbiased

$\mathbb{E}[M_{Lap,\mathcal{A},\epsilon}({\vec{x}})]={\vec{x}},$ (6)

and, if $\|\cdot\|$ is $\ell_{2}$-norm, there exists a constant $K>0$ so that
for all ${\vec{x}}\in\mathbb{N}^{d}$, and large enough $t>0$,

$\Pr[\|M_{Lap,\mathcal{A},\epsilon}({\vec{x}})-{\vec{x}}\|_{2}\geq t]\leq
Kt^{d-k}\exp(-\epsilon t).$ (7)

Moreover, $K$ can be $\frac{4\cdot
2^{d-k}V_{d-k}}{\sqrt{\det({\bm{C}}_{\mathcal{A}}^{\top}{\bm{C}}_{\mathcal{A}})}}$
where $V_{d-k}$ is the volume of the $(d-k)$ dimensional unit sphere.

Note that as the dimension $d-k$ increases $2^{d-k}V_{d-k}\to 0$, so $K$ is
determined by the counting constraints
$1/\sqrt{\det({\bm{C}}_{\mathcal{A}}^{\top}{\bm{C}}_{\mathcal{A}})}$.

###  3.3 Generalized Gaussian mechanism

In this section, we propose generalized Gaussian mechanisms that relate to
Lattice-based cryptography discussed in Section 4. We first define Gaussian
random variables on lattices. Given a lattice $\Lambda\subset\mathbb{R}^{d}$,
and ${\vec{c}}\in\mathbb{R}^{d}$, the spherical _Gaussian random variable_
${\vec{Z}}$ on $\Lambda$ with variance $\sigma$ and center ${\vec{c}}$ is

$\Pr[{\vec{Z}}={\vec{z}}]\propto\exp\left(-\frac{1}{2\sigma^{2}}\|{\vec{z}}-{\vec{c}}\|_{2}^{2}\right),\forall{\vec{z}}\in\Lambda.$
(8)

Note that the generalized Gaussian mechanism utilizes an $\ell_{2}$ norm to
measure the noise scale.

######  Definition 13 (label = def:gaussian).

Given a collection of counting invariant condition $\mathcal{A}$ of full rank,
$\epsilon>0$, and $\delta>0$, the _generalized Gaussian mechanism_ is defined
as $M_{G,\mathcal{A},\epsilon,\delta}({\vec{x}})={\vec{x}}+{\vec{z}}$ where
${\vec{z}}$ is the Gaussian random variable on $\Lambda_{\mathcal{A}}$ with
center at ${\mathbf{0}}$ and variance
$\sigma_{\epsilon,\delta}^{2}=\frac{2c_{\mathcal{A}}\ln
1/\delta}{\epsilon^{2}}$ for some constant
$c_{\mathcal{A}}=\mathcal{O}(\max\\{(d-k)\ln(d-k),\ln K\\})$ that only depends
on the dimension and the collection of counting constraints where $K$ is
defined in LABEL:thm:lap_acc.

While the variance scale in the original Gaussian mechanisms is independent of
the dimension, generalized Gaussian mechanism’s $\sigma_{\epsilon,\delta}$ in
LABEL:def:gaussian depends on the dimension. This may be due to a lack of
symmetry of the lattice space. Recall that the proof of the original Gaussian
mechanism reduces the high dimensional case to the one dimensional setting,
thanks to the Gaussian being spherically symmetric. However, our lattice space
$\Lambda_{\mathcal{A}}$ is generally not spherically symmetric, and simple
reduction does not work.

######  Theorem 14 (label = thm:gaussian).

Given full rank $\mathcal{A}$, $\epsilon$ and $\delta$ with
$0<\delta<\epsilon<1/e$, mechanism
$M_{G,\mathcal{A},\epsilon,\delta}:\mathbb{N}^{d}\to\mathbb{Z}^{d}$ in
LABEL:def:gaussian is $\mathcal{A}$-induced integer subspace differentially
private with $(\epsilon,\delta)$.

After the privacy guarantee, we turn to show the utility of the generalized
Gaussian mechanism. First, since the distribution function of the discrete
Gaussian is an even function, the additive errors of the discrete Gaussian
mechanism is unbiased Eq. 9. Furthermore, similar to Gaussians on the
Euclidean space, Gaussians on lattices are sub-Gaussian.222This is nontrivial,
because the support of Gaussian in Eq. 8 is not uniform. For instance, suppose
$\Lambda$ is not a lattice space, but has $2^{r^{2}}$ points in each integer
radius $r$ shell, the a Gaussian on such a space is not sub-Gaussian.

######  Theorem 15 (name = unbiasedness and accuracy).

Given the condition in LABEL:thm:gaussian, $M_{G,\mathcal{A},\epsilon,\delta}$
in LABEL:def:gaussian is unbiased,

$\mathbb{E}[M_{G,\mathcal{A},\epsilon,\delta}({\vec{x}})]={\vec{x}}\text{, for
all }{\vec{x}}\in\mathbb{N}^{d}.$ (9)

Additionally, there exists a constant $K>0$ defined in LABEL:thm:lap_acc so
that for all ${\vec{x}}\in\mathbb{N}^{d}$, and large enough $t>0$,

$\Pr[\|M_{G,\mathcal{A},\epsilon,\delta}({\vec{x}})-{\vec{x}}\|_{2}\geq t]\leq
Kt^{d-k}e^{-\frac{t^{2}}{2\sigma_{\epsilon,\delta}^{2}}}.$ (10)

Last, we remark that previous work has used discrete Gaussian mechanisms
(Canonne et al.,, 2020)for integer valued noises. But these mechanisms do not
satisfy (non-trivial) invariants constraints. Note that when there is no
invariant constraint, the resulting lattice space becomes an integer grid and
our generalized Gaussian mechanism reduces to the discrete Gaussian mechanism.

##  4 Implementation

We discuss implementation issues of discrete Gaussian and generalized Laplace
in lattice spaces. In general the design of exact samplers for complex privacy
mechanisms (such as the Exponential mechanism, to which our generalized
Laplace mechanism is a special case) is a widely acknowledged challenge in the
literature. Section 3 of Abowd et al. (Abowd et al.,, 2022) specifically
discussed this challenge in the context of invariants and integer constraints
and concluded that direct sampling from the Exponential mechanism is
“infeasible” for their TopDown algorithm. We address this problem in two ways.
For certain discrete Gaussian, we can use exact samplers; for generalized
Laplace we develop an efficient and practical solution using Markov chain
Monte Carlo (MCMC). Our experiments all use the MCMC sampler for practicality.

Sampling discrete Gaussian on a lattice. There are several efficient
algorithms to sample discrete variables within negligible statistical distance
of any discrete Gaussian distribution whose scales are not exceedingly narrow
(Gentry et al.,, 2008; Peikert,, 2010). Moreover, (Brakerski et al.,, 2013)
provide an exact Gaussian sampler.

######  Theorem 16 (name = Lemma 2.3 in Brakerski et al., (2013)).

There is a probabilistic polynomial-time algorithm that takes as input a basis
${\bm{B}}=({\vec{b}}_{1},\dots,{\vec{b}}_{m})$ for a lattice
$\Lambda({\bm{B}})\subset\mathbb{R}^{d}$, a center
${\vec{c}}\in\mathbb{R}^{d}$ and parameter
$\sigma=\Omega(\max_{i}\|{\vec{b}}_{i}\|_{2}\ln d)$ and outputs a vector that
is distributed exactly as Eq. 8.

We note that, however, it is believed to be computationally difficult to
sample efficiently on general lattice spaces if the target distribution is not
sufficiently dispersed. Most lattice-based cryptographic schemes (Folláth,,
2014) are based on the hardness assumption of solving the short integer
solution (SIS) problem and the learning with errors (LWE) problem (Brakerski
et al.,, 2013), and the ability of an efficient sampler of a narrow
distribution on general lattice efficiently would enable us to solve short
integer solution and break lattice-based cryptographic schemes.

MCMC sampling for generalized Laplace on a lattice. Since a discrete Gaussian
mechanism only provides $({\epsilon},\delta)$-DP, it is important to develop
pure-DP mechanisms using the generalized Laplace mechanism on a lattice. We
devise an MCMC sampling scheme to instantiate the generalized Laplace
mechanism (LABEL:def:lap). As a practical solution, MCMC has been
theoretically studied in the literature; e.g. Wang et al., (2015) for
stochastic gradient Monte Carlo and Ganesh and Talwar, (2020) for discretized
Langevin MCMC. The challenge to designing such a sampling scheme is twofold.
First, the invariants render the target state space highly constrained. To
achieve sampling efficiency requires efficient proposals, which ideally
satisfy the invariants themselves in order to maintain a reasonable acceptance
rate (c.f. Gong and Meng,, 2020). Second, MCMC algorithms are generally not
exact samplers, and may fail to converge to the target distribution in finite
time, thus incur additional cost to privacy that can be difficult to quantify
for a given application. To this end, we 1) target only the additive noise
component of the mechanism, ensuring that any privacy leakage due to sampling
is independent of the underlying confidential data; and 2) devise an empirical
convergence assessment for the upper bound on the total variation (TV)
distance between the chain’s marginal distribution the target, based on the
$L$-lag coupling method proposed by (Biswas et al.,, 2019).

Proposal design. We wish to sample the additive noise ${\vec{z}}$ employed by
the generalized Laplace mechanism. For an incidence matrix
${\bm{A}}\in\mathbb{Z}^{k\times d}$ specifying the invariants, recall its
Smith normal form ${\bm{D}}={\bm{U}}{\bm{A}}{\bm{V}}$ in Eq. 4. Consider a
jumping distribution with probability mass function
$g\left({\vec{u}}\right)=g_{{\bm{A}}}\left({\bm{V}}^{-1}{\vec{u}}\right),$
where

$g_{{\bm{A}}}\left({\vec{e}}\right)=\prod_{j=1}^{k}{\bf
1}\left\\{e_{j}=0\right\\}\prod_{j=k+1}^{d}\eta_{j}\left(e_{j}\right),$ (11)

for ${\vec{e}}=\left(e_{1},\ldots,e_{d}\right)$ a _pre-jump_ object whose
first $k$ entries are exactly zero, and ${\vec{u}}={\bm{V}}{\vec{e}}$ the
proposed jump. The $\eta_{j}$’s in Eq. 11 can be set to any interger-valued
univariate probabilities that are unbiased and symmetric around zero. A
straightforward choice is the double geometric distribution, with mass
function $\eta_{j}\left(e\right)=\frac{1-a}{1+a}a^{\left\|e\right\|_{1}}$ and
parameter $a$. Note that the jumping distribution $g$ connects to $g_{\bm{A}}$
in Eq. 11 without a Jacobian because ${\vec{u}}$ and ${\vec{e}}$ are integer-
valued and ${\bm{V}}$ is invertible. The proposed jump ${\vec{u}}\sim g$ is
unbiased and symmetric, i.e. $\mathbb{E}_{g}({\vec{u}})={\bf 0}$ and
$g({\vec{u}})=g(-{\vec{u}})$, and always respects the desired invariants
specified by ${\bm{A}}$:

${\bm{A}}{\vec{u}}={\bm{U}}^{-1}{\bm{U}}{\bm{A}}{\bm{V}}{\vec{e}}={\bm{U}}^{-1}{\bm{D}}{\vec{e}}={\bm{U}}^{-1}{\bf
0}={\bf 0}.$

A Metropolis sampling scheme. The target distribution
$q_{\epsilon}\left({\vec{z}}\right)$ is given in Eq. 5 and is known only up to
a normalizing constant. Algorithm 1 in Appendix D presents a Gibbs-within-
Metropolis sampler that produces a sequences of dependent draws
$\left({\vec{z}}^{\left(l\right)}\right)_{0\leq l\leq\text{nsim}}$ from the
target distribution $q_{\epsilon}$ in Eq. 5 known only up to a normalizing
constant. We use an additive jumping distribution whose element-wise
construction is described in Eq. 11. The algorithm incurs a transition kernel
that dictates how the chain moves from an existing state to the next one:
${\vec{z}}^{\left(l\right)}\sim
K\left({\vec{z}}^{\left(l-1\right)},\cdot\right)$. The initial distribution
$\pi_{0}$ may be chosen simply as $g$ for convenience. The choice of
$\eta_{j}$ is a tuning decision for the algorithm, and should be made to
encourage fast mixing of the chain. We discuss our choices for examples in the
Experiments section. Note that steps 5 through 7 of Algorithm 1 updates the
proposed jump in a Gibbs sweep (i.e. one dimension at a time), utilizing the
fact that the pre-jump object ${\vec{e}}$ is element-wise independent under
$g_{\bm{A}}$. Doing so facilitates the $L$-lag coupling, to be discussed next,
for assessing empirical convergence of the chain. In practice, these updates
may be performed simultaneously.

Empirical convergence assessment using $L$-lag coupling. We perform empirical
convergence assessment of the proposed algorithm using $L$-lag coupling.
Construct a joint transition kernel $\tilde{K}$ of two Markov chains, each
having the same target distribution $q_{\epsilon}$ induced by the marginal
transition kernel $K$ as defined in Algorithm 1. The joint kernel $\tilde{K}$,
given by Algorithm 2 in Appendix D, is a maximal coupling kernel such that the
$L$-th lag of the two chains will couple in finite time with probability one.
The random $L$-lag meeting time, which Algorithm 3 in Appendix D samples,
provides an estimate of the upper bound on the TV distance,

$d_{TV}(q_{\epsilon}^{(l)},q_{\epsilon})\leq\mathbb{E}\left[\max\left(0,\left\lceil\left(\tau^{(L)}-L-l\right)/L\right\rceil\right)\right]$

between the target distribution $q_{\epsilon}$ and the marginal distribution
of the chain at time $l$, denoted $q_{\epsilon}^{(l)}$ (Biswas et al.,, 2019,
Theorem 2.5). The upper bounds are obtained as empirical averages over
independent runs of coupled Markov chains. The lag $L>0$ would need to be set,
and we discuss its choice in the next section and Appendix D. Compared to
earlier work which uses MCMC to instantiate privacy mechanisms, our use of
$L$-lag coupling provides real-time (rather than asymptotic) assessment on the
MCMC convergence behavior. The number of iterations needed to ensure
convergence is empirically assessed via the upper bound on the total variation
(note that this is still an estimate). In general, there is a tradeoff between
the number of iterations and an extra privacy loss budget $\delta^{\prime}$,
in the sense that Eq. Eq. 5 can be absorbed as another additive error in the
DP guarantee.

##  5 Experiments

Figure 1: Estimated TV bound from target

County | Low | Medium | High | Very High | Total
---|---|---|---|---|---
Alpha | -1 | -1 | 0 | 2 | 0
Beta | 5 | 4 | -2 | -7 | 0
Gamma | -5 | -3 | 5 | 3 | 0
Delta | 1 | 0 | -3 | 2 | 0
Total | 0 | 0 | 0 | 0 | 0
Table 1: Generalized Laplace additive noise ($\epsilon=0.25$, $\ell_{1}$-norm
target distribution) using Algorithm 1, which preserves row and column
margins. Figure 2: Privacy noise from the generalized Laplace mechanism (
LABEL:def:lap) via Algorithm 1 (blue squares: one instance; boxplot: 1000
instances) for county populations of Illinois in increasing county sizes.
State population total is invariant. The proposed noises are integer-valued
and unbiased. For comparison are DAS errors from the Nov 2020 vintage 2010
Census demonstration data (Van Riper et al.,, 2020) (red dots).

Intersecting counting constraints. Consider a synthetic example in which a set
of three counting constraints $\mathcal{A}=(A_{1},A_{2},A_{3})$ are defined
over 14 records, where the intersection of all subsets of constrains are
nonempty. The constraints are schematically depicted by Fig. 3 in Appendix E,
and may be encoded by an incident matrix ${\bm{A}}\in\mathbb{Z}^{3\times 14}$,
with each row corresponding to one of $A_{1},A_{2}$ and $A_{3}$, and each
element being $1$ at indices corresponding to the record within that
constraint and $0$ otherwise. We apply Algorithm 1 to instantiate the
generalized Laplace mechanism, with both $\ell_{1}$-norm and $\ell_{2}$-norm
targets, to privatize a data product that conforms to the constraint
${\bm{A}}$. To ensure adequate dispersion of the target distribution, we set
$\epsilon=0.25$, a value on the smaller end within the range of meaningful
privacy protection (e.g. Dwork,, 2011). The pre-jump proposal distributions
$\eta_{j}$ are double geometric distributions, with parameter $a=\exp{(-1)}$
for the $\ell_{1}$-norm target and $a=\exp{(-1.5)}$ for the $\ell_{2}$-norm
target. Details on the tuning, the noise distribution and convergence
assessment can be found in Appendix E. In particular, Fig. 4 in Appendix E
shows that the chains empirically converge at around $10^{5}$ iterations for
the $\ell_{1}$-norm target and $10^{6}$ iterations for the $\ell_{2}$-norm
target.

Delinquent children by county and household head education level. The Federal
Committee on Statistical Methodology published a fictitious dataset concerning
delinquent children in the form of a $4\times 4$ contingency table, tabulated
across four counties by education level of household head (Table 4 in Federal
Committee on Statistical Methodology,, 2005, reproduced in Table 2 of Appendix
E), to illustrate various traditional SDL techniques (Slavković and Lee,,
2010).

The charge is to extend privacy protection to the sensitive individual records
while preserving the margins of the contingency table for data publication. To
do so, we apply Algorithm 1 to instantiate the generalized Laplace mechanism
with both $\ell_{1}$\- and $\ell_{2}$-norm targets and with $\epsilon=0.25$.
Fig. 1 shows the evolution of the TV upper bound on the chain’s marginal
distributions to the $\ell_{1}$-norm target, estimated with 200 independent
coupled chains, which appear to converge after about $10^{4}$ iterations and
are stable at various choices of $L$. Table 1 shows one instance of the
proposed additive noise, obtained after the chain achieves empirical
convergence. They are integer-valued, with zero row and column totals which
would preserve the margins of the confidential table. Convergence assessment
for the $\ell_{2}$-norm target and discussions on the choice of proposal are
in Appendix E.

2010 U.S. Census county-level population data. We consider the publication of
county-level population counts subject to the invariant of state population
size, and compare with the privacy-protected demonstration files produced by
preliminary versions of the 2020 Census DAS. The confidential values are the
2010 Census Summary Files (CSF), curated by IPUMS NHGIS and are publicly
available (Van Riper et al.,, 2020). Employed for this example are the
November 2020 vintage demonstration data, protected by pure differential
privacy with $\epsilon=0.192=4\text{ (total)}\times 0.16\text{ (county
level)}\times 0.3\text{ (population query)}$.

We demonstrate Algorithm 2 using the generalized Laplace mechanism under
$\ell_{1}$ norm, with $\epsilon$ set to accord to the Census Bureau’s
specification and $a=\exp(-2.5)$. Fig. 2 showcases the proposed county-level
errors applied to the population of Illinois. The $x$-axis is arranged in
increasing true county population sizes. Blue squares shows an instance of the
proposed noise vector. The boxplots summarize $1000$ proposed noise vectors
(thinned at $0.01\%$), with whiskers indicating $1.5$ times the interquartile
range and cross-marks indicating extreme realizations. The proposed errors are
integer-valued and centered around zero, confirming their marginal
unbiasedness. Note that the published DAS privacy errors (red dots) show a
clear negative bias as a function of increasing underlying county population
size. The bias is likely due to optimization-based post-processing during the
estimation phase which imposes non-negativity on the constituent geographic
areas with small population counts. Similar observations were made for other
states (Appendix E).

##  6 Discussion and Future Work

This paper provides solutions for sanitizing private data that are required to
simultaneously observe pre-specified invariants and integral characteristics.
The proposed _integer subspace differential privacy_ scheme allows for
rigorous statements of privacy guarantees while maintaining good statistical
properties, including unbiasedness, accuracy, and probabilistic transparency
of the privatized output. An efficient MCMC scheme is devised to instantiate
the proposed mechanism, alongside tools to assess empirical convergence.

One direction for future work is to design exact sampling algorithms for the
proposed generalized discrete Laplace and Gaussian mechanisms. In the case
that the normalizing constant to the target distribution $q_{\epsilon}$ is
known exactly, it may be feasible to design, for example, efficient rejection
sampling based on intelligent choices of proposal distributions. It is also
conceivable that a perfect sampling scheme, such as _coupling from the past_
(Propp and Wilson,, 1996), may be designed. The recent work of Seeman et al.,
(2021) proposes an exact sampling scheme for Exponential mechanisms using
_artificial atoms_. However, the technique requires the state space be
compact, which is not the case for the constrained yet unbounded target
distribution $q_{\epsilon}$ considered in this work. Another direction for
future work is to extend invariant-respecting privacy protection data products
that are required to be _binary_ (i.e., exhibiting or not exhibiting a private
attribute), as well as to obey _inequality_ constraints (e.g. non-negativity
of count data). The challenge there is again on the design of efficient
methods to generate noise vectors in an extremely constrained discrete space.

##  7 Acknowledgements

Prathamesh Dharangutte and Jie Gao would like to acknowledge support from NSF
through grants CCF-2118953, IIS-2207440, CCF-2208663 and CNS-2137245. Ruobin
Gong would like to acknowledge support from NSF through grant DMS-1916002.

## References

  * Abowd et al., (2022) Abowd, J. M., Ashmead, R., Cumings-Menon, R., Garfinkel, S., Heineck, M., Heiss, C., Johns, R., Kifer, D., Leclerc, P., Machanavajjhala, A., Moran, B., Sexton, W., Spence, M., and Zhuravlev, P. (2022).  The census topdown algorithm.  Harvard Data Science Review.
  * Ashmead et al., (2019) Ashmead, R., Kifer, D., Leclerc, P., Machanavajjhala, A., and Sexton, W. (2019).  Effective privacy after adjusting for invariants with applications to the 2020 census.  Technical report, US Census Bureau.
  * Barak et al., (2007) Barak, B., Chaudhuri, K., Dwork, C., Kale, S., McSherry, F., and Talwar, K. (2007).  Privacy, accuracy, and consistency too: a holistic solution to contingency table release.  In Proceedings of the twenty-sixth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems, pages 273–282.
  * Barber and Duchi, (2014) Barber, R. F. and Duchi, J. C. (2014).  Privacy and statistical risk: Formalisms and minimax bounds.  arXiv 1412.4451.
  * Bentkus and Gotze, (1999) Bentkus, V. and Gotze, F. (1999).  Lattice point problems and distribution of values of quadratic forms.  Annals of Mathematics, 150(3):977–1027.
  * Biswas et al., (2019) Biswas, N., Jacob, P. E., and Vanetti, P. (2019).  Estimating convergence of markov chains with l-lag couplings.  Advances in Neural Information Processing Systems, 32.
  * Brakerski et al., (2013) Brakerski, Z., Langlois, A., Peikert, C., Regev, O., and Stehlé, D. (2013).  Classical hardness of learning with errors.  page 575–584.
  * Canonne et al., (2020) Canonne, C. L., Kamath, G., and Steinke, T. (2020).  The discrete gaussian for differential privacy.  Advances in Neural Information Processing Systems, 33:15676–15688.
  * Chatzikokolakis et al., (2013) Chatzikokolakis, K., Andrés, M. E., Bordenabe, N. E., and Palamidessi, C. (2013).  Broadening the scope of differential privacy using metrics.  In Cristofaro, E. D. and Wright, M. K., editors, Privacy Enhancing Technologies - 13th International Symposium, PETS 2013, Bloomington, IN, USA, July 10-12, 2013. Proceedings, volume 7981 of Lecture Notes in Computer Science, pages 82–102. Springer.
  * Cormode et al., (2017) Cormode, G., Kulkarni, T., and Srivastava, D. (2017).  Constrained differential privacy for count data.  arXiv preprint arXiv:1710.00608.
  * Desfontaines and Pejó, (2019) Desfontaines, D. and Pejó, B. (2019).  Sok: Differential privacies.  CoRR, abs/1906.01337.
  * Dwork, (2011) Dwork, C. (2011).  A firm foundation for private data analysis.  Communications of the ACM, 54(1):86–95.
  * Dwork et al., (2006) Dwork, C., McSherry, F., Nissim, K., and Smith, A. (2006).  Calibrating noise to sensitivity in private data analysis.  In Theory of cryptography conference, pages 265–284. Springer.
  * Federal Committee on Statistical Methodology, (2005) Federal Committee on Statistical Methodology (2005).  Report on statistical disclosure limitation methodology.  Technical report, Statistical Policy Working Paper 22 (Second Version).  https://www.hhs.gov/sites/default/files/spwp22.pdf [Accessed: 05-15-2022].
  * Folláth, (2014) Folláth, J. (2014).  Gaussian sampling in lattice based cryptography.  Tatra Mountains Mathematical Publications, 60(1):1–23.
  * Ganesh and Talwar, (2020) Ganesh and Talwar (2020).  Faster differentially private samplers via rényi divergence analysis of discretized langevin MCMC.  In Adv. Neural Inf. Process. Syst., pages 7222–7233.
  * Gao et al., (2022) Gao, J., Gong, R., and Yu, F.-Y. (2022).  Subspace differential privacy.  In Proceedings of the The Thirty-Sixth AAAI Conference on Artificial Intelligence (AAAI-22), pages 3986–3995.  arXiv:2108.11527.
  * Gentry et al., (2008) Gentry, C., Peikert, C., and Vaikuntanathan, V. (2008).  Trapdoors for hard lattices and new cryptographic constructions.  In Proceedings of the fortieth annual ACM symposium on Theory of computing, pages 197–206.
  * Gilbert and Pathria, (1990) Gilbert, W. J. and Pathria, A. (1990).  Linear diophantine equations.  preprint.
  * Gong, (2022) Gong, R. (2022).  Transparent privacy is principled privacy.  Harvard Data Science Review (Special Issue 2).  doi:10.1162/99608f92.b5d3faaa.
  * Gong and Meng, (2020) Gong, R. and Meng, X.-L. (2020).  Congenial differential privacy under mandated disclosure.  In Proceedings of the 2020 ACM-IMS on Foundations of Data Science Conference, pages 59–70.
  * Greenberg, (1971) Greenberg, H. (1971).  Integer programming.  Academic Press.
  * Hay et al., (2009) Hay, M., Rastogi, V., Miklau, G., and Suciu, D. (2009).  Boosting the accuracy of differentially-private histograms through consistency.  arXiv preprint arXiv:0904.0942.
  * He et al., (2014) He, X., Machanavajjhala, A., and Ding, B. (2014).  Blowfish privacy: Tuning privacy-utility trade-offs using policies.  In Proceedings of the 2014 ACM SIGMOD international conference on Management of data, pages 1447–1458.
  * Hotz and Salvo, (2020) Hotz, J. V. and Salvo, J. (2020).  Assessing the use of differential privacy for the 2020 Census: Summary of what we learned from the CNSTAT workshop.  Technical report, National Academies Committee on National Statistics (CNSTAT).  https://www.amstat.org/asa/files/pdfs/POL-CNSTAT_CensusDP_WorkshopLessonsLearnedSummary.pdf [Accessed: 04-08-2020].
  * Kifer and Machanavajjhala, (2014) Kifer, D. and Machanavajjhala, A. (2014).  Pufferfish: A framework for mathematical privacy definitions.  ACM Transactions on Database Systems (TODS), 39(1):1–36.
  * Landau, (1915) Landau, E. (1915).  Zur analytischen zahlentheorie der definiten quadratischen formen (uber die gitterpunkte in einem mehrdimensionalen ellipsoid), sitzungsber. koniglich breissischen.  Akad. Wiss., 31:458–476.
  * Peikert, (2010) Peikert, C. (2010).  An efficient and parallel gaussian sampler for lattices.  In Annual Cryptology Conference, pages 80–97. Springer.
  * Propp and Wilson, (1996) Propp, J. G. and Wilson, D. B. (1996).  Exact sampling with coupled markov chains and applications to statistical mechanics.  Random Structures & Algorithms, 9(1-2):223–252.
  * Schrijver, (1998) Schrijver, A. (1998).  Theory of linear and integer programming.  John Wiley & Sons.
  * Seeman et al., (2021) Seeman, J., Reimherr, M., and Slavković, A. (2021).  Exact privacy guarantees for markov chain implementations of the exponential mechanism with artificial atoms.  Advances in Neural Information Processing Systems, 34.
  * Seeman et al., (2022) Seeman, J., Slavkovic, A., and Reimherr, M. (2022).  A formal privacy framework for partially private data.  arXiv preprint arXiv:2204.01102.
  * Slavković and Lee, (2010) Slavković, A. B. and Lee, J. (2010).  Synthetic two-way contingency tables that preserve conditional frequencies.  Statistical Methodology, 7(3):225–239.
  * Tsuji, (1953) Tsuji, M. (1953).  On lattice points in an n-dimensional ellipsoid.  Journal of the Mathematical Society of Japan, 5(3-4):295–306.
  * U.S. Census Bureau, (2019) U.S. Census Bureau (2019).  Memorandum 2019.25: 2010 demonstration data products - design parameters and global privacy-loss budget.  https://www.census.gov/programs-surveys/decennial-census/decade/2020/planning-management/plan/memo-series/2020-memo-2019_25.html [Accessed: 04-18-2022].
  * U.S. Census Bureau, (2021) U.S. Census Bureau (2021).  Privacy-loss budget allocation 2021-06-08.  https://www2.census.gov/programs-surveys/decennial/2020/program-management/data-product-planning/2010-demonstration-data-products/ppmf20210608/2021-06-08-privacy-loss_budgetallocation.pdf [Accessed: 02-15-2022].
  * Van Riper et al., (2020) Van Riper, D., Kugler, T., and Schroeder, J. (2020).  IPUMS NHGIS privacy-protected 2010 Census demonstration data, version 20201116.  Minneapolis, MN: IPUMS.
  * Wang et al., (2015) Wang, Y.-X., Fienberg, S. E., and Smola, A. (2015).  Privacy for free: Posterior sampling and stochastic gradient monte carlo.  In Proceedings of the International Conference on Machine Learning, pages 2493–2502.
  * Zhu et al., (2021) Zhu, K., Van Hentenryck, P., and Fioretto, F. (2021).  Bias and variance of post-processing in differential privacy.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 11177–11184.

##  Appendix A Proofs of Composition and Post-Processing

###### Proof of Proposition 6.

Let
$c_{1}=\frac{e^{\epsilon_{1}\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon_{1}}-1}\delta_{1}$,
$c_{2}=\frac{e^{\epsilon_{2}\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon_{2}}-1}\delta_{2}$
and $x\wedge y:=\min(x,y)$. For any fixed output $(y_{1},y_{2})$ we have

$\begin{split}&\Pr[M_{1,2}({\vec{x}})=(y_{1},y_{2})]\\\
=&\Pr[M_{1}({\vec{x}})=y_{1}]\Pr[M_{2}({\vec{x}})=y_{2}|M_{1}({\vec{x}})=y_{1}]\\\
\leq&\Pr[M_{1}({\vec{x}})=y_{1}]\\\ &\hskip
20.0pt\cdot(e^{\epsilon_{2}\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M_{2}({\vec{x}}^{\prime})=y_{2}|M_{1}({\vec{x}})=y_{1}]\wedge
1+c_{2})\\\
\leq&\big{(}e^{\epsilon_{2}\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M_{2}({\vec{x}}^{\prime})=y_{2}|M_{1}({\vec{x}})=y_{1}]\wedge
1\big{)}\\\ &\hskip
20.0pt\cdot(e^{\epsilon_{1}\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M_{1}({\vec{x}}^{\prime})=y_{1}]+c_{1})+c_{2}\\\
\leq&e^{(\epsilon_{1}+\epsilon_{2})\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M_{1,2}({\vec{x}}^{\prime})=(y_{1},y_{2})]+(c_{1}+c_{2})\end{split}$

which completes the proof. ∎

###### Proof of Proposition 7.

For $M:\mathcal{X}\to\mathcal{Y}$ which is $(\epsilon,\delta)$-differentially
private on $\equiv$, given a output space $\mathcal{Z}$, and an stochastic
mapping $F$ from $\mathcal{Y}$ to $\mathcal{Z}$, for any event
$E\subset\mathcal{Z}$, we have

$\begin{split}&\Pr[F(M({\vec{x}}))\in E]\\\ =&\mathbb{E}_{f\sim
F}[\Pr[f(M({\vec{x}}))\in E]]\\\ =&\mathbb{E}_{f\sim F}[\Pr[M({\vec{x}})\in
f^{-1}(E)]]\\\ \leq&\mathbb{E}_{f\sim
F}\Big{[}e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[M({\vec{x}}^{\prime})\in
f^{-1}(E)]+\frac{e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon}-1}\delta\Big{]}\\\
=&e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}\Pr[F(M({\vec{x}}^{\prime}))\in
E]+\frac{e^{\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|}-1}{e^{\epsilon}-1}\delta\end{split}$

where the second last inequality follows form Definition 2. ∎

##  Appendix B Proofs for Laplace mechanisms

###### Proof of LABEL:thm:lap.

By definition, the mechanism is $\mathcal{A}$-invariant. Now we show the
mechanism is $\epsilon$-$\mathcal{A}$-induced integer subspace differentially
private on $\mathcal{A}$. For all
${\vec{x}}\equiv_{\mathcal{A}}{\vec{x}}^{\prime}$ and a feasible outcome
${\vec{y}}$ so that ${\vec{y}}\equiv_{\mathcal{A}}{\vec{x}}$, by Proposition 8
${\vec{y}}\in{\vec{x}}+\Lambda_{\mathcal{A}}={\vec{x}}^{\prime}+\Lambda_{\mathcal{A}}$.
The privacy loss is

$\displaystyle\ln\frac{\Pr[M({\vec{x}})={\vec{y}}]}{\Pr[M({\vec{x}}^{\prime})={\vec{y}}]}=$
$\displaystyle\ln\frac{q_{\epsilon}(\|{\vec{y}}-{\vec{x}}\|)}{q_{\epsilon}(\|{\vec{y}}-{\vec{x}}^{\prime}\|)}$
$\displaystyle\cdot\ln\frac{\sum_{{\vec{w}}\in{\vec{x}}+\Lambda_{\mathcal{A}}}q_{\epsilon}(\|{\vec{w}}-{\vec{x}}\|)}{\sum_{{\vec{w}}\in{\vec{x}}^{\prime}+\Lambda_{\mathcal{A}}}q_{\epsilon}(\|{\vec{w}}-{\vec{x}}^{\prime}\|)}$

For the first term,

$\ln\frac{q_{\epsilon}(\|{\vec{y}}-{\vec{x}}\|)}{q_{\epsilon}(\|{\vec{w}}-{\vec{x}}^{\prime}\|)}=-\epsilon\left(\|{\vec{y}}-{\vec{x}}\|-\|{\vec{y}}-{\vec{x}}^{\prime}\|\right)\leq\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|,$

by the triangle inequality of $\|\cdot\|$. For the second term, ${\vec{w}}$ is
feasible if and only if ${\vec{w}}+{\vec{x}}-{\vec{x}}^{\prime}$ is still
feasible because ${\vec{x}}\equiv_{\mathcal{A}}{\vec{x}}^{\prime}$. Therefore

$\displaystyle\sum_{{\vec{w}}\in{\vec{x}}+\Lambda_{\mathcal{A}}}q_{\epsilon}(\|{\vec{w}}-{\vec{x}}\|)$
$\displaystyle=\sum_{{\vec{w}}+{\vec{x}}-{\vec{x}}^{\prime}\in{\vec{x}}+\Lambda_{\mathcal{A}}}q_{\epsilon}(\|{\vec{w}}+{\vec{x}}-{\vec{x}}^{\prime}-{\vec{x}}\|)$
$\displaystyle=\sum_{{\vec{w}}\in{\vec{x}}^{\prime}+\Lambda_{\mathcal{A}}}q_{\epsilon}(\|{\vec{w}}-{\vec{x}}^{\prime}\|)$

and the second term is zero. Combining these two we have

$\ln\frac{\Pr[M({\vec{x}})={\vec{y}}]}{\Pr[M({\vec{x}}^{\prime})={\vec{y}}]}\leq\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|$

which completes the proof. ∎

Now we show the unbiasedness and accuracy guarantee in LABEL:thm:lap_acc. To
control the error’s tail bound, we first show that lattice spaces are similar
to Euclidean spaces which are “regular” and have bounded dimensions. We use
the idea of doubling dimension that bounds the grow rate of the number of
lattice points in spheres as the radius increases in LABEL:lem:doubling
Formally, given a positive definite matrix ${\bm{Q}}\in\mathbb{R}^{l\times l}$
and $s>0$, we define the ellipsoid as a set

$E_{\bm{Q}}(s)=\\{{\vec{v}}\in\mathbb{R}^{l}:{\vec{v}}^{\top}{\bm{Q}}{\vec{v}}\leq
s^{2}\\}.$

We use $\operatorname{Vol}(E_{\bm{Q}}(s))$ as the Lebesgue measure of
$E_{\bm{Q}}(s)$ and
$\operatorname{Vol}_{\mathbb{Z}}(E_{\bm{Q}}(s)):=|\mathbb{Z}^{l}\cap
E_{\bm{Q}}(s)|$ as the number of points in $E_{\bm{Q}}(s)$ with integer
coordinates. There is a long line of work dating back to (Landau,, 1915;
Tsuji,, 1953; Bentkus and Gotze,, 1999) showing that the number of lattice
points in $E_{\bm{Q}}(s)$ is approximated by the volume of $E_{\bm{Q}}(s)$,
when $s$ is large enough.

######  Theorem 17.

(Landau,, 1915) For all $\ell\geq 2$ and ${\bm{Q}}$, when $s>0$ is large
enough, we have

$\frac{|\operatorname{Vol}(E_{\bm{Q}}(s))-\operatorname{Vol}_{\mathbb{Z}}(E_{\bm{Q}}(s))|}{\operatorname{Vol}(E_{\bm{Q}}(s))}=\mathcal{O}\left(s^{-2+\frac{2}{1+\ell}}\right)$

Additionally, let the volume of a unit $l$-dimensional sphere be
$V_{l}=\operatorname{Vol}(\mathbb{I})=\frac{\pi^{l/2}}{\Gamma(l/2+1)}$ where
$\Gamma$ is the gamma function. The volume of a general ellipsoid
$E_{\bm{Q}}(s)$ of dimension $l$ is
$V_{l}\left(\frac{s^{2}}{\det({\bm{Q}})}\right)^{l/2}$ where $\det({\bm{Q}})$
is the determinant of ${\bm{Q}}$. Therefore, we have for all positive definite
${\bm{Q}}$ and $\rho$ there exists $s_{\rho}$ so that for all $s\geq s_{\rho}$

$\frac{1}{\rho}s^{l}\leq\frac{\operatorname{Vol}_{\mathbb{Z}}(E_{\bm{Q}}(s))}{\frac{V_{l}}{\sqrt{\det({\bm{Q}})}}}\leq\rho
s^{l}.$ (12)

Informally, as the radius of the ellipsoid doubles, the number of integer
points in the ellipsoid increases by a factor $2^{l}$, where $l$ is the
dimension. Now we are ready to prove LABEL:lem:doubling and LABEL:thm:lap_acc.

###### Proof of LABEL:lem:doubling.

Observe that

$\displaystyle|\\{{\vec{z}}\in\Lambda({\bm{B}}):\|{\vec{z}}\|\leq r\\}|=$
$\displaystyle|\\{{\vec{v}}\in\mathbb{Z}^{m}:\|{\bm{B}}{\vec{v}}\|\leq r\\}|$
$\displaystyle=$
$\displaystyle|\\{{\vec{v}}\in\mathbb{Z}^{m}:{\vec{v}}^{\top}({\bm{B}}^{\top}{\bm{B}}){\vec{v}}\leq
r^{2}\\}|$ $\displaystyle=$
$\displaystyle\operatorname{Vol}_{\mathbb{Z}}\left(E_{{\bm{B}}^{\top}{\bm{B}}}(r)\right)$

Because ${\bm{B}}^{\top}{\bm{B}}$ is full rank, by Eq. 12, there exists
$r^{*}>0$ so that for all $r\geq r^{*}$,

$|B_{\mathcal{A}}(r)|\leq\frac{2V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}r^{m},$

which completes the proof. ∎

###### Proof of LABEL:thm:lap_acc.

Let ${\vec{Z}}=M_{Lap,\mathcal{A},\epsilon}({\vec{x}})-{\vec{x}}$ in
LABEL:def:lap be a random vector.

First it is easy to show that the noise of the generalized Laplace mechanism
is unbiased, because the density
$\exp\left(-\epsilon\|{\vec{z}}\|\right)\mathbf{1}[{\vec{z}}\in\Lambda_{\mathcal{A}}]$
is an even function, for all ${\vec{z}}\in\mathbb{Z}^{d}$

$\exp\left(-\epsilon\|-{\vec{z}}\|\right)\mathbf{1}[-{\vec{z}}\in\Lambda_{\mathcal{A}}]=\exp\left(-\epsilon\|{\vec{z}}\|\right)\mathbf{1}[{\vec{z}}\in\Lambda_{\mathcal{A}}].$

Given $\Lambda_{\mathcal{A}}\subseteq\mathbb{R}^{d}$, let
$D_{\mathcal{A}}(r)=\\{{\vec{z}}\in\Lambda_{\mathcal{A}}:\|{\vec{z}}\|\leq
r\\}$. The tail probability of ${\vec{z}}$ is

$\Pr[\|{\vec{z}}\|\geq
t]=\frac{\sum_{{\vec{z}}\in\Lambda_{\mathcal{A}}\setminus
D_{\mathcal{A}}(t)}\exp(-\epsilon\|{\vec{z}}\|)}{\sum_{{\vec{z}}\in\Lambda_{\mathcal{A}}}\exp(-\epsilon\|{\vec{z}}\|)}.$

We can upper bound it by lower bounding the denominator and upper bounding the
numerator.

For the denominator,

$\sum_{{\vec{z}}\in\Lambda_{\mathcal{A}}}\exp(-\epsilon\|{\vec{z}}\|)\geq\exp(-\epsilon\|\mathbf{0}\|)=1.$

For the numerator, by LABEL:lem:doubling, because ${\bm{C}}_{\mathcal{A}}$ has
full rank, there exists a constant
$L=\frac{2V_{d-k}}{\sqrt{\det({\bm{C}}_{\mathcal{A}}^{\top}{\bm{C}}_{\mathcal{A}})}}$,

$|D_{\mathcal{A}}(r)|\leq Lr^{d-k},\text{ for all large enough }r>0.$ (13)

Hence, if $t$ is large enough, we have

$\displaystyle\sum_{{\vec{z}}\in\Lambda_{\mathcal{A}}\setminus
D_{\mathcal{A}}(t)}\exp(-\epsilon\|{\vec{z}}\|)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{{\vec{z}}\in
D_{\mathcal{A}}(2^{l+1}t)\setminus
D_{\mathcal{A}}(2^{l}t)}\exp(-\epsilon\|{\vec{z}}\|)$ $\displaystyle\leq$
$\displaystyle\sum_{l=0}^{\infty}\sum_{{\vec{z}}\in
D_{\mathcal{A}}(2^{l+1}t)\setminus D_{\mathcal{A}}(2^{l}t)}\exp(-\epsilon
2^{l}t)$ $\displaystyle\leq$
$\displaystyle\sum_{l=0}^{\infty}|D_{\mathcal{A}}(2^{l+1}t)|\exp(-\epsilon
2^{l}t)$ $\displaystyle\leq$
$\displaystyle\sum_{l=0}^{\infty}L(2^{l+1}t)^{d-k}\exp(-\epsilon 2^{l}t)$ (by
Eq. 13 and large enough $t$) $\displaystyle=$ $\displaystyle
L(2t)^{d-k}\sum_{l=0}^{\infty}2^{l(d-k)}\exp(-\epsilon 2^{l}t)$
$\displaystyle\leq$ $\displaystyle 2L(2t)^{d-k}\exp(-\epsilon t)$ (by ratio
test and large enough $t$)

Therefore, we complete the proof by taking $K=2\cdot 2^{d-k}L=\frac{4\cdot
2^{d-k}V_{d-k}}{\sqrt{\det({\bm{C}}_{\mathcal{A}}^{\top}{\bm{C}}_{\mathcal{A}})}}$.
∎

##  Appendix C Proofs for Gaussian mechanisms

Similar to the Gaussian distribution on Euclidean space, a Gaussian
distribution on lattices is sub-Gaussian (LABEL:prop:gaussian), with proof
below. This is non-trivial, because the support of Gaussian in Eq. 8 is not
uniform. For instance, suppose $\Lambda$ is not a lattice space, but has
$2^{r^{2}}$ points in each integer radius $r$ shell. Then the probability of
Gaussian on such space is not sub-Gaussian.

######  Proposition 18 (label = prop:gaussian).

Let ${\bm{B}}\in\mathbb{R}^{d\times m}$ be a rank $m\leq d$ matrix. For any
Gaussian ${\vec{Z}}$ on $\Lambda=\Lambda({\bm{B}})$ with
${\vec{c}}={\mathbf{0}}$ and $\sigma>0$, we have
$\mathbb{E}[{\vec{Z}}]={\mathbf{0}}$ and

$\Pr[\|{\vec{Z}}\|_{2}\geq t]\leq
Kt^{m}\exp\left(-\frac{t^{2}}{2\sigma^{2}}\right)$

when $t>0$ is large enough and $K$ is defined in LABEL:thm:lap_acc.

###### Proof of LABEL:prop:gaussian.

The proof is mostly identical to LABEL:thm:lap_acc. The unbiasedness is
trivial. Given $\Lambda\subseteq\mathbb{R}^{d}$, let
$D(r)=\\{{\vec{z}}\in\Lambda:\|{\vec{z}}\|\leq r\\}$. The tail probability of
${\vec{z}}$ is

$\Pr[\|{\vec{z}}\|\geq t]=\frac{\sum_{{\vec{z}}\in\Lambda\setminus
D(t)}\exp(-\frac{1}{2\sigma^{2}}\|{\vec{z}}\|^{2})}{\sum_{{\vec{z}}\in\Lambda}\exp(-\frac{1}{2\sigma^{2}}\|{\vec{z}}\|^{2})}.$

We will upper bound it by lower bounding the denominator and upper bounding
the numerator.

For the denominator,

$\sum_{{\vec{z}}\in\Lambda}\exp(-\frac{1}{2\sigma^{2}}\|{\vec{z}}\|^{2})\geq\exp(-\epsilon\|\mathbf{0}\|^{2})=1.$

For the numerator, by LABEL:lem:doubling, for all large enough $r>0$,

$|D(r)|\leq\frac{2V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}r^{m}.$ (14)

Hence, if $t$ is large enough, we have

$\displaystyle\sum_{{\vec{z}}\in\Lambda\setminus
D(t)}\exp(-\frac{1}{2\sigma^{2}}\|{\vec{z}}\|^{2})$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{{\vec{z}}\in D(2^{l+1}t)\setminus
D(2^{l}t)}\exp(-\frac{1}{2\sigma^{2}}\|{\vec{z}}\|^{2})$ $\displaystyle\leq$
$\displaystyle\sum_{l=0}^{\infty}\sum_{{\vec{z}}\in D(2^{l+1}t)\setminus
D(2^{l}t)}\exp(-\frac{1}{2\sigma^{2}}2^{2l}t^{2})$ $\displaystyle\leq$
$\displaystyle\sum_{l=0}^{\infty}|D(2^{l+1}t)|\exp\left(-\frac{2^{2l}t^{2}}{2\sigma^{2}}\right)$
$\displaystyle\leq$
$\displaystyle\frac{2V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}\sum_{l=0}^{\infty}(2^{l+1}t)^{m}\exp\left(-\frac{2^{2l}t^{2}}{2\sigma^{2}}\right)$
(by Eq. 14 and large enough $t$) $\displaystyle=$
$\displaystyle\frac{2(2t)^{m}V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}\sum_{l=0}^{\infty}2^{lm}\exp\left(-\frac{2^{2l}t^{2}}{2\sigma^{2}}\right)$
$\displaystyle\leq$
$\displaystyle\frac{2^{m+1}V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}t^{m}\exp\left(-\frac{t^{2}}{2\sigma^{2}}\right)$
(by ratio test and large enough $t$)

Therefore, we complete the proof by taking
$K=\frac{2^{m+1}V_{m}}{\sqrt{\det({\bm{B}}^{\top}{\bm{B}})}}$. ∎

###### Proof of LABEL:thm:gaussian.

By definition, the mechanism is $\mathcal{A}$-invariant. Now we show the
mechanism is $(\epsilon,\delta)$ $\mathcal{A}$-induced integer subspace
differentially private. For all distinct
${\vec{x}}\equiv_{\mathcal{A}}{\vec{x}}^{\prime}$ and a feasible outcome
${\vec{y}}$, because the normalizing constants are identical, the privacy loss
is

$\displaystyle\ln\frac{\Pr[M({\vec{x}})={\vec{y}}]}{\Pr[M({\vec{x}}^{\prime})={\vec{y}}]}=$
$\displaystyle\ln\frac{\exp\left(-\frac{1}{2\sigma_{\epsilon,\delta}^{2}}\|{\vec{y}}-{\vec{x}}\|^{2}\right)}{\exp\left(-\frac{1}{2\sigma_{\epsilon,\delta}^{2}}\|{\vec{y}}-{\vec{x}}^{\prime}\|^{2}\right)}$
$\displaystyle=\frac{1}{2\sigma_{\epsilon,\delta}^{2}}\left(\|{\vec{y}}-{\vec{x}}^{\prime}\|^{2}-\|{\vec{y}}-{\vec{x}}\|^{2}\right),$

If we set ${\vec{z}}:={\vec{y}}-{\vec{x}}$, by cosine theorem the privacy loss
is

$\ln\frac{\Pr[M({\vec{x}})={\vec{y}}]}{\Pr[M({\vec{x}}^{\prime})={\vec{y}}]}\leq\frac{1}{2\sigma_{\epsilon,\delta}^{2}}\left(\|{\vec{x}}-{\vec{x}}^{\prime}\|^{2}+2\|{\vec{x}}-{\vec{x}}^{\prime}\|\|{\vec{z}}\|\right).$

Therefore, the privacy loss is bounded by
$\epsilon\|{\vec{x}}-{\vec{x}}^{\prime}\|$ when ${\vec{z}}$ satisfies

$\|{\vec{z}}\|\leq\sigma_{\epsilon,\delta}^{2}\epsilon-\frac{1}{2}\|{\vec{x}}-{\vec{x}}^{\prime}\|.$

If the datasets ${\vec{x}}$ and ${\vec{x}}^{\prime}$ are too far apart so that
$\|{\vec{x}}-{\vec{x}}^{\prime}\|\geq
2\sigma_{\epsilon,\delta}^{2}\epsilon=\frac{4c_{\mathcal{A}}\ln
1/\delta}{\epsilon}$, we cannot have any privacy guarantee as discussed in
Definition 2. On the other hand, when $\|{\vec{x}}-{\vec{x}}^{\prime}\|\leq
2\sigma_{\epsilon,\delta}^{2}\epsilon$, let
$r=\|{\vec{x}}-{\vec{x}}^{\prime}\|$ and
$t=\sigma_{\epsilon,\delta}^{2}\epsilon-\frac{1}{2}r$. If we can show

$Kt^{d-k}\exp\left(-\frac{t^{2}}{2\sigma_{\epsilon,\delta}^{2}}\right)\leq
e^{\epsilon(r-1)}\delta\leq\frac{e^{\epsilon r}-1}{e^{\epsilon}-1}\delta,$
(15)

we completes the proof by LABEL:prop:gaussian.

We first take log on both side of Eq. 15 and the above inequality is
equivalent to $\epsilon(r-1)-\ln
1/\delta+\frac{t^{2}}{2\sigma_{\epsilon,\delta}^{2}}-(d-k)\ln t-\ln K\geq 0$.

$\displaystyle\frac{t^{2}}{2\sigma_{\epsilon,\delta}^{2}}+\epsilon(r-1)-\ln
1/\delta-(d-k)\ln t-\ln K$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma_{\epsilon,\delta}^{2}\epsilon^{2}-\frac{1}{2}\epsilon
r+\frac{r^{2}}{8\sigma_{\epsilon,\delta}^{2}}+\epsilon(r-1)-\ln 1/\delta$ (16)
$\displaystyle-(d-k)\ln t-\ln K$ (replace
$t=\sigma_{\epsilon,\delta}^{2}\epsilon-\frac{1}{2}r$) $\displaystyle\geq$
$\displaystyle\frac{1}{2}\sigma_{\epsilon,\delta}^{2}\epsilon^{2}-\epsilon-\ln
1/\delta-(d-k)\ln t-\ln K$ $\displaystyle=$ $\displaystyle c_{\mathcal{A}}\ln
1/\delta-\epsilon-\ln 1/\delta-(d-k)\ln t-\ln K$ (17)

Now we show Eq. 17 is greater than zero for some
$c_{\mathcal{A}}=\mathcal{O}(\max\\{(d-k)\ln(d-k),\ln K\\})$ that prove Eq.
15. We show this in three parts. First because
$\delta<\epsilon<e^{-\epsilon}$, $\ln 1/\delta-\epsilon>0$. Second, for some
$c_{\mathcal{A}}=\mathcal{O}(\max\\{(d-k)\ln(d-k),\ln K\\})$, we have
$\frac{1}{3}c_{\mathcal{A}}\ln 1/\delta\geq\frac{c_{\mathcal{A}}}{3}>\ln K$.
Finally,

$\displaystyle\left(\frac{2}{3}c_{\mathcal{A}}-1\right)\ln 1/\delta-(d-k)\ln
k$ $\displaystyle\geq$
$\displaystyle\left(\frac{2}{3}c_{\mathcal{A}}-1\right)\ln
1/\delta-(d-k)\ln\sigma_{\epsilon,\delta}^{2}\epsilon$ (because
$t\leq\sigma_{\epsilon,\delta}^{2}\epsilon$) $\displaystyle=$
$\displaystyle\left(\frac{2}{3}c_{\mathcal{A}}-1\right)\ln
1/\delta-(d-k)\left(\ln\ln 1/\delta+\ln 1/\epsilon+\ln c_{\mathcal{A}}\right)$
$\displaystyle\geq$
$\displaystyle\left(\frac{2}{3}c_{\mathcal{A}}-1-2(d-k)\right)\ln
1/\delta-(d-k)\ln c_{\mathcal{A}}$ ($\delta<\epsilon<1/e$) $\displaystyle\geq$
$\displaystyle 0$

The last inequality holds for some $c_{\mathcal{A}}=O((d-k)\ln(d-k))$ that
completes the proof. ∎

##  Appendix D Details on Implementation

###  D.1 The Gibbs-within-Metropolis sampler

Algorithm 1 presents a Gibbs-within-Metropolis sampler, which produces a
sequences of dependent draws $\left({\vec{z}}^{\left(l\right)}\right)_{0\leq
l\leq\text{nsim}}$ from the target distribution $q_{\epsilon}$, using an
additive jumping distribution whose element-wise construction is as described
in Eq. 11.

Algorithm 1 A Gibbs-within-Metropolis sampler

1: input: initial distribution $\pi_{0}$, pre-jump distribution
$\\{\eta_{j}\\}_{k+1\leq j\leq n}$, target distribution $q_{\epsilon}$,
${\bm{V}}\in\mathbb{Z}^{d\times d}$ from the Smith normal form of constraint
${\bm{A}}\in\mathbb{Z}^{k\times d}$ of rank $k$;

2: Initialize ${\vec{v}}^{\left(0\right)}\sim\pi_{0}$;

3: for $l=1,\ldots,\text{nsim}$ do

4: set ${\vec{v}}^{*}={\vec{v}}^{\left(l-1\right)}$;

5: for $j=k+1,\ldots,n$ do

6: sample ${e}_{j}^{*}\sim\eta_{j}$ and update
${\vec{v}}^{*}_{[j]}={\vec{v}}^{\left(l-1\right)}_{[j]}+e^{*}_{j}$;

7: end for

8: sample $r\sim U[0,1]$;

9: if
$r\leq{q_{\epsilon}\left({\bm{V}}{\vec{v}}^{*}\right)}/{q_{\epsilon}\left({\bm{V}}{\vec{v}}^{\left(l-1\right)}\right)}$
then set ${\vec{v}}^{\left(l\right)}={\vec{v}}^{*}$, otherwise set
${\vec{v}}^{\left(l\right)}={\vec{v}}^{\left(l-1\right)}$;

10: set ${\vec{z}}^{\left(l\right)}={\bm{V}}{\vec{v}}^{\left(l\right)}$;

11: end for

12: return: a chain $\left({\vec{z}}^{\left(l\right)}\right)_{0\leq
l\leq\text{nsim}}$ following the single transition kernel
${\vec{z}}^{\left(l\right)}\sim
K\left({\vec{z}}^{\left(l-1\right)},\cdot\right)$ with target distribution
$q_{\epsilon}$.

The choice of $\eta_{j}$ is a tuning decision for the algorithm. If using the
double geometric proposal as described in Section 4, the scale parameter $a$
should be set to encourage fast mixing of the chain. For the examples in
Section 5, we report the specific choices as they occur, as well as in
Appendix E. The initial distribution $\pi_{0}$ may be chosen simply as $g$ for
convenience.

As noted in Section 4, steps 5 through 7 of Algorithm 1 updates the proposed
jump in a Gibbs sweep (i.e. one dimension at a time), utilizing the fact that
the pre-jump object ${\vec{e}}$ is element-wise independent under
$g_{\bm{A}}$. Doing so facilitates the $L$-lag coupling for assessing
empirical convergence of the chain. In practice, these updates may also be
performed simultaneously.

###  D.2 $L$-lag maximal coupling algorithm and meeting time sampler.

We adapt the method proposed by (Biswas et al.,, 2019) to perform empirical
convergence assessment of the proposed Metropolis algorithm using $L$-lag
coupling. As Section 4 explains, we construct a joint transition kernel
$\tilde{K}$ of two Markov chains, each having the same target distribution
$q_{\epsilon}$ induced by the marginal transition kernel $K$ as defined in
Algorithm 1. The joint kernel $\tilde{K}$ is a maximal coupling kernel and is
given in Algorithm 2. It is designed in such a way that the $L$-th lag of the
two chains will couple in finite time with probability one. This $L$-lag
meeting time is a random variable, and is sampled via Algorithm 3, which
employs Algorithm 1 and 2 as subroutines.

Algorithm 2 A maximal $L$-lag Gibbs-within-Metropolis coupling

1: input: lag $L\geq 1$, iteration $l>L$, states
$\left({\vec{z}}^{\left(l-1\right)},{\vec{z}}_{\circ}^{\left(l-L-1\right)}\right)$

2: Initialize
${\vec{v}}^{\left(l-1\right)}={\bm{V}}^{-1}{\vec{z}}^{\left(l-1\right)}$,
${\vec{w}}^{\left(l-L-1\right)}={\bm{V}}^{-1}{\vec{z}}_{\circ}^{\left(l-L-1\right)}$

3: for $j=k+1,\ldots,n$ do

4: set ${\vec{v}}^{*}={\vec{v}}^{\left(l-1\right)}$ and
${\vec{w}}^{*}={\vec{w}}^{\left(l-L-1\right)}$;

5: sample ${e}_{j}^{*}\sim\eta_{j}$ and update
${\vec{v}}^{*}_{[j]}={\vec{v}}^{\left(l-1\right)}_{[j]}+e^{*}_{j}$;

6: sample $s\sim U[0,1]$;

7: if
$s\eta_{j}\left(e^{*}_{j}\right)\leq\eta_{j}\left({\vec{v}}^{*}_{[j]}-{\vec{w}}^{(l-L-1)}_{[j]}\right)$
then update ${\vec{w}}^{*}_{[j]}={\vec{v}}^{*}_{[j]}$; else sample
$\tilde{e}_{j}\sim\eta_{j}$, $\tilde{s}\sim U[0,1]$, and set
$\tilde{w}_{j}={\vec{w}}^{(l-L-1)}_{[j]}+\tilde{e}_{j}$ until
$\tilde{s}\eta_{j}\left(\tilde{e}_{j}\right)>\eta_{j}\left(\tilde{w}_{j}-{\vec{v}}^{(l-1)}_{[j]}\right)$;
update ${\vec{w}}^{*}_{[j]}=\tilde{w}_{j}$.

8: end for

9: sample $r\sim U[0,1]$;

10: if
$r\leq{q_{\epsilon}\left({\bm{V}}{\vec{v}}^{*}\right)}/{q_{\epsilon}\left({\bm{V}}{\vec{v}}^{\left(l-1\right)}\right)}$
then set ${\vec{v}}^{\left(l\right)}={\vec{v}}^{*}$, else set
${\vec{v}}^{\left(l\right)}={\vec{v}}^{\left(l-1\right)}$;

11: if
$r\leq{q_{\epsilon}\left({\bm{V}}{\vec{w}}^{*}\right)}/{q_{\epsilon}\left({\bm{V}}{\vec{w}}^{\left(l-L-1\right)}\right)}$
then set ${\vec{w}}^{\left(l\right)}={\vec{w}}^{*}$, else set
${\vec{w}}^{\left(l-L\right)}={\vec{w}}^{\left(l-L-1\right)}$;

12: set ${\vec{z}}^{\left(l\right)}={\bm{V}}{\vec{v}}^{\left(l\right)}$,
${\vec{z}}_{\circ}^{\left(l-L\right)}={\bm{V}}{\vec{w}}^{\left(l-L\right)^{\prime}}$.

13: return: a pair of draws
$\left({\vec{z}}^{\left(l\right)},{\vec{z}}_{\circ}^{\left(l-L\right)}\right)$
following the joint transition kernel
$\tilde{K}\left(\left({\vec{z}}^{\left(l-1\right)},{\vec{z}}_{\circ}^{\left(l-L-1\right)}\right),\cdot\right)$,
each with marginal target distribution $q_{\epsilon}$.

The output of Algorithm 3 is used to construct an estimate of the upper bound
on the total variation distance between the target distribution $q_{\epsilon}$
and the marginal distribution of the chain at time $l$, denoted
$q_{\epsilon}^{(l)}$ (Biswas et al.,, 2019, Theorem 2.5). The upper bounds are
obtained as empirical averages over independent runs of coupled Markov chains.

Algorithm 3 $L$-lag meeting time sampler (adapted from (Biswas et al.,, 2019))

1: input: lag $L\geq 1$, initial distribution $\pi_{0}$, single kernel $K$ and
joint kernel $\tilde{K}$

2: output: meeting time $\tau^{(L)}$, chains
$\left({\vec{z}}^{(l)}\right)_{0\leq l\leq\tau^{(L)}}$,
$\left({\vec{z}}_{\circ}^{(l)}\right)_{0\leq l\leq\tau^{(L)}-L}$

3: Initialize ${\vec{z}}^{(0)}\sim\pi_{0}$,
${\vec{z}}^{\left(l+1\right)}\mid{\vec{z}}^{\left(l\right)}\sim
K\left({\vec{z}}^{\left(l\right)},\cdot\right)$ for $1\leq l\leq L$, and
${\vec{z}}_{\circ}^{(0)}\sim\pi_{0}$

4: for $l>L$ do

5: sample
$\left({\vec{z}}^{\left(l\right)},{\vec{z}}_{\circ}^{\left(l-L\right)}\right)\sim\tilde{K}\left(\left({\vec{z}}^{\left(l-1\right)},{\vec{z}}_{\circ}^{\left(l-L-1\right)}\right),\cdot\right)$

6: if ${\vec{z}}^{\left(l\right)}={\vec{z}}_{\circ}^{\left(l-L\right)}$ then
return $\tau^{(L)}:=l$ and chains $\left({\vec{z}}^{(l)}\right)_{0\leq
l\leq\tau^{(L)}}$ and $\left({\vec{z}}_{\circ}^{(l)}\right)_{0\leq
l\leq\tau^{(L)}-L}$

7: end for

For Algorithm 2, one needs to choose a lag $L>0$. As (Biswas et al.,, 2019)
discusses, smaller $L$ encourages faster coupling of the two chains, however
the estimated total variation distance as a upper bound is a poorer one. On
the other hand, larger $L$ produces a a tighter bound at the expense of a
heavier computational burden. (Biswas et al.,, 2019) recommends starting with
$L=1$, and increasing it to a point where the estimated TV upper bound
$d_{TV}(q_{\epsilon}^{(0)},q_{\epsilon})$ is close to 1 (i.e. non-vacuous).
(Biswas et al.,, 2019) also advises against increasing $L$ arbitrarily,
because the TV upper bound itself is not an optimal one thus the benefit of
sharpness is limited. In our experiments in Section 5, we employ $L$ values
such that the distribution of the $L$-lag coupling time, $\tau^{(L)}-L$,
appears to be empirically stable in $L$. The choices of $L$ are reported as
they occur, as well as in Appendix E.

##  Appendix E Additional experimental results and details

###  E.1 Additional details

The code and data for reproducing plots is included in the supplementary
material. All experiments were run on a Macbook with M1 processor and 16GB
memory. The packages and software used are all available freely.

###  E.2 Intersecting counting constraints.

Consider a synthetic example in which a set of three counting constraints
$\mathcal{A}=(A_{1},A_{2},A_{3})$ are defined over 14 records, where the
intersection of all subsets of constrains are nonempty. The constraints are
schematically depicted by Fig. 3. The corresponding incident matrix is
${\bm{A}}\in\mathbb{Z}^{3\times 14}$, with each row corresponding to one of
$A_{1},A_{2}$ and $A_{3}$, and each element being $1$ at indices corresponding
to the record within that constraint and $0$ otherwise. We apply Algorithm 1
to instantiate the generalized Laplace mechanism with both $\ell_{1}$-norm and
$\ell_{2}$-norm targets, to privatize a data product that conforms to the
constraint ${\bm{A}}$ with $\epsilon=0.25$. The pre-jump proposal
distributions $\eta_{j}$ are double geometric distributions, with parameter
$a=\exp{(-1)}$ for the $\ell_{1}$-norm target and $a=\exp{(-1.5)}$ for the
$\ell_{2}$-norm target. Fig. 4 shows the estimated upper bound total variation
distance using $L$-lag coupling. The chains appear to empirically converge at
around $10^{5}$ iterations for the $\ell_{1}$-norm target distribution and
$10^{6}$ iterations for the $\ell_{2}$-norm target distribution.

Figure 3: Schematic depiction of the intersecting counting constraints
$\mathcal{A}=(A_{1},A_{2},A_{3})$.

Figure 4: Intersecting constraints: estimated upper bound on the chain’s total
variation distance from target distribution using $L$-lag coupling (Biswas et
al.,, 2019). Left: $\ell_{1}$-norm target; right: $\ell_{2}$-norm target by
the generalized Laplace mechanism with $\epsilon=0.25$.

Fig. 5 shows the boxplot of 1000 realizations of the proposed entry-wise
errors respectively under the $\ell_{1}$-norm and $\ell_{2}$-norm target
distributions. The proposed additive noise are integer-valued and are entry-
wise unbiased.

Figure 5: Intersecting counting constraints: boxplots of 1000 realizations of
proposed entrywise errors; $\ell_{1}$-norm target distribution (top) and
$\ell_{2}$-norm target distribution (bottom), $\epsilon=0.25$.

###  E.3 Delinquent Children by County and Household Head Education Level.

In the _report on statistical disclosure limitation (SDL) methodology_ , the
Federal Committee on Statistical Methodology published a fictitious dataset
concerning delinquent children in the form of a $4\times 4$ contingency table,
tabulated across four counties by education level of household head (Federal
Committee on Statistical Methodology,, 2005, Table 4). Reproduced in Table 2,
the dataset was used to illustrate various traditional SDL techniques, and was
employed again by (Slavković and Lee,, 2010) in a comparative study on the
effect of swapping versus synthetic table generation on disclosure risk and on
the validity of downstream statistical analysis.

Table 2: $4\times 4$ contingency table: delinquent children by county and household head education level (reproduced from Table 4 of (Federal Committee on Statistical Methodology,, 2005)) County | Low | Medium | High | Very High | Total
---|---|---|---|---|---
Alpha | 15 | 1 | 3 | 1 | 20
Beta | 20 | 10 | 10 | 15 | 55
Gamma | 3 | 10 | 10 | 2 | 35
Delta | 12 | 14 | 7 | 2 | 35
Total | 50 | 35 | 30 | 20 | 135

Consider the privatization of the data while preserving the margins of the
contingency table. We apply Algorithm 1 to instantiate the generalized Laplace
mechanism with both $\ell_{1}$\- and $\ell_{2}$-norm targets and with
$\epsilon=0.25$. The pre-jump proposal distributions $\eta_{j}$ are double
geometric distributions, with parameter $a=\exp(-1)$ for the $\ell_{1}$-norm
target distribution, and $a=\exp(-2)$ for the $\ell_{2}$-norm target
distribution.

Fig. 1 in Section 5 show the evolution of the TV upper bounds on the chain’s
marginal distributions to the $\ell_{1}$-norm target distribution, each
estimated with 200 independent coupled chains, and Table 1 shows one
realization of the noises obtained after the chain has achieved empirical
convergence under the $\ell_{1}$ norm target. Table 3 and Fig. 6 presented
here are analogs of Table 1 and Fig. 1, respectively, but for the
$\ell_{2}$-norm target distribution. Under the $\ell_{2}$-norm target, the
chains appear to converge after about $10^{5}$ iterations, and are stable at
various choices of $L$ as shown in the figures.

County | Low | Medium | High | Very High | Total
---|---|---|---|---|---
Alpha | -11 | 1 | 10 | 0 | 0
Beta | -3 | -6 | 1 | 8 | 0
Gamma | -6 | -5 | 13 | -2 | 0
Delta | 20 | 10 | -24 | -6 | 0
Total | 0 | 0 | 0 | 0 | 0

Table 3:  Proposed generalized Laplace additive noise ($\epsilon=0.25$,
$\ell_{2}$-norm target distribution) Algorithm 1, which preserves row and
column margins.

Figure 6: Estimated upper bound: TV distance from target ($\ell_{2}$ norm)

In addition, Fig. 7 shows the boxplot of 1000 realizations (thinned at 0.01%)
of the proposed errors under both the $\ell_{1}$-norm and $\ell_{1}$-norm
target distributions. The proposed additive noise are integer-valued and are
cellwise unbiased.

Figure 7: $4\times 4$ contingency table: boxplot of 1000 realizations of
proposed cell-wise errors; $\ell_{1}$-norm target distribution (top) and
$\ell_{2}$-norm target distribution (bottom), $\epsilon=0.25$.

###  E.4  2010 U.S. Census county-level population data.

We consider the publication of county-level population counts subject to the
invariant of state population size. The U.S. Census Bureau published
iterations of privacy-protected demonstration files produced by preliminary
versions of its 2020 Disclosure Avoidance System (DAS), which treats the 2010
Census Summary Files (CSF) as the confidential values. These data have been
curated by IPUMS NHGIS and are publicly available (Van Riper et al.,, 2020).
Employed here are the November 2020 vintage demonstration data, protected by
pure differential privacy with $\epsilon=0.192=4\text{ (total)}\times
0.16\text{ (county level)}\times 0.3\text{ (population query)}$.

We demonstrate Algorithm 2 using the generalized Laplace mechanism under
$\ell_{1}$ norm, with $\epsilon$ set to accord to the Census Bureau’s
specification. The states that are used for demonstration in this work are
identified by (Gao et al.,, 2022) as those for which the Census DAS produced
statistically significantly biased errors (at the $\alpha=0.01$ level) for the
county populations. There are 11 states in total. The result for the state of
Illinois is presented in Fig. 2 in Section 1. Fig. 8 below presents results
for 10 additional states.

To assess the empirical convergence of the sampler, for each state we
initialize four independent chains from over-dispersed starting values. To
obtain these starting values, independent chains are run for $10^{6}$
iterations with an over-dispersed target distribution, specifically
$q_{\epsilon}$ with $\epsilon=0.1$. The pre-jump proposal parameter is set to
$a=\exp(-2.5)$ and is the same for all the states. With those starting values,
each independent chains again expends $10^{6}$ steps of burn-in with the
desired target distribution, $q_{\epsilon}$ with $\epsilon=0.192$. The
potential scale reduction factors calculated across independent chains are
observed to reach $<1.01$ for each county of each state.

Figure 8:  Privacy noise from DAS demonstration data (red dots) (11/2020
vintage; Van Riper et al.,, 2020) vs. generalized Laplace mechanism (
LABEL:def:lap) via Algorithm 1 (blue squares: one instance; boxplot: 1000
instances) for county populations within a state. The $x$-axes are arranged in
increasing true county populations. State population total is invariant. The
states shown here (plus Illinois in Fig. 2) are those identified by (Gao et
al.,, 2022) for which the DAS produced statistically significantly biased
errors (at the $\alpha=0.01$ level), as can be seen by the downwardly biased
trends among the red. In contrast, the proposed noises are integer-valued and
unbiased.

Figure 9: continued.

#### Comparison with projected discrete Laplace mechanism

For publishing just the county level data under invariants, the most basic
approach would be to add noise which ensures $\epsilon$-DP and then projecting
onto space satisfying the constraints. As the latter portion of this procedure
constitutes post-processing, the output is still private. One drawback with
this approach is the procedure is now less transparent with regards to the
randomness involved and the utility is now unclear due to optimization
procedure required to project onto constraint satisfying subspace. Integer
subspace differential privacy address both these concerns and though
empirically noise magnitudes are similar we have a better understanding of the
process compared to adding noise and then post-processing.

To compare, we consider databases that are at distance 1, i.e.
$\|x-x^{\prime}\|_{1}=1$ and compare generalized Laplace mechanism to
$\epsilon-$DP mechanism for sum constraint (the sum of all entries is held
invariant). The $\epsilon-$DP mechanism is obtained by adding Double geometric
noise and then solving linear program to find closest vector in $\ell_{2}$
norm that satisfies the sum constraint. We then consider another vector closes
in $\ell_{1}$ norm but with additional constraint that variables are now
integers. We do not consider non-negativity constraints. Fig. 10 depicts in
green - box plots per county for 1000 independent simulations of the above
mentioned procedure whereas the setup for our mechanism is same as before. The
box-plots majorly overlap indicating the magnitude of noise is comparable for
both approaches.

#### Comparison with TopDown Algorithm (Abowd et al.,, 2022)

The TopDown algorithm, in addition to invariants and integer values, also
performs post-processing to ensure that all released counts are non-negative.
This requirement of non-negativity introduces systematic bias into the
privatized tabulations (Zhu et al.,, 2021; Gao et al.,, 2022). Hotz and Salvo,
(2020) note that earlier versions of the DAS introduces inaccuracy in small-
area data that compromised existing statistical analyses. An example may be
seen from Fig. 2: when the state population total is fixed, positive errors
tend to associate with smaller county population counts, and negative errors
with larger counts. Through six rounds of iterations spanning 2019-2021, the
bureau increased the privacy loss budget for the PL. 94-171 files from
$\epsilon=6$ (U.S. Census Bureau,, 2019) to
$(\epsilon,\delta)=(19.71,10^{-10})$ (U.S. Census Bureau,, 2021).

When non-negativity is enforced, bias is theoretically unavoidable. However,
bias comes in a matter of degrees. The DAS demonstration dataset that we
showcase here is tabulated at the county level, where the underlying microdata
are dense and the confidential counts are far from zero. As a consequence, the
non-negativity post-processing of counts privatized with moderately sized
unbiased noise should not induce visible systematic bias as shown in Fig. 2.
To be precise, if we infuse the county-level counts with errors from the
proposed generalized Laplace mechanism (shown in Fig. 2 and 8) and then apply
post-processing to truncate negative counts, none of the 1000 instantiations
triggers the truncation in any county of any state. This is not surprising,
since the smallest counties of most states have population sizes in the
thousands (with the least populous county of the U.S. being Loving, Texas with
a population of 82), whereas the overwhelming majority of our proposed noise
instantiations are within [-30, 30] using the Bureau’s privacy loss budget
allocation for the county population table. Thus, even after nonnegativity
post-processing is applied to our proposed solution, the resulting privatized
counts are still effectively unbiased, because the post-processing applied to
county-level tables is effectively vacuous.

A key design feature of the TopDown algorithm is that privatized counts at
different geographic levels are produced separately, rather than only at the
block level and aggregated all the way up. As Section 6.5 of Abowd et al.,
(2022) explains, this takes advantage of the hierarchical structure of the
Census geographies and facilitates fine-grained control over the Bureau’s own
fitness-for-use accuracy targets. While maintaining unbiasedness and
nonnegativity is difficult at the block level, it is achievable at coarser
levels as we demonstrate here. When the underlying geography is dense in
counts, our proposal contributes a solution that delivers privatized
tabulations that are effectively unbiased even if nonnegativity is furthermore
imposed.

Figure 10:  Fig. 8 superimposed with additional post-processing on privacy
noise generated from the generalized Laplace mechanism (green squares: one
instance; boxplot: 1000 instances) to enforce non-negativity of the privatized
county-level population counts. In this case, post-processing did not alter
the distribution of the original noise from the generalized Laplace mechanism
( LABEL:def:lap via Algorithm 1; blue squares: one instance; boxplot: 1000
instances) in any significant way, and did not induce a downward bias as seen
from the privacy noise from DAS demonstration data (red dots) (11/2020
vintage; Van Riper et al.,, 2020).

Figure 11: continued.

Figure 12: continued.