# CytOpT: Optimal Transport with Domain Adaptation for Interpreting Flow
Cytometry data

Paul Freulon 111correspondence<EMAIL_ADDRESS>Université de
Bordeaux, Bordeaux, 33000, France. Institut de Mathématiques de Bordeaux et
CNRS (UMR 5251), 33400 Talence, France.  Jérémie Bigot  Université de
Bordeaux, Bordeaux, 33000, France. Institut de Mathématiques de Bordeaux et
CNRS (UMR 5251), 33400 Talence, France.  Boris P. Hejblum  Université de
Bordeaux, Bordeaux, 33000, France. Bordeaux Population Health Research Center
Inserm U1219, Inria SISTM, 33000 Bordeaux, France.  Vaccine Research Institute
(VRI), 94010 Créteil, France.

###### Abstract

The automated analysis of flow cytometry measurements is an active research
field. We introduce a new algorithm, referred to as CytOpT, using regularized
optimal transport to directly estimate the different cell population
proportions from a biological sample characterized with flow cytometry
measurements. We rely on the regularized Wasserstein metric to compare
cytometry measurements from different samples, thus accounting for possible
mis-alignment of a given cell population across samples (due to technical
variability from the technology of measurements). In this work, we rely on a
supervised learning technique based on the Wasserstein metric that is used to
estimate an optimal re-weighting of class proportions in a mixture model from
a source distribution (with known segmentation into cell sub-populations) to
fit a target distribution with unknown segmentation. Due to the high-
dimensionality of flow cytometry data, we use stochastic algorithms to
approximate the regularized Wasserstein metric to solve the optimization
problem involved in the estimation of optimal weights representing the cell
population proportions in the target distribution. Several flow cytometry data
sets are used to illustrate the performances of CytOpT that are also compared
to those of existing algorithms for automatic gating based on supervised
learning.

Keywords: Automatic gating; flow cytometry; Optimal Transport; Stochastic
Optimization

##  1 Introduction

Flow cytometry is a high-throughput biotechnology used to characterize a large
amount of cells from a biological sample. Flow cytometry is paramount to many
biological and immunological research with applications, for instance, in the
monitoring of the immune system of HIV patients by counting the number of CD4
cells.

The first step to characterize cells from a biological sample with flow
cytometry is to stain those cells. Specifically, cells are stained with
multiple fluorescently-conjugated monoclonal antibodies directed to the
cellular markers of interest. Then, the cells flow one by one through the
cytometer laser beam. The scattered light is characteristic to the biological
markers of the cells N.Aghaeepour et al. (2013). Thus, from a biological
sample analyzed by a flow cytometer, we get a data set $X_{1},...,X_{I}$ where
each observation $X_{i}$ corresponds to a single cell crossing the laser beam.
For an observation $X_{i}\in\mathbb{R}^{d}$, the coordinate $X_{i}^{(m)}$
corresponds to the light intensity emitted by the fluorescent antibody
attached to the biological marker $m$. Interestingly, such a data set may also
be considered as a discrete probability distribution
$\frac{1}{I}\sum_{i=1}^{I}\delta_{X_{i}}$ with support in $\mathbb{R}^{d}$,
this is the point of view taken in this paper.

One can thus assess individual characteristics at a cellular level for samples
composed of hundreds of thousands of cells. The development of this technology
now leads to state-of-the-art cytometers that can measure up to 50 biological
markers at once on a single cell Y.Saeys et al. (2016). The constant increase
of the number of measured markers allows to identify more precisely the
different cell populations and subtypes present in the biological sample.

The analysis of cytometry data is generally done manually, by drawing
geometric shapes (referred to as “gates”) around populations of interest in a
sequence of two-dimensional data projections. Such manual analysis features
several drawbacks: i) it is extremely time-consuming ; ii) manual gating lacks
reproducibility across different operators N.Aghaeepour et al. (2013). To
overcome these shortcomings, several automated methods have been proposed
N.Aghaeepour et al. (2013). Those automated approaches aim at a clustering of
the flow cytometry data to derive the proportions of the cell populations that
are in the biological sample. Some methods follow an unsupervised approach.
For instance, FlowMeans Aghaeepour et al. (2011) is an automated method based
on the K-means designed for flow cytometry data. We also mention Cytometree
D.Commenges et al. (2018), an algorithm based on the construction of a binary
tree. Other unsupervised methods which perform model based clustering have
been proposed, see e.g. Y.Ge and S.C.Sealfon (2012); B.P.Hejblum et al. (2019)
To improve the accuracy of the classification, supervised machine learning
techniques have been applied to flow cytometry data analysis. Among those
techniques, one may cite DeepCyTOF introduced in H.Li et al. (2017) which is
based on deep-learning algorithms to gate cytometry data. In M.Lux et al.
(2018), the authors introduced flowlearn, a method that uses manually gated
samples to predict gates on other samples. We also want to mention a new
supervised approach named OptimalFlow developed by Barrio et al. (2019). This
method drew our attention as OptimalFlow introduces the Wasserstein distance
to quantify the discrepancy between cytometry data sets.

In spite of the numerous efforts to automate cytometry data analysis, manual
gating remains the gold-standard for benchmarking. The fact that automated
analysis has not outdone manual gating can be explained, at least in part, by
the significant variability of flow cytometry data. The variability of the
flow cytometry is twofold. First, variability is induced by biological
heterogeneity across the samples analyzed. For instance, we cannot expect the
same cell subtypes proportions within the biological sample of a healthy
patient and the biological sample of a sick patient. In addition to this
variability of the biological phenomena of interest, there is a technical
variability. This variability appears during the process of flow cytometry
analysis. For instance, differences in the staining procedure, in the data
acquisition settings or cytometers performances are very likely to happen and
lead to undesirable variability between flow cytometry data. This variability
of flow cytometry data makes the development of automated methods a
challenging task.

###  1.1 Illustrative data sets

As an illustrative example, we shall analyze in this paper the flow cytometry
data from the T-cell panel of the Human Immunology Project Consortium (HIPC) –
publicly available on ImmuneSpace R.Gottardo et al. (2014). Seven laboratories
stained three replicates (denoted A, B, and C) of three cryo-preserved
biological samples denoted patient 1, 2, and 3 (e.g. cytometry measurements
from the Stanford laboratory for replicate C from patient 1 will be denoted as
the “Stanford1C” data set). After performing cytometry measurements in each
center, the resulting FCS files were manually gated centrally for quantifying
10 cell populations: CD4 Effector (CD4 E), CD4 Naive (CD4 N), CD4 Central
memory (CD4 CM), CD4 Effector memory (CD4 EM), CD4 Activated (CD4 A), CD8
Effector (CD8 E), CD8 Naive (CD8 N), CD8 Central memory (CD8 CM), CD8 Effector
memory (CD8 EM) and CD8 Activated (CD8 A). Hence, for these data sets, a
manual clustering is at our disposal to evaluate the performances of automatic
gating methods. Those 10 sub populations have a size that range from 15 554
observations for the smallest data set to 112 318 observations for the largest
data set. For each cell we have the measurement of seven markers that thus
lead to cytometry observations $X_{i}$ that belong to $\mathbb{R}^{d}$ with
$d=7$. A 3D projection using three markers is displayed in Figure 1 for the
“Stanford1A” data set with the corresponding manual gating into 10 clusters.
Additionally, we also benchmarked our method on the flow cytometry data sets
used in Barrio et al. (2019).

Figure 1: Example flow cytometry data set of a biological sample from Stanford
patient 1 replicate A. Left: Manual clustering of the cytometry data. Existing
automated methods target a clustering in order to derive the class
proportions. Right: Class proportion derived from the manual gating. Our
method CytOpT aims at going straight to the estimation of class proportions
without clustering the cytometry data.

###  1.2 Main contributions

We propose a new supervised method to estimate the relative proportions of the
different cell sub-types in a biological sample analyzed by a cytometer. Our
approach aims at finding an optimal re-weighting of class proportions in a
mixture model between a source data set (with known segmentation into cell
sub-populations) to fit a target data set with unknown segmentation. This
estimation of the class proportions is done without any preliminary clustering
of the target cytometry data. To the best of our knowledge, all the previous
automated methods are based on the classification of the cells of the
biological sample to deduce a class proportions estimation. However, from a
clinical perspective, the relevant information is the class proportions, while
the clustering of individual cells is simply a mean to an end to get there.
Our method can nonetheless be extended in order to obtain a clustering of
cytometry data.

We believe that going straight to the estimation of class proportions is an
original approach to tackle the issue of flow cytometry data analysis. Hence,
the result of our algorithm is not a clustering of the cells analyzed by the
cytometers but a vector $h=(h_{1},...,h_{K})$, where each coefficient $h_{k}$
amounts to the percentage of the $k^{\text{th}}$ population cell among all the
cells analyzed in the target sample. To do so, we rely on tools derived from
optimal transport and regularized Wasserstein distance between probability
measures.

Optimal transport has recently gained interest in machine learning and
statistics. Indeed, the introduction of approximate solvers that can tackle
large dimension problems allowed to move beyond the high computational cost of
optimal transport. Thus, optimal transport has found various applications in
machine learning for regression H.Janati et al. (2018), classification
R.Flamary et al. (2018) and generative modeling M.Arjovsky et al. (2017).
Those computational progress in optimal transport have also allowed its use in
imaging sciences J.Solomon et al. (2015). Recent progress in applied optimal
transport has been fueled by the development of efficient, large-scale
optimization algorithms for problems in this domain. In particular, our method
rely on stochastic algorithms  A.Genevay et al. (2016); B.Bercu and J.Bigot
(2020) to handle the computational cost of optimal transport.

Then, to estimate the proportions $h=(h_{1},...,h_{K})$ of sub-populations
cell in a target data set (whose gating is unknown) by minimizing the
regularized Wasserstein distance between the target distribution and a re-
weighted vector of class proportions in the source distribution (from a source
data set whose gating into $K$ sub-populations cell of interest is known).
Since the Wasserstein distance is able to capture the underlying geometry of
the measurements, it has the ability to make meaningful comparisons between
distributions whose supports do not overlap and to quantify potential spatial
shifts. Our work demonstrates the benefits of using Wasserstein distance for
the analysis of the cytometry data as it allows to handle the technical
variability induced by different settings of measurements.

###  1.3 Organization of the paper

In Section 2 we present the Wasserstein distance and its regularized version.
We also introduce the stochastic algorithm that is used to compute the
regularized Wasserstein distance. Section 3 details our estimation method that
yields an estimate of the class proportions in an unsegmented data set. Then,
in Section 4, we test our approach on simulated data sets. Finally, in Section
5 we demonstrate the performance of CytOpT on the T-cell panel proposed by the
Human Immunology Project Consortium (HIPC) described previously. In this
Section 5, we also compare the performance of CytOpT, with the performances of
existing methods.

##  2 Mathematical framework

Optimal transport allows the definition of a metric between two probability
distributions $\alpha$ and $\beta$ supported on $\mathbb{R}^{d}$. This metric
is informally defined as the lowest cost to move the mass from one probability
measure, the source measure $\alpha$, onto the other, the target measure
$\beta$.

As optimal transport handles probability distributions, we must describe a
cytometry data set $X_{1},...,X_{I}$ where each $X_{i}$ belongs to
$\mathbb{R}^{d}$ as a probability distribution in $\mathbb{R}^{d}$. A natural
choice to move from observations to probability measure is to use the
empirical measure. Therefore, for a data set $X_{1},...,X_{I}$, we consider
its empirical measure defined as $\frac{1}{I}\sum_{i=1}^{I}\delta_{X_{i}}.$ In
this work, we shall consider optimal transport between discrete measures on
$\mathbb{R}^{d}$, and we denote by $\alpha=\sum_{i=1}^{I}a_{i}\delta_{X_{i}}$
and $\beta=\sum_{j=1}^{J}b_{j}\delta_{Y_{j}}$ two such measures. Note that
$a_{i}$ is not necessarily equal to $1/I$ as a re-weighting of the source
measure will be considered. Optimal transport then seeks a plan that moves the
probability measure $\alpha$ to the probability measure $\beta$ with a minimum
cost. The transportation cost is encoded by a function
$c:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$ where $c(x,y)$
represents the cost to move one unit of mass from $x$ to $y$. In the discrete
setting the transportation cost is encoded by a matrix
$C\in\mathbb{R}^{I\times J}$ where $C_{i,j}=c(X_{i},Y_{j})$.

By denoting $\Pi(\alpha,\beta)=\\{P\in\mathbb{R}^{I\times
J}:P1_{J}=a\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
P^{T}1_{I}=b\\}$ all the coupling matrices from $\alpha$ to $\beta$,
minimizing the transportation cost from $\alpha$ to $\beta$ boils down to the
following optimization problem:

$W_{c}(\alpha,\beta)=\inf_{\pi\in\Pi(\alpha,\beta)}\sum_{i=1}^{I}\sum_{j=1}^{J}C_{i,j}P_{i,j}=\inf_{\pi\in\Pi(\alpha,\beta)}\langle
C,P\rangle$ (2.1)

All along this work the coefficient $C_{i,j}$ of the cost matrix is defined as
the squared Euclidean distance between $x_{i}$ and $y_{j}$, i.e.
$C_{i,j}=||x_{i}-y_{j}||_{2}^{2}$, and we use the notation
$W(\alpha,\beta)=W_{2}^{2}(\alpha,\beta)$. However, the cost to evaluate the
Wasserstein distance between two discrete probability distributions with
support of equal size $N$ is generally of order $O(N^{3}\log(N))$. To allow
the evaluation of the Wasserstein distance between large data sets M.Cuturi
(2013) has proposed to add an entropic regularization term to the linear
optimization problem (2.1) in order to accelerate the computation of an
approximate solution. The regularized problem is

$W^{\varepsilon}(\alpha,\beta)=\inf_{\pi\in\Pi(\alpha,\beta)}\langle
C,P\rangle+\varepsilon H(P),$ (2.2)

where $H:\mathbb{R}_{+}^{I\times J}\rightarrow\mathbb{R}$ is defined for
$P\in\mathbb{R}_{+}^{I\times J}$ by
$H(P)=\sum_{i,j}\left(\log\left(P_{i,j}\right)\right.$
$\left.-1\right)P_{i,j}$. The entropic regularization of the Kantorovich
problem (2.2) leads to an approximation of the Wasserstein distance which can
be calculated in $O(N^{2}\log(N))$ operations if the two distributions have a
support of size $N$ G.Peyré et al. (2019).

###  2.1 Calculation of regularized Wasserstein distances

In this work, the regularized Wasserstein distance is calculated with a
statistical procedure based on the Robbins-Monro algorithm for stochastic
optimization. This way to calculate the Wasserstein distance is investigated
in A.Genevay et al. (2016) and B.Bercu and J.Bigot (2020). The keystone of
this approach is that the regularized Wasserstein problem can be written as
the following stochastic optimization problem:

$W^{\varepsilon}(\alpha,\beta)=\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u)],$
(2.3)

where $Y$ is a random variable with distribution $\beta$, and
$g_{\varepsilon}$ is defined as:

$g_{\varepsilon}(y,u)=\sum_{i=1}^{I}u_{i}a_{i}+u_{c,\varepsilon}(y)-\varepsilon,$
(2.4)

where
$u_{c,\varepsilon}(y_{j})=\varepsilon\left(\log(b_{j})-\log\left(\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)\right)\right).$
We point out that the solution $P_{\varepsilon}$ of the primal problem (2.2)
is recovered from any $u^{*}\in\mathbb{R}^{I}$ solution of the semi-dual
problem (2.3) as:

$(P_{\varepsilon})_{i,j}=\exp\left(\frac{u_{i}^{*}+u_{c,\varepsilon}^{*}(y_{j})-c(x_{i},y_{j})}{\varepsilon}\right)$
(2.5)

Formulation (2.3) and the fact that for all $y\in\mathbb{R}^{d}$, the function
$g(y,.)$ is concave lead us to estimate the vector $u^{*}$ by the Robbins-
Monro algorithm H.Robbins and S.Monro (1951) given, for all $n\geq 0$, by:

$\widehat{U}_{n+1}=\widehat{U}_{n}+\gamma_{n+1}\nabla_{u}g_{\varepsilon}(Y_{n+1},\widehat{U}_{n}),$
(2.6)

where the initial value $\hat{U}_{0}$ is a random vector which can be
arbitrarily chosen, $Y_{1},...,Y_{n+1}$ is an i.i.d. sequence of random
variables sampled from the distribution $\beta$, and $(\gamma_{n})_{n\geq 0}$
is a positive sequence of real numbers decreasing toward zero satisfying

$\sum_{n=1}^{\infty}\gamma_{n}=+\infty\qquad\text{and}\qquad\sum_{n=1}^{\infty}\gamma_{n}^{2}<+\infty.$
(2.7)

It follows from M.Cuturi (2013) that the maximizer $u^{*}$ of (2.3) is unique
up to a scalar translation of the form $u^{*}-t1_{I}$ for any
$t\in\mathbb{R}$. Throughout this paper, we shall denote by $u^{*}$ the
maximizer of (2.3) satisfying $\langle u^{*},1_{I}\rangle=0$ which means that
$u^{*}$ belongs to $\langle 1_{J}\rangle^{\perp}$ where $\langle 1_{J}\rangle$
is the one-dimensional subspace of $\mathbb{R}^{I}$ spanned by $1_{J}$. As
shown in B.Bercu and J.Bigot (2020), to obtain a consistent estimator of
$u^{*}$, one must slightly modify the Robbins-Monro algorithm by requiring
that $\widehat{U}_{0}$ belongs to $\langle 1_{I}\rangle^{\perp}$. To
approximate the regularized Wasserstein distance
$W^{\varepsilon}(\alpha,\beta)$ B.Bercu and J.Bigot (2020) proposed the
recursive estimator:

$\widehat{W}_{n}=\frac{1}{n}\sum_{k=1}^{n}g_{\varepsilon}(Y_{k},\widehat{U}_{k-1}).$
(2.8)

Throughout this work, we will use this estimator.

###  2.2 Statistical model

Let us consider $X_{1}^{s},...,X_{I}^{s}$ the cytometry measurements from a
first biological sample that is referred to as the source sample or source
observations. The distribution of the source sample is modeled by a mixture of
$K$ distributions: $\alpha=\sum_{k=1}^{K}\rho_{k}\alpha_{k}$, where each term
of the mixture corresponds to the cytometry measurements of one type of cell.
For instance, $\alpha_{k}$ represents the underlying distribution behind the
cytometry measurements of the CD4 effector T cells that are in the biological
sample. The weights $\rho\in\Sigma_{K}$ are coefficients lying in the
probability simplex $\Sigma_{K}=\\{h\in\mathbb{R}_{+}^{K}\leavevmode\nobreak\
:\leavevmode\nobreak\ \sum_{k=1}^{K}h_{k}=1\\}$, and $\rho_{k}$ represents the
proportion of one cell sub-type among all the cells found in the sample.

We also consider a second set of cytometry measurements
$X_{1}^{t},...,X_{J}^{t}$ that is referred to as the target sample or source
observations, and that corresponds to a second biological sample that may come
from another patient or a cytometry analysis performed in a different
laboratory. We assume that the underlying distribution of the the target
observations is another mixture of $K$ distributions:
$\beta=\sum_{k=1}^{K}\pi_{k}\beta_{k}$, where $\pi\in\Sigma_{K}$ represents
the unknown class proportions and $\beta_{k}$ corresponds to the distribution
of the $k^{\text{th}}$ cell population. Note that, in the framework of flow
cytometry, we cannot make the assumption that $\alpha_{k}=\beta_{k}$, neither
that $\rho_{k}=\pi_{k}$. Indeed, $\rho$ could differ from $\pi$ due to
biological differences. For instance if the source sample comes from a healthy
patient and the target sample comes from a sick patient, it is likely that
$\pi$ will feature significant difference from $\rho$. Moreover, between each
component $\beta_{k}$ and $\alpha_{k}$ there could be a difference (e.g.
location shift) due to the technical variability of the cytometry
measurements.

These potential differences between the distributions of two flow cytometry
measurements make the development of supervised methods to infer the unknown
proportions $\pi$ in the target data from those of the source data a difficult
task. In this work, we rely on the geometric properties of the Wasserstein
distance to handle the differences between samples that are due to technical
reasons or inter-variability between healthy and sick patients.

##  3 Class proportion estimation

We now detail a supervised algorithm that estimates the class proportions
$\pi$ in the target distribution. While state-of-the-art automated methods in
cytometry data analysis attempt to classify the observations
$X_{1}^{t},...,X_{J}^{t}$, the method that we propose will go straight to the
estimation of the class proportions in the unsegmented target data set. To do
so, we borrow ideas from the domain adaptation technique proposed in I.Redko
et al. (2018) where it is proposed to re-weight the source observations by
searching for weights minimizing the Wasserstein distance between the re-
weighted source measure and the target measure. Contrary to I.Redko et al.
(2018), we handle the Wasserstein distance and the minimization problem by a
stochastic approach. Indeed the large number of observations produced by flow
cytometry makes stochastic techniques competitive to handle the high-
dimensionality of such data sets.

###  3.1 A new estimator of the class proportions

We now specify the definition of the estimator of the class proportions in the
target data set. For the target sample $X_{1}^{t},...,X_{J}^{t}$, the
segmentation into various cell sub-types is not available, hence we define the
empirical target measure as
$\hat{\beta}=\frac{1}{J}\sum_{j=1}^{J}\delta_{X_{j}^{t}}$. From the source
observations $X_{1}^{s},...,X_{I}^{s}$, we define the empirical source measure
as $\hat{\alpha}=\frac{1}{I}\sum_{i=1}^{I}\delta_{X_{i}^{s}}$

Then, the knowledge of the gating of the source data allows us to re-write the
measure $\hat{\alpha}$ as a mixture of probability measures where each
component corresponds to a known sub-population of cells in the data

$\hat{\alpha}=\sum_{k=1}^{K}\frac{n_{k}}{I}\left(\sum_{i:X_{i}^{s}\in
C_{k}}\frac{1}{n_{k}}\delta_{X_{i}^{s}}\right)=\sum_{k=1}^{K}\frac{n_{k}}{I}\hat{\alpha}_{k},$
(3.1)

where $n_{k}=\\#C_{k}$ and $\hat{\alpha}_{k}=\sum_{i:X_{i}^{s}\in
C_{k}}\frac{1}{n_{k}}\delta_{X_{i}^{s}}$. Namely, the component
$\hat{\alpha}_{k}$ is the empirical measure of the observations that belong to
the sub-population $C_{k}$ (known class).

Then, instead of only considering the true class proportions
$(n_{1}/I,...,n_{K}/I)$ in the source data set, we can re-weight the clusters
$C_{k}$ in the empirical distribution as desired. Indeed, for a probability
vector $h=(h_{1},...,h_{K})\in\Sigma_{K}$ we can define the measure
$\hat{\alpha(h)}$ that corresponds to the re-weighted measure $\hat{\alpha}$
such that for all $k\in\\{1,...,K\\}$ the component $\hat{\alpha}_{k}$ amounts
for $h_{k}$ in the measure $\hat{\alpha}(h)$. In mathematical terms, the
measure $\hat{\alpha}(h)$ is thus defined by:

$\hat{\alpha}(h)=\sum_{k=1}^{K}h_{k}\hat{\alpha}_{k}.$ (3.2)

Then, to derive the class proportions in the target data, we minimize the
regularized Wasserstein distance (2.2) between the re-weighted source
empirical distribution (3.2) and the target empirical distribution. The main
idea is that the source distribution will get closer to the target
distribution as the class proportions in its re-weighted version get closer to
the class proportions of the target distribution. Thus, we propose to estimate
the weights $\pi=(\pi_{1},...,\pi_{K})\in\Sigma_{K}$ of the underlying
distribution $\beta$ behind the observations $X_{1}^{t},...,X_{J}^{t}$ of the
unlabelled target data set by

$\hat{\pi}\in\operatorname*{arg\,min}_{h\in\Sigma_{K}}W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})$
(3.3)

The combination of the estimator $\hat{\pi}$ and the algorithms described in
Subsection 3.4 to solve the associated minimization problem (3.3) will be
referred to as CytOpT.

###  3.2 Extension of CytOpT to a soft assignment method

While our method aims at estimating the class proportions and not to classify
the target observations, regularized optimal transport offers a natural soft
assignment method which can be used to derive a soft classification of the
target data set, as illustrated in Figure 7. Assuming that we have access to
the optimal transport plan $P_{\varepsilon}$ with respect to the regularized
problem (2.2), the coefficient $(P_{\varepsilon})_{i,j}/b_{j}$ can be
interpreted as the probability that $X_{i}^{s}$ is assigned to $X_{j}^{t}$.
Here, $b_{j}$ is the weight associated to the observation $X_{j}^{t}$. Thus,
the probability $\gamma_{j}^{(k)}$ that $X_{j}^{t}$ belongs to the class
$C_{k}$ is $\gamma_{j}^{(k)}=\frac{1}{b_{j}}\sum_{i=1}^{I}1_{X_{i}^{s}\in
C_{k}}(P_{\varepsilon})_{i,j}$ By choosing the class with highest probability,
we can derive a classification for the observation $X_{j}^{t}$. One can
compute an approximate $\hat{P}$ of $P_{\varepsilon}$ by plugging
$\widehat{U}$ the Robbins-Monro approximation of the optimal dual vector
$u^{*}$ in formula $\eqref{eq:dual_to_primal}$.

Thus, we can derive an automatic gating of the target data thanks to optimal
transport. However, the transfer of the classification from the source data
set toward the target data set requires to compute all the columns of the
optimal transport plan. Therefore, obtaining a classification in the target
data calls for additional calculations. Nevertheless, in flow cytometry, the
information that matters from a clinical point of view is the relative
proportions of the different cell sub-types. Hence, we focus on the estimation
of class proportions in this work.

###  3.3 Measure of performance

To evaluate our approach, the benchmark class proportions
$\pi=(\pi_{1},...,\pi_{K})\in\Sigma_{K}$, will be defined thanks to the manual
segmentation of the target observations. We recall that our algorithm does not
make use of the segmentation of the target data set, and that it is only used
to evaluate our method.

In flow cytometry data analysis the F-measure is a popular tool to asses the
performances of the clustering methods. As our method does not yield a
clustering but an estimate $\hat{\pi}$ of the class proportions $\pi$, we
cannot rely on this measure. Indeed, CytOpT yields a probability vector
$\hat{\pi}\in\Sigma_{K}$ that we wish to be the closest to the class
proportions $\pi\in\Sigma_{K}$ in the target data set. Hence, a natural way to
measure the discrepancy between $\hat{\pi}$ and $\pi$, is the Kullback-
Leilbler divergence. Thus, with an estimation
$\hat{\pi}=(\hat{\pi}_{1},...,\hat{\pi}_{K})$ and a benchmark
$\pi=(\pi_{1},...,\pi_{K})$, to assess the quality of the estimator
$\hat{\pi}$, we compute the Kullback-Leibler divergence $\operatorname*{KL}$
defined as
$\operatorname*{KL}(\hat{\pi}|\pi)=\sum_{k=1}^{K}\hat{\pi}_{k}\log\left(\frac{\hat{\pi}_{k}}{\pi_{k}}\right)$.

We also relied on graphical diagnoses from Bland-Altman plots (JM.Bland and
DG.Altman (1986)) to visually assess and compare the performance of CytOpT for
class proportions estimation. Those provide an overview of CytOpT behavior and
performance on a collection of several cytometry data sets at once. For one
target data set, a Bland-Altman plot compares the estimation
$\hat{\pi}=(\hat{\pi}_{1},...,\hat{\pi}_{K})$ with the benchmark
$\pi=(\pi_{1},...,\pi_{K})$ by plotting the difference $\hat{\pi}-\pi$ against
the mean $(\hat{\pi}+\pi)/2$. To simultaneously visualize the results from the
analysis of two data sets, one just adds on to the graph the points defined by
$\hat{\pi}^{(2)}-\pi^{(2)}$ for the y-axis and $(\hat{\pi}^{(2)}+\pi^{(2)})/2$
for the x-axis (with $\hat{\pi}^{(2)}$ denoting the estimation on the second
data set and $\pi^{(2)}$ the benchmark on that data set). This way, one can
actually represent the results from the analysis of several data sets in a
single plot. As we aim for $\hat{\pi}$ to be as close as possible to $\pi$,
the closer to the $x$-axis the points are, the better the performance.

###  3.4 A brief overview of the minimization procedures

The optimization problem (3.3) leading to the estimator $\hat{\pi}$ of the
class proportions does not have a solution in a closed form expression, and a
numerical approximation is needed. We propose two methods to solve this
optimization problem. In this section, we only give the mains ideas of these
two methods without focusing on the technical details. For a more thorough
presentation, we refer the reader to Appendix A.

####  3.4.1 Descent-Ascent procedure

First, we re-write problem (3.3) as

$\min_{h\in\Sigma_{K}}W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})=\min_{h\in\Sigma_{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,h)],$
(3.4)

where $Y$ is a random variable with distribution $\hat{\beta}$, and
$g_{\varepsilon}$ is the function defined by (2.4) that can be easily
computed. Second, to avoid projecting $h$ on the simplex $\Sigma_{K}$ at each
iteration, we re-parameterize problem $\eqref{eq:min_max_expectation_one}$
with a soft max function $\sigma$ whose complete expression is given in (A.5).
Hence, our new minimization problem is

$\min_{z\in\mathbb{R}^{K}}W^{\varepsilon}(\hat{\alpha}(\sigma(z)),\hat{\beta})=\min_{z\in\mathbb{R}^{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,\sigma(z))]$
(3.5)

The gradient of the function $F:z\mapsto
W^{\varepsilon}(\hat{\alpha}(\sigma(z)),\hat{\beta})$ can be expressed for
$z\in\mathbb{R}^{K}$ as $\nabla_{z}F(z)=(\Gamma J_{\sigma}(z))^{T}u_{z}^{*}$
where $u_{z}^{*}$ is a maximizer of the function
$u\mapsto\mathbb{E}[g_{\varepsilon}(Y,u,\sigma(z))]$. The expressions of
$\Gamma$ and $J_{\sigma}$ can be found in (A.1) and (A.8) respectively. Hence,
a natural way to solve (3.5) is to use a gradient descent ascent algorithm
where we alternate between a gradient descent on the variable $z$ and a
gradient ascent on the variable $u$. The approximate of the variable
$u_{z}^{*}$ allows us to compute the gradient in the outer loop. An
interesting feature of our algorithm is that to compute an approximate $U_{z}$
of the vector $u_{z}^{*}$ we rely on the Robbins-Monro algorithm (2.6). Then,
by computing $\hat{\omega}(z)=(\Gamma J_{\sigma}(z))^{T}U_{z}$ we get a
stochastic approximation of $\nabla_{z}F(z)$ to proceed with the gradient
ascent.

$z\leftarrow 1_{K}$

for  _$l\leftarrow 1$ to $n_{out}$_ do

$U\leftarrow$ an arbitrary vector in $\mathbb{R}^{I}$

for  _$k\leftarrow 1$ to $n_{in}$_ do

$Y\sim\hat{\beta}$

$U\leftarrow U+\gamma_{k}\nabla_{u}g_{\varepsilon}(Y,U,\sigma(z))$

end for

/* Approximation of the gradient of $z\mapsto
W^{\varepsilon}(\hat{\alpha}(\sigma(z)),\hat{\beta})$ */

$\hat{\omega}(z)\leftarrow(\Gamma J_{\sigma}(z))^{T}U$

$z\leftarrow z-\eta\hat{\omega}(z)$

end for

/* Computation of the estimator of the class proportion $\hat{\pi}$ */

return _$\hat{\pi}=\sigma(z)$_

Algorithm 1 Solving
$\min_{z\in\mathbb{R}^{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}_{Y\sim\hat{\beta}}[g_{\varepsilon}(Y,u,\sigma(z))]$

####  3.4.2 Minmax swapping procedure

To avoid the use of an inner-loop algorithm, we borrow ideas developed in
M.Ballu et al. (2020). By adding an entropic penalization $\varphi$ defined in
(A.13) to problem $\eqref{eq:estimator_proportion}$, we obtain the following
optimization problem:

$\min_{h\in\Sigma_{K}}W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})+\lambda\varphi(h)=\min_{h\in\Sigma_{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,h)]+\lambda\varphi(h).$
(3.6)

The entropic penalization allows to swap the min and the max. Then, for a
given $u\in\mathbb{R}^{I}$, we can compute an explicit solution
$h(u)\in\Sigma_{K}$ of the problem
$\min_{h\in\Sigma_{K}}\mathbb{E}[g_{\varepsilon}(Y,u,h)]+\lambda\varphi(h)$.
This solution $h(u)$ is defined by, $\forall
k\in\\{1,...,K\\},\leavevmode\nobreak\
(h(u))_{k}=\frac{\exp\left(-\frac{(\Gamma^{T}u)_{k}}{\lambda}\right)}{\sum_{l=1}^{K}\exp\left(-\frac{(\Gamma^{T}u)_{l}}{\lambda}\right)}$,
where $\Gamma$ the linear operator defined in (A.1). Thus, adding an entropy
penalization, induces a natural parameterization of the simplex by the soft-
max function. It follows from this result that problem (3.6) boils down to
$\max_{u\in\mathbb{R}^{I}}\mathbb{E}_{Y\sim\hat{\beta}}[f_{\varepsilon,\lambda}(Y,u)]$,
where the complete expression of $f_{\varepsilon,\lambda}$ can be found in
(A.20). Once again, we solve this problem using a Robbins-Monro algorithm,
similarly with (2.6).

$U\leftarrow 0_{I}$

for  _$l\leftarrow 1$ to $n$_ do

$Y\sim\hat{\beta}$

$U\leftarrow U+\gamma_{l}\nabla_{u}f_{\varepsilon,\lambda}(Y,U)$

end for

/* Computation of the estimator of the class proportion $\hat{\pi}$ */

$\hat{\pi}\leftarrow h(U)$

return _$\hat{\pi}$_

Algorithm 2 Solving
$\max_{u\in\mathbb{R}^{I}}\mathbb{E}_{Y\sim\hat{\beta}}[f_{\varepsilon,\lambda}(Y,u)]$

####  3.4.3 Comparison of the minimization procedures

As said in Subsection 3.3, we assess the performance of our estimator
$\hat{\pi}$ by computing the Kullback-Leibler divergence between $\hat{\pi}$
and $\pi$ the benchmark class proportions. Figure 2 displays the evolution of
the Kullback-Leibler divergence along the iterations of the two minimization
procedures. For the descent-ascent procedure, the iterations displayed
correspond to the iterations of the outer loop. This comparison has been
realized when the source data set is Stanford1A segmented into 10 classes. The
target data set is Stanford3A where we try to estimate the proportions of the
10 cell sub-populations.

Both algorithms have been implemented in Python. From our numerical
experiments, we have found that the sequence of estimates $\hat{\pi}$ produced
by the second algorithm described in 3.4.2 levels off approximately five times
faster than the sequence of estimates produced with the descent-ascent
procedure 3.4.1. This computational time gain can be accounted by the simple
loop complexity of this second algorithm. However, in practice, the descent-
ascent procedure seems to yield an estimate that is slightly more accurate
than the estimate produced by the second procedure. Indeed, after enough
iterations, the estimates of the descent-ascent procedure seem to get closer
(according to the Kullback-Leibler divergence) to the benchmark $\pi$ than the
estimates of the min-max swapping procedure. Hence, for all the experiments on
flow cytometry data presented below, we report results produced by the
descent-ascent procedure.

Figure 2: Evolution of the Kullback-Leibler divergence between the estimated
proportions $\hat{\pi}$ and the manual gating benchmark $\pi$.

##  4 Simulation Study

To assess the performance of our method we generate two data sets from two
different mixtures of K = 10 components each. The first data set is drawn from
a Gaussian mixture distribution in dimension $d=10$ with a vector of
proportions $\rho\in\Sigma_{K}$. Thus,
$\alpha=\sum_{k=1}^{10}\rho_{k}\alpha_{k}$ where the $\alpha_{k}$ are Gaussian
distributions. The second data set is drawn from the same Gaussian components
$\alpha_{k}$ but with different mixture proportions $\pi\neq\rho$. Then, we
add a non-linear mapping $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ mimicking domain-
shift. Thus, $\beta=\sum_{k=1}^{10}\pi_{k}\beta_{k}$ where
$\beta_{k}=T\\#\alpha_{k}$ (that is the push-forward of the measure
$\alpha_{k}$ by the mapping $T$). In this simulation study, we apply the same
pre-processing of the data as the one applied to the real flow cytometry data.
That is, we threshold the negative values at zero and re-scale the
observations such that every observation $X$ belongs to $[0,1]^{10}$. See
Figure 3 for an illustration of a generated data-set.

We compare our method with three state-of-the-art classification methods: one
unsupervised method – the $k$-means ; and two supervised methods – quadratic
discriminant analysis (QDA), and random forest (RF). After a training step on
the source data set, each supervised method yields a classification of the
target data set. Hence, after this classification step by either QDA or RF, we
can derive an estimation of the class proportions in the target data set. We
indicate, that in this framework (namely with a single learning data set), the
QDA algorithm corresponds to the OptimalFlow method Barrio et al. (2019) that
is a state-of-the-art method designed to cluster flow cytometry data. As the
$k$-means is an unsupervised method, we apply it straight to the target data
set. Since the $k$-means does not yield a labelled classification, it is
necessary to add an annotation step for each cluster $C_{k}$ returned by the
$k$-means algorithm. To annotate a cluster $C_{k}$, we retrieve the majority
population according to the true labels and then we label all the observations
that belong to $C_{k}$ with this label. Therefore, we provide a significant
upper hand to this method with this additional information.

To study the performance and the stability of these various methods, we
sampled 100 source data sets with respect to the source model (with a fixed
number $I=115,783$ of observations), and sampled 100 target data sets with
respect to the target model (with a fixed number $J=68,981$ of observations).
Then, for each pair of source and target data sets we apply the supervised
methods, that is CytOpT, QDA and RF, and on each target data set, we apply the
$k$-means algorithm followed by the annotation step described above. Hence,
for each method we obtain 100 estimates of the class proportions. Figure 3
displays an example of a pair of source and target data sets with their
corresponding class proportions.

Figure 4 displays the median class proportions retrieved with theses 100
estimates and the variability of these estimates. In this scenario, CytOpT is
the method which offers the best trade-off between accuracy and stability.

Figure 3: Simulated source data set and Simulated target data set. Using the
segmentation of the source data set (Left) CytOpt estimates the class
proportions in a target data set (Right) without making use of the
classification of the target data set. Figure 4: Performances of the automated
methods tested to estimate the class proportions. For each method, one boxplot
displays the median and the variation of the estimation of the weight of one
component.

##  5 Application to real flow cytometry data analysis

###  5.1 Illustration of the methodology with a two classes example

In this section, we test and illustrate our proposed method in the setting
where the cytometry data from HIPC described in Section 1.1 are divided only
among two broad classes: the CD4 cells and the CD8 cells. For ease of
visualization, we use only two markers: the CD4 marker and the CD8 marker –
that is $d=2$. This basic case, where two-dimensional data are divided into
two classes, is a first illustration of our method in a favourable situation.
We aim to demonstrate that our approach may reach a good approximation of the
CD4 proportion and the CD8 proportion in a target cytometry data set with
unknown gating.

We consider the two data sets from the HIPC T-cell panel that are displayed in
Figure 5. The first data set is a series of cytometry measurements done in a
Stanford laboratory on a biological sample that comes from a patient
identified as patient 1. The Stanford data set is chosen as the segmented
source data. The second data set is a series of cytometry measurements done in
the same laboratory on a biological sample that comes from an other patient
identified as patient 3. This second data set will be the target data set. The
manual gating classification, and thus the class proportions are available for
both samples. Nevertheless, we only use the classification available for the
source (i.e. Stanford1A) data set in order to estimate the class proportions
in the target (i.e. Stanford3A) data set. The manual gating classification for
the target data set is only used as a benchmark to evaluate our method.

Figure 5: Illustration of the CytOpT framework. CytOpt estimate the class
proportions in an unclassified data set (Stanford3A, right) from one
classified data set (Stanford1A, left) without clustering the observations of
Stanford3A.

First, to assess the relevance of our method, we present in Figure 6 the
evolution of the regularized Wasserstein distance as a function of the weights
associated to each class in the source distribution. To this end, we evaluate
the function $F:h_{1}\mapsto W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})$,
where $h=(h_{1},1-h_{1})$, on a finite grid
$\mathcal{H}=\\{h^{(1)},...,h^{(m)}\\}$. For $h_{1}\in\mathcal{H}$ we
approximate $W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})$ by the estimator
$\widehat{W}_{n}$ defined in equation $\eqref{eq:estimate_wasserstein}$. It
can be observed that the regularized Wasserstein distance decreases as the
class proportions of the source data set get closer to the class proportions
of the target data set.

Figure 6: Approximation of the function $h_{1}\mapsto
W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})$, for $h=(h_{1},1-h_{1})$,
$h_{1}$ represents the weight associated to the CD8 cells. The approximation
of $W^{\varepsilon}(\hat{\alpha}(h),\hat{\beta})$ is produced using the
estimator (2.8).

Using the segmentation of a cytometry data set, CytOpT adequately retrieves
the true class proportions of an unlabelled cytometry data set. Even if the
two classes situation is somewhat an easy scenario, one can observe the
significant gap in the class proportions between the source data set and the
target data set. In the Stanford1A data set the CD4 cells constitute 45.1 % of
the cells and the CD8 cells 54.9 % of the cells, whereas in the Stanford3A
data set the CD4 cells constitute 73.9 % (estimated at 73.3 % by CytOpT) of
the cells and the CD8 cells 26.1 % (estimated at 26.7 % by CytOpT).

Figure 7: Optimal transport plan $P_{\varepsilon}$ between the source and
target distribution. A green line between $X_{i}^{s}$ and $X_{j}^{t}$
indicates that the optimal transport plan $P_{\varepsilon}$ moves some mass
from $X_{i}^{s}$ to $X_{j}^{t}$. To facilitate readability, the target data
set has been shifted, and only 500 coefficients of $P_{\varepsilon}$ have been
represented. Without re-weighting, mass from the CD8 class in the source data
set is sent toward the CD4 class in the target data set. With re-weighting,
the mass from one class in the source data set is sent to the corresponding
class in the target data set.

The results presented in Figure 7 and Figure 8 represent the use of the soft-
clustering and classification methods described in Section 3.2, and they
outline the importance of re-weighting the source data. Without re-weighting,
the classification obtained by optimal transport is really unlike the manual
gating classification. However, one can re-weight the observations with the
transformation $a(\pi)=\Gamma\pi$ in order to match the class proportions
$\pi$ in the target data set. Here, $\Gamma$ is the linear operator defined in
(A.1). Thus, when we re-weight our observations with the estimated class
proportions $\hat{\pi}$ we get a much more interesting classification. Indeed,
once the source data are re-weighted, the classification by regularized
transport is very similar to the manual gating classification.

Figure 8: Results of the soft assignment method and comparison with the manual
clustering on the data set Stanford3A. From left to right: Label transfer
without reweighting - Label transfer with reweighting - Manual gating
benchmark

###  5.2 Estimating 10 cell sub-type proportions from 7 cellular markers

We apply our method to the full T-cell panel of the HIPC flow cytometry data
described in Section 1.1. Contrary to the illustrative results presented in
Section 5.1, we now use all of the d=7 markers at hand to estimate the
proportions of K = 10 cell populations. This represents a more involved
classification task than the one considered in Section 5.1 with $K=2$. We
arbitrarily chose the data set ”Stanford1A” as the source measure to derive
the class proportions in all the other data sets as targets. A comprehensive
evaluation of CytOpT performance is provided in Figure 9, which features a
Bland-Altman plot JM.Bland and DG.Altman (1986) displaying all the cell sub-
types proportions estimations by CytOpT when targeting all available data sets
across all the seven different centers available. In this example, we solely
used the information from one reference classification from manual gating of
one data set (namely Stanford1A) to estimate the class proportions in each of
the 61 unsegmented data sets targeted. In more than $90\%$ of the cases the
absolute difference between the estimated proportion and the manual gating
gold-standard proportion is below $5\%$. And in more than $99\%$ of the cases
the absolute error is no more than $10\%$. Due to the stochastic nature of our
algorithm, a new call to CytOpT would lead to a slightly different estimate
$\hat{\pi}$ compare to a former estimation. However, the results displayed in
Figure 9 and Figure 12 are representative of the general quality of the
estimation produced by CytOpT in such settings.

Finally, note that when applying our method to estimate class proportions on
real flow cytometry data, a simple pre-processing is required. First, the
signal processing of the cytometer can induce some contrived negative values
of light intensity Tung et al. (2007). To undo this effect, we merely
threshold those few negative values at zero. Second, to set the parameters of
our algorithms, in particular $\varepsilon$, we need to bound the displacement
cost. To do so, we scale the data such that: $\forall
i\in\\{1,...,I\\},\leavevmode\nobreak\ X_{i}^{s}\in[0,1]^{d}$ and $\forall
j\in\\{1,...,J\\},\leavevmode\nobreak\ X_{j}^{t}\in[0,1]^{d}$.

####  5.2.1 Comparison with other methods

On the HIPC data, we compare CytOpT with 3 state-of-the-art automated-gating
approaches specifically designed to analyze flow cytometry data. First,
flowMeans Aghaeepour et al. (2011) which is a variation of the $k$-Means
algorithm designed to cluster flow cytometry data. As flowMeans is an
unsupervised method, we will label every observation of a cluster retrieved by
flowMeans with the major population according to the reference manual
clustering (Note that this in theory should provide a substantial advantage
over the other methods). We also apply Cytometree to cluster the HIPC data.
This method offers an automated annotation of its estimated clustering.
Therefore, with Cytometree, we can directly retrieve an estimation of the
class proportions. Finally, we compare the performances of our method with
OptimalFlow, recently developed by Barrio et al. (2019) and which performs
supervised automatic gating. It is the first method that make use of the
Wasserstein distance in the context of automatic-gating of flow cytometry data
analysis. We stress that, even though both our method and OptimalFlow are
supervised approaches relying on the Wasserstein distance, they greatly
differ. First, they have different objectives: OptimalFlow aims at classifying
the cells, while CytOpT directly aims at estimating the cellular sub-types
proportions. Second, OptimalFlow uses the Wasserstein distance to cluster a
compendium of pre-gated cytometry data sets and uses the notion of Wasserstein
barycenter to produce a prototype data set for each cluster. Then, using the
Wasserstein distance, OptimalFlow browses among the prototypes to extract the
most relevant prototype to classify a test data set. But ultimately, the
classification is performed with tools such as tclust B.Dost et al. (2010), or
Quadratic Discriminant Analysis T.Hastie et al. (2009) that do not belong to
the field of optimal transport. On the contrary, CytOpT only requires one pre-
gated data set and the estimator of the class proportions is fully based on
regularized optimal transport. In this comparison, OptimalFlow had a learning
database composed of the nine Stanford cytometry data sets, and it was set to
produce one prototype from the database (Note that OptimalFlow thus leverages
the information from multiple source data sets, not just one). For the rest of
the method, we applied OptimalFlow with the default settings proposed in the
vignette222 at https://github.com/HristoInouzhe/optimalFlow.

Figure 9 displays a Bland-Altman plot for each of the compared method.
flowMeans and Cytometree clearly stand-out with less accurate class
proportions estimation that cannot compete with the good performance of CytOpT
and OptimalFlow. Although visually CytOpT and OptimalFlow seem to yield close
results, one can notice that OptimalFlow estimations are systematically biased
as it underestimates the rarest cell populations proportions while
overestimating the most represented cell populations frequencies.

To get a more thorough comparison between CytOpT and Optimalflow, we used the
Kullback-Leibler divergence to assess the discrepancy between the gold-
standard class proportions derived manually and the estimation for the two
methods. Figure 10 shows that on almost every data set the quality of
estimation is rather close. However when the target data set comes from Miami,
the Kullback-Leibler divergence between the benchmark and CytOpT estimations
is significantly lower than the Kullback-Leibler divergence between the
benchmark and OptimalFlow estimations. This means that in theses cases CytOpT
performs significantly better than OptimalFlow. We account for theses results
by the spatial shift between the source and the target. This spatial shifting
is displayed on Figure 11. To end this comparison, we present in Figure 12 a
comparison of the class proportions estimation between CytOpT and OptimalFlow
in two cases. The first case is the worst estimation for Optimalflow, namely,
when the target data set is Miami3C. In this case, OptimalFlow overestimates
very much the proportion of CD4 effector and the proportion of CD4 activated.
The second case corresponds to the worst estimation for CytOpT, namely when
the target data set is Nhlbi3B. If OptimalFlow performs slightly better in
this case, CytOpT still offers an acceptable and competitive estimation.

Figure 9: Comparison of the proportions $\hat{\pi}$ estimated by four
automated methods and the manual benchmark $\pi$ on the HIPC database. By
plotting the difference against the mean of the results of two different
methods, the Bland-Altman plot allows to assess the agreement between an
automated method and the manual benchmark. Figure 10: Comparison between our
algorithm CytOpT and one of the state of the art automated method for
cytometry data analysis: OptimalFlow. The comparison was performed on the data
sets of the HIPC panel. Figure 11: 2D projection of two cytometry data sets.
Left: measurements performed in Stanford. Right: measurements performed in
Miami. On this 2D projection one can assess the spatial shift induced by the
technical variability of cytometry. Figure 12: Comparison of the estimated
proportions $\hat{\pi}$ by CytOpT and OptimalFlow with the manual gating
benchmark $\pi$. Left: The target data set is Miami3C. Right: The target data
set is Nhlbi3B. In both cases, the data set Stanford1A has been used as a
source data set for CytOpT.

###  5.3 Application to the OptimalFlow data

A second evaluation of our method was performed on the cytometry data used in
Barrio et al. (2019). This second cytometry panel is available on GitHub 333
at https://github.com/HristoInouzhe/optimalFlowData, and we will refer to this
panel as the OptimalFlow data sets. First, in order to run a sensible
comparison between CytOpT, OptimalFlow and FlowMeans we selected nine cell
populations that were present in the first 21 data sets. As FlowMeans is an
unsupervised method, it was directly applied to the 21 data sets. For
OptimalFlow, it was set to use the 3 first data sets of the database, to build
one template, and then to apply the classification step. Hence OptimalFlow
yielded a class proportions estimation for the 18 remaining data sets. For
CytOpT, we used the first data set as a source data set and we estimated the
class proportions in the 20 remaining data sets. Figure 13 shows the Bland-
Altman plots of the results for those three methods. One can notice that
CytOpT is the method where the estimation of the class proportions
$\hat{\pi}_{k}$ is overall the closest to the manual gold-standard $\pi_{k}$.

Figure 13: Comparison of the proportions $\hat{\pi}$ estimated with CytOpT and
the manual benchmark $\pi$ on the OptimalFlow database.

##  6 Discussion

We have introduced CytOpT, a supervised method that directly estimates the
cell class proportions, without a preliminary clustering or classification
step. While obtaining a classification of the observation is often a
mathematical tool for achieving automatic gating, we emphasize that from a
clinical perspective it is the proportions of the different cell populations
that matter. Indeed, this is the key parameter which informs on a patient’s
biological and health condition.

To further improve the performance of CytOpT, an additional pre-processing
step to handle outliers could make CytOpT more robust to extreme cellular
observations that have little biological meaning. Indeed, due to the needed
scaling standardization of the data, one outlier observation can change the
layout of the source data in comparison with the layout of the target data. A
data driven strategy to choose the regularization parameter $\varepsilon$
could also help to tackle this outlier issue.

Following the supervised approach in Barrio et al. (2019) one could also think
of an extension of CytOpT that would leverage the information from several
labelled source data sets. A potential approach would be to select the most
relevant labelled data set to estimate the class proportions in an
unclassified target data set. To be more specific, let us assume that we have
$M$ classified data sets $X^{1},...,X^{M}$. For each data set $X^{m}$, we can
define the re-weighted empirical measure $\alpha^{m}(h)$ such that the weight
of the $k^{th}$ class equals $h_{k}$ for all $k\in\\{1,...,K\\}$, which would
yield an optimal re-weighting $\hat{\pi}^{m}$ that corresponds to an
estimation of the class proportions $\pi$ in the target data set. Finally, to
select the most relevant estimator $\hat{\pi}^{m}$, we choose the one leading
to the smallest Wasserstein cost between the re-weighted empirical measure
$\hat{\alpha}^{m}(\hat{\pi}^{m})$ and the empirical target distribution
$\hat{\beta}$ that is :

$\hat{\pi}=\operatorname*{arg\,min}_{\hat{\pi}^{1},...,\hat{\pi}^{M}}\\{W^{\varepsilon}(\hat{\alpha}^{1}(\hat{\pi}^{1}),\hat{\beta}),...,W^{\varepsilon}(\hat{\alpha}^{M}(\hat{\pi}^{M}),\hat{\beta})\\}$
(6.1)

Another extension of this work could also consider the Sinkhorn divergence
$S^{\varepsilon}(\alpha,\beta)=$
$W^{\varepsilon}(\alpha,\beta)-\frac{1}{2}W^{\varepsilon}(\alpha,\alpha)-\frac{1}{2}W^{\varepsilon}(\beta,\beta)$,
instead of $W^{\varepsilon}(\alpha,\beta)$. Indeed, the authors of J.Feydy et
al. (2018) showed that $S^{\varepsilon}$ has good theoretical properties that
make it a suitable loss for applications in machine learning. Future research
may also include a theoretical study of the convergence rates of the
stochastic algorithms described in Section 3 that aim at solving optimization
problem (3.3). Indeed, an improved optimization strategy could lead to
practical advantages and alleviate the computational burden of CytOpT.
Similarly, an appropriate adaptive decreasing step policy could lessen the
variations of the results induced by the stochastic nature of our approach.

In conclusion, the proposed algorithm is very promising in terms of its
performances in difficult situations where the technical variability of flow
cytometry induces spatial shifts between samples. We think this new strategy
is relevant for clinical applications, given the constant increase of the
number of cellular markers used in cytometry analysis. This raises the
question of the impact of such higher-dimensions, both for manual gating
(dramatically increasing the resources necessary) and automated approaches
(impacting the computational cost and efficiency). Regarding CytOpT, the
increase in the number of markers does not extend the computational cost by
much as the dimension only impacts the initial computation of the distance
matrix between the observations. We argue that with higher-dimensional data,
our method could perform even better as the Wasserstein distance would take
advantage of all the dimensions available at once.

## Codes and data availability

The illustrative HIPC data sets described in Subsection 1.1, as well as Python
Notebooks to reproduce the Figures presented in this paper are available at
https://github.com/Paul-Freulon/CytOpt.

## References

  * N.Aghaeepour et al. [2013] N.Aghaeepour, G.Finak, H.Hoos, T.R.Mosmann, R.Brinkman, R.Gottardo, R.H.Scheuermann, FlowCAP Consortium, Dream Consortium, et al.  Critical assessment of automated flow cytometry data analysis techniques.  _Nature methods_ , 10(3):228, 2013.
  * Y.Saeys et al. [2016] Y.Saeys, S.Van Gassen, and B.N.Lambrecht.  Computational flow cytometry: helping to make sense of high-dimensional immunology data.  _Nature Reviews Immunology_ , 16(7):449, 2016\.
  * Aghaeepour et al. [2011] Nima Aghaeepour, Radina Nikolic, Holger H Hoos, and Ryan R Brinkman.  Rapid cell population identification in flow cytometry data.  _Cytometry Part A_ , 79(1):6–13, 2011.
  * D.Commenges et al. [2018] D.Commenges, C.Alkhassim, R.Gottardo, B.P.Hejblum, and R.Thiébaut.  cytometree: a binary tree algorithm for automatic gating in cytometry analysis.  _Cytometry Part A_ , 93(11):1132–1140, 2018.
  * Y.Ge and S.C.Sealfon [2012] Y.Ge and S.C.Sealfon.  flowPeaks: a fast unsupervised clustering for flow cytometry data via K-means and density peak finding.  _Bioinformatics_ , 28(15):2052–2058, 05 2012\.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/bts300.  URL https://doi.org/10.1093/bioinformatics/bts300.
  * B.P.Hejblum et al. [2019] B.P.Hejblum, C.Alkhassim, R.Gottardo, F.Caron, R.Thiébaut, et al.  Sequential dirichlet process mixtures of multivariate skew $t$-distributions for model-based clustering of flow cytometry data.  _The Annals of Applied Statistics_ , 13(1):638–660, 2019.
  * H.Li et al. [2017] H.Li, U.Shaham, K.P.Stanton, Y.Yao, R.R.Montgomery, and Y.Kluger.  Gating mass cytometry data by deep learning.  _Bioinformatics_ , 33(21):3423–3430, 07 2017\.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btx448.  URL https://doi.org/10.1093/bioinformatics/btx448.
  * M.Lux et al. [2018] M.Lux, R.R.Brinkman, C.Chauve, A.Laing, A.Lorenc, L.Abeler-Dörner, and B.Hammer.  flowLearn: fast and precise identification and quality checking of cell populations in flow cytometry.  _Bioinformatics_ , 34(13):2245–2253, 02 2018\.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/bty082.  URL https://doi.org/10.1093/bioinformatics/bty082.
  * Barrio et al. [2019] E.del Barrio, H.Inouzhe, JM.Loubes, C.Matrán, and A.Mayo-Íscar.  optimalflow: Optimal-transport approach to flow cytometry gating and population matching.  _arXiv preprint arXiv:1907.08006_ , 2019.
  * R.Gottardo et al. [2014] V.Brusic R.Gottardo, S.H.Kleinstein, M.M.Davis, D.A.Hafler, H.Quill, A.K.Palucka, G.A.Poland, B.Pulendran, E.L.Reinherz, et al.  Computational resources for high-dimensional immune analysis from the human immunology project consortium.  _Nature biotechnology_ , 32(2):146, 2014.
  * H.Janati et al. [2018] H.Janati, M.Cuturi, and A.Gramfort.  Wasserstein regularization for sparse multi-task regression.  _arXiv preprint arXiv:1805.07833_ , 2018.
  * R.Flamary et al. [2018] R.Flamary, M.Cuturi, N.Courty, and A.Rakotomamonjy.  Wasserstein Discriminant Analysis.  _Machine Learning_ , 107(12):1923–1945, December 2018.  doi: 10.1007/s10994-018-5717-1.  URL https://hal.archives-ouvertes.fr/hal-01377528.
  * M.Arjovsky et al. [2017] M.Arjovsky, S.Chintala, and L.Bottou.  Wasserstein gan.  _arXiv preprint arXiv:1701.07875_ , 2017.
  * J.Solomon et al. [2015] J.Solomon, F.de Goes, G.Peyre, M.Cuturi, A.Butscher, A.Nguyen, T.Du, and L.Guibas.  Convolutional wasserstein distances: Efficient optimal transportation on geometric domains.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–11, 2015.
  * A.Genevay et al. [2016] A.Genevay, M.Cuturi, G.Peyré, and F.Bach.  Stochastic optimization for large-scale optimal transport.  In _Advances in neural information processing systems_ , pages 3440–3448, 2016.
  * B.Bercu and J.Bigot [2020] B.Bercu and J.Bigot.  Asymptotic distribution and convergence rates of stochastic algorithms for entropic optimal transportation between probability measures.  _Annals of Statistics_ , To be published, 2020.
  * M.Cuturi [2013] M.Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in neural information processing systems_ , pages 2292–2300, 2013.
  * G.Peyré et al. [2019] G.Peyré, M.Cuturi, et al.  Computational optimal transport.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * H.Robbins and S.Monro [1951] H.Robbins and S.Monro.  A stochastic approximation method.  _The annals of mathematical statistics_ , pages 400–407, 1951.
  * I.Redko et al. [2018] I.Redko, N.Courty, R.Flamary, and D.Tuia.  Optimal transport for multi-source domain adaptation under target shift.  _arXiv preprint arXiv:1803.04899_ , 2018.
  * JM.Bland and DG.Altman [1986] JM.Bland and DG.Altman.  Statistical methods for assessing agreement between two methods of clinical measurement.  _The lancet_ , 327(8476):307–310, 1986.
  * M.Ballu et al. [2020] M.Ballu, Q.Berthet, and F.Bach.  Stochastic optimization for regularized wasserstein estimators.  _arXiv preprint arXiv:2002.08695_ , 2020.
  * Tung et al. [2007] James W Tung, Kartoosh Heydari, Rabin Tirouvanziam, Bita Sahaf, David R Parks, Leonard A Herzenberg, and Leonore A Herzenberg.  Modern flow cytometry: a practical approach.  _Clinics in laboratory medicine_ , 27(3):453–468, 2007.
  * B.Dost et al. [2010] B.Dost, C.Wu, A.Su, and V.Bafna.  Tclust: A fast method for clustering genome-scale expression data.  _IEEE/ACM transactions on computational biology and bioinformatics_ , 8(3):808–818, 2010.
  * T.Hastie et al. [2009] T.Hastie, R.Tibshirani, and J.Friedman.  _The elements of statistical learning: data mining, inference, and prediction_.  Springer Science & Business Media, 2009.
  * J.Feydy et al. [2018] J.Feydy, T.Séjourné, FX.Vialard, SI.Amari, A.Trouvé, and G.Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  _arXiv preprint arXiv:1810.08278_ , 2018.
  * J.M.Borwein and D.Zhuang [1986] J.M.Borwein and D.Zhuang.  On fan’s minimax theorem.  _Mathematical programming_ , 34(2):232–234, 1986\.
  * C.Boyer et al. [2014] C.Boyer, P.Weiss, and J.Bigot.  An algorithm for variable density sampling with block-constrained acquisition.  _SIAM Journal on Imaging Sciences_ , 7(2):1080–1107, 2014.

##  Appendix A Two different numerical optimization schemes

We will now present the technical details of the optimization strategies
discussed in Section 3.4.

###  A.1 Descent Ascent procedure

In the optimization problem (3.3) that yields an estimator $\hat{\pi}$ of the
class proportions, the weights of the source distribution are modified, but
the support is fixed. Therefore, the parameterized measure $\alpha(h)$ is
entirely represented through the mapping $a:h\mapsto a(h)$. Hence, we define
the linear operator $\Gamma\in\mathbb{R}^{I\times K}$ that allows to express
the transformation from the class proportions to the weight vector $a(h)$ as:

$\forall(i,k)\in\\{1,...,I\\}\times\\{1,...,K\\},\leavevmode\nobreak\
\Gamma_{i,k}=\left\\{\begin{array}[]{ll}\frac{1}{n_{k}}&\mbox{if }X_{i}^{s}\in
C_{k},\\\ 0&\mbox{otherwise}.\end{array}\right.$ (A.1)

Thanks to this transformation, the vector $a(h)=\Gamma h$ is such that the
class proportions in the discrete measure
$\sum_{i=1}^{I}a_{i}(h)\delta_{X_{i}^{s}}$ equals the vector $h$. To present
the ideas of our first way to numerically solve
$\eqref{eq:estimator_proportion}$, we re-write this optimization problem as:

$\min_{h\in\Sigma_{K}}W^{\varepsilon}(\alpha(h),\beta)=\min_{h\in\Sigma_{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,h)],$
(A.2)

where $Y$ is a random variable with distribution $\beta$, and for
$y_{j}\in\mathbb{R}^{d}$ an observation of the random variable $Y$,

$g_{\varepsilon}(y_{j},u,h)=\sum_{i=1}^{I}u_{i}a_{i}(h)+\varepsilon\left(\log(b_{j})-\log\left(\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)\right)\right)-\varepsilon$
(A.3)

We propose a gradient descent ascent algorithm in order to approximate
$\hat{\pi}$ a solution of the optimization problem
$\eqref{eq:min_max_expectation_appendix}$. At each iteration, our method runs
multiple stochastic gradient ascent steps to estimate the solution of the
inner maximization problem. Thanks to this estimation, we can approximate the
gradient of the function $h\mapsto W(\alpha(h),\beta)$ that we wish to
minimize.

We will now detail the optimization procedure. First, in order to avoid
projecting $h_{n}$ onto the simplex $\Sigma_{K}$ at each step of the
procedure, we re-parameterized our problem with a soft-max function. Our new
optimization problem is now:

$\min_{z\in\mathbb{R}^{K}}W^{\varepsilon}(\Gamma(\sigma(z)),b)=\min_{z\in\mathbb{R}^{K}}F(z),$
(A.4)

where $F(z)=W^{\varepsilon}(\Gamma(\sigma(z)),b)$, and
$\sigma:\mathbb{R}^{K}\rightarrow\Sigma_{K}$ is the soft-max function defined
as

$\sigma(z)_{l}=\frac{\exp(z_{l})}{\sum_{k=1}^{K}\exp(z_{k})}.$ (A.5)

Then, if we denote $\hat{z}$ a minimizer of (A.4), we will derive an estimator
of $\pi$ as $\hat{\pi}=\sigma(\hat{z}).$

Our approach uses the fact that the gradient of the function $a\mapsto
W^{\varepsilon}(a,b)$ is given by:

$\frac{\partial}{\partial a}W^{\varepsilon}(a,b)=u^{*}$ (A.6)

where $u^{*}$ is the unique solution of (2.3) centered such that
$\sum_{i=1}^{I}u_{i}^{*}=0$. This result is established in Proposition 4.6 of
G.Peyré et al. [2019]. Hence, by applying the chain rule of differentiation,
one obtains that:

$\nabla_{z}W^{\varepsilon}(\Gamma(\sigma(z)),b)=(\Gamma
J_{\sigma}(z))^{T}u_{z}^{*},$ (A.7)

where $\Gamma$ is the linear operator defined in (A.1), $u_{z}^{*}$ denotes
the maximiser of (2.3) when the weights of the distribution $\alpha$ equal
$a=\Gamma\sigma(z)$ and $J_{\sigma}(z)$ is the Jacobian matrix of $\sigma$,
defined for $z\in\mathbb{R}^{K}$ by:

$\forall(i,j)\in\\{1,...,K\\}^{2},\leavevmode\nobreak\
J_{\sigma}(z)_{i,j}=\left\\{\begin{array}[]{lll}\frac{\exp(x_{j})\left(\sum_{k=1}^{K}\exp(x_{k})\right)-\exp(2x_{j})}{\left(\sum_{k=1}^{K}\exp(x_{k})\right)^{2}}&\mbox{if
}i=j,\\\ \\\
-\frac{\exp(x_{i}+x_{j})}{\left(\sum_{k=1}^{K}\exp(x_{k})\right)}&\mbox{otherwise.}\end{array}\right.$
(A.8)

The stochastic algorithm that we consider to find a minimizer $z^{*}$ of
${F:z\mapsto W^{\varepsilon}(\Gamma\sigma(z),b)}$ is then given by the
recursive procedure

$\hat{z}_{n+1}=\hat{z}_{n}-\eta\hat{\omega}(\hat{z}_{n}),$ (A.9)

where $\hat{\omega}(\hat{z}_{n})$ stands for a stochastic approximation of
$\nabla_{z}W^{\varepsilon}(\Gamma(\sigma(\hat{z}_{n})),b)$. We propose to
estimate $\nabla_{z}W^{\varepsilon}(\Gamma(\sigma(\hat{z}_{n})),b)$ by:

$\hat{\omega}(\hat{z}_{n})=(\Gamma
J_{\sigma}(\hat{z}_{n}))^{T}\widehat{U}_{m+1}^{(n+1)}.$ (A.10)

That is, we plug an estimate of $u_{\hat{z}_{n}}^{*}$ in the formula (A.7) of
the gradient of the objective function when the weights of the source
distribution equal $\Gamma(\sigma(\hat{z}_{n}))$.

The estimate $\widehat{U}_{m+1}^{(n+1)}$ of $u_{\hat{z}_{n}}^{*}$ is
calculated by the Robbins-Monro Algorithm H.Robbins and S.Monro [1951] which
is defined through the recurrence:

$\widehat{U}_{k+1}^{(n+1)}=\widehat{U}_{k}^{(n+1)}+\gamma_{k+1}\nabla_{u}g_{\varepsilon}(Y_{k+1}^{(n+1)},\widehat{U}_{k}^{(n+1)},\hat{z}_{n}),$
(A.11)

where $\widehat{U}_{0}^{(n+1)}$ is arbitrarily chosen such that
$\langle\widehat{U}_{0}^{(n+1)},1_{I}\rangle=0$, and
$Y_{1}^{(n+1)},...,Y_{m+1}^{(n+1)}$ are i.i.d random variables sampled from
$\beta$ at the iteration to pass from $\hat{z}_{n}$ to $\hat{z}_{n+1}$, and
$(\gamma_{n})_{n\geq 0}$ is a positive sequence of real numbers decreasing
toward zero satisfying condition (2.7).

Once our minimization procedure is over and we have computed $\hat{z}$ a
minimizer of (A.4), we compute an estimate of the class proportions in the
target observations. To do so, we set

$\hat{\pi}=\sigma(\hat{z}).$ (A.12)

##### Choice of the parameters for the descent-ascent procedure

For this Descent-Ascent procedure, we need to set several parameters. For the
step size policy $(\gamma_{k})_{k>0}$ of the inner loop, we rely on the
recommendation done in B.Bercu and J.Bigot [2020]. Therefore, we chose
$\gamma_{k}=\gamma/n^{c}$ where $c=0.51$ and $\gamma=J\varepsilon/1.9$ with
$J$ the cardinal of $\beta$. We set the other parameters experimentally. Thus,
we chose $n_{out}=10000$, $n_{in}=10$, $\eta=10$ and $\varepsilon=0.0005$.

###  A.2 Minmax swapping procedure

In this section we use the ideas of M.Ballu et al. [2020]. To propose an
alternative scheme to solve the minimization problem
$\eqref{eq:estimator_proportion}$ we slightly modify this problem by adding an
entropic term

$\varphi(h)=\sum_{k=1}^{K}h_{k}\log(h_{k}).$ (A.13)

Thus, our new problem is now:

$\min_{h\in\Sigma}W^{\varepsilon}(\alpha(h),\beta)+\lambda\varphi(h)=\min_{h\in\Sigma_{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,h)]+\lambda\varphi(h),$
(A.14)

where $Y$ is a random variable with distribution $\beta$, and
$g_{\varepsilon}$ is defined in (A.3). By swapping the minimum and the maximum
according to Fan’s minimax theorem J.M.Borwein and D.Zhuang [1986] we get

$\begin{split}\min_{h\in\Sigma_{K}}\max_{u\in\mathbb{R}^{I}}\mathbb{E}[g_{\varepsilon}(Y,u,h)]+\lambda\varphi(h)&=\max_{u\in\mathbb{R}^{I}}\min_{h\in\Sigma_{K}}\mathbb{E}[g_{\varepsilon}(Y,u,h)]+\lambda\varphi(h)\\\
&=\max_{u\in\mathbb{R}^{I}}\min_{h\in\Sigma_{K}}\sum_{i=1}^{I}u_{i}(\Gamma
h)_{i}+A(u)-\varepsilon+\lambda\varphi(h)\\\
&=\max_{u\in\mathbb{R}^{I}}A(u)+\min_{h\in\Sigma_{K}}\sum_{i=1}^{I}u_{i}(\Gamma
h)_{i}+\lambda\varphi(h)-\varepsilon\\\
&=\max_{u\in\mathbb{R}^{I}}A(u)+J_{\lambda}(u)-\varepsilon,\end{split}$ (A.15)

where
$A(u)=\varepsilon\sum_{j=1}^{J}\left(\log(b_{j})-\log\left(\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)\right)\right)b_{j}$
and

$J_{\lambda}(u)=\min_{h\in\Sigma_{K}}u_{i}(\Gamma h)_{i}+\lambda\varphi(h).$
(A.16)

Arguing e.g. as in the proof of [C.Boyer et al., 2014, Proposition 4.1], the
solution $h(u)$ of (A.16) is defined by

$\forall k\in\\{1,...,K\\},\leavevmode\nobreak\
(h(u))_{k}=\frac{\exp\left(-\frac{(\Gamma^{T}u)_{k}}{\lambda}\right)}{\sum_{l=1}^{K}\exp\left(-\frac{(\Gamma^{T}u)_{l}}{\lambda}\right)}.$
(A.17)

By plugging (A.17) in (A.15), Problem (A.14) boils down to the following
maximization problem with respect to $u\in\mathbb{R}^{d}$:

$\max_{u\in\mathbb{R}^{I}}\varepsilon\sum_{j=1}^{J}\left(\log(b_{j})-\log\left(\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)\right)\right)b_{j}-\lambda\log\left(\sum_{l=1}^{K}\exp\left(-\frac{(\Gamma^{T}u)_{l}}{\lambda}\right)\right)-\varepsilon,$
(A.18)

which can be rewritten

$\max_{u\in\mathbb{R}^{I}}\mathbb{E}[f_{\varepsilon,\lambda}(Y,u)].$ (A.19)

Here, $Y$ is a random variable with distribution $\beta$, and for
$y_{j}\in\mathbb{R}^{d}$ an observation of $Y$, and $u\in\mathbb{R}^{I}$,

$f_{\varepsilon,\lambda}(y_{j},u)=\varepsilon\left(\log(b_{j})-\log\left(\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)\right)\right)-\lambda\log\left(\sum_{l=1}^{K}\exp\left(-\frac{(\Gamma^{T}u)_{l}}{\lambda}\right)\right)-\varepsilon.$
(A.20)

The fact that for all $y\in\mathbb{R}^{d}$ the function
$f_{\varepsilon,\lambda}(y,.)$ is concave, and the expectation form (A.19) of
problem (A.18) lead us to estimate the optimal vector $u^{*}$ by the Robbins-
Monro algorithm given, for all $n\geq 0$ by

$\widehat{U}_{n+1}=\widehat{U}_{n}+\gamma_{n+1}\nabla_{u}f_{\varepsilon,\lambda}(Y_{n+1},\widehat{U}_{n}),$
(A.21)

where the initial value $\hat{U}_{0}$ is a random vector which can be
arbitrarily chosen, $Y_{1},...,Y_{n+1}$ are i.i.d random variables sampled
from $\beta$ ,and $(\gamma_{n})_{n\geq 0}$ is a positive sequence of real
numbers decreasing toward zero satisfying condition (2.7).

For $y_{j}\in\mathbb{R}^{d}$ an observation of $Y$, and $u\in\mathbb{R}^{I}$,
the gradient of $\nabla_{u}f_{\varepsilon,\lambda}(y_{j},u)$ is defined by
$\forall i_{0}\in\\{1,...,I\\}$,

$(\nabla_{u}f_{\varepsilon,\lambda}(y_{j},u))_{i_{0}}=\frac{\sum_{l=1}^{K}\Gamma_{i_{0},l}\exp\left(-\frac{(\Gamma^{T}u)_{l}}{\lambda}\right)}{\sum_{k=1}^{K}\exp\left(-\frac{(\Gamma^{T}u)_{k}}{\lambda}\right)}-\frac{\exp\left(\frac{u_{i_{0}}-c(x_{i_{0}},y_{j})}{\varepsilon}\right)}{\sum_{i=1}^{I}\exp\left(\frac{u_{i}-c(x_{i},y_{j})}{\varepsilon}\right)}.$
(A.22)

Once the algorithm has converged, and we get $\hat{U}$, a satisfactory
approximation of a maximizer of problem A.19, one can compute an estimate of
the class proportions by setting

$\hat{\pi}=h(\hat{U}),$ (A.23)

where $h(u)$ is defined in A.17.

##### Choice of the parameters for the min-max swapping procedure

For this second procedure we have to set several parameters. To satisfy
condition (2.7) for the step size policy $(\gamma_{n})_{n>0}$ of the Robbins-
Monro algorithm, we chose $\gamma_{n}=\gamma/n^{c}$ where $\gamma=5$ and
$c=0.99$. The other parameters have been set empirically. We chose
$\lambda=0.0001$ and $\varepsilon=0.0001$. The number of iterations for this
stochastic ascent algorithm has been set to $n=10000$.