${\widetilde{e}}^{(s)}_{n}(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n})=\frac{1}{c_{n}^{3}}(\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}-c_{n}^{2}\epsilon_{n}-\frac{c_{n}^{3}\widetilde{q}_{0,sn}}{\sqrt{n_{s}}})$.
Then
$\Pi_{n}\left\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:{\widetilde{e}}^{(s)}_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)>C\tau_{n}\mid\mathcal{D}_{n}\right\\}\to
0$ in $\mathbb{P}_{0n}$-probability for a sufficiently large $M$ and a
constant $C>0$, if for every sufficiently large $j\in\mathbb{N}$ the following
conditions hold

  1. (I’)

For some $C>0$,
$\Pi_{n}\left\\{B_{n,0}(\mbox{\boldmath$\Theta$}_{0n},\tau_{n})\right\\}\geq
e^{-Cn\tau_{n}^{2}}$.

  2. (II’)

Define the set
$\widetilde{G}_{j,n}=\left\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}<{\widetilde{e}}^{(s)}_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)\leq
2j\tau_{n}\right\\}$. There exists tests $\varphi_{sn}$ such that for some
$K>0$,
$\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{0n}}\varphi_{sn}=0$
and
$\sup_{\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq\exp(-Knj^{2}\tau_{n}^{2})$.

Proof of the above lemma follows same line of arguments as Lemma 2.
Additionally condition (I’) is identical to (I) from Lemma 2. We show the
existence of a sequence of test functions satisfying condition (II’) in Lemma
5. This implies that

$\Pi_{n}\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:{\widetilde{e}}^{(s)}_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)>C\tau_{n}\mid\mathcal{D}_{n}\\}\to
0\text{ in $\mathbb{P}_{0n}$-probability}.$ (S.17)

Note that
${\widetilde{e}}^{(s)}_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)>C\tau_{n}\Leftrightarrow\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}>c_{n}^{2}\epsilon_{n}+\frac{c_{n}^{3}\widetilde{q}_{0,sn}}{\sqrt{n_{s}}}+Cc_{n}^{3}\tau_{n}$.
From respective definitions of $\epsilon_{n}$ and $\varepsilon_{n}$,
$c_{n}^{2}\epsilon_{n}\asymp\varepsilon_{n}$ and from (D1)
$\varepsilon_{n}\succ\max\left(\frac{c_{n}^{3}\widetilde{q}_{0,sn}}{\sqrt{n_{s}}},c_{n}^{3}\tau_{n}\right)$.
Thus (S.17) implies that for some large enough $M$,
$\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\to
0$ in $\mathbb{P}_{0n}$-probability. ∎

######  Lemma 5.

Recalling the definitions from Lemma 4, we let the set
$\widetilde{G}_{j,n}=\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}\leq{\widetilde{e}}^{(s)}_{n}(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n})\leq
2j\tau_{n}\\}$ for positive integer
$j=o\left[{\sqrt{n}}\left\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{-\frac{1}{2}}\right]$.
Based on observed data $\mathcal{D}_{n}$ consider the following hypothesis
testing problem
$H_{0}:\mbox{\boldmath$\Theta$}_{n}=\mbox{\boldmath$\Theta$}_{0n}$ versus
$H_{1}:\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}$. We use
$\varphi_{sn}=\mathbbm{1}\left(\norm{\frac{1}{n_{s}}\mbox{\boldmath$\Xi$}_{sn}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}>\epsilon_{n}\right)$
from Theorem 4 as the test function. Then for some absolute constant $K>0$

$\mathbb{E}_{H_{0}}(\varphi_{sn})\to
0;\qquad\sup_{\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq\exp(-Knj^{2}\tau_{n}^{2}).$

###### Proof.

Type-I error:

Please refer to corresponding section in the proof of Theorem 4 where it has
already been established that
$\lim_{n\to\infty}\mathbb{E}_{H_{0}}(\varphi_{sn})=0$.

Type-II error:

From equation (S.14),

$\norm{\bm{\Sigma}_{\eta,0sn}^{-1}\bm{\Sigma}_{\eta,sn}}_{2}\leq\norm{\bm{\Sigma}_{sn}}_{2}\leq\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}+\norm{\bm{\Sigma}_{sn}}_{2}<\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}+c_{n}.$
(S.18)

We obtain the last inequality in the previous display using (C3) and (C4).
Using (S.13)

$\norm{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}}{\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}(1+\norm{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}}_{2})}\geq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}}{c_{n}^{2}}.$
(S.19)

Combining (S.9), (S.18) and (S.19) we get

$\displaystyle 1-\varphi_{sn}$ $\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\left(\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}+c_{n}\right)\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}}{c_{n}^{2}}-\epsilon_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}-\epsilon_{n}c_{n}^{2}}{c_{n}^{2}\left(\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}+c_{n}\right)}\right\\}.$
(S.20)

Using
$j=o\left[{\sqrt{n}}\left\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{-\frac{1}{2}}\right]$
and similar arguments as in the proof of Lemma 10, it can be straightforwardly
seen that
$\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}-\epsilon_{n}c_{n}^{2}}{c_{n}^{2}\left(\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}+c_{n}\right)}>\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{1}{\sqrt{n_{s}}}$
for some absolute constant $\widetilde{C}>0$ if
$\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}$. Further applying
[Equation 5.26]vershynin_book, for all
$\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}$ we get

$\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq\Pr\left\\{\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{1}{\sqrt{n_{s}}}\right\\}\leq
e^{-Kj^{2}n_{s}\tau_{n}^{2}}.$

Using (C1) we conclude
$\sup_{\mbox{\boldmath$\Theta$}_{n}\in\widetilde{G}_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq
e^{-K^{\prime}j^{2}n\tau_{n}^{2}}$ for some $K^{\prime}>0$.

∎

###### Proof of Lemma 3.

Each HMC step comprises calculating the gradients of the $\log$-posterior
followed by evaluating the density at the proposed value. We separately
compute the numerical complexities of each.

##### Gradient computation:

As steps 2 and 2 are simply parallelized sums of objects with size at most
$qd$, it suffices to focus attention on steps 2-2. These are parallelized
across studies $s=1,\dots,S$, so we analyze the cost for each study as
follows.

Step 2:

Computing the Cholesky factor ${\mathbf{C}}_{s}^{\frac{1}{2}}$ is order
$\mathcal{O}(q^{3})$ since ${\mathbf{C}}_{s}$ has a low-rank plus diagonal
decomposition, and Cholesky decomposition of a $q\times q$ matrix is
$\mathcal{O}(q^{3})$ in the worst case. Next, the matrix product
$\widetilde{\bm{\Lambda}}_{s}=\bm{\Lambda}{\mathbf{C}}_{s}^{\frac{1}{2}}$ is
$\mathcal{O}(q^{2}d)$. Finally, owing to the low-rank and diagonal
decomposition structure, the inversion $\bm{\Sigma}_{s}^{-1}$ has complexity
$\mathcal{O}(q^{2}d)$. Therefore the overall complexity of this step is
$\mathcal{O}\\{\max(q^{2}d,q^{3})\\}=\mathcal{O}(q^{2}d)$ since in our setting
$q\ll d$.

Step 2:

Noting that
$\bm{\Sigma}_{s}^{-1}=(\widetilde{\bm{\Lambda}}_{s}\widetilde{\bm{\Lambda}}_{s}^{\rm
T}+\bm{\Delta})^{-1}=\bm{\Delta}^{-1}$
$-\bm{\Delta}^{-1}\widetilde{\bm{\Lambda}}_{s}({\mathbf{I}}_{q}+\widetilde{\bm{\Lambda}}_{s}^{\rm
T}\bm{\Delta}^{-1}\widetilde{\bm{\Lambda}}_{s})^{-1}\widetilde{\bm{\Lambda}}_{s}^{\rm
T}\bm{\Delta}^{-1}$ reveals that computing ${\mathbf{G}}_{s}$ is
$\mathcal{O}(qd^{2})$.

Steps 2 and 2:

Because the factor ${\mathbf{G}}_{s}$ is cached from the preceding step,
obtaining the derivative (following the flowchart in Figure 6) only requires
simple matrix multiplies, with complexity no more than $\mathcal{O}(qd^{2})$.

Hence, the dominant complexity of steps 2-2 is $\mathcal{O}(qd^{2})$, so that
the combined computational complexity of all leapfrog steps is
$\mathcal{O}(Lqd^{2})$.

#####  $\log$-posterior computation:

From equation (8) it can be seen that evaluating ${\cal L}$ at a proposal is
the expensive step involving high-dimensional matrix operations while the
prior densities are quite simple to evaluate. Therefore, we focus on the
complexity of calculating ${\cal L}$: owing to the structure of
$\bm{\Sigma}_{s}$, $\absolutevalue{\bm{\Sigma}_{s}}$ can be computed in
$\mathcal{O}(q^{2}d)$. Using the already cached
$\bm{\Sigma}_{s}^{-1}{\mathbf{W}}_{s}$ from step 2,
$\mathrm{trace}(\bm{\Sigma}_{s}^{-1}{\mathbf{W}}_{s})$ can be computed in
$\mathcal{O}(d)$. Both of the aforementioned operations are done in parallel
across $s$, so evaluating the $\log$-posterior has complexity no more than
$\mathcal{O}(q^{2}d)$.

Hence, the runtime complexity of each HMC step is $\mathcal{O}(Lqd^{2})$. ∎

###  Technical Lemmas and Complete Proofs

###### Proof of Lemma 2.

Lemma 2 closely resembles Theorem 8.22 from ghosal_book. In the discussion
following equation (8.22) in the book, the authors note that for the iid case,
simpler theorems are obtained by using an absolute lower bound on the prior
mass and by replacing the local entropy by the global entropy. In particular,
(I) implies Theorem 8.19(i) in the book. Similarly, (II) is the same condition
in Theorem 8.20, and so the result follows. ∎

######  Lemma 6.

For any two matrices ${\mathbf{A}}$ and ${\mathbf{B}}$,

  1. (i)

$s_{\min}\left({\mathbf{A}}\right)\norm{{\mathbf{B}}}_{F}\leq\norm{{\mathbf{A}}{\mathbf{B}}}_{F}\leq\norm{{\mathbf{A}}}_{2}\norm{{\mathbf{B}}}_{F}$.

  2. (ii)

$s_{\min}\left({\mathbf{A}}\right)\norm{{\mathbf{B}}}_{2}\leq\norm{{\mathbf{A}}{\mathbf{B}}}_{2}\leq\norm{{\mathbf{A}}}_{2}\norm{{\mathbf{B}}}_{2}$.

  3. (iii)

$s_{\min}\left({\mathbf{A}}\right)s_{\min}\left({\mathbf{B}}\right)\leq
s_{\min}\left({\mathbf{A}}{\mathbf{B}}\right)\leq\norm{{\mathbf{A}}}_{2}s_{\min}\left({\mathbf{B}}\right)$.

###### Proof.

We borrow Lemma 6 from pati2014 (Lemma 1.1. from the supplement). ∎

######  Lemma 7.

Let ${\mathbf{x}}=(x_{1},\dots,x_{n})^{\rm T}$ be a random vector with
$x_{i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(0,\sigma^{2})$
and ${\mathbf{x}}_{0}=(x_{01},\dots,x_{0n})^{\rm T}$ be a fixed vector. Then
for some absolute constant $C>0$,
$\Pr(\norm{{\mathbf{x}}-{\mathbf{x}}_{0}}\leq\tau)\geq
e^{-C\max\\{n\log\frac{\sqrt{\sigma}}{\tau},n\log
n,\frac{1}{\sigma^{2}}(\norm{{\mathbf{x}}_{0}}+\tau)^{2}\\}}$.

###### Proof.

Let us first define $v_{n}(r)$ to be the $n$-dimensional Euclidean ball of
radius $r$ centered at zero and $\absolutevalue{v_{n}(r)}$ denote its volume.
For the sake of brevity, denote $v_{n}=\absolutevalue{v_{n}(1)}$, so that
$\absolutevalue{v_{n}(r)}=r^{n}v_{n}$. Using [Lemma 5.3]castillo2012needles,
$v_{n}\asymp(2\pi e)^{n/2}n^{-(n+1)/2}$.

Note that
$\Pr(\norm{{\mathbf{x}}-{\mathbf{x}}_{0}}\leq\tau)\geq\absolutevalue{v_{n}(\tau)}\inf_{{\mathbf{z}}:\norm{{\mathbf{z}}-{\mathbf{x}}_{0}}\leq\tau}\mathrm{N}_{n}({\mathbf{z}};{\mathbf{0}},\sigma^{2}{\mathbf{I}}_{n})$
where $\mathrm{N}_{n}({\mathbf{z}};\mbox{\boldmath$\mu$},\bm{\Sigma})$ denotes
the density of a $n$ dimensional multivariate normal distribution with mean
$\mu$ and covariance matrix $\bm{\Sigma}$ at the point ${\mathbf{z}}$. Now the
greatest distance between the origin of the $n$-dimensional Euclidean space
and the $n$-dimensional sphere of radius $\tau$ centered at ${\mathbf{x}}_{0}$
is bounded above by $\norm{{\mathbf{x}}_{0}}+\tau$. Since the density of
$\mathrm{N}_{n}({\mathbf{0}},\sigma^{2}{\mathbf{I}}_{n})$ at a point
${\mathbf{z}}$ monotonically decreases with $\norm{{\mathbf{z}}}$, we have
$\inf_{{\mathbf{z}}:\norm{{\mathbf{z}}-{\mathbf{x}}_{0}}\leq\tau}\mathrm{N}_{n}({\mathbf{z}};{\mathbf{0}},\sigma^{2}{\mathbf{I}}_{n})\geq\frac{1}{(\sqrt{2\pi}\sigma)^{n}}e^{-\frac{1}{2\sigma^{2}}(\norm{{\mathbf{x}}_{0}}+\tau)^{2}}$.
Hence for absolute constants $C,C^{\prime},C^{\prime\prime}>0$

$\displaystyle\Pr(\norm{{\mathbf{x}}-{\mathbf{x}}_{0}}\leq\tau)\geq\absolutevalue{v_{n}(\tau)}\inf_{{\mathbf{z}}:\norm{{\mathbf{z}}-{\mathbf{x}}_{0}}\leq\tau}\mathrm{N}_{n}({\mathbf{z}};{\mathbf{0}},\sigma^{2}{\mathbf{I}}_{n})$
$\displaystyle\geq$ $\displaystyle C\tau^{n}(2\pi
e)^{n/2}n^{-{(n+1)}/2}\times\frac{1}{(\sqrt{2\pi}\sigma)^{n}}e^{-\frac{1}{2\sigma^{2}}(\norm{{\mathbf{x}}_{0}}+\tau)^{2}}\geq
e^{n\log\frac{\tau}{\sqrt{\sigma}}-C^{\prime}n(\log
n+C^{\prime\prime})-\frac{1}{2\sigma^{2}}(\norm{{\mathbf{x}}_{0}}+\tau)^{2}}$

which concludes the proof. ∎

######  Lemma 8.

$2\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)=\sum_{s=1}^{S}n_{s}\left\\{\mathrm{trace}(\bm{\Sigma}_{sn}^{-1}\bm{\Sigma}_{0sn}-{\mathbf{I}}_{d_{n}})-\log\absolutevalue{\bm{\Sigma}_{sn}^{-1}\bm{\Sigma}_{0sn}}\right\\}\leq\sum_{s=1}^{S}\frac{n_{s}\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}}{{\delta}_{\min}^{2}s_{\min}\left(\bm{\Sigma}_{sn}\right)}.$

###### Proof.

Recall that $\mathbb{P}_{n}$ and $\mathbb{P}_{0n}$ denote the joint
distributions of $\mathcal{D}_{n}$ under $\mbox{\boldmath$\Theta$}_{n}$ and
the true value $\mbox{\boldmath$\Theta$}_{0n}$, respectively. Let
$f_{\mathbb{P}}({\mathbf{Y}}_{s,i})$ denote the marginal distribution of
${\mathbf{Y}}_{s,i}$ under a measure $\mathbb{P}$. Since
$f_{\mathbb{P}_{n}}({\mathbf{Y}}_{s,i})\equiv\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Sigma}_{sn})$
and
$f_{\mathbb{P}_{0n}}({\mathbf{Y}}_{s,i})\equiv\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Sigma}_{0sn})$
for all $i$, and the ${\mathbf{Y}}_{s,i}$s are independent, we have

$\displaystyle\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)$
$\displaystyle=\sum_{s=1}^{S}\sum_{i=1}^{n_{s}}\mathbb{KL}\left(f_{\mathbb{P}_{0n}}({\mathbf{Y}}_{s,i})\parallel
f_{\mathbb{P}_{n}}({\mathbf{Y}}_{s,i})\right)$
$\displaystyle=\sum_{s=1}^{S}\frac{n_{s}}{2}\left\\{\mathrm{trace}(\bm{\Sigma}_{sn}^{-1}\bm{\Sigma}_{0sn}-{\mathbf{I}}_{d_{n}})-\log\absolutevalue{\bm{\Sigma}_{sn}^{-1}\bm{\Sigma}_{0sn}}\right\\}.$
(S.21)

Let
${\mathbf{H}}=\bm{\Sigma}_{sn}^{-\frac{1}{2}}\bm{\Sigma}_{0sn}\bm{\Sigma}_{sn}^{-\frac{1}{2}}$
and
$\mathbb{KL}_{s}=-\log\absolutevalue{\bm{\Sigma}_{sn}^{-1}\bm{\Sigma}_{0sn}}+\mathrm{trace}(\bm{\Sigma}_{0sn}\bm{\Sigma}_{sn}^{-1}-{\mathbf{I}}_{d_{n}})$.
Letting $\psi_{1},\dots,\psi_{d_{n}}$ be the eigenvalues of ${\mathbf{H}}$, we
note that
$\mathbb{KL}_{s}=\sum_{j=1}^{d_{n}}\left\\{(\psi_{j}-1)-\log\psi_{j}\right\\}$.
Observe that for any $x>-1$, $\log(1+x)\geq\frac{x}{1+x}$; additionally
$\psi_{j}>0$ for all $j=1,\dots,d_{n}$ since these are eigenvalues of the
positive definite matrix ${\mathbf{H}}$. Hence $\log\psi_{j}\geq
1-\frac{1}{\psi_{j}}$ implying that
$(\psi_{j}-1)-\log\psi_{j}\leq\frac{(\psi_{j}-1)^{2}}{\psi_{j}}$. Therefore,

$\textstyle{\mathbb{KL}_{s}\leq\sum_{j=1}^{d_{n}}\frac{(\psi_{j}-1)^{2}}{\psi_{j}}=\norm{({\mathbf{H}}-{\mathbf{I}}_{d_{n}})^{2}{\mathbf{H}}^{-1}}_{F}=\norm{(\bm{\Sigma}_{0sn}-\bm{\Sigma}_{sn})^{2}\bm{\Sigma}_{sn}^{-1}{\mathbf{H}}^{-1}\bm{\Sigma}_{sn}^{-1}}_{F}}.$
(S.22)

Now
$\bm{\Sigma}_{sn}{\mathbf{H}}\bm{\Sigma}_{sn}=\bm{\Sigma}_{sn}^{\frac{1}{2}}\bm{\Sigma}_{0sn}\bm{\Sigma}_{sn}^{\frac{1}{2}}$.
Hence, from (S.22) and using Lemma 6 we arrive at

$\mathbb{KL}_{s}\leq\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}\norm{(\bm{\Sigma}_{sn}^{\frac{1}{2}}\bm{\Sigma}_{0sn}\bm{\Sigma}_{sn}^{\frac{1}{2}})^{-1}}_{2}\leq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}}{s_{\min}\left(\bm{\Sigma}_{0sn}\right)s_{\min}\left(\bm{\Sigma}_{sn}\right)}\leq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}}{{\delta}_{\min}^{2}s_{\min}\left(\bm{\Sigma}_{sn}\right)}.$

Combining the above display with (S.21) we conclude the proof. ∎

######  Lemma 9.

Assume that
$\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}<\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}$,
$\max_{s}\norm{{\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn}}_{F}<\frac{C\tau_{n}}{c_{n}\sqrt{q_{0n}}}$,
$\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq
C\tau_{n}$ and $s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)>\nu$ where
$C=\nu{\delta}_{\min}^{2}/5$ and $0<\nu<{\delta}_{\min}^{2}$. Then
$2\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq
n\tau_{n}^{2}$.

###### Proof.

Using Lemma 8
$\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq\sum_{s=1}^{S}\frac{n_{s}\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}}{2{\delta}_{\min}^{2}s_{\min}\left(\bm{\Sigma}_{sn}\right)}$.
Now
$\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}\leq\norm{\bm{\Lambda}_{n}\bm{\Lambda}_{n}^{\rm
T}-\bm{\Lambda}_{0n}\bm{\Lambda}_{0n}^{\rm
T}}_{F}+\norm{\bm{\Lambda}_{sn}\bm{\Lambda}_{sn}^{\rm
T}-\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}}_{F}+\norm{\mbox{\boldmath$\Delta$}_{0n}-\mbox{\boldmath$\Delta$}_{n}}_{F}$
where $\bm{\Lambda}_{sn}:=\bm{\Lambda}_{n}{\mathbf{A}}_{sn}$ and
$\bm{\Lambda}_{0sn}:=\bm{\Lambda}_{0n}{\mathbf{A}}_{0sn}$. Recalling the
notations defined in the beginning of Section S.2 we have

$\displaystyle\norm{\bm{\Lambda}_{sn}\bm{\Lambda}_{sn}^{\rm
T}-\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}}_{F}=\norm{\bm{\Lambda}_{sn}\bm{\Lambda}_{sn}^{\rm
T}-\widetilde{\bm{\Lambda}}_{0sn}\widetilde{\bm{\Lambda}}_{0sn}^{\rm T}}_{F}$
$\displaystyle=$
$\displaystyle\norm{(\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn})\bm{\Lambda}_{sn}^{\rm
T}+\widetilde{\bm{\Lambda}}_{0sn}(\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn})^{\rm
T}}_{F}$ $\displaystyle\leq$
$\displaystyle\norm{\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn}}_{F}\times\left(\norm{\bm{\Lambda}_{sn}}_{2}+\norm{\widetilde{\bm{\Lambda}}_{0sn}}_{2}\right).$
(S.23)

We define $\bm{\Lambda}_{sn}=\bm{\Lambda}_{n}{\mathbf{A}}_{sn}$ and
$\widetilde{\bm{\Lambda}}_{0sn}=\widetilde{\bm{\Lambda}}_{0n}\widetilde{{\mathbf{A}}}_{0sn}$.
Then using Lemma 6

$\displaystyle\norm{\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn}}_{F}=\norm{\bm{\Lambda}_{n}{\mathbf{A}}_{sn}-\widetilde{\bm{\Lambda}}_{0n}\widetilde{{\mathbf{A}}}_{0sn}}_{F}=\norm{\bm{\Lambda}_{n}({\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn})+(\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n})\widetilde{{\mathbf{A}}}_{0sn}}_{F}$
$\displaystyle\leq$
$\displaystyle\norm{\bm{\Lambda}_{n}}_{2}\norm{{\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn}}_{F}+\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}\norm{\widetilde{{\mathbf{A}}}_{0sn}}_{2}.$
(S.24)

Note that
$\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}<\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}\Rightarrow\norm{\bm{\Lambda}_{n}}_{2}\leq\norm{\bm{\Lambda}_{n}}_{F}\leq\norm{\widetilde{\bm{\Lambda}}_{0n}}_{F}+\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}<\sqrt{c_{n}}$.
The last display holds since
$\ref{ass3}\Rightarrow\norm{\bm{\Lambda}_{0n}}_{2}=o(\sqrt{c_{n}/q_{0n}})\Rightarrow\norm{\bm{\Lambda}_{0n}}_{F}=o(\sqrt{c_{n}})$.
Combining the condition $\max_{1\leq s\leq
S}\norm{{\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn}}_{F}\leq\frac{C\tau_{n}}{c_{n}\sqrt{q_{0n}}}$
and (S.24) we get

$\max_{1\leq s\leq
S}\norm{\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn}}_{F}\leq\sqrt{c_{n}}\times\frac{C\tau_{n}}{c_{n}\sqrt{q_{0n}}}+\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}\times\sqrt{q_{0n}}<\frac{2C\tau_{n}}{\sqrt{c_{n}q_{0n}}}.$
(S.25)

Now
$\norm{\bm{\Lambda}_{sn}-\widetilde{\bm{\Lambda}}_{0sn}}_{F}<\frac{C\tau_{n}}{\sqrt{c_{n}q_{0n}}}\Rightarrow\norm{\bm{\Lambda}_{sn}}_{2}\leq\norm{\bm{\Lambda}_{sn}}_{F}<\norm{\widetilde{\bm{\Lambda}}_{0sn}}_{F}+\frac{C\tau_{n}}{\sqrt{c_{n}q_{0n}}}$.
Note that
$\norm{\widetilde{\bm{\Lambda}}_{0sn}}_{F}\leq\sqrt{q_{0sn}}\times\norm{\widetilde{\bm{\Lambda}}_{0sn}}_{2}\leq\sqrt{q_{0sn}}\times\norm{\bm{\Lambda}_{0n}}_{2}\norm{{\mathbf{A}}_{0sn}}_{2}<\sqrt{c_{n}q_{0n}}$
implying that
$\norm{\bm{\Lambda}_{sn}}_{2}+\norm{\widetilde{\bm{\Lambda}}_{0sn}}_{2}<2\sqrt{c_{n}q_{0n}}$.
Combining the last display with (S.23) and (S.25) we obtain
$\norm{\bm{\Lambda}_{sn}\bm{\Lambda}_{sn}^{\rm
T}-\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm T}}_{F}\leq 2C\tau_{n}$.
Similarly $\norm{\bm{\Lambda}_{n}\bm{\Lambda}_{n}^{\rm
T}-\bm{\Lambda}_{0n}\bm{\Lambda}_{0n}^{\rm
T}}_{F}\leq\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}(\norm{\bm{\Lambda}_{n}}_{F}+\norm{\widetilde{\bm{\Lambda}}_{0n}}_{F})\leq
2C\tau_{n}$. Using the additional condition
$\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq
C\tau_{n}$ and $s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)>\nu$, we get
$\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}\leq 5C\tau_{n}$. Furthermore
$s_{\min}\left(\bm{\Sigma}_{sn}\right)\geq
s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)$. Hence,
$2\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq\sum_{s=1}^{S}\frac{n_{s}\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{F}^{2}}{{\delta}_{\min}^{2}s_{\min}\left(\bm{\Sigma}_{sn}\right)}\leq
n\tau_{n}^{2}$. ∎

######  Lemma 10.

Recall the notations in Theorem 4. If $\mbox{\boldmath$\Theta$}_{n}\in
G_{j,n}$, then
$\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}-\epsilon_{n}c_{n}^{2}}{c_{n}^{2}\left(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}+c_{n}\right)\left(1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\right)}>\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{2}}{\sqrt{n_{s}}}$
for some absolute constant $\widetilde{C}>0$.

###### Proof.

For all $\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}$,
$1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\leq Cn\tau_{n}^{2}$ for some $C>0$ for all
$s=1,\dots,S$. For brevity of notations, we define $\omega=Cn\tau_{n}^{2}$,
$x=\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}$ and $a=\epsilon_{n}c_{n}^{2}$.
Hence observe that for $\mbox{\boldmath$\Theta$}_{n}\in G_{j,n}$

$\displaystyle\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}-\epsilon_{n}c_{n}^{2}}{c_{n}^{2}{\left(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}+c_{n}\right)\left(1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\right)}}>\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{1}{\sqrt{n_{s}}}$
(S.26) $\displaystyle\Leftrightarrow$
$\displaystyle\frac{x-a}{(x+c_{n})w}>\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{2}}{\sqrt{n_{s}}}\Leftarrow
x>\frac{a+\omega\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{3}}{\sqrt{n_{s}}}}{1-\omega\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{2}}{\sqrt{n_{s}}}}.$

The specifics on the postulated model from Section 3 imply
$\widetilde{q}_{0,sn}\leq q_{n}+q_{sn}=o(n\tau_{n}^{2})$. As
$j=o\left[\frac{\sqrt{n}}{c_{n}^{5}}\left\\{s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{-\frac{3}{2}}\right]$,
$\omega\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{2}}{\sqrt{n_{s}}}\precsim\frac{c_{n}^{5}}{\sqrt{n}}{\left\\{s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{\frac{3}{2}}}=o(1)$
which can be seen from (D1). Therefore,
$\frac{a+\omega\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{3}}{\sqrt{n_{s}}}}{1-\omega\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{c_{n}^{2}}{\sqrt{n_{s}}}}<a+j\times
Cc_{n}^{3}\sqrt{\frac{(n\tau_{n}^{2})^{3}}{n_{s}}}$ for some absolute constant
$C>0$ and large enough $n$. Hence, (S.26) holds if $x>a+j\times
Cc_{n}^{3}\sqrt{\frac{(n\tau_{n}^{2})^{3}}{n_{s}}}$ or equivalently if
$j\tau_{n}\leq
e_{n}(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n})$,
establishing the result. ∎

##  Appendix S.3 Hamiltonian Monte Carlo Algorithm

Let $\Pi(\mbox{\boldmath$\Theta$})$ be the target density where $\Pi(\cdot)$
is differentiable with respect to $\mbox{\boldmath$\Theta$}\in\mathbb{R}^{D}$.
In HMC, a dynamical system is considered in which auxiliary “momentum”
variables ${\mathbf{p}}\in\mathbb{R}^{D}$ are introduced and the uncertain
parameters $\Theta$ in the target distribution are treated as the variables
for the displacement. The total energy (Hamiltonian function) of the dynamical
system is defined by
$H(\mbox{\boldmath$\Theta$},{\mathbf{p}})=V(\mbox{\boldmath$\Theta$})+K({\mathbf{p}})/2$,
where its potential energy
$V(\mbox{\boldmath$\Theta$})=-\log\Pi(\mbox{\boldmath$\Theta$})$ and its
kinetic energy $K({\mathbf{p}})={\mathbf{p}}^{\rm
T}{\mathbf{M}}^{-1}{\mathbf{p}}$ depends only on ${\mathbf{p}}$ and some
chosen positive definite “mass” matrix ${\mathbf{M}}\in\mathbb{R}^{D\times
D}$. The Hamiltonian dynamics will preserve the distribution
$e^{-H(\mbox{\boldmath$\Theta$},{\mathbf{p}})}$ and the invariant distribution
will have $\Pi(\cdot)$ to be the marginal distribution of $\Theta$. Using
Hamilton’s equations, the evolution of $\Theta$, ${\mathbf{p}}$ through “time”
$t$ is given by

$\derivative{{\mathbf{p}}}{t}=-\partialderivative{H}{\mbox{\boldmath$\Theta$}}=-\nabla
V(\mbox{\boldmath$\Theta$}),\qquad\derivative{\mbox{\boldmath$\Theta$}}{t}=-\partialderivative{H}{{\mathbf{p}}}=-{\mathbf{M}}^{-1}{\mathbf{p}}.$
(S.27)

If we start with $\mbox{\boldmath$\Theta$}(0)$ and draw a sample
${\mathbf{p}}(0)$ from $\mathrm{N}({\mathbf{0}},{\mathbf{M}})$, the final
values $\mbox{\boldmath$\Theta$}(t)$, ${\mathbf{p}}(t)$ will provide an
independent sample $\mbox{\boldmath$\Theta$}(t)$ from $\Pi$. The leapfrog
algorithm DUANE_hmc is popularly used to approximately solve the differential
equations in (S.27). For time step $\delta t$, we have

$\displaystyle\mbox{\boldmath$\Theta$}(t+\delta
t)=\mbox{\boldmath$\Theta$}(t)+\delta
t{\mathbf{M}}^{-1}\left[{\mathbf{p}}(t)-\frac{\delta t}{2}\nabla
V\left\\{\mbox{\boldmath$\Theta$}(t)\right\\}\right],$ (S.28)
$\displaystyle{\mathbf{p}}(t+\delta t)={\mathbf{p}}(t)-\frac{\delta
t}{2}\left[\nabla V\left\\{\mbox{\boldmath$\Theta$}(t)\right\\}+\nabla
V\left\\{\mbox{\boldmath$\Theta$}(t+\delta t)\right\\}\right].$ (S.29)

The complete HMC sampler is summarized as follows for some ${\mathbf{M}}$,
$\delta t$, and $L$ leapfrog steps.

1

2Initialize $\Theta$ and simulate
${\mathbf{p}}(0)\sim\mathrm{N}({\mathbf{0}},{\mathbf{M}})$.

3 for  _$i=1,\dots,N$_ do

4 In iteration $i$, let the most recent sample be
$(\mbox{\boldmath$\Theta$}_{i-1},{\mathbf{p}}_{i-1})$, then do the following
to simulate a new sample $(\mbox{\boldmath$\Theta$}_{i},{\mathbf{p}}_{i})$:

  1. (i)

Randomly draw a new momentum vector ${\mathbf{p}}^{\prime}$ from
$\mathrm{N}({\mathbf{0}},{\mathbf{M}})$.

  2. 5(ii)

Initiate the leapfrog algorithm with
$\\{\mbox{\boldmath$\Theta$}(0),{\mathbf{p}}(0)\\}=(\mbox{\boldmath$\Theta$}_{i-1},{\mathbf{p}}^{\prime})$
and run the algorithm (S.28)-(S.29)

for $L$ time steps to obtain a new candidate sample
$(\mbox{\boldmath$\Theta$}^{\prime\prime},{\mathbf{p}}^{\prime\prime})=(\mbox{\boldmath$\Theta$}(t+L\delta
t),{\mathbf{p}}(t+L\delta t))$.

* 6(iii)

Set
$(\mbox{\boldmath$\Theta$}_{i},{\mathbf{p}}_{i})=(\mbox{\boldmath$\Theta$}^{\prime\prime},{\mathbf{p}}^{\prime\prime})$
with probability
$\min\left\\{1,e^{H\left(\mbox{\boldmath$\Theta$}^{\prime\prime},{\mathbf{p}}^{\prime\prime}\right)-H\left(\mbox{\boldmath$\Theta$}_{i-1},{\mathbf{p}}^{\prime}\right)}\right\\}$.

7 end for

Algorithm 3 Hamiltonian Monte Carlo

Thus we obtain the MCMC samples
$\mbox{\boldmath$\Theta$}_{1},\dots,\mbox{\boldmath$\Theta$}_{N}$.

###  Choice of Tuning Parameters

Practical performances of HMC samplers can be sensitive to the choice of
hyperparameters $\delta t$, $L$ and ${\mathbf{M}}$. To reduce sensitivity on
particular choices of the hyperparameters, in each MCMC iteration we randomly
sample $\delta t\sim\mathrm{Unif}(0,0.01)$, $L\sim\mathrm{Poisson}(5)$
truncated in the range $[1,10]$ and following usual practices set
${\mathbf{M}}$ to be the identity matrix.

##  Appendix S.4 Extended Simulation Studies

###  True Simulation Settings

##### Details on shared $\bm{\Lambda}$:

Scenario 1:

We fill a consecutive 25% of the elements in each column of $\bm{\Lambda}$
with independent samples from $\mathrm{Unif}(-2,2)$ and set the rest as zero.
We choose the starting position of the consecutive non-zero elements randomly.
To avoid having rows with all zero elements, we fill a randomly chosen 5
elements from each null row with iid samples from $\mathrm{Unif}(-2,2)$. A
typical example for a $200\times 20$ $\bm{\Lambda}$ is shown in the leftmost
panel of Figure S.1.

Scenario 2:

We separately fill the first and second consecutive $d/2$ values in each
column of $\bm{\Lambda}$. For each half we follow the same strategy as
scenario 1 and eliminate having entirely null rows using the aforementioned
scheme. A typical example for a $200\times 20$ $\bm{\Lambda}$ is shown in the
middle panel of Figure S.1. From the correlation structures shown under the
“True” panel of Figure 2 in the main paper, it can be seen that the above two
scenarios induce very different correlation structures in the marginal
distribution.

Scenario 3:

We randomly choose and fill 25% of elements in each column of $\bm{\Lambda}$
with independent samples from $\mathrm{Unif}(-2,2)$ and set the rest as zero.
We eliminate null rows using the same strategy of scenario 1. A typical
example for a $200\times 20$ $\bm{\Lambda}$ is shown in the rightmost panel of
Figure S.1.

Figure S.1: Simulation truths of $\bm{\Lambda}$ across different scenarios
under consideration. The gray color encodes exact 0 values.

##### Details on study-specific loadings:

Slight misspecification:

We set the study-specific loading
$\bm{\Phi}_{s}=\bm{\Lambda}{\mathbf{A}}_{s}+{\mathbf{E}}_{s}$ where we
generate ${\mathbf{A}}_{s}=((a_{s,j,h}))$ as
$a_{s,j,h}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(0,0.25^{2})$
and ${\mathbf{E}}_{s}=((e_{s,j,h}))$ as
$e_{s,j,h}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(0,0.10^{2})$.
Note that the SUFA model in (1) assumes that study-specific loading matrices
exactly admit the form $\bm{\Lambda}{\mathbf{A}}_{s}$, i.e., they are in
$\mathbb{C}\left(\bm{\Lambda}\right)$. Since we slightly deviate from that
assumption in this simulation setting, we refer to this as the “slight
misspecification” case.

Complete misspecification:

We generate $\bm{\Phi}_{s}$’s from the null-space of $\bm{\Lambda}$ for all
$s$ using the MASS package in R. We additionally take $\bm{\Phi}_{s}^{\rm
T}\bm{\Phi}_{s^{\prime}}={\mathbf{0}}$ for all $s\neq s^{\prime}$ to ensure
that the column spaces of study-specific loading matrices are non-overlapping.
We refer to this as the “complete misspecification” case since
$\mathbb{C}\left(\bm{\Lambda}\right)\cap\mathbb{C}\left(\bm{\Phi}_{s}\right)=\phi$
for all $s$.

Note that the shared and study-specific signals are characterized by the
magnitudes of the elements in the respective factor loading matrices. To
investigate the signal ratios in shared versus study-specific covariance
structures in the simulation truths, for each loading matrix we compute the
mean of squares for each column separately and compare the medians of these
mean squares across shared and study-specific loading matrices. In Figure S.2
we plot the histogram of these medians across independent replicates. The
figure indicates that the elements in the shared and study-specific loading
matrices are comparable in magnitude. This implies that despite having
similarities via the shared loading matrix $\bm{\Lambda}$, the marginal
covariance structures across studies have substantial differences.

Figure S.2: Investigation of the signal-to-noise ratios in shared versus
study-specific loadings in the simulation truths

##### Details on idiosyncratic variance $\bm{\Delta}$:

We set the diagonals of $\bm{\Delta}$ as $0.50$.

###  Recovering Shared Loading Matrix

In addition to the discussion in Section 5.1 in the main paper, here we
include an extended simulation study showcasing recovery of the shared factor
loading matrix $\bm{\Lambda}$. We apply the post-processing strategy discussed
in Section 2.1.2 on the MCMC samples to align them with respect to a common
orthogonal rotation and then use the mean of the post-processed MCMC samples
as a point estimate. The uncertainty quantification via the posterior samples
guide inference on whether entries of $\bm{\Lambda}=((\lambda_{j,h}))$ are $0$
or not; following common practice in Bayesian inferences on sparsity patterns
in matrices KsheeraSagar2021precision, we set $\lambda_{j,h}=0$ if its 95%
posterior credible interval includes 0.

Figure S.3: The true and estimated shared factor loading matrices: The rows
correspond to the three different scenarios, viz., FM1, FM2 and FM3 under
consideration. In each row, the true factor loading matrix is shown in the
leftmost plot; the middle and right plots are the point estimates in the
slight and complete misspecified cases, respectively.

In Figure S.3, we plot the shared factor loading matrices $\bm{\Lambda}$
recovered by SUFA. In each setting, as a representative of the 100 independent
replicates, we choose the one with the median Frobenius norm from the truth
with respect to the shared covariance matrix
$\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\rm T}+\bm{\Delta}$. We use the
approach elaborated in Section 2.1.2 of the main paper on the MCMC samples of
$\bm{\Lambda}$ to obtain a point estimate. For sake of brevity, we only show
results for $d=200$ in Figure S.3. The heatplots under the “True” panel of
Figure S.3 show the simulation truths of the three different correlation
structures; the “Slight” and “Complete” misspecification panels show the
recovered $\bm{\Lambda}$ by the SUFA model for the two model misspecification
types under consideration.

##### Comparison with alternate approaches:

Figure S.4: Coefficient of determination $R^{2}$ from regressing the columns
of the true $\bm{\Lambda}$ on the estimated $\widehat{\bm{\Lambda}}$ across
all simulations settings and alternate approaches under consideration.

Although the post-processing scheme aligns the MCMC samples of $\bm{\Lambda}$
with respect to a common orthogonal rotation, from Figure S.3 it is not
immediately obvious whether the estimate $\widehat{\bm{\Lambda}}$ is well-
aligned with the simulation truth. To assess this, we further regress the
columns of the true $\bm{\Lambda}$ on $\widehat{\bm{\Lambda}}$, and quantify
accuracy in terms of the coefficient of determination $R^{2}$.

To this end, we regress $\mbox{\boldmath$\lambda$}_{j\cdot}$ on
$\widehat{\bm{\Lambda}}$ for all $j=1,\dots,q$, where
$\mbox{\boldmath$\lambda$}_{j\cdot}$ denotes the $j^{th}$ column of the ground
truth. Considering $\widehat{\bm{\Lambda}}$ as the predictor matrix in this
fashion, we expect $R^{2}$ should be close to 1 if the true $\bm{\Lambda}$ is
well-learned up to orthogonal rotations, since the coefficient of
determination is invariant under (non-singular) matrix multiplication. We
report the median over the $R^{2}$ values obtained from the $q$ fitted
regressions, and show the boxplots of the median $R^{2}$ values from the
independent replicates in Figure S.4, across all simulation settings together
with estimates under peer methods. Note that the post-processing strategy in
Section 2.1.2 does not depend on any prior specification or knowledge of
$\bm{\Lambda}$; we use the same approach on samples from peer methods to
obtain point estimates of $\bm{\Lambda}$.

In Figure S.4, we see that even in the completely misspecified case,
performance of our proposed SUFA is comparable with the other approaches, and
evidently becomes better in higher dimensions. This is coherent with the
observations reported in Section 5.1 in the main paper. The performances of
BFR and PFA are not as favorable as their peers. Notably, neither BFR nor PFA
has separate components for inferences on study-specific correlation
structures between the variables. The better performance of SUFA compared to
B-MSFA is most likely due to nullification of the information switching issue
allowing improved learning of the shared vs study-specific structures.

###  Recovering Study-specific Loading Matrices

In this section we study how well the study-specific loading matrices are
recovered by the different methods. As BFR and PFA do not have study-specific
factor components we compare B-MSFA and SUFA only. From the posterior MCMC
samples of the study-specific loading matrices we obtain a point estimate in
the same manner described in Section S.4.2 and use the coefficient of
determination as the performance-metric.

Figure S.5: Coefficient of determination $R^{2}$ from regressing the columns
of the true study-specific loading matrices on their estimated counterparts
across all simulations settings for B-MSFA and SUFA.

B-MSFA recovers the study-specific loading matrices uniformly better than SUFA
in the complete misspecification case. This is expected as the simulation
truths of the study-specific loading matrices are from the null space of
$\bm{\Lambda}$, see Section S.4.1 for details. SUFA assumes that the study-
specific loadings are in the column space of the shared loading matrix
$\bm{\Lambda}$ and therefore it can not recover the study-specific loadings in
this setup.

On the contrary, there are overlaps between the column spaces of
$\bm{\Lambda}$ and the study-specific loading matrices in the slight
misspecification case. Performance of SUFA is more favorable in this case as
B-MSFA is susceptible to the information switching issue.

###  Scalability Analysis

In this section we study the scalability of the SUFA model with increasing
sample sizes. We generate data from the “complete misspecification” case in
“scenario 1”, see Section S.4.1 for details. We vary $d=50,200$ and $450$. For
each $d$ we set $n_{s}=\mathrm{Poisson}(m\times\frac{d}{S})$, for each
$s=1,\dots,S$ and repeat the experiments for $m\in\\{1,10,25\\}$. We generate
$q_{s}\sim\mathrm{Poisson}(\frac{q}{S})$ where $q=10$ for the $d=50$ and
$q=20$ otherwise. We fit SUFA using the HMC-within-Gibbs sampler in Algorithm
1 and report the boxplots of execution times (in minutes) across independent
replicates in Figure S.6. We used a system with $13^{th}$ Gen Intel(R)
Core(TM) i9-13900K CPU and 128GB RAM.

Figure S.6: Execution times (in minutes) of SUFA implementation with
increasing sample sizes and dimensions.

As expected, we see that the execution time is increasing with dimension.
However, it remains stable for increasing sample sizes for fixed dimensions.
This is because our proposed sampler uses the sufficient statistics calculated
before the MCMC iterations and becomes free of the sample size. This indicates
that our proposed approach is scalable to massive sample sizes.

##  Appendix S.5 Supplementary to Section 5.2

###  Heatmaps of the Gene Expressions

Figure S.7: Heatmaps of two microarray and a bulk RNAseq dataset on a common
set of genes. The X and Y axes represent genes and cells across studies,
respectively.

###  Pre-processing the Data

The normalized dataset has more than $21,000$ gene expressions from $628$ and
$146$ immune cells in the two microarray datasets, respectively, and more than
$49,000$ genes from $156$ cells in the bulk RNAseq dataset. We made a
$\log_{2}$ transformation of the data and filtered the top $5\%$ genes with
highest variances using the genefilter R package genefilter from each of the
datasets and considered the intersection of the filtered genes in our
analysis, resulting in a $d=474$ dimensional problem. Since different cell-
types exhibit very different gene expression profiles, we centered the gene-
expressions separately within each cell type.

###  Discussion on the Gene Network

In addition to corroborating qualitative results in the prior literature with
a principled statistical analysis, our study also reveals new insights via
integrating the data. Bcl2a1b and Bcl2a1d genes–two functional isoforms of the
B cell leukemia 2 family member bcl2 – show strong positive correlations. We
find strong positive associations between genes Bub1, Ccna2, Ccnb2, Top2a,
which corroborates findings from a study by tcellleukemia on adult human T
cell leukemia. Interestingly, we find this group of genes strongly negatively
correlated with those in the Gimap family, which are involved in lymphocyte
development and play important roles in immune system homeostasis. On the
other hand, the analysis reveals that the Gimap class features positive
within-group associations, supporting recent studies that suggest Gimap
proteins may interact with each other in roles such as moving cellular cargo
along the cytoskeletal network gimap. Up-regulation of the Ccr2 gene has been
found to be associated with cancer advancement, metastasis and relapse ccr2,
whereas Ccr5 inhibitors exhibit negative association with lung metastasis of
human breast cancer cell lines ccr5. Coherent with these studies, we find
strong positive correlation between Ccr2 and Ccr5. We also observe Ccr2 to be
highly positively correlated with Il18rap, supported by genome-wide
association studies identifying their susceptibility with coeliac diseases
il18rap. Additionally we see strong positive association between the
structurally related genes Il18rap and Il18r1 within the Il18 receptor complex
il18. Interestingly, strong negative correlation is observed between the Il18
family of genes and the gene Cd81, which is required for multiple normal
physiological functions cd81.

###  Interpretation of Factor Loading Matrices

Note that genes Il18r1, Il18rap, Klrk1, Ccl5, Slamf7, Tbx21, Ccr2, Ccr5, etc.
located in the bottom of Figure 5(a) have negative loadings corresponding to
the first latent factor in the shared space as well as in the bulkRNASeq space
but negative loadings for Microarray 2. To interpret the loading values we
primarily focus on the two bottom-most genes Il18rap and Il18r1 with
$\lambda_{g,h}$ and $\lambda_{g^{\prime},h}$ being their respective shared
loadings corresponding to factor $h$. From Figure 5(a) we see that
$\lambda_{g,h},\lambda_{g^{\prime},h}\approx 0$ for $h>1$. As covariances in
factor analysis models are entirely characterized by the dot-product between
the row-vectors of $\bm{\Lambda}$, for the sake of argument we assume that the
shared covariance between the two aforementioned genes is
$\mathrm{cov}_{g,g^{\prime}}\approx\lambda_{g,1}\lambda_{g^{\prime},1}$. For
characterizing $\mathrm{cov}_{g,g^{\prime}}$ the relative signs of
$\lambda_{g,1}$ and $\lambda_{g^{\prime},1}$ are meaningful, i.e.,
${\mathrm{sign}}(\lambda_{g,1})={\mathrm{sign}}(\lambda_{g^{\prime},1})$
implies $\mathrm{cov}_{g,g^{\prime}}>0$. Accordingly we observe strong
positive correlation between Il18rap and Il18r1 in the shared covariance
matrix.

Now let $\lambda^{\prime}_{g,1}$ and $\lambda^{\prime}_{g^{\prime},1}$ be the
loadings of Il18rap and Il18r1 respectively, corresponding to bulkRNASeq-
specific factor $1$. From Figure 5(a) we see that
${\mathrm{sign}}(\lambda^{\prime}_{g,1})=-{\mathrm{sign}}(\lambda_{g,1})$ and
${\mathrm{sign}}(\lambda^{\prime}_{g^{\prime},1})=-{\mathrm{sign}}(\lambda_{g^{\prime},1})$,
i.e., opposite signs in the shared and study-specific loadings. However, we
have
${\mathrm{sign}}(\lambda_{g,1}\lambda_{g^{\prime},1})={\mathrm{sign}}(\lambda^{\prime}_{g,1}\lambda^{\prime}_{g^{\prime},1})$
implying that the two genes exhibit positive correlation in the bulkRNASeq-
specific covariace part as well.

Similar phenomena can also be observed between genes Fosb, Nr4a1, Egr1, Dusp1,
Hspa1a, etc. located towards the top of Figure 5(a) that exhibit positive
loadings corresponding to latent factor 5 in the shared space but negative
loadings corresponding to the bulkRNASeq-specific latent factor 2. However,
the relative signs of the loadings within the shared space are the same
implying that they are positively correlated. The same is observed for the
bulkRNASeq-specific covariance as well. Therefore, we conclude that the signs
of the loadings must be interpreted together and not individually.

The magnitudes of the loadings indicate the absolute value of the respective
variances and covariances. For example the genes Klrk1 and Ccl5 (third and
fourth from the bottom in Figure 5(a)) have higher shared loadings compared to
most other genes and therefore they are expected to exhibit higher shared
covariance and variances.

natbib main


# Reliable Probability Intervals For Classification Using Inductive Venn
Predictors Based on Distance Learning

Dimitrios Boursinos and Xenofon Koutsoukos
Institute for Software Integrated Systems
Vanderbilt University
Nashville TN, USA
{dimitrios.boursinos<EMAIL_ADDRESS>
###### Abstract

Deep neural networks are frequently used by autonomous systems for their
ability to learn complex, non-linear data patterns and make accurate
predictions in dynamic environments. However, their use as black boxes
introduces risks as the confidence in each prediction is unknown. Different
frameworks have been proposed to compute accurate confidence measures along
with the predictions but at the same time introduce a number of limitations
like execution time overhead or inability to be used with high-dimensional
data. In this paper, we use the Inductive Venn Predictors framework for
computing probability intervals regarding the correctness of each prediction
in real-time. We propose taxonomies based on distance metric learning to
compute informative probability intervals in applications involving high-
dimensional inputs. Empirical evaluation on image classification and botnet
attacks detection in Internet-of-Things (IoT) applications demonstrates
improved accuracy and calibration. The proposed method is computationally
efficient, and therefore, can be used in real-time.

###### Index Terms:

deep neural networks, assurance monitoring, inductive venn predictors,
probability intervals

††publicationid: pubid: 978-1-6654-3156-9/21/$31.00 ©2021 IEEE

##  I Introduction

Modern Deep Neural Network (DNN) architectures have the capacity to be trained
using high-dimensional data and make accurate decisions in dynamic and
uncertain environments. This ability makes them a common choice for many
autonomous system applications. However, when DNNs are used as black boxes in
safety-critical systems, they may result in disastrous consequences if it is
not possible to reason about their predictions.

The training of a Learning Enabled Component (LEC) requires specification of
the task, performance measure for evaluating how well the task is performed,
and experience in the form of training and testing data. An LEC, such as a
DNN, during system operation exhibits some nonzero error rate and the true
error rate is unknown and can only be approximated during design time using
the available data. Confidence values, such as the softmax probabilities which
are used by most DNNs for classification, are usually greater than the actual
posterior probability that the prediction is correct. Important factors that
make modern DNNs overconfident are the depth, width, and techniques like
weight decay, and batch normalization [5].

Our objective is to complement the predictions made by DNNs with a computation
of confidence. The confidence can be expressed as probability intervals to
characterize the correctness of the DNN prediction. An efficient and robust
approach must ensure that the actual accuracy of a DNN is contained in the
computed intervals and the width of the intervals is small. We focus on
computationally efficient algorithms that can be used in real-time. The
proposed approach is based on the Inductive Venn Predictors (IVP) framework
[28]. IVP computes the probability intervals for an unknown input leveraging
knowledge it has acquired from previous predictions on labeled data. Most of
the IVP or Venn Predictors (VP) applications in the literature are evaluated
on low-dimensional data [29, 28, 13, 19, 14, 20].

The estimation of reliable predictive uncertainty has become an important part
of many modern machine learning components used in safety-critical
applications. Even though many of the proposed methods produce well-calibrated
models, their application in the real world is challenging. In [22, 12], new
training algorithms and loss functions are proposed to achieve well-calibrated
DNNs. These approaches require training DNN models from scratch and cannot be
used with pre-trained ones. Another category of calibration methods like the
Platt’s scaling [23] and temperature scaling [5] proposes ways of post-
processing the outputs of already trained models to produce calibrated
confidence measures. In [11, 18], it is shown that these methods are not as
well-calibrated as it is reported especially when the validation data are not
independent and identically distributed (IID) and in the presence of
distribution shifts. The Conformal Prediction (CP) framework is developed to
compute prediction sets to satisfy a desired significance level [28, 26, 1].
The confidence value assigned to each possible class is in the form of
$p$-values which is less intuitive than estimating the confidence as
probabilities. Another way of obtaining confidence information about
predictions is by using algorithms based on the Bayesian framework. The use of
this framework, however, require some prior knowledge about the distribution
generating the data. In the real world, this distribution is unknown and it
has to be chosen arbitrarily. In [21], it is shown that the predictive regions
produced by Gaussian Processes, a popular Bayesian machine learning approach,
may be incorrect and misleading when the correct prior is not known.

The main contribution of our work is that we compute low-dimensional,
appropriate, embedding representations of the original inputs in a space where
the Euclidean distance is a measure of similarity between the original inputs,
in order to handle high-dimensional inputs in real-time. Then, we implement
four different taxonomies that split the low-dimensional data into categories
based on their similarity. Last, we present an empirical evaluation of the
approach using two datasets for image classification problems with a large
number of classes as well as detection of botnet attacks in an IoT device. The
underlying models are chosen according to the input size and shape keeping
into account the low-latency and low-power properties to meet the resource
constraints of the variety of use cases [8].

##  II Problem

A perception component in an autonomous system aims to observe and interpret
the environment in order to provide information for decision-making. For
example, a DNN can be used for classifying traffic signs in autonomous
vehicles. The problem is to complement the prediction of the DNN with a
computation of confidence. An efficient and robust approach must ensure a
small and well-calibrated error rate to enable real-time operation. The
approach must ensure a bounded small error rate while limiting the number of
inputs for which an accurate prediction cannot be made.

During system operation, for each new input a prediction is made, usually by a
LEC and the objective is to compute a valid measure of the prediction’s
confidence. The objective is twofold: (1) provide guarantees for the error
rate of the prediction and (2) limit the number of input examples for which a
confident prediction cannot be made. Well-calibrated confidence in terms of
probabilities can be used for decision-making, for example, by generating
warnings when human intervention is required.

The Venn Prediction (VP) framework can produce predictions with well-
calibrated confidence intervals that guarantee to include the true
probabilities for each class output to occur[28]. The confidence intervals for
a test input are generated by considering the class distribution of labeled
inputs assigned to the same category that are collected offline and are
available to the system. In the literature, VP implementations use Support
Vector Machines (SVMs) or DNN classifiers to create categories of labeled
data[28, 13, 20]. The additional problem we are considering is the computation
of appropriate embedding representations that can lead to more efficient VPs.
The main idea is to use distance metric learning and enable DNNs to learn a
lower-dimensional representation for each input on an embedding space where
the Euclidean distance between the input representations is a measure of
similarity between the original inputs themselves. Using such representations
we define taxonomies to form categories of similar input data. This not only
reduces the memory requirements but is also more efficient in producing more
informative intervals.

##  III Probability Intervals based on Distance Metric Learning

Venn Predictors is a machine learning framework that can be combined with
existing classifier architectures for producing well-calibrated multi-
probability predictions under the IID assumption [1, 28]. This means that the
confidence assigned to a prediction is a probability distribution which in
effect defines lower and upper bounds regarding the probability of correctness
for all possible classes. VPs are well-calibrated and the probability bounds
asymptotically contain the corresponding true conditional probabilities (proof
in [28]). However the framework is computationally inefficient as it requires
training the underlying algorithm after every new test input. Computational
efficiency can be addressed using the Inductive Venn Predictors [14, 13], an
extension of the VP framework.

Central to the VP and IVP frameworks is the definition of a Venn taxonomy.
This is a way of clustering data points into a number of categories according
to their similarity and is based on an underlying algorithm. For example a
taxonomy can be defined to put in the same category examples that are
classified in the same class by a DNN. The main idea of our approach is that
the taxonomy can be defined efficiently by learning embedding representations
of the inputs for which the Euclidean distance is a measure of similarity. To
compute the embedding representations of the inputs we train a siamese network
using contrastive loss [6, 9].

We consider the training examples, $z_{1},\dots,z_{l}$ from $\bm{Z}$, where
each $z_{i}$ is a pair $(x_{i},y_{i})$ with $x_{i}$ the feature vector and
$y_{i}$ the corresponding label. We also consider a test input $x_{l+1}$ which
we wish to classify. IVP assumes that all the examples $z_{1},\dots,z_{l+1}$
are independent and identically distributed (IID) generated from the same but
usually unknown probability distribution. The available training examples are
split into two parts: the proper training set with $q$ examples and the
calibration set with $l-q$ examples. The examples in the proper training set
are used to train the siamese network which is used to define different Venn
taxonomies. The roll of the taxonomy is to divide the $l-q$ calibration
examples into a number of categories based on their similarity. This process
takes place during the design time.

After placing the calibration data into categories using the underlying
algorithm for the taxonomy, during execution time we consider a test input
$x_{l+1}$ and place it in a category $k_{l+1}$. The true class $y_{l+1}$ is
unknown and IVP computes a lower and an upper probability
$[L(Y_{j}),U(Y_{j})]$ for every possible class $j=1,\dots,c$ based on the
number of samples of each class in $k_{l+1}$. The predicted class for the
classification is computed as:

$j_{best}=\text{arg}\max_{j=1,\dots,c}\overline{p(Y_{j})}$ (1)

where $\overline{p(Y_{j})}$ is the mean of the probability interval assigned
to $Y_{j}$. Along with the class $Y_{j_{best}}$ the IVP framework outputs the
probability interval $[L(Y_{j_{best}}),U(Y_{j_{best}})]$. The steps taking
place during execution are illustrated in Fig. 1.

##  IV Distance-based Taxonomies

Environment Siamese Network Distance-based Taxonomy Probability Distributions
Set Calibration Data Proper Training Data $x_{l+1}$$r_{l+1}$$k_{l+1}$
$Y_{j_{best}}$ $[L(Y_{j_{best}}),U(Y_{j_{best}})]$ Figure 1: IVP classifier
based on distance metric learning

As proved in[28] the probability intervals assigned to each classification by
the VP are well-calibrated regardless of the choice of the Venn taxonomy and
this holds in practice for IVP as well[13]. However, the choice of the
taxonomy affects the efficiency of the IVP. The probability intervals are
desirable to be relatively narrow to minimize the uncertainty in the
probability of correctness as well as create better separation between the
probabilities of each class. We propose four different Venn taxonomies based
on distance metric learning. The first two taxonomies are based on a
$k$-Nearest Neighbors classifier. The naive approach, that we call k-NN
$V_{1}$, trains a $k$-NN classifier using the embedding representations of the
proper training set. Then the calibration data, as well as each new test
input, are placed to a category that is defined by the $k$-NN prediction using
the computed embedding representations. That is, for a data point $x_{l+1}$
that needs to be placed into a category, its embedding representation is
computed using the siamese network, $r_{l+1}=\text{Net}(x_{l+1})$ and its $k$
nearest training data are found. Depending on the class $\hat{y}_{l+1}$ that
most neighbors belong to, the data point is assigned to the category

$k_{l+1}=\hat{y}_{l+1}.$ (2)

This taxonomy leads to a number of categories that is equal to the number of
classes in the dataset. An extension of the previous taxonomy, k-NN $V_{2}$,
is also based on a $k$-Nearest Neighbors classifier. However, we attempt to
more accurately split the data into categories by taking into account how many
of the $k$ nearest training data points are labeled different than the
predicted class. For a data point $x_{l+1}$ with embedding representation
$r_{l+1}$ that needs to be placed into a category we compute the $k$-nearest
neighbors in the training set and store their labels in a multi-set $\Omega$.
The category where $x_{l+1}$ is placed is computed as:

$k_{l+1}=\hat{y}_{l+1}\times(k-\left\lfloor{\frac{k}{c}}\right\rfloor)+|i\in\Omega:i\neq\hat{y}_{l+1}|$
(3)

where $\hat{y}_{l+1}$ is the $k$-NN classification of $r_{l+1}$, $k$ is the
number of nearest neighbors and $c$ is the number of different classes. This
taxonomy aims at further improving the similarity of the data in each category
leveraging the classifier’s confidence. It is expected that the more similar
labeled neighbor training data points, the higher the chance of the
corresponding class being the correct one. That way each category of k-NN
$V_{1}$ is further split into $k-\left\lfloor{\frac{k}{c}}\right\rfloor$ new
categories.

The ability of siamese networks to create clusters of similar data can be used
to further reduce the Venn taxonomy computational requirements when there is a
large amount of training data. Each class cluster $i$ corresponding to class
$Y_{i}$, $i=1\dots,c$ can then be represented by its centroid
$\mu_{i}=\frac{\sum_{j=1}^{n_{i}}r_{j}^{i}}{n_{i}}$, where $r_{j}^{i}$ is the
embedding representation of the $j^{th}$ training example from class $Y_{i}$
and $n_{i}$ is the number of training examples labeled as $Y_{i}$. We propose
another family of taxonomies based on the Nearest Centroids. The NC $V_{1}$
places the calibration data as well as each new test input to a category that
is the same as the class assigned to their nearest centroid. The category
where an example $x_{l+1}$ is placed is computed as:

$k_{l+1}=\text{arg}\min_{j=1,\dots,c}d(r_{l+1},\mu_{j})$ (4)

where $d$ the Euclidean distance. This leads to a number of categories that is
equal to the number of classes in the dataset. An extension of this taxonomy,
the NC $V_{2}$, attempts to form more accurate categories by taking into
account the classification confidence. We expect data points of the same class
to be more similar to each other when their embedding representations are
placed at similar distances to their class centroid. That way each category of
NC $V_{1}$ is further split into two categories based on how close an example
$x_{l+1}$ is to its nearest centroid:

$k_{l+1}=2\times\text{arg}\min_{j=1,\dots,c}d(r_{l+1},\mu_{j})+h,$ (5)

$h=\begin{cases}0,&\text{if }d(r_{l+1},\mu_{\text{min}})\leq\theta\\\
1,&\text{otherwise}\end{cases}$

where $\mu_{\text{min}}=\text{arg}\min_{j=1,\dots,c}d(r_{l+1},\mu_{j})$ is the
distance to the nearest centroid and $\theta$ a chosen distance threshold.

##  V Evaluation Metrics

The performance of IVP based on the proposed taxonomies is evaluated regarding
the accuracy, calibration and efficiency. The objective is for the computed
probability intervals to contain the true probability of correctness for each
prediction. The probability interval for a given input $x$ with predicted
class $\hat{y}$ is $[L(\hat{y}),U(\hat{y})]$. Equivalently, the probability
that $\hat{y}$ is not the correct classification will be in the complimentary
interval $[1-U(\hat{y}),1-L(\hat{y})]$, called error probability interval. The
true probability of correctness for a single prediction is unknown so the
correctness of the computed intervals is evaluated over a number of samples.
To do this we use the following metrics:

  * •

cumulative errors

$E_{n}=\sum_{i=1}^{n}err_{i},$ (6) $err_{i}=\begin{cases}1,&\text{if
}\text{classification }\hat{y}_{i}\text{ is incorrect}\\\
0,&\text{otherwise}\end{cases}$

  * •

cumulative lower and upper error probabilities

$\displaystyle LEP_{n}=\sum_{i=1}^{n}[1-U(\hat{y})],$ $\displaystyle
UEP_{n}=\sum_{i=1}^{n}[1-L(\hat{y})]$ (7)

To compare the IVP implementations based on our proposed taxonomies with the
baseline taxonomies, scalar metrics are used that represent the performance
regarding accuracy, calibration, and efficiency. Unlike the NN classifiers
that produce a single softmax probability for each class, the IVP framework
produces probability intervals. For the computation of the evaluation metrics
the probability assigned to a class $Y_{j}$ will be $\overline{p(Y_{j})}$ like
in (1). The accuracy of an IVP implementation is evaluated as the number of
correct classifications over the number of attempted classifications and it is
computed as

$accuracy=1-\dfrac{E_{n}}{n}.$ (8)

An efficient, or informative, IVP is one that makes predictions with small
diameter probability intervals and their median is as close to zero or one.
The most popular quality metrics for probability assessments are the negative
log-likelihood (NLL) and the Brier score (BS) [4]. NLL is the simplest out of
the two and only considers the probability assigned to the predicted class in
(1). It is computed as

$NLL=-\sum_{i=1}^{n}\sum_{j=1}^{c}t_{i}^{j}\log(o_{i}^{j}),$ (9)

where $o_{i}^{j}=\overline{p(Y_{j})}$ of example $i$ and $t_{i}^{j}$ the one-
hot representation of the ground truth classification label $y_{i}$ of example
$i$, that is

$t_{i}^{j}=\begin{cases}1,&\text{if }\text{classification }y_{i}=Y_{j}\\\
0,&\text{otherwise}\end{cases}$

This metric is minimized by producing intervals that are narrow and have
median probability close to one assigned to the correct class. Computational
issues may occur as the log score explodes if we observe an event that the
classifier considers impossible. BS is computed as

$BS=\dfrac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{c}(o_{i}^{j}-t_{i}^{j})^{2}$ (10)

This is, in effect, the mean squared error of the predictions. Unlike NLL, BS
considers the probabilities assigned to all possible classes and will penalize
probability intervals assigned to incorrect classes that are not close to
zero. There are different views in the literature regarding which scoring rule
is more appropriate. [2] emphasizes in the importance of the locality
property, meaning, the scoring rule should only depend on the probability of
events that actually occur and only NLL satisfies this. On the other hand,
[25] states that a scoring rule should be symmetric and only BS satisfies
this. This means that if the true class probability is $p$ and the predicted
probability is $\hat{p}$, then the score should be equal to the case where the
true probability is $\hat{p}$ and the predicted probability is $p$. However,
we think that both metrics produce useful insights in probability assessment
so both are reported in our experiment results. The interval size has a
significant role on how informative and interpretable a prediction is. We
evaluate the size of the probability intervals by computing the average
interval diameter as

$D=\dfrac{\sum_{i=1}^{n}U(\hat{y})-\sum_{i=1}^{n}L(\hat{y})}{n}$ (11)

A well-calibrated IVP computes probability intervals that are representative
of the true correctness likelihood. Formally a model is well-calibrated when

$\mathbb{P}\left(\hat{y}=Y|\hat{p}=p\right)=p,\quad\quad\forall p\in[0,1]$
(12)

However, $\hat{p}$ is a continuous random variable so the probability in (12)
cannot be approximated using finitely many samples. According to (12) a
measure of miscalibration can be expressed as
$\underset{\hat{p}}{\mathbb{E}}\left[\left|\mathbb{P}\left(\hat{y}=y|\hat{p}=p\right)-p\right|\right]$.
The Expected Calibration Error (ECE) [17] computes an approximation of this
expected value across bins:

$\text{ECE}=\sum_{m=1}^{M}\dfrac{|B_{m}|}{n}\left|\text{acc}(B_{m})-\text{conf}(B_{m})\right|$
(13)

where $|B_{m}|$ is the number of samples in bin $B_{m}$, $n$ is the total
number of samples and acc$(B_{m})$ and conf$(B_{m})$ are the accuracy and
confidence of bin $B_{m}$ respectively as defined in [17]. Many times in
safety critical applications it is more useful to compute the maximum
miscalibration of a model than the mean value. This metric is called Maximum
Calibration Error (MCE) [17] and is computed as:

$\text{MCE}=\max_{m\in\\{1,\dots,M\\}}\left|\text{acc}(B_{m})-\text{conf}(B_{m})\right|$
(14)

##  VI Evaluation

In this section, we evaluate the IVPs that use distance-based taxonomies with
regard to accuracy, calibration, and efficiency. Additionally, for the
evaluation of our proposed taxonomies, we use metrics regarding the
performance of the siamese network in clustering similar input data, the
execution time of the framework, and the required memory.

###  VI-A Experimental Setup

The embedding representation computations, part of our proposed taxonomies,
are not application-specific and can improve the performance of IVP in cases
where inputs are high-dimensional. We evaluate the performance of IVP with
distance-learning in two different classification problems. First, we have two
case studies in image classification. The German Traffic Sign Recognition
Benchmark (GTSRB) dataset is a collection of traffic sign images to be
classified in 43 classes [27]. The labeled sign images are of various sizes
between 15x15 to 250x250 pixels depending on the observed distance. We convert
all the images to a fixed shape of 96x96 pixels. The second dataset is the
Fruits360 [16]. This dataset contains images of 131 different kinds of fruits
and vegetables. The input data are used in their original size, 100x100
pixels.

The second classification problem we consider is the detection of botnet
attacks in IoT devices. As part of the evaluation in [15], authors made
available data regarding network traffic while infecting different common IoT
devices two families of botnets. Mirai and BASHLITE are two common IoT-based
botnets and their harmful capabilities are presented in [10]. In the dataset
there are data for the following ten attacks:

  * •

BASHLITE Attacks

    1. 1.

Scan: Scanning the network for vulnerable devices

    2. 2.

Junk: Sending spam data

    3. 3.

UDP: UDP flooding

    4. 4.

TCP: TCP flooding

    5. 5.

COMBO: Sending spam data and opening a connection to a specified IP address
and port

  * •

Mirai Attacks

    1. 1.

Scan: Scanning the network for vulnerable devices

    2. 2.

Ack: Ack flooding

    3. 3.

Syn: Syn flooding

    4. 4.

UDP: UDP flooding

    5. 5.

UDPplain: UDP flooding with fewer options, optimized for higher PPS

Including the benign network traffic we approach this as a classification
problem with eleven classes. The available data are in the form of 115
statistical features extracted from the raw network traffic. The same 23
features, presented in [15], are extracted from five time windows of the most
recent 100ms, 500ms, 1.5sec, 10sec, 1min. The features summarize the traffic
in each of these time windows that has (1) the same source IP address, (2) the
same source IP and MAC address, (3) been sent between the source and
destination IP address, (4) been sent between the source and destination
TCP/UDP sockets. These features are computed incrementally and in real-time.

The available data are used throughout the evaluation process the same way in
every dataset. 10% of the data are taken out to be used for testing and the
rest is the training set. The training set is then split into the proper
training set and the calibration set. The proper training set is randomly
chosen as 80% of the training set and is used to train the underlying models
and for the computation of the categories. The calibration set is the
remaining 20% of the available training data is used only to form the
categories during the design time. The reported evaluation results are
computed on the separate test set. All the experiments run in a desktop
computer equipped with and Intel(R) Core(TM) i9-9900K CPU and 32 GB RAM and a
Geforce RTX 2080 GPU with 8 GB memory.

###  VI-B Baseline

To understand the effect of the distance metric learning in IVP we compare it
with approaches that use DNN classifiers as underlying algorithms. A variety
of Venn taxonomy definitions based on DNNs is proposed in [20]. $V_{1}$
assigns two examples to the same category if their maximum softmax outputs
correspond to the same class. $V_{2}$, divides the examples in the categories
defined by $V_{1}$ into two smaller categories based on the value of their
maximum softmax output. Their chosen threshold for the maximum output to
create the two smaller categories is $0.75$. $V_{3}$ divides the examples in
the categories defined by $V_{1}$ into two smaller categories but this time
based on the second highest softmax output. Their chosen threshold for the
second-highest output is $0.25$. $V_{4}$ divides each category of taxonomy
$V_{1}$ in two, based on the difference between the highest and second-highest
softmax outputs. The threshold for this difference is $0.5$. In the same
paper, they proposed a fifth taxonomy that creates the categories based on
which classes have softmax outputs above a certain threshold. This taxonomy
creates $2^{C}$ number of categories making its use infeasible in our
evaluation datasets.

###  VI-C Evaluation Results

The difficulty to assign an input to a category and the memory demands
increase as the size and complexity of the inputs increases. Our goal is to
evaluate our method using general-purpose and lightweight DNNs. For the image
classification problems, we use the MobileNet architecture for both the
embedding representation computation as well as the classifier used for the
baseline taxonomies for its low latency and low memory requirements. The
trade-off between accuracy and latency is configured by the hyperparameter
$\alpha$. We set $\alpha=0.5$ in the case of GTSRB and $\alpha=1$ for the
Fruit360. In both cases the embedding representation vectors are of size 128.
In the case of the botnet attacks detection, the input data are arranged in
vectors of 115 values so we use a fully connected DNN with two hidden layers,
the first has 10 units, and the second which produces the embedding
representations has 32 units.

After training the siamese network and before it is used as part of the
taxonomies we need to evaluate how well it performs in clustering similar
inputs. For comparison, we use the embedding space produced by the penultimate
layer of the DNN classifier [7]. A commonly used metric of the separation
between class clusters is the silhouette coefficient [24]. This metric
evaluates how close together samples from the same class are, and far from
samples of different classes and takes values in [-1,1]. The results on the
silhouette analysis for the test inputs from both datasets are shown in Table
I. The siamese network produces representations that are well clustered based
on their similarity and better than the representations produced by the
classifier DNN. This is important for constructing efficient categories using
our proposed distance-based taxonomies.

Table I: Silhouette Coefficient Comparison | Classifier Embeddings | Siamese Embeddings
---|---|---
GTSRB | 0.56 | 0.98
Fruits360 | 0.52 | 0.85
Ecobee Thermostat | 0.27 | 0.46

For illustration, the cumulative upper and lower error probabilities as well
as the cumulative error are plotted on the same axis in Fig. 2 using the NC
$V_{2}$ taxonomy and test data from the GTSRB dataset. The computed
probability intervals successfully bound the true error-rate.

$0$$500$$1{,}000$$1{,}500$$2{,}000$$2{,}500$$0$$10$$20$$30$$40$sample #Cumulative errorLEPUEPCE Figure 2: Illustration of cumulative metrics in the GTSRB dataset using the NC $V_{2}$ taxonomy Table II: Evaluation metrics results Dataset | Taxonomy | Accuracy | NLL | BS | D | ECE | MCE | Time | Memory
---|---|---|---|---|---|---|---|---|---
GTSRB | $V_{1}$ | 0.994 | 111.835 | 0.013 | 3.29e-06 | 0.005 | 0.005 | 3.6ms | 11.2MB
$V_{2}$ | 0.992 | 58.104 | 0.055 | 5.30e-06 | 0.011 | 0.583 | 6.6ms | 11.2MB
$V_{3}$ | 0.993 | 75.394 | 0.038 | 4.48e-06 | 0.008 | 0.750 | 3.7ms | 11.2MB
$V_{4}$ | 0.991 | 70.279 | 0.053 | 5.06e-06 | 0.009 | 0.750 | 2.9ms | 11.2MB
$k$-nn $V_{1}$ | 0.998 | 41.575 | 0.005 | 3.30e-06 | 0.004 | 0.004 | 3.2ms | 19MB
$k$-nn $V_{2}$ | 0.998 | 41.126 | 0.005 | 3.30e-06 | 0.004 | 0.004 | 3.6ms | 19.8MB
NC $V_{1}$ | 0.998 | 41.575 | 0.005 | 3.30e-06 | 0.004 | 0.004 | 2.9ms | 3.9MB
NC $V_{2}$ | 0.996 | 38.444 | 0.046 | 6.21e-06 | 0.007 | 0.500 | 3ms | 3.9MB
Fruits360 | $V_{1}$ | 0.983 | 1089.938 | 0.043 | 1.19e-06 | 0.008 | 0.113 | 4.4ms | 41MB
$V_{2}$ | 0.986 | 816.893 | 0.144 | 1.57e-06 | 0.013 | 0.407 | 4ms | 41.2MB
$V_{3}$ | 0.985 | 870.470 | 0.154 | 1.55e-06 | 0.012 | 0.392 | 4.6ms | 41.2MB
$V_{4}$ | 0.985 | 836.295 | 0.159 | 1.58e-06 | 0.013 | 0.384 | 2.7ms | 41.2MB
$k$-nn $V_{1}$ | 0.993 | 532.314 | 0.025 | 1.19e-06 | 0.010 | 0.073 | 3.3ms | 127.5MB
$k$-nn $V_{2}$ | 0.993 | 466.311 | 0.088 | 1.42e-06 | 0.011 | 0.243 | 3.7ms | 128.1MB
NC $V_{1}$ | 0.991 | 605.087 | 0.027 | 1.19e-06 | 0.010 | 0.045 | 3.6ms | 14MB
NC $V_{2}$ | 0.988 | 725.556 | 0.208 | 2.22e-06 | 0.018 | 0.500 | 3.5ms | 14.2MB
Ecobee Thermostat | $V_{1}$ | 0.823 | 4732.483 | 0.218 | 2.67e-08 | 0.003 | 0.009 | 0.7ms | 52.2kB
$V_{2}$ | 0.830 | 4310.008 | 0.200 | 4.20e-08 | 0.003 | 0.014 | 0.7ms | 53.2kB
$V_{3}$ | 0.830 | 4311.460 | 0.200 | 4.33e-08 | 0.002 | 0.015 | 0.7ms | 53.2kB
$V_{4}$ | 0.830 | 4306.791 | 0.200 | 4.26e-08 | 0.003 | 0.040 | 0.6ms | 53.2kB
$k$-nn $V_{1}$ | 0.935 | 2872.725 | 0.113 | 2.94e-08 | 0.001 | 0.003 | 1.9ms | 43.8MB
$k$-nn $V_{2}$ | 0.935 | 2299.023 | 0.096 | 1.33e-07 | 0.006 | 0.375 | 2.4ms | 43.8MB
NC $V_{1}$ | 0.794 | 5550.013 | 0.255 | 2.78e-08 | 0.006 | 0.017 | 1ms | 24kB
NC $V_{2}$ | 0.794 | 5541.171 | 0.255 | 4.02e-08 | 0.006 | 0.023 | 0.9ms | 25kB

The evaluation results are shown in Table II. For both datasets, we observe
that using the proposed distance-based taxonomies, IVP produces more accurate
classifications. Even though the baseline $V_{1}$ taxonomy produces
probability intervals that are as narrow as the intervals produced by some of
the proposed taxonomies, the proposed taxonomies produce better quality
intervals by keeping the intervals assigned to the correct class close to 1
and the intervals of the incorrect classes close to 0, as shown by the NLL and
BS metrics. The differences in ECE are not significant but most of the
proposed taxonomies produce probabilities that are better calibrated in the
whole probability space $[0,1]$ with no areas of miscalibration as indicated
by MCE.

The times required for the computation of a classification and the probability
intervals when a new input arrives are similar in both the baseline and our
proposed taxonomies and indicate they can be used for real-time operation. The
speed bottleneck is the computations by the DNNs for either the
classifications or the representation mapping. The $k$-NN computation step in
the low-dimensional embedding representation space adds minimal overhead in
the execution time. The memory requirements have two main parts: the memory
required to store the DNN weights and the memory required to store the
categories after calibration. The proposed taxonomies have the additional
requirement to store either the embedding representations of the training data
to be used by the $k$-NN or the centroid of each class. The representations of
the training data are stored in a $k-d$ tree [3] for fast $k$-NN computation.
With the use of low-dimensional representations, the additional memory
required for the nearest centroid based taxonomies is small compared to the
underlying DNN size.

##  VII Conclusion

Although DNNs offer advanced capabilities, they must be complemented by
engineering methods and practices for them to provide accurate measures of
prediction confidence. For classification tasks, the IVP framework computes
probability intervals that contain the probability of the prediction’s
correctness by examining the underlying model’s accuracy on similar data. We
presented computationally efficient algorithms based on appropriate embedding
representations learned by siamese networks that make it possible for IVP to
be used with high-dimensional data for real-time applications. The evaluation
results demonstrate that the IVP framework using distance-based taxonomies
produces high accuracy and probability intervals that are efficient and well-
calibrated. Our choice of lightweight DNNs and small embedding representation
size make the approach computationally efficient and can be used in real-time.
A direction for future extension of this work is to improve the probability
intervals, regarding their efficiency, during execution time.

## Acknowledgment

The material presented in this paper is based upon work supported by the
Defense Advanced Research Projects Agency (DARPA) through contract number
FA8750-18-C-0089. The views and conclusions contained herein are those of the
authors and should not be interpreted as necessarily representing the official
policies or endorsements, either expressed or implied, of DARPA.

## References

  * [1] V. Balasubramanian, S.-S. Ho, and V. Vovk.  Conformal Prediction for Reliable Machine Learning: Theory, Adaptations and Applications.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1st edition, 2014\.
  * [2] R. Benedetti.  Scoring rules for forecast verification.  Monthly Weather Review, 138(1):203–211, 2010.
  * [3] J. L. Bentley.  Multidimensional binary search trees used for associative searching.  Commun. ACM, 18(9):509–517, Sept. 1975.
  * [4] G. W. Brier.  Verification of forecasts expressed in terms of probability.  Monthly Weather Review, 78(1):1 – 3, 01 Jan. 1950.
  * [5] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger.  On calibration of modern neural networks.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, pages 1321–1330. JMLR.org, 2017.
  * [6] R. Hadsell, S. Chopra, and Y. LeCun.  Dimensionality reduction by learning an invariant mapping.  In 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), volume 2, pages 1735–1742. IEEE, 2006.
  * [7] G. E. Hinton.  Learning multiple layers of representation.  Trends in cognitive sciences, 11(10):428–434, 2007.
  * [8] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861, 2017.
  * [9] G. Koch, R. Zemel, and R. Salakhutdinov.  Siamese neural networks for one-shot image recognition.  In ICML deep learning workshop, volume 2. Lille, 2015.
  * [10] C. Kolias, G. Kambourakis, A. Stavrou, and J. Voas.  Ddos in the iot: Mirai and other botnets.  Computer, 50(7):80–84, 2017.
  * [11] A. Kumar, P. S. Liang, and T. Ma.  Verified uncertainty calibration.  In Advances in Neural Information Processing Systems, pages 3792–3803, 2019.
  * [12] A. Kumar, S. Sarawagi, and U. Jain.  Trainable calibration measures for neural networks from kernel mean embeddings.  In J. Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2805–2814, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018. PMLR.
  * [13] A. Lambrou, I. Nouretdinov, and H. Papadopoulos.  Inductive venn prediction.  Annals of Mathematics and Artificial Intelligence, 74(1-2):181–201, 2015.
  * [14] A. Lambrou, H. Papadopoulos, I. Nouretdinov, and A. Gammerman.  Reliable probability estimates based on support vector machines for large multiclass datasets.  In IFIP International Conference on Artificial Intelligence Applications and Innovations, pages 182–191. Springer, 2012.
  * [15] Y. Meidan, M. Bohadana, Y. Mathov, Y. Mirsky, A. Shabtai, D. Breitenbacher, and Y. Elovici.  N-baiot—network-based detection of iot botnet attacks using deep autoencoders.  IEEE Pervasive Computing, 17(3):12–22, 2018.
  * [16] H. Mureşan and M. Oltean.  Fruit recognition from images using deep learning.  arXiv preprint arXiv:1712.00580, 2017.
  * [17] M. P. Naeini, G. Cooper, and M. Hauskrecht.  Obtaining well calibrated probabilities using bayesian binning.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 29, 2015.
  * [18] Y. Ovadia, E. Fertig, J. Ren, Z. Nado, D. Sculley, S. Nowozin, J. Dillon, B. Lakshminarayanan, and J. Snoek.  Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift.  In Advances in Neural Information Processing Systems, pages 13991–14002, 2019.
  * [19] H. Papadopoulos.  Reliable probabilistic prediction for medical decision support.  In Artificial Intelligence Applications and Innovations, pages 265–274. Springer, 2011.
  * [20] H. Papadopoulos.  Reliable probabilistic classification with neural networks.  Neurocomputing, 107:59–68, 2013.
  * [21] H. Papadopoulos, V. Vovk, and A. Gammerman.  Regression conformal prediction with nearest neighbours.  J. Artif. Int. Res., 40(1):815–840, Jan. 2011.
  * [22] G. Pereyra, G. Tucker, J. Chorowski, Ł. Kaiser, and G. Hinton.  Regularizing neural networks by penalizing confident output distributions.  arXiv preprint arXiv:1701.06548, 2017.
  * [23] J. C. Platt.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  In ADVANCES IN LARGE MARGIN CLASSIFIERS, pages 61–74. MIT Press, 1999.
  * [24] P. J. Rousseeuw.  Silhouettes: A graphical aid to the interpretation and validation of cluster analysis.  Journal of Computational and Applied Mathematics, 20:53 – 65, 1987\.
  * [25] R. Selten.  Axiomatic characterization of the quadratic scoring rule.  Experimental Economics, 1(1):43–61, 1998.
  * [26] G. Shafer and V. Vovk.  A tutorial on conformal prediction.  J. Mach. Learn. Res., 9:371–421, June 2008.
  * [27] J. Stallkamp, M. Schlipsing, J. Salmen, and C. Igel.  Man vs. computer: Benchmarking machine learning algorithms for traffic sign recognition.  Neural Networks, 32:323 – 332, 2012.  Selected Papers from IJCNN 2011.
  * [28] V. Vovk, A. Gammerman, and G. Shafer.  Algorithmic Learning in a Random World.  Springer-Verlag, Berlin, Heidelberg, 2005.
  * [29] V. Vovk, G. Shafer, and I. Nouretdinov.  Self-calibrating probability forecasting.  In NIPS, pages 1133–1140, 2003.