\section{Introduction}
An essential focus in computer vision, information-theoretic learning, and computational neuroscience has been centered on the question: how can a learning system discover the latent structure of the real world through sensors? This paper proposes a novel and principled learning framework built upon the bi-dimensional nature of random processes (r.p.): the probability measure of the random variables (r.v.) and their statistical dependence across the functional space.

\textcolor{black}{By definition, an r.p. is a sequence of Borel measurable random variables $(\Omega)$ over a functional space $( \mathcal{F})$, indexed by a continuous or discrete variable $t$ such as time instances of a signal or spatial coordinates of an image. For the discrete-time stochastic process $\mathbf{x} = \{\mathbf{x}(t, \omega)\}_{t=1}^{T}$ where $\omega$ is a subset of the common sample space $\Omega$}, there are two different components in the r.p.: the ensemble of realizations $x_1, x_2, \cdots, x_N$ of the r.p. sampled from the joint probability density function (\textit{pdf}), and the temporal or spatial domain $x_n(1), x_n(2), \cdots, x_n(T)$ for every realization $x_n$, containing a rich latent functional structure. For simplicity, we will omit the dependence on $\omega$ for the remaining presentation. Statistical dependence has been a fundamental approach to characterize the r.p. functional structure. Our approach is closely related to the Wiener filtering theory (WF), Gaussian processes (GP), and the alternating conditional expectation (ACE) algorithm.

The r.p. covariance is the fundamental statistical dependence measure utilized by Parzen, among many other researchers, for statistical inference under the linear model and Gaussian statistics (\cite{parzen1961approach, parzen1999stochastic}). The WF, as a supervised modeling framework, uses a linear model to minimize the mean squared error (MSE) or the error projected on a space defined by the covariance of the r.p. (\cite{wiener1949extrapolation, hagenblad1999aspects, hagenblad2000maximum, wills2013identification}). However, the WF requires a stringent wide-sense stationarity assumption, and there are no known closed-form solutions to find nonlinear WF extensions. Alternatively, discrete statistical signal processing solutions based on search techniques in reproducing kernel Hilbert spaces (RKHS) give rise to the kernel adaptive filters (KAF) (\cite{Principe2010information, principe2011kernel}). The KAF uses nonlinear kernel functions to estimate the mean of the joint \textit{pdf} over pairs of samples at two time intervals. Multilayer neural networks (either convolution neural networks (CNN) for spatial r.p. or recurrent neural networks (RNN) for temporal r.p.) can also be substituted for KAF. \textcolor{black}{An active research area is estimating the inverse covariance matrix for high-dimensional r.p. by applying optimization and sparse matrices techniques (\cite{yuan2010high, ravikumar2011high, pourahmadi2013high, bollhofer2019large}). However, the linearity of the covariance is still a shortcoming. The log-determinant of the covariance plays a central role in these approaches as a matrix rank approximator (\cite{fazel2003log}), which is also used recently for explaining neural networks (\cite{chen2018harnessing, yu2020learning, chan2022redunet}).}

The GP establishes a powerful Bayesian framework by hypothesizing a fixed Gaussian prior with a joint $pdf$ over a functional space, providing a variational interpretation of the kernel function (\cite{williams1998prediction}). Similar to other Bayesian inference frameworks, the GP implements maximum a posteriori (MAP) estimation by minimizing the negative log-likelihood (NLL) rather than constructing an error r.v.. The theoretical importance of GP resides in the eigenanalysis of the Gaussian kernel function (\cite{zhu1997gaussian, williams2006gaussian}). The Bayesian system contains the $pdf$ of the r.p., and the eigenanalysis of the kernel function characterizes the functional structure of the r.p.. However, there are two practical issues with GP when modeling non-stationary and non-Gaussian real-world r.p.. First, a stationary assumption is imposed when the prior is fixed. Second, the autocorrelation function that evaluates the posterior over the functional space, which corresponds to correntropy (\cite{Principe2010information}), only evaluates the joint distribution mean along the bisector of the joint space in a region specified by the kernel bandwidth. Moreover, the Nystr{\"o}m method (\cite{baker1977numerical}) to compute the eigenfunctions is still restricted. More specific quantifiers of the joint density require higher order moments of the $\textit{pdf}$ and are very computationally intensive. Therefore, the GP cannot provide a complete description of statistical dependence. \textcolor{black}{Many recent methods are proposed to improve GP kernels, including learning kernels with neural networks (\cite{wilson2016deep}) and selecting kernels with Bayesian variational methods (\cite{wilson2016stochastic, malkomes2016bayesian}). }

\textcolor{black}{An alternative to quantify higher order statistics of the input data is to substitute the MSE cost function by Shannon's mutual information (MI) (\cite{shannon1948mathematical}), which is currently the most widely used methodology to quantify statistical dependence. Methods for quantifying MI include estimating directly the \textit{pdf} from discrete samples (\cite{parzen1962estimation}; \cite{Liamestimate}; \cite{Principe2010information}; \cite{hu2021training}), employing statistical embeddings by projecting discrete samples to an RKHS (\cite{bach2002kernel}; \cite{gretton2007kernel}; \cite{sriperumbudur2010hilbert}), and maximizing the variational lower bound with a universal functional approximator (\cite{nguyen2010estimating}; \cite{belghazi2018mutual}; \cite{birrell2020f}; \cite{sreekumar2022neural}). Among these methods, kernel independent component analysis (KICA) (\cite{bach2002kernel}) is one of the first approaches to establish the relationship between the covariance functional in RKHS and the mutual information for characteristic kernels. Mutual information neural estimator (MINE) (\cite{belghazi2018mutual}) uses a neural network as a function approximator for maximizing the variational lower bound.} However, unlike KAF and GP, these approaches are not adequate for high-dimensional and non-stationary r.p., as they provide only a scalar-valued quantification of statistical dependence. As we can expect, a single scalar is unable to quantify the complex dependence structure in the functional space.

\textcolor{black}{Instead of using MI, the maximal correlation functional proposed by Alfred R\'enyi (\cite{renyi1959measures}) is preferable because it extends to r.p. naturally.} Given any two r.v. $\mathbf{x}$ and $\mathbf{u}$, the maximal correlation is the supremum of the correlation coefficient between $f(\mathbf{x})$ and $g(\mathbf{u})$ over all feasible scalar-valued functions $f$ and $g$. The alternating conditional expectations (ACE) (\cite{breiman1985estimating}) is one of the first algorithms implementing the maximal correlation using bidirectional recursive updates. At each iteration of the updates, the function $f(x)$ is updated to match the conditional mean $\mathbb{E}_{\mathbf{u}}[g(\mathbf{u})|x]$, and then $g(u)$ is updated to match $\mathbb{E}_{\mathbf{x}}[f(\mathbf{x})|u]$. The recursive updates reach the equilibrium as the two functions become the solutions for the maximal correlation. \textcolor{black}{However, ACE applies only to r.v. as it uses only the correlation coefficient and two scalar-valued functions. Moreover, the ACE applies only for discrete r.v., and fails to generalize to continuous r.v. because of the function approximator's limitations. One fundamental work discusses the probabilistic relations between the maximal correlation and the covariance functional of multivariate Gaussian distributions (\cite{bach2002kernel, bach2005probabilistic}). More recently, deep networks are used to approximate the maximal correlation~(\cite{andrew2013deep, hu2021mimo}). Similar works also include learning the Gaussian universal features with the maximal correlation (\cite{huang2018gaussian, huang2019universal}). The common limitation of these methods is that they are still limited to discrete variables. As we will see, kernel analysis applies naturally to continuous variables.}

To extend the theory of GP and ACE to r.p. and variational inference, our approach focuses on a Bayesian framework based on the maximal correlation and provides a totally new interpretation because it is based on the joint \textit{pdf} without imposing Gaussianity. Some definitions are necessary and will constitute the first theoretical part of the paper. \textcolor{black}{Given two r.p. $\mathbf{x} = \{\mathbf{x}(t)\}_{t=1}^{T_1}$ and $\mathbf{u} = \{\mathbf{u}(t)\}_{t=1}^{T_2}$, the former is called the primary r.p. (the input to a learning system), and the latter is called the reference r.p., normally class information in classification or target information in regression.} Let us construct two internal latent variables $\mathbf{x}'$ and $\mathbf{u}'$ to help characterize the statistical dependence between $\mathbf{x}$ and $\mathbf{u}$. The creation of the internal variables leads to a novel Bayesian variational system, named a bidirectional recursion. Assuming the joint density $p(\mathbf{x}=x, \mathbf{u}=u)$ exists between the two external r.p., the conditional densities can naturally be applied in two directions: the forward recursion $p(\mathbf{u}=u|\mathbf{x}=x)$ and the backward recursion $p(\mathbf{x} = x|\mathbf{u}=u)$. The internal variable $\mathbf{x}'$ is obtained by applying both recursions and marginalizing over $\mathcal{U}$ as $p(\mathbf{x} = x, \mathbf{x}' = x') = \int_{\mathcal{U}} p(\mathbf{x}' = x'|\mathbf{u}=u) p(\mathbf{u}=u|\mathbf{x}=x) p(\mathbf{x}=x) du$ to yield the joint density $p(\mathbf{x} = x, \mathbf{x}' = x')$ between the external and latent internal variables. Similarly, the internal variable $\mathbf{u}'$ can be obtained by marginalizing over $\mathcal{X}$ to yield the joint density $p(\mathbf{u} = u, \mathbf{u}' = u')$. This procedure forms a bidirectional variational system involving four variables $\mathbf{x}$, $\mathbf{x}'$, $\mathbf{u}$, $\mathbf{u}'$. Similar to GP, each of the joint $pdf$ $p(\mathbf{x} = x, \mathbf{x}' = x')$ or $p(\mathbf{u} = u, \mathbf{u}' = u')$ of the variational system defines a density kernel function, named the cross density kernel function (CDKF). The eigenspectrum $\{ \lambda_i \}_{i=1}^\infty$ of the CDKF is a novel statistical dependence measure applicable to r.p.. The eigenfunctions $\{\phi_i(x)\}_{i=1}^\infty$ and $\{\psi_i(x)\}_{i=1}^\infty$ of each CDKF characterize the functional structure of the respective r.p.. The bidirectional recursion can be seen as a major improvement over GP. Instead of using a fixed Gaussian prior, the internal variables $\mathbf{x}'$ and $\mathbf{u}'$ are learned from the data and have the same marginal distribution of $\mathbf{x}$ and $\mathbf{u}$ respectively, even when the data is no longer stationary. Moreover, unlike any kernel functions used in GP, our proposed CDKF characterizes the multivariate \textit{pdf} of each r.p. instead of the autocorrelation.

\textcolor{black}{In the second section, a learning framework is presented to implement the bidirectional recursion from r.p. realizations with the goal of estimating the latent internal variables and solving the CDKF eigenproblem.} Our implementation is very similar to the recursive updates in the ACE but no longer limited to two scalar-valued functions. Instead, two multivariate neural networks are employed to produce two families of functions $\{f_i\}_{i=1}^K$ and $\{g_i\}_{i=1}^L$ for approximating the eigenfunctions of the two CDKF. Instead of employing the correlation coefficient score as the ACE, the two networks are updated with recursive updates to satisfy the Wiener solution, bidirectionally. We name this new algorithm the functional maximal correlation algorithm (FMCA). We will theoretically show that when the recursive updates of FMCA converge, we obtain the eigenvalues and the eigenfunctions of the CDKF. \textcolor{black}{The \textit{pdf} component of the r.p. is contained in the joint density in CDKF, and the functional structure of the r.p. is contained in the eigenfunctions.} This cannot be achieved by either the ACE or by approximating MI.

The existence of \textcolor{black}{two external variables} implies a supervised framework for estimating statistical dependence and the dual architecture creates a latent space of a user-defined dimensionality by projecting nonlinearly the primary r.p. $\mathbf{x} = \{\mathbf{x}(t)\}_{t=1}^{T_1}$ concatenated with the reference user selected r.p. $\mathbf{u} = \{\mathbf{u}(t)\}_{t=1}^{T_2}$. Examples of the reference r.p. in supervised learning are multivariate labels for classification, target time series for regression and adaptive filtering. However, when we select the reference r.p. as non-structured white noise, the learning problem is no longer supervised. This case yields the factorial codes (FC) for the input data and can be used directly for variational inference with an external head network. Section~\ref{section_fmca} presents the properties and the analysis of FMCA, where a bidirectional neural network architecture is required to implement FMCA from realizations. Section~\ref{validation_results} and Section~\ref{real_world_applications} present a set of experiments with synthetic and real-world datasets to exemplify the performance of the technique in two problem domains: classification and variational inference. We compare FMCA with alternatives being used in the machine learning literature for medium-size datasets to show its scaling-up capability. The manuscript will end with a conclusion.

\section{A Novel Variational Approach to Measure Statistical Dependence}
\label{bidirectional_system}

The framework is presented in two parts. In this section, the focus is on the theoretical foundations of the variational system and the density kernel. Assuming that all \textit{pdfs} of the given r.p. are available, how can the statistical dependence within and between two r.p. be properly quantified? Unlike current methodologies, the proposed approach is to construct a Bayesian system and create latent internal variables that best reflect the statistical dependence.

\subsection{Definition of the bidirectional recursion and the cross density kernel function}

Suppose two r.p. $\mathbf{x} = \{\mathbf{x}(t)\}_{t=1}^{T_1}$ and $\mathbf{u} = \{\mathbf{u}(t)\}_{t=1}^{T_2}$ are given with a fixed and finite length $T_1$ and $T_2$. For the proofs, the joint densities and the conditional densities associated with the two r.p. are assumed to exist and are Lebesgue measurable. No stationarity condition is necessary. The variational system is defined as follows.

\begin{figure}[t]
\captionsetup{type=figure}
\centering
{\includegraphics[width=1\textwidth]{{./used_figures/pipelinefigure}.pdf}}\\
\captionof{figure}{\small Diagram of the bidirectional recursion process. (a) Given two r.p. $\mathbf{x}$ and $\mathbf{u}$ with a joint density, we view the two conditional densities as bidirectional iterative functions. By recursively applying the two conditional mappings, we create two internal variables $\mathbf{x}'$ and $\mathbf{u}'$. (b) The joint densities $p(x, x')$ and $p(u, u')$ obtained by applying the marginalization evaluate the level of independence in the space of the internal variables, drawn in red and blue regions. (c) Each joint density induces the density kernel function CDKF. By applying Mercer's theorem, the eigenvalues for the statistical dependence measure and the eigenfunctions for the approximation are obtained.}
\label{figureprocess}
\end{figure}

\begin{definition}
(Bidirectional recursion) Suppose two r.p. $\mathbf{x} = \{\mathbf{x}(t)\}_{t=1}^{T_1}$, $\mathbf{u} = \{\mathbf{u}(t)\}_{t=1}^{T_2}$ and their joint density $p(\mathbf{x} = x, \mathbf{u} = u)$ are given. The two conditional densities have the forms
\begin{equation}
\begin{aligned}
p(\mathbf{x} = x|\mathbf{u} = u) = p\Big(\mathbf{x}(1) = x(1), \cdots, \mathbf{x}(T_1) = x(T_1)\big|\mathbf{u}(1) = u(1), \cdots, \mathbf{u}(T_2) = u(T_2)\Big), \\
p(\mathbf{u} = u|\mathbf{x} = x) = p\Big(\mathbf{u}(1) = u(1), \cdots, \mathbf{u}(T_2) = u(T_2)\big|\mathbf{x}(1) = x(1), \cdots, \mathbf{x}(T_1) = x(T_1)\Big).
\end{aligned}
\label{forward_backward}
\end{equation}
To simplify the explanation we attribute different directions to each. The conditional density $p(\mathbf{u} = u|\mathbf{x} = x)$ is interpreted as an iterative function in the forward direction, and  $p(\mathbf{x} = x|\mathbf{u} = u)$ as an iterative function in the backward direction. The bidirectional recursion is a variational system constructed by recursively applying both conditional mappings. Starting with the variable $\mathbf{x}$ (or $\mathbf{u}$), define two internal variables $\mathbf{x}'$ or $\mathbf{u}'$ by computing the joint $pdf$
\begin{equation}
\begin{aligned}
p(\mathbf{x}=x, \mathbf{x}'=x')  &= \int_{\mathcal{U}} p(\mathbf{x}' = x'|\mathbf{u} = u) p(\mathbf{u} = u| \mathbf{x} = x) p(\mathbf{x} = x) pu,\\
p(\mathbf{u}=u, \mathbf{u}'=u') &= \int_{\mathcal{X}} p(\mathbf{u}' = u'|\mathbf{x} = x) p(\mathbf{x} = x| \mathbf{u} = u) p(\mathbf{u} = u) px.
\end{aligned}
\label{cross_densities}
\end{equation}
These two joint densities~\eqref{cross_densities} are named the cross densities. The four variables $\mathbf{x}, \mathbf{u}, \mathbf{x}', \mathbf{u}'$ form the variational system, which is named the bidirectional recursion. For convenience, the \textcolor{black}{following notations are employed}
\vspace{5pt}

(i) Marginal densities: $p(x):=p(\mathbf{x}=x),\; p(u):=p(\mathbf{u}=u);$ \vspace{5pt}

(ii) Conditional densities: $p(u|x):=p(\mathbf{u} = u|\mathbf{x} = x),\; p(x'|u):=p(\mathbf{x} = x'|\mathbf{u} = u) $; \vspace{5pt}

(iii) \textcolor{black}{The joint and the conditional cross densities}:
\begin{equation}
\begin{gathered}
p(x, x'):= p(\mathbf{x}=x, \mathbf{x}'=x'), p(u, u') := p(\mathbf{u}=u, \mathbf{u}'=u'),\\
p(x'|x) := p(\mathbf{x}'=x'|\mathbf{x}=x), \,\,p(u'| u) := p(\mathbf{u}'=u'|\mathbf{u}=u).\vspace{5pt}
\end{gathered}
\end{equation}
Each cross density $p(x, x')$ and $p(u, u')$ induces a density kernel function.
\label{definition_density}
\end{definition}

Figure~\ref{definition_density} illustrates this variational procedure and answers why this structure is essential for characterizing statistical dependence. The internal variables are created in two steps. First, a variational system with three variables $\mathbf{x}\rightarrow\mathbf{u}\rightarrow\mathbf{x}'$ or $\mathbf{u}\rightarrow\mathbf{x}'\rightarrow\mathbf{u}'$ is assumed, starting with any one of the r.p. $\mathbf{x}$ or $\mathbf{u}$. The second step applies marginalization \eqref{cross_densities} to obtain the joint density $p(x, x')$ or $p(u, u')$. The blue and red regions in Figure~\ref{definition_density} illustrate in the space of internal variables the domain where statistical independence with the marginalized variable exists. The cross density is a canonical quantity defined for the proposed statistical mapping that involves computation in both directions. Suppose $\mathbf{x}$ and $\mathbf{u}$ are strictly dependent. Given a realization $x$ of $\mathbf{x}$, if there is a one-to-one mapping between $x$ and another realization $u$ of $\mathbf{u}$, there will also be a one-to-one mapping between $x$ and a realization $x'$ of the variable $\mathbf{x}'$. In this case, the cross density $p(x, x')$ is nonzero only when $x' = x$. On the other hand, \textcolor{black}{if $\mathbf{x}$ and $\mathbf{u}$ are not strictly dependent}, there will not be only one realization $x'$ but a region with a positive cross density. In the \textcolor{black}{extreme} case when $\mathbf{x}$ and $\mathbf{u}$ are statistically independent, the blue and red regions should coincide with the domain of $\mathbf{x}$ and $\mathbf{u}$. \vspace{3pt}

\noindent \textbf{Self-evaluating internal variables:} \textcolor{black}{Comparing to GP, the bidirectional recursion no longer assumes stationarity or Gaussian priors. Instead, the internal variables $\mathbf{x}'$ and $\mathbf{u}'$ follow a special self-evaluating principle, sharing the same marginal as the external variables $\mathbf{x}$ and $\mathbf{u}$.} Moreover, after the marginalization, the internal variable $\mathbf{x}'$ depends only on $\mathbf{x}$ with a joint density $p(x, x')$, and $\mathbf{u}'$ depends only on $\mathbf{u}$ with $p(u, u')$. This procedure creates two systems that are self-evaluating, with each system defined only on one marginal of the joint space. This characteristic brings tremendous advantages because the interactions between individual systems provide an optimization framework to evaluate statistical dependence.

\vspace{3pt}
\noindent \textbf{Relations to the ACE:} The ACE implicitly implements the bidirectional recursion with realizations through recursive updates of two functions. As ACE converges, the two functions reach an equilibrium, representing an agreement between two self-evaluating systems. We should not forget that $p(x, x')$ and $p(u, u')$ may be joint densities of high-dimensional r.p., impossible to approximate with only two functions. The proposed approach instead seeks to approximate the exact cross density to characterize each of the internal variables with a kernel function, named the cross density kernel function (CDKF). As we will show, the equilibrium of the ACE expresses only two eigenfunctions of the CDKF. The definitions of the two CDKF $K(x, x')$ and $K(u, u')$ are presented in Theorem~\ref{lemma_2_kernel}.

\begin{theorem}
(CDKF definition) The functions $K(x, x') = \frac{p(x, x')}{p^{\frac{1}{2}}(x) p^{\frac{1}{2}}(x')}$ and $K(u, u') = \frac{p(u, u')}{p^{\frac{1}{2}}(u) p^{\frac{1}{2}}(u')}$ are positive-definite kernel functions.
\label{lemma_2_kernel}
\end{theorem}
\begin{proof}
We first show that $K(x, x')$ is symmetrical. \textcolor{black}
The denominator of $K(x, x')$ is symmetrical because $\mathbf{x}'$ is self-evaluating and has the same marginal density as $\mathbf{x}$}, i.e.
\begin{equation}
\begin{aligned}
p(\mathbf{x}' = x') = \int_{\mathcal{U}} p(\mathbf{x}' = x' |\mathbf{u}=u) p(\mathbf{u}=u)du = \int_{\mathcal{U}} p(\mathbf{x} = x' |\mathbf{u}=u) p(\mathbf{u}=u)du = p(\mathbf{x} = x').
\end{aligned}
\end{equation}
\vspace{-5pt}The numerator of $K(x, x')$ is also symmetrical because
\begin{equation}
\begin{aligned}
p(x, x') = p(\mathbf{x} = x, \mathbf{x}' = x') &= \int_\mathcal{U}p(\mathbf{x}' =x'|\mathbf{u} = u)p(\mathbf{u} = u|\mathbf{x} = x)p(\mathbf{x} = x) du\\
& = \int_\mathcal{U}p(\mathbf{x} =x'|\mathbf{u} = u)p(\mathbf{u} = u|\mathbf{x}' = x)p(\mathbf{x}' = x) du \\
&= p(\mathbf{x} = x', \mathbf{x}' = x) = p(x', x).
\end{aligned}
\end{equation}
Therefore $K(x, x')$ is symmetrical. Next, we show that $K(x, x')$ is positive-definite. Suppose we are given a sequence of possible realizations $\{x_i\}_{i=1}^N$ and a set of scalars $\{\alpha_i\}_{i=1}^N$, we have
\begin{equation}
\begin{aligned}
\sum_{i=1}^N \sum_{j=1}^N \alpha_i \alpha_j K(x_i, x_j) &= \sum_{i=1}^N \sum_{j=1}^N \alpha_i \alpha_j \frac{\int_{\mathcal{U}} p(x_j|u) p(u|x_i) p(x_i) du}{p^{\frac{1}{2}}(x_i) p^{\frac{1}{2}}(x_j)} \\
&= \int_{\mathcal{U}} \sum_{i=1}^N \sum_{j=1}^N \frac{  \alpha_i \alpha_j  p(x_i, u) p(x_j, u)}{p^{\frac{1}{2}}(x_i) p^{\frac{1}{2}}(x_j) p(u)} du \\
&= \int_{\mathcal{U}} \left( \sum_{i=1}^N \frac{\alpha_i p(x_i, u)}{p^{\frac{1}{2}}(x_i)}  \right)^2 \frac{1}{p(u)} du \geq 0.
\end{aligned}
\end{equation}
\textcolor{black}{Therefore $K(x, x')$ is also positive-definite. Since $K(x, x')$ is both symmetrical and positive-definite, the function is a positive-definite kernel function.} Following the same argument, we can also show that $K(u, u') = \frac{p(u, u')}{p^{\frac{1}{2}}(u) p^{\frac{1}{2}}(u')}$ is also a positive-definite kernel.
\end{proof}
Since an obvious descriptor of a kernel function is its eigenfunction expansion, we \textcolor{black}{immediately} apply Mercer's theorem to the linear operator of the CDKF kernel and show that its eigenvalues provide a numerical quantification of statistical dependence.

\begin{corollary}
\textcolor{black}{By applying} Mercer's theorem, the linear operator
\begin{equation}
\begin{aligned}
T \varphi (x) = \int_{\mathcal{X}}K(x, x') \varphi (x') dx'
\end{aligned}
\label{linear_operator}
\end{equation}
induces an orthonormal eigenfunction set $\{\phi_i \}_{i=1}^\infty$ and eigenvalues $\{ \lambda_i \}_{i=1}^\infty$ such that \begin{equation}
\begin{gathered}
\int_{\mathcal{X}}\phi^2_i(x) dx = 1,\; T\phi_i(x) = \lambda_i \phi_i,\; i=1,2,\cdots, \infty,\\
K(x, x') = \sum_{i=1}^\infty \lambda_i \phi_i(x) \phi_i(x').
\end{gathered}
\end{equation}
The CDKF eigenspectrum $\{ \lambda_i \}_{i=1}^\infty$ is a measure of statistical dependence. The kernel $K(u, u')$ also induces an orthonormal eigenfunction set $\{\psi_i \}_{i=1}^\infty$. Moreover, the two CDKF $K(x, x')$ and $K(u, u')$ share the same eigenvalues $\{ \lambda_i \}_{i=1}^\infty$.
\label{eiproperty}
\end{corollary}
In this paper, without loss of generality, we assume that the eigenvalues are ordered, as $\lambda_1 \geq \lambda_2 \geq \cdots $. The eigenfunctions $\phi_1, \phi_2, \cdots$ and $\psi_1, \psi_2, \cdots$ are also ordered by their eigenvalues. The most important property of CDKF is that its eigenspectrum $\{ \lambda_i \}_{i=1}^\infty$ is a proper statistical dependence measure.
\begin{lemma}
For any given r.p., the largest eigenvalue of CDKF satisfies $\lambda_1 = 1$. All the other eigenvalues of CDKF are zero $\lambda_2 = \lambda_3 = \cdots = 0$ if and only if two r.p. are strictly independent.

\label{dependence}
\end{lemma}
\begin{proof}
We first show $\lambda_1 = 1$ for any given r.p.. The strategy is to show there is an eigenfunction $\phi(x) = p^{\frac{1}{2}}(x)$ with an eigenvalue $\lambda_1 = 1$. By choosing $\phi_1(x) = p^{\frac{1}{2}}(x)$, it follows easily $\int \phi_1^2(x) dx = \int p(x) dx = 1$, thus the norm of the function is $1$. Next we apply the linear operator to $\phi_1(x) = p^{\frac{1}{2}}(x)$, which gives
\begin{equation}
\begin{gathered}
T\phi_1(x) = \int_{\mathcal{X}} \frac{p(x, x')}{p^{\frac{1}{2}}(x)p^{\frac{1}{2}}(x')} p^{\frac{1}{2}}(x') dx' = p^{\frac{1}{2}}(x) = \phi_1(x).
\end{gathered}
\end{equation}
Thus $\phi_1(x) = p^{\frac{1}{2}}(x)$ is an eigenfunction with the maximal eigenvalue $\lambda_1 = 1$. Now suppose the two r.p. are strictly independent, \textcolor{black}{the CDKF will have the form} $K(x, x') = p^{\frac{1}{2}}(x)p^{\frac{1}{2}}(x') = \phi_1(x) \phi_1(x')$. This implies that $K(x, x')$ only has one positive eigenvalue and one eigenfunction. Conversely, if $K(x, x')$ can be written in the form $K(x, x') = \phi_1(x) \phi_1(x') = p^{\frac{1}{2}}(x)p^{\frac{1}{2}}(x')$, it can be shown that
\begin{equation}
\begin{gathered}
K(x, x') = \frac{p(x, x')}{p^{\frac{1}{2}}(x)p^{\frac{1}{2}}(x')} = p^{\frac{1}{2}}(x)p^{\frac{1}{2}}(x') \Rightarrow p(x, x') = p(x)p(x'),
\end{gathered}
\end{equation}
which is satisfied only when two r.p. $\mathbf{x}$ and $\mathbf{u}$ are strictly independent.
\end{proof}
Lemma~\ref{dependence} shows that the eigenvalues of CDKF satisfy one of the important postulates for any measure of statistical dependence on two r.v.~(\cite{renyi1959measures}), which states that the value of the measure must vanish when two r.v. are strictly independent. \textcolor{black}{Lemma~\ref{dependence} is an extension of this postulate} to characterize multidimensional statistical independence or dependence by employing the set of eigenvalues from the full eigenspectrum. Now with Lemma~\ref{dependence} as the theoretical foundation, we can claim Definition~\ref{definition_eigenspectrum} for the eigenspectrum.
\begin{definition}
\textcolor{black}{The eigenspectrum $\{ \lambda_i \}_{i=1}^\infty$ from each of the two CDKF $K(x, x')$ and $K(u, u')$ is a novel multivariate statistical dependence measure. If a bandwidth $K$ is specified, the leading $K$ eigenvalues $\{ \lambda_i \}_{i=1}^K$ are taken as the statistical dependence measure, which also becomes the learning objective. }
\label{definition_eigenspectrum}
\end{definition}

The CDKF and its eigenspectrum have tremendous advantages compared to conventional methods and cost functions employed in the neural network community for statistical dependence estimation. \textcolor{black}{Many measures or cost functions are either one-dimensional or with a dimension limited to the number of classes or the target signal dimension.} Instead, CDKF quantifies all eigenvalues over the spatial or time structure in the functional space. For cost functions implementing proxies of statistical dependence, it is impossible to fully characterize the error signal \textit{pdf}. CDKF also extends the ACE. Instead of finding only two scalar-valued functions, CDKF characterizes the full joint densities \textcolor{black}{by solving the eigenproblem}.
\subsection{Learning framework with r.p. realizations}
 \label{design_principle_rp}
The remaining problem is that the integration of the linear operator is with respect to the Lebesgue measure, and normally we are given realizations of r.p.. The standard variational technique of changing the probability measures has been used in the literature of the GP (\cite{williams2006gaussian}), \textcolor{black}{which can be applied here as follows.}
\begin{lemma}
Suppose $\{ \lambda_i \}_{i=1}^\infty$ are \textcolor{black}{the eigenvalues of CDKF, there exists a set of functions} $\{\widehat{\phi_i}(x)\}_{i=1}^\infty$ such that
\begin{equation}
\begin{gathered}
\int_{\mathcal{X}} \widehat{\phi_i}(x) \widehat{\phi_j}(x) p(x) dx = \begin{cases} 1, \;i=j& \\ 0, \;i\neq j \end{cases}\hspace{-8pt},\;\; i, j=1,2,\cdots, \infty,
\end{gathered}
\label{Orthonormal}
\end{equation}
\begin{equation}
\begin{gathered}
\mathbb{E}_{\mathbf{x}'}[\widehat{\phi_i} (\mathbf{x}')|x] = \int_{\mathcal{X}} p(x'|x) \widehat{\phi_i} (x') dx' = \lambda_i \widehat{\phi_i} (x)\;\; a.e. p(x), \;\; i=1,2,\cdots, \infty.\\
\end{gathered}
\label{conditional_mean}
\end{equation}
Equation ~\eqref{Orthonormal} is named the orthonormal condition and~\eqref{conditional_mean} the equilibrium condition. Moreover, the set of functions $\{\widehat{\phi_i}(x)\}_{i=1}^\infty$ approximates a special probabilistic quantity which is called the cross density ratio (CDR)
\begin{equation}
\begin{gathered}
 \frac{p(x, x')}{p(x) p(x')} = \sum_{i=1}^\infty \lambda_i \widehat{\phi_i}(x) \widehat{\phi_i}(x).
\end{gathered}
\label{final_equation}
\end{equation}
The CDR in Equation~\eqref{final_equation} is denoted as $\rho (x, x'):=\frac{p(x, x')}{p(x) p(x')}$. Similarly, there exists a set of functions $\{\widehat{\psi_i}(u)\}_{i=1}^\infty$ that satisfy the two conditions for the r.p. $\mathbf{u}$, and approximate the CDR \textcolor{black}{$\rho(u, u') := \frac{p(u, u')}{p(u) p(u')} = \sum_{i=1}^\infty \lambda_i \widehat{\psi_i}(u) \widehat{\psi_i}(u)$.}
\label{theorem1}
\end{lemma}
As this variational technique is a standard approach, the details of the proof appear in Appendix~\ref{CDR_appendix}. Lemma~\ref{theorem1} also states that the variational eigenproblem shares the same eigenvalues as the CDKF. From the statistical dependence measure perspective, the decomposition differs from its original formulation only by $p^{-\frac{1}{2}}(x)$.  Most importantly, the equilibrium condition involves only empirical quantities and leads directly to an optimization problem from data using neural networks or other functional approximators.

The CDR (or CDKF) defines the cost function and the optimal criterion for the learning system. Whenever the CDR has a large value, it implies that there is a high level of statistical dependence in the space of the marginal. For supervised learning, $\rho(x, x')$ should display a large value iff $x$ and $x'$ are in the same class. For unsupervised learning,  $\rho(x, x')$ should have a large value iff $x$ and $x'$ are the same sample. We will also show how $\rho(x, x')$ is associated with the theory of factorial codes~(\cite{barlow1989finding}).

\section{Solving the variational eigenproblem with the functional maximal correlation algorithm}
\label{section_fmca}

Solving the variational eigenproblem requires implementing the bidirectional recursion with realizations and finding the eigenfunctions with universal function approximators. The ACE and the maximal correlation algorithm (MCA)~(\cite{hu2021mimo}) are two  \textcolor{black}{implementations closely} related to our proposed approach. ACE implements the bidirectional recursion to find an equilibrium through recursive updates. MCA uses two neural networks to optimize a correlation ratio cost function between a multidimensional nonlinear function of the input data and the scalar target or desired data for the purpose of regression and system identification.

Here we introduce a new algorithm based on the proposed variational system, named the functional maximal correlation (FMCA). \textcolor{black}{FMCA is also based on bidirectional recursive updates. The difference between FMCA versus previous attempts is to find an equilibrium over two functional families, using the autocorrelation function (ACF), the crosscorrelation function (CCF), and the Wiener filters concept. At the practical level, the implementation requires two neural networks, as shown in Figure~\ref{figure_fmca_pipeline}.} Given two r.p., \textcolor{black}{each neural network will take one r.p. as the input and produce a set of functions (i.e. approximate eigenfunctions). The goal for each neural network at each iteration is to predict the outputs of the other neural network until a good approximation of eigenvalues is found}. The most important theorem in this section is Theorem~\ref{main_theorem}, where we show that neural networks learned by FMCA \textcolor{black}{solve} the variational eigenproblem, approximate the CDR, and find the leading eigenvalues in the eigenspectrum. The strategy is to show that the trained neural networks will satisfy the \textcolor{black}{orthonormality} and the equilibrium requirements of Lemma~\ref{theorem1}.

\begin{figure}[t]
\captionsetup{type=figure}
\centering
{\includegraphics[width=1\textwidth]{{./used_figures/pipelinefigure_2}.pdf}}\\
\captionof{figure}{\small Diagram describing the architecture and the FMCA cost function to find the leading eigenvalues of the CDKF. The pipeline has three parts:  constructs the input data latent space and the reference r.p., optimizes the two neural networks, and obtains the functional decomposition by applying two normalization schemes. (a) Suppose two r.p. are given, FMCA is applied to quantify the statistical dependence between them. For the reference r.p., we will follow design rules introduced in Section~\ref{IMPLEMENTATIONS} for both supervised and unsupervised learning. (b) The optimization of the two neural networks follows a special bidirectional optimization scheme, where the output of each network is used as the target for the other in alternation. (c) After the networks are trained, we apply two normalization schemes to obtain the CDKF eigenvalues and the CDR decomposition. }
\label{figure_fmca_pipeline}
\vspace{-10pt}
\end{figure}

Suppose the data from the real world contains a high-dimensional r.p. with rich and complex functional structures, such as an image dataset. This r.p. is the primary r.p.. Correspondingly, a \textcolor{black}{reference} r.p. can also be available that has lower dimensions and less complexity for the purpose of achieving a specific task, such as the labels of the images for classification, a target r.p. for regression or simply uniform, unstructured random noise for variational applications. Therefore, depending upon the design of the reference r.p., various CDR forms are possible. The latent structure comprising the two neural network outputs can be considered a coding mechanism, which unifies supervised or unsupervised learning tasks. Among the proposed choices of the reference r.p., the unstructured random noise is an unsupervised mechanism that addresses the bi-dimensional property of an r.p.: the probability measure and its functional structure, and will lead to the implementation of factorial
\textcolor{black}{codes}.

\subsection{The maximal correlation is the second largest eigenvalue of CDKF}
\label{single_variate_representation}
Here the connection between R\'enyi's definition of maximal correlation (\cite{renyi1959measures}) and CDKF is discussed. \textcolor{black}{Lemma~\ref{lemma_10} shows that the maximal correlation is the second largest eigenvalue $\lambda_2$. The solutions of the maximal correlation are CDR's two eigenfunctions $\widehat{\phi_2}(x)$ and $\widehat{\psi_2}(u)$.} Suppose two r.v. $\mathbf{x}\in \mathcal{X}$ and $\mathbf{u}\in\mathcal{U}$ are given. The maximal correlation uses two scalar-valued functions $f:\mathcal{X}\rightarrow \mathbb{R}$ and $g:\mathcal{U}\rightarrow \mathbb{R}$ to find the supremum of the correlation coefficient $\textbf{cc}^*$ between $f(\mathbf{x})$ and $g(\mathbf{u})$, defined as
\begin{equation}
\begin{aligned}
\textbf{cc}(f(\mathbf{x}), g(\mathbf{u})) = \frac{\mathbb{E}[f(\mathbf{x}) g(\mathbf{u})] - \mathbb{E}[f(\mathbf{x})] \mathbb{E}[g(\mathbf{u})]}{\textbf{std}(f(\mathbf{x}))\, \textbf{std}(g(\mathbf{u}))},\;\;\textbf{cc}^*(\mathbf{x}, \mathbf{u}) = \sup_{f, g} \; \textbf{cc}(f(\mathbf{x}), g(\mathbf{u})),
\end{aligned}
\label{form1}
\end{equation}
where $\textbf{std}$ is the standard deviation of the random variables.

\vspace{5pt}
\noindent \textbf{Recursive updates of MCA:} MCA substitutes two functions $f$ and $g$ by two neural networks $f_\theta$ and $g_\omega$, then the supremum $\textbf{cc}^*$~\eqref{form1} is approximated by training the neural networks by gradient descent. Our previous work~(\cite{hu2021mimo}) follows this procedure. ACE implements the optimization problem~\eqref{form1} with bidirectional recursive updates, but the networks are restricted to have scalar outputs. For a fixed neural network $f_\theta$, the network $g_\omega$ is updated to match $f_\theta$ until \textcolor{black}{$g_\omega(u) \leftarrow \frac{\mathbb{E}_\mathbf{x}[f_\theta(\mathbf{x})|u]}{\text{std}(\mathbb{E}_\mathbf{x}[f_\theta(\mathbf{x})|\mathbf{u}])}$}, improving the overall correlation coefficient score. Conversely, for a fixed function $g_\omega$, the network $f_\theta$ will be updated to match $g_\omega$ until \textcolor{black}{$f_\theta(x) \leftarrow \frac{\mathbb{E}_\mathbf{u}[g_\omega(\mathbf{u})|x]}{\text{std}(\mathbb{E}_\mathbf{u}[g_\omega(\mathbf{u})|\mathbf{x}])}$.} \textcolor{black}{Each recursive update takes a step on the bidirectional recursion shown in Figure~\ref{figureprocess}. This relation is discussed mathematically in Appendix~\ref{reslationship_ace_recursion}.} The recursive updates continue until the two neural networks reach an equilibrium point.  The analysis of the variational eigenproblem answers the question of what this equilibrium point is. The lemma along with the important properties of the maximal correlation are presented below.

\begin{remark} The single-variate maximal correlation satisfies the following property:
\vspace{5pt}

(i) (Upper and lower bound) The maximal correlation is bounded as $0 \leq \mathbf{cc}^* \leq 1$;\vspace{5pt}

(ii) (Defining a metric) $r^* = \log(1-{\mathbf{cc}^*}^2)$ is a metric between $f^*$ and $g^*$;\vspace{5pt}

(iii) (Scale and shift-invariance)  Suppose $f^*$ and $g^*$ are the optimal solutions to the maximal correlation, then for any given scalars $a, b, c, d\in \mathbb{R}\backslash\{0\}$, the functions $a f^*+b$ and $c g^*+d$ are also the optimal solutions to the maximal correlation; \vspace{5pt}

(iv) (Vanishing condition) $\mathbf{cc}^*= 0$ if and only if two r.v. $\mathbf{x}$ and $\mathbf{u}$ are strictly independent;\vspace{5pt}

(v) (Invariance to bijective transformations) Applying any bijective functions to $\mathbf{x}$ or $\mathbf{u}$ does not change the value of the maximal correlation.
\label{remark}
\end{remark}

\begin{lemma}
Assuming $f^*$ and $g^*$ are the optimal solutions for the maximal correlation, then construct $\overline{f^*}(x) = \frac{f^*(x) - \mathbb{E}(f^*(\mathbf{x}))}{\mathbf{std}(f^*(\mathbf{x}))}$ and $\overline{g^*}(u) = \frac{g^*(u) - \mathbb{E}(g^*(\mathbf{u}))}{\mathbf{std}(g^*(\mathbf{u}))}$. The maximal correlation finds the second largest eigenvalue and the corresponding function of the variational eigenproblem, as $(\mathbf{cc}^*)^2 = \lambda_2$, $\overline{f^*}(x) = \widehat{\phi}_2(x)$ $a.e.p(x)$, and $\overline{g^*}(u) = \widehat{\psi}_2(u)$ $a.e.p(u)$.
\label{lemma_10}
\end{lemma}
The proof of Lemma~\ref{lemma_10} applies the Cauchy-Schwarz inequality and can be found in~\cite{renyi1959measures}. The details of the proof are included in the appendix. However, finding the maximal correlation is not adequate for many real-world applications, as it applies only to r.v. To illustrate this shortcoming, let us analyze the following example.
\begin{example}
Quantification of statistical dependence of an r.p. requires more than the second largest eigenvalue of the CDKF. Let $\mathbf{x}= \{\mathbf{x}_t\}_{t=1}^{T_1}$ and $\mathbf{u}=\{\mathbf{u}_t\}_{t=1}^{T_2}$ be given. Suppose there is a one-to-one mapping (strict dependence) between the first element of each process, \textcolor{black}{i.e.} $\mathbf{x}(1)$ and $\mathbf{u}(1)$, the coordinate maps $\pi_x(\mathbf{x}) = \mathbf{x}(1)$ and $\pi_u(\mathbf{u}) = \mathbf{u}(1)$ satisfy $\mathbf{cc}(\pi_x, \pi_u) = 1$. Regardless of the length $T_1$, $T_2$ and the functional structure in $\{\mathbf{x}_t\}_{t=2}^{T_1}$ and $\{\mathbf{u}_t\}_{t=2}^{T_2}$, the maximal correlation will always be $1$. Thus the maximal correlation overestimates the functional structure of the r.p..
\end{example}
Another justification is that the CDR approximation is not accurate with only a single scalar-valued function, regardless of its model complexity or representation power. As CDR further defines the joint density of the internal variables, it means that ACE or MCA only produces suboptimal internal variables. Notice that this shortcoming happens too often in neural networks where the cost is a scalar variable. This is exactly what FMCA improves.

\subsection{The functional maximal correlation algorithm}
\label{FMCA}

Instead of using only two scalar-valued functions, FMCA requires two \textcolor{black}{functional families} $\{f_i\}_{i=1}^K$ and $\{g_i\}_{i=1}^L$. We denote $\mathbf{f}=[f_1, \cdots, f_K]^\intercal:\mathcal{X} \rightarrow \mathbb{R}^{K}$ and $\mathbf{g}=[g_1, \cdots, g_L]^\intercal:\mathcal{U} \rightarrow \mathbb{R}^{L}$. Since the correlation coefficient score only evaluates the distance between two scalar-valued \textcolor{black}{functions}, a new quantity is needed that evaluates the overall score amongst the two families of functions match. The goal of the new score function is to find the maximal multivariate statistical dependence so it should include the top eigenvalues, where the actual number of top eigenvalues is a hyperparameter under the control of the user. In practice, the score function is estimated from the autocorrelation function (ACF) and the crosscorrelation function (CCF) between the outputs of the two networks. The following abbreviations will be used.
\begin{definition} Given $\mathbf{f}(\mathbf{x})$ and $\mathbf{g}(\mathbf{u})$ define:
\vspace{5pt}

(i) The ACF of $\mathbf{f}(\mathbf{x}):$ $\mathbf{R}_F := \mathbf{R}(\mathbf{f}(\mathbf{x})) = \mathbb{E}[\mathbf{f}(\mathbf{x})  \mathbf{f}^\intercal(\mathbf{x}) ] = \int_{\mathcal{X}} \mathbf{f}(x) \mathbf{f}^\intercal(x) p(x) dx;$ \vspace{5pt}

(ii) The ACF of $\mathbf{g}(\mathbf{u}):$ $\mathbf{R}_G := \mathbf{R}(\mathbf{g}(\mathbf{u})) = \mathbb{E}[\mathbf{g}(\mathbf{u})  \mathbf{g}^\intercal(\mathbf{u}) ]= \int_{\mathcal{U}} \mathbf{g}(u) \mathbf{g}^\intercal(u) p(u) du ;$ \vspace{5pt}

(iii) The CCF between $\mathbf{f}(x)$ and $\mathbf{g}(u):$
\begin{equation}
\begin{gathered}
\mathbf{P}_{FG} := \mathbf{P}(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u})) = \mathbb{E}[\mathbf{f}(\mathbf{x})  \mathbf{g}^\intercal(\mathbf{u}) ] = \int_{\mathcal{X}}\int_\mathcal{U} \mathbf{f}(x) \mathbf{g}^\intercal(u) p(x, u) dx du;
\end{gathered}
\end{equation}

(iv) The joint ACF combining $\mathbf{f}(x)$ and $\mathbf{g}(u):$
$\mathbf{R}_{FG} := \mathbf{R}(f(\mathbf{x}), g(\mathbf{u})) = \begin{bmatrix}\mathbf{R}_F & \mathbf{P}_{FG} \\
\mathbf{P}^\intercal_{FG} & \mathbf{R}_G
\end{bmatrix}.
$\vspace{5pt}

\label{definition_11}
\end{definition}
Based on Definition~\ref{definition_11}, the FMCA score function becomes
\begin{equation}
\begin{gathered}
r(\mathbf{f}(\mathbf{x}), g(\mathbf{u})) = \log\det \mathbf{R}(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u}))-\log \det \mathbf{R}(\mathbf{f}(\mathbf{x})) - \log\det \mathbf{R}(\mathbf{g}(\mathbf{u})).
\end{gathered}
\label{fmca_cost}
\end{equation}
FMCA seeks the infimum of the score function~\eqref{fmca_cost}. Such a problem can be written as
\begin{equation}
\begin{gathered}
r^*(\mathbf{x}, \mathbf{u}) = \inf_{\mathbf{f}, \mathbf{g}} r(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u})).
\end{gathered}
\label{fmca_maximal}
\end{equation}
For the purpose of regularization, a diagonal matrix with small values $\epsilon I$ is added to $\mathbf{R}(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u}))$, $\mathbf{R}(\mathbf{f}(\mathbf{x}))$ and $\mathbf{R}(\mathbf{g}(\mathbf{u}))$ whenever the log-determinant is computed. The maximal correlation~\eqref{form1} is a special case of FMCA when $K=L=2$. The two families of functions $\mathbf{f}$ and $\mathbf{g}$ can also be substituted by two multiple-output neural networks $\mathbf{f}_\theta$ and $\mathbf{g}_\omega$ for optimization.

\vspace{8pt}
\noindent \textbf{Recursive updates of FMCA:} The FMCA uses \textcolor{black}{bidirectional recursive updates} to obtain the infimum~\eqref{fmca_maximal}. This can be shown by applying the Schur complement to $\mathbf{R}_{FG}$ as
\begin{equation}
\begin{aligned}
\log\det \begin{bmatrix}{\mathbf{R}}_F & {\mathbf{P}}_{FG} \\
{\mathbf{P}}^\intercal_{FG} & {\mathbf{R}}_G
\end{bmatrix} &= \log\det {\mathbf{R}}_F + \log\det ({\mathbf{R}}_G  - {\mathbf{P}}^\intercal_{FG} {\mathbf{R}}_F^{-1} {\mathbf{P}}_{FG})\\
&= \log\det {\mathbf{R}}_G + \log\det ({\mathbf{R}}_F  - {\mathbf{P}}_{FG} {\mathbf{R}}_G^{-1} {\mathbf{P}}^\intercal_{FG}).
\end{aligned}
\label{equation_RFG}
\end{equation}
Next, these two forms are applied to the score function $r(\mathbf{f}(\mathbf{x}), g(\mathbf{u}))$, which yields
\begin{equation}
\begin{aligned}
r(\mathbf{f}(\mathbf{x}), g(\mathbf{u})) &= \log\det ({\mathbf{R}}_F  - {\mathbf{P}}_{FG} {\mathbf{R}}_G^{-1} {\mathbf{P}}^\intercal_{FG}) - \log\det {\mathbf{R}}_G \\
&= \log\det ({\mathbf{R}}_G  - {\mathbf{P}}^\intercal_{FG} {\mathbf{R}}_F^{-1} {\mathbf{P}}_{FG}) - \log\det {\mathbf{R}}_F.
\end{aligned}
\label{wiener_solution_schuml}
\end{equation}
For those familiar with Wiener filter theories, the new form~\eqref{wiener_solution_schuml} is the log determinant of the prediction error for two optimal Wiener filters, added with a regularizer that ensures invertible ACF. The two equalities in Equation~\eqref{wiener_solution_schuml}, \textcolor{black}{depending} upon how the Schur complement is applied, \textcolor{black}{represent} the two optimal Wiener filters with two directions.
The next subsection shows that during adaptation, the two \textcolor{black}{multivariate neural networks} reach the eigenfunctions and the leading eigenvalues of the variational eigenproblem, upon the application of a normalization scheme. \textcolor{black}{Still, one unsolved difficulty is the nonconvexity in the optimization. To implement the exact bidirectional recursion, the system has to attain the infimum of the score function at each update, which is obviously impossible for neural networks because of the enormous amount of local minima. Thus each FMCA update may not take the exact step over the bidirectional recursion. The error may also accumulate with recursive updates. How to apply convex learning machines to solve this issue is still under investigation. }

\subsection{FMCA solves the variational eigenproblem}
\label{CORE-THEOREM}
Our discussion will show that the functions found by FMCA satisfy the orthonormal condition~\eqref{Orthonormal} and the equilibrium requirements~\eqref{conditional_mean} in Lemma~\ref{theorem1}. The two neural networks will obtain the eigenspectrum as well as the approximation of CDR. Since CDR is defined by the joint densities of internal variables, upon convergence, FMCA will obtain the optimal internal variables for the bidirectional recursion. The essential part of the proof is to construct two normalization schemes that can be applied to two functionals, as expressed in  Equation~\eqref{normalize_inverse} in Lemma~\ref{lemma_invariance} and Equation~\eqref{second_normalization} in Theorem~\ref{main_theorem}.

There are two important properties of FMCA that constitute the proof: First, the score function evaluates the distance between two multivariate functions and it is upper bounded by $0$, which is shown in Lemma~\ref{lemma_less_than_0}; Second, the optimal solutions of FMCA is invariant to any affine transformation, which is shown in Lemma~\ref{lemma_invariance}. The second property extends the shift and scale-invariant properties of maximal correlation. After proving these two properties, we will present the main theorem of this section.
\begin{lemma}
(Score function's upper bound) For any given $\mathbf{f}$ and $\mathbf{g}$ such that $\mathbf{R}(\mathbf{f}(\mathbf{x}))$ and $\mathbf{R}(\mathbf{g}(\mathbf{u}))$ are invertible, the FMCA score function is upper bounded by $0$, i.e.
\begin{equation}
\begin{aligned}
r(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u})) \leq 0.
\end{aligned}
\end{equation}
And equivalently the optimal score function $r^*(\mathbf{x}, \mathbf{u}) \leq 0$.
\label{lemma_less_than_0}
\end{lemma}

\begin{proof}
This lemma follows the proof of Fischer's inequality. The proof starts by diagonalizing $\mathbf{R}_{FG}$ by normalizing $\mathbf{f}$ and $\mathbf{g}$. Assuming that $\mathbf{f}$ and $\mathbf{g}$ are such that $\mathbf{R}_F$ and $\mathbf{R}_G$ are full rank, represent the negative half power of $\mathbf{R}_F$ and $\mathbf{R}_G$ as $\mathbf{R}_F^{-\frac{1}{2}}$ and $\mathbf{R}_G^{-\frac{1}{2}}$, then diagonalize $\mathbf{R}_{FG}$ as
\begin{equation}
\begin{aligned}
\begin{bmatrix}\mathbf{R}^{-\frac{1}{2}}_F & \mathbf{0} \\
\mathbf{0} & \mathbf{R}^{-\frac{1}{2}}_G
\end{bmatrix} \begin{bmatrix}\mathbf{R}_F & \mathbf{P}_{FG} \\
\mathbf{P}^\intercal_{FG} & \mathbf{R}_G.
\end{bmatrix} \begin{bmatrix}\mathbf{R}^{-\frac{1}{2}}_F & \mathbf{0} \\
\mathbf{0} & \mathbf{R}^{-\frac{1}{2}}_G
\end{bmatrix} = \begin{bmatrix} \mathbf{I}_K & \mathbf{R}^{-\frac{1}{2}}_F \mathbf{P}_{FG} \mathbf{R}^{-\frac{1}{2}}_G \\
\mathbf{R}^{-\frac{1}{2}}_G \mathbf{P}^\intercal_{FG} \mathbf{R}^{-\frac{1}{2}}_F & \mathbf{I}_L
\end{bmatrix}:= \overline{\mathbf{R}_{FG}}.
\end{aligned}
\label{normalized}
\end{equation}
Since the diagonal of $\overline{\mathbf{R}_{FG}}$ are $1$, its trace $\text{Tr}(\overline{\mathbf{R}_{FG}}) = K+L$ is a constant for any feasible $\mathbf{f}$ and $\mathbf{g}$. Denoting the eigenvalues of $\overline{\mathbf{R}_{FG}}$ as $\{\lambda_i(\overline{\mathbf{R}_{FG}})\}_{i=1}^{K+L}$ yields
\begin{equation}
\begin{aligned}
\text{Tr}(\overline{\mathbf{R}_{FG}}) = \sum_{i=1}^{K+L} \lambda_i(\overline{\mathbf{R}_{FG}}) = K+L.
\end{aligned}
\label{determinant_inequality}
\end{equation}
The next step is to apply the inequality of arithmetic and geometric means to the eigenvalues $\{\lambda_i(\overline{\mathbf{R}_{FG}})\}_{i=1}^{K+L}$, as
\begin{equation}
\begin{aligned}
\det \overline{\mathbf{R}_{FG}} = \prod_{i=1}^{K+L} \lambda_i(\overline{\mathbf{R}_{FG}})
\leq \left(\frac{1}{K+L}\sum_{i=1}^{K+L} \lambda_i(\overline{\mathbf{R}_{FG}}) \right)^{K+L}= 1,
\end{aligned}
\label{inequality}
\end{equation}
Notice that if $K=L=1$, the inequality~\eqref{inequality} becomes Cauchy-Schwarz inequality. Observe that the cost function $r(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u}))$ can be written as the log determinant of $\overline{\mathbf{R}_{FG}}$, as
\begin{equation}
\begin{aligned}
\log\det \overline{\mathbf{R}_{FG}} &= \log\det \begin{bmatrix}\mathbf{R}_F & \mathbf{P}_{FG} \\
\mathbf{P}^\intercal_{FG} & \mathbf{R}_G
\end{bmatrix} + 2\log\det \begin{bmatrix}\mathbf{R}^{-\frac{1}{2}}_F & \mathbf{0} \\
\mathbf{0} & \mathbf{R}^{-\frac{1}{2}}_G
\end{bmatrix} \\
&= \log\det \mathbf{R}_{FG} - \log\det \mathbf{R}_F -  \log\det \mathbf{R}_G = r(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u})).
\end{aligned}
\label{equation_31}
\end{equation}
Applying ~\eqref{inequality} to~\eqref{equation_31} yields $r(\mathbf{f}(\mathbf{x}), \mathbf{g}(\mathbf{u}))\leq 0$. Since $r^*(\mathbf{x}, \mathbf{u})$ is the infimum, the optimal score $r^*(\mathbf{x}, \mathbf{u})\leq 0$.
\end{proof}
Lemma~\ref{lemma_less_than_0} not only shows the upper bound for the score function but also proposes an approach to normalize $\mathbf{R}_{FG}$ by the half-powered matrices $\mathbf{R}_F^{-\frac{1}{2}}$ and $\mathbf{R}_G^{-\frac{1}{2}}$. Furthermore, it can be shown that the upper bound is tight if and only if the two r.p. are strictly independent.

The next step is to show that the optimal FMCA solution is invariant to affine transformations (the second property). The notation $^*$ represents any optimal solution such that the infimum $r^*$ is attained. The corresponding ACFs and CCFs of $\mathbf{f}^*$ and $\mathbf{g}^*$ are denoted as $\mathbf{R}^*_F$, $\mathbf{R}^*_G$, $\mathbf{P}^*_{FG}$ and $\mathbf{R}^*_{FG}$. The lemma is presented as follows:
\begin{lemma}
(Invariance to affine transformations) Suppose two functions $\mathbf{f}^*$ and $\mathbf{g}^*$ are optimal solutions such that the infimum $r^*$ is attained. The following normalizations
\begin{equation}
\begin{aligned}
\overline{\mathbf{f}^*} = {\mathbf{R}^*_F}^{-\frac{1}{2}} \mathbf{f}^*, \;\; \overline{\mathbf{g}^*} = {\mathbf{R}^*_G}^{-\frac{1}{2}} \mathbf{g}^*.
\end{aligned}
\label{normalize_inverse}
\end{equation}
are also optimal solutions of FMCA, i.e.,
\begin{equation}
\begin{aligned}
r(\overline{\mathbf{f}^*}(\mathbf{x}), \overline{\mathbf{g}^*}(\mathbf{u})) = r({\mathbf{f}^*}(\mathbf{x}), {\mathbf{g}^*}(\mathbf{u})) =  r^*(\mathbf{x}, \mathbf{u}).
\end{aligned}
\label{maximalcorr}
\end{equation}
The ACFs and the CCFs of $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$ have the following forms
\begin{equation}
\begin{gathered}
\overline{\mathbf{P}^*_{FG}}:=\mathbf{P}(\overline{\mathbf{f}^*}, \overline{\mathbf{g}^*}) = {\mathbf{R}^*_F}^{-\frac{1}{2}} \mathbf{P}_{FG}^* {\mathbf{R}^*_F}^{-\frac{1}{2}}, \;\; \mathbf{R}(\overline{\mathbf{f}^*}) = \mathbf{I}_K, \;\; \mathbf{R}(\overline{\mathbf{g}^*}) = \mathbf{I}_L,\\
r({\mathbf{f}^*}(\mathbf{x}), {\mathbf{g}^*}(\mathbf{u})) = \log\det \overline{\mathbf{R}_{FG}^*}.
\end{gathered}
\label{r_star_rfg}
\end{equation}
Given any two invertible matrices $\mathbf{A} \in \mathbb{R}^{K \times K}$ and $\mathbf{B} \in \mathbb{R}^{L \times L}$, the functions $\mathbf{A} \overline{\mathbf{f}^*}$ and $\mathbf{B} \overline{\mathbf{g}^*}$ are also optimal solutions for the functional maximal correlation. Moreover, any two functionals from the sets
\begin{equation}
\begin{aligned}
\mathcal{F}^* &= \{\mathbf{A} \overline{\mathbf{f}^*} : \mathbf{A} \in \mathbb{R}^{K \times K} \hspace{-4pt}\text{ and full rank}\},\\
\mathcal{G}^* &= \{\mathbf{B} \overline{\mathbf{g}^*} : \mathbf{B} \in \mathbb{R}^{L \times L} \text{ and full rank}\}.
\end{aligned}
\label{solution_set}
\end{equation}
create a pair of optimal solutions for FMCA.
\label{lemma_invariance}
\end{lemma}
\begin{proof}
Since the functions are normalized, their ACFs are diagonalized as $\mathbf{R}(\overline{\mathbf{f}^*}(\mathbf{x})) = \mathbb{E}[\overline{\mathbf{f}^*} \;\overline{\mathbf{f}^*}^\intercal ] = {\mathbf{R}^*_F}^{-\frac{1}{2}}\mathbb{E}[\mathbf{f}^* {\mathbf{f}^*}^\intercal ]{\mathbf{R}^*_F}^{-\frac{1}{2}} = \mathbf{I}_K$. Next, we prove that applying the affine transformations with matrices $\mathbf{A}$ and $\mathbf{B}$ preserves the optimal FMCA solutions. Observe that the ACFs and the CCFs of $\mathbf{A} \overline{\mathbf{f}^*}$ and $\mathbf{B} \overline{\mathbf{g}^*}$ can be written as
\begin{equation}
\begin{gathered}
\mathbf{R}(\mathbf{A} \overline{\mathbf{f}^*}) = \mathbb{E}[\mathbf{A} \overline{\mathbf{f}^*}\;\overline{\mathbf{f}^*}^\intercal \mathbf{A}^\intercal] = \mathbf{A} \mathbf{I}_K \mathbf{A}^\intercal = \mathbf{A} \mathbf{A}^\intercal,\;\;
\mathbf{R}(\mathbf{B} \overline{\mathbf{g}^*})  = \mathbb{E}[\mathbf{B} \overline{\mathbf{g}^*}\;\overline{\mathbf{g}^*}^\intercal \mathbf{B}^\intercal] = \mathbf{B} \mathbf{I}_L \mathbf{B}^\intercal = \mathbf{B} \mathbf{B}^\intercal, \vspace{3pt}\\
\mathbf{P}(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) = \mathbb{E}[\mathbf{A} \overline{\mathbf{f}^*}\;\overline{\mathbf{g}^*}^\intercal \mathbf{B}^\intercal]= \mathbf{A} \overline{\mathbf{P}^*_{FG}} \mathbf{B}^\intercal, \vspace{3pt}\\
\mathbf{R}(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) = \begin{bmatrix}\mathbf{R}(\mathbf{A} \overline{\mathbf{f}^*}) & \mathbf{P}(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) \vspace{3pt}\\
\mathbf{P}^\intercal(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) & \mathbf{R}(\mathbf{B} \overline{\mathbf{g}^*})
\end{bmatrix} = \begin{bmatrix}\mathbf{A} & \mathbf{0} \\
\mathbf{0} & \mathbf{B}
\end{bmatrix} \overline{\mathbf{R}_{FG}^*} \begin{bmatrix}\mathbf{A}^\intercal & \mathbf{0} \\
\mathbf{0} & \mathbf{B}^\intercal
\end{bmatrix}.
\end{gathered}
\label{quantities_affine}
\end{equation}
Now we can compute $r(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*})$ with~\eqref{quantities_affine}, as
\begin{equation}
\begin{aligned}
r(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) &= \log\det \mathbf{R}(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) - \mathbf{R}(\mathbf{A} \overline{\mathbf{f}^*}) - \mathbf{R}(\mathbf{B} \overline{\mathbf{g}^*}) \\
&= \log\det \overline{\mathbf{R}_{FG}^*} + 2\log\det \mathbf{A} + 2\log\det\mathbf{B} - 2\log\det \mathbf{A} - 2\log\det \mathbf{B}\\
&= \log\det \overline{\mathbf{R}_{FG}^*} = r(\overline{\mathbf{f}^*}, \overline{\mathbf{g}^*}) = r^*.
\end{aligned}
\end{equation}
It can be seen that matrices $\mathbf{A}$ and $\mathbf{B}$ are canceled in the deviation, and $r(\mathbf{A} \overline{\mathbf{f}^*}, \mathbf{B} \overline{\mathbf{g}^*}) = r(\overline{\mathbf{f}^*}, \overline{\mathbf{g}^*}) = r^*$. Thus $\mathbf{A} \overline{\mathbf{f}^*}$ and $\mathbf{B} \overline{\mathbf{g}^*}$ are also the optimal solutions of FMCA for any invertible matrices $\mathbf{A}$ and $\mathbf{B}$.
\end{proof}
This property is important to better apply our previous work (\cite{hu2021mimo}) of exploiting the MCA invariance property to demonstrate that neural networks can be trained layerwise. First, the optimal solutions for $r^*$ form a solution set~\eqref{solution_set}. Any two solutions in the set~\eqref{solution_set} are equivalent, which greatly simplifies the optimization and the design of universal mappers. Second, Lemma~\ref{lemma_invariance} introduces the first normalization scheme~\eqref{normalize_inverse}, which indicates that any of the optimal solutions are equivalent to $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$ with diagonalized ACFs. Therefore, optimal solutions are specified by $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$, as all other solutions can be normalized with Equation~\eqref{normalize_inverse}.

Next, the main theorem of the section demonstrates that FMCA finds the leading eigenvalues and the decomposition of the CDR shown in Lemma~\ref{theorem1}. The initial claim is that the optimization reaches the equilibrium when the two \textcolor{black}{families} of functions satisfy the Wiener solution bidirectionally, as shown in Theorem~\ref{theorem_wiener}. By applying Theorem~\ref{theorem_wiener} over the bidirectional recursion, two fixed-point equations shown in Lemma~\ref{lemma_P} can be derived. Moreover, Lemma~\ref{lemma_P} also introduces the second normalization scheme~\eqref{second_normalization}. By applying both normalization schemes \eqref{normalize_inverse} and \eqref{second_normalization}, we prove in Theorem~\ref{main_theorem} that the optimal solutions of FMCA satisfy both the orthonormal condition~\eqref{Orthonormal} and the equilibrium condition~\eqref{conditional_mean}.
\begin{theorem}
(The FMCA equilibrium is a Wiener solution) The functions $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$ satisfy the Wiener solutions
\begin{equation}
\begin{aligned}
\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|x] &= \int_{\mathcal{U}} \overline{\mathbf{g}^*}(u) p(u|x) du = \overline{\mathbf{P}^*_{FG}}  \; \overline{\mathbf{f}^*}(x), \\
\mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|u] &= \int_{\mathcal{X}} \overline{\mathbf{f}^*}(x) p(x|u) du = \overline{\mathbf{P}^*_{FG}}^\intercal  \; \overline{\mathbf{g}^*}(u). \\
\end{aligned}
\label{P_fg}
\end{equation}
\label{theorem_wiener}
\end{theorem}

\begin{lemma}
For simplicity, denote $\overline{\mathbf{P}^*} := \overline{\mathbf{P}^*_{FG}}$. Given the bidirectional recursion in Definition~\ref{definition_density}, the functionals $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$ have the following relations
\begin{equation}
\begin{gathered}
(\overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal})^{-1} \int_{\mathcal{X}} \overline{{f}^*_i}(x) P(x'|x) dx' = \overline{{f}^*_i}(x), \;\; i=1,2,\cdots, K,\\
(\overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*})^{-1} \int_{\mathcal{U}} \overline{{g}^*_i}(x) P(u'|u) du' = \overline{{g}^*_i}(u), \;\; i=1,2,\cdots, L.
\end{gathered}
\end{equation}
\label{lemma_P}
\end{lemma}

\begin{theorem} (Main theorem) Suppose $K\leq L$, the functional maximal correlation finds the leading $K$ eigenvalues $\{\lambda_i\}_{i=1}^K$, the eigenfunctions $\{\widehat{\phi_i}\}_{i=1}^K$, and the eigenfunctions $\{\widehat{\psi_i}\}_{i=1}^K$ of CDR. Suppose the eigendecomposition of $\overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal}$ and $\overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}$ are given by
\begin{equation}
\begin{gathered}
\overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}= \mathbf{Q}_G^* \mathbf{\Sigma}^* {\mathbf{Q}_G^*}^\intercal,\;\; \overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal} = \mathbf{Q}_F^* \mathbf{\Sigma}^* {\mathbf{Q}_F^*}^\intercal,\;\; \mathbf{\Sigma}^* = \begin{bmatrix}
    \sigma_{1}^* &  \\
    & \ddots  \\
    & & \sigma_{K}^*
\end{bmatrix}.
\end{gathered}
\end{equation}
Construct the following functions $\widehat{\mathbf{f}^*}$ and $\widehat{\mathbf{g}^*}$ as
\begin{equation}
\begin{gathered}
\widehat{\mathbf{f}^*} = {\mathbf{Q}^*_F}^\intercal \overline{\mathbf{f}^*} = {\mathbf{Q}^*_F}^\intercal {\mathbf{R}^*_F}^{-\frac{1}{2}} \mathbf{f}^* :=  [\widehat{f^*_1}, \widehat{f^*_2}, \cdots, \widehat{f^*_K}]^\intercal,\\
\widehat{\mathbf{g}^*} = {\mathbf{Q}^*_G}^\intercal \overline{\mathbf{g}^*} = {\mathbf{Q}^*_G}^\intercal {\mathbf{R}^*_G}^{-\frac{1}{2}} \mathbf{g}^* :=  [\widehat{g^*_1}, \widehat{g^*_2}, \cdots, \widehat{g^*_L}]^\intercal.
\end{gathered}
\label{second_normalization}
\end{equation}
\textcolor{black}{It follows that $\widehat{\mathbf{f}^*}$ and $\widehat{\mathbf{g}^*}$ are eigenfunctions of CDR}, as
\begin{equation}
\begin{gathered}
\sigma^*_i = \lambda_i, \;\; i=1,2,\cdots, K \\
\widehat{f^*_i}(x) = \widehat{\phi_i}(x) \text{ a.e. } p(x), \;\; i=1,2,\cdots, K\\
\widehat{g^*_i}(u) = \widehat{\psi_i}(u) \text{ a.e. } p(u), \;\; i=1,2,\cdots, K
\end{gathered}
\label{equation_equl}
\end{equation}
\label{main_theorem}
\end{theorem}

\begin{proof}
Since $r^*(\mathbf{x}, \mathbf{u}) = \log\det \overline{\mathbf{R}_{FG}^*} = \log\det \begin{bmatrix} \mathbf{I}_K & \overline{\mathbf{P}^*} \\
\overline{\mathbf{P}^*}^\intercal  & \mathbf{I}_L
\end{bmatrix}$, the application of the Schur complement yields
\begin{equation}
\begin{gathered}
r^* = \log\det \overline{\mathbf{R}_{FG}^*} = \log \det(\mathbf{I}_K - \overline{\mathbf{P}^*}\; \overline{\mathbf{P}^*}^{\intercal}) = \log \det(\mathbf{I}_L - \overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}).
\end{gathered}
\label{Schur_complement}
\end{equation}
\textcolor{black}{Similar} to maximal correlation, let us first bound $r^*$ \eqref{Schur_complement} from below. Since $r^*$ is the infimum, the lower bound must be tight. To show this, rewrite the CCF $\overline{\mathbf{P}^*}$ in the form of the conditional operator, as
\begin{equation}
\begin{aligned}
\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|x] &= \int_{\mathcal{U}} \overline{\mathbf{g}^*}(u) p(u|x) du: \mathcal{X} \rightarrow \mathbb{R}^L, \;\;\mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|u] &= \int_{\mathcal{X}} \overline{\mathbf{f}^*}(x) p(x|u) du : \mathcal{U} \rightarrow \mathbb{R}^K,
\end{aligned}
\end{equation}
\vspace{-15pt}
\begin{equation}
\begin{aligned}
\overline{\mathbf{P}^*} = \int_{\mathcal{X}} \int_{\mathcal{U}} \overline{\mathbf{f}^*}(x)\overline{\mathbf{g}^*}(u) p(x, u)du dx = \int_{\mathcal{X}}  \overline{\mathbf{f}^*}(x) \mathbb{E}^\intercal_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|x] p(x)dx = \mathbf{P}(\overline{\mathbf{f}^*}(\mathbf{x}), \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]).
\end{aligned}
\label{CCF_conditional}
\end{equation}
Equation~\eqref{CCF_conditional} implies that the CCF between $\overline{\mathbf{f}^*}(\mathbf{x})$ and ${\mathbf{g}^*}(\mathbf{u})$ is equal to the CCF between $\overline{\mathbf{f}^*}(\mathbf{x})$ and the conditional mean $\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]$. This further implies that the ACF by combining $\overline{\mathbf{f}^*}(\mathbf{x})$ and $\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]$ has the following special form
\begin{equation}
\begin{aligned}
\mathbf{R}(\overline{\mathbf{f}^*}(\mathbf{x}), \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) =
\begin{bmatrix}
    \mathbf{I}_K &  \overline{\mathbf{P}^*}\\
    \overline{\mathbf{P}^*}^\intercal & \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}])
\end{bmatrix}.
\end{aligned}
\label{ACF_combined}
\end{equation}
Since $\mathbf{R}(\overline{\mathbf{f}^*}(\mathbf{x}), \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}])$ itself is an ACF, the properties of the ACF will be used for the proof. Since any ACF is Hermitian and positive semi-definite, the eigenvalues for any ACF are positive. The application of this property and the Schur complement to~\eqref{ACF_combined} yields
\begin{equation}
\begin{aligned}
\log\det \mathbf{R}(\overline{\mathbf{f}^*}(\mathbf{x}), \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) = \log\det \big( \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) - \overline{\mathbf{P}^*} \;\overline{\mathbf{P}^*}^\intercal \big) \geq 0.
\end{aligned}
\label{ACF_equation}
\end{equation}
Next, the concavity of the determinant leads to the following inequality for $r^*$
\begin{equation}
\begin{aligned}
\exp(r^*)= \det( \mathbf{I}_K - \overline{\mathbf{P}^*}\; \overline{\mathbf{P}^*}^{\intercal} ) &= \det \big(\mathbf{I}_K - \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) + \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) - \overline{\mathbf{P}^*} \;\overline{\mathbf{P}^*}^\intercal\big)\\
& \geq \det \big(\mathbf{I}_K - \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}])\big) + \det\big(\mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) - \overline{\mathbf{P}^*} \;\overline{\mathbf{P}^*}^\intercal\big)\\
&\geq \det \big(\mathbf{I}_K - \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}])\big).
\end{aligned}
\label{concavity}
\end{equation}
The second line of~\eqref{concavity} is known as the Minkowski determinant theorem. The proof is based on the concavity of the determinant and can be found in~\cite{marcus1992survey}. The third line of~\eqref{concavity} uses the lower bound $\det\big(\mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) - \overline{\mathbf{P}^*}\big)\geq 0$, proved in~\eqref{ACF_equation}.

Similarly to maximal correlation, since $r^*$ is obtained as the infimum, the inequality must become equality when the optimal condition is achieved, which implies that
\begin{equation}
\resizebox{1\linewidth}{!}{
$\begin{aligned}
\det(\mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) - \overline{\mathbf{P}^*} \;\overline{\mathbf{P}^*}^\intercal) = 0 \Rightarrow \mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) = \overline{\mathbf{P}^*} \; \overline{\mathbf{P}^*}^\intercal, \;\; \mathbf{R}( \mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|\mathbf{u}]) = \overline{\mathbf{P}^*}^\intercal \overline{\mathbf{P}^*}.
\end{aligned}$}
\label{equality_PP}
\end{equation}
Equation~\eqref{equality_PP} is well known as the Wiener solution that uses ACFs and CCFs to find the optimal linear projection. Since both $\overline{\mathbf{f}^*}$ and $\overline{\mathbf{g}^*}$ are normalized to be orthonormal, the Wiener solution is given by their CCF. It then follows that the infimum $r^*$ is obtained \textcolor{black}{iff} the Wiener solution~\eqref{Wiener_solution} is satisfied bidirectionally, as
\begin{equation}
\begin{aligned}
\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|x] = \int_{\mathcal{U}} \overline{\mathbf{g}^*}(u) p(u|x) du = \overline{\mathbf{P}^*}  \; \overline{\mathbf{f}^*}(x),\;\;\mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|u] = \int_{\mathcal{U}} \overline{\mathbf{f}^*}(u) p(x|u) du = \overline{\mathbf{P}^*}^\intercal \overline{\mathbf{g}^*}(u).
\end{aligned}
\label{Wiener_solution}
\end{equation}\vspace{-15pt}
\begin{equation}
\begin{aligned}
\Longleftrightarrow\;\; &\mathbf{R}( \mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}^*}(\mathbf{u})|\mathbf{x}]) = \int_{\mathcal{X}} \overline{\mathbf{P}^*}  \; \overline{\mathbf{f}^*}(x) \overline{\mathbf{f}^*}^\intercal (x) \overline{\mathbf{P}^*}^\intercal p(x) dx = \overline{\mathbf{P}^*} \mathbf{I}_K \overline{\mathbf{P}^*}^\intercal = \overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^\intercal, \\
 &\mathbf{R}( \mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|\mathbf{u}]) = \int_{\mathcal{U}} \overline{\mathbf{P}^*}^\intercal   \overline{\mathbf{g}^*}(u) \overline{\mathbf{g}^*}^\intercal (u) \overline{\mathbf{P}^*}\, p(u) du = \overline{\mathbf{P}^*}^\intercal \mathbf{I}_L \overline{\mathbf{P}^*} = \overline{\mathbf{P}^*}^\intercal \overline{\mathbf{P}^*}.
\end{aligned}
\label{satisfied}
\end{equation}

\end{proof}

\begin{proof}
Lemma~\ref{lemma_P} follows directly from Theorem~\ref{theorem_wiener} by applying Equation~\eqref{Wiener_solution} twice using the definition of the bidirectional recursion internal variables. For functions $\overline{\mathbf{f}^*}$, observe that
\begin{equation}
\begin{aligned}
\int_{\mathcal{X}} \overline{\mathbf{f}^*}(x') p(x'|x)dx' &= \int_{\mathcal{U}} \int_{\mathcal{X}} \overline{\mathbf{f}^*}(x') p(x'|u) dx' p(u|x) du \\
&= \int_{\mathcal{U}} \mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}^*}(\mathbf{x})|u] p(u|x) du \\
&= \overline{\mathbf{P}^*}\int_{\mathcal{U}} \overline{\mathbf{g}^*} (u) p(u|x) du \\
&= \overline{\mathbf{P}^*} \mathbb{E}_{\mathbf{u}} [\overline{\mathbf{g}^*} (\mathbf{u})|x] = \overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{f}^*}(x).
\end{aligned}
\end{equation}
The same proof applies to $\overline{\mathbf{g}^*}(u)$. Therefore we obtained Lemma~\ref{lemma_P}.
\end{proof}

\begin{proof}
The proof of Theorem~\ref{main_theorem} uses the conclusion from Lemma~\ref{theorem1} that states: if two sets of functions form a decomposition of the CDR, the functions must satisfy both the orthonormal condition~\eqref{Orthonormal} and the equilibrium condition~\eqref{conditional_mean}.

First, construct the functions $\widehat{\mathbf{f}^*}$ and $\widehat{\mathbf{g}^*}$ to satisfy these two conditions, which is the second normalization scheme in FMCA. From Lemma~\ref{lemma_P}, its eigendecomposition yields $\overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal} = \mathbf{Q}_F^* \mathbf{\Sigma}^* {\mathbf{Q}_F^*}^\intercal$ and $\overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}= \mathbf{Q}_G^* \mathbf{\Sigma}^* {\mathbf{Q}_G^*}^\intercal$, where the eigenvalues are denoted as $\sigma_{1}^*, \cdots, \sigma_{K}^*$. Using the property that the eigenvectors matrix $\mathbf{Q}_F^*$ is unitary, it follows that
\begin{equation}
\begin{aligned}
\int_{\mathcal{X}} \overline{\mathbf{f}^*}(x') p(x'|x)dx' = \mathbf{Q}_F^* \mathbf{\Sigma}^* {\mathbf{Q}_F^*}^\intercal \overline{\mathbf{f}^*}(x) \Rightarrow \int_{\mathcal{X}} {\mathbf{Q}_F^*}^\intercal \overline{\mathbf{f}^*}(x') p(x'|x)dx' = \mathbf{\Sigma}^* {\mathbf{Q}_F^*}^\intercal \overline{\mathbf{f}^*}(x).
\end{aligned}
\label{equation_eigen_QF}
\end{equation}
Now construct new functions $\widehat{\mathbf{f}^*} = {\mathbf{Q}^*_F}^\intercal \overline{\mathbf{f}^*} =  [\widehat{f^*_1}, \widehat{f^*_2}, \cdots, \widehat{f^*_K}]^\intercal$ and rewrite Equation~\eqref{equation_eigen_QF} as
\begin{equation}
\begin{aligned}
\int_{\mathcal{X}} \widehat{f^*_i}(x') p(x'|x)dx'  = \sigma_i \widehat{f^*_i}(x), \;\; i=1,2,\cdots, K,
\end{aligned}
\label{equlibrium_condition_f}
\end{equation}
which proves the equilibrium condition~\eqref{conditional_mean}. However, the functional family $\{\widehat{f^*_i}\}_{i=1}^K$ must satisfy the orthonormal condition $\int_{\mathcal{X}} \widehat{f^*_i}(x) \widehat{f^*_j}(x) p(x) dx = \begin{cases} 1, \;i=j& \\ 0, \;i\neq j \end{cases}$\hspace{-5pt} in~\eqref{Orthonormal}. To show this, observe that the ACF of $\widehat{\mathbf{f}^*}$ satisfies
\begin{equation}
\begin{aligned}
\int_{\mathcal{X}} \widehat{\mathbf{f}^*}(x) \widehat{\mathbf{f}^*}^\intercal(x) p(x)dx &= {\mathbf{Q}^*_F}^\intercal  {\mathbf{R}^*_F}^{-\frac{1}{2}} \int_{\mathcal{X}}\mathbf{f}^*(x) {\mathbf{f}^*}^\intercal(x) p(x)dx  {\mathbf{R}^*_F}^{-\frac{1}{2}} {\mathbf{Q}^*_F} \\
&= {\mathbf{Q}^*_F}^\intercal {\mathbf{R}^*_F}^{-\frac{1}{2}} {\mathbf{R}^*_F} {\mathbf{R}^*_F}^{-\frac{1}{2}} {\mathbf{Q}^*_F} = {\mathbf{Q}^*_F}^\intercal \mathbf{I_K} {\mathbf{Q}^*_F} = \mathbf{I_K}.
\end{aligned}
\label{orthonormal_main_theorem}
\end{equation}
Therefore the ACF of $\widehat{\mathbf{f}^*}(x)$ is an identity matrix, and the orthonormality is preserved through normalization. The direction for $\widehat{\mathbf{g}^*}(u)$ can be proved by symmetry. Now~\eqref{equlibrium_condition_f} implies the equilibrium condition~\eqref{conditional_mean}, and~\eqref{orthonormal_main_theorem} implies the orthonormal condition~\eqref{Orthonormal}. Therefore $\{\widehat{f^*_i}\}_{i=1}^K \subseteq \{\widehat{\phi_i}\}_{i=1}^\infty$ and $\{\sigma_i^*\}_{i=1}^K \subseteq \{\lambda_i\}_{i=1}^\infty$. To see that $\{\sigma_i^*\}_{i=1}^K$ are the leading eigenvalues, observe that
\begin{equation}
\begin{aligned}
r^*(\mathbf{x}, \mathbf{u}) = \det(\mathbf{I}_K - \overline{\mathbf{P}^*}\; \overline{\mathbf{P}^*}^{\intercal}) = \det(\mathbf{I}_L - \overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}) = \sum_{i=1}^K \log(1 - \sigma_i^*).
\end{aligned}
\end{equation}
Thus maximizing $r$ finds the leading eigenvalues, and the proof is complete.
\end{proof}
One simple corollary can be enunciated as follows.
\begin{corollary}
Given the bidirectional recursion defined on two r.p. $\mathbf{x}$ and $\mathbf{u}$, the two CDKF $K(x, x')$ and $K(u, u')$ share the same eigenspectrum $\{\lambda_i\}_{i=1}^\infty$.
\end{corollary}
This can be easily proved since $\overline{\mathbf{P}^*}\, \overline{\mathbf{P}^*}^{\intercal}$ and $\overline{\mathbf{P}^*}^{\intercal} \overline{\mathbf{P}^*}$ have the same eigenvalues. This corollary answers the question of how to choose the dimensionality $K$ and $L$. Since the number of the captured eigenvalues is limited by $\min(K, L)$, for ideal universal approximators, there will be no loss of generality assuming $K=L$. However, due to the limitation of \textcolor{black}{neural network size and} the adaptation, the total number of the captured eigenvalues may not reach the theoretical values. In practice, the dimensionality $K$ and $L$ should be chosen according to the complexity in the functional structure of \textcolor{black}{each} r.p..

\subsection{Implementing FMCA with neural networks}
\label{IMPLEMENTATIONS}

So far we have proved theoretical aspects of FMCA. This section focuses on the implementation of the learning system consisting of neural network mappers, the cost function, and the optimization methodology, shown in Figure~\ref{figure_fmca_pipeline}. Suppose we are given $N$ realizations of the r.p. $(x_n)_{n=1}^N$ and $(u_n)_{n=1}^N$ sampled from a joint density. Let two neural networks $\mathbf{f}_\theta:\mathcal{X} \rightarrow \mathbb{R}^K$ and $\mathbf{g}_\omega: \mathcal{U} \rightarrow \mathbb{R}^L$ be given, parameterized by $\theta$ and $\omega$ respectively. FMCA learns a set of eigenvalues $\{\sigma_i\}_{k=1}^K$, produces a decomposition to the functional space
\begin{equation}
\begin{aligned}
x_n(1), x_n(2),\cdots, x_n(T_1) &\Rightarrow (\widehat{\mathbf{f}_\theta})_1 (x_n), (\widehat{\mathbf{f}_\theta})_2(x_n), \cdots, (\widehat{\mathbf{f}_\theta})_K(x_n),\\
u_n(1), u_n(2),\cdots, u_n(T_2) &\Rightarrow (\widehat{\mathbf{g}_\omega})_1 (u_n), (\widehat{\mathbf{g}_\omega})_2(u_n), \cdots, (\widehat{\mathbf{g}_\omega})_L(u_n),
\end{aligned}
\label{functional_representation_fg}
\end{equation}
and approximates the CDR $\widetilde{\rho}(x, x') \approx \frac{p(x, x')}{p(x)p(x')}$ and $\widetilde{\rho}(u, u') \approx \frac{p(u, u')}{p(u)p(u')}$. Because realizations are employed, the ACFs and CCFs of neural networks are empirical estimations that will be denoted as $\widetilde{\mathbf{R}}_F$, $\widetilde{\mathbf{R}}_G$ and $\widetilde{\mathbf{P}}_{FG}$, written as
\begin{equation}
\resizebox{1\linewidth}{!}{
$\begin{gathered}
\widetilde{\mathbf{R}}_F = \frac{1}{N} \sum_{n=1}^N \mathbf{f}_\theta(x_n)  \mathbf{f}^\intercal_\theta(x_n) +\epsilon I_{K}, \;\; \widetilde{\mathbf{R}}_G = \frac{1}{N} \sum_{n=1}^N \mathbf{g}_\omega(x_n)  \mathbf{g}^\intercal_\omega(x_n) +\epsilon I_{L}, \;\; \widetilde{\mathbf{P}}_{FG} = \sum_{n=1}^N  \mathbf{f}_\theta(x_n) \mathbf{g}_\omega^\intercal(x_n).
\end{gathered}$}
\label{ACFCCFNN}
\end{equation}
For numerical stability, small diagonal matrices $\epsilon I_{K}$ and $\epsilon I_{L}$ are added to the ACFs. The estimation of the score function and the minimization problem are denoted as
\begin{equation}
\begin{gathered}
\widetilde{r}((\mathbf{f}_\theta(x_n), \mathbf{g}_\omega(u_n))_{n=1}^N) = \log\det \begin{bmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bmatrix}  - \log\det \widetilde{\mathbf{R}}_F - \log\det \widetilde{\mathbf{R}}_G,\\
\text{minimize}_{\theta, \omega} \;\; \widetilde{r}((\mathbf{f}_\theta(x_n), \mathbf{g}_\omega(u_n))_{n=1}^N).
\end{gathered}
\label{optimization_problem}
\end{equation}
Using the properties of the log determinant, the gradient of $\widetilde{r}$ has the following formula
\begin{equation}
\begin{gathered}
\frac{\partial \widetilde{r}}{\partial \theta} = \text{Tr} \left( \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}^{-1}  \frac{\partial \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}}{\partial \theta} \right)  - \text{Tr}\left( {\widetilde{\mathbf{R}}_F^{-1}} \frac{\partial \widetilde{\mathbf{R}}_F}{\partial \theta} \right), \\
\frac{\partial \widetilde{r}}{\partial \omega} = \text{Tr} \left( \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}^{-1}  \frac{\partial \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}}{\partial \omega} \right)  - \text{Tr}\left( {\widetilde{\mathbf{R}}_G^{-1}} \frac{\partial \widetilde{\mathbf{R}}_G}{\partial \omega} \right).\vspace{-5pt}
\label{gradient_form}
\end{gathered}
\end{equation}
During adaptation, one network is fixed, and the gradient steps~\eqref{gradient_form} \textcolor{black}{apply} only to the other network until $\widetilde{r}$ convergence. Then the process is switched to the other network and recursively repeated until the minimum from both directions is reached. A simplified approach is to take gradient updates~\eqref{gradient_form} of both networks at each iteration. Although this approach is easier to implement and converges faster, the recursive estimators are necessary to produce a stable and low-variance training trajectory. An optimizer such as Adam~\eqref{gradient_form} is adequate to estimate the ACFs and CCFs~\eqref{ACFCCFNN}. Observe that the gradients~\eqref{gradient_form} contain terms such as $\begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}^{-1}$, ${\widetilde{\mathbf{R}}_G^{-1}}$ and ${\widetilde{\mathbf{R}}_F^{-1}}$ that require taking the inverse. We always use the estimated quantities when we compute the inverse. The details of applying the recursive estimators can be found in Algorithm~\ref{algorithm_1}.

\begin{algorithm}[t]
\small
\caption{Pipeline of the functional maximal correlation algorithm (FMCA)}
\begin{algorithmic}[1]
\State Initialize neural networks $\mathbf{f}_\theta:\mathcal{X}\rightarrow \mathbb{R}^K$ and $\mathbf{g}_\omega: \mathcal{U} \rightarrow \mathbb{R}^L$. Initialize the ACF matrices $\mathbf{R}_F: K\times K$ and $\mathbf{R}_G: L\times L$ with identity matrices. Initialize $\mathbf{P}_{FG}$ with zeros.
\State Given the input random process $\{x_n(t):t=1,\cdots,T_1\}_{n=1}^N$; Construct the reference random process $\{u_n(t):t=1,\cdots,T_2\}_{n=1}^N$ following Section~\ref{reference_process}.
\While{the equilibrium is not reached}
\State $k \leftarrow k+1$
\State Sample $M$ realizations $\{(x_m. u_m)\}_{m=1}^M$ from the joint random process.
\State Calculate ${\mathbf{R}}_F = \frac{1}{M} \sum_{m=1}^M \mathbf{f}_\theta(x_m)  \mathbf{f}^\intercal_\theta(x_m) +\delta I_{K}$, \; ${\mathbf{R}}_G = \frac{1}{M} \sum_{m=1}^M \mathbf{g}_\theta(x_m)  \mathbf{g}^\intercal_\theta(x_m) +\delta I_{L}$,

${\mathbf{P}}_{FG} = \sum_{m=1}^M  \mathbf{f}_\theta(x_m) \mathbf{g}_\theta^\intercal(x_m)$.

\State Apply recursive estimators ${\mathbf{R}}_F^{(k+1)} \leftarrow \beta {\mathbf{R}}_F^{(k)} +  (1-\beta) \mathbf{R}_F$,  ${\mathbf{R}}_G^{(k+1)} \leftarrow \beta {\mathbf{R}}_G^{(k)} +  (1-\beta) \mathbf{R}_G$,

${\mathbf{P}}_{FG}^{(k+1)} \leftarrow \beta {\mathbf{P}}_{FG}^{(k)} +  (1-\beta) \mathbf{P}_{FG}$.

\State Calculate $\widetilde{\mathbf{R}}_F = {\mathbf{R}}_F/(1-\beta^k)$, $\widetilde{\mathbf{R}}_G = {\mathbf{R}}_G/(1-\beta^k)$, $\widetilde{\mathbf{P}}_{FG} = {\mathbf{P}}_{FG}/(1-\beta^k)$.

\State Calculate $\frac{\partial \widetilde{r}}{\partial \theta} = \text{Tr} ( \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}^{-1}  \frac{\partial \begin{bsmallmatrix}{\mathbf{R}}_F & {\mathbf{P}}_{FG} \\
{\mathbf{P}}^\intercal_{FG} & {\mathbf{R}}_G
\end{bsmallmatrix}}{\partial \theta} )  - \text{Tr}( {\widetilde{\mathbf{R}}_F^{-1}} \frac{\partial {\mathbf{R}}_F}{\partial \theta} )$

$\frac{\partial \widetilde{r}}{\partial \omega} = \text{Tr} ( \begin{bsmallmatrix}\widetilde{\mathbf{R}}_F & \widetilde{\mathbf{P}}_{FG} \\
\widetilde{\mathbf{P}}^\intercal_{FG} & \widetilde{\mathbf{R}}_G
\end{bsmallmatrix}^{-1}  \frac{\partial \begin{bsmallmatrix}{\mathbf{R}}_F & {\mathbf{P}}_{FG} \\
{\mathbf{P}}^\intercal_{FG} & {\mathbf{R}}_G
\end{bsmallmatrix}}{\partial \omega} )  - \text{Tr}( {\widetilde{\mathbf{R}}_G^{-1}} \frac{\partial {\mathbf{R}}_G}{\partial \omega} )$.

\State Gradient descent with an optimizer $\theta \leftarrow \text{optimizer}(\theta, \frac{\partial \widetilde{r}}{\partial \theta})$, $\omega \leftarrow \text{optimizer}(\omega, \frac{\partial \widetilde{r}}{\partial \omega})$.

\If{the eigenspectrum is needed}
\State Run eigendecomposition: ${\widetilde{\mathbf{P}}}_{FG}' =  \widetilde{\mathbf{R}}_F^{-\frac{1}{2}} \widetilde{\mathbf{R}}_{FG} \widetilde{\mathbf{R}}_G^{-\frac{1}{2}}, \;
{\widetilde{\mathbf{P}}}_{FG}'^{\intercal} {\widetilde{\mathbf{P}}}_{FG}' = \mathbf{Q}_G \mathbf{\Sigma} {\mathbf{Q}_G}^\intercal$,

\;\;\;\;\;${\widetilde{\mathbf{P}}}_{FG}' {\widetilde{\mathbf{P}}}_{FG}'^{\intercal} = \mathbf{Q}_F \mathbf{\Sigma} {\mathbf{Q}_F}^\intercal$.
\State Obtain the eigenvalues: $\mathbf{\Sigma} \Rightarrow \{\sigma_1, \cdots, \sigma_K\}$.
\State Normalize the functional by the ACF: $\overline{\mathbf{f}_\theta} = \widetilde{\mathbf{R}}_F^{-\frac{1}{2}} \mathbf{f}_\theta, \;\;  \overline{\mathbf{g}_\omega} = \widetilde{\mathbf{R}}_G^{-\frac{1}{2}} \mathbf{g}_\omega$.
\State Normalize the functional by the eigenvectors: $\widehat{\mathbf{f}_\theta} = \mathbf{Q}_F^T \overline{\mathbf{f}_\theta},\;\; \widehat{\mathbf{g}_\omega} = \mathbf{Q}_G^T \overline{\mathbf{g}_\omega}$.
\EndIf
\EndWhile
\State Given any $(x(t))_{t=1}^{T_1}$, obtain its functional representation $(x(t))_{t=1}^{T_1} \Rightarrow ((\widehat{\mathbf{f}_\theta})_k (x))_{i=1}^K.$
\State Given any $(x(t))_{t=1}^{T_1}$ and $(x'(t))_{t=1}^{T_1}$, obtain the density ratio $\widetilde{\rho}(x, x') = \widehat{\mathbf{f}_\theta}^\intercal(x) \mathbf{\Sigma} \widehat{\mathbf{f}_\theta}(x')$.
\end{algorithmic}
\label{algorithm_1}
\end{algorithm}%

After the networks are trained, the two normalization schemes discussed in Section~\ref{CORE-THEOREM} are applied to obtain the eigenvalues and the approximation of CDR. The first normalization step follows as
\begin{equation}
\begin{aligned}
\overline{\mathbf{f}_\theta} = \widetilde{\mathbf{R}}_F^{-\frac{1}{2}} \mathbf{f}_\theta, \;\;  \overline{\mathbf{g}_\omega} = \widetilde{\mathbf{R}}_G^{-\frac{1}{2}} \mathbf{g}_\omega.
\end{aligned}
\label{normalization_1}
\end{equation}
This step ensures that the obtained functions are orthonormal with respect to the probability measure. For the second step, denote the CCF between $\overline{\mathbf{f}_\theta}$ and $\overline{\mathbf{g}_\theta}$ as ${\widetilde{\mathbf{P}}}_{FG}'$. Then the eigendecomposition is performed to obtain the eigenvectors ${{\mathbf{Q}}}_F$, ${{\mathbf{Q}}}_G$ and the eigenspectrum ${{\mathbf{\Sigma}}}$, written as
\begin{equation}
\resizebox{1\linewidth}{!}{
$\begin{gathered}
{\widetilde{\mathbf{P}}}_{FG}' =  \sum_{n=1}^N  \overline{\mathbf{f}_\theta}(x_n) \overline{\mathbf{g}_\omega}^\intercal(x_n), \;\; {\widetilde{\mathbf{P}}}_{FG}' {\widetilde{\mathbf{P}}}_{FG}'^{\intercal} = \mathbf{Q}_F \mathbf{\Sigma} {\mathbf{Q}_F}^\intercal,\;\;
{\widetilde{\mathbf{P}}}_{FG}'^{\intercal} {\widetilde{\mathbf{P}}}_{FG}' = \mathbf{Q}_G \mathbf{\Sigma} {\mathbf{Q}_G}^\intercal,\;\;  \mathbf{\Sigma} = \begin{bmatrix}
    \sigma_{1} &  \\
    & \ddots  \\
    & & \sigma_{K}
\end{bmatrix}.
\end{gathered}$}
\label{eigenvalues_implementation}
\end{equation}
Next, employ $\mathbf{Q}_F$ and $\mathbf{Q}_G$ to obtain the eigenfunctions $\widehat{\mathbf{f}_\theta}$ and $\widehat{\mathbf{g}_\omega}$ of the CDR, as
\begin{equation}
\begin{gathered}
\widehat{\mathbf{f}_\theta} = \mathbf{Q}_F^T \overline{\mathbf{f}_\theta},\;\; \widehat{\mathbf{g}_\omega} = \mathbf{Q}_G^T \overline{\mathbf{g}_\omega}.
\end{gathered}
\label{normalization_2}
\end{equation}
The obtained eigenvalues $\{\sigma_1, \cdots, \sigma_K\}$~\eqref{eigenvalues_implementation} are the approximation of the leading eigenvalues of CDKF. The obtained functions $\widehat{\mathbf{f}_\theta} = \{\widehat{f_\theta^{(i)}}\}_{i=1}^K$ and $\widehat{\mathbf{g}_\omega} = \{\widehat{g_\omega^{(i)}}\}_{i=1}^L$ \textcolor{black}{are the eigenfunctions} of the r.p.. Finally, the CDR approximation becomes
\begin{equation}
\begin{gathered}
\widetilde{\rho}(x, x') = \sum_{i=1}^K \sigma_i \widehat{{f}_\theta^{(i)}}(x) \widehat{{f}_\theta^{(i)}}(x') \approx \frac{p(x, x')}{p(x) p(x')}, \;\; \widetilde{\rho}(u, u') = \sum_{i=1}^L \sigma_i \widehat{{g}_\omega^{(i)}}(u) \widehat{{g}_\omega^{(i)}}(u') \approx  \frac{p(u, u')}{p(u) p(u')}.
\end{gathered}
\end{equation}
The CDR is defined by the joint densities of the internal latent variables $\mathbf{x}'$ and $\mathbf{u}'$, and it is the perfect tool to achieve different machine learning tasks from the perspective of statistical dependence. Intuitively, the value of CDR implies the separability between two realizations $x$ and $x'$ of the r.p. $\mathbf{x}$, \textcolor{black}{under the condition that} the r.p. $\mathbf{u}$ is given. In the next section, we discuss how CDR can be utilized to design learning systems for various tasks. Algorithm~\ref{algorithm_1} presents the pseudo-code and summarizes this section.

\subsection{Designing coding mechanisms with FMCA}
\label{reference_process}

The learning system design starts from selecting the CDR $\widetilde{\rho}(x, x')$ according to the goal. We introduce four mechanisms for real-world datasets. For finite-horizon time series, the dataset contains realizations sampled from a 1D random process $\mathbf{x} = \{\mathbf{x}(t)\}_{t=1}^{T_1}$. For real-world images, the dataset contains the realizations from a 2D random process $\mathbf{x} = \{\mathbf{x}(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$. Real-world images are highly non-stationary. The marginal density of each given pixel and the joint densities between pairs or groups of pixels change drastically in 2D image space across the training set. This real-world r.p. is assigned as the primary r.p. and it is paired with a reference r.p. $\mathbf{u}$. The goal is to design the reference r.p. such that the CDR solves the desired task. Each form of the CDR is effectively a coding mechanism for the real-world r.p.. But notice that the goal does not change: Given any realization $x$, the CDR $\widetilde{\rho}(x, x')$ always provides an evaluation of the separability between $x$ and all realizations $x'$ in the training set. This property can obviously be used for both supervised and unsupervised learning. In the rest of the section, we propose four straightforward coding mechanisms, but many others may exist \textcolor{black}{for further investigations}. \vspace{5pt}

\noindent \textbf{Class coding:} Suppose labels of images are given for classification. A straightforward choice takes the labels as the reference r.p.. Each realization $u_n = [u_n(1), \cdots, u_n(T_2)]^\intercal$ is a one-hot vector sampled from a categorical distribution, where $T_2$ is the number of the classes. It is easy to verify that if $x$ and $x'$ are from two different classes, we have $\widetilde{\rho}(x, x') = 0$. \textcolor{black}{If} $x$ and $x'$ are from the same class and the class prior is $p$, it can be expected that $\widetilde{\rho}(x, x') = \frac{1}{p}$. \textcolor{black}{The CDR not only describes which class the sample is affiliated with but also reveals its prior probability in the dataset. Classes with fewer samples will have a higher CDR. }\vspace{5pt}

\noindent \textbf{Similarity coding:} Again, assume that the labels for the images are given, but instead of using the labels as the reference r.p. to train the network $\mathbf{g}_\omega$, at each iteration, for each image $\{x_n(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$, another image from the same class is sampled as the reference. It can be verified that the similarity coding achieves the same form of CDR as the class coding. But notice that this may require much fewer examples for good performance. \vspace{5pt}

\noindent \textcolor{black}{\textbf{Temporal recurrence coding:} FMCA can also be easily adapted for regression of time series. An infinite-horizon time series $\mathbf{x}(1), \mathbf{x}(2), \cdots$ as the input signal is paired with the target signal $\mathbf{d}(1), \mathbf{d}(2), \cdots$. The goal is to obtain a prediction $\mathbf{y}(t)$ of the target signal $\mathbf{d}(t)$ using the past samples of the input signal. This requires learning a filter with a user-specified order $K$, which takes $K$ samples of the input signal $\mathbf{x}(t-K+1), \;\mathbf{x}(t-K+2),\;\cdots,\; \mathbf{x}(t)$ and produces the prediction $\mathbf{y}(t)$. FMCA can follow this procedure and produce an iterated adaptive system. At each time step, the primary network $\mathbf{f}_\theta$ takes the past input signal of length $K$, and $\mathbf{g}_\omega$ take the target $\mathbf{d}(t)$ as the reference r.p. The two networks become two nonlinear filters of order $K$ and produce a special class of time-invariant eigenfunctions that generalize temporally to the entire time series. The FMCA eigenfunctions are a description of the dynamical system. To further obtain the regressor, we simply find the optimal Wiener solution between the CDR and the target signal.} \vspace{5pt}

\noindent \textbf{Maximal dependence coding:} When no label information is available for the dataset (unsupervised learning), the proved lower bound of the FMCA cost can be exploited. For each given image $\{x_n(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$ used as the primary r.p, the exact same image $\{x_n(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$ is selected as the reference r.p. sample. This procedure measures the statistical dependence between the image random process realization and itself, which obtains the maximal statistical dependence. In this case, the optimality condition of the CDR is that $\widetilde{\rho}(x, x')$ is positive (activated) if and only if $x=x'$. For any two different samples $x\neq x'$, the CDR is zero (deactivated), as $\widetilde{\rho}(x, x') = 0$. If the optimal CDR is reached, the learning system can capture and distinguish all realizations in the training set of the primary r.p., and the network outputs form the decomposition as a representation of the r.p.. However in experiments, we find that the maximal dependence coding is more inefficient and inaccurate than the factorial coding, which will be introduced next. \vspace{5pt}

\noindent \textbf{Factorial coding:} Factorial coding is an interesting case that uses uniform noise as the reference r.p.. Before the training starts, for each realization $\{x_n(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$, a series of codes from a uniform distribution are sampled as $\{u_n(t)\}_{t=1}^{T_2}$ such that the output of the reference network learns no spatial structure. Each realization is assigned a unique code of length $T_2$. Both the factorial coding and the maximal dependence coding have the same goal of capturing all statistical dependence in the realizations of the random process by setting the CDR $\widetilde{\rho}(x, x')$ to be positive if and only if $x=x'$, which captures the dissimilarities between all realizations.
\vspace{5pt}

\subsection{Importance of factorial codes for learning}
\label{Section_factorial_codes}
The concept of factorial codes was proposed in the neuroscience community (\cite{barlow1989finding, barlow1989unsupervised}) as an attempt to explain how the visual cortex and the brain can represent efficiently the world. The goal of factorial coding is to find invertible internal representations, which require that the occurrence of one code symbol is independent of all the other codes. Hence, the probability of occurrence of a particular input is simply the product of the probability of the corresponding codes. For this reason, the codes are called factorial. Factorial codes will be an excellent property for any machine learning algorithm because internal representations are indeed codes for the input data. An optimal system described by factorial codes does not follow the maximal mutual information or the smallest bandwidth principles but requires that a particular input image triggers a unique pattern of internal activation. Unfortunately, designing factorial codes is not easy. \textcolor{black}{\cite{barlow1989finding}} suggested minimizing the sum of the bit entropies of binary code components, which turns out to be a difficult task. \cite{schmidhuber1992learning} proposed a general unsupervised principle to extract distributed nonredundant internal representations of inputs by implementing two competitive forces: first, each internal unit should minimize its predictability given the input; second, for each representation unit there is a predictor that attempts to predict the response of the other units. This is a complex requirement for parameter tuning. Later attempts linked factorial coding to independent component analysis~(\cite{bach2002kernel}) and sparsity (\cite{olshausen1996emergence}), but there are no satisfactory solutions so far.

From a theoretical perspective, learning factorial codes is a natural goal for learning theory because they are associated with a generative model for r.p. e.g. image databases: starting with a uniform noise prior for each pixel, there is a ``law'' that can paint each pixel of each image with the appropriate color and intensity (the image generation model) to create a specific image (realization). Since learning is an inverse model of this generation process, the natural goal of a generative learning system is to produce outputs that are unstructured noise from any image in the training set, because in this case all the statistical dependence in the input imagery will be captured during learning.

\textcolor{black}{The CDKF and CDR} provide a perfect tool for implementing factorial codes. As the factorial code concept \textcolor{black}{require} that a particular latent pattern should appear if and only if a particular sample occurs, we translate this requirement into the condition that the CDR $\widetilde{\rho}(x, x')$ is positive if and only if $x=x'$. When $x$ and $x'$ are different samples, the CDR stays at $0$, to guarantee statistical independence of latent representations. If such an optimal CDR is learned by the system, it suffices to say that the learning system captures all dependence in the realizations. So, contradicting~\cite{schmidhuber1992learning}, there is a polynomial algorithm to learn factorial codes.

\section{Experiments for solving the eigenproblem with FMCA}
\label{validation_results}

The experimental results are divided in two sections: the current section tests theoretical predictions in some simple datasets. The focus is on the interpretation of the neural networks outputs as eigenvalues, which is also indirectly testing how close the neural network functional mapping approximates the eigenfunctions for the problem. We show that the training is very smooth and that network outputs actually approximate the eigenvalues. We also see that the eigenfunctions converge one at a time, as in adaptive PCA~(\cite{fancourt1998competitive}). This coupling is detrimental to fast adaptation and would benefit from a fixed-point update. Nevertheless, the bidirectional architecture is simpler to adapt than conventional neural networks because of the intrinsic orthogonality of the outputs, and should generalize better. We will verify that the eigenvalues and eigenfunctions reach the theoretical values. We first present the two estimated quantities of interest that will be checked in the examples:

\vspace{5pt}
\noindent\textbf{Learning eigenvalues:} FMCA solves the eigenproblem and approximates the following quantities from the trained neural nets.
\begin{equation}
\begin{gathered}
\widetilde{\sigma}_i = \lambda_i({\widetilde{\mathbf{P}}}_{FG}'^{\intercal} {\widetilde{\mathbf{P}}}_{FG}') = \lambda_i({\widetilde{\mathbf{P}}}_{FG}' {\widetilde{\mathbf{P}}}_{FG}'^{\intercal}) \approx {\lambda_i},\;\; i=1,2,\cdots, K\\
\widehat{f_\theta^{(i)}}(x) \approx \widehat{\phi_i}(x) \text{ a.e. } p(x), \;\; \widehat{g_\omega^{(i)}}(u) \approx \widehat{\psi_i}(u) \text{ a.e. } p(u),\;\; i=1,2,\cdots, K\\
\widetilde{\rho}(x, x') = \sum_{i=1}^K \widetilde{\sigma}_i \widehat{{f}_\theta^{(i)}}(x) \widehat{{f}_\theta^{(i)}}(x') \approx \frac{p(x, x')}{p(x) p(x')}, \;\; \widetilde{\rho}(u, u') = \sum_{i=1}^L \widetilde{\sigma}_i \widehat{{g}_\omega^{(i)}}(u) \widehat{{g}_\omega^{(i)}}(u') \approx  \frac{p(u, u')}{p(u) p(u')}.
\end{gathered}
\label{equivalence}
\end{equation}
Since $\widehat{f_\theta^{(i)}}(x)$ approximates $\widehat{\phi_i}(x)$ almost everywhere $p(x)$, for the continuous r.v. or r.p., we multiply $\widehat{f_\theta^{(i)}}(x)$ by $p(x)$ and visualize the approximation of the joint density
\begin{equation}
\begin{gathered}
\widehat{f_\theta^{(i)}}(x)p(x) \approx \widehat{\phi_i}(x)p(x) \text{ for } i=1,2,\cdots, K, \;\; \widetilde{\rho}(x, x')p(x)p(x') \approx p(x, x').
\end{gathered}
\label{almost_everywhere}
\end{equation}
For convenience, we show the comparison mainly between  $\widehat{f_\theta^{(i)}}$ and $\widehat{\phi_i}(x)$.
\vspace{5pt}

\noindent\textbf{Wiener solution:} The FMCA optimization reaches the equilibrium as the two neural networks satisfy the Wiener solution bidirectionally
\begin{equation}
\resizebox{1\linewidth}{!}{
$\begin{gathered}
\mathbb{E}_{\mathbf{u}}[\overline{\mathbf{g}_\omega}(\mathbf{u})|x] = \int_{\mathcal{U}} \overline{\mathbf{g}_\omega}(u) p(u|x) du = {\widetilde{\mathbf{P}}}_{FG}'  \; \overline{\mathbf{f}_\theta}(x),\;\;\mathbb{E}_{\mathbf{x}}[\overline{\mathbf{f}_\theta}(\mathbf{x})|u] = \int_{\mathcal{X}} \overline{\mathbf{f}_\theta}(x) p(x|u) du = {\widetilde{\mathbf{P}}}_{FG}'^{\intercal} \; \overline{\mathbf{g}_\omega}(u).
\end{gathered}$}
\label{Wiener_networks}
\end{equation}
In the synthetic experiments, the integration in~\eqref{Wiener_networks} is computed with the known $p(u|x)$ and $p(x|u)$ to calculate the fitting.
\vspace{5pt}

\subsection{ Finding the leading eigenvalues of CDKF with FMCA}
\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=0.247\textwidth]{{./used_figures/gauss_eig_save_figure_8_v2}.pdf}}\hspace{-5pt}%
\subfloat[]{\includegraphics[width=0.247\textwidth]{{./used_figures/gauss_eig_save_figure_13_v2}.pdf}}\hspace{-5pt}%
\subfloat[]{\includegraphics[width=0.257\textwidth]{{./used_figures/mkv_save_figure_4_v2}.pdf}}\hspace{-5pt}%
\subfloat[]{\includegraphics[width=0.257\textwidth]{{./used_figures/mkv_save_figure_1_v2}.pdf}}\\
\vspace{-33pt}\subfloat[]{\includegraphics[width=0.5\textwidth]{{./used_figures/mkv_spectrum_v3}.pdf}}%
\subfloat[]{\includegraphics[width=0.5\textwidth]{{./used_figures/gauss_spectrum_v3}.pdf}}\\
\vspace{-24pt}\captionof{figure}{\small Results of applying FMCA to GAUSS and MC. The first row shows the learning curves of the eigenvalues, where each color is a different eigenvalue. The theoretical values are marked with dotted lines. We show that the eigenvalues are sequentially learned from the largest to the smallest with roughly the same convergence speed. The second row shows the learned eigenvalues w.r.t. the level of statistical dependence, parameterized by $p_\alpha$ and $\rho_\alpha$. }
\label{learning_curve_eigenvalue}
\vspace{-10pt}
\end{figure}

In this section, FMCA is applied to the multiple discrete state space Markov chains (MC) where finding the eigenvalues is a standard eigenvectors expansion problem, and to the multiple multivariate Gaussian distributions (GAUSS), which can be computed in closed form using Hermite polynomials (\cite{williams2006gaussian}, \cite{zhu1997gaussian}). The hyperparameters for MC and GAUSS are respectively the diagonal of the transitional probability $p_\alpha$ and the determinant of the covariance matrix $\rho_\alpha$. As $p_\alpha$ or $\rho_\alpha$ approaches $1$, there will be a strict dependence in the joint space. As $p_\alpha$ or $\rho_\alpha$ approaches $0$, there will be strict independence in the joint space.
Figure~\ref{learning_curve_eigenvalue} shows the learning dynamics of eigenvalues and the comparison with their theoretical values, which shows that FMCA is both fast (around 5k iterations to converge), smooth-converging, and highly accurate. We also plot in Figure~\ref{learning_curve_eigenvalue} the eigenspectrum as a function of the hyperparameter $p_\alpha$ and $\rho_\alpha$, which further shows the accuracy of FMCA.

The eigenvalues converge sequentially from the largest eigenvalue to the smallest eigenvalue. FMCA convergence time depends upon the number of eigenvalues, which is similar to a Gram-Schmidt process. Moreover, the convergence speed for each eigenvalue is almost the same. This characteristic is useful for designing stopping and module training schemes since the user can choose how many eigenvalues are adequate according to the application. As expected from the theory, there is always a constant eigenvalue $\lambda_1 = 1$ captured by FMCA.

\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=.175\textwidth]{{./used_figures/gauss_eigenfunction_0}.pdf}}\hspace{-7.5pt}%
\subfloat[]{\includegraphics[width=.180\textwidth]{{./used_figures/gauss_eigenfunction_1}.pdf}}\hspace{-7.5pt}%
\subfloat[]{\includegraphics[width=.190\textwidth]{{./used_figures/gauss_eigenfunction_9}.pdf}}\hspace{-8pt}%
\subfloat[]{\includegraphics[width=.175\textwidth]{{./used_figures/mkv_eigenfunction_0_v2}.pdf}}\hspace{-7.5pt}%
\subfloat[]{\includegraphics[width=.175\textwidth]{{./used_figures/mkv_eigenfunction_1_v2}.pdf}}\hspace{-7.5pt}%
\subfloat[]{\includegraphics[width=.175\textwidth]{{./used_figures/mkv_eigenfunction_9_v2}.pdf}}\\
\vspace{-27pt}\subfloat[]{\includegraphics[width=.47\textwidth]{{./used_figures/ratio_plot_2}.pdf}}\hspace{0.001\textwidth}
\subfloat[]{\includegraphics[width=.49\textwidth]{{./used_figures/combine_plot_v2_2}.pdf}}
\vspace{-20pt}\captionof{figure}{\small Visualization of neural network outputs. The first row compares the eigenfunctions learned by neural nets (blue solid lines) and the theoretical solutions (red dotted line). The CDRs $\Tilde{\rho}(x, x')$ are also plotted for the Gauss and MC (bottom left), along with the squared distance between the model approximation and the theoretical values. We also show on the right bottom quadrant that the outputs of the network reach the Wiener solution of the Gauss and MC, by calculating the squared distance between  $\overline{\mathbf{f}_\theta}(x)$ and the conditional expectation $\mathbb{E}_\mathbf{u}[\overline{\mathbf{g}_\omega}(\mathbf{u})|x]$. This figure shows the high accuracy of the approximation using FMCA, for Theorem~\ref{theorem_wiener} and Theorem~\ref{main_theorem}.}
\label{eigenfunction}
\vspace{-10pt}
\end{figure}

Figure~\ref{eigenfunction} shows that the outputs of the neural networks after the two normalization schemes approximate and decompose the CDR $\Tilde{\rho}(x, x')$ with high accuracy. FMCA captures the Hermite polynomial for GAUSS, and the eigenvectors for MC. Figure~\ref{eigenfunction} also shows that the Wiener solution~\eqref{Wiener_networks} is satisfied for the training networks.

\subsection{CDKF eigenvalues as a novel multivariate statistical dependence measure}
\label{comparison_section_with_mine}

\begin{figure}
\captionsetup{type=table}
\captionsetup[subfigure]{labelformat=empty}
\centering
\resizebox{1\linewidth}{!}{
\begin{tabular}{rl|cccccccc|c|cc}
\toprule[3pt]
\multicolumn{2}{c|}{DATASET}&\multicolumn{9}{c|}{FMCA}&\multicolumn{2}{c}{MINE}\\
\cmidrule{3-13}
\multicolumn{2}{c|}{} &$\widetilde{\sigma}_2$&$\widetilde{\sigma}_3$&$\widetilde{\sigma}_4$&$\widetilde{\sigma}_5$&$\widetilde{\sigma}_6$&$\widetilde{\sigma}_6$&$\widetilde{\sigma}_7$&$\widetilde{\sigma}_8$&TOTAL&KL&RY\\
\midrule[2pt]
GAUSS&$\rho_\alpha=$ 0.214&0.041&0.001&0.0&0.0&0.0&0.0&0.0&0.0&0.002&0.024&0.024\\
&$\rho_\alpha=$ 0.357&0.124&0.015&0.001&0.0&0.0&0.0&0.0&0.0&0.007&0.07&0.067\\
&$\rho_\alpha=$ 0.5&0.245&0.057&0.011&0.002&0.0&0.0&0.0&0.0&0.018&0.141&0.14\\
&$\rho_\alpha=$ 0.643&0.411&0.16&0.054&0.013&0.002&0.0&0.0&0.0&0.039&0.259&0.258\\
&$\rho_\alpha=$ 0.786&0.608&0.36&0.196&0.097&0.039&0.013&0.003&0.001&0.088&0.479&0.443\\
&$\rho_\alpha=$ 0.928&0.856&0.722&0.586&0.469&0.349&0.253&0.171&0.107&0.293&0.992&0.770\\
\midrule[1.5pt]
MC&$p_\alpha=$ 0.143&0.579&0.499&0.494&0.447&0.402&0.371&0.315&0.196&0.231&0.931&0.722\\
&$p_\alpha=$ 0.357&0.773&0.742&0.731&0.683&0.68&0.678&0.62&0.523&0.492&1.46&0.963\\
&$p_\alpha=$ 0.571&0.896&0.878&0.873&0.866&0.837&0.802&0.776&0.758&0.804&1.81&1.05\\
&$p_\alpha=$ 0.786&0.962&0.953&0.949&0.944&0.937&0.931&0.909&0.894&1.2&2.06&1.12\\
&$p_\alpha=$ 1.0&0.997&0.996&0.996&0.996&0.995&0.995&0.995&0.994&2.41&2.24&1.15\\
\midrule[1.5pt]
GMM&EXP1&0.984&0.97&0.913&0.892&0.791&0.664&0.323&0.297&0.804&1.65&0.981\\
&EXP2&0.967&0.955&0.911&0.776&0.722&0.552&0.361&0.267&0.742&1.53&0.999\\
&EXP3&0.996&0.987&0.86&0.751&0.538&0.36&0.196&0.051&0.738&1.41&0.885\\
&EXP4&0.995&0.987&0.806&0.796&0.553&0.351&0.174&0.152&0.739&1.43&0.929\\
&EXP5&0.994&0.992&0.977&0.904&0.87&0.661&0.335&0.26&1.027
&1.73&1.01\\
\midrule[1.5pt]
OTHER& TM &0.432&0.425&0.344&0.171&0.042&0.029&0.007&0.004&0.091&0.704&0.463\\
&SPIRAL &0.773&0.73&0.698&0.662&0.598&0.576&0.564&0.549&0.804&1.120&0.668\\
&CHECK &0.93&0.003&0.002&0.001&0.001&0.001&0.0&0.0&0.129&0.59&0.342\\
\bottomrule[3pt]
\end{tabular}}
\caption{\small Comparison of the approximations by FMCA ($\tilde{\sigma}_1, \tilde{\sigma}_2, \cdots$) and MINE (KL and RY) as the baselines. The total power produced by FMCA is also shown in the table as TOTAL. Although both methods can be used to quantify the total power of statistical dependence, only FMCA provides a component-wise quantification of the dependence. A visualization of the results in polar plots is shown in Figure~\ref{CIRCLE_PLOT}.}
\label{table_giant}

\begin{minipage}[h]{1.0\linewidth}
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=.265\textwidth]{{./used_figures/GMM_1_eigencompare}.pdf}}\hspace{-11pt}%
\subfloat[]{\includegraphics[width=.265\textwidth]{{./used_figures/GMM_2_eigencompare}.pdf}}\hspace{-11pt}%
\subfloat[]{\includegraphics[width=.265\textwidth]{{./used_figures/SPIRAL_2_eigencompare}.pdf}}\hspace{-15pt}
\subfloat[]{\includegraphics[width=.265\textwidth]{{./used_figures/CHECK_2_eigencompare}.pdf}}
\vspace{-20pt}\captionof{figure}{\small Visualization of Table~\ref{table_giant} in polar plots. While the total power is very close for different datasets, the component-wise evaluations reveal complex hidden dependence structures. For example, the top eigenvalues of GMM-II are larger than GMM-I, but the last few eigenvalues of GMM-II are smaller than GMM-I. SPIRAL has more uniformly distributed eigenvalues that are more difficult to learn. CHECK only has one positive eigenvalue excluding the constant eigenvalue at $1$. }
\label{CIRCLE_PLOT}
\end{minipage}
\end{figure}

\begin{figure}
\captionsetup{type=table}
\captionsetup[subfigure]{labelformat=empty}
\centering
\resizebox{.9\linewidth}{!}{
\begin{tabular}{rl|cccccc}
\toprule[3pt]
\multicolumn{2}{c|}{DATASET}&\multicolumn{6}{c}{FMCA APPROX ERRORS}\\
\cmidrule{3-8}
\multicolumn{2}{c|}{} &$\widetilde{\sigma}_2$&$\widetilde{\sigma}_3$&$\widetilde{\sigma}_4$&$\widetilde{\sigma}_5$&$\widetilde{\sigma}_6$&$\widetilde{\sigma}_7$\\
\midrule[2pt]
GAUSS&$\rho_\alpha=$ 0.928&3.6$\times 10^{-5}$& 1.9$\times 10^{-5}$& 1.6$\times 10^{-5}$& 4.4$\times 10^{-5}$& 1.2$\times 10^{-5}$& 3.6$\times 10^{-5}$\\
\midrule[1.5pt]
MC&$p_\alpha=$ 0.786&9.0$\times 10^{-6}$& 1.6$\times 10^{-5}$& 1.6$\times 10^{-5}$& 1.6$\times 10^{-5}$& 2.5$\times 10^{-5}$& 2.5$\times 10^{-5}$\\
\midrule[1.5pt]
GMM&EXP1&4.0$\times 10^{-6}$& 9.0$\times 10^{-6}$& 4.0$\times 10^{-6}$& 6.2$\times 10^{-5}$& 9.0$\times 10^{-6}$& 5.7$\times 10^{-5}$\\
&EXP2&8.1$\times 10^{-5}$& 9.0$\times 10^{-6}$& 3.6$\times 10^{-5}$& 2.8$\times 10^{-5}$& 4.0$\times 10^{-6}$& 4.0$\times 10^{-5}$\\
&EXP3&$1.7 \times 10^{-7}$& 1.6$\times 10^{-5}$& $3.3 \times 10^{-7}$& 2.5$\times 10^{-5}$& 4.8$\times 10^{-5}$&1.2$\times 10^{-5}$\\
&EXP4&6.4$\times 10^{-5}$& 2.6$\times 10^{-5}$& 4.9$\times 10^{-5}$& 2.5$\times 10^{-5}$& 1.6$\times 10^{-5}$&  4.8$\times 10^{-5}$\\
&EXP5&1.2$\times 10^{-5}$&4.4$\times 10^{-5}$&2.9$\times 10^{-5}$&4.9$\times 10^{-5}$&2.5$\times 10^{-5}$&1.0$\times 10^{-5}$\\
\bottomrule[3pt]
\end{tabular}}
\caption{\small Performance of FMCA by showing the squared distance between the eigenvalues approximation ($\tilde{\sigma}_1, \tilde{\sigma}_2, \cdots$) and the ground truth (${\sigma}^*_1, {\sigma}^*_2, \cdots$). It can be seen that the approximation is almost highly accurate with a very small squared error.}
\label{table_2_error_plot}
\vspace{-10pt}
\end{figure}

\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=.34\textwidth]{{./used_figures/mine_lr}.pdf}}\hspace{-5pt}%
\subfloat[]{\includegraphics[width=.34\textwidth]{{./used_figures/fmca_figure}.pdf}}\hspace{-5pt}%
\subfloat[]{\includegraphics[width=.34\textwidth]{{./used_figures/fmca_eigen_forth}.pdf}}\\
\vspace{-20pt}\captionof{figure}{\small Comparison between FMCA and MINE. We increase $\rho_\alpha$ from $0.9$ to $0.9999$ and compare the learning curves. As $\rho_\alpha$ approaches $1$, MINE tends to have a very high variance, even with recursive estimators. On the contrary, the learning curves of FMCA are much more stable. The learning curves of the fourth FMCA eigenvalue are also plotted to show that the smaller eigenvalue tends to have a higher variance.}
\label{MINE_COMPARE}
\vspace{-10pt}
\end{figure}

Next, we investigate the CDKF eigenvalues as a novel multivariate statistical dependence measure. Since there is a growing interest in applying mutual information neural estimators (MINE) (\cite{belghazi2018mutual}, \cite{birrell2020f}, \cite{sreekumar2022neural}), the comparison focuses on the implementation of MINE with the Kullback–Leibler divergence (KL) and the second-order R\'enyi's divergence (RY), which are both single-variate measures that produce only one value for the quantification. We show the advantage of CDKF from two perspectives: the functional structure and the stability of FMCA.

\vspace{5pt}
\noindent\textbf{CDKF quantifies functional structures:} To show the advantages of a multivariate statistical dependence measure, we apply FMCA to various datasets such as Gaussian mixtures (GMM), two moons (TM), spirals (SPIRAL) and check dataset (CHECK), and compare with MINE, shown in Table~\ref{table_giant}. Figure~\ref{CIRCLE_PLOT} further visualizes the results in a polar axis. CDKF reveals a more complete hidden functional structure that a single-variate measure cannot quantify. While the mutual information provides the total statistical dependence, CDKF evaluates the statistical dependence across the eigenfunctions of the multivariate space, which when added provides the total estimated dependence. The form $-\frac{1}{20}\sum_{i=1}^\infty \log(1- \widetilde{\sigma}
_i)$ is used in the table and can be directly compared with MINE. However, the hidden structures revealed by the component-wise description can vary for different datasets. For example, CHECK has only one major component and one effective eigenvalue. Applying FMCA to GMM reveals the statistical dependence in each Gaussian component. We also apply techniques from GP to find the ground truth of the eigenvalues for GMM. By comparing the approximation $\Tilde{\sigma}_1, \Tilde{\sigma}_2, \cdots$ with the ground truth ${\sigma}^*_1, {\sigma}^*_2, \cdots$ in Table~\ref{table_2_error_plot}, we further show the high accuracy of FMCA.

\vspace{5pt}

\noindent\textbf{CDKF provides training stability:} A well-known disadvantage of MINE is its instability (\cite{hjelm2018learning}, \cite{hu2021information}) especially near the boundary of strict dependence. In our previous work (\cite{hu2021information}), we described that this issue is caused by the breakdown of absolute continuity due to approximation in the cost function. Moreover, mutual information estimators based on differential entropy for continuous r.v. or r.p. do not have an upper bound, thus will have unacceptably high variance. CDKF instead, has all eigenvalues bounded by $1$ because of normalization techniques. We apply both methods to GAUSS from $\rho_\alpha = 0.9$ to $\rho_\alpha = 0.9999$ shown in Figure~\ref{MINE_COMPARE}. The figure shows that FMCA is much more stable at the boundary of high dependence.

\section{Experiments for modeling real-world images as r.p.}
\label{real_world_applications}
The second set of experiments applies the concept of CDKF and the FMCA to MNIST, CIFAR10, CIFAR100 and the Street View House Numbers (SVHN), which are real-world high-dimensional r.p. with complex and rich functional structures. We will first apply FMCA to classification tasks and investigate the properties and applications of learned eigenvalues and eigenfunctions. Next, factorial coding will be presented as an example of the bi-dimensional nature of the r.p. for generative inference.

\subsection{Learning dynamics of CDKF eigenvalues in classification}

\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=.5\textwidth]{{./used_figures/mnist_curve_v2}.pdf}}%
\subfloat[]{\includegraphics[width=.5\textwidth]{{./used_figures/cifar_curve_v2}.pdf}}\\
\vspace{-30pt}\subfloat[]{\includegraphics[width=.5\textwidth]{{./used_figures/mnist_curve_2_v2}.pdf}}%
\subfloat[]{\includegraphics[width=.5\textwidth]{{./used_figures/cifar_curve_2_v2}.pdf}}
\vspace{-20pt}\captionof{figure}{\small Testing FMCA-CC and FMCA-SC on real-world datasets for classification. Each curve in the plot is a learning curve of a different eigenvalue. We plot the curves with a log-scaled x-axis for better visualization. The figure demonstrates that for real-world datasets, the eigenvalues are also learned sequentially with almost the same convergence speed. It can also be seen that FMCA-CC is more efficient than FMCA-SC.}
\label{real-world-data}
\vspace{-8pt}
\end{figure}

\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\begin{minipage}[h]{1.0\linewidth}
\center
\subfloat[]{\includegraphics[width=.95\textwidth]{{./used_figures/combine_plot_v4_2}.pdf}}%
\end{minipage}\vspace{-25pt}\\
\captionof{figure}{\small Visualizations of FMCA-CC learning dynamics on CIFAR10. The first row shows in the initial stage, the first learned eigenfunction classifies the general machine class against the animal class. The second row shows in the final stage of the training, the CDR classifies the cat class against the dog class. As the CDR $\rho(x, x')$ is an argument between two variables, we fix one sample (marked in red) and calculate all other samples (marked in blue and green). The final row shows the CDR as a matrix over the training.}
\label{CDR_visualization_CLASS}
\vspace{-12pt}
\end{figure}

Applying CDKF and FMCA for classification is implicitly solving the task of quantifying statistical dependence between a given image dataset, which contains realizations of the r.p. and the label r.p.. We compare the class coding (CC) and the similarity coding (SC) introduced in Section~\ref{reference_process} to train the mapping parameters. As a reminder, the CC uses directly the labels as the input to the network $g_\theta$. The SC requires pairs of high-dimensional images from the same class, sampled from every class, as inputs. The learning curves for MNIST and CIFAR10 are shown in Figure~\ref{real-world-data}. In testing, the CDR is constructed  $\Tilde{\rho}(x, x')$ between test set sample $x'$ and samples $x$ from each class in the training sets. The class with the highest CDR is taken as the class prediction for $x'$. This procedure easily produces state-of-the-art classification results, with a test accuracy of $98.1\%$ for MNIST, $95.2\%$ for CIFAR10, $97.3\%$ for SVHN, and $82.1\%$ for CIFAR100. We can observe that CC is much less noisy than SC, as we can expect. In fact, SC resembles clustering because it is based on a similarity measure, instead of a distance measure as most clustering techniques.

The learning dynamics of FMCA shown in Figure~\ref{real-world-data} provides new insight into the classification task under the FMCA criterion. Recall that for simple data sets, FMCA learns the eigenvalues sequentially from the largest to the smallest, with almost the same speed. For real-world datasets, this property means that FMCA during training will classify the classes in the order of their classification difficulty. Taking CIFAR10 as an example, we analyzed the class prediction and the learned eigenfunctions during training in Figure~\ref{CDR_visualization_CLASS}. When the first learned eigenfunction stabilizes, the system separates the images into two subgroups, which are the animals (such as cats, dogs, deer) and the machines (such as vehicles, trucks, planes). Then the classifier further separates the two groups into more subgroups and classes for the remaining eigenfunctions i.e., the system separates the cats and dogs class, birds and deer class, and eventually captures all classes of CIFAR10. This observation is visualized in Figure~\ref{CDR_visualization_CLASS}, which is a unique behavior of FMCA and provides a new interpretation of the classification task, resembling human perception.

\subsection{Implementing factorial codes with FMCA}

\begin{figure}[t]
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\begin{minipage}[h]{1.0\linewidth}
\centering
\subfloat[]{\includegraphics[width=1\textwidth]{{./used_figures/combine_plot_v5}.pdf}}\\
\vspace{-8pt}{\small (a) Visualizations of the CDR $\Tilde{\rho}(x, x')$ with FMCA-MDC.}
\vspace{5pt}
\end{minipage}
\begin{minipage}[h]{1.0\linewidth}
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[]{\includegraphics[width=1\textwidth]{{./used_figures/combine_plot_v6}.pdf}}\\
\vspace{-10pt}{\small (b) Visualizations of the CDR $\Tilde{\rho}(x, x')$ with FMCA-FC.}
\vspace{-5pt}
\end{minipage}%
\captionof{figure}{\small We visualize the CDR $\rho(x, x') $learned by (a) FMCA-MDC and (b) FMCA-FC on MNIST. We vary the number of eigenfunctions for the given sample $x$. The goal of FC is to map image samples to delta functions. When the number of eigenfunctions increases, the delta functions become more accurate. FC performs better than MDC. The codes produced by FC do not contain any similarities between samples and capture the realization component of the r.p., which cannot be achieved by MDC. }
\label{CDR_visualization}
\vspace{-10pt}
\end{figure}

\begin{figure}
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\begin{minipage}[h]{1.0\linewidth}
\subfloat[]{\includegraphics[width=1\textwidth]{{./used_figures/combine_plot_v7}.pdf}}%
\vspace{-20pt}\captionof{figure}{\small Since the CDR has two arguments, we can visualize the CDR between all image pairs across MNIST. The construction with the top $20$ and $200$ eigenfunctions is shown. It can be seen that while both methods try to generate a diagonal matrix, FC is much more accurate and efficient. MDC instead still contains similarities between samples even with $200$ eigenfunctions. }
\label{MDC_MATRIX}
\end{minipage}
\center
\resizebox{0.9\linewidth}{!}{
\begin{tabular}{rl|cc|cc|cc|cc}
\toprule[3pt]
\multicolumn{2}{c|}{NUM OF EIGENV}&\multicolumn{2}{c|}{MNIST}&\multicolumn{2}{c|}{CIFAR10}&\multicolumn{2}{c|}{CIFAR100}&\multicolumn{2}{c}{SVHN}\\
\cmidrule{3-10}
\multicolumn{2}{c|}{}&MDC&FC&MDC&FC&MDC&FC&MDC&FC\\
\midrule[2pt]
$K=20$&$\widehat{\mathbf{f}_\theta}(1:5)$&0.0131&0.0114&0.0193&0.0175&0.0168&0.0165&0.00751&0.00622\\
&$\widehat{\mathbf{f}_\theta}(1:10)$&0.00424&0.00157&0.0113&0.00674&0.01&0.00612&0.00514&0.00204\\
&$\widehat{\mathbf{f}_\theta}(1:15)$&0.00244&0.000898&0.00905&0.00528&0.00824&0.00485&0.00414&0.00129\\
&$\widehat{\mathbf{f}_\theta}(1:20)$&0.0017&0.000678&0.00787&0.00476&0.00719&0.0043&0.00364&0.0012\\
\midrule[1.5pt]
$K=200$&$\widehat{\mathbf{f}_\theta}(1:20)$&0.000787&0.00698&0.0056&0.00426&0.00621&0.00507&0.00136&0.00097\\
&$\widehat{\mathbf{f}_\theta}(1:50)$&0.000459&0.000362&0.00423&0.00315&0.00434&0.00388&0.000862&0.000547\\
&$\widehat{\mathbf{f}_\theta}(1:100)$&0.000337&0.000216&0.00339&0.00263&0.00334&0.00336&0.000652&0.000524\\
&$\widehat{\mathbf{f}_\theta}(1:200)$&0.000255&0.000195&0.00298&0.00206&0.00289&0.00279&0.000546&0.000412\\
\bottomrule[3pt]
\end{tabular}}
\captionof{table}{\small We evaluate the efficiency of MDC and FC numerically by training another neural network to find the inverse mapping from the eigenfunctions to the images, by minimizing the reconstruction error. We train the two networks in FMCA with $K=20$ or $K=200$, then change the number of eigenfunctions for the reconstruction. FC always produces a lower MSE than MDC. As the number of eigenfunctions increases, the reconstruction becomes progressively better. The reconstructed images are shown in Figure~\ref{RECONSTRUCTION_IMAGES}. }
\label{FMCA_UNSUPERVISED_CODING}
\vspace{-5pt}
\captionsetup{type=figure}
\captionsetup[subfigure]{labelformat=empty}
\centering
\hfill\subfloat[]{\includegraphics[width=.25\textwidth]{{./used_figures/FC_TOP_5_RECONSTRUCTION}.pdf}}\hspace{5pt}
\subfloat[]{\includegraphics[width=.26\textwidth]{{./used_figures/FC_TOP_10_RECONSTRUCTION}.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=.27\textwidth]{{./used_figures/FC_TOP_200_RECONSTRUCTION}.pdf}}\hfill
\vspace{-20pt}\captionof{figure}{\small The reconstructed samples with FC for CIFAR10. With only 10 eigenfunctions, the reconstructed images are quite similar to the original.}
\label{RECONSTRUCTION_IMAGES}
\vspace{-10pt}
\end{figure}

The second experiment implements unsupervised learning coding schemes shown in Section~\ref{reference_process} including maximal dependence coding (MDC) and factorial coding (FC). We focus on showing how FMCA is the natural solution to implement factorial codes for feature representation, as well as how it can be employed for generative modeling.

As discussed in Section~\ref{Section_factorial_codes}, an image dataset can be viewed as an r.p. $\mathbf{x} := \{\mathbf{x}(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$ with two coordinates. The r.p. naturally has a bi-dimensional structure consisting of the \textit{pdf} structure (realizations of the r.v. over time) $\{x_n\}_{n=1}^N$, and the spatial structure $\{x_n(t_1, t_2)\}_{t_1, t_2=1}^{T_1}$ for every realization $x_n$. The spatial structure is captured with the eigenfunctions learned by the network $\mathbf{f}_\theta$, as
\begin{equation}
\begin{gathered}
x_n(1, 1), \cdots, x_n(T_1, T_2) \Rightarrow \widehat{f_\theta^{(1)}}(x_n), \cdots, \widehat{f_\theta^{(K)}}(x_n), \; n=1,\cdots, N.
\end{gathered}
\end{equation}
The $\textit{pdf}$ structure is instead captured by the CDR by evaluating every realization $x_n$ against all the others.
\begin{equation}
\begin{gathered}
x_1, x_2, \cdots, x_N \Rightarrow \tilde{\rho}(x_1, \cdot), \tilde{\rho}(x_2, \cdot), \cdots, \tilde{\rho}(x_N, \cdot).
\end{gathered}
\label{CDR_transformation}
\end{equation}
Since the factorial codes require a unique code pattern for every different r.v. realization, our final goal is to make $\tilde{\rho}(x_n, x)$ in Equation~\eqref{CDR_transformation} positive if and only if $x=x_n$.

The comparison results between MDC and FC are shown through visualization in Figure~\ref{CDR_visualization}. As a reminder, MDC employs the same image as the input to both networks $\mathbf{f}_\theta$ and $\mathbf{g}_\omega$, which attains the lower bound of $\Tilde{r}$. However, FC uses presampled non-structured uniform noise as the input of $\mathbf{g}_\omega$ to ensure that each image is assigned to a unique code. Figure~\ref{CDR_visualization} depicts for both methods the CDR approximation using the top 20 and the top 200 eigenfunctions. As depicted, both methods use the full network capacity to assign a unique code for each different realization i.e., every different image sample is mapped to a "delta function". As the number of eigenfunctions increases, the peakiness increases and approximates better the delta function.

Figure~\ref{CDR_visualization} also depicts the codes produced by FC, and they have a much better discriminability than MDC. The MDC codes always include residual similarities between samples. For a given sample (marked in red), the CDR is activated not only for the selected sample, but also for other samples that look similar to it. Even with a large number of eigenfunctions, MDC is significantly less accurate than FC because of this reason. Instead, because FC uses non-structured noise as the input to $\mathbf{g}_\omega$, the codes no longer contain the similarities between samples. To generate perfect factorial codes, FC captures the dissimilarities between samples. And this preserves the independence between realizations of the r.p., which cannot be achieved by MDC. Since the CDR is a quantity between two arguments, we further visualize the CDR as a matrix for all samples in the dataset, shown in Figure~\ref{MDC_MATRIX}. A factorial code is reflected by a diagonal matrix, and Figure~\ref{MDC_MATRIX} further demonstrates that FC produces high-quality factorial codes.

Finally, we compare numerically the two coding schemes for various datasets using the following procedure: By applying FMCA, we learn a mapping from the samples to the eigenfunctions. For reconstruction, we want to find the inverse mapping from the eigenfunctions to the samples, as
\begin{equation}
\begin{gathered}
\widehat{f_\theta^{(1)}}(x_n), \cdots, \widehat{f_\theta^{(K)}}(x_n) \Rightarrow x_n(1, 1), \cdots, x_n(T_1, T_2), \; n=1,\cdots, N.
\end{gathered}
\label{inverse_mapping_reconstruction}
\end{equation}
After training the two networks $\mathbf{f}_\theta$ and $\mathbf{g}_\omega$ with FMCA, we train another neural network using the outputs of $\mathbf{f}_\theta$ by minimizing the reconstruction error to obtain a decoder. The reconstruction network and the networks of FMCA are trained separately. Table~\ref{FMCA_UNSUPERVISED_CODING} shows the comparisons of the reconstruction error. For all image datasets, FC produces a low reconstruction error. By finding the inverse mapping between the codes and the samples, we further demonstrate that FC captures both the realization component and the spatial structure of the r.p.. In Figure~\ref{RECONSTRUCTION_IMAGES} we show the reconstruction results with different numbers of eigenmodes. As expected the accuracy of the reconstruction improves with the number of eigenfunctions, but it is remarkable that with only 10 functions the picture content can be recognized, which shows that the eigenfunctions are very efficient in capturing statistical dependence. This also attests to the fact that the intrinsic image dimensionality is far smaller than the number of pixels we currently employ for image processing. We submit that the eigenvalues of FC can be employed to estimate intrinsic image dimensionality. The next obvious step is to apply the FMCA for generative models of images, but this topic will not be pursued in this paper.

\section{Discussion and Conclusion}

This paper proposes a new cross density kernel function (CDKF) to quantify statistical dependence between two r.p.. The CDKF is associated with a joint $pdf$ of internal variables created by a novel variational Bayesian system named the bidirectional recursion. The CDKF is a major improvement to Wiener filter theory and Gaussian processes because of its nonlinearity and the $pdf$ structure. The CDKF eigenspectrum is the first multivariate measure of statistical dependence between two r.p. and it extends both the maximal correlation and conventional mutual information.

The CDKF eigenproblem can also be solved with realizations of real-world data, using a neural architecture with two networks. The optimization algorithm, FMCA is shown to solve the CDKF variational eigenproblem. Hence, the adaptation of the neural networks extends the recursive updates of the alternating conditional expectation algorithm to a multivariate scenario with an orthogonality and Wiener solution requirement. FMCA is also the first implementation of factorial codes directly from real-world data.

While this paper focuses on the multivariate characteristic of the proposed statistical dependence measure, the framework opens many new doors and can have tremendous applications for information-theoretical learning. A list of potential future applications follows:

\begin{itemize}[leftmargin=*]
\vspace{-5pt}
\item The bidirectional recursion provides a whole new way of describing the relationship between data and internal variables. It exceeds the conventional variational inference methods that consider only a single direction from the latent variable to the data. Based on the bidirectional recursion, it is possible to design a theoretically justified generative model with novel properties that other variational inference methods cannot achieve because of the limitations on the architecture and cost functions employed. Furthermore, in generative adversarial networks, the CDKF eigenspectrum has the capability of providing much more control over the quality of the generated data and avoiding issues like mode dropping or mode collapse.

\item One obvious application of FMCA is to train networks for optimal feature extraction. Currently, there is a lot of interest in low dimensional "universal" features for high dimensional data, using eigendecomposition techniques and information theoretic principles including maximal correlation. The FMCA and the proposed neural architecture of two networks introduced in this paper allow learning these features directly from real data, instead of evaluating them from networks trained with other cost functions.

\item It has been historically shown that Reproducing Kernel Hilbert spaces (RKHS) are natural descriptors of r.p.. The CDKF as a density kernel function motivates broad applications to time series and image analysis. The CDKF is special because it creates an RKHS with the $pdf$ of the data, unlike conventional kernels using predesigned functions e.g., the Gaussian. We should utilize the reproducing and pair-wise properties of the RKHS defined by the CDKF to bring transparency and geometrical intuition for the study of r.p.. Moreover, as this paper uses neural networks as function approximators, we should further use kernels as universal function approximators and perform the optimization in the RKHS with the potential advantage of convexity.

\item The FMCA algorithm can be applied to any problem that requires describing the relationship between multiple signal sources or different modalities. This can be particularly important for neuroscience where cross-domain scenarios are very common.

\item Finally, the perspective of eigenanalysis is extremely important for interpreting the internal representation of neural networks. One major contribution of this paper is to treat the eigenspectrum of the data as a neural network cost function. FMCA is the first practical neural network algorithm to learn eigenfunctions as objectives. Applications of kernel-based learning methods or deep neural networks in physics and chemistry are showing the importance of a more effective alternative to traditional solutions of differential equations for physical processes. However, to be really useful in practice, a particular focus will lie on the challenge of interpreting their solutions. Eigenanalysis of the CDKF appears as a very principled alternative to control the quality of the learned features and also to help interpret what the deep learning models are actually learning.

\end{itemize}

\section{Acknowledgment}
The authors gratefully acknowledge that this work was partially supported by the Office of Naval Research (N00014-18-1-2306) and (N00014-21-1-2345). The author also gratefully acknowledges Dr.Yuheng Bu, Assistant Professor with the Department of Electrical \& Computer Engineering (ECE) at the University of Florida, for constructive discussions and suggestions to improve the paper readability.

\vskip 0.2in