Reactive Motion Generation on Learned
Riemannian Manifolds

The International Journal of Robotics
Research
XX(X):1–22©The Author(s) 2021
Reprints and permission:
sagepub.co.uk/journalsPermissions.nav
DOI: 10.1177/ToBeAssigned
www.sagepub.com/

SAGE

Hadi Beik-Mohammadi12 , Søren Hauberg3 , Georgios Arvanitidis3 , Gerhard Neumann2 , andLeonel Rozo1

Abstract
In recent decades, advancements in motion learning have enabled robots to acquire new skills and adapt to unseen
conditions in both structured and unstructured environments. In practice, motion learning methods capture relevant
patterns and adjust them to new conditions such as dynamic obstacle avoidance or variable targets. In this paper, weinvestigate the robot motion learning paradigm from a Riemannian manifold perspective. We argue that Riemannian
manifolds may be learned via human demonstrations in which geodesics are natural motion skills. The geodesics are
generated using a learned Riemannian metric produced by our novel variational autoencoder (VAE), which is intended
to recover full-pose end-effector states and joint space configurations. In addition, we propose a technique for facilitating
on-the-fly end-effector/multiple-limb obstacle avoidance by reshaping the learned manifold using an obstacle-aware
ambient metric. The motion generated using these geodesics may naturally result in multiple-solution tasks that havenot been explicitly demonstrated previously. We extensively tested our approach in task space and joint space scenarios
using a 7-DoF robotic manipulator. We demonstrate that our method is capable of learning and generating motion skills
based on complicated motion patterns demonstrated by a human operator. Additionally, we assess several obstacle
avoidance strategies and generate trajectories in multiple-mode settings.
Keywords
Robot Motion Learning, Manifold Learning, Riemannian Manifolds, Geodesic Motion Skills

1

Introduction

Robot learning has gained interest recently because of
its potential to endow robots with a repertoire of motion
skills to operate autonomously or assist humans in
repetitive, precise, and risk-averse tasks. In this paper,
we tackle the robot motion generation problem from arobot learning perspective, which involves modeling andreproducing motion skills such as picking, placing, etc.
In practice, we may broadly classify motion generation
techniques into two categories: Classic motion planners andmovement primitives. The former generates obstacle-freecontinuous trajectories between start and goal configurations
using search or sampling-based algorithms (Elbanhawi and
Simic 2014; Mohanan and Salgoankar 2018). Movement
primitives, in contrast, are a modular abstraction of robot
movements, where a primitive represents an “atomic action”
or an “elementary movement”, which are often designed
using robot learning techniques (Schaal et al. 2003). A
relevant distinction between the aforementioned approaches
is that motion planners require a precise description of
the robot workspace to generate a continuous path to thegoal, while movement primitives rely on observed motion
patterns to encapsulate the desired trajectory without an
explicit model of the environment. If an unseen obstacle
shows up in the robot workspace, motion planners often
need to replan the entire trajectory from scratch, while most
movement primitives require to be retrained, unless they
are provided with via-point passing features. We aim for
Prepared using sagej.cls [Version: 2017/01/17 v1.20]

a trade-off between these two categories: A learning-based
reactive motion generation that generates robot movements
resembling the observed motion patterns while avoiding
obstacles on the fly.
Human movements often involve following geodesic
paths, which are the most efficient and natural-looking paths.
This is particularly true when completing tasks that require
hand skills, such as reaching for an object. In these cases,
individuals tend to take the path that minimizes the distanceand energy required for the movement. This has been
observed in different studies (Moon et al. 2012; Klein et al.
2022). Therefore, we resort to learning-from-demonstration(LfD) to build a model of robot motions, a technique in which
the skill model is learned by encapsulating motion patterns
from human demonstrations (Osa et al. 2018). The current
LfD approaches do not need an environment model andsome can quickly adapt to settings with dynamic targets (Osa

1 Bosch Center for Artificial Intelligence (BCAI), Renningen, Germany
2 Autonomous Learning Robots Lab, Karlsruhe Institute of Technology(KIT), Karlsruhe, Germany
3 Section for Cognitive Systems, Technical University of Denmark (DTU),
Lyngby, Denmark
Corresponding author:
Hadi Beik-Mohammadi, Bosch Center for Artificial Intelligence (BCAI)
Reinforcement learning and planning (CR/PJ-AI-R31), 70049 Stuttgart,
GERMANY
Email: hadi.beik-mohammadi@de.bosch.com

2

et al. 2018; Ijspeert et al. 2013; Calinon 2016). We canbroadly identify three main groups of LfD approaches thatuse motion primitives as their building blocks. The firstgroup explicitly considers the motion dynamics (Ijspeert
et al. 2013; Manschitz et al. 2018). The second groupbuilds on probabilistic approaches that exploit demonstration
variability for motion adaptation and control (Huang et al.
2019; Calinon 2016; Paraschos et al. 2018). The last groupuses neural networks as their skill model and focuses on
robot motion generalization (Seker et al. 2019; Osa andIkemoto 2019). Despite their significant contributions (seeSection 2), several challenges are still open: Encoding
and reproducing full-pose end-effector movements, skill
adaptation to unseen or dynamic obstacles at both task
and joint spaces, handling multiple-solution (a.k.a. multiplemode) tasks in the aforementioned spaces, and generalization
to unseen situations (e.g. new target poses), among others.
In our previous paper (Beik-Mohammadi et al. 2021),
we provided an LfD approach that addressed several of theforegoing problems by leveraging a Riemannian perspective
for learning robot motions using task space demonstrations.
We employed a Riemannian formulation to represent a
motion skill, in which human demonstrations were assumed
to form a Riemannian manifold (i.e. a smooth surface),
which could be learned in task space R3 × S3 . This approach
showed to have a number of advantages: Strong capabilities
to model the non-linearity of motion data, adaptability tounseen conditions (e.g. obstacles), and generation of noveltask solutions. The core of our previous approach was built
on a variational autoencoder (VAE) that was exploited tocompute a Riemannian metric describing the underlying
structure of these skill manifolds. This metric allowed us
to measure distances on Riemannian manifolds, a property
that is used to compute geodesics (i.e., the shortest pathson Riemannian manifolds), which was then leveraged as our
motion generation mechanism. In addition, we designed anambient space metric to reshape the skill manifold on the
fly in order to perform obstacle avoidance. Furthermore, thelearned skill manifold was able to encode multiple-solution
tasks, which naturally resulted in novel hybrid solutions
based on the synergy of a subset of the demonstrations.
In this paper, we extend our previous work in two
different ways. First, inspired by the need of avoiding
obstacles at any location of the robot body, we propose a
new reactive motion generation method that also leverages
the Riemannian approach proposed in (Beik-Mohammadi
et al. 2021) for joint space skills. To do so, we develop
a new VAE architecture that integrates the robot’s forward
kinematics to access task space information at any point on
the robot body. This new approach makes joint space skill
adaptation possible, allowing the robot to simultaneously
avoid unseen and dynamic obstacles and handle multiplesolution (a.k.a. multiple-mode) tasks in both task and joint
spaces. Second, we describe both methodologies with simple
examples and extensively evaluate them on real tasks using
a 7-DoF robot manipulator. The experiments showcase how
our techniques allow a robot to learn and reconstruct motion
skills of varying complexity. In particular, we show how
our methods enable obstacle avoidance at task and joint
space levels. Moreover, we provide experiments where therobot exploits multiple-solution tasks to effectively generate
Prepared using sagej.cls

The International Journal of Robotics Research XX(X)

Figure 1. From demonstrations we learn a variational
autoencoder that spans a random Riemannian manifold.
Geodesics on this manifold are viewed as motion skills.

hybrid solutions without model retraining. Furthermore, we
experimentally analyze how the latent space dimensionality
impacts the quality of the geodesic motion generator,and how the robot behaves when controlled by geodesics
crossing the manifold boundaries.
In summary, we propose a new Riemannian learning
framework for reactive motion skills that lies on the
middle ground between classical motion planning andmotion primitives. Similar to motion primitives, we exploit
human demonstrations to learn a model of a skill withthe assumption that these demonstrations belong to asmooth surface characterized as a Riemannian manifold.
Additionally, our method, like motion planners, derives
a time-independent reference trajectory from a geodesic
motion generator while incorporating obstacle avoidance
capabilities.
In Section 2, we review relevant work on learningfrom demonstration, variational auto-encoders, Riemannian
manifolds, their connections, and ambient space metrics.
In Section 3, we examine how a skill manifold may be
trained and used to reconstruct robot motion when thedemonstrations are provided in the task and joint spaces.
Section 4 investigates geodesic computation and their ability
to avoid unseen obstacles on the fly. In Section 5, we show
our method capabilities via a series of real-world robot
experiments. Finally, in Section 6, we discuss some of thelimitations of our approach and provide possible solutions to
be explored in future work.

2

Background and Related WorkDue to the strong connections between our work and thelearning from demonstration (LfD) technique, we begin
this section by discussing pertinent work on learning andsynthesizing joint and task space motion skills via LfD.
We look into how these approaches reconstruct full-pose

Beik-Mohammadi et al.

trajectories, provide obstacle avoidance capabilities, andreproduce multiple-solution (multiple-mode) tasks. Second,
we review Riemannian geometry as well as variational
auto-encoders (VAEs), which are essential concepts forour approach. Finally, we explain recent advances in theconnection between VAEs and Riemannian geometry, whichis the backbone of our methods.

2.1

Learning from Demonstration (LfD)

LfD is a robot programming technique that leverages
human demonstrations, recorded via kinesthetic teachingor teleoperation, to learn a model of a task (Ravichandar
et al. 2020). End-effector positions and orientations, joint
configurations, linear or angular velocities, and accelerations
are all examples of data that can be used in LfD. The methods
for dealing with motion dynamics are outside the scopeof this study; Instead, we focus on techniques for learningand synthesizing robot skills built on joint and task space
trajectories.
We identify three key lines of research and their
particular features (e.g. obstacle avoidance) for robot motion
learning. The first group focuses on motion dynamics
learning (Ijspeert et al. 2013; Manschitz et al. 2018),
where demonstrations are considered solutions to specificdynamical systems. These techniques are well-behaved
when confronted with changes in the environment due totheir autonomous systems formulation. The second set of
approaches builds on probabilistic methods that exploit
data variability and model uncertainty (Huang et al. 2019;
Calinon 2016; Paraschos et al. 2018). Their probabilistic
formulation allows robots to execute the skill using alarge diversity of trajectories sampled from the skill model.
The last category includes approaches that use neuralnetworks for increasing generalization in robot motion
learning (Seker et al. 2019; Osa and Ikemoto 2019).
Furthermore, there exist methods that combine dynamicalsystems and neural networks (Bahl et al. 2020; Rana et al.
2020), or dynamical systems and probabilistic models (Ugurand Girgin 2020; Khansari-Zadeh and Billard 2011). Allof the aforementioned methods belong to the categoryof movement primitives (MPs), which is an alternative
approach to classic motion planning (Elbanhawi and Simic
2014) for robot motion generation.
In contrast, we propose a reactive motion generation
technique that lies on a middle ground between movement
primitives and motion planning. Similar to motion primitives, we exploit human demonstrations to learn a modelof a skill with the assumption that these demonstrations
belong to a smooth surface characterized as a Riemannian
manifold. Additionally, our method, like motion planners,
derives a time-independent reference trajectory generated
from geodesic curves, which can be locally deformed to
avoid unseen obstacles. Specifically, our method leverages
a neural network (VAE) to learn a Riemannian metric thatincorporates the network uncertainty. This metric allows us
to generate motions that resemble the demonstrations via
geodesics. The decoded geodesics are then used as reference
trajectories on the robot to reproduce motion trajectories that
resemble the demonstrations.
Complex robot movements may involve sophisticatedfull-pose end-effector trajectories, making it imperative to
Prepared using sagej.cls

3

have a learning framework capable of encoding full-pose
motion patterns. The main challenge is then how to properly
encode and reproduce orientation movements. Despite most
LfD approaches have overlooked this problem (Paraschos
et al. 2018; Seker et al. 2019; Calinon 2016; Huang et al.
2019), recent works have addressed it using probabilistic
models (Rozo and Dave 2021; Zeestraten 2018). In
our previous work (Beik-Mohammadi et al. 2021), we
proposed a VAE architecture capable of encoding full-pose
trajectories, which is here exploited for learning a variety ofreal robotic tasks.
Obstacle avoidance is another feature that a reactive
motion generation should offer. While several approaches
rely solely on obstacles information given before learning (Prescott and Mayhew 1992; Aljalbout et al. 2020), these
are ineffective in dynamic environments. Other techniques
exploit via-points (Paraschos et al. 2018; Seker et al. 2019;
Huang et al. 2019) to tackle this problem in an indirect
fashion, but they do not require retraining the skill model.
A different perspective on the obstacle avoidance problem
is to see obstacles as costs in an optimization frameworkthat seeks to generate optimal and obstacle-free trajectories (Urain et al. 2021; Ratliff et al. 2018).
Meanwhile, our method tackles the obstacle avoidance
problem from a metric-reshaping viewpoint. We design
obstacle-aware ambient space metrics to reshape the learned
Riemannian metric. The combination of these two metrics
yields a new metric that is exploited to generate trajectories
that simultaneously follow the demonstrations while also
avoiding obstacles. Note that the ambient space metricalways represents a notion of distance to the obstacles
in task space. It is worth mentioning that the choice of
demonstration ambient space defines how this metric canbe designed to avoid obstacles. For example, when using
joint space demonstrations, the ambient space metric canincorporate information about the distance from any pointof the robot to the obstacles, resulting in a multiple-limb
obstacle avoidance capability. On the contrary, when using
task space end-effector demonstrations, this metric only
provides obstacle avoidance at the level of the robot endeffector.
Our obstacle avoidance technique employs a techniquesimilar to that of CHOMP (Ratliff et al. 2009), which
involves the use of a simplified geometric representationof the robot and obstacles in order to construct a uniform
grid in the task space. This grid is used to determine thepotential for collision between the trajectory of the robot andany obstacles present. Additionally, our method incorporates
this information into the geodesic energy leading to optimal
paths that adapt accordingly. This approach is conceptuallycomparable to CHOMP, which also utilizes a cost functionthat considers the presence of obstacles.
Human demonstrations can naturally show various
solutions to a single motion skill (Rozo et al. 2011; Seker
et al. 2019), which is typically addressed using hierarchicaltechniques (Konidaris et al. 2012; Ewerton et al. 2015). Our
approach enables the encoding of multiple-solution tasks on
the learned Riemannian manifold, which is then exploited toreplicate the demonstrated skill and generate novel hybrid
solutions based on a synergy of a (sub)set of demonstrations.
These hybrid solutions naturally emerge from our geodesic

4

The International Journal of Robotics Research XX(X)motion generator on the learned Riemannian manifold. Notethat the aforementioned robot motion learning techniques
only provide trajectories that strictly follow the demonstrated
patterns without the ability to generate hybrid solutions.

2.2

Variational auto-encoders (VAEs)

VAEs are generative models (Kingma and Welling 2014)
that learn and reconstruct data by encapsulating their density
into a lower-dimensional latent space Z. Specifically, VAEs
encode the training data density p(x) in an ambient space X
through a low-dimensional latent variable z. For simplicity,
we consider the generative process of a Gaussian VAE
defined as,p(z) = N (z|0, Id ) ,pϕ (x|z) = N

2z|µϕ (z), ID σϕ
(z)



,

z ∈ Z;(1)x ∈ X.

(2)

where µϕ : Z → X and σϕ : Z → RD
+ are deep neuralnetworks with parameters ϕ estimating the mean and thevariance of the posterior distribution pϕ (x|z), and ID and
Id are identity matrices of size D and d, respectively. Since
the exact inference of the generative process is in general
intractable, a variational approximation of the evidence(marginal likelihood) can be used,
LELBO = Eqξ (z|x) [log(pϕ (x|z))]− KL (qξ (z|x)||p(z)) ,(3)qξ (z|x) = N(x|µξ (x), Id σξ2 (x))

whereapproximates theposterior distribution p(z|x) by two deep neural networksµξ (x) : X → Z and σξ (x)) : X → Rd+ . The posterior distribution pξ (z|x) is called inference or encoder distribution,
while the generative distribution pϕ (x|z) is known as thegenerator or decoder. In the next subsection, we use a VAE
to learn a skill-specific Riemannian manifold from humandemonstrations.

2.3

Riemannian Manifolds

In differential geometry, Riemannian manifolds are referred
to as curved d-dimensional continuous and differentiablesurfaces characterized by a Riemannian metric (Lee 2018).
This metric is characterized by a family of smoothly varying
positive-definite inner products acting on the tangent spacesof the manifold, which locally resembles the Euclidean space
Rd . In this paper, we use the mapping function f to represent
a manifold M immersed in the ambient space X defined as,
M = f (Z)with f : Z → X,

(4)where Z and X are open subsets of Euclidean spaces withdim Z < dim X.
An important operation on Riemannian manifolds is thecomputation of the length of a smooth curve c : [0, 1] → Z,
defined as,
Z 1
Lc =∥∂t f (c(t))∥dt.
(5)0This length can be reformulated using the chain rule as,
Z 1qċ(t)T M (c(t))ċ(t)dt,
(6)
Lc =
0

Prepared using sagej.cls

where M and ċt = ∂t ct are the Riemannian metric and curve
derivative, respectively. Note that the Riemannian metriccorresponds to,
M (z) = Jf (z)T Jf (z).

(7)

Here, Jf (z) is the Jacobian of the mapping function f . This
metric can be used to measure local distances in Z. The
shortest path on the manifold, also known as the geodesic,
can be computed given the curve length in Eq. (6). Geodesics
on Riemannian manifolds can be seen as a generalizationof straight lines in Euclidean space. However, geodesics
might not be unique, e.g. great circles on the spheremanifold. Later, we demonstrate that calculating geodesics
on a learned Riemannian manifold can be leveraged torecover demonstrated motion patterns. It should be notedthat geodesics have recently been utilized as solutions oftrajectory optimizers for quadrotor control Scannell et al.
(2021).

2.4

Learning Riemannian Manifolds with VAEs

In this subsection, we examine the link between VAEs andRiemannian geometry. To begin, we first define the VAE
generative process of Eq. (2) as a stochastic function,fϕ (z) = µϕ (z) + diag(ϵ)σϕ (z),ϵ ∼ N(0, ID ).

(8)

where µϕ (z) and σϕ (z) are decoder mean and varianceneural networks, respectively. Also, diag(·) is a diagonal
matrix, and ID is a D × D identity matrix. The above
formulation is referred to as the reparameterization
trick (Kingma and Welling 2014), which can be interpreted
as samples generated out of a random projection of amanifold jointly spanned by µϕ and σϕ , as depicted in Fig. 2.
Riemannian manifolds may arise from mapping functions
between two spaces as in Eq. (4). As a result, Eq. (8)may be seen as a stochastic version of the mapping
function of Eq. (4), which in turn defines a Riemannian
manifold (Hauberg 2019). We can now write the stochastic
form of the Riemannian metric of Eq. (7). To do so, we firstrecast the stochastic function Eq. (8) as follows (Eklund andHauberg 2019), µϕ (z)fϕ (z) = ID , diag(ϵ)= P g(z), (9)σϕ (z)
where P is a random matrix, and g(z) is the concatenationof µϕ (z) and σϕ (z). Therefore, the VAE can be seen as arandom projection of a deterministic manifold spanned by g.
Given that this stochastic mapping function is defined by a
combination of mean µϕ (z) and variance σϕ (z), the metricis likewise based on a mixture of both as follows,M̄ (z) = Jµϕ (z)T Jµϕ (z) + Jσϕ (z)T Jσϕ (z).

(10)

where Jµϕ (z), Jσϕ (z), Jµϕ (z)T , and Jσϕ (z)T are
respectively the Jacobian of µϕ (z) and σϕ (z) and theircorresponding transpose evaluated at z ∈ Z, with Z being
the VAE low-dimensional latent space.
Notably, the decoder variance network σϕ (z) approximates the data uncertainty, which plays a critical role in themetric Eq. (10) by associating low values to regions with aBeik-Mohammadi et al.

5

Figure 2. In a Gaussian VAE,
samples are generated by a random
projection of the manifold jointly
spanned by µ and σ .

high number of data points and vice-versa. Indeed, omitting
this element results in a low-curvature manifold geometry (Hauberg 2019). For example, Shao et al. (2018) suggests a similar technique that learns Riemannian manifolds
using generative models that do not model data uncertainty,
resulting in low-curvature manifolds and often straight lines
as geodesics. As explained earlier, the Riemannian metricis required to compute the geodesics, which conform to thegeometry of the training data (Arvanitidis et al. 2018).
In summary, we exploit the link between VAEs
and Riemannian metrics for robot motion generation.
Specifically, we learn a Riemannian metric that describes
the motion patterns observed during the demonstrations.
These demonstrations may take place in two different
ambient spaces: Task and joint space, which define the VAE
architecture needed to learn the manifold of interest. The
geodesic curves generated on the learned manifold produce
robot movements that mimic the given demonstrations in theambient space.

is then computed as follows,M̄ (z) = Jµϕ (z)T MX (µϕ (z))Jµϕ (z)+ Jσϕ (z)T MX (µϕ (z))Jσϕ (z).

(12)

Given this metric, we can compute geodesics that are
repelled from certain regions of the ambient space X
by increasing the value of the ambient metric MX . We
demonstrate how this metric reshaping method is leveraged
to generate obstacle-free robot motions in Section 5.

3

Riemannian manifold learning

In this section, we describe how learning complex
robot motion skills from demonstrations can be treatedfrom a Riemannian manifold perspective. Unlike previous
works Havoutis and Ramamoorthy (2013); Li et al. (2018),
where skill manifolds are built from locally smooth manifold
learning Dollár et al. (2007), we leverage a Riemannian
formulation. We develop a model that has enough capacity
to learn and synthesize the relevant patterns of a motion
2.5 Ambient space metricwhile being flexible enough to adapt to new conditions(e.g. dynamic obstacles). In this context, related approaches
According to the preceding section, we leverage a VAE tosuch as geometric control methods build on the geometriclearn a Riemannian metric. However, there are situations
properties of a system in the design of control laws. These
where this metric needs to be changed. For example, if our
methods are particularly useful in situations where themetric encapsulates the main patterns of the robot motion
dynamics of the system are complex and difficult to model
demonstrations, we may be interested in endowing the robot
accurately (Bullo and Lewis 2004).
with obstacle-avoidance capabilities. This means that we
On a related note, our approach shares some conceptual
now require to reshape the manifold to take these obstacles
aspects
with optimal control theory. Broadly speaking,into account. To do so, we need to reshape the previouslyoptimal
control aims at finding optimal control inputs
learned metric so that the new geodesics lead to obstacle-freethatminimizea specific cost functional (Kirk 1970). The
robot motions.
connection between inverse optimal control and Riemannian
A naı̈ve approach would entail retraining the VAE model
manifold learning is that the learned Riemannian metricwith new data, which is time-consuming and data-intensive,can be understood as the cost function that the humanand can only be executed offline. We propose to reshapedemonstrator is optimizing. More specifically, the optimal
the learned metric by considering problem-specific ambient
expert demonstrations of a motion are used to learn ametrics. Note that although the definition of curve length
Riemannian metric that defines the optimal or shortest
relies on the Euclidean metric of X, this is not a strict
path (i.e., geodesics) that minimizes the functional inrequirement. Indeed, Arvanitidis et al. (2021) argued thatEq. (6). In comparison, our method differs conceptually
there is value in giving the ambient space a manually-definedfrom the aforementioned control techniques, as it leverages
Riemannian metric and including that into the definition
geodesic on a learned Riemannian manifold for robot motionof curve length. The resulting metric can then be used togeneration. While all approaches aim to achieve precise andcompute the curve length as,
efficient control or planning for mechanical systems, they doso through different perspectives. We next describe how we
Z 1qTTċ(t) Jfϕ (c(t)) MX (fϕ (c(t)))Jfϕ (c(t))ċ(t)dt, use VAEs to access a low-dimensional learned manifold of
Lc =the demonstrations to learn an ambient space Riemannian0(11) metric.
As mentioned previously, this metric is exploited towhere MX is the ambient metric, which can now vary reconstruct robot motions in both task space R3 × S3 andsmoothly across X. The reshaped Riemannian metric in Z joint space Rη . In addition, we discuss the design of the
Prepared using sagej.cls

6

The International Journal of Robotics Research XX(X)corresponding VAE for each ambient space as well as theformulation of the corresponding Riemannian metric. In
Section 4, we explain how we exploit these learned metrics
to generate robot motion trajectories using geodesics.

3.1

Task space R3 × S3

To begin, we focus on learning motion skills characterized
by full-pose end-effector trajectories, where each poseis represented in R3 × S3 . Before exploiting the VAE
to compute the Riemannian metric, we must ensure its
capability to properly learn and reconstruct full-pose endeffector states, i.e. position x ∈ R3 and orientation q ∈ S3 ,while accounting for specific properties of the data, suchas quaternions antipodality. To do so, we propose a VAE
architecture that models the joint density of the robot endeffector state. Our model retains the usual Gaussian priorp(z) = N(z|0, Id ), but modifies the generative distribution
pϕ,ψ (x, q|z). Specifically, we assume that position andorientation are conditionally independent,
pϕ,ψ (x, q|z) = pϕ (x|z)pψ (q|z),(13)where the latent variable z captures the correlation between
position and quaternion data. Next, we describe how eachconditional distribution is parameterized and learned.
3.1.1

where CD is the normalization constantD

CD (κ) =Position encoding in R3 : To model the conditional

distribution of end-effector positions x, we opt for simplicity
and choose this to be Gaussian,2
pϕ (x|z) = N(x|µϕ (z), I3 σϕ
(z)),

(14)

where µϕ and σϕ are neural networks parametrized by ϕ.
This is a somewhat simplistic model as the Gaussian model
will assign probability mass outside the workspace of therobot, i.e. physically infeasible robot states. With simplicity
in mind, we disregard this concern as we only optimize theassociated likelihood, and therefore never sample infeasible
states. However, if we were to use the VAE as a fully
generative model, this likelihood should be replaced by amore elaborated model that accounts for physical properties.
3.1.2

antipodal identification that q and −q correspond to thesame orientation. Formally, a unit quaternion q lying on
the surface of a 3-sphere S3 can be represented using a
4-dimensional unit vector q = [qw , qx , qy , qz ], where thescalar qw and vector (qx , qy , qz ) represent the real andimaginary parts of the quaternion, respectively. To copewith antipodality, one could model q as a point in a
projective space, but for simplicity we let q live on theunit sphere S3 . We then choose a generative distributionpψ (q|z) such that pψ (q|z) = pψ (−q|z). In other words,
the quaternions q and −q are considered to be antipodal: theylie on diametrically opposite points on the 3-sphere while
representing the same orientation.
To formulate a suitable distribution pψ (q|z) over S3 , weleverage the von Mises-Fischer (vMF) distribution, which ismerely an isotropic Gaussian constrained to lie on the unit
sphere (Sra 2018). This distribution is described by a mean
direction µ with ∥µ∥ = 1, and a concentration parameterκ ≥ 0. The vMF density function is defined as,
vMF(q|µ, κ) = CD (κ) exp κµT q ,∥µ∥= 1, (15)

Quaternion encoding in S3 : On a robot motion

trajectory, each position is paired with an orientation, andtogether they define the full pose of the end-effector.
There are several representations for the end-effector
orientation, for example, Euler angles, rotation matrices,and unit quaternions. Euler angles and rotation matrices arewidely used for their simplicity and intuitiveness, however,the former suffer from gimbal lock (Hemingway andO’Reilly 2018) which makes them an inadequate orientation
parametrization, and the latter are a redundant representation
requiring a high number of parameters.
Unit quaternions, on the other hand, are a convenient
way to represent orientations since they are compact,
not redundant, and prevent gimbal locks. Also, they
provide strong stability guarantees in closed-loop orientation
control (Campa and Camarillo 2008), and they havebeen recently exploited in complicated robotic tasks
learning (Rozo et al. 2020), and for data-efficient robot
control tuning (Jaquier et al. 2020) using Riemannianmanifold formulations. We choose to represent orientations
q as a unit quaternion, such that q ∈ S3 with the additional
Prepared using sagej.cls

κ 2 −1D(2π) 2 I D −1 (κ)

,

(16)

2

with I D −1 (κ) being the modified Bessel function of the first2
kind. Like the Gaussian, from which the distribution wasconstructed, the von Mises-Fischer distribution is unimodal.
To build a distribution that is antipodal symmetric, i.e.
pψ (q|z) = pψ (−q|z), we define a mixture of antipodalvMF distributions (Hauberg et al. 2016),
1
vMF(q|µψ (z), κψ (z))
2
1
+ vMF(q|−µψ (z), κψ (z)),
2pψ (q|z) =(17)

where µ and κ are parametrized as neural networks.
This mixture model is conceptually similar to a Bingham
distribution (Sra 2018), but is easier to implement
numerically.
Variational inference: To train the VAE, we maximize an adapted evidence lower bound (ELBO) Eq. (3),defined as

3.1.3

LELBO = β1 Lx + β2 Lq − KL (qξ (z|x)||p(z)) ,(18)

Lx = Eqξ (z|x) [log pϕ (x|z)] ,(19)

Lq = Eqξ (z|x) [log pψ (q|z)] ,(20)

where x ∈ R3 and q ∈ S3 represent the position andorientation of the end-effector, respectively. The scaling
factors β1 > 0 and β2 > 0 balance the log-likelihood of
position and orientation components. This approach isinspired by the β-VAE method (Higgins et al. 2017),
which, despite not yielding a valid lower bound on the
ELBO objective function, has been shown to produce highquality results comparable to methods that provide a valid
lower bound, thus it is commonly used in practice. Due tothe quaternions antipodality, raw demonstration data may
contain negative or positive values for the same orientation.

Beik-Mohammadi et al.

7

So, we avoid any pre-processing step of the data byconsidering two vMF distributions that encode the sameorientation at both sides of the hypersphere. Practically,
we double the training data, by including qn and −qn for
all observations qn . Note that as the Riemannian manifold
is learned using task space data, the model is kinematics
agnostic, which means that the generated motion may be
used across different robots as long as the trajectory isreachable.
Induced Riemannian metric: Our generative process is parametrized by a set of neural networks. Specifically,µϕ and σϕ are position mean and variance neural networks
parameterized by ϕ, while µψ and κψ are neural networks
parameterized by ψ that represent the mean and concentration of the quaternion distribution. Following Section 2.4, theJacobians of these functions govern the induced Riemannian
metric as,

3.1.4

M (z) = Mµx (z) + Mσx (z) + Mµq (z) + Mκq (z), (21)demonstrations. Similarly, a new VAE architecture can be
designed to compute a Riemannian metric from joint spacedemonstrations. We use this metric to compute geodesics andgenerate robot movements that resemble the demonstrations
in joint space. This joint space approach also allows us
to endow the robot with whole-body obstacle avoidance
capabilities. By using ambient metrics, we can again
reshape the learned metric to make the robot move awayfrom obstacles in an online fashion. The ambient metrics
exclusively use task space information of the obstacles andthe robot body, in contrast to classical motion planningthat often works in the configuration space. Note that thedata manifold learned using joint space demonstrations iskinematics-dependent, meaning that the generated motion
cannot be directly transferred to other robots with different
kinematics.
Variational inference: To train the joint space VAE,
we maximize a modified version of the evidence lower bound(ELBO) Eq. (3), defined as
3.2.1

LELBO = Lθ − KL (qξ (zi |θi )||p(z)) ,with
Mµx (z) = Jµϕ (z)T Jµϕ (z),(22)

Mσx (z) = Jσϕ (z)T Jσϕ (z),(23)

Mµq (z) = Jµψ (z)T Jµψ (z),Mκq (z) = Jκψ (z)T Jκψ (z),

(24)
(25)

where Jµϕ , Jσϕ , Jµψ , Jκψ are the Jacobians of the functions
representing the position mean and variance, as well as thequaternion mean and concentration, respectively.
In practice, we want this Riemannian metric M (z) totake large values in regions with little or no data, so thatgeodesics avoid passing through them. We achieve this byusing radial basis function (RBF) networks as our variance
representation, whose kernels reliably extrapolate over thewhole space (Arvanitidis et al. 2018). Since one of themain differences between Gaussian and von Mises-Fischer
distributions is the representation of the data dispersion,
the RBF network should consider a reciprocal behavior
when estimating variance for positions. In summary, the data
uncertainty is encoded by the RBF networks representing−1σϕ(z) and κψ (z), which affect the Riemannian metricthrough their corresponding Jacobians as in Eq. (21).

3.2

Joint space Rη

The joint space Rη , also known as configuration space, isanother space to represent robot motion trajectories* . In thisspace, each trajectory point is represented as the vector θ =[θ1 , θ2 , . . . , θη ] ∈ Rη , where η is the number of degrees
of freedom of the robot. Learning motion skills in thisspace is known to be challenging as it is less intuitive toprovide joint-level demonstrations. However, being able tolearn and generate joint space movements is relevant as sometasks may demand specific robot postures. Moreover, joint
space skills can be extended to provide whole-body obstacle
avoidance. In this context, we formulate a Riemannian robot
motion learning approach to generate collision-free joint
space movements.
Previously, we computed a Riemannian metric in thelatent space using the VAE decoder trained on task space
Prepared using sagej.cls

(26)

with,
Lθ = Eqξ (zi |θi ) [pX (fFK (θ)|zi )] ,= Eqξ (zi |θi ) [log(pΘ (θ|zi )) − log(V)] ,where pΘ (θ|zi ) and pX (fFK (θ)|zi ) are the estimatedconditional densities in the joint space Θ and task space X,
respectively. Also, V is the volume measure defined asq(27)
V = (det(JfTFK (µα (zi ))JfFK (µα (zi ))),
where JfFK is the Jacobian of the forward kinematics
fFK given the joint configuration estimated by the
VAE decoder µα . Furthermore, the generative distributionpΘ (θ|zi ) = N(µα (zi ), Iη σα (zi )2 ) is parameterized by the
VAE decoder mean µα (zi ) and variance σα (zi ) networks.
Note that the new ELBO formulation in Eq. (26) leverages
the change of variable theorem (Deisenroth et al. 2020) totransform probability densities from joint to task space. As aresult, the VAE is still trained using task space information,
while the given demonstration trajectories are defined in joint
space. This is motivated by the fact that most robot skills may
still depend on task space variables (e.g. the manipulated
objects pose), despite the same skill is also required to imitate
particular robot postures.
As we are interested in whole-body obstacle avoidance,
we can leverage the forward kinematics model to access the
Cartesian position of different points on the robot (e.g., joint
locations). Therefore, we use a set of M forward kinematic
1:M
1:nfunctions fFK(µ1:η
α ), where µα is n = 1, . . . , η elementsof the joint-values vector µα (z), and M is the number ofconsidered points on the robot. Note that for certain points
on the robot structure, the forward kinematics only needs asubset of the joint values. For simplicity, we consider M tobe equal to the number of robot joints plus the end-effector∗ We did not model the robot joint space as a high-dimensional torus for

simplicity. However, we showed that our approach can easily encode dataon Riemannian manifolds as in the task space case.

8

The International Journal of Robotics Research XX(X)(i.e. M = η + 1). Then, the full forward kinematic function
fFK is defined as,

Riemannian-metric-based Geodesic
Euclidean-metric-based Geodesic
Encoded Demonstrations
12

 1

M
ee
1:η T
(µ1α ), . . . , fFK(µ1:η
fFK (µα ) = fFK,α ), fFK (µα )T= p1 , . . . , pM , pee , qee ,10where given the joint value vector µ1:n
α as input, all thefunctions compute the corresponding position pm of the meeth point on the robot, except the last function fFK
which also
provides both the position pee and the orientation qee of theend-effector.
Furthermore, the volume measure V in Eq. (27) uses theJacobian of the full forward kinematics function, which isdefined as,

8

6

4

2


T
JfFK (µα ) = Jp1 , . . . , JpM , Jpee , Jqee ,where Jpi and Jqi are the linear and angular components ofthe corresponding Jacobians.

22

Induced Riemannian metric: With the new integrated forward kinematic layer, we can calculate a pullback
metric that directly uses task space information. This adds anadditional step in the formulation of the Riemannian metric,
which now requires the Jacobian of the forward kinematics
JfFK as well as the Jacobians of the VAE decoder Jµα
and Jσα , computed from the mean and variance decoder
networks. Using these two Jacobians the metric can be
defined as,
3.2.2

M θ (z) = Mµθα (z) + Mσθα (z)TMµθα (z) = (JfFK (µα (z))Jµα (z)) (JfFK (µα (z))Jµα (z)) ,T

Mσθα (z) = (JfFK (µα (z))Jσα (z)) (JfFK (µα (z))Jσα (z)) .
Similarly to our Riemannian metric in task space, this new
metric M θ (z) takes large values in regions with little
or no data so that geodesics avoid passing through them.
Therefore, geodesic curves generated via Eq. (28) allow us
to reproduce joint space robot skills.

Geodesic motion skills

As explained previously, geodesics follow the trend of thedata, and they are here exploited to reconstruct motion
skills that resemble human demonstrations. In this section,
we describe geodesic computation in both settings, namely,
where the VAE is trained on task space or joint space
trajectories. Moreover, we explain how new geodesic paths,
that avoid obstacles on the fly, can be obtained by metricreshaping. In particular, we exploit ambient space metrics
defined as a function of the obstacle’s configuration to locally
deform the original learned Riemannian metric. kinematic
Last but not least, our approach can encode multiplesolution skills, from which new hybrid trajectories (notpreviously shown to the robot) can be synthesized. To
elaborate, a multiple-solution task requires the robot to
combine multiple demonstrated solutions, i.e. trajectories,
in order to complete it. By being able to combine multiple
solutions to generate a mixed novel trajectory, the robot
Prepared using sagej.cls

18

16

(28)

with,

4

20

Figure 3. Top: The variance measure, bottom: The
magnification factor of the Riemannian manifold learned fromtrajectories based on J and C English alphabet charactersdefined in R2 × S2 . The semi-transparent white points depict
the encoded training set, and the yellow curve depicts thegeodesic in the latent space. The resulting manifold iscomposed of two similar clusters due to the antipodal encodingof the quaternions, where each cluster represents one side of
the hyper-sphere. The yellow and red curves show thegeodesics computed based on Riemannian and Euclidean
metrics, respectively.

can find a novel solution to a task despite that not a single
demonstration provides a complete solution on its own.
This can be particularly useful in complex or dynamic
environments where the robot needs to be able to adapt andrespond to changing conditions. We elaborate on each ofthese features in the sequel.

4.1

Generating motion:

Robot motion generation techniques that leverage demonstrations aim to replicate the demonstrated movement patterns with a high degree of accuracy, reliability, and efficiency. This is achieved by exploiting the demonstrations
as a source of information about the range of (manifold of)possible movements that the robot should perform. The robot
uses this information to generate similar movements on its
own.
One way to represent the observed movement patterns
mathematically is through the use of Riemannian manifolds.

Beik-Mohammadi et al.

Riemannian manifolds allow us to describe data sub-spaces(defined by the demonstrations) within a larger spaceand provide analytical tools for calculating distances andangles within these sub-spaces. In the context of robot
motion generation, a Riemannian manifold can be used tocharacterize the region of the state space covered by thedemonstrations.
In this paper, the Riemannian metrics that describe thestructure of these manifolds are highly influenced by the data
uncertainty and they are used to find the path between two
points on the manifold, the so-called geodesic. In the context
of Riemannian manifolds, geodesics are the curves thatminimize the distance between two points on the manifold.
It is important to keep in mind that the distance measurement
is done on the manifold, using the distance metric of thephysical space where the robot moves. This means thatgeodesics on a Riemannian manifold represent the mostefficient paths through the space represented by the manifold.
Given that Riemannian manifolds represent the sub-spacewithin which the demonstrations reside, it follows that thegeodesics are equivalent to the optimal motion trajectories.
Specifically, the Riemannian metrics Eq. (10) and Eq. (28)tell us that geodesics are penalized for crossing throughregions where the VAE predictive uncertainty grows. Thisimplies that if a set of demonstrations follows a circularmotion pattern, geodesics starting from arbitrary points on
the learned manifold will also generate a circular motion (seeFig. 1). This behavior is due to the way that the metric M isdefined: Our Riemannian metric M is characterized by low
values where data uncertainty is low (and vice-versa). Since
the geodesics minimize the energy of the curve between two
points on M, which is a function of M , they tend to stay on
the learned manifold and avoid outside regions.
This property follows the common characteristics ofmotion generation techniques to reproduce motion learned
from demonstrations and, makes us suggest that geodesics
form a natural motion generation mechanism. Note thatwhen using a Euclidean metric (i.e., an identity matrix),
geodesics correspond to straight lines. Such geodesics
certainly neglect the data manifold geometry.
Noted that geodesics do not typically follow a closedform equation on these learned manifolds, and numericalapproximations are required. This can be done by direct
minimization of curve length (Shao et al. 2018; Kalatzis
et al. 2020), A∗ search (Chen et al. 2019), integration of
the associated ODE (Arvanitidis et al. 2019), or various
heuristics (Chen et al. 2018). In this paper, we compute
geodesics on M by approximating them by cubic splinesc ≈ ωλ (zc ), with zc = {zc0 , . . . , zcK }, where zck ∈ Z isa vector defining a control point of the spline over the latent
space Z. Given K control points, K − 1 cubic polynomials
ωλi with coefficients λi,0 , λi,1 , λi,2 , λi,3 have to be estimatedto minimize its Riemannian length,
Z 1p⟨ω̇λ (zc ), M (ωλ (zc ))ω̇λ (zc )⟩dt. (29)Lωλ (zc ) =0

The resulting geodesic c computed in Z is used to generate
the robot motion by decoding it through the VAE networksµψ and µψ or µα depending on the ambient space setting.
The obtained trajectory is then executed on the robot arm
to reproduce the required skill. In the task space setting,
Prepared using sagej.cls

9

Figure 4. As an illustration, we consider synthetic data thatbelong to R2 × S2 . The left panel depicts the J-shaped position
data in R2 and the right panel shows the C-shaped orientation
data on S2 . The yellow and red curves show the decoded
geodesics depicted in Fig. 3, computed according to the
Riemannian and Euclidean metrics, respectively.

the decoded geodesics can be deployed on the robot using a
Cartesian impedance controller or inverse kinematics. In the
joint space setting, the decoded geodesics can be employed
directly on the robot as a joint trajectory reference to be
tracked by joint position or impedance controllers.

4.2

Geodesics in task space R3 × S3 :In this section, we investigate the geodesic motion generationin task space. To illustrate the motion generation mechanism,
we consider a simple experiment where the demonstration
data at each time point is confined to R2 × S2 , i.e. only
two-dimensional positions and orientations are considered.
We artificially create position data that follows a J shapeand orientation data that follows a C shape projected on thesphere (see Fig. 4). We fit our VAE model to this dataset,and visualize the corresponding latent space in Fig. 3, where
the top panel shows the latent mean embeddings of the
training data with a background color corresponding to thepredictive uncertainty. We see low uncertainty near the data,and high otherwise. The bottom panel of Fig. 3 shows thesame embedding
but with a background color proportional√
to log det M . This quantity, known as the magnification
factor (Bishop et al. 1997), generally takes large values inregions where distances are large, implying that geodesics
will try to avoid such regions. In Fig. 3, we notice thatthe magnification factor is generally low, except on the‘boundary’ of the data manifold, i.e. in regions where thepredictive variance grows. Consequently, we observe thatRiemannian geodesics (yellow curves in the figure) staywithin the ‘boundary’ and hence resemble the training data
patterns. In contrast, Euclidean geodesics (red curves in the
figure) fail to stay in the data manifold. Our proposal is to use
these geodesics on the learned manifolds as our robot motion
generation mechanism.
Note that both panels in Fig. 3 depict two distinct
horseshoe-like clusters, which is a result of the antipodalityof the data in S2 . More precisely, the bimodal distributionof the antipodal quaternion data is encapsulated by these twoclusters in the latent space. In practical settings, as long as thegeodesic curve does not cross across clusters (both start and10

The International Journal of Robotics Research XX(X)

goal points belong to the same cluster), the quaternion signis unchanged. We experimentally examine this in Section 5.

4.3

Geodesics in joint space Rη :In this section, we investigate the geodesic computation
for joint space movements. We use a toy example wherea 2-DOF robot arm follows an S-shaped trajectory in task
space using two different joint configurations (i.e., twodifferent inverse-kinematics solutions), as shown in Fig. 5–
left. We observe two sets of demonstrations in joint spacethat reproduce the same end-effector movements whenapplied to the robot. We can also see in the middle andright panels of Fig. 5 the geodesics computed using the
Riemannian and Euclidean metrics, depicted as blue and red
curves, respectively. The background of the middle panel
illustrates the predictive uncertainty over the latent space
Z, where we again see low uncertainty near the data, andhigh otherwise. For completeness, the background in theright panel illustrates the magnification factor. The latentmean embedding of the training data is depicted as semitransparent white points. Similar to the previous section, thegeodesics generated using the Riemannian metric stay withinthe ‘boundary’ near the training data.
Furthermore, it is easy to note that the learned manifold
comprises two clusters, but unlike the previous task spaceexample, these clusters arise from the two different joint
space solutions provided in the training data. This indicatesthat the clusters in the learned manifold encapsulate theprovided solutions in the demonstrations. When the numberof clusters grows, the geodesic has a higher chance to travel
among them to find a path with minimal energy as the highenergy regions may become narrow. However, unnecessary
frequent switching among these clusters may often lead tojerky geodesics, therefore negatively impacting the geodesics
quality, particularly in robots with a high degree of freedom(e.g. DOF ≥ 7). Later in Section 5, we experimentally showthat this issue can be alleviated by increasing the latent spacedimensionality.

4.4

Obstacle avoidance using ambient space
metrics

Often human demonstrations do not include any notion of
obstacles in the environment. Therefore, obstacle avoidanceis usually treated as a separate problem when generating
robot motions in unstructured environments. A possible
solution to integrate both problems is to provide obstacle–aware demonstrations, where the robot is explicitly taught
how to avoid known obstacles. The main drawback is thatthe robot is still unable to avoid unseen obstacles on the fly.
Our Riemannian approach provides a natural solution tothis problem. The learned metrics in latent space Eq. (10)and Eq. (28) measure the length of a geodesic curve under
the Euclidean space of the ambient space X. We can easily
modify this to account for unseen and dynamic obstacles.
Intuitively, we can increase the length of curves that intersect
obstacles, such that geodesics are less likely to go near them.
Next, we explain how obstacle avoidance can be achieved forboth task space and joint space settings.
4.4.1

Obstacle avoidance in task space R3 × S3 : Here

we explain how we can reshape the learned metric to avoid
Prepared using sagej.cls

obstacles in the task space setting, where only the robot endeffector is considered. Formally, we propose to define theambient metric of the end-effector position to be−∥x − o∥2x
I3 , x ∈ R 3 ,MX (x) = 1 + ζ exp
2r2(30)
where ζ > 0 scales the cost, o ∈ R3 and r > 0 represent
the position and radius of the obstacle, respectively. For
the orientation component, we assume a flat ambient metric
MXq (x) = I4 . Under this new ambient metric, geodesics will
generally avoid the obstacle, though we emphasize this is
only a soft constraint. This approach is similar in spirit to
CHOMP (Ratliff et al. 2009) except our formulation works
along a low-dimensional learned manifold, whose solution isthen decoded to the task space of the robot.
Under this ambient metric, the associated (reshaped)
Riemannian metric of the latent space Z becomes,
M (z) = Mµx (z) + Mσx (z) + Mµq (z) + Mκq (z), (31)with Mµx (z) = Jµϕ (z)T MXx (µϕ (z))Jµϕ (z),Mσx (z) = Jσϕ (z)T MXx (µϕ (z))Jσϕ (z),Mµq (z) = Jµψ (z)T MXq (µψ (z))Jµψ (z),Mκq (z) = Jκψ (z)T MXq (µψ (z))Jκψ (z),where MXx and MXq represent the position and orientation
components of the obstacle-avoidance metric MX , respectively. We emphasize that as the object changes position,
the VAE does not need to be re-trained as the change isonly in the ambient metric. As stated before, obstacles can
be avoided only by the end-effector under this task spacesetting.
In the multiple-limb obstacle avoidance setting, rather
than just using the end effector, the entire body of therobot is taken into account. This technique takes advantageof redundant solutions in joint space to find a robot
configuration that avoids collisions. This is useful insituations where the end effector may be safe from an
obstacle, but other parts of the robot’s body (such as its links)
are in danger of colliding with the obstacle. By providing
multiple solutions in joint space during the demonstration
phase, the robot can choose a configuration that avoids the
obstacle and continue with its task. In cases where thedemonstrations do not provide joint space solutions, a varietyof obstacle avoidance techniques can be used to move therobot away from obstacles.
There are a variety of methods to incorporate obstacle
avoidance into robotic systems. One commonly employedmethod is to utilize an off-the-shelf obstacle avoidancetechnique as post-processing of the generated trajectory.
This approach involves first generating a trajectory, andthen applying an obstacle avoidance technique to prevent
collisions with obstacles present in the environment. Thismethod can be relatively straightforward to implement andmay effectively avoid obstacles. However, it can also lead todeviation from the intended motion and may not accuratelyreflect the demonstrations. It is worth noting that obstacle
avoidance techniques that operate in joint space are often

Beik-Mohammadi et al.

11

Figure 5. Illustration of joint space motion generation via geodesics. Left: S-shaped joint space demonstrations. Middle: theresulting variance measure. Right: The magnification factor of the learned Riemannian manifold. The semi-transparent white points
depict the encoded training set and the blue curve represents the geodesic in the latent space. The resulting manifold is composedof two similar clusters due to the two different inverse-kinematics solutions for the task. The blue and red curves show thegeodesics computed based on Riemannian and Euclidean metrics, respectively.

The technique that will be outlined in the following
section offers a capability for efficient and effective obstacle
avoidance for multiple limbs, without reliance on additional
obstacle avoidance techniques.

Encoder

Decoder
𝜇𝜙

𝓍

𝜇𝜉

𝑞

𝜎𝜉

𝓍

Latent
Space

𝜎𝜙

𝑧

Joint
space

𝜇𝜓

𝑞
𝜅𝜓

Variational Autoencoder
Encoder

Decoder

Joint
space

ForwardKinematic(Untrainable)

Task
space

Obstacle avoidance in joint space Rη : Avoiding

obstacles at the robot link level while performing motion
skills requires considering the whole robot’s kinematic
structure. Classical motion planning methods model thegeometry of the obstacles into the configuration spaceand later compute an obstacle-free path via sampling
methods (Elbanhawi and Simic 2014). In contrast, we takeadvantage of the forward kinematics layer (see Fig. 6bottom), which provides us with task space poses of any
point on the robot body, to compute an obstacle-avoidance
ambient metric. Similar to the task space formulation
presented previously, this ambient metric is then exploited toreshape the learned metric and generate modified geodesic
curves that produce collision-free robot movements.
Specifically, we need to define a collection of points on
the robot body p1 , . . . , pM with pm ∈ R3 . These points arethen used to compute the ambient space metric for obstacleavoidance purposes. Therefore, a larger collection of points
provides a more robust obstacle-avoidance performance at
the cost of higher computational complexity. Given the setof points of interest, we compute an associated ambient
metric following Eq. (30) with x = pm . Similar to thetask space setting, since the orientation of obstacles is notconsidered, the corresponding ambient space metric is an
identity matrix. Finally, we form the whole ambient
metric
as MX = blockdiag( MXp1 , MXp2 , · · · , MXpM ), which isthen used to reshape the learned metric of Eq. (28) as,
M (z) = Mµθ (z) + Mσθ (z),Prepared using sagej.cls

(32)

𝜇𝜉

𝓍

𝜇𝛼

Latent
Space

𝜃

𝜃

𝑧

𝑞
𝜎𝛼

𝜎𝜉
Pullback
Metric

4.4.2

Task
space

Variational Autoencoder

Task
space

Pullback
Metriccomputationally intensive. Furthermore, this method cannot
be applied in the latent space, as the obstacle avoidancetechnique must have knowledge of geometry and be able togenerate trajectories on the manifold.

Figure 6. Top: The VAE architecture under the task spacesetting. The blue, orange, green, and gray blocks correspond toambient spaces, functions with trainable parameters, latentvariables, and functions with fixed parameters. The arrows
indicate the direction in which the data flows during the query.
Bottom: The architecture of the VAE under the joint spacesetting.

with,TMµθ (z) = (JfFK (µα )Jµα (z)) MX (JfFK (µα )Jµα (z)) ,TMσθ (z) = (JfFK (µα )Jσα (z)) MX (JfFK (µα )Jσα (z)) .
Generating geodesics on discrete manifolds: To
generate robot motion, our approach requires computing
geodesics, which can be done in several ways (Peyré et al.
2010). While this can be computed by solving an ordinary
differential equation, this approach may not be suitable forreal-time robotic applications due to its computational cost,
as highlighted in previous studies (Arvanitidis et al. 2019).
A commonly used alternative is to employ gradient descent
to minimize the curve length (Noakes and Zhang 2022).
While this approach is straightforward and effective in some
4.4.3

12

The International Journal of Robotics Research XX(X)cases, it may not be optimal when operating in real-time.
In our experiments, we found that the gradient descentapproach often gets trapped in local minima, depending on a
non-trivial initialization. As a consequence, multiple restarts
may be required, which increases computational cost. We,
therefore, consider an alternative discrete approach.
As we work with low-dimensional latent spaces, wehere propose to simply discretizing them on a regular gridand use a graph-based algorithm for computing shortest
paths. If this grid is sufficiently dense, this approach isglobally optimal and computationally inexpensive. However,this approach is only feasible in low-dimensional latent
spaces due to the curse of dimensionality. In particular,
the memory requirements of the graph grow exponentiallywith dimension, and the short path search grows with thegraph size. Our experiments are, however, concerned with
low-dimensional latent spaces, where these limitations are
negligible.
In detail, we create a uniform grid over the latent space Z,
and assign a weight to each edge in the graph corresponding
to the Riemannian distance between neighboring nodes(see Fig. 8). Geodesics are then found using Dijkstra’salgorithm (Cormen et al. 2009). This algorithm selects a setof graph nodes,
Gc = {g0 , g1 , . . . , gS−1 , gS } ,g s ∈ RD ,where g0 and gS represent the start and the target of thegeodesic in the graph, respectively. To select these points,
the shortest path on the graph is calculated by minimizing
the accumulated weight (cost) of each edge connecting two
nodes, computed as in Eq. (6). To ensure a smooth trajectory,
we fit a cubic spline ωλ to the resulting set Gc by minimizing
the mean square error. The spline computed in Z is finallyused to generate the robot moves through the mean decoder
networks: µθ and µψ or µα . The resulting trajectory can be
executed on the robot arm to reproduce the required skill.
In order to verify the validity of the discrete geodesic
approximation (graph-based) in comparison to the continuous approximation (gradient-based), a toy example scenario
was devised for evaluating their respective accuracies. Threedifferent experimental setups were compared, one using
gradient descent and two using a graph-based approach. In
the graph-based setting, the manifold was discretized, andthe Dijkstra algorithm was used to find the shortest path. The
discrete setups were also evaluated using two different resolutions to assess the effect of the grid density on the accuracy
and efficiency. We employed the same pre-trained modelfrom the toy example presented in Section 4.2. Figure 7
displays the resulting geodesics. As observed, all methods
were successful in approximating a geodesic that follows thetrend of the data in the latent space. However, the geodesic
generated using the 50 × 50 grid, depicted in yellow, exhibits
an unnecessary deviation from the data, likely due to therelatively low resolution of the grid. However, the geodesics
generated using the 100 × 100 grid and gradient descentmethod exhibit a comparable level of accuracy.
Moreover, one issue with the graph-based approach is thatdynamic obstacles imply that geodesic distances between
nodes may change dynamically. To avoid recomputing alledge weights every time an obstacle moves we do as follows.
Since the learned manifold does not change, we can keep a
Prepared using sagej.cls

22

20

18

16

50 x 50
100 x 100
Gradient Decent
Encoded training set

Figure 7. The geodesics generated under continuous anddiscrete setups. The results show that the geodesic computed
using the 50 × 50 grid exhibits an unnecessary deviation fromthe data, likely due to the relatively low resolution of the grid. In
contrast, the geodesic generated using the 100 × 100 grid has
a level of accuracy that is comparable to the (continuous)
gradient descent method.

decoded version of the latent graph in memory (Fig. 8). Thisway we avoid querying the decoders at run-time. We can thenfind points on the decoded graph that are near obstacles andrescale their weights to take the ambient metric into account.
Once the obstacle moves, we can reset the weights of points
that are now sufficiently far from the obstacle.

5

Experiments

In this section, we showcase the capabilities of our method
in different experiments where human demonstrations areprovided in both task and joint spaces. We provide afull description of the experimental setup, which includes
the design of the VAE networks, Riemannian manifold
learning, real-time geodesic computation, and calculationof the ambient metric to avoid obstacles. We evaluate theperformance of our approach using real-robot experiments,
namely, Reach-to-grasp and Pouring tasks. The Pouring
tasks are primarily designed to demonstrate the model’s
capability to generate multiple-solution trajectories andavoid obstacles under the task space setting. The Reach-tograsp task is intended to demonstrate the aforementioned
features under a joint space setting. Furthermore, thevideo is available at: https://sites.google.com/view/motiongeneration-on-manifolds/home.

5.1

Setup description

We consider a set of demonstrations involving a 7-DOF
Franka Emika Panda robot arm endowed with a twofinger gripper. The demonstrations were recorded usingkinesthetic teaching at a frequency of 10Hz. We calculate
geodesics on 100 × 100 and 50 × 50 × 50 graphs undertask space and joint space settings, respectively. All thelatent space graphs were bounded by different experimental
ranges based on the specific settings of the experiment.
For example, our experiments showed that the range of[−10.0, 10.0] is appropriate for the toy example setting forBeik-Mohammadi et al.
Latent graph

13

with Mµx,q (z) = Jµϕ (z)T MXQ (z)Jµϕ (z),MX (µϕ (z))
0
MXQ (z) =,0
MQ (µψ (z))

Decoded graphMσx (z) = Jσϕ (z)T MX (µxϕ (z))Jσϕ (z),Mκq (z) = Jκψ (z)T MQ (µψ (z))Jκψ (z).
Obstacle

Figure 8. Concept drawing. Left: The latent space isdiscretized to form a grid graph consisting of linearly spaced
nodes with edge weights matching Riemannian distances.
Right: To efficiently handle obstacles, the graph is decoded,
such that obstacles can easily be mapped to latent space.

all axes. Our Python implementation runs at 100Hz on
ordinary PC hardware. The approach readily runs in realtime. Additionally, obstacles are simulated along with adigital twin of the robot in a simulated environment toprovide real-time obstacle information (we did not integrate
obstacle localization systems in our setups).

where Jµϕ ∈ R(DX +DQ )×d is the Jacobian of the joint
decoder mean network (position and quaternion), andJσϕ ∈ RDX ×d and Jκψ ∈ RDQ ×d are the Jacobians of
the decoder variance and concentration RBF networks,
respectively. Since the position and quaternion share thesame decoder mean network, the output vector is split into
two parts, accordingly. The quaternion part of the decoder
mean is projected to the S3 to then define the corresponding
von Mises-Fischer distribution as in Eq. (15). The yellow
arrow in Figure 6-top shows the flow of the data fromambient space back to the latent space in order to computethe pullback metric.
The ELBO parameters β1 and β2 in Eq. (18) are foundexperimentally to guarantee good reconstruction of both
position and quaternion data. It is worth pointing outthat we manually provided antipodal quaternions during
training, which leads to better latent space structures andreconstruction accuracy. The same architecture is used for
all the experiments.
5.2.2

5.2

Architecture

In this section, we describe the VAE architecture underboth R3 × S3 and R7 settings. The VAE architectures are
implemented using PyTorch (Paszke et al. 2019).

5.2.1

VAE architecture in R3 × S3 : This particular VAE

network is designed to reconstruct end-effector poses in
R3 × S3 . The overall architecture is depicted in Figure 6top with the different components that are required tocorrectly reproduce end-effector poses. In this architecture,both the decoder and encoder networks have two hidden
layers with 200 and 100 neuron units (depicted as beige
boxes) which output the mean and variance vectors (depicted
as orange boxes). Furthermore, as previously explained,
the variance and concentration parameters for position andquaternion data are estimated using RBF decoder networks
with 500 kernels calculated by k-means over the training
dataset (Arvanitidis et al. 2018) with predefined bandwidth.
In this particular setup, the VAE uses a 2-dimensional
latent space (depicted by the green box) to encode the 7dimensional input vector (depicted as blue boxes on the left).
Moreover, the VAE reconstructs the encoded input back to
the ambient space (depicted as blue boxes on the right) using
the decoder networks.
In this setting, the VAE is implemented using a single
neural network as the decoder mean for the position andorientation, while the variance and concentration networks
are implemented separately. As a result, the learned metric isdefined as,
M (z) = Mµx,q (z) + Mσx (z) + Mκq (z),Prepared using sagej.cls(33)

VAE architecture in Rη : Here, we describe our VAE

network that reconstructs joint space movements in R7 . The
overall architecture is depicted in Figure 6–bottom withdifferent components. The input vector (depicted as the bluebox on the left) is a joint-value vector representing a single
configuration of the robot arm on a trajectory. This vector is
fed to the encoder network with two hidden layers of 200
and 100 neuron units (depicted as beige boxes) which are themean and variance vectors for the latent variable (depicted asorange boxes). Moreover, the variance RBF decoder network
uses 500 kernels calculated by k-means over the training
with predefined bandwidth.
Under the joint space setting, the VAE uses a 3dimensional latent space Z to encode the input vectors θ. As
usual, the decoder network reconstructs the encoded inputs
back to the joint space. However, in order to access thetask space information, necessary for whole-body obstacle
avoidance, our architecture is integrated with a forward
kinematics layer (depicted as the gray box). Note that
this layer is predefined based on the robot arm kinematic
model and does not change during training. We leverage
this layer to compute task space information regarding
multiple points on the robot (and not just the end-effector)
given the input configuration vector θ. To implement thiscomponent we used the Python Kinematic and Dynamic
Library PyKDL (Orocos 2021).
It is worth noting that as this VAE architecture uses
the forward kinematics layer during training, singularkinematic configurations need special attention. The main
problem arises in the formulation of our ELBO in
Eq. (26), which uses the volume measure computed as afunction of the determinant of the Jacobian of the forward
kinematics function. We can detect singularities whendet(JfFK (µα (zi ))) = 0, which may occur due to random
initialization of the VAE. To circumvent this issue and14

The International Journal of Robotics Research XX(X)

Figure 9. Left: The yellow
curves in the right cluster show
geodesics starting from thesame points and ending up at
random targets, and the blue
curves in the left cluster connect
random points on the manifold.
The background depicts themagnification factor derived
from the learned manifold, andthe semi-transparent white
points show the encoded
training dataset. Right: The
decoded geodesic employed on
the robot.

Figure 10. The evolution of the individual quaternion
components of the green geodesic as it passes throughhigh-energy regions.

guarantee that the learning process is not disrupted, we
simply add a small regularization term to the kinematic
Jacobian.
Finally, we evaluate the obstacle avoidance capabilities indifferent scenarios where the obstacles partially or entirely
obstruct the solutions in joint space. In these experiments,
the ambient metric of Eq.(32) is formulated by considering
all the joints on the robot, in addition to the end-effector.
This ensures that the robot will avoid obstacles as long as
they obstruct the solution for one or more joints or the endeffector. In other words, we do not consider points on the
robot links lying between joints, since it was not necessary
in our experiments. However, extra points on the robot links
can be easily added to guarantee a more robust obstacle
avoidance using the whole robot body.

5.3

Experiments in R3 × S3

The first set of experiments focuses on tasks where only the
robot end-effector motion is relevant for the task, therefore
the demonstrations are recorded in R3 × S3 .
Reach-to-grasp: The first set of experiments isbased on a dataset collected while an operator performs
kinesthetic demonstrations of a grasping skill. The grasping
motion includes a 90° rotation when approaching the object
for performing a side grasp (see Fig. 9). The demonstration
trajectories start from the same end-effector pose, and theyreach the same target position in the task space. To reproduce
this grasping skill, we computed a geodesic in Z which is
decoded to generate a continuous trajectory in X, whichclosely reproduces the rotation pattern observed during
demonstrations. Figure 9–left depicts the magnification
factor related to the learned manifold. The semi-transparent
white points correspond to the latent representation of the
training set, and the yellow and blue curves are geodesics
between points assigned from the start and endpoints of the
5.3.1

Prepared using sagej.cls

30

25

20

15

10

demonstrations. The left panel in Fig. 9 shows geodesics in
two different scenarios: The yellow geodesics in the rightcluster start from the same pose and end up at different
targets, while the blue geodesics in the left cluster start andend in different random targets. The results show that themethod can successfully generate geodesics that respect thegeometry of the manifold learned from the demonstration.
As expected, the magnification factor (Fig. 9–left) showsthat the learned manifold is composed of two similar clusters,
similar to the illustrative example in Fig. 3. We observedthat this behavior emerged due to the antipodal encodingof the quaternions, where each cluster represents one sideof the hyper-sphere. To confirm this, we generated a newgeodesic, depicted in green in the right panel of Fig. 9, which
is designed to cross the boundaries of the cluster (start andend locations belong to different clusters). Figure 10 depicts
the evolution of the quaternion components corresponding
to the decoded green geodesic. As this geodesic curve
crosses the clusters, the sign of the end-effector quaternion
flips (highlighted by the black rectangle). It is worthemphasizing that by staying on the manifold and avoiding
these boundaries, no post-processing of raw quaternion data
is required during training or reconstruction. This can beeasily guaranteed by monitoring the energy along the curves,which indicates when geodesics approach a high-energy
region. For instance, the average energy of the blue andyellow geodesics are 7.50 and 10.51, respectively, while thatof the green geodesic is 2.49 × 109 . As a result, we cansimply identify and avoid these geodesics.
Figure 9–right shows the reconstructed geodesic executedby the robot, where the overlapping images display the time
evolution of the skill. It is easy to observe that the desired
motion is successfully generated by our model. Note how theend-effector orientation evolves on the decoded geodesic in
the ambient space, showing that the 90° rotation is properly
encoded and reproduced using our approach.
Pouring: To evaluate our model on a more complicated scenario, we collected a dataset of pouring task
demonstrations. The task involves grasping 3 bottles from
3 different positions and then pouring 3 cups placed at 3
different locations. The demonstrated trajectories cross eachother, therefore providing a multiple-solution setting. As aresult, with 3 sets of demonstrations, all 9 permutations forgrasping any bottle from the table and then pouring at any
5.3.2

Beik-Mohammadi et al.

15
Figure 11. Left: Red and yellow
curves depict geodesics at two
different time frames performing
dynamic obstacle avoidance withobstacles depicted as red and yellow
circles. These geodesics also takeadvantage of a different group ofdemonstrations to generate a hybrid
solution. The background depicts themagnification factor derived from thelearned manifold with encodeddemonstration sets depicted as dot
clusters in white, green, and blue.
Right: The decoded geodesics
performed on the robot in different
time steps.

30.0

27.5

25.0

22.5

20.0

17.5

15.0

12.5

cup are feasible, despite only a small subset of them being
demonstrated.
The first feature we want to test in this setting is the
obstacle avoidance capabilities via metric reshaping. To
do so, we compute the ambient space metric in Eq. (30)based on a spherical obstacle that partially blocks the lowenergy regions that the geodesics exploit to find a solution.
This way the geodesics are forced to either use the lowenergy regions that the individual demonstrations provide
or improvise and find a hybrid novel path based on asubset of demonstrations. Figure 11–left shows the geodesic
performing obstacle avoidance around a moving obstacle
while following the geometry of the manifold. Two time
instances of the obstacle are depicted as red and yellow
circles in the latent space for illustration purposes. The
red and yellow curves represent geodesics avoiding the redand yellow obstacles, respectively. These curves correspondto one time frame of the adapted geodesic, showing how
our method can deal with dynamic obstacles. Fig. 11–rightshows the decoded geodesics performed on the robot, wheretransparent robot arms show the temporal evolution of theskill. In order to correctly perform obstacle avoidance, theparameter x in ambient space metric in Eq. (30) represents
the position of the bottle when grasped and the obstacle
radius r is modified to account for the bottle. To do so, we
simplify the bottle geometry by approximating it using a
sphere and add its radius rbot to the radius of the obstacle
as r = robs + rbot , where robs is the obstacle radius. This
prevents the bottle from colliding with the yellow and red
spheres that represent the obstacle.
Figure 12–left shows the ability to leverage multiplesolution tasks to generate novel movements emerging as
combinations of the observed demonstrations. Specifically,
we generate a combination of three geodesics that can beused to pour all the cups with a single bottle. The geodesic(depicted as a yellow curve) starts from an initial point(depicted in green) in the white demonstrations group at thetop, grasps the bottle, pours the first cup using the green
demonstration group, then goes to the second cup using theblue-demonstrations group, and finally pours the third cup
using the white-demonstration group again. Figure 12–rightuses the same visualization technique to show the decoded
geodesics employed on the robot by superimposing images
from different time frames. The task begins by pouring theright-side cup, then the center, and lastly the left-side cup. As
Prepared using sagej.cls

observed, our approach generates the geodesics that properly
switch among the demonstrated solutions to reproduce novel
movements in ambient space and successfully pour all thecups in the given order.
Multiple-limb obstacle avoidance: Multiple-limb
obstacle avoidance is a technique that aims to avoid obstacles
by considering the entire robot body, rather than just theend-effector. This approach offers a more comprehensive
way of obstacle avoidance as compared to techniques that
focus exclusively on end-effector motion. To implement thismethod, we leverage redundant solutions in joint space fromdemonstrations. Therefore, by operating in the joint space,
we can achieve redundancy at the joint level, which isessential for the success of multiple-limb obstacle avoidance.
When working under the task space setting, it is notpossible to achieve the redundancy in the joint level required
for multiple-limb obstacle avoidance. Thus, attempts toavoid obstacles using multiple limbs may be unsuccessful.
However, an off-the-shelf obstacle avoidance technique canbe used to navigate the robot away from obstacles andcomplete the task under the task space setting.
To evaluate this, we used the trained model from Sec. 5.3.1
to show how likely an obstacle avoidance technique fails toprovide an obstacle-free joint configuration. To do so, we
decoded a 150 × 150 equidistant latent grid to access theircorresponding input space states which represent the endeffector poses. Then, we used each pose to compute 100
different joint configurations using inverse kinematics. The
obstacle was placed in a way that was close to the robot bodybut did not block any demonstration in the task space, as
depicted in Fig. 18. Therefore, this setting requires multiplelimb obstacle avoidance in order to successfully performthe task. Each configuration was scored 1.0 even if a single
point on the robot collided with the obstacle, otherwise, thescore was 0.0. In Fig. 13, the background shows the sumof these scores for each latent state over the latent grid, andthe white points depict the encoded demonstrations in thelatent space. The results show that in the regions close todemonstrations, there is up to 90 percent chance to generatea joint configuration that leads to a collision.
5.3.3

5.4

Experiments in the joint space Rη

In this section, we focus on tasks where joint-level motion
patterns are relevant, and therefore the human teacher

16

The International Journal of Robotics Research XX(X)

Figure 12. Left: The geodesic
shown as the yellow curve combines
the blue, green, and whitedemonstration groups to form ahybrid solution. This experiment uses
three geodesics to pour all the cups
using a single bottle. Right: The
decoded geodesics performed on the
robot are depicted by superimposing
images from different time frames.
The transparent robot arms depict
the trace of the motion trajectory,
which begins with pouring the rightcup, the middle, and lastly the left
cup.
100

80

60

40

20

0

Figure 13. The background illustrates the frequency of failuresby an obstacle avoidance technique to provide obstacle-freetrajectories for the entire robot body. The white points depict theencoded demonstrations. The results indicate that the robot
configurations generated using inverse kinematics in the regions
close to the demonstrations are likely to result in collisions withsome regions exceeding 90 collisions in 100 trials.

provides kinesthetic demonstrations in Rη , where η is thenumber of DOF of the robot. When learning Riemannian
metrics in this setting, we initially designed a 2-dimensional
latent space for our VAE, which proved to be insufficient
to encode the demonstrated motion patterns. Specifically,
we analyze the capacity of the latent space to encode theskill manifold experimentally. To do so, we investigated theswitching behavior of geodesics in Z. In other words, overlap
between low-energy regions in Z, representing two or moredifferent demonstration sets, may lead to unnecessary andfrequent switches between these solutions when computing
a geodesic. To put it differently, when computing a geodesic,
switching between two demonstration sets is unavoidable
when the total energy of the geodesic switching between
them is less than the total energy without the switch.
However, switching behaviors may be avoided by having
high-energy regions among demonstration clusters in Z.
Furthermore, specifically under the joint space setting, thefrequent switching in geodesics may lead to jittery motion intask space when applied to the robot. To illustrate this issue,
we designed a simple experiment in which a robotic arm
follows a circular pattern with its end-effector. The start andPrepared using sagej.cls

target configuration of the geodesic is selected from the samedemonstrated trajectory, and the objective is to evaluate if thegeodesic stays on the low-energy regions corresponding to
the same demonstrated trajectory.
Fig. 14–left shows a geodesic curve computed in the
2-dimensional latent space depicted as the yellow curve.
This geodesic exhibits unnecessary switches among different
solutions (i.e. circular white demonstrations). Therefore,
when the decoded geodesic is deployed on the robot,
it results in jerky movements and undesirable back-andforth motions. To solve this issue, we evaluated the same
experiment using a 3-dimensional latent space. Figure 14–
right shows the magnification factor of the metric learned
using the same training data but in a 3–dimensional latent
space. This magnification factor shows several torus-likeclusters, representing separate demonstration groups insteadof collapsing them into a plane, as in the 2-dimensional case.
Moreover, the resulting geodesic (depicted as the yellow
curve) does not switch among solutions, which provides
stable and smooth robot end-effector movements when
decoded.
To provide further details on the learned manifold in the
3-dimensional latent space, we create an illustration shownin Fig. 15, where the learned metric is shown as yellow
hollow tubes. Their inner part encodes low-energy regions
which are surrounded by high-energy boundaries. Figure 15
displays a horizontal cut to show the inner part of the learned
metric. To illustrate the obstacles, the right hollow tube ispartially blocked by the red sphere representing the ambient
space metric in the latent space Z, which is a solid highenergy region. Additionally, the figure depicts a geodesic
curve traveling successfully along both tubes, showcasing acollision-free trajectory at the right-side plot.
Reach-to-grasp: Similar to the task space setting,
we used the reach-to-grasp task to evaluate our approach
under the joint space setting. In this case, the demonstrations
are quite similar at the end-effector level but differ at thejoint space, as we exploited the kinematic robot redundancy
to provide different joint trajectories. Figure 16–left shows
the magnification factor in a 3-dimensional latent space,
where it can be seen that the learned metric corresponds toseveral connected and separated hollow tubes. As mentioned
previously, these hollow tubes have low-energy inner regions
surrounded by high-energy boundaries, analogous to thelearned metric in a 2-dimensional latent space. As shownin the figure, the generated geodesics stay inside the tubes
5.4.1

Beik-Mohammadi et al.

17

Figure 14. Left: The
geodesics calculated in
2–dimensional latent space,depicted as the yellow curve,
reveal several unnecessarytransitions between different
solutions when connecting two
points in the samedemonstration. Right:
Geodesics computed in
3–dimensional latent space,shown by the yellow curve,
shows that the geodesic does
not switch between clusters.

Geodesic
Encoded Demonstrations

16

14

12

10

8

Figure 15. Concept drawing. Left: The hollow tube represents
the metric, with low-energy regions inside surrounded byhigh-energy boundaries. The geodesic depicted as the blue
curve travels through the low-energy regions. Right: The hollow
tube is partially blocked by a solid high-energy region
corresponding to the ambient space metric representing the
obstacle. The geodesic depicted as the blue curve successfully
travels through the low-energy regions meanwhile avoiding the
obstacle.

and avoid crossing the boundaries. Figure 17–left shows therobot executing the decoded geodesic by applying the joint
values directly on the robot using a joint position controller.
We can observe that the decoded geodesic is able to generate
the demonstrated grasp motion with 90° rotation during theapproaching part. Note that the start and end points of thegeodesics are extracted from the demonstrations.
Figure 16–middle displays the magnification factor wherean obstacle is placed in such a way that all possible solutions
in joint space are blocked (e.g. obstacle placed on thecommon target of all the demonstrations). Note that the
obstacle introduces a high energy zone (depicted in red)that passes through all of the hollow tubes representing
the learned manifold, therefore, none of the geodesics arefeasible. The feasibility of the geodesics is assessed byevaluating the total energy of the curve. If a geodesic
passes through an obstacle, its energy will suddenly increase,providing an indication of a potential collision. Additionally,
Fig. 16–right illustrates the same learned metric but reshapedusing a different ambient space metric. We can now see thatthe obstacle partially obstructs the possible solutions, andtherefore some few geodesics (depicted as blue curves) arestill able to successfully generate obstacle-free movements.
Figure 17–right shows the robot executing the decoded
geodesic using a joint position controller while avoiding the
obstacle.
Prepared using sagej.cls

We designed another experiment to showcase how themultiple-solution capabilities can be leveraged to generate
collision-free movements. If the different demonstrated joint
space trajectories sufficiently overlap, the learned manifold
will be characterized by several overlapping low-energy
regions (i.e. hollow tubes), which geodesic curves cantravel through. As a result, if the obstacles partially block
the learned manifold, the geodesic may still travel among
solutions to generate new movements out of combinationsof the provided demonstrations. To show this behavior, weused a different demonstration dataset where we accounted
for several overlapping solutions starting from the same
joint configuration. Figure 18–left shows the magnification
factor of the learned Riemannian metric in the 3–dimensional
latent space. The learned manifold and the ambient metriccan be distinguished visually based on their energy values.
The ambient metric, which represents the obstacle regions(depicted in red), encodes high-energy zones. The learned
manifold (depicted in yellow) is characterized by several
entangled hollow tubes. The geodesic, shown as the blue
curve, begins from the common start configuration on theleft and successfully navigates around the obstacle to reachthe target on the right side. Figure 18–right displays the
decoded geodesic deployed directly on the robot arm usinga joint position controller. Note how the robot arm avoids
the obstacle by carefully maneuvering around it while
successfully performing the grasping skill.
Pouring: To evaluate the method in a more complex
scenario, we also performed the pouring task under thejoint space setting. Similar to the task space experiment, thedemonstrations also overlap in joint space. For this specificexperiment, we mainly focus on evaluating the multiplesolution trajectories. Figure 19–left shows the magnification
factor in 3-dimensional latent space showing entangled
hollow tubes, which represent the learned Riemannian
metric. Each hollow tube encodes low-energy regions
surrounded by high-energy boundaries, the latter outlined
by green, red, and blue solid lines. Each color represents
one group of demonstrations. Since the magnification factor
learned with the original dataset of 5 demonstrations pergroup (corresponding to the bottle’s initial position) is very
hard to visualize, we opt for simplicity and used a subset
of 2 demonstrations per bottle. As shown, the geodesic(depicted as the black curve) starts from a point in thesecond demonstration group (green border) and switches tothe first group (red border) to reach the target located in
5.4.2

18

The International Journal of Robotics Research XX(X)

Figure 16. left: The magnification factor in 3–dimensional latent space contains hollow tubes representing the learned Riemannian
metric. These hollow tubes contain low-energy regions surrounded by high-energy boundaries. Several geodesics are calculated
and visualized in blue. The manifold’s low-energy areas are rendered transparent. middle: The same magnification factor when an
obstacle is introduced in such a manner that all viable solutions in the ambient space are blocked. Due to the fact that the obstacle
introduces a high energy zone (in red) that goes through all of these hollow tubes, none of the geodesics are feasible, therefore no
geodesic is shown in this plot. right: The same magnification factor when an obstacle is introduced in such a manner that partiallyobstructs the solutions in the ambient space. Several geodesics (blue curves) were left outside of the obstacle region on the left
side of the panel.

Figure 17. Left: The decoded
geodesic is employed on the robot arm
when no obstacle is present in theenvironment. Right: The decoded
geodesic employed on the robot arm
with an obstacle partially obstructing
the solutions. The obstacle is depicted
as the green sphere.

the third group (blue boundary). Figure 19–right shows thatthe decoded geodesic employed on the robot successfully
performs the pouring task. It is also evident the geodesic uses
a multiple-solution strategy to reproduce a new trajectorythat was not explicitly demonstrated to the robot in the
training phase.

6

Discussion

We presented a novel LfD approach that learns Riemannian
manifolds from human demonstrations, over which geodesic
curves are computed and used as a motion generator.
These geodesics are capable of recovering and synthesizing
robot movements similar to the provided demonstrations.
Therefore, our technique enables robot motion planning
between two arbitrary points on the learned manifold via
geodesics. This is achieved through a variational autoencoder(VAE) from which a pullback Riemannian metric can becomputed in the VAE latent space and used to generate
movements in both task and joint spaces. Additionally,
we proposed reshaping the learned Riemannian metricusing ambient metrics to avoid static or dynamic obstacles
on the fly. This metric reshaping along with a forward
kinematics layer enabled multiple-limb obstacle avoidance
capabilities, which allowed us to achieve reactive motions.
Prepared using sagej.cls

In order to guarantee fast geodesic computation and motion
generation, we developed a graph-based technique built on
the Dijkstra algorithm to provide real-time motion synthesis
and adaptation. We extensively tested our approach in bothjoint and task space settings, showing that our geodesic
motion generation performs effectively in both settings andin a variety of tasks such as grasping and pouring. We would
like to draw attention to the fact that the developed system
has certain limitations and potential areas for improvement.
In this regard, we would like to present a selection of suchaspects for discussion.
Geodesics generalization: During our experiments, weobserved that when the geodesic curves are forced to leave
the learned manifold, thus crossing high-energy boundaries,
the parts of the geodesic that lie outside the manifold
may cause inaccurate and undesirable motions. This is apotential issue when the given start or target points areplaced outside the data manifold, or when the obstacles
fully blocked it, which forces the geodesics to leave the
manifold to comply with the desired task specifications. We
believe this is a consequence of using VAEs to learn theskill manifold where the data lying outside the training data
support may be arbitrarily misrepresented in the latent space
Z. We hypothesize that this problem may be addressed by

Beik-Mohammadi et al.

19

Figure 18. Left: The geodesic depicted
as the blue curve avoids the obstacleby traveling through different
demonstrations. Right: The decoded
geodesic employed on the robot arm.
Images from different time steps aresuperimposed to show the trace of the
motion. The green sphere represents
the obstacle in the ambient space.

Figure 19. Left: The magnification
factor in 3–dimensional latent space
representing the learned Riemannian
metric. This metric contains low-energy
regions surrounded by high-energy
boundaries depicted as red hollow
tubes. The geodesic curve is depictedin black. The outline color of each tube
(green, red, or blue) indicates which
demonstration group it belongs to.
Right: The decoded geodesic
employed on the robot.

learning a bijective mapping between old demonstrations andnew conditions, and then using this function to transformthe learned manifold (e.g., by expanding or rotating) to fit
another region of the space.
Obstacle avoidance: Regarding our obstacle avoidance
approach, the ambient metric used to reshape the learned
metric enforces a “soft constraint” rather than a “hardconstraint”. This means, under certain situations, thegeodesic might still cross the obstacle instead of avoiding
it. It is worth noting that we were unable to replicate this
scenario since the generated geodesics tended to avoid the
obstacle by abandoning the manifold rather than crossing
the obstacle, which can be easily detected and prevented if
necessary. To address this potential issue, the graph-based
geodesic can be further exploited by removing the nodes
near the obstacle from the graph instead of re-weighting
the edges. While this strategy may reduce the computational
overhead and work in practice, it is not theoreticallygrounded. Finally, our obstacle avoidance formulation only
considered simple obstacles, but the strategy can be extended
to multiple dynamic obstacles. Instead of working withsingle balls, one can imagine extending the approach tocomplex obstacle shapes represented as point clouds. This
may increase the implementation demands in order to remain
real-time, but such an extension seems reasonable.
Latent space topology: As explained in Section 5, we

observed that the latent space dimensionality plays a criticalrole when learning joint space motions. As reported, a 3dimensional latent space was necessary to calculate smooth
Prepared using sagej.cls

geodesic curves for joint space tasks. This resulted inincreased computational complexity when calculating theambient metric and geodesics in the latent space. Thisimplies that this issue may exacerbate when working
with high-DOF robots, requiring higher-dimensional latent
spaces, and therefore raising the need to design more efficient
geodesic computation methods. Note that this phenomenon
may be related to the latent space topology, which is here
assumed to be Euclidean in accordance with the used VAE
Gaussian prior. We plan to theoretically and experimentally
analyze the effects of the latent space topology when learning
Riemannian manifolds for robot motion generation.

Scalability to high-dimensional latent spaces: In ourcurrent setup, we used task- and joint-space trajectories,
and therefore the latent space did not encode any external
perceptual information. In this setup, low-dimensional
latent spaces are particularly preferred due to their
reduced computational demands for the calculation of the
Riemannian metric and subsequent geodesics. However,higher dimensional latent spaces may also be of interest, as
they have a larger capacity to encode trajectories with highdimensional perceptual data (e.g. images). One question thatarises in this context is what the learned manifold or metricis representing, given that the latent space may combine
perceptual and motion data. Additionally, it is important
to consider whether disentangled latent representations arenecessary in such cases. Thus, the potential extension of ourwork involves learning a latent representation of task/joint

20

space trajectories and exploring the integration of perceptual
information into high-dimensional latent spaces.
The effective application of such extensions in higher
dimensional latent spaces hinges on the careful selection ofan appropriate real-time geodesic approximation and predictive variances modeling techniques. Our work has exploredvarious approaches for approximating geodesics, including
graph-based and gradient descent-based methods, which arepotential candidates for achieving this objective. In particular, our gradient-based approach exhibits scalability to highdimensional spaces, and existing literature provides evidence
supporting the scalability of graph-based approaches up to
20 dimensions (Chen et al. 2019). Moreover, our model
relies on RBF networks for modeling predictive variances,which are known to scale poorly to high-dimensional settings (Arvanitidis et al. 2018). There is thus a need for both
improving uncertainty estimation in generative models, aswell as geodesic computations, in order to robustly scale tohigh-dimensional latent spaces.
Latent space dimensionality: This paper presents our useof a 2-dimensional latent space in a task-space setting anda 3-dimensional latent space in a joint-space setting. The
choice of dimensionality for the latent space was basedon empirical observations, and as such, may vary acrossdifferent tasks and spaces. A more analytical approach may
involve determining the minimum number of dimensions
required to encode the data accurately. For instance, there
could exist a correlation between a system’s relative degree
and the dimensionality of the latent space. Specifically, asystem with a low relative degree is generally less nonlinear
and may be represented with fewer dimensions, while asystem with a high relative degree is more nonlinear andmay require a higher-dimensional latent space to accurately
capture its behavior. Nonetheless, the dimensionality of thelatent space is not solely dependent on a system’s relative
degree, as other factors such as task complexity, required
accuracy, and precision, amongst others, may also play a
relevant role.
Choice of distributions: Using an isotropic Gaussian kernel,
in which the covariance matrix is just a multiple ofthe identity matrix, is computationally simpler and more
efficient than using an anisotropic Gaussian kernel. This isbecause the isotropic version only requires the specificationof a single variance parameter, whereas an anisotropic
Gaussian kernel needs the specification of a full covariancematrix, which is more computationally costly and requires
more parameters to optimize. Additionally, the isotropic
Gaussian kernel may be more interpretable, as it has asimpler functional form and only depends on a single
variance parameter. However, the use of an isotropic
Gaussian kernel might compromise the model’s accuracy,
as it does not allow for capturing correlations between
different dimensions of the input space. Moreover, it isworth noting that the von Mises-Fischer distribution has been
leveraged as an approximation for data distribution on thehypersphere when encoding quaternion data. However, it
should be noted that there exist more expressive and accurate
distributions, such as the Bingham distribution, which may
provide more accurate modeling of the true distribution of
the observed quaternion data. In the future, we plan to
Prepared using sagej.cls

The International Journal of Robotics Research XX(X)

investigate the impact of using different distributions in theposition, orientation, and joint spaces.
Acknowledgements
This work was supported by a research grant (42062) from VILLUM
FONDEN . This project has also received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement 757360). Our
work was funded in part by the Novo Nordisk Foundation throughthe Center for Basic Machine Learning Research in Life Sciences(MLSS) with reference number NNF20OC0062606.

References
Aljalbout E, Chen J, Ritt K, Ulmer M and Haddadin S(2020) Learning vision-based reactive policies for obstacle avoidance.
In: Conference on Robot Learning(CoRL). URL https://proceedings.mlr.press/
v155/aljalbout21a.html.
Arvanitidis G, Hansen LK and Hauberg S (2018) Latent
space oddity: on the curvature of deep generative models.
In: International Conference on Learning Representations(ICLR). URL https://openreview.net/forum?id=
SJzRZ-WCZ.
Arvanitidis G, Hauberg S, Hennig P and Schober M (2019) Fast androbust shortest paths on manifolds learned from data. In: International Conference on Artificial Intelligence and Statistics(AISTATS). pp. 1506–1515. URL http://proceedings.
mlr.press/v89/arvanitidis19a.html.
Arvanitidis G, Hauberg S and Schölkopf B (2021) Geometricallyenriched latent spaces.
In: International Conferenceon Artificial Intelligence and Statistics (AISTATS).
URL
https://proceedings.mlr.press/v130/
arvanitidis21a.html.
Bahl S, Mukadam M, Gupta A and Pathak D (2020)
Neural dynamic policies for end-to-end sensorimotor
learning. In: Advances in Neural Information Processing
Systems (NeurIPS). pp. 5058–5069.
URL https://proceedings.neurips.cc/paper/2020/file/354ac345fd8c6d7ef634d9a8e3d47b83-Paper.
pdf.
Beik-Mohammadi H, Hauberg S, Arvanitidis G, Neumann G
and Rozo L (2021) Learning Riemannian manifolds for
geodesic motion skills. In: Robotics: Science and Systems(R:SS). URL https://roboticsconference.org/
2021/program/papers/082/index.html.
Bishop CM, Svensén M and Williams CK (1997) Magnification
factors for the SOM and GTM algorithms. In: Workshop
on Self-Organizing Maps. URL https://research.
aston.ac.uk/files/113377/NCRG_97_008.pdf.
Bullo F and Lewis AD (2004) Geometric Control of Mechanical
Systems: Modeling, Analysis, and Design for Simple Mechanical Control Systems. Texts in Applied Mathematics. Springer
New York. URL https://books.google.de/books?
id=jg7VumxOoe4C.
Calinon S (2016) A tutorial on task-parameterized movement learning and retrieval. Intelligent Service Robotics 9(1): 1–29. URL
https://doi.org/10.1007/s11370-015-0187-9.
Campa R and Camarillo K (2008) Unit quaternions: A
mathematical tool for modeling, path planning and control of

Beik-Mohammadi et al.

robot manipulators. In: Robot Manipulators, chapter 2. Rijeka:
IntechOpen. URL https://doi.org/10.5772/6197.
Chen N, Ferroni F, Klushyn A, Paraschos A, Bayer J andvan der Smagt P (2019) Fast approximate geodesics fordeep generative models.
In: International Conferenceon Artificial Neural Networks (ICANN). pp. 554–566.
URL https://link.springer.com/content/pdf/
10.10072F978-3-030-30484-3_45.pdf.
Chen N, Klushyn A, Kurle R, Jiang X, Bayer J and Smagt P(2018) Metrics for deep generative models. In: International
Conference on Artificial Intelligence and Statistics (AISTATS).
pp. 1540–1550.
URL http://proceedings.mlr.
press/v84/chen18e.html.
Cormen TH, Leiserson CE, Rivest RL and Stein C(2009) Introduction to algorithms.
MIT press.
URL
https://mitpress.mit.edu/books/
introduction-algorithms-third-edition.
Deisenroth MP, Faisal AA and Ong CS (2020) Mathematics for
Machine Learning. Cambridge University Press. URL
https://doi.org/10.1017/9781108679930.
Dollár P, Rabaud V and Belongie S (2007) Non-isometricmanifold learning: Analysis and an algorithm.
In:
International Conference on Machine Learning (ICML).
pp. 241–248.
DOI:10.1145/1273496.1273527.
URL
https://pdollar.github.io/files/papers/
DollarICML07manifold.pdf.
Eklund D and Hauberg S (2019) Expected path length on random
manifolds. In: arXiv preprint. URL https://arxiv.org/
abs/1908.07377.
Elbanhawi M and Simic M (2014) Sampling-based robot motion
planning: A review. IEEE Access 2: 56–77. URL https://ieeexplore.ieee.org/document/6722915.
Ewerton M, Neumann G, Lioutikov R, Ben Amor H, Peters
J and Maeda G (2015) Learning multiple collaborative
tasks with a mixture of interaction primitives. In: IEEE
International Conference on Robotics and Automation (ICRA).
pp. 1535–1542. URL https://doi.org/10.1109/
ICRA.2015.7139393.
Hauberg S (2019) Only bayes should learn a manifold. URL
https://arxiv.org/abs/1806.04994.
Hauberg S, Feragen A, Enficiaud R and Black MJ (2016)
Scalable robust principal component analysis using Grassmann
averages. IEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI) 38(11): 2298–2311. URL https://doi.org/10.1109/TPAMI.2015.2511743.
Havoutis I and Ramamoorthy S (2013) Motion planning andreactive control on learnt skill manifolds. International Journalof Robotics Research (IJRR) 32(9-10): 1120–1150. DOI:10.
1177/0278364913482016. URL https://doi.org/10.
1177/0278364913482016.
Hemingway EG and O’Reilly OM (2018) Perspectives
on Euler angle singularities, gimbal lock, and theorthogonality of applied forces and applied moments.
Multibody System Dynamics 44(1): 31–56.
DOI:
10.1007/s11044-018-9620-0.
URL https://doi.
org/10.1007/s11044-018-9620-0.
Higgins I, Matthey L, Pal A, Burgess CP, Glorot X, Botvinick MM,
Mohamed S and Lerchner A (2017) beta-vae: Learning basic
visual concepts with a constrained variational framework. In:
International Conference on Learning Representations. URL

Prepared using sagej.cls

21

https://openreview.net/forum?id=Sy2fzU9gl.
Huang Y, Rozo L, Silvério J and Caldwell DG (2019) Kernelized
movement primitives.
International Journal of Robotics
Research (IJRR) 38(7): 833–852. URL https://doi.
org/10.1177/0278364919846363.
Ijspeert AJ, Nakanishi J, Hoffmann H, Pastor P and Schaal S(2013) Dynamical movement primitives: Learning attractor
models for motor behaviors.
Neural Computation 25:
328–373.
URL https://homes.cs.washington.
edu/˜todorov/courses/amath579/reading/
DynamicPrimitives.pdf.
Jaquier N, Rozo L, Calinon S and Bürger M (2020) Bayesian
optimization meets Riemannian manifolds in robot learning.
In: Conference on Robot Learning (CoRL). pp. 233–
246. URL http://proceedings.mlr.press/v100/
jaquier20a.html.
Kalatzis D, Eklund D, Arvanitidis G and Hauberg S (2020)
Variational autoencoders with Riemannian Brownian
motion priors.
In: International Conference on
Machine Learning (ICML), volume 119. pp. 6789–6799.
URL
http://proceedings.mlr.press/v119/
kalatzis20a/kalatzis20a.pdf.
Khansari-Zadeh SM and Billard A (2011) Learning stable nonlinear
dynamical systems with Gaussian mixture models. IEEE
Transactions on Robotics 27(5): 943–957. URL https://ieeexplore.ieee.org/document/5953529.
Kingma DP and Welling M (2014) Auto-encoding variational
Bayes. In: International Conference on Learning Representations (ICLR). URL https://openreview.net/forum?
id=33X9fd2-9FyZd.
Kirk D (1970) Optimal Control Theory: An Introduction.
Networks series. Prentice-Hall.
ISBN 9780136380986.
URL https://store.doverpublications.com/
0486434842.html.
Klein H, Jaquier N, Meixner A and Asfour T (2022) A
Riemannian take on human motion analysis and retargeting.
In: IEEE/RSJ International Conference on Intelligent Robots
and Systems, IROS. IEEE, pp. 5210–5217. DOI:10.1109/
IROS47612.2022.9982127. URL https://doi.org/10.
1109/IROS47612.2022.9982127.
Konidaris G, Kuindersma S, Grupen R and Barto A (2012)
Robot learning from demonstration by constructing skill
trees. The International Journal of Robotics Research (IJRR)
31(3): 360–375. URL https://doi.org/10.1177/
0278364911428653.
Lee J (2018) Introduction to Riemannian Manifolds. 2 edition.
Springer. URL https://sites.math.washington.
edu/˜lee/Books/RM/.
Li M, Tahara K and Billard A (2018) Learning task manifolds
for constrained object manipulation. Autonomous Robots
42: 159–174.
URL https://doi.org/10.1007/
s10514-017-9643-z.
Manschitz S, Gienger M, Kober J and Peters J (2018) Mixtureof attractors: A novel movement primitive representation forlearning motor skills from demonstrations. IEEE Robotics
and Automation Letters (RA-L) 3(2): 926–933. URL https://ieeexplore.ieee.org/document/8255585.
Mohanan M and Salgoankar A (2018) A survey of robotic
motion planning in dynamic environments.
Robotics
and Autonomous Systems 100: 171–185.
URL

22

https://www.sciencedirect.com/science/
article/pii/S0921889017300313.
Moon H, Hoáng N, Robson N and Langari R (2012) Human arm
motion planning against a joint constraint. pp. 401–406. DOI:
10.1109/BioRob.2012.6290923. URL https://doi.org/
10.1109/BioRob.2012.6290923.
Noakes L and Zhang E (2022) Finding geodesics joining given
points. Advances in Computational Mathematics 48(4). URL
https://rdcu.be/dba9z.
Orocos (2021) orocos/orocos˙kinematics˙dynamics: Orocos kinematics and dynamics c library. URL https://github.
com/orocos/orocos_kinematics_dynamics.
Osa T and Ikemoto S (2019) Variational autoencoder trajectoryprimitives with continuous and discrete latent codes. CoRR
abs/1912.04063. URL https://link.springer.com/
article/10.1007/s42979-020-00324-7.
Osa T, Pajarinen J, Neumann G, Bagnell JA, Abbeel P and Peters
J (2018) An algorithmic perspective on imitation learning.
Foundations and Trends in Robotics 7(1–2): 1–179. DOI:10.
1561/2300000053. URL https://ieeexplore.ieee.
org/document/8620668.
Paraschos A, Daniel C, Peters J and Neumann G (2018) Using
probabilistic movement primitives in robotics. Autonomous
Robots 42: 529–551. URL https://link.springer.
com/article/10.1007/s10514-017-9648-7.
Paszke A, Gross S, Massa F, Lerer A, Bradbury J, ChananG, Killeen T, Lin Z, Gimelshein N, Antiga L, Desmaison
A, Kopf A, Yang E, DeVito Z, Raison M, Tejani A,
Chilamkurthy S, Steiner B, Fang L, Bai J and Chintala
S (2019) Pytorch: An imperative style, high-performance
deep learning library. In: Advances in Neural Information
Processing Systems (NeurIPS). pp. 8024–8035. URL https://proceedings.neurips.cc/paper/2019/file/
bdbca288fee7f92f2bfa9f7012727740-Paper.
pdf.
Peyré G, Péchaud M, Keriven R and Cohen LD (2010) Geodesic
methods in computer vision and graphics. Foundations andTrends in Computer Graphics and Vision 5(3–4): 197–397.
URL http://dx.doi.org/10.1561/0600000029.
Prescott T and Mayhew J (1992) Obstacle avoidance
through reinforcement learning.
In: Advances in
Neural Information Processing Systems (NeurIPS),
volume 4. Morgan-Kaufmann.
URL https://proceedings.neurips.cc/paper/1991/file/
9431c87f273e507e6040fcb07dcb4509-Paper.
pdf.
Rana MA, Li A, Fox D, Boots B, Ramos F and Ratliff N(2020) Euclideanizing flows: Diffeomorphic reduction forlearning stable dynamical systems. In: Proceedings of the
2nd Conference on Learning for Dynamics and Control,
Proceedings of Machine Learning Research, volume 120.
PMLR, pp. 630–639. URL https://proceedings.
mlr.press/v120/rana20a.html.
Ratliff N, Zucker M, Bagnell JAD and Srinivasa S (2009) CHOMP:
Gradient optimization techniques for efficient motion planning.
In: International Conference on Robotics and Automation(ICRA). pp. 489–494.
URL https://doi.org/10.
1109/ROBOT.2009.5152817.
Ratliff ND, Issac J, Kappler D, Birchfield S and Fox D (2018)
Riemannian motion policies. URL https://arxiv.org/

Prepared using sagej.cls

The International Journal of Robotics Research XX(X)abs/1801.02854.
Ravichandar H, Polydoros AS, Chernova S and Billard A (2020)
Recent Advances in Robot Learning from Demonstration.
Annual Review of Control, Robotics, and Autonomous Systems
3(1): 297–330.
URL https://doi.org/10.1146/
annurev-control-100819-063206.
Rozo L and Dave V (2021) Orientation probabilistic movement
primitives on Riemannian manifold. In: Conference on Robot
Learning (CoRL). URL https://proceedings.mlr.
press/v164/rozo22a.html.
Rozo L, Guo M, Kupcsik AG, Todescato M, Schillinger P,
Giftthaler M, Ochs M, Spies M, Waniek N, Kesper P
and Bürger M (2020) Learning and sequencing of objectcentric manipulation skills for industrial tasks. In: IEEE/RSJ
International Conference on Intelligent Robots and Systems(IROS). pp. 9072–9079. URL https://doi.org/10.
1109/IROS45743.2020.9341570.
Rozo L, Jimenez P and Torras C (2011) Robot learning fromdemonstration of force-based tasks with multiple solution
trajectories.
In: International Conference on Advanced
Robotics (ICAR). pp. 124–129. URL https://doi.org/
10.1109/ICAR.2011.6088633.
Scannell A, Ek CH and Richards A (2021) Trajectory optimisationin learned multimodal dynamical systems via latent-ode
collocation. In: IEEE International Conference on Robotics
and Automation (ICRA). pp. 12745–12751. URL https://ieeexplore.ieee.org/document/9561362.
Schaal S, Ijspeert A and Billard A (2003) Computational
approaches to motor learning by imitation. Phil. Trans. R.
Soc. Lond. B 358: 537–547. URL http://doi.org/10.
1098/rstb.2002.1258.
Seker MY, Imre M, Piater JH and Ugur E (2019) Conditional
neural movement primitives. In: Robotics: Science and Systems(R:SS). URL http://www.roboticsproceedings.
org/rss15/p71.html.
Shao H, Kumar A and Fletcher PT (2018) The Riemannian
geometry of deep generative models.
In: IEEE/CVF
Conference on Computer Vision and Pattern Recognition
Workshops (CVPRW). pp. 428–4288. URL https://doi.
org/10.1109/CVPRW.2018.00071.
Sra S (2018) Directional statistics in machine learning: A
brief review. In: Ley C and Verdebout T (eds.) Applied
Directional Statistics. URL https://arxiv.org/abs/
1605.00316. 1st edition.
Ugur E and Girgin H (2020) Compliant parametric dynamic
movement primitives. Robotica 38(3): 457–474. URL
https://www.cmpe.boun.edu.tr/˜emre/
papers/Ugur-2019-ROB.pdf.
Urain J, Anqi L, Puze L, D’eramo C and Peters J (2021)
Composable energy policies for reactive motion generation andreinforcement learning. In: Robotics: Science and Systems(R:SS). URL http://www.roboticsproceedings.
org/rss17/p052.html.
Zeestraten M (2018) Programming by Demonstration on Riemannian Manifolds. PhD thesis, University of Genova, Italy.
URL https://iris.unige.it/retrieve/handle/
11567/930621/245803/phdunige_4054447.pdf.