text:
\subsubsection*{\bibname}}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\usepackage{multirow}
\usepackage{graphicx}
\usepackage{epstopdf}

\usepackage{dblfloatfix}

\input{sections/definitions}

\begin{document}

\twocolumn[

\aistatstitle{Sample-Efficient Reinforcement Learning using Deep Gaussian Processes}

\aistatsauthor{ Charles Gadd\textsuperscript{1} \And Markus Heinonen\textsuperscript{1} \And Harri L\"{a}hdesm\"{a}ki\textsuperscript{1} \And Samuel Kaski\textsuperscript{1,2}}

\aistatsaddress{ \textsuperscript{1}Aalto University, Finland \And \textsuperscript{2}University of Manchester, UK}  ]

\begin{abstract}
\input{sections/S0_abstract}
\end{abstract}

\section{INTRODUCTION}
\input{sections/S1_introduction}

\section{MODEL-BASED REINFORCEMENT LEARNING}
\input{sections/S2_MBRL}

\section{PLANNING WITH MODEL PREDICTIVE CONTROL}
\input{sections/S3_planning}

\section{EXPERIMENTAL RESULTS}
\input{sections/S4_examples}

\section{DISCUSSION}
\input{sections/S5_discussion}

\newpage
\bibliographystyle{plainnat}

\subsubsection*{\bibname}}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\usepackage[dvips]{graphicx}
\input{sections/definitions}
\usepackage{multirow}

\usepackage{mathtools}
\DeclarePairedDelimiter\floor{\lfloor}{\rfloor}
\usepackage{booktabs}
\usepackage{pifont}
\newcommand*\rot{\rotatebox{90}}

\begin{document}

\onecolumn
\aistatstitle{Supplementary Material: Sample-Efficient Reinforcement Learning using Deep Gaussian Processes}

\section{Trajectory sampling - posterior predictive updates}
\input{sections/sup1_TS}

\section{Cross Entropy Method}
\input{sections/sup1_CEM}

\newpage
\section{Reproducing environments and experiments}
\input{sections/sup2_setting}

\newpage
\section{Additional experiments}
\input{sections/sup2_additional_experiments}

\newpage
\subsection{Prior information. Dynamics of contact points.}
\input{sections/sup2_kernel_choice_dynamics}

\bibliographystyle{plainnat}

\subsection{Modeling dynamics with deep Gaussian processes}

Given observed transition tuples $\left\{\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right\}_t$ of the dynamical process, we approximate the transition function with a deep Gaussian process (DGP) \citep{damianou2013deep}, a composition of $L$ Gaussian process layers,

\begin{align} \label{eq:dgp}
    \mathbf{s}_{t+1} &= \mathbf{s}_t + \underbrace{\mathbf{f}_L}_{\mathcal{S} \times \mathcal{A} \mapsto \mathcal{S}} \circ \underbrace{\mathbf{f}_{L-1}  \cdots \circ \mathbf{f}_1(\mathbf{s}_t,\mathbf{a}_t)}_{\mathcal{S} \times \mathcal{A} \mapsto \mathcal{S} \times \mathcal{A}} + \epsilon,
\end{align}
where
\begin{align*}
    \mathbf{f}_\ell &\sim \mathcal{GP}(\mathbf{m}\left(\cdot\right), k_{\pmb{\theta}_\ell}\left(\cdot, \cdot\right)
    ),
\end{align*}
with $\mathcal{GP}$ denoting a Gaussian process (GP) prior over function space, fully determined by a mean function $\mathbf{m}$ and kernel function $k$ with hyper-parameters $\pmb{\theta}_\ell$ \citep{williams2006gaussian}. A GP defines a distribution over functions $\mathbf{f}_\ell(\cdot)$ such that for any subset of inputs $\{ (\mathbf{s},\mathbf{a}) \}$ the function outputs $\{ \mathbf{f}_\ell(\mathbf{s},\mathbf{a}) \}$ follow a Gaussian distribution with a mean and covariance
\begin{align}
    \mathbb{E}[\mathbf{f}_\ell(\mathbf{s},\mathbf{a})] &= \mathbf{m}(\mathbf{s},\mathbf{a}) \\
    \mathbf{cov}[\mathbf{f}_\ell(\mathbf{s},\mathbf{a}), \mathbf{f}_\ell(\mathbf{s}',\mathbf{a}')] &= k_{\pmb{\theta}_\ell}((\mathbf{s},\mathbf{a}), (\mathbf{s}',\mathbf{a}')).
\end{align}
In addition, $\epsilon$ is Gaussian additive noise with zero mean and a precision hyper-parameter $\beta$. This is included in our statistical model as a nugget on the last layer.

In Equation \eqref{eq:dgp} the first $L-1$ layers $\mathbf{f}_\ell : \mathcal{S} \times \mathcal{A} \mapsto \mathcal{S} \times \mathcal{A}$ \emph{input-warp} the state-action pair $(\mathbf{s}_t,\mathbf{a}_t)$, and the final layer Gaussian process $\mathbf{f}_L : \mathcal{S} \times \mathcal{A} \mapsto \mathcal{S}$ maps into state transitions. The input-warping adds to our model its capability to estimate non-Gaussian uncertainty \citep{salimbeni2017doubly} over non-stationary transitions \citep{shen2019learning}.

\subsection{Sparse deep Gaussian process inference}

Following the sparse inducing point constructions of \citep{quinonero2005unifying,salimbeni2017doubly} for scalability, the Gaussian process transitions are guided by $M$ added inducing \emph{input} variables $\mathbf{Z}_\ell = \{ \mathbf{z}_{\ell,m} \}_{m=1}^M$ per layer consisting of state-action pair support vectors $\mathbf{z} \in \mathbb{R}^{S+A}$. Corresponding inducing \emph{outputs} $\mathbf{U}_{\ell} = \{ \mathbf{f}_\ell(\mathbf{z}_{\ell,m}) \}_{m=1}^M$ are also introduced that are in the state-action space for warping layers and in only the state space for final layer $L$.

The joint distribution of the inducing DGP states and random variables at time $t$ follows
\begin{align} \label{eq:dgpjoint}
    p( \mathbf{s}_{t+1}, &\mathbf{f}, \mathbf{U} | \mathbf{s}_t, \mathbf{a}_t, \mathbf{Z}, \boldsymbol{\theta} ) \\
    &= p(\mathbf{s}_{t+1} | \mathbf{f}_L, \beta) \prod_{\ell=0}^{L-1} p(\mathbf{f}_{\ell+1}|\mathbf{f}_\ell, \mathbf{U}_\ell, \mathbf{Z}_\ell, \boldsymbol{\theta}_\ell) p(\mathbf{U}_\ell),
\end{align}

with $\mathbf{f}_0 \triangleq (\mathbf{s}_t,\mathbf{a}_t)$,
\begin{align*}
     p(\mathbf{f}_{\ell+1}|\mathbf{f}_\ell, &\mathbf{U}_\ell, \mathbf{Z}_\ell, \boldsymbol{\theta}_\ell) \\
    &= \mathcal{N}( \mathbf{f}_{\ell+1} | \mathbf{K}_{\mathbf{f}\mathbf{z}} \mathbf{K}_{\mathbf{z}\z}^{-1} \mathrm{vec} \mathbf{U}_\ell, \mathbf{K}_{\mathbf{f}\f} - \mathbf{K}_{\mathbf{f}\mathbf{z}} \mathbf{K}_{\mathbf{z}\z}^{-1} \mathbf{K}_{\mathbf{z}\mathbf{f}})
\end{align*}
and
\begin{equation*}
    p(\mathbf{U}_\ell) = \mathcal{N}(\mathrm{vec} \mathbf{U}_\ell | \mathbf{m}_{\mathbf{z}}, \mathbf{K}_{\mathbf{z}\z}),
\end{equation*}
where we use shorthands $\mathbf{f} = \{ \mathbf{f}_\ell \}$, $\mathbf{U} = \{ \mathbf{U}_\ell \}$, $\mathbf{Z} = \{\mathbf{Z}_\ell \}$, $\boldsymbol{\theta} = \{ \boldsymbol{\theta}_\ell \}$.
The kernel matrices $\mathbf{K}_{\mathbf{f}\f}, \mathbf{K}_{\mathbf{z}\mathbf{f}}, \mathbf{K}_{\mathbf{z}\z}$ are evaluated at all pairs of corresponding sets $\mathbf{f}_\ell$ and $\mathbf{Z}_\ell$ with kernel parameters $\boldsymbol{\theta}_\ell$. We use the ARD technique to estimate separate kernel bandwidths for each input dimension \citep{williams2006gaussian}. The first term of Equation~\eqref{eq:dgpjoint} is the Gaussian likelihood, while the second term is the kernel-based Gaussian process conditional of the next state. The final term is the inducing smoothness prior. Each output dimension is assumed to be independent, conditional on shared hyper-parameters.

We follow the approach of \citep{havasi2018inference} to obtain posterior samples of $(\mathbf{s}_{1:H}, \mathbf{f}, \mathbf{U})$ conditioned on observed trajectories with stochastic gradient Hamiltonian Monte Carlo (SG-HMC) inference, whilst estimating empirical Bayes values for $\boldsymbol{\theta}$ and $\beta$ with a moving window approach. Further trajectories can be sampled efficiently from the posterior estimates of $(\mathbf{f},\mathbf{U})$ to obtain estimates of the trajectory uncertainty under the parameter uncertainty. We add jitter to all kernel inverses for numerical stability, and use separable, diagonal matrix-valued kernels \citep{micchelli2005learning}.

\subsection{Trajectory sampling}

We require a scheme in which we use the model to sample the posterior trajectory given a sequence of proposed actions,  Equation. We desire an approach which accurately captures model uncertainty and is not bound by restrictive distributional assumptions.

\citep{deisenroth2011pilco} used moment matching approximation, which maintains Gaussianity of the state $\mathbf{s}_t$ over time, which prevents estimating multimodal or skewed distributions. In \citep{chua2018deep} transitions are sampled from a bootstrapped model resulting in ensemble reward estimates. We propose an asymptotically exact particle-based trajectory sampling (TS) approach from the Bayesian model. We sample $P$ particles recursively and independently through the DGP-MDP
\begin{align}
    \mathbf{s}_{t+1}^{(p)} \sim \mathbf{s}_t^{(p)} + \tilde{\mathbf{f}}_L \circ \tilde{\mathbf{f}}_{L-1} \circ \cdots \circ \tilde{\mathbf{f}}_1( \mathbf{s}_t^{(p)}, \mathbf{a}_t),
\end{align}
where $p=1,\dots,P$ is the particle index and the $\tilde{\mathbf{f}}_\ell$ indicate the predictive posterior of each layer of our model. The predicted new state, from $\mathbf{s}_t$ under action $\mathbf{a}_t$ is then sampled from our model as $\mathbf{s}_{t+1}^{(p)}$. This procedure is carried out along the planning horizon, where each particle consists of recursive samples, and particles are independent of each other given the initial known state.

These full-uncertainty sample trajectories are then pushed to the deterministic reward function. The long-term reward is the sum of rewards along this horizon, and the expectation under the model is used with the cross-entropy method to optimise the action sequence, as outlined in Algorithm \eqref{alg:ver2}.

\paragraph{Joint predictive posterior samples.} When planning, predictions for each particle are made using the joint predictive posterior distribution across all proposed action sequences simultaneously. The upshot of this is that, within each independent particle, trajectories remain correlated by state and actions. Consequently this smooths our reward objective. If two sequences are similar, or lead to a similar point in the state space then we expect the trajectories to follow consistent dynamics with each other.

\begin{algorithm}[t]
\KwData{Initial data $\mathcal{D}_0 = \left\{\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right\}_t$ obtained by performing random actions. Initial action proposal distribution defined by $\mathbf{m}$ and $\mathbf{v}$.
}
\For{each episode $i$}{
  Train model: optimise $\left\{\mathbf{Z}, \theta, \beta\right\}$ and sample $\mathbf{U}$\\
  \For{each $t = 1, \ldots,$ TaskHorizon}{
    \For{each control optimisation iter}{
    Sample $K=300$ action sequences $\mathbf{a}^{(k)} \sim \mathcal{N}(\mathbf{m}, \mathrm{diag} \, \mathbf{v})$ \\
    For each sequence, sample $P=5$ trajectories under the model\\ \hspace{1mm} $\mathbf{s}^{(p,k)}|\mathbf{a}^{(k)}$  \\
    For each sequence, calculate the expected long term reward over horizon $H$ \\ \hspace{1mm} $\sum_{\tau = t}^{t+H} \frac{1}{P}\sum_{p=1}^P r(\mathbf{s}_{\tau+1}^{\left(p,k\right)},\mathbf{a}_\tau^{(k)})$ \\
    Update $\mathbf{m}$ and $\mathbf{v}$ using the mean and variance of the top performing action sequences \\
    }
    Perform first action of the sequence $\mathbf{a}_t$ \\
    Append observed tuple $\left(\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right)$ to data set.\\
    Re-sample inducing outputs $\pmb{U}$ \\
  }
 }
 \caption{Our DGP-MPC  algorithm.}
  \label{alg:ver2}
  \end{algorithm}

\subsection{Benchmark environments}

We consider three environments upon which to test our model. The first example is a modified cartpole, in which we demonstrate the advantages of improved model capacity even in simple settings. We then consider two additional environments: the reacher, in which the goal is to control a robot arm to reach a target; and the half cheetah environment, in which the goal to teach a cheetah to run. To the best of our knowledge this latter environment has previously posed insurmountable challenges for Gaussian process based models due to the contact points leading to significantly non-stationary dynamics. We demonstrate that increased model capacity allows us to tackle this problem. Additionally, we report significant sample efficiency improvements in the early stages of learning.

\begin{figure}
    \centering
    \subfigure[Cartpole]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/pendulum.png}
        \label{fig:env_cartpole}
    }
    \subfigure[Reacher]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/PR2.png}
        \label{fig:env_reacher}
    }
    \subfigure[Half-cheetah]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/cheetah.png}
        \label{fig:env_cheetah}
    }
\caption{Environments}
\end{figure}

\paragraph{Non-stationary cartpole.}
We consider the simple swinging cartpole task, as depicted in figure~\ref{fig:env_cartpole}, modified such that each episode begins with the cart at the end of the rail. Similar modifications have been considered in previous works with varying motivations,~\citep{kamthe2017data}. Within this task, the objective is to balance the pole above the cart, with a modification that the reward's target lies near the boundary.  Consequently, the true underlying dynamics are then non-stationary near the contact point at the end of the rail. In all cases reported in this example, where relevant we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles. Each episode takes $200$ control iterations.

\paragraph{Reacher.}
We next consider the reacher environment, in which the objective is to direct a robot arm's end effector to reach a target ball, as depicted in figure~\ref{fig:env_reacher}.
For this example we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles.  Each episode takes $150$ control iterations

\paragraph{Half-cheetah.}
Lastly we consider the half-cheetah environment, depicted in figure~\ref{fig:env_cheetah}. Here our task is to make the cheetah run, where more reward is accumulated with greater speed. This environment poses a significant challenge for Gaussian process based models due to the increased number of dimensions and the highly non-stationary dynamics as a result of the contact points between the cheetah's feet and the floor. To the best of our knowledge no Gaussian process based dynamics model has yet been successfully applied to this environment in an RL setting. For this environment we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles. We apply risk-sensitivity to control by penalising large back angles. Each episode takes $1000$ control iterations.

\begin{table*}[!b]
\centering
\begin{tabular}{ |p{3.5cm}||p{3.5cm}||p{2.5cm}|p{2.5cm}|p{2.5cm}|  }
 \hline
\multirow{2}{5em}{Environment}  & \multirow{2}{3em}{Kernel}  & \multicolumn{3}{|c|}{Average episode reward}   \\
& & $L=1$ & $L=2$ & $L=3$ \\
 \hline \hline
\multirow{4}{8em}{Modified Cartpole \hfill over $15$ episodes} &   Squared exponential   & $90.71 \pm 15.43$ & $132.73 \pm 7.10$ & $147.53 \pm 4.01$\\
 &   Mat\'{e}rn-$5/2$   &  $85.38 \pm 15.38$ & $133.36 \pm 6.61$  & $141.34 \pm 3.39$ \\
 &   Mat\'{e}rn-$3/2$ &  $95.97 \pm 15.59$ & $140.81 \pm 4.35$ & $149.81 \pm 3.76$ \\
  &   Mat\'{e}rn-$1/2$ &  $98.33 \pm 9.80$ & $141.98 \pm 4.86$ & $147.67 \pm 4.26$ \\
 \hline
 \multirow{4}{8em}{Reacher \qquad\qquad over $15$ episodes} &   Squared exponential  &  $-50.88 \pm 3.81$ & $-44.62 \pm 2.78$ & $-49.04 \pm 6.38$\\
 &   Mat\'{e}rn-$5/2$    & $-47.42 \pm 3.57$ & $-46.74 \pm 3.97$ & $-44.85 \pm 2.37$ \\
 &   Mat\'{e}rn-$3/2$    & $-45.09 \pm 3.53$ & $-44.39 \pm 1.90$ & $-44.68 \pm 3.14 $ \\
  &   Mat\'{e}rn-$1/2$   &  $-39.89 \pm 0.57 $ & $-40.82 \pm 1.26$ & $-41.28 \pm 2.39$ \\
 \hline
  \multirow{4}{8em}{Half cheetah \qquad over $10$ episodes} &   Squared exponential & $359 \pm 124.5$ & $1256 \pm 133.5$ &  $1498 \pm 70.4$\\
 &   Mat\'{e}rn-$5/2$    & $474 \pm 123.9$ & $1603 \pm 42.7$  & $1612 \pm 112.4$\\
 &   Mat\'{e}rn-$3/2$    & $671 \pm 83.9$ & $1540 \pm 139.5$ & $1680 \pm 83.3$\\
  &   Mat\'{e}rn-$1/2$   & $584 \pm 68.0$ & $1264 \pm 85.9$ & $1400 \pm 116.7$\\
 \hline
\end{tabular}
\label{tab:results}
\caption{Average episode reward for each environment under different kernel choices and model depth. Values reported are the mean and deviation over $10$ independent trials.}
\end{table*}

\subsection{Prior knowledge through functional priors.}

Through choice of kernel and depth, our Bayesian model leverages prior knowledge on the smoothness and non-stationarity of dynamical systems. In this section we study the effect of these model choices on performance.

\paragraph{Modelling contact dynamics.}

We first begin by motivating the use of deep Gaussian processes as an approach which accurately captures uncertainty, whilst incorporating prior information through kernel choice and increasing model capacity through depth. Here we consider the first modified cartpole environment. We train our DGP model on $50$ sequences of length $10$, with each sequence beginning at a random initial state. Sequences of positive actions are  sampled uniformly over $[0,3]$ and the resulting transitions are used to train our DGP model. Contact with the edge of the rail is then observed in all sequences. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=1$ from the initial state mean, and comparing to the ground truth.

We consider the effect of increasing depth and choice of kernel in figure~\ref{fig:kernels_modified_cartpole_dynamics_crash}. Here posterior trajectory samples of Equation~\eqref{eqn:jointtransition} under our DGP model are shown in blue, whilst the true trajectory is shown in red. We clearly observe both improved predictive performance and more accurate uncertainty bounds  when using the once differentiable Mat\'{e}rn-$\frac{3}{2}$ kernel over the infinitely differentiable squared exponential (SEXP). Additionally, increased depth leads to further improvements as this allows the GP model to account for non-stationarity between regions of smooth and contact point dynamics. A more detailed study, including different scenarios, kernel choices, and numbers of layers is given in the supplementary material.

\begin{figure}
    \centering
    \subfigure[SEXP with $L=1$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/SEXP_l1y3.png}
}
\subfigure[SEXP with $L=2$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/SEXP_l2y3.png}
}
\subfigure[Mat\'{e}rn-$\frac{3}{2}$ with $L=1$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/Matern32_l1y3.png}
}
\subfigure[Mat\'{e}rn-$\frac{3}{2}$ with $L=2$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/Matern32_l2y3.png}
}
\caption{Open loop trajectory samples for contact dynamics of the modified Cartpole environment. The change in state of the angular velocity is plotted against the time index of the sequence. Trajectory samples are shown in blue, whilst the true underlying dynamics are shown in red.}
    \label{fig:kernels_modified_cartpole_dynamics_crash}
\end{figure}

\paragraph{Prior knowledge with control.}

We now study the effect of increasing the number of layers in the case of applying control. In figure~\ref{fig:cartpole_layers} we consider the modified cartpole environment, plotting the reward observed with increasing numbers of layers. Each experiment using our approach is repeated to account for stochastic effects. We observe a clear improvement as the number of layers increases. This is due to the increased ability to model the uncertainty of the non-stationary system whilst retaining accurate Bayesian uncertainty estimates. In addition, as the number of layers decreases, the model becomes increasingly miss-specified leading to larger predictive bias, increased variation across seeds, and consequently a less consistently performing controller. Similar results are seen across each environment.

\begin{figure*}
     \centering
     \includegraphics[width=\textwidth, height=4cm]{Images/examples/cartpole/cartpole_layers.png}
     \caption{Modified cartpole reward under our approach with increasing number of layers $L$ and squared exponential kernels. Lines in the background  depict the reward observed across each independent experiment, whilst those in the foreground depict the mean observed reward across experiments.}
     \label{fig:cartpole_layers}
 \end{figure*}

Table~\ref{tab:results} reports the average obtained reward and deviation over independent trials. We choose this metric as it values both final performance and learning rate. Here, we consider the effect of incorporating prior knowledge on both depth and kernel smoothness. We observe that increasing depth greatly benefits those environments with non-stationary dynamics, whilst choosing an appropriately smooth kernel improves performance and improves stability.

We note that whilst not used here, task-specific kernels have been shown to notably improve model specification even further~\citep{calandra2016manifold}. The ability to leverage this prior information within the kernels is a distinct advantage to our approach.

\subsection{Benchmark comparison}

We compare our approach to a number of competing methods. Firstly we consider \textit{PETS}, where the dynamics are modelled with ensembles of probabilistic neural networks and control is performed with planning via MPC. Two different trajectory sampling approaches are used: \textit{TS1}, in which a particle's bootstrap is continuously re-sampled; and \textit{TSinf}, in which a particle's boot strap does not change during a trial~\citep{chua2018deep}. To ensure a fair comparison throughout, all PETS experiments use comparable control parameters to our approach.  We also consider further Gaussian process approaches: \textit{GP-MM}, which approximately propagates the trajectory using moment matching; and \textit{GP-DS}, which approximately propagates the trajectory using distributional sampling. For these later GP experiments we again kept parameters comparable with our own choices.

In figure~\ref{fig:rewards} the maximum seen rewards are plotted against the number of control steps. We report the result for our approach with $L=3$ and for each benchmark we report the mean over $10$ independent runs.
When comparing our approach to each benchmark we clearly observe both a faster learning rate, and in many cases an improved final reward. In the half cheetah task we observe an exceedingly fast initial learning rate, with reward being accumulated from the very beginning given only a single episode of random interactions. In contrast GP-MM and GP-DS were unable to accumulate any reward, whilst PETS requires many more interactions before accumulating positive reward. As expected,  PETS does eventually surpass our performance, but requires many times more samples to do so. Plots showing performance with more episodes are given in the supplementary material, but lie beyond the scope of our sample-efficiency setting.

\begin{figure*}[t]
     \centering
     \includegraphics[width=\textwidth, height=4cm]{Images/examples/results.png}
     \caption{Maximum seen episodic reward, averaged over random seeds. Our approach using a DGP compared against a number of benchmarks.}
\label{fig:rewards}
 \end{figure*}

\begin{figure*}[t]
\centering
\subfigure[Cartpole]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/cartpole/cartpole_var.png}
    \label{fig:cartpole_variance}
}
\subfigure[Reacher]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/reacher/reacher_var.png}
\label{fig:reacher_variance}
}
\subfigure[Half-cheetah]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/cheetah/cheetah_var.png}
    \label{fig:cheetah_variance}
}
\caption{Mean and deviation of the average reward obtained over the first $10$ episodes across independent trials. The larger the average reward, the quicker our agent converged to a strong optima.
}
\label{fig:variances}
\end{figure*}

We also consider the consistency of the learning process. Whilst the maximum seen reward gives a good depiction of the optimal performance of a controller through the exploration process, we are also concerned with the learning curves as depicted in figure~\ref{fig:cartpole_layers}, including their variability.
 For this we consider the variability across seeds of the mean reward across episodes.  The larger this mean reward, the faster a method led to a strong controller; the smaller the variability, the more consistently the model performs. This is shown in figure~\ref{fig:variances}. Variability of this estimator over seeds then also encapsulates erratic differences in rewards obtained between episodes, i.e. due to a controller which performs inconsistently. Our method's mean performance is clearly better than the alternatives in all cases, and variation is at most as large as the others, except for the half-cheetah where the others did not perform at all comparably.

\subsection{Unmodified cartpole}

We include here additional experiments for the unmodified cartpole example, in which the starting point remains at the center of the rail. Figure~\ref{fig:unmodified_cartpole} gives the maximum seen reward plots, in which the squared exponential kernel is used on every layer. We observe that we retain high sample efficiency in the simplified case.

\begin{figure}[ht]
\centering
    \includegraphics[width=\textwidth]{Images/supp/plots/unmodified_cartpole_layer_analysis_maxFalse.png}
\caption{The unmodified cartpole reward with increasing numbers of layers. Each case used the square exponential kernel across all layers.
}
\label{fig:unmodified_cartpole}
\end{figure}

\subsection{Half cheetah - beyond small sample sizes}

We now consider the performance of our method when applied on the half cheetah, beyond the sample-efficiency domain. In Figure~\ref{fig:cheetah_longer} we present the reward obtained using our DGP-MPC approach for the half cheetah environment over a larger number of environment observations.

\begin{figure}[ht]
\centering
\includegraphics[width=\textwidth]{Images/supp/plots/half_cheetah_longer.png}
\caption{The half cheetah reward with increasing numbers of layers. Each case used the Mat\'{e}rn-3/2 kernel across all layers..
}
\label{fig:cheetah_longer}
\end{figure}

\subsubsection*{Modified cartpole} See Figure~\ref{fig:kernels_modified_cartpole}

\begin{figure}[ht]
\centering
\subfigure[Squared exponential kernel on all layers.]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/modified_cartpole_sexp_sexp.png}
    \label{fig:}
}
\subfigure[Squared exponential on first layers, matern $\frac{5}{2}$ kernel on final layer.]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/modified_cartpole_sexp_m52.png}
    \label{fig:}
}
\caption{Modified cartpole environment 'as seen' reward.
}
\label{fig:kernels_modified_cartpole}
\end{figure}

\newpage
\subsubsection*{Reacher}

\newpage
\paragraph{Half cheetah} The case of no reward shaping with different kernel choices is given in~ Figure~\ref{fig:kernels_halfcheetah_noshaping}. The results when adding reward shaping are shown in Figure~\ref{fig:kernels_halfcheetah_shaping}.

\begin{figure}[ht]
\centering
\subfigure[Squared exponential kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_sexp.png}
    \label{fig:}
}
\subfigure[Matern $\frac{5}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern52.png}
    \label{fig:}
}
\subfigure[Matern $\frac{3}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern32.png}
    \label{fig:}
}
\subfigure[Matern $\frac{1}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, with no reward shaping. \textbf{Different kernel} choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_noshaping}
\end{figure}

\begin{figure}[ht]
\centering
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, with reward shaping]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/halfcheetah_shaped_matern12.png}
    \label{fig:}
}
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, using the lower credible bound (minus one deviation) of the pushforward measure on the reward]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/halfcheetah_LCB_m12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, and with \textbf{reward shaping}. Different kernel choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_shaping}
\end{figure}

\begin{figure}[ht]
\centering
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, without reward shaping but with noisy observations]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/half_cheetah_noisy_matern12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, and without shaping and with \textbf{noise added to observations}. Different kernel choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_shaping}
\end{figure}

\subsubsection*{Modified cartpole} We study the effect of deep Gaussian process depth and kernel choice in the modified cartpole environment. Towards this end we consider three situations: a) an isolated non-stationary case of the cart hitting the wall b) an isolated stationary case of the cart travelling along the rail c) a hybrid case containing transitions of both. In each case the trajectory plotted is an open loop, where at each index the next sample is conditional on the previous sampled state, with no feedback. Here, steps are made on a finer time-scale to reduce discretisation errors.

\paragraph{Case a} We train our DGP model on $50$ sequences of length $10$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of positive actions are then uniformly sampled in $[0,3]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=0.75$ from the initial state mean, and comparing to the ground truth.  See Figure~\ref{fig:kernels_modified_cartpole_dynamics_crash}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case a}. Analysis on the effect of kernel choice and depth when modelling the non-stationary contact dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_crash}
\end{figure}

\paragraph{Case b} We train our DGP model on $50$ sequences of length $10$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of negative actions are then uniformly sampled in $[-3,0]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=-1$ from the initial state mean, and comparing to the ground truth. See Figure~\ref{fig:kernels_modified_cartpole_dynamics_no_crash}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case b}. Analysis on the effect of kernel choice and depth when modelling the stationary dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_no_crash}
\end{figure}

\paragraph{Case c} We train our DGP model on $50$ sequences of length $20$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of actions are then uniformly sampled in $[-3,3]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=0.75$ from the initial state mean, and comparing to the ground truth. See Figure~\ref{fig:kernels_modified_cartpole_dynamics_mixed}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case c}. Analysis on the effect of kernel choice and depth when modelling both the stationary and non-stationary dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_mixed}
\end{figure}

\subsubsection*{\bibname}}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\usepackage{multirow}
\usepackage{graphicx}
\usepackage{epstopdf}

\usepackage{dblfloatfix}

\input{sections/definitions}

\begin{document}

\twocolumn[

\aistatstitle{Sample-Efficient Reinforcement Learning using Deep Gaussian Processes}

\aistatsauthor{ Charles Gadd\textsuperscript{1} \And Markus Heinonen\textsuperscript{1} \And Harri L\"{a}hdesm\"{a}ki\textsuperscript{1} \And Samuel Kaski\textsuperscript{1,2}}

\aistatsaddress{ \textsuperscript{1}Aalto University, Finland \And \textsuperscript{2}University of Manchester, UK}  ]

\begin{abstract}
\input{sections/S0_abstract}
\end{abstract}

\section{INTRODUCTION}
\input{sections/S1_introduction}

\section{MODEL-BASED REINFORCEMENT LEARNING}
\input{sections/S2_MBRL}

\section{PLANNING WITH MODEL PREDICTIVE CONTROL}
\input{sections/S3_planning}

\section{EXPERIMENTAL RESULTS}
\input{sections/S4_examples}

\section{DISCUSSION}
\input{sections/S5_discussion}

\newpage
\bibliographystyle{plainnat}

\subsubsection*{\bibname}}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\usepackage[dvips]{graphicx}
\input{sections/definitions}
\usepackage{multirow}

\usepackage{mathtools}
\DeclarePairedDelimiter\floor{\lfloor}{\rfloor}
\usepackage{booktabs}
\usepackage{pifont}
\newcommand*\rot{\rotatebox{90}}

\begin{document}

\onecolumn
\aistatstitle{Supplementary Material: Sample-Efficient Reinforcement Learning using Deep Gaussian Processes}

\section{Trajectory sampling - posterior predictive updates}
\input{sections/sup1_TS}

\section{Cross Entropy Method}
\input{sections/sup1_CEM}

\newpage
\section{Reproducing environments and experiments}
\input{sections/sup2_setting}

\newpage
\section{Additional experiments}
\input{sections/sup2_additional_experiments}

\newpage
\subsection{Prior information. Dynamics of contact points.}
\input{sections/sup2_kernel_choice_dynamics}

\bibliographystyle{plainnat}

\subsection{Modeling dynamics with deep Gaussian processes}

Given observed transition tuples $\left\{\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right\}_t$ of the dynamical process, we approximate the transition function with a deep Gaussian process (DGP) \citep{damianou2013deep}, a composition of $L$ Gaussian process layers,

\begin{align} \label{eq:dgp}
    \mathbf{s}_{t+1} &= \mathbf{s}_t + \underbrace{\mathbf{f}_L}_{\mathcal{S} \times \mathcal{A} \mapsto \mathcal{S}} \circ \underbrace{\mathbf{f}_{L-1}  \cdots \circ \mathbf{f}_1(\mathbf{s}_t,\mathbf{a}_t)}_{\mathcal{S} \times \mathcal{A} \mapsto \mathcal{S} \times \mathcal{A}} + \epsilon,
\end{align}
where
\begin{align*}
    \mathbf{f}_\ell &\sim \mathcal{GP}(\mathbf{m}\left(\cdot\right), k_{\pmb{\theta}_\ell}\left(\cdot, \cdot\right)
    ),
\end{align*}
with $\mathcal{GP}$ denoting a Gaussian process (GP) prior over function space, fully determined by a mean function $\mathbf{m}$ and kernel function $k$ with hyper-parameters $\pmb{\theta}_\ell$ \citep{williams2006gaussian}. A GP defines a distribution over functions $\mathbf{f}_\ell(\cdot)$ such that for any subset of inputs $\{ (\mathbf{s},\mathbf{a}) \}$ the function outputs $\{ \mathbf{f}_\ell(\mathbf{s},\mathbf{a}) \}$ follow a Gaussian distribution with a mean and covariance
\begin{align}
    \mathbb{E}[\mathbf{f}_\ell(\mathbf{s},\mathbf{a})] &= \mathbf{m}(\mathbf{s},\mathbf{a}) \\
    \mathbf{cov}[\mathbf{f}_\ell(\mathbf{s},\mathbf{a}), \mathbf{f}_\ell(\mathbf{s}',\mathbf{a}')] &= k_{\pmb{\theta}_\ell}((\mathbf{s},\mathbf{a}), (\mathbf{s}',\mathbf{a}')).
\end{align}
In addition, $\epsilon$ is Gaussian additive noise with zero mean and a precision hyper-parameter $\beta$. This is included in our statistical model as a nugget on the last layer.

In Equation \eqref{eq:dgp} the first $L-1$ layers $\mathbf{f}_\ell : \mathcal{S} \times \mathcal{A} \mapsto \mathcal{S} \times \mathcal{A}$ \emph{input-warp} the state-action pair $(\mathbf{s}_t,\mathbf{a}_t)$, and the final layer Gaussian process $\mathbf{f}_L : \mathcal{S} \times \mathcal{A} \mapsto \mathcal{S}$ maps into state transitions. The input-warping adds to our model its capability to estimate non-Gaussian uncertainty \citep{salimbeni2017doubly} over non-stationary transitions \citep{shen2019learning}.

\subsection{Sparse deep Gaussian process inference}

Following the sparse inducing point constructions of \citep{quinonero2005unifying,salimbeni2017doubly} for scalability, the Gaussian process transitions are guided by $M$ added inducing \emph{input} variables $\mathbf{Z}_\ell = \{ \mathbf{z}_{\ell,m} \}_{m=1}^M$ per layer consisting of state-action pair support vectors $\mathbf{z} \in \mathbb{R}^{S+A}$. Corresponding inducing \emph{outputs} $\mathbf{U}_{\ell} = \{ \mathbf{f}_\ell(\mathbf{z}_{\ell,m}) \}_{m=1}^M$ are also introduced that are in the state-action space for warping layers and in only the state space for final layer $L$.

The joint distribution of the inducing DGP states and random variables at time $t$ follows
\begin{align} \label{eq:dgpjoint}
    p( \mathbf{s}_{t+1}, &\mathbf{f}, \mathbf{U} | \mathbf{s}_t, \mathbf{a}_t, \mathbf{Z}, \boldsymbol{\theta} ) \\
    &= p(\mathbf{s}_{t+1} | \mathbf{f}_L, \beta) \prod_{\ell=0}^{L-1} p(\mathbf{f}_{\ell+1}|\mathbf{f}_\ell, \mathbf{U}_\ell, \mathbf{Z}_\ell, \boldsymbol{\theta}_\ell) p(\mathbf{U}_\ell),
\end{align}

with $\mathbf{f}_0 \triangleq (\mathbf{s}_t,\mathbf{a}_t)$,
\begin{align*}
     p(\mathbf{f}_{\ell+1}|\mathbf{f}_\ell, &\mathbf{U}_\ell, \mathbf{Z}_\ell, \boldsymbol{\theta}_\ell) \\
    &= \mathcal{N}( \mathbf{f}_{\ell+1} | \mathbf{K}_{\mathbf{f}\mathbf{z}} \mathbf{K}_{\mathbf{z}\z}^{-1} \mathrm{vec} \mathbf{U}_\ell, \mathbf{K}_{\mathbf{f}\f} - \mathbf{K}_{\mathbf{f}\mathbf{z}} \mathbf{K}_{\mathbf{z}\z}^{-1} \mathbf{K}_{\mathbf{z}\mathbf{f}})
\end{align*}
and
\begin{equation*}
    p(\mathbf{U}_\ell) = \mathcal{N}(\mathrm{vec} \mathbf{U}_\ell | \mathbf{m}_{\mathbf{z}}, \mathbf{K}_{\mathbf{z}\z}),
\end{equation*}
where we use shorthands $\mathbf{f} = \{ \mathbf{f}_\ell \}$, $\mathbf{U} = \{ \mathbf{U}_\ell \}$, $\mathbf{Z} = \{\mathbf{Z}_\ell \}$, $\boldsymbol{\theta} = \{ \boldsymbol{\theta}_\ell \}$.
The kernel matrices $\mathbf{K}_{\mathbf{f}\f}, \mathbf{K}_{\mathbf{z}\mathbf{f}}, \mathbf{K}_{\mathbf{z}\z}$ are evaluated at all pairs of corresponding sets $\mathbf{f}_\ell$ and $\mathbf{Z}_\ell$ with kernel parameters $\boldsymbol{\theta}_\ell$. We use the ARD technique to estimate separate kernel bandwidths for each input dimension \citep{williams2006gaussian}. The first term of Equation~\eqref{eq:dgpjoint} is the Gaussian likelihood, while the second term is the kernel-based Gaussian process conditional of the next state. The final term is the inducing smoothness prior. Each output dimension is assumed to be independent, conditional on shared hyper-parameters.

We follow the approach of \citep{havasi2018inference} to obtain posterior samples of $(\mathbf{s}_{1:H}, \mathbf{f}, \mathbf{U})$ conditioned on observed trajectories with stochastic gradient Hamiltonian Monte Carlo (SG-HMC) inference, whilst estimating empirical Bayes values for $\boldsymbol{\theta}$ and $\beta$ with a moving window approach. Further trajectories can be sampled efficiently from the posterior estimates of $(\mathbf{f},\mathbf{U})$ to obtain estimates of the trajectory uncertainty under the parameter uncertainty. We add jitter to all kernel inverses for numerical stability, and use separable, diagonal matrix-valued kernels \citep{micchelli2005learning}.

\subsection{Trajectory sampling}

We require a scheme in which we use the model to sample the posterior trajectory given a sequence of proposed actions,  Equation. We desire an approach which accurately captures model uncertainty and is not bound by restrictive distributional assumptions.

\citep{deisenroth2011pilco} used moment matching approximation, which maintains Gaussianity of the state $\mathbf{s}_t$ over time, which prevents estimating multimodal or skewed distributions. In \citep{chua2018deep} transitions are sampled from a bootstrapped model resulting in ensemble reward estimates. We propose an asymptotically exact particle-based trajectory sampling (TS) approach from the Bayesian model. We sample $P$ particles recursively and independently through the DGP-MDP
\begin{align}
    \mathbf{s}_{t+1}^{(p)} \sim \mathbf{s}_t^{(p)} + \tilde{\mathbf{f}}_L \circ \tilde{\mathbf{f}}_{L-1} \circ \cdots \circ \tilde{\mathbf{f}}_1( \mathbf{s}_t^{(p)}, \mathbf{a}_t),
\end{align}
where $p=1,\dots,P$ is the particle index and the $\tilde{\mathbf{f}}_\ell$ indicate the predictive posterior of each layer of our model. The predicted new state, from $\mathbf{s}_t$ under action $\mathbf{a}_t$ is then sampled from our model as $\mathbf{s}_{t+1}^{(p)}$. This procedure is carried out along the planning horizon, where each particle consists of recursive samples, and particles are independent of each other given the initial known state.

These full-uncertainty sample trajectories are then pushed to the deterministic reward function. The long-term reward is the sum of rewards along this horizon, and the expectation under the model is used with the cross-entropy method to optimise the action sequence, as outlined in Algorithm \eqref{alg:ver2}.

\paragraph{Joint predictive posterior samples.} When planning, predictions for each particle are made using the joint predictive posterior distribution across all proposed action sequences simultaneously. The upshot of this is that, within each independent particle, trajectories remain correlated by state and actions. Consequently this smooths our reward objective. If two sequences are similar, or lead to a similar point in the state space then we expect the trajectories to follow consistent dynamics with each other.

\begin{algorithm}[t]
\KwData{Initial data $\mathcal{D}_0 = \left\{\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right\}_t$ obtained by performing random actions. Initial action proposal distribution defined by $\mathbf{m}$ and $\mathbf{v}$.
}
\For{each episode $i$}{
  Train model: optimise $\left\{\mathbf{Z}, \theta, \beta\right\}$ and sample $\mathbf{U}$\\
  \For{each $t = 1, \ldots,$ TaskHorizon}{
    \For{each control optimisation iter}{
    Sample $K=300$ action sequences $\mathbf{a}^{(k)} \sim \mathcal{N}(\mathbf{m}, \mathrm{diag} \, \mathbf{v})$ \\
    For each sequence, sample $P=5$ trajectories under the model\\ \hspace{1mm} $\mathbf{s}^{(p,k)}|\mathbf{a}^{(k)}$  \\
    For each sequence, calculate the expected long term reward over horizon $H$ \\ \hspace{1mm} $\sum_{\tau = t}^{t+H} \frac{1}{P}\sum_{p=1}^P r(\mathbf{s}_{\tau+1}^{\left(p,k\right)},\mathbf{a}_\tau^{(k)})$ \\
    Update $\mathbf{m}$ and $\mathbf{v}$ using the mean and variance of the top performing action sequences \\
    }
    Perform first action of the sequence $\mathbf{a}_t$ \\
    Append observed tuple $\left(\mathbf{s}_t, \mathbf{a}_t, \mathbf{s}_{t+1}\right)$ to data set.\\
    Re-sample inducing outputs $\pmb{U}$ \\
  }
 }
 \caption{Our DGP-MPC  algorithm.}
  \label{alg:ver2}
  \end{algorithm}

\subsection{Benchmark environments}

We consider three environments upon which to test our model. The first example is a modified cartpole, in which we demonstrate the advantages of improved model capacity even in simple settings. We then consider two additional environments: the reacher, in which the goal is to control a robot arm to reach a target; and the half cheetah environment, in which the goal to teach a cheetah to run. To the best of our knowledge this latter environment has previously posed insurmountable challenges for Gaussian process based models due to the contact points leading to significantly non-stationary dynamics. We demonstrate that increased model capacity allows us to tackle this problem. Additionally, we report significant sample efficiency improvements in the early stages of learning.

\begin{figure}
    \centering
    \subfigure[Cartpole]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/pendulum.png}
        \label{fig:env_cartpole}
    }
    \subfigure[Reacher]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/PR2.png}
        \label{fig:env_reacher}
    }
    \subfigure[Half-cheetah]{
        \includegraphics[width=0.26\columnwidth]{Images/environments/cheetah.png}
        \label{fig:env_cheetah}
    }
\caption{Environments}
\end{figure}

\paragraph{Non-stationary cartpole.}
We consider the simple swinging cartpole task, as depicted in figure~\ref{fig:env_cartpole}, modified such that each episode begins with the cart at the end of the rail. Similar modifications have been considered in previous works with varying motivations,~\citep{kamthe2017data}. Within this task, the objective is to balance the pole above the cart, with a modification that the reward's target lies near the boundary.  Consequently, the true underlying dynamics are then non-stationary near the contact point at the end of the rail. In all cases reported in this example, where relevant we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles. Each episode takes $200$ control iterations.

\paragraph{Reacher.}
We next consider the reacher environment, in which the objective is to direct a robot arm's end effector to reach a target ball, as depicted in figure~\ref{fig:env_reacher}.
For this example we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles.  Each episode takes $150$ control iterations

\paragraph{Half-cheetah.}
Lastly we consider the half-cheetah environment, depicted in figure~\ref{fig:env_cheetah}. Here our task is to make the cheetah run, where more reward is accumulated with greater speed. This environment poses a significant challenge for Gaussian process based models due to the increased number of dimensions and the highly non-stationary dynamics as a result of the contact points between the cheetah's feet and the floor. To the best of our knowledge no Gaussian process based dynamics model has yet been successfully applied to this environment in an RL setting. For this environment we sample $K=300$ action sequences, use $200$ inducing points, and use $P=5$ particles. We apply risk-sensitivity to control by penalising large back angles. Each episode takes $1000$ control iterations.

\begin{table*}[!b]
\centering
\begin{tabular}{ |p{3.5cm}||p{3.5cm}||p{2.5cm}|p{2.5cm}|p{2.5cm}|  }
 \hline
\multirow{2}{5em}{Environment}  & \multirow{2}{3em}{Kernel}  & \multicolumn{3}{|c|}{Average episode reward}   \\
& & $L=1$ & $L=2$ & $L=3$ \\
 \hline \hline
\multirow{4}{8em}{Modified Cartpole \hfill over $15$ episodes} &   Squared exponential   & $90.71 \pm 15.43$ & $132.73 \pm 7.10$ & $147.53 \pm 4.01$\\
 &   Mat\'{e}rn-$5/2$   &  $85.38 \pm 15.38$ & $133.36 \pm 6.61$  & $141.34 \pm 3.39$ \\
 &   Mat\'{e}rn-$3/2$ &  $95.97 \pm 15.59$ & $140.81 \pm 4.35$ & $149.81 \pm 3.76$ \\
  &   Mat\'{e}rn-$1/2$ &  $98.33 \pm 9.80$ & $141.98 \pm 4.86$ & $147.67 \pm 4.26$ \\
 \hline
 \multirow{4}{8em}{Reacher \qquad\qquad over $15$ episodes} &   Squared exponential  &  $-50.88 \pm 3.81$ & $-44.62 \pm 2.78$ & $-49.04 \pm 6.38$\\
 &   Mat\'{e}rn-$5/2$    & $-47.42 \pm 3.57$ & $-46.74 \pm 3.97$ & $-44.85 \pm 2.37$ \\
 &   Mat\'{e}rn-$3/2$    & $-45.09 \pm 3.53$ & $-44.39 \pm 1.90$ & $-44.68 \pm 3.14 $ \\
  &   Mat\'{e}rn-$1/2$   &  $-39.89 \pm 0.57 $ & $-40.82 \pm 1.26$ & $-41.28 \pm 2.39$ \\
 \hline
  \multirow{4}{8em}{Half cheetah \qquad over $10$ episodes} &   Squared exponential & $359 \pm 124.5$ & $1256 \pm 133.5$ &  $1498 \pm 70.4$\\
 &   Mat\'{e}rn-$5/2$    & $474 \pm 123.9$ & $1603 \pm 42.7$  & $1612 \pm 112.4$\\
 &   Mat\'{e}rn-$3/2$    & $671 \pm 83.9$ & $1540 \pm 139.5$ & $1680 \pm 83.3$\\
  &   Mat\'{e}rn-$1/2$   & $584 \pm 68.0$ & $1264 \pm 85.9$ & $1400 \pm 116.7$\\
 \hline
\end{tabular}
\label{tab:results}
\caption{Average episode reward for each environment under different kernel choices and model depth. Values reported are the mean and deviation over $10$ independent trials.}
\end{table*}

\subsection{Prior knowledge through functional priors.}

Through choice of kernel and depth, our Bayesian model leverages prior knowledge on the smoothness and non-stationarity of dynamical systems. In this section we study the effect of these model choices on performance.

\paragraph{Modelling contact dynamics.}

We first begin by motivating the use of deep Gaussian processes as an approach which accurately captures uncertainty, whilst incorporating prior information through kernel choice and increasing model capacity through depth. Here we consider the first modified cartpole environment. We train our DGP model on $50$ sequences of length $10$, with each sequence beginning at a random initial state. Sequences of positive actions are  sampled uniformly over $[0,3]$ and the resulting transitions are used to train our DGP model. Contact with the edge of the rail is then observed in all sequences. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=1$ from the initial state mean, and comparing to the ground truth.

We consider the effect of increasing depth and choice of kernel in figure~\ref{fig:kernels_modified_cartpole_dynamics_crash}. Here posterior trajectory samples of Equation~\eqref{eqn:jointtransition} under our DGP model are shown in blue, whilst the true trajectory is shown in red. We clearly observe both improved predictive performance and more accurate uncertainty bounds  when using the once differentiable Mat\'{e}rn-$\frac{3}{2}$ kernel over the infinitely differentiable squared exponential (SEXP). Additionally, increased depth leads to further improvements as this allows the GP model to account for non-stationarity between regions of smooth and contact point dynamics. A more detailed study, including different scenarios, kernel choices, and numbers of layers is given in the supplementary material.

\begin{figure}
    \centering
    \subfigure[SEXP with $L=1$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/SEXP_l1y3.png}
}
\subfigure[SEXP with $L=2$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/SEXP_l2y3.png}
}
\subfigure[Mat\'{e}rn-$\frac{3}{2}$ with $L=1$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/Matern32_l1y3.png}
}
\subfigure[Mat\'{e}rn-$\frac{3}{2}$ with $L=2$.]{
    \includegraphics[width=0.45\columnwidth]{Images/supp/plots/dynamics/crash/Matern32_l2y3.png}
}
\caption{Open loop trajectory samples for contact dynamics of the modified Cartpole environment. The change in state of the angular velocity is plotted against the time index of the sequence. Trajectory samples are shown in blue, whilst the true underlying dynamics are shown in red.}
    \label{fig:kernels_modified_cartpole_dynamics_crash}
\end{figure}

\paragraph{Prior knowledge with control.}

We now study the effect of increasing the number of layers in the case of applying control. In figure~\ref{fig:cartpole_layers} we consider the modified cartpole environment, plotting the reward observed with increasing numbers of layers. Each experiment using our approach is repeated to account for stochastic effects. We observe a clear improvement as the number of layers increases. This is due to the increased ability to model the uncertainty of the non-stationary system whilst retaining accurate Bayesian uncertainty estimates. In addition, as the number of layers decreases, the model becomes increasingly miss-specified leading to larger predictive bias, increased variation across seeds, and consequently a less consistently performing controller. Similar results are seen across each environment.

\begin{figure*}
     \centering
     \includegraphics[width=\textwidth, height=4cm]{Images/examples/cartpole/cartpole_layers.png}
     \caption{Modified cartpole reward under our approach with increasing number of layers $L$ and squared exponential kernels. Lines in the background  depict the reward observed across each independent experiment, whilst those in the foreground depict the mean observed reward across experiments.}
     \label{fig:cartpole_layers}
 \end{figure*}

Table~\ref{tab:results} reports the average obtained reward and deviation over independent trials. We choose this metric as it values both final performance and learning rate. Here, we consider the effect of incorporating prior knowledge on both depth and kernel smoothness. We observe that increasing depth greatly benefits those environments with non-stationary dynamics, whilst choosing an appropriately smooth kernel improves performance and improves stability.

We note that whilst not used here, task-specific kernels have been shown to notably improve model specification even further~\citep{calandra2016manifold}. The ability to leverage this prior information within the kernels is a distinct advantage to our approach.

\subsection{Benchmark comparison}

We compare our approach to a number of competing methods. Firstly we consider \textit{PETS}, where the dynamics are modelled with ensembles of probabilistic neural networks and control is performed with planning via MPC. Two different trajectory sampling approaches are used: \textit{TS1}, in which a particle's bootstrap is continuously re-sampled; and \textit{TSinf}, in which a particle's boot strap does not change during a trial~\citep{chua2018deep}. To ensure a fair comparison throughout, all PETS experiments use comparable control parameters to our approach.  We also consider further Gaussian process approaches: \textit{GP-MM}, which approximately propagates the trajectory using moment matching; and \textit{GP-DS}, which approximately propagates the trajectory using distributional sampling. For these later GP experiments we again kept parameters comparable with our own choices.

In figure~\ref{fig:rewards} the maximum seen rewards are plotted against the number of control steps. We report the result for our approach with $L=3$ and for each benchmark we report the mean over $10$ independent runs.
When comparing our approach to each benchmark we clearly observe both a faster learning rate, and in many cases an improved final reward. In the half cheetah task we observe an exceedingly fast initial learning rate, with reward being accumulated from the very beginning given only a single episode of random interactions. In contrast GP-MM and GP-DS were unable to accumulate any reward, whilst PETS requires many more interactions before accumulating positive reward. As expected,  PETS does eventually surpass our performance, but requires many times more samples to do so. Plots showing performance with more episodes are given in the supplementary material, but lie beyond the scope of our sample-efficiency setting.

\begin{figure*}[t]
     \centering
     \includegraphics[width=\textwidth, height=4cm]{Images/examples/results.png}
     \caption{Maximum seen episodic reward, averaged over random seeds. Our approach using a DGP compared against a number of benchmarks.}
\label{fig:rewards}
 \end{figure*}

\begin{figure*}[t]
\centering
\subfigure[Cartpole]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/cartpole/cartpole_var.png}
    \label{fig:cartpole_variance}
}
\subfigure[Reacher]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/reacher/reacher_var.png}
\label{fig:reacher_variance}
}
\subfigure[Half-cheetah]{
    \includegraphics[width=0.31\textwidth, height=4.25cm]{Images/examples/cheetah/cheetah_var.png}
    \label{fig:cheetah_variance}
}
\caption{Mean and deviation of the average reward obtained over the first $10$ episodes across independent trials. The larger the average reward, the quicker our agent converged to a strong optima.
}
\label{fig:variances}
\end{figure*}

We also consider the consistency of the learning process. Whilst the maximum seen reward gives a good depiction of the optimal performance of a controller through the exploration process, we are also concerned with the learning curves as depicted in figure~\ref{fig:cartpole_layers}, including their variability.
 For this we consider the variability across seeds of the mean reward across episodes.  The larger this mean reward, the faster a method led to a strong controller; the smaller the variability, the more consistently the model performs. This is shown in figure~\ref{fig:variances}. Variability of this estimator over seeds then also encapsulates erratic differences in rewards obtained between episodes, i.e. due to a controller which performs inconsistently. Our method's mean performance is clearly better than the alternatives in all cases, and variation is at most as large as the others, except for the half-cheetah where the others did not perform at all comparably.

\subsection{Unmodified cartpole}

We include here additional experiments for the unmodified cartpole example, in which the starting point remains at the center of the rail. Figure~\ref{fig:unmodified_cartpole} gives the maximum seen reward plots, in which the squared exponential kernel is used on every layer. We observe that we retain high sample efficiency in the simplified case.

\begin{figure}[ht]
\centering
    \includegraphics[width=\textwidth]{Images/supp/plots/unmodified_cartpole_layer_analysis_maxFalse.png}
\caption{The unmodified cartpole reward with increasing numbers of layers. Each case used the square exponential kernel across all layers.
}
\label{fig:unmodified_cartpole}
\end{figure}

\subsection{Half cheetah - beyond small sample sizes}

We now consider the performance of our method when applied on the half cheetah, beyond the sample-efficiency domain. In Figure~\ref{fig:cheetah_longer} we present the reward obtained using our DGP-MPC approach for the half cheetah environment over a larger number of environment observations.

\begin{figure}[ht]
\centering
\includegraphics[width=\textwidth]{Images/supp/plots/half_cheetah_longer.png}
\caption{The half cheetah reward with increasing numbers of layers. Each case used the Mat\'{e}rn-3/2 kernel across all layers..
}
\label{fig:cheetah_longer}
\end{figure}

\subsubsection*{Modified cartpole} See Figure~\ref{fig:kernels_modified_cartpole}

\begin{figure}[ht]
\centering
\subfigure[Squared exponential kernel on all layers.]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/modified_cartpole_sexp_sexp.png}
    \label{fig:}
}
\subfigure[Squared exponential on first layers, matern $\frac{5}{2}$ kernel on final layer.]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/modified_cartpole_sexp_m52.png}
    \label{fig:}
}
\caption{Modified cartpole environment 'as seen' reward.
}
\label{fig:kernels_modified_cartpole}
\end{figure}

\newpage
\subsubsection*{Reacher}

\newpage
\paragraph{Half cheetah} The case of no reward shaping with different kernel choices is given in~ Figure~\ref{fig:kernels_halfcheetah_noshaping}. The results when adding reward shaping are shown in Figure~\ref{fig:kernels_halfcheetah_shaping}.

\begin{figure}[ht]
\centering
\subfigure[Squared exponential kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_sexp.png}
    \label{fig:}
}
\subfigure[Matern $\frac{5}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern52.png}
    \label{fig:}
}
\subfigure[Matern $\frac{3}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern32.png}
    \label{fig:}
}
\subfigure[Matern $\frac{1}{2}$ kernel on all layers.]{
    \includegraphics[width=0.8\textwidth]{neurips/supp/plots/halfcheetah_matern12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, with no reward shaping. \textbf{Different kernel} choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_noshaping}
\end{figure}

\begin{figure}[ht]
\centering
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, with reward shaping]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/halfcheetah_shaped_matern12.png}
    \label{fig:}
}
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, using the lower credible bound (minus one deviation) of the pushforward measure on the reward]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/halfcheetah_LCB_m12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, and with \textbf{reward shaping}. Different kernel choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_shaping}
\end{figure}

\begin{figure}[ht]
\centering
\subfigure[Matern $\frac{1}{2}$ kernel on all layers, without reward shaping but with noisy observations]{
    \includegraphics[width=\textwidth]{neurips/supp/plots/half_cheetah_noisy_matern12.png}
    \label{fig:}
}
\caption{Half-cheetah environment 'as seen' reward, and without shaping and with \textbf{noise added to observations}. Different kernel choices encode prior knowledge on the dynamic system's smoothness.
}
\label{fig:kernels_halfcheetah_shaping}
\end{figure}

\subsubsection*{Modified cartpole} We study the effect of deep Gaussian process depth and kernel choice in the modified cartpole environment. Towards this end we consider three situations: a) an isolated non-stationary case of the cart hitting the wall b) an isolated stationary case of the cart travelling along the rail c) a hybrid case containing transitions of both. In each case the trajectory plotted is an open loop, where at each index the next sample is conditional on the previous sampled state, with no feedback. Here, steps are made on a finer time-scale to reduce discretisation errors.

\paragraph{Case a} We train our DGP model on $50$ sequences of length $10$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of positive actions are then uniformly sampled in $[0,3]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=0.75$ from the initial state mean, and comparing to the ground truth.  See Figure~\ref{fig:kernels_modified_cartpole_dynamics_crash}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth, height=3cm]{Images/supp/plots/dynamics/crash/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case a}. Analysis on the effect of kernel choice and depth when modelling the non-stationary contact dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_crash}
\end{figure}

\paragraph{Case b} We train our DGP model on $50$ sequences of length $10$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of negative actions are then uniformly sampled in $[-3,0]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=-1$ from the initial state mean, and comparing to the ground truth. See Figure~\ref{fig:kernels_modified_cartpole_dynamics_no_crash}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/no_crash/new/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case b}. Analysis on the effect of kernel choice and depth when modelling the stationary dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_no_crash}
\end{figure}

\paragraph{Case c} We train our DGP model on $50$ sequences of length $20$. Each sequence begins at a random initial point, under the same distribution used in our control experiments. Sequences of actions are then uniformly sampled in $[-3,3]$ and the resulting transitions are used to train our DGP model. This is then tested by sampling the trajectory under a sequence of fixed actions with $a=0.75$ from the initial state mean, and comparing to the ground truth. See Figure~\ref{fig:kernels_modified_cartpole_dynamics_mixed}.

\begin{figure}[ht!]
\centering
\subfigure[SEXP kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l1y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l2y3.png}
    \label{fig:}
}
\subfigure[SEXP kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/SEXP_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern52 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l1y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l2y3.png}
    \label{fig:}
}
\subfigure[Matern52 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern52_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern32 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l1y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l2y3.png}
    \label{fig:}
}
\subfigure[Matern32 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern32_l3y3.png}
    \label{fig:}
} \\
\subfigure[Matern12 kernel, $l=1$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l1y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=2$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l2y3.png}
    \label{fig:}
}
\subfigure[Matern12 kernel, $l=3$.]{
    \includegraphics[width=0.3\textwidth]{Images/supp/plots/dynamics/mixed/new/Matern12_l3y3.png}
    \label{fig:}
} \\
\caption{\textbf{Case c}. Analysis on the effect of kernel choice and depth when modelling both the stationary and non-stationary dynamics in the modified cartpole example.
}
\label{fig:kernels_modified_cartpole_dynamics_mixed}
\end{figure}

redpajama_set_name: RedPajamaArXiv


text:
Anti-poverty advocates call for BC Hydro to implement an electricity affordability program for BC's poor
admin August 15, 2016 June 14, 2017 BC Hydro, BC Utilities Commission, Energy Poverty, Energy Utility, News, Press Release
For Immediate Release | BCPIAC
August 15, 2016 (Vancouver) Starting on August 16, 2016, the BC Utilities Commission (BCUC) will hold a public hearing to review BC Hydro's Rate Design Application (RDA). In this process, the BCUC will hear evidence and submissions from BC Hydro and intervener groups and determine rate structures and terms and conditions of service for residential, business and industrial customers.
The BC Public Interest Advocacy Centre (BCPIAC) is intervening on behalf of seven organizations in this rare opportunity to ask the BCUC to order BC Hydro to implement programs for low income residential ratepayers including a discounted rate for electricity, low income customer rules, and a crisis intervention fund.
BCPIAC is seeking these programs because low income customers are struggling to pay BC Hydro's endlessly increasing electricity rates. In support of these proposals, BCPIAC has presented evidence from two expert witnesses, and a number of advocates and low income ratepayers.
BC Hydro has increased residential electricity rates by 50% in the last 10 years, and is on track to increase them by over 30% in the next eight years. Rates are projected to continue to rise significantly in future years as BC Hydro proceeds with multi-billion dollar projects such as the Site C dam which have been exempted from a full public review by the BCUC, and is currently estimated to cost over $9 billion.
BC Hydro's rate increases have grossly outpaced increases in income for low income British Columbians. For example, BC social assistance rates have been frozen since 2007 at $610 per month for basic assistance and $906 for disability assistance, and in the last 10 years the BC general minimum wage has only increased by $2.45 an hour.
Keith Simmonds, a Minister with Duncan United Church in Duncan BC, who provided evidence in support of BCPIAC's proposals, has seen an increasing number of people who are either facing disconnection or have already been cut off because they cannot afford to pay their BC Hydro bills. "Inability to pay rising Hydro costs has a huge impact on low income people, especially seniors on fixed pensions and people who rely on the very low social assistance rates that have been frozen since 2007" Simmonds said. "Many people have to choose between paying rent, purchasing food, and paying their electricity bill."
"About 10% of BC Hydro residential customers live below Statistics Canada's Low Income Cut-off, which translates to about 170,000 customers", said Sarah Khan, one of the lawyers at BCPIAC who is bringing this issue to the BCUC. Khan added that "BC Hydro's rates are increasing much faster than incomes for low income people, and our proposals will help to mitigate the impact of BC Hydro's never-ending rate increases".
BC Hydro offers no rates or terms and conditions that specifically apply to low income customers. The only programs available to these customers are energy saving kits and in more limited cases, energy efficiency audits and certain home upgrades. While these programs are important, they are not offsetting BC Hydro's rate increases.
BCPIAC's low income proposals were developed by Roger Colton, an expert in low income rate design from the United States, who will testify as an expert witness during the public hearing. Mr. Colton will testify that measures such as a discounted electricity rate can simultaneously address affordability concerns, improve cost reflectivity in rates, and improve the efficiency of BC Hydro's operations.
BCPIAC will also be calling Seth Klein, Director of the Canadian Centre for Policy Alternatives, BC Office, to testify about the extent and profile of poverty in BC and the difficulties low income people in BC have paying for the basic necessities of life including residential electricity. Mr. Klein's evidence will provide context for the need for mitigating measures like BCPIAC's proposed bill affordability programs.
BCPIAC is representing the following groups in this proceeding: Active Support Against Poverty, BC Old Age Pensioners' Organization, BC Poverty Reduction Coalition, Council of Senior Citizens' Organizations of BC, Disability Alliance BC, Together Against Poverty Society, and Tenant Resource & Advisory Centre.
The hearing will take place on August 16-18 and 23-24, 2016 in downtown Vancouver.
Tagged on: Electricity Affordability Erin Sarah
← BC woman launches challenge for right to legal representation in Mental Health Act detentions
BC man wins important human rights victory for Deaf people in need of welfare →
BC Hydro
BC Utilities Commission
Fortis BC
Racialized & Immigrant Communities

redpajama_set_name: RedPajamaCommonCrawl