# Multi-View Non-negative Matrix Factorization Discriminant Learning via Cross
Entropy Loss

Jian-wei Liu, Yuan-fang Wang, Run-kun Lu, Xionglin Luo  Department of
Automation, College of Information Science and Engineering, China University
of Petroleum , Beijing, Beijing, China

###### Abstract.

Multi-view learning accomplishes the task objectives of classification by
leverag-ing the relationships between different views of the same object. Most
existing methods usually focus on consistency and complementarity between
multiple views. But not all of this information is useful for classification
tasks. Instead, it is the specific discriminating information that plays an
important role. Zhong Zhang et al. explore the discriminative and non-
discriminative information exist-ing in common and view-specific parts among
different views via joint non-negative matrix factorization. In this paper, we
improve this algorithm on this ba-sis by using the cross entropy loss function
to constrain the objective function better. At last, we implement better
classification effect than original on the same data sets and show its
superiority over many state-of-the-art algorithms.

Multi-View learning, matrix factorization, cross entropy loss, discriminative
learning, classification

††ccs: Computing methodologies Dimensionality reduction and manifold
learning††ccs: Computing methodologies Non-negative matrix factorization

##  1\. Introduction

In many scientific data analysis tasks, data are often collected through
different measuring methods, such as various feature extractors or different
sensors, as usually the single particular measuring method cannot
comprehensively extract all the information of the object or entities. In this
case, the features of each example can be naturally divided into groups, each
of which can be considered as a view. For instance, for images and videos,
color features and texture features can be regarded as two views. It is
important to make good use of the information from different views. Multi-view
learning is a branch of machine learning which studies and utilizes the
information and relationship between views.

There are also many entities represented with different views in the real
world, such as web pages (Blum and Mitchell, 1998; Zhou and Burges, 2007a),
multi-lingual news (Cai et al., 2010; Guan et al., 2015; Kumar et al., 2011),
and neuroimaging (Shao et al., 2018, 2012a; Shao et al., 2012b). And most of
relevant literatures of multi-view learning regard the consistency and
complementarity as two main underlying properties of the multi-view data (Xu
et al., 2013), which like the bridges to link all views together (Chaudhuri et
al., 2009; Liu et al., 2013; Ye et al., 2015). Consistency assumptions
indicate that all views share consistent information. Obviously, using only
consistent information to take advantage of multi-view data is not enough,
because each view also contains additional knowledge that other views do not
have (Blum and Mitchell, 1998; Gupta et al., 2010a; Liu et al., 2014a).
Therefore, the complementarity of views is another important properties of
learning multi-view data, which is worth developing and utilizing.

However, it is a doubt that whether the consistent and complementary
information really always support a better classification performance. Because
some experiences of previous experiments show that the predictive performance
of multi-view data may even be worse than only using single-view data in some
real data sets. Zhang Zhong et al. hold point of view that the multi-view data
contains parts that are helpful and unhelpful for classification, which are
called discriminate and non-discriminate information and the consistent or
complementary information does not learn discriminative information directly
(Zhang et al., 2018). The classifier constructed by multi-view data may give
an even worse classification performance if the learned consistent and
complementary information contains no clear discriminative information.

In this paper, we propose multi-view non-negative matrix factorization
discriminant learning via cross entropy loss. we distinguish discriminative
and non-discriminative information existing in the consistent and
complementary parts, and only use discriminative information for
classification. More specifically, multi-view data is factorized into common
part shared across views and view-specific parts existing within each view in
the usual multi-view learning manner. Besides, for both common part and view-
specific part, they are further factorized into discriminative part and non-
discriminative part. In this situation, each view of data is factorized into
four parts: the common discriminative part, the common non-discriminative
part, the specific discriminative part and the specific non-discriminative
part.

In previous research for multi-view learning, mean square error is habitually
used as a loss function for predicting labels, but the truth is not the case.
As the training results of gradient descent show, using mean squared error
term to predict the output labels has a certain degree of defects. Its partial
derivative value of mean squared error term is very small when the probability
of output labels is close to 0 or 1, which could cause the value of partial
derivative to almost vanish at the beginning of training. In order to speed up
the training processes With the decrease of error between predicted label and
actual label values, we choose the cross-entropy loss function instead of mean
squared error. Cross-entropy loss function is usually used to measure the
performance of a classification model whose input is a probability value
between 0 and 1. And cross-entropy loss increases as the predicted probability
diverging from the actual label. Moreover, a supervised constraint is added
objective function to guide the joint NMF factorization to obtain better
discriminative parts.

To find the optimal decomposition, we follow the block coordinate descent
framework (Kim et al., 2014) to solve our proposed objective function. And
only the derived discriminative parts from common part and view-specific parts
are used to construct a classifier. Finally, experimental results on seven
real-world data sets verify the effectiveness of our proposed method.

In conclusion, our contributions are summarized as follows:

(1)We propose a new multi-view non-negative matrix factorization discriminant
learning algorithm, we dub Multi-View non-negative matrix factorization
Discriminant Learning via Cross Entropy Loss(MV-DLCSL), which utilizes the
discriminative information of data to classify, and outperforms many state-of-
the-art algorithms on seven real-world data sets.

(2)To our best knowledge, we introduce cross-entropy loss into the non-
negative matrix factorization framework for the first time, which is more
reasonable objective function to measure predicting output label errors.

(3)To find the optimal decomposition, we follow the block coordinate descent
framework to solve our proposed optimization problem. However, the cross-
entropy loss function contains the nonlinear softmax function, which
normalized the predicting output label probability. For better preforming the
block coordinate descent, we devise approach to obtain the derivative of the
cross-entropy loss function for each matrix vectors.

(4) In order to visually compare the difference between cross entropy loss and
mean squared error, we use a visualization tool, t-SNT, which is very suitable
for high-dimensional data dimensionality reduction to 2D or 3D for
visualization. According to the result, our proposed method that use cross
entropy as loss function is able to find the feature matrix more accurately
than the mean squared error function.

##  2\. Related Work

Two representative works of multi-view learning in the early days are
canonical correlation analysis (CCA) (Hotelling 1936) (Hotelling, 1992) and
co-training (Blum and Mitchell 1998) (Blum and Mitchell, 1998). They represent
two core ideas for dealing with multi-view problems through consistency and
complementarity.

Studies in exploiting consistency generally looks for commonalities between
multiple views, which have minimum disagreement. Canonical Correlation
Analysis related algorithms (Chaudhuri et al., 2009; Farquhar et al., 2006;
Hardoon et al., 2004; Kan et al., 2015; Sharma et al., 2012) project two or
more views into latent subspaces by maximizing the correlations among
projected views. Spectral methods (Kumar et al., 2011; De Sa, 2005; Nie et
al., 2016; Xia et al., 2010; Zhou and Burges, 2007b) use a weighted summation
to merge graph Laplacian matrices from different views into one optimal map
for further clustering or embedding. Matrix factorization based methods(Guan
et al., 2015; Liu et al., 2013; Wang et al., 2016) jointly factorize multi-
view data into a common centroid representation by minimizing the overall
reconstruction loss of different views. In addition, multiple kernel learning
(MKL) (Gönen and Alpaydın, 2011) can also be viewed as taking advantage of the
consistency of different views, where each view is mapped to the new space
using kernel tricks, and then all kernel matrices are combined to unify the
kernel by minimizing predefined objective functions.

Another method is to explicitly preserve complementary information of
different views. The co-training algorithms (Blum and Mitchell, 1998; Kumar
and Daumé, 2011; Wang and Zhou, 2010; Zhang and Zhou, 2011) treat each view as
complementary. Generally, it iteratively trains two classifiers on two
different views, and each classifier generates its complementary information
to help other classifiers train in the next iteration.

In summary, most existing multi-view learning algorithms mainly focus on
learning consistency and complementarity from multi-view data. Zhong Zhang et
al. (Zhang et al., 2018) break away from convention and explores the
discriminative and non-discriminative information existing in common and view-
specific parts among different views via joint non-negative matrix
factorization, which provides novel ideas for multi-view learning.

##  3\. The Proposed Method

###  3.1. Non-negative Matrix Factorization

Given a non-negative matrix $X\in\mathbb{R}_{+}^{m\times n}$, which represents
$n$ examples with features. NMF aims to find non-negative matrix factors
$W\in\mathbb{R}_{+}^{m\times k}$ and $H\in\mathbb{R}_{+}^{n\times k}$such
that:

(1) $X\approx W{H^{T}}$

Then the objective function can be written as follows:

(2)
$\displaystyle\underset{W,H}{\mathop{\arg\min}}\,\left\|X-W{{H}^{T}}\right\|_{F}^{2}$
$\displaystyle s.t.\quad W,H\geq 0$

where $|\cdot|{|_{F}}$ denotes the Frobenius norm. Note that the original data
matrix is a linear combination of all column vectors in $W$ with weight of
corresponding column vec-tors in $H$. Therefore, $W$ and $H$ are usually
called the basis matrix and the coefficient matrix respectively. For multi-
view data, the objective function learning representation of NMF-based
approaches is as follows:

(3)
$\displaystyle\mathop{\arg\min}\limits_{W,H}\sum\limits_{v=1}^{{n_{v}}}{\left\|{{X^{(v)}}-W{H^{T}}}\right\|_{F}^{2}+\Phi(W,H)}$
$\displaystyle s.t.\quad W,H\geq 0$

where ${{n}_{v}}$ denotes the number of views, and ${{X}^{(v)}}$ denotes the
data matrix of $v$-th view. $\Phi(\cdot)$ is a regularization term for $W$ and
$H$.

###  3.2. Multi-view Learning via DICS

There are consistency and complementary among different multiple views.
Therefore, DICS algorithm decomposes the multi-view data matrix into two
parts: common part and view-specific parts. As other approaches (Gupta et al.,
2010b; Gupta et al., 2013; Kim et al., 2015; Liu et al., 2014b), define
${{W}_{C}}$ represents the common subspace shared by all views and
$W_{S}^{(v)}$ represents the specific subspace for the $v$-st view. Therefore,
the data matrix of each view can be written as
${{X}^{(v)}}={{W}_{C}}H_{C}^{T}+W_{S}^{(v)}H{{{}_{S}^{(v)}}^{T}}$. In order to
learn the discriminate information, not only the data matrix is divided into
the common and view-specific parts, but also each part of data matrix is
divided into the discriminative and non-discriminative part as follows:

(4) $\begin{array}[]{l}W=[{W_{CD}}\,{W_{CN}}\,W_{SD}^{(v)}\,W_{SN}^{(v)}]\\\
H=[{H_{CD}}\,{H_{CN}}\,H_{SD}^{(v)}\,H_{SN}^{(v)}]\end{array}$

where ${{W}_{CD}}$ and ${{H}_{CD}}$represent the common discriminate matrixes
respectively, while${{W}_{CN}}$ and ${{H}_{CN}}$ separately represent the
common non-discriminate matrixes. Similarly, $W_{SD}^{(v)}$ and $H_{SD}^{(v)}$
respectively indicate the view-specific discriminate parts, while
$W_{SN}^{(v)}$ and $H_{SN}^{(v)}$ separately indicate the view-specific non-
discriminate parts. Afterwards, the discriminate matrixes is used to predict
the output label through a linear projection matrix
$B=[{{B}_{CD}}\,B_{SD}^{(v)}]$, which is corresponding to ${{H}_{CD}}$ and
$H_{SD}^{(v)}$. Therefore, the objection function of DICS is further
reformulated as follows:

(5)
$\displaystyle\mathop{\arg\min}\limits_{W,H,B}\sum\limits_{v=1}^{{n_{v}}}{\left\|{{X^{(v)}}-W{H^{T}}}\right\|_{F}^{2}+\Phi(W,H)}$
$\displaystyle+\gamma\left\|{Y-\left[{{B_{CD}}\,B_{SD}^{(v)}}\right]\left[{\begin{array}[]{*{20}{c}}{H_{CD}^{T}}{H_{SD}^{(v)T}}\end{array}}\right]}\right\|_{F}^{2}$
$\displaystyle s.t.\quad W,H\geq 0$

where $Y\in{{\mathbb{R}}^{c\times n}}$ is the label matrix, c is the number of
classes and n is the number of data instances. Moreover, ${{y}_{i,j}}=1$if the
$j$-st instance belong to class $i$and otherwise ${{y}_{i,j}}=0$.

###  3.3. Improved Objective Function via Cross-Entropy Loss

However, as training results of gradient descent show in section 4 , the mean
squared error term for predicting the output label has a certain degree of
deficiencies. Its partial derivative value is very small when the output
probability value is close to 0 or 1, which could suffer from the partial
derivative value to almost disappear at the beginning of training. In order to
improve the training convergence speed in time with decrease of the error
between predicted and actual values, the cross-entropy loss function is
selected instead of the mean squared error. Cross-entropy loss function is
usually used to measure the distance between two probability distributions
whose input is a probability value between 0 and 1. And cross-entropy loss
increases as the predicted probability diverging from the actual label. For
multi-classification tasks, the form of the cross-entropy loss function is as
follows:

(6) $J=-\sum\limits_{i=1}^{c}{{y_{i}}\log({p_{i}})}$

Where $c$ denotes the number of classes, ${{y}_{i}}$ is an probability
distribution of real labels, if the instance belong to the class i,
${{y}_{i}}$ is 1 otherwise 0. And ${{p}_{i}}$ is the predicted probability
that the observed example belongs to class i. Therefore, the improved
objective function could be written as follows:

(7)
$\displaystyle\mathop{\arg\min}\limits_{W,H,B}\sum\limits_{v=1}^{{n_{v}}}{\left\|{{X^{(v)}}-W{H^{T}}}\right\|_{F}^{2}+\alpha{{\left\|{W_{D}^{T}{W_{D}}}\right\|}_{1,1}}+\beta{{\left\|{{H_{D}}}\right\|}_{1,1}}}$
$\displaystyle-\gamma\sum\limits_{j=1}^{n}{\sum\limits_{i=1}^{c}{{y_{ij}}\ln{p_{ij}}}}$

Where $\alpha,\beta,\gamma$are non-negative parameters to balance the
regularization term and ${{n}_{v}}$ denotes the number of views, $n$ denotes
the number of instances, $c$ denotes the number of classes. the
${{\left\|\,\cdot\,\right\|}_{1,1}}$ is a ${{L}_{1,1}}$ norm constraint on the
discriminative matrix ${{W}_{D}}$, this second term in objective function can
be factorized into two parts:

(8)
${\left\|{W_{D}^{T}{W_{D}}}\right\|_{1,1}}=\sum\nolimits_{i}{w_{Di}^{T}{w_{Di}}}+\sum\nolimits_{i\neq
j}{w_{Di}^{T}{w_{Dj}}}$

the second term in objective function encourage basis vectors to be as
orthogonal as possible, and which reduces the redundancy of discriminative
bases. the third term in objective function is used to prevent overfitting,
What’s more, using ${{L}_{1,1}}$ norm constraint on ${{H}_{D}}$ makes the
discriminate coefficients sparse. The reason is that data points of different
classes should not possess identical basis vectors. ${{p}_{ij}}$ is a
predictive label normalized by softmax function. Specifically, its form is as
follows:

(9)
${p_{ij}}=\frac{{{e^{\sum\nolimits_{k=1}^{{k_{1}}}{{b_{C{D_{i,k}}}}{\rm{\cdot}}{h_{C{D_{j,k}}}}}+\sum\nolimits_{k=1}^{{k_{3}}}{b_{S{D_{i,k}}}^{(v)}{\rm{\cdot}}h_{{}_{S{D_{j,k}}}}^{(v)}}}}}}{{\sum\nolimits_{t=1}^{n}{{e^{\sum\nolimits_{k=1}^{{k_{1}}}{{b_{C{D_{t,k}}}}{\rm{\cdot}}{h_{C{D_{j,k}}}}}+\sum\nolimits_{k=1}^{{k_{3}}}{b_{S{D_{t,k}}}^{(v)}{\rm{\cdot}}h_{{}_{S{D_{j,k}}}}^{(v)}}}}}}}$

Thus the objective function can be written as follows in the form of element-
wise and vectors:

(10)
$\displaystyle\begin{array}[]{l}f(W,H,B)=\sum\limits_{v=1}^{{n_{v}}}{{\rm{||}}{X^{(v)}}-\sum\limits_{i=1}^{k1}{{w_{C{D_{i}}}}h_{C{D_{i}}}^{T}}-\sum\limits_{i=1}^{k2}{{w_{C{N_{i}}}}h_{C{N_{i}}}^{T}}}\\\
\quad\quad\quad\quad-\sum\limits_{i=1}^{k3}{w_{S{D_{i}}}^{(v)}h_{S{D_{i}}}^{(v)T}}-\sum\limits_{i=1}^{k4}{w_{S{N_{i}}}^{(v)}h_{S{N_{i}}}^{(v)T}}{\rm{||}}_{F}^{2}\\\
\quad\quad\quad\quad{\rm{+}}\alpha(\sum\limits_{i=1}^{k1}{\sum\limits_{j=1}^{k1}{w_{C{D_{i}}}^{T}{w_{C{D_{j}}}}}}+\sum\limits_{i=1}^{k3}{\sum\limits_{j=1}^{k3}{w_{S{D_{i}}}^{(v)T}w_{S{D_{j}}}^{(v)}}}\\\
\quad\quad\quad\quad+2\sum\limits_{i=1}^{k1}{\sum\limits_{j=1}^{k3}{w_{C{D_{i}}}^{T}w_{S{D_{j}}}^{(v)}}})\\\
\quad\quad\quad\quad+\beta{1_{1\times
n}}\left({\sum\limits_{i=1}^{k1}{{h_{C{D_{i}}}}}+\sum\limits_{i=1}^{k3}{h_{S{D_{i}}}^{(v)}}}\right)-\gamma\sum\limits_{j=1}^{n}{\sum\limits_{i=1}^{c}{{y_{ij}}\ln{p_{ij}}}}\end{array}$

where the subscript $i$denotes the $i$-th row and the subscript $j$denotes the
$j$-th column of the corresponding matrix. And ${{1}_{1\times n}}$ is a row
vector of length $n$ with all elements 1. By fixing all column vectors except
the one we want to update, we can obtain the convex sub-problem respect to it,
then solve it based on the block coordinate descent framework. In order to
facilitate the derivation, the objective function is divided into two parts:
the nonlinear part of cross-entropy loss function and the other linear parts.
The nonlinear part is written as follows:

(11)
$\begin{array}[]{l}L=-\gamma\sum\limits_{j=1}^{n}{\sum\limits_{i=1}^{c}{{y_{ij}}\ln{p_{ij}}}}\\\
{p_{ij}}=\frac{{{e^{{z_{ij}}}}}}{{\sum\nolimits_{t=1}^{n}{{e^{{z_{tj}}}}}}}\\\
{z_{ij}}=\sum\limits_{k=1}^{{k_{1}}}{{b_{C{D_{i,k}}}}{\rm{\cdot}}{h_{C{D_{j,k}}}}}+\sum\limits_{k=1}^{{k_{3}}}{b_{S{D_{i,k}}}^{(v)}{\rm{\cdot}}h_{{}_{C{D_{j,k}}}}^{(v)}}\end{array}$

Take the partial derivative of L to ${{h}_{C{{D}_{i,j}}}}$as an example:

(12) $\frac{{\partial L}}{{\partial{h_{C{D_{ij}}}}}}=\frac{{\partial
L}}{{\partial{p_{tj}}}}\cdot\frac{{\partial{p_{tj}}}}{{\partial{z_{ij}}}}\cdot\frac{{\partial{z_{ij}}}}{{\partial{h_{C{D_{ij}}}}}}$

and the same process is also suitable to the other parameter, given space
limitations, no further details are given here. As the Fig 1, since the
softmax function contains the output of all previous layers, a partial
derivative is obtained for all n layers.

Figure 1. Schematic diagram of the softmax function

The result of first term is as follows:

(13) $\frac{{\partial
L}}{{\partial{p_{tj}}}}=-\sum\limits_{j=1}^{n}{\sum\limits_{i,t=1}^{c}{\frac{{{y_{ij}}}}{{{p_{tj}}}}}}$

The second term of deriving the softmax function
$\frac{{\partial{p_{tj}}}}{{\partial{z_{ij}}}}$ is more involved and needs to
consider two situations:

(1) Deriving the current node when$t=i$:

(14)
$\begin{array}[]{l}\frac{{\partial{p_{tj}}}}{{\partial{z_{ij}}}}=\frac{{\partial\left({\frac{{{e^{{z_{ij}}}}}}{{\sum\nolimits_{t=1}^{c}{{e^{{z_{tj}}}}}}}}\right)}}{{\partial{z_{ij}}}}=\frac{{{e^{{z_{ij}}}}\cdot\sum\nolimits_{t=1}^{c}{{e^{{z_{tj}}}}}-{e^{{z_{ij}}}}\cdot{e^{{z_{ij}}}}}}{{{{(\sum\nolimits_{t=1}^{c}{{e^{{z_{tj}}}}})}^{2}}}}\\\
\quad\quad={p_{ij}}\cdot(1-{p_{ij}})\end{array}$

(2) Deriving the other nodes when$t\neq i$:

(15)
$\frac{\partial{{p}_{tj}}}{\partial{{z}_{ij}}}=\frac{\partial\left(\frac{{{e}^{{{z}_{tj}}}}}{\sum\nolimits_{t=1}^{c}{{{e}^{{{z}_{tj}}}}}}\right)}{\partial{{z}_{ij}}}=-\frac{{{e}^{{{z}_{tj}}}}}{\sum\nolimits_{t=1}^{c}{{{e}^{{{z}_{tj}}}}}{{}^{2}}}\cdot{{e}^{{{z}_{ij}}}}=-{{p}_{tj}}\cdot{{p}_{ij}}$

We take an example for deriving ${{z}_{2j}}$ as the Fig 2 to explain the two
situations.

Figure 2. Two situations for deriving ${{z}_{2j}}$

And the third term is ${{b}_{C{{D}_{i,j}}}}$. Then combine the results of the
above three items, we have:

(16) $\begin{array}[]{l}\frac{{\partial
L}}{{\partial{h_{C{D_{ij}}}}}}=-\sum\limits_{j=1}^{n}{\sum\limits_{i,t=1}^{c}{{y_{ij}}\cdot\frac{1}{{{p_{tj}}}}\cdot\frac{{\partial{p_{tj}}}}{{\partial{z_{ij}}}}\cdot{b_{C{D_{i,j}}}}}}\\\
=\sum\limits_{j=1}^{n}{\sum\limits_{i,t=1}^{c}{[-\frac{{{y_{ij}}}}{{{p_{ij}}}}\cdot{p_{ij}}\cdot(1-{p_{ij}})+\sum\limits_{t\neq
i}{\frac{{{y_{tj}}}}{{{p_{tj}}}}\cdot{p_{tj}}\cdot{p_{ij}}}]\cdot{b_{C{D_{i,j}}}}}}\\\
=\sum\limits_{j=1}^{n}{\sum\limits_{i,t=1}^{c}{(-{y_{ij}}+{p_{ij}}\sum\limits_{t}{{y_{tj}}})\cdot{b_{C{D_{i,j}}}}}}\end{array}$

Because $Y$ is a label matrix with 0 or 1 elements, And for every sample,
$Y$satisfies the constraint$\sum\limits_{t}{{{y}_{tj}}}=1$. So the partial
derivative of L to ${{h}_{C{{D}_{i,j}}}}$can be written as follows:

(17) $\begin{array}[]{l}\frac{{\partial
L}}{{\partial{h_{C{D_{ij}}}}}}=\sum\limits_{j=1}^{n}{\sum\limits_{i=1}^{c}{({p_{ij}}-{y_{ij}})\cdot{b_{C{D_{i,j}}}}}}\\\
=\sum\limits_{j=1}^{n}{\sum\limits_{i=1}^{c}{(\frac{{{e^{\sum\nolimits_{k=1}^{{k_{1}}}{{b_{C{D_{i,k}}}}{\rm{\cdot}}{h_{C{D_{j,k}}}}}+\sum\nolimits_{k=1}^{{k_{3}}}{b_{S{D_{i,k}}}^{(v)}{\rm{\cdot}}h_{{}_{S{D_{j,k}}}}^{(v)}}}}}}{{\sum\nolimits_{t=1}^{c}{{e^{\sum\nolimits_{k=1}^{{k_{1}}}{{b_{C{D_{t,k}}}}{\rm{\cdot}}{h_{C{D_{j,k}}}}}+\sum\nolimits_{k=1}^{{k_{3}}}{b_{S{D_{t,k}}}^{(v)}{\rm{\cdot}}h_{{}_{S{D_{j,k}}}}^{(v)}}}}}}}}}\\\
\quad\quad\quad\quad-{y_{ij}})\cdot{b_{C{D_{i,j}}}}\end{array}$

Then the partial derivative of the column vector ${{h}_{C{{D}_{i}}}}$ is

(18) $\frac{{\partial
L}}{{\partial{h_{C{D_{i}}}}}}=\sum\limits_{i=1}^{n}{{{({p_{i}}-{y_{i}})}^{T}}\cdot{b_{C{D_{i}}}}}$

where ${{p}_{i}}$ is calculated by the softmax function. Let ${{[\,]}_{+}}$ to
denote $max\left(0,\cdot\right)$, which projects the negative value to the
boundary of feasible region of zero and guarantees the non-negative nature of
the matrix. And we get the update process for each matrix with learning rate
$\eta$ as follows:

(19)
$\begin{array}[]{l}{w_{C{D_{i}}}}={w_{C{D_{i}}}}+\eta[\sum\nolimits_{v=1}^{{n_{v}}}{({R^{(v)}}{h_{C{D_{i}}}}-\alpha({W_{CD}}{1_{k1\times
1}}}\\\ \quad\quad\quad\quad\quad\quad\quad\quad+W_{SD}^{(v)}{1_{k3\times
1}}){]_{+}}\end{array}$ (20)
${w_{C{N_{i}}}}={w_{C{N_{i}}}}+\eta{\left[{\sum\nolimits_{v=1}^{{n_{v}}}{{R^{(v)}}{h_{C{N_{i}}}}}}\right]_{+}}$
(21)
$w_{S{D_{i}}}^{(v)}=w_{S{D_{i}}}^{(v)}+\eta{\left[{{R^{(v)}}h_{S{D_{i}}}^{(v)}-\alpha({W_{CD}}{1_{k1\times
1}}+W_{SD}^{(v)}{1_{k3\times 1}})}\right]_{+}}$ (22)
$w_{S{N_{i}}}^{(v)}=w_{S{N_{i}}}^{(v)}+\eta{\left[{{R^{(v)}}h_{S{N_{i}}}^{(v)}}\right]_{+}}$
(23)
$\begin{array}[]{l}{h_{C{D_{i}}}}={h_{C{D_{i}}}}+\eta[\sum\nolimits_{v=1}^{{n_{v}}}{({R^{(v)}}{w_{C{D_{i}}}}-\frac{\beta}{{\rm{2}}}{1_{n\times
1}})}\\\
\quad\quad\quad\quad\quad\quad\quad-\frac{\gamma}{{\rm{2}}}{Q^{(v)}}^{T}{b_{C{D_{i}}}}){]_{+}}\end{array}$
(24)
${h_{C{N_{i}}}}={h_{C{N_{i}}}}+\eta{\left[{\sum\nolimits_{v}^{{n_{v}}}{{R^{(v)T}}{w_{C{N_{i}}}}}}\right]_{+}}$
(25)
$h_{S{D_{i}}}^{(v)}=h_{S{D_{i}}}^{(v)}+\eta{\left[{{R^{(v)T}}w_{S{D_{i}}}^{(v)}-\frac{\beta}{2}{1_{n\times
1}}-\frac{\gamma}{2}{Q^{(v)T}}b_{SD}^{(v)}}\right]_{+}}$ (26)
$h_{S{N_{i}}}^{(v)}=h_{S{N_{i}}}^{(v)}+\eta{\left[{{R^{(v)T}}w_{S{N_{i}}}^{(v)}}\right]_{+}}$

Where ${{R}^{(v)}}$ and ${{Q}^{(v)}}$ are as follows:

(27)
${R^{(v)}}{\rm{=}}{X^{(v)}}-{W_{CD}}H_{CD}^{T}-{W_{CD}}H_{CD}^{T}-W_{SD}^{(v)}H_{SD}^{(v)T}-W_{SN}^{(v)}H_{SN}^{(v)T}$
(28)
${Q^{(v)}}={\rm{soft}}\max({B_{CD}}H_{CD}^{T}+B_{SD}^{(v)}H_{SD}^{(v)T})-Y$

Where $\eta$ denotes the learning rate, and the second derivative of the
objective function can be used in the actual calculation process, which is
closest to the best step size. Furthermore, when the other variables are
fixed, the projection matrix ${{B}_{CD}}$ and $B_{SD}^{(v)}$ could be
calculated. Let the cross-entropy loss function be zero, and we can get the
projection matrix ${{B}_{CD}}$ and $B_{SD}^{(v)}$ as follows:

(29)
${B_{CD}}=\frac{1}{{{n_{v}}}}\sum\nolimits_{v=1}^{{n_{v}}}{\left({Y-B_{SD}^{(v)}H_{SD}^{(v)T}}\right){H_{CD}}}{\left({H_{CD}^{T}{H_{CD}}+\lambda
I}\right)^{-1}}$ (30)
$B_{SD}^{(v)}=\left({Y-{B_{CD}}H_{CD}^{T}}\right)H_{SD}^{(v)}{\left({H_{SD}^{(v)T}H_{SD}^{(v)}+\lambda
I}\right)^{-1}}$

Where $I$is the identity matrix and $\lambda$ is a small positive number. The
specific processes of the block coordinate descent are depicted in Algorithm
1.

Algorithm 1 MV-DLCSL

0: Multi-view data matrix ${X^{(1)}},{X^{(2)}},\ldots,{X^{({n_{v}})}}$, label
matrix $Y$, parameters $\alpha,\beta,\gamma$, number of latent factors
$k1,k2,k3,k4$;

0: Basis matrixes
$W=\left\\{{{W_{CD}},{W_{CN}},W_{SD}^{(v)},W_{SN}^{(v)}}\right\\}$;Coefficient
matrixes
$H=\left\\{{{H_{CD}},{H_{CN}},H_{SD}^{(v)},H_{SN}^{(v)}}\right\\}$;Projection
matrixes $B=\left\\{{{B_{CD}},B_{SD}^{(v)}}\right\\}$;

1: Initialize matrixes $W,H$ and $B$;

2: repeat

3: Update each column of ${W_{CD}}$ using Eq.(19);

4: Update each column of ${W_{CN}}$ using Eq.(20);

5: for $v=1$ to ${n_{v}}$ do

6: Update each column of $W_{SD}^{(v)}$ using Eq.(21);

7: end for

8: for $v=1$ to ${n_{v}}$ do

9: Update each column of $W_{SN}^{(v)}$ using Eq.(22);

10: end for

11: Update each column of ${H_{CD}}$ using Eq.(23);

12: Update each column of ${H_{CN}}$ using Eq.(24);

13: for $v=1$ to ${n_{v}}$ do

14: Update each column of $H_{SD}^{(v)}$ using Eq.(25);

15: end for

16: for $v=1$ to ${n_{v}}$ do

17: Update each column of $H_{SN}^{(v)}$ using Eq.(26);

18: end for

19: Update ${B_{CD}}$ using Eq.(29);

20: for $v=1$ to ${n_{v}}$ do

21: Update $B_{SD}^{(v)}$ using Eq.(30);

22: end for

23: until convergence or reach max iterations

##  4\. Experiment

In this section, we experimentally evaluate the proposed MV-DLCSL algorithm in
classification task on seven real-world multi-view data sets, and analyze the
convergence of our proposed block coordinate descent algorithm.

###  4.1. Datasets

In this paper, we use seven real-world multi-view data sets to verify the
performance of the proposed algorithm, including Reuters, YaleFace, BBC,
Cornell, Texas, Washington, and Wisconsin datasets. And Cornell, Texas,
Washington, and Wisconsin dataset are four subset of data sets selected from
WebKB data sets. The properties of data sets are summarized in Table 1.

Table 1. Characteristics of the datasets Data Set | Characteristics
---|---
| The numbers of
---
Instances
Views | Classes | Dimensions
Reuters | 1200 | 5 | 6 | 2000 for all
YaleFace | 256 | 2 | 8 | 2016 for all
BBC | 685 | 4 | 5 | 4659/4633/4665/4686
Cornell | 195 | 2 | 5 | 1703/585
Texas | 187 | 2 | 5 | 1703/561
Washington | 230 | 2 | 5 | 1703/690
Wisconsin | 265 | 2 | 8 | 2703/795
Table 2. Accuracy of different methods Method | ACC(%)
---|---
Reuters | YaleFace | BBC | Cornell | Texas | Washington | Wisconsin
GNMF | 40.8$\pm$1.2 | 50.0$\pm$2.5 | 38.0$\pm$1.5 | 41.0$\pm$1.8 | 57.9$\pm$1.8 | 69.6$\pm$2.2 | 52.8$\pm$1.4
MultiNMF | 52.7$\pm$0.2 | 64.2$\pm$4.2 | 73.1$\pm$0.2 | 49.7$\pm$7.7 | 68.7$\pm$3.4 | 59.3$\pm$2.6 | 50.3$\pm$3.5
MVCC | 54.4$\pm$1.9 | 33.3$\pm$6.9 | 95.8$\pm$2.6 | 60.8$\pm$5.0 | 64.7$\pm$5.5 | 62.8$\pm$3.8 | 64.3$\pm$2.7
DICS | 70.3$\pm$4.0 | 89.1$\pm$3.2 | 90.2$\pm$2.4 | 72.8$\pm$6.1 | 81.6$\pm$4.0 | 77.4$\pm$6.0 | 85.1$\pm$4.5
MV-DLCSL | 70.5$\pm$1.3 | 90.2$\pm$1.5 | 91.7$\pm$1.2 | 75.3$\pm$2.2 | 83.8$\pm$2.1 | 83.5$\pm$2.2 | 83.6$\pm$1.9

(a) BBC (b) Reuters (c) YaleFace

Figure 3. Iterative curves on all data sets

###  4.2. Comparison Algorithms

We compare our algorithm with several single-view and multi-view algorithms to
show its effectiveness, including KNN, NMF, SSNMF, GNMF, multiNMF, MVCC, MCL,
and DICS. And the parameters of all algorithms are selected within the range
that the author suggested.

KNN(set $k=1$, i.e., 1-nearest neighbor classifier.) is regarded as the
baseline algorithm and we apply KNN on all single views and report the best
performance on the view. Also we apply the KNN algorithm on the concatenated
feature vector(KNNcat) which is extracted by baseline unsupervised comparison
algorithms .

NMF is applied on each of the single view data and the concatenated feature
vector (i.e. NMFcat) extracted by NMF on the single view data, which is
regarded as inputs another baseline classification algorithm, such as1-nearest
neighbor classifier.

SSNMF is a supervised NMF variant (Lee et al., 2009), which incorporates a
linear classifier to encode the supervised information. We select the
regularization parameter $\lambda$ within the range of [0.5:0.5:3].

GNMF is a manifold regularized version of NMF (Cai et al., 2010), which
preserves the local similarity by imposing a graph Laplacian regularization.
We use the normalized dot product (cosine similarity) to construct the
affinity graph, and select the regularization parameter $\lambda$ within the
set of $\\{{10^{0}},{10^{1}},{10^{2}},{10^{3}},$ ${10^{4}}\\}$.

MultiNMF is an NMF-based multi-view clustering algorithm (Liu et al., 2013),
which can get compatible clustering results across multiple views. We select
the regularization parameter $\lambda$ within the set of $\\{10^{0},10^{-1},$
$10^{-2},10^{-3}\\}$.

MVCC is a novel multi-view clustering method based on concept factorization
with local manifold regularization (Wang et al., 2016), which drives a common
consensus representation for multiple views. We set parameter $\alpha$ to 100,
and select $\beta$ and $\gamma$ within the set of $\\{50,100,200,500,1000\\}$.

MCL is a semi-supervised multi-view NMF variant with graph regularized
constraint (Guan et al., 2015). We select parameter $\alpha$ within the range
of [100:50:250], $\beta$ within the set of ${\rm{\\{}}0.01,{\rm{}}0.02,$
$0.03{\rm{\\}}}$, and set $\gamma$ to 0.005 as author suggested.

DICS is the prototype model of our MV-DLCSL algorithm, which is an NMF-based
multi-view learning algorithm, by exploring the discriminative and non-
discriminative information existing in common and view-specific parts among
different views via joint non-negative matrix factorization, and produce
discriminative and non-discriminative feature from all subspaces. What’s more,
discriminative and non-discriminative features are further used to produce
classification results. We select parameters $\alpha$ and $\beta$ within a
small range of $\left[0,1\right]$, and set parameter $\gamma$ to 1.

###  4.3. Result

For all algorithms, we first perform a five-folds cross validation to select
the parame-ters that has the best accuracies and generalization performance.
Due to randomness, we run all algorithms 10 times on each dataset and report
the mean values and standard deviation of accuracies. All the classification
results of seven multi-view datasets are summarized in Table 2, and the best
result on each dataset is highlighted in boldface. As we can see, the proposed
algorithm achieves better accuracy on most of the datasets, and slightly worse
than other algorithms on BBC and Wisconsin datasets, and It is worth
mentioning that, the standard deviations of accuracies for MV-DLCSL are much
lower than DICS method.

(a) DICS on BBC data sets (left panel) and our MV-DLCSL (right panel) (b) DICS
on Cornell data sets (left panel) and our MV-DLCSL (right panel) (c) DICS on
Reuters data sets (left panel) and our MV-DLCSL (right panel) (d) DICS on
Texas data sets (left panel) and our MV-DLCSL (right panel)

(e) DICS on Washington data sets (left panel) and our MV-DLCSL (right panel)
(f) DICS on Wisconsin data sets (left panel) and our MV-DLCSL (right panel)
(g) DICS on YaleFace data sets (left panel) and our MV-DLCSL (right panel)

Figure 4. t-SNE scatter plot of the feature matrix on all data sets

###  4.4. Convergence Analysis

In order to empirically investigate the convergence property of our algorithm,
we plot the iterative curves of objective function on five typical data sets
in Fig 3. From Fig 3, we can observe that the objective function values drop
sharply and then iterative curves begin to grow/decrease mildly, then it
converges eventually. Usually, the algo-rithm will converge in no more than
100 iterations.

###  4.5. Discriminant Matrix Visualization

t-SNE (t-distributed stochastic neighbour embedding) is a machine learning
algorithm used for dimensionality reduction (Maaten and Hinton, 2008). It was
proposed by Laurens van der Maaten and Geoffrey Hinton in 2008. In addition,
t-SNE is a nonlinear dimensionality reduction algorithm, which is very
suitable for high-dimensional data dimensionality reduction to 2D or 3D for
visualization. We utilize DICS and our MV-DLCSL to obtain the basis matrix
${{W}_{D}}$ and the coefficient matrix ${{H}_{D}}$ through iterative updating
processes, which is used to generate the feature matrix of raw data matrix$X$.
By using the t-SNE for dimension reduction processing of the feature matrix,
it is possible to visually see the distinguishable feature matrix. Therefore,
we compare DICS with our MV-DLCSL method by two-dimensional visualized feature
matrix to discern the classification capacities of feature matrices extracting
by two methods. In the scatter plot generated by t-SNE, each point represents
an instance and different colors indicate different class labels. Therefore,
most points with the same color are clustered into clusters in the scatter
plot. From Fig 4, we could see the t-SNE scatter plot of feature matrices that
our MV-DLCSL method extract having better clustering results, and there are
fewer instances assigned to the wrong cluster. What’s more, the clusters in
our MV-DLCSL method are more concentrated than DICS. Thus, our method is able
to find the feature matrix more accurately than DICS.

##  5\. Conclusions and Future Work

In this paper, we propose a novel multi-view network MV-DLCSL. The proposed
algorithm explores the discriminative and non-discriminative information
existing in common and view-specific parts among different views via joint
non-negative matrix factorization, and use the cross entropy loss to constrain
the objective function, which shows better classification results than the
mean square error. The experimental results on seven real-world data sets have
demonstrated the effectiveness of our proposed algorithm. For future studies,
we plan to consider the distance farther between the view-specific matrixes to
grasp the information related to the common and view-specific part, which is
benefit to precisely classify instances through properties between different
views.

## References

  * (1)
  * Blum and Mitchell (1998) Avrim Blum and Tom Mitchell. 1998.  Combining labeled and unlabeled data with co-training. In Proceedings of the eleventh annual conference on Computational learning theory. ACM, 92–100.
  * Cai et al. (2010) Deng Cai, Xiaofei He, Jiawei Han, and Thomas S Huang. 2010\.  Graph regularized nonnegative matrix factorization for data representation.  IEEE transactions on pattern analysis and machine intelligence 33, 8 (2010), 1548–1560.
  * Chaudhuri et al. (2009) Kamalika Chaudhuri, Sham M Kakade, Karen Livescu, and Karthik Sridharan. 2009.  Multi-view clustering via canonical correlation analysis. In Proceedings of the 26th annual international conference on machine learning. ACM, 129–136.
  * De Sa (2005) Virginia R De Sa. 2005\.  Spectral clustering with two views. In ICML workshop on learning with multiple views. 20–27.
  * Farquhar et al. (2006) Jason Farquhar, David Hardoon, Hongying Meng, John S Shawe-taylor, and Sandor Szedmak. 2006\.  Two view learning: SVM-2K, theory and practice. In Advances in neural information processing systems. 355–362.
  * Gönen and Alpaydın (2011) Mehmet Gönen and Ethem Alpaydın. 2011.  Multiple kernel learning algorithms.  Journal of machine learning research 12, Jul (2011), 2211–2268.
  * Guan et al. (2015) Ziyu Guan, Lijun Zhang, Jinye Peng, and Jianping Fan. 2015\.  Multi-view concept learning for data representation.  IEEE Transactions on Knowledge and Data Engineering 27, 11 (2015), 3016–3028.
  * Gupta et al. (2010a) Sunil Kumar Gupta, Dinh Phung, Brett Adams, Truyen Tran, and Svetha Venkatesh. 2010a.  Nonnegative shared subspace learning and its application to social media retrieval. In Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 1169–1178.
  * Gupta et al. (2010b) Sunil Kumar Gupta, Dinh Phung, Brett Adams, Truyen Tran, and Svetha Venkatesh. 2010b.  Nonnegative shared subspace learning and its application to social media retrieval. In Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 1169–1178.
  * Gupta et al. (2013) Sunil Kumar Gupta, Dinh Phung, Brett Adams, and Svetha Venkatesh. 2013.  Regularized nonnegative shared subspace learning.  Data mining and knowledge discovery 26, 1 (2013), 57–97.
  * Hardoon et al. (2004) David R Hardoon, Sandor Szedmak, and John Shawe-Taylor. 2004\.  Canonical correlation analysis: An overview with application to learning methods.  Neural computation 16, 12 (2004), 2639–2664.
  * Hotelling (1992) Harold Hotelling. 1992\.  Relations between two sets of variates.  In Breakthroughs in statistics. Springer, 162–190.
  * Kan et al. (2015) Meina Kan, Shiguang Shan, Haihong Zhang, Shihong Lao, and Xilin Chen. 2015.  Multi-view discriminant analysis.  IEEE transactions on pattern analysis and machine intelligence 38, 1 (2015), 188–194.
  * Kim et al. (2015) Hannah Kim, Jaegul Choo, Jingu Kim, Chandan K Reddy, and Haesun Park. 2015.  Simultaneous discovery of common and discriminative topics via joint nonnegative matrix factorization. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 567–576.
  * Kim et al. (2014) Jingu Kim, Yunlong He, and Haesun Park. 2014.  Algorithms for nonnegative matrix and tensor factorizations: A unified view based on block coordinate descent framework.  Journal of Global Optimization 58, 2 (2014), 285–319.
  * Kumar and Daumé (2011) Abhishek Kumar and Hal Daumé. 2011.  A co-training approach for multi-view spectral clustering. In Proceedings of the 28th International Conference on Machine Learning (ICML-11). 393–400.
  * Kumar et al. (2011) Abhishek Kumar, Piyush Rai, and Hal Daume. 2011.  Co-regularized multi-view spectral clustering. In Advances in neural information processing systems. 1413–1421.
  * Lee et al. (2009) Hyekyoung Lee, Jiho Yoo, and Seungjin Choi. 2009.  Semi-supervised nonnegative matrix factorization.  IEEE Signal Processing Letters 17, 1 (2009), 4–7.
  * Liu et al. (2014a) Jing Liu, Yu Jiang, Zechao Li, Zhi-Hua Zhou, and Hanqing Lu. 2014a.  Partially shared latent factor learning with multiview data.  IEEE transactions on neural networks and learning systems 26, 6 (2014), 1233–1246.
  * Liu et al. (2014b) Jing Liu, Yu Jiang, Zechao Li, Zhi-Hua Zhou, and Hanqing Lu. 2014b.  Partially shared latent factor learning with multiview data.  IEEE transactions on neural networks and learning systems 26, 6 (2014), 1233–1246.
  * Liu et al. (2013) Jialu Liu, Chi Wang, Jing Gao, and Jiawei Han. 2013\.  Multi-view clustering via joint nonnegative matrix factorization. In Proceedings of the 2013 SIAM International Conference on Data Mining. SIAM, 252–260.
  * Maaten and Hinton (2008) Laurens van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-SNE.  Journal of machine learning research 9, Nov (2008), 2579–2605.
  * Nie et al. (2016) Feiping Nie, Jing Li, Xuelong Li, and others. 2016\.  Parameter-Free Auto-Weighted Multiple Graph Learning: A Framework for Multiview Clustering and Semi-Supervised Classification.. In IJCAI. 1881–1887.
  * Shao et al. (2018) Junming Shao, Chun Meng, Masoud Tahmasian, Felix Brandl, Qinli Yang, Guangchun Luo, Cheng Luo, Dezhong Yao, Lianli Gao, Valentin Riedl, and others. 2018.  Common and distinct changes of default mode and salience network in schizophrenia and major depression.  Brain imaging and behavior (2018), 1–12.
  * Shao et al. (2012a) Junming Shao, Nicholas Myers, Qinli Yang, Jing Feng, Claudia Plant, Christian Böhm, Hans Förstl, Alexander Kurz, Claus Zimmer, Chun Meng, and others. 2012a.  Prediction of Alzheimer’s disease using individual structural connectivity networks.  Neurobiology of aging 33, 12 (2012), 2756–2765.
  * Shao et al. (2012b) Junming Shao, Qinli Yang, Afra Wohlschlaeger, and Christian Sorg. 2012b.  Discovering aberrant patterns of human connectome in alzheimer’s disease via subgraph mining. In 2012 IEEE 12th International Conference on Data Mining Workshops. IEEE, 86–93.
  * Sharma et al. (2012) Abhishek Sharma, Abhishek Kumar, Hal Daume, and David W Jacobs. 2012\.  Generalized multiview analysis: A discriminative latent space. In 2012 IEEE Conference on Computer Vision and Pattern Recognition. IEEE, 2160–2167.
  * Wang et al. (2016) Hao Wang, Yan Yang, and Tianrui Li. 2016.  Multi-view clustering via concept factorization with local manifold regularization. In 2016 IEEE 16th International Conference on Data Mining (ICDM). IEEE, 1245–1250.
  * Wang and Zhou (2010) Wei Wang and Zhi-Hua Zhou. 2010.  A New Analysis of Co-Training.. In ICML, Vol. 2. 3\.
  * Xia et al. (2010) Tian Xia, Dacheng Tao, Tao Mei, and Yongdong Zhang. 2010\.  Multiview spectral embedding.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 40, 6 (2010), 1438–1446.
  * Xu et al. (2013) Chang Xu, Dacheng Tao, and Chao Xu. 2013.  A survey on multi-view learning.  arXiv preprint arXiv:1304.5634 (2013).
  * Ye et al. (2015) Han-Jia Ye, De-Chuan Zhan, Yuan Miao, Yuan Jiang, and Zhi-Hua Zhou. 2015.  Rank consistency based multi-view learning: a privacy-preserving approach. In Proceedings of the 24th ACM International on Conference on Information and Knowledge Management. ACM, 991–1000.
  * Zhang and Zhou (2011) Min-Ling Zhang and Zhi-Hua Zhou. 2011.  CoTrade: confident co-training with data editing.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 41, 6 (2011), 1612–1626.
  * Zhang et al. (2018) Zhong Zhang, Zhili Qin, Peiyan Li, Qinli Yang, and Junming Shao. 2018.  Multi-view Discriminative Learning via Joint Non-negative Matrix Factorization. In International Conference on Database Systems for Advanced Applications. Springer, 542–557.
  * Zhou and Burges (2007a) Dengyong Zhou and Christopher JC Burges. 2007a.  Spectral clustering and transductive learning with multiple views. In Proceedings of the 24th international conference on Machine learning. ACM, 1159–1166.
  * Zhou and Burges (2007b) Dengyong Zhou and Christopher JC Burges. 2007b.  Spectral clustering and transductive learning with multiple views. In Proceedings of the 24th international conference on Machine learning. ACM, 1159–1166.


# The role of solvent interfacial structural ordering in maintaining stable
graphene dispersions

Urvesh Patil1 and Nuala M. Caffrey2,3 1 School of Physics & CRANN Institute,
Trinity College Dublin, Dublin 2, Ireland  2 School of Physics, University
College Dublin, Dublin 4, Ireland  3 Centre for Quantum Engineering, Science,
and Technology, University College Dublin, Dublin 4, Ireland
<EMAIL_ADDRESS>
###### Abstract

Liquid phase exfoliation (LPE) is the most promising method for the low-cost,
scalable production of two-dimensional nanosheets from their bulk
counterparts. Extensive exfoliation occurs in most solvents due to the huge
amount of energy introduced by sonication or shear mixing. However, the
subsequent dispersion is not always stable, with extensive reaggregation
occurring in some solvents. Identifying the optimal solvent for a particular
layered material is difficult and requires a fundamental understanding of the
mechanism involved in maintaining a stable dispersion. Here, we use molecular
dynamics calculations to show that when graphene is immersed in a solvent,
distinct solvation layers are formed irrespective of the choice of solvent and
their formation is energetically favourable for all considered solvents.
However, energetic considerations such as these do not explain the
experimental solvent-dependence of the dispersion concentration. Instead, we
find that solvents with high diffusion coefficients parallel to the graphene
layer result in the lowest experimental concentration of graphene in solution.
This can be explained by the enhanced ease of reaggregation in these solvents.
Solvents with smaller diffusion coefficients result in higher experimental
graphene concentrations as reaggregation is prevented. In the low diffusion
limit, however, this relationship breaks down. We suggest that here the
concentration of graphene in solution depends primarily on the separation
efficiency of the initial exfoliation step. Based on this, we predict that the
concentration of exfoliated graphene in solvents such as benzaldehyde and
quinoline, which have low diffusion constants, can be increased dramatically
by careful tuning of the experimental sonication parameters.

liquid phase exfoliation; graphene, molecular dynamics; diffusion

##  1 Introduction

The liquid-phase exfoliation (LPE) of layered materials is a scalable top-down
method to produce industrial quantities of monolayer or few-layer two-
dimensional (2D) nanosheets at a reasonable cost [1, 2, 3, 4]. Solvent-
dispersed 2D nanosheets can find applications as inks for printed flexible
electronics [5, 6], energy generation and storage [7], sensing [8] and fire
resistant coatings [9]. The ability of LPE to produce large quantities
compensates for relatively small lateral size flakes and low yield achieved
compared to other exfoliation methods [10]. However, the quality and quantity
of the exfoliated layers depends dramatically on the solvent used [11, 12]. A
careful choice of both the solvent and the sonication parameters can improve
the monolayer yield and increase the lateral dimensions achieved, but a
fundamental lack of understanding of the mechanism driving the exfoliation and
stabilization processes impedes this optimisation procedure [13, 14].

These two processes, exfoliation and stabilization, can be considered
separately even though they occur effectively simultaneously experimentally.
Exfoliation occurs when the van der Waals interactions binding the layers
together are overcome via large shear forces that are introduced in the
presence of a solvent through either sonication, high-shear mixing or wet-ball
milling [10]. The solvent has been suggested to play a relatively minor role
in this initial separation of layers, with extensive exfoliation occurring in
most solvents due to the huge amounts of energy involved [15, 16]. The solvent
then stabilizes the exfoliated nanosheets preventing reaggregation and
sedimentation. The choice of solvent now plays a critical role. Solvent
exchange methods, where different solvents are used in the exfoliation and
stabilization steps, find that those solvents which do not perform well in the
initial exfoliation step may perform very well in maintaining a stable
dispersion of already exfoliated monolayers [17, 18, 19]. Evidently, the
solvent plays a different role in the two steps but as they occur effectively
simultaneously, their exact role is difficult to disentangle.

The majority of investigations in this regard concern the exfoliation of
graphite. Hernandez et al. compared the performance of a wide selection of
solvents and found that the maximum graphene concentration was achieved by
exfoliation in cyclopentanone, resulting in 8.5 $\pm$ 1.2 $\mu$g/ml [20].
Exfoliation in non-polar solvents, including toluene, heptane, hexane and
pentane, resulted in much lower graphene concentrations of 0.8 $\pm$ 0.4
$\mu$g/ml, 0.3 $\pm$ 0.4 $\mu$g/ml, 0.2 $\pm$ 0.1 $\mu$g/ml and 0.16 $\pm$
0.05 $\mu$g/ml, respectively. However, some highly polar molecules such as
water and formamide also perform poorly, with the best performing solvents
having a slightly polar nature (although not all slightly polar solvents
perform well). The effectiveness of LPE for other layered materials, such as
WS2, MoS2 and h-BN, is also strongly dependent on the choice of solvent [2,
21].

Clearly, the inability to identify the optimal solvent for the exfoliation of
a particular material without resorting to expensive trial-and-error is a
limitation for its extensive industrial use. Understanding the nature of the
interactions between the layered material and the solvent would allow for the
preemptive screening of effective solvents for the exfoliation of layered
materials. As well as aiming to maximize monolayer yield, there are other
restrictions on the choice of solvent. For example, many of the better
solvents for graphite exfoliation, including NMP, DMA and DMF, are toxic and
facing restrictions by the European Chemicals Agency, rendering them
unsustainable for future industrial use [22]. Devising a screening tool to
identify non-toxic, low-boiling point solvents which do not comprise on
performance would be extremely useful [23, 24]. Ideally these screening tools
would be generalizable to other two-dimensional materials beyond graphene.

One of the first attempts to screen for effective solvents found that the
concentration of dispersed graphene was maximized for those solvents with a
surface energy or tension similar to that of graphite, following the
traditional surface wettability argument [12]. Although this method has had
some successes, there are also examples where it breaks down [20]. Hansen
solubility parameters, characterized by three intermolecular interactions,
namely hydrogen bonding, polar and nonpolar (dispersive), have also been used
as screening parameters. Graphene was found to have Hansen dispersion,
hydrogen bonding and polar solubility parameters of $\delta_{D}$ = 18 MPa1/2 ,
$\delta_{P}$ = 9.3 MPa1/2 , and $\delta_{H}$ = 7.7 MPa1/2, respectively [20].
Searching for solvents with solubility parameters close to those of graphene
led to the prediction and confirmation of cyclopenanone and cyclohexanone as
excellent solvents of graphene. However, here again, there are some examples
where this screening technique fails – dimethyl phthalate has very similar
Hansen parameters to cyclopentanone, yet is a poor graphene solvent [20]. A
peculiarity of the Hansen solubility parameters attributed to graphene is the
non-zero polar component ($\delta_{P}=9.3$ MPa1/2). The mechanism which could
lead to this is so-far unknown. Contribution from edge-sites or the
unintentional chemical functionalisation of the basal plane were touted as a
potential physical mechanisms, yet there is little experimental evidence for
either [25, 26].

All of the above suggests that empirical solubility models, which consider
only macroscopic solution thermodynamics, are blunt tools, missing some
information about the stabilization mechanism [27, 28]. These methods assume
that kinetic effects driving reaggregation play a minor role. There is growing
acceptance that this is unlikely to be true and that explicit structural and
electronic interactions between the solvent molecules and the solute are
important in the stabilization of the exfoliated monolayers [29, 30, 31, 32].

Here, we systematically study the origin of the solvent dependence of the
stabilization of exfoliated monolayers. We address several open questions: Is
it sufficient to consider only energetic effects via the enthalpy of mixing
when choosing screening descriptors, what is the physical mechanism leading to
the experimental finding of a non-zero polar component of the Hansen
solubility parameter of graphene and, and how important are kinetic effects,
such as reaggregation, to the stabilization of graphene in a particular
solvent?

##  2 Methods

###  2.1 Choice of solvents

The twelve solvents included in this investigation are listed in Table 1. They
were chosen to include those reported by Hernandez et al. [33] to result in a
wide range of graphene concentrations after exfoliation and include
cyclopentanone (best performing solvent, polar aprotic), N-methyl-pyrrolidone
(commonly used effective solvent, polar aprotic), benzaldehyde (poor
performance, polar aprotic), bromobenzene (intermediate performance, slightly
non-polar) and toluene (poor performance, non-polar). All others have
intermediate exfoliation capabilities or polar properties. Note that while the
concentration of exfoliated graphene in individual solvents has been increased
considerably since the initial work of Hernandez et al, it is, to the best of
our knowledge, the only available comparative study of the effectiveness of a
wide range of solvents for graphite exfoliation.

###  2.2 Molecular Dynamics

Molecular dynamics calculations are performed using the Large-scale Atomic /
Molecular Massively Parallel Simulator (LAMMPS) [34]. As the interaction
between graphene and isolated solvent molecules was found to be primarily
attractive dispersion forces [35, 36], it can be modelled accurately using
classical van der Waals (vdW) force fields. The graphene monolayer is modeled
as uncharged vdW spheres and are described by interaction parameters
originally reported by Steele et al. [37] and later used to describe the
interaction between graphene and molecular adsorbants [38, 30, 39, 40]. The
molecules are described by the All Atom Optimized Potentials for Liquid
Simulations (OPLS-AA) potentials, which were obtained from LigPargen web
server [41, 42, 43]. The pair interaction coefficients between graphene and
the molecules are obtained using Lorentz-Berthelot rules [44, 45]. Binding
energies calculated using these potentials agree with experiment for small
organic molecules adsorbed on graphene [35]. The initial solvent structure was
created using packmol [46] in a box size of 105.91 Å $\times$ 106.65 Å
$\times$ 120.0 Å to reproduce the experimental room-temperature densities (see
Table 3). The graphene monolayer was then included in the simulation box by
increasing the box size in the direction normal to the plane by (2 $\times$
3.5) Å. The python package pymatgen was used to generate the LAMMPS structural
data file [47]. The systems were annealed from 600 K to 300 K for 50 fs with
an integration timestep of 0.01 fs and then equilibriated at 300 K for 1 ns
with an integration timestep of 1 fs. The system was further simulated with
the NVE ensemble for 2 ns, of which 2000 samples from the last 1 ns were used
to generate statistical distribution functions. For the duration of the
complete simulation the graphene layer is held fixed.

###  2.3 Thermodynamics of Mixing

The Helmholtz free energy of mixing, $\Delta A_{\mathrm{mix}}$, is an
indication of the favorable formation of a solvation structure. As $\Delta
A_{\mathrm{mix}}=\Delta H_{\mathrm{mix}}-T\Delta S_{\mathrm{mix}}$, where
$\Delta H_{\mathrm{mix}}$ and $\Delta S_{\mathrm{mix}}$ are the enthalpy and
entropy of mixing, respectively, a requirement of mixing is that $\Delta
H_{\mathrm{mix}}<T\Delta S_{\mathrm{mix}}$. Note that the PV term of the
enthalpy is negligible for liquids and solids so it can be neglected. In this
case, the internal energy of mixing at 0 K will be equal to the enthalpy of
mixing at 0 K. $\Delta A_{\mathrm{mix}}$ was determined from molecular
dynamics simulations using achemical free energy methods [48]. Such methods
simulate a series of non-physical intermediates to calculate the free energy
of transferring a solute (here graphene) from vacuum (state 0) to solution
(state 1). A continuous variable $\lambda$ parameterizes the path between
state $0$ and state $1$. $\Delta A_{\mathrm{mix}}$ is calculated here using
the finite difference thermodynamic integration (FDTI) method [49, 50, 51]. MD
simulations are carried out for a discrete set of $\lambda$ values and the
free energy of mixing is calculated as:

$\Delta A=\int_{\lambda=0}^{\lambda=1}\biggl{<}\frac{\partial
U}{\partial\lambda}\biggr{>}_{\lambda}d\lambda\approx\sum_{i=0}^{N-1}w_{i}\biggl{<}\frac{U(\lambda_{i}+\delta)-U(\lambda_{i})}{\delta}\biggr{>}$

where $\langle\dots\rangle$ indicates an ensemble average, $w_{i}$ are
integral weights and $\delta$ is the finite difference parameter which
satisfies $(\lambda_{i}-\lambda_{i-1})>>\delta$. To avoid numerical
instabilities the Leonard-Jones potential is modified to be a function of
$\lambda$. The system is simulated for 10 ns in the NVT ensemble at 300 K. The
parameter $\lambda$ is incremented in 40 steps (every 250 ps) with the finite
difference parameter $\delta=0.001$. For each $\lambda$, the first 100 ps are
discarded for equilibration and the ensemble average is calculated by sampling
every 20 fs from the last 150 ps. The enthalpy and entropy contributions to
the total free energy are extracted from the slope and $y$-intercept of
$\Delta A/T$ versus $1/T$, respectively, where the free energy is calculated
at 290 K, 300 K and 310 K.

###  2.4 Diffusion Coefficients

The diffusion coefficient of solvent molecules confined in a region $[a,b]$
parallel to the graphene sheet, $D_{xx}$, was determined from MD trajectories
using the method described in Liu et al. [52] and Agosta et al. [53]. $D_{xx}$
is calculated as:

$D_{xx}(\\{a,b\\})=\lim_{\tau\rightarrow\infty}\frac{\langle\Delta
x(\tau)^{2}\rangle_{\\{a,b\\}}}{2\,\tau\,P(\tau)}$

where the terms on the right hand side can all be calculated using standard
molecular dynamics calculations. $\langle\Delta
x(\tau)^{2}\rangle_{\\{a,b\\}}$ is the mean square displacement of particles
in the direction parallel to the surface in a region $\\{a,b\\}$ in a time
window $\tau$ and $P(\tau)$ is the probability that a particle remains in the
region $\\{a,b\\}$ during that time, i.e, the survival probability.
$\langle\Delta x(\tau)^{2}\rangle_{\\{a,b\\}}$ is extracted from MD
trajectories via:

$\langle\Delta
x(\tau)^{2}\rangle_{a,b}=\frac{1}{T}\sum_{t=1}^{T}\frac{1}{N(t)}\sum_{i\in\zeta(t,t+\tau)}(x_{i}(t+\tau)-x_{i}(t))^{2}$

where $\zeta(t,t+\tau)$ is the set of all particles that stay in the volume
$a<z<b$ from time $t$ to $t+\tau$, $N(t)$ is the number of particles in the
region $\\{a,b\\}$ at time $t$, and $T$ is the number of time steps to
average. The survival probability is calculated as:

$P(\tau)=\frac{1}{T}\sum_{t}^{T}\frac{N(t,t+\tau)}{N(t)}$

where $N(t,t+\tau)$ is the number of particles that stay within the region
$\\{a,b\\}$ between times $t$ and $t+\tau$ . Note that the parameters in the
calculation of diffusion coefficient are $\tau$, $T$, and the definition of
the region between $a$ and $b$. Here, the equilibrated structure was simulated
with the NVE ensemble for $T=150$ ps and the center of mass of each molecule
recorded at each time step. A time window ($\tau$) of 35 ps is used to
calculate the diffusion coefficient. $a$ and $b$ are chosen to encompass the
entire first solvation shell of each system.

##  3 Results

In Section 3.1, we show that the free energy of mixing is always negative
regardless of the solvent due to a large enthalphic contribution for solvating
graphene. A significant entropic penalty, however, suggests that the solvent
must undergo significant ordering in the presence of graphene. Section 3.2
discusses this structural and orientational ordering. Finally, Section 3.3
presents the diffusion coefficients of solvent molecules in the first
solvation shell and discusses how these molecules play a major role in
preventing the reaggregation of graphene sheets.

###  3.1 Free Energy of Mixing

It was previously shown, using _ab initio_ calculations, that the interaction
between an isolated gas phase solvent molecule and a pristine graphene
monolayer is van der Waals in nature, with no charge transfer involved. This
is independent of the type of molecule, and the binding strength of a solvent
molecule on graphene is not correlated with its ability to exfoliate graphite
[36]. To determine if instead the collective behaviour of the molecules near
the surface of a graphene monolayer is the determining quantity, we calculate
the Helmholtz free energy of mixing, $\Delta A_{\mathrm{mix}}$, of graphene in
a solvent. If $\Delta A_{\mathrm{mix}}<0$, then the mixture is
thermodynamically more stable than each of the two components separately. The
Helmholtz free energies of mixing for graphene at room temperature are shown
in the first column of Table 1 for a variety of solvents.

Solvent |  $\Delta A_{\mathrm{mix}}$ (MJ/mol) |  $\Delta H_{\mathrm{mix}}$ (MJ/mol) |  $T\Delta S_{\mathrm{mix}}$ (MJ/mol)
---|---|---|---
Quinoline | -11.08 | -28.32 | -17.13
Bromobenzene | -10.41 | -9.57 | 0.83
NMP | -10.172 | -28.15 | -18.01
Benzaldehyde | -9.90 | -26.57 | -16.64
Chlorobenzene | -9.65 | -23.10 | -13.45
1,3-dioxolane | -9.22 | -25.23 | -16.0
Toluene | -8.99 | -21.10 | -12.07
Benzonitrile | -8.94 | -21.98 | -13.06
Cyclopentanone | -8.33 | -24.07 | -15.74
Cyclohexanone | -7.57 | -22.26 | -14.69
Table 1:  The Helmholtz free energy of mixing, $\Delta A_{\mathrm{mix}}$, the
corresponding enthalpy, $\Delta H_{\mathrm{mix}}$ and the entropy term
calculated at room temperature (300 K), $T\Delta S_{\mathrm{mix}}$. These
values are determined for a graphene area of $112.953$ nm2.

At 300 K, $\Delta A_{\mathrm{mix}}$ is always negative, irrespective of the
solvent, suggesting the favorable mixing of monolayer graphene in each of the
solvents. It is most favourable in quinoline ($\Delta A_{\mathrm{mix}}=-11.08$
MJ/mol) followed by bromobenzene ($\Delta A_{\mathrm{mix}}=-11.41$ MJ/mol). It
is least favourable in cyclohexanone, at $-7.57$ MJ/mol. The magnitudes of
$\Delta A_{\mathrm{mix}}$ reported here are consistent with those reported by
Oyer et al. for the free energy of mixing of graphene in benzene and
hexaflourene [54] and by Mukhopadhyay et al. for the free energy of mixing of
graphene in water and DMSO [55].

The decomposition of the Helmholtz energy into the enthalpic and entropic
contributions are given in the second and third column of Table 1,
respectively. In almost all cases there is a significant entropic penalty for
solvating graphene, with $T\Delta S_{\mathrm{mix}}$ lying between $-17.13$
MJ/mol and $-12.07$ MJ/mol. This entropic penalty is compensated by a large
enthalpy change ranging between $-28.32$ MJ/mol in quinoline to $-21.10$
MJ/mol in toluene, driven by strong van der Waals interactions. Bromobenzene
is an exception to this: $\Delta H_{mix}$ is much smaller in magnitude than
for the other solvents at $-9.57$ MJ/mol, while the entropy contribution is
positive at 0.83 MJ/mol. However, this unusual behavior is not evident in the
overall Helmoltz free energy of mixing which is similar to that of the other
considered solvents, at $-10.41$ MJ/mol.

In general, no correlation can be found between the magnitude of $\Delta
A_{\mathrm{mix}}$ and the experimentally-determined graphene concentrations
after exfoliation. As such, the Helmoltz free energy of mixing cannot be used
as a mechanism of screening for more effective solvents. Instead, the
calculated free energies indicate that the solvent-dependent exfoliation is
not energetically bound and other dynamical effects should be considered.

###  3.2 Solvation Shell Formation

The entropic penalty suggests that there must be some solvent reconstruction
or ordering that occurs due to the presence of graphene. To understand this,
the location and orientation of solvent molecules close to the graphene layer
are extracted from molecular dynamics calculations. The position distribution
function $g(z)$ and the angle distribution function $\alpha(z)$ are shown in
Fig. 1. The position distribution function $g(z)$ is given relative to the
bulk solvent number density so that a value greater than 1 indicates an
accumulation of molecules relative to the bulk solvent, and a value less than
1 indicates a depletion of molecules relative to the bulk. $\alpha(z)$ is the
angle made by the normal to the plane of the molecule with respect to the $z$
axis, where the plane is defined as containing at least 3 atoms of the
molecular ring system.

Figure 1: Pair distribution functions. (Top panel) $g(z)$ is the position
distribution of the center of mass of the solvent molecules. $g(z)$ is given
relative to the bulk solvent number density. (Bottom panel) $\alpha(z)$ is the
probability distribution function of the angle made by the normal to the plane
of molecule with the axis perpendicular to the plane of graphene. The graphene
monolayer is located at $z=0$ Å.

Irrespective of the polar nature of the molecule, distinct solvation shells
are formed next to the graphene layer as a result of surface confinement,
molecule - molecule interactions and molecule - graphene interactions. The
first solvation shell (grey shaded area) can be recognized as a peak in $g(z)$
followed by a deep trough near the surface of graphene. The peak associated
with the first solvation shell is also the sharpest, indicating that the
highest degree of transverse ordering occurs for those closest to graphene.
The first solvation shell has a complex structure in most cases due to the
adsorption geometry of the molecules on the graphene surface. This is
particularly evident for cyclopentanone (Fig 1(a)) and NMP (Fig 1(f)) where
the first solvation shell is composed of two distinct peaks. The reason for
this can be found by looking at $\alpha(z)$: the solvent adsorbs on graphene
with two distinct angular orientations – one with the molecular plane lying
approximately parallel to the surface ($\alpha=0^{\circ}$) and one with the
plane almost perpendicular ($\alpha=90^{\circ}$). As the distribution function
is based on the centre of mass of the molecules, this manifests as a double
peak structure in $g(z)$. As a result, the dipole moments of the molecules
also show two distinct polar angles $(\Theta(z))$, whereas the azimuthal angle
$(\phi(z))$ shows no distinct orientation preference indicating no in-plane
ordering of the dipoles (see Fig 5). This behavior was previously found in
molecular dynamic calculations involving NMP interacting with both graphene
and carbon nanotubes [56, 39, 57]. In other cases, the distinction between the
two peaks is not as clear. For example, the first solvation shell of
cyclohexanone (Fig 1(b)) has a small peak corresponding to those molecules
located closest to graphene and orientated parallel to it, but the remainder
of the molecules in the first solvation shell do not exhibit an angular
preference. This is also the case for NFP, bromobenzene, chlorobenzene,
dioxolane, quinoline, acetone and toluene. For the case of benzonitrile and
benzaldehyde the first solvation shell is comprised almost entirely of
molecules lying parallel to the graphene layer, with no second peak visible.
There is some experimental evidence for the confinement of solvent molecules
near the surface of graphene - Arunachalam et al. showed that the NMP
molecules near the surface of graphene show a reduction in rotational degrees
of freedom using rotating frame Overhauser effect spectroscopy – nuclear
magnetic resonance technique [58].

Figure 2: A snapshot of molecules in first solvation shell of benzonitrile
(left) and NMP (right) taken from the molecular dynamics simulation. Red
shading indicates molecules belonging to the first peak in the solvation shell
and green shading indicates the rest of the molecules.

To visualize this orientational ordering, Fig 2 shows a snapshot from the
molecular dynamics simulation showing only the first solvation shell of
benzonitrile and NMP. All molecules belonging to the very first peak are
shaded in red, and all others are shaded in green. As expected in the case of
benzonitrile, it can be seen that essentially all molecules in the first peak
are orientated parallel to the plane of graphene $(\alpha\approx 0^{\circ})$
with very few molecules adsorbed perpendicularly. In contrast, a small but
significant number of NMP molecules are adsorbed perpendicular to the sheet.
These molecules appear to fill all the gaps that exist in between the flat-
lying molecules.

To determine if any long-range lateral order is present, we calculated the
pair correlation functions for the geometric centers for the molecules in the
first solvation shell. These are shown in Fig 3. In all cases, there is a
local ordering around the center (0,0) due to the presence of an exclusion
zone around the molecule. This ordering then decays rapidly moving away from
the center. The strongest ordering is seen in case of cyclopentanone,
cyclohexanone and 1,3-dioxolane. Intermediate ordering behaviour is seen in
case of NFP, benzonitrile, NMP, quinoline, acetone and toluene, with no
ordering beyond the first shell seen in case of bromobenzene, chlorobenzene
and benzaldehyde. This is in agreement with the work of Terrones et al. who
found that while a monolayer of NMP molecules confined between two graphene
sheets exhibit long-range hexagonal ordering, this ordering is lost as the
thickness of the molecular layer increases [57].

Figure 3: Pair correlation function for the center of geometry of the solvent
molecules present in the first solvation shell.

In response to the formation of the first molecular layer, the rest of the
solvent then reorganizes to form further solvation shells, as seen in Fig 1.
As expected, the amplitude of the solvation shell peaks decay as the influence
of the solute wanes, finally approaching a constant value of 1, indicating
bulk solvent behaviour. The distance between graphene and the nearest solvent
atoms ($\mathrm{h}_{\mathrm{min}}$), the number of solvation shells
$(\mathrm{N}_{\mathrm{shells}})$, the depth to which the solvation shells
extend away from graphene, and the average distance between two consecutive
solvation shells are given in Table 2. As the distance between graphene and
the atom of the molecule nearest it is determined by the vdW interaction, it
is very similar across all solvents regardless of the molecular orientations,
varying from 1.99 Å for bromobenzene to 2.07 Å for benzaldehyde. The number of
solvation shells varies from a minimum of 3 in the cases of bromobenzene,
benzonitrile, chlorobenzene and benzaldehyde to a maximum of 6 for
cyclohexanone and 1,3-dioxolane, where in counting the total number of
solvation shells we have defined the last shell as the final peak that
deviates from bulk by at most 5%. The average distance between adjacent
solvation shells ranges from 4.25 Å for benzonitrile to 5.25 Å for quinoline.

Molecules | hmin | Nshells | Depth (Å) | Period (Å)
---|---|---|---|---
Cyclopentanone | 2.03 | 5 | 23.34 | 4.97
Cyclohexanone | 2.02 | 6 | 29.31 | 5.02
NFP | 2.0 | 5 | 24.67 | 5.14
Bromobenzene | 1.99 | 3 | 13.11 | 4.95
Benzonitrile | 2.03 | 3 | 12.16 | 4.25
NMP | 2.0 | 4 | 18.7 | 5.08
Chlorobenzene | 2.02 | 3 | 12.92 | 4.64
1,3-dioxolane | 2.05 | 6 | 25.94 | 4.41
Quinoline | 2.0 | 5 | 24.86 | 5.25
Benzaldehyde | 2.07 | 3 | 12.41 | 4.38
Acetone | 2.02 | 4 | 17.37 | 4.7
Toluene | 2.01 | 4 | 18.51 | 4.95
Table 2:  Data extracted from the pair distribution function $g(z)$ shown in
Fig. 1. hmin is the distance between graphene and the solvent atoms closest to
it. Nshells indicates the number of solvation shells present where the final
solvation shell is defined to that peak which deviates from the bulk value by
at most 5%. Depth is the distance of last solvation shell from graphene.
Period indicates the average periodicity of the solvation shells.

Despite these differences, it is clear from Fig 1 that the solvation structure
cannot be used to explain the observed phenomenon of solvent-dependent
graphene stabilization. For example, there is no appreciable difference
between the solvation structure of NMP and toluene, one of the best and worst
solvents for the exfoliation of graphite, respectively. However, the existence
of these solvation layers can explain the origin of the non-polar Hansen
solubility parameter attributed to graphene. Rather than to graphene alone,
the solubility parameters can be attributed to an effective solute comprised
of graphene and its first few solvation shells. Due to partial charges on the
atoms of the first solvation shell, the effective solute appears to have
surface charges, resulting in a designation of a non-zero polar solubility
parameter.

###  3.3 Molecular Diffusion

In the analysis up until now, we have assumed that the process of
stabilization is thermostatically driven, and kinetic effects such as
sedimentation play a minor role. This may not be necessarily true, and in
certain conditions may be the most important consideration. If the solvent-
dependence of graphene stability in a dispersion is determined by
reaggregation effects, then screening descriptors based on energetic
considerations alone will not be not appropriate. Long-term stability can only
be obtained if the dispersion is kinetically stabilized against reaggregation
when layer collision occurs. If the layers are far apart, this will be
influenced by the viscosity of the solvent [24, 59]. When the layers are close
together, the diffusion coefficient of molecules confined close to the
graphene sheet will determine the ease at they can be ejected to facilitate
reaggregation. Indeed, the addition of surfactants or polymers to graphene
dispersions to prevent reaggregation operate on this principle [60, 61].

Figure 4: The diffusion coefficient of solvent molecules in the first
solvation shell in the direction parallel to the surface of graphene as a
function of the experimental concentration of graphene $C_{G}$
($\mathrm{\mu}$g/ml$)$ in that solvent. Experimental data is taken from
Herndandez et al. [33].

To determine the extent to which solvent molecules can kinetically block
reaggregation, we calculate the diffusion coefficient in the direction
parallel to the graphene surface of the molecules in the first solvation
shell. The results are shown in Fig. 4 as a function of the experimental
graphene concentrations taken from Ref. [33]. The highest parallel diffusion
coefficient is found for acetone, at 0.35 Å2/ps. This is followed by toluene
(0.18 Å2/ps), dioxolane (0.16 Å2/ps) and chlorobenzene (0.16 Å2/ps). It is
notable that all of these are poor solvents for the exfoliation of graphite,
with the worst solvent corresponding to that with the largest parallel
diffusion coefficient. As the diffusion constant decreases, the experimental
concentration increases, with the best solvent (cyclopentanone) having a
parallel diffusion constant of 0.085 Å2/ps. However, in the very-low diffusion
regime the concentration does not follow this upward trend. As the molecular
hindrance to re-aggregation is higher in this regime, one would expect that
the exfoliated layers would remain in the dispersion resulting in an increased
monolayer concentration. Instead, the concentration decreases again.

While this requires further investigation, we hypothesis that the
concentration of graphene is now dependent on the ability of that solvent to
faciliate the exfoliation process itself, as opposed to the stabiliziation
process. When graphene is dispersed in a solvent with a low surface diffusion
coefficient, such as benzaldehyde, the reaggregation rate will be low and the
concentration of graphene layers will be maintained at its initial
concentration. In contrast, when the diffusion coefficient is high, the
graphene concentration will always be low due to a high rate of reaggregation.
Such solvents will result in a low concentration dispersion even if a large
number of monolayers are initially exfoliated due to efficient reaggregation.
As a consequence then, solvents with low surface diffusion coefficients, such
as quinoline and benzaldehyde, are good candidate for the optimization of the
initial separation process. If the initial concentration of graphene in these
solvents can be increased, the low diffusion rate in the first solvation shell
means that reaggregation will be hindered and the high concentration
maintained.

##  4 Conclusion

We find that there is a molecular-level structural and orientational ordering
of solvent molecules at the graphene-solvent interface which extends up to 30
Å away from the graphene layer. These solvation layers form in all cases and
behaves similarly irrespective of the solvent polarity. We speculate that the
non-polar Hansen parameter obtained by fitting the experimental data can be
attributed to the formation of the first solvation shell as the molecules in
the first solvation shell and the graphene form an effective solute. The
formation of these solvation layers is spontaneous as indicated by the free
energy and its decomposition in enthalpy and entropy contributions. This
independence of the solvation structure on the nature of the solvent molecule
suggests that non energetic contributions are responsible for the observed
solvent-dependence of graphene stabilization, and specifically the kinetics of
reaggregation. We propose that the diffusion coefficient of molecules in the
first solvation shell is an appropriate property to determine the probability
of reaggregation in a solvent. We note that in the low-diffusion limit, this
appears to break-down, at least for the solvents investigated here. This
origin of this requires further work, but could be related to the role of the
solvent in the initial separation step. We predict that the concentration of
graphene exfoliated in solvents with low diffusion coefficients parallel to
the graphene plane, such as benzaldehyde and quinoline, can be enhanced
dramatically by careful tuning of the experimental sonication parameters, as
any produced monolayers will be prevented from reaggregation by steric effects
due to the slowly diffusing solvent molecules on the monolayer surface.

This work was supported by a Science Foundation Ireland Starting Investigator
Research Grant (15/SIRG/3314). Computational resources were provided by the
supercomputer facilities at the Trinity Center for High Performance Computing
(TCHPC) and at the Irish Center for High-End Computing (ICHEC).

##  Appendix A Polar and azimuthal angle distribution functions

The position distribution function $g(z)$, the polar angle distribution
function $\theta(z)$ and the azimuthal angle distribution function $\phi(z)$
of the molecular dipoles close to graphene are shown in Fig5. The graphene
layer is located at 0 Å. The polar angle is measured with respect to the axis
aligned to the positive direction (towards the right). $g(z)$ is given
relative to the bulk solvent number density so that a value greater (less)
than 1 indicates an accumulation (depletion) of molecules relative to the bulk
solvent. $\theta(z)$ is the probability distribution function of the polar
angle of the molecular dipole, where 90∘ corresponds to the molecular dipole
in the plane parallel to graphene and 0∘ (or 180∘) corresponds to the
molecular dipole in the plane perpendicular to the plane of graphene.
$\phi(z)$ is the probability distribution function of the azimuthal angle of
the molecular dipole, indicating the orientation of the dipole in the plane
parallel to the plane of graphene, with 0∘ (or 180∘) corresponding to the
orientation along positive (negative) x-axis, 90∘ (-90∘) corresponding to the
orientation along the positive (negative) y-axis.

Figure 5: Pair distribution functions. (Top panel) $g(z)$ a histogram of the
center of mass of the solvent molecules normalised by $(\frac{N}{L}dz)$ where
$N$ is the number of molecules in the unit cell, $L$ is the length of the unit
cell and $dz$ is the bin width. The graphene monolayer is located at $z=0$ Å.
(This is the same as presented in Fig. 1). $\theta(z)$ is the probability
distribution function of the polar angle of the molecular dipole, where
$90^{\circ}$ corresponding to the molecular dipole in the plane parallel to
graphene and $0^{\circ}$ ($180^{\circ}$) corresponds to the molecular dipole
in the plane perpendicular to the plane of graphene. (Bottom panel) $\phi(z)$
is the probability distribution function of the azimuthal angle of the
molecular dipole, indicating the orientation of the dipole in the plane
parallel to the plane of graphene, with $0^{\circ}$ ($180^{\circ}$)
corresponding to the orientation along positive (negative) $x$-axis,
$90^{\circ}$ ($-90^{\circ}$) corresponding to the orientation along the
positive (negative) $y$-axis.

##  Appendix B Solvent Densities and Dipole Moments

Molecules | density (g/ml) | dipole moment (D)
---|---|---
Cyclopentanone | 0.95 | 3.28
Cyclohexanone | 0.95 | 2.9
NFP | 1.02 | -
Bromobenzene | 1.5 | 1.74
Benzonitrile | 1.01 | 3.2
NMP | 1.03 | 4.1
Chlorobenzene | 1.11 | 1.55
1,3-dioxolane | 1.06 | 1.19
Quinoline | 1.09 | 2.0
Benzaldehyde | 1.04 | 2.89
Acetone | 0.78 | 2.9
Toluene | 0.87 | 0.31
Table 3:  Experimental densities and dipole moments of the solvent molecules
considered here extracted from the on-line databases [62, 63].

## References

  * [1] Nicolosi V, Chhowalla M, Kanatzidis M G, Strano M S and Coleman J N 2013 Science 340
  * [2] Coleman J N 2012 Accounts of chemical research 46 14–22
  * [3] Lin S, Chui Y, Li Y and Lau S P 2017 FlatChem 2 15–37
  * [4] Raza A, Hassan J Z, Ikram M, Ali S, Farooq U, Khan Q and Maqbool M 2021 Advanced Materials Interfaces 8 2002205
  * [5] Secor E B, Ahn B Y, Gao T Z, Lewis J A and Hersam M C 2015 Advanced Materials 27 6683–6688
  * [6] Torrisi F and Coleman J N 2014 Nature nanotechnology 9 738–739
  * [7] Worsley R, Pimpolari L, McManus D, Ge N, Ionescu R, Wittkopf J A, Alieva A, Basso G, Macucci M, Iannaccone G et al. 2018 ACS nano 13 54–60
  * [8] Boland C S, Khan U, Ryan G, Barwich S, Charifou R, Harvey A, Backes C, Li Z, Ferreira M S, Möbius M E et al. 2016 Science 354 1257–1260
  * [9] Davesne A L, Lazar S, Bellayer S, Qin S, Grunlan J C, Bourbigot S and Jimenez M 2019 ACS Applied Nano Materials 2 5450–5459
  * [10] Li Z, Young R J, Backes C, Zhao W, Zhang X, Zhukov A A, Tillotson E, Conlan A P, Ding F, Haigh S J et al. 2020 ACS nano 14 10976–10985
  * [11] Backes C, Higgins T M, Kelly A, Boland C, Harvey A, Hanlon D and Coleman J N 2017 Chemistry of materials 29 243–255
  * [12] Hernandez Y, Nicolosi V, Lotya M, Blighe F M, Sun Z, De S, McGovern I, Holland B, Byrne M, Gun’Ko Y K et al. 2008 Nature nanotechnology 3 563–568
  * [13] Tung T T, Yoo J, Alotaibi F K, Nine M J, Karunagaran R, Krebsz M, Nguyen G T, Tran D N, Feller J F and Losic D 2016 ACS applied materials & interfaces 8 16521–16532
  * [14] Khan U, O’Neill A, Porwal H, May P, Nawaz K and Coleman J N 2012 Carbon 50 470–475
  * [15] Backes C, Campi D, Szydlowska B M, Synnatschke K, Ojala E, Rashvand F, Harvey A, Griffin A, Sofer Z, Marzari N et al. 2019 ACS nano 13 7050–7061
  * [16] Xiang K, Huang S, Song H, Bazhenov V, Bellucci V, Birnsteinova S, de Wijn R, Koliyadu J C, Koua F H, Round A et al. 2023 arXiv preprint arXiv:2305.08538
  * [17] Zhang X, Coleman A C, Katsonis N, Browne W R, Van Wees B J and Feringa B L 2010 Chemical Communications 46 7539–7541
  * [18] Liang Y T and Hersam M C 2010 Journal of the American Chemical Society 132 17661–17663
  * [19] Li J, Ye F, Vaziri S, Muhammed M, Lemme M C and Östling M 2012 Carbon 50 3113–3116
  * [20] Hernandez Y, Lotya M, Rickard D, Bergin S D and Coleman J N 2009 Langmuir 26 3208–3213
  * [21] Coleman J N, Lotya M, O’Neill A, Bergin S D, King P J, Khan U, Young K, Gaucher A, De S, Smith R J et al. 2011 Science 331 568–571
  * [22] ECHA 2020 Candidate list of substances of very high concern for authorisation https://echa.europa.eu/candidate-list-table
  * [23] Morton J A, Kaur A, Khavari M, Tyurnina A V, Priyadarshi A, Eskin D G, Mi J, Porfyrakis K, Prentice P and Tzanakis I 2023 Carbon 204 434–446
  * [24] Salavagione H J, Sherwood J, Budarin V, Ellis G, Clark J, Shuttleworth P et al. 2017 Green Chemistry 19 2550–2560
  * [25] Backes C, Paton K R, Hanlon D, Yuan S, Katsnelson M I, Houston J, Smith R J, McCloskey D, Donegan J F and Coleman J N 2016 Nanoscale 8 4311–4323
  * [26] Tao H, Zhang Y, Gao Y, Sun Z, Yan C and Texter J 2017 Physical Chemistry Chemical Physics 19 921–960
  * [27] Liang K Y and Yang W D 2018 AIP Advances 8 015018
  * [28] Shen J, He Y, Wu J, Gao C, Keyshar K, Zhang X, Yang Y, Ye M, Vajtai R, Lou J et al. 2015 Nano letters 15 5449–5454
  * [29] Goldie S J, Degiacomi M T, Jiang S, Clark S J, Erastova V and Coleman K S 2022 ACS nano 16 16109–16117
  * [30] Lin S, Shih C J, Strano M S and Blankschtein D 2011 Journal of the American Chemical Society 133 12810–12823
  * [31] Mukhopadhyay T K and Datta A 2017 The Journal of Physical Chemistry C 121 10210–10223
  * [32] Mukhopadhyay T K and Datta A 2017 the Journal of Physical Chemistry C 121 811–822
  * [33] Hernandez Y, Lotya M, Rickard D, Bergin S D and Coleman J N 2010 Langmuir 26 3208–3213
  * [34] Plimpton S 1995 Journal of computational physics 117 1–19
  * [35] Lazar P, Karlicky F, Jurecka P, Kocman M, Otyepkova E, Safarova K and Otyepka M 2013 Journal of the American Chemical Society 135 6372–6377
  * [36] Patil U and Caffrey N M 2018 The Journal of Chemical Physics 149 094702
  * [37] Cheng A and Steele W 1990 The Journal of chemical physics 92 3858–3866
  * [38] Hardy A, Dix J, Williams C D, Siperstein F R, Carbone P and Bock H 2018 ACS nano 12 1043–1049
  * [39] Fu C and Yang X 2013 Carbon 55 350–360
  * [40] Kamath G and Baker G A 2012 Physical Chemistry Chemical Physics 14 7929–7933
  * [41] Jorgensen W L and Tirado-Rives J 2005 Proceedings of the National Academy of Sciences 102 6665–6670
  * [42] Dodda L S, Vilseck J Z, Tirado-Rives J and Jorgensen W L 2017 The Journal of Physical Chemistry B 121 3864–3870
  * [43] Dodda L S, Cabeza de Vaca I, Tirado-Rives J and Jorgensen W L 2017 Nucleic acids research 45 W331–W336
  * [44] Lorentz H 1881 Annalen der physik 248 127–136
  * [45] Berthelot D 1898 Compt. Rendus 126 1703–1706
  * [46] Martínez L, Andrade R, Birgin E G and Martínez J M 2009 Journal of computational chemistry 30 2157–2164
  * [47] Ong S P, Richards W D, Jain A, Hautier G, Kocher M, Cholia S, Gunter D, Chevrier V L, Persson K A and Ceder G 2013 Computational Materials Science 68 314–319
  * [48] Duarte Ramos Matos G, Kyu D Y, Loeffler H H, Chodera J D, Shirts M R and Mobley D L 2017 Journal of Chemical & Engineering Data 62 1559–1569
  * [49] Mezei M 1987 The Journal of chemical physics 86 7084–7088
  * [50] Kirkwood J G 1935 The Journal of Chemical Physics 3 300–313
  * [51] Jorgensen W L and Ravimohan C 1985 The Journal of chemical physics 83 3050–3054
  * [52] Liu P, Harder E and Berne B 2004 The Journal of Physical Chemistry B 108 6595–6602
  * [53] Agosta L, Brandt E G and Lyubartsev A P 2017 The Journal of chemical physics 147 024704
  * [54] Oyer A J, Carrillo J M Y, Hire C C, Schniepp H C, Asandei A D, Dobrynin A V and Adamson D H 2012 Journal of the American Chemical Society 134 5018–5021
  * [55] Mukhopadhyay T K and Datta A 2019 JOURNAL OF THE INDIAN CHEMICAL SOCIETY 96 753–766
  * [56] Frolov A I, Arif R N, Kolar M, Romanova A O, Fedorov M V and Rozhin A G 2012 Chemical science 3 541–548
  * [57] Terrones J, Kiley P J and Elliott J A 2016 Scientific reports 6 1–11
  * [58] Arunachalam V and Vasudevan S 2018 The Journal of Physical Chemistry C 122 1881–1888
  * [59] Gilliam M S, Yousaf A, Guo Y, Li D O, Momenah A, Wang Q H and Green A A 2021 Langmuir 37 1194–1205
  * [60] Griffin A, Nisi K, Pepper J, Harvey A, Szydłowska B M, Coleman J N and Backes C 2020 Chemistry of Materials 32 2852–2862
  * [61] Hu C X, Shin Y, Read O and Casiraghi C 2021 Nanoscale 13 460–484
  * [62] Saturated liquid density URL http://ddbonline.ddbst.de/DIPPR105DensityCalculation/DIPPR105CalculationCGI.exe?component=Acetone
  * [63] Kim S, Chen J, Cheng T, Gindulyte A, He J, He S, Li Q, Shoemaker B A, Thiessen P A, Yu B et al. 2019 Nucleic acids research 47 D1102–D1109