# On Effective Scheduling of Model-based Reinforcement Learning

Hang Lai 1, Jian Shen∗1, Weinan Zhang†1, Yimin Huang2,
Xing Zhang2, Ruiming Tang2, Yong Yu1, Zhenguo Li2
1Shanghai Jiao Tong University, 2Huawei Noah’s Ark Lab
{laihang<EMAIL_ADDRESS>
equal contribution. †Weinan Zhang is the corresponding author.

###### Abstract

Model-based reinforcement learning has attracted wide attention due to its
superior sample efficiency. Despite its impressive success so far, it is still
unclear how to appropriately schedule the important hyperparameters to achieve
adequate performance, such as the real data ratio for policy optimization in
Dyna-style model-based algorithms. In this paper, we first theoretically
analyze the role of real data in policy training, which suggests that
gradually increasing the ratio of real data yields better performance.
Inspired by the analysis, we propose a framework named AutoMBPO to
automatically schedule the real data ratio as well as other hyperparameters in
training model-based policy optimization (MBPO) algorithm, a representative
running case of model-based methods. On several continuous control tasks, the
MBPO instance trained with hyperparameters scheduled by AutoMBPO can
significantly surpass the original one, and the real data ratio schedule found
by AutoMBPO shows consistency with our theoretical analysis.

##  1 Introduction

Deep model-free reinforcement learning (MFRL) has achieved great successes in
complex decision-making problems such as Go [26] and robotic control [10], to
name a few. Although MFRL can achieve high asymptotic performance, a
tremendous number of samples collected in interactions with the environment
are required. In contrast, model-based reinforcement learning (MBRL), which
alternately learns a model of the environment and derives an agent with the
help of current model estimation, is considered to be more sample efficient
than MFRL [27].

MBRL methods generally fall into different categories according to the
specific usage of the learned model [25, 33]. Among them, Dyna-style
algorithms [28] adopt some off-the-shelf MFRL methods to train a policy with
both real data from environment and imaginary data generated by the model.
Since Dyna-style algorithms can seamlessly take advantage of innovations in
MFRL literature and have recently shown impressive performance [12, 30], this
paper mainly focuses on Dyna-style algorithms. Although some theoretical
analysis [17, 12, 22] and thorough empirical evaluation [30] have been
conducted in the previous literature, it remains unclear how to appropriately
schedule the hyperparameters to achieve optimum performance when training a
Dyna-style MBRL algorithm. In practice, many important hyperparameters may
primarily affect performance.

Firstly, since Dyna-style MBRL algorithms consist of model learning and policy
optimization, how to balance the alternate optimization of these two parts is
a crucial problem [22]. Intuitively, insufficient model training may lead to
inaccurate estimation while excessive model training will cause overfitting,
and a small amount of policy training may not utilize the model adequately
while exhaustive policy training will exploit the model’s deficiencies.
Moreover, the rollout length of the imaginary trajectory is also critical [12]
since too short rollout length fails to sufficiently leverage the model to
plan forward, while too long rollout length may bring disastrous compounding
error [1]. Janner et al. [12] manually design a schedule that linearly
increases the rollout length across epochs. Finally, when using both real
samples and imaginary samples to train the policy, how to control the ratio of
the two datasets remains unclear. Janner et al. [12] fix the ratio of real
data as $5\%$, while Kalweit and Boedecker [13] use the uncertainty to
adaptively choose the ratio, which tends to use more imaginary samples
initially and gradually use more real samples afterward. According to these
existing works, the optimal hyperparameters in model-based methods may be
dynamic during the whole training process, which further strengthens the
burden of manually scheduling the hyperparameters.

Based on these considerations, in this work, we aim to investigate how to
appropriately schedule these hyperparameters, i.e., real data ratio, model
training frequency, policy training iteration, and rollout length, to achieve
optimal performance of Dyna-style MBRL algorithms. Although real data ratio is
an essential factor empirically [13, 12], it has not yet been studied
thoroughly in theory. To bridge this gap, we first derive a return discrepancy
upper bound for model-based value iteration, which reveals that _gradually
increasing the ratio of real data yields a tighter bound than choosing a fixed
value_. Inspired by the analysis, considering the complex interplay between
real data ratio and other hyperparameters in practice, we develop AutoMBPO to
automatically determine the joint schedule of the above hyperparameters in
model-based policy optimization (MBPO), a representative running case of Dyna-
style MBRL algorithms. Specifically, AutoMBPO introduces a parametric hyper-
controller to sequentially choose the value of hyperparameters in the whole
optimization process to maximize the performance of MBPO. We apply AutoMBPO to
several continuous control tasks. On all these tasks, the MBPO instance
trained with hyperparameters scheduled by AutoMBPO can significantly surpass
the one with original configuration [12]. Furthermore, the hyperparameter
schedule found by AutoMBPO is consistent with our theoretical analysis, which
yields a better understanding of MBPO and provides insights for the design of
other MBRL methods.

##  2 Related Work

Model-based reinforcement learning (MBRL) has been widely studied due to its
high sample efficiency. MBRL methods can be roughly categorized into four
types according to different model usage: (i) Dyna-style algorithms [28, 17,
4, 14, 12, 15] leverage the model to generate imaginary data and adopt some
off-the-shelf MFRL algorithms to train a policy using both real data and
imaginary data; (ii) shooting algorithms [19, 3] use model predictive control
(MPC) to plan directly without explicit policy; (iii) analytic-gradient
algorithms [7, 16, 5] search policies with back-propagation through time by
exploiting the model derivatives; (iv) model-augmented value expansion
algorithms [9, 2] utilize model rollouts to improve the target value for
temporal difference (TD) updates. This paper mainly focuses on the first
category.

Dyna-style algorithms [28] typically alternate between sampling real data,
learning dynamics models, generating imaginary data, and optimizing policy.
Many efforts have been made to improve one or more of these steps. For model
learning, the deep ensemble technique [21, 3, 14] has been leveraged to resist
overfitting. As for imaginary data generation, it is suggested to generate
short rollouts to reduce compounding model error [12], and the model’s
uncertainty is further incorporated to choose reliable imaginary data
dynamically [20]. When optimizing the policy, we can use merely imaginary data
[17, 14] or a mixture of real data and imaginary one in a fixed ratio [12].
Moreover, the ratio can also be dynamically adjusted according to the
estimated uncertainty of the Q-function [13]. Besides, Rajeswaran et al. [22]
present a game-theoretic framework for MBRL by formulating the optimization of
model and policy as a two-player game.

As for automatic hyperparameter optimization for MBRL, Zhang et al. [32]
utilize population based training (PBT) to tune the hyperparameters of a
shooting algorithm dynamically. The reinforcement on reinforcement (RoR)
framework proposed by Dong et al. [8] applies an RL algorithm to control the
sampling and training process of an MBRL method. Our work differs from theirs
in four major aspects: (1) We provide a theoretical analysis for the
hyperparameter being scheduled, which is not included in their work. (2) In
terms of problem settings, RoR assumes one can directly access an initial
state without further interaction from that state, which is impractical for
most RL problems. (3) As for the method design, our formulated hyper-MDP in
Section 5.1 is totally different from RoR’s, including state, action and
reward. In detail, our state definition contains more training information,
and the reward design is based on the return rather than a single-step reward
in the target-MDP. (4) Finally, in the experiments in Section 6, our method
achieves better performance and generalization in all environments, and the
empirical findings also show differences.

##  3 Preliminaries

We first briefly introduce the RL problem with notations used throughout this
paper and the previous MBRL method on which our proposed framework is based.

###  3.1 Reinforcement Learning

In reinforcement learning (RL), there originally exists a Markov decision
process (MDP), which we will call _target-MDP_ throughout the paper. The
target-MDP is defined by the tuple ($\mathcal{S},\mathcal{A},T,r,\gamma$).
$\mathcal{S}$ and $\mathcal{A}$ are the state and action spaces, respectively,
and $T(s^{\prime}\mid s,a)$ is the transition density of state $s^{\prime}$
given action $a$ taken under state $s$. The reward function is denoted as
$r(s,a)$, and $\gamma\in(0,1)$ is the discount factor. The goal of RL is to
find the optimal policy $\pi^{*}$ that maximizes the expected return, denoted
by $\eta$:

$\small\pi^{*}\vcentcolon=\mathop{\arg\max}_{\pi}\eta[\pi]=\mathop{\arg\max}_{\pi}\mathbb{E}_{\pi}\Big{[}\sum_{t=0}^{\infty}\gamma^{t}r(s_{t},a_{t})\Big{]},$
(1)

where $s_{t+1}\sim T(s\mid s_{t},a_{t})$ and $a_{t}\sim\pi(a\mid s_{t})$.
Generally, the ground truth transition $T$ is unknown, and MBRL methods aim to
construct a model $\hat{T}$ of the transition dynamics to help improve the
policy.

To evaluate a reinforcement learning method (no matter whether it is model-
free or model-based), we mainly care about two metrics: (1) asymptotic
performance: the expected return of the converged policy in real environments;
(2) sample efficiency: the number of samples $(s,a,s^{\prime},r)$ collected in
target-MDP to achieve some performance. Therefore, an effective RL method
means it achieves higher expected returns with fewer real samples.

###  3.2 Model-Based Policy Optimization

Model-based policy optimization (MBPO) [12] is a state-of-the-art MBRL method
which will be used as the running case in our framework. MBPO learns a
bootstrapped ensemble of probabilistic dynamics models
$\\{\hat{T}_{\theta}^{1},...,\hat{T}_{\theta}^{E}\\}$ where $E$ is the
ensemble size. Each individual dynamics model in the ensemble is a
probabilistic neural network that outputs a Gaussian distribution with
diagonal covariance
$\hat{T}_{\theta}^{i}\left(s^{\prime}|s,a\right)=\mathcal{N}\left({\mu}_{\theta}^{i}\left({s},{a}\right),{\Sigma}_{\theta}^{i}\left({s},{a}\right)\right)$.
These models are trained independently via maximum likelihood with different
initializations and training batches. The corresponding loss function is

$\displaystyle\mathcal{L}_{\hat{T}}({\theta})=$
$\displaystyle\sum_{n=1}^{N}\left[\mu_{\theta}\left({s}_{n},{a}_{n}\right)-{s}_{n+1}\right]^{\top}{\Sigma}_{\theta}^{-1}\left({s}_{n},{a}_{n}\right)\left[\mu_{\theta}\left({s}_{n},{a}_{n}\right)-{s}_{n+1}\right]+\log\operatorname{det}{\Sigma}_{\theta}\left({s}_{n},{a}_{n}\right).$

In practical implementation, an early stopping trick is adopted in training
the model ensemble. To be more specific, when training each individual model,
a hold-out dataset will be created, and the training will early stop if the
loss evaluated on the hold-out data does not decrease.

Every time the ensemble models are trained, they are then used to generate
$k$-length imaginary rollouts, which begin from states sampled from real data
and follow the current policy. Each time the policy interacts with the target-
MDP, it will be trained for $G$ gradient updates using both real data and
imaginary data in a certain ratio. To make it clear, in each batch for policy
training, we denote the proportion of real data to the whole batch as $\beta$
and the proportion of imaginary data as $1-\beta$. We call $\beta$ the _real
ratio_ and call $G$ the _policy training iteration_ throughout the paper. The
policy optimization algorithm is SAC (Soft Actor-Critic) [11], and the loss
functions for the actor $\pi_{\omega}$ and the critic $Q_{\phi}$ are

$\displaystyle\mathcal{L}_{\pi}(\omega)=~{}$
$\displaystyle\mathbb{E}_{s_{t}\sim
D}\Big{[}\mathbb{E}_{a_{t}\sim\pi_{\omega}}[\alpha\log(\pi_{\omega}(a_{t}|s_{t}))-Q_{\phi}(s_{t},a_{t})]\Big{]},$
(2) $\displaystyle\mathcal{L}_{Q}(\phi)=~{}$
$\displaystyle\mathbb{E}_{(s_{t},a_{t})\sim
D}\Big{[}\frac{1}{2}\big{(}Q_{\phi}(s_{t},a_{t})-(r(s_{t},a_{t})$ (3)
$\displaystyle+\gamma\mathbb{E}_{s_{t+1},a_{t+1}}[Q_{\hat{\phi}}(s_{t+1},a_{t+1})-\alpha\log\pi_{\omega}(a_{t+1}|s_{t+1})])\big{)}^{2}\Big{]},$

where $Q_{\hat{\phi}}$ is the target Q-function.

##  4 Analysis of real ratio schedule

Previous theoretical findings in MBRL mainly focus on the optimization of
model and policy [17, 22] or the rollout length [12, 15], while the usage of
real data in policy optimization remains unclear. Kalweit and Boedecker [13]
only use real data near convergence but accept noisier imaginary data
initially, which seems confusing since the model may not be accurate enough at
the beginning. In this section, we provide a theoretical analysis of the real
ratio schedule by deriving a return discrepancy upper bound of sampling-based
Fitted Value Iteration (FVI) in model-based settings. The analysis mainly
follows Munos and Szepesvári [18], where all the training data are sampled
from the underlying dynamics. We extend the framework to a more general case:
the data can either be sampled from the underlying dynamics w.p. (with
probability) $\beta$ or from a learned model w.p. $1-\beta$, which we call
_$\beta$ -mixture sampling-based FVI_. Besides, we consider deterministic
environments here, i.e., given $(s,a)$, the next state $s^{\prime}=T(s,a)$ is
unique, just as in all the environments in Section 6.

For the sake of specificity, we first present a detailed description of the
$\beta$-mixture sampling-based FVI algorithm. Let $\mathcal{F}$ be the value
function space, and $V_{0}\in\mathcal{F}$ be the initial value function. The
FVI algorithm produces a sequence of functions $\\{V_{k}\\}_{0\leq k\leq
K}\subset\mathcal{F}$ iteratively. Given $V_{k}$ and $N$ sampled states
$\\{s_{1},...,s_{N}\\}$, the value function of the next iteration $V_{k+1}$ is
computed as

$\displaystyle\hat{V}_{k}(s_{i})$
$\displaystyle=\max_{a\in\mathcal{A}}\big{(}r_{i}+\gamma
V_{k}\left(s_{i}^{\prime}\right)\big{)},i=1,2,\ldots,N,$ (4) $\displaystyle
V_{k+1}$
$\displaystyle=\underset{f\in\mathcal{F}}{\operatorname{argmin}}\sum_{i=1}^{N}\left|f\left(s_{i}\right)-\hat{V}_{k}(s_{i})\right|^{p},$
(5)

where $s_{i}$ is uniformly sampled from a certain state distribution $\mu$,
and for each possible action $a$, the next state $s_{i}^{\prime}$ is sampled
from the underlying dynamics w.p. $\beta$ or from the learned model w.p.
$1-\beta$. As final preparatory steps, we define the _Bellman operator_ for
deterministic MDPs $B:\mathcal{F}\rightarrow\mathcal{F}$ as

$BV(s)=\max_{a\in\mathcal{A}}\left\\{r(s,a)+\gamma V(s^{\prime})\right\\},$

where $s^{\prime}=T(s,a)$; and for function $g$ over $\mathcal{X}$,
$\|g\|_{p,\mu}$ is defined as $\|g\|_{p,\mu}^{p}=\int|g(x)|^{p}\mu(dx)$.
Following Munos and Szepesvári [18], we begin with the error bound in a single
iteration.

######  Lemma 4.1.

(Single Iteration Error Bound) Let $V_{k}$ and $V_{k+1}$ be the value
functions of iteration $k$ and $k+1$, and $V_{max}=r_{max}/(1-\gamma)$. For
$p\geq 1$ and a certain state distribution $\mu$, let the inherent Bellman
error of the value function space $\mathcal{F}$ be defined by
$d_{p,\mu}(B\mathcal{F},\mathcal{F})=\sup_{V\in\mathcal{F}}\inf_{f\in\mathcal{F}}\|f-BV\|_{p,\mu}$.
Assume the value functions are $K_{V}$-$Lipschitz$ continuous, i.e., for any k
and states pair $(s_{i},s_{j})$, it holds that
$|V_{k}(s_{i})-V_{k}(s_{j})|\leq K_{V}\|s_{i}-s_{j}\|_{2}$. And assume the
$L^{2}$-norm model error between the real next state $s^{\prime}$ and the
predicted ones $\hat{s^{\prime}}$ obeys a half-normal distribution, i.e., the
probability density function
$f(\|\hat{s^{\prime}}-s^{\prime}\|_{2})=\frac{2}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\|\hat{s^{\prime}}-s^{\prime}\|_{2}^{2}}{2\sigma^{2}}\right)$.
Let
$\mathcal{N}_{0}(N)=\mathcal{N}\left(\frac{1}{8}\left(\frac{\epsilon_{0}}{4}\right)^{p},\mathcal{F},N,\mu\right)$
denotes the covering number defined in Lemma A.1. Then for any
$\epsilon_{0},\delta_{0}>0$, the inequality

$\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}\leq
d_{p,\mu}(B\mathcal{F},\mathcal{F})+\epsilon_{0}$

holds w.p. at least $1-\delta_{0}$ provided that

$N>128\left(\frac{8V_{\max}}{\epsilon_{0}}\right)^{2p}\Big{(}\log(1/\delta_{0})+\log\left(32\mathcal{N}_{0}(N)\right)\Big{)}$
(6)

and

$\sigma<\frac{\epsilon_{0}}{4\gamma
K_{V}\Phi^{-1}\Big{(}1-\delta_{0}/\big{(}8N|A|(1-\beta)\big{)}\Big{)}}.$ (7)

###### Proof.

See _Appendix A_, _Lemma A.2_. ∎

According to the model error assumption (more details in Appendix B), the
model will give precise predictions if $\sigma$ is small enough. Therefore,
Lemma 4.1 states that with an accurate model and sufficient large $N$, the
errors between $BV_{k}$ and $V_{k+1}$ in each iteration can be bounded with
high probability. Now let us turn to bound the return discrepancy between the
greedy policy w.r.t. $V_{K}$ and the optimal policy.

######  Theorem 4.1.

($\beta$-mixture sampling-based FVI bound) Assume that the pseudo-dimension of
the value function space is finite as in Proposition A.1, under the
concentrability assumption and the same assumptions of Lemma 4.1, let $\rho$
be an arbitrary state distribution, $C_{\rho,\mu}$ be the discounted-average
concentrability coefficient and $\pi_{K}$ be the greedy policy w.r.t. $V_{K}$.
Let $V^{\pi_{K}}$ and $V^{*}$ be the expected return of executing $\pi_{K}$
and the optimal policy $\pi^{*}$ in real environment, respectively. Define
$N_{\rm{real}}=N\cdot|\mathcal{A}|\cdot\beta$ as the expected number of real
samples. Then the following bound holds w.p. at least $1-\delta$:

$\displaystyle\|V^{*}-V^{\pi_{K}}\|_{p,\rho}\leq$
$\displaystyle\frac{2\gamma}{(1-\gamma)^{2}}C_{\rho,\mu}^{1/p}d_{p,\mu}(B\mathcal{F},\mathcal{F})+O\left(\Big{(}\frac{\beta|\mathcal{A}|}{N_{\rm{real}}}\big{(}\log(\frac{N_{\rm{real}}}{\beta|\mathcal{A}|})+\log(\frac{K}{\delta})\big{)}\Big{)}^{\frac{1}{2p}}\right)$
(8)
$\displaystyle+O\Big{(}\Phi^{-1}\big{(}1-\frac{\beta\delta}{8KN_{\rm{real}}(1-\beta)}\big{)}\sigma\Big{)}+O\left(\gamma^{K/p}V_{\max}\right).$

###### Proof.

See _Appendix A_, _Theorem A.1_. ∎

Remark. The definition of the concentrability assumption is provided in
Appendix A. Theorem 4.1 gives an upper bound on the return discrepancy. To
achieve small return discrepancy, we hope to minimize the upper bound. For the
first term in the bound, the approximation error
$d_{p,\mu}(B\mathcal{F},\mathcal{F})$ can be made small by selecting
$\mathcal{F}$ to be large enough, and the last term arises due to the finite
number of iterations. For the remaining two terms, as $\beta/N_{\rm{real}}$
increases, the second term increases while the third term decreases.
Therefore, there exists an optimal value for $\beta/N_{\rm{real}}$ to trade
off these two terms. In online MBRL, $N_{\rm{real}}$ increases as real samples
are continuously collected from the environment throughout the training
process, so that gradually increasing the real ratio $\beta$ is promising to
achieve good performance according to the upper bound.

From another perspective, the second term corresponds to the estimation error
of using limited data to learn the value functions. By utilizing a learned
model, we can generate as much data as we want, which, however, tends to
suffer from the model error represented by the third term. The hyperparameter
$\beta$ serves as a trade-off between these two kinds of errors. Intuitively,
with more real data available, the pressure of limited data is alleviated,
which means we can increase $\beta$ accordingly to reduce the reliance on
inaccurate models.

##  5 The AutoMBPO Framework

The analysis in the previous section shows that the optimal real ratio $\beta$
during training should gradually increase according to the current
optimization states, such as the number of real samples and model error. These
states may also be affected by other hyperparameters, e.g., model training
frequency and rollout length, which makes the hyperparameter scheduling non-
trivial in practice. In order to investigate the optimal joint scheduling of
these hyperparameters, we propose a practical framework named AutoMBPO to
schedule these hyperparameters automatically when training an MBRL algorithm.
The high-level idea of AutoMBPO is to formulate hyperparameter scheduling as a
sequential decision-making problem and then adopt some RL algorithm to solve
it. Here we choose MBPO [12] as a representative running case since MBPO has
achieved remarkable success in continuous control benchmark tasks and can
naturally take all the hyperparameters we care about into consideration (a
comparison of the hyperparameters that can be scheduled in different MBRL
algorithms is provided in Appendix C). We illustrate the overview framework of
AutoMBPO in Figure 1 and detail the algorithm in Algorithm 5.1.

###  5.1 Hyper-MDP Formulation

We formulate the problem of determining the hyperparameter schedule as a new
sequential decision-making problem, which we will call _hyper-MDP_. In the
hyper-MDP, the hyper-controller needs to decide the hyperparameter values
throughout one complete run of the MBPO learning process. Below we elaborate
on several critical aspects of the hyper-MDP formulation used in this paper.

Episode. An episode of the hyper-MDP consists of the whole process of training
an MBPO instance from scratch, which contains multiple episodes of the target-
MDP. We will terminate an episode of the hyper-MDP if $N$ samples have been
collected in the target-MDP. To be more specific, if the episode length of the
target-MDP is $H$, and we only train MBPO for $m$ episodes, then we have
$N=H\cdot m$.

Figure 1: Illustration of AutoMBPO framework. The hyper-controller (left)
treats the whole learning process of MBPO (right) as its hyper-MDP and takes
action according to the state containing all the information of MBPO training.
The action is to adjust the hyperparameters of MBPO, e.g., real ratio, model
training frequency, policy training iteration, and rollout length. Afterward,
a reward is given to the hyper-controller based on the performance of MBPO
policy in target-MDP.

State. The state should be defined to include all the information of MBPO
training to help hyper-controller make decisions accordingly. We use a
concatenation of several features to represent the state. As suggested by
Theorem 4.1, the state should contain $N_{\rm{real}}$, the number of real
samples so far, and the model error, which can be estimated by the model loss
$\mathcal{L}_{\hat{T}}$ on a validation set. Besides, we also add other
features to better capture the current optimization situation: (1) critic
loss: $\mathcal{L}_{Q}$, which to some extent represents how well the value
function fits the training data; (2) policy change:
$\epsilon_{\pi}=E_{(s,a)\sim\pi_{D}}[|\pi(a|s)-\pi_{D}(a|s)|]$, which
evaluates the distance of current policy $\pi$ and the data-collecting policy
$\pi_{D}$ on sampled data. In practice, we calculate $\epsilon_{\pi}$ using
recently collected data to reflect the updating speed of policy roughly; (3)
validation metrics: we use the average return on the target-MDP to represent
the performance of the policy; (4) current hyperparameters: the hyperparameter
values used currently by MBPO. All these features are normalized into $[0,1]$
to make the learning process stable.

Action. The action set in the hyper-MDP should contain all the hyperparameters
we care about. Firstly, we would like to adjust the ratio of real data when
optimizing the policy. Secondly, we also consider investigating the model
training and policy learning frequency to balance their alternate optimization
[22]. Lastly, the rollout length when generating imaginary rollouts is crucial
in MBPO, but it remains unclear how to choose its value [12]. Besides, since
there exists an optimal value for $\beta/N_{\text{real}}$ according to the
theoretical analysis and $N_{\text{real}}$ gradually increases during training
in practice, for the action space, adjusting the real ratio by multiplying a
constant is enough to get an effective schedule. To sum up, the action set and
the corresponding space of each action are: (1) adjusting the real ratio
$\beta$: $\\{\times c^{-1},\times 1,\times c\\}$, where c is a constant larger
than $1$; (2) deciding whether to train the model: $\\{0,1\\}$. Note that in
MBPO, the model training iterations are implicitly decided by an early
stopping trick, so we only decide whether to train the model; (3) adjusting
the policy training iteration $G$: $\\{-1,0,+1\\}$; (4) adjusting the model
rollout length $k$: $\\{-1,0,+1\\}$. In practice, the actions are taken every
$\tau$ real samples are collected where $\tau$ is the hyperparameter of the
hyper-controller.

Reward. One direct way is to define a reward at the end of an episode in
hyper-MDP, which, however, will lead to the sparse reward problem. We define
our reward as follows: every $H$ real samples, i.e., one episode in target-
MDP, we evaluate the MBPO policy on target-MDP and define the reward as the
average evaluated return.

[t]  AutoMBPO Initialize a hyper-controller MC

2 repeat

3 Initialize a MBPO instance; initialize the hyper-controller sample count
$i=0$

4 for  _$m$ target-MDP episodes_ do

5 for  _$h$ = $0:H-1$ real timesteps (an episode in target-MDP)_ do

6 Interact with the target-MDP using the MBPO policy; add the real samples to
$\mathcal{D}_{\mathrm{env}}$

7 if  _$h$ mod $\tau=0$_ then

8 Calculate the current state $\mathbb{S}_{i}$ for hyper-MDP; take an action
based on $\mathbb{S}_{i}$:
$\mathbb{A}_{i}\leftarrow\textbf{MC}(\mathbb{S}_{i})$

9 Decide whether to train the model, and adjust the hyperparameters
$k,G,\beta$ according to $\mathbb{A}_{i}$

10

11 for  _F model rollouts_ do

12 Sample $s_{t}$ uniformly from $\mathcal{D}_{\mathrm{env}}$; perform
$k$-step rollout from $s_{t}$; add the imaginary data to
$\mathcal{D}_{\mathrm{model}}$

13

14 for  _$G$ gradient updates_ do

15 Train the MBPO policy using SAC on $\mathcal{D}_{\mathrm{env}}$ and
$\mathcal{D}_{\mathrm{model}}$ in the ratio $\beta$

16

17 if  _$h=H-1$ _ then

18 Evaluate MBPO and construct the reward as the average return:
$\mathbb{R}_{i}=Avg(\eta)$; update $i=i+1$

19

20 else if  _$h$ mod $\tau=0$ and $h\not=H-\tau$_ then

21 Construct the reward $\mathbb{R}_{i}=0$; update $i=i+1$

22

23

24 Train the hyper-controller MC on the data
$\\{(\mathbb{S}_{i},\mathbb{A}_{i},\mathbb{R}_{i})\\}_{i=0}^{m\cdot H/\tau}$

25until  _acceptable performance is achieved_

###  5.2 Hyper-Controller Learning

After formulating the hyper-MDP, we are now ready to learn a hyper-controller
in the hyper-MDP. In fact, we can use any off-the-shelf reinforcement learning
algorithm, and we choose Proximal Policy Optimization (PPO) [24] as the hyper-
controller algorithm since PPO has achieved considerable success in the
optimization schedule problem [31]. To be more specific, we maximize the
following clipped PPO objective:

$\mathcal{J}(\Omega)=\hat{\mathbb{E}}_{t}\left[\min\big{(}r_{t}(\Omega)\hat{A}_{t},\operatorname{clip}\left(r_{t}(\Omega),1-\epsilon,1+\epsilon\right)\hat{A}_{t}\big{)}\right],$

where
$r_{t}(\Omega)=\textbf{MC}_{\Omega}(a_{t}|s_{t})/\textbf{MC}_{\Omega_{\text{old}}}(a_{t}|s_{t})$
and $\textbf{MC}_{\Omega}$ is the hyper-controller policy parameterized by
$\Omega$.

There are multiple choices of the advantage function [23], and we use the
baseline version of the Monte-Carlo returns to reduce the variance:

$\hat{A}_{t}=\Sigma_{i=t}^{m\cdot H/\tau}(R_{i}-R^{\prime}_{i}).$ (9)

Here, $R^{\prime}_{i}$ is the average return of MBPO policy trained in advance
with the same amount of real data using the original hyperparameters, and we
treat $\Sigma_{i=t}^{m\cdot H/\tau}R^{\prime}_{i}$ as the baseline. The
advantage function can also be regarded as the performance improvement
compared with the original MBPO configuration.

##  6 Experiment

In this section, we conduct experiments to verify the effectiveness of
AutoMBPO and provide comprehensive empirical analysis. Our experiments aim to
answer the following three questions: i) How does the MBPO instance learned by
AutoMBPO perform compared to the original configuration? ii) Do the learned
hyper-controllers in different target-MDPs share similar schedules, and are
the schedules consistent with previous theoretical analysis? iii) Which
hyperparameter in MBPO is more important/sensitive to the final performance?
More experimental results are provided in Appendix D.

Figure 2: Learning curves of AutoMBPO and baselines on different continuous
control environments. The solid lines indicate the mean, and shaded areas
indicate the standard deviation of 10 trials over different random seeds. Each
trial is evaluated every 1000 environment steps, where each evaluation reports
the average return over 10 episodes.

###  6.1 Comparative Evaluation

We mainly compare the MBPO instance found by AutoMBPO to the MBPO with
original hyperparameters [12] since AutoMBPO is directly based on MBPO, and
there is no need to compare to other model-based methods. Besides, we also
compare to Soft Actor-Critic (SAC) [11], the state-of-the-art model-free
algorithm in terms of sample efficiency and asymptotic performance. Actually,
SAC can be viewed as an extreme variant of MBPO where the real ratio in
training policy is 1. To make a fair comparison, we also provide the result of
SAC where the policy training iteration is enforced to be the same as in MBPO.
To be more specific, the original SAC, denoted as SAC(1), trains the policy
one time per step in the real environment, while in MBPO, the policy is
trained 20 times per real step, so we add the baseline SAC(20). We also
compare our method to the PBT [32] and RoR [8] algorithm for a comprehensive
comparison. Note that the hyper-controllers in AutoMBPO and RoR are trained in
advance and are only evaluated here. The whole training process of AutoMBPO is
provided in Appendix D.3.

The compared methods are evaluated on three MuJoCo [29] (Hopper, Ant,
Humanoid), and three PyBullet [6] (HopperBullet, Walker2dBullet,
HalcheetahBullet) continuous control tasks as our target-MDPs111Code is
available at: https://github.com/hanglai/autombpo. In practice, the three
hyperparameters (real ratio, model training frequency, and policy training
iteration) are scheduled on all tasks, while the rollout length is not
scheduled on the PyBullet tasks since preliminary experiments showed that even
2-step rollouts degrade the performance on the PyBullet environments.
Moreover, due to the high computational cost, the hyper-controller is only
trained in previous $m$ episodes of MBPO and will be used to continuously
control MBPO instances running for $M$ episodes, which is usually $2-3$ times
larger than $m$ (e.g., $m=50,M=150$ on Ant). Furthermore, we give a penalty
for the model training to speed up the training process, and we find it does
not influence the results much. More experimental details and hyperparameter
configuration for hyper-controller can be found in Appendix F and H.

The comparison results are shown in Figure 2. We can observe that: i) the MBPO
instances learned by AutoMBPO significantly outperform MBPO with original
configuration in all tasks, highlighting the strength of automatic scheduling.
ii) Compared to other hyperparameter optimization baselines (i.e., PBT and
RoR), our approach AutoMBPO achieves better performance and generalization in
different environments, verifying the effectiveness of our proposed framework.
iii) MFRL methods with exhaustive policy training may overfit to a small
amount of real data since SAC(20) achieves worse performance than SAC(1) on
some complex tasks like Ant and Walker2dBullet. The overfitting problem is
much less severe in MBPO since we can generate plentiful data using the model.
iv) Though only trained in previous episodes, the hyper-controller can still
find reasonable hyperparameter schedules and can be extended to the longer
training process of MBPO.

###  6.2 Hyperparameter Schedule Visualization

We visualize the schedules of the hyperparameters found by AutoMBPO in Figure
3. From the results, the hyperparameter schedules on different tasks show both
similarities and differences: i) A similar increasing schedule of the real
ratio can be observed in all six environments, which is consistent with our
derivation in Section 4 and the empirical results in Kalweit and Boedecker
[13]. ii) The hyper-controller tends to train the model more frequently in
complex environments, such as Ant and Humanoid. iii) To a certain extent,
increasing the policy training iteration can better exploit the model.
However, it may also increase the risk of instability [22], and we observe
that the hyper-controller decreases the policy training iteration in the
subsequent training phase of Hopper when the policy is close to convergence.
iv) As for the rollout length, the hyper-controller adopts a nearly linear
increasing schedule, which is close to the manually designed one used in the
original MBPO [12].

Figure 3: Hyperparameter schedule found by AutoMBPO on different tasks. The
dashed lines indicate the original hyperparameter configuration in MBPO. For
better visualization, model training frequency is scaled by the frequency of
the original MBPO, and the training process (x-axis) is scaled to [0,1] since
the max episodes vary on different tasks.

###  6.3 Importance of Hyperparameters

In this section, we investigate the importance of different hyperparameters.
Specifically, we utilize the hyper-controller trained previously to schedule
one certain hyperparameter respectively, i.e., real ratio (AutoMBPO-R), model
training frequency (AutoMBPO-M), policy training iteration (AutoMBPO-P), and
rollout length (AutoMBPO-L), while fixing other hyperparameters to the same as
in the original MBPO. As shown in Figure 4, the importance of various
hyperparameters changes across different tasks: using the hyper-controller to
schedule the policy training iteration is helpful on Hopper but fails on Ant,
while the situation reverses for model training frequency. However, using the
hyper-controller to schedule the real ratio retains much of the benefit of
AutoMBPO on almost all tasks, which further verifies the importance of real
ratio.

Figure 4: Performance of the instances with only one specific hyperparameter
scheduled by hyper-controller while other hyperparameters fixed to the same as
in MBPO. The complete figures of all six environments are provided in Appendix
D.1.

##  7 Conclusion

In this paper, we conduct a theoretical analysis to justify real data usage
for policy optimization in MBRL algorithms. Inspired by the analysis, we
present a framework called AutoMBPO to automatically schedule the
hyperparameters in the MBPO algorithm. Extensive experiments on complex
continuous control tasks demonstrate that the derived MBPO instance with
hyperparameters scheduled by AutoMBPO can achieve significantly better
performance than the MBPO instance of the original configuration. Moreover,
the real ratio schedule discovered by AutoMBPO is consistent with our
theoretical analysis. We hope this insight can help design more effective MBRL
algorithms. For future work, we plan to apply our hyper-controller training
scheme to other MBRL methods.

## Acknowledgments

The SJTU team is supported by “New Generation of AI 2030” Major Project
(2018AAA0100900), Shanghai Municipal Science and Technology Major Project
(2021SHZDZX0102) and National Natural Science Foundation of China (62076161,
61772333, 61632017). The work is also sponsored by Huawei Innovation Research
Program.

## References

  * Asadi et al. [2018] Kavosh Asadi, Dipendra Misra, and Michael Littman.  Lipschitz continuity in model-based reinforcement learning.  In _International Conference on Machine Learning_ , pages 264–273. PMLR, 2018.
  * Buckman et al. [2018] Jacob Buckman, Danijar Hafner, George Tucker, Eugene Brevdo, and Honglak Lee.  Sample-efficient reinforcement learning with stochastic ensemble value expansion.  In _Advances in Neural Information Processing Systems_ , pages 8224–8234, 2018.
  * Chua et al. [2018] Kurtland Chua, Roberto Calandra, Rowan McAllister, and Sergey Levine.  Deep reinforcement learning in a handful of trials using probabilistic dynamics models.  In _Advances in Neural Information Processing Systems_ , pages 4754–4765, 2018.
  * Clavera et al. [2018] Ignasi Clavera, Jonas Rothfuss, John Schulman, Yasuhiro Fujita, Tamim Asfour, and Pieter Abbeel.  Model-based reinforcement learning via meta-policy optimization.  In _Conference on Robot Learning_ , pages 617–629, 2018.
  * Clavera et al. [2020] Ignasi Clavera, Violet Fu, and Pieter Abbeel.  Model-augmented actor-critic: Backpropagating through paths.  _arXiv preprint arXiv:2005.08068_ , 2020.
  * Coumans and Bai [2016] Erwin Coumans and Yunfei Bai.  Pybullet, a python module for physics simulation for games, robotics and machine learning.  2016\.
  * Deisenroth and Rasmussen [2011] Marc Deisenroth and Carl E Rasmussen.  Pilco: A model-based and data-efficient approach to policy search.  In _Proceedings of the 28th International Conference on machine learning (ICML-11)_ , pages 465–472, 2011.
  * Dong et al. [2020] Linsen Dong, Yuanlong Li, Xin Zhou, Yonggang Wen, and Kyle Guan.  Intelligent trainer for dyna-style model-based deep reinforcement learning.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Feinberg et al. [2018] Vladimir Feinberg, Alvin Wan, Ion Stoica, Michael I Jordan, Joseph E Gonzalez, and Sergey Levine.  Model-based value estimation for efficient model-free reinforcement learning.  _arXiv preprint arXiv:1803.00101_ , 2018.
  * Gu et al. [2017] Shixiang Gu, Ethan Holly, Timothy Lillicrap, and Sergey Levine.  Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates.  In _2017 IEEE international conference on robotics and automation (ICRA)_ , pages 3389–3396. IEEE, 2017.
  * Haarnoja et al. [2018] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  _arXiv preprint arXiv:1801.01290_ , 2018.
  * Janner et al. [2019] Michael Janner, Justin Fu, Marvin Zhang, and Sergey Levine.  When to trust your model: Model-based policy optimization.  In _Advances in Neural Information Processing Systems_ , pages 12519–12530, 2019.
  * Kalweit and Boedecker [2017] Gabriel Kalweit and Joschka Boedecker.  Uncertainty-driven imagination for continuous deep reinforcement learning.  In _Conference on Robot Learning_ , pages 195–206, 2017.
  * Kurutach et al. [2018] Thanard Kurutach, Ignasi Clavera, Yan Duan, Aviv Tamar, and Pieter Abbeel.  Model-ensemble trust-region policy optimization.  In _International Conference on Learning Representations_ , 2018.
  * Lai et al. [2020] Hang Lai, Jian Shen, Weinan Zhang, and Yong Yu.  Bidirectional model-based policy optimization.  In _International Conference on Machine Learning_ , pages 5618–5627. PMLR, 2020.
  * Levine and Koltun [2013] Sergey Levine and Vladlen Koltun.  Guided policy search.  In _International Conference on Machine Learning_ , pages 1–9, 2013\.
  * Luo et al. [2018] Yuping Luo, Huazhe Xu, Yuanzhi Li, Yuandong Tian, Trevor Darrell, and Tengyu Ma.  Algorithmic framework for model-based deep reinforcement learning with theoretical guarantees.  In _International Conference on Learning Representations_ , 2018.
  * Munos and Szepesvári [2008] Rémi Munos and Csaba Szepesvári.  Finite-time bounds for fitted value iteration.  _Journal of Machine Learning Research_ , 9(May):815–857, 2008.
  * Nagabandi et al. [2018] Anusha Nagabandi, Gregory Kahn, Ronald S Fearing, and Sergey Levine.  Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 7559–7566. IEEE, 2018.
  * Pan et al. [2020] Feiyang Pan, Jia He, Dandan Tu, and Qing He.  Trust the model when it is confident: Masked model-based actor-critic.  _arXiv preprint arXiv:2010.04893_ , 2020.
  * Rajeswaran et al. [2016] Aravind Rajeswaran, Sarvjeet Ghotra, Balaraman Ravindran, and Sergey Levine.  Epopt: Learning robust neural network policies using model ensembles.  _arXiv preprint arXiv:1610.01283_ , 2016.
  * Rajeswaran et al. [2020] Aravind Rajeswaran, Igor Mordatch, and Vikash Kumar.  A game theoretic framework for model based reinforcement learning.  _arXiv preprint arXiv:2004.07804_ , 2020.
  * Schulman et al. [2015] John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel.  High-dimensional continuous control using generalized advantage estimation.  _arXiv preprint arXiv:1506.02438_ , 2015.
  * Schulman et al. [2017] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ , 2017.
  * Shen et al. [2020] Jian Shen, Han Zhao, Weinan Zhang, and Yong Yu.  Model-based policy optimization with unsupervised model adaptation.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Silver et al. [2016] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  _nature_ , 529(7587):484–489, 2016.
  * Sun et al. [2019] Wen Sun, Nan Jiang, Akshay Krishnamurthy, Alekh Agarwal, and John Langford.  Model-based rl in contextual decision processes: Pac bounds and exponential improvements over model-free approaches.  In _Proceedings of the Thirty-Second Conference on Learning Theory_ , pages 2898–2933, 2019.
  * Sutton [1990] Richard S Sutton.  Integrated architectures for learning, planning, and reacting based on approximating dynamic programming.  In _Machine learning proceedings 1990_ , pages 216–224. Elsevier, 1990.
  * Todorov et al. [2012] Emanuel Todorov, Tom Erez, and Yuval Tassa.  Mujoco: A physics engine for model-based control.  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 5026–5033. IEEE, 2012.
  * Wang et al. [2019] Tingwu Wang, Xuchan Bao, Ignasi Clavera, Jerrick Hoang, Yeming Wen, Eric Langlois, Shunshi Zhang, Guodong Zhang, Pieter Abbeel, and Jimmy Ba.  Benchmarking model-based reinforcement learning.  _arXiv preprint arXiv:1907.02057_ , 2019.
  * Xu et al. [2018] Haowen Xu, Hao Zhang, Zhiting Hu, Xiaodan Liang, Ruslan Salakhutdinov, and Eric Xing.  Autoloss: Learning discrete schedule for alternate optimization.  In _International Conference on Learning Representations_ , 2018.
  * Zhang et al. [2021] Baohe Zhang, Raghu Rajan, Luis Pineda, Nathan Lambert, André Biedenkapp, Kurtland Chua, Frank Hutter, and Roberto Calandra.  On the importance of hyperparameter optimization for model-based reinforcement learning.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4015–4023. PMLR, 2021.
  * Zhu et al. [2020] Guangxiang Zhu, Minghao Zhang, Honglak Lee, and Chongjie Zhang.  Bridging imagination and reality for model-based deep reinforcement learning.  _Advances in Neural Information Processing Systems_ , 33, 2020.

## Appendix for: On Effective Scheduling of
Model-based Reinforcement Learning

##  Appendix A Proofs

######  Assumption A.1.

(Discounted-average concentrability of future-state distributions) Given
$\rho,\mu,m\geq 1$, and an arbitrary sequence of stationary policies
$\\{\pi_{m}\\}_{m\geq 1}$, assume that the future-state distribution $\rho
T^{\pi_{1}}T^{\pi_{2}}\ldots T^{\pi_{m}}$ is absolutely continuous w.r.t.
$\mu$. Assume that

$c(m)\stackrel{{\scriptstyle\rm
def}}{{=}}\sup_{\pi_{1},\ldots,\pi_{m}}\left\|\frac{d\left(\rho
T^{\pi_{1}}T^{\pi_{2}}\ldots T^{\pi_{m}}\right)}{d\mu}\right\|_{\infty}$

satisfies

$C_{\rho,\mu}\stackrel{{\scriptstyle\rm def}}{{=}}(1-\gamma)^{2}\sum_{m\geq
1}m\gamma^{m-1}c(m)<+\infty.$

We call $c(m)$ the _$m$ -step concentrability_ of a future-state distribution
and call $C_{\rho,\mu}$ the _discounted-average concentrability coefficient_
of the future-state distributions. The class of MDPs that satisfies this
concentrability assumption is quite large, which is further discussed in Munos
and Szepesvári [18].

######  Lemma A.1.

(Pollard, 1984) Let $\mathcal{F}$ be a set of measurable functions
$f:\mathcal{X}\rightarrow[0,K]$ and let $\epsilon>0$, $N$ be arbitrary. If
$X_{i}$, $i=1,\ldots,N$ is an i.i.d. sequence taking values in the space
$\mathcal{X}$ then

$\mathbb{P}\left(\sup_{f\in\mathcal{F}}\left|\frac{1}{N}\sum_{i=1}^{N}f\left(X_{i}\right)-\mathbb{E}\left[f\left(X_{1}\right)\right]\right|>\epsilon\right)\leq
8\mathbb{E}\left[\mathcal{N}\left(\epsilon/8,\mathcal{F}\left(X^{1:N}\right)\right)\right]e^{-\frac{N\epsilon^{2}}{128K^{2}}},$

where $\mathcal{N}_{q}\left(\epsilon,\mathcal{F}\left(X^{1:N}\right)\right)$
denotes the ($\epsilon$,q)-covering number of the set
$\mathcal{F}(X^{1:N})=\left\\{\left(f\left(X_{1}\right),\ldots,f\left(X_{N}\right)\right)\mid
f\in\mathcal{F}\right\\}$, i.e., the smallest integer $m$ such that
$\mathcal{F}(X^{1:N})$ can be covered by $m$ balls of the normed-space
$\left(\mathbb{R}^{N},\|\cdot\|_{q}\right)$ with centers in
$\mathcal{F}(X^{1:N})$ and radius $N^{1/q}\epsilon$. When $X^{1:N}$ are i.i.d
with the distribution $\mu$, we denote
$\mathbb{E}\left[\mathcal{N}_{q}\left(\epsilon,\mathcal{F}\left(X^{1:N}\right)\right)\right]$
as $\mathcal{N}_{q}(\epsilon,\mathcal{F},N,\mu)$. And when $q=1$,
$\mathcal{N}$ is used instead of $\mathcal{N}_{1}$.

######  Proposition A.1.

(Haussler $(1995)$, Corollary 3$)$ For any set $X$, any points $x^{1:N}\in
X^{N}$, any class $\mathcal{F}$ of functions on $\mathcal{X}$ taking values in
$[0,L]$ with pseudo-dimension $V_{\mathcal{F}^{+}}<\infty$, and any
$\epsilon>0$,

$\mathcal{N}\left(\epsilon,\mathcal{F}\left(x^{1:N}\right)\right)\leq
e\left(V_{\mathcal{F}+}+1\right)\left(\frac{2eL}{\epsilon}\right)^{V_{\mathcal{F}+}}.$

######  Lemma A.2.

(Single Iteration Error Bound) Let $V_{k}$ and $V_{k+1}$ be the value
functions of iteration $k$ and $k+1$, and $V_{max}=r_{max}/(1-\gamma)$. For
$p\geq 1$ and a certain state distribution $\mu$, let the inherent Bellman
error of the value function space $\mathcal{F}$ be defined by
$d_{p,\mu}(B\mathcal{F},\mathcal{F})=\sup_{V\in\mathcal{F}}\inf_{f\in\mathcal{F}}\|f-BV\|_{p,\mu}$.
Assume the value functions are $K_{V}$-$Lipschitz$ continuous, i.e., for any k
and states pair $(s_{i},s_{j})$, it holds that
$|V_{k}(s_{i})-V_{k}(s_{j})|\leq K_{V}\|s_{i}-s_{j}\|_{2}$. And assume the
$L^{2}$-norm model error between the real next state $s^{\prime}$ and the
predicted ones $\hat{s^{\prime}}$ obeys a half-normal distribution, i.e., the
probability density function
$f(\|\hat{s^{\prime}}-s^{\prime}\|_{2})=\frac{2}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\|\hat{s^{\prime}}-s^{\prime}\|_{2}^{2}}{2\sigma^{2}}\right)$.
Let
$\mathcal{N}_{0}(N)=\mathcal{N}\left(\frac{1}{8}\left(\frac{\epsilon_{0}}{4}\right)^{p},\mathcal{F},N,\mu\right)$
denotes the covering number defined in Lemma A.1. Then for any
$\epsilon_{0},\delta_{0}>0$, the inequality

$\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}\leq
d_{p,\mu}(B\mathcal{F},\mathcal{F})+\epsilon_{0}$

holds w.p. at least $1-\delta_{0}$ provided that

$N>128\left(\frac{8V_{\max}}{\epsilon_{0}}\right)^{2p}\Big{(}\log(1/\delta_{0})+\log\left(32\mathcal{N}_{0}(N)\right)\Big{)}$
(10)

and

$\sigma<\frac{\epsilon_{0}}{4\gamma
K_{V}\Phi^{-1}\Big{(}1-\delta_{0}/\big{(}8N|A|(1-\beta)\big{)}\Big{)}}.$ (11)

###### Proof.

Let $f^{*}$ be the best fit to $BV_{k}$ in $\mathcal{F}$:
$f^{*}=\Pi_{\mathcal{F}}BV_{k}$, and let $\hat{\mu}$ denotes the distribution
of the $N$ sampled states $\\{s_{i}\\}_{i=1}^{N}$. Define
$\|\cdot\|_{p,\hat{\mu}}$ as:

$\|f\|^{p}_{p,\hat{\mu}}=\frac{1}{N}\sum_{i=1}^{N}\left|f\left(s_{i}\right)\right|^{p}.$

If the following inequalities hold simultaneously w.p. not smaller than
$1-\delta_{0}$, then the Lemma can be proved by choosing
$\epsilon_{0}^{\prime}=\epsilon_{0}/4$:

$\displaystyle\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}$
$\displaystyle\leq\left\|V_{k+1}-BV_{k}\right\|_{p,\hat{\mu}}+\epsilon_{0}^{\prime}$
(12)
$\displaystyle\leq\|V_{k+1}-\hat{V}_{k}\|_{p,\hat{\mu}}+2\epsilon_{0}^{\prime}$
(13)
$\displaystyle\leq\|f^{*}-\hat{V}_{k}\|_{p,\hat{\mu}}+2\epsilon_{0}^{\prime}$
(14)
$\displaystyle\leq\left\|f^{*}-BV_{k}\right\|_{p,\hat{\mu}}+3\epsilon_{0}^{\prime}$
(15)
$\displaystyle\leq\left\|f^{*}-BV_{k}\right\|_{p,\mu}+4\epsilon_{0}^{\prime}$
(16) $\displaystyle\leq
d_{p,\mu}(B\mathcal{F},\mathcal{F})+4\epsilon_{0}^{\prime}$

Since $V_{k+1}$ is the best fit to $\hat{V}_{k}$ in $\mathcal{F}$,
$(\ref{single iter inequality 3})$ holds w.p. 1. So we only need to prove that
$(\ref{single iter inequality 1})$, $(\ref{single iter inequality 2})$,
$(\ref{single iter inequality 4})$, $(\ref{single iter inequality 5})$ hold
w.p. at least $1-\delta_{0}^{\prime}$ where
$\delta_{0}^{\prime}=\delta_{0}/4$. The proof of $(\ref{single iter inequality
1})$ and $(\ref{single iter inequality 5})$ is the same as that in (Munos and
Szepesvári [18], Lemma1), and we provide a quick proof here for completeness.
Let

$Q=\max\left(\left|\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}-\left\|V_{k+1}-BV_{k}\right\|_{p,\hat{\mu}}\right|,\left|\left\|f^{*}-BV_{k}\right\|_{p,\mu}-\left\|f^{*}-BV_{k}\right\|_{p,\hat{\mu}}\right|\right).$

So $(\ref{single iter inequality 1})$ and $(\ref{single iter inequality 5})$
will follow if

$\mathbb{P}\left(Q>\epsilon_{0}^{\prime}\right)\leq\delta_{0}^{\prime}.$ (17)

And due to the inequality:

$Q\leq\sup_{f\in\mathcal{F}}\Big{|}\|f-BV_{k}\|_{p,\mu}-\|f-BV_{k}\|_{p,\hat{\mu}}\Big{|},$
(18)

we have:

$\mathbb{P}\left(Q>\epsilon_{0}^{\prime}\right)\leq\mathbb{P}\left(\sup_{f\in\mathcal{F}}\Big{|}\|f-BV_{k}\|_{p,\mu}-\|f-BV_{k}\|_{p,\hat{\mu}}\Big{|}>\epsilon_{0}^{\prime}\right).$
(19)

For any event $\omega$ such that

$\sup_{f\in\mathcal{F}}\Big{|}\|f-BV_{k}\|_{p,\mu}-\|f-BV_{k}\|_{p,\hat{\mu}}\Big{|}>\epsilon_{0}^{\prime},$
(20)

there exist a function $f^{\prime}\in\mathcal{F}$ such that

$\Big{|}\|f^{\prime}-BV_{k}\|_{p,\mu}-\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}\Big{|}>\epsilon_{0}^{\prime}.$
(21)

First assume that
$\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}\leq\|f^{\prime}-BV_{k}\|_{p,\mu}$. Hence,
$\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}+\epsilon_{0}^{\prime}\leq\|f^{\prime}-BV_{k}\|_{p,\mu}$.
Since the elementary inequality $x^{p}+y^{p}\leq(x+y)^{p}$ holds for $p\geq 1$
and any non-negative x, y, we can get
$\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}^{p}+(\epsilon_{0}^{\prime})^{p}\leq(\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}+\epsilon_{0}^{\prime})^{p}\leq\|f^{\prime}-BV_{k}\|_{p,\mu}^{p}$.
And thus

$\left|\left\|f^{\prime}-BV_{k}\right\|_{p,{\mu}}^{p}-\left\|f^{\prime}-BV_{k}\right\|_{p,\hat{\mu}}^{p}\right|>(\epsilon_{0}^{\prime})^{p}.$
(22)

This inequality holds for an analogous reason when
$\|f^{\prime}-BV_{k}\|_{p,\hat{\mu}}>\|f^{\prime}-BV_{k}\|_{p,\mu}$. And since

$\sup_{f\in\mathcal{F}}\left|\|f-BV_{k}\|_{p,\mu}^{p}-\|f-BV_{k}\|_{p,\hat{\mu}}^{p}\right|\geq\left|\left\|f^{\prime}-BV_{k}\right\|_{p,\mu}^{p}-\left\|f^{\prime}-BV_{k}\right\|_{p,\hat{\mu}}^{p}\right|,$
(23)

we can get

$\mathbb{P}\left(\sup_{f\in\mathcal{F}}\Big{|}\|f\\!-\\!BV_{k}\|_{p,\mu}\\!-\\!\|f\\!-\\!BV_{k}\|_{p,\hat{\mu}}\Big{|}\\!>\\!\epsilon_{0}^{\prime}\right)\\!\leq\\!\mathbb{P}\left(\sup_{f\in\mathcal{F}}\left|\|f\\!-\\!BV_{k}\|_{p,\mu}^{p}\\!-\\!\|f\\!-\\!BV_{k}\|_{p,\hat{\mu}}^{p}\right|\\!>\\!\left(\epsilon_{0}^{\prime}\right)^{p}\right).$
(24)

Observe that
$\|f-BV_{k}\|_{p,\mu}^{p}=\mathbb{E}\left[\left|\left(f\left(s_{1}\right)-BV_{k}\left(s_{1}\right)\right)\right|^{p}\right]$,
and $\|f-BV_{k}\|_{p,\hat{\mu}}^{p}$ is just the sample average approximation
of $\|f-BV_{k}\|_{p,\mu}^{p}$. Using Lemma A.1, we can get

$\displaystyle\mathbb{P}\left(\sup_{f\in\mathcal{F}}\left|\|f\\!-\\!BV_{k}\|_{p,\mu}^{p}\\!-\\!\|f\\!-\\!BV_{k}\|_{p,\hat{\mu}}^{p}\right|\\!>\\!\left(\epsilon_{0}^{\prime}\right)^{p}\right)$
$\displaystyle\leq\\!8\mathbb{E}\\!\left[\mathcal{N}\Big{(}\frac{\left(\epsilon_{0}^{\prime}\right)^{p}}{8},\mathcal{F}\left(X^{1:N}\right)\\!\Big{)}\\!\right]e^{-\frac{N}{2}\\!\left(\\!\frac{1}{8}\\!\big{(}\frac{\epsilon_{0}^{\prime}}{2V_{max}}\\!\big{)}^{p}\right)^{2}}$
(25)
$\displaystyle=8\mathcal{N}_{0}(N)e^{-\frac{N}{2}\left(\frac{1}{8}\big{(}\frac{\epsilon_{0}^{\prime}}{2V_{max}}\big{)}^{p}\right)^{2}}.$

Making the right-hand side upper bounded by $\delta_{0}^{\prime}=\delta_{0}/4$
yields a lower bound of $N$, displayed in (10). Then (17) can be proved by
combining (19), (24), and (25).

Now we turn to prove (13) and (15). For an arbitrary function
$f\in\mathcal{F}$, using the triangle inequality, we have

$\left|\|f-BV_{k}\|_{p,\hat{\mu}}-\|f-\hat{V}_{k}\|_{p,\hat{\mu}}\right|\leq\|BV_{k}-\hat{V}_{k}\|_{p,\hat{\mu}}.$
(26)

So if we show that
$\|BV_{k}-\hat{V}_{k}\|_{p,\hat{\mu}}\leq\epsilon_{0}^{\prime}$ holds w.p.
$1-\delta_{0}^{\prime}$, then (13) and (15) can be proved by choosing
$f=V_{k+1}$ and $f=f^{*}$, respectively. To prove
$\|BV_{k}-\hat{V}_{k}\|_{p,\hat{\mu}}\leq\epsilon_{0}^{\prime}$, Recall that
the data $(s,a,\hat{s^{\prime}},\hat{r})$ used to calculate $\hat{V}_{k}$ is
sampled from real environment w.p. $\beta$ (case 1), or from the learned model
w.p. $1-\beta$ (case 2).

For case 1, it satisfies that $\hat{s^{\prime}}=s^{\prime}$ and $\hat{r}=r$.
So we have

$\mathbb{P}\left(|\hat{r}+\gamma V_{k}(\hat{s^{\prime}})-r-\gamma
V_{k}(s^{\prime})|>\epsilon_{0}^{\prime}\right)=0.$ (27)

And for case 2, like in previous literature [3, 17], assume the reward
function is known, i.e., $\hat{r}=r$. We have

$\displaystyle\mathbb{P}\left(|\hat{r}+\gamma V_{k}(\hat{s^{\prime}})-r-\gamma
V_{k}(s^{\prime})|>\epsilon_{0}^{\prime}\right)$
$\displaystyle=\mathbb{P}\left(\gamma|(V_{k}(\hat{s^{\prime}})-V_{k}(s^{\prime}))|>\epsilon_{0}^{\prime}\right)$
$\displaystyle\leq\mathbb{P}\left(\gamma
K_{V}\|(\hat{s^{\prime}}-s^{\prime})\|_{2}>\epsilon_{0}^{\prime}\right)$
$\displaystyle=\mathbb{P}\left(\|\hat{s^{\prime}}-s^{\prime}\|_{2}>\frac{\epsilon_{0}^{\prime}}{\gamma
K_{V}}\right).$ (28)

By using the model error assumption:
$f(\|\hat{s^{\prime}}-s^{\prime}\|_{2})=\frac{2}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\|\hat{s^{\prime}}-s^{\prime}\|_{2}^{2}}{2\sigma^{2}}\right)$,
we can write that:

$\displaystyle\mathbb{P}\left(|\hat{r}+\gamma V_{k}(\hat{s^{\prime}})-r-\gamma
V_{k}(s^{\prime})|>\epsilon_{0}^{\prime}\right)$
$\displaystyle\leq\mathbb{P}\left(\|\hat{s^{\prime}}-s^{\prime}\|_{2}>\frac{\epsilon_{0}^{\prime}}{\gamma
K_{V}}\right)$
$\displaystyle=\mathbb{P}\left(\frac{\|\hat{s^{\prime}}-s^{\prime}\|}{\sigma}>\frac{\epsilon_{0}^{\prime}}{\gamma
K_{V}\sigma}\right)$
$\displaystyle=2\Big{(}1-\Phi(\frac{\epsilon_{0}^{\prime}}{\gamma
K_{V}\sigma})\Big{)},$ (29)

where $\Phi$ is the Cumulative Distribution Function (CDF) of the standard
normal distribution. Combine these two cases, we can get

$\mathbb{P}\left(\Big{|}\hat{r}+\gamma V_{k}(\hat{s^{\prime}})-r-\gamma
V_{k}(s^{\prime})\Big{|}>\epsilon_{0}^{\prime}\right)\leq
2(1-\beta)\Big{(}1-\Phi(\frac{\epsilon_{0}^{\prime}}{\gamma
K_{V}\sigma})\Big{)}.$ (30)

Making the right-hand side upper bounded by
$\delta_{0}^{\prime}/(N|\mathcal{A}|)$ yields an upper bound of $\sigma$,
displayed in (11). And for each sampled state $s_{i}$, since

$\left|BV_{k}\left(s_{i}\right)-\hat{V}_{k}\left(s_{i}\right)\right|\leq\max_{a\in\mathcal{A}}\Big{|}\hat{r}+\gamma
V_{k}(\hat{s^{\prime}})-r-\gamma V_{k}(s^{\prime})\Big{|},$

by using a union bounding argument, we can get

$\mathbb{P}\left(\left|BV_{k}\left(s_{i}\right)-\hat{V}_{k}\left(s_{i}\right)\right|>\epsilon_{0}^{\prime}\right)\leq\delta_{0}^{\prime}/N.$
(31)

And by another union bounding argument, we have

$\mathbb{P}\left(\max_{i=1,\ldots,N}\left|BV_{k}\left(s_{i}\right)-\hat{V}_{k}\left(s_{i}\right)\right|>\epsilon_{0}^{\prime}\right)\leq\delta_{0}^{\prime}.$
(32)

And therefore,

$\mathbb{P}\left(\frac{1}{N}\left|BV_{k}\left(s_{i}\right)-\hat{V}_{k}\left(s_{i}\right)\right|>\epsilon_{0}^{\prime}\right)\leq\delta_{0}^{\prime}.$
(33)

Hence, we have proved that
$\|BV_{k}-\hat{V}_{k}\|_{p,\hat{\mu}}\leq\epsilon_{0}^{\prime}$ holds w.p. at
least $1-\delta_{0}^{\prime}$. This completes the whole proof. ∎

######  Lemma A.3.

(Munos and Szepesvári [18], Theorem 2) Under the concentrability assumption
(Assumption A.1), let $\rho$ be an arbitrary state distribution, and
$C_{\rho,\mu}$ be the discounted-average concentrability coefficient. If each
iteration error can be bounded as $\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}\leq
d_{p,\mu}(B\mathcal{F},\mathcal{F})+(1-\gamma)^{2}\epsilon/(4\gamma
C_{\rho,\mu}^{1/p})$ w.p. at least $1-\delta/K$ for $0\leq k<K$, then the loss
due to using $\pi_{K}$ instead of the optimal policy $\pi^{*}$ satisfies that
w.p. at least $1-\delta$,

$\left\|V^{*}-V^{\pi_{K}}\right\|_{p,\rho}\leq\frac{2\gamma}{(1-\gamma)^{2}}C_{\rho,\mu}^{1/p}d_{p,\mu}(B\mathcal{F},\mathcal{F})+\epsilon,$
(34)

provided that

$\gamma^{K}<\left[\frac{(1-\gamma)^{2}}{8\gamma V_{\max}}\epsilon\right]^{p}.$
(35)

Lemma A.3 shows that if the error in each iteration can be bounded with high
probability, then the return discrepancy between $\pi_{K}$ and $\pi^{*}$ can
be bounded with high probability when the number of iterations $K$ is large
enough. Moreover, we can set $\rho$ to the distribution of the states we care
about more, e.g., the initial states we start to execute $\pi_{K}$.

######  Theorem A.1.

($\beta$-mixture sampling-based FVI bound) Assume that the pseudo-dimension of
the value function space is finite as in Proposition A.1, under the
concentrability assumption (Assumption A.1) and the same assumptions of Lemma
4.1, let $\rho$ be an arbitrary state distribution, $C_{\rho,\mu}$ be the
discounted-average concentrability coefficient and $\pi_{K}$ be the greedy
policy w.r.t. $V_{K}$. Let $V^{\pi_{K}}$ and $V^{*}$ be the expected return of
executing $\pi_{K}$ and the optimal policy $\pi^{*}$ in real environment,
respectively. Define $N_{\rm{real}}=N\cdot|\mathcal{A}|\cdot\beta$ as the
expected number of real samples. Then the following bound holds w.p. at least
$1-\delta$:

$\displaystyle\|V^{*}-V^{\pi_{K}}\|_{p,\rho}\leq$
$\displaystyle\frac{2\gamma}{(1-\gamma)^{2}}C_{\rho,\mu}^{1/p}d_{p,\mu}(B\mathcal{F},\mathcal{F})+O\left(\Big{(}\frac{\beta|\mathcal{A}|}{N_{\rm{real}}}\big{(}\log(\frac{N_{\rm{real}}}{\beta|\mathcal{A}|})+\log(\frac{K}{\delta})\big{)}\Big{)}^{\frac{1}{2p}}\right)$
(36)
$\displaystyle+O\Big{(}\Phi^{-1}\big{(}1-\frac{\beta\delta}{8KN_{\rm{real}}(1-\beta)}\big{)}\sigma\Big{)}+O\left(\gamma^{K/p}V_{\max}\right).$

###### Proof.

Lemma A.3 implies that we can bound
$\left\|V^{*}-V^{\pi_{K}}\right\|_{p,\rho}$ w.p. at least $1-\delta$ as

$\left\|V^{*}-V^{\pi_{K}}\right\|_{p,\rho}\leq\frac{2\gamma}{(1-\gamma)^{2}}C_{\rho,\mu}^{1/p}d_{p,\mu}(B\mathcal{F},\mathcal{F})+\epsilon,$
(37)

via bounding $\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}$ w.p. at least
$1-\delta/K$ as

$\left\|V_{k+1}-BV_{k}\right\|_{p,\mu}\leq
d_{p,\mu}(B\mathcal{F},\mathcal{F})+(1-\gamma)^{2}\epsilon/(4\gamma
C_{\rho,\mu}^{1/p}).$ (38)

Using Lemma A.2 to bound (38) by setting
$\epsilon_{0}=(1-\gamma)^{2}\epsilon/(4\gamma C_{\rho,\mu}^{1/p})$ and
$\delta_{0}=\delta/K$. We can get the corresponding bounds of $N$ and $\sigma$
w.r.t. $\epsilon$. Combined with (35), we can write $\epsilon$ as

$\epsilon=O\left(\Big{(}\frac{1}{N}\big{(}\log(\mathcal{N}_{0}(N))+\log(\frac{K}{\delta})\big{)}\Big{)}^{\frac{1}{2p}}\right)+O\Big{(}\Phi^{-1}\big{(}1-\frac{\delta}{8KN|\mathcal{A}|(1-\beta)}\big{)}\sigma\Big{)}+O\left(\gamma^{K/p}V_{\max}\right).$
(39)

Then using Proposition A.1 to bound $\mathcal{N}_{0}(N)$, we can get:

$\displaystyle\epsilon$ $\displaystyle\leq
O\left(\Big{(}\frac{1}{N}\big{(}\log(1/\epsilon_{0})+\log(\frac{K}{\delta})\big{)}\Big{)}^{\frac{1}{2p}}\right)+O\Big{(}\Phi^{-1}\big{(}1-\frac{\delta}{8KN|\mathcal{A}|(1-\beta)}\big{)}\sigma\Big{)}+O\left(\gamma^{K/p}V_{\max}\right)$
(40) $\displaystyle\leq
O\left(\Big{(}\frac{1}{N}\big{(}\log(N)+\log(\frac{K}{\delta})\big{)}\Big{)}^{\frac{1}{2p}}\right)+O\Big{(}\Phi^{-1}\big{(}1-\frac{\delta}{8KN|\mathcal{A}|(1-\beta)}\big{)}\sigma\Big{)}+O\left(\gamma^{K/p}V_{\max}\right),$

where the last inequality comes from inequality 10. Plugging (40) into (37)
and substituting $N$ with $N_{\rm{real}}/(\beta|\mathcal{A}|)$ complete the
proof. ∎

##  Appendix B Model Error Assumption

In Lemma 4.1 and Theorem 4.1, we assume that the model error between the real
next state $s^{\prime}$ and the predicted ones $\hat{s^{\prime}}$ obeys a
half-normal distribution, i.e., the probability density function
$f(\|\hat{s^{\prime}}-s^{\prime}\|_{2})=\frac{2}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\|\hat{s^{\prime}}-s^{\prime}\|_{2}^{2}}{2\sigma^{2}}\right)$.
In this section, we give an empirical analysis of this model error assumption.
To be more specific, we test the trained dynamics model in MBPO on newly
collected data (typically 10k transitions) and plot the frequency of different
model error ranges in Figure 5 to approximate the distribution. We can observe
that the approximated distribution is close to the half-normal distribution,
i.e., the probability of small model error is higher than that of large model
error, which empirically shows the rationality of the assumption.

Figure 5: Frequency distribution histogram of the model error
$\|\hat{s^{\prime}}-s^{\prime}\|_{2}$. The y-axis indicates the frequency of
the corresponding range.

##  Appendix C Algorithm Comparison

We compare the hyperparameters that can be scheduled in different Dyna-style
MBRL algorithms in Table 1. From the comparison, real ratio and rollout length
can only be scheduled in MBPO and MA-BDDPG since other algorithms merely use
the imaginary data to train the policy and need to generate model rollouts
from the initial state to the end of the episode. Considering that MBPO is
much more effective than MA-BDDPG in continuous control benchmark tasks [12,
13], we finally choose MBPO as a representative running case.

Table 1: Comparison of the hyperparameters that can be scheduled in different Dyna-style MBRL algorithms. | Real Ratio |  | Policy Training
---
Iteration
| Model Training
---
Frequency
Rollout Length
MBPO [12] | ✓ | ✓ | ✓ | ✓
MA-BDDPG [13] | ✓ | ✓ | ✓ | ✓
SLBO [17] | $\times$ | ✓ | ✓ | $\times$
ME-TRPO [14] | $\times$ | ✓ | ✓ | $\times$
MB-MPO [4] | $\times$ | ✓ | ✓ | $\times$
PAL/MAL [22] | $\times$ | ✓ | ✓ | $\times$

##  Appendix D More Experimental Results

###  D.1 Hyperparameter Importance

Figure 6: Complete figures of the hyperparameter importance experiments in
Section 6.3.

Due to the page limit, we only show the results of hyperparameter importance
on three environments in Section 6.3. For a more comprehensive analysis, we
plot the results on all the six environments in Figure 6. The conclusion is
the same as discussed in Section 6.3, i.e., using the hyper-controller to
schedule the real ratio retains much of the advantage of AutoMBPO.

###  D.2 Controller Transfer

From the results in Figure 3, we observe that the hyperparameter schedules on
the three PyBullet environments are similar, especially for the real ratio
$\beta$. Then we further conduct experiments to test whether the hyper-
controller trained on these three tasks can transfer to others without
additional fine-tuning. Results are shown in Figure 7. Though the performance
of the transferred hyper-controller is slightly worse than the original ones,
they all surpass the MBPO baseline with a considerable margin, which shows the
potential of the hyper-controller to capture the commonality to generalize to
different tasks.

Figure 7: The transfer results of the hyper-controller among PyBullet
environments with $A_{3}^{2}=6$ transferring cases in total. Each figure’s
title represents the target task, and each specific line indicates the source
task. For example, the green line in the middle figure represents transferring
the hyper-controller learned on HopperBullet to Walker2dBullet.

###  D.3 Effectiveness of AutoMBPO

In Algorithm 5.1, we train the hyper-controller until it achieves acceptable
performance. One concern is whether the hyper-controller improves
monotonically during training. We plot the performance of the hyper-controller
in different training phases in Figure 8. For example, suppose we train the
hyper-controller with 200 hyper-MDP episodes, i.e., 200 MBPO instances in
total, then the purple line ($1\%$-$20\%$) represents the average return of
the $1$-$40$th MBPO instances, and so on. Monotonic improvement of hyper-
controller can be observed during its training phase, which further
demonstrates the effectiveness of our algorithm.

Another concern may be the extensive computational requirements for AutoMBPO.
According to the computational time table in Appendix G, take Hopper as an
example: our algorithm takes about 90 hours to train the hyper-controller,
about 5-6 times the time to train a complete MBPO instance, i.e., equivalent
to 5-6 trials of MBPO hyperparameters. It is almost impossible to manually
find a suitable configuration for these hyperparameters within such few
trials.

Figure 8: Performance of the hyper-controller in different training phases.
Results are the average on six trials of hyper-controller training over
different random seeds.

###  D.4 Hyperparameter Study

Figure 9: Top: Performance comparison of AutoMBPO(0.5, 10, 20), AutoMBPO(0.05,
20, 40) and the original AutoMBPO(0.05, 10, 20). Middle and bottom: the
corresponding schedule of real ratio and policy training iteration.

Since AutoMBPO utilizes the hyper-controller to adjust the hyperparameters of
MBPO dynamically, e.g., $\pm 1$, the initial values of these hyperparameters
become the hyperparameters of AutoMBPO. We also want to investigate the effect
of these initial values on the learning process. In the experiments of Section
6, we initialize the real ratio to $0.05$, the policy training iteration to
$10$, and limit the policy training iteration to $[1,20]$ for computation
efficiency. We denote the original one as AutoMBPO(0.05, 10, 20). In this
section, we further conduct experiments of AutoMBPO(0.5, 10, 20) and
AutoMBPO(0.05, 20, 40) on three PyBullet environments. The performance and the
hyperparameter schedules are shown in Figure 9.

From the comparison, we can find that changing the initial value of real ratio
or policy training iteration does not influence the final performance much,
which shows the robustness of AutoMBPO to the hyperparameter initialization.
Notice that the policy training iteration of AutoMBPO(0.05, 20, 40) is much
larger than that of AutoMBPO(0.05, 10, 20), while the performance does not
improve much. This finding is consistent with our conclusion in Section 6.2,
i.e., increasing the policy training iteration is not always a good choice.

###  D.5 State Feature Ablation

In section 5.1, inspired by Theorem 4.1, we include the number of real samples
$N_{\rm{real}}$ and the model loss $\mathcal{L}_{\hat{T}}$ into the state
formulation of hyper-MDP. In this section, we ablate these two features and
only use other features as the state to train the controller, denoted as
AutoMBPO-SA. Results are shown in Figure 10. We find that ablating these two
features degrades the performance, which highlights the importance of
elaborate state design for hyper-controller learning.

Figure 10: Results of the state feature ablation experiment. AutoMBPO-SA
denotes the AutoMBPO variant of excluding real samples number and model loss
from the hyper-MDP state.

###  D.6 Statistical Significance

Since the shaded areas in Figure 2 overlap, we further add a t-test to the
results of the original MBPO and AutoMBPO. We use the average return to
perform t-test and list the p-values in Table 2. From the result, the p-values
are all less than 0.05, confirming the statistical significance of the
results.

Table 2: t-test to the average returns of the original MBPO and AutoMBPO. | Hopper | Ant | Humanoid |  | Hopper
---
Bullet
| Walker2d
---
Bullet
| HalfCheetah
---
Bullet
p-value | 0.0078 | 0.0031 | 6e-5 | 0.0052 | 9e-5 | 0.0276

##  Appendix E Baseline Implementation Details

PBT. The implementation of PBT mainly follows Zhang et al. [32]. Specifically,
we train 10 MBPO instances in parallel, each with randomly initialized
hyperparameters. After each episode, 20$\%$ of the MBPO instances with low
returns are replaced by the top 20$\%$ (both hyperparameters, network
parameters, and data buffer) or re-initialized (only hyperparameters) with a
certain probability.

RoR. For RoR, we utilize the same Hyper-MDP definition as AutoMBPO since the
original Hyper-MDP definition in Dong et al. [8] is not suitable for MBPO.
Moreover, we run RoR 6 times, each using the same amount of data as our
method. In the result of Figure 2, we compare to the best performing RoR for a
fair comparison.

##  Appendix F Experimental Settings

We first provide the details of the experimental environments in Table 3.
Among them, Hopper, Ant, and Humanoid are the same version used in Janner et
al. [12]. For the remaining three PyBullet environments, we modify the reward
function to be similar to that of Mujoco to facilitate training since we found
that the default reward setting is too complicated for MBPO to solve in
preliminary experiments. Then, we present the experimental settings in
different environments in Table 4.

Table 3: Environment settings in our experiments. $\theta_{t}$ denotes the joint angle, $x_{t}$ denotes the position in x-direction, $a_{t}$ denotes the action control input, and $z_{t}$ denotes the height. |  | State
---
Dimension
| Action
---
Dimension
Reward Function |  | Termination States
---
Condition
Hopper | 11 | 3 | $\dot{x}_{t}-0.001\left\|a_{t}\right\|_{2}^{2}+1$ |  $z_{t}\leq 0.7$ or $\theta_{t}\geq 0.2$
Ant | 27 | 8 | $\dot{x}_{t}-0.5\left\|a_{t}\right\|_{2}^{2}+1$ |  $z_{t}\leq 0.2$ or $z_{t}\geq 1.0$
Humanoid | 45 | 17 | 0.25$\dot{x}_{t}-0.1\left\|a_{t}\right\|_{2}^{2}+5$ |  $z_{t}\leq 1.0$ or $z_{t}\geq 2.0$
| Hopper
---
Bullet
15 | 3 | 5$\dot{x}_{t}-0.001\left\|a_{t}\right\|_{2}^{2}+1$ |  $z_{t}\leq 0.8$ or $|\theta_{t}|\geq 1.0$
| Walker2d
---
Bullet
22 | 6 | 5$\dot{x}_{t}-0.001\left\|a_{t}\right\|_{2}^{2}+1$ |  $z_{t}\leq 0.8$ or $|\theta_{t}|\geq 1.0$
| HalfCheetah
---
Bullet
26 | 6 | 5$\dot{x}_{t}-0.001\left\|a_{t}\right\|_{2}^{2}$ | None
Table 4: Experimental settings in different environments. Specifically, a hyper-MDP episode consists of m target-MDP episodes, i.e., the whole training process of an MBPO instance, and a target-MDP episode consists of $H$ timesteps in the environments. |  | Hopper | Ant | Humanoid |  | Hopper
---
Bullet
| Walker2d
---
Bullet
| HalfCheetah
---
Bullet
|  | hyper-MDP episodes
---
100 | 100 | 200 | 200 | 200 | 100
$m$ |  | target-MDP epi-
---
sodes for training
50 | 50 | 100 | 30 | 50 | 20
$M$ |  | target-MDP epi-
---
sodes for evaluation
100 | 150 | 300 | 100 | 100 | 50
$H$ |  | timesteps per
---
target-MDP episode
1000

##  Appendix G Computing Infrastructure

We present the computing infrastructure and the corresponding computational
time used to train the hyper-controller in Table 5.

Table 5: Computing infrastructure and the corresponding computational time. | Hopper | Ant | Humanoid |  | Hopper
---
Bullet
| Walker2d
---
Bullet
| HalfCheetah
---
Bullet
CPU | 32 cores | 16 cores
GPU | RTX2080TI$\times$2 | V100 $\times$2
| computation
---
time in hours
90.92 | 95.05 | 245.33 | 56.91 | 149.88 | 30.48

##  Appendix H Hyperparameters

Table 6 lists the hyperparameters used in training the hyper-controller. Other
hyperparameters of MBPO not scheduled by the hyper-controller are the same as
the original one [12]. Note that the hyperparameter $\tau$ in Humanoid varies
from that in other environments since the original MBPO configuration of
Humanoid is different. The original MBPO trains the model per 1000 real
timesteps in Humanoid but per 250 real timesteps in other environments. So we
set $\tau$ to half of the original interval to keep the maximum model training
frequency two times of the original configuration.

Table 6: Hyperparameter settings for hyper-controller. |  | Hopper | Ant | Humanoid |  | Hopper
---
Bullet
| Walker2d
---
Bullet
| HalfCheetah
---
Bullet
$\tau$ |  | real timesteps
---
interval per action
125 | 500 | 125
|  | policy network
---
architecture
MLP with one hidden layer of size 256
| learning rate | $3\cdot 10^{-4}$
| batch size | 64
|  | policy updates per
---
hyper-MDP episode
30
|  | initial value
---
of real ratio
0.05
|  | initial value of poli-
---
cy training iteration
10
|  | initial value of
---
rollout length
1
$\epsilon$ |  | PPO clip
---
constant
0.2
$c$ |  | real ratio
---
change constant
1.2
|  | penalty for each
---
model training
0.1


# Skeletal Point Representations with Geometric Deep Learning

###### Abstract

Skeletonization has been a popular shape analysis technique that models both
the interior and exterior of an object. Existing template-based calculations
of skeletal models from anatomical structures are a time-consuming manual
process. Recently, learning-based methods have been used to extract skeletons
from 3D shapes. In this work, we propose novel additional geometric terms for
calculating skeletal structures of objects. The results are similar to
traditional fitted s-reps but but are produced much more quickly. Evaluation
on real clinical data shows that the learned model predicts accurate skeletal
representations and shows the impact of proposed geometric losses along with
using s-reps as weak supervision.

Index Terms—  Geometric learning, Skeletal representations, Shape analysis,
Point clouds.

##  1 Introduction

Skeletonization has been a popular approach for modeling anatomical structures
due to the ability of skeletal representations to capture both the boundary
and interior of an object. This is in contrast with simpler models such as
densely-sampled boundary landmarks. The most widely used approach to define an
object’s skeleton is through a medial axis transform [1] (MAT)(fig. 1, left).
The MAT consists of a set of points and associated radii, called spokes, which
form the set of maximally inscribed spheres inside the shape. MAT-based models
have been used for a wide range of applications [2] such as segmentation,
registration, finite element modeling, and statistics of object shape.

The main limitation of the MAT is that it has the tendency to amplify small-
scale noise on an object’s boundary into large deviations in skeleton location
and topology across a population. This makes it unsuited for uses which
require consistent structure such as statistical analysis.

This limitation has led to modifications [3, 4, 5] to the MAT to address these
shortcomings. In particular, s-reps [6](fig. 1, right) are a class of discrete
skeletal representations related to the MAT which have a fixed topology and
consistent sampling across a population. This is done by fitting a template
s-rep to an object via an optimization [7] rather than direct computation from
the object’s boundary. Having a fixed template which is then adapted to each
individual in a population yields improved consistency and correspondence.

The optimization is constrained in such a way as to have the final object be
nearly medial in that points on the skeleton are approximately equidistant
from the top and bottom of the object’s surface and the radii emanating from
these points are nearly orthogonal to the boundary. This allows for s-reps to
have the stability required for statistical analysis while still taking
advantage of some of the beneficial properties of a purely medial skeleton. A
downside to this approach is that the optimization can be slow and often
requires manual template generation and parameter tuning when applied to a new
dataset.

Learning-based methods for skeletonizing images and shapes are a relatively
recent line of research which have shown promise in robustly computing
skeletal representations. Earlier learning-based methods primarily focused on
the task of extracting skeletons from 2D natural images [8]. They typically
involved use of fully convolutional architectures and pose the problem in a
similar fashion to that of image segmentation as seen in works based on the
popular U-Net architecture [9, 10]. In contrast, there is less work on
learning-based skeletonization of 3D objects, partially due to their increased
complexity and partially due to the lack of a benchmark dataset for training.
This led to the development of point-based methods like Point2Skeleton [11]
which utilizes PointNet++ [12] as a point encoder and tries to predict weights
on input points to generate the skeleton as a convex combination of the inputs
in a manner similar to [13].

Fig. 1: (Left) Medial axis for a 2D shape, (Right) s-rep for a hippocampus
surface with yellow lines as the spoke vectors. Fig. 2: Model overview and
geometric terms used for training.

One of the drawbacks of the point based methods is that they are not optimal
for use in shape analysis tasks like s-rep fitting due to inconsistent
sampling density and sometimes large differences in the MAT of similar shapes.
Such deviations are harder to tolerate for applications in shape analysis
where a consistent pattern in the skeletons across a population is necessary.
In this paper we build upon these point-based learning methods by introducing
a novel way to utilize s-reps as a form of supervision (instead of computing
MATs) and introducing two additional geometric loss terms, one for enforcing
the medial approximation of the skeletal points and a second regularization
term to encourage spread between them. The proposed changes show good results
on synthetic data as well as promising application to two clinical problems of
modeling hippocampi and tricuspid valve leaflets.

##  2 Methodology

In this section we present the base model that defines a learning based s-rep
to then motivate the three geometric terms for better skeletal point
prediction.

###  2.1 Base Model

Our basic definition of a learning-based s-rep is built upon the previous
Point2Skeleton [11, 13] method. As seen in fig. 2, the input for our model is
a segmented shape’s surface point cloud. We then use the PointNet++ layers as
encoders and have a shared multilayer perceptron (MLP) across all input points
to predict weights over them. Assuming there are N input points $P$ after the
encoding and M skeletal points $S$ to be predicted, the model will try to
predict a weight matrix $W\in\mathcal{R}^{N\times M}$. The weights are
constrained to form a convex combination over the input points and hence we
obtain the skeletal points $S=W^{T}P$. The radii are predicted in a similar
fashion by first computing the distance $d_{i}=min_{s_{j}\in
S}||p_{i}-s_{j}||_{2},\;\forall p_{i}\in P$. Then the radius for each skeletal
point is predicted using $r(s_{j})=\sum_{i=1}^{N}W_{ij}d_{i}$ using the convex
combination over the input points. The drawback with this approach is that the
model does not learn to well approximate a smooth skeletal sheet due to the
lack of supervision. The predictions are scattered around the interior as seen
on the right of fig. 3.

Fig. 3: Effect of s-rep supervision: White dots are the G.T skeletal points,
black are the predictions with supervision and red are without supervision.

###  2.2 Novel Geometric Losses

####  2.2.1 MAT supervision via s-reps

We propose to use supervision to allow the network to find skeletal points
which better approximate a smooth sheet. Using the MAT for supervision will
yield a network susceptible to producing skeletons with the same undesirable
properties. We instead propose to use pre-computed s-reps as a form of
supervision to allow the model to better approximate actual s-reps. The
supervision with s-reps as ground truth is enforced through the means of a
symmetric Chamfer distance loss ($L_{CD}$) with the predicted skeletal points
as seen in Eqn. 1, where $S1$ and $S2$ are the predicted and ground truth
skeletal sheets.

$\begin{split}L_{CD}(S_{1},S_{2})=\frac{1}{|S_{1}|}\sum_{a\in S_{1}}\min_{b\in
S_{2}}||a-b||_{2}^{2}\\\ +\frac{1}{|S_{2}|}\sum_{b\in S_{2}}\min_{a\in
S_{1}}||b-a||_{2}^{2}\end{split}$ (1)

####  2.2.2 Spread regularization

Another aspect of the base model is that there is no regularization on the
skeletal points in terms of their spread. Ideally we want the skeletal points
to cover the interior of the object as well as possible. Without this, some
parts of the object interior are not well represented which in turn affects
the quality of the surface reconstructed from the skeleton. Thus we impose a
spread regularization term $L_{spread}$ shown in Eqn. 2 to discourage the
clustering of skeletal points by penalizing the average pairwise distances
between them.

$L_{spread}=-\frac{1}{M^{2}}\sum_{i=1}^{N}\sum_{j=1}^{M}||s_{i}-s_{j}||_{2}$
(2)

####  2.2.3 Medial Enforcement

Recall that the skeletal points in a MAT have more than one closest point on
the surface. The base model incorporates this loss through a point to sphere
loss [11] with two terms: one on the input points forcing them to lie on the
closest skeletal surface and another on the skeletal points forcing them to
touch the closest input point. We modify the second to encourage a medial
approximation by having each skeletal point’s sphere touch its nearest three
points following the definition of MAT. This loss $L_{medial}$, shown in Eqn.
3 will help train the network to adhere to a stronger medial approximation.

$L_{medial}=\sum_{s\in S}\sum_{p\in C_{3}(s)}||p-s||_{2}-r(s)$ (3)

Here the set $C_{3}(s)\subset P$ represents the closest three input surface
points to the skeletal point $s$ and $r(s)$ represents the radius for the
skeletal sphere at $s$.

####  2.2.4 Spoke vectors

To ensure a fair comparison with s-reps and obtain a surface reconstruction we
also leveraged a simple way to predict the spoke vectors toward the surface
for each skeletal point from the weight matrix. Intuitively the surface points
closer to the skeletal points have a higher weight and hence hold clues to
which direction the nearest surface location lies for a skeletal point.
Therefore the spoke vectors are predicted to be the direction vectors to the
highest weight input point (for each skeletal point) as shown in Eqn. 4

$spoke(s_{j})=\frac{(p_{i^{*}_{j}}-s_{j})}{||p_{i^{*}_{j}}-s_{j}||_{2}}\;\text{where
}i^{*}_{j}=\operatorname*{arg\,min}_{i=1,\dots,N}W_{ij}$ (4)

Note that all surface reconstruction results utilize the notion of spoke
vectors which is a bit different from earlier works where either (a) uniform
sampling is used on the skeletal sphere, or (b) the merged mesh from all
skeletal spheres is used to reconstruct the surface.

##  3 Experiments

###  3.1 Implementation

We used PyTorch version 1.12 for our model implementation running on CUDA
11.6. The experiments were performed on a machine with a single NVIDIA RTX
A5000 GPU having 24GB vRAM. The Adam optimizer was used for training with a
learning rate of 5e-4 and a batch size of 6. We followed a similar training
regime as that of [11] for enabling the geometric losses on the model. For
both training and inference, 1000 points were randomly sampled from the
normalized input (zero mean, scaled down to unit sphere). We set 100 skeletal
points as the desired output from the model. The s-reps for the hippocampi and
leaflets data were obtained using the skeletal representation module in
SlicerSALT [14], a toolbox of shape analysis software. For the quantitative
evaluation metrics, we utilized the symmetric versions of the Chamfer and
Hausdorff distances between two point sets abbreviated as C.D and H.D
respectively in Table 1 and Table 2.

In our implementation, predicting an s-rep takes approximately 1.7 seconds.
The time needed for the s-rep optimization in SlicerSALT varies depending on
how much refinement of the initial s-rep fit is needed. Typical runs take
between 5 and 30 minutes but runs over an hour occasionally occur.

###  3.2 Synthetic Dataset

As stated before, there is no standard 3D dataset for the task of skeletal
point prediction. Taking inspiration from current s-rep fitting approaches
which rely on deforming ellipsoids to match an object, we generated a
synthetic dataset of deformed ellipsoids along with their s-reps. Starting
from a base ellipsoid and s-rep which are axis aligned, we first apply random
scale factors sampled from the normal distribution $\mathcal{N}(1,0.15)$ to
each axis independently. We then deform the ellipsoid by bending the long axis
by angles sampled from $\mathcal{N}(\frac{\pi}{3},\frac{\pi}{8})$ and twisting
by angles sampled from $\mathcal{N}(\frac{\pi}{6},\frac{\pi}{8})$. Fig. 4
shows an example of an ellipsoid pre and post deformation.

Fig. 4:  (Top) The left image shows an ellipsoid (black) and a version
deformed into a leaflet (white) while the right image shows the ground truth
skeletal points (blue sheet). (Bottom) Sample images (front/side view) of
leaflets and hippocampi.

###  3.3 Ablation Study

We conducted an ablation study of the proposed geometric terms on 5000
deformed ellipsoids generated as described in 3.2. Table 1 shows the
performance of the individual and combined effects on held out data. We refer
to using the base model with the supervised loss as the baseline method. The
ablation experiments evaluated two aspects of the model, the first measuring
how close the predicted skeletal points match the s-rep skeletal points and
the second measuring the effectiveness of surface reconstruction using the
spoke vectors.

Table 1: Ablation study results on different configurations. Loss Terms | G.T. Skel Points | Surf Recon
---|---|---
Spread | Medial | C.D | H.D | C.D | H.D
✗ | ✗ | 0.0085 | 0.3528 | 0.0395 | 0.4122
✓ | ✗ | 0.0089 | 0.2082 | 0.0297 | 0.4157
✗ | ✓ | 0.0097 | 0.3527 | 0.0385 | 0.4396
✓ | ✓ | 0.0073 | 0.1903 | 0.0247 | 0.4366

As we can see in table 1, the combined effect of both the terms tends to give
the best result. The effect of the spread regularization is more noticeable
than that of the medial term, perhaps due to the fact that the base model also
enforces the medial terms albeit less strongly. The surface reconstruction
results with the spoke vectors under the H.D metric are weaker than those of
the skeleton due to limitations with predicted spoke vector orientation.

###  3.4 Clinical Data Evaluation

We evaluated the model on shape data from hippocampi and heart valve leaflets
with boundaries computed via SPHARM-PDM. Since the real world shape data for
such anatomical structures is quite limited, we utilized a fine-tuning
strategy. First we trained the model using the synthetic deformed ellipsoid
data to enable learning on a wide variety of shapes. Using this pre-trained
model, we fine-tuned using the limited amount of hippocampi or leaflet data.
Similar to the ablation study, we compared the skeletal point predictions with
the s-rep skeleton points along with the surface reconstruction quality.
Additionally, we also evaluated how similar is the surface reconstruction to
that derived from the s-rep.

Table 2: Results on Hippocampi and Leaflet data Data | Skel. Points | Surf Recon | Srep Recon
---|---|---|---
| C.D | H.D | C.D | H.D | C.D | H.D
Hipp. | 0.004 | 0.097 | 0.014 | 0.238 | 0.023 | 0.231
Leaf. | 0.007 | 0.146 | 0.025 | 0.097 | 0.044 | 0.275

Table. 2 summarizes the results on the held out samples for each dataset. As
we can see, the predicted skeletons are close to what we would expect after
the s-rep fitting process while also being able to reconstruct the surface
fairly well even though the spoke vectors are not refined. Fig. 5 shows
qualitative results of the skeletal point prediction on the two datasets. The
figure shows well spread out skeletal points along with them lying a sheet-
like surface.

###  3.5 Discussion

The experiments shows promising results in obtaining skeletal representations
for slabular structures. Due to its learning-based nature, its runtime has an
advantage over the s-reps fitting process. It is important to note that there
is a trade-off since s-reps may be more accurate for surface reconstruction
due to principled computation of spoke vectors. This shows one of the
limitations of point based models since we do not account for the surface
normal information, which could help align the spoke vectors to the surface.
Adapting the model to work with surface meshes may also refine the results
after taking into account other geometric features like curvature.

Fig. 5: Qualitative results on the hippocampi (left) and leaflets (right)
showing front and side views for both. The white dots are the predicted
skeletal points.

A limitation of the point based models is the requirement to specify the
number of output skeletal points while training. This could result in the the
skeletal surface having a lesser density of predicted points. Spline-based
interpolation methods could help here to obtain a dense and continuous
representation for the skeleton. Finally, in the context of real world shape
analysis applications, it is important to note that all skeletonization
methods assume a segmented boundary surface. For anatomical structures, it is
likely that the surface is the result of an algorithm that will be uncertain
on certain parts of the surface. Therefore, future work could involve
incorporating the uncertainty into the skeletal point prediction, specifically
while assigning weights to input points.

##  4 Conclusion

In this paper we presented novel technique to leverage recent point-based
learning methods for 3D shape skeletonization. We showed the relative benefits
of utilizing s-reps for supervision and the impact of two new geometric losses
for training more appropriate skeletons. The results on real world data show
promise in terms of both qualitative and quantitative measures and also point
towards future research in terms of adapting to input modality and surface
reconstruction.

##  5 Compliance with Ethical Standards

Heart leaflets from the tricuspid valve were obtained from transthoracic 3DE
images. This modality is part of a standard clinical echo lab protocol for
children with congenital heart disease at Children’s Hospital of Philadelphia.

The hippocampus data was provided by Martin Styner, UNC Neuro Image Analysis
Laboratory, and is publicly available as part of the SlicerSALT and SPHARM-PDM
tool distributions. Original data acquisition was funded by the Stanley
Foundation.

Both studies were performed according to a protocol approved by the
institutional review board at the relevant institutions.

##  6 Acknowledgments

This research was supported by the NIH under project R01EB021391. We thank
Matthew A. Jolley of Children’s Hospital of Philadelphia for lending the
clinical heart valve leaflet data for validation of the algorithms presented
here. Ninad Khargonkar is currently affiliated to the University of Texas at
Dallas but the work contained in this manuscript was performed while he was an
intern at Kitware. The authors have no other financial or non-financial
interests to disclose.

## References

  * [1] Harry Blum,  “A transformation for extracting new descriptions of shape,”  Models for the perception of speech and visual form, pp. 362–380, 1967.
  * [2] Andrea Tagliasacchi, Thomas Delame, Michela Spagnuolo, Nina Amenta, and Alexandru Telea,  “3d skeletons: A state-of-the-art report,”  in Computer Graphics Forum. Wiley Online Library, 2016, vol. 35, pp. 573–597.
  * [3] Yajie Yan, David Letscher, and Tao Ju,  “Voxel cores: Efficient, robust, and provably good approximation of 3d medial axes,”  ACM Transactions on Graphics (TOG), vol. 37, no. 4, pp. 1–13, 2018\.
  * [4] Shihao Wu, Hui Huang, Minglun Gong, Matthias Zwicker, and Daniel Cohen-Or,  “Deep points consolidation,”  ACM Transactions on Graphics (ToG), vol. 34, no. 6, pp. 1–13, 2015\.
  * [5] Hui Huang, Shihao Wu, Daniel Cohen-Or, Minglun Gong, Hao Zhang, Guiqing Li, and Baoquan Chen,  “L1-medial skeleton of point cloud.,”  ACM Trans. Graph., vol. 32, no. 4, pp. 65–1, 2013.
  * [6] Stephen M Pizer, Junpyo Hong, Jared Vicory, Zhiyuan Liu, JS Marron, Hyo-young Choi, James Damon, Sungkyu Jung, Beatriz Paniagua, Jörn Schulz, et al.,  “Object shape representation via skeletal models (s-reps) and statistical analysis,”  in Riemannian Geometric Statistics in Medical Image Analysis, pp. 233–271. Elsevier, 2020.
  * [7] Zhiyuan Liu, Junpyo Hong, Jared Vicory, James N Damon, and Stephen M Pizer,  “Fitting unbranching skeletal structures to objects,”  Medical Image Analysis, vol. 70, pp. 102020, 2021.
  * [8] Wei Shen, Kai Zhao, Yuan Jiang, Yan Wang, Xiang Bai, and Alan Yuille,  “Deepskeleton: Learning multi-task scale-associated deep side outputs for object skeleton extraction in natural images,”  IEEE Transactions on Image Processing, vol. 26, no. 11, pp. 5298–5311, 2017.
  * [9] Oleg Panichev and Alona Voloshyna,  “U-net based convolutional neural network for skeleton extraction,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2019, pp. 0–0.
  * [10] Nam Hoang Nguyen,  “U-net based skeletonization and bag of tricks,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 2105–2109.
  * [11] Cheng Lin, Changjian Li, Yuan Liu, Nenglun Chen, Yi-King Choi, and Wenping Wang,  “Point2skeleton: Learning skeletal representations from point clouds,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 4277–4286.
  * [12] Charles Ruizhongtai Qi, Li Yi, Hao Su, and Leonidas J Guibas,  “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,”  Advances in neural information processing systems, vol. 30, 2017\.
  * [13] Nenglun Chen, Lingjie Liu, Zhiming Cui, Runnan Chen, Duygu Ceylan, Changhe Tu, and Wenping Wang,  “Unsupervised learning of intrinsic structural representation points,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 9121–9130.
  * [14] Jared Vicory, Laura Pascal, Pablo Hernandez, James Fishbaugh, Juan Prieto, Mahmoud Mostapha, Chao Huang, Hina Shah, Junpyo Hong, Zhiyuan Liu, et al.,  “Slicersalt: Shape analysis toolbox,”  in International Workshop on Shape in Medical Imaging. Springer, 2018, pp. 65–72.