# Identifying Hosts of Families of Viruses:
A Machine Learning Approach

Anil Raj  Department of Applied Physics and Applied Mathematics
Columbia University, New York<EMAIL_ADDRESS>Michael Dewar  Department
of Applied Physics and Applied Mathematics
Columbia University, New York  Gustavo Palacios  Center for Infection and
Immunity
Mailman School of Public Health
Columbia University, New York  Raul Rabadan  Department of Biomedical
Informatics
Center for Computational Biology and Bioinformatics
Columbia University College of Physicians and Surgeons, New York  Chris H.
Wiggins  Department of Applied Physics and Applied Mathematics
Center for Computational Biology and Bioinformatics
Columbia University, New York

###### Abstract

## Abstract

Identifying viral pathogens and characterizing their transmission is essential
to developing effective public health measures in response to a pandemic.
Phylogenetics, though currently the most popular tool used to characterize the
likely host of a virus, can be ambiguous when studying species very distant to
known species and when there is very little reliable sequence information
available in the early stages of the pandemic. Motivated by an existing
framework for representing biological sequence information, we learn sparse,
tree-structured models, built from decision rules based on subsequences, to
predict viral hosts from protein sequence data using popular discriminative
machine learning tools. Furthermore, the predictive motifs robustly selected
by the learning algorithm are found to show strong host-specificity and occur
in highly conserved regions of the viral proteome.

viral host, machine learning, adaboost, alternating decision tree, mismatch
k-mers

##  1 Introduction

Emerging pathogens constitute a continuous threat to our society, as it is
notoriously difficult to perform a realistic assessment of optimal public
health measures when little information on the pathogen is available. Recent
outbreaks include the West Nile virus in New York (1999); SARS coronavirus in
Hong Kong (2002-2003); LUJO virus in Lusaka (2008); H1N1 influenza pandemic
virus in Mexico and the US (2009); and cholera in Haiti (2010). In all these
cases, an outbreak of unusual clinical diagnoses triggered a rapid response,
and an essential part of this response is the accurate identification and
characterization of the pathogen.

Sequencing is becoming the most common and reliable technique to identify
novel organisms. For instance, LUJO was identified as a novel, very distinct
virus after the sequence of its genome was compared to other arenaviruses
Briese2009 . The genome of an organism is a unique fingerprint that reveals
many of its properties and past history. For instance, arenaviruses are
zoonotic agents usually transmitted from rodents.

Another promising area of research is metagenomics, in which DNA and RNA
samples from different environments are sequenced using shotgun approaches.
Metagenomics is providing an unbiased understanding of the different species
that inhabit a particular niche. Examples include the human microbiome and
virome, and the Ocean metagenomics collection Williamson2008 . It has been
estimated that there are more than 600 bacterial species living in the mouth
but that only 20% have been characterized.

Pathogen identification and metagenomic analysis point to an extremely rich
diversity of unknown species, where partial genomic sequence is the only
information available. The main aim of this work is to develop approaches that
can help infer characteristics of an organism from subsequences of its genomic
sequence where primary sequence information analysis does not allow us to
identify its origin. In particular, our work will focus on predicting the host
of a virus from the viral genome.

The most common approach to deduce a likely host of a virus from the viral
genome is sequence / phylogenetic similarity (i.e., the most likely host of a
particular virus is the one that is infected by related viral species).
However, similarity measures based on genomic / protein sequence or protein
structure could be misleading when dealing with species very distant to known,
annotated species. Other approaches are based on the fact that viruses undergo
mutational and evolutionary pressures from the host. For instance, viruses
could adapt their codon bias for a more efficient interaction with the host
translational machinery or they could be under pressure of deaminating enzymes
(e.g. APOBEC3G or HIV infection). All these factors imprint characteristic
signatures in the viral genome. Several techniques have been developed to
extract these patterns (e.g., nucleotide and dinucleotide compositional
biases, and frequency analysis techniques Touchon2008 ). Although most of
these techniques could reveal an underlying biological mechanism, they lack
sufficient accuracy to provide reliable assessments. A relatively similar
approach to the one presented here is DNA barcoding. Genetic barcoding
identifies conserved genomic structures that contain the necessary information
for classification.

Using contemporary machine learning techinques, we present an approach to
prediciting the hosts of unseen viruses, based on the amino acid sequences of
proteins of viruses whose hosts are well known. Using sequence and host
information of known viruses, we learn a multi-class classifier composed of
simple sequence-motif based questions (e.g., does the viral sequence contain
the motif ‘DALMWLPD’?) that achieves high prediction accuracies on held-out
data. Prediction accuracy of the classifier is measured by the area under the
ROC curve, and is compared to a straightforward nearest-neighbour classifier.
Importantly (and quite surprisingly), a post-processing study of the highly
predictive sequence-motifs selected by the algorithm identifies strongly
conserved regions of the viral genome, facilitating biological interpretation.

##  2 Methods

Our overall aim is to discover aspects of the relationship between a virus and
its host. Our approach is to develop a model that is able to predict the host
of a virus given its sequence; those features of the sequence that prove most
useful are then assumed to have a special biological significance. Hence, an
ideal model is one that is parsimonious and easy to interpret, whilst
incorporating combinations of biologically relevant features. In addition, the
interpretability of the results is improved if we have a simple learning
algorithm which can be straightforwardly verified.

Formally, for a given virus family, we learn a function
$g:\mathcal{S}\rightarrow\mathcal{H}$, where $\mathcal{S}$ is the space of
viral sequences and $\mathcal{H}$ is the space of viral hosts. The space of
viral sequences $\mathcal{S}$ is generated by an alphabet $\mathcal{A}$ where,
$|\mathcal{A}|=4$ (genome sequence) or $|\mathcal{A}|=20$ (primary protein
sequence).

Defining a function on a sequence requires representation of the sequence in
some feature space. Below, we specify a representation
$\phi:\mathcal{S}\rightarrow\mathcal{X}$, where a sequence $s\in\mathcal{S}$
is mapped to a vector of counts of subsequences
$x\in\mathcal{X}\subset\mathbb{N}_{0}^{D}$. Given this representation, we have
the well-posed problem of finding a function
${f}:\mathcal{X}\rightarrow\mathcal{H}$ built from a space of simple binary-
valued functions.

###  2.1 Collected Data

The collected data consist of $N$ genome sequences or primary protein
sequences, denoted $s_{1}\ldots s_{N}$, of viruses whose host class, denoted
$h_{1}\ldots h_{N}$ is known. For example, these could be ‘plant’,
‘vertebrate’ and ‘invertebrate’. The label for each virus is represented
numerically as $\mathbf{y}\in\mathcal{Y}=\\{0,1\\}^{L}$ where
$[\mathbf{y}]_{l}=1$ if the index of the host class of the virus is $l$, and
where $L$ denotes the number of host classes. Note that this representation
allows for a virus to have multiple host classes. Here and below we use
boldface variables to indicate vectors and square brackets to denote the
selection of a specific element in the vector, e.g., $[\mathbf{y}_{n}]_{l}$ is
the $l^{\mathrm{th}}$ element of the $n^{\mathrm{th}}$ label vector.

###  2.2 Mismatch Feature Space

A possible feature space representation of a viral sequence is the vector of
counts of exact matches of all possible $k$-length subsequences ($k$-mers).
However, due to the high mutation rate of viral genomes Duffy2008 ; Pybus2009
, a predictive function learned using this simple representation of counts
would fail to generalize well to new viruses. Instead, motivated by Leslie2004
, we count not just the presence of an individual $k$-mer but also the
presence of subsequences within $m$ mismatches from that $k$-mer. The
mismatch- or $m$-neighborhood of a $k$-mer $\alpha$, denoted
$\mathcal{N}^{m}_{\alpha}$, is the set of all $k$-mers with a Hamming distance
Hamming1950  at most $m$ from it, as shown in Table 1. Let
$\delta_{\mathcal{N}^{m}_{\alpha}}$ denote the indicator function of the
$m$-neigbourhood of $\alpha$ such that

$\delta_{\mathcal{N}^{m}_{\alpha}}(\beta)=\left\\{\begin{array}[]{rl}1&\mathrm{if}~{}\beta\in\mathcal{N}^{m}_{\alpha}\\\
0&\mathrm{otherwise}.\end{array}\right.$ (1)

Table 1: Mismatch feature space representation. The mismatch feature space representation of a segment of a protein sequence …AQGPRIYDDTCQHPSWWMNFEYRGSP… $m=0$ | $m=1$ | $m=2$
---|---|---
kmer | count | kmer | count | kmer | count
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
DQGPS | 0 | DQGPS | 0 | DQGPS | 1
CQGPS | 0 | CQGPS | 1 | CQGPS | 1
CQHPS | 1 | CQHPS | 1 | CQHPS | 1
CQIPS | 0 | CQIPS | 1 | CQIPS | 1
DQIPS | 0 | DQIPS | 0 | DQIPS | 1
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
APGPQ | 0 | APGPQ | 0 | APGPQ | 1
AQGPQ | 0 | AQGPQ | 1 | AQGPQ | 1
AQGPR | 1 | AQGPR | 1 | AQGPR | 1
AQGPS | 0 | AQGPS | 1 | AQGPS | 1
ASGPS | 0 | ASGPS | 0 | ASGPS | 1
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
ARGMP | 0 | ARGMP | 0 | ARGMP | 1
ARGSP | 0 | ARGSP | 1 | ARGSP | 1
YRGSP | 1 | YRGSP | 1 | YRGSP | 1
WRGSP | 0 | WRGSP | 1 | WRGSP | 1
WRGNP | 0 | WRGNP | 0 | WRGNP | 1
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$

We can then define, for any possible $k$-mer $\beta$, the mapping $\phi$ from
the sequence $s$ onto a count of the elements in $\beta$’s $m$-neighbourhood
as

$\phi_{k,m}(s,\beta)=\sum_{\begin{subarray}{c}\alpha\in s\\\
|\alpha|=k\end{subarray}}\delta_{\mathcal{N}^{m}_{\alpha}}(\beta).$ (2)

Finally, the $d^{\mathrm{th}}$ element of the feature vector for a given
sequence is then defined elementwise as

$[\mathbf{x}]_{d}=\phi_{k,m}(s,\beta_{d})$ (3)

for every possible $k$-mer $\beta_{d}\in\mathcal{A}^{k}$, where $d=1\dots D$
and $D=|\mathcal{A}^{k}|$.

Note that when $m=0$, $\phi_{k,0}$ exactly captures the simple count
representation described earlier. This biologically realistic relaxation
allows us to learn discriminative functions that better capture rapidly
mutating and yet functionally conserved regions in the viral genome
facilitating generalization to new viruses.

###  2.3 Alternating Decision Trees

Given this representation of the data, we aim to learn a discriminative
function that maps features $\mathbf{x}$ onto host class labels $\mathbf{y}$,
given some training data
$\\{(\mathbf{x}_{1},\mathbf{y}_{1}),\dots,(\mathbf{x}_{N},\mathbf{y}_{N})\\}$.
We want the discriminative function to output a measure of “confidence”
Schapire1999  in addition to a predicted host class label. To this end, we
learn on a class of functions
$\mathbf{f}:\mathcal{X}\rightarrow\mathbb{R}^{L}$, where the indices of
positive elements of $\mathbf{f}(\mathbf{x})$ can be interpreted as the
predicted labels to be assigned to $\mathbf{x}$ and the magnitudes of these
elements to be the confidence in the predictions.

A simple class of such real-valued discriminative functions can be constructed
from the linear combination of simple binary-valued functions
$\psi:\mathcal{X}\rightarrow\\{0,1\\}$. The functions $\psi$ can, in general,
be a combination of single-feature decision rules or their negations:

$\displaystyle\mathbf{f}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\sum_{p=1}^{P}\mathbf{a}_{p}\psi_{p}(\mathbf{x})$ (4)
$\displaystyle\psi_{p}(\mathbf{x})$ $\displaystyle=$ $\displaystyle\prod_{d\in
S_{p}}\mathbb{I}(x_{d}\geq\theta_{d})$ (5)

where $\mathbf{a}_{p}\in\mathbb{R}^{L}$, $P$ is the number of binary-valued
functions, $\mathbb{I}(\cdot)$ is 1 if its argument is true, and zero
otherwise, $\theta\in\\{0,1,\dots,\Theta\\}$, where
$\Theta=\max_{d,n}[\mathbf{x}_{n}]_{d}$, and $S_{p}$ is a subset of feature
indices. This formulation allows functions to be constructed using
combinations of simple rules. For example, we could define a function $\psi$
as the following

$\psi(\mathbf{x})=\mathbb{I}(x_{5}\geq 2)\times\lnot\mathbb{I}(x_{11}\geq
1)\times\mathbb{I}(x_{1}\geq 4)$ (6)

where $\lnot\mathbb{I}(\cdot)=1-\mathbb{I}(\cdot)$.

Alternatively, we can view each function $\psi_{p}$ to be parameterized by a
vector of thresholds $\boldsymbol{\theta}_{p}\in\\{0,1,\dots,\Theta\\}^{D}$,
where $[\boldsymbol{\theta}_{p}]_{d}=0$ indicates $\psi_{p}$ is not a function
of the $d^{\mathrm{th}}$ feature $[\mathbf{x}]_{d}$. In addition, following
Busa2009 , we can decompose the weights
$\mathbf{a}_{p}=\alpha_{p}\mathbf{v}_{p}$ into a vote vector
$\mathbf{v}\in\\{+1,-1\\}^{L}$ and a scalar weight $\alpha\in\mathbb{R}_{+}$.
The discriminative model, then, can be written as

$\displaystyle\mathbf{f}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\sum_{p=1}^{P}\alpha_{p}\mathbf{v}_{p}\psi_{\boldsymbol{\theta}_{p}}(\mathbf{x}),$
(7) $\displaystyle\psi(\mathbf{x};\boldsymbol{\theta}_{p})$ $\displaystyle=$
$\displaystyle\prod_{d=1}^{D}\mathbb{I}(x_{d}\geq[\boldsymbol{\theta}_{p}]_{d}).$
(8)

Figure 1: Alternating Decision Tree. An example of an ADT where rectangles are
decision nodes, circles are output nodes and, in each decision node,
$[\beta]=\phi_{k,m}(s,\beta)$ is the feature associated with the $k$-mer
$\beta$ in sequence $s$. The output nodes connected to each decision node are
associated with a pair of binary-valued functions $(\psi,\tilde{\psi})$. The
binary-valued function corresponding to the highlighted path is given as
$\tilde{\psi}(\mathbf{x};\boldsymbol{\theta}_{3})=\mathbb{I}([\mathrm{AKNELSID}]\geq
2)\times\lnot\mathbb{I}([\mathrm{AAALASTM}]\geq 1)$ and the associated
$\tilde{\alpha}=0.3$.

Every function in this class of models can be concisely represented as an
Alternating Decision Tree (ADT) Freund1999 . Similar to ordinary decision
trees, ADTs have two kinds of nodes: decision nodes and output nodes. Every
decision node is associated with a single-feature decision rule, the
attributes of the node being the relevant feature and corresponding threshold.
Each decision node is connected to two output nodes corresponding to the
associated decision rule and its negation. Thus, binary-valued functions in
the model come in pairs $(\psi,\tilde{\psi})$; each pair is associated with
the the pair of output nodes for a given decision node in the tree (see Figure
1). Note that $\psi$ and $\tilde{\psi}$ share the same threshold vector
$\boldsymbol{\theta}$ and only differ in whether they contain the associated
decision rule or its negation. The attributes of the output node pair are the
vote vectors $(\mathbf{v},\tilde{\mathbf{v}})$ and the scalar weights
$(\alpha,\tilde{\alpha})$ associated with the corresponding functions
$(\psi,\tilde{\psi})$.

Each function $\psi$ has a one-to-one correspondence with a path from the root
node to its associated output node in the tree; the single-feature decision
rules in $\psi$ being the same as those rules associated with decision nodes
in the path, with negations applied appropriately. Combinatorial features can,
thus, be incorporated into the model by allowing for trees of depth greater
than 1. Including a new function $\psi$ in the model is, then, equivalent to
either adding a new path of decision and output nodes at the root node in the
tree or growing an existing path at one of the existing output nodes. This
tree-structured representation of the model will play an important role in
specifying how Adaboost, the learning algorithm, greedily searches over an
exponentially large space of binary-valued functions. It is important to note
that, unlike ordinary decision trees, each example runs down an ADT through
every path originating from the root node.

###  2.4 Multi-class Adaboost

Having specified a representation for the data and the model, we now briefly
describe Adaboost, a large-margin supervised learning algorithm which we use
to learn an ADT given a data set. Ideally, a supervised learning algorithm
learns a discriminative function $\mathbf{f}^{*}(\mathbf{x})$ that minimizes
the number of mistakes on the training data, known as the Hamming loss
Hamming1950 :

$\mathbf{f}^{*}(\mathbf{x})=\arg\min_{\mathbf{f}}\mathcal{L}_{h}(\mathbf{f})=\sum_{\begin{subarray}{c}1\leq
n\leq N\\\ 1\leq l\leq
L\end{subarray}}\mathbb{I}\left(H([\mathbf{f}(\mathbf{x}_{n})]_{l})\neq[\mathbf{y}_{n}]_{l}\right)$
(9)

where $H(.)$ denotes the Heaviside function. The Hamming loss, however, is
discontinuous and non-convex, making optimization intractable for large-scale
problems.

Adaboost is the unconstrained minimization of the exponential loss, a smooth,
convex upper-bound to the Hamming loss, using a coordinate descent algorithm.

$\tilde{\mathbf{f}}^{*}(\mathbf{x})=\arg\min_{\mathbf{f}}\mathcal{L}_{e}(\mathbf{f})=\sum_{\begin{subarray}{c}1\leq
n\leq N\\\ 1\leq l\leq
L\end{subarray}}\exp\left(-[\mathbf{y}_{n}]_{l}[f_{l}(\mathbf{x}_{n})]_{l}\right).$
(10)

Adaboost learns a discriminative function $\mathbf{f}(\mathbf{x})$ by
iteratively selecting the $\psi$ that maximally decreases the exponential
loss. Since each $\psi$ is parameterized by a $D$-dimensional vector of
thresholds $\boldsymbol{\theta}$, the space of functions $\psi$ is of size
$O((\Theta+1)^{D})$, where $\Theta$ is the largest $k$-mer count observed in
the data, making an exhaustive search at each iteration intractable for high-
dimensional problems.

To avoid this problem, at each iteration, we only allow the ADT to grow by
adding one decision node to one of the existing output nodes. To formalize
this, let us define $Z(\boldsymbol{\theta})=\\{d:[\boldsymbol{\theta}]_{d}\neq
0\\}$ to be the set of active features corresponding to a function $\psi$. At
the $t^{\mathrm{th}}$ iteration of boosting, the search space of possible
threshold vectors is then given as
$\\{\boldsymbol{\theta}:\exists\tau<t,Z(\boldsymbol{\theta})\supset
Z(\boldsymbol{\theta}_{\tau}),|Z(\boldsymbol{\theta})|-|Z(\boldsymbol{\theta}_{\tau})|=1\\}$.
In this case, the search space of thresholds at the $t^{\mathrm{th}}$
iteration is of size $O(t\Theta D)$ and grows linearly in a greedy fashion at
each iteration (see Figure 1). Note, however, that this greedy search,
enforced to make the algorithm tractable, is not relevant when the class of
models are constrained to belong to ADTs of depth 1.

In order to pick the best function $\psi$, we need to compute the decrease in
exponential loss admitted by each function in the search space, given the
model at the current iteration. Formally, given the model at the
$t^{\mathrm{th}}$ iteration, denoted $\mathbf{f}^{t}(\mathbf{x})$, the
exponential loss upon inclusion of a new decision node, and hence the creation
of two new paths
$(\psi_{\boldsymbol{\theta}},\tilde{\psi}_{\boldsymbol{\theta}})$, into the
model can be written as

$\displaystyle\mathcal{L}_{e}(\mathbf{f}^{t+1})$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{l=1}^{L}\exp\left(-[\mathbf{y}_{n}]_{l}[\mathbf{f}^{t}(\mathbf{x}_{n})+\alpha\mathbf{v}\psi_{\boldsymbol{\theta}}(\mathbf{x}_{n})+\tilde{\alpha}\tilde{\mathbf{v}}\tilde{\psi}_{\boldsymbol{\theta}}(\mathbf{x}_{n})]_{l}\right)$
(11) $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{l=1}^{L}w^{t}_{nl}\exp\left(-[\mathbf{y}_{n}]_{l}[\alpha\mathbf{v}\psi_{\boldsymbol{\theta}}(\mathbf{x}_{n})+\tilde{\alpha}\tilde{\mathbf{v}}\tilde{\psi}_{\boldsymbol{\theta}}(x_{n})]_{l}\right)$
(12)

where
$w^{t}_{nl}=\exp\left(-[\mathbf{y}_{n}]_{l}[\mathbf{f}^{t}(\mathbf{x}_{n})]_{l}\right)$.
Here $w_{nl}^{t}$ is interpreted as the weight on each sample, for each label,
at boosting round $t$. If, at boosting round $t-1$, the model disagrees with
the true label $l$ for sample $n$, then $w_{nl}^{t}$ is large. If the model
agrees with the label then the weight is small. This ensures that the boosting
algorithm chooses a decision rule at round $t$, preferentially discriminating
those examples with a large weight, as this will lead to the largest reduction
in $\mathcal{L}_{e}$.

For every possible new decision node that can be introduced to the tree,
Adaboost finds the ($\alpha$,$\mathbf{v}$) pair that minimizes the exponential
loss on the training data. These optima can be derived as

$[\mathbf{v}^{*}]_{l}=\begin{cases}1&\text{if
}\omega^{t}_{+,l}\geq\omega^{t}_{-,l}\\\ -1&\text{otherwise}\end{cases}$ (13)
$\alpha^{*}=\frac{1}{2}\ln\frac{W^{t}_{+}}{W^{t}_{-}}$ (14)

where for each new path $\psi_{n}$ associated with each new decision node

$\displaystyle\omega^{t}_{\pm,l}$ $\displaystyle=$
$\displaystyle\sum_{n:\psi_{n}y_{nl}=\pm 1}w^{t}_{nl}$ (15) $\displaystyle
W^{t}_{\pm}$ $\displaystyle=$ $\displaystyle\sum_{n,l:v_{l}\psi_{n}y_{nl}=\pm
1}w^{t}_{nl}.$ (16)

Corresponding equations for the ($\tilde{\alpha}$,$\tilde{\mathbf{v}}$) pair
can be written in terms of $\tilde{W}^{t}_{\pm,l}$ and $\tilde{W}^{t}_{\pm}$
obtained by replacing $\psi_{n}$ with $\tilde{\psi}_{n}$ in the equations
above. The minimum loss function for the threshold $\boldsymbol{\theta}$ is
then given as

$\mathcal{L}_{e}(\mathbf{f}^{t+1})=2\sqrt{W^{t}_{+}W^{t}_{-}}+2\sqrt{\tilde{W}^{t}_{+}\tilde{W}^{t}_{-}}+W^{t}_{o}$
(17)

where $W^{t}_{o}=\sum_{n,l:\psi_{n}=\tilde{\psi}_{n}=0}w^{t}_{nl}$. Based on
these model update equations, each iteration of the Adaboost algorithm
involves building the set of possible binary-valued functions to search over,
selecting the one for which the loss function given by Eq. 17 and computing
the associated $(\alpha,\mathbf{v})$ pair using Eq. 13 and Eq. 14.

##  3 Results

###  3.1 Data specifications

We aim to learn a predictive model to identify hosts of viruses belonging to a
specific family; we show results for _Picornaviridae_ and _Rhabdoviridae_.
_Picornaviridae_ is a family of viruses that contain a single stranded,
positive sense RNA. The viral genome usually contains about 1-2 Open Reading
Frames (ORF), each coding for protein sequences about 2000-3000 amino acids
long. _Rhabdoviridae_ is a family of negative sense single stranded RNA
viruses whose genomes typically code for five different proteins: large
protein (L), nucleoprotein (N), phosphoprotein (P), glycoprotein (G), and
matrix protein (M). The data consist of 148 viruses in the _Picornaviridae_
family and 50 viruses in the _Rhabdoviridae_ family. For some choice of $k$
and $m$, we represent each virus as a vector of counts of all possible
$k$-mers, up to $m$-mismatches, generated from the amino-acid alphabet. Each
virus is also assigned a label depending on its host: vertebrate /
invertebrate / plant in the case of _Picornaviridae_ , and animal / plant in
the case of _Rhabdoviridae_ (see Table S1 for the names and label assignments
of viruses). Using multiclass Adaboost, we learn an ADT classifier on training
data drawn from the set of labeled viruses and test the model on the held-out
viruses.

###  3.2 BLAST Classifier accuracy

Given whole protein sequences, a straightforward classifier is given by a
nearest-neighbour approach based on the Basic Local Alignment Search Tool
(BLAST) Altschul1990 . We can use BLAST score (or $P$-value) as a measure of
the distances between the unknown virus and a set of viruses with known hosts.
The nearest neighbor approach to classification then assigns the host of the
closest virus to the unknown virus. Intuitively, as this approach uses the
whole protein to perform the classification, we expect the accuracy to be very
high. This is indeed the case – BLAST, along with a $1$-nearest neighbor
classifier, successfully classifies all viruses in the _Rhabdoviridae_ family,
and all but 3 viruses in the _Picornaviridae_ family. What is missing from
this approach, however, is the ability to ascertain and interpret host
relevant motifs.

###  3.3 ADT Classifier accuracy

The accuracy of the ADT model, at each round of boosting, is evaluated using a
multi-class extension of the Area Under the Curve (AUC). Here the ‘curve’ is
the Receiver Operating Characteristic (ROC) which traces a measure of the
classification accuracy of the ADT for each value of a real-valued
discrimination threshold. As this threshold is varied, a virus is considered a
true (or false) positive if the prediction of the ADT model for the true class
of that protein is greater (or less) than the threshold value. The ROC curve
is then traced out in True Positive Rate – False Positive Rate space by
changing the threshold value and the AUC score is defined as the area under
this ROC curve.

The ADT is trained using 10-fold cross validation, calculating the AUC, at
each round of boosting, for each fold using the held-out data. The mean AUC
and standard deviation over all folds is plotted against boosting round in
Figures 2, 3. Note that the ‘smoothing effect’ introduced by using the
mismatch feature space allows for improved prediction accuracy for $m>0$. For
_Picornaviridae_ , the best accuracy is achieved at $m=5$, for a choice of
$k=12$; this degree of ‘smoothing’ is optimal for the algorithm to capture
predictive amino-acid subsequences present, up to a certain mismatch, in
rapidly mutating viral protein sequences. For _Rhabdoviridae_ , near perfect
accuracy is achieved with merely one decision rule, i.e., _Rhabdoviridae_ with
plant or animal hosts can be distinguished based on the presence or absence of
one highly conserved region in the L protein.

Figure 2: Prediction accuracy for _Picornaviridae_. A plot of (a) mean AUC vs
boosting round, and (b) 95% confidence interval vs boosting round. The mean
and standard deviation were computed over 10-folds of held-out data, for
_Picornaviridae_ , where $k=12$. Figure 3: Prediction accuracy for
_Rhabdoviridae_. A plot of (a) mean AUC vs boosting round, and (b) 95%
confidence interval vs boosting round. The mean and standard deviation were
computed over 5-folds of held-out data, for _Rhabdoviridae_ , where $k=10$.
The relatively higher uncertainty for this virus family was likely due to very
small sample sizes. Note that the cyan curve lies on top of the red curve.

###  3.4 Predictive subsequences are conserved within hosts

Having learned a highly predictive model, we would like to locate where the
selected $k$-mers occur in the viral proteomes. We visualize the $k$-mer
subsequences selected in a specific ADT by indicating elements of the mismatch
neighborhood of each selected subsequence on the virus protein sequences. In
Figure 4, the virus proteomes are grouped vertically by their label with their
lengths scaled to $[0,1]$. Quite surprisingly, the predictive $k$-mers occur
in regions that are strongly conserved among viruses sharing a specific host.
Note that the representation we used for viral sequences retained no
information regarding the location of each $k$-mer on the virus protein.
Furthermore, these selected $k$-mers are significant as they are robustly
selected by Adaboost for different choices of train / test split of the data,
as shown in Figure 5.

Figure 4: Visualizing predictive subsequences. A visualization of the mismatch
neighborhood of the first 7 $k$-mers selected in an ADT for _Picornaviridae_ ,
where $k=12,m=5$. The virus proteomes are grouped vertically by their label
with their lengths scaled to $[0,1]$. Regions containing elements of the
mismatch neighborhood of each $k$-mer are then indicated on the virus
proteome. Note that the proteomes are not aligned along the selected $k$-mers
but merely stacked vertically with their lengths normalized. Figure 5:
Visualizing predictive regions of protein sequences. A visualization of the
mismatch neighborhood of the first 7 $k$-mers, selected in all ADTs over
10-fold cross validation, for _Picornaviridae_ , where $k=12,m=5$. Regions
containing elements of the mismatch neighborhood of each selected $k$-mer are
indicated on the virus proteome, with the grayscale intensity on the plot
being inversely proportional to the number of cross-validation folds in which
some $k$-mer in that region was selected by Adaboost. Thus, darker spots
indicate that some $k$-mer in that part of the proteome was robustly selected
by Adaboost. Furthermore, a vertical cluster of dark spots indicate that
region, selected by Adaboost to be predictive, is also strongly conserved
among viruses sharing a common host type.

##  4 Discussion

We have presented a supervised learning algorithm that learns a model to
classify viruses according to their host and identifies a set of highly
discriminative oligopeptide motifs. As expected, the $k$-mers selected in the
ADT for _Picornaviridae_ (Figure 4, 5) and _Rhabdoviridae_ (Figure S.1, S.2)
are mostly selected in areas corresponding to the replicase motifs of the
polymerase – one of the most conserved parts of the viral genome. Thus, given
that partial genomic sequence is normally the only information available, we
could achieve quicker bioinformatic characterization by focusing on the
selection and amplification of these highly predictive regions of the genome,
instead of full genomic characterization and contiguing. Moreover, in contrast
with generic approaches currently under use, such a targeted amplification
approach might also speed up the process of sample preparation and improve the
sensitivity for viral discovery.

Over representation of highly similar viruses within the data used for
learning is an important source of overfitting that should be kept in mind
when using this technique. Specifically, if the data largely consist of
nearly-similar viral sequences (e.g. different sequence reads from the same
virus), the learned ADT model would overfit to insignificant variations within
the data (even if 10-fold cross validation were employed), making
generalization to new subfamilies of viruses extremely poor. To check for
this, we hold out viruses corresponding to a particular subfamily (see Table
S1 for subfamily annotation of the viruses used), run 10-fold cross validation
on the remaining data and compute the expected fraction of misclassified
viruses in the held-out subfamily, averaged over the learned ADT models. For
_Picornaviridae_ , viruses belonging to the subfamilies _Parechovirus_ (0.47),
_Tremovirus_ (0.8), _Sequivirus_ (0.5), and _Cripavirus_ (1.0) were poorly
classified with misclassification rates indicated in parentheses. Note that
the _Picornaviridae_ data used consist mostly of Cripaviruses; thus, the high
misclassification rate could also be attributed to a significantly lower
sample size used in learning. For _Rhabdoviridae_ , viruses belonging to
_Novirhabdovirus_ (0.75) and _Cytorhabdovirus_ (0.77) were poorly classified.
The poorly classified subfamilies, however, contain a very small number of
viruses, showing that the method is strongly generalizable on average.

Other applications for this technique include identification of novel
pathogens using genomic data, analysis of the most informative fingerprints
that determine host specificity, and classification of metagenomic data using
genomic information. For example, an alternative application of our approach
would be the automatic discovery of multi-locus barcoding genes. Multi-locus
barcoding is the use of a set of genes which are discriminative between
species, in order to identify known specimens and to flag possible new species
Seberg2009 . While we have focused on virus host in this work, ADTs could be
applied straightforwardly to the barcoding problem, replacing the host label
with a species label. Additional constraints on the loss function would have
to be introduced to capture the desire for suitable flanking sites of each
selected $k$-mer in order to develop the universal PCR primers important for a
wide application of the discovered barcode Kress2008 .

## Acknowledgments

The authors thank Vladimir Trifonov and Joseph Chan for interesting
suggestions and discussions.

## References

  * [1] Thomas Briese, Janusz T Paweska, Laura K McMullan, Stephen K Hutchison, Craig Street, Gustavo Palacios, Marina L Khristova, Jacqueline Weyer, Robert Swanepoel, Michael Egholm, Stuart T Nichol, and W Ian Lipkin.  Genetic detection and characterization of Lujo virus, a new hemorrhagic fever-associated arenavirus from southern Africa.  PLoS pathogens, 5(5):e1000455, May 2009.
  * [2] Shannon J Williamson, Douglas B Rusch, Shibu Yooseph, Aaron L Halpern, Karla B Heidelberg, John I Glass, Cynthia Andrews-Pfannkoch, Douglas Fadrosh, Christopher S Miller, Granger Sutton, Marvin Frazier, and J Craig Venter.  The Sorcerer II Global Ocean Sampling Expedition: metagenomic characterization of viruses within aquatic microbial samples.  PloS one, 3(1):e1456, January 2008.
  * [3] Marie Touchon and Eduardo P C Rocha.  From GC skews to wavelets: a gentle guide to the analysis of compositional asymmetries in genomic data.  Biochimie, 90(4):648–59, 2008.
  * [4] Siobain Duffy, Laura A Shackelton, and Edward C Holmes.  Rates of evolutionary change in viruses: patterns and determinants.  Nature Reviews Genetics, 9(4):267–276, 2008.
  * [5] Oliver G Pybus and Andrew Rambaut.  Evolutionary analysis of the dynamics of viral infectious disease.  Nature reviews. Genetics, 10(8):540–50, August 2009.
  * [6] C S Leslie, E Eskin, A Cohen, J Weston, and W S Noble.  Mismatch string kernels for discriminative protein classification.  Bioinformatics, 20(4):467–476, 2004.
  * [7] R W Hamming.  Error detecting and error correcting codes.  Bell System Technical Journal, 29:147—-160, 1950.
  * [8] R E Schapire and Y Singer.  Improved boosting algorithms using confidence-rated predictions.  Machine Learning, 37(3):297–336, 1999.
  * [9] R Busa-Fekete and B Kegl.  Accelerating AdaBoost using UCB.  In JMLR: Workshop and Conference Proceedings, KDD cup 2009, volume 7, pages 111–122, 2009.
  * [10] Y Freund and L Mason.  The Alternating Decision Tree Algorithm.  In Proceedings of the 16th International Conference on Machine Learning, pages 124–133, 1999.
  * [11] Stephen F Altschul, Warren Gish, Webb Miller, Eugene W Myers, and David J Lipman.  Basic Local Alignment Search Tool.  Journal of Molecular Biology, 215:403–410, 1990.
  * [12] Ole Seberg and Gitte Petersen.  How many loci does it take to DNA barcode a crocus?  PloS one, 4(2):e4598, January 2009.
  * [13] W John Kress and David L Erickson.  DNA barcodes: genes, genomics, and bioinformatics.  Proceedings of the National Academy of Sciences of the United States of America, 105(8):2761–2762, 2008.

## Supplementary Figures

Figure S.1: Visualizing predictive subsequences for _Rhabdoviridae_. A
visualization of the mismatch neighborhood of the $k$-mer selected in an ADT
for _Rhabdoviridae_ , where $k=10,m=2$. The virus proteomes are grouped
vertically by their label with their lengths scaled to $[0,1]$. Regions
containing elements of the mismatch neighborhood of each $k$-mer are then
indicated on the virus proteome. Figure S.2: Visualizing predictive regions
for _Rhabdoviridae_. A visualization of the mismatch neighborhood of the
$k$-mers selected in an ADT for _Rhabdoviridae_ , where $k=10,m=2$. The virus
proteomes are grouped vertically by their label with their lengths scaled to
$[0,1]$. Regions containing elements of the mismatch neighborhood of each
$k$-mer are then indicated on the virus proteome.


abstractOblivious transfer (OT) protocol is an essential tool in cryptography that provides a wide range of applications such as secure multi-party computation, private information retrieval, private set intersection, contract signing, and privacy-preserving location-based services. The OT protocol has different variants such as one-out-of-2, one-out-of-n, k-out-of-n, and OT extension. In the OT (one-out-of-2, one-out-of-n, and OT extension) protocol, the sender has a set of messages, whereas the receiver has a key. The receiver sends that key to the sender in a secure way; the sender cannot get any information about the received key. The sender encrypts every message by operating on every message using the received key and sends all the encrypted messages to the receiver. The receiver is able to extract only the required message using his key. However, in the k-out-of-n OT protocol, the receiver sends a set of k keys to the sender, and in replay, the sender sends all the encrypted messages. The receiver uses his keys and extracts the required messages, but it cannot gain any information about the messages that it has not requested. Generally, the OT protocol requires high communication and computation cost if we transfer millions of oblivious messages. The OT extension protocol provides a solution for this, where the receiver transfers a set of keys to the sender by executing a few numbers of OT protocols. Then, the sender encrypts all the messages using cheap symmetric key cryptography with the help of a received set of keys and transfers millions of oblivious messages to the receiver. In this work, we present different variants of OT protocols such as one-out-of-2, one-out-of-n, k-out-of-n, and OT extension. Furthermore, we cover various aspects of theoretical security guarantees such as semi-honest and malicious adversaries, universally composable, used techniques, computation, and communication efficiency aspects. From the analysis, we found that the semi-honest adversary-based OT protocols required low communication and computation costs as compared to malicious adversary-based OT protocols.
journaltitleACM Computing Surveys
shortjournalACM Comput. Surv.
titleA Survey of Oblivious Transfer Protocol
multiparty computation,oblivious transfer extension protocol,Oblivious transfer protocol,pseudorandom generator,random oracle model

[title=Appendix References,filter=appendixOnlyFilter]

§ PROOFS OF STATEMENTS IN PRELIMINARIES

§ PROOFS OF SECURITY OF THE BIT OT PROTOCOL

§ PROOF OF THE ZKOQS AND $K$-OUT-OF-$N$ STRING OT PROTOCOLS

§ PROOF OF THE COMPOSABILITY OF <CIT.>


# Domain Generalization on Medical Imaging Classification using Episodic
Training with Task Augmentation

Chenxin Li  Qi Qi  Xinghao Ding  Yue Huang  Dong Liang  and Yizhou Yu  Fellow,
IEEE This work was partially supported by National Key Research and
Development Program of China (No.2020YFC2003900) and the National Natural
Science Foundation of China under Grants U19B2031, 61971369, in part by
Fundamental Research Funds for the Central Universities 20720200003.Chenxin
Li, Qi Qi, Xinghao Ding, Yue Huang are with the School of Informatics, Xiamen
University, Xiamen 361005, China. (e-mail:huangyue05@gmail)Dong Liang is with
the Paul C. Lauterbur Research Center for Biomedical Imaging, Shenzhen
Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen
518055, China. (dong.liang@siat.ac.cn)Yizhou Yu are with the Deepwise AI
Laboratory, Beijing 100125, China. (e-mail<EMAIL_ADDRESS>
###### Abstract

Medical imaging datasets usually exhibit domain shift due to the variations of
scanner vendors, imaging protocols, etc. This raises the concern about the
generalization capacity of machine learning models. Domain generalization
(DG), which aims to learn a model from multiple source domains such that it
can be directly generalized to unseen test domains, seems particularly
promising to medical imaging community. To address DG, recent model-agnostic
meta-learning (MAML) has been introduced, which transfers the knowledge from
previous training tasks to facilitate the learning of novel testing tasks.
However, in clinical practice, there are usually only a few annotated source
domains available, which decreases the capacity of training task generation
and thus increases the risk of overfitting to training tasks in the paradigm.
In this paper, we propose a novel DG scheme of episodic training with task
augmentation on medical imaging classification. Based on meta-learning, we
develop the paradigm of episodic training to construct the knowledge transfer
from episodic training-task simulation to the real testing task of DG.
Motivated by the limited number of source domains in real-world medical
deployment, we consider the unique task-level overfitting and we propose task
augmentation to enhance the variety during training task generation to
alleviate it. With the established learning framework, we further exploit a
novel meta-objective to regularize the deep embedding of training domains. To
validate the effectiveness of the proposed method, we perform experiments on
histopathological images and abdominal CT images.

Medical imaging classification, domain generalization, model-agnostic meta-
learning, episodic training, task-level overfitting

##  1 Introduction

Medical imaging analysis usually provides crucial clinical information for
diagnosis. For instance, in oncology studies, the classification of epithelium
and stroma for histopathological images is a prior step to calculate
epithelium-stroma (ES) ratio. Then this ratio can be used to perform patient
stratification in many different tumor types, especially breast cancer [1].
Another example is abdominal computed tomography (CT) images. The accurate
segmentation of objective organ from the surrounding regions is crucial for
the radiation therapy planning and follow-up evaluation [2].

In clinical practice, manual annotations are tedious, error-prone and time-
consuming. They are even more difficult to acquire in some conditions where
specialized biomedical expert knowledge is required. The computer-assisted
automatic algorithms are thus developed, including the classical machine
learning methods [3] and recent deep convolutional neural networks (CNN) [4].
However, these methods are based on IID assumption, i.e., training and testing
set are independent and identically distributed [5], and they would suffer
performance degradation when IID is not satisfied [6]. Moreover, this non-IID
scenario (also called as domain shift problem) is common in clinical practice,
where the medical images are usually acquired by different institutions with
various types of scanner vendors, patient populations, disease severity, etc.
To cope with it, a straightforward way is to fine-tune [7], but it requires
extra annotated data from the test domain. Without requiring the annotation on
test domain, unsupervised domain adaptation (UDA) [8, 9, 10, 11] seems more
appealing. However, it still needs to access the test-domain data during the
training procedure. Then, for every new test domain, we have to repeatedly
collect the training images from the test domain, which limits real-world
medical deployment.

Figure 1: Pipeline of domain generalization (DG) on medical images. Trained on
multiple annotated source domains
$\\{\mathcal{D}_{1},\mathcal{D}_{2},...,\mathcal{D}_{K}\\}$, the model can
directly generalize to an unseen test domain $\mathcal{D}_{K+1}$, i.e.,
without accessing any information of $\mathcal{D}_{K+1}$ during the training
phase.

In contrast, it has a broad prospect for clinical practice to train a
generalizable model that can be directly applied to a new domain (cf. Figure
1), i.e., domain generalization (DG) [12, 13, 14, 15, 16, 17]. Without
accessing the information of test domain, the principle of DG is to capture
the representation that is insensitive to the changes of domain-specific
statistics. To achieve this goal, some works explore data augmentation or
self-supervised learning signals [13, 14] to promote the model to learn a
domain-invariant embedding. However, these designs are somehow heuristic and
highly rely on the characteristics of tasks. The others aim to align the
feature embedding of multiple source domains [12, 18], which has shown
effectiveness on domain adaptation (DA). However, embedding alignment between
source domains does not necessarily lead to the alignment between source and
unseen test domains.

Recently, model-agnostic meta-learning (MAML) [19], which aims to train the
model that can be fast adapted to new tasks, has been promising in many
applications like few-shot learning. Then this spirit is introduced to address
DG with a paradigm of episodic training [15, 16, 17]. A key difference between
episodic training and conventional supervised training is that the data point
in episodic training is a task. Specifically, a meta-task (virtual DG task) is
generated episodically, with the domain shift simulation from splitting the
available source domains into meta-train (virtual training) and meta-test
(virtual testing) domains. This scheme implies the knowledge transfer of
coping with domain shift from virtual training tasks to real testing tasks,
based on the task-level similarity between them.

However, since the test domain is unseen during training, the real domain
shift is also unobserved. Thus it inevitably exists a shift between the
training task with simulated domain shift and the testing task with real
domain shift. Furthermore, in clinical practice, it is common that there are
only several labeled source domains available. This makes the meta-task
generation as well as domain shift simulation restricted into just a few
patterns, which further limits the capacity of training task generation. Thus,
the risk of overfitting to training tasks in the episodic paradigm is
increased, under which the model would suffer the performance drop of DG from
simulated training task to real testing task.

In this work, we aim to alleviate this special risk of task-level overfitting.
We propose task augmentation to enhance the variety for the meta-task
generation, which promotes the model to extract the task-invariant knowledge
against various domain shift simulations. We further incorporate a novel meta-
update objective in the proposed episodic learning framework to perform the
regularization for deep semantic embedding.

Our main contributions are summarized as follows.

  * •

We propose a novel DG method for generalizable medical image classification on
unseen test domains, through the paradigm of episodic training.

  * •

We consider the risk of overfitting to the simulated training tasks in the
existing scheme of episodic training, which is motivated by the limited number
of available annotated source domains in clinical practice.

  * •

We develop the task augmentation to alleviate this unique task-level
overfitting, using our proposed mixed task sampling (MTS) strategy to enhance
the variety for the training task simulation.

  * •

We exploit the proposed learning framework in the deep semantic embedding with
a novel meta-objective, which regularizes the cross-domain alignment from both
sample-wise and prototype-wise.

  * •

We conduct extensive experiments on two typical tasks, including epithelium-
stroma classification on histopathological images, and liver segmentation on
abdominal CT images. The results demonstrate the superiority of our method
over the state-of-the-art DG methods.

##  2 Related Work

In this section, we review medical imaging classification, mainly on the
practice of two typical tasks. Then we survey some methods of domain
adaptation (DA). Finally, we review the recent domain generalization (DG)
methods.

###  2.1 Medical Imaging Classification

Medical imaging classification is critical in many clinical applications [4].
In this paper, we investigate two common tasks, i.e., epithelium-stroma
classification and liver segmentation (which can be treated as pixel-wise
classification). In oncology studies, epithelium and stroma are basic animal
tissue types in histopathological images. The ratio of epithelium-stroma (ES)
can be used to perform patient stratification and follow-up [1], and the
classification of epithelium-stroma is a prior step for calculating it. In the
early machine learning methods, hand-crafted features are extracted according
to the characteristics of histopathological images [20, 21]. Recently, deep
learning approaches achieve remarkable performance [22, 23], and some methods
further investigate the generalization ability of deep models on epithelium-
stroma classification [24, 25].

Moreover, in computer-aided hepatic disease diagnosis and treatment planning,
the accurate segmentation of liver organ is a crucial prerequisite.
Statistical deformable models utilize shape priors, intensity distributions as
well as boundary and region information [26, 27]. Learning-based methods have
also been explored to seek powerful features [28]. Recently, convolutional
neural networks (CNNs) have been widely utilized in the liver segmentation
[29, 30]. But the variety of imaging protocols, scanner etc. further raises
the concern on generalization ability of these models.

###  2.2 Domain Adaptation

Domain adaptation (DA) plays a crucial role in cross-dataset evaluation
especially when the training and testing set have domain shift. Generally, it
is performed by mitigating the distribution mismatch between the source and
target domains[31]. In medical image analysis, many works focus on feature
distribution alignment by minimizing the distance between the feature
distributions of source and target domain [8, 9, 10, 11], using maximum mean
discrepancy, adversarial learning, or regularization term. Furthermore, some
methods consider to conduct the domain alignment and refine the task-specific
decision boundary at the same time [32, 33]. This motivates us to incorporate
the paradigm of episodic learning in the deep embedding alignment from both
sample-wise and prototype-wise perspectives.

###  2.3 Domain Generalization

As a recently emerged concept, domain generalization (DG) aims to learn a
model from multiple distributed source domains that can generalize to unseen
test domains. To achieve this, some methods try to learn domain-invariant
representation by minimizing the domain discrepancy across source domains,
using maximum mean discrepancy (MMD) [34], contrastive loss [35], adversarial
learning [18, 12], autoencoders [36]. The others design heuristic schemes,
such as manipulating deep neural network architectures [37], leveraging self-
supervision signals [13], designing novel training paradigm [38, 39], or
conducting data augmentation [40, 14]. Very recently, with the rise of model-
agnostic meta-learning (MAML) [19], an paradigm of episodic learning is
introduced to address DG [41, 42, 43, 16, 17]. The principle behind this is to
simulate virtual meta-task episodically so that the knowledge of maintaining
generalization under domain shift can be transferred to the real task of DG.
However, there is still a shift between virtual training task and real testing
task, which raises the risk of overfitting to training tasks in the paradigm.

Figure 2: Overview of the proposed domain generalization (DG) method using
episodic training with task augmentation. The meta-task is simulated from
training domains episodically as a virtual training task. The task
augmentation is performed to enhance the variety of meta-task using mixed task
sampling (MTS). The overall proposed training framework is incorporated in the
cross-domain deep embedded alignment from sample-wise $\mathcal{L}_{sa}$ and
prototype222Prototype [44, 45] means the centroid of the embedding for same-
class samples. It is usually used in conjunction with cosine distance to
perform a decision boundary or classifier.-wise $\mathcal{L}_{pa}$.

##  3 Method

In this paper, following the common setting of domain generalization (DG)
[46], we denote input and label space by $X$ and $Y$, and the $K$ source
domains available during training as
$\\{\mathcal{D}_{1},\mathcal{D}_{2},...,\mathcal{D}_{K}\\}$ with different
distributions on the joint space $X\times Y$. $Y$ is the discrete set
$\\{1,2,...C\\}$ shared by all domains and $C$ denotes the number of classes.
Samples are drawn from the domain dataset
$\mathcal{D}_{k}=\\{(x_{n},y_{n})\\}_{n=1}^{N_{k}}$, where $N_{k}$ is the
number of samples in $k$ domain. The goal of DG is to use multiple source
domains $\\{\mathcal{D}_{1},\mathcal{D}_{2},...,\mathcal{D}_{K}\\}$ to learn a
mapping $X\rightarrow Y$ that can generalize well on an unseen test domain
$\mathcal{D}_{K+1}$, i.e., without accessing any information of
$\mathcal{D}_{K+1}$ during the training phase.

The overview of the proposed method is presented in Fig. 2. The remainders of
this section are organized as follows. First, we describe the model-agnostic
episodic training backbone. Then, we introduce task augmentation, which is
proposed to enhance task-level variety via a mixed task sampling (MTS)
strategy. Finally, we present the meta-objective to regularize the proposed
learning procedure in the deep semantic embedding.

###  3.1 Episodic Training Backbone for DG

Generally, we consider a basic model in medical image classification,
consisting of a feature extractor $\mathcal{F}_{\phi}:X\rightarrow Z$ and a
classifier network $\mathcal{C}_{\theta}:Z\rightarrow Y$, where $Z$ is the
lower-dimensional feature space.333For segmentation task, Z is the high-
dimensional space for feature maps. In standard supervised learning, the
parameters $(\phi,\theta)$ are usually optimized on a task-specific loss
$\mathcal{L}_{task}$, e.g. cross-entropy loss. Although the minimization of
$\mathcal{L}_{task}$ may produce highly discriminative features
$z=\mathcal{F}_{\phi}(x)$ as well as a well-behaved classifier on training
domains, nothing in this process emphasizes the generalization performance of
the deep model to unseen test domains. In contrast, based on model-agnostic
meta-learning (MAML) [19], the paradigm of episodic training optimizes the
model according to the performance on the virtual test domain, which reflects
the generalization ability.

In particular, the virtual training task of DG, i.e., meta-task
$\mathcal{T}_{i}$ is episodically sampled from the set
${\mathcal{T}=\\{\mathcal{T}_{1},\mathcal{T}_{2},...,\mathcal{T}_{I}}\\}$, and
the available source domains are split into meta-train
$\mathcal{D}^{\mathcal{T}_{i}}_{tr}$ and meta-test
$\mathcal{D}^{\mathcal{T}_{i}}_{te}$ domains that depends on
$\mathcal{T}_{i}$. Then two loops of optimization are applied on the model
training procedure. The first one is inner optimization, where the model of
feature net $\mathcal{F}_{\phi}$ and classifier $\mathcal{C}_{\theta}$ is
updated with one or more steps of gradient descent on
$\mathcal{D}^{\mathcal{T}_{i}}_{tr}$ with respect to a task loss
$\mathcal{L}_{task}$:

$(\phi^{\prime}_{i},\theta^{\prime}_{i})\leftarrow(\phi,\theta)-\alpha\nabla_{(\phi,\theta)}\mathcal{L}_{task}(D^{\mathcal{T}_{i}}_{tr};\phi,\theta)$
(1)

where $\alpha$ denotes the learning rate w.r.t. the inner optimization. This
leads to a temporary predictive model
$\mathcal{C}_{\theta^{\prime}_{i}}\circ\mathcal{F}_{\phi^{\prime}_{i}}$ with
the high task accuracy on the meta-train domains
$\mathcal{D}^{\mathcal{T}_{i}}_{tr}$. To further evaluate the generalization
performance, the second outer optimization is conducted on the generated
temporary model with weights of $(\phi^{\prime}_{i},\theta^{\prime}_{i})$.
Specifically, besides the task-specific loss, a meta-learning objective is
utilized to emphasize some certain properties that we desire the temporary
predictive model to exhibit. Specially, the meta-objective is calculated on
the temporary model with weights of $(\phi^{\prime}_{i},\theta^{\prime}_{i})$
while the gradients are computed towards the origin model with weights
of$(\phi,\theta)$:

$\displaystyle(\phi,\theta)\leftarrow(\phi,\theta)-\beta\nabla_{(\phi,\theta)}\\{$
$\displaystyle\mathcal{L}_{task}(D^{\mathcal{T}_{i}}_{tr};\phi,\theta)$ (2)
$\displaystyle+\mathcal{L}_{meta}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})\\},$

where $\beta$ is the learning rate w.r.t. outer optimization.

###  3.2 Task Augmentation via Mixed Task Sampling

The paradigm of episodic learning optimizes the model to behave well on the
virtual training task of DG such that it can also behave well on the real
testing task, based on the similarity between training and testing tasks.
However, the risk of overfitting to training tasks emerges due to the
inevitable shift between training and testing tasks. This risk could be
increased in clinical practice where it is often seen that only several
labeled source domains are available, which further limits the capacity of
training task generation. To alleviate this special overfitting, inspired by
the art of data augmentation, we propose the novel task augmentation to
enhance the variety during training meta-task simulation.

Specifically, in the existing paradigm of episodic training, the meta-task is
generated by the task sampling (TS) strategy that $K$ available source domains
are split into $K-1$ meta-train and the held-out meta-test to simulate domain
shift. Intuitively, the number of meta-task patterns is $K$ (i.e.,
$C_{K}^{1}$). This is limited when $K$ is small (e.g., $K=3$), which is a
common phenomenon in real-world medical deployment. In comparison, we perform
the task augmentation to enhance the capacity of meta-task simualtion. We
present the mixed task sampling (MTS) strategy, where the meta-train domains
are generated in the same way while the meta-test domain is acquired by the
interpolation among all the source domains:

$\mathcal{D}_{te}=\bigcup_{k=1}^{K}Sample(\mathcal{D}_{k},\;r_{k}\cdot
N_{te}),\quad s.t.\sum^{K}_{k=1}r_{k}=1$ (3)

where $Sample(\mathcal{D},N)$ is the sampling operator that we define to
randomly select $N$ samples from domain $\mathcal{D}$. $N_{te}$ denotes the
number of samples for $D_{te}$ and $r_{k}$ represents the interpolation ratio
for $k$-$th$ source domain. In particular, the existing dominant TS strategy
is the special case of our MTS approach when $r_{K}=1$. Furthermore, compared
to the previous TS approach, the proposed MTS increases the capacity of meta-
task patterns largely from original $K$, which highlights the variety
enhancement during training task generation.

###  3.3 Meta-Objective to Regularize Deep Semantic Embedding

Here we further discuss the design of meta-learning objective
$\mathcal{L}_{meta}$ in Eq. 2. A straightforward approach is to encourage the
high task accuracy on meta-test domain, i.e.,
$\mathcal{L}_{meta}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})=\mathcal{L}_{task}(D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})$,
which exhibits the generalization performance of the model. Recently, the
alignment of semantic embedding across training domains has been developed in
many tasks like DA and DG [8, 9, 12, 18, 16, 17], which promotes the model to
learn class-discriminative whilst domain-invariant representation. In light of
this, we incorporate the meta-objective of episodic learning in the learned
deep embedding, with the regularization of cross-domain alignment from sample-
wise and prototype-wise, respectively.

Specially, we use sample-wise alignment to encourage clustering the features
of the same class while keeping the features of different class far apart.
Equipped with contrastive learning [47] and cosine distance based decision
boundary [45, 44], the loss of sample-wise alignment is:

$\displaystyle\mathcal{L}_{sa}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})=$
(4)
$\displaystyle\sum_{x_{j}^{c}\in(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te})}\\{(1-\cos(\mathcal{F}(x_{j}^{c}),\mu^{g}_{c}))+\sum_{d=1}^{C}\mathbbm{1}(d\neq
c)\cos(\mathcal{F}(x_{j}^{c}),\mu^{g}_{d})\\}$

where $\mathbbm{1}(\cdot)$ is an indicator function.
$\mu^{g}=[\mu^{g}_{1},\mu^{g}_{2},...,\mu^{g}_{C}]$ denotes the set of domain-
general prototypes. They are actually the weight vectors of classifier
$\mathcal{C}$, as
$\cos(\mathcal{F}(x_{j}^{c}),\mu^{g}_{c})=\frac{\mathcal{F}(x_{j}^{c})}{||\mathcal{F}(x_{j}^{c})||}\cdot\frac{\mu^{g}_{c}}{||\mu^{g}_{c}||}$.
Moreover, we introduce the domain-specific prototypes, which are dependent on
the respective statistics of training domains. Specifically, for domain $k$,
we calculate the centroid of class $c$ as
$\mu^{k}_{c}=\frac{1}{n_{c}^{k}}\sum_{(x_{j}^{k},y_{j}^{k}=c)\in
D_{k}}\mathcal{F}(x_{j}^{k})$, and the prototypes of all classes are
$\mu^{k}=[\mu^{k}_{1},\mu^{k}_{2},...,\mu^{k}_{C}]$. With these prototypes,
the prediction scores of an input sample can be obtained by calculating the
cosine distance between feature embedding and prototypes. In particular, in
the meta-task $\mathcal{T}_{i}$, assume that the meta-train is composed of the
$K-1$ domains $\\{\mathcal{D}_{1},\mathcal{D}_{2},...,\mathcal{D}_{K-1}\\}$ in
all $K$ source domains. Given a query sample $x_{j}$ from meta-test, we
calculate K prediction scores in total based on the $K-1$ domain-specific
prototypes $\\{\mu^{1},\mu^{2},...,\mu^{{K-1}}\\}$ and the domain-general
prototypes $\mu^{g}$. We denote the prediction scores for sample $x_{j}$ as
$\\{P^{1}_{j},P^{2}_{j},...,P^{{K-1}}_{j},P^{g}_{j}\\}$. Then the prototype-
wise alignment is performed to encourage the prediction scores across
different prototypes to be consistent with each other, i.e., the loss of
prototype-wise alignment is:

$\mathcal{L}_{pa}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})=\sum_{a\neq
b}\sum_{x_{j}\in
D_{te}}\frac{1}{2}(KL(P_{j}^{a}||P_{j}^{b})+KL(P_{j}^{b}||P_{j}^{a}))$ (5)

where  $(a,b)\in\\{1,2,...,K-1,g\\}$,
$KL(P_{j}^{a}||P_{j}^{b})=\sum_{c=1}^{C}P^{a}_{jc}*\log{\frac{P^{a}_{jc}}{P^{b}_{kc}}}$.
Finally, the overall meta-objective is used to regularize the proposed
episodic training framework in the deep embedding by integrating the cross-
domain alignment from both sample-wise and prototype-wise:

$\displaystyle\mathcal{L}_{meta}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{tr};\phi^{\prime}_{i},\theta^{\prime}_{i})=\gamma_{1}\mathcal{L}_{fa}+\gamma_{2}\mathcal{L}_{ca}$
(6)

where $\gamma_{1}$ and $\gamma_{2}$ are the trade-off parameters to balance
the two alignment losses. The algorithm is shown in Algorithm 1.

1

Input : Source training domains $\\{D_{k}\\}^{K}_{k=1}$;

Hyperparameters $\alpha,\beta,\gamma_{1},\gamma_{2}>0$;

Maximum number of iterations $I$;

Output : Feature extractor $\mathcal{F}_{\phi}$, classifier network
$\mathcal{C}_{\theta}$,

2 for  _$i=1:I$_ do

3 Generate a meta-task $\mathcal{T}_{i}$, with meta-train $\mathcal{D}_{tr}$
and meta-test $\mathcal{D}_{te}$ sampled from source domains $\mathcal{D}$ via
MTS ;

4 // Eq.3

5
$(\phi^{\prime}_{i},\theta^{\prime}_{i})\leftarrow(\phi,\theta)-\alpha\nabla_{(\phi,\theta)}\mathcal{L}_{task}(D^{\mathcal{T}_{i}}_{tr};\phi,\theta)$
;

6 // Inner optim, Eq.1

7 Compute sample-wise alignment loss:
$\mathcal{L}_{sa}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})$
;

8 // Eq.4

9 Compute prototype-wise alignment loss:
$\mathcal{L}_{pa}(D^{\mathcal{T}_{i}}_{tr},D^{\mathcal{T}_{i}}_{te};\phi^{\prime}_{i},\theta^{\prime}_{i})$
;

10 // Eq.5

11
$\mathcal{L}_{meta}=\gamma_{1}*\mathcal{L}_{sa}+\gamma_{2}*\mathcal{L}_{pa}$ ;

12 // Eq.6

13
$(\phi,\theta)\leftarrow(\phi,\theta)-\beta\nabla_{(\phi,\theta)}(\mathcal{L}_{task}+\mathcal{L}_{meta})$
;

14 // Outer optim, Eq.2

15

16 end for

Algorithm 1 The proposed paradigm of episodic training with task augmentation
for DG. Table 1: Statistics of the datasets used in our experiments.

Histopathological Images | Abdominal CT Scans
---|---
Dataset | Staining method | Organ | Scanner | Location | # Image | Dataset | # Subject | # Image
VGH | H&E | Breast | Aperio | Netherlands | 5920 | BTCV | 30 | 5488
NKI | H&E | Breast | Aperio | Vancouver | 8337 | CHAOS | 20 | 2341
IHC | IHC | Coloretal | Mirax Scan | Finland | 1376 | IRCAD | 20 | 2074
NCH | H&E | Coloretal | Unknown | Germany | 26417 | LITS | 40 | 5785

##  4 Experimental Results

###  4.1 Dataset and Evaluation Metrics

We evaluate the proposed method on epithelium-stroma classification, using the
histopathological images from four common publicly available datasets,
including VGH [48], NKI [48], IHC [49], and NCH [50]. Among them, each image
is labeled as the class of epithelium or stroma. They are collected from
different staining methods, human organs, scanners as well as institutions,
which causes the domain shift among them. The detailed information is shown in
Figure 3 and Table 1. The input images are resize to 224$\times$224, following
previous practice [22, 24, 23, 25].

We also validate our method on the abdominal CT images from four public data
sources, including BTCV [51], CHAOS [52], IRCAD [53], LiTS [54]. Although some
of these datasets are labeled in a multi-organ manner, the organ that is
commonly labeled in the four datasets is only liver. So we conduct DG of liver
segmentation task on these datasets. They exhibit heterogeneous distributions
due to the varying imaging conditions such as different scanners, clinical
centers and subjects. The detailed information of each dataset is shown in
Figure 3 and Table 1. Please note the details of data acquisition for the
abdominal CT datasets are not present because they are not provided in their
official dataset descriptions. Meanwhile, we use only the first forty subjects
in LiTS such that the number of used subjects is in a similar scale compared
to other datasets. For pre-processing, the HU-value of CT scans are truncated
to the range $[-200,250]$ to remove the irrelevant information, and then they
are normalized to $[0,1]$ in intensity value. All the images are resized to
128$\times$128 in the axial plane. Simple data augmentation as random flipping
is applied.

Following the common practice in DG [13, 16, 17], we use leave-one-domain-out
cross-validation, i.e., training on multiple source domains and testing on the
remaining unseen one, with randomly dividing the samples of each domain into
70% training data and 30% testing data. To evaluate the performance on
epithelium-stroma classification, we adopt the accuracy metric as
$ACC=\frac{N_{c}}{N_{t}}$, where $N_{c}$ and $N_{t}$ represents the number of
correctly classified images and total number of images, respectively. To
evaluate the performance on liver segmentation, we adopt the most commonly-
used metric, i.e., Dice coefficient (DC). It is calculated as
$DC=\frac{{2|P\cap G|}}{|P|+|G|}$, where P and G denote the set of pixel-wise
binary prediction and ground-truth, respectively.

Figure 3: Histopathological images and abdominal CT scans from diverse
datasets. The example images and the intensity histograms for every dataset
are present.

###  4.2 Implementation Details

For image classification of epithelium-stroma, we adopt the AlexNet [55] pre-
trained on ImageNet as the model backbone. The feature extractor $\mathcal{F}$
comprises the top layer of the AlexNet up to the pool7 layer and an additional
bottleneck layer fc8 with 512 units. Accordingly, the classifier $\mathcal{C}$
is composed of the fully connected layer with weight size as $512\times 2$.
For the trade-off parameters in $\mathcal{L}_{fa}$ and $\mathcal{L}_{ca}$, we
set $\gamma_{1}=1.0$ and $\gamma_{2}=0.5$ such that it is in a similar scale
to $\mathcal{L}_{task}$. For the outer optimization to obtain
$(\phi^{\prime}_{i},\theta^{\prime}_{i})$, we use Adam optimizer with learning
rate $\beta=5\times 10^{-5}$. For the inner optimization to update
$(\phi,\theta)$, we use plain, non-adaptive gradient descent with learning
rate $\alpha=5\times 10^{-5}$, and clip the gradients by norm (threshold by
2.0) to prevent them from exploding. To construct the meta-task during each
iteration, the samples with mini-batch size = 120 are selected from each meta-
train and meta-test domains. The total number of training iterations is 10000.
Our experiments are implemented in Tensorflow with an Nvidia TITAN Xp 12 GB
GPU. On all the experiments, we report the average and standard deviation over
three independent runs. For liver segmentation, we modify the network backbone
as U-Net and reduce the mini-batch size to 60 due to computational overhead.
The rest settings remain the same.

Table 2: Quantitative results of DG on epithelium-stroma classification
(Acc[%]).

Source | Target | MLDG | Epi-FCR | MetaReg | JiGen | MASF | DeepAll | ETTA-SE
---|---|---|---|---|---|---|---|---
[12] | [41] | [38] | [13] | [16] | (Baseline) | (Ours)
NKI, IHC, NCH | VGH | 91.13$\pm$0.08 | 91.49$\pm$0.13 | 91.74$\pm$0.09 | 92.05$\pm$0.11 | 92.43$\pm$0.13 | 90.43$\pm$0.21 | 93.51$\pm$0.19
VGH, IHC, NCH | NKI | 89.98$\pm$0.05 | 90.15$\pm$0.12 | 90.33$\pm$0.17 | 90.46$\pm$0.09 | 90.32$\pm$0.06 | 89.62$\pm$0.19 | 91.95$\pm$0.17
VGH, NKI, NCH | IHC | 86.68$\pm$0.07 | 87.56$\pm$0.08 | 87.77$\pm$0.19 | 88.29$\pm$0.13 | 88.37$\pm$0.16 | 86.48$\pm$0.08 | 90.32$\pm$0.12
VGH, NKI, IHC | NCH | 87.69$\pm$0.16 | 88.02$\pm$0.31 | 88.57$\pm$0.21 | 89.05$\pm$0.10 | 89.52$\pm$0.06 | 87.17$\pm$0.13 | 91.28$\pm$0.18
Average | 88.85 | 89.31 | 89.60 | 89.96 | 90.16 | 88.43 | 91.77

###  4.3 Experiments on the epithelium-stroma Classification in
Histopathological Images

####  4.3.1 Comparison with Baseline Model

We merge the images from all the source domains and use them to train the DG
baseline in a supervised way, as the DeepAll model. For clarity, the overall
proposed method is called as episodic training with task augmentation to
regularize the deep semantic embedding (ETTA-SE). Table 2 presents the
quantitative results on epithelium-stroma classification by them. It can be
observed that the performance of our method is better than the baseline in all
the settings, raising the accuracy from 88.43% to 91.77% on average. Moreover,
t-SNE [56] is used to analyze the extracted features from $\mathcal{F}_{\phi}$
by DeepAll and our ETTA-SE, as shown in Figure 4. The color type (i.e., blue,
green, gold, red) represents the domain, and the color shade (i.e., dark,
light) denotes the class. It appears that our method produces a better
separation of classes and learns the more domain-invariant embedding, compared
to the DeepAll baseline. In particular, for the unseen test domain (blue
color), some samples of epithelium class (dark blue) are pushed apart from the
cluster of stroma class (light blue) by our method.

(a) DeepAll baseline.

(b) Our proposed ETTA-SE.

Figure 4: The t-SNE visualization of feature embedding learned by (a) DeepAll
baseline and (b) our proposed ETTA-SE. The color type indicates domain
(’blue’-unseen test domain; ’green’,’golden’,’red’-source domains) and color
shade indicates class (’dark’-epithelium class; ’light’-stroma class).

####  4.3.2 Comparison with Other DG methods

We compare with recent state-of-the-art DG methods, which are described as
follows. MLDG [41]: A meta-learning method for domain generalization, which
contains a model-agnostic training procedure to simulates train/test domain
shift during training. Epi-FCR [38]: A single deep network with an episodic
training procedure, which exposes the model to domain shift. MetaReg [42]: A
method that uses episodic training with regularization functions to encode the
notion of domain generalization. JiGen [13]: A model trained by the gradient
signals from the supervised learning for semantic labels and the self-
supervised learning for a puzzle problem. MASF [16]: A recent state-of-the-art
method that utilizes the paradigm of model-agnostic episodic learning for
semantic features. These methods are implemented on the task of epithelium-
stroma classification using the public codes provided by their authors. From
the quantitative results list in Table 2, it is observed that MLDG improves
slightly compared to DeepAll baseline. Epi-FCR and MetaReg further enhance the
performance. The not bad results of MLDG, Epi-FCR, MetaReg indicate that
exposing the model to domain shift simulation episodically benefits the
generalization performance of deep models to unseen domains. JiGen and MASF
improve 1.53% and 1.73% at average accuracy compared to the DeepAll baseline.
The performance achieved by MASF shows the effectiveness of introducing the
episodic training scheme in the learning of semantic features. Finally, the
proposed ETTA-SE, which incorporates the task augmentation as well as deep
embedded regularization, outperforms the state-of-the-art DG methods on Across
all the settings on different target domains, it appears that our method
consistently outperforms the state-of-the-art DG methods.

Table 3: Ablation study on key components of our method on epithelium-stroma
classification (Acc[%]).

ET | TA | $\mathcal{L}_{sa}$ | $\mathcal{L}_{pa}$ | NKI,IHC,NCH-VGH | VGH,IHC,NCH-NKI | VGH,NKI,NCH-IHC | VGH,NKI,IHC-IHC | Average
---|---|---|---|---|---|---|---|---
- | - | - | - | 90.43$\pm$0.21 | 89.62$\pm$0.19 | 86.48$\pm$0.08 | 87.17$\pm$0.13 | 88.43
✓ | - | - | - | 91.27$\pm$0.23 | 90.67$\pm$0.15 | 86.17$\pm$0.25 | 89.81$\pm$0.09 | 89.48
✓ | ✓ | - | - | 91.33$\pm$0.27 | 90.97$\pm$0.18 | 88.76$\pm$0.14 | 90.90$\pm$0.10 | 90.49
✓ | ✓ | ✓ | - | 91.96$\pm$0.24 | 90.93$\pm$0.28 | 89.84$\pm$0.10 | 91.18$\pm$0.19 | 90.98
✓ | ✓ | - | ✓ | 92.78$\pm$0.19 | 91.15$\pm$0.24 | 90.12$\pm$0.13 | 91.16$\pm$0.24 | 91.27
✓ | ✓ | ✓ | ✓ | 93.51$\pm$0.19 | 91.95$\pm$0.17 | 90.32$\pm$0.12 | 91.28$\pm$0.18 | 91.77

####  4.3.3 Ablation Study

We conduct an ablation study to investigate five key issues in our ETTA-SE
framework on the task of epithelium-stroma classification : 1) contribution of
each component, 2) effect of task augmentation (TA) against task-level
overfitting, 3) ablation on the mixed degree of mixed task sampling (MTS)
strategy, 4) impact of the domain alignment of semantic embedding (SE) between
source and unseen domains, 5) stability analysis of our overall model.

Contribution of each component. We validate the effect of the four key
components in our method, i.e., episodic training (ET), task augmentation
(TA), the loss of sample-wise alignment $\mathcal{L}_{sa}$ and prototype-wise
alignment $\mathcal{L}_{pa}$, respectively, by an ablation study shown in
Table 3. The 1st row corresponds to the DeepAll baseline, where the standard
training is performed by merging all the data from source domains. The last
row corresponds to the proposed ETTA-SE. Note that when neither
$\mathcal{L}_{fa}$ nor $\mathcal{L}_{ca}$ is added (i.e., ETTA, without SE), a
task loss calculated on meta-test is used as the basic $\mathcal{L}_{meta}$.
It can be observed that each component plays its own role. Specifically,
compared to DeepAll, ET improves 1.05% at average Acc, showing the effect of
simulating domain shift to solve DG (2nd row). TA further improves 1.01%,
which demonstrates it is necessary to introduce task augmentation (3rd row).
Meanwhile, the two loss terms for regularizing deep embedded alignment are
both beneficial to the model performance. More importantly, we observe that
the performance gain achieved by using them together is greater than the sum
of the gain by using them respectively. Thus, these two perspectives for deep
embedded alignment appear to be complementary to each other.

(a) ET w/o TA

(b) ET w/ TA

Figure 5: The exhibition of the task-level overfitting. The area of shaded
region is calculated to measure the performance gap between training and
testing tasks (i.e., overfitting degree). ’ET’-episodic training, ’TA’-task
augmentation.

Analysis of task augmentation against overfitting. We conduct the ablation
analysis on the paradigm of episodic training (ET) w/ and w/o our task
augmentation (TA). In standard supervised learning, a common approach to
observe overfitting is by comparing the task loss on the training and testing
sets. Accordingly, in episodic training, we observe the loss of DG performance
on the training and testing tasks, respectively. Concretely, we use the task
loss on the meta-test domain as the former and that on the unseen test domain
as the latter. As shown in Figure 5, the model can fit to the training tasks
well during ET whether with TA or not. However, the obvious gap (shaded
region) between the loss curves exhibits the overfitting phenomenon to
training tasks. We further calculate the area of shaded region as an
approximate observation for overfitting degree. When equipping the ET paradigm
with TA, the area decreases from 936.28 to 808.59, which shows that the
mitigation of the task-level overfitting is achieved by the proposed TA.

Figure 6: Ablation on the mixed degree of MTS.

Figure 7: Effect of the domain alignment on SE.

Ablation on the mixed degree of mixed task sampling. We investigate how the
enhanced variety by task augmentation affects the performance of DG, by an
ablation study on the mixed degree of mixed task. Specifically, given $K$
source domains, $K-1$ domains are randomly selected as meta-train, and the
held-out domain denote as (denoted as $\mathcal{D}_{ho}$) as meta-test, using
the task sampling (TS) strategy in existing episodic training paradigm. In
comparison, in the proposed mixed task sampling (MTS), the meta-test domain is
acquired by the interpolation from all K source domain (cf. Eq. 3). To
facilitate illustration, we denote the interpolation ratio of held-out domain
$\mathcal{D}_{ho}$ as $r_{ho}$. As shown in Figure 7, we conduct experiments
of MTS with $r_{ho}$ in different situations, including dynamic values with
the range of $[0.0,0.5]$, $[0.5,0.1]$, $[0.0,1.0]$, and fixed values at $0.0$,
$1.0$. $r_{ho}=0$ corresponds to $\mathcal{D}_{ho}$ excluded in the
interpolation process for meta-test generation, which means the meta-test
domain is just the combination of meta-train domains. It appears to behave not
well due to the signal of domain shift simulation is too weak in such a
situation. $r_{ho}=1$ corresponds to TS and seems to outperforms $r_{ho}=0$,
which shows the effectivenss of episodic simulation for virtual training tasks
of DG. However, in TS, the risk of overfitting to virtual training tasks
prevents the model from achieving the good performance on the real testing DG
task. As a comparison, we observe the proposed MTS outperforms TS under most
variations of $r_{ho}$, which demonstrates the validity of enhancing variety
during the meta-task sampling process. An exception is $r_{ho}\in[0,0.5]$ that
achieves the lower but very color results with TS, which is probably due to
the still weak domain shift simulation. Finally, we see the proposed MTS
yields the best results under $r_{ho}\in[0,1]$.

Semantic embedded alignment between source and unseen domains. The proposed
meta-objective performs the alignment regularization on the semantic embedding
across source domains episodically. Here we analyze the domain alignment of
the deep embedding between source and unseen test domains, which reflects the
effect of meta-objective on the real DG task. As shown in Figure 7, the loss
curves of $L_{sa}$ and $L_{pa}$ are drawn to exhibit the effect of domain
alignment. We see the alignment losses between source domains (light-colored
curves) falls normally during the training phase. At the same time, the losses
between source and unseen domains (dark-colored curves) drops following the
trends of those between source domains. This demonstrates although the
proposed meta-objective for domain alignment is applied on the semantic
embedding between source domains, the regularization effect can be generalized
to the real DG task, which promotes the deep embedded alignment between source
and unseen test domains.

Table 4: Quantitative results of DG on liver segmentation (Dice[%]).

Source | Target | MLDG | Epi-FCR | MetaReg | JiGen | MASF | DeepAll | ETTA-SE
---|---|---|---|---|---|---|---|---
[12] | [41] | [38] | [13] | [16] | (Baseline) | (Ours)
CHAOS, IRCAD, LITS | BTCV | 80.40$\pm$0.24 | 82.50$\pm$0.11 | 84.83$\pm$0.23 | 85.90$\pm$0.20 | 86.29$\pm$0.23 | 81.26$\pm$0.18 | 86.70$\pm$0.24
BTCV, IRCAD, LITS | CHAOS | 87.22$\pm$0.14 | 90.42$\pm$0.08 | 91.08$\pm$0.19 | 90.93$\pm$0.16 | 91.18$\pm$0.16 | 84.38$\pm$0.15 | 91.93$\pm$0.19
BTCV, CHAOS, LITS | IRCAD | 89.17$\pm$0.12 | 89.26$\pm$0.13 | 89.17$\pm$0.09 | 91.44$\pm$0.24 | 90.89$\pm$0.10 | 90.10$\pm$0.13 | 92.14$\pm$0.09
BTCV, CHAOS, IRCAD | LITS | 88.67$\pm$0.08 | 88.48$\pm$0.21 | 89.41$\pm$0.15 | 89.13$\pm$0.21 | 90.11$\pm$0.18 | 88.19$\pm$0.24 | 89.73$\pm$0.26
Average | 86.37 | 87.67 | 88.63 | 89.16 | 89.54 | 85.98 | 90.20

Figure 8: Training Loss.

Figure 9: Sensitivity study on$\gamma_{1}$, $\gamma_{2}$.

Training stability and parameter sensitivity. To evaluate the training
stability of the proposed ETTA-SE framework, the curves of the loss terms
involved in the training procedure are all drawn in Figure 9. We observe that
they drop in a systematic way as the training goes, which shows the stability
of our system. We also see that $\mathcal{L}_{task}$ decreases in the early
training stage and then oscillates around low values. It demonstrates that the
meta-objective for deep embedded alignment, i.e., $\mathcal{L}_{fa}$ and
$\mathcal{L}_{ca}$, does perform the regularization effect on the model. That
is, the model is promoted to perform the appropriately well rather than
extremely well on the training set, which is more in line with the goal of
generalizing to unseen domains.

In addition, we evaluate the sensitivity of our ETTA-SE with respect to the
trade-off parameters of meta-objective, i.e., $\gamma_{1}$ and $\gamma_{2}$.
Figure 9 shows the performance of our proposed model when $\gamma_{1}$
($\gamma_{2}$) varies and the other parameter $\gamma_{2}$ ($\gamma_{1}$) is
fixed. It appears that the performance of our model is not sensitive to the
variations of $\gamma_{1}$ ($\gamma_{2}$) especially when $\gamma_{1}$
($\gamma_{2}$) varies from 0.25 (0.25) to 1.5 (0.75).

Figure 10: Visual results of DG on liver segmentation.

###  4.4 Experiments on the liver segmentation in CT images

####  4.4.1 Comparison with Baseline Model

We also evaluate our DG method on the liver segmentation in CT images. Table 4
presents the results of the DeepAll baseline and our proposed ETTA-SE method.
Although the domain shift is not easy to observe directly from the grayscale
CT scans, it makes the cross-domain performance of liver segmentation drops a
lot. In addition, the degree of performance degradation varies under different
unseen test domains, which shows the diverse difficulties of these DG
missions. Compared to the DeepAll baseline, our proposed ETTA-SE achieves the
improvement as 4.22% Dice score on average. As shown in Figure 10, the DeepAll
baseline fails to accurately segment the liver regions in some cases. In
comparison, the results yielded by our method are more refined, and exhibit
the greater visual similarity to the ground-truth.

####  4.4.2 Comparison with Other DG methods

We further report the quantitative results compared with the prior DG methods,
as shown in Table 2. We implement these approaches on the task of liver
segmentation based on the official public codes with some necessary
modifications. When the source domains are BTCV, CHAOS, IRCAD and the unseen
test domain is LITS, the proposed ETTA-SE method does not outperforms MASF
[16]. However, our method can still achieve the state-of-the-art performance
on all the remaining settings. The proposed method also surpasses other
advanced methods on the average Dice score. As shown in Figure 10, some over-
segmented and under-segmented results produced by prior methods can be refined
by our method.

##  5 Conclusion

In this paper, we propose the novel episodic training with task augmentation
method for domain generalization (DG) on medical imaging classification. The
episodic training paradigm transfers the knowledge of DG from meta-task
simulation to the real task of DG. Inspired by the limited number of available
source domains in clinical practice, we consider the special risk of
overfitting to training meta-tasks. To alleviate it, task augmentation is
developed to enhance the variety during training task generation, using mixed
task sampling (MTS). The meta-objective is integrated into the established
learning framework to regularize the deep semantic embedding to align across
training domains. The results of extensive experiments demonstrate our state-
of-the-art performance for DG.

## References

  * [1] C. Downey, S. Simpkins, J. White, D. Holliday, J. Jones, L. Jordan, J. Kulka, S. Pollock, S. Rajan, H. Thygesen, et al., “The prognostic significance of tumour–stroma ratio in oestrogen receptor-positive breast cancer,” British journal of cancer, vol. 110, no. 7, pp. 1744–1747, 2014\.
  * [2] T. Heimann, B. Van Ginneken, M. A. Styner, Y. Arzhaeva, V. Aurich, C. Bauer, A. Beck, C. Becker, R. Beichel, G. Bekes, et al., “Comparison and evaluation of methods for liver segmentation from ct datasets,” IEEE transactions on medical imaging, vol. 28, no. 8, pp. 1251–1265, 2009.
  * [3] B. J. Erickson, P. Korfiatis, Z. Akkus, and T. L. Kline, “Machine learning for medical imaging,” Radiographics, vol. 37, no. 2, pp. 505–515, 2017.
  * [4] D. Shen, G. Wu, and H.-I. Suk, “Deep learning in medical image analysis,” Annual review of biomedical engineering, vol. 19, pp. 221–248, 2017.
  * [5] L. G. Epstein and M. Schneider, “Iid: independently and indistinguishably distributed,” Journal of Economic Theory, vol. 113, no. 1, pp. 32–50, 2003\.
  * [6] M. Dundar, B. Krishnapuram, J. Bi, and R. B. Rao, “Learning classifiers when the training data is not iid.,” in IJCAI, vol. 2007, pp. 756–61, 2007\.
  * [7] M. Ghafoorian, A. Mehrtash, T. Kapur, N. Karssemeijer, E. Marchiori, M. Pesteie, C. R. Guttmann, F.-E. de Leeuw, C. M. Tempany, B. Van Ginneken, et al., “Transfer learning for domain adaptation in mri: Application in brain lesion segmentation,” in International conference on medical image computing and computer-assisted intervention, pp. 516–524, Springer, 2017\.
  * [8] Q. Dou, C. Ouyang, C. Chen, H. Chen, and P.-A. Heng, “Unsupervised cross-modality domain adaptation of convnets for biomedical image segmentations with adversarial loss,” arXiv preprint arXiv:1804.10916, 2018\.
  * [9] Y. Zhang, S. Miao, T. Mansi, and R. Liao, “Task driven generative modeling for unsupervised domain adaptation: Application to x-ray image segmentation,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 599–607, Springer, 2018.
  * [10] C. Li, Y. Zhang, Z. Liang, W. Ma, Y. Huang, and X. Ding, “Consistent posterior distributions under vessel-mixing: A regularization for cross-domain retinal artery/vein classification,” arXiv preprint arXiv:2103.09097, 2021.
  * [11] C. Chen, Q. Dou, H. Chen, J. Qin, and P.-A. Heng, “Synergistic image and feature adaptation: Towards cross-modality domain adaptation for medical image segmentation,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 33, pp. 865–872, 2019.
  * [12] H. Li, S. J. Pan, S. Wang, and A. C. Kot, “Domain generalization with adversarial feature learning,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5400–5409, 2018.
  * [13] F. M. Carlucci, A. D’Innocente, S. Bucci, B. Caputo, and T. Tommasi, “Domain generalization by solving jigsaw puzzles,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2229–2238, 2019.
  * [14] L. Zhang, X. Wang, D. Yang, T. Sanford, S. Harmon, B. Turkbey, B. J. Wood, H. Roth, A. Myronenko, D. Xu, et al., “Generalizing deep learning for medical image segmentation to unseen domains via deep stacked transformation,” IEEE transactions on medical imaging, vol. 39, no. 7, pp. 2531–2540, 2020.
  * [15] C. Yoon, G. Hamarneh, and R. Garbi, “Generalizable feature learning in the presence of data bias and domain class imbalance with application to skin lesion classification,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 365–373, Springer, 2019.
  * [16] Q. Dou, D. C. de Castro, K. Kamnitsas, and B. Glocker, “Domain generalization via model-agnostic learning of semantic features,” in Advances in Neural Information Processing Systems, pp. 6447–6458, 2019.
  * [17] Q. Liu, C. Chen, J. Qin, Q. Dou, and P.-A. Heng, “Feddg: Federated domain generalization on medical image segmentation via episodic learning in continuous frequency space,” arXiv preprint arXiv:2103.06030, 2021.
  * [18] Y. Li, X. Tian, M. Gong, Y. Liu, T. Liu, K. Zhang, and D. Tao, “Deep domain generalization via conditional invariant adversarial networks,” in Proceedings of the European Conference on Computer Vision (ECCV), pp. 624–639, 2018.
  * [19] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” in Proceedings of the 34th International Conference on Machine Learning-Volume 70, pp. 1126–1135, JMLR. org, 2017.
  * [20] F. Bianconi, A. Álvarez-Larrán, and A. Fernández, “Discrimination between tumour epithelium and stroma via perception-based features,” Neurocomputing, vol. 154, pp. 119–126, 2015.
  * [21] R. Nava, G. González, J. Kybic, and B. Escalante-Ramírez, “Classification of tumor epithelium and stroma in colorectal cancer based on discrete tchebichef moments,” in Workshop on Clinical Image-Based Procedures, pp. 79–87, Springer, 2015.
  * [22] Y. Huang, H. Zheng, C. Liu, G. Rohde, D. Zeng, J. Wang, and X. Ding, “Epithelium-stroma classification in histopathological images via convolutional neural networks and self-taught learning,” in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1073–1077, IEEE, 2017.
  * [23] Q. Qi, Y. Li, J. Wang, H. Zheng, Y. Huang, X. Ding, and G. K. Rohde, “Label-efficient breast cancer histopathological image classification,” IEEE journal of biomedical and health informatics, vol. 23, no. 5, pp. 2108–2116, 2018.
  * [24] Y. Huang, H. Zheng, C. Liu, X. Ding, and G. K. Rohde, “Epithelium-stroma classification via convolutional neural networks and unsupervised domain adaptation in histopathological images,” IEEE journal of biomedical and health informatics, vol. 21, no. 6, pp. 1625–1632, 2017.
  * [25] Q. Qi, X. Lin, C. Chen, W. Xie, Y. Huang, X. Ding, X. Liu, and Y. Yu, “Curriculum feature alignment domain adaptation for epithelium-stroma classification in histopathological images,” IEEE Journal of Biomedical and Health Informatics, 2020.
  * [26] D. Kainmüller, T. Lange, and H. Lamecker, “Shape constrained automatic segmentation of the liver based on a heuristic intensity model,” in Proc. MICCAI Workshop 3D Segmentation in the Clinic: A Grand Challenge, vol. 109, p. 116, 2007.
  * [27] A. Wimmer, G. Soza, and J. Hornegger, “A generic probabilistic active shape model for organ segmentation,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 26–33, Springer, 2009\.
  * [28] S. D. S. Al-Shaikhli, M. Y. Yang, and B. Rosenhahn, “Automatic 3d liver segmentation using sparse representation of global and local image information via level set formulation,” arXiv preprint arXiv:1508.01521, 2015.
  * [29] Q. Dou, H. Chen, Y. Jin, L. Yu, J. Qin, and P.-A. Heng, “3d deeply supervised network for automatic liver segmentation from ct volumes,” in International conference on medical image computing and computer-assisted intervention, pp. 149–157, Springer, 2016.
  * [30] D. Yang, D. Xu, S. K. Zhou, B. Georgescu, M. Chen, S. Grbic, D. Metaxas, and D. Comaniciu, “Automatic liver segmentation using an adversarial image-to-image network,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 507–515, Springer, 2017.
  * [31] M. Long, H. Zhu, J. Wang, and M. I. Jordan, “Unsupervised domain adaptation with residual transfer networks,” in Advances in neural information processing systems, pp. 136–144, 2016.
  * [32] K. Saito, K. Watanabe, Y. Ushiku, and T. Harada, “Maximum classifier discrepancy for unsupervised domain adaptation,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 3723–3732, 2018\.
  * [33] Y. Zhu, F. Zhuang, and D. Wang, “Aligning domain-specific distribution and classifier for cross-domain classification from multiple sources,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 33, pp. 5989–5996, 2019.
  * [34] K. Muandet, D. Balduzzi, and B. Schölkopf, “Domain generalization via invariant feature representation,” in International Conference on Machine Learning, pp. 10–18, PMLR, 2013.
  * [35] S. Motiian, M. Piccirilli, D. A. Adjeroh, and G. Doretto, “Unified deep supervised domain adaptation and generalization,” in 2017 IEEE International Conference on Computer Vision (ICCV), pp. 5716–5726, Oct 2017\.
  * [36] M. Ilse, J. M. Tomczak, C. Louizos, and M. Welling, “Diva: Domain invariant variational autoencoders,” arXiv preprint arXiv:1905.10427, 2019.
  * [37] D. Li, Y. Yang, Y.-Z. Song, and T. M. Hospedales, “Deeper, broader and artier domain generalization,” in Proceedings of the IEEE international conference on computer vision, pp. 5542–5550, 2017.
  * [38] D. Li, J. Zhang, Y. Yang, C. Liu, Y.-Z. Song, and T. M. Hospedales, “Episodic training for domain generalization,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 1446–1455, 2019.
  * [39] S. Wang, L. Yu, K. Li, X. Yang, C.-W. Fu, and P.-A. Heng, “Dofe: Domain-oriented feature embedding for generalizable fundus image segmentation on unseen datasets,” IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 4237–4248, 2020.
  * [40] R. Volpi, H. Namkoong, O. Sener, J. Duchi, V. Murino, and S. Savarese, “Generalizing to unseen domains via adversarial data augmentation,” arXiv preprint arXiv:1805.12018, 2018.
  * [41] D. Li, Y. Yang, Y.-Z. Song, and T. M. Hospedales, “Learning to generalize: Meta-learning for domain generalization,” in Thirty-Second AAAI Conference on Artificial Intelligence, 2018.
  * [42] Y. Balaji, S. Sankaranarayanan, and R. Chellappa, “Metareg: Towards domain generalization using meta-regularization,” in Advances in Neural Information Processing Systems, pp. 998–1008, 2018.
  * [43] Y. Li, Y. Yang, W. Zhou, and T. M. Hospedales, “Feature-critic networks for heterogeneous domain generalization,” arXiv preprint arXiv:1901.11448, 2019\.
  * [44] J. Snell, K. Swersky, and R. Zemel, “Prototypical networks for few-shot learning,” in Advances in neural information processing systems, pp. 4077–4087, 2017.
  * [45] W.-Y. Chen, Y.-C. Liu, Z. Kira, Y.-C. F. Wang, and J.-B. Huang, “A closer look at few-shot classification,” arXiv preprint arXiv:1904.04232, 2019.
  * [46] S. Gidaris and N. Komodakis, “Dynamic few-shot visual learning without forgetting,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 4367–4375, 2018.
  * [47] R. Hadsell, S. Chopra, and Y. LeCun, “Dimensionality reduction by learning an invariant mapping,” in 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), vol. 2, pp. 1735–1742, IEEE, 2006.
  * [48] A. H. Beck, A. R. Sangoi, S. Leung, R. J. Marinelli, T. O. Nielsen, M. J. Van De Vijver, R. B. West, M. Van De Rijn, and D. Koller, “Systematic analysis of breast cancer morphology uncovers stromal features associated with survival,” Science translational medicine, vol. 3, no. 108, pp. 108ra113–108ra113, 2011.
  * [49] N. Linder, J. Konsti, R. Turkki, E. Rahtu, M. Lundin, S. Nordling, C. Haglund, T. Ahonen, M. Pietikäinen, and J. Lundin, “Identification of tumor epithelium and stroma in tissue microarrays using texture analysis,” Diagnostic pathology, vol. 7, no. 1, pp. 1–11, 2012.
  * [50] J. N. Kather, J. Krisam, P. Charoentong, T. Luedde, E. Herpel, C.-A. Weis, T. Gaiser, A. Marx, N. A. Valous, D. Ferber, et al., “Predicting survival from colorectal cancer histology slides using deep learning: A retrospective multicenter study,” PLoS medicine, vol. 16, no. 1, p. e1002730, 2019.
  * [51] B. Landman, Z. Xu, J. Igelsias, M. Styner, T. Langerak, and A. Klein, “Miccai multi-atlas labeling beyond the cranial vault–workshop and challenge,” in Proc. MICCAI Multi-Atlas Labeling Beyond Cranial Vault—Workshop Challenge, 2015.
  * [52] A. E. Kavur, M. A. Selver, O. Dicle, M. Barış, and N. S. Gezer, “CHAOS - Combined (CT-MR) Healthy Abdominal Organ Segmentation Challenge Data,” Apr. 2019\.
  * [53] L. Soler, A. Hostettler, V. Agnus, A. Charnoz, J. Fasquel, J. Moreau, A. Osswald, M. Bouhadjar, and J. Marescaux, “3d image reconstruction for comparison of algorithm database: A patient specific anatomical and medical image database,” IRCAD, Strasbourg, France, Tech. Rep, 2010.
  * [54] P. Bilic, P. F. Christ, E. Vorontsov, G. Chlebus, H. Chen, Q. Dou, C.-W. Fu, X. Han, P.-A. Heng, J. Hesser, et al., “The liver tumor segmentation benchmark (lits),” arXiv preprint arXiv:1901.04056, 2019.
  * [55] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in Advances in neural information processing systems, pp. 1097–1105, 2012.
  * [56] L. v. d. Maaten and G. Hinton, “Visualizing data using t-sne,” Journal of machine learning research, vol. 9, no. Nov, pp. 2579–2605, 2008.