# Learning Coordinated Terrain-Adaptive Locomotion by Imitating a Centroidal
Dynamics Planner

Philémon Brakel, Steven Bohez, Leonard Hasenclever, Nicolas Heess,
Konstantinos Bousmalis Deepmind
<EMAIL_ADDRESS>https://youtu.be/dvFwkvef2Ww

###### Abstract

Dynamic quadruped locomotion over challenging terrains with precise foot
placements is a hard problem for both optimal control methods and
Reinforcement Learning (RL). Non-linear solvers can produce coordinated
constraint satisfying motions, but often take too long to converge for online
application. RL methods can learn dynamic reactive controllers but require
carefully tuned shaping rewards to produce good gaits and can have trouble
discovering precise coordinated movements. Imitation learning circumvents this
problem and has been used with motion capture data to extract quadruped gaits
for flat terrains. However, it would be costly to acquire motion capture data
for a very large variety of terrains with height differences. In this work, we
combine the advantages of trajectory optimization and learning methods and
show that terrain adaptive controllers can be obtained by training policies to
imitate trajectories that have been planned over procedural terrains by a non-
linear solver. We show that the learned policies transfer to unseen terrains
and can be fine-tuned to dynamically traverse challenging terrains that
require precise foot placements and are very hard to solve with standard RL.

##  I Introduction

Robots with legs can access locations that are out of reach for platforms with
wheels and can carry loads that would be infeasible for airborne vehicles.
This makes legged locomotion an application area with practical impact. On top
of that, it is also a challenging research domain with complex control and
perception problems for which we know that solutions exist in nature.

While recent hardware and algorithmic developments have led to impressive
legged locomotion demonstrations of both learned controllers [1, 2], model-
based planning approaches [3, 4] and combinations of the two [5, 6], there
appears to be a trade-off between the overall agility of the systems and the
ability to traverse complex terrains that require precise foot placements.
Broadly speaking, methods that involve planning tend to struggle more with
dynamic movements while learned controllers struggle more with producing
precise coordinated foot placements and desirable gaits.

For many planning-based locomotion controllers, the lack of dynamic movements
is attributable to the necessity for the planner to provide trajectory plans
fast enough to run online. This requirement often rules out general non-linear
solvers and forces the search space to be constrained to more conservative
motions or the foothold planning and the dynamics planning to happen
separately. When footholds are planned independently of the robot dynamics,
they also have to be more conservative to guarantee feasibility, especially on
rough terrains.

Reinforcement learning (RL) methods generate closed-loop control policies
which can produce very dynamic and robust gaits that have been successfully
transferred to hardware [2, 7]. However, it often requires significant tuning
and the use of many shaping rewards to produce desirable gaits. For trajectory
optimization this is less of an issue as constraints can easily be added and
tuned relatively quickly. Furthermore, it can be hard for learning algorithms
to discover very precise movements over long time horizons via random
exploration. For these reasons, many successful learning methods employ
planning methods to provide more effective exploration [8, 9, 6]. A downside
of such hybrid approaches is that even if the planner doesn’t need to run
online during hardware deployment, the rate at which it can generate
trajectories can quickly become a bottleneck in the learning iteration. This
still limits most use cases of planners in learning systems to conservative
trajectory spaces or relatively short time horizons compared to state-of-the-
art trajectory optimization.

One way to circumvent exploration problems and also prevent undesirable gaits,
is through imitation learning. Motion capture data from a dog has been shown
to be a useful source of behaviors from which natural looking gaits were
successfully transferred to a quadruped robot [10]. However, obtaining motion
capture data for a wide variety of terrains would be very costly. It is also
not clear how much data would be needed and whether imitation learning
provides the required generalization to traverse unseen terrains.

In this paper, we demonstrate that imitation learning from planner data can be
used to efficiently create general terrain-adaptive controllers with precise
foot placements. We first use a non-linear centroidal dynamics planning
formulation in an offline fashion to generate a dataset of trajectories that
traverse a diverse set of procedurally generated terrains. Subsequently, we
use reinforcement learning to obtain policies that imitate these trajectories.
This setup has two advantages over conventional optimal control and RL
methods: first, the trajectory optimizer can take as much time as is needed to
find long time horizon trajectories over complex terrains because it doesn’t
need to operate in real time; second, the reinforcement learning agent does
not need to discover the coordinated behaviors and constrained gaits through
exploration with shaping rewards. Concretely, we make the following
contributions:

  1. 1.

We provide a recipe for generating varied terrain traversal datasets using a
centroidal dynamics planner.

  2. 2.

We demonstrate how to use RL to effectively imitate the generated
trajectories.

  3. 3.

We show empirically that the learned policies display reasonable zero-shot
performance on unseen terrains and that they can be fine-tuned to traverse
complex stepping stones and mixed procedural terrains that are hard to learn
from scratch.

##  II Related Work

###  II-A RL for Terrain-Aware Locomotion

RL has been used to learn terrain-adaptive policies for various simulated
characters/robots [11, 12]. In recent work [13], a multi-task setup was used
to learn a quadruped controller that was able to traverse stairs and various
other rough terrains in simulation. This controller used a pattern generator
modulating action space and it is not clear if the system could learn the
precise movements needed for the environments in our work. Similarly, in work
concurrent to ours, evolved pattern generators were combined with RL to
traverse challenging terrains [14]. In the ALLSTEPS approach [15], RL was used
to train policies that were able to walk over stepping stones requiring
precise foot placements. The behaviors were learned using a curriculum, highly
task specific rewards and feature observations based on step-by-step target
foothold locations. Other types of curricula have been used together with
shaping rewards to learn blind [1] and very recently also depth perceptive
[16] terrain adaptive controllers via RL. It is not clear whether such
curricula can easily be designed for any type of terrain.

###  II-B Combining RL and Planning for Locomotion

Various systems combine ideas from RL and optimal control for locomotion.
DeepGait [5], for example, is a hierarchical system in which a higher level
policy produces goal poses that are judged for feasibility by a centroidal
dynamics planner. This system can traverse challenging terrains but its
motions tend to be conservative and fairly static. This is probably due to the
simplified (convex) formulation that was used to make the planner fast enough.
In MPC-Net [6], trajectories are generated using an MPC controller and
optimized with a loss that minimizes the original control Hamiltonian. One
advantage of this method is that it takes the original optimal control
constraints into account. However, the learned policy is limited by what the
teacher MPC controller can do and the teacher still needs to employ an optimal
control solver in its inner-loop. RLOC [17] combines RL and optimal control by
learning a terrain-aware foothold planner that interfaces with a whole-body
controller (WBC). This systems was shown to traverse a wide variety of
terrains. Unlike our approach, RLOC still uses a WBC which may limit the types
of movements that can be selected by the perceptive foothold setting policy.
GLiDE [18] learns an RL policy for the centroidal dynamics directly. This
system was shown to traverse narrow beams and small stepping stones. However,
the foothold planning was done with terrain specific hand-designed heuristics.
Pre-planned trajectories have been used to guide learning of policies [19] to
aid exploration and constraint satisfaction but so far this method has only
been applied to flat terrains. Finally, learning methods have also been
combined with motion planning to speed up the planning process by training
neural networks to provide good initializations for the solver [20, 21]. The
obtained plans still need to be executed by a WBC but such methods are
complementary to our approach and could be used to speed up the dataset
construction.

###  II-C Planner Imitation for Legged Robots

In very recent work [22] that is perhaps most similar to ours, a single
trajectory was generated with a motion planner and RL was used to imitate it
in a timing independent way and learn a bounding or hopping behavior which was
fine-tuned and evaluated on a robot. Our work is different in that we generate
an entire dataset of trajectories and use it for learning terrain adaptive
policies that can generalize and be fine-tuned for a variety of challenging
terrains.

###  II-D Character Animation and Mocap Tracking

Our work is inspired by recent methods for animated character control via RL-
based imitation of motion capture data [23, 24, 25]. Motion capture data has
also been used to transfer dog behaviors to a real quadruped platform [10] or
to provide adversarially learned motion priors [26]. Holden et al. [27] use
terrains that were fitted to existing motion capture trajectories for
kinematic replay. This is similar in spirit to our method for distorting
terrains but used in a different context.

##  III Method

As Fig. 1 shows, our method can be divided into three phases: data generation,
imitation learning and evaluation/fine-tuning. To generate the data, we first
generated procedural terrains as detailed in Sec. III-A, and defined a
feasible planner formulation which is described in Sec. III-B. In Sec. III-C,
we detail how this formulation was used to produce diverse trajectories that
traverse the terrains. To explain how imitation learning was used to reproduce
the trajectories with a feedback policy, we describe the RL tracking task in
Sec. III-D, the learning algorithm in Sec. III-E, the observations and network
architecture in Sec. III-F and finally the terrain distortions in Sec. III-G.
In Sec. III-H we describe how evaluation and fine-tuning were done for new
terrains.

Figure 1: A schematic overview of our method that can be separated in a _data
generation_ step in which the terrains and planned trajectories are generated,
an _imitation learning_ step in which the terrains are distorted and policies
are trained to imitate the trajectories and an _evaluation/fine-tuning_ step
where the input commands are changed to produce straight walking behavior and
there is potential further training for new types of terrains.

###  III-A Procedural Terrains

To obtain a diverse set of trajectories with variations in contact positions
and orientations, we implemented a system for procedural terrain generation.
We represented elevation data as a triangulated grid. We aimed to have
terrains with both flat and tilted surfaces. To do this, we created random
grids of tiles with differing heights and subsequently scaled these terrains
by multiplying with Gaussian functions representing round bumps. The height
fields were scaled down to always have a maximum height of
$30\text{\,}\mathrm{c}\mathrm{m}$. See Fig. 2 for examples of height fields
generated with this procedure and Appendix A for a more detailed description.

Figure 2: Examples of procedural terrains provided to the planner and used in
some of the evaluation environments.

###  III-B Planner Formulation

Our motion planner formulation extends the centroidal dynamics planner
proposed by Winkler et al. [28] and is implemented on top of their open-source
library (TOWR) which also provides an interface for the IPOpt solver [29] that
we used for optimization. The default formulation in TOWR is a nonlinear
program (NLP) of the _simultaneous_ type where both the states and controls
are optimization variables and the dynamics are represented as constraints.
The dynamics constraints describe the motion of a single rigid body (the
center of mass of the robot) and are a function of the contact forces and
paths of the end-effectors (feet). The motions and paths are represented as
splines of cubic Hermite polynomials. Optimization is done over the parameters
of these polynomials, which can be interpreted as the positions and velocities
of the spline nodes and the duration of each spline segment.

More precisely, the formulation defines a search over the CoM positions
$\mathbf{r}(t)\in\mathbb{R}^{3}$, the base orientations
$\mathbf{\theta}(t)\in\mathbb{R}^{3}$, the end-effector positions
$\mathbf{p}_{i}\in\mathbb{R}^{3}$, force profiles $\mathbf{f}_{i}$, and the
phase durations for each of the end-effectors $\Delta T_{i,j}$, where $i$
indexes the end-effector and $j$ the phase. The forces are only non-zero
during the contact phases. The optimization has a trivial cost but the
solution must satisfy a set of constraints that can be grouped in the
following way:

  * •

Dynamics of the single body system enforced at fixed time intervals.

  * •

Contact force constraints to ensure positive and bounded values in the terrain
surface normal direction and adherence to the friction pyramids.

  * •

Swing height constraints that ensure that the position spline nodes are above
the terrain.

  * •

Kinematic limit constraints that restrict the end-effector positions to
bounding boxes at fixed offsets of the center of mass.

We refer to the work by Winkler et al. [28] for the details of the formulation
and will only elaborate on our changes.

Similar to previous work [20], we modified the formulation to improve the
physical feasibility of the optimized trajectories. Since the force
constraints are by default only enforced at the nodes of the end-effector
force splines, it is possible for the end-effectors in contact phase to
generate forces that are large or even negative in the normal direction in-
between two node values. Similarly, the swing height constraints on the end-
effector positions are only enforced at the node positions and we observed
end-effectors passing through the terrains. To address this, we also enforced
the force and swing height constraints at fixed time intervals ($\delta
t=0.08s$ and $\delta t=0.04s$ respectively). These additional constraints made
the optimization problem harder and we increased the number of nodes per
spline segment for the end-effector position trajectories from 2 to 3 to
compensate and provide additional flexibility. This is different from the
approach by Melon et al. [20], where different constraints and also _cost_
terms where added to minimize ground-reaction forces and to encourage higher
swing-phase foot positions. Since our goal was to generate a _diverse_ dataset
and not to directly execute the trajectories with a WBC, our priorities were
different from Melon et al. and we opted for a more permissive formulation
without any cost terms. We think that cost terms are more likely to reduce the
diversity of the optimization results compared to a formulation where the
optimization landscape is flat within the region where the constraints are
satisfied.

###  III-C Data Generation

To create a dataset of trajectories and terrains, the terrain generation
routine from Section III-A and the trajectory solver are called repeatedly.
While the randomness of the terrains should lead to a diverse set of
trajectories, we also added some randomization to the initialization of the
solver. The end-effector and CoM spline node positions were initialized using
the procedure from TOWR using one of the predefined gait patterns (in TOWR
denoted as ‘C1’, which defines a fly trot pattern111A fly trot is probably not
the best gait for scaling complex terrains but it appeared to be a good
initial setting for the solver.), to which we added normally distributed noise
with a standard deviation of $\sigma=$10\text{\,}\mathrm{c}\mathrm{m}$$.
Similarly, we added noise with a standard deviation of
$\sigma=$5\text{\,}\mathrm{N}$$ to the initial values of the spline nodes of
the force variables. We set the goal location to be $2.3\text{\,}\mathrm{m}$
in front of the robot (based on the terrain size and optimization time) and
set the amount of time allowed to reach the goal to $4.6\text{\,}\mathrm{s}$
to promote an average forward velocity of
$0.5\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$. To give the solver enough time
for difficult trajectories but to also prevent problematic optimization
regions from halting progress, we set the maximum CPU time of the solver to 30
minutes. The result of the optimization are the splines describing the CoM and
end-effector motions and the durations of the swing and stance phases. These
were sampled into discrete values at a rate of
$100\text{\,}\mathrm{H}\mathrm{z}$.

To compute joint states (which were used for initialization) from the end-
effector and CoM trajectories, we used a geometric approach based on the
method implemented in the ANYbotics codebase. To obtain a consistent mapping,
we constrained the knees to always bend inwards. For the ANYmal B robot model
that we used, this corresponds to the so called x-configuration of the legs.
See [30] for an example of how to compute these quantities for a similar
quadruped.

The final dataset has $8926$ clips of $4.6\text{\,}\mathrm{s}$ each,
corresponding to about $11.4$ hours in total. Each clip $C_{n}$ was
represented by a tuple
$(\mathcal{M}^{n},\\{\mathbf{p}^{n}_{\text{com}}(t),\mathbf{v}^{n}(t),\mathbf{\omega}^{n}(t),\mathbf{p}_{i}(t),g^{n}_{i}(t),\mathbf{q}(t),\mathbf{\dot{q}}(t)\colon
0\leq t\leq T\\})$ combining a terrain image and discretized trajectory
information, where $\mathcal{M}$ is the $16\times 16$ terrain image,
$\mathbf{p}^{n}_{\text{com}}(t)$ the center-of-mass position,
$\mathbf{v}^{n}(t)$ the linear CoM velocity, $\mathbf{\omega}^{n}(t)$ the
angular CoM velocity, $\mathbf{p}_{i}(t)$ the end-effector positions,
$g^{n}_{i}(t)$ is a binary indicator of whether each end-effector $i$ was in
contact with the ground or not and $\mathbf{q}(t)$ and $\mathbf{\dot{q}}(t)$
are the joint positions and velocities obtained by inverse kinematics and
finite differences. Fig. 3(a) shows the distribution of the end-effector
contacts in the world frame $xy$-plane. The positions are clearly spread out
and indicative of a diverse set of trajectories. The forward velocity of the
trajectories varied significantly over time and among different trajectories
as can be seen in Fig. 3(b).

(a)

(b)

Figure 3:  Distribution plots of the dataset. Fig. 3(a) shows the planar end-
effector contact distribution data in the world frame, where $x$ is aligned
with the robot facing forward, up to $1.75$ meters (the end-position was
always the same so we left it out of the distribution). Fig. 3(b) shows the
velocity of the center of mass projected on the $x$-axis for 20 randomly
chosen trajectories as a function of time.

###  III-D Tracking Task

To train a policy to track the planned trajectories, we used setup similar to
the multi-clip motion capture tracking task from Hasenclever et al. [25]. Each
episode, a clip was chosen from the dataset and a random starting frame was
selected to initialize the simulated robot state. The episode ran until the
end of the clip was reached or a termination criterion had been triggered. The
agent was rewarded based on the similarity of the simulated trajectory to the
reference clip.

We used the same truncation criterion as Hasenclever et al. [25] to terminate
episodes when the robot would stray too far from the reference trajectory.
This criterion is defined as

$\displaystyle\epsilon$
$\displaystyle=\frac{1}{3|\mathcal{B}|}\sum_{i\in\mathcal{B}}\|\mathbf{p}_{i}-\mathbf{p}^{\text{ref}}_{i}\|_{1}+\frac{1}{|\mathcal{J}|}\sum_{j\in\mathcal{J}}\|q_{j}-q^{\text{ref}}_{j}\|_{1},$
(1) $\displaystyle r_{\text{trunc}}$ $\displaystyle=1-\frac{\epsilon}{\tau},$
(2)

where the vectors $\mathbf{p}_{i}$ represent the positions of all bodies and
the values $q_{j}$ are the joint positions with indices $\mathcal{J}$. An
episode was terminated when $r_{\text{trunc}}>\tau$, with $\tau=0.5$ for all
our experiments. Episodes were also terminated on contacts that were not
between the end-effectors and the environment.

The tracking reward consists of multiple objectives. The total reward used to
monitor tracking performance consists of five terms defined as

$\displaystyle r_{\text{com}}$
$\displaystyle=0.2\cdot\exp(-80\|\mathbf{p}_{\text{com}}-\mathbf{p}^{\text{ref}}_{\text{com}}\|^{2}),$
(3) $\displaystyle r_{\text{ee}}$
$\displaystyle=0.2\cdot\exp\left(-80\sum_{i}\|\mathbf{p}_{\text{ee,i}}-\mathbf{p}^{\text{ref}}_{\text{ee,i}}\|^{2}\right),$
(4) $\displaystyle r_{\text{linvel}}$
$\displaystyle=0.2\cdot\exp(-10\|\mathbf{\dot{p}}_{\text{com}}-\mathbf{\dot{p}}^{\text{ref}}_{\text{com}}\|^{2}),$
(5) $\displaystyle r_{\text{angvel}}$
$\displaystyle=0.2\cdot\exp(-10\|\mathbf{\omega}_{\text{com}}-\mathbf{\omega}^{\text{ref}}_{\text{com}}\|^{2}),$
(6) $\displaystyle r_{\text{quat}}$
$\displaystyle=0.2\cdot\exp(-2\|\mathbf{q}_{\text{com}}\ominus\mathbf{q}^{\text{ref}}_{\text{com}}\|^{2}),$
(7)

where $\mathbf{p}_{\text{com}}$ and $\mathbf{p}_{\text{ee,i}}$ are the
positions of the center of mass and the end-effectors,
$\mathbf{\dot{p}}_{\text{com}}$ and $\mathbf{\omega}_{\text{com}}$ are the
linear and angular velocities of the center of mass and
$\mathbf{q}_{\text{com}}$ is the quaternion representing the base orientation
at the center of mass.

###  III-E Learning Algorithm

Our learning problem can be described as a multi-objective POMPD where there
is more than one reward function and we have only access to observations
rather than complete states. To optimize the set of objectives, we used the
Multi-Objective variant of the VMPO algorithm (MO-VMPO) [31, 11], which has
been shown to work well for motion capture tracking. MO-VMPO is an on-policy
actor-critic policy improvement method, suitable for learning with neural
networks. Like other actor-critic methods, there is a policy (actor) network
$\pi(a_{t}|o_{t})$ that predicts actions $a_{t}$ given observations $o_{t}$
and a value function network (critic) $V_{l}(o_{t})$ which maps observations
to estimates of the expected discounted future returns for each objective
indexed by $l$. The MO-VMPO algorithm was configured to optimize the reward
terms in Equations 3 to 7 separately using individual critic outputs for value
estimation. The policies were trained at a control rate of
$100\text{\,}\mathrm{H}\mathrm{z}$, with a discount factor of $0.99$, using a
batch size of $1024$, a trajectory unroll length of 20, MO-VMPO
$\epsilon_{k}=0.002$ for the $k$ objectives and a learning rate of $10^{-4}$.
Optimization of the networks was done with Adam [32]. See Appendix D for more
hyperparameter settings.

###  III-F Observations and Architecture

At each time step, the policy observed a tuple of features

$\displaystyle\phi(t)=[$
$\displaystyle\mathbf{M}(t),\mathbf{q}(t),\mathbf{\dot{q}}(t),\mathbf{v}(t),\mathbf{\omega}(t),\mathbf{p}_{i}(t),h(t),$
$\displaystyle\mathbf{c}(t),\mathbf{X}(t),\mathbf{\hat{X}}(t),a_{t-1},p^{y}_{\text{com}},\mathbf{z}(t)],$
(8)

where $\mathbf{M}$ is an image with local height information ($32\times 32$
pixels representing the
$$1.7\text{\,}\mathrm{m}$\times$1.7\text{\,}\mathrm{m}$$ square centered
$40\text{\,}\mathrm{c}\mathrm{m}$ in front of the CoM, updated at
$10\text{\,}\mathrm{H}\mathrm{z}$), $\mathbf{q}$ and $\mathbf{\dot{q}}$ are
the positions and velocities of the joints, $\mathbf{v}$ and $\mathbf{\omega}$
the linear and angular velocities of the base, $\mathbf{p}$ the end-effector
positions, $h$ the height of the base and $\mathbf{X}$ the orientation of the
base represented as a rotation matrix. The features $\mathbf{c}$ and
$\mathbf{\hat{X}}$ contain the relative position and orientation of the base
with respect to the pose of the robot at the last time step that the local
height image was updated. The previously taken action was observed as well.
The absolute position of the center of mass on the $y$-axis
$p^{y}_{\text{com}}$ was added to allow policies to track a straight path
during evaluation. Finally, a position ‘command’ $\mathbf{z}(t)$ was provided
which consists of the pair
$(p^{x,\text{ref}}_{\text{com}}(t+10)-p^{x}_{\text{com}}(t),p^{y,\text{ref}}_{\text{com}}(t+10))$,
i.e., the relative $x$-coordinate of the reference center of mass and its
absolute $y$-coordinate $10$ timesteps in the future. These commands were
added to the observation to help the controller disambiguate between different
trajectories to track. The choice to use absolute $y$-coordinates was again to
easily switch to straight line following in downstream tasks by setting this
$y$ command to zero.

The actor is a conditional normal distribution parameterized by a deep neural
network. Since the observations don’t provide full state information, we added
memory to the network by incorporating an LSTM [33]. Convolutional layers were
used to process the local terrain information. See Fig. 4 for a schematic
representation of the actor network architecture. The critic network received
the same observations as the actor but also an embedding vector representing
the clip number and a set of future reference body positions and quaternions.
See Appendix E for more details about the networks architectures used.
Finally, the critic parameters were kept during fine-tuning but it’s
additional observations (see Sec. III-F) were replaced with zeros. We didn’t
find that resetting the critic parameters lead to very different results (see
Appendix F).

Figure 4: Network architecture of the actor network up to the penultimate
layer. The final layer maps densely to the parameters of the Gaussian action
distribution.

###  III-G Terrain Distortions

For half the seeds of the tracking experiments, the terrains were _distorted_
to provide more diverse height map observations while respecting the locations
and orientations of the planned end-effector contacts. The $16\times 16$
height map of a clip was first embedded in the center of a $46\times 46$
matrix. Subsequently 8 sets of 4 integers representing the corner indices of
rectangles where sampled uniformly. The heights within sampled rectangles from
the middle region with the original height map or the flat region in front of
it were scaled with a factor sampled uniformly from $(0.7,1.0)$. If the
rectangle was in one of the other regions, the scaling was sampled from the
interval $(0.8,1.5)$. Finally, all terrain surface triangles overlapping with
$10\text{\,}\mathrm{c}\mathrm{m}$ squares centered on the end-effector contact
locations were reset to their original heights. See Fig. 5 for examples of
distorted terrains.

Figure 5: Examples of distorted terrains after embedding the height maps in a
larger area and rescaling the height values at random rectangular regions.

###  III-H Evaluation and Fine-Tuning

During zero-shot evaluation and fine-tuning experiments, the position commands
$\mathbf{z}(t)$ were fixed to $(0.05,0.0)$. During fine-tuning, the reward
function was also adapted. The squared difference in the reward term
$r_{\text{com}}$ was replaced with
$(v_{x}-0.5)^{2}+(p^{y,\text{ref}}_{\text{com}})^{2}$ to promote straight
walking at a fixed velocity and the other rewards terms were discarded. This
also means that no updates were done for the value heads of the critic
representing these rewards and that these value heads were ignored in the MO-
VMPO policy improvement step. Note that we replaced a position target here
with a velocity target but in a way that should produce a maximum rewards for
the same base velocity. The parameters of both the actor and the critic
networks were initialized with those of the tracking policy. Finally, the
additional privileged inputs of the critic were set to zero as they were
specific to the clip tracking task.

(a) Stairs

(b) Procedural

(c) Wavy Steps

(d) Mixed

Figure 6: Evaluation and fine-tuning environments used in the experiments.

##  IV Experiments

We used the Mujoco library [34] to simulate the dynamics of the full ANYmal
model. The actuators were modelled with a neural network trained on real data,
similar to what was done in other recent work on the ANYmal [7] to make the
simulation more realistic. Our qualitiative results can be seen in the
accompanying video.

###  IV-A Environments

To evaluate and fine-tune our policies, we used a variety of environments,
which are shown in Fig. 6. The _stairs_ environment consists of steps of
randomly varying heights and randomly positioned boxes that serve as
obstacles. The _procedural_ terrain environment contains height maps that were
created by the same process as was used to generate the planner trajectories.
The _wavy steps_ environment consists of a track of pairs of stepping stones
with randomized height differences and gaps between them. The pitch and roll
rotations of the stepping stones are also randomized and the overall elevation
pattern follows a sine wave. Finally, we used a _mixed_ environment which
combines the stairs, procedural terrains and wavy steps while adding two
additional types of terrains: flat regions with randomly spaced and sized
slits and another type of procedural terrain based on a higher resolution
height map and Perlin noise. See Appendix B for more detailed descriptions of
the environments and Appendix C for additional screenshots in a larger format.

###  IV-B Baselines

As baselines, we ran multiple seeds _from scratch_ on each environment using
the same architecture and learning algorithm as in the fine-tuning experiments
but with randomly initialized network weights. We ran three baseline seeds for
each environment, except for the wavy steps for which we observed more
variance and ran ten. We report the average performance over 100 episodes of
the best performing baseline seed. We also fine-tuned the baselines trained on
the procedural terrains on the wavy steps track and added these policies to
the seeds to take the maximum over for this environment. This was done to
disentangle curriculum learning benefits of first training on procedural
terrains from potential benefits due to properties of the imitated
trajectories. All baselines were trained for 20 billion environment steps but
performance typically stopped improving much earlier.

###  IV-C Tracking Performance and Zero-Shot Transfer

To find out how much data is needed to achieve good generalization, we trained
tracking policies using 8 different train set sizes. To monitor tracking
performance, we used the sum of all the reward terms defined in Equations 3 to
7. As performance kept increasing for the distorted terrains, we trained for
50 billion steps. As is clear from Fig. 7(a) and Fig. 7(b), held-out and train
performances keep approaching each other as the size of the dataset increases,
showing that sufficiently many trajectories are needed for good
generalization. The terrain distortions don’t seem to affect this pattern. As
shown in the accompanying video, the policies trained on the full dataset can
track the planned trajectories over challenging terrains that require precise
foot placements.

(a)

(b)

(c)

(d)

Figure 7: Convergence plots for the tracking task with distortion (a) and
without distortion (b) and fine-tuning convergence on the wavy steps (c) and
the mixed terrains (d).

We evaluated the tracking policies without fine-tuning on the _procedural_ and
_stairs_ environments. Episodes are $60\text{\,}\mathrm{s}$ long and we report
the 100 episode average distance traversed in Fig. 8(a) and 8(b). As Fig. 8(a)
shows, zero-shot performance on the procedural environment is better for
larger datasets. The distortions seem to help a bit for the larger datasets
but hurt performance for smaller sets. The policies trained without
distortions traverse about $18\text{\,}\mathrm{m}$ of the terrains on average.
This shows, together with the tracking performance on held-out data, that the
imitation learning setup can generalize to new terrains. While this zero-shot
performance is on par with the baseline (which was trained on this terrain for
20 billion steps), the behaviors of the tracking policies look qualitatively
very different. The performance on the stairs environment (not seen during
training) is also directly related to the size of the datasets, as can be seen
in Fig. 8(b) and the best performing policies are able to traverse more than
$6\text{\,}\mathrm{m}$ on average. The task specific baseline performed much
better on this task with an average of over $23\text{\,}\mathrm{m}$ but we
consider the zero-shot performance promising.

###  IV-D Fine-tuning

We fine-tuned the tracking policies on the wavy steps and mixed environments
using 15 billion steps. The former was constructed to be hard to learn with
low-level trial-and-error exploration while the latter was constructed to test
generality. Fig. 7(c), 7(d) and 8(e) show that convergence speed is generally
correlated with the size of the tracking dataset. As Fig. 8(c) shows, the wavy
steps were hard to learn and policies would either get stuck at about
$2.5\text{\,}\mathrm{m}$ like all our baseline seeds or learn the task and
traverse about $20\text{\,}\mathrm{m}$ on average. There seems to be a
relation between dataset size and final performance. Distortions seem to hurt
performance on this environment. As Fig. 8(d) shows, all the tracking policies
trained on all but the smallest dataset performed better than the baseline for
the mixed terrains. The effect of the distortions on final performance is not
clear. However, as shown in Fig. 8(e), the number of steps needed to fine-tune
past a reward of 500 (an arbitrary number that all policies achieved at some
point) was lower for the policies trained with distortions and larger
datasets. As the accompanying video shows, our policies can traverse both the
wavy stepping stones parts and the rough procedural parts of the terrain. A
limitation is that the gaits diverge more from the planned trajectories as
fine-tuning continues.

To evaluate how promising the behaviors are for transfer to the real platform,
we looked at the distributions of the joint torques and velocities. As shown
in Fig. 8(f), there are no significant peaks above the suggested maximum
torque of $40\text{\,}\mathrm{N}\mathrm{m}$ [35]. The largest observed torque
was $54.78\text{\,}\mathrm{N}\mathrm{m}$ and only $1.28\%$ of the values were
above $40\text{\,}\mathrm{N}\mathrm{m}$. The velocity distribution looks well-
behaved too. The largest observed value was
$21\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$ and only
$1.68\%$ of the observed values surpassed the suggested maximum of
$12\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$. We think that
these values do not provide reason to suspect that the behaviors would be
infeasible on the real robot.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: (a-d) Meters traversed for the procedural and stairs environments
(zero shot) and the wavy steps and mixed terrains (after fine tuning). (d)
Number of steps needed to obtain more than 500 reward on the mixed terrains
(rolling average with a window size corresponding to $10\%$ of the data). (f)
Torque and velocity histograms measured for a policy fine-tuned on the wavy
steps terrain.

##  V Conclusion and Discussion

We showed that imitation from a diverse set of planned trajectories can be
used to learn dynamic terrain-adaptive controllers. After fine-tuning, the
policies can traverse environments that are very hard to solve with a standard
RL setup. Our data augmentation method showed mixed results but aided fine-
tuning convergence speed. While we only presented results in simulation, we
used a realistic robot model and observations that can all at least be
estimated on a real platform. The fine-tuning stage doesn’t prevent the
policies from drifting away towards gaits that don’t resemble the tracking
policy anymore. This can be partially remedied by early stopping but a more
principled method would be to add a penalty term that enforces the policies to
remain similar to the tracking policy or to the planner trajectories. A
possible way to do this would be via adversarial motion priors [26]. In future
work, we also want to add domain randomization and test the controllers on a
real robot via sim2real. We hope that our method for learning dynamic terrain-
adaptive controllers will allow robot platforms to move more dynamically over
more challenging terrains with less engineering effort and we believe that it
provides a useful starting point for research on terrain-adaptive skill reuse
and adaptation.

## Acknowledgment

We’d like to thank Michael Neunert, Josh Merel, Fereshteh Sadeghi, Markus
Wulfmeier, Yuval Tassa, Michael Bloesch, Saran Tunyasuvunakool, Raia Hadsell
and Emilio Parisotto for their feedback and useful discussions. We also thank
Piotr Trochim for technical support.

## References

  * [1] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science robotics_ , vol. 5, no. 47, 2020.
  * [2] J. Siekmann, K. Green, J. Warila, A. Fern, and J. Hurst, “Blind bipedal stair traversal via sim-to-real reinforcement learning,” _arXiv preprint arXiv:2105.08328_ , 2021.
  * [3] M. Kalakrishnan, J. Buchli, P. Pastor, M. Mistry, and S. Schaal, “Fast, robust quadruped locomotion over challenging terrain,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 2665–2670.
  * [4] C. D. Bellicoso, F. Jenelten, C. Gehring, and M. Hutter, “Dynamic locomotion through online nonlinear motion optimization for quadrupedal robots,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2261–2268, 2018\.
  * [5] V. Tsounis, M. Alge, J. Lee, F. Farshidian, and M. Hutter, “Deepgait: Planning and control of quadrupedal gaits using deep reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3699–3706, 2020\.
  * [6] J. Carius, F. Farshidian, and M. Hutter, “Mpc-net: A first principles guided policy search,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2897–2904, 2020.
  * [7] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, 2019.
  * [8] K. Lowrey, A. Rajeswaran, S. Kakade, E. Todorov, and I. Mordatch, “Plan online, learn offline: Efficient learning and exploration via model-based control,” _arXiv preprint arXiv:1811.01848_ , 2018.
  * [9] S. Levine and V. Koltun, “Guided policy search,” in _Proceedings of the 30th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, S. Dasgupta and D. McAllester, Eds., vol. 28, no. 3. Atlanta, Georgia, USA: PMLR, 17–19 Jun 2013, pp. 1–9. [Online]. Available: https://proceedings.mlr.press/v28/levine13.html
  * [10] X. B. Peng, E. Coumans, T. Zhang, T.-W. Lee, J. Tan, and S. Levine, “Learning agile robotic locomotion skills by imitating animals,” _arXiv preprint arXiv:2004.00784_ , 2020.
  * [11] H. F. Song, A. Abdolmaleki, J. T. Springenberg, A. Clark, H. Soyer, J. W. Rae, S. Noury, A. Ahuja, S. Liu, D. Tirumala, _et al._ , “V-mpo: On-policy maximum a posteriori policy optimization for discrete and continuous control,” _arXiv preprint arXiv:1909.12238_ , 2019.
  * [12] N. Heess, D. TB, S. Sriram, J. Lemmon, J. Merel, G. Wayne, Y. Tassa, T. Erez, Z. Wang, S. Eslami, _et al._ , “Emergence of locomotion behaviours in rich environments,” _arXiv preprint arXiv:1707.02286_ , 2017.
  * [13] A. Escontrela, G. Yu, P. Xu, A. Iscen, and J. Tan, “Zero-shot terrain generalization for visual locomotion policies,” _arXiv preprint arXiv:2011.05513_ , 2020.
  * [14] H. Shi, B. Zhou, H. Zeng, F. Wang, Y. Dong, J. Li, K. Wang, H. Tian, and M. Q.-H. Meng, “Reinforcement learning with evolutionary trajectory generator: A general approach for quadrupedal locomotion,” _arXiv preprint arXiv:2109.06409_ , 2021.
  * [15] Z. Xie, H. Y. Ling, N. H. Kim, and M. van de Panne, “Allsteps: Curriculum-driven learning of stepping stone skills,” in _Computer Graphics Forum_ , vol. 39, no. 8. Wiley Online Library, 2020, pp. 213–224.
  * [16] N. Rudin, D. Hoeller, P. Reist, and M. Hutter, “Learning to walk in minutes using massively parallel deep reinforcement learning,” _arXiv preprint arXiv:2109.11978_ , 2021.
  * [17] S. Gangapurwala, M. Geisert, R. Orsolino, M. Fallon, and I. Havoutis, “Rloc: Terrain-aware legged locomotion using reinforcement learning and optimal control,” _arXiv preprint arXiv:2012.03094_ , 2020.
  * [18] Z. Xie, X. Da, B. Babich, A. Garg, and M. van de Panne, “Glide: Generalizable quadrupedal locomotion in diverse environments with a centroidal model,” _arXiv preprint arXiv:2104.09771_ , 2021.
  * [19] S. Gangapurwala, A. Mitchell, and I. Havoutis, “Guided constrained policy optimization for dynamic quadrupedal robot locomotion,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3642–3649, 2020.
  * [20] O. Melon, M. Geisert, D. Surovik, I. Havoutis, and M. Fallon, “Reliable trajectories for dynamic quadrupeds using analytical costs and learned initializations,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 1410–1416.
  * [21] T. Kwon, Y. Lee, and M. Van De Panne, “Fast and flexible multilegged locomotion using learned centroidal dynamics,” _ACM Transactions on Graphics (TOG)_ , vol. 39, no. 4, pp. 46–1, 2020.
  * [22] M. Bogdanovic, M. Khadiv, and L. Righetti, “Model-free reinforcement learning for robust locomotion using trajectory optimization for exploration,” _arXiv preprint arXiv:2107.06629_ , 2021.
  * [23] X. B. Peng, P. Abbeel, S. Levine, and M. van de Panne, “Deepmimic: Example-guided deep reinforcement learning of physics-based character skills,” _ACM Transactions on Graphics (TOG)_ , vol. 37, no. 4, pp. 1–14, 2018.
  * [24] J. Merel, L. Hasenclever, A. Galashov, A. Ahuja, V. Pham, G. Wayne, Y. W. Teh, and N. Heess, “Neural probabilistic motor primitives for humanoid control,” _arXiv preprint arXiv:1811.11711_ , 2018.
  * [25] L. Hasenclever, F. Pardo, R. Hadsell, N. Heess, and J. Merel, “Comic: Complementary task learning & mimicry for reusable skills,” in _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , ser. Proceedings of Machine Learning Research, vol. 119. PMLR, 2020, pp. 4105–4115. [Online]. Available: http://proceedings.mlr.press/v119/hasenclever20a.html
  * [26] X. B. Peng, Z. Ma, P. Abbeel, S. Levine, and A. Kanazawa, “Amp: Adversarial motion priors for stylized physics-based character control,” _arXiv preprint arXiv:2104.02180_ , 2021.
  * [27] D. Holden, T. Komura, and J. Saito, “Phase-functioned neural networks for character control,” _ACM Transactions on Graphics (TOG)_ , vol. 36, no. 4, pp. 1–13, 2017.
  * [28] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through Phase-Based End-Effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, July 2018.
  * [29] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006.
  * [30] M. A. Sen, V. Bakircioglu, and M. Kalyoncu, “Inverse kinematic analysis of a quadruped robot,” _International journal of scientific & technology research_, vol. 6, no. 9, pp. 285–289, 2017.
  * [31] A. Abdolmaleki, S. Huang, L. Hasenclever, M. Neunert, F. Song, M. Zambelli, M. Martins, N. Heess, R. Hadsell, and M. Riedmiller, “A distributional view on multi-objective policy optimization,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 11–22.
  * [32] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [33] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [34] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 5026–5033.
  * [35] M. Hutter, C. Gehring, D. Jud, A. Lauber, C. D. Bellicoso, V. Tsounis, J. Hwangbo, K. Bodie, P. Fankhauser, M. Bloesch, _et al._ , “Anymal-a highly mobile and dynamic quadrupedal robot,” in _2016 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2016, pp. 38–44.
  * [36] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” _arXiv preprint arXiv:1511.07289_ , 2015.
  * [37] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in _Icml_ , 2010.

##  Appendix A Terrain Generation Procedure

For each terrain, we first generated an $8\times 8$ matrix of coefficients
sampled uniformly from $[0,0.2]$ and took the Kronecker product with a
$2\times 2$ matrix of ones to create a pattern of tiles. Subsequently, a
random shift was applied to the columns to prevent the tiles from always being
at the same locations. Finally, three $x-y$ locations were sampled uniformly
from the region covered by the height map, which served as centers
$\mathbf{c}_{i}$ of radial basis functions
$f_{i}(\mathbf{p})=\exp(-\|\mathbf{c}_{i}-\mathbf{p}\|/h_{i})$, where
$\mathbf{p}$ is a point on the $xy$-plane and $h_{i}$ is a bandwidth parameter
which was chosen as $h_{i}=\exp(-2\eta_{i})$, where $\eta_{i}$ is uniformly
sampled from $[0,1]$. These radial basis functions were used to scale down all
the values in the terrain grid based on their distance of the centers. The
height fields were scaled down to always have a maximum height of
$30\text{\,}\mathrm{c}\mathrm{m}$ and represent an area of
$$2\text{\,}\mathrm{m}$\times$2\text{\,}\mathrm{m}$$.

##  Appendix B Environment Details

###  B-A Procedural

A series of height fields generated with the procedure in Sec. A are generated
(15 in total) and connected together by platforms of
$$3\text{\,}\mathrm{m}$\times$2\text{\,}\mathrm{m}$$.

###  B-B Stairs

The step heights were uniformly sampled between
$0\text{\,}\mathrm{c}\mathrm{m}$ and $10\text{\,}\mathrm{c}\mathrm{m}$. The
step lengths were sampled between $10\text{\,}\mathrm{c}\mathrm{m}$ and
$75\text{\,}\mathrm{c}\mathrm{m}$. The steps were $2\text{\,}\mathrm{m}$ wide.
Each step also had three boxes that with heights sampled from the same
distribution as the steps heights. The horizontal positions of the boxes were
sampled to be within some maximum distance of each other that was sampled from
$($0\text{\,}\mathrm{m}$,$1\text{\,}\mathrm{m}$0)$ every episode. This would
cause the boxes to be somewhat bundled together. Whenever two boxes overlapped
with each other, their positioning along the length of the stairs ($x$-axis)
was offset by $5\text{\,}\mathrm{c}\mathrm{m}$.

###  B-C Wavy Steps

The _wavy steps_ environment was constructed out of a series of pairs of boxes
for the robot to step on. To make the task more challenging, various
properties of the track were randomized. All distributions over continuous
intervals were uniform. All units are centimeters. The spaces separating each
pair of steps laterally were sampled from $(10,15)$, the vertical offset and
longitudinal offset of each step from $(-5,5)$, the length of the gaps in the
longitudinal direction from $(15,20)$ and the length and the width of the
steps of each pair from $(20,40)$ and $(30,60)$ respectively. The steps were
$72\text{\,}\mathrm{c}\mathrm{m}$ high and we applied random rotations on the
pitch and roll axes with values sampled from $(-0.15,0.15)$ radians. Note that
these rotations could significantly widen the gaps between the steps. Finally,
the function $\sin((x-6)/3)$ was added to the heights of the steps, where $x$
is the longitudinal distance in meters.

###  B-D Mixed

The mixed terrains contained five segments which where randomly sampled from
to create different terrain combinations at each episode. The procedural and
stairs terrains were identical to the descriptions above. A procedural segment
corresponded to a single terrain segment of $4\text{\,}\mathrm{m}$. A stairs
segment corresponded to five steps. The third terrain segment was
$4\text{\,}\mathrm{m}$ of the wavy steps track but without the sine function
being added to the heights. The forth terrain was a segment of
$3\text{\,}\mathrm{m}$ of flat platforms separated by slits/gaps. The platform
lengths were sampled from $($0.1\text{\,}\mathrm{m}$,$1\text{\,}\mathrm{m}$)$
and the lengths of the gaps from
$($0.15\text{\,}\mathrm{m}$,$0.2\text{\,}\mathrm{m}$)$. Finally the fifth
terrain was a $4\text{\,}\mathrm{m}$ long height map with Perlin noise (6
harmonics) of a maximum height of $0.5\text{\,}\mathrm{m}$ and a resolution of
65 vertices per meter. All segments were $2\text{\,}\mathrm{m}$ wide and the
total track combined 15 of them.

##  Appendix C Additional Images of the Environments

See Fig. 9.

(a) Stairs

(b) Procedural

(c) Wavy Steps

(d) Wavy Steps (side view)

(e) Mixed

(f) Mixed (different sample)

Figure 9: Additional and larger images of the environments used in the
experiments.

##  Appendix D Learning Algorithm Hyperparameters

See Tabel I for an overview of the hyperparameter settings used in the
experiments. We refer to [31] and [11] for the details of the MO-VMPO
algorithm but will explain the relevant hyperparameters here. The MO-VMPO
algorithm consists of an E-step and M-step which both optimize objectives
under KL constraints using Lagrange multipliers. The E-step constraints are
defined for each of the $k$ objectives via hyperparameters $\epsilon_{l}$ and
corresponding Lagrange multpliers $\eta_{l}$. The M-step only has a single KL-
divergence score but since this is a KL-divergence between Gaussian
distributions, the mean and covariance terms can be separated and are
constrained using their own individual hyperparameter and lagrange multiplier
pairs $(\epsilon_{\mu},\alpha_{\mu})$ and
$(\epsilon_{\Sigma},\alpha_{\Sigma})$, respectively. The Lagrange multipliers
are updated automatically but their initial values can strongly influence
learning and are therefore relevant hyperparameters as well.

TABLE I: Learning hyperparameters. Hyperparameter | Value
---|---
training |
Adam learning rate | $10^{-4}$
target network update period | $100$
discount | $0.99$
batch size | $1024$
MO-VMPO |
variational distribution KL constraints, $\epsilon_{l}$ | $0.002$
initial variational dist. Lagrange multipliers, $\eta_{l,0}$ | $10$
mean KL constraint, $\epsilon_{\mu}$ | $0.1$
covariance KL constraint, $\epsilon_{\Sigma}$ | $10^{-5}$
initial mean lagrange multiplier, $\alpha_{\mu,0}$ | $0.1$
initial covariance Lagrange multiplier, $\alpha_{\Sigma,0}$ | $100$

##  Appendix E Architecture Details

The actor is a conditional normal distribution parameterized by a deep neural
network. Since the observations don’t provide full state information, we added
memory to the network by incorporating an LSTM [33]. More precisely, the local
terrain information $\mathbf{M}(t)$ was first processed by two convolutional
layers without pooling with respectively $16$ filters of dimension $5\times 5$
and $32$ filters of dimension $3\times 3$. Each layer was followed by an ELU
non-linearity [36]. The final outputs were passed through a fully connected
layer of 128 dimensions and concatenated with the remaining observation
features. Subsequently, the concatenated features were passed through a 256
unit fully connected tanh layer, followed by the 512 unit large LSTM. The LSTM
output was passed through by one more 512 unit fully connected tanh layer and
finally projected to the parameters of a multivariate normal distribution with
diagonal covariance.

The critic received the same observations as the actor but also an embedding
vector representing the clip number and a set of future reference body
positions and quaternions $\\{\xi_{i}(t+k)\colon
i\in\mathcal{B},k\in\\{2,4,6,8,10\\}\\}$, where $\mathcal{B}$ is the set of
body indices. The terrain and proprioceptive features went through the same
computations as in the actor (but with separate parameters) up to the
concatenation of the image and proprioceptive paths. This was followed by
three fully connected layers (two tanh and one ReLU [37], an unintentional
inconsistency) and finally a linear transformation to each of the value
predictions. The critic parameters were kept during fine-tuning but it’s
additional observations (see Sec. III-F) were replaced with zeros.

##  Appendix F Comparison of Preserving or Resetting the Critic Parameters

To see if resetting the critic instead of preserving it would lead to
substantially different results, we ran an additional fine-tuning experiment
on the wavy steps task. As can be seen in Fig. 10, the overall pattern of
performance looks the same for both preserving or resetting the critic
weights. The most notable difference is that without distortions, the run in
which the critic was reset failed on the largest dataset in the sweep. We
suspect that this is just an unlucky seed and that more runs are needed to
verify this result. However, overall, these results don’t suggest that
rerunning the full experiments while resetting the critic would have lead to
substantially different results.

(a)

(b)

Figure 10: Comparisons of the average rewards over the first 8 billion steps
for fine-tuning on wavy steps when preserving both the actor and critic
weights or only those of the actor.


# Poisoning $\times$ Evasion: Symbiotic Adversarial Robustness for Graph
Neural Networks

Ege Erdogan, Simon Geisler, Stephan Günnemann
{ege.erdogan, s.geisler<EMAIL_ADDRESS>
Department of Computer Science
Technical University of Munich

###### Abstract

It is well-known that deep learning models are vulnerable w.r.t. small input
perturbations. Such perturbed instances are called adversarial examples.
Adversarial examples are commonly crafted to fool a model either at training
time (poisoning) or test time (evasion). In this work, we study the symbiosis
of poisoning and evasion. We show that combining both threat models can
substantially improve the devastating efficacy of adversarial attacks.
Specifically, we study the robustness of Graph Neural Networks (GNNs) under
structure perturbations and devise a memory-efficient adaptive end-to-end
attack for the novel threat model using first-order optimization.

##  1 Introduction

Figure 1: Perturbed accuracies (with standard errors) after evasion,
poisoning, and symbiotic (ours) attacks on different models on PubMed.

Graph neural networks (GNNs) are increasingly used in different domains such
as product recommendations (Hao et al.,, 2020) and drug discovery (Guo et
al.,, 2021). Nevertheless, GNNs are vulnerable to adversarial attacks across
many tasks such as node classification (Zügner et al.,, 2018; Dai et al.,,
2018; Geisler et al.,, 2021), graph classification (Dai et al.,, 2018; Wang et
al.,, 2023), link prediction (Chen et al.,, 2020), and node embeddings
(Bojchevski and Günnemann,, 2019; Zhang et al.,, 2019). With attacks being
able to scale to very large graphs (Geisler et al.,, 2021), studying the
adversarial robustness of GNNs is of growing importance. GNNs can be attacked
during test time (evasion) or train time (poisoning); yet, a threat model
combining evasion and poisoning has not been considered in the literature. It
is nevertheless a reasonable threat model considering e.g. publicly available
graphs, or graphs extracted from sources such as social media sites.

The Problem. We consider node classification tasks, and an adversary able to
change the structure of the graph (i.e. insert/remove edges), with both train
and test time access to the graph. The adversary’s end goal is to minimize the
classification accuracy on the test set. It is constrained by a global budget
and attacks entire graph at once rather than targeting specific nodes. We call
such attacks combining evasion and poisoning in this threat model symbiotic
attacks.

Contributions. We initiate the study of this threat model and compare it with
plain poisoning and evasion adversaries by adapting the previous PR-BCD attack
(Geisler et al.,, 2021) to the symbiotic threat model, resulting in memory-
efficient attacks that can scale to large graphs. Our main findings are:

  * •

The symbiotic attacks are consistently stronger than a poisoning attack,
indicating that it is beneficial to allocate part of the resources of the
poisoning attack for the evasion objective.

  * •

Evasion attacks are constrained by the number of test nodes, with larger test
sets making an evasion attack harder. The symbiotic attacks are affected less
significantly by the size of the test set since they can also manipulate the
graph during training, leading to almost-zero accuracy in cases such as when
the share of labeled train is low and test nodes high.

Overall, the significance of the potential improvement the symbiotic threat
model provides indicates that it requires further study, and we outline paths
of future work in our conclusion.

##  2 Preliminaries

Notation. Throughout we denote a graph by $\mathcal{G}$ with $n$ nodes,
adjacency matrix $\mathbf{A}\in\\{0,1\\}^{n\times n}$ and feature matrix
$\mathbf{X}\in\mathrm{R}^{n\times d}$, and a GNN applied to the graph as
$f_{\theta}(\mathcal{G})$ with parameters $\theta$. $\Phi(\mathcal{G})$ is the
set of admissible adversarial graphs resulting from $\mathcal{G}$, and
$\mathcal{L}_{\text{atk}}$ and $\mathcal{L}_{\text{train}}$ denote the
adversarial and training objectives.

###  2.1 Adversarial Robustness of GNNs

An adversarial attack on a GNN can change the structure of the graph by
inserting/removing edges and nodes, or can modify the node features. We focus
on node classification and edge-level structural perturbations.

Attacks can be divided into two categories: evasion and poisoning. An evasion
attack targets a fixed GNN with $\theta$ obtained on a clean graph, and thus
tries to solve the optimization problem

$\max_{\hat{\mathcal{G}}\in\Phi(\mathcal{G})}\mathcal{L}_{\text{atk}}(f_{\theta}(\hat{\mathcal{G}})),$
(1)

while a poisoning attack is performed before training but aims to degrade
performance after training:

$\max_{\hat{\mathcal{G}}\in\Phi(\mathcal{G})}\mathcal{L}_{\text{atk}}(f_{\theta^{*}}(\hat{\mathcal{G}}))\quad\text{where}\quad\theta^{*}=\operatorname*{argmin}_{\theta}\mathcal{L}_{\text{train}}(f_{\theta}(\hat{\mathcal{G}})).$
(2)

A poisoning attack is admittedly more challenging. Previous work has tried
using evasion perturbations as poisoning perturbations (Zügner et al.,, 2018),
or unrolling the training procedure as part of the optimization to compute
meta-gradients (gradients w.r.t. hyperparameters) of
$\mathcal{L}_{\text{atk}}$ w.r.t. $\mathbf{A}$ (Zügner et al.,, 2020).

Finally, since we only consider changes to the binary adjacency matrix, we
define $\Phi(\mathcal{G})$ to include graphs reachable from $\mathcal{G}$
after at most $\Delta$ edge perturbations (e.g. for undirected $\mathcal{G}$,
$\Phi(\mathcal{G})=\\{\tilde{\mathcal{G}}\,|\,\|\tilde{\mathbf{A}}-\mathbf{A}\|_{0}\leq
2\Delta\,\land\,\tilde{\mathbf{A}}^{\top}=\tilde{\mathbf{A}}\\}$), though it
would also be straightforward to use metrics such as the graph’s degree
distribution or diameter.

####  2.1.1 PR-BCD

Our work builds on the Projected Randomized Block Coordinate Descent (PR-BCD)
attack proposed in Geisler et al., (2021). Similar to the Projected Gradient
Descent (PGD) attack (Xu et al.,, 2019), the adjacency matrix is relaxed to
$\mathbf{P}\in[0,1]^{n\times n}$ to enable continuous gradient updates, and
each entry denotes the probability of flipping that edge with the final
perturbations sampled from $\text{Bernoulli}(\mathbf{P})$. However, since the
adjacency matrix grows quadratically with the number of nodes, scalability of
plain PGD becomes a challenge on larger graphs.

PR-BCD employs Randomized Block Coordinate Descent (R-BCD) (Nesterov,, 2012)
and updates a block of size $b$ of $\mathbf{P}$ in each iteration. The
projection step ensures that the budget is enforced in expectation; i.e.
$\mathrm{E}[\text{Bernoulli}(\mathbf{P})]=\sum\mathbf{P}\leq\Delta$ and that
$\mathbf{P}\in[0,1]^{n\times n}$. After each iteration, rather than sampling
the entire block again, the promising entries of the block are kept and only
the rest is resampled. In other words, $\mathbf{P}$ is kept sparse in a
survival-of-the-fittest manner. Multiple samples are drawn at the end and the
best-performing one is returned as the final perturbed graph.

PGD can also be applied for a poisoning attack as the Meta-PGD attack in
Mujkanovic et al., (2022). In our attacks, we employ the same principle with
PR-BCD to scale better to larger graphs.

While we only consider a single global budget $\Delta$, it is straightforward
to include more sophisticated constraints when beneficial for the application
at hand (Gosch et al., 2023b, ; Gosch et al., 2023a, ).

##  3 Symbiotic Attacks

The Symbiotic Objective. The problem of a symbiotic attack has a similar form
to the bi-level optimization of Equation 2, but the main objective is
conditioned on the evasion graph $\mathcal{G}^{*}$ in addition to the
parameters $\theta^{*}$:

$\displaystyle\begin{split}\max_{\hat{\mathcal{G}}\in\Phi(\mathcal{G})}\mathcal{L}_{\text{pois}}(f_{\theta^{*}}(\mathcal{G}^{*}))\quad&\text{where}\quad\theta^{*}=\operatorname*{argmin}_{\theta}\mathcal{L}_{\text{train}}(f_{\theta}(\hat{\mathcal{G}})),\\\
&\text{and}\quad\mathcal{G}^{*}=\operatorname*{argmax}_{\tilde{\mathcal{G}}\in\Phi(\hat{\mathcal{G}})}\mathcal{L}_{\text{ev}}(f_{\theta^{*}}(\tilde{\mathcal{G}}))\end{split}$
(3)

where for clarity we separate the poisoning and evasion objectives
$\mathcal{L}_{\text{pois}}$ and $\mathcal{L}_{\text{ev}}$ although they might
be the same function.

Threat Model. We model an attacker with the goal of degrading a model’s
performance on node classification tasks. Combining poisoning and evasion
adversaries, our attacker has full access to the graph both at train and test
times, knows the model architecture being used so that it can create surrogate
models, but can access the trained model only as a black-box. Finally, our
attacker is limited by a global budget of edge insertions/removals.

With the symbiotic objective and the threat model, we propose two attacks as
approximations to the optimal solution.

The Sequential Attack. A simple way of launching a symbiotic attack is to
split the budget into two, and then to launch an evasion attack with the
second part after poisoning the model using the first part. In this case, the
poisoning attack has no knowledge of a future evasion, but it can can still
help the evasion attack by causing some nodes to be misclassified and reducing
the classification margin of some nodes, making them easier to be
misclassified during an evasion attack.

The Joint Attack. The poisoning attack can also be designed to “fit" the
future evasion graph by including the evasion attack within the poisoning
loss. In this joint attack, an evasion attack is launched within each
poisoning attack iteration, and the poisoning loss is computed with the
poisoned model over the evasion graph, so that the poisoning attack not only
reduces the model’s accuracy but also makes it more vulnerable to evasion.

Although the joint and sequential attacks can be instantiated using different
evasion/poisoning attacks, we choose to build on PR-BCD since it scales well
to larger graphs. This is a significant consideration since especially for the
joint attack the inner evasion attack has to performed many times. Also note
that the sequential attack is in fact a special case of the joint attack with
zero iterations per inner evasion attack.

##  4 Evaluation

Table 1: Numbers of nodes, edges, and classes in the datasets we include in
our evaluations.
Dataset | Nodes | Edges | Classes
---|---|---|---
Cora | 2,708 | 10,556 | 7
CiteSeer | 3,327 | 9,104 | 6
PubMed | 19,717 | 88,648 | 3

###  4.1 Setup

In this section, we compare the symbiotic threat model with evasion and
poisoning attacks, instantiated through the PR-BCD evasion (Geisler et al.,,
2021) and poisoning attacks (Mujkanovic et al.,, 2022) on Cora, CiteSeer
(McCallum et al.,, 2000), and PubMed (Sen et al.,, 2008). We study the
robustness of GCN (Kipf and Welling,, 2016), GAT (Veličković et al.,, 2018),
APPNP (Gasteiger et al.,, 2018), and GPRGNN (Chien et al.,, 2020). We also
consider R-GCN (Zhu et al.,, 2019) and Jaccard purification (Wu et al.,, 2019)
as potential defensive mechanisms. For each dataset, we allocate 20 nodes of
each class for the labeled training set and allocate 10% of the nodes as the
test set. For transductive tasks, the test nodes are also included as
unlabeled train nodes during training, whereas for inductive tasks the test
nodes are added to the graph after training. Table 1 displays the numbers of
nodes, edges, and classes in each of our benchmark datasets.

Figure 2: Perturbed accuracies and standard errors after the four attacks with
different test set sizes and a fixed 5% global budget, using a GCN on the
three benchmark datasets.

###  4.2 Results

Table 2 displays the perturbed accuracy values on the test set (10% of nodes)
for our benchmark datasets/models averaged over 10 runs, along with the
standard error of the mean. We limit the attacker with a budget 5% the number
of edges, with the budget split equally between poisoning and evasion for the
symbiotic attacks. To better focus on the threat model, we report the better-
performing of the two attacks for the symbiotic threat model and leave a
comparison between the attacks to the appendix.

Symbiotic attacks are stronger than poisoning across all tasks, and stronger
than plain evasion especially against GCN and GPRGNN models. The effect of the
symbiotic threat model is most evident on the larger PubMed graph, with the
accuracy dropping to almost zero e.g. against a GCN, especially as the evasion
attack’s performance drops against larger test sets as we discuss next.

Table 2: Average (± standard error) perturbed accuracies for the evasion, poisoning, and symbiotic attacks with a 5% budget. The -J suffix indicates the graph has been pre-processed with Jaccard purification (Wu et al.,, 2019) and (ind.) stands for inductive learning. The strongest (lowest accuracy) results for each setup are written in bold. Model | Dataset | Clean | Evasion | Poisoning | Symbiotic
---|---|---|---|---|---
GCN | CiteSeer | 0.68 ± 0.01 | 0.41 ± 0.01 | 0.4 ± 0.01 | 0.38 ± 0.01
CiteSeer (ind.) | 0.67 ± 0.01 | 0.41 ± 0.01 | 0.62 ± 0.01 | 0.33 ± 0.01
CiteSeer-J | 0.68 ± 0.01 | 0.41 ± 0.01 | 0.41 ± 0.02 | 0.38 ± 0.01
| Cora | 0.78 ± 0.01 | 0.41 ± 0.01 | 0.46 ± 0.02 | 0.35 ± 0.01
| Cora (ind.) | 0.75 ± 0.02 | 0.42 ± 0.01 | 0.68 ± 0.03 | 0.3 ± 0.01
| Cora-J | 0.74 ± 0.01 | 0.39 ± 0.01 | 0.43 ± 0.02 | 0.36 ± 0.01
| PubMed | 0.78 ± 0.01 | 0.41 ± 0.01 | 0.12 ± 0.02 | 0.03 ± 0.01
| PubMed-J | 0.77 ± 0.01 | 0.41 ± 0.01 | 0.11 ± 0.01 | 0.02 ± 0.0
GAT | CiteSeer | 0.62 ± 0.02 | 0.27 ± 0.02 | 0.41 ± 0.02 | 0.3 ± 0.03
CiteSeer (ind.) | 0.68 ± 0.01 | 0.37 ± 0.01 | 0.64 ± 0.02 | 0.56 ± 0.02
CiteSeer-J | 0.64 ± 0.01 | 0.32 ± 0.03 | 0.41 ± 0.03 | 0.3 ± 0.03
| Cora | 0.69 ± 0.02 | 0.22 ± 0.02 | 0.48 ± 0.03 | 0.29 ± 0.02
| Cora (ind.) | 0.77 ± 0.01 | 0.21 ± 0.01 | 0.61 ± 0.04 | 0.35 ± 0.03
| Cora-J | 0.67 ± 0.01 | 0.23 ± 0.02 | 0.45 ± 0.02 | 0.28 ± 0.02
| PubMed | 0.73 ± 0.01 | 0.38 ± 0.04 | 0.41 ± 0.01 | 0.2 ± 0.03
| PubMed-J | 0.74 ± 0.01 | 0.34 ± 0.04 | 0.38 ± 0.04 | 0.19 ± 0.02
APPNP | CiteSeer | 0.69 ± 0.01 | 0.45 ± 0.01 | 0.56 ± 0.01 | 0.47 ± 0.01
CiteSeer (ind.) | 0.71 ± 0.01 | 0.47 ± 0.01 | 0.66 ± 0.02 | 0.4 ± 0.01
CiteSeer-J | 0.68 ± 0.01 | 0.43 ± 0.01 | 0.52 ± 0.02 | 0.45 ± 0.02
| Cora | 0.82 ± 0.02 | 0.48 ± 0.03 | 0.64 ± 0.02 | 0.51 ± 0.04
| Cora (ind.) | 0.82 ± 0.02 | 0.53 ± 0.02 | 0.78 ± 0.01 | 0.37 ± 0.01
| Cora-J | 0.82 ± 0.01 | 0.5 ± 0.01 | 0.67 ± 0.01 | 0.54 ± 0.01
| PubMed | 0.79 ± 0.0 | 0.46 ± 0.01 | 0.21 ± 0.02 | 0.09 ± 0.01
| PubMed-J | 0.77 ± 0.01 | 0.45 ± 0.01 | 0.19 ± 0.03 | 0.1 ± 0.02
GPRGNN | CiteSeer | 0.66 ± 0.01 | 0.34 ± 0.01 | 0.44 ± 0.02 | 0.33 ± 0.01
CiteSeer (ind.) | 0.67 ± 0.01 | 0.37 ± 0.01 | 0.56 ± 0.01 | 0.34 ± 0.01
CiteSeer-J | 0.65 ± 0.01 | 0.35 ± 0.01 | 0.44 ± 0.01 | 0.35 ± 0.01
| Cora | 0.82 ± 0.01 | 0.46 ± 0.01 | 0.53 ± 0.01 | 0.4 ± 0.01
| Cora (ind.) | 0.8 ± 0.02 | 0.44 ± 0.01 | 0.74 ± 0.01 | 0.35 ± 0.01
| Cora-J | 0.79 ± 0.01 | 0.44 ± 0.01 | 0.54 ± 0.01 | 0.4 ± 0.01
| PubMed | 0.78 ± 0.01 | 0.42 ± 0.01 | 0.28 ± 0.03 | 0.08 ± 0.02
| PubMed-J | 0.78 ± 0.01 | 0.42 ± 0.01 | 0.38 ± 0.04 | 0.15 ± 0.04
RGCN | CiteSeer | 0.63 ± 0.01 | 0.39 ± 0.01 | 0.59 ± 0.02 | 0.47 ± 0.01
Cora | 0.74 ± 0.02 | 0.44 ± 0.01 | 0.74 ± 0.01 | 0.52 ± 0.02
PubMed | 0.77 ± 0.01 | 0.43 ± 0.01 | 0.42 ± 0.04 | 0.15 ± 0.03
Figure 3: Classification accuracy after the four attacks on a GCN, with the
varying block sizes for the PR-BCD optimization displayed along the x-axis.

###  4.3 Effect of the Number of Test Nodes

To illustrate a fundamental point of difference between poisoning and evasion
objectives, Figure 2 displays the perturbed accuracies for evasion, poisoning,
and symbiotic attacks over varying fractions of test nodes using a GCN with a
5% global budget. Attacking a very small number of nodes is easy as all
attacks can obtain zero accuracy since a few edge perturbations should be
enough to manipulate a small number of nodes.

As the number of test nodes grows, evasion becomes considerably more difficult
across all datasets. Although poisoning and symbiotic attacks become more
difficult as well with more test nodes, especially on PubMed they are much
more robust than the evasion attack. Thus the degrading performance cannot be
explained only by the attacks having to use the same budget to target a larger
number of nodes. The poisoning attack is more lightly affected as it can also
manipulate the flow of information across the graph during training; e.g. the
PubMed in our setup has the lowest share of labeled train nodes (20 nodes per
class), which leads to poisoning attacks having devastating impact.

Similarly for the symbiotic attacks, poisoning helps in two ways: by reducing
the base accuracy evasion starts with, and possibly changing the structure of
the graph in a more effective way to block the flow of label information from
the labeled train nodes. This results in symbiotic attacks being both more
robust against larger test sets than plain evasion, while also being stronger
than poisoning alone.

Figure 4: Perturbed accuracy of a GCN after the four attacks. The x-axis shows
the global budget as a fraction of the number of edges.

###  4.4 Hyperparameters

Block size. Figure 3 shows the results of the four attacks with varying block
sizes, a fixed 5% budget, and 125 attack iterations against a GCN. We observe
that for very small block sizes, the attacks are less effective since the PR-
BCD optimization can cover only a small part of the adjacency matrix. However,
the marginal benefit of larger block sizes decreases once a large part of the
adjacency matrix can be covered.

Budget. Considering different fractions of the number of edges as the global
budget, Figure 4 indicates that all four attacks follow a similar trend with
an increasing budget. Especially on PubMed where the share of labeled train
nodes is the lowest, changing 5% of the edges is enough to achieve close to
zero accuracy under the symbiotic threat model. This highlights the
potentially devastating effect of the joint attacks especially in large graphs
with few labeled train nodes.

##  5 Conclusion / Future Work

In this work, we considered the symbiotic threat model of a combined test and
train time adversary aiming to degrade the overall performance of a GNN on
node classification tasks. We proposed two methods of obtaining perturbations,
one a special case of the other, and demonstrated the potentially devastating
impact such an attacker can have compared to a plain poisoning or evasion
attacker. We now conclude by outlining potential lines of future work
regarding the symbiotic threat model.

Different Attacks as the Inner Evasion Step. The joint attack we described is
not limited to PR-BCD, and can be used with other evasion attacks, or attacks
designed specifically for the symbiotic setting, to obtain perturbations
within the poisoning attack.

Inductive Tasks, Local Budgets, Targeted Attacks. The combined threat model is
also applicable to inductive tasks in which the test graph is different from
the train graph. The poisoning half of the joint attack can be performed by
targeting the validation nodes rather than the test nodes, and the evasion
half is performed on the new test graph. Similarly, the attacks can also be
applied with small changes under per-node local budgets, or target specific
nodes. We plan to include further evaluations on these settings as our next
step.

New Evasion-Aware Poisoning Attacks. Finally, further novel poisoning attacks
can also be developed which utilize the knowledge of a future evasion in
different ways.

## Acknowledgments and Disclosure of Funding

This research was supported by the Helmholtz Association under the joint
research school “Munich School for Data Science - MUDS“.

## References

  * Bojchevski and Günnemann, (2019) Bojchevski, A. and Günnemann, S. (2019).  Adversarial Attacks on Node Embeddings via Graph Poisoning.
  * Chen et al., (2020) Chen, J., Lin, X., Shi, Z., and Liu, Y. (2020).  Link prediction adversarial attack via iterative gradient attack.  IEEE Transactions on Computational Social Systems, 7(4):1081–1094.
  * Chien et al., (2020) Chien, E., Peng, J., Li, P., and Milenkovic, O. (2020).  Adaptive Universal Generalized PageRank Graph Neural Network.  In International Conference on Learning Representations.
  * Dai et al., (2018) Dai, H., Li, H., Tian, T., Huang, X., Wang, L., Zhu, J., and Song, L. (2018).  Adversarial Attack on Graph Structured Data.  In Proceedings of the 35th International Conference on Machine Learning, pages 1115–1124. PMLR.
  * Gasteiger et al., (2018) Gasteiger, J., Bojchevski, A., and Günnemann, S. (2018).  Predict then Propagate: Graph Neural Networks meet Personalized PageRank.  In International Conference on Learning Representations.
  * Geisler et al., (2021) Geisler, S., Schmidt, T., Şirin, H., Zügner, D., Bojchevski, A., and Günnemann, S. (2021).  Robustness of Graph Neural Networks at Scale.  In Advances in Neural Information Processing Systems, volume 34, pages 7637–7649. Curran Associates, Inc.
  * (7) Gosch, L., Geisler, S., Sturm, D., Charpentier, B., Zügner, D., and Günnemann, S. (2023a).  Adversarial training for graph neural networks.  In Neural Information Processing Systems, NeurIPS.
  * (8) Gosch, L., Sturm, D., Geisler, S., and Günnemann, S. (2023b).  Revisiting robustness in graph machine learning.  In 11th International Conference on Learning Representations, ICLR.
  * Guo et al., (2021) Guo, J., Janet, J. P., Bauer, M. R., Nittinger, E., Giblin, K. A., Papadopoulos, K., Voronov, A., Patronov, A., Engkvist, O., and Margreitter, C. (2021).  Dockstream: a docking wrapper to enhance de novo molecular design.  Journal of cheminformatics, 13(1):1–21.
  * Hao et al., (2020) Hao, J., Zhao, T., Li, J., Dong, X. L., Faloutsos, C., Sun, Y., and Wang, W. (2020).  P-companion: A principled framework for diversified complementary product recommendation.  In Proceedings of the 29th ACM International Conference on Information & Knowledge Management, pages 2517–2524.
  * Kipf and Welling, (2016) Kipf, T. N. and Welling, M. (2016).  Semi-Supervised Classification with Graph Convolutional Networks.  In International Conference on Learning Representations.
  * McCallum et al., (2000) McCallum, A. K., Nigam, K., Rennie, J., and Seymore, K. (2000).  Automating the construction of internet portals with machine learning.  Information Retrieval, 3:127–163.
  * Mujkanovic et al., (2022) Mujkanovic, F., Geisler, S., Günnemann, S., and Bojchevski, A. (2022).  Are Defenses for Graph Neural Networks Robust?  Advances in Neural Information Processing Systems, 35:8954–8968.
  * Nesterov, (2012) Nesterov, Yu. (2012).  Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems.  SIAM Journal on Optimization, 22(2):341–362.
  * Sen et al., (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Galligher, B., and Eliassi-Rad, T. (2008).  Collective classification in network data.  AI magazine, 29(3):93–93.
  * Veličković et al., (2018) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y. (2018).  Graph Attention Networks.
  * Wang et al., (2023) Wang, X., Chang, H., Xie, B., Bian, T., Zhou, S., Wang, D., Zhang, Z., and Zhu, W. (2023).  Revisiting adversarial attacks on graph neural networks for graph classification.  IEEE Transactions on Knowledge and Data Engineering, pages 1–12.
  * Wu et al., (2019) Wu, H., Wang, C., Tyshetskiy, Y., Docherty, A., Lu, K., and Zhu, L. (2019).  Adversarial Examples on Graph Data: Deep Insights into Attack and Defense.
  * Xu et al., (2019) Xu, K., Chen, H., Liu, S., Chen, P.-Y., Weng, T.-W., Hong, M., and Lin, X. (2019).  Topology Attack and Defense for Graph Neural Networks: An Optimization Perspective.
  * Zhang et al., (2019) Zhang, H., Zheng, T., Gao, J., Miao, C., Su, L., Li, Y., and Ren, K. (2019).  Data poisoning attack against knowledge graph embedding.  In Proceedings of the 28th International Joint Conference on Artificial Intelligence, IJCAI’19, pages 4853–4859, Macao, China. AAAI Press.
  * Zhu et al., (2019) Zhu, D., Zhang, Z., Cui, P., and Zhu, W. (2019).  Robust Graph Convolutional Networks Against Adversarial Attacks.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1399–1407, Anchorage AK USA. ACM.
  * Zügner et al., (2018) Zügner, D., Akbarnejad, A., and Günnemann, S. (2018).  Adversarial Attacks on Neural Networks for Graph Data.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD ’18, pages 2847–2856, New York, NY, USA. Association for Computing Machinery.
  * Zügner et al., (2020) Zügner, D., Borchert, O., Akbarnejad, A., and Günnemann, S. (2020).  Adversarial Attacks on Graph Neural Networks: Perturbations and their Patterns.  ACM Transactions on Knowledge Discovery from Data, 14(5):57:1–57:31.

Table 3: Perturbed accuracies (± standard error) of the joint and sequential attacks under the symbiotic threat model with a 5% global budget. The -J suffix indicates the graph has been pre-processed with Jaccard purification (Wu et al.,, 2019). Model | Dataset | Clean | Sequential | Joint
---|---|---|---|---
GCN | CiteSeer | 0.68 ± 0.01 | 0.41 ± 0.01 | 0.38 ± 0.01
CiteSeer-J | 0.68 ± 0.01 | 0.4 ± 0.01 | 0.38 ± 0.01
Cora | 0.78 ± 0.01 | 0.37 ± 0.02 | 0.35 ± 0.01
| Cora-J | 0.74 ± 0.01 | 0.36 ± 0.01 | 0.36 ± 0.02
| PubMed | 0.78 ± 0.01 | 0.05 ± 0.01 | 0.03 ± 0.01
| PubMed-J | 0.77 ± 0.01 | 0.04 ± 0.01 | 0.02 ± 0.0
GAT | CiteSeer | 0.62 ± 0.02 | 0.3 ± 0.03 | 0.38 ± 0.02
CiteSeer-J | 0.64 ± 0.01 | 0.3 ± 0.03 | 0.36 ± 0.02
Cora | 0.69 ± 0.02 | 0.29 ± 0.02 | 0.32 ± 0.02
| Cora-J | 0.67 ± 0.01 | 0.28 ± 0.02 | 0.3 ± 0.03
| PubMed | 0.73 ± 0.01 | 0.24 ± 0.02 | 0.2 ± 0.03
| PubMed-J | 0.74 ± 0.01 | 0.27 ± 0.04 | 0.19 ± 0.02
APPNP | CiteSeer | 0.69 ± 0.01 | 0.47 ± 0.01 | 0.48 ± 0.01
CiteSeer-J | 0.68 ± 0.01 | 0.45 ± 0.02 | 0.45 ± 0.02
Cora | 0.82 ± 0.02 | 0.54 ± 0.02 | 0.51 ± 0.04
| Cora-J | 0.82 ± 0.01 | 0.57 ± 0.01 | 0.54 ± 0.01
| PubMed | 0.79 ± 0.0 | 0.09 ± 0.02 | 0.09 ± 0.01
| PubMed-J | 0.77 ± 0.01 | 0.1 ± 0.02 | 0.12 ± 0.02
GPRGNN | CiteSeer | 0.66 ± 0.01 | 0.34 ± 0.01 | 0.33 ± 0.01
CiteSeer-J | 0.65 ± 0.01 | 0.35 ± 0.01 | 0.35 ± 0.01
Cora | 0.82 ± 0.01 | 0.41 ± 0.01 | 0.4 ± 0.01
| Cora-J | 0.79 ± 0.01 | 0.42 ± 0.01 | 0.4 ± 0.01
| PubMed | 0.78 ± 0.01 | 0.08 ± 0.02 | 0.11 ± 0.03
| PubMed-J | 0.78 ± 0.01 | 0.16 ± 0.05 | 0.15 ± 0.04
RGCN | CiteSeer | 0.63 ± 0.01 | 0.47 ± 0.01 | 0.47 ± 0.01
Cora | 0.74 ± 0.02 | 0.56 ± 0.01 | 0.52 ± 0.02
PubMed | 0.77 ± 0.01 | 0.28 ± 0.04 | 0.15 ± 0.03
Figure 5: Perturbed accuracies and standard error after the joint attack on a
GCN, with the x-axis indicating the number of iterations for the inner evasion
part of the joint attack.

##  Appendix A Comparing the Joint and Sequential Attacks

Extending on Table 2, Table 3 further compares the joint and sequential
attacks of the symbiotic threat model. Although the joint attack is stronger
than the sequential attack in a higher number of cases, they are often within
one standard error of each other so it is difficult to argue for a
statistically significant difference. This might indicate the difficulty of
estimating the future evasion perturbations from within the poisoning attack
due to the very large search space, and highlight the space for potential
developments of stronger attacks.

##  Appendix B Inner Evasion Iterations.

The number of iterations the evasion attack runs for within the joint attack
is also configurable. Figure 5 displays the accuracy of a GCN on our three
benchmark graphs after the joint attack with different number of inner-evasion
iterations. We observe, contrary to our expectations, that increasing the
number of iterations for the inner-evasion attack has no significant influence
on the perturbed accuracy. This perhaps highlights a potential point of
improvement for future symbiotic attacks.

Table 4: Inductive Model | Dataset | Clean | Evasion | Poisoning | Symbiotic
---|---|---|---|---|---
GCN | CiteSeer | 0.67 ± 0.01 | 0.41 ± 0.01 | 0.62 ± 0.01 | 0.33 ± 0.01
Cora | 0.75 ± 0.02 | 0.42 ± 0.01 | 0.68 ± 0.03 | 0.3 ± 0.01
GAT | CiteSeer | 0.68 ± 0.01 | 0.37 ± 0.01 | 0.64 ± 0.02 | 0.56 ± 0.02
Cora | 0.77 ± 0.01 | 0.21 ± 0.01 | 0.61 ± 0.04 | 0.35 ± 0.03
APPNP | CiteSeer | 0.71 ± 0.01 | 0.47 ± 0.01 | 0.66 ± 0.02 | 0.4 ± 0.01
Cora | 0.82 ± 0.02 | 0.53 ± 0.02 | 0.78 ± 0.01 | 0.37 ± 0.01
GPRGNN | CiteSeer | 0.67 ± 0.01 | 0.37 ± 0.01 | 0.56 ± 0.01 | 0.34 ± 0.01
Cora | 0.8 ± 0.02 | 0.44 ± 0.01 | 0.74 ± 0.01 | 0.35 ± 0.01