Deep Learning for Physical Processes: Incorporating Prior Scientific Knowledge

Introduction

A physical process is a sustained phenomenon marked by gradual changes through a series of states occurring in the physical world. Physicists and environmental scientists attempt to model these processes in a principled way through analytic descriptions of the scientist’s prior knowledge of the underlying processes. Conservation laws, physical principles or phenomenological behaviors are generally formalized using differential equations. This physical paradigm has been, and still is the main framework for modeling complex natural phenomena like e.g. those involved in climate. With the availability of very large datasets captured via different types of sensors, this physical modeling paradigm is being challenged by the statistical Machine Learning (ML) paradigm, which offers a prior-agnostic approach. However, despite impressive successes in a variety of domains as demonstrated by the deployment of Deep Learning methods in fields such as vision, language, speech, etc, the statistical approach is not yet ready to challenge the physical paradigm for modeling complex natural phenomena, or at least it has not demonstrated how to. This is a new challenge for this field and an emerging research direction in the ML community. We believe that knowledge and techniques accumulated for modeling physical processes in well developed fields such as maths or physics could be useful as a guideline to design efficient learning systems and conversely, that the ML paradigm could open new directions for modeling such complex phenomena. In this paper we then raise two issues: 1) are modern ML techniques ready to be used to model complex physical phenomena, and 2) how general knowledge gained from the physical modeling paradigm could help designing efficient ML models.

In this work, we tackle these questions by considering a specific physical modeling problem: forecasting sea surface temperature (SST). SST plays a significant role in analyzing and assessing the dynamics of weather and other biological systems. Accurately modeling and predicting its dynamics is critical in various applications such as weather forecasting, or planning of coastal activities. Since 1982, weather satellites have made huge quantities of very high resolution SST data available (Bernstein, 1982). Standard physical methods for forecasting SST use coupled ocean-atmosphere prediction systems, based on the Navier Stokes equations. These models rely on multiple physical hypotheses and do not optimally exploit the information available in the data. On the other hand, despite the availability of large amounts of data, direct applications of ML methods do not lead to competitive state of the art results, as will be seen in section 4. We use SST as a typical and representative problem of intermediate complexity. Our goal is not to offer one more solution to this problem, but to use it as an illustration for advancing on the challenges mentioned above. The way we handle this problem is general enough to be transfered to a more general class of transport problems.

We propose a Deep Neural Network (NN) model, inspired from general physical motivations which offers a new approach for solving this family of problems. We first motivate our approach by introducing in section 2 the solution of a general class of partial differential equations (PDE) which is a core component of a large family of transport and propagation phenomena in physics. This general solution is used as a guideline for introducing a Deep Learning architecture for SST prediction which is described in section 3. Experiments and comparison with a series of baselines is introduced in section 4. A review of related work is finally presented in section 5.

The main contributions of this work are: 1) an example showing how to incorporate general physical background for designing a NN aimed at modeling a relatively complex prediction task. We believe the approach to be general enough to be used for a family of transport problems obeying general advection-diffusion principles. 2) formal links between our model’s prediction and the solution of a general advection diffusion PDE 3) an unsupervised model for estimating motion fields, given a sequence of images. 4) a proof, on a relatively complex physical modeling problem, that full data intensive approaches based on deep architectures can be competitive with state of the art dedicated numerical methods.

Physical Motivation

Forecasting consists in predicting future temperature maps using past records. Temperatures are acquired via satellite imagery. If we focus on a specific area, we can formulate the problem as prediction of future temperature images of this area using past images. The classical approach to forecasting SST consists in using numerical models representing prior knowledge on the conservation laws and physical principles, which take the form of PDEs. These models are then coupled with SST data using assimilation techniques in order to adjust to initial conditions. It is then integrated forward in time to predict SST evolution. For the sea surface, temperature variation is related to a fluid transport problem. In fluids, transport occurs through the combination of two principles: advection and diffusion. During advection, a fluid transports some conserved quantity $I$ (the temperature for SST) or material via bulk motion, i.e.for small variations $\Delta x$, $\Delta t$ conservation is expressed as:

applying a first order approximation of the right hand side and moving the resulting terms to the left hand side of equation 1, we obtain the advection equation, known also as the Brightness Constancy Constraint Equation (BCCE):

where $\nabla$ denotes the gradient operator, and $w$ the motion vector $\frac{\Delta x}{\Delta t}$. This equation describes the temporal evolution of quantity $I$ for displacement $w$. Note that this equation is also the basis for many variational methods for Optical Flow. To retrieve the motion, numerical schemes are applied, and the resulting system of equations, along with a an additional constraint on $w$ is solved for $w$. This motion can then be used to forecast the future value of $I$

Advection alone is not sufficient to explain the evolution of many physical processes (including SST). Diffusion corresponds to the movement which spreads out the quantity $I$ from areas of high concentration to areas of low concentration. Both advection and diffusion should be considered together. The following equation describes the transport of quantity $I$ through advection and diffusion:

$\nabla^{2}$ denotes the Laplacian operator and $D$ the diffusion coefficient. Note that when $D\!\rightarrow\!0$, we recover the advection equation 2. This equation describes a large family of physical processes (e.g. fluid dynamics, heat conduction, wind dynamics, etc). Let us now state a result, characterizing the general solutions of equation 3.

and $k(u,v)=\frac{1}{4\pi Dt}e^{-\frac{1}{4Dt}\left\|u-v\right\|^{2}}$ is a radial basis function kernel, or alternatively, a Gaussian probability density with mean $x-w$ and variance $2Dt$ in its second argument.

Equation 4 provides a principled way to calculate $I(x,t)$ for any time $t$ using the initial condition $I_{0}$, provided the motion $w$ and diffusion coefficient $D$ are known. It states that quantity $I(x,t)$ can be computed from the initial condition $I_{0}$ via a convolution with a Gaussian probability density function. In other words, if $I$ was used as a model for the evolution of the SST and the surface’s underlying advecting mechanisms were known, future surface temperatures could be predicted from previous ones. Unfortunately neither the initial conditions, the motion vectors nor the diffusion coefficient are known. They have to be estimated from the data. Inspired from the general form of solution 4, we propose a ML method, expressed as a Deep Learning architecture for predicting SST. This model will learn to predict a motion field analog to the $w$ in equation 3, which will be used to predict future images.

Model

The model consists of two main components, as illustrated in Figure 1. One predicts the motion field from a sequence of past input images, this is convolutional-deconvolutional (CDNN) module on the top of figure 1), and the other warps the last input image using the motion field from the first component, in order to produce an image forecast. The entire system is trained in an end-to-end fashion, using only the supervision from the target SST image. By doing so, we are able to produce an interpretable latent state which corresponds in our problem to the velocity field advecting the temperatures.

Generally, and this is the case for our problem, we do not have a direct supervision on the motion vector field, since the target motion is usually not available. Using the warping scheme introduced below, we will nonetheless be able to (weakly) supervise $\hat{w}$, based on discrepancy of the warped version of $I_{t}$ image and the target image $I_{t+1}$.

2 Warping Scheme

Discretizing the solution of the advection-diffusion equation in section 2 by replacing the integral with a sum, and setting image $I_{t}$ as the initial condition, we obtain a method to calculate the future image, based on the motion field estimate $\hat{w}$. The latter is used as a warping scheme:

where $k(x-\hat{w},y)=\frac{1}{4\pi D\Delta t}e^{-\frac{1}{4D\Delta t}\left\|x-\hat{w}-y\right\|^{2}}$ is a radial basis function kernel, as in equation 4, parameterized by the diffusion coefficient $D$ and the time step value $\Delta t$ between $t$ and $t+1$. To calculate the temperature for time $t+1$ at position $x$, we compute the scalar product between $k(x-\hat{w},.)$, a Gaussian centered in $x-\hat{w}$, and the previous image $I_{t}$. Simply put, it is a weighted average of the temperatures $I_{t}$, where the weight values are larger when the pixel’s positions are closer to $x-\hat{w}$. Informally, $x-\hat{w}$ corresponds to the pixel’s previous position at time $t$. See figure 3.

As seen by the relation with the solution of the advection-diffusion equation, the proposed warping mechanism is then clearly adapted to the modeling of phenomena governed by the advection-diffusion equation. SST forecasting is a particular case, but the proposed scheme can be used for any problem problems in which advection and diffusion is occurring. Moreover, this warping scheme is entirely differentiable, allowing backpropagation of the error signal to the motion estimating module.

This warping mechanism has been inspired by the Spatial Transformer Network (STN) (Jaderberg et al., 2015), originally designed to be incorporated as a layer in a convolutional neural network architecture in order to gain invariance under geometric transformations. Using the notations in (Jaderberg et al., 2015), when the inverse geometric transformation $\mathcal{T}_{\theta}$ of the grid generator step is set to $\mathcal{T}_{\theta}(x)=x-\hat{w}(x)$, and the kernels $k(\,.\,;\Phi_{x})$ and $k(\,.\,;\Phi_{y})$ in the sampling step are rbf kernels, we recover our warping scheme. The latter can be seen as a specific case of the STN, without the localization step. This result theoretically grounds the use of the STN for Optical Flow in many recent articles (Zhu et al., 2017), (Yu et al., 2016), (Patraucean et al., 2015), (Finn et al., 2016): in equation 3, when $D\to 0$, we recover the brightness constancy constraint equation, used in the latter.

For training, supervision is provided at the output of the warping module. It consists in minimizing the discrepancy between the warped image $\hat{I}_{t+1}$ and the target image $I_{t+1}$. The loss is measured via a differentiable function and the gradient is back propagated through the warping function in order to adjust the parameters of the convolutional-deconvolutional module generating the vector field. This is detailed in the next section.

3 Loss function

The proposed NN model has been designed according to the intuition gained from general background knowledge of a physical phenomenon, here advection-diffusion equations. Additional prior knowledge – expressed as partial differential equations, or through constraints – can be easily incorporated in our model, by adding penalty terms in the loss function. As the displacement $w$ is explicitly part of our model, one strength of our model is its capacity to apply some regularization term directly on the motion field. In our experiments, we tested the influence of different terms: divergence $\nabla.\,w_{t}(x)^{2}$, magnitude $(\left\|w_{t}(x)\right\|)^{2}$ and smoothness $\left\|\nabla w_{t}(x)\right\|^{2}$. Data assimilation techniques sometimes use these weighted penalties to control the rotation and divergence fields, for example. We evaluate the influence of these terms in the experiments section. The objective function we used to train the model may be written as:

Experiments

Since 1982, high resolution SST data has been made available by the NOAA6 weather satellite (Bernstein, 1982). Dealing directly with these data requires a lot of preprocessing (e.g. some regions are not available due to clouds hindering temperature acquisition). In order to avoid such complications which are beyond the scope of this work, we used synthetic but realistic SST data of the Atlantic ocean generated by a sophisticated simulation engine: NEMO (Nucleus for European Modeling of the Ocean) engine NEMO data are available at http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=GLOBAL_ANALYSIS_FORECAST_PHY_001_024 (Madec, 2008). NEMO is a state-of-the-art modelling framework of ocean related engines. It is a primitive equation model adapted to the regional and global ocean circulation problems. Historical data is accumulated in the model to generate a synthesized estimate of the states of the system using data reanalysis, a specific data assimilation scheme, which means that the data does follow the true temperatures. The resulting dataset is constituted of daily temperature acquisitions of 481 by 781 pixels, from 2006-12-28 to 2017-04-05 (3734 acquisitions).

We extract 64 by 64 pixel sized sub-regions as indicated in figure 4. We use data from years 2006 to 2015 for training and validation (94743 training examples), and years 2016 to 2017 for testing. We withhold 20% of the training data for validation, selected uniformly at random at the beginning of each experiment. For the tests we used sub-regions enumerated 17 to 20 in figure 4, where the interactions between hot and cold waters make the dynamics interesting to study. All the regions numbered in figure 4, from 2006 to 2015 where used for training non numbered regions correspond to land and not sea on the figure. Each sequence of images used for training or for evaluation corresponds to a specific numbered sub-region. We make the simplifying hypothesis that the data in a single sub-region contains enough information to forecast the future of the sub-region. As the forecast is for a small temporal horizon we can assume that the influence from outside the region is small enough.

We standardise the daily SST acquisitions of each sub region using the mean and the standard deviation of all the SST data of the sub region acquired on the same day of the year, i.e. the SST acquisition of sub region 2 on date September 8th 2017 is standardized using the data of all the September 8th available in the dataset, for each sub-region. This removes the seasonal component from SST data.

2 Baseline Comparison

We compare our model with several baselines. Each model is evaluated with a mean square error metric, forecasting images on a horizon of 6 (we forecast from $I_{t+1}$ to $I_{t+6}$ and then average the MSE). The hyperparameters are tuned using the validation set. Neural network based models are run on a Titan Xp GPU, and runtime is given for comparison purpose.

Concerning the constraints on the vector field $w$ (equation 6). the regularization coefficients selected via validation are $\lambda_{\text{div}}=1$, $\lambda_{\text{magn}}=-0.03$ and $\lambda_{\text{grad}}=0.4$. We also compare the results with the model without any regularization.

Our reference model for forecasting is (Béréziat and Herlin, 2015), a numerical assimilation model which relies on data assimilation. In (Béréziat and Herlin, 2015), the ocean’s dynamics are modeled using shallow water equations (Vallis, 2017) and the initial conditions, along with other terms, are estimated using assimilation techniques (Trémolet, 2006). This is a state of the art assimilation model for predicting ocean dynamics, here SST.

The other baselines are 1) an autoregressive convolutional-deconvolutional NN (ACNN), with an architecture similar to our CDNN module, but trained to predict the future image directly, without explicitly representing the motion vector field. Each past observation is used as an input channel, and the output is used as new input for multi step forecasting, 2) a ConvLSTM model (Shi et al., 2015), which uses convolutional transitions in the inner LSTM module, and 3) the model in (Mathieu et al., 2015) which is a multi-scale ACNN trained as a Generative Adversial Network (GAN). We have used a non-official code for (Mathieu et al., 2015), which is made available at https://github.com/dyelax/Adversarial_Video_Generation. For (Béréziat and Herlin, 2015), the code has been provided by the authors of the paper. We have implemented the ACNN and ConvLSTM models ourselves. The code for our models, along with these baselines will be made available.

3 Quantitative Results

Quantitatively, our model performs well. The MSE score is better than any of the baselines. The closest NN baseline is (Mathieu et al., 2015) which regularizes a regression convolution-deconvolution model with a GAN. The performance is however clearly below the proposed model and it does not allow to easily incorporate prior constraints inspired from the physics of the phenomenon. ACNN is a direct predictor of the image sequence, implemented via a CDNN module identical to the one used in our model. Its performance is poor. Clearly, a straightforward use of prediction models is not adapted to the complexity of the phenomenon. ConvLSTM performs better: as opposed to the ACNN, it seems to be able to capture a dynamic, although not very accurately. Overall, direct prediction models are not able to capture the complex underlying dynamics and produce blurry sequences of images. The GAN explicitly forces the network output to eliminate the blurring effect and then makes it able to capture short term dynamics. The state of the art numerical model (Béréziat and Herlin, 2015), performs well but has a lower performance than our regularized model, although it incorporates more prior constraints. This shows that pure ML models, when conceived adequately and when trained with enough data, can be competitive with state of the art dedicated models. Regularizing the motion vector $w$ notably increases the performance w.r.t. to the unregularized model. The choice of the constraints (divergence, magnitude and smoothness) inspired here by physical background correspond to relevant priors on the dynamics of the model.

As for the running time, the proposed model is extremely fast, being just above the ConvLSTM model of (Shi et al., 2015). The running time of (Béréziat and Herlin, 2015)’s model is not comparable to the others. It was run on a CPU (no GPU code) when all the others were run on Titan Xp GPU. However, an optimization procedure is required to estimate the motion field, and it is clearly slower than the straightforward NN predictions. Moreover, in order to prevent the numerical scheme from diverging, multiple intermediate forecasts are required.

Besides MSE, we need to analyze the prediction samples qualitatively. Figure 5 shows the predictions obtained by the different models. The top row is the ground truth for a sequence of 4 temperature images corresponding to time $t$, $t+1$, $t+3$ and $t+6$. The second row corresponds to our regularized model prediction at times $t+1$, $t+3$ and $t+6$ (time $t$ corresponds to the last input image, it is repeated on each row). The model seems to conserve temperatures. The prediction is close to the target for $t+1$, $t+3$ and starts to move away at time $t+6$. The third row shows the motion flow estimated by the model. Each color in the flow images corresponds to a motion vector. There is clearly a strong evolving dynamic captured for this sequence. Row 4 is the numerical assimilation model of (Béréziat and Herlin, 2015). It also clearly captures some dynamics and shows interesting patterns, but it tends to diverge when the prediction horizon increases. The ACNN model (row 5) rapidly produces blurry images; it does not preserve the temperatures and does not seem to capture any dynamics. Row 6 shows the predictions of the ConvLSTM model. Temperature is not preserved and although a dynamic is captured, it does not correspond to the target. Overall, the proposed model seems to forecast SST quite accurately, while retrieving a coherent motion vector field.

Related Work

ML for Physical modeling Close to this work is the field of spatio-temporal statistics. In their reference book (Cressie and Wikle, 2015) also advocate the use of physical background knowledge to build statistical models. They show how the design of statistical models can be inspired from partial differential equations linked to an observed physical phenomenon. They mainly consider auto-regressive models within a hierarchical Bayesian framework. Another interesting research direction is the use of NNs for reducing the complexity of numerical simulation for physical processes. Generally, in these approaches statistical models are used in place of a computational demanding component of the numerical simulation process. For example in the domain of fluid dynamics, (Tompson et al., 2017) and (Ladický et al., 2015) propose to use regressors for simulating fluid and smoke animation. (Ladický et al., 2015) use a random forest to compute particle location and (Tompson et al., 2017) use a CNN to approximate part of a numerical PDE scheme. In these approaches, ML is only a component of a numerical simulation scheme whereas we aim at modeling the whole physical process via a Deep Learning approach. Farther to our objective, (Rudy et al., 2017) make use of a sparse regression method for discovering the governing partial differential equation(s) of a given system by time series measurements in the spatial domain.

Our work is also related to recent developments in computer vision, in the related but distinct fields of video prediction and motion estimation in videos. Our goal and the domain of application are clearly different from video modeling, but since our solution involves predicting a motion field and the next SST image, the solutions share some similarities. Motion estimation and video predictions by deep architectures have motivated a series of work over the last two years. We briefly review them below and outline the differences.

Optical Flow Optical flow consists in retrieving the apparent motion of objects, surfaces, or particles between two consecutive frames of a video. The extracted motion can be used in many areas such as object detection, object tracking, movement detection, robot navigation and visual odometry. In the vision community, this is considered as a problem by itself and several papers are dedicated to this topic. Classical methods rely on the brightness constancy constrain equation (BCCE) (equation 2), derived from the observation that surfaces usually persist over time and hence the intensity value of a small region remains the same despite its position change (Sun et al., 2008). Since using BCCE directly leads to complicated optimizing issues, classic approaches – namely differential methods – approximate BCCE using a first order Taylor expansion and develop variational methods.

As an alternative to these methods, Deep Learning models have been recently proposed for estimating the optical flow between 2 images. (Fischer et al., 2015) formulate optical flow as a supervised regression problem, using a CNN to predict motion. (Ilg et al., 2016) build on this approach and propose to use an ensemble of these CNN architectures. They assess results on par with state of the art methods for optical flow, while maintaining a small computational overhead. The difficulty here is that these methods require a notable quantity of target data, i.e. optical flow images, while because of the complexity of manually annotating flow images, there are only a few small annotated collections available. (Fischer et al., 2015) and (Ilg et al., 2016) chose to pretrain their model on a synthetic dataset made of computer animations and their associated motion and show that this information transfers well to real videos. (Yu et al., 2016) demonstrate that it is possible to predict the optical flow between two input images in an unsupervised way using a CNN and a Spatial Transformer Network. This is however not extensible for prediction as is done in our setting since this requires the two images $I_{t}$ and $I_{t+1}$ as input while $I_{t+1}$ is not available at inference time for prediction.

It is only very recently that video prediction emerged as a task in the Deep Learning community. For this task, people are generally interested at predicting accurately the displacement/ emergence/ disappearing of objects in the video. In our application, the goal is clearly different since we are interested into modeling the whole dynamics behind image changes and not at following moving objects. Let us first introduce some methods that perform prediction by computing optical flow or a similar transformation. Both (Patraucean et al., 2015) and (Finn et al., 2016) use some form of motion flow estimation. For next frame prediction (Patraucean et al., 2015) introduce a STN module at the hidden layer of a LSTM in order do estimate a motion field in this latent space. The resulting image is then decoded in the original image space for prediction. This approach clearly does not allow introducing prior knowledge on the field vector as this has been done in our work. (Finn et al., 2016) learn affine transformations on image parts in order to predict object displacement and (Van Amersfoort et al., 2017) proposed a similar model.

Let us now consider models that directly attempt to predict the next frame without estimating a motion field. As shown in the experimental section, plain application of autoregressive models produces blurred images. (Mathieu et al., 2015), one of our baseline proposed to use different loss functions and a GAN regularization of a CDNN predictor which led to sharper and higher quality predictions. Significant improvements have been obtained with the Video Pixel Network of (Kalchbrenner et al., 2016), which is a sophisticated architecture composed of resolution preserving CNN encoders, LSTM and PixelCNN decoders which form a conditional Spatio-temporal video autoencoder with differentiable memory. This model is probably state of the art today for video prediction, They reach a high accuracy on moving MNIST and good performance on a robot video dataset. A drawback is the complexity of the model and the number of parameters: they are using respectively 20 M and 1 M frames on these two datasets. We did not test this model since up to our knowledge no code was available.

Conclusion

The data intensive paradigm offers alternative directions to the classical physical approaches for modeling complex natural processes. Our belief is that cross fertilization of both paradigms is essential for pushing further the frontier of complex data modeling. By using as an example application a problem of intermediate complexity concerning ocean dynamics, we proposed a principled way to design Deep Learning models using inspiration from the physics. The proposed approach can be easily generalized to a class of problems for which the underlying dynamics follow advection-diffusion principles. We have compared the proposed approach to a series of baselines. It is able to reach performance comparable to a state of the art numerical model and clearly outperforms alternative NN models used as baselines.

This work was partially funded by ANR project LOCUST - ANR-15-CE23-0027 and by CLEAR Lab.

References

Appendix A Proof of the theorem in section 2 1

We apply the Fourier Transform $\mathcal{F}$ to both sides of 3. As consequence of the linearity of the Fourier transform, we can calculate decompose the Fourier transform of the left hand side in the sum of the transforms of each term. We have three terms: $\frac{\partial I}{\partial t}$, $(w.\nabla)I$ and $-D\nabla^{2}I$.

Regrouping all three previously calculated terms, we obtain

This is a first order ordinary differential equation of the form $f^{\prime}(t)+af(t)=0$, which admits a known solution $f(t)=f(0)e^{-at}$. Thus, the solution of 11 is

where $\mathcal{F}(I)_{0}$ denotes the initial condition of the advection diffusion equation in the frequency domain. In order to obtain a solution of 3 in the spatial domain, we calculate the inverse Fourier Transform $\mathcal{F}^{-1}$ of 12. The multiplication of two functions in the frequency domain is equivalent to their convolution in the spatial domain, i.e. $\mathcal{F}(f*g)=\mathcal{F}(f)\mathcal{F}(g)$. Hence, the inverse of both terms $\mathcal{F}(I)_{0}\,e^{-i<\xi,w>t}$ and $e^{-Dt\left\|\xi\right\|^{2}}$ can be calculated separately:

For the second term, we use the fact that the Fourier Transform of a Gaussian function also is a Gaussian function, i.e. $\mathcal{F}(\frac{1}{2\pi\sigma^{2}}e^{-\frac{1}{2\sigma^{2}}\left\|\mathbf{x}\right\|^{2}})=e^{-\frac{1}{2}\sigma^{2}\left\|\xi\right\|^{2}}$. Identifying $\sigma^{2}$ with $2Dt$, we have:

As has been stated above, the solution is a convolution of both previously calculated terms:

Appendix B Additional samples from our model