# Probabilistic thermal stability prediction through sparsity promoting
transformer representation

Yevgen Zainchkovskyy
DTU Compute & Novo Nordisk A/S
<EMAIL_ADDRESS>Jesper Ferkinghoff-Borg
Novo Nordisk A/S
<EMAIL_ADDRESS>Anja Bennett
Novo Nordisk A/S
<EMAIL_ADDRESS>Thomas Egebjerg
Novo Nordisk A/S
<EMAIL_ADDRESS>Nikolai Lorenzen
Novo Nordisk A/S
<EMAIL_ADDRESS>Per Jr. Greisen
Novo Nordisk A/S
<EMAIL_ADDRESS>Søren Hauberg
DTU Compute
<EMAIL_ADDRESS>Carsten Stahlhut
Novo Nordisk A/S
<EMAIL_ADDRESS>Corresponding author

###### Abstract

Pre-trained protein language models have demonstrated significant
applicability in different protein engineering task [1, 2]. A general usage of
these pre-trained transformer models producing latent representation is to use
a mean pool across residue positions to reduce the feature dimensions to
further downstream tasks such as predicting bio-physical properties or other
functional behaviours. In this paper we provide a two-fold contribution to
machine learning (ML) driven drug design. Firstly, we demonstrate the power of
sparsity by promoting penalization of pre-trained transformer models to secure
more robust and accurate melting temperature (Tm) prediction of single-chain
variable fragments with a mean absolute error of $0.23^{\circ}$C. Secondly, we
demonstrate the power of framing our prediction problem in a probabilistic
framework. Specifically, we advocate for the need of adopting probabilistic
frameworks especially in the context of ML driven drug design.

##  1 Introduction

Peptide and protein engineering is the process of optimizing peptide and
proteins towards desired and valuable features for technological or medical
applications [3]. In protein engineering, we seek to optimize the function of
a protein with respect to e.g. its expression level, solubility, or thermal
stability. Their functional behavior is directly determined by their amino
acid sequence. Thus, to develop new or optimize desired properties for e.g.,
biomedical applications require to invert the relationship of the function
given the sequence [4], also generally known in statistics and machine
learning as the inverse problem. However, existing design methods have serious
problems in distinguishing the functional levels of closely related proteins
[5, 6]. While both protein engineering and design is a NP-hard problem [7], a
direct search in the protein space simply becomes an overwhelming and
intractable approach in linear time. Directed evolution has successfully
demonstrated its applicability of mapping peptide and protein sequencing to
functional behavior. However, it is highly limited by the fact that even high-
throughput techniques only can sample a minor fraction of sequences
constructed from diversification methods [8]. Among others, Bedbrook and co-
workers have demonstrated the direct applicability of utilizing machine
learning for optimizing a property that would not have been possible to
engineer through directed evolution alone, [6, 9, 10]. On the other hand,
machine learning models are heavily dependent on learning from data - a
crucial part in designing a machine learning driven drug design pipeline is
therefore the accessibility of relevant functional data for the task at hand.

In this contribution, we demonstrate and discuss classical challenges in both
designing compounds with dedicated properties from a minimal set of
observations and data sets with quite scarce diversity. While we at one hand
wish to provide as diverse molecules to maximize the coverage of the chemical
search space and on the other hand seek to ensure optimized properties within
minimum number of design rounds (experiments) - we are facing a typical active
learning problem balancing explore vs exploit steps through the usage of model
uncertainty estimates. In section 2, we provide a unified probabilistic
framework for integrating compact latent pre-trained transformer features with
Gaussian Process (GP) regression models, [11]. Through careful variant design
train, development (dev), and test splits we demonstrate the applicability of
uncertainty estimates to assess the models own notion of what it does not
know. We examine the effect of training data with 1-5 mutations away from a
wild-type sequence and the models ability to reason of its own predictive
power to generalize to multiple mutations. While we in this paper limit
ourselves to the quantification of predictive performance of the models, the
GPs can be utilized as the surrogate model in a Bayesian Optimization
framework for optimizing and searching the sequence space.

##  2 Background

Here, we motivate our problem and provide a brief overview of the core
architecture of our probabilistic models utilizing a transformer architecture
as input to our downstream regression models. Significant improvement and
applicability of protein language models have been demonstrated over the last
years, where among others the UniRep [12], Evolutionary Scale Modeling (ESM)
[2], ProtBert [13] models can be mentioned. In [12], they utilize pretrained
language model representation UniRep to generalize the representation to
unseen regions of sequence space. Furthermore, Vig & Rao argues for attention
in the transformer models corresponds to known biological properties like
structure and binding sites that can enable contact prediction [14, 15]. In a
drug design setting, we are especially interested in enabling pretrained
representations in our protein engineering tasks for designing improved drugs
as we are highly limited by the number of experiments we can conduct relative
to the enormous sequence space at hand, e.g. $10^{130}$ for proteins of 100
amino-acids length. Thus, searching the space intelligently is needed even
when utilizing high-throughput experimental setups.

While designing or engineering proteins, we are faced with the problem of
optimizing towards specific functions of the molecules and effectively only
interested in a tiny subspace of sequence space. The main challenge is
naturally how can we utilize the general protein representation for a direct
fine tuning to the downstream functional optimization task at hand. In this
contribution, we seek to build a model for predicting the thermal stability of
the antibody format single chain variable fragment (scFv). Due to the small
sizes and the stranded nature of scFvs, these are commonly used as building
blocks to construct recombinant multi-specific antibody formats, [16].
Unfortunately, the scFvs has been reported to be less termostable than larger
antibody formats and thus more likely to lead to undesired aggregation and low
Tms when utilized in a multi-specific format, [17]. To improve biophysical
behavior of the scFV, our goal is to build a predictive model of the
experimental measured Tm values determined by nano differential scanning
fluorimetry (nanoDSF) [18]. Having an accurate model for prediction the
melting temperature is needed to assess which variants to test experimentally
for increased thermal stability. To quantify the accurate of our predictions,
we follow a probabilistic approach where we not only obtain our mean
predictions but just as importantly can provide uncertainty estimates on the
predictions. Uncertainty estimates is critically needed for providing
quantitative and directed search strategies balancing both exploitation and
exploration.

###  2.1 Transformer-based models

Transformers are revolutionising NLP, have recently been repurposed to model
biological sequences.In the core of a transformer architecture is the
attention mechanism allowing to capture long-range dependencies between
positions in a sequence. Originating as a solution to classic sequence-to-
sequence (seq2seq) models, attention mechanism shows better performance and
scaling characteristics than traditional RNNs or LSTMs. Common to those
architectures is a context vector comprising of a hidden state of the network
being carried through subsequent propagation, resulting in degraded
performance with increased sequence lengths. On the other hand, transformers
utilize self-attention, which allows processing of the whole sequence while
still focusing on specific parts of it.

In the context of biological sequence modelling, the hidden state of a
Transformer model corresponds to individual amino acid residues and represents
the given amino acid in its context as a point in a high dimensional space
(embedding). Thus, similar sequences are assigned similar representations by
the network and are mapped to nearby points in space.

In this contribution, we use the ESM1-b variant of a Transformer protein
language model from Facebook AI Research [2] encoding each of our sequences to
an embedding $\mathbf{x}\in\mathbb{R}^{250\times 1280}$.

###  2.2 Gaussian process regression

A Gaussian Processes (GP) is a powerful probabilistic framework enabling
nonparametric, nonlinear Bayesian models [11]. A GP defines a prior
distribution over the set of function $f\left(\mathbf{x}\right)$ mapping the
relation between our $M$-dimensional feature representation of our protein
sequences to the target property $y=f\left(\mathbf{x}\right)+\epsilon$. Here
$\epsilon$ represents additive observation noise. Using a standard zero-mean
GP prior we obtain

$p\left(f\left(\mathbf{X}\right)\right)=p\left(\mathbf{f}_{X}\right)={{\mathcal{N}}\left(0,\mathbf{K}\right)},$
(1)

where $\mathbf{K}$ is the covariance matrix between our training input
features $\mathbf{X}=[\mathbf{x}_{1},\cdots,\mathbf{x}_{N}]$ such that
$K_{ij}=k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ defines the covariance
function between input $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. We utilize one
of the typically applied kernels for GP regression, Matern $\frac{5}{2}$
covariance function, with shared length-scale parameters for each input
dimension $\sigma_{l}$, yielding

$k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)=\sigma_{f}^{2}\left(1+\frac{\sqrt{3}r}{\sigma_{l}}\right)\exp\left(-\frac{\sqrt{3}r}{\sigma_{l}}\right)\quad\mathrm{where}\quad
r=\sqrt{\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)^{\top}\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)}$
(2)

Assuming additive independent identically distributed Gaussian noise with
variance $\sigma^{2}_{\epsilon}$ our predictive distribution for our new test
proteins $\mathbf{Z}$ reads, [11],

$\displaystyle p\left(f_{Z}|\mathbf{X},\mathbf{y},\mathbf{Z}\right)$
$\displaystyle=$
$\displaystyle\mathcal{N}\left(\mathbf{\mu}_{z},\mathbf{\Sigma}_{z}\right),\quad\mathrm{where}\quad\mathbf{\mu}_{z}=k\left(\mathbf{Z},\mathbf{X}\right)\left(\mathbf{K}+\sigma^{2}_{\epsilon}\mathbf{I}\right)^{-1}\mathbf{y}$
$\displaystyle\mathbf{\Sigma}_{z}$ $\displaystyle=$ $\displaystyle
k\left(\mathbf{Z},\mathbf{Z}\right)-k\left(\mathbf{Z},\mathbf{X}\right)\left(\mathbf{K}+\sigma^{2}_{\epsilon}\mathbf{I}\right)^{-1}k\left(\mathbf{X},\mathbf{Z}\right).$
(3)

##  3 Driving sparsity through learned masks

Having extracted an embedding for each residue in a sequence of length $P$, a
typical approach used to represent a complete protein as a single vector
$\hat{\mathbf{x}}$ is by averaging across the transformer’s hidden
representation $\mathbf{x}$ at each sequence position $p$ (mean pooling):

$\hat{\mathbf{x}}=\frac{1}{P}\sum\nolimits_{p=1}^{P}{\mathbf{x}_{p}}$ (4)

While significantly reducing the overall dimensionality of the embedding and
allowing to represent proteins of different lengths, this approach inevitably
results in loss of information. Intuitively, averaging assigns equal weight to
all residues in the sequence, while in reality, only a handful of positions
might influence the target of interest.

In this work, we investigate 3 different positional-weighted approaches to the
typical averaging: a positively constrained Learned mask $\mathbf{w_{l}}$,
sparsity promoting Sigmoid-transformed mask $\mathbf{w_{s}}$ and a Half-Cauchy
Prior based mask $\mathbf{w_{p}}$:

$\mathbf{\hat{x}_{l}}=\frac{\sum_{p=1}^{P}{\exp({\mathbf{w_{l}}_{p}})\mathbf{x}_{p}}}{\sum{\exp({\mathbf{w_{l}})}}}$
(5)

$\mathbf{\hat{x}_{s}}=\frac{\sum_{p=1}^{P}{\mathrm{S}({\mathbf{w_{s}}_{p}})\mathbf{x}_{p}}}{\sum{\mathrm{S}({\mathbf{w_{s}})}}}\quad\textrm{where}\quad\mathrm{S}(x)=\frac{1}{1+\exp(-x)}$
(6)

$\mathbf{\hat{x}_{p}}=\frac{\sum_{p=1}^{P}{{\mathbf{w_{p}}_{p}}\mathbf{x}_{p}}}{\sum{{\mathbf{w_{p}}}}}\quad\textrm{where}\quad\mathbf{w_{p}}\sim\textrm{Half-
Cauchy}(0,\sigma)$ (7)

From a practical perspective, for all three approaches, we learn the masks as
a part of the standard gradient based GP Maximum Log-Likelihood maximization
procedure.

##  4 Results

Using Mean Absolute Error (MAE) as our metric, and partitioning the data set,
we ran 3 sets of experiments corresponding to 3 different splits. First, we
tested on the subset of the training set where single-site mutations were used
as the validation (1MUT), next we used a random sample of the training set as
validation (Uniform Shuffle) and finally evaluated the proposed methods on the
hold-out test set itself. The reason for this partitioning is the spread of
the positions of mutations in the wild-type sequence and relative sizes of
training and validation set. Naturally, having learned a specific mask on the
smaller training set, performance will degrade if mask does not reflect the
mutated positions in the bigger test-set. This effect is simulated for the
"Uniform Shuffle" split where the size of validation set was 24 samples (vs.
10 samples for the 1-MUT split).

Table 1: Evaluation of the Baseline and the proposed methods. For each split and method we report mean $\pm$ std over 64 trials. Bold values denote statistical significance against Baseline ($p<0.05).$ | 1-MUT Shuffle | Uniform Shuffle | Test Set (Hold out)
---|---|---|---
Baseline | 1.216 $\pm$ 0.306 | 0.819 $\pm$ 0.157 | 0.273 $\pm$ 0.006
Learned mask | $\mathbf{1.110\pm 0.338}$ | 0.836 $\pm$ 0.231 | $\mathbf{0.227\pm 0.010}$
Learned mask (Sigmoid) | $\mathbf{1.150\pm 0.347}$ | 0.813 $\pm$ 0.161 | $\mathbf{0.227\pm 0.008}$
Learned mask (Prior) | $\mathbf{1.090\pm 0.314}$ | 0.841 $\pm$ 0.158 | $\mathbf{0.226\pm 0.010}$

Gross metrics are reported in Table 1. Here, we see that the Learned Mask
(prior) proposed method outperforms the standard averaging approach referenced
as Baseline in the cases of balanced train/validation splits. Naturally, this
method also depends on the actual value of the prior: $\sigma$ \- the second
moment of Half-Cauchy distribution. In our case, a prior of $\sigma=0.15$ was
chosen as a result of a parameter sweep on the 1-MUT Shuffle resulting in
lowest MAE. That prior is then reused for all Learned Mask (prior) runs.

Figure 1: Masks learned on the hold out set. The number of zero entries
(values below $1e^{-5}$ threshold) are: 208 for Learned mask, 203 for Learned
mask (Sigmoid) and 9 for Learned mask (Prior). Unused entries are marked with
a cross $\times$.

Learned masks are shown in Figure 1. In terms of the importance on melting
temperature they emphasize roughly same positions in the sequence.
Interestingly, the addition of Half-Cauchy prior, results in a more _dense_
mask, as seen on the number of zero-entries. This appear to help
generalization.

We summarize the individual test-set predictions and their corresponding
uncertainties in Figure 2 (full version in Appendix, Figure 6).

Figure 2: A subset of melting temperature prediction and corresponding
uncertainties on the test-set.

Additionally, we train a model on a subset of the training-data comprising of
only 80 samples all being single site mutations. While this restricted model
does not perform particularly well in terms of our metric (resulting
$\textrm{MAE}=1.33$), it does great job in terms of uncertainties as shown in
Figure 3 (full version in Appendix, Figure 7). Note that this limited model is
consistently underestimating the target, as it is unable to capture additive
effects of mutations.

Figure 3: A subset of melting temperature prediction and corresponding
uncertainties on the test-set when model is trained only on the sigle-
mutations $N{\textrm{train}}=80$

##  5 Conclusions

Utilization of machine learning driven techniques for navigating the drug
design process towards optimized properties requires good compact
representation of the molecule space of interest. We have demonstrated that
sparsity promoting concentration of the larger pre-trained latent space
provided by the protein language model, ESM-1b, leads to more robust estimate
of a dedicated thermal stability optimization task for scFvs. We have proposed
three variants of sparsity promoting effects through GP regression models
integrating learned masks. In general all three models leads to improved
predictive performance relative to a standard mean pooled feature
representation. Even though our sparsity promoting models outperform our
baseline model without mask on our final test data, our validation data
indicates that the learned masks are sensitive to too aggressive sparsity when
validation data is out-of-distribution of the training data. In fact, this
makes sense as the sparse models exactly will seek to favor sparse
representations given the training data at hand. Thus, if we seek to utilize
the models solely for optimization in regions (residue positions) outside the
support of previous seen data, care should be taken in utilizing the masks.
From a Bayesian Optimization perspective, this type of evaluation would
correspond to the explorative evaluation and thus the mean prediction
evaluation is not suitable here. Instead of providing the mean prediction for
evaluation we would seek opportunities to enrich the model support towards new
regions e.g. through the upper-confidence-bound as acquisition function
evaluation. Future work will examine the applicability sparsity promoting
models in the context of steering explorative searches.

## Acknowledgments and Disclosure of Funding

This project was supported in part by the Innovation Fund Denmark via the
Industrial PhD Programme (grant no. 0153-00007B), Novo Nordisk R&D STAR
Fellowship Programme and the Novo Nordisk Foundation (grant no.
NNF20OC0062606) via the Center for Basic Machine Learning Research in Life
Science (MLLS, %\newpagehttps://www.mlls.dk).

## References

  * [1] Roshan Rao, Nicholas Bhattacharya, Neil Thomas, Yan Duan, Peter Chen, John Canny, Pieter Abbeel, and Yun Song.  Evaluating protein transfer learning with tape.  Advances in neural information processing systems, 32, 2019.
  * [2] Alexander Rives, Joshua Meier, Tom Sercu, Siddharth Goyal, Zeming Lin, Jason Liu, Demi Guo, Myle Ott, C. Lawrence Zitnick, Jerry Ma, and Rob Fergus.  Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences.  bioRxiv, 2019.
  * [3] Jesse D Bloom, Sy T Labthavikul, Christopher R Otey, and Frances H Arnold.  Protein stability promotes evolvability.  Proceedings of the National Academy of Sciences, 103(15):5869–5874, 2006.
  * [4] Simona Cocco, Christoph Feinauer, Matteo Figliuzzi, Rémi Monasson, and Martin Weigt.  Inverse statistical physics of protein sequences: a key issues review.  Reports on Progress in Physics, 81(3):032601, 2018.
  * [5] Jiayi Dou, Lindsey Doyle, Per Jr Greisen, Alberto Schena, Hahnbeom Park, Kai Johnsson, Barry L Stoddard, and David Baker.  Sampling and energy evaluation challenges in ligand binding protein design.  Protein Science, 26(12):2426–2437, 2017.
  * [6] Kevin K Yang, Zachary Wu, and Frances H Arnold.  Machine-learning-guided directed evolution for protein engineering.  Nature methods, 16(8):687–694, 2019.
  * [7] Niles A Pierce and Erik Winfree.  Protein design is np-hard.  Protein engineering, 15(10):779–782, 2002.
  * [8] Nobuhiko Tokuriki and Dan S Tawfik.  Stability effects of mutations and protein evolvability.  Current opinion in structural biology, 19(5):596–604, 2009.
  * [9] Claire N Bedbrook, Kevin K Yang, J Elliott Robinson, Elisha D Mackey, Viviana Gradinaru, and Frances H Arnold.  Machine learning-guided channelrhodopsin engineering enables minimally invasive optogenetics.  Nature methods, 16(11):1176–1184, 2019.
  * [10] Florian Richter, Andrew Leaver-Fay, Sagar D. Khare, Sinisa Bjelic, and David Baker.  De novo enzyme design using rosetta3.  PLOS ONE, 6(5):1–12, 05 2011.
  * [11] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer, 2003\.
  * [12] Ethan C Alley, Grigory Khimulya, Surojit Biswas, Mohammed AlQuraishi, and George M Church.  Unified rational protein engineering with sequence-based deep representation learning.  Nature methods, 16(12):1315–1322, 2019.
  * [13] Ahmed Elnaggar, Michael Heinzinger, Christian Dallago, Ghalia Rehawi, Yu Wang, Llion Jones, Tom Gibbs, Tamas Feher, Christoph Angerer, Martin Steinegger, et al.  Prottrans: towards cracking the language of lifes code through self-supervised deep learning and high performance computing.  IEEE transactions on pattern analysis and machine intelligence, 2021\.
  * [14] Jesse Vig, Ali Madani, Lav R Varshney, Caiming Xiong, Richard Socher, and Nazneen Fatema Rajani.  Bertology meets biology: interpreting attention in protein language models.  arXiv preprint arXiv:2006.15222, 2020.
  * [15] Roshan Rao, Joshua Meier, Tom Sercu, Sergey Ovchinnikov, and Alexander Rives.  Transformer protein language models are unsupervised structure learners.  Biorxiv, 2020.
  * [16] Robert E Bird, Karl D Hardman, James W Jacobson, Syd Johnson, Bennett M Kaufman, Shwu-Maan Lee, Timothy Lee, Sharon H Pope, Gary S Riordan, and Marc Whitlow.  Single-chain antigen-binding proteins.  Science, 242(4877):423–426, 1988.
  * [17] Anthony L Fink.  Protein aggregation: folding aggregates, inclusion bodies and amyloid.  Folding and design, 3(1):R9–R23, 1998.
  * [18] Wyatt Strutz.  Exploring protein stability by nanodsf.  Biophysical Journal, 110(3):393a, 2016.

##  Appendix A Appendix

Figure 4: Histogram of the mutations occurring at the respective positions in
training and test set. Figure 5: Histogram of the single site mutations of the
training set. Figure 6: Complete test-set melting temperature prediction and
corresponding uncertainties. Figure 7: Complete test-set melting temperature
prediction and corresponding uncertainties.


page1

[ Supplementary Material

] Evaluation metrics AUROC is the measure of the area under the plot of true
positive rate (TPR) vs. false positive rate (FPR). It measures the probability
that an OOD example is assigned a lower score than an ID example. TPR is
computed as $\mathrm{TPR}=\mathrm{TP}/(\mathrm{TP}+\mathrm{FN})$ and FPR is
computed as $\mathrm{FPR}=\mathrm{FP}/(\mathrm{FP}+\mathrm{TN})$, where TP,
FP, TN, and FN denote True Positive, False Positive, True Negative, and False
Negative, respectively.

AUPR is the Area under the Precision-Recall curve. The PR curve is a graph
showing the $precision=\mathrm{TP}/(\mathrm{TP}+\mathrm{FP})$ and
$recall=\mathrm{TP}/(\mathrm{TP}+\mathrm{FN})$ against each other. The metric
AUPR-In and AUPR-Out are the area under the precision-recall curve where ID
and OOD samples are specified as positives, respectively.

Detection Error $D_{e}$ measures the misclassification probability when TPR is
95%. The formula to calculate $D_{e}=0.5(1-\mathrm{TPR})+0.5\mathrm{FPR}$,
where we assume that both positive and negative examples have equal
probability of appearing in the test set.

FPR @ 95% TPR is the probability that an OOD (negative) example is correctly
identified when the TPR is 95%. Additional results

###  Complete results

tab:full_results_all shows the additional results with all evaluation metrics
and additional layers for all common OOD objects (i.e., OOD objects that are
detected by both the base model and the contrastive model for a given run) for
all OOD datasets. Results separated for each OOD dataset are in
tab:full_results_carla,tab:full_results_ignored,tab:full_results_fp,tab:full_results_waymo.
We highlight two additional observations. For Mahalanobis distance and
normalizing flows, backbone layer provides best overall performance. However,
the optimal layer can vary for different type of OOD objects. For instance,
for KITTI ignored and Waymo objects
(tab:full_results_ignored,tab:full_results_waymo), using normalizing flows
with conv4x features gives better results than with backbone features, whereas
for KITTI FP objects (tab:full_results_fp), backbone layer has 20% higher
AUROC and AUPR compared to conv4x layer. For Mahalanobis distance, the
variation is smaller. Conv2x layer provides better performance for a few
metrics than backbone layer for Carla and KITTI FP objects
(tab:full_results_carla,tab:full_results_fp). For the contrastive model, we
also utilize the contrastive layer features (labeled ‘ctrst’ in the table) to
train Mahalanobis distance, OC-SVM, and normalizing flows. The contrastive
layer does not provide a consistent overall improvement for any of the methods
over all types of OOD objects (tab:full_results_all). Interestingly, for KITTI
ignored objects (tab:full_results_ignored), the contrastive layer is the best
performing layer for all three methods, outperforming other layers by a large
margin. This suggests that the contrastive features are highly biased towards
certain types of OOD objects. It would be interesting to investigate whether
or not this issue also exists for image-based object detection and
classification.

###  Performance impact of OOD detection

To demonstrate how OOD detection would impact the object detector’s
performance during deployment, we evaluate the mAP performance of the base
model on KITTI moderate objects after removing predictions with OOD scores
higher than a threshold. This is done across multiple OOD score thresholds. We
choose normalizing flows (RealNVP) with the backbone layer as our OOD
detector, which has the best overall performance for the base model. We use
one set of experiments with the dataset augmented with Carla, KITTI FP, KITTI
Ignored, and Waymo objects. fig:performance_impact shows the mAP, number of
FP, and number of OOD across multiple OOD score thresholds. With a lower OOD
score threshold, more predictions are marked as OOD objects, which naturally
eliminates more FP and OOD objects but introduces more FN at the same time,
leading to a decrease in mAP performance. From fig:performance_impact, we see
that with 2% mAP performance decrease, over 50% of the FP can be eliminated.
However, the OOD recall varies significantly for different types of OOD
objects. Overall, 35.3% of the objects labeled as OOD can be successfully
identified, which shows that the OOD detection method is a practical addition
to a 3D object detector in deployment.

[b]0.49 [width=]fig/performance_impact_carla.pdf [b]0.49
[width=]fig/performance_impact_ignored.pdf [b]0.49
[width=]fig/performance_impact_fp.pdf [b]0.49
[width=]fig/performance_impact_waymo.pdf

Figure : Carla Figure : KITTI Ignored Figure : KITTI FP Figure : Waymo Figure
: Performance impact of OOD detection.

###  Qualitative results

fig:qualitative shows the qualitative results for OOD detection with max-
softmax score and normalizing flows OOD detection methods. We use the base
model without contrastive learning and normalizing flows with the backbone
layer, which has the best overall performance for the base model
(tab:full_results_all). For each type of OOD objects, we show three examples:
successful detection, normalizing flows (RealNVP) failure, and max-softmax
failure. Successful detections are examples where the OOD object is assigned
low max-softmax score and high OOD score; RealNVP failures are examples where
the OOD score from the normalizing flows for the OOD object is lower than ID
objects; Max-softmax failures are examples where the max-softmax score for OOD
object is higher than ID objects.

[width=0.8]fig/qualitative.pdf

Figure : Qualitative results for OOD detection with max-softmax and OOD scores
from normalizing flows (RealNVP) with backbone layer. Red boxes represents
ground truth ID objects, yellow boxes are OOD objects, and green boxes are
predictions from the base model. The text labels for RealNVP and max-softmax
failures are red.

##  1 Dataset description

We present alg:cap to generate an OOD dataset for a model from an OOD object
database and an ID dataset. It takes an ID dataset, an OOD object database, an
object detection model, and two parameters: maximum number of trials
$\gamma_{max}$ and maximum number of objects $\zeta_{max}$ as inputs. The
output is an OOD dataset that includes ID objects from the ID dataset and OOD
objects from the OOD object database. We also provide a diagram in
fig:ood_gen_diagram describing the OOD dataset generation process. [H] OOD
dataset generation [1] ID dataset ($I$), OOD object database ($O$), detection
model ($M$), maximum number of attempt $\gamma_{max}$, maximum number of OOD
objects ($\zeta_{max}$) OOD dataset ($I$) OODGenerator
$I$,$O$,$M$,$\gamma_{max}$, $\zeta_{max}$ $c\in classes(O$) $\zeta\leftarrow
0$ $f\in I$ and $\zeta\leq\zeta_{max}$ $P\leftarrow M(f)$ obtain predictions
$\gamma\leftarrow 0$ $\gamma\leq\gamma_{max}$ $o\leftarrow rand\\_obj(c)$
select an object $\in c$ $a\leftarrow rand\\_angle()$select an relative angle
$f^{\prime}\leftarrow insert(f,o,a)$ add object to frame IoU($o$, GT($f$),
$P$) $\neq 0$ or $o$ is not in the FOV range
GT returns ground truth objects $\gamma\leftarrow\gamma+1$ continue
$P^{\prime}\leftarrow\textit{M}(f^{\prime})$ obtain predictions $o\notin
P^{\prime}$ $\gamma\leftarrow\gamma+1$ continue $I\leftarrow
update(f^{\prime},I)$ update the dataset $\gamma\leftarrow\gamma+1$
$\zeta\leftarrow\zeta+1$ break return $I$ We utilize different sources to
gather OOD objects to be inserted into the KITTI dataset
[geiger2013kittidataset]. For synthetic objects, we use the Carla simulation
[dosovitskiy2017carla]. We run LiDAR simulation to obtain the point cloud of
objects at various angles and distances. For real objects, we use the KITTI
ignored objects, the KITTI False Positive (FP) objects, and weird vehicle
objects from the Waymo dataset [waymo]. For Carla objects, we set the
intensity values for the points to the median of the KITTI intensities. For
Waymo objects, we transform the original intensities using $\tanh$ then adjust
the intensity values so that the mean and variance under log scale is the same
as the KITTI intensities under log scale. The KITTI FP objects are background
objects classified as Pedestrians. We manually label and categorize them into
classes such as potted plant, bike rack, low traffic sign, sidewalk sign, and
thin sign. Furthermore, we verify if KITTI FP objects appear in the KITTI
training dataset. The tab:fp-in-training shows that except bike rack, other
KITTI FP objects do appear relatively frequently in the KITTI training
dataset. Therefore, KITTI FP objects are a good candidate for the type
\small5⃝ OOD objects. We manually identify objects among the Waymo vehicle
class that are not FG in KITTI, including motorcycle, scooter, digger, and
excavator. To ease the annotation procedure, we first cluster all collected
vehicles in the Waymo validation dataset into either a main cluster or an
outlier cluster using the DBSCAN algorithm and only annotate the vehicles
clustered into the outlier cluster. We run four separate DBSCAN clusterings on
the following four sets of features:
$\text{PCA}([l,w,h,l/w,l/h,w/h,lw,lh,wh],\text{ndim}=2)$, $[l,w]$, $[l,h]$,
and $[w,h]$, where $l$, $w$, and $h$ correspond to the length, width, and
height of the collected vehicle, respectively. We then define the global
outlier cluster as
$\\{\text{vehicles}v|v\text{isanoutlierinatleastoneclustering}\\}$. Next, we
manually collect the vehicles in the global outlier cluster that are likely to
be OOD objects for the KITTI dataset and assign a class label to each
collected out-of-distribution vehicle. Finally, we collect the LiDAR point
cloud of each annotated outlier Waymo vehicle across multiple frames (every
ten frames) in the scene where the vehicle exists. In tab:object_statistics,
we list the number of objects per class from each OOD database that we select
to create OOD datasets. We select these classes because, based on our
experiment, they are likely to be misclassified as foreground objects when we
insert them into the KITTI validation dataset. For each class, we insert 300
objects into an OOD dataset. The number of attempted frames is the number of
frames the algorithm must go through to insert 300 objects. In each frame, the
maximum number of trials is 100. If the number of trials exceeds the maximum
number of trials, we abandon the frame and move on to the next frame.
tab:object_statistics also shows the statistics of inserting objects from each
database to create OOD datasets. Injection failures and detection failures are
the numbers of failed trials averaged over the number of attempted frames. The
higher injection failures, the higher chance that an object overlaps with
existing objects in a frame. The higher the detection failures, the lower
likelihood that an object is detected when we insert it into a frame. For
example, in Carla OOD datasets, bench and swing couch have high injection
failures because they are large objects. Thus, they easily overlap with
existing objects in a frame. These objects also have high detection failures
because a model cannot easily detect them. Furthermore, by comparing the
statistics among datasets, we can conclude that it is easier to insert Waymo
and KITTI FP objects than to insert Carla objects. It makes sense because
LiDAR sensors in both Waymo and KITTI datasets have a similar specification.
In contrast, the Carla objects are synthetic data from a simulation. In
addition, it is easier to insert Carla objects using the contrastive model
than the base model. The observation is opposite for the Waymo and KITTI FP
objects.


# Inverse soft limit construction of QCD amplitude

Andriniaina Narindra Rasoanaivo Sciences Expérimentales et des Mathématiques,
Ecole Normale Supérieure Ampefiloha,
Antananarivo - BP 881, Université d’Antananarivo, Madagascar

###### Abstract

Abstract: QCD amplitudes are one of the most important ingredients for the
understanding of the early universe. In this work we present how the knowledge
of the asymptotic states can be used to calculate the scattering amplitude of
the underline QCD events at high energy. To do so, we show how the
understanding of the soft factorization and the soft structure of amplitudes
can provide us all the necessary information to fix the kinematics of
amplitudes using inverse soft limit techniques.

##  I Introduction

Quantum field theory is known to be the main framework to study the behaviour
of elementary particles and to describe the fundamental interactions between
particles. The computation of scattering amplitudes is the central focus to
connect such theoretical descriptions to the experimental measurements. The
standard method to compute scattering amplitudes is the Feynman diagrammatic
approach, however with this technique the computation tends to be very
complicated as the number of particles involved in the process get large. In
quantum chromodynamic (QCD), the number of gluon tends to grow very fast in
the process which makes the computation of amplitudes so challenging to match
with the high precision experimental results.

Many progresses have been made toward a more efficient computation, the focus
point of such progress is to understand the mathematical structure of
scattering amplitudes[]. The standard approach of Feynman is known to make the
locality and unitarity structure of the amplitude manifest at every step of
the computation. The idea to make both locality and unitarity to manifest at
same time is known to be source of all the complexity in the standard
computation. A modern approach, like the BCFW recursion of Brito Cachazo Fen
Witten Britto:2004ap , only make the unitarity manifest along side the on-
shell property of the particles involved in the process. Despite the non
conventional variables (onshell variables) used in the recursive approach, the
computation turns out very efficient to compute a large number point
amplitudes. Modern approaches also make various properties and mathematical
structure of amplitudes manifest and easy to study.

From our previous work Rasoanaivo:2022kkm , we showed how scattering
amplitudes can be decomposed around the softness of one of its particle which
is an interesting structure around the soft limit theorem Weinberg:1965nx ;
Jackiw:1968zza ; Elvang:2016qvq . The two components of such decomposition are
first the soft core, which is the main subject to the universal Weinberg soft
factorisation, and second the soft shell that needs to be computed in order to
fully reconstruct a lager point amplitude from inverse soft limit (ISL) of a
lower point amplitudes. Many works have been done toward building full
amplitudes from ISL, and such construction can be fully fixed from the
asymptotic symmetry of the individual particles involved in process.

In this work, our main objective is to investigate the possibility of
construction QCD amplitudes from ISL. We will show from our experiences of
soft decomposition of known amplitudes an algorithmic patterns to reconstruct
bipolar amplitudes, a configuration class of amplitudes, from ISL. Such
algorithm will be presented as graph, oriented graph, showing the different
operation to reconstruct full amplitudes.

##  II Inverse soft limit for MHV amplitudes

Before we talk about the inverse soft limit, it is important to talk about the
Weinberg theorem Weinberg:1965nx . This theorem shows that a scattering
amplitude has an universal factorisation behaviour as the momentum of a
massless boson tends to be near zero, Rasoanaivo:2020yii ; Campoleoni:2017mbt
:

$\lim_{k\to
0}A_{n+1}(p_{1},p_{2},\ldots,p_{n},k)=\hat{S}(k)A_{n}(p_{1},p_{2},\ldots,p_{n}),$
(1)

where $A_{n}$ is an $n$-point amplitude and $\hat{S}(k)$ is the soft operator
related to the soft momentum $k$. This factorisation is known to be universal,
and its allow us to connects an $(n+1)$-points amplitude to an $(n)$-point. As
shown in Rasoanaivo:2020yii  the soft operator $\hat{S}$ could be derived
independently from the amplitude by resolving the following equation

$\left[\hat{H}_{i},\hat{S}(k_{j})\right]=h_{i}\hat{S}(k_{j})\delta_{ij}.$ (2)

In this equation, $\hat{S}(k_{j})$ is the soft operator associated to $j$-th
particle taken to be soft, $h_{i}$ is the helicity of the $i$-th particle,
$\delta_{i}j$ is the Kronecker symbol, and $\hat{H}_{i}$ is the helicity
operator associated to $i$-th particle where with the spinor helicity
variables Rasoanaivo:2020yii

$\hat{H}_{i}=-\frac{1}{2}\left(\lambda^{a}_{i}\frac{\partial\;}{\partial\lambda^{a}_{i}}-\bar{\lambda}^{\dot{a}}_{i}\frac{\partial\;}{\partial\bar{\lambda}^{\dot{a}}_{i}}\right).$
(3)

The inverse soft limit is the opposite action where a soft operator is applied
on lower point amplitude while its momentum is taken to be hard so the the
$n$-point will give rise to an $(n+1)$-point. However as shown in
Rasoanaivo:2022kkm , such process is not possible with one single ISL. The
soft limit can be represented as

$A_{n+1}\xrightarrow[SL]{\;\text{soft limit}\;}A_{n}$ (4)

while the inverse soft

$A_{n}\xrightarrow[ISL]{\;\text{invert soft
limit}\;}A_{n+1}^{ISL}=A_{n+1}-R_{n+1},$ (5)

where $R_{n+1}$ is the missing part of the higher point amplitude that can be
completed from an ISL of other soft sector provided by a different particle.
In this soft decomposition, the main objective of the ISL reconstruction is to
find an algorithm to calculate $R_{n}$ so we can complete the amplitude.
Around this idea of soft reconstruction, it is worth explore the decomposition
of an amplitude around the soft behaviour of a given particle. Let us consider
the $i$-th particle in the process, an amplitude can be decompose as the core
of the soft theorem of the $i$-th particle $A^{[i]}_{n}$(the part that can be
recovered by ISL), and the soft shell which is lost in this soft limit
$R_{n}^{[i]}$,

$A_{n}=A^{[i]}_{n}+R^{[i]}_{n}\Longleftrightarrow\left\\{\begin{aligned}
&A_{n}\xrightarrow[]{\;k_{i}\to 0\;}A_{n-1}\xrightarrow{\;ISL\>}A_{n}^{[i]}\\\
&R^{[i]}_{n}\xrightarrow[]{\;k_{i}\to 0\;}0.\end{aligned}\right.$ (6)

One of the main problem of the standard approach is that amplitudes are
computed in a more generalised way where the helicity of particles are only
specified at the very last step of the computation. In the modern approaches,
especially after Park and Taylor introduce formula for some $n$-points
amplitude Parke:1986gb , the computation of amplitudes simplified depending on
them helicity configuration.

Figure 1: Six-point amplitude from one MHV and one $\overline{\text{MHV}}$
amplitudes

In order to understand better these classification let us only consider pure
gluon amplitudes were all the particles are all outgoing. Depending on the
helicity configuration, for any $n$-point amplitudes, here are the
classification of gluon amplitudes

  * •

Vanishing amplitudes: are mplitude vanishes when all the particles have the
same helicity or there is only one negative of positive helicity.

  * •

MHV amplitudes: these amplitudes are the first non vanishing amplitudes in
most non conserved way of the helicity in which we only have two negative
helicity particles.

  * •

$\overline{\text{MHV}}$ amplitudes: the MHV bar amplitudes are similar to the
MHV amplitudes in which we only have two positive helicity particles.

  * •

N${}^{k}\\!$MHV amplitudes: are amplitudes that have $k+2$ negative helicity
particles, $k$-th next the MHV amplitudes.

The well understood among these amplitudes are the MHV and the
$\overline{\text{MHV}}$ amplitudes. Them expression can be derived directly
from the Park-Taylor formula Parke:1986gb . And the amazing part is that in a
soft limit where $A_{n+1}^{MHV}$ amplitudes is factorized into a soft operator
and a MHV sub-amplitude $A_{n}^{MHV}$ a single inverse soft limit is enough to
recover the original amplitude,

$A_{n+1}^{MHV}\xrightarrow{\;\text{ soft limit
}\;}A_{n}^{MHV}\xrightarrow{\;\text{ inverse soft limit }\;}A_{n+1}^{MHV}.$
(7)

This stability of the MHV and $\overline{\text{MHV}}$ amplitudes and them
simplicity are the main reasons we choose MHV amplitude to be central part of
amplitudes. For example the 6-point amplitude
$A_{6}(1^{-}2^{-}3^{-}4^{+}5^{+}6^{+})$ [2205.13873], it can be constructed
for two inverse soft limits as shown by the graph in the Fig. 1.

##  III Bipolar amplitudes

Figure 2: Seven-point bipolar amplitude reconstruction.

In our experiments, we analysed the soft decomposition of many different
helicity amplitudes. From these decomposition, we explore the possibility of
inverting the process using ISL calculations. It is important to mention that
we are in the search of pattern to generalise the ISL reconstruction of any
point amplitude.

Out of these experiments, we found that there seems a pattern in some class of
helicity configurations. These amplitude have all of its positive helicity
grouped in one region. A typical form of that is when all the positive
helicity is grouped to the right and all the negative to left such as
$A_{6}(1^{-}2^{-}3^{+}4^{+}5^{+})$, $A_{6}(1^{-}2^{-}3^{-}4^{+}5^{+}6^{+})$,
$A_{6}(1^{-}2^{-}3^{-}4^{+}5^{+}6^{+}7^{+})$, …Here we can see that a bipolar
amplitude can ever MHV or any of the NkMHV amplitudes since we do not care
about the number the negative helicity but more on how are they grouped as
bipolar objects.

For simplicity we are calling these amplitude bipolar amplitude. In the soft
decomposition experiments we found that any $n$-point bipolar amplitude can be
constructed directly from two ISL of some other $(n-1)$-points bipolar. This
finding allow to generate an ISL algorithm to reconstruct those amplitude from
MHV and/or $\overline{\text{MHV}}$ amplitudes. The algorithm can be
represented in a bi-parted graph in which:

  * •

the root: is the bipolar amplitude that we aim to reconstruct.

  * •

the leaves: are MHV/$\overline{\text{MHV}}$ amplitudes which are the seeds of
the construction.

  * •

the oriented edges: are the ISL operations that connect the $n$-points to the
$(n+1)$-points.

  * •

the intermediate vertices: are the intermediate bipolar amplitude that
connects the MHV’s to the root.

The Fig. 2 and Fig. 3 are two examples of the calculation we did using the bi-
parted algorithm above. It is worth to mention that the ISL operation made at
every edge in the graph consist to multiply the amplitude with the
corresponding soft operator of the particle labelled on the arrow and to use a
momentum shift as presented in Boucher-Veronneau:2011rwd ; Nandan:2012rk

Figure 3: Eight-point bipolar amplitude reconstruction.

##  IV Conclusion

In conclusion, we aim to construct QCD amplitudes from the asymptotic
behaviour of the individual particles using inverse soft limit. And we found
that such process is well establish for the case of MHV amplitudes since they
have a simplicity in the soft decomposition. Such decomposition is main
ingredient of this reconstruction, and the analysis of decomposition of many
different amplitudes allow us to observe a pattern in the soft structure
certain amplitudes we called bipolar amplitudes. In the present work presented
at the HEPMAD22, we only gave the ISL algorithm, the formal proof will be
presented in a separate paper of our coming works. The extension of the method
can be considered for any helicity configuration, the only complication of
some helicity configuration is the apparition of loops in the graph associated
to the ISL algorithm.

## References

  * (1)

## REFERENCES

* (2) R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B 715, 499-522 (2005) [arXiv:hep-th/0412308 [hep-th]].
* (3) A. N. Rasoanaivo, [arXiv:2205.13873 [hep-th]].
* (4) S. Weinberg, Phys. Rev. 140, B516-B524 (1965)
* (5) R. Jackiw, Phys. Rev. 168, 1623-1633 (1968)
* (6) H. Elvang, C. R. T. Jones and S. G. Naculich, Phys. Rev. Lett. 118, no.23, 231601 (2017) [arXiv:1611.07534 [hep-th]].
* (7) A. N. Rasoanaivo, [arXiv:2002.02120 [hep-th]].
* (8) A. Campoleoni, D. Francia and C. Heissenberg, JHEP 05, 120 (2017) [arXiv:1703.01351 [hep-th]].
* (9) S. J. Parke and T. R. Taylor, Phys. Rev. Lett. 56, 2459 (1986)
* (10) C. Boucher-Veronneau and A. J. Larkoski, JHEP 09, 130 (2011) [arXiv:1108.5385 [hep-th]].
* (11) D. Nandan and C. Wen, JHEP 08, 040 (2012) [arXiv:1204.4841 [hep-th]].