\section{Introduction}
    \label{sec:intro}
    \input{sections/1_introduction}

\section{Multivariate Gaussianization}
    \label{sec:methods}
    \input{sections/3_methods}

\section{Experiments}
    \label{sec:experiments}
    \input{sections/4_experiments}

\section{Conclusions}
    \label{sec:concl}
    \input{sections/6_conclusions}

\section{Acknowledgements}

This research was funded by the European Research Council (ERC) under the ERC-Consolidator Grant 2014 Statistical Learning for Earth Observation Data Analysis. project (grant agreement 647423). J.E.J. thanks the European Space Agency (ESA) for support via the Early Adopter Call of the Earth System Data Lab project; M.D.M. thanks the ESA for the long-term support of this initiative. Additional support was provided by the Project RTI2018-096765-A-100 (MCIU/AEI/FEDER, UE).

\bibliographystyle{IEEEtran}

\subsection{Main Points}
\begin{itemize}
    \item ESDC
    \item Importance of Spatial-Temporal
    \item How do we choose which spatial-temporal decomposition?
    \item Try all models is too expensive...
    \item maybe measure similarity?
    \item how: IT measures and PDF estimation
\end{itemize}

\subsection{Main Contributions}
\begin{itemize}
    \item Spatial-Temporal Data structures
    \item Compare Multivariate Linear and Multivariate PDF Estimation
    \item Information theoretic perspective for measuring
\end{itemize}

\subsection{Similarity}
\commentE{
\begin{itemize}
    \item Briefly introduce similarity - what is it, what do we expect. reference Kornblith similarity paper
    \item introduce IT as the ultimate measure that satisfies some criteria.
    \item \TODO{TODO: shorten a lot...no need for the table}
\end{itemize}}
In this section, we will define what we mean by similarity and what we can expect when using a measure of similarity. Similar to what was presented in \cite{kornblith19a}, we take the perspective of the Earth science applications and argue which properties one would consider important in this setting.
\newline\newline
\textbf{Invariance to Transformations}\newline We prefer similarity measures that are not effected by key transformations. The simplest transform to be invariant to is scaling but a harder transformation to capture is an orthogonal transformation. \newline
1) Scaling are transformations that apply a scalar factor as either a single constant for all features (isotropic scaling, $\alpha \in \mathbb{R}^{+}$) or a constant per feature (non-isotropic scaling, $\alpha \in \mathbb{R}^{d}$). A similarity measure should not be dependent upon isotropic scaling as that would require a lot of preprocessing steps to ensure that the features are all scaled appropriately. However, being independent to non-isotropic scaling is a lot harder to comply with as this is essentially a non-linear transformation and can be useful in some situations. \newline
2) Orthogonal transformations are transformations that preserve scalar products and Euclidean distances between samples. This invariance is important for dealing with ESDC, e.g. there are cases where we need to do dimensionality reduction to reduce the features \cite{Mahecha2010} in order to extract relevant spatial-temporal patterns. A similarity measure will allow one to preserve the key similarity components within two datasets despite applying orthogonal transformation.
\newline\newline
\textbf{Multivariate Methods}\newline We limit ourselves to multivariate methods which are necessary to analyze multivariate data. A common occurrence is to use methods designed for univariate datasets and then do tricks in order to make them work for high-dimensional data. For example, the Pearson correlation coefficient can never be used because as-is, it cannot summarize multivariate information without unraveling the features. However, the congruence-coefficient or the $\rho$-Vector coefficient is a non-linear extension of the standard Pearson correlation coefficient. This class of algorithms usually has a matrix that contain all pairwise points and a matrix-norm is used to summarize this matrix. We will explain more in detail about some individual algorithm characteristics in the next section.
\newline\newline
\textbf{Non-Linear Methods}\newline Linear methods are powerful because they are simple and easy to understand but we argue for the need to use non-linear methods to be able to capture non-linear relationships. The relationship between spatial-temporal data is not always linear and with many features, this relationship can be more difficult to capture using simple non-linear functions. A trade-off we come across often is that non-linear methods tend to come with many parameters that are not always interpretable. For example, we could use non-linear models to measure the $R^2$ values between the input $\mathbf{X}$ and the output $\mathbf{Y}$. The final relationship may also be very sensitive to which parameters are chosen. However, if chosen correctly, one could mitigate the trade-off to reap the best results.

\begin{table}[h!]
\caption{A general breakdown of different methods and their compliance with the requirements for similarity measures.}
\renewcommand{\arraystretch}{1.5}
\resizebox{\linewidth}{!}{%
\begin{tabular}{|l|c|c|c|}
\hline\hline
\multicolumn{1}{|c|}{\textbf{Method}} & \textbf{Non-Linear} & \textbf{Multivariate} & \textbf{Transformations} \\ \hline\hline
\textbf{Linear}                       & No                  & Yes                   & Yes                      \\ \hline
\textbf{Distance-Based}               & Yes                 & Yes                   & Some                  \\ \hline
\textbf{Kernel-Based}                 & Yes                 & Yes                   & Some                  \\ \hline
\textbf{Information Theoretic}        & Yes                 & Yes                   & Yes                      \\ \hline
\end{tabular}%
}
\label{tab:my-table}
\end{table}

\newline\newline
\textbf{Information Theoretic Methods}\newline
Another way to estimate similarity is to estimate the probability distribution function (PDF) directly. We could estimate a PDF for $\mathcal{X}$ and $\mathcal{Y}$ respectively and then compare the samples. There are also summary statistics called Information Theory (IT) metrics which allow us to effectively summarize the information contained within each data distribution. We assume Shannon information: the amount of "surprise", i.e. the expectation of low probability events. Shannon states this as $i({\mathbf x})=-\log p({\mathbf x})$, as any event that is low is inversely proportional to the log of that probability.
To summarize this, we can calculate the Entropy which is a measure of uncertainty or the amount of information within a distribution. This metric is the expected amount of surprise we can expect to see from a given distribution.
\textbf{\begin{equation}\label{eq:ent}
    H({\mathbf x}) = -\mathbf{E}_{\mathbf x} \left[ \log p({\mathbf x}) \right]
\end{equation}}
This measure is analogous to standard deviation.
Alternatively, given two random variables, we can...the amount of information of one variable given another.
\begin{equation}\label{eq:mi}
    I({\mathbf x},{\vect y}) = H({\mathbf x}) + H({\vect y}) - H([{\mathbf x},{\vect y}])
\end{equation}
Like the metrics we mentioned above, if $I=0$, then there is not shared information between $\mathcal{X},\mathcal{Y}$ and they are independent variables. Mutual information has been shown to be a good way for comparing two multidimensional distributions.
This measure is unbounded and to make this comparable, we need a normalized variant. There are many perspectives about how to normalize mutual information. Considering the simplest formulation by \cite{Strehl2002ClusterE}, we can take the analogous method to the Pearson correlation coefficient:
\begin{equation}
    nI(\mathbf{X}, \mathbf{Y}) = \frac{I(\mathbf{X},\mathbf{Y})}{\sqrt{H(\mathbf{X})H(\mathbf{Y})}}
\end{equation}
However, it's not straightforward with the continuous entropy because it can be negative which would affect the squareroot term in the denominator. One strategy is to bias the entropy estimate such that all values are above zero. Although we have added a factor, the trends for Entropy and thus Mutual Information would remain the same. This is the strategy we've chosen in this paper.

\fi

\subsection{PDF Estimation}\label{sec:pdf}

\iffalse
\commentE{
\begin{itemize}
 \item I introduce PDF estimation, what is it?
 \item Why it's important
 \item define the PDF estimation problem somewhat explicitly
\end{itemize}
}
\fi

Most problems in signal and image processing, information theory and machine learning involve the challenging task of multidimensional probability density function (PDF) estimation.
A probability density function or simply a density $p(\cdot)$ takes an input ${\mathbf x}\in{\mathcal X}$ and outputs a density, which follows the properties that 1) $p({\mathbf x}) \geq 0, \forall {\mathbf x}\in\mathbb{R}^{D}$, and 2) it has to sum to one, $\int_{\mathcal X} p({\mathbf x})\text{d}{\mathbf x}=1$.
In practice, we usually do not have access to the PDF $p(\cdot)$ but we do have a set of (multivariate) samples drawn from the generating process ${\mathbf x} = \left\{ {\mathbf x}_1, {\mathbf x}_2, \ldots, {\mathbf x}_N \right\}$ to estimate the PDF from. An accurate PDF estimation is important because it allows us to: 1) calculate the probability of any arbitrary input data point, which accounts for the relative likelihood that the value of the random variable (r.v.) would equal the sample;
2) generate samples ${\mathbf x}'\sim p({\mathbf x})$ from this distribution thus allowing data synthesis, background and support estimation, as well as anomaly detection; and
3) calculate expectations for functions (or transformations) of arbitrary form $f({\mathbf x})$ given $p({\mathbf x})$, i.e. ${\mathbb E}_{\mathbf x}[f({\mathbf x})]$, which allows us to e.g. characterize the system.

Having access to all of these properties gives us the ability to tackle long-standing problems in machine learning and statistics. With accurate PDF estimates, one can model conditional densities of data generated from a prior distribution, develop accurate and efficient compression schemes, and use principled objective functions such as maximum likelihood.
In addition, having access to an accurate density estimator can be useful in many hybrid applications to deal with out-of-sample or out-of-distribution problems too \cite{Nalisnick2019HybridMW}. The problem is therefore to estimate the density $p({\mathbf x})$ given a set of samples from $\mathcal{X}$.

\iffalse
\commentE{
\begin{itemize}
 \item I want to briefly mention parametric models
 \item I want to stress that these methods are good on univariate data
\end{itemize}
}
\fi

The simplest approach to PDF estimation assumes the density has a {\em parametric} functional form defined by a fixed number of tuneable parameters. The Gaussian assumption is the most widely adopted for unimodal distributions, which comes parameterized by a mean $\boldsymbol{\mu}$ and a covariance function $\boldsymbol{\Sigma}$.
If more than one mode is assumed, then a mixture of Gaussians
generally leads to better fits.
\iffalse
\commentE{
\begin{itemize}
 \item I want to mention the core PDF estimators
 \item how histograms have a \textbf{fixed} bandwidth, kernels are \textbf{smooth} and kNN are \textbf{adaptive}...because I like the way it sounds.
 \item I want to stress that they suck with high-dimensional data
 \item I want to stress that they're good on univariate data
\end{itemize}
}
\fi
However, finding a parametric form for the distribution that fits properly to a particular data is really difficult in most cases.

The alternative approach comes from {\em non-parametric models}, which do not assume a specific form for the distribution and are learned from data. The simplest non-parametric method estimates the PDF
by partitioning the data space in non-overlapping bins whereby the density is estimated as the fraction of data points in the bin divided by the volume of the bin. This histogram-based PDF estimation method poorly copes with dimensionality, typically leads to either overfitting or underfitting, and selecting an appropriate number of bins per dimension is a challenge in itself.
Alternative parametric estimates for these methods following likelihood-estimation schemes for the optimal bin width determined by the maximum likelihood have been introduced \cite{Papamakarios2019NeuralDE}. However, they are very {\em rigid}
approaches
and lead to very rough density functions. To achieve {\em smoother} PDF estimates, the kernel density estimation (KDE) method is another popular non-parametric method.
KDE places a non-linear kernel function with a varying bandwidth parameter to control the degree of smoothness on top of each example. Unfortunately, a bias-variance trade-off
will result in over/underfitting the PDF, especially in moderate-to-high dimensional problems.
In the previous approaches, the bandwidth is typically fixed {\em a priori} following heuristics in the literature~\cite{Bishop2007PatternRA}, and rarely take into account the concentration of points, i.e. that smaller bins should be placed in regions with a higher concentration of points, in a form of adaptive bit-allocation scheme. This can be addressed by using $k$-nearest neighbors (kNN), which has an {\em adaptive} bandwidth per location and depends on the number of training points available. However, all of the density estimators above suffer from the curse of dimensionality: as the dimensionality increases, the space becomes sparser and density estimates are unreliable.

\begin{table*}[t!]
    \caption{Summary of all components of the Gaussianization algorithm.}
    \vspace{-0.0cm}
    \label{tab:rbigparts}
    \centering
\begin{tabular}{|c|c|m{3.5cm}|p{1.5cm}|m{1.1cm}|m{1.1cm}|}
\hline \hline
\multicolumn{1}{c}{\textbf{Description}}
& \multicolumn{1}{c}{\textbf{Notation}}
& \multicolumn{1}{c}{\textbf{Transformation}}
& \multicolumn{1}{c}{\textbf{Domain}}
& \multicolumn{1}{c}{\textbf{Before}}
& \multicolumn{1}{c}{\textbf{After}}  \\ \hline \hline
     Marginal Uniformization  & $\text{U}$                                                  & Histogram \cite{Laparra2011RBIG}, \newline Kernel Density Estimation \cite{Meng2020GaussianizationF},\newline Lambert\cite{goerg2012lambert}, Splines\cite{NSF}, Box-Cox \cite{boxcox}
     & $\mathbb{R} \rightarrow \mathbb{R}^{[0,1]}$       &  \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x_u.png}        \\ \hline
     Inverse CDF              & $\text{CDF}^{-1}$
     & Inverse Gaussian CDF, Logit,\newline Inverse Cauchy CDF
     & $\mathbb{R}^{[0,1]} \rightarrow \mathbb{R}$
      &
     \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x_u.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x_g.png}  \\ \hline
     Marginal Gaussianization & $\Psi = \text{CDF}^{-1}\circ \text{U}$                      & Marginal Uniformization + Inverse CDF                                                                                &
     $\mathbb{R} \rightarrow \mathbb{R}$ &
     \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/hist_x_g.png}\\ \hline
     Rotation                 & $\mathbf{R}$
     & Principal Components Analysis \cite{Laparra2011RBIG},\newline Independent Components Analysis \cite{Chen2000Gaussianization}, \newline Random Rotations \cite{Laparra2011RBIG}
     &
     $\mathbb{R}^d \rightarrow \mathbb{R}^d$ &
     \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x_g.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x_l1.png}\\ \hline
     Gaussianization Block    & $G_\ell = \mathbf{R} \left[ \Psi^1 \cdots \Psi^d  \right]$
     & Composition of \newline Rotation + Marginal Gaussianization
     & $\mathbb{R}^d \rightarrow \mathbb{R}^d$
     &
     \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x_l1.png}\\ \hline \hline
     Gaussianization transform
     & $G= \left[G_{1} \circ  \cdots \circ G_{L}  \right]$
     & Composition of \newline Gaussianization Blocks
     &
     $\mathbb{R}^d \rightarrow \mathbb{R}^d$
     &\includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x.png}            &      \includegraphics[width=1.2cm,height=1.2cm]{figures/rbig_components/joint_x_l50.png} \\ \hline\hline
\end{tabular}
\end{table*}
\begin{figure*}
    \centering
    \vspace{3mm} \includegraphics[width = 16cm]{figures/demo/rbig_flow.png}
    \caption{A demonstration of a complete Gaussianization of a noisy sine wave to a marginally and jointly Gaussian distributed. We use {PCA} for the rotation matrix and  {histogram} CDF estimator for the marginal transformation.}
    \label{fig:rbigflow}
\end{figure*}

\subsection{Gaussianization for PDF estimation}

\iffalse
\commentE{
\begin{itemize}
 \item Introduce Gaussianization - original (chen), random rotations (laparra), fully parameterized (meng)
 \item redo the cost function with KL divergence related to the above section
 \item stress on how that representation relates to the IT measures
 \item talk about some of the details (rotations, marginal gaussianization options, computational cost)
\end{itemize}}
\fi

An alternative way to estimate a PDF from observational data is through a data transformation to a {\em convenient} domain, instead of working explicitly in the high-dimensional input domain. The question of what is a convenient domain is a long-standing one, yet ideally this should be a domain with independent components so one can work in each dimension independently to get rid of the curse of dimensionality, one that allows to perform operations and compute quantities therein, and one that is invertible so that one can express these quantities in meaningful units of the input domain.

The Gaussian distribution has desirable properties of showing independent components and being mathematically tractable and is thus a good candidate for density estimation. A class of Gaussianization methods \cite{Meng2020GaussianizationF, Laparra2011RBIG} look for transforms to a multivariate Gaussian domain.
These transforms are related to projection pursuit transformations originally introduced in \cite{Friedman1987ExploratoryPP} and seek to transform a multivariate distribution $p({\mathbf x})$, where ${\mathbf x}\in\mathbb R^d$, into a standardized multivariate Gaussian distribution  \cite{Chen2000Gaussianization,Laparra2011RBIG}:
\begin{eqnarray}
\begin{array}{llll}
      G_{\boldsymbol \theta }: & {\mathbf x}\in \mathbb R^d & \mapsto & \vect{z}\in \mathbb R^d \\
      &   \sim p({\mathbf x}) & &  \sim {\mathcal N}({\bf 0},{\bf I}_d),
\end{array}
\end{eqnarray}
where $\boldsymbol \theta $ are the parameters learned to Gaussianize the data ${\mathbf x}$, $\bf{0}$ is a vector of zeros (for the means) and $\bf{I}_d$ is the identity matrix (for the covariance). By construction, the Gaussianization transform is a parameterized function $G_{\boldsymbol \theta }$ consisting of a sequence of $L$ iterations (or layers), each one of them performing an orthogonal rotation of the data and a marginal Gaussianization transformation to each feature.

The transformation $G_{\boldsymbol \theta }$ in each iteration $\ell$ is defined as:
\begin{equation*}
G_{\boldsymbol \theta }: {\mathbf x}_{\ell+1} = {\bf R}_\ell \boldsymbol{\Psi}_\ell({\mathbf x}_{\ell}),~~~\ell=1,\ldots,L
\end{equation*}
where ${\mathbf x}_{0}$ corresponds to the original data ${\mathbf x}$, $\boldsymbol{\Psi}_\ell$ is the marginal Gaussianization of each dimension of ${\mathbf x}^\ell$ for the iteration $\ell$, and ${\bf R}_\ell$ is a rotation matrix for the marginally Gaussianized variable $\boldsymbol{\Psi}_\ell({\mathbf x}^{\ell})$. After convergence in $L$ iterations, the transformation contains all the needed information to transform data coming from the original density into a multivariate Gaussian. \textcolor{black}{$\boldsymbol \theta $ collectively group all parameters in the method: those from the rotation matrix $\mathbf{R}$ and the marginal transformation $\boldsymbol \Psi $. For example, one could use a principal components analysis (PCA) transformation for the rotation matrix $\mathbf{R}$ and a histogram transformation for the marginal Gaussianization transformation $\Psi$. Then, the eigenvectors obtained from PCA describing $\mathbf{R}$ and the parameterizations of $\boldsymbol \Psi $ would define $\boldsymbol \theta $. See Table \ref{tab:rbigparts} for more details on the decomposition of this formula and Fig. \ref{fig:rbigflow} for a full decomposition of a toy dataset.}

We can use the change of variables formula to calculate the PDF of ${\mathbf x}$ as
\begin{equation}\label{eq:changeofvariables}
     p_{{\mathbf x}}({\mathbf x}) = p_{\vect{z}}\left(G_{\boldsymbol \theta }({\mathbf x}) \right)\left| \nabla_{\mathbf x} G_{\boldsymbol \theta ({\mathbf x})} \right|,
\end{equation}
where $\left| \nabla_{{\mathbf x}} G_{\boldsymbol \theta }({\mathbf x}) \right|$ is the determinant of the Jacobian of $G_{\boldsymbol \theta }$ w.r.t. ${\mathbf x}$. Generally, any unknown PDF of ${\mathbf x}$ can be estimated as long as we have the transformation $G_{\boldsymbol \theta }$ along with its Jacobian.
Intuitively, this transformation essentially destroys the density of $\mathcal{X}$ into unstructured noise (often Gaussian) \cite{Inouye2018DeepDD}. There is no limit to the amount of composite transformations $G_{\mathbf{\Theta}}=G_{\boldsymbol \theta _1} \circ G_{\boldsymbol \theta _2} \circ \cdots \circ G_{\boldsymbol \theta _L}$ that can be used in order to sufficiently converge to the Gaussian distribution. In addition, because $G_{\mathbf{\Theta}}$ is invertible, we can sample points in the original domain ${\mathbf x}'\in\mathcal{X}$ by generating samples in the transformed Gaussian domain and propagating this through the inverse transformation $G_{\mathbf{\Theta}}^{-1}$. Because the transform is a product of linear and marginal operations, both the Jacobian and the inverse transform can be computed easily~\cite{Laparra2011RBIG, sosjaini19a}.

\begin{wrapfigure}{r}{4.4cm}
\vspace{-0.5cm}
\begin{mdframed}[backgroundcolor=gray!20]
The transformation is invertible and allows for density estimation, synthesis and information quantification in high-dimensional Earth data problems
\end{mdframed}
\vspace{-0.5cm}
\end{wrapfigure}

\textcolor{black}{
The original Gaussianization algorithm \cite{Chen2000Gaussianization} worked by applying an orthogonal rotation matrix via independent components analysis (ICA) and then a mixture of Gaussians (MOGs) for the marginal Gaussian transformation. After enough repetitions $L$, it was shown that this converged to a multivariate Gaussian distribution \cite{Chen2000Gaussianization}. In \cite{Laparra2011RBIG} we extended Gaussianization by realizing that the method will converge with any orthogonal rotation matrix $\mathbf{R}$ and we named the algorithm \textit{Rotation-Based Iterative Gaussianization} (RBIG). This allowed for more simpler and faster algorithms such as the Principal Components Analysis (PCA) and even randomly generated orthogonal rotation matrices. In addition, much simpler univariate estimators like the histogram was used to speed up the algorithm significantly.
Meng et. al. \cite{Meng2020GaussianizationF} coined the term {\em Gaussianization Flows} and extended the iterative algorithm to be fully parameterized and trainable by incorporating a mixture of logistics as the marginal Gaussianization layer and a sequence of Householder Flows \cite{Liu2019GMMHF, Tomczak2016HF} as the rotation layer. They also proved this is a universal approximator and showed convincing results that Gaussianization is comparable to some other classes of methods specifically designed for density estimation or sampling \cite{Meng2020GaussianizationF}. All transformations and example variants can be found in table \ref{tab:rbigparts}.}

Irregardless of the method chose, in order to find the parameters $\boldsymbol \theta $ for the transformation $G_{\boldsymbol \theta }$, we minimize the following cost function w.r.t. $\boldsymbol \theta $:
\begin{equation}\label{eq:loss_ng}
 \mathcal{L}(\boldsymbol \theta ) = \text{D}_{\text{KL}} \left[ p_\vect{z}\left(G_{\boldsymbol \theta }({\mathbf x}) \right) || \mathcal{N}({\bf 0},{\bf I}_D) \right],
\end{equation}
which is the Kullback-Leibler (KL) divergence between the estimated Gaussian distribution and the true multivariate Gaussian distribution of mean $\mathbf{0}$ and covariance $\mathbf{I}$; in other words, this is a measure of how much {\em non-Gaussian} our distribution is after transformation. This reveals a direct relationship with information theoretic concepts and measures.
Chen \cite{Chen2000Gaussianization,Cardoso03} showed that \eqref{eq:loss_ng} can be decomposed as
\begin{equation}\label{eq:loss_it}
 \mathcal{L}(\theta) = T({\mathbf x}) + J_m({\mathbf x}),
\end{equation}
where $T({\mathbf x})$ is the Total Correlation (T) (a.k.a. Multi-Information, Multivariate Mutual Information) between all of the marginal distributions, and $J_m({\mathbf x})$ is the KL divergence between the marginal distributions and the standard Gaussian normal distribution. Intuitively, this cost function is trying to minimize the shared information between each of the marginal distributions and ensuring that they follow a standard normal Gaussian distribution.
We want to highlight here that RBIG vastly transforms and simplifies the PDF estimation problem: from estimating the density of the high-dimensional multivariate distribution in $\mathcal{X}$ directly, to doing it indirectly through a transformation to a Gaussian domain. All this by using a series of marginal transformations, which are straightforward and fast.

\begin{figure}[t!]
\small
\begin{center}
\setlength{\tabcolsep}{2pt}
\begin{tabular}{ccc}
\Large$\mathcal X$ & \Large$\hat \vect{z} = \mathbf{G}_\theta({\mathbf x})$ & \Large${\mathbf x} = \mathbf{G}_\theta^{-1} (\hat \vect{z})$\\[0mm]
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/original.png} &
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/approx_base.png} &
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/inverse.png} \\
\Large$\mathcal Z$ & \Large$\hat {\mathbf x} = \mathbf{G}_\theta^{-1} (\vect{z})$ & \Large$\Delta T$, Non-Gaussianity \\[0mm]
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/base_gen.png} &
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/generated.png} &
\hspace{-3mm}
\includegraphics[width=4cm,height=4cm]{figures/demo/rbig/losses.png}
\end{tabular}
\vspace{-0.0cm}
\caption{Density estimation of sinusoid with heteroscedastic noise using RBIG.
Top: The original data distribution $\mathcal{X}$ is mapped to a Gaussian domain $\mathcal{Z}$ with transform $\mathbf{G}_\theta$ parameterized by a set of rotations and marginal Gaussianizations collectively denoted as $\theta$, which has an analytic inverse transformation, ${\mathbf x} = \mathbf D_\theta^{-1} (\hat z)$ to recover the original data.
Bottom: One can sample random data from the Gaussian in domain $\mathcal{Z}$ and use the inverse transformation of $\vect{z}$ to $\hat {\mathbf x}$ for data synthesis. We also demonstrate the losses; the equivalence of the change in total correlation between layers $\Delta T$ and the KL-Divergence between transformed data and a multivariate Gaussian (Non-Gaussianity)}
\label{fig:rbigtoy}
\end{center}
\end{figure}

An illustrative example of how RBIG works on a simple 2D toy dataset is shown in Figure~\ref{fig:rbigtoy}. We transform a non-Gaussian 2D dataset into a 2D marginal and jointly Gaussian distribution along with the inverse transformation (first row). The second row showcases how we can use RBIG to synthesize points in the data domain using the inverse transformation. The bottom right figure shows evolution through iterations of the final total correlation (as a measure of redundancy) and the Non-Gaussianity (as a measure of distance to a Gaussian). Please see

\centerline{\fbox{RBIG site: \url{https://ipl-uv.github.io/rbig/}}}
\vspace{0.25cm}

\noindent for a working implementation of the RBIG algorithm in Python and MATLAB.

\subsection{Information Theory Measures using the RBIG Transform}

\iffalse
\commentE{
\begin{itemize}
 \item Talk about how RBIG calculates the IT measures we mentioned
 \item Information - density estimation
 \item Entropy - for free through the density estimation
 \item Mutual Info - 2 RBIGs and 1 RBIG with the residual entropy of the transformed variables
 \item Probably need a diagram for how RBIG does all of these similar to the ICML poster. It's \textbf{much} easier to understand that way.
 \item Demo Toy Dataset - something similar to the circles, IT measures demo from Valero compared to RV Coefficient.
\end{itemize}}
\fi

RBIG was designed for density estimation, but was inspired by and had connections to information theory~\cite{Cover06}. The series of transformations learned by RBIG lead to a Gaussian domain so features are statistically  independent. This reduction in redundancy is achieved iteratively and can be explicitly computed by summing up all the layer redundancy reductions. This metric is known as the total correlation and computing this metric subsequently allows us to derive {\em some information-theoretic measures} from the data.

\subsubsection{Information}
\label{sec:shannon_info}

Shannon information $I$~\cite{Shannon1948TheMT} is based on the idea that a sample, ${\mathbf x}_i$, is more interesting (carries more information) if it is less probable. The formal definition of information is:
\begin{equation}
    I({\mathbf x}_i) = -\log(p_{\mathbf x}({\mathbf x}_i)).
    \label{eq:information}
\end{equation}
It can be used for instance to highlight regions of more interest in a dataset. Information can be computed for each sample in our dataset by using RBIG and Eq.~\eqref{eq:changeofvariables}.

The expected value of the provided information by a complete dataset, ${\mathbf x}$, is called entropy:
\begin{equation}
    H({\mathbf x}) = \mathbb{E}_{{\mathbf x}} [-\log(p_{\mathbf x}({\mathbf x}))].
    \label{eq:entropy}
\end{equation}
While the entropy could be computed by estimating the information of each sample in the dataset using Eq.~\eqref{eq:information} and averaging, it is more convenient to compute it by using the ability of RBIG to compute Total Correlation as we will see in the following section.

\subsubsection{Total Correlation}

The total correlation, $T$, accounts for the information shared among the dimensions of a multidimensional random variable \cite{Watanabe1960TC,Studeny1998MI}. Details on how to compute $T$ using RBIG can be found in \cite{Laparra2011RBIG}, here we sketch the main idea. Given data  ${\mathbf x} \in \mathbb R^D$, we first learn the Gaussianization transform with $L$ iterations, and compute the cumulative reduction in total correlation in each iteration as:
\begin{equation}\label{eq:T}
 T({\mathbf x}) = \sum_{\ell=1}^{L} \bigg(D~H(\mathcal{N}(0,1)) - \sum_{d=1}^D H({\mathbf x}^{\ell}_d) \bigg).
\end{equation}
The number of layers $L$ will be determined by the reduction in total correlation with each transformation. If there is no change in total correlation after some threshold number of layers, we can assume that $x_d$ are completely independent. It is important to note that all entropy calculations only involved marginal operations which are simple and fast which allows RBIG to be used on large datasets with a high number of dimensions.

\subsubsection{Joint entropy}

While the concept of information is attached to a a particular sample, entropy is used in different fields to characterize how unpredictable a complete process is.
Entropy can be easily computed from the learned RBIG transformation by
\begin{equation}
 H({\mathbf x}) = \sum_{d=1}^D H({\mathbf x}_i) - T({\mathbf x}),
\end{equation}
where $\sum_{d=1}^D H({\mathbf x}_i)$ are marginal entropy estimations and $T({\mathbf x})$ also involves marginal estimations, cf.~\eqref{eq:T}.

\subsubsection{Multivariate Mutual Information}

The multivariate mutual information (MI) accounts for the information shared by two datasets \cite{Cover06}. Estimating MI can be very challenging when working with high dimensional data. Our approach is based on the invariance property of mutual information to reparameterize the space of each variable \cite{Kraskov2004}. Therefore, we essentially Gaussianize the two datasets, $\mathbf{X}$ and $\mathbf{Y}$, with corresponding transforms that remove their total correlations. Then, the total correlation remaining among both Gaussianized datasets is equivalent to the mutual information between the original datasets:
\begin{equation}
MI(\mathbf{X}, \mathbf{Y}) = T([\mathbf{G}_{\theta_{\mathbf x}}(\mathbf{X}), \mathbf{G}_{\theta_{\vect y}}(\mathbf{Y})]),
\end{equation}
which again implies only marginal operations, cf.
eq.\eqref{eq:T}.

Figure~\ref{fig:it_measures3} shows a venn diagram illustrating the different information theory measures used in this paper and Table~\ref{tab:itms} shows a visual demonstration of how they compare to the popular Pearson correlation coefficient for different toy datasets.

\begin{figure}[h!]
\begin{center}
\vspace{0.5cm}
\begin{tabular}{c}
\includegraphics[width=10.0cm]{figures/demo/it_measures/it_measures.png}
\end{tabular}
\vspace{-0.15cm}
\caption{A venn diagram showing the relationships of all information theory measures used in this paper. The solid coloured circles represent marginal variables and the intersection regions with bold lines represent regions for information theory measures like mutual information, $MI$ and total correlation, $T$}
\label{fig:it_measures3}
\end{center}
\end{figure}

\ra{1.2}

\begin{table}[h!]
\caption{\textcolor{black}{Demonstration showing how different information theory measures discussed compare to the popular Pearson correlation coefficient, $\rho$. This table is also a visual demonstration of how to interpret Mutual Information and how it's related to marginal entropy and the joint entropy; $MI(\mathbf{x}, \mathbf{y})=H(\mathbf{x}) + H(\mathbf{y}) - H(\mathbf{x}, \mathbf{y})$}}
\label{tab:itms}
\begin{tabular}{m{1.5cm}ccccccc}
&
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_random_joint.png}
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_parabola_joint.png}
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_noise_joint.png}
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_circle_joint.png}
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_nlinear_joint.png}
& \includegraphics[width=1.8cm, height=1.9cm]{figures/demo/it_measures/itm_linear_joint.png} \\ \hline\hline
\textbf{Correlation}  & $\rho(\mathbf{x}, \mathbf{y})$
& Low  & Medium & Low    & Low  & Medium & High \\
\textbf{Mutual Information} & $MI(\mathbf{x}, \mathbf{y})$
& Low  & Medium & High   & High & High   & High \\
\textbf{Marginal Entropy} & $H(\mathbf{x}), H(\mathbf{y})$
& High & High   & High   & High & High   & High \\
\textbf{Joint Entropy} & $H(\mathbf{x}, \mathbf{y})$
& High & Medium & Medium & Low  & Low    & Low \\ \hline\hline
\end{tabular}
\end{table}
\subsection{Experiment 1: Gaussianization in remote sensing data}

This first set of experiments considers the use of RBIG for standard remote sensing image processing. We will show the performance of RBIG in hyperspectral, very high resolution and radar imagery, and for several applications: joint (multivariate) Gaussianization, data synthesis and information estimation.

\subsubsection{Gaussianization of radar images}

The first part of the experiment focuses on analysis of radar imagery. Data used here was collected in the Urban Expansion Monitoring (UrbEx) \href{http://dup.esrin.esa.int/ionia/projects/summaryp30.asp)}{ESA-ESRIN DUP} project \cite{Chova06}. Results from UrbEx project were used to perform the analysis of the selected test sites and for validation purposes. We consider an ERS-2 SAR pair selected with perpendicular baselines between 20 and 150 m in order to obtain the interferometric coherence from each complex SAR image pair. The corresponding pair ($I_1$, $I_2$) of the SAR backscattering intensities (0-35 days) were stacked for analysis, Figure~\ref{fig:sar}[left]. The relation between the intensity features is strongly nonlinear and non-Gaussian and shows a large dispersion, see Figure~\ref{fig:sar}[a]. The total correlation, $T$ is computed with RBIG for the original domain is $T=0.0929$ bits.
A standard approach in SAR image (pre)processing consists of noise removal and marginal Gaussianization, which can address these problems only partially. This marginal Gaussianization cannot deal with the saturation for high and low signal values, Figure~\ref{fig:sar}[b].
A multivariate Gaussianization leads to a fully Gaussian density, Figure~\ref{fig:sar}[c]. This is confirmed by the estimated total correlation $T=0.0095$ bits as it is less than the marginally Gaussianized data.

\begin{wrapfigure}{r}{4.4cm}
\vspace{-0.5cm}
\begin{mdframed}[backgroundcolor=gray!20]
Multispectral, radar and hyperspectral data exhibit complex nonlinear relations in high dimensional spaces. RBIG maps the data to a convenient Gaussian domain that allows to synthesize new data and quantify information
\end{mdframed}
\vspace{-0.5cm}
\end{wrapfigure}

\subsubsection{Synthesizing hyperspectral images}

\begin{figure*}[t!]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{cc}
\includegraphics[height=3.5cm]{figures/experiments/remote_sensing/sarimage1.png} &
\includegraphics[height=3.5cm]{figures/experiments/remote_sensing/sarimage2.png}\\
\end{tabular}
\qquad
\begin{tabular}{ccc}
(a) & (b) & (c)\\
\includegraphics[height=5cm]{figures/experiments/remote_sensing/sar_joint_original.png} &
\includegraphics[height=5cm]{figures/experiments/remote_sensing/sar_joint_gaussianized_marginally.png}&
\includegraphics[height=5cm]{figures/experiments/remote_sensing/sar_joint_gaussianized_jointly.png}  \\
\end{tabular}
\end{center}
\vspace{-0.5cm}
   \caption{Radar image processing. We illustrate the Gaussinization of 2D radar data comprised of a pair (I1, I2) of ERS-2 SAR backscattering intensities. The joint distribution is non-Gaussian (a) and preprocessing before applying any algorithm is generally convenient. The standard marginal Gaussiaization (b) does not achieve a full spherical (joint) Gaussian, unlike the RBIG transformation (c).
   }
\label{fig:sar}
\end{figure*}

To show the capabilities of the method to deal with high dimensional data here we consider hyperspectral image processing. We took the standard AVIRIS Indian Pines data set~\cite{PURR1947},
where the data has
spectrally redundancy and complex joint distrubtions. The image contains 200 spectral channels, which are considered here as the (very high) input dimensionality. We learned a Gaussianization transform leading to a multivariate Gaussian domain of 200 dimensions spectral bands. Then we sampled from a multivariate Gaussian $n=10^6$ samples of 200 dimensions, and inverted them back to the spectral domain. RBIG can be used this way to generate synthetic spectra easily. Figure~\ref{fig:indianpines} (a)
shows the original and the synthesized spectra. This shows how the proposed method allows us to generate/synthesize seemingly spectral distributions even in such a high-dimensional setting.

\textcolor{black}{Figure~\ref{fig:indianpines} (a) shows the original and the synthesized spectra. This shows how the proposed method allows us to generate/synthesize seemingly spectral distributions even in such a high-dimensional setting. In addition, figure~\ref{fig:indianpines} (b) and (c) show corner plots illustrating the joint distributions between various spectral bands (10, 20, 50, 100, and 150). We see that the marginal and joint distributions for the generated spectra by RBIG in (c) are very similar to the real data in (b) across all pairwise band combinations. This is important to highlight that some of the most widely used methods such as PCA would be able to replicate figure~\ref{fig:indianpines}(a) with a good approximate mean and standard deviation but they would not be able to replicate figure~\ref{fig:indianpines}(d) where {\em all joint distributions} are approximately Gaussian.}

\begin{figure*}[h!]
    \ra{1.2}
\begin{center}
\begin{tabular}{ccc}
(a) Generated spectra & (b) Real & (c) RBIG \\
\includegraphics[height=4.5cm]{figures/experiments/hsi/gen_data_original.png}&
\includegraphics[height=4.5cm]{figures/experiments/hsi/corner_gen_real.png}&
\includegraphics[height=4.5cm]{figures/experiments/hsi/corner_gen_rbig.png}
\end{tabular}
\end{center}
\caption{\textcolor{black}{Gaussianization and synthesis of hyperspectral data using RBIG. In (a) we show the mean and standard deviation spectrum for the 21000 real pixels (mean = black, standard deviation = darker shade) and the 1 million pixels generated synthetically (mean = red, standard deviation = lighter shade) using RBIG.  In (b) and (c) we show the marginal and joint distributions of 10, 20, 50, 100, 150 spectral bands for the real data and data generated with RBIG respectively.}}
\label{fig:indianpines}
\end{figure*}

\subsubsection{Information and redundancy in high resolution images}\label{sec:ucmerced}

Very high resolution images are constantly acquired with the new generation of sensors, both on airborne and spaceborne platforms. A systematic analysis of the images is necessary. Machine learning, and deep learning in particular, has led to an important leap in classification accuracy. However, owing to the wealth of data and the diversity, it becomes necessary to design algorithms that exploit most of the information content of images, both in terms of relevant features and examples.

\iffalse
\begin{figure*}[h!]
\small
\begin{center}
\centerline{\includegraphics[width=18cm]{figures/experiments/remote_sensing/ucmercedpanel.png}}
\setlength{\tabcolsep}{0pt}
\vspace{0.15cm}
\begin{tabular}{cc}
(a) Averaged convergence & (b) Class-specific $T$ \\%& (c) KLD \\
\includegraphics[height=7cm]{figures/experiments/remote_sensing/ucmerced_mi_convergence.png} &  \includegraphics[height=7cm]{figures/experiments/remote_sensing/ucmerced_mi.pdf}
\end{tabular}
\end{center}
\caption{
\textcolor{blue}{Estimation of total correlation, $T$ in aerial imagery. Top: We show four illustrative images for each of the 21 classes in the database.
Bottom: (a) Average total correlation computed iteratively for the different 21 class-specific RBIG models over 50 iterations, with the mean $T$ (solid) and the $T$ standard deviation (shaded) over all models. Convergence is achieved very rapidly for all classes (note the log-scale). (b) ranked $T$ per class computed from the RBIG models.}
}
\label{fig:ucmerced}
\end{figure*}
\fi

\begin{figure*}[h!]
\small
\begin{center}
\scriptsize
\setlength{\tabcolsep}{1pt}
\begin{tabular}{ccccccccccccccccccccc}
forest & chap. & agric.& park & sparse & golf & airplane & beach & river & harbor & tennis & parking & medium & tanks & dense & baseball & inters. & build & freeway & pass & runway\\
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c8_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c6_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c1_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c14_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c19_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c10_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c2_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c4_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c17_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c11_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c21_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c16_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c13_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c20_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c7_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c3_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c12_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c5_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c9_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c15_im1.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c18_im1.png} \\
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c8_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c6_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c1_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c14_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c19_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c10_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c2_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c4_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c17_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c11_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c21_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c16_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c13_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c20_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c7_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c3_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c12_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c5_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c9_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c15_im2.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c18_im2.png} \\
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c8_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c6_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c1_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c14_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c19_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c10_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c2_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c4_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c17_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c11_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c21_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c16_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c13_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c20_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c7_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c3_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c12_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c5_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c9_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c15_im3.png} &
\includegraphics[height=0.80cm]{figures/experiments/remote_sensing/ucmerced_5_examples/ucmerced_c18_im3.png} \\

\end{tabular}

\setlength{\tabcolsep}{0pt}
\vspace{0.5cm}
\begin{tabular}{ccc}
(a) Data partition & (c) Averaged convergence & (d) Class-specific $T$ \\
 \includegraphics[height=2.5cm]{figures/experiments/remote_sensing/FIG_patches_expl.png} & &

 \\
 (b) Measured magnitude & &  \\
 \includegraphics[height=2.5cm]{figures/experiments/remote_sensing/FIG_patches_info.png} & \multirow{-1}{*}[5cm]{\includegraphics[height=5cm]{figures/experiments/remote_sensing/ucmerced_mi_convergence.pdf}} & \multirow{-1}{*}[5cm]{\includegraphics[height=5cm]{figures/experiments/remote_sensing/ucmerced_mi.pdf}}
\end{tabular}
\end{center}
\caption{Estimation of total correlation, $T$ in very high resolution aerial imagery. Top: Three illustrative images for each of the 21 classes in the database, ranked according to their estimated $T$.
Bottom: \textcolor{black}{(a) Each image is decomposed in 3 $\times$ 3 patches with three channels (rgb), making samples of 27 dimensions. (b) We measure how much overlapping there is between the information content (i.e. the total correlation) of these 27 dimensions for each class. We show a Venn diagram to illustrate the measured information following the same criteria as in Fig.\ref{fig:it_measures3}.} (c) Average total correlation computed iteratively for the different 21 class-specific RBIG models over 50 iterations, with the mean $T$ (solid) and the $T$ standard deviation (shaded) over all models. Convergence is achieved very rapidly for all classes (note the log-scale). (d) ranked $T$ per class computed from the RBIG models.
}
\label{fig:ucmerced}
\end{figure*}

Here we validate RBIG to estimate total correlation (multi-information)
in a set of aerial scenes collected in the \href{http://vision.ucmerced.edu/datasets/landuse.html}{UCMerced} data set~\cite{yang2010bag}. The data set contains manually extracted images from the USGS National Map Urban Area Imagery collection from 21 aerial scene categories, with 1-ft/pixel resolution. \textcolor{black}{The data set contains highly overlapping classes and has 100 images per class, see some examples per class in Fig.~\ref{fig:ucmerced}[top].
We extracted color patches of size $3\times 3 \times 3$ from each image, which yielded a total of 6499950 27-dimensional feature vectors per class. Then, we developed a Gaussianization transformation for each class individually and computed the (spatio-spectral) $T$ using RBIG, see Fig.~\ref{fig:ucmerced}[bottom left]. We show in Fig.~\ref{fig:ucmerced}[bottom right] the average and standard deviation of the $T$ evolution through 50 iterations for the 21 classes (note the log-scale) and the total correlation per class. More textured classes like runaways, freeway, buildings and intersections lead to higher $T$, while rather homogeneous/flat classes like chaparral, agricultural or forests reveal low information content.}

\iffalse
\begin{figure*}[t!]
\centerline{
\includegraphics[width=18cm]{figures/experiments/remotesensingpanel.png}}
\caption{
{\sf \bf A} Hyperspectral imagery. Left: The AVIRIS Indian Pines data set 3D representation, where high spatial and spectral redundancy is present. Right: We Gaussianize this high dimensional image (220 spectral bands, dimensions), and then sampled from a multivariate Gaussian $n=10^6$ pixels, inverted them back to the spectral domain and showed the spectra distribution. We show that the proposed method generates/synthesizes seemingly spectral distribution even in such high-dimensional setting.
{\sf \bf B} Radar image processing. Left: Synthetic aperture radar (SAR) data is typically processed before application in remote sensing problems, such as classification or parameter estimation. We show ERS-2 SAR intensity images (I1, I2) and their filtered versions (F1, F2) using specific filters based on median and Wiener filters to highlight urban classes~\cite{Chova06}. Right: The marginal distributions (a) show distinct and non-Gaussian features, which can lead to problems for model development and analysis. Marginal Gaussianization (b) is a standard procedure in radar signal processing, but the joint multivariate Gaussianization noticeably improves results (c).
{\sf \bf C} Vegetation products derived from multispectral imagery are typically used to monitor the terrestrial biosphere. We show fAPAR, FVC and LAI products obtained from SPOT Vegetation (2013-05-13). They are typically highly cross-correlated products, and when analyzed globally or at continental scales lead to highly bimodal, skewed and heavy-tailed distributions. Gaussianization also works well in such situations, and can be a very valuable tool for model and product inter-comparison.
{\sf \bf D} Very high resolution aerial images. Top: We validate RBIG to estimate multi-information (total correlation) in a set of aerial scenes collected in the UCMerced data set. The data set contains manually extracted images from the USGS National Map Urban Area Imagery collection (http://vision.ucmerced.edu/datasets/landuse.html) from 21 aerial scene categories, with 1-ft/pixel resolution. The data set contains highly overlapping classes and has 100 images per class.
Bottom: We developed a Gaussianization transformation for each class individually and computed the MI as a cumulative reduction in MI through iterations. We show the average and standard deviation of the MI evolution through 50 iterations for the 21 classes (bottom left, note the log-scale) and the total correlation per class (bottom right). More textured classes like airplanes or tennis courts lead to higher entropy, while homogeneous or random distributions like chaparral or golf course reveal low information content. The implications of this analysis of the data can help in designing classifier losses, sampling strategies and feature ranking.}
    \label{fig:remotesensingexamples}
\end{figure*}
\fi

\subsection{Experiment 2: Information Quantification of Terrestrial Biosphere Variables in Time}
\label{sec:exps}

\begin{figure}[h!]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{lc}
\multirow{2}{*}[4cm]{\includegraphics[width=0.34\textwidth]{figures/experiments/droughts/map_California_with_star}} &
\includegraphics[width=0.45\textwidth]{figures/experiments/droughts/time_series_SMADI}\\
 & \includegraphics[width=0.46\textwidth]{figures/experiments/droughts/timeSeries_California_USDM}
\end{tabular}
\vspace{-0.15cm}
\caption{Left: distribution of croplands within California, according to the MODIS IGBP land cover classification. Top right: time series of normalized VOD, NDVI, SM, and LST, as well as the SMADI index \cite{smadidata} obtained at the selected pixel (blue star). SMADI extreme drought category is marked with an orange horizontal line. Bottom right: percent area of California in U.S. drought monitor~\cite{usdrought} categories
.}
\label{fig:cali}
\end{center}
\end{figure}

According to climate projections, extreme events are likely to intensify and become more frequent over the coming years \cite{Zscheischler2013}. The effects of extreme events (such as droughts) are prevalent, not only in the biosphere and atmosphere, but in the anthroposphere too. Drought is a major cause of limited agricultural productivity which accounts for a large proportion of crop losses and annual yield variations throughout the world \cite{Boyer1982}. Droughts are also being currently placed as direct contributors to social conflicts, migration, and political unrest (e.g. \cite{Kelley2015}).

There are many studies showing the value of incorporating Earth observation (EO) data for global agricultural systems and applications \cite{Fritz2019,Weiss2020}. Variables such as land surface temperature (LST) and the normalized difference vegetation index (NDVI) derived from optical satellites and, more recently, soil moisture (SM) and vegetation optical depth (VOD) derived from passive microwave sensors are just a few of the many features that can potentially be key for the early detection of drought events \cite{Sanchez2018,Sadri2018,Fernandez_Moran_2017}. The Soil Moisture Agricultural Drought Index (SMADI) was proposed in \cite{Sanchez2016} to integrate SM with LST and NDVI, showing good agreement with other drought indices and with documented events of drought world-wide~\cite{Sanchez2018}.

In this experiment, we quantify the information in and between LST, NDVI, SM and VOD variables for the study area of California (only agricultural fields), see figure~\ref{fig:cali}. LST and NDVI are descriptors of the surface temperature and vegetation chlorophyll content, whereas SM and VOD characterize the water content in soils and vegetation \cite{Sanchez2016,Fernandez_Moran_2017}. We will also use information measures as a means to evaluate whether it would be worthwhile to include VOD as an additional variable in the SMADI ensemble to characterize droughts. Prior to the analyses, variables were resampled into a common 0.05$^\circ$ grid and biweekly temporal resolution. Details on the data sets are provided in Table \ref{tab:expsummary}. Measures are conducted for years 2010-2011 and years 2014-2016 separately, which are representative of non-drought and drought conditions in the study region (see Figure~\ref{fig:cali}).

We focus here in computing multivariate information theory measures in a temporal feature setting, in which previous time steps are included as input features. So for example, 1 input feature includes the current time stamp, 2 input features includes the current time stamp and the time stamp 14 days previously, and so on. This allows us to investigate the temporal scales that maximize the shared information among the remotely-sensed variables. This is particularly relevant for droughts, since there is a time lag between soil/climatic conditions (e.g. represented here by SM, LST) and the plant response (e.g. described by NDVI and VOD), which varies in the literature from two or three weeks up to three months \cite{Petropoulos2017}.

\begin{figure}[h!]
\small
\begin{center}
\setlength{\tabcolsep}{-0pt}
\begin{tabular}{cc}
(a)  \\
\includegraphics[height=5.5cm]{figures/experiments/droughts/H_individual_both.png} & \multirow{0}{*}[5.5cm]{\includegraphics[width=0.25\columnwidth]{figures/experiments/droughts/FIG_VOD_marginal_entropy.png}}\\
(b) \\
\includegraphics[height=5.5cm]{figures/experiments/droughts/H_grouped_both.png} &
\includegraphics[height=5.5cm]{figures/experiments/droughts/FIG_VOD_SMADI_vs.png}\\
\end{tabular}
\vspace{+0.15cm}
\caption{(a) Compares the Entropy for VOD, LST, NDVI and SM individually against the number of temporal dimensions considered. (b) Compares the contribution in entropy of VOD on the joint multidimensional variables integrated in SMADI [LST, NDVI, SM] and SMADI+ [LST, NDVI, SM, VOD] and how it changes as we include more temporal dimensions. The solid lines are mean estimates and the shaded regions are the variance estimates for the non-drought years (2010-2011) and the drought years (2014-2015). \textcolor{black}{Next to each graphic we show a Venn diagram to illustrate the measured information for 3 temporal dimensions as an example, following the same criteria as in Fig.\ref{fig:it_measures3}.}}
\label{fig:droughtH}
\end{center}
\end{figure}

\begin{wrapfigure}{r}{4.4cm}
\vspace{-0.5cm}
\begin{mdframed}[backgroundcolor=gray!20]
How much information is adding a particular variable for drought monitoring? What temporal scales are the most informative? RBIG can answer these questions explicitly in bits
\end{mdframed}
\vspace{-0.5cm}
\end{wrapfigure}
The amount of expected information $H$ for each of the four variables and how it changes as we include more temporal dimensions is analyzed in Figure~\ref{fig:droughtH} (a). Entropy will always increase with more features. So the entropy shown here has been normalized by the total amount of features present which allows us to quantify the amount of {\em entropy per feature}. It can be seen that the amount of entropy for VOD is the highest in all temporal settings, closely followed by LST. All variables decrease in entropy as we add more temporal features. NDVI saturates at around 1.5 bits whereas the other variables have a steady smooth decline. We can also see that LST and VOD show the largest difference between drought and non-drought years, and that the difference is largest as we increase the temporal dimension. This result suggests that LST and VOD observed during longer periods could be more useful in detecting droughts.
Figure~\ref{fig:droughtH} (b) shows that VOD increases the amount of expected information when added to the SMADI variable ensemble in all the temporal settings considered, suggesting that it would be worthwhile to include VOD in agricultural drought studies. Results indicate operational settings of vegetation monitoring could benefit from synergistic approaches that allow including multi-sensor multi-dimensional variables, in particular under stress and disturbances such as agricultural droughts.

The MI of every pair of multidimensional variables was analyzed to investigate the relation and redundancy between them as well as the optimal time scales to combine them. Note standard measures for pairwise comparison such as Pearson's correlation are restricted to one temporal dimension and hence do not allow exploring these scales. The MI scores obtained for LST relations are shown in Figure~\ref{fig:droughtMI}.
Interestingly, it shows that LST-NDVI and LST-VOD show an increase in mutual information up to about 2-4 temporal dimensions and then it saturates. This result suggests that a period of about 1-2 months is needed to capture the soil-plant status with the remotely-sensed variables analyzed in our study region. The curves are relatively similar irregardless of whether it is a drought year or not, and the spread of values for the drought years is considerably reduced for all variables and especially for VOD. This could be related to a reduced variability (limited range of values) under drought episodes, but further studies are needed to confirm this. We also observed that MI is consistently low between SM and all variables with any number of temporal dimensions, and is also low between NDVI and VOD, highlighting the value of combining optical and microwave variables for vegetation/land monitoring.

\begin{figure}[h!]
\small
\begin{center}
\setlength{\tabcolsep}{-0pt}
\begin{tabular}{cc}
\includegraphics[height=5.5cm]{figures/experiments/droughts/LST_I_norm_both.png} &
\multirow{0}{*}[5.5cm]{\includegraphics[width=0.25\columnwidth]{figures/experiments/droughts/FIG_VOD_MI.png}}\\
\end{tabular}
\vspace{-0.3cm}
\caption{Mutual information between pairs of multidimensional variables: LST-VOD, LST-NDVI, and LST-SM.  The solid lines are mean estimates and the shaded regions are the variance estimates for the non-drought years (2010-2011) and the drought years (2014-2015). \textcolor{black}{The Venn diagram illustrates the measured information for 3 temporal dimensions as an example, following the same criteria as in Fig.\ref{fig:it_measures3}.}}
\label{fig:droughtMI}
\end{center}
\end{figure}

\subsection{Experiment 3: Information in Spatial-Temporal Earth data}

\iffalse
\begin{itemize}
\item More data description + ESDL
\item Mega-maps of information (triangle + anomalies)
\item isolated areas, (zoomed in regions) + bars per biome
\item Entropy to summarize, H
\end{itemize}
\fi

\subsubsection{Data}

For our experiments we used observational and model simulated variables from the Earth Science Data Lab (ESDL)~\cite{Mahecha19esdc}\footnote{\url{https://www.earthsystemdatalab.net/}}, which is a platform that provides an opportunity for data-centric processing methodologies.
The analysis-ready data-cube contains and harmonizes more than 40 variables relevant to monitor key processes of the terrestrial land-surface and atmosphere.
Data exhibit clear spatial-temporal relations, which need to be taken into account to properly convey and quantify information. Figure~\ref{fig:ESDC_data} illustrates how we represent this spatial-temporal relations as inputs given a single variable. Here we focus on three key
land-surface variables:
precipitation,
sensible heat and evaporation, which are outlined below:
\begin{itemize}
\item {\em Precipitation, $Precip$.} This is a fundamental variable in land-atmosphere processes. The collected data comprises the period 1980--2015, and comes from the Global Precipitation Climatology Project (GPCP)~\cite{ADLERetal2003,HUFFMANetal2009}.

\item {\em Sensible heat, $SH$.} These data comprises 2001--2012, and was generated by training an ensemble of machine learning algorithms with eddy covariance data from FLUXNET and satellite observations in a cross-validation approach, regressions from these observations to different kinds of carbon and energy fluxes were established and used to generate data sets with a spatial resolution of 5 arc-minutes and a temporal resolution of 8 days. The H resembles the sensible heat flux from the surface and is expressed in [W~m$^{-2}$]~\cite{TRAMONTANAetal2016}.

\item {\em Evaporation, $E$.} These data covers 2001--2011, and builds on the Global Land Evaporation Amsterdam Model (GLEAM), which consists on a set of algorithms that separately estimate the different components of land evaporation using input forcing data sets from reanalyses, optical and microwave satellites and other merged sources. The model itself consists of four modules: potential evaporation (Priestley and Taylor equation), interception (Gash analytical model), soil (multilayer soil model plus data assimilation) and stress (semi-empirical). The data are sampled on a grid of 0.25$^\circ$ and have a daily temporal coverage~\cite{MARTENSetal2017,MIRALLESetal2011}.
\end{itemize}
The data is organized in a 4-dimensional data cube ${\bf x}(u,v,t,k)$ involving (latitude, longitude) spatial coordinates $(u,v)$, time sampling $t$, and the variable $k$. The available data is provided at two spatial resolutions (0.083$^o$ and 0.25$^o$) and at a temporal resolution of 8 days, spanning the years 2001-2011. In our experiments, we focus on the lower resolution products and on the period 2008-2010.

\subsubsection{Spatial-Temporal Analysis}

The considered variables (precipitation, sensible heat, evaporation) are fully coupled.
Moisture and precipitation interactions are vastly modulated by both land-atmosphere exchanges and large-scale atmospheric circulation.
Nevertheless, before understanding variable relations, it is important to {\em identify when and where} individual variables are expressive. This may help in assessing the coupling
mechanisms between variables and improve Earth system models.

\begin{wrapfigure}{r}{4.4cm}
\vspace{-0.5cm}
\begin{mdframed}[backgroundcolor=gray!20]
What are the optimal spatial and temporal scales to exploit each variable in the coupled land-atmosphere system? RBIG sheds light on uncertainty quantification and information in multidimensional spatio-temporal Earth data cubes
\end{mdframed}
\vspace{-0.5cm}
\end{wrapfigure}
The question we want to address in this experiment is what are the optimal (in information terms) spatial and temporal scales to exploit each variable's information. Using RBIG, we show here that the ratio of spatial-temporal neighbouring pixels giving the most amount of information can be explicitly calculated. We used RBIG to calculate the entropy $H$ for the aforementioned variables under different spatial-temporal configurations (fully temporal, spatio-temporal and fully spatial) as well as the corresponding information $I({\mathbf x})$ for each time pixel and variable.

Figure~\ref{fig:entropypanel} shows the entropy for the different variables and configurations following the same procedure as in \cite{Laparra2015} (and used in experiment \S\ref{sec:ucmerced} too in the spatial domain only).
Essentially we formed cubes with the same dimensionality but different spatio-temporal configuration and computed the entropy values for each of them.
We chose several configurations ranging from a ratio of purely spatial (ratio=0) up to purely temporal (ratio=1). We also looked at different configurations for the amount of spatial-temporal dimensions used, e.g. a maximum of 4 dimensions up to a maximum of 49 dimensions (temporally, this is approximately 1 year).
Notice how each variable has a different spatial-temporal relationship with entropy, but in general temporal configurations (ratio=1) convey more information than purely spatial (ratio=0) for all the considered variables.
Trends are clear in particular for precipitation, where incorporating temporal information for any amount of dimensions has higher expected information.
For sensible heat and evaporation the entropy paths are similar, and reveal a fast increase in entropy for particular spatio-temporal configurations (ratio$\sim$0.8).
These results suggest different {\em optimal} (in information terms) time and space scales for different variables, which may have implications in further analyses and applications.

Using the same data configurations we have computed the information content of each sample following the procedure described in section \ref{sec:shannon_info}. This help us to visualize the regions with more and less information. We show in Fig. \ref{fig:spatiotemporalpanel} the results of a spatio-temporal analysis of the information content of all three variables.
In regions where we expect pronounced seasonal patterns, the information (complexity) is apparently high in fully temporal configurations as the seasonal cycle controls ecosystem dynamics.
Actually, seasonal (temporal) modes are of lower informative content in the spatial domain, as they are mainly driven by solar forcing.
The information values tend to be higher in tropical regions, whereas arid regions show low-complexity (low-information) patterns. Let us now look in deeper detail at the different spatio-temporal configurations and their information patterns.

Global patterns in rainfall are traditionally related to a strong seasonality, dominated by the position of the Inter-Tropical Convergence Zone (ITCZ) in the tropics, and the El Ni\~no-La Ni\~na cycles, which occur irregularly at intervals of 2-7 years. Spatial information generally dominates with high probability in the Amazonia and the tropics and with low information in desertic areas (e.g. California, Arabian peninsula and central Australia). As we quantify information in spatio-temporal configurations, more clear patterns of low information (e.g. Australia) and high probability (e.g. east-west US gradient) emerge \cite{tuttle2016empirical}. Studying precipitation in the fully temporal configuration translates into a clear ruling of the winter season in Amazonia, Indonesia, as well as northern Europe. Yet, a comparison of temporal vs. spatial information in Fig.~\ref{fig:spatiotemporalpanel}[bottom row] reveals that spatial information dominates in desertic areas (e.g. Australia, Iberian peninsula, Sahara, Mexico) which are reasonably independent of time, and temporal information dominates in Sahel (Savanna), northern latitudes and SW China, which are generally characterized by high rain factors, seasons and moisture.

Transfer of sensible heat $SH$ into the air is dependent on the temperature gradient between the surface and the air above. Patterns of the information of sensible heat $SH$ stand out clearly. While the (fully) spatial information dominates in the Northern hemisphere, the (fully) temporal information patterns appear in the tropics where rainfall is present over larger regions and seasons. The global spatial distribution of $SH$ information shows the largest values in subtropical dry regions where available energy is preferentially partitioned to sensible heat rather than latent heat~\cite{jung2011global}, and seem to be anti-correlated with the amplitude of the mean seasonal cycle.
These results reveal a maximum information of $SH$ in the tropical and subtropical deserts, where the high surface temperature conducts much heat into the air above, and the lowest near the poles where the surface temperatures are much lower.
Information is mainly concentrated in the tropics too, and show similar patterns to precipitation with the exception of clear spatial information in the Indian peninsula.
Evaporation maps of information reveal that the spatial information dominates in deserts and dry regions where evaporation is limited, while temporal information (more interannual variability) resides in Northern latitudes. This is mainly due to the low temperatures and radiation which relates to little evaporation all year round. The temperate areas show increased evaporation information in both purely spatial and temporal configurations, coinciding with increasing temperatures over ground moistened by winter rains. Cooler winter temperatures in Southern hemisphere reduce evaporation, which is also captured in the spatial-vs-temporal divergent maps. Note that in very dry regions information is higher (lower evaporative fraction), conversely for very humid regions, in agreement with~\cite{jung2011global}.

\iffalse
\begin{figure*}[t!]
    \centering
    \small
    \setlength{\tabcolsep}{1pt}
    \begin{tabular}{c|cccc}
    & precipitation & root moisture & sensible heat & evaporation \\
    \rotatebox{90}{\hspace{5mm} Entropy} &
    \includegraphics[width=3.85cm]{figures/experiments/spatial_temporal/temp/entropy/esdc_precipitation.png} &
    \includegraphics[width=3.85cm]{figures/experiments/spatial_temporal/temp/entropy/esdc_root_moisture.png} &
    \includegraphics[width=3.85cm]{figures/experiments/spatial_temporal/temp/entropy/esdc_sensible_heat.png} &
    \includegraphics[width=3.85cm]{figures/experiments/spatial_temporal/temp/entropy/esdc_evaporation.png} \\
    \rotatebox{90}{\hspace{10mm} $7\times 7 \times 1$} &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var1_config1.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var3_config1.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var4_config1.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var5_config1.png}  \\
    \rotatebox{90}{\hspace{5mm} $3\times 3 \times 6$} &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var1_config2.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var3_config2.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var4_config2.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var5_config2.png}  \\
    \rotatebox{90}{\hspace{5mm} $1\times 1 \times 46$} &        \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var1_config3.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var3_config3.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var4_config3.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/var5_config3.png}  \\
    \rotatebox{90}{\hspace{5mm}  Spat-temp. $p(x)$} &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/compo_var1.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/compo_var3.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/compo_var4.png}  &
    \includegraphics[width=4.33cm]{figures/experiments/spatial_temporal_2/compo_var5.png}  \\

    \end{tabular}
    \caption{ (A) The entropy information theory measure for different ratios of spatial-temporal dimensions (ratio=0 meaning fully spatial while ratio=1 means fully temporal information). (B) Information content maps for several variables (precipitation, root zone soil moisture, sensible heat, and evaporation) using different spatial-temporal cube configurations (a) $7\times 7$ spatial width, $1$ temporal length, (b) $3\times 3$ spatial width, $6$ temporal length and (c) $1\times 1$ spatial width, $46$ temporal length. (C) map of the trade-off between the spatial ($7\times 7\times 1$) and fully temporal ($1\times 1\times 46$) information.
    \label{fig:spatiotemporalpanel}
\end{figure*}
\fi

\subsection{Entropy}

\begin{figure}[h]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{m{45mm}m{45mm}}
H (Root Mosture) & Normalized H (Root Moisture) \\
\includegraphics[width=\linewidth]{figures/prelim/rm_entropy.png} &
\includegraphics[width=\linewidth]{figures/prelim/rm_n_entropy.png} \\
H (LST) & Normalized H (LST) \\
\includegraphics[width=\linewidth]{figures/prelim/lst_entropy.png} &
\includegraphics[width=\linewidth]{figures/prelim/lst_n_entropy.png} \\
\end{tabular}
\caption{Root Moisture and Soil Moisture for different spatial-temporal values. We show the entropy and normalized entropy} \label{fig:spatemp_h}
\vspace{0.15cm}
\end{center}
\end{figure}

\subsection{Correlation Measures}

\begin{figure}[h]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{m{45mm}m{45mm}}
Mutual Info & (b) Normalized Mutual Info  \\
\includegraphics[width=\linewidth]{figures/prelim/rm_mi.png} &
\includegraphics[width=\linewidth]{figures/prelim/rm_nmi.png} \\
RV Coefficient & (b) CKA  \\
\includegraphics[width=\linewidth]{figures/prelim/rm_rv.png} &
\includegraphics[width=\linewidth]{figures/prelim/rm_cka.png} \\
\end{tabular}
\caption{Spatial-Temporal results for \textbf{Root Moisture}.} \label{fig:test_stat}
\vspace{0.15cm}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{m{45mm}m{45mm}}
Mutual Info & (b) Normalized Mutual Info  \\
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_mi.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_nmi.png} \\
RV Coefficient & (b) CKA  \\
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_rv.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_cka.png} \\
\end{tabular}
\caption{Spatial-Temporal results for \textbf{Gross Primary Productivity}.} \label{fig:test_stat}
\vspace{0.15cm}
\end{center}
\end{figure}

\subsection{Distnace Measures}

\begin{figure}[h]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{m{45mm}m{45mm}}
Random Forest, RMSE) & Linear Regression, RMSE) \\
\includegraphics[width=\linewidth]{figures/prelim/models/sm_rf_rmse.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/sm_rlr_rmse.png} \\
Variation Information & RV RMSE \\
\includegraphics[width=\linewidth]{figures/prelim/models/sm_vi.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/sm_rv.png} \\
\\
\end{tabular}
\caption{Soil Moisture for different spatial-temporal relations. } \label{fig:dist_sm}
\vspace{0.15cm}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\setlength{\tabcolsep}{0pt}
\begin{tabular}{m{45mm}m{45mm}}
Random Forest, RMSE) & Linear Regression, RMSE) \\
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_rf_rmse.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_rlr_rmse.png} \\
Variation Information & RV RMSE \\
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_vi.png} &
\includegraphics[width=\linewidth]{figures/prelim/models/gpp_rv.png} \\
\\
\end{tabular}
\caption{Gross Primary Productivity for different spatial-temporal relations. } \label{fig:dist_gpp}
\vspace{0.15cm}
\end{center}
\end{figure}