Title: LEARNING AND POLICY SEARCH IN STOCHASTIC DYNAMICAL SYSTEMS WITH BAYESIAN NEURAL NETWORKS

Abstract: We present an algorithm for policy search in stochastic dynamical systems using model-based reinforcement learning. The system dynamics are described with Bayesian neural networks (BNNs) that include stochastic input variables. These input variables allow us to capture complex statistical patterns in the transition dynamics (e.g. multi-modality and heteroskedasticity), which are usually missed by alternative modeling approaches. After learning the dynamics, our BNNs are then fed into an algorithm that performs random roll-outs and uses stochastic optimization for policy learning. We train our BNNs by minimizing α-divergences with α = 0.5, which usually produces better results than other techniques such as variational Bayes. We illustrate the performance of our method by solving a challenging problem where model-based approaches usually fail and by obtaining promising results in real-world scenarios including the control of a gas turbine and an industrial benchmark.

1. INTRODUCTION:

In model-based reinforcement learning, an agent uses its experience to first learn a model of the environment and then uses that model to reason about what action to take next. We consider the case in which the agent observes the current state s t , takes some action a, and then observes the next state s t+1 . The problem of learning the model corresponds then to learning a stochastic transition function p(s t+1 |s t , a) specifying the conditional distribution of s t+1 given s t and a. Most classic control theory texts, e.g. Bertsekas (2002) , will start with the most general model of dynamical systems:

s t+1 = f (s t , a, z, W)

where f is some deterministic function parameterized by weights W that takes as input the current state s t , the control signal a, and some stochastic disturbance z.

However, to date, we have not been able to robustly learn dynamical system models to such a level of generality. Popular modes for transition functions include Gaussian processes (Rasmussen et al., 2003; Ko et al., 2007; Deisenroth & Rasmussen, 2011) , fixed bases such as Laguerre functions (Wahlberg, 1991) , and adaptive basis functions or neural networks (Draeger et al., 1995) . All of these methods assume deterministic transition functions, perhaps with some addition of Gaussian observation noise. Thus, they are severely limited in the kinds of stochasticity-or transition noise-they can express. In many real-world scenarios stochasticity may often arise due to some unobserved environmental feature that can affect the dynamics in complex ways (such as unmeasured gusts of wind on a boat).

In this work we use Bayesian neural networks (BNNs) in conjunction with a random input noise source z to express stochastic dynamics. We take advantage of a very recent inference advance based on α-divergence minimization (Hernández-Lobato et al., 2016) , with α = 0.5, to learn with high accuracy BNN transition functions that are both scalable and expressive in terms of stochastic patterns. Previous work achieved one but not both of these two characteristics.

We focus our evaluation on the off-policy batch reinforcement learning scenario, in which we are given an initial batch of data from an already-running system and are asked to find a better (ideally near-optimal) policy. Such scenarios are common in real-world industry settings such as turbine control, where exploration is restricted to avoid possible damage to the system. We propose an algorithm that uses random roll-outs and stochastic optimization for learning an optimal policy from the predictions of BNNs. This method produces (to our knowledge) the first model-based solution of a 20-year-old benchmark problem: the Wet-Chicken (Tresp, 1994) . We also obtain very promising results on a real-world application on controlling gas turbines and on an industrial benchmark.

2. BACKGROUND:



2.1. MODEL-BASED REINFORCEMENT LEARNING:

We consider reinforcement learning problems in which an agent acts in a stochastic environment by sequentially choosing actions over a sequence of time steps, in order to minimize a cumulative cost. We assume that our environment has some true dynamics T true (s t+1 |s, a), and we are given a cost function c(s t ). In the model-based reinforcement learning setting, our goal is to learn an approximation T approx (s t+1 |s, a) for the true dynamics based on collected samples (s t , a, s t+1 ). The agent then tries to solve the control problem in which T approx is assumed to be the true dynamics.

2.2. BAYESIAN NEURAL NETWORKS WITH STOCHASTIC INPUTS:

Given data D = {x n , y n } N n=1 , formed by feature vectors x n ∈ R D and targets y n ∈ R K , we assume that y n = f (x n , z n ; W) + n , where f (•, •; W) is the output of a neural network with weights W. The network receives as input the feature vector x n and the random disturbance z n ∼ N (0, γ). The activation functions for the hidden layers are rectifiers: ϕ(x) = max(x, 0). The activation functions for the output layers are the identity function: ϕ(x) = x. The network output is corrupted by the additive noise variable n ∼ N (0, Σ) with diagonal covariance matrix Σ. The role of the noise disturbance z n is to capture unobserved stochastic features that can affect the network's output in complex ways. Without z n , randomness is only given by the additive Gaussian observation noise n , which can only describe limited stochastic patterns. The network has L layers, with V l hidden units in layer l, and W = {W l } L l=1 is the collection of V l × (V l-1 + 1) weight matrices. The +1 is introduced here to account for the additional per-layer biases.

One could argue why n is needed at all when we are already using the more flexible stochastic model based on z n . The reason for this is that, in practice, we make predictions with the above model by averaging over a finite number of samples of z n and W. By using n , we obtain a predictive distribution whose density is well defined and given by a mixture of Gaussians. If we eliminate n , the predictive density is degenerate and given by a mixture of delta functions.

Let Y be an N × K matrix with the targets y n and X be an N × D matrix of feature vectors x n . We denote by z the N -dimensional vector with the values of the random disturbances z 1 , . . . , z N that were used to generate the data. The likelihood function is

p(Y | W, z, X) = N n=1 p(y n | W, z, x n ) = N n=1 K k=1 N (y n,k | f (x n , z n ; W), Σ) .

(1)

The prior for each entry in z is N (0, γ). We also specify a Gaussian prior distribution for each entry in each of the weight matrices in W. That is,

EQUATION

where w ij,l is the entry in the i-th row and j-th column of W l and γ and λ are a prior variances. The posterior distribution for the weights W and the random disturbances z is given by Bayes' rule:

EQUATION

Variational Bayes 0.5 1 0 q tends to fit a mode of p q tends to fit p globally Figure 1 : Solution for the minimization of the α-divergence between the posterior p (in blue) and the Gaussian approximation q (in red and unnormalized). Figure source Minka (2005) .

Given a new input vector x , we can then make predictions for y using the predictive distribution

EQUATION

Unfortunately, the exact computation of ( 4) is intractable and we have to use approximations.

2.3. α-DIVERGENCE MINIMIZATION:

We approximate the exact posterior distribution p(W, z | D) with the factorized Gaussian distribution

EQUATION

The parameters m w ij,l , v w ij,l and m z n , v z n are determined by minimizing a divergence between p(W, z | D) and the approximation q. After fitting q, we make predictions by replacing p(W, z | D) with q in (4) and approximating the integrals in (4) with empirical averages over samples of W ∼ q.

We aim to adjust the parameters of (5) by minimizing the α-divergence between p(W, z | D) and q(W, z) (Minka, 2005) :

EQUATION

which includes a parameter α ∈ R that controls the properties of the optimal q. Figure 1 illustrates these properties for the one-dimensional case. When α ≥ 1, q tends to cover the whole posterior distribution p. When α ≤ 0, q tends to fit a local mode in p. The value α = 0.5 is expected to achieve a balance between these two tendencies. Importantly, when α → 0, the solution obtained is the same as with variational Bayes (VB) (Wainwright & Jordan, 2008) .

The direct minimization of ( 6) is infeasible in practice for arbitrary α. Instead, we follow Hernández-Lobato et al. (2016) and optimize an energy function whose minimizer corresponds to a local minimization of α-divergences, with one α-divergence for each of the N likelihood factors in (1). Since q is Gaussian and the priors p(W) and p(z) are also Gaussian, we represent q as q(W, z)

EQUATION

where f (W) is a Gaussian factor that approximates the geometric mean of the N likelihood factors in (1) as a function of W. Each f n (z n ) is also a Gaussian factor that approximates the n-th likelihood factor in (1) as a function of z n . We adjust f (W) and the f n (z n ) by minimizing local α-divergences.

In particular, we minimize the energy function -Lobato et al., 2016) , where f (W) and f n (z n ) are in exponential Gaussian form and parameterized in terms of the parameters of q and the priors p(W) and p(z n ), that is,

EQUATION

EQUATION

EQUATION

and log Z q is the logarithm of the normalization constant of the exponential Gaussian form of q:

EQUATION

The scalable optimization of ( 8) is done in practice by using stochastic gradient descent. For this, we subsample the sums for n = 1, . . . , N in ( 8) and ( 11) using mini-batches and approximate the expectations over q in (8) with an average over K samples drawn from q. We can then use the reparametrization trick (Kingma et al., 2015) to obtain gradients from the resulting stochastic approximator to (8). The hyper-parameters Σ, λ and γ can also be tuned by minimizing ( 8). In practice we only tune Σ and keep λ = 1 and γ = d. The latter means that the prior scale of each z n grows with the data dimensionality. This guarantees that, a priori, the effect of each z n in the neural network's output does not diminish when more and more features are available.

Minimizing (8) when α → 0 is equivalent to running the method VB (Hernández-Lobato et al., 2016) , which has recently been used to train Bayesian neural networks in reinforcement learning problems (Blundell et al., 2015; Houthooft et al., 2016; Gal et al., 2016) . However, we propose to minimize (8) using α = 0.5, which often results in better test log-likelihood values.

We have also observed α = 0.5 to be more robust than VB when q(z) is not fully optimized. In particular, α = 0.5 can still capture complex stochastic patterns even when we do not learn q(z) and instead keep it fixed to the prior p(z). By contrast, VB fails completely in this case (see Appendix A).

3. POLICY SEARCH USING BNNS WITH STOCHASTIC INPUTS:

We now describe a gradient-based policy search algorithm that uses the BNNs with stochastic disturbances from the previous section. The motivation for our approach lies in its applicability to industrial systems: we wish to estimate a policy in parametric form, using only an available batch of state transitions obtained from an already-running system. We assume that the true dynamics present stochastic patterns that arise due to some unobserved process affecting the system in complex ways.

Model-based policy search methods include two key parts (Deisenroth et al., 2013) . The first part consists in learning a dynamics model from data in the form of state transitions (s t , a t , s t+1 ), where s t denotes the current state, a t is the action applied and s t+1 is the resulting state. The second part consists in learning the parameters W π of a deterministic policy function π that returns the optimal action a t = π(s t ; W π ) as function of the current state s t . The policy function can be a neural network with deterministic weights given by W π .

The first part in the aforementioned procedure is a standard regression task, which we solve by using the modeling approach from the previous section. We assume the dynamics to be stochastic with the following true transition model:

s t = f true (s t-1 , a t-1 , z t ; W true ) , z t ∼ N (0, γ) .

(12) where the input disturbances z t ∼ N (0, γ) account for the stochasticity in the dynamics. When the Markov state s t is hidden and we are given only observations o t , we can use the time embedding theorem using a suitable window of length n and approximate:

ŝ(t) = [o t-n , • • • , o t ] .

(13) The transition model in equation 12 specifies a probability distribution p(s t |s t-1 , a t-1 ) that we approximate using a BNN with stochastic inputs:

EQUATION

Algorithm 1 Model-based policy search using Bayesian neural networks with stochastic inputs. 1: Input: D = {s n , a n , ∆ n } for n ∈ 1..N 2: Fit q(W) and Σ by optimizing (8).

3: function UNFOLD(s 0 ) 4:

sample{W 1 , .., W K } from q(W)

5:

C ← 0 6:

for k = 1 : K do 7:

for t = 0 : T do 8:

z k t+1 ∼ N (0, γ) 9: ∆ t ← f (s t , π(s t ; W π ), z k t+1 ; W k ) 10: k t+1 ∼ N (0, Σ)

11: where the feature vectors in our BNN are now s t-1 and a t-1 and the targets are given by s t . In this expression, the integration with respect to W accounts for stochasticity arising from lack of knowledge of the model parameters, while the integration with respect to z t accounts for stochasticity arising from unobserved processes that cannot be modeled. In practice, these integrals are approximated by an average over samples of z t ∼ N (0, γ) and W ∼ q.

s t+1 ← s t + ∆ t + k t+1 12: C ← C + c(s t+1 ) 13: return C/K 14: Fit W π by optimizing 1 N N n=1 UNFOLD(s n )

In the second part of our model-based policy search algorithm, we optimize the parameters W π of a policy that minimizes the sum of expected cost over a finite horizon T with respect to our belief q(W). This expected cost is obtained by averaging over multiple virtual roll-outs. For each roll-out we sample W i ∼ q and then simulate state trajectories using the model s t+1 = f (s t , a t , z t ; W i ) + t+1 with policy a t = π(s t ; W π ), input noise z t ∼ N (0, γ) and additive noise t+1 ∼ N (0, Σ). This procedure allows us to obtain estimates of the policy's expected cost for any particular cost function.

If model, policy and cost function are differentiable, we are then able to tune W π by stochastic gradient descent over the roll-out average.

Given a cost function c(s t ), the objective to be optimized by our policy search algorithm is

EQUATION

We approximate (15) by using ( 14), replacing a t with π(s t ; W π ) and using sampling to approximate the expectations:

EQUATION

)

The first line in ( 16) is obtained by using the assumption that the dynamics are Markovian with respect to the current state and the current action and by replacing p(s t |s t-1 , a t-1 ) with the right-hand side of ( 14). In the second line, s W,{z1,...,zt},{ 1,..., t},Wπ t

is the state that is obtained at time t in a roll-out generated by using a policy with parameters W π , a transition function parameterized by W and input noise z 1 , . . . , z t , with additive noise values 1 , . . . , t . In the last line we have approximated the integration with respect to W, z 1 , . . . , z T , 1 , . . . , T and s 0 by averaging over K samples of these variables. To sample s 0 , we draw this variable uniformly from the available transitions (s t , a t , s t+1 ).

The expected cost (15) can then be optimized by stochastic gradient descent using the gradients of the Monte Carlo approximation given by the last line of ( 16). Algorithm 1 computes this Monte Table 1 : Policy performances over different benchmarks. Printed are average values over 5 runs with respective standard errors. Bottom row is the average rank over all 5 × 3 runs.

Carlo approximation. The gradients can then be obtained using automatic differentiation tools such as Theano (Theano Development Team). Note that Algorithm 1 uses the BNNs to make predictions for the change in the state ∆ t = s t+1 -s t instead of for the next state s t+1 since this approach often performs better in practice (Deisenroth & Rasmussen, 2011) .

4. EXPERIMENTS:

We now evaluate the performance of our algorithm for policy search in different benchmark problems. These problems are chosen based on two reasons. First, they contain complex stochastic dynamics and second, they represent real-world applications common in industrial settings. A theano implementation of algorithm 1 is available online1 . See the appendix B for a short introduction to all methods we compare to and appendix C for the hyper-parameters used.

4.1. WET-CHICKEN BENCHMARK:

The Wet-Chicken benchmark (Tresp, 1994) is a challenging problem for model-based policy search that presents both bi-modal and heteroskedastic transition dynamics. We use the two-dimensional version of the problem (Hans & Udluft, 2009) and extend it to the continuous case.

In this problem, a canoeist is paddling on a two-dimensional river. The canoeist's position at time t is (x t , y t ). The river has width w = 5 and length l = 5 with a waterfall at the end, that is, at y t = l. The canoeist wants to move as close to the waterfall as possible because at time t he gets reward r t = -(l -y t ). However, going beyond the waterfall boundary makes the canoeist fall down, having to start back again at the origin (0, 0). At time t the canoeist can choose an action (a t,x , a t,y ) ∈ [-1, 1] 2 that represents the direction and magnitude of his paddling. The river dynamics have stochastic turbulences s t and drift v t that depend on the canoeist's position on the x axis. The larger x t , the larger the drift and the smaller x t , the larger the turbulences. The underlying dynamics are given by the following system of equations. The drift and the turbulence magnitude are given by v t = 3x t w -1 and s t = 3.5 -v t , respectively. The new location (x t+1 , y t+1 ) is given by the current 

EQUATION

)

where ŷt+1 = y t + (a t,y -1) + v t + s t τ t and τ t ∼ Unif([-1, 1]) is a random variable that represents the current turbulence. These dynamics result in rich transition distributions depending on the position as illustrated by the plots in Figure 2 . As the canoeist moves closer to the waterfall, the distribution for the next state becomes increasingly bi-modal (see Figure 1c ) because when he is close to the waterfall, the change in the current location can be large if the canoeist falls down the waterfall and starts again at (0, 0). The distribution may also be truncated uniform for states close to the borders (see Figure 1d ). Furthermore the system has heteroskedastic noise, the smaller the value of x t the higher the noise variance (compare Figure 1a with 1b). Because of these properties, the Wet-Chicken problem is especially difficult for model-based reinforcement learning methods. To our knowledge it has only been solved using model-free approaches after a discretization of the state and action sets (Hans & Udluft, 2009) . For model training we use a batch 2500 random state transitions.

The predictive distributions of different models for y t+1 are shown in Figure 2 for specific choices of (x t , y t ) and (a x,t , a y,t ). These plots show that BNNs with α = 0.5 are very close to the ground-truth.

While it is expected that Gaussian processes fail to model multi-modalities in Figure 1c , the FTIC approximation allows them to model the heteroskedasticity to an extent. VB captures the stochastic patterns on a global level, but often under or over-estimates the true probability density in specific regions. The test-loglikelihood and test MSE in y-dimension are reported in Table 2 for all methods.

(the transitions for x are deterministic given y).

After fitting the models, we train policies using Algorithm 1 with a horizon of size T = 5. Table 1 shows the average reward obtained by each method. BNNs with α = 0.5 perform best and produce policies that are very close to the optimal upper bound, as indicated by the performance of the particle swarm optimization policy (PSO-P). In this problem VB seems to lack robustness and has much larger empirical variance across experiment repetitions than α = 0.5 or α = 1.0.

Figure 3 shows three example policies, π VB ,π α=0.5 and π GP (Figure 3a ,3b and 3c, respectively). The policies obtained by BNNs with random inputs (VB and α = 0.5) show a richer selection of actions. The biggest differences are in the middle-right regions of the plots, where the drift towards the waterfall is large and the bi-modal transition for y (missed by the GP) is more important.

4.2. INDUSTRIAL APPLICATIONS:

We now present results on two industrial cases. First, we focus on data generated by a real gas turbine and second, we consider a recently introduced simulator called the "industrial benchmark", with code publicly available2 (Hein et al., 2016b) . According to the authors: "The "industrial benchmark" aims at being realistic in the sense, that it includes a variety of aspects that we found to be vital in industrial applications." c(t) that together form the reward signal R(t) = -(3f (t) + c(t)). Also visible is the setpoint S, a constant hyper-parameter of the benchmark that indicates the complexity of the dynamics.

For each setpoint S ∈ {10, 20, • • • , 100} we generate 7 trajectories of length 1000 using random exploration. This batch with 70, 000 state transitions forms the training set. We use 30, 000 state transitions, consisting of 3 trajectories for each setpoint, as test set.

For data preprocessing, in addition to the standard normalization process, we apply a log transformation to the reward variable. Because the reward is bounded in the interval [0, R max ], we use a logit transformation to map this interval into the real line. We define the functional form for the dynamics as

R t = f (A t-15 , • • • , A t , R t-15 , • • • , R t-1 ).

The test errors and log-likelihood are given in Table 2 . We see that BNNs with α = 0.5 and α = 1.0 perform best here, whereas Gaussian processes or the MLP obtain rather poor results.

Each row in Figure 4 visualizes long term predictions of the MLP and BNNs trained with VB and α = 0.5 in two specific cases. In the top row we see that while all three methods produce wrong predictions in expectation (compare dark blue curve to red curve). However, BNNs trained with V B and with α = 0.5 exhibit a bi-modal distribution of predicted trajectories, with one mode following the ground-truth very closely. By contrast, the MLP misses the upper mode completely. The bottom row shows that the VB and α = 0.5 also produce more tight confident bands in other settings.

Next, we learn policies using the trained models. Here we use a relatively long horizon of T = 75 steps. Table 1 shows average rewards obtained when applying the policies to the real dynamics.

Because both benchmark and models have an autoregressive component, we do an initial warm-up phase using random exploration before we apply the policies to the system and start to measure rewards.

We observe that GPs perform very poorly in this benchmark. We believe the reason for this is the long search horizon, which makes the uncertainties in the predictive distributions of the GPs become very large. Tighter confidence bands, as illustrated in Figure 4 seem to be key for learning good policies. Overall, α = 1.0 performs best with α = 0.5 being very close.

5. RELATED WORK:

There has been relatively little attention to using Bayesian neural networks for reinforcement learning.

In Blundell et al. (2015) a Thompson sampling approach is used for a contextual bandits problem; the focus is tackling the exploration-exploitation trade-off, while the work in Watter et al. (2015) combines variational auto-encoder with stochastic optimal control for visual data. Compared to our approach the first of these contributions focusses on the exploration/exploitation dilemma, while the second one uses a stochastic optimal control approach to solve the learning problem. By contrast, our work seeks to find an optimal parameterized policy.

Policy gradient techniques are a prominent class of policy search algorithms (Peters & Schaal, 2008) . While model-based approaches were often used in discrete spaces (Wang & Dietterich, 2003) , model-free approaches tended to be more popular in continuous spaces (e.g. Peters & Schaal (2006) ).

Our work can be seen as a Monte-Carlo model-based policy gradient technique in continuous stochastic systems. Similar work was done using Gaussian processes (Deisenroth & Rasmussen, 2011) and with recurrent neural networks (Schaefer et al., 2007) . The Gaussian process approach, while restricted to a Gaussian state distribution, allows propagating beliefs over the roll-out procedure.

More recently Gu et al. (2016) augment a model-free learning procedure with data generated from model-based roll-outs.

6. CONCLUSION AND FUTURE WORK:

We have extended the standard Bayesian neural network (BNN) model with the addition of a random input noise source z. This enables principled Bayesian inference over complex stochastic functions.

We have shown that our BNNs with random inputs can be trained with high accuracy by minimizing α-divergences, with α = 0.5, which often produces better results than variational Bayes. We have also presented an algorithm that uses random roll-outs and stochastic optimization for learning a parameterized policy in a batch scenario. This algorithm particular suited for industry domains.

Our BNNs with random inputs have allowed us to solve a challenging benchmark problem where model-based approaches usually fail. They have also shown promising results on industry benchmarks including real-world data from a gas turbine. In particular, our experiments indicate that a BNN trained with α = 0.5 as divergence measure in conjunction with the presented algorithm for policy optimization is a powerful black-box tool for policy search.

As future work we will consider safety and exploration. For safety, we believe having uncertainty over the underlaying stochastic functions will allows us to optimize policies by focusing on worst case results instead of on average performance. For exploration, having uncertainty on the stochastic functions will be useful for efficient data collection. A ROBUSTNESS OF α = 0.5 AND α = 1.0 WHEN q(z) IS NOT LEARNED

We evaluate the accuracy of the predictive distributions generated by BNNs with stochastic inputs trained by minimizing (8) for different BNNs parameterized by α in two simple regression problems. The first one is characterized by a bimodal predictive distribution. The second is characterized by a heteroskedastic predictive distribution. In the latter case the magnitude of the noise in the targets changes as a function of the input features.

In the first problem x ∈ [-2, 2] and y is obtained as y = 10 sin(x) + with probability 0.5 and y = 10 cos(x) + , otherwise, where ∼ N (0, 1) and is independent of x. The plot in the top of the 1st column in Figure 5 shows a training dataset obtained by sampling 2500 values of x uniformly at random. The plot clearly shows that the distribution of y for a particular x is bimodal. In the second problem x ∈ [-4, 4] and y is obtained as y = 7 sin(x) + 3| cos(x/2)| . The plot in the bottom of the 1st column in Figure 5 shows a training dataset obtained with 1000 values of x uniformly at random. The plot clearly shows that the distribution of y is heteroskedastic, with a noise variance that is a function of x.

We evaluated the predictive performance obtained by minimizing (8) using α = 0.5 and α = 1.0 and also by running VB. However, we do not learn q(z) and keeping it instead fixed to the prior p(z).

We fitted a neural network with 2 hidden layers and 50 hidden units per layer using Adam with its default parameter values, with a learning rate of 0.01 in the first problem and 0.002 in the second problem. We used mini-batches of size 250 and 1000 training epochs. To approximate the expectations in 8, we draw K = 50 samples from q.

The plots in the 3rd and 4th columns of Figure 5 show the predictions obtained with α = 0.5 and α = 1.0, respectively. In these cases, the predictive distribution is able to capture the bimodality in the first problem and the heteroskedasticity pattern in the second problem in both cases The plots in the 2nd column of Figure 5 show the predictions obtained with VB, which converges to suboptimal solutions in which the predictive distribution has a single mode (in the first problem) or is homoskedastic (in the second problem). Tables 3 and 4 show the average test RMSE and log-likelihood obtained by each method on each problem.

These results show that Bayesian neural networks trained with α = 0.5 or α = 1.0 are more robust than VB and can still model complex predictive distributions, which may be multimodal and heteroskedastic, even when q(z) is not learned and is instead kept fixed to the prior p(z). By contrast, VB fails to capture complex stochastic patterns in this setting.

:

https://github.com/siemens/policy_search_bb-alpha

http://github.com/siemens/industrialbenchmark

Figures and Tables:

Figure 2: Predictive distribution of y t given by different methods in four different scenarios. Ground truth (red) is obtained by sampling from the real dynamics.

Figure 3: Visualization of three policies in state space. Waterfall is indicated by top black bar. Left: policy π V B obtained with a BNN trained with VB. Avg. reward is -2.53. Middle: policy π α=0.5 obtained with a BNN trained with α = 0.5. Avg. reward is -2.31. Right: policy π GP obtained by using a Gaussian process model. Avg. reward is -2.94. Color and arrow indicate direction of paddling of policy when in state s t , arrow length indicates action magnitude. Best viewed in color. Dataset MLP VB α=0.5 α=1.0 GP PSO-P Wetchicken -2.71±0.09 -2.67±0.10 -2.37±0.01 -2.42±0.01 -3.05±0.06 -2.34 Turbine -0.65±0.14 -0.45±0.02 -0.41±0.03 -0.55±0.08 -0.64±0.18 NA Industrial -183.5±4.1 -180.2±0.6 -174.2±1.1 -171.1±2.1 -285.2±20.5 -145.5 Avg. Rank 3.6±0.3 3.1±0.2 1.5±0.2 2.3±0.3 4.5±0.3

Figure 5: Ground truth and predictive distributions for two toy problems introduced in main text. Top: bi-modal prediction problem, Bottom: heteroskedastic prediction problem. Left column: Training data (blue points) and ground truth functions (red). Columns 2-4: predictions generated with VB, α = 0.5 and α = 1.0, respectively.

Model test error and test log-likelihood for different benchmarks. Printed are average values over 5 runs with respective standard errors. Bottom row is the average rank over all 5 × 3 runs. location (x t , y t ) and current action (a t,x , a t,y ) using

Table content: <table><tr><td>Dataset</td><td>MLP</td><td>VB</td><td>α=0.5</td><td>α=1.0</td><td>GP</td></tr><tr><td>MSE</td><td/><td/><td/><td/><td/></tr><tr><td>WetChicken</td><td>1.289±0.013</td><td>1.347±0.015</td><td>1.347±0.008</td><td>1.359±0.017</td><td>1.359±0.017</td></tr><tr><td>Turbine</td><td>0.16±0.001</td><td>0.21±0.003</td><td>0.192±0.002</td><td>0.237±0.004</td><td>0.492±0.026</td></tr><tr><td>Industrial</td><td colspan="5">0.0186±0.0052 0.0182±0.0052 0.017±0.0046 0.0171±0.0047 0.0233±0.0049</td></tr><tr><td>Avg. Rank</td><td>2.0±0.34</td><td>3.1±0.24</td><td>2.4±0.23</td><td>2.9±0.36</td><td>4.6±0.23</td></tr><tr><td>Log-Likelihood</td><td/><td/><td/><td/><td/></tr><tr><td>WetChicken</td><td colspan="5">-1.755±0.003 -1.140±0.033 -1.057±0.014 -1.070±0.011 -1.722±0.011</td></tr><tr><td>Turbine</td><td colspan="5">-0.868±0.007 -0.775±0.004 -0.746±0.013 -0.774±0.015 -2.663±0.131</td></tr><tr><td>Industrial</td><td>0.767±0.047</td><td>1.132±0.064</td><td>1.328±0.108</td><td>1.326±0.098</td><td>0.724±0.04</td></tr><tr><td>Avg. Rank</td><td>4.3±0.12</td><td>2.6±0.16</td><td>1.3±0.15</td><td>2.1±0.18</td><td>4.7±0.12</td></tr></table>

Peer review (iclr2017):
This paper introduces an approach for model-based control of stochastic dynamical systems with policy search, based on (1) learning the stochastic dynamics of the underlying system with a Bayesian deep neural network (BNN) that allows some of its inputs to be stochastic, and (2) a policy optimization method based on simulated rollouts from the learned dynamics. BNN training is carried out using alpha-divergence minimization, the specific form of which was introduced in previous work by the authors. Validation and comparison of the approach is undertaken on a simulated domain, as well as real-world scenarios. The paper is tightly written, and easy to follow. Its approach to fitting Bayesian neural networks with alpha divergence is interesting and appears novel in this context. The resulting application to model-based control appears to have significant practical impact, particularly in light of the explainability that a system model can bring to specific decisions made by the policy. As such, I think that the paper brings a valuable contribution to the literature. That said, I have a few questions and suggestions: 1) In section 2.2, it should be explained how the random z_n input is used by the neural network: is it just concatenated to the other inputs and used as-is, or is there a special treatment? 2) Moreover, much case is made for the need to have stochastic inputs, but only a scalar input seems to be provided throughout. Is this enough? How computationally difficult would providing stochastic inputs of higher dimensionality be? 3) How important is the normality assumption in z_n? How is the variance gamma established? 4) It is mentioned that the hidden layers of the neural network are made of rectifiers, but no further utilization of this fact is made in the paper. Is this assumption somehow important in the optimization of the alpha-divergence (beyond what we know about rectifiers to mitigate the vanishing gradient problem) ? 5) Equation (3), denominator mathbf{y} should be mathbf{Y} ? 6) Section 2.3: it would be helpful to have an overview or discussion of the computational complexity of training BNNs, to understand whether and when they can practicably be used. 7) Between eq (12) and (13), a citation to the statement of the time embedding theorem would be helpful, as well as an indication of how the embedding dimension should be chosen. 8) Figure 1: the subplots should have the letters by which they are referenced in the text on p. 7. 9) In section 4.2.1, it is not clear if the gas turbine data is publicly available, and if so where. In addition more details should be provided, such as the dimensionality of the variables E_t, N_t and A_t. 10) Perhaps the comparisons with Gaussian processes should include variants that support stochastic inputs, such as Girard et al. (2003), to provide some of the same modelling capabilities as what’s made use of here. At least, this strand of work should be mentioned in Section 5. References: Girard, A., Rasmussen, C. E., Quiñonero Candela, J., & Murray Smith, R. (2003). Gaussian process priors with uncertain inputs-application to multiple-step ahead time series forecasting. Advances in Neural Information Processing Systems, 545-552.