Title: LEARNING TEMPORAL EVOLUTION OF PROBABILITY DISTRIBUTION WITH RECURRENT NEURAL NETWORK Anonymous

Abstract: We propose to tackle a time series regression problem by computing temporal evolution of a probability density function to provide a probabilistic forecast. A Recurrent Neural Network (RNN) based model is employed to learn a nonlinear operator for temporal evolution of a probability density function. We use a softmax layer for a numerical discretization of a smooth probability density functions, which transforms a function approximation problem to a classification task. Explicit and implicit regularization strategies are introduced to impose a smoothness condition on the estimated probability distribution. A Monte Carlo procedure to compute the temporal evolution of the distribution for a multiple-step forecast is presented. The evaluation of the proposed algorithm on three synthetic and two real data sets shows advantage over the compared baselines.

1. INTRODUCTION:

Application of the deep learning for manufacturing processes has attracted a great attention as one of the core technologies in Industry 4.0 (Lasi et al., 2014) . In many manufacturing processes, e.g. blast furnace, smelter, and milling, the complexity of the overall system makes it almost impossible or impractical to develop a simulation model from the first principles. Hence, system identification from sensor observations has been a long-standing research topic (Wang et al., 2016) . Still, when the observation is noisy and there is no prior knowledge on the underlying dynamics, there is only a very limited number of methods for the reconstruction of nonlinear dynamics.

In this work, we consider the following class of problems, where the system is driven by a complex underlying dynamical system, e.g., ∂y ∂t = F(y(t), y(tτ ), u(t)).

(1)

Here, y(t) is a continuous process, F is a nonlinear operator, τ is a delay-time parameter, and u(t) is an exogenous forcing, such as control parameters. At time step t, we then observe a noisy measurement of y(t) which can be defined by the following noise model

EQUATION

where ν t is a multiplicative and t is an additive noise process. In (1) and ( 2), we place no assumption on function F, do not assume any distributional properties of noises ν t and t , but assume the knowledge of the control parameters u(t).

Since the noise components, ν t and t , are stochastic processes, the observation ŷt is a random variable. In this work, we are interested in computing temporal evolution of the probability density function (PDF) of ŷ, given the observations up to time step t, i.e., p(ŷ t+n | Y 0:t , U 0:t+n-1 ) for n ≥ 1, where Y 0:t = (ŷ 0 , • • • , ŷt ) is a trajectory of the past observations and U 0:t+n-1 = (u 0 , • • • , u t+n-1 ) consists of the history of the known control actions, U 0:t-1 , and a future control scenario, U t:t+n-1 . We show, in Section 3, a class of problems, where simple regression problem of forecasting the value of ŷt+n is not sufficient or not possible, e.g., chaotic systems. Note that the computation of time evolution of a PDF has been a long-standing topic in statistical physics. For a simple Markov process, there are well-established theories based on the Fokker-Planck equation. However, it is very difficult to extend those theories to a more general problem, such as delay-time dynamical systems, or apply it to complex nonlinear systems.

Modeling of the system (1) has been extensively studied in the past, in particular, under the linearity assumptions on F and certain noise models, e.g., Gaussian t and ν t = 1 in (2). The approaches based on auto-regressive processes (Lütkepohl, 2005) and Kalman filter (Harvey, 1990) are good examples. Although these methods do estimate the predictive probability distribution and enable the computation of the forecast uncertainty, the assumptions on the noise and linearity in many cases make it challenging to model real nonlinear dynamical systems.

Recently, a nonlinear state-space model based on the Gaussian process, called the Gaussian Process State Space Model (GPSSM), has been extended for the identification of nonlinear system (Frigola et al., 2013; Eleftheriadis et al., 2017) . GPSSM is capable of representing a nonlinear system and is particularly advantageous when the size of the data set is relatively small that it is difficult to train a deep learning model. However, the joint Gaussian assumption of GPSSM may restrict the representation capability for a complex non-Gaussian noise.

A recent success of deep learning created a flurry of new approaches for time series modeling and prediction. The ability of deep neural networks, such as RNN, to learn complex nonlinear spatiotemporal relationships in the data enabled these methods to outperform the classical time series approaches. For example, in the recent works of Qin et al. (2017) ; Hsu (2017) ; Dasgupta & Osogami (2017) , the authors proposed different variants of the RNN-based algorithms to perform time series predictions and showed their advantage over the traditional methods. Although encouraging, these approaches lack the ability to estimate the probability distribution of the predictions since RNN is a deterministic model and unable to fully capture the stochastic nature of the data.

To enable RNN to model the stochastic properties of the data, Chung et al. (2015) augmented RNN with a latent random variable included in the hidden state and proposed to estimate the resulting model using variational inference. In a similar vein, the works of Archer et al. (2015) ; Krishnan et al. (2017) extend the traditional Kalman filter to handle nonlinear dynamics when the inference becomes intractable. Their approach is based on formulating the variational lower bound and optimizing it under the assumption of Gaussian posterior.

Another recent line of works enabled stochasticity in the RNN-based models by drawing a connection between Bayesian variation inference and a dropout technique. In particular, Gal & Ghahramani (2016) showed that the model parameter uncertainty (which then leads to uncertainty in model predictions), that traditionally was estimated using variational inference, can be approximated using a dropout method (a random removal of some connections in the network structure). The prediction uncertainty is then estimated by evaluating the model outputs at different realizations of the dropout weights. Following the ideas of Gal & Ghahramani (2016) , Zhu & Laptev (2017) proposed additional ways (besides modeling the parameter uncertainty) to quantify the forecast uncertainty in RNN, which included the model mis-specification error and the inherent noise of the data.

1.1. OVERVIEW OF THE PROPOSED WORK:

We propose an RNN-model to compute the temporal evolution of a PDF, p(ŷ t+n | Y 0:t , U 0:t+n-1 ).

To avoid the difficulties in directly estimating the continuous function, we use a numerical discretization technique, which converts the function approximation problem to a classification task (see Section 2.2). We note that the use of the traditional cross-entropy (CE) loss in our formulated classification problem can be problematic since it is oblivious to the class ordering. We empirically show that the proposed modeling approach enables us to represent a continuous PDF of any arbitrary shape, including the ability to handle the multiplicative data noises in (2). Since the probability distribution is computed, the RNN-model can also be used for a regression task by computing the expectation (see Section 2.4). Hereafter, we use DE-RNN for the proposed RNN model, considering the similarity with the density-estimation task.

In summary, the contributions of this work are as follows: (i) formulate the classical regression problem for time series prediction as a predictive density-estimation problem, which can be solved by a classification task (ii) propose an approach to compute the time evolution of probability distribution using SMC on the distributions from DE-RNN (iii) proposed two regularizations for CE loss to capture the ordering of the classes in the discretized PDF. We evaluated the proposed algorithm on three synthetic and two real datasets, showing its advantage over the baselines. Note that DE-RNN has a direct relevance to a wide range of problems in physics and engineering, in particular, for uncertainty quantification and propagation (Zhang & Karniadakis, 2017) .

2. LSTM FOR NOISY DYNAMICAL SYSTEM:

In this Section we present the details of the proposed algorithm using a specific form of RNN, called Long Short-Term Memory (LSTM) network. Although in the following presentation and experiments we used LSTM, other networks, e.g., GRU (Chung et al., 2014) , can be used instead.

2.1. REVIEW OF LONG SHORT-TERM MEMORY NETWORK:

The Long Short-Term Memory network (LSTM) (Hochreiter & Schmidhuber, 1997; Gers et al., 2000) consists of a set of nonlinear transformations of input variables z t ∈ R m ;

Gating functions:

G i,f,o = ϕ S • L(z t ), (3) Internal state: s t = (1 -G f ) s t-1 + G i (ϕ T • L(z t )) , (4) Output: h t = G o s t .

(5)

Here, ϕ S and ϕ T , respectively, denote the sigmoid and hyperbolic tangent functions, L is a linear layer, which includes a bias, s t ∈ R Nc is the internal state, h t ∈ R Nc is the output of the LSTM network, N c is the number of the LSTM units, and a b denote a component-wise multiplication.

Interesting observation can be made about equation ( 4). We can re-write equation (4) as

EQUATION

for some functions f and g. With a simple re-scaling, this equation can be interpreted as a first-order Euler scheme for a linear dynamical system, ds dt = -f (z)s + g(z).

Thus, LSTM can be understood as a series expansion, where a complex nonlinear dynamical system is approximated by a combination of many simpler dynamical systems.

Usually, LSTM network is supplemented by feed-forward neural networks, e.g.,

EQUATION

in which x t is the input feature. Using (5), we can denote by Ψ e and Ψ d a collection of the operators from input to internal state (encoder) and from internal state to the output P (decoder):

EQUATION

2.2. DISCRETE APPROXIMATION OF PROBABILITY DENSITY FUNCTION:

In this Section we first consider the problem of modeling the conditional PDF, p(ŷ t+1 | Y 0:t , U 0:t ).

Although ŷt+1 has a dependence on the past trajectories of both ŷ and u, using the "state space" LSTM model argument in Section 2.1, the conditional PDF can be modeled as a Markov process

EQUATION

Hence, to simplify the problem, we consider a task of estimating the PDF of a random variable ŷ, given an input x, i.e., p(ŷ|x). The obtained results can then be directly applied to the original problem of estimating p(ŷ t+1 |s t ).

Let α = (α 0 , • • • , α K ) denote a set of real numbers, such that α i-1 < α i for i = 1, • • • , K, which defines K disjoint intervals, I i = (α i-1 , α i ).

Then, a discrete probability distribution can be defined

EQUATION

where it is clear that p(k|x) is a numerical discretization of the continuous PDF, p(ŷ|x). Using the LSTM from Section 2.1, the discrete probability p(k|x) can be modeled by the softmax layer (P ) as an output of Ψ d in ( 9) such that

EQUATION

Thus, the original problem of estimating a smooth function, p(ŷ|x), is transformed into a classification problem of estimating p(k|x) in a discrete space. Obviously, the size of the bin, |I j |, affects the fidelity of the approximation. The effects of the bin size are presented in Section 3.1. There is a similarity between the discretization and the idea of Lin et al. (2007) . However, it should be noted that the same discretization technique, often called "finite volume method", has been widely used in the numerical simulations of partial differential equations for a long time.

The discretization naturally leads to the conventional cross-entropy (CE) minimization. Suppose we have a data set, D R = {(ŷ i , x i ); ŷi ∈ R, x i ∈ R, and i = 1, . . . , N }. We can define a mapping

C : R → N + such that C(ŷ) = k, if y ∈ I k . Then, D R can be easily converted to a new data set for target labels, D C = {(c i , ŷi , x i ); c i ∈ N + , ŷi ∈ R, x i ∈ R, and i = 1, . . . , N }, where c i = C(ŷ i ).

D C provides a training data set for the following CE minimization problem,

EQUATION

Note, however, that the CE minimization does not explicitly guarantee the smoothness of the estimated distribution. Since CE loss function depends only on P i of a correct label, δ cnk , as a result, in the optimization problem every element P i , except for the one corresponding to the correct label, P cn , is penalized in the same way, which is natural in the conventional classification tasks where a geometric proximity between the classes is not relevant. In the present study, however, the softmax layer, or class probability, is used as a discrete approximation to a smooth function. Hence, it is expected that P cn and P cn±1 (i.e., the nearby classes) should be close to each other. To address this issue, in the following Sections 2.2.1 and 2.2.2, we propose two types of regularization to impose the class proximity structure in the CE loss.

2.2.1. EXPLICIT REGULARIZATION OF CROSS-ENTROPY LOSS:

To explicitly impose the smoothness between the classes, we propose to use a regularized crossentropy (RCE) minimization, defined by the following loss function

EQUATION

)

where λ is a penalty parameter and the Laplacian matrix

EQUATION

)

RCE is analogous to the penalized maximum likelihood solution for density estimation (Silverman, 1986) . Assuming a uniform bin size,

|I 0 | = • • • = |I K | = δy,

the Laplacian of a distribution can be approximated by a Taylor expansion p (ŷ|x)

| y=α i-1/2 (P i-1 -2P i + P i+1 )/δy 2 , where α i-1/2 = 0.5(α i-1 + α i ).

Then, it is clear that

EQUATION

In other words, RCE aims to smooth out the distribution by penalizing local minima or maxima.

2.2.2. IMPLICIT REGULARIZATION OF CROSS-ENTROPY LOSS:

Alternative to adding an explicit regularization to CE, the smoothness can be achieved by enforcing a spatial correlation in the network output. Here, we use an one-dimensional convolution layer to enforce smoothness. Let o ∈ R K denote the last layer of DE-RNN, which was the input to the softmax layer. We can add a convolution layer, o ∈ R K , on top of o, such that

EQUATION

where the parameter h determines the smoothness of the estimated distribution. Then, o is supplied to the softmax layer. Using ( 17), DE-RNN can now be trained by the standard CE. The implicit regularization, here we call convolution CE (CCE), is analogous to a kernel density estimation.

2.3. MULTIVARIATE TIME SERIES:

In the modeling of multivariate time series, it is usually assumed that the noise is independent, i.e., the covariance of the noise is a diagonal matrix. In this case, it is straightforward to extend DE-RNN, by using multiple softmax layers as the output of DE-RNN. However, such an independent noise assumption significantly limits the representative capability of an RNN. Here, we propose to use a set of independently trained DE-RNNs to compute the joint PDF of a multivariate time series.

Let ŷt be a l-dimensional multivariate time series; ŷt = (ŷ

(1) t , • • • , ŷ(l) t ).

The joint PDF can be represented by a product rule,

p(ŷ t+1 ) = p ŷ(l) t+1 |ŷ (l-1) t+1 , • • • , ŷ(1) t+1 p ŷ(l-1) t+1 |ŷ (l-2) t+1 , • • • , ŷ(1) t+1 • • • p ŷ(2) t+1 |ŷ

(1)

t+1 p ŷ(1) t+1 ,

where the dependency on the past trajectory ( Y 0:t , U 0:t ) is omitted in the notation for simplicity. Directly learning the joint PDF, p(ŷ t+1 | Y 0:t , U 0:t ), in a tensor product space is not scalable. 

2.4. COMPUTING TIME EVOLUTION OF PROBABILITY DISTRIBUTION:

The inputs to a DE-RNN are (ŷ t , u t ), and the output is the probability distribution,

P t+1 = Ψ d (s t ) = Ψ d • Ψ e (ŷ t , u t , s t-1 ).

Note that D C is used only in the training stage. Then, the moments of the predictive distribution can be easily evaluated, e.g.,

EQUATION

) where α 1/2 = (α 1/2 , α 1+1/2 , • • • , α K-1/2 ) T , α 2 1/2 = α 1/2 α 1/2

, and α i-1/2 = 0.5(α i-1 + α i ). Next, we consider a multiple-step forecast, which corresponds to computing a temporal evolution of the probability distribution, i.e., p(ŷ t+n | Y 0:t , U 0:t+n-1 ) for n > 1. For simplicity, the multiple-step forecast is shown only for a univariate time series. An extension to a multivariate time series is straightforward (Appendix A).

Applying the results of Section 2.2, once the distribution of ŷt+1 in (10) is computed, the distribution of ŷt+2 can be similarly obtained as p(ŷ t+2 |s t+1 ). Observe that s t+1 is computed from a deterministic function of s t , u t+1 , and ŷt+1 , i.e., s t+1 = Ψ e (ŷ t+1 , u t+1 , s t ).

Here, u t+1 and s t are already known, while ŷt+1 is a random variable, whose distribution p(ŷ t+1 |s t ) is computed from the deterministic function Ψ d (s t ). Then, s t+1 is also a random variable. The distribution, p(s t+1 |s t , u t+1 ), can be obtained by applying a change of variables on p(ŷ t+1 |s t ) with a nonlinear mapping Ψ e . Repeating this process, the multiple-step-ahead predictive distribution can therefore be computed as

EQUATION

)

Since the high dimensional integration in ( 19) is intractable, we propose to approximate it by a sequential Monte Carlo method. The Monte Carlo procedure is outlined in Algorithm 1.

Algorithm 1 Sequential Monte Carlo method for LSTM multi-step-ahead prediction

Input: Y 0:t , U 0:t , number of Monte Carlo samples, N s , and forecast horizon n Output: p(ŷ t+n | Y 0:t , U 0:t+n-1 ) (density estimation from ŷt+n ) Initialization: Set LSTM states to s 0 = h 0 = 0 Perform a sequential update of LSTM up to time t from the noisy observations ( Y 0:t ).

s i = Ψ e (ŷ i , u i , s i-1 ) for i = 1, • • • , t.

Make N s replicas of the internal state,

s 1 t = • • • = s Ns t = s t . repeat

Compute the predictive distribution of ŷi t+1 for each sample

P i t+1 = Ψ d (s i t ), for i = 1, • • • , N s .

Sample the target variable at t + 1, ŷi t+1 , from the distribution 1. Sample the class label from the discrete distribution:

c i ∼ P i t+1 2. Sample ŷi t+1 in I c i : ŷi t+1 ∼ U(α c i -1 , α c i ) Update the internal state of LSTM s i t+1 = Ψ e (ŷ i t+1 , u t+1 , s i t )

. until (all ŷt+n are sampled)

3. EXPERIMENTS:

In this section, DE-RNN is tested against three synthetic and two real data sets. The LSTM architecture used in all of the numerical experiments is identical. Two feed-forward networks are used before and after the LSTM;

EQUATION

in which ϕ SP and ϕ SM denote the softplus and softmax functions, respectively. The size of the bins is kept uniform, i.e.,

|I 1 | = • • • = |I K | = δy.

The LSTM is trained by using ADAM (Kingma & Ba, 2015) with a minibath size of 20 and a learning rate of η = 10 -3 .

3.1. COX-INGERSOLL-ROSS PROCESS:

First, we consider a modified Cox-Ingersoll-Ross (CIR) process, which is represented by the following stochastic differential equation,

EQUATION

in which W is the Weiner process. The original CIR process is used to model the valuation of interest rate derivatives in finance. Equation ( 21) is solved by the forward Euler method with the time step size δt = 0.1. The simulation is performed for T = (0, 160000]δt to generate the training data and T = (160000, 162000]δt is used for the testing. Note that the noise component of CIR is multiplicative, which depends on y(t). The experiments are performed for two different bin sizes, dy = 0.08 and 0.04. The DE-RNN has 64 LSTM cells. Figure 1 shows the errors in the expectation and the standard deviation with respect to the analytical solution;

EQUATION

Here, p T denotes the true distribution of the CIR process. The normalized root mean-square errors (NRMSE) are defined as

EQUATION

)

EQUATION

in which • denotes an average over the testing data, p L is the distribution from the LSTM, and sd[y] denotes the standard deviation of the data. The error in the expectation is normalized against a zeroth-order prediction, which assumes y t+1 = y t .

In Figure 1 , it is clearly shown that the prediction results are improved when a regularization is used to impose a smoothness condition. Comparing Figures 1 (a ) and (b), for RCE, e µ and e σ become smaller when a smaller δy is used. As expected, e σ increases when λ is large. But, for the smaller bin size, δy = 0.04, both e µ and e σ are not so sensitive to λ. Similar to RCE, e µ and e σ for CCE decrease at first as the penalty parameter h increases. However, in general, RCE provides a better prediction compared to CCE.

NRMEs are listed in Table 1 . For a comparison, the predictions by AR(1) and KF are shown. The CIR process is essentially a first-order autoregressive process. So, it is not surprising to see that AR(1) and KF, which are designed for the first-order AR process, outperforms DE-RNN for the prediction of the expectation. However, e σ of AR(1) and KF are much larger than that of DE-RNN, because those models assume an additive noise. The Gaussian process (GP) model has a relatively large e µ . But, GP outperforms AR(1) and KF in the prediction of the noise (e σ ). Still, e σ of RCE and CCE are less than 4%, while that of GP is about 10%, indicating that DE-RNN can model the complex noise process much better.

In Figure 2 The SMC forecast can accurately capture the behavior. Such kind of behavior can not be represented if a simple additive noise is assumed.

3.2. MACKEY-GLASS TIME SERIES:

For the next test, we applied DE-RNN for a time series generated from the Mackey-Galss equation (Mackey & Glass, 1977) ;

EQUATION

)

We use the parameters adopted from Gers ( 2001), α = 0.2, β = 10, γ = 0.1, and τ = 17.

The Mackey-Glass equation is solved by using a third-order Adams-Bashforth method with a time step size of 0.02. The time series is generated by down-sampling, such that the time interval between consecutive data is δt = 1. A noisy observation is made by adding a white noise;

ŷt = y t + t .

t is a zero-mean Gaussian random variable with the noise level sd[ t ] = 0.3sd [y] . A time series of the length 1.6 × 10 5 δt is generated for the model trainig and another time series of length 2 × 10 3 δt is made for the validation. DE-RNN is trained for δy = 0.04sd[y] and consists of 128 LSTM cells. The prediction errors are shown in table 2. NRMSEs are defined as,

EQUATION

NRMSEs are computed with respect to the ground truth. Again, e µ compares the prediction error to the zeroth-order prediction. In this example, the errors are not so sensitive to the regularization parameters. The best result is achieved by CCE. DE-RNN can make a very good estimation of the noise. The error in the noise component, e σ , is only 2% ∼ 5%. Unlike the CIR process, NRMSEs from KF and ARIMA are much larger than those of DE-RNN. Because the underlying process is a delay-time nonlinear dynamical system, those linear models can not accurately approximate the complex dynamics. Since GP is capable of representing a nonlinear behavior of data, GP outperforms KF and ARIMA both in e µ and e σ . In particular, e σ of GP is similar to that of DE-RNN.

However, e µ of GP is about 1.5 times larger than DE-RNN.

A multiple-step forecast of the Mackey-Glass time series is shown in Figure 4 . In the validation time series, the observations in t ∈ [1, 100]δt are supplied to the DE-RNN to develop the internal state, and a 500-step forecast is made for t ∈ [101, 600]δt. In Figure 4 (a), it is shown that a multiple-step forecast by a standard regression LSTM approximates y(t) very well initially, e.g, for t < 80δt, but eventually diverges for larger t. Because of the Mackey-Glass time series is chaotic, theoretically it is impossible to make a long time forecast. But, in the DE-RNN forecast, y(t) is bounded by the 95%-confidence interval (CI) even for the 500-step forecast. Note that the uncertainty, denoted by 95%-CI grows only at a very mild rate in time. In fact, it is observed that CI is not a monotonic function of time. In DE-RNN, the 95%-CI may grow or decrease following the dynamics of the system, while for the conventional time series models, such as ARIMA and KF, the uncertainty is a non-decreasing function of time.

3.3. MAUNA LOA CO 2 OBSERVATION:

In this experiments, DE-RNN is tested against the atmospheric CO 2 observation at Mauna Loa Observatory, Hawaii (Keeling et al., 2001) . The CO 2 data set consists of weekly-average atmospheric CO 2 concentration from Mar-29-1958 to Sep-23-2017 (Figure 5 a) . The data from Mar-29-1958 to Apr-01-2000 is used to train DE-RNN and a 17-year forecast is made from Apr-01-2000 to Sep-23-2017. This CO 2 data has been used in previous studies (Gal & Ghahramani, 2016; Rasmussen & Williams, 2006) . In DE-RNN, 64 LSTM cells and δy = 0.1sd[dy t ], in which dy t = y t+1y t , are used. The 17-year DE-RNN forecast, with 1,000 MC samples, is shown in Figure 5 (b). DE-RNN well represents the growing trend and the oscillatory patten of the CO 2 data. The CO 2 data is nonstationary, where the rate of increase of CO 2 is an increasing function of time. Since DE-RNN is trained against the history data, where the rate of CO 2 increase is smaller than current, it is expected that the forecast will underestimate the future CO 2 . E[ŷ n+s |ŷ s ] agrees very well with the observation for the first 200 weeks, but eventually underestimates CO 2 concentration. It is interesting to observe that the upper bound of the 95%-CI grows more rapidly than the expectation and provides a good approximation of the observation. For a comparison, the forecast by a regression LSTM is also shown. Similar to the chaotic Mackey-Glass time series, the regression LSTM makes a good prediction for a short time, e.g., t < 100 weeks, but eventually diverges from the observation. Note that the lower bound of 95%-CI encompasses the regression LSTM. Figure 5 (c) shows a forecast made by GP, following setup suggested by Rasmussen & Williams (2005) . For a shortterm forecast (< 50 weeks), GP provides a sharper estimate of the uncertainty, i.e., a smaller 95%-CI interval. However, even for a mid-term forecast, 100 ∼ 600 weeks, the ground truth is close or slightly above the upper bound of 95%-CI. Because of the different behaviors, it is difficult to conclude which method, DE-RNN or GP, provides a better forecast. But, it should be noted that the GP is based on a combination of handcrafted kernels specifically designed for this particular problem, while such careful tuning is not required for DE-RNN.

3.4. CPU TEMPERATURE FORECAST:

In the last experiment, IBM Power System S822LC and NAS Parallel Benchmark (NPB) are used to generate the temperature trace. Figure 6 (a) shows the temperature of a CPU. The temperature sensor generates a discrete data, which has a resolution of 1 • C. The CPU temperature is controlled by three major parameters; CPU frequency, CPU utilization, and cooling fan speed. In this experiment, we have randomized the CPU frequencies and job arrival time to mimic the real workload behavior, while the fan speed is fixed to 3300RPM. The time step size is δt = 2 seconds. Accurate forecast of CPU temperature for a future workload scenario is essential in developing an energy-efficient control strategy for the thermal management of a cloud system. Figure 6 (c) and (d) show multiple-step forecasts of the CPU temperature by RCE and a regression LSTM, respectively. The bin size is δy = 0.18 • C, which is smaller than the sensor resolution. In the forecast, the future control parameters are given to DE-RNN. In other words, DE-RNN predicts the probability distribution of future temperature with respect to a control scenario, i.e., p(ŷ t+n | Y 0:t , U 0:t+n-1 ). The forecast is made by using 5,000 Monte Carlo samples. Here, 1800-step forecast is made, t = 0 ∼ 3, 600 sec. and only the results in t ∈ (50, 1800) sec. is shown. While the regression LSTM makes a very noisy prediction near the local peak temperature at t 500, RCE provides a much more stable forecast. Table 3 shows the l ∞ -errors, i.e., maximum absolute Figure 4 shows normalized errors for the next-step prediction. The normalized errors are defined as

e µ = 1 3 3 i=1 (E[ŷ (i) t+1 ] -y (i) t+1 ) 2 (ŷ (i) t -y (i) t+1 ) 2 , e σ = 1 3 3 i=1 V ar[ŷ (i) t+1 ] σ 2 i -1, e Σij = Cov[ŷ (i) t+1 , ŷ(j) t+1 ] Σ ij -1.

(29) The moments of the joint PDF are computed by the Monte Carlo method with a sample size of 5 × 10 3 . It is shown that DE-RNN makes a very good prediction of both expectations and covariances. The error in the covariance is less than 4% except for Σ 12 . It is shown that DE-RNN is able to make correct predictions of not only the magnitude, but also the signs of the covariance. A vector autoregressive model (VAR) also predicts the signs of the covariance correctly. But, the errors in the expectation and covariances are much larger than DE-RNN. The GP used in the experiment assumes an independent noise, i.e. Σ = ρ 2 I. Hence, e Σij is not evaluated for GP. Similar to the Mackey-Glass time series, GP outperforms VAR, but the errors are larger than DE-RNN. The error in e σ is about 10 times larger than DE-RNN, while e µ of GP is about 3 times larger than DE-RNN.

Figure 7 shows 300-step forecasts of the Lorenz time series. The expectations from DE-RNNs make a good prediction of the ground truth up to ts < 1.5. Then, the expectations start to diverge from the ground truth, which is accompanied by a sudden increase in the 95%-CI. For a longer-time forecast, e.g., ts > 4, the 95%-CI exhibits an oscillatory patten depending on the dynamics of the Lorenz system and well captures the oscillations in the ground truth.

4. CONCLUDING REMARKS:

We present DE-RNN to compute the time evolution of a probability distribution for complex time series data. DE-RNN employs LSTM to learn multiscale, nonlinear dynamics from the noisy observations, which is supplemented by a softmax layer to approximate a probability density function. To assign probability to the softmax output, we use a mapping from R to N + , which leads to a cross-entropy minimization problem. To impose a geometric structure in the distribution, two regularization strategies are proposed. The regularized cross-entropy method is analogous to the penalized maximum likelihood estimate for the density estimation, while the convolution cross-entropy method is motivated by the kernel density estimation. The proposed algorithm is validated against three synthetic data set, for which we can compare with the analytical solutions, and two real data sets. 

:



Figures and Tables:

Instead, a set of DE-RNN is trained to represent the conditional PDFs shown on the right hand side of the above expression. Then, the joint PDF can be computed by a product of the Softmax outputs of the DE-RNNs. Note that, although it requires training l DE-RNNs to compute the full joint PDF, there is no dependency between the DE-RNNs in the training phase. So, the set of DE-RNNs can be trained in parallel, which can greatly reduce the training time. The details of the multivariate DE-RNN are explained in Appendix A.

Figure 1: NRMSE of the next-step prediction of the CIR process by RCE (a,b) and CCE (c,d). The bin size is (a,c) δy = 0.08 and (b,d) 0.04. The hollow circles (•) denote NRMSE in the expectation (e µ ) and the solid circles (•) are the standard deviation (e σ ).

, a 200-step forecast by DE-RNN is compared with a Monte-Carlo solution of equation (21). DE-RNN is trained with δy = 0.04 and λ = 200. For the DE-RNN forecast, the testing data is supplied to DE-RNN for the first 100 time steps, i.e., for t = -10 to t = 0, and the SMC multiple-step forecast is performed for the next 200 time steps with 20,000 samples. It is shown that the multiple-step forecast by DE-RNN agrees very well with the MC solution of the CIR process. Note that, in Figure2(b), the noise process, as reflected in sd[y t ], is a function of y t , and hence the multi-step forecast of the noise increases rapidly first and then decreases before reaching a plateau.

Figure 2: 200-step forecast of (a) expectation and (b) standard deviation of the CIR process. The circles denote the solution of Eqn (21) from a Monte Carlo method with 10 7 samples.

Figure 3: (a) Next-step prediction (•) and the noisy observation (•) for the Mackey-Glass equation. The solid line denotes the ground truth, y(t). (b) The next-step probability distribution, p(ŷ t+1 |ŷ t ), from the standard CE (•) and CCE (•) with h = 5.

Figure 3 (a) shows the noisy observation and the expectation of the next-step prediction, E[ŷ t+1 |ŷ t ], in a phase space. It is shown that DE-RNN can filter out the noise and reconstruct the original dynamics accurately. Even though the noisy data are used as an input, E[ŷ t+1 |ŷ t ] accurately represents the original attractor of the chaotic system, indicating a strong de-noising capability of DE-RNN.The estimated probability distribution is shown in Figure3 (b). Without a regularization, the standard CE results in a noisy distribution, while the distribution from CCE shows a smooth Gaussian shape.

Figure 5: (a) Mauna Loa CO 2 observation. The vertical line denotes the boundary of the training and testing data. (b) 17-year forecast of the CO 2 concentration by DE-RNN: Apr-01-2000 Sep-23-2017. (c) 17-year forecast by Gaussian processes. The solid and dashed lines denote the expectation and 95%-CI, respectively. The observation is shown as the solid circles (•) and a regression LSTM is shown as the hollow circles (•) in (b). The time unit is a week.

Figure 7: Multiple-step forecast of the Lorenz time series by DE-RNN trained with CCE. The solid and dashed lines denote the expectation and 95%-CI, respectively. The ground truth is shown as sold circles.

Figure 8: Architecture of the DE-RNN algorithm.

Figure 9: Details of the computation for the one-step-ahead predictions. At a given step the model computes a discrete approximation of the expectation for the next observation.

Figure 10: Visualization of the multi-step-ahead forecast made by DE-RNN. The shown illustration is for computing 3-step-ahead predictive probability distribution p(ŷ t+3 | 0:t , U 0:t+2 ). N s independent DE-RNN instances, i = 1, . . . , N s , are initialized identically with the hidden states from the time step of the last observation. Each DE-RNN instance i uses a two-stage sampling procedure to propagate its own value of ŷi in time. At the end, N s samples of ŷi t+3 are used to compute a discrete approximation of the predictive distribution p(ŷ t+3 | Y 0:t , U 0:t+2 ).

NRMSE of the mean and standard deviation of the next-step prediction. The DE-RNN results are compared with the first-order autoregressive model (AR), Kalman filter (KF), and Gaussian process (GP).

Table content: <table><tr><td>CE</td><td>RCE</td><td>CCE AR(1)</td><td>KF</td><td>GP</td></tr><tr><td colspan="5">e µ 0.238 0.0549 0.149 0.029 0.029 0.831</td></tr><tr><td colspan="5">e σ 0.066 0.017 0.038 0.228 0.228 0.095</td></tr></table>

l ∞ error for RCE, CCE and regression LSTM in • C. RCE CCE LSTM E[ŷ n+s |ŷ s ] -ŷn+s ∞ 4.61 2.70 9.48

Table content: <table/>

Normalized errors of the Lorenz time series. DE-RNN results are compared with vector autoregressive (VAR) and Gaussian process (GP) models.

Table content: <table><tr><td/><td>e µ</td><td>e σ</td><td>e Σ12</td><td>e Σ13</td><td>e Σ23</td></tr><tr><td>CE</td><td colspan="5">0.19 0.027 0.146 -0.012 0.024</td></tr><tr><td colspan="6">RCE 0.20 0.030 0.155 -0.010 0.026</td></tr><tr><td colspan="6">CCE 0.20 0.018 0.123 -0.009 0.022</td></tr><tr><td colspan="6">VAR 1.00 0.318 1.485 0.483 0.649</td></tr><tr><td>GP</td><td colspan="2">0.58 0.296</td><td>-</td><td>-</td><td>-</td></tr></table>

Peer review (iclr2018):
 Interesting ideas that extend LSTM to produce probabilistic forecasts for univariate time series, experiments are okay. Unclear if this would work at all in higher-dimensional time series. It is also unclear to me what are the sources of the uncertainties captured. The author proposed to incorporate 2 different discretisation techniques into LSTM, in order to produce probabilistic forecasts of univariate time series. The proposed approach deviates from the Bayesian framework where there are well-defined priors on the model, and the parameter uncertainties are subsequently updated to incorporate information from the observed data, and propagated to the forecasts. Instead, the conditional density p(y_t|y_{1:t-1|, .theta}) was discretised by 1 of the 2 proposed schemes and parameterised by a LSTM. The LSTM was trained using discretised data and cross-entropy loss with regularisations to account for ordering of the discretised labels. Therefore, the uncertainties produced by the model appear to be a black-box. It is probably unlikely that the discretisation method can be generalised to high-dimensional setting? Quality: The experiments with synthetic data sufficiently showed that the model can produce good forecasts and predictive standard deviations that agree with the ground truth. In the experiments with real data, it*s unclear how good the uncertainties produced by the model are. It may be useful to compare to the uncertainty produced by a GP with suitable kernels. In Fig 6c, the 95pct CI looks more or less constant over time. Is there an explanation for that? Clarity: The paper is well-written. The presentations of the ideas are pretty clear. Originality: Above average. I think the regularisation techniques proposed to preserve the ordering of the discretised class label are quite clever. Significance: Average. It would be excellent if the authors can extend this to higher dimensional time series. I*m unsure about the correctness of Algorithm 1 as I don*t have knowledge in SMC.