# A normal approximation for joint frequency estimatation under Local
Differential Privacy

Thomas Carette  Sony R&D

###### Abstract

In the recent years, Local Differential Privacy (LDP) has been one of the
corner stone of privacy preserving data analysis. However, many challenges
still opposes its widespread application. One of these problems is the
scalability of LDP to high dimensional data, in particular for estimating
joint-distributions. In this paper, we develop an approximate estimator for
frequency joint-distribution estimation under so-called pure LDP protocols
[1].

##  1 Introduction

With the advances of data-driven technologies in the last decades, privacy has
become a hot topic at all levels of society. Correspondingly, the litterature
on methods that allow statistical analysis and inference while protecting the
anonymity of the data subjects has exploded.

One such technique is Differential Privacy [2] and its extension, Local
Differential Privacy (LDP) [3]. The key contribution of (Local) Differential
Privacy is to quantify, with parameters $\epsilon$ and $\delta$, the degree of
anonymization of disclosed data. The parameters $\epsilon$ and $\delta$ are
usually chosen depending on the problem at hand. One of its features is to
ensure _plausible deniability_ which is at the root of the invention of the
_Randomized Response_ [4] method in 1965. Randomized Response (RR), applied on
binary data, and its generalization to categorical data (GRR), are
prototypical LDP protocols.

Since the introduction of the LDP formalism by Dwork and co-workers, a vast
literature has derived from marrying its methodology with (G)RR, as attested
by two recent comprehensive reviews [5, 6]. While the privacy preserving
collection of a single discrete information is relatively well understood,
generic solution to the problem of accessing multi-dimensional statistics can
still be improved.

One of the corner stone of recent research stems from the observation that
splitting the privacy budget $\epsilon$ for concurrent question leads to an
exponentially growing variance, but splitting the population so that each user
answers a single question only leads to a linear growth of the variance [1].
Of course, the downside of the latter approach is that the correlation between
the different answers is lost. Various algorithms have been proposed to
collect joint-distribution using LDP, based on expectation maximization [7, 8]
or other approaches [9]. Both methods suffer from scalability issues in
general. On the question of $k$-way marginal release, which is core to a
series of statistical methods, the CALM method [10] was introduced to try to
strike a balance between the two strategies. However this balance remains
dataset- and task-dependent, meaning higher level algorithms need to be
designed.

Another question, and opportunity, lies in how to efficiently encode a complex
dataset. Complexity can stem from large or unknown dictionary sizes as well as
from high dimensionality and low entropy. Any categorical dataset can be
turned one-dimensional at the cost of an exponentially large dictionary.

When using binary encoding of the dictionary, the randomized response can be
applied on each bit independently. This approach has been explored
extensively. Various method achieves good or optimal variance-utility
tradeoffs for frequency estimation [11, 1], but often rely on data transfers
increasing exponentially with the number of features. Once again, some level
of optimization can be done by randomizing lossy encoding over users, thereby
reducing communication and focussing the use of the privacy budget for each
user [12, 13, 1], but by limiting the amount of information each user
contribute, large populations are required to obtain useful signal/noise ratio
on very small marginals.

In order to push back the limit on the number of features reported by single
users, Ren  _et al._ [8] proposed a multi-demensional version of RAPPOR [12,
7] by essentially concatenating the Bloom filters of each feature rather than
considering the full row as an element of a dictionary. The advantage of this
approach is that the size of bit array to communicate grows only linearly with
the number of feature instead of exponentially. Also, this decomposition of
the information allows easy access to lower-dimensional marginals without the
need of suming over all the noisy cells of the full contingency table. The
issue with this method is that reconstructing the final joint-distribution can
be rather expensive. Indeed, the only method proposed so far was to obtain it
numerically using Expectation Maximization (EM). Some acceleration or
stabilisation of the convergence can be achieved by using clever initial
guesses [8, 14], but this approach is still unsatisfactory in most cases.

Another issue with using numerically estimated joint distribution is that it
constitutes a dead-end for analytical research of new algorithms. As an
illustration, the research performed on analytical description of the
distribution of LDP frequency estimations was critical to investiate theoretic
optimal solutions [15, 16, 1]. However, to our knowledge, no such description
exist that efficiently describes the distribution of _marginals_ in a multi-
dimensional dataset. Many LDP statistical learning algorithms, some of them
iterative, rely on the knowledge of conditional probabilities in the dataset
[5, 6] but because those conditional probabilities are not known analytically,
researchers must rely on costly simulations and intuition.

In this paper, we tackle the problem of analytical approximation of the
posterior joint distribution of a discrete dataset after perturbation. We
focus on pure LDP protocols proposed by [1]. By solving this problem we hope
to facilitate the development of new algorithms leveraging the analytical
properties of the deduced equations.

##  2 Background and notations

In this paper, we will use boldface font to denote random variables or
functions, and normal font for values.

$\varepsilon,\delta$ | traditional DP parameters | definition 2.1
---|---|---
A | random function defining a LDP protocol | definition 2.1
$D,D^{\prime}$ | representations of datasets | definition 2.1
$\mathcal{X},\mathcal{Y}$ | input and output spaces of a LDP protocol | definition 2.1
$p,q$ | traditional randomized response paraeters | definition 2.3
$S$ | Support of a random function A | definition 2.3
$\mathbb{S}$ | Support operator corresponding to $S$ | equation 8
t | ground truth set of all joint frequencies of values in a dataset $D$ |
${\bf\hat{X}}$ |  _evidence matrix_ of perturbed dataset | section 2.4
c | set of all marginal frequencies of ${\bf\hat{X}}$ (PMDH) | equations 16-18
t | estimator of ground truth computed with c (TMDH) | section 2.4
$\mathcal{I},\mathcal{J},\mathcal{S}$ | Sets of columns in ${\bf\hat{X}}$ | section 3.1
$\alpha,\beta$ | size of column sets | section 3.1
$\mathbb{K}^{(\alpha\beta)}(\mathcal{J},\mathcal{S})$ | Subset incidence matrix for $\mathcal{S}\subseteq\mathcal{J}$ | equation 26
Table 1: Notation

###  2.1 Differential Privacy

Differential privacy is defined in the context of a trusted curator that
process the users data and releases anonymized reports. In differential
privacy, the level of anonymity is defined using a statistical framework [2].

######  Definition 2.1.

An random function A, mapping elements of its domain $\mathcal{X}$ to its
image $\mathcal{Y}$, satisfies $(\varepsilon,\delta)$-differential privacy
with $\varepsilon>0$, if and only if for any two datasets
$D,D^{\prime}\in\mathbb{N}^{|\mathcal{X}|}$ satisfying $||D-D^{\prime}||=1$,
_i.e._ that differ only in one element, and for all $y\in\mathcal{Y}$

$P({\bf A}(D)=y)\leq e^{\varepsilon}P({\bf A}(D^{\prime})=y)+\delta$

This definition establishes a framework to _measure_ privacy leakage in the
parameter $\varepsilon$ with a certain slack $\delta$. In order to lift the
requirement to rely on a trusted curator, Kasiviswanathan  _et al._ [3]
introduced a stricter definition, _Local Differential Privacy_. In this case,
rather than ensuring that the random algorithm has an output distribution that
is close for neighboring databases, it is imposed that any two elements have a
finite probability ratio to have the same image

######  Definition 2.2.

An random function A, mapping elements of its domain $\mathcal{X}$ to its
image $\mathcal{Y}$, satisfies $(\varepsilon,\delta)$-local differential
privacy with $\varepsilon>0$, if and only if for any two inputs
$x,x^{\prime}\in\mathcal{X}$ we have for all $y\in\mathcal{Y}$

$P({\bf A}(x)=y)\leq e^{\varepsilon}P({\bf A}(x^{\prime})=y)+\delta$

We will focus on the cases where $\delta=0$, _i.e._ strictly no accidental
privacy leakage is allowed.

###  2.2 Randomized response

The randomized response protocol was introduced by Warner [4] to guarantee
anonymity in surveys. It can be defined formally on binary answers as

$\displaystyle P({\mathbf{R}}_{p,q}(1)=1)$ $\displaystyle=$ $\displaystyle q$
(1) $\displaystyle P({\mathbf{R}}_{p,q}(0)=1)$ $\displaystyle=$ $\displaystyle
p$ (2)

where ${\mathbf{R}}_{p,q}$ represents the randomized response function.

Note that when $0<p<q<1$, the space of ${\mathbf{R}}_{p,q}$ form a monoid with
the composition operation $\circ$, i.e. it forms an algebra which has all the
properties of a group except for the inversion.

$\displaystyle{\mathbf{R}}_{p(q^{\prime}-p^{\prime})+p^{\prime},q(q^{\prime}-p^{\prime})+p^{\prime}}$
$\displaystyle=$
$\displaystyle{\mathbf{R}}_{p,q}\circ{\mathbf{R}}_{p^{\prime},q^{\prime}}\qquad\textrm{(complete)}$
(3) $\displaystyle{\mathbf{R}}_{p,q}$ $\displaystyle=$
$\displaystyle{\mathbf{R}}_{p,q}\circ{\mathbf{R}}_{0,1}\qquad\textrm{(identity)}$
(4)
$\displaystyle{\mathbf{R}}_{p,q}\circ({\mathbf{R}}_{p^{\prime},q^{\prime}}\circ{\mathbf{R}}_{p^{\prime\prime},q^{\prime\prime}})$
$\displaystyle=$
$\displaystyle({\mathbf{R}}_{p,q}\circ{\mathbf{R}}_{p^{\prime},q^{\prime}})\circ{\mathbf{R}}_{p^{\prime\prime},q^{\prime\prime}}\qquad\textrm{(associative)}$
(5)

The identity is a direct consequence of the relation of completeness and the
composition operation is always associative. Note also that we need both
$p=1-q$ and $p^{\prime}=1-q^{\prime}$ to have commutativity between the
composition of two randomized response functions.

###  2.3 Pure LDP protocols

Wang  _et al._ [1] introduced a family of randomized-response-like protocols,
dubbed _pure LDP protocols_.

######  Definition 2.3.

An encoding or hashing function
$A_{e}:\mathcal{X}\rightarrow\mathcal{Y}_{e}\subseteq\mathcal{Y}$, a random
function $\mathbf{A}_{p}:\mathcal{Y}_{e}\rightarrow\mathcal{Y}$, and a support
function $S:\mathcal{Y}\rightarrow\bar{\mathcal{X}}$ with $\bar{\mathcal{X}}$
the set of sets over $\mathcal{X}$, form a _Pure LDP protocol_ if we have two
probability values $p<q$ such that for all $x,x^{\prime}\in\mathcal{X}$,
$x\neq x^{\prime}$

$\displaystyle P(\mathbf{A}(x)\in\\{y|x\in S(y)\\})$ $\displaystyle=$
$\displaystyle q$ (6) $\displaystyle P(\mathbf{A}(x^{\prime})\in\\{y|x\in
S(y)\\})$ $\displaystyle=$ $\displaystyle p$ (7)

where $\mathbf{A}=\mathbf{A}_{p}\circ A_{e}$

$\mathcal{X}$ is the dataset dictionary of values, while $\mathcal{Y}_{e}$ is
the dictionary of encoded values (which can be smaller than the original
dictionary). The space $\mathcal{Y}$ can be outside of the $\mathcal{Y}_{e}$
dictionary space so that the support function cannot be considered as a
decoder.

The above definition states that components of an encoding of an input value
are conserved with a probability $q$, and that there is a probability $p$ than
any other input value would share the same component111Unfortunately, p and q
are used interchangeably in the litterature. We adopt the convention of the
RAPPOR paper [12] while Wang  _et al._ [1] use $p$ for the probability of
conserving the original value, and $q$ for the probability of polution.. It
can be seen as a generalization of the Randomized Response protocol, which is
itself a pure LDP protocol. As noted by the original authors, it is a direct
consequence of the definition 2.3 that the use of Bloom filters as encoding
function, as in RAPPOR, does not lead to a pure protocol as soon as any
hashing collision $h_{i}(x)=h^{\prime}_{i}(x^{\prime})$ occurs, $h$ and
$h^{\prime}$ being any hashing function of the filter.

The support operator maps perturbed vectors in the Unary Encoding space

$\mathbb{S}_{x}(y)=1\textrm{ if }x\in S(y)\textrm{ and }0\textrm{ otherwise}$
(8)

Wang  _et al._ [1] prove that

$\mathbf{t}_{x}=\frac{\sum_{i}\mathbb{S}_{x}(\mathbf{A}(x_{i}))-np}{q-p}$ (9)

is an unbiased estimator for the frequency of any element $i$, where
$x,\\{x_{i}\\}\in\mathcal{X}$, and $i$ spans all the elements of the dataset.
Note that the random variable $\mathbb{S}_{x}(\mathbf{A}(x_{i}))$ follow
Bernouilli distributions

$\mathbb{S}_{x}(\mathbf{A}(x_{i}))\sim\left\\{\begin{array}[]{ll}\textrm{Bernouilli}(q)&\textrm{if
}x=x_{i}\\\ \textrm{Bernouilli}(p)&\textrm{if }x\neq x_{i}\\\
\end{array}\right.$ (10)

Two variables $\mathbb{S}_{x}(\mathbf{A}(x_{i}))$ and
$\mathbb{S}_{x}(\mathbf{A}(x_{j}))$ obviously follow independent
distributions. However, for a single input they are in general not
independent, _i.e._
$\mathbb{S}_{x}(\mathbf{y})\centernot{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}\mathbb{S}_{x^{\prime}}(\mathbf{y})$.

###  2.4 Multi-dimensionality, marginals and pure LDP protocols

Wang  _et al._ [1] focussed on instances of pure LDP protocols applied to
element frequency estimation, but it is not restricted to this case. The
advantage of considering explicitely the multi-dimensionality of the dataset
is to be able to avoid the exponential growth of the dictionary with the
number of categories, which poses problems for communication costs,
information loss during encoding, or marginal retrieval.

A first step to fight ths curse of dimensionality is to follow in the
footsteps of Ren  _et al._ [8] and split the dataset row domain into
independent sub-domains for each feature. The encoding function $A_{E}$
encodes each feature independently and concatenate the results. Similarly, the
perturbation function $\mathbf{A}_{P}$ can be a combination of randomization
steps each acting in its own feature sub-domain. Finally the support function
$S$ as well can be defined as mapping each feature encoding space into the
corresponding dictionary. In this setting, any arbitrary combination of
feature-specific pure LDP protocol will form a global pure LDP-protocol.

If we have $f$ features, then we have for the dataset value space
$\mathcal{X}$ as well as for the encoded value space $\mathcal{Y}$

$\displaystyle\mathcal{X}$ $\displaystyle=$
$\displaystyle\mathcal{X}^{(1)}\oplus\cdots\oplus\mathcal{X}^{(f)}$ (11)
$\displaystyle\mathcal{Y}$ $\displaystyle=$
$\displaystyle\mathcal{Y}^{(1)}\oplus\cdots\oplus\mathcal{Y}^{(f)}$ (12)

and

$\mathbb{S}_{x}(\mathbf{A}_{p}(A_{e}(x_{i})))=\mathbb{S}^{(1)}_{x}(\mathbf{A}^{(1)}_{p}(A^{(1)}_{e}(x^{(1)}_{i})))\oplus\cdots\oplus\mathbb{S}^{(f)}_{x}(\mathbf{A}^{(f)}_{p}(A^{(f)}_{e}(x^{(f)}_{i})))$
(14)

where the superscripts indicate the feature sub-domain and the $\oplus$ is the
direct sum of the resulting elements. In the following, we will assume that a
basis set is built on such a decomposition of $\mathcal{X}$ and that any index
$i$ indicates a specific value $x^{(k)}$ belonging to a specific feature space
$\mathcal{X}^{(k)}$, thereby allowing us to unambiguously omit the feature
superscript in the notation. We define the _evidence matrix_ random variable

${\bf\hat{X}}_{ni}=\mathbb{S}_{i}(\mathbf{A}_{p}(A_{e}(x_{n})))$ (15)

For a specific $i$, it forms a sequence of independent Bernouilli experiments
so that the marginal sum

${\mathbf{c}}_{i}=\sum_{n=0}^{N}{\bf\hat{X}}_{ni}$ (16)

follows the distribution of sum of Binomial random variables, one of
probability $p$ and the other of probability $q$. We can also construct higher
order perturbed marginal counts

$\displaystyle{\mathbf{c}}_{ij}$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{N}{\bf\hat{X}}_{ni}{\bf\hat{X}}_{nj}$ (17)
$\displaystyle{\mathbf{c}}_{ijk}$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{N}{\bf\hat{X}}_{ni}{\bf\hat{X}}_{nj}{\bf\hat{X}}_{nk}$
(18) $\displaystyle\cdots$

Those counts also follow the distribution of sum of Binomial random variables.
Considering that every possible row in ${\bf\hat{X}}$ occurs with a certain
probability, it is straightforward that a change of basis set in the histogram
space reduces the joint distribution of all c to a sum of multinomial
distributions. This change of basis set is explicited in appendix A and
mulinomial properties are used in the lemma of appendix B.3. However, it
cannot be used to derive an analytical posterior for the joint frequency
distribution.

In this paper, we perform an analytical derivation of multivariate normal
parameters for the joint distribution
$\\{{\mathbf{c}}_{i},{\mathbf{c}}_{ij},\dots\\}$, later referred to as the
Perturbed Multi-Dimensional Histogram (PMDH). This allows to use Bayesian
Inference to deduce the Truth Multi-Dimensional Histogram (TMDH)

$P(A|B)\propto P(A)P(B|A)$ (19)

with
$A=(\\{{\mathbf{t}}_{i},{\mathbf{t}}_{ij},\dots\\}=\\{t_{i},t_{ij},\dots\\})$
and
$B=()\\{{\mathbf{c}}_{i},{\mathbf{c}}_{ij},\dots\\}=\\{c_{i},c_{ij},\dots\\})$.
We demonstrate the use of this model by computing analytically the
distribution of the LDP correlation estimator.

Some generalization of our results could be performed using the properties of
the normal distribution under linear basis transformations. Furthermore, the
developments presented in the appendix of this paper can serve as a basis for
dealing with more complex cases. A relatively straightforward but tedious work
would be to adapt them to LH or GRR of each independent feature space.

The decomposition of c in independent distributions can easily be done if
$\mathbb{S}_{x}\circ\mathbf{A}_{p}\circ A_{e}$ is block diagonal and that the
index $i,j,k,\dots$ belong to different blocks, _i.e._ if
${\bf\hat{X}}_{ni}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}{\bf\hat{X}}_{nj}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}{\bf\hat{X}}_{nk}$.
In general, $\mathbf{(}A(x_{n}))$ is composed of multiple independent random
variables. For instance one for GRR, and $|\mathcal{X}|$ for UE. If
${\bf\hat{X}}_{ni}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}{\bf\hat{X}}_{nj}$,
then the set of random variables they depend on is disjoint, and they compound
to the probabilities $p$ or $q$ according to 2.3. If
${\bf\hat{X}}_{ni}\centernot{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}{\bf\hat{X}}_{nj}$
it means they depend on common random variables. For instance, using Bloom
filters over set values data (to avoid collisions), two values $x$ and
$x^{\prime}$ can share a bit in their support function, the other bits of the
filter serving to disambiguate the values. In LH and GRR, all dictionary
values are supported by the same underlying multinomial random variable.

We list a few examples of correlated co-occurrence probabilities in Table 2.
In GRR, all co-occurrence counts ${\mathbf{c}}_{ij},{\mathbf{c}}_{ijk},\dots$
are simply 0 if any pair of indices belong to the same block, while in LH they
are merely modulated depending on the size of the hashing dictionary.

Besides standard hashing and encoding techniques, $A_{e}$ can also implement
various linear transformations like for instance in the Hadamard Randomized
Response [17] which also lead to dependencies between evidence matrix
elements. A careful investigation of all those aspects is beyond the scope of
the present work. In the following, we assume
${\bf\hat{X}}_{ni}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}{\bf\hat{X}}_{nj}$.

Table 2:  Co-occurrence probabilities. $p_{i}=E({\bf\hat{X}}_{ni})$. The
probability is the success probability of $\prod_{i=1}^{k}{\bf\hat{X}}_{ni}$.
LH: $d$ is the dimension of the hash space
Protocol | probability
---|---
UE | $\prod_{i}p_{i}$
GRR | 0
LH | $p_{1}/d^{k}$

###  2.5 Previous work

Different aspects of joint-distribution have been considered, but few concern
the analytical development of marginal estimators. Fanti  _et al._ [7]
proposed the use an EM-algorithm to reconstruct the joint-distribution of a
dataset encoded row by row. Ren  _et al._ [8] also considered using EM,
optimizing the solution by opting for a column-wise encoding scheme and
speeding up convergence using an initial guess obtained with Lasso. Wang  _et
al._ [14] proposed an analytical, iterative formula for the co-occurrences
frequency estimation, but they didn’t derive the co-variance matrix between
estimators so that EM was still necessary to obtain the joint-distribution.
Cormode  _et al._ [9] showed the superiority of Hadamard (or Fourier)
transforms to collects marginals of binary datasets. The CALM framework of
Zhang  _et al._ [10] further improves numerical derivation of marginals by
offering a systematic way to neglect specific correlations, and using maximum
entropy to extrapolate between lower dimensional marginals. These lower-
dimensional marginals are still constructed using established methods of full
row value frequency estimation.

In this work, we leverage the concept of pure LDP protocols [1] (as discussed
in section 2.4), inspired by the LoPub approach to multi-dimensional datasets
[8], and develop those following the empirical estimation methodology [15,
16].

##  3 Normal parameters for Multi-Dimensional Histograms

According to the Central Limit Theorem (CLT), the distribution
$\mathbf{C}=\\{{\mathbf{c}}_{i},{\mathbf{c}}_{ij},{\mathbf{c}}_{ijk},\dots\\}$
converges to a normal distribution for large $N$. We have

$\displaystyle\lim_{N\rightarrow\infty}\mathbf{C}/N$ $\displaystyle\sim$
$\displaystyle\mathcal{P}(c)\
{\mathcal{N}}\left(\hat{\mu},\hat{\Sigma};t\right)$ (20)
$\displaystyle\lim_{N\rightarrow\infty}\mathbf{T}/N$ $\displaystyle\sim$
$\displaystyle\mathcal{P}(t)\ {\mathcal{N}}\left(\mu,\Sigma;c\right)$ (21)

where $\mathcal{P}(t)$ and $\mathcal{P}(c)$ are priors. Without any specific
prior knowledge, a so-called “non-informative prior” is chosen over the TMDH
domain. This domain is constrained in general by the fact that the individual
frequencies must be comprised between 0 and 1, that counts are discrete, and
that co-occurrences cannot be larger than the corresponding frequencies. Other
constrains are specific to the encoding function which, for instance, could
guarantee that categories in a single dimension are mutually exclusive.

In general, the dimension of this representation makes it computationally
intractable, but in practice it is possible to truncate it to low-order co-
occurrences as in most cases high-order co-occurrences are extremely unlikely
to be statistically significant. It is however not our purpose to propose
scalable algorithms using the PMDH and TMDH distributions, but rather to lay
the mathematical ground work to be built upon.

In section 3.1, we derive analytical formulas for $\hat{\mu}$ and
$\hat{\Sigma}$ for arbitrary co-occurrence order. In section 3.2, we derive
the associated inferred parameters $\mu$ and $\Sigma$.

For defining arbitrary co-occurrences, we will be using curved letters
$\mathcal{I}$, $\mathcal{J}$, $\mathcal{S}$,…as sets of indices for which bits
corresponding to a single row are positive, either before or after the
randomization process depending on the context. The number of elements in such
a set, _i.e._ the number of positive bits in the row, will be denoted
$|\mathcal{I}|$.

###  3.1 PMDH normal parameters

####  3.1.1 Expected values

######  Theorem 3.1.

The expected value of a co-occurrence of order $|\mathcal{I}|$ of indices
$\mathcal{I}=[i_{1},\dots,i_{k}]$ can be developed as

$E({\mathbf{c}}_{\mathcal{I}})=\sum_{\mathcal{S}\subset\mathcal{I}}(q-p)^{|\mathcal{S}|}p^{|\mathcal{I}|-|\mathcal{S}|}t_{\mathcal{S}}$
(22)

where $t_{\mathcal{\emptyset}}=N$, the number of elements in the dataset.

The proof can be found in Appendix B.2. Because the multi-variate expected
value $E(\mathbf{C})$ is built out of the expected value of the
marginals$\\{E({\mathbf{c}}_{i}),E({\mathbf{c}}_{ij}),\dots\\}$ , we can write
the above theorem in matrix form. Our basis set for the multi-dimensional
histogram is split in the order of co-occurrence considered: $\alpha=1$
corresponds to basic counts ${\mathbf{c}}_{i}$, $\alpha=2$ to pair-wise co-
occurrences ${\mathbf{c}}_{ij}$, $\alpha=3$, 3rd order co-occurrences
${\mathbf{c}}_{ijk}$, etc. Except stated otherwise, we will use greek letters
$\alpha$, $\beta$, $\gamma$ to denote indices of _blocks_ in the histrogram
vector representation, _i.e._ range of rows and columns corresponding to co-
occurrence counts of a given order. From theorem 3.1 we have

$\hat{\mu}=M\cdot T+PN$ (23)

defining

$\displaystyle M$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{D}\sum_{\beta=1}^{\alpha}p^{\alpha-\beta}(q-p)^{\beta}\mathbb{K}^{(\alpha\beta)}$
(24) $\displaystyle P$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{D}p^{\alpha}\mathbb{V^{(\alpha)}}$ (25)

where $\mathbb{V^{(\alpha)}}$ the vector of ones in the $\alpha$ block and 0
elsewhere. The matrix $M$ is triangular. The matrices
$\mathbb{K}^{(\alpha\beta)}$ are binary matrices defined to map all sets of
indices $\mathcal{S}$ of size $\beta$ to the set of indices $\mathcal{J}$ of
size $\alpha$. It is defined element-wise as

$\mathbb{K}^{(\alpha\beta)}(\mathcal{J},\mathcal{S})=1\ \mathrm{iif}\
\mathcal{S}\subseteq\mathcal{J}$ (26)

An example for three columns is given in appendix B.5. Note that
$\mathbb{K}^{\alpha\alpha}$ is the identity matrix in the $\alpha$ block and 0
elsewhere. The result above relies on the assumption that we know the dataset
size, _i.e._ $t_{\emptyset}=N$ ($\alpha=0$), $V(N)=0$.

Anticipating on subsequent development, it is also possible to include it in
the vector representation, hence dropping the $PN$ term of equation 23 and to
treat the knowledge of $N$ as a prior in the inference process. Alternatively,
we can account for overdetermined encoding strategies. For instance, with
unary encoding, we can account for a contrain $\sum_{i\in\mathcal{I}}t_{i}=N$
for a pertinent subset of columns $\mathcal{I}$. For a single category, the
modified $M$ takes the form $\tilde{M}=M+b\mathbb{J}$ with $\mathbb{J}_{ij}=1$
$\forall$ $i,j$. The inverse of that matrix is easily obtained from $M^{-1}$
using the Sherman-Morrison formula

$\tilde{M}^{-1}=M^{-1}-\frac{b\
M^{-1}\mathbb{J}M^{-1}}{1+b{\bar{1}}^{T}M^{-1}{\bar{1}}}$ (27)

where ${\bar{1}}_{i}=1$ $\forall i$. It is trivial to extend this to any
linear constrain on the $t_{i}$’s.

####  3.1.2 Covariance matrix

######  Theorem 3.2.

The variance of a co-occurrence of order $|\mathcal{I}|$, of indices
$\mathcal{I}=[i_{1},\dots,i_{k}]$ can be developed as

$\displaystyle V({\mathbf{c}}_{I})$ $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq\mathcal{I}}t_{\mathcal{S}}\
(q-p)^{|\mathcal{S}|}p^{|\mathcal{I}|-|\mathcal{S}|}(1-(q+p)^{|\mathcal{S}|}p^{|\mathcal{I}|-|\mathcal{S}|})$
(28)

The proof can be found in Appendix B.4.1.

######  Theorem 3.3.

Consider a PMDH built on a dataset of $N$ rows and of columns $\mathcal{I}$
with corresponding underlying ground truth
$T=\\{t_{\mathcal{J}}|\mathcal{J}\subseteq\mathcal{I}\\}$. Consider two
arbitrary column sets $\mathcal{J}$ and $\mathcal{J}^{\prime}$. The covariance
between the co-occurrences counts ${\mathbf{c}}_{\mathcal{J}}$ and
${\mathbf{c}}_{\mathcal{J}}^{\prime}$ is

$V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})=\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
(q-p)^{|\mathcal{S}|}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left(1-(q+p)^{|\mathcal{S}\cap\mathcal{J}_{\cap}|}p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}\right)$
(29)

where $\mathcal{J}_{\cup}=\mathcal{J}\cup\mathcal{J}^{\prime}$ and
$\mathcal{J}_{\cap}=\mathcal{J}\cap\mathcal{J}^{\prime}$.

The proof can be found in Appendix B.4.2.

We can see directly that that Theorem 3.2 is a special case of Theorem 3.3
with $\mathcal{J}_{\cap}=\mathcal{J}_{\cup}=\mathcal{J}=\mathcal{J}^{\prime}$.
Similarly, if $\mathcal{J}_{\cap}=\emptyset$ the correlation vanishes, as
expected. Examples of low-order cases are

$\displaystyle V({\mathbf{c}}_{i},{\mathbf{c}}_{ij})$ $\displaystyle=$
$\displaystyle N\ p^{2}(1-p)+t_{j}\ p(q-p)(1-p)$ (30)
$\displaystyle\qquad\qquad\qquad+t_{i}\ p(q-p)[1-q-p]+t_{ij}\
(q-p)^{2}(1-q-p)$ $\displaystyle V({\mathbf{c}}_{ij},{\mathbf{c}}_{jk})$
$\displaystyle=$ $\displaystyle N\ p^{3}(1-p)+(t_{i}+t_{k})\
p^{2}(q-p)(1-p)+t_{j}\ p^{2}(q-p)(1-q-p)$ (31)
$\displaystyle\qquad\qquad\qquad\ (t_{ij}+t_{jk})\ p(q-p)^{2}(1-q-p)+t_{ik}\
p(q-p)^{2}(1-p)$ $\displaystyle\qquad\qquad\qquad\ t_{ijk}(q-p)^{3}(1-q-p)$

###  3.2 TMDH normal parameters

The likelihood of a certain TMDH value after observing one PMDH
$C=\\{c_{i},c_{ij},\dots\\}$ is obtained in the normal approximation 21 using
equation 23

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle M^{-1}(C-NP)$ (32)
$\displaystyle\Sigma$ $\displaystyle=$ $\displaystyle
M^{-1}\hat{\Sigma}(M^{-1})^{T}$ (33)

Here we must realize that $\hat{\Sigma}$ depends on the ground truth $T$ and
hence cannot be used as-is to compute the likelihood. The simplest way is to
inject the expected value $\mu$ in the computation of $\Sigma$. However, a
more careful approach is sometimes necessary in the presence of prior
information.

In order to compute mean and variance of the TMDH, we first deduce the form of
$M^{-1}$

######  Lemma 3.4.

Given a matrix $M$ of the form 24, its inverse takes the form

$M^{-1}=\sum_{\alpha=1}^{D}(q-p)^{-\alpha}\sum_{\beta=1}^{\alpha}(-p)^{\alpha-\beta}\
\mathbb{K}^{(\alpha\beta)}$ (34)

The proof can be found in appendix C.1. An example for three columns is given
in appendix B.5.

The formulas for the mean and covariance matrix are provided in sections 3.2.1
and 3.2.2 respectively and demonstrated in appendix C.2.

####  3.2.1 Expected values

######  Theorem 3.5.

In estimated true co-occurrence $\mu_{\mathcal{I}}(C)$ on a set of columns
$\mathcal{I}$ depending on an observed PMDH is

$\displaystyle\mu_{\mathcal{I}}(C)$ $\displaystyle=$
$\displaystyle({q-p})^{-|\mathcal{I}|}\sum_{\mathcal{S}\subseteq\mathcal{I}}(-p)^{|\mathcal{I}|-|\mathcal{S}|}c_{\mathcal{S}}$
(35)

For $|\mathcal{I}|=1$ we recover the well known frequency estimator of pure
LDP protocols, and randomized response in particular.

####  3.2.2 Covariance matrix

The law of total variance dictates

$\displaystyle V(\mathbf{T})$ $\displaystyle=$ $\displaystyle
V_{\mathbf{T}}(E_{\mathbf{C}}(\mathbf{T}))+E_{\mathbf{T}}(V_{\mathbf{C}}(\mathbf{T}))$
(36) $\displaystyle=$ $\displaystyle
V_{\mathbf{T}}(\mu)+E_{\mathbf{T}}(\Sigma)$ (37) $\displaystyle=$
$\displaystyle 0+M^{-1}\hat{\Sigma}(E(\mathbf{T}))(M^{-1})^{T}$ (38)
$\displaystyle=$ $\displaystyle M^{-1}\hat{\Sigma}(\mu)(M^{-1})^{T}$ (39)

where $V_{\mathbf{C}}$ and $V_{\mathbf{T}}$ are the variance at fixed
$\mathbf{T}$ and $\mathbf{C}$ respectively, and $E_{\mathbf{C}}$ and
$E_{\mathbf{T}}$ are similarly defined expected values. The last two steps are
possible thanks to the linearity of the expected value operator. The first
term vanishes because the $\mu$ estimator does not depend on ground truth.

######  Theorem 3.6.

Consider a TMDH inferred from a randomized dataset $\hat{X}$ of $N$ rows and
of columns $\mathcal{I}$ characterized by a PMDH
$C=\\{c_{\mathcal{J}}|\mathcal{J}\subseteq\mathcal{I}\\}$. Consider two
arbitrary column sets $\mathcal{J}$ and $\mathcal{J}^{\prime}$. The covariance
between the co-occurrences counts ${\mathbf{t}}_{\mathcal{J}}$ and
${\mathbf{t}}_{\mathcal{J}}^{\prime}$ is

$\displaystyle
V({\mathbf{t}}_{\mathcal{J}},{\mathbf{t}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle\left(\frac{1}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J^{\prime}}|}\\!\\!\sum_{\mathcal{S}\subseteq\mathcal{J}_{\cup}}c_{\mathcal{S}}(-p)^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left((-p)^{|\mathcal{J}_{\cap}\setminus\mathcal{S}|}(1-2p)^{|\mathcal{J}_{\cap}\cap\mathcal{S}|}-(q-p)^{|\mathcal{J_{\cap}}|}\right)$

where $\mathcal{J}_{\cup}=\mathcal{J}\cup\mathcal{J}^{\prime}$ and
$\mathcal{J}_{\cap}=\mathcal{J}\cap\mathcal{J}^{\prime}$.

We see that $\mathcal{J}_{\cap}=\emptyset$ lead trivially to a vanishing
correlation. We can derive the following examples of variances

$\displaystyle V({\mathbf{t}}_{i})$ $\displaystyle=$
$\displaystyle(q-p)^{-2}\left(qpN+(1-p-q)c_{i}\right)$ (41) $\displaystyle
V({\mathbf{t}}_{ij})$ $\displaystyle=$
$\displaystyle(q-p)^{-4}\left(\left((1-2p)^{2}-(q-p)^{2}\right)c_{ij}+p\left(p(1-2p)+(p-q)^{2}\right)(c_{i}+c_{j})-p^{2}q(q-2p)N\right)$
(42)

and co-variances

$\displaystyle V({\mathbf{t}}_{ij},{\mathbf{t}}_{i})$ $\displaystyle=$
$\displaystyle(q-p)^{-3}\left((1-q-p)(c_{ij}-p\ c_{i})+qp(c_{j}-pN)\right)$
(43) $\displaystyle V({\mathbf{t}}_{ij},{\mathbf{t}}_{jk})$ $\displaystyle=$
$\displaystyle(q-p)^{-4}\big{(}(1-p-q)(c_{ijk}-p(c_{ij}+c_{jk})+p^{2}c_{j})+pq(c_{ik}-p(c_{i}+c_{k})+p^{2}N)\big{)}$
(44)

##  4 Validation and discussion

We provide the code for the numerical validation of the normal approximation
in the repository [code]. In this paper, we summarize the key findings of our
numerical validation. These were obtained using an artificial dataset of 4
independent binary variables for which we studies the full joint distribution
under a pure LDP protocol with $q=0.8$ and $p=0.1$. The dataset size was
varied from $\sim 10$ to $\sim 10^{5}$. This arguably corresponds to higher
$\epsilon$ than disirable ($\sim 2$) and smaller dataset size than typical.
However we trust this choice is reasonable as it shows both the transition to
well determined TMDH and the issues that one can encounter with the normal
approximation for very low frequencies. The user is free to download our code
to test other configurations.

###  4.1 Validation of the PMDH

There are two main approximations in the equations 3.1 and 3.2 related to
omitting high order moments in the distribution : we neglect the distortion
close to domain boundaries (_i.e._ when a mean frequency is close to 0 or 1,
or when a co-occurrence is close to the limits imposed by the individual
frequencies) and we assume a linear error (_i.e._ the variance of a marginal
frequency is independent of the error on other marginal frequencies). The
latter can be illustrated by considering the conditional probabilities over
the exact multinomial and comparing them to the normal approximation. It can
be proven that for partial sum following a multinomial
$\\{{\mathbf{c}}_{i},{\mathbf{c}}_{j},{\mathbf{c}}_{ij}\\}$,
$({\mathbf{c}}_{ij}|{\mathbf{c}}_{i}=x,{\mathbf{c}}_{j}=y)\sim\mathcal{H}(N,x,y)$
where $\mathcal{H}$ is a hypergeometric function. This leads to

$E({\mathbf{c}}_{ij}|{\mathbf{c}}_{i}=x,{\mathbf{c}}_{j}=y)=\frac{xy}{N}\quad;\qquad
V({\mathbf{c}}_{ij}|{\mathbf{c}}_{i}=x,{\mathbf{c}}_{j}=y)=\frac{xy(N-x)(N-y)}{N^{2}(N-1)}$
(45)

while te normal normal approximation yields

$E({\mathbf{c}}_{ij}|{\mathbf{c}}_{i}=x,{\mathbf{c}}_{j}=y)\approx
xp_{j}+yp_{i}-Np_{j}p_{i}\quad;\qquad
V({\mathbf{c}}_{ij}|{\mathbf{c}}_{i}=x,{\mathbf{c}}_{j}=y)\approx
Np_{i}p_{j}(1-p_{i})(1-p_{j})$ (46)

where $p_{i},p_{j}$ are the probabilities of success of ${\bf\hat{X}}_{i}$ and
${\bf\hat{X}}_{j}$ respectively. This issue should be kept in mind when trying
to force priors as conditions on the posterior distribution directly.

We try 2 different simulation techniques : full simulation, where the binary
dataset is built and perturbed explicitely, and a more efficient approach
where we sample from the exact sum of multinomial distributions. As expected,
both converge to the same PMDH distribution. We therefore use the latter for
our tests.

Figure 1 shows that, as expected, the calculated normal parameters agree with
the simulations within the numerical accuracy. For metrics emphasizing
divergence to normal, we observe a larger error which diminishes with
increasing dataset size (see Figure 2).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 1: Examples of mean and covariance matrix elements computed
analytically vs. simulated. Bracket on top of figures signify the set of
columns $\mathcal{I}$. (a-c) $E({\mathbf{c}}_{\mathcal{I}})$, (d-f)
$\left(V({\mathbf{c}}_{\mathcal{I}})\right)^{1/2}$, (g-i)
$V({\mathbf{c}}_{\mathcal{I}},{\mathbf{c}}_{\mathcal{J}})$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Examples of metrics showing divergence to normal behavior. Bracket
on top of figures signify the set of columns $\mathcal{I}$ and $\mathcal{J}$.
(a-c) relative divergence of median to mean $1-mean/median$, (d-f)
$Skew({\mathbf{c}}_{\mathcal{I}})$

###  4.2 Validation of the TMDH

To verify the quality of our estimates as a funciton of the dataset size, we
sample PMDH from their exact distribution and compute the corresponding
expected values for the true distribution. As in the case of the PMDH, we
observe that the expected value of our mean TMDH estimator defined in 3.5 is
within numerical uncertainty of the true frequencies. Furthermore, the
variance of this estimator is also closely matched by the TMDH variance
estimator defined in 3.6. As shown in figure 3(a), we observe errors above 1%
only for low frequency ${\mathbf{t}}_{\mathcal{I}}$ in smaller datasets
($<10^{3}$). In those cases, the uncertainty on the variance estimator can
exceed 5% (see figure 3(b)). We do not attribute this error on simulation
noise, but rather on the effect of our approximation close to the boundaries
of the exact distribution. This is issue is clearly visible in the figure 3(f)
where the discretization of the support is visible in the reconstructed
distribution. Metrics of deviation from the normal as a function of dataset
size are displayed in figures 3g and 3h.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: Reconstruction of true frequencies based on PMDH. (a-c) Errors on
normal parameters. (a-b) see text. (c) Difference between covariance of
estimated means and theoretical covariance between ${1,2}$ and various sets
$\mathcal{I}$ (d-f) comparison of the distribution of reconstructed
$E({\mathbf{t}}_{i})$ (blue histogram) and a normal distribution of same mean
and variance (orange curve) for a dataset of 200 datapoints. (g-h) metrics of
deviation to the normal of the distribution of reconstructed
$E({\mathbf{t}}_{i})$ : mean-median deviation (g) and skew (h).

It would be interesting to compare the theoretical distribution of the TMDH to
a numerical solution, e.g. obtained with Expectation Maximization (EM).
However, our purpose here is to demonstrate that we propose a useful estimator
of the joint frequencies, as well as robust error estimates on those
frequencies. Therefore we consider that comparing our method with numerical
approaches is beyond the scope of this paper. Another point we do not study
here is the impact of priors on the estimates. Simple solutions like
truncating and/or rounding the normal distribution will have either a positive
or little impact on our conclusions, while more involved solutions are beyond
the scope of this work.

##  5 Conclusion

In this work, we propose an normal approximation for the joint estimation of
category frequencies of arbitrary dimensionality under pure LDP protocols. We
beleive this can be useful both to further the research on this difficult
problem as well as help speedup current numerical algorithms.

We show numerically that the normal approximation of the TMDH has a wide
domain of applicability, loosing accuracy only close to the histogram
boundaries. Future work should compare our solution to numerical solutions. We
expect using the analytical TMDH could significantly improve the performance
of EM-based algorithms.

We also briefly highlight potential opportunities and caveats of the normal
approximation to account for prior information and constrains on the
distribution support. In section 3.1, we discuss some constrains in the
encoding space which can be handled using equation of the form 27. However,
more subtle priors could be hard to handle. In principle, working with
multivariate gaussian distributions offer all sorts of possibilities to build
conditional probabilities. In practice however, one can easily face numerical
instabilities when working directly on the posterior distribution with sharp
constrains. In section 4.1, we highlight how error can grow when constrains
stray too far from the distribution center. Such issues should be studied more
carefully to enable a wider use of the results reported here.

##  Appendix A Decomposition of the randomized matrix

We denote ${\bf\hat{X}}[\mathcal{I}]$, with $\mathcal{I}$ a column set
$[ij\dots]$, the slice of ${\bf\hat{X}}$ constituted of the columns
$i,j,\dots$ only. We denote ${\bf\hat{X}}^{b}[\mathcal{I}]$ the horizontal
slice of ${\bf\hat{X}}[\mathcal{I}]$ that corresponds to ground truth rows
equal to the binary vector $b=[b_{i},b_{j},\dots]$ defined on $\mathcal{I}$.

######  Definition A.1.

A Synchronous Bernoulli Sequence (SBS) is a series of finite length (denoted
$N$) of $D$ distinguishable synchronous Bernoulli experiments with fixed
(permanent) success probabilities $[p_{1},p_{2},\dots,p_{D}]$ respectively.

The matrix ${\bf\hat{X}}^{b}[\mathcal{I}]$ is analogous to a SBS where each
column is a sequence, each row represent synchronous experiments, and the
perturbed binary value in each cell is a single experiment output. We will
denote $L({\bf\hat{X}})$ the number of rows in a matrix ${\bf\hat{X}}$.

######  Lemma A.1.

For any column subset $\mathcal{I}$, and corresponding row-set defined by the
bit array $b_{\mathcal{J}}=[b_{i}b_{j}\dots]$, denoting
$\mathcal{J}=\\{i|b_{i}=1\\}$, then

$L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{I}])=\sum_{\mathcal{J}\subseteq\mathcal{S}\subseteq\mathcal{I}}(-1)^{|\mathcal{S}|-|\mathcal{J}|}t_{\mathcal{S}}$
(47)

###### Proof.

We will prove the result using recursion.

By definition,
$L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{J}])=t_{\mathcal{J}}$, so that the
result is immediate when $\mathcal{J}=\mathcal{I}$.

When $|\mathcal{I}|-|\mathcal{J}|=1$, _i.e._ there is one and only one zero in
$b_{\mathcal{J}}$, we have

$L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{I}])=t_{\mathcal{J}}-t_{\mathcal{I}}$
(48)

In words, that means that the number of co-occurrences on $\mathcal{J}$ where
the remaining bit of $\mathcal{I}$ is 0 is the number of co-occurrences on
$\mathcal{J}$ minus the number of co-occurrences on $\mathcal{I}$.

If we consider an arbitrary $\mathcal{J}$ with $|\mathcal{J}|<|\mathcal{I}|$,
we know that we have to remove from $t_{\mathcal{J}}$ all the rows
corresponding to bit vectors $b^{\prime}$ in which at least one 0 has been
flipped up compared to $b_{\mathcal{J}}$. We have

$\displaystyle L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{I}])$
$\displaystyle=$ $\displaystyle
t_{\mathcal{J}}-\sum_{b^{\prime}}L({\bf\hat{X}}^{b^{\prime}}[\mathcal{I}])$
(49) $\displaystyle=$ $\displaystyle
t_{\mathcal{J}}-\sum_{\mathcal{J}\subset\mathcal{B}\subseteq\mathcal{I}}L({\bf\hat{X}}^{b^{\prime}}[\mathcal{I}])$
(50)

where $\mathcal{B}=\\{i|b^{\prime}_{i}=1\\}$. If we assume the lemma holds for
all $b^{\prime}$

$\displaystyle L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{I}])$
$\displaystyle=$ $\displaystyle t_{\mathcal{J}}-\
\sum_{\mathcal{J}\subset\mathcal{B}\subseteq\mathcal{I}}\
\sum_{\mathcal{B}\subseteq\mathcal{S}\subseteq\mathcal{I}}(-1)^{|\mathcal{S}|-|\mathcal{B}|}t_{\mathcal{S}}$
(51)

Let us consider the contribution of a single count $t_{\mathcal{S}}$ in this
sum. It is the same as asking, for each $l$, how many distinct $\mathcal{B}$
are subsets of $\mathcal{S}$. It is equivalent to choosing $k-l$ elements
($=|\mathcal{S}\setminus\mathcal{B}|$) from $\mathcal{S}\setminus\mathcal{I}$
to remove from $\mathcal{S}$. Since
$|\mathcal{S}\setminus\mathcal{J}|=k-|\mathcal{J}|$, we have

$\displaystyle L({\bf\hat{X}}^{b_{\mathcal{J}}}[\mathcal{I}])$
$\displaystyle=$ $\displaystyle
t_{\mathcal{J}}-\sum_{l=|\mathcal{J}|+1}^{|\mathcal{I}|}\
\sum_{k=l}^{|\mathcal{I}|}(-1)^{k-l}\binom{k-|\mathcal{J}|}{k-l}\
\sum_{\begin{subarray}{c}\mathcal{J}\subset\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (52) $\displaystyle=$
$\displaystyle t_{\mathcal{J}}+\sum_{k=|\mathcal{J}|+1}^{|\mathcal{I}|}\
\sum_{l=|\mathcal{J}|+1}^{k}(-1)^{k-l-1}\binom{k-|\mathcal{J}|}{k-l}\
\sum_{\begin{subarray}{c}\mathcal{J}\subset\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (53) $\displaystyle=$
$\displaystyle t_{\mathcal{J}}+\
\sum_{k=|\mathcal{J}|+1}^{|\mathcal{I}|}(-1)^{k-|\mathcal{J}|}\
\sum_{\begin{subarray}{c}\mathcal{J}\subset\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (54)

where the sum on $l$ is performed using the binomial theorem. We conclude that
the lemma is also valid for $b_{\mathcal{J}}$ with the assumption on all
parent vectors $b^{\prime}$. By recursion we prove the lemma. ∎

In the remainder of the appendix, we will refer to row-type and co-occurrence
representations as being the histogram of distinct rows and the histogram of
co-occurrences, respectively. Both form a complete representation of the joint
distribution of a binary matrix.

##  Appendix B PMDH distribution

###  B.1 Marginal distribution of co-occurrences

######  Lemma B.1.

Consider a set of columns $\mathcal{I}$ of which we count the number of co-
occurrences ${\mathbf{c}}_{\mathcal{I}}$ after randomization. We have

${\mathbf{c}}_{\mathcal{I}}=\sum_{b=0}^{|\mathcal{I}|}{\mathbf{B}}\left(\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\binom{k}{b}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}},q^{b}p^{|\mathcal{I}|-b}\right)$
(55)

where ${\mathbf{B}}(n,p)$ signifies a Binomial random variable of parameters
$n$ and $p$.

###### Proof.

By virtue of the independence of the randomization on each row, each segment
of ${\bf\hat{X}}$ with the same probability of yielding a co-occurrence will
contribute one Bernoulli experiment in a binomial random variable. This
experiment will be positive if all positive bits remain positive (probability
$=q^{|\mathcal{J}|}$) and, independently, all the remaining bits are flipped
up by the randomization (probability=$p^{|\mathcal{I}|-|\mathcal{J}|}$). The
total probability of a row to contribute to ${\mathbf{c}}_{\mathcal{I}}$ is
therefore $q^{b}p^{|\mathcal{I}|-b}$. A given $b\in[0,|\mathcal{I}|]$
corresponds to several sets $\mathcal{J}$ of positive bits in a row, all with
$|\mathcal{J}|=b$. Using the lemma A.1 for each of those yields

$\displaystyle\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{I}\\\
|\mathcal{J}|=b\end{subarray}}\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ $\displaystyle=$
$\displaystyle\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{S}\\\
|\mathcal{J}|=b\end{subarray}}1$ (56) $\displaystyle=$
$\displaystyle\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}\binom{k}{b}$ (57)

∎

###  B.2 Proof of theorem 3.1

###### Proof.

By the linearity of the expected value, and using $E({\mathbf{B}}(n,p))=np$

$\displaystyle E({\mathbf{c}}_{\mathcal{I}})$ $\displaystyle=$
$\displaystyle\sum_{b=0}^{|\mathcal{I}|}\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\binom{k}{b}q^{b}p^{|\mathcal{I}|-b}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (58) $\displaystyle=$
$\displaystyle\sum_{k=0}^{|\mathcal{I}|}\sum_{b=0}^{k}(-1)^{k-b}\binom{k}{b}q^{b}p^{|\mathcal{I}|-b}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (59) $\displaystyle=$
$\displaystyle\sum_{k=0}^{|\mathcal{I}|}\left(\sum_{b=0}^{k}\binom{k}{b}q^{b}(-p)^{k-b}\right)p^{|\mathcal{I}|-k}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (60) $\displaystyle=$
$\displaystyle\sum_{k=0}^{|\mathcal{I}|}(q-p)^{k}p^{|\mathcal{I}|-k}\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (61)

∎

###  B.3 Variance between arbitrary co-occurrences in SBS

######  Lemma B.2.

Consider a SBS of composed of $N$ experiments and $D=|\mathcal{I}|$ parallel
sequences $\mathcal{I}$ of corresponding success probabilities $\\{p_{i}\ |\
i\in\mathcal{I}\\}$. Consider two arbitrary column sets $\mathcal{J}$ and
$\mathcal{J}^{\prime}$. The covariance between the co-occurrences counts
${\mathbf{c}}_{\mathcal{J}}$ and ${\mathbf{c}}_{\mathcal{J}}^{\prime}$ is

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$ $\displaystyle
Np_{\mathcal{J}\cup\mathcal{J}^{\prime}}(1-p_{\mathcal{J}_{\cap}})$ (62)

with

$p_{\mathcal{S}}=\prod_{j\in\mathcal{S}}p_{j}$ (63)

###### Proof.

The sum of successes and co-occurring successes in SBS can be modeled using a
multinomial distribution in the row-type representation. Indeed, lemme A.1
allows to compute how many rows of a given type is present in a dataset
characterized by a co-occurrence histogram. Here the former representation is
interpreted as a given output for set of distinguishable synchronous Bernoulli
experiments (success set $\mathcal{J}$), and the latter is the partial PMDH on
the corresponding SBS segment ${\bf\hat{X}}^{b}[{\mathcal{I}}]$. We therefore
define the change of distribution representation

$\tilde{\mathbf{c}}_{\mathcal{J}}=\sum_{k=|\mathcal{J}|}^{|{\mathcal{I}}|}(-1)^{k-|\mathcal{J}|}\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{S}\subseteq{\mathcal{I}}\\\
|\mathcal{S}|=k\end{subarray}}{\mathbf{c}}_{\mathcal{S}}$ (64)

The probability that the SBS succeeds on the indexes $\mathcal{J}$ and fails
elsewhere is

$\tilde{p}_{\mathcal{J}}=\prod_{j\in\mathcal{J}}p_{j}\prod_{i\in{\mathcal{I}}\setminus\mathcal{J}}(1-p_{i})$
(65)

Because each row can be considered as one event and that each row is
independently randomized, we have that the joint probability distribution to
observe any given row follows a multinomial distribution

$P(\tilde{\mathbf{c}}_{\mathcal{J}}=\tilde{c}_{\mathcal{J}}\ \forall\
\mathcal{J}\subseteq{\mathcal{I}})\sim\mathrm{Mutinomial}(\\{\tilde{c}_{\mathcal{J}}\
\forall\ \mathcal{J}\subseteq{\mathcal{I}}\\};\\{\tilde{p}_{\mathcal{J}}\
\forall\ \mathcal{J}\subseteq{\mathcal{I}}\\})$ (66)

We can then compute the full covariance matrix of
$\mathbf{\tilde{C}}=\\{\tilde{\mathbf{c}}_{\mathcal{J}}\ \forall\
\mathcal{J}\subseteq{\mathcal{I}}\\}$ with

$V(\mathbf{\tilde{C}})=N\left(\mathrm{diag}(\tilde{P})-\tilde{P}\tilde{P}^{T}\right)$
(67)

where $\tilde{P}=\\{\tilde{p}_{\mathcal{J}}\ \forall\
\mathcal{J}\subseteq{\mathcal{I}}\\}$. The invert transformation from the row-
type form back to the co-occurrence form is

${\mathbf{c}}_{\mathcal{J}}=\sum_{\mathcal{J}\subseteq\mathcal{S}\subseteq{\mathcal{I}}}\tilde{\mathbf{c}}_{\mathcal{S}}$
(68)

so that by linearity of the covariance operator we can write

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{S}\subseteq{\mathcal{I}}\\\
\mathcal{J}^{\prime}\subseteq\mathcal{S}^{\prime}\subseteq{\mathcal{I}}\end{subarray}}V(\tilde{\mathbf{c}}_{\mathcal{S}},\tilde{\mathbf{c}}_{\mathcal{S}^{\prime}})$
(69)

which become after substituting with the result of lemma B.2

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}}){}$
$\displaystyle=$ $\displaystyle
N\sum_{\begin{subarray}{c}\mathcal{J}\subseteq\mathcal{S}\subseteq{\mathcal{I}}\\\
\mathcal{J}^{\prime}\subseteq\mathcal{S}^{\prime}\subseteq{\mathcal{I}}\end{subarray}}\tilde{p}_{\mathcal{S}}\left(\delta(\mathcal{S},\mathcal{S}^{\prime})-\tilde{p}_{\mathcal{S}^{\prime}}\right)$
(70)

The probabilities 65 behave as a vector in the row-type to co-occurrence
transformation 68

$p_{\mathcal{J}}=\sum_{\mathcal{J}\subseteq\mathcal{S}\subseteq{\mathcal{I}}}\tilde{p}_{\mathcal{S}}$
(71)

This can be understood easily by the fact that the sum of all
$\tilde{p}_{\mathcal{S}}/p_{\mathcal{J}}$, representing the probability of any
variation of columns outside of $\mathcal{J}$, sum to 1. We obtain

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})/N$
$\displaystyle=$ $\displaystyle-
p_{\mathcal{J}}p_{\mathcal{J}^{\prime}}+\sum_{\mathcal{J}\cup\mathcal{J}^{\prime}\subseteq\mathcal{S}\subseteq{\mathcal{I}}}\tilde{p}_{\mathcal{S}}$
(72)

The final result is a consequence of equation 71 and the event independence

$p_{\mathcal{J}}p_{\mathcal{J}^{\prime}}=p_{\mathcal{J\cup
J}^{\prime}}(1-p_{\mathcal{J\cap J}^{\prime}}).$ (74)

∎

###  B.4 Covariance matrix of PMDH

####  B.4.1 Proof of theorem 3.2

###### Proof.

By the independence of each row in ${\bf\hat{X}}$, the linearity of the
variance, and using $V({\mathbf{B}}(n,p))=np(1-p)$ we have

$\displaystyle V({\mathbf{c}}_{\mathcal{I}})$ $\displaystyle=$
$\displaystyle\sum_{b=0}^{|\mathcal{I}|}\sum_{k=b}^{|\mathcal{I}|}(-1)^{k-b}\binom{k}{b}q^{b}p^{|\mathcal{I}|-b}(1-q^{b}p^{|\mathcal{I}|-b})\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (75) $\displaystyle=$
$\displaystyle\sum_{k=0}^{|\mathcal{I}|}\sum_{b=0}^{k}\left(\binom{k}{b}q^{b}(-p)^{k-b}p^{|\mathcal{I}|-k}-q^{2b}(-p^{2})^{k-b}p^{2(|\mathcal{I}|-k)}\right)\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (76) $\displaystyle=$
$\displaystyle\sum_{k=0}^{|\mathcal{I}|}(q-p)^{k}p^{|\mathcal{I}|-k}(1-(q+p)^{k}p^{|\mathcal{I}|-k})\sum_{\begin{subarray}{c}\mathcal{S}\subseteq\mathcal{I}\\\
|\mathcal{S}|=k\end{subarray}}t_{\mathcal{S}}$ (77)

which leads directly to the result. ∎

####  B.4.2 Proof of theorem 3.3

###### Proof.

For convenience, we introduce the notation
$\mathcal{J}_{\cup}=\mathcal{J}\cup\mathcal{J}^{\prime}$. Theorem B.2 provides
an expression for the covariance matrix in a segment
${\bf\hat{X}}^{b}[{\mathcal{J}_{\cup}}]$ of ${\bf\hat{X}}$. By the
independence of each row in ${\bf\hat{X}}$, and thanks to the covariance
linearity, we have

$V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})=\sum_{\mathcal{B}\subseteq\mathcal{J}_{\cup}}V({\mathbf{c}}^{(b)}_{\mathcal{J}},{\mathbf{c}}^{(b)}_{\mathcal{J}^{\prime}})$
(78)

where the sum is performed on all possible $b$ ($\sim\mathcal{B}$) and
${\mathbf{c}}^{(b)}$ refers to the partial PMDH computed on the segment
defined by $b$. We have using $L$ to count the number of elements in
${\bf\hat{X}}^{b}[{\mathcal{J}_{\cup}}]$

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{B}\subseteq\mathcal{J}_{\cup}}L({\bf\hat{X}}^{b}[{\mathcal{J}_{\cup}}])\
p_{\mathcal{J}\cup\mathcal{J}^{\prime}}(1-p_{\mathcal{J}_{\cap}})$ (79)

We pause to interpret the above equation. Choosing $b$ is equivalent to
choosing a set $\mathcal{B}$ defining the ground truth from which a row has
been generated. $\mathcal{J}$ and $\mathcal{J}^{\prime}$ however represent the
_output_ of the randomization. They define the actual row in ${\bf\hat{X}}$,
and $\tilde{p}_{\mathcal{J}}$ is the probability to obtain it from a ground
truth $b$. Therefore $\tilde{p}_{\mathcal{J}}$ depends implicitly on $b$

$p_{\mathcal{J}}(\mathcal{B})=q^{|\mathcal{J}\cap\mathcal{B}|}p^{|\mathcal{J}\setminus\mathcal{B}|}$
(80)

We use lemma A.1 and introduce the shorthand notation
$\mathcal{J}_{\cap}=\mathcal{J}\cap\mathcal{J}^{\prime}$ to obtain

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{B}\subseteq{\mathcal{J}_{\cup}}}p_{\mathcal{J}_{\cup}}(\mathcal{B})(1-p_{\mathcal{J}_{\cap}}(\mathcal{B}))\sum_{\mathcal{B}\subseteq\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}(-)^{|\mathcal{S}|-|\mathcal{B}|}t_{\mathcal{S}}$
(81) $\displaystyle=$
$\displaystyle\sum_{\mathcal{B}\subseteq{\mathcal{J}_{\cup}}}q^{|\mathcal{B}|}p^{|\mathcal{J}_{\cup}|-|\mathcal{B}|}(1-p_{\mathcal{J}_{\cap}}(\mathcal{B}))\sum_{\mathcal{B}\subseteq\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}(-)^{|\mathcal{S}|-|\mathcal{B}|}t_{\mathcal{S}}$
(82) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\sum_{\mathcal{B}\subseteq\mathcal{S}}q^{|\mathcal{B}|}(-p)^{|\mathcal{S}|-|\mathcal{B}|}(1-p_{\mathcal{J}_{\cap}}(\mathcal{B}))$
(83) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left((q-p)^{|\mathcal{S}|}-\sum_{\mathcal{B}\subseteq\mathcal{S}}q^{|\mathcal{B}|+|\mathcal{J}_{\cap}\cap\mathcal{B}|}(-p)^{|\mathcal{S}|-|\mathcal{B}|}p^{|\mathcal{J}_{\cap}\setminus\mathcal{B}|}\right)$
(84) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left((q-p)^{|\mathcal{S}|}-\sum_{\mathcal{B}^{\prime}\subseteq(\mathcal{S}\cap\mathcal{J}_{\cap})}q^{|\mathcal{B}^{\prime}|}p^{|\mathcal{J}_{\cap}|-|\mathcal{B}^{\prime}|}\sum_{\mathcal{B}^{\prime\prime}\subseteq(\mathcal{S}\setminus\mathcal{J}_{\cap})}q^{|\mathcal{B}|}(-p)^{|\mathcal{S}|-|\mathcal{B}|}\right)$
(85)

where in the last step we split $\mathcal{B}$ into the part intersecting
$\mathcal{B}$, $\mathcal{B}^{\prime}$ and the part excluding $\mathcal{B}$,
$\mathcal{B}^{\prime\prime}$. We have
$\mathcal{B}=\mathcal{B}^{\prime}+\mathcal{B}^{\prime\prime}$ and
$|\mathcal{B}|=|\mathcal{B}^{\prime}|+|\mathcal{B}^{\prime\prime}|$. Using the
binomial theorem twice, we obtain the result

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}},{\mathbf{c}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}$ (87)
$\displaystyle\left((q-p)^{|\mathcal{S}|}-p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}\sum_{\mathcal{B}^{\prime}\subseteq(\mathcal{S}\cap\mathcal{J}_{\cap})}q^{2|\mathcal{B}^{\prime}|}(-p^{2})^{|\mathcal{S}\cap\mathcal{J}_{\cap}|-|\mathcal{B}^{\prime}|}\sum_{\mathcal{B}^{\prime\prime}\subseteq(\mathcal{S}\setminus\mathcal{J}_{\cap})}q^{|\mathcal{B}^{\prime\prime}|}(-p)^{|\mathcal{S}\setminus\mathcal{J}_{\cap}|-|\mathcal{B}^{\prime\prime}|}\right)$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left((q-p)^{|\mathcal{S}|}-p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}(q^{2}-p^{2})^{|\mathcal{S}\cap\mathcal{J}_{\cap}|}(q-p)^{|\mathcal{S}\setminus\mathcal{J}_{\cap}|}\right)$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}t_{\mathcal{S}}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\
(q-p)^{|\mathcal{S}|}\left(1-p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}(q+p)^{|\mathcal{S}\cap\mathcal{J}_{\cap}|}\right)$
(88)

∎

###  B.5 Example of M matrix of PMDH for 3 columns and it’s inverse

For three columns, the maximum co-occurrence order is 3. The first block is
composed of the three column counts ${\mathbf{c}}_{1},{\mathbf{c}}_{2},$ and
${\mathbf{c}}_{3}$. The second block is also three-dimensional, with
${\mathbf{c}}_{12},{\mathbf{c}}_{13}$, and ${\mathbf{c}}_{23}$. The last block
has only one dimension, ${\mathbf{c}}_{123}$ We then have

$M=\begin{bmatrix}(q-p)\mathbb{I}&0&0\\\
p(q-p)\mathbb{K}^{(21)}&(q-p)^{2}\mathbb{I}&0\\\ p^{2}(q-p)\
\mathbb{K}^{(31)}&p(q-p)^{2}\ \mathbb{K}^{(32)}&(q-p)^{3}\\\ \end{bmatrix}$
(89)

where here we denote the 3 by 3 identity matrix $\mathbb{I}$ for simplicity.
We also have

$\mathbb{K}^{(21)}=\begin{bmatrix}1&1&0\\\ 1&0&1\\\
0&1&1\end{bmatrix}\qquad\rm{and}\qquad\mathbb{K}^{(31)}=\mathbb{K}^{(32)}=[1,1,1].$
(90)

The inverse $M^{-1}$ is

$M^{-1}=\frac{1}{(q-p)^{3}}\left(\begin{bmatrix}(q-p)^{2}\mathbb{I}&0&0\\\
-p(q-p)\mathbb{K}^{(21)}&(q-p)\mathbb{I}&0\\\ p^{2}\ \mathbb{K}^{(31)}&-p\
\mathbb{K}^{(32)}&1\\\ \end{bmatrix}\right)$ (91)

As an example, we can verify that the $(31)$ element of the product $MM^{-1}$
vanishes

$\displaystyle p^{2}(q-p)\ \mathbb{K}^{(31)}(q-p)^{-1}\mathbb{I}-p(q-p)^{2}\
\mathbb{K}^{(32)}p(q-p)^{-2}\mathbb{K}^{(21)}+(q-p)^{3}p^{2}(q-p)^{-3}\
\mathbb{K}^{(31)}$ $\displaystyle=$ $\displaystyle 0$ (92) $\displaystyle
p^{2}\ \mathbb{K}^{(31)}-p^{2}\ \mathbb{K}^{(32)}\mathbb{K}^{(21)}+p^{2}\
\mathbb{K}^{(31)}$ $\displaystyle=$ $\displaystyle 0$ (93)

which is obtained because we have

$\mathbb{K}^{(32)}\mathbb{K}^{(21)}=2\mathbb{K}^{(31)}$ (94)

as demonstrated in the next section (lemma C.1).

##  Appendix C TMDH distribution

###  C.1 Proof of theorem 3.4

######  Lemma C.1.

For any $\alpha>\beta>\gamma$, we have

$\mathbb{K}^{(\alpha\beta)}\mathbb{K}^{(\beta\gamma)}=\binom{\alpha-\gamma}{\alpha-\beta}\mathbb{K}^{(\alpha\gamma)}$
(95)

###### Proof.

We introduce the arbitrary sets $\mathcal{A}$ and $\mathcal{C}$, of size
$\alpha$ and $\gamma$ respectively. We have

$\displaystyle(\mathbb{K}^{(\alpha\beta)}\mathbb{K}^{(\beta\gamma)})(\mathcal{A},\mathcal{C})$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{B},|\mathcal{B}|=\beta}\mathbb{K}^{(\alpha\beta)}(\mathcal{A},\mathcal{B})\mathbb{K}^{(\beta\gamma)}(\mathcal{B},\mathcal{C})$
(96) $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}\mathcal{C}\subset\mathcal{B}\subset\mathcal{A}\\\
|\mathcal{B}|=\beta\end{subarray}}1$ (97) $\displaystyle=$
$\displaystyle\binom{\alpha-\gamma}{\alpha-\beta}\mathbb{K}^{(\alpha\gamma)}(\mathcal{A},\mathcal{C})$
(98)

∎

We now turn to the proof of theorem 3.4.

###### Proof.

In order to obtain the TMDH normal parameters, we first need to invert the
matrix $M$. This can be done using the Neumann series:

$(A-B)^{-1}=(\mathbb{I}-A^{-1}B)^{-1}A^{-1}=\left(\sum_{k=0}^{\infty}(A^{-1}B)^{k}\right)A^{-1}$
(99)

In our case, we have

$A^{-1}=\sum_{\alpha=1}^{D}(q-p)^{-\alpha}\mathbb{I}^{(\alpha)}$ (100)

and

$A^{-1}B=-\sum_{\alpha=2}^{D}\sum_{\beta=1}^{\alpha-1}\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}$
(101)

This is a lower triangular matrix so that we are assured that the series 99
terminates. Using the lemma C.1 we obtain

$\displaystyle\left(\sum_{\alpha=2}^{D}\sum_{\beta=1}^{\alpha-1}\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}\right)^{2}$
$\displaystyle=$
$\displaystyle\sum_{\alpha=3}^{D}\sum_{\beta=2}^{\alpha-1}\sum_{\gamma=1}^{\beta-1}\left(\frac{p}{q-p}\right)^{\alpha-\gamma}\binom{\alpha-\gamma}{\alpha-\beta}\mathbb{K}^{(\alpha\gamma)}$
(102) $\displaystyle=$
$\displaystyle\sum_{\alpha=3}^{D}\sum_{\gamma=1}^{\alpha-2}\left(\frac{p}{q-p}\right)^{\alpha-\gamma}\mathbb{K}^{(\alpha\gamma)}\sum_{\beta^{\prime}=1}^{\alpha-\gamma-1}\binom{\alpha-\gamma}{\beta^{\prime}}$
(103) $\displaystyle=$
$\displaystyle\sum_{\alpha=3}^{D}\sum_{\beta=1}^{\alpha-2}(2^{\alpha-\beta}-2)\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}$
(104)

where we rename $\gamma$ to $\beta$ in the last step. Similarly

$\displaystyle\left(\sum_{\alpha=2}^{D}\sum_{\beta=1}^{\alpha-1}\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}\right)^{3}$
$\displaystyle=$
$\displaystyle\sum_{\alpha=4}^{D}\sum_{\beta=2}^{\alpha-2}\sum_{\gamma=1}^{\beta-1}(2^{\alpha-\beta}-2)\left(\frac{p}{q-p}\right)^{\alpha-\gamma}\binom{\alpha-\gamma}{\alpha-\beta}\mathbb{K}^{(\alpha\gamma)}$
(105) $\displaystyle=$
$\displaystyle\sum_{\alpha=4}^{D}\sum_{\gamma=1}^{\alpha-3}\left(\frac{p}{q-p}\right)^{\alpha-\gamma}\mathbb{K}^{(\alpha\gamma)}\sum_{\beta^{\prime}=1}^{\alpha-\gamma-2}(2^{\alpha-\gamma-\beta^{\prime}}-2)\binom{\alpha-\gamma}{\beta^{\prime}}$
(106)

By repeating the process we deduce that in general

$\left(\sum_{\alpha=2}^{D}\sum_{\beta=1}^{\alpha-1}\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}\right)^{k}=\sum_{\alpha=k+1}^{D}\sum_{\gamma=1}^{\alpha-k}\kappa^{(k)}_{\alpha-\gamma}\left(\frac{p}{q-p}\right)^{\alpha-\gamma}\mathbb{K}^{(\alpha\gamma)}$
(107)

where $\kappa^{(k)}_{\alpha\gamma}$ are combinatorial coefficients that can be
computed recursively

$\kappa^{(k+1)}_{\delta}=\sum_{\beta=1}^{\delta-k}\kappa^{(k)}_{\delta-\beta}\binom{\delta}{\beta}=\sum_{\beta^{\prime}=k}^{\delta-1}\kappa^{(k)}_{\beta^{\prime}}\binom{\delta}{\beta}$
(108)

Therefore $\kappa^{(k)}_{\delta}=0$ if $k>\delta$ and the highest non-
vanishing $k$ is $D-1$. Let us investigate those coefficients with
combinatorics reasoning.

The coefficients $\kappa^{(k)}_{\delta}$ appear also in the exponentiation of
the full binary matrix $\mathbb{K}$ composed of all the blocks
$\mathbb{K}^{(\alpha\gamma)}$ without factors. We can interpret
$\mathbb{K}^{k}$ as the matrix counting all the ways to connect any to subset
of the range $[0,D]$ using $k$ steps. If we consider
$\mathbb{K}^{k}(\mathcal{A},\mathcal{B})$, then the corresponding weight is
$\kappa^{(k)}_{\alpha-\beta}$ if $\mathcal{B}\subset\mathcal{A}$. This count
represents all the ways $\mathcal{A}\setminus\mathcal{B}$ can be divided into
$k$ partitions in a given order. Bare the ordering, this corresponds to the
Stirling numbers of the second kind $S(n,k)$ defined as

$S(n,k)=\frac{1}{k!}\sum_{i=0}^{k}(-)^{i}\binom{k}{i}(k-i)^{n}$ (109)

In our case, we have to multiply those factors to account for the permutations
of the partitions

$\kappa^{(k)}_{\delta}=k!S(\delta,k)$ (110)

Those numbers are listed as the Triangle of numbers in ref. [18]. They obey
the triangular relation

$\kappa^{(k)}_{\delta}=k\left(\kappa^{(k-1)}_{\delta-1}+\kappa^{(k)}_{\delta-1}\right)$
(111)

Substituting equation 107 in the series 99 and replacing $\gamma$ by $\beta$
we have

$\displaystyle M^{-1}$ $\displaystyle=$
$\displaystyle\left(\sum_{k=0}^{D-1}(-)^{k}\sum_{\alpha=k+1}^{D}\sum_{\beta=1}^{\alpha-k}\kappa^{(k)}_{\alpha-\beta}\left(\frac{p}{q-p}\right)^{\alpha-\beta}\mathbb{K}^{(\alpha\beta)}\right)A^{-1}$
(112) $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{D}(q-p)^{-\alpha}\left(\mathbb{I}^{(\alpha)}+\sum_{\beta=1}^{\alpha-1}\textrm{\large$\varkappa$}_{\alpha-\beta}\
p^{\alpha-\beta}\ \mathbb{K}^{(\alpha\beta)}\right)$ (113)

where

$\displaystyle\textrm{\large$\varkappa$}_{\delta}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\delta}(-)^{k}\kappa^{(k)}_{\delta}$ (114)

We can verify manually that the first $\textrm{\large$\varkappa$}_{\delta}$
are equal to $(-)^{\delta}$. Using the triangle relation 111 we have

$\displaystyle\textrm{\large$\varkappa$}_{\delta}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\delta}(-)^{k}k\left(\kappa^{(k-1)}_{\delta-1}+\kappa^{(k)}_{\delta-1}\right)$
(116) $\displaystyle=$
$\displaystyle(-)^{\delta}\delta\kappa^{\delta}_{\delta-1}+\kappa^{0}_{\delta-1}+\sum_{k=1}^{\delta-1}(-)^{k}\left(-(k+1)\kappa^{(k)}_{\delta-1}+k\kappa^{(k)}_{\delta-1}\right)$
(117) $\displaystyle=$ $\displaystyle-\textrm{\large$\varkappa$}_{\delta-1}$
(118)

with $\kappa^{\delta}_{\delta-1}=\kappa^{0}_{\delta-1}=0$

∎

###  C.2 Proof of theorems 3.5 and 3.6

###### proof of theorem 3.5.

$\displaystyle\mu_{\mathcal{I}}(C)$ $\displaystyle=$
$\displaystyle\left(M^{-1}(C-NP)\right)_{\mathcal{I}}$ (119) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}}M_{\mathcal{I}\mathcal{S}}^{-1}\left(c_{\mathcal{S}}-p^{|\mathcal{S}|}N\right)$
(120) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}}(q-p)^{-|\mathcal{I}|}(-p)^{|\mathcal{I}|-|\mathcal{S}|}\mathbb{K}^{(|\mathcal{I}||\mathcal{S}|)}(\mathcal{I},\mathcal{S})\left(c_{\mathcal{S}}-p^{|\mathcal{S}|}N\right)$
(121) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq\mathcal{I}}\frac{(-p)^{|\mathcal{I}|-|\mathcal{S}|}}{(q-p)^{|\mathcal{I}|}}\left(c_{\mathcal{S}}-p^{|\mathcal{S}|}N\right)$
(122) $\displaystyle=$
$\displaystyle({q-p})^{-|\mathcal{I}|}\sum_{\mathcal{S}\subseteq\mathcal{I}}(-p)^{|\mathcal{I}|-|\mathcal{S}|}c_{\mathcal{S}}$
(123)

where we use the binomial formula together with ${\mathbf{c}}_{\emptyset}=N$
in the last step. ∎

###### proof of theorem 3.6.

For any matrix element of $\hat{\Sigma}(\mu)$ we have using theorems 3.3 and
3.5

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}}(\mu),{\mathbf{c}}_{\mathcal{J}^{\prime}}(\mu))$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}\mu_{\mathcal{S}}\
(q-p)^{|\mathcal{S}|}\
p^{|\mathcal{J}_{\cup}|-|\mathcal{S}|}\left(1-(q+p)^{|\mathcal{S}\cap\mathcal{J}_{\cap}|}p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}\right)$
(124) $\displaystyle=$
$\displaystyle\sum_{\mathcal{S^{\prime}}\subseteq\mathcal{S}\subseteq{\mathcal{J}_{\cup}}}c_{\mathcal{S^{\prime}}}(-)^{|\mathcal{S}|-|\mathcal{S^{\prime}}|}p^{|\mathcal{J}_{\cup}|-|\mathcal{S^{\prime}}|}\left(1-(q+p)^{|\mathcal{S}\cap\mathcal{J}_{\cap}|}p^{|\mathcal{J}_{\cap}|-|\mathcal{S}\cap\mathcal{J}_{\cap}|}\right)$
(125)

Each term depends on $\mathcal{S}$ only through the sign $(-)^{|\mathcal{S}|}$
and the size of its overlap with the intersection $\mathcal{J}_{\cap}$
($|\mathcal{S}\cap\mathcal{J}_{\cap}|$). The summation on $\mathcal{S}$ can be
simplified using
$\mathcal{J}_{\Delta}=\mathcal{J}_{\cup}\setminus\mathcal{J}_{\cap}$

$\displaystyle\sum_{s_{\cap}=0}^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}\sum_{s_{\Delta}=0}^{|\mathcal{J}_{\Delta}\setminus\mathcal{S^{\prime}}|}(-)^{|\mathcal{S^{\prime}}|+s_{\cap}+s_{\Delta}}\binom{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}{s_{\cap}}\binom{|\mathcal{J}_{\Delta}\setminus\mathcal{S^{\prime}}|}{s_{\Delta}}\left(1-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|+s_{\cap}}p^{|\mathcal{J}_{\cap}|-|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|-s_{\cap}}\right)$
(127) $\displaystyle=$
$\displaystyle(-)^{|\mathcal{S^{\prime}}|}\mathbb{I}(\mathcal{J}_{\Delta}\setminus\mathcal{S^{\prime}},\emptyset)\sum_{s_{\cap}=0}^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}(-)^{s_{\cap}}\binom{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}{s_{\cap}}\left(1-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|+s_{\cap}}p^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|-s_{\cap}}\right)$
$\displaystyle=$
$\displaystyle(-)^{|\mathcal{S^{\prime}}|}\mathbb{I}(\mathcal{J}_{\Delta}\setminus\mathcal{S^{\prime}},\emptyset)\left(\mathbb{I}(\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}},\emptyset)-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}(-q)^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}\right)$
(128)

so that

$\displaystyle
V({\mathbf{c}}_{\mathcal{J}}(\mu),{\mathbf{c}}_{\mathcal{J}^{\prime}}(\mu))$
$\displaystyle=$
$\displaystyle\sum_{\mathcal{S^{\prime}}\subseteq{\mathcal{J}_{\cup}}}c_{\mathcal{S^{\prime}}}p^{|\mathcal{J}_{\cup}|-|\mathcal{S^{\prime}}|}\mathbb{I}(\mathcal{J}_{\Delta}\setminus\mathcal{S^{\prime}},\emptyset)\left(\mathbb{I}(\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}},\emptyset)-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}(-q)^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}\right)$
(130) $\displaystyle=$
$\displaystyle\sum_{\mathcal{J}_{\Delta}\subseteq\mathcal{S^{\prime}}\subseteq{\mathcal{J}_{\cup}}}c_{\mathcal{S^{\prime}}}p^{|\mathcal{J}_{\cup}|-|\mathcal{S^{\prime}}|}\left(\mathbb{I}(\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}},\emptyset)-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}(-q)^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}\right)$

We can then develop equation 39 and inject the identity deduced above

$\displaystyle
V({\mathbf{t}}_{\mathcal{J}},{\mathbf{t}}_{\mathcal{J}^{\prime}})$
$\displaystyle=$
$\displaystyle(q-p)^{-(|\mathcal{J}|+|\mathcal{J}^{\prime}|)}\sum_{\begin{subarray}{c}\mathcal{B}\subseteq\mathcal{J}\\\
\mathcal{B}^{\prime}\subseteq\mathcal{J}^{\prime}\end{subarray}}(-p)^{|\mathcal{J}|+|\mathcal{J}^{\prime}|-|\mathcal{B}|-|\mathcal{B}^{\prime}|}V({\mathbf{c}}_{\mathcal{B}},{\mathbf{c}}_{\mathcal{B}^{\prime}})$
(131) $\displaystyle=$
$\displaystyle\left(\frac{-p}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J}^{\prime}|}\\!\\!\sum_{\begin{subarray}{c}\mathcal{B}\subseteq\mathcal{J}\\\
\mathcal{B}^{\prime}\subseteq\mathcal{J}^{\prime}\\\
\mathcal{B}_{\Delta}\subseteq\mathcal{S^{\prime}}\subseteq\mathcal{B}_{\cup}\end{subarray}}c_{\mathcal{S^{\prime}}}\
(-)^{|\mathcal{B}|+|\mathcal{B^{\prime}}|}\frac{\mathbb{I}(\mathcal{B}_{\cap}\setminus\mathcal{S^{\prime}},\emptyset)-(q+p)^{|\mathcal{B}_{\cap}\cap\mathcal{S^{\prime}}|}(-q)^{|\mathcal{B}_{\cap}\setminus\mathcal{S^{\prime}}|}}{p^{2|\mathcal{B}_{\cup}|-|\mathcal{B}_{\Delta}|-|\mathcal{B}_{\cup}|+|\mathcal{S^{\prime}}|}}$
(134) $\displaystyle=$
$\displaystyle\left(\frac{-p}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J}^{\prime}|}\\!\\!\sum_{\mathcal{S^{\prime}}\subseteq\mathcal{J}_{\cup}}c_{\mathcal{S^{\prime}}}\sum_{x=0}^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}\sum_{y=0}^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}(-)^{y+|\mathcal{S^{\prime}}\setminus\mathcal{J}_{\cap}|}2^{y}$
$\displaystyle\qquad\qquad\qquad\quad\binom{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}{x}\binom{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}{y}\frac{\delta^{0}_{x}-(q+p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|-y}(-q)^{x}}{p^{|\mathcal{S^{\prime}}|+x+|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|-y}}$
$\displaystyle=$
$\displaystyle\left(\frac{-p}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J}^{\prime}|}\\!\\!\sum_{\mathcal{S^{\prime}}\subseteq\mathcal{J}_{\cup}}c_{\mathcal{S^{\prime}}}\sum_{x=0}^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}(-)^{|\mathcal{S^{\prime}}\setminus\mathcal{J}_{\cap}|}\binom{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}{x}\frac{\delta^{0}_{x}(1-2p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}-(q-p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}(-q)^{x}}{p^{|\mathcal{S^{\prime}}|+x+|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}}$
$\displaystyle=$
$\displaystyle\left(\frac{-p}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J}^{\prime}|}\\!\\!\sum_{\mathcal{S^{\prime}}\subseteq\mathcal{J}_{\cup}}c_{\mathcal{S^{\prime}}}(-)^{|\mathcal{S^{\prime}}\setminus\mathcal{J}_{\cap}|}\frac{(1-2p)^{|\mathcal{S^{\prime}}\cap\mathcal{J}_{\cap}|}-(-p)^{-|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}(q-p)^{|\mathcal{J}_{\cap}\cap\mathcal{S^{\prime}}|+|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}}{p^{|\mathcal{S^{\prime}}|+|\mathcal{J}_{\cap}\cap\mathcal{S^{\prime}}|}}$
(135) $\displaystyle=$
$\displaystyle\left(\frac{1}{q-p}\right)^{|\mathcal{J}|+|\mathcal{J^{\prime}}|}\\!\\!\sum_{\mathcal{S^{\prime}}\subseteq\mathcal{J}_{\cup}}c_{\mathcal{S^{\prime}}}(-p)^{|\mathcal{J}_{\cup}|-|\mathcal{S^{\prime}}|}\left((-p)^{|\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}|}(1-2p)^{|\mathcal{J}_{\cap}\cap\mathcal{S^{\prime}}|}-(q-p)^{|\mathcal{J_{\cap}}|}\right)$
(136)

using that

$|\mathcal{B}|+|\mathcal{B}^{\prime}|=|\mathcal{B}_{\cup}|+|\mathcal{B}_{\cap}|=2|\mathcal{B}_{\cup}|-|\mathcal{B}_{\Delta}|$
(137)

that
$\mathcal{B}_{\cup}\setminus\mathcal{J}_{\cap}\subseteq\mathcal{B}_{\Delta}$,
and that
$\mathcal{J}_{\cup}\setminus\mathcal{S^{\prime}}=\mathcal{J}_{\cap}\setminus\mathcal{S^{\prime}}$.
∎

## References

  * [1] Tianhao Wang, Jeremiah Blocki, Ninghui Li, and Somesh Jha.  Locally differentially private protocols for frequency estimation.  In 26th USENIX Security Symposium (USENIX Security 17), pages 729–745, Vancouver, BC, August 2017. USENIX Association.
  * [2] C. Dwork.  The algorithmic foundations of differential privacy.  Theoretical Computer Science, 9(3-4):211–407, 2014.
  * [3] Shiva Prasad Kasiviswanathan, Homin K. Lee, Kobbi Nissim, Sofya Raskhodnikova, and Adam Smith.  What can we learn privately?  SIAM Journal on Computing, 40(3):793–826, 2011.
  * [4] Stanley L. Warner.  Randomized response: A survey technique for eliminating evasive answer bias.  Journal of the American Statistical Association, 60(309):63–69, 1965.  PMID: 12261830.
  * [5] Xingxing Xiong, Shubo Liu, Dan Li, Zhaohui Cai, and Xiaoguang Niu.  A comprehensive survey on local differential privacy.  Security and Communication Networks, 2020(8829523), 08 Oct. 2020\.
  * [6] Mengmeng Yang, Lingjuan Lyu, Jun Zhao, Tianqing Zhu, and Kwok-Yan Lam.  Local differential privacy and its applications: A comprehensive survey, 2020.
  * [7] Giulia Fanti, Vasyl Pihur, and Úlfar Erlingsson.  Building a rappor with the unknown: Privacy-preserving learning of associations and data dictionaries.  Proceedings on Privacy Enhancing Technologies, 2016(3):41 – 61, 01 Jul. 2016.
  * [8] Xuebin Ren, Chia-Mu Yu, Weiren Yu, Shusen Yang, Xinyu Yang, Julie A. McCann, and Philip S. Yu.  LoPub: High-Dimensional Crowdsourced Data Publication with Local Differential Privacy.  arXiv e-prints, page arXiv:1612.04350, December 2016.
  * [9] Graham Cormode, Tejas Kulkarni, and Divesh Srivastava.  Marginal release under local differential privacy.  In Proceedings of the 2018 International Conference on Management of Data, SIGMOD ’18, page 131–146, New York, NY, USA, 2018. Association for Computing Machinery.
  * [10] Zhikun Zhang, Tianhao Wang, Ninghui Li, Shibo He, and Jiming Chen.  Calm: Consistent adaptive local marginal for marginal release under local differential privacy.  In Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security, CCS ’18, page 212–229, New York, NY, USA, 2018. Association for Computing Machinery.
  * [11] Shaowei Wang, Liusheng Huang, Pengzhan Wang, Yiwen Nie, Hongli Xu, Wei Yang, Xiang-Yang Li, and Chunming Qiao.  Mutual information optimally local private discrete distribution estimation, 2016.
  * [12] Úlfar Erlingsson, Vasyl Pihur, and Aleksandra Korolova.  Rappor: Randomized aggregatable privacy-preserving ordinal response.  In Proceedings of the 2014 ACM SIGSAC Conference on Computer and Communications Security, CCS ’14, page 1054–1067, New York, NY, USA, 2014. Association for Computing Machinery.
  * [13] Raef Bassily and Adam Smith.  Local, private, efficient protocols for succinct histograms.  In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing, STOC ’15, page 127–135, New York, NY, USA, 2015. Association for Computing Machinery.
  * [14] Tianhao Wang, Ninghui Li, and Somesh Jha.  Locally differentially private heavy hitter identification, 2017.
  * [15] N. Holohan, D. J. Leith, and O. Mason.  Optimal differentially private mechanisms for randomised response.  IEEE Transactions on Information Forensics and Security, 12(11):2726–2735, 2017.
  * [16] Peter Kairouz, Keith Bonawitz, and Daniel Ramage.  Discrete distribution estimation under local privacy.  In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 2436–2444, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * [17] Jayadev Acharya, Ziteng Sun, and Huanyu Zhang.  Hadamard response: Estimating distributions privately, efficiently, and with little communication.  In Kamalika Chaudhuri and Masashi Sugiyama, editors, Proceedings of Machine Learning Research, volume 89 of Proceedings of Machine Learning Research, pages 1120–1129. PMLR, 16–18 Apr 2019.
  * [18] OEIS Foundation Inc.  The on-line encyclopedia of integer sequences.  http://oeis.org/A019538.