In the proof of Theorem <ref>, the differential in quadratic means assumption (Assumption <ref>(c)) is only used to establish Equation (<ref>). Assumption <ref>(e) states that each $\log p_\theta^i$ is twice continuously differentiable with nonsingular Hessian matrices. Since $f_t^j(\theta)$ is a positive linear combination of $\log p_\theta^i$, $f_t^j(\theta)$ is also twice continuously differentiable with a nonsingular Hessian matrix.  Therefore, Equation (<ref>) holds when we replace Assumption <ref>(c) with Assumption <ref>(e).

Denote $\langle \cdot, \cdot \rangle$ as the tensor inner product, and $\otimes$ as the tensor product. Let Assumption <ref>(f) hold.

We have a third - order Taylor expansion of $f_t^j(\theta)$ around $\hat \theta_t^j$ with nonsingular quadratic terms:

\begin{equation*}
f_t^j(\theta) = f_t^j(\hat \theta_t^j) + \frac{1}{2} (\theta - \theta_t^j)^T \nabla_\theta^2 f_t^j(\hat \theta_t^j) (\theta - \theta_t^j) + \frac{1}{6} \langle \nabla_\theta^3 f_t^j(\tilde \theta), (\theta - \theta_t^j) \otimes (\theta - \theta_t^j) \otimes (\theta - \theta_t^j) \rangle,
\end{equation*}

for some $\tilde \theta$ between $\theta_0$ and $\hat \theta_t^j$.

Under Assumption <ref>(f), we have $|\nabla^3_\theta \log p_\theta^j(x)| \leq \phi^j(x)$ for integrable $\phi$ in a neighborhood $\mathcal N^j$ of $\theta_0$.  Define $\mathcal N = \cap_{j = 1}^m \mathcal N^j$. The consistency assumption implies that $\bP_0(\hat \theta_t^j \in \mathcal N) \to 1$. For all $\theta \in \mathcal N$, we have

\begin{equation*}
\bP_0 |\nabla_\theta^3 f_t^j(\tilde \theta)| \leq \bP_0 \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] |\nabla^3 \log p_\theta^j(X_k^j)| \leq \bP_0 \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \phi^j(X_k^j) < \infty.
\end{equation*}

By the uniform large of large number (Theorem 1.3.3, [Ghosh and Ramamoorthi, 2003]) and Lemma <ref>,  we have

\begin{equation*}
\sup_{\theta \in \mathcal N} \left\|\nabla_\theta^3 f_t^j(\tilde \theta) - \frac{1}{m}  \sum_{i = 1}^m \bP_0 \nabla^3 \log p_\theta^i \right\| \to 0.
\end{equation*}

In the limit, we have $\frac{1}{m} \sum_{i = 1}^m \bP_0 |\nabla^3 \log p_\theta^i| \leq \frac{1}{m} \sum_{i = 1}^m \bP_0 \phi^i$, hence $\nabla_\theta^3 f_t^j(\tilde \theta)$ is asymptotically tight and the cubic term is $o_p(1)$.

We also have that

\begin{equation*}
\nabla_\theta^2 f_t^j(\hat \theta_t^j) \leq  \nabla_\theta^2 f_t^j(\theta_0) + \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] |\langle \phi^j(X_k^j) ,  \hat \theta_t^j - \theta_0 \rangle|,
\end{equation*}

hence $\nabla_\theta^2 f_t^j(\hat \theta_t^j) = \nabla_\theta^2 f_t^j(\theta_0) + o_p(1)$.

Finally, we have $f_t^j(\hat \theta_t^j) = f_t^j(\theta_0) + o_p(1)$ by continuous mapping theorem (Theorem 18.11, [Van der Vaart, 2000]).

Putting the terms together, we conclude that

\begin{equation*}
f_t^j(\theta) = f_t^j(\theta_0) + \frac{1}{2} (\theta - \theta_t^j)^T \nabla_\theta^2 f_t^j(\theta_0) (\theta - \theta_t^j) + o_p(1).
\end{equation*}

Thus, Equation (<ref>) holds we substitute Assumption <ref>(c) and <ref>(d).

By Theorem 7.2 from [Van der Vaart, 2000], Differentiability in quadratic means implies that $\bP_0 \nabla \log \bp_{\theta_0}^j = 0$ and the Fisher information matrix $V_{\theta_0}^j = \bP_0 \nabla \log \bp_{\theta_0}^j (\nabla \log \bp_{\theta_0}^j)^T$ exists and is non-singular.

By the central limit theorem, we have

\begin{equation*}
\frac{1}{\sqrt{t}} \sum_{k = 1}^t \nabla \log \bp_{\theta_0}^j(X_k^j) \dto N(0, V_{\theta_0}^j)
\end{equation*}

Let's define $V_{\theta_0}$ as the average Fisher information matrix, i.e., $V_{\theta_0} = \frac{1}{m} \sum_{j = 1}^m V_{\theta_0}^j$.

If we choose $\epsilon_t^j = \frac{1}{\sqrt{t}}$, then by Taylor expansion,

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

Define $\Delta_{t, \theta_0}^j$ and $V_t$ as:

\begin{align*}
&\Delta_{t, \theta_0}^j =\frac{1}{\sqrt{t}} V_{\theta_0}^{-1} \left(\sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]\nabla \log \bp_{\theta_0}^i(X_k^i)\right), \quad \text{and} \\
& V_t^j = \frac{1}{t}\left( \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}]\nabla \log \bp_{\theta_0}^i(X_k^i) {\nabla \log \bp_{\theta_0}^i(X_k^i)}^T \right).
\end{align*}

Substituting these terms in the Taylor expansion, we have:

\begin{equation*}
t f_t^j(\theta_0 + \epsilon_t^j h)  =  t  f_t^j(\theta_0) - h^T V_{\theta_0} \Delta_{t, \theta_0}^j  + \frac{1}{2} h^T V_t^j h + o_p(|h^2|).
\end{equation*}

By Lemma <ref> (distributed central limit theorem), we have

\begin{equation*}
\Delta_{t, \theta_0}^j \dto N(0, (m V_{\theta_0})^{-1}),
\end{equation*}

thus $\Delta_{t, \theta_0}^j$ is bounded in $[\bP_0]$-probability by Prokhorov's theorem (Theorem 2.4, [Van der Vaart, 2000]).

By Lemma <ref> (distributed law of large number),

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

This establishes the stochastic LAN:

\begin{equation*}
\sup_{h \in K} \left| - t f_t^j(\theta_0 + \epsilon_t^j h) + t f_t^j(\theta_0) -h^T V_{\theta_0} \Delta_{t, \theta_0}^j + \frac{1}{2} h^T V_{\theta_0} h \right| = \sup_{h \in K}o_p(|h|) \to 0.
\end{equation*}

Since $P_t^j\left(A \right)$ is uniformly bounded by $1$ in $L_\infty$, the set $\{P_t^j(A) \}_{t \in \N}$ is uniformly integrable. Then assumption (<ref>) implies that

\begin{equation*}
\bP_0 P_t^j\left(\theta: d(\theta,  \theta_0) \geq M_t \epsilon_t^j \right) \to 0.
\end{equation*}

The rest of the proof is split into two parts: in the first part, we prove the assertion conditional on an arbitrary compact set $K \subset \R^p$, and in the second part, we use this to prove the asymptotic normality convergence. Throughout the proof, we denote the posterior for $h_t^j = (\theta - \theta_0)/\epsilon_t^j$ when $\theta \sim P_t^j$ by $Q_t^j$. For all Borel sets $B$, the posterior for $h_t^j$ follows from that for $\theta$ by

\begin{equation*}
Q_t^j(h_t^j \in B) = P_t^j\left(\theta - \theta_0 \in \epsilon_t^j B \right).
\end{equation*}

The assumption that (<ref>) holds at $\theta_0$ means that for all $j \in [m]$, there exists non-singular $V_{\theta_0}^j$ and asymptotically tight sequence $\Delta_{t, \theta_0}^j$ such that for every compact set $K \subset \R^p$, as $t \to \infty$,

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

We denote the normal distribution $N(\Delta_{t, \theta}^j, V_{\theta_0}^{-1})$ by $\Phi_t^j$. For a compact subset $K \subset \R^p$ such that $Q_t^j(h_t^j \in K) > 0$, we define a conditional measure $Q_{K,t}^j$ of $Q_t^j$ by $Q_{K,t}^j(B) = Q_t^j(B \cap K)/Q_t^j(K)$. Similarly, we define a conditional measure $\Phi_{K,t}^j$ corresponding to $\Phi_t^j$.

Let $K \subset \R^p$ be a compact subset of $\R^p$. For every open neighborhood $U \subset \Theta$ of $\theta_0$, $\theta_0 + K \epsilon_t^j \subset U$ for large enough $t$. Since $\theta_0$ is an interior point of $\Theta$, for large enough $t$, the random function $\psi_t^j: K \times K \mapsto \R$

\begin{equation*}
\psi_t^j(h, h') = \left\lvert 1 - \frac{\phi_t(h')}{\phi_t(h)}
\frac{q_t^j(h')}{q_t^j(h)} \exp\left(t f_t^j(\theta_0 + \epsilon_t^j h) - t f_t^j(\theta_0 + \epsilon_t^j h')\right)\right\rvert_+
\end{equation*}

is well-defined,  where $\phi_t: K \rightarrow \R$ is the Lebesgue density of the distribution $N(\Delta_{t, \theta_0}^j, V_{\theta_0}^{-1})$, $q_t^j: K \rightarrow \R$ is the Lebesgue density of the prior for the centered and rescaled parameter $h_t^j$, and $f_t^j: K \rightarrow \R$ is the random functions defined in (<ref>).

By the distributed LAN assumption (<ref>), we have for every $h \in K$,

\begin{equation*}
t f_t^j(\theta_0 + \epsilon_t^j h) - t f_t^j(\theta_0)= h^T V_{\theta_0} \Delta_{t, \theta_0}^j - \frac{1}{2} h^T V_{\theta_0} h + o_{\bP_0}(1).
\end{equation*}

For any two sequences $h_t, h'_t \in K$, $\lim_{t \to \infty}q_t^j(h_t)/q_t^j(h'_t) =  1$ implies that

\begin{align*}
&\log \left(\frac{\phi_t(h'_t)}{\phi_t(h_t)}  \frac{q_t^j(h_t)}{q_t^j(h'_t)} \right) + t f_t^j(\theta_0 + \epsilon_t^j h_t) - t f_t^j(\theta_0 + \epsilon_t^j h'_t) \\
&= (h'_t - h_t)^T V_{\theta_0} \Delta_{t, \theta_0}^j - \frac{1}{2} h_t^T V_{\theta_0} h_t + \frac{1}{2} {h'_t}^T V_{\theta_0} h'_t + o_{\bP_0}(1)  \\
& + \frac{1}{2} (h_t' - \Delta_{t, \theta_0}^j )^T V_{\theta_0}(h_t' - \Delta_{t, \theta_0}^j )- \frac{1}{2} (h_t - \Delta_{t, \theta_0}^j )^T V_{\theta_0}(h_t - \Delta_{t, \theta_0}^j )
\\
&= o_{\bP_0}(1),
\end{align*}

as $t \to \infty$.

Since all functions $\psi_t^j$ depend continuously on $(h, h')$ and $K \times K$ is compact, we conclude that,

\begin{equation} \label{eqn:psi-uniform-convergence}
\sup_{h,h' \in K} \psi_t^j(h, h') \pto 0, \text{ as } t \to \infty,
\end{equation}

where the convergence is in outer probability.

Assume that $K$ contains a neighborhood of $0$ (to guarantee that $\Phi_t^j(K) > 0$) and let $\Xi_t^j$ denote the event that $Q_t^j(K) > 0$. Let $\eta > 0$ be given, and based on that, define the events:

\begin{equation*}
\Omega_t^j = \left\{ \sup_{h,h' \in K} \psi_t^j(h, h') \leq \eta \right\}.
\end{equation*}

Consider the inequality (recall that the total-variation norm $||\cdot||_{TV}$ is bounded by $2$):

\begin{align}
\bP_0 \| Q_{K,t}^j - \Phi_{K,t}^j \|_{TV} 1_{\Xi_t^j}
&\leq \bP_0 \| Q_{K,t}^j - \Phi_{K,t}^j \|_{TV} 1_{\Omega_t^j \cap \Xi_t^j} + 2\bP_0(\Xi_t^j \setminus \Omega_t^j). \label{eqn:inequality}
\end{align}

As a result of Equation (<ref>), the second term is $o(1)$ as $t \to \infty$. The first term on the r.h.s. is calculated as follows:

\begin{align*}
\frac{1}{2}\bP_0 \| Q_{K,t}^j - \Phi_{K,t}^j \|_{TV} 1_{\Omega_t^j \cap \Xi_t^j}
&= \bP_0 \int \left(1 - \frac{d\Phi_{K,t}^j}{dQ_{K,t}^j}\right)_+ dQ_{K,t}^j 1_{\Omega_t^j \cap \Xi_t^j} \\
&= \bP_0 \int_K \left(1 - \int_K\frac{ \phi_{K,t}^j(h) q_t^j(h') s_t^j(h')}{\phi_{K,t}^j(h') q_t^j(h)s_t^j(h)} d\Phi_{K,t}^j(h') \right)_+  dQ_{K,t}^j(h) 1_{\Omega_t^j \cap \Xi_t^j},
\end{align*}

where $s_t^j: K \rightarrow \R$ is the function defined as $s_t^j(h) = \exp(t f_t^j(\theta_0 + \epsilon_t^j h))$, and $q_t^j: K \rightarrow \R$ is the Lebesgue density of the prior for the centered and rescaled parameter $h_t^j$.

Note that for all $h, h' \in K$, $\phi_{K,t}^j(h)/\phi_{K,t}^j(h') = \phi_t(h)/\phi_t(h')$, since on $K$, $\phi_{K,t}^j$ differs from $\phi_t$ only by a normalization factor. Using Jensen's inequality (with respect to the $\Phi_{K,t}^j$-expectation) for the convex function $x \mapsto (1 - x)_+$, we derive:

\begin{align*}
&\frac{1}{2} \bP_0 \| Q_{K,t}^j - \Phi_{K,t}^j \|_{TV} 1_{\Omega_t^j \cap \Xi_t^j} \\
&\leq \bP_0 \int \left(1 - \frac{s_t^j(h') q_t^j(h') \phi_t(h')}{s_t^j(h') q_t^j(h') \phi_t(h')}\right)_+ d\Phi_{K,t}^j(h') dQ_{K,t}^j(h') 1_{\Omega_t^j \cap \Xi_t^j} \\
&\leq \bP_0 \int_{\Omega_t^j \cap \Xi_t^j}\sup_{h,h' \in K} \psi_t^j(h, h')  d\Phi_{K,t}^j(h') dQ_{K,t}^j(h') \\
&\leq \eta.
\end{align*}

Combination with (<ref>) shows that for all compact $K \subset \R^p$ containing a neighborhood of $0$,

\begin{equation} \label{eqn:tv-l1-convergence}
\bP_0 \| Q_{K,t}^j - \Phi_{K,t}^j \|_{TV} 1_{\Xi_t^j} \to 0.
\end{equation}

Now, let $(K_n)$ be a sequence of balls centered at $0$ with radii $M_n \to \infty$. For each $n \geq 1$, the Equation (<ref>) holds. Hence, the intuition is that if we choose a sequence of balls $(K_t)$ that traverses the sequence $K_n$ slowly enough, the convergence of the expected total variational distance to zero can still be guaranteed. Moreover, the corresponding events $\Xi_t^j = \{Q_t^j(K_t) > 0\}$ satisfy $\bP_0(\Xi_t^j) \to 1$ as a result of assumption (<ref>).

We conclude that there exists a sequence of radii $(M_t)$ such that $M_t \to \infty$ and $\bP_0 \| Q_{K_t,t}^j - \Phi_{K_t,t}^j \|_{TV} \pto 0$ (where the conditional probabilities on the l.h.s. are well-defined on sets of probability growing to one). The total variation distance between a measure and its conditional version on a set $K$ satisfies $\|Q - Q_K\|_{TV} \leq 2Q(K^c)$. Combining this with assumption (<ref>) and the Hellinger integral test, we conclude that

\begin{equation*}
\bP_0 \| Q_t^j - \Phi_t^j \|_{TV} \to 0,
\end{equation*}

which implies the main result.

§.§ Proofs of Results in Section <ref>

For any function $f \geq 0$ and two probability measures $P, Q$, we have

\begin{equation*}
\int f(x) dQ(x) \leq D_{KL}(Q \parallel P) + \log \int \exp(f(x)) dP(x).
\end{equation*}

By the definition of KL divergence,

\begin{align*}
D_{KL}(Q \parallel P) + \log \exp(f(x)) d P(x) &= \int \log\left(\frac{dQ(x) \int \exp(f(y)) dP(y)}{dP(x)} dQ(x) \right) \\
&= \int \log\left(\frac{dQ(x) \int \exp(f(y)) dP(y)}{\exp(f(x)) dP(x)} dQ(x) \right)  + \int f(x) dQ(x) \\
&= D(Q \parallel P') + \int f(x) dQ(x) \geq \int f(x) dQ(x)
\end{align*}

where $P'$ is given by

\begin{equation*}
dP'(x) = \frac{\exp(f(x)) dP(x)}{\int \exp(f(y) dP(y)}.
\end{equation*}

For $P_t^j$ defined in (<ref>) and $P_t$ defined in (<ref>), we have the following upper bound on the expected loss under $P_t^j$ and $P_0$.

\begin{equation*}
\bP_0 P_t^j d(\theta, \theta_0) \leq  \inf_{a > 0} \frac{1}{a} \left[\frac{1}{mt} \left(\bP_0 D_{KL}(P_t^j \parallel P_t) + \log \bP_{\theta_0} P_te^{a mt d(\theta, \theta_0) }\right) + \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j) \right].
\end{equation*}

By Lemma <ref>, we have

\begin{equation*}
amtP_t^jd(\theta, \theta_0)  \leq D_{KL}(P_t^j \parallel P_t) + \log P_t\exp(amt d(\theta, \theta_0)),
\end{equation*}

for all $a > 0$.
Taking expectations on both sides, we have

\begin{equation*}
amt\bP_0 P_t^jd(\theta, \theta_0)  \leq \bP_0 D_{KL}(P_t^j \parallel P_t) + \bP_0 \log P_t\exp(amt d(\theta, \theta_0)),
\end{equation*}

By Lemma <ref> and the tenderization property of KL divergence, we have

\begin{align*}
\bP_0 \log P_t\exp(amtd(\theta, \theta_0)) &\leq mt \left( \frac{1}{m} \sum_{j = 1}^m D_{Kl}(P_0 \parallel P_{\theta_0}^j) \right) + \log \bP_{\theta_0} P_t\exp(amt d(\theta, \theta_0)).
\end{align*}

Putting it all together, we have

\begin{equation*}
\bP_0 P_t^jd(\theta, \theta_0) \leq  \inf_{a > 0} \frac{1}{a mt} \left(\bP_0 D_{KL}(P_t^j \parallel P_t) + \log \bP_{\theta_0} P_t\exp(a mt d(\theta, \theta_0) ) \right) + \frac{1}{am}  \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j).
\end{equation*}

Let the assumptions of Theorem <ref> hold. We have

\begin{equation*}
\bP_{\theta_0} P_t(d(\theta, \theta_0) > C_1 \epsilon^2) \leq  \exp(-C_0 t \epsilon^2) + \exp(-\lambda t \epsilon^2) + 2 \exp(-t \epsilon^2),
\end{equation*}

for all $\epsilon \geq \epsilon_{m,t}^2$, where $\lambda = \rho -1$ in assumption (<ref>).

We define the sets

\begin{equation*}
U = \{\theta: d(\theta, \theta_0)  > C_1 \epsilon^2\}, \qquad K_t = \{\theta: \frac{1}{m}\sum_{j = 1}^m D_{1 + \lambda}(\bP_{\theta_0}^j \parallel \bP_\theta^j) \leq C_3 \epsilon_{m,t}^2\}.
\end{equation*}

Let $\tilde \Pi$ be a probability measure defined as $\tilde \Pi(B) = \frac{\Pi(B \cap K_t)}{\Pi(K_t)}$.  Define the event

\begin{equation*}
A_t = \{X^{(mt)}: \int_\Theta \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \tilde \Pi(\theta) \leq \exp(-(C_3 + 1) t \epsilon^2)\},
\end{equation*}

Let $\Pi_t^j$ and $\phi_t$ be the set and testing function in (<ref>). We bound $\bP_{\theta_0}  P_t (U)$ by

\begin{align*}
\bP_{\theta_0}  P_t (U)& \leq \bP_{\theta_0} \phi_t(X^{(mt)}) + \bP_{\theta_0}(A_t) + \bP_{\theta_0} (1 - \phi_t(X^{(mt)})) P_t(U) I_{(A_t)^c} \\
&= \bP_{\theta_0} \phi_t(X^{(mt)})
+ \bP_{\theta_0}(A_t)
+ \bP_{\theta_0} (1 - \phi_t(X^{(mt)})) I_{(A_t)^c} \frac{\int_{U} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)}{\int_{\Theta} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)}.
\end{align*}

By (<ref>), we bound the first term by

\begin{equation*}
\bP_{\theta_0} \phi_t(X^{(mt)})  \leq \exp(-C_0 t \epsilon^2).
\end{equation*}

By Jensen's inequality, we have

\begin{align*}
D_\rho(\bP_{\theta_0} \parallel \bP_\theta) &= \sum_{i = 1}^m \frac{1}{\rho - 1} \log \int [dP_{\theta_0}^i]^{\frac{\rho}{m}} [dP_{\theta_0}^i]^{\frac{1}{m} - \frac{\rho}{m}} d \mu \\
&\leq \frac{1}{m} \sum_{i = 1}^m \frac{1}{\rho - 1} \log \int [dP_{\theta_0}^i]^\rho [dP_{\theta_0}^i]^{1-\rho} d \mu  = \frac{1}{m} \sum_{i = 1}^m  D_\rho(\bP_{\theta_0}^i \parallel \bP_\theta^i)
\end{align*}

Using the definitions of event $A_t$ and Markov inequality, we have

\begin{align*}
\bP_{\theta_0} (A_t) &=  \bP_{\theta_0}\left(\left(\int_\Theta \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \tilde \Pi(\theta) \right)^{-\lambda } < \exp(\lambda(C_3 + 1) t \epsilon^2 \right)  \\
&\leq \exp(-\lambda (C_3 + 1) t \epsilon^2)  \bP_{\theta_0}\left(\int_\Theta \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \tilde \Pi(\theta) \right)^{-\lambda } \\
&\leq \exp(-\lambda(C_3 + 1) t \epsilon^2) \int_\Theta \left( \int_{\X^{mt}}\frac{(\bp_{\theta_0}(x^{(mt)}))^{\lambda }}{(\bp_\theta(x^{(mt)}))^{\lambda } }  d \bP_{\theta_0}(x^{(mt)})\right) d \tilde \Pi(\theta) \\
&= \exp(-\lambda(C_3 + 1) t \epsilon^2) \int_\Theta \exp(\lambda  D_{1 + \lambda}(\prod_{k = 1}^t  \bP_{\theta_0} \parallel \prod_{k = 1}^t\bP_\theta)) d \tilde \Pi(\theta) \\
&= \exp(-\lambda(C_3 + 1) t \epsilon^2) \int_{K_t} \exp(\lambda  t  D_{1 + \lambda}(  \bP_{\theta_0} \parallel \bP_\theta)) d \tilde \Pi(\theta) \\
&\leq \exp(-\lambda(C_3 + 1) t \epsilon^2) \int_{K_t} \exp(\lambda t \frac{1}{m} \sum_{j = 1}^m D_{1 + \lambda}(\bP_{\theta_0}^j \parallel \bP_\theta^j)) d \tilde \Pi(\theta) \\
&\leq \exp(- \lambda(C_3 + 1) t \epsilon^2 + \lambda C_3 t \epsilon_{m,t}^2) \\
&\leq \exp(-\lambda t \epsilon^2).
\end{align*}

Let's analyze the third term. Conditioned on $(A_t)^c$, we have

\begin{equation*}
\int_{\Theta} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta) \geq \Pi(K_t)  \int_{\Theta} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \tilde \Pi(\theta)  \geq \exp(-(C_2 + C_3 + 1) t \epsilon^2),
\end{equation*}

where the last inequality follows from <ref>. It follows that

\begin{align*}
&\bP_{\theta_0} (1 - \phi_t(X^{(mt)})) I_{(A_t)^c} \frac{\int_{U} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)}{\int_{\Theta} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)} \\
&\leq \exp((C_2 + C_3 + 1) t \epsilon^2)\bP_{\theta_0} \int_{U} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta) \\
&\leq \exp((C_2 + C_3 + 1) t \epsilon^2)\left[\int_{U \cap \Theta_t(\epsilon)} \bP_\theta (1- \phi_t (X^{(mt)})) d \Pi(\theta) + \Pi\left(\Theta_t(\epsilon)^c\right) \right] \\
& \leq \exp((C_2 + C_3 + 1) t \epsilon^2)(\exp(-C t \epsilon^2) + \exp(-C t \epsilon^2)),
\end{align*}

where the last inequality follows from <ref> and <ref>. Since $C_0 > C_3 + C_2 + 2$, we obtain the upper bound for the third term.

\begin{equation*}
\bP_{\theta_0} (1 - \phi_t(X^{(mt)})) I_{(A_t)^c} \frac{\int_{U} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)}{\int_{\Theta} \frac{\bp_\theta}{\bp_{\theta_0}}(X^{(mt)}) d \Pi(\theta)}  \leq 2 \exp(- t \epsilon^2).
\end{equation*}

Putting the bounds all together, we have

\begin{equation*}
\bP_{\theta_0}  P_t (U) \leq \exp(-C_0 t \epsilon^2) + \exp(-\lambda t \epsilon^2) + 2 \exp(-t \epsilon^2)
\end{equation*}

Let the random variable $X$ satisfies

\begin{equation*}
\bP(X \geq \delta) \leq c_1 \exp(-c_2 \delta),
\end{equation*}

for all $\delta \geq \delta_0 > 0$. Then, for any $0 < a \leq \frac{1}{2} c_2$, we have

\begin{equation*}
\E \exp(a X) \leq \exp(a \delta_0) + c_1,
\end{equation*}

Define $Y = \exp(aX)$ for some $0 < a \leq \frac{1}{2}c_2$. For any $C > 0$,

\begin{align*}
\E Y \leq C + \int_C^\infty \bP(Y \geq y) dy = C + \int_C^\infty \bP(X \geq \frac{1}{a} \log y) dy \leq C + c_1 \int_C^\infty y^{-c_2/a} dy.
\end{align*}

Take $C = \exp(a \delta_0)$. SInce $a \leq \frac{1}{2} c_2$, we have

\begin{equation*}
\E Y \leq \exp(a_t) + c_1 \exp(-a\delta_0) \leq \exp(a\delta_0) + c_1.
\end{equation*}

By Lemma <ref>, for all $\delta \geq \delta_0$, we have

\begin{equation*}
\bP_{\theta_0} P_t(d(\theta, \theta_0) > \delta) \leq c_1 \exp(-c_2 \delta).
\end{equation*}

Take $c_1 = 4, c_2 = mt\min(\lambda, 1)/C_1$ as $C_0 > C_1 + C_2 + 2 > 1$ and $\delta_0 = C_1 m \epsilon_{m,t}^2$. By Lemma <ref>, we have

\begin{equation*}
\bP_{\theta_0} P_t\exp(a mtd(\theta, \theta_0))  \leq \exp(a mt C_1\epsilon_{m,t}^2) + 4,
\end{equation*}

for all $a \leq \min(\lambda, 1)/(2 C_1)$.

Take $a = \min(\lambda, 1)/(2 C_1)$. By Lemma <ref>, we have

\begin{align*}
\bP_0  P_t^j d(\theta, \theta_0)  &\leq \frac{mt \gamma_{j,m,t}^2 + \log(4 + \exp(a C_1 mt \epsilon_{m,t}^2)) + t \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j)}{a mt}  \\
&\leq \frac{\gamma_{j,m,t}^2}{a} + C_1 \epsilon_{m,t}^2 + \frac{4 \exp(-a C_1 mt \epsilon_{m,t}^2)}{amt} + \frac{ \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j)}{a m}\\
&\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{align*}

for some $C$ that depends only on $C_0, C_1, \lambda$.

The first result is a consequence of Markov's inequality.

\begin{align*}
&\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) \\
&\leq \frac{\bP_0 P_t^j  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))} \leq \frac{C}{M_t}  \to 0
\end{align*}

The second result follows from Jensen's inequality

\begin{equation*}
\bP_0 d(P_t^j \theta, \theta_0) \leq \bP_0 P_t 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*}

We have

\begin{align*}
\bP_0 D_{KL}(P_t^j \parallel P_t) &=  \bP_0 P_t^j \left( \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}] \log \bP_\theta^i (X_k^i) -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \log \bP_\theta^i (X_k^i) \right) \\
&=P_t^j \bP_0 \left( \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}] \log \bP_\theta^i (X_k^i) -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \log \bP_\theta^i (X_k^i) \right) \\
&= P_t^j  \left( \sum_{k = 1}^t \sum_{i = 1}^m [A^{t - k}_{ij}] \bP_0 \log \bP_\theta^i -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \bP_0 \log \bP_\theta^i \right) \\
&= \sum_{k = 1}^t \sum_{i = 1}^m
\left\{  [A^{t - k}_{ij}] -  \frac{1}{m} \right\} P_t^j  \bP_0 \log \bP_\theta^i \\
&= \sum_{k = 1}^t \sum_{i = 1}^m
\left\{  [A^{t - k}_{ij}] -  \frac{1}{m} \right\} \left\{\bP_0\log \bp_0  - P_t^j D_{KL}(\bP_0 \parallel \bP_\theta^i)\right\}.
\end{align*}

For each $j$, we have

\begin{align*}
\bP_0 D_{KL}(P_t^j \parallel P_t) &\leq  \sum_{k = 1}^t \sum_{i = 1}^m | [A^{t - k}_{ij}] -  \frac{1}{m}| \left|\bP_0\log \bp_0 - \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right| \\
& \leq  \sum_{k = 1}^t \sum_{i = 1}^m | [A^{t - k}_{ij}] -  \frac{1}{m}| \left\{|\bP_0\log \bp_0| + \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right\}.
\end{align*}

By Lemma <ref>, we get

\begin{equation*}
\bP_0 D_{KL}(P_t^j \parallel P_t)  \leq \frac{16 m^2 \log m}{\nu} \left(|\bP_0\log \bp_0| + \frac{1}{m}\sum_{i = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^i)\right).
\end{equation*}

By Assumption <ref>(a), 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*}

§.§ Proofs of Results in Section <ref>

Define $\mathcal T \subseteq [t]$ as the set of $\tau$ where $G_\tau = G$.  This allows us to break the left-hand side of the inequality into two parts:

\begin{align*}
\sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}- \frac{\lambda}{m} \right|  &=
\sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau \in [k, t-1]\cap \mathcal T} A \right]_{ij}- \frac{\lambda}{m} \right| =  \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[ A^{|[k, t-1]\cap \mathcal T|} \right]_{ij}- \frac{\lambda}{m} \right|.
\end{align*}

As $t \to \infty$, $|\mathcal T|/t \asto \lambda$ and $|[k, t-1]\cap \mathcal T| \asto \lambda (t - k)$.  Then with probability $1$, we have

\begin{align*}
\limsup_{t \to \infty} \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij}- \frac{\lambda}{m} \right|
&=  \limsup_{t \to \infty}  \sum_{k = 1}^t \sum_{j = 1}^m \left|\left[ A^{\lambda(t - k)} \right]_{ij}- \frac{\lambda}{m} \right|.
\end{align*}

Define $\delta = \max(|\lambda_2(A)|, |\lambda_m(A)|)$. The spectral radius of $A$ is $1$, thus $\delta <1$ by Perron-Frobenius theorem.  Under the assumptions <ref>, by the standard property of stochastic matrices (see e.g. [Rosenthal, 1995]), the diagonalizable matrix $A$ satisfies

\begin{equation} \label{eqn:graph_conv-1}
\|e_i^T A^{\lambda t} - \frac{1}{m} \ones \|_1 \leq m \delta^{\lambda t}.
\end{equation}

for any $i \in [m]$, using the fact that $\frac{1}{m} \ones$ is the stationary distribution of the Markov chain with transition matrix $A$.

Assume that $\lambda \geq \frac{c}{m}$.  For any $t - k \geq \tilde t = \frac{\log \frac{m}{c} + \log \lambda}{- \lambda \log \delta}$,

\begin{equation*}
m \delta^{\lambda(t-k)} \leq  m \delta^{\lambda \tilde t} \leq \frac{m}{\frac{m}{c} \lambda} \leq \frac{c}{\lambda}.
\end{equation*}

Since $ \|e_i^T A^{\lambda t} - \frac{1}{m} \ones \|_1 \leq 2$ by the double stochasticity of $A$, we use (<ref>) to break the quantity of interest $\sum_{k = 1}^t \sum_{j = 1}^m |[A^{\lambda(t-k)}_{ij}] - \frac{1}{m}|$ into two parts. For any $i \in [m]$,

\begin{align*}
\sum_{k = 1}^t \sum_{j = 1}^m |[A^{\lambda(t-k)}_{ij}] - \frac{1}{m}|  &=  \sum_{k = 1}^t \| e_i^T A^{\lambda(t-k)}- \frac{1}{m} \ones \|_1  \\
&=  \sum_{k = 1}^{t - \tilde t} \| e_i^T A^{\lambda(t-k)} - \frac{1}{m} \ones \|_1 + \sum_{k > t - \tilde t}^{t} \| e_i^T A^{\lambda(t-k)} - \frac{1}{m} \ones \|_1 \\
&\leq \sum_{k = 1}^{t - \tilde t}  m \delta^{\lambda(t-k)} + 2 \tilde t \leq \frac{m \delta^{\lambda \tilde t}}{1 - \delta} + 2 \tilde t \leq \frac{c}{\lambda(1 - \delta)} + \frac{2 \log \frac{m}{c}+2 \log \lambda}{-\lambda \log \delta}.
\end{align*}

Noting that $1 - \delta \leq - \log \delta$ and $\delta = 1 - \frac{\nu}{4 m^2}$ for any doubly stochastic matrix $A$,  we get

\begin{align*}
\sum_{k = 1}^t \sum_{j = 1}^m |\left[ A^{\lambda(t - k)} \right]_{ij} - \frac{1}{m}| & \leq  \frac{c + 2 \log \frac{m}{c} + \log \lambda}{\lambda(1 - \delta)}  \leq  \frac{4m^2 c + 8 m^2 \log \frac{m}{c} + 8 m^2\log \lambda}{\lambda \nu}.
\end{align*}

When $\lambda \geq \frac{2}{m}$, the optimal bound is achieved when $c = \frac{1}{2}$,

\begin{align*}
\limsup_{t \to \infty}  \sum_{k = 1}^t \sum_{j = 1}^m |\left[ A^{\lambda(t - k)} \right]_{ij} - \frac{1}{m}| &\leq   \frac{8m^2 (1 + \log \frac{m}{2}) + 8 m^2\log \lambda}{\lambda \nu} \leq \frac{16m^2 \log m + 8 m^2\log \lambda}{\lambda \nu}.
\end{align*}

When $\lambda < \frac{2}{m}$,  the optimal bound is achieved when $c = \lambda $,

\begin{align*}
\limsup_{t \to \infty}  \sum_{k = 1}^t \sum_{j = 1}^m |\left[ A^{\lambda(t - k)} \right]_{ij} - \frac{1}{m}| &\leq  \frac{4m^3}{\nu}.
\end{align*}

Finally, when $\lambda = 0$,

\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{\lambda}{m} \right|
=  \limsup_{t \to \infty} \sum_{k = 1}^t  \frac{m-1}{m} = \infty.
\end{equation*}

We have

\begin{align*}
\bP_0 D_{KL}(P_t^j \parallel P_t) &=  \bP_0 P_t^j \left( \sum_{k = 1}^t \sum_{i = 1}^m \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} \log \bP_\theta^i (X_k^i) -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \log \bP_\theta^i (X_k^i) \right) \\
&=P_t^j \bP_0 \left( \sum_{k = 1}^t \sum_{i = 1}^m \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} \log \bP_\theta^i (X_k^i) -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \log \bP_\theta^i (X_k^i) \right) \\
&= P_t^j  \left( \sum_{k = 1}^t \sum_{i = 1}^m \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} \bP_0 \log \bP_\theta^i -  \frac{1}{m}\sum_{k = 1}^t \sum_{i = 1}^m  \bP_0 \log \bP_\theta^i \right) \\
&= \sum_{k = 1}^t \sum_{i = 1}^m
\left\{  \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} -  \frac{1}{m} \right\} P_t^j  \bP_0 \log \bP_\theta^i \\
&= \sum_{k = 1}^t \sum_{i = 1}^m
\left\{  \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} -  \frac{1}{m} \right\} \left\{\bP_0\log \bp_0  - P_t^j D_{KL}(\bP_0 \parallel \bP_\theta^i)\right\}.
\end{align*}

For each $j$, we have

\begin{align*}
\bP_0 D_{KL}(P_t^j \parallel P_t) &\leq  \sum_{k = 1}^t \sum_{i = 1}^m | \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} -  \frac{1}{m}| \left|\bP_0\log \bp_0 - \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right| \\
& \leq  \sum_{k = 1}^t \sum_{i = 1}^m | \left[\prod_{\tau = k}^{t-1}A_\tau \right]_{ij} -  \frac{1}{m}| \left\{|\bP_0\log \bp_0| + \max_{i \in [m]}\inf_{\theta \in \Theta} D_{KL}(\bP_0 \parallel \bP_\theta^i) \right\}.
\end{align*}

This implies that

\begin{equation*}
\bP_0 D_{KL}(P_t^j \parallel P_t)  \leq    \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| \left(|\bP_0\log \bp_0| + \frac{1}{m}\sum_{i = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^i)\right).
\end{equation*}

Applying Proposition <ref> gives us the desired results.

Since $P_t$ does not depend on the communication graph, Theorem <ref> still holds.  This implies the following upper bound on the contraction rate:

\begin{equation} \label{eqn-proofcor-time-varying2:1}
\bP_0 P_t^j d(\theta, \theta_0) \lesssim \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}

The only term depending on the communication graph is $\gamma_{j,m,t}^2$. By Corollary <ref>, we have the following upper bounds with probability $1$ (with respect to the probability measure of $A_1, A_2, \cdots$).

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

\begin{equation*}
\gamma_{j,m,t}^2 \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*}
\gamma_{j,m,t}^2 \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*}

Combining the upper bounds with Equation (<ref>), we have:

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

\begin{equation*}
\bP_0 P_t^j d(\theta, \theta_0) \lesssim  \frac{1}{t} + \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) + \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{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)+ \frac{1}{m} \sum_{j = 1}^m D_{KL}(\bP_0 \parallel \bP_{\theta_0}^j).
\end{equation*}

This is the desired result after we remove the constants.

§.§ Proofs of Results in Section <ref>

By Definition (<ref>),

\begin{equation*}
\log p_t^j(\theta) = t f_t^j(\theta)+ \pi(\theta) + const.
\end{equation*}

for $f_t^j$ defined in Equation (<ref>).

We can compute $f_t^j(\theta)$ explicitly ,

\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 \sum_{i = 1}^m [A^{t-k}_{ji}] \left[ \langle \theta, T^i(X_k^i) \rangle - \psi^i(\theta)  \right] + const  \\
&= -\frac{1}{t}\left[ \langle \theta,  \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  T^i(X_k^i) \rangle -  \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  \psi^i(\theta)  \right] + const.
\end{align*}

This implies that

\begin{align*}
\log p_t^j(\theta) &= -t f_t^j(\theta)+ \pi(\theta)  + const \\
&=  \langle \theta,   \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  T^i(X_k^i) \rangle - \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  \psi^i(\theta)   + \langle \theta, u\rangle - \psi^0(\theta) + const\\
&=  \langle \theta,   \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  T^i(X_k^i) + u\rangle -  \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}]  \psi^i(\theta) - \psi^0(\theta) + const.
\end{align*}

This shows that $p_t^j$ is a member of the exponential family of the form (<ref>) with sufficient statistic

\begin{equation*}
\chi_t^j + u,
\end{equation*}

and log-partition function

\begin{equation*}
B_t^j(\theta)  + \psi^0(\theta),
\end{equation*}

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

By Lemma <ref>, $P_t^j$ is a member of the exponential family given by Equation (<ref>).
The full rankness of $P_t^j$ follows directly from the definition of full-rank exponential family and the compact representation (<ref>). By Lemma 4.5 of [Van der Vaart, 2000], the log-partition function $B_t^j: \Theta \to R$ is infinitely times differentiable with the gradient being the mean parameter of $\chi_t^j$ and the Hessian is the covariance matrix of $\chi_t^j$. Since $P_t^j$ is full-rank, it is equivalent to the nonsingularity of the covariance matrix for $\chi_t^j$ [Van der Vaart, 2000]. This implies that the gradient $\nabla_\theta B_t^j$ is one-to-one in the interior of $\Theta$, thus, $\nabla_\theta B_t^j$ is invertible.

By Lemma (<ref>), the sequence of estimators $\hat \theta_t^j$ is well-defined with $[\bP_0]-$probability tending to $1$. By Theorem 4.6 of [Van der Vaart, 2000], the sequence of estimators $\hat \theta_t^j$ has limiting distribution

\begin{equation*}
\sqrt{t} (\hat \theta_t^j - \theta_0) \dto N(0, I_{\theta_0}^{j^{-1}}),
\end{equation*}

where $I_{\theta_0}^j$ is the Fisher information for $\bP_\theta^j$ at $\theta_0$.  This corresponds to the Sandwich covariance matrix.

\begin{equation*}
I_{\theta_0}^j = \left[ \phi^{'^{-1}}_{\theta_0} \text{Cov}_{\theta_0}(\chi_t^j) \phi^{-1}_{\theta_0} \right]^{-1}  = \text{Cov}_{\theta_0}(\chi_t^j).
\end{equation*}

Since $B_t^j$ is infinitely times differentiable in a neighborhood of $\theta_0$, the Fisher information $I_{\theta}^j$ must be uniformly bounded. Thus, $\nabla \log \bp_\theta^j$ is Lipchistz around $\theta_0$ and Assumption <ref>(d) is satisfied.

When $\text{int}(\Theta) \neq \emptyset$, the exponential family $\bP_\Theta^i$ satisfies the differential in the quadratic mean Assumption <ref>(c) by Example 7.7 of [Van der Vaart, 2000]. Given that $\theta_0 \in \text{int}(\Theta)$ and $\pi$ have support on $\Theta$, the prior $\pi$ puts positive mass on every neighborhood of $\theta_0$, thus Assumption <ref>(b) is satisfied.

Under the full-rank assumptions on $\bP_\theta^j$,  the Hessian $\nabla^2_\theta \log \bp_\theta^j$ is negative definite for all $\theta$. Thus, the Hessian of $f_t^j$ is negative definite as a positive linear combination of $\nabla^2_\theta \log \bp_\theta^j$. We have that $f_t^j$ is strictly concave on $\Theta$ with a unique maximizer at $\tilde \theta_t^j$. The strict concavity implies 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*}

By Lemma <ref>, $\tilde \theta_t^j$ is consistent at $\theta_0$. Then Assumption <ref>(a) is satisfied as a consequence of  Lemma <ref>.

The assumptions of Theorem <ref> are all satisfied. Considering the sequence of random variables $\theta \sim P_t^j$ centered at the moment estimator $\hat \theta_t^j$, we have

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

where $q_t^j$ is the density of $\sqrt{t}(\theta -\hat \theta_t^j)$, and $V_{\theta_0} = \frac{1}{m} \sum_{i = 1}^m \text{Cov}(T^i)$.

The prior $\pi$ is assumed to have full support over $\R^p$, which trivially satisfies Assumption <ref>(a).

Define functions $f_t^j$, $f_t$, and $f$ on $\Theta$ as per Equations (<ref>)–(<ref>). By Assumption (<ref>)(a), the function $f$ exists and we note that $f_t^j(\theta) \pto f(\theta)$.

The estimators $\hat \theta_t^j$ and the true parameter $\theta_0$ are respectively given by:

\begin{equation*}
\hat \theta_t^j = \argmin_{\theta \in \Theta} f_t^j(\theta), \quad \theta_0 = \argmin_{\theta \in \Theta} f(\theta).
\end{equation*}

The gradient of $f_t^j$ is given by

\begin{equation*}
\nabla  f_t^j(\theta) = -\frac{1}{t} T_t^j + \frac{1}{t} \sum_{i = 1}^m \sum_{k = 1}^t \frac{ [A^{t-k}_{ji}] X_k^j e^{\langle \theta,  X_k^j \rangle}}{1 + e^{\langle \theta,  X_k^j \rangle}},
\end{equation*}

and the Hessian is

\begin{equation*}
%  \label{eqn:hessian-ftj}
\nabla^2 f_t^j(\theta) = \frac{1}{t} \sum_{i = 1}^m \sum_{k = 1}^t \frac{ [A^{t-k}_{ji}] X_k^j X_k^{j^T} e^{\langle \theta,  X_k^j \rangle}}{\left(1 + e^{\langle \theta,  X_k^j \rangle} \right)^2}.
\end{equation*}

Assumption i) implies that
Given that $\nabla^2 f_t^j(\theta) > 0$, the function $f_t^j$ is strictly convex in $\theta$. This implies that the minimizer $\hat \theta_t^j$ is unique, and we have $\hat \theta_t^j \in \text{int}(\Theta)$.

The gradient of $f$ is given by

\begin{equation*}
\nabla  f(\theta) = \bP_{\theta_0}\left[\frac{X_1^1 e^{\langle \theta,  X_1^1 \rangle}}{1 + e^{\langle \theta,  X_1^1 \rangle}} - X_1^1 Y_1^1 \right] =  \bP_{\theta_0}\left[X_1^1 \left( \frac{e^{\langle \theta,  X_1^1 \rangle}}{1 + e^{\langle \theta,  X_1^1 \rangle}} - \frac{e^{\langle \theta_0,  X_1^1 \rangle}}{1 + e^{\langle \theta_0,  X_1^1 \rangle}}\right) \right],
\end{equation*}

The differentiation and expectation are interchanged by the Lebesgue-dominated convergence theorem since $ |\nabla f(\theta)| \leq |X_1^1|$ and $\bP_0 |X_1^1| < \infty$ by Assumption ii).

Since $\nabla f(\theta)$ is strictly increasing in $\theta$, we have $\nabla f(\theta - \epsilon) < 0 < \nabla f(\theta + \epsilon)$ in each coordinate. Hence, Assumption <ref>(b) is satisfied. By Lemma <ref>, we have $\hat \theta_t^j \pto \theta_0 $.

Let $\mathcal{E} = \{\eta \in \mathbb{R} : |\sigma(\eta)| < \infty \}$. The set $\mathcal{E}$ is open and nonempty. Additionally, $\eta$ is identifiable, as $\sigma(\eta)$ is a one-to-one function. Trivially, the set $\Theta$ is open and convex, and $\theta^T x \in \mathcal{E}$ for all $\theta \in \Theta$ and $x \in \X$.

For all $\eta \in \mathcal{E}$, we have $0 < \sigma(\eta) < 1$ and

\begin{equation*}
|\sigma'''(\eta)| = |\frac{\sigma(1 - \sigma(\eta))(1 - 2\sigma(\eta))^2 - 2\sigma^2(1-\sigma(\eta))^2}{(1 - \sigma(\eta))^2}| \leq 3.
\end{equation*}

After algebraic manipulations, we have

\begin{equation*}
\nabla_\theta^3 \log p_\theta^j(. \mid x_k^j)_{a,b, c} = \sigma^{'''}(\langle \theta, x_k^j \rangle) x_{k,a}^j x_{k,b}^j x_{k,c}^j \leq 3 x_{k,a}^j x_{k,b}^j x_{k,c}^j,
\end{equation*}

for all $\theta \in \Theta, x_k^j \in \X$ and $a, b,c \in [p]$.  Thus, $\nabla_\theta^3 \log p_\theta^j$ is uniformly bounded by an integrable function.

The Fisher information for agent $j$ is given by

\begin{equation*}
V_{\theta_0}^j = - \bP_0 \nabla_\theta^2 \log p_\theta^j(. \mid X_k^j) = \bP_0  \left[\frac{X_k^{j^T} e^{\langle \theta,  X_k^j \rangle} X_k^j}{\left(1 + e^{\langle \theta,  X_k^j \rangle} \right)^2} \right].
\end{equation*}

By Assumption i), $V_{\theta_0}^j$ exists and is nonsingular. Hence, Assumption <ref>(f) is satisfied.

Assumption <ref>(c), <ref>(a) are satisfied with the same argument as Proposition <ref>.  By Corollary <ref> and Slusky's theorem, we have

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

where $\hat V_{\theta_0}$ is given by as the finite - sample version of the $\frac{1}{m} \sum_{j = 1}^m V_{\theta_0}^j$. Specifically,

\begin{equation*}
\hat V_\theta = \frac{1}{m} \sum_{i = 1}^m {X^i}^T  \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) X^i.
\end{equation*}

By Lemma <ref>, $f_t^j(\theta) \pto f(\theta)$ for each $\theta \in \Theta$.  For all $\theta \in [0,1]^2$, $f_t^j$ is uniformly bounded.

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

The gradient $\nabla f_t^j(\theta)$ is uniformly bounded over $[0,1]^2$. To see this, consider

\begin{align*}
\nabla f_t^j(\theta) = A + B,
\end{align*}

where $A, B$ are given by

\begin{align*}
A &=  \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \left[\frac{\theta - Z^j}{\sigma^{j^2}} - \frac{|X_t^j - Z^j| (\theta - Z^j)}{|\theta - Z^j|\sigma^{j^2}}\right], \\
B &=  \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \frac{(\theta - Z^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]}{\sigma^j |\theta - Z^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]}.
\end{align*}

Applying the following upper bounds on $A$ and $B$,

\begin{align*}
A &\leq \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \frac{\sqrt{2}}{\sigma^{j^2}}, \\
B &\leq \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] \frac{\phi\left(\frac{|Z^j| + \frac{1}{2}}{\sigma^j}\right) - \phi\left(\frac{-\sqrt{2}}{\sigma^j}\right)}{\sigma^j\left[\Phi\left(\frac{|Z^j| + \frac{1}{2} - \sqrt{2}}{\sigma^j}\right) - \Phi\left(\frac{-\sqrt{2}}{\sigma^j}\right)\right]},
\end{align*}

we obtain

\begin{equation*}
\|\nabla f_t^j(\theta)\| = \|A + B\| \leq \left[ \frac{\sqrt{2}}{ \sigma^{j^2}} + \frac{\phi(\frac{|Z^j|+ \frac{1}{2}}{\sigma^j}) - \phi(\frac{- \sqrt{2}}{\sigma^j}) }{\sigma^j\left[\Phi(\frac{|Z^j|+ \frac{1}{2} - \sqrt{2}}{\sigma^j}) - \Phi(\frac{- \sqrt{2}}{\sigma^j}) \right]} \right]\left(1 + \frac{16 m^2 \log m}{\nu} \right).
\end{equation*}

Then $f_t^j$ is uniformly equicontinuous in $t$, as $f_t^j$ is Lipschitz.  By Lemma <ref>, we have $\|f_t^j - f\|_\infty \pto 0$.

Note that $\hat{\theta}_t^j = \argmin_{\theta \in [0,1]^2} f_t^j(\theta)$ and $\theta_0 = \argmin_{\theta \in [0,1]^2} f(\theta)$. By the Argmax theorem (Theorem 3.2.2, [Vaart and Wellner, 2023]), we have $\hat{\theta}_t^j \pto \theta_0$.

Since $\theta_0 \in \Theta$, for small enough $\epsilon$, we have

\begin{align*}
\bP_0(|\hat \theta_t^j - \theta_0| > \epsilon) &= \bP_0(|\hat \theta_t^j  - \theta_0| > \epsilon, \hat \theta_t^j  \notin \Theta) + \bP_0(|\hat \theta_t^j - \theta_0| > \epsilon, \hat \theta_t^j \in \Theta)  \\
&= \bP_0(\hat \theta_t^j \notin \Theta) + \bP_0(|\hat \theta_t^j - \theta_0| > \epsilon, \hat \theta_t^j \in \Theta) \to 0,
\end{align*}

thus $\bP_0(\hat \theta_t^j \notin \Theta) \to 0$.

The prior $\pi$ is assumed to have full support over $\Theta$, which trivially satisfies Assumption <ref>(a).

From the proof of Lemma <ref>, we proved that $\nabla f_t^j$ is uniformly bounded for all $t$, thus $f_t^j$ is uniformly equicontinuous. By definition, there exists $\delta > 0$ such that $\|\theta - \hat \theta_t^j\| > \delta$ implies that $|f_t^j(\theta) - f_t^j(\hat \theta_t^j)| > \epsilon$. This verifies Assumption <ref>(a).

The gradient of $\log p_\theta^j$ is given by

\begin{align*}
\nabla \log p_\theta^j(X_t^j) = \frac{\theta - Z^j}{|\theta - Z^j| \sigma^{j^2}} \left[ |\theta - Z^j| - |X_t^j - Z^j| - \frac{\phi\left(\frac{|Z^j| + \frac{1}{2} - |\theta - Z^j|}{\sigma^j}\right) - \phi\left(\frac{- |\theta - Z^j|}{\sigma^j}\right)}{\Phi\left(\frac{|Z^j| + \frac{1}{2} - |\theta - Z^j|}{\sigma^j}\right) - \Phi\left(\frac{- |\theta - Z^j|}{\sigma^j}\right)} \right].
\end{align*}

This allows us to derive the Fisher information of $p_\theta^j$.

\begin{align*}
V_\theta^j &=  \bP_0[\nabla_\theta \log p_\theta^j \nabla_\theta \log p_\theta^{j^T}] \\
&=  \frac{(\theta - Z^j) (\theta - Z^j)^T}{\sigma^{j^4} |\theta - Z^j|^2} \left[ |\theta - Z^j| - |\theta_0 - Z^j| - \frac{\phi(\frac{|Z^j|+ \frac{1}{2} - |\theta - Z^j|}{\sigma^j}) - \phi(\frac{- |\theta - Z^j|}{\sigma^j})}{ \Phi(\frac{|Z^j|+ \frac{1}{2} - |\theta - Z^j|}{\sigma^j}) - \Phi(\frac{- |\theta - Z^j|}{\sigma^j}) } \right]^2.
\end{align*}

At $\theta_0$, the Hessian simplifies to $V_{\theta_0}^j = \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$ which is nonsingular. Thus, Assumption <ref>(e) is satisfied.

Reusing the argument in the proof of Lemma <ref>, we have

\begin{align*}
\sup_{\theta \in [0,1]^2} \|\nabla \log p_{\theta}^j\| \leq \frac{\sqrt{2}}{\sigma^{j^2}} + \frac{\phi(\frac{|Z^j|+ \frac{1}{2}}{\sigma^j}) - \phi(\frac{- \sqrt{2}}{\sigma^j}) }{\sigma^j\left[\Phi(\frac{|Z^j|+ \frac{1}{2} - \sqrt{2}}{\sigma^j}) - \Phi(\frac{- \sqrt{2}}{\sigma^j}) \right]}.
\end{align*}

This implies that $\nabla \log p_\theta^j$ is a continuously differentiable function over a compact domain $[0,1]^2$. Thus, it is Lipschitz continuous. Assumption <ref>(d) is satisfied.

The Hessian for $f_t^j$ is simply the average of agent Fisher information:

\begin{equation*}
\nabla^2 f_t^j(\theta) = \frac{1}{t} \sum_{k = 1}^t \sum_{i = 1}^m [A^{t-k}_{ji}] V_\theta^i.
\end{equation*}

Recall that in the proof of
Since $f_t^j$ is uniformly equicontinuous,

We define the event $E_t$ as

\begin{equation*}
E = \{X^{(mt)}:\hat \theta_t^j \in \Theta\}.
\end{equation*}

Conditioning on $E_t$, by Corollary <ref> and Slusky's theorem, we have

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

for  $q_t^j$ defined as the density of $\sqrt{t}(\theta - \hat{\theta}_t^j)$ and $V_{\theta_0}$ defined in Equation (<ref>).

By Lemma <ref>, $\lim_{t \to \infty}\bP_0(E_t) = 1$. This completes our proof.

§ APPENDIX C. SUPPORTING RESULTS

Define $\delta = \max(|\lambda_2(A)|, |\lambda_m(A)|)$. The spectral radius of $A$ is $1$, thus $\delta <1$ by Perron-Frobenius theorem.  Under the assumptions <ref>, by the standard property of stochastic matrices (see e.g. [Rosenthal, 1995]), the diagonalizable matrix $A$ satisfies

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

for any $i \in [m]$, using the fact that $\frac{1}{m} \ones$ is the stationary distribution of the Markov chain with transition matrix $A$.

For any $t - k \geq \tilde t = \frac{\log \frac{m}{2}}{- \log \delta}$,

\begin{equation*}
m \delta^{t-k} \leq  m \delta^{\tilde t} \leq \frac{m}{\frac{m}{2}} \leq 2
\end{equation*}

Since $ \|e_i^T A^t - \frac{1}{m} \ones \|_1 \leq 2$ by the double stochasticity of $A$, we use (<ref>) to break the quantity of interest $\sum_{k = 1}^t \sum_{j = 1}^m |[A^{t - k}_{ij}] - \frac{1}{m}|$ into two parts, that is, for any $i \in [m]$,

\begin{align*}
\sum_{k = 1}^t \sum_{j = 1}^m |[A^{t - k}_{ij}] - \frac{1}{m}|  &=  \sum_{k = 1}^t \| e_i^T A^{t - k} - \frac{1}{m} \ones \|_1  \\
&=  \sum_{k = 1}^{t - \tilde t} \| e_i^T A^{t - k} - \frac{1}{m} \ones \|_1 + \sum_{k > t - \tilde t}^{t} \| e_i^T A^{t - k} - \frac{1}{m} \ones \|_1 \\
&\leq \sum_{k = 1}^{t - \tilde t}  m \delta^{t-k} + 2 \tilde t \leq \frac{m \delta^{\tilde t}}{1 - \delta} + 2 \tilde t \leq \frac{2}{1 - \delta} + \frac{2 \log \frac{m}{2}}{- \log \delta}
\end{align*}

Noting that $1 - \delta \leq - \log \delta$ and $\delta = 1 - \frac{\nu}{4 m^2}$ for any doubly stochastic matrix $A$,  we get

\begin{equation*}
\sum_{k = 1}^t \sum_{j = 1}^m |[A^{t - k}_{ij}] - \frac{1}{m}|  \leq  \frac{2 + 2 \log \frac{m}{2}}{1 - \delta} \leq  \frac{4 \log m}{1 - \delta} \leq  \frac{16 m^2 \log m}{\nu}
\end{equation*}

Assume that $S_t^i, t \geq 1$ are i.i.d. random variables with $\E[|S_t^i|]$ exists and are finite, for all $ i \in [m]$. Under Assumption <ref>, for any $j \in [m]$, the random variables

\begin{equation*}
Z_t^j = \frac{1}{t} \sum_{i = 1}^m\sum_{k = 1}^t [A_{ij}^{t-k}] S_t^i
\end{equation*}

converge in probability to $\frac{1}{m} \sum_{i = 1}^m \E[S_1^i]$ as $t \to \infty$.

Let $Z_t= [Z_t^1, \cdots, Z_t^m]^T$ and $S_t= [S_t^1, \cdots, S_t^m]^T$.
Since $\E[|S_t^i|] < \infty$ exists and is finite, we have

\begin{equation*}
Z_t = \frac{1}{t} \sum_{k = 1}^t A^{t-k} S_k.
\end{equation*}

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\|Z_t  - \sum_{k = 1}^t \frac{1}{m} \ones \ones^T S_k \right\| = \left\|\frac{1}{t} \sum_{k = 1}^t A^{t-k} S_k - \frac{1}{m} \ones \ones^T S_k \right\| \\
&\leq \frac{1}{t} \sum_{k = 1}^t \|(A^{t-k} -\frac{1}{m} \ones \ones^T) S_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]}|S_k^i|, \quad \text{by Cauchy - Schwarz inequality} \\
&\leq \frac{1}{t} \sum_{k = 1}^t  m \delta^{t-k} \max_{i \in [m]}|S_k^i| \\
&\leq \frac{m}{t} \sum_{k = 1}^t  \sum_{i = 1}^m  \delta^{t-k} |S_k^i| \\
&\leq \frac{m}{t} \sum_{k = t - c + 1}^t  \sum_{i = 1}^m  \delta^{t-k} |S_k^i|  +  \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m |S_k^i|,
\end{align*}

for some fixed constant $c$.

By the strong law of large number,  $\left\|\sum_{k = 1}^t \frac{1}{m} \ones \ones^T S_k - \frac{1}{m}  \sum_{i = 1}^m\E[S_1^i]\right\| = o_d(1)$.  Then

\begin{align*}
\left\|Z_t  - \frac{1}{m}  \sum_{i = 1}^m\E[S_1^i]\right\|  &\leq  \left\|Z_t  - \sum_{k = 1}^t \frac{1}{m} \ones \ones^T S_k \right\|  + \left\|\sum_{k = 1}^t \frac{1}{m} \ones \ones^T S_k  - \frac{1}{m}  \sum_{i = 1}^m\E[S_1^i]\right\| \\
&=   \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m   |S_k^i| +\frac{m}{t} \sum_{k = t - c + 1}^t  \sum_{i = 1}^m  \delta^{t-k} |S_k^i| +o_d(1) \\
&\leq   \frac{m}{t} \delta^c \sum_{k = 1}^{t - c}  \sum_{i = 1}^m   |S_k^i| +o_d(1) \\
&\leq   m \delta^c  \sum_{i = 1}^m\E[S_1^i]+o_d(1).
\end{align*}

Since this for arbitrarily large $c$, we conclude that $ \left|Z_t  - \frac{1}{m}  \sum_{i = 1}^m\E[S_1^i]\right| = o_d(1)$, that is

\begin{equation*}
Z_t \to \frac{1}{m}  \sum_{i = 1}^m\E[S_1^i] \quad \text{in probability}.
\end{equation*}

Assume that $S_t^i, t \geq 1$ are i.i.d. random variables with $\E[S_t^i] = \mu^i$ and $cov[S_t^i] = \Sigma^i$ exists and are finite, for all $ i \in [m]$. Under Assumption <ref>, for any $j \in [m]$, the random variables

\begin{equation*}
Z_t^j = \sqrt{\frac{m}{t}}\left(\frac{1}{m}\sum_{i = 1}^m \Sigma^i \right)^{-1/2} \sum_{i = 1}^m\sum_{k = 1}^t [A_{ij}^{t-k}](S_t^i - \mu^i)
\end{equation*}

converge in distribution to a standard normal distribution as $t \to \infty$.

The proof applies the Linderberg central limit theorem. Without loss of generality, assume that $\mu^i = 0$. Let $\Sigma_{m,t}$ be the sum of covariance. We have

\begin{align*}
\Sigma_{m,t} &= \sum_{i = 1}^m \sum_{k = 1}^t \text{Cov}([A_{ij}^{t-k}]S_t^i) =  \sum_{i = 1}^m \sum_{k = 1}^t[A_{ij}^{t-k}]^2 \Sigma^i .
\end{align*}

The matrix $\Sigma_{m,t}$ is asymptotically equivalent to the scaling factor $\Sigma_{m,t}^* = \frac{t}{m^2} \sum_{i = 1}^m  \Sigma^i$.

\begin{align*}
\left(\Sigma_{m,t}^* \right)^{-1}\Sigma_{m,t} &= \frac{m}{t}  \left(\frac{1}{m} \sum_{i = 1}^m  \Sigma^i\right)^{-1} \sum_{i = 1}^m \sum_{k = 1}^t[A_{ij}^{t-k}]^2 \Sigma^i  \\
&=m  \left(\frac{1}{m} \sum_{i = 1}^m  \Sigma^i\right)^{-1} \sum_{i = 1}^m \Sigma^i \frac{\sum_{k = 1}^t[A_{ij}^{t-k}]^2 }{t} \\
&= m  \left(\frac{1}{m} \sum_{i = 1}^m  \Sigma^i\right)^{-1} \sum_{i = 1}^m \Sigma^i \frac{\sum_{k = 1}^t[A_{ij}^k]^2 }{t}
\end{align*}

By Assumption <ref>, $\lim_{t \to \infty}A_{ij}^k \to \frac{1}{m}$, thus $\lim_{t \to \infty}[A_{ij}^k]^2 \to \frac{1}{m^2}$.Then we have

\begin{align*}
\lim_{t \to \infty}   \left(\Sigma_{m,t}^* \right)^{-1}\Sigma_{m,t} &= m  \left(\frac{1}{m} \sum_{i = 1}^m  \Sigma^i\right)^{-1} \sum_{i = 1}^m \Sigma^i  \lim_{t \to \infty} \frac{\sum_{k = 1}^t[A_{ij}^k]^2 }{t} \\
&=  m  \left(\frac{1}{m} \sum_{i = 1}^m  \Sigma^i\right)^{-1} \sum_{i = 1}^m \Sigma^i  \frac{1}{m^2} = I_m.
\end{align*}

It remains to verify the Linderberg condition

\begin{equation*}
\lim_{t \to \infty} \sum_{i = 1}^m \sum_{k = 1}^t \E\left[\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i \|^2 I_{\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i\| > \epsilon}\right]  = 0, \quad \forall \epsilon > 0.
\end{equation*}

The condition is equivalent to

\begin{equation*}
\lim_{t \to \infty} \sum_{i = 1}^m \sum_{k = 1}^t \E\left[\|\Sigma_{m,t}^{*^{-1/2}} [A_{ij}^{t-k}] S_t^i\|^2 I_{\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i\| > \epsilon}\right]  = 0, \quad \forall \epsilon > 0.
\end{equation*}

Since $\lim_{t \to \infty}\Sigma^*_{m,t} =  \infty$, $\lim_{t \to \infty}\Sigma_{m,t} = \infty$ and $ I_{\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i\| > \epsilon}$ converges to $0$ almost surely.

Because $\max_{i \in [m]} \sup_t [A_{ij}^t]^2 < \infty$, we have

\begin{align*}
\max_{i \in [m]}\sup_{t} t\E\left[ \|\Sigma_{m,t}^{*^{-1/2}} [A_{ij}^{t-k}] S_t^i\|^2\right] &= \max_{i \in [m]} \sup_{t} \E\left[\|m(\frac{1}{m}\sum_{i = 1}^m  \Sigma^i)^{-1} [A_{ij}^{t-k}] S_t^i\|^2\right]  \\
&\leq m (\max_{i \in [m]} \sup_t [A_{ij}^t]^2 ) \|(\frac{1}{m}\sum_{i = 1}^m  \Sigma^i)^{-1}\|_{op} \max_{i \in [m]} \|\Sigma^i \|_{op}
\end{align*}

By Lebesgue dominated convergence theorem, we have

\begin{equation*}
\lim_{t \to \infty} t\E\left[\|\Sigma_{m,t}^{*^{-1/2}} [A_{ij}^{t-k}] S_t^i\|^2 I_{\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i\| > \epsilon}\right]  \to 0.
\end{equation*}

The Cesaro sums also converge to $0$:

\begin{equation*}
\lim_{t \to \infty} \sum_{i = 1}^m \frac{1}{t}\sum_{k = 1}^t t \E\left[\|\Sigma_{m,t}^{*^{-1/2}} [A_{ij}^{t-k}] S_t^i\|^2 I_{\|\Sigma_{m,t}^{-1/2} [A_{ij}^{t-k}] S_t^i\| > \epsilon}\right]  = 0, \quad \forall \epsilon > 0.
\end{equation*}

This completes the proof

Suppose that $h_n: E \to F$ for $n \in \N$, where $E$ is a totally bounded space and $F$ is a normed space. If $h_n$ converges pointwise and is equicontinuous, then it converges uniformly.

Suppose that the posterior distribution $\Pi_n(\cdot \mid X(n))$ is consistent (or strongly consistent) at $\theta_0$ relative to the metric $d$ on $\Theta$. Then $\hat{\theta}_n$, defined as the center of a (nearly) smallest ball that contains posterior mass at least $\delta$ for $\delta > 0$, satisfies

\begin{equation*}
d(\hat{\theta}_n, \theta_0) \rightarrow 0,
\end{equation*}

in $[P_{\theta_0}]-$probability or $[P_{\theta_0}]-$almost surely, respectively.


# LRGNet: Learnable Region Growing for Class-Agnostic Point Cloud Segmentation

Jingdao Chen1 , Zsolt Kira1and Yong K. Cho1 1 Institute for Robotics and
Intelligent Machines, Georgia Institute of Technology, Atlanta, GA 30332, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

3D point cloud segmentation is an important function that helps robots
understand the layout of their surrounding environment and perform tasks such
as grasping objects, avoiding obstacles, and finding landmarks. Current
segmentation methods are mostly class-specific, many of which are tuned to
work with specific object categories and may not be generalizable to different
types of scenes. This research proposes a learnable region growing method for
class-agnostic point cloud segmentation, specifically for the task of instance
label prediction. The proposed method is able to segment any class of objects
using a single deep neural network without any assumptions about their shapes
and sizes. The deep neural network is trained to predict how to add or remove
points from a point cloud region to morph it into incrementally more complete
regions of an object instance. Segmentation results on the S3DIS and ScanNet
datasets show that the proposed method outperforms competing methods by 1%-9%
on 6 different evaluation metrics.

##  I INTRODUCTION

A key step in parsing a 3D point cloud scene is point cloud segmentation.
Segmentation is useful for many robotic applications such as obstacle
detection [1], localization [2], and object recognition [3][4]. In this
research, point cloud segmentation is defined as subdividing a point cloud
scene into individual object segments where each point in the point cloud is
assigned to a unique instance label (i.e. class-agnostic segmentation). One
popular type of approach is to perform region growing to iteratively create
different point cloud segments [5][6][7]. However, this approach is highly
sensitive to specific threshold hyperparameters defined over the raw feature
values. Another popular type of approach is to directly learn the segmentation
labels using a deep neural network [8][9][10]. These approaches are data-
driven and are usually more accurate and robust, but require discretizing the
point cloud scene into 1m x 1m blocks [4][8] and need class annotations during
training. We propose to extend these two types of approaches, combining the
region growing approach with data-driven deep learning techniques in order to
achieve the benefits of both types of approaches.

Thus, this research proposes a point cloud segmentation method where the
region growing process itself is learned from data. Given a raw point cloud
scene, the proposed method tries to form point cloud segments that each
correspond to a complete object by making a series of incremental predictions.
Specifically, a deep neural network is trained to predict at each decision
point: (i) which points to add to grow the current region, (ii) which
unnecessary points to remove from the current region, conditioned on the
current set of points. The proposed method is evaluated on the S3DIS [11] and
ScanNet [12] datasets with a focus on generalization performance across
datasets. The proposed method has the following advantages:

  * •

able to segment any class of objects using a single neural network

  * •

able to segment objects without assumptions about their shapes and sizes

  * •

can be trained without semantic labels

  * •

does not require block division of the point cloud scene

  * •

more accurate and generalizable compared to existing methods

In addition, the code and data have been made publicly available.
111https://github.com/jingdao/learn_region_grow

##  II LITERATURE REVIEW

Point cloud segmentation can take the form of: (i) semantic label prediction -
the task of predicting a class label (e.g. wall, window) for each point of a
given point cloud scene, or (ii) instance label prediction, the task of
assigning points to unique objects instances in a given point cloud scene.
This research will specifically focus instance label prediction, which is
especially challenging in real-world point cloud scenes because they are often
cluttered and contain objects with a wide variety of shapes and sizes.

###  II-A Class-agnostic segmentation

Class-agnostic methods generally perform segmentation without making any
assumptions about the class of objects that are in the point cloud scene.
Class-agnostic methods are generally instance-based segmentation methods. In
the robotics domain, these methods are commonly implemented as region growing,
where groups of points with similar features such as surface normal and color
are incrementally merged together [2][6]. These methods face issues with
oversegmentation or undersegmentation in cluttered scenes due to the reliance
on fixed thresholds or heuristics to control the region growing process. Some
works attempt to use reinforcement learning [13][14] or evolutionary methods
[15] to more effectively tune the threshold hyperparameters, but do not
directly optimize the region growing process itself. An alternative to region
growing is to perform segmentation as a global operation by using clustering
algorithms such as K-means [16] or DBSCAN [17]. The problem with these
clustering algorithms is that the results are sensitive to the hyperparameters
such as number of segments and the neighborhood range. Other ideas for class-
agnostic segmentation include using a network to predict object seeds [18] or
searching over the space of candidate segmentations and returning one which
scores well according to an “objectness” model [19]. However, the former
method was only evaluated for the tabletop setting where most objects are
roughly the same size whereas the latter method was only applied in outdoor
environments where objects do not physically overlap with one another.

###  II-B Class-specific segmentation

Class-specific segmentation methods try to separate objects from the
background based on prior information about the geometry of particular classes
of objects. Class-specific methods are generally semantic segmentation
methods. These methods are usually data-driven and make use of deep neural
networks. A popular approach is to use 3D bounding box regression to predict
the location and size of objects in the scene, used by methods such as MV3D
[20], PointRCNN [10], PointPillars [21], and 3D-BoNet [22]. The downside of
this approach is the need to define class-specific anchor boxes [20] or train
multiple networks separately for each class [23], which is not scalable for a
large number of classes. In addition, bounding boxes only provide a coarse
outline and do not cleanly separate object segments. A different strategy is
to directly predict labels for each point in the point cloud using neural
networks based on PointNet [4] or PointNet++ [24]. These network are usually
used to perform semantic segmentation or part segmentation and do not directly
predict the instance labels. Later methods combine the semantic segmentation
results with grouping matrix prediction [8][25], centroid prediction [26][27],
or clustering on a learnt embedding [9][28] in order to obtain instance
segmentation results. However, a common problem is that they require a pre-
processing step of subdividing the point cloud into fixed-size blocks and a
post-processing step of BlockMerging [8] to combine segments that span
different blocks together, which leads to information loss. Another weakness
is that they rely on having training data that is labeled with class
information and ignore objects that are not within the predefined set of
classes (usually designated as ”clutter”). In contrast to the proposed method
which is class-agnostic, these class-specific segmentation methods are
difficult to generalize across datasets which have different classes.

##  III METHODOLOGY

This section describes the proposed LRGNet method for class-agnostic point
cloud segmentation including the formal definition, network architecture,
training process, inference time algorithm, and local search optimization.

Figure 1: LRGNet network architecture to perform an intermediate step of
region growing. The network takes as input two sets of points: the inlier set
and the neighbor set. Features are computed from both sets of points and then
pooled and concatenated to form a global feature vector. This feature vector
is concatenated back into the point-level features to predict which points in
the inlier set to remove (i.e. remove mask prediction) and which points in the
neighbor set to add (i.e. add mask prediction)

###  III-A Formal definition

The point cloud segmentation problem is defined as follows. The input is a
point cloud P, represented as a $N\times F$ matrix of real numbers where the
input point cloud contains $N$ points with each point having $F$ features
(e.g. 3D position, color, etc.). The task is to output the instance labels L,
represented as a $N$-dimensional vector of integers, which contains one
integer label for each point indicating which object instance it should belong
to.

The segmentation problem can be decomposed into several region-growing
subproblems. The region-growing procedure takes in an input seed point
$p_{seed}\in\textbf{P}$, and outputs a set of points $P^{*}\subset P$, that
all belong to the same object instance as $p_{seed}$. The instance label
vector, L, is then updated such that all points in $P^{*}$ have the same
instance label. The idea is that this region-growing procedure can be carried
out repeatedly from any remaining unlabeled seed points in the point cloud P
until all points have been assigned an instance label.

In this study, the region-growing subproblem is further decomposed into a
series of steps, where each step takes in an input set of points $Q_{k}$ and
outputs a set of points $Q_{k+1}$, where $Q_{k},Q_{k+1}\subset\textbf{P}$.
This can be represented as $f(Q_{k})=Q_{k+1}$, where $f$ is a learnable
function. The region-growing process starts with $Q_{0}=\\{p_{seed}\\}$ and
computes $Q_{1},Q_{2},Q_{3},...$ etc. until $Q\rightarrow P^{*}$, the set of
all points which belong to the same object instance as $p_{seed}$. This study
aims to learn the function $f$ which dictates how a point cloud region grows
from one set of points to another. The key novelty of this formulation is to
decompose the region growing problem and use a learnable function to carry it
out, thus achieving a ”learnable” region growing method.

###  III-B Network architecture

The learnable function $f(Q_{k})=Q_{k+1}$ takes the form of a deep neural
network named LRGNet, as shown in Figure 1. The input to the network consists
of two sets of points, the inlier set (represented as an $I\times F$ matrix)
and the neighbor set (represented as a $J\times F$ matrix). The inlier set is
constructed by taking $I$ points from $Q_{k}$. Note that the full set of
points in the current region, $Q_{k}$, is maintained as a separate global
variable. The neighbor set is constructed by taking all the unlabeled points
that are within a predetermined $\delta$ distance of any point in $Q_{k}$ and
selecting $J$ points. When there are more or fewer than $I$ or $J$ points in
the subsets, random resampling with replacement is used to obtain exactly $I$
or $J$ points.

The network uses an encoder-decoder structure. First, the input inlier set and
neighbor set are separately processed by individual network branches each
containing 5 point-wise feature convolution layers. The features across all
the points are then aggregated using a max-pooling layer to obtain a global
feature vector incorporating 512 features from the inlier set and 512 features
from the neighbor set. This global feature vector is then separately
concatenated back into two decoder branches, one processing the inlier set and
one processing the neighbor set. These two branches are then further processed
by point-wise feature convolution layers. This network design is largely based
off of PointNet [4], but with two input branches and two output branches
instead of a single input branch and a single output branch.

The upper branch finally outputs an $I$-dimensional vector which constitutes
the remove mask probability prediction output. This is then used to determine
the remove mask, a binary mask which determines which points in the inlier set
should be removed from $Q_{k}$. On the other hand, the lower branch finally
outputs a $J$-dimensional vector which constitutes the add mask probability
prediction output. This is then used to determine the add mask, a binary mask
which determines which points in the neighbor set should be added to $Q_{k}$.
The reason for having two separate masks is because the remove mask prediction
allows the network to correct for previous errors and makes it more robust to
noise. Finally, the new set of points $Q_{k+1}$ is determined by applying the
remove mask to the inlier set and removing the corresponding points from
$Q_{k}$, then applying the add mask to the neighbor set and adding those
points to $Q_{k}$.

This study uses a total of $F=13$ input features for each point, which are the
local XYZ-coordinates, room-normalized XYZ-coordinates, RGB-color, normal
vector, and curvature. The room-normalized XYZ-coordinates are obtained by
dividing the X,Y,Z coordinates by the length, width, height of the room
respectively to obtain values ranging from 0 to 1. The local XYZ-coordinates
describe a point in relation to its object center whereas the room-normalized
XYZ-coordinates describe a point in relation to the global room scene. The
normal vector and curvature are computed using the Principal Component
Analysis method from [6] whereas the XYZ features are computed as in [4].
These point features (except the room-normalized XYZ-coordinates) are
normalized column-by-column by subtracting the median over all the points in
the inlier set. On the other hand, the parameters $\delta$ (neighbor threshold
distance), $I$ (inlier set size), $J$ (neighbor set size) are set to 0.1m,
512, 512 respectively. These parameters are tuned based on an ablation study
on the S3DIS dataset.

The network is trained on the binary cross entropy loss on both the add mask
prediction and the remove mask prediction. This is shown in (1), where
$\hat{x_{i}}$ denotes the predicted probability for removing point $i$ in the
inlier set and $x_{i}$ is the corresponding ground truth label, whereas
$\hat{y_{j}}$ denotes the predicted probability for adding point $j$ in the
neighbor set and $y_{j}$ is the corresponding ground truth label.

$\begin{split}\mathcal{L}=-\frac{1}{I}\sum^{I}_{i}{[x_{i}\log{\hat{x_{i}}}+(1-x_{i})\log{(1-\hat{x_{i}})}]}\\\
-\frac{1}{J}\sum^{J}_{j}{[y_{j}\log{\hat{y_{j}}}+(1-y_{j})\log{(1-\hat{y_{j}})}]}\end{split}$
(1)

Figure 2: Visualization of four different intermediate steps of the region
growing process starting from a seed point (points in the inlier set are shown
in green whereas points in the neighbor set are shown in blue).

###  III-C Region growing simulation and training

Input: PointCloud, P

Output: InstanceLabels, L

$\phi$ $\equiv$ empty label ;

$L=\\{\phi,\forall p\in P\\}$ ;

$obj\\_id=1$;

while  _$\phi\notin L$_ do

$p_{seed}$ = unlabeled point with min curvature;

$k=0$;

$Q_{k}=\\{p_{seed}\\}$;

while  _true_ do

$inliers=sample(Q_{k},count=I)$;

$neighbors=sample(\\{p:|p-q|<\delta,p\in P,q\in Q_{k},L[p]=\phi\\},count=J)$;

$add,remove=net(inliers,neighbors)$;

if  _termination condition_ then

$L[Q_{k}]=obj\\_id$;

$obj\\_id=obj\\_id+1$;

break;

else

$Q_{k+1}=Q_{k}-inliers[remove]+neighbors[add]$;

$k=k+1$;

Algorithm 1 Region Growing

In this research, training data is obtained by carrying out a series of region
growing simulations from labeled point cloud datasets, S3DIS [11] and ScanNet
[12]. The region growing process for each object instance in a point cloud
scene is simulated in a series of steps, starting from a randomly selected
seed point. At each step, the desired next point set, $Q_{k+1}$, is created
from the current point set, $Q_{k}$, by adding all the points that are within
$\delta$ distance of any point in $Q_{k}$ and also have the same ground truth
instance label. To prevent the network from overfitting to only clean point
cloud samples, noise is artificially introduced into the region growing
process by randomly adding and removing points to $Q_{k+1}$ according to a
mistake probability, $\alpha$. Within one region growing simulation, the
mistake probability starts at a high value and is then gradually decreased by
0.01 each step for the current region to ensure that region growing process
will eventually converge to a complete object instance.

Data augmentation is carried out in several ways. First, each point cloud
scene is randomly flipped to swap the x and y axes. Next, the point cloud
scene is randomly rotated in increments of 90∘. Also, the mistake probability,
$\alpha$, for each object instance is randomly assigned an initial value
between 0.2 and 0.4 to simulate different levels of noise in the region
growing process. For example, after data augmentation on the S3DIS dataset,
the training set consists of 1406516 point clouds originating from 63047
object instances.

The acquired training data is used to train LRGNet as described in Section
IIIB, implemented in Tensorflow. Note that the training process only uses the
ground truth instance labels and does not require semantic labels. The network
is trained with the ADAM optimizer until convergence at around 40 epochs. The
learning rate used is 0.001 and the batch size is 100.

###  III-D Region growing at inference time

At inference time, segmentation of a point cloud scene is carried out by
continuously performing region growing until each point in the scene is
assigned an instance label. Starting from a seed point, the trained neural
network is applied to predict the next region to grow to using the add mask
and remove mask output. This process is repeated until one of the following
termination conditions is met:

  * •

No unassigned neighbor points remaining

  * •

The set of points to be added is predicted to be empty

  * •

Region does not expand for two consecutive steps

The last termination condition is necessary to avoid situations where
different points are alternatingly added then removed on consecutive steps in
an alternating manner, which could occur near edges. Upon reaching the
termination condition, all the points in the current region areis assigned the
same instance label and the region growing process starts over again from a
new seed point. Note that the region points may change over the course of
region growing for a single object instance but the assigned instance label
will be permanent once the region growing process switches to a new object.

Since the region growing result could be variable due to the choice of the
seed point, this study uses the strategy of picking the remaining point with
minimum curvature as the seed point. This strategy is also used in [29]. In
addition, as the region growing progresses, certain groups of points may
become isolated. This frequently occurs for edges in between two object
instances. To overcome this, if the resulting point cloud segment is too small
(e.g. less than 10 points), then a new instance is not created. Instead, the
points in that isolated segment are assigned instance labels according to the
nearest neighbor among points that have already been assigned instance labels.
A similar strategy is used in [4][30].

Algorithm 1 shows a brief pseudocode summary of the region growing process for
a given point cloud scene. To further illustrate the process, Figure 2 shows a
visualization of region growing for a wall object. As shown in Figure 2, there
tends to be a higher number of added and removed points during the middle of
the process and a lower number of added and removed points during the start
and end of the process. On average, the percentage of added points and removed
points per step are 70% and 59% respectively for a typical point cloud.

###  III-E Local search to optimize region growing

Since the learned function $f(Q_{k})=Q_{k+1}$ is not perfect and also depends
on randomly sampling the inlier set and neighbor set, the region growing
process may not converge on the perfect point cloud segment that contains all
the points from one object instance. To overcome this, local search techniques
can be applied to search through different sequences of applying the learned
region growing function to find the best final region. This study considers 5
different methods of arriving at the final region from a given seed point: (i)
greedy: the add and remove outputs are greedily applied as binary masks with a
cutoff threshold of 0.5, (ii) random restart - Maximum Likelihood (ML): the
region growing process is restarted from the same seed point multiple times,
each time probabilistically sampling points to be added and removed and
keeping track of the sum of the log-likelihoods from each step; the best final
region is determined using maximum likelihood as the optimality criterion,
(iii) random restart - Number of Points (NP): same as the previous method but
with the number of points in the final region as the optimality criterion
(based on the assumption that regions with more points are more complete),
(iv) beam search - ML: multiple states for the current region are
simultaneously tracked and expanded [31]; after each expansion step, the
number of states is reduced to the $K$ best states using maximum likelihood as
the optimality criterion, (v) beam search - NP: same as the previous method
but with the number of points as the optimality criterion. A comparison of the
segmentation performance between different local search methods will be given
in the Results section.

##  IV RESULTS

###  IV-A Segmentation in indoor scenes

TABLE I: Comparison of segmentation performance on S3DIS and ScanNet by testing generalization across datasets Method | NMI | AMI | ARI | PRC | RCL | mIOU
---|---|---|---|---|---|---
ScanNet (train) $\rightarrow$ S3DIS (test)
Region growing | 0.71 | 0.70 | 0.59 | 0.19 | 0.34 | 0.38
Rabbani et al. [32] | 0.72 | 0.71 | 0.62 | 0.17 | 0.31 | 0.36
FPFH [33] | 0.62 | 0.60 | 0.39 | 0.14 | 0.25 | 0.32
PointNet [4] | 0.58 | 0.48 | 0.38 | 0.18 | 0.17 | 0.25
PointNet++ [24] | 0.62 | 0.56 | 0.40 | 0.15 | 0.22 | 0.31
JSIS3D [9] | 0.74 | 0.73 | 0.63 | 0.28 | 0.29 | 0.36
3D-BoNet [10] | 0.75 | 0.72 | 0.68 | 0.20 | 0.29 | 0.35
LRGNet | 0.75 | 0.74 | 0.67 | 0.25 | 0.41 | 0.43
LRGNet + local search | 0.76 | 0.75 | 0.68 | 0.34 | 0.44 | 0.45
S3DIS (train) $\rightarrow$ ScanNet (test)
Region growing | 0.62 | 0.60 | 0.44 | 0.17 | 0.23 | 0.30
Rabbani et al. [32] | 0.64 | 0.62 | 0.49 | 0.13 | 0.24 | 0.32
FPFH [33] | 0.53 | 0.51 | 0.28 | 0.10 | 0.14 | 0.26
PointNet [4] | 0.57 | 0.51 | 0.40 | 0.08 | 0.13 | 0.26
PointNet++ [24] | 0.63 | 0.57 | 0.47 | 0.15 | 0.21 | 0.32
JSIS3D [9] | 0.57 | 0.56 | 0.31 | 0.15 | 0.13 | 0.22
3D-BoNet [10] | 0.59 | 0.54 | 0.34 | 0.10 | 0.13 | 0.24
LRGNet | 0.69 | 0.67 | 0.54 | 0.25 | 0.33 | 0.39
LRGNet + local search | 0.69 | 0.68 | 0.56 | 0.31 | 0.33 | 0.38

The segmentation performance in indoor scenes was evaluated on the S3DIS
dataset [11] and the ScanNet dataset [12] alternatingly with one dataset used
as training data while the other dataset is held out as test data. The first
three evaluation metrics used are from the clustering literature, namely
normalized mutual information (NMI), adjusted mutual information (AMI), and
adjusted rand index (ARI), as defined in [34]. The other three evaluation
metrics used are from the object detection literature, namely precision (PRC),
recall (RCL), and mean intersection-over-union (mIOU). For the purposes of
calculating precision and recall, true positive detections are defined as
cases where the IOU between a predicted point cloud segment and a ground truth
point cloud segment is greater than 50% (as established in [11]). Each
evaluation metric is calculated per room of the test dataset and averaged
across all rooms (note that this setup means that they differ from the metrics
used in [9][10] which are averaged across classes). The proposed region
growing method is compared against several class-agnostic competing methods
including (i) conventional region growing which uses the raw features directly
with thresholding, (ii) Rabbani et al. which uses smoothness constraints [35],
(iii) segmentation using the Fast Point Feature Histogram (FPFH) descriptor
[33]. The class-specific competing methods are (iv) PointNet [4] (v)
PointNet++ [24] (vi) JSIS3D [9] and (vii) 3D-BONET [10]. The latter four
methods are implemented based on publicly-available open-source code.

Table I shows a quantitative comparison of each evaluation metric between the
different segmentation methods. Results show that the proposed LRGNet achieves
the overall best results across all metrics whether tested on the S3DIS
dataset or the ScanNet dataset. In addition, the results also show that
applying a local search method (random restart with number of points as the
optimality criterion) can further improve the results by optimizing the region
growing process of LRGNet. Note that PointNet and PointNet++ do not directly
compute the instance labels. Instead, instance labels are obtained by
clustering neighboring points that have the same semantic label, as described
in [4]. This process could cause inaccuracies which lead to worse performance
in Table 1 when compared to other methods such as region growing or FPFH.
Furthermore, this generalization study also demonstrates the importance of
having a class-agnostic point cloud segmentation method. ScanNet is collected
in mostly home environments whereas S3DIS is collected in mostly office
environments which contain different classes of objects. As a result, most
methods that are trained with class information on the S3DIS dataset do not
perform as well on the ScanNet dataset. On the other hand, class-agnostic
methods such as the proposed LRGNet are able to better generalize to the new
point cloud scenes.

Figures 3 and 4 show a qualitative comparison of the segmentation results from
an overhead perspective, with the ceiling removed to improve visualization.
Results show that competing methods have more cases of undersegmentation or
oversegmentation, whereas the proposed method is able to perform better since
it is directly trained on optimizing the region growing segmentation process.

Figure 3: Visualization of segmentation results on the S3DIS dataset: (from
left to right) (i) original RGB point cloud (ii) ground truth (iii) PointNet++
(iv) 3D-BoNet (v) proposed method  Figure 4: Visualization of segmentation
results on the ScanNet dataset: (from left to right) (i) original RGB point
cloud (ii) ground truth (iii) PointNet++ (iv) 3D-BoNet (v) proposed method

###  IV-B Segmentation in outdoor scenes

The segmentation performance in outdoor scenes is evaluated on the Semantic
KITTI dataset [36]. To be consistent with the experimental setup in indoor
scenes, the original Semantic KITTI dataset is modified as follows: (i) RGB
channels are projected from the left camera images to the LiDAR scans to get
color-mapped point clouds, (ii) the semantic labels are discarded and only
instance labels are used during training, (iii) points without a valid RGB
color and instance label are filtered out, (iv) scan data from 20 consecutive
LiDAR scans are merged to get a dense point cloud. 10 driving sequences are
used as training data whereas 1 driving sequence is held out as test data.
Figure 5 shows a qualitative comparison between the proposed method and other
class-agnostic segmentation methods such as simple region growing and FPFH.
Results show that the proposed method is able to better separate out the road
from sidewalk areas compared to other methods which tend to undersegment those
areas. The proposed method is also able to generate more complete segments of
building and vegetation compared to other methods which tend to oversegment
those objects.

Figure 5: Visualization of segmentation results on the Semantic KITTI dataset:
(from left to right) (i) original RGB point cloud (ii) ground truth (iii)
simple region growing (iv) FPFH (v) proposed method

###  IV-C Ablation studies

Table II shows a series of ablation studies of the segmentation performance on
the S3DIS dataset for different architecture design choices. In these ablation
studies, Areas 1,2,3,4,6 of S3DIS are used as the training set whereas Area 5
of S3DIS is used as the validation set. These ablation studies include using
the original method but (i) without using the remove mask prediction to remove
outlier points from the current region, (ii) without using the minimum
curvature criteria when selecting seed points, instead selecting them
randomly, (iii) without normalizing the input point features, (iv) only using
XYZ features, (v) only using XYZ+RGB features, (vi) using $I=128$, $J=128$, or
(vii) using $I=256$, $J=256$. Results show that these variants lead to a 5-10%
change in performance but the complete method has the best performance
overall.

TABLE II: Ablation study of different design choices on S3DIS dataset Method | NMI | AMI | ARI | PRC | RCL | mIOU
---|---|---|---|---|---|---
no remove mask | 0.78 | 0.75 | 0.73 | 0.51 | 0.41 | 0.42
no seed selection | 0.81 | 0.77 | 0.76 | 0.40 | 0.55 | 0.53
no feature normalization | 0.81 | 0.76 | 0.75 | 0.39 | 0.53 | 0.52
only XYZ | 0.76 | 0.75 | 0.61 | 0.24 | 0.42 | 0.44
only XYZ+RGB | 0.78 | 0.78 | 0.70 | 0.33 | 0.53 | 0.51
$I=128$, $J=128$ | 0.79 | 0.79 | 0.69 | 0.32 | 0.57 | 0.55
$I=256$, $J=256$ | 0.80 | 0.80 | 0.74 | 0.35 | 0.56 | 0.54
complete method | 0.81 | 0.78 | 0.77 | 0.43 | 0.56 | 0.54

Table III shows the segmentation performance on the S3DIS dataset when
comparing different local search methods based on the quality (NMI, AMI, ARS
metrics) as well as efficiency (average number of inference steps needed to
segment one object instance). The NMI, AMI, and ARS metrics are calculated as
mean $\pm$ standard deviation across all rooms of the S3DIS dataset. Random
restart is applied with a total of 10 restarts whereas beam search is applied
with 3 candidates each expanded 3 times. This setting is so that random
restart and beam search will have roughly the same number of computation
steps. Note that the non-greedy methods use a noisier randomized selection
criteria compared to the greedy method, so the number of computational steps
is not exactly proportional to the number of restarts (or beam width).
Overall, both random restart and beam search improved the segmentation
performance compared to the greedy method, with random restart - NP having the
best performance overall. However, this comes with the computational tradeoff
of having many more passes through the prediction network for each object
instance.

TABLE III: Segmentation performance on S3DIS dataset with different local search methods Method | NMI(%) | AMI(%) | ARS(%) | Avg steps
---|---|---|---|---
greedy | 82$\pm$4 | 73$\pm$7 | 74$\pm$12 | 13.38
random restart - ML | 82$\pm$5 | 78$\pm$7 | 77$\pm$12 | 190.06
random restart - NP | 82$\pm$4 | 79$\pm$6 | 77$\pm$10 | 188.70
beam search - ML | 82$\pm$5 | 78$\pm$7 | 77$\pm$12 | 159.79
beam search - NP | 82$\pm$4 | 78$\pm$6 | 77$\pm$10 | 175.50

###  IV-D Computation time analysis

Table IV shows a comparison of the computation time involved to completely
segment a point cloud scene, evaluated on a randomly selected subset of 50
scenes from the S3DIS dataset. The computation time is measured by the authors
on an Intel Xeon E3-1200 CPU with a NVIDIA GTX1080 GPU. Compared to other
methods, LRGNet lies on the higher end of the spectrum in terms of computation
cost due to having to run a large number of region growing iterations to
segment each scene. In future work, using parallel processing methods or an
adaptive region growing step size could potentially improve the computational
speed.

TABLE IV: Computation time comparison for segmentation per scene Method | Computation time (s)
---|---
| Min | Mean | Max
Region growing | 0.4 | 4.8 | 18.6
Rabbani et al. | 0.3 | 4.6 | 18.3
FPFH | 0.5 | 4.6 | 17.8
Pointnet | 0.1 | 0.6 | 2.5
PointNet++ | 0.1 | 0.9 | 3.5
JSIS3D | 1.0 | 539.2 | 16713.9
3D-BoNet | 1.5 | 14.1 | 69.3
LRGNet | 0.8 | 64.9 | 620.9

##  V CONCLUSIONS

3D point cloud segmentation is important for robots to be able to understand
the layout and geometry of objects in the surrounding environment. To achieve
this, this paper proposes a learnable region growing method that can perform
class-agnostic point cloud segmentation. The proposed method works by
decomposing the segmentation problem into a series of region growing
subproblems. A deep neural network is trained to predict how to add or remove
points from a point cloud region to morph it into incrementally more complete
regions of an object instance. The advantages of the proposed method are that
is can be used to segment objects of arbitrary shape, size, or class, and that
it is more accurate and generalizable compared to existing methods.

## ACKNOWLEDGMENT

The work reported herein was supported by the United States Air Force Office
of Scientific Research (Award No. FA2386-17-1-4655) and by a grant
(18CTAP-C144787-01) funded by the Ministry of Land, Infrastructure, and
Transport (MOLIT) of Korea Agency for Infrastructure Technology Advancement
(KAIA). Any opinions, findings, and conclusions or recommendations expressed
in this material are those of the authors and do not necessarily reflect the
views of the United States Air Force and MOLIT. Special thanks to John Yi Seon
Keun for assisting with the segmentation experiments.

## References

  * [1] J. Chen, P. Kim, Y. K. Cho, and J. Ueda, “Object-sensitive potential fields for mobile robot navigation and mapping in indoor environments,” in _2018 15th International Conference on Ubiquitous Robots (UR)_ , June 2018\.
  * [2] R. Dube, M. G. Gollub, H. Sommer, I. Gilitschenski, R. Siegwart, C. Cadena, and J. Nieto, “Incremental-segment-based localization in 3-d point clouds,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, July 2018.
  * [3] P. Kim, J. Chen, and Y. K. Cho, “Robotic sensing and object recognition from thermal-mapped point clouds,” _International Journal of Intelligent Robotics and Applications_ , vol. 1, no. 3, Sep 2017.
  * [4] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” _Proc. Computer Vision and Pattern Recognition (CVPR), IEEE_ , 2017.
  * [5] A. C. C. Law, N. Southon, N. Senin, P. Stavroulakis, R. Leach, R. D. Goodridge, and Z. Kong, “Curvature-based segmentation of powder bed point clouds for inprocess monitoring,” in _Proceedings of the 29th Annual International Solid Freeform Fabrication Symposium_ , 2018.
  * [6] R. B. Rusu, “Semantic 3d object maps for everyday manipulation in human living environments,” _KI - Künstliche Intelligenz_ , vol. 24, no. 4, 2010.
  * [7] X. Wang, L. Zou, X. Shen, Y. Ren, and Y. Qin, “A region-growing approach for automatic outcrop fracture extraction from a three-dimensional point cloud,” _Computers and Geosciences_ , vol. 99, 2017.
  * [8] W. Wang, R. Yu, Q. Huang, and U. Neumann, “Sgpn: Similarity group proposal network for 3d point cloud instance segmentation,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018.
  * [9] Q.-H. Pham, T. Nguyen, B.-S. Hua, G. Roig, and S.-K. Yeung, “Jsis3d: Joint semantic-instance segmentation of 3d point clouds with multi-task pointwise networks and multi-value conditional random fields,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * [10] B. Yang, J. Wang, R. Clark, Q. Hu, S. Wang, A. Markham, and N. Trigoni, “Learning object bounding boxes for 3d instance segmentation on point clouds,” in _Advances in Neural Information Processing Systems 32_ , 2019\.
  * [11] I. Armeni, O. Sener, A. R. Zamir, H. Jiang, I. Brilakis, M. Fischer, and S. Savarese, “3d semantic parsing of large-scale indoor spaces,” in _Proceedings of the IEEE International Conference on Computer Vision and Pattern Recognition_ , 2016.
  * [12] A. Dai, A. X. Chang, M. Savva, M. Halber, T. Funkhouser, and M. Nießner, “Scannet: Richly-annotated 3d reconstructions of indoor scenes,” in _Proc. Computer Vision and Pattern Recognition (CVPR), IEEE_ , 2017.
  * [13] F. Sahba, H. R. Tizhoosh, and M. M. A. Salama, “A reinforcement learning framework for medical image segmentation,” in _The 2006 IEEE International Joint Conference on Neural Network Proceedings_ , July 2006.
  * [14] M. Bhardwaj, S. Choudhury, and S. A. Scherer, “Learning heuristic search via imitation,” in _1st Annual Conference on Robot Learning, CoRL 2017, Mountain View, California, USA, November 13-15, 2017, Proceedings_ , ser. Proceedings of Machine Learning Research, vol. 78. PMLR, 2017.
  * [15] S. Derivaux, G. Forestier, C. Wemmert, and S. Lefèvre, “Supervised image segmentation using watershed transform, fuzzy classification and evolutionary computation,” _Pattern Recognition Letters_ , vol. 31, no. 15, 2010.
  * [16] X. X. Zhu and M. Shahzad, “Facade reconstruction using multiview spaceborne tomosar point clouds,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 52, no. 6, June 2014.
  * [17] S. M. Ahmed and C. M. Chew, “Density-based clustering for 3d object detection in point clouds,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020.
  * [18] C. Xie, Y. Xiang, A. Mousavian, and D. Fox, “Unseen object instance segmentation for robotic environments,” in _arXiv:2007.08073_ , 2020.
  * [19] P. Hu, D. Held, and D. Ramanan, “Learning to optimally segment point clouds,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, April 2020\.
  * [20] X. Chen, H. Ma, J. Wan, B. Li, and T. Xia, “Multi-view 3d object detection network for autonomous driving,” in _IEEE CVPR_ , 2017.
  * [21] A. H. Lang, S. Vora, H. Caesar, L. Zhou, J. Yang, and O. Beijbom, “Pointpillars: Fast encoders for object detection from point clouds,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * [22] S. Shi, X. Wang, and H. Li, “Pointrcnn: 3d object proposal generation and detection from point cloud,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * [23] F. Liu, S. Li, L. Zhang, C. Zhou, R. Ye, Y. Wang, and J. Lu, “3dcnn-dqn-rnn: A deep reinforcement learning framework for semantic parsing of large-scale 3d point clouds,” in _The IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017.
  * [24] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,” _arXiv preprint arXiv:1706.02413_ , 2017.
  * [25] H. Jiang, F. Yan, J. Cai, J. Zheng, and J. Xiao, “End-to-end 3d point cloud instance segmentation without detection,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020\.
  * [26] L. Jiang, H. Zhao, S. Shi, S. Liu, C.-W. Fu, and J. Jia, “Pointgroup: Dual-set point grouping for 3d instance segmentation,” _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [27] F. Engelmann, M. Bokeloh, A. Fathi, B. Leibe, and M. Niessner, “3d-mpa: Multi-proposal aggregation for 3d semantic instance segmentation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020.
  * [28] T. H. He, Y. Liu, C. Shen, W. Xinlong, and C. Sun, “Instance-aware embedding for point cloud instance segmentation,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2020.
  * [29] A. Dimitrov and M. Golparvar-Fard, “Segmentation of building point cloud models including detailed architectural/structural features and mep systems,” _Automation in Construction_ , vol. 51, 2015.
  * [30] J. Chen, Z. Kira, and Y. K. Cho, “Deep learning approach to point cloud scene understanding for automated scan to 3d reconstruction,” _Journal of Computing in Civil Engineering_ , vol. 33, no. 4, 2019.
  * [31] R. Haeb-Umbach and H. Ney, “Improvements in beam search for 10000-word continuous-speech recognition,” _IEEE Transactions on Speech and Audio Processing_ , vol. 2, no. 2, April 1994.
  * [32] T. Rabbani, F. van den Heuvel, and G. Vosselman, “Segmentation of point clouds using smoothness constraints,” in _ISPRS 2006 : Proceedings of the ISPRS commission V symposium Vol. 35, part 6 : image engineering and vision metrology, Dresden, Germany 25-27 September 2006_ , H. Maas and D. Schneider, Eds., vol. 35. International Society for Photogrammetry and Remote Sensing (ISPRS), 2006.
  * [33] R. B. Rusu, N. Blodow, and M. Beetz, “Fast point feature histograms (fpfh) for 3d registration,” in _2009 IEEE International Conference on Robotics and Automation_ , 2009.
  * [34] N. X. Vinh, J. Epps, and J. Bailey, “Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance,” _J. Mach. Learn. Res._ , vol. 11, Dec. 2010.
  * [35] Q. Zhan, Y. Liang, and Y. Xiao, “Color-based segmentation of point clouds,” in _IAPRS XXXVIII Part 3/W8_ , 2009.
  * [36] J. Behley, M. Garbade, A. Milioto, J. Quenzel, S. Behnke, C. Stachniss, and J. Gall, “SemanticKITTI: A Dataset for Semantic Scene Understanding of LiDAR Sequences,” in _Proc. of the IEEE/CVF International Conf. on Computer Vision (ICCV)_ , 2019.