# Open-Ended Learning Strategies for Learning Complex Locomotion Skills

Fangqin Zhou
<EMAIL_ADDRESS>
&Joaquin Vanschoren
<EMAIL_ADDRESS>

###### Abstract

Teaching robots to learn diverse locomotion skills under complex three-
dimensional environmental settings via Reinforcement Learning (RL) is still
challenging. It has been shown that training agents in simple settings before
moving them on to complex settings improves the training process, but so far
only in the context of relatively simple locomotion skills. In this work, we
adapt the Enhanced Paired Open-Ended Trailblazer (ePOET) approach to train
more complex agents to walk efficiently on complex three-dimensional terrains.
First, to generate more rugged and diverse three-dimensional training terrains
with increasing complexity, we extend the Compositional Pattern Producing
Networks - Neuroevolution of Augmenting Topologies (CPPN-NEAT) approach and
include randomized shapes. Second, we combine ePOET with Soft Actor-Critic
off-policy optimization, yielding ePOET-SAC, to ensure that the agent could
learn more diverse skills to solve more challenging tasks. Our experimental
results show that the newly generated three-dimensional terrains have
sufficient diversity and complexity to guide learning, that ePOET successfully
learns complex locomotion skills on these terrains, and that our proposed
ePOET-SAC approach slightly improves upon ePOET.

Keywords:  Reinforcement Learning, POET, ePOET, Locomotion Skills, Legged
Robots, CPPN-NEAT

##  1 Introduction

In recent years, Reinforcement Learning (RL) and Deep RL (DRL) have achieved
remarkable successes in the area of legged robot locomotion, especially in
controlling robots or agents to successfully walk on flat terrains [15, 16,
12]. However, many of those agents usually fail to maintain well-balanced
behaviors on rugged terrains [17, 9], because locomotion on uneven terrains
requires the ability to perceive the environment. Therefore, dynamically
generating diverse terrains for the robots to interact with is one of the key
challenges in controlling robots on complex terrains. To address this
challenge, Wang et al. [25] introduced an approach that automatically
generates diverse environments while optimizing the policy, called the Paired
Open-Ended Trailblazer (POET). Moreover, an improved version, called Enhanced
POET (ePOET) [26], can generate more diverse and complex challenges by
leveraging compositional pattern producing networks (CPPNs) [19] as a terrain
encoding method, coupled with the Neuroevolution of Augmenting Topologies
(NEAT) [20, 21] algorithm to evolve increasingly rugged environments. Despite
achieving impressing results in training a simple bipedal walker walking over
two-dimensional terrains, the question remains whether this approach also
works for complex agents in three-dimensional settings. Modeling gait
transitions of legged robots is more complex and much harder compared to biped
walkers. Additionally, generating three-dimensional terrains with gradually
increasing complexity also requires a more careful and complex design in
comparison to two-dimensional worlds. By addressing these issues, we aim to
make this work more applicable to real-world cases.

To this end, this work aims to make the following contributions: 1)
automatically generating diverse three-dimensional complex terrains while
optimizing policies; 2) analyzing the effectiveness of ePOET in training a
complex hexapod robot to adapt suitable gaits to these terrains; 3) improve
upon ePOET so that the agent learns more diverse skills to solve more
challenging tasks.

This paper is organized as follows. First, section 2 presents a detailed
literature review of previous work on complex locomotion, as well as the POET
and the ePOET approach. We then describe our proposed ePOET-SAC approach in
section 3. Section 4 describes how we leverage CPPN-NEAT [20, 21] to generate
complex terrains. Our conducted experiments and experimental results are
presented in section 5. Lastly, research directions are proposed for future
study in section 6.

##  2 Related Work

Complex Locomotion in simulation: Levine et al. [11] applied direct policy
search methods to produce motions for a bipedal walker running on uneven
terrains. Particularly, their approach succeeded in learning a push recovery
behavior. Heess et al. [10] applied stochastic policy gradient methods in
learning continuous control policies for locomotion agents to address two
potential limitations (reliance on planning and restriction to deterministic
models) caused by value gradient methods. A 3-dimensional physics-based
locomotion controller is presented by Mordatch et al. [14]. Their experiments
demonstrated that by optimizing a low-dimensional physical model, their
controller can successfully traverse not only flat terrains but also many
types of uneven and constrained terrain. However, their controllers failed to
consistently handle certain types of motion such as generating $180^{o}$ turns
when the character is running quickly. Heess et al. [9] used Proximal Policy
Optimization (PPO) and distributed PPO to explore complex locomotion
behaviours over a wide range of environmental conditions. Furthermore,
imitation learning has been used for producing high-quality locomotion through
imitating well-defined experts’ behaviors, such as work done by Chentanez et
al. [3]. Peng et al. [16] proposed a mixture of actor-critic experts (MACE)
model, to work directly with high-dimensional character and terrain state
descriptions without requiring feature engineering. More recently, Luo et al.
[13] applied Generative Adversarial Networks (GAN) to enable a high-level
gating network to approximate previously learned natural action distributions.
Learning locomotion on various terrains can be treated as a problem of
learning several skills for a variety of terrains and choosing the right skill
for a certain type of terrain. From this point of view, a two-layer recurrent
policy combining with PPO was introduced by Azayev et al. [1] to train a
hexapod walker to adapt to different terrains. To improve sample efficiency
and generalization, a hierarchical RL structure that combines an off-policy
Soft Actor-Critic method with a model-based planning approach is proposed by
Li et al. [12]. Their approach helped a Daisy robot to reach goals up to 12
meters away from its start point and to follow the waypoints defined by a
user.

Limitations: Despite the outstanding results that previous studies have
achieved, to our knowledge, they have shown two primary limitations. First of
all, their experiments were performed with mostly manufactured or randomly
generated terrains that either lack complexity and diversity or have uneven
difficulty levels [1, 9, 16, 11, 14, 15, 12]. With very simple terrains, an
agent could only acquire some very basic skills, while with overly challenging
terrains, an agent could fail to learn at the very beginning. Therefore,
terrains with well-balanced and gradually increasing difficulty levels are
essential for an agent to learn diverse locomotion skills. Secondly, some
algorithms are unstable for training legged robots’ locomotion tasks and are
highly sensitive to hyper-parameter settings, such as the actor-critic
architecture [12, 16]. An unstable training could cause the trained agent to
behave unnaturally and inefficiently, especially when the agents face uneven
training environments.

Paired Open-Ended Trailblazer (POET), proposed by Wang et al. [25], aims to
directly confront open-endedness, i.e., an unbounded invention of learning
environments and their solutions. This is done by evolving a set of diverse
and increasingly complex environmental challenges while collectively
optimizing their solutions. These environmental challenges and solutions
together form a class of environment-agent (EA) pairs. When agents prove
successful in one environment, they are transferred to another, usually more
complex environment. This exploits the opportunity to transfer high-quality
solutions from one objective to another. Each EA pair is optimized with
Evolution Strategies (ES) [8, 18, 28, 27]. Enhanced POET (ePOET) [26] improves
upon POET in several ways. Firstly, a domain-general Environment
Characterization (EC), called the Performance of All Transferred Agents EC
(PATA-EC), is used to evaluate how all agents perform in each newly generated
environment instead of relying on domain-specific information. Secondly, ePOET
simplifies the transfer mechanism of the original POET by introducing a more
stringent threshold to save computing resources. Finally, another environment
encoding mechanism, based on compositional pattern producing networks (CPPNs)
[19] is applied to generate environments with increasing complexity and
express any possible landscape at any conceivable resolution or size. More
details about CPPN are described in section 4.1.

Although both POET and ePOET have shown exciting results, they have only been
evaluated within a simple 2D bipedal walker environment. In three-dimensional
settings, both the observation space and action space could have more
dimensions, so training could also be much harder. Hence, it is unclear
whether this approach generalizes to more complex settings, such as 3D agents
walking over complex three-dimensional terrains. To address this question, we
modify the ePOET algorithm to perform experiments in 3D space and extend the
CPPN approach to generate more rugged terrains. Moreover, to address the
problems of local convergence and mutation-sensitivity when using evolution
strategies, we combine ePOET with the Soft Actor-Critic approach to encourage
more random behaviors while maximizing returns.

##  3 Methodology

###  3.1 Preliminaries

We define an RL model as
$M=\left(\mathcal{S},\mathcal{A},\mathcal{T},R,\gamma\right)$, in which
$\mathcal{S}$ represents a set of states, $\mathcal{A}$ denotes a set of
actions, $\mathcal{T}$ is a transition probability distribution, $R$ is the
reward function, and $\lambda\in[0,1]$ is a discount factor. Under this
definition, at each environment-step $t$, the agent observes a state $s_{t}$
from its state space $\mathcal{S}$, then takes an action $a_{t}$ from its
action space $\mathcal{A}$ according to a policy $\pi(a=a_{t}\mid s=s_{t})$,
and then receives a corresponding reward $r(s=s_{t},a=a_{t})$ calculated with
the taken action and its reward function $R$. Meanwhile, the transition
probability
$P_{a}\left(s,s^{\prime}\right)=\operatorname{Pr}\left(s_{t+1}=s^{\prime}\mid
s_{t}=s,a_{t}=a\right)$ is calculated and the corresponding transition
$T(s_{t},a_{t},s_{t+1})$ is added into the transition probability distribution
$\mathcal{T}$. This indicates the probability of changing the state from the
current $s$ to the next state $s^{\prime}$ by taking action $a$ at
environment-step $t$. The goal of an RL model is to seek an optimal policy
$\pi^{*}$ that maximizes the expected cumulative reward $J(\pi)$, i.e.,
$\pi^{*}=\underset{\pi}{\arg\max}\hskip 2.84544ptJ(\pi)$.

###  3.2 Motivation

The evolution strategies applied in ePOET have good scalability but suffer
from local convergence [18, 4], while SAC has less stability but encourages
more diverse motions through maximizing entropy. Hence, combining ES and SAC
promises to produce more diverse motions that enable the agent to overcome
more complex challenges while retaining the scalability of ES. This idea is
inspired by Suri et al. [22], who introduced the Evolution-based Soft Actor-
Critic (ESAC) method. This method adds a SAC actor into a population of actors
optimized using evolutions to exploit gradient-based knowledge. They
demonstrated that it combines high sample efficiency and scalability. Since
each agent-environment pair of ePOET is optimized with ES, it could be
similarly combined with SAC to encourage a more diverse exploration of each
agent. However, different from ESAC, which optimizes a single ES optimizer, in
our method we add the SAC actor to the ePOET agent pool to form a new ES
optimizer, and then perform crossover between the newly added optimizer and
all other ES optimizers in the pool.

We name this combination ePOET-SAC. The workflow of this approach is
illustrated in Figure 1, and detailed in Algorithm 1. The steps highlighted in
red are additional steps in comparison to ePOET. Initially, similar to ePOET,
an environment is generated together with a randomly initialized agent
(parameter vector), and the agent is to be optimized under the environment via
an ES optimizer. Simultaneously, the SAC model is initialized with one policy
network and two Q-networks (as two Q-networks can significantly speed up
training [7]), as well as a global empty replay buffer for storing old
experiences. Meanwhile, an environment is randomly sampled from the
environment pool to train the SAC model. If there is only one initial
environment in the pool, then the same environment is used for the SAC model.
At each iteration, the ePOET and the SAC are trained in parallel, and the SAC
networks are updated every environment step as explained in Appendix A.1.

Figure 1: Illustration of Algorithm 1. ePOET and SAC are trained
simultaneously. Trained SAC agents are added to the ePOET agent pool, and
crossover is performed with all active ePOET agents. The original agents are
replaced when they are outperformed by the newly generated agents.

Input: initial environment $E^{\text{init }}(\cdot)$, its paired agent’s
policy parameter vector $\theta^{\text{init }}$, learning rate $\alpha$, noise
standard deviation $\sigma$, iterations $T$, mutation interval
$N^{\text{mutate }}$, transfer interval $N^{\text{transfer }}$

Initialize _Set EA $\\_$list empty _

Initialize _SAC, an empty replay buffer R (global variable)_

Add $\left(E^{\text{init }}(\cdot),\theta^{\text{init }}\right)$ to
EA$\\_$list

for  _$t=0$ to $T-1$ _ do

if  _$t >0$ and $t$ mod $N^{\text{mutate }}=0$ _ then

EA$\\_$list = MUTATE$\\_$ENVS(EA$\\_$list) $\\#$new environments created by
mutation

end if

$M=\operatorname{len}(\mathrm{EA}_{-}$list$)$

for  _$m=1$ to $M$ _ do

$E^{m}(\cdot),\theta_{t}^{m}=\mathrm{EA}_{-}$list$[\mathrm{m}]$

$\theta_{t+1}^{m}=\theta_{t}^{m}+$ ES$\\_$STEP
$\left(\theta_{t}^{m},E^{m}(\cdot),\alpha,\sigma,\right)\quad$ $\\#$ each
agent independently optimized

end for

Train$\\_$SAC()  $\\#$ Algorithm 2 (see appendix)

for  _$m=1$ to $M$ _ do

if  _$M >1$ and $t$ mod $N^{\text{transfer }}=0$ _ then

$\theta^{\text{top }}=$ EVALUATE$\\_$CANDIDATES
$\left(\theta_{t+1}^{1},\ldots,\theta_{t+1}^{m-1},\theta_{t+1}^{m+1},\ldots,\theta_{t+1}^{M},E^{m}(\cdot),\alpha,\sigma\right)$

if  _$E^{m}\left(\theta^{\text{top }}\right)
>E^{m}\left(\theta_{t+1}^{m}\right)$ _ then

$\theta_{t+1}^{m}=\theta^{\text{top }}$ $\\#$ transfer attempts

end if

end if

EA$\\_$list$[m]=\left(E^{m}(\cdot),\theta_{t+1}^{m}\right)$

end for

if  _$M >1$ and $t$ mod $N^{\text{transfer }}=0$ _ then

if  _Score(SAC $\\_$Actor) $<$ Score($\theta^{\text{top }}$)_ then

SAC$\\_$Actor = $\theta^{\text{top }}$ $\\#$update SAC actor with the best
agent in the pool

end if

end if

if  _$t >0$ and $t$ mod ($N^{\text{transfer }}*4)=N^{\text{transfer }}*4-2$ _
then

Add SAC$\\_$Actor to the agent pool as $\theta_{t}^{sac}$

Crossover(agents list, $\theta_{t}^{sac}$)  $\\#$ Algorithm 3 (see appendix)

end if

end for

Algorithm 1 ePOET-SAC Main Loop

After training for a user-specified number of iterations, the updated SAC
Actor is added into the agent pool and crossover is performed with the active
agents in the pool. The crossover subroutine is illustrated in Algorithm 3 of
Appendix A.1. In essence, during crossover, we randomly exchange the weights
and biases of two agents (two input parameter vectors $\theta^{m}$ and
$E^{m}$): for each node of each layer, we flip a coin and randomly replace the
weight of $\theta^{m}$ with the weight of $E^{m}$ from the same position. We
then perform the same operation for all biases. Before performing every single
crossover, we evaluate the candidate agent $\theta^{m}$ and save its score.
Then, we evaluate the agent again after performing the single crossover. If
the evaluation score of an agent after crossover is higher than its original
score, then it is used to replace the original. Otherwise, we keep the
original agent and discard the other.

Simultaneously, at every environment step, the transitions
$<s,s^{\prime},a,r,d>$ of all active agents in the pool are pushed into the
replay buffer from which the SAC samples previous experiences. Hence, the
replay buffer contains not only transitions from the SAC model itself but also
old experiences from the agents in the pool. This is because, instead of
reusing only the SAC actor’s experiences, learning from more agents in the
pool could help improve the SAC actor’s performance and speed up training.
Each transition pushed into the replay buffer consists of the current state
$s$, the next state $s^{\prime}$, the current action $a$, the reward $r$ and
the terminating indicator $d$ after taking action $a$. Furthermore, after
every transfer step, the best parameter vector (best agent in the pool) is
also transferred to the SAC actor if it outperforms the SAC actor. This is to
speed up the training of SAC.

##  4 Environment Evolution

This section explains how we combine the Compositional Pattern Producing
Networks (CPPNs) encoding with the Neuroevolution of Augmenting Topologies
(NEAT) algorithm to produce diverse three-dimensional terrains with gradually
increasing complexity.

###  4.1 CPPN-NEAT

Compositional Pattern Producing Networks (CPPNs) are an indirect encoding
method that abstracts the process of natural development without requiring the
simulation of diffusing chemicals but through a composition of functions [19].
These functions are similar to the activation functions of Artificial Neural
Networks (ANNs), however, unlike ANNs, CPPNs could allow each node to select a
unique activation function so that the output of each is patterned,
symmetrical and unique. Due to the structural similarity with ANNs, CPPNs
could be evolved through neuroevolution algorithms such as the Neuroevolution
of Augmenting Topologies (NEAT) algorithm [20, 21], together called CPPN-NEAT.
The NEAT algorithm is an evolutionary algorithm that could evolve increasingly
complex networks over generations by adding nodes and more connections, or
deleting existing connections in the population. In short, the idea behind
CPPN-NEAT is that the CPPNs take geometric coordinates as inputs and output
expression patterns that describe the phenotypes, and through evolving the
networks over generations by the NEAT algorithm, increasingly complex
phenotype expression patterns can be produced. More specifically, in the
beginning, the CPPNs are initialized with random simple structures without
hidden nodes. By evolving through the NEAT method over generations, extra
nodes and connections are added into the networks, and new offspring
(networks) are generated. Hence, the population of CPPNs becomes more complex
as evolution continues and topological mutations are applied.

###  4.2 Generating three-dimensional terrains with CPPN-NEAT

In ePOET [26], the authors used CPPNs to produce y-coordinates for each given
x-coordinate and then render the corresponding (x,y) coordinates into a two-
dimensional terrain. The CPPN was updated by NEAT every user-specified number
of iterations. Similarly, in this paper, our goal is to adapt this idea to
generate more complex rich three-dimensional terrains. Figure 2 illustrates
the overall idea of generating three-dimensional terrains with CPPN-NEAT. The
terrain is generated with a heightmap, which consists of two parts, namely a
random bowl-shaped part and the CPPN-NEAT part. The random bowl shape is
created through sampling points from a uniform distribution and applying a
cosine function to these sampled points, and then subtracting the cosine value
with a gradually increased threshold (the maximum height). An example terrain
generated from the random bowl shape method can be found in Figure 4 of
Appendix A.2.

Figure 2: An overview of generating three-dimensional terrains. It contains
two parts: the random bowl shape (left) and the CPPN-NEAT part (right), with a
weight of 0.3 and 0.7, respectively.

The second part is to generate more diverse terrains through CPPN encoding and
to evolve both the weights and the architecture of CPPN through the NEAT
algorithm. We define two hidden nodes for the initial network to make it
neither over-complex nor over-simple. The $(x,y)$ coordinates in the three-
dimensional plane are taken as the inputs of the CPPN. The output of each node
is determined by formula:
$activation(bias+weight*(response*aggregation(inputs)))$, in which the
activation function is chosen from
$\\{sin,sigmoid,square,tanh,identity,gauss\\}$, the aggregation function is
set to $sum$, and the $bias,weight,response$ parameters are updated through
the NEAT algorithm. More settings are listed in Table 6 of Appendix A.5. With
these settings, the network produces varied height values for each point
$(x,y)$ in the plane, and the produced heights together with the heights
generated from the random bowl shape form different height maps. The height
maps are used to generate terrains in MuJoCo [23]. Importantly, the fitness
function and the fitness threshold of the NEAT algorithm are not used in this
approach because we applied the PATA-EC method also used in ePOET to guarantee
the novelty and complexity of every newly generated terrain. The key points
that enable the generated terrains to have sufficient diversity and gradually
increasing complexity are illustrated in Appendix A.2, and some generated
example terrains can be seen in Figure 5 of Appendix A.2.

##  5 Experimental Results

We evaluate our method using a Hexapod walker created by Azayev et al [1] with
18 Degrees-of-Freedom (DoF), which allows it to be more flexible on uneven
terrains and to have a suitable training complexity. The agent and its
implementation are detailed in Appendix A.3 and A.4, respectively. We trained
the hexapod agent with PPO, SAC, VMPO, ePOET, and ePOET-SAC approaches to
allow a direct comparison. All the experiments are modeled with the MuJoCo
[23] physics simulator. The ePOET and our proposed ePOET-SAC ran on 10 CPU
workers (Intel Xeon Broadwell-EP 2683v4 @ 2.1GHz), while PPO, VMPO, and SAC
ran on GPU (ASUS Turbo GeForce GTX 1080 Ti, CUDA 10.2). The hyper-parameter
settings of each algorithm are presented in Appendix A.5.

The average training returns of ePOET, PPO, SAC, and VMPO are shown in Figure
3 (a), showing that ePOET performs well compared to other methods. In this
figure, ePOET$\\_$gen0 denotes the initial agent of ePOET, and
ePOET$\\_$gen$i$ denotes the paired agents of reproduced environments from
mutating, and the number $i$ indicates their generation order. For the
training terrains, because the CPPN-NEAT encoding method cannot be used in
PPO, SAC, and VMPO, their training terrains contain only the random bowl shape
with a maximum height of 1. These are actually easier than the terrains
generated for ePOET because CPPN-NEAT generates terrains with larger height
variation. In addition, the terrain’s complexity of agent ePOET$\\_$gen$i$
increases with their generation order. Furthermore, Figure 3 (b) shows all
generations of ePOET-SAC in 24 hours running. Because we set the interval of
adding SAC actors to the agent pool as $iterations\%(15*4)=15*4-2$, the ePOET-
SAC$\\_$gen1 (yellow) in Figure 3 represents the SAC actor added at iteration
58 and the ePOET-SAC$\\_$gen2 (green) represents the SAC actor added at
iteration 118. Through the transfer across generations, we see that ePOET-
SAC$\\_$gen1 (yellow) outperforms ePOET-SAC$\\_$gen0 (blue) after some
iterations and ePOET-SAC$\\_$gen0 (blue) is also improved dramatically by
transferring knowledge from ePOET-SAC$\\_$gen1 (yellow) or ePOET-SAC$\\_$gen2
(green).

(a) The training of ePOET, SAC, PPO, VMPO.

(b) The training of ePOET-SAC.

Figure 3: (a) Comparisons among ePOET, SAC, PPO, and VMPO showing the average
training return over the number of iterations. ePOET$\\_$gen$i$ refers to the
ePOET agents of the $i$th generation. (b) All generated agents of ePOET-SAC
running in 24 hours with 10 workers.

To measure how well the trained agents perform when encountering unseen
terrains, we evaluate the hexapod walker trained by PPO, SAC, ePOET, and
ePOET-SAC with the same 32 environments. These 32 terrains were generated by
ePOET from a run that is different from the above ePOET agents, which means
that the terrains are not only new for PPO and SAC but also new for the ePOET
and ePOET-SAC agents. Importantly, in our settings, an eligible parent
environment generates at most 8 child environments and only one of them is
admitted, which is to prevent generating overly challenging or overly easy
terrains (as explained in section 3.2). Therefore, only 4 environments among
these 32 environments are ensured to be neither very difficult nor very easy,
and some environments could be overly challenging. We took all these 32
environments to compare the performance on not only easy terrains but also
challenging ones. Furthermore, since the agent’s travel route of each run
could be different, we run the agent 5 times on each terrain to get a better
performance estimate. Thus, in total, each trained agent runs on $32*5=160$
terrains. We consider an environment ’solved’ if the agent yields a score of
2000 or higher.

Table 1: Success rate (percentage of environments reaching a score of 2000 or higher) across 32 generated environments. Each environment is run for 5 times by each agent. Algorithm | Easy | Medium difficulty | Hard | Total
---|---|---|---|---
PPO | 72% | 26% | 0% | 36%
SAC | 78% | 40% | 0% | 43%
ePOET |  91% | 79% | 0% | 64%
ePOET-SAC |  91% |  82% |  3% |  66%

The comparison results are shown in Table 1. As can be seen, ePOET passed more
terrains than PPO and SAC in both easy and middle terrains. In general, ePOET
performed better than PPO and SAC, and none of these approaches can solve hard
terrains. In contrast, ePOET-SAC did solve a few of the hard-level terrains,
and it solved more terrains than ePOET in total. Interestingly, the hexapod
trained with these approaches behaves very differently under the same
environment, which can be seen in the video.111Video comparing learned
locomotion skills: https://youtu.be/YZ2rbdW29r0 The hexapod trained with SAC
can solve an environment that the hexapod trained with PPO cannot solve.
However, the locomotion learned with SAC is often less efficient than the one
learned with PPO. In comparison to the agent trained with ePOET, the agent
trained with ePOET-SAC shows the best balanced and efficient walking styles on
different terrains. This comparison is shown in a second video.222Video
comparing ePOET and ePOET-SAC: https://youtu.be/i8QM1-FoAto To gain more
insight in their behavioural differences, we evaluate the trained agents on 16
easy and medium difficulty environments. Table 2 summarizes their average
evaluation scores over 5 runs on the 16 unknown terrains. As can be seen,
ePOET-SAC obtained both the highest average and maximum returns, yet very
close to ePOET. More behavior analysis are shown in Appendix A.6.

Table 2: Average Evaluation over 16 environments. Algorithm | Average Return | Maximum Return
---|---|---
PPO | 564.63±1669.03 | 3144.96
SAC | 1657.44±2981.49 | 3387.95
ePOET | 4030.22±246.29 | 4421.11
ePOET-SAC | 4070.78±248.37 | 4473.74

##  6 Conclusions & Future Scope

#### Contributions

The contributions of this work center around three key aspects. First, our
experimental results have shown that the ePOET approach can indeed outperform
some classic RL algorithms such as PPO, SAC, and VMPO in acquiring diverse
locomotion skills in complex three-dimensional environmental settings. Second,
our adapted approach of combining the CPPN-NEAT approach with random bowl
shapes can generate diverse three-dimensional terrains with gradually
increasing complexity. Third, our proposed ePOET-SAC approach slightly
outperforms the ePOET algorithm, especially on hard terrains, by combining
ePOET with SAC so that the trained agent could learn more diverse locomotion
skills and overcome more challenging terrains.

#### Limitations

There are still some limitations in our work. First, the diversity of
generated terrains could be further improved. For instance, adding small deep
gaps, continuous stairs, or obstacles. This might be achieved with a better
design or more careful fine-tuning of CPPN-NEAT. Secondly, due to the hidden
layer shape differences between SAC ($[256,256]$) and ES ($[40,40]$) of ePOET
(as described in Appendix A.5), we randomly choose parameters to reshape the
parameter vectors when adding the SAC actor into the ES agent pool. This
random selection could weaken the effectiveness of the SAC actor to the ES
agents when performing crossover. In our analysis, we found that the
improvement of our ePOET-SAC compared to the ePOET mainly benefits from adding
SAC actors into the pool rather than performing crossover. This is because the
replacement of the existing agents with the new ones generated from crossover
happens only when crossover outperforms the original agents. To address this
problem, we attempted to use a middle hidden layer shape of $[128,128]$, but
this led to degraded performance. Further research is needed to address these
challenges.

#### Future Scope

The potential future scope would focus on addressing the above limitations.
Firstly, more careful fine-tuning of the hyper-parameters of CPPN-NEAT would
be desired, and a method of monitoring the difficulty levels of generated
terrains would be helpful for providing more insight into the generated
environments. Secondly, based on the fact that our ePOET-SAC does not
massively outperform ePOET, tuning the entropy parameter of SAC could help to
improve the exploration process. Furthermore, the training algorithm could be
further improved with respect to two possible aspects. One is to pre-train the
agent with a Diversity is All You Need (DIAYN) [6] approach to learn useful
skills and then to adapt the learned skills to different terrains. Because
learning locomotion skills is highly sensitive to reward function design, and
human-designed reward signals could bring about potential limitations,
learning without reward functions would allow us to address this problem.
Another possible future direction is to make use of additional meta-learning
approaches. Because the agent encounters new terrains every iteration, a meta-
learner (running in an outer loop) could help to first identify the current
terrain type (e.g. roughness, steepness, depth, height, etc) and then apply
suitable locomotion gaits to the corresponding terrain. Moreover, meta-
learning can also be used to learn good initial policy weights or hyper-
parameter values over a distribution of environments. Finally, using the CPPN-
NEAT approach to evolve the agent so that the agent can flexibly adapt its
body to different terrains might also be an interesting and possible research
direction.

## Acknowledgements

This research was partially supported by TAILOR, a project funded by EU
Horizon 2020 research and innovation programme under GA No 952215.

## References

  * [1] Teymur Azayev and Karel Zimmerman.  Blind Hexapod Locomotion in Complex Terrain with Gait Adaptation Using Deep Reinforcement Learning and Classification.  Journal of Intelligent & Robotic Systems, 99, 09 2020.
  * [2] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  Openai gym, 2016.
  * [3] Nuttapong Chentanez, Matthias Müller, Miles Macklin, Viktor Makoviychuk, and Stefan Jeschke.  Physics-based motion capture imitation with deep reinforcement learning.  In Proceedings of the 11th Annual International Conference on Motion, Interaction, and Games, MIG ’18, New York, NY, USA, 2018. Association for Computing Machinery.
  * [4] Edoardo Conti, Vashisht Madhavan, Felipe Petroski Such, Joel Lehman, Kenneth O. Stanley, and Jeff Clune.  Improving Exploration in Evolution Strategies for Deep Reinforcement Learning via a Population of Novelty-Seeking Agents.  CoRR, abs/1712.06560, 2017.
  * [5] Prafulla Dhariwal, Christopher Hesse, Oleg Klimov, Alex Nichol, Matthias Plappert, Alec Radford, John Schulman, Szymon Sidor, Yuhuai Wu, and Peter Zhokhov.  Openai baselines.  https://github.com/openai/baselines, 2017.
  * [6] Benjamin Eysenbach, Abhishek Gupta, Julian Ibarz, and Sergey Levine.  Diversity is all you need: Learning skills without a reward function.  CoRR, abs/1802.06070, 2018.
  * [7] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  CoRR, abs/1801.01290, 2018.
  * [8] Nikolaus Hansen, Sibylle Müller, and Petros Koumoutsakos.  Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES).  Evolutionary computation, 11:1–18, 02 2003.
  * [9] Nicolas Heess, Dhruva TB, Srinivasan Sriram, Jay Lemmon, Josh Merel, Greg Wayne, Yuval Tassa, Tom Erez, Ziyu Wang, S. M. Ali Eslami, Martin A. Riedmiller, and David Silver.  Emergence of Locomotion Behaviours in Rich Environments.  CoRR, abs/1707.02286, 2017.
  * [10] Nicolas Heess, Greg Wayne, David Silver, Timothy P. Lillicrap, Yuval Tassa, and Tom Erez.  Learning Continuous Control Policies by Stochastic Value Gradients.  CoRR, abs/1510.09142, 2015.
  * [11] Sergey Levine and Vladlen Koltun.  Learning complex neural network policies with trajectory optimization.  In Eric P. Xing and Tony Jebara, editors, Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 829–837, Bejing, China, 22–24 Jun 2014. PMLR.
  * [12] Tianyu Li, Nathan O. Lambert, Roberto Calandra, Franziska Meier, and Akshara Rai.  Learning generalizable locomotion skills with hierarchical reinforcement learning.  CoRR, abs/1909.12324, 2019.
  * [13] Ying-Sheng Luo, Jonathan Hans Soeseno, Trista Pei-Chun Chen, and Wei-Chao Chen.  CARL: controllable agent with reinforcement learning for quadruped locomotion.  CoRR, abs/2005.03288, 2020.
  * [14] Igor Mordatch, Martin de Lasa, and Aaron Hertzmann.  Robust physics-based locomotion using low-dimensional planning.  In ACM SIGGRAPH 2010 Papers, SIGGRAPH ’10, New York, NY, USA, 2010\. Association for Computing Machinery.
  * [15] J. Morimoto, G. Cheng, C.G. Atkeson, and G. Zeglin.  A simple reinforcement learning algorithm for biped walking.  In IEEE International Conference on Robotics and Automation, 2004\. Proceedings. ICRA ’04. 2004, volume 3, pages 3030–3035 Vol.3, 2004.
  * [16] Xue Bin Peng, Glen Berseth, and Michiel van de Panne.  Terrain-adaptive locomotion skills using deep reinforcement learning.  ACM Trans. Graph., 35(4), July 2016.
  * [17] Branislav Pilňan.  Exploring dynamic locomotion of a quadruped robot: a study of reinforcement learning for the anymal robot.  2019\.
  * [18] Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever.  Evolution strategies as a scalable alternative to reinforcement learning, 2017.
  * [19] Kenneth Stanley.  Compositional pattern producing networks: A novel abstraction of development.  Genetic Programming and Evolvable Machines, 8:131–162, 06 2007\.
  * [20] Kenneth O. Stanley.  The hybercube-based neuroevolution of augmenting topologies: http://eplex.cs.ucf.edu/hyperneatpage/hyperneat.html.
  * [21] Kenneth O. Stanley.  The neuroevolution of augmenting topologies: https://www.cs.ucf.edu/ kstanley/neat.html.
  * [22] Karush Suri, Xiao Qi Shi, Konstantinos N. Plataniotis, and Yuri A. Lawryshyn.  Evolve To Control: Evolution-based Soft Actor-Critic for Scalable Reinforcement Learning.  CoRR, abs/2007.13690, 2020.
  * [23] Emanuel Todorov, Tom Erez, and Yuval Tassa.  Mujoco: A physics engine for model-based control.  In Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on, pages 5026–5033. IEEE, 2012.
  * [24] Laurens van der Maaten and Geoffrey Hinton.  Visualizing data using t-SNE.  Journal of Machine Learning Research, 9(86):2579–2605, 2008.
  * [25] Rui Wang, Joel Lehman, Jeff Clune, and Kenneth O. Stanley.  Paired open-ended trailblazer (POET): endlessly generating increasingly complex and diverse learning environments and their solutions.  CoRR, abs/1901.01753, 2019.
  * [26] Rui Wang, Joel Lehman, Aditya Rawal, Jiale Zhi, Yulun Li, Jeff Clune, and Kenneth O. Stanley.  Enhanced POET: Open-Ended Reinforcement Learning through Unbounded Invention of Learning Challenges and their Solutions, 2020.
  * [27] Daan Wierstra, Tom Schaul, Tobias Glasmachers, Yi Sun, Jan Peters, and Jürgen Schmidhuber.  Natural evolution strategies.  J. Mach. Learn. Res., 15(1):949–980, January 2014.
  * [28] Daan Wierstra, Tom Schaul, Jan Peters, and Jürgen Schmidhuber.  Natural evolution strategies.  pages 3381–3387, 06 2008.
  * [29] Rchal Yang.  Pytorch implementation of reinforcement learning methods : https://github.com/rchalyang/torchrl, 2021.
  * [30] Jiale Zhi, Rui Wang, Jeff Clune, and Kenneth O. Stanley.  Fiber: A platform for efficient development and distributed training for reinforcement learning and population-based methods, 2020.

##  Appendix A Appendix

###  A.1 Algorithms

The algorithm of training SAC for one iteration is shown in algorithm 2, in
which $\psi,\theta,\phi,\bar{\psi}$ represent the parameters of the value
network, the Q-networks, the policy network, and the target value networks,
respectively. The gradient of the soft value function $J_{V}(\psi)$ can be
estimated with an unbiased estimator:

$\hat{\nabla}_{\psi}J_{V}(\psi)=\nabla_{\psi}V_{\psi}\left(\mathbf{s}_{t}\right)\left(V_{\psi}\left(\mathbf{s}_{t}\right)-Q_{\theta}\left(\mathbf{s}_{t},\mathbf{a}_{t}\right)+\log\pi_{\phi}\left(\mathbf{a}_{t}\mid\mathbf{s}_{t}\right)\right).$

The soft Q-values are approximated by two parameterized functions
$Q_{\theta_{1}}(s,a)$ and $Q_{\theta_{2}}(s,a)$, and they can be optimized
with stochastic gradients:

$\hat{\nabla}_{\theta_{i}}J_{Q}(\theta_{i})=\nabla_{\theta_{i}}Q_{\theta_{i}}\left(\mathbf{a}_{t},\mathbf{s}_{t}\right)\left(Q_{\theta_{i}}\left(\mathbf{s}_{t},\mathbf{a}_{t}\right)-r\left(\mathbf{s}_{t},\mathbf{a}_{t}\right)-\gamma
V_{\bar{\psi_{i}}}\left(\mathbf{s}_{t+1}\right)\right)\hskip
5.69046ptfor\hskip 2.84544pti\in\\{1,2\\},$

in which $V_{\bar{\psi_{i}}}$ represents corresponding target value networks.
The target networks are updated through:

$\bar{\psi_{i}}\leftarrow\tau\psi+(1-\tau)\bar{\psi_{i}}\hskip
5.69046ptfor\hskip 2.84544pti\in\\{1,2\\}.$

One can refer to the original SAC paper [7] for more proof about the above
equations. One difference from the paper is that the replay buffer contains
not only the old experiences of SAC itself but also experiences from the
agents trained with ES. Additionally, we update the networks every environment
step as the authors [7] set both their target update interval and gradient
steps as 1 in their practice.

Input: global replay buffer $\mathcal{R}$

Initialize _parameter vectors $\psi,\bar{\psi},\theta,\phi$_

Initialize _learning rate
$\lambda_{V},\lambda_{Q},\lambda_{\pi},\lambda_{\alpha},\tau$_

for  _each environment step t_ do

select action
$\mathbf{a}_{t}\sim\pi_{\phi}\left(\mathbf{a}_{t}\mid\mathbf{s}_{t}\right)$

observe new state $\mathbf{s}_{t+1}\sim
p\left(\mathbf{s}_{t+1}\mid\mathbf{s}_{t},\mathbf{a}_{t}\right)$

$\mathcal{R}\leftarrow\mathcal{R}\cup\left\\{\left(\mathbf{s}_{t},\mathbf{a}_{t},r\left(\mathbf{s}_{t},\mathbf{a}_{t}\right),\mathbf{s}_{t+1},\mathbf{d}_{t}\right)\right\\}$

update value network
$\psi\leftarrow\psi-\lambda_{V}\hat{\nabla}_{\psi}J_{V}(\psi)$

update action value Q-networks
$\theta_{i}\leftarrow\theta_{i}-\lambda_{Q}\hat{\nabla}_{\theta_{i}}J_{Q}\left(\theta_{i}\right)$
for $i\in\\{1,2\\}$

update policy network
$\phi\leftarrow\phi-\lambda_{\pi}\hat{\nabla}_{\phi}J_{\pi}(\phi)$

update target value networks
$\bar{\psi_{i}}\leftarrow\tau\psi+(1-\tau)\bar{\psi_{i}}$ for $i\in\\{1,2\\}$

update temperature
$\alpha\leftarrow\alpha-\lambda_{\alpha}\nabla_{\alpha}J(\alpha)$

end for

Algorithm 2 Train$\\_$SAC

Input: a list of active parameter vectors
$\left(\theta^{1},\ldots,\theta^{m-1},\theta^{m+1},\ldots,\theta^{M}\right)$,
a list of active environments
$\left(E^{1},\ldots,E^{m-1},E^{m+1},\ldots,E^{M}\right)$, $\theta_{t}^{sac}$

for  _$m=1$ to $M$ _ do

old$\\_$return$\\_$m = rollout($\theta^{m},E^{m}$)

$new\\_\theta^{m}$ = crossover$\\_$single($\theta^{m}$, $\theta_{t}^{sac}$)

new$\\_$return$\\_$m = rollout($new\\_\theta^{m},E^{m}$)

if  _old $\\_$return$\\_$m $<$ new$\\_$return$\\_$m _ then

archive $\theta^{m}$

$\theta^{m}$ = $new\\_\theta^{m}$

end if

end for

Algorithm 3 Crossover

###  A.2 Generated Environments

Following the ePOET, crossover between CPPNs is also not performed in our
approach of generating three-dimensional terrains. Behind this approach, five
key points enable the generated terrains to have sufficient diversity and
gradually increasing complexity:

  * •

The CPPN-NEAT is set to be updated (produce new generations) after a user-
specified number of iterations (e.g. 75 is set in our experiments) and when
the agent has achieved a user-defined score threshold (e.g. 2000 is set in our
experiments) under the current encoded terrain. This is to ensure that the
agent has solved the current challenge before producing a new challenge level.
It also prevents the CPPN from being updated too often and generating overly
complex terrains.

  * •

The domain-general PATA-EC method, which is described section
LABEL:subsec:existingapp, is used to measure the novelty of every newly
generated terrain so that the newly added terrains are distinguished from the
existing ones in terms of diversity and complexity.

  * •

The part of the random bowl-shape is re-generated every new iteration, to
ensure that the agent is always trained on changed terrains every new
iteration even if the CPPN-NEAT is not updated.

  * •

The weight for the two components is set as 0.3 and 0.7, respectively, which
makes a good balance for the generated terrains based on our experiments.
Because the major complexity and diversity should be contributed by the CPPN-
NEAT part, whereas the random bowl shape part with a relatively small portion
(0.3) introduces small changes when the CPPN-NEAT is not updated.

  * •

The maximum terrain height ($elevation\\_z$ in MuJoCo [23]) is increased by
0.01 every iteration, which also ensures that the complexity of the generated
terrains gradually increase.

Figure 4 shows a terrain generated from the random bowl shape method with the
maximum height ($elevation\\_z$ in MuJoCo [23]) of 1.5.

Figure 4: A three-dimensional terrain that is generated from the random bowl
shape method.

Some generated example terrains can be seen in Figure 5. They are sorted into
five levels (easy, challenging, solvable, hard, and extremely hard) according
to their average height variances. The number of solved terrains by the
hexapod walker among all generated terrains is explained in section 5.
Moreover, Figure 6 illustrates that the complexity of created terrains over
generations gradually increases. These environments were generated by ePOET
and ePOET-SAC. The number corresponds to the generated order, for instance,
Env0 and Env11 correspond to the initial environment and the 11th environment
respectively. To some extent, the complexity of the generated environments
gradually increases as their generated orders.

Figure 5: Examples of generated terrains through our approach. The difficulty
level is decided by the average height variance.

(a) Env0 generated by ePOET

(b) Env1 generated by ePOET

(c) Env2 generated by ePOET

(d) Env4 generated by ePOET

(e) Env9 generated by ePOET

(f) Env10 generated by ePOET

(g) Env9 generated by ePOET-SAC

(h) Env11 generated by ePOET-SAC

Figure 6: Some generated environments by ePOET and ePOET-SAC.

###  A.3 Implementation

All our implementation is done with the Python language as it is widely used
in Machine Learning and Reinforcement Learning (RL) communities. Moreover,
Python provides many common packages for neural networks, such as the most
popular Tensorflow and Pytorch packages. There are also many RL baselines
implemented with Python, such as the baselines provided by OpenAI [5]. MuJoCo
simulator and mujoco_py API are used to simulate the agents, as well as OpenAI
Gym [2] benchmarks. For the RL baseline algorithms (PPO, SAC, VMPO), we used
Yang’s Pytorch implementation [29] and made small changes to adapt them to our
experiments. For the implementation of our approaches, we used the framework
of the ePOET that uses Fiber (a distributed computing library developed by
OpenAI) [30] for parallel processing. Fiber is similar to Python
Multiprocessing but more powerful. For instance, it has better ability in
terms of error handling. To generate diverse complex terrains, we applied the
Neat-python API to evolve the CPPNs. The list below shows some primary
packages and corresponding versions that are used in our implementation. The
source code and trained models can be found at https://github.com/ml-
tue/ePOET$\\_$3D.git.

  * •

Python 3.8

  * •

Fiber 0.2.1

  * •

Pytorch 1.8.1

  * •

MuJoCo 200

  * •

mujoco_py 2.0.2.13

  * •

Gym 0.18.0

  * •

PyOpenGL 3.1.5

  * •

Neat-python 0.92

###  A.4 The agent

Our chosen agent is a Hexapod walker created by Azayev et al [1] as it has 18
Degrees-of-Freedom (DoF), which allows it to be more flexible on uneven
terrains and to have a suitable training complexity. As can be seen in Figure
7 (left), it has 6 legs and each actuator is defined with a feedback gain of
40. The highest and the lowest torso height is shown in Figure 7 (right-top)
and (right-bottom) respectively. The observation space and the action space of
this agent is 53 and 18, respectively. Since we did not add any sensors to the
hexapod, the observations of each time-step have to be informative enough for
the agent to perform the task. To this sense, the observations consist of the
torso positions, joint velocities, contacts information, as well as some
sampled heights of the agent’s current position. Borrowing the formula from
Azayev et al. [1], observations at time-step $t$ are defined as:
$o_{t}=\left\\{j_{t}^{1},..,j_{t}^{n},j_{t}^{1},..,j_{t}^{in},c_{t}^{1},..,c_{t}^{m},q_{t}^{w},q_{t}^{x},q_{t}^{y},q_{t}^{z}\right\\}$,
in which $j^{i}_{t}$ is the $i_{th}$ joint angle, $j_{t}^{n}$ is the
$i_{\text{th }}$ joint velocity, $c_{t}^{k}$ is a binary value representing
the $k_{th}$ contact of leg $k$ with the ground, and
$q_{t}^{w},q_{t}^{x},q_{t}^{y},q_{t}^{z}$ are the quaternion of rotation of
the torso.

Figure 7: A hexapod walker with 18 joints. The right-top image shows the
maximum body height that the hexapod can achieve, and the right-bottom image
shows its minimum body height.

The reward signal of the hexapod agent at time-step $t$ is determined as:

$R_{t}=w_{v}*R^{v}_{t}+w_{\theta}*R^{\theta}_{t}-R^{c}_{t}$ (1)

where $R^{v}_{t}$ represents the velocity reward that motivates the agent to
move forward, $R^{\theta}_{t}$ indicates the correcting heading error,
$R^{c}_{t}$ consists of a variety of costs, and $w_{i}$ indicates their
corresponding weights. The velocity reward $R^{v}_{t}$ in Equation 1 is
defined as:

$R^{v}_{t}=\left(\frac{1}{\left|x_{t}-v_{tar}\right|+1}-\frac{1}{v_{tar}+1}\right)*\frac{1}{\left(1+30*y_{t}^{2}\right)},$

where $x_{t},y_{t}$ represents the velocity in x- and y-direction,
respectively. The $v_{tar}$ represents the maximum target velocity to prevent
a very fast gait and erratic jumping behavior. It is important to find a
suitable $v_{tar}$ because a large $v_{tar}$ could cause erratic jumping
behaviors while a small $v_{tar}$ could result in slow movements. It is set to
0.4 in our experiments, which is a finetuned value by Azayev et al [1]. The
correcting heading error $R^{\theta}_{t}$ in Equation 1 is defined as:

$R^{\theta}_{t}=\left(\left|deviation_{\theta-1}\right|-\left|deviation_{\theta}\right|\right)+\left(\left|y\\_deviation_{\theta-1}\right|-\left|y\\_deviation_{\theta}\right|\right).$

This term is to motivate the agent to correct its head forward to the target
direction before starting to walk. According to Azayev et al. [1], they found
that correcting heading errors works significantly better than penalizing
heading deviations. The penalty term $R^{c}_{t}$ sums up the torso angle,
acceleration penalization, y-axis penalization, velocity in the z-axis, and
control costs (external forces of the body), with different weights. Finally,
the weights $w_{v}$ and $w_{\theta}$ is set to 6 and 10, respectively.

The reward shaping requires a lot of fine-tuning since it is not the main
focus of our study, we reused the hexapod agent created by Azayev et al [1]
and made two small changes to adapt it to our experiments instead of creating
an entirely new agent. One change is that we add a bonus reward if the agent
reaches the target position (the edge of the terrain along the positive
x-axis), and we assign it a $finished$ signal. This is to encourage the agent
to move towards the positive x-axis and the $finished$ signal is used to count
the number of terrains that the agent passed. The other change is that we add
an indicator to continuously track the agent’s direction along the x-axis. If
the agent continuously moves toward the negative direction along the x-axis
for more than 500 steps, then the current iteration is terminated. This is to
save training time.

###  A.5 Hyper-parameters

The settings shown in Table 4 are applied for both SAC experiments and ePOET-
SAC experiments. Note that the hidden layer shape of SAC differs from the
hidden layer shape of ES in the ePOET, which is $[256,256]$ as shown in Table
4 and $[40,40]$ as shown in Table 5, respectively. This is because the SAC
with the hidden layer shape of $[256,256]$ outperforms the hidden layer shape
of $[40,40]$, while the ES performs oppositely with these two hidden layer
shapes for the hexapod according to our experiments. Thus, we use a different
hidden layer size for SAC and ES, and randomly select the parameters when
performing crossover in our ePOET-SAC. We have also attempted to use the t-SNE
[24] dimension reduction technique to keep the distance relation while
reducing the shape from $[256,256]$ to $[40,40]$, but the resulted parameter
vectors had worse performance than the randomly selected parameter vectors.
Furthermore, to ensure that the SAC actor indeed helps improve the ES agents,
we pre-train the SAC for 500 iterations so that the pre-trained SAC actor
obtains a certain score before performing crossover with the ES agents.

It is also important to mention that we set the lower bound and upper bound of
PATA-EC (explained in section LABEL:subsec:existingapp) as 500 and 3000
respectively in our experiments, as can be seen in Table 5. These bounds are
to prevent that a newly generated terrains from being over-simple and over-
complex. Moreover, the threshold of reproducing new environments is set as
2000, which means that the agent needs to acquire a score of higher than 2000
for the current terrain in order to reproduce child terrains. According to our
experiments, a score of 2000 or higher indicates that the agent can somewhat
pass the current challenge. In the ePOET and ePOET-SAC experiments, every 75
iterations, a parent agent generates its child agent-environment pair if the
agent reaches a score of 2000. During every reproducing process, the parent
environment produces 8 child environments via CPPN-NEAT, but only one
environment is admitted and added to the pool with its paired agent. For the
experiments with RL methods (PPO, SAC, VMPO), since the algorithms do not have
CPPN-NEAT mechanism, their training terrains only contain the random bowl
shape. Furthermore, one training episode terminates when 2,000 environment
steps have elapsed, when the agent continuously heads to the opposite
direction for more than 500 steps, or when the agent arrives at the finish
line. An environment is considered solved when the agent reaches the finish
line and obtains a score of 2000 or above.

Table 3: Hyper-parameter settings of PPO Hyper-parameter | Value
---|---
Activation | Tanh
batch_size | 64
hidden_shape | [64,64]
trajectory_length | 2000
eval_episodes | 1
policy_net_learning_rate | 3e-4
value_net_learning_rate | 3e-4
tau | 0.95
entropy_coeff | 0.005
discount | 0.99
clip_para | 0.2
reward_scale | 1
opt_epochs | 10
replay_buffer_size | 1e-6

Table 4: Hyper-parameter settings of SAC Hyper-parameter | Value
---|---
Activation | ReLU
batch_size | 256
hidden_shape | [256,256]
trajectory_length | 2000
eval_episodes | 1
policy_net_learning_rate | 3e-4
value_net_learning_rate | 3e-4
tau | 0.95
discount | 0.99
pretrain_epoch | 1
reward_scale | 5
opt_epochs | 10
replay_buffer_size | 1e-6
target_hard_update_period | 1000
automatic_entropy_tuning | True
reparameterization | True
use_soft_update | True

Table 5: Hyper-parameter settings of ePOET and ePOET-SAC Hyper-parameter | Value
---|---
Activation | tanh
number of sample points for each ES step | 500
batch_size | 1
batches_per_chunk | 256
hidden_shape | [40,40]
trajectory_length | 2000
max_num_envs | 40
num_workers | 10
weight update method | Adam
initial learning rate | 0.01
lower bound of learning rate | 0.001
decay factor of learning rate per ES step | 0.9999
initial noise standard deviation for ES | 0.02
lower bound of noise standard deviation | 0.01
decay factor of noise standard deviation per ES step | 0.999
iterations_before_transfer | 15
adjust_interval | 5
mc_lower | 500
mc_upper | 3000
repro_threshold | 2000
Table 6: Hyper-parameters for instantiating and evolving CPPNs. Hyper-parameter | Value
---|---
initial connection | full_nodirect
activation default | random
activation_mutate_rate | 0.5
activation options | identity sin sigmoid square tanh gauss
aggregation default | sum
aggregation_mutate_rate | 0.1
aggregation_options | sum
bias init stdev | 0.1
bias init type | gaussian
bias max value | 10.0
bias min value | -10.0
bias mutate power | 0.1
bias mutate rate | 0.75
compatibility disjoint coefficient | 1.0
compatibility weight coefficient | 0.5
enabled default | True
feed forward | True
node add prob | 0.1
node delete prob | 0.075
num inputs | 2
num_hidden | 2
num outputs | 1
response init mean | 1.0
response init type | gaussian
response max value | 10.0
response min value | -10.0
response mutate power | 0.2
single structural mutation | True
structural mutation surer | default
weight init stdev | 0.25
weight init type | gaussian
weight max value | 10.0
weight min value | -10.0
weight mutate power | 0.1
weight mutate rate | 0.75

###  A.6 Behavior Analysis

To better understand the behavioral differences of the hexapod trained with
different approaches, we visualize the agent’s action representations with a
dimensionality reduction technique known as t-Distributed Stochastic Neighbor
Embedding (t-SNE) [24]. Figure 8 shows the visualizations of the agent’s
action representations while walking over three terrains (corresponding to
three different colors) with t-SNE. Each subplot corresponds to the hexapod
trained with different approaches. A more dispersed scatter plotting indicates
a larger distance among actions, thus, represents more diverse locomotion
behaviors. As can be seen in Figure 8, the hexapod trained with SAC (a) has
the most erratic behaviors while the hexapod trained with ePOET (c) has the
least behavioral diversity. This is not surprising as SAC encourages random
exploration while optimizing the rewards through a trade-off between
maximizing rewards and maximizing entropy. In comparison to ePOET, it seems
that ePOET-SAC (d) helped the agent to explore slightly more diverse
locomotion behaviors. A possible reason is that ePOET-SAC benefits from the
randomness of SAC.

(a) SAC

(b) PPO

(c) ePOET

(d) ePOET-SAC

Figure 8: Visualizations of the agent’s action representations while walking
over three terrains (corresponding to three different colors) with t-SNE. The
visualization from (a) to (e) corresponds to the hexapod trained with
different methods (SAC, PPO, ePOET, ePOET-SAC).

Through rendering the interactions between the agent and the environments,
which can be found in the video,333Video: https://youtu.be/i8QM1-FoAto we
found that the hexapod trained with ePOET-SAC shows the best balanced and
efficient walking styles on different terrains. Behavioral comparisons are
shown in Figure 9 (trained with ePOET-SAC) and Figure 10 (trained with ePOET).
The behaviors shown in these two figures are captured with the same terrain.
In Figure 9, the agent trained with ePOET-SAC shows well-balanced movements
and a wide range of joint angles. Particularly, the left-front leg is well
balanced with the right-back leg in order to stretch to a maximum horizontal
length. Moreover, flexibly making use of the inner joint of each leg with a
wide-angle range prevents it from slipping down when encountering a steep
slope, as shown with the orange markers in Figure 9. By contrast, although the
agent trained with ePOET also shows quite natural and well-balanced locomotion
behaviors, its weakness could be exposed when facing challenging tasks. As
shown in Figure 10, due to using a narrow joint angle of the left-front leg
when stretching the right-back leg, the stability (holding power) is weakened
so that it easily slips when facing a steep slope. Another reason that the
agent failed to overcome the slope is that it did not control the torso angle
in a suitable range that resulted in an excessive body height (straightly
standing pose). Therefore, it fell down from the slope. More performance
comparisons can be found in video3.

Figure 9: Behavioral visualization of the agent trained with ePOET-SAC. The
agent applies a wide angle range of each joint and has a good balance between
legs (left-front-leg is paired with right-back-leg, and right-front-leg is
paired with left-back-leg to reach the longest distance per step). Figure 10:
Behavioral visualization of the agent trained with ePOET. The angle range is
smaller than the above agent, which results in a shorter reach distance per
step and an unstable position when facing a steep slope.


# On harvesting physical predictions
from asymptotically safe quantum field theories

Frank Saueressig<EMAIL_ADDRESS>Agustín Silva
<EMAIL_ADDRESS>High Energy Physics Department, Institute for Mathematics,
Astrophysics, and Particle Physics, Radboud University, Nijmegen, The
Netherlands.

###### Abstract

Asymptotic safety is a powerful mechanism for obtaining a consistent and
predictive quantum field theory beyond the realm of perturbation theory. It
hinges on an interacting fixed point of the Wilsonian renormalization group
flow which controls the microscopic dynamics. Connecting the fixed point to
observations requires constructing the set of effective actions compatible
with this microscopic dynamics. Technically, this information is stored in the
UV-critical surface of the fixed point. In this work, we describe a novel
approach for extracting this information based on analytical and pseudo-
spectral methods. Our construction is illustrated at the level of the two-
dimensional Ising model and easily generalizes to any asymptotically safe
quantum field theory. It also constitutes an important step towards setting up
a well-founded swampland program within the gravitational asymptotic safety
program.

##  I Introduction

The Wilsonian renormalization group (RG) [1] provides a powerful tool for
understanding the impact of statistical or quantum fluctuations in a given
physical system [2, 3, 4, 5]. It organizes fluctuations in momentum shells and
integrates out the modes, starting from the most energetic ones and
subsequently moving to lower energies. This procedure creates RG trajectories
which connect effective descriptions of the same system at different values of
the coarse-graining scale.

A major success of this approach is a comprehensive picture of critical
phenomena and universality which are readily explained in terms of RG fixed
points which control the theory’s infrared (IR) behavior. Typically, one has
to tune a small number of parameters so that a RG trajectory is dragged into
the fixed point in the limit where all fluctuations are integrated out. These
trajectories then “forget” their microscopic origin and physical quantities
like correlation functions are completely dictated by the properties of the IR
fixed point.

From a physics perspective, one may also be interested in situations where a
RG fixed point provides the microscopic description of the system. This
situation applies to asymptotically free quantum field theories like quantum
chromodynamics where the relevant fixed point is the free theory. Asymptotic
safety, first proposed in [6, 7], generalizes this construction to interacting
fixed points, also called non-Gaussian fixed points (NGFPs). In both cases,
the fixed point equips the construction with predictive power. By definition,
RG trajectories whose ultraviolet (UV) completion is provided by the fixed
point span its UV-critical surface $\mathcal{S}_{\rm UV}$. In the vicinity of
the fixed point, this surface can be found by linearizing the RG flow and
identifying the UV-attractive directions. Typically, $\mathcal{S}_{\rm UV}$ is
finite-dimensional and embedded in a larger space called the theory space
$\mathcal{T}$. This allows to predict values for renormalized couplings based
on the free parameters specifying a RG trajectory within $\mathcal{S}_{\rm
UV}$. These predictions are obtained from the end-point of the RG trajectory
where all fluctuations have been accounted for. In the case of an IR fixed
point, this limit is controlled by the fixed point itself, so that its
properties are directly related to observables. For an UV fixed point, this
link is highly non-trival though and requires the construction of
$\mathcal{S}_{\rm UV}$. Depending on the dimensions of $\mathcal{S}_{\rm UV}$
and $\mathcal{T}$ this becomes cumbersome rather quickly. A prototypical
example is provided by asymptotically safe quantum gravity supplemented by the
matter degrees of freedom of the standard model, where dim($\mathcal{S}_{\rm
UV}$) is expected to be of order twenty [8]. This surface should then be
embedded in an even larger space also comprising couplings associated with
beyond the standard model physics for making predictions. It is clear that the
canonical approach of mapping out $\mathcal{S}_{\rm UV}$ based on a shooting-
method is intractable in such cases. Nevertheless, a solid knowledge about the
effective actions within $\mathcal{S}_{\rm UV}$ is crucial for making
predictions and potentially falsifying the asymptotic safety hypothesis.

The goal of this letter is to introduce novel strategies which allow to map
out $\mathcal{S}_{\rm UV}$ in a computationally efficient way using pseudo-
spectral methods [9, 10] which are readily boosted by machine-learning
algorithms to improve convergence. We illustrate the generic algorithm based
on the RG flow of a scalar field theory in less than four dimensions which has
already been extensively studied, e.g., in [11, 12, 13, 14, 15, 16, 17, 18].
Owed to the flexibility of the RG [3], our algorithm is applicable to a much
wider range of settings, including the study of phase transitions in
statistical physics [19], estimates of the rate of spontaneous nucleation
[20], strongly interacting and high temperature theories [21], up to quantum
gravity building on the gravitational asymptotic safety program [22, 23].

##  II RG flows, fixed points
and free parameters

A primary tool for computing Wilsonian RG flows is the Wetterich equation [24,
25, 26, 27],111Our construction is readily adapted to other RG equations like
the Polchinski equation [28]. Essentially, it applies to any situation where a
vector field is used to generate a constraint surface also including the case
of IR-repulsive surfaces of IR-fixed points.

$k\partial_{k}\Gamma_{k}=\frac{1}{2}\mathbf{Tr}\left[\left(\Gamma_{k}^{(2)}+\mathcal{R}_{k}\right)^{-1}k\partial_{k}\mathcal{R}_{k}\right],$
(1)

which encodes the dependence of the effective average action $\Gamma_{k}$ on
the coarse-graining scale $k$, see [29, 30, 31] for more details. The
effective average action lives on theory space $\mathcal{T}$. By definition,
$\mathcal{T}$ contains all action functionals which can be constructed from
the field content of the theory and are compatible with its symmetry
requirements. For the purpose of this work, we are interested in complete
(approximate) solutions of (1), $k\mapsto\Gamma_{k}$, which interpolate
between a NGFP in the limit $k\rightarrow\infty$ (asymptotic safety) and the
effective action $\Gamma\equiv\lim_{k\rightarrow 0}\Gamma_{k}$ associated with
observables.

Introducing a basis $\\{\mathcal{O}_{n}\\}$ on $\mathcal{T}$, $\Gamma_{k}$ can
be expanded as

$\Gamma_{k}=\sum_{n}\;\tilde{g}^{n}_{k}\,\mathcal{O}_{n}\,.$ (2)

Here the $\tilde{g}^{n}_{k}$ are the dimensionful couplings of the theory
which depend on $k$. Denoting the mass-dimension of $\tilde{g}^{n}_{k}$ by
$d_{n}$, these couplings are conveniently traded for their dimensionless
counterparts

$g^{n}_{k}\equiv\tilde{g}^{n}_{k}\,k^{d_{n}}\,.$ (3)

In practice, the expansion (2) retains a finite number of operators, $N$,
only. The couplings $g^{n}$, $n=1,\cdots,N$ can then be read as coordinates on
$\mathbb{R}^{N}$. Substituting (2) into (1), the beta functions capturing the
$k$-dependence of these couplings $g_{k}^{n}$ can be read off as the
coefficients multiplying the basis elements $\\{\mathcal{O}_{n}\\}$:

$k\partial_{k}g_{k}^{n}=\beta^{n}(g^{m})\,.$ (4)

In the sequel, we will assume that the beta functions have been computed and
admit a NGFP suitable for asymptotic safety.

By definition, fixed points $\\{g^{n}_{*}\\}$ satisfy
$\beta^{n}(g^{m}_{*})=0,\forall n$. The dimension of $\mathcal{S}_{\rm UV}$
associated with the fixed point is readily inferred by linearizing the beta
functions

$k\partial_{k}g^{n}_{k}\simeq\sum_{m}B^{n}{}_{m}(g^{m}_{k}-g^{m}_{*})\,,\quad
B^{n}{}_{m}\equiv\left.\frac{\partial\beta^{n}}{\partial
g^{m}}\right|_{g=g_{*}}\,.$ (5)

These equations are readily solved in terms of the stability coefficients
$\theta_{I}$ and right-eigenvectors $V_{I}$, satisfying
$BV_{I}=-\theta_{I}V_{I}$. Separating the stability coefficients
$\\{\theta_{I}\\}\rightarrow\\{\theta_{\alpha},\tilde{\theta}_{\mu}\\}$ where
$\alpha=1,\cdots,\dim(\mathcal{S}_{\rm UV})$ and
$\mu=1,\cdots,N-\dim(\mathcal{S}_{\rm UV})$ enumerate the stability
coefficients with Re$(\theta_{\alpha})>0$ and Re$(\tilde{\theta}_{\mu})<0$,
respectively, one has

$g^{n}_{k}\approx
g^{n}_{*}+\sum_{\alpha}C_{\alpha}\,V_{\alpha}^{n}\left(\frac{k_{0}}{k}\right)^{\theta_{\alpha}}+\sum_{\mu}\tilde{C}_{\mu}\,\tilde{V}_{\mu}^{n}\left(\frac{k_{0}}{k}\right)^{\tilde{\theta}_{\mu}}\,.$
(6)

Here $k_{0}$ is a reference scale and $C_{\alpha},\tilde{C}_{\mu}$ are
integration constants labelling the specific solutions. The eigendirections
$V_{\alpha}$ $(\tilde{V}_{\mu})$ are attractive (repulsive) in the limit
$k\rightarrow\infty$. Thus asymptotic safety fixes $\tilde{C}_{\mu}=0$ while
the ${C}_{\alpha}$ are free parameters. RG trajectories where
$\tilde{C}_{\mu}=0$ span $\mathcal{S}_{\rm UV}$.

Subsequently, we reformulate the conditions on the integration constants in
terms of the dimensionless couplings. For this purpose, we split

$g^{n}\mapsto\\{u^{\alpha}\,,\,v^{\mu}\\}\,,\quad\beta^{n}\mapsto\\{\beta^{\alpha}\,,\,\beta^{\mu}\\}\,.$
(7)

The basic idea illustrated in Fig. 1 is to use the $u^{\alpha}$ as coordinates
on $\mathcal{S}_{\rm UV}$.222Depending on the structure of $\mathcal{S}_{\rm
UV}$, the initial coordinate system may not cover the entire surface. In this
case, $\mathcal{S}_{\rm UV}$ has to be covered with multiple coordinate
patches. A simple illustration is provided by the surface $x^{2}+y^{2}-1=0$
embedded into $\mathbb{R}^{2}$ where there is no single-valued function $y(x)$
that covers the entire surface. In the context of the RG such coordinate
changes have recently been discussed in [32]. Its embedding into
$\mathbb{R}^{N}$ is then given by a set of generating functions

$F_{\nu}(u,v)=0\,,\qquad\nu=1,\cdots,N-\dim(\mathcal{S}_{UV})\,.$ (8)

The generating functions are combinations of couplings whose values are
conserved along the RG flow. We adopt the convention that $\mathcal{S}_{\rm
UV}$ is generated by setting the conserved quantity to zero. Solving these
relations for $v^{\mu}$ then determines the values of the couplings $v^{\mu}$
in terms of the coordinates on the surface

$\left.v^{\mu}\right|_{S_{\rm UV}}=v^{\mu}(u^{\alpha})\,.$ (9)

Figure 1: Using the choice of couplings introduced below eq. (12), the $k\to
0$ limit of the couplings $u^{\alpha}$ is $\\{u^{\alpha}\\}=0$. Every complete
RG trajectory emanating from the NGFP gives rise to an effective action whose
parameters $v^{\mu}$ are read off from the intersect of $\mathcal{S}_{\rm UV}$
with the axes $\\{u^{\alpha}\\}=0$.

In the vicinity of the NGFP, the functions $F_{\nu}$ can be found as follows.
The vectors $\\{V_{\alpha},\tilde{V}_{\mu}\\}$ span the tangent space at
$\\{g_{*}^{n}\\}$. The vectors $V_{\alpha}$ are tangent to $\mathcal{S}_{\rm
UV}$ while, in general, the $\tilde{V}_{\mu}$ are not normal to the surface. A
set of normal vectors can then be constructed in two steps. First, we project
the $V_{\alpha}$ into a orthogonal basis $t_{\alpha}$, carrying out a Gram-
Schmidt process. Starting from the vectors $\tilde{V}_{\mu}$, the set of
normals is obtained as

$n_{\mu}=\tilde{V}_{\mu}-\sum_{\alpha}\,t_{\alpha}\frac{\langle
t_{\alpha}|\tilde{V}_{\mu}\rangle}{\langle t_{\alpha}|t_{\alpha}\rangle}\,,$
(10)

where $\langle x|y\rangle$ is the standard vector product on $\mathbb{R}^{N}$.
At the linearized level, the constrained equations (8) are then given by

$F_{\mu}^{\ast}(g^{n})=\sum_{n}\,n^{n}_{\mu}\,\left(g^{n}-g^{n}_{*}\right)\,.$
(11)

Here $\ast$ emphasizes that this the expansion of the full $F_{\mu}$ to linear
order at the NGFP. One may then single out couplings $u^{\alpha}$ according to
(7) and solve these linear relations to express the $v^{\mu}$ in terms of
these parameters.

##  III Constructing $\mathcal{S}_{\rm UV}$

Finding the IR-endpoints of asymptotically safe RG trajectories requires the
construction of $\mathcal{S}_{\rm UV}$ beyond the linear approximation (11).
Technically, it is then convenient to work with a compact subregion of
$\mathbb{R}^{N}$ which contains the NGFP and the endpoints. In order to
achieve this, it is instructive to discuss the $k\rightarrow 0$-limit of eq.
(3). Keeping $\lim_{k\rightarrow 0}\tilde{g}^{n}_{k}$ fixed and finite, which
is a natural assumption for couplings appearing in the effective action, the
dimensionless couplings with non-zero mass-dimension $d_{n}$ either approach
zero or diverge

$\begin{array}[]{ll}\lim_{k\rightarrow
0}g^{n}_{k}=0\,,&d_{n}<0\,,\\\\[5.16663pt] \lim_{k\rightarrow
0}g^{n}_{k}=\pm\infty\,,&d_{n}>0\,.\end{array}$ (12)

Couplings with $d_{n}=0$ take the finite value of the coupling appearing in
the effective action. Based on these insights, we choose the couplings
$\tilde{u}^{\alpha}$ to have a negative mass dimension. This can always be
achieved by making the identification $\tilde{u}^{n}=1/\tilde{g}^{n}$ if the
initial coupling comes with a positive mass dimension. This has the advantage
that the IR-endpoints are located at $\\{u^{\mu}\\}=0$. Furthermore, we can
take the $v^{\mu}$ to be dimensionless by constructing appropriate ratios of
the couplings, see (20) for an explicit example. The physical predictions for
the couplings $v^{\mu}$ are then read off at $\\{u^{\alpha}\\}=0$.

###  III.1 Master equation and explicit form

We now derive our master equation encoding the structure of $\mathcal{S}_{\rm
UV}$. Taking the $k$-derivative of the generating functions $F_{\mu}(u,v)$ and
using the definition of the beta functions (4) leads to

$\sum_{\alpha}\,\frac{\partial F_{\mu}}{\partial
u^{\alpha}}\,\beta^{\alpha}+\sum_{\nu}\frac{\partial F_{\mu}}{\partial
v^{\nu}}\,\beta^{\nu}=0\,.$ (13)

At this point we assume that $\det(\frac{\partial F_{\mu}}{\partial
v^{\nu}})\not=0$, so that the implicit function theorem guarantees that the
system $F_{\mu}=0$ has a (local) solution in the form (9).

The master equation can be turned into an system of non-linear partial
differential equations for the functions (9). Substituting
$F_{\mu}(u,v)=v^{\mu}-v^{\mu}(u)|_{\mathcal{S}_{\rm UV}}$, which obviously
satisfies $F_{\mu}|_{\mathcal{S}_{\rm UV}}=0$, leads to the linear partial
differential equation

$\sum_{\alpha}\,\beta^{\alpha}\,\frac{\partial v^{\mu}}{\partial
u^{\alpha}}=\beta^{\mu}\,.$ (14)

The boundary conditions for the system are given in the linearized regime, eq.
(11). Solutions can then be obtained either by solving the system (13) or (14)
using numerical methods or extending the power series (11) to include higher-
order terms. The latter strategy requires checking whether the endpoint of
interest is within the radius of convergence of the series.

###  III.2 Implicit solutions

The key strength of the master equation (13) is that it is a multi-linear,
first order system of partial differential equations for the generating
functions $F_{\mu}(u,v)$. Finding the implicit form of $\mathcal{S}_{\rm UV}$
by solving this system leads to a significant reduction in complexity as
compared to solving the explicit system (14). In particular, the linear nature
of the master equation makes the use of pseudo-spectral methods and
collocation techniques [33, 34] a highly efficient tool for obtaining
solutions.333For earlier applications of pseudo-spectral methods in the
context of the functional renormalization group see [35, 36, 37, 38].

The basic idea is to select a dense set of functions $\\{\psi_{i}(u,v)\\}$.
Typical examples used in practice include Chebyshev polynomials [9] or radial
basis functions [39]. The generating functions are then approximated by
expanding in this basis

$F_{\mu}(u,v)\simeq\sum_{i=1}^{N_{p}}p_{\mu,i}\,\psi_{i}(u,v)\,,$ (15)

i.e., the functions are determined by the $N_{p}\times(N-\dim(\mathcal{S}_{\rm
UV}))$ free parameters $p_{\mu,i}$. Since (13) contains derivatives of
$F_{\mu}(u,v)$ only, it determines solutions up to an additive constant only.
Following (11), we fix this freedom by demanding that

$F_{\mu}(u_{*},v_{*})=0\,.$ (16)

Secondly, the first derivatives of $F_{\mu}$ evaluated at the fixed point are
taken to agree with the linearized approximation

$\left.\frac{\partial F_{\mu}}{\partial
u^{\alpha}}\right|_{g^{n}_{*}}=\frac{\partial F^{\ast}_{\mu}}{\partial
u^{\alpha}}\,,\qquad\left.\frac{\partial F_{\mu}}{\partial
v^{\nu}}\right|_{g^{n}_{*}}=\frac{\partial F^{\ast}_{\mu}}{\partial
v^{\nu}}\,.$ (17)

This enforces that some of the expansion coefficients in (15) must take non-
vanishing values thereby eliminating the trivial solution.

Subsequently, one chooses $N_{p}$ collocation points. It is natural to take
the NGFP as one of these points and construct a grid covering the parameter
space of interest. Substituting the expansion (15) into (13) and evaluating
the resulting expressions at the collocation points then leads to a system of
algebraic equations which determines the coefficients $p_{\mu,i}$. The
predictions $v^{\mu}$ compatible with asymptotic safety are then found by
evaluating the approximate solution at $\\{u^{\alpha}\\}=0$.

_Systematic improvements._ In the simplest case, one may arrange the
collocation points $(u^{\alpha}_{c,i},v^{\mu}_{c,i})$ in a regular lattice and
expand the generating functions in terms of Multivariate Cauchy Distributions
(MCDs)

$\psi_{i}(u,v)=\left(1+\sum_{\alpha}\frac{(u^{\alpha}-u^{\alpha}_{c,i})^{2}}{\sigma^{2}}+\sum_{\mu}\frac{(v^{\mu}-v^{\mu}_{c,i})^{2}}{\sigma^{2}}\right)^{-1}\,.$
(18)

Here $\sigma$ is a smoothness parameter whose typical value is chosen to be of
the order of the average separation between the collocation points. Since the
substitution of (15) into the system (13) results in a linear equation, taking
one collocation point for each MCD leads to a unique solution for the
parameters $p_{\mu,i}$.

Depending on the concrete application, this “basic” setup may be improved by
working with less parameters than collocation points, and/or with base
functions more suitable for describing sharp edges, periodicity, or optimized
convergence [40]. Moreover, the positioning of the collocation points may be
optimized near interesting regions. In the case of using less parameters than
collocation points, the construction is overdetermined and finding solutions
turns into an optimization problem, where one could use machine learning [41]
techniques, such as gradient descent methods [42]. Since an exhaustive survey
of all these possibilities is beyond the scope of this letter, we just limit
the discussion to the simplest setting, which may then be seen as a starting
point for improvements.

##  IV Example: scalar field in two dimensions

We illustrate the working of our general algorithm based on a scalar field
theory in a two-dimensional Euclidean spacetime. The RG flow possesses a NGFP,
corresponding to the Ising universality class. This fixed point is already
seen in simple approximations and comes with a single UV-attractive
eigendirection. For our purpose, it is then sufficient to work with the local
potential approximation (LPA) [24]. In this case, the effective average action
is approximated by

$\Gamma_{k}[\phi]\simeq\int
d^{2}x\left[\frac{1}{2}(\partial\phi)^{2}+U_{k}[\phi]\right]\,,$ (19)

where $\phi(x)$ is a real scalar field with mass-dimension zero. Implementing
our convention on the couplings, we expand the potential according to

$U_{k}[\phi]=\frac{1}{4!\,\tilde{u}_{k}}\phi^{4}+\sum_{0<n\,,\,n\neq
2}^{N_{\rm max}}\frac{v^{n}_{k}}{(2n)!\,\tilde{u}_{k}}\,(\phi^{2})^{n}\,.$
(20)

This singles out the coupling $\tilde{u}_{k}$ (taken with a negative mass
dimension) associated with the $\phi^{4}$-term as the free parameter allowed
by asymptotic safety. The IR-value of the dimensionless couplings
$v^{n}=(v^{1},v^{3},v^{4},\cdots)$ will be predicted based on the UV-critical
surface of the NGFP.

$\Big{.}N_{\rm max}$ | $u_{*}$ | $v^{1}_{*}$ | $v^{3}_{*}$ | $v^{4}_{*}$ | $v^{5}_{*}$ | $\theta_{1}$ | $\tilde{\theta}_{1}$ | $\tilde{\theta}_{2}$ | $\tilde{\theta}_{3}$ | $\tilde{\theta}_{4}$ | $v^{1}$ | $v^{3}$ | $v^{4}$ | $v^{5}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\Big{.}4$ | $0.3790$ | $-0.0541$ |  |  |  | $1.852$ | $-2.519$ |  |  |  | $0.15$ |  |  |
$\Big{.}6$ | $0.2765$ | $-0.0746$ | $16.3346$ |  |  | $1.699$ | $-2.875$ | $-25.381$ |  |  | $0.08$ | $31$ |  |
$\Big{.}8$ | $0.2695$ | $-0.0961$ | $24.1984$ | $941.6$ |  | $1.573$ | $-3.043$ | $-22.469$ | $-78.936$ |  | $0.07$ | $46$ | $3209$ |
$\Big{.}10$ | $0.2809$ | $-0.1174$ | $28.1641$ | $1561.8$ | $110325$ | $1.468$ | $-3.134$ | $-21.111$ | $-66.552$ | $-173.426$ | $0.066$ | $55$ | $5538$ | $701866$
Table 1: Properties of the NGFP associated with the Ising universality class
obtained from the LPA, eq. (20), up to order $N_{\rm max}=10$. The first and
second block give the position and stability coefficients of the fixed point,
respectively. The values $v^{\mu}\equiv\lim_{k\rightarrow 0}v_{k}^{\mu}$
predicted by asymptotic safety are summarized in the third block.

In Table 1 we report results up to order $N_{\rm max}=10$.444These results
serve as an illustration of the method and should not be mistaken as an
attempt of a precision computation. The latter may be achieved by extending
the approximation along the lines of the derivative expansion used for the
Wilson-Fisher fixed point [18, 43]. For simplicity, we limit our exposition to
the simplest case where $N_{\rm max}=2$, i.e., we are dealing with two
couplings $(u,v)$ and refer to the supplementary notebook for the detailed
implementation of this case.

We start from the beta functions given in [23]. Adapting to our
parameterization of the potential and setting $N_{\rm max}=2$ and $d=2$, we
have

$\beta_{u}=2u-\frac{3u{}^{3}}{2\pi\left(v+u\right){}^{3}},\quad\beta_{v}=-\frac{u{}^{2}\left(7v+u\right)}{4\pi\left(v+u\right){}^{3}}.$
(21)

The NGFP of the system is located at
$(u_{*},v_{*})=(\frac{343}{288\pi},-\frac{49}{288\pi})$. The stability
coefficients and right-eigenvectors entering into the linearized solution (6)
are

$\begin{split}\theta_{1}=\frac{1}{3}(1+\sqrt{43})\,,\quad&V_{1}=\left(\frac{1}{21}(8+\sqrt{43})\,,\,1\right)^{\rm
T}\,,\\\
\tilde{\theta}_{1}=\frac{1}{3}(1-\sqrt{43})\,,\quad&\tilde{V}_{1}=\left(\frac{1}{21}(8-\sqrt{43})\,,\,1\right)^{\rm
T}\,.\end{split}$ (22)

Thus we are dealing with a saddle-point with one UV-attractive and one UV-
repulsive eigendirection. The UV-critical surface is one-dimensional and we
expect a prediction for $v$ from asymptotic safety. Following the procedure
leading to (11), the generating function for the UV-critical surface is

$\begin{split}F^{\ast}(u,v)=&\left(172-137\sqrt{43}\right)u+\left(44\sqrt{43}+215\right)v\\\
&+\frac{49\left(1003\sqrt{43}-989\right)}{288\pi}\,.\end{split}$ (23)

Based on this input we compute the prediction for $v$ using 1) the traditional
shooting method, 2) a solution of the explicit system (14) based on series
expansions of orders 10 to 15, and 3) a spectral solution evaluating the
implicit equation (13) for a linear combination of 500 MCDs on a regular
lattice centered on the NGFP. All methods yield

$v\equiv\lim_{k\rightarrow 0}v_{k}\simeq 0.15\,,$ (24)

with the numerical difference appearing in the third relevant digit. Fig. 2
illustrates the structures underlying this prediction.

Figure 2:  Phase diagram generated by the beta functions (21) with arrows
pointing towards lower values of $k$. The NGFP and the endpoint of
$\mathcal{S}_{\rm UV}$ giving rise to the prediction (24) are marked by $*$.
The UV-critical surface obtained from the shooting method is marked by hollow
green dots, the solution of the implicit method by brown dashes, and the one
of the explicit method by a solid magenta line. For completeness we added the
RG trajectories generated by $\tilde{V}_{1}$ as solid orange dots.

##  V Summary and Outlook

The Wilsonian renormalization group (RG) is a powerful tool to study the
effect of fluctuations in statistical and quantum systems. Typically, such
investigations proceed along the following steps. 1) one computes the beta
functions of the theory in a suitable approximation. On this basis, one
determines the RG fixed points and their stability properties. 2) In the case
of asymptotic safety, where the UV-completion is provided by an interacting
fixed point, one constructs the UV-critical surface of the fixed point in
order to obtain the set of effective actions compatible with this UV-
completion. 3) The output of step 2) is used to construct physical observables
which allow to confront the predictions with observations.

Our work proposes an efficient algorithm for carrying out the second step in
this procedure, by solving multi-linear partial differential equations for the
generating functions encoding the embedding of the UV-critical surface. This
constitutes a significant improvement compared to the usual shooting method
(employed e.g., in [44, 45, 46, 47, 48]) which tracing individual RG
trajectories. We illustrated the working of our algorithm based on the
interacting fixed point corresponding to the Ising universality class.
Interpreting this fixed point as the UV-completion of the RG flow, we
predicted the couplings of the effective potential based on the single free
parameter allowed by asymptotic safety (see Table 1).

We stress that the methods described in this letter are applicable in a much
broader context. Potential applications beyond the scope of the present work
are the gravitational asymptotic safety program [22, 23], including its
extension by matter degrees of freedom [49, 8], asymptotically safe gauge
theories [50] and their supersymmetric extensions [51], and asymptotically
safe quantum electrodynamics [52]. The authors have already tested its
applicability for asymptotically safe scalar-tensor theories [53, 54, 55, 56,
57, 58, 59, 60] and the resulting shapes of the effective potentials will be
reported in [61].

## Acknowledgements

We thank M. Becker, C. Laporte, R. Loll, and J. Wang for useful discussions
during the process of developing these ideas.

## References

  * Wilson and Kogut [1974] K. G. Wilson and J. B. Kogut, The Renormalization group and the epsilon expansion, Phys. Rept. 12, 75 (1974).
  * Berges  _et al._ [2002] J. Berges, N. Tetradis, and C. Wetterich, Non-perturbative renormalization flow in quantum field theory and statistical physics, Physics Reports 363, 223 (2002).
  * Dupuis  _et al._ [2021] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, The nonperturbative functional renormalization group and its applications, Phys. Rept. 910, 1 (2021), arXiv:2006.04853 [cond-mat.stat-mech] .
  * Reuter and Saueressig [2012] M. Reuter and F. Saueressig, Quantum Einstein Gravity, New J. Phys. 14, 055022 (2012), arXiv:1202.2274 [hep-th] .
  * Zinn-Justin [2021] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , Vol. 171 (Oxford university press, 2021).
  * Weinberg [1978] S. Weinberg, Critical phenomena for field theorists, in  _Understanding the fundamental constituents of matter_ (Springer, 1978) pp. 1–52.
  * Weinberg [1979] S. Weinberg, Ultraviolet divergences in quantum theories of gravitation, in  _General relativity_ (1979).
  * Pastor-Gutiérrez  _et al._ [2023] A. Pastor-Gutiérrez, J. M. Pawlowski, and M. Reichert, The Asymptotically Safe Standard Model: From quantum gravity to dynamical chiral symmetry breaking, SciPost Phys. 15, 105 (2023), arXiv:2207.09817 [hep-th] .
  * Boyd [2013] J. Boyd,  _Chebyshev and Fourier Spectral Methods: Second Revised Edition_, Dover Books on Mathematics (Dover Publications, 2013).
  * Shen  _et al._ [2011] J. Shen, T. Tang, and L.-L. Wang,  _Spectral methods: algorithms, analysis and applications_ , Vol. 41 (Springer Science & Business Media, 2011).
  * Tetradis and Wetterich [1994a] N. Tetradis and C. Wetterich, Critical exponents from effective average action, Nucl. Phys. B 422, 541 (1994a), arXiv:hep-ph/9308214 .
  * Morris [1994a] T. R. Morris, On truncations of the exact renormalization group, Phys. Lett. B 334, 355 (1994a), arXiv:hep-th/9405190 .
  * Adams  _et al._ [1995] J. A. Adams, J. Berges, S. Bornholdt, F. Freire, N. Tetradis, and C. Wetterich, Solving nonperturbative flow equations, Mod. Phys. Lett. A 10, 2367 (1995), arXiv:hep-th/9507093 .
  * Litim [2002] D. F. Litim, Critical exponents from optimized renormalization group flows, Nucl. Phys. B 631, 128 (2002), arXiv:hep-th/0203006 .
  * Canet  _et al._ [2003] L. Canet, B. Delamotte, D. Mouhanna, and J. Vidal, Nonperturbative renormalization group approach to the Ising model: A Derivative expansion at order partial**4, Phys. Rev. B 68, 064421 (2003), arXiv:hep-th/0302227 .
  * Jüttner  _et al._ [2017] A. Jüttner, D. F. Litim, and E. Marchais, Global Wilson–Fisher fixed points, Nucl. Phys. B 921, 769 (2017), arXiv:1701.05168 [hep-th] .
  * De Polsi  _et al._ [2020] G. De Polsi, I. Balog, M. Tissier, and N. Wschebor, Precision calculation of critical exponents in the $O(N)$ universality classes with the nonperturbative renormalization group, Phys. Rev. E 101, 042113 (2020), arXiv:2001.07525 [cond-mat.stat-mech] .
  * Balog  _et al._ [2019] I. Balog, H. Chaté, B. Delamotte, M. Marohnic, and N. Wschebor, Convergence of Nonperturbative Approximations to the Renormalization Group, Phys. Rev. Lett. 123, 240604 (2019), arXiv:1907.01829 [cond-mat.stat-mech] .
  * Berges and Wetterich [1997] J. Berges and C. Wetterich, Equation of state and coarse grained free energy for matrix models, Nuclear Physics B 487, 675 (1997).
  * Strumia  _et al._ [1999] A. Strumia, N. Tetradis, and C. Wetterich, The region of validity of homogeneous nucleation theory, Physics Letters B 467, 279 (1999).
  * Tetradis and Wetterich [1994b] N. Tetradis and C. Wetterich, High temperature phase transitions without infrared divergences, International Journal of Modern Physics A 9, 4029 (1994b).
  * Percacci [2017] R. Percacci,  _An introduction to covariant quantum gravity and asymptotic safety_ , Vol. 3 (World Scientific, 2017).
  * Reuter and Saueressig [2019] M. Reuter and F. Saueressig,  _Quantum gravity and the functional renormalization group: the road towards asymptotic safety_ (Cambridge University Press, 2019).
  * Wetterich [1993] C. Wetterich, Exact evolution equation for the effective potential, Physics Letters B 301, 90 (1993).
  * Morris [1994b] T. R. Morris, The exact renormalization group and approximate solutions, International Journal of Modern Physics A 9, 2411 (1994b).
  * Reuter and Wetterich [1994] M. Reuter and C. Wetterich, Effective average action for gauge theories and exact evolution equations, Nucl. Phys. B 417, 181 (1994).
  * Reuter [1998] M. Reuter, Nonperturbative evolution equation for quantum gravity, Phys. Rev. D 57, 971 (1998), arXiv:hep-th/9605030 .
  * Polchinski [1984] J. Polchinski, Renormalization and effective lagrangians, Nuclear Physics B 231, 269 (1984).
  * Pawlowski [2007] J. M. Pawlowski, Aspects of the functional renormalisation group, Annals Phys. 322, 2831 (2007), arXiv:hep-th/0512261 .
  * Rosten [2012] O. J. Rosten, Fundamentals of the Exact Renormalization Group, Phys. Rept. 511, 177 (2012), arXiv:1003.1366 [hep-th] .
  * Saueressig [2023] F. Saueressig, The Functional Renormalization Group in Quantum Gravity (2023) arXiv:2302.14152 [hep-th] .
  * Buccio and Percacci [2022] D. Buccio and R. Percacci, Renormalization group flows between Gaussian fixed points, JHEP 10, 113, arXiv:2207.10596 [hep-th] .
  * Fornberg and Sloan [1994] B. Fornberg and D. M. Sloan, A review of pseudospectral methods for solving partial differential equations, Acta numerica 3, 203 (1994).
  * Hussaini  _et al._ [1983] M. Y. Hussaini, C. L. Streett, and T. A. Zang,  _Spectral methods for partial differential equations_ , Tech. Rep. (1983).
  * Borchardt and Knorr [2016] J. Borchardt and B. Knorr, Solving functional flow equations with pseudo-spectral methods, Phys. Rev. D 94, 025027 (2016), arXiv:1603.06726 [hep-th] .
  * Borchardt and Knorr [2015] J. Borchardt and B. Knorr, Global solutions of functional fixed point equations via pseudospectral methods, Phys. Rev. D 91, 105011 (2015), [Erratum: Phys.Rev.D 93, 089904 (2016)], arXiv:1502.07511 [hep-th] .
  * Draper  _et al._ [2020] T. Draper, B. Knorr, C. Ripken, and F. Saueressig, Finite Quantum Gravity Amplitudes: No Strings Attached, Phys. Rev. Lett. 125, 181301 (2020), arXiv:2007.00733 [hep-th] .
  * Fehre  _et al._ [2023] J. Fehre, D. F. Litim, J. M. Pawlowski, and M. Reichert, Lorentzian Quantum Gravity and the Graviton Spectral Function, Phys. Rev. Lett. 130, 081501 (2023), arXiv:2111.13232 [hep-th] .
  * Buhmann [2000] M. D. Buhmann, Radial basis functions, Acta numerica 9, 1 (2000).
  * Pachón and Trefethen [2009] R. Pachón and L. N. Trefethen, Barycentric-remez algorithms for best polynomial approximation in the chebfun system, BIT Numerical Mathematics 49, 721 (2009).
  * Powell [1981] M. J. D. Powell,  _Approximation theory and methods_ (Cambridge university press, 1981).
  * Bottou [2010] L. Bottou, Large-scale machine learning with stochastic gradient descent, in  _Proceedings of COMPSTAT’2010: 19th International Conference on Computational StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers_ (Springer, 2010) pp. 177–186.
  * Delamotte  _et al._ [2024] B. Delamotte, G. De Polsi, M. Tissier, and N. Wschebor, Conformal invariance and composite operators: A strategy for improving the derivative expansion of the nonperturbative renormalization group, (2024), arXiv:2401.02517 [cond-mat.stat-mech] .
  * Morris [1997] T. R. Morris, Three-dimensional massive scalar field theory and the derivative expansion of the renormalization group, Nucl. Phys. B 495, 477 (1997), arXiv:hep-th/9612117 .
  * Bervillier  _et al._ [2008] C. Bervillier, B. Boisseau, and H. Giacomini, Analytical approximation schemes for solving exact renormalization group equations in the local potential approximation, Nucl. Phys. B 789, 525 (2008), arXiv:0706.0990 [hep-th] .
  * Reuter and Weyer [2004] M. Reuter and H. Weyer, Renormalization group improved gravitational actions: A Brans-Dicke approach, Phys. Rev. D 69, 104022 (2004), arXiv:hep-th/0311196 .
  * Shaposhnikov [2012] M. Shaposhnikov, Asymptotic safety of gravity and the Higgs-boson mass, Theor. Math. Phys. 170, 229 (2012).
  * Eichhorn and Held [2018] A. Eichhorn and A. Held, Top mass from asymptotic safety, Phys. Lett. B 777, 217 (2018), arXiv:1707.01107 [hep-th] .
  * Eichhorn and Schiffer [2022] A. Eichhorn and M. Schiffer, Asymptotic safety of gravity with matter, (2022), arXiv:2212.07456 [hep-th] .
  * Litim and Sannino [2014] D. F. Litim and F. Sannino, Asymptotic safety guaranteed, JHEP 12, 178, arXiv:1406.2337 [hep-th] .
  * Bond and Litim [2017] A. D. Bond and D. F. Litim, Asymptotic safety guaranteed in supersymmetry, Phys. Rev. Lett. 119, 211601 (2017), arXiv:1709.06953 [hep-th] .
  * Gies and Ziebell [2020] H. Gies and J. Ziebell, Asymptotically Safe QED, Eur. Phys. J. C 80, 607 (2020), arXiv:2005.07586 [hep-th] .
  * Narain and Percacci [2010] G. Narain and R. Percacci, Renormalization Group Flow in Scalar-Tensor Theories. I, Class. Quant. Grav. 27, 075001 (2010), arXiv:0911.0386 [hep-th] .
  * Henz  _et al._ [2013] T. Henz, J. M. Pawlowski, A. Rodigast, and C. Wetterich, Dilaton Quantum Gravity, Phys. Lett. B 727, 298 (2013), arXiv:1304.7743 [hep-th] .
  * Percacci and Vacca [2015] R. Percacci and G. P. Vacca, Search of scaling solutions in scalar-tensor gravity, Eur. Phys. J. C 75, 188 (2015), arXiv:1501.00888 [hep-th] .
  * Henz  _et al._ [2017] T. Henz, J. M. Pawlowski, and C. Wetterich, Scaling solutions for Dilaton Quantum Gravity, Phys. Lett. B 769, 105 (2017), arXiv:1605.01858 [hep-th] .
  * Pawlowski  _et al._ [2019] J. M. Pawlowski, M. Reichert, C. Wetterich, and M. Yamada, Higgs scalar potential in asymptotically safe quantum gravity, Phys. Rev. D 99, 086010 (2019), arXiv:1811.11706 [hep-th] .
  * Laporte  _et al._ [2021] C. Laporte, A. D. Pereira, F. Saueressig, and J. Wang, Scalar-tensor theories within Asymptotic Safety, JHEP 12, 001, arXiv:2110.09566 [hep-th] .
  * Laporte  _et al._ [2023] C. Laporte, N. Locht, A. D. Pereira, and F. Saueressig, Evidence for a novel shift-symmetric universality class from the functional renormalization group, Phys. Lett. B 838, 137666 (2023), arXiv:2207.06749 [hep-th] .
  * de Brito  _et al._ [2023] G. P. de Brito, B. Knorr, and M. Schiffer, On the weak-gravity bound for a shift-symmetric scalar field, Phys. Rev. D 108, 026004 (2023), arXiv:2302.10989 [hep-th] .
  * [61] A. Silva, Emergence of inflaton potential from asymptotically safe gravity, In preparation .