\section{Data Pre-processing}\label{sec:preproc}
The main issue of applying optimization techniques on geographical data refers to the high data dimensionality, especially whether vectorial and matrix computations are involved.\\
In the Machine Learning domain, this is often addressed as \textit{curse of dimensionality}, i.e. the impossibility of efficiently learning the space where data samples are distributed, which becomes exponentially sparse as the dimensions of the feature space increase. For this reason, here comes the necessity to introduce some pre-processing techniques, able to create a suitable low-dimensional encoding, both minimizing the loss of information, and maximizing the performance of the optimization process.
\subsection{Empirical Orthogonal Function (EOF) Analysis}
Empirical Orthogonal Functions (EOFs) \parencite{eof} represents a standard reduction technique in geo-science-related experiments. Firstly adopted to reduce the feature space to few relevant variables (sharing a common concept with PCA, Section \ref{subsubsec:pca}), it is mainly used to extract individual modes of variability (i.e. pre-dominant patterns, often sharing common structures with the main regimes). There exist several variants of the algorithm, according if a spatial-time cross-correlation is kept into account; in this section, we refer to the standard algorithm, where the components are considered to be totally de-correlated (orthogonal). Although we know that the real physical components share some dependencies, we believe this is a suitable trade-off leading to a minor computational complexity.\\
Given the time-indexed geo-dataset, we can indicate it as:
$$ \mathbf{X}_{t,y,x},\quad t=1,\dots, N$$
i.e., indexed by time $t$, latitude $y$, longitude $x$, respectively.\\
The EOF procedure requires a 2D-flattened version of this dataset, where the spatial indexes are converted from a grid to a vectorial format, leading to:
$$ \mathbf{X}_{t,p},\quad t=1,\dots, N,\ p=(y_i, x_i), \ i = 1,\dots, P$$
Since data is provided in a regular grid format, there could be the necessity to introduce a weighting scheme prior to analyse data and fed it to EOF analysis. Indeed, the density given by the un-weighted grid, would be major poleward, than along the latitude. This attempt requires then the following mathematical passage:
$$ \mathbf{X}^w = \mathbf{X} D_\theta, \quad D_\theta = \mathrm{diag}[\cos\theta_1, \dots, \cos\theta_P]$$
Finally, the Empirical Orthogonal Functions are computed solving the \textit{Singular Value Decomposition} (SVD):
$$ \mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T$$
factorizing the dataset in:
\begin{itemize}
    \item $\mathbf{U} \in \mathbb{R}^{n \times n}$, whose columns are called \textit{left singular vectors}, representing the eignevectors of the matrix $\mathbf{X}\mathbf{X}^T$, as well as the EOFs along the spatial dimension
    \item $\mathbf{\Sigma} \in \mathbb{R}^{n \times p}$: the matrix of the singular values
    \item $\mathbf{V} \in \mathbb{R}^{p \times p}$, whose columns are called \textit{right singular vectors}, representing the eignevectors of the matrix $\mathbf{X}^T\mathbf{X}$, as well as the Principal Components (PCs) along the time dimension
\end{itemize}
By construction, EOFs are stationary (i.e. they do not evolve over time), only changing in the sign and value of their amplitude, thus resembling the state of the atmosphere.
\begin{figure}[!htbp]
    \centering
    \includegraphics[width = .65\textwidth]{imgs/eofs.png}
    \caption{Visualization of the four main EOFs (Empirical Orthogonal Functions) in the North-Atlantic/European zone. It is possible to identify: a NAO+ resembling component (top left); an anti-Scandinavian Block SB- (top right); an anti-Atlantic Ridge AR- component (bottom left)}
    \label{fig:eofs}
\end{figure}
\begin{figure}[!htbp]
    \centering
    \begin{minipage}{.68\hsize}
       \includegraphics[width = \textwidth]{imgs/PC1.png}
    \end{minipage}
    \begin{minipage}{.31\hsize}
       \includegraphics[width = \textwidth]{imgs/density.png}
    \end{minipage}
    %\includegraphics[width = .85\textwidth]{imgs/PC1.png}
    \caption{Bar Visualization (left) of the main Principal Component (PC-1) of the Time Series, in the North-Atlantic/European zone, grouped by the January-February-December months of the same solar year. PC is here standardized to have zero mean and unit variance. When the values are out of the range $\pm 1$, they are generally associated to more extreme meteorological events. Associated to it, the NAO Index \parencite{naoindex} evolution during winter months: it is evident the high correlation with PC-1. On the right, a density visualization of the first two PCs}
    \label{fig:pc1}
\end{figure}

\subsection{Principal Component Analysis (PCA)}\label{subsubsec:pca}
PCA represents a standard pre-processing technique for dimensionality reduction in Machine Learning problems. As for EOF analysis, the intent is to find a low-dimensional (latent) representation $\mathbf{Z}$ of data, through a projection matrix $\mathbf{W}$, such that:
$$ \mathbf{Z} = \mathbf{XW}$$
by maximizing the variance (and thus, the explainability) of it. It is an optimization problem, described by the objective:
$$ \underset{w:\ \|w \| <1}{\mathrm{argmax} \ } \| \mathbf{X}w\|_2^2$$
where we introduce a constraint on the projection weights $w$, in order to bound the problem and regularize the magnitude of them.\\
As a result, it translates into a closed-form analytical solution, coinciding with the eigen-decomposition problem. The projection matrix $\mathbf{W}$ is thus the matrix of eigen-vectors of $\mathbf{X}$. By picking the top-$n$ eigenvectors, associated to the largest-magnituded eigenvalues, we can control the quantity of explained variance associated to the latent variables $z$.
\subsection{Variational Autoencoder (VAE)}
As clearly proved, the PCA approach performs a linear transformation of data into an encoded space. This feature represents both the strength and weakness of the method, able to computationally scale with large datasets, but often inferring a trivial mapping.\\
For this reason, we also provide a more advanced dimensionality reduction technique, leveraging the approximation power of neural networks.\\
Autoencoders (AEs) are a Deep Learning architecture, composed of two stages (Figure \ref{fig:ae}): an \textit{encoder}, mapping original data into a low-dimensional representation, called \textit{bottleneck} (in this case, we deal with \textit{undercomplete} AEs); a \textit{decoder}, learning to reconstruct the input samples from the latent space.
\begin{figure}[!htbp]
    \centering
    \includegraphics[width = .7\textwidth]{imgs/autoencoder.png}
    \caption{Visualization of an Autoencoder architecture\parencite{ae_image}}
    \label{fig:ae}
\end{figure}
The simplest architecture we can think of, is a 3-layers network. In the absence of non-linear activation functions, this architecture should converge to the PCA, being the bottleneck obtained through a matrix-to-matrix multiplication with the input.\\
In our study, we decide to tackle the problem in a more sophisticated way, directly inferring a mapping on the anomaly maps \parencite{dim_red_AEs}, thus employing techniques from the computer vision domain \parencite{DBLP:journals/corr/abs-2005-04613}.\\
We switch to a probabilistic approach, employing a Variational Autoencoder (VAE) architecture. By definition, a VAE is an Autoencoder where the latent space summarizes a multi-variate distribution, learning its mean and covariance parameters (Figure \ref{fig:vae}).
\begin{figure}[!htbp]
    \centering
    \includegraphics[width = .85\textwidth]{imgs/vae.jpeg}
    \caption{Visualization of a Variational Autoencoder (VAE) architecture\parencite{vae_image}}
    \label{fig:vae}
\end{figure}
The ``variational'' addressing is referred to the statistical approximation performed by the two networks, optimizing the so-called \textit{variational objective}. Indeed, under a mathematical point of view, a VAE tries to find the distribution $p(z|x)$ which best explains data $x$ starting from its low-dimensional representation $z$. Recalling the Bayes' rule we can fragment this conditional \textit{posterior} by means of a \textit{likelihood} $p(x|z)$ on the reconstructed data, a \textit{prior} $p(z)$ on the latent space, and a \textit{marginal likelihood} (often addressed as model evidence) $p(x)$. Therefore, according to the Bayes' theorem it follows that:
$$ p(z|x) = \frac{p(x|z) p(z)}{p(x)}$$
Under a more practical point of view and with a slight change of notation, the decoder, parametrized by $\theta$ parameters, approximates the likelihood distribution $p_\phi(x|z)$, while the encoder ($\phi$-parametrized) is in charge of approximating the posterior $p_\phi(z|x)$. Finally, the prior on the latent space is typically chosen to be as more general as possible, with a standard multivariate Gaussian representing the usual choice.\\
The framework presents an innermost issue, though. The marginal likelihood term $p(x)$ is indeed tough to be estimated analytically, leading to an intractable analytical solution. Hence, we resort to \textit{variational inference}, introducing an approximating distribution $q(z)$, which we optimize to be as closer as possible to the true posterior distribution.\\
The objective statement then aims at reducing the gap between these two distributions, by minimizing their \textit{Kullback-Leibler (KL-)Divergence}:
$$ \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ }\mathbb{KL}(q_\phi(z) \| p_\phi(z|x))$$
By definition, we can write the KL-term as:
$$ \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ }\mathbb{KL}(q_\phi(z) \| p_\phi(z|x)) = $$ $$ = \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ }\int q_\phi(z) \log \frac{q_\phi(z)}{p_\phi(z|x)} = \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ }\mathbb{E}_{q_\phi(z)}\Bigg[\log \frac{q_\phi(z)}{p_\phi(z|x)}\Bigg]$$
Then, exploiting the logarithm properties:
$$ \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ } \mathbb{E}_{q_\phi(z)}[\log q_\phi(z)] - \underbrace{\mathbb{E}_{q_\phi(z)} [\log p_\phi(z|x)]}_{\text{probabilistic term}}$$
Focusing on the probabilistic term, we can apply the Bayes' rule to simplify a bit the notation:
$$ \mathbb{E}_{q_\phi(z)} [\log p_\phi(z|x)] = \int q_\phi(z) \log \Bigg[ \frac{p_\theta(x|z) p(z)}{p(x)}\Bigg]$$
$$ = \mathbb{E}_{q_\phi(z)} \log p_\theta(x|z) + \mathbb{E}_{q_\phi(z)} p(z) - \mathbb{E}_{q_\phi(z)} p(x)$$
Recalling the initial objective, where we intentionally omit the $\mathrm{argmin}$ notation for sake of brevity:
$$ \mathbb{KL}(q_\phi(z) \| p_\phi(z|x)) = \mathbb{E}_{q_\phi(z)}[\log q_\phi(z)] - (\mathbb{E}_{q_\phi(z)} \log p_\theta(x|z) + \mathbb{E}_{q_\phi(z)} p(z) - \mathbb{E}_{q_\phi(z)} p(x))$$
With a bit of re-arrangement, we can define the variational objective, also known as \textit{Evidence Lower Bound (ELBO)}:
$$ p(x) - \mathbb{KL}(q_\phi(z) \| p_\phi(z|x)) = 
\mathbb{E}_{q_\phi(z)} \log p_\theta(x|z) - \underbrace{(\mathbb{E}_{q_\phi(z)}[\log q_\phi(z)] - \mathbb{E}_{q_\phi(z)} p(z))}_{= \mathbb{KL}(q_\phi(z) \| p(z))}$$
Given the change of sign, this new objective is maximized, and it is equivalent to the maximization of the marginal likelihood metric, representing a lower bound on it. Indeed, by minimizing the KL-divergence between the approximating and true posteriors, we try to make the lower bound tighter. At the same time, the log-likelihood on the reconstructed data is maximized. The last term represents a \textit{regularization} term, which explains the intended evolution of the optimization process: at the beginning the approximating distribution $q_\phi(z)$ is thought to be very close to the uninformative prior of the latent space $p(z)$. The scope is then to gradually shift the approximating distribution towards the more complex posterior.\\
Under a more practical point of view, it is often hard to find the correct balance to reach the desired performances. Indeed, the risk is that over-parameterized model could lead to an excessive weight of the last term, often resulting in trivial or weakly-distinguished encodings. For this reason, over the years several variants of VAEs have been proposed, with the scope of introducing more dynamic objectives, with twofold purposes:
\begin{itemize}
    \item leading to better reconstructions
    \item obtaining the so-called \textit{disentangled} representations \parencite{DBLP:journals/corr/abs-1811-12359, DBLP:journals/corr/abs-1905-13662, DBLP:journals/corr/abs-1905-01258,pmlr-v119-locatello20a,pmlr-v139-trauble21a,DBLP:journals/corr/abs-2010-14407, mita2021identifiable}
\end{itemize}
While the former is often more related to image-reconstruction problems, where standard VAEs tend to produce blurry reconstructions, the latter has been at the center of many studies. We define disentagled representations as a suitably small set of latent features, uncorrelated between each other, and thought of being sufficiently informative for the generation process. As a standard example, characteristics like the geometry, shape, size, colour, orientation of an object in an image could be tought as disentangled features. Being able to retrieve them means finding a very compressed feature mapping, enhancing the overall data and model interpretabilities.\\
A drawback of this approach is that it often requires supervised information for at least few data samples, turning the problem into a supervised/weakly-supervised learning problem. In our study, this represents a big obstacle since we do not have at disposal any labeling. Also, there could not be the possibility to artificially re-create a pseudo-labeling, since each anomaly map could fall into multiple classes with different probabilities, making an hard-classification impossible.\\
However, we are able to find a trade-off, employing two well-noted VAE variants:
\begin{itemize}
    \item $\beta$-VAE \parencite{Higgins2017betaVAELB}: this variant introduces the hyper-parameter $\beta$ in the variational objective, weighting the KL-term. It results then in:
    $$ \underset{z,\ \phi, \ \theta}{\mathrm{argmin}\ } \log p_\theta(x|z) - \mathbf{\beta} \cdot \mathbb{KL}(q_\phi(z) \| p(z))$$
    The adoption of a $\beta > 1$ should leverage the learning of disentangled representation, by forcing the latent variables in a more constrained/regularized space. As a trade-off we lose in reconstruction capabilities.\\
    In our experiments, although we encounter this effect of better decoupled latent features, the adoption of the $\beta$-VAE results in decreased performances of the subsequent clustering task in the pipeline. Indeed, the strong imposed regularization projects the points in a much more restricted region of the space, having as a side effect a distortion of the clustering validation measures based on the points dispersion
    %good but high regularization given by the $\beta$ term leads to feature mappings in a small region of space, too close 
    \item $\sigma$-VAE \parencite{rybkin2020sigmavae}, an architecture introducing the variance term of the likelihood as a trainable parameter  to \textit{calibrate} (in statistical sense) the decoder distribution. The variational objective is thus modified as follows:
    $$ D \log \sigma + \frac 1{\sigma^2} \log p_\theta(x|z) - \mathbb{KL}(q_\phi(z) \| p(z))$$
    where $D$ is the input-data dimensionality.\\
    Conceptually, the decoder distribution (typically either a multivariate Gaussian or a Multinoulli) is set as:
    $$ p_\theta(x|z) = \mathrm{Dist}(\mu_\theta(z), \sigma^2_\theta(z))$$
    even though, in practice, imposing a common shared diagonal covariance $\mathbf{\Sigma} = \sigma^2 \mathbf{I}$ works better.\\
    Furthermore, it allows to resort to the $\beta$-VAE objective, by setting the variance-term $\sigma^2$ constant and equal to $\frac{\beta}2$. In this way, there is no need to manually tune $\beta$.
\end{itemize}
We reach a noticeable improvement with respect to the PCA reduction, managing to encode the anomaly maps in a 5-dimensional space, ensuring major interpretability and cheaper computations in the subsequent modeling task.

To assess the improvement on the encoded space we rely on a twofold proof. First, the results we obtain for modeling weather regimes (Section \ref{sec:results}) with the VAE algorithm are satisfactory and well in line with the literature, despite a much denser feature space. Thus, it suggests that the VAE algorithm is more efficient to separate the points belonging to different clusters. Secondly, we employ the t-SNE (t-distributed Stochastic Neighbor Embedding) algorithm to find a 2D representation of the high-dimensional manifold of the reduced space, visually proving the enhancements we reach with the VAE (Figure \ref{fig:tsne}). Each point in the chart corresponds to a different anomaly map in our historical dataset, associated to a strong signal of a regime (i.e. where an hard-assignment to a singular regime would be a faithful approximation). The representation is clearly less sparse in terms of intra-cluster variance for the VAE reduction scheme than under the PCA projection. Further, we highlight also how in both the two projections the points associated to AR are typically mixed within the other regimes (more on this in the subsequent Sections).\\

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_PCA10.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_PCA30.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_PCA50.png}
    \end{subfigure}
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_VAE10.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_VAE30.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.32\textwidth}
       \centering
       \includegraphics[width = \linewidth]{imgs/tsne_VAE50.png}
    \end{subfigure}
    \caption{t-SNE applied on the PCA-reduced (top row) and VAE-reduced (bottom row) spaces at different levels of perplxities: the VAE-encoded features leads to a better separability of the clusters, as further proved by the better minimized KL-divergence objective}
    \label{fig:tsne}
\end{figure}

In order to guarantee a truthfulness of these considerations, such visualization is computed at different levels of \textit{perplexity}, representing a parameter of the t-SNE algorithm determining the number of neighboring points used to estimate the projection of the manifold in the bi-dimensional space. To do so, the algorithm minimizes a KL-Divergence objective, whose value we use to further assess the better projection entailed by the VAE reduction.
