We only show the Type II modes with positive $\omega_{rm}$ in Figure 16. The
negative counterpart would overlap with the shown Type II curves in the
$\gamma_{m}$–$k$ plots and would appear as mirrored images of the
positive-$\omega_{rm}$ Type II curves about the $y$-axis in the
$\gamma$–$\omega_{rm}$ plots. This reconciles with the mathematical notion
that a bifurcation point should have two inbound and two outbound branches.

c) Very sharp bumps, e.g., $A=1.4,\Delta w/H=1.2$. The two types of modes
again have independent curves. In the $\gamma_{m}$–$k$ plot, the Type I curve
still has a maximum followed by a long tail, whereas the low-$k$ flank of the
Type II curve becomes concave up (showing a positive second derivative).

The smooth shift between the two types of instabilities implies a smooth
transition of the instability mechanism. In a transition from Type II to Type
I, the eigenmodal phase velocity $\omega_{rm}/k$ gradually reduces to zero,
and the co-rotation radius moves from outside to the center of the dust bump.
Also, cross-correlation calculations similar to those in Section 5.4.2 show
that the phase difference between $q_{1}$ and $S_{\rm bar1}$ rapidly increases
from a small angle to well over $90^{\circ}$; meanwhile, the phase difference
between the vortensity perturbation and advection diminishes to zero. This
indicates that, during the transition, advection takes over the role from the
baroclinic source in driving the instability in the dust bump, as might be
expected. Further studies of the mechanism behind the bifurcation would be
more complex and lie beyond the scope of this paper.

##  Appendix E Tentative propagation mechanism of the Type II DRWI

Figure 17: Schematic illustration of the tentative propagation mechanism of
the Type II DRWI in the rotating frame with $\omega_{rm}>0$. Top: vortensity
perturbation, advection and source patterns on the $x$-$y$ plane, extracted
from Figure 5. The patterns are separated by grey dotted lines into three
regions, namely, {1}⃝, {2}⃝, and {3}⃝ from left to right. Bottom: velocity
decomposition in each region. The phase and shear velocities combine to give
the modal azimuthal velocity $\omega_{rm}/k$. The two panels share the
$x$-axis.

One main feature of the Type II DRWI is a non-zero phase velocity. We have
explained the propagation of perturbation patterns in the dust bump in Section
5.4, but the Rossby waves are also important components of this instability.
Here, we tentatively explain the propagation of the Type II DRWI with the
entire mode in mind, as illustrated in Figure 17. With a simplified physical
picture, we focus on how the different components share one modal phase
velocity despite (or perhaps as a result of) the asymmetric eigenfunction
patterns.

We first apply the interpretation in Ono et al. (2016) to explain the
propagation of the gas perturbation patterns. In the following, we will use
the words “upper” and “lower” to describe location with higher or lower
$y$-coordinate. The leftmost region (labelled as {1}⃝) in the top panel in
Figure 17, where $f_{d}\ll 1$ and the background vortensity satisfies
$dq_{0}/dx<0$, illustrates a Rossby wave with a negative phase velocity
$v_{\rm phase}$ in $\hat{y}$. For sake of explanation, we temporarily ignore
any interaction between regions {1}⃝ and {2}⃝. Positive vortensity
perturbations imply clockwise loops of gas flow, while negative
counterclockwise. Such flows drive $\hat{x}$-direction advection of vortensity
along the $q_{0}$ gradient, which induces $y$-direction migration of the loops
themselves. For example, the flows at $ky\sim\pi/2$ and $ky\sim 3\pi/2$ in
region ① increase and decrease local vortensity respectively, the former
extending the lower side of the $+q_{1}$ loop and the latter paring the upper
side. The result is a downward shift of the perturbation patterns. However,
here the Rossby wave travels on a background with shear velocity $v_{\rm
shear}=v^{\prime}_{0y}-3\Omega_{0}x/2$, and thus the patterns overall move
upwards in the frame co-rotating with $x=0$.

The decomposition of the motion of the Rossby wave into $v_{\rm phase}$ and
$v_{\rm shear}$ is essentially built in the vortensity equation (34).
Formally, the phase velocity of the Rossby wave can be defined as

$v_{\rm phase}={\rm
Re}\left[i\frac{-(dq_{0}/dx)v^{\prime}_{1x}}{kq_{1}}\right]\ ,$ (94)

where the numerator appeared in Equation (34) as the radial advection term.
Replacing it with the full vortensity advection term would recover
$\omega_{rm}/k$ assuming conservation of vortensity, while replacing it with
the azimuthal advection term would give $v_{\rm shear}$. Equation (94) implies
that, when vortensity is conserved, $v_{\rm phase}$ will depend on the
strength of the radial vortensity advection relative to that of the vortensity
perturbation.

Application of the argument above on both regions {1}⃝ and {3}⃝ gives two
Rossby waves with $v_{\rm phase}$ in the opposite direction. In region {1}⃝
resides what we term the “major Rossby wave” (solid) with much more
conspicuous perturbation compared to the “minor Rossby wave” (dashed) in
region {3}⃝. However, the $x$-direction vortensity advection has comparable
strengths in the two regions, hence more influential relative to the minor
Rossby wave in terms of shifting the perturbation patterns. This leads to a
faster phase velocity in region ③, as illustrated in the bottom panel in
Figure 17. Then, taking the similar but opposite shearing velocity into
account, one observes that $|v_{\rm phase,{3}⃝}|>|v_{\rm
shear,{3}⃝}|\simeq|v_{\rm shear,{1}⃝}|>|v_{\rm phase,{1}⃝}|$.This explains why
the two regions share an overall modal phase velocity $\omega_{rm}/k=v_{\rm
phase,{3}⃝}-v_{\rm shear,{3}⃝}=v_{\rm shear,{1}⃝}-v_{\rm phase,{1}⃝}$. Our
explanation suggests a close link between the non-zero real part of the
eigenvalue and the asymmetry of the eigenfunction.

The propagation of $q_{1}$ patterns in region ② also fits in this picture. As
shown in Section 5.4, here the vortensity advection governs the propagation
with a $\sim 90^{\circ}$ spatial phase lead over $q_{1}$, thus inducing the
$\pm q_{1}$ loops in region ② to travel upwards. This process is similar to
that of the Rossby wave, except that the background $v_{\rm shear}$ is now in
the same direction as $v_{\rm phase}$. Here, the material flow advects
vortensity across different regions (e.g., the inclined arrows between region
① and ② in Figure 17), indicating a close interaction between the Rossby wave
and the dust bump.


# Introducing instance label correlation in multiple instance learning.
Application to cancer detection on histopathological images

Pablo Morales-Álvarez111Work done during the Margarita Salas fellowship
outgoing phase at Cambridge University, UK<EMAIL_ADDRESS>Arne
Schmidt  José Miguel Hernández-Lobato  Rafael Molina  Department of Computer
Science and AI, University of Granada, Spain  Department of Engineering,
University of Cambridge, UK

###### Abstract

In the last years, the weakly supervised paradigm of multiple instance
learning (MIL) has become very popular in many different areas. A paradigmatic
example is computational pathology, where the lack of patch-level labels for
whole-slide images prevents the application of supervised models.
Probabilistic MIL methods based on Gaussian Processes (GPs) have obtained
promising results due to their excellent uncertainty estimation capabilities.
However, these are general-purpose MIL methods that do not take into account
one important fact: in (histopathological) images, the labels of neighboring
patches are expected to be correlated. In this work, we extend a state-of-the-
art GP-based MIL method, which is called VGPMIL-PR, to exploit such
correlation. To do so, we develop a novel coupling term inspired by the
statistical physics Ising model. We use variational inference to estimate all
the model parameters. Interestingly, the VGPMIL-PR formulation is recovered
when the weight that regulates the strength of the Ising term vanishes. The
performance of the proposed method is assessed in two real-world problems of
prostate cancer detection. We show that our model achieves better results than
other state-of-the-art probabilistic MIL methods. We also provide different
visualizations and analysis to gain insights into the influence of the novel
Ising term. These insights are expected to facilitate the application of the
proposed model to other research areas.

###### keywords:

Multiple Instance Learning, Gaussian Processes , Ising model , Variational
Inference , Whole Slide Images , Histopathology

##  1 Introduction

Multiple instance learning (MIL) has caught great attention in fields where
there is a challenging lack of labelled data. Although it has been applied in
many different areas [1], we will focus on the case of computational
pathology. In the last years, thanks to the increasing digitalization of
whole-slide images (WSIs), the field of computational pathology is developing
computer-aided diagnosis systems based on machine learning for cancer
detection [2, 3]. The goal of computational pathology is to provide a fast and
reliable diagnosis for the most prototypical cases, letting the pathologists
focus on the most challenging ones. Ultimately, this will enable a much wider
access to early cancer diagnosis [4].

In order to make accurate predictions, machine learning classification methods
need to be trained using a labelled set of instances [5]. In the case of
computational pathology, these instances are typically patches from the WSIs
(and not the complete images themselves) [6, 7, 8]. The reason for this is
twofold: i) it is useful to have predictions at patch level in order to know
where exactly in the image the cancer is located, and ii) WSIs are extremely
large and cannot be directly fed to a classifier. As a consequence, notice
that expert pathologists must label _every single patch_ in the training data
as cancerous or not (we will consider the binary problem cancer/no-cancer
throughout this work). Given the large number of patches and the limited
availability of pathologists, this becomes a daunting task in real practice
[9].

To address this problem, different weakly supervised learning paradigms have
been proposed in recent years. Here we focus on MIL, which has become very
popular in the medical domain [10, 8]. The idea in MIL is that instances are
grouped in bags, and only bag labels are needed for training. In the case of
WSIs, all the patches coming from the same image are considered a bag.
Therefore, the labelling workload on pathologists decreases enormously: from
labelling every single patch, to only labelling the complete WSI as cancerous
or not.

Different machine learning algorithms have been developed to learn under the
MIL setting. Notice that dealing with uncertainty is essential in MIL models,
since instance-level labels are unknown. To deal with uncertainties, different
probabilistic methods have been developed, such as Dirichlet Process Mixture
Models [11], Markov chain [12], Monte-Carlo chain [13, 14] and Gaussian
Processes (GPs) [15, 16]. In particular, GPs have attracted plenty of
attention in the last years, due to their expressive power and their capacity
to handle uncertainty in a principled manner. Moreover, we are interested in
this type of probabilistic models, since they will allow for introducing
correlations in a theoretically sound way.

Among GP-based MIL methods, we will focus on the two most successful ones:
VGPMIL and VGPMIL-PR. VGPMIL [17] was proposed in 2017 to overcome the
limitations of two earlier formulations [15, 16] (namely, the use of the
inefficient Laplace approximation and the impossibility to obtain instance-
level predictions, respectively). In short, VGPMIL relies on variational
inference and allows for closed-form updates of its parameters. However, the
use of the logistic function implies that VGPMIL needs to resort to a
theoretical approximation during inference (namely, the Jaakola bound [17, Eq.
(10)]). As shown in [18], such approximation hurts predictive performance in
practice. As an alternative, the authors of [18] propose the utilization of
the probit function, which removes the need for the aforementioned
approximation. This method, which will be referred to as VGPMIL-PR, is
considered the current state of the art among probabilistic MIL approaches.

Methods such as VGPMIL and VGPMIL-PR are general-purpose MIL models that can
be used in any MIL problem (that is, whenever the label is known only at bag
level, see different use-cases in [17, 19]). However, the underlying MIL
assumption that the labels of the instances in a bag are independent of each
other is unrealistic in many real problems. For example, in the particular
case of WSI images (and in many image-related MIL problems), the labels of
neighboring patches are expected to be correlated [20]. We hypothesize that
the predictive performance of MIL methods can be enhanced by incorporating
this type of prior knowledge into the model.

In this work, we introduce a novel GP-based MIL algorithm that takes into
account the correlation between the labels of neighboring patches, and we
apply it to the real-world problem of prostate cancer detection on
histopathological images. We model the correlation through a coupling term
inspired by the Ising model [5, Section 19.4.1], an statistical physics method
that has found several applications in computer vision [20, 21]. Our GP-MIL
modeling builds on VGPMIL-PR, so our method will be referred to as VGPMIL-PR-I
(Ising). In VGPMIL-PR-I, a hyperparameter $\lambda$ regulates the influence of
the Ising-inspired terms. Variational inference is used to estimate the model
parameters, and the update formulas of VGPMIL-PR are recovered when
$\lambda\to 0$ (that is, when the influence of the coupling term vanishes). In
the experimental section, we show that VGPMIL-PR-I outperforms the state-of-
the-art GP-based MIL approaches VGPMIL and VGPMIL-PR when predicting at both
instance and bag levels, while keeping an analogous computational cost.
Moreover, to gain insights into the influence of the new coupling term, we
analyze the role of $\lambda$, and provide several visualizations for the
predictions.

The rest of the paper is organized as follows. Section 2 presents the
probabilistic model and inference for the novel VGPMIL-PR-I. Closely related
methods such as VGPMIL and VGPMIL-PR are also discussed in this section.
Section 3 focuses on the empirical evaluation of the model, including the data
description, the experimental framework, and the discussion of results.
Section 4 provides the main conclusions and some future outlook.

##  2 Probabilistic model and inference

In this section we present the theoretical description for VGPMIL-PR-I.
Specifically, Section 2.1 explains the problem formulation and the main
notation. Section 2.2 explains the closely-related methods VGPMIL and VGPMIL-
PR, which are at the base of our formulation. Section 2.3 introduces the novel
coupling term that accounts for patch label correlation, which is used to
define VGPMIL-PR-I. Section 2.4 shows how to perform variational inference to
estimate the parameters in VGPMIL-PR-I. Section 2.5 explains the procedure to
make predictions at both instance and bag levels.

###  2.1 Notation and problem formulation

Our notation follows the state-of-the-art work [22]. The training data is
given by a set of bags ${\mathbf{X}}=\\{{\mathbf{X}}_{b}\\}_{b\in\mathcal{B}}$
and their corresponding labels ${\mathbf{y}}=\\{y_{b}\\}_{b\in\mathcal{B}}$.
We deal with a binary problem, i.e. $y_{b}\in\\{0,1\\}$. Each bag
${\mathbf{X}}_{b}=\\{{\mathbf{x}}_{i}\\}_{i\in b}$ contains $|b|$ instances,
i.e. $b=\\{i_{1},\dots,i_{|b|}\\}\subseteq[N]$ ($N$ is the total amount of
instances). Notice that different bags may have different amounts of
instances. Each instance ${\mathbf{x}}_{i}$ is given by a vector in
$\mathbb{R}^{D}$. In the MIL setting, one assumes that each instance has its
(unknown) label $h_{i}\in\\{0,1\\}$. We write ${\mathbf{h}}_{b}$ for the
labels of all the instances belonging to bag $b$. The MIL labelling assumption
dictates that a bag is considered positive (class 1) if at least one of its
instances is positive. Mathematically, this is

${\mathrm{p}}(y_{b}|{\mathbf{h}}_{b})=1[y_{b}=\max_{i\in b}h_{i}],$ (1)

where $1[\cdot]$ is the indicator function (i.e. it equals one when its
argument is true and zero otherwise). Finally, we will collectively denote
${\mathbf{h}}=\\{{\mathbf{h}}_{b}\\}_{b\in\mathcal{B}}$.

In the case of WSIs, each ${\mathbf{X}}_{b}$ is an image, which is composed of
its patches $\\{{\mathbf{x}}_{i}\\}_{i\in b}$. Each patch has an unknown label
$h_{i}$ (0 for non-cancerous and 1 for cancerous), and we only have access to
the bag label $y_{b}$ (whether the image is cancerous or not, i.e. whether it
contains at least one patch that is cancerous).

The goal in MIL is to train a model based only on bag labels
$\\{y_{b}\\}_{b\in\mathcal{B}}$. And such model must be able to predict at
both instance and bag levels. That is, given a previously unseen instance
${\mathbf{x}}^{\star}\in\mathbb{R}^{D}$, we are interested in the probability
${\mathrm{p}}(h^{\star}=1)$. Likewise, given a previously unseen complete bag
${\mathbf{X}}^{\star}$, we are interested in ${\mathrm{p}}(y^{\star}=1)$.

###  2.2 Background: VGPMIL and VGPMIL-PR

As mentioned in the introduction, our model is inspired by two closely related
methods: VGPMIL and VGPMIL-PR. To understand our contribution, it is essential
to fully understand their formulations, which we explain next.

####  2.2.1 VGPMIL formulation

VGPMIL was introduced in [17]. The idea is to consider a sparse GP
classification model [23] to describe the relationship between instance
features ${\mathbf{X}}$ and their (unknown) labels ${\mathbf{h}}$. Then, an
additional bag likelihood must be considered to model the (observed) bag
labels ${\mathbf{y}}$ given the instance labels ${\mathbf{h}}$. Both
components are described next.

The sparse GP classification model. Instances ${\mathbf{x}}_{i}$ are
associated latent variables $f_{i}\in\mathbb{R}$ which are modelled through a
GP, $f\sim\mathcal{GP}(0,\kappa)$. We write $\kappa$ for the GP kernel, which
encodes the properties of the considered functions. Then, the instance labels
$h_{i}$ are defined from $f_{i}$ through a classification likelihood $\nu$:

${\mathrm{p}}(h_{i}|f_{i})=\nu(f_{i})^{h_{i}}(1-\nu(f_{i}))^{1-h_{i}}.$ (2)

Specifically, VGPMIL uses the logistic function $\nu(x)=(1+e^{-x})^{-1}$.
Intuitively, a large (resp. low) value of $f_{i}$ implies that the class is
likely to be one (resp. zero). Moreover, since standard GPs scale poorly with
the number of training instances $N$, VGPMIL makes use of sparse GPs [23],
which summarize the training data through $M\ll N$ inducing points. These
inducing points ${\mathbf{u}}=\\{u_{1},\dots,u_{M}\\}$ represent the value of
the GP at some inducing points locations
${\mathbf{Z}}=\\{{\mathbf{z}}_{1},\dots,{\mathbf{z}}_{M}\\}$ (just like
${\mathbf{f}}=\\{f_{1},\dots,f_{N}\\}$ are the GP values at
${\mathbf{X}}=\\{{\mathbf{x}}_{1},\dots,{\mathbf{x}}_{N}\\}$). Therefore, the
distributions of ${\mathbf{u}}$ and ${\mathbf{f}}|{\mathbf{u}}$ are:

$\displaystyle{\mathrm{p}}({\mathbf{u}})$
$\displaystyle=\mathcal{N}({\mathbf{u}}|\mathbf{0},{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}),$
(3) $\displaystyle{\mathrm{p}}({\mathbf{f}}|{\mathbf{u}})$
$\displaystyle=\mathcal{N}({\mathbf{f}}|{\mathbf{K}}_{{\mathbf{X}}{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{u}},\mathcal{K}).$
(4)

The expression of $\mathcal{K}$ is given by the particular sparse GP approach
used in VGPMIL, which is FITC [23], so we have
$\mathcal{K}=\operatorname*{diag}({\mathbf{K}}_{{\mathbf{X}}{\mathbf{X}}}-{\mathbf{K}}_{{\mathbf{X}}{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{X}}})$.
As is standard in GP literature, we are writing
${\mathbf{K}}_{{\mathbf{X}}{\mathbf{X}}}$ for the $N\times N$ covariance
matrix
${\mathbf{K}}_{{\mathbf{X}}{\mathbf{X}}}=(\kappa({\mathbf{x}}_{i},{\mathbf{x}}_{j}))_{1\leq
i,j\leq N}$. The definitions for ${\mathbf{K}}_{{\mathbf{X}}{\mathbf{Z}}}$ and
${\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}$ are analogous.

The bag likelihood. VGPMIL introduces the following parameterization to model
the bag labels from the instance labels:

${\mathrm{p}}(y_{b}|{\mathbf{h}}_{b})=\frac{H^{G_{b}}}{H+1},$ (5)

where $G_{b}:=1[y_{b}=\max_{i\in b}h_{i}]$ and $H$ is a large and fixed value
(in their examples, they use $H=100$). Eq. (5) approximates the MIL assumption
introduced in eq. (1): if some instance label $h_{i}$ is one, then the bag
label $y_{b}$ is one with very high probability (namely, with probability
$\frac{H}{H+1}$). Otherwise (that is, if all instance labels $h_{i}$ are
zero), the bag label is one with very low probability (namely,
$\frac{1}{H+1}$).

In summary, VGPMIL is given by eqs. (3), (4), (2) and (5). For additional
details, the interested reader is referred to the original work [17].

####  2.2.2 VGPMIL-PR formulation

(a) VGPMIL-PR | $u_{m}$$f_{i}$$m_{i}$$y_{b}$$m=1,\dots,M$$i\in b$$b\in\mathcal{B}$
---|---
(b) VGPMIL-PR-I | $u_{m}$$f_{i}$$m_{i}$$y_{b}$$m=1,\dots,M$$i\in b$$i\in b$$b\in\mathcal{B}$
Figure 1: Probabilistic graphical model for VGPMIL-PR (a) and VGPMIL-PR-I (b).
Gray nodes are observed variables, and white ones are latent variables to be
estimated. The main difference is that VGPMIL-PR-I introduces correlation
between instances in the same bag. Therefore, the distribution of
${\mathbf{m}}_{b}$ given ${\mathbf{f}}_{b}$ does not factorize across
instances. The correlation is introduced through a novel term inspired by the
Ising model, see Section 2.3 for details.

VGPMIL-PR was recently proposed in [18] as an improvement over VGPMIL. Namely,
the logistic function used by VGPMIL in eq. (2) is not conjugate with the
Gaussian distribution coming from the GP, recall eqs. (3)–(4). This means
that, in order to achieve mathematical tractability, VGPMIL needs to resort to
the Jaakola bound [17, Eq. (10)]. However, the use of this bound introduces an
_approximation_ in the training objective. As shown in [18], such
approximation damages the predictive performance in practice. Consequently,
the authors of [18] introduce an alternative formulation based on the probit
function, VGPMIL-PR. They show that, via a variable augmentation approach,
VGPMIL-PR allows for directly optimizing the training objective (without
approximations).

More specifically, VGPMIL-PR uses the probit function
$\nu(x)=\int_{-\infty}^{x}\mathcal{N}(t|0,1)\mathrm{d}t$ in eq. (2). Also, the
bag likelihood is given by eq. (1) (instead of eq. (5)). Then, to circumvent
the need for approximations, VGPMIL-PR leverages a variable augmentation
approach [24]. Namely, for each instance we introduce a new variable
$m_{i}\in\mathbb{R}$ between $f_{i}$ and $h_{i}$, which is defined as
$m_{i}\sim\mathcal{N}(f_{i},1)$. Since we are using a probit likelihood, we
have that $h_{i}=1[m_{i}>0]$. Analogously to the rest of variables, we write
${\mathbf{m}}_{b}=\\{m_{i}\\}_{i\in b}$ for all the $m_{i}$’s inside bag $b$,
and we use ${\mathbf{m}}=\\{{\mathbf{m}}_{b}\\}_{b\in\mathcal{B}}$ to
collectively denote all the $m_{i}$’s in the model. Then, by marginalizing out
${\mathbf{h}}$, we have:

$\displaystyle{\mathrm{p}}({\mathbf{m}}|{\mathbf{f}})$
$\displaystyle=\prod_{b}{\mathrm{p}}({\mathbf{m}}_{b}|{\mathbf{f}}_{b})=\prod_{b}\mathcal{N}({\mathbf{m}}_{b}|{\mathbf{f}}_{b},I),$
(6) $\displaystyle{\mathrm{p}}({\mathbf{y}}|{\mathbf{m}})$
$\displaystyle=\prod_{b}{\mathrm{p}}(y_{b}|{\mathbf{m}}_{b}),$ (7)

where $I$ is the identity matrix (of size $|b|$) and
${\mathrm{p}}(y_{b}=0|{\mathbf{m}}_{b})=\prod_{i\in b}1[m_{i}<0]$.
Importantly, these augmented variables ${\mathbf{m}}$ will prove extremely
helpful to introduce the Ising correlation in the next section.

In summary, VGPMIL-PR is given by eqs. (3), (4), (6) and (7). Figure 1(a)
shows the probabilistic graphical model for VGPMIL-PR.

###  2.3 Correlating patch labels: VGPMIL-PR-I

$P_{1}$$P_{2}$$P_{3}$$P_{4}$$P_{5}$ Figure 2: A simplified illustration of an
image composed by five patches: $P1,\dots,P5$.

As mentioned in the introduction, VGPMIL and VGPMIL-PR are general MIL
approaches that can be used in any MIL problem. Indeed, there are plenty of
applications where MIL methods can be used. For instance, think of a
recommendation system where a reviewer has not evaluated every single item in
the database, but has reviewed “groups” of them (e.g., he/she likes science-
fiction movies, although he/she has not rated individual movies). Consider
also a task of anomaly detection in which we do not have labels for individual
transactions, but we only know whether there was some anomalous behavior in a
certain period of time (which contains many different transactions).

Here we focus in the particular use-case of images, where bags are images and
their instances are their patches. In this case, there exists very valuable
information coming from the structure of the image itself, which can be
exploited in the model. For example, it is natural to think that neighboring
patches in the same image are likely to have similar labels. The main goal of
this work is to introduce such correlation into the VGPMIL-PR formulation. To
do so, we are inspired by the Ising model.

The novel coupling term. The Ising model arose from statistical physics to
describe the behavior of magnets. In some magnets, called ferro-magnets,
neighboring spins tend to line up in the same direction, whereas in other
kinds of magnets, called anti-ferromagnets, the spins are repelled from their
neighbors [5]. This type of interactions based on the Ising model have been
used previously in machine learning and computer vision to describe
relationships between pixels of an image, see e.g. [20, 21, 5]. However, to
the best of our knowledge, they have never been used in the context of MIL.

Our first idea was to consider an Ising model over the patch labels of each
image, $\\{h_{i}\\}_{i\in b}$. However, inference proved very challenging in
this case, due to the non-conjugacy of the Ising model and the GP-based MIL
formulation. As an alternative, we considered a continuous counterpart of the
Ising model over the variables ${\mathbf{m}}_{b}=\\{m_{i}\\}_{i\in b}$
introduced in VGPMIL-PR, which are directly related to the patch labels
(recall from Section 2.2.2 that $h_{i}=1[m_{i}>0]$). Notice that such
continuous version of the Ising model corresponds to the well-known
Conditional Autoregression (CAR) [25]. Importantly, as we will see in the rest
of this section, this alternative formulation yields a Gaussian distribution
on ${\mathbf{m}}_{b}$, which can be treated analytically together with the GP-
based MIL model.

Specifically, we consider the following coupling term for each image, which is
defined over the augmented variables ${\mathbf{m}}_{b}$, recall Section 2.2.2:

$\displaystyle\mathcal{C}({\mathbf{m}}_{b})=$
$\displaystyle\exp\left(-\frac{\lambda}{2}\cdot\sum_{\begin{subarray}{c}i,j\in
b\\\ i<j\end{subarray}}\mathbf{1}[i,j\textrm{ are
contiguous}]\cdot(m_{i}-m_{j})^{2}\right)=$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\lambda}{2}{\mathbf{m}}_{b}^{\intercal}{\mathbf{C}}_{b}{\mathbf{m}}_{b}\right).$
(8)

Notice that $\mathcal{C}({\mathbf{m}}_{b})$ is always in the range $[0,1]$,
and it becomes close to zero when the value of $m$ is very different for
neighboring patches. For the second equality in eq. (2.3), notice that the sum
only produces quadratic terms in $m$, so it can be written as
${\mathbf{m}}_{b}^{\intercal}{\mathbf{C}}_{b}{\mathbf{m}}_{b}$ for some
positive semidefinite matrix ${\mathbf{C}}_{b}$.

Since $m$ determines the label of each patch (recall from Section 2.2.2 that
$h_{i}=1[m_{i}>0]$), the term $\mathcal{C}({\mathbf{m}}_{b})$ can be used to
favor “smoothness” in the labels associated to the different patches. Also,
the hyperparameter $\lambda$, which can be set to any non-negative value,
regulates the strength of the coupling term: the larger $\lambda$, the more
importance is given to differences in $m$. For example, when $\lambda=0$,
$\mathcal{C}({\mathbf{m}}_{b})$ becomes constant and it does not account for
correlation between patch labels.

An example of $\mathcal{C}({\mathbf{m}}_{b})$. To illustrate the proposed
coupling term, consider an image with five patches $P_{1},\dots,P_{5}$
distributed as in Figure 2. In this case, the quadratic terms of
${\mathbf{m}}_{b}$ are:

$(m_{1}-m_{2})^{2}+(m_{1}-m_{3})^{2}+(m_{2}-m_{4})^{2}+(m_{3}-m_{4})^{2}+(m_{4}-m_{5})^{2},$
(9)

and therefore we have:

${\mathbf{C}}_{b}=\begin{pmatrix}2&-1&-1&0&0\\\ -1&2&0&-1&0\\\ -1&0&2&-1&0\\\
0&-1&-1&3&-1\\\ 0&0&0&-1&1\\\ \end{pmatrix}.$ (10)

In general, it is easy to compute the matrix ${\mathbf{C}}_{b}$ for any given
image. Notice that it is always a positive semidefinite matrix, and thus it is
associated to a (singular) normal distribution.

The VGPMIL-PR-I formulation. To introduce the new coupling term
$\mathcal{C}({\mathbf{m}}_{b})$ in the MIL formulation, we modify eq. (6) and
define:

${\mathrm{p}}({\mathbf{m}}|{\mathbf{f}})\propto\prod_{b}\mathcal{C}({\mathbf{m}}_{b})\cdot\mathcal{N}({\mathbf{m}}_{b}|{\mathbf{f}}_{b},I).$
(11)

Notice that the probability of a configuration ${\mathbf{m}}_{b}$ is
proportional to the coupling term $\mathcal{C}({\mathbf{m}}_{b})$, which
favors smoothness across labels of neighboring patches. Decisively, since both
$\mathcal{C}({\mathbf{m}}_{b})$ and
$\mathcal{N}({\mathbf{m}}_{b}|{\mathbf{f}}_{b},I)$ only contain (the
exponential of) quadratic terms in ${\mathbf{m}}_{b}$, the new distribution
can be written as a Gaussian:

${\mathrm{p}}({\mathbf{m}}|{\mathbf{f}})=\prod_{b}\mathcal{N}({\mathbf{m}}_{b}|{\boldsymbol{\Sigma}}_{b}{\mathbf{f}}_{b},{\boldsymbol{\Sigma}}_{b}),$
(12)

with ${\boldsymbol{\Sigma}}_{b}=(\lambda{\mathbf{C}}_{b}+I)^{-1}$. Notice that
this new formulation provides a generalization of VGPMIL-PR. Namely, when
$\lambda\to 0$, we have that ${\boldsymbol{\Sigma}}_{b}\to I$ and we recover
eq. (6).

In summary, the proposed model is given by eqs. (3), (4), (12), and (7).
Notice also that, instead of FITC, in eq. (4) we leverage the more recent
sparse GP approach introduced in [26]. Basically, this means that
$\mathcal{K}={\mathbf{K}}_{{\mathbf{X}}{\mathbf{X}}}-{\mathbf{K}}_{{\mathbf{X}}{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{X}}}$
in eq. (4). Our method will be referred to as VGPMIL-PR-I (I denotes Ising).
The probabilistic graphical model is depicted in Figure 1(b).

###  2.4 Variational inference

In order to make inference in the proposed model, we need to compute the
posterior distribution
${\mathrm{p}}({\mathbf{u}},{\mathbf{f}},{\mathbf{m}}|{\mathbf{y}})$. However,
this is not analytically tractable due to the definition of the bag likelihood
in eq. (7), which depends on the sign of the $m_{i}$’s. Following [17] and
[18], we leverage standard mean-field variational inference (VI) theory [27,
Section 10.1.1] to calculate an approximate posterior distribution that
factorizes as

${\mathrm{q}}({\mathbf{u}},{\mathbf{f}},{\mathbf{m}})={\mathrm{q}}({\mathbf{u}}){\mathrm{p}}({\mathbf{f}}|{\mathbf{u}}){\mathrm{q}}({\mathbf{m}}).$
(13)

Applying the well-known mean-field VI update equation [27, Eq. (10.9)], we
have that ${\mathrm{q}}({\mathbf{u}})$ and ${\mathrm{q}}({\mathbf{m}})$ can be
iteratively computed as

$\displaystyle{\mathrm{q}}({\mathbf{u}})$
$\displaystyle=\mathcal{N}({\mathbf{u}}|{\boldsymbol{\mu}}^{u},{\boldsymbol{\Sigma}}^{u}),$
(14) $\displaystyle{\mathrm{q}}({\mathbf{m}})$
$\displaystyle\propto\prod_{b}{\mathrm{p}}(y_{b}|{\mathbf{m}}_{b})\mathcal{N}({\mathbf{m}}_{b}|{\boldsymbol{\mu}}^{{\mathbf{m}}_{b}},{\boldsymbol{\Sigma}}_{b}),$
(15)

where

$\displaystyle{\boldsymbol{\Sigma}}^{u}$
$\displaystyle=\left({\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}+{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{X}}}{\boldsymbol{\Sigma}}{\mathbf{K}}_{{\mathbf{X}}{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}\right)^{-1},$
(16) $\displaystyle{\boldsymbol{\mu}}^{u}$
$\displaystyle={\boldsymbol{\Sigma}}^{u}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{X}}}\mathbb{E}_{{\mathrm{q}}({\mathbf{m}})}({\mathbf{m}}),$
(17)

and

${\boldsymbol{\mu}}^{{\mathbf{m}}_{b}}={\boldsymbol{\Sigma}}_{b}{\mathbf{K}}_{b{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\boldsymbol{\mu}}^{u}.$
(18)

Here we are writing ${\boldsymbol{\Sigma}}$ for the $N\times N$ block-diagonal
matrix that contains all the ${\boldsymbol{\Sigma}}_{b}$’s, $b\in\mathcal{B}$.
Also, we are writing ${\mathbf{K}}_{b{\mathbf{Z}}}$ for the $|b|\times M$
matrix of covariances between ${\mathbf{X}}_{b}$ and ${\mathbf{Z}}$. Very
importantly, notice that these update rules generalize those derived in [18]
for VGPMIL-PR. Namely, when $\lambda\to 0$, we have that
${\boldsymbol{\Sigma}}_{b},{\boldsymbol{\Sigma}}\to I$, and then eqs.
(14)–(18) match eqs.(15)–(19) in [18].

All the computations involved in eqs. (14)–(18) are straightforward, except
for $\mathbb{E}_{{\mathrm{q}}({\mathbf{m}})}({\mathbf{m}})$. Indeed, each
${\mathrm{q}}({\mathbf{m}}_{b})$ is a multivariate Gaussian truncated to
$(-\infty,0)^{|b|}$ (or $\mathbb{R}^{|b|}\setminus(-\infty,0)^{|b|}$,
depending on whether $y_{b}=0$ or $y_{b}=1$, respectively). It is well-known
that the expectation of a truncated multivariate Gaussian cannot be obtained
in closed-form [28]. Notice that this is not an issue for VGPMIL-PR [18],
where the absence of Ising terms implies dealing with _univariate_ Gaussians,
whose expectations can be analytically computed.

To overcome the problem, we first tried to leverage numerical methods proposed
in [29] to approximate the expectation for the multivariate truncated case.
However, these methods proved computationally too expensive to be integrated
within our iterative calculation of ${\mathrm{q}}({\mathbf{u}})$ and
${\mathrm{q}}({\mathbf{m}})$. Therefore, we decided to approximate the
multivariate Gaussian
$\mathcal{N}({\mathbf{m}}_{b}|{\boldsymbol{\mu}}^{{\mathbf{m}}_{b}},{\boldsymbol{\Sigma}}_{b})$
by the factorized
$\mathcal{N}({\mathbf{m}}_{b}|{\boldsymbol{\mu}}^{{\mathbf{m}}_{b}},\mathrm{diag}({\boldsymbol{\Sigma}}_{b}))$
and utilize the expression for one-dimensional truncated Gaussians. Although
such approximation reduces the influence of the Ising correlation at this
specific computation, notice that the coupling terms, which are included in
${\boldsymbol{\Sigma}}_{b}$, affect the update equations in more places across
eqs. (16)–(18).

Specifically, the expression for each
$\mathbb{E}_{{\mathrm{q}}({\mathbf{m}}_{b})}({\mathbf{m}}_{b})$ is as follows.
For bags with $y_{b}=0$, we have that each ${\mathrm{q}}(m_{i})$, $i\in b$, is
a univariate normal distribution
$\mathcal{N}(({\boldsymbol{\mu}}^{{\mathbf{m}}_{b}})_{i},({\boldsymbol{\Sigma}}_{b})_{ii})$
truncated to $(-\infty,0)$. The expectation of such a distribution is well-
known and can be obtained in closed-form [30]:

$E_{i}=\mu_{i}-\frac{\phi(\mu_{i}/\sigma_{i})}{1-\Phi(\mu_{i}/\sigma_{i})}\sigma_{i},$
(19)

where $\phi$ and $\Phi$ are, respectively, the density and cumulative
distribution functions of a standard Gaussian $\mathcal{N}(0,1)$ (recall that
both are efficiently implemented in standard software packages such as
Python’s Scipy). We have also abbreviated
$\mu_{i}=({\boldsymbol{\mu}}^{{\mathbf{m}}_{b}})_{i}$ and
$\sigma_{i}=\sqrt{({\boldsymbol{\Sigma}}_{b})_{ii}}$. For bags with $y_{b}=1$,
we proceed analogously to [19] to obtain the normalization constant $Z$ of the
distribution of interest (that is, the factorized Gaussian
$\mathcal{N}({\mathbf{m}}_{b}|{\boldsymbol{\mu}}^{{\mathbf{m}}_{b}},\mathrm{diag}({\boldsymbol{\Sigma}}_{b}))$
truncated to $\mathbb{R}^{|b|}-(-\infty,0)^{|b|}$). Then, the expectation of
each ${\mathrm{q}}(m_{i})$, $i\in b$, is given by:

$\mathbb{E}_{{\mathrm{q}}(m_{i})}(m_{i})=\frac{\mu_{i}-(1-Z)E_{i}}{Z},$ (20)

where $Z=1-\prod_{i\in b}(1-\Phi(\mu_{i}/\sigma_{i}))$, $E_{i}$ is given by
eq. (19), and we are again abbreviating
$\mu_{i}=({\boldsymbol{\mu}}^{{\mathbf{m}}_{b}})_{i}$ and
$\sigma_{i}=\sqrt{({\boldsymbol{\Sigma}}_{b})_{ii}}$.

The full training algorithm is summarized in Algorithm 1. It is an iterative
process that alternates the updates between ${\mathrm{q}}({\mathbf{u}})$ and
${\mathrm{q}}({\mathbf{m}})$. Details on the GP kernel and initializations
used in this work are provided in Section 3.1. The code for the proposed
method will be publicly available upon acceptance of the paper.

Input :

Bags ${\mathbf{X}}=\\{{\mathbf{X}}_{b}\\}_{b\in\mathcal{B}}$ and bag labels
${\mathbf{y}}=\\{y_{b}\\}_{b\in\mathcal{B}}$.

Calculate the matrices ${\mathbf{C}}_{b}$ that account for the instance
correlation inside each bag $b\in\mathcal{B}$, recall eq. (2.3) and the
example at eq. (10).

Initialize GP kernel parameters and inducing points locations, as well as the
posterior distributions ${\mathrm{q}}({\mathbf{u}})$ and
${\mathrm{q}}({\mathbf{m}})$. Details on initializations in the text.

foreach  _$\mathrm{iteration}\ t=1,\dots,T$_ do

Update ${\mathrm{q}}({\mathbf{u}})$, using eqs. (14), (16) and (17). Update
${\mathrm{q}}({\mathbf{m}})$ and obtain
$\mathbb{E}_{{\mathrm{q}}({\mathbf{m}})}({\mathbf{m}})$, using eqs. (15),
(18), (19), and (20).

Output : Posterior distributions ${\mathrm{q}}({\mathbf{u}})$ and
${\mathrm{q}}({\mathbf{m}})$.

Algorithm 1 Training procedure for VGPMIL-PR-I.

###  2.5 Making predictions

Suppose we are given a new bag
${\mathbf{X}}^{\star}=\\{{\mathbf{x}}^{\star}_{i}\\}_{i\in b_{\star}}$. As
explained at the end of section 2.1, we are interested in both instance-level
and bag-level predictions. For this, we first need to compute the predictive
distributions over ${\mathbf{m}}^{\star}$.

By using the learned posterior ${\mathrm{q}}({\mathbf{u}})$ along with
${\mathrm{p}}({\mathbf{f}}|{\mathbf{u}})$, we can obtain the joint
distribution over ${\mathbf{f}}^{\star}$:

${\mathrm{p}}({\mathbf{f}}^{\star})=\int{\mathrm{p}}({\mathbf{f}}^{\star}|{\mathbf{u}}){\mathrm{p}}({\mathbf{u}}){\mathrm{d}}{\mathbf{u}}=\mathcal{N}({\mathbf{f}}^{\star}|{\boldsymbol{\mu}}^{\star},{\mathbf{S}}^{\star}),$
(21)

with ${\boldsymbol{\mu}}^{\star}$ and ${\mathbf{S}}^{\star}$ given by the
standard sparse GP predictions:

${\boldsymbol{\mu}}^{\star}={\mathbf{K}}_{\star{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\boldsymbol{\mu}}^{u},\quad{\mathbf{S}}^{\star}={\mathbf{K}}_{\star\star}-{\mathbf{K}}_{\star{\mathbf{Z}}}{\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}({\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}-{\boldsymbol{\Sigma}}^{u}){\mathbf{K}}_{{\mathbf{Z}}{\mathbf{Z}}}^{-1}{\mathbf{K}}_{{\mathbf{Z}}\star}.$
(22)

Here, ${\boldsymbol{\mu}}^{u}$ and ${\boldsymbol{\Sigma}}^{u}$ are the
parameters learned during training, recall eq. (16) and (17). Naturally, the
subscript $\star$ in the kernel matrices ${\mathbf{K}}$ indicates that we are
using the new bag ${\mathbf{X}}^{\star}$. Then, since the distribution
${\mathrm{p}}({\mathbf{m}}|{\mathbf{f}})$ is also Gaussian, recall eq. (12),
we can compute the joint distribution over ${\mathbf{m}}^{\star}$ in closed-
form:

${\mathrm{p}}({\mathbf{m}}^{\star})=\int{\mathrm{p}}({\mathbf{m}}^{\star}|{\mathbf{f}}^{\star}){\mathrm{p}}({\mathbf{f}}^{\star}){\mathrm{d}}{\mathbf{f}}^{\star}=\mathcal{N}({\mathbf{m}}^{\star}|{\boldsymbol{\mu}}_{m}^{\star},{\mathbf{S}}_{m}^{\star}),$
(23)

with ${\boldsymbol{\mu}}_{m}^{\star}$ and ${\mathbf{S}}_{m}^{\star}$ given by

${\boldsymbol{\mu}}_{m}^{\star}={\boldsymbol{\Sigma}}_{\star}\cdot{\boldsymbol{\mu}}^{\star},\quad{\mathbf{S}}_{m}^{\star}={\boldsymbol{\Sigma}}_{\star}+{\boldsymbol{\Sigma}}_{\star}\cdot{\mathbf{S}}^{\star}\cdot{\boldsymbol{\Sigma}}_{\star}^{\intercal}.$
(24)

Here, ${\boldsymbol{\mu}}^{\star}$ and ${\mathbf{S}}^{\star}$ are given by eq.
(22), and ${\boldsymbol{\Sigma}}_{\star}$ is the matrix that accounts for
correlation among instances in the test bag ${\mathbf{X}}^{\star}$, which is
defined analogously to the training case, recall matrix
${\boldsymbol{\Sigma}}_{b}$ in eq. (12).

Once we have the joint distribution over ${\mathbf{m}}^{\star}$, the instance-
level and bag-level predictions are given as:

$\displaystyle{\mathrm{p}}(h_{i}^{\star}=1)$
$\displaystyle={\mathrm{p}}(m_{i}^{\star}>0)=\Phi\left(({\boldsymbol{\mu}}_{m}^{\star})_{i}/\sqrt{({\mathbf{S}}_{m}^{\star})_{ii}}\right),$
(25) $\displaystyle{\mathrm{p}}(y^{\star}=1)$
$\displaystyle=1-\int_{{\mathbf{m}}^{\star}\in(-\infty,0)^{|b_{\star}|}}{\mathrm{p}}({\mathbf{m}}^{\star}){\mathrm{d}}{\mathbf{m}}^{\star}.$
(26)

Notice that the integral in eq. (26) can be computed efficiently with the
cumulative distribution function of a multivariate Gaussian, which is also
available in most standard statistical packages, such as Python’s Scipy.

Interestingly, these predictions generalize those obtained in VGPMIL-PR [18].
Indeed, if we do not consider correlation among instances in the test bag,
i.e. ${\boldsymbol{\Sigma}}_{\star}={\mathbf{I}}$, then eqs. (25) and (26)
match those in [18] (last two equations before Section 3.5). Finally, although
here we have detailed how to make predictions for a complete previously unseen
bag ${\mathbf{X}}^{\star}$, the same process can be applied to make
predictions on previously unseen individual instances ${\mathbf{x}}^{\star}$
(patches).

##  3 Experiments

In this section we thoroughly evaluate VGPMIL-PR-I in a real-world problem of
prostate cancer detection. The experimental framework, including data, metrics
and baselines, is explained in Section 3.1. The results are discussed in
Section 3.2. Finally, in Section 3.3 we evaluate our method in a much larger
prostate cancer detection dataset: the well-known PANDA challenge.

|  |  |
---|---|---|---
Cancerous | Non-cancerous
Figure 3: Two examples of cancerous (left) and non-cancerous (right) patches
in the SICAPv2 test set.

###  3.1 Experimental framework

Data description. In this paper we focus on the problem of prostate cancer
detection. However, notice that the algorithm can be applied for any other
type of cancer (and more generally, for any other type of image). Prostate
cancer is the most commonly occurring cancer in men, and the second most
commonly occurring cancer overall, according to the latest 2020 statistics on
age-standardized incidence rate from the World Health Organisation (WHO)
Global Cancer Observatory [31]. We will use the prostate cancer database
presented in [32], which is called SICAPv2 and is publicly available. Although
this database includes information on the Gleason score, which is used to
evaluate the severity of the disease, in this work we will focus on the binary
task of presence/absence of cancer.

We use the original partition of the dataset, which contains 95 training and
31 test WSIs, respectively. These very large images are split in 512x512
patches, resulting in a total amount of 15132 patches for training and 5246
for testing. Following the MIL paradigm, for the training set we only use
binary labels benign/malign at the level of images (bags), but we do not have
information at the level of patches (instances). In order to evaluate the
predictive performance at instance-level, we do have labels for the patches in
the test set. The amount of cancerous (resp. non-cancerous) images for the
train set is 70 (resp. 25). For the test set, it is 25 (resp. 6). For
illustration purposes, a couple of cancerous and non-cancerous patches are
shown in Figure 3. In order to train our model, each patch is represented
through a 128-dimensional feature vector extracted in previous work [7].

Baselines and metrics. Since our model is framed in the field of probabilistic
GP-based MIL methods, we compare with the two most popular approaches VGPMIL
[17] and VGPMIL-PR [18], which were reviewed in Section 2.2. For a fair
comparison, the parameters used for the baselines are the same as those used
for our method (described in next paragraph). For those parameters that do not
have an analogous in our method (e.g. the initialization of
${\mathrm{q}}({\mathbf{y}})$ in VGPMIL), we use the default values proposed in
the original papers. To evaluate the performance of the compared methods we
use four metrics: accuracy, precision, recall and F1-score (which provides a
trade-off between precision and recall). For all the metrics, we use the
standard implementations in the popular Python scikit-learn library [33].

Experimental details. For the underlying GPs, in this work we use the well-
known squared exponential kernel [34], i.e.
$\kappa({\mathbf{x}},{\mathbf{y}})=\gamma\cdot\exp\left(-||{\mathbf{x}}-{\mathbf{y}}||^{2}/(2\ell^{2})\right)$.
Following [17] and [18], we use standard values for the kernel
hyperparameters, i.e. $\gamma=1$ and $\ell$ equals the square root of the
number of features of ${\mathbf{x}},{\mathbf{y}}$ (in this work we set
$\ell=11\approx\sqrt{128}$). The number of inducing points is set to $M=200$,
and their locations are initialized through K-means clustering as in previous
work [17, 18] (namely, $100$ of them are obtained by doing clustering on the
patches that belong to the positive images, and the other $100$ on the patches
that belong to the negative ones). The number of iterations is set to $T=200$,
which was enough to achieve convergence in practice. The expectation of the
posterior distribution $\mathbb{E}_{{\mathrm{q}}({\mathbf{m}})}({\mathbf{m}})$
is initialized with a standard Gaussian for each instance independently.
Notice that the initialization of ${\mathrm{q}}({\mathbf{u}})$ is irrelevant
since it gets updated first in Algorithm 1. As for the value of $\lambda$,
which regulates the strength of the Ising correlation (recall eq. (2.3)), we
will analyze five different values in the experiments,
$\lambda\in\\{0.1,0.5,1.0,5.0,10.0\\}$. This will allow us to empirically
illustrate the effect of $\lambda$.

###  3.2 Experimental results

| $\lambda$ | Accuracy | Precision | Recall | F1-score
---|---|---|---|---|---
VGPMIL | - |  92.22$\pm$0.00 |  96.40$\pm$0.00 |  92.29$\pm$0.00 |  94.30$\pm$0.00
VGPMIL-PR | - |  92.38$\pm$0.03 |  96.44$\pm$0.06 |  92.48$\pm$0.06 |  94.42$\pm$0.02
VGPMIL-PR-I | 0.1 |  92.94$\pm$0.05 |  96.24$\pm$0.15 |  93.52$\pm$0.12 |  94.86$\pm$0.04
0.5 | 93.85$\pm$0.04 | 97.17$\pm$0.19 | 93.90$\pm$0.14 | 95.51$\pm$0.02
1.0 | 94.58$\pm$0.03 | 97.32$\pm$0.06 | 94.83$\pm$0.04 | 96.06$\pm$0.02
5.0 | 95.11$\pm$0.06 | 97.74$\pm$0.14 | 95.18$\pm$0.09 | 96.44$\pm$0.04
10.0 | 95.03$\pm$0.14 | 97.72$\pm$0.16 | 95.09$\pm$0.06 | 96.39$\pm$0.10
Table 1: Predictive performance at the level of patches (instances). In bold, we highlight the values of $\lambda$ for which VGPMIL-PR-I gets better (or equal) performance than both baselines in all the metrics. The results are the mean and standard deviation over five independent runs. | $\lambda$ | Accuracy | Precision | Recall | F1-score
---|---|---|---|---|---
VGPMIL | - |  83.87$\pm$0.00 |  83.33$\pm$0.00 |  100.00$\pm$0.00 |  90.91$\pm$0.00
VGPMIL-PR | - |  90.32$\pm$0.00 |  89.29$\pm$0.00 |  100.00$\pm$0.00 |  94.34$\pm$0.00
VGPMIL-PR-I | 0.1 | 93.55$\pm$0.00 | 92.59$\pm$0.00 | 100.00$\pm$0.00 | 96.15$\pm$0.00
0.5 | 93.55$\pm$0.00 | 92.59$\pm$0.00 | 100.00$\pm$0.00 | 96.15$\pm$0.00
1.0 |  90.32$\pm$0.00 |  92.31$\pm$0.00 |  96.00$\pm$0.00 |  94.12$\pm$0.00
5.0 |  87.10$\pm$0.00 |  95.65$\pm$0.00 |  88.00$\pm$0.00 |  91.67$\pm$0.00
10.0 |  83.87$\pm$0.00 |  95.45$\pm$0.00 |  84.00$\pm$0.00 |  89.36$\pm$0.00
Table 2: Predictive performance at the level of images (bags). In bold, we
highlight the values of $\lambda$ for which VGPMIL-PR-I gets better (or equal)
performance than both baselines in all the metrics. The results are the mean
and standard deviation over five independent runs.

|  | VGPMIL | VGPMIL-PR
---|---|---|---
|  | Neg. | Pos. | Neg. | Pos.
Actual | Neg. | 1 | 5 | 3 | 3
Pos. | 0 | 25 | 0 | 25

|  | VGPMIL-PR-I
---|---|---
|  | $\lambda=0.1$ | $\lambda=0.5$ | $\lambda=1$ | $\lambda=5$ | $\lambda=10$
|  | Neg. | Pos. | Neg. | Pos. | Neg. | Pos. | Neg. | Pos. | Neg. | Pos.
Actual | Neg. | 4 | 2 | 4 | 2 | 4 | 2 | 5 | 1 | 5 | 1
Pos. | 0 | 25 | 0 | 25 | 1 | 24 | 3 | 22 | 4 | 21

Table 3: Confusion matrices obtained at the level of images for the compared
methods.

Figure 4: Patch-level predictions inside each one of the six negative (non-cancerous) WSIs in the test set. Each column is an image (the header is the image identifier in the SICAPv2 dataset). The rows refers to the three compared methods. Each subplot has red/green axis depending on whether the image is correctly classified or not by that method. The blue bars inside the subplots represent the probability of cancer for the different patches inside the image (for ease of visualization, they are sorted in increasing order).  |
---|---
(a) Original image | (b) VGPMIL predictions
|
(c) VGPMIL-PR predictions | (d) VGPMIL-PR-I predictions
Figure 5: Patch level predictions obtained by the compared methods for image 16B0028138, which is non-cancerous. The original image is shown in (a). For predictions (b)–(d), the brightness of the patch is proportional to the probability of cancer (the brighter, the more probability). Figure 6: Variability in the patch-level predictions inside images. In VGPMIL-PR-I, the variability decreases as the strength of the Ising terms, given by $\lambda$, increases. VGPMIL-PR, which does not include Ising terms (i.e. $\lambda=0$), gets larger variability. As explained in the text, the variability inside an image is measured as the standard deviation of the probability of cancer for all the patches inside that image. | $\lambda$ | Training time | Testing time
---|---|---|---
VGPMIL | - | 15.45$\pm$0.53 | 3.14$\pm$0.16
VGPMIL-PR | - | 11.78$\pm$0.75 | 2.63$\pm$0.12
VGPMIL-PR-I | 0.1 | 12.21$\pm$0.26 | 2.23$\pm$0.09
0.5 | 12.18$\pm$0.57 | 2.21$\pm$0.07
1.0 | 11.79$\pm$0.44 | 2.24$\pm$0.15
5.0 | 11.80$\pm$0.74 | 2.22$\pm$0.08
10.0 | 11.64$\pm$0.91 | 2.26$\pm$0.09
Table 4: Computational cost for training and testing the compared methods (in
seconds). We are using 200 iterations in all cases, recall the experimental
details in Section 3.1. The results are the mean and standard deviation over
five independent runs.

In this section we evaluate the performance of the compared methods on the
aforementioned prostate cancer problem. We analyze eight different research
questions, which are discussed in the following paragraphs.

Predictions at the level of instances (patches). Although they only use bag
labels for training, the compared methods can make predictions at the level of
instances, recall Section 2.5. This is important to determine more precisely
in which region (patch) the cancer is present. Table 1 shows the results when
making predictions at patch level. We observe that VGPMIL-PR-I outperforms
both baselines in all the metrics for four out of five values of $\lambda$
(and for $\lambda=0.1$, the baselines are only better in terms of precision).
We also appreciate that VGPMIL-PR-I results are robust across different runs,
obtaining low values of standard deviation. This stability is important for
real-world applications, where one wants to avoid high sensitivity to random
initializations. Finally, notice that, as argued in [19], we also observe that
VGPMIL-PR (slightly) outperforms VGPMIL.

Predictions at the level of bags (images). Table 2 shows the results when
making predictions at the level of images. We observe that VGPMIL-PR-I
outperforms both baselines in all the metrics when $\lambda\in\\{0.1,0.5\\}$.
However, when $\lambda$ becomes larger, the results of VGPMIL-PR-I get worse.
This fact can be explained theoretically because, whereas having low-to-
moderate correlation among patches can be helpful, having strong ones tends to
make the predictions too homogeneous, damaging the bag-level prediction (which
takes into account the correlation among patches). Indeed, in the next
research question we analyze with greater detail how $\lambda$ is affecting
the predictions on cancerous and non-cancerous images separately, which will
provide additional insights. Finally, similar to the patch-level results, we
observe that VGPMIL-PR obtains better results than VGPMIL, as expected. Also,
the predictions at the level of images are very stable across runs (notice the
zero standard deviation).

Analyzing the confusion matrices at bag level. Here we analyze more in detail
the results presented in the previous paragraph (i.e. at the level of images).
Table 3 shows the confusion matrices for the compared methods. Notice that
both baselines classify all the 25 positive (cancerous) images correctly.
However, the difficulties arise at the non-cancerous images. This happens
because of the nature of the MIL problem: as soon as a few patches obtain a
non-negligible probability of cancer, the image will be likely predicted as
cancerous (recall that the MIL formulation establishes that a bag has positive
class as soon as one instance inside the bag has positive class).

Interestingly, the Ising model can help to avoid isolated false positive
predictions on the patch-level (which lead to false positive bag predictions)
by using the patch-level correlation. This is reflected in Table 3, where we
observe that increasingly more negative images are predicted correctly as
$\lambda$ gets higher. In contrast, notice that strong correlation damage the
performance in the positive class, as they penalize the appearance of positive
patches (which would break the homogeneity of the bag, where most patches do
not contain cancer). Therefore, we conclude that instance correlation is
beneficial when used with a low-to-medium intensity. For instance, in this
application $\lambda=0.5$ is the best performing value, and it will be the one
used by default in the sequel.

Analyzing the instance-level predictions for negative bags. In the previous
paragraph, we have explained that negative images are incorrectly classified
because a few patches inside them get classified as positive. Here we provide
a visualization to support this. Figure 4 shows how the patch-level
predictions are distributed inside the six negative images available in the
test set. We observe that the amount of patches with a non-negligible
probability of cancer gets reduced as we move from VGPMIL to VGPMIL-PR, and
then to VGPMIL-PR-I. This translates into better performance at bag-level
(observe that the amount of green-axis subplots increases in the same sequence
VGPMIL $\rightarrow$ VGPMIL-PR $\rightarrow$ VGPMIL-PR-I). Notice that the
improvement from VGPMIL to VGPMIL-PR is larger than from VGPMIL-PR to VGPMIL-
PR-I. This may be due to the simplification that was introduced when computing
the expectation of the truncated multivariate Gaussian in VGPMIL-PR-I, recall
the third-to-last paragraph in Section 2.4.

Visualizing the predictions. In the last paragraph, we have analyzed how the
patch-level predictions are distributed inside non-cancerous images
_quantitatively_. Indeed, Figure 4 represents each patch through a bar.
However, this hampers the _qualitative_ visualization from a medical
viewpoint. Here we focus on such qualitative assessment by visualizing the
predictions obtained for image 16B0028138, see Figure 5. We have chosen this
image because it illustrates best the effect of the coupling term. Notice
that, thanks to these terms, the proposed VGPMIL-PR-I manages to keep all
patches with a probability closer to zero than VGPMIL and VGPMIL-PR. As a
consequence, VGPMIL-PR-I is the only method that correctly classifies this
image as non-cancerous, recall Figure 4. For the other two methods, there are
some patches that trigger the image prediction to be cancerous.

An explicit analysis on the role of $\lambda$. The hyperparameter $\lambda$ is
at the core of the novel VGPMIL-PR-I. It was introduced in the probabilistic
model to regulate the strength of the coupling term, recall eq. (2.3). This
role has been confirmed _indirectly_ in Table 3: when $\lambda$ gets higher,
the predictive performance improves for negative bags and degrades for
positive ones. This can be explained because a higher $\lambda$ homogenizes
the patch-level predictions and difficulties the appearance of positive
patches. Here we perform a more _direct_ measure to gain insights into the
role of $\lambda$. Specifically, we define the “variability inside a bag” as
the standard deviation of the probability of cancer for all the patches inside
that bag. Therefore, this metric measures the dispersion in the patch-level
predictions obtained inside a bag. Figure 6 shows the evolution of this metric
for VPGMIL-PR-I as $\lambda$ increases. As theoretically expected, the metric
decreases as $\lambda$ gets higher. Also, notice that the metric value for
VGPMIL-PR is higher. This is explained because VGPMIL-PR does not incorporate
Ising correlation, i.e. $\lambda=0$. The value for VPGMIL is even higher,
$0.30$, and it is not included in Figure 6 for ease of visualization. This
greater value is probably due to the additional approximations that VGPMIL
involves, which deepens the independence among patches.

Computational cost. Finally, we report the computational training and testing
time for the compared methods, see Table 4. The results are in the same order
of magnitude in all cases, which justifies the practical utility of the novel
VGPMIL-PR-I, which obtained better predictive performance, recall Tables 1 and
2. In fact, VGPMIL-PR-I is slightly faster than VGPMIL, since the Jaakola
bound approximation leveraged in the latter introduces additional parameters
$\xi$ to be estimated. As theoretically expected, the computational cost of
VGPMIL-PR-I and VGPMIL-PR is analogous, since the update equations for the
former are just a generalization of those for the latter, recall Sections 2.4
and 2.5. Finally, notice that the value of $\lambda$ does not affect the
computational cost of VGPMIL-PR-I, as $\lambda$ only regulates the intensity
of the Ising terms (but it does not introduce any additional computation).

Comparison to other related MIL approaches. So far we have focused on the
comparison of VGPMIL-PR-I with VGPMIL and VGPMIL-PR. Since VGPMIL-PR-I builds
on the same type of GP-based modeling, this is the most meaningful comparison
in order to evaluate our main contribution (the Ising term to account for
correlations among patches). However, to provide a wider perspective, it is
interesting to compare the novel VGPMIL-PR-I to other state-of-the-art and
popular families of MIL methods. We consider three families: attention-based
methods, where the two algorithms proposed in [35] are the most popular
approaches; MIL methods based on pseudo-labels such as the recent [7]; and
classical pooling/aggregation methods such us the mean aggregation [36]. These
will be referred to as Att-MIL, Gated-Att-MIL, PS-MIL and Mean-Agg,
respectively.

Let us discuss the results both at instance (patch) and bag (image) levels.
For the former, notice that the formulation of attention-based methods (Att-
MIL and Gated-Att-MIL) and classical aggregation methods (Mean-Agg) do not
allow for making predictions at instance level in a natural way. Namely, the
instance-level labels are not modelled explicitly in this type of methods, and
this is precisely one of their main limitations [37]. Compared to PS-MIL,
which does model instance labels explicitly, the novel VGPMIL-PR-I achieves
higher predictive performance (95.11 vs 85.01 in accuracy and 96.44 vs 88.05
in F1-Score). Regarding bag-level performance, the results are shown in Table
5. We observe that the best results are obtained by attention-based methods
Att-MIL and Gated-Att-MIL, followed by the novel VGPMIL-PR-I.

In conclusion, we observe that the results are quite different depending on
the nature of the model and the user requirements. If one is interested in
predictions at patch level, then the novel VGPMIL-PR-I is the best choice.
However, if one is only interested in image level performance, then attention-
based approaches are the best option for this data. Indeed, we hypothesise
that the performance of attention-based approaches could be even enhanced by
leveraging correlation among patches in a similar way to VGPMIL-PR-I. This is
a interesting line of future research, see Section 4.

| Accuracy | Precision | Recall | F1-score
---|---|---|---|---
VGPMIL-PR-I |  93.55$\pm$0.00 |  92.59$\pm$0.00 |  100.00$\pm$0.00 |  96.15$\pm$0.00
Att-MIL |  96.80$\pm$0.00 |  96.20$\pm$0.00 |  100.00$\pm$0.00 |  98.00$\pm$0.00
Gated-Att-MIL |  96.80$\pm$ 0.00 |  96.20$\pm$0.00 |  100.00$\pm$0.00 |  98.00$\pm$0.00
Mean-Agg |  87.10$\pm$0.00 |  95.70$\pm$0.00 |  88.00$\pm$0.00 |  91.70$\pm$0.00
PS-MIL |  90.32$\pm$NA |  89.28$\pm$NA |  100.00$\pm$NA |  94.33$\pm$NA
Table 5: Comparison with other related MIL methods which are not based on the
GP modeling. The predictive performance at the level of images (bags) is
shown. The results are the mean and standard deviation over five independent
runs. The algorithm PS-MIL was run only once because of its high computational
training cost.

###  3.3 Evaluation on a larger dataset: PANDA

The SICAPv2 dataset used so far is of medium-size (total amount of 126 WSIs,
leading to 20378 patches; recall Section 3.1). This has allowed us to carry
out a very detailed analysis of the results. In this section we show that the
novel VGPMIL-PR-I also performs well on larger datasets, such us the well-
known PANDA set. Although scalability is not an issue from a theoretical
perspective, since the model is based on sparse GPs, it is important to verify
it in practice.

PANDA also tackles the problem of prostate cancer detection, and was presented
at the MICCAI 2020 conference as a challenge222https://panda.grand-
challenge.org/. Since the test set of PANDA is not publicly available, we use
the train/test split proposed in [38], where each split follows the overall
class proportions. Namely, the dataset used here features a total amount of
10503 WSIs, which leads to 1107931 patches. Notice that this is much larger
than SICAPv2 (83 times larger in terms of WSIs).

Table 6 shows the predictive performance at image level, an aspect where
attention-based methods stood out in the previous dataset. In this case, we
observe that VGPMIL-PR-I obtains consistently better results. Additionally, as
outlined in the last research question in Section 3.2, VGPMIL-PR-I is able to
provide instance-level predictions, which is not the case for attention-based
models. We conclude that, for the PANDA dataset, the proposed method is the
best choice in comparison to the other tested approaches.

| Accuracy | Precision | Recall | F1-score
---|---|---|---|---
VGPMIL |  74.87$\pm$0.08 |  74.25$\pm$0.06 | 99.77$\pm$0.00 |  85.13$\pm$0.04
VGPMIL-PR |  90.64$\pm$0.06 |  90.97$\pm$0.05 |  96.60$\pm$0.04 |  93.70$\pm$0.04
VGPMIL-PR-I | 92.57$\pm$0.14 | 95.42$\pm$0.24 |  94.22$\pm$0.11 | 94.82$\pm$0.09
Att-MIL |  90.92$\pm$0.01 |  94.17$\pm$0.01 |  93.43$\pm$0.01 |  93.79$\pm$0.01
Gated-Att-MIL |  91.70$\pm$0.01 |  94.31$\pm$0.01 |  94.40$\pm$0.01 |  94.34$\pm$0.01
Mean-Agg |  88.09$\pm$0.00 |  91.57$\pm$0.01 |  92.27$\pm$0.01 |  91.91$\pm$0.00
PS-MIL |  88.36$\pm$NA |  87.99$\pm$NA |  97.11$\pm$NA |  92.33$\pm$NA
Table 6: Predictive performance at the level of images (bags) in the PANDA
dataset. The results are the mean and standard deviation over five independent
runs. The algorithm PS-MIL was run only once because of its high computational
training cost.

##  4 Conclusions, limitations and future work

In this work we have introduced VGPMIL-PR-I, a novel MIL methodology that
incorporates instance label correlation through a coupling term inspired by
the Ising model. VGPMIL-PR-I is a generalization of another probabilistic MIL
method, whose formulation is theoretically recovered when the influence of the
Ising term converges to zero. In the experimental section, we have shown that
VGPMIL-PR-I outperforms other related state-of-the-art probabilistic MIL
approaches in two real-world problems of prostate cancer detection,
effectively reducing false positive bag predictions and providing instance-
level predictions. We have also provided different visualizations to better
understand the behavior of the proposed model, specially the influence of the
new coupling term.

As discussed along the paper, our model presents several limitations which we
summarize next. Firstly, we needed to introduce a diagonal approximation to
compute the expectation of the truncated mulivariate Gaussian in VGPMIL-PR-I,
recall Section 2.4. This is probably reflected in the empirical performance,
as the improvement when moving from VGPMIL to VGPMIL-PR is generally larger
than when moving from VGPMIL-PR to VGPMIL-PR-I. Secondly, we have observed
that the behaviour of VGPMIL-PR-I depends on the value of $\lambda$, which
regulates the strength of the coupling term. Although we have discussed the
role of $\lambda$ and tested different values, it remains a hyperparameter
that has to be found empirically using the validation set. We believe that its
value (or even distribution over it) could be estimated from the data by
introducing $\lambda$ in the probabilistic modeling. Even more, $\lambda$
could be estimated per image, since the level of correlation could be image-
dependent. Thirdly, we have observed that the image-level performance of
VGPMIL-PR-I is not generally better than that of attention-based methods. This
is probably due to the different nature of the models. Indeed, the explicit
modeling of instance label in GP-based models, which allows them to provide
instance-level predictions, may come at the cost of less accurate bag-level
predictions.

In addition to the aforementioned ideas, this work opens other future research
lines. First, seeing the performance boost obtained in GP-based methods
through the novel coupling term, and taking into account the good results of
attention-based methods in bag-level prediction, it is very interesting to
explore the modeling of instance label correlations in the context of
attention-based methods. Second, notice that we are using mean-field
variational inference to estimate the model parameters in VGPMIL-PR-I. A
promising alternative is to estimate them by directly optimizing the evidence
lower bound (ELBO). Finally, although we have focused on modeling correlation
between neighboring patches in histopathogical images, we expect that the
ideas behind our proposal can boost further research in MIL, by exploiting the
particular structure of the data used in different applications.

## Acknowledgements

This work has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska Curie grant agreement No
860627 (CLARIFY Project), from the Spanish Ministry of Science and Innovation
under project PID2019-105142RB-C22, and by the University of Granada and
FEDER/Junta de Andalucía under project B-TIC-324-UGR20 (Proyectos de I+D+i en
el marco del Programa Operativo FEDER Andalucía). PMA has been supported by
the Margarita Salas postdoctoral fellowship (Spanish Ministry of Universities
with Next Generation EU funds).

## References

  * [1] M.-A. Carbonneau, V. Cheplygina, E. Granger, G. Gagnon, Multiple instance learning: A survey of problem characteristics and applications, Pattern Recognition 77 (2018) 329–353.
  * [2] M. Cui, D. Y. Zhang, Artificial intelligence and computational pathology, Laboratory Investigation 101 (4) (2021) 412–422.
  * [3] S. Huang, Z. Liu, W. Jin, Y. Mu, Bag dissimilarity regularized multi-instance learning, Pattern Recognition 126 (2022) 108583.
  * [4] M. K. K. Niazi, A. V. Parwani, M. N. Gurcan, Digital pathology and artificial intelligence, The lancet oncology 20 (5) (2019) e253–e261.
  * [5] K. P. Murphy, Machine learning: a probabilistic perspective, MIT press, 2012.
  * [6] M. López-Pérez, M. Amgad, P. Morales-Álvarez, P. Ruiz, L. A. Cooper, R. Molina, A. K. Katsaggelos, Learning from crowds in digital pathology using scalable variational gaussian processes, Scientific reports 11 (1) (2021) 1–9.
  * [7] A. Schmidt, J. Silva-Rodríguez, R. Molina, V. Naranjo, Efficient cancer classification by coupling semi supervised and multiple instance learning, IEEE Access 10 (2022) 9763–9773.
  * [8] G. Campanella, M. G. Hanna, L. Geneslaw, A. Miraflor, V. Werneck Krauss Silva, K. J. Busam, E. Brogi, V. E. Reuter, D. S. Klimstra, T. J. Fuchs, Clinical-grade computational pathology using weakly supervised deep learning on whole slide images, Nature medicine 25 (8) (2019) 1301–1309.
  * [9] A. Schmidt, J. Silva-Rodríguez, R. Molina, V. Naranjo, Coupling semi-supervised and multiple instance learning for histopathological image classification, IEEE Access (2022).
  * [10] Y. Wu, A. Schmidt, E. Hernández-Sánchez, R. Molina, A. K. Katsaggelos, Combining attention-based multiple instance learning and gaussian processes for ct hemorrhage detection, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2021, pp. 582–591.
  * [11] M. Kandemir, F. A. Hamprecht, et al., Instance label prediction by dirichlet process multiple instance learning., in: UAI, 2014, pp. 380–389.
  * [12] J. Read, L. Martino, J. Hollmén, Multi-label methods for prediction with sequential data, Pattern Recognition 63 (2017) 45–55.
  * [13] J. Read, L. Martino, D. Luengo, Efficient monte carlo methods for multi-dimensional learning with classifier chains, Pattern Recognition 47 (3) (2014) 1535–1546.
  * [14] J. Read, L. Martino, Probabilistic regressor chains with monte carlo methods, Neurocomputing 413 (2019) 471–486.
  * [15] M. Kim, F. De la Torre, Gaussian processes multiple instance learning, in: ICML, 2010.
  * [16] M. Kandemir, M. Haußmann, F. Diego, K. T. Rajamani, J. Van Der Laak, F. A. Hamprecht, Variational weakly supervised gaussian processes, in: E. R. H. Richard C. Wilson, W. A. P. Smith (Eds.), Proceedings of the British Machine Vision Conference (BMVC), BMVA Press, 2016, pp. 71.1–71.12.
  * [17] M. Haußmann, F. A. Hamprecht, M. Kandemir, Variational bayesian multiple instance learning with gaussian processes, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 6570–6579.
  * [18] F. Wang, A. Pinar, The multiple instance learning gaussian process probit model, in: International Conference on Artificial Intelligence and Statistics, PMLR, 2021, pp. 3034–3042.
  * [19] Z. Wang, A. Lamb, E. Saveliev, P. Cameron, J. Zaykov, J. M. Hernandez-Lobato, R. E. Turner, R. G. Baraniuk, C. Barton, S. P. Jones, et al., Results and insights from diagnostic questions: The neurips 2020 education challenge, in: NeurIPS 2020 Competition and Demonstration Track, PMLR, 2021, pp. 191–205.
  * [20] N. Ding, J. Deng, K. P. Murphy, H. Neven, Probabilistic label relation graphs with ising models, in: Proceedings of the IEEE International Conference on Computer Vision, 2015, pp. 1161–1169.
  * [21] P. Qin, J. Zhao, A polynomial-time algorithm for image segmentation using ising models, in: 2011 Seventh International Conference on Natural Computation, Vol. 2, 2011, pp. 932–935.
  * [22] F. Wang, A. Pinar, The multiple instance learning gaussian process probit model, in: A. Banerjee, K. Fukumizu (Eds.), Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Vol. 130 of Proceedings of Machine Learning Research, PMLR, 2021, pp. 3034–3042.
  * [23] E. Snelson, Z. Ghahramani, Sparse gaussian processes using pseudo-inputs, in: Advances in Neural Information Processing Systems, Vol. 18, 2006.
  * [24] M. Girolami, S. Rogers, Variational bayesian multinomial probit regression with gaussian process priors, Neural Computation 18 (8) (2006) 1790–1817.
  * [25] B. D. Ripley, Spatial statistics, John Wiley & Sons, 2005.
  * [26] J. Hensman, A. De G. Matthews, Z. Ghahramani, Scalable variational Gaussian process classification, in: International conference on artificial intelligence and statistics, 2015, pp. 351–360.
  * [27] C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics), Springer-Verlag, Berlin, Heidelberg, 2006.
  * [28] S. Wilhelm, B. Manjunath, tmvtnorm: A package for the truncated multivariate normal distribution, R Journal 2 (2) (2010) 1–25.
  * [29] Y. Li, S. K. Ghosh, Efficient sampling methods for truncated multivariate normal and student-t distributions subject to linear inequality constraints, Journal of Statistical Theory and Practice 9 (4) (2015) 712–732.
  * [30] N. L. Johnson, S. Kotz, N. Balakrishnan, Continuous univariate distributions, volume 2, Vol. 289, John wiley & sons, 1995.
  * [31] Global cancer observatory, world health organisation, https://gco.iarc.fr/, accessed: 2010-03-15.
  * [32] J. Silva-Rodríguez, A. Colomer, M. A. Sales, R. Molina, V. Naranjo, Going deeper through the gleason scoring scale: An automatic end-to-end system for histology prostate grading and cribriform pattern detection, Computer Methods and Programs in Biomedicine 195 (2020) 105637.
  * [33] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
  * [34] C. Williams, C. Rasmussen, Gaussian processes for machine learning, Vol. 2, MIT press Cambridge, MA, 2006.
  * [35] M. Ilse, J. Tomczak, M. Welling, Attention-based deep multiple instance learning, in: International conference on machine learning, PMLR, 2018, pp. 2127–2136.
  * [36] W. Zhu, Q. Lou, Y. S. Vang, X. Xie, Deep multi-instance networks with sparse label assignment for whole mammogram classification, in: International conference on medical image computing and computer-assisted intervention, Springer, 2017, pp. 603–611.
  * [37] A. Schmidt, P. Morales-Álvarez, R. Molina, Probabilistic attention based on gaussian processes for deep multiple instance learning, IEEE Transactions on Neural Networks and Learning Systems (2023) 1–14doi:10.1109/TNNLS.2023.3245329.
  * [38] J. Silva-Rodriguez, A. Colomer, J. Dolz, V. Naranjo, Self-learning for weakly supervised gleason grading of local patterns, IEEE journal of biomedical and health informatics 25 (8) (2021) 3094–3104.