obtained via the following marginalisation,

$p(\tilde{\mathbf{x}}_{t}|y_{t}=k,\mathcal{D}_{l})=\int\int
p(\tilde{\mathbf{x}}_{t}|\bm{\mu}_{k},\Sigma_{k})p(\bm{\mu}_{k},\Sigma_{k}|y_{t}=k,\mathcal{D}_{l})d\bm{\mu}_{k}d\Sigma_{k}$
(A.13)

resulting in the Student-t distribution [46],

$(\tilde{\mathbf{x}}_{t}|y_{t}=k,\mathcal{D}_{l})\sim\mathcal{T}\bigg{(}\bm{m}_{n},\frac{\kappa_{n}+1}{\kappa_{n}(v_{n}-D+1)}S_{n},v_{n}-D+1\bigg{)}$
(A.14)

where $\bm{m}_{n}$, $\kappa_{n}$, $v_{n}$, $S_{n}$ are the updated
hyperparameters and $D$ is the dimensionality of the feature space. Here, the
first two parameters of the Student-t distribution correspond to the mean and
scale, respectively. The third parameter specifies the degrees of freedom and
relates to the probability mass found in the tails of the distribution. The
full functional form of the Student-t distribution can be found in [46].

Via a similar procedure, the posterior predictive distribution over the labels
is obtained via the following marginalisation,

$p(\tilde{y}_{t}|\mathcal{D}_{l})=\int
p(\tilde{y}_{t}|\bm{\lambda})p(\bm{\lambda}|\mathcal{D}_{l})d\bm{\lambda}$
(A.15)

resulting in,

$p(\tilde{y}_{t}=k|\mathcal{D}_{l})=\frac{n_{k}+\alpha_{k}}{n+\alpha_{0}}$
(A.16)

where $n=\sum_{k=1}^{K}n_{k}$ and $\alpha_{0}=\sum_{k=1}^{K}\alpha_{k}$.

Finally, the predictive distribution for the class labels given a new
unlabelled observation $\tilde{\mathbf{x}}_{t}$ can be obtained using the
posterior predictive distribution and applying Bayes’ rule [7],

$p(\tilde{y}_{t}=k|\tilde{\mathbf{x}}_{t},\mathcal{D}_{l})=\frac{p(\tilde{\mathbf{x}}_{t}|\tilde{y}_{t}=k,\mathcal{D}_{l})p(\tilde{y}_{t}=k|\mathcal{D}_{l})}{p(\tilde{\mathbf{x}}_{t}|\mathcal{D}_{l})}$
(A.17)

##  Appendix B Expectation-Maximisation for Gaussian Mixture Models

The aim of semi-supervised learning via the EM algorithm is to infer updated
distribution parameters $\bm{\Theta}$ from
$\mathcal{D}=\mathcal{D}_{l}\cup\mathcal{D}_{u}$. The maximum a posteriori MAP
estimate of these updated parameters is specified as follows,

$\bm{\Theta}|\mathcal{D}=\text{argmax}_{\bm{\Theta}}\biggl{\\{}\frac{p(\mathcal{D}|\bm{\Theta})p(\bm{\Theta})}{p(\mathcal{D})}\biggr{\\}}=\text{argmax}_{\bm{\Theta}}\biggl{\\{}\frac{p(\mathcal{D}_{u}|\bm{\Theta})p(\mathcal{D}_{l}|\bm{\Theta})p(\bm{\Theta})}{p(\mathcal{D}_{u},\mathcal{D}_{l})}\biggr{\\}}$
(B.1)

Implicit in the factorisation of $p(\mathcal{D}|\bm{\Theta})$ in equation
(B.1), is the assumption that $\mathcal{D}_{l}$ and $\mathcal{D}_{u}$ are
conditionally independent. This assumption holds for random querying, as
samples selected in this manner are independent and identically distributed
(iid). Unfortunately, active learning violates this assumption, as data in
$\mathcal{D}_{l}$ are not iid because of the preferential querying process and
iterative model updating [20, 11]. Convenient assumptions such as this are
frequently relied upon in statistical and engineering analyses (particularly
in active learning contexts [7]). As such, the assumption is embraced for the
current case study in order to demonstrate the decision-making performance
that one may achieve despite the violation.

To circumvent numerical instabilities, the MAP estimate in equation (B.1) is
formulated as a maximisation of the expected joint log-likelihood across
$\mathcal{D}$ [26],

$\mathcal{L}(\bm{\Theta}|\mathcal{D}_{l},\mathcal{D}_{u})\propto\sum_{t=1}^{m}\log\sum_{k=1}^{K}p(\tilde{\mathbf{x}}_{t}|y_{t}=k,\bm{\Theta})p(y_{t}=k|\bm{\Theta})\\\
+\sum_{t=1}^{n}\log[p(\mathbf{x}|y_{t}=k)p(y_{t}=k|\bm{\Theta})]+\log
p(\bm{\Theta})$ (B.2)

where the first term corresponds to the log-likelihood of the model over
$\mathcal{D}_{u}$, the second term corresponds to the log-likelihood of the
model over $\mathcal{D}_{l}$ and the final term is the log-prior-likelihood of
the model parameters. It should be noted that the first term contains a
summation over the label space; this marginalises out $y_{t}$, which is
considered a latent variable for data in $\mathcal{D}_{u}$.

During each E-step of the EM algorithm, the unlabelled observations are
classified using the current estimate of the model parameters. During the
M-step, updates for parameters $\bm{\Theta}$ are found using the predicted
labels determined via the E-step, in addition to the acquired labelled data.

The E-step and M-step are more formally defined as follows [29].

E-step: A responsibility matrix $R$ is computed for the unlabelled data
corresponding to the posterior label prediction given by equation (A.17),

$R[t,k]=r_{tk}=p(\tilde{y}_{t}=k|\tilde{\mathbf{x}}_{t},\bm{\Theta}),\quad\forall\tilde{\mathbf{x}}_{t}\in\mathcal{D}_{u}$
(B.3)

The posterior label predictions for observations in $\mathcal{D}_{l}$ are
given by the known labels $y_{t}$, and can be represented using discrete delta
functions [48],

$p(y_{t}=k|\mathbf{x}_{t})=\delta_{k,y_{t}},\quad\forall(\mathbf{x}_{t},y_{t})\in\mathcal{D}_{l}$
(B.4)

where $\delta_{k,y_{t}}$ is the Kronecker delta function – equal to 1 when
$y_{t}=k$, and 0 otherwise.

Whereas, for $\mathcal{D}_{l}$, the number of observations corresponding to
each class $n_{k}$ are known, these values are indeterminate for data in
$\mathcal{D}_{u}$. As such, it is convenient to define for each class their
effective counts $r_{k}$ in $D_{u}$; these may be calculated from the
responsibility matrix as follows [46],

$r_{k}=\sum_{t=1}^{m}r_{tk}$ (B.5)

The total (effective) counts per class over $\mathcal{D}$ can be summarised
as,

$N_{k}=n_{k}+r_{k}$ (B.6)

M-step: Updates for $\bm{\Theta}$ are computed via modified versions of
equations (A.8) to (A.11) and equation (A.16). The mean and covariance
parameters are updated as follows,

$\bm{m}_{n}=\frac{\kappa_{0}}{\kappa_{0}+N_{k}}\bm{m}_{0}+\frac{N_{k}}{\kappa_{0}+N_{k}}\bar{\mathbf{x}}_{k}$
(B.7)

$\bar{\mathbf{x}}_{k}\overset{\Delta}{=}\frac{\sum_{t=1}^{n}\delta_{k,y_{t}}\mathbf{x}_{t}+\sum_{t=1}^{m}r_{tk}\tilde{\mathbf{x}}_{t}}{N_{k}}$
(B.8)

$\kappa_{n}=\kappa_{0}+N_{k}$ (B.9)

$v_{n}=v_{0}+N_{k}$ (B.10)

$S_{n}=S_{0}+S_{k}+\kappa_{0}\bm{m}_{0}\bm{m}_{0}^{\top}-\kappa_{n}\bm{m}_{n}\bm{m}_{n}^{\top}$
(B.11)

$S_{k}\overset{\Delta}{=}\sum_{t=1}^{n}\delta_{k,y_{t}}\mathbf{x}_{t}\mathbf{x}_{t}^{\top}+\sum_{t=1}^{m}r_{tk}\tilde{\mathbf{x}}_{t}\tilde{\mathbf{x}}_{t}^{\top}$
(B.12)

which lead to MAP estimates given by,

$\hat{\bm{\mu}}_{k}=\bm{m}_{n}$ (B.13)

$\hat{\Sigma}_{k}=\frac{S_{n}}{v_{n}+D+2}$ (B.14)

The mixing proportions $\bm{\lambda}$ are updated as follows,

$\frac{\alpha_{k}+N_{k}-1}{\alpha_{0}+N-K}$ (B.15)

where $N=|\mathcal{D}|=n+m$.

##  Appendix C Multiclass Relevance Vector Machines

For the multiclass classification of an unlabelled data point
$\tilde{\mathbf{x}}_{t}$, mRVMs employ a fundamentally (generalised) linear
model-form as a foundation,

$\mathbf{f}_{t}=\mathbf{W}^{\top}\mathbf{k}(\mathbf{x}_{l},\tilde{\mathbf{x}}_{t})$
(C.1)

where $\mathbf{f}_{t}=\\{f_{1},\ldots,f_{K}\\}^{\top}$ is a vector of $K$
auxiliary variables that provide a ranking system by which the class
membership of an unlabelled data point may be assessed.
$\mathbf{k}(\mathbf{x}_{l},\tilde{\mathbf{x}}_{t})$ is an $n\times 1$ vector
for which the $i^{\text{th}}$ element is specified by the kernel function
$k(\mathbf{x}_{i},\tilde{\mathbf{x}}_{t})$. The kernel function specifies a
set of basis functions that reflect the similarity between
$\tilde{\mathbf{x}}_{t}$ and training inputs
$\mathbf{x}_{l}=\\{\mathbf{x}_{i}|(\mathbf{x}_{i},y_{i})\in\mathcal{D}_{l}\\}^{n}_{i=1}$;
nonlinearity can be introduced into the RVM by selecting a nonlinear kernel
function. For compactness, the vector
$\mathbf{k}(\mathbf{x}_{l},\tilde{\mathbf{x}}_{t})$ is denoted as
$\mathbf{k}_{t}$ herein.
$\mathbf{W}=\\{\mathbf{w}_{1},\ldots,\mathbf{w}_{K}\\}$ is an $n\times K$
matrix of tunable parameters referred to as weights, and where
$\mathbf{w}_{k}=\\{w_{1,k},\ldots,w_{n,k}\\}^{\top}$. These weights act as a
voting system that indicate which data in $\mathcal{D}_{l}$ are important, or
‘relevant’, for discriminating between classes.

In accordance with [12], the auxiliary variables in $\mathbf{f}_{t}$ are
assumed to adhere to a standardised noise model:

$f_{k}|\mathbf{w}_{k},\mathbf{k}_{t}\sim\mathcal{N}(\mathbf{w}_{k}^{\top}\mathbf{k}_{t},1)$
(C.2)

Predicted class labels $\hat{y}_{t}$ are assigned to otherwise unlabelled data
via the auxiliary variables $\mathbf{f}_{t}$ by using a criterion specified by
the multinomial probit link,

$\hat{y}_{t}=k\iff f_{k}>f_{j}\text{ }\forall j\neq k$ (C.3)

Moreover, a probabilistic representation of class membership can be obtained
via the following marginalisation,

$p(\tilde{y}_{t}=k|\mathbf{W},\mathbf{k}_{t})=\int
p(\tilde{y}_{t}=k|\mathbf{f}_{t})p(\mathbf{f}_{t}|\mathbf{W},\mathbf{k}_{t})d\mathbf{f}_{t}$
(C.4)

where $p(\tilde{y}_{t}=k|\mathbf{f}_{t})=\delta_{k,\hat{y}_{t}}$. This
marginalisation yields the multinomial probit likelihood [25],

$p(\tilde{y}_{t}=k|\mathbf{W},\mathbf{k}_{t})=\mathbb{E}_{p(u)}\Biggl{[}\prod_{j\neq
k}\Phi(u+(\mathbf{w}_{k}-\mathbf{w}_{j})^{\top}\mathbf{k}_{t})\Biggr{]}$ (C.5)

where $u\sim\mathcal{N}(0,1)$ and $\Phi$ denotes the Gaussian cumulative
distribution function. Here, it is worth noting that, in practice, equation
(C.5) cannot be computed analytically and is instead approximated via Gauss-
Hermite quadrature [12].

Having specified the predictive model, the relevant supervised learning
problem for mRVMs can be expressed as the Bayesian inference of weights
$\mathbf{W}$ from $\mathcal{D}_{l}$. In order to conduct Bayesian inference, a
prior must first be placed upon the model parameters $\mathbf{W}$; in this
case,

$w_{i,k}\sim\mathcal{N}(0,\alpha_{i,k}^{-1})$ (C.6)

where the scale parameters $\alpha_{i,k}$ can be summarised in an $n\times K$
matrix $\mathbf{A}$ and are assigned the following hyperprior,

$\alpha_{i,k}\sim\Gamma(\tau,\nu)$ (C.7)

where $\tau$ and $\nu$ are hyperparameters related to the shape and scale,
respectively. Given small values for the hyperparameters $\tau$ and $\nu$, the
prior and hyperprior presented above result in a Student-$t$ distribution with
zero-mean and small variance over the weights in $\mathbf{W}$. This
restrictive distribution causes few weights to be non-zero, thereby inducing
sparsity in the model [34, 12].

The training of the mRVM model is accomplished by updating the model
parameters via an EM algorithm [29]. For a detailed exposition of the learning
procedure, the reader is directed to [35]. Here, the EM steps are provided.

The model weight parameters are updated as the MAP estimate of the posterior
distribution over the weights; given by
$\hat{\mathbf{W}}=\operatorname*{argmax}_{\mathbf{W}}p(\mathbf{W}|\mathbf{F}_{l},\mathbf{K}_{l},\mathbf{A},\mathcal{D}_{l})$,
where $\mathbf{F}_{l}=\\{\mathbf{f}_{1},\ldots,\mathbf{f}_{K}\\}^{\top}$ are
auxiliary variables for data in $\mathcal{D}_{l}$, and $\mathbf{K}_{l}$
denotes $\mathbf{K}(\mathbf{x}_{l},\mathbf{x}_{l})$. It follows that, for a
given class, the update for the weights across data in $\mathcal{D}_{l}$ is
given by,

$\hat{\mathbf{w}}_{k}=(\mathbf{K}_{l}\mathbf{K}_{l}^{\top}+\mathbf{A}_{k})^{-1}\mathbf{K}_{l}\mathbf{f}_{k}^{\top}$
(C.8)

where $\mathbf{A}_{k}$ is a diagonal matrix formed from the $k^{\text{th}}$
column of $\mathbf{A}$, i.e.
$\mathbf{A}_{k}=\text{diag}(\alpha_{1,k},\ldots,\alpha_{n,k})$.

Following the formulation in [12, 35], the posterior distribution over the
auxiliary variables is derived to be a product of conically-truncated Gaussian
distributions. For a given class $j$, the auxiliary variables can be updated
$\forall k\neq j$ as follows,

$\tilde{f}_{k,i}\leftarrow\hat{\mathbf{w}}_{k}^{\top}\mathbf{k}_{i}-\frac{\mathbb{E}_{p(u)}\bigl{[}\mathcal{N}(\hat{\mathbf{w}}_{k}^{\top}\mathbf{k}_{i}-\hat{\mathbf{w}}_{j}^{\top}\mathbf{k}_{i},1)\prod_{\kappa\neq
j,k}\Phi(u+\hat{\mathbf{w}}_{j}^{\top}\mathbf{k}_{i}-\hat{\mathbf{w}}_{\kappa}^{\top}\mathbf{k}_{i})\bigr{]}}{\mathbb{E}_{p(u)}\bigl{[}\Phi(u+\hat{\mathbf{w}}_{j}^{\top}\mathbf{k}_{i}-\hat{\mathbf{w}}_{k}^{\top}\mathbf{k}_{i})\prod_{\kappa\neq
j,k}\Phi(u+\hat{\mathbf{w}}_{j}^{\top}\mathbf{k}_{i}-\hat{\mathbf{w}}_{\kappa}^{\top}\mathbf{k}_{i})\bigr{]}}$
(C.9)

and for class $j$,

$\tilde{f}_{j,i}\leftarrow\hat{\mathbf{w}}_{j}^{\top}\mathbf{k}_{i}-\sum_{\kappa\neq
j}(\tilde{f}_{\kappa,i}-\hat{\mathbf{w}}_{\kappa}^{\top}\mathbf{k}_{i})$
(C.10)

where $\tilde{f}_{k,i}$ denotes the mean of the distribution over the latent
auxiliary variable corresponding to class $k$ for data point $i$ in
$\mathcal{D}_{l}$.

By continuing to follow [35], one can derive the update for the hyperpriors
$\alpha_{i,k}$ as the mean of the posterior distribution
$p(\mathbf{A}|\mathbf{W},\mathcal{D}_{l})$, given by,

$\tilde{\alpha}_{i,k}=\frac{2\tau+1}{w_{i,k}^{2}+2\nu}$ (C.11)

The iterative learning process requires the repeated application of the
updates given in equations (C.11), (C.8), (C.9) and (C.10), until a
convergence criterion is met. Suitable convergence criteria are provided in
[12].

As discussed in Section 4.2.1, RVMs utilise a subset of ‘relevant’ samples to
construct basis functions, with mRVM1 adopting a bottom-up approach and mRVM2
adopting a top-down approach.

For mRVM1, $\mathcal{A}$ is initialised as an empty set with samples
subsequently added or removed from $\mathcal{A}$ based upon their contribution
to the objective function given by the marginal log-likelihood
$\mathcal{L}(\mathbf{A})=\log
p(\mathbf{F}_{l}|\mathbf{K_{l}},\mathbf{A})=\log\int
p(\mathbf{F}_{l}|\mathbf{K}_{l},\mathbf{W})p(\mathbf{W}|\mathbf{A})d\mathbf{W}$.
In accordance with [12], to ensure $\mathcal{L}(\mathbf{A})$ is differentiable
such that the fast type-II Maximum-Likelihood approach [49], can be employed,
it is assumed that the scale parameter for each sample is common across
classes, i.e. $\forall k\in\\{1,\ldots,K\\},\>\alpha_{i,k}=\alpha_{i}$. By
following the procedure detailed in [49], one arrives at the following
marginal likelihood decomposition,

$\mathcal{L}(\mathbf{A})=\mathcal{L}(\mathbf{A}_{-i})+\ell(\alpha_{i})$ (C.12)

where,

$\mathcal{L}(\mathbf{A}_{-i})=\sum_{k=1}^{K}-\frac{1}{2}\biggl{[}N\log
2\pi+\log|\mathcal{C}_{-i}|+\mathbf{f}_{k}^{\top}\mathcal{C}_{-i}^{-1}\mathbf{f}_{k}\biggr{]}$
(C.13)

and,

$\ell(\alpha_{i})=\sum_{k=1}^{K}\frac{1}{2}\biggl{[}\log\alpha_{i}-\log(\alpha_{i}+s_{i})+\frac{q_{k,i}^{2}}{\alpha_{i}+s_{i}}\biggr{]}$
(C.14)

with
$\mathcal{C}_{-i}=\mathbf{I}_{n-1}+\mathbf{K}_{l,-i}^{\top}\mathbf{A}_{-i}^{-1}\mathbf{K}_{l,-i}$,
where $\mathbf{I}_{n-1}$ denotes the identity matrix of size $n-1$, and where
the subscript $-i$ is used to indicate matrices with entries corresponding to
the $i^{\text{th}}$ data point in $\mathcal{D}_{l}$ removed. The quantities
$s_{i}$ and $q_{k,i}$, introduced in equation (C.14) are given by [49],

$s_{i}=\mathbf{k}_{i}^{\top}\mathcal{C}_{-i}^{-1}\mathbf{k}_{i}\quad\text{and}\quad
q_{k,i}=\mathbf{k}_{i}^{\top}\mathcal{C}_{-i}^{-1}\mathbf{f}_{k}$ (C.15)

and can be interpreted as a ‘sparsity factor’ and ‘quality factor’,
respectively. The sparsity factor indicates how much the descriptive
information provided by the $i^{\text{th}}$ data point is already provided by
the existing samples. The quality factor provides a measure of the
$i^{\text{th}}$ sample’s ability to describe class $k$. Maximising
$\mathcal{L}(\mathbf{A})$ with respect to $\alpha_{i}$ by following the
procedure presented in [49, 12], one can quantify the contribution of data
point $i$ to the objective function as,

$\theta_{i}=\sum_{k=1}^{K}q_{k,i}^{2}-Ks_{i}$ (C.16)

For each iteration in the learning procedure, the contribution $\theta_{i}$ is
then used to construct $\mathcal{A}$ with samples satisfying $\theta_{i}>0$
included, and other samples excluded. The model parameter update in equation
(C.8) can be expressed utilising the sparse subset of data $\mathcal{A}$ as
follows:

$\hat{\mathbf{W}}_{\ast}=(\mathbf{K}_{\ast}\mathbf{K}_{\ast}^{\top}+\mathbf{A}_{\ast})^{-1}\mathbf{K}_{\ast}\tilde{\mathbf{F}}^{\top}$
(C.17)

where $\mathbf{K}_{\ast}$ denotes $\mathbf{K}(\mathbf{x}_{l},\mathcal{A})$ and
is $n^{\ast}\times n$, and $\mathbf{A}_{\ast}$ is $n^{\ast}\times n^{\ast}$.
Finally, the update given in equation (C.11) becomes,

$\alpha_{i}=\frac{Ks_{i}^{2}}{\theta_{i}}$ (C.18)

As discussed in Section 4.2.1, mRVM2 prunes samples from $\mathcal{A}$. This
is achieved by excluding samples with scales $\alpha_{i,k}$ sufficiently large
that $w_{i,k}$ is negligible. The $i^{\text{th}}$ data point can be considered
insignificant and can be removed from $\mathcal{A}$ when
$\alpha_{i,k}>10^{5}\;\forall k\in\\{1,\ldots,K\\}$.