# Schema matching using Gaussian mixture models with Wasserstein distance

Mateusz Przyborowski
University of Warsaw, Warsaw, Poland
QED Software, Warsaw, Poland
<EMAIL_ADDRESS>
Mateusz Pabiś
QED Software, Warsaw, Poland
<EMAIL_ADDRESS>
Andrzej Janusz
University of Warsaw, Warsaw, Poland
QED Software, Warsaw, Poland
<EMAIL_ADDRESS>
Dominik Ślęzak
University of Warsaw, Warsaw, Poland
QED Software, Warsaw, Poland
<EMAIL_ADDRESS>

###### Abstract

Gaussian mixture models find their place as a powerful tool, mostly in the
clustering problem, but with proper preparation also in feature extraction,
pattern recognition, image segmentation and in general machine learning. When
faced with the problem of schema matching, different mixture models computed
on different pieces of data can maintain crucial information about the
structure of the dataset. In order to measure or compare results from mixture
models, the Wasserstein distance can be very useful, however it is not easy to
calculate for mixture distributions. In this paper we derive one of possible
approximations for the Wasserstein distance between Gaussian mixture models
and reduce it to linear problem. Furthermore, application examples concerning
real world data are shown.

_K_ eywords Gaussian mixture models $\cdot$ Wasserstein distance $\cdot$
schema matching $\cdot$ big data

##  1 Introduction

Mixture model is a probabilistic model that is able to infer subpopulations
from total population without additional information (within the paradigm of
unsupervised learning). Mixture models closely correspond to the mixture
distributions of the probabilistic distributions of observations. In general,
in the structure of mixture model, we make assumptions over latent variables
that evaluate membership of each observation. Given the dataset, we can assume
that it is a sample and then mixture model can estimate the parameters of the
probability distributions that created points of this dataset, as well as
assign each observation vector of probabilities indicating the original
distribution.
Comparing different mixture models can be considered a generalization of the
problem of comparing different distributions. From the viewpoint of optimal
transport theory, the Wasserstein distance is an important method for
measuring similarities and the maintenance of the explainable nature of
mixture models.
In this paper we derive one of possible approximations of Wasserstein
distances computed between mixture models, which may be reduced to linear
optimization problem, and we present examples of usage

##  2 Related work

Gaussian mixture models with Wasserstein distance find their place in many
areas of machine learning. In case of generative networks[3], the use of
Wasserstein distance has been proved to model more complex distributions.
Autoencoder architectures equipped with Wasserstein distance (WAE), unlike
variational autonencoders (VAE), allow to use deterministic mapping to a
latent space[7]. In image processing, Gaussian mixture models equipped with
Wasserstein distance proved to be useful in tasks of color transfer and
texture synthesis[2]. When dealing with heterogeneous data, mixture models
have the advantage of simplicity and Wasserstein distance provides a suitable
convergence rate[5]. Moreover, Wasserstein distance holds an important place
in optimal transport theory[6][8].

##  3 Problem formulation

Let $p(x|z)$ be a probability distribution of the given data with unknown
vector of parameters $z$. Modeling the data using statistics and machine
learning comes to modeling probability distribution. In real world
applications, data is usually composed of multiple different probability
distributions. Hence comes an elementary idea of modeling the data using
mixture model, where each observation is assigned a probability of originating
from the given probability distribution. The problem of choosing types of
probability distributions for each component is usually skipped by assuming
normality (Gaussian) of individual components, as normal distribution has
important probabilistic properties. This approach is focused on a general
summary of the very origin of the data, therefore its applications are
widespread:

  1. 1.

in cluster analysis, Gaussian mixture models (GMM) may be seen as an extension
to K-means algorithm, yielding additional information about given
observations;

  2. 2.

in supervised learning, associating a type of labels from training data with
one or more components may give us a similarity function between observations,
based on whether they originate from the same probability distribution;

  3. 3.

in natural language processing, distribution of words in documents can be
modelled as mixture of different categorical distributions.

###  3.1 Big data

Nowadays dealing with the big data is a popular issue. While focusing on a big
volume of moderately dimensional data, mixture models can help with summary of
the most common type of observations. Suppose that size of the data makes it
impractical to repetitively perform calculations using the entire data. If we
could summarize the data by creating representations, which allow to maintain
most important features of the data, as well as to perform calculations
yielding approximate but much faster solutions, we would save a lot of
computing power and time in practical applications. Mixture models may be
considered as one of such approaches, in which the data representation is made
of the components understood as parameters of probability distributions. The
mixture model of a given dataset is itself an approximation of the underlying
probability distribution. While it gives a way to compare different
observations from the same dataset, one may think about comparing different
representations, i.e. different mixture models. Suppose that we split a
labeled dataset into datasets based on label, then compute mixture model for
each of such datasets. Under the assumption that different labels indicate a
different distribution of features, comparing the mixture models allows us to
conclude that two datasets originate from similar sources. This problem is
more widely known as schema matching problem and is a common task in data
integration and database management.

###  3.2 Comparing mixture models

Mixture models, by the very way they are calculated, are based on the values
of many observations. The only difference between the resulting models must be
a manifestation of the different values of the respective observations. This
interpretation yields a corollary that the difference between models could be
measured by how much and how many observations making up one mixture model
must be transformed in order for the resulting model to be more similar to the
one with which it is compared. This intuition is realized in the Wasserstein
metric, where the distance between two probability distributions is the amount
of ”work” that needs to be done in order to transform one distribution into
another. Further explanation is provided in the following sections.
Gaussian mixture models allow us to summarize large datasets, while
Wasserstein distance makes a tool for comparing different representations.

##  4 Gaussian mixture models

Henceforth we will focus on the Gaussian mixture models, i.e. mixture models
only with normal components.

######  Definition 1.

Let $n\in\mathbb{N}_{+}$, $w=[w_{1},...,w_{n}]$ s.t. $\forall i\leq n\quad
0\leq w_{i}\leq 1$ and $\sum_{i=1}^{n}w_{i}=1$. A Gaussian mixture model of
size $n$ components is a probability distribution $p$ defined as:

$p=\sum_{i=1}^{n}w_{i}\mathcal{N}(\mu_{i},\Sigma_{i})$

where $\mathcal{N}(\mu_{i},\Sigma_{i})$ is a normal distribution with
$\mu_{i}$the vector of mean and $\Sigma_{i}$ the covariance matrix as
parameters corresponding to the $i$-th component.

Thereafter, if not stated otherwise, we will assume that distributions are
defined over $\mathbb{R}^{d}$ with the dimension $d$. For the sake of
simplicity, if it is not necessary we will omit the dimension.
Fitting Gaussian mixture model to a given data is a task of finding
appropriate values for parameters $\\{w_{i},\mu_{i},\Sigma_{i}\\}_{i=1}^{n}$
s.t. the resulting model describes the dataset. Using maximum likelihood
estimation (MLE) let $X=[X_{1},...,X_{k}]$ be a vector of $k$ observations
from our dataset. The joint probability distribution is then defined as:

$p(X_{1}=x_{1},...,X_{k}=x_{k})=\prod_{j=1}^{k}\sum_{i=1}^{n}w_{i}\mathcal{N}(x_{j}|\mu_{i},\Sigma_{i})$

Likelihood function is defined as:

$L(w,\mu,\Sigma|X)=\prod_{j=1}^{k}\sum_{i=1}^{n}w_{i}\mathcal{N}(x_{j}|\mu_{i},\Sigma_{i})$
$\log(L(w,\mu,\Sigma|X))=\sum_{j=1}^{k}\log(\sum_{i=1}^{n}w_{i}\mathcal{N}(x_{j}|\mu_{i},\Sigma_{i}))$

Unfortunately, differentiation and comparing to $0$ will not allow us to
analytically solve this equation. In order to help with this, we will
introduce a latent variables $z$ that explain which component generated given
observation. Then:

$p(X)=\sum_{i=1}^{n}p(x,z=i)=\sum_{i=1}^{n}p(z=i)p(x|z=i)=\sum_{i=1}^{n}w_{i}\mathcal{N}(\mu_{i},\Sigma_{i})$
$p(z_{j}=i_{0}|X_{j})=\frac{p(X_{j}|z_{j}=i_{0})p(z_{j}=i_{0})}{p(X_{j})}=\frac{w_{i_{0}}\mathcal{N}(\mu_{i_{0}},\Sigma_{i_{0}})}{\sum_{i=1}^{n}w_{i}\mathcal{N}(\mu_{i},\Sigma_{i})}=\gamma_{i_{0}}(z_{j})$

###  4.1 EM algorithm

We can notice that knowing either parameters or $\gamma$ allows us to compute
the missing part. Furthermore, having a random guess about parameters, we can
evaluate $\gamma$ probabilities and then estimate new parameters. Repeating
this process, as well as measuring progress with log-likelihood, is a sketch
of an iterative method known as the expectation–maximization (EM) algorithm.
Let $\theta=(w,\mu,\Sigma)$, then:

$Q(\theta,\theta_{0})=\mathbb{E}\log
p(X,Z|\theta)=\sum_{z}p(z|X,\theta^{0})\log p(X,z|\theta)$

We can simplify to this form:

$Q(\theta,\theta_{0})=\sum_{z}\gamma_{(}z)\log p(X,z|\theta)$

And since:

$\log
p(X,Z|\theta)=\log\prod_{j=1}^{k}\prod_{i=1}^{n}w_{i}^{z_{i,j}}\mathcal{N}(X|\mu_{i},\Sigma_{i})^{z_{i,j}}=\sum_{j=1}^{k}\sum_{i=1}^{n}z_{i,j}(\log(w_{i})+\log(\mathcal{N}(X|\mu_{i},\Sigma_{i})))$

In summary, we have:

$Q(\theta,\theta_{0})=\sum_{j=1}^{k}\sum_{i=1}^{n}\gamma_{i}(z_{j})(\log(w_{i})+\log(\mathcal{N}(X|\mu_{i},\Sigma_{i})))$

As for expectation phase (E), we evaluate $Q(\theta,\theta_{0})$ given initial
$\theta_{0}$. Maximization step (M) consists of solving $\theta_{1}=\arg
max_{\theta}Q(\theta,\theta_{0})$. Solving for $Q$ is performed using Lagrange
multipliers. These steps are repeated until the stop conditions are met.

Algorithm 1 EM algorithm

$\theta,\gamma(Z)$

$\theta_{0}\leftarrow$ initial guess

repeat

calculate $Q(\theta,\theta_{0})$

$\theta_{0}\leftarrow\arg max_{\theta}Q(\theta,\theta_{0})$

until stop condition satisfied

###  4.2 Bayesian Gaussian mixture models

Finding parameters for GMM with EM algorithm does not include particular
hyperparameters, as e.g. number of components. One can imagine that having $N$
observations and $N$ components in the form of Dirac delta functions would
perfectly model a dataset, yet it would not be useful. The number of
components can also be a very important parameter for the regularization of
overfitting, a phenomenon in which the model may not be able to generalize
outside of the training set. Bayesian interpretation allows us to use the
prior probability distribution (Dirichlet distribution) to model the parameter
space. Estimating the approximate posterior distribution over the parameters
of a Gaussian mixture distribution yields the number of components from the
dataset.

##  5 Wasserstein distance

######  Definition 2.

Let $p$ and $q$ be two $d$-dimensional probability distributions and let
$\Gamma(p,q)$ be a set of probability distributions whose marginals are $p$
and $q$ as the first and second factors respectively. Let $n\geq 1$. The
$n$-th Wasserstein distance $W_{n}$ between $p$ and $q$ is defined as:

$W_{n}(p,q)=\inf_{\nu\in\Gamma(p,q)}(\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}||x-y||^{n}d\nu(x,y))^{1/n}=\inf(\mathbb{E}(||X-Y||^{n}))^{1/n}$

We may notice that looking at the measures corresponding to distributions $p$
and $q$, the set $\Gamma(p,q)$ is compact in the sense of weak convergence,
therefore infimum is achievable.

######  Lemma 1.

Let $\gamma\in\Gamma(p,q)$ and $x,y\in supp(\gamma)$ be $1$-dimensional
($d=1$) elements from support of $\gamma$ s.t. $x_{1}<y_{1}$ and
$y_{2}<x_{2}$. Then the following inequality holds:

$|x_{1}-y_{2}|^{n}+|y_{1}-x_{2}|^{n}<|x_{1}-x_{2}|^{n}+|y_{1}-y_{2}|^{n}$

###### Proof.

We can notice that:

$\exists t\in(0,1)\quad tx_{1}+(1-t)x_{2}=ty_{1}+(1-t)y_{2}$

Since points $x_{1},y_{2},tx_{1}+(1-t)x_{2}$ lie on the same line:

$|x_{1}-y_{2}|=|x_{1}-tx_{1}-(1-t)x_{2}|+|tx_{1}+(1-t)x_{2}-y_{2}|=(1-t)|x_{1}-x_{2}|+t|y_{1}-y_{2}|$

From Jensen’s inequality we have:

$|x_{1}-y_{2}|^{n}\leq(1-t)^{n}|x_{1}-x_{2}|+t^{n}|y_{1}-y_{2}|$

Symmetrically we get the result:

$|x_{2}-y_{1}|^{n}\leq t^{n}|x_{1}-x_{2}|+(1-t)^{n}|y_{1}-y_{2}|$

Summing gives us the inequality. ∎

######  Theorem 1.

Let $P,Q$ be cumulative distribution functions of distributions $p$ and $q$,
by $P^{-1}$ and $Q^{-1}$ we mean inverse cumulative distribution functions or
quantile functions; for $d=1$ we have:

$W_{n}(p,q)=(\int_{0}^{1}|P^{-1}(t)-Q^{-1}(t)|^{n}dt)^{1/n}$

###### Proof.

If $\gamma\in\Gamma(p,q)$ satisfies infimum in the definition of Wasserstein
distance, then $x_{2}\leq y_{2}$. Otherwise, by lemma, it would mean that
there exists a better fit, where swapping $y_{1}$ and $y_{2}$ gives smaller
value.
Let $x\in supp(p),y\in supp(q)$; then we can notice that $(x,y)\in
supp(\gamma)\iff P(x)=Q(y)$. Indeed,

$(x,y)\in
supp(\gamma)\iff\gamma(\mathbb{R},(-\infty,y])=\gamma((-\infty,x],(-\infty,y])=\gamma((-\infty,x],\mathbb{R})$

Therefore we conclude that:

$\int_{\mathbb{R}\times\mathbb{R}}|x-y|^{n}d\gamma(x,y)=\int_{supp(\gamma)}|x-y|^{n}d\gamma(x,y)=\int_{0}^{1}|P^{-1}(t)-Q^{-1}(t)|^{n}dt$

∎

Thereafter, if not stated otherwise, we will consider $W:=W_{2}$, i.e.
Wasserstein distance for $n=2$.

###  5.1 Connections with transportation theory

While considering probability as a mass over some space, Wasserstein distance
realises the optimal transport problem for transforming one probability
distribution into another. Suppose we have a cost function $c$ and probability
distributions $p,q$. A transport plan is a function $\gamma$ s.t.
$\gamma(x,y)$ is a volume of mass that needs to be moved from $x$ to $y$. Cost
of a transport plan $\gamma$ is:

$\int\int\gamma(x,y)c(x,y)dxdy=\int c(x,y)d\gamma(x,y)$

Depending on the selection of the function $c$, going with infimum over
possible plans yields us cost of optimal transport.

##  6 Wasserstein distance between two Gaussian mixture models

In order to calculate Wasserstein distance between Gaussian mixture models, we
would need to calculate an inverse cumulative distribution function for
mixture of normal distributions. Since it is analytically impossible, a
similar idea is adopted.

######  Theorem 2.

[2]

$W_{2}^{2}(\mathcal{N}(\mu_{1},\Sigma_{1}),\mathcal{N}(\mu_{2},\Sigma_{2}))=||\mu_{1}-\mu_{2}||^{2}+tr(\Sigma_{1}+\Sigma_{2}-2(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2})^{1/2})$

######  Definition 3.

111Similar definition is proposed in [2], but we only consider the finite
case.

Let
$p_{1}=\sum_{i=1}^{n}w_{1i}\mathcal{N}(\mu_{1i},\Sigma_{1i}),p_{2}=\sum_{i=1}^{m}w_{2i}\mathcal{N}(\mu_{2i},\Sigma_{2i})$
be two GMMs. We define approximate Wasserstein distance between $p_{1}$ and
$p_{2}$ in the following way:

$\hat{W}(p_{1},q_{1})=min\\{\sum_{i=1}^{n}\sum_{j=1}^{m}t_{ij}W(\mathcal{N}(\mu_{1i},\Sigma_{1i}),\mathcal{N}(\mu_{2j},\Sigma_{2j}))\;|\;{t_{i,j}\geq
0:\forall i\sum_{j}t_{i,j}=w_{1i}\land\forall j\sum_{i}t_{i,j}=w_{2j}}\\}$

Proposed Wasserstein distance between mixture models is a straightforward
extension of intuitions lying behind original Wasserstein distance. From the
transport point of view, we are looking for the best assignment between
corresponding mixtures. This can be extended to infinite dimensional form,
i.e. when the number of components is not finite. The main difference is that
here we do not seek the best transportation plan understood as a function or a
measure, but a matrix of size $n\times m$.

###  6.1 Dual problem

Let us consider a more general problem; let $a_{i},b_{j},c_{ij},d_{ij}$ be
given nonnegative integers, the problem is following:

$minimize\quad\sum_{i=1}^{n}\sum_{j=1}^{m}d_{ij}x_{ij}$
$subject\>to:\quad\forall i=1,...,n\sum_{j=1}^{m}x_{ij}=a_{i}\land\forall
j=1,...,m\sum_{i=1}^{n}x_{ij}=b_{j}$ $\forall i,j\quad 0\leq x_{ij}\leq
c_{ij}$

In the case of Wasserstein distance $c_{ij}\equiv 1$, $d_{ij}$ is a
Wasserstein distance between $i$-th component from first mixture and $j$-th
component from second mixture, ${a_{i}}_{i=1}^{n}$ are weights of components
from first mixture and ${b_{j}}_{j=1}^{m}$ are weights of components from
second mixture.
The dual problem has the following form:

$maximize\quad\sum_{i=1}^{n}a_{i}\alpha_{i}+\sum_{j=1}^{m}b_{j}\beta_{j}+\sum_{i=1}^{n}\sum_{j=1}^{m}c_{ij}\gamma_{ij}$
$subject\>to:\quad\forall i=1,...,n\forall j=1,...,m\quad
d_{ij}\geq\alpha_{i}+\beta_{j}+\gamma_{ij}\land\gamma_{ij}\leq 0$

It is worth to notice that the dual form immediately yields us a possible
solution: setting all $\alpha_{i}=\beta_{j}=\gamma_{ij}=0$.

###  6.2 Solving with linear programming

Finding Wasserstein distance between two mixture models comes down to solving
particular transport problem. Therefore we can use the notations of graph
theory: with a directed complete bipartite graph we have a cost over each edge
being a Wasserstein distance between the given components, a capacity at each
edge corresponding to the weight of the component and amount of flow, i.e. the
value sought. We can use the network simplex algorithm to solve such a
problem.

###  6.3 GMM with Wasserstein distance as a classifier

We present the algorithm for classification problem using Gaussian mixture
models and Wasserstein distance.

Algorithm 2 GMM with Wasserstein distance

$\\{(x_{i},y_{i})\\}_{i=1}^{k}$ \- training dataset with $m$ different classes

$\\{U_{i}\\}_{i=1}^{n}$ \- $n$ test datasets; within dataset each observation
has the same label

1\. Split training dataset into $m$ sets $\\{Z_{j}\\}_{j=1}^{,}$ based on the
label

2\. Fit Gaussian mixture model $p_{j}$ for each $Z_{j}$

3\. Fit Gaussian mixture model $q_{i}$ for each $U_{i}$

4\. Compute Wasserstein between $p_{j}$ and $q_{i}$ for each $i$, $j$

5\. Label the set $U_{i}$ with a label of the set $Z_{j_{0}}$, where
$j_{0}=argmin_{j}\hat{W}(p_{j},q_{i})$

##  7 Experiments

###  7.1 STL-10 dataset

In the first experiment we used features extracted from an autoencoder neural
network, which was trained in the case of image recognition. Original dataset
is the STL-10 dataset[1], which consists of total $13000$ images labeled as
one of ten possible classes. Extracted representation has dimensionality of
$512$, therefore during experiment we randomly choose some subset of
dimensions. Received results has been compares with $L2$ distance and
quadratic Jensen-Rényi divergence.

Figure 1: OX axis indicates number of selected components. At each step
experiment was repeated $20$ times. Solid lines stands for mean results, while
shaded area indicates $\pm 1$ standard deviation from the mean.

The task of matching data types using proposed method, i.e. Gaussian mixture
models with Wasserstein distance, relies only on applied preprocessing. While
autoencoder representations can be summaries themselves[4], given the large
volume of the data, our method is far more practical.

###  7.2 Text data

In the second experiment we operated on a large volume ($13000000$) of short
text data, divided into chunks of examples with the same label. The task was
to predict label of entire chunk. Preprocessing consists of transforming
characters into features based on length and frequency of occurrence of given
letters and signs. During experiment we performed $5$ $2$-fold cross-
validations. Results were compared with a different approach (KNN algorithm)
and similar approach with different distance function ($L2$).

Figure 2: Mean results from $5$ $2$-fold cross-validations. Standard
deviations are indicated.

Results conclude that proposed framework works better than compared methods.

##  8 Conclusions and future work

We derived the approximate easy easy-to-calculate version of the Wasserstein
distance between Gaussian mixture models, that may find many applications in
various fields of machine learning. In the case of big data, the greatest
advantage is the avoidance of multiple calculations over the entire dataset,
as the obtained summary allow for the estimation of similarity based only on
the compacted data representations.
The future work may involve a statistical analysis of the properties extracted
by Gaussian mixture models from a dataset, e.g. selecting important
observations that may have had the greatest impact on the parameters.

## References

  * [1] A. Coates, H. Lee, A. Y. Ng, ”An Analysis of Single Layer Networks in Unsupervised Feature Learning”, AISTATS, 2011.
  * [2] J. Delon and A. Desolneux, ”A Wasserstein-type distance in the space of Gaussian Mixture Models”, arXiv:1907.05254v4, 2019.
  * [3] B. Gaujac, I. Feige, and D. Barber, ”Gaussian mixture models with Wasserstein distance,” arXiv:1806.04465v1, 2018.
  * [4] M. Przyborowski, T. Tajmajer, Ł. Grad, A. Janusz, P. Biczyk, D. Ślęzak, ”Toward Machine Learning on Granulated Data - a Case of Compact Autoencoder-based Representations of Satellite Images”, IEEE BigData 2018: pp. 2657-2662, 2018.
  * [5] S. Ozkan and G. B. Akar, ”mproved deep spectral convolution network for hyperspectral unmixing with multinomial mixture kernel and endmember uncertainty,” arXiv:1808.01104v1, 2018.
  * [6] G. Peyré and M. Cuturi, ”Computational Optimal Transport”, arXiv:1803.00567v4, 2020.
  * [7] M. Śmieja, M. Wołczyk, J. Tabol, B. C. Geiger, ”SeGMA: Semi-Supervised Gaussian Mixture Auto-Encoder”, arXiv:1906.09333v2, 2019.
  * [8] A. Takatsu and T. Yokota, ”Cone structure of $L2$-Wasserstein spaces”, arXiv:0812.2752v3, 2009.


# Prediction of Progression to Alzheimer’s disease with Deep InfoMax ††thanks:
Corresponding author: Alex Fedorov<EMAIL_ADDRESS>**Data used in the
preparation of this article were obtained from the Alzheimer’s Disease
Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the
investigators within the ADNI contributed to the design and implementation of
ADNI and/or provided data but did not participate in analysis or writing of
this report.

Alex Fedorov12, R Devon Hjelm345, Anees Abrol12, Zening Fu1, Yuhui Du16,
Sergey Plis1, Vince D. Calhoun12
and for the Alzheimer’s Disease Neuroimaging Initiative  1The Mind Research
Network, Albuquerque, New Mexico, USA  2Department of Electrical and Computer
Engineering, University of New Mexico, Albuquerque, New Mexico, USA
3Microsoft Research, Montreal, Canada  4Montreal Institute for Learning
Algorithms, Montreal, Canada  5Department of Computer Science and Operations
Research, University of Montreal, Montreal, Canada  6School of Computer &
Information Technology, Shanxi University, Taiyuan, China

###### Abstract

Arguably, unsupervised learning plays a crucial role in the majority of
algorithms for processing brain imaging. A recently introduced unsupervised
approach Deep InfoMax (DIM) is a promising tool for exploring brain structure
in a flexible non-linear way. In this paper, we investigate the use of
variants of DIM in a setting of progression to Alzheimer’s disease in
comparison with supervised AlexNet and ResNet inspired convolutional neural
networks. As a benchmark, we use a classification task between four groups:
patients with stable, and progressive mild cognitive impairment (MCI), with
Alzheimer’s disease, and healthy controls. Our dataset is comprised of 828
subjects from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database.
Our experiments highlight encouraging evidence of the high potential utility
of DIM in future neuroimaging studies.

###### Index Terms:

CNN, MRI, Deep InfoMax, classification, unsupervised

##  I Introduction

According to [1], the economic costs of mental disorders have the highest
impact on economic growth, direct and indirect costs and the statistical value
of life. One essential tool for better understanding mental illness is to use
noninvasive neuroimaging (e.g., structural magnetic resonance imaging (MRI)
images) along with machine learning to learn brain structure.

Deep Learning has been integral to the successes of machine learning for
numerous demanding real-world applications, e.g., state-of-the-art image
classification [2] and self-driving cars [3]. While many of Deep Learning’s
successes involve supervised learning, supervised approaches can fail when
data annotation (e.g., labels) is limited or unavailable. When there is
sufficient data, supervised models can not only perform well on holdout sets
but provide representations that generalize well to other supervised settings
[4]. However, when there is insufficient data, a supervised learner tends to
discriminate on low-level (e.g., pixel-level, _trivial_) information, which
hurts generalization performance. A model that generalizes well needs to
extract meaningful high-level information (e.g., a collection of important
features at the input level). In order to address this, many successful
applications of machine learning to neuroscience rely on unsupervised learning
[5, 6, 7, 8] to extract representations of brain imaging data. These
representations are then used as input to an off-the-shelf classifier (i.e.,
semi-supervised learning).

However, prior work on unsupervised learning of brain imaging data is either
linear or weakly nonlinear [5, 6] or are highly restrictive in
parameterization [7], and do not represent flexible methodology for learning
representations.

In this work, we explore using DIM [9] to learn deep non-linear
representations of neuroimaging data as an output of a convolutional neural
network. DIM works by maximizing the mutual information between a high-level
feature vector and low-level feature maps of a highly flexible convolutional
_encoder_ network by training a second neural network that maximizes a lower
bound on a divergence (probabilistic measure of difference) between the joint
or the product of marginals of the encoder input and output. The estimates
provided by this second network can be used to maximize the mutual information
of the features in the encoder with the input. Unlike other popular
unsupervised auto-encoding approaches such as VAE [10], DIM doesn’t require a
decoder. Hence it significantly reduces memory requirements of the model for
volumetric data.

We evaluate DIM by performing a downstream classification task between four
groups: patients with stable and progressive MCI, with Alzheimer’s disease and
healthy controls, using only the resulting representation from DIM as input to
the classifier. We compare DIM to two convolutional networks with AlexNet [11]
and ResNet [12] inspired architectures trained with supervised learning. On
strict evaluation, we show comparable performance to supervised methods and to
previously reported [13, 14, 15, 16] classification performance.

##  II Materials and Methods

###  II-A Deep InfoMax

Let $\mathbf{X}:=\\{x^{(i)}\in\mathcal{X}\\}$ and
$\mathbf{Z}:=\\{z^{(i)}\in\mathcal{Z}\\}$ be the input and output variables of
a neural network encoder, $E_{\phi}:\mathcal{X}\rightarrow\mathcal{Z}$ with
parameters $\phi$, where $\mathcal{X}$, and $\mathcal{Z}$ are its domain and
range. We wish to find the parameters that maximize the following objective:

$(\hat{\phi},\hat{\theta})_{G}=\operatorname*{\arg\max}_{\phi,\theta}\hat{\mathcal{I}}_{\theta}(X;Z),$
(1)

where $\hat{\mathcal{I}}_{\theta}$ is the mutual information estimate provided
by a different network with parameters $\theta$, and $Z=E_{\phi}(X)$ is the
output of the encoder.

A parametric estimator for the mutual information can be found by training a
_statistics network_ to maximize a lower bound based on the Fenchel-dual [17]
or the Donsker-Varadhan representation [18, 19] of the Kullback–Leibler
divergence $D_{KL}$. The Donsker-Varadhan-based estimator is a consistent,
asymptotically unbiased estimator has been shown to outperform nonparametric
estimators, and can also be used to improve deep generative models [19].
However, $D_{KL}$ is unbounded, which can be problematic if the above
estimators are used for training deterministic neural network encoders. [9]
showed that using an estimator based on the Jensen-Shannon divergence (JSD)
(i.e., simple binary cross-entropy) is more stable and works well in practice,
and it has been shown that this estimator also yields a good estimator for
mutual information [9, 20]:

$\displaystyle\begin{split}\hat{\mathcal{I}}_{\phi,\theta}^{\text{(JSD)}}(X;E_{\phi}(X)):=&\mathbb{E}_{\mathbb{P}_{X}}[-\text{sp}(-T_{\theta}(X,E_{\phi}(X)))]-\\\
&\mathbb{E}_{\mathbb{P}_{X}\otimes\mathbb{P}_{X}}[\text{sp}(T_{\theta}(X^{\prime},E_{\phi}(X)))],\end{split}$
(2)

where $T_{\theta}$ is a statistics network with parameters $\theta$,
$\text{sp}=\log(1+e^{z})$ (softplus function) and $X^{\prime}$ is another
input sampled from the data distribution independently from $X$. In addition,
the Noise-Contrastive variant of the estimator (NCE) [21] was shown to work
well in practice [9]:

$\displaystyle\begin{split}&\hat{\mathcal{I}}_{\phi,\theta}^{\text{(NCE)}}(X;E_{\phi}(X)):=\\\
&\mathbb{E}_{\mathbb{P}}\left[T_{\theta}(X,E_{\phi}(X))-\log\sum_{X^{\prime}\in\mathcal{X}_{b}}e^{T_{\theta}(X^{\prime},E_{\phi}(X))}\right].\end{split}$
(3)

Here, $\mathcal{X}_{b}=\\{X\\}\bigcup\mathcal{X}_{n}$ are a set of samples
where $\mathcal{X}_{n}$ are a set of _negative samples_ drawn from the data
distribution, such that there is exactly one positive example in
$\mathcal{X}_{b}$ ($X$ occurs exactly once).

[9] showed that maximizing the mutual information between the _complete_ input
and output of an encoder are insufficient for learning good representations
for downstream classification tasks, as this approach can still focus on
lower-level “trivial” localized details. Instead, they show that maximizing
the mutual information between the high-level representation, $Z=E_{\phi}(X)$
and _patches_ of an input image can achieve highly competitive results. The
intuition is that this approach encourages the high-level representation to
learn information that is _shared_ across the input. It is suitable for many
classification tasks, as we expect that class-discriminative features should
be evident across many spatial locations of the input. For a convolutional
encoder $E$, the _local_ DIM objective can be written in a compact form:

$(\hat{\phi},\hat{\theta})_{L}=\operatorname*{\arg\max}_{\phi,\theta}\frac{1}{M^{2}}\sum_{i=1}^{M^{2}}\hat{\mathcal{I}}_{\phi,\theta}(C^{(i)}_{\phi}(X);E_{\phi}(X)),$
(4)

where $C^{(i)}_{\phi}(X)$ is a feature map location from encoder (with a
limited receptive field corresponding to an input patch with size $M$) at some
intermediate layer of the network.

Due to stronger performance of AlexNet architecture (Section (II-B)) in our
experiments (see Section (IV)) we used it as an encoder for DIM method. Last
linear layer of AlexNet we changed with a layer for $64$-dimensional output
representation.

To estimate mutual information using eq. (4) we used the encode-and-dot-
product architecture (Fig. $6$ from [9]). First, patches $C^{(i)}_{\phi}(X)$
taken from third convolutional layer of AlexNet were mapped using
convolutional encoder-and-dot architecture (Tab. $9$ from [9]) with $512$
units and their representation $Z=E_{\phi}(X)$ — linear encoder-and-dot
architecture (Tab. $8$ from [9]). Then flattened encoded mappings of patches
and representations were combined using the dot product to create real and
fake samples efficiently. The real sample is a dot product of a “local” patch
and its “global” representation mappings, while fake — between mapping of some
“local” patch with global representation coming from an unrelated input.
Eventually we estimated JSD based loss eq. (2) and NCE — eq. (3) using these
samples. Since NCE needs to have more negative samples to be competitive with
JSD [9], all possible combinations between the patch and representation
mappings were used a similar way to create negative samples.

To evaluate the performance of the learned representation by DIM, we trained
three additional neural networks using as input features output from last
convolutional layer with size $128\times 2\times 2\times 2$, the first fully
connected layer with $1024$ units, and final fully connected layer with
$64$-dimensional representation, which we call as Conv, FC, and $Z$. The
classifiers are composed of one fully-connected layer with $200$ hidden units,
dropout [22] with $p=0.1$, batch normalization [23] and a ReLU [24]
activation.

TABLE I: AlexNet and ResNet architectures AlexNet
---
3D Conv $(1,64,5,2,0)$ \- BN 3D - ReLU - MP 3D $(3,3)$
3D Conv $(64,128,3,1,0)$ \- BN 3D - ReLU - MP 3D $(3,3)$
3D Conv $(128,192,3,1,1)$ \- BN 3D - ReLU
3D Conv $(192,192,3,1,1)$ \- BN 3D - ReLU
3D Conv $(192,128,3,1,1)$ \- BN 3D - ReLU - MP 3D $(3,3)$
Linear $(1024,1024)$ \- BN 1D - ReLU
Linear $(1024,4)$ \- SoftMax - ArgMax
ResNet
3D Conv $(1,64,3,2,0)$ \- BN 3D - ReLU - MP 3D $(3,3)$
Residual Layer 1
BB0 - 2 x (3D Conv $(64,64,3,1,1)$ \- BN 3D - ReLU)
BB1 - 2 x (3D Conv $(64,64,3,1,1)$ \- BN 3D - ReLU)
Residual Layer 2
BB0 - 3D Conv $(64,128,3,2,1)$ \- BN 3D - ReLU
BB0 - 3D Conv $(128,128,3,1,1)$ \- BN 3D - ReLU
BB0 downsample - 3D Conv $(64,128,3,2,1)$ \- BN 3D
BB1 - 2 x (3D Conv $(128,128,3,2,1)$ \- BN 3D - ReLU)
Residual Layer 3
BB0 - 3D Conv $(128,256,3,2,1)$ \- BN 3D - ReLU
BB0 - 3D Conv $(256,256,3,1,1)$ \- BN 3D - ReLU
BB0 downsample - 3D Conv $(128,256,3,2,1)$ \- BN 3D
BB1 - 2 x (3D Conv $(256,256,3,2,1)$ \- BN 3D - ReLU)
MaxPool 3D $(3,3)$
Linear $(2048,1024)$ \- BN 1D - ReLU
Linear $(1024,4)$ \- SoftMax - ArgMax

###  II-B Supervised baselines

As baselines we have considered supervised methods — two convolutional
networks, one based on a simplified AlexNet [11] architecture and the other a
ResNet [12] architecture. Both networks use convolutions and max pooling with
volumetric kernels, batch normalization, ReLU and two fully connected layers
in the end (see Tab. (I) for details). The notations in Tab. (I) denotes: _BN_
for batch normalization, _BB_ — a basic block, _MP $(k,s)$_ — max pooling with
kernel size $k$ and stride $s$, for convolutions _$(i,o,k,s,p)$_ — a number of
input and output channels, a kernel size, a stride and a padding
respectively). Cross-entropy loss used as a training objective.

###  II-C Regularization

For small datasets, it is common to penalize the number of the model
parameters by driving most of them to zero using $L_{1}$ regularization.
Formally, this penalty is defined as:

$L_{1}(\omega)=\lambda||\omega||_{1}=\lambda\sum_{i}|\omega_{i}|,$ (5)

where $\omega$ is parameter vector of the model and $\lambda$ — coefficient.
$L_{1}$ regularization imposes a sparse solution. This penalty is added to
JSD, NCE and cross-entropy losses in different setting. For our experiments we
used $\lambda=1$.

##  III Experiments

###  III-A Datasets and preprocessing

For the downstream classification task, the data was obtained from the ADNI
database adni.loni.usc.edu (for up-to-date information, see www.adni-
info.org). We use T1w MRI images of $830$ subjects with four different groups:
patients with stable, and progressive MCI, Alzheimer’s disease and healthy
controls.

Structural MRI (sMRI) data was pre-processed to grey matter volume (modulated)
maps using SPM12 toolbox. To segment grey matter, the MRI images were
spatially normalized and smoothed by 6 mm full width at half maximum (FWHM) 3D
Gaussian kernel. After quality control, two subjects from ADNI dataset were
excluded. The final dataset consisted of $828$ subjects with a volume size of
$[121,145,121]$.

###  III-B Experimental setup

####  III-B1 Data

The dataset was divided in approximately $93\%$ and $7\%$ subjects for cross-
validation and hold-out test sets using a stratified split. Then, $93\%$
subjects were split into five stratified folds.

For AlexNet and ResNet architectures, we used simple data augmentation of the
training dataset to reduce overfitting to the small number of annotated
samples available. Our augmentation consisted of zero padding and random
cropping to size $128$ along all dimensions along with randomly flipping the
input with probability $0.5$ for each axis. The whole brain was included in
the crop.

For DIM, we didn’t use data augmentation, but we used zero padding to make
sure that input size is equal to $128$ along all dimensions.

####  III-B2 Training

The models were trained using the AMSGrad [25] optimizer with learning rate
$0.001$ for CNN models and $0.0001$ for DIM using a batch size of $8$ but
dropping the incomplete last batch. The training of the supervised
architectures was performed for $500$ epochs, DIM — for $1000$ epochs as pre-
training and for $1000$ epochs for training the classifiers on top of frozen
features from the encoder.

####  III-B3 Evaluation

Since the dataset is not completely balanced, the evaluation was performed
using balanced accuracy [26], defined as the average of recall of each class
(implementation in scikit-learn [27]).

####  III-B4 Implementation and hardware

The implementation was written using Deep Learning frameworks PyTorch [28] and
Cortex [29]. The DIM code is based on openly available DIM implementation
[30]. The experiments were performed on NVIDIA GeForce Titan X Pascal and 1080
Ti and 8 CPU threads.

Figure 1: Performance of the models TABLE II: Performance

Model Balanced Accuracy Stratified 5-Fold Balanced Accuracy Hold-out Mean gap
Wilcoxon test Stat $p$-value AlexNet $47.31\pm 2.69$ $50.36\pm 4.62$ $3.05$
6.5 $0.290$ AlexNet Aug $49.82\pm 3.18$ $52.14\pm 7.41$ $2.32$ 7.0 $0.554$
Sparse AlexNet Aug $51.85\pm 5.14$ $51.07\pm 3.91$ $0.78$ N/A N/A ResNet
$47.9\pm 4.28$ $43.57\pm 8.71$ $4.33$ 10.0 0.034 ResNet Aug $50.72\pm 3.8$
$47.14\pm 6.51$ $3.58$ 14.0 0.039 Sparse ResNet Aug $50.07\pm 3.09$ $43.93\pm
6.26$ $6.14$ 10.0 0.033 JSD Conv $47.82\pm 2.11$ $45.36\pm 7.74$ $2.46$ 13.5
$0.052$ JSD Conv SS $48.83\pm 2.79$ $47.5\pm 3.24$ $1.33$ 13.5 $0.052$ Sparse
JSD Conv $49.61\pm 2.35$ $44.29\pm 5.27$ $5.32$ 15.0 0.022 Sparse JSD Sparse
Conv $49.29\pm 4.22$ $48.57\pm 5.27$ $0.72$ 6.0 $0.054$ NCE Conv $45.82\pm
2.82$ $40.71\pm 5.98$ $5.11$ 15.0 0.022 NCE Conv SS $45.08\pm 1.94$ $43.57\pm
2.99$ $1.51$ 10.0 0.034 Sparse NCE Conv $47.59\pm 3.21$ $45.0\pm 3.19$ $2.59$
15.0 0.022 Sparse NCE Sparse Conv $47.23\pm 2.69$ $39.64\pm 2.93$ $7.59$ 15.0
0.022 JSD FC $44.18\pm 2.05$ $30.71\pm 6.11$ $13.47$ 15.0 0.022 JSD FC SS
$49.77\pm 2.73$ $45.36\pm 5.73$ $4.41$ 14.0 0.040 Sparse JSD FC $45.42\pm
4.03$ $37.5\pm 5.05$ $7.92$ 15.0 0.022 Sparse JSD Sparse FC $46.6\pm 3.92$
$45.0\pm 3.87$ $1.6$ 15.0 0.022 NCE FC $44.89\pm 2.93$ $32.86\pm 4.48$ $12.03$
15.0 0.022 NCE FC SS $46.53\pm 2.59$ $40.36\pm 4.66$ $6.17$ 15.0 0.020 Sparse
NCE FC $46.2\pm 2.64$ $37.86\pm 6.49$ $8.34$ 15.0 0.022 Sparse NCE Sparse FC
$47.01\pm 3.12$ $36.43\pm 7.84$ $10.58$ 15.0 0.022 JSD Z $43.12\pm 1.1$
$31.43\pm 5.14$ $11.69$ 15.0 0.021 JSD Z SS $48.44\pm 2.79$ $44.64\pm 4.19$
$3.8$ 15.0 0.022 Sparse JSD Sparse Z $47.74\pm 5.12$ $40.36\pm 6.39$ $7.38$
15.0 0.022 Sparse JSD Z $45.27\pm 3.38$ $40.0\pm 6.98$ $5.27$ 15.0 0.022 NCE Z
$44.14\pm 3.84$ $36.07\pm 2.93$ $8.07$ 15.0 0.022 NCE Z SS $45.28\pm 3.55$
$40.0\pm 7.84$ $5.28$ 14.0 0.040 Sparse NCE Sparse Z $46.63\pm 2.69$ $37.14\pm
4.45$ $9.49$ 15.0 0.022 Sparse NCE Z $45.27\pm 4.85$ $43.21\pm 5.56$ $2.06$
15.0 0.021

##  IV Results

The final trained models used further to evaluate the performance were
selected based on the best-balanced accuracy but from a checkpoint where the
validation score was lower than the training score. We gave the model a burn-
in period before applying this rule to deal with initial stochasticity. The
models notations are as follows: _Aug_ denotes augmentation of the training
dataset, the first _sparse_ — a model trained with $L_{1}$ regularization, the
second — a classifier on top of the frozen features from encoders trained
using $L_{1}$ regularization, _SS_ — stands for training an unsupervised model
with an additional supervised loss from $Z$-classifier.

Table II reports the balanced accuracy rates including mean, standard
deviation values, and the gap between mean values on cross-validation and
hold-out. The bold text distinguishes the best scores and the name of the
models. The last column shows $p$-value and statistic for the one-sided
Wilcoxon test. The bold $p$-values indicate acceptance of the null hypothesis.
The test was performed to compare each method with the best model (_Sparse
AlexNet Aug_) based on the five values of balanced accuracy on hold-out. An
alternative hypothesis is that the model Sparse AlexNet Aug is better. Fig. 1
highlights the distributions of the performance.

With all modifications, _ResNet_ shows a lower performance on hold-out (at
most $47.14\pm 6.51$) than _AlexNet_. It is reasonable since the capacity of
the ResNet architecture is larger and the dataset is small. For $\alpha=0.05$
Wilcoxon test also rejects H0 supporting the worse performance of ResNet.
Performance of _JSD Conv_ , _JSD Conv SS_ , _Sparse JSD Sparse Conv_ ,
_AlexNet_ , _AlexNet Aug_ is statistically indistinguishable from that of
_Sparse AlexNet Aug_. Follows that unsupervised DIM has comparable performance
to supervised methods.

Among DIM variants, JSD has higher scores than NCE. Lower scores of NCE can be
explained by its requirement of a large number of negative samples during
training to be competitive with JSD. Our dataset is not large enough to
support the needed level of negative sampling.

The best score with convolutional features—$48.57\pm 5.27$—was obtained by an
encoder and classifier trained with $L_{1}$ regularization which is the
_Sparse JSD Sparse Conv_ model. For features from the fully-connected layer —
_JSD FC SS_ model with $45.36\pm 5.73$ using semi-supervised loss was the
best. However, _Sparse JSD Sparse FC_ has similar results $45.0\pm 3.87$ and a
smaller mean gap $1.6$ but it has a lower mean cross-validation score by
$3.17\%$ . For the smallest $64$-dimensional representation, semi-supervised
model _JSD Z SS_ gives the best performance $44.64\pm 4.19$, but similar
result $43.21\pm 5.56$ were obtained by _Sparse NCE Z_ model. Semi-supervised
loss and $L_{1}$ regularization improved models’ generalization by reducing
the gap between cross-validation and hold-out scores. The observed degradation
in performance between _Conv_ , _FC_ , and $Z$ can be explained by the reduced
capacity of the features. $L_{1}$ regularization and dropout could also be
adjusted. However, a more compact input representation can be of independent
use, for example, for dimensionality reduction.

In previous studies, the best reported accuracy for the ResNet architecture in
a 4-class sMRI classification task was $54\%$ [13], while stacked autoencoders
(SAE) [15] reached for sMRI only $46.30\pm 4.24$ and for sMRI+PET $53.79\pm
4.76$, and DW-S2MTL [16] — for sMRI $47.83$ or for sMRI+PET+CSF $53.72$. Our
values can’t be completely comparable since the evaluation is different.
Reproduced ResNet can be used as a proxy to estimate performance relative to
this prior work. Note, however, it is not one of the best-performing methods
in our study.

##  V Conclusions

This work proposes an unsupervised method DIM for learning representations
from structural neuroimaging data. The evaluation of the prediction of
progression to Alzheimer’s disease demonstrates results comparable to
supervised methods. In the future, we will scale up our experiments with
increased sample size and address the cases of other diseases. Our future
efforts will also be focused on the multi-modal fusion of brain imaging data
[31] to increase the predictive strength of the model.

## Acknowledgement

This study is supported by NIH grants R01EB020407, R01EB006841, P20GM103472,
P30GM122734.

Data collection and sharing for this project was funded by the Alzheimer’s
Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant
U01 AG024904) and DOD ADNI (Department of Defense award number
W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the
National Institute of Biomedical Imaging and Bioengineering, and through
generous contributions from the following: AbbVie, Alzheimer’s Association;
Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.;
Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.;
Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La
Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE
Healthcare; IXICO Ltd.;Janssen Alzheimer Immunotherapy Research & Development,
LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity;
Lundbeck; Merck & Co., Inc.;Meso Scale Diagnostics, LLC.; NeuroRx Research;
Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.;
Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition
Therapeutics. The Canadian Institutes of Health Research is providing funds to
support ADNI clinical sites in Canada. Private sector contributions are
facilitated by the Foundation for the National Institutes of Health
(www.fnih.org). The grantee organization is the Northern California Institute
for Research and Education, and the study is coordinated by the Alzheimer’s
Therapeutic Research Institute at the University of Southern California. ADNI
data are disseminated by the Laboratory for Neuro Imaging at the University of
Southern California.

## References

  * [1] S. Trautmann, J. Rehm, and H.-U. Wittchen, “The economic costs of mental disorders: Do our societies react appropriately to the burden of mental disorders?” _EMBO reports_ , p. e201642951, 2016.
  * [2] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2017, pp. 2261–2269.
  * [3] Y. Xiong, R. Liao, H. Zhao, R. Hu, M. Bai, E. Yumer, and R. Urtasun, “Upsnet: A unified panoptic segmentation network,” _arXiv preprint arXiv:1901.03784_ , 2019.
  * [4] L. Zhang, T. Xiang, and S. Gong, “Learning a deep embedding model for zero-shot learning,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2017, pp. 3010–3019.
  * [5] V. D. Calhoun, T. Adali, G. D. Pearlson, and J. Pekar, “A method for making group inferences from functional mri data using independent component analysis,” _Human brain mapping_ , vol. 14, no. 3, pp. 140–151, 2001.
  * [6] R. D. Hjelm, V. D. Calhoun, R. Salakhutdinov, E. A. Allen, T. Adali, and S. M. Plis, “Restricted boltzmann machines for neuroimaging: an application in identifying intrinsic networks,” _NeuroImage_ , vol. 96, pp. 245–260, 2014\.
  * [7] E. Castro, R. D. Hjelm, S. M. Plis, L. Dinh, J. A. Turner, and V. D. Calhoun, “Deep independence network analysis of structural brain imaging: application to schizophrenia,” _IEEE transactions on medical imaging_ , vol. 35, no. 7, pp. 1729–1740, 2016.
  * [8] S. M. Plis, D. R. Hjelm, R. Salakhutdinov, E. A. Allen, H. J. Bockholt, J. D. Long, H. J. Johnson, J. S. Paulsen, J. A. Turner, and V. D. Calhoun, “Deep learning for neuroimaging: a validation study,” _Frontiers in neuroscience_ , vol. 8, p. 229, 2014.
  * [9] R. D. Hjelm, A. Fedorov, S. Lavoie-Marchildon, K. Grewal, A. Trischler, and Y. Bengio, “Learning deep representations by mutual information estimation and maximization,” _arXiv preprint arXiv:1808.06670_ , 2018.
  * [10] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [11] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in Neural Information Processing Systems 25_ , F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, Eds. Curran Associates, Inc., 2012, pp. 1097–1105.
  * [12] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [13] A. Abrol, M. Bhattarai, A. Fedorov, Y. Du, S. Plis, and V. Calhoun, “Deep residual learning for neuroimaging: An application to predict progression to alzheimer’s disease,” _bioRxiv_ , p. 470252, 2018.
  * [14] S. Vieira, W. H. Pinaya, and A. Mechelli, “Using deep learning to investigate the neuroimaging correlates of psychiatric and neurological disorders: Methods and applications,” _Neuroscience & Biobehavioral Reviews_, vol. 74, pp. 58–75, 2017.
  * [15] S. Liu, S. Liu, W. Cai, H. Che, S. Pujol, R. Kikinis, D. Feng, M. J. Fulham, and ADNI, “Multimodal neuroimaging feature learning for multiclass diagnosis of alzheimer’s disease,” _IEEE Transactions on Biomedical Engineering_ , vol. 62, no. 4, pp. 1132–1140, April 2015.
  * [16] H.-I. Suk, S.-W. Lee, D. Shen, A. D. N. Initiative _et al._ , “Deep sparse multi-task learning for feature selection in alzheimer’s disease diagnosis,” _Brain Structure and Function_ , vol. 221, no. 5, pp. 2569–2587, 2016.
  * [17] S. Nowozin, B. Cseke, and R. Tomioka, “f-gan: Training generative neural samplers using variational divergence minimization,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 271–279.
  * [18] M. D. Donsker and S. S. Varadhan, “Asymptotic evaluation of certain markov process expectations for large time. iv,” _Communications on Pure and Applied Mathematics_ , vol. 36, no. 2, pp. 183–212, 1983.
  * [19] M. I. Belghazi, A. Baratin, S. Rajeswar, S. Ozair, Y. Bengio, A. Courville, and R. D. Hjelm, “Mine: mutual information neural estimation,” _arXiv preprint arXiv:1801.04062_ , 2018.
  * [20] B. Poole, S. Ozair, A. van den Oord, A. A. Alemi, and G. Tucker, “On variational lower bounds of mutual information.”
  * [21] A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _arXiv preprint arXiv:1807.03748_ , 2018\.
  * [22] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” _The Journal of Machine Learning Research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [23] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” _arXiv preprint arXiv:1502.03167_ , 2015.
  * [24] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in _Proceedings of the 27th international conference on machine learning (ICML-10)_ , 2010, pp. 807–814.
  * [25] S. J. Reddi, S. Kale, and S. Kumar, “On the convergence of adam and beyond,” _arXiv preprint arXiv:1904.09237_ , 2019.
  * [26] K. H. Brodersen, C. S. Ong, K. E. Stephan, and J. M. Buhmann, “The balanced accuracy and its posterior distribution,” in _Pattern recognition (ICPR), 2010 20th international conference on_. IEEE, 2010, pp. 3121–3124.
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [28] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “Automatic differentiation in pytorch,” 2017.
  * [29] rdevon/cortex: A machine learning library for pytorch. [Online]. Available: https://github.com/rdevon/cortex
  * [30] rdevon/dim: Deep infomax (dim), or ”learning deep representations by mutual information estimation and maximization”. [Online]. Available: https://github.com/rdevon/DIM
  * [31] V. D. Calhoun and J. Sui, “Multimodal fusion of brain imaging data: a key to finding the missing link (s) in complex mental illness,” _Biological psychiatry: cognitive neuroscience and neuroimaging_ , vol. 1, no. 3, pp. 230–244, 2016.


#  Time evolution of the von Neumann entropy in open quantum systems

Kohei Kobayashi1 1Global Research Center for Quantum Information Science,
National Institute of Informatics, 2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo
101-8340, Japan

###### Abstract

Control of open quantum dynamics is of great interest for realizing quantum
technologies. Therefore, it is an important task to quantify and characterize
the entropy for open quantum systems under decoherence. In this paper, we
study the time evolution of the von Neumann entropy for open quantum systems
described by the Lindblad master equation. Note that, in particular, when the
decoherence corresponds to the measurement for the observable in the system,
the von Neumann entropy tends to monotonically increases as the variance
becomes larger. Furthermore, we present a lower bound of the von Neumann
entropy in the long-time limit. This lower bound has advantages of being
straightforwardly calculated and applicable to a general Markovian open
quantum system.

von Neumann entropy, open quantum system, Lindblad master equation

##  I Introduction

Control of quantum dynamics is of significant interest for realizing quantum
technologies. Because real quantum systems are subjected to the interaction
with uncontrollable environments, an open quantum system offers a powerful
tool for modeling actual phenomena in quantum mechanics Nielsen . Due to the
interaction between system and environments, the loss of energy and phase
informationm, so-called decoherence, is inevitably caused. Therefore, it is
important to characterize how quantum information vary in an open quantum
system.

The calculation of entropy in open quantum systems is generally an essential
task. Obtaining a limit on the rate of entropy change or quantifyig the
entropy itself enable us to understand the dynamical behavior of
nonequilibrium quantum states, as a result, which leads to ideal control
toward desired states.

Towards this problem, the von Neumann entropy in open quantum systems has
already been discussed almost 30 years ago Benatti . In recent years, several
results for studying the entropy in open quantum dynamics have been reported
Audenaert ; Ou ; Abe ; Bakh ; Mohan ; for instance, the time derivative of the
von Neumann entropy in a time-dependent harmonic oscillator was studied Ou ,
the lower bound on the Re’nyi entropy rate was given Abe , and the relation of
the von Neumann entropy and the quantum speed limit was investigated Mohan .
However, these results remain as a formal mathematical expression and physical
intuition is hard to extract, and also, there has been no result giving any
estimate on the entropy itself.

Motivated by the above background, we study the time evolution of the von
Neumann entropy for open quantum systems descrived by the Lindblad master
equation. From an obtained inequality for the time derivative of the von
Neumann entropy, we derive a condition that the entropy in the quantum system
monotonically increases under decoherence. Note that, in particular, when the
decoherence is the measurement for the observable of the system, this
condition is characterized by the variance of it. Moreover, we present a lower
bound of the von Neumann entropy of an open quantum system in the long-time
limit. This bound gives a limit on how low the entropy in the quantum system
can maintain in the presence of decoherence. In addition, as a mathematical
tool, this bound has advantages of being easy to calculate and applicable to a
general Markovian open quantum system.

##  II Entropy bound

###  II.1 Controlled quantum dynamics

Let $\mathcal{H}$ denote the Hilbert space, where the quantum state is
defined. The state of a quantum system is represented by the density operator
$\rho$, which satisfies following properties: $\rho=\rho^{\dagger}$, $\rho\geq
0$, and ${\rm Tr}(\rho)=1$. The identity operator is denoted by $I$. Let us
assume that the qauntum state $\rho$ is Markov (i.e., the system is weakly
coupled to the environment and the dynamics of the environment changes much
faster than the system). Therefore, the information only flows from the system
to the environment and there is no information flowing back. In this case, the
time evolution of the Markovian open quantum system $\rho_{t}$ is described by
the following equation Lindblad :

$\frac{d\rho_{t}}{dt}=-i[H,\rho_{t}]+\mathcal{D}[L]\rho_{t},$ (1)

which is referred as the Lindblad master equation. $H$ is the system
Hamiltonian, and $L$ is the Lindblad operator representing the decoherence
process, and thus $\mathcal{D}[A]\rho=A\rho
A^{\dagger}-A^{\dagger}A\rho/2-\rho A^{\dagger}A/2$ (we set $\hbar=1$ in the
following). The Lindblad master equation is known as the most general
generator of the quantum dynamical semigroup Heinz .

###  II.2 Main result

We consider the von Neumann entropy $S_{t}$:

$S_{t}=-{\rm Tr}\left(\rho_{t}\ln\rho_{t}\right),\ \ \ 0\leq S_{t}\leq\ln d,$
(2)

where $d$ is the rank of the quantum system. The von Neumann entropy
quantifies the amount of information contained in a state when many identical
and independent copies of the state are available. $S_{t}$ takes zero when
$\rho_{t}$ is pure, and the maximum is achieved only when the system is
maximally mixed $\rho_{t}=I/d$. Because $\rho_{t}$ can be diagonalized by a
certain orthonormal basis $\\{|u_{j}\rangle\\}_{j}$ as
$\rho_{t}=\sum_{j}\lambda_{j}|u_{j}\rangle\langle u_{j}|$, the von Neumann
entropy is also written as

$S_{t}=-\sum_{j}\lambda_{j}\ln\lambda_{j},$ (3)

which is analogous to the Shannon entropy in classical information theory
Shannon .

The time evolution of the von Neumann entropy is calculated as follows:

$\displaystyle\frac{dS_{t}}{dt}$ $\displaystyle=$ $\displaystyle-{\rm
Tr}\left(\frac{d\rho_{t}}{dt}\ln\rho_{t}\right)$ (4) $\displaystyle=$
$\displaystyle-{\rm
Tr}\left\\{(-i[H,\rho_{t}]+\mathcal{D}[L]\rho_{t})\ln\rho_{t}\right\\}$
$\displaystyle=$ $\displaystyle-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\ln\rho_{t}\right)+\frac{1}{2}{\rm
Tr}\left(L^{\dagger}L\rho_{t}\ln\rho_{t}\right)+\frac{1}{2}{\rm
Tr}\left(\rho_{t}L^{\dagger}L\ln\rho_{t}\right)$ $\displaystyle=$
$\displaystyle{\rm Tr}\left(L^{\dagger}L\rho_{t}\ln\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\ln\rho_{t}\right),$

where we used ${\rm Tr}(d\rho_{t}/dt)=0$, $[\rho_{t},\ln\rho_{t}]=0$ and ${\rm
Tr}\left([H,\rho_{t}]\ln\rho_{t}\right)=0$.

Here, due to the fact that $-\rho_{t}\ln\rho_{t}$ is positive semidefinite, by
using the inequality ${\rm Tr}(AB)\leq{\rm Tr}(A){\rm Tr}(B)$ for $A$, $B\geq
0$, the first term of Eq. (4) can be bounded as follows:

$\displaystyle{\rm Tr}\left(L^{\dagger}L\rho_{t}\ln\rho_{t}\right)$
$\displaystyle=$ $\displaystyle-{\rm
Tr}\left[L^{\dagger}L(-\rho_{t}\ln\rho_{t})\right]$ (5) $\displaystyle\geq$
$\displaystyle-{\rm Tr}(L^{\dagger}L)S_{t}=\|L\|^{2}_{\rm F}S_{t},$

where $\|A\|^{2}_{\rm F}:={\rm Tr}(A^{\dagger}A)^{1/2}$ is the Frobenius norm.
Next, to calculate the lower bound of Eq. (4), we introduce the following
inequality:

$-\ln\rho\geq I-\rho.$ (6)

We give the proof of this inequality in Appendix A. Using it, the second term
of Eq. (4) can be bounded as follows:

$\displaystyle{\rm Tr}\left[L\rho_{t}L^{\dagger}(-\ln\rho_{t})\right]$
$\displaystyle\geq$ $\displaystyle{\rm
Tr}\left[L\rho_{t}L^{\dagger}(I-\rho_{t})\right]$ (7) $\displaystyle\geq$
$\displaystyle{\rm Tr}\left(L\rho_{t}L^{\dagger}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right)\geq 0.$

Therefore, we have the inequality:

$\frac{dS_{t}}{dt}\geq-\|L\|^{2}_{\rm F}S_{t}+{\rm
Tr}\left(L^{\dagger}L\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right).$ (8)

Therefore, it is guaranteed that the time evolution of $S_{t}$ is always
positive when

$S_{t}\leq\frac{{\rm Tr}\left(L^{\dagger}L\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right)}{\|L\|^{2}_{\rm F}}.$ (9)

This upper bound is less than one, because

$\frac{{\rm Tr}\left(L^{\dagger}L\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right)}{\|L\|^{2}_{\rm F}}\leq\frac{{\rm
Tr}\left(L^{\dagger}L\right){\rm Tr}\left(\rho_{t}\right)}{\|L\|^{2}_{\rm
F}}\leq 1.$ (10)

Because $S_{t}\leq\ln d$, the inequality (9) gives a meaningful result for
$d\geq 3$.

Furthermore, let us assume that $L$ is self-adjoint (i.e., $L=L^{\dagger}$);
$L$ represents the measurement for an observable in the system. In this case,
the numerator of Eq. (9) is bounded from below as follows:

$\displaystyle{\rm Tr}\left(L^{\dagger}L\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right)$ $\displaystyle=$
$\displaystyle{\rm Tr}\left(L^{2}\rho_{t}\right)-{\rm
Tr}\left[(\rho^{1/2}_{t}L\rho^{1/2}_{t})^{2}\right]$ (11) $\displaystyle\geq$
$\displaystyle{\rm Tr}\left(L^{2}\rho_{t}\right)-{\rm Tr}(L\rho_{t})^{2},$

where ${\rm Tr}(A^{2})\leq{\rm Tr}(A)^{2}$ for $A\geq 0$ was used. It is
identical to the variance of the observable $L$ in $\rho_{t}$. Therefore,
denoting ${\rm Var}[L]:={\rm Tr}(L^{2}\rho_{t})-{\rm Tr}(L\rho_{t})^{2}$, Eq.
(9) is rewritten as

$S_{t}\leq\frac{{\rm Var}[L]}{\|L\|^{2}_{\rm F}}.$ (12)

This result is reasonable because the more the measurement outcomes
fluctuates, the more easier the conditions for the von Neumann entropy
monotonically increasing. On the other hand, by suppressing the observable
variance small, it may be possible to prevent the entropy from monotonically
increasing.

Next, we show that a fundamental lower bound of $S_{t}$ can be derived in the
limit $t\to\infty$ from the inequality (8). To obtain it, let us consider the
function:

$f(x)=-\|L\|^{2}_{\rm F}x+{\rm Tr}\left(L^{\dagger}L\rho_{t}\right)-{\rm
Tr}\left(L\rho_{t}L^{\dagger}\rho_{t}\right).$ (13)

$f(0)\geq 0$, and when $x=\ln d$ ( $\rho_{t}=I/d$), then

$f(\ln d)=\|L\|^{2}_{\rm F}\left(\frac{1}{d}-\frac{1}{d^{2}}-\ln d\right)<0.$
(14)

Then, from the above properties of $f(x)$, $f(S_{*})=0$ has a unique solution
$S_{*}$ in $[0,\ln d]$. If $S_{T}<S_{*}$ at a given evolution time $T$, $S$ is
monotonically increasing. On the other hand, for the range such that
$S_{T}\geq S_{*}$, the inequality (7) does not say anything about the time
evolution of $S_{T}$ for $t\geq T$. As a result, in the long-time limit
$t\to\infty$, we obtain the lower bound:

$S_{\infty}\geq S_{*}:=\frac{\beta}{\|L\|^{2}_{\rm F}},$ (15)

where $\beta:={\rm Tr}\left(L^{\dagger}L\rho_{\infty}\right)-{\rm
Tr}\left(L\rho_{\infty}L^{\dagger}\rho_{\infty}\right)$. That is, $S_{*}$
gives a limit on how low the entropy in quantum system can maintain under
decoherence.

Here we would like to emphasize the performance of $S_{*}$ as a tool for
quantifying the entropy. This lower bound is applicable to a general Markovian
open quantum system without imposing any assumptions on the decoherence $L$.
Also it can be directly calculated once $\rho_{\infty}$ and $L$ are specified,
thanks to which, there is no need for solving any equation or doing integral
calculation. When the steady state is maximally mixed $\rho_{\infty}=I/d$,
$S_{*}$ is given as the function of $d$ as

$\displaystyle S_{*}\left(I/d\right)=\frac{d-1}{d^{2}}.$ (16)

Then, the lower bound takes the maximum when $d=2$, $S_{*}=1/4$ and
$S_{\infty}=\ln 2\approx 0.69$, meanwhile as $d$ increases, its tightness
becomes weak.

Finally, the bound (15) can be generalized to the case where the system is
sujected to multiple decoherence channels $\sum_{j}\mathcal{D}[L_{j}]\rho_{t}$
as follows:

$\displaystyle S_{*}=\frac{\sum_{j}\beta_{j}}{\sum_{j}\|L_{j}\|^{2}_{\rm F}}.$
(17)

##  III Conclusion

In this paper, we have derived an inequality for the time derivative of the
von Neumann entropy in an open quantum system under decoherence. Based on it,
we have characterized the condition that the von Neumann entropy is
monotonically increasing. In particular, it is noted that when the decoherence
corresponds to the measurement for an observable of the system, this condition
is upper bounded by its variance. Moreover, we derived a lower bound on the
von Neumann entropy itself in the long-time limit. This lower bound is
applicable to a general Markovian open system described by the Lindblad master
equation and can be straightforwardly computed. The results obtained in this
paper are expected to be useful for studying the dynamics of quantum
properties in the presence of decoherence. A remaining work is to further
generalize the present discussion to an important dynamics, such as non-
Markovian dynamics or continuously monitored dynamics Handel ; Geremia .

This work was supported by MEXT Quantum Leap Flagship Program Grant
JPMXS0120351339.

##  Appendix A Proof of Eq. (6)

To show the inequality Eq. (6), we use the diagonalization of the density
matrix

$\rho=U\Lambda U^{\dagger},$ (18)

where $U$ is an appropri ate unitary matrix and $\Lambda={\rm
diag}\\{\lambda_{1},\cdots,\lambda_{d}\\}$ with the eigenvalues
$1\geq\lambda_{1}\geq\cdots\geq\lambda_{d}\geq 0$. Then,

$\displaystyle-\ln\rho$ $\displaystyle=$ $\displaystyle-U\ln\lambda
U^{\dagger}$ (19) $\displaystyle=$ $\displaystyle U{\rm
diag}\\{-\ln\lambda_{1},\cdots,-\ln\lambda_{d}\\}U^{\dagger}$
$\displaystyle\geq$ $\displaystyle U{\rm
diag}\\{1-\lambda_{1},\cdots,1-\lambda_{d}\\}U^{\dagger}$ $\displaystyle=$
$\displaystyle U\left(I-\Lambda\right)U^{\dagger}=I-\rho,$

where we used $-\ln x\geq 1-x$ for $x>0$.

## References

  * (1) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information. Cambridge University Press, Cambridge (2010).
  * (2) F. Benatti and H. Narnhofer, Lett. Math. Phys. 15, 325 (1988).
  * (3) K. M. R. Audenaert, J. Math. Phys. A: Math. Theor. 40 (28), 8127-8136 (2007).
  * (4) C. Ou, R.V. Chamberlin, and S. Abe, Physica A 466, 450 (2017).
  * (5) S. Abe, Phys. Rev. E 94, 022106 (2016).
  * (6) F. Bakhshinezhad and A. Rezakhani, arXiv:1808.08305 [quant-ph] (2018).
  * (7) B. Mohan, S. Das, and A. K. Pati, New J. Phys. 24 065003 (2022).
  * (8) G. Lindblad, Commun. Math. Phys. 48 119 (1976).
  * (9) H.-P. Breuer and Francesco Petruccione, The Theory of Open Quantum Systems. Oxford University Press (2002).
  * (10) C. Shannon, Bell System Technical Journal 27, 379-423, 623-656 (1948).
  * (11) R. van Handel, JK. Stockton, and H. Mabuchi, Phys. Rev. A 70, 022106 (2004).
  * (12) J.M. Geremia, Phys. Rev. Lett. 97, 073601 (2006).