# Benchmarking Reinforcement Learning Techniques for Autonomous Navigation

Zifan Xu1, Bo Liu1, Xuesu Xiao2,3, Anirudh Nair1, and Peter Stone1,4
1Department of Computer Science, University of Texas at Austin 2Department of
Computer Science, George Mason University 3Everyday Robots 4Sony AI.

###### Abstract

Deep reinforcement learning (RL) has brought many successes for autonomous
robot navigation. However, there still exists important limitations that
prevent _real-world_ use of RL-based navigation systems. For example, most
learning approaches lack safety guarantees; and learned navigation systems may
not generalize well to unseen environments. Despite a variety of recent
learning techniques to tackle these challenges in general, a lack of an open-
source benchmark and reproducible learning methods specifically for autonomous
navigation makes it difficult for roboticists to choose what learning methods
to use for their mobile robots and for learning researchers to identify
current shortcomings of general learning methods for autonomous navigation. In
this paper, we identify four major desiderata of applying deep RL approaches
for autonomous navigation: (D1) reasoning under uncertainty, (D2) safety, (D3)
learning from limited trial-and-error data, and (D4) generalization to diverse
and novel environments. Then, we explore four major classes of learning
techniques with the purpose of achieving one or more of the four desiderata:
memory-based neural network architectures (D1), safe RL (D2), model-based RL
(D2, D3), and domain randomization (D4). By deploying these learning
techniques in a new open-source large-scale navigation benchmark and real-
world environments, we perform a comprehensive study aimed at establishing to
what extent can these techniques achieve these desiderata for RL-based
navigation systems.

##  I Introduction

Autonomous robot navigation, i.e., moving a robot from one point to another
without colliding with any obstacle, has been studied by the robotics
community for decades. Classical navigation systems [1, 2] can successfully
solve such navigation problem in many real-world scenarios, e.g., handling
noisy, partially observable sensory input but still providing verifiable
collision-free safety guarantees. However, these systems require extensive
engineering effort, and can still be brittle in challenging scenarios, e.g.,
in highly constrained environments. This is reflected by a recent competition
(The BARN Challenge [3]) held in ICRA 2022, which suggests that even
experienced roboticists tend to underestimate how difficult navigation
scenarios are for real robots. Recently, data-driven approaches have also been
used to tackle the navigation problem [4] thanks to advances in the machine
learning community. In particular, Reinforcement Learning (RL), i.e., learning
from self-supervised trial-and-error data, has achieved tremendous progress on
multiple fronts, including safety [5, 6, 7], generalizability [8, 9, 10, 11],
sample efficiency [12, 13], and addressing temporal data [14, 15, 16]. For the
problem of navigation, learned navigation systems from RL [17] have the
potential to relieve roboticists from extensive engineering efforts [18, 19,
20, 21, 22] spent on developing and fine-tuning classical systems. Moreover, a
simple case study conducted in five randomly generated obstacle courses where
classical navigation systems often fail shows that a RL-based navigation has
the potential to achieve superior behaviors in terms of successful collision
avoidance and goal reaching (Fig. 1 left).

Despite such promising advantages, learning-based navigation systems are far
from finding their way into real-world robotics use cases, which currently
still rely heavily on their classical counterparts. Such reluctance in
adopting learning-based systems in the real world stems from a series of
fundamental limitations of learning methods, e.g., lack of safety,
explainability, and generalizability. To make things even worse, a lack of
well-established comparison metrics and reproducible learning methods further
obfuscates the effects of different learning approaches on navigation across
both the robotics and learning communities, making it difficult to assess the
state of the art and therefore to adopt learned navigation systems in the real
world.

To facilitate research in developing RL-based navigation systems with the goal
of deploying them in real-world scenarios, we introduce a new open-source
large-scale navigation benchmark with a variety of challenging, highly
constrained obstacle courses to evaluate different learning approaches, along
with the implementation of several state-of-the-art RL algorithms. The
obstacle courses resemble highly-constraint real-world navigation environments
(Fig. 1 right), and present major challenges to existing classical navigation
systems, while RL-based navigation systems have the potential to perform well
in them (Fig. 1 left).

Figure 1: Left: Success rate of two classical navigation systems, DWA [2]
(red) and E-band [1] (blue), and vanilla end-to-end RL-based (green)
navigation systems (individually trained) in five randomly generated difficult
obstacle courses. The insets at the top show top-down views of the five
obstacle courses. Right: Navigation environments in the real world (left) and
the proposed benchmark (right) are similar to the robot perception system
(e.g., white/red laser scans and cyan/purple costmaps).

We identify four major desiderata that ought to be fulfilled by any learning-
based system that is to be deployed: (D1) reasoning under uncertainty of
partially observed sensory inputs, (D2) safety, (D3) learning from limited
trial-and-error data, and (D4) generalization to diverse and novel
environments. By deploying four major classes of learning techniques: memory-
based neural network architectures, safe RL, model-based RL, and domain
randomization, we perform extensive experiments and empirically compare a
large range of RL-based methods based on the degree to which they achieve each
of these desiderata. Each trial of these methods runs for about 8 hours and
uses 100 CPU cores on a computational cluster for distributed rollout of the
RL agents. The policies are updated on a local machine with one A40 GPU. The
whole set of experiments includes 120 independent trials which account for
about 100000 hours of CPU time and 800 hours of GPU time in total. To avoid
any requirement on having access to a high-performance compute cluster, the
provided code base collects the trajectories on one single local machine,
which uses about 16 CPU cores to finish one run in 48 hours. Moreover, by
deploying six selected navigation systems in three qualitatively different
real-world navigation environments, we investigate to what degree the
conclusions drawn from the benchmark can be applied to the real world.
Supplementary videos and material for this work are available on the project
webpage.111https://cs.gmu.edu/ xiao/Research/RLNavBenchmark/

##  II Desiderata for Learning-based Navigation

In this section, we introduce four desiderata for learning-based autonomous
navigation systems and briefly discuss the learning techniques as their
corresponding solutions.

(D1) reasoning under uncertainty of partially observed sensory inputs.
Autonomous navigation without explicit mapping and localization is usually
formalized as a Partially Observable Markov Decision Process (POMDP), where
the agent produces the motion of the robot only based on limited sensory
inputs that are usually not sufficient to recover the full state of the
navigation environment. Most RL approaches solve POMDPs by maintaining a
history of past observations and actions [14, 15]. Then, neural network
architectures like Recurrent Neural Networks (RNNs) that process sequential
data are employed to encode history and address partial observability. In this
study, we investigate various design choices of history-dependent
architectures.

(D2) safety. Even though in some cases deep RL methods achieve comparable
performance to classical navigation, they still suffer from poor
explainability and do not guarantee collision-free navigation. The lack of
safety guarantee is a major challenge preventing RL-based navigation from
being used in the real world. Prior works have addressed this challenge by
formalizing the navigation as a multi-objective problem that treats collision
avoidance as a separate objective from reaching the goal and solving it with
Lagrangian or Lyapunov-based methods [5]. For simplicity, we only explore
Lagrangian method and investigate whether explicitly treat safety as a
separate objective leads to safer and smoother learned navigation behavior.

(D3) learning from limited trial-and-error data. Although deep RL approaches
can alleviate roboticists from extensive engineering effort, a large amount of
data is still required to train a typical deep RL agent. However, autonomous
navigation data is usually expensive to collect in the real world. Therefore,
data collection is usually conducted in simulation, e.g., in the Robot
Operating System (ROS) Gazebo simulator, which provides an easy interface with
real-world robots. However, simulating a full navigation stack from perception
to actuation is more computationally expensive compared to other RL domains,
e.g., MuJuCo or Atari games [23, 24], which presents a high requirement for
sample efficiency. Most prior works have used off-policy RL algorithms to
improve sample efficiency with experience replay [25, 26]. In addition, model-
based RL methods can explicitly improve sample efficiency, and are widely used
in robot control problems. In this study, we compare two common classes of
model-based RL method [12, 13] combined with an off-policy RL algorithm, and
empirically study to what extent model-based approaches improve sample
efficiency when provided with different amounts of data.

(D4) generalization to diverse and novel environments. The ultimate goal of
deep RL approaches for autonomous navigation is to learn a generalizable
policy for all kinds of navigation environments in the real world. A simple
strategy is to train the agent in as many diverse navigation environments as
possible or domain randomization, but it is unclear what is the necessary
amount of training environments to efficiently achieve good generalization.
Utilizing the large-scale navigation benchmark proposed in this paper, we
empirically study the dependence of generalization on the number of training
environments.

##  III Navigation Benchmark

This section details the proposed navigation benchmark for RL-based navigation
systems, which aims to provide a unified and comprehensive testbed for future
autonomous navigation research. First, Sec. III-A discusses the difference
between the proposed benchmark and existing navigation benchmarks. In Sec.
III-B and III-C, the navigation task is formally defined and formulated as a
POMDP. More detailed background of MDP and POMDP can be found on the project
webpage. Finally, Sec. III-D introduces simulated and real-world environments
that benchmark different aspects of navigation performance.

###  III-A Existing Navigation Benchmarks

Our proposed benchmark differs from existing benchmarks in three aspects: (1)
high-fidelity physics: the navigation tasks are simulated by Gazebo [27],
which is based on realistic physical dynamics and therefore tests motion
planners that directly produce low-level motion commands, i.e., linear and
angular velocies, in contrast to high-level instructions such as turn left,
turn right, move forward [28, 29]. In other words, we focus on “how to
navigate” (motion planning), instead of “where to navigate” (path planning);
(2) ROS integration: our benchmark is based on ROS [30], which allows seamless
transfer of a navigation method developed and benchmarked in simulation
directly onto a physical robot with little (if any) effort; and (3) collision-
free navigation: the benchmark includes both static and dynamic environments,
and requires _collision-free_ navigation, whereas other benchmarks assume that
either collisions are possible [29] or collision-avoidance will be addressed
by other low-level controllers out of the scope of the benchmark [28]. A
special case is the photo-realistic interactive Gibson benchmark by Xia et.
al. [31], which intentionally allows physical interaction with objects (e.g.,
pushing) and therefore pose no challenges to the collision-avoidance system.

###  III-B Navigation Problem Definition

######  Definition 1 (Robot Navigation Problem).

Situated within a navigation environment $e$ which includes information of all
the obstacle locations at any time $t$, a start location $(x_{i},y_{i})$, a
start orientation $\theta_{i}$, and a goal location $(x_{g},y_{g})$, the
navigation problem $\mathcal{T}_{e}$ is to maximize the probability $p$ of a
mobile robot reaching the goal location from the start location and
orientation under a constraint on the number of collisions with any obstacle
$C<1$ and a time limit $t<T_{max}$.

A navigation problem can be formally defined as above. Given the current
location $(x_{t},y_{t})$, the robot is considered to have reached the goal
location if and only if its distance to the goal location is smaller than a
threshold, $d_{t}<d_{s}$, where $d_{t}$ is the Euclidean distance between
$(x_{t},y_{t})$ and $(x_{g},y_{g})$, and $d_{s}$ is a constant threshold.

###  III-C POMDP Formulation

A navigation task $\mathcal{T}_{e}$ can be formulated as a POMDP conditioned
on a navigation environment $e$, which can be represented by a 7-tuple
$(S_{e},A_{e},O_{e},T_{e},\gamma_{e},R_{e},Z_{e})$ . In this POMDP, the state
$s_{t}\in S_{e}$ is a 5-tuple $(x_{t},y_{t},\theta_{t},c_{t},e)$ with
$x_{t},y_{t},\theta_{t}$ the two-dimensional coordinates and the orientation
of the robot at time step $t$, $c_{t}$ a binary indicator of whether a
collision has occurred since the last time step $t-1$, and $e$ the navigation
environment. The action $a_{t}=(v_{t},\omega_{t})\in A_{e}$ is a two-
dimensional continuous vector that encodes the robot’s linear and angular
velocity. The observation $o_{t}=(\chi_{t},\bar{x_{t}},\bar{y_{t}})\in O_{e}$
is a 3-tuple composed of the sensory input $\chi_{t}$ from LiDAR scans and the
relative goal position $(\bar{x_{t}},\bar{y_{t}})$ in the robot frame. The
observation model $Z:S\rightarrow O$ maps the state to the observation. The
reward function for this POMDP is defined as follows:

$R_{e}(s_{t},a_{t})=+b_{f}\cdot\mathbb{1}(d_{t}<d_{s})+b_{p}\cdot(d_{t-1}-d_{t})-b_{c}\cdot
c_{t},$ (1)

where $\mathbb{1}(d_{t}<d_{s})$ is the indicator function of reaching the goal
location, $d_{t}$ is the Euclidean distance to the goal location, and $b_{f}$,
$b_{p}$, $b_{c}$ are the coefficient constants. In this reward function, the
first term is the true reward function that assigns a positive constant
$b_{f}$ for the success of an agent, which matches with the objective of the
navigation task in Definition 1. The second and third terms are auxiliary
rewards that facilitate the training by encouraging local progress and
penalizing collisions.

We perform a grid search over different values of the coefficients in this
reward function, and the result shows that the auxiliary reward term
$(d_{t-1}-d_{t})$ is necessary for successful training, and a much smaller
coefficient $b_{p}$ relative to $b_{f}$ can lead to a better asymptotic
performance. The agent can learn without the penalty reward for collision
($b_{c}=0$), but a moderate value of $b_{c}$ can improve the asymptotic
performance and speed up training. For all the experiments in this paper, we
fix the coefficients as $b_{f}=20$, $b_{p}=1$ and $b_{c}=4$.

In our experiments, the RL algorithm solves a multi-task RL problem where the
tasks are randomly sampled from a task distribution $\mathcal{T}_{e}\sim
p(\mathcal{T}_{e})$. Here the task distribution
$p(\mathcal{T}_{e}):=U(\\{e_{i}\\}_{i=1}^{N})$ is a uniform distribution on a
set of $N$ navigation environments $\\{e_{i}\\}_{i=1}^{N}$. The overall
objective of this multi-task RL problem is to find an optimal policy
$\pi^{*}=\max_{\pi}\mathbb{E}_{\mathcal{T}_{e}\sim
p(\mathcal{T}_{e}),\tau_{t}\sim\pi}\big{[}\sum_{t=0}^{\infty}\gamma^{t}R_{e}(s_{t},a_{t})\big{]}$.

###  III-D Navigation Environments

Figure 2: Three types of navigation environments: static (left), dynamic box
(middle), and dynamic-wall (right). The red squares mark the obstacle fields,
and the yellow arrows mark the direction of navigation. In dynamic-wall, the
green (blue) arrows indicate the case when the two walls are moving apart
(together). In dynamic box, the red arrows indicate the velocities of
obstacles.

The navigation is performed by a ClearPath Jackal differential-drive ground
robot in simulated by the Gazebo simulator. More details of the robot and
simulation can be found on the project webpage. Each environment in this
benchmark will have a navigation system navigating the robot through a 10m
navigation path that passes through a highly constrained obstacle course.
Walls are placed at three edges of a square so that passing through the
obstacle field is the only path to the goal location (see Fig. 2). The
benchmark includes 300 static environments, 100 dynamic-box environments, and
100 dynamic-wall environments. The static environments contains a diverse set
of obstacle course covering a large range of difficulty levels from easy to
hard. A dynamic-box environment has small boxes with random shapes and
velocities to test the system’s immediate reactions to small moving obstacles.
A dynamic-wall has two walls moving oppositely that requires the system to
make a longer-term decision of whether to pass or wait. The detailed
procedures of generating these environments can be found on the project
webpage. We randomly select 50 environments from each type as the test sets,
which are denoted as static-test, dynamic-box-test, and dynamic-wall-test. The
remaining environments are denoted as static-train, dynamic-box-train, and
dynamic-wall-train respectively. To study the effect of randomization, static-
train is further separated as static-train-5, static-train-10, static-
train-50, static-train-100, and static-train-250 by randomly sampling 5, 10,
50, 100, and all 250 environments from static-train.

To test the sim-to-real transferability of the policies learning with
different techniques, the navigation systems are deployed in three
qualitatively different static navigation environments including a benchmark-
like environment (Fig. 3 left), an indoor highly-constrained environment (Fig.
3 right), and a large-scale environment of 30 meters in length. We denote them
as real-world-1, real-world-2, and real-world-3 respectively.

Figure 3: Real-world benchmark-like (top-right), in-door highly-constrained
(top-left), and large-scale (bottom) environments. The yellow curves mark the
paths of navigation.

##  IV Experiments

In this section, we present experimental results of each studied technique to
achieve the proposed desiderata in Sec. II. We implement distributed training
pipelines (similar to [21]) of different RL algorithms including TD3[32],
SAC[33], and DDPG[34]. They perform similarly in the study of different neural
network architectures. For simplicity, all the experiments mentioned in this
section use TD3 combined with the corresponding techniques, and all the data
points are averaged over three independent runs.

Success rate (%) ($\uparrow$) | static env. | dynamic-box env. | dynamic-wall env.
---|---|---|---
$H=1$ | $H=4$ | $H=8$ | $H=1$ | $H=4$ | $H=8$ | $H=1$ | $H=4$ | $H=8$
MLP | $65\pm 4$ | $57\pm 7$ | $42\pm 2$ | $50\pm 5$ | $35\pm 2$ | $46\pm 3$ | $67\pm 7$ | $72\pm 1$ | $69\pm 4$
GRU | - | $51\pm 2$ | $43\pm 4$ | - | $48\pm 4$ | $45\pm 1$ | - | $\boldsymbol{82\pm 4}$ | $78\pm 5$
CNN | - | $55\pm 4$ | $45\pm 5$ | - | $42\pm 5$ | $40\pm 1$ | - | $63\pm 3$ | $43\pm 3$
Transformer | - | $\boldsymbol{68\pm 2}$ | $46\pm 3$ | - | $\boldsymbol{52\pm 1}$ | $44\pm 4$ | - | $33\pm 28$ | $15\pm 13$

TABLE I: (D1) Success rate (%) ($\uparrow$) of policies trained with different
neural network architectures and history lengths. $H$ is the history length of
the memory. Bold font indicates the best success rate for each type of
environment.

Methods | Baseline (model-free) | Lagrangian method | MPC (model-based) | DWA | TEB
---|---|---|---|---|---
Success rate (%) ($\uparrow$) | $65\pm 4$ | $74\pm 2$ | $70\pm 3$ | $\boldsymbol{82}$ | 70
Survival time (s) ($\uparrow$) | $8.0\pm 1.5$ | $16.2\pm 2.5$ | $55.7\pm 4.9$ | $\boldsymbol{62.7}$ | 26.9
Traversal time (s) ($\downarrow$) | $\boldsymbol{7.5\pm 0.3}$ | $8.6\pm 0.2$ | $24.7\pm 2.0$ | $35.6$ | 26.9

TABLE II: (D2) Success rate ($\uparrow$), survival time ($\uparrow$), and
traversal time ($\downarrow$) of policies trained with Lagrangian method, MPC
with probabilistic transition model, and DWA. The bold font indicates the best
number achieved for each type of metric.

###  IV-A Memory-based Neural Network Architectures (D1)

To benchmark the performance of different neural network (NN) architectures,
deep RL policies represented by architectures of Multilayer Perceptron (MLP),
One-dimensional Convolutional Neural Network (CNN), Gated Recurrent Units
(GRU), and Transformer with history length of 4 and 8 are trained in static-
train-50, and the two types of dynamic environments dynamic-box-train and
dynamic-wall-train from Sec. III-D. After training, the policies are tested in
their corresponding test sets. In addition, MLP with history length of one is
added as a memory-less baseline. Table II shows the success rates of policies
with different architectures and history lengths evaluated in static-test (
left), dynamic-wall-test (middle) and dynamic-box-test respectively.

Memory-based NNs only marginally improve navigation performance in static
environments. In Table II, the policy represented by Transformer with a
history length of 4 shows the best success rate of $68\%$, with a slightly
worse success rate of $65\%$ achieved by the baseline MLP. Additionally, a
monotonic decrease in success rate with increasing history length is observed
in each tested NN architecture. For example, a $32\%$ drop in the success rate
of Transformer is shown by increasing the history length from 4 to 8. One
possible explanation is that, if only few past observations are useful to make
the decision, including more history will make it more difficult to learn a
generalized policy in this very diverse training set.

Memory is essential when possible catastrophic failures will happen by making
the wrong long-term decisions. Memory usually matters for dynamic environments
when a single time frame is not sufficient to estimate the motion of
obstacles. Surprisingly, in dynamic-box where the dynamic obstacles are
completely random, the memory-based NN architectures do not outperform the
memory-less baseline. On the other hand, in dynamic-wall with a manually
designed dynamic challenge, the best success rate of $82\%$ is observed in GRU
with a history length of 4, which improves about $15\%$ over the non-memory
baseline. During our deployment of the policies, we observe that, in dynamic-
box even though the memory-less agent does not estimate the motion and adjust
its plan in advance, it tends to perform safely and avoids the obstacles when
they get close enough. This simple strategy works surprisingly well and
achieves similar success rate as the memory-based policies. However, this
strategy does not work in the manually designed dynamic challenges like
dynamic-wall where the agent has to estimate the motion of the obstacles to
pass safely.

###  IV-B Safe RL (D2)

To investigate to what extent safe RL methods can help to improve safety, a
TD3 agent with the Lagrangian-based safe RL method is trained in static-
train-50, and then tested in static-test. The policy is represented by a MLP
with its input containing only one history length. Table II shows the success
rate, average survival time, and average traversal time of the safe RL agent
trained with Lagrangian method and a baseline MLP agent tested in static-test.
We define survival time as the time cost of an unsuccessful episode (collision
or exceeding a time limit of 80s). Traversal time, instead, is the time cost
of a successful episode. With the same level of success rate, a longer
survival time means that the agent tends to, at least, avoid collisions if it
cannot succeed. To compare the safe RL method with classical navigation
systems which are believed to have better safety, we also add evaluation
metrics from a classical navigation stack with the Dynamic Window Approach
(DWA) [2] local planner.

Lagrangian method reduces the gap between training and test environments. When
deployed in the training environments, both the baseline MLP and the safe RL
method achieves about $80\%$ success rate. However, in the test environments,
the Lagrangian method has a better success rate of $74\%$ compare to $65\%$ by
the baseline MLP. We hypothesize that the safety constraint applied by the
safe RL methods forms a way of regularization, and therefore, improves the
generalization to unseen environments.

Lagrangian method increases the average survival time in failed episodes. As
expected, the Lagrangian method increases the average survival time by $8.2s$
compared to the baseline MLP at a cost of $1.1s$ longer average traversal
time. However, such improved safety are still worse than the classical
navigation systems given the best survival time of $88.6s$ achieved by DWA.

###  IV-C Model-based RL (D2 and D3)

To explore how the model-based approaches help with the autonomous navigation
tasks, we implement Dyna-style, MPC, and MBPO, and evaluate the methods in
static environments. The transition models are either represented by a
deterministic NN or a probabilistic NN that predicts the mean and variance of
the next state. During the training in static-train-50, the policies are saved
when 100k, 500k and 2000k transition samples are collected, then tested in
static-test. The success rates of these policies are reported in Table IV.

Model-based methods do not improve sample efficiency. As shown in the second
and third columns in Table IV, better success rates of $13\%$ and $58\%$ are
achieved by the baseline MLP method provided by limited 100k and 500k
transition samples respectively. In addition, Higher success rates at 500k
transition samples are observed in probabilistic models compared to their
deterministic counterparts, which indicates a more efficient learning with
probabilistic transition models. Notice that MBPO exploits more heavily on the
model compared to the Dyna-style method, which leads to much worse asymptotic
performance (about 20% success rate in the end).

Model-based methods with probabilistic dynamic models improve the asymptotic
performance. In the last column of Table IV, both Dyna-style and MPC with
probabilistic dynamic models achieve slightly better success rates of $70\%$
compared to $65\%$ in the baseline MLP method when sufficient transition
samples of 2000k are given to the learning agent.

The MPC policy performs conservatively when deployed in unseen test
environments and shows a better safety performance. The safety performances of
MPC policies with probabilistic dynamic models are also tested (see Table II).
We observe that the agents with MPC policies navigate very conservatively with
an average traversal time of $24.7s$, which is about two times more than the
MLP baseline. In the meantime, MPC policies achieve improved safety with the
best survival time of $55.7s$ among the RL-based methods.

Figure 4: (D4) Success rate (%) of policies trained with different number of
training environments.

| H | # envs | real-world-1 | real-world-2 | real-world-3
---|---|---|---|---|---
| traversal time ($\downarrow$) (# successful trials ($\uparrow$) / total #
trials)
MLP | 1 | 50 |  $6.9$ (1/3) |  $10.6$ (1/3) | N (0/3)
MLP | 1 | 250 | $4.6\pm 0.8~{}(3/3)$ | $\boldsymbol{6.6\pm 0.6~{}(3/3)}$ | $\boldsymbol{22.6\pm 0.5~{}(3/3)}$
Transformer | 4 | 50 |  $6.1\pm 0.4$ (3/3) |  $6.1\pm 0.1$ (2/3) |  $20.5\pm 2$ (2/3)
Lagrangian | 1 | 50 | $\boldsymbol{4.4\pm 0.6~{}(3/3)}$ |  $7.1\pm 0.1$ (2/3) |  $26.2$ (1/3)
MPC | 1 | 50 |  $13.2\pm 0.7$ (3/3) |  $24.8\pm 3.7$ (3/3) | N (0/3)
DWA | - | - |  $16.2\pm 0.7$ (3/3) |  $35.2\pm 8.2$ (2/3) |  $66.9\pm 0.6$ (3/3)

TABLE III: Physical experiments. The table shows the traversal time (s)
($\downarrow$) and the number of successful trials ($\uparrow$) of 5 RL-based
navigation systems and a classical navigation system (DWA) evaluated in three
real-world environments. The bold font indicates the best traversal time when
all three trials are successful.

###  IV-D Domain Randomization (D4)

To explore how model generalization depends on the degree of randomness in the
training environments, baseline MLP policies with one history length are
trained in the environment sets with 5, 10, 50, 100, and 250 training
environments. The trained policies are tested in the same static-test. To
investigate the performance gap between training and test, the policies
trained with 50, 100, and 250 environments are also tested on static-train-50,
which is part of their training sets. Fig 4 shows the success rate of policies
trained with different number of training environments.

The generalization to unseen environments improves with increasing number of
training environments. As shown in Fig. 4, the performances on the unseen test
environments monotonously increase from 43% to 74% with the number of training
environments increasing from 5 to 250. Moreover, the gaps between training and
test environments gradually shrink by adding more training environments
provided by that the polices are robust enough to maintain similar
performances of about 80% on the training environments.

Transition samples | 100k | 500k | 2000k
---|---|---|---
MLP | $\boldsymbol{13\pm 7}$ | $\boldsymbol{58\pm 2}$ | $65\pm 4$
Dyna-style deterministic | $8\pm 2$ | $30\pm 10$ | $66\pm 5$
MPC deterministic | $0\pm 0$ | $21\pm 10$ | $62\pm 3$
Dyna-style probabilistic | $0\pm 0$ | $48\pm 4$ | $\boldsymbol{70\pm 1}$
MPC probabilistic | $0\pm 0$ | $45\pm 4$ | $\boldsymbol{70\pm 3}$
MBPO | ${{\color[rgb]{0,0,0}0\pm 0}}$ | ${{\color[rgb]{0,0,0}0\pm 0}}$ | ${{\color[rgb]{0,0,0}21.9\pm 3}}$

TABLE IV: (D3) Success rate (%) ($\uparrow$) of policies trained with
different model-based methods and different number of transition samples. The
bold font indicates the best success rate for each number of transition
samples.

###  IV-E Physical experiments

To study the consistency of the above observations in simulation and the real
world, we deploy one baseline MLP policy, one best policy for each studied
desideratum, and one classical navigation system (DWA [2]) in the three real-
world environments introduced in Sec. III-D. Each deployment is repeated three
times, and the average traversal time and the number of successful trials are
reported in Table. III. Even though the best memory-based policy, transformer
architecture with 4 history length, was only marginally better than the
baseline MLP in simulation, in the real world it can navigate very smoothly
and fails only once in real-world-2 and real-world-3, while baseline MLP fails
most of the trials in all the environments including the benchmark-like
environment. Similarly, MLP policy trained with 250 environments can
successfully navigate in all the environments without any failures, while
baseline MLP trained with 50 environments fails most of the trials. Safe RL
improves the chances of success in all the environments and can navigate more
safely by performing backups and small adjustments of robots’ poses. Similar
to the simulation, MPC navigates very conservatively and succeeds in all the
trials in real-world-1 and real-world-2, but has much more difficulty
generalizing to large-scale real-world-3.

##  V Conclusion

In this section, we discuss the conclusions we draw from these benchmark
experiments. We organize these conclusions by the desiderata as follows:

(D1) reasoning under uncertainty of partially observed sensory inputs does not
obviously benefit from adding memory in simulated static environments and very
random dynamic (dynamic-box) environments, but much more significant
improvements were observed in the _real world_ and in more _challenging
dynamic environments_ (dynamic-wall).

(D2) safety is improved by both safe RL and model-based MPC methods. However,
classical navigation systems still achieve the best safety performance at a
cost of very long traversal time. Whether RL-based navigation systems can
achieve similar safety guarantees as classical navigation systems and whether
safety can be improved without significantly sacrificing the traversal time
are still open questions.

(D3) the ability to learn from limited trial-and-error data is not improved by
the evaluated model-based methods. Currently, we observe that model-based RL
methods indeed improve sample-efficiency, but only when the number of
imaginary rollouts from the learned model is large (e.g. $\geq 2000k$) and
when they are sampled with randomness. We therefore hypothesize that the
improvement comes from the robustness brought by learning on more data sampled
from the learned model. Hence, this result motivates not only more accurate
model learning for reducing the number of imaginary rollouts, but also
theoretical understanding of how the model helps improve the robustness or
even safety of navigation.

(D4) the generalization to diverse and novel environments is improved by
increasing the randomness of training environments. However, a noticeable gap
of about $5\%$ between training and test environments is not eliminated by
further increasing the number of training environments to 250. This reflects
the limitation of simple _domain randomization_ to increase the
generalization, which is, however, widely used by the community.

In summary, although the proposed benchmark is not intended to represent every
real-world navigation scenario, it serves as a simple yet comprehensive
testbed for RL-based navigation methods. We observed that for every
desideratum, no method can achieve 100% success rate on all _training_
environments. Even though we ensured that we have made sure that every
environment is indeed individually solvable. This alone indicates that there
exists an optimization and generalization challenge when we have a large
number of training environments as in our proposed benchmark.

## References

  * [1] S. Quinlan and O. Khatib, “Elastic bands: Connecting path planning and control,” in _[1993] Proceedings IEEE International Conference on Robotics and Automation_. IEEE, 1993, pp. 802–807.
  * [2] D. Fox, W. Burgard, and S. Thrun, “The dynamic window approach to collision avoidance,” _IEEE Robotics & Automation Magazine_, vol. 4, no. 1, pp. 23–33, 1997.
  * [3] X. Xiao, Z. Xu, Z. Wang, Y. Song, G. Warnell, P. Stone, T. Zhang, S. Ravi, G. Wang, H. Karnan _et al._ , “Autonomous ground navigation in highly constrained spaces: Lessons learned from the barn challenge at icra 2022,” _arXiv preprint arXiv:2208.10473_ , 2022.
  * [4] X. Xiao, B. Liu, G. Warnell, and P. Stone, “Motion planning and control for mobile robot navigation using machine learning: a survey,” _Autonomous Robots_ , pp. 1–29, 2022.
  * [5] Y. Chow, O. Nachum, A. Faust, M. Ghavamzadeh, and E. A. Duéñez-Guzmán, “Lyapunov-based safe policy optimization for continuous control,” _CoRR_ , vol. abs/1901.10031, 2019. [Online]. Available: http://arxiv.org/abs/1901.10031
  * [6] G. Thomas, Y. Luo, and T. Ma, “Safe reinforcement learning by imagining the near future,” 2022.
  * [7] E. Rodríguez-Seda, D. Stipanovic, and M. Spong, “Lyapunov-based cooperative avoidance control for multiple lagrangian systems with bounded sensing uncertainties,” in _2011 50th IEEE Conference on Decision and Control and European Control Conference, CDC-ECC 2011_ , ser. Proceedings of the IEEE Conference on Decision and Control, Dec. 2011, pp. 4207–4213, 2011 50th IEEE Conference on Decision and Control and European Control Conference, CDC-ECC 2011 ; Conference date: 12-12-2011 Through 15-12-2011.
  * [8] K. Cobbe, O. Klimov, C. Hesse, T. Kim, and J. Schulman, “Quantifying generalization in reinforcement learning,” in _ICML_ , 2019.
  * [9] K. Cobbe, C. Hesse, J. Hilton, and J. Schulman, “Leveraging procedural generation to benchmark reinforcement learning,” _arXiv preprint arXiv:1912.01588_ , 2019.
  * [10] N. Justesen, R. R. Torrado, P. Bontrager, A. Khalifa, J. Togelius, and S. Risi, “Illuminating generalization in deep reinforcement learning through procedural level generation,” _arXiv: Learning_ , 2018.
  * [11] J. Tobin, R. Fong, A. Ray, J. Schneider, W. Zaremba, and P. Abbeel, “Domain randomization for transferring deep neural networks from simulation to the real world,” in _2017 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2017, pp. 23–30.
  * [12] R. S. Sutton, “Dyna, an integrated architecture for learning, planning, and reacting,” _SIGART Bull._ , vol. 2, no. 4, p. 160–163, jul 1991. [Online]. Available: https://doi.org/10.1145/122344.122377
  * [13] A. Nagabandi, G. Kahn, R. S. Fearing, and S. Levine, “Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 7559–7566.
  * [14] M. J. Hausknecht and P. Stone, “Deep recurrent q-learning for partially observable mdps,” in _AAAI Fall Symposia_ , 2015.
  * [15] D. Wierstra, A. Förster, J. Peters, and J. Schmidhuber, “Solving deep memory pomdps with recurrent policy gradients,” in _ICANN_ , 2007.
  * [16] K. Chua, R. Calandra, R. McAllister, and S. Levine, “Deep reinforcement learning in a handful of trials using probabilistic dynamics models,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [17] H.-T. L. Chiang, A. Faust, M. Fiser, and A. Francis, “Learning navigation behaviors end-to-end with autorl,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 2007–2014, 2019.
  * [18] X. Xiao, B. Liu, G. Warnell, J. Fink, and P. Stone, “Appld: Adaptive planner parameter learning from demonstration,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4541–4547, 2020.
  * [19] Z. Wang, X. Xiao, B. Liu, G. Warnell, and P. Stone, “APPLI: Adaptive planner parameter learning from interventions,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021.
  * [20] Z. Wang, X. Xiao, G. Warnell, and P. Stone, “Apple: Adaptive planner parameter learning from evaluative feedback,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 7744–7749, 2021.
  * [21] Z. Xu, G. Dhamankar, A. Nair, X. Xiao, G. Warnell, B. Liu, Z. Wang, and P. Stone, “APPLR: Adaptive planner parameter learning from reinforcement,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021.
  * [22] X. Xiao, Z. Wang, Z. Xu, B. Liu, G. Warnell, G. Dhamankar, A. Nair, and P. Stone, “Appl: Adaptive planner parameter learning,” _Robotics and Autonomous Systems_ , vol. 154, p. 104132, 2022.
  * [23] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033.
  * [24] M. G. Bellemare, Y. Naddaf, J. Veness, and M. Bowling, “The arcade learning environment: An evaluation platform for general agents,” _Journal of Artificial Intelligence Research_ , vol. 47, pp. 253–279, jun 2013.
  * [25] H.-T. L. Chiang, A. Faust, M. Fiser, and A. Francis, “Learning navigation behaviors end-to-end with autorl,” _IEEE Robotics and Automation Letters_ , vol. 4, pp. 2007–2014, 2019.
  * [26] A. Wahid, A. Toshev, M. Fiser, and T.-W. E. Lee, “Long range neural navigation policies for the real world,” _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 82–89, 2019.
  * [27] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)(IEEE Cat. No. 04CH37566)_ , vol. 3. IEEE, 2004, pp. 2149–2154.
  * [28] Y. Zhu, R. Mottaghi, E. Kolve, J. J. Lim, A. Gupta, L. Fei-Fei, and A. Farhadi, “Target-driven visual navigation in indoor scenes using deep reinforcement learning,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 3357–3364.
  * [29] L. Harries, S. Lee, J. Rzepecki, K. Hofmann, and S. Devlin, “Mazeexplorer: A customisable 3d benchmark for assessing generalisation in reinforcement learning,” in _2019 IEEE Conference on Games (CoG)_. IEEE, 2019, pp. 1–4.
  * [30] Stanford Artificial Intelligence Laboratory et al., “Robotic operating system.” [Online]. Available: https://www.ros.org
  * [31] F. Xia, W. B. Shen, C. Li, P. Kasimbeg, M. E. Tchapmi, A. Toshev, R. Martín-Martín, and S. Savarese, “Interactive gibson benchmark: A benchmark for interactive navigation in cluttered environments,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 713–720, 2020\.
  * [32] S. Fujimoto, H. van Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” 2018.
  * [33] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International conference on machine learning_. PMLR, 2018, pp. 1861–1870.
  * [34] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.


# Discovering Sparse Representations of Lie Groups
with Machine Learning

Roy T. Forestano<EMAIL_ADDRESS>Institute for Fundamental Theory,
Physics Department, University of Florida, Gainesville, FL 32611, USA
Konstantin T. Matchev<EMAIL_ADDRESS>Institute for Fundamental Theory,
Physics Department, University of Florida, Gainesville, FL 32611, USA  Katia
Matcheva<EMAIL_ADDRESS>Institute for Fundamental Theory, Physics
Department, University of Florida, Gainesville, FL 32611, USA  Alexander Roman
<EMAIL_ADDRESS>Institute for Fundamental Theory, Physics Department,
University of Florida, Gainesville, FL 32611, USA  Eyup B. Unlu
<EMAIL_ADDRESS>Institute for Fundamental Theory, Physics Department,
University of Florida, Gainesville, FL 32611, USA  Sarunas Verner
<EMAIL_ADDRESS>Institute for Fundamental Theory, Physics Department,
University of Florida, Gainesville, FL 32611, USA

(February 10, 2023)

###### Abstract

Recent work has used deep learning to derive symmetry transformations, which
preserve conserved quantities, and to obtain the corresponding algebras of
generators. In this letter, we extend this technique to derive sparse
representations of arbitrary Lie algebras. We show that our method reproduces
the canonical (sparse) representations of the generators of the Lorentz group,
as well as the $U(n)$ and $SU(n)$ families of Lie groups. This approach is
completely general and can be used to find the infinitesimal generators for
any Lie group.

##  I Introduction

The transformative advances in theoretical physics at the turn of the 20th
century came with the realization that symmetry is a primary feature of nature
that constrains the allowable dynamical laws Gross (1996). By Noether’s
theorem Noether (1918), the presence of a continuous symmetry in the problem
implies a conservation law which is universally applicable and indispensable
in understanding the system’s behavior and evolution. Symmetries can be found
at all scales — from the microscopic description of subatomic particles in the
Standard Model to the large-scale structure of the Universe. In particle
physics, symmetries provide an organizing principle behind the observed
particle zoo and its interactions, and guide model-builders in the search for
viable extensions of the Standard Model Csáki  _et al._ (2018).

The formal study of symmetries has historically been done within the framework
of group theory Wigner  _et al._ (1959). While virtually all types of
classical Lie groups find some applications in various physics subfields, the
most commonly encountered Lie groups in particle physics are the special
orthogonal groups $SO(n)$ and the special unitary groups $SU(n)$. In
particular, the rotation group $SO(3)$ reflects the isotropy of space and all
fundamental laws of motion, the Lorentz group $SO(1,3)$ is the symmetry group
of spacetime in special relativity, and $SU(2)$ and $SU(3)$ are the Standard
Model gauge groups responsible for the weak and strong interactions,
respectively. An advanced course in group theory is a necessary tool in the
arsenal of the aspiring physics graduate student.

Recently there has been increased interest in the application of machine
learning (ML) for the discovery and identification of symmetries in data Iten
_et al._ (2020); Wetzel  _et al._ (2020); Krippendorf and Syvaeri (2020); Liu
and Tegmark (2021); Barenboim  _et al._ (2021); Dillon  _et al._ (2022); Liu
and Tegmark (2022); Desai  _et al._ (2022); Craven  _et al._ (2022); Moskalev
_et al._ (2022); Forestano  _et al._ (2023); Roman  _et al._ (2023). Relevant
applications of ML to group theory include computing tensor products and
branching rules of irreducible representations of Lie groups Chen  _et al._
(2021), and testing and/or deriving Lie group generators of a symmetry present
in the data Liu and Tegmark (2022); Moskalev  _et al._ (2022); Forestano  _et
al._ (2023); Roman  _et al._ (2023). In this letter we extend the method of
Forestano  _et al._ (2023); Roman  _et al._ (2023) to derive sparse
representations of arbitrary Lie algebras. We shall demonstrate that, with the
appropriate loss function modifications discussed in Section II, one can use
machine learning to derive from first principles the canonical (sparse) form
of the generators for the Lorentz group (Section III), as well as for the
$U(n)$ and $SU(n)$ families (Section IV). As emphasized in the summary Section
V, the described method is completely general and can be applied to any other
classical Lie group.

##  II Method

We begin by briefly summarizing the symmetry finding procedure outlined in
Ref. Forestano  _et al._ (2023). The classical groups are the linear groups of
transformations over the reals $\mathbb{R}$, the complex numbers $\mathbb{C}$,
and quaternions $\mathbb{H}$. Therefore, a symmetry transformation acts on a
feature vector ${\mathbf{x}}\equiv\\{x^{(1)},x^{(2)},\ldots,x^{(n)}\\}$, where
${\mathbf{x}}\in\mathbb{R}^{n}$ and ${\mathbf{x}}\in\mathbb{C}^{n}$ in our
$SO(n)$ and $SU(n)$ examples, respectively. In order to capture the effect of
a group transformation on ${\mathbb{R}}^{n}$ or ${\mathbb{C}}^{n}$, we
consider a representative set of $m$ points
$\left\\{\mathbf{x}\right\\}\equiv\left\\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{m}\right\\}$
sampled from some finite domain. The choice of a sampling distribution and/or
the size and location of the domain is inconsequential. For definiteness we
use a standard normal distribution with $m=100$.

The classical groups are defined in terms of polynomial invariants over their
respective fields. For example, $O(n)$ preserves the values of the polynomial
oracle

$\varphi_{O}(\mathbf{x})\;\equiv\;|\mathbf{x}|^{2}=\sum_{j=1}^{n}[x^{(j)}]^{2},\quad
x^{(j)}\in\mathbb{R},$ (1)

the Lorentz group in $n=4$ dimensions preserves

$\varphi_{L}(\mathbf{x})\;\equiv\;\bigl{(}x^{(1)}\bigr{)}^{2}-\bigl{(}x^{(2)}\bigr{)}^{2}-\bigl{(}x^{(3)}\bigr{)}^{2}-\bigl{(}x^{(4)}\bigr{)}^{2},\quad
x^{(j)}\in\mathbb{R},$ (2)

while $U(n)$ preserves

$\varphi_{U}(\mathbf{x})\;\equiv\;\sum_{j=1}^{n}\bigl{(}x^{(j)}\bigl{)}^{\ast}x^{(j)},\quad
x^{(j)}\in\mathbb{C}.$ (3)

A symmetry transformation $\mathbf{f}$ is a map
$\mathbf{x}^{\prime}=\mathbf{f}(\mathbf{x})$ which preserves the respective
oracle (1-3) everywhere, or in our case, for each of the sampled $m$ points:

$\varphi(\mathbf{x}^{\prime}_{i})\;\equiv\;\varphi(\mathbf{f}(\mathbf{x}_{i}))=\varphi(\mathbf{x}_{i}),\quad\forall
i\,=\,1,2,\ldots,m\,.$ (4)

In order to focus on the generators of the group of symmetry transformations,
$\mathbf{f}$ is linearized by considering infinitesimal transformations
$\delta{\mathbf{f}}$ in the vicinity of the identity transformation
$\mathbb{I}$:

$\delta{\mathbf{f}}\;\equiv\;\mathbb{I}+\varepsilon\,{\mathbb{G}}\,,$ (5)

where $\varepsilon$ is an infinitesimal parameter and ${\mathbb{G}}$ is a
$n\times n$ matrix. Upon training with a suitable loss function
$L({\mathbb{G}},\\{\mathbf{x}\\})$, the components of ${\mathbb{G}}$ are
driven to their trained values, thus producing one symmetry generator
Forestano  _et al._ (2023)

$\mathbb{J}\;\equiv\;\operatorname*{arg\,min}_{\mathbb{G}}\Bigl{(}L({\mathbb{G}},\\{\mathbf{x}\\})\Bigr{)}.$
(6)

By repeating the above procedure $N_{g}$ times under different initial
conditions ${\mathbb{G}}_{0}$ for the candidate generator components, or for
different random seeds and hyperparameters, one obtains a whole set of $N_{g}$
generators $\mathbb{J}_{\alpha}$, $\alpha=1,2,\ldots,N_{g}$. The loss function
is chosen to ensure that the obtained generators have the following
properties:

  1. 1.

Invariance, i.e., preserving the oracle values for all sampled datapoints
$\left\\{\mathbf{x}\right\\}$ under the set $\\{\mathbb{G}\\}$ of all
candidate transformations $\mathbb{G}_{\alpha}$, $\alpha=1,2,\ldots,N_{g}$:

$L_{\text{inv}}\bigl{(}\\{\mathbb{G}\\},\\{\mathbf{x}\\}\bigl{)}=\frac{1}{m\varepsilon^{2}}\sum_{\alpha=1}^{N_{g}}\sum_{i=1}^{m}\Bigl{[}\varphi\bigl{(}\mathbf{x}_{i}+\varepsilon{\mathbb{G}_{\alpha}}\cdot\mathbf{x}_{i}\bigr{)}-\varphi(\mathbf{x}_{i})\Bigr{]}^{2}\,,$
(7)

where “$\cdot$” denotes ordinary tensor multiplication.

  2. 2.

Orthogonality, which guarantees that the obtained generators are distinct:

$L_{\text{ortho}}\bigl{(}\\{\mathbb{G}\\}\bigr{)}=\Bigl{[}\sum_{\alpha=1}^{N_{g}-1}\sum_{\begin{subarray}{c}\beta=\alpha+1\end{subarray}}^{N_{g}}\text{Tr}\bigl{(}{\mathbb{G}}_{\alpha}\cdot{\mathbb{G}}_{\beta}^{T}\bigr{)}\Bigr{]}^{2}.$
(8)

  3. 3.

Normalization, ensuring that the transformation is not trivial:

$L_{\text{norm}}\bigl{(}\\{\mathbb{G}\\}\bigr{)}=\sum_{\alpha=1}^{N_{g}}\Bigl{[}\text{Tr}\Bigl{(}{\mathbb{G}}_{\alpha}\cdot{\mathbb{G}}^{T}_{\alpha}\Bigr{)}-2\Bigr{]}^{2}.$
(9)

  4. 4.

Closure of the algebra. This tests whether the set of candidate generators
$\\{\mathbb{G}\\}$ forms a closed algebra
$\bigl{[}\mathbb{G}_{\alpha},\mathbb{G}_{\beta}\bigr{]}=\sum_{\gamma=1}^{N_{g}}a_{\alpha\beta}^{\gamma}\mathbb{G}_{\gamma}$
with some structure constants $a_{\alpha\beta}^{\gamma}$:

$L_{\text{closure}}\bigl{(}\\{\mathbb{G}\\},a_{\alpha\beta}^{\gamma}\bigr{)}=\sum_{\alpha=1}^{N_{g}-1}\sum_{\begin{subarray}{c}\beta=\alpha+1\end{subarray}}^{N_{g}}\Tr\left(\mathbb{C}_{\alpha\beta}\cdot\mathbb{C}_{\alpha\beta}^{T}\right),$
(10)

where the closure mismatch is defined by

$\mathbb{C}_{\alpha\beta}\bigl{(}\\{\mathbb{G}\\},a_{\alpha\beta}^{\gamma}\bigr{)}\equiv\bigl{[}\mathbb{G}_{\alpha},\mathbb{G}_{\beta}\bigr{]}-\sum_{\gamma=1}^{N_{g}}a_{\alpha\beta}^{\gamma}\mathbb{G}_{\gamma}.$
(11)

The loss function in Forestano  _et al._ (2023) was the sum of those four
terms with suitable weight factors (hyperparameters) $h_{\text{ortho}}$,
$h_{\text{norm}}$ and $h_{\text{closure}}$:

$L_{4}=L_{\text{inv}}+h_{\text{ortho}}\,L_{\text{ortho}}+h_{\text{norm}}\,L_{\text{norm}}+h_{\text{closure}}\,L_{\text{closure}}.$
(12)

The analysis in Forestano  _et al._ (2023) demonstrated that this method is
capable of learning the symmetry groups (and their subgroups) generated by the
oracles (1) and (2) without any prior assumptions, in a fully agnostic
fashion. Since the method is completely general, however, the representations
learned in Forestano  _et al._ (2023) were not conveniently aligned with the
axes in the feature space, thus differing from the nice conventional results
found in the textbooks. To rectify this problem, here we introduce an
additional loss term designed to encourage sparsity:

$L_{\text{sparsity}}\bigl{(}\\{\mathbb{G}\\}\bigr{)}=\sum_{\alpha=1}^{N_{g}}\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{\begin{subarray}{c}j^{\prime}=1\end{subarray}}^{n}\sum_{\begin{subarray}{c}k^{\prime}=1\end{subarray}}^{n}\left|{\mathbb{G}}_{\alpha}^{(jk)}{\mathbb{G}}_{\alpha}^{(j^{\prime}k^{\prime})}\right|\Bigl{(}1-\delta_{jj^{\prime}}\delta_{kk^{\prime}}\Bigr{)},$
(13)

where ${\mathbb{G}}_{\alpha}^{(jk)}$ denotes the $jk$-component of
$\mathbb{G}_{\alpha}$. The full loss function is therefore

$L\bigl{(}\\{\mathbb{G}\\},a_{\alpha\beta}^{\gamma}\bigr{)}=L_{4}\bigl{(}\\{\mathbb{G}\\},a_{\alpha\beta}^{\gamma}\bigr{)}+h_{\text{sparsity}}\,L_{\text{sparsity}}\bigl{(}\\{\mathbb{G}\\}\bigr{)}$
(14)

and the minimization is performed over all generator components
${\mathbb{G}}_{\alpha}^{(jk)}$ and structure constants
$a_{\alpha\beta}^{\gamma}$, using the Adam optimizer with a learning rate
between $0.0001$ and $0.01$. Unless otherwise specified, the hyperparameters
$h$ were held fixed to 1.

##  III The Lorentz group $SO(1,3)$

As a reference benchmark, the left column of Figure 1 shows the result from
Forestano  _et al._ (2023) for the $N_{g}=6$ generators of the Lorentz group,
using the loss function (12) and the oracle (2). In this and all subsequent
such figures, each $n\times n$ panel represents a learned generator
${\mathbb{J}}_{\alpha}$ in matrix form. The values of the individual matrix
elements are indicated by the color bar. In the bottom panels of Figure 1,
each row (labeled $\alpha\beta=12,13,\ldots,56$) represents one of the 15
possible unique commutators $[\mathbb{J}_{\alpha},\mathbb{J}_{\beta}]$,
whereas the columns (labeled $\gamma=1,2,\ldots,6$) represent the found
generators $\mathbb{J}_{\gamma}$. Each cell then represents a structure
constant $a_{\alpha\beta}^{\gamma}$, with its value given by the color bar.

Figure 1: Left: the Lorentz group generators found in Ref. Forestano  _et al._
(2023) by using the loss function (12) and the oracle (2) (top panels) and the
corresponding structure constants (bottom panel). Right: the corresponding
results after inclusion of the sparsity term (13) with
$h_{\text{sparsity}}=1$.

While the six generators found in Forestano  _et al._ (2023) do satisfy the
four criteria 1-4 discussed in the previous section, they are not easily
interpretable, since the generators end up being generic linear combinations
of the familiar generators of boosts ${\mathbb{K}}_{k}$ and rotations
${\mathbb{L}}_{k}$, $k=1,2,3$. In other words, the representations in the top
left panels in Figure 1 are not sparse. After including the sparsity term (13)
in the loss function (14), we obtain the corresponding result shown in the
right panels of Figure 1. We easily recognize
$\mathbb{J}_{1}=-\mathbb{K}_{2}$, $\mathbb{J}_{3}=-\mathbb{K}_{3}$ and
$\mathbb{J}_{4}=\mathbb{K}_{1}$ as the canonical boost generators and
$\mathbb{J}_{2}=\mathbb{L}_{3}$, $\mathbb{J}_{5}=-\mathbb{L}_{2}$ and
$\mathbb{J}_{6}=-\mathbb{L}_{1}$ as the canonical generators of rotations.

Figure 2: The learned generators of the $N_{g}=4$ subalgebra of the Lorentz
group for different values of the $h_{\text{sparsity}}$ hyperparameter:
$h_{\text{sparsity}}=0,0.01,0.1,1,10$ (from top to bottom). Each training
session ran for 3,000 epochs.

Like any ML approach, our method relies on the minimization of a loss
function, which involves some hyperparameter tuning. For optimal results, we
do need to successfully balance the relative weights of the different terms in
the loss function (14). Figure 2 illustrates the dependence of the learned
generators for the $N_{g}=4$ subalgebra of the Lorentz group on the sparsity
hyperparameter $h_{\text{sparsity}}$, which is scanned from 0 to 10. As the
sparsity term is gradually turned on, the learned representations become
sparser, and for a certain range of values of $h_{\text{sparsity}}$, we
recover a canonical set $\mathbb{J}_{1}=\mathbb{L}_{3}$,
$\mathbb{J}_{2}=\mathbb{K}_{2}+\mathbb{L}_{1}$,
$\mathbb{J}_{3}=\mathbb{K}_{1}-\mathbb{L}_{2}$, and
$\mathbb{J}_{4}=\mathbb{K}_{3}$. However, increasing $h_{\text{sparsity}}$
further leads to solutions which are very sparse, but only at the cost of
violating one of the other conditions, in this case normalization — note that
$\mathbb{J}_{4}$ in the last row of Figure 2 is essentially trivial.

##  IV The $U(n)$ and $SU(n)$ family

We now turn our attention to the $U(n)$ and $SU(n)$ families which were not
considered in Forestano  _et al._ (2023) or elsewhere in the current context.
The two families are related as $U(n)=SU(n)\times U(1)$, where the $U(1)$
factor accounts for the $\det U=1$ condition. Since the invariant polynomial
(3) is defined over the field of complex numbers, we need to make a few
adjustments to the procedure in Section II: i) create a dataset
$\\{\mathbf{x}\\}$ which now consists of $m$ $n$-dimensional complex vectors;
ii) use the complex form of the functions defined in the loss function by
replacing any products of $\mathbb{G}\mathbb{G}$ with
$\mathbb{G}\mathbb{G}^{\ast}$; and iii) insert the conventional imaginary unit
factor in front of $\varepsilon$ in (5) as
$\delta{\mathbf{f}}=\mathbb{I}+i\varepsilon\,{\mathbb{G}}$. Upon training with
$n^{2}$ generators, we will find the algebra of $U(n)$. In order to find the
algebra of $SU(n)$, we take advantage of the ability of our method to find
sets of generators forming closed algebras. Therefore, all we need to do is
reduce the number of generators to $n^{2}-1$ and proceed as before.

###  IV.1 The case of $U(2)$ and $SU(2)$

Figure 3: The learned sparse generators for the case of $U(2)$ (left) and
$SU(2)$ (right).

Our results for the case of $n=2$ are shown in Figure 3. Since the components
of the generator matrices are now complex, we show separate panels for their
real ($\mathfrak{R}$) and imaginary ($\mathfrak{I}$) parts. The generators of
$U(2)$ are comprised of the generators of $SU(2)$, usually taken as the Pauli
matrices

$\sigma_{1}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\qquad\sigma_{2}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\qquad\sigma_{3}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix},$
(15)

supplemented with the $2\times 2$ identity matrix
$\sigma_{0}={\mathbb{I}}_{2}$. We see that for the $N_{g}=4$ case of $U(2)$
displayed on the left side of Figure 3, the method correctly finds
$\mathbb{J}_{1}=-\sigma_{1}$, $\mathbb{J}_{3}=-\sigma_{2}$,
$\mathbb{J}_{2}=(\sigma_{0}+\sigma_{3})/\sqrt{2}$,
$\mathbb{J}_{4}=(\sigma_{0}-\sigma_{3})/\sqrt{2}$. The last two linear
combinations are chosen by the algorithm because they are sparser than the
canonical set $\sigma_{0}$ and $\sigma_{3}$. Note that when a generator matrix
has a single nonvanishing element, its value is fixed to $\pm\sqrt{2}$ by the
normalization condition (9). When repeating the exercise with $N_{g}=3$, as is
the case of $SU(2)$, we find the result on the right in Figure 3. We recognize
$\mathbb{J}_{1}=\sigma_{3}$, $\mathbb{J}_{2}=\sigma_{1}$,
$\mathbb{J}_{3}=-\sigma_{2}$, which is precisely the canonical $SU(2)$ algebra
up to signs and relabelling.

###  IV.2 The case of $U(3)$ and $SU(3)$

Figure 4: The same as Figure 3, but for the case of $U(3)$ (top panels) and
$SU(3)$ (bottom panels).

We now move on to the case of $n=3$. The generators of $SU(3)$ are the Gell-
Mann matrices which are the $SU(3)$ analogues of the Pauli matrices (15)

$\displaystyle\lambda_{1}$ $\displaystyle=\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix},\qquad\lambda_{2}=\begin{pmatrix}0&-i&0\\\ i&0&0\\\
0&0&0\end{pmatrix},\qquad\lambda_{3}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&0\end{pmatrix},$ $\displaystyle\lambda_{4}$
$\displaystyle=\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix},\qquad\lambda_{5}=\begin{pmatrix}0&0&-i\\\ 0&0&0\\\
i&0&0\end{pmatrix},$ (16) $\displaystyle\lambda_{6}$
$\displaystyle=\begin{pmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{pmatrix},\qquad\lambda_{7}=\begin{pmatrix}0&0&0\\\ 0&0&-i\\\
0&i&0\end{pmatrix},\qquad\lambda_{8}=\frac{1}{\sqrt{3}}\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&-2\end{pmatrix}.$

Our results are shown in Figure 4. Upon inspection, we see that the generators
are successfully reproduced. In particular, for $SU(3)$, the algorithm found
$\mathbb{J}_{1}=-\lambda_{4}$,
$\mathbb{J}_{2}=(\lambda_{8}-\sqrt{3}\lambda_{3})/2$,
$\mathbb{J}_{3}=\lambda_{6}$, $\mathbb{J}_{4}=\lambda_{2}$,
$\mathbb{J}_{5}=\lambda_{7}$, $\mathbb{J}_{6}=\lambda_{1}$,
$\mathbb{J}_{7}=-\lambda_{5}$,
$\mathbb{J}_{8}=-(\lambda_{3}+\sqrt{3}\lambda_{8})/2$.

###  IV.3 The case of $U(4)$ and $SU(4)$

Figure 5: The same as Figure 4, but for the case of $U(4)$ (top panels) and
$SU(4)$ (bottom panels).

Figure 5 shows the results for the next most complicated case, $n=4$. The
generators of $SU(4)$ are the corresponding generalizations of the Gell-Mann
matrices (16) and consist of 6 symmetric matrices with real off-diagonal
components, 6 anti-symmetric matrices with purely imaginary off-diagonal
components, and 3 matrices with real diagonal components. This is precisely
what we observe in the lower panels of Figure 5. The $U(4)$ algebra has one
additional diagonal generator, which can be taken as the identity matrix
$\mathbb{I}_{4}$. Once again, the sparsity term in the loss function forces
the $U(4)$ generators with real diagonal components to have a single non-zero
entry, see $\mathbb{J}_{1}$, $\mathbb{J}_{4}$, $\mathbb{J}_{7}$, and
$\mathbb{J}_{10}$.

###  IV.4 The case of $SU(5)$ and beyond

Figure 6: The learned sparse generators for the case of $SU(5)$ (top panels)
and $SU(6)$ (bottom panels).

We now turn to the cases of $n=5$ and $n=6$. The learned generators for
$SU(5)$ and $SU(6)$ are shown in Figure 6. In the case of $SU(5)$, we have 4
generators with real diagonal entries ($\mathbb{J}_{1}$, $\mathbb{J}_{2}$,
$\mathbb{J}_{6}$ and $\mathbb{J}_{20}$), 10 symmetric matrices with real off-
diagonal components ($\mathbb{J}_{\alpha}$ with
$\alpha\in\\{3,4,5,9,10,13,15,17,22,24\\}$), and 10 anti-symmetric matrices
with imaginary off-diagonal components ($\mathbb{J}_{\alpha}$ with
$\alpha\in\\{7,8,11,12,14,16,18,19,21,23\\}$). In the case of $SU(6)$, out of
the 35 total generators, we have 5 generators with real diagonal entries
($\mathbb{J}_{\alpha}$ with $\alpha\in\\{2,8,10,13,24\\}$), 15 symmetric
matrices with real off-diagonal components ($\mathbb{J}_{\alpha}$ with
$\alpha\in\\{1,4,5,7,9,11,15,21,23,25,28,29,30,31,34\\}$), and 15 anti-
symmetric matrices with imaginary off-diagonal components
($\mathbb{J}_{\alpha}$ with
$\alpha\in\\{3,6,12,14,16,17,18,19,20,22,26,27,32,33,35\\}$). Having
demonstrated the successful discovery of the $SU(5)$ and $SU(6)$ generators in
Figure 6, we believe that our procedure has proven to be fully generalizable
for any higher dimensional groups $SU(n)$ with $n>6$.

##  V Summary and Outlook

In this letter we demonstrated how to use machine learning to derive sparse
representations of the generators of an arbitrary Lie algebra. The key
observation is that a Lie group preserves a certain invariant polynomial under
the group action on the vector space (for the examples considered here, those
were the oracle functions (1-3)). As shown in Forestano  _et al._ (2023), one
can design suitable loss functions for the generator components, so that
minimization of the loss function leads to a set of valid generators. In this
letter we extended the approach of Forestano  _et al._ (2023) by adding terms
to the loss function that encourage finding sparse representations. We
validated our method for the Lorentz group (in Section III) and for the $U(n)$
and $SU(n)$ family of Lie groups (in Section IV). Although not shown here, we
have also verified the procedure on all of the $SO(n)$ examples considered in
Forestano  _et al._ (2023).

We should note that all of our numerical work was done on a personal laptop,
and did not use any high performance facilities. All models were trained for
5,000 epochs and the typical training times ranged from a few minutes to a
couple of hours for the case of $SU(6)$. While the method is applicable to
arbitrarily large groups, in practice the training in such cases may become
relatively slow and could benefit from parallelization, further code
optimization, further hyperparameter tuning, and better understanding of the
loss landscape. Examples of much larger groups, like the exceptional Lie
groups Ramond (1976), will be explored in future work.

###### Acknowledgements.

We thank A. Davis, S. Gleyzer, R. Houtz, K. Kong, S. Mrenna, H. Prosper and P.
Shyamsundar for useful discussions. We thank P. Ramond for group theory
insights and inspiration. This work is supported in part by the U.S.
Department of Energy award number DE-SC0022148.

## References

  * Gross (1996) David J. Gross, “The role of symmetry in fundamental physics,” PNAS 93, 14256–14259 (1996).
  * Noether (1918) E. Noether, “Invariante variationsprobleme,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1918, 235–257 (1918).
  * Csáki  _et al._ (2018) Csaba Csáki, Salvator Lombardo, and Ofri Telem, “TASI Lectures on Non-supersymmetric BSM Models,” in  _Proceedings, Theoretical Advanced Study Institute in Elementary Particle Physics : Anticipating the Next Discoveries in Particle Physics (TASI 2016): Boulder, CO, USA, June 6-July 1, 2016_, edited by Rouven Essig and Ian Low (WSP, 2018) pp. 501–570, arXiv:1811.04279 [hep-ph] .
  * Wigner  _et al._ (1959) E.P. Wigner, J.J. Griffin, and J.J. Griffin,  _Group Theory and Its Application to the Quantum Mechanics of Atomic Spectra_ , Pure and Applied Physics : a series of monographs and textbooks (Academic Press, 1959).
  * Iten  _et al._ (2020) Raban Iten, Tony Metger, Henrik Wilming, Lídia del Rio, and Renato Renner, “Discovering Physical Concepts with Neural Networks,” Phys. Rev. Lett. 124, 010508 (2020), arXiv:1807.10300 [quant-ph] .
  * Wetzel  _et al._ (2020) Sebastian J. Wetzel, Roger G. Melko, Joseph Scott, Maysum Panju, and Vijay Ganesh, “Discovering Symmetry Invariants and Conserved Quantities by Interpreting Siamese Neural Networks,” Phys. Rev. Res. 2, 033499 (2020), arXiv:2003.04299 [physics.comp-ph] .
  * Krippendorf and Syvaeri (2020) Sven Krippendorf and Marc Syvaeri, “Detecting Symmetries with Neural Networks,” (2020), arXiv:2003.13679 [physics.comp-ph] .
  * Liu and Tegmark (2021) Ziming Liu and Max Tegmark, “Machine Learning Conservation Laws from Trajectories,” Phys. Rev. Lett. 126, 180604 (2021), arXiv:2011.04698 [cs.LG] .
  * Barenboim  _et al._ (2021) Gabriela Barenboim, Johannes Hirn, and Veronica Sanz, “Symmetry meets AI,” SciPost Phys. 11, 014 (2021), arXiv:2103.06115 [cs.LG] .
  * Dillon  _et al._ (2022) Barry M. Dillon, Gregor Kasieczka, Hans Olischlager, Tilman Plehn, Peter Sorrenson, and Lorenz Vogel, “Symmetries, safety, and self-supervision,” SciPost Phys. 12, 188 (2022), arXiv:2108.04253 [hep-ph] .
  * Liu and Tegmark (2022) Ziming Liu and Max Tegmark, “Machine Learning Hidden Symmetries,” Phys. Rev. Lett. 128, 180201 (2022), arXiv:2109.09721 [cs.LG] .
  * Desai  _et al._ (2022) Krish Desai, Benjamin Nachman, and Jesse Thaler, “Symmetry discovery with deep learning,” Phys. Rev. D 105, 096031 (2022), arXiv:2112.05722 [hep-ph] .
  * Craven  _et al._ (2022) Sean Craven, Djuna Croon, Daniel Cutting, and Rachel Houtz, “Machine learning a manifold,” Phys. Rev. D 105, 096030 (2022), arXiv:2112.07673 [hep-ph] .
  * Moskalev  _et al._ (2022) Artem Moskalev, Anna Sepliarskaia, Ivan Sosnovik, and Arnold Smeulders, “Liegg: Studying learned lie group generators,”  (2022), arXiv:2210.04345 [cs.LG] .
  * Forestano  _et al._ (2023) Roy T. Forestano, Konstantin T. Matchev, Katia Matcheva, Alexander Roman, Eyup Unlu, and Sarunas Verner, “Deep Learning Symmetries and Their Lie Groups, Algebras, and Subalgebras from First Principles,” (2023), arXiv:2301.05638 [hep-ph] .
  * Roman  _et al._ (2023) Alexander Roman, Roy T. Forestano, Konstantin T. Matchev, Katia Matcheva, and Eyup B. Unlu, “Oracle-Preserving Latent Flows,” (2023), arXiv:2302.00806 [cs.LG] .
  * Chen  _et al._ (2021) Heng-Yu Chen, Yang-Hui He, Shailesh Lal, and Suvajit Majumder, “Machine learning Lie structures & applications to physics,” Phys. Lett. B 817, 136297 (2021), arXiv:2011.00871 [hep-th] .
  * Ramond (1976) Pierre Ramond,  _Introduction to Exceptional Lie Groups and Algebras_ , Tech. Rep. CALT-68-577 (California Institute of Technology, Pasadena, California 91125, 1976).


# Physics-inspired Ising Computing with Ring Oscillator Activated p-bits

‡ Navid Anjum Aadit∗, ‡ Andrea Grimaldi†, Giovanni Finocchio†,∘
and Kerem Y. Camsari∗,♮ $\textsection$‡ Equally contributing authors
†Department of Mathematical and Computer Sciences, Physical Sciences and Earth
Sciences,
University of Messina, Messina, Italy
∗Department of Electrical and Computer Engineering, University of California,
Santa Barbara,
Santa Barbara, CA, 93106, USA

###### Abstract

The nearing end of Moore’s Law has been driving the development of domain-
specific hardware tailored to solve a special set of problems. Along these
lines, probabilistic computing with inherently stochastic building blocks
(p-bits) have shown significant promise, particularly in the context of hard
optimization and statistical sampling problems. p-bits have been proposed and
demonstrated in different hardware substrates ranging from small-scale
stochastic magnetic tunnel junctions (sMTJs) in asynchronous architectures to
large-scale CMOS in synchronous architectures. Here, we design and implement a
truly asynchronous and medium-scale p-computer (with $\approx$ 800 p-bits)
that closely emulates the asynchronous dynamics of sMTJs in Field Programmable
Gate Arrays (FPGAs). Using hard instances of the planted Ising glass problem
on the Chimera lattice, we evaluate the performance of the asynchronous
architecture against an ideal, synchronous design that performs parallelized
(chromatic) exact Gibbs sampling. We find that despite the lack of any careful
synchronization, the asynchronous design achieves parallelism with comparable
algorithmic scaling in the ideal, carefully tuned and parallelized synchronous
design. Our results highlight the promise of massively scaled p-computers with
millions of free-running p-bits made out of nanoscale building blocks such as
stochastic magnetic tunnel junctions.

###### Index Terms:

p-bits, combinatorial optimization, planted Ising, Chimera lattice,
asynchronous computing, massive parallelism, magnetic tunnel junctions

§§footnotetext: Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  I Introduction

With the nearing end of Moore’s Law, domain-specific hardware and
architectures are growing rapidly. The notion of performing some tasks more
efficiently (area, speed and/or energy) rather than improving performance for
general purpose computing has led to the proliferation of special-purpose
accelerators. With their widespread use, hard optimization problems have been
a primary target of this approach and a variety of different domain-specific
hardware architectures have emerged (see, Ref. [1] for a general and recent
review).

As an example of this growing trend, probabilistic bits or p-bits were
introduced [2] as a building block which can accelerate a broad family of
algorithms including Monte Carlo, Markov Chain Monte Carlo [3], Quantum Monte
Carlo, statistical sampling for Bayesian inference and Boltzmann machine
learning [4] methods. p-bits have been shown to be compatible with powerful
optimization techniques such as parallel tempering [5] with competitive
performance relative to all other Ising machines (classical and quantum) in
select problems such as integer factorization and Boolean satisfiability [6].
Their combination with sophisticated algorithms [7] could yield further
advantages.

A natural advantage of the p-bit model is its native mapping to the Ising
model and to the natural generalization of Ising Models. This ensures that
coupled p-bits can systematically probe the exact Boltzmann distribution
through Gibbs or Metropolis sampling without any approximations or reductions,
often necessary in alternative, non-bistable abstractions of the Ising spin.

One particularly promising small-scale demonstration of p-bits in an
asynchronously operating mode was performed in Ref. [8]. Combined with key
breakthrough experiments demonstrating nanosecond fluctuations in suitably
designed low barrier magnetic tunnel junctions (MTJ) [9, 10], these results
suggest the intriguing possibility of designing $>$ million bit probabilistic
computers [11] in light of the remarkable advances in the magnetic memory chip
industry reaching gigabit densities [12, 13]. Even though large scale p-bit
emulators have been designed and tested in FPGAs or ASICs, [14, 15, 11, 6, 3],
virtually all of these implementations have been on synchronous hardware where
a global clock controlled the information flow.

FIG. 1: Upper panel: Physics-inspired analogy between asynchronously
interacting bodies and p-bits: both systems are asynchronous, local (sparse
connectivity) and massively parallel. (a) Asynchronous computer architecture:
The local field (Eq. (1)) for each p-bit is computed combinationally. Each
p-bit has a different clock with asynchronous activation (b) Asynchronous ring
oscillator-based clocks (ROSC) where we draw an analogy between stochastic
Magnetic Tunnel Junction (sMTJ)-based p-bits and ROSCs. Large rings correspond
to large energy-barrier sMTJs with higher retention times. (c) Measured
frequency range for ROSCs with different ring sizes in the FPGA
implementation.

In this paper, we make a first attempt in designing and building a physics-
inspired, truly asynchronous architecture, closely emulating the dynamics of
interacting nanodevice-based p-bits. This physics-inspired architecture bears
similarities to locally interacting (sparsely connected) and asynchronous
bodies with probabilistic dynamics (FIG. 1, upper panel). We achieve the
design by an unconventional use of FPGAs where individual p-bits are activated
by decoupled ring oscillators and can have overlapping and out-of-phase clocks
with different frequencies. Considering how variations may influence
individual p-bit behavior in magnetic tunnel junction based designs [16] the
behavior of asynchronous p-computers with built-in variations is worth
investigating.

To compare the performance of the truly asynchronous p-computer in the FPGA,
we choose the planted Ising model where a hard optimization problem is
generated with a planted solution [17, 18], allowing a reliable evaluation of
the asynchronous design with respect to exact Gibbs sampling.

##  II Physics-inspired Architecture

The main equations of the p-bit model (FIG. 1a) involves stochastic activation
and a local field (synapse) calculation, given by:

$m_{i}=\mathrm{sgn}(\mathrm{tanh}(\beta I_{i})-r_{U})\quad I_{i}=\sum
J_{ij}m_{j}+h_{i}$ (1)

where $m_{i}$ represents the bipolar p-bit state ($\pm 1$), $r_{U}$ is a
uniform random number between $(-1,+1)$ and $[J],\\{h\\}$ are the weights and
biases for a given problem and $\beta$ is the inverse temperature.

Standard Gibbs sampling iterates Eq. (1) to reach the Boltzmann distribution
defined by the weights and typically involves a serialized update procedure
with nested for loops. One way to avoid this serial for loop is to perform
block updates between unconnected p-bits. This approach when applied in
software is named “chromatic sampling” [19] and a low-level hardware
realization of it was recently reported in Ref. [6]. However, this design also
involves carefully designed and equally phase shifted synchronous clocks so
that multiple blocks do not update simultaneously.

In this work, inspired by truly asynchronous small-scale implementations of
p-computers with nanodevices (based on stochastic MTJs [8, 4]), we implemented
a physics-inspired, truly asynchronous Ising Computer where different p-bits
receive clocks with different frequencies with random phases. In contrast to
synchronous designs, no careful engineering between the clocks of asynchronous
p-bits were made. Moreover, unavoidable variations of sMTJs in highly scaled
p-computers with nanodevices would make such engineering extremely difficult
if not impossible. We found that despite the deliberate randomization of p-bit
clocks and unavoidable collisions breaking exact Gibbs sampling, the physics-
inspired design exhibited massive parallelism observed in carefully tuned
synchronous designs, not observed in standard CPU-based Gibbs sampling (FIG.
3).

FIG. 2: Planted Ising problems. (a) Chimera lattice with different sizes,
depending on number of tiles. Different tiles were chosen in the same hardware
to change problem sizes, similar to Ref.[17] but all problems used a fixed
800-spin Chimera in our hardware. (b) Loop/close generation process shown for
$c_{l}=[6,8]$. (c) Example graph with three loops (the colored hatchings) and
a planted solution (the colors of the nodes). For each loop, the weights of
the couplings between neighbors are assigned to be either ferromagnetic (FM)
interactions or antiferromagnetic (AFM), based on the plant. One randomly
selected coupling in each loop is flipped (spark symbol).

Ring Oscillator Generation: A ROSC clock consists of an odd number of looped
NOT gates. In our FPGA (Xilinx, VCU118), we attach controllable delays to our
inverters to make logical delays comparable to wire delays (FIG. 1b). We
designed the delay unit as a flip flop with a very fast master clock ($300$
MHz) compared to the ROSC frequencies that essentially acts as a combinational
delay unit. In this way, we were able to obtain highly regular ROSC clocks as
a function of ring sizes whose frequencies were measured by specially designed
counters (FIG. 1c). In our experiments, we used 10 ROSCs to drive 800 p-bits
in a Chimera lattice. Each p-bit has a pseudorandom number generator, which is
a 32-bit Linear Feedback Shift Register (LFSR). The ROSCs activate the LFSRs
of the p-bits randomly based on the frequencies. In this work, we have
distributed the clocks evenly among the p-bits between 5 and 17 MHz. However,
different distributions for the clocks, e.g., Gaussian, could be used. Since
the Chimera graph is bi-partite, we did not assign the same clock to two
p-bits that are on the different partitions to avoid systematic parallel
updates between connected p-bits. Future work will consider dynamic clocking
schemes where each p-bit can have a different “retention time” much like MTJ-
based p-bits.

##  III Planted Ising Model

An important class of hard optimization problems are those with “planted”
ground states that allow effective evaluation of performance. We construct
frustrated spin glasses with planted solutions [17] on a 800-spin Chimera
graph where we changed the number of tiles for different problem sizes (FIG.
2a). A Hamiltonian generated by this process is the sum of several local
frustrated Hamiltonians which we will call “clauses”. A planted solution will
be used to define these clauses so that it will be the ground energy of the
final Hamiltonian. Every instance is characterized by two parameters: the
clause density $\alpha$, defined as $n_{c}/n_{n}$, where $n_{c}$ is the number
of clauses and $n_{n}$ is the number of nodes of the graph, and the length of
possible loops that form the clauses, $c_{l}=[l_{min},l_{max}]$, where
$l_{min/max}$ is the min/max loop length, respectively. In this work, we chose
$\alpha=0.4$ and $c_{l}=[4,8]$ for all our instances used in this paper, that
run on the same 800-spin Chimera lattice in our hardware.

Clause generation: A total of $n_{c}=\alpha n_{n}$ clauses is generated. Each
clause is an ordered sequence of nodes that creates a loop of acceptable
length in the graph. To obtain one, following Ref. [17], we pick a random node
and start a non-backtracking random walk of at most $l_{max}$ steps. If the
walker lands on an already visited node, it means that a loop was formed and
the node can be considered its initial point. If the length of the loop is
$>l_{min}$ the clause is accepted, if it is not or if the maximum number of
steps is reached without closing the loop, the process is repeated. FIG. 2b
shows a few examples of this process. A planted solution $s$ is generated by
creating a random array of $-1$s and $+1$s of length $n_{n}$. A clause can be
defined as $c_{m}=\\{n_{1},n_{2},...,n_{k},n_{k+1}\\}$, with
$k\in[l_{min},l_{max}]$, where $n_{k+1}=n_{1}$, representing the closing of
the loop. Now, $\forall i\in[1,k]:i\neq j$ we increase $J_{n_{i},n_{i+1}}$ by
$s_{i}s_{i+1}$, while for $j\in[1,k]$, picked at random, we increase
$J_{n_{j},n_{j+1}}$ by $-s_{j}s_{j+1}$. This last step serves to create a
frustrated loop. Once this is done for all clauses, the final $J$ is
calculated by summing $J$ and $J^{T}$ and by normalizing so that all $J$ lie
between $[-1,+1]$.

##  IV Performance Comparison

We follow the time-to-solution formulation[17, 18] to measure performance of
the physics-inspired asynchronous architecture.

$\mathrm{TTS}(\tau,p_{R})=\tau N_{r}(\tau,p_{R})$ (2)

where $N_{r}$ is the expected number of repetitions we need to perform an
annealing schedule of time $\tau$ to the energy ground state at least once
with probability $p_{R}$. $N_{r}$ is defined as:

$N_{r}(\tau,p_{R})=\frac{\ln{(1-p_{R})}}{\ln{(1-p_{S}(\tau))}}$ (3)

FIG. 3: Performance comparison of synchronous (graph colored-tuned),
asynchronous (physics-inspired, not colored, not tuned) and standard Gibbs
(CPU) samplers on the planted Ising problem.

where $p_{S}$ is the probability of success in finding the ground state in one
annealing process of length $\tau$.

To evaluate the performance of our asynchronous architecture, we compared it
to serialized Gibbs sampling on CPU ($2.6$ GHz) and to synchronous colored
Gibbs sampling on FPGA. We investigated the scaling difficulty of planted
Ising instances with fixed $c_{l}$ and $\alpha$ across several Chimera graphs
increasing in size by changing the number of tiles used in a Chimera, as
illustrated in FIG. 2. For each set of tiles, we generated 100 planted Ising
instances, performing simulated annealing trials to estimate the success
probability, $p_{S}(\tau)$, for each instance (FIG. 3). On FPGA, we performed
500 trials per instance while on CPU we only performed 50 trials because of
the exceedingly long run-times.

After estimating a $p_{S}$ for every instance, the average $\mathrm{TTS}$ for
a point is obtained by simply calculating the average of the $p_{S}$ for
instances of that size and applying Eq. (2) with the appropriate annealing
time $\tau$. The reference probability $p_{R}$ was set to $99\%$. The error
bars are obtained through bootstrapping with 95% confidence intervals and
$10^{4}$ samples.

We present the $\mathrm{TTS}$ for solving 100 instances of the planted Ising
problems of 9 different sizes in FIG. 3. The standard Gibbs (CPU) was
implemented in Python using optimized libraries for matrix calculations. The
final two points were not computed because of time limitations. We solved the
exact same instances on the FPGA programmed with the asynchronous ROSC
activated 800 p-bits with a fixed point representation using 10-bits. As a
reference, the synchronous solver performing chromatic Gibbs sampling solves
the same instances on the same FPGA where careful phase shifting ensures no
simultaneous or incorrect updates (where $I_{i}$ calculation is not complete)
occur between neighboring p-bits (as in Ref. [6]). On the other hand, the
asynchronous solver is expected to take samples with both of those errors when
p-bit clock edges happen to be closely separated. Our experiment investigates
the usefulness of such samples.

We defined a common linear simulated annealing schedule for all the
architectures with $\beta=0.5$ to $7.0$ with a step of $0.5$ where, at each
$\beta$, a total of $937$ sweeps (attempted flips of all p-bits) are executed.
In the FPGA the annealing time was fixed to $\tau=1.4$ ms for each trial. To
obtain a configuration comparable to the asynchronous architecture, the
synchronous architecture was set up with two stable and oppositely phase
shifted clocks with the average frequency (9.375 MHz) of the 10 ROSC clocks.
We believe this arrangement made the two designs equivalent beyond the
asynchronous and inexact dynamics of the ROSC since both designs approximately
take the same amount of samples within the fixed annealing time $\tau$. The
key result we obtained is shown in FIG. 3. We observe a clear scaling
difference between the CPU implementation of standard serialized Gibbs
sampling and the massively parallel FPGA implementations which (ideally)
obtain a scaling factor of $\approx N$ in their flips/second due to their
massively parallel architecture. Both solvers provide a roughly 5-orders of
magnitude prefactor improvement over the CPU. Intriguingly, the scaling of the
synchronous and asynchronous FPGA remain similar, despite the possibility of
many collisions (parallel or incorrect updates) in the asynchronous design.
Indeed, the carefully tuned synchronous design performs strictly better than
the asynchronous one in all instances. Nevertheless, it is encouraging to
observe that the asynchronous design without any carefully engineered clocks
or tuning performs nearly as well, leading to the promising possibility of
truly asynchronous, million bit p-computers with stochastic MTJs or other
nanodevices.

## Acknowledgment

K.Y.C. and N.A.A. acknowledge support through National Science Foundation (CCF
2106260) and K.Y.C. through the Samsung GRO program. A.G. and G.F. were
supported under the project PRIN 2020LWPKH7 funded by the Italian Ministry of
University and Research and by the PETASPIN Association (www.petaspin.com).

## References

  * [1] N. Mohseni, P. L. McMahon, and T. Byrnes, “Ising machines as hardware solvers of combinatorial optimization problems,” _arXiv preprint arXiv:2204.00276_ , 2022.
  * [2] K. Y. Camsari et al., “Stochastic p-bits for invertible logic,” _Physical Review X_ , vol. 7, no. 3, p. 031014, 2017.
  * [3] J. Kaiser et al., “Benchmarking a probabilistic coprocessor,” _arXiv preprint arXiv:2109.14801_ , 2021.
  * [4] J. Kaiser et al., “Hardware-aware in situ learning based on stochastic magnetic tunnel junctions,” _Physical Review Applied_ , vol. 17, no. 1, p. 014016, 2022.
  * [5] A. Grimaldi et al., “Spintronics-compatible approach to solving maximum-satisfiability problems with probabilistic computing, invertible logic, and parallel tempering,” _Physical Review Applied_ , vol. 17, no. 2, p. 024052, 2022.
  * [6] N. A. Aadit et al., “Massively parallel probabilistic computing with sparse ising machines,” _arXiv preprint arXiv:2110.02481_ , 2021.
  * [7] M. Mohseni et al., “Nonequilibrium monte carlo for unfreezing variables in hard combinatorial optimization,” _arXiv preprint arXiv:2111.13628_ , 2021\.
  * [8] W. A. Borders et al., “Integer factorization using stochastic magnetic tunnel junctions,” _Nature_ , 2019.
  * [9] K. Hayakawa et al., “Nanosecond random telegraph noise in in-plane magnetic tunnel junctions,” _Physical Review Letters_ , vol. 126, no. 11, p. 117202, 2021.
  * [10] C. Safranski et al., “Demonstration of nanosecond operation in stochastic magnetic tunnel junctions,” _Nano Letters_ , vol. 21, no. 5, pp. 2040–2045, 2021.
  * [11] B. Sutton et al., “Autonomous probabilistic coprocessing with petaflips per second,” _IEEE Access_ , vol. 8, pp. 157 238–157 252, 2020.
  * [12] G. Finocchio et al., “The promise of spintronics for unconventional computing,” _Journal of Magnetism and Magnetic Materials_ , vol. 521, p. 167506, 2021.
  * [13] S. Bhatti et al., “Spintronics based random access memory: a review,” _Materials Today_ , vol. 20, no. 9, pp. 530–548, 2017.
  * [14] A. Z. Pervaiz et al., “Weighted p-bits for fpga implementation of probabilistic circuits,” _IEEE transactions on neural networks and learning systems_ , 2018.
  * [15] S. C. Smithson et al., “Efficient cmos invertible logic using stochastic computing,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 66, no. 6, pp. 2263–2274, 2019.
  * [16] R. Rahman and S. Bandyopadhyay, “Variability of binary stochastic neurons employing low energy barrier nanomagnets with in-plane anisotropy,” _arXiv preprint arXiv:2108.04319_ , 2021.
  * [17] I. Hen et al., “Probing for quantum speedup in spin-glass problems with planted solutions,” _Physical Review A_ , vol. 92, no. 4, p. 042325, 2015\.
  * [18] T. Albash and D. A. Lidar, “Demonstration of a scaling advantage for a quantum annealer over simulated annealing,” _Physical Review X_ , vol. 8, no. 3, p. 031016, 2018.
  * [19] J. Gonzalez et al., “Parallel gibbs sampling: From colored fields to thin junction trees,” in _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics_. JMLR Workshop and Conference Proceedings, 2011, pp. 324–332.