# Controlling Steering with Energy-Based Models

Mikita Balesni
Autonomous Driving Lab
University of Tartu
<EMAIL_ADDRESS>
&Ardi Tampuu
Autonomous Driving Lab
University of Tartu
<EMAIL_ADDRESS>
Tambet Matiisen
Autonomous Driving Lab
University of Tartu
<EMAIL_ADDRESS>

###### Abstract

So-called implicit behavioral cloning with energy-based models has shown
promising results in robotic manipulation tasks. We tested if the method’s
advantages carry on to controlling the steering of a real self-driving car
with an end-to-end driving model. We performed an extensive comparison of the
implicit behavioral cloning approach with explicit baseline approaches, all
sharing the same neural network backbone architecture. Baseline explicit
models were trained with regression (MAE) loss, classification loss (softmax
and cross-entropy on a discretization), or as mixture density networks (MDN).
While models using the energy-based formulation performed comparably to
baseline approaches in terms of safety driver interventions, they had a higher
whiteness measure, indicating higher jerk. To alleviate this, we show two
methods that can be used to improve the smoothness of steering. We confirmed
that energy-based models handle multimodalities slightly better than simple
regression, but this did not translate to significantly better driving
ability. We argue that the steering-only road-following task has too few
multimodalities to benefit from energy-based models. This shows that applying
implicit behavioral cloning to real-world tasks can be challenging, and
further investigation is needed to bring out the theoretical advantages of
energy-based models.

##  1 Introduction

Implicit behavioral cloning[14] with energy-based models [22] has shown a lot
of promise in robotic manipulation tasks. The theoretical advantages of
energy-based models include increased data efficiency and the ability to model
discontinuities and multimodalities in the output action distribution [14].
Here, we set out to evaluate energy-based models for controlling the steering
of a self-driving vehicle using end-to-end driving models [27, 5]. We work
with steering-only models as the usefulness of multimodality-handling is
evident in this output modality, as illustrated by the theoretical and
practical failure cases of unimodal models on Figure 1. Adding longitudinal
control would complicate the task and demand more training data while not
necessarily adding more multimodal situations.

Figure 1: Left: If the experts have passed the tree from left and right with
equal frequency in the training data, behavioral cloning with a unimodal
policy and no high-level navigation commands would average the training
trajectories and drive straight into the tree. Right: In practice, we have
experienced such behavior at locations where side roads enter the main road -
unimodal regression models tend to swerve slightly (the red trajectory)
towards the side road. These swerves are minor because side roads are rarely
taken in our training data, and keeping straight is the dominant behavior.

To validate the theoretical advantages empirically, we compare a simple
energy-based model (EBM) with several baseline approaches in a road-following
task in the real world. The main experiments were also repeated in the
VISTA[2] simulator. The explicit baseline models are based on the same neural
network architecture as the EBM, with only the necessary modifications. They
are trained with the MAE loss, classification loss (softmax with cross-
entropy), or as mixture density networks [4] (MDN). During on-policy testing,
the models controlled only the car’s steering; the location- and direction-
relevant velocity was taken from a previously recorded expert trajectory. The
evaluation was performed on a WRC 2022 Rally Estonia track designed to be
challenging for humans, which was not included in the training set.

According to the main evaluation metric, safety-driver intervention count,
energy-based models performed comparably to the baseline methods but had
noticeably jerkier steering. To alleviate this, we proposed two methods:
temporal smoothing of predicted steering angles and spatially-aware soft
targets for cross-entropy loss. Still, unimodal explicit behavioral cloning
performed best in terms of safety driver interventions and jerk. We argue that
the inductive bias enforced by unimodal losses makes them more data-efficient
in simpler road situations, but more data is needed to model situations
requiring true multimodalities.

The main contributions of the paper are as follows:

1. We show that controlling the steering of an autonomous vehicle with energy-based models in the real world performs comparably to the baseline explicit behavioral cloning approaches.

2. We propose two methods that effectively reduce the steering jerk of energy-based models: temporal smoothing of predicted steering angles and soft targets for the cross-entropy loss.

3. To our surprise, we find that energy-based models do not outperform any of the similar-architecture explicit behavioral cloning baseline approaches in the real-world road-following task and that representing multimodalities does not translate into better driving.

##  2 Background

#### End-to-end driving

End-to-end driving attempts to replace the classical modular self-driving
pipeline with a single neural network model [27]. In the purest form, an end-
to-end self-driving model takes in raw sensor data and yields actionable
commands such as steering angle, throttle, and brake values. Such models are
commonly trained with behavioral cloning [26, 5] to imitate human expert
commands in the same situation.

One of the popular end-to-end driving models is NVIDIA PilotNet [5]. In our
work, we use the PilotNet network architecture as the backbone of all models
because it is relatively fast to train and sufficient for the road-following
task we aimed for. We do not use conditioning on high-level commands [8] or
more complex network architectures [9] to keep the setup simple - our goal is
to compare similar-capacity implicit and explicit models, not to aim for the
best performance. We use steering angle as the network output, rather than
trajectory[6, 15] or costmap[29], to be in line with prior implicit learning
work in robotics [14] where the network predicted raw actuator signal.

Figure 2: Controlling steering with implicit behavioral cloning. Left:
explicit models output the predicted angle directly. Right: implicit energy-
based models return an energy value per observation-action pair. The action
yielding the lowest energy for the current observation is chosen via argmin.

A significant amount of end-to-end driving research is done in
simulations[10], which through repeatability and access to the state of the
world, allow benchmarking with more complex measures of driving quality [9,
1]. However, models created in simulations cannot be deployed in the real
world without adaptations [3]. Also, some problems related to steering delays
and passenger comfort only become apparent in the real world. Here, we have
sufficient data for learning the relatively simple road-following task, so we
choose to work in the real world.

While the data-driven approach to autonomous driving is viewed as the most
promising path to full autonomy by some authors [18, 16], significant problems
remain to be solved, most prominently in generalization, in the explainability
of decisions, and in providing safety guarantees [9, 24, 19].

#### Energy-based models

An energy function, described by Lecun et al. [22], is any continuous function
that measures "goodness" between two sets of variables, where "good" pairs
have a low energy value. Following Florence et al. [14] who coined the term,
we call behavioral cloning policies "implicit" when they are composed of
argmin and a continuous energy function $E$, such that:

$\mathbf{\hat{a}}=\operatorname*{arg\,min}_{\mathbf{a}\in\mathcal{A}}E(\mathbf{o},\mathbf{a}),$
(1)

where $\mathbf{o}$ is an observation, e.g., a camera image, and $a$ is an
action, e.g., a steering wheel angle. In the present work, $E$ is implemented
by a neural network with PilotNet architecture with minor modifications
(discussed below). The classical approach of a model directly computing the
action based on an observation is in this context called "explicit" behavioral
cloning (Figure 2).

Implicit behavioral cloning with energy-based models promises the following
three advantages compared to classical explicit behavioral cloning: the
ability to represent discontinuities sharply, the ability to represent
multimodal action distributions, and better generalization with improved data
efficiency. In this work, we mainly focus on the ability to model
multimodalities.

#### Evaluation of driving models

With behavioral cloning, models are optimized to make momentary decisions on
data originating from the distribution resulting from human driving. When
deployed, however, the solutions face a sequential decision-making task on
data originating from a distribution caused by their own driving. Off-policy
metrics computed on held-out datasets of expert driving, such as mean absolute
error (MAE), measure only the predictive ability of the models. However, such
measures are insufficient for predicting success at the sequential decision-
making task when deployed [7]. Despite modest correlations with driving
ability, we use these metrics for model selection, as is often done in related
works.

Among the on-policy metrics measured during model deployment, the number of
safety driver interventions, the mean distance between interventions (DBI),
and the amount of time or distance traveled autonomously are the most popular
metrics [27]. In our work, we chose the number of interventions as the main
metric, as the distance traveled was fixed.

Beyond just completing routes safely, the comfort of passengers matters. The
smoothness of driving has been related directly to passenger comfort and
perceived safety[17, 11]. We follow multiple previous works [12, 13, 28] that
quantify the smoothness of steering with whiteness, defined as:

$W=\sqrt{\frac{1}{D}\sum_{i=1}^{D}\Big{(}\frac{\delta P_{i}}{\delta
t}\Big{)}^{2}},$ (2)

where $\delta P_{i}$ is the steering angle change, $D$ is dataset size, and
$\delta t$ is the time between decisions.

##  3 Methods

Figure 3: Energy-based model training and inference procedures. Left: feeding
the observation and different action values (one action at a time) through the
energy-computing network results in a vector of energy values, which is then
optimized via CE loss to make the ground truth action have the lowest energy.
Right: Observation is fed through the network with different candidate action
values, producing an energy value per candidate; the lowest-energy action is
chosen.

#### Baseline models

We adapt the PilotNet architecture’s output layer to produce a variety of
baseline models. The regression baseline has just one output node optimized to
produce the steering value via MAE loss. The classification baseline has the
action space divided into N bins and optimizes predicting the right bin via
cross-entropy loss. Mixture density networks output between 1 to 5 triplets of
mean, standard deviation, and relevancy scores $\alpha_{i}$ ($\alpha_{i}\geq
0,\sum_{i}\alpha_{i}=1$), a linear combination of which produces a Gaussian
mixture model over action values. We use the mean of the most likely Gaussian
during deployment.

#### EBM Training and Inference

We adapt the training and inference algorithms from energy-based model
literature [14] with a few modifications (see Figure 3). As a first
modification, we use a constant grid of linearly-spaced steering angles durig
training and inference instead of sampling uniformly. Second, we do not use
inference-time optimization to improve the initial candidate actions. The
candidate action values cover the steering angle range densely, and further
optimization yielded no gains. Furthermore, a fixed set of values is required
by one of the proposed EBM modifications and helps to make a cleaner
comparison with the classification baseline. Early offline experiments (see
Appendix Figure 7) showed that these changes resulted in at least as good
performance on steering prediction as random sampling and inference-time
optimization.

This results in the following loss function:

$L_{EBM}=\text{CE}(\text{softmax}(\mathbf{-e}),\mathbf{y})=-\sum_{i=1}^{n+1}(y_{i}\cdot
log(-\frac{e^{e_{i}}}{\sum_{j=1}^{n+1}e^{e_{j}}})),$ (3)

where $\mathbf{e}$ is the vector of energy values produced by neural network
outputs $e_{i}=E(\mathbf{o},a_{i})$ and $\mathbf{y}$ is one-hot vector having
1 at the position of ground truth action. Both $\mathbf{e}$ and $\mathbf{y}$
contain $n+1$ elements: $n$ sampled values and one ground truth. To get the
action energy vector $\mathbf{e}$ for a single observation, the neural network
$E$ is run on a batch of samples with the same repeated observation
$\mathbf{o}$ and different actions $a_{i}$. In practice, the convolutional
part of the network runs on an image only once, with actions fused into the
model before the MLP head to reduce the time and memory usage.

Our initial implementation of EBM demonstrated a high whiteness score, i.e.,
high lateral jerk. This characteristic did not depend on the amount of
training data (see Appendix Figure 6). We explored two changes to the EBM to
combat this undesired characteristic.

#### EBM with Temporal Smoothing

If one could reduce a model’s sensitivity to slight differences in subsequent
camera frames, one would achieve temporally smoother predictions. An obvious
choice in the case of energy-based models is to minimize the difference
between predicted energy distributions at subsequent frames. So, we propose
adding a temporal smoothness loss term, defined as:

$L_{temp}=\alpha\Big{\|}\mathbf{e}_{t}-\mathbf{e}_{t+1}\Big{\|},$ (4)

where $\mathbf{e}_{t}$ stands for the vector of predicted energy values at
timestep $t$, and $\alpha$ is the smoothing strength. EBMs take actions as
input, so $\mathbf{e}_{t}$ and $\mathbf{e}_{t+1}$ have to be computed with the
same steering angle inputs, which motivated our use of a constant action grid
instead of random sampling. However, to stick with the conventional sampling,
one could also draw a random sample once per pair of frames. Since the ground
truth steering angle is often different for consecutive frames, its energy
values are masked in this loss term. A range of well-performing smoothing
strengths was found empirically. We use $\alpha=1.0$ for the temporally
smoothed EBM in the final experiments.

#### EBM with Soft Targets

We hypothesized that using one-hot targets in the cross-entropy loss is a
major contribution to the higher whiteness of EBMs. Forcing nearby steering
values to have drastically different energy is likely to make learning less
efficient and the energy landscape noisier. This can lead to higher variance
when choosing the best action via _argmin_.

Hence, we investigated a simple fix: use soft targets for the cross-entropy
loss. Whereas soft targets have been widely used in neural networks with the
purpose of regularization and better calibration [25, 23], our use case is a
bit different. Unlike usual classification targets, our outputs are ordinal,
and we aim to enforce spatial smoothness. We replace the one-hot ground-truth
vector with a vector assigning some of the probability to actions a few
degrees away from the ground truth (see Figure 4). Target probabilities are
computed as:

$\mathbf{p^{*}}=\text{softmax}\Big{(}\frac{-(\mathbf{a}-a_{GT})^{2}}{T}\Big{)},$
(5)

where $\mathbf{p^{*}}$ is a vector of target probabilities for cross-entropy
loss, $\mathbf{a}$ is the vector of input candidate steering values, $a_{GT}$
is the ground truth steering angle, and $T$ is the softmax temperature
($2.5*10^{-3}$ in all reported tests with soft targets). We picked the
temperature value such that $99.9\%$ of the probability mass was on $\pm 5$
degrees around the ground truth.

Figure 4: (a) One-hot (standard) cross-entropy pushes the probability of
negative samples down, even if they are indistinguishably close to the target.
(b) Soft targets give positive weights to the samples around the target
proportionally to the distance. Note that the horizontal axis is trimmed to
focus on the GT value; the full range of input actions we use is from
$\sim$-250 to $\sim$250 degrees.

##  4 Experimental Setup

#### Dataset and Training Pipeline

We use the training dataset by Tampuu et al. [28], which consists of 540 km of
human driving on WRC Rally Estonia tracks. These are usually very low-traffic
gravel roads. The recordings that make up the dataset are broken into training
(460 km) and evaluation (80 km), with the evaluation recordings used for off-
policy metric calculation and early stopping. The original dataset includes
camera and LiDAR images, but only camera frames are used in our experiments.
Image pre-processing and training details are specified in the supplementary
materials.

#### On-Policy Evaluation

The on-policy evaluation was performed on a 4.3 km section of the WRC Estonia
2022 SS10+14 Elva track 111https://www.rally-maps.com/Rally-Estonia-2022/Elva,
driven in both directions. No recordings from this track were in the training
set of the models. The speed was set to 80% of the speed a human driver used
in the same location and direction, extracted from a prior recording. In
practice, this meant a speed of up to 40 km/h. Setting the speed to 100% was
attempted but felt too dangerous with certain models. The testing was
performed over multiple weeks of September 2022.

The evaluation track is narrow, and driving off the road edge is hazardous for
the car, so the safety driver was free to intervene when they perceived
danger. An intervention was counted when the driver applied force to turn the
steering wheel. If the model turned the steering wheel simultaneously in the
same direction as the safety driver, it would not cause an intervention since
no force was applied.

Alongside the intervention count, we also report the whiteness of steering as
an on-policy metric. Here, _command whiteness_ $W_{cmd}$ stands for the
whiteness of the predicted steering commands during on-policy evaluation, and
_effective whiteness_ $W_{eff}$ refers to the resulting actual whiteness of
the front wheel angles as measured by the sensors. Command whiteness is
usually higher than effective whiteness due to the smoothing effect of the
real-car actuators. No matter the force (i.e., the angular acceleration of the
steering wheel), it takes time to reach the target value.

First, around twenty test runs were completed to select the best
representative of each model type over several hyperparameters and random
seeds. For the final experiments, the six most promising models were chosen:
an EBM with 512 linearly-spaced candidate values (standard, with temporal
smoothing, or with soft targets), a classification model with 512 bins, MDN
with 5 Gaussians, and a regression model with MAE loss. We performed four
evaluation runs per model across four days. The worst run for each model was
discarded to account for out-of-distribution weather (excessive sun or rain
drops on camera) or safety driver variance.

#### Evaluation in VISTA

Evaluating driving models in the real world can make it harder for other
researchers to replicate results. To aid reproducibility, we additionally run
the main experiments in the VISTA Driving Simulator[2]. VISTA is a data-driven
simulator that allows replaying recordings of real-world drives interactively
by reprojecting the viewpoint as desired. Thus, a simulator can be used for
on-policy, closed-loop evaluation, allowing fast and reproducible model
evaluation (as in standard model-based simulators) while staying visually
close to the real-world data distribution. To this end, we release our
evaluation code and the recording we used for evaluation in the VISTA
format.222https://github.com/UT-ADL/vista-evaluation We used a recording
produced by the strongest model completing the track without interventions.
This recording had the highest correlation with our results on all models (see
Appendix Table 3), due to being most in-distribution with the weather and
vegetation at the time of the real-world tests.

Absent a safety driver, we define crashes in VISTA as moments when the car
drives more than 2 meters away from the expert-driven trajectory. After a
crash, we restart the car two seconds further down the road. This evaluation
scheme has obvious limitations, for example, as the expert does not always
drive in the center of the road, and 2 meters would be too late or too soon to
cause the safety driver to disengage in reality. Yet, empirical results from
the evaluation in VISTA support the findings from our real-world experiments,
suggesting that VISTA can be used to reproduce our key results.

##  5 Results

Table 1: Generalization results, with three real-world and three virtual driving sessions per model. | Real world | VISTA
---|---|---
Model | Interventions | $W_{eff}$ | $W_{cmd}$ | Crashes | $W_{cmd}$
EBM | 4 | 35.25°/s | 176.93°/s | 2 | 114.33°/s
1 | 32.34°/s | 96.94°/s | 1 | 121.57°/s
2 | 28.57°/s | 223.59°/s | 2 | 121.67°/s
mean: | 2.33 | 32.05°/s | 165.82°/s | 1.67 | 119.19°/s
EBM Temp. Smoothing | 5 | 49.92°/s | 119.39°/s | 3 | 58.70°/s
2 | 38.96°/s | 137.22°/s | 2 | 60.37°/s
3 | 34.21°/s | 77.28°/s | 2 | 48.86°/s
mean: | 3.33 | 41.03°/s | 111.30°/s | 2.33 | 55.98°/s
EBM Soft Targets | 5 | 27.80°/s | 56.33°/s | 3 | 85.72°/s
5 | 46.83°/s | 57.15°/s | 3 | 74.97°/s
4 | 33.72°/s | 56.86°/s | 3 | 81.87°/s
mean: | 4.66 | 36.12°/s | 56.78°/s | 3 | 80.85°/s
Regression (MAE) | 2 | 26.69°/s | 37.84°/s | 0 | 24.39°/s
2 | 29.65°/s | 75.34°/s | 0 | 24.75°/s
1 | 26.28°/s | 33.10°/s | 0 | 24.25°/s
mean: | 1.66 | 27.54°/s | 48.76°/s | 0 | 24.47°/s
Classification | 1 | 41.05°/s | 182.39°/s | 1 | 123.69°/s
7 | 62.17°/s | 287.14°/s | 1 | 105.13°/s
1 | 34.11°/s | 162.27°/s | 1 | 104.31°/s
mean: | 3.00 | 45.77°/s | 210.60°/s | 1 | 111.04°/s
MDN | 1 | 25.32°/s | 33.62°/s | 3 | 37.22°/s
5 | 24.82°/s | 35.46°/s | 3 | 35.74°/s
5 | 26.66°/s | 37.39°/s | 3 | 35.84°/s
mean: | 3.66 | 25.59°/s | 35.49°/s | 3 | 36.27°/s

The results of the final test runs are presented in Table 4. Each row in the
table corresponds to approximately 20 minutes of driving, so each model’s
total intervention count is produced by an hour of driving. We believe this is
enough time to reveal noteworthy differences in performance on a simple task
such as road following. Real-world results are supported by VISTA showing
similar trends with lower variance (Pearson R = 87.5%, Spearman R = 92.7% for
real-world interventions in final experiments vs VISTA crashes).

We observe no clear benefit of using energy-based implicit behavioral cloning
models over the explicit baseline models. The baseline regression model
resulted in the least interventions and had also the most smooth driving
according to the whiteness measure. The other two explicit baseline models:
classification and mixture density networks, resulted in more interventions
across the three runs than the EBM. Supported by VISTA, the results reveal
that no one solution stands out clearly, and hence that using the EBM
formulation did not significantly improve the performance.

To bring the whiteness values of EBMs closer to the values of explicit models,
we implemented two variations to the EBM. Both approaches significantly
reduced the whiteness of the model predictions during deployment ($W_{cmd}$ in
Table 4). However, this reduction did not translate into a reduction of
effective whiteness, while both of these approaches resulted in more safety
driver interventions.

The increased temporal stability of modified EBM approaches compared to the
naive EBM is also visible in Figure 5. This illustration visualizes the
outputs of different models computed off-policy on a recording of the vehicle
passing an intersection. Furthermore, we see the classifier model would swerve
to the left if given control of the car near the intersection.

Figure 5: Outputs of different models during 5 seconds of human driving (a-f).
The period corresponds to the car passing an intersection, with Y-axis
representing the time and X-axis the predicted steering angle. The camera
image 2 seconds into the recording is given on the right (g). The red dashed
line represents the predicted steering values, grey dashed line is the ground
truth.

EBM models exhibit slight multimodality when passing the intersection with an
area of lower energy values towards the left, corresponding to a left turn
(Figure 5 (d-f)). The standard EBM also shows an alternative hypothesis of
turning to the right around 3 seconds into the recording. As reported above,
however, this ability to represent alternative hypotheses in the energy
landscape does not translate into improved performance.

Table 2: Handling multimodalities. Swerve rate shows the percentage of
challenging road sections where the model slightly swerved towards the side
road. There are three such road sections, with three on-policy runs per model
type. Slight swerves are difficult to set a threshold for and are counted as
half a swerve. A lower number is better.
Model | Swerve rate
---|---
EBM | 44%
EBM Temp. Smoothing | 66%
EBM Soft Targets | 39%
Regression | 89%
Classification | 44%
MDN | 33%

Supposedly, a model having an explicit representation of alternative paths can
select the most likely and not produce intermediate behavior that is average
of multiple options. At intersections, such average behavior would show as
swerving towards the side road. When quantified (see Table 2), classical
regression models were shown to produce most swerving towards side roads. This
observation aligns with the unimodal nature of regression models trained with
the MAE loss. In contrast, other models swerved less, which can be attributed
to their richer representations.

##  6 Discussion

This project investigated if the reported benefits of using the energy-based
model formulation for behavioral cloning carry over to the task of real-world
road following. We hypothesized that the claimed better generalization and
handling of multimodalities could be useful in this task. However, the results
show no improvements in overall driving ability. We believe the chosen task
has too few multimodalities to make EBMs stand out. Prior work on implicit
behavior cloning[14] used tasks where actions are distributed less normally
than in road-following steering control (see Appendix Figure 8). Past tasks
also had at least two-dimensional action spaces, which may have more
multimodalities. However, our models only predicted steering, and representing
different steering hypotheses was useful only at a few intersections along the
route.

We did observe some improvement in handling situations that presumably require
modeling multimodal distributions, such as intersections and ignoring side
roads. The only unimodal baseline, MAE-based regression, swerved towards side
roads more frequently than other models. However, less frequent swerving of
multimodal models did not result in fewer interventions overall. We attribute
this to the low proportion of the task requiring a multimodal policy.
Conversely, the unimodal loss seems to introduce an inductive bias, increasing
the data efficiency in learning to handle simpler (unimodal) road situations
that dominate the task.

We observed higher lateral jerk for multimodality-representing models. Our
proposed modifications to the EBM training process significantly reduced the
jerk of the model-predicted commands. However, this did not improve the
effective whiteness of the car’s front wheels. We attribute this to the car’s
actuators acting as a low-pass filter on the noise in the command sequence.
This prevents even a large change in command whiteness from translating to a
drop in front-wheel whiteness, motivating research on more powerful smoothing
techniques. The temporal loss term was computed on consecutive images only 33
ms apart. However, effective whiteness seems to be caused by output
variability on a slightly higher time scale. In future work, smoothing actions
across a slightly longer timescale should be attempted. Soft targets proved
surprisingly effective in reducing command-sequence whiteness, given that they
do not directly enforce similarity across time.

##  7 Conclusion

We tested implicit behavioral cloning with energy-based models for controlling
the steering of a real self-driving car. We showed that energy-based models
perform comparably to classical explicit behavioral cloning baselines in terms
of safety driver interventions but have higher jerk that reduces the comfort
of the drive. We show two methods for reducing the steering jerk, measured as
a whiteness score. Even though these methods greatly reduce the whiteness of
predicted steering angles, it does not translate into improved whiteness of
real steering, as the actuator delays in a real car smooth out radical
steering movements anyway.

In our experiments, the simple regression-based explicit behavioral cloning
baseline was the best in terms of interventions and jerk of the drive.
However, the regression approach tended to swerve towards side roads, which
comes from the unimodal nature of its loss function. We show that
multimodality-capable models handle the situation with side roads better and
do fewer swerves but do not eliminate the problem. Based on the analysis of
action distributions in our task and in prior work where EBMs outperformed
explicit models, we conclude that the lateral control in the road-following
task has too few multimodalities to make EBMs useful. Altogether, this shows
that while energy-based models have a number of theoretical advantages, it can
be challenging to bring those out in real-world scenarios, and more research
is needed to make efficient use of them.

## Acknowledgments and Disclosure of Funding

This work was supported by the collaboration project LLTAT21278 with Bolt
Technologies, and by the Estonian Research Council grant PRG1604. Mikita
Balesni was funded by the Estonian Research Council grant for Ukrainian
Researchers.

## References

  * car [2022] Carla leaderboard - evaluation and metrics, 2022.  URL https://leaderboard.carla.org/#evaluation-and-metrics.
  * Amini et al. [2022] A. Amini, T.-H. Wang, I. Gilitschenski, W. Schwarting, Z. Liu, S. Han, S. Karaman, and D. Rus.  Vista 2.0: An open, data-driven simulator for multimodal sensing and policy learning for autonomous vehicles.  In _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022.
  * Bewley et al. [2019] A. Bewley, J. Rigley, Y. Liu, J. Hawke, R. Shen, V.-D. Lam, and A. Kendall.  Learning to drive from simulation without real world labels.  In _2019 International conference on robotics and automation (ICRA)_ , pages 4818–4824. IEEE, 2019.
  * Bishop [1994] C. M. Bishop.  Mixture density networks.  1994\.
  * Bojarski et al. [2016] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, et al.  End to end learning for self-driving cars.  _arXiv preprint arXiv:1604.07316_ , 2016.
  * Bojarski et al. [2020] M. Bojarski, C. Chen, J. Daw, A. Değirmenci, J. Deri, B. Firner, B. Flepp, S. Gogri, J. Hong, L. Jackel, et al.  The nvidia pilotnet experiments.  _arXiv preprint arXiv:2010.08776_ , 2020.
  * Codevilla et al. [2018a] F. Codevilla, A. M. López, V. Koltun, and A. Dosovitskiy.  On offline evaluation of vision-based driving models.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 236–251, 2018a.
  * Codevilla et al. [2018b] F. Codevilla, M. Miiller, A. López, V. Koltun, and A. Dosovitskiy.  End-to-end driving via conditional imitation learning.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1–9. IEEE, 2018b.
  * Codevilla et al. [2019] F. Codevilla, E. Santana, A. M. López, and A. Gaidon.  Exploring the limitations of behavior cloning for autonomous driving.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 9329–9338, 2019.
  * Dosovitskiy et al. [2017] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun.  CARLA: An open urban driving simulator.  _arXiv preprint arXiv:1711.03938_ , 2017.
  * Elbanhawi et al. [2015] M. Elbanhawi, M. Simic, and R. Jazar.  In the passenger seat: investigating ride comfort measures in autonomous cars.  _IEEE Intelligent transportation systems magazine_ , 7(3):4–17, 2015.
  * Eraqi et al. [2017] H. M. Eraqi, M. N. Moustafa, and J. Honer.  End-to-end deep learning for steering autonomous vehicles considering temporal dependencies, 2017.  URL https://arxiv.org/abs/1710.03804.
  * Fernandez [2018] N. Fernandez.  Two-stream convolutional networks for end-to-end learning of self-driving cars.  2018\.  doi: 10.48550/ARXIV.1811.05785.  URL https://arxiv.org/abs/1811.05785.
  * Florence et al. [2021] P. Florence, C. Lynch, A. Zeng, O. Ramirez, A. Wahid, L. Downs, A. Wong, J. Lee, I. Mordatch, and J. Tompson.  Implicit behavioral cloning.  _Conference on Robot Learning (CoRL)_ , 2021.
  * Hawke et al. [2020] J. Hawke, R. Shen, C. Gurau, S. Sharma, D. Reda, N. Nikolov, P. Mazur, S. Micklethwaite, N. Griffiths, A. Shah, et al.  Urban driving with conditional imitation learning.  In _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 251–257. IEEE, 2020.
  * Hawke et al. [2021] J. Hawke, V. Badrinarayanan, A. Kendall, et al.  Reimagining an autonomous vehicle.  _arXiv preprint arXiv:2108.05805_ , 2021.
  * Hecker et al. [2020] S. Hecker, D. Dai, A. Liniger, M. Hahner, and L. Van Gool.  Learning accurate and human-like driving using semantic maps and attention.  In _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 2346–2353. IEEE, 2020.
  * Jain et al. [2021] A. Jain, L. Del Pero, H. Grimmett, and P. Ondruska.  Autonomy 2.0: Why is self-driving always 5 years away?  _arXiv preprint arXiv:2107.08142_ , 2021.
  * Kalra and Paddock [2016] N. Kalra and S. M. Paddock.  Driving to safety: How many miles of driving would it take to demonstrate autonomous vehicle reliability?  _Transportation Research Part A: Policy and Practice_ , 94:182–193, 2016.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2014.  URL https://arxiv.org/abs/1412.6980.
  * Krogh and Hertz [1991] A. Krogh and J. A. Hertz.  A simple weight decay can improve generalization.  In _Proceedings of the 4th International Conference on Neural Information Processing Systems_ , NIPS’91, page 950–957, San Francisco, CA, USA, 1991. Morgan Kaufmann Publishers Inc.  ISBN 1558602224.
  * Lecun et al. [2006] Y. Lecun, S. Chopra, and R. Hadsell.  _A tutorial on energy-based learning_.  01 2006.
  * Müller et al. [2019] R. Müller, S. Kornblith, and G. E. Hinton.  When does label smoothing help?  _Advances in neural information processing systems_ , 32, 2019.
  * Nassi et al. [2020] B. Nassi, D. Nassi, R. Ben-Netanel, Y. Mirsky, O. Drokin, and Y. Elovici.  Phantom of the adas: Phantom attacks on driver-assistance systems.  https://eprint.iacr.org/2020/085.pdf, 2020.
  * Pereyra et al. [2017] G. Pereyra, G. Tucker, J. Chorowski, Ł. Kaiser, and G. Hinton.  Regularizing neural networks by penalizing confident output distributions.  _arXiv preprint arXiv:1701.06548_ , 2017.
  * Pomerleau [1989] D. A. Pomerleau.  Alvinn: An autonomous land vehicle in a neural network.  In _Advances in neural information processing systems_ , pages 305–313, 1989.
  * Tampuu et al. [2020] A. Tampuu, T. Matiisen, M. Semikin, D. Fishman, and N. Muhammad.  A survey of end-to-end driving: Architectures and training methods.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Tampuu et al. [2022] A. Tampuu, R. Aidla, J. A. van Gent, and T. Matiisen.  Lidar-as-camera for end-to-end driving, 2022.  URL https://arxiv.org/abs/2206.15170.
  * Zeng et al. [2019] W. Zeng, W. Luo, S. Suo, A. Sadat, B. Yang, S. Casas, and R. Urtasun.  End-to-end interpretable neural motion planner.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8660–8669, 2019.

##  Appendix A Appendix

###  A.1 Pre-processing, training, and hardware

#### Pre-processing

The frames are cropped to remove the car’s hood and everything beyond the
horizon and to limit the view to 90 degrees of the front center. The resulting
frames of shape 264x68x3 are then min-max normalized and fed into the model.

#### Training

When training, a mini-batch is created by sampling uniformly from all
recordings. For experiments with temporal smoothing, a different sampling
approach is used, where sequences of two consecutive frames are sampled
instead. The sequence dimension is flattened such that a mini-batch of
sequences becomes a mini-batch of frames. The target labels correspond to the
steering wheel angles of the human drivers.

The Adam [20] optimizer is used with default hyperparameters (learning rate
$1*10^{-3}$, betas $0.9$ and $0.999$) and $1*10^{-2}$ weight decay [21].
Finally, early stopping is used on validation MAE with a patience of 10
epochs.

#### Car Hardware and Software Stack

We perform the experiments with Lexus RX 450h fitted with a PACMod v3 drive-
by-wire system. The following sensors are used: a NovAtel PwrPak7D-E2 GNSS
device and a Sekonix SF3324 120-degree FOV camera. The car computer is
equipped with a GeForce GTX 2080 GPU. The camera works at 30 Hz, but our end-
to-end stack is slower ($\sim$12 Hz). To accommodate for the differing
processing speeds all but the latest frame in the queue are dropped.

###  A.2 Validation loss and whiteness with different data amounts

Figure 6: Left: Varying the dataset size intuitively changes the predictive
accuracy of all model formulations. Right: In contrast, whiteness does not
seem to be influenced by the amount of data when the amount was varied by 5
times.

###  A.3 How do our EBM simplifications affect modeling performance?

Figure 7: Left: Using a constant grid of actions results in at least as good
MAE as the more common uniformly random sample on each decision. Right:
Derivative-free optimization does not improve performance beyond a one-shot
argmin; to our surprise, it even hurts validation MAE.

###  A.4 Is road following just too unimodal?

(a)

(b)

(c)

(d)

Figure 8: Action distributions from three randomly-picked example tasks from
prior work on IBC[14] where EBMs outperformed explicit models (a-c). Action
distribution for our road-following task (d) looks much more gaussian, which
is a hint for a lower number of possible multimodalities.

###  A.5 VISTA agreement with real-world results

In the main text, we report VISTA’s agreement with the results in the main
experiments. In Table 3, we report the results for all models we tested
throughout the project. When run on the same track, VISTA seems to have a very
high agreement with real-world results.

Table 3: Per-model mean metrics correlations for VISTA vs reality (n=17 models). Measure | Interventions | $W_{cmd}$
---|---|---
Pearson | 83% | 89%
Spearman | 84% | 86%

###  A.6 Notepad

Table 4: Generalization results, with three real-world and three virtual driving sessions per model. | Real world | VISTA
---|---|---
Model | Interventions | Whiteness | Crashes | Whiteness
EBM | 4 | 176.93°/s | 2 | 114.33°/s
1 | 96.94°/s | 1 | 121.57°/s
2 | 223.59°/s | 2 | 121.67°/s
mean: | 2.33 | 165.82°/s | 1.67 | 119.19°/s
EBM Temp. Smoothing | 5 | 119.39°/s | 3 | 58.70°/s
2 | 137.22°/s | 2 | 60.37°/s
3 | 77.28°/s | 2 | 48.86°/s
mean: | 3.33 | 111.30°/s | 2.33 | 55.98°/s
EBM Soft Targets | 5 | 56.33°/s | 3 | 85.72°/s
5 | 57.15°/s | 3 | 74.97°/s
4 | 56.86°/s | 3 | 81.87°/s
mean: | 4.66 | 56.78°/s | 3 | 80.85°/s
Regression (MAE) | 2 | 37.84°/s | 0 | 24.39°/s
2 | 75.34°/s | 0 | 24.75°/s
1 | 33.10°/s | 0 | 24.25°/s
mean: | 1.66 | 48.76°/s | 0 | 24.47°/s
Classification | 1 | 182.39°/s | 1 | 123.69°/s
7 | 287.14°/s | 1 | 105.13°/s
1 | 162.27°/s | 1 | 104.31°/s
mean: | 3.00 | 210.60°/s | 1 | 111.04°/s
MDN | 1 | 33.62°/s | 3 | 37.22°/s
5 | 35.46°/s | 3 | 35.74°/s
5 | 37.39°/s | 3 | 35.84°/s
mean: | 3.66 | 35.49°/s | 3 | 36.27°/s


# Fast evaluation of spherical harmonics with sphericart

Filippo Bigi  Guillaume Fraux  Laboratory of Computational Science and
Modelling, Institute of Materials, École Polytechnique Fédérale de Lausanne,
Lausanne 1015, Switzerland  National Centre for Computational Design and
Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de
Lausanne, 1015 Lausanne, Switzerland  Nicholas J. Browning  Swiss National
Supercomputing Centre (CSCS), 6900, Lugano, Switzerland  Michele Ceriotti
Laboratory of Computational Science and Modelling, Institute of Materials,
École Polytechnique Fédérale de Lausanne, Lausanne 1015, Switzerland  National
Centre for Computational Design and Discovery of Novel Materials (MARVEL),
École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
<EMAIL_ADDRESS>
###### Abstract

Spherical harmonics provide a smooth, orthogonal, and symmetry-adapted basis
to expand functions on a sphere, and they are used routinely in physical and
theoretical chemistry as well as in different fields of science and
technology, from geology and atmospheric sciences to signal processing and
computer graphics. More recently, they have become a key component of
rotationally equivariant models in geometric machine learning, including
applications to atomic-scale modeling of molecules and materials. We present
an elegant and efficient algorithm for the evaluation of the real-valued
spherical harmonics. Our construction features many of the desirable
properties of existing schemes and allows to compute Cartesian derivatives in
a numerically stable and computationally efficient manner. To facilitate
usage, we implement this algorithm in sphericart, a fast C++ library which
also provides C bindings, a Python API, and a PyTorch implementation that
includes a GPU kernel.

##  I Introduction

The spherical harmonics $Y^{m}_{l}$ are basis functions for the irreducible
representations of the SO(3) group,1 which makes them a key tool in
understanding physical phenomena that exhibit rotational symmetries, and to
design practical algorithms to model them on a computer. Examples include the
distribution of charge in atoms2, the behavior of gravitational3 and
magnetic4; 5 fields, and the propagation of light6 and sound7 in the
atmosphere. Moreover, they provide a complete set of smooth, orthogonal
functions defined on the surface of a sphere, and in this sense they are
widely used in many fields including computer graphics8; 9, quantum 10; 11;
12; 13 and physical 14; 15; 16 chemistry, and signal processing17; 18. More
recently, they have become an essential tool in the context of geometric deep
learning19; 20, as a structural descriptor needed in SO(3)-, O(3)-, and
E(3)-equivariant machine learning models21; 22; 23; 24 and more specifically
in the construction of symmetry-adapted descriptors of atomic structures in
chemical machine learning25; 26; 27; 28; 29. Derivatives of the spherical
harmonics are also very often used in these applications. An example is that
of the calculation of forces on a configuration of atoms by an E(3)-invariant
machine learning model.30 Likewise, spherical harmonics gradients have been
used in computer graphics for mid-range illumination, irradiance volumes, and
optimization methods.9

In many of these contexts, the spherical harmonics are used together with an
expansion in the radial direction, to compute expressions of the form

$c_{nlm}=\sum_{k}R_{nl}(r_{k})Y^{m}_{l}(\hat{\mathbf{r}}_{k}),$ (1)

where $R_{nl}$ indicates a radial basis, and
$\mathbf{r}_{k}=r_{k}\hat{\mathbf{r}}_{k}$ a set of points that correspond to
either particles, or to the positions of a Cartesian mesh. In others, such as
message-passing neural networks,31 directional terms that depend on the
orientation of interatomic separation vectors are multiplied by continuous
filters that are a function of the radial distance. For this reason, although
the spherical harmonics are most often defined as complex functions in
spherical coordinates, in practical implementations it is often preferred to
use their real-valued combinations, and to express their value and derivatives
directly in terms of the Cartesian coordinates of the points at which they are
evaluated. With these applications in mind, we derive compact, efficient
expressions to compute the real spherical harmonics and their derivatives of
arbitrary order as polynomials of the Cartesian coordinates, we discuss the
computational implications of this formulation, and we present a simple yet
efficient implementation that can be used both as a C and C++ library and as a
Python module.

##  II Analytical expressions

The real-valued spherical harmonics can be defined in spherical coordinates
$(\theta,\phi)$ as

$Y_{l}^{m}(\theta,\phi)=F_{l}^{|m|}*P_{l}^{|m|}(\cos{\theta})*\begin{cases}\sin{(|m|\phi)}&\text{if
$m<0$}\\\ 1/\sqrt{2}&\text{if $m=0$}\\\ \cos{(m\phi)}&\text{if
$m>0$}\end{cases}$ (2)

where $P_{l}^{m}$ is an associated Legendre polynomial and $F_{l}^{m}$ is a
prefactor which takes the form

$F_{l}^{m}=(-1)^{m}\sqrt{\frac{2l+1}{2\pi}\frac{(l-m)!}{(l+m)!}}.$ (3)

Possible strategies for a stable and efficient computation of the $F_{l}^{m}$
are discussed in Appendix C.

Traditionally32, the real spherical harmonics of points in 3D space have been
calculated by first converting the Cartesian coordinates $x$, $y$, and $z$
into the spherical coordinates $\theta$ and $\phi$, and then using Eq. 2 to
calculate the spherical harmonics, most often via the standard recurrence
relations for $P_{l}^{m}(t)$:

$\begin{split}&P_{0}^{0}=1,\quad
P_{l}^{l}=-\,(2l-1)\,\sqrt{1-t^{2}}\,P_{l-1}^{l-1},\\\
&P_{l}^{l-1}=(2l-1)\,t\,P_{l-1}^{l-1},\\\
&P_{l}^{m}=[(2l-1)\,t\,P_{l-1}^{m}-(l+m-1)\,P_{l-2}^{m}]/(l-m).\end{split}$
(4)

It should be noted how naive differentiation of these recursions introduces
poles on the z axis where $t=\cos{\theta}=\pm 1$, as well as numerical
instabilities in the calculation of the derivatives for points close to the z
axis.

In contrast, our algorithm aims at calculating
$\tilde{Y}_{l}^{m}=r^{l}Y_{l}^{m}$, where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. These
are the so-called _solid_ harmonics, which consist of simple homogeneous
polynomials of the Cartesian coordinates. This choice avoids the need to
normalize $r$, and it leads to simple and stable iterations to compute
$\tilde{Y}_{l}^{m}$, as well as very compact expressions for their derivatives
with respect to the Cartesian coordinates which re-use the same factors needed
to evaluate $\tilde{Y}_{l}^{m}$. In most applications that require spherical
harmonics in Cartesian coordinates, the radial direction is dealt with by a
separate expansion (cf. Eq. (1)), and the $r^{l}$ factor that is included in
the scaled $\tilde{Y}_{l}^{m}$ can be compensated for at little to no
additional cost by a corresponding $r^{-l}$ factor in the radial term. If,
instead, the spherical harmonics are needed in their conventional version,
they can be recovered easily from the radially scaled (or solid) kind.

As shown in Appendix A, the scaled Cartesian harmonics $\tilde{Y}_{l}^{m}$ can
be computed as

$\tilde{Y}_{l}^{m}(x,y,z)=F_{l}^{|m|}Q_{l}^{|m|}(z,r)\times\begin{cases}s_{|m|}(x,y)&\text{if
$m<0$}\\\ 1/\sqrt{2}&\text{if $m=0$}\\\ c_{m}(x,y)&\text{if $m>0$}\end{cases}$
(5)

where we define

$\begin{split}Q_{l}^{m}=&r^{l}r_{xy}^{-m}\,P_{l}^{m},\\\
s_{m}=&r_{xy}^{m}\,\sin{(m\phi)},\\\
c_{m}=&r_{xy}^{m}\,\cos{(m\phi)},\end{split}$ (6)

with $r_{xy}=\sqrt{x^{2}+y^{2}}$. Similar (but not equivalent) definitions
have been used often in the literature, e.g. in Refs. 33; 14; 34. The
quantities $Q_{l}^{m}$, $c_{m}$, and $s_{m}$ can be evaluated very efficiently
by recursion in Cartesian coordinates. For example, the recursions for
$Q_{l}^{m}$ follow almost immediately (see Appendix A) from those for
$P_{l}^{m}$ (Eq. (4)) and the definition of $Q_{l}^{m}$ (Eq. (6)):

$\begin{split}&Q_{0}^{0}=1,\quad Q_{l}^{l}=-\,(2l-1)\,Q_{l-1}^{l-1},\\\
&Q_{l}^{l-1}=(2l-1)\,z\,Q_{l-1}^{l-1}=-zQ_{l}^{l},\\\
&Q_{l}^{m}=[(2l-1)\,z\,Q_{l-1}^{m}-(l+m-1)\,r^{2}\,Q_{l-2}^{m}]/(l-m).\end{split}$
(7)

Many other recursive expressions can be derived based on analogous, well-known
relations for $P_{l}^{m}$, e.g.

$\begin{split}Q_{l}^{m}=\frac{2(m+1)\,z\,Q_{l}^{m+1}+r_{xy}^{2}\,Q_{l}^{m+2}}{(l+m+1)(l-m)}.\end{split}$
(8)

which can be used to iterate _down_ from $Q_{l}^{l}$, avoiding the pole at
$r_{xy}=0$ that is present for the similar recursion for $P_{l}^{m}$.

Similarly, recursive relations for $s_{m}$ and $c_{m}$ can be derived (see
Appendix A) as

$\begin{split}&s_{0}=0,\quad c_{0}=1,\\\ &s_{m}=xs_{m-1}+yc_{m-1},\quad
c_{m}=-ys_{m-1}+xc_{m-1}.\end{split}$ (9)

Once the $Q_{l}^{m}$, $c_{m}$, and $s_{m}$ quantities are known, their
Cartesian derivatives also follow in an extremely compact form:

$\begin{split}&\frac{\partial Q_{l}^{m}}{\partial x}\\!=\\!xQ_{l-1}^{m+1},\
\frac{\partial Q_{l}^{m}}{\partial y}\\!=\\!yQ_{l-1}^{m+1},\ \frac{\partial
Q_{l}^{m}}{\partial z}\\!=\\!(l+m)Q_{l-1}^{m},\\\ &\frac{\partial
s_{m}}{\partial x}=m\,s_{m-1},\ \frac{\partial s_{m}}{\partial y}=m\,c_{m-1},\
\\\ &\frac{\partial c_{m}}{\partial x}=m\,c_{m-1},\ \frac{\partial
c_{m}}{\partial y}=-\,m\,s_{m-1}.\end{split}$ (10)

These relationships (proven in Appendix B) are much simpler than the standard
recurrence relations for the derivatives of $P_{l}^{m}$ (see e.g. Ref. 35),
and they do not lead to poles or instabilities when combined to compute the
derivatives of $\tilde{Y}^{m}_{l}$. The simplicity of the expressions in Eq.
10 opens the door to the efficient calculation of Cartesian derivatives of
$\tilde{Y}^{m}_{l}$ of any order. In addition, one can find several
expressions that directly link _some_ of the derivatives of the Cartesian
spherical harmonics to other $(l,m)$ values, e.g.

$\begin{split}\frac{\partial\tilde{Y}^{m}_{l}}{\partial
z}=({F_{l}^{|m|}}/{F_{l-1}^{|m|}})(l+m)\tilde{Y}_{m}^{l-1},\\\
\frac{\partial\tilde{Y}_{l}^{l}}{\partial
x}=-l(2l-1)({F_{l}^{l}}/{F_{l-1}^{l-1}})\tilde{Y}_{l-1}^{l-1}\end{split}$ (11)

that simplify even further the calculation of the derivatives of the spherical
harmonics. For $l\leq 6$, we use a computer algebra system to automatically
find the expressions that provide $\tilde{Y}^{m}_{l}$ and their derivatives in
terms of lower-$l$ values with the smallest number of multiplications, and we
use them to generate hard-coded implementations, as discussed in the following
Section.

##  III Computer implementation and benchmarking

Most computational applications require the evaluation of all the spherical
harmonics, and possibly their derivatives, for the values of $l$ up to a
maximum degree $l_{\mathrm{max}}$. In many cases, the spherical harmonics have
to be computed for many points simultaneously, e.g. for the interatomic
separation vectors of all neighbors of a selected atom.

Even though we recommend to use the scaled spherical harmonics in
applications, accounting for the fact that accompanying radial terms have to
compensate for the scaling, we also provide an implementation of “normalized”
spherical harmonics, by simply evaluating $\tilde{Y}_{l}^{m}(x/r,y/r,z/r)$ and
applying the chain rule to the derivatives. These additional operations
typically result in an overhead of 5-10% for the calculation of the spherical
harmonics and their derivatives.

We implement a general routine that takes a list of Cartesian coordinates and
evaluates the scaled $\tilde{Y}_{l}^{m}$ and, optionally, their derivatives.
We use C++ for the implementation, using templates to exploit compile-time
knowledge of the maximum angular momentum, the need to compute derivatives,
etc. We then define a pure C API covering the typical use cases on top of this
C++ API. Since almost all programming languages have a way to call C
functions, this enables using our code from most languages used in a
scientific context. We also provide a Python package with a high-level
interface to the C library.

Furthermore, we use the C++ API to provide an implementation compatible with
PyTorch36, using a custom backward function to compute gradients using the
derivatives evaluated in the forward pass, and making sure the code is
compatible with TorchScript, allowing to use models without a Python
interpreter. Finally, we implement a GPU-accelerated version of the PyTorch
code for NVIDIA GPUs, using the CUDA language.

###  III.1 CPU implementation details

We apply a number of trivial (and a few less obvious) optimizations. For
example, we pre-compute the factors $F_{l}^{m}$ (Eq. (3)) to minimize the
number of operations in the inner loop of the iterative algorithm discussed
above. In order to further accelerate the evaluation of low-$l$ Cartesian
harmonics, we use a computer algebra system to derive expressions that
evaluate the full $\tilde{Y}^{m}_{l}$ and their derivatives with hard-coded
expressions using the smallest possible number of multiplications. As shown in
Table 1, and as noted in previous implementations33; 24, there is considerable
scope for optimization by using ad hoc expressions. However, the speedup is
less remarkable when including the calculation of the derivatives through Eq.
(10), which re-uses quantities that have already been computed when evaluating
$\tilde{Y}_{l}^{m}$.

| general-purpose | hard-coded
---|---|---
|  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$
$l_{\text{max}}=1$ |  3.49 |  11.0 |  1.33 |  7.85
$l_{\text{max}}=2$ |  7.16 |  21.6 |  4.21 |  16.4
$l_{\text{max}}=3$ |  13.3 |  36.1 |  7.26 |  28.4
$l_{\text{max}}=4$ |  20.3 |  55.3 |  11.8 |  46.3
$l_{\text{max}}=5$ |  29.8 |  81.9 |  16.3 |  68.3
$l_{\text{max}}=6$ |  42.8 |  121 |  22.2 |  95.2
Table 1: Serial execution time (in ns/point) for computing Cartesian spherical
harmonics and their derivatives up to the indicated value of $l_{\text{max}}$
in double precision on an Intel Xeon Gold 6226R CPU, averaged over 1000 calls
for 10 000 points, comparing our general purpose recursive algorithm to an
hard-coded implementation. Figure 1: Scaling of wall-clock timing to evaluate
the Cartesian spherical harmonics and their derivatives for 10 000 points
using different $l_{\text{max}}$ and numbers of OpenMP threads. Results refer
to up to 64 cores of up to two Intel Xeon Platinum 8360Y (2.4 GHz) CPUs.

As a compromise between the convenience of a function that works for arbitrary
$l_{\text{max}}$ and the efficiency of optimized expressions, we provide an
interface for the hard-coded implementation up to $l_{\text{max}}=6$, and use
by default a hybrid implementation that applies the hard-coded version for
small $l$, and then switches to the general expression, using the alternative
recursion (8) to avoid computing low-$l$ values of the modified Legendre
polynomials $Q_{l}^{m}$. All our functions can be applied to many 3D points at
once, and they are trivially parallelized over the sample index using
OpenMP.37 Despite the simplicity of our design, we achieve good parallel
scaling, particularly when using a relatively-large $l_{\text{max}}$ so that
there is a substantial amount of computation for each thread, as illustrated
in Fig. 1. The curves show some degree of irregular behavior, suggesting that
there might be room for further optimization. As shown in Table 2, for large
$l_{\text{max}}$ and numbers of data points, there is a substantial
performance gain by using single-precision floating-point arithmetics.
However, for smaller amounts of computation there is very small advantage, and
in some corner cases ($l_{\text{max}}=8$, $n_{\text{samples}}=10,000$,
$n_{\text{threads}}=16$) the single-precision version can be _slower_ than
that using 64-bit floating-point values.

| 1 OpenMP thread | 16 OpenMP threads
---|---|---
|  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$
| Single Precision
$l_{\text{max}}=1$ |  1.08 |  2.00 |  0.572 |  0.680
$l_{\text{max}}=2$ |  3.74 |  12.4 |  0.564 |  1.10
$l_{\text{max}}=4$ |  9.47 |  35.8 |  0.922 |  2.58
$l_{\text{max}}=8$ |  56.4 |  169 |  17.5 |  54.2
$l_{\text{max}}=16$ |  240 |  1093 |  28.9 |  126
$l_{\text{max}}=32$ |  1128 |  3101 |  109 |  366
| Double Precision
$l_{\text{max}}=1$ |  1.06 |  6.47 |  0.451 |  0.573
$l_{\text{max}}=2$ |  4.02 |  15.9 |  0.578 |  1.28
$l_{\text{max}}=4$ |  11.6 |  47.3 |  1.00 |  3.73
$l_{\text{max}}=8$ |  57.0 |  280 |  5.15 |  21.7
$l_{\text{max}}=16$ |  252 |  1602 |  18.3 |  164
$l_{\text{max}}=32$ |  1385 |  4651 |  148 |  440
Table 2: CPU-parallel execution time (in ns/point) for computing Cartesian
spherical harmonics (or Cartesian spherical harmonics and their derivatives)
up to the indicated value of $l_{\text{max}}$. Timings refer to 1 and 16
OpenMP threads respectively on Intel Xeon Gold 6226R CPU, averaged over 1000
calls for 10 000 points.

###  III.2 CUDA implementation details

We adapt the implementation to a custom CUDA kernel, that allows efficient
execution on GPUs. Similarly to the CPU implementation, we use hard-coded
spherical harmonics and derivatives but up to a reduced hard-coded
$l_{\text{max}}=3$ to lower shared memory requirements, and we use the general
expression for the remaining terms. We parallelize the computation with a two-
dimensional thread-block of $16\times 8$ threads, using a grid dimension of
$(n_{\text{samples}}/16)$ blocks. The first dimension in the thread-block
parallelizes over samples, while the first thread in the second dimension is
responsible for performing the computation, and the remaining perform
coalesced writing of the temporary buffers to global memory. We store the
intermediary work buffers for the spherical harmonics and derivatives in fast
shared memory, which defaults to 48KB for most NVIDIA cards. Memory accesses
are designed such that each sample-dimension thread is mapped to shared memory
bank contiguously, eliminating possible bank conflicts. The CUDA wrapper
automatically re-allocates the necessary amount of shared memory beyond the
default 48KB on the first compute call, or attempts to reduce the number of
threads launched in the CUDA kernel if the required allocation is too large.
As a consequence, for accurate timings we call the forward step once to
perform the initialization and then measure timings thereafter. We note that
using the recursion (7) would require storing temporary values for all $l\leq
l_{\text{max}}$, so that the default 48KB shared memory allocation would be
filled for $l_{\text{max}}\approx 8$ when computing derivatives in 64-bit
floating-point format. Using the alternative recursion (8) allows us to
evaluate one $l$ channel at a time, so values up to $l_{\text{max}}\approx 30$
can be computed without adjusting the shared memory allocation or reducing the
number of sample-dimension threads. For GPUs which support more than 48KB
shared memory per block, for example the A100 (164KB), the recursion (8)
allows us to evaluate values up to $l_{\text{max}}\approx 82$.

Table 3 shows timings for GPU-accelerated computations on an A100 SXM4 (80GB)
card. In general, using single-precision floating-point arithmetic results in
half the computation time than double-precision arithmetic. We note that the
GPU per-sample timings reduce significantly upon increasing number of samples,
as the GPU is not fully saturated for $n_{\text{samples}}=10\,000$. For
example for $n_{\text{samples}}=100\,000$, the per-sample timings reduce by a
factor of 1.5 to 4. Although the parallel CPU implementation is faster than
the GPU for $l_{\text{max}}\leq 4$ when using $n_{\text{samples}}=10\,000$,
this trend reverses for higher values of $l_{\text{max}}$. For example, the
64-bit floating-point computations with $n_{\text{samples}}=10\,000$ and
$l_{\text{max}}=16$ show the GPU outperforming the CPU by a factor of 8.

We note in closing that, even though precise timings depend on the hardware
and the maximum angular momentum considered, sphericart is between 10 and 40
times faster than the spherical harmonics implementation in e3nn38, a widely
used library for equivariant neural networks. We provide wrappers that use the
same conventions as e3nn, to simplify integration of our implementation in
existing models and to accelerate computational frameworks that are limited by
the evaluation of $Y_{l}^{m}$.

| 10k points | 100k points
---|---|---
|  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m}$ |  $\tilde{Y}_{l}^{m},\gradient\tilde{Y}_{l}^{m}$
| Single Precision
$l_{\text{max}}=1$ |  2.1 |  2.4 |  0.3 |  0.4
$l_{\text{max}}=2$ |  2.2 |  2.4 |  0.3 |  0.4
$l_{\text{max}}=4$ |  2.4 |  3.1 |  0.5 |  0.8
$l_{\text{max}}=8$ |  3.1 |  4.8 |  1.0 |  3.3
$l_{\text{max}}=16$ |  5.6 |  15.1 |  2.8 |  10.9
$l_{\text{max}}=32$ |  15.5 |  43.6 |  9.6 |  33.5
| Double Precision
$l_{\text{max}}=1$ |  2.2 |  2.6 |  0.4 |  0.6
$l_{\text{max}}=2$ |  2.3 |  2.9 |  0.4 |  0.8
$l_{\text{max}}=4$ |  2.7 |  4.0 |  0.7 |  1.7
$l_{\text{max}}=8$ |  4.2 |  8.5 |  1.7 |  4.6
$l_{\text{max}}=16$ |  9.3 |  22.2 |  5.2 |  16.0
$l_{\text{max}}=32$ |  28.8 |  75.0 |  20.6 |  65.1
Table 3: GPU-parallel execution time (in ns/point) for computing Cartesian
spherical harmonics (or Cartesian spherical harmonics and their derivatives)
up to the indicated value of $l_{\text{max}}$. Timings refer to a single A100
SXM4/80GB GPU, and are averaged over 10 000 calls on 10 000 and 100 000 points
respectively.

##  IV Conclusions and perspectives

Spherical harmonics are ubiquitous in the computational sciences, and their
efficient calculation has become even more important in light of the
widespread adoption of equivariant machine-learning models in chemistry.
Reformulating the calculation of spherical harmonics in a scaled form that
corresponds to real-valued polynomials of the Cartesian coordinates provides
simple expressions for their recursive evaluation. Derivatives can be obtained
with little overhead, re-using the same factors that enter the definition of
the scaled $\tilde{Y}_{l}^{m}$, and without numerical instabilities. The
conventional form of the spherical harmonics can be recovered easily by
computing $\tilde{Y}_{l}^{m}$ at $(x/r,y/r,z/r)$ and applying the
corresponding correction to the derivatives. In many applications, this
normalization may be unnecessary, as the scaling can be incorporated
(explicitly or through regression) into additional radial terms. We provide an
efficient implementation as a C++ library, complete with a C API, Python and
PyTorch bindings, that tackles the most common use case in scientific
computing and geometric machine learning, i.e., the evaluation of real-valued
spherical harmonics for all angular momentum channels up to a given cutoff
$l_{\text{max}}$ and possibly for many 3D points at once. The function call is
parallelized over the sample direction, and it uses more efficient hard-coded
expressions for low-$l$ terms. Future efforts will focus on the extension of
the software library to different programming languages and frameworks, as
well as on further optimization on new hardware platforms and improved
parallelism.

## Data and software availability

The _sphericart_ library can be freely downloaded under the Apache License
version 2.0 from the public git repository https://github.com/lab-
cosmo/sphericart. A Python package that provides a convenient class to compute
spherical harmonics and derivatives from a _NumPy_39 or a _PyTorch_36 array of
3D coordinates is available on PyPI at https://pypi.org/project/sphericart/.
The timing data can be generated using the benchmarking code that is included
in the software distribution.

## AUTHOR CONTRIBUTIONS

Filippo Bigi and Michele Ceriotti conceived the project and performed
analytical derivations. All authors contributed to software development and to
the writing of the paper.

## Acknowledgements

The Authors would like to thank Kevin Kazuki Huguenin-Dumittan for useful
discussions, Philip Loche for help with the sphericart library, and Prashanth
Kanduri for help with the continuous integration framework. MC acknowledges
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No
101001890-FIAMMA). GF acknowledges support from the Platform for Advanced
Scientific Computing (PASC). MC and FB acknowledge support from the NCCR
MARVEL, funded by the Swiss National Science Foundation (grant number 182892).

## References

  * Müller (1966) C. Müller,  _Spherical harmonics_ (Springer, 1966).
  * Steiner (1963) E. Steiner, The Journal of Chemical Physics 39, 2365 (1963).
  * Rexer  _et al._ (2016) M. Rexer, C. Hirt, S. Claessens, and R. Tenzer, Surveys in Geophysics 37, 1035 (2016).
  * Knaack and Stenflo (2005) R. Knaack and J. O. Stenflo, Astronomy and Astrophysics 438, 349 (2005).
  * Morschhauser  _et al._ (2014) A. Morschhauser, V. Lesur, and M. Grott, Journal of Geophysical Research: Planets 119, 1162 (2014).
  * Evans (1998) K. F. Evans, Journal of the Atmospheric Sciences 55, 429 (1998).
  * poletti (2005) m. a. poletti, journal of the audio engineering society 53, 1004 (2005).
  * Max and Getzoff (1988) N. L. Max and E. D. Getzoff, IEEE Computer Graphics and Applications 8, 42 (1988).
  * Sloan (2008) P.-P. Sloan, in  _Game developers conference_ , Vol. 9 (2008) p. 42.
  * Schlegel and Frisch (1995) H. B. Schlegel and M. J. Frisch, International Journal of Quantum Chemistry 54, 83 (1995).
  * Varganov  _et al._ (2008) S. A. Varganov, A. T. Gilbert, E. Deplazes, and P. M. Gill, The Journal of chemical physics 128, 201104 (2008).
  * Gill and Gilbert (2009) P. M. Gill and A. T. Gilbert, Chemical Physics 356, 86 (2009).
  * Maintz  _et al._ (2016) S. Maintz, M. Esser, and R. Dronskowski, Acta Physica Polonica B 47 (2016).
  * Pérez-Jordá and Yang (1996) J. M. Pérez-Jordá and W. Yang, The Journal of Chemical Physics 104, 8003 (1996).
  * Choi  _et al._ (1999) C. H. Choi, J. Ivanic, M. S. Gordon, and K. Ruedenberg, The Journal of Chemical Physics 111, 8825 (1999).
  * Ding  _et al._ (2017) L. Ding, M. Levesque, D. Borgis, and L. Belloni, The Journal of Chemical Physics 147, 094107 (2017).
  * Zotkin  _et al._ (2009) D. N. Zotkin, R. Duraiswami, and N. A. Gumerov, in  _2009 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics_ (IEEE, 2009) pp. 257–260.
  * Li  _et al._ (2011) X. Li, S. Yan, X. Ma, and C. Hou, Applied Acoustics 72, 646 (2011).
  * Cohen and Welling (2016) T. Cohen and M. Welling, in  _Int. Conf. Mach. Learn._ (PMLR, 2016) pp. 2990–2999.
  * Bronstein  _et al._ (2021) M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković, arxiv:2104.13478 (2021).
  * Thomas  _et al._ (2018) N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, arxiv:1802.08219 (2018).
  * Anderson  _et al._ (2019) B. Anderson, T. S. Hy, and R. Kondor, in  _NeurIPS_ (2019) p. 10.
  * Klicpera  _et al._ (2021) J. Klicpera, F. Becker, and S. Günnemann, arxiv:2106.08903 (2021).
  * Geiger and Smidt (2022) M. Geiger and T. Smidt, arXiv preprint arXiv:2207.09453 (2022).
  * Bartók  _et al._ (2013) A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Thompson  _et al._ (2015) A. Thompson, L. Swiler, C. Trott, S. Foiles, and G. Tucker, Journal of Computational Physics 285, 316 (2015).
  * Willatt  _et al._ (2019) M. J. Willatt, F. Musil, and M. Ceriotti, J. Chem. Phys. 150, 154110 (2019).
  * Drautz (2019) R. Drautz, Phys. Rev. B 99, 014104 (2019).
  * Musil  _et al._ (2021) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Chem. Rev. 121, 9759 (2021).
  * Christensen and Von Lilienfeld (2020) A. S. Christensen and O. A. Von Lilienfeld, Machine Learning: Science and Technology 1, 045018 (2020).
  * Schütt  _et al._ (2021) K. Schütt, O. Unke, and M. Gastegger, in  _Int. Conf. Mach. Learn._ (PMLR, 2021) pp. 9377–9388.
  * Press (2007) W. H. Press,  _Numerical Recipes: The Art of Scientific Computing_ (Cambridge University Press, 2007).
  * Sloan (2013) P.-P. Sloan, J. Comput. Graph. Tech. JCGT 2, 84 (2013).
  * Drautz (2020) R. Drautz, Phys. Rev. B 102, 024104 (2020).
  * Alken (2022) P. Alken, (2022).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Dagum and Menon (1998) L. Dagum and R. Menon, IEEE computational science and engineering 5, 46 (1998).
  * Geiger  _et al._ (2022) M. Geiger, T. Smidt, A. M., B. K. Miller, W. Boomsma, B. Dice, K. Lapchevskyi, M. Weiler, M. Tyszkiewicz, S. Batzner, D. Madisetti, M. Uhrin, J. Frellsen, N. Jung, S. Sanborn, M. Wen, J. Rackers, M. Rød, and M. Bailey, “Euclidean neural networks: e3nn,”  (2022).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, _et al._ , Nature 585, 357 (2020).
  * Bosch (2000) W. Bosch, Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy 25, 655 (2000).

##  Appendix A Recurrence relations for the Cartesian spherical harmonics

In this Appendix, we derive the Cartesian form of the recurrence relations
presented in the main text. Starting from Eq. 2, it is sufficient to multiply
both sides by $r^{l}$ and multiply and divide the right-hand side by
$r_{xy}^{|m|}$ to obtain

$r^{l}Y_{l}^{m}=F_{l}^{|m|}*r^{l}r_{xy}^{-|m|}P_{l}^{|m|}*\begin{cases}r_{xy}^{|m|}\sin{(|m|\phi)}&\text{if
$m<0$}\\\ r_{xy}^{0}/\sqrt{2}&\text{if $m=0$}\\\
r_{xy}^{m}\cos{(m\phi)}&\text{if $m>0$}\end{cases}$ (12)

Then, using the definitions in Eq. 6, Eq. 5 follows.

Let us now derive Eq. 7. These relationships can be obtained as a combination
of Eqs. 4 and the definiton of $Q_{l}^{m}$ in Eq. 8, with the additional
observation that $t=\cos{\theta}=z/r$:

$Q_{0}^{0}=r^{0}r_{xy}^{0}P_{0}^{0}=P_{0}^{0}=1$ (13)
$Q_{m}^{m}=r^{m}r_{xy}^{-m}P_{m}^{m}=-\,(2m-1)\,\sqrt{1-t^{2}}\,r^{m}r_{xy}^{-m}P_{m-1}^{m-1}=\\\
-\,(2m-1)\,\sqrt{1-(z/r)^{2}}\,(r/r_{xy})\,Q_{m-1}^{m-1}=-\,(2m-1)\,Q_{m-1}^{m-1}$
(14) $Q_{m}^{m-1}=r^{m}r_{xy}^{-(m-1)}P_{m}^{m-1}=\\\
(2m-1)\,rt\,r^{-(m-1)}r_{xy}^{-(m-1)}P_{m-1}^{m-1}=(2m-1)\,z\,Q_{m-1}^{m-1}$
(15) $Q_{l}^{m}=r^{l}r_{xy}^{-m}P_{l}^{m}=\\\
r^{l}r_{xy}^{-m}((2l-1)\,t\,P_{l-1}^{m}+(l+m-1)\,P_{l-2}^{m})/(l-m)=\\\
((2l-1)\,rt\,r^{l-1}r_{xy}^{-m}P_{l-1}^{m}+(l+m-1)\,r^{2}r^{l-2}r_{xy}^{-m}P_{l-2}^{m})/(l-m)=\\\
((2l-1)\,z\,Q_{l-1}^{m}+(l+m-1)\,r^{2}Q_{l-2}^{m})/(l-m)$ (16)

And, finally, we can derive the recurrence relations for $s_{m}$ and $c_{m}$
in the following way:

$s_{m}=r_{xy}^{m}\sin{m\phi}=r_{xy}^{m}\sin{((m-1)\phi+\phi)}=\\\
r_{xy}^{m}(\sin{((m-1)\phi)}\cos{\phi}+\cos{((m-1)\phi)}\sin{\phi})=\\\
s_{m-1}r_{xy}\cos{\phi}+c_{m-1}r_{xy}\sin{\phi}=xs_{m-1}+yc_{m-1}$ (17)
$c_{m}=r_{xy}^{m}\cos{m\phi}=r_{xy}^{m}\cos{((m-1)\phi+\phi)}=\\\
r_{xy}^{m}(\cos{((m-1)\phi)}\cos{\phi}-\sin{((m-1)\phi)}\sin{\phi})=\\\
c_{m-1}r_{xy}\cos{\phi}-s_{m-1}r_{xy}\sin{\phi}=xc_{m-1}-ys_{m-1}$ (18)

where we have used some well-known trigonometric relations and the fact that
$\sin{\phi}=y/r_{xy}$, $\cos{\phi}=x/r_{xy}$.

##  Appendix B Cartesian derivatives of the scaled spherical harmonics

In this Appendix, we prove the derivative formulas presented in the main text,
i.e., Eq. 10. Let us start from the derivatives of $s_{m}$ and $c_{m}$. To
this end, we can rewrite Eq. 9 in matrix form:

$\begin{pmatrix}s_{m}\\\ c_{m}\end{pmatrix}=\begin{pmatrix}x&y\\\
-y&x\end{pmatrix}\begin{pmatrix}s_{m-1}\\\ c_{m-1}\end{pmatrix},$ (19)

from which it is easy to see that

$\begin{pmatrix}s_{m}\\\ c_{m}\end{pmatrix}=\begin{pmatrix}x&y\\\
-y&x\end{pmatrix}^{m}\begin{pmatrix}s_{0}\\\ c_{0}\end{pmatrix}.$ (20)

Differentiation with respect to $x$ yields

$\frac{\partial}{\partial x}\begin{pmatrix}s_{m}\\\
c_{m}\end{pmatrix}=m\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\begin{pmatrix}x&y\\\
-y&x\end{pmatrix}^{m-1}\begin{pmatrix}s_{0}\\\ c_{0}\end{pmatrix}=\\\
m\begin{pmatrix}s_{m-1}\\\ c_{m-1}\end{pmatrix},$ (21)

which proves the $x$ derivatives of $s_{m}$ and $c_{m}$ in Eq. 10, while
differentiation with respect to $y$ gives

$\frac{\partial}{\partial y}\begin{pmatrix}s_{m}\\\
c_{m}\end{pmatrix}=m\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\begin{pmatrix}x&y\\\
-y&x\end{pmatrix}^{m-1}\begin{pmatrix}s_{0}\\\ c_{0}\end{pmatrix}=\\\
m\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}\begin{pmatrix}s_{m-1}\\\
c_{m-1}\end{pmatrix}=m\begin{pmatrix}c_{m-1}\\\ -s_{m-1}\end{pmatrix},$ (22)

which results in the $y$ derivatives of $s_{m}$ and $c_{m}$ in Eq. 10.

We can now turn to the $\partial Q_{l}^{m}/\partial z$ derivative in Eq. 10,
which we prove by induction over the $l$ variable. The base case ($l=m+1$) can
easily be checked for any $m$ by considering the third equality in Eq. 7,
which implies

$Q_{m+1}^{m}=(2m+1)\,z\,Q_{m}^{m}.$ (23)

Since $Q_{m}^{m}$ does not depend on $z$ (see Eq. 10), differentiating with
respect to $z$ gives

$\frac{\partial Q_{m+1}^{m}}{\partial z}=(2m+1)\,Q_{m}^{m}.$ (24)

This is indeed exactly the $\partial Q_{l}^{m}/\partial z$ derivative in Eq.
10 for $l=m+1$. In the induction step, we prove the $l$ case from the $l-1$
and $l-2$ cases, so it is also necessary to prove the $l=m+2$ case in advance.
From the last line in Eq. 7, we can write

$Q_{m+2}^{m}=\frac{1}{2}((2m+3)\,z\,Q_{m+1}^{m}-(2m+1)\,r^{2}\,Q_{m}^{m}).$
(25)

Differentiation leads to

$\frac{\partial{Q_{m+2}^{m}}}{\partial
z}=\frac{1}{2}((2m+3)\,(Q_{m+1}^{m}+z\,\frac{\partial{Q_{m+1}^{m}}}{\partial
z})\\\ -(2m+1)(2z\,Q_{m}^{m}+r^{2}\,\frac{\partial{Q_{m}^{m}}}{\partial z})).$
(26)

Using $\partial{Q_{m}^{m}}/\partial z=0$ (see Eq. 7) and
$\partial{Q_{m+1}^{m}}/\partial z=(2m+1)\,Q_{m}^{m}$ (which we have just
proved), we obtain

$\frac{\partial{Q_{m+2}^{m}}}{\partial
z}=\frac{1}{2}((2m+3)\,(Q_{m+1}^{m}+z\,(2m+1)\,Q_{m}^{m})-\\\
(2m+1)\,2z\,Q_{m}^{m})=(2m+2)\,Q_{m+1}^{m},$ (27)

where we have used $(2m+1)\,z\,Q_{m}^{m}=Q_{m+1}^{m}$ from the third equality
of Eq. 7 and hidden some elementary algebra. This corresponds to the $\partial
Q_{l}^{m}/\partial z$ derivative in Eq. 10 for $l=m+2$. For the induction
step, we can now assume $\partial Q_{l-1}^{m}/\partial z=(l+m-1)\,Q_{l-2}^{m}$
and $\partial Q_{l-2}^{m}/\partial z=(l+m-2)\,Q_{l-3}^{m}$. Differentiating
the last equality of Eq. 7 with respect to $z$ on both sides results in

$\frac{\partial Q_{l}^{m}}{\partial
z}=((2l-1)\,Q_{l-1}^{m}+(2l-1)\,z\,\frac{\partial Q_{l-1}^{m}}{\partial z}-\\\
(l+m-1)\,2z\,Q_{l-2}^{m}-(l+m-1)\,r^{2}\,\frac{\partial Q_{l-2}^{m}}{\partial
z})/(l-m).$ (28)

We can now insert the assumptions of the induction step into the right-hand-
side expression to obtain

$\frac{\partial Q_{l}^{m}}{\partial
z}=((2l-1)\,Q_{l-1}^{m}+(2l-1)\,z\,(l+m-1)\,Q_{l-2}^{m}-\\\
(l+m-1)\,2z\,Q_{l-2}^{m}-(l+m-1)\,r^{2}\,(l+m-2)\,Q_{l-3}^{m})/(l-m).$ (29)

An elementary manipulation of the $z\,Q_{l-2}^{m}$ terms affords

$\frac{\partial Q_{l}^{m}}{\partial
z}=((2l-1)\,Q_{l-1}^{m}+(l+m-1)(2l-3)\,z\,Q_{l-2}^{m}-\\\
(l+m-1)\,(l+m-2)\,r^{2}\,Q_{l-3}^{m})/(l-m).$ (30)

Finally, application of the last equality of Eq. 7 absorbs the $Q_{l-2}^{m}$
and $Q_{l-3}^{m}$ terms into

$\frac{\partial Q_{l}^{m}}{\partial
z}=((2l-1)\,Q_{l-1}^{m}+(l+m-1)^{2}\,Q_{l-1}^{m})/(l-m),$ (31)

which simplifies into

$\frac{\partial Q_{l}^{m}}{\partial z}=(l+m)\,Q_{l-1}^{m},$ (32)

as required.

To prove the expressions for $\partial Q_{l}^{m}/\partial x$ and $\partial
Q_{l}^{m}/\partial y$ in Eq. 10, we first note that, since $Q_{l}^{m}$ is
formally only a function of $z$ and $r$, any $x$\- or $y$-dependence in
$Q_{l}^{m}$ must come from the involvement of the $r$ variable, so that

$\frac{\partial Q_{l}^{m}}{\partial x}=\frac{\partial Q_{l}^{m}}{\partial
r}\frac{\partial r}{\partial x}=\frac{\partial Q_{l}^{m}}{\partial
r}\frac{x}{r}$ (33)

and

$\frac{\partial Q_{l}^{m}}{\partial y}=\frac{\partial Q_{l}^{m}}{\partial
r}\frac{\partial r}{\partial y}=\frac{\partial Q_{l}^{m}}{\partial
r}\frac{y}{r}.$ (34)

Hence, in order to justify the $x$ and $y$ derivatives of $Q_{l}^{m}$ in Eq.
10, we simply need to prove that

$\frac{\partial Q_{l}^{m}}{\partial r}=r\,Q_{l-1}^{m+1}.$ (35)

Before doing so, we borrow Eq. 1 from Ref. 40. It is important to note how
Ref. 40 follows a different convention, and it does not include a $(-1)^{m}$
factor in the definition of the associated Legendre polynomials. Hence,
compared to Ref. 40, we change the sign of the $P_{l-1}^{m-1}$ term to obtain

$P_{l}^{m}=-\,(2l-1)\,\sin{\theta}P_{l-1}^{m-1}+P_{l-2}^{m},$ (36)

which is now consistent with our conventions. Noting that
$\sin{\theta}=r_{xy}/r$ and multiplying by $r^{l}r_{xy}^{-m}$ on both sides
gives

$r^{l}r_{xy}^{-m}P_{l}^{m}=\\\
-\,(2l-1)(r_{xy}/r)\,r^{l}r_{xy}^{-m}P_{l-1}^{m-1}+r^{l}r_{xy}^{-m}P_{l-2}^{m},$
(37)

which, thanks to the definition of $Q_{l}^{m}$ in Eq. 6, translates to

$Q_{l}^{m}=-\,(2l-1)\,Q_{l-1}^{m-1}+r^{2}Q_{l-2}^{m}.$ (38)

Eq. 38 and the last line of Eq. 7 provide two expressions for $Q_{l}^{m}$.
Imposing their equality results in

$((2l-1)\,z\,Q_{l-1}^{m}+(l+m-1)\,r^{2}\,Q_{l-2}^{m})/(l-m)=\\\
-\,(2l-1)\,Q_{l-1}^{m-1}+r^{2}Q_{l-2}^{m}.$ (39)

This can be rearranged into

$z\,Q_{l-1}^{m}=-\,(l-m)\,Q_{l-1}^{m-1}+r^{2}Q_{l-2}^{m}.$ (40)

From this equation, $Q_{l-1}^{m-1}$ can be extracted as

$Q_{l-1}^{m-1}=(-\,z\,Q_{l-1}^{m}+r^{2}Q_{l-2}^{m})/(l-m).$ (41)

Let us now go back to Eq. 35. We prove it by induction over the $l$ variable,
similar to what we did for the $z$-derivative. Given the $m\leq l$ constraint,
the base case reads $\partial Q_{m+1}^{m-1}/\partial r=r\,Q_{m}^{m}.$ From the
third equality in Eq. 7, we have $Q_{m}^{m-1}=(2m-1)\,z\,Q_{m-1}^{m-1}$. Now,
using the last line of Eq. 7 with $l=m+1$ and (abusing notation) $m=m-1$, we
can express $Q_{m+1}^{m-1}$ as

$Q_{m+1}^{m-1}=((2m+1)\,z\,Q_{m}^{m-1}-(2m-1)\,r^{2}\,Q_{m-1}^{m-1})/2=\\\
((2m+1)(2m-1)\,z^{2}\,Q_{m-1}^{m-1}-(2m-1)\,r^{2}\,Q_{m-1}^{m-1})/2=\\\
(2m-1)\,Q_{m-1}^{m-1}\,((2m+1)\,z^{2}-r^{2})/2=\\\
-\,Q_{m}^{m}\,((2m+1)\,z^{2}-r^{2})/2,$ (42)

where we have used the second equality of Eq. 7 in the last line. Since
$Q_{m}^{m}$ does not depend on $r$ (see Eq. 7), differentiating both sides of
Eq. 42 with respect to $r$ yields

$\frac{\partial Q_{m+1}^{m-1}}{\partial r}=r\,Q_{m}^{m},$ (43)

as required. As in the case of the $z$ derivative, the base case also includes
a further equality: $\partial Q_{m+2}^{m-1}/\partial r=r\,Q_{m}^{m+1}$. Using
the last line of Eq. 7, $Q_{m+2}^{m-1}$ can be written as

$Q_{m+2}^{m-1}=\frac{1}{3}((2m+3)\,z\,Q_{m+1}^{m-1}-2m\,r^{2}\,Q_{m}^{m-1}).$
(44)

partial differentiation of the above with respect to $r$ yields

$\frac{\partial Q_{m+2}^{m-1}}{\partial
r}=\frac{1}{3}((2m+3)\,z\,\frac{\partial Q_{m+1}^{m-1}}{\partial
r}-2m\,(2r\,Q_{m}^{m-1}+r^{2}\,\frac{\partial Q_{m}^{m-1}}{\partial r})).$
(45)

However, we have proved $\partial Q_{m+1}^{m-1}/\partial r=r\,Q_{m}^{m}$, and
we have $\partial Q_{m}^{m-1}/\partial r=0$ from the first three equalities of
Eq. 7. Hence,

$\frac{\partial Q_{m+2}^{m-1}}{\partial
r}=\frac{1}{3}((2m+3)\,z\,r\,Q_{m}^{m}-4mr\,Q_{m}^{m-1})=\\\
\frac{1}{3}((2m+3)\,z\,r\,Q_{m}^{m}+4mrz\,Q_{m}^{m})=(2m+1)\,zr\,Q_{m}^{m}=r\,Q_{m}^{m+1},$
(46)

where we have used the third equality of 7 twice. This concludes the proof of
the base cases.

To prove the induction step, we start by differentiating both sides of the
last equality of Eq. 7 with respect to $r$. We obtain

$\frac{\partial Q_{l}^{m}}{\partial r}=((2l-1)\,z\frac{\partial
Q_{l-1}^{m}}{\partial r}-2(l+m-1)\,r\,Q_{l-2}^{m}-\\\
(l+m-1)\,r^{2}\,\frac{\partial Q_{l-2}^{m}}{\partial r})/(l-m).$ (47)

We can now assume that $\partial Q_{l-1}^{m}/\partial r=r\,Q_{l-2}^{m+1}$ and
$\partial Q_{l-2}^{m}/\partial r=r\,Q_{l-3}^{m+1}$ to get

$\frac{\partial Q_{l}^{m}}{\partial
r}=((2l-1)\,zr\,Q_{l-2}^{m+1}-2(l+m-1)\,r\,Q_{l-2}^{m}-\\\
(l+m-1)\,r^{3}\,Q_{l-3}^{m+1})/(l-m).$ (48)

Now, Eq. 41 implies
$Q_{l-2}^{m}=(-\,z\,Q_{l-2}^{m+1}+r^{2}Q_{l-3}^{m+1})/(l-m-2)$. Substituting
this into Eq. 48 leads to

$\frac{\partial Q_{l}^{m}}{\partial r}=((2l-1)\,zr\,Q_{l-2}^{m+1}-\\\
\frac{2(l+m-1)}{l-m-2}\,(-\,zr\,Q_{l-2}^{m+1}+r^{3}Q_{l-3}^{m+1})-\\\
(l+m-1)\,r^{3}\,Q_{l-3}^{m+1})/(l-m).$ (49)

Adding like terms (i.e., those in $zr\,Q_{l-2}^{m+1}$ and those in
$r^{3}\,Q_{l-3}^{m+1}$) results in

$\frac{\partial Q_{l}^{m}}{\partial
r}=\frac{1}{(l-m)(l-m-2)}((2l-3)(l-m)\,zr\,Q_{l-2}^{m+1}-\\\
(l+m-1)(l-m)\,r^{3}\,Q_{l-3}^{m+1})=\\\
\frac{1}{l-m-2}((2l-3)\,zr\,Q_{l-2}^{m+1}-(l+m-1)\,r^{3}\,Q_{l-3}^{m+1})=\\\
r\frac{1}{l-m-2}((2l-3)\,z\,Q_{l-2}^{m+1}-(l+m-1)\,r^{2}\,Q_{l-3}^{m+1}).$
(50)

However, by considering the last line of Eq. 7 with $l$ decreased by one and m
increased by one, the last expression can be simplified into

$\frac{\partial Q_{l}^{m}}{\partial r}=r\,Q_{l-1}^{m+1},$ (51)

which concludes the proof.

##  Appendix C Considerations on the prefactor

The prefactors $F_{l}^{m}$ contain a ratio of factorials that can lead to
numerical instabilities in a naïve implementation. It is however easy to see
that one can compute them iteratively as

$F_{l}^{0}=\sqrt{\frac{2l+1}{2\pi}},\quad
F_{l}^{m}=-\frac{F_{l}^{m-1}}{\sqrt{(l+m)(l+1-m)}}.$ (52)

It is also possible to incorporate the prefactors in the definition of the
modified associated Legendre polynomials, defining
$\tilde{Q}_{l}^{m}=F_{l}^{m}{Q}_{l}^{m}$. This simplifies somehow the
construction of $\tilde{Y}_{l}^{m}$ and avoids possible instabilities
connected with the fact that $F_{l}^{m}$ become very small for large $m\approx
l$, at the price of complicating slightly the expressions for the recursion
and derivatives of $Q_{l}^{m}$, e.g.

$\frac{\partial\tilde{Q}_{l}^{m}}{\partial
z}=\sqrt{l^{2}-m^{2}}\,\tilde{Q}_{l-1}^{m}.$ (53)