# Adaptive Online Bayesian Estimation of Frequency Distributions with Local
Differential Privacy

Soner Aydın  Sinan Yıldırım

###### Abstract

We propose a novel Bayesian approach for the adaptive and online estimation of
the frequency distribution of a finite number of categories under the local
differential privacy (LDP) framework. The proposed algorithm performs Bayesian
parameter estimation via posterior sampling and adapts the randomization
mechanism for LDP based on the obtained posterior samples. We propose a
randomized mechanism for LDP which uses a subset of categories as an input and
whose performance depends on the selected subset and the true frequency
distribution. By using the posterior sample as an estimate of the frequency
distribution, the algorithm performs a computationally tractable subset
selection step to maximize the utility of the privatized response of the next
user. We propose several utility functions related to well-known information
metrics, such as (but not limited to) Fisher information matrix, total
variation distance, and information entropy. We compare each of these utility
metrics in terms of their computational complexity. We employ stochastic
gradient Langevin dynamics for posterior sampling, a computationally efficient
approximate Markov chain Monte Carlo method. We provide a theoretical analysis
showing that (i) the posterior distribution targeted by the algorithm
converges to the true parameter even for approximate posterior sampling, and
(ii) the algorithm selects the optimal subset with high probability if
posterior sampling is performed exactly. We also provide numerical results
that empirically demonstrate the estimation accuracy of our algorithm where we
compare it with nonadaptive and semi-adaptive approaches under experimental
settings with various combinations of privacy parameters and population
distribution parameters.

Keywords: Local differential privacy, frequency distribution, adaptive
estimation, posterior sampling, stochastic gradient Langevin dynamics

###### Contents

  1. 1 Introduction
  2. 2 Related Literature
  3. 3 Problem definition and general framework
  4. 4 Constructing informative randomized response mechanisms
    1. 4.1 The randomly restricted randomized response (RRRR) mechanism
    2. 4.2 Subset selection for RRRR
      1. 4.2.1 Fisher information matrix
      2. 4.2.2 Entropy of randomized response
      3. 4.2.3 Total variation distance
      4. 4.2.4 Expected mean squared error
      5. 4.2.5 Probability of honest response
      6. 4.2.6 Semi-adaptive approach
    3. 4.3 Computational complexity of utility functions
  5. 5 Posterior sampling
    1. 5.1 Stochastic gradient Langevin dynamics
    2. 5.2 Gibbs sampling
  6. 6 Theoretical analysis
    1. 6.1 Convergence of the posterior distribution
    2. 6.2 Selecting the best subset
  7. 7 Numerical results
  8. 8 Conclusion
  9. A Appendix
    1. A.1 Proofs for LDP of RRRR
    2. A.2 Proofs about utility functions
    3. A.3 Proof for SGLD update
    4. A.4 Proofs for convergence and consistency results
      1. A.4.1 Preliminary results
      2. A.4.2 Convergence of the posterior distribution
      3. A.4.3 Convergence of the expected frequency

##  1 Introduction

Frequency distributions are encountered in many applications that involve
personal and private data. Suppose a type of sensitive information is
represented as a categorical random variable $X$ which has a population
distribution denoted by $\text{Cat}(\theta)$, where $\theta$ is a
$K$-dimensional probability vector.

In this paper, we propose an adaptive and online algorithm to estimate
$\theta$ in a Local Differential Privacy (LDP) framework, when $X$ is
unobserved and instead, we have access to a randomized response $Y$ derived
from $X$. In the LDP framework, a central aggregator receives
privatized/perturbed data from each user and uses this noisy data for
inferential tasks. In that sense, LDP differs from standard DP framework
(Dwork, , 2006) where the aggregator privatizes operations on the sensitive
dataset after it collects the sensitive data without noise. Hence LDP can be
said to provide a stricter form of privacy and is used in cases where the
aggregator may not be trustable (Kasiviswanathan et al., , 2011). This feature
of LDP makes it especially useful for users of social media platforms and
cloud computing services where the central aggregator is more prone to leakage
of personal information with or without the knowledge of the company that
aggregates the data. Below we give a more formal definition of $\epsilon$-LDP
as a property that concerns a randomized mechanism.

######  Definition 1 (Local differential privacy).

A randomized mechanism $\mathcal{M}:\mathcal{X}\mapsto\mathcal{Y}$ satisfies
$\epsilon$-LDP if the following inequality holds for any pairs of inputs
$x,x^{\prime}\in\mathcal{X}$, and for any output (response) $y\in\mathcal{Y}$:

$e^{-\epsilon}\leq\frac{\mathbb{P}(\mathcal{M}(x)=y)}{\mathbb{P}(\mathcal{M}(x^{\prime})=y)}\leq
e^{\epsilon}.$

The definition of LDP is almost the same as that of global DP. The main
difference is that, in the global DP, neighboring inputs $x,x^{\prime}$ are
two datasets that differ in only one record; whereas in LDP, $x,x^{\prime}$
pairs are two different data points from the same domain. In Definition 1,
$\epsilon\geq 0$ is the privacy parameter. A smaller $\epsilon$ value provides
stronger privacy.

One main challenge in most differential privacy settings is to decide on the
randomized mechanism. In the case of LDP, this is how an individual data point
$X$ should be randomized. For a given randomized algorithm, too little
randomization may not guarantee the privacy of individual users; whereas too
severe randomization deteriorates the utility of the output of the randomized
algorithm. Balancing these conflicting objectives (privacy vs utility) is the
main goal of the research on estimation under privacy constraints.

While sequential observation of the individual points in a dynamic/online
setting may make the estimation task under the LDP constraint harder, it may
offer an opportunity to adapt the randomized mechanism in time to improve the
estimation quality. Motivated by that, in this paper, we address the problem
of online Bayesian estimation of a frequency distribution ($\theta$) under
$\epsilon$-LDP, while at the same time choosing the randomization mechanism
adaptively so that the utility is improved continually in time.

##### Contribution:

This paper presents AdOBEst-LDP: a novel methodological framework for Adaptive
Online Bayesian Estimation with LDP. A flowchart diagram of AdOBEst-LDP is
given in Figure 1 to expose the reader to the main idea of the framework.
AdOBEst-LDP seeks to adapt its randomized response mechanism to past data to
generate the ‘most informative’ responses about $\theta$ while maintaining
$\epsilon$-LDP for all users. AdOBEst-LDP uses a randomized mechanism, called
_Randomly Restricted Randomized Response_ (RRRR), to generate randomized
responses under $\epsilon$-LDP. The RRRR mechanism is a generalization of the
standard randomized response mechanism in that it restricts the response to a
_subset_ of the sample space of $X$. The subset used in RRRR is decided such
that the sensitive information $X$ from the next individual is likely
contained in that subset. To ensure this, the subset selection step uses two
inputs: (i) a sample for $\theta$ drawn from the posterior distribution of
$\theta$ conditional on the past data, (ii) a utility function that scores the
informativeness of the randomized response obtained from RRRR when it is run
with a given subset. To that end, we propose several utility functions to
score the informativeness of the randomized response. The utility functions
are based on well-known tools from probability and statistics. The first
utility function is based on the Fisher information of the privatized response
(Alparslan and Yıldırım, , 2022; Lopuhaä-Zwakenberg et al., , 2022;
Steinberger, , 2024; Yıldırım, , 2024). We propose a second utility function
based on the entropy of the randomized response. We also propose two utility
functions based on total variation distance. A fifth utility function is
proposed based on the expected squared error of the Bayesian estimator of the
sensitive data given the randomized response. Finally, the last utility
function considered is the probability of honest response, i.e., $Y=X$. We
provide some insight into those utility functions both theoretically and
numerically. Moreover, we also provide a computational complexity analysis for
the proposed utility functions.

Figure 1: Flowchart diagram of AdOBest-LDP

We equip AdOBEst-LDP with an online posterior sampling for parameter
estimation. Standard MCMC methods, such as Gibbs sampling, have a computation
complexity quadratic in the number of individuals whose data have been
collected. As a remedy to this, similar to Mazumdar et al.,  (2020) we propose
a stochastic gradient Langevin dynamics (SGLD)-based algorithm to obtain
approximate posterior samples (Welling and Teh, , 2011). By working on subsets
of data, the SGLD scales in time.

The numerical experiments show that AdOBEst-LPD outperforms its non-adaptive
counterpart when run with SGLD for posterior sampling. The results also
suggest that the utility functions considered in this paper are robust and
perform well. The MATLAB code at https://github.com/soneraydin/AdOBEst_LDP can
be used to reproduce the results obtained in this paper.

Finally, we provide a theoretical analysis of AdOBEst-LDP. We show two key
results:

  1. (i)

The targeted posterior distribution conditional on the generated observations
by the adaptive scheme converges to the true parameter in probability in the
number of observations, $n$. This convergence result owes mainly to the
smoothness and a special form of concavity of the marginal log-likelihood
function of the randomized responses. Another key factor is that the second
moment of the sum up to time $n$ of the gradient of this log-marginal
likelihood is $\mathcal{O}(n)$.

  2. (ii)

If perfect sampling is used for the posterior sampling stage, the expected
frequency of the algorithm choosing the best subset (according to the utility
function) tends to $1$ as $n$ goes to $\infty$.

The theoretical results depend on fairly weak, realistic, and verifiable
assumptions.

##### Outline:

In Section 2, we discuss the earlier work related to ours. Section 3 presents
LDP and the frequency estimation problem and introduces AdOBEst-LDP as a
general framework. In Section 4, we delve deeper into the details of AdOBEst-
LDP by first presenting RRRR, the proposed randomized response mechanism, then
explaining how it chooses an ‘optimal’ subset of categories adaptively at each
iteration. Section 4 also presents the utility metrics considered for choosing
these subsets in this paper. In Section 5, we provide the details of the
posterior sampling methods considered in this paper, particularly SGLD. The
theoretical analysis of AdOBEst-LDP is provided in Section 6. Section 7
contains the numerical experiments. Finally, Section 8 provides some
concluding remarks. All the proofs of the theoretical results are given in the
appendices.

##  2 Related Literature

Estimating frequency distributions under the LDP setting has been an
increasingly popular research area in recent years. Along with its basic
application (estimation of discrete probabilities from locally privatized
data), it is also used for a wide range of other estimation and learning
purposes; to name a few, estimation of confidence intervals and confidence
sets for a population mean (Waudby-Smith et al., , 2023), estimation or
identification of heavy hitters (Zhu et al., , 2024) (Wang et al., 2024b, )
(Jia and Gong, , 2019), estimation of quantiles (Cormode and Bharadwaj, ,
2022), frequent itemset mining (Zhao et al., , 2023), estimation of degree
distribution in social networks (Wang et al., 2024a, ), distributed training
of graph neural networks with categorical features and labels (Bhaila et al.,
, 2024). The methods that are proposed for $\epsilon$-LDP frequency estimation
also form the basis of more complex inferential tasks (with some modifications
on these methods), such as the release of ‘marginals’ (contingency tables)
between multiple categorical features and their correlations, as in the work
of Cormode et al.,  (2018).

AdOBEst-LDP employs a modified version of the Generalized Randomized Response
mechanism (GRR) to generate the randomized responses. Known also as
$k$-randomized response and direct encoding, this LDP mechanism (or protocol)
is the generalization of the binary randomized response mechanism for the
privatization of categorical data when the number of categories is $K>2$. The
basic version of GRR is based on simply flipping the category of a given
response with probability $\frac{1}{e^{\epsilon}+K-1}$ or reporting its true
category with probability $\frac{e^{\epsilon}}{e^{\epsilon}+K-1}$. This is a
well-studied mechanism in the DP literature, and the statistical properties of
its basic version (such as its estimation variance) can be found in the works
by (Wang et al., , 2017) and (Wang et al., , 2020). For convenience and unity
of terminology, we will call this basic approach as standard randomized
response (SRR) mechanism throughout the paper. Unlike the basic version of
SRR, the perturbation probability of responses in our algorithm changes
adaptively, depending on the cardinality of the selected subset of categories
(which we explain in detail in Section 4) for the privatization of real $X$,
and the cardinality of its complementary set.

The use of information metrics as utility functions in LDP protocols has been
an active line of research in recent years. In the work of Kairouz et al.,
(2016), information metrics like $f$-divergence and mutual information are
used for selecting optimal LDP protocols. In the same vein, Steinberger,
(2024) uses Fisher Information as the utility metric for finding a nearly
optimal LDP protocol for the frequency estimation problem, and Lopuhaä-
Zwakenberg et al.,  (2022) uses it for comparing the utility of various LDP
protocols for estimation of frequency distributions and finding the optimal
one. In these works, the mentioned information metrics are used statically,
i.e., to choose a protocol once and for all, for a given estimation task. The
approaches in these works suffer from computational complexity for large
values of $K$, because the search space for optimal protocols there grows in
the order of $2^{K}$. In some other works, such as Wang et al.,  (2016), a
randomly sampled subset of size $k\leq K$ is used to improve the efficiency of
this task, where the optimal $k$ is determined by maximizing the mutual
information between real data and the privatized data. However, this approach
is also static, the optimal subset size $k$ is selected only once, and the
optimization procedure only determines $k$ and not the subset itself. Unlike
those static approaches, AdOBEst-LDP dynamically encompasses the information
metric (such as the Fisher Information matrix and the other alternatives in
Section 4.2) to select the optimal subset for each time step. In addition, in
the subset step of AdOBEst-LDP, only $K$ candidate subsets are compared in
terms of their utilities at each iteration, enabling computational
traceability. This way of tackling the problem requires computing the given
information metric for only $K$ times at each iteration. We will provide
further details of this approach and elaborate on its computational advantages
in Section 4.2.

Another use of the Fisher Information in the LDP literature is for bounding
the estimation error for a given LDP protocol. For example, Barnes et al.,
(2020) uses Fisher Information inside van Trees inequality, the Bayesian
version of the Cramér-Rao bound (Gill and Levit, , 1995), for bounding the
estimation error of various LDP protocols for Gaussian mean estimation and
frequency estimation. Again, their work provides rules for choosing optimal
protocols for a given $\epsilon$ in a static way. As a similar example,
Acharya et al.,  (2023) derives a general information contraction bound for
parameter estimation problems under LDP and shows its relation to van Trees
inequality as its special case. To our knowledge, our approach is the first
one that adaptively uses a utility metric to dynamically update the inner
workings of an LDP protocol for estimating frequency distributions.

The idea of building adaptive mechanisms for improved estimation under the LDP
has been studied in the literature, although the focus and methodology of
those works differ from ours. For example, Joseph et al.,  (2019) proposed a
two-step adaptive method to estimate the unknown mean parameter of data from
Gaussian distribution. In this method, the users are split into two groups, an
initial mean estimate is obtained from the perturbed data of the first group,
and the data from the second group is transformed adaptively according to that
initial estimate. Similarly, Wei et al.,  (2024) proposed another two-step
adaptive method for the mean estimation problem, in which, the aggregator
first computes a rough distribution estimate from the noisy data of a small
sample of users, which is then used for adaptively adjusting the amount of
perturbation for the data of remaining users. While Joseph et al.,  (2019);
Wei et al.,  (2024) consider a two-stage method, AdOBEst-LDP seeks to adapt
continually by updating its LDP mechanism each time an individual’s
information is collected. Similar to our work, Yıldırım,  (2024) has recently
proposed an adaptive LDP mechanism for online parameter estimation for
continuous distributions. The LDP mechanism of Yıldırım,  (2024) contains a
truncation step with boundaries adapted to the estimate from the past data
according to a utility function based on the Fisher information.
Unfortunately, the parameter estimation step of Yıldırım,  (2024) does not
scale well in time. Differently from Yıldırım,  (2024), AdOBEst-LDP focuses on
frequency distributions, considers several other utility functions to update
its LDP mechanism, employs a scalable parameter estimation step, and its
performance is backed up with theoretical results.

To avoid potential confusion, we remark that other notions of ‘adaptiveness’
are also used in the DP literature. For example, Du et al.,  (2021) used an
adaptive scheme to partition the domain of non-categorical features for a
given dataset, to frugally use the privacy budget for responding to range
queries. Similarly, Qardaji et al.,  (2013) proposed an adaptive method to
partition geospatial data into grids, to answer range queries. As another
example of using a different notion of adaptiveness, da Costa Filho and
Machado,  (2023) proposed an adaptive selection scheme to choose one of the
two well-known frequency estimation protocols, namely GRR and OLH (optimal
local hash), as a step of their main algorithm. Similarly, Zhang et al.,
(2018) proposed an adaptive scheme to choose between GRR and OUE (optimized
unary encoding) protocols for a given task of marginal release. Li and Miklau,
(2012) proposed an adaptive strategy to distribute the tasks of answering
batch queries (which they call workload) in $(\epsilon,\delta)$-DP setting, to
use a given privacy budget effectively. McKenna et al.,  (2020) also uses the
same notion of adaptiveness for answering batch queries, but under the
$\epsilon$-LDP setting.

##  3 Problem definition and general framework

Suppose we are interested in a discrete probability distribution $\mathcal{P}$
of a certain form of sensitive categorical information
$X\in\mathcal{X}:=\\{1,\ldots,K\\}$ of individuals in a population. Hence,
$\mathcal{P}$ is a categorical distribution $\text{Cat}(\theta^{\ast})$ with a
probability vector

$\theta^{\ast}:=(\theta^{\ast}_{1},\ldots,\theta^{\ast}_{K})\in\Delta,$

where $\Delta$ is the $(K-1)$-dimensional probability simplex,

$\Delta:=\left\\{\theta\in\mathbb{R}^{K}:\sum_{k=1}^{K}\theta_{k}=1\text{ and
}\theta_{k}\geq 0\text{ for }k\in\\{1,\ldots,K\\}\right\\}.$

We assume a setting where the individuals’ sensitive data are collected
_privately_ and _sequentially in time_. The privatization is performed via a
randomized algorithm that, upon taking a category index in $\\{1,\ldots,K\\}$
as an input, returns a random category index in $\\{1,\ldots,K\\}$ such that
the whole data collection process is $\epsilon$-LDP. (See Definition 1 above.)
Let $X_{t}$ and $Y_{t}$ be the private information and randomized responses of
individual $t$, respectively. According to Definition 1 for LDP, the following
inequality must be satisfied for all triples
$(x,x^{\prime},y)\in\mathcal{X}^{3}$ for the randomized mechanism to be
$\epsilon$-LDP.

$\mathbb{P}(Y_{t}=y|X_{t}=x)\leq
e^{\epsilon}\mathbb{P}(Y_{t}=y|X_{t}=x^{\prime}).$ (1)

The inferential goal is to estimate $\theta^{\ast}$ sequentially based on the
responses $Y_{1},Y_{2},\ldots$, and the mechanisms
$\mathcal{M}_{1},\mathcal{M}_{2},\ldots$ that are used to generate those
responses. Specifically, Bayesian estimation is considered, whereby the target
is the posterior distribution, denoted by
$\Pi(\mathrm{d}\theta|Y_{1:n},\mathcal{M}_{1:n})$, given a prior probability
distribution with pdf $\eta(\theta)$ on $\Delta$.

This paper concerns online Bayesian estimation of $\theta$ while adapting the
randomized mechanism to improve the estimation utility continually. We propose
a general framework called Adaptive Online Bayesian Estimation with LDP
(AdOBEst-LDP), in which the randomized mechanism at time $t+1$ is adapted to
the data collected until time $t$. AdOBEst-LDP is outlined in Algorithm 1.

1 Initialization: Start with an initial estimator
$\Theta_{0}=\theta_{\text{init}}$.

2for  _$t=1,2,\ldots$_ do

3 Step 1: Adapting the LDP mechanism: Based on $\Theta_{t-1}$, determine the
$\epsilon$-LDP mechanism $\mathcal{M}_{t}$ for the next individual according
to a utility metric.

4 Step 2: LDP response generation The sensitive information $X_{t}$ of
individual $t$ is shared as $Y_{t}$ using the $\epsilon$-LDP mechanism
$\mathcal{M}_{t}$.

5 Step 3: Draw a sample (approximately) from the posterior distribution
$\Theta_{t}\sim\Pi(\cdot|Y_{1:t},\mathcal{M}_{1:t}).$

Algorithm 1 AdOBEst-LDP: Adaptive Online Bayesian Estimation with LDP

Algorithm 1 is fairly general and it does not describe how to choose the
$\epsilon$-LDP mechanism $\mathcal{M}_{t}$ at time $t$, nor does it provide
the details of the posterior sampling. However, it is still worth making some
critical observations about the nature of the algorithm. Firstly, note that
the selection of the $\epsilon$-LDP mechanism in Step 1 relies on the
posterior sample $\Theta_{t}$, which serves as an _estimator_ of the true
parameter $\theta^{\ast}$. As we shall see in Section 4, at Step 1 the ‘best’
$\epsilon$-LDP mechanism is chosen from a set of candidate LDP mechanisms
according to a utility function. This step is relevant only when $\Theta_{t}$
is a reliable estimator of $\theta^{\ast}$. In other words, Step 1 ‘exploits’
the estimator $\Theta_{t}$. Moreover, the random nature of posterior sampling
prevents having too much confidence in the current estimator $\Theta_{t}$ and
enables a certain degree of ‘exploration’. In conclusion, Algorithm 1 utilizes
an ‘exploration-exploitation’ approach reminiscent of reinforcement learning.
In particular, posterior sampling in Step 3 suggests a strong parallelism
between AdOBEst-LDP and the well-known exploration-exploitation approach
called Thompson sampling (Russo et al., , 2018).

The details of Steps 1-2 and Step 3 of Algorithm 1 are given in Sections 4 and
5, respectively.

##  4 Constructing informative randomized response mechanisms

In this section, we describe Step 1 of AdOBEst-LDP in Algorithm 1 where the
$\epsilon$-LDP mechanism $\mathcal{M}_{t}$ is selected at time $t$ based on
the posterior sample $\Theta_{t}$. For ease of exposition, we will drop the
time index $t$ throughout the section and let $\Theta_{t}=\theta$.

Recall that an $\epsilon$-LDP randomized mechanism is associated with a
conditional probability distribution that satisfies (1). An $\epsilon$-LDP
mechanism is not unique. One such mechanism, which we call Standard Randomized
Response (SRR), returns $Y=X$ with probability
$e^{\epsilon}/(e^{\epsilon}+K-1)$, otherwise returns a uniformly random number
on $\\{1,\ldots,K\\}/\\{X\\}$. The transition law corresponding to the SRR
from $X$ to $Y$ is

$\mathbb{P}_{\text{SRR}}(Y=y|X=x)=\frac{e^{\epsilon}}{e^{\epsilon}+K-1}\mathbb{I}(x=y)+\frac{1}{e^{\epsilon}+K-1}\mathbb{I}(x\neq
y).$

For subsequent use, we define $\texttt{SRR}(\epsilon,\mathcal{X})$ to be the
SRR mechanism that operates on a sample space $\mathcal{X}$ with DP parameter
$\epsilon$, and $\texttt{SRR}(x|K,\epsilon)$ to be its random output when it
takes $X=x$ as its input.

We aim to develop an alternative randomized mechanism such that the generated
randomized response $Y$ is more informative about $\theta^{\ast}$ than the one
generated by SRR. The main idea is as follows. Supposing that the posterior
sample $\Theta_{t}=\theta$ is an accurate estimate of $\theta^{\ast}$, it is
reasonable to aim for the ‘best’ $\epsilon$-LDP mechanism (among a range of
candidates) which would maximize the (estimation) utility of $Y$ if the true
parameter were $\theta^{\ast}=\theta$. We follow this main idea to develop the
proposed $\epsilon$-LDP mechanism.

###  4.1 The randomly restricted randomized response (RRRR) mechanism

Given $\Theta_{t}=\theta\in\Delta$, an informative randomized response
mechanism can be constructed by considering a _high-probability set_
$S\subset\\{1,\ldots,K\\}$ for $X$, and a _low-probability set_
$S^{c}=\\{1,\ldots,K\\}/S$ for $X$. Let
$\sigma_{\theta}=\sigma_{\theta}(1),\ldots,\sigma_{\theta}(K)$ be a
permutation vector for $\theta$ so that
$\theta_{\sigma_{\theta}(1)}\geq\ldots\geq\theta_{\sigma_{\theta}(K)}$. Then,
a high-probability set $S$ can be constructed as

$S=\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k)\\}$

for some $0\leq k<K$. Since $S$ contains the indices of the $k$ highest-valued
components of $\theta^{\ast}$, it is expected to cover a large portion of the
total probability for $X$. This can be the case even for a small value of $k$
relative to $K$ when the components of $\theta^{\ast}$ are not evenly
distributed. Hence, a sensible alternative to SRR would be to confine the
randomized response to the set $S$ (unioned by a random element from $S^{c}$
to remain LDP). The expected benefit of this approach is due to (i) using less
amount of randomization since $k<K$, and (ii) having an informative response
since $X$ belongs $S$ with a high probability.

Based on this approach, we propose the Randomly Restricted Randomized Response
(RRRR) mechanism, given in Algorithm 2. RRRR has three algorithmic parameters:
a subset $S$ of $\\{1,\ldots,K\\}$ and two privacy parameters $\epsilon_{1}$
and $\epsilon_{2}$ which operates on $S$ and $S^{c}$, respectively. Note that
when $|S|=0$ and $\epsilon_{2}=\epsilon$, RRRR reduces to SRR. Theorem 1
states the necessary conditions for $\epsilon_{1},\epsilon$ for RRRR to be
$\epsilon$-LDP. A proof of Theorem 1 is given in Appendix A.1.

Input: Sample space size $K$, a subset $S\subset\\{1,\ldots,K\\}$,
randomization parameters $\epsilon_{1},\epsilon_{2}>0$, input
$X\in\\{1,\ldots,K\\}$

Output: Randomized response $Y\in\\{1,\ldots,K\\}$

1

2Step 1: Generate $U\sim\textup{Uniform}(S^{c})$ and set
$R=U\cdot\mathbb{I}(X\in
S)+\texttt{SRR}(X|\epsilon_{2},S^{c})\cdot\mathbb{I}(X\notin S).$

3Step 2: Set $Z=X\cdot\mathbb{I}(X\in S)+R\cdot\mathbb{I}(X\in S^{c})$.

4Step 3: Set $Y=\texttt{SRR}(Z|\epsilon_{1},S\cup\\{R\\})$.

return _$Y$_

Algorithm 2 RRRR: Randomly restricted randomized response mechanism

######  Theorem 1.

RRRR is $\epsilon$-DP if $\epsilon_{1}\leq\epsilon$ and

$\epsilon_{2}=\begin{cases}\min\left\\{\epsilon,\ln\frac{|S^{c}|-1}{e^{\epsilon_{1}-\epsilon}|S^{c}|-1}\right\\}&\text{for
}\epsilon-\epsilon_{1}<\ln|S^{c}|\text{ and }|S|>0\\\
\epsilon&\text{else}\end{cases}.$ (2)

We elaborate on the behavior of RRRR and the choice of $\epsilon_{1}$ and
$\epsilon_{2}$ in the light of Theorem 1. The theorem implies that there is a
degree of freedom to choose $\epsilon_{1}<\epsilon$, while $\epsilon_{2}$ in
(2) is the largest possible choice given $\epsilon$ and $\epsilon_{1}$. In
RRRR, the probability of an honest response, i.e., $X=Y$, given $X\in S$ is

$\mathbb{P}_{\text{RRRR}}(Y=X|X\in
S)=\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+|S|},$

which should be contrasted to $\mathbb{P}_{\text{SRR}}(X=Y|X\in
S)=\frac{e^{\epsilon}}{e^{\epsilon}+K-1}$ for SRR. Anticipating that $\\{X\in
S\\}$ is likely, one should at least aim for a higher probability for $Y=X$
than that of SRR when $X\in S$ for RRRR to be relevant. That is, we would like
to choose $\epsilon_{1}$ to satisfy $\mathbb{P}_{\text{RRRR}}(X=Y|X\in
S)\geq\mathbb{P}_{\text{SRR}}(X=Y|X\in S)$, which is equivalent to

$\epsilon_{1}\geq\epsilon+\ln|S|-\ln(K-1).$ (3)

Taking into account also the constraint that $\epsilon_{1}\leq\epsilon$ (by
Theorem 1), we consider $\epsilon_{1}=\kappa\epsilon$, where $\kappa\in[0,1]$.
It is possible to choose $\kappa=1$; however, again by Theorem 1, this
requires that $\epsilon_{2}=0$, which renders $Y$ completely uninformative
when $X\notin S$. We suggest taking $\kappa$ a number close to $1$, such as
$0.9$, to ensure (3) with a significant margin. In Section 7, we discuss the
choice of $\kappa$ in more detail.

######  Example 1 (Numerical illustration).

We illustrate the benefit of RRRR and the role of $k$ with an example. We
considered $\theta$ values such that $\theta_{i}/\theta_{i+1}$ is constant for
$i=1,\ldots,K-1$. The ratio $\theta_{i}/\theta_{i+1}$ controls the degree of
‘evenness’ in $\theta$: The smaller ratio indicates a more evenly distributed
$\theta$. Note that $\theta$ is already ordered in this example; hence, we
consider using $S=\\{1,\ldots,k\\}$ which has the $k$ most likely values for
$X$ according to $\theta$. Also, for a given $\epsilon$, we fix
$\epsilon_{1}=0.9\epsilon$ and set $\epsilon_{2}$ according to (2).

Figure 2 shows, for a fixed $\epsilon$ and $K=20$, and various values of $k$,
the probability of the randomized response being equal to the sensitive
information, i.e., $\mathbb{P}_{\textup{RRRR}}(Y=X|\theta)$, vs
$\theta_{i}/\theta_{i+1}$ when $S=\\{1,\ldots,k\\}$ in RRRR. (This probability
is considered one of the utility functions presented in Sections 4.2.1-4.2.2.)
Comparing this probability with $\mathbb{P}_{\textup{SRR}}(Y=X|\theta)$, it
can be observed that RRRR can do significantly better than SRR if $k$ can be
chosen suitably. The plots demonstrate that the “suitable” $k$ depends on
$\theta$: While the best $k$ tends to be larger for more even $\theta$, small
$k$ becomes the better choice for non-even $\theta$ (large
$\theta_{i}/\theta_{i+1}$). This is because, when $\theta_{i}/\theta_{i+1}$ is
large, the probability is concentrated on just a few components, and $S$ with
a small $k$ captures most of the probability. Moreover, the plots for
$\epsilon=1$ and $\epsilon=5$ also show the effect of the level of privacy. In
more challenging scenarios where $\epsilon$ is smaller, the gain obtained by
RRRR compared to SRR is bigger.

Figure 2: $\mathbb{P}_{\textup{RRRR}}(Y=X|\theta)$ vs
$\theta_{i}/\theta_{i+1}$ for all $i=1,\ldots,K-1$ with $K=20$. Left:
$\epsilon=1$, Right: $\epsilon=5$.

###  4.2 Subset selection for RRRR

RRRR in Algorithm 2 can be used for any $S$ by choosing
$\epsilon_{1},\epsilon_{2}$ according to $|S|$, as suggested by Theorem 1. In
this paper, we assume a fixed $\kappa\in(0,1)$ and set
$\epsilon_{1}=\kappa\epsilon$; and we focus on the selection of $S$. (Note
that one can also seek to optimize $\epsilon_{1}$, too. However, this comes
with the expense of additional computation due to the requirement to explore a
larger hyperparameter space.)

To formalize the discussion, let $\texttt{RRRR}(S,\epsilon)$ be the RRRR
mechanism that achieves $\epsilon$-LDP by using the subset $S$
$\epsilon_{1}=\kappa\epsilon$ and let $\texttt{RRRR}(x|S,\epsilon)$ be its
random output when the input is $X=x$. Furthermore, let $U(\theta,S,\epsilon)$
be the ‘utility’ (concerning inference) when the true parameter is
$\theta^{\ast}=\theta$ and $Y=\texttt{RRRR}(X|S,\epsilon)$. One would like to
choose $S$ such that randomized response $Y$ maximizes $U(\theta,S,\epsilon)$.
However, since there are $2^{K}-1$ feasible choices for $S$, one must confine
the search space for $S$ in practice. As suggested in Example 1, RRRR becomes
most relevant when the set $S$ is chosen to have the indices of the $k$
_largest_ components of $\theta$. Therefore, for a given $\theta$, we confine
the choices for $S$ to

$S_{\theta,k}:=\\{\sigma_{\theta}(1),\sigma_{\theta}(2),\ldots,\sigma_{\theta}(k)\\},\quad
k=1,\ldots,K.$ (4)

Then the subset selection problem can be formulated as finding

$k^{\ast}=\arg\max_{k\in\\{0,\ldots,K-1\\}}U(\theta,S_{k,\theta},\epsilon).$
(5)

Note that the set of solutions considered in (5) covers SRR, obtained with
$k=0$ (leading to $S=\emptyset$ and $\epsilon_{2}=\epsilon$).

In the subsequent sections, we present six different utility functions and
justify their relevance to estimation; the usefulness of the proposed
functions is also demonstrated in the numerical experiments.

####  4.2.1 Fisher information matrix

The first utility function under consideration is based on the Fisher
information matrix at $\theta$ according to the distribution of $Y$ given
$\theta$. It is well-known that the inverse of the Fisher information matrix
sets the Cramer-Rao lower bound for the variance of an unbiased estimator
(Casella and Berger, , 2021). Hence the Fisher information can be regarded as
a reasonable metric to quantify the information contained in $Y$ about
$\theta$. This approach is adopted in Lopuhaä-Zwakenberg et al.,  (2022);
Steinberger,  (2024) for LDP applications for estimating discrete
distributions, and Alparslan and Yıldırım,  (2022); Yıldırım,  (2024) in
similar problems involving parametric continuous distributions.

For a given $\theta\in\Delta$, let $F(\theta;S,\epsilon)$ be the Fisher
information matrix evaluated at $\theta$ when $X\sim\text{Cat}(\theta)$ and
$Y=\texttt{RRRR}(X|S,\epsilon)$. Let $g_{S,\epsilon}$ be the conditional
probability mass function of the distribution of $Y$ given $X$ according to
$Y=\texttt{RRRR}(X|S,\epsilon)$, that is,

$g_{S,\epsilon}(y|x):=\mathbb{P}_{\textup{RRRR}}(Y=y|X=x)$

when $Y=\texttt{RRRR}(X|S,\epsilon)$. The following result states
$F(\theta;S,\epsilon)$ in terms of $g_{S,\epsilon}$ and $\theta$. The result
is derived in Lopuhaä-Zwakenberg et al.,  (2022); we also give a simple proof
in Appendix A.2. Note that $F(\theta;S,\epsilon)$ is $(K-1)\times(K-1)$ since
$\theta$ has $K-1$ free components and
$\theta_{K}=1-\sum_{i=1}^{K-1}\theta_{i}$.

######  Proposition 1.

The Fisher information matrix for RRRR is given by

$F(\theta;S,\epsilon)=A_{S,\epsilon}^{\top}D_{\theta}^{-1}A_{S,\epsilon},$ (6)

where $A_{S,\epsilon}$ is a $K\times(K-1)$ matrix whose entries are
$A_{S,\epsilon}(i,j):=g_{S,\epsilon}(i|j)-g_{S,\epsilon}(i|K)$ and
$D_{\theta}$ is a $K\times K$ diagonal matrix with elements
$D_{\theta}(i,i):=\sum_{j=1}^{K}g_{S,\epsilon}(i|j)\theta_{j}$.

We define the following utility function based on the Fisher information

$U_{1}(\theta,S,\epsilon):=-\textup{Tr}\left[F^{-1}(\theta;S,\epsilon)\right].$
(7)

This utility function depends on the Fisher information differently from
Lopuhaä-Zwakenberg et al.,  (2022); Steinberger,  (2024), who considered the
determinant of the FIM as the utility function. The rationale behind (7) is
that the for an unbiased estimator $\hat{\theta}(Y)$ of $\theta^{\ast}$ based
on $Y=\texttt{RRRR}(X|S,\epsilon)$ and $X\sim\text{Cat}(\theta^{\ast})$, the
expected mean squared error is bounded by
$E_{\theta^{\ast}}[\|\hat{\theta}(Y)-\theta^{\ast}\|^{2}]\leq\textup{Tr}\left[F^{-1}(\theta^{\ast};S,\epsilon)\right]$.
For the utility function in (7) to be well defined, the FIM needs to be
invertible. Proposition 2, proven in Appendix A.2, states that this is indeed
the case.

######  Proposition 2.

$F(\theta;S,\epsilon)$ in (6) is invertible for all $\theta\in\Delta$,
$S\subset\\{1,\ldots,K\\}$, and $\epsilon_{1},\epsilon_{2}>0$.

####  4.2.2 Entropy of randomized response

For discrete distributions, entropy measures the uniformity. Hence in the LDP
framework, a lower entropy for the randomized response $Y$ implies a more
informative $Y$. Based on that observation, a utility function can be defined
as the negative entropy of the marginal distribution of $Y$.

$U_{2}(\theta,S,\epsilon):=\sum_{y=1}^{K}\ln
h_{S,\epsilon}(y|\theta)h_{S,\epsilon}(y|\theta),$

where $h_{S,\epsilon}(y|\theta)$ is the marginal probability of $Y=y$ given
$\theta$,

$h_{S,\epsilon}(y|\theta):=\sum_{x=1}^{K}g_{S,\epsilon}(y|x)\theta_{x}.$

####  4.2.3 Total variation distance

The TV distance between two discrete probability distributions $\mu,\nu$ on
$\mathcal{X}$ is given by

$\textup{TV}(\mu,\nu):=\frac{1}{2}\sum_{x\in\mathcal{X}}|\mu(x)-\nu(x)|.$

We consider two utility functions based on TV distance. The first function
arises from the observation that a more informative response $Y$ generally
leads to a larger change in the posterior distribution of $X$ given
$Y,\theta$,

$p_{S,\epsilon}(x|y,\theta):=\frac{\theta_{x}\cdot
g_{S,\epsilon}(y|x)}{h_{S,\epsilon}(y|\theta)},\quad x=1,\ldots,K,$ (8)

relative to its prior $\text{Cat}(\theta)$. The expected amount of change can
be formulated as the expectation of the TV distance between the prior and
posterior distributions with respect to the marginal distribution of $Y$ given
$\theta$. Then a utility function can be defined as

$\displaystyle U_{3}(\theta,S,\epsilon)$
$\displaystyle:=\mathbb{E}_{\theta}\left[\text{TV}(p_{S,\epsilon}(\cdot|Y,\theta),\text{Cat}(\theta))\right]$
$\displaystyle=\frac{1}{2}\sum_{x=1}^{K}\sum_{y=1}^{K}\left|g_{S,\epsilon}(y|x)\theta_{x}-h_{S,\epsilon}(y|\theta)\theta_{x}\right|.$

Another utility function is related to the TV distance between the marginal
probability distributions of $X$ given $\theta$ and $Y$ given $\theta$. Since
$X$ is more informative about $\theta$ than the randomized response $Y$, the
mentioned TV distance is desired to be as small as possible. Hence, a utility
function may be formulated as

$\displaystyle U_{4}(\theta,S,\epsilon)$
$\displaystyle:=-\text{TV}(h_{S,\epsilon}(\cdot|\theta),\text{Cat}(\theta))$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{K}|h_{S,\epsilon}(i|\theta)-\theta_{i}|.$

####  4.2.4 Expected mean squared error

One can also wish to choose $S$ such that the Bayesian estimator of $X$ given
$Y$ has the lowest expected squared error. Specifically, given
$k\in\\{1,\ldots,K\\}$ let $e_{k}$ be a $K\times 1$ vector of $0$s except that
$e_{k}=1$. A utility function can be defined based on that as

$U_{5}(\theta,S,\epsilon):=-\arg\min_{\widehat{e_{X}}}\mathbb{E}_{\theta}\left[\|e_{X}-\widehat{e_{X}}(Y)\|^{2}\right],$
(9)

where $\mathbb{E}_{\theta}\left[\|e_{X}-\widehat{e_{X}}(Y)\|^{2}\right]$ is
the mean squared error for the estimator $\widehat{e_{X}}$ of $e_{X}$ given
$Y$ when $X\sim\text{Cat}(\theta)$ and $Y=\texttt{RRRR}(X|S,\epsilon)$, which
is known to be minimized when $\widehat{e_{X}}$ is the Bayesian estimator of
$e_{X}$. Proposition 3 provides an explicit formula for this utility function.
A proof is given in Appendix A.2.

######  Proposition 3.

For the utility function in (9), we have

$U_{5}(\theta,S,\epsilon)=\sum_{y=1}^{K}\sum_{x=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)}-1.$

####  4.2.5 Probability of honest response

Our last alternative for the utility function is a simple yet intuitive one,
which is the probability of an honest response, i.e.,

$U_{6}(\theta,S,\epsilon):=\mathbb{P}(Y=X|S,\theta).$ (10)

This probability is explicitly given by

$\displaystyle\mathbb{P}_{\theta}(Y=X|S)$ $\displaystyle=\mathbb{P}(Y=X|X\in
S)\mathbb{P}_{\theta}(X\in S)+\mathbb{P}(Y=X|X\notin
S)\mathbb{P}_{\theta}(X\notin S)$
$\displaystyle=\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+|S|}\left(\sum_{i\in
S}\theta_{i}+\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-|S|-1}\sum_{i\notin
S}\theta_{i}\right).$

Recall that, for computational tractability, we confined the possible sets for
$S$ to the subsets $\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k)\\}$,
$k=0,\ldots,K-1$ and select $S$ by solving the maximization problem in (5).
Remarkably, if $U_{6}(\theta,S,\epsilon)$ is used for the utility function,
the restricted maximization (5) is equivalent to _global_ maximization, i.e.,
finding the best $S$ among all the $2^{K}$ possible subsets $S$. We state this
as a theorem and prove it in Appendix A.2.

######  Theorem 2.

For the utility function $U_{6}(\theta,S,\epsilon)$ in (10) and
$S_{k,\theta}$s in (4), we have

$\max_{k=0,\ldots,K-1}U_{6}(\theta,S_{k,\theta},\epsilon)=\max_{S\subset\\{1,\ldots,K\\}}U_{6}(\theta,S,\epsilon).$

####  4.2.6 Semi-adaptive approach

We also consider a semi-adaptive approach which uses a fixed parameter
$\alpha\in(0,1)$ to select the smallest $S_{k,\theta}$ in (4) such that
$\mathbb{P}_{\theta}(X\in S_{k,\theta})\geq\alpha$, that is,
$S=\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k^{\ast})\\}$ is taken such
that

$\mathbb{P}_{\theta}(X\in\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k^{\ast}-1)\\})<\alpha\text{
and
}\mathbb{P}_{\theta}(X\in\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k^{\ast})\\})\geq\alpha.$

Again, the idea is to randomize the most likely values of $X$ with a high
accuracy. The approach forms the subset $S$ by including values for $X$ in
descending order of their probabilities (given by $\theta$) until the
cumulative probability exceeds $\alpha$. In that way, it is expected to have
set $S$ that is small-sized (especially when $\theta$ is unbalanced) and
captures the most likely values of $X$. The resulting $S$ has varying
cardinality depending on the sampled $\theta$ at the current time step.

We call this approach “semi-adaptive” because, while it still adapts to
$\theta$, it uses the fixed parameter $\alpha$ which is not updated during
AdOBEst-LDP. As we will see in Section 7, the best $\alpha$ depends on various
parameters such as $\epsilon$, $K$, and the degree of evenness in $\theta$.

###  4.3 Computational complexity of utility functions

We now provide the computational complexity analysis of the utility metrics
presented in Section 4.2.1-4.2.5, and that of the semi-adaptive approach in
Section 4.2.6, as a function of $K$. The first row of Table 1 shows the
computational complexities of calculating the utility function for a fixed
$S$, and the second row shows the complexities of choosing the best $S$
according to (5). To find (5), the utility function generally needs to be
calculated $K$ times, which explains the additional $K$ factor in the
computational complexities in the second row.

The least demanding utility function is $U_{6}$, that is based on
$\mathbb{P}_{\theta}(Y=X)$, whose complexity is $\mathcal{O}(K)$. Moreover,
finding the best $S$ can also be done in $\mathcal{O}(K)$ time because one can
compute this utility metric for all $k=0,\ldots,K-1$ by starting with
$S=\emptyset$ and populating it incrementally. Also note that the semi-
adaptive approach does not use a utility metric and finding $k^{\ast}$ can be
done in $\mathcal{O}(K)$ time by summing the components of $\theta$ from
largest to smallest until the cumulative sum exceeds the given $\alpha$
parameter. So, its complexity is $\mathcal{O}(K)$.

For all these approaches, it is additionally required to sort $\theta$
beforehand, which is an $\mathcal{O}(K\ln K)$ operation with an efficient
sorting algorithm like merge sort.

| Fisher | Entropy | $\text{TV}_{1}$ | $\text{TV}_{2}$ | MSE | $\mathbb{P}_{\theta}(Y=X)$ | Semi-adaptive
---|---|---|---|---|---|---|---
Computing utility | $\mathcal{O}(K^{3})$ | $\mathcal{O}(K^{2})$ | $\mathcal{O}(K^{2})$ | $\mathcal{O}(K^{2})$ | $\mathcal{O}(K^{2})$ | $\mathcal{O}(K)$ | NA
Choosing $S$ | $\mathcal{O}(K^{4})$ | $\mathcal{O}(K^{3})$ | $\mathcal{O}(K^{3})$ | $\mathcal{O}(K^{3})$ | $\mathcal{O}(K^{3})$ | $\mathcal{O}(K)$ | $\mathcal{O}(K)$
Table 1: Computational complexity of utility functions and choosing $S$

##  5 Posterior sampling

Step 3 of Algorithm 1 requires sampling from the posterior distribution
$\Pi(\cdot|Y_{1:n},S_{1:n})$ of $\theta$ given $Y_{1:n}$ and $S_{1:n}$, where
$S_{t}$ is the subset selected at time $t$ to generate $Y_{t}$ from $X_{t}$.
Let $\pi(\theta|Y_{1:n},S_{1:n})$ denote the pdf of
$\Pi(\cdot|Y_{1:n},S_{1:n})$. Given $Y_{1:n}=y_{1:n}$ and $S_{1:n}=s_{1:n}$,
the posterior density can be written as

$\displaystyle\pi(\theta|y_{1:n},s_{1:n})$
$\displaystyle\propto\eta(\theta)\prod_{t=1}^{n}h_{s_{t},\epsilon}(y_{t}|\theta).$
(11)

Note that the right-hand side does not include a transition probability for
$S_{t}$’s because the sampling procedure of $S_{t}$ given $Y_{1:t-1}$ and
$S_{1:t-1}$ does _not_ depend on $\theta$. Furthermore, we assume that the
prior distribution $\eta(\theta)$ is a Dirichlet distribution
$\theta\sim\textup{Dir}(\rho_{1},\ldots,\rho_{K})$ with prior hyperparameters
$\rho_{k}>0$, for $k=1,\ldots,K$.

Unfortunately, exact sampling from the posterior in (11) is difficult.
Therefore, we resort to approximate sampling approaches using MCMC. Below we
present two MCMC sampling methods. The first method is based on SGLD, which
works with random small-sized subsamples and therefore is computationally
scalable in $n$. SGLD is an asymptotically exact MCMC sampling approach,
targeting the posterior distribution approximately. The second method is a
Gibbs sampler, an exact MCMC method unlike the SGLD. However, the
computational cost of the Gibbs sampler for $n$ time steps is
$\mathcal{O}(n^{2}K)$, making the Gibbs sampler affordable when $n$ is small
to moderate, but too expensive for large $n$. In our numerical experiments, we
preferred SGLD due to its computational advantage over the Gibbs sampler.

###  5.1 Stochastic gradient Langevin dynamics

SGLD is an asymptotically exact MCMC sampling approach that enables the use of
subsamples of size $m\ll t$. A direct application of SGLD to generate samples
for $\theta$ from the posterior distribution in (11) is difficult. This is
because $\theta$ lives in the probability simplex $\Delta$, which makes the
task of keeping the iterates for $\theta$ inside $\Delta$ challenging. We
overcome this problem by defining the surrogate variables
$\phi_{1},\ldots,\phi_{K}$ with

$\phi_{k}\overset{\text{ind.}}{\sim}\text{Gamma}(\rho_{k},1),\quad
k=1,\ldots,K,$

and the mapping from $\phi$ to $\theta$ as

$\theta(\phi)_{k}:=\frac{\phi_{k}}{\sum_{j=1}^{K}\phi_{j}},\quad
k=1,\ldots,K.$ (12)

It is well-known that the resulting $(\theta_{1},\ldots,\theta_{K})$ has a
Dirichlet distribution $\text{Dir}(\rho_{1},\ldots,\rho_{K})$, which is
exactly the prior distribution $\eta(\theta)$. Therefore, this change of
variables preserves the originally constructed probabilistic model. Moreover,
since $\phi=(\phi_{1},\ldots,\phi_{K})$ takes values in $[0,\infty)^{K}$, we
run SGLD for $\phi$, where the $j$’th update is

$\phi^{(j)}=\left|\phi^{(j-1)}+\frac{\gamma}{2}\left(\nabla_{\phi}\ln
p(\phi^{(j-1)})+\frac{t}{m}\sum_{i=1}^{m}\nabla_{\phi}\ln
p(y_{u_{i}}|\phi^{(j-1)})\right)+\gamma W_{t}\right|,\quad
W_{t}\sim\mathcal{N}(0,I_{K}).$ (13)

where $u=(u_{1},\ldots,u_{m})$ is a random subsample of $\\{1,\ldots,t\\}$. In
(13), the ‘new’ prior and likelihood functions are

$p(\phi)=\prod_{k=1}^{K}\text{Gamma}(\phi_{k};\alpha_{i},1),\quad
p(y|\phi)=h_{S,\epsilon}(y|\theta(\phi)).$ (14)

The reflection in (13) via taking the component-wise absolute value is
necessary because each $\phi_{i}^{(j)}$ must be positive. Step 3 of Algorithm
1 can be approximated by running SGLD for some $M>0$ iterations. To exploit
the SGLD updates from the previous time, one should start the updates at time
$t$ by setting the initial value for $\phi$ to the last SGLD iterate at time
$t-1$.

The next proposition provides the explicit formulae for the gradients of the
log-prior and the log-likelihood of $\phi$ in (13). A proof is given in
Appendix A.3.

######  Proposition 4.

For $p(\phi)$ and in $p(y|\phi)$ in (14), we have

$[\nabla_{\phi}\ln
p(\phi)]_{i}=\frac{\alpha_{i}-1}{\phi_{i}}-1,\quad[\nabla_{\phi}\ln
p(y|\phi)]_{i}=\sum_{k=1}^{K-1}J(i,k)\frac{g_{S,\epsilon}(y|k)-g_{S,\epsilon}(y|K)}{h_{S,\epsilon}(y|\theta(\phi))},$

where $J$ is a $K\times(K-1)$ Jacobian matrix whose $(i,j)$th element is

$J(i,j)=\mathbb{I}(i=j)\frac{1}{\sum_{k=1}^{K}\phi_{k}}-\frac{\phi_{j}}{\left(\sum_{k=1}^{K}\phi_{k}\right)^{2}}.$

###  5.2 Gibbs sampling

An alternative to SGLD is the Gibbs sampler, which operates on the joint
posterior distribution of $\theta$ and $X_{1:n}$ given $Y_{1:n}=y_{1:n}$ and
$S_{1:t}=s_{1:n}$,

$p(\theta,x_{1:n}|y_{1:n},s_{1:n})\propto\eta(\theta)\left[\prod_{t=1}^{n}\theta_{x_{t}}g_{s_{t},\epsilon}(y_{t}|x_{t})\right].$

The full conditional distributions of $X_{1:n}$ and $\theta$ are tractable.
Specifically, for $X_{1:n}$, we have

$p(x_{1:n}|y_{1:n},s_{1:n},\theta)=\prod_{t=1}^{n}p_{s_{t},\epsilon}(x_{t}|y_{t},\theta),$
(15)

where $p_{s_{t},\epsilon}(x_{t}|y_{t},\theta)$ is defined in (8). Therefore,
(15) is a product of $n$ categorical distributions each with support
$\\{1,\ldots,K\\}$. Furthermore, the full conditional distribution of $\theta$
is a Dirichlet distribution due to the conjugacy between the categorical and
the Dirichlet distributions. Specifically,

$p(\theta|x_{1:n},y_{1:n},s_{1:n})=\text{Dir}(\theta|\rho^{\text{post}}_{1},\ldots,\rho^{\text{post}}_{K}),$

where the hyperparameters of the posterior distribution are given by
$\rho^{\text{post}}_{k}:=\rho_{k}+\sum_{t=1}^{n}\mathbb{I}(x_{t}=k)$ for
$k=1,\ldots,K$.

Computational load at time $t$ of sampling from $t$ distributions in (15), is
proportional to $tK$, which renders the computational complexity of Gibbs
sampling $\mathcal{O}(n^{2}K)$ after $n$ time steps. This can be
computationally prohibitive when $n$ gets large.

##  6 Theoretical analysis

We address two questions concerning AdOBEst-LDP in Algorithm 1 when run with
RRRR. (i) Does the targeted posterior distribution based on the observations
generated by Algorithm 1 converge to the true value $\theta^{\ast}$? (ii) How
frequently does Algorithm 1 with RRRR select the optimum subset $S$ according
to the chosen utility function?

###  6.1 Convergence of the posterior distribution

We begin by developing the joint probability distribution of the random
variables involved in AdOBEst-LDP.

  * •

Given $Y_{1:n}$ and $S_{1:n}$, the posterior distribution
$\Pi(\cdot|Y_{1:n},S_{1:n})$ is defined such that for any measurable set
$A\subseteq\Delta$, the posterior probability of $\\{\theta\in A\\}$ is given
by

$\Pi(A|Y_{1:n},S_{1:n}):=\frac{\int_{A}\eta(\theta)\prod_{t=1}^{n}h_{S_{t},\epsilon}(y_{t}|\theta)\mathrm{d}\theta}{\int_{\Delta}\eta(\theta)\prod_{t=1}^{n}h_{S_{t},\epsilon}(y_{t}|\theta)\mathrm{d}\theta}.$
(16)

  * •

Let $Q(\cdot|Y_{1:n},S_{1:n},\Theta_{n-1})$ be the probability distribution
corresponding to the posterior sampling process for $\Theta_{n}$. Note that,
if perfect posterior sampling were used, we would have
$Q(A|Y_{1:n},S_{1:n},\Theta_{n-1})=\Pi(A|Y_{1:n},S_{1:n})$; however, when
approximate sampling techniques are used to target $\Pi$, such as SGLD or
Gibbs sampling, the equality does not hold in general.

  * •

For $\theta\in\Delta$, let

$S^{\ast}_{\theta}:=\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k_{\theta}^{\ast})\\},\quad
k^{\ast}_{\theta}:=\arg\max_{k\in\\{0,\ldots,K-1\\}}U(\theta,S_{k,\theta},\epsilon)$

be the best subset according to $\theta$, where
$S_{k,\theta}=\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k)\\}$ is defined
in (4). Given $\Theta_{1:t-1}$ and $Y_{1:t}$, $S_{t}$ depends only on
$\Theta_{t-1}$ and it is given by $S_{t}=S^{\ast}_{\Theta_{t-1}}$.

Combining all, the joint law of $S_{1:n},Y_{1:n}$ can be expressed as

$P_{\theta^{\ast}}(S_{1:n},Y_{1:n}):=\prod_{t=1}^{n}h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\left[\int_{\Delta}\mathbb{I}(S_{t}=S_{k^{\ast},\theta_{t-1}})Q(\mathrm{d}\theta_{t-1}|Y_{1:t-1},S_{1:t-1},\theta_{t-2})\right],$
(17)

where we use the convention that
$Q(\mathrm{d}\theta_{0}|Y_{1:0},S_{1:0},\theta_{-1})=\delta_{\theta_{\text{init}}}(\mathrm{d}\theta_{0})$
for an initial value $\theta_{\text{init}}$.

The posterior probability in (16) is a random variable with respect to
$P_{\theta^{\ast}}$ defined in (17). Theorem 3 establishes that under the
fairly mild Assumption 1 on the prior, the $\Pi(\cdot|Y_{1:n},S_{1:n})$
converges to $\theta^{\ast}$ regardless of the choice of $Q$ for posterior
sampling.

######  Assumption 1.

There exist finite positive constants $d>0$ and $B>0$ such that
$\eta(\theta)/\eta(\theta^{\prime})<B$ for all
$\theta,\theta^{\prime}\in\Delta$ whenever
$\|\theta^{\prime}-\theta^{\ast}\|<d$.

######  Theorem 3.

Under Assumption 1, there exists a constant $c>0$ such that, for any $0<a<1$
and the sequence of sets

$\Omega_{n}=\\{\theta\in\Delta:\|\theta-\theta^{\ast}\|^{2}\leq cn^{-a}\\},$

the sequence of probabilities

$\lim_{n\rightarrow\infty}\Pi(\Omega_{n}|Y_{1:n},S_{1:n})\overset{P_{\theta^{\ast}}}{\rightarrow}1,$

regardless of the choice of $Q$.

A proof is given in Appendix A.4.2, where the constant $c$ in the sets
$\Omega_{n}$ is explicitly given.

###  6.2 Selecting the best subset

Let $S^{\ast}:=S^{\ast}_{\theta^{\ast}}$ be the best subset at
$\theta^{\ast}$. In this part, we prove that, if posterior sampling is
performed exactly, the best subset is chosen with an expected long-run
frequency of $1$. Our result relies on some mild assumptions.

######  Assumption 2.

The components of $\theta^{\ast}$ are strictly ordered, that is,
$\theta_{\sigma_{\theta^{\ast}}(1)}>\ldots>\theta_{\sigma_{\theta^{\ast}}(K)}$.

######  Assumption 3.

Given any $S\subset\\{1,\ldots,K\\}$ and $\epsilon>0$, $U(\theta,S,\epsilon)$
is a continuous function of $\theta$ with respect to the $L_{2}$-norm.

######  Assumption 4.

The solution of (5) is unique at $\theta^{\ast}$.

Assumption 2 is required to avoid technical issues regarding the uniqueness of
$S^{\ast}$. Assumptions 3 and 4 impose a certain form of regularity on the
utility function.

######  Theorem 4.

Suppose Assumptions 2-4 hold and $\Theta_{t}$s are generated by perfect
sampling, that is, $Q(A|Y_{1:t},S_{1:t})=\Pi(A|Y_{1:t},S_{1:t})$ for all
measurable $A\subseteq\Delta$. Then,

$P_{\theta^{\ast}}(S_{t}=S^{\ast})\rightarrow 1.$ (18)

As a corollary, $S^{\ast}$ is selected with an expected long-run frequency of
$1$, that is,

$\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{t=1}^{n}E_{\theta^{\ast}}\left[\mathbb{I}(S_{t}=S^{\ast})\right]=1.$
(19)

The result in (19) can be likened to sublinear regret from the reinforcement
learning theory.

##  7 Numerical results

We tested the performance of AdOBEst-LDP run with RRRR for randomization when
the subset $S$ in RRRR is determined according to a utility function in
Section 4.2.111The MATLAB code at https://github.com/soneraydin/AdOBEst_LDP
can be used to reproduce the results obtained in this paper. We compared
AdOBEst-LDP when combined with each of the utility functions defined in
Sections 4.2.1-4.2.5 with its non-adaptive counterpart when SRR is used to
generate $Y_{t}$ at all steps. We also included the semi-adaptive subset
selection method in Section 4.2.6 into the comparison. For the semi-adaptive
approach, we obtained results for 5 different values of its $\alpha$
parameter, namely $\alpha\in\\{0.2,0.6,0.8,0.9,0.95\\}$.

We ran each method for 50 Monte Carlo runs. Each run contained $T=500K$ time
steps. For each run, the sensitive information is generated as
$X_{t}\overset{\text{i.i.d.}}{\sim}\text{Cat}(\theta^{\ast})$ where
$\theta^{\ast}$ itself was randomly drawn from
$\text{Dirichlet}(\rho,\ldots,\rho)$. Here, the parameter $\rho$ was used to
control the unevenness among the components of $\theta^{\ast}$. (Smaller
$\rho$ leads to more uneven components in general). At each time step, Step 3
of Algorithm 1 was performed by running $M=20$ updates of an SGLD-based MCMC
kernel as described Section 5.1. In SGLD, we took the subsample size $m=50$
and the step-size parameter $a=\frac{0.5}{T}$. Prior hyperparameters for the
gamma distribution were taken $\rho_{0}=1_{K}$. The posterior sample
$\Theta_{t}$ was taken as the last iterate of those SGLD updates. Only for the
last time step, $t=T$, the number of MCMC iterations was taken $2000$ to
reliably calculate the final estimate $\hat{\theta}$ of $\theta$ by averaging
the last $1000$ of those $2000$ iterates. (This average is the MCMC
approximation of the posterior mean of $\theta$ given $Y_{1:T}$ and
$S_{1:T}$.) We compared the mean posterior estimate of $\theta$ and the true
value and the performance measure was taken as the TV distance between
$\text{Cat}(\theta^{\ast})$ and $\text{Cat}(\hat{\theta})$, that is,

$\frac{1}{2}\sum_{i=1}^{K}|\hat{\theta}_{i}-\theta_{i}|.$ (20)

Finally, the comparison among the methods was repeated for all the
combinations ($K\times\epsilon\times\epsilon_{1}\times\rho$) of
$K\in\\{10,20\\}$, $\epsilon\in\\{0.5,1,5\\}$,
$\epsilon_{1}\in\\{0.8\epsilon,0.9\epsilon\\}$, and $\rho\in\\{0.01,0.1,1\\}$.

The accuracy results for the methods in comparison are summarized in Figures 3
and 4 in terms of the error given in (20). The box plots are centered at the
error median and the whiskers stretch from the minimum to the maximum over the
50 MC runs excluding the outliers. When the medians are compared, the fully
adaptive algorithms, which use a utility function to select $S_{t}$, yield
comparable results to the best semi-adaptive approach in both figures. As one
may expect, the non-adaptive approach yielded the worst results in general,
especially in the high-privacy regimes (smaller $\epsilon$) and uneven
$\theta^{\ast}$ (smaller $\rho$). We also observe that, while most utility
metrics are generally robust, the one based on FIM seems sensitive to the
choice of $\epsilon_{1}$ parameter. This can be attributed to the fact that
the FIM approaches singularity when $\epsilon_{2}$ is too small, which is the
case if $\epsilon_{1}$ is chosen too close to $\epsilon$. Supporting this, we
see that when $\epsilon_{1}=0.8\epsilon$, the utility metric based on FIM
becomes more robust. Another remarkable observation is that the utility
function based on the probability of honest response, $U_{6}$, has competitive
performance despite being the lightest utility metric in computational
complexity. Finally, while the semi-adaptive approach is computationally less
demanding than most fully adaptive versions, the results show it can
dramatically fail if its $\alpha$ hyperparameter is not tuned properly. In
contrast, the fully adaptive approaches adapt well to $\epsilon$ or $\rho$ and
do not need additional tuning.

In addition to the error graphs, the heat maps in Figures 5 and 6 show the
effect of parameters $\rho$ and $\epsilon$ on the average cardinality of the
subsets $S$ chosen by each algorithm (again, averaged over 50 Monte Carlo
runs). According to these figures, increasing the value of $\rho$ causes an
increase in the cardinalities of subsets chosen by each algorithm (except the
nonadaptive one, since it uses all $K$ categories, rather than a smaller
subset). This is expected, since higher $\rho$ values cause
$\text{Cat}(\theta^{\ast})$ to be closer to the uniform distribution, thus
causing $X$ to be more evenly distributed among the categories. Moreover, for
small $\rho$, increasing the value of $\epsilon$ causes a decrease in the
cardinalities of these subsets, which can be attributed to a higher $\epsilon$
leading to a more accurate estimation. When we compare the utility functions
for the adaptive approach among themselves, we observe that for
$\epsilon_{1}=0.8\epsilon$, the third utility function (TV1) uses the subsets
with the largest cardinality (on average). However, when we increase the
$\epsilon_{1}$ value to $\epsilon_{1}=0.9\epsilon$, the second utility
function (FIM) uses the subsets with the largest cardinality. This might be
due to the sensitivity of the FIM-based utility function to the choice of
$\epsilon_{1}$ parameter that we mentioned before, which affects the
invertibility of the Fisher information matrix when $\epsilon_{1}$ is too
close to $\epsilon$.

Figure 3: TV distance in (20) for $K\in\\{10,20\\}$,
$\epsilon_{1}=0.8\epsilon$

Figure 4: TV distance in (20) for $K\in\\{10,20\\}$,
$\epsilon_{1}=0.9\epsilon$

Figure 5: Average cardinalities of the subsets selected by each method, for
$K\in\\{10,20\\}$, $\epsilon_{1}=0.8\epsilon$

Figure 6: Average cardinalities of the subsets selected by each method, for
$K\in\\{10,20\\}$, $\epsilon_{1}=0.9\epsilon$

##  8 Conclusion

In this paper, we proposed a novel adaptive framework, AdOBEst-LDP, for online
estimation of the frequency distribution of categorical data under the
$\epsilon$-LDP constraint. AdOBEst-LDP, run with RRRR for randomization,
encompasses both the privatization/perturbation of the given data and the
accurate estimation of population parameters from privatized data, in a
dynamic way, using Bayesian methodology. Our privatization mechanism (RRRR) is
distinguished from the baseline approach (SRR) in a way that it operates on a
smaller subset of the categories rather than using all categories which is
often suboptimal in terms of accuracy. Finding the optimal subset for this
task in a static way is a computationally inefficient approach whose
complexity is in the order of $2^{K}$. For that reason, we employed an
adaptive approach to dynamically update the subsets at each iteration, given
the current sample of the population parameter and the privacy parameter.
Optimal selection of these subsets was guided by various alternative utility
functions that we used throughout the paper, and in each of them, their
‘utility’ was based on optimizing a different information metric to select the
most informative subset among all possible combinations, in a computationally
tractable way.

For the posterior sampling of $\theta$ at each iteration, we employed an
efficient SGLD-based sampling scheme on a constrained region, namely the
$K$-dimensional probability simplex. We distinguished this scheme from Gibbs
sampling which uses all of the historical data, which is not scalable to large
datasets.

In the numerical experiments, we demonstrated that AdOBEst-LDP can estimate
the population distribution more accurately than the non-adaptive approach
under experimental settings with various privacy levels $\epsilon$ and degrees
of evenness among the components of $\theta^{\ast}$. While the performance of
AdOBEst-LDP is generally robust for all the utility functions considered in
the paper, the utility function based on the probability of honest response
can be preferred due to its much lower computational complexity than the other
utility functions. Our experiments also showed that the accuracy of the
adaptive approach is comparable to that of the semi-adaptive approach.
However, the semi-adaptive approach requires determining an optimal parameter
($\alpha$) in a continuous space, so the adaptive approach is practically more
useful in real-life applications.

Via theoretical analysis, we showed that, regardless of whether the posterior
sampling is conducted exactly or approximately, the target posterior
distribution of AdOBEst-LDP converges to the true population parameter
$\theta^{\ast}$. We also showed that, under exact posterior sampling, the best
subset given utility function is selected with probability $1$ in the long
run.

Note that the observations $\\{Y_{t}\\}_{t\geq 1}$ generated by AdOBEst-LDP
are dependent. Hence, the theoretical work in Section 6 can also be regarded
as a contribution to the literature on the convergence of posterior
distributions under dependent observations. Moreover, we have already pointed
out a form of analogy between AdOBEst-LDP and Thompson sampling (Russo et al.,
, 2018). Observe that posterior sampling is present in both methods and the
subset selection step can be likened to the action step of a reinforcement
learning scheme. In that sense, we also believe that the theoretical results
can inspire studies on the convergence of dynamic reinforcement learning
algorithms, particularly Thompson sampling.

## References

  * Acharya et al., (2023) Acharya, J., Canonne, C. L., Sun, Z., and Tyagi, H. (2023).  Unified lower bounds for interactive high-dimensional estimation under information constraints.  In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S., editors, Advances in Neural Information Processing Systems, volume 36, pages 51133–51165, New Orleans, US. Curran Associates, Inc.
  * Alparslan and Yıldırım, (2022) Alparslan, B. and Yıldırım, S. (2022).  Statistic selection and MCMC for differentially private Bayesian estimation.  Statistics and Computing, 32(5):66.
  * Barnes et al., (2020) Barnes, L. P., Chen, W.-N., and Özgür, A. (2020).  Fisher information under local differential privacy.  IEEE Journal on Selected Areas in Information Theory, 1(3):645–659.
  * Bhaila et al., (2024) Bhaila, K., Huang, W., Wu, Y., and Wu, X. (2024).  Local differential privacy in graph neural networks: a reconstruction approach.  In Proceedings of the 2024 SIAM International Conference on Data Mining (SDM), pages 1–9, Texas, US. SIAM, SIAM.
  * Casella and Berger, (2021) Casella, G. and Berger, R. L. (2021).  Statistical inference.  Cengage Learning.
  * Cormode and Bharadwaj, (2022) Cormode, G. and Bharadwaj, A. (2022).  Sample-and-threshold differential privacy: Histograms and applications.  In International Conference on Artificial Intelligence and Statistics, pages 1420–1431, Valencia, Spain. PMLR.
  * Cormode et al., (2018) Cormode, G., Kulkarni, T., and Srivastava, D. (2018).  Marginal release under local differential privacy.  In Proceedings of the 2018 International Conference on Management of Data, pages 131–146.
  * da Costa Filho and Machado, (2023) da Costa Filho, J. S. and Machado, J. C. (2023).  Felip: A local differentially private approach to frequency estimation on multidimensional datasets.  In EDBT, pages 671–683.
  * Du et al., (2021) Du, L., Zhang, Z., Bai, S., Liu, C., Ji, S., Cheng, P., and Chen, J. (2021).  Ahead: adaptive hierarchical decomposition for range query under local differential privacy.  In Proceedings of the 2021 ACM SIGSAC Conference on Computer and Communications Security, pages 1266–1288.
  * Dwork, (2006) Dwork, C. (2006).  Differential privacy.  In International colloquium on automata, languages, and programming, pages 1–12. Springer.
  * Gill and Levit, (1995) Gill, R. D. and Levit, B. Y. (1995).  Applications of the van Trees Inequality: A Bayesian Cramér-Rao Bound.  Bernoulli, 1(1/2):59–79.
  * Jia and Gong, (2019) Jia, J. and Gong, N. Z. (2019).  Calibrate: Frequency estimation and heavy hitter identification with local differential privacy via incorporating prior knowledge.  In IEEE INFOCOM 2019-IEEE Conference on Computer Communications, pages 2008–2016. IEEE.
  * Joseph et al., (2019) Joseph, M., Kulkarni, J., Mao, J., and Wu, S. Z. (2019).  Locally private Gaussian estimation.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
  * Kairouz et al., (2016) Kairouz, P., Oh, S., and Viswanath, P. (2016).  Extremal mechanisms for local differential privacy.  The Journal of Machine Learning Research, 17(1):492–542.
  * Kasiviswanathan et al., (2011) Kasiviswanathan, S. P., Lee, H. K., Nissim, K., Raskhodnikova, S., and Smith, A. (2011).  What can we learn privately?  SIAM Journal on Computing, 40(3):793–826.
  * Li and Miklau, (2012) Li, C. and Miklau, G. (2012).  An adaptive mechanism for accurate query answering under differential privacy.  Proc. VLDB Endow., 5(6):514–525.
  * Lopuhaä-Zwakenberg et al., (2022) Lopuhaä-Zwakenberg, M., Škorić, B., and Li, N. (2022).  Fisher information as a utility metric for frequency estimation under local differential privacy.  In Proceedings of the 21st Workshop on Privacy in the Electronic Society, pages 41–53.
  * Mazumdar et al., (2020) Mazumdar, E., Pacchiano, A., Ma, Y.-A., Bartlett, P. L., and Jordan, M. I. (2020).  On approximate Thompson sampling with Langevin algorithms.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org.
  * McKenna et al., (2020) McKenna, R., Maity, R. K., Mazumdar, A., and Miklau, G. (2020).  A workload-adaptive mechanism for linear queries under local differential privacy.  Proc. VLDB Endow., 13(12):1905–1918.
  * Pelekis and Ramon, (2017) Pelekis, C. and Ramon, J. (2017).  Hoeffding’s inequality for sums of dependent random variables.  Mediterranean Journal of Mathematics, 14(6):243.
  * Qardaji et al., (2013) Qardaji, W., Yang, W., and Li, N. (2013).  Differentially private grids for geospatial data.  In 2013 IEEE 29th international conference on data engineering (ICDE), pages 757–768. IEEE.
  * Russo et al., (2018) Russo, D., Roy, B. V., Kazerouni, A., Osband, I., and Wen, Z. (2018).  A tutorial on thompson sampling.  Foundations and Trends in Machine Learning, 11(1):1–96.
  * Steinberger, (2024) Steinberger, L. (2024).  Efficiency in local differential privacy.
  * (24) Wang, M., Jiang, H., Peng, P., and Li, Y. (2024a).  Accurately estimating frequencies of relations with relation privacy preserving in decentralized networks.  IEEE Transactions on Mobile Computing, 23(05):6408–6422.
  * Wang et al., (2016) Wang, S., Huang, L., Wang, P., Nie, Y., Xu, H., Yang, W., Li, X.-Y., and Qiao, C. (2016).  Mutual information optimally local private discrete distribution estimation.
  * (26) Wang, S., Li, Y., Zhong, Y., Chen, K., Wang, X., Zhou, Z., Peng, F., Qian, Y., Du, J., and Yang, W. (2024b).  Locally private set-valued data analyses: Distribution and heavy hitters estimation.  IEEE Transactions on Mobile Computing [preprint], pages 1–14.
  * Wang et al., (2017) Wang, T., Blocki, J., Li, N., and Jha, S. (2017).  Locally differentially private protocols for frequency estimation.  In 26th USENIX Security Symposium (USENIX Security 17), pages 729–745.
  * Wang et al., (2020) Wang, T., Lopuhaä-Zwakenberg, M., Li, Z., Skoric, B., and Li, N. (2020).  Locally differentially private frequency estimation with consistency.  In 27th Annual Network and Distributed System Security Symposium, NDSS 2020.  Cited by: 33; All Open Access, Bronze Open Access, Green Open Access.
  * Waudby-Smith et al., (2023) Waudby-Smith, I., Wu, S., and Ramdas, A. (2023).  Nonparametric extensions of randomized response for private confidence sets.  In International Conference on Machine Learning, pages 36748–36789. PMLR.
  * Wei et al., (2024) Wei, F., Bao, E., Xiao, X., Yang, Y., and Ding, B. (2024).  Aaa: an adaptive mechanism for locally differential private mean estimation.
  * Welling and Teh, (2011) Welling, M. and Teh, Y. W. (2011).  Bayesian learning via stochastic gradient langevin dynamics.  In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, page 681–688, Madison, WI, USA. Omnipress.
  * Yıldırım, (2024) Yıldırım, S. (2024).  Differentially private online bayesian estimation with adaptive truncation.  Turkish Journal of Electrical Engineering and Computer Sciences, 32(2):34–50.
  * Zhang et al., (2018) Zhang, Z., Wang, T., Li, N., He, S., and Chen, J. (2018).  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, pages 212–229.
  * Zhao et al., (2023) Zhao, D., Zhao, S.-Y., Chen, H., Liu, R.-X., Li, C.-P., and Zhang, X.-Y. (2023).  Hadamard encoding based frequent itemset mining under local differential privacy.  Journal of Computer Science and Technology, 38(6):1403–1422.
  * Zhu et al., (2024) Zhu, Y., Cao, Y., Xue, Q., Wu, Q., and Zhang, Y. (2024).  Heavy hitter identification over large-domain set-valued data with local differential privacy.  IEEE Transactions on Information Forensics and Security, 19:414–426.

##  Appendix A Appendix

###  A.1 Proofs for LDP of RRRR

###### Proof of Theorem 1.

Let $k=|S|$. We can write as

$g_{S,\epsilon}(y|x)=\begin{cases}\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}&x\in
S,y\in S,x=y\\\ \frac{1}{e^{\epsilon_{1}}+k}&x\in S,y\in S,x\neq y\\\
\frac{1}{K-k}\frac{1}{e^{\epsilon_{1}}+k}&x\in S,y\notin S\\\
\frac{1}{e^{\epsilon_{1}}+k}&x\notin S,y\in S\\\
\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}&x\notin
S,y\notin S,x=y\\\
\frac{1}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}&x\notin
S,y\notin S,x\neq y\\\ \end{cases}.$ (21)

We will show that when $\epsilon_{1},\epsilon_{2}$ are chosen according to the
theorem,

$e^{-\epsilon}\leq\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}\leq
e^{\epsilon}$ (22)

for all possible $x,x^{\prime},y\in\\{1,\ldots,K\\}$. When $S=\emptyset$, the
proof is trivial; we focus on the non-trivial case $S\neq\emptyset$. For the
non-trivial case, the transition probability $g_{S,\epsilon}(y|x)$ requires
checking the ratio in (22) in $10$ different cases for $x,x^{\prime},y$
concerning their interrelation.

  1. (C1)

$x\in S$, $x^{\prime}\notin S$, $y\in S$, $y=x$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}}{\frac{1}{e^{\epsilon_{1}}+k}}=e^{\epsilon_{1}}.$

Since $\epsilon_{1}<\epsilon$, (22) holds.

  2. (C2)

$x\in S$, $x^{\prime}\notin S$, $y\in S$, $y\neq x$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{\frac{1}{e^{\epsilon_{1}}+k}}{\frac{1}{e^{\epsilon_{1}}+k}}=1,$

which trivially implies (22).

  3. (C3)

$x\in S$, $x^{\prime}\notin S$, $y\notin S$, $y=x^{\prime}$. We need

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{\frac{1}{K-k}\frac{1}{e^{\epsilon_{1}}+k}}{\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}}=\frac{e^{\epsilon_{2}}+K-k-1}{(K-k)e^{\epsilon_{1}+\epsilon_{2}}}.$

We can show that $\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}\leq
1\leq e^{\epsilon}$ already holds since

$\frac{e^{\epsilon_{2}}+K-k-1}{(K-k)e^{\epsilon_{1}+\epsilon_{2}}}=\frac{(K-k-1)+e^{\epsilon_{2}}}{(K-k-1)e^{\epsilon_{1}+\epsilon_{2}}+e^{\epsilon_{1}+\epsilon_{2}}},$

and the first and the second terms in the numerator are smaller than those in
the denominator, respectively. For the other side of the inequality,

$\frac{e^{\epsilon_{2}}+K-k-1}{(K-k)e^{\epsilon_{1}+\epsilon_{2}}}\geq
e^{-\epsilon}$

implies

$e^{\epsilon_{2}}\leq\frac{K-k-1}{e^{\epsilon_{1}-\epsilon}(K-k)-1}$

whenever $e^{\epsilon_{1}-\epsilon}(K-k)-1>0$, which is the condition given in
the theorem.

  4. (C4)

$x\in S$, $x^{\prime}\notin S$, $y\notin S$, $y\neq x^{\prime}$. We need

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{\frac{1}{K-k}\frac{1}{e^{\epsilon_{1}}+k}}{\frac{1}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}}=\frac{e^{\epsilon_{2}}+(K-k)-1}{(K-k)e^{\epsilon_{1}}}.$

Since $\epsilon_{2}\leq\epsilon$, we have

$e^{\epsilon_{2}}\leq(K-k)(e^{\epsilon+\epsilon_{1}}-1)+1.$

Hence

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}\leq\frac{(K-k)(e^{\epsilon+\epsilon_{1}}-1)+1+K-k-1}{(K-k)e^{\epsilon_{1}}}\leq
e^{\epsilon},$

hence we proved the right-hand side inequality. For the left-hand side, we
have

$\frac{e^{\epsilon_{2}}+K-k-1}{(K-k)e^{\epsilon_{1}}}=\frac{(e^{\epsilon_{2}}-1)+K-k}{(K-k)e^{\epsilon_{1}}}\geq\frac{K-k}{(K-k)e^{\epsilon_{1}}}=e^{-\epsilon_{1}}\geq
e^{-\epsilon}$

since $\epsilon_{2}\geq 0$ and $\epsilon_{1}\leq\epsilon$.

  5. (C5)

$x,x^{\prime}\in S$, $y\in S$, $y=x$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{e^{\epsilon_{1}}/(e^{\epsilon_{1}}+k)}{1/(e^{\epsilon_{1}}+k)}=e^{\epsilon_{1}}.$

Since $\epsilon_{1}<\epsilon$, (22) holds.

  6. (C6)

$x,x^{\prime}\in S$, $y\in S$, $y\neq x$ and $y\neq x^{\prime}$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{1/(e^{\epsilon_{1}}+k)}{1/(e^{\epsilon_{1}}+k)}=1.$

So (22), trivially holds.

  7. (C7)

$x,x^{\prime}\in S$, $y\notin S$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{\frac{1}{K-k}\frac{1}{e^{\epsilon_{1}}+k}}{\frac{1}{K-k}\frac{1}{e^{\epsilon_{1}}+k}}=1.$

So, (22) trivially holds.

  8. (C8)

$x,x^{\prime}\notin S$, $y\notin S$, $y=x$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{e^{\epsilon_{2}}/(e^{\epsilon_{2}}+K-k-1)e^{\epsilon_{1}}/(e^{\epsilon_{1}}+k)}{1/(e^{\epsilon_{2}}+K-k-1)e^{\epsilon_{1}}/(e^{\epsilon_{1}}+k)}=e^{\epsilon_{2}}.$

Since $\epsilon_{2}\leq\epsilon$, (22) holds.

  9. (C9)

$x,x^{\prime}\notin S$, $y\notin S$, $y\neq x$, $y\neq x^{\prime}$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{1/(e^{\epsilon_{2}}+K-k-1)e^{\epsilon_{1}}/(e^{\epsilon_{1}}+k)}{1/(e^{\epsilon_{2}}+K-k-1)e^{\epsilon_{1}}/(e^{\epsilon_{1}}+k)}=1.$

So, (22) trivially holds.

  10. (C10)

$x,x^{\prime}\notin S$, $y\in S$. We have

$\frac{g_{S,\epsilon}(y|x)}{g_{S,\epsilon}(y|x^{\prime})}=\frac{1/(e^{\epsilon_{1}}+|S|)}{1/(e^{\epsilon_{1}}+|S|)}=1.$

So (22), trivially holds.

We conclude the proof by noting that any other case left out is symmetric in
$(x,x^{\prime})$ to one of the covered cases and therefore does not need to be
checked separately. ∎

###  A.2 Proofs about utility functions

###### Proof of Proposition 1.

Given $\theta\in\Delta$, let $\vartheta$ be the $(K-1)\times 1$ column vector
such that $\vartheta_{i}=\theta_{i}$ for $i=1,\ldots,K-1$. We can write the
Fisher information matrix in terms of the score vector as follows.

$F(\theta;S,\epsilon)=\mathbb{E}_{Y}\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y|\theta)\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y|\theta)^{\top}\right]=\sum_{y=1}^{K}h_{S,\epsilon}(y|\theta)\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}\right].$

Noting that

$h_{S,\epsilon}(y|\theta)=\sum_{k=1}^{K-1}g_{S,\epsilon}(y|k)\vartheta_{k}+g_{S,\epsilon}(y|K)\left(1-\sum_{k=1}^{K-1}\vartheta_{k}\right),$

the score vector can be derived as

$[\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)]_{k}=\frac{g_{S,\epsilon}(y|k)-g_{S,\epsilon}(y|K)}{h_{S,\epsilon}(y|\theta)},\quad
k=1,\ldots,K-1.$ (23)

As the $K\times(K-1)$ matrix $A_{S,\epsilon}$ defined as
$A(i,j)=g(i|j)-g(i|K)$, we can rewrite (23) as $[\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)]_{k}=A_{S,\epsilon}(y,k)/h_{S,\epsilon}(y|\theta)$.
Let $a_{y}$ be the $y$’th row of $A_{S,\epsilon}$, and recall that
$D_{\theta}$ is defined as a diagonal matrix with $1/h_{S,\epsilon}(j|\theta)$
being the $j$’th element in the diagonal. Then, the Fisher information matrix
is

$F(\theta;S,\epsilon)=\sum_{y=1}^{K}\frac{a_{y}^{\top}}{h_{S,\epsilon}(y|\theta)}\frac{a_{y}}{h_{S,\epsilon}(y|\theta)}h_{S,\epsilon}(y|\theta)=\sum_{y=1}^{K}a_{y}^{\top}\frac{1}{h_{S,\epsilon}(y|\theta)}a_{y}=A_{S,\epsilon}^{\top}D_{\theta}A_{S,\epsilon}$

as claimed. ∎

Next, we prove that $F(\theta;S,\epsilon)$ is invertible. Let $G_{S,\epsilon}$
be the $K\times K$ matrix whose elements are

$G_{S,\epsilon}(i,j)=g_{S,\epsilon}(i|j),\quad i,j=1,\ldots,K.$ (24)

To prove that $F(\theta;S,\epsilon)$ is invertible, we first prove the
intermediate result that $G_{S,\epsilon}$ is invertible.

######  Lemma 1.

$G_{S,\epsilon}$ is invertible for all $S\subset\\{1,\ldots,K\\}$ and
$\epsilon>0$.

###### Proof.

It suffices to prove that of $G_{S,\epsilon}$ is invertible for
$S=\\{1,2,\ldots,k\\}$ and for all $k\in\\{0,\ldots,K-1\\}$. For other $S$,
$G_{S,\epsilon}$ can be obtained by permutation. Fix $k$ and let
$S=\\{1,2,\ldots,k\\}$. It can be verified by inspection that $G_{S,\epsilon}$
is a block matrix as

$G_{S,\epsilon}=\begin{bmatrix}a_{1}I_{k}+a_{2}1_{k}1_{k}^{\top}&b1_{k}1_{K-k}^{\top}\\\
c1_{K-k}1_{k}^{\top}&d_{1}I_{K-k}+d_{2}1_{K-k}1_{K-k}^{\top}\end{bmatrix},$

where $I_{n}$ is the identity matrix of size $n$ and $1_{n}$ is the column
vector of $1$’s of size $n$. The constants $a_{1},a_{2},b,c,d_{1},d_{2}$ are
given as

$\displaystyle a_{1}=\frac{e^{\epsilon_{1}}}{k+e^{\epsilon_{1}}}-a_{2},\quad
a_{2}=\frac{1}{k+e^{\epsilon_{1}}},\quad b=\frac{1}{k+e^{\epsilon_{1}}},\quad
c=\frac{1}{K-k}a_{2}$ $\displaystyle
d_{1}=\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{k+e^{\epsilon_{1}}}-d_{2},\quad
d_{2}=\frac{1}{e^{\epsilon_{2}}+K-k-1}\frac{e^{\epsilon_{1}}}{k+\epsilon^{1}}.$

Also, note that since $\epsilon_{1}>0$ and $\epsilon_{2}>0$, $a_{1}$ and
$d_{1}$ (whenever it is defined) are strictly positive.

The case $k=0$ is trivial since then
$G_{S,\epsilon}=d_{1}I_{K}+d_{2}1_{K}1_{K}^{\top}$ is invertible. Hence, we
focus on the case $0<k<K$. For this case, firstly note that the matrices on
the diagonal are invertible. So, by Weinstein–Aronszajn identity, for
$G_{S,\epsilon}$ to be invertible, it suffices to show that the matrix

$M=a_{1}I_{k}+a_{2}1_{k}1_{k}^{\top}-b1_{k}1_{K-k}^{\top}(d_{1}I_{K-k}+d_{2}1_{K-k}1_{K-k}^{\top})^{-1}c1_{K-k}1_{k}^{\top}$

is invertible. Using the Woodbury matrix identity, the matrix $M$ can be
expanded as

$\displaystyle M$
$\displaystyle=a_{1}I_{k}+a_{2}1_{k}1_{k}^{\top}-b1_{k}1_{K-k}^{\top}\left(\frac{I_{K-k}}{d_{1}}-\frac{1}{d_{1}}1_{K-k}\left(\frac{1}{d_{2}}+1_{K-k}^{\top}\frac{1}{d_{1}}1_{K-k}\right)^{-1}1_{K-k}^{\top}\frac{1}{d_{1}}\right)c1_{K-k}1_{k}^{\top}$
$\displaystyle=a_{1}I_{k}+a_{2}1_{k}1_{k}^{\top}-\frac{bc}{d_{1}}1_{k}1_{K-k}^{\top}1_{K-k}1_{k}^{\top}+\frac{bc}{d_{1}^{2}}\left(\frac{1}{d_{2}}+\frac{K-k}{d_{1}}\right)^{-1}1_{k}1_{K-k}^{\top}1_{K-k}1_{K-k}^{\top}1_{K-k}1_{k}^{\top}$
$\displaystyle=a_{1}I_{k}+\left[a_{2}-\frac{(K-k)bc}{d_{1}}+\frac{bc}{d_{1}^{2}}\left(\frac{1}{d_{2}}+\frac{K(K-k)}{d_{1}}\right)^{-1}\right]1_{k}1_{k}^{\top}.$

Inside the square brackets is a scalar, therefore $M$ in question is the sum
of an identity matrix and a rank-$1$ matrix, which is invertible. Hence,
$G_{S,\epsilon}$ is invertible. ∎

###### Proof of Proposition 2.

Note that $A_{S,\epsilon}=G_{S,\epsilon}J$, where the $K\times(K-1)$ matrix
$J$ satisfies $J(i,i)=1$ and $J(i,K)=-1$ for $i=1,\ldots,K$, and $J(i,j)=0$
otherwise. Since $G_{S,\epsilon}$ is invertible, it is full rank. Also, the
columns of $A_{S,\epsilon}$, denoted by $c^{A}_{i}$, $i=1,\ldots,K-1$ are
given by

$c^{A}_{1}=c^{G}_{1}-c^{G}_{K},\quad\ldots,\quad
c^{A}_{K-1}=c^{G}_{K-1}-c^{G}_{K},$

where $c^{G}_{i}$ is the $i$’th column of $G_{S,\epsilon}$ for $i=1,\ldots,K$.
Observe that $c^{A}_{i}$, $i=1,\ldots,K-1$ are linearly independent, since any
linear combination of those columns is in the form of

$\sum_{i=1}^{K-1}a_{i}c^{A}_{i}=\sum_{i=1}^{K-1}a_{i}c^{G}_{i}-\left(\sum_{i=1}^{K-1}a_{i}\right)c^{G}_{K}.$

Since the columns of $G_{S,\epsilon}$ are linearly independent, the linear
combination above becomes $0$ only if $a_{1}=\ldots=a_{K-1}=0$. This shows
that the columns of $A_{S,\epsilon}$ are also linearly independent, thus we
conclude that $A_{S,\epsilon}$ has rank $K-1$. Finally, since $D_{\theta}$ is
diagonal with positive diagonal entries,
$A_{S,\epsilon}^{\top}D_{\theta}A_{S,\epsilon}=A_{S,\epsilon}^{\top}D_{\theta}^{1/2}D_{\theta}^{1/2}A_{S,\epsilon}$
is positive definite, hence invertible. ∎

The following proof contains a derivation of the utility function based on the
MSE of the Bayesian estimator of $X$ given $Y$.

###### Proof of Proposition 3.

It is well-known that the expectation in (9) is minimized when
$\widehat{e_{X}}=\nu(Y):=\mathbb{E}_{\theta}[e_{X}|Y]$, i.e. the posterior
expectation of $e_{X}$ given $Y$. That is,

$\min_{\widehat{e_{X}}}\mathbb{E}_{\theta}\left[\|e_{X}-\widehat{e_{X}}(Y)\|^{2}\right]=\mathbb{E}_{\theta}\left[\|e_{X}-\nu(Y)\|^{2}\right].$

For the squared norm inside the expectation, we have

$\displaystyle\|e_{X}-\nu(Y)\|^{2}$
$\displaystyle=(1-v(Y)_{X})^{2}+\sum_{k\neq X}v(Y)_{k}^{2}$
$\displaystyle=1+\nu(Y)_{X}^{2}-2\nu(Y)_{X}+\sum_{k\neq X}v(Y)_{k}^{2}$
$\displaystyle=1-2\nu(Y)_{X}+\sum_{k=1}^{K}v(Y)_{k}^{2}.$ (25)

The expectation of the last term in (25) is

$\displaystyle\mathbb{E}_{\theta}\left[\sum_{k=1}^{K}v(Y)_{k}^{2}\right]$
$\displaystyle=\sum_{y}h_{S,\epsilon}(y|\theta)\sum_{x=1}^{K}p_{S,\epsilon}(x|y,\theta)^{2}$
$\displaystyle=\sum_{y=1}^{K}\sum_{x=1}^{K}h_{S,\epsilon}(y|\theta)p_{S,\epsilon}(x|y,\theta)^{2}$
$\displaystyle=\sum_{y=1}^{K}\sum_{x=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)}.$
(26)

For the expectation of the second term in (25), we have

$\mathbb{E}_{\theta}\left[\nu(Y)_{X}\right]=\sum_{x,y}p_{S,\epsilon}(x|y,\theta)p_{S,\epsilon}(x,y|\theta),$

where $p(x,y|\theta)$ denotes the joint probability of $X,Y$ given $\theta$.
Substituting
$p(x,y|\theta)=p_{S,\epsilon}(x|y,\theta)h_{S,\epsilon}(y|\theta)$ into the
equation above, we get

$\displaystyle\mathbb{E}_{\theta}\left[\nu(Y)_{X}\right]$
$\displaystyle=\sum_{x=1}^{K}\sum_{y=1}^{K}h_{S,\epsilon}(y|\theta)p_{S,\epsilon}(x|y,\theta)^{2}.$
$\displaystyle=\sum_{x=1}^{K}\sum_{y=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)},$
(27)

which is equal to what we get in (26). Substituting (26) and (27) into (25),
we obtain

$\displaystyle\mathbb{E}_{\theta}\left[\|e_{X}-\nu(Y)\|^{2}\right]$
$\displaystyle=1-2\sum_{x=1}^{K}\sum_{y=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)}+\sum_{x=1}^{K}\sum_{y=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)}$
$\displaystyle=1-\sum_{x=1}^{K}\sum_{y=1}^{K}\frac{g_{S,\epsilon}(y|x)^{2}\theta_{x}^{2}}{h_{S,\epsilon}(y|\theta)}.$

Finally, using the definition
$U_{5}(\theta,S,\epsilon)=-\min_{\widehat{e_{X}}}\mathbb{E}_{\theta}\left[\|e_{X}-\widehat{e_{X}}(Y)\|^{2}\right]$,
we conclude the proof. ∎

###### Proof of Theorem 2.

The global maximization of $U_{6}$ over all $S\subset\\{1,\ldots,K\\}$ can be
decomposed as

$\displaystyle\max_{S\subset\\{1,\ldots,K\\}}U_{6}(\theta,S,\epsilon)=\max_{k\in\\{0,\ldots,K-1\\}}\left\\{\max_{S\subset\\{1,\ldots,K\\}:|S|=k}U_{6}(\theta,S,\epsilon)\right\\}.$
(28)

This inner maximization is equivalent to fixing the cardinality of $S$ to $k$
and finding the best $S$ with cardinality $k$. Now, the utility function can
be written as

$\displaystyle U_{6}(\theta,S,\epsilon)$
$\displaystyle=\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}\left(\sum_{i\in
S}\theta_{i}+\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-k-1}\sum_{i\notin
S}\theta_{i}\right)$
$\displaystyle=\left(\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}\right)\sum_{i\in
S}\theta_{i}+\left(\frac{e^{\epsilon_{1}}}{e^{\epsilon_{1}}+k}\right)\left(\frac{e^{\epsilon_{2}}}{e^{\epsilon_{2}}+K-k-1}\right)\sum_{i\notin
S}\theta_{i},$

where $k$ appears in the first line since $|S|=k$. Note that $\sum_{i\in
S}\theta_{i}$ and $\sum_{i\notin S}\theta_{i}$ sum to $1$ and the constants in
front of the first sum is larger than that of the second. Hence, we seek to
maximize an expression in the form of

$ax+b(1-x)$

over a variable $x>0$ when $a>b>0$. This is maximized when $x>0$ is taken as
large as possible. Therefore, $U_{6}(\theta,S,\epsilon)$ is maximized when
$\sum_{i\in S}\theta_{i}$ is made as large as possible under the constraint
that $|S|=k$. Under this constraint this sum is maximized when $S$ has the
indices of the $k$ largest components of $\theta$, that is, when
$S=S_{k,\theta}=\\{\sigma_{\theta}(1),\ldots,\sigma_{\theta}(k)\\}$.

Then, (28) reduces to
$\max_{k=1,\ldots,K}U_{6}(\theta,S_{k,\theta},\epsilon)$. Hence, we conclude.
∎

###  A.3 Proof for SGLD update

###### Proof of Proposition 4.

For the prior component of the gradient, recall that we have
$\phi=(\phi_{1},\ldots,\phi_{K})$, where

$\phi_{i}\overset{\text{iid}}{\sim}\text{Gamma}(\rho_{i},1),\quad
i=1,\ldots,K.$

Then, the marginal pdf of $\phi_{i}$ satisfies

$\ln p(\phi_{i})=(\rho_{i}-1)\ln\phi_{i}-\ln\Gamma(\rho_{i})-\phi_{i},\quad
i=1,\ldots,K.$

Taking the partial derivatives of $\ln p(\phi_{i})$ with respect to
$\phi_{i}$, we have

$[\nabla_{\phi}\ln p(\phi)]_{i}=\frac{\rho_{i}-1}{\phi_{i}}-1,\quad
i=1,\ldots,K.$

For the likelihood component, given $\theta\in\Delta$, let the $(K-1)\times 1$
vector $\vartheta$ be the reparametrization of $\theta$ such that
$\vartheta_{i}=\theta_{i}$ for $i=1,\ldots,K-1$. Then, according to (12),

$\vartheta_{k}=\frac{\phi_{k}}{\sum_{j=1}^{K}\phi_{j}},\quad k=1,\ldots,K-1.$
(29)

Using the chain rule, we can write the gradient of the log-likelihood with
respect to $\phi$ as

$\displaystyle\nabla_{\phi}\ln p(y|\phi)$
$\displaystyle=J\cdot\nabla_{\vartheta}\ln h_{S,\epsilon}(y|\theta).$

where $J$ is the $K\times(K-1)$ Jacobian matrix for the mapping from $\phi$ to
$\vartheta$ in (29), whose $(i,j)$th element can be derived as

$\displaystyle J(i,j)$
$\displaystyle=\frac{\partial\vartheta_{j}}{\partial\phi_{i}}=\frac{\partial}{\partial\phi_{i}}\frac{\phi_{j}}{\sum_{k=1}^{K}\phi_{k}}$
$\displaystyle=\mathbb{I}(i=j)\frac{1}{\sum_{k=1}^{K}\phi_{k}}-\frac{\phi_{j}}{\left(\sum_{k=1}^{K}\phi_{k}\right)^{2}}.$

Using (23) for $\nabla_{\vartheta}\ln h_{S,\epsilon}(y|\theta)$, we complete
the proof. ∎

###  A.4 Proofs for convergence and consistency results

####  A.4.1 Preliminary results

######  Lemma 2.

Given $\epsilon\geq 0$, there exists constants $0<c\leq C<\infty$ such that
for all $\theta\in\Delta$ and all $S\subset\\{1,\ldots,K\\}$, we have

$c\leq g_{S,\epsilon}(y|x)\leq C,\quad c\leq h_{S,\epsilon}(y|\theta)\leq C.$

###### Proof.

The bounds for $g_{S,\epsilon}(y|x)$ can directly be verified from (21).
Moreover,

$c\leq\min_{i=1,\ldots,K}g_{S,\epsilon}(y|i)\leq
h_{S,\epsilon}(y|\theta)=\sum_{i=1}^{K}g_{S,\epsilon}(y|i)\theta_{i}\leq\max_{i=1,\ldots,K}g_{S,\epsilon}(y|i)\leq
C.$

Hence, we conclude. ∎

######  Remark 1.

For the symbol $\theta$, which is used for $K\times 1$ probability vectors in
$\Delta$, we will associate the symbol $\vartheta$ such that $\vartheta$
denotes the shortened vector of the first $K-1$ elements of $\theta$.
Accordingly, we will use $\vartheta$, $\vartheta^{\prime}$,
$\vartheta^{\ast}$, etc, to denote the shortened versions of $\theta$,
$\theta^{\prime}$, $\theta^{\ast}$, etc.

######  Lemma 3.

For any $\theta,\theta^{\prime}\in\Delta$, $y\in\\{1,\ldots,K\\}$, and
$S\subset\\{1,\ldots,K\\}$, we have

$\displaystyle\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})$
$\displaystyle=\frac{h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}.$
(30)

###### Proof.

Recall from (23) that

$[\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)]_{i}=\frac{g_{S,\epsilon}(y|i)-g_{S,\epsilon}(y|K)}{h_{S,\epsilon}(y|\theta)},\quad
i=1,\ldots,K-1.$

Hence,

$\displaystyle\nabla_{\vartheta}\ln h_{S,\epsilon}(y|\theta)^{\top}\vartheta$
$\displaystyle=\frac{1}{h_{S,\epsilon}(y|\theta)}\sum_{i=1}^{K-1}(g_{S,\epsilon}(y|i)-g_{S,\epsilon}(y|K))\vartheta_{i}$
$\displaystyle=\frac{1}{h_{S,\epsilon}(y|\theta)}\left[\sum_{i=1}^{K-1}g_{S,\epsilon}(y|i)\vartheta_{i}-g_{S,\epsilon}(y|K)\sum_{i=1}^{K-1}\vartheta_{i}\right]$
$\displaystyle=\frac{1}{h_{S,\epsilon}(y|\theta)}\left[\sum_{i=1}^{K-1}g_{S,\epsilon}(y,i)\theta_{i}-g_{S,\epsilon}(y|K)(1-\theta_{K})\right]$
$\displaystyle=\frac{1}{h_{S,\epsilon}(y|\theta)}\left[h_{S,\epsilon}(y|\theta)-g_{S,\epsilon}(y|K)\right],$

and likewise $\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}\vartheta^{\prime}=\frac{1}{h_{S,\epsilon}(y|\theta)}\left[h_{S,\epsilon}(y|\theta^{\prime})-g_{S,\epsilon}(y|K)\right]$.
Taking the difference between $\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}\vartheta$ and $\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}\vartheta^{\prime}$, we arrive at the result. ∎

##### Concavity of $\ln h_{S,\epsilon}(y|\theta)$:

The following lemmas help with proving the concavity of
$h_{S,\epsilon}(y|\theta)$ as a function of $\theta$.

######  Lemma 4.

For $0<b\leq a\leq 1$, we have
$\ln\frac{a}{b}\geq\frac{a-b}{a}+\frac{(a-b)^{2}}{2}$.

###### Proof.

For $z>0$, using the series based on the inverse hyperbolic tangent function,
we can write

$\ln z=2\sum_{k=0}^{\infty}\frac{1}{2k+1}\left(\frac{z-1}{z+1}\right)^{2k+1}.$

Apply the expansion to $z=a/b$ when $a\geq b$. Noting that
$(z-1)/(z+1)=(a-b)/(a+b)$,

$\displaystyle\ln\frac{a}{b}$ $\displaystyle\geq
2\frac{a-b}{a+b}\geq\frac{2(a-b)}{2a}=\frac{a-b}{a}$

where the difference is

$\frac{2(a-b)}{a+b}-\frac{a-b}{a}=\frac{2a^{2}-2ab-a^{2}+b^{2}}{a(a+b)}=\frac{(a-b)^{2}}{(a+b)a}\geq\frac{(a-b)^{2}}{2}$

since $0<a,b\leq 1$. ∎

######  Lemma 5.

Let $0<\alpha\leq 1$. For $x\geq 1$, we have
$\frac{1}{\alpha}(x^{\alpha}-1)\leq(x-1)$.

###### Proof.

Consider $f(x)=\frac{1}{\alpha}(x^{\alpha}-1)-(x-1)$. We have $f(1)=0$ and
$f^{\prime}(x)=x^{\alpha-1}-1\leq 0$ for $x\geq 1$. Hence we conclude. ∎

######  Lemma 6.

Let $0<a\leq b\leq 1$ and $0<\alpha\leq 1$. Then,
$b^{\alpha}-a^{\alpha}\geq\alpha(b-a)$.

###### Proof.

Fix $a$ and let $b=a+x$ for $0\leq x\leq 1-a$. Consider the function
$f(x)=(a+x)^{\alpha}-a^{\alpha}-\alpha x$. We have $f(0)=0$ and
$f^{\prime}(x)=\alpha(a+x)^{\alpha-1}-\alpha\geq 0$ over $0\leq x\leq(1-a)$
since $a+x\leq 1$ and $\alpha-1\leq 0$. Hence, we conclude. ∎

######  Lemma 7.

Given $\epsilon>0$, there exists $m_{0}>0$ such that, for all
$S\subset\\{1,\ldots,K\\}$ and $y\in\\{1,\ldots,K\\}$, $\ln
h_{S,\epsilon}(y|\theta)$ is a concave function of $\vartheta$ that satisfies

$\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})\geq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+m_{0}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2}$

for all $\theta,\theta^{\prime}\in\Delta$.

###### Proof.

We will look at the cases $h_{S,\epsilon}(y|\theta)\geq
h_{S,\epsilon}(y|\theta^{\prime})$ and $h_{S,\epsilon}(y|\theta)\leq
h_{S,\epsilon}(y|\theta^{\prime})$ separately.

  1. 1.

Assume $h_{S,\epsilon}(y|\theta)\geq h_{S,\epsilon}(y|\theta^{\prime})$. Using
Lemma 4, we have

$\displaystyle\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})$
$\displaystyle\geq\frac{h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}+\frac{1}{2}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2}$
$\displaystyle=\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{1}{2}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2},$

where the last line follows from Lemma 3.

  2. 2.

Assume $h_{S,\epsilon}(y|\theta)\leq h_{S,\epsilon}(y|\theta^{\prime})$. Let
$a=h_{S,\epsilon}(y|\theta)$ and $b=h_{S,\epsilon}(y|\theta^{\prime})$. Let

$\alpha=\min\left\\{1,\frac{\ln 2}{\ln(C/c)}\right\\},$

where $c$ and $C$ are given in Lemma 2. This $\alpha$ ensures that
$0<\alpha\leq 1$ and $b^{\alpha}/a^{\alpha}\leq 2$, so that we can use
Taylor’s expansion of $\ln\frac{b^{\alpha}}{a^{\alpha}}$ around $1$ and have

$\ln\frac{b}{a}=\frac{1}{\alpha}\ln\frac{b^{\alpha}}{a^{\alpha}}=\frac{1}{\alpha}\left[\sum_{k=1}^{\infty}(-1)^{k+1}\frac{1}{k}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{k}\right].$

Approximating the expansion up to its third term, we have the following upper
bound on $\ln\frac{b}{a}$ as

$\ln\frac{b}{a}\leq\frac{1}{\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)-\frac{1}{2\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{2}+\frac{1}{3\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{3}.$
(31)

For the third term, we have

$\frac{1}{3\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{3}=\frac{1}{3\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{2}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)\leq\frac{1}{3\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{2},$

since $1\leq\frac{b^{\alpha}}{a^{\alpha}}\leq 2$. Substituting this into (31),
the inequality can be continued as

$\displaystyle\ln\frac{b}{a}\leq\frac{1}{\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)-\frac{1}{6\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{2}.$
(32)

Using Lemma 5 to bound the first term in (32), we have

$\displaystyle\ln\frac{b}{a}\leq\left(\frac{b}{a}-1\right)-\frac{1}{6\alpha}\left(\frac{b^{\alpha}}{a^{\alpha}}-1\right)^{2}.$
(33)

Finally, using Lemma 6 we can lower-bound the second term in (32) as

$\frac{b^{\alpha}}{a^{\alpha}}-1=\frac{b^{\alpha}-a^{\alpha}}{a^{\alpha}}\geq\frac{\alpha(b-a)}{a^{\alpha}}\geq\alpha(b-a),$

where the last inequality follows from $a\leq 1$ and $\alpha>0$. We end up
with

$\displaystyle\ln\frac{b}{a}\leq\left(\frac{b}{a}-1\right)-\frac{\alpha}{6}\left(b-a\right)^{2}.$

Referring to the definitions of $a$ and $b$, we have

$\ln h_{S,\epsilon}(y|\theta^{\prime})-\ln
h_{S,\epsilon}(y|\theta)\leq\left(\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}-1\right)-\frac{\alpha}{6}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2},$

or, reversing the inequality,

$\displaystyle\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})\geq\left(1-\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}\right)+\frac{\alpha}{6}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2}.$
(34)

Using (30) in Lemma 3, we rewrite (34) as

$\displaystyle\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})\geq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{\alpha}{6}(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime}))^{2},$
(35)

which is the inequality we look for.

To cover both cases, take $m_{0}=\min\\{\frac{\alpha}{6},\frac{1}{2}\\}$. So
the proof is complete. ∎

Recalling that $S_{t}$ is the selected subset at time $t$, define

$V_{t}(\theta,\theta^{\prime}):=(h_{S_{t},\epsilon}(Y_{i}|\theta)-h_{S_{t},\epsilon}(Y_{t}|\theta^{\prime}))^{2}.$
(36)

The proof of Theorem 3 requires a probabilistic bound for
$\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})$, which we provide next.

######  Lemma 8.

For all $\theta,\theta^{\prime}\in\Delta$ and $t\geq 0$,
$V_{t}(\theta,\theta^{\prime})\leq\|\theta-\theta^{\prime}\|^{2}$.

###### Proof.

Let the $1\times K$ vector $r_{i}$ be the $i$th row of $G_{S_{t},\epsilon}$.
Then, we obtain

$\displaystyle
V_{t}(\theta,\theta^{\prime})=[h_{S_{t},\epsilon}(y|\theta)-h_{S_{t},\epsilon}(y|\theta^{\prime})]^{2}$
$\displaystyle=(r_{i}(\theta-\theta^{\prime}))^{2}$
$\displaystyle=(\theta-\theta^{\prime})^{\top}r_{i}^{\top}r_{i}(\theta-\theta^{\prime})$
$\displaystyle\leq\|\theta-\theta^{\prime}\|^{2},$

since every element of $r_{i}$ is at most $1$. ∎

######  Lemma 9.

For all $\theta,\theta^{\prime}\in\Delta$ and $t\geq 1$, there exists a
constant $c_{u}>0$ such that

$E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})]\geq
c_{u}\|\theta-\theta^{\prime}\|^{2},$

where $E_{\theta^{\ast}}$ is the expectation operator with respect to
$P_{\theta^{\ast}}$ defined in (17), $\lambda_{\textup{min}}(A)$ is the
minimum eigenvalue of the square matrix $A$ and the matrix $G_{S,\epsilon}$ is
defined in (24).

###### Proof.

The overall expectation can be written as

$\displaystyle E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})]$
$\displaystyle=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})|S_{t}=S]$

where the conditional expectation can be bounded as

$\displaystyle E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})|S_{t}=S]$
$\displaystyle=\sum_{i=1}^{K}(h_{S,\epsilon}(i|\theta)-h_{S,\epsilon}(i|\theta^{\prime}))^{2}h_{S,\epsilon}(i|\theta^{\ast})$
$\displaystyle\geq\sum_{i=1}^{K}(h_{S,\epsilon}(i|\theta)-h_{S,\epsilon}(i|\theta^{\prime}))^{2}c,$
(37)

where the second line is due to Lemma 2. Further, let the $1\times K$ vector
$r_{i}$ be the $i$th row of $G_{S,\epsilon}$. Then,

$\displaystyle\sum_{i=1}^{K}(h_{S,\epsilon}(i|\theta)-h_{S,\epsilon}(i|\theta^{\prime}))^{2}$
$\displaystyle=\sum_{i=1}^{K}(r_{i}(\theta-\theta^{\prime}))^{2}$
$\displaystyle=\sum_{i=1}^{K}(\theta-\theta^{\prime})^{\top}r_{i}^{\top}r_{i}(\theta-\theta^{\prime})$
$\displaystyle=(\theta-\theta^{\prime})^{\top}G_{S,\epsilon}^{\top}G_{S,\epsilon}(\theta-\theta^{\prime})$
$\displaystyle=\|G_{S,\epsilon}(\theta-\theta^{\prime})\|^{2}$
$\displaystyle\geq\lambda_{\min}(G_{S,\epsilon}^{\top}G_{S,\epsilon})\|\theta-\theta^{\prime}\|^{2}.$
(38)

Combining (37) and (38) and letting
$c_{u}:=c\min_{S\subset\\{1,\ldots,K\\}}\lambda_{\min}(G_{S,\epsilon}^{\top}G_{S,\epsilon})$,
we have

$E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})|S_{t}=S]\geq
c_{u}\|\theta-\theta^{\prime}\|^{2}.$ (39)

for all $S$, which directly implies
$E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})]\geq
c_{u}\|\theta-\theta^{\prime}\|^{2}$ for the overall expectation. Finally,
$c_{u}>0$ since, by Lemma 1, every $G_{S,\epsilon}$ is invertible. ∎

Further, define

$W_{t}(\theta,\theta^{\prime}):=1-\frac{V_{t}(\theta,\theta^{\prime})}{\|\theta-\theta^{\prime}\|^{2}}.$

######  Lemma 10.

For any $0<t_{1}<\ldots<t_{k}$, we have
$E_{\theta^{\ast}}\left[\prod_{i=1}^{k}W_{t_{i}}(\theta,\theta^{\prime})\right]\leq(1-c_{u})^{k}$.

###### Proof.

For simplicity, we drop $(\theta,\theta^{\prime})$ from the notation and
denote the random variables in question as $W_{t_{1}},\ldots,W_{t_{k}}$. We
can write

$\displaystyle E_{\theta^{\ast}}\left(\prod_{i=1}^{k}W_{t_{i}}\right)$
$\displaystyle=E_{\theta^{\ast}}\left[E_{\theta^{\ast}}\left(\left.\prod_{i=1}^{k}W_{t_{i}}\right|W_{t_{1}},\ldots,W_{t_{k-1}}\right)\right]$
$\displaystyle=E_{\theta^{\ast}}\left[\left(\prod_{i=1}^{k-1}W_{t_{i}}\right)E_{\theta^{\ast}}\left(W_{t_{k}}|W_{t_{1}},\ldots,W_{t_{k-1}}\right)\right].$
(40)

By construction of $W_{i}(\theta,\theta^{\prime})$, we have
$E_{\theta^{\ast}}(W_{i}|S_{i}=S)\leq 1-c_{u}$, which follows from (39). Using
this, the inner conditional expectation can be bounded as

$\displaystyle
E_{\theta^{\ast}}\left(W_{t_{k}}|W_{t_{1}},\ldots,W_{t_{k-1}}\right)$
$\displaystyle=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t_{k}}=S|W_{t_{1}},\ldots,W_{t_{k-1}})E_{\theta^{\ast}}(W_{t_{k}}|S_{t_{k}}=S)$
$\displaystyle\leq\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t_{k}}=S|W_{t_{1}},\ldots,W_{t_{k-1}})\left(1-c_{u}\right)$
$\displaystyle=\left(1-c_{u}\right).$ (41)

Combining (40) and (41), we have

$E_{\theta^{\ast}}\left(\prod_{i=1}^{k}W_{t_{i}}\right)\leq\left(1-c_{u}\right)E_{\theta^{\ast}}\left(\prod_{i=1}^{k-1}W_{t_{i}}\right).$
(42)

By Lemmas 8 and 9, we have
$V_{t}(\theta,\theta^{\prime})\leq\|\theta-\theta^{\prime}\|^{2}$ and
$E_{\theta^{\ast}}[V_{t}(\theta,\theta^{\prime})]\geq
c_{u}\|\theta-\theta^{\prime}\|^{2}$. Thus we necessarily have $c_{u}<1$.
Therefore, the recursion in (42) can be used until $k=1$ to obtain the desired
result. ∎

Note that $W_{t}(\theta,\theta^{\prime})$ is bounded as $0\leq
W_{i}(\theta,\theta^{\prime})\leq 1$. We now quote a critical theorem from
Pelekis and Ramon,  (2017, Theorem 3.2) regarding the sum of dependent and
bounded random variables, which will be useful for bounding
$\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})$.

######  Theorem 5.

(Pelekis and Ramon, , 2017, Theorem 3.2) Let $W_{1},\ldots,W_{n}$ be random
variables, such that $0\leq W_{t}\leq 1$, for $t=1,\ldots,n$. Fix a real
number $\tau\in(0,n)$ and let $k$ be any positive integer, such that
$0<k<\tau$. Then

$\displaystyle\mathbb{P}\left(\sum_{t=1}^{n}W_{t}\geq\tau\right)\leq\frac{1}{\binom{\tau}{k}}\sum_{A\subseteq\\{1,\ldots,n\\}:|A|=k}\mathbb{E}\left[\prod_{i\in
A}W_{i}\right],$ (43)

where $\binom{\tau}{k}=\frac{\tau(\tau-1)\ldots(\tau-k+1)}{k!}$.

In the following, we apply Theorem 5 for
$\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})$.

######  Lemma 11.

For every $\theta,\theta^{\prime}\in\Delta$ and $a\in(0,1)$,

$\lim_{n\rightarrow\infty}P_{\theta^{\ast}}\left(\frac{1}{n}\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})\leq
ac_{u}\|\theta-\theta^{\prime}\|^{2}\right)=0.$

###### Proof.

For any integer $k<n(1-ac_{u})<n$, we have

$\displaystyle
P_{\theta^{\ast}}\left(\frac{1}{n}\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})\leq
ac_{u}\|\theta-\theta^{\prime}\|^{2}\right)$
$\displaystyle=P_{\theta^{\ast}}\left(\|\theta-\theta^{\prime}\|^{2}\frac{1}{n}\sum_{t=1}^{n}(1-W_{t}(\theta,\theta^{\prime}))\leq
ac_{u}\|\theta-\theta^{\prime}\|^{2}\right)$
$\displaystyle=P_{\theta^{\ast}}\left(\frac{1}{n}\sum_{t=1}^{n}(1-W_{t}(\theta,\theta^{\prime}))\leq
ac_{u}\right)$
$\displaystyle=P_{\theta^{\ast}}\left(\sum_{t=1}^{n}W_{t}(\theta,\theta^{\prime})\geq
n\left(1-ac_{u}\right)\right)$
$\displaystyle=P_{\theta^{\ast}}\left(\sum_{t=1}^{n}W_{t}(\theta,\theta^{\prime})\geq
n\left(1-ac_{u}\right)\right)$
$\displaystyle\leq\frac{1}{\binom{n\left(1-ac_{u}\right)}{k}}\sum_{A\subseteq\\{1,\ldots,n\\}:|A|=k}E_{\theta^{\ast}}\left[\prod_{i\in
A}W_{i}(\theta,\theta^{\prime})\right]$
$\displaystyle\leq\frac{1}{\binom{n\left(1-ac_{u}\right)}{k}}\binom{n}{k}\left(1-c_{u}\right)^{k},$

where the last two lines follow from Theorem 5 and Lemma 10, respectively.
Select $k=k^{\ast}=\lceil n(1-a)\rceil$ and note that when $n>1/(a-ac_{u})$
one always has $k^{\ast}<n(1-ac_{u})$. Then, for $n>1/(a-ac_{u})$,

$\displaystyle
P_{\theta^{\ast}}\left(\frac{1}{n}\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime})\leq
ac_{u}\|\theta-\theta^{\prime}\|^{2}\right)$
$\displaystyle\leq\frac{1}{\binom{n\left(1-ac_{u}\right)}{k^{\ast}}}{\binom{n}{k^{\ast}}}\left(1-c_{u}\right)^{k^{\ast}}$
$\displaystyle=\left(1-c_{u}\right)^{k^{\ast}}\prod_{i=1}^{k^{\ast}}\frac{n-i+1}{n(1-ac_{u})-i+1}.$
(44)

The right-hand side does not depend on $\theta,\theta^{\prime}$ and converges
to $0$. ∎

##### Smoothness of $\ln h_{S,\epsilon}(y|\theta)$:

Next, we establish the $L$-smoothness of $h_{S,\epsilon}(y|\theta)$ as a
function of $\theta$ for any $y\in\\{1,\ldots,K\\}$ and
$S\subset\\{1,\ldots,K\\}$. Some technical lemmas are needed first.

######  Lemma 12.

For $x>0$, $\ln(1+x)\geq x-\frac{1}{2}x^{2}$.

###### Proof.

The function $\ln(1+x)-x+0.5x^{2}$ is $0$ at $x=0$ and its derivative
$1/(1+x)-1+x=\frac{1}{1+x}-(1-x)=\frac{1-(1-x^{2})}{1+x}=x^{2}/(1+x)>0$ when
$x>0$. ∎

######  Lemma 13.

For $x>0$, $\ln(x+1)\leq 1-\frac{1}{x+1}+\frac{1}{2}x^{2}$.

###### Proof.

The function $\ln(x+1)-1+\frac{1}{x+1}-\frac{1}{2}x^{2}$ is $0$ at $x=0$ and
its derivative,

$\frac{1}{x+1}-\frac{1}{(x+1)^{2}}-x=\frac{x-x(x+1)^{2}}{(x+1)^{2}}=\frac{x(1-(x+1)^{2})}{(x+1)^{2}},$

is negative for $x>0$. ∎

######  Lemma 14.

There exists an $L_{0}>0$ such that, for all $S\subset\\{1,\ldots,K\\}$ and
$y\in\\{1,\ldots,K\\}$, the function $\ln h_{S,\epsilon}(y|\theta)$ is an
$L_{0}$-smooth function of $\theta$. That is for all
$\theta,\theta^{\prime}\in\Delta$, we have

$\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})\leq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+L_{0}\|\theta-\theta^{\prime}\|^{2}.$

###### Proof.

Assume $h_{S,\epsilon}(y|\theta)\geq h_{S,\epsilon}(y|\theta^{\prime})$. Using
Lemma 13 with
$x=\frac{h_{S,\epsilon}(y|\theta)}{h_{S,\epsilon}(y|\theta^{\prime})}-1\geq
0$, we have

$\displaystyle\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})$ $\displaystyle\leq
1-\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}+\frac{1}{2}\left(\frac{h_{S,\epsilon}(y|\theta)}{h_{S,\epsilon}(y|\theta^{\prime})}-1\right)^{2}$
$\displaystyle=\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{1}{2}\left(\frac{h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta^{\prime})}\right)^{2}$
$\displaystyle\leq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{1}{2c^{2}}\left(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime})\right)^{2}$
(45) $\displaystyle\leq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{1}{2c^{2}}\|\theta-\theta^{\prime}\|^{2},$

where the first term in the second line follows from Lemma 3, the third line
follows from Lemma 2, and the last line follows from Lemma 8.

Now, assume $h_{S,\epsilon}(y|\theta)\leq h_{S,\epsilon}(y|\theta^{\prime})$.
By Lemma 12 with
$x=\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}-1\geq
0$,

$\ln h_{S,\epsilon}(y|\theta^{\prime})-\ln
h_{S,\epsilon}(y|\theta)\geq\left(\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}-1\right)-\frac{1}{2}\left(\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}-1\right)^{2},$

or, reversing the sign of inequality,

$\displaystyle\ln h_{S,\epsilon}(y|\theta)-\ln
h_{S,\epsilon}(y|\theta^{\prime})$
$\displaystyle\leq\left(1-\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}\right)+\frac{1}{2}\left(\frac{h_{S,\epsilon}(y|\theta^{\prime})}{h_{S,\epsilon}(y|\theta)}-1\right)^{2}$
$\displaystyle\leq\nabla_{\vartheta}\ln
h_{S,\epsilon}(y|\theta)^{\top}(\vartheta-\vartheta^{\prime})+\frac{1}{2c^{2}}\left(h_{S,\epsilon}(y|\theta)-h_{S,\epsilon}(y|\theta^{\prime})\right)^{2},$

where the third line follows from Lemma 2. Hence we have arrived at the same
inequality as (A.4.1). Hence, Lemma 14 holds with $L_{0}=\frac{1}{2c^{2}}$. ∎

##### Second moment of the gradient at $\theta^{\ast}$:

Let the average log-marginal likelihoods be defined as

$\Phi_{n}(\theta):=\frac{1}{n}\sum_{t=1}^{n}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta),\quad n\geq 1.$ (46)

The following bound on the second moment of this average at $\theta^{\ast}$
will be useful.

######  Lemma 15.

For $\Phi_{n}(\theta)$ defined in (46), we have

$\displaystyle
E_{\theta^{\ast}}\left[\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}\right]\leq\frac{1}{n}\max_{S\subset\\{1,\ldots,K\\}}\textup{Tr}\left[F(\theta^{\ast};S,\epsilon)\right].$

###### Proof.

First, we evaluate the mean at $\theta=\theta^{\ast}$.

$\displaystyle
E_{\theta^{\ast}}\left[\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\right]=\frac{1}{n}\sum_{t=1}^{n}E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\right].$

Focusing on a single term,

$\displaystyle E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\right]$
$\displaystyle=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y_{t}|\theta^{\ast})|S_{t}=S\right].$

Each term in the sum is equal to $0$, since

$E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y_{t}|\theta^{\ast})|S_{t}=S\right]=\sum_{k=1}^{K}\nabla_{\vartheta}\ln
h_{S,\epsilon}(k|\theta^{\ast})h_{S,\epsilon}(k|\theta^{\ast})=0.$ (47)

For the second moment at $\theta=\theta^{\ast}$,

$\displaystyle
E_{\theta^{\ast}}\left[\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}\right]=$
$\displaystyle\frac{1}{n^{2}}\sum_{t=1}^{n}E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})^{\top}\right]$
$\displaystyle+\frac{2}{n^{2}}\sum_{t=1}^{n}\sum_{t^{\prime}=1}^{t-1}E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t^{\prime}}}(Y_{t^{\prime}}|\theta^{\ast})^{\top}\right].$

For the diagonal terms, for all $t=1,\ldots,n$, we have

$\displaystyle E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})^{\top}\right]$
$\displaystyle\quad\quad=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y_{t}|\theta^{\ast})^{\top}|S_{t}=S\right]$
$\displaystyle\quad\quad=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)F(\theta^{\ast};S,\epsilon).$

For the cross terms, for $1\leq t^{\prime}<t\leq n$,

$\displaystyle E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t^{\prime}}}(Y_{t^{\prime}}|\theta^{\ast})^{\top}\right]$
$\displaystyle=E_{\theta^{\ast}}\left\\{E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t^{\prime}}}(Y_{t^{\prime}}|\theta^{\ast})^{\top}|Y_{t^{\prime}},S_{t^{\prime}}\right]\right\\}$
$\displaystyle=E_{\theta^{\ast}}\left\\{E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})|Y_{t^{\prime}},S_{t^{\prime}}\right]\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t^{\prime}}|\theta^{\ast})^{\top}\right\\}.$

The conditional expectation inside is zero, since, by (47),

$\displaystyle E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})|Y_{t^{\prime}},A_{t^{\prime}}\right]=\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S|Y_{t^{\prime}},S_{t^{\prime}})E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S,\epsilon}(Y_{t}|\theta^{\ast})|S_{t}=S\right]=0.$

Therefore, all the cross terms are zero,

$E_{\theta^{\ast}}\left[\nabla_{\vartheta}\ln
h_{S_{t},\epsilon}(Y_{t}|\theta^{\ast})\nabla_{\vartheta}\ln
h_{S_{t^{\prime}}}(Y_{t^{\prime}}|\theta^{\ast})^{\top}\right]=0,$

and hence we arrive at

$\displaystyle
E_{\theta^{\ast}}\left[\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}\right]=\frac{1}{n^{2}}\sum_{t=1}^{n}\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)\left[F(\theta^{\ast};S,\epsilon)\right]$
(48)

for the second moment of the gradient of $\Phi_{n}(\theta^{\ast})$. Therefore,

$\displaystyle
E_{\theta^{\ast}}\left[\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}\right]$
$\displaystyle=E_{\theta^{\ast}}\left[\textup{Tr}\left(\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}\right)\right]$
$\displaystyle=\textup{Tr}\left(E_{\theta^{\ast}}\left[\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}\right]\right)\
$
$\displaystyle=\frac{1}{n^{2}}\sum_{t=1}^{n}\sum_{S\subset\\{1,\ldots,K\\}}P_{\theta^{\ast}}(S_{t}=S)\textup{Tr}(F(\theta^{\ast};S,\epsilon))$
$\displaystyle\leq\frac{1}{n}\max_{S\subset\\{1,\ldots,K\\}}\textup{Tr}\left(F(\theta^{\ast};S,\epsilon)\right),$

which concludes the proof. ∎

####  A.4.2 Convergence of the posterior distribution

Let $\mu\in(0,1)$ and, for $\theta,\theta^{\prime}\in\Delta$, define

$\mathcal{E}^{\mu}_{n}(\theta,\theta^{\prime}):=\mu
c_{u}\|\theta-\theta^{\prime}\|^{2}-\frac{1}{n}\sum_{t=1}^{n}V_{t}(\theta,\theta^{\prime}),$
(49)

where $V_{t}(\theta,\theta^{\prime})$ was defined in (36) and $c_{u}>0$ was
defined in the proof of Lemma 9, respectively. The proof of Theorem 3 requires
the following lemma concerning
$\mathcal{E}^{\mu}_{n}(\theta,\theta^{\prime})$.

######  Lemma 16.

There exists $\mu\in(0,1)$ such that, for any $\varepsilon>0$, we have

$\lim_{n\rightarrow\infty}P_{\theta^{\ast}}\left(\int_{\Delta}e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\mathrm{d}\theta>e^{\varepsilon}\right)=0.$

###### Proof.

Define the product measure

$P^{\otimes}(\mathrm{d}(\theta,\cdot)):=\frac{\mathrm{d}\theta}{|\Delta|}\times
dP_{\theta^{\ast}}(\cdot)$ (50)

for random variables
$(\Theta\in\Delta,\\{S_{t}\subset\\{1,\ldots,K\\},Y_{t}\in\\{1,\ldots,K\\}\\}_{t\geq
1})$, where $\mathrm{d}\theta$ is the Lebesgue measure for $\vartheta$
restricted to $\Delta$ and $|\Delta|:=\int_{\Delta}\mathrm{d}\theta$. We will
show that the parameter $\mu$ in (49) can be chosen such that the collection
of random variables

$\mathcal{C}:=\\{f_{n}:=\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\Theta,\theta^{\ast})}\\}:n\geq
1\\}$

is uniformly integrable with respect to $P^{\otimes}$. For uniform
integrability, we need to show that for any $\varepsilon>0$, there exists a
$K>0$ such that

$E^{\otimes}[|f_{n}|\cdot\mathbb{I}(f_{n}>K)]<\varepsilon,\quad\forall n\geq
1.$

For any $K>1$ and $n\geq 1$, we have

$\displaystyle
E^{\otimes}[e^{nm_{0}\mathcal{E}^{\mu}_{n}(\Theta,\theta^{\ast})}\mathbb{I}(f_{n}>K)]$
$\displaystyle\leq\sup_{\theta\in
A_{n}}e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}P^{\otimes}(f_{n}>K)$
$\displaystyle=\sup_{\theta\in
A_{n}}e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\int_{\Delta}P_{\theta^{\ast}}\left(\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})>\frac{\ln
K}{nm_{0}}\right)\frac{\mathrm{d}\theta}{|\Delta|}$ $\displaystyle\leq e^{n\mu
m_{0}c_{u}}P_{\theta^{\ast}}\left(\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})>0\right),$

where the last line follows from
$\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\leq\mu
c_{u}\|\theta-\theta^{\ast}\|^{2}\leq\mu c_{u}$. Using (44), the last
expression can be upper-bounded as

$\displaystyle e^{nm_{0}\mu
c_{u}}P_{\theta^{\ast}}\left(\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})>0\right)$
$\displaystyle=e^{n\mu
m_{0}c_{u}}P_{\theta^{\ast}}\left(\frac{1}{n}\sum_{t=1}^{t}V_{t}(\theta,\theta^{\ast})<\mu
c_{u}\|\theta-\theta^{\prime}\|^{2}\right)$ $\displaystyle\leq e^{n\mu
m_{0}c_{u}}\left(1-c_{u}\right)^{\lceil n(1-\mu)\rceil}\prod_{i=1}^{\lceil
n(1-\mu)\rceil}\frac{n-i+1}{n(1-\mu c_{u})-i+1}.$ (51)

The parameter $\mu$ can be arranged such that (51) converges to $0$. For such
$\mu$, we have that for any $\varepsilon$ there exists a $N_{\varepsilon}>0$
such that for all $n>N_{\varepsilon}$,
$E^{\otimes}[|f_{n}|\cdot\mathbb{I}(f_{n}>K)]<\varepsilon$ for any $K>0$.
Finally, choose $K_{\varepsilon}=\sup_{\theta\in
A_{N_{\varepsilon}}}e^{N_{\varepsilon}m_{0}\mu c_{u}}$ so that
$E^{\otimes}[|f_{n}|\cdot\mathbb{I}(f_{n}>K_{\varepsilon})]<\varepsilon$ for
any $n\geq 1$. Hence, $\mathcal{C}$ is uniformly integrable for a suitable
choice of $\mu$.

Next, we show that each $f_{n}$ in $\mathcal{C}$ converges in probability to
$1$. The convergence is implied by the fact that for every $\theta\in A_{n}$
and $\varepsilon>0$, we have

$P_{\theta^{\ast}}(\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\\}>e^{\varepsilon})=P_{\theta^{\ast}}(\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})>\varepsilon)\rightarrow
0.$

by Lemma 11. Since $\mathcal{C}$ is uniformly integrable, the Vitali
convergence theorem ensures that $f_{n}$ converges in distribution (with
respect to $P^{\otimes}$) to $1$, i.e.,
$\lim_{n\rightarrow\infty}E^{\otimes}(f_{n})=1$. Since
$P^{\otimes}=\frac{\mathrm{d}\theta}{|\Delta|}\times
dP_{\theta^{\ast}}(\cdot)$ is a product measure as defined in (50), the stated
limit implies that

$E_{\theta^{\ast}}\left[\int_{\Delta}\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\\}\mathrm{d}\theta\right]\rightarrow
1,$

that is, the sequence
$\int_{A_{n}}\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\\}\mathrm{d}\theta$
converges to $1$ in distribution with respect to $P_{\theta^{\ast}}$. Since
convergence in distribution to a constant implies convergence in probability,
we have

$\int_{\Delta}\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\\}\mathrm{d}\theta\overset{P_{\theta^{\ast}}}{\rightarrow}1.$
(52)

Finally, since we have

$\int_{\Delta}e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\mathrm{d}\theta\leq\int_{\Delta}\max\\{1,e^{nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})}\\}\mathrm{d}\theta,$

and the right-hand side converges in probability to 1, we conclude. ∎

###### Proof of Theorem 3.

Writing down Lemma 7 with $\theta^{\ast}$ and any $\theta\in\Delta$ separately
for $t=1,\ldots,n$, summing the inequalities and dividing by $n$, we obtain

$\displaystyle\Phi_{n}(\theta^{\ast})-\Phi_{n}(\theta)$
$\displaystyle\geq\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}(\vartheta^{\ast}-\vartheta)+m_{0}\sum_{t=1}^{n}V_{t}(\theta,\theta^{\ast})$
$\displaystyle=\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}(\vartheta^{\ast}-\vartheta)+m\|\theta-\theta^{\ast}\|^{2}-m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$

where $m:=\mu m_{0}c_{u}$. Reversing the sign,

$\displaystyle\Phi_{n}(\theta)-\Phi_{n}(\theta^{\ast})\leq\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}(\vartheta-\vartheta^{\ast})-m\|\theta^{\ast}-\theta\|^{2}+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$

Using Cauchy-Schwarz inequality for the first term on the right-hand side, we
get

$\displaystyle\Phi_{n}(\theta)-\Phi_{n}(\theta^{\ast})\leq\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|\|\theta-\theta^{\ast}\|-m\|\theta^{\ast}-\theta\|^{2}+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$

Using Young’s inequality $uv\leq\frac{u^{2}}{2\kappa}+\frac{v^{2}\kappa}{2}$
for the second term with $u=\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|$,
$v=\|\theta^{\ast}-\theta\|$, and $\kappa=m$, we get

$\Phi_{n}(\theta)-\Phi_{n}(\theta^{\ast})\leq\frac{\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}}{2m}-\frac{m}{2}\|\theta^{\ast}-\theta\|^{2}+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$
(53)

Similarly, using Lemma 14 with $\theta^{\ast}$ and any
$\theta^{\prime}\in\Delta$ for $t=1,\ldots,n$, summing the inequalities and
dividing by $n$, we obtain

$\Phi_{n}(\theta^{\ast})-\Phi_{n}(\theta^{\prime})\leq\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})^{\top}(\vartheta^{\ast}-\vartheta^{\prime})+L_{0}\|\theta^{\ast}-\theta^{\prime}\|^{2}.$

Again, using Cauchy-Schwarz inequality and Young’s inequality
$uv\leq\frac{u^{2}}{2\kappa}+\frac{v^{2}\kappa}{2}$ with
$u=\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|$,
$v=\|\theta^{\ast}-\theta\|$, and $\kappa=2L_{0}$, we get

$\displaystyle\Phi_{n}(\theta^{\ast})-\Phi_{n}(\theta^{\prime})$
$\displaystyle\leq\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|\|\theta^{\ast}-\theta^{\prime}\|+L_{0}\|\theta^{\ast}-\theta^{\prime}\|^{2}$
$\displaystyle\leq\frac{\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}}{4L_{0}}+2L_{0}\|\theta^{\ast}-\theta^{\prime}\|^{2}$
$\displaystyle=\frac{\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}}{2L}+L\|\theta^{\ast}-\theta^{\prime}\|^{2},$
(54)

where we let $L:=2L_{0}$. Summing the inequalities in (53) and (54), we obtain

$\Phi_{n}(\theta)-\Phi_{n}(\theta^{\prime})\leq\left(\frac{1}{2L}+\frac{1}{2m}\right)\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-\frac{m}{2}\|\theta^{\ast}-\theta\|^{2}+\frac{L}{2}\|\theta^{\ast}-\theta^{\prime}\|^{2}+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$
(55)

Let $a\in(0,1)$ be a constant and define the sequences

$\displaystyle\Omega_{n}$
$\displaystyle:=\\{\theta\in\Delta:\|\theta-\theta^{\ast}\|^{2}<\max\\{4/m,4/L\\}n^{-a}\\},\quad
n\geq 1.$ $\displaystyle A_{n}$
$\displaystyle:=\\{\theta\in\Delta:\|\theta-\theta^{\ast}\|^{2}>\max\\{4/m,4/L\\}n^{-a}\\},\quad
n\geq 1.$ $\displaystyle B_{n}$
$\displaystyle:=\\{\theta\in\Delta:\|\theta-\theta^{\ast}\|^{2}\leq\min\\{2/m,2/L\\}n^{-a}\\},\quad
n\geq 1.$

For $\theta\in A_{n}$ and $\theta^{\prime}\in B_{n}$, (55) can be used to
obtain

$\Phi_{n}(\theta)-\Phi_{n}(\theta^{\prime})\leq\left(\frac{1}{2L}+\frac{1}{2m}\right)\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{-a}\left(\max\left\\{2,\frac{2m}{L}\right\\}-\min\left\\{\frac{L}{m},1\right\\}\right)+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}).$

Noting that
$\max\left\\{2,\frac{2m}{L}\right\\}-\min\left\\{\frac{L}{m},1\right\\}\geq
1$, we have

$\Phi_{n}(\theta)-\Phi_{n}(\theta^{\prime})\leq\left(\frac{1}{2L}+\frac{1}{2m}\right)\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{-a}+m_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast}),\quad\theta\in
A_{n};\theta^{\prime}\in B_{n}.$ (56)

Multiplying (56) with $n$, exponentiating, and multiplying the ratio of the
priors, we get

$\displaystyle\frac{\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}}{\eta(\theta^{\prime})\exp\\{n\Phi_{n}(\theta^{\prime})\\}}$
$\displaystyle\leq\frac{\eta(\theta)}{\eta(\theta^{\prime})}\exp\left[\left(\frac{1}{2L}+\frac{1}{2m}\right)n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{1-a}+nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]$
$\displaystyle\leq
C_{\eta,n}\exp\left[C_{1}n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{1-a}+nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]$
(57)

for all $\theta\in A_{n}$ and $\theta^{\prime}\in B_{n}$, where
$C_{1}:=\frac{1}{2L}+\frac{1}{2m}$ and $C_{\eta,n}:=\sup_{\theta\in
A_{n},\theta^{\prime}\in B_{n}}\frac{\eta(\theta)}{\eta(\theta^{\prime})}$.
The bound in (57) can be used to bound the ratio between the posterior
probabilities $\Pi(A_{n}|Y_{1:n},S_{1:n})$ and $\Pi(B_{n}|Y_{1:n},S_{1:n})$,
since

$\displaystyle\frac{\Pi(A_{n}|Y_{1:n},S_{1:n})}{\Pi(B_{n}|Y_{1:n},S_{1:n})}$
$\displaystyle=\frac{\int_{A_{n}}\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}\mathrm{d}\theta}{\int_{B_{n}}\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}\mathrm{d}\theta}$
$\displaystyle=\frac{\int_{A_{n}}\frac{\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}}{\inf_{\theta^{\prime}\in
B_{n}}\eta(\theta^{\prime})\exp\\{n\Phi_{n}(\theta^{\prime})\\}}\mathrm{d}\theta}{\int_{B_{n}}\frac{\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}}{\inf_{\theta^{\prime}\in
B_{n}}\eta(\theta^{\prime})\exp\\{n\Phi_{n}(\theta^{\prime})\\}}\mathrm{d}\theta}$
$\displaystyle\leq\frac{\int_{A_{n}}\frac{\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}}{\inf_{\theta^{\prime}\in
B_{n}}\eta(\theta^{\prime})\exp\\{n\Phi_{n}(\theta^{\prime})\\}}\mathrm{d}\theta}{\int_{B_{n}}\frac{\eta(\theta)\exp\\{n\Phi_{n}(\theta)\\}}{\inf_{\theta^{\prime}\in
B_{n}}\eta(\theta^{\prime})\exp\\{n\Phi_{n}(\theta^{\prime})\\}}\mathrm{d}\theta}$
$\displaystyle\leq\frac{\int_{A_{n}}C_{\eta,n}\exp\left[C_{1}n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{1-a}+nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]\mathrm{d}\theta}{\int_{B_{n}}1\mathrm{d}\theta}$
$\displaystyle=\frac{1}{\text{Vol}(B_{n})}C_{\eta,n}\exp\left[C_{1}n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}-n^{1-a}\right]\int_{A_{n}}\exp\left[nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]\mathrm{d}\theta,$

where $\text{Vol}(B_{n}):=\int_{B_{n}\cap\Delta}\mathrm{d}\theta$. Note that
$B_{n}$ shrinks with $n$, so there exists a $N_{B}$ such that for $n>N_{B}$,
the volume $B_{n}$ can be lower-bounded as

$B_{n}\geq\frac{1}{2(K-2)!}\frac{(\sqrt{\pi}\min\\{2/m,2/L\\}n^{-a})^{(K-1)}}{\Gamma((K-1)/2+1)},$

where the factor $\frac{1}{2(K-2)!}$ corresponds to the worst-case situation
where $\theta^{\ast}$ is on one of the corners of $\Delta$, such as
$\theta^{\ast}=(1,0,\ldots,0)^{\top}$, and the rest is the volume of a $K-1$
dimensional sphere with radius $\min\\{2/m,2/L\\}n^{-\alpha}$. The lower bound
is the volume of the intersection of a simplex with a sphere centered at one
of the sharpest corners of the simplex. Therefore, for $n>N_{B}$, ratio can
further be bounded as

$\frac{\Pi(A_{n}|Y_{1:n},S_{1:n})}{\Pi(B_{n}|Y_{1:n},S_{1:n})}\leq
C_{2}C_{\eta,n}\exp\left[C_{1}n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}+(K-1)a\ln
n-n^{1-a}\right]\int_{A_{n}}\exp\left[nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]\mathrm{d}\theta,$

where $C_{2}:=\frac{\Gamma((K-1)/2+1)}{\min\\{2/m,2/L\\}^{K-1}\pi^{(K-1)/2}}$
does not depend on $n$.

Next, we prove that the sequence of random variables

$Z_{n}:=C_{\eta,n}\exp\left[C_{1}n\|\nabla_{\vartheta}\Phi_{n}(\theta^{\ast})\|^{2}+(K-1)a\ln
n-n^{1-a}\right]\int_{A_{n}}\exp\left[nm_{0}\mathcal{E}^{\mu}_{n}(\theta,\theta^{\ast})\right]\mathrm{d}\theta$

converges to $0$ in probability, which in turn proves the convergence of
$\frac{\Pi(A_{n}|Y_{1:n},S_{1:n})}{\Pi(B_{n}|Y_{1:n},S_{1:n})}$ in probability
to $0$. To do that, we need to prove that for each $\varepsilon>0$ and
$\delta>0$ there exists a $N>0$ such that for all $n>N$ we have