\section{Introduction} \label{sect:Intro}

Many real-world datasets can be described from multiple ``viewpoints", such as pictures taken from different angles of the same object, different language expressions of the same semantic, texts and images on the same web page, etc. The representations from different perspectives can be treated as different views. The essence of multi-view learning (MVL) is to exploit the consensual and complementary information between different views \cite{Tang2018,Xu2013} to achieve better learning performance.

MVL approaches can be divided into three major categories \cite{Xue2017,Baltrusaitis2019}:
\begin{enumerate}
\item \emph{Co-training} \cite{Blum1998,Abney2002}, which exchanges discriminative information between two views by training the two models alternately.
\item \emph{Multi-kernel learning} \cite{Gehler2009,Gonen2011}, which maps data to different feature spaces with different kernels, and then combines those projected features from all spaces.
\item \emph{Subspace learning} \cite{Hotelling1936,Diethe2008,Xia2010}, which assumes all views are generated from a latent common space where shared information of all views can be exploited.
\end{enumerate}

Canonical correlation analysis (CCA), first proposed by Hotelling \cite{Hotelling1936} in 1936, is a typical subspace learning approach. Its main idea is to find pairs of projections for different views so that the correlations between them are maximized. One example is illustrated in Fig.~\ref{fig:CCA}. Since CCA takes the relationship between different feature sets into account, which is consistent with the idea of MVL, it is widely used in MVL \cite{Zhao2017,Sun2013}, including multi-view dimensionality reduction \cite{Foster2008}, multi-view clustering \cite{Kamalika2009,Blaschko2008}, multi-view regression \cite{Kakade2007}, and so on.

\begin{figure}[htpb] \centering
\includegraphics[width=.98\linewidth,clip]{CCA}
\caption{Multi-view data and CCA-based subspace learning.} \label{fig:CCA}
\end{figure}

Traditional CCA has the following limitations:
\begin{enumerate}
\item It cannot handle more than two views.
\item It can only calculate the linear correlation between two views, whereas in many real-world applications the true relationship between the views may be nonlinear.
\item In supervised classification, labels are available; however, CCA, as an unsupervised algorithm, completely ignores the labels, and hence wastes information.
\end{enumerate}
Many extensions have been proposed in the past few decades \cite{Akaho2006,Sun2007,Wang2015DCCAE,Sun2008,Andrew2013,Benton2016,Elmadany2016,Luo2015,Nielsen2002} to accommodate these limitations. Some representative ones are briefly introduced next.

In 1961, Horst \cite{Horst1961} first proposed generalized canonical correlation analysis (GCCA) to estimate the pairwise correlations of multiple views. He provided two formulations: the sum of correlation (SUMCOR), and the maximum variance (MAXVAR). Carroll \cite{Carroll1968} in 1968 proposed to find a shared latent correlated space, which was shown to be identical to MAXVAR GCCA. In 1971, Kettenring \cite{Kettenring1971} added three new formulations to GCCA, which maximizes the sum of the squared correlation (SSQCOR), minimizes the smallest eigenvalue (MINVAR), and minimizes the determinant of the correlation matrix (GENVAR), respectively. Two decades later, Nielsen \cite{Nielsen2002} summarized four constraints for these five formulations, forming a total of 20 combinations. In 2007, Via \cite{Via2007} proposed least squares based CCA (LS-CCA), and showed that it is essentially identical to MAXVAR. Luo \emph{et al.} \cite{Luo2015} proposed in 2015 that the correlation of multiple views can be directly maximized by analyzing the high-order covariance tensor, namely tensor canonical correlation analysis (TCCA). In 2017, Benton \emph{et al.} \cite{Benton2016} proposed DNN-based deep generalized canonical correlation analysis (DGCCA), which was a nonlinear extension of MAXVAR GCCA.

Generally, there are three approaches to deal with complex nonlinear relationship between two views \cite{Sun2007}. The most common one is to project data onto a higher dimensional space using the ``kernel trick", e.g., kernel canonical correlation analysis (KCCA) \cite{Akaho2006}. However, global kernelization suffers from high computational complexity, and it is not easy to choose an optimal kernel function. The second approach tries to preserve the locality of data. Inspired by the graph model, Sun and Chen  \cite{Sun2007} proposed locality preserving canonical correlation analysis (LPCCA), which aimed to reduce the global nonlinear dimensionality while preserving the local linear structure of data. Since the distances between neighbors need to be calculated, it is time-consuming when the sample size is large. The third approach is based on deep neural networks (DNN), which can give a very complex mapping between data. Andrew \emph{et al.} \cite{Andrew2013} first proposed deep canonical correlation analysis (DCCA) in 2013. Inspired by the autoencoder, Wang \emph{et al.} \cite{Wang2015DCCAE} proposed deep canonically correlated autoencoders (DCCAE) in 2015. However, DNN models have poor interpretability, and require a large amount of data to fit.

The label information is critical for classification problems. In order to make full use of the discriminant information, Sun \emph{et al.} \cite{Sun2008} proposed discriminant canonical correlation analysis (DisCCA) in 2008, by taking the inter-class and intra-class similarities of different views into consideration. Similar to DisCCA, Sun \emph{et al.} \cite{Sun2016} further proposed multi-view linear discriminant analysis (MLDA), which combined CCA and linear discriminant analysis (LDA) \cite{Fisher1936}. Elmadany \emph{et al.} \cite{Elmadany2016} integrated neural networks into DisCCA to obtain a nonlinear supervised model. Benton \emph{et al.} \cite{Benton2016} made use of the discriminant information by treating the one-hot encoding matrix of the labels as an additional view.

This paper provides an overview of many representative CCA approaches. A summary of their main characteristics is shown in Table~\ref{tab:comparation}.

\begin{table}[htpb] \centering
\caption{Comparison of some typical CCA-based approaches.} \label{tab:comparation}
\begin{tabular}{l| c c c c} \toprule[0.08em]
\textbf{Approaches}&\textbf{Multiple views}&\textbf{Supervised}&\textbf{Nonlinear}&\textbf{DNN} \\ \hline
    CCA         &  &  &   &    \\
    sCCA  &  &  &   &    \\
    KCCA        &  &  & \checkmark  &    \\
    RCCA        &  &  & \checkmark  &    \\
    LPCCA       &  &  & \checkmark  &    \\
    DisDCCA    &  & \checkmark &   &    \\
    MLDA        &  & \checkmark &   &    \\
    MULDA       &  & \checkmark &   &    \\
    DCCA        &  &  &  \checkmark & \checkmark   \\
    DisDCCA    &  & \checkmark & \checkmark  &  \checkmark  \\
    DCCAE       &  &  & \checkmark  & \checkmark   \\
    DisDCCAE   &  & \checkmark & \checkmark  &   \checkmark   \\
    VCCA        &  &  & \checkmark  &  \checkmark   \\
    GCCA        & \checkmark &  &   &    \\
    MCCA        & \checkmark &  &   &    \\
    LS-CCA      & \checkmark &  &   &    \\
    TCCA        & \checkmark &  &   &    \\
    DGCCA       & \checkmark & \checkmark &  \checkmark &   \checkmark   \\
\bottomrule
\end{tabular}
\end{table}

The remainder of this paper is organized as follows: Sections~\ref{sec:TwoViewCCA} and \ref{sec:MultiViewCCA}  introduce CCA-based MVL approaches for two views and more than two views, respectively. Section~\ref{sec:application} provides their applications. Section~\ref{sec:Conclusion} draws conclusions.

\section{CCA-Based MVL Approaches for Two Views} \label{sec:TwoViewCCA}

This section first reviews the tradition CCA, and then introduces some representative nonlinear, sparse, and/or supervised extensions. Table~\ref{tab:notation} summarizes the notations used in this paper.

\begin{table}[htpb] \centering
\caption{List of notations.} \label{tab:notation}
\begin{tabular}{l l}
\toprule[0.08em]
\textbf{Notation}&\textbf{Description} \\ \hline
    $X$, $Y$ & Data matrices from two different views $x$ and $y$ \\
    $\left(\mathbf{x}_{i},\mathbf{y}_{i}\right)$ & The $i$-th paired instances from Views $x$ and $y$ \\
    $N$ &	Number of instances in each view\\
    $c$ &   Number of classes \\
    $n_{i}$& Number of instances in Class $i$ \\
    $J$ &   Number of views\\
    $K$ &   Number of canonical vectors \\
    $d_x$ & Feature dimensionality of View $x$ \\
    $\mathbf{w}_x$ &	A canonical vector of View $x$ \\
    $W_x$ &	Canonical matrix of View $x$ \\
    $\Sigma_{xy}$ & Covariance matrix of View $x$ and View $y$\\
    $\hat{\Sigma}_{xy}$ & Regularized covariance matrix of View $x$ and View $y$\\
    $r_x$ & Regularization coefficients of View $x$ \\
    $m$ & Dimension of low rank approximation  \\
    $\cdot^{T}$ & Transpose of a matrix or vector \\
    $I$ & The identity matrix\\
    $\|\cdot\|_{2}$ & 2-norm of a vector \\
    $\|\cdot\|_{F}$ & The Frobenius-norm of a matrix\\
\bottomrule
\end{tabular}
\end{table}

\subsection{Canonical Correlation Analysis (CCA) }\label{subsection:CCA}

Let $X=\left[\mathbf{x}_{1}, \mathbf{x}_{2}, \dots, \mathbf{x}_{N}\right] \in \mathbb{R}^{d_{x} \times N}$ and $Y =\left[\mathbf{y}_{1}, \mathbf{y}_{2}, \dots, \mathbf{y}_{N}\right]\in \mathbb{R}^{d_{y} \times N}$ be two mean-zero data matrices with $N$ instances and $d_x$ and $d_y$ features, respectively. CCA aims to find $K$ pairs of linear projections ${W}_{x}=\left[\mathbf{w}_{x,1}, \mathbf{w}_{x,2}, \dots, \mathbf{w}_{x,K}\right]\in \mathbb{R}^{d_{x}\times K}$ and ${W}_{y}=\left[\mathbf{w}_{y,1}, \mathbf{w}_{y,2}, \dots, \mathbf{w}_{y,K}\right]\in \mathbb{R}^{d_{x}\times K}$, called canonical vectors, so that the correlations between ${W}_{x}^{T}X$ and ${W}_{y}^{T}Y$ are maximized.

Take a canonical vector $\mathbf{w}_{x}\in \mathbb{R}^{d_{x}\times 1}$ for $X$ and a canonical vector $\mathbf{w}_{y}\in \mathbb{R}^{d_{y}\times 1}$ for $Y$ for example. CCA maximizes the correlation coefficient $\rho$ between $X^{T}\mathbf{w}_{x}$ and $Y^{T}\mathbf{w}_{y}$ \cite{Hardoon2004}, i.e.,
\begin{align}
\rho\left(X^{T}\mathbf{w}_{x}, Y^{T}\mathbf{w}_{y}\right)=\frac{\mathbf{w}_{x}^{T}X Y^{T} \mathbf{w}_{y}}{\sqrt{\left(\mathbf{w}_{x}^{T}XX^{T} \mathbf{w}_{x}\right)\left(\mathbf{w}_{y}^{T} Y Y^{T}\mathbf{w}_{y}\right)}}. \label{equ:CCA}
\end{align}

Since (\ref{equ:CCA}) is invariant to the scaling of $\mathbf{w}_{x}$ and $\mathbf{w}_{y}$, it can be transformed into the following constrained form:
\begin{align}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}}\ & \mathbf{w}_{x}^{T} XY^{T} \mathbf{w}_{y}\\ \nonumber
s.t.\ & \mathbf{w}_{x}^{T}X X^{T} \mathbf{w}_{x}=1, \mathbf{w}_{y}^{T}YY^{T} \mathbf{w}_{y}=1.
\end{align}

When the feature dimensionality is high, especially when $d_x>N$ (or $d_y>N$), the covariance matrix $XX^T$ (or $YY^T$) is singular, and hence the optimization problem is under-determined. Regularizations can be added to the covariance matrices to remedy this problem \cite{Bickel2008,Tijl2003,Warton2008}, by introducing:
\begin{align}
\hat{\Sigma}_{xx}&=\frac{1}{N} X X^{T}+ r_{x} I,\\
\hat{\Sigma}_{yy}&=\frac{1}{N} Y Y^{T}+r_{y} I,
\end{align}
where $r_{x}$ and $r_{y}$ are non-negative regularization coefficients.

There are two approaches for computing $W_x$ and $W_y$ directly. The first is to solve the following generalized eigenvalue decomposition problem \cite{Hardoon2004}:
\begin{align}
\left[ \begin{array}{cc}{\mathbf{0}} & {{\Sigma}_{xy}} \\ {{\Sigma}_{yx}} & {\mathbf{0}}\end{array}\right] \left[ \begin{array}{c}{\mathbf{w}_{\mathbf{x}}} \\ {\mathbf{w}_{\mathbf{y}}}\end{array}\right]=\lambda \left[ \begin{array}{cc}{\hat{\Sigma}_{xx}} & {\mathbf{0}} \\ {\mathbf{0}} & {\hat{\Sigma}_{yy}}\end{array}\right] \left[ \begin{array}{c}{\mathbf{w}_{\mathbf{x}}} \\ {\mathbf{w}_{\mathbf{y}}}\end{array}\right],
\end{align}
where
\begin{align}
{\Sigma}_{xy}&=\frac{1}{N} X Y^{T},\\
{\Sigma}_{yx}&=\frac{1}{N} Y X^{T}.
\end{align}
$\{[\mathbf{w}_{x,k};\mathbf{w}_{y,k}]\}_{k=1}^K$ are then the $K$ leading generalized eigenvectors. The correlation $\rho\left(\mathbf{w}_{x,k}^{T}X, \mathbf{w}_{y,k}^{T}Y\right)$ is equal to the $k$-th leading  generalized eigenvalue.

The second solution \cite{Andrew2013} performs singular value decomposition (SVD) on matrix $T=\hat{\Sigma}_{xx}^{-1 / 2} {\Sigma}_{xy} \hat{\Sigma}_{yy}^{-1 / 2}$. Let $\tilde{W}_{x}$ and $\tilde{W}_{y}$ be the $K$ leading left and right singular vectors of $T$. Then, the canonical matrices are $W_{x}= \hat{\Sigma}_{xx}^{-1 / 2} \tilde{W}_{x}$ and $W_{y}= \hat{\Sigma}_{yy}^{-1 / 2} \tilde{W}_{y}$. The correlation $\rho\left(\mathbf{w}_{x,k}^{T}X, \mathbf{w}_{y,k}^{T}Y\right)$ is equal to the $k$-th leading  singular value of $T$.

Once $W_x$ and $W_y$ are obtained, the projected new features, called canonical variables, are computed by
$Z_{x}=W_{x}^{T} X$ and $Z_{y}=W_{y}^{T} Y$.

\subsection{Sparse CCA (sCCA)}

Canonical matrices $W_{x}$ and $W_{y}$ calculated by CCA is dense. Based on penalized matrix decomposition (PMD) \cite{Witten2009}, Witten \emph{et al.}  \cite{Witten2009Extensions} proposed penalized CCA, called sCCA in this paper, to obtain spare canonical vectors. sCCA also remedies the problem that canonical vectors are not unique when $N$ is smaller than $d_{x}$ and/or $d_{y}$.

PMD \cite{Eckart1936} performs rank-$m$ approximation of an arbitrary matrix $M \in \mathbb{R}^{d \times N}$ by ${\hat{M}}=\sum_{i=1}^{m} d_{i} \mathbf{u}_{i} \mathbf{v}_{i}^{T}$, using SVD, where $d_{i}$ is the $i$-th leading singular value, and $\mathbf{u}_{i}\in \mathbb{R}^{d\times 1}$ and $\mathbf{v}_{i}\in \mathbb{R}^{N\times 1}$ are the corresponding left and right singular vectors, respectively. The sparsity of $\hat{M}$ can be achieved by imposing LASSO-constraints on $\mathbf{u}_{i}$ and $\mathbf{v}_{i}$. When $m=1$, minimizing $\|M-\hat{M}\|_{F}^{2}$ is equivalent to maximizing the following objective function:
\begin{align}
\max_ {\mathbf{u}, \mathbf{v}}\ & \mathbf{u}^{T} M \mathbf{v} \label{eq:PMD}\\
s.t.\ &\|\mathbf{u}\|_{2}^{2}=1,\ \|\mathbf{v}\|_{2}^{2}=1,\ p_{1}(\mathbf{u}) \leqslant c_{1},\ p_{2}(\mathbf{v}) \leqslant c_{2}, \nonumber
\end{align}
where $p_{1}(\mathbf{u})=\sum_{i=1}^d\left|u_{i}\right|$ and $p_2(\mathbf{v})=\sum_{i=1}^N\left|v_{i}\right|$ represent the LASSO penalties, and parameters $c_{1}$ and $c_{2}$ control the degree of sparsity.

sCCA can be solved by replacing $M$ in (\ref{eq:PMD}) with the cross covariance matrix $XY^T$, i.e.,
\begin{align}
\max_ {\mathbf{w}_{x}, \mathbf{w}_{y}}\ & \mathbf{w}_{x}^{T} X Y^T \mathbf{w}_{y}\\ \nonumber
s.t.\ & \|\mathbf{w}_{x}\|_{2}^{2}=1, \|\mathbf{w}_{y}\|_{2}^{2}=1, p_{1}(\mathbf{w}_{x}) \leqslant c_{1}, p_{2}(\mathbf{w}_{y}) \leqslant c_{2}.
\end{align}

$K$ pairs of canonical vectors, $\{\left(\mathbf{w}_{x,k},\mathbf{w}_{x,k}\right)\}_{k = 1}^K$, can be obtained by solving a multi-factor PMD \cite{Witten2009}.

\subsection{Kernel CCA (KCCA)}

The traditional CCA cannot be applied when the correlation between different views is nonlinear. Kernel CCA (KCCA) uses (nonlinear) kernels to project data onto a higher dimensional space for correlation analysis.

Let the projections be $\phi_{x}$ and $\phi_{y}$, and the projected views in the high-dimensional space be $\Phi_{x}=\left[\phi_{x}\left(\mathbf{x}_{1}\right), \phi_{x}\left(\mathbf{x}_{2}\right), \ldots, \phi_{x}\left(\mathbf{x}_{N}\right)\right]$ and $\Phi_{\mathrm{y}}=\left[\phi_{y}\left(\mathbf{y}_{1}\right), \phi_{y}\left(\mathbf{y}_{2}\right), \ldots, \phi_{y}\left(\mathbf{y}_{N}\right)\right]$, respectively. The objective function of KCCA is:
\begin{align}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}}\ & \mathbf{w}_{x}^{T} \Phi_{x}\Phi_{y}^{T} \mathbf{w}_{y} \label{eq:KCCA} \\
s.t.\ &\mathbf{w}_{x}^{T} \Phi_{x} \Phi_{x}^{T} \mathbf{w}_{x}=1, \mathbf{w}_{y}^{T} \Phi_{y} \Phi_{y}^{T} \mathbf{w}_{y}=1.\nonumber
\end{align}

Expressing $\mathbf{w}_{x}$ and $\mathbf{w}_{y}$ as linear combinations of the columns of ${\Phi}_{x}$ and ${\Phi}_{y}$, respectively \cite{Akaho2006}:
\begin{align}
\mathbf{w}_{x}&=\sum_{i=1}^{N} a^{i} \phi_{x}\left(\mathbf{x}_{i}\right)={\Phi}_{x} \mathbf{a},\\
\mathbf{w}_{y}&=\sum_{i=1}^{N} b^{i} \phi_{y}\left(\mathbf{y}_{i}\right)={\Phi}_{y} \mathbf{b},
\end{align}
where $\mathbf{a}=\left[a^{1},\cdots,a^{N}\right]^{T}$ and $\mathbf{b}=\left[b^{1},\cdots,b^{N}\right]^{T}$ are linear coefficients. (\ref{eq:KCCA}) can then be rewritten as:
\begin{align}
&\max _{\mathbf{a}, \mathbf{b}}\ \mathbf{a}^{T} {\Phi}_{x}^{T} {\Phi}_{x} {\Phi}_{y}^{T} {\Phi}_{y} \mathbf{b} \label{eq:KCCA2}\\
&s.t.\ \mathbf{a}^{T} {\Phi}_{x}^{T} {\Phi}_{x} {\Phi}_{x}^{T} {\Phi}_{x} \mathbf{a}=1,\
\mathbf{b}^{T} {\Phi}_{y}^{T} {\Phi}_{y} {\Phi}_{y}^{T} {\Phi}_{y} \mathbf{b}=1.\nonumber
\end{align}

Let $K_{x}={\Phi}_{x}^{T} {\Phi}_{x}$ and $K_y={\Phi}_{y}^{T} {\Phi}_{y}$ be kernel matrices, i.e., $({K}_{x})_{i j}={\kappa}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)$, where $\mathbf{\kappa}$ is a kernel function, such as a radial basis function (RBF). Then,  (\ref{eq:KCCA2}) can be further simplified to:
\begin{align}
\max _{\mathbf{a}, \mathbf{b}}\ & \mathbf{a}^{T} {K}_{x} {K}_{y} \mathbf{b}\\
s.t.\ &\mathbf{a}^{T} {K}_{x} {K}_{x} \mathbf{a}=1, \mathbf{b}^{T} {K}_{y} {K}_{y} \mathbf{b}=1.\nonumber
\end{align}

We can further add regularizations to the kernel matrices to make them numerically more stable \cite{Hardoon2004}, e.g., replace ${K}_{x}{K}_{x}$ with ${K}_{x}{K}_{x}+r_{x} {K}_{x}$, and ${K}_{y}{K}_{y}$ with ${K}_{y}{K}_{y}+r_{y} {K}_{y}$. The linear coefficient vectors $\mathbf{a}$ and $\mathbf{b}$ can then be solved by the following generalized eigen decomposition problem:
\begin{align}
&\left[\begin{array}{cc}
      \mathbf{0}   & K_{x}K_{y} \\
      K_{y}K_{x} &\mathbf{0}
      \end{array}\right]\left[\begin{array}{c}
                                \mathbf{a} \\
                                \mathbf{b}
                              \end{array}\right]\nonumber\\
=&      \lambda\left[\begin{array}{cc}
                     K_xK_x+r_xK_x & \mathbf{0} \\
                     \mathbf{0} & K_yK_y+r_yK_y
                   \end{array}\right]
                   \left[\begin{array}{c}
                                \mathbf{a} \\
                                \mathbf{b}
                              \end{array}\right]
\end{align}

\subsection{Randomized Nonlinear CCA (RCCA)}

KCCA has very high computational complexity. Lopez-Paz \emph{et al.} \cite{Lopez-Paz2014} handled the nonlinear characteristics of data by performing random nonlinear projection, which greatly reduced the computational difficulty, with little scarification of performance:
\begin{align}
\text{RCCA}\left(X,Y\right) :=\text{CCA}(f_x(X),f_y(Y))\approx \text{KCCA}(X,Y),
\end{align}
where $f_{x}(\cdot)$ represents a mapping function, whose parameters $Q_{x} = \left(\boldsymbol{q}_{x,1}, \boldsymbol{q}_{x,2}, \ldots, \boldsymbol{q}_{x,D}\right) \in \mathbb{R}^{d_{x} \times D}$ are randomly sampled from a given data-independent distribution $p(\boldsymbol{q})$, e.g., a Gaussian distribution. The $D$-dimensional nonlinear random features $Z_{x}$ can be obtained by:
\begin{align}
\mathbf{z}_{x,d} :=&\left[\cos \left(\boldsymbol{q}_{x,d}^{T} \mathbf{x}_{1}+b_{d}\right), \ldots, \cos \left(\boldsymbol{q}_{x,d}^{T} \mathbf{x}_{N}+b_{d}\right)\right]^{T} \in \mathbb{R}^{N\times 1}, \label{eq:zxd} \\
Z_{x}=&\left[\mathbf{z}_{x,1}, \ldots, \mathbf{z}_{x,D} \right]^{T}\in \mathbb{R}^{D\times N},\nonumber
\end{align}
where $\mathbf{b}=[b_{1}, b_{2}, \ldots, b_{D}]^T\in \mathbb{R}^{D\times 1}$ is a randomly generated bias term. The cosine mapping function in (\ref{eq:zxd}) can also be replaced by other nonlinear functions, e.g., sine. As in CCA, the computational complexity of RCCA also increases linearly with the sample size.

\subsection{Locality Preserving CCA (LPCCA)}

In addition to the computational complexity, the global kernelization strategy of KCCA ignores the local linear structure of complex data. Locality preserving CCA (LPCCA) \cite{Sun2007} is another nonlinear CCA extension, which preserves the local linear structure of the data while performing global nonlinear dimensionality reduction.

LPCCA assumes that the corresponding instances of different views should be as close as possible in the common latent space, so it can be expressed in the following equivalent form \cite{Hoegaerts2005}:
\begin{align}
\min _{\mathbf{w}_{x}, \mathbf{w}_{y}}& \sum_{i=1}^{N}\left\|\mathbf{w}_{x}^{T}\left(\mathbf{x}_{i}-\overline{\mathbf{x}}\right)
-\mathbf{w}_{y}^{T}\left(\mathbf{y}_{i}-\overline{\mathbf{y}}\right)\right\|_{2}^{2},\label{eq:wx0} \\
s.t.& \sum_{i=1}^{N}\left\|\mathbf{w}_{x}^{T}\left(\mathbf{x}_{i}-\overline{\mathbf{x}}\right)\right\|^{2}=1, \sum_{i=1}^{n}\left\|\mathbf{w}_{y}^{T}\left(\mathbf{y}_{i}-\overline{\mathbf{y}}\right)\right\|^{2}=1 \nonumber
\end{align}
where $\overline{\mathbf{x}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i}$ and $\overline{\mathbf{y}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{y}_i$ represent the mean vectors of ${X}$ and ${Y}$, respectively. After some algebraic manipulations, (\ref{eq:wx0}) can be re-expressed as:
\begin{align}
\max_{\mathbf{w}_{x}, \mathbf{w}_{y}} &\mathbf{w}_{x}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)^{T} \cdot \mathbf{w}_{y} \label{equ:CCA2}\\
s.t. \ &\mathbf{w}_{x}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)^{T} \cdot \mathbf{w}_{x}=1, \nonumber \\
&\mathbf{w}_{y}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)^{T} \cdot \mathbf{w}_{y}=1. \nonumber
\end{align}

Let $\text{Nei} \left(\mathbf{x}_i\right)$ be the neighbor set of $\mathbf{x}_{i}$. Define a similarity matrix ${S}_{x}$, whose $ij$-th element is:
\begin{align}
S_{x,i j}=\left\{\begin{array}{ll}
\exp \left(-\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|_{2}^{2} / t_{x}\right),& \mathbf{x}_{j} \in \text{Nei}\left(\mathbf{x}_{i}\right) \\
0, &\text { otherwise }\end{array}\right. ,
\end{align}
where $t_{x}=\sum_{i=1}^{N}\sum_{j=1}^{N} \frac{ 2\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|_{2}^{2}}{N(N-1)}$ represents the mean squared distance between all instances. The similarity matrix ${S}_{y}$ of ${Y}$ can be computed in a similar manner.

Substituting $S_x$ and $S_y$ into (\ref{equ:CCA2}), the global correlation is decomposed into many local linear correlations between the neighboring instances. The objective of LPCCA is:
\begin{align}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}} &\mathbf{w}_{x}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}S_{x,i j}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)S_{y,ij}\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)^{T} \cdot \mathbf{w}_{y}\\  \nonumber
s.t. \ &\mathbf{w}_{x}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}S_{x,i j}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)S_{x,ij}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)^{T} \cdot \mathbf{w}_{x}=1, \\ \nonumber
&\mathbf{w}_{y}^{T} \cdot \sum_{i=1}^{N} \sum_{j=1}^{N}S_{y,i j}\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)S_{y,ij}\left(\mathbf{y}_{i}-\mathbf{y}_{j}\right)^{T} \cdot \mathbf{w}_{y}=1,
\end{align}
which can be transformed into a generalized eigen decomposition problem:
\begin{align}
&\left[
\begin{array}{cc}
  \mathbf{0} &XS_{xy}Y^T \\
  YS_{yx}X^T & \mathbf{0}
\end{array}\right]\left[\begin{array}{c}
                          \mathbf{w}_x \\
                          \mathbf{w}_y
                        \end{array}\right] \nonumber \\
=&\lambda\left[\begin{array}{cc}
                 XS_{xx}X^T & \mathbf{0} \\
                 \mathbf{0} & YS_{yy}Y^T
               \end{array}\right]
               \left[\begin{array}{c}
                          \mathbf{w}_x \\
                          \mathbf{w}_y
                        \end{array}\right],
\end{align}
where
\begin{align}
{S}_{xy}&={D}_{xy}-{S}_{x}\odot{S}_y,\\
{S}_{yx}&={D}_{yx}-{S}_{y}\odot{S}_x,\\
{S}_{xx}&={D}_{xx}-{S}_{x}\odot{S}_x,\\
{S}_{yy}&={D}_{yy}-{S}_{y}\odot{S}_y,
\end{align}
in which $\odot$ denotes the element-wise product operation, and $D_{xy}$ is a diagonal matrix, with $(D_{xy})_{ii}=\sum_{j=1}^{N}({S}_{x}\odot{S}_y)_{ij}$. $D_{yx}$, $D_{xx}$ and $D_{yy}$ are defined similarly.

\subsection{Discriminative CCA (DisCCA)}

Traditional CCA is unsupervised. In supervised classification, we have label information, which should be taken into consideration to help extract more discriminative features.

Discriminative CCA (DisCCA) is one such approach. It maximizes the within-class similarity and minimizes the between-class similarity. Rearrange $X$ and $Y$ according to the classes:
\begin{align}
&\hat{{X}}=\left[\mathbf{x}_{1}^{(1)}, \ldots, \mathbf{x}_{n_{1}}^{(1)}, \ldots \ldots, \mathbf{x}_{1}^{(c)}, \dots, \mathbf{x}_{n_{c}}^{(c)}\right],\label{equ:Xhat}\\
&\hat{{Y}}=\left[\mathbf{y}_{1}^{(1)}, \ldots, \mathbf{y}_{n_{1}}^{(1)}, \ldots \ldots, \mathbf{y}_{1}^{(c)}, \ldots, \mathbf{y}_{n_{c}}^{(c)}\right],
\end{align}
where $\mathbf{x}_{j}^{(i)}$ and $\mathbf{y}_{j}^{(i)}$ are the $j$-th instance of Class $i$ from the two views.  The objective function of DisCCA is:
\begin{align}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}}\ &\mathbf{w}_{x}^{T} {C}_{\mathrm{w}} \mathbf{\mathrm{w}}_{y}-\eta \cdot \mathbf{w}_{x}^{T} {C}_{\mathrm{b}} \mathbf{w}_{y} \label{eq:wx}\\
s.t.\ & \mathbf{w}_{x}^{T} \hat{X} \hat{X}^{T} \mathbf{w}_{x}=1, \mathbf{w}_{y}^{T} \hat{Y} \hat{Y}^{T} \mathbf{w}_{y}=1,\nonumber
\end{align}
where $\eta$ is a trade-off parameter, and the within-class similarity matrix ${C}_{\mathrm{w}}$ and between-class similarity matrix ${C}_{\mathrm{b}}$ are defined as:
\begin{align}
&{C}_{\mathrm{w}}=\sum_{i=1}^{c} \sum_{k=1}^{n_{i}} \sum_{l=1}^{n_{j}} \mathbf{x}_{k}^{(i)} \left(\mathbf{y}_{l}^{(j)}\right)^T=\hat{X} {A} \hat{Y}^{T},\\
&{C}_{\mathrm{b}}=\sum_{i=1}^{c} \sum_{j=1 \atop j \neq i}^{c} \sum_{k=1}^{n_{i}} \sum_{l=1}^{n_{j}} \mathbf{x}_{k}^{(i)} \left(\mathbf{y}_{l}^{(j)}\right)^ T=-\hat{X} {A} \hat{Y}^{T},
\end{align}
where $A$ is a block-diagonal matrix:
\begin{align}\label{equ:A}
{A} = \left[\begin{array}{ccccc}
             \mathbf{1}_{n_{1} \times n_{1}} & \quad & \quad & \quad & \quad \\
             \quad & \ddots & \quad & \quad & \quad \\
             \quad & \quad & \mathbf{1}_{n_{i} \times n_{i}} & \quad & \quad \\
             \quad & \quad & \quad & \ddots & \quad \\
             \quad & \quad & \quad & \quad & \mathbf{1}_{n_{c} \times n_{c}}
           \end{array}\right].
\end{align}

Since ${C}_{\mathrm{w}}=-{C}_{\mathrm{b}}$, the optimization is independent of $\eta$, and hence (\ref{eq:wx}) can be rewritten as:
\begin{align}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}}\ &\mathbf{w}_{x}^{T} \hat{X} {A} \hat{Y}^{T} \mathbf{w}_{y}\\   \nonumber
s.t.\ &\mathbf{w}_{x}^{T} \hat{X} \hat{X}^{T} \mathbf{w}_{x}=1, \mathbf{w}_{y}^{T} \hat{Y} \hat{Y}^{T} \mathbf{w}_{y}=1,
\end{align}
whose solution is similar to (\ref{equ:CCA}) and can also be expressed as a generalized eigen decomposition problem:
\begin{align}
\begin{bmatrix}{\mathbf{0}}&{\tilde{\Sigma}_{xy}} \\ {\tilde{\Sigma}_{yx}}&{\mathbf{0}}\end{bmatrix}
\begin{bmatrix}{\mathbf{w}_{\mathrm{x}}} \\ {\mathbf{w}_{\mathrm{y}}}\end{bmatrix}=\lambda
\begin{bmatrix}{\hat{\Sigma}_{xx}} & {\mathbf{0}} \\ {\mathbf{0}} & {\hat{\Sigma}_{yy}}\end{bmatrix}
\begin{bmatrix}{\mathbf{w}_{x}} \\ {\mathbf{w}_{y}}\end{bmatrix},
\end{align}
where $\tilde{\Sigma}_{xy}=\frac{1}{N}\hat{X}{A} \hat{Y}^{T}$ and $\tilde{\Sigma}_{yx}=\frac{1}{N}\hat{Y} {A} \hat{X}^{T}$.

\subsection{Multi-view Linear Discriminant Analysis (MLDA)}

LDA \cite{Fisher1936} is a supervised algorithm for a single view that minimizes the within-class variance and maximizes the between-class variance. Multi-view linear discriminant analysis (MLDA) \cite{Sun2016} combines LDA and CCA, which not only ensures the discriminative ability within a single view, but also maximizes the correlation between different views.

Define a between-class scatter matrix $S_{b}=\frac{1}{N} \hat{{X}} {W} \hat{X}^{T}$, where $\hat{X}$ is given in (\ref{equ:Xhat}), $W = diag({W}_{1},{W}_2,\dots, {W}_{c})$, with all elements in ${W}_{i}\in \mathbb{R}^{n_{i}\times n_{i}}$ equal $\frac{1}{n_i}$. Then, for a single view, the objective function of LDA is:
\begin{align}
\max _{\mathbf{w}}\ & \mathbf{w}^{T} {S}_{\mathrm{b},x} \mathbf{w}\\ \nonumber
s.t.\ & \mathbf{w}^{T}\hat{\Sigma}_{xx} \mathbf{w}=1.
\end{align}

Integrating with (\ref{equ:CCA}), the objective function of MLDA is:
\begin{align}\label{equ:MLDA}
\max _{\mathbf{w}_{x}, \mathbf{w}_{y}}\ & \mathbf{w}_{x}^{T} {S}_{\mathrm{b},x} \mathbf{w}_{x}+\mathbf{w}_{y}^{T} {S}_{\mathrm{b},y} \mathbf{w}_{y}+\eta \mathbf{w}_{x}^{T} \Sigma_{xy} \mathbf{w}_{y} \\  \nonumber
s.t.\ & \mathbf{w}_{x}^{T} \hat{\Sigma}_{xx} \mathbf{w}_{x}+\sigma\mathbf{w}_{y}^{T} \hat{\Sigma}_{yy} \mathbf{w}_{y}=1,
\end{align}
where $\eta$ is a trade-off parameter, and $\sigma=\frac{tr(\hat{\Sigma}_{xx})}{tr(\hat{\Sigma}_{yy})}$. The constraint in (\ref{equ:MLDA}) is in the form of summation to ensure a closed-form solution can be obtained. Using the Lagrangian multiplier, (\ref{equ:MLDA}) can also be solved by generalized eigen decomposition:
\begin{align}
\begin{bmatrix}{S_{\mathrm{b},x}} & {\eta {\Sigma}_{x y}} \\ {\eta {\Sigma}_{y x}} & {S_{\mathrm{b},y}}\end{bmatrix}
\begin{bmatrix}{\mathbf{w}_{\mathrm{x}}} \\ {\mathbf{w}_{\mathrm{y}}}\end{bmatrix}=\lambda
\begin{bmatrix}{\hat{\Sigma}_{xx}} & {\mathbf{0}} \\ {\mathbf{0}} & {\hat{\Sigma}_{yy}}\end{bmatrix}
\begin{bmatrix}{\mathbf{w}_{\mathbf{x}}} \\ {\mathbf{w}_{\mathbf{y}}}\end{bmatrix}.
\end{align}

\subsection{Multi-view Uncorrelated Linear Discriminant Analysis (MULDA)}

Canonical vectors $\{\left(\mathbf{w}_{x,k}, \mathbf{w}_{y,k}\right)\}_{k=1}^K$, calculated by MLDA, are correlated, which means the projected canonical variables contain redundant information. Multi-view uncorrelated linear discriminant analysis (MULDA) \cite{Sun2016} aims to eliminate the information redundancy of the canonical variables using uncorrelated linear discriminant analysis (ULDA) \cite{Jin2001}.

ULDA imposes orthogonal constraints to LDA, and is solved recursively. Each time a projection vector $\mathbf{w}_k$ is calculated to be orthogonal to the $k-1$ projection vectors that have been obtained:
\begin{align}
\max _{\mathbf{w}} \ & \mathbf{w}^{T} {S}_{\mathrm{b},x} \mathbf{w}\\   \nonumber
s.t.\ & \mathbf{w}_{k}^{T} \hat{\Sigma}_{xx} \mathbf{w}_{i}=0,\quad i=1,2, \ldots, k-1.
\end{align}

Integrating ULDA with MLDA, the objective function of MULDA can be written as:
\begin{align}
\max _{\mathbf{w}_{x,k}, \mathbf{w}_{y,k}}\ & \mathbf{w}_{x,k}^{T} {S}_{\mathrm{b},x} \mathbf{w}_{x,k}+\mathbf{w}_{y,k}^{{T}} {S}_{\mathrm{b},y} \mathbf{w}_{y,k}+\eta \mathbf{w}_{x,k}^{T} \Sigma_{x y} \mathbf{w}_{y,k} \\ \nonumber
s.t.\ & \mathbf{w}_{x,k}^{T} \hat{\Sigma}_{xx} \mathbf{w}_{x,k}+\sigma\mathbf{w}_{y,k}^{T} \hat{\Sigma}_{yy} \mathbf{w}_{y,k}=1, \\ \nonumber
&\mathbf{w}_{x,k}^{T} \hat{\Sigma}_{xx} \mathbf{w}_{x,i}=0, \mathbf{w}_{y,k}^{T} \hat{\Sigma}_{yy} \mathbf{w}_{y,i}=0, \\  \nonumber
& i=1,\ldots,k-1.
\end{align}
The $k$-th pair of canonical vectors, $\left(\mathbf{w}_{x,k},\mathbf{w}_{y,k}\right)$, is the leading eigenvector of
\begin{align}
\begin{bmatrix} {P_{x}}& {\mathbf{0}} \\ {\mathbf{0}} & {P_{y}} \end{bmatrix}
\begin{bmatrix} {S_{\mathrm{b},x}} & {\eta {\Sigma}_{x y}} \\ {\eta {\Sigma}_{y x}} & {S_{\mathrm{b},y}} \end{bmatrix}
\begin{bmatrix} {\mathbf{w}_{x,k}} \\ {\mathbf{w}_{y,k}}\end{bmatrix}=\lambda
\begin{bmatrix} {\hat{\Sigma}_{xx}} & {\mathbf{0}} \\ {\mathbf{0}} & {\hat{\Sigma}_{yy}}\end{bmatrix}
\begin{bmatrix} {\mathbf{w}_{x,k}} \\ {\mathbf{w}_{y,k}}\end{bmatrix},
\end{align}
where
\begin{align}
&{P}_{x}={I}-\hat{\Sigma}_{xx} {D}_{x}^{T}\left({D}_{x} \hat{\Sigma}_{xx} {D}_{x}^{T}\right)^{-1} {D}_{x}, \\
&{D}_{x}=\left[\mathbf{w}_{x,1}, \mathbf{w}_{x,2}, \cdots, \mathbf{w}_{x,k-1}\right]^{T}.
\end{align}
${P}_{y}$ and ${D}_{y}$ are defined similarly.

Unlike DisCCA, MLDA and MULDA only focus on the within-view class scatters, without considering the between-view class scatter. Sun \emph{et al.} \cite{Sun2016} proposed that CCA in MLDA and MULDA can be replaced by DisCCA, to also consider the between-view class scatter. Moreover, MLDA and MULDA can also be extended using kernels to consider the nonlinear relationship between different views.

\subsection{Deep CCA (DCCA) and Discriminative DCCA (DisDCCA)}

Deep canonical correlation analysis (DCCA) \cite{Andrew2013}, shown in Fig.~\ref{fig:DCCA}, first extracts nonlinear features through deep neural networks (DNNs), and then uses linear CCA to calculate the canonical matrices.

Let $\mathbf{f}_x$ and $\mathbf{f}_y$ be two DNNs, and ${H}_{x}=\mathbf{f}_{x}(X)$ and ${H}_{y}=\mathbf{f}_{y}(Y)$ be their outputs, respectively. We already know that the total correlation of the $K$ canonical variables equals the sum of the first $K$ singular values of matrix ${T}=\hat{\Sigma}_{x x}^{-1/2}{\Sigma}_{x y}\hat{\Sigma}_{y y}^{-1/2}$ from \ref{subsection:CCA}. Define
\begin{align} \label{equ:T}
\hat{T}=&\left(\frac{1}{N}H_{x}H_{x}^{T}+r_{x}I\right)^{-\frac{1}{2}}\nonumber\\
& \left(\frac{1}{N}H_{x}H_{y}^{T}\right)\left(\frac{1}{N}H_{y}H_{y}^{T}+r_{y}I\right)^{-\frac{1}{2}}.
\end{align}
The objective function of DCCA is:
\begin{align}
\max_{\mathbf{f}_{x}, \mathbf{f}_{y}, \mathbf{w}_{x}, \mathbf{w_{y}}}&\sum_{k=1}^{K} \sigma_{k}(\hat{T})\label{equ:DCCA} \\
s.t.\ &\mathbf{w}_{x}{\left(\frac{1}{N}H_{x}H_{x}^{T}+r_{x}I\right)}\mathbf{w}_x^{T}=1,\nonumber \\
& \mathbf{w}_{y}{\left(\frac{1}{N}H_{y}H_{y}^{T}+r_{y}I\right)}\mathbf{w}_y^{T}=1,\nonumber
\end{align}
where $\sigma_{k}(\hat{T})$ denotes the $k$-th largest singular value of $\hat{T}$.

Andrew \emph{et al.} \cite{Andrew2013} employed a full-batch optimization algorithm (L-BFGS) to optimize (\ref{equ:DCCA}), because the covariance matrix of the entire training set needs to be computed. If the dataset is divided into small batches, then the covariance matrix computed from each batch may not be accurate. However, full-batch optimization is both memory-hungry and time-consuming, especially when the dataset is large. Wang \emph{et al.} \cite{Wang2015Unsupervised} proposed that with a large batch size, the instances in each batch can be sufficient for estimating the covariance matrix, and the efficient mini-batch stochastic gradient descent approach can be used.

\begin{figure}[htbp] \centering
\includegraphics[height=.7\linewidth,clip]{DCCA}
\caption{Schematic diagram of DCCA. Two deep neural networks $\mathbf{f}_{x}$ and $\mathbf{f}_{y}$ first extract nonlinear features from $X$ and $Y$, respectively, and then a linear CCA is applied for correlation analysis. } \label{fig:DCCA}
\end{figure}

In experiments, the network outputs on the training set were first obtained, and then a linear CCA was used to obtain the canonical matrices ${W}_{x}$ and ${W}_y$. On the test set, ${W}_{x}^{T}\mathbf{f}_{x}(\cdot)$ and ${W}_{y}^{T}\mathbf{f}_{y}(\cdot)$ were the final inputs to a classifier.

Similar to DCCA, discriminative deep canonical correlation analsis (DisDCCA) \cite{Elmadany2016} is a DNN-based extension of discriminative CCA. There are two major differences in their implementations. First, when calculating the loss of each batch, the instances are first rearranged according to their classes, and then $\frac{1}{N}H_{x}H_y^T$ in (\ref{equ:T}) is replaced by $\frac{1}{N} {H}_{x}{A} {H}_{y}^{T}$, where $A$ is defined in (\ref{equ:A}). Second, after the model is trained, DisCCA is used to obtain the canonical matrices ${W}_{x}$ and ${W}_{y}$.

\subsection{Deep Canonically Correlated Autoencoders (DCCAE) and Discriminative DCCAE (DisDCCAE)}

Deep canonically correlated autoencoders (DCCAE) \cite{Wang2015DCCAE} improves DCCA, by using autoencoders to make sure the information captured by $\mathbf{f}_x$ and $\mathbf{f}_y$ can also accurately reconstruct the original $X$ and $Y$. Fig.~\ref{fig:DCCAE} shows the DCCAE model structure, where $\mathbf{g}_{x}$ and $\mathbf{g}_{y}$ represent the decoder networks for reconstructing $X$ and $Y$, respectively. Therefore, two reconstruction errors are incorporated into (\ref{equ:DCCA}):
\begin{align}
\max_{\mathbf{f}_{x}, \mathbf{f}_{y},\mathbf{g}_{x}, \mathbf{g}_{y}, \mathbf{w}_{x}, \mathbf{w_{y}}}\ &\sum_{k=1}^{K} \sigma_{k}\left(\hat{T}\right)
-\frac{\lambda}{N}\left(\left\|{X}-\mathbf{g}_{x}\left(\mathbf{f}_{x}\left({X}\right)\right)\right\|^{2}_{F} \right. \nonumber\\  &
\left.+ \left\|{Y}-\mathbf{g}_{y}\left(\mathbf{f}_{y}\left({Y}\right)\right)\right\|^{2}_{F}\right) \label{equ:DCCAE}\\
s.t.\ &\mathbf{w}_{x}{(\frac{1}{N}H_{x}H_{x}^{T}+r_{x}I)}\mathbf{w}_x^{T}=1,\nonumber \\
& \mathbf{w}_{y}{(\frac{1}{N}H_{y}H_{y}^{T}+r_{y}I)}\mathbf{w}_y^{T}=1.\nonumber
\end{align}

\begin{figure}[htbp] \centering
\includegraphics[height=.7\linewidth,clip]{DCCAE}
\caption{Schematic diagram of DCCAE. Two deep neural networks $\mathbf{f}_{x}$ and $\mathbf{f}_{y}$ first extract nonlinear features from $X$ and $Y$, respectively, and then a linear CCA is applied for correlation analysis. $\mathbf{g}_{x}$ and $\mathbf{g}_{y}$ are two decoders for reconstructing View $X$ and View $Y$, respectively. }\label{fig:DCCAE}
\end{figure}

The test procedure is similar to DCCA, except that ${W}_{x}^{T}\mathbf{f}_{x}(\cdot)$ and ${W}_{y}^{T}\mathbf{f}_{y}(\cdot)$ are used as inputs to a classifier. The advantage of DCCAE is that the autoencoders can alleviate the over-fitting of the model.

Since the first term of (\ref{equ:DCCAE}) is identical to that in (\ref{equ:DCCA}), extending DCCAE to discriminative DCCAE (DisDCCAE) is similar to extending DCCA to DisDCCA. However, to our knowledge, no one has done that before. We included both DCCAE and DisDCCAE in our experiments.

\subsection{Variational CCA (VCCA) and VCCA-private}

Bach and Jordan \cite{Bach2005} explained CCA from the perspective of a probabilistic latent variable model. It assumes that the instances ${\mathbf{x}}$ and ${\mathbf{y}}$ from two views are independently conditioned on the multivariate latent variable ${\mathbf{z}}\in \mathbb{R}^{d_{z}\times 1}$, and CCA aims to maximize the joint probability distribution of $\mathbf{x}$ and $\mathbf{y}$:
\begin{align} p(\mathbf{x}, \mathbf{y}, \mathbf{z}) &=p(\mathbf{z}) p(\mathbf{x} | \mathbf{z}) p(\mathbf{y} | \mathbf{z}), \\ p(\mathbf{x}, \mathbf{y}) &=\int p(\mathbf{x}, \mathbf{y}, \mathbf{z}) d \mathbf{z}.
\end{align}

Deep variational canonical correlation analysis (VCCA) \cite{Wang2016}, shown in Fig.~\ref{fig:DVCCA}, was extended from variational autoencoders (VAE) \cite{Kingma2013}. First, the mean vector $\boldsymbol{\mu}_{i}$ and the diagonal covariance matrix $\boldsymbol{\Sigma}_{i}$ of $\mathbf{x}_i$ are calculated by the encoder $\mathbf{f}_{x}$. Then, $L$ instances, $\{\mathbf{z}_{i}^{(l)}\}_{l=1}^L$, are randomly sampled from the distribution $\mathcal{N}\left(\boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)$. Finally, decoders $\mathbf{g}_{x}$ and $\mathbf{g}_{y}$ reconstruct instances $\mathbf{x}$ and $\mathbf{y}$, respectively.

The objective function of VCCA is:
\begin{align}
\min _{\mathbf{f}_x,\mathbf{g}_x,\mathbf{g}_y}& \frac{1}{N} \sum_{i=1}^{N} D_{K L}\left(q( \mathbf{z}_i | \mathbf{x}_i)| | p\left(\mathbf{z}_{i}\right) \right)+ \\ \nonumber
&\frac{\lambda}{N L} \sum_{i=1}^{N} \sum_{l=1}^{L} \left(\log p\left(\mathbf{x}_{i} | \mathbf{z}_{i}^{(l)}\right)+\log p\left(\mathbf{y}_{i} | \mathbf{z}_{i}^{(l)}\right)\right),
\end{align}
where the first term denotes the Kullback-Leibler (KL) divergence between the posterior distributions $q(\mathbf{z}_{i}|\mathbf{x}_{i})$ and $ p\left(\mathbf{z}_{i}\right)$, and the latter two terms denote the expectations of the log-likelihood under the approximate posterior distributions, which are equivalent to the reconstruction error. In testing, the mean vector $\boldsymbol{\mu}_{i}$  is used as the input features.

VCCA only takes the common latent variables $\mathbf{z}$ into consideration, which may not be sufficient to describe all information in $X$ and $Y$. In order to more adequately capture the private information of the views, Wang \emph{et al.} \cite{Wang2016} further proposed VCCA-private, shown in Fig.~\ref{fig:DVCCAP}. Its probabilistic model is defined as:
\begin{align}
p(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{h}_{x}, \mathbf{h}_{y}) &=p(\mathbf{z})p(\mathbf{h}_{x})p(\mathbf{h}_{y}) p(\mathbf{x} | \mathbf{z},\mathbf{h}_{x}) p(\mathbf{y} | \mathbf{z}, \mathbf{h}_{y}), \\
p(\mathbf{x}, \mathbf{y}) &=\iiint p(\mathbf{x}, \mathbf{y}, \mathbf{z},\mathbf{h}_{x}, \mathbf{h}_{y} ) d \mathbf{z} d \mathbf{h}_{x} d \mathbf{h}_{y}.
\end{align}

\begin{figure*}[htbp] \centering
\subfigure[VCCA]{\label{fig:DVCCA} \includegraphics[height=.35\linewidth,trim = 120 5 120 10, clip]{DVCCA}}
\subfigure[VCCA-private]{\label{fig:DVCCAP} \includegraphics[height=.35\linewidth,trim = 100 5 100 10, clip]{DVCCAP}}\centering
\caption{Schematic diagrams of DVCCA and DVCCA-private. The encoders $q(\cdot |x)$ and $q(\cdot |y)$ are used to fit the mean and covariance of the input, and the $p(x|\cdot)$ and $p(y| \cdot)$ are decoders for reconstructing input data. }
\end{figure*}

The corresponding objective function of VCCA-private is:
\begin{align}
\min_{\mathbf{f}_{x,h},\mathbf{f}_{x,z},\mathbf{f}_{y,h}, \mathbf{g}_{x}, \mathbf{g}_{y}} &\frac{1}{N} \sum_{i=1}^{N}
(D_{K L}\left(q\left( \mathbf{z}_{i} | \mathbf{x}_{i}\right)| | p\left(\mathbf{z}_{i}\right) \right) \\
&+D_{K L}\left(q\left( \mathbf{h}_{x,i} | \mathbf{x}_{i}\right)| | p\left(\mathbf{h}_{x,i}\right) \right)\nonumber \\
&+D_{K L}\left(q\left( \mathbf{h}_{y,i} | \mathbf{y}_{i}\right)| | p\left(\mathbf{h}_{y,i}\right) \right)
)\nonumber \\
&+ \frac{\lambda}{N L} \sum_{i=1}^{N} \sum_{l=1}^{L} \left(\log p\left(\mathbf{x}_{i} | \mathbf{z}_{i}^{(l)},\mathbf{h}_{x,i}^{(l)}\right)\right.\nonumber \\
&+\left.\log p\left(\mathbf{y}_{i} | \mathbf{z}_{i}^{(l)},\mathbf{h}_{y,i}^{(l)}\right)\right).\nonumber
\end{align}

In testing, the mean vectors output by the three encoders are stacked as new features. According to \cite{Wang2016}, VCCA-private is more effective than VCCA.

\subsection{Summary}

This section has introduced the traditional CCA and its six nonlinear extensions, three supervised extensions, and one sparse extension. Three ideas were used in the nonlinear extensions, namely, kernel-based KCCA, locality preserving LPCCA, and DNN-based DCCA (DCCAE, DVCCA). The first two are more suitable for small datasets. The supervised extensions (DisCCA, MLDA and MULDA) are closely related to LDA, because they all take the between-class and within-class scatter matrices into consideration.

\section{CCA-Based MVL for More than Two Views}\label{sec:MultiViewCCA}

This section summarizes five representative CCA approaches for more than two views, which are SUMCOR-GCCA, MAXVAR-GCCA, LS-CCA, TCCA, and DGCCA.

\subsection{SUMCOR-GCCA}

CCA maximizes the correlation between two views. For more than two views, a natural extension is to maximize the sum of the pairwise correlations \cite{Horst1961} \cite{Kettenring1971}, which is the idea of SUMCOR-GCCA.

Let $\left\{{X}_{j} \in \mathbb{R}^{d_{j} \times N}\right\}_{j=1}^J$ be a dataset containing $J$ mean-zero views, where $d_{j}$ is the feature dimensionality of View $j$. The objective function of SUMCOR-GCCA is:
\begin{align}
\max _{\left\{\mathbf{w}_{i}\right\}_{i=1}^{J}}\ & \sum_{i=1}^{J}\sum_{j=1}^{J}\mathbf{w}_{i}^{T}{X}_{i}{X}_{j}^{T}\mathbf{w}_{j}\\ \nonumber
s.t.\ &\mathbf{w}_{j}^{T}X_{j}X_{j}^{T}\mathbf{w}_{j}=1,\quad j=1,\cdots,J.
\end{align}

There is no closed-form solution for $\mathbf{w}_{j}$. Therefore, we followed \cite{Asendorf2015} and solved it using the Manopt \cite{Manopt} package\footnote{www.manopt.org}.

\subsection{MAXVAR-GCCA}

GCCA \cite{Carroll1968} assumes that each view can be generated from a set of multivariate latent variables ${G}= [\mathbf{g}_{1}, \mathbf{g}_{2},\ldots, \mathbf{g}_{N}]^{T} \in \mathbb{R}^{N \times K}$, in which all views are correlated. Its objective function is:
\begin{align}
\min_{G, \left\{W_{j}\right\}_{j=1}^{J}} \sum_{j=1}^{J}\left\|G-X_{j}^{T} W_{j}\right\|_{F}^{2} \quad  s.t.\ G^{T}G=I,
\end{align}
where ${W}_{j}\in \mathbb{R}^{d_{j} \times K}$ is the canonical matrix of View $j$, and the constraint guarantees that $G$ is a unit orthogonal matrix.

This objective is similar to a least squares problem. In order to avoid memory problems for large datasets, following \cite{Rastogi2015,Arora2012}, an SVD is first used to obtain a rank-$m$ approximation of each view:
\begin{align}
{X}_{j} \approx {U}_{j}{S}_{j}{V}_{j}^{T},\quad  j=1,\cdots,J, \label{eq:Xj}
\end{align}
where ${S}_{j}\in \mathbb{R}^{m \times m}$ is a diagonal matrix composed of the $m$ largest singular values, and ${U}_{j}\in {R}^{d_{j} \times m}$ and ${V}_{j}\in \mathbb{R}^{N \times m}$ are the corresponding left and right singular matrices, respectively.

$G$ then consists of the $K$ leading eigenvectors of matrix ${M}=\tilde{M}\tilde{M}^{T}$, where
\begin{align}
\tilde{M}&=\left[V_{1}T_{1},\cdots,V_{J}T_{J}\right]\in \mathbb{R}^{N\times mJ}.
\end{align}
It has been shown \cite{Rastogi2015} that the diagonal matrix $T_{j}$ satisfies $T_{j}T_{j}^{T}=S_{j}\left(S_{j}^{T} S_{j}+r_{j} I\right)^{-1} S_{j}^{T}$. Hence, given $S_j$ in (\ref{eq:Xj}), $T_j$ can be easily computed.

Once $G$ is obtained, $W_j$ can be computed as:
\begin{align}
W_{j}&=\left(X_{j} X_{j}^{T}+r_{j}I\right)^{-1} {X}_{j} G.
\end{align}

Benton \emph{et al.} \cite{Benton2016} also proposed view-weighted GCCA to consider the varying importance of different views.

\subsection{Least Squares based Generalized CCA (LS-CCA)}

LS-CCA \cite{Via2007} aims to minimize the distances among the canonical variables, so that different views are maximally overlapping with each other after mapping. Its objective function is:
\begin{align}
\min  _{\left\{\mathbf{w}_{j}\right\}_{j=1}^{J}} \ & \frac{1}{2 J(J-1)}\sum_{i=1}^{J}\sum_{j=1}^{J}\left\|X_{i}^{T} \mathbf{w}_{i}-X_{j}^{T} \mathbf{w}_{j}\right\|_{2}^{2} \\ \nonumber
s.t.\ &\frac{1}{J} \sum_{j=1}^{J} \mathbf{w}_{j}^{T}XX^{T} \mathbf{w}_{j}=1.
\end{align}

The $k$-th stacked canonical vector $\mathbf{w}^{\left(k\right)}=[\mathbf{w}_{1}^{\left(k\right)};\cdots; \mathbf{w}_{J}^{\left(k\right)}]$ can be obtained by performing a generalized eigen decomposition:
\begin{align}
\frac{1}{J-1}(R-D) \mathbf{w}^{(k)}=\lambda^{(k)} D \mathbf{w}^{(k)},
\end{align}
where
\begin{align}
R=
\begin{bmatrix}
\Sigma_{11} & \cdots & \Sigma_{1J} \\
\vdots & \ddots & \vdots \\
\Sigma_{J1} & \cdots & \Sigma_{JJ}
\end{bmatrix},\quad
D=
\begin{bmatrix}
\hat{\Sigma}_{11} & \cdots & \mathbf{0} \\
\vdots & \ddots & \vdots \\
\mathbf{0} & \cdots & \hat{\Sigma}_{JJ}
\end{bmatrix}.
\end{align}

Since the computational cost of this solution is high, LS-CCA can also be solved in an iterative manner, e.g., using partial least squares (PLS) \cite{Via2007}, so that the canonical vectors are obtained directly without using low rank approximation. Via \emph{et al.} \cite{Via2007} showed that LS-CCA and MAXVAR-GCCA are equivalent. Moreover, when $J=2$, LS-CCA degrades to CCA.

\subsection{Tensor CCA (TCCA)}

Traditional CCA optimizes the pairwise view correlation only, and ignores high-order statistics. TCCA \cite{Luo2015} directly maximizes the correlation of all views, using a high-order covariance tensor. Its objective function is:
\begin{align}
\max _{\left\{\mathbf{w}_{j}\right\}_{j=1}^{J}}\ & \rho\left(\mathbf{w}_{1}^{T}X_{1}, \mathbf{w}_{2}^{T}X_{2}, \ldots, \mathbf{w}_{J}^{T}X_{J}\right)\\ \nonumber
s.t.\  & \mathbf{w}_{j}^{T}XX^{T}\mathbf{w}_{j}=1, \ j=1, \ldots, J,
\end{align}
where $\rho$ denotes the correlation among the $J$ views. Luo \emph{et al.} \cite{Luo2015} showed that
\begin{align}\label{equ:TCCA}
\rho\left(\mathbf{w}_{1}^{T}X_{1}, \ldots, \mathbf{w}_{J}^{T}X_{J}\right)=\mathcal{C}_{1,2, \ldots ,J} \times_{1} \mathbf{w}_{1}^{T} \ldots \times_{J} \mathbf{w}_{J}^{T},
\end{align}
where $\times_{j}$ is the $j$-mode product, and $\mathcal{C}_{1,2, \ldots, J}$ denotes a $d_{1} \times d_{2} \times \cdots \times d_{J}$ covariance tensor. Let $\circ$ be the tensor product. Then,
\begin{align}
\mathcal{C}_{1,2, \ldots, J}=\frac{1}{N} \sum_{i=1}^{N} \mathbf{x}_{1 i} \circ \mathbf{x}_{2 i} \circ \ldots \circ \mathbf{x}_{J i},
\end{align}
where $\mathbf{x}_{ji}$ ($j=1,...,J$) denotes the $i$-th instance of View $j$. Let $\mathbf{u}_{j}=\hat{\Sigma}_{j j}^{1 / 2} \mathbf{w}_{j}$, and
\begin{align}
\mathcal{M}=\mathcal{C}_{1,2, \ldots, J} \times_{1} \hat{\Sigma}_{11}^{-1 / 2} \times_{2} \hat{\Sigma}_{22}^{-1 / 2} \ldots \times_{J} \hat{\Sigma}_{J J}^{-1 / 2}.
\end{align}
Then, (\ref{equ:TCCA}) can be rewritten as:
\begin{align}
\max_ {\left\{\mathbf{u}_{j}\right\}_{j=1}^{J}}\ & \mathcal{M} \times_{1} \mathbf{u}_{1}^{T} \times_{2} \mathbf{u}_{2}^{T} \ldots \times_{J} \mathbf{u}_{J}^{T} \\ \nonumber
s.t.\ & \mathbf{u}_{j}^{T} \mathbf{u}_{j}=1, \quad j=1, \ldots, J.
\end{align}

This problem can be solved by decomposing $\mathcal{M}$ into $K$ best rank-1 approximations $\left[\mathbf{u}_{1}^{k}, \mathbf{u}_{2}^{k}, \cdots, \mathbf{u}_{J}^{k}\right]_{k=1}^{K}$ with alternating least squares, and then stacking the corresponding vectors to obtain the canonical vectors $W_{j}=\left[\hat{\Sigma}_{jj}^{-1/2}\mathbf{u}_{j}^{1}, \cdots, \hat{\Sigma}_{jj}^{-1/2}\mathbf{u}_{j}^{K}\right]$ of View $j$.

\subsection{Deep Generalized CCA (DGCCA)}

DGCCA \cite{Benton2017}, shown in Fig.~\ref{fig:DGCCA}, is a DNN-based nonlinear extension of GCCA. Each view first passes through a multi-layer perceptron neural network to obtain nonlinear features $H_{j}=\mathbf{f}_{j}\left(X_{j}\right)$, $j=1,\cdots,J$. The objective of DGCCA is defined as the sum of the reconstruction errors between the canonical variables $X_{j}^{T}W_{j} $ and the common representation $G$:
\begin{align}
\min_{\left\{\mathbf{f}_{j}\right\}_{j=1}^{J}} \ \sum_{j=1}^{J}\left\|G-\mathbf{f}_{j}\left(X_{j}\right)^{T} {W}_{j}\right\|_{F}^{2} \quad  s.t.\ G^{T}G=I.
\end{align}

\begin{figure}[htpb] \centering
\includegraphics[width=.95\linewidth,clip]{DGCCA}
\caption{A Schematic diagram of DGCCA for $J$ views. Each view $X_{j}, j=1,\cdots,J$, first passes through a multi-layer neural network $\mathbf{f}_{j}, j=1,\cdots,J$, to obtain nonlinear features, and then GCCA is applied to project them onto a common correlated space.}
\label{fig:DGCCA}
\end{figure}

In training, the network's outputs for a mini-batch of instances are used as the inputs to GCCA, and then the reconstruction error between each view and the common representation is calculated to update the corresponding network parameters. In testing, the common representation $G$ or the stacked canonical variables $\left[{X}_{1}^{T}{W}_{1},\cdots,X_{J}^{T}W_{J}\right]$ can be used as the input to a classifier.

DGCCA is unsupervised. To utilize the label information, Benton \emph{et al.} \cite{Benton2017} added the labels' one-hot coding matrix as an additional view.

\subsection{Summary}

This section has introduced five representative CCA approaches for more than two views. To our knowledge, TCCA is currently the only approach to directly calculate the correlation among more than two views. All other approaches perform the computation indirectly.

\section{Application}\label{sec:application}

As a well-known approach for analysing correlations between two sets of variables, CCA and its extensions have been used in a variety of machine learning tasks. According to the type of problems, they can be divided into three main categories: multi-view recognition tasks, cross-modal retrieval/classificaion, and multi-view embedding.

\subsection{Multi-view Recognition Tasks}

In these CCA-based MVL applications, recognition tasks generally refers to supervised or unsupervised classification or regression. These methods have a wide range of applications in computer vision. To name a few, face recognition were achieved by treating the Daubechies wavelet transformed low-frequent images and the original images as two different views \cite{Sun2008, Sun2016}. Zheng \emph{et al.} \cite{Zheng2006} recognized expression of facial images by using the psychologist's semantic rating as the second view. To improve the performance of human gender recognition, Shan \emph{et al.} \cite{Shan2008} fused gait and face modalities at feature level with CCA.  Another challenging problem in computer version is pose estimation, which is try to compute the pose parameters of an image in the three-dimensional scene. A representative approach is first to construct connections between training images with pose parameters by using CCA-based MVL approaches \cite{Sun2007}, and then the pose parameters were estimated by calculating the neighbors of test images in projected space. Gait recognition \cite{Xing2016} across various viewing angles was accomplished by first finding a unified feature space, and then a nearest neighbor classifier was used.

In addition to computer vision applications, many other areas also used CCA-based MVL approaches. Take phonetic recognition as example, it has been shown that the phonetic recognition accuracy can be improved by learning both acoustic and articulatory features with CCA \cite{Wang2016, Wang2015Unsupervised, Arora2013}. There are many other similar MVL problems, include handwritten digit recognition \cite{Elmadany2016, Sun2008, Shin2011, Sun2016}, biometric structure prediction \cite{Luo2015}, text categorization \cite{Sun2008}, advertisement classification \cite{Luo2015}, and so on.

Generally, in a variety of multi-view recognition tasks, the CCA-based MVL approaches are first used to learn different new feature representations for each view, namely canonical variables. Once the canonical variables are obtained, there are two representative approaches for the subsequent learning. Take classification as example:
\begin{enumerate}
\item Directly put concatenated canonical variables of different views $[Z_{x}; Z_{y}]$ into SVM, KNN or other classifiers \cite{Luo2015,Shin2011,Shan2008}.
\item Take canonical variables of the first view $Z_x$ as input to a classifier \cite{Wang2016, Sun2008, Wang2015DCCAE}. There are two cases to choose this approach. The first is that in a real-world situation, only one view data is provided in the test. The second is to use target as the second view when applying the CCA-based approaches, so the prediction of test instance is accomplished by calculating the nearest neighbors in the first projected view.
\end{enumerate}

\subsection{Cross-modal Retrieval and Classification}

Due to the rapid growth of multimedia data on the Internet, cross-media retrieval has become a research hotspot in the fields of
artificial intelligence, multimedia information retrieval and computer vision. The main challenge of cross-media retrieval is the ``heterogeneity gap", which means that the data from different media types is inconsistent. This gap makes it diffcult to measure the cross-media similarity of instances \cite{Chi2017}.For example, a document containing images and text, and our goal is to retrieve text articles in response to query images and vice-versa. One of the traditional solutions is to find a common representation for different views with CCA \cite{Rasiwasia2010, Peng2017, Wu2006}. Firstly, taking images and text as two different views, the $k$-dimensional canonical variables of each instances are obtained with CCA-based MVL methods. Consider that all canonical variables from different views as coordinates in a common correlated space, as is shown in \ref{fig:CCA}. When given a query image, after mapping onto the common correlated space with canonical vectors,  the most closely matches can be obtained by computing its nearest neighbors in the range of all text canonical variables.

Cross-model classification \cite{Chandar2016} has the same data distribution assumption that the paired canonical variables of two different views are close to each other in the common correlated space. Therefore, a classifier trained on the first view can perform well in the second view.

\subsection{Multi-view Embedding}

One of the important applications of CCA-based MVL approaches is word embedding in many natural language processing tasks \cite{Dhillon2011, Lu2015, Faruqui2014}. Faruqui and Dyer \cite{Faruqui2014} proposed to learn words embeddings by adding multilingual context with DCCA, and found that the resulting embeddings could improve the performance on word similarity tasks. Benton \emph{et al.} \cite{Benton2016} learned multi-view embeddings of social media users to provided a solution for friend recommendation with GCCA \cite{Benton2017}.

In generally, multi-view embedding have two benefits. The first is to capture information from multiple views to obtain a representation with a better measure of similarity. The second is to learn a low-dimensional representation for high-dimensional multi-view data, which is critical to reducing the complexity of algorithms.

\subsection{Summary}
This section has introduced three main application of CCA-based MVL approaches, including multi-view recognition tasks, cross-modal retrieval/classificaion, and multi-view embedding. In addition to the applications described above, there are some other applications, including phonetic transcription \cite{Benton2017} and super-resolution images reconstruction \cite{An2014, Huang2010}, and so on.

\section{Conclusion} \label{sec:Conclusion}

CCA is widely used in MVL. However, traditional CCA has some limitations: 1) it can only handle two views; 2)  it can only optimize the linear correlation between different views; 2) it is unsupervised, and hence cannot make use of the label information. Many nonlinear, supervised, or generalized CCA extensions have been proposed for remedy. This paper has presented a survey of representative CCA-based MVL approaches.


\section{Introduction}

Shadowed by a century of neuronal recordings, astrocytes, the electrically silent non-neuronal cells in the brain, have remained absent from most studies on the biological principles of intelligence and from efforts to translate this understanding to computational primitives of artificial intelligence. Recent advancements in $Ca^{2+}$ imaging \cite{volterra2014astrocyte} and selective stimulation \cite{pereaoptogenetic} have signified the active signaling that occurs among astrocytes as well as between astrocytes and neurons \cite{volterra2005}. Astrocytes and other glial cells are now known to do far more than just providing nutritional and structural support to neurons \cite{volterra2005}. Among the many roles attributed to astrocytes in brain function and dysfunction \cite{barres2008}, a striking one is their neuromodulatory ability that originates at the cellular level, where synaptic plasticity \cite{HaydonCSH, Arlington} and neuronal synchronization \cite{Fellin, Knoxville} take place, and extends to the network level, where brain rhythms are observed \cite{LeoRhythms, fries2015}, which is a major component of behavioral functions including memory \cite{LeoMemory, adamsky2018astrocytic} and cognition \cite{Perea2014, han2013, fries2015}.

Astrocytic neuromodulation takes place in several brain regions \cite{araque2014gliotransmitters, Letellier, pereaoptogenetic} and has elevated the synapse to a "tripartite" unit (Fig. 1A), where neurons and astrocytes process and learn information independently and collaboratively \cite{perea2009tripartite}. The current hypothesis is that, by forming a large number of tripartite synapses, an astrocyte acts as a spatiotemporal integrator of the synaptic activity \cite{deemyad2018}. We have recently developed computational models of astrocytic cells that replicate their experimentally reported intra- and inter-cellular activity. Briefly, our models suggest intracellular mechanisms that astrocytes may use to encode and modulate both the excitatory \cite{Arlington, Knoxville} and the inhibitory \cite{Chicago} synaptic activity. The astrocytic models sense the neurotransmitters' spillover from the presynaptic neurons and release gliotransmitters to the postsynaptic neurons, suggesting a spatially constrained cascade mechanism for neuromodulation. Upon binding to astrocytic receptors, the presynaptic neurotransmitters produce an intracellular messenger molecule, Inositol 1,4,5-Triphosphate ($IP_3$); When the $IP_3$ concentration exceeds a threshold, an intracellular $Ca^{2+}$ wave is generated, having a similar temporal and spatial profile to the one observed in \textit{in-vitro} studies \cite{shigetomi2016probing}. This $Ca^{2+}$ wave triggers the release of gliotransmitters which in turn induce a slow inward current ($I_{SIC}$) into the postsynaptic neuron \cite{parpura2000physiological}.

Following the prevailing dogma that brain computation equals neuronal processing, neuromorphic processors are built on the premise that neurons are the \textit{only} computing unit. To overcome the intrinsic limitations of the connectionist models of neuronal computation, multiple neuromorphic chips have now become available \cite{davies2018loihi,merolla2014million,schemmel2010wafer,furber2014spinnaker}. These chips employ a time-dependent computational formalism, spiking neural networks (SNN), where asynchronous computing units are simulated as spiking neurons and memory is distributed at the neuronal synapses. Notably, Intel's Loihi processor is the most advanced large-scale neuromorphic hardware, offering multi-compartmental processing, on-chip learning and an unconstrained scalability \cite{davies2018loihi}. A handful of Loihi applications have already demonstrated unparalleled power efficiency in batch-mode image processing \cite{lin2018programming,davies2018loihi} and speech recognition \cite{blouw2018benchmarking} as well as in real-time control of a mobile robot, demonstrated by our Lab \cite{tang2019spiking}.

Alongside efforts to implement backprop algorithms in SNN \cite{shrestha2018slayer}, most neuromorphic algorithms are either restricted enough to be dictated by the underlying connectome associated with the targeted function \cite{tang2019spiking}, or simple enough to be trained via variations of spike-time-dependent plasticity (STDP), a biologically relevant yet computationally weak Hebbian-type rule \cite{diehl2015unsupervised,bing2018end}. By correlating pre- and post-synaptic activities, STDP rules contribute to the asynchronous parallelism of neuromorphic chips but they rather oversimplify other types of global, localized or temporal learning rules in the brain \cite{Fitzpatrick}. For example, STDP considers the synaptic activity solely within a single synapse and within a milliseconds-scale, contrary to what is happening in the brain. This inhibits SNN from taking place in regional or global tasks such as synchronization, a well-studied yet poorly-understood brain principle believed to convey information in space and time \cite{fries2015}. Further, by learning short-term correlations, STDP ignores events at larger time-scales, which are crucial for real-time applications. The asynchronous neuromorphic chips lack an activity-dependent mechanism to entrain their SNN into rhythmic activities, much like the brain does to efficiently and timely compute. Also, SNN inherited from the ANN the assumption that learning takes place only in the connection strengths between neurons, although structural plasticity rules that span beyond the single synapse have been known for many years to increase the computational capacity of neurons \cite{poirazi2001}. In fact, learning and memory are optimally balanced in networks that macroscopically operate at the edge of chaos \cite{legenstein2007edge}, a narrow dynamical regime largely ignored in learning algorithms, yet also exhibited by the brain \cite{friedman2012universal}. Overall, neuromorphic solutions seem to be missing a \textit{mesoscopic} learning mechanism that can combine the computational efficiency of having a global network goal with the versatility of a local, activity-dependent, plastic mechanism.

In this paper, we introduce spiking Neuronal-Astrocytic Networks (SNAN) and demonstrate some of the distinct computational and learning abilitities that astrocytes may add to neuromorphic hardware. We first describe our Loihi module\footnote{The module is available at https://github.com/combra-lab/combra\_loihi}, where astrocytes are used as an information processing unit, communicating with neurons in the tripartite synapses and with each other via intra- and inter-cellular $Ca^{2+}$ waves \cite{scemes2006astrocyte}. We then exhibit three example uses of our Loihi module, where astrocytes can extend learning at the meso-scale: Specifically, a) we show how astrocytes may modulate the neuronal component of the SNAN, by using their $Ca^{2+}$ dependent SIC to impose an on-demand synchronization across neuronal areas; b) We present how they can be used for single-shot pattern memorization via astrocyte-reinforced STDP; c) We finally show how a single astrocytic cell may be used to continuously monitor the neuronal component for its transition from order to chaos. With SNAN already showing their potential in a growing number of applications as reinforcement to neuronal networks \cite{porto2011artificial,valenza2013novel,LeoMemory}, the proposed Loihi module can become a computational framework for exploring and exploiting the unique computational principles that astrocytes are now known to exhibit.

\section{The Loihi Astrocytic Module (LAM)}

\begin{figure}
\includegraphics[scale=1]{paperfig1.png}
\caption{A. Schematic of the biological tripartite synapse. B. Its implementation on our Loihi Astrocytic Module, preserving the communication pathways between the astrocytic and the neuronal components of the SNAN.}
\end{figure}

\subsection{Module Overview}

\begin{figure*}
\includegraphics[scale=1]{paperfig2.png}
\caption{A feedforward SNAN where an astrocyte senses and imposes synchronous neuronal activity. A. SNAN Architecture. B. Activity of the astrocyte's main compartments: Integrating synaptic activity induces a SIC injected to postsynaptic neurons.}
\end{figure*}

Our Loihi Astrocytic Module (LAM) emulates the basic communication pathways between neurons and astrocytes in the tripartite synapse as well as the sub-cellular processes that astrocytes use to generate calcium waves in response to synaptic activity. The module's building block is the astrocytic cell that is created as an assembly of 4 Loihi compartments that approximate in hardware the leaky-integrate-and-fire neuronal model. We connected these compartments and tweaked their internal dynamics to represent the slower astrocytic processes (Fig. 1B). This ensures that the astrocytic component can be seamlessly integrated into any spiking neural network (SNN), by only declaring the connections between any of its spiking neurons and astrocytes.

\subsection{Module Implementation}
The LAM has 2 main components: (i) The input integration component that senses the neuronal activity by emulating the astrocytic receptors and $IP_3$ production, and (ii) the output generation component that modulates the neuronal component of the SNAN. The input integration component is comprised of 2 spiking Loihi compartments, the input Spike Receiver (SR) and the $IP_3$ integrator (IP3). The SR compartment integrates synaptic activity into IP3, resembling the behavior of astrocytic receptors. The IP3 dynamics are much slower than the neuronal dynamics, with their time-scale ranging from hundreds of milliseconds to tens of seconds. This allows for a temporal integration of the synaptic activity. When IP3 levels reach a threshold, the output component generates bursting spikes to either the postsynaptic neurons or the Loihi's reinforcement channel of neuronal connections. The output generation component has a dendritic tree where the non-spiking SIC generator compartment (SIC) feeds its voltage into the burst Spike Generator compartment (SG). The SIC produces a continuous voltage signal similar in shape to the biological $I_{SIC}$. The SG discretizes the signal into bursting spikes sent to the user-defined destinations.

\subsection{Developing a feedforward SNAN using the LAM}
To demonstrate the internal dynamics of the astrocyte compartments as the mechanisms that an astrocyte uses to impose synchronous activity in postsynaptic neurons \cite{Fellin}, we present a simple 2-layer feedforward SNAN (Fig. 2A). In this network, presynaptic neurons with Poisson spike trains were randomly connected to the same number of postsynaptic neurons. The astrocyte integrated the presynaptic activity and generated an IP3 spike at 6 seconds, which, in turn, triggered the SG bursting spikes and generated synchronous activity in the postsynaptic neurons for 400 milliseconds (Fig. 2B). The astrocytic dynamics are configurable in our astrocyte SDK that is described below.

\section{Astrocyte SDK}

\begin{figure}[b]
\includegraphics[scale=1]{sec3fig1.png}
\caption{Overall Astrocyte SDK Structure.}
\end{figure}

\subsection{SDK Architecture}

The Astrocyte SDK builds on the existing NxSDK Python-based programming model for Loihi \cite{lin2018programming} and it is designed to facilitate the integration of astrocytes with SNNs. The SDK's API provides access to 3 main classes, depicted in Figure 3:
\begin{enumerate}
  \item \textbf{AstrocytePrototype}: The prototype class for defining astrocytes allows users to configure the properties of each of the four compartments within the astrocyte, including the spiking thresholds for the IP3 compartment, the current decay of the SIC compartment and the weights between the compartments.
  \item \textbf{Astrocyte}: An astrocyte instance describes the internal compartments of the astrocyte model. This class provides the functions for specifying the presynaptic (input) neurons and the postsynaptic (output) neurons, to create tripartite synapses. Each instance of the Astrocyte class requires an AstrocytePrototype; If none is provided, a default prototype is created.
  \item \textbf{AstrocyteGroup}: A wrapper class that instantiates multiple astrocytes. The constructor accepts lists of prototypes, the size of the group of astrocytes and other configurations for mapping astrocyte instances to prototypes and Loihi's logical cores.
\end{enumerate}

\subsection{Automatic Parameter Setup}

Our Astrocyte SDK allows the user to define astrocyte behavior by controlling the following three parameters in the AstrocytePrototype class:

\begin{enumerate}
    \item \textbf{ip3\_sensitivity}: the weight between the SR and the IP3.
    \item \textbf{sic\_amplitude}: the maximum firing rate of the SG.
    \item \textbf{sic\_window}: time duration in milliseconds between the first and last of the spikes generated by the SG.
\end{enumerate}

These SDK level parameters are automatically mapped to Loihi parameters using a precalculated configuration table, generated using the $sic\_data\_table.py$ file located within the utils folder of the module. This script performs a brute force search of Loihi parameter configurations, saving all configurations that yield a minimal discretized SIC output, or a single burst spike. Finally, the closest matching configuration is selected based on the Euclidean distance metric defined as

\begin{equation*}
    \begin{split}
        cost = &(sic\_amplitude_{target} - sic\_amplitude_{config})^2 \\
        &- (sic\_window_{target} - sic\_window_{config})^2.
    \end{split}
\end{equation*}

Users can define their own range of values for the weights or current decays by regenerating the configurations table using the $sic\_data\_table.py$ file.

\section{Astrocytes induce Neuronal Synchronization}

\begin{figure}
\includegraphics[scale=1]{sec3fig2.png}
\caption{An example Astrocyte SDK program for integrating 2 astrocytes into an SNAN using the AstrocyteGroup.}
\label{fig:example_astrocyte_program}
\end{figure}

\begin{figure}
\includegraphics[scale=1]{sec3fig3.png}
\caption{An astrocyte group inducing synchronized activities in postsynaptic neurons.}
\end{figure}

We demonstrate how a group of astrocytes can be instantiated automatically, as described in Section 3.2, to synchronize sub-networks of neurons in an SNAN. We provide the code to create a single-layer feedforward SNN with 20 input neurons and 20 output neurons (Fig. 4). In this example, the synapses between the input and output neurons were randomly created with a probability of 8\% and a weight of 3. Briefly, to form the SNAN, we created the astrocyte group by defining an AstrocytePrototype (Fig. \ref{fig:example_astrocyte_program}a). We then created a group of astrocytes (Fig. \ref{fig:example_astrocyte_program}b) by specifying the NxNet instance it belonged to, the prototype(s) to use for instantiation and the number of astrocytes to create within the AstrocyteGroup. We also demonstrate how to form tripartite synapses based on the connection masks, indicating which of the possible neurons to connect to, and the weights for those connections (Figs. \ref{fig:example_astrocyte_program}c, d). Replicating experimental studies as well as our computational results \cite{Fellin, Knoxville}, the two Loihi astrocytes imposed distinct synchronization patterns in the neuronal component, forming two synchronized groups (Fig. 5).

\section{Astrocytes expand Hebbian Learning in Space and Time}

\subsection{Experimental and Computational Evidence}

Astrocytes are now known to take active part in several brain functions, from improving cognitive tasks and memory \cite{han2013} to suppressing futile unsuccessful behavior \cite{zebra}; Interestingly enough, they have also been found to not only be necessary but also sufficient for new memory formation \cite{adamsky2018astrocytic}. Astrocytes modify behavior by changing the network's structure in ways that cannot be captured by Hebbian-type rules. In fact, what distinguishes astrocytes from neurons is that they integrate synaptic activity over a much wider spatial scale and much longer temporal scale \cite{araque2014gliotransmitters}. This allowed us to explore possible mechanisms that astrocytes may use to increase or decrease synaptic weights, individually from neurons.

The newly identified mechanisms in astrocytes can be used to expand two main characteristics of Hebbian-type learning. First, astrocytes can be used to introduce regional learning: Hebbian learning is localized and lacks any knowledge of the neuronal activity in the vicinity of the synapse. Any change in synaptic weights is driven by the spiking activity of the pre- and post-synaptic terminals that form a single synapse. However, the spatial scale of astrocyte-induced plasticity spans from neighboring spines \cite{Zhang} to synapses found hundreds of micrometers away from the active synapse \cite{purinergic}, extending the spatial reach of plasticity. Second, astrocytes can integrate synaptic activity over time, ranging from hundreds of milliseconds to seconds, extending plasticity beyond the short-term correlations between pre- and post-synaptic activities. Therefore, a change of synaptic weights that currently depends only on very short temporal dependencies, can now be influenced by the activity in a much larger time window, allowing, e.g., for temporally delayed plasticity effects \cite{Serrano, Letellier}. Such computational primitives step away from the typical Hebbian learning rules and can endow astrocytes with unique computational and learning abilities that are applicable to SNN.

\subsection{Astrocytic-induced Heterosynaptic Depression (HSD)}

Heterosynaptic plasticity offers a straightforward astrocyte - implicated extension of Hebbian learning, in time and space. For example, astrocytes are found to detect the high (tetanic) activity of a pathway and respond to it by inducing presynaptic depression on the neighboring inactive pathways \cite{Serrano}. This astrocyte-driven Heterosynaptic Depression (HSD) may be used to complement the associative nature of STDP by counterbalancing the homosynaptic alterations, making the heterosynaptic modifications more prominent. In the next session, we show how HSD can use its spatial information to augment STDP, introducing single-shot memories.

\begin{figure*}
\includegraphics[scale=1]{paperfig4.png}
\caption{Single-shot pattern memorization and retrieval. A. Architecture of the SNAN with astrocyte-reinforced STDP. B. Activity of the sensory neurons and the memory neuron during learning and retrieval. C. Dynamics of the synaptic weights potentiated by the astrocyte: Switching STDP on and off is based on the temporal persistence of the learnable input.}
\end{figure*}

\subsection{Example Implementation of HSD in LAM to Introduce Single-Shot Memory}

We now describe how LAM supports astrocyte-driven plasticity mechanisms, to incorporate space and time in synaptic weight modification. To showcase this functionality, we created a 2-layer feedforward SNAN (Fig. 6A), with the astrocyte monitoring the activity of the sensory neurons. LAM uses the reward channels of Loihi's learning engine to implement the astrocyte-induced HSD on top of the regular neuronal STDP component.

The astrocyte-reinforced learning rule is defined as a sum of trace products \cite{lin2018programming},
\begin{equation}
    dw = \sum_i S_i \prod_j T_{ij},
\end{equation}
where $T_{ij}$ is the trace and $S_i$ is the learning rate.

The regular STDP component implemented on Loihi is defined as
\begin{equation}
    dw_{stdp} = a*x_1*y_0-b*x_0*y_1,
\end{equation}
where $x_0$ and $y_0$ represent the presence of pre- and post-synaptic spikes, respectively, and $x_1$ and $y_1$ are the pre- and post-synaptic spike traces, decaying exponentially.

The astrocyte-induced HSD component is defined as
\begin{equation}
    \begin{split}
        dw_{hsd} & = -(c*y_0 - d*x_0)*r_1 \\
        & = -c*y_0*r_1 + d*x_0*r_1,
    \end{split}
\end{equation}
where the reward spike trace $r_1$ represents the activity of the astrocyte. The HSD component decreased the weights of all neuronal synapses when the astrocyte was active. For more prominent heterosynaptic changes, the component induced slower weight decreases for sensory neurons with stronger activity and faster decreases for those with weaker activity.

We finally define the HSD-reinforced STDP on Loihi as
\begin{equation}
    dw = dw_{stdp} + dw_{hsd}
\end{equation}

\subsubsection{Experiment}

We demonstrate how LAM can be used to enhance the regular STDP learning rule by superimposing an HSD-reinforced component in a feedforward SNAN (Fig. 6A). The network consisted of 9 sensory neurons, each of which encoded one block of a 3x3 input grid into a Poisson spike train. The input neurons were connected to one output (memory) neuron. An astrocyte monitored the synaptic connections and modulated the memory neuron by injecting a SIC. We used 5 input patterns, with each pattern having 3 neurons active. All patterns had an overlap of one active neuron. The input neurons had a baseline firing rate of 5 Hz (white blocks) and an active firing rate of 100 Hz (black blocks). The Loihi parameters used in this implementation were as follows: $a=2^{-5}$, $b=2^{-6}$, $c=2^{-2}$, $d=2^{-1}$, the impulses for $x_1$ and $y_1$ were 16 whereas for $r_1$ was 8, and the decay time constants for all traces were 2. First, the network was trained on a desired pattern, which was presented for 2 seconds (memory learning). Then, the different overlapping patterns were presented to the network for 0.2 seconds each (memory retrieval).

\subsubsection{Results}
By virtue of its own slow dynamics, the astrocyte switched on and off STDP learning, controlling when the network memorized the pattern. Specifically, when an input pattern was presented to the network long enough (in the order of seconds), this persistent neuronal activity triggered the astrocytic SIC. Further, as it integrated inputs from regional synapses, the astrocyte became active only when the overall input was strong enough. When the astrocyte was activated, all synapses encoding inactive blocks decreased their weights to negative values through HSD, thereby inhibiting the memory neuron (Fig. 6C). This created a distinction between regular STDP and the astrocyte-reinforced STDP: With regular STDP, any active block that overlaps among patterns (e.g. the central block in the first two and the fourth patterns - Fig. 6A) increases the weight of its synapse every time a pattern is presented. In the presence of HSD, when active blocks overlapped across patterns, the memory neuron received inhibitory input from the sensory neurons corresponding to active grid blocks forming a pattern different than the one learned. During memory retrieval, the memory neuron fired maximally when the learned pattern was presented (Fig. 6B). Neurons that were irrelevant to the learned pattern had negative weights and inhibited postsynaptic firing. With no astrocytic contribution, this learning rule downgraded to a regular STDP (Fig. 6C) that does not change the synaptic weights when the postsynaptic activity is low. In other words, the astrocyte controlled \textit{when} the network memorized a new pattern, based on the input's persistence. Overall, LAM enables an astrocytic gating mechanism for neuronal learning that expands the regular STDP to space and time.

\section{Astrocytes Learn to Monitor Transition into Chaos}

\subsection{LAM Support of Astrocytic Learning}

LAM enables the creation of astrocytes with custom input plasticity mechanisms that are limited only by the set of operations supported by Loihi. It is known that astrocytes perform large-scale spatial and temporal integration of synaptic activity \cite{perea2009tripartite,araque2014gliotransmitters} and exhibit bidirectional homeostatic plasticity (BHP), where astrocytic connections continuously increase/decrease sensitivity to synapses with low/high activity over the course of several hours \cite{xie2012}. We incorporated rate-dependent astrocytic learning in LAM as follows:
\begin{equation}
    \begin{split}
        dt &= f_{r}(x_0)\\
        dw &= f_{w}(w,t),
    \end{split}
\end{equation}
where $f_{r}$ dynamically integrated synaptic input spikes sensed by the astrocyte, resulting in a spike rate approximation variable $t$ for each tripartite synapse, and $f_{w}(w,t)$ implemented astrocytic weight learning defined in the astrocyte input compartment. In the next session, we demonstrate how the proposed astrocytic model with rate-based learning may detect the transition of the neuronal network, from order to chaos.

\subsection{An Astrocytic Model For Detecting Network Chaos}

SNNs fall under the general category of complex nonlinear dynamical systems, which maximize their computational capacity when they operate at the edge of chaos, as shown by studies ranging from cellular automata to boolean networks \cite{packard1988adaptation,ribeiro2008}. It is at the critical point between ordered and chaotic dynamics where balance of robustness and versatility emerges not only in brain-imitating networks \cite{legenstein2007edge} but also in the brain itself \cite{friedman2012universal}. Interestingly, the design of learning algorithms largely does not take into account this general computing principle that would enable SNNs to operate at the edge of chaos. In our ongoing effort to develop self-tuning near-critical SNNs, we proposed and developed a homeostatically plastic, astrocyte model capable of detecting and signaling the approach of network chaos, detailed in our modeling paper \cite{vlad2019unpublished}, and abstracted here as,
\begin{equation}
f_{astro}(r_{i}, {\hat{r}}_{i}) = g \Bigg( \bigg( \frac{1}{N} \bigg) \sum_{i=1}^{N} \frac{r_{i}}{r_{max}} \log{\Bigg(   \frac{1}{\frac{{\hat{r}}_{i}}{r_{max}}}\Bigg)} \Bigg),
\end{equation}
where  $r_i$ is the short-term synaptic firing rate, ${\hat{r}}_{i}$ is the long-term synaptic firing rate, $r_{max}$ is the maximum synaptic firing rate, $g$ is a nonlinear activation function, and $N$ is the number of synaptic inputs into the model.  Intuitively, this model compares the run-time, short-term synaptic firing rates, $r_i$, to the memorized, long-term synaptic firing rates, ${\hat{r}}_{i}$, and aggregates the extent of the individual synaptic differences in the overall astrocytic output frequency, $f_{astro}$.

\begin{figure*}
\includegraphics[scale=1]{paperfig5.png}
\caption{Homeostatic astrocyte on Loihi sensing transition from ordered to chaotic network state. A. SNAN architecture showing training and testing phases. B. Learned weights for all the neuronal-astrocytic connections. C. (upper row) Ordered and chaotic neuronal network state (middle row) number of active neurons in both states showing similar neuronal activity (bottom row) astrocytic activity during the transition from order to chaos.}
\end{figure*}

\subsection{Implementation of Astrocytic Learning in LAM to Detect Chaos}

We implemented an SNAN with a single-layer of neurons connected to an astrocyte endowed with the BHP learning mechanism (Fig. 7A). To approximate the BHP rule on Loihi, we decomposed Equation (5) into a weighted sum,
\begin{equation}
    \begin{split}
        f_{astro}(r_{i}, {\hat{r}}_{i}) = g \Bigg(\eta \sum_{i=1}^{N} W_i*N_i\Bigg) \\
        W_i  = \log{\Bigg(   \frac{1}{\frac{{\hat{r}}_{i}}{r_{max}}}\Bigg)} ,\ N_i  = \frac{r_{i}}{r_{max}}.
    \end{split}
\end{equation}

The input $N_i$ was represented as the integrated current of real-time presynaptic spikes $I_i$ and the parameter $W_i$ was the weight of the neural-astrocytic connection $w_i$, described by
\begin{equation}
Activity_{astro} = Astrocyte(\sum_{i=1}^{N}w_i*I_i).
\end{equation}

Following the general form in Equation (5), the weight $w_i$ was learned based on the BHP learning rule, defined as,
\begin{equation}
    \begin{split}
        dt &= a*x_0\\
        dw &= b*u_k*(w_{max}-w)-c*u_k*t,
    \end{split}
\end{equation}
where $t$ represents the estimated long-term presynaptic neuron firing rate, and $w$ approximates the inverse relationship between the neuronal-astrocytic weight and the presynaptic firing rate described by $W_i$. Intuitively, the learning rule kept decreasing/increasing $w$ for high/low $t$ values in a predefined learning time window, where $t$ did not saturate (Fig 7B.).

\subsubsection{Experiment}

The astrocyte was trained using the BHP rule (Eq. (8)) on a single layer of spiking neurons driven by activity that was generated off-chip. To represent the neuronal activity, we used a non-isotropic Ising model, which is described in \cite{vlad2019unpublished}. Briefly, we modeled synaptic activity as 256x256 magnetic spins positioned on a 2-dimensional lattice with non-uniform, clustered couplings between spins resulting in stationary spin state patterns. We used the Markov Chain Monte Carlo (MCMC) algorithm to simulate the evolution of spin states by minimizing the system's overall energy. We determined ordered and chaotic Ising states (Fig. 7C) using the magnetic susceptibility measure \cite{yeomans1992statistical}. To drive the neurons on Loihi, we downsampled 42x42 spins from the original Ising system, evenly distributed throughout the Ising lattice, and converted each binary Ising spin state to a corresponding neuronal spike state. All Ising spins were updated and transformed into neuronal spikes every 5 milliseconds, keeping the maximum neuron firing rate to 200 Hz (Fig. 7C). The astrocyte's receptor weights were trained for 25 seconds on neuronal dynamics driven by ordered Ising spin evolution. Then, the model's activity was tested with respect to both states for 25 seconds each. The parameters for the BHP rule were $a=2^{-6}$, $b=2^{-2}$, $c=2^{-3}$, $w_{max}=16$, and $k=4$ in $u_k$ which controlled the $w$ to only update every $2^4$ learning epochs.

\subsubsection{Results}

The frequency of the astrocytic calcium wave increased when the neuronal component transitioned from an ordered to a chaotic state of neuronal firing (Fig. 7C). This nearly 200\% increase in frequency resulted from the weights that were learned through the BHP rule, since network activity levels remained unchanged throughout the transition (middle row, Fig. 7C). Specifically, the learned weights generated by the BHP inversely correlated with the long-term mean firing rates of the Ising neurons during training (Fig. 7B), imitating the functional form of $W_i$ in Equation (6). As the minimally active neurons increased their firing rates with the transition to chaos (first row, Fig. 7C), their higher weights induced increased input current to the astrocyte, which resulted to increased astrocyte activity. That is how the transition between states was monitored by a single astrocytic cell.

\section{Discussion and Conclusion}

Here, we introduced LAM, a module that enables the integration of astrocytes to Loihi, endowing to SNN the computational primitives recently identified in these long-neglected non-neuronal cells. Adding astrocytes to neuromorphic processors, which are built on the premise that the neuron is the \textit{only} computational unit in the brain, has a great potential to advance neuromorphic computing. To demonstrate this potential, we demonstrated three example uses of astrocytic computation, based on several pathways for astrocyte-neuron interaction and neuromodulation.

Following biological evidence that astrocytes exhibit homosynaptic \cite{Perea2014} and heterosynaptic \cite{Serrano} plasticity mechanisms, at both presynaptic \cite{prePlast} and postsynaptic \cite{postPlast2} sites, LAM supports (i) HSD, shown here to be implicated in single shot pattern memorization, and (ii) BHP, demonstrated as they key mechanism for monitoring the transition of the neuronal network to chaos. These astrocytic learning primitives showcased how network learning can extend beyond updating the weights of the neuronal synapses and onto the wider operating time and space of astrocytes.

In addition, the astrocytes in LAM were used to induce neuronal synchronization, replicating both experimental and computational studies \cite{Knoxville,tewari2013possible}. Astrocytes are long-implicated in network frequency modulation \cite{tewari2013possible,LeoRhythms}, improving network polychronicity \cite{valenza2013novel} and enabling memory transition in attractor networks \cite{LeoMemory}. In our example, time and space became learnable parameters beyond the reach of Hebbian-type learning, allowing the introduction of “on-demand” synchronization or binding, principles in the brain that are long-believed to process information efficiently but are missing in asynchronous neuromorphic hardware.

Finally, we showed here how astrocytes in LAM may induce heterosynaptic depression constrained within their microdomain. By introducing spatial processing, astrocytes can further induce local plasticity and, through it, synaptic clustering, which is a computational primitive for robust brain computation that is also missing in neuromorphic solutions.

In taking our next steps to expand LAM, it has not escaped our attention that neuromorphic chips may be used as a computational framework enabling astrocytic functions to be both explored in brain hypotheses and exploited for computing purposes. Currently, LAM develops astrocytes as dimensionless units, having no inter-unit communication. A natural extension for our module is to extend astrocytic computation to (i) the subcellular level and (ii) network level. Both imaging \cite{DiCastro} and stimulation techniques \cite{Zhiping} provide evidence that the majority of astrocytic activity is spatially confined to functionally independent subcellular astrocytic compartments, rarely spanning the whole cellular domain. Adding dimensions in our astrocytic cells and networks would enable them to operate, process and learn, on a semi-local, "regional" or long-range scale. In fact, astrocytes, coupled through inter-astrocytic gap junctions, form a network syncytium, thereby enabling long distance communication on top of neuronal networks \cite{Bell}. Enabling both multi-compartment processing and inter-astrocytic communication would further bridge the computation capabilities of LAM with those of biological astrocytes.

Overall, although they probably require a lot more small insights before they can be fully integrated to SNN, the fact that our Loihi-run SNAN already emulate biological principles associated with robustness and versatility, indicates that the integration of astrocytes to neuromorphic chips is a direction worth exploring.

\begin{acks}
This research is funded by Intel's NRC Research Award.
\end{acks}