Xiang Zhang 11institutetext: University of Innsbruck, Technikerstr. 21a,
Innsbruck,
11email<EMAIL_ADDRESS>
WWW home page: https://iis.uibk.ac.at/

# Learning Descriptor of Constrained Task from Demonstration

Xiang Zhang 11 Matteo Saveriano 11 Justus Piater 11

###### Abstract

Constrained objects, such as doors and drawers are often complex and share a
similar structure in the human environment. A robot needs to interact
accurately with constrained objects to safely and successfully complete a
task. Learning from Demonstration offers an appropriate path to learn the
object structure of the demonstration for unknown objects for unknown tasks.
There is work that extracts the kinematic model from motion. However, the gap
remains when the robot faces a new object with a similar model but different
contexts, e.g. size, appearance, etc. In this paper, we propose a framework
that integrates all the information needed to learn a constrained motion from
a depth camera into a descriptor of the constrained task. The descriptor
consists of object information, grasping point model, constrained model, and
reference frame model. By associating constrained learning and reference frame
with the constrained object, we demonstrate that the robot can learn the book
opening model and parameter of the constraints from demonstration and
generalize to novel books.

###### keywords:

learning from demonstration, articulation model

##  1 INTRODUCTION

In our daily life, manipulating daily objects is complex tasks often have a
certain degree of constraints for utilities and space efficiency. For
instance, a door is rotating around an axle; drawers can be only pulled out in
one direction; a robot needs to slide a card in a card reader with a specific
orientation to get the card properly read. These tasks share a consistent
mathematical structure.

Nowadays, robots are capable of navigating in a complex environment and plan
their manipulation with obstacles. In the factory, the company is exploring
the possibilities of teaching a robot by human demonstration or coaching. [1]
The main technical obstacle is how robots could conveniently capture,
reproduce and generalize the complex manipulation on objects especially with
physical constraints

For a robot executing tasks, the generalization is applied in two-fold.
Firstly, the robot should be able to adapt itself to different environments,
e.g. obstacles. Secondly, the robot should be able to manipulate similar
objects with different states (e.g. poses) and sizes. More loosely speaking,
the robot should have generalizable skills for a dynamic environment.

The robot can capture and learn physical constraints from human demonstration
via two means. On the one hand, a robot should learn the property of the
object together with its affordance [2] with which the robot plans its
behavior during execution [3]. One the other hand, researchers approach this
problem by directly encoding the motion of human manipulating constrained
objects. And then adapt those problems by adapting hyperparameters, namely
task parameters. Traditional methods merely learn the trajectory from the
demonstration with the hope that the constraints are implicitly handled by the
model ([4]).

In one branch of research, the robot only learns the trajectory from
demonstration, the generalization depends on the task parameters which is
determined by the human or some hyper-parameter. In contrast, there are
algorithms designed only to obtain the object model, the returned model only
predicts the motion of the objects [5], or even the motion of the objects from
images [6]. Both kinds of methods are not sufficient for efficient
generalization. As for generalization, given a similar object, the object
model should not only detect the constraints but also return them as suitable
task parameters to the robot for executing the learned policy.

In this article, we introduce a framework with a defined task descriptor to
learn constrained motion from the demonstration. The learned task descriptor
consists of object recognition, grasping point, task parameter detection
model, constrained model, reference frame, and task parameters. The role of
the reference frame is chosen where the parameters of the learned object model
can be conclusively turned into executable end-effector trajectory. For
instance, if a robot is learning how to open a door, the descriptor is modeled
by a rotational constrained model. The rotation center is the center of the
reference frame. The robot needs to learn how to identify the center of the
reference frame and the parameter of the constrained motion, namely the
rotation axle and radius. As a result, the motion can be generalized for any
kind of doors.

This framework closes the gap between the learning object model and robot
execution from the demonstration. We propose a standard procedure that
integrates all necessary information, i.e. object recognition, grasping point
detection, constrained motions, and reference frame, for learning and
generalizing a constrained task model from the demonstration. By the
integrated information from the vision system, the ambiguity of the
demonstrated trajectory can be clarified by the element defined in our task
descriptor.

We take color and depth frames from human demonstration and learn the
descriptor of the constrained tasks. All of the components we are proposing is
for demonstrating the framework and can be replaced by more dedicated and
advanced methods in the future. We specify each of the components in the
following sections.

##  2 Related work

Learning from Demonstration has extensively focused on reproducing solely
demonstrated trajectories [7]. The constraints in the task are specified by
task parameters, such as coupling terms in DMPs [8]. Other methods learn
constraints implicitly by putting a restriction on the region with less
variation across multiple demonstrations.

The current research of vision assisted programming by demonstration ([9])
only focus on trajectory learning, human gesture learning, object recognition,
and robot navigation. Linking the vision, constraints and the resulting
trajectory is still an open question.

Modeling and estimating the motion of the constrained objects attract a lot of
research interest. Klingbeil [10] proposes a method that identifies a door
handle for dealing with this particular type of articulated object. More
sophisticated models are proposed describing and estimating the articulated
kinematic map from observed motions [5]. Following this method, the algorithm
is integrated into a vision processing algorithm so that the constrained model
is automatically extracted solely from a visual recording of the manipulation
articulated object [6]. Furthermore, an articulation model change detection
algorithm is proposed to handle a more complex and realistic situation, such
as moving a stapler and using a stapler. Huang 2012. As for faulty
articulation models, Huang 2015 proposes methods that use particle filter to
eliminate the wrong models and refine constrained object models.

Instead of solely rely on the perception, the constrained object model can be
obtained or refined by robot interactive perception. Nekum [11] aims to
extract a constrained object model in 3D space given the observation with
occlusion which can be clarified by robot interaction. Other methods discover
and learn the object model learning totally from scratch by interaction [12].
Martin [13] extends the idea to the extent that the algorithm only needs
visual input from the interaction.

There are also methods that cultivate robot applications. At the basic level,
the robotics programmer needs to handle the constrained motion with an
efficient task kinematic description and framework [14]. Furthermore,
controllers are designed based on the articulation models which update the
articulation models parameters online while interacting with objects [15].
Besides adaptive controllers, Dang Proposes a framework for manipulation
primitives, namely a position and a direction, where they define a task
descriptors for constrained object manipulations [16]. Otte integrates the
articulation object research into a framework that guides the robot to proper
exploration points and facilitates complex robot exploration [17]. Claudia
proposes a method that takes a constrained parameter and its reference frame
as the basis for planning multi-step manipulation tasks.

Some of the methods rely on other sources of information. Nagahama [18]
proposes a framework that takes grasping and manipulation into consideration.
With the provided constraints parameters from users, the objects from vision
went through a pipeline ending up with a robot manipulation. Daniele [19] uses
the visual demonstration to learn an articulated manipulation using language
processing and description for model selection.

##  3 Methodology

We argue that the manipulation skill should be defined and restricted to its
manipulated objects. Therefore, for every manipulation skill, the algorithm
should contain the object information (namely object identification and
segmentation). For manipulation, affordance, constrained model 3.1 and
reference frame determination are also necessary 3.2. The task descriptor of
the constrained task is as follows:

  * •

Identified Object: the name and location of the manipulated object

  * •

Grasping point and pose: grasping point for manipulation

  * •

Segmented objects: bounding box of the objects

  * •

Constrained model: the type of constraints and its corresponding parameters.

  * •

Reference frame: the reference frame describing the motion, depending on the
type of the constrained model. The reference frame is learned and obtained by
a model that map the image of the object to a rotation center and an axle in
the image.

We will show how the proposed framework works with simple algorithms as its
elements. The details about how each of the element is constructed from
demonstration and how they are recall and applied in the novel scene and
objects are illustrated in fig.1 and explained in the following sections.

Figure 1: The overall scheme

###  3.1 Kinematic Models Acquisition

Our framework depends on a depth camera. The learning procedure is shown on
the left side of Fig.1. The input from the camera is shown in Fig. 2. We then
use YOLO-net to identity potential objects in the scene. Then we obtain the
grasping point with a QR code. The manipulated object can be then labeled with
an object name. Consequently, we can segment the object out of the image with
the bounding box rendered by YoLo-net.

In order to learn the model of the object from the depth image, we need to get
the motion of the grasping point. Given the depth frames of the demonstration,
the robot tracked the position of the grasping points in the color/depth
($p_{x},p_{y},d$) image which is de-projected to a point cloud ($x,y,z$) for a
constrained motion.

The reason why we choose to work within the color image instead of point cloud
for rotation axle is that the point cloud is not always stable under different
light conditions, whereas project pixels from 2 camera mainly depends on the
camera intrinsic and relative stability to the environmental conditions. In
addition, a lot of patterns on a surface cannot be detected easily from the
point cloud but pretty obvious in the color image.

Subsequently, the constrained object can be described by the best kinematic
graph model $\hat{G}$ [5] given a sequence of pose observations $\mathcal{D}$.
The algorithm maximizes the posterior probability of observing the pose
$\mathcal{D}$ by $\hat{G}$. Our application keeps the assumption of a
kinematic tree.

$\hat{G}=\operatorname*{arg\,max}_{G}p(G|\mathcal{D})$ (1)

Under our assumption, the model is restricted to five types of articulated,
namely, resolute model, prismatic model, parallel model, and rigid motion
which indicate the demonstration is not a motion constrained on objects. The
optimization of $G$ consists of two-part, namely selecting the model type $M$
and learning the parameter $\theta$. Each type of model has a specific form of
a parameter vector that is used to construct a rotation matrix for
optimization. The optimization is done by MLESAC (maximum likelihood
consensus) [20].

###  3.2 Reference frame Model

The reference frame is crucial to learning and generalization of the
demonstrated motion. With the correct object model, the manipulation can be
transferred to a novel object by tuning parameter. In order to represent the
motion in a convenient and intuitive manner, we believe that the robot should
be able to determine the reference frame in which the motion to be tracked by
the robot is described and generated, for instance, the reference the rotation
center will reduce the motion to an equation determined by radius. We only
need to transform the motion back to the robot reference frame. Therefore our
framework needs to learn a model that maps the initial object information of
the manipulated objects to the location of the reference frame.

For instance, and the rotation center for opening a door is Based on these
structural constraints, it is possible to generalize the learned kinematic
model to the other similar objects.

The reference frame is defined by a center $c$ and a vector $v$. The key
component of this framework is to get the constrained object model through
visions and consequently obtain the reference frame for generalization and
execution. (Green Block in fig. 1) The parameters in the kinematic models have
the physical meaning which naturally corresponds to the shape of the observed
objects. The robots are able to learn the reference frame from kinematic
models from a demonstration by projecting the constrained model from section
3.1 to the image, which will be used to train mapping from the raw vision into
a rotation axle in the pixel coordinates. The reason why we choose the pixel
frame is that the 3D pose of the object is estimated from a 2D depth image. In
addition, some of the objects may not be detectable in the 3D spaces with low-
cost sensors, such as a door that is nicely sealed or a thin book on the
table.

We need to have a model taking the initial state of object images and grasping
as input and the rotation axle as a vector in images as output. For
demonstrations of the framework, from the initial state of a rotation object,
we propose a simple heuristic algorithm. For objects with a revolute joint,
the rotation axle $v$ is chosen either on the longer or shorter edge of the
object that is further away from the grasping point. The rotation center is
the perpendicular intersection point on the rotation axle through the grasping
point.

As for the prismatic constraints, the position and orientation of the grasping
point and the directional vector is sufficient to place the reference frame
during generalization. Similar heuristics can be defined for reference frame
placement. But this is not the focus of introducing the framework.

Figure 2: The input image (left) and depth frame (right) from our sensor. The
bounding box on the left is the result of object recognition YOLOnet-v3 [21]

###  3.3 Manipulation

The generalization of opening a door is demonstrated in the framework
illustrated in Fig. 1. The information is extracted in the same way as in
learning. The task parameter, i.e. rotation axle and center, is first
predicted from the image of the door and grasping point in the image which is
de-projected to the real world for robot execution.

As mentioned, the geometric information alone is not conducive to robot task
execution. In general, a motion-defined in the robot sensorimotor system can
be transformed into an end-effector in the robot world coordinates. Therefore
we define each type of articulation motion with respect to its identified
reference frame with the following rules.

For articulated rotation, we place the reference frame at the location of the
detected rotation center. The orientation of the reference frame is defined in
such a way that the Y axle is aligned with the rotation axle and the z axle is
aligned with the surface norm of the manipulated object. Similarly, the
reference frame of the prismatic motions has centered at the grasping point
and the sliding director is aligned with the frame orientation in the same
way. Having the reference frame, grasping point, and the type of the
constrained motion, the end-effector trajectory is calculated and transformed
into a robot execution frame.

As a result, through this framework, the perceived articulated model in the
demonstration generalized to novel books. The generalization heuristics can be
replaced by more advanced approaches that directly taking object image as
input.

##  4 Experiment

In this section, we illustrate our framework with a book-open use case. The
robot learns how to open a book from the demonstration and apply this learned
skill to novel books with different sizes. We do not use experiments on
prismatic motions since it is simpler than rotational motion and trivial for
the topic of this paper. We use Intel® RealSense depth camera to get visual
input. The camera is able to return the aligned color image, depth image, and
the resulting point cloud as shown in Figure 2. The sampling rate is 20 Hz.

For object information, we use YOLO-net-v3 [21] to obtain the name of the
manipulated objects in the scene during the demonstration. In this article, we
use pre-define QR code to indicate the position and orientation of the
grasping points since it is not the focus of our research. The encoded
information in the QR code is the orientation of the grasp which is similar to
the output of the grasp planning.

###  4.1 Learning from Demonstration

During the demonstration, the algorithm first recognizes the object and the
grasping point from the color image. The object where the grasping point is
attached is recognized by object classification as shown in Fig. 2. The robot
then adds the object name from YOLO-net into the object information. After
that, the initial object is obtained from the bounding box and saved for
learning later.

The trajectory of the grasping point is recorded in the 3D camera coordinates.
We apply the articulation model detection algorithm described in 3.1 on the
recorded trajectory and obtain rotation center and axle as shown in Figure 3.

Figure 3: Learning: The rotation axle obtained from the deprojected point
cloud. The blue dots in the image is the blue dots each of which is connected
to the rotation center with a red line and visualizes a red plane. The blue
line starts from the rotation center is the rotation axle.

After that, the rotation axle is projected back to the color frame, which is
shown in Fig. 4. According to the set heuristic rules in 3.2, the robot check
whether the rotational axle is on the longer edge and shorter edge. And its
distance to the grasping point. The robot adds them to the task descriptor for
generalization. In addition, a grasping model may be learned nevertheless in
our article it is replaced by QR code detection. With this information above
the robot can generate the reference frame for execution.

Figure 4: Learning: The rotation center points and the rotation axle projected
back to the initial state of the object from the point cloud. The blue dot is
the rotation center and (the line between two dots are the rotation axle)

###  4.2 Generalization

When the robot receives command for opening a book, the robot first look for
’books’ from object classification shown in Fig. 5.

Figure 5: The robot first find the book (the green bounding box) to be
manipulated in a scattered environment.

Then the grasping point and object image are sent to generate the reference
frame model, in this case, the robot looking for the long edge on the opposite
side of the grasping point. The axle is shown as the green line in Fig. 6.

Figure 6: Within the bounding box, the robot looks for the rotational axle and
grasping points and consequently the radius required for execution.

The predicted axle and rotation center in Fig. 6 is de-projected back to the
3D camera coordinates. The rotation center is the origin of the reference
frame and the Y axle is aligned with the rotation axle. Finally, the
trajectory for the new book to be tracked is obtained by applying circle
motion with the radius from grasping point to rotation axle. The result is
similar to as shown in Fig. 3 but with a different radius.

##  5 CONCLUSIONS

In this paper, we propose a framework for robot learning constrained-object
manipulation, which is based on the articulation learning model. Our framework
defines a descriptor of the constrained tasks that integrate all essential
information required to learn and execute constrained motion from the
demonstration.

The framework only on prismatic and rotational motions. The constrained motion
is imposed by the manipulated objects or tasks per se, the obstacle avoidance
should be applied after the constrained motion is determined. We also ignore
the constraints in tool handling because the tool is designed for humans, and
tools such as pliers should be provided to a robot as accessories.

The methods for the modules in our framework are simplified for the purpose of
illustration. The object recognition and grasping point detection can be
replaced with any method in their field [2]. In the future, we will replace
the heuristic method for identifying parameters from vision with an advanced
machine learning method that can be applied to more general situations.

## References

  * [1] M. Ikeda, S. Maddukuri, M. Hofmann, A. Pichler, X. Zhang, A. Polydoros, J. Piater, K. Winkler, K. Brenner, I. Harton, et al., “Flexrop-flexible, assistive robots for customized production,” in Proceedings of Austrian robotics workshop, pp. 53–58, 2018.
  * [2] P. Zech, S. Haller, S. R. Lakani, B. Ridge, E. Ugur, and J. Piater, “Computational models of affordance in robotics: a taxonomy and systematic classification,” Adaptive Behavior, vol. 25, no. 5, pp. 235–271, 2017\.
  * [3] C. Perez-D’Arpino and J. A. Shah, “C-LEARN: Learning geometric constraints from demonstrations for multi-step manipulation in shared autonomy,” Proceedings - IEEE International Conference on Robotics and Automation, vol. 3, no. c, pp. 4058–4065, 2017.
  * [4] A. J. Ijspeert, “Dynamical Movement Primitives: Learning Attractor Models for Motor Behaviors Auke,” Neural computation, vol. 25, no. 2, pp. 328–373, 2013.
  * [5] J. Sturm, C. Stachniss, and W. Burgard, “A probabilistic framework for learning kinematic models of articulated objects,” Journal of Artificial Intelligence Research, vol. 41, pp. 477–526, 2011.
  * [6] S. Pillai, M. Walter, and S. Teller, “Learning Articulated Motions From Visual Demonstration,” 2015.
  * [7] A. Billard, S. Calinon, R. Dillmann, and S. Schaal, “Robot Programming by Demonstration,” Robotics handbooks, pp. 1371–1394, 2008.
  * [8] Y. Zhou and T. Asfour, “Task-oriented generalization of dynamic movement primitive,” IEEE International Conference on Intelligent Robots and Systems, vol. 2017-Septe, pp. 3202–3209, 2017.
  * [9] H. Zhang, X. Lv, W. Leng, and X. Ma, “Recent Advances on Vision-Based Robot Learning by Demonstration,” Recent Patents on Mechanical Engineering, vol. 11, no. 4, pp. 269–284, 2018.
  * [10] E. Klingbeil, A. Saxena, and A. Y. Ng, “Learning to open new doors,” IEEE/RSJ 2010 International Conference on Intelligent Robots and Systems, IROS 2010 - Conference Proceedings, pp. 2751–2757, 2010.
  * [11] K. Hausman, S. Niekum, S. Osentoski, and G. S. Sukhatme, “Active articulation model estimation through interactive perception,” Proceedings - IEEE International Conference on Robotics and Automation, vol. 2015-June, no. June, pp. 3305–3312, 2015.
  * [12] D. Katz, A. Orthey, and O. Brock, “Interactive perception of articulated objects,” Springer Tracts in Advanced Robotics, vol. 79, pp. 301–315, 2014.
  * [13] R. Martin Martin and O. Brock, “Online interactive perception of articulated objects with multi-level recursive estimation based on task-specific priors,” IEEE International Conference on Intelligent Robots and Systems, pp. 2494–2501, 2014.
  * [14] J. De Schutter, T. De Laet, J. Rutgeerts, W. Decré, R. Smits, E. Aertbeliën, K. Claes, and H. Bruyninckx, “Constraint-based task specification and estimation for sensor-based robot systems in the presence of geometric uncertainty,” International Journal of Robotics Research, vol. 26, no. 5, pp. 433–455, 2007.
  * [15] A. Jain and C. C. Kemp, “Pulling open doors and drawers: Coordinating an omni-directional base and a compliant arm with equilibrium point control,” Proceedings - IEEE International Conference on Robotics and Automation, pp. 1807–1814, 2010.
  * [16] H. Dang and P. K. Allen, “Robot learning of everyday object manipulations via human demonstration,” IEEE/RSJ 2010 International Conference on Intelligent Robots and Systems, IROS 2010 - Conference Proceedings, pp. 1284–1289, 2010.
  * [17] S. Otte, J. Kulick, M. Toussaint, and O. Brock, “Entropy-based strategies for physical exploration of the environment’s degrees of freedom,” IEEE International Conference on Intelligent Robots and Systems, no. Iros, pp. 615–622, 2014.
  * [18] K. Nagahama, K. Takeshita, H. Yaguchi, K. Yamazaki, T. Yamamoto, and M. Inaba, “Estimating Door Shape and Manipulation Model for Daily Assistive Robots Based on the Integration of Visual and Touch Information,” IEEE International Conference on Intelligent Robots and Systems, no. 2, pp. 7660–7666, 2018.
  * [19] A. F. Daniele, T. M. Howard, and M. R. Walter, “A Multiview Approach to Learning Articulated Motion Models,” pp. 371–386, 2020.
  * [20] P. H. Torr and A. Zisserman, “Mlesac: A new robust estimator with application to estimating image geometry,” Computer vision and image understanding, vol. 78, no. 1, pp. 138–156, 2000.
  * [21] J. Redmon and A. Farhadi, “Yolo9000: better, faster, stronger,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 7263–7271, 2017.


[a]Tuan Nguyen

# Riemannian manifold hybrid Monte Carlo in lattice QCD

Peter Boyle  Norman Christ  Yong-Chull Jang  Chulwoo Jung

###### Abstract

Critical slowing down presents a critical obstacle to lattice QCD calculation
at the smaller lattice spacings made possible by Exascale computers. Inspired
by the concept of Fourier acceleration, we study a version of the Riemannian
Manifold HMC (RMHMC) algorithm in which the canonical mass term of the HMC
algorithm is replaced by a rational function of the SU(3) gauge covariant
Laplacian. We have developed a suite of tools using Chebyshev filters based on
the SU(3) gauge covariant Laplacian that provides the power spectra of both
the gauge and fermion forces and determines the spectral dependence of the
resulting RMHMC evolution of long- and short-distance QCD observables. These
tools can be used to optimize the RMHMC mass term and to monitor the resulting
acceleration mode-wise.

##  1 Introduction

Overcoming critical slowing down is a crucial objective of modern lattice
computations. The problem of critical slowing down stems from the increasing
range of distance or energy scales that must be accommodated in the usual
Markov-chain gauge-evolution algorithms as the continuum limit is approached.
This problem is made worse by the inherent non-renormalizability of the
favored Hybrid Monte Carlo (HMC) algorithm [1]. For free fields, we expect for
fixed physical lattice size but with decreasing lattice spacing $a$ that non-
local observables will decorrelate after integrating over a physical
trajectory length [2]. Because of the higher frequency modes the trajectory
step size must scale with the lattice spacing $a$, meaning the naive
autocorrelation time scales as $a^{-1}$. However, Lüscher and Schaefer present
evidence that the true autocorrelation time of these observables scales as
$a^{-2}$ [3] and suggest that this Langevin-like behavior results from the
non-renormalizability of the HMC algorithm.

The goal of Fourier acceleration is to modify the kinetic energy term in the
molecular dynamics Hamiltonian with a coordinate-dependent mass term such that
the low modes move at a higher MD velocity. In certain systems, a coordinate-
dependent mass term can be chosen such that all modes move at the same
velocity and critical slowing down is eliminated. For example, consider a
lattice pure $U(1)$ gauge theory in the continuum limit. Substituting the
gauge links with gauge fields via $U_{\mu}(x)=\exp(iA_{\mu}(x)/\sqrt{2\beta})$
and changing the MD mass term to be
$\mathcal{M}=(1-\kappa)-\kappa\partial^{2}$ the MD Hamiltonian becomes

$\mathcal{H}=\frac{1}{2}P_{\mu}\quantity[(1-\kappa)-\kappa\partial^{2}]^{-1}P_{\mu}+\frac{1}{2}A_{\mu}(\partial_{\mu}\partial_{\nu}-\delta_{\mu\nu}\partial^{2})A_{\nu}+\mathcal{O}(1/\beta)\,.$
(1)

Since the gauge fields are abelian, MD updates on the gauge links commute with
gauge transformations, and thus we are free to fix the system to Lorenz gauge.
Therefore the gauge modes are simply plane waves that evolve as
$\derivative*[2]{\tilde{A}_{\mu}(k)}{\tau}=-\omega^{2}(k)\tilde{A}_{\mu}(k)$
with characteristic frequency

$\omega^{2}(k)=\frac{k^{2}}{1-\kappa+\kappa k^{2}}\,.$ (2)

In the limit $\kappa\to 1$, all modes move with the same velocity and the
modified HMC is fully accelerated. [4] This same choice of mass term may also
accelerate $\phi^{4}$ theories in the perturbative limit since the eigenmodes
of the system are still approximately plane waves, meaning one can tune the
mode-dependence of the mass term to achieve approximate Fourier acceleration.
[5]

The situation is more complicated when we try to accelerate the HMC for
lattice QCD since the theory is non-abelian and non-perturbative. We can
parameterize the coordinate dependent mass term as a function of the SU(3)
gauge-covariant Laplace operator in order to maintain the gauge invariance of
the MD Hamiltonian. However, in order to use the eigenmodes of the covariant
Laplace operator to tune the mass term for Fourier acceleration we must show
that the eigenmodes of the theory have sufficient overlap with the covariant
Laplace eigenmodes. As we will show, when projected onto the subspace of
Laplace eigenvectors whose eigenvalues lie in a narrow interval
$[\lambda-\Delta\lambda/2,\lambda+\Delta\lambda/2]$ the size of the HMC force
that results from the gauge action is nearly linear with respect to $\lambda$.
Additionally, we will show that an extension of Fourier acceleration for non-
abelian pure gauge theories, the Riemannian Manifold HMC (RMHMC) algorithm, is
capable of moving certain long-distance observables faster than the naive HMC,
indicating some level of Fourier acceleration of the theory.

##  2 The HMC Gauge Force Power Spectrum

Figure 1: The power spectrum per mode of the gauge force for various lattice
actions is plotted with respect to the eigenvalues corresponding to each
eigenmode of the normalized covariant Laplace operator. The inverse lattice
spacings of each action are $a^{-1}=2\ \text{GeV}$ for the DBW2 action,
$a^{-1}=270\ \text{GeV}$ for the Wilson $\beta=10$ action, and $a^{-1}=3.6\
\text{GeV}$ for the Wilson $\beta=6.4$ action. We cut off the plot for very
low and very high eigenvalues of the normalized covariant Laplace operator
because there are too few covariant Laplace modes corresponding to these
eigenvalues to accurately measure the power spectrum.

In order to show that the eigenmodes of the covariant Laplace operator can be
used to Fourier accelerate the HMC algorithm, we calculate the spectral
dependence of the naive HMC gauge force per mode with respect to the
eigenvalues of the normalized covariant Laplace operator
$-\laplacian_{G}/(4d)$ where $\laplacian_{G}$ is defined in the adjoint
representation as

$\laplacian_{G}\phi(x)=\sum_{\mu}^{d}\quantity[U_{\mu}(x)\phi(x+\mu)U_{\mu}^{\dagger}(x)+U_{\mu}(x-\mu)^{\dagger}\phi(x-\mu)U_{\mu}^{\dagger}(x-\mu)-2\phi(x)]\,.$
(3)

We note that the covariant Laplace operator explicitly depends on the gauge
links, which can cause subtle complications that will be addressed later when
we present the Riemannian Manifold HMC algorithm.

We define a bandpass filter $\mathcal{B}(\lambda)\approx
P_{\lambda,\Delta\lambda}\quantity[-\laplacian_{G}/16]$ that approximates a
projection operator onto a narrow range of eigenvalues of $-\laplacian_{G}/16$
with a sum of Chebyshev polynomials: $\mathcal{B}(\lambda)=\sum
c_{i}T_{i}[-\laplacian_{G}/16]$. Using this bandpass filter, we project the
gauge force $F_{G}=-\partialderivative*{S[U]}{U_{\mu}(x)}$ (where
$\partialderivative*{U_{\mu}}$ represents the conventional derivative with
respect to the group parameters resulting in the left multiplication of the
link matrix $U_{\mu}$ by the group generators) onto a narrow interval of
eigenvalues near $\lambda$ to extract the contributions from covariant Laplace
eigenmodes corresponding to those eigenvalues. We calculate the power spectrum
of the gauge force per mode as

$\matrixelement{F_{G}}{\mathcal{B}(\lambda)}{F_{G}}\Big{/}\Tr[\mathcal{B}(\lambda)]\,.$
(4)

The gauge force power spectrum per mode is plotted in Figure 1 for various
lattice actions. The power spectrum is a relative measure of how quickly an
eigenmode of the covariant Laplace operator would move after integrating a
small trajectory time step. For the pure $U(1)$ gauge theory for example, the
power spectrum is a linear function of the wavenumber squared. From Figure 1
we see that power spectrum is surprisingly nearly linear with respect to the
eigenvalue of $-\laplacian_{G}/16$ for all actions, which confirms that the
eigenmodes of the gauge action are approximately covariant Laplace eigenmodes.

##  3 The RMHMC Algorithm

The Riemannian manifold HMC (RMHMC) is a class of Fourier accelerated HMC
algorithms where we add a coordinate dependent mass term that depends on the
covariant Laplace operator. The full Hamiltonian is given by

$\mathcal{H}=\dfrac{1}{2}\Tr\quantity{P_{\mu}^{\dagger}(x)\quantity[\mathcal{M}\quantity(-\laplacian_{G}/16)]^{-1}P_{\mu}(x)}+S\quantity[U_{\mu}(x)]+\log\absolutevalue{\mathcal{M}}$
(5)

The addition of the $\log\absolutevalue{\mathcal{M}}$ term is required in
order to cancel the additional force
$F_{M}=-\Tr\quantity{P^{\dagger}[\partialderivative*{\mathcal{M}^{-1}}{U_{\mu}(x)}]P}$
caused by the dependence of the mass term on the gauge links. This additional
term in the Hamiltonian will be realized by adding auxiliary fields and
auxiliary conjugate momenta. A further complication is the nonseparability of
the Hamiltonian, which requires an implicit leap-frog integration scheme.[6]

In Duane et al., they consider a simple mass term that is linear in the gauge-
invariant Laplacian. [4] Our initial studies of this simple mass term show
little if any reduction in autocorrelation. We therefore consider more general
sums of rational functions of the covariant Laplace operator to parameterize
our mass term:

$\mathcal{M}=\quantity[c+\sum_{i=1}^{N}\dfrac{a_{0,i}+a_{1,i}\quantity(-\laplacian_{G}/16)}{b_{0,i}+b_{1,i}\quantity(-\laplacian_{G}/16)+\quantity(-\laplacian_{G}/16)^{2}}]^{-2}\,.$
(6)

Figure 2:  The change in Wilson flowed energy (flow time of 8) after a single
HMC leapfrog integration step of $10^{-4}$ MD time units with injected
conjugate momenta projected onto specific covariant Laplace eigenmodes. All
algorithms have the largest change in Wilson flowed energy after injecting low
modes. All calculations are done with a DBW2 action with $a^{-1}=2\
\text{GeV}$ on a $16^{4}$ lattice.

For our current study we did not pursue a full-length autocorrelation study
between the naive HMC and the RMHMC. Instead, we investigate whether or not
the RMHMC is able to move long-distance observables further than the naive
HMC, allowing us to more quickly probe whether we have indeed Fourier
accelerated the theory. We chose to study changes in the Wilson flowed energy
between the two algorithms, since Wilson flowed observables are sensitive to
length-scales on the order of the square root of the flow time. [7]

As a preliminary test of the correlation between the Laplace eigenmodes with
small eigenvalues and the long-distance physics measured by the Wilson flowed
energy at large flow time, we examine how sensitive the Wilson flowed energy
is to perturbations to the gauge links by specific covariant Laplace modes. We
compare the naive HMC algorithm to an RMHMC algorithm with a mass term that is
hand-tuned to approximate the spectral dependence of the gauge force term
while remaining easy to invert. For the Wilson action, we consider an
additional RMHMC mass term that is fitted as a sum of many rational functions
to match the spectral dependence of the gauge force, which should
theoretically give optimal acceleration.

For each algorithm, we generate random conjugate momenta as if initiating an
HMC trajectory with a heatbath, but then use the bandpass filters from before
to project these momenta onto a narrow range of covariant Laplace eigenvalues.
Since the HMC heatbath depends on the mass term for the generation of random
conjugate momenta, the RMHMC algorithm will generate larger momenta
corresponding to low covariant Laplace eigenmodes. We then inject these
filtered momenta into a short HMC trajectory and measure the change in the
Wilson flowed energy. As shown in Figure 2, the Wilson flowed energy is pushed
further by the RMHMC algorithms compared to the naive HMC during a single,
short HMC integration step after injecting momenta corresponding to low modes,
which indicates that the covariant Laplace eigenmodes with small eigenvalues
are highly correlated with long-distance scales.

Figure 3: The change in Wilson flowed energy from that of the initial
configuration after a single HMC or RMHMC trajectory versus the number of MD
integration steps in that trajectory. The left plot is run with a DBW2 action
with $a^{-1}=2\ \text{GeV}$ on a $16^{4}$ lattice, while the right plot is run
with a Wilson action with $\beta=6.4$ and $a^{-1}=3.6\ \text{GeV}$ on a
$32^{4}$ lattice. Due to working on a larger lattice for $\beta=6.4$, we
double the flow time to 16 for the right plot.

Finally, we plot the change in the Wilson flowed energy after a single
trajectory as a function of the number of MD integration steps in that
trajectory for each algorithm to see if the RMHMC algorithms indeed change the
long-distance stucture of the gauge configurations faster than the naive HMC.
For each of our algorithms, we tune the MD trajectory step size so that all
algorithms give us similar accept/reject probabilities. Our results are shown
in Figure 3 for two actions, a DBW2 action with $a^{-1}=2\ \text{GeV}$ on a
$16^{4}$ lattice and a Wilson action with $\beta=6.4$ and $a^{-1}=3.6\
\text{GeV}$ on a $32^{4}$ lattice. For both actions, we see that the RMHMC
algorithms are able to move the Wilson flowed energy about twice as fast in
the first few MD evolution steps than the naive HMC, before slowing down to
about the same rate of change as the naive HMC. Thus, we conclude that we are
able to accelerate the low modes of the theory for the actions considered and
that the RMHMC is effective at moving long-distance observables faster than
the HMC at similar computation cost, i.e. with a similar number of integration
steps at fixed acceptance.

##  4 Conclusion and Next Steps

Figure 4: The power spectrum per mode of the fermion force for various Iwasaki
actions of different lattice sizes and coupling constants. Axes are the same
as in Figure 1

In this work, we propose the RMHMC class of algorithms for Fourier
accelerating the HMC algorithm by modifying the mass term in the MD
Hamiltonian such that the low modes of the theory can be moved more quickly.
By measuring the Wilson flowed energy, an observable that is sensitive to the
long-distance structure of the gauge configurations, we have shown that the
RMHMC algorithm is able to move these low modes in fewer evolutions steps than
the naive HMC. Ultimately, we are interested in using the RMHMC with fermions.
Applying the same techniques as before to measure the power spectrum of the
fermion force per mode, we show in Figure 4 that the fermion force has a
milder spectral dependence than the gauge force, and it already moves low
modes further than high modes. This is unsurprising since the fermion
determinant is a highly non-local operator. We expect that the RMHMC can be
applied to a hierarchical integration scheme where we separate the fermion
force and the gauge force onto different integration levels, where we would
only need the implicit integration scheme with the lower, gauge-only levels.

In future studies, we are interested in running a full-length study of the
autocorrelation time of some long-distance observables to show definitive
reduction in autocorrelation time. In addition, we would like to tune and run
the RMHMC algorithm with fermions. One potential first step in this direction
would be to investigate the spectral dependence of fermion forces from
individual Hasenbusch ratios.

## Acknowledgements

We thank our RBC and UKQCD Collaboration colleagues, in particular Joseph
Karpie and Robert Mawhinney, for helpful discussions and ideas. This research
was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative
effort of the U.S. Department of Energy Office of Science and the National
Nuclear Security Administration. We are grateful for the computational
resources provided by Argonne Leadership Computing Facility, which is a DOE
Office of Science User Facility supported under Contract DE-AC02-06CH11357 and
the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department
of Energy under Contract No. DE-AC05-00OR22725. In addition NHC was supported
by U.S. DOE grant #DE-SC0011941.

## References

  * [1] M. Lüscher and S. Schaefer, _Non-renormalizability of the HMC algorithm_ , _Journal of High Energy Physics_ 2011 (2011) 104.
  * [2] A.D. Kennedy and B. Pendleton, _Cost of the generalized hybrid Monte Carlo algorithm for free field theory_ , _Nucl. Phys. B_ 607 (2001) 456 [hep-lat/0008020].
  * [3] M. Lüscher and S. Schaefer, _Lattice qcd without topology barriers_ , _Journal of High Energy Physics_ 2011 (2011) .
  * [4] S. Duane, R. Kenway, B.J. Pendleton and D. Roweth, _Acceleration of gauge field dynamics_ , _Physics Letters B_ 176 (1986) 143.
  * [5] N.H. Christ and E.W. Wickenden, _Fourier acceleration, the hmc algorithm and renormalizability_ , 2018.
  * [6] M. Girolami and B. Calderhead, _Riemann manifold langevin and hamiltonian monte carlo methods_ , _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73 (2011) 123 [https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9868.2010.00765.x].
  * [7] M. Lüscher, _Properties and uses of the wilson flow in lattice qcd_ , _Journal of High Energy Physics_ 2010 (2010) 71.