Plugging (30) in (31), we have with probability at least $1-\delta$,

$\displaystyle\int_{\Theta}f(\theta)\pi_{n}^{D}(\theta)\mathrm{d}\theta$
$\displaystyle\leq\inf_{\rho\in\mathcal{P}(\Theta)}\left\\{\int_{\Theta}f(\theta)\rho(\theta)\mathrm{d}\theta+\frac{1}{n}\text{KL}(\rho\|\pi)\right\\}+\frac{8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta\sqrt{n}}.$
(32)

Plugging the trivial bound $f(\theta)\leq
f(\theta_{*})+|f(\theta)-f(\theta_{*})|$ into (32), we have

$\displaystyle\eqref{eq:pc_eq4}\leq
f(\theta_{*})+\inf_{\rho\in\mathcal{P}(\Theta)}\bigg{\\{}\int_{\Theta}\left|f(\theta)-f(\theta_{*})\right|\rho(\theta)\mathrm{d}\theta+\frac{1}{n}\text{KL}(\rho\|\pi)\bigg{\\}}+\frac{8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta\sqrt{n}}.$

Notice that the infimum term can be upper bounded by any choice of
$\rho\in\mathcal{P}(\Theta)$. Letting
$\Pi(B_{n}):=\int_{B_{n}}\pi(\theta)\mathrm{d}\theta$, we take
$\rho(\theta)=\pi(\theta)/\Pi(B_{n})$ for $\theta\in B_{n}$ and
$\rho(\theta)=0$ for $\theta\not\in B_{n}$. Then Assumption 4 part (2) ensures
that
$\int_{B_{n}}|f(\theta)-f(\theta_{*})|\rho(\theta)\mathrm{d}\theta\leq\alpha_{1}/\sqrt{n}$
and that
$\mathrm{KL}(\rho\|\pi)=\int_{\Theta}\log\left(\rho(\theta)/\pi(\theta)\right)\rho(\theta)\mathrm{d}\theta=\int_{B_{n}}-\log(\Pi(B_{n}))\pi(\theta)\mathrm{d}\theta/\Pi(B_{n})=-\log\Pi(B_{n})\leq\alpha_{2}\sqrt{n}$.
Thus

$\displaystyle\int_{\Theta}f(\theta)\pi_{n}^{D}(\theta)\mathrm{d}\theta\leq
f(\theta_{*})+\left(\alpha_{1}+\alpha_{2}+\frac{8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta}\right)\frac{1}{\sqrt{n}},$
(33)

with probability at least $1-\delta$, as claimed. ∎

Now we turn to the proof of Theorem 1:

###### Proof of Theorem 1.

Since $\theta_{*}$ uniquely minimise $f$,

$\displaystyle f(\theta)-f(\theta_{*})\geq
0,\quad\forall\theta\in\Theta\qquad\Longrightarrow\qquad\int_{\Theta}f(\theta)\pi_{n}^{D}(\theta)\mathrm{d}\theta-f(\theta_{*})\geq
0.$

Thus, from Lemma 8,

$\displaystyle\mathbb{P}\left(\left|\int_{\Theta}f(\theta)\pi_{n}^{D}(\theta)\mathrm{d}\theta-f(\theta_{*})\right|\leq\left(\alpha_{1}+\alpha_{2}+\frac{8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta}\right)\frac{1}{\sqrt{n}}\right)\geq
1-\delta.$

Applying the simplifying upper bound

$\alpha_{1}+\alpha_{2}+\frac{8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta}\leq\frac{\alpha_{1}+\alpha_{2}+8\sup_{\theta\in\Theta}\sigma(\theta)}{\delta},$

taking complement of the probability and performing a change of variables, we
obtain the stated result. ∎

###  B.5 Proof of Theorem 2 (Bernstein–von Mises)

In this section we define the notation
$f_{n}(\theta):=\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})$
and $f(\theta):=\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})$.
Similarly, denote $H_{n}:=\nabla_{\theta}^{2}f_{n}(\theta_{n})$ and
$H_{*}=\nabla_{\theta}^{2}f(\theta_{*})$. Our aim is to verify the conditions
of Theorem 4 in Miller (2021). The following technical lemma lists and
establishes the conditions that are required:

######  Lemma 9.

Suppose that Assumption 2 ($r_{\max}=3$), Assumption 3, and part (1) of
Assumption 4 hold. Assume that $H$ is nonsingular. Then the following
statements almost surely hold:

  1. 1.

the prior density $\pi$ is continuous at $\theta_{*}$ and $\pi(\theta_{*})>0$,

  2. 2.

$\theta_{n}\to\theta_{*}$,

  3. 3.

the Taylor expansion
$f_{n}(\theta)=f_{n}(\theta_{n})+\frac{1}{2}(\theta-\theta_{n})\cdot
H_{n}(\theta-\theta_{n})+r_{n}(\theta-\theta_{n})$ holds on $\Theta$, where
the remainder $r_{n}$ satisfies $|r_{n}(\vartheta)|\leq
C\|\vartheta\|_{2}^{3}$ for all $\|\vartheta\|_{2}\leq\epsilon$, all
sufficiently large $n$ and some $C$ and $\epsilon>0$,

  4. 4.

$H_{n}\to H_{*}$, where $H_{n}$ is symmetric and $H_{*}$ is positive definite,

  5. 5.

$\liminf_{n\to\infty}(\inf_{\\{\theta\in\Theta:\|\theta-\theta_{n}\|_{2}\geq\epsilon\\}}f_{n}(\theta)-f_{n}(\theta_{n}))>0$
for any $\epsilon>0$.

###### Proof.

We sequentially prove each statement in the list.

Part (1): Directly assumed in Assumption 4 part (1).

Part (2): Assumption 2 ($r_{\max}=3$) and 3 are sufficient for Lemma 3 and
hence $\theta_{n}\overset{a.s.}{\to}\theta_{*}$.

Part (3): From Assumption 2 ($r_{\max}=3$), for all $h\in\mathcal{H}$ and
$x\in\mathcal{X}$ the map
$\theta\mapsto\mathcal{S}_{\mathbb{P}_{\theta}}[h](x)$ is three times
continuously differentiable, meaning that $f_{n}$ is three times continuously
differentiable on $\Theta$. Hence a second order Taylor expansion gives

$\displaystyle f_{n}(\theta)=f_{n}(\theta_{n})+\nabla
f_{n}(\theta_{n})(\theta-\theta_{n})+\frac{1}{2}(\theta-\theta_{n})\cdot
H_{n}(\theta-\theta_{n})+r_{n}(\theta-\theta_{n})$

where, for all sufficiently large $n$, $\nabla f_{n}(\theta_{n})=0$ was
assumed and the mean value form of the remainder term $r_{n}$ in the Taylor
expansion provides a bound

$\displaystyle|r_{n}(\theta-\theta_{n})|\leq\sup_{\theta\in\Theta}\|\nabla^{3}f_{n}(\theta)\|_{2}\|\theta-\theta_{n}\|_{2}^{3}.$

Finally,
$\limsup_{n\to\infty}\sup_{\theta\in\Theta}\|\nabla^{3}f_{n}(\theta)\|_{2}<\infty$
a.s. by the auxiliary Lemma 14 in Appendix C.

Part (4): $H_{n}$ is symmetric since the assumed regularity of $f_{n}$ allows
the mixed second order partial derivatives of $f_{n}$ to be interchanged. The
auxiliary Lemma 15 in in Appendix C establishes that
$H_{n}\overset{a.s.}{\to}H_{*}$ where $H_{*}$ is positive semi-definite. Thus,
since we assumed $H_{*}$ is nonsingular, it follows that $H_{*}$ is positive
definite.

Part (5): The inequality
$\liminf_{n\to\infty}(a_{n}+b_{n})\geq\liminf_{n\to\infty}a_{n}+\liminf_{n\to\infty}b_{n}$
holds for any sequences of $a_{n},b_{n}\in\mathbb{R}$. Combining the property
$\liminf_{n\to\infty}(-b_{n})=-\limsup_{n\to\infty}b_{n}$, we have that
$\liminf_{n\to\infty}(a_{n}-b_{n})\geq\liminf_{n\to\infty}a_{n}-\limsup_{n\to\infty}b_{n}$.
Applying this inequality,

$\displaystyle\liminf_{n\to\infty}\left(\inf_{\\{\theta\in\Theta:\|\theta-\theta_{n}\|_{2}\geq\epsilon\\}}f_{n}(\theta)-f_{n}(\theta_{n})\right)\geq\liminf_{n\to\infty}\inf_{\\{\theta\in\Theta:\|\theta-\theta_{n}\|_{2}\geq\epsilon\\}}f_{n}(\theta)-\limsup_{n\to\infty}f_{n}(\theta_{n})=:(*)$

Since $f_{n}(\cdot)\overset{a.s.}{\to}f(\cdot)$ uniformly on $\Theta$ by Lemma
2 and $\theta_{n}\overset{a.s.}{\to}\theta_{*}$ by Lemma 3,

$\displaystyle(*)$
$\displaystyle\overset{a.s.}{=}\inf_{\\{\theta\in\Theta:\|\theta-\theta_{*}\|_{2}\geq\epsilon\\}}f(\theta)-f(\theta_{*})>0$

where the last inequality follows from Assumption 3. ∎

Now we turn to the main proof:

###### Proof of Theorem 2.

Our aim is to verify the conditions of Theorem 4 in Miller (2021). Note that
this result in Miller (2021) views $\\{f_{n}\\}_{n=1}^{\infty}$ as a
deterministic sequence; we therefore aim to show that the conditions of
Theorem 4 in Miller (2021) are a.s. satisfied by our random sequence
$\\{f_{n}\\}_{n=1}^{\infty}$.

Recall that the generalised posterior has p.d.f.
$\pi_{n}^{D}(\theta)\propto\exp\left(-nf_{n}(\theta)\right)\pi(\theta)$
defined on $\Theta\subset\mathbb{R}^{p}$. This p.d.f. can be trivially
extended to a p.d.f. on $\mathbb{R}^{p}$ by defining $\pi(\theta)=0$ and
(e.g.) $f_{n}(\theta)=\inf_{\theta\in\Theta}f_{n}(\theta)+1$ for all
$\theta\in\mathbb{R}^{p}\setminus\Theta$. This brings us into the setting of
Miller (2021). The assumptions of Miller (2021, Theorem 4) are precisely the
list in the statement of Lemma 9, and the conclusion is that

$\displaystyle\int_{\mathbb{R}^{p}}\left|\hat{\pi}_{n}^{D}(\theta)-\frac{1}{\det(2\pi
H_{*}^{-1})^{1/2}}\exp\left(-\frac{1}{2}\theta\cdot
H_{*}\theta\right)\right|\mathrm{d}\theta\rightarrow 0.$ (34)

Thus, since from Lemma 9 the assumptions of Miller (2021, Theorem 4) are a.s.
satisfied, the conclusion in Equation 34 a.s. holds, as claimed. ∎

###  B.6 Proof of Robustness Results

####  B.6.1 Proof of Lemma 5

###### Proof.

First of all, (17) of Ghosh and Basu (2016) demonstrates that

$\displaystyle\operatorname{PIF}(y,\theta,\mathbb{P}_{n})$
$\displaystyle=\beta
n\pi_{n}^{L}(\theta)\left(-\operatorname{D}L(y,\theta,\mathbb{P}_{n})+\int_{\Theta}\operatorname{D}L(y,\theta^{\prime},\mathbb{P}_{n})\pi_{n}^{L}(\theta^{\prime})\mathrm{d}\theta^{\prime}\right).$

By Jensen’s inequality, we have an upper bounded

$\displaystyle\sup_{\theta\in\Theta}\sup_{y\in\mathcal{X}}|\operatorname{PIF}(y,\theta,\mathbb{P}_{n})|$
$\displaystyle\leq\beta
n\sup_{\theta\in\Theta}\pi_{n}^{L}(\theta)\left(\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|+\int_{\Theta}\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta^{\prime},\mathbb{P}_{n})\right|\pi_{n}^{L}(\theta^{\prime})\mathrm{d}\theta^{\prime}\right).$

Recall that $\pi_{n}^{L}(\theta)=\pi(\theta)\exp(-\beta
nL(\theta;\mathbb{P}_{n}))/Z$ where $0<Z<\infty$ is the normalising constant.
Thus we can obtain an upper bound
$\pi_{n}^{L}(\theta)\leq\pi(\theta)\exp(-\beta
n\inf_{\theta\in\Theta}L_{n}(\theta;\mathbb{P}_{n}))/Z=:C\pi(\theta)$ for some
constant $0<C<\infty$, since $L_{n}(\theta;\mathbb{P}_{n})$ is lower bounded
by assumption and $n$ is fixed. From this upper bound, we have

$\displaystyle\sup_{\theta\in\Theta}\sup_{y\in\mathcal{X}}|\operatorname{PIF}(y,\theta,\mathbb{P}_{n})|$
$\displaystyle\leq\beta
nC\sup_{\theta\in\Theta}\pi(\theta)\left(\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|+C\int_{\Theta}\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta^{\prime},\mathbb{P}_{n})\right|\pi(\theta^{\prime})\mathrm{d}\theta^{\prime}\right)$
$\displaystyle\leq\beta
nC\sup_{\theta\in\Theta}\left(\pi(\theta)\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|\right)+$
$\displaystyle\hskip 90.0pt\beta
nC^{2}\left(\sup_{\theta\in\Theta}\pi(\theta)\right)\int_{\Theta}\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta^{\prime},\mathbb{P}_{n})\right|\pi(\theta^{\prime})\mathrm{d}\theta^{\prime}.$

Since $\sup_{\theta\in\Theta}\pi(\theta)<\infty$ by assumption in the
statement of Lemma 5, it follows that

$\displaystyle\sup_{\theta\in\Theta}\left(\pi(\theta)\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|\right)<\infty\quad\text{
and
}\quad\int_{\Theta}\pi(\theta)\sup_{y\in\mathcal{X}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|\mathrm{d}\theta<\infty$

are sufficient conditions for
$\sup_{\theta\in\Theta}\sup_{y\in\mathcal{X}}|\operatorname{PIF}(y,\theta,\mathbb{P}_{n})|<\infty$,
as claimed. ∎

####  B.6.2 The Form of $\operatorname{D}L(y,\theta,\mathbb{P}_{n})$ for KSD

The following lemma clarifies the form of
$\operatorname{D}L(y,\theta,\mathbb{P}_{n})$ for KSD:

######  Lemma 10.

For
$L(\theta;\mathbb{P}_{n,\epsilon,y})=\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n,\epsilon,y})$,
we have

$\displaystyle\operatorname{D}L(y,\theta,\mathbb{P}_{n})=2\mathbb{E}_{X\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,y)\big{]}-2\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})\big{]}.$
(35)

###### Proof.

From the definition of the $\epsilon$-contamination model as a mixture model,
and using the symmetry of $K$, we have

$\displaystyle\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n,\epsilon,y})$
$\displaystyle={\mathbb{E}}_{X,X^{\prime}\sim\mathbb{P}_{n,\epsilon,y}}\left[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})\right]$
$\displaystyle=(1-\epsilon)^{2}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}_{n}}\left[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})\right]+2(1-\epsilon)\epsilon\mathbb{E}_{X\sim\mathbb{P}_{n}}\left[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,y)\right]$
$\displaystyle\hskip
40.0pt+\epsilon^{2}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y).$

Direct differentiation then yields

$\displaystyle\operatorname{D}L(y,\theta,\mathbb{P}_{n})$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n,\epsilon,y})\bigg{|}_{\epsilon=0}$
$\displaystyle=2\mathbb{E}_{X\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,y)\big{]}-2\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})\big{]},$

as claimed. ∎

####  B.6.3 Proof of Theorem 3

###### Proof.

From Lemma 5 with $\mathcal{X}=\mathbb{R}^{d}$, it is sufficient to show that

$\displaystyle\text{(i)}\
\sup_{\theta\in\Theta}\left(\pi(\theta)\sup_{y\in\mathbb{R}^{d}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|\right)<\infty\quad\text{
and }\quad\text{(ii)}\
\int_{\Theta}\sup_{y\in\mathbb{R}^{d}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|\pi(\theta)\mathrm{d}\theta<\infty.$

To establish (i) and (ii) we exploit the expression for
$\operatorname{D}L(y,\theta,\mathbb{P}_{n})$ in Lemma 10. This furnishes us
with the bound

$\displaystyle\big{|}\operatorname{D}L(y,\theta,\mathbb{P}_{n})\big{|}\leq
2\mathbb{E}_{X\sim\mathbb{P}_{n}}\big{[}\underbrace{|\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,y)|}_{=:(*_{1})}\big{]}+2\underbrace{\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})\big{]}}_{=:(*_{2})}.$
(36)

From Lemma 6,
$(*_{1})\leq\sqrt{\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)}\sqrt{\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)}$
and
$(*_{2})\leq\mathbb{E}_{X\sim\mathbb{P}_{n}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]$.
Plugging these bounds into (36) and using Jensen’s inequality gives

$\displaystyle\eqref{eq:bias-robust_GD_upb}$ $\displaystyle\leq
2\sqrt{\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)}\sqrt{\mathbb{E}_{X\sim\mathbb{P}_{n}}\left[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)\right]}+2\mathbb{E}_{X\sim\mathbb{P}_{n}}\left[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)\right].$
(37)

Now, observing that

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}_{n}}\big{[}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)\big{]}\leq\mathbb{E}_{X\sim\mathbb{P}_{n}}\Big{[}\sup_{y\in\mathbb{R}^{d}}(\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y))\Big{]}=\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)$
(38)

and taking a supremum over $y$ in (37), we obtain the bound

$\displaystyle\sup_{y\in\mathbb{R}^{d}}\left|\operatorname{D}L(y,\theta,\mathbb{P}_{n})\right|$
$\displaystyle\leq
4\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y).$
(39)

Therefore, from (39), it suffices to verify the conditions

$\displaystyle\text{(I)}\
\sup_{\theta\in\Theta}\left(\pi(\theta)\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)\right)<\infty\quad\text{
and }\quad\text{(II)}\
\int_{\Theta}\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)\pi(\theta)\mathrm{d}\theta<\infty,$

which imply the original conditions (i) and (ii). To this end, in the
remainder we (a) exploit the specific form of
$\mathcal{S}_{\mathbb{P}_{\theta}}$ to derive the an explicit upper bound on
$\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)$,
then (b) verify the conditions (I) and (II) based on this upper bound.

Part (a): By the reproducing property of $K$, the definition of the diffusion
Stein operator $\mathcal{S}_{\mathbb{P}_{\theta}}$, and the fact
$(a_{1}+a_{2})^{2}\leq 2(a_{1}^{2}+a_{2}^{2})$ for $a_{1},a_{2}\in\mathbb{R}$,
we have the bound

$\displaystyle\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)=\|\mathcal{S}_{\mathbb{P}_{\theta}}K(y,\cdot)\|_{\mathcal{H}}^{2}$
$\displaystyle=\left\|\nabla_{y}\log p_{\theta}(y)\cdot
K(y,\cdot)+\nabla_{y}\cdot K(y,\cdot)\right\|_{\mathcal{H}}^{2}$
$\displaystyle\leq 2\|\nabla_{y}\log p_{\theta}(y)\cdot
K(y,\cdot)\|_{\mathcal{H}}^{2}+2\left\|\nabla_{y}\cdot
K(y,\cdot)\right\|_{\mathcal{H}}^{2}.$

For the first term, the reproducing property of $K$ gives that

$\displaystyle\|\nabla_{y}\log p_{\theta}(y)\cdot
K(y,\cdot)\|_{\mathcal{H}}^{2}$ $\displaystyle=\nabla_{y}\log
p_{\theta}(y)\cdot K(y,y)\nabla_{y}\log p_{\theta}(y),$

while for the second term, and letting
$R(x,x^{\prime}):=\nabla_{x}\cdot(\nabla_{x^{\prime}}\cdot K(x,x^{\prime}))$,
the reproducing property gives that

$\displaystyle\left\|\nabla_{y}K(y,\cdot)\right\|_{\mathcal{H}}^{2}$
$\displaystyle=\big{\langle}\nabla_{y}\cdot K(y,\cdot),\nabla_{y}\cdot
K(y,\cdot)\big{\rangle}_{\mathcal{H}}=R(y,y).$

Thus, taking the supremum with respect to $y\in\mathbb{R}^{d}$ yields the
upper bound,

$\displaystyle\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)\leq
2\sup_{y\in\mathbb{R}^{d}}\big{(}\nabla_{y}\log p_{\theta}(y)\cdot
K(y,y)\nabla_{y}\log p_{\theta}(y)\big{)}+2\sup_{y\in\mathbb{R}^{d}}R(y,y).$

Since $K\in C_{b}^{1\times 1}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ by
assumption, it follows that $C_{MK}:=\sup_{y\in\mathbb{R}^{d}}R(y,y)<\infty$.
Thus we have arrived at

$\displaystyle\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)$
$\displaystyle\leq 2\gamma(\theta)+2C_{MK},$ (40)

where $\gamma(\theta)$ was defined in the statement of Theorem 3.

Part (b): Now we are in a position to verify conditions (I) and (II). For
condition (I), we use (40) to obtain

$\displaystyle\sup_{\theta\in\Theta}\Big{(}\pi(\theta)\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)\Big{)}\leq
2\sup_{\theta\in\Theta}\pi(\theta)\gamma(\theta)+2C_{MK}\sup_{\theta\in\Theta}\pi(\theta)$

which is finite by assumption. Similarly, for condition (II), we use (40) to
obtain

$\displaystyle\int_{\Theta}\sup_{y\in\mathbb{R}^{d}}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(y,y)\pi(\theta)\mathrm{d}\theta\leq
2\int_{\Theta}\pi(\theta)\gamma(\theta)\mathrm{d}\theta+2C_{MK}\int_{\Theta}\pi(\theta)\mathrm{d}\theta,$

which is also finite by assumption. This completes the proof. ∎

###  B.7 Verifying Assumptions 3, 2 and 4

In this appendix we demonstrate how Assumptions 3, 2 and 4 can be verified for
the exponential family model when the Langevin Stein operator is employed. For
simplicity, consider the case where the data dimension is $d=1$, the parameter
dimension is $p=1$, and the conjugate prior
$\pi(\theta)\propto\exp(-\theta^{2}/2)$ is used. From (9), a canonical
exponential family model with $\eta(\theta)=\theta$ and
$\mathcal{X}=\mathbb{R}$ is given by

$\displaystyle p_{\theta}(x)=\exp(\theta\cdot t(x)-a(\theta)+b(x))$

where $t:\mathbb{R}\to\mathbb{R}$, $a:\Theta\to\mathbb{R}$ and
$b:\mathbb{R}\to\mathbb{R}$. Accordingly, the log derivative is given by
$\nabla\log p_{\theta}(x)=\nabla t(x)\theta+\nabla b(x)$. Identical
calculations to Proposition 2 show that the KSD of the exponential family
model with the Langevin Stein operator takes a quadratic form

$\displaystyle\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})=C_{1,n}\theta^{2}+C_{2,n}\theta+C_{3,n}\quad\text{
and
}\quad\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})=C_{1}\theta^{2}+C_{2}\theta+C_{3}.$

where $C_{i,n}=(1/n^{2})\sum_{i,j=1}^{n}c_{i}(x_{i},x_{j})$ and
$C_{i}=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[c_{i}(X,X^{\prime})]$ and

$\displaystyle c_{1}(x,x^{\prime})$ $\displaystyle:=\nabla t(x)\cdot
K(x,x^{\prime})\nabla t(x^{\prime})$ $\displaystyle c_{2}(x,x^{\prime})$
$\displaystyle:=\nabla t(x)\cdot\big{(}\nabla_{x^{\prime}}\cdot
K(x,x^{\prime})\big{)}+\nabla t(x^{\prime})\cdot\big{(}\nabla_{x}\cdot
K(x,x^{\prime})\big{)}+2\nabla t(x)\cdot K(x,x^{\prime})\nabla b(x^{\prime})$
$\displaystyle c_{3}(x,x^{\prime})$ $\displaystyle:=b(x)\cdot
K(x,x^{\prime})b(x^{\prime})+\nabla_{x}\cdot(\nabla_{x^{\prime}}\cdot
K(x,x^{\prime}))$ $\displaystyle\hskip
160.0pt+b(x)\cdot(\nabla_{x^{\prime}}\cdot
K(x,x^{\prime}))+b(x^{\prime})\cdot(\nabla_{x}\cdot K(x,x^{\prime})).$

Note that $C_{1,n}>0$ and $C_{1}>0$ if a positive definite kernel $K$ is used.

Verifying Assumption 2 ($r_{\max}=3$): First, note that
$H_{*}=\nabla_{\theta}^{2}\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})|_{\theta=\theta_{*}}$
is non-singular since
$\nabla_{\theta}\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})=\nabla_{\theta}^{2}(C_{1}\theta^{2}+C_{2}\theta+C_{3})=2C_{1}>0$.
Now, as demonstrated in Section 4.1, when $\mathcal{S}_{\mathbb{P}_{\theta}}$
is the Langevin Stein operator, we have
$(\partial^{r}\mathcal{S}_{\mathbb{P}_{\theta}})[h](x)=(\partial^{r}\nabla_{x}\log
p_{\theta}(x))\cdot h(x)$ and
$h\mapsto(\partial^{r}\mathcal{S}_{\mathbb{P}_{\theta}})[h](x)$ is a
continuous linear functional on $\mathcal{H}$ for each fixed
$x\in\mathcal{X}$. In the exponential family case, the map
$\theta\mapsto\nabla_{x}\log p_{\theta}(x)$ is infinitely differentiable over
$\Theta$ since it is polynomial, leading to

$\displaystyle\partial^{1}\nabla_{x}\log p_{\theta}(x)=\nabla
t(x),\qquad\partial^{2}\nabla_{x}\log
p_{\theta}(x)=0,\qquad\partial^{3}\nabla_{x}\log p_{\theta}(x)=0.$

It is then clear that
$\mathbb{E}_{X\sim\mathbb{P}}[\sup_{\theta\in\Theta}((\partial^{r}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial^{r}\mathcal{S}_{\mathbb{P}_{\theta}})K(X,X))]<\infty$
for $r=2,3$. For $r=1$,

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}\left[\sup_{\theta\in\Theta}\big{(}(\partial^{1}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(X,X)\big{)}\right]=\mathbb{E}_{X\sim\mathbb{P}}\left[\nabla
t(X)\cdot K(X,X)\nabla t(X)\right].$ (41)

For the remaining term in Assumption 2, by essentially same calculations as
Proposition 2,

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}\left[\sup_{\theta\in\Theta}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)\big{)}\right]$
$\displaystyle=\mathbb{E}_{X\sim\mathbb{P}}\left[\sup_{\theta\in\Theta}\big{(}c_{1}(X,X)\theta^{2}+c_{2}(X,X)\theta+c_{3}(X,X)\big{)}\right]$
$\displaystyle\leq\mathbb{E}_{X\sim\mathbb{P}}\left[c_{1}(X,X)\sup_{\theta\in\Theta}\theta^{2}+c_{2}(X,X)\sup_{\theta\in\Theta}\theta+c_{3}(X,X)\right]$
(42)

Since $\Theta$ is a bounded set in $\mathbb{R}$, it is clear that
$\sup_{\theta\in\Theta}\theta<\infty$. The finiteness of (41) and (42) can
therefore be interpreted as finite moment conditions involving $t$, $b$, $K$
and $\mathbb{P}$.

Verifying Assumption 3: If both
$\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})$ and
$\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})$ are of quadratic
form with $C_{1,n}>0$ and $C_{1}>0$, the estimator $\theta_{n}$ exists and the
minimiser $\theta_{*}$ is unique over $\mathbb{R}$. It depends on
$C_{1},C_{2},C_{3}$ whether $\theta_{*}$ is contained in $\Theta$, but if we
are free to select $\Theta$ then we may select it such that
$\theta^{*}\in\Theta$. Since $C_{1}>0$, the _well-separated_ property of
$\theta_{*}$ is automatically satisfied; i.e.
$\operatorname{KSD}(\mathbb{P}_{\theta_{*}},\mathbb{P})<\inf_{\\{\theta\in\Theta:\|\theta-\theta_{*}\|_{2}\geq\epsilon\\}}\operatorname{KSD}(\mathbb{P}_{\theta},\mathbb{P})$
for all $\epsilon>0$.

Verifying Assumption 4: Part (1) is immediately satisfied since the prior
density $\pi(\theta)\propto\exp(-\theta^{2}/2)$ is continuous and positive on
$\Theta$. For part (2), we first have

$\displaystyle|\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})-\operatorname{KSD}^{2}(\mathbb{P}_{\theta_{*}}\|\mathbb{P})|=|C_{1}\theta^{2}+C_{2}\theta-
C_{1}\theta_{*}^{2}-C_{2}\theta_{*}|=C_{1}\left|(\theta+Z_{2})^{2}-Z_{1}\right|$

where $Z_{1}:=C_{2}^{2}/(4C_{1}^{2})+\theta_{*}^{2}+(C_{2}/C_{1})\theta_{*}$
and $Z_{2}:=C_{2}/(2C_{1})$ by completing the square. By the simple
calculation, the set
$B_{n}(\alpha_{1})=\\{\theta\in\Theta:|\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})-\operatorname{KSD}^{2}(\mathbb{P}_{\theta_{*}}\|\mathbb{P})|\leq\alpha_{1}/\sqrt{n}\\}$
is then given by

$\displaystyle
B_{n}(\alpha_{1})=\left\\{\theta\in\Theta:-\left(\frac{\alpha_{1}}{C_{1}\sqrt{n}}+Z_{1}\right)^{\frac{1}{2}}-Z_{2}\leq\theta\leq\left(\frac{\alpha_{1}}{C_{1}\sqrt{n}}+Z_{1}\right)^{\frac{1}{2}}-Z_{2}\right\\}$

While it is difficult to derive an explicit inequality between $\Pi(B_{n})$
and $\exp(-\alpha_{2}\sqrt{n})$, since it requires division into cases
according to the values of $C_{1},C_{2},C_{3}$, $\theta_{*}$, and the set
$\Theta$, the explicit form of $B_{n}$ renders it straightforward to
numerically determine which values for $\alpha_{1}>0$ and $\alpha_{2}>0$
ensure that $\Pi(B_{n})\geq\exp(-\alpha_{2}\sqrt{n})$ holds for all
$n\in\mathbb{N}$.

Quantities $S_{n}(x,\theta)$ and $J_{n}$: Here we provide the explicit form of
$S_{n}(x,\theta)$ and $J_{n}$ used to determine the value of $\beta$ for
exponential family model. From the definition,

$\displaystyle S_{n}(x,\theta)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x_{i})\big{)}=2\left(\frac{1}{n}\sum_{i=1}^{n}c_{1}(x,x_{i})\right)\theta+\left(\frac{1}{n}\sum_{i=1}^{n}c_{2}(x,x_{i})\right).$

Let $c_{1,n}(x):=(1/n)\sum_{i=1}^{n}c_{1}(x,x_{i})$ and
$c_{2,n}(x):=(1/n)\sum_{i=1}^{n}c_{2}(x,x_{i})$. From the definition,

$\displaystyle J_{n}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}S_{n}(x_{i},\theta_{n})S_{n}(x_{i},\theta_{n})^{\top}=\frac{1}{n}\sum_{i=1}^{n}\big{(}2c_{1,n}(x_{i})\theta_{n}+c_{2,n}(x_{i})\big{)}\big{(}2c_{1,n}(x_{i})\theta_{n}+c_{2,n}(x_{i})\big{)}^{\top}.$

Together with $H_{n}=C_{1,n}$, the default choice of $\beta$ is given by (14)
in Section 5.

##  Appendix C Auxiliary Theoretical Results

In Appendix B we exploited a number of auxiliary results, the details of which
are now provided. Recall that Standing Assumptions 1 and 2 continue to hold
throughout.

###  C.1 Derivative Bounds

Our auxiliary results mainly concern moments of derivative quantities, and the
aim of Section C.1 is to establish the main bounds that will be used. Recall
that $\partial^{1}$, $\partial^{2}$ and $\partial^{3}$ denote the partial
derivatives $(\partial/\partial\theta_{h})$,
$(\partial^{2}/\partial\theta_{h}\partial\theta_{k})$ and
$(\partial^{3}/\partial\theta_{h}\partial\theta_{k}\partial\theta_{l})$
respectively. For the proofs in Section C.1, we make the index explicit by re-
writing them as $\partial_{(h)}^{1}$, $\partial_{(h,k)}^{2}$ and
$\partial_{(h,k,l)}^{3}$. For $x\in\mathcal{X}$ and
$(h,k,l)\in\\{1,\ldots,p\\}^{3}$, we define

$\displaystyle m^{0}(x)$
$\displaystyle:=\sup_{\theta\in\Theta}\sqrt{\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x)},\qquad
m^{1}(x):=\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,x)},$
$\displaystyle m^{2}(x)$
$\displaystyle:=\sup_{\theta\in\Theta}\sqrt{\sum_{h,k=1}^{p}(\partial_{(h,k)}^{2}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h,k)}^{2}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,x)},$
$\displaystyle m^{3}(x)$
$\displaystyle:=\sup_{\theta\in\Theta}\sqrt{\sum_{h,k,l=1}^{p}(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,x)}.$

where we continue to use the convention that the first and second operator in
expressions such as
$(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,x^{\prime})$
are respectively applied to the first and second argument of $K$. Further
define

$\displaystyle M^{1}(x,x^{\prime})$
$\displaystyle:=m^{1}(x)m^{0}(x^{\prime})+m^{0}(x)m^{1}(x^{\prime}),$
$\displaystyle M^{2}(x,x^{\prime})$
$\displaystyle:=m^{2}(x)m^{0}(x^{\prime})+2m^{1}(x)m^{1}(x^{\prime})+m^{0}(x)m^{2}(x^{\prime}),$
$\displaystyle M^{3}(x,x^{\prime})$
$\displaystyle:=m^{3}(x)m^{0}(x^{\prime})+3m^{2}(x)m^{1}(x^{\prime})+3m^{1}(x)m^{2}(x^{\prime})+m^{0}(x)m^{3}(x^{\prime}).$

Based on these quantities, we now provide three technical results, Lemma 11,
Lemma 13 and Lemma 12.

######  Lemma 11.

Suppose Assumption 2 ($r_{max}=3$) holds. For each $r=1,2,3$, and for any
$x,x^{\prime}\in\mathcal{X}$,

$\displaystyle\sup_{\theta\in\Theta}\big{\|}\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\big{\|}_{2}$
$\displaystyle\leq M^{r}(x,x^{\prime}).$ (43)

If instead Assumption 2 ($r_{max}=1$) holds, then (43) holds for $r=1$.

###### Proof.

We first derive the upper bound for $r=1$ and then apply the same argument for
the remaining upper bound for $r=2$ and $r=3$. By the definition of
$\nabla_{\theta}$,

$\displaystyle\sup_{\theta\in\Theta}\Big{\|}\nabla_{\theta}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\Big{\|}_{2}$
$\displaystyle=\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}\left(\partial_{(h)}^{1}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\right)^{2}}.$
(44)

By Lemma 6 and Standing Assumption 2, we have
$\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot)\in\mathcal{H}$ for any
$x\in\mathcal{X}$ and

$\displaystyle(*_{1}):=\partial_{(h)}^{1}\Big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\Big{)}$
$\displaystyle=\partial_{(h)}^{1}\Big{(}\left\langle\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot),\mathcal{S}_{\mathbb{P}_{\theta}}K(x^{\prime},\cdot)\right\rangle_{\mathcal{H}}\Big{)}.$
(45)

From Assumption 2 ($r_{max}=1$), the operator
$(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})$ exists over $\Theta$
and satisfies the preconditions of Lemma 6. Hence, by setting
$\mathcal{S}_{\mathbb{Q}}=(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})$
in Lemma 6, we have that
$(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\in\mathcal{H}$
for each $x\in\mathcal{X}$. Let
$f_{\theta}(\cdot)=\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot)$ and
$g_{\theta}(\cdot)=\mathcal{S}_{\mathbb{P}_{\theta}}K(x^{\prime},\cdot)$. Then
the following product rule holds:

$\displaystyle\partial_{(h)}^{1}\langle
f_{\theta},g_{\theta}\rangle_{\mathcal{H}}$
$\displaystyle=\langle\partial_{(h)}^{1}f_{\theta},g_{\theta}\rangle_{\mathcal{H}}+\langle
f_{\theta},\partial_{(h)}^{1}g_{\theta}\rangle_{\mathcal{H}},$ (46)

which is verified from definition of differentiation as a limit and continuity
of the inner product. Note that
$\partial_{(h)}f_{\theta}(\cdot)=(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\in\mathcal{H}$
and
$\partial_{(h)}g_{\theta}(\cdot)=(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x^{\prime},\cdot)\in\mathcal{H}$.
Therefore by (46) and the Cauchy–Schwarz inequality,

$\displaystyle(*_{1})$
$\displaystyle=\left\langle\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot),\mathcal{S}_{\mathbb{P}_{\theta}}K(x^{\prime},\cdot)\right\rangle_{\mathcal{H}}+\left\langle\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot),\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}}K(x^{\prime},\cdot)\right\rangle_{\mathcal{H}}$
$\displaystyle\leq\underbrace{\left\|(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\right\|_{\mathcal{H}}}_{(*_{a})}\underbrace{\vphantom{\left|\partial_{(h)}^{1}\right\|_{\mathcal{H}}}\left\|\mathcal{S}_{\mathbb{P}_{\theta}}K(x^{\prime},\cdot)\right\|_{\mathcal{H}}}_{(*_{b})}+\underbrace{\vphantom{\left|\partial_{(h)}^{1}\right\|_{\mathcal{H}}}\left\|\mathcal{S}_{\mathbb{P}_{\theta}}K(x,\cdot)\right\|_{\mathcal{H}}}_{(*_{c})}\underbrace{\left\|(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x^{\prime},\cdot)\right\|_{\mathcal{H}}}_{(*_{d})}.$

For the original term (44), by the triangle inequality,

$\displaystyle\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}(*_{1})^{2}}$
$\displaystyle\leq\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}\Big{(}(*_{a})(*_{b})+(*_{c})(*_{d})\Big{)}^{2}}\leq\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}(*_{a})^{2}(*_{b})^{2}}+\sup_{\theta\in\Theta}\sqrt{\sum_{h=1}^{p}(*_{c})^{2}(*_{d})^{2}}.$

For the term $(*_{a})$, expanding the norm yields that

$\displaystyle(*_{a})^{2}$
$\displaystyle=\left\langle(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot),(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\right\rangle_{\mathcal{H}}=(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h)}^{1}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,x).$

A similar argument applied to $(*_{b})^{2}$, $(*_{c})^{2}$ and $(*_{d})^{2}$
leads to the overall bound

$\displaystyle\sup_{\theta\in\Theta}\big{\|}\nabla_{\theta}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\big{\|}_{2}$
$\displaystyle\leq
m^{1}(x)m^{0}(x^{\prime})+m^{0}(x)m^{1}(x^{\prime})=M^{1}(x,x^{\prime}).$

The upper bounds for $r=2$ and $r=3$ are obtained by an analogous argument.
Indeed, from the definition of $\nabla_{\theta}^{2}$ and
$\nabla_{\theta}^{3}$,

$\displaystyle\sup_{\theta\in\Theta}\left\|\nabla_{\theta}^{2}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\right\|_{2}$
$\displaystyle=\sup_{\theta\in\Theta}\sqrt{\sum_{h,k=1}^{p}\left(\partial_{(h,k)}^{2}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\right)^{2}}=:(*^{\prime\prime}),$
$\displaystyle\sup_{\theta\in\Theta}\left\|\nabla_{\theta}^{3}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\right\|_{2}$
$\displaystyle=\sup_{\theta\in\Theta}\sqrt{\sum_{h,k,l=1}^{p}\left(\partial_{(h,k,l)}^{3}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\right)^{2}}=:(*^{\prime\prime\prime}).$

From Assumption 2 ($r_{max}=3$), the operators
$(\partial_{(h,k)}^{2}\mathcal{S}_{\mathbb{P}_{\theta}})$ and
$(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})$ exist over
$\Theta$ and satisfy the preconditions of Lemma 6. Hence from Lemma 6,
$\partial_{(h,k)}^{2}f_{\theta}(\cdot)=(\partial_{(h,k)}^{2}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\in\mathcal{H}$
and
$\partial_{(h,k,l)}^{3}f_{\theta}(\cdot)=(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})K(x,\cdot)\in\mathcal{H}$
for any $x\in\mathcal{X}$, and in turn
$\partial_{(h,k)}^{2}g_{\theta}(\cdot)\in\mathcal{H}$ and
$\partial_{(h,k,l)}^{3}g_{\theta}(\cdot)\in\mathcal{H}$. Repeated application
of the product rule (46) gives that

$\displaystyle\partial_{(h,k)}^{2}\langle
f_{\theta},g_{\theta}\rangle_{\mathcal{H}}=$
$\displaystyle\langle\partial_{(h,k)}^{2}f_{\theta},g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(h)}^{1}f_{\theta},\partial_{(k)}^{1}g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(k)}^{1}f_{\theta},\partial_{(h)}^{1}g_{\theta}\rangle_{\mathcal{H}}+\langle
f_{\theta},\partial_{(h,k)}^{2}g_{\theta}\rangle_{\mathcal{H}},$
$\displaystyle\partial_{(h,k,l)}^{3}\langle
f_{\theta},g_{\theta}\rangle_{\mathcal{H}}=$
$\displaystyle\langle\partial_{(h,k,l)}^{3}f_{\theta},g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(h,k)}^{2}f_{\theta},\partial_{(l)}^{1}g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(h,l)}^{2}f_{\theta},\partial_{(k)}^{1}g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(k,l)}^{2}f_{\theta},\partial_{(h)}g_{\theta}\rangle_{\mathcal{H}}$
$\displaystyle+\langle\partial_{(h)}^{1}f_{\theta},\partial_{(k,l)}^{2}g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(k)}^{1}f_{\theta},\partial_{(h,l)}^{2}g_{\theta}\rangle_{\mathcal{H}}+\langle\partial_{(l)}^{1}f_{\theta},\partial_{(h,k)}^{2}g_{\theta}\rangle_{\mathcal{H}}+\langle
f_{\theta},\partial_{(h,k,l)}^{3}g_{\theta}\rangle_{\mathcal{H}}.$

Following the same argument as the preceding upper bound for $r=1$, the
triangle inequality and Cauchy–Schwarz imply that

$\displaystyle(*^{\prime\prime})$ $\displaystyle\leq
m^{2}(x)m^{0}(x^{\prime})+m^{1}(x)m^{1}(x^{\prime})+m^{1}(x)m^{1}(x^{\prime})+m^{0}(x)m^{2}(x^{\prime})$
$\displaystyle=m^{2}(x)m^{0}(x^{\prime})+2m^{1}(x)m^{1}(x^{\prime})+m^{0}(x)m^{2}(x^{\prime})=M^{2}(x,x^{\prime}),$
$\displaystyle(*^{\prime\prime\prime})$ $\displaystyle\leq
m^{3}(x)m^{0}(x^{\prime})+m^{2}(x)m^{1}(x^{\prime})+m^{2}(x)m^{1}(x^{\prime})+m^{2}(x)m^{1}(x^{\prime})$
$\displaystyle\qquad+m^{1}(x)m^{2}(x^{\prime})+m^{1}(x)m^{2}(x^{\prime})+m^{1}(x)m^{2}(x^{\prime})+m^{0}(x)m^{3}(x^{\prime})$
$\displaystyle=m^{3}(x)m^{0}(x^{\prime})+3m^{2}(x)m^{1}(x^{\prime})+3m^{1}(x)m^{2}(x^{\prime})+m^{0}(x)m^{3}(x^{\prime})=M^{3}(x,x^{\prime}),$

which are the claimed upper bounds for the cases $r=2$ and $r=3$. ∎

######  Lemma 12.

Suppose Assumption 2 ($r_{max}=3$) holds. For $r=0,1,2,3$,
$\mathbb{E}_{X\sim\mathbb{P}}[|m^{r}(X)|]<\infty$ and
$\mathbb{E}_{X\sim\mathbb{P}}[|m^{r}(X)|^{2}]<\infty$. For $r=1,2,3$,
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[|M^{r}(X,X^{\prime})|]<\infty$ and
$\mathbb{E}_{X\sim\mathbb{P}}[|M^{r}(X,X)|]<\infty$. If instead Assumption 2
($r_{max}=1$) holds, these results hold for $0\leq r\leq 1$.

###### Proof.

First, note that positivity of $m^{r}(\cdot)$ and $M^{r}(\cdot)$ implies that
the absolute value signs can be neglected. Moreover, from Jensen’s inequality
$(\mathbb{E}_{X\sim\mathbb{P}}[m^{r}(X)])^{2}\leq\mathbb{E}_{X\sim\mathbb{P}}[m^{r}(X)^{2}]$.
Thus it is sufficient to show that (a)
$\mathbb{E}_{X\sim\mathbb{P}}[m^{r}(X)^{2}]<\infty$, (b)
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]<\infty$ and (c)
$\mathbb{E}_{X\sim\mathbb{P}}[M^{r}(X,X)]<\infty$.

Part (a): The argument is analogous for each $r=0,1,2,3$ and we present it
with $r=3$. The bounded follows from Jensen’s inequality and the triangle
inequality:

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}\left[m^{3}(X)^{2}\right]$
$\displaystyle\leq\mathbb{E}_{X\sim\mathbb{P}}\left[\sup_{\theta\in\Theta}\sum_{h,k,l=1}^{p}(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})K(X,X)\right]$
$\displaystyle\leq\sum_{h,k,l=1}^{p}\mathbb{E}_{X\sim\mathbb{P}}\left[\sup_{\theta\in\Theta}\big{(}(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})(\partial_{(h,k,l)}^{3}\mathcal{S}_{\mathbb{P}_{\theta}})K(X,X)\big{)}\right]$

where the terms in the sum are finite by Assumption 2 ($r_{max}=3$).

Part (b): Since $X,X^{\prime}$ are independent in the expectation
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]$, it is clear
from the definition of $M^{r}$ that
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]$ exists if the
expectation of each term $m^{s}(X)$, $s\leq r$, exists. Thus by part (a),
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]<\infty$ for
$r=1,2,3$.

Part (c): From the definition of $M^{r}(x,x)$ for $r=1,2,3$,

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{1}(X,X)]$
$\displaystyle=2\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)m^{0}(X)],$
$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{2}(X,X)]$
$\displaystyle=2\mathbb{E}_{X\sim\mathbb{P}}[m^{2}(X)m^{0}(X)]+2\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)m^{1}(X)],$
$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{3}(X,X)]$
$\displaystyle=2\mathbb{E}_{X\sim\mathbb{P}}[m^{3}(X)m^{0}(X)]+6\mathbb{E}_{X\sim\mathbb{P}}[m^{2}(X)m^{1}(X)].$

Applying the Cauchy Schwartz inequality for each term

$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{1}(X,X)]$ $\displaystyle\leq
2\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)^{2}]}\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{0}(X)^{2}]},$
$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{2}(X,X)]$ $\displaystyle\leq
2\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{2}(X)^{2}]}\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{0}(X)^{2}]}+2\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)^{2}]}\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)^{2}]},$
$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[M^{3}(X,X)]$ $\displaystyle\leq
2\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{3}(X)^{2}]}\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{0}(X)^{2}]}+6\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{2}(X)^{2}]}\sqrt{\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)^{2}]}.$

Since each of the latter expectations is finite by part (a),
$\mathbb{E}_{X\sim\mathbb{P}}[M^{r}(X,X)]<\infty$ for $r=1,2,3$.

Inspection of the proof reveals that these results hold for $r=0,1$ if instead
Assumption 2 ($r_{max}=1$) holds. ∎

######  Lemma 13.

Suppose Assumption 2 ($r_{max}=3$) holds. Then, for $r=1,2,3$,

$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{r}(x_{i},x_{j})\overset{a.s}{\longrightarrow}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]<\infty.$
(47)

If instead Assumption 2 ($r_{max}=1$) holds, then (47) holds for $r=1$.

###### Proof.

The proof is based on the strong law of large numbers, the sufficient
conditions for which are provided by Lemma 12, which shows that
$\mathbb{E}_{X\sim\mathbb{P}}\left[|m^{r}(X)|\right]<\infty$ for $r=0,1,2,3$
under Assumption 2 ($r_{max}=3$). Then the strong law of large numbers
(Durrett, 2010, Theorem 2.5.10) yields that
$(1/n)\sum_{i=1}^{n}m^{r}(x_{i})\overset{a.s.}{\to}\mathbb{E}_{X\sim\mathbb{P}}\left[m^{r}(X)\right]=:(*_{r})$
for $r=0,1,2,3$. Then, from the definition of $M^{1}$,

$\displaystyle\lim_{n\to\infty}\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{1}(x_{i},x_{j})=\lim_{n\to\infty}\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\Big{(}m^{1}(x_{i})m^{0}(x_{j})+m^{0}(x_{i})m^{1}(x_{j})\Big{)}$
$\displaystyle\hskip
5.0pt=\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}m^{1}(x_{i})\times\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}m^{0}(x_{j})+\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}m^{0}(x_{i})\times\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}m^{1}(x_{j}).$

Since each limit in the right hand side converges a.s. to either $(*_{0})$ or
$(*_{1})$, so that

$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{1}(x_{i},x_{j})\overset{a.s.}{\longrightarrow}$
$\displaystyle\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)]\times\mathbb{E}_{X\sim\mathbb{P}}[m^{0}(X)]+\mathbb{E}_{X\sim\mathbb{P}}[m^{0}(X)]\times\mathbb{E}_{X\sim\mathbb{P}}[m^{1}(X)]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[m^{1}(X)m^{0}(X^{\prime})+m^{0}(X)m^{1}(X^{\prime})]=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{1}(X,X^{\prime})],$

where $X,X^{\prime}$ are independent. An analogous argument holds for
$M^{2}(x_{i},x_{j})$ and $M^{3}(x_{i},x_{j})$, giving that

$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{2}(x_{i},x_{j})$
$\displaystyle\overset{a.s.}{\longrightarrow}(*_{2})(*_{0})+2(*_{1})(*_{1})+(*_{0})(*_{2})=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{2}(X,X^{\prime})],$
$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{3}(x_{i},x_{j})$
$\displaystyle\overset{a.s.}{\longrightarrow}(*_{3})(*_{0})+3(*_{2})(*_{1})+3(*_{1})(*_{2})+(*_{0})(*_{3})=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{3}(X,X^{\prime})].$

Inspection of the proof reveals that (47) still holds for $r=1$ if Assumption
2 ($r_{max}=1$) holds instead. ∎

###  C.2 Proof of Auxiliary Results

Throughout this section we let
$f_{n}(\theta):=\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})$
and $f(\theta):=\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})$.
Similarly to $\nabla_{\theta}^{2}$, we let
$\nabla_{\theta}^{3}:=\nabla_{\theta}\otimes\nabla_{\theta}\otimes\nabla_{\theta}$
denote the tensor product $\otimes$ where each component is given by
$\partial_{h,k,l}^{3}$. For a matrix $a\in\mathbb{R}^{p\times p}$ and tensor
$b\in\mathbb{R}^{p\times p\times p}$, denote their Euclidean norms by
$\|a\|_{2}$ and $\|b\|_{2}$.

######  Lemma 14 (Derivatives a.s. Bounded).

Suppose Assumption 2 ($r_{max}=3$) holds. Then
$\limsup_{n\to\infty}\sup_{\theta\in\Theta}\left\|\nabla_{\theta}^{r}f_{n}(\theta)\right\|_{2}<\infty$
a.s. for $r=1,2,3$. If instead Assumption 2 ($r_{max}=1$) holds, then the
result holds for $r=1$.

###### Proof.

First of all, for finite $n$ we have

$\displaystyle\nabla_{\theta}^{r}f_{n}(\theta)=\nabla_{\theta}^{r}\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})\big{)}.$

From the triangle inequality and Lemma 11, we further have

$\displaystyle\sup_{\theta\in\Theta}\|\nabla_{\theta}^{r}f_{n}(\theta)\|_{2}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sup_{\theta\in\Theta}\left\|\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})\big{)}\right\|_{2}\leq\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}M^{r}(x_{i},x_{j}).$

It follows from Lemma 13 that
$(1/n^{2})\sum_{i=1}^{n}\sum_{j=1}^{n}M^{r}(x_{i},x_{j})\overset{a.s.}{\longrightarrow}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[M^{r}(X,X^{\prime})]<\infty$.
Therefore, a.s.
$\limsup_{n\to\infty}\sup_{\theta\in\Theta}\left\|\nabla_{\theta}^{r}f_{n}(\theta)\right\|_{2}<\infty$.
Inspection of the proof reveals that the argument still holds for $r=1$ if
Assumption 2 ($r_{max}=1$) holds instead. ∎

######  Lemma 15 (A.S. Convergence of Derivatives).

Suppose Assumption 2 ($r_{max}=3$) and 3 hold. Then we have
$\nabla_{\theta}^{r}f_{n}(\theta_{*})\overset{a.s.}{\to}\nabla_{\theta}^{r}f(\theta_{*})$
for $r=1,2,3$. Let $H_{n}:=\nabla_{\theta}^{2}f_{n}(\theta_{n})$ and
$H_{*}:=\nabla_{\theta}^{2}f(\theta_{*})$. We further have
$H_{n}\overset{a.s.}{\to}H_{*}$, where $H_{n}$ and $H_{*}$ are symmetric and
$H_{*}$ is semi positive definite.

###### Proof.

The proof is structured as follows: First we show (a)
$\nabla_{\theta}^{r}f_{n}(\theta_{*})\overset{a.s.}{\to}\nabla_{\theta}^{r}f(\theta_{*})$
for $r=1,2,3$. Then we show (b) $H_{n}\overset{a.s.}{\to}H_{*}$. Finally we
show (c) $H_{n}$ is symmetric and $H_{*}$ is semi-positive definite.

Part (a): The argument here is analogous to that used to prove Lemma 1, based
on the decomposition

$\displaystyle\nabla_{\theta}^{r}f_{n}(\theta)=\nabla_{\theta}^{r}\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})\big{)}.$

Let
$F(x,x^{\prime}):=\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}$
to see that

$\displaystyle\nabla_{\theta}^{r}f_{n}(\theta)$
$\displaystyle=\frac{1}{n}\underbrace{\frac{1}{n}\sum_{i=1}^{n}F(x_{i},x_{i})}_{(*_{1})}+\frac{n-1}{n}\underbrace{\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\neq
i}^{n}F(x_{i},x_{j})}_{(*_{2})}.$

It follows from the strong law of large number (Durrett, 2010, Theorem 2.5.10)
that $(*_{1})\overset{a.s.}{\to}\mathbb{E}_{X\sim\mathbb{P}}[F(X,X)]$ provided
$E_{X\sim\mathbb{P}}[\|F(X,X)\|_{2}]<\infty$. Similarly, it follows from the
strong law of large number for U-statistics (Hoeffding, 1961) that
$(*_{2})\overset{a.s.}{\to}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[F(X,X^{\prime})]$
provided $E_{X,X^{\prime}\sim\mathbb{P}}[\|F(X,X^{\prime})\|_{2}]<\infty$.
Both the required conditions holds by Lemma 12 and the fact that
$\|F(x,x^{\prime})\|_{2}\leq\sup_{\theta\in\Theta}\|\nabla_{\theta}^{r}(\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime}))\|_{2}\leq
M^{r}(x,x^{\prime})$ from Lemma 11. Thus

$\displaystyle\nabla_{\theta}^{r}f_{n}(\theta)$
$\displaystyle\overset{a.s.}{\longrightarrow}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[F(X,X^{\prime})]=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[\nabla_{\theta}^{r}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})\big{)}].$

Since
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[\|F(X,X^{\prime})\|_{2}]<\infty$, we
may apply the dominated convergence theorem to interchange expectation and
differentiation:

$\displaystyle\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[\nabla_{\theta}^{r}(\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime}))]=\nabla_{\theta}^{r}\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X^{\prime})]=\nabla_{\theta}^{r}f(\theta).$

Therefore, setting $\theta=\theta_{*}$, we conclude that
$\nabla_{\theta}^{r}f_{n}(\theta_{*})\overset{a.s.}{\to}\nabla_{\theta}^{r}f(\theta_{*})$.

Part (b): First of all, by the triangle inequality,

$\displaystyle\left\|\nabla_{\theta}^{2}f_{n}(\theta_{n})-\nabla_{\theta}^{2}f(\theta_{*})\right\|_{2}\leq\underbrace{\left\|\nabla_{\theta}^{2}f_{n}(\theta_{n})-\nabla_{\theta}^{2}f_{n}(\theta_{*})\right\|_{2}}_{(**_{1})}+\underbrace{\left\|\nabla_{\theta}^{2}f_{n}(\theta_{*})-\nabla_{\theta}^{2}f(\theta_{*})\right\|_{2}}_{(**_{2})}.$

By the mean value theorem applied to $(**_{1})$ and Lemma 14 (i.e.
$\lim_{n\to\infty}\sup_{\theta\in\Theta}\|\nabla_{\theta}^{3}f_{n}(\theta)\|_{2}<\infty$
a.s.), there a.s. exists a constant $0<C<\infty$ s.t., for all sufficiently
large $n$,

$\displaystyle(**_{1})=\left\|\nabla_{\theta}^{2}f_{n}(\theta_{n})-\nabla_{\theta}^{2}f_{n}(\theta_{*})\right\|_{2}$
$\displaystyle\leq\sup_{\theta\in\Theta}\|\nabla_{\theta}^{3}f_{n}(\theta)\|_{2}\|\theta_{n}-\theta_{*}\|_{2}\leq
C\|\theta_{n}-\theta_{*}\|_{2}.$

Then applying Lemma 3 (i.e.
$\|\theta_{n}-\theta_{*}\|_{2}\overset{a.s.}{\to}0$), we have
$(**_{1})\overset{a.s.}{\to}0$. Further the preceding part (a) implied that
$(**_{2})\overset{a.s.}{\to}0$. Therefore, we conclude that
$\nabla_{\theta}^{2}f_{n}(\theta_{n})\overset{a.s.}{\to}\nabla_{\theta}^{2}f(\theta_{*})$.

Part (c): Since $f_{n}$ is twice continuously differentiable over $\Theta$ by
assumption, commutation of two partial derivatives
$\partial_{(h)}\partial_{(k)}f_{n}(\theta)=\partial_{(k)}\partial_{(h)}f_{n}(\theta)$
holds over $\Theta$ by the Clairaut’s theorem. Therefore the $(h,k)$-th entry
and $(k,h)$-th entry of $H_{n}=\nabla_{\theta}^{2}f_{n}(\theta_{n})$ are
equal. An analogous argument applies to
$H_{*}=\nabla_{\theta}^{2}f(\theta_{*})$, so that both $H_{n}$ and $H_{*}$ are
symmetric. Furthermore, the Hessian $H_{*}$ is semi positive definite since
$\theta_{*}$ is the minimiser of $f$ from Assumption 3. ∎

######  Lemma 16 (Moment Condition for Asymptotic Normality).

Suppose that Assumption 2 ($r_{max}=3$) holds. Let
$F(x,x^{\prime}):=\nabla_{\theta}(\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime}))$
for any fixed $\theta\in\Theta$. Then we have
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}\left[\left\|F(X,X^{\prime})\right\|_{2}^{2}\right]<\infty$
and
$\mathbb{E}_{X\sim\mathbb{P}}\left[\left\|F(X,X)\right\|_{2}\right]<\infty$.

###### Proof.

First of all, it follows from Lemma 11 that for any
$x,x^{\prime}\in\mathcal{X}$,

$\displaystyle\|F(x,x^{\prime})\|_{2}\leq\sup_{\theta\in\Theta}\|\nabla_{\theta}\big{(}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})\big{)}\|_{2}\leq
M^{1}(x,x^{\prime}).$

Thus for the first moment we have
$\mathbb{E}_{X\sim\mathbb{P}}\left[\|F(X,X)\|_{2}\right]\leq\mathbb{E}_{X\sim\mathbb{P}}[M^{1}(X,X)]<\infty$
from Lemma 12. For the second moment,
$\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}\left[\left\|F(X,X^{\prime})\right\|_{2}^{2}\right]\leq\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}\left[M^{1}(X,X^{\prime})^{2}\right]=:(*)$.
By definition,

$\displaystyle(*)$
$\displaystyle=\mathbb{E}_{X,X^{\prime}\sim\mathbb{P}}\left[\big{(}m^{1}(X)m^{0}(X^{\prime})+m^{0}(X)m^{1}(X^{\prime})\big{)}^{2}\right]=4\mathbb{E}_{X\sim\mathbb{P}}\left[m^{1}(X)^{2}\right]\mathbb{E}_{X\sim\mathbb{P}}\left[m^{0}(X)^{2}\right].$

Each of these latter expectations is finite by Lemma 12, which completes the
proof. ∎

######  Theorem 4 (Concentration Inequality for KSD).

Let
$\sigma(\theta):={\mathbb{E}}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]$.
Then

$\displaystyle\mathbb{P}\left(\left|f_{n}(\theta)-f(\theta)\right|\geq\delta\right)\leq\frac{4\sigma(\theta)}{\delta\sqrt{n}},\quad\forall\theta\in\Theta,$

where the probability is with respect to realisations of the dataset
$\\{x_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathbb{P}$.

###### Proof.

Since $|a^{2}-b^{2}|=|(a+b)(a-b)|=(a+b)|a-b|$ for all $a,b\in[0,\infty)$, we
have the bound

$\displaystyle\underbrace{\left|\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})-\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})\right|}_{=:(*)}=\underbrace{\left(\operatorname{KSD}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})+\operatorname{KSD}(\mathbb{P}_{\theta}\|\mathbb{P})\right)}_{=:(*_{1})}\underbrace{\left|\operatorname{KSD}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})-\operatorname{KSD}(\mathbb{P}_{\theta}\|\mathbb{P})\right|}_{=:(*_{2})}.$

In what follows we use $\mathbb{E}$ to denote an expectation with respect to
the dataset $\\{x_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathbb{P}$. Applying Markov’s inequality followed by
Cauchy–Schwarz, we have

$\displaystyle\mathbb{P}((*)\geq\delta)\leq\frac{1}{\delta}\mathbb{E}[(*)]=\frac{1}{\delta}\mathbb{E}[(*_{1})(*_{2})]\leq\frac{1}{\delta}\sqrt{\mathbb{E}[(*_{1})^{2}]}\sqrt{\mathbb{E}[(*_{2})^{2}]}.$
(48)

To conclude the proof, we bound the two expectations one the right hand side.

Bounding $\mathbb{E}[(*_{1})^{2}]$: From the fact that $(a+b)^{2}\leq
2(a^{2}+b^{2})$ for $a,b\in\mathbb{R}$,

$\displaystyle\mathbb{E}[(*_{1})^{2}]\leq
2\mathbb{E}\left[\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})+\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})\right]=2\Big{(}\mathbb{E}\left[\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})\right]+\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})\Big{)}.$

The preconditions of Lemma 6 holds due to Standing Assumption 2. Thus from
Lemma 6 part (iii), together with Jensen’s inequality, we have the two bounds
$\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P}_{n})\leq(1/n)\sum_{i=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{i})$
and
$\operatorname{KSD}^{2}(\mathbb{P}_{\theta}\|\mathbb{P})\leq{\mathbb{E}}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]$.
Plugging these into the previous inequality, and exploiting independence of
$x_{i}$ and $x_{j}$ whenever $i\neq j$, we have

$\displaystyle\mathbb{E}[(*_{1})^{2}]$ $\displaystyle\leq
2\left(\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{i})\right]+{\mathbb{E}}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]\right)$
$\displaystyle=2\Big{(}\mathbb{E}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]+\mathbb{E}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]\Big{)}=4\sigma(\theta),$

where existence of $\sigma(\theta)$ for all $\theta\in\Theta$ is ensured by
Standing Assumption 2.

Bounding $\mathbb{E}[(*_{2})^{2}]$: From the fact
$|\sup_{x}|f(x)|-\sup_{y}|g(y)||\leq\sup_{x}|f(x)-g(x)|$ for functions $f$ and
$g$, the term $(*_{2})$ is upper bounded by

$\displaystyle(*_{2})$ $\displaystyle=\left|\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}[h](x_{i})\right|-\sup_{\|h\|_{\mathcal{H}}\leq
1}\Bigg{|}\mathbb{E}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}[h](X)]\Bigg{|}\right|$
$\displaystyle\leq\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}[h](x_{i})-\mathbb{E}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}[u](X)]\right|=\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}f(x_{i})-\mathbb{E}_{X\sim\mathbb{P}}[f(X)]\right|.$

where
$\mathcal{F}:=\\{\mathcal{S}_{\mathbb{P}_{\theta}}[h]\mid\|h\|_{\mathcal{H}}\leq
1\\}$. We can see from this expression that standard arguments in the context
of Rademacher complexity theory can be applied. Noting that $|\cdot|^{2}$ is a
convex function, Proposition 4.11 in Wainwright (2019) gives that

$\displaystyle\mathbb{E}\left[(*_{2})^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\left(\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}f(x_{i})-\mathbb{E}_{X\sim\mathbb{P}}[f(X)]\right|\right)^{2}\right]\leq\mathbb{E}\mathbb{E}_{\epsilon}\left[2^{2}\left(\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}f(x_{i})\right|\right)^{2}\right]$

where $\\{\epsilon_{i}\\}_{i=1}^{n}$ are independent random variables taking
values in $\\{-1,+1\\}$ with equiprobability $1/2$ and $\mathbb{E}_{\epsilon}$
is the expectation over $\\{\epsilon_{i}\\}_{i=1}^{n}$. From the essentially
same derivation as Proposition 1, the following equality holds:

$\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}f(x_{i})\right|=\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\mathcal{S}_{\mathbb{P}_{\theta}}[h](x_{i})\right|=\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\left\langle
h,\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},\cdot)\right\rangle_{\mathcal{H}}\right|\\\
=\left\|\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},\cdot)\right\|_{\mathcal{H}}=\sqrt{\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\epsilon_{i}\epsilon_{j}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})}.$

Plugging this equality into the upper bound of
$\mathbb{E}\left[(*_{2})^{2}\right]$, we have

$\displaystyle\mathbb{E}\left[(*_{2})^{2}\right]$ $\displaystyle\leq
4\mathbb{E}\mathbb{E}_{\epsilon}\left[\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\epsilon_{i}\epsilon_{j}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})\right]$
$\displaystyle=4\mathbb{E}\left[\frac{1}{n^{2}}\sum_{i=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}[K(X_{i},X_{i})]\right]=\frac{4}{n}{\mathbb{E}}_{X\sim\mathbb{P}}[\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(X,X)]=\frac{4\sigma(\theta)}{n}.$

Bounding $\mathbb{E}[(*)^{2}]$: Returning to (48), we have the overall bound

$\displaystyle\mathrm{P}((*)\geq\delta)\leq\frac{\sqrt{4\sigma(\theta)}\sqrt{4\sigma(\theta)}}{\delta\sqrt{n}}\leq\frac{4\sigma(\theta)}{\delta\sqrt{n}}$

as claimed. ∎

##  Appendix D Additional Empirical Results

This appendix contains additional empirical results referred to in the main
text. Section D.1 investigates the sensitivity of the generalised posterior to
the choice of parameters employed in the kernel $K$. Section D.2 investigates
the sampling distribution of $\beta$, controlling the scale of the generalised
posterior, when estimated using the approach proposed in Section 5.2. An
extended discussion of the choice of weighting function, $M$, and the
associated trade-off between statistical efficiency and robustness, is
contained in Section D.3. A comparison of KSD-Bayes with other generalised
Bayesian procedures developed for tractable likelihood is presented in Section
D.4. Finally, the use of KSD-Bayes in the context of discrete state spaces is
demonstrated in Section D.5.

###  D.1 Sensitivity to Kernel Parameters

Figure 7:  Sensitivity to kernel parameters: Kernels of the form (49), with
length-scale parameter $\sigma$ and exponent $\gamma$, are considered in the
context of the normal location model in Section 6.1. The settings
$\sigma\approx 1$, $\gamma=0.5$ (central panel) were used in the main text.
The true parameter value is $\theta=1$, while a proportion $\epsilon$ of the
data were contaminated by noise of the form $\mathcal{N}(y,1)$. Here $y=10$ is
fixed and $\epsilon\in\\{0,0.1,0.2\\}$ are considered.

The kernel $K$ that we recommend as a default in Section 5.1 has no degrees of
freedom to be specified (with the exception of the weighting function $M$,
whose choice is further explored in Section D.3). Nevertheless, it is
interesting to ask whether the generalised posterior is sensitive to our
recommended choice of kernel. To this end, we considered the family of kernels
of the form

$\displaystyle
K(x,x^{\prime})=\left(1+\sigma^{-2}\|x-x^{\prime}\|_{2}^{2}\right)^{-\gamma}\times
I_{d}$ (49)

where $\sigma>0$ and $\gamma\in(0,1)$. Our recommended kernel sets $\sigma$
equal to a regularised version of the sample standard deviation of the dataset
and $\gamma=1/2$. To investigate how the generalised KSD-Bayes posterior
depends on the choice of $\sigma$ and $\gamma$, we re-ran the normal location
model experiment from Section 6.1 using values $\sigma\in\\{0.5,1,2\\}$ and
$\gamma\in\\{0.1,0.5,0.9\\}$. To limit scope, we consider the performance of
the robust version of KSD-Bayes from Section 6.1, with weight function
$M(x)=(1+x^{2})^{-1/2}$, in the case where the contaminant is fixed to $y=10$
and the proportion of contamination is varied in $\epsilon\in\\{0,0.1,0.2\\}$.
Results in Figure 7 indicate that the generalised posterior is insensitive to
$\sigma$, with almost identical output for each value of $\sigma$ considered.
The results for $\gamma\in\\{0.5,0.9\\}$ were almost identical, but the
generalised posterior appeared to be less robust to contamination when
$\gamma=0.1$. These results support the default choices recommended in the
main text ($\sigma\approx 1$, $\gamma=0.5$) and provide reassurance that the
generalised posterior is not overly sensitive to how these values are
specified.

###  D.2 Sampling Distribution of $\beta$

An important component of the KSD-Bayes method is the use of a data-adaptive
$\beta$, as specified in Section 5.2. In this appendix the sampling
distribution of this data-adaptive $\beta$ is investigated. Of particular
interest are (1) the extent to which $\beta$ varies at small sample sizes, and
(2) how the behaviour of $\beta$ changes when the data-generating model is
mis-specified. To investigate, we considered multiple independent realisations
of the dataset in the context of the normal location model from Section 6.1,
collecting the corresponding estimates of $\beta$ together into box plots, so
that the sampling distribution of $\beta$ can be visualised. To limit scope,
we consider the performance of the standard version of KSD-Bayes from Section
6.1 (i.e. with weight function $M(x)=1$), in the case where the contaminant is
fixed to $y=10$ and the proportion of contamination is varied in
$\epsilon\in\\{0,0.1,0.2\\}$. The dataset sizes $n\in\\{10,50,100\\}$ were
considered. Results in Figure 8 show that, in the case $\epsilon=0$ where the
model is well-specified, the value $\beta=1$ is typically selected. This value
ensures that the scale of the KSD-Bayes posterior matches that of the standard
posterior in this example, so that the approach used to select $\beta$ can be
considered successful. In the mis-specified regimes
$\epsilon\in\\{0.1,0.2\\}$, with small $n$ the estimation of an appropriate
weight $\beta$ is expected to be difficult and indeed the default choice of
$\beta=1$ in (14) is automatically adopted. At larger values of $n$ it is
possible to reliably estimate a weight $\beta<1$ and this weight is seen to be
smaller on average when data are more contaminated. These results support our
recommended approach to selecting $\beta$ in (14).

Figure 8:  Sampling distribution of $\beta$: Box plots are used to summarise
the sampling distribution of $\beta$ in the context of the normal location
model in Section 6.1. The sample size $n$ and the contamination proportion
$\epsilon$ were each varied.

###  D.3 Efficiency/Robustness Trade-Off

There is a well-known trade-off between statistical efficiency and robustness
to model mis-specification, as exemplified by the data-agnostic statistician
who is robust by not learning from data. Minimum distance estimation, which
can be considered the frequentist analogue of generalised Bayesian inference,
can strike an attractive balance between these competing goals (see e.g.
Lindsay, 1994, Basu et al., 2019). In Section 4.3 it was demonstrated that
global bias-robustness can be achieved using KSD-Bayes through the inclusion
of an appropriate weighting function $M$ in the kernel, and in Section 6 it
was demonstrated that KSD-Bayes can learn from data whilst being bias-robust.
However, it remains to investigate the extent to which statistical efficiency
is lost in KSD-Bayes, compared to standard Bayesian inference, in the case
where the data-generating model is correctly specified. In this appendix we
return to the normal location model of Section 6.1 and explore the effect of
the choice of weighting function $M$ on the efficiency of the inferences that
are produced.

Recall from Theorem 3 that KSD-Bayes is globally bias-robust if there is a
function $\gamma:\Theta\rightarrow\mathbb{R}$ such that

$\displaystyle\sup_{y\in\mathbb{R}^{d}}\Big{(}\nabla_{y}\log
p_{\theta}(y)\cdot K(y,y)\nabla_{y}\log
p_{\theta}(y)\Big{)}\leq\gamma(\theta)$ (50)

where $\sup_{\theta\in\Theta}|\pi(\theta)\gamma(\theta)|<\infty$ and
$\int_{\Theta}\pi(\theta)\gamma(\theta)\mathrm{d}\theta<\infty$. For our
recommended kernel $K$ in (13), the expression on the left hand side of (50)
reduces to

$\sup_{y\in\mathbb{R}^{d}}\|M(y)^{\top}\nabla_{y}\log
p_{\theta}(y)\|_{2}^{2}.$

For the normal location model in Section 6.1 we have $\nabla_{y}\log
p_{\theta}(y)=\theta-y$ and thus, with our recommended kernel from Equation
13, we have

$\displaystyle\|M(y)^{\top}\nabla_{y}\log
p_{\theta}(y)\|_{2}^{2}=(y-\theta)^{2}M(y)^{2}.$ (51)

In order that (51) is bounded over $y\in\mathbb{R}$ we require $M(y)$ to decay
at the rate $\mathcal{O}(|y|^{-1})$ as $|y|\rightarrow\infty$. This decay is
achieved, for example, by functions of the form

$\displaystyle M(y)=\left(\frac{a^{2}}{a^{2}+(y-b)^{2}}\right)^{c/2}$ (52)

for any $a\neq 0$, $b\in\mathbb{R}$ and any $c\geq 1$, although of course
there are infinitely many other such functions that could be considered. The
particular value $c=1$, which we considered in Section 6.1 of the main text
and consider here in the sequel, represents the smallest value of $c$ for
which (51) is bounded over $y\in\mathbb{R}$. For this choice we have that (51)
is maximised by $y=\theta\pm\sqrt{a^{2}+(\theta-b)^{2}}$ and

$\displaystyle\sup_{y\in\mathbb{R}^{d}}(y-\theta)^{2}M(y)^{2}=\frac{[a^{2}+(\theta-b)^{2}]a^{2}}{a^{2}+[\theta-b\pm\sqrt{a^{2}+(\theta-b)^{2}}]^{2}}\leq
a^{2}+(\theta-b)^{2}=:\gamma(\theta).$

For this bound $\gamma(\theta)$, all conditions of Theorem 3 are satisfied.
The aim in what follows is to investigate how the performance of KSD-Bayes
depends on the specific choices of $a$ and $b$ and in (52).

To limit scope, we consider performance in the case where the contaminant is
fixed to $y=10$ and the proportion of contamination is varied in
$\epsilon\in\\{0,0.1,0.2\\}$. The dataset sizes was fixed at $n=100$ as per
the main text. Recall from Section 6.1 of the main text that the choices
$a=1$, $b=0$ lead to statistical efficiency comparable to that of standard
Bayesian inference. Results in Figure 9 show that $a=0.1$ led to almost total
robustness to contamination at the expense of inefficient estimation, with the
spread of the generalised posterior approximately twice as large as the case
where $a=1$. The setting $a=10$ causes the generalised posterior to
approximate the non-robust KSD-Bayes approach with $M\equiv 1$, as would be
expected from inspection of (52). The generalised posterior was somewhat
insensitive to $b$, though we note that the choice $b=-5$ conferred additional
robustness at the expense of efficiency, while the choice $b=5$ sacrificed
both robustness and efficiency, in both cases relative to $b=0$. These results
broadly support the choices of $a=1$ and $b=0$ for this inference problem, as
we considered in the main text.

Figure 9:  Efficiency/robustness trade-off: Weight functions of the form (52),
with length-scale parameter $a$ and location parameter $b$, are considered in
the context of the normal location model in Section 6.1. The settings $a=1$,
$b=0$ (central panel) were used in the main text. The true parameter value is
$\theta=1$, while a proportion $\epsilon$ of the data were contaminated by
noise of the form $\mathcal{N}(y,1)$. Here $y=10$ is fixed and
$\epsilon\in\\{0,0.1,0.2\\}$ are considered.

###  D.4 Comparison with Robust Generalised Bayesian Procedures

This paper presented a generalised Bayesian approach to inference for models
that involve an intractable likelihood. However, several generalised Bayesian
approaches exist for tractable likelihood and it is interesting to ask how the
performance of KSD-Bayes compares to these existing approaches in the case of
a tractable likelihood. To this end, we return to the normal location model of
Section 6.1, which has a tractable likelihood, and consider two distinct
generalised Bayesian procedures that have been developed in this context; the
power posterior approach of Holmes and Walker (2017) and the MMD-Bayes
approach of Cherief-Abdellatif and Alquier (2020). These approaches are
representative of two of the main classes of robust statistical methodology;
data-adaptive scaling parameters $\beta$ and minimum discrepancy methods. Both
approaches are briefly recalled:

##### Power Posteriors

Motivated by the coherence argument of Bissiri et al. (2016), the authors
Holmes and Walker (2017) consider a generalised posterior of the form, for
some $\beta>0$,

$\displaystyle\pi_{n}(\theta)\propto\pi(\theta)\exp\left\\{\beta\sum_{i=1}^{n}\log
p_{\theta}(x_{i})\right\\},$

which we call a power posterior (e.g. following Friel and Pettitt, 2008). To
select an appropriate value for $\beta$, with the intention to “allow for
Bayesian learning under model misspecification”, the authors first introduce
the function

$\displaystyle\Delta(x)=\int_{\Theta}\pi(\theta)\;\|\partial^{1}\log
p_{\theta}(x)\|_{2}^{2}\;\mathrm{d}\theta,$

where we recall that, in our notation,
$\partial^{1}=(\partial_{\theta_{1}},\dots,\partial_{\theta_{p}})$. Then the
authors set

$\displaystyle\beta=\left\\{\frac{\int_{\mathcal{X}}p_{\hat{\theta}_{n}}(x)\Delta(x)\mathrm{d}x}{\frac{1}{n}\sum_{i=1}^{n}\Delta(x_{i})}\right\\}^{\frac{1}{2}},$
(53)

where $\hat{\theta}_{n}$ is a maximiser of the likelihood. The motivation for
(53) is quite involved, so we refer the reader to Holmes and Walker (2017) for
further background. The authors prove that $\beta\rightarrow 1$ in probability
when the model is well-specified (Holmes and Walker, 2017, Lemma 2.1), and
present empirical evidence of robustness when the model is mis-specified.

For the normal location model of Section 6.1 we can compute $\partial^{1}\log
p_{\theta}(x)=x-\theta$, $\Delta(x)=1+x^{2}$,
$\hat{\theta}_{n}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$, and
$\int_{\mathcal{X}}p_{\hat{\theta}_{n}}(x)\Delta(x)\mathrm{d}x=2+(\hat{\theta}_{n})^{2}$,
leading to the recommended weight

$\displaystyle\beta=\left\\{\frac{2+\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}\right)^{2}}{1+\frac{1}{n}\sum_{i=1}^{n}x_{i}^{2}}\right\\}^{\frac{1}{2}}$

and an associated generalised posterior that is again Gaussian with mean
$(\frac{\beta n}{1+\beta n})(\frac{1}{n}\sum_{i=1}^{n}x_{i})$ and variance
$\frac{1}{1+\beta n}$.

##### MMD-Bayes

An analogue of KSD-Bayes for tractable likelihood is provided by the MMD-Bayes
approach of Cherief-Abdellatif and Alquier (2020), where a maximum mean
discrepancy (MMD) is employed in place of KSD. In identical notation to that
used in Section 3.1, the MMD-Bayes generalised posterior is defined, for some
$\beta>0$, as

$\displaystyle\pi_{n}^{D}(\theta)\propto\pi(\theta)\exp\left\\{-\beta
n\operatorname{MMD}^{2}(\mathbb{P}_{\theta},\mathbb{P}_{n})\right\\}$ (54)

where, for a given reproducing kernel Hilbert space $\mathcal{H}$ with
reproducing kernel $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, the
MMD between distributions $\mathbb{P}$ and $\mathbb{Q}$ on $\mathcal{X}$ is
defined as

$\displaystyle\operatorname{MMD}(\mathbb{P},\mathbb{Q})=\|\mu_{\mathbb{P}}-\mu_{\mathbb{Q}}\|_{\mathcal{H}},$

where the Bochner intergals
$\mu_{\mathbb{P}}(\cdot)=\int_{\mathcal{X}}k(\cdot,x)\mathrm{d}\mathbb{P}(x)$
and
$\mu_{\mathbb{Q}}(\cdot)=\int_{\mathcal{X}}k(\cdot,x)\mathrm{d}\mathbb{Q}(x)$
are the kernel mean embeddings of $\mathbb{P}$ and $\mathbb{Q}$ in
$\mathcal{H}$. The authors prove a generalisation bound for MMD-Bayes
(Cherief-Abdellatif and Alquier, 2020, Theorem 1), which they interpret as
showing “the MMD-Bayes posterior distribution is robust to misspecification”.
The authors do not recommend a default choice of $\beta$ in the main text888
Cherief-Abdellatif and Alquier (2020) absorbed the $n$ factor in (54) into
their definition of $\beta$, but for convenience of the reader we have
adjusted the presentation of MMD-Bayes to match that used for KSD-Bayes in the
main text., but in private correspondence they recommend $\beta=O(1)$, and we
use $\beta=1$ as a default. The kernel $k(x,y)=\exp(-\|x-y\|_{2}^{2}/d)$ was
used in our experiment, following Appendix F in Cherief-Abdellatif and Alquier
(2020).

For the normal location model of Section 6.1 we can compute the kernel mean
embeddings
$\mu_{\mathbb{P}_{\theta}}(x)=\sqrt{\frac{1}{3}}\exp(-\frac{1}{3}(x-\theta)^{2})$,
$\mu_{\mathbb{P}_{n}}(x)=\frac{1}{n}\exp(-(x-x_{i})^{2})$, obtaining an
overall expression for the MMD:

$\displaystyle\operatorname{MMD}(\mathbb{P}_{\theta},\mathbb{P}_{n})^{2}=\frac{1}{3}\exp\left(-\frac{\theta^{2}}{6}\right)-\frac{2}{n}\sum_{i=1}^{n}\sqrt{\frac{1}{3}}\exp\left(-\frac{(\theta-
x_{i})^{2}}{3}\right)+\frac{1}{n^{2}}\sum_{i,j=1}^{n}\exp\left(-(x_{i}-x_{j})^{2}\right)$

The un-normalised density associated with this generalised posterior can be
pointwise evaluated; we do this over a fine grid to approximate the
normalisation constant in the experiments that we report.

##### Results

The experiment of Section 6.1 was conducted using the power posterior and MMD-
Bayes methods just described, with results shown in Figure 10. Power
posteriors exhibited similar performance to (non-robust) KSD-Bayes (i.e. with
$M\equiv 1$; see Figure 2 in the main text), and was therefore less robust to
contamination compared with robust KSD-Bayes (i.e. with
$M(x)=(1+x^{2})^{-1/2}$). MMD-Bayes generalised posteriors provided similar
performance to robust KSD-Bayes in this experiment, albeit exhibiting greater
spread. The spread of the MMD-Bayes generalised posterior might be improved if
a data-adaptive learning rate $\beta$ is used, but such an approach was not
proposed in Cherief-Abdellatif and Alquier (2020).

Figure 10:  Comparison with robust generalised Bayesian procedures: Robust
KSD-Bayes (this paper), power posterior (Holmes and Walker, 2017) and MMD-
Bayes (Cherief-Abdellatif and Alquier, 2020) approaches are considered in the
context of the normal location model in Section 6.1. The true parameter value
is $\theta=1$, while a proportion $\epsilon$ of the data were contaminated by
noise of the form $\mathcal{N}(y,1)$. In the top row $y=10$ is fixed and
$\epsilon\in\\{0,0.1,0.2\\}$ are considered, while in the bottom row
$\epsilon=0.1$ is fixed and $y\in\\{1,10,20\\}$ are considered.

###  D.5 Application to Discrete Data

This section illustrates how KSD-Bayes may be applied to an intractable
discrete-space model; note that the theoretical results in Section 4.1 and
Section 4.2 cover both the discrete and continuous data context. For
demonstration purposes we consider a simple Ising model $\mathbb{P}_{\theta}$
on a vectorised $10\times 10$ lattice $\mathcal{X}=\\{-1,1\\}^{100}$, with a
temperature parameter $\theta\in(0,\infty)$, whose density is

$\displaystyle p_{\theta}(x)$
$\displaystyle\propto\exp\Big{(}\theta^{-1}\sum_{(i,j)\in
E}x_{(i)}x_{(j)}\Big{)},$ (55)

where $x=\\{x_{(1)},\dots,x_{(100)}\\}$ and $E$ is a index set containing all
pairs of adjacent states in the lattice $\mathcal{X}$. The classical Ising
model describes the statistical mechanics of molecular magnetic dipoles, with
$\theta$ controlling the intensity of interaction between each adjacent
magnetic dipole. To construct a Stein operator in this setting we follow Yang
et al. (2018), defining the _difference_ operators $\nabla^{+}$ and
$\nabla^{-}$ for a function $h:\\{-1,1\\}^{d}\to\mathbb{R}$ as

$\displaystyle\nabla^{+}h(x)=\left[\begin{array}[]{c}h(x^{+,1})-h(x)\\\
\vdots\\\ h(x^{+,d})-h(x)\\\ \end{array}\right]\quad\text{ and
}\quad\nabla^{-}h(x)=\left[\begin{array}[]{c}h(x^{-,1})-h(x)\\\ \vdots\\\
h(x^{-,d})-h(x)\\\ \end{array}\right]$

where $x^{(+,i)}$ and $x^{(-,i)}$ are vectors whose $i$-th coordinate is $1$
if $x_{(i)}=-1$ and $-1$ if $x_{(i)}=1$, with all other coordinates identical
to their values in $x$. The difference operators can be extend to act element-
wise on vector-valued functions $h:\\{0,1\\}^{d}\to\mathbb{R}^{d}$, so that
$\nabla^{+}h(x)$ and $\nabla^{-}h(x)$ are $d\times d$ matrices whose $i$-th
columns are given, respectively, by $\nabla^{+}h_{i}(x)$ and
$\nabla^{-}h_{i}(x)$. Further, for a vector-valued function
$h:\\{0,1\\}^{d}\to\mathbb{R}^{d}$ we let $\nabla^{-}\cdot
h(x)=\sum_{i=1}^{d}h_{i}(x^{(-,i)})-h_{i}(x)=\operatorname{Tr}(\nabla^{-}h(x))\in\mathbb{R}$.
The operator $\nabla^{-}\cdot$ will be applied to a matrix-valued kernel $K$;
in the same manner as the divergence operator in continuous domain,
$\nabla^{-}\cdot h(x)$ takes a value in $\mathbb{R}^{d}$ for a matrix-valued
function $h:\\{0,1\\}^{d}\to\mathbb{R}^{d\times d}$ where $[\nabla^{-}\cdot
h(x)]_{i}=\sum_{j=1}^{d}h_{ij}(x^{(-,j)})-h_{ij}(x)$. The Stein operator
$\mathcal{S}_{\mathbb{P}_{\theta}}$ we consider in this example is as follows:

$\displaystyle\mathcal{S}_{\mathbb{P}_{\theta}}h(x)=\frac{\nabla^{+}p_{\theta}(x)}{p_{\theta}(x)}\cdot
h(x)+\nabla^{-}\cdot h(x)$

For an empirical distribution $\mathbb{P}_{n}$ associated to a dataset
$\\{x_{i}\\}_{i=1}^{n}$, the corresponding $\operatorname{KSD}$ i.e.
$\operatorname{KSD}^{2}(\mathbb{P}_{\theta},\mathbb{P}_{n})=(1/n^{2})\sum_{i=1}^{n}\sum_{j=1}^{n}\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x_{i},x_{j})$
is based on

$\mathcal{S}_{\mathbb{P}_{\theta}}\mathcal{S}_{\mathbb{P}_{\theta}}K(x,x^{\prime})=\frac{\nabla^{+}p_{\theta}(x)}{p_{\theta}(x)}\cdot
K(x,x^{\prime})\frac{\nabla^{+}p_{\theta}(x^{\prime})}{p_{\theta}(x^{\prime})}+\operatorname{tr}\big{(}\nabla^{-}_{x}\cdot\nabla^{-}_{x^{\prime}}\cdot
K(x,x^{\prime})\big{)}\\\
+\frac{\nabla^{+}p_{\theta}(x)}{p_{\theta}(x)}\cdot\big{(}\nabla^{-}_{x^{\prime}}\cdot
K(x,x^{\prime})\big{)}+\big{(}\nabla^{-}_{x}\cdot
K(x,x^{\prime})\big{)}\cdot\frac{\nabla^{+}p_{\theta}(x^{\prime})}{p_{\theta}(x^{\prime})}.$

where $\nabla^{-}_{x^{\prime}}$ denotes an action of the operator
$\nabla^{-}\cdot$ with respect to the argument $x^{\prime}$ and likewise for
$\nabla^{-}_{x}$. See Yang et al. (2018) for further detail. The availability
of a discrete KSD enables the application of our KSD-Bayes methodology to the
Ising model.

Figure 11: Application to discrete data: Left: Samples from the Ising model
(55) at temperature $\theta=5$, defined on a $10\times 10$ lattice, where
white cells represent $+1$ and black cells represent $-1$. Right: The KSD-
Bayes and robust KSD-Bayes generalised posteriors obtained from an
uncontaminated ($\epsilon=0$) and contaminated ($\epsilon=0.1$) dataset.

As an empirical demonstration, we consider the same setting as Yang et al.
(2018); we approximately draw 1000 samples $\\{x_{i}\\}_{i=1}^{1000}$ from
$\mathbb{P}_{\theta}$ with $\theta=5$ using thinned MCMC (see Figure 11,
left). The prior $\pi$ was taken to be a half-normal distribution over
$\Theta=[0,\infty)$ with the scale hyper-parameter $3.0$. Our focus is on
robustness of the generalised posterior, and for the contamination model we
replaced a proportion $\epsilon$ of the data with the vector $(1,1,\cdots,1)$,
corresponding to the all-white lattice (a configuration more typically
observed at low values of the temperature parameter $\theta$). For KSD-Bayes,
the kernel in Yang et al. (2018) was used in combination with a weighting
function $M(x)$, i.e. our kernel is

$\displaystyle
K(x,x^{\prime})=M(x)\exp\left(-\frac{1}{2d}\sum_{i=1}^{d}|x_{(i)}-x_{(i)}^{\prime}|\right)M(x)^{\top}$

where $d=100$. For the weighting function $M(x)$, we examined two choices: (i)
$M(x)=I_{d}$ and (ii) $M(x)=\mathbbm{1}\\{|\sum_{i}x_{(i)}|\leq 90\\}\times
I_{d}$. The kernel in case (i) coincides with the one used in Yang et al.
(2018). The weighting function in case (ii) is designed to limit the influence
of data whose coordinates are almost all equal. The generalised posterior in
cases (i) and (ii) will be called, respectively, the KSD-Bayes posterior and
the robust KSD-Bayes posterior. The KSD-Bayes and robust KSD-Bayes posteriors
were approximated using Hamiltonian Monte Carlo. For simplicity, the weight
$\beta=1$ was fixed in this experiment. Results in Figure 11 (right) present
the generalised posteriors for a uncontaminated ($\epsilon=0.0$) and
contaminated ($\epsilon=0.1$) dataset. It can be observed that both the KSD-
Bayes and robust KSD-Bayes posteriors place their mass near the true parameter
$\theta=5$ when there is no contamination $\epsilon=0$. Furthermore, when
contamination is present, the robust KSD-Bayes posterior is not strongly
affected.

The computational challenge associated with discrete intractable likelihoods,
as exemplified by the Ising model, continues to attract attention (e.g. Kim et
al., 2021). Perhaps as a consequence, there has been little consideration of
robust estimation in this context. The nature of data contamination in
discrete spaces, and the extent to which this can be mitigated by careful
selection of the weighting function in KSD-Bayes, requires further careful
examination and will be addressed in a sequel. However, these preliminary
results are an encouraging proof-of-concept.