Frequentist Guarantees of Distributed (Non)-Bayesian Inference

Bohan Wu<EMAIL_ADDRESS>
Department of Statistics

Columbia University

New York, NY 10027, USA

César A. Uribe<EMAIL_ADDRESS>
Department of Electrical and Computer Engineering

Rice University

Houston, TX 77005, USA

We establish Frequentist properties, i.e., posterior consistency, asymptotic normality, and posterior contraction rates, for the distributed (non-)Bayes Inference problem for a set of agents connected over a network. These results are motivated by the need to analyze large, decentralized datasets, where distributed (non)-Bayesian inference has become a critical research area across multiple fields, including statistics, machine learning, and economics. Our results show that, under appropriate assumptions on the communication graph, distributed (non)-Bayesian inference retains parametric efficiency while enhancing robustness in uncertainty quantification.  We also explore the trade-off between statistical efficiency and communication efficiency by examining how the design and size of the communication graph impact the posterior contraction rate. Furthermore, we extend our analysis to time-varying graphs and apply our results to exponential family models, distributed logistic regression, and decentralized detection models.

Distributed inference, Bayesian theory, Berstein-von Mises, communication efficiency, estimation over networks

§ INTRODUCTION

Modern datasets are frequently generated and stored by distributed systems, including social media, sensor networks, blockchain, and cloud-based databases. However, transmission costs make analyzing these datasets on a centralized machine prohibitively expensive and, in some cases, infeasible. To address this challenge, researchers have turned to distributed algorithms that enable decentralized data-driven decision-making under communication constraints [Borkar and Varaiya, 1982, Tsitsiklis and Athans, 1984, Gubner, 1993]. In such systems, a set of agents operates within a communication network structure, where each agent can only share information locally with its neighbors. The agents sequentially analyze the data, with each agent performing inference independently and sharing the results through edges defined by the network structure, which may vary over time [Nedić et al., 2017, Uribe et al., 2022].

Decentralized or distributed Bayesian inference originates in statistics [DeGroot, 1974, Gilardoni and Clayton, 1993]. However, it wasn't until the massive advances in computing power in the past decade that the ideas of distributed inference started regaining interest in the statistics community [Uribe et al., 2022, Uribe et al., 2022]. There is a growing line of works on Distributed Bayesian inference, which aims to develop scalable and efficient algorithms for posterior computation on large datasets [Jordan et al., 2018]. One of the main challenges in this area is to design data-parallel procedures that can handle massive datasets by breaking them into smaller blocks that can be processed independently on individual machines. Much of the current literature focuses on “one-shot” or “embarrassingly parallel” approaches, which involve only one round of communication between local machines and a central node. In these approaches, estimators or posterior samples are obtained in parallel on local machines, communicated to a central node, and combined to form a global estimator or approximation to the posterior distribution.

From the Markov chain Monte Carlo (MCMC) perspective, there have been several developments in parallel MCMC methods for distributed Bayesian inference [Neiswanger et al., 2013, Wang and Dunson, 2013, Minsker et al., 2014, Wang et al., 2015, Rabinovich et al., 2015, Scott et al., 2016, Li et al., 2017, Minsker et al., 2017]. These methods draw samples from the subset posterior in parallel agents and combine the samples to obtain an approximation to the posterior measure for the complete data.

From the variational Bayes perspective, algorithms such as stochastic variational inference [Hoffman et al., 2013] have been proposed for distributed Bayesian inference. These algorithms distribute the data across machines, implement the local variational updates in parallel through stochastic gradient descent (SGD), and update the global variational parameters as a weighted average of local optimums. The variational interpretation of the Bayes rule [Walker, 2006] allows the representation between the variational optimization problem and posterior to go both ways.

In parallel to the success in the statistics community, microeconomics has also seen a surge of interest in Distributed Bayesian inference, often referred to as “non-Bayesian social learning.” Notable works in this area have focused on its axiomatic foundations [Epstein et al., 2010], conditions for achieving consensus [Acemoglu et al., 2011], various learning rules [Golub and Sadler, 2017], and the effects of information aggregation [Molavi et al., 2018]. This cross-disciplinary interest underscores the broad relevance and applicability of Distributed Bayesian inference techniques. Here, the agents represent individuals seeking to learn about an underlying state of the world $\theta$. Unlike traditional Bayesian approaches, where each agent makes inferences based on the full data, the “non-Bayesian learning” (as economists call it) model captures how individuals make inferences in the presence of other decision-makers, often with limited access to information and interaction with social networks. The Bayesian distributed learning framework offers a promising solution that retains the desirable properties of Bayesian learning, such as ease of uncertainty quantification and flexibility while incorporating a form of information aggregation that aligns with realistic behavioral assumptions. Indeed, the distributed Bayes rule has been shown to reflect reasonable assumptions about individual behavior in society [Jadbabaie et al., 2012]. The distributed framework is analytically tractable under certain distributional assumptions and computationally feasible in general.  For a more comprehensive literature review, see
[Molavi et al., 2018].

Studies have also investigated distributed Bayesian frameworks for building large-scale decentralized machine-learning systems. These algorithms generally involve an aggregation step consisting of a weighted geometric or arithmetic average of the received beliefs, followed by a Bayesian update using locally available data. In this setup, the set of hypotheses is assumed to be common across the network, but individual observations may come from different distributions. Recent works, including [Lalitha et al., 2014, Qipeng et al., 2015, Shahrampour et al., 2015, Rahimian et al., 2015], proposed variations of this approach and established consistent, geometric, and non-asymptotic convergence rates for a general class of distributed algorithms, covering asymptotic analysis, non-asymptotic bounds, time-varying directed graphs, and adversarial agents, transmission, and node failures.

Distributed Bayesian procedures have attracted substantial interest across disciplines such as electrical engineering, statistics, and economics. Yet, the broad adoption of these methods in the statistical community has been hindered by a lack of rigorous analysis of their statistical properties. Moreover, understanding these properties is key to deepening our knowledge of the consensus behavior of agents within varying communication patterns, a topic of interest to the electrical engineering and economics communities. In this paper, we investigate the distributed Bayes procedures that arise from applying the stochastic mirror descent (SMD) algorithm to statistical estimation problems [Uribe et al., 2022, Uribe et al., 2022]. Our work fills a crucial gap by establishing the Frequentist properties of such distributed Bayes procedures, including posterior consistency, asymptotic normality, and posterior contraction rates. We also explore the tradeoff between statistical efficiency and communication cost by investigating the relationship between the posterior contraction rate and the structure of the communication graph. Furthermore, we illustrate practical applications of the Bernstein von - Mises results to emphasize their utility in uncertainty quantification. Ultimately, we hope to stimulate further interest in distributed Bayes methods within the statistical community and to establish a solid theoretical foundation for distributed Bayesian inference in fields such as economics and electrical engineering.

Table of Notation

$D_{KL}$     $\triangleq$     Kullback–Leibler divergence,

$L$     $\triangleq$    loss function,

$\langle.,.\rangle$     $\triangleq$     $L_2$ inner product, as in $\langle f,g \rangle = \int_\R f(x)g(x) dx$.

3cProbability Distributions:

$\bP, \bp$     $\triangleq$     data generating measure, data generating density,

$\bP f(X)$     $\triangleq$    expectation of $f(X)$ when $X \sim \bP$, same as $\E_\bP f(X)$,

$\Pi, \pi$     $\triangleq$    prior measure, prior density.

$P, p$     $\triangleq$    posterior measure, posterior density,

$P_t^j, p_t^j$     $\triangleq$     posterior measure, density for agent $j$ at iterant (time) $t$,

$G,A, A_{ij}$     $\triangleq$     Graph, adjacency matrix and its $(i,j)^{th}$ entry,

$\lambda_j(A)$     $\triangleq$     $j^{th}$ largest eigenvalue of $A$,

$\ones$     $\triangleq$     the vector of all ones.

The rest of the paper is outlined as follows: Section <ref> introduces the distributed Bayesian inference problem from an optimization-centric viewpoint and rigorously defines the distributed Bayes posterior. In Section <ref>, we outline sufficient conditions for the consistency and asymptotic normality of the distributed Bayes posterior. Section <ref> establishes the consistency of distributed Bayes posterior (Theorem <ref>). Section <ref> establishes the Bernstein-von Mises theorems under both correct (Theorem <ref>) and incorrect model specifications(Theorem <ref>). Section <ref> provides both the abstract and concrete upper bounds on the posterior contraction rate (Theorem <ref> and Theorem <ref>), with an emphasis on model misspecification. Our analysis is extended to time-varying graphs in Section <ref>, where we establish posterior contraction rates under various communication frequency regimes (Theorem <ref>). In Section <ref>, we demonstrate the practical use of our findings by establishing Bernstein-von Mises results for three statistical models, including exponential family (Proposition <ref>), logistic regression models(Proposition <ref>), and the distributed detection problem (Proposition <ref>)- a canonical problem in electrical engineering. The paper concludes in Section <ref> with a discussion on future research directions following our findings.

§ PRELIMINARIES

Suppose we observe a sequence of i.i.d random variables $X_1, X_2, \cdots$ all taking values in a probability space $(\mathcal{X}, \mathbb{B}, \bP_0)$ where the true distribution $\bP_0$ is unknown. Moreover, let a family of probability distribution be given in the form of $\{\bP_\theta: \theta \in \Theta\}$ where $(\Theta, \mathcal{A})$ is a measurable space. Each $\bP_\theta$ is a probability measure defined on $(\mathcal{X}, \mathbb{B})$ and the mapping $\theta \to \bP_\theta(B)$ is measurable for every $B \in \mathbb{B}$. We refer to $\{\bP_\theta: \theta \in \Theta\}$ as the statistical model. The parameter space $\Theta$ is typically taken as a subset of the Euclidean or Hilbert space to avoid any measurability issue with $\{\bP_\theta: \theta \in \Theta\}$ [Ghosal and Van der Vaart, 2017].

The centralized statistical estimation problem is to fund a subset $\Theta_0 \subseteq \Theta$ such that for $\theta_0 \in \Theta_0$, $\bP_{\theta_0}$ is the “closest” to $\bP_0$ with respect to a metric.  Geometrically, the goal is to find a point in the subset of the probability measures $\{\bP_\theta: \theta \in \Theta\}$ closest to $\bP_0$ under a given topology. The topology is often defined by divergence on the space of probabilities, such as the Kullback-Leibler (KL) divergence, Rényi divergence, etc. The definition of KL, Rényi divergence, and other divergence functions are reviewed in Section <ref> of the Appendix.

This article uses the KL divergence to define the objective function for statistical estimation problems. This is a natural choice because the KL minimization problem is equivalent to the maximum likelihood estimation in statistics. In short, consider the following estimation problem,

\begin{equation}
\label{eqn: simple-KL}
\min_{\theta \in \Theta} D_{KL}(\bP_0\parallel \bP_\theta).
\end{equation}

Recall the definition of KL divergence,

\begin{equation*}
D_{KL}(\bP_0\parallel \bP_\theta) = - \bP_0 \log \frac{d \bP_\theta}{d \bP_0}.
\end{equation*}

where $\frac{d\bp_\theta(x)}{d\bp_0(x)}$ is the Radon-Nikodyn derivative of $\bP_\theta$ with respect to $\bP_0$.

There are a couple of challenges to solving (<ref>). First, $\bP_0$ is unknown, and we want to estimate $\bP_0$ using the realizations of $X$. For example, suppose the true distribution $\bP_0 = N(0,1)$, and $\bP_\theta = \{N(\theta, 1);  \theta \in [-1,1]\}$ where $N(\mu, \sigma^2)$ is the normal distribution with mean $\mu$ and $\sigma^2$.  Then, $\theta_0$ is estimated to be $0$ under the KL divergence. Secondly, first-order stochastic methods such as the gradient descent algorithm are the default approach to solving this optimization problem. But the set $\Theta$ might be discrete or non-smooth.

Instead of minimizing over points $\Theta$, we “smooth” the problem by minimizing over the probability distributions on $\Theta$ as described next.

Let $\Delta_\Theta$ be the space of probability density functions defined on $\Theta$. It follows that

\begin{equation}\label{eqn: simple-KL-Bayes}
\min_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta) = \min_{p  \in \Delta_\Theta} \int_\Theta p(\theta) D_{KL}(\bP_0\parallel \bP_\theta) d \theta,
\end{equation}

where $\Delta_\Theta$ is the (hypothetical) space of all probability distributions over $\Theta$, and , with $p(\theta) \geq 0$. The problem transforms from choosing the appropriate parameter to choosing the appropriate prior. If $\theta_0$ solves (<ref>), $p^* = \delta_{\theta_0}$ solves (<ref>), therefore the two problems are equivalent. The advantage is that formulation (<ref>) is a continuous optimization problem regardless of the topology in which $\Theta$ is defined, i.e., $\Theta$ might be a finite discrete set.

There is an equivalent formulation of Problem (<ref>) as a minimization problem over the space of probability measures on $\Theta$. Given the linearity of the expectation functional, the Riez representation theorem implies the existence of an inner product $\langle., .\rangle$ that characterizes the expectation over $\Theta$. We then reformulate the KL-minimization problem in terms of linear stochastic optimization as follows:

\begin{equation}
\begin{aligned}
\min_{\theta \in \Theta} D_{KL}(\bP_0\parallel \bP_\theta)  &= \min_{p \in \Delta_\Theta}  \E_{\bP_0}\langle p,- \log \bp_\theta(x) \rangle.
\end{aligned}
\end{equation}

Note that

\begin{align*}
\min_{\theta \in \Theta} D_{KL}(\bP_0\parallel \bP_\theta) &= \min_{p \in \Delta_\Theta} \int_\Theta p(\theta) D_{KL}(\bP_0\parallel \bP_\theta) d \theta  \\
&= \min_{p \in \Delta_\Theta} \int_\Theta p(\theta) \int_{\mathcal X} \bp_0(x) \log \frac{\bp_0(x)}{\bp_\theta(x)} dx d \theta \\
& = \min_{p \in \Delta_\Theta} \int_{\mathcal X} \bp_0(x) \int_\Theta p(\theta) \log \frac{\bp_0(x)}{\bp_\theta(x)} d \theta dx, \qquad \text{by Fubini's theorem} \\
&= \min_{p \in \Delta_\Theta} \int_{\mathcal X} \bp_0(x) \int_\Theta p(\theta) [\log \bp_0(x) - \log \bp_\theta(x) ] d \theta dx \\
&= \min_{p \in \Delta_\Theta} - \int_{\mathcal X} \bp_0(x) \int_\Theta p(\theta) \log \bp_\theta(x)  d \theta dx.
\end{align*}

In this reformulation, $\Delta_\Theta$ is the probability simplex over $\Theta$, and $\langle p,- \bp_\theta(x) \rangle$ represents the inner product between the simplex vector $p$ and the negative log-likelihood $\bp_\theta(x)$, taken with respect to the expectation under $\bP_0$.

This representation enables efficient computational methods from the stochastic optimization literature, and a proven solution is the stochastic mirror descent (SMD) algorithm [Uribe et al., 2022, Uribe et al., 2022]. Algorithm 6 in  [Uribe et al., 2022, Uribe et al., 2022] iterates through a KL-regularized optimization problem

\begin{equation}
\label{eqn:smd}
p_{t+1} = \argmin_{p \in \Delta_\Theta} \left\{-\langle p,  \log \bP_\theta(x_{t+1})\rangle + \alpha D_{KL}(p \parallel p_t) \right\}.
\end{equation}

At time $t$, the agent's goal is to maximize the expected log-likelihood given the belief $p_t$, i.e., minimizing the divergence. The trade-off between the two behaviors is tuned with the learning rate $\alpha$.

Equation (<ref>) is simultaneously the variational representation of a generalized Bayesian inference. We unpack this as follows:

\begin{align*}
p_{t+1} &= \argmin_{p \in \Delta_\Theta} \left\{-\langle p,  \log \bP_\theta(x_{t+1})\rangle + \alpha D_{KL}(p \parallel p_t) \right\}\\
&= \argmin_{p \in \Delta_\Theta}\left\{ - \E_p \log \bP_\theta(x_{t+1}) + \alpha \E_p \log p(\theta) -  \alpha \E_p \log p_t(\theta) \right\} \\
&= \argmin_{p \in \Delta_\Theta}\left\{ - \frac{1}{\alpha}\E_p \log \bP_\theta(X_{t+1}) + \E_p \log p(\theta) -  \E_p \log p_t(\theta) \right\}\\
&= \argmin_{p \in \Delta_\Theta} \left\{\E_p \log \frac{p(\theta)}{(\bP_\theta(x_{t+1}))^{\frac{1}{\alpha}} p_t(\theta)} \right\}\\
&=  \argmin_{p \in \Delta_\Theta} \left\{ D_{KL}(p \parallel \tilde p) \right\}, \quad \tilde p \propto (\bP_\theta(x_{t+1}))^{\frac{1}{\alpha}} p_t(\theta).
\end{align*}

A generalized posterior is a probability distribution $p(\theta \mid X_1, \ldots, X_t) \propto \pi(\theta) \tilde L_t(\theta)$, where $\pi$ is a prior and $\tilde L_t$ is a surrogate likelihood [Miller, 2021]. From this lens, the solution to stochastic mirror descent is a generalized Bayes posterior:

\begin{equation}
\label{eqn:simple-bayes}
p_{t+1}(\theta) \propto (\bP_\theta(x_{t+1}))^{\frac{1}{\alpha}} p_t(\theta).
\end{equation}

The parameter $\alpha$ is simultaneously the step size of the SMD update and a temperature parameter in the generalized Bayes' rule. When $\alpha$ is set to $1$, we recover the standard Bayes' theorem.

The result above suggests a fascinating interplay between stochastic optimization and Bayesian inference. While this may seem somewhat unexpected, it aligns well with existing literature that characterizes Bayesian rules and their generalizations as equivalent optimization problems [Knoblauch et al., 2022]. The connection between Bayesian methods and Stochastic Mirror Descent (SMD) is particularly fruitful to explore, as it enables us to formulate the distributed Bayesian rule due to the distributed SMD algorithm.

§.§ Distributed Bayesian Inference

Now, we present the distributed Bayesian inference problem.
Consider a set of $m$ independent agents, represented by the index $j = 1, ..., m$. Each agent independently observes a sequence of random variables at discrete time steps $t = 0,1, 2, ...$.  The random variable observed by agent $j$ at time $t$ is denoted as $X_t^j$. The random process observed by agent $j$ is denoted as $X^j$. The collection of random variables all agents observe at time $t$ is denoted as $X_t$. In general, each $X^j$ may be endowed with a different probability space with a common parameter space. For simplicity, we assume that all random variables are defined on the measure space $(\mathcal{X}, \mathbb{B})$, and follow a true probability distribution $\P$ with density $p$.

Given a  common set of parameters $\Theta$, the private parametric statistical model that agent $j$ can access is defined as $\mathcal{P}_\Theta^j = \{\bP_\theta^j: \theta \in \Theta\}$. Each model $\mathcal{P}_\Theta^j$ has the same support for all $j \in [m]$.

At time $t$, the agents interact through an undirected communication graph $G_t = (V, E_t)$, where $V = [m]$ denotes the set of agents. An edge $(j, i) \in E_t$ implies that agent $j$ can communicate with agent $i$ at time $t$.

The weighted adjacency matrix associated with $G_t$ is denoted by $A_t$.  We assume that $A_t$ is a doubly stochastic matrix obtained by normalizing the matrix $A_t'$, where $[A_t']_{ij} = 1$ if there exists a communication link between agent $i$ and agent $j$ and $[A_t']_{ij} = 0$ otherwise. For undirected graphs, one can always construct a doubly stochastic matrix by normalizing the adjacency matrix, so $A_t$ is guaranteed to exist.

In our framework, agents share information by communicating their beliefs, represented by $p_t^j(\theta)$. This notation represents the posterior distribution of $\theta$ as perceived by agent $j$ at time $t$.

Analogous to the centralized case, we formulate distributed statistical inference as a distributed optimization problem:

\begin{equation}
\label{eqn: dKL}
\min_{\theta \in \Theta} \sum_{j = 1}^m D_{KL}(\bP_0\parallel \bP_\theta^j).
\end{equation}

For each agent $j \in [m]$, let $\pi^j$ represent their initial belief or prior. Therefore, we have $p_0^j = \pi^j$. Similar to the statistical estimation for a single agent, the distributed statistical inference problem admits a reformulation via the distributed stochastic mirror descent algorithm [Uribe et al., 2022, Uribe et al., 2022].  The belief of the agent $j$ at time $t + 1$ is obtained through the following mirror descent update.

\begin{equation}
\label{eqn:dSMD}
p_{t+1}^j = \argmin_{p \in \Delta_\Theta} \left\{  -\langle \log \bp_\theta^j(x_{t+1}^j), p \rangle + \sum_{i  = 1}^m [A_t]_{ij} D_{KL}(p \parallel p_t^i) \right\},
\end{equation}

which is solved via the update

\begin{equation}
\label{eqn:dBr}
p_{t+1}^j(\theta) \propto \bp_\theta^j(x_{t+1}^j)\prod_{i = 1}^m (p_t^i(\theta))^{[A_t]_{ij}}.
\end{equation}

The derivation from (<ref>) to (<ref>) is identical to that of the case of a single agent (<ref>)–(<ref>).

Equation (<ref>) is aptly referred to as the distributed Bayes rule, as it generalizes the Bayesian rule to a distributed setting. This distributed Bayes rule induces a distributed Bayes posterior $p_{t+1}^j$, which might be viewed as the belief of each agent after updating its prior belief based on the new data and information received from its neighbors.

It's important to clarify that despite its name, the distributed Bayes posterior is not a posterior distribution in the conventional sense. The classic Bayes rule is a special case of this general distributed Bayes rule, where the communication graph has no edges, denoted as $A_t = I_m$. In this scenario, each agent operates independently, with no influence from others

\begin{equation*}
\begin{aligned}
p_{t + 1}^j(\theta)
&\propto  \bp_\theta(x_{t + 1}^j)p_t^j(\theta).
\end{aligned}
\end{equation*}

In the special case where the communication graph $G_t$ is fully connected, each agent can communicate with every other agent. This leads to the distributed Bayes posterior effectively acting as a weighted Bayesian update rule, with equal weights assigned to all agents.

The distributed Bayes posterior becomes a product of the likelihood of the incoming data and a “tempered” posterior with power $1/m$

\begin{equation*}
\begin{aligned}
p_{t + 1}^j(\theta)
&\propto  \bp_\theta(x_{t+1}^j) \prod_{i = 1}^m \prod_{k = 1}^t (\bp_\theta(x_k^i))^{1/m} \prod_{i = 1}^m \pi^i(\theta)^{1/m}.
\end{aligned}
\end{equation*}

These examples show that distributed Bayesian posterior can adapt to the underlying communication structure, balancing individualistic updates (when there are no communications) and collective updates (when the communication happens according to the network).

Analyzing $p_{t+1}^j(\theta)$ from the perspective of generalized posterior distributions is useful. In the next Section, we provide the measure-theoretic definition of the distributed Bayes posterior as a generalized Bayes posterior.

§.§ Distributed Bayes Posterior

This section provides a rigorous definition of the distributed Bayes posterior we study throughout this paper.

Define the prior measure $\Pi$ as follows:

\begin{align*}
\Pi(B) = \frac{\int_B \prod_{j = 1}^m \pi^j( \theta)^{1/m} d\theta}{\int_\Theta \prod_{j = 1}^m \pi^j( \theta)^{1/m} d\theta} .
\end{align*}

For each $t$ and $j$, denote $z_t^j = \int_\Theta \bp_\theta^j(x_t^j) \prod_{k = 1}^{t - 1} \prod_{i = 1}^m \bp_\theta^i(x_k^i)^{\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}}  \Pi(d\theta)$ and assume $z_t^j < \infty$.

The distributed Bayes posterior is defined as the probability measure $P_t^j$ such that for $B \in \B$,

\begin{equation} \label{def-db-posterior}
P_{t}^j(B) = \frac{1}{z_t^j}\int_B \bp_\theta^j(x_t^j) \prod_{k = 1}^{t - 1} \prod_{i = 1}^m \bp_\theta^i(x_k^i)^{\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}}  \Pi(d\theta).
\end{equation}

The measure-theoretic definition of distributed Bayes posterior is equivalent to the recursive formulation (<ref>). The definition shows that the distributed Bayes posterior is the product of the prior distribution and a surrogate likelihood.

The asymptotic properties of posteriors with surrogate likelihood have been studied in [Miller, 2021], which establish sufficient assumptions for these posteriors to exhibit posterior concentration, asymptotic normality (Bernstein von-Mises), and asymptotically correct Frequentist coverage.  However, the one key difference between the distributed Bayes posterior and generalized posterior defined in [Miller, 2021] is that the surrogate likelihood examined in [Miller, 2021] depends solely on the sample size $t$. In contrast, the surrogate likelihood used in distributed Bayes posterior depends on the sample size $t$, the statistical model $\bP_\Theta^j$, and the communication graph structure $G$. Therefore, the theoretical analysis of distributed Bayes posterior requires careful adjustments of the standard results for generalized posteriors for new constraints and challenges.

We define the following surrogate loss functions $f_t^j, f_t, f$ on $\Theta$.

\begin{align}
\label{def:f_t^j}
f_t^j(\theta) &= -\frac{1}{t} \log \bp_\theta^j(x_t^j)
- \frac{1}{t} \sum_{k = 1}^{t-1} \sum_{i = 1}^m \left[ \prod_{\tau = k}^{t-1} A_\tau \right]_{ij} \log \bp_\theta^i(X_k^i), \quad j \in [m], \\
\label{def:f_t}
f_t(\theta) &= -\frac{1}{mt} \sum_{k = 1}^t \sum_{i = 1}^m \log \bp_\theta^i(X_k^i), \\
\label{def:f}
f(\theta) &= - \frac{1}{m} \sum_{i = 1}^m \bP_0 \log \bp^i_\theta.
\end{align}

The function  $f_t^j$ corresponds to the surrogate likelihood for the distributed Bayes posterior,  $f_t$ corresponds to the empirical mean of $\log p^j_\theta(X)$ and a special case of $f_t^j$ when the adjacency matrix has weights of all $1/m$. Informally, $f_t^j$ and $f_t$ have the same asymptotic properties under mild assumptions, but $f_t$ has much nicer statistical behaviors, so one should expect distributed Bayes posterior to behave like a generalized posterior with likelihood given by $f_t$.

§ ASSUMPTIONS

Throughout the paper, we impose the following assumptions:

* $\Theta$ is an open subset of $\R^p$ with standard Euclidean metric $d$ .

* There exists a unique parameter $\theta_0 \in \Theta$ that minimizes Problem (<ref>), that is,

\begin{equation}\label{def:theta_0}
\theta_0 = \argmin_{\theta \in \Theta}\sum_{j = 1}^m D_{KL}(\bP_0\parallel \bP_\theta^j).
\end{equation}

The remaining assumptions fall into four categories: communication graph structure, regularity of private statistical models, prior mass assumptions, and consistent testing assumptions. The latter three are standard in Bayes posterior studies with a well-established history (see Appendix Section <ref>). The assumptions on the communication graph are relaxed in Section <ref> to allow for more complex situations.

We assume that the communication graph $G_t, t \geq 1$ is static and undirected, i.e., $G_t = G = (V, E)$.
[Graph Assumptions]
The communication graph $G_t$ is static and undirected, i.e., $G_t = G = (V, E), t \geq 1$. Moreover, $G$ and the adjacency matrix $A$ satisfy

* $A$ is symmetric and row stochastic.

* $A$ has positive diagonal entries, i.e. $a_{ii} > 0$ for all $i \in [m]$.

* $G$ is connected, i.e., a directed path exists from any agent $i$ to any agent $j$.

The assumption of connectedness is standard in the literature on network communication to ensure information flow between agents [Shahrampour et al., 2015, Nedić et al., 2017]. The three assumptions ensure that the Markov chain with transition matrix $A$ is irreducible and aperiodic.

We proceed to establish a result for the convergence of $A^t$ to $\frac{1}{m}\ones \ones^T$ and an upper bound on the convergence rate of $\|A^t - \frac{1}{m}\ones \ones^T\|_1$.

[[Acemoglu et al., 2011, Nedić et al., 2017]]
Let Assumption <ref> hold. The matrix $A$ satisfies

\begin{equation*}
\limsup_{t \to \infty} \sum_{k = 1}^t \sum_{j = 1}^m |[A^{t - k}_{ij}] - \frac{1}{m}| \leq \frac{16m^2 \log m}{\nu}, \qquad i \in [m]
\end{equation*}

where $\nu$ is the smallest positive entry in $A$.

Lemma <ref> states that given fixed $A$, the sum $\sum_{k = 1}^t \sum_{j = 1}^m |[A^{t - k}_{ij}] - \frac{1}{m}|$ is uniformly bounded in $t$, and scales in the order of $m^2 \log m$. The argument for Lemma <ref> relies crucially on the geometric convergence rate of an irreducible, aperiodic Markov chain.  See Section <ref> in the Appendix for the proof.

The next group of assumptions concerns the smoothness of the log-likelihood of the private statistical models.
[Regularity Assumptions]
The regularity assumptions on the statistical model $\bP_\theta^j = \{\bP_\theta^j,  \theta \in \Theta\}$ involves, for every $j \in [m]$ and almost surely $[\bP_0]$,

* For every $\theta \in \Theta$,  $\bP_0 |\log \bp_\theta^j| < \infty$, i.e. the expectation of $|\log \bp_\theta^j(X)|$ is finite under $X \sim \bP_0$.

* The mapping $\theta \mapsto \log \bp_\theta^j$ is convex for every $x$ in a neighborhood of $\theta_0$, and $\nabla f(\theta_0 - \epsilon) < 0 < \nabla f(\theta_0 + \epsilon)$ in every coordinate.

* The statistical model $\bP_\theta^j$ is differentiable in quadratic mean at $\theta_0$ with nonsingular Fisher information matrix $V_{\theta_0}^j$.

* there exists a measurable function $s^j$ with $\bP_0 [s^j(X)]^2 < \infty$ such that, for every $\theta_1, \theta_2$ in a neighborhood of $\theta_0$,

\begin{equation*}
|\nabla \log \bp_{\theta_1}^j - \nabla \log \bp_{\theta_2}^j| \leq s^j(x) \|\theta_1 - \theta_2 \|
\end{equation*}

* the mapping $\theta \mapsto \log \bp_\theta^j$ is twice continuously differentiable for every $x$ in a neighborhood of $\theta_0$, and the Fisher information matrix $V_{\theta_0}^j$ exists and is nonsingular.

* the mapping $\theta \mapsto \log \bp_\theta^j$ is three-times continuously differentiable for every $x$ in a neighborhood of $\theta_0$, the Fisher information matrix $V_{\theta_0}^j$ exists and is nonsingular, and the third-order partial derivatives uniformly bounded by an integrable function in a neighborhood of $\theta_0$.

The regularity conditions described in Assumption <ref> are commonly required prerequisites for establishing the asymptotic properties of posterior distributions [Van der Vaart, 2000]. The underlying intuition is that the differentiability or convexity of the log-likelihood ensures the consistency and asymptotic normality of the maximum likelihood estimators. Asymptotically, the posterior distributions behave like the asymptotic distribution of the MLEs.

The next group of assumptions is typically referred to as the prior mass or prior thickness assumption.
[Prior Assumptions]
For every $j \in [m]$, the prior distribution $\Pi_j$ satisfies

* $\Pi_j(\U_\epsilon) > 0$ for all $\epsilon > 0$, where $\U_\epsilon = \{\theta \in \Theta: \frac{1}{m}\sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_\theta^j)<  \epsilon\}$.

* the density $\pi_j$ is continuous and positive at $\theta_0$.

The first assumption states that prior put a sufficient amount of mass in a KL neighborhood of the target distribution $\bP_0$. An extra continuity assumption of the density $\pi$ at $\theta_0$ is sometimes assumed to connect the first assumption to a statement about the neighborhood of $\theta_0$ [Van der Vaart, 2000].

Uniform, consistent testing assumptions enable $\theta_0$ to be distinguished with a sequence of test functions.  This assumption ensures that asymptotically negligible mass is placed outside a neighborhood of $\theta_0$.
[Uniform Consistent Testing]
Let $\Theta \subseteq \R^p$.  For each $j \in [m]$,

* For all $\epsilon > 0$, there exists $\delta > 0$ such that $\lim_{t \to \infty} \bP_0(\inf_{\|\theta - \theta_0\| > \delta}|f_t^j(\theta) - f_t^j(\theta_0)| \geq \epsilon) = 1$.

* There exists a sequence of test functions $\phi_t$ such that $\bP_{\theta_0}^j (1 - \phi_t(X^{(mt)})) \to 0$ and $\sup_{\|\theta - \theta_0\| \geq \epsilon} \bP_\theta^j (1 - \phi_t(X^{(mt)})) \to 0$ for every $\epsilon$.

Assumption <ref>(b) is attributed to Theorem 10.1 in [Van der Vaart, 2000], and Assumption <ref>(a) is given first in [Ghosh and Ramamoorthi, 2003]. The two assumptions in <ref> play an identical role in the theory but are slightly different. We mainly use assumptions <ref>(a).

§ POSTERIOR CONSISTENCY

This section establishes the posterior consistency for the distributed Bayes posterior.  Theorem <ref>, to be presented next, provides a general concentration result for the distributed Bayes posterior $P_t^j$ over the measurable space $(\Theta, \mathcal{A})$. The proof of Theorem <ref> takes its cues from Schwartz's theorem [Schwartz, 1965]. This fundamental theorem implies two supporting findings. Corollary <ref> shows that Theorem <ref> holds under model misspecification and alternative definitions of the neighborhood around the true parameter. Moreover, Lemma <ref> lays down more user - friend sufficient assumptions for the regularity assumption on the log-likelihood (Assumption <ref>(a)) in Theorem <ref>.

First, we introduce a supporting lemma to demonstrate that the surrogate loss functions $f_t^j$ and $f_t$ are asymptotically equivalent, and both converge to the population loss function $f$. This lemma suggests that under the conditions of a connected communication graph and a first-moment condition satisfied by the statistical model, the distributed Bayes posterior $P_t^j$ becomes asymptotically equivalent to a posterior derived from $f_t$.

Let Assumptions <ref> and <ref>(a) hold. Then $f_t^j$ and $f_t$ converges to $f$ on $\Theta$ in $[\bP_0]-$probability.

The proof can be bound in Section <ref> of the Appendix. The argument is based on a distributed version of the law of large numbers and geometric convergence of the adjacency matrices described in Lemma <ref>.

We are ready to state the main result of this section.

Let $\theta_0 \in \Theta$ be defined in (<ref>), $\epsilon>0$, and denote $\U_\epsilon = \{\theta \in \Theta: \frac{1}{m}\sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_\theta^j)<  \epsilon\}$. Moreover, let Assumptions <ref>, <ref>(a), <ref>(a) hold.  Then, the distributed Bayes posterior $P_t^j$ defined in (<ref>) has the following property:

\begin{equation*}
\lim_{t \to \infty} P_t^j (\U_\epsilon) = 1,
\end{equation*}

in $[\bP_0]-$probability for each $\epsilon > 0$ and $j \in [m]$.

See Section <ref> for the proof. The proof uses the standard technique based on Schwartz's theorem [Miller, 2021]. The primary insight is that, given the defined regularity assumptions on private statistical models and the presence of a connected communication graph, the distributed Bayes posterior concentrates at the same point as the $\frac{1}{m}-$geometric average of the individual Bayes posteriors.

Theorem <ref> does not assume correct model specification. When the model is misspecified, the distributed Bayes posterior concentrates around the unique minimizer $\theta_0$ of Problem (<ref>), which is assumed to exist.

Under additional assumptions,  the neighborhoods on the space of probability measures used in stating Theorem <ref> can be substituted with neighborhoods on $\Theta$.

Let $(\Theta, d)$ be metric space, $\theta_0 \in \Theta$, $\epsilon>0$, and denote $\mathcal N_\epsilon = \{\theta \in \Theta:d(\theta, \theta_0) < \epsilon\}$.   Moreover, let Assumptions <ref>, <ref>(a), <ref>(b), <ref>(c) hold.  Then the distributed Bayes posterior $P_t^j$ defined in (<ref>) has the following property

\begin{equation*}
\lim_{t \to \infty} P_t^j(\mathcal N_\epsilon) = 1,
\end{equation*}

in $[\bP_0]-$probability for each $\epsilon > 0$ and $j \in [m]$.

The result follows directly from Theorem 3 of [Miller, 2021]. Hence, the proof is omitted.

In practice, Assumption <ref>(a) is challenging to verify due to the lack of information on the moments of the unknown distribution. The expectation may not exist under a heavy-tailed distribution like the Cauchy distribution. In the result below, we provide an information-theoretic equivalent of Assumption <ref>(a) that is easier to check.

Assumption <ref>(a) is equivalent to either of the following assumptions:

(a). $\bP_0$ is absolutely continuous with respect to $\bP_\theta^j$ for all $\theta \in \Theta$,

(b). $D_{KL}(\bP_0 \parallel \bP_\theta^j) < \infty$.

When $\bP_0$ has a density $\bp_0$, the absolute continuity assumption is equivalent to the assumption that if $\bp_\theta^j(x) > 0$, then $\bp_0(x) > 0$  for every $x$ [Cover, 1999]. The second assumption is the minimum assumption for the problem of distributed statistical estimation to be tractable. Both assumptions are easy to check, as illustrated by the following example.

Assume that $\bP_\theta^j = \{N(\theta, \sigma^2_j), \theta \in \Theta\}$ for known $\sigma^2_j > 0$ and $\bP_0 = N(\theta_0, \sigma_0^2)$ for unknown mean $\theta_0$ and variance $\sigma_0^2$. The Gaussian distribution with positive variance has support on the whole real line. Therefore, $\bP_0$ is absolutely continuous with respect to $\bP_\theta^j$ for all $j$.  For the second assumption, the KL divergence between $\bP_\theta^j$ and $\bP_0$ has an explicit formula:
The Kullback-Leibler (KL) divergence between two Gaussian distributions $\mathcal{N}(\mu_1, \sigma_1^2)$ and $\mathcal{N}(\mu_2, \sigma_2^2)$ is given by:

\begin{equation*}
D_{KL}(\bP_\theta^j \parallel \bP_0) = \frac{1}{2} \left( \frac{(\theta - \theta_0)^2}{\sigma_0^2} + \frac{\sigma_j^2}{\sigma_0^2} - 1 - \log \frac{\sigma_j^2}{\sigma_0^2} \right).
\end{equation*}

Since $\sigma_j^2, \sigma_0^2, \theta, \theta_0$ are finite, the KL divergence is bounded.

§ ASYMPTOTIC NORMALITY

The Bernstein-Von Mises (BvM) theorem states that under certain conditions on the prior, the posterior distribution approximates a Gaussian distribution centered at a consistent estimator, such as the maximum likelihood estimator, as data increases. This theorem is pivotal in Bayesian statistics for at least two reasons. First, it provides a quantitative description of how the posterior contracts to the truth. Second,  this theorem justifies Bayesian credible sets as valid Frequentist confidence sets, i.e., sets of posterior probability $1 - \alpha$ contain the true parameter at the confidence level $1 - \alpha$.

In this section,  we establish Bernstein von - Mises theorems for the distributed Bayes posterior defined in Equation (<ref>). Our results address the asymptotic normality of the distributed Bayes posterior under both correct and incorrect model specifications. We consider scenarios where all agents accurately specify the model, and some agents gather observations from a 'true distribution' that does not belong to their statistical models.

Theorem <ref> provides sufficient assumptions for the distributed posterior to converge to a normal distribution centered around a sequence of M-estimators $\hat \theta_t^j$. Theorem <ref> generates an analogous result under classical assumptions, with the normal approximation centered around $\theta_0$ (defined in Equation (<ref>)). The supporting lemmas are provided preceding the main results. Given the abstract nature of the assumptions, Theorem <ref> is supplemented by Corollary <ref> and Corollary <ref>, which outlines more user-friendly assumptions that guarantee the same Bernstein von-Mises (BvM) results.

The Bernstein von - Mises argument relies critically on the following sequence of M-estimators:

\begin{equation} \label{eqn:m-estimator}
\hat \theta_t^j = \argmin_{\theta \in \Theta} f_t^j(\theta).
\end{equation}

The existence and consistency of $\hat \theta_t^j$ is not always guaranteed; one set of sufficient assumptions is provided as follows.

Let Assumptions  <ref>, <ref>(a), <ref>(f), <ref>(b) hold.  Then the probability that the equation $\nabla f_t^j(\hat \theta_t^j) = 0$ has at least one solution converges to $1$ as $t \to \infty$, and there exists a sequence of solutions $\hat \theta_t^j$ such that $\hat \theta_t^j \pto \theta_0$.

This is a direct consequence of Theorem 5.42 of [Van der Vaart, 2000]; thus, the proof is omitted. Lemma <ref> states that the sequence of M estimators exists under third-order smoothness assumptions on the log - log-likelihood.  From now on, we assume the existence of M estimators $\hat \theta_t^j$ that satisfies (<ref>) for every $t, j$.

The high-order smoothness assumptions on the private log-likelihoods (Assumption (f)) is a strong assumption for guaranteeing consistency. It can be replaced with a more relaxed and amenable convexity assumption (Assumption <ref>(b)).

Let Assumptions  <ref>, <ref>(a), <ref>(b) hold. Then $\hat \theta_t^j \pto \theta_0$.

The next two lemmas are useful in the proof of Theorem <ref>.

Let Assumptions <ref>, <ref>(c), <ref>(a) hold and $\hat \theta_t^j \pto \theta_0$. Then for every $\epsilon > 0$, there exists $\delta > 0$ such that

\begin{equation*}
\lim_{t \to \infty} \bP_0(\inf_{\|\theta - \hat\theta_t^j\| > \delta}|f_t^j(\theta) - f_t^j(\hat\theta_t^j)| \geq \epsilon) = 1
\end{equation*}

The result states that if a sequence of estimators $\hat \theta_t^j$ is consistent at $\theta_0$, then the existence of uniformly consistent tests at $\theta_0$ (Assumption <ref>(a)) implies analogous result at $\hat \theta_t^j$.

Let $\hat \theta_t \in \R^p$ such that $\hat \theta_t \pto \theta_0$ for some $\theta_0 \in \R^p$, $\pi_t$ be a density with respect to Lebesgue measure on $\R^p$. Suppose $q_t$ is the density of $a_t (\theta - \hat \theta_t)$ where $\theta \sim \pi_t$ and $a_t^{-1} = o(1)$. If $\int |q_t(x) - q(x)|dx \to 0$ in $[\bP_0]-$probability for some probability density $q$, then for every $\epsilon > 0$,

\begin{equation*}
lim_{t \to \infty}\Pi_t(\mathcal N_\epsilon) = 1,
\end{equation*}

where $N_\epsilon$ is the neighborhood of $\theta_0$ defined in Corollary <ref>.

The proofs of Lemma <ref>,  Lemma <ref>,  Lemma <ref>, and Lemma <ref> can be found in Section <ref> of the Appendix. Lemma <ref>  and <ref> relies on the usual consistency argument of M - estimators. Lemma <ref> and Lemma <ref> use measure-theoretic arguments.

Our first main result of the section provides general sufficient assumptions under which a distributed Bayes posterior exhibits asymptotic normality and an asymptotically correct Laplace approximation, along
with the posterior concentration at $\theta_0$.

Let $\Theta$ be an open subset of $\R^p$. Assume that there exists $\theta_0 \in \Theta$ such that $P_0 = P_{\theta_0}^j$ for every $j \in [m]$. Moreover, let Assumptions <ref>,  <ref>(a), <ref>(c), <ref>(d), <ref>(b), <ref>(a)  hold.  If the sequence $\hat \theta_t^j$ defined in (<ref>) satisfies that $\hat \theta_t^j \pto \theta_0$, then

\begin{equation} \label{eqn-bvm-1}
f_t^j(\theta)  = f_t^j(\hat \theta_t^j) - \frac{1}{2}(\theta - \hat \theta_t^j)^T \hat V_t^j (\theta -\hat \theta_t^j) + r_t^j(\theta -\hat \theta_t^j),
\end{equation}

where $\hat V_t^j$ is a sequence of matrices that converges in probability to the average Fisher information matrix $V_{\theta_0} = \frac{1}{m} \sum_{i = 1}^m V_{\theta_0}^i$, and $|r_t^j(h)| = O(|h|^3)$ for large enough $t$.

Let Equation (<ref>) and Assumption <ref>(b), <ref>(a) hold. As $t \to \infty$, we have

\begin{equation} \label{eqn-bvm-2}
\int_{B_\epsilon(\theta_0)} p_t^j(\theta) d\theta \pto 1 \quad \forall \epsilon >0.
\end{equation}

that is, the distributed Bayes posterior $p_t^j$ is weakly consistent around $\theta_0$.

Let $q_t^j$ be the density of $\sqrt{t}(\theta -\hat \theta_t^j)$ when $\theta \sim P_t^j$. Then,

\begin{equation} \label{eqn-bvm-3}
\int_\Theta |q_t^j(x) - N(0, V_{\theta_0}^{-1})| dx \pto 0.
\end{equation}

i.e., the total variational distance between $q_t^j$ and $N(0, V_{\theta_0}^{-1})$ vanishes in $[\bP_0]-$probability.

See Section <ref> for the proof.
The assumptions of Theorem <ref> involve structural and statistical assumptions that are standard within the Bernstein von Mises literature. These include a connected communication graph (Assumption <ref>), bounded entropy condition of the private statistical models (Assumption <ref>(a)), a distributed version of the differentiable in quadratic means (DQM) assumption (Assumption <ref>(c)) and the Lipschitz gradient regularity assumption on the log-likelihood around $\theta_0$ (Assumption <ref>(d)). The M-estimators, $\hat \theta_t^j$, are assumed to be consistent, with sufficient conditions outlined in Lemmas <ref> and <ref>. These regularity assumptions are fairly mild and applicable to, for example, most exponential family models. They serve as the foundation for Equation (<ref>), which mirrors the local asymptotic normality assumption in classical BvM theory [Cam et al., 2000]. Additionally, we specify a prior mass assumption (Assumption <ref>(b)) and an assumption related to uniform, consistent testing (Assumption <ref>(a)). The prior mass assumption (Assumption <ref>(b)) is slightly different from the one used in Theorem <ref>, but they are equivalent when $f$ is continuous at $\theta_0$. These, together with Equation (<ref>), are the main sufficient assumptions for the BvM results in Equations (<ref>) and (<ref>).

The result shows a sampling rate of $\sqrt{t}$, meaning no agent loses parametric efficiency when performing distributed inference in the fully connected network compared to the single-node network. This helps ensure the scalability of the approach. However, the main benefit of having communication between agents is the robustness of the uncertainty quantification. In Theorem <ref>, the normal approximation to the distributed Bayes posterior uses the inverse of the average Fisher information across all agents as a covariance matrix. The intuition is that the individual agent's uncertainty about the ground truth eventually aligns with the network's average level of uncertainty, achieved irrespective of the agent's initial uncertainty or the particular statistical model they employ. Notably, this attribute is a protective mechanism in scenarios where a few agents might be subject to anomalous levels of uncertainty, potentially due to adversarial activities or other forms of data corruption.

The differentiability in quadratic means (DQM) assumption (Assumption <ref>(c)) in Theorem <ref> may be difficult to verify in practical settings. We replace the abstract DQM assumptions with a second-order smoothness assumption on the private log-likelihoods Assumption <ref>(e).

Theorem <ref> holds if  Assumption <ref>(c) is replaced with Assumption <ref>(e).

Both differentiability in quadratic means (DQM) assumption (Assumption <ref>(c)) and the Lipschitz gradient assumption (Assumption <ref>(d)) can be replaced with a third-order smoothness assumption (Assumption <ref>(f)) which is often called the classical condition for asymptotic normality of M - estimators [Van der Vaart, 2000].

Theorem <ref> holds if  Assumption <ref>(c) and <ref>(d) are replaced with Assumption <ref>(f).

The proofs of Corollary <ref> and Corollary <ref> are based on bounding the second and third-order terms in the Taylor expansion of $f_t^j(\theta)$, respectively. We can substitute more model regularity assumptions by bounding higher-order terms in the Taylor expansion.

The Bernstein-von Mises theorem, as traditionally stated [Van der Vaart, 2000], relies on a stochastic rendition of Local Asymptotic Normality (LAN). Here, we introduce a distributed version of stochastic LAN. A distributed family of statistical models $\left(\{\bP_\theta^j\}_{j \in [m]}, G \right)$ is stochastic LAN  at $\theta \in \Theta$ if, for every $j \in [m]$, and in relation to a non-singular scaling factor $\epsilon_t^j \to 0$, there exists a random vector $\Delta_{t, \theta}^j$ and a non-singular matrix $V_\theta$ such that $\Delta_{t, \theta}^j$ is bounded in probability and for every compact subset $K \subset \R^p$: as $t \to \infty$,

\begin{equation} \label{eqn:s-lan}
\sup_{h \in K} \left| - t f_t^j(\theta + \epsilon_t^j h) + t f_t^j(\theta)  -h^T V_\theta \Delta_{t, \theta}^j   + \frac{1}{2} h^T V_\theta h \right| \pto 0.
\end{equation}

The random vector $\Delta_{t, \theta}^j$ is often called the “local sufficient statistics,” and $V_\theta$ typically corresponds to the average of Fisher information. It's worth noting that the $V_\theta$ matrix must be the same across all agents. The scaling factor $\epsilon_t^j$ is chosen to ensure $\Delta_{t, \theta}^j$ is asymptotically normal and $V_\theta$ converges to a nonsingular matrix. The default choice is $\epsilon_t^j = \frac{1}{\sqrt{t}}$. However, an agent-specific choice is available if any private statistical model has a different parametric concentration rate.

A sufficient assumption for the distributed LAN is a distributed version of the DQM condition.

Let Assumptions <ref>, <ref>(a), <ref>(c) hold. The distributed family of statistical models, denoted as $\left(\{\bP_\Theta^j\}_{j \in [m]}, G \right)$, is stochastically locally asymptotically normal (LAN) at $\theta_0$.

We now relax the assumption of the correct model specification to explore the more realistic scenario where the data-generating processes are not perfectly aligned with the assumed private statistical models. Specifically, we extend our previous results to the case where the true data-generating distribution $\bP_0 \neq \bP_{\theta_0}^j$ for some $j \in [m]$.

Let $\theta_0 \in (\Theta, d)$ be defined in (<ref>). Moreover, let Assumptions <ref>(b), <ref>(b) hold, and assume that the distributed stochastic LAN (<ref>) holds at $\theta_0$. Let a sequence of constants $\epsilon_t^j$ satisfy that for any $M_t \to \infty$,

\begin{equation} \label{assumption:bvm2}
P_t^j\left(\theta: d(\theta, \theta_0) \geq M_t \epsilon_t^j \right) \pto 0.
\end{equation}

If $q_t^j$ is the density of $(\theta - \theta_0)/\epsilon_t^j$ when $\theta \sim P_t^j$, then

\begin{equation*}\label{eqn:bvm2-main-result}
\int_\Theta |q_t^j(x) - N(0, V_{\theta_0}^{-1})| dx \pto 0,
\end{equation*}

for $V_{\theta_0}$ provided in (<ref>).

Since the covariance matrix $V_{\theta_0}$ in the misspecified Bernstein-Von Mises theorem fails to match the average of sandwich covariance matrices, the posterior credible sets derived from $P_t^j$ have valid Frequentist coverage. While these sets may be properly centered at $\hat \theta_t^j$, their width may be inaccurate, and they don't typically correspond to confidence sets with level $1 - \alpha$.

The sequence $\epsilon_t^j$ that satisfies assumption (<ref>) is called a posterior contraction rate of the distributed Bayes posterior $P_t^j$. This quantity determines the convergence rate of $P_t^j$ to the unknown parameter $\theta_0$. Results to control the contraction rates are provided in the next section.

§ CONTRACTION RATE

Contraction rates quantify the speed at which a posterior distribution approaches the true parameter of the data-generating distribution. Controlling the contraction rates not only refines our understanding of posterior consistency but also helps control the sampling complexity in the misspecified Bernstein von - Mises Theorem (Theorem <ref>). Unlike the previous sections, which focus on asymptotic results, we provide non-asymptotic bounds and scaling laws of the posterior contraction rates in this section. Our results involve the sample size $t$, dimension $p$, and the number of agents $m$.

For two positive sequences $x_t$ and $y_t$, we use $x_t \lesssim y_t$ to denote the existence of a constant $c$, independent of $n$, such that $x_t \leq c y_t$. Furthermore, we write $x_t \asymp y_t$ when $x_t \lesssim y_t$ and $y_t \lesssim x_t$.

Let $(\Theta, d)$ be a metric space. A sequence of constants $\epsilon_t^j$ is a posterior contraction rate for $P_t^j$ at the parameter $\theta_0$ if, for every $M_t \to \infty$,

\begin{equation} \label{def:contraction-rate}
P_t^j\left(\theta: d(\theta, \theta_0) \geq M_t \epsilon_t^j \right) \pto 0,
\end{equation}

The posterior contraction rate is not a unique quantity. Any rate slower than a contraction rate is also a contraction rate. Although the fastest rate is desirable, they may be hard to find. The natural goal is to establish a rate that is a close upper bound for the “optimal” rate. We refer to this upper bound as the contraction rate.

We define a probability measure $\bP_\theta$ on the product space $\X^m$. For a measurable set $A \subseteq \X^m$, we have

\begin{equation*}
\bP_\theta(A) = \frac{\int_A  \prod_{j = 1}^m [\bp_{\theta_0}^j(x^j)]^{\frac{1}{m}} dx^1 \ldots d x^m}{\int_\Theta  \prod_{j = 1}^m [\bp_{\theta_0}^j(x^j)]^{\frac{1}{m}} dx^1 \ldots d x^m}.
\end{equation*}

The density is given by $\bp_\theta(x) = \prod_{j = 1}^m [\bp_{\theta_0}^j(x^j)]^{\frac{1}{m}}$ for $x \in \X^m$. Given the prior measure $\Pi$,  we define the distributed ideal posterior $P_t$ as the posterior measure corresponding to $\bP_\theta$ and $\Pi$. For a measurable set $B \subseteq \Theta$, we have

\begin{equation} \label{def-ideal-posterior}
P_t(B) =\frac{\int_B \prod_{k = 1}^t \bp_\theta(x_k) \Pi(d\theta)}{\int_\Theta \prod_{k = 1}^t \bp_\theta(x_k) \Pi(d\theta)} = \frac{\int_B \prod_{k = 1}^t \prod_{j = 1}^m [\bp_\theta^j(x_k^j)]^{\frac{1}{m}}\Pi(d\theta)}{\int_\Theta \prod_{k = 1}^t \prod_{j = 1}^m [\bp_\theta^j(x_k^j)]^{\frac{1}{m}}\Pi(d\theta)},
\end{equation}

Let $D_{\rho}(p \parallel q)$ denote the $\rho-$Rényi divergence, as defined in Section (<ref>) of the Appendix.
We establish that the contraction rate of the distributed Bayes posterior$P_t^j$ is given by:

\begin{equation} \label{eqn:contraction-rate}
\epsilon_{m,t}^2 + \frac{1}{mt} \bP_0 D_{\text{KL}}(P_t^j \parallel P_t) + \frac{1}{m} \sum_{i = 1}^m D_{\text{KL}}(\bP_0 \parallel \bP_{\theta_0}^i).
\end{equation}

Here $\epsilon_{m,t}$ is the contraction rate of the distributed ideal posterior $P_t$. The second term quantifies the approximation error when the distributed Bayes posterior is approximated by the corresponding distributed ideal posterior under the true distribution $\bP_0$. The last term captures the average divergence between the model specifications from individual agents and the true underlying distribution.

The main theorem of the section is stated under the “prior mass and testing” framework. The assumptions are sub-exponential refinements of the assumptions for Bernstein von - Mises theorems: (a) The prior is required to put a minimal amount of mass in a neighborhood of the true parameter. (b) Restricted to a subset of the parameter space, there exists a test function that can distinguish the truth from the complement of its neighborhood; (c) The prior is essentially supported on the subset described in (b).

Suppose $\epsilon_{m,t}$ is a sequence such that $mt\epsilon_{m,t}^2 \geq 1$. Let $C_0$, $C_1$, $C_2$, $C_3 > 0$ be constants such that $C_0 > C_2 + C_3 + 2$. Let the following assumptions hold:

* For any $\epsilon > \epsilon_{m,t}$, there exists a set $\Theta_t(\epsilon)$ and a testing function $\phi_t$ such that:

\begin{align}
\bP_{\theta_0} \phi_t(X^{(mt)}) + \sup_{\substack{\theta \in \Theta_t(\epsilon)\\ d(\theta, \theta_0) \geq C_1 \epsilon^2}} \bP_\theta (1 - \phi_t(X^{(mt)})) &\leq \exp\left(-C_0 t \epsilon^2\right). \tag{C1} \label{assumption:C1}
\end{align}

* For any $\epsilon > \epsilon_{m,t}$, the set $\Theta_t(\epsilon)$ above satisfies:

\begin{equation}
\Pi\left(\Theta_t(\epsilon)^c\right) \leq \exp\left(-C_0 t \epsilon^2\right). \tag{C2} \label{assumption:C2}
\end{equation}

* For some constant $\rho > 1$:

\begin{equation}
\Pi\left(\theta \in \Theta,  \frac{1}{m} \sum_{j = 1}^m D_\rho(\bP_{\theta_0}^j \parallel \bP_\theta^j) \leq C_3 \epsilon_{m,t}^2 \right) \geq \exp\left(-C_2 t \epsilon_{m,t}^2\right). \tag{C3} \label{assumption:C3}
\end{equation}

Then, for the distributed Bayes posterior $P_t^j$ defined in (<ref>), we have:

\begin{equation}\label{eqn:contraction11}
\bP_0 P_t^j d(\theta, \theta_0) \leq C (\epsilon_{m,t}^2 + \gamma_{j,m,t}^2+ \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j)),
\end{equation}

for some constant $C$ depending on $C_0, C_1$, where the quantity $\gamma_{j,m,t}^2$ is defined as:

\begin{equation*}
\gamma_{j,m,t}^2 =  \frac{1}{mt} \bP_0 D_{KL}(P_t^j \parallel P_t).
\end{equation*}

The proof of Theorem <ref> can be found in Section <ref> of the Appendix  The proof of Theorem <ref> is a consequence of four support lemmas which make use of the Gibbs variational representation and the subexponential decay of sub-exponential decay of $d(\theta, \theta_0)$ under $P_t$.

Assumption (<ref>) and  (<ref>) is a refinement of assumptions <ref> for the uniform consistent testing and states that there is a sequence of tests such that the sum of Type I and Type II errors decrease exponentially with sample size, where the alternative hypothesis is taken in a large enough set under the prior. Assumption (<ref>) refines the prior mass assumptions <ref> by stating that the prior mass decreases exponentially away from a $\rho-$Rényi neighborhood of the true distribution. This assumption is slightly stronger than the equivalent assumption stated with the KL neighborhood because $D_\rho(P \parallel Q) > D_{KL}(P \parallel Q)$ for $\rho >1 $.

The contraction rate is the sum of two terms, $\epsilon_{m,t}$ and $\gamma_{j,m,t}^2$. The first term is the contraction rate of the distributed ideal posterior. The second term $\gamma_{j,m,t}^2$ characterizes the distance between the distributed Bayes posterior $P_t^j$ and the distributed ideal posterior $ P_t$. A larger or less connected communication graph means potentially more deviation between the two distributions, and thus, the rate $\gamma_{j,m,t}^2$ is larger.

We can show by Markov inequality that the upper bound on $ \bP_0 P_t^j d(\theta, \theta_0)$ is indeed the contraction rate for the distributed Bayes posterior $P_t^j$. We obtain a point estimator $\hat \theta$ that converges to the KL minimizer $\theta_0$ for convex loss functions at the same rate. The following result elaborates on this.

Under the assumptions of Theorem <ref>, for any diverging sequence $M_t \to \infty$, we have

\begin{equation*}
\bP_0 P_t^j \left(d(\theta, \theta_0) > M_t (\epsilon_{m,t}^2 + \gamma_{j,m,t}^2 +  \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j)) \right) \to 0
\end{equation*}

Furthermore, if $d(\theta, \theta_0)$ is convex in $\theta$, the distributed posterior posterior mean $\hat \theta = \int_\Theta \theta dP_t^j(\theta)$ satisfies

\begin{equation*}
\bP_0  d(\hat \theta, \theta_0) \leq C(\epsilon_{m,t}^2 + \gamma_{j,m,t}^2 +  \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j))
\end{equation*}

where $C$ is the same constant in (<ref>).

The contraction rate described in Theorem <ref> remains somewhat abstract because it is characterized by terms such as $\epsilon_{m,t}^2$ and $\gamma_{j,m,t}^2$ which do not directly inform the design of the underlying communication network. For practical purposes, it's preferable to express the rate in terms of parameters such as \(m\), \(t\), and \(\nu\), which can guide the design of a statistically efficient network structure. To this end, the following section offers more concrete upper bounds on the terms \(\epsilon_{m,t}^2\) and \(\gamma_{j,m,t}^2\).

The upper bounds for $\epsilon_{m,t}^2$ can be directly borrowed from the existing theory on posterior contraction rates under model misspecification [Kleijn and van der Vaart, 2012]. For finite-dimensional models denoted by $\{ \mathbb{P}_\theta, \theta \in \Theta \}$, the optimal contraction rate is given by $t^{-1/2}$. However, this result does not follow from Theorem <ref> because it requires a more restrictive metric entropy assumption involving the $2^{nd}$-order KL divergence. See, for example, Theorem 2.2 of [Kleijn and van der Vaart, 2012].

The general theory for posterior contraction rates typically combines a prior mass assumption, often in the form of <ref>, with either a model entropy assumption or a consistent testing assumption in the form of <ref> and <ref>. For a review of the theory of posterior contraction rates, see Chapter 8 of [Ghosal and Van der Vaart, 2017] and the references therein.

As touched upon in Section <ref>, one should expect the distributed Bayes posterior to be well approximated by the distributed ideal posterior as the sample size increases. The next result provides a uniform bound on the approximation error $\gamma_{j,m,t}^2$.

Let Assumptions <ref> and <ref>(a) hold. For $P_t^j$ defined in (<ref>) and $P_t$ defined in (<ref>), we have

\begin{equation*}
\gamma_{j,m,t}^2 \leq  \frac{16m \log m}{\nu t} \left(|\bP_0\log \bp_0| +  \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right)
\end{equation*}

The upper bounds on $\gamma_{j,m,t}^2$ consist of two terms: the first term depends on the structure of the graph, and the second term depends on the graph and the worst-case model misspecification error, denoted by $\max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i)$. The intuition is that there is a tradeoff between the size of the communication network and statistical efficiency, and the communication cost could be much higher if one or more agents in the network use misspecified models. If all models are correctly specified, the contraction rate degrades with $m$ at a rate of $m \log m$. If the term $\max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) $ scales at a rate of $m^2$, then the contraction rate decreases at a much higher rate of $m^3 \log m$.

Suppose we ignore the constants and focus on the scaling law. In that case, the contraction rate of the distributed Bayes posterior is a function of the sample size, the number of agents, the smallest positive adjacency weight (spectral gap), the worst-case model misspecification error, and the average model misspecification error. We formalize this in the following result.

Let Assumptions <ref>,  <ref>(a), and the assumptions of Theorem (<ref>) hold. Let the distributed ideal posterior $P_t$ satisfies a contraction rate of $\epsilon_{m,t}^2 \lesssim t^{-1}$. For the distributed Bayes posterior $P_t^j$ defined in (<ref>), we have:

\begin{equation*}
\bP_0 P_t^j d(\theta, \theta_0) \lesssim \frac{1}{t} + \frac{m \log m}{\nu t} \left(1 + \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right)  + \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j).
\end{equation*}

The result follows directly from Theorem <ref> and <ref>. Thus, the proof is omitted.

The last theorem outlines how the design ingredients determine the contraction rate of a distributed Bayes posterior. Specifically, the contraction rate slows as the number of agents increases by a factor of $m \log m$. The rate is further hindered linearly by the worst-case misspecification error. Unlike the other two factors, the average error term remains constant and does not vanish as the sample size increases. It sets a 'baseline' level of deviation between the asymptotic behavior of the posterior and the true parameter value. The three factors collectively determine the statistical efficiency of the distributed Bayes posterior, which highlights the complexities involved in its sampling performance.

§ EXTENSION TO TIME-VARYING GRAPHS

This section extends the theory to sequences of time-varying graphs, represented as $G_t = (V, E_t), t \geq 1$. While the node set $V$ remains constant, the edge set $E_t$ undergoes temporal changes. The principal aim of this exploration is to understand the statistical attributes of the distributed Bayes posterior in a fluctuating communication landscape.  We focus on one scenario where agents communicate within a fully connected graph with a fixed probability and have no communication otherwise. This setup provides a structured yet flexible framework for analyzing the implications of time-variant connectivity on the statistical efficiency of distributed Bayesian inference.

the identity matrix $I_m$.
Let Assumption <ref> be satisfied for graph $G$ with the adjacency matrix $A$ We assume that $G_t$ and $A_t$ are independent, random graphs and matrices such that $G_t = G$ and $A_t = A$ with probability $\lambda \in [0, 1]$, respectively, and $A_t = I_m$ otherwise.

Assumption <ref> proposes a setting useful for various reasons. Theoretically, this setting leads to a rigorous study of the tradeoff between statistical efficiency and communication cost. It allows us to explore what constitutes an "optimal" level of communication, given a targeted level of statistical efficiency. On the practical side, this setting can be conceptualized as a network experiencing complete failure with probability $\lambda$. In such scenarios, the primary objective is to quantify how this level of intermittent connectivity impacts the overall learning quality within the system.

We provide an analogous result to Lemma <ref> in the Assumption <ref> setting.

Let Assumption <ref>, hold and $m \geq 2$. Then the sequence of adjacency matrices $(A_t)$ satisfies the following scaling law with probability $1$:

If $\lambda \geq \frac{2}{m}$, then

\begin{equation*}
\limsup_{t \to \infty} \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}- \frac{1}{m} \right| \leq  \frac{16 m^2 \log m + 8m^2 \log \lambda}{\lambda \nu}, \qquad i \in [m]
\end{equation*}

If $0<\lambda < \frac{2}{m}$, then

\begin{equation*}
\limsup_{t \to \infty} \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}- \frac{1}{m} \right| \leq  \frac{4m^3}{\nu}, \qquad i \in [m]
\end{equation*}

If $\lambda = 0$, then

\begin{equation*}
\limsup_{t \to \infty} \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}- \frac{1}{m} \right|= \infty, \qquad i \in [m]
\end{equation*}

In these inequalities, $\nu$ denotes the smallest positive entry in $A$.

Proposition <ref> provides the scaling laws for three regimes of communication frequency.  Consensus can be reached as long as the communication occurs for a non-zero proportion of the time. However, suppose the goal is to achieve optimal scaling behavior relative to the size of the communication graph ($m$). In that case, the communications need to happen with a minimum frequency of $\frac{2}{m}$, with a higher frequency being more desirable. On the other hand, if there is a need to adhere to the low-frequency regime ($0<\lambda < \frac{2}{m}$) due to constraints such as high communication costs, the optimal strategy becomes one of minimal communication. Notably, within this regime, a decrease in communication frequency does not affect the rate of attaining consensus among the agents.

Recall that the contraction rate of $P_t^j$ is decomposed into three pieces: the contraction rate of the distributed ideal posterior, the approximation error, and the average misspecification error across the network.

\begin{equation*}
\epsilon_{m,t}^2 + \frac{1}{mt} \bP_0 D_{\text{KL}}(P_t^j \parallel P_t) + \frac{1}{m} \sum_{i = 1}^m D_{\text{KL}}(\bP_0 \parallel \bP_{\theta_0}^i).
\end{equation*}

The structure of the communication network affects the contraction rate through the second term. The following Corollary illustrates this.

Let Assumptions <ref> and <ref>(a) hold. For $P_t^j$ defined in (<ref>) and $P_t$ defined in (<ref>), the following holds with probability $1$:

If $\lambda \geq \frac{2}{m}$, then

\begin{equation*}
\frac{1}{mt} \bP_0 D_{\text{KL}}(P_t^j \parallel P_t) \leq  \frac{16 m \log m + 8m \log \lambda}{\lambda \nu t} \left(|\bP_0\log \bp_0| +  \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right).
\end{equation*}

If $0<\lambda < \frac{2}{m}$, then

\begin{equation*}
\frac{1}{mt} \bP_0 D_{\text{KL}}(P_t^j \parallel P_t) \leq  \frac{4m^2}{\nu t}\left(|\bP_0\log \bp_0| +  \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right).
\end{equation*}

In Corollary <ref>, the term "with probability \(1\)" specifies the likelihood under the measure defined over the sequence \(A_1, A_2, \ldots\), where each \(A_t\) is an independent random matrix that toggles between two pre-determined values with a probability of \(\lambda\). Importantly, this is the same probability measure that underlies Proposition <ref>, and we continue to employ this measure in the forthcoming results.

Different regimes of the communication frequency $\lambda$ lead to a difference in the statistical performance. In the following result, we demonstrate the impact of $\lambda$ on the contraction rate of the distributed Bayes posterior.

Let Assumptions <ref>,  <ref>(a), and the assumptions of Theorem <ref> hold.Let $\epsilon_{m,t}^2 \lesssim t^{-1}$. Then for the distributed Bayes posterior $P_t^j$ defined in (<ref>), the following holds with probability $1$:

If $\lambda \geq \frac{2}{m}$, then

\begin{equation*}
\bP_0 P_t^j d(\theta, \theta_0) \lesssim  \frac{1}{t} + \frac{m \log m + m \log \lambda}{\lambda \nu t} \left(1 +  \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right) + \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j).
\end{equation*}

If $0<\lambda < \frac{2}{m}$, then

\begin{equation*}
\bP_0 P_t^j d(\theta, \theta_0) \lesssim  \frac{1}{t} +\frac{m^2}{\nu t} \left(1 +  \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right) + \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j).
\end{equation*}

The distributed ideal posterior $P_t$ does not depend on the communication graph, thus Theorem <ref> holds and the assumption that $\epsilon_{m,t}^2 \lesssim t^{-1}$ is still valid. Since the only term in the contraction rate that depends on the communication graph is $\gamma_{j,m,t}^2$, the argument for Corollary <ref> then directly builds on Theorem <ref> and Corollary <ref>.  See Section <ref> in the Appendix for the proof.

This is the first result we know of on the impact of time-varying communication networks on the efficiency of distributed statistical inference. The implication deviates from the conventional wisdom about the relationship between communication cost and statistical efficiency.  Naturally, one might expect a communication-statistical tradeoff because increasing the communication frequency leads to faster posterior contraction rates. While higher communication frequencies ($\lambda \geq \frac{2}{m}$) do lead to accelerated contraction rates, scaled by a factor of $\frac{\log \lambda}{\lambda}$, the story is different when the frequency is below $\frac{2}{m}$. In this regime, reducing communication costs ceases to have a linear effect on the contraction rate. Instead, the contraction rate scales quadratically with the number of agents $m$—at a rate of $m^2$ rather than $m \log m$.

This suggests the possibility of a “phase transition” phenomenon at $\lambda = \frac{2}{m}$, worth exploring for future research. In practical terms, if the goal is to balance communication efficiency and statistical performance in a large network, a communication frequency of $\lambda^* = \frac{2}{m}$ appears to be the optimal choice.

§ ILLUSTRATIVE EXAMPLES

§.§ Exponential Family Distributions

Let the distributed statistical models $\left(\{\bP_\Theta^j\}_{j \in [m]}, G \right)$ be well - specified, and the Assumption <ref> be satisfied for $G$. Let $\eta: \Theta \to \R^p$ and $T: \mathcal{X} \to \R^p$ be some sufficient statistics and let $\psi^j: \Theta \to \R, h: \X \to \R$ be normalizing functions. We assume that $\bP_\theta^j$ is a member of the canonical exponential family

\begin{equation} \label{eqn-EF-1}
p_\theta^j(x) = h(x) \exp\left(\langle \theta, T^j(x)\rangle - \psi^j(\theta) \right).
\end{equation}

The exponential family includes commonly used Gaussian, exponential, gamma, chi-square, Beta, Dirichlet, Bernoulli, categorical, Poisson, Wishart, inverse Wishart, and geometric distributions. In this section, we only consider the canonical exponential family and interchangeably refer to the exponential and canonical exponential families.

Exponential family distributions are helpful for Bayesian inference because of the conjugate properties, i.e., the posterior is a member of the exponential family if the prior and likelihood are members of the exponential family. This property is preserved in the distributed setting.

Let us consider a scenario where the beliefs of all agents at time $t$ belong to the natural exponential family. In this setting, the belief of agent$i$ concerning the parameters$\theta$ can be described as follows:

\begin{equation*}
p_t^i(\theta) \propto \exp\left( \langle \theta,  \chi_t^i \rangle - B^i(\theta) \right),
\end{equation*}

where$\chi_t^i$ represents the sufficient statistic for agent$i$ at time$t$, and$A^i(\theta)$ is the log-partition function.

We update the one-step-ahead posterior $p_{t+1}^j$ with the distributed Bayes rule,

\begin{equation*}
\begin{aligned}
p_{t+1}^j (\theta) &\propto \log \bp_\theta(X_{t+1}^j) \prod_{i = 1}^m  p_t^i(\theta) \\
&\propto \exp(\langle \theta, T(X_{t+1}^j) \rangle - \psi^j(\theta))\exp(\langle \theta,\sum_{i = 1}^m A_{ij} \chi_t^i\rangle - \sum_{i = 1}^m B^i(\theta)) \\
&\propto  \exp\left(\langle \theta,T^j(X_{t+1}^j)  +  \sum_{i = 1}^m A_{ij} \chi_t^i \rangle - \psi^j(\theta) - \sum_{i = 1}^m B^i(\theta) \right).
\end{aligned}
\end{equation*}

The distribution $p_{t+1}^j$ is a member of the exponential family with sufficient statistic $T^j(X_{t+1}^j)  +  \sum_{i = 1}^m A_{ij} \chi_t^i$ and the log-partition function $\psi^j(\theta) + \sum_{i = 1}^m B^i(\theta)$. This provides an easy-to-implement algorithm to learn the distributed Bayes posterior.

Let the prior be a member of the natural exponential family.

\begin{equation} \label{eqn-EF-2}
\pi(\theta) = g(u) \exp\left(\langle \theta, u\rangle - \psi^0(\theta) \right).
\end{equation}

Leveraging the conjugate property of the exponential family, we derive a closed-form expression for the density of the distributed Bayes posterior $P_t^j$.

Let assumption (<ref>) hold. Assume that the likelihood $P_\theta^i$ has an exponential family form given by Equation (<ref>). Assume that the prior $\Pi$ has an exponential family form given by Equation (<ref>).  Then the distributed Bayes posterior $P_t^j$ defined in (<ref>) is given by the following formula:

\begin{equation}\label{eqn-EF-3}
p_t^j (\theta) = h(X^{(mt)}) \exp(\langle \theta, \chi_t^j + u\rangle - B_t^j(\theta) - \psi^0(\theta)),
\end{equation}

\begin{equation}\label{eqn-EF-4}
\chi_t^j =   \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] T^i(X_k^i),  \quad B_t^j(\theta) = \sum_{i = 1}^m \sum_{k = 1}^t [A^{t-k}_{ji}] \psi^i(\theta) ,
\end{equation}

are the sufficient statistic and the log-partition function, respectively.

An exponential family is full-rank if no linear combination of the sufficient statistic is constant. For example, we say the distribution $p_\theta^i$ in Equation (<ref>) is full-rank if no linear combination of the $D-$dimensional sufficient statistic $T^i(X) = [T^i_1(X), \cdots, T^i_D(X)]$ leads to a constant. This is a mild assumption on the form of the exponential family.

Assume that $\text{int}(\Theta) \neq \emptyset$. Assume that the likelihood $P_\theta^i$  given by Equation (<ref>)  and prior $\Pi$ given by Equation (<ref>) belong to exponential families with full rank. Then $P_t^j$ belongs to an exponential family of full rank. Moreover, the gradient of the log-partition function $\nabla B_t^j$ is invertible on the interior of $\Theta$.

Lemma <ref> establishes the conditions under which the distributed Bayes posterior belongs to the full-rank exponential family. Critically, the invertibility of \(\nabla B_t^j\) is instrumental for constructing a sequence of consistent M-estimators for \(\theta_0\). These estimators serve as the centering sequence for the Laplace approximation, which is key in proving our Bernstein–von Mises result.

Let $\hat \theta_t^j$ be the M-estimators corresponding to $f_t^j$, as defined in Equation (<ref>).

\begin{equation*}
\hat \theta_t^j = \argmax_{\theta \in \R^p}\langle \theta, \chi_t^j \rangle - B_t^j(\theta).
\end{equation*}

The closed-form expression for $\hat \theta_t^j$ can be obtained as follows:

\begin{equation}\label{eqn-EF-5}
\hat \theta_t^j = (\nabla_\theta B_t^j)^{-1}(\chi_t^j),
\end{equation}

where $\chi_t^j, B_t^j$ are defined in Equation (<ref>).

We now state the asymptotic property of the distributed Bayes posteriors for exponential family distributions.

Let $\Theta$ be an open subset of $\R^p$. Assume that there exists $\theta_0 \in \Theta$ such that $P_0 = P_{\theta_0}^j$ for every $j \in [m]$. Moreover, let Assumptions <ref>, <ref>(a) hold. Let the sequence of estimators $\hat \theta_t^j$ be defined in Equation (<ref>). Let $q_t^j$ be the density of $\sqrt{t}(\theta -\hat \theta_t^j)$ when $\theta \sim P_t^j$ Then,

\begin{equation*}
\int_\Theta |q_t^j(x) - N(0, V_{\theta_0}^{-1})| dx \pto 0,
\end{equation*}

where $V_{\theta_0}$ is the average of the covariance of $T^i$ evaluated at $\theta_0$, i.e. $V_{\theta_0} = \frac{1}{m} \sum_{i = 1}^m \text{Cov}(T^i)$.

The BvM result suggests a sample-based Laplace approximation to the distributed Bayes posterior that belongs to the canonical exponential family. One example of such an approximation is the normal distribution with mean given by the moment estimator $\hat{\theta}_{t}^j$ and covariance given by as$\left(\frac{1}{m} \sum_{i=1}^{m} \hat V(T^{i})\right)^{-1}$, where $\hat V(T^{i})$ is the bootstrapped sample covariance of $T^{i}$. By Proposition <ref> and Slusky's theorem, we obtain that the total variational distance between $P_{t}^j$ and the normal distribution $N\left(\hat{\theta}_{t}^j, \left(\frac{1}{m} \sum_{i=1}^{m} \hat V(T^{i})\right)^{-1}\right)$ asymptotically converges to zero.

One advantage of the Laplace approximation is that it directly calculates an asymptotically valid credible region.  Let $\chi_{\alpha, p}^2$ be the critical value from a $\chi^2$ distribution with $p$ degrees of freedom.  If $\theta \sim N\left(\hat{\theta}_{t}^j, V_{\theta_0}^{-1}\right)$, a credible region for $\theta$ at level $1 - \alpha$ is given by the set of $\theta$ that satisfies

\begin{equation*}
(\hat{\theta}_{t}^j - \theta)^T \left(\frac{1}{m} \sum_{i=1}^{m} \hat V(T^{i})\right)  (\hat{\theta}_{t}^j - \theta)^T \leq \frac{\chi_{\alpha, p}^2}{t},
\end{equation*}

where $\hat V(T^{i})$ can be replaced by any consistent estimator for $\text{Cov}(T^i)$. By Proposition <ref> and Slusky's theorem, the credible region has asymptotic coverage of $1 - \alpha$ under $P_t^j$; hence it is asymptotic valid.

§.§ Distributed Logistic Regression

We consider i.i.d. observations $D_k^j = (X_k^j, Y_k^j)$ for $k \in [t] $ and $j \in [m]$ where covariates $X_k^j \in \X$ for $\X \subseteq \R^p$ and responses $Y_k^j \in \{0, 1\}$. Let $\theta_0 \in \R^p$ be the true and unknown parameter. A logistic regression model generates the data.

\begin{equation}\label{eqn-logistic-1}
Y_k^j \sim \text{Ber}(\nu_k^j), \quad \log \left( \frac{\nu_k^j}{1 - \nu_k^j} \right) = \theta_0^T X_k^j.
\end{equation}

For a logistic regression model with coefficient $\theta$, the conditional distribution $p_\theta(y_k^j \mid x_k^j)$ after marginalizing out $\nu_k^j$ is expressed as:

\begin{equation}\label{eqn-logistic-3}
p_\theta(y_k^j \mid x_k^j) = \exp\left(\langle \theta, x_k^j y_k^j \rangle - \sigma(\langle \theta, x_k^j \rangle)\right),
\end{equation}

where $\sigma(\eta) = \log(1 + e^\eta)$. The function $\sigma(\eta)$ is strictly convex, since $\sigma^{''}(\eta) = \frac{e^\eta}{(1 + e^\eta)^2} > 0$.

Conditioning on the covariates, the model (<ref>) aligns with the canonical form of the exponential family. This allows us to apply the results from Section <ref>.
Although the conditional distribution $p_\theta(y_k^j \mid x_k^j)$ is not a member of the exponential family in $\theta$, we can reparametrize it to fit the canonical form of the exponential family by introducing $\eta_k^j = \theta^T x_k^j$. Under this reparametrization, the distribution takes the form:

\begin{equation}
\label{eqn-logistic-EF-1}
p_{\eta_t^j}(y_k^j) = \exp\left(\eta_k^j y_k^j - \sigma(\eta_k^j)\right),
\end{equation}

For each agent $j$, we let the statistical model $\P_\Theta^j$ be a logistic regression model (<ref>) with prior $\pi(\theta)$ supported on $\R^p$.   By aggregating the log-likelihood, we obtain the closed-form expression for the distributed Bayes posterior:

\begin{equation} \label{eqn-logistic-3.5}
p_t^j (\theta) = \exp\left(\langle \theta, T_t^j \rangle - B_t^j(\theta) - \log \pi(\theta)\right),
\end{equation}

\begin{align*}
T_t^j =   \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] X_k^i Y_k^i,  \quad B_t^j(\theta) = \sum_{i = 1}^m \sum_{k = 1}^t [A^{t-k}_{ji}] \sigma(\langle \theta,  X_k^j \rangle).
\end{align*}

It is difficult to sample from the distribution (<ref>) due to the complex form of $\sigma(\theta)$. As a result, the sample-based Laplace approximation characterized by the Bernstein-von Mises theorem becomes especially helpful. Now, we proceed to state this result.

Let $\hat \theta_t^j$ be the M-estimators corresponding to $f_t^j$:

\begin{equation}\label{eqn-logistic-4}
\hat \theta_t^j = \argmax_{\theta \in \R^p}\langle \theta, T_t^j \rangle - B_t^j(\theta).
\end{equation}

The objective function in (<ref>) is strictly concave, guaranteeing the uniqueness of $\hat{\theta}_t^j$. This also enables efficient numerical methods, such as Newton-Raphson, for finding the value of $\hat{\theta}_t^j$.

Let $\Theta = \R^p$.  Let $D_k^j = (X_k^j, Y_k^j), k \in [t], j \in [m]$ be a sequence of i.i.d. paired random variables generated according to the model (<ref>). Let Assumptions <ref>, <ref>(a) hold.  Under these conditions, the estimator sequence $\hat{\theta}_t^j$ defined in Equation (<ref>) converges to $\theta_0$ in $[\bP_0]-$probability.

Define $q_t^j$ as the density of $\sqrt{t}(\theta - \hat{\theta}_t^j)$, where $\theta \sim P_t^j$. If the following assumptions are satisfied:

* $\bP_{\theta_0} X_1^1 X_1^{1^T}$ exists and is finite and nonsingular.

* $\bP_{\theta_0}|X_{k, a}^1 X_{k, b}^1 X_{k, c}^1| < \infty$ holds for all $a, b, c \in [p]$,

then we have

\begin{equation}
\label{eqn-logistic-5}
\int_\Theta \left|q_t^j(x) - N(0, \hat V_{\theta_0}^{-1})\right| \, dx \pto 0.
\end{equation}

The matrix $\hat V_\theta$ is computed as

\begin{equation} \label{eqn-logistic-6}
\hat V_\theta = \frac{1}{m} \sum_{i = 1}^m {X^i}^T W^i(\theta) X^i,
\end{equation}

where $W^i(\theta)$ is the diagonal matrix defined by

\begin{equation*}
W^i(\theta) = \text{diag}\left(
\frac{e^{\sum_{j=0}^{p} \theta_j x_{1j}^i}}{\left(1 + e^{\sum_{j=0}^{p} \theta_j x_{1j}^i}\right)^2},
\ldots,
\frac{e^{\sum_{j=0}^{p} \theta_j x_{tj}^i}}{\left(1 + e^{\sum_{j=0}^{p} \theta_j x_{tj}^i}\right)^2}
\right).
\end{equation*}

The proposition relies on two model-specific assumptions. Assumption i) focuses on the regularity of the private Fisher information. It is a standard prerequisite in Generalized Linear Models (GLMs) to ensure that the model parameter $\theta$ is identifiable (Example 16.8, [Van der Vaart, 2000]). Assumption ii) imposes a more stringent condition requiring a bounded third moment for the covariates. This assumption is used to verify Assumption <ref>(f) and is generally reasonable in most applications.

Let $\chi_{\alpha, p}^2$ be the critical value from a $\chi^2$ distribution with $p$ degrees of freedom. Proposition <ref> specifies an asymptotically valid credible region for the parameter $\theta$ under the probability measure $P_t^j$ at level $1 - \alpha$. The region is given by the set of $\theta$ that satisfies

\begin{equation*}
(\hat{\theta}_{t}^j - \theta)^T \hat V_{\theta_0} (\hat{\theta}_{t}^j - \theta)^T \leq \frac{\chi_{\alpha, p}^2}{t},
\end{equation*}

where $\hat{\theta}_{t}^j$ is the estimator in Equation (<ref>) and $\hat V_{\theta_0}$ is defined in Equation (<ref>).

Proposition <ref> strengthens the robustness of distributed logistic regression models. The averaging Fisher information in the approximate covariance ensures that extreme covariate values do not overly influence the approximate uncertainty. The model-specific assumptions i) and ii) also serve as robustness checks for when the distributed logistic regression model attains the desired asymptotic property. In practice, checking the assumptions and using Laplace approximation provides an efficient, reliable approach to statistical inference in distributed logistic regression models, making it more resilient to various data irregularities.

§.§ Distributed Detection

In this example, we extend the distributed source location scenario from [Nedić et al., 2017] to demonstrate our theoretical results in a real-world context. Consider a communication network spread over a unit square $[0,1] \times [0,1]$, comprised of $m$ sensor agents located at points $Z^j, j \in [m]$. The agents are interconnected via a connected undirected graph $G$, represented by the adjacency matrix $A$. The goal is to locate a target at an unknown position $\theta_0 \in (0,1) \times (0,1)$. This target generates data $X_t^1, \ldots, X_t^m \sim \delta_{\theta_0}$ at time $t$, as seen by the agents. However, agents do not receive the data $X_t^j$ directly; instead, the $j^{th}$ agent receives a noisy version of $|X_t^j - Z^j|$, which is the Euclidean distance between the data point $X_t^j$ and the agent's location $Z^j$.

The goal is for the agents to identify the unknown location parameter, $\theta_0$ collectively. Due to noise and data corruption during transmission, the $j^{th}$ agent employs a statistical model to describe the observed signal $|X_t^j - Z^j|$. Specifically, this signal is modeled as a  normal distribution $N(|\theta - Z^j|, \sigma^{j^2}),\sigma^j > 0$ truncated to the interval $[0, |Z^j| + \frac{1}{2}]$.

Let the parameter space be $\Theta = (0,1)^2$. For each $\theta \in \Theta$, the $j^{th}$ agent's statistical model is given by:

\begin{equation} \label{eqn-detection-1}
p_\theta^j(X_k^j) = \frac{\phi\left(\frac{|X_t^j - Z^j| - |\theta - Z^j|}{\sigma^j}\right)}{\sigma^j \left[\Phi\left(\frac{|Z^j| + \frac{1}{2} - |\theta - Z^j|}{\sigma^j}\right) - \Phi\left(\frac{- |\theta - Z^j|}{\sigma^j}\right)\right]}  I(0 \leq |X_t^j - Z^j| \leq |Z^j| + \frac{1}{2}).
\end{equation}

Before collecting any data, the agents collectively decide on a uniform prior for the unknown location parameter $\theta$, defined as:

\begin{equation}\label{eqn-detection-2}
\pi(\theta) = I(\theta \in [0,1]^2).
\end{equation}

The statistical models from Equation (<ref>) are based on truncated normal distributions with fixed support, which belongs to the exponential family. The sufficient statistic for $|\theta - Z^j|$ is $|X_t^j - Z^j|$.

We can express the generalized likelihood function $f_t^j(\theta)$, up to a constant, as:

\begin{align*}
f_t^j(\theta) &= -\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \left[\frac{(|X_t^j - Z^j| - |\theta - Z^j|)^2}{2 \sigma^{j^2}} \right. \\
& \qquad \left. + \log \left(\Phi\left(\frac{|Z^j| + \frac{1}{2} - |\theta - Z^j|}{\sigma^j}\right) - \Phi\left(\frac{- |\theta - Z^j|}{\sigma^j}\right)\right) \right].
\end{align*}

The M-estimators for the loss function $f_t^j$ are denoted by $\hat{\theta}_t^j$, and they minimize $f_t^j(\theta)$ over the parameter space $[0,1]^2$:

\begin{equation}\label{eqn-detection-3}
\hat{\theta}_t^j = \argmin_{\theta \in [0,1]^2} f_t^j(\theta).
\end{equation}

Due to the continuity and strict concavity of the objective function $f_t^j$, there exists a unique $\hat{\theta}_t^j$, which enables efficient numerical methods like the Newton-Raphson algorithm for its calculation.

Care is needed when $f_t^j$ achieves its minimum at the corners of $[0,1]^2$, namely at the points $(0,0), (0,1), (1,0), (1,1)$. In these scenarios, $\hat{\theta}_t^j$ does not belong to $\Theta$. Including this technicality is primarily to ensure that our definition of $\hat{\theta}_t^j$ satisfies the conditions of the Argmax Theorem (Theorem 3.2.2 in [Vaart and Wellner, 2023]).

Subsequent results confirm that the sequence of estimators $\hat{\theta}_t^j$ is consistent with $\theta_0$, allowing us to disregard the concern about the extreme corner cases with probability one.

Let $\theta_0 \in \Theta$ be given.  Under Assumption <ref>(a), the sequence of estimators $\hat \theta_t^j$ defined in Equation (<ref>) converges to $\theta_0$ in $[\bP_0]-$probability. Moreover,  $\lim_{t \to \infty}\bP_0(\hat \theta_t^j \in \Theta) = 1$.

We can use the Bernstein von - Mises theorem for the distributed location detection problem.

Let the private statistical models be given by Equation (<ref>) and the prior be given by Equation (<ref>). Let Assumption <ref>, <ref>(a) hold.  For the sequence of estimators $\hat \theta_t^j$ defined in Equation (<ref>), define $q_t^j$ as the density of $\sqrt{t}(\theta - \hat{\theta}_t^j)$, where $\theta \sim P_t^j$. Then we have

\begin{equation}
\int_\Theta \left|q_t^j(x) - N(0, V_{\theta_0}^{-1})\right| \, dx \pto 0,
\end{equation}

where $V_{\theta_0}$ is given by

\begin{equation} \label{eqn-detection-4}
V_{\theta_0} = \frac{1}{m} \sum_{j = 1}^m \frac{(\theta_0 - Z^j) (\theta_0 - Z^j)^T}{\sigma^{j^4} |\theta_0 - Z^j|^2} \left[ \frac{\phi(\frac{|Z^j|+ \frac{1}{2} - |\theta_0 - Z^j|}{\sigma^j}) - \phi(\frac{- |\theta_0 - Z^j|}{\sigma^j})}{ \Phi(\frac{|Z^j|+ \frac{1}{2} - |\theta_0 - Z^j|}{\sigma^j}) - \Phi(\frac{- |\theta_0 - Z^j|}{\sigma^j}) } \right]^2.
\end{equation}

Let $\chi_{\alpha, p}^2$ denote the critical value from a chi-squared distribution with $p$ degrees of freedom. Proposition <ref> delineates an asymptotically valid credible region for the parameter $\theta$ under the probability measure $P_t^j$ at a confidence level of $1 - \alpha$. Specifically, this credible region encompasses all values of $\theta$ that satisfy the following inequality:

\begin{equation*}
(\hat{\theta}_{t}^j - \theta)^\top \hat{V}_{\theta_0} (\hat{\theta}_{t}^j - \theta) \leq \frac{\chi_{\alpha, p}^2}{t},
\end{equation*}

where $\hat{\theta}_{t}^j$ represents the estimator defined in Equation (<ref>). Furthermore, $\hat{V}_{\theta_0}$ serves as the empirical analogue of the covariance matrix $V_{\theta_0}$, as defined in Equation (<ref>), with the true parameter $\theta_0$ replaced by its estimate $\hat{\theta}_{t}^j$.

§ DISCUSSION AND FUTURE DIRECTIONS

We have studied the Frequentist statistical properties of the distributed (Non-Bayesian) Bayesian inference in terms of “posterior” consistency, Bernstein von - Mises theorems, and “posterior” contraction rates. Our results provide the first rigorous insights into the statistical efficiency of distributed Bayesian inference and its dependence on the design parameters of the underlying communication network, such as the number of agents and the network topology. The promising results offer several avenues for future research.

Future work should investigate the Frequentist statistical properties of the distributed Bayes posterior under more complex time-varying network structures, such as networks with single link failures [Shahrampour et al., 2015] or under other random graph structures such as Erdős–Rényi graphs [Erdős et al., 1960]. Understanding how distributed Bayesian inference adapts to such scenarios will be crucial for applying theoretical insights in unpredictable or adversarial settings.

Another compelling direction is to explore the Frequentist coverage properties of the distributed Bayes posterior. This could include analyses of confidence intervals or hypothesis tests when the underlying model is well-specified and misspecified. Understanding how the distributed Bayes posterior aligns with Frequentist criteria can offer additional validation and robustness checks for the model.

There is also potential for bridging our work to the theory of non-Bayesian social learning in econometrics, such as the results in [Molavi et al., 2018]. Insights from non-Bayesian frameworks like variants of the DeGroot model  [DeGroot, 1974, Acemoglu et al., 2011] could shed light on the convergence properties of distributed Bayesian methods, especially in settings where agents have heterogeneous prior beliefs or are influenced by external signals.

Our work finds resonance in the literature on Frequentist distributed inference that focuses on communication efficiency. Future studies could integrate insights from this body of work, including communication-efficient methods [Jordan et al., 2018] and high-dimensional distributed statistical inference [Battey et al., 2015], to further enhance the efficacy and scalability of distributed Bayesian approaches.

Lastly, our theoretical findings could inform the design of communication networks in practical applications. Specifically, an analytical understanding of how the number of agents and communication costs interact could guide the construction of more efficient and robust distributed Bayesian systems.

[Acemoglu et al., 2011]
Daron Acemoglu, Munther A Dahleh, Ilan Lobel, and Asuman Ozdaglar.
Bayesian learning in social networks.
The Review of Economic Studies, 780 (4):0 1201–1236, 2011.

[Battey et al., 2015]
Heather Battey, Jianqing Fan, Han Liu, Junwei Lu, and Ziwei Zhu.
Distributed estimation and inference with statistical guarantees.
arXiv preprint arXiv:1509.05457, 2015.

[Bickel and Kleijn, 2012]
Peter J Bickel and Bas J K Kleijn.
The semiparametric bernstein–von mises theorem.
The Annals of Statistics, 40:0 206–237, 2012.
ISSN 0090-5364.

[Bickel and Yahav, 1969]
Peter J Bickel and Joseph A Yahav.
Some contributions to the asymptotic theory of bayes solutions.
Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 11:0 257–276, 1969.
ISSN 1432-2064.

[Bochkina and Green, 2014]
Natalia A Bochkina and Peter J Green.
The bernstein–von mises theorem and nonregular models.
The Annals of Statistics, 420 (5):0 1850–1878, 2014.

[Borkar and Varaiya, 1982]
Vivek Borkar and Pravin Varaiya.
Asymptotic agreement in distributed estimation.
IEEE transactions on automatic control, 270 (3):0 650–655, 1982.

[Cam et al., 2000]
Lucien Le Cam, Lucien Marie LeCam, and Grace Lo Yang.
Asymptotics in statistics: some basic concepts.
Springer Science & Business Media, 2000.
ISBN 0387950362.

[Castillo and Rousseau, 2015]
Ismaël Castillo and Judith Rousseau.
A bernstein–von mises theorem for smooth functionals in semiparametric models.
The Annals of Statistics, 43:0 2353–2383, 2015.
ISSN 0090-5364.

[Cover, 1999]
Thomas M Cover.
Elements of information theory.
John Wiley & Sons, 1999.

[DeGroot, 1974]
Morris H DeGroot.
Reaching a consensus.
Journal of the American Statistical association, 690 (345):0 118–121, 1974.

[Epstein et al., 2010]
Larry G Epstein, Jawwad Noor, and Alvaro Sandroni.
Non-bayesian learning.
The BE Journal of Theoretical Economics, 10:0 1–20, 2010.

[Erdős et al., 1960]
Paul Erdős, Alfréd Rényi, et al.
On the evolution of random graphs.
Publ. math. inst. hung. acad. sci, 50 (1):0 17–60, 1960.

[Ghosal and Van der Vaart, 2017]
Subhashis Ghosal and Aad Van der Vaart.
Fundamentals of nonparametric Bayesian inference, volume 44.
Cambridge University Press, 2017.

[Ghosh and Ramamoorthi, 2003]
Jayanta K. Ghosh and R.V. Ramamoorthi.
Bayesian Nonparametrics.
Springer-Verlag, 2003.
URL <https://doi.org/10.1007%2Fb97842>.

[Gilardoni and Clayton, 1993]
Gustavo L Gilardoni and Murray K Clayton.
On reaching a consensus using degroot's iterative pooling.
The Annals of Statistics, 210 (1):0 391–401, 1993.

[Golub and Sadler, 2017]
Benjamin Golub and Evan Sadler.
Learning in social networks.
Available at SSRN 2919146, 2017.

[Gubner, 1993]
John A Gubner.
Distributed estimation and quantization.
IEEE Transactions on Information Theory, 390 (4):0 1456–1459, 1993.

[Hoffman et al., 2013]
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.
Stochastic variational inference.
Journal of Machine Learning Research, 2013.

[Jadbabaie et al., 2012]
Ali Jadbabaie, Pooya Molavi, Alvaro Sandroni, and Alireza Tahbaz-Salehi.
Non-bayesian social learning.
Games and Economic Behavior, 76:0 210–225, 2012.
ISSN 0899-8256.

[Jordan et al., 2018]
Michael I Jordan, Jason D Lee, and Yun Yang.
Communication-efficient distributed statistical inference.
Journal of the American Statistical Association, 2018.

[Katsevich, 2023]
Anya Katsevich.
Improved scaling with dimension in the bernstein-von mises theorem for two statistical models.
arXiv preprint arXiv:2308.06899, 2023.

[Kleijn and van der Vaart, 2012]
Bas J K Kleijn and Aad W van der Vaart.
The bernstein-von-mises theorem under misspecification.
Electronic Journal of Statistics, 6:0 354–381, 2012.
ISSN 1935-7524.

[Knoblauch et al., 2022]
Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas.
An optimization-centric view on bayes’ rule: Reviewing and generalizing variational inference.
Journal of Machine Learning Research, 230 (132):0 1–109, 2022.

[Lalitha et al., 2014]
Anusha Lalitha, Anand Sarwate, and Tara Javidi.
Social learning and distributed hypothesis testing.
In 2014 IEEE International Symposium on Information Theory, pages 551–555. IEEE, 2014.

[Laplace, 1809]
Pierre-Simon Laplace.
Supplément au mémoire sur les approximations des formules qui sont fonction de très grands nombres.
Mémoires de l’Académie Royale des sciences de Paris, 1810:0 559–565, 1809.

[LeCam, 1953]
Lucien LeCam.
On some asymptotic properties of maximum likelihood estimates and related bayes estimates.
Univ. California Pub. Statist., 1:0 277–330, 1953.

[Li et al., 2017]
Cheng Li, Sanvesh Srivastava, and David B Dunson.
Simple, scalable and accurate posterior interval estimation.
Biometrika, 104:0 665–680, 2017.
ISSN 0006-3444.

[Medina et al., 2022]
Marco Avella Medina, José Luis Montiel Olea, Cynthia Rush, and Amilcar Velez.
On the robustness to misspecification of $\alpha$-posteriors and their variational approximations.
Journal of Machine Learning Research, 230 (147):0 1–51, 2022.

[Miller, 2021]
Jeffrey W Miller.
Asymptotic normality, concentration, and coverage of generalized posteriors.
J. Mach. Learn. Res., 22:0 168–1, 2021.

[Minsker et al., 2014]
Stanislav Minsker, Sanvesh Srivastava, Lizhen Lin, and David Dunson.
Scalable and robust bayesian inference via the median posterior.
In International conference on machine learning, pages 1656–1664. PMLR, 2014.

[Minsker et al., 2017]
Stanislav Minsker, Sanvesh Srivastava, Lizhen Lin, and David B Dunson.
Robust and scalable bayes via a median of subset posterior measures.
The Journal of Machine Learning Research, 18:0 4488–4527, 2017.
ISSN 1532-4435.

[Molavi et al., 2018]
Pooya Molavi, Alireza Tahbaz-Salehi, and Ali Jadbabaie.
A theory of non-bayesian social learning.
Econometrica, 860 (2):0 445–490, 2018a.

[Molavi et al., 2018]
Pooya Molavi, Alireza Tahbaz‐Salehi, and Ali Jadbabaie.
A theory of non‐bayesian social learning.
Econometrica, 86:0 445–490, 2018b.
ISSN 0012-9682.

[Nedić et al., 2017]
Angelia Nedić, Alex Olshevsky, and César A Uribe.
Fast convergence rates for distributed non-bayesian learning.
IEEE Transactions on Automatic Control, 62:0 5538–5553, 2017.
ISSN 0018-9286.

[Neiswanger et al., 2013]
Willie Neiswanger, Chong Wang, and Eric Xing.
Asymptotically exact, embarrassingly parallel mcmc.
arXiv preprint arXiv:1311.4780, 2013.

[Panov and Spokoiny, 2015]
Maxim Panov and Vladimir Spokoiny.
Finite sample bernstein–von mises theorem for semiparametric problems.
Bayesian Analysis, 100 (3):0 665–710, 2015.

[Qipeng et al., 2015]
Liu Qipeng, Zhao Jiuhua, and Wang Xiaofan.
Distributed detection via bayesian updates and consensus.
In 2015 34th Chinese Control Conference (CCC), pages 6992–6997. IEEE, 2015.

[Rabinovich et al., 2015]
Maxim Rabinovich, Elaine Angelino, and Michael I Jordan.
Variational consensus monte carlo.
Advances in Neural Information Processing Systems, 28, 2015.

[Rahimian et al., 2015]
Mohammad Amin Rahimian, Shahin Shahrampour, and Ali Jadbabaie.
Learning without recall by random walks on directed graphs.
In 2015 54th IEEE Conference on Decision and Control (CDC), pages 5538–5543. IEEE, 2015.

[Rosenthal, 1995]
Jeffrey S Rosenthal.
Convergence rates for markov chains.
Siam Review, 370 (3):0 387–405, 1995.

[Schwartz, 1965]
Lorraine Schwartz.
On bayes procedures.
Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 4:0 10–26, 1965.
ISSN 1432-2064.

[Scott et al., 2016]
Steven L Scott, Alexander W Blocker, Fernando V Bonassi, Hugh A Chipman, Edward I George, and Robert E McCulloch.
Bayes and big data: The consensus monte carlo algorithm.
International Journal of Management Science and Engineering Management, 11:0 78–88, 2016.
ISSN 1750-9653.

[Shahrampour et al., 2015]
Shahin Shahrampour, Alexander Rakhlin, and Ali Jadbabaie.
Distributed detection: Finite-time analysis and impact of network topology.
IEEE Transactions on Automatic Control, 61:0 3256–3268, 2015.
ISSN 0018-9286.

[Shen, 2002]
Xiaotong Shen.
Asymptotic normality of semiparametric and nonparametric posterior distributions.
Journal of the American Statistical Association, 97:0 222–235, 2002.
ISSN 0162-1459.

[Tsitsiklis and Athans, 1984]
John Tsitsiklis and Michael Athans.
Convergence and asymptotic agreement in distributed decision problems.
IEEE Transactions on Automatic Control, 290 (1):0 42–50, 1984.

[Uribe et al., 2022]
Cesar A Uribe, Alexander Olshevsky, and Angelia Nedich.
Non-asymptotic concentration rates in cooperative learning part i: Variational non-bayesian social learning.
IEEE Transactions on Control of Network Systems, 2022a.
ISSN 2325-5870.

[Uribe et al., 2022]
Cesar A Uribe, Alexander Olshevsky, and Angelia Nedich.
Non-asymptotic concentration rates in cooperative learning part ii: Inference on compact hypothesis sets.
IEEE Transactions on Control of Network Systems, 2022b.
ISSN 2325-5870.

[Vaart and Wellner, 2023]
AW van der Vaart and Jon A Wellner.
Empirical processes.
In Weak Convergence and Empirical Processes: With Applications to Statistics, pages 127–384. Springer, 2023.

[Van der Vaart, 2000]
Aad W. Van der Vaart.
Asymptotic Statistics, volume 3.
Cambridge University Press, 2000.
ISBN 0521784506.

[Walker, 2006]
Stephen G Walker.
Bayesian inference via a minimization rule.
Sankhyā: The Indian Journal of Statistics, pages 542–553, 2006.
ISSN 0972-7671.

[Wang and Dunson, 2013]
Xiangyu Wang and David B Dunson.
Parallelizing mcmc via weierstrass sampler.
arXiv preprint arXiv:1312.4605, 2013.

[Wang et al., 2015]
Xiangyu Wang, Fangjian Guo, Katherine A Heller, and David B Dunson.
Parallelizing mcmc with random partition trees.
Advances in neural information processing systems, 28, 2015.

[Wang and Blei, 2019]
Yixin Wang and David M Blei.
Frequentist consistency of variational bayes.
Journal of the American Statistical Association, 1140 (527):0 1147–1161, 2019.

§ REVIEW OF CONCEPTS

§.§ Distances and Divergences

This section reviews a class of divergence functions.

Let $\rho > 0$ and $\rho \neq 1$. The $\rho$-Rényi divergence between two probability measures $P_1$ and $P_2$ is defined as

\begin{equation*}
D_\rho(P_1 \parallel P_2) =
\begin{cases}
\frac{1}{\rho - 1} \log \int \left(\frac{dP_1}{dP_2}\right)^{\rho -1} dP_1 & \text{if } P_1 \ll P_2, \\
+\infty & \text{otherwise.}
\end{cases}
\end{equation*}

The relations between the Rényi divergence and other divergence functions are summarized below:

* When $\rho \to 1$, the Rényi divergence converges to the Kullback–Leibler divergence, defined as

\begin{equation*}
D_{KL}(P_1 \parallel P_2) =
\begin{cases}
\int \log \left(\frac{dP_1}{dP_2}\right) dP_1 & \text{if } P_1 \ll P_2, \\
+\infty & \text{otherwise.}
\end{cases}
\end{equation*}

* When $\rho = 1/2$, the Rényi divergence is related to the Hellinger distance by

\begin{equation*}
D_{1/2}(P_1 \parallel P_2) = -2 \log \left(1 - H(P_1,P_2)\right),
\end{equation*}

and the Hellinger distance is defined as

\begin{equation*}
H(P_1,P_2) = \sqrt{\frac{1}{2} \int \left(\sqrt{dP_1} - \sqrt{dP_2}\right)^2}.
\end{equation*}

* When $\rho = 2$, the Rényi divergence is related to the $\chi^2$-divergence by

\begin{equation*}
D_2(P_1 \parallel P_2) = \log \left(1 + \chi^2(P_1 \parallel P_2)\right),
\end{equation*}

and the $\chi^2$-divergence is defined as

\begin{equation*}
\chi^2(P_1 \parallel P_2) = \int \left(\frac{(dP_1)^2}{dP_2} - 1\right).
\end{equation*}

The total variation distance between two probability measures $P_1$ and $P_2$ is defined as

\begin{equation*}
TV(P_1, P_2) = \frac{1}{2} \int |dP_1 - dP_2|.
\end{equation*}

§.§ An (incomplete) Summary of Previous BvM Results

The history of the Bernstein von Mises (BvM) theorem is marked by steady evolution, expanding its reach to wider contexts with weaker assumptions.  The origins of the Berstein von-Mises (BvM) theorem trace back to [Laplace, 1809].  Earlier works on BvM theorem are restricted to well-specified, i.i.d statistical models with fixed, finite-dimensional parameters and use assumptions that involve up to fourth-order derivatives of the log-likelihood [LeCam, 1953, Bickel and Yahav, 1969, Cam et al., 2000]).

This classical Bernstein von - Mises theorem for i.i.d., finite-dimensional data was brilliantly summarized in <cit.>. The Theorem relies on three assumptions: the existence of a sequence of consistent estimators, local asymptotic normality (LAN), and uniform consistent testing. One typical sufficient assumption for LAN is differentiable in quadratic means (DQM), which involves only first-order derivatives of the log-likelihood to achieve a desirable quadratic expansion. In addition to the LAN assumption, Schwartz [Schwartz, 1965] proposed the concept of uniformly consistent test assumptions around the true parameter $\theta_0$.

The contemporary wave of research has extended BvM theorems beyond the canonical settings. These studies have broached areas such as model misspecification [Kleijn and van der Vaart, 2012], semi-parametric models [Shen, 2002, Bickel and Kleijn, 2012, Panov and Spokoiny, 2015, Castillo and Rousseau, 2015], and parameters positioned on the boundary of the parameter space [Bochkina and Green, 2014]. Among these explorations, a particularly relevant line of inquiry has focused on BvM for non-standard Bayes procedures [Knoblauch et al., 2022], including generalized posteriors [Miller, 2021], Bayes rules derived from optimization perspectives, and variational Bayes posteriors [Wang and Blei, 2019, Medina et al., 2022]. Another emerging thread of work establishes BvM theorems for high dimensional models, with most recent results allowing the dimension to grow at the order $d^2 \lesssim n$ [Panov and Spokoiny, 2015, Katsevich, 2023].

§ PROOFS

§.§ Proofs of Results in Section <ref>

Define vector - valued functions $\phi_t(\theta) = \left[f_t^1(\theta), \cdots, f_t^m(\theta) \right]^T$ and $\log \bp_\theta(X_t) = \left[\log \bP_\theta^1(X_t^1), \cdots, \log \bP_\theta^m (X_t^m)\right]^T$.

Since $\bP_0 |\log \bP_\theta^i| = \bP_0 (\log \bP_\theta^i)_{+} + \bP_0 (\log \bP_\theta^i)_{-} < \infty$, we have

\begin{align*}
f_t^j(\theta) &= -\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \log \bp_\theta^i(X_t^i) = \frac{1}{t} \sum_{k = 1}^t \left< [A^{t-k}_{j,.}] ,\log \bP_\theta\right>. \\
\mbox{ and } \phi_t(\theta) &= \frac{1}{t} \sum_{k = 1}^t A^{t-k} \log \bp_\theta(X_k)
\end{align*}

Using the property of stochastic matrices [Rosenthal, 1995], the diagonalizable matrix $A$ satisfies

\begin{equation*}
\|e_i^T A^t - \frac{1}{m} \ones \|_1 \leq m \delta^t
\end{equation*}

for all $i \in [m]$ and some $\delta < 1$.

This implies hat

\begin{align*}
&\left|\phi_t(\theta)  - \sum_{k = 1}^t \frac{1}{m} \ones \ones^T \log \bp_\theta(X_k) \right| = \left|\frac{1}{t} \sum_{k = 1}^t A^{t-k} \log \bp_\theta(X_k) - \frac{1}{m} \ones \ones^T \log \bp_\theta(X_k) \right| \\
&\leq \frac{1}{t} \sum_{k = 1}^t \left|A^{t-k} -\frac{1}{m} \ones \ones^T \right| |\log \bp_\theta(X_k) | \\
&\leq \frac{1}{t} \sum_{k = 1}^t \max_{i \in [m]} \left|e_i^T A^{t-k} - \frac{1}{m} \ones \right| \max_{i \in [m]}|\log p^i_\theta(X_k^i) | \\
&\leq \frac{1}{t} \sum_{k = 1}^t  m \delta^{t-k} \max_{i \in [m]}|\log p^i_\theta(X_k^i) | \\
&\leq \frac{m}{t} \sum_{k = 1}^t  \sum_{i = 1}^m  \delta^{t-k} |\log p^i_\theta(X_k^i)| \\
&\leq \frac{m}{t} \sum_{k = t - c + 1}^t  \sum_{i = 1}^m  \delta^{t-k} |\log p^i_\theta(X_k^i)|  +  \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m   |\log p^i_\theta(X_k^i)|, \qquad \text{for fixed }c\\
\end{align*}

By the strong law of large number,  $\left|\sum_{k = 1}^t \frac{1}{m} \ones \ones^T \log \bP_\theta^k - \frac{1}{m}  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i \right| = o_d(1)$.  Then

\begin{align*}
\left|\phi_t(\theta)  - \frac{1}{m}  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i \right|  &\leq  \left|\phi_t(\theta)  - \sum_{k = 1}^t \frac{1}{m} \ones \ones^T \log \bp_\theta(X_k) \right|  + \left|\sum_{k = 1}^t \frac{1}{m} \ones \ones^T \log \bp_\theta(X_k)  - \frac{1}{m}  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i \right| \\
&=   \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m   |\log p^i_\theta(X_k^i)| +\frac{m}{t} \sum_{k = t - c + 1}^t  \sum_{i = 1}^m  \delta^{t-k} |\log p^i_\theta(X_k^i)| +o_d(1) \\
&\leq   \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m   |\log p^i_\theta(X_k^i)| +o_d(1) \\
&\leq   m \delta^c  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i +o_d(1)
\end{align*}

Since this for arbitrarily large $c$, we conclude that $ \left|\phi_t(\theta)  - \frac{1}{m}  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i \right| = o_d(1)$, i.e.
By Lemma <ref> [Distributed law of large numbers],

\begin{equation*}
\phi_t(\theta) \pto \frac{1}{m}  \sum_{i = 1}^m \bP_0 \log \bP_\theta^i =  f(\theta) \ones.
\end{equation*}

which implies the coordinate-wise convergence $f_t^j(\theta) \pto f(\theta)$ for every $j \in [m]$.

Let $\epsilon > 0$ be given. If we define $\mu_t^j(B) = \int_B \exp(-t f_t^j(\theta)) \Pi(d\theta) $, then $P_t^j(B) = \frac{\mu_t^j(B)}{\mu_t^j(\Theta)}$ and $\mu_t^j(\Theta)  = z_t^j < \infty$.

In the proof, we consider the $\Theta = \Theta \backslash U_{null}$ where $\pi(\theta) = 0$ for every $\theta \in U_{null}$. This null set should be of little importance because the set has zero weight under $P_t^j$.

We now show that there exists $\alpha > 0$ such that $\inf_{\theta \in \U_\epsilon^c} f_t^j (\theta) \geq  f(\theta_0) + \alpha$ for all $t$ large enough.  By lemma <ref>, we get that there exists $T$ such that $f_t^j (\theta) - f (\theta) < \frac{\epsilon}{2}$ for $t > T$. Given that $\inf_{\theta \in \U_\epsilon^c} f_t^j (\theta) -f(\theta_0) = \inf_{\theta \in \U_\epsilon^c} [f_t^j (\theta) - f (\theta)] + [f(\theta) -f(\theta_0)]$, $\inf_{\theta \in \U_\epsilon^c} f_t^j (\theta) -f(\theta_0) > \frac{\epsilon}{2}$ for $t > T$ because for $\theta \in \U_\epsilon^c, f(\theta) -f(\theta_0) > \epsilon$.

The result shows that for $t > T$, for all $\theta \in \U_\epsilon^c$, $exp(-t(f_t^j(\theta) - f(\theta_0) - \alpha)) \leq 1$. Therefore, for $t > T$,

\begin{equation*}
\exp(t(f(\theta_0) + \alpha)) \mu_t^j (\U_\epsilon^c) = \int_{\U_\epsilon^c} \exp(-t(f_t^j(\theta) - f(\theta_0) - \alpha)) \Pi(d\theta) \leq \int_{\U_\epsilon^c} \Pi(d\theta) \leq 1
\end{equation*}

For any $\theta \in A_{\alpha/2}$, $\lim_{t \to \infty} f_t^j(\theta) - f(\theta_0) - \alpha = f(\theta) - f(\theta_0) - \alpha < -\alpha /2 < 0$. Thus, $exp(-t(f_t^j(\theta) - f(\theta_0) - \alpha)) \to \infty$ as $t \to \infty$. By Fatou's lemma, since $\Pi(A_{\alpha/2}) > 0$ by assumption 2,

\begin{equation*}
\liminf_{t} e^{t(f(\theta_0) + \alpha)} \mu_t^j(A_{\alpha/2}) =  \liminf_{t} \int_{A_{\alpha/2}} \exp(-t(f_t^j(\theta) - f(\theta_0) - \alpha)) \Pi(d\theta) = \infty.
\end{equation*}

Since $\mu_t^j(\Theta) \geq \mu_t^j(A_{\alpha/2})$, we obtain $e^{t(f(\theta_0) + \epsilon)} \mu_t^j(\Theta) \to \infty$.

Combining the two results,

\begin{equation*}
1 - P_t^j(\U_\epsilon) = P_t^j(\U_\epsilon^c) =  \frac{\mu_t^j(\U_\epsilon^c)}{\mu_t^j(\Theta)} = \frac{\exp(t(f(\theta_0) + \alpha))\mu_t^j(\U_\epsilon^c)}{\exp(t(f(\theta_0) + \alpha))\mu_t^j(\Theta)} \to 0,
\end{equation*}

as $t \to \infty$.

Let $j$ and $\theta$ be fixed. If $\bP_0 \ll \bP_\theta^j$, then $D_{KL}(\bP_0 \parallel \bP_\theta^j) < \infty$. This implies that $\bP_0 \log \bp_\theta^j > - \infty$. Likewise, the expectation of $\log \bP_0$ under $\bP_0$ is also finite, i.e., $\bP_0 \log \bP_0 < \infty$.

Since the KL divergence is non-negative, $D_{KL}(\bP_0 \parallel \bP_\theta^j) \geq 0$, it follows that the expectation of $\log \bp_\theta^j$ under $\bP_0$ is less than or equal to the expectation of $\log \bP_0$ under $\bP_0$. Therefore, $\bP_0 \log \bp_\theta^j \leq \bP_0 \log \bp_0 < \infty$, and this completes the proof.

§.§ Proofs of Results in Section <ref>

Suppose $\U$ is a neighborhood of $\theta_0$, within which $\theta \mapsto \log \bP_\theta^i$ is convex for all $i \in [m]$. Being a linear combination of convex mappings, $f_t^j(\theta)$ retains convexity in $\U$. As a convex function, $f_t^j$ is Lebesgue - almost surely differentiable on $\U$. According to Lemma <ref>, $f_t^j \pto f$, and since $f$ is also convex, it remains differentiable a.s. on $\U$.  Given that $\nabla f_t^j(\theta)$ is non-decreasing a.s. on $\U$ for every coordinate, $\nabla f_t^j(\theta_0 - \epsilon) < - \eta$ and $\hat \theta_t^j \leq \theta_0 - \epsilon$ implies $\nabla f_t^j(\hat \theta_t^j) < - \eta$, which has probability tending to $0$ for every $\eta > 0$ if $f_t^j(\hat \theta_t^j) = o_d(1)$. This shows that for every $\epsilon, \eta > 0$,

\begin{equation*}
P( \|\hat \theta_t^j - \theta_0 \| > \epsilon) + o(1) \geq  P( \nabla f_t^j(\theta_0 - \epsilon) < -\eta,\nabla f_t^j(\theta_0 + \epsilon) > \eta )
\end{equation*}

Since $\nabla f(\theta_0 - \epsilon) < 0 < \nabla f(\theta_0 + \epsilon)$, taking $2 \eta$ to be the smallest coordinates of $\nabla f(\theta_0 - \epsilon)$ and $\nabla f(\theta_0 + \epsilon)$ makes the right hand side converge to $1$.

Using the consistency of $\hat \theta_t^j$ along with Lemma <ref>, given the scaling factor $\epsilon_t^j = \frac{1}{\sqrt{t}} $, we have

\begin{equation*}
f_t^j(\hat\theta_t^j) = f_t^j(\theta_0) + \frac{1}{\sqrt{t}} (\hat \theta_t^j - \theta_0)^T V_{\theta_0} \Delta_{t, \theta}^j + \frac{1}{2} (\hat \theta_t^j - \theta_0)^T V_{\theta_0} (\hat \theta_t^j - \theta_0) + o_d(1) = f_t^j(\theta_0) + o_d(1),
\end{equation*}

by applications of Slusky's theorem.

By Lemma <ref>, $f_t^j(\theta_0)$ converges to $f(\theta_0)$ in probability so $f_t^j(\hat\theta_t^j) = f(\theta_0) + o_d(1)$.

Assumption <ref>(a) assures the existence of $\delta > 0$ such that:

\begin{align*}
\inf_{\theta \in B_\epsilon(\hat\theta_t^j)^c} (f_t^j(\theta) - f_t^j(\hat\theta_t^j)) &= \inf_{\theta \in B_{2\epsilon}(\theta_0)^c} (f_t^j(\theta) - f_t^j(\hat\theta_t^j)) \\
&=  \inf_{\theta \in B_{2\epsilon}(\theta_0)^c} (f_t^j(\theta) - f(\theta_0)) + o_d(1) \\
& > 2\delta + o_d(1) \overset{\text{$t$ large}}{>} \delta.
\end{align*}

hence $\inf_{\theta \in B_\epsilon(\hat\theta_t^j)^c} (f_t^j(\theta) - f_t^j(\hat\theta_t^j)) > \delta$ with probability tending to $1$ under $[\bP_0]$.

For any $\epsilon >0$ and $\delta > 0$,

\begin{equation*}
Q_t(B_\delta(0)) = \Pi_t (B_{\delta/a_t}(\theta_t)) \leq \Pi_t(B_\epsilon(\theta_0)),
\end{equation*}

for all $t$ sufficiently large. Therefore, since $Q_t \to Q$ in total variational distance, we obtain

\begin{equation*}
Q(B_\delta(0)) = \lim_{t \to \infty}   Q_t(B_\delta(0)) \leq \liminf_t \Pi_t(B_\epsilon(\theta_0)).
\end{equation*}

Take $\delta$ arbitrarily large shows that $\lim_{t \to \infty} \Pi_t(B_\epsilon(\theta_0)) = 1$.

For each $j \in [m]$,  the sequence of estimators $\hat \theta_t^j$ are are chosen to satisfy the equation

\begin{equation*}
\sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}] \nabla \log p_{\hat \theta_t^j}^i (X_k^i) = 0.
\end{equation*}

The function $f_t^j(\theta)$ admits the following quadratic expansion.

\begin{align*}
&f_t^j(\theta) =  -\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]\log \bp_\theta^i(X_k^i)  \\
&=  -\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}] \{\log p_{\hat \theta_t^j}^i(X_k^i)  + \nabla \log p_{\hat \theta_t^j}^i(X_k^i)(\theta -\hat \theta_t^j) \\
&- \frac{1}{2}(\theta -\hat \theta_t^j)  \nabla \log p_{\hat \theta_t^j}^i(X_k^i) \left(\nabla \log p_{\hat \theta_t^j}^i(X_k^i)\right)^T (\theta -\hat \theta_t^j)^T  + O_p(|\theta -\hat \theta_t^j|^3)\}\\
&= f_t^j(\hat\theta_t^j) + \frac{1}{2}(\theta -\hat \theta_t^j) \left\{\underbrace{\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]     \nabla \log p_{\hat \theta_t^j}^i(X_k^i)   \left( \nabla\log p_{\hat \theta_t^j}^i(X_k^i)\right)^T}_{\hat V_t^j} \right\}(\theta -\hat \theta_t^j)^T - \underbrace{O_p(|\theta -\hat \theta_t^j|^3)}_{r_t^j(\theta -\hat \theta_t^j)}.
\end{align*}

The sequence of local estimates $\hat \theta_t^j \in \Theta$ satisfies the consistency requirement.  By denoting the remainder term as $r_t^j(\theta -\hat \theta_t^j)$, it immediately holds that

\begin{equation*}
r_t^j(h) = O_p(|h|^3).
\end{equation*}

By Theorem 7.2 of [Van der Vaart, 2000],  Assumption <ref>(c) implies the existence of the nonsingular Fisher information matrix $V_{\theta_0}^i$ for any $i \in [m]$.  The average of Fisher information matrix, defined as $V_{\theta_0} = \frac{1}{m} \sum_{i = 1}^m V_{\theta_0}^i$, is also nonsingular and positive definite. It remains to show that $\hat V_t^j \pto V_{\theta_0}$.

We represent the sequence of matrices $\hat V_t^j$ as:

\begin{align*}
\hat V_t^j &= \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]     \nabla \log p_{\hat \theta_t^j}^i(X_k^i)   \left(\nabla \log p_{\hat \theta_t^j}^i(X_k^i)\right)^T .
\end{align*}

Using Lipschitz continuity of the gradient (Assumption <ref>(d)), we can decompose $\hat V_t^j$ as:

\begin{align*}
\hat V_t^j   &=  \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]     \left( \nabla \log \bp_{\theta_0}^i(X_k^i)  +s^i(X_k^i)(\hat \theta_t^j - \theta_0) \right) \left( \nabla \log \bp_{\theta_0}^i(X_k^i)  +s^i(X_k^i)(\hat \theta_t^j - \theta_0) \right)^T .
\end{align*}

Applying the distributive property, we can decompose this into three terms:

\begin{align*}
\hat V_t^j &=\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  \nabla \log \bp_{\theta_0}^i(X_k^i) \left(\nabla \log \bp_{\theta_0}^i(X_k^i) \right)^T \\
&+ (\hat \theta_t^j - \theta_0) \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  s^i(X_k^i) \nabla \log \bp_{\theta_0}^i(X_k^i)^T  \\
&+ (\hat \theta_t^j - \theta_0)  (\hat \theta_t^j - \theta_0)^T \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  [s^i(X_k^i)]^2.
\end{align*}

By Lemma <ref>, we have

\begin{equation*}
\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  \nabla \log \bp_{\theta_0}^i(X_k^i) \left(\nabla \log \bp_{\theta_0}^i(X_k^i) \right)^T  \pto V_{\theta_0},
\end{equation*}

\begin{equation*}
\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  s^i(X_k^i) \nabla \log \bp_{\theta_0}^i(X_k^i)^T  \pto \bP_0[s^i \nabla \log \bp_{\theta_0}^i],
\end{equation*}

By the Cauchy-Schwarz inequality, we have:

\begin{equation*}
\bP_0[s^i \nabla \log \bp_{\theta_0}^i] \leq \sqrt{\bP_0 [s^i]^2} \sqrt{V_{\theta_0}^i} < \infty,
\end{equation*}

thus the second term is $o(1)$.

\begin{equation*}
\frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]  [s^i(X_k^i)]^2 \pto \bP_0 [s^i(X_k^i)]^2.
\end{equation*}

Finally, since $\hat \theta_t^j - \theta_0 = o_p(1)$, the second and the third term are $o_p(1)$. Combining all these results, we have

\begin{equation*}
\hat V_t^j \pto V_{\theta_0}.
\end{equation*}

Note that $q_t^j(x) = p_t^j(\hat\theta_t^j + x/\sqrt{t}) t^{-p/2}$. Define

\begin{align*}
g_t^j(x) &= \exp(-t[f_t^j(\hat\theta_t^j+ x/\sqrt{t}) - f_t^j(\hat\theta_t^j)]) \pi(\hat\theta_t^j + x/\sqrt{t}) \\
&= \exp(t f_t^j(\hat\theta_t^j)) q_t^j (x) z_t^j t^{p/2}.
\end{align*}

recalling that $z_t^j < \infty$ by assumption and define

\begin{equation*}
g_0(x) = \exp(-\frac{1}{2} x^T V_{\theta_0} x) \pi(\theta_0)
\end{equation*}

Let $\alpha \in (0, \lambda)$, where $\lambda$ is the smallest eigenvalue of $V_{\theta_0}$. Let $\epsilon >0$ small enough that $\epsilon < \alpha/ (2 c_0), \epsilon < \epsilon_0$ and $\pi(\theta) \leq 2 \pi(\theta_0)$ for all $\theta \in B_{2 \epsilon}(\theta_0)$. Let $\alpha = \liminf_t \inf_{\theta \in B_\epsilon(\hat\theta_t^j)^c} (f_t^j(\theta) -f_t^j(\hat\theta_t^j))$, noting that $\delta > 0$ with asymptotic probability $1$ by lemma <ref>. Let $A_t^j = V_t^j - \alpha I$ and $A_0 = V_{\theta_0} - \alpha I$, define

\begin{align*}
h_t^j(x) &= \begin{cases} \exp(-\frac{1}{2} x^T A_t^j x) 2 \pi(\theta_0) & \mbox{if }|x| < \epsilon \sqrt{t} \\
\exp(-\frac{1}{2} t \delta) 2\pi(\hat\theta_t^j + x/\sqrt{t}) & \mbox{if }|x| \geq \epsilon \sqrt{t} \end{cases} \\
h_0(x) &= \exp(-\frac{1}{2}x^T A_0 x) 2\pi(\theta_0)
\end{align*}

We show that

* $g_t^j \to g_0$ and $h_t^j \to h_0$ almost surely

* $\int h_t^j$ converges to $\int h_0$

* $g_t^j = |g_t^j| \leq h_t^j$ with asymptotic probability $1$, and

* $g_t^j, g_0, h_t^j, h_0 \in L^1(dx)$ with asymptotic probability $1$.

By the generalized dominated convergence theorem, these results imply that $\int g_t^j \pto \int g_0$ and $\int |(g_t^j){1/m} - g_0| dx \pto 0$ with $[\bP_0]-$probability. Assuming these results are true, we want to show how asymptotic normality follows. Since $\int q_t^j = 1$, the definition of $g_t^j$ implies that

\begin{equation*}
\exp(t f_t^j(\hat\theta_t^j)) t^{p/2} z_t^j = \int g_t^j \to \int g_0 = \pi(\theta_0) \frac{(2\pi)^{p/2}}{|m V_{\theta_0}|^{1/2}}.
\end{equation*}

where $|V_{\theta_0}| = |det(V_{\theta_0})|$ and therefore

\begin{equation*}
z_t^j \approx \frac{exp(-t f_t^j(\hat\theta_t^j)) \pi(\theta_0)}{|m V_{\theta_0}|^{1/2}} (\frac{2\pi}{t})^{p/2}.
\end{equation*}

as $t \to \infty$.
For any $a_t^j \to a \in \R$, we have that $\int |a_t^j  g_t^j - ag_0| \to 0$ since

\begin{equation*}
\int |a_t^j  g_t^j - a g_0| \leq \int |a_t^j  g_t^j - a_t^j  g_0| + \int |a_t^j  g_0 - a g_0| = |a_t^j | \int |g_t^j - g_0| + |a_t^j  - a| \int |g_0| \to 0
\end{equation*}

Therefore, if we let $a_t^j  = (exp(tf_t^j(\hat\theta_t^j)) t^{p/2} z_t^j)^{-1}$ and $a = (\pi(\theta_0) \frac{(2\pi)^{p/2}}{|V_{\theta_0}|^{1/2}})^{-1}$, we have shown that $a_t^j \to a$ and thus

\begin{equation*}
\int |q_t^j(x) - \frac{|V_{\theta_0}|^{1/2}}{(2\pi)^{p/2}} \exp(-\frac{1}{2}x^T V_{\theta_0} x)| dx \pto 0.
\end{equation*}

This shows asymptotic normality. Finally, the posterior consistency at $\theta_0$ follows from Lemma <ref>. It remains to show statements 1 to 4 above.

1) Fix $x \in \R^p$. For $t$ sufficiently large, $|x| \leq \epsilon \sqrt{t}$. It follows that

\begin{equation*}
h_t^j(x) = \exp(-\frac{1}{2}x^T A_t^j x) 2\pi(\theta) \to \exp(-\frac{1}{2}x^T V_{\theta_0} x) 2 \pi(\theta_0) = h_0(x).
\end{equation*}

since $A_t^j \to A_0$ almost surely. For $g_t^j$, first note that $\pi(\hat\theta_t^j+ x/\sqrt{t}) \to \pi(\theta_0)$ since $\pi$ is continuous at $\theta_0$ and $\theta_t^j\to \theta_0$ and $x/\sqrt{t} \to 0$. By the uniform convergence of $f_t^j$ to $f_t$,

\begin{align*}
t(f_t^j(\hat\theta_t^j+ x/\sqrt{t}) - f_t^j(\hat\theta_t^j))
&= -\frac{1}{2} x^T V_t^j x +  t r_t^j(x/\sqrt{t}) \to -\frac{1}{2} x^T V_{\theta_0} x.
\end{align*}

because $V_t^j \to V_{\theta_0}$ and $|x/\sqrt{t}| < \epsilon_0$ for all $t$ sufficiently large. We then infer a bound on $r_t^j$,

\begin{equation*}
|t r_t^j(x/\sqrt{t})| \leq t c_0 |x/\sqrt{t}|^3 = c_0 |x|^3 / \sqrt{t} \to 0.
\end{equation*}

as $t \to \infty$. Therefore, $g_t^j(x) \to g_0(x)$.

2) By the definition of $h_t^j$, letting $B_t = B_{\epsilon \sqrt{t}}(0)$,

\begin{equation*}
\int h_t^j = \int_{B_t} \exp(-\frac{1}{2} x^T A_t^j x) 2 \pi(\theta_0) dx + \int_{B_t^c}exp(-\frac{1}{2} n \delta) 2\pi(\hat\theta_t^j + x/\sqrt{t}) dx.
\end{equation*}

Since $A_t^j \to A_0$ almost surely and $A_0$ is positive definite, for all $t$ sufficiently large, $A_t^j$ is also positive definite and the first term equals

\begin{equation*}
2\pi(\theta_0) \frac{(2\pi)^{p/2}}{|A_t^j|^{1/2}} P(|(A_t^j)^{-1/2}Z | < \epsilon \sqrt{t}) \asto   2\pi(\theta_0) \frac{(2\pi)^{p/2}}{|A_0|^{1/2}} = \int h_0 .
\end{equation*}

where $Z \sim N(0,1)$. The second integral converges to $0$, since it is nonnegative and has an upper bound as

\begin{equation*}
\int_{\R^p}\exp(-\frac{1}{2} t \delta) 2\pi(\hat\theta_t^j + x/\sqrt{t}) dx = \exp^{-\frac{1}{2} n \delta} n^{p/2} \to 0.
\end{equation*}

using the fact that $\pi(\hat\theta_t^j+ x/\sqrt{t}) t^{-p/2}$ is the density of $\sqrt{t}(\theta -\hat \theta_t^j)$ where $\theta \sim \pi$.

3) With $\bP_0-$probability going to $1$,  $|\theta_t^j - \theta_0| < \epsilon$ thus the bound on $r_t^j$ applies, $\inf_{\theta \in B_\epsilon(\hat\theta_t^j)^c} (f_t^j(\theta) - f_t^j(\hat\theta_t^j)) > \delta$. Let $x \in \R^p$ be given again. If $|x| > \epsilon \sqrt{t}$, then $f_t^j(\hat\theta_t^j+ x/\sqrt{t}) - f_t^j(\hat\theta_t^j) > \delta/2$, and thus,

\begin{equation*}
g_t^j(x) \leq \exp(-\frac{1}{2} t \delta) \pi(\hat\theta_t^j+ x/\sqrt{t}) = h_t^j(x).
\end{equation*}

In the meantime, if $|x| < \epsilon \sqrt{t}$, then $\pi(\hat\theta_t^j+ x/\sqrt{t}) \leq 2 \pi(\theta_0)$ by our choice of $\epsilon$, and

\begin{align*}
t(f_t^j(\hat\theta_t^j + x/\sqrt{t}) - f_t^j(\hat\theta_t^j))
&\geq \frac{1}{2} x^T V_t^j x + t r_t^j(x/\sqrt{t}) \\
&\geq  \frac{1}{2} x^T V_t^j x - \frac{1}{2}\alpha xT x = \frac{1}{2}x^T A_t^j x.
\end{align*}

since $|t r_t^j (x/\sqrt{t})| \leq 1/2 \alpha |x|^2$, by the fact that $|x/\sqrt{t}| < \epsilon < \epsilon_0$ and $\epsilon < \alpha/(2c_0)$. Therefore,

\begin{equation*}
\lim_{t \to \infty}\bP_0\left(g_t^j (x) \leq h_t^j(x)\right) = 1.
\end{equation*}

4) Since $V_{\theta_0}$ and $A_0$ are positive definite, $\int g_0$ and $\int h_0$ are finite. Since $\int h_t^j \to \int h_0$, we have $\lim_{t \to \infty} \bP_0(\int g_t^j \leq \int h_t^j < \infty)  = 1$. The measurability of $g_t^j, h_t^j$ with respect to $\bP_0$ follows from the measurability of $f_t^j$ and $\pi$.