---
abstract: 'The goal of reinforcement learning (RL) is to let an agent learn an optimal control policy in an unknown environment so that future expected rewards are maximized. The model-free RL approach directly learns the policy based on data samples. Although using many samples tends to improve the accuracy of policy learning, collecting a large number of samples is often expensive in practice. On the other hand, the model-based RL approach first estimates the transition model of the environment and then learns the policy based on the estimated transition model. Thus, if the transition model is accurately learned from a small amount of data, the model-based approach can perform better than the model-free approach. In this paper, we propose a novel model-based RL method by combining a recently proposed model-free policy search method called *policy gradients with parameter-based exploration* and the state-of-the-art transition model estimator called *least-squares conditional density estimation*. Through experiments, we demonstrate the practical usefulness of the proposed method.'
author:
- |
    Syogo Mori$^1$, Voot Tangkaratt$^1$, Tingting Zhao$^1$\
    Jun Morimoto$^2$, and Masashi Sugiyama$^1$\
    $^1$Tokyo Institute of Technology, Japan.\
    {mori@sg., voot@sg., tingting@sg., sugi@}cs.titech.ac.jp\
    $^2$ATR Computational Neuroscience Labs, Japan\
    xmorimo@atr.jp
title: |
     Model-Based Policy Gradients\
    with Parameter-Based Exploration\
    by Least-Squares Conditional Density Estimation 
---

Introduction {#sec:introduction}
============

*Reinforcement learning* (RL) is a framework to let an agent learn an optimal control policy in an unknown environment so that expected future rewards are maximized [@JAIR:Kaelbling+Littman+Moore:1996]. The RL methods developed so far can be categorized into two types: *Policy iteration* where policies are learned based on value function approximation [@book:Sutton+Barto:1998; @JMLR:Lagoudakis+Parr:2003] and *policy search* where policies are learned directly to maximize expected future rewards [@mach:Williams:1992; @nc:Dayan+Hinton:1997; @SutMcASinMan00; @nips02-CN11; @NN:Sehnke+etal:2010; @NC:Zhao+etal:2013].

Policy Iteration Vs. Policy Search
----------------------------------

A value function represents expected future rewards as a function of a state or a state and an action. In the policy iteration framework, approximation of the value function for the current policy and improvement of the policy based on the learned value function are iteratively performed until an optimal policy is found. Thus, accurately approximating the value function is a challenge in the value function based approach. So far, various machine learning techniques have been employed for better value function approximation, such as least-squares approximation [@JMLR:Lagoudakis+Parr:2003], manifold learning [@AutoRobo:Sugiyama+etal:2008], efficient sample reuse [@NN:Hachiya+etal:2009], active learning [@NN:Akiyama+etal:2010], and robust learning [@IEICE:Sugiyama+etal:2010b].

However, because policy functions are learned indirectly via value functions in policy iteration, improving the quality of value function approximation does not necessarily yield a better policy function. Furthermore, because a small change in value functions can cause a big change in policy functions, it is not safe to use the value function based approach for controlling expensive dynamic systems such as a humanoid robot. Another weakness of the value function approach is that it is difficult to handle continuous actions because a maximizer of the value function with respect to an action needs to be found for policy improvement.

On the other hand, in the policy search approach, policy functions are determined so that expected future rewards are directly maximized. A popular policy search method is to update policy functions via gradient ascent. However, a classic policy gradient method called REINFORCE [@mach:Williams:1992] tends to produce gradient estimates with large variance, which results in unreliable policy improvement [@IROS:Peters+Schaal:2006]. More theoretically, it was shown that the variance of policy gradients can be proportional to the length of an agent’s trajectory, due to the stochasticity of policies [@NN:Zhao+etal:2012]. This can be a critical limitation in RL problems with long trajectories.

To cope with this problem, a novel policy gradient method called *policy gradients with parameter-based exploration* (PGPE) was proposed [@NN:Sehnke+etal:2010]. In PGPE, deterministic policies are used to suppress irrelevant randomness and useful stochasticity is introduced by drawing policy parameters from a prior distribution. Then, instead of policy parameters, hyper-parameters included in the prior distribution are learned from data. Thanks to this prior-based formulation, the variance of gradient estimates in PGPE is independent of the length of an agent’s trajectory [@NN:Zhao+etal:2012]. However, PGPE still suffers from an instability problem in small sample cases. To further improve the practical performance of PGPE, an efficient sample reuse method called *importance-weighted PGPE* (IW-PGPE) was proposed recently and demonstrated to achieve the state-of-the-art performance [@NC:Zhao+etal:2013].

Model-Based Vs. Model-Free
--------------------------

The RL methods reviewed above are categorized into the *model-free* approach, where policies are learned without explicitly modeling the unknown environment (i.e., the transition probability of the agent in the environment). On the other hand, an alternative approach called the *model-based* approach explicitly models the environment in advance and uses the learned environment model for policy learning [@Wang:2003; @Deisenroth:2011]. In the model-based approach, no additional sampling cost is necessary to generate artificial samples from the learned environment model.

The model-based approach is particularly advantageous in the policy search scenario. For example, given a fixed budget for data collection, IW-PGPE requires us to determine the *sampling schedule* in advance. More specifically, we need to decide, e.g., whether we gather many samples in the beginning or only a small batch of samples are collected for a longer period. However, optimizing the sampling schedule in advance is not possible without strong prior knowledge. Thus, we need to just blindly design the sampling schedule in practice, which can cause significant performance degradation. On the other hand, the model-based approach does not suffer from this problem because we can draw as many trajectory samples as we want from the learned transition model without additional sampling costs.

Another advantage of the model-based approach lies in *baseline subtraction*. In the gradient-based policy search methods such as REINFORCE and PGPE, subtraction of a baseline from a gradient estimate is a vital technique to reduce the estimation variance of policy gradients [@IROS:Peters+Schaal:2006; @NC:Zhao+etal:2013]. If the baseline is estimated from samples that are statistically independent of samples used for the estimation of policy gradients, variance reduction can be carried out without increasing the estimation bias. However, such independent samples are not available in practice (if available, they should be used for policy gradient estimation), and thus variance reduction by baseline subtraction is practically performed at the expense of bias increase. On the other hand, in the model-based scenario, we can draw as many trajectory samples as we want from the learned transition model without additional sampling costs. Therefore, two statistically independent sets of samples can be generated and they can be separately used for policy gradient estimation and baseline estimation.

Transition Model Learning by Least-Squares Conditional Density Estimation
-------------------------------------------------------------------------

If the unknown environment is accurately approximated, the model-based approach can fully enjoy all the above advantages. However, accurately estimating the transition model from a limited amount of trajectory data in multi-dimensional continuous state and action spaces is highly challenging. Although the model-based method that does not require an accurate transition model was developed [@Abbeel:2006], it is only applicable to deterministic environments, which significantly limits its range of applications in practice. On the other hand, a recently proposed model-based policy search method called PILCO [@Deisenroth:2011] learns a probabilistic transition model by the Gaussian process (GP) [@book:Rasmussen+Williams:2006], and explicitly incorporates long-term model uncertainty. However, PILCO requires states and actions to follow Gaussian distributions and the reward function to be a particular exponential form to ensure that the policy evaluation is performed in a closed form and policy gradients are computed analytically for policy improvement. These strong requirements make PILCO practically restrictive. To overcome such limitations of existing approaches, we propose a highly practical policy-search algorithm by extending the model-free PGPE method to the model-based scenario. In the proposed model-based PGPE (M-PGPE) method, the transition model is learned by the state-of-the-art non-parametric conditional density estimator called *least-squares conditional density estimation* (LSCDE) [@IEICE:Sugiyama+etal:2010a], which has various superior properties: It can directly handle multi-dimensional inputs and outputs, it was proved to achieve the optimal convergence rate [@ML:Kanamori+etal:2012], it has high numerical stability [@ML:Kanamori+etal:2012b], it is robust against outliers [@AISM:Sugiyama+etal:2012], its solution can be analytically and efficiently computed just by solving a system of linear equations [@JMLR:Kanamori+etal:2009], and generating samples from the learned conditional density is straightforward. Through experiments, we demonstrate that the proposed M-PGPE method is a promising approach. The rest of this paper is structured as follows. In Section \[sec:MF\], we formulate the RL problem and review model-free RL methods including PGPE. We then propose the model-based PGPE method in Section \[sec:MB\], and experimentally demonstrate its usefulness in Section \[sec:experiments\]. Finally, we conclude in Section \[sec:conclusion\].

Problem Formulation and Model-Free Policy Search {#sec:MF}
================================================

In this section, we first formulate our RL problem and review existing model-free policy search methods.

Formulation {#subsec:formulation}
-----------

Let us consider a Markov decision problem consisting of the following elements:

-   ${{\mathcal S}}$: A set of continuous states.

-   ${{\mathcal A}}$: A set of continuous actions.

-   ${p({\boldsymbol{s}})}$: The (unknown) probability density of initial states.

-   ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$: The (unknown) conditional probability density of visiting state ${\boldsymbol{s}}'$ from state ${\boldsymbol{s}}$ by action ${a}$.

-   ${{R}({\boldsymbol{s}},{a},{\boldsymbol{s}}')}$: The immediate reward function for the transition from ${\boldsymbol{s}}$ to ${\boldsymbol{s}}'$ by ${a}$.

Let ${p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}$ be a policy of an agent parameterized by ${\boldsymbol{\theta}}$, which is the conditional probability density of taking action ${a}$ at state ${\boldsymbol{s}}$. Let $$\begin{aligned}
  {h}:=[{\boldsymbol{s}}_1,{a}_1,\ldots,{\boldsymbol{s}}_{T},{a}_{T},{{\boldsymbol{s}}_{{T}+1}}]\end{aligned}$$ be a history, which is a sequence of states and actions with finite length ${T}$ generated as follows: First, the initial state ${\boldsymbol{s}}_1$ is determined following the initial-state probability density ${p({\boldsymbol{s}})}$. Then action ${a}_1$ is chosen following policy ${p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}$, and next state ${\boldsymbol{s}}_2$ is determined following the transition probability density ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$. This process is repeated ${T}$ times.

Let ${{r}({h})}$ be the return for history ${h}$, which is the discounted sum of future rewards the agent can obtain: $$\begin{aligned}
{{r}({h})}:=\sum^{{T}}_{{t}=1}{\gamma}^{{t}-1}{{R}({{\boldsymbol{s}}_{{t}}},{{a}_{{t}}},{{\boldsymbol{s}}_{{t}+1}})},
$$ where ${\gamma}\in(0,1]$ is a discount factor. The expected return is given by $$\begin{aligned}
{J({\boldsymbol{\theta}})}:=\int{{r}({h})}{p({h}|{\boldsymbol{\theta}})}\mathrm{d}{h},\end{aligned}$$ where ${p({h}|{\boldsymbol{\theta}})}$ is the probability density of observing history ${h}$: $$\begin{aligned}
{p({h}|{\boldsymbol{\theta}})}={p({{\boldsymbol{s}}_{1}})}\prod_{{t}=1}^{{T}}
{p({{a}_{{t}}}|{{\boldsymbol{s}}_{{t}}},{\boldsymbol{\theta}})}{p({{\boldsymbol{s}}_{{t}+1}}|{{\boldsymbol{s}}_{{t}}},{{a}_{{t}}})}.
$$

The goal of RL is to find optimal policy parameter ${{\boldsymbol{\theta}}^{*}}$ that maximizes the expected return ${J({\boldsymbol{\theta}})}$: $$\begin{aligned}
{{\boldsymbol{\theta}}^{*}}:={\mathop{\mathrm{argmax}}}_{{\boldsymbol{\theta}}}{J({\boldsymbol{\theta}})}.\end{aligned}$$

REINFORCE {#subsec:reinforce}
---------

REINFORCE [@mach:Williams:1992] is a classic method for learning the policy parameter ${\boldsymbol{\theta}}$ via gradient ascent: $$\begin{aligned}
{\boldsymbol{\theta}}\leftarrow{\boldsymbol{\theta}}+{\varepsilon}{\nabla_{{\boldsymbol{\theta}}}}{J({\boldsymbol{\theta}})},\end{aligned}$$ where ${\varepsilon}>0$ denotes the learning rate and ${\nabla_{{\boldsymbol{\theta}}}}{J({\boldsymbol{\theta}})}$ denotes the gradient of ${J({\boldsymbol{\theta}})}$ with respect to ${\boldsymbol{\theta}}$.

The gradient ${\nabla_{{\boldsymbol{\theta}}}}{J({\boldsymbol{\theta}})}$ can be expressed as $$\begin{aligned}
{\nabla_{{\boldsymbol{\theta}}}}{J({\boldsymbol{\theta}})}&=\int{{r}({h})}{\nabla_{{\boldsymbol{\theta}}}}{p({h}|{\boldsymbol{\theta}})}\mathrm{d}{h}\\
&=\int{{r}(h)}{p({h}|{\boldsymbol{\theta}})}{\nabla_{{\boldsymbol{\theta}}}}\log{p({h}|{\boldsymbol{\theta}})}\mathrm{d}{h}\\
&=\int{{r}(h)}{p({h}|{\boldsymbol{\theta}})}\sum_{{t}=1}^{{T}}
{\nabla_{{\boldsymbol{\theta}}}}\log{p({{a}_{{t}}}|{{\boldsymbol{s}}_{{t}}},{\boldsymbol{\theta}})}\mathrm{d}{{h}},\end{aligned}$$ where we used $$\begin{aligned}
{\nabla_{{\boldsymbol{\theta}}}}{p({h}|{\boldsymbol{\theta}})}={p({h}|{\boldsymbol{\theta}})}{\nabla_{{\boldsymbol{\theta}}}}\log{p({h}|{\boldsymbol{\theta}})}.\end{aligned}$$

In the above expression, the probability density of histories, ${p({h}|{\boldsymbol{\theta}})}$, is unknown. Suppose that we are given ${N}$ roll-out samples $\{{{h}_{{n}}}\}_{{n}=1}^{{N}}$ for the current policy, where $$\begin{aligned}
  {{h}_{{n}}}=[{{\boldsymbol{s}}^{{n}}_{1}},{{a}^{{n}}_{1}},\ldots,{{\boldsymbol{s}}^{{n}}_{{T}}},
  {{a}^{{n}}_{{T}}},{{\boldsymbol{s}}^{{n}}_{{T}+1}}].\end{aligned}$$ Then the expectation over ${p({h}|{\boldsymbol{\theta}})}$ can be approximated by the empirical average over the samples $\{{{h}_{{n}}}\}_{{n}=1}^{{N}}$, i.e., an empirical approximation of the gradient ${\nabla_{{\boldsymbol{\theta}}}}{J({\boldsymbol{\theta}})}$ is given by $$\begin{aligned}
  {\nabla_{{\boldsymbol{\theta}}}}{\widehat{J}({\boldsymbol{\theta}})}:=\frac{1}{{N}}\sum_{{n}=1}^{{N}}{{r}({{h}_{{n}}})}
  \sum_{{t}=1}^{{T}}{\nabla_{{\boldsymbol{\theta}}}}\log{p({{a}^{{n}}_{{t}}}|{{\boldsymbol{s}}^{{n}}_{{t}}},{\boldsymbol{\theta}})}.\end{aligned}$$ It is known [@IROS:Peters+Schaal:2006] that the variance of the above gradient estimator can be reduced by subtracting the baseline ${b}$: $$\begin{aligned}
  {\nabla_{{\boldsymbol{\theta}}}}{\widehat{J}^{{b}}({\boldsymbol{\theta}})}:=\frac{1}{{N}}\sum_{{n}=1}^{{N}}({{r}({{h}_{{n}}})}-{b})
  \sum_{{t}=1}^{{T}}{\nabla_{{\boldsymbol{\theta}}}}\log{p({{a}^{{n}}_{{t}}}|{{\boldsymbol{s}}^{{n}}_{{t}}},{\boldsymbol{\theta}})},\end{aligned}$$ where $$\begin{aligned}
  {b}&=\frac{\frac{1}{{N}}\sum_{{n}=1}^{{N}}{{{r}({{h}_{{n}}})}
      {{\bigl\|{\nabla_{{\boldsymbol{\theta}}}}\log{p({{h}_{{n}}}|{\boldsymbol{\theta}})}\bigr\|}^2}}}
  {\frac{1}{{N}}\sum_{{n}=1}^{{N}}{{{\bigl\|{\nabla_{{\boldsymbol{\theta}}}}\log{p({{h}_{{n}}}|{\boldsymbol{\theta}})}\bigr\|}^2}}}.\end{aligned}$$

Let us consider the following Gaussian policy model with policy parameter ${\boldsymbol{\theta}}=({\boldsymbol{\mu}},{\sigma}^2)$: $$\begin{aligned}
  {p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}=\frac{1}{\sqrt{2\pi{\sigma}^2}}
    \exp\left(-\frac{({a}-{\boldsymbol{\mu}}^\top{\boldsymbol{\phi}}({\boldsymbol{s}}))^2}{2{\sigma}^2}\right),\end{aligned}$$ where $^\top$ denotes the transpose, ${\boldsymbol{\mu}}$ is the Gaussian mean, ${\sigma}$ is the Gaussian standard deviation, and ${\boldsymbol{\phi}}({\boldsymbol{s}})$ is the basis function vector. Then the policy gradients are explicitly expressed as $$\begin{aligned}
  {\nabla_{{\boldsymbol{\mu}}}}\log{p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}&=\frac{{a}-{\boldsymbol{\mu}}^\top{\boldsymbol{\phi}}({\boldsymbol{s}})}{{\sigma}^2}{\boldsymbol{\phi}}({\boldsymbol{s}}),\\
  {\nabla_{{\sigma}}}\log{p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}&=\frac{(a-{\boldsymbol{\mu}}^\top{\boldsymbol{\phi}}({\boldsymbol{s}}))^2-{\sigma}^2}{{\sigma}^3}.\end{aligned}$$

REINFORCE is a simple policy-search algorithm that directly updates policies to increase the expected return. However, gradient estimates tend to have large variance even if it is combined with variance reduction by baseline subtraction. For this reason, policy update by REINFORCE tends to be unreliable [@IROS:Peters+Schaal:2006]. In particular, the variance of gradient estimates in REINFORCE can be proportional to the length of the history, ${T}$, due to the stochasticity of policies [@NN:Zhao+etal:2012]. This can be a critical limitation when the history is long.

Policy Gradients with Parameter-Based Exploration (PGPE) {#subsec:PGPE}
--------------------------------------------------------

To overcome the above limitation of REINFORCE, a novel policy-search method called *policy gradients with parameter-based exploration* (PGPE) was proposed recently [@NN:Sehnke+etal:2010]. In PGPE, a deterministic policy (such as the linear policy) is adopted, and the stochasticity for exploration is introduced by drawing the policy parameter ${\boldsymbol{\theta}}$ from a prior distribution ${p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}$ with hyper-parameter ${\boldsymbol{\rho}}$. Thanks to this per-trajectory formulation, the variance of gradient estimates can be drastically reduced.

In the PGPE formulation, the expected return is represented as a function of the hyper-parameter ${\boldsymbol{\rho}}$: $$\begin{aligned}
  {J({\boldsymbol{\rho}})}&:=\iint{{r}({h})}{p({h}|{\boldsymbol{\theta}})}{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}\mathrm{d}{h}\mathrm{d}{\boldsymbol{\theta}}.\end{aligned}$$ Differentiating this with respect to ${\boldsymbol{\rho}}$, we have $$\begin{aligned}
  {\nabla_{{\boldsymbol{\rho}}}}{J({\boldsymbol{\rho}})}&=\iint{{r}({h})}{p({h}|{\boldsymbol{\theta}})}{\nabla_{{\boldsymbol{\rho}}}}{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}
  \mathrm{d}{h}\mathrm{d}{\boldsymbol{\theta}}\\
  &=\iint{{r}(h)}{p({h}|{\boldsymbol{\theta}})}{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}{\nabla_{{\boldsymbol{\rho}}}}\log{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}
  \mathrm{d}{h}\mathrm{d}{\boldsymbol{\theta}}.\end{aligned}$$

Because of the per-trajectory formulation, roll-out samples in the PGPE framework are accompanied with policy parameters, i.e., $\{({{h}_{{n}}},{{\boldsymbol{\theta}}_{{n}}})\}_{{n}=1}^{{N}}$. Based on these paired samples, an empirical estimator of the above gradient (with baseline subtraction) is given as follows [@NN:Zhao+etal:2012]: $$\begin{aligned}
  {\nabla_{{\boldsymbol{\rho}}}}{\widehat{J}^{{b}}({\boldsymbol{\rho}})}:=\frac{1}{{N}}\sum_{{n}=1}^{{N}}
  ({{r}({{h}_{{n}}})}-{b}){\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}|{\boldsymbol{\rho}})},\end{aligned}$$ where $$\begin{aligned}
  {b}&=\frac{\frac{1}{{N}}\sum_{{n}=1}^{{N}}{{{r}({{h}_{{n}}})}
      {{\bigl\|{\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}|{\boldsymbol{\rho}})}\bigr\|}^2}}}
  {\frac{1}{{N}}\sum_{{n}=1}^{{N}}{
      {{\bigl\|{\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}|{\boldsymbol{\rho}})}\bigr\|}^2}}}.\end{aligned}$$

Let us employ the linear deterministic policy, i.e., action ${a}$ is chosen as ${\boldsymbol{\theta}}^\top{\boldsymbol{\phi}}({\boldsymbol{s}})$ for some basis function ${\boldsymbol{\phi}}$. The parameter vector ${\boldsymbol{\theta}}$ is drawn from the Gaussian prior distribution with hyper-parameter ${\boldsymbol{\rho}}=({\boldsymbol{\eta}},{\boldsymbol{\tau}})$. Here ${\boldsymbol{\eta}}$ denotes the Gaussian mean vector and ${\boldsymbol{\tau}}$ denotes the vector consisting of the Gaussian standard deviation in each element: $$\begin{aligned}
  {p(\theta_i|\rho_i)}=\frac{1}{\sqrt{2\pi\tau_i^2}}
  \exp\left(-\frac{(\theta_i-\eta_i)^2}{2\tau_i^2}\right),\end{aligned}$$ where $\theta_i$, $\rho_i$, $\eta_i$, and $\tau_i$ are the $i$-th elements of ${\boldsymbol{\theta}}$, ${\boldsymbol{\rho}}$, ${\boldsymbol{\eta}}$, and ${\boldsymbol{\tau}}$, respectively. Then the derivatives of $\log {p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}$ with respect to $\eta_i$ and $\tau_i$ are given as follows: $$\begin{aligned}
  {\nabla_{\eta_i}}\log {p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}&=\frac{\theta_i-\eta_i}{\tau_i^2},\\
  {\nabla_{\tau_i}} \log {p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}&=\frac{(\theta_i-\eta_i)^2-\tau_i^2}{\tau_i^3}.\end{aligned}$$

Importance-Weighted PGPE (IW-PGPE) {#subsec:IWPGPE}
----------------------------------

A popular idea to further improve the performance of RL methods is to reuse previously collected samples [@book:Sutton+Barto:1998; @NN:Hachiya+etal:2009]. Such a sample-reuse strategy is particularly useful when data sampling costs is high (e.g., robot control).

*Importance-weighted* PGPE (IW-PGPE) [@NC:Zhao+etal:2013] combines the sample-reuse idea with PGPE. Technically, IW-PGPE can be regarded as an off-policy extension of PGPE, where data collecting policies are different from the current policy. In the PGPE formulation, such a off-policy scenario can be regarded as the situation where data collecting policies and the current policy are drawn from different prior distributions (more specifically, different hyper-parameters). Let ${\boldsymbol{\rho}}$ be the hyper-parameter for the current policy and ${\boldsymbol{\rho}}'$ be the hyper-parameter for a data collecting policy. Let us denote data samples collected with hyper-parameter ${\boldsymbol{\rho}}'$ as $\{({{\boldsymbol{\theta}}_{{n}}}',{{h}_{{n}}}')\}_{{n}=1}^{{N}'}$.

When the data collecting policy is different from the current policy, *importance sampling* is a useful technique to correct the estimation bias caused by differing distributions [@book:Sugiyama+Kawanabe:2012]. More specifically, the gradient is estimated as $$\begin{aligned}
{\nabla_{{\boldsymbol{\rho}}}}{\widehat{J}_{{b}}^{w}({\boldsymbol{\rho}})}{}=
\frac{1}{{N}'}\sum_{{n}=1}^{{N}'}w({{\boldsymbol{\theta}}_{{n}}}')({{r}({{h}_{{n}}}')}-{b})
{\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}'|{\boldsymbol{\rho}})},\end{aligned}$$ where $w({\boldsymbol{\theta}})$ is the *importance weight* defined as $$\begin{aligned}
w({\boldsymbol{\theta}}):=\frac{{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}}{{p({\boldsymbol{\theta}}|{\boldsymbol{\rho}}')}},\end{aligned}$$ and ${b}$ is the baseline given by $$\begin{aligned}
{b}&=\frac{\frac{1}{{N}'}\sum_{{n}=1}^{{N}'}
{{r}({{h}_{{n}}}')}w^2({{\boldsymbol{\theta}}_{{n}}}'){{\bigl\|{\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}'|{\boldsymbol{\rho}})}\bigr\|}^2}}
{\frac{1}{{N}'}\sum_{{n}=1}^{{N}'}
w^2({{\boldsymbol{\theta}}_{{n}}}'){{\bigl\|{\nabla_{{\boldsymbol{\rho}}}}\log{p({{\boldsymbol{\theta}}_{{n}}}'|{\boldsymbol{\rho}})}\bigr\|}^2}}.\end{aligned}$$ Through experiments, the IW-PGPE method was demonstrated to be the best performing algorithm in model-free RL approaches [@NC:Zhao+etal:2013].

The purpose of this paper is to develop a model-based counterpart of PGPE.

Model-Based Policy Search {#sec:MB}
=========================

Model-based RL first estimates the transition model and then learns a policy based on the estimated transition model. Because one can draw as many trajectory samples as one wants from the learned transition model without additional sampling costs, the model-based approach can work well if the transition model is accurately estimated [@Wang:2003; @Deisenroth:2011]. In this section, we extend PGPE to a model-based scenario. We first review an existing model estimation method based on the Gaussian process (GP) [@book:Rasmussen+Williams:2006] and point out its limitations. Then we propose to use the state-of-the-art conditional density estimator called *least-squares conditional density estimation* (LSCDE) [@IEICE:Sugiyama+etal:2010a] in the model-based PGPE method.

Model-Based PGPE (M-PGPE) {#subsec:prop}
-------------------------

PGPE can be extended to a model-based scenario as follows.

1.  Collect transition samples $\{({\boldsymbol{s}}_m,{a}_m,{\boldsymbol{s}}'_m)\}_{m=1}^M$.

2.  Obtain transition model ${\widehat{p}({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ by a model estimation method from $\{({\boldsymbol{s}}_m,{a}_m,{\boldsymbol{s}}'_m)\}_{m=1}^M$.

3.  Initialize hyper-parameter ${\boldsymbol{\rho}}$.

4.  Draw policy parameter ${\boldsymbol{\theta}}$ from prior distribution ${p({\boldsymbol{\theta}}|{\boldsymbol{\rho}})}$.

5.  Generate many samples $\{\widetilde{h}_{n}\}_{n=1}^{\widetilde{N}}$ from ${\widehat{p}({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ and current policy ${p({a}|{\boldsymbol{s}},{\boldsymbol{\theta}})}$.

6.  Estimate baseline ${b}$ and gradient ${\nabla_{{\boldsymbol{\rho}}}}{\widehat{J}^{{b}}({\boldsymbol{\rho}})}$ from disjoint subsets of $\{\widetilde{h}_{n}\}_{n=1}^{\widetilde{N}}$.

7.  Update hyper-parameter as ${\boldsymbol{\rho}}\leftarrow{\boldsymbol{\rho}}+{\varepsilon}{\nabla_{{\boldsymbol{\rho}}}}{\widehat{J}^{{b}}({\boldsymbol{\rho}})}$, where ${\varepsilon}>0$ denotes the learning rate.

8.  Repeat Steps 4–7 until ${\boldsymbol{\rho}}$ converges.

Below, we consider the problem of approximating the transition probability ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ from samples $\{({\boldsymbol{s}}_m,{a}_m,{\boldsymbol{s}}'_m)\}_{m=1}^M$, and review transition model estimation methods.

Gaussian Process (GP)
---------------------

Here we review a transition model estimation method based on GP.

In the GP framework, the problem of transition probability estimation is formulated as the regression problem of predicting output ${\boldsymbol{s}}'$ given input ${\boldsymbol{s}}$ and ${a}$ under Gaussian noise: $${\boldsymbol{s}}' = f({\boldsymbol{s}}, {a}) + \varepsilon,$$ where $f$ is an unknown regression function and $\varepsilon\sim \mathcal{N}(0,\sigma_{\varepsilon}^2)$ is independent Gaussian noise. Then, the GP estimate of the transition probability density ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ for an arbitrary test input ${\boldsymbol{s}}$ and ${a}$ is given by the Gaussian distribution with mean and variance given by $$\begin{aligned}
\boldsymbol{k}^\top(\boldsymbol{K}+\sigma_{\varepsilon}^2 \boldsymbol{I}_M)^{-1}\boldsymbol{y}
~~~\mbox{and}~~~
k({\boldsymbol{s}},{a},{\boldsymbol{s}},{a})-\boldsymbol{k}^\top(\boldsymbol{K}+\sigma_{\varepsilon}^2 \boldsymbol{I}_M)^{-1}\boldsymbol{k},\end{aligned}$$ respectively. Here, ${\boldsymbol{I}_{M}}$ denotes the $M$-dimensional identity matrix. $\boldsymbol{k}$ is the $M$-dimensional vector and $\boldsymbol{K}$ is the $M\times M$ Gram matrix defined by $$\begin{aligned}
  \boldsymbol{k}&=
  \begin{pmatrix}
    k({\boldsymbol{s}}_1,{a}_1,{\boldsymbol{s}},{a})\\
    \vdots\\
    k({\boldsymbol{s}}_M,{a}_M,{\boldsymbol{s}},{a})
  \end{pmatrix}
  ~~~\mbox{and}~~~
  \boldsymbol{K}=
  \begin{pmatrix}
 k({\boldsymbol{s}}_1,{a}_1,{\boldsymbol{s}}_1,{a}_1)&\ldots&k({\boldsymbol{s}}_M,{a}_M,{\boldsymbol{s}}_1,{a}_1)\\
 \vdots&\ddots&\vdots\\
 k({\boldsymbol{s}}_1,{a}_1,{\boldsymbol{s}}_M,{a}_M)&\ldots&k({\boldsymbol{s}}_M,{a}_M,{\boldsymbol{s}}_M,{a}_M)\\
  \end{pmatrix}.\end{aligned}$$ $k({\boldsymbol{s}},{a},{\boldsymbol{s}}',{a}')$ denotes the covariance function, which is, e.g., defined by $$k({\boldsymbol{s}},{a},{\boldsymbol{s}}',{a}') = \theta \exp\left(-([{\boldsymbol{s}}^\top,{a}]-[{\boldsymbol{s}}'^\top,{a}'])\boldsymbol{\Theta}([{\boldsymbol{s}}^\top,{a}]-[{\boldsymbol{s}}'^\top,{a}'])^\top\right).$$ Here, $\theta$ and $\boldsymbol{\Theta}$ are hyperparameters, and together with the noise variance $\sigma_{\varepsilon}^2$, the hyperparameters are determined by evidence maximization [@book:Rasmussen+Williams:2006].

As shown above, the GP-based model estimation method requires the strong assumption that the transition probability density ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ is Gaussian. That is, GP is non-parametric as a regression method of estimating the conditional mean, it is parametric (Gaussian) as a conditional density estimator. Such a conditional Gaussian assumption is highly restrictive in RL problems.

Least-Squares Conditional Density Estimation (LSCDE) {#subsec:LSCDE}
----------------------------------------------------

To overcome the restriction of the GP-based model estimation method, we propose to use LSCDE.

Let us model the transition probability ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}$ by the following linear-in-parameter model: $$\begin{aligned}
{q({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}&:={{\boldsymbol{\alpha}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')\nonumber\\
&\phantom{:}=
\sum_{m=1}^M
\alpha_m\exp\left(-\frac{\|{\boldsymbol{s}}-{\boldsymbol{s}}_m\|^2}{2\kappa^2}\right)
\exp\left(-\frac{({a}-{a}_m)^2}{2\kappa^2}\right)
\exp\left(-\frac{\|{\boldsymbol{s}}'-{\boldsymbol{s}}'_m\|^2}{2\kappa^2}\right),
\label{LSCDE-model}\end{aligned}$$ where ${\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')$ is the $M$-dimensional basis function vector and ${\boldsymbol{\alpha}}$ is the $M$-dimensional parameter vector. If $M$ is too large, we may reduce the number of basis functions by only using a subset of samples as Gaussian centers. We may use different Gaussian widths for ${\boldsymbol{s}}$ and ${a}$ if necessary.

The parameter ${\boldsymbol{\alpha}}$ in the model is learned so that the following squared error is minimized: $$\begin{aligned}
{G_0({\boldsymbol{\alpha}})}:=\frac{1}{2}\iiint\Big({q({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}-{p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}\Big)^2
p({\boldsymbol{s}},{a})\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'.\end{aligned}$$ This can be expressed as $$\begin{aligned}
{G_0({\boldsymbol{\alpha}})}
&=\frac{1}{2}\iiint{q({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}^2p({\boldsymbol{s}},{a})\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'\\
&\phantom{=}
-\iiint{q({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}p({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'+C\\
&=\frac{1}{2}\iiint\Big({{\boldsymbol{\alpha}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')\Big)^2
p({\boldsymbol{s}},{a})\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'\\
&\phantom{=}
-\iiint{{\boldsymbol{\alpha}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')p({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'+C,\end{aligned}$$ where we used ${p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}=p({\boldsymbol{s}},{a},{\boldsymbol{s}}')/p({\boldsymbol{s}},{a})$ in the second term and $$\begin{aligned}
C:=\frac{1}{2}\iiint{p({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}p({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'.\end{aligned}$$ Because $C$ is constant, we only consider the first two terms from here on: $$\begin{aligned}
{G({\boldsymbol{\alpha}})}&:={G_0({\boldsymbol{\alpha}})}-C\\
&\phantom{:}=\frac{1}{2}{{\boldsymbol{\alpha}}^\top}{\boldsymbol{H}}{\boldsymbol{\alpha}}-{{\boldsymbol{\alpha}}^\top}{\boldsymbol{h}},\end{aligned}$$ where $$\begin{aligned}
{\boldsymbol{H}}&:=\iint{\overline{\boldsymbol{\Phi}}}({\boldsymbol{s}},{a})p({\boldsymbol{s}},{a})\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\in\mathbb{R}^{M\times M},\\
{\boldsymbol{h}}&:=\iiint{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')p({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}\mathrm{d}{a}\mathrm{d}{\boldsymbol{s}}'
\in\mathbb{R}^{M},\\
{\overline{\boldsymbol{\Phi}}}({\boldsymbol{s}},{a})&:=\int{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}'){{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')^\top}\mathrm{d}{\boldsymbol{s}}'
\in\mathbb{R}^{M\times M}.\end{aligned}$$ Note that, for the Gaussian model , the $(m,m')$-th element of ${\overline{\boldsymbol{\Phi}}}({\boldsymbol{s}},{a})$ can be computed analytically as $$\begin{aligned}
{\overline{\boldsymbol{\Phi}}}_{m,m'}({\boldsymbol{s}},{a})
&=(\sqrt{\pi}\kappa)^{\mathrm{dim}({\boldsymbol{s}}')}
\exp\left(-\frac{\|{\boldsymbol{s}}-{\boldsymbol{s}}_m\|^2+\|{\boldsymbol{s}}-{\boldsymbol{s}}_{m'}\|^2}{2\kappa^2}\right)\\
&\phantom{=}\times
\exp\left(-\frac{({a}-{a}_m)^2+({a}-{a}_{m'})^2}{2\kappa^2}\right)
\exp\left(-\frac{\|{\boldsymbol{s}}'_m-{\boldsymbol{s}}'_{m'}\|^2}{4\kappa^2}\right).\end{aligned}$$

Because ${\boldsymbol{H}}$ and ${\boldsymbol{h}}$ included in ${G({\boldsymbol{\alpha}})}$ contain the expectations over unknown densities $p({\boldsymbol{s}},{a})$ and $p({\boldsymbol{s}},{a},{\boldsymbol{s}}')$, they are approximated by sample averages. Then we have $$\begin{aligned}
{\widehat{G}({\boldsymbol{\alpha}})}:=\frac{1}{2}{{\boldsymbol{\alpha}}^\top}{\widehat{{\boldsymbol{H}}}}{\boldsymbol{\alpha}}-{{\widehat{{\boldsymbol{h}}}}^\top}{\boldsymbol{\alpha}},\end{aligned}$$ where $$\begin{aligned}
{\widehat{{\boldsymbol{H}}}}&:=\frac{1}{M}\sum^{M}_{m=1}{\overline{\boldsymbol{\Phi}}}({\boldsymbol{s}}_m,{a}_m)\in\mathbb{R}^{M\times M},\\
{\widehat{{\boldsymbol{h}}}}&:=\frac{1}{M}\sum^{M}_{m=1}{\boldsymbol{\phi}}({\boldsymbol{s}}_m,{a}_m,{\boldsymbol{s}}'_m)\in\mathbb{R}^{M}.\end{aligned}$$

By adding an $\ell_2$-regularizer to ${\widehat{G}({\boldsymbol{\alpha}})}$ to avoid overfitting, the LSCDE optimization criterion is given as $$\begin{aligned}
{\widetilde{{\boldsymbol{\alpha}}}}:={\mathop{\mathrm{argmin}}}_{{\boldsymbol{\alpha}}\in\mathbb{R}^M}\left[{\widehat{G}({\boldsymbol{\alpha}})}+\frac{\lambda}{2}\|{\boldsymbol{\alpha}}\|^2\right],\end{aligned}$$ where $\lambda$ $(\ge0)$ is the regularization parameter. Taking the derivative of the above objective function and equating it to zero, we can see that the solution ${\widetilde{{\boldsymbol{\alpha}}}}$ can be obtained just by solving the following system of linear equations: $$\begin{aligned}
({\widehat{{\boldsymbol{H}}}}+\lambda{\boldsymbol{I}_{M}}){\boldsymbol{\alpha}}={\widehat{{\boldsymbol{h}}}},\end{aligned}$$ where ${\boldsymbol{I}_{M}}$ denotes the $M$-dimensional identity matrix. Thus, the solution ${\widetilde{{\boldsymbol{\alpha}}}}$ is given analytically as $$\begin{aligned}
{\widetilde{{\boldsymbol{\alpha}}}}=({\widehat{{\boldsymbol{H}}}}+\lambda{\boldsymbol{I}_{M}})^{-1}{\widehat{{\boldsymbol{h}}}}.\end{aligned}$$ Because conditional probability densities are non-negative by definition, we modify the solution ${\widetilde{{\boldsymbol{\alpha}}}}$ as $$\begin{aligned}
{\widehat{{\boldsymbol{\alpha}}}}:={\mathop{\mathrm{max}}}(\boldsymbol{0}_M,{\widetilde{{\boldsymbol{\alpha}}}}),\end{aligned}$$ where $\boldsymbol{0}_M$ denotes the $M$-dimensional zero vector and ‘${\mathop{\mathrm{max}}}$’ for vectors are applied in the element-wise manner.

Finally, we renormalize the solution in the test phase. More specifically, given a test input point $({\boldsymbol{s}},{a})$, the final LSCDE solution is given as $$\begin{aligned}
{\widehat{p}({\boldsymbol{s}}'|{\boldsymbol{s}},{a})}=\frac{{{\widehat{{\boldsymbol{\alpha}}}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')}
{\int{{\widehat{{\boldsymbol{\alpha}}}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}'},
\label{eq:est-p}\end{aligned}$$ where, for the Gaussian model , the denominator in Eq. can be analytically computed as $$\begin{aligned}
  \int{{\widehat{{\boldsymbol{\alpha}}}}^\top}{\boldsymbol{\phi}}({\boldsymbol{s}},{a},{\boldsymbol{s}}')\mathrm{d}{\boldsymbol{s}}'
=(\sqrt{2\pi}\kappa)^{\mathrm{dim}({\boldsymbol{s}}')}
\sum_{m=1}^M\alpha_m
\exp\left(-\frac{\|{\boldsymbol{s}}-{\boldsymbol{s}}_m\|^2+({a}-{a}_m)^2}{2\kappa^2}\right).\end{aligned}$$ LSCDE was proved to achieve the optimal non-parametric convergence rate to the true conditional density in the mini-max sense [@IEICE:Sugiyama+etal:2010a], meaning that no method can outperform this simple LSCDE method asymptotically.

Model selection of the Gaussian width $\kappa$ and the regularization parameter $\lambda$ is possible by cross-validation. A MATLAB$^\text{\textregistered}$ implementation of LSCDE is available from

‘<http://sugiyama-www.cs.titech.ac.jp/~sugi/software/LSCDE/>’.

Experiments {#sec:experiments}
===========

In this section, we demonstrate the usefulness of the proposed method through experiments.

Continuous Chain Walk {#subsec:chain}
---------------------

For illustration purposes, let us first consider a simple continuous chain-walk task (Figure \[fig:cwalk-view\]).

### Setup

Let $$\begin{aligned}
{{\mathcal S}}&=[0, 10],\\
{{\mathcal A}}&=[-5, 5],\\
{{R}(s,{a},s')}&=
\begin{cases}
1&(4<s'<6),\\
0&(\mathrm{otherwise}).\\
\end{cases}\end{aligned}$$ That is, the agent receives positive reward $+1$ at the center of the state space. We set the episode length at ${T}=10$, the discount factor at ${\gamma}=0.99$, and the learning rate at ${\varepsilon}=0.1$. We use the following linear-in-parameter policy model: $$\begin{aligned}
\sum_{i=1}^6\theta_i\exp\left(-\frac{(s-c_i)^2}{2}\right),\end{aligned}$$ where $(c_1,\ldots,c_6)=(0,2,4,6,8,10)$.

![Illustration of continuous chain walk.[]{data-label="fig:cwalk-view"}](graph/cwalk-view.eps){width="60.00000%"}

As transition dynamics, we consider two setups:

Gaussian:

:   The true transition dynamics is given by $$\begin{aligned}
    s_{t+1}=s_{t}+{a}_{t}+\epsilon_t,\end{aligned}$$ where $\epsilon_t$ is the Gaussian noise with mean $0$ and standard deviation $0.3$.

Bimodal:

:   The true transition dynamics is given by $$\begin{aligned}
    s_{t+1}=s_{t}\pm{a}_{t}+\epsilon_t,\end{aligned}$$ where $\epsilon_t$ is the Gaussian noise with mean $0$ and standard deviation $0.3$, and the sign of $a_{t}$ is randomly chosen with probability $1/2$.

We compare the following three policy search methods:

M-PGPE(LSCDE):

:   The model-based PGPE method with transition model estimated by LSCDE.

M-PGPE(GP):

:   The model-based PGPE method with transition model estimated by GP.

IW-PGPE:

:   The model-free PGPE method with sample reuse by importance weighting[^1] [@NC:Zhao+etal:2013].

Below, we consider the situation where the budget for data collection is limited to ${N}=20$ episodic samples.

### LSCDE Vs. GP {#subsubsec:model}

When the transition model is learned in the M-PGPE methods, all ${N}=20$ samples are gathered randomly in the beginning at once. More specifically, the initial state $s_1$ and the action ${a}_1$ are chosen from the uniform distributions over ${{\mathcal S}}$ and ${{\mathcal A}}$, respectively. Then the next state $s_2$ and the immediate reward $r_1$ are obtained. Then the action ${a}_2$ is chosen from the uniform distribution over ${{\mathcal A}}$, and the next state $s_3$ and the immediate reward $r_2$ are obtained. This process is repeated until we obtain $r_{{T}}$. This gives a trajectory sample, and we repeat this data generation process ${N}$ times to obtain ${N}$ trajectory samples.

Figure \[fig:trans\_pure\] and Figure \[fig:trans\_bi\] illustrate the true transition dynamics and its estimates obtained by LSCDE and GP in the Gaussian and bimodal cases. Figure \[fig:trans\_pure\] shows that both LSCDE and GP can learn the entire profile of the true transition dynamics well in the Gaussian case. On the other hand, Figure \[fig:trans\_bi\] shows that LSCDE can still successfully capture the entire profile of the true transition model well even in the bimodal case, but GP fails to capture the bimodal structure.

Based on the estimated transition models, we learn policies by the M-PGPE method. We generate $1000$ artificial samples for policy gradient estimation and another $1000$ artificial samples for baseline estimation from the learned transition model. Then policy is updated based on these artificial samples. We repeat this policy update step $20$ times. For evaluating the return of a learned policy, we use $100$ additional test episodic samples which are not used for policy learning. Figure \[fig:cwalk-results-pure\] and Figure \[fig:cwalk-results-bi\] depict the average performance of learned policies over $100$ runs. As expected, the GP-based method performs very well in the Gaussian case, but LSCDE still exhibits reasonably good performance. In the bimodal case, GP performs poorly and LSCDE gives much better policies than GP. This illustrates the high flexibility of LSCDE.

### Model-Based Vs. Model-Free

Next, we compare the performance of M-PGPE with the model-free IW-PGPE method.

For the IW-PGPE method, we need to determine the schedule of collecting $20$ samples under the fixed budget scenario. First, we illustrate how the choice of sampling schedules affects the performance of IW-PGPE. Figure \[fig:cwalk-schedule-pure\] and Figure \[fig:cwalk-schedule-bi\] show expected returns averaged over $100$ runs under the sampling schedule that a batch of $k$ samples are gathered $20/k$ times for different $k$ values. In our implementation of IW-PGPE, policy update is performed $100$ times after observing each batch of $k$ samples, because we empirically observed that this performs better than performing policy update only once. Figure \[fig:cwalk-schedule-pure\] shows that the performance of IW-PGPE depends heavily on the sampling schedule, and gathering $k=20$ samples at once is shown to be the best choice in the Gaussian case. Figure \[fig:cwalk-schedule-bi\] shows that gathering $k=20$ samples at once is also the best choice in the bimodal case.

Although the best sampling schedule is not accessible in practice, we use this optimal sampling schedule for IW-PGPE. Figure \[fig:cwalk-results-pure\] and Figure \[fig:cwalk-results-bi\] also include returns of IW-PGPE averaged over $100$ runs as functions of the sampling steps. These graphs show that IW-PGPE can improve the policies only in the beginning, because all samples are gathered at once in the beginning. The performance of IW-PGPE may be further improved if it is possible to gather more samples, but this is prohibited under the fixed budget scenario. On the other hand, return values of M-PGPE constantly increase throughout iterations, because artificial samples can be kept generated without additional sampling costs. This illustrates a potential advantage of model-based RL methods.

\

![Returns obtained by IW-PGPE averaged over $100$ runs for Gaussian transition dynamics with different sampling schedules (e.g., $5\times4$ means gathering $k=5$ samples $4$ times).[]{data-label="fig:cwalk-schedule-pure"}](graph/performance_pure.eps){width="100.00000%"}

  

![Returns obtained by IW-PGPE averaged over $100$ runs for Gaussian transition dynamics with different sampling schedules (e.g., $5\times4$ means gathering $k=5$ samples $4$ times).[]{data-label="fig:cwalk-schedule-pure"}](graph/sampling_pure.eps){width="100.00000%"}

\

![Returns obtained by IW-PGPE averaged over $100$ runs for bimodal transition with different sampling schedules (e.g., $5\times4$ means gathering $k=5$ samples $4$ times).[]{data-label="fig:cwalk-schedule-bi"}](graph/performance_bi.eps){width="\textwidth"}

  

![Returns obtained by IW-PGPE averaged over $100$ runs for bimodal transition with different sampling schedules (e.g., $5\times4$ means gathering $k=5$ samples $4$ times).[]{data-label="fig:cwalk-schedule-bi"}](graph/sampling_cw_ting.eps){width="\textwidth"}

Humanoid Robot Control {#subsec:humanoid}
----------------------

Finally, we evaluate the performance of M-PGPE on a practical control problem of a simulated upper-body model of the humanoid robot *CB-i* [@Cheng:2007] (see Figure \[fig:cbi\]). We use its simulator for experiments (see Figure \[fig:simulator\]). The goal of the control problem is to lead the end-effector of the right arm (right hand) to the target object.

### Setup

The simulator is based on the upper-body of the CB-i humanoid robot, which has 9 joints for shoulder pitch, shoulder roll, elbow pitch of the right arm, shoulder pitch, shoulder roll, elbow pitch of the left arm, waist yaw, torso roll, and torso pitch.

At each time step, the controller receives a state vector from the system and sends out an action vector. The state vector is 18-dimensional and real-valued, which corresponds to the current angle in degree and the current angular velocity for each joint. The action vector is 9-dimensional and real-valued, which corresponds to the target angle of each joint in degree.

We simulate a noisy control system by perturbing action vectors with independent bimodal Gaussian noise. More specifically, for each action element, we add Gaussian noise with mean $0$ and standard deviation $3$ with probability $0.6$, and Gaussian noise with mean $-5$ and standard deviation $3$ with probability $0.4$.

The initial posture of the robot is fixed to standing up straight with arms down. The the target object is located in front-above of the right hand which is reachable by using the controllable joints. The reward function at each time step is defined as $$r_t=\exp(-10 d_t) - 0.000005 \min\{c_t,1000000\},$$ where $d_t$ is the distance between the right hand and target object at time step $t$, and $c_t$ is the sum of control costs for each joint. The coefficient $0.000005$ is multiplied to keep the values of the two terms in the same order of magnitude. The deterministic policy model used in PGPE is defined as $a=\boldsymbol{\theta}^\top \boldsymbol{\phi(s)}$ with the basis function $\boldsymbol{\phi(s)} = \boldsymbol{s}$. We set the episode length at $T = 100$, the discount factor at $\gamma = 0.9$, and the learning rate at ${\varepsilon}=0.1/\|\nabla_{\rho} \hat{J}(\bm{\rho})\|$.

### Experiment with 2 Joints

First, we only use 2 joints among the 9 joints, i.e., we allow only the right shoulder pitch and right elbow pitch to be controlled, while the other joints remain still at each time step (no control signal is sent to these joints). Therefore, the dimensionality of state vector $\boldsymbol{s}$ and action vector $\boldsymbol{a}$ is $4$ and $2$, respectively. Under this simplified setup, we compare the performance of M-PGPE(LSCDE), M-PGPE(GP), and IW-PGPE.

We suppose that the budget for data collection is limited to $N = 50$ episodic samples. For the M-PGPE methods, all samples are collected at first using the uniformly random initial states and policy. More specifically, the initial state is chosen from the uniform distributions over $\mathcal{S}$. At each time step, the $i$-th element of action vector $a_i$ is chosen from the uniform distribution on $[s_i-5,s_i+5]$. In total, we have 5000 transition samples for model estimation. Then, we generate 1000 artificial samples for policy gradient estimation and another 1000 artificial samples for baseline estimation from the learned transition model, and update the control policy based on these artificial samples. For the IW-PGPE method, we performed preliminary experiments to determine the optimal sampling schedule (Figure \[fig:simu\_schedule2j\]), showing that collecting $k=5$ samples $50/k$ times yields the highest average return. We use this sampling schedule for performance comparison with the M-PGPE methods. Returns obtained by each method averaged over 10 runs are plotted in Figure \[fig:simu\_result2j\], showing that M-PGPE(LSCDE) tends to outperform both M-PGPE(GP) and IW-PGPE. Figure \[fig:simu\_traject2j\] illustrates an example of the reaching motion with 2-joints obtained by M-PGPE(LSCDE) at the $60$th iteration policy. This shows that the learned policy successfully leads the right hand to the target object within only $13$ steps in this noisy control system.

![ Averages and standard errors of obtained returns over 10 runs in 2-joint humanoid robot simulator. All methods use 50 samples for policy learning. In M-PGPE(LSCDE) and M-PGPE(GP), all 50 samples are gathered in the beginning and the environment model is learned; then 2000 artificial samples are generated in each update iteration. In IW-PGPE, a batch of 5 samples are gathered for 10 iterations, which was shown to be the best sampling scheduling (see Figure \[fig:simu\_schedule2j\]). Note that policy update is performed 100 times after observing each batch of samples, which we confirmed to perform well. The IW-PGPE curve is elongated to have the same horizontal scale as others. []{data-label="fig:simu_result2j"}](graph/simu_schedule2j.eps){width="60.00000%"}

![ Averages and standard errors of obtained returns over 10 runs in 2-joint humanoid robot simulator. All methods use 50 samples for policy learning. In M-PGPE(LSCDE) and M-PGPE(GP), all 50 samples are gathered in the beginning and the environment model is learned; then 2000 artificial samples are generated in each update iteration. In IW-PGPE, a batch of 5 samples are gathered for 10 iterations, which was shown to be the best sampling scheduling (see Figure \[fig:simu\_schedule2j\]). Note that policy update is performed 100 times after observing each batch of samples, which we confirmed to perform well. The IW-PGPE curve is elongated to have the same horizontal scale as others. []{data-label="fig:simu_result2j"}](graph/simu_result2j.eps){width="60.00000%"}

### Experiment with 9 Joints

Finally, we evaluate the performance of the proposed method on the reaching task with 9 joints, i.e., all joints are allowed to move. In this experiment, we compare learning performance between M-PGPE(LSCDE) and IW-PGPE. We do not include M-PGPE(GP) since it is outperformed by M-PGPE(LSCDE) on the previous 2-joints experiments, and furthermore the GP-based method requires an enormous amount of computation time. The experimental setup is essentially the same as the 2-joints experiments, but we have a budget for gathering $N = 1000$ samples for this complex and high-dimensional task. The position of the target object is moved to far left, which is not reachable by using just 2-joints. Thus, the robot is required to move other joints to reach the object with right hand. We randomly choose 5000 samples for Gaussian centers for M-PGPE(LSCDE). The sampling schedule for IW-PGPE was set to 1000 samples at once, which is the best sampling schedule according to Figure \[fig:simu\_schedule9j\]. The returns obtained by M-PGPE(LSCDE) and IW-PGPE averaged over 30 runs are plotted in Figure \[fig:simu\_result9j\], showing that M-PGPE(LSCDE) tends to outperform the state-of-the-art IW-PGPE method in this challenging robot control task.

Figure \[fig:simu-traject9j\] shows a typical example of the reaching motion with 9 joints obtained by M-PGPE(LSCDE) at the 1000th iteration. The images show that the policy learned by M-PGPE(LSCDE) leads the right hand to the distant object successfully within 14 steps.

![ Returns obtained by IW-PGPE averaged over 30 runs in humanoid robot simulator with 9 joints for different sampling schedules (e.g., $100 \times 10$ means gathering $k=100$ samples 10 times). []{data-label="fig:simu_schedule9j"}](graph/simu_schedule9j.eps){width="60.00000%"}

![ Averages and standard errors of obtained returns over 30 runs in humanoid robot simulator with 9 joints. Both methods use 1000 samples for policy learning. In M-PGPE(LSCDE), all 1000 samples are gathered in the beginning and the environment model is learned; then 2000 artificial samples are generated in each update iteration. In IW-PGPE, a batch of 1000 samples are gathered at once, which was shown to be the best scheduling (see Figure \[fig:simu\_schedule9j\]). Note that policy update is performed 100 times after observing each batch of samples. The IW-PGPE curve is elongated to have the same horizontal scale as M-PGPE(LSCDE). []{data-label="fig:simu_result9j"}](graph/simu_result9j.eps){width="60.00000%"}

Overall, the proposed M-PGPE(LSCDE) method is shown to be promising in the noisy and high-dimensional humanoid robot arm reaching task.

Conclusion {#sec:conclusion}
==========

We extended the model-free PGPE method to a model-based scenario, and proposed to combine it with a model estimator called LSCDE. Under the fixed sampling budget, appropriately designing a sampling schedule is critical for the model-free IW-PGPE method, while this is not a problem for the proposed model-based PGPE method. Through experiments, we confirmed that GP-based model estimation is not as flexible as the LSCDE-based method when the transition model is not Gaussian, and the proposed model-based PGPE based on LSCDE was overall demonstrated to be promising.

Acknowledgments {#acknowledgments .unnumbered}
===============

VT was supported by the JASSO scholarship, TZ was supported by the MEXT scholarship, JM was supported by MEXT KAKENHI 23120004, and MS was supported by the FIRST project.

[10]{}

P. Abbeel, M. Quigley, and A. Y. Ng. Using inaccurate models in reinforcement learning. , pages 1–8, 2006.

T. Akiyama, H. Hachiya, and M. Sugiyama. Efficient exploration through active learning for value function approximation in reinforcement learning. , 23(5):639–648, 2010.

G. Cheng, S. Hyon, J. Morimoto, A. Ude, G.H. Joshua, Glenn Colvin, Wayco Scroggin, and C.J. Stephen. Cb: A humanoid research platform for exploring neuroscience. , 21(10):1097–1114, 2007.

P. Dayan and G. E. Hinton. Using expectation-maximization for reinforcement learning. , 9(2):271–278, 1997.

M. P. Deisenroth and C. E. Rasmussen. Pilco: A model-based and data-efficient approach to policy search. , pages 465–473, 2011.

H. Hachiya, T. Akiyama, M. Sugiyama, and J. Peters. Adaptive importance sampling for value function approximation in off-policy reinforcement learning. , 22(10):1399–1410, 2009.

L. P. Kaelbling, M. L. Littman, and A. W. Moore. Reinforcement learning: A survey. , 4:237–285, 1996.

S. Kakade. A natural policy gradient. In T. G. Dietterich, S. Becker, and Z. Ghahramani, editors, [ *Advances in Neural Information Processing Systems 14*]{}, pages 1531–1538, Cambridge, MA, 2002. MIT Press.

T. Kanamori, S. Hido, and M. Sugiyama. A least-squares approach to direct importance estimation. , 10:1391–1445, Jul. 2009.

T. Kanamori, T. Suzuki, and M. Sugiyama. Computational complexity of kernel-based density-ratio estimation: [A]{} condition number analysis. , 2012. to appear.

T. Kanamori, T. Suzuki, and M. Sugiyama. Statistical analysis of kernel-based least-squares density-ratio estimation. , 86(3):335–367, 2012.

M. G. Lagoudakis and R. Parr. Least-squares policy iteration. , 4:1107–1149, 2003.

J. Peters and S. Schaal. Policy gradient methods for robotics. In [*Processing of the IEEE/RSJ Internatinal Conference on Inatelligent Robots and Systems (IROS)*]{}, pages 2219–2225, 2006.

C. E. Rasmussen and C. K. I. Williams. . Press, Cambridge, MA, USA, 2006.

F. Sehnke, C. Osendorfer, T. Rückstiess, A. Graves, J. Peters, and J. Schmidhuber. Parameter-exploring policy gradients. , 23(4):551–559, 2010.

M. Sugiyama, H. Hachiya, H. Kashima, and T. Morimura. Least absolute policy iteration—[A]{} robust approach to value function approximation. , E93-D(9):2555–2565, 2010.

M. Sugiyama, H. Hachiya, C. Towell, and S. Vijayakumar. Geodesic [G]{}aussian kernels for value function approximation. , 25(3):287–304, 2008.

M. Sugiyama and M. Kawanabe. . MIT Press, Cambridge, Massachusetts, USA, 2012.

M. Sugiyama, T. Suzuki, and T. Kanamori. Density ratio matching under the [B]{}regman divergence: A unified framework of density ratio estimation. , 64(5):1009–1044, 2012.

M. Sugiyama, I. Takeuchi, T. Suzuki, T. Kanamori, H. Hachiya, and D. Okanohara. Least-squares conditional density estimation. , E93-D(3):583–594, 2010.

R. S. Sutton and G. A. Barto. . Press, Cambridge, MA, USA, 1998.

R. S. Sutton, D. McAllester, S. Singh, and Y. Mansour. Policy gradient methods for reinforcement learning with function approximation. In S.A. Solla, T.K. Leen, and K.-R. M[ü]{}ller, editors, [ *Advances in Neural Information Processing Systems 12*]{}, pages 1057–1063. MIT Press, 2000.

X. Wang and T. G. Dietterich. Model-based policy gradient reinforcement learning. , pages 776–783, 2003.

Ronald J. Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. , 8:229, 1992.

T. Zhao, H. Hachiya, G. Niu, and M. Sugiyama. Analysis and improvement of policy gradient estimation. , 26:118–129, 2012.

T. Zhao, H. Hachiya, V. Tangkaratt, J. Morimoto, and M. Sugiyama. Efficient sample reuse in policy gradients with parameter-based exploration. , 25:1512–1547, 2013.

[^1]: We have also tested the plain PGPE method without importance weighting. However, this did not perform well in our preliminary experiments, and thus we decided to omit the results.
