# Spiking Neural Networks with Dynamic Time Steps for Vision Transformers

Gourav Datta∗, Zeyu Liu∗, Anni Li, and Peter A. Beerel Ming Hsieh Dept. of
Electrical and Computer Engineering, University of Southern California, Los
Angeles, USA
{gdatta, liuzeyu, annili<EMAIL_ADDRESS>
∗Equally contributing authors.†This work was supported by a gift funding from
Intel.

###### Abstract

Spiking Neural Networks (SNNs) have emerged as a popular spatio-temporal
computing paradigm for complex vision tasks. Recently proposed SNN training
algorithms have significantly reduced the number of time steps (down to 1) for
improved latency and energy efficiency, however, they target only
convolutional neural networks (CNN). These algorithms, when applied on the
recently spotlighted vision transformers (ViT), either require a large number
of time steps or fail to converge. Based on analysis of the histograms of the
ANN and SNN activation maps, we hypothesize that each ViT block has a
different sensitivity to the number of time steps. We propose a novel training
framework that dynamically allocates the number of time steps to each ViT
module depending on a trainable score assigned to each timestep. In
particular, we generate a scalar binary time step mask that filters spikes
emitted by each neuron in a leaky-integrate-and-fire (LIF) layer. The
resulting SNNs have high activation sparsity and require only accumulate
operations (AC), except for the input embedding layer, in contrast to
expensive multiply-and-accumulates (MAC) needed in traditional ViTs. This
yields significant improvements in energy efficiency. We evaluate our training
framework and resulting SNNs on image recognition tasks including CIFAR10,
CIFAR100, and ImageNet with different ViT architectures. We obtain a test
accuracy of 95.97% with 4.97 time steps with direct encoding on CIFAR10.

##  I Introduction

Spiking neural networks (SNN) have emerged as an attractive edge-based
computing paradigm, thanks to their high activation sparsity and use of cheap
accumulation (AC) instead of energy-expensive multiply-and-accumulation (MAC)
[1, 2]. Early research on SNNs suggested that they can achieve good
performance in a wide range of modern vision tasks [2], but not without a
noticeable accuracy drop compared to the state-of-the-art (SOTA), particularly
when restricted to energy-efficient low-latency SNNs. This motivated a
plethora of research on supervised SNN training algorithms, that can generally
be grouped into i) ANN-to-SNN conversion, and ii) direct training. ANN-to-SNN
conversion leverages standard backpropagation for training in the ANN domain,
and yields SNNs by approximating the non-linear activation functions (e.g.
ReLU in CNNs) with the firing rate of spiking neurons [3, 2, 4, 5]. In
contrast, direct training algorithms typically leverage surrogate gradient
learning that approximates the discontinuous spiking neuron functionality with
a continuous differentiable model [6, 7, 8, 9]. While direct training
algorithms require the gradients to be integrated over all the SNN time steps,
thereby increasing the training complexity, it requires far less number of
time steps compared to ANN-to-SNN conversion. To obtain the best of both
worlds, a few research works [10, 11, 12, 13] attempted to perform ANN-to-SNN
conversion, followed by fine-tuning in the SNN domain. The latest SNN
literature [14, 15, 16] suggests that direct training is possible with only a
few time steps ($1$ or $2$) and with a slight increase in training complexity
compared to non-spiking networks. These works significantly decrease latency
and energy of SNNs while not compromising accuracy, thereby enabling the
deployment of SNNs at the extreme edge.

In spite of these tremendous improvements in SNN training algorithms, most (if
not all) existing SNN models are based on CNN backbones, such as VGG [17] and
ResNet [18]. On the other hand, the last two years have witnessed the
dominance of vision transformers (ViT) [19] that have outperformed SOTA CNNs
in complex vision tasks. There is little research on SNNs for ViTs,
particularly for low time steps and edge applications, where SNNs are the most
attractive candidates for edge efficiency. Existing works propose several
architectural modifications in ViTs for SNNs that incur a large number of time
steps ($4$ and above) to obtain similar accuracy as the non-spiking
counterparts [20, 21]. There has not been any thorough analysis on the spiking
activation functions for ViTs (in particular the self-attention block that is
unique to ViTs) or training optimizations that reduce the number of time
steps.

With this motivation, this paper first shows that the existing CNN-based low-
latency SNN training approaches cannot be ported to ViTs, without a
significant drop in accuracy, particularly at low time steps, implying there
are large gaps in both the accuracy and efficiency advantages of the
traditional vision and SNN models. In order to mitigate these gaps, we
investigate the key sources of information loss in ViT-based SNNs and propose
a novel training algorithm that allocates time steps dynamically to each ViT
module based on its precision sensitivity. Our algorithm yields SNN models
that can obtain close to SOTA accuracy with an average number of time steps as
low as ${\sim}2.73$ on the CIFAR100 dataset.

##  II Preliminaries

###  II-A LIF Model & Surrogate gradient

In this work, we adopt the popular leaky-integrate-and-fire (LIF) neuron [22]
models to capture the computation dynamics of an SNN, where each neuron in the
$l^{th}$ layer has an internal state called its membrane potential $U_{l}^{t}$
($t$ denotes the time step). It captures the sum of the weight $W_{l}$
modulated incoming spikes $O_{l{-}1}^{t}$ from each neuron in the previous
layer $(l{-}1)$ as shown below.

$\displaystyle H_{l}^{t}$
$\displaystyle=\lambda_{l}U_{l}^{t{-}1}+W_{l}O_{l{-}1}^{t}\ \ O_{l}^{t}$
$\displaystyle=\begin{cases}V_{l}^{th},&\text{if }H_{l}^{t}>V_{l}^{th}\\\
0,&\text{otherwise}\\\ \end{cases}$ (1) $\displaystyle U_{l}^{t}$
$\displaystyle=H_{l}^{t}-O_{l}^{t}$ (2)

Note that $U_{l}^{t}$ leaks with a fixed time constant, denoted as
$\lambda_{l}$. Since the spiking neuron function $O_{l}^{t}$ is discontinuous
and non-differentiable, it is difficult to implement gradient descent based
backpropagation in SNNs. Hence, previous works [7, 23] approximate the spiking
function with a continuous differentiable function, which helps back-propagate
non-zero gradients known as surrogate gradients. The resulting weight update
in the $l^{th}$ hidden layer in the SNN is calculated as

$\displaystyle\Delta{W_{l}}{=}\sum_{t}\frac{\partial\mathcal{L}}{\partial
W_{l}}{=}\sum_{t}\frac{\partial\mathcal{L}}{\partial\mbox{$O$}_{l}^{t}}\frac{\partial\mbox{$O$}_{l}^{t}}{\partial\mbox{$H$}_{l}^{t}}\frac{\partial\mbox{$H$}_{l}^{t}}{\partial
W_{l}}{=}\sum_{t}\frac{\partial\mathcal{L}}{\partial\mbox{$O$}_{l}^{t}}\frac{\partial\mbox{$O$}_{l}^{t}}{\partial\mbox{$H$}_{l}^{t}}\mbox{$O$}_{l-1}^{t}$

where $\frac{\partial{O}_{l}^{t}}{\partial{U}_{l}^{t}}$ is the non-
differentiable gradient which can be approximated with the surrogate gradient
$\frac{\gamma}{V_{l}^{th}}\cdot{max(0,1-\text{abs}(\mbox{$\frac{U_{l}^{t}}{V_{l}^{th}}-1$))}}$,
where $\gamma$ is a hyperparameter denoting the maximum gradient value.

###  II-B Architectural Modifications for Spiking ViT

Inspired by [20], we adopt some architectural modifications to the standard
ViT backbone to simplify the training in the SNN domain. Firstly, we remove
the Dropout and Droppath present in standard ViTs. Secondly, we remove the
layer norm before each self-attention and multi-layer perceptron (MLP) block,
and instead, add batch norm (BN) after each linear layer. Thirdly, we remove
the softmax layer in the multi-head self-attention (MSA) block since the
attention maps computed from the spiking query and key outputs are always
positive. This is because we use a spiking LIF layer before every dot-product
operation incurred in convolutional/linear/MSA layers so that we require only
accumulate operations for energy-efficiency, as promised in SNNs. Details of
our spiking ViT architecture are further illustrated in Fig. 1.

##  III Proposed Method

Based on the LIF model, we propose a Dynamic Time Step Spiking layer (DTSS
layer) that adaptively assigns different time steps to different LIF layers of
the spiking ViT. Intuitively, this scheme can skip the computations and
latency incurred in processing timesteps that can be potentially unnecessary
in some layers. Furthermore, we analyzed the distribution of SNN activation
maps for each layer of the Spikformer [20] at different timesteps. This
analysis shows that the sensitivity of each transformer layer varies with the
different time steps. Based on this intuition and analysis, we propose a
dynamic time step spiking ViT training framework, that achieves close to SOTA
test accuracy with lower latency.

(a)

(b)

Figure 1: (a) Proposed network architecture of the spiking ViT used in this
work, where PS denotes patch splitting and (b) Proposed dynamic time step
spiking block and the MSA block (dot-product w/ scale) without the softmax
layer.

###  III-A Dynamic Time Steps Spiking Layer

In this subsection, we present our DTSS layer based on the traditional multi-
time-step LIF model. Specifically, we first set a maximum number of time steps
$T_{max}$ for our SNN model, which means that the number of time steps
assigned to each layer is less than or equal to $T_{max}$. We then define a
trainable tensor of length $T_{max}$, called the parameter tensor $TP_{l}$,
and an all-one lower triangular matrix of size $T_{max}{\times}T_{max}$ called
the coefficient matrix $\it{coeff}$. The parameter tensor $TP_{l}$ is
initialized with all non-negative values, so that $TS_{l}$, which we denote as
the score tensor, and which is the product of $TP_{l}$ and $\it{coeff}$ (see
Eq. 5 below), is a non-increasing tensor. This helps ensure that the score
tensor $TS_{l}$ represents a valid time step mask such that we do not skip any
intermediate time steps i.e., we only skip time steps towards the end.

To convert the score tensor at a particular time step $t$ ($TS_{l}^{t}$) to a
binary mask tensor denoted by $TM_{l}^{t}$, we propose a simple spiking layer
as illustrated in Eq.(5). Since it is not possible to update $TP_{l}$ with the
real gradients backpropagated from $TM_{l}^{t}$ (gradients are zero almost
everywhere), we use surrogate gradients. Note that this spiking layer is
different from the LIF layers used in the spiking ViT network architecture. As
illustrated in Fig. 1, $TM_{l}^{t}$ dynamically retains the output of the LIF
neurons of a particular layer in some specific time steps while setting the
output in the remaining time steps to zero. This helps in shutting off the
computations incurred in the subsequent layer for the remaining time steps,
which increases the compute efficiency of the SNN. Moreover, during inference
when $TM_{l}^{t}$ is fixed, this also reduces the number of time steps of the
subsequent layer, which with appropriate spike routing support from the
underlying hardware compiler, can also increase the SNN latency efficiency. To
quantify this latency reduction, we define the total number of time steps of
our SNN as

$T_{avg}=\frac{\sum_{l=1}^{L}T_{l}}{L},l\in[1,2,3,...,L]$ (3)

where $T_{l}$ is the total number of time steps of the $l^{th}$ layer, and $L$
is the total number of layers.

Without loss of generality, let us illustrate the function of our DTSS layer
for the case when $T_{max}{=}4$. Denoting ${O}_{l}^{t}$ as the input and
${W}_{l}$ as the weights of the convolutional/linear layer $l$, we first
obtain the temporary spiking outputs $Y_{l}^{t}$ following the LIF model
illustrated in Eq. (1-3). Let us denote this model as $SN$. We then apply our
mask tensor $TM_{l}^{t}$ to filter out the final spiking outputs ${O}_{l}^{t}$
that would be fed to layer $l{+}1$ only when all the elements of ${O}_{l}^{t}$
are not zero.

$\small{coeff{=}\begin{bmatrix}1&0&0&0\\\ 1&1&0&0\\\ 1&1&1&0\\\
1&1&1&1\end{bmatrix}\;TS_{l}{=}TP_{l}{\cdot}coeff{=}[TS_{l}^{0},TS_{l}^{1},TS_{l}^{2},TS_{l}^{3}]}$
(4)
${TM}_{l}^{t}=\Theta({TS}_{l}^{t}){=}\begin{cases}1,&\text{if}{TS}_{l}^{t}\geq
1;\\\ 0,&\text{otherwise}\end{cases}$ (5)
$\frac{\partial({TM}_{l}^{t})}{\partial({TS}_{l}^{t})}{=}\begin{cases}1,&\text{if
}0<{TS}_{l}^{t}<2;\\\ 0,&\text{otherwise}\end{cases}$ (6)

${Y}_{l}^{t}=SN({W}_{l}{O}_{l{-}1}^{t})\;\;\;\;{O}^{t}_{l}={Y}_{l}^{t}*TM_{l}^{t}$
(7)

###  III-B Distribution of SNN Membrane Potentials

Apparently, the number of time steps in a spiking ViT layer $l$ is the sum of
$TM_{l}^{t}$ across the time dimension $t$, denoted as $TM_{l}$. Therefore, to
reduce the number of time steps, we define the mask loss, which is the sum of
$TM_{l}$ across all the layers. However, our experimental results indicate
that naively applying this loss to each layer, incurs a large drop in the SNN
test accuracy. To mitigate this issue, we analyze the input distribution of
each LIF layer for the spiking ViT model at different time steps. For example,
in Fig. 2, we visualize the distribution map of some representative spiking
ViT layers at the first and the sixth time step with a total of six time
steps.

For any multi-time-steps SNN, each LIF layer distribution map can change
either slightly or significantly across different time steps. We can interpret
the distribution map at the first time step as the baseline, which can change
in the subsequent time steps based on the sensitivity (that represents the
distribution similarity) of the layer with the different time steps. For
example, if the distribution map of layer $l$ in the first time step is
similar to its distribution in all the remaining time steps, it is likely that
the layer $l$ can be assigned only a single time step without a significant
accuracy drop. Furthermore, considering each distribution map as a probability
density vector, we use the cosine similarity metric to quantitatively measure
the sensitivity. Empirically, we observe that the input distribution of the
LIF layers in the MSA, combined across the different time steps, is similar
between different blocks. Based on that, we further assume that the
sensitivity to the different time steps should be linked to certain groups of
layers instead of individual layers. Therefore, we bin all the spiking ViT
layers into four categories, the ”sps” part, the ”qkv” part, the ”attn” part,
and the ”mlp” part, as shown in Fig 2. The ”sps” part represents the layers in
the Spiking Patch Splitting (SPS) module [20]. The ”qkv” part represents the
layers that compute the query, key, and value tensors from the inputs. The
”attn” part represents the QKV matrix dot-product followed by the scaling
operation and the subsequent linear operation after the ”qkv” part in the MSA
block. The ”mlp” part indicates the multi-layer perception block. Each part
has varying sensitivites to the different time steps as shown in Fig. 2, which
controls the decision to skip a particular time step differently.

Figure 2: The distribution of activation maps of the spiking models, the
orange color means the distribution of the one-time-step model while the blue
color means the distribution of the six-time-steps model. The label in each
plot shows the name of the layer. The first five figures are the distributions
of the layers in the spiking patch splitting (SPS) module, the next five are
the distributions of the layers in the spiking self-attention (SSA) block, and
the last two are the distributions of the layers in the multi-layer perception
(MLP) block. To keep it simplified, we omit the distributions of the remaining
layers.

###  III-C Training Loss

Based on the above findings, we compute the loss function of our SNN as the
sum of the cross entropy loss (computed from surrogate gradient descent) and
the mask loss, as illustrated in Eq.(8), where, $l_{s}$ denotes the index of
the layers that are sufficiently sensitive to the different time steps and the
hyperparameter $\lambda_{m}$ helps to control the accuracy-efficiency trade-
off. During the surrogate gradient descent, the weights of the SNN are updated
based on the gradients to the normal SNNs but with a $TP^{t}_{l}$, which is
illustrated in Eq.(9). However, computing the gradients of $TP^{t}_{l}$ is not
straight forward. As illustrated in Eq.(10), its gradient flows from both
$L_{CE}$ and $L_{mask}$ in two separate directions.

$L_{total}=L_{CE}+L_{mask}=L_{CE}+\lambda_{m}\sum_{l_{s}}\sum_{t}{TM}_{l_{s}}^{t}$
(8) $\frac{\partial L_{total}}{\partial W_{l}}=\frac{\partial L_{CE}}{\partial
O_{l}^{t}}\frac{\partial O_{l}^{t}}{\partial W_{l}}={TM}^{t}_{l}\frac{\partial
L_{CE}}{\partial O_{l}^{t}}\frac{\partial Y_{l}^{t}}{W_{l}}$ (9)
$\frac{\partial L_{total}}{\partial TP^{t}_{l}}==Y_{l}^{t}\frac{\partial
L_{CE}}{\partial O_{l}^{t}}\frac{\partial TM_{l}^{t}}{\partial
TP_{l}^{t}}+\frac{\partial\lVert
TM^{t}_{l}\rVert_{1}}{TM^{t}_{l}}\frac{\partial TM^{t}_{l}}{\partial
TP^{t}_{l}}$ (10) $\frac{\partial TM^{t}_{l}}{\partial
TP^{t}_{l}}=\frac{\partial TM^{t}_{l}}{\partial TS^{t}_{l}}\frac{\partial
TS^{t}_{l}}{\partial TP^{t}_{l}}$ (11)

##  IV Experimental Results

###  IV-A Datasets & Models

Our models are based on the compact ViT architecture [24] in which we replace
all the linear layers with convolutional layers and remove their bias
parameters. We performed object recognition experiments on the CIFAR10/100
[25] and ImageNet [26] datasets. The architecture name ”Spikformer-N-D”
indicates that the model is based on the spikformer architecture [20] and has
N self-attention encoder blocks, and the feature embedding dimension is D. For
the CIFAR10/100 dataset, we scale the image resolution to $32{\times}32$ and
use the spikformer-4-384 architecture. For training, we use Adamw [27]
optimizer with an initial learning rate of 0.004, and a batch size of 256. On
the other hand, for the ImageNet dataset, we scale the image resolution to
$224{\times}224$ and use the spikformer-8-512 architecture to evaluate our
approach on large-size models. The initial learning rate is 0.001, the total
batch size is 168, and we train for a total of 300 epochs.

TABLE I: Comparison of our dynamic time-step spiking ViT models to the original spikformer model with different time steps.∗ denotes results based on our replicated setup. Reference | Architecture | Acc. (%) | Time steps
---|---|---|---
Dataset : CIFAR10
CCT[24] | Spikformer-4-384 | 96.61 | 1
Spikformer[20] | Spikformer-4-384 | 94.13 | 1
95.29 | 2
95.66 | 4
95.94 | 6
This work | Spikformer-4-384 | 94.39 | 1.15
95.97 | 4.79
Dataset : CIFAR100
CCT[24] | Spikformer-4-384 | 80.19 | 1
Spikformer[20] | Spikformer-4-384 | 75.05 | 1
77.10 | 2
77.61 | 3
78.14 | 4
78.85 | 6
This work | Spikformer-4-384 | 78.09 | 2.73
78.90 | 3.81
Dataset : ImageNet
Spikformer[20] | Spikformer-8-512 | 63.12* | 1
66.03* | 2
73.38 | 4
This work | Spikformer-8-512 | 68.04 | 1.30

###  IV-B Comparison with SOTA ViT-based SNNs

Since our spiking ViT network architecture is based on Spikformer[20], we
exhaustively compare the SNN accuracy obtained by our framework with
Spikformer for different number of time steps. We observe that Spikformer
incurs a drastic drop in test accuracy as the total number of time steps
decreases from $2$ to $1$. Therefore, it is important to yield a SNN that can
optimize this trade-off between the test accuracy and time steps. As shown in
Table 1, our framework can improve the accuracy by 2.01% in the ImageNet
dataset with a decrease in the total number of time steps by 46% compared to
the one-time-step model. Our SNNs can thus obtain a balanced and controlled
(with the hyperparameter $\lambda_{m}$ shown in Eq. 10) accuracy-efficiency
trade-off compared to the spikformer models with both one and two time steps.
Furthermore, as also shown in Table 1, for the CIFAR10/100 dataset, our
framework can achieve similar accuracy as Spikformer with 6 time steps while
reducing the latency ($T_{avg}$) by about 20% (from 6 to 4.79) and 36.5% (from
6 to 3.81).

###  IV-C Comparison with SOTA CNN-based SNNs

As shown in Table II, the SNNs obtained by our framework outperforms the
existing CNN-based SNNs in terms of latency-accuracy trade-off. Although,
unlike ViTs, the CNN-based SNNs can achieve accuracies similar to the non-
spiking CNNs even at low latencies ($1{-}4$ time steps), they are still
relatively far from the SOTA that ViTs can provide. Our framework can help
achieve this SOTA while still keeping the number of time steps low.

As shown in Table 2, applying the SOTA CNN-based Hoyer regularized training
method [9] (that achieves impressive SNN accuracies with only one time step)
on ViTs, lead to a significant accuracy drop in CIFAR10/100. This is because
it cannot cleanly separate the input distribution of the LIF layers in the MSA
block around the threshold, causing training convergence issues. Other
training techniques, including TET [14], also cannot achieve close to SOTA
accuracies, when applied to ViTs.

TABLE II: Comparison of our dynamic time-step spiking ViT models to existing CNN and ViT-based low-latency SNNs. ∗ denotes results are based on our replicated setup. Reference | Architecture | Acc. (%) | Time steps
---|---|---|---
Dataset : CIFAR10
Diet-SNN[10] | VGG16 | 93.44 | 10
STBP-tdBN[28] | ResNet-19 | 93.16 | 6
Hoyer-reg[9] | VGG16 | 93.44 | 1
DSR[29] | PreAct-ResNet-18 | 95.40 | 20
TET[14] | ResNet-19 | 94.44 | 4
Hoyer-reg[9] | Spikformer-4-384 | 93.69* | 2
TET[14] | Spikformer-4-384 | 94.08* | 2
This work | Spikformer-4-384 | 94.39 | 1.15
95.97 | 4.79
Dataset : CIFAR100
Diet-SNN[10] | VGG16 | 69.67 | 5
STBP-tdBN[28] | ResNet-19 | 70.86 | 4
Hoyer-reg[9] | Spikformer-4-384 | 71.55 | 1
DSR[29] | PreAct-ResNet-18 | 78.50 | 20
TET[14] | ResNet-19 | 74.47 | 4
Hoyer-reg[9] | Spikformer-4-384 | 71.29* | 2
TET[14] | Spikformer-4-384 | 74.88* | 2
This work | Spikformer-4-384 | 78.09 | 2.73
78.90 | 3.81

###  IV-D Ablation study

In this subsection, we demonstrate the impact of mask loss on the accuracy and
the FLOPs of our SNNs through two groups of controlled experiments. One group
is related to the count of different components in the mask loss, while the
other is related to the weight of the mask loss.

The results in Table V validate our hypothesis described in the Section III-B.
All the experiments are trained and tested on the CIFAR100 dataset. The
initial time step is $2$, and the weight of the mask loss is $1e{-}6$. As
shown in the first row of Table 3, naively incorporating the mask loss for all
the layers can significantly reduce the total number of time steps and hence,
the FLOPs, but the accuracy also suffers significantly; the last row indicates
that, without the regularization provided by the mask loss, the model still
has a large number of time steps after training. Considering the specificity
of the SPS module, when we remove the ”sps loss” (see second row), the
accuracy remains similar, however, the total number of time steps reduces by a
factor of $2$. Furthermore, based on our analysis of the sensitivity of layers
to the different time steps, we observe that omitting the mask loss of the
self-attention layers (along with SPS layers) yields the highest accuracy (see
third row), and the number of time steps is reasonably small.

TABLE III: Ablation study results on different mask loss combinations, where TS denotes the total number of time steps, and SA denotes spiking activity. qkv loss |  attn loss |  mlp loss |  sps loss | Acc. (%) | TS | SA(%) | SOPs (M)
---|---|---|---|---|---|---|---
$\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | 74.64 | 1.03 | 2.90 | 101.24
$\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | 77.16 | 1.96 | 4.87 | 181.50
$\checkmark$ | $\times$ | $\checkmark$ | $\times$ | 78.09 | 2.73 | 5.56 | 185.31
$\times$ | $\times$ | $\times$ | $\times$ | 77.10 | 3.93 | 5.33 | 182.69

The experiments in Table IV use the mask loss consisting of ”qkv” loss and
”mlp” loss since this configuration provided the best accuracy. When
initialized with a small initial time step, a smaller weight of mask loss
$\lambda_{m}$ incurs significantly higher spiking activity and FLOPs (see
first two rows). However, when the initial time step increases, the effect of
the $\lambda_{m}$ becomes less pronounced. Note that a too-small initial time
step fails to get high accuracy, and a too-large initial time step fails to
reduce the total number of time steps and hence, the FLOPs significantly. When
both the initial time step and the maximum time step are equal (see the
seventh row), although the total number of time steps do not change after
training, the accuracy improves significantly. This indicates that our mask
loss may help the model find a better local optimization. We empirically
observe initializing the time step to half of the maximum value yields the
best accuracy-efficiency trade-off.

TABLE IV: Ablation study results on the different weights of mask loss and different initial time steps, where $\lambda_{m}$ denotes the weight of mask loss, TS denotes the number of time steps, and SA denotes spiking activity. Init. TS |  Max. TS | $\lambda_{m}$ | Acc. (%) | TS | SA(%) | SOPs (M)
---|---|---|---|---|---|---
1 | 4 | 1e-4 | 74.45 | 2.45 | 3.73 | 172.54
1 | 4 | 1e-6 | 73.90 | 3.0 | 20.97 | 429.18
2 | 4 | 1e-4 | 76.96 | 2.69 | 5.41 | 183.65
2 | 4 | 1e-6 | 78.09 | 2.73 | 5.33 | 185.31
3 | 4 | 1e-4 | 78.16 | 3.36 | 7.70 | 254.39
3 | 4 | 1e-6 | 78.07 | 3.36 | 7.71 | 255.78
4 | 4 | 1e-4 | 78.69 | 4.0 | 9.85 | 319.75
4 | 6 | 1e-4 | 78.90 | 3.81 | 6.54 | 319.51

###  IV-E Energy Efficiency of our Proposed Approach

We compare the energy-efficiency of our dynamic time step SNNs with non-
spiking ViTs and existing spikformer-based SNNs in Table 5. The compute-
efficiency of SNNs stems from two factors: 1) sparsity, that reduces the
number of floating point operations (FLOPs) in convolutional, linear, and MSA
layers compared to non-spiking DNNs according to $SNN^{flops}_{l}=S_{l}\times
DNN^{flops}_{l}$, where $S_{l}$ denotes the average number of spikes per
neuron per inference over all timesteps in layer $l$. Note that the sparsity
induces a small overhead of checking whether the 1-bit activation is zero,
which consumes 0.05pJ in 28nm Kintex-7 FPGA platform according to our post
place-and-route simulations. 2) Use of only AC (1.8pJ) operations that consume
$7.4\times$ lower compared to each MAC (13.32pJ) operation in 45nm CMOS
technology [30] in our FPGA setup for floating-point (FP) representation. The
binary activations can replace the FP multiplications with logical operations,
i.e., conditional assignment to 0 with a bank of AND gates.

TABLE V: Comparison of the normalized compute (CE), memory (ME), and total (TE) energy between our dynamic time step SNN (DTS SNN), SOTA spikformer-based SNN (Spik. SNN) and non-spiking ViT (NS ViT). The number of time steps (TS) for the dynamic time step and spikformer-based SNN are chosen at iso-accuracy. Dataset | Type | TS |  Norm. CE |  Norm. ME |  Norm. TE
---|---|---|---|---|---
CIFAR10 | NS ViT | - | 22.57$\times$ | 0.84$\times$ | 16.12$\times$
CIFAR10 | Spik. SNN | 4 | 2.79$\times$ | 3.68$\times$ | 2.98$\times$
CIFAR10 | DTS SNN | 1.15 | 1$\times$ | 1$\times$ | 1$\times$
CIFAR100 | NS ViT | - | 25.7$\times$ | 0.88$\times$ | 17.15$\times$
CIFAR100 | Spik. SNN | 6 | 1.70$\times$ | 2.32$\times$ | 1.89$\times$
CIFAR100 | DTS SNN | 3.81 | 1$\times$ | 1$\times$ | 1$\times$

Since our dynamic time step SNNs skip time steps, and consequently, skip more
FLOPs compared to existing spikformer-based SNNs (and even more FLOPs compared
to non-spiking ViTs that do not have sparsity), as shown in Table 1, our
proposed approach leads to higher compute efficiency. These lower FLOPs,
coupled with the $7.4\times$ reduction for AC operations leads to a
$22.57\times$ and $25.7\times$ reduction in compute energy on CIFAR10 and
CIFAR100 respectively compared to non-spiking ViTs with Spikformer-4-384.

Additionally, our SNNs also enjoy superior memory-efficiency compared to
existing SNNs since the latter requires less number of membrane potential
reads and writes to and from the on-chip memory. We can avoid the potential
access for the skipped time step using zero gating logic, that can be
supported in neuromorphic chips, such as Loihi. This potential access
dominates the memory footprint of the SNNs during inference [9, 31] since it
is not influenced by the SNN sparsity since each potential is the sum of a
large number of (typically a few 100s) weight-modulated spikes, and so it is
almost impossible that all these spike values are zero for the membrane
potential to be kept unchanged at a particular time step. Note that the number
of weight read and write accesses can be reduced with the spike sparsity, and
thus typically do not dominate the memory footprint of the SNN. However, the
exact memory savings would depend on the data re-use and underlying hardware.
Assuming a simple weight stationary scheme [32], our SNNs (1.15 time steps)
yield $3.68\times$ and $2.32\times$ reduced memory energy compared to
Spikformer (4 time steps) at iso-accuracy on CIFAR10 and CIFAR100
respectively. The total (compute+memory) energy savings become $16.12\times$
and $2.98\times$ compared to non-spiking ViT and Spikingformer on CIFAR10. The
detailed energy models used to evaluate these numbers are in [33, 34], while
the energy per operation numbers are obtained from our FPGA simulations.

##  V Discussions

Since the introduction of deep SNNs, many researchers focused on reducing the
total number of time steps while not sacrificing the accuracy compared to the
non-spiking counterparts. Existing research on low-time-step SNNs have
assigned the same number of time steps to all the layers and obtained close to
SOTA accuracy. However, we empirically show that this approach is only limited
to CNNs (we are the first to the best of our knowledge), and fails when
applied to ViTs. We conduct a comprehensive study of the sensitivity of the
layers in the model to the number of time steps in order to guide our decision
to select the number of time steps for each layer. Based on this study, we
propose a training framework with dynamic time step spiking layers as well as
a mask loss. The spiking ViT models obtained by our proposed framework
outperform the SOTA SNN models with $20\%\sim 30\%$ lower latency with no drop
in test accuracy on CIFAR10/100 datasets.

Our dynamic time step approach is related to the recently proposed work, SEENN
[16], that also adjusts the number of time steps dyanmically. However, SEENN
conditions the time steps based on the difficulty of the input sample, and not
sensitivity of the SNN block, and is limited to CNNs111The implementation of
SEENN is not made public yet for us to evaluate it on ViTs. We provide
insights on how the sensitivity of the activation distributions in ViT can be
leveraged to adjust the time steps independent of the sample difficulty. SEENN
can potentially be applied to our work to further improve the accuracy-
timestep trade-off.

## References

  * [1] M. Pfeiffer _et al._ , “Deep learning with spiking neurons: Opportunities and challenges,” _Frontiers in Neuroscience_ , vol. 12, p. 774, 2018.
  * [2] Y. Cao _et al._ , “Spiking deep convolutional neural networks for energy-efficient object recognition,” _International Journal of Computer Vision_ , vol. 113, pp. 54–66, 05 2015.
  * [3] A. Sengupta _et al._ , “Going deeper in spiking neural networks: VGG and residual architectures,” _Frontiers in Neuroscience_ , vol. 13, p. 95, 2019.
  * [4] S. Deng and S. Gu, “Optimal conversion of conventional artificial neural networks to spiking neural networks,” in _International Conference on Learning Representations_ , 2021. [Online]. Available: https://openreview.net/forum?id=FZ1oTwcXchK
  * [5] G. Datta and P. A. Beerel, “Can deep neural networks be converted to ultra low-latency spiking neural networks?” in _2022 Design, Automation & Test in Europe Conference & Exhibition (DATE)_, vol. 1, no. 1, 2022, pp. 718–723.
  * [6] Y. Wu _et al._ , “Direct training for spiking neural networks: Faster, larger, better,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, pp. 1311–1318, Jul. 2019.
  * [7] P. Panda _et al._ , “Toward scalable, efficient, and accurate deep spiking neural networks with backward residual connections, stochastic softmax, and hybridization,” _Frontiers in Neuroscience_ , vol. 14, 2020.
  * [8] G. Datta, S. Kundu, and P. A. Beerel, “Training energy-efficient deep spiking neural networks with single-spike hybrid input encoding,” in _2021 International Joint Conference on Neural Networks (IJCNN)_ , vol. 1, no. 1, 2021, pp. 1–8.
  * [9] G. Datta _et al._ , “Hoyer regularizer is all you need for ultra low-latency spiking neural networks,” _arXiv preprint arXiv:2212.10170_ , 2022.
  * [10] N. Rathi and K. Roy, “Diet-snn: A low-latency spiking neural network with direct input encoding and leakage and threshold optimization,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 1, no. 1, pp. 1–9, 2021.
  * [11] N. Rathi, G. Srinivasan, P. Panda, and K. Roy, “Enabling deep spiking neural networks with hybrid conversion and spike timing dependent backpropagation,” in _International Conference on Learning Representations_ , 2020. [Online]. Available: https://openreview.net/forum?id=B1xSperKvH
  * [12] G. Datta, S. Kundu, A. R. Jaiswal, and P. A. Beerel, “ACE-SNN: Algorithm-Hardware co-design of energy-efficient & low-latency deep spiking neural networks for 3D image recognition,” _Frontiers in Neuroscience_ , vol. 16, 2022. [Online]. Available: https://www.frontiersin.org/articles/10.3389/fnins.2022.815258
  * [13] G. Datta, H. Deng, R. Aviles, Z. Liu, and P. A. Beerel, “Bridging the gap between spiking neural networks & lstms for latency & energy efficiency,” in _2023 IEEE/ACM International Symposium on Low Power Electronics and Design (ISLPED)_ , vol. 1, no. 1, 2023, pp. 1–6.
  * [14] S. Deng _et al._ , “Temporal efficient training of spiking neural network via gradient re-weighting,” in _International Conference on Learning Representations_ , 2022. [Online]. Available: https://openreview.net/forum?id=_XNtisL32jv
  * [15] S. S. Chowdhury _et al._ , “Towards ultra low latency spiking neural networks for vision and sequential tasks using temporal pruning,” in _Computer Vision – ECCV 2022_. Springer Nature Switzerland, 2022, pp. 709–726.
  * [16] Y. Li _et al._ , “Seenn: Towards temporal spiking early-exit neural networks,” _arXiv preprint arXiv:2304.01230_ , 2023.
  * [17] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [18] K. He _et al._ , “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [19] A. Dosovitskiy _et al._ , “An image is worth 16x16 words: Transformers for image recognition at scale,” in _International Conference on Learning Representations_ , 2021. [Online]. Available: https://openreview.net/forum?id=YicbFdNTTy
  * [20] Z. Zhou _et al._ , “Spikformer: When spiking neural network meets transformer,” _arXiv preprint arXiv:2209.15425_ , 2022.
  * [21] Y. Li _et al._ , “Spikeformer: A novel architecture for training high-performance low-latency spiking neural network,” _arXiv preprint arXiv:2211.10686_ , 2022.
  * [22] C. Lee _et al._ , “Enabling spike-based backpropagation for training deep neural network architectures,” _Frontiers in Neuroscience_ , vol. 14, 2020.
  * [23] E. O. Neftci _et al._ , “Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks,” _IEEE Signal Processing Magazine_ , vol. 36, no. 6, pp. 51–63, 2019.
  * [24] A. Hassani _et al._ , “Escaping the big data paradigm with compact transformers,” _arXiv preprint arXiv:2104.05704_ , 2021.
  * [25] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * [26] A. Krizhevsky _et al._ , “Imagenet classification with deep convolutional neural networks,” in _Advances in Neural Information Processing Systems 25_ , 2012, pp. 1097–1105.
  * [27] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” _arXiv preprint arXiv:1711.05101_ , 2017.
  * [28] H. Zheng _et al._ , “Going deeper with directly-trained larger spiking neural networks,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 12, pp. 11 062–11 070, May 2021.
  * [29] Q. Meng _et al._ , “Training high-performance low-latency spiking neural networks by differentiation on spike representation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 12 444–12 453.
  * [30] M. Horowitz, “Computing’s energy problem (and what we can do about it),” in _2014 IEEE International Solid-State Circuits Conference Digest of Technical Papers (ISSCC)_ , 2014, pp. 10–14.
  * [31] R. Yin, A. Moitra, A. Bhattacharjee, Y. Kim, and P. Panda, “Sata: Sparsity-aware training accelerator for spiking neural networks,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 1, no. 1, pp. 1–1, 2022.
  * [32] V. Sze, Y.-H. Chen, T.-J. Yang, and J. Emer, “Efficient processing of deep neural networks: A tutorial and survey,” _arXiv preprint arXiv:1703.09039_ , 2017.
  * [33] S. Kundu, G. Datta, M. Pedram, and P. A. Beerel, “Spike-thrift: Towards energy-efficient deep spiking neural networks by limiting spiking activity via attention-guided compression,” in _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV)_ , January 2021, pp. 3953–3962.
  * [34] F. Ottati, C. Gao, Q. Chen, G. Brignone, M. R. Casu, J. K. Eshraghian, and L. Lavagno, “To spike or not to spike: A digital hardware perspective on deep learning acceleration,” _arXiv preprint arXiv:2306.15749_ , 2023.


11institutetext: Department of Physics, Swansea University, Swansea SA2 8PP,
United Kingdom 22institutetext: European Centre for Theoretical Studies in
Nuclear Physics and Related Areas (ECT*) & Fondazione Bruno Kessler Strada
delle Tabarelle 286, 38123 Villazzano (TN), Italy 33institutetext: Department
of Theoretical Physics, National University of Ireland Maynooth, County
Kildare, Ireland 44institutetext: CP3-Origins & Danish IAS, Department of
Mathematics and Computer Science, University of Southern Denmark, 5230 Odense
M, Denmark 55institutetext: Department of Physics, Sejong University, Seoul
143-747, Korea 66institutetext: INFN, Sezione di Firenze, 50019 Sesto
Fiorentino (FI), Italy 77institutetext: School of Mathematics, Trinity
College, Dublin, Ireland

# A comparison of spectral reconstruction methods applied to non-zero
temperature NRQCD meson correlation functions111Talk presented at A Virtual
Tribute to Quark Confinement and the Hadron Spectrum, August 2-6 2021,
University of Stavanger, Stavanger, Norway.

Thomas Spriggs 11<EMAIL_ADDRESS>Gert Aarts 1122 Chris Allton 11 Timothy
Burns 11 Rachel Horohan D’Arcy 33 Benjamin Jäger 44 Seyong Kim 55 Maria-Paola
Lombardo 66 Sam Offler 11 Ben Page 11 Sinead M. Ryan 77 Jon-Ivar Skullerud 33

###### Abstract

We present results from the fastsum collaboration’s programme to determine the
spectrum of the bottomonium system as a function of temperature. Three
different methods of extracting spectral information are discussed: a Maximum
Likelihood approach using a Gaussian spectral function for the ground state,
the Backus Gilbert method, and the Kernel Ridge Regression machine learning
procedure. We employ the fastsum anisotropic lattices with 2+1 dynamical quark
flavours, with temperatures ranging from 47 to 375 MeV.

##  1 Introduction

There has been a great deal of interest in onia systems in the context of
heavy-ion collision experiments, particularly since the proposal they behave
as a proxy for the temperature Matsui:1986dk . While initial work concentrated
on the charmonium system, interest turned to bottomonium for several reasons:
these mesons are produced copiously in LHC heavy-ion collision experiments,
they act as probes of the quark-gluon plasma, and results from the CMS
Collaboration indicate sequential suppression in this system CMS:2011all .

The fastsum collaboration has had a long programme of studying the bottomonium
spectrum at non-zero temperature using the NRQCD method Aarts:2014cda ,
principally using the Maximum Entropy Method Asakawa:2000tr . We have
determined both S- and P-wave masses and determined upper bounds for the
state’s widths. We extend this work here to include three new analysis
techniques to determine the bottomonium spectrum, including the widths of the
states. These approaches are: a Maximum Likelihood approach using a Gaussian
spectral function for the ground state (see sec. 3), the Backus Gilbert method
(sec. 4), and the Kernel Ridge Regression machine learning procedure (sec. 5).
This contribution contains a progress report on this project, in which the
long-term aim is to test and compare these and other methods, using the same
lattice dataset in order to extract the best estimates of the bottomonium
spectrum at finite temperature.

##  2 Lattice Method

Because the $b$-quark’s mass is larger than any other mass scale, it can be
approximated as a non-relativistic particle and the NRQCD effective theory can
be used for its dynamics. Note that this effective theory also holds in the
case of thermal QCD since the temperature reached in heavy-ion collisions is
less than the $b$-quark mass. We use an ${\cal O}(v^{4})$ lattice
implementation of NRQCD where $v=|{\bf p}|/M$ is the velocity of the $b$-quark
in the bottomonium’s rest frame. Details of the action and the tuning of the
$b$-quark parameters can be found in Aarts:2014cda .

NRQCD has simpler time evolution properties compared to relativistic theories,
it is an “initial value” problem with the quark propagating forward in time
only. The missing backward movers lead to a simpler spectral decomposition.
The meson correlator, $G(\tau)$, can be expressed in terms of the spectral
function, $\rho(\omega)$, via

$G(\tau)=\int^{\omega_{max}}_{\omega_{min}}\frac{\mathrm{d}\omega}{2\pi}K(\tau,\omega)\rho(\omega),\qquad\qquad
K(\tau,\omega)=e^{-\omega\tau}.$ (1)

In NRQCD the kernel $K(\tau,\omega)$ takes the simple exponential form,
whereas in the relativistic case it is
$\cosh[\omega(\tau-1/2T)]/\sinh(\omega/2T)$.

The challenge is to reconstruct $\rho(\omega)$ from $G(\tau)$. Eq. (1)
illustrates the “ill-posed” nature of this problem: $G(\tau)$ is typically
known at ${\cal O}(10-100)$ data points, whereas $\rho(\omega)$ is a
continuous function which requires ${\cal O}(1000)$ points to correctly
represent it. This is a well-known problem and many techniques have been
developed across several research fields to solve this problem. In the lattice
context, it is important to compare and test these methods to ascertain which
one(s) are most applicable.

fastsum’s “Generation 2L” anisotropic lattice ensembles are used, with a
physical size of (32$a_{s})^{3}\times(N_{\tau}a_{\tau})$ where
$a_{\tau}^{-1}=5.997(34)$ GeV and $a_{s}/a_{\tau}=3.453(6)$ Edwards:2008ja ;
HadronSpectrum:2008xlg ; Cheung:2016bym . The simulation is performed with 2+1
dynamical Wilson-clover quark flavours where the pion mass is $M_{\pi}=236(2)$
MeV. There are ${\cal O}(1000)$ configurations for all temperatures and the
pseudocritical temperature is $162(1)$ MeV, as measured from the inflection
point of the renormalised chiral condensate Aarts:2020vyb . The temperatures
and corresponding $N_{\tau}$ values are shows in Table 1.

$N_{\tau}$ | 16 | 20 | 24 | 28 | 32 | 36 | 40 | 48 | 56 | 64 | 128
---|---|---|---|---|---|---|---|---|---|---|---
T [MeV] | 375 | 300 | 250 | 214 | 187 | 167 | 150 | 125 | 107 | 94 | 47
Table 1: Temporal lattice sizes and corresponding temperatures for the fastsum
Generation 2L ensembles Aarts:2020vyb . The $N_{\tau}=128$ ensemble was kindly
provided by the Hadron Spectrum Collaboration Edwards:2008ja ;
HadronSpectrum:2008xlg ; Cheung:2016bym .

To start the discussion we show in Fig.1 the effective mass in lattice units
for the $\Upsilon$ meson, defined via

$M_{\text{eff}}(\tau)=\log\left(\frac{G(\tau)}{G(\tau+1)}\right),$ (2)

for all the temperatures considered. Some evidence of thermal effects can be
observed.

Figure 1: Effective mass plot in the $\Upsilon$ channel, for all temperatures
considered.

##  3 Gaussian Maximum Likelihood

In this approach the bottomonium spectral function is parametrised by a
Gaussian ground state (see also Larsen:2019bwy ) and a single
$\delta$-function to account for all spectral weight beyond the ground state,

$\rho_{\text{ansatz}}(\omega)=A_{\text{ground}}\exp\left(\frac{(\omega-
M_{\text{ground}})^{2}}{2\sigma^{2}}\right)+A_{\text{excited}}\;\delta(\omega-
M_{\text{excited}}).$ (3)

The Gaussian ansatz for the ground state is chosen to allow for a finite
width, which is expected in a thermal medium. The spectral weight beyond the
ground state could be modelled with a function more sophisticated than a
$\delta$-function, but the extra parameters in any such function would reduce
the method’s predictability. Inserting Eq. (3) into Eq. (1) leads to a closed
form expression for the correlation function $G_{\text{ansatz}}(\tau)$. We
then apply the standard Maximum Likelihood method to determine the best fit
parameters $\\{M_{\text{ground/excited}},A_{\text{ground/excited}},\sigma\\}$.

In the Maximum Likelihood method, we need to chose a time window
$\tau\in[\tau_{1},\tau_{2}]$ to perform the fits. This introduces a systematic
effect which we study in Fig. 2 (left). Here the full width at half maximum
(FWHM) is plotted against $1/\tau_{2}$ for the $\Upsilon$ channel at $T=47$
MeV with $\tau_{1}=8$ throughout. Note that we use temporal lattice units
here, $a_{\tau}\equiv 1$, such that the $\tau$’s are integers. We expect the
best estimate of the width to be given as $1/\tau_{2}\rightarrow 0$, since
this isolates the ground state. However, at non-zero temperature the temporal
extent of the lattice is naturally restricted. As can be seen, the width
heavily depends on the time window. We consider two extrapolations
$1/\tau_{2}\rightarrow 0$; the blue line is a linear fit to the eight leftmost
(blue) data points, whereas the red line is a linear fit to all (blue and red)
points. Taking these extrapolations at face-value yields a width estimate of
$\sim 10$ MeV.

The effect of varying the time window $\tau\in[8,\tau_{2}]$ is shown in Fig. 2
(right) for all temperatures. Note that results from the same time window are
plotted with the same colour. The variation in the width as the fit parameter
$\tau_{2}$ is varied is substantial. Further details can be found elsewhere
tom_lat21 , including a more detailed analysis of the systematic effects. As a
preliminary result we note from the right lower pane the indication that the
width increases above $T_{\rm pc}$ when results from the same time window but
at different temperatures are compared.

Figure 2: Left: Full width at half maximum (FWHM) for the $\Upsilon$ at $T=47$
MeV as a function of the inverse time window parameter $1/\tau_{2}$, with
$\tau_{1}=8$ fixed. The red (blue) line is a linear extrapolation on all
(eight leftmost) data points. Right: FWHM of the $\Upsilon$ state for all
temperatures for various $\tau_{2}$ values. The pseudo-critical temperature,
$T_{\rm pc}=162(1)$ MeV, is shown by vertical dashed line. The lower pane
contains a close-up using $\tau_{1}=8$ and $\tau_{2}=19$ for all temperatures.

##  4 Backus Gilbert

The Backus-Gilbert method was introduced originally to solve an ill-posed
problem relevant to geology Backus1968 . It estimates a solution to Eq. (1),
$\hat{\rho}(\omega)$, in a point-wise manner by sampling the target spectrum,
$\rho(\omega)$, using a basis of resolution functions, $A(\omega,\omega_{0})$,
that are peaked around some $\omega_{0}\in[\omega_{min},\omega_{max}]$,

$\hat{\rho}(\omega_{0})=\int_{\omega_{min}}^{\omega_{max}}A(\omega,\omega_{0})\rho(\omega)d\omega.$
(4)

Ideally, these resolution functions are as close to the delta function
$\delta(\omega-\omega_{0})$ as possible. They are expressed as a linear
combination of the kernel function,

$A(\omega,\omega_{0})=\sum_{\tau}c_{\tau}(\omega_{0})K(\omega,\tau),$ (5)

where the coefficients, $c_{\tau}(\omega_{0})$ are to be determined. This
means that, by combining Eqs. (4) and (5), the solution estimate
$\hat{\rho}(\omega_{0})$ is obtained from a linear combination of the original
correlation functions, $G(\tau)$,

$\hat{\rho}(\omega_{0})=\sum_{\tau}c_{\tau}(\omega_{0})G(\tau).$ (6)

The coefficients $c_{\tau}(\omega_{0})$ can be determined by the Dirichlet
least-squares criterion of minimising the distance between the resolution
functions and the delta function,

$J(\omega_{0})=\int_{\omega_{min}}^{\omega_{max}}[A(\omega,\omega_{0})-\delta(\omega-\omega_{0})]^{2}d\omega,$
(7)

leading to

$\mathcal{K}_{\tau\tau^{\prime}}\cdot
c_{\tau^{\prime}}(\omega_{0})=K(\omega_{0},\tau)\hskip
28.45274pt\text{where}\hskip
5.69054pt\mathcal{K}_{\tau\tau^{\prime}}=\int_{\omega_{min}}^{\omega_{max}}K(\omega,\tau)K(\omega,\tau^{\prime})d\omega.$
(8)

The variance in the solution is also a function of $c_{\tau}(\omega_{0})$ and
is given by

$\text{Var}[\hat{\rho}(\omega_{0})]=\sum_{\tau,\tau^{\prime}}c_{\tau}\Sigma_{\tau,\tau^{\prime}}c_{\tau^{\prime}},$
(9)

where $\Sigma_{\tau,\tau^{\prime}}$ is the covariance in $G(\tau)$.

Figure 3: Left: Reconstructed spectral functions of $\Upsilon$ from the
Backus-Gilbert method for a sample of temperatures. Right: The ground state’s
FWHM against temperature with the vertical dashed line indicating $T_{\rm
pc}$. The grey band shows the minimum resolution, and its error band, that the
resolution functions could achieve at a given temperature.

Reconstructed spectral functions for a sample of temperatures are shown in
Fig. 3 (left). From this we can see little change in the spectrum between 47
and 97 MeV, but then an increase in mass and a large apparent increase in
width at larger temperatures. We note that the features below 9 GeV are
artefacts of the method. In Fig. 3 (right) the FWHM is plotted against
temperature where the error bars are statistical only. At first sight, there
appears a clear increase in the FWHM with temperature. However, also plotted
in the grey band is the FWHM of the resolution function $A(\omega,\omega_{0})$
and its error bound calculated at $\omega_{0}=M_{\text{ground}}$. This
represents the lower bound in the width that the method can resolve – the
method is incapable of predicting widths below this bound, and only widths
which lie above this bound can be taken as predictions. Since the bottomonium
widths calculated using this method follow this lower bound, the Backus
Gilbert method is currently incapable of resolving the bottomonium width.
Further discussion can be found in ben_lat21 ; inpreparation .

##  5 Kernel Ridge Regression

As is common in machine learning paradigms, kernel ridge regression (KRR)
infers a prediction of a quantity based on previously observed training data.
In the case of solving the inversion problem for NRQCD, mock data with known
spectral functions are generated as training data.

Details of the KRR approach in the context of NRQCD can be found in sam-old .
Here we note that we modified our method how the input data, the set of
correlators $G_{i}(\tau)$ (where $i=1,2,\ldots,N_{\rm train})$, is combined
into the matrix (or kernel) $\mathbf{C}$ before the regression procedure. In
contrast to sam-old , we use here the following encoding for the matrix
elements of $\mathbf{C}$,

$C_{ij}=\exp\left(-\gamma\sum_{n}\left[\frac{G_{i}(\tau_{n})-G_{j}(\tau_{n})}{\overline{G}(\tau_{n})}\right]^{2}\right),$
(10)

where the sum $n$ runs over all time slices (starting at $n=4$), $\gamma$ is a
hyperparameter that sets the correlation length in the space of correlators,
and $\overline{G}(\tau_{n})=1/N_{\rm train}\sum_{i}G_{i}(\tau_{n})$ is the
average correlator of the training data at time $\tau_{n}$.

Figure 4: Left: Frequency of central values generated to form the set of mock
correlation functions overlaid with the actual correlation functions, the data
is constrained to lie within the dashed lines and be most frequent near the
true data. Right: Reconstructed spectral functions at a sample of temperatures
with a zoom in of the ground state peak highlighted to observe the temperature
dependence.

We continue with some remarks on the generation of mock data, which has
applications beyond KRR. KRR works by inferring a spectral function from a
training set of mock spectral functions weighted by how much the concomitant
correlation function resembles the true lattice correlation function. Clearly
then, the closer the training correlation functions are to the actual lattice
correlation functions obtained via numerical simulations, the more relevant
training spectral functions the KRR has to draw upon to infer a spectral
function. Building on sam-old , work has been done to increase the overlap of
the training data and the lattice data. The overlap between the new training
data and the actual lattice correlators is shown in Fig. 4 (left). Here, the
solid lines are the lattice correlators obtained in the simulations at a given
temperature ($N_{\tau}=128$ in this case). The heatmap shows the distribution
of the training data around these values and the dashed lines enclose the
region inside which all training data lies. As can be seen, the mock data
tracks the actual correlators well. More details can be found in sam_lat21 .

Resulting spectral functions are shown in Fig. 4 (right). We observe ground
state peaks at all four temperatures shown, with some temperature dependence
visible in the left pane. It is noted that the mass of both the ground state
and first excited state lie below the experimental values. Further discussion
is presented in sam_lat21 ; inpreparation .

##  6 Conclusion

This work presents a progress report of the fastsum collaboration’s programme
to apply a number of methods to reconstruct the bottomonium spectral functions
from NRQCD lattice correlation functions. The three methods used are a Maximum
Likelihood approach using a Gaussian spectral function for the ground state;
the Backus Gilbert method; and the Kernel Ridge Regression machine learning
procedure. The results presented are preliminary and indicate the substantial
systematic uncertainties.

Further work to analyse these systematic effects is in progress inpreparation
.

## Acknowledgements

This work is supported by STFC grant ST/T000813/1. RHD has been supported by a
Maynooth University SPUR scholarship. SK is supported by the National Research
Foundation of Korea under grant NRF-2021R1A2C1092701 funded by the Korean
government (MEST). BP has been supported by a Swansea University Research
Excellence Scholarship (SURES). This work used the DiRAC Extreme Scaling
service at the University of Edinburgh, operated by the Edinburgh Parallel
Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).
This equipment was funded by BEIS capital funding via STFC capital grant
ST/R00238X/1 and STFC DiRAC Operations grant ST/R001006/1. DiRAC is part of
the National e-Infrastructure. This work was performed using PRACE resources
at Cineca via grants 2015133079 and 2018194714. We acknowledge the support of
the Supercomputing Wales project, which is part-funded by the European
Regional Development Fund (ERDF) via Welsh Government, and the University of
Southern Denmark for use of computing facilities. We are grateful to the
Hadron Spectrum Collaboration for the use of their zero temperature ensemble.

## References

  * (1) T. Matsui and H. Satz, Phys. Lett. B 178 (1986), 416-422.
  * (2) S. Chatrchyan et al. [CMS], Phys. Rev. Lett. 107 (2011), 052302 [arXiv:1105.4894 [nucl-ex]].
  * (3) G. Aarts, C. Allton, T. Harris, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, JHEP 07 (2014), 097 [arXiv:1402.6210 [hep-lat]].
  * (4) M. Asakawa, T. Hatsuda and Y. Nakahara, Prog. Part. Nucl. Phys. 46 (2001), 459-508 [arXiv:hep-lat/0011040 [hep-lat]].
  * (5) R. G. Edwards, B. Joo and H. W. Lin, Phys. Rev. D 78 (2008), 054501 doi:10.1103/PhysRevD.78.054501 [arXiv:0803.3960 [hep-lat]].
  * (6) H. W. Lin et al. [Hadron Spectrum], Phys. Rev. D 79 (2009), 034502 doi:10.1103/PhysRevD.79.034502 [arXiv:0810.3588 [hep-lat]].
  * (7) G. K. C. Cheung et al. [Hadron Spectrum], JHEP 12 (2016), 089 doi:10.1007/JHEP12(2016)089 [arXiv:1610.01073 [hep-lat]].
  * (8) G. Aarts, C. Allton, J. Glesaaen, S. Hands, B. Jäger, S. Kim, M. P. Lombardo, A. A. Nikolaev, S. M. Ryan and J. I. Skullerud, et al. [arXiv:2007.04188 [hep-lat]].
  * (9) R. Larsen, S. Meinel, S. Mukherjee and P. Petreczky, Phys. Rev. D 100 (2019) no.7, 074506 [arXiv:1908.08437 [hep-lat]].
  * (10) T. Spriggs et al., PoS LATTICE2021 (2021), 077 [arXiv:2112.01599 [hep-lat]].
  * (11) Fastsum, in preparation.
  * (12) G. Backus and F. Gilbert. Geophysical Journal International, 16(2) (1968) 169-205.
  * (13) B. Page et al., PoS LATTICE2021 (2021), 134 [arXiv:2112.02075 [hep-lat]].
  * (14) D. W. Oldenburg, "An introduction to linear inverse theory," IEEE Transactions on Geoscience and Remote Sensing, GE-22 (1984) 665-674.
  * (15) S. Offler, G. Aarts, C. Allton, J. Glesaaen, B. Jäger, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, PoS LATTICE2019 (2019), 076 [arXiv:1912.12900 [hep-lat]].
  * (16) S. Offler et al., PoS LATTICE2021 (2021), 509 [arXiv:2112.02116 [hep-lat]].