PILCO, PETS, Ensembles, and MBPO › Probabilistic dynamics and data efficiency

The motivating problem is data efficiency. On some benchmark tasks, model-free policy-gradient and actor-critic baselines require millions of environment steps, as illustrated by the MBPO benchmark comparison. A model-based method adds a supervised learning problem: each observed $(s,a,s')$ transition trains a predictor of local dynamics. Under an appropriate smoothness prior or other inductive bias, that observation can inform predictions at nearby state-action pairs. Once the model is useful, planning or policy learning can consume synthetic transitions without further environment interaction. PILCO was an influential early demonstration that this strategy could achieve exceptional data efficiency on continuous-control benchmarks, as a later probabilistic-MPC comparison also notes.

The contrast is not that model-free methods reduce every transition to a few bits or ignore its next state. Value and policy algorithms may use states, actions, rewards, next states, temporal-difference errors, advantages, or distributional targets. The distinction is that a dynamics learner uses the transition as direct supervision for a reusable transition model. Generalization from one transition to a neighborhood comes from the model class and its fitted inductive bias, not from the tuple alone. Queries to the learned model then produce synthetic transitions, not additional real experience. Suitable priors or smoothness assumptions may improve data efficiency by helping a dynamics model generalize from limited data. Propagating uncertainty and choosing an appropriate objective can limit exposure to model bias when their assumptions hold. None of these choices guarantees safe extrapolation.

State uncertainty and dynamics uncertainty are also distinct rather than mutually exclusive. A state estimator may maintain a belief over the current hidden state while using known, estimated, or jointly inferred dynamics. PILCO learns uncertain dynamics and propagates uncertain state distributions through them. In both cases dispersion records uncertainty conditional on the model assumptions, and new observations update the relevant belief. A fuller belief-space controller can maintain uncertainty over state and model parameters or functions at the same time.

Suppose you fit a deterministic model $\hat{f}(s,a)$ and optimize a controller against its predicted trajectories. Away from observed data, the predictions are weakly constrained by evidence and are determined largely by architecture, regularization, and optimization. A sufficiently aggressive controller search can then select actions that benefit from those off-support errors. This is a risk of optimizing one fitted hypothesis, not proof that every point-estimate controller must fail.

The failure can arise structurally when three conditions coincide: the model is inaccurate in reachable regions, the optimizer can search those regions, and the objective rewards the resulting prediction errors. Better coverage, regularization, constraints, robust objectives, or benign task structure can prevent it. Representing model uncertainty can help, but only if the decision rule uses that information appropriately; posterior averaging alone is not a universal safety mechanism.

Formally, let the true dynamics be the conditional distribution $p(s_{t+1} \mid s_t, a_t)$. A probabilistic model assigns a distribution $q_\phi(s_{t+1} \mid s_t, a_t)$ parameterized by $\phi$. In the Bayesian view, $\phi$ itself is uncertain, giving a posterior $p(\phi \mid \mathcal{D})$ after seeing data $\mathcal{D}$. To predict the next state, we average over all models weighted by how plausible each is under the data, which means marginalizing over the uncertain model parameters:

$$q(s_{t+1} \mid s_t, a_t, \mathcal{D}) = \int q_\phi(s_{t+1} \mid s_t, a_t) \, p(\phi \mid \mathcal{D}) \, d\phi$$

where:

• $s_t$: the state at time $t$
• $a_t$: the action taken at time $t$
• $\phi$: the parameters of the probabilistic dynamics model
• $q_\phi(s_{t+1} \mid s_t, a_t)$: the model's predicted distribution over the next state given the current state and action
• $p(\phi \mid \mathcal{D})$: the posterior distribution over model parameters after observing data $\mathcal{D}$
• $\mathcal{D}$: the dataset of observed transitions

The integral averages predictions from all parameter settings according to their posterior plausibility. Under an identifiable, well-specified model with an appropriate prior and likelihood, informative data can concentrate the posterior and reduce latent predictive uncertainty. In data-sparse regions, uncertainty may remain larger. Neither behavior is automatic under misspecification, and integrating over parameters is not itself a risk penalty: it defines a posterior-predictive distribution.

Consider a dynamics model with one unknown gain. Data near a well-explored operating point may narrow the posterior over that gain, while predictions elsewhere may remain sensitive to the prior and model assumptions. Bayesian averaging prevents the planner from evaluating only the single most flattering gain, but a risk-neutral expected return can still favor an uncertain action if its posterior mean payoff is high. A controller that should avoid optimistic tails needs a suitable nonlinear utility, a lower-confidence or worst-case objective, CVaR, or another explicit risk criterion.

An operational test requires controlled conditions. Repeating the same fully specified physical transition in a stationary regime can reveal intrinsic process variation, although unmodeled hidden state can confound the comparison. If the prediction target is the latent physical transition, stochastic sensor noise can look like process noise; for a measured-outcome target, residual stochastic sensor noise is an aleatoric component, while systematic sensor bias is not automatically one. Disagreement among models retrained with different samples or optimization randomness is a practical epistemic proxy, not a pure measurement: it contains sampling and training variability. Additional informative local data may narrow disagreement when the model class and training procedure turn those observations into more similar local fits, but the change is not guaranteed or monotone. The two uncertainties call for different responses. Irreducible variability may require expected-outcome or robust planning; reducible ignorance may justify targeted data collection. Models that expose only a collapsed total predictive variance make that distinction harder to use.

Let $s\in\mathbb{R}^{d_s}$ and $a\in\mathbb{R}^{d_a}$, and set $x=(s,a)\in\mathbb{R}^d$ with $d=d_s+d_a$. Let $f(s,a)$ denote one scalar component of a latent dynamics function. Assume $f\sim\mathcal{GP}(0,k)$ and observations $y_i=f(x_i)+\epsilon_i$, where the $\epsilon_i\sim\mathcal{N}(0,\sigma_n^2)$ are mutually independent and independent of $f$. A positive-semidefinite kernel $k(x,x')$ defines the covariance of any finite collection of function values. Assume $\sigma_n^2>0$, so $K+\sigma_n^2I$ is invertible even if training inputs coincide. Given $n$ training inputs $X$ and scalar-output targets $y\in\mathbb{R}^n$, the posterior distribution of the latent function value at $x_*$ is Gaussian with

$$\begin{aligned} \mu_f(x_*) &= k_*^\top (K+\sigma_n^2I)^{-1}y, \\ v_f(x_*) &= k_{**}-k_*^\top (K+\sigma_n^2I)^{-1}k_*. \end{aligned}$$

Here $K\in\mathbb{R}^{n\times n}$ has entries $K_{ij}=k(x_i,x_j)$, $k_*\in\mathbb{R}^n$ contains the covariances between the query and all training inputs, $I$ is the $n\times n$ identity, and $k_{**}=k(x_*,x_*)$. The displayed $v_f$ is latent-function variance. The predictive variance of a new noisy observation is $v_y(x_*)=v_f(x_*)+\sigma_n^2$ under this homoscedastic likelihood.

For fixed inputs and fixed kernel and noise hyperparameters, the posterior mean is linear in the targets, and its effective weights are the entries of $(K+\sigma_n^2I)^{-1}k_*$. They are obtained by solving the full kernel system, generally depend jointly on the training inputs, and can even be negative. A training point's query-to-point similarity alone therefore does not determine its weight; the full Gram matrix, query covariance vector, and noise variance matter. With a local squared-exponential kernel, nearby isolated observations often exert more influence, but that is an intuition rather than a coefficient-by-coefficient guarantee.

The latent variance depends on the input locations and fitted kernel and noise parameters, not on the observed target values. It generally decreases near informative observations, but with noisy data it need not reach zero even at a training input. For a stationary kernel whose covariance decays with distance, it approaches the prior variance when all query-to-training covariances vanish. These statements remain conditional on the kernel, likelihood, hyperparameters, and model specification: a misspecified GP can still be confidently wrong. With PILCO's squared-exponential kernel and differentiable moment, policy, and cost expressions, the resulting approximate predictive moments are differentiable in the quantities used for policy optimization.

The kernel encodes assumptions about smoothness and characteristic variation scales. PILCO uses a squared-exponential kernel,

$$k(x,x')=\sigma_f^2\exp\!\left(-\frac{1}{2}(x-x')^\top\Lambda^{-1}(x-x')\right),$$

where $x$ and $x'$ concatenate state and action, $\sigma_f^2$ is the signal variance, and

$$\Lambda=\operatorname{diag}(\ell_1^2,\ldots,\ell_d^2)$$

is a $d\times d$ matrix containing the squared characteristic length scales. A shorter $\ell_j$ permits covariance to decay more rapidly along input dimension $j$; a longer one favors slower variation. PILCO fits these and other GP hyperparameters with gradient-based log-marginal-likelihood optimization. Its released GP trainer also adds numerical-stability penalties for extreme length scales and signal-to-noise ratios.

The squared-exponential prior favors extremely smooth functions: its sample paths are infinitely differentiable under the usual construction, and nearby inputs receive high prior covariance. The length scales control typical variation, but they do not impose a deterministic bound on a function's derivative; derivative distributions still have unbounded Gaussian support. This strong smoothness prior can work well for smoothly varying mechanics, but stationary squared-exponential dynamics can fit contact, stiction, and discrete switching poorly because those regimes conflict with the prior's smoothness and stationarity assumptions.

Suppose the state at time $t$ is approximated by $s_t\sim\mathcal{N}(\mu_t,\Sigma_t)$ and the deterministic feedback policy is $a_t=\pi(s_t;\theta)$. A generic additive-noise transition can be written as

$$s_{t+1}=f(s_t,\pi(s_t;\theta))+\epsilon_t.$$

PILCO specifically fits independent GPs to state increments $\Delta_t=s_{t+1}-s_t$ and reconstructs the next state by adding the predicted increment to $s_t$. The fitted target noise and the GP's latent-function uncertainty must therefore be distinguished from any physical process-noise interpretation.

Passing an uncertain Gaussian state through a nonlinear policy and GP predictive model generally produces a non-Gaussian distribution. Non-Gaussianity alone does not make moments intractable, but the required integrals are not available in closed form for arbitrary kernels and policies. PILCO exploits compatible squared-exponential kernels and policy forms to compute the needed first and second moments analytically, then projects the result back to a Gaussian. Thus the integrals used for each update are analytic while the repeated Gaussian moment-matching inference is approximate. The PILCO authors also note that their rollout does not retain temporal correlation in model uncertainty: it treats that uncertainty similarly to temporally uncorrelated noise, which can underestimate it. This caveat does not erase the benefit of uncertainty-aware evaluation, but it rules out reading the rollout as an exact joint distribution over a fixed uncertain dynamics function.

The resulting $\mu_{t+1}$ and $\Sigma_{t+1}$ are differentiable functions of $\mu_t$, $\Sigma_t$, and $\theta$, so gradients of predicted cost can flow back to the policy parameters. Many alternative kernels or policy nonlinearities would require new moment derivations or numerical approximation; the limitation is tractability for this pipeline, not a claim that every other kernel has intractable moments.

With the moment-matching machinery in place, we can state the control objective: choose policy parameters that minimize the total cost the agent expects to incur along a finite horizon. Let $\ell(s)$ be a cost function (PILCO uses a saturating cost; more on that shortly). The expected cost over a finite horizon of length $T$ is

$$J(\theta) = \sum_{t=1}^{T} \mathbb{E}_{s_t \sim \mathcal{N}(\mu_t, \Sigma_t)}\!\big[\ell(s_t)\big]$$

where:

• $J(\theta)$: the expected cost over the horizon, as a function of the policy parameters $\theta$
• $T$: the finite planning horizon
• $\ell(s)$: the per-state cost function
• $\mathcal{N}(\mu_t, \Sigma_t)$: the Gaussian predictive distribution over states at time $t$, obtained by moment matching
• $\mu_t, \Sigma_t$: the mean and covariance of that predictive distribution

Because the matched moments are differentiable in $\theta$, PILCO computes $\nabla_\theta J$ without Monte Carlo trajectory gradients. The expected cost incorporates the GP's predictive uncertainty and can mitigate model exploitation, but approximate inference and model misspecification prevent an unconditional guarantee.

The objective is a finite, undiscounted cost sum. PILCO does not optimize one fixed open-loop action sequence: it propagates the state distribution induced by the same parameterized feedback policy over the horizon. Each step depends on the previous predictive moments and on $\theta$, so the gradient follows the same dependency pattern as backpropagation through time. The recurrent operation is an analytic moment update under an approximate Gaussian state distribution rather than a deterministic update of one state.

One more design choice deserves attention. A positive-definite quadratic distance cost grows without bound as the state moves away from the target. Under a Gaussian predictive state distribution, its expectation contains both squared mean error and a term proportional to predictive covariance. Large uncertainty can therefore strongly influence the objective, although unstable gradients or literal tail dominance do not follow automatically. For normalized or otherwise commensurate state coordinates, an isotropic PILCO-style bounded saturating cost is

$$\ell(s)=1-\exp\!\left(-\frac{\|s-s_{\text{target}}\|^2}{2\sigma_c^2}\right),$$

where $\sigma_c$ sets the distance scale. The exponential term is largest at the target and decays with distance, so the cost moves smoothly from low cost near the target to a high but bounded cost far away. Specifically, $\ell(s)=0$ at the target and approaches $1$ as the distance grows.

The bounded cost limits how much very distant predictive states can contribute. Once the distance is large relative to $\sigma_c$, moving still farther changes the cost only slightly. Compared with an unbounded quadratic, this reduces the influence of remote uncertain states and emphasizes distinctions near the task-relevant target scale. It does not prove that those remote states are impossible or irrelevant; it is an objective-design choice that controls their influence on policy optimization.

The PILCO experiments reported control learning in a small number of trials on several continuous-control systems. The result made the method an influential example of data-efficient policy search. Its restrictions matter just as much:

• Computational cost scales poorly. Standard exact GP training with a dense, unstructured covariance matrix costs $O(n^3)$ time per Cholesky factorization and $O(n^2)$ memory, and moment propagation is expensive. The workable point count is implementation- and problem-dependent; hundreds to low thousands is a rough historical scale, not a hard limit.
• The prior and analytic forms are restrictive. Stationary squared-exponential kernels are poorly matched to discontinuities and switching, while analytic moment calculations require compatible kernels and policy forms.
• Repeated Gaussian projection can degrade. Long-horizon propagation repeatedly compresses a potentially non-Gaussian distribution into its first two moments.

PETS directly addresses GP scalability by using probabilistic neural ensembles and Monte Carlo propagation. MBPO addresses a different problem: bias from using learned dynamics to generate policy-training data over long horizons. It would be misleading to present MBPO's short branches as a direct remedy for PILCO's particular approximation errors.