# Variational Inference for Semiparametric Bayesian Novelty Detection in Large
Datasets

Luca Benedetti  Eric Boniardi  Leonardo Chiani
Jacopo Ghirri  Marta Mastropietro  Andrea Cappozzo
Politecnico di Milano
Francesco<EMAIL_ADDRESS>
Università Cattolica del Sacro Cuore, Milan

###### Abstract

After being trained on a fully-labeled training set, where the observations
are grouped into a certain number of known classes, novelty detection methods
aim to classify the instances of an unlabeled test set while allowing for the
presence of previously unseen classes. These models are valuable in many
areas, ranging from social network and food adulteration analyses to biology,
where an evolving population may be present. In this paper, we focus on a two-
stage Bayesian semiparametric novelty detector, also known as Brand, recently
introduced in the literature. Leveraging on a model-based mixture
representation, Brand allows clustering the test observations into known
training terms or a single novelty term. Furthermore, the novelty term is
modeled with a Dirichlet Process mixture model to flexibly capture any
departure from the known patterns. Brand was originally estimated using MCMC
schemes, which are prohibitively costly when applied to high-dimensional data.
To scale up Brand applicability to large datasets, we propose to resort to a
variational Bayes approach, providing an efficient algorithm for posterior
approximation. We demonstrate a significant gain in efficiency and excellent
classification performance with thorough simulation studies. Finally, to
showcase its applicability, we perform a novelty detection analysis using the
openly-available Statlog dataset, a large collection of satellite imaging
spectra, to search for novel soil types.

##  1 Introduction

Statistical methods for novelty detection are becoming increasingly popular in
recent literature. Similar to standard supervised classifiers, these models
are trained on a fully labeled dataset and subsequently employed to group
unlabeled data. In addition, novelty detectors allow for the discovery and
consequent classification of samples showcasing patterns not previously
observed in the learning phase. For example, novelty detection models have
been successfully employed in uncovering unknown adulterants in food
authenticity studies (Cappozzo et al., 2020; Fop et al., 2022), collective
anomalies in high energy particle physics (Vatanen et al., 2012), and novel
communities in social network analyses (Bouveyron, 2014).
Formally, we define a _novelty detector_ as a classifier trained on a set (the
training set) characterized by a specific number of classes that is used to
predict the labels of a second set (the test set). For a detailed account of
the topic, the interested reader is referred to the reviews of Markou and
Singh (2003b) and Markou and Singh (2003a), where the former is devoted
explicitly to statistical methods for novelty detection.

Within the statistical approach, different methodologies have been proposed to
construct novelty detectors. Recently, Denti et al. (2021) introduced a
Bayesian Robust Adaptive Novelty Detector (Brand). Brand is a semiparametric
classifier divided into two stages, focused on training and test datasets,
respectively. In the first stage, a robust estimator is applied to each class
of the labeled training set. By doing so, Brand recovers the representative
traits of each known class. The second step consists in fitting a
semiparametric nested mixture to the test set: the hierarchical structure of
the model specifies a convex combination between terms already observed in the
training set and a novelty term, where the latter is decomposed into a
potentially unbounded number of novel classes. At this point, Brand uses the
information extracted from the first phase to elicit reliable priors for the
known components.

In their article, Denti et al. (2021) devised an MCMC algorithm for posterior
estimation. The algorithm is based on a variation of the slice sampler (Kalli
et al., 2011) for nonparametric mixtures, which avoids the ex-ante
specification of the number of previously unseen components, reflecting the
expected ignorance about the structure of the novelty term. As a result, the
algorithm obtains good performance as it targets the true posterior without
resorting to any truncation. Nonetheless, as is often the case when full MCMC
inference is performed, it becomes excessively slow when applied to large
multidimensional datasets.
In this paper, we aim to scale up the applicability of Brand by adopting a
variational inference approach, vastly improving its computational efficiency.
Variational inference is an estimating technique that approximates a complex
probability distribution by resorting to optimization (Jordan et al., 1999;
Ormerod and Wand, 2010), which has received significant attention within the
statistical literature in the past decade (Blei et al., 2017). In the Bayesian
set-up, the distribution of interest is, of course, the posterior
distribution. Variational inference applied to Bayesian problems (also known
as variational Bayes, VB from now on) induces a paradigm shift in the
approximation of a complicated posterior: we switch from a simulation problem
(MCMC) to an optimization one. Broadly speaking, the VB approach starts with
the specification of a class of simple distributions called the variational
family. Then, within this specified family, one looks for the member that
minimizes a suitable distributional divergence (e.g., the Kullbak-Lielber
divergence) from the actual posterior. By dealing with a minimization problem
instead of a simulation one, we can considerably scale up the applicability of
Brand, obtaining results for datasets with thousands of observations measured
in high-dimension in a fraction of the time needed by MCMC techniques.
Variational techniques have showcased the potential to enhance the relevance
of the Bayesian framework even to large datasets, sparking interest in its
theoretical properties (e.g., Wang and Titterington, 2012; Nieman et al.,
2022; Ray and Szabó, 2022), and its applicability to, for instance, logistic
regression (Jaakkola and Jordan, 2000; Rigon, 2019), network analysis
(Aliverti and Russo, 2022), and, more in general, non-conjugate (Wang, 2012)
and advanced models (Zhang et al., 2019).

The present paper is structured as follows. In Section 2, we review Brand and
provide a summary of the variational Bayes approach. In Section 3, we discuss
the algorithm and the hyperparameters needed for the VB version of the
considered model. Section 4 reports extensive simulation studies that examine
the efficiency and robustness of our method. Then, in Section 5, we present an
application to the Statlog dataset, openly available from the UCI dataset
repository. In this context, novelty detection is used to discover novel types
of soil starting from satellite images. We remark that this analysis would
have been prohibitive with a classical MCMC approach. Lastly, Section 6
concludes the manuscript.

##  2 Background: Bayesian novelty detection and variational Bayes

This section introduces the two stages in which Brand articulates and briefly
presents the core concepts of the variational Bayes approach.

###  2.1 The Brand model

In what follows, we review the two-stage procedure proposed in Denti et al.
(2021). The first stage centers around a fully labeled learning set from which
we extract robust information to set up a Bayesian semiparametric model in the
second stage. More specifically, consider a labeled training dataset with $n$
observations grouped into $J$ classes. This paper will focus on classes
distributed as multivariate Gaussians, but one can readily extend the model
using different distributions. To this regard, we will write
$\mathcal{N}\left(\cdot\mid\boldsymbol{\Theta}\right)$ to indicate a
Multivariate Gaussian density with generic location and scale parameters
$\boldsymbol{\Theta}=\left(\boldsymbol{\mu},\Sigma\right)$. Within each of the
$J$ training classes, we separately apply the Minimum Regularized Covariance
Determinant (MRCD) estimator (Boudt et al., 2020) to retrieve robust estimates
for the associated mean vector and covariance matrix. In detail, the MRCD
searches for a subset of observations with fixed size whose regularized sample
covariance has the lowest possible determinant. MRCD estimates are then
defined as the multivariate location and regularized scatter based on the so-
identified subset. Compared to the popular Minimum Covariance Determinant
(Rousseeuw and Driessen, 1999; Hubert et al., 2018), MRCD takes advantage of
regularization with a pre-specified positive definite target matrix to ensure
the existence of the solution even when the data dimension exceeds the sample
size. This feature implies that the robustness properties of the original
procedure are preserved whilst being applicable also to “$p>n$” problems. Such
a characteristic is paramount in the context considered in the present paper,
where we specifically aim to scale up the applicability of Brand to high-
dimensional scenarios. We will denote the estimates obtained for class $j$,
$j=1,\ldots,J,$ with
$\bar{\Theta}^{MRCD}_{j}=\left(\hat{\boldsymbol{\mu}}^{MRCD}_{j},\hat{\boldsymbol{\Sigma}}^{MRCD}_{j}\right)$.
These quantities will be employed to elicit informative priors for the $J$
known components in the second stage. Specifically, in so doing, outliers and
label noise that might be present in the training set will not hamper the
Bayesian mixture model specification hereafter reported.

In the second stage, we estimate a semiparametric Bayesian classifier on a
test set with $M$ unlabeled observations. We want to build a novelty detector,
i.e., a model that can discern between “known” units - which follow a pattern
already observed in the training set - or “novelties”. At this point, the
likelihood for each observation is a simple two-group mixture between a
generic novelty component $f_{nov}$ and a density reflecting the behavior
found in the training set $f_{obs}$. Given stage 1, it is immediate to specify
$f_{obs}$ as a mixture of $J$ multivariate Gaussians, whose hyperpriors
elicitation will be guided by each of the robust estimates
$\bar{\Theta}^{MRCD}_{j}$. Therefore, we can now write the likelihood for the
generic observation $\boldsymbol{y}_{m}$, where $m=1,\ldots,M$, as a mixture
of $J+1$ distributions:

$\mathcal{L}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}^{obs},\boldsymbol{\pi}\right)=\pi_{0}f_{nov}+\sum_{j=1}^{J}\pi_{j}\mathcal{N}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}^{obs}_{j}\right).$
(1)

Estimating this model allows us to either allocate each of the $M$
observations into one of the previously observed $J$ Gaussian classes or flag
it as a novelty generated from an unknown distribution $f_{nov}$.

The specification of $f_{nov}$ is more delicate. To specify a flexible
distribution that would reflect our ignorance regarding the novelty term, we
employ a Dirichlet process mixture model with Gaussian kernels (Escobar and
West, 1995). It is a mixture model where the mixing distribution is sampled
from a Dirichlet process (Ferguson, 1973), characterized by concentration
parameter $\gamma$ and base measure $H$. In other words, we have that
$f_{nov}(\boldsymbol{y}_{m})=\sum_{k=1}^{\infty}\omega_{k}\mathcal{N}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}^{nov}_{k}\right)$,
where the atoms are sampled from the base measure, i.e.,
$\boldsymbol{\Theta}^{nov}\sim H$, and the mixture weights follow the so-
called stick-breaking construction (Sethuraman, 1994), where
$w_{k}=v_{k}\prod_{l<k}(1-v_{l})$ and $v_{l}\sim Beta(1,\gamma)$. To indicate
this process, we will write $\boldsymbol{\omega}\sim SB(\gamma)$. Thus, the
likelihood in (1) becomes, for $m=1,\ldots,M$,

$\mathcal{L}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta},\boldsymbol{\pi}\right)=\pi_{0}\left[\sum_{k=1}^{\infty}\omega_{k}\mathcal{N}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}^{nov}_{k}\right)\right]+\sum_{j=1}^{J}\pi_{j}\mathcal{N}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{j}\right).$

This nested-mixture expression of the likelihood highlights a two-fold
advantage. First, it is highly flexible, effectively capturing departures from
the known patterns and flagging them as novelties. Second, the mixture nature
of $f_{nov}$ allows to automatically cluster the novelties, capturing
potential patterns that may arise. Furthermore, clusters in the novelty terms
characterized by very small sizes can be interpreted as simple outliers.

Notice that the previous nested-mixture likelihood can be conveniently re-
expressed as:

$\displaystyle\mathcal{L}\left(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}^{obs},\boldsymbol{\pi}\right)$
$\displaystyle=\sum_{k=1}^{\infty}\tilde{\pi}_{k}\mathcal{N}\left(\boldsymbol{y}_{m}\mid\tilde{\boldsymbol{\Theta}}_{k}\right),$
(2)

where, for $k\in\mathbb{N}$, we define
$\tilde{\pi}_{k}={\pi_{k}}^{\mathbbm{1}{\\{0<k\leq
J\\}}}(\pi_{0}\omega_{k-J})^{\mathbbm{1}{\\{k\geq J\\}}}$ and
$\tilde{\boldsymbol{\Theta}}_{k}=(\boldsymbol{\Theta}_{k}^{obs})^{\mathbbm{1}{\\{0<k\leq
J\\}}}(\boldsymbol{\Theta}_{k}^{nov})^{\mathbbm{1}{\\{k\geq J\\}}}$. Note
that, without loss of generality, we regard the first $J$ mixture components
as the known components and all the remaining ones as novelties.

Finally, as customary in mixture models, to ease the computation, we augment
the model by introducing the auxiliary variables $\xi_{m}\in\mathbb{N}$, for
$m=1,\ldots,M$, where $\xi_{m}=l$ means that the $m$-th observation has been
assigned to the $l$-th component. Therefore, the model in (2) simplifies to

$\displaystyle\boldsymbol{y}_{m}\mid\xi_{m},\tilde{\boldsymbol{\Theta}},\boldsymbol{\pi}\sim$
$\displaystyle\mathcal{N}\left(\boldsymbol{y}_{m}\mid\tilde{\boldsymbol{\Theta}}_{\xi_{m}}\right),\quad\quad{\xi_{m}}\mid{\tilde{\boldsymbol{\pi}}\overset{iid}{\sim}\sum_{k=1}^{\infty}{\tilde{\pi}_{k}\delta_{k}(\cdot)}}.$
(3)

We complete our Bayesian model with the following prior specifications for the
weights and the atoms:

$\displaystyle\boldsymbol{\pi}$
$\displaystyle\sim{Dirichlet}(\alpha_{0},\alpha_{1},...,\alpha_{J}),\quad\quad\boldsymbol{\omega}\sim{SB}(\gamma),$
(4) $\displaystyle\boldsymbol{\Theta}_{k}^{obs}$
$\displaystyle\sim{\mathcal{NIW}}(\boldsymbol{\mu}^{{obs}}_{k},\nu^{{obs}}_{k},\lambda^{{obs}}_{k},\boldsymbol{\Psi}^{{obs}}_{k}),\quad\quad
k=1,\ldots,J,$ $\displaystyle\boldsymbol{\Theta}_{k}^{nov}$
$\displaystyle\sim{\mathcal{NIW}}({\boldsymbol{\mu}_{0}^{nov}},\nu_{0}^{{nov}},\lambda_{0}^{{nov}},\boldsymbol{\Psi}_{0}^{nov}),\quad\quad
k=J+1,\ldots,\infty,$

where $\mathcal{NIW}$ indicates a Normal Inverse-Wishart distribution. To ease
the notation, let
$\boldsymbol{\Theta}=\\{\boldsymbol{\Theta}_{k}^{obs}\\}_{k=1}^{J}\cup\\{\boldsymbol{\Theta}_{k}^{nov}\\}_{k=J+1}^{\infty}$,
$\boldsymbol{\varrho}_{k}=(\boldsymbol{\mu}^{{obs}}_{k},\nu^{{obs}}_{k},\lambda^{{obs}}_{k},\boldsymbol{\Psi}^{{obs}}_{k})$
and
$\boldsymbol{\varrho}_{0}=(\boldsymbol{\mu}^{nov}_{0},\nu^{nov}_{0},\lambda^{nov}_{0},\boldsymbol{\Psi}^{nov}_{0})$.
We remark that the values for the hyperparameters
$\\{\boldsymbol{\varrho}_{k}^{obs}\\}_{k=1}^{J}$ are defined according to the
robust estimates obtained in the first stage. See Denti et al. (2021) for more
details about the hyperparameters specification.

Given the model, it is easy to derive the joint law of the data and the
parameters, which is proportional to the posterior distribution of interest
$p\left({\boldsymbol{\Theta}},\boldsymbol{\pi},\boldsymbol{\xi}\mid\boldsymbol{y}\right)$.
Therefore, the posterior we target for inference is proportional to:

$\displaystyle
p\left({\boldsymbol{\Theta}},\boldsymbol{\pi},\boldsymbol{\xi}\mid\boldsymbol{y}\right)\propto$
$\displaystyle\prod_{m=1}^{M}\left[\prod_{k=1}^{J}\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{obs})^{\mathbbm{1}_{\xi_{m}=k}}\prod_{k=J+1}^{\infty}\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{nov})^{\mathbbm{1}_{\xi_{m}=k}}\right]\times$
$\displaystyle\prod_{k=1}^{J}\mathcal{\mathcal{NIW}}(\boldsymbol{\Theta}_{k}^{obs}\mid\boldsymbol{\varrho}_{k})\prod_{k=J+1}^{\infty}{\mathcal{NIW}}(\boldsymbol{\Theta}_{k}^{nov}\mid\boldsymbol{\varrho}_{0})\times$
$\displaystyle\prod_{m=1}^{M}\left[\prod_{k=1}^{J}(\pi_{k})^{\mathbbm{1}_{\xi_{m}=k}}\prod_{k=J+1}^{\infty}\pi_{0}\left(v_{k-J}\prod_{h=1}^{k-J-1}(1-v_{h})\right)^{\mathbbm{1}_{\xi_{m}=k}}\right]\times$
$\displaystyle\prod_{k=0}^{J}(\pi_{k})^{\alpha_{k}-1}\prod_{k=1}^{\infty}(1-v_{k})^{\gamma-1}.$
(5)

The left panel of Figure 1 contains a diagram that summarizes how Brand works.
In the following, we will devise a VB approach to approximate (5) in a timely
and efficient manner. The following subsection briefly outlines the general
strategy underlying a mean-field variational approach, while the thorough
derivation of the algorithm used to estimate Brand is deferred to Section 3.

Figure 1: Diagrams depicting how Brand (left panel) and variational inference
(right panel) work, respectively.

###  2.2 A short summary of mean-field variational Bayes

Working in a Bayesian setting, we ultimately seek to estimate the posterior
distribution to draw inference. Unfortunately, the expression in (5) is not
available in closed form, and therefore we need to rely on approximations.
MCMC algorithms, which simulate draws from (5), could be prohibitively slow
when applied to large datasets. To this extent, we resort to variational
inference to recasts the approximation of (5) into an optimization problem.
For notational simplicity, while we present the basic ideas of VB, let us
denote a generic target distribution with
$p(\boldsymbol{\theta}\mid\boldsymbol{y})\propto
p(\boldsymbol{\theta},\boldsymbol{y})$.

As in Blei et al. (2017), we focus on a mean-field variational family
$\mathcal{Q}$, a set containing distributions for the parameters of interest
$\boldsymbol{\theta}$ that are all mutually independent:
$\mathcal{Q}=\\{q_{\boldsymbol{\zeta}}(\boldsymbol{\theta}):q_{\boldsymbol{\zeta}}\boldsymbol{\theta}=\prod_{j}q_{\zeta_{j}}({\theta}_{j})\\}$.
The postulated independence dramatically simplifies the problem at hand.
Notice that each member of $\mathcal{Q}$ depends on a set of variational
parameters denoted by $\boldsymbol{\zeta}$.

We seek, among the members of this family, the candidate that provides the
best approximation of our posterior distribution
$p(\boldsymbol{\theta}\mid\boldsymbol{y})$. Herein, the goodness of the
approximation is quantified by the Kullback-Leibler (KL) divergence. Thus, we
aim to find the member of the variational family $\mathcal{Q}$ that minimizes
the KL divergence between the actual posterior distribution and the
variational approximation. The KL divergence ${D}_{KL}(\cdot\mid\mid\cdot)$
can be written as ${D}_{KL}(p(\boldsymbol{\theta}\mid\boldsymbol{y})\mid\mid
q_{\boldsymbol{\zeta}}(\boldsymbol{\theta})))=\mathbb{E}_{q}[\log{q_{\boldsymbol{\zeta}}(\boldsymbol{\theta})}]-\mathbb{E}_{q}[\log{p(\boldsymbol{\theta},\boldsymbol{y})}]+\log
p(\boldsymbol{y}).$ Unfortunately, we cannot easily compute the _evidence_
$p(\boldsymbol{y})$. However, the evidence does not depend on any variational
parameter and can be treated as fixed w.r.t. $\boldsymbol{\theta}$ during the
optimization process. We then re-formulate the problem into an equivalent one,
the maximization of the Evidence Lower Bound (ELBO), which is fully
computable:

$ELBO(q)=\mathbb{E}_{q}[\log{p(\boldsymbol{z},\boldsymbol{y})}]-\mathbb{E}_{q}[\log{q_{\boldsymbol{\zeta}}(\boldsymbol{z})}].$
(6)

Maximizing (6) is equivalent to minimize the aforementioned KL divergence.
To detect the optimal member $q^{\star}\in\mathcal{Q}$, we employ a widely
used algorithm called _Coordinate Ascent Variational Inference_ (CAVI). It
consists of a one-variable-at-a-time optimization procedure. Indeed,
exploiting the independence postulated by the mean-field property, one can
show that the optimal variational distribution for the parameter $\theta_{j}$
is given by:

$q^{\star}_{\zeta_{j}}(\theta_{j})\propto\exp\\{\mathbb{E}_{-j}\left[\log
p(\theta_{j},\boldsymbol{\theta}_{-j},\boldsymbol{y})\right]\\},$ (7)

where $\boldsymbol{\theta}_{-j}=\\{\theta_{l}\\}_{l\neq j}$ and the expected
value is taken w.r.t. the densities of $\boldsymbol{\theta}_{-j}$. The CAVI
algorithm iteratively computes (7) for every $j$ until the ELBO does not
register any significant improvement. The basic idea behind the CAVI algorithm
is depicted in the right-half of Figure 1.
In the next subsections, we will derive the CAVI updates for the variational
approximation of our model, along with the corresponding expression of the
ELBO. We call the resulting algorithm the variational Brand, or VBrand.

##  3 Variational Bayes for the Brand model

In this section, we tailor the generic variational inference algorithm to our
specific case. First, we write our variational approximation, highlighting the
dependence of each distribution on its specific variational parameters,
collected in
$\boldsymbol{\zeta}=\left(\boldsymbol{\eta},\boldsymbol{\varphi},\boldsymbol{a},\boldsymbol{b},\boldsymbol{\rho}^{obs},\boldsymbol{\rho}^{nov}\right)$.
The factorized form we adopt is:

$\begin{split}q_{\boldsymbol{\zeta}}(\boldsymbol{\xi},\boldsymbol{\pi},\boldsymbol{v},\boldsymbol{\Theta}^{obs},\boldsymbol{\Theta}^{nov})=&q_{\boldsymbol{\eta}}(\boldsymbol{\pi})\prod_{m=1}^{M}q_{\boldsymbol{\varphi}^{(m)}}(\xi^{(m)})\prod_{k=1}^{T-1}q_{a_{k},b_{k}}(v_{k})\times\\\
&\prod_{k=1}^{J}q_{\boldsymbol{\rho}_{k}^{obs}}(\boldsymbol{\Theta}^{obs}_{k})\prod_{k=J+1}^{J+T}q_{\boldsymbol{\rho}_{k}^{nov}}(\boldsymbol{\Theta}^{nov}_{k}).\end{split}$
(8)

In Equation (8), we truncated the stick-breaking representation of the
Dirichlet Process on the novelty term at a pre-specified threshold $T$, as
suggested in Blei and Jordan (2006). This implies that $q(v_{T}=1)=1$ and that
all the variational mixture weights indexed by $t>T$ are equal to zero.

Then, we can exploit a key property of VB. Note that all the full conditionals
of the Brand model have closed-form expressions (c.f.r. Section S1 of the
Supplementary Material) and belong to the exponential family. This feature
greatly simplifies the search for the variational solution. Indeed, it ensures
that the corresponding optimal variational distributions belong to the same
family of the corresponding full conditional, with properly updated
variational parameters. Therefore, we can already state that
$q_{\boldsymbol{\eta}}(\boldsymbol{\pi})$ is the density function of a
Dirichlet distribution, $q_{\boldsymbol{\varphi}^{(m)}}(\xi^{(m)})$ follows a
Categorical distribution, each $q_{a_{k},b_{k}}(v_{k})$ follows a Beta
distributions, and
$q_{\boldsymbol{\rho}_{k}^{obs}}(\boldsymbol{\Theta}^{obs}_{k})$ and
$q_{\boldsymbol{\rho}_{k}^{nov}}(\boldsymbol{\Theta}^{nov}_{k})$ are both
distributed according to a Normal Inverse-Wishart.

###  3.1 The CAVI parameters update

Once the parametric expressions for the members of the variational family are
obtained, we can derive the explicit formulas to optimize the parameters via
the CAVI algorithm. In this subsection, we state the updating rules that have
to be iteratively computed to fit VBrand. As we will observe, the set of
responsibilities $\\{\varphi_{k}^{(m)}\\}_{k=1}^{J+T}$, i.e., the variational
probabilities $\varphi^{(m)}_{k}=q(\xi^{(m)}=k)$, will play a central role in
all the steps. In detail, the CAVI steps are as follows:

  1. 1.

$q^{\star}_{\boldsymbol{\eta}}(\boldsymbol{\pi})$ is the density of a
$Dirichlet(\eta_{0},\eta_{1},\ldots,\eta_{J})$, where

$\eta_{0}=\alpha_{0}+\sum_{m=1}^{M}\sum_{l=J+1}^{J+T}\varphi_{l}^{(m)},\quad\;\eta_{j}=\alpha_{j}+\sum_{m=1}^{M}\varphi_{j}^{(m)},\;\text{
for }j=1,\ldots,J.\\\ $ (9)

Here, each hyperparameter linked to a known component is updated with the sum
of the responsibilities of the data belonging to the same specific known
component. Likewise, the variational novelty probability hyperparameter
$\eta_{0}$ contains the sum of the responsibilities of all data belonging to
all the novelty terms.

  2. 2.

Each $q^{\star}_{a_{k},b_{k}}(v_{k})$, for $k=1,\ldots,T-1$, is the density of
a $Beta(a_{k},b_{k})$. The update for these variational parameters is given
by:

$a_{k}=1+\sum_{m=1}^{M}\varphi_{J+k}^{(m)},\quad
b_{k}=\gamma+\sum_{l=k+1}^{T}\sum_{m=1}^{M}\varphi_{J+l}^{(m)}.$ (10)

Here, as expected from the stick-breaking construction, the first parameter of
the variational $Beta$ distribution is updated with the sum of the
probabilities of each point belonging to a specific novelty cluster $k$. At
the same time, the second parameter is updated with the sum of variational
probabilities of belonging to one of the next novelty components.

  3. 3.

Both $q^{\star}_{\boldsymbol{\rho}_{k}^{nov}}(\boldsymbol{\Theta}^{nov}_{k})$
and $q^{\star}_{\boldsymbol{\rho}_{k}^{obs}}(\boldsymbol{\Theta}^{obs}_{k})$
are Nomal Inverse-Wishart densities. Let us start with the updating rules of
the known components. Note that each variational parameter
$\boldsymbol{\rho}_{k}^{obs}$ is a shorthand for
$\left(\boldsymbol{m}^{obs}_{k},\ell^{obs}_{k},u^{obs}_{k},\boldsymbol{S}^{obs}_{k}\right)$.
These parameters have the same interpretation as the parameters of (4),
contained in $\boldsymbol{\varrho}_{k}$. So we have

$\displaystyle\boldsymbol{m}^{obs}_{k}$
$\displaystyle=\frac{1}{\lambda^{obs}_{k}+\sum_{m=1}^{M}\varphi_{k}^{(m)}}\left(\lambda^{obs}_{k}\boldsymbol{\mu}^{obs}_{k}+\sum_{m=1}^{M}\boldsymbol{y}_{m}\varphi_{k}^{(m)}\right),$
(11) $\displaystyle\ell_{k}$
$\displaystyle=\lambda_{k}^{obs}+\sum_{m=1}^{M}\varphi_{k}^{(m)},\quad
u^{obs}_{k}=\nu_{K}^{obs}+\sum_{m=1}^{M}\varphi_{k}^{(m)},$
$\displaystyle\boldsymbol{S}^{obs}_{k}$
$\displaystyle=\boldsymbol{\Psi}^{obs}_{k}+\sum_{m=1}^{M}\hat{\boldsymbol{\Sigma}}_{k}^{(m)}+\frac{\lambda^{obs}_{k}\sum_{m=1}^{M}\varphi_{k}^{(m)}}{\lambda^{obs}_{k}+\sum_{m=1}^{M}\varphi_{k}^{(m)}}(\overline{\boldsymbol{y}}_{k}-\boldsymbol{\mu}^{obs}_{k})^{T}(\overline{\boldsymbol{y}}_{k}-\boldsymbol{\mu}^{obs}_{k}),$

where we defined
$\hat{\boldsymbol{\Sigma}}_{k}^{(m)}=(\boldsymbol{y}_{m}-\overline{\boldsymbol{y}}_{k})(\boldsymbol{y}_{m}-\overline{\boldsymbol{y}}_{k})^{T}\varphi_{k}^{(m)}$
and
$\overline{\boldsymbol{y}}_{k}=\sum_{m=1}^{M}\boldsymbol{y}_{m}\varphi_{k}^{(m)}/\sum_{m=1}^{M}\varphi_{k}^{(m)}.$
The update for the parameters in $\boldsymbol{\rho}_{k}^{nov}$ follows the
same structure, with the hyperprior parameters in $\boldsymbol{\varrho}^{nov}$
carefully substituted to $\boldsymbol{\varrho}^{obs}$.

  4. 4.

Updating the responsibilities $\\{\varphi_{k}^{(m)}\\}_{k=1}^{J+T}$ for
$m=1,\ldots,M,$ is the most challenging step of the algorithm, given the
nested nature of the mixture in (2). We recall that, for a given $m$, the
distribution $q_{\boldsymbol{\varphi}^{(m)}}(\xi^{(m)})$ is categorical with
$J+T$ levels. Thus, we need to compute the values for the $J+T$ corresponding
probabilities. For the known classes $k=1,\ldots,J$, we have

$\log\varphi_{k}^{(m)}\propto\mathbb{E}[\log{\pi_{k}}]+\mathbb{E}[\log{\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{obs})}],$
(12)

while for the novelty terms $k=J+1,\ldots,J+T$, we have

$\displaystyle\log\varphi_{k}^{(m)}\propto$
$\displaystyle\mathbb{E}[\log{\pi_{0}}]+\mathbb{E}[\log{v_{k-J}}]+\sum_{l=1}^{k-J-1}\mathbb{E}[\log\left(1-v_{l}\right)]+$
(13)
$\displaystyle\mathbb{E}[\log{\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{nov})}].$

For the sake of conciseness, we report the explicit expression for all the
terms of (12) and (13) in Section S2 of the Supplementary Material.

This means that the probability of datum $\boldsymbol{y}_{m}$ to belong to
cluster $k$ depends on the likelihood of $\boldsymbol{y}_{m}$ under that same
cluster and on the overall relevance of $k^{th}$ cluster. Such relevance is
determined as the expected value of the relative component of the
${Dirichlet}$ distributed $\boldsymbol{\pi}$ and, for novelties, the stick-
breaking weight, here unrolled in its ${Beta}$-distributed components.

###  3.2 The expression of the ELBO for VBrand

In this subsection, we report the terms that need to be derived to obtain the
ELBO in Equation (6) for VBrand model. We start by computing the first term of
Equation (6), which takes the following form:

$\begin{split}\mathbb{E}[\log{p}]=&\sum_{m=1}^{M}\left(\sum_{k=1}^{J}f_{1}^{(m,k)}+\sum_{k=J+1}^{J+T}f_{2}^{(k)}\right)+\sum_{k=1}^{J}f_{3}^{(k)}+\sum_{k=J+1}^{J+T}f_{4}^{(k)}+\\\
&\sum_{m=1}^{M}\left(\sum_{k=1}^{J}f_{5}^{(m,k)}+\sum_{k=J+1}^{J+T}f_{6}^{(m,k)}\right)+\sum_{k=0}^{J}f_{7}^{(k)}+\sum_{l=1}^{T}f_{8}^{(l)}+const,\end{split}$

where the quantities $\\{f_{k}\\}_{k=1}^{8}$ have the following expressions
(note we have suppressed the superscripts to ease the notation, and that
$\boldsymbol{\psi}(\cdot)$ indicates the digamma function):

$\displaystyle f_{1}$
$\displaystyle=\varphi_{k}^{(m)}\mathbb{E}[\log{{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{obs})}],$
$\displaystyle\quad f_{2}$
$\displaystyle=\varphi_{k}^{(m)}\mathbb{E}[\log{{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{nov})}],$
$\displaystyle f_{3}$
$\displaystyle=\mathbb{E}[\log{{\mathcal{NIW}}(\boldsymbol{\Theta}_{k}^{obs}\mid\boldsymbol{\varrho}_{k})}],$
$\displaystyle\quad f_{4}$
$\displaystyle=\mathbb{E}[\log{{\mathcal{NIW}}(\boldsymbol{\Theta}_{k}^{nov}\mid\boldsymbol{\varrho}_{0}})].$

Moreover, we have
$f_{5}=\varphi_{k}^{(m)}\left(\boldsymbol{\psi}\left(\eta_{k}\right)-\boldsymbol{\psi}\left(\sum_{j=0}^{J}\eta_{j}\right)\right)$,

$\displaystyle
f_{6}=\varphi_{k}^{(m)}\left[\boldsymbol{\psi}(\eta_{0})-\boldsymbol{\psi}(\sum_{j=0}^{J}\eta_{j})+\boldsymbol{\psi}(a_{k-J}+b_{k-J})+\sum_{h=1}^{k-J-1}\boldsymbol{\psi}(b_{h})-\boldsymbol{\psi}(a_{h}+b_{h})\right],$

and, lastly,
$f_{7}=(\alpha_{k}-1)(\boldsymbol{\psi}(\eta_{k})-\boldsymbol{\psi}(\sum_{j=0}^{J}\eta_{j})),$
and
$f_{8}=(\gamma-1)(\boldsymbol{\psi}(b_{l})-\boldsymbol{\psi}(a_{l}+b_{l}))$.
The second term of Equation (6) can be written as

$\mathbb{E}[\log{q}]=\sum_{m=1}^{M}\left(h_{1}^{(m)}+h_{2}^{(m)}+\sum_{k=1}^{T}h_{3}^{(m,k)}+\sum_{k=1}^{J}h_{4}^{(m,k)}+\sum_{k=1}^{T}h_{5}^{(m,k)}\right),$

with
$h_{1}=\sum_{k=1}^{J+T}\varphi_{k}^{(m)}\ln{\varphi_{k}^{(m)}}-\ln{\sum_{k=0}^{J+T}\varphi_{k}^{(m)}}$,
$h_{4}=\mathbb{E}[\log{\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{obs}\mid\boldsymbol{\rho}_{k}^{obs})}]$,
$h_{5}^{(m,k)}=\mathbb{E}[\log{\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{nov}\mid\boldsymbol{\rho}_{k}^{nov})}]$,
and, finally,

$\displaystyle h_{2}=$
$\displaystyle\sum_{j=0}^{J}(\eta_{j}-1)\left(\boldsymbol{\psi}(\eta_{j})-\boldsymbol{\psi}(\sum_{j=0}^{J}\eta_{j})\right)+\ln{\Gamma\left(\sum_{j=0}^{J}{\eta_{j}}\right)}-\sum_{j=0}^{J}\ln{\Gamma(\eta_{j})},$
$\displaystyle h_{3}=$
$\displaystyle(a_{k}-1)(\boldsymbol{\psi}(a_{k})-\boldsymbol{\psi}(a_{k}+b_{k}))+$
$\displaystyle(b_{k}-1)(\boldsymbol{\psi}(b_{k})-\boldsymbol{\psi}(a_{k}+b_{k}))-\ln{\left(\frac{\Gamma(a_{k})\Gamma(b_{k})}{\Gamma(a_{k}+b_{k})}\right)}.$

Additional details are deferred to Section S2 of the Supplementary Material.

At each step of the algorithm, one needs to update the parameters according to
the rules (9)-(13) and to evaluate all the terms in $\\{f_{k}\\}_{k=1}^{8}$
and $\\{h_{k}\\}_{k=1}^{5}$. Albeit faster than MCMC, this articulates over
numerous steps. To further reduce the overall computing time, an R package
relying on an efficient C++ implementation has been implemented. The package
is openly available at the GitHub repository JacopoGhirri/VarBRAND. In the
same repository, the interested reader can find all the R scripts written to
run the simulation studies and real data application that we will discuss in
the following sections.

##  4 Simulation studies

In this section, we report the performance of our variational algorithm on a
range of simulated datasets. In particular, our simulation study is
articulated into three different experiments. Each experiment investigates a
different aspect of the model while altogether providing a multi-faceted
description of its performance.

The first experiment focuses on the scaling capabilities of our proposal; the
second compares the results and efficiency of VBrand with the original slice
sampler, while the last one assesses the sensitivity of the recovered
partition to the hyperprior specification. For the MCMC algorithm, we follow
the default hyperprior setting suggested in Denti et al. (2021). Moreover, we
run the slice sampler for 20,000 iterations, discarding the first 10,000 as
burn-in. As for VBrand, we use the same hyperprior specifications used for the
MCMC (unless otherwise stated), while we use $k$-means estimation to
initialize the novelty terms’ means. We set a threshold $\varepsilon=10^{-9}$
as stopping rule.

SS1 | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}^{*}$ | $C_{5}^{*}$ | $C_{6}^{*}$ | $C_{7}^{*}$
---|---|---|---|---|---|---|---
$\boldsymbol{\mu}_{k}$ | (-5,-5) | (-4,-4) | (4,4) | (0,0) | (5,-10) | (5,-10) | (-10,-10)
$\sigma^{2}$ | 1 | 2 | 2 | 1 | 1 | 1 | 0.1
$\rho$ | 0.9 | 0 | 0 | -0.75 | 0.9 | 0.9 | 0
$n_{Tr}$ | 300 | 300 | 300 | – | – | – | –
$n_{Te}$ | 200 | 200 | 250 | 90 | 100 | 100 | 60
SS2 | $C_{1}$ | $C_{2}$ | $C_{3}^{*}$ | $C_{4}^{*}$ | $C_{5}^{*}$ | – | –
$\boldsymbol{\mu}_{k}$ | (2,2) | (-2,-2) | (2,-2) | (-2,2) | (0,0) | – | –
$n_{Tr}$ | 500 | 500 | – | – | – | – | –
$n_{Te}$ | 200 | 200 | 200 | 200 | 200 | – | –
SS3 | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}^{*}$ | $C_{5}^{*}$ | $C_{6}^{*}$ | –
$\boldsymbol{\mu}_{k}$ | (-5,5) | (5,5) | (5,-5) | (0,0) | (0,0) | (-5,5) | –
$\sigma^{2}$ † | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 1.5 | –
$\rho$ | 0 | 0 | 0 | 0.8 | -0.8 | 0 | –
$n_{Tr}$ | 50 | 50 | 50 | – | – | – | –
$n_{Te}$ | 1950 | 1950 | 1950 | 2000 | 2000 | 100 | –
Table 1: Characteristics of the synthetic datasets used in the three
simulation studies (SS1$-$SS3). The components flagged with ∗ are novelties. †
the variances reported in this table refer to the low overlap scenario. For
the high overlap scenario, we consider $\sigma^{2\star}=6\sigma^{2}$.

###  4.1 Classification performance

We test the classification performance of our proposal by applying VBrand to a
sequence of increasingly complex variations of a synthetic dataset. We monitor
the computation time and different clustering metrics to provide a complete
picture of the overall performance. In particular, we compute the Adjusted
Rand Index (ARI, Hubert and Arabie (1985)), the Adjusted Mutual Information
(AMI, Vinh et al. (2009)), and the Fowlkes-Mallows Index (FMI, Fowlkes and
Mallows (1983)). While the first two metrics correct the effect of agreement
solely due to chance, the latter performs well also if noise is added to an
existing partition. Thereupon, the joint inspection of these three quantities
aim to provide a complete picture of the results.

The considered data generating process (DGP) for this experiment is based on a
mixture of 7 bivariate Normals. In detail, the first three components
represent the known classes, appearing in both training and test sets, while
the remaining $4$ are deemed to be novelties, present only in the test set.
The components are characterized by different mean vectors
$\boldsymbol{\mu}_{k}$, correlation matrices $(\sigma_{k}^{2},\rho_{k})$, and
cardinalities in the training ($n_{Tr}$) and test ($n_{Te}$). The main
attributes of this DGP are summarised in the first block, named SS1, of Table
1.

Starting from this basic mechanism, we subsequently increment the difficulty
of the classification tasks. Specifically, we increase both the data dimension
and the sample size as follows:

  * •

_Sample size_ : while keeping unaltered the mixture proportions, we consider
the sample sizes $\tilde{n}_{k}=q\cdot n_{k}$, with multiplicative factor
$q\in\\{0.5,1,2.5,5,10\\}$ in both the training and test sets;

  * •

_Data dimensionality_ : we augment the dimensionality $p$ of the problem by
considering $p\in\\{2,3,5,7,10\\}$. Each added dimension (above the second)
comprises independent realization from a standard Gaussian. Note that the
resulting datasets define a particularly challenging discrimination task: all
the information needed to distinguish the different components is contained in
the first two dimensions. In contrast, the remaining ones only display
overlapping noise.

For each combination of sample size and dataset dimension, we generate 50
replications of each simulated dataset and summarize the results by computing
the means and the standard errors of the chosen metrics.

Figure 2: Performance metrics and elapsed time obtained by the VBrand
algorithm stratified by number of variables $p$ and size scaling $n$. The dots
represent the averages obtained over 50 replicates of the simulated
experiment, while the vertical bars display the associated standard errors.

Results for this experiment are reported in Figure 2. We immediately notice
that the clustering performances deteriorate as the dimensionality of the
problem increases. This trend is expected, especially given the induced
overlap in the added dimensions. However, the classification abilities of our
method remain consistent and satisfactory across all metrics. Indeed, as we
can see from the three panels at the top, ARI, AMI, and FMI are all strictly
above 70% across all scenarios. This outcome indicates that not only are the
known classes correctly identified and clustered as such, but the flexible
nonparametric component effectively captures the novelty term. The computation
time (bottom panel) grows exponentially as a function of the test set
cardinality. Interestingly, the increment of data dimensionality does not
significantly impact the computational costs, suggesting an effective
scalability of our proposal to high dimensional problems. Indeed, even when
the test size is in the order of tens of thousands, the devised CAVI algorithm
always reaches convergence in less than half a minute. Lastly, we remark that
the time needed for convergence is sensitive to the different initialization
provided to the algorithm, explaining the high variance that characterized the
computational costs.

###  4.2 Comparison with MCMC

We now compare the variational and the MCMC approaches for approximating the
Brand posterior. The MCMC algorithm we consider is the modified slice-sampler
introduced in Denti et al. (2021). We compare the estimating approaches
leveraging on the same DGP highlighted in the previous section and slightly
modified accordingly. In detail, we consider five spherical bivariate Gaussian
to generate the classes, out of which three are deemed as novelties: the basic
details of the resulting DGP can be found in the second block of Table 1
(SS2). We consider 50 different scenarios resulting from the interactions of
the levels of following three attributes:

  * •

_Simple vs. complex scenarios_ : we set the variance of all the mixture
components to either $\sigma_{S}^{2}=0.2$ or $\sigma_{C}^{2}=0.75$ (the only
exception being $\sigma_{C_{5}^{*}}^{2}=0.375$). The former value implies
clear separation among the elements in the simple scenario. In contrast, the
latter variance defines the complex case, where we induce some overlap that
may hinder the classification. A descriptive plot is displayed in Section S3
of the Supplementary Material;

  * •

_Sample size_ : we modify the default sample size $n=1000$ by considering
different multiplicative factors in $q\in\\{0.5,1,2.5,5,10\\}$, thus obtaining
datasets ranging from $500$ to $10000$ observations.

  * •

_Data dimensionality_ : we augment the dimensionality $p$ of the problem by
considering $p\in\\{2,3,5,7,10\\}$. The dimensionality augmentation is carried
out as described in Section 4.1.

We assess the classification performances with the same metrics previously
introduced. For each of the 50 simulated scenarios, we perform 50 Monte Carlo
replicates to assess the variation in the performances. A summary of the
classification results under the complex scenario are reported in Figure 3.
The panels show that the slice sampler always outperforms the VB
implementation in the two-dimensional case. However, as the dimensionality of
the dataset increases the MCMC performance rapidly drops, while VBrand always
obtains good clustering recovery, irrespective of the data dimensionality and
the sample size. Similar results are obtained under the simple scenario, for
which a summarizing plot can be found in Section S3 of the Supplementary
Material.

Figure 3: Performance metrics (MCMC in red, variational Bayes - VB - in blue)
stratified by number of variables $p$ and sample size scaling factor $q$ under
the complex scenario. The dots represent the averages obtained over 50
replicates of the simulated experiment, while the vertical bars display the
associated standard errors.

Figure 4 compares the algorithms in terms of computation time. As expected,
VBrand provides results in just a fraction of the time required by the MCMC
approach, being approximately two orders of magnitude faster.

Figure 4: Computational time in seconds (MCMC in red, variational Bayes - VB -
in blue) grouped by number of variables $p$, sample size scaling factor $q$,
and type of scenario. The dots represent the averages obtained over 50
replicates of the simulated experiment, while the vertical bars display the
associated standard errors.

These results cast light not only on the apparent gain in computational speed
when using the variational approach, which is expected by any means, but also
on the superior recovering of the underlying partition in the test set in more
complex scenarios.

###  4.3 Sensitivity Analysis

Finally, we investigate the sensitivity of the model classification to
different hyperprior specifications. Two pairs of crucial hyperparameters may
considerably affect the clustering results. The first two are
$\boldsymbol{\alpha}$ and $\gamma$, which drive the parametric and
nonparametric mixture weights, respectively. The second pair is given by the
$\mathcal{NIW}$ precision $\lambda^{nov}_{0}$ and degrees of freedom
$\nu^{nov}_{0}$ of the novelty components. To assess their impact, we devise a
sensitivity analysis considering each possible combination of the following
hyperparameters: $\boldsymbol{\alpha}\in\\{0.1,0.55,1\\}$,
$\gamma\in\\{1,5.5,10\\}$, $\lambda^{nov}_{0}\in\\{1,5,10\\}$, and
$\nu^{nov}_{0}\in\\{4,52,100\\}$, thus defining 81 scenarios. We fit VBrand to
a dataset composed of five Normals in a bivariate domain, considering a fixed
sample size for both training and test sets. We chose a small sample size for
the training set to limit the informativeness of the robust estimation
procedure. Moreover, for each combination of the hyperparameters, we consider
both low and high values for the variances of the mixing components, obtaining
scenarios with low overlap (LOV) and high overlap (HOV), respectively.
Additional details about the data-generating process can be found in the third
block of Table 1 (SS3). For this experiment, we compare the retrieved
partitions in terms of ARI, as done in Sections 4.1 and 4.2, and by monitoring
the $F1$ score, i.e., the harmonic mean of precision and recall. The results
are reported in Figure 5. We immediately notice by inspecting the panels that
for the LOV case the method performs well regardless of the combination of
hyper-parameters chosen for the prior specification. In the HOV scenario the
recovery of the underlying true data partition is less effective, as it is
nevertheless expected. In particular, it seems that setting a high value for
the degrees of freedom $\nu^{nov}_{0}$ produces a slight drop in the ARI
metric. This behavior is due to the extra flexibility allowed to the novelty
component, by which some of the units belonging to the known groups are
incorrectly captured by the novelty term.

Figure 5: Classification metrics obtained over 50 replicates for 81
combinations of the hyperparameters $\boldsymbol{\alpha}$, $\gamma$,
$\lambda^{nov}_{0}$, and $\nu^{nov}_{0}$ under the low overlap (LOV) and high
overlap (HOV) cases. The dots represent the averages obtained over 50
replicates of the simulated experiment, while the vertical bars display the
associated standard errors.

In summary, our proposal increases the scalability of the approach introduced
in Denti et al. (2021), allowing for the original novelty detector to be
successfully applied to complex scenarios in a more timely manner, without
being it affected by the hyper-parameters specification.

##  5 Application to novel soil type detection

For our application, we consider the Statlog (Landsat Satellite) Data Set,
publicly obtained from the UCI machine-learning
repository222https://archive.ics.uci.edu/ml/datasets/Statlog+%28Landsat+Satellite%29.
It consists of a collection of observations from satellite images of different
soils. Each image contains four spectral measurements over a 3x3 grid,
recorded to classify the soil type captured in the picture. There are a total
of six different soil types recorded: Red Soil (RS), Grey Soil (GS), Damp Grey
Soil (DGS), Very Damp Grey Soil (VDGS), Cotton Crop (CC), and Soil with
Vegetation Stubble (SVS). We frame the original classification problem in a
novelty detection task by removing the images of CC and SVS from the training
set, leaving these groups in the test set to be detected as novelties.

Even when performing a simple classification, a method that can account for
the possible presence of previously unseen scenarios can be of paramount
utility in many fields. For example, new plants (Christenhusz and Byng, 2016),
animals (Camilo et al., 2011), and viruses (Woolhouse et al., 2012) are
progressively discovered every year. Likewise, related to our application,
landscapes present novelties at increasing rates (Finsinger et al., 2017).
Moreover, a scalable model that can discern and separate outlying observations
is necessary when dealing with real-world data, allowing the results to be
robust to outliers or otherwise irregular observations. Once these
observations are flagged, they can be the objective of future investigations.
Thus, our novelty detection application to the Statlog dataset is a nontrivial
example that may inspire future applications of our method.

The original data are already split into training and test sets. After
removing the CC and SVS classes from the training set, we obtain a training
set of 3486 observations. The test set instead contains 2000 instances. Each
observation includes the four spectral values recorded over the 9-pixel grids.
Therefore, we will model these data with a semiparametric mixture of
36-dimensional multivariate Normals. Given the large dimensionality of the
dataset, the application of the MCMC estimation approach is problematic in
terms of both required memory and computational time. Indeeed, estimating the
model via slice sampler becomes unfeasible for most laptop computers.
Moreover, we recall that the MCMC approach showed some numerical instabilities
in our simulation studies when applied to large dimensional datasets.

We apply VBrand adopting a mixture with full covariance matrices to capture
the potential dependence across the different pixels. Being primarily
interested in clustering, we first rescale the data to avoid numerical
problems. In detail, we divided all the values in both training and test by
4.5 to reduce the dispersion. Indeed, before the correction, the variabilities
within groups ranged from 25.40 to 228.04. After, the within-group variability
ranges from 1.25 to 11.26, significantly improving the stability of the
algorithm.

Since the variational techniques are likely to find a locally optimal
solution, we run the CAVI algorithm 200 times adopting different
initializations. For each run, we obtain different random starting point as
follows:

  * •

we set the centers for the novelty NIWs equal to the centers returned by a
$k$-means algorithm performed over the whole test set, with $k$ being equal to
the chosen truncation;

  * •

the Dirichlet parameters, $\boldsymbol{\nu}^{nov}$ and
$\boldsymbol{\lambda}^{nov}$ are randomly selected. In particular, we sample
the Dirichlet parameters from $\left(0.1,1\right)$ and the $\mathcal{NIW}$
parameters from $\left(1,10\right)$ through Latin Hypercube sampling;

The other variational hyperparameters are assumed fixed, equal to the
corresponding prior hyperparameters. In Section S4 of the Supplementary
Material, we report a detailed list of the hyperparameter specifications that
we adopted for this analysis.

As a final result, we select the run with the highest ELBO value, being the
one whose variational distribution has the lowest KL divergence from the true
posterior. The top-left panel of Figure 6 shows the ELBO trends for all the
run we performed.

Figure 6: Top-left panel: collection of ELBO trajectories obtained via CAVI
updates starting from 200 different random intializations; the trajectory
providing the highest ELBO is highlighted in blue (the $y$ axis is truncated
for improved visualization). Top-right panel: projection of the test dataset
onto a two-dimensional space via the tSNE algorithm. Bottom panel: heatmap of
the resulting confusion matrix.

The bottom panel of Figure 6 reports the resulting confusion matrix: we
observe that the algorithm successfully detected both novelties, achieving a
satisfactory classification performance of the previously observed soil types.
However, the model struggles with classifying the DGS instances, often
mistaken for GS or VDGS. Such difficulty is explained by the overlap between
these groups, as shown by the visualization of the test set obtained via the
tSNE projection (Hinton and van der Maaten, 2008), reported in the top-right
panel of Figure 6: from the plot, we see that it is not straightforward to
establish clear boundaries between GS, DGS, and VDGS soil types.

Overall, VBrand captures the main traits of the data and flags some
observations as outliers (e.g., Novelty clusters 4 and 5), which may warrant
further investigation. All in all, our variational approach provides a good
clustering solution in a few seconds and it is fast enough to allow for a
brute-force search for a better, albeit only locally optimal, solution
employing multiple initializations.

##  6 Discussion and conclusions

In this paper we introduced VBrand, a variational Bayes algorithm for novelty
detection, to classify instances of a test set that may conceal classes not
observed in a training set. We showed how VBrand outperforms the previously
proposed slice sampler implementation in terms of both computational time and
robustness of the estimates. The application to soil data provides an example
of the versatility of our method in a context where the MCMC algorithm fails
because of the large dimensionality of the problem.
Our results pave the way for many possible extensions. First, given recent
developments in Bayesian nonparametric literature, the variational algorithm
can be enriched by adding a hyperprior distribution for the concentration
parameter of the novelty DP (Ascolani et al., 2022). While in practice VBrand
already obtains very good classification performance, this addition would lead
to a consistent model for the number of true clusters. Second, we can consider
different likelihood specifications and develop variational inference novelty
detectors for, but not limited to, functional or graphical data. Third, at the
expense of efficiency, we can explore more complex specifications for the
variational distributions, as in structured variational inference. Albeit
potentially slower, this choice would lead to an algorithm that could better
capture the complex structure of the posterior distribution we are targeting.
Finally, we can resort to stochastic variational inference algorithms (Hoffman
et al., 2003), to scale up VBrand’s applicability to massive datasets that
could benefit from novelty detection techniques.

##  7 Acknowledgments

The authors are grateful to Fabio Lavezzo for his help in the early-stage
research that led to the development of this article.

## References

  * Aliverti and Russo [2022] Emanule Aliverti and Massimiliano Russo.  Stratified Stochastic Variational Inference for High-Dimensional Network Factor Model.  _Journal of Computational and Graphical Statistics_ , 31(2):502–511, 2022.  ISSN 15372715.  10.1080/10618600.2021.1984929.
  * Ascolani et al. [2022] F Ascolani, A Lijoi, G Rebaudo, and G Zanella.  Clustering consistency with Dirichlet process mixtures.  _Biometrika_ , 2022.  ISSN 0006-3444.  10.1093/biomet/asac051.
  * Bishop [2006] Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer New York, NY, 2006.  ISBN 978-0-387-31073-2.
  * Blei and Jordan [2006] David M Blei and Michael I Jordan.  Variational Inference for Dirichlet Process Mixtures.  _Bayesian Analysis_ , 1(1):121–144, 2006.  URL http://www.cs.berkeley.edu/$∼$blei/.
  * Blei et al. [2017] David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe.  Variational Inference: A Review for Statisticians.  _Journal of the American Statistical Association_ , 112(518):859–877, 2017.  ISSN 1537274X.  10.1080/01621459.2017.1285773.  URL https://doi.org/10.1080/01621459.2017.1285773.
  * Boudt et al. [2020] Kris Boudt, Peter J. Rousseeuw, Steven Vanduffel, and Tim Verdonck.  The minimum regularized covariance determinant estimator.  _Statistics and Computing_ , 30(1):113–128, February 2020.  ISSN 0960-3174.  10.1007/s11222-019-09869-x.  URL http://link.springer.com/10.1007/s11222-019-09869-x.
  * Bouveyron [2014] Charles Bouveyron.  Adaptive mixture discriminant analysis for supervised learning with unobserved classes.  _Journal of Classification_ , 31(1):49–84, 2014\.  ISSN 14321343.  10.1007/s00357-014-9147-x.  URL https://link.springer.com/content/pdf/10.1007/s00357-014-9147-x.pdf.
  * Camilo et al. [2011] Mora Camilo, P. Tittensor Derek, Adl Sina, G. B. Simpson Alastair, and Worm Boris.  How Many Species Are There on Earth and in the Ocean?  _PLoS Biology_ , 9(8):e1001127, 2011.
  * Cappozzo et al. [2020] Andrea Cappozzo, Francesca Greselin, and Thomas Brendan Murphy.  Anomaly and Novelty detection for robust semi-supervised learning.  _Statistics and Computing_ , 30(5):1545–1571, September 2020.  ISSN 0960-3174.  10.1007/s11222-020-09959-1.  URL http://arxiv.org/abs/1911.08381http://link.springer.com/10.1007/s11222-020-09959-1.
  * Christenhusz and Byng [2016] Maarten J.M. Christenhusz and James W. Byng.  The number of known plants species in the world and its annual increase.  _Phytotaxa_ , 261(3):201–217, 2016.  ISSN 11793163.  10.11646/phytotaxa.261.3.1.
  * Denti et al. [2021] Francesco Denti, Andrea Cappozzo, and Francesca Greselin.  A two-stage bayesian semiparametric model for novelty detection with robust prior information.  _Statistics and Computing_ , 31(4):1–19, 2021\.
  * Escobar and West [1995] Michael D. Escobar and Mike West.  Bayesian Density Estimation and Inference Using Mixtures.  _Journal of the American Statistical Association_ , 90(430):577–588, June 1995.  ISSN 0162-1459.  10.1080/01621459.1995.10476550.  URL http://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476550.
  * Ferguson [1973] Thomas S. Ferguson.  A Bayesian Analysis of Some Nonparametric Problems.  _The Annals of Statistics_ , 1(2):209–230, 1973\.  ISSN 0090-5364.  10.1214/aos/1176342360.
  * Finsinger et al. [2017] Walter Finsinger, Thomas Giesecke, Simon Brewer, and Michelle Leydet.  Emergence patterns of novelty in European vegetation assemblages over the past 15 000 years.  _Ecology Letters_ , 20(3):336–346, 2017.  ISSN 14610248.  10.1111/ele.12731.
  * Fop et al. [2022] Michael Fop, Pierre-Alexandre Mattei, Charles Bouveyron, and Thomas Brendan Murphy.  Unobserved classes and extra variables in high-dimensional discriminant analysis.  _Advances in Data Analysis and Classification_ , 16(1):55–92, mar 2022.  ISSN 1862-5347.  10.1007/s11634-021-00474-3.  URL https://doi.org/10.1007/s11634-021-00474-3https://link.springer.com/10.1007/s11634-021-00474-3.
  * Fowlkes and Mallows [1983] E. B. Fowlkes and C. L. Mallows.  A method for comparing two hierarchical clusterings.  _Journal of the American Statistical Association_ , 78(383):553–569, 1983.  ISSN 1537274X.  10.1080/01621459.1983.10478008.
  * Hinton and van der Maaten [2008] Geoffrey Hinton and Laurens van der Maaten.  Visualizing Data using t-SNE.  _Journal of Machine Learning Research_ , 9:2579–2605, 2008\.  URL http://jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf%0Ahttp://www.jmlr.org/papers/v9/vandermaaten08a.html%5Cnfile:///Files/63/63E4B948-D809-4073-8CE0-E56194C96FD8.pdf.
  * Hoffman et al. [2003] Matt Hoffman, Chong Wang, and John Paisley.  Stochastic Variational Inference.  _Journal ofMachine Learning Research_ , pages 1–52, 2003.
  * Hubert and Arabie [1985] Lawrence Hubert and Phipps Arabie.  Comparing partitions.  _Journal of Classification_ , 2(1):193–218, 1985\.  ISSN 01764268.  10.1007/BF01908075.
  * Hubert et al. [2018] Mia Hubert, Michiel Debruyne, and Peter J. Rousseeuw.  Minimum covariance determinant and extensions.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 10(3):1–11, 2018.  ISSN 19390068.  10.1002/wics.1421.
  * Jaakkola and Jordan [2000] Tommi S. Jaakkola and Michael I. Jordan.  Bayesian parameter estimation via variational methods.  _Statistics and Computing_ , 10(1):25–37, 2000\.  ISSN 09603174.  10.1023/A:1008932416310.
  * Jordan et al. [1999] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul.  Introduction to variational methods for graphical models.  _Machine Learning_ , 37(2):183–233, 1999.  ISSN 08856125.  10.1023/A:1007665907178.
  * Kalli et al. [2011] Maria Kalli, Jim E Griffin, and Stephen G Walker.  Slice sampling mixture models.  _Statistics and Computing_ , 21(1):93–105, 2011\.  ISSN 09603174.  10.1007/s11222-009-9150-y.  URL https://link-springer-com.proxy.unimib.it/content/pdf/10.1007{%}2Fs11222-009-9150-y.pdf.
  * Markou and Singh [2003a] Markos Markou and Sameer Singh.  Novelty detection: a review—part 2:.  _Signal Processing_ , 83(12):2499–2521, December 2003a.  ISSN 01651684.  10.1016/j.sigpro.2003.07.019.  URL https://linkinghub.elsevier.com/retrieve/pii/S0165168403002032.
  * Markou and Singh [2003b] Markos Markou and Sameer Singh.  Novelty detection: a review—part 1: statistical approaches.  _Signal Processing_ , 83(12):2481–2497, December 2003b.  ISSN 01651684.  10.1016/j.sigpro.2003.07.018.  URL https://linkinghub.elsevier.com/retrieve/pii/S0165168403002020.
  * Nieman et al. [2022] Dennis Nieman, Botond Szabo, and Harry van Zanten.  Contraction rates for sparse variational approximations in Gaussian process regression.  _Journal of Machine Learning Research_ , 23:1–26, 2022\.  URL http://arxiv.org/abs/2109.10755.
  * Ormerod and Wand [2010] J. T. Ormerod and M. P. Wand.  Explaining variational approximations.  _American Statistician_ , 64(2):140–153, 2010\.  ISSN 00031305.  10.1198/tast.2010.09058.
  * Ray and Szabó [2022] Kolyan Ray and Botond Szabó.  Variational Bayes for High-Dimensional Linear Regression With Sparse Priors.  _Journal of the American Statistical Association_ , 117(539):1270–1281, 2022.  ISSN 1537274X.  10.1080/01621459.2020.1847121.
  * Rigon [2019] Tommaso Rigon.  An enriched mixture model for functional clustering.  _arXiv_ , 2019.  ISSN 23318422.  URL http://arxiv.org/abs/1907.02493.
  * Rousseeuw and Driessen [1999] Peter J. Rousseeuw and Katrien Van Driessen.  A fast algorithm for the minimum covariance determinant estimator.  _Technometrics_ , 41(3):212–223, August 1999\.  ISSN 0040-1706.  10.1080/00401706.1999.10485670.  URL http://www.tandfonline.com/doi/abs/10.1080/00401706.1999.10485670.
  * Sethuraman [1994] Jayaram Sethuraman.  A constructive definition of Dirichlet Process prior.  _Statistica Sinica_ , 4(2):639–650, 1994.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24305538.
  * Vatanen et al. [2012] Tommi Vatanen, Mikael Kuusela, Eric Malmi, Tapani Raiko, Timo Aaltonen, and Yoshikazu Nagai.  Semi-supervised detection of collective anomalies with an application in high energy particle physics.  In _Proceedings of the International Joint Conference on Neural Networks_ , 2012.  ISBN 9781467314909.  10.1109/IJCNN.2012.6252712.  URL http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.223.215{&}rep=rep1{&}type=pdf.
  * Vinh et al. [2009] Nguyen Xuan Vinh, Julien Epps, and James Bailey.  Information theoretic measures for clusterings comparison: Is a correction for chance necessary?  _ACM International Conference Proceeding Series_ , 382:1–8, 2009.  10.1145/1553374.1553511.
  * Wang and Titterington [2012] Bo Wang and D. Titterington.  Convergence and asymptotic normality of variational Bayesian approximations for exponential family models with missing values.  _ArXiv preprint_ , 2012.  URL http://arxiv.org/abs/1207.4159.
  * Wang [2012] Chong Wang.  Variational Inference in Nonconjugate Models.  _Journal of Machine Learning Research_ , 2012.
  * Woolhouse et al. [2012] Mark Woolhouse, Fiona Scott, Zoe Hudson, Richard Howey, and Margo Chase-Topping.  Human viruses: Discovery and emeraence.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ , 367(1604):2864–2871, 2012.  ISSN 14712970.  10.1098/rstb.2011.0354.
  * Zhang et al. [2019] Cheng Zhang, Judith Butepage, Hedvig Kjellstrom, and Stephan Mandt.  Advances in Variational Inference.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 41(8):2008–2026, 2019.  ISSN 19393539.  10.1109/TPAMI.2018.2889774.

## Supplementary Material

This Supplementary Material contains the formulas to implement the full
conditionals of the Brand Gibbs sampler in Section S1. In Section S2, we
report additional details that we used to derive the updating rules of the
CAVI algorithm and the value of the ELBO. Section S3 contains additional
figures regarding the simulation studies presented in the main paper. Finally,
we report the hyperparameter specifications we adopted for the application to
the Statlog dataset in Section S4.

##  S1 Full conditionals for Brand with multivariate Gaussian components

Let $X\mid\cdots$ denote the full conditional distribution for the random
variable $X$. Then, the full conditionals for the parameters of the model
introduced in Denti et al. [2021] are as follows:

  * •

The memberships labels are updated according to

$\mathbb{P}\left[\xi_{m}=k\mid\cdots\right]=\begin{cases}\pi_{k}\cdot\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{obs})\cdot\mathbbm{1}_{\xi_{m}=k}&\quad\;k=1,\;...\;,J,\\\
\pi_{0}\cdot
v_{k}\left(\prod_{l=1}^{k-J-1}(1-v_{l})\right)\cdot\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{nov})\cdot\mathbbm{1}_{\xi_{m}=k},&\quad
k\geq J+1.\\\ \end{cases}$

  * •

The weights of the parametric mixtures and the stick-breaking variables are
sampled from

$\boldsymbol{\pi}\mid\cdots\overset{}{\sim}Dirichlet\left(\\{\alpha_{k}+n_{k}\\}_{k=0}^{J}\right),\quad\quad
v_{k}\mid\cdots\overset{}{\sim}Beta\left(n_{J+k}+1,\gamma+\sum_{l=k+1}^{\infty}n_{J+l}\right),$

  * •

The mixture locations are updated according to
$\boldsymbol{\Theta}_{k}^{q}\mid\cdots\overset{}{\sim}\mathcal{\mathcal{missing}}{NIW}(\boldsymbol{\mu}_{n},\lambda_{n},\nu_{n},\boldsymbol{\Psi}_{n})$,
with $q\in\\{obs,nov\\}$, where

$\displaystyle n_{k}=\begin{cases}\\#[\xi_{m}=k]&\text{if $k$ $\geq$ 1 },\\\
\\#[\xi_{m}>J]&\text{if $k$ = 0 },\\\
\end{cases}\quad\quad\boldsymbol{\mu}_{n}=\frac{\lambda_{0}\cdot\boldsymbol{\mu}_{0}+\sum_{m=1}^{M}\boldsymbol{y}_{m}\cdot\mathbbm{1}_{\xi_{m}=k}}{\lambda_{0}+n_{k}},$
$\displaystyle\lambda_{n}=\lambda_{0}+n_{k},\quad\quad\nu_{n}=\nu_{0}+n_{k},$
$\displaystyle\boldsymbol{\Psi}_{n}=\boldsymbol{\Psi}_{0}+\boldsymbol{S}+\frac{\lambda_{0}\cdot
n_{k}}{\lambda_{0}+n_{k}}\left(\frac{\sum_{m=1}^{M}\boldsymbol{y}_{m}\cdot\mathbbm{1}_{\xi_{m}=k}}{n_{k}}-\boldsymbol{\mu}_{0}\right)\left(\frac{\sum_{m=1}^{M}\boldsymbol{y}_{m}\cdot\mathbbm{1}_{\xi_{m}=k}}{n_{k}}-\boldsymbol{\mu}_{0}\right)^{T},$
$\displaystyle\boldsymbol{S}=\sum_{m=1}^{M}(\boldsymbol{y}_{m}-\overline{\boldsymbol{y}}^{k})\cdot(\boldsymbol{y}_{m}-\overline{\boldsymbol{y}}^{k})^{T}\cdot\mathbbm{1}_{\xi_{m}=k},\quad\quad\overline{\boldsymbol{y}}_{k}=\frac{\sum_{m=1}^{M}\boldsymbol{y}_{m}\cdot\mathbbm{1}_{\xi_{m}=k}}{\sum_{m=1}^{M}\mathbbm{1}_{\xi_{m}=k}}.$

Here, ($\boldsymbol{\mu}_{0}$, $\lambda_{0}$, $\nu_{0}$,
$\boldsymbol{\Psi}_{0}$) indicate the hyperparameters of the relative ${NIW}$
distribution, according to which atom $\boldsymbol{\Theta}_{k}^{q}$ is being
considered.

##  S2 Additional details regarding the VBrand algorithm

###  S2.1 Quantities useful to devise the CAVI updating rules

Here, we report some additional formulas to clarify the derivation of the CAVI
algorithm for VBrand. In our implementation, we used the following well-known
expected values:

  * •

$\mathbb{E}[\log{\pi_{k}}]=\boldsymbol{\psi}(\eta_{k})-\boldsymbol{\psi}(\overline{\eta}),$
with $\overline{\eta}=\sum_{k}\eta_{k}$

  * •

$\mathbb{E}[\log{v_{k-J}}]=\boldsymbol{\psi}(a_{k-J})-\boldsymbol{\psi}(a_{k-J}+b_{k-J})$
and
$\mathbb{E}[\log{(1-v_{l})}]=\boldsymbol{\psi}(b_{l})-\boldsymbol{\psi}(a_{l}+b_{l})$,

  * •

If $\Sigma^{-1}\sim\textit{Wishart}$, we have that:

$\mathbb{E}[\log{\mid\Sigma^{-1}\mid}]=-p\log{2}-\log{\mid\boldsymbol{\Psi}_{k}^{-1}\mid}-\sum_{l=1}^{p}\boldsymbol{\psi}\left(\frac{\nu_{k}-l+1}{2}\right).$

Given the previous results and the updating rules derived in Section S1, we
can easily derive the value for the negative entropy of a multivariate Normal
distribution for $q\in\\{obs,nov\\}$:

$\begin{split}\mathbb{E}[\log{\mathcal{N}(\boldsymbol{y}_{m}\mid\boldsymbol{\Theta}_{k}^{q})}]&=\frac{1}{2}\left(-p\log{2\pi}+\sum_{l=1}^{p}\boldsymbol{\psi}\left(\frac{\nu_{k}-l+1}{2}\right)+p\log{2}+\log{\mid\boldsymbol{\Psi}_{k}^{-1}\mid}\right)+\\\
&-\frac{1}{2}\left(\frac{p}{\lambda_{k}}+\nu_{k}(\boldsymbol{y}_{m}-\boldsymbol{\mu}_{k})^{T}\boldsymbol{\Psi}_{k}^{-1}(\boldsymbol{y}_{m}-\boldsymbol{\mu}_{k})\right).\end{split}$

###  S2.2 Quantities useful to derive the Elbo components

In the following formulas, $Tr(\cdot)$ denotes the trace of a matrix.
Moreover, the expression of the function $B(\cdot,\cdot)$ can be found in
Bishop [2006], Appendix B, p.693.

$\begin{split}\mathbb{E}[\log(\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{obs}\mid\boldsymbol{\varrho}_{k}))]=&\frac{1}{2}\big{(}p\log{\frac{\lambda_{k}^{obs}}{2\pi}}+\sum_{i=1}^{p}\psi(\frac{u_{k}^{obs}-i+1}{2})+\log{\mid\boldsymbol{S}_{k}^{obs^{-1}}\mid}+\\\
&-p\frac{\lambda_{k}^{obs}}{l_{k}^{obs}}-\lambda_{k}^{obs}u_{k}^{obs}(\boldsymbol{m}_{k}^{obs}-\boldsymbol{\mu}_{k}^{obs})^{T}\boldsymbol{S}_{k}^{obs^{-1}}(\boldsymbol{m}_{k}^{obs}-\boldsymbol{\mu}_{k}^{obs})\big{)}+\\\
&+\log{B(\boldsymbol{\Psi}_{k}^{obs^{-1}},\nu_{k}^{obs})}-\frac{1}{2}u_{k}^{obs}Tr(\boldsymbol{\Psi}_{k}^{obs}\boldsymbol{S}_{k}^{obs^{-1}})+const,\end{split}$

$\begin{split}\mathbb{E}[\log(\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{nov}\mid\boldsymbol{\varrho}_{0}))]=&\frac{1}{2}\big{(}p\log{\frac{\lambda_{0}^{nov}}{2\pi}}+\sum_{i=1}^{p}\psi(\frac{u_{k}^{nov}-i+1}{2})+\log{\mid\boldsymbol{S}_{k}^{nov^{-1}}\mid}+\\\
&-p\frac{\lambda_{0}^{nov}}{l_{k}^{nov}}-\lambda_{0}^{nov}u_{k}^{nov}(\boldsymbol{m}_{k}^{nov}-\boldsymbol{\mu}_{0}^{nov})^{T}\boldsymbol{S}_{k}^{nov^{-1}}(\boldsymbol{m}_{k}^{nov}-\boldsymbol{\mu}_{0}^{nov})\big{)}+\\\
&+\log{B(\boldsymbol{\Psi}_{0}^{nov^{-1}},\nu_{0}^{nov})}-\frac{1}{2}u_{0}^{nov}Tr(\boldsymbol{\Psi}_{0}^{nov}\boldsymbol{S}_{k}^{nov^{-1}})+const,\end{split}$

$\displaystyle\mathbb{E}[\log(\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{obs}\mid\boldsymbol{\rho}_{k}^{obs}))]\
=\frac{u_{k}^{obs}}{2}\log{\mid\boldsymbol{S}_{k}^{obs}\mid}-\frac{u_{k}^{obs}p}{2}\log{2}-\log\Gamma_{p}\left(\frac{u_{k}^{obs}}{2}\right)+\frac{p}{2}\log{l_{k}^{obs}}+$
$\displaystyle+\frac{u_{k}^{obs}+p+2}{2}\left(\log{\mid\boldsymbol{S}_{k}^{obs^{-1}}\mid}+\sum_{i=1}^{p}\boldsymbol{\psi}\left(\frac{u_{k}^{obs}-i+1}{2}\right)\right)-\frac{u_{k}^{obs}}{2}p,$

$\displaystyle\mathbb{E}[\log(\mathcal{NIW}(\boldsymbol{\Theta}_{k}^{nov}\mid\boldsymbol{\rho}_{k}^{nov}))]\
=\frac{u_{k}^{nov}}{2}\log{\mid\boldsymbol{S}_{k}^{nov}\mid}-\frac{u_{k}^{nov}p}{2}\log{2}-\log\Gamma_{p}\left(\frac{u_{k}^{nov}}{2}\right)+\frac{p}{2}\log{l_{k}^{nov}}+$
$\displaystyle+\frac{u_{k}^{nov}+p+2}{2}\left(\log{\mid\boldsymbol{S}_{k}^{nov^{-1}}\mid}+\sum_{i=1}^{p}\boldsymbol{\psi}\left(\frac{u_{k}^{nov}-i+1}{2}\right)\right)-\frac{u_{k}^{nov}}{2}p.$

##  S3 Additional figures

Figure 7: Example of the first two dimensions of the generated data under the
simulation study described in Section 4.2 of the main paper. Figure 8: Results
of the simulation study described in Section 4.2 of the main paper under the
simple scenario.

##  S4 Hyperparameter initializations for the novel soil detection analysis

We report a detailed list of the hyperparameters initialization of the VBrand
model used for the novelty detection of soil type in Section 5 of the main
paper.

  1. 1.

The truncation is set at $T=10$.

  2. 2.

The concentration parameter of the stick-breaking process is set to
$\gamma=5$.

  3. 3.

The hyperparameters of the Dirichlet distribution are all set to
$\alpha_{j}=0.1\;\;\forall j$, inducing the same probability of being in each
known component and of being a novelty.

  4. 4.

The $\mathcal{NIW}$ prior for the novelty locations is parameterized as
follows: $\boldsymbol{\mu}_{0}^{nov}$ is set equal to the grand mean of the
training set, to allow an hyperprior specification that does not rely on the
test set; the precision is given as ${\lambda}_{0}^{nov}=0.1$,the degrees of
freedom ${\nu}_{0}^{nov}=p+2$, where $p$ is the dimensionality of the dataset.
This is the minimum integer value that is allowed to have so that the
$\mathcal{IW}$ distribution has a finite mean. Finally, the variance is set to
be $(p+1)$ times the training sample covariance matrix. This choice implies
that the expected variance sampled from the $\mathcal{NIW}$ is an inflated
version of the overall sample covariance matrix of the training set.

  5. 5.

For the known components, the mean vector, and the variance-covariance
matrices are estimated via the MRCD procedure. To put high mass on these
estimates, for all the known groups, we set $\lambda_{k}^{obs}=200$ and
$\nu^{obs}_{k}=p+1+200$.