---
abstract: 'We explore the use of deep neural networks for nonlinear dimensionality reduction in climate applications. We train convolutional autoencoders (CAEs) to encode two temperature field datasets from pre-industrial control runs in the CMIP5 first ensemble, obtained with the CCSM4 model and the IPSL-CM5A-LR model, respectively. With the later dataset, consisting of 36500 96$\times$96 surface temperature fields, the CAE out-performs PCA in terms of mean squared error of the reconstruction from a 40 dimensional encoding. Moreover, the noise in the filters of the convolutional layers in the autoencoders suggests that the CAE can be trained to produce better results. Our results indicate that convolutional autoencoders may provide an effective platform for the construction of surrogate climate models.'
author:
- 'Juan A. Saenz$^{1}$[^1] [^2], Nicholas Lubbers$^{1,2}$, Nathan M. Urban$^{1}$'
title: 'Dimensionality-Reduction of Climate Data using Deep Autoencoders'
---

Introduction
============

Uncertainty quantification of the response of the Earth system to greenhouse-gas emission scenarios is important for evaluating the impacts of climate change on infrastructure, agriculture, and the environment, among other areas. However, simulations using global, coupled earth system models are computationally expensive, making it impossible to produce large ensembles needed for statistical uncertainty quantification. To overcome this, surrogate models–simplified models that emulate more complex climate models–are built and trained [@little_etal_2013; @little_etal_2013b]. The computational cost of running these climate model emulators is much lower than their full complexity counterparts. As a result, large ensembles of simulations (order tens-hundreds of thousands) can be produced and used to carry out uncertainty quantification.

Emulators can be devised to dynamically evolve the state of the climate on a dimensionally reduced manifold [@ross_2009]. An important requirement of such emulators is that the state in physical dimensions can be recovered. Linear dimensionality reduction via principal component analysis (PCA) is well-known in the climate science community, however, nonlinear methods have not been fully explored. Ross [@ross_2009] investigated nonlinear dimensionality reduction methods for a different climate aplication: identifying low-dimensional nonlinear dynamics in El Nino variability. Errors from reconstructions using nonlinear methods were not significantly better than using linear PCA. Methods include nonlinear PCA (autoencoders), Isomap, and Hessian locally linear embedding. Here, motivated by the success of deep autoencoders for dimensionality reduction [@hinton_salakhutdinov_2006] and convolutional neural networks for image processing [@alexnet; @vggnet], we present work in progress using convolutional autoencoders [@convautoencoder] to reduce the dimensionality of data from climate models. In section \[s:methods\] we describe our methods, and in section \[s:results\] we present results on two pre-industrial climate model simulation datasets: the CCSM4-T31 temperature at the surface dataset, and the IPSL-CM5A-LR temperature at the surface data. We end with a brief discussion in section \[s:discussion\].

Methods {#s:methods}
=======

Principal Components Analysis {#s:pca}
-----------------------------

Consider a dataset of dimensionality $M$ with $N$ datapoints collected into a $N \times M$ data matrix $\mathbf{X}$. PCA constructs a rank $m$ reduced matrix $ \hat{\mathbf{X}}_{pca}$ by projecting $\tilde{ \mathbf{X} }$ (obtained by centering and normalizing $\mathbf{X}$ using the global mean and standard deviation) onto the first $m$ principal components which maximizes the covariance of the data, thus minimizing the mean reconstruction error $$\label{eq:mse}
MSE = \frac{1}{NM} \sum_n{(\mathbf{x}_n-\mathbf{\hat{x}}_n)^2} =\frac{1}{NM} {{\lVert\mathbf{X}-\hat{\mathbf{X}}\rVert}_2^2}.$$ This can be obtained by singular value decomposition of $\mathbf{X}$ = $\mathbf{U}\mathbf{\Sigma}\mathbf{W}^T$ with $\mathbf{\Sigma}$ the diagonal matrix of singular values. The data covariance is given by $\mathbf{X}^T \mathbf{X}$ = $\mathbf{W} \mathbf{\Sigma}^2 \mathbf{W}^T$, and so the first $m$ principal component vectors are the columns of $\mathbf{W}$ associated with the $m$ largest singular values.

Convolutional Autoencoder {#s:convauto}
-------------------------

The autoencoder provides an alternative method for dimensionality reduction of $\mathbf{X}$. Data is fed though a series of neural network layers, i.e. an affine transform followed by an elementwise nonlinearity $f$, to create activations $\mathbf{X}_l$ at layer $l$: $$\mathbf{X}_l = f(\mathbf{W}_l \mathbf{X}_{l-1}+\mathbf{b}_l).$$ The trainable parameters of the network are the weights $\mathbf{W}_l$ and biases $\mathbf{b}_l$ of each layer. The waist layer of the autoencoder is constrained to a number of neurons $m$, so that the activations at that layer can be used as an $m$ length code for each image. This is followed by decoding layers, and in the final layer of the autoencoder, $\mathbf{\hat{X}}$ is constructed to have the same dimensionality as the input $\mathbf{X}$. The parameters of the network are then trained to minimize the same reconstruction error $MSE$ (Eq. \[eq:mse\]) as PCA.

The convolutional autoencoder (CAE) is an extension to the autoencoder which faciltates the analysis of data on regular grids, such as images. Here, each data point can be indexed by two pixel positions. In the convolutional and deconvolutional layers, the weights consist of small image filter kernels, and the product of the weights and the data consist of spatial convolutions. Denoting the feature $k$ of the pixel indexed by $i$ and $j$ by $\mathbf{x}_{i,j}^k$, the preactivations $\mathbf{W}\mathbf{x}$ are computed using $$(\mathbf{W} \mathbf{x})_{i,j}^l = \sum_{k,l,m}{\mathbf{W}^{lk}_{{i-l},{j-m}} \mathbf{x}^k_{i,j}} \,.$$ Convolutional layers operate as a set of local image filters with the capacity to extract patterns that increase in complexity with depth [@convpatterns]. They have the benefit of greatly reducing the number of learnable parameters per layer, and faciliate training by seeking out features which encode the structure of local image patches [@lecun98].

Finally, CAEs have pooling layers that coarse-grain the image plane after each convolution. We use $2\times2$ max pooling, with unpooling layers which use piece-wise constant $2\times2$ upsampling. We use two regularization strategies to improve training regularity and smoothness of learned filters. In the first, we use image flipping and weight decay with strength $\beta$. In the second, we use noise injection (denoising autoencoder), applying pixel-wise Gaussian noise with mean $0$ and standard deviation $\gamma$ to the images.

Autoencoders presented here are trained using stochastic gradient descent for 1000 epochs with a learning rate of 0.01, which is updated using a Nesterov scheme with a momentum of 0.975, and a batch size of 128. Convolutional layers use linear activation functions (we have not been able to successfully train non-linear activations). Fully connected layers use rectified linear (ReLU) activation functions. Decoding layers have the reverse structure of the encoding layers, but we do not tie the weights between the encoding and decoding layers. Convolutional boundary conditions use valid convolutions (no padding), with decoding convolutional layer dimensions computed to produce the correct size output. As preprocessing, the temperature fields are normalized to mean $0$ and standard deviation $1$ using the global mean and standard deviation.

\[t:architectures\]

   label    encoding architecture     $\beta$   $\gamma$
  ------- -------------------------- --------- ----------
    A1          CL7,32-PL-FC40        0.00025     N/A
    A2     CL5,32-PL-CL5,64-PL-FC40   0.00025     N/A
    B1     CL5,32-PL-CL5,64-PL-FC40     N/A       N/A
    B2     CL5,32-PL-CL5,64-PL-FC40     N/A       0.1
    B3     CL5,32-PL-CL5,64-PL-FC40     N/A       0.5

  : Architectures of Convoluational Autoencoders

CL$n$,$m$ is a convolutional layer with an $n\times n$ receptive field and $m$ features, PL is a 2$\times$2 pooling layer, and FC$m$ is a fully connected layer with $m$ neurons.

Results {#s:results}
=======

CCSM4 {#s:CCSM4_T31}
-----

We train encoders with data from the Community Climate System Model - version 4 (CCSM4) [@gent_etal_2011]. We employ 150 years of surface temperature ($T_s$) monthly climatology data from the first ensemble, pre-industrial control run of the 5th version of the Climate Model Intercomparison Project (CMIP5). The data was regridded to a 3.75$\times$3.75 degree grid (T31 grid). The resulting dataset has 1800 samples of 48$\times$96 $T_s$ fields.

The trained architectures (A1 and A2, table \[t:architectures\]) do a good job at recovering the global structure of the temperature field $T_s$. There are differences in local features in some regions (not shown). In table \[t:mse\] we compare the $MSE$ (Eq. \[eq:mse\]) of reconstructions using CAEs A1 and A2 and PCA. Architecture A1 performs better than A2, but PCA performed better than both autoencoders. The weights of the first convolutional layer indicated that this type of regularization was not effective (Fig. \[f:weights\], upper left). There is some repetitive structure between weights, and they are noisy, indicating that the neural networks are not well trained. One possible reason for this is that there are not enough samples ($N$=1800). To investigate this, we now use a dataset with a much larger number of samples.

![Weights from the first convolutional layer ($\mathbf{W}^1$) in the trained architectures A1 (top left) and B1 (top right), B2 (bottom left) and B3 (bottom right). Brightness scale is arbitrary, only the patterns visible in the plots are meaningful here.[]{data-label="f:weights"}](weigths_T31_l12.pdf "fig:"){width="0.3\hsize"} ![Weights from the first convolutional layer ($\mathbf{W}^1$) in the trained architectures A1 (top left) and B1 (top right), B2 (bottom left) and B3 (bottom right). Brightness scale is arbitrary, only the patterns visible in the plots are meaningful here.[]{data-label="f:weights"}](IPSL_l25_weights_conv1.pdf "fig:"){width="0.3\hsize"}\
![Weights from the first convolutional layer ($\mathbf{W}^1$) in the trained architectures A1 (top left) and B1 (top right), B2 (bottom left) and B3 (bottom right). Brightness scale is arbitrary, only the patterns visible in the plots are meaningful here.[]{data-label="f:weights"}]({IPSL_l26_01_weights_conv1}.pdf "fig:"){width="0.3\hsize"} ![Weights from the first convolutional layer ($\mathbf{W}^1$) in the trained architectures A1 (top left) and B1 (top right), B2 (bottom left) and B3 (bottom right). Brightness scale is arbitrary, only the patterns visible in the plots are meaningful here.[]{data-label="f:weights"}]({IPSL_l26_02_weights_conv1}.pdf "fig:"){width="0.3\hsize"}\

IPSL-CM5A-LR {#s:IPSL}
------------

Here, we use near-surface air temperature $T_{as}$ from the low resolution pre-industrial control run produced by the IPSL-CM5A-LR model, part of the CMIP5 first ensemble [@dufresne_etal_2013]. The data is on a 1.9$\times$3.75 degree grid, 96$\times$96, and we use daily output for 100 years between 1800 and 1900 ($N$=36500 samples), although 600 years of data are available (219000 samples). As a result, the dataset has 36500 samples of 96$\times$96 $T_{as}$ fields, resulting in $\mathbf{X}$ with shape $(N=36500) \times (M=9216)$.

We implement an autoencoder architecutre, B1, which is similar to A2, but without any type of regularization (table \[t:architectures\]). The MSE using B1, shown in table \[t:mse\], is smaller than the error obtained with PCA. However, the weights are still noisy (Figure \[f:weights\], upper right). To remedy this, we explored architectures B2 and B3 (table \[t:architectures\]) which are regularized using injected noise. The weights become smoother with more noise (Fig \[f:weights\], lower left and right), but the errors are larger, as shown in table \[t:mse\].

In figure \[f:IPSL\] we show the reconstructed temperature fields, which are very similar to the temperature in the original dataset. Large scale features of the global temperature patterns are preserved. Smaller scale features in regions such as over the Antarctic peninsula, the North Atlantic and the South Pacific are filtered out. Some small scale features above high elevation topography, such as the Andes and the Himalayas, appear to be well preserved.

\[t:mse\]

  -- -------- -------- -------- --------
       PCA       A1       A2    
      0.8022   1.5495   2.1608  
       PCA       B1       B2       B3
      4.9014   4.2415   4.3806   4.5243
  -- -------- -------- -------- --------

  : Mean squared error of reconstructions \[K$^2$\].

\

Discussion {#s:discussion}
==========

The results of this work in progress indicate that there is potential to devise deep autoencoders for dimensionality reduction of climate data. Noise in the filters of the trained networks indicate that finding effective representations is dependent on regularization and the availability of data.

Future developments will focus on testing networks with data not used for training/validating, using larger datasets, using nonlinear convolutional activation functions, implementing other regularization methods (e.g. dropout), and using deeper networks, with the aim of improving the reconstruction of small scale features. Future analysis will include the investigation of patterns extracted by the convolutions.

Acknowledgments {#acknowledgments .unnumbered}
===============

This work was supported by the Office of Science (BER), U. S. Department of Energy. Autoencoders are built using Python libraries Theano [@theano_2016], Lasagne [@sander_dieleman_2015_27878], and nolearn [@nolearn]. We thank A. Jonko for providing the CCSM4-T31 data. We acknowledge the WCRP-WGCM, which is responsible for CMIP. We thank the DOE and NCAR for developing the CCSM4, producing and making available their model output.

[1]{}

C. M. Little, M. Oppenheimer, and N. M. Urban, “Upper bounds on twenty-first-century antarctic ice loss assessed using a probabilistic framework,” [*Nature Clim. Change*]{}, vol. 3, pp. 654–659, 07 2013.

C. M. Little, N. M. Urban, and M. Oppenheimer, “Probabilistic framework for assessing the ice sheet contribution to sea level change,” [*Proceedings of the National Academy of Sciences*]{}, vol. 110, no. 9, pp. 3264–3269, 2013.

I. [Ross]{}, [*[Nonlinear Dimensionality Reduction Methods in Climate Data Analysis]{}*]{}. PhD thesis, University of Bristol, Jan. 2009.

G. E. Hinton and R. R. Salakhutdinov, “Reducing the dimensionality of data with neural networks,” [*Science*]{}, vol. 313, no. 5786, pp. 504–507, 2006.

A. Krizhevsky, I. Sutskever, and G. E. Hinton. “ImageNet classification with deep convolutional neural networks”. *Adv. Neural Inf. Process. Syst.*, pages 1097–1105, 2012.

Simonyan, Karen, and Andrew Zisserman. “Very deep convolutional networks for large-scale image recognition." *arXiv preprint* arXiv:1409.1556 (2014).

Masci, Jonathan, et al. “Stacked convolutional auto-encoders for hierarchical feature extraction." *International Conference on Artificial Neural Networks.* Springer Berlin Heidelberg, 2011.

Lee, Honglak, et al. “Convolutional deep belief networks for scalable unsupervised learning of hierarchical representations." *Proceedings of the 26th annual international conference on machine learning.* ACM, 2009.

LeCun, Yann, et al. “Gradient-based learning applied to document recognition." *Proceedings of the IEEE* 86.11 (1998): 2278-2324.

P.R. Gent, G. Danabasoglu, L.J. Donner, M.M. Holland, E.C. Hunke, S.R. Jayne, D.M. Lawrence, R.B. Neale, P.J. Rasch, M. Vertenstein, P.H. Worley, Z. Yang, and M. Zhang, “The Community Climate System Model Version 4”. [*J. Climate*]{}, vol. 24, pp. 4973–4991.

J.-L. Dufresne, M.-A. Foujols, S. Denvil, A. Caubel, O. Marti, O. Aumont, Y. Balkanski, S. Bekki, H. Bellenger, R. Benshila, S. Bony, L. Bopp, P. Braconnot, P. Brockmann, P. Cadule, F. Cheruy, F. Codron, A. Cozic, D. Cugnet, N. de Noblet, J.-P. Duvel, C. Eth[é]{}, L. Fairhead, T. Fichefet, S. Flavoni, P. Friedlingstein, J.-Y. Grandpeix, L. Guez, E. Guilyardi, D. Hauglustaine, F. Hourdin, A. Idelkadi, J. Ghattas, S. Joussaume, M. Kageyama, G. Krinner, S. Labetoulle, A. Lahellec, M.-P. Lefebvre, F. Lefevre, C. Levy, Z. X. Li, J. Lloyd, F. Lott, G. Madec, M. Mancip, M. Marchand, S. Masson, Y. Meurdesoif, J. Mignot, I. Musat, S. Parouty, J. Polcher, C. Rio, M. Schulz, D. Swingedouw, S. Szopa, C. Talandier, P. Terray, N. Viovy, and N. Vuichard, “Climate change projections using the ipsl-cm5 earth system model: from cmip3 to cmip5,” [*Climate Dynamics*]{}, vol. 40, no. 9, pp. 2123–2165, 2013.

R. Al[-]{}Rfou, G. Alain, A. Almahairi et al. 2016 “Theano: A Python framework for fast computation of mathematical expressions” [*CoRR*]{}, abs/1605.02688.

S. Dieleman, J. Schlüter, C. Raffel, E. Olson, S.K. Sønderby, D. Nouri, D. Maturana, M. Thoma, E. Battenberg, J. Kelly, J. De Fauw, M. Heilman, dogo149, B. McFee, H. Weideman, takacsg84, peterderivaz, Jon, instagibbs, K. Rasul, CongLiu, Britefury, J. Degrave. 2015 “Lasagne: First release” doi: 10.5281/zenodo.27878

D. Nouri. 2014. “nolearn: scikit-learn compatible neural network library”. https://github.com/dnouri/nolearn

[^1]: Corresponding author: J.A. Saenz, juan.saenz@lanl.gov; N. Lubbers, nlubbers@bu.edu; N.M. Urban, nurban@lanl.gov.

[^2]: $^1$Los Alamos National Laboratory, $^2$Boston University.
