\section{Introduction}
\label{sec:introduction}
\if0
\begin{figure}[tb]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=\columnwidth]{experiment/comparison-scaling-factor-cnn6-flops.pdf}}
\caption{Comparison among layer-wise penalty terms. The $x$-axis represents the layer number,
and the $y$-axis represents the ratio of some layer-wise metric or penalty term to some total metric or penalty term.
The network architecture is CNN7 (App.~\ref{app:detail-of-network-architecture}), and the computation details are described in  Sec.~\ref{sec:method}.
Our penalty term (blue) is precisely proportional to the ground truth (gray).}
\label{fig:comparison-scaling-factor-cnn6-flops}
\end{center}
\vskip -0.2in
\end{figure}
\fi
With the rapid growth and spread of neural networks, realizing energy-efficient neural networks is an urgent mission for sustainable development.
One such model is the spiking neural network (SNN), which is also known to be more biologically plausible than ordinary artificial neural networks (ANNs).
SNNs are energy-efficiently driven on neuromorphic chips~\citep{akopyan2015truenorth,davies2018loihi} or certain field-programmable gate arrays (FPGAs)~\citep{maguire2007challenges,misra2010artificial} by asynchronously processing spike signals.
However, as the spike firing rate of an SNN increases, the energy consumption does as well, and thus, the advantage of the SNN diminishes.
Therefore, in addition to the shift from ANNs to SNNs,
it is advantageous to adopt training methods that reduce energy consumption in the inference phase.
At the same time, such a training method should be independent of the network architecture to avoid limitations in the application.
That is, our goal is to develop a training method that realizes energy-efficient SNNs without any constraint on the network architecture.

There are various approaches toward energy-efficient SNNs, such as pruning, quantization, and knowledge distillation~\citep{kundu2021spike,chowdhury2021spatio,lee2021energy}, which are widely-used approaches also in ANNs.
Further, there are SNN-specific approaches sparsifying the spiking activity related to the energy consumption~\citep{lee2020enabling,kim2021revisiting,naya2021spiking}, to which our method belongs.
In particular, the methods that penalize the spike firing rate in the training phase are close to our aforementioned goal~\citep{esser2016cover,sorbaro2020optimizing,pellegrini2021low}.
However, they indirectly reduce the energy consumption by arbitrarily reducing the spike firing rate, where there is no strict positive correlation between them.
Hence, reducing energy consumption while maintaining accuracy is difficult.


\begin{figure}[tb]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=\columnwidth]{experiment/comparison-scaling-factor-cnn6-flops.pdf}}
\caption{Comparison among layer-wise penalty terms. The $x$-axis represents the layer number,
and the $y$-axis represents the ratio of some layer-wise metric or penalty term to some total metric or penalty term.
The network architecture is CNN7 (App.~\ref{app:detail-of-network-architecture}), and the computation details are described in  Sec.~\ref{sec:method}.
Our penalty term (blue) is precisely proportional to the ground truth (gray).}
\label{fig:comparison-scaling-factor-cnn6-flops}
\end{center}
\vskip -0.2in
\end{figure}

\paragraph{Principle and Idea}
Our principle is that---we should optimize the metric as it is in the training phase.
In this spirit, we propose to introduce a proper penalty term for the spiking activity---a \textit{spiking synaptic penalty}---into the objective function.
It is derived so that its expected value is precisely proportional to the energy consumption metric for the SNN (see Fig.~\ref{fig:comparison-scaling-factor-cnn6-flops}).
Although the difference between the proposed and existing methods is only at the scaling factor for each spike, we demonstrate that this minor correction causes significant improvement.


\paragraph{Main Contributions}
\begin{itemize}
    \item We derived a novel penalty term that can directly optimize the metric for the total energy consumption of an SNN without modifying the network architecture.
    \item We demonstrated that the proposed method is compatible with the weight decay, which imposes implicit sparsity on the network~\citep{yaguchi2018adam}, and that the proposed method creates a higher sparsification effect than the weight decay.
    \item We demonstrated that the proposed method can reduce the energy consumption more than other methods while maintaining the accuracy for image classification tasks, which mitigates the dilemma of the energy--accuracy trade-off.
\end{itemize}


\section{Related Work}
\label{sec:related-work}

\subsection{Spike Sparsification in Direct SNN Training}
\label{sec:related-work_direct-training}
The most relevant approaches to our proposal introduce the penalty term for the spike firing rate and directly train SNNs by the surrogate gradient method~\citep{esser2016cover, pellegrini2021low}.
It is a straightforward idea to penalize the spike firing rate to obtain energy-efficient SNNs because it appears in the SNN energy consumption metric~\citep{lee2020enabling,kim2021revisiting}.
We refer to the reduction in spike firing rate as spike sparsification.
Although the spike firing rate cannot be optimized by the ordinal backpropagation method owing to non-differentiability, it can be optimized by the surrogate gradient method, which is the same technique as training spiking neurons of SNNs~\citep{zenke2018superspike,shrestha2018slayer}.
However, neither of these penalty terms (\citet{esser2016cover,pellegrini2021low}) precisely matches the energy consumption metric.
As opposed to them, our spiking synaptic penalty resolves this limitation.

\subsection{Spike Sparsification via Conversion from ANN}
\label{sec:related-work_conversion}
Other approaches introduce the penalty term for corresponding ReLU networks (ANNs with ReLU activations) and convert them to SNNs~\citep{sorbaro2020optimizing,narduzzi2022optimizing}.
Although there is no guarantee that the penalty terms for ReLU networks contribute to the reduction of the energy consumption for converted SNNs, ReLU networks can be optimized by the ordinal backpropagation method.
Note that the same synaptic scaling factor for the penalty term as ours is proposed to reduce the energy consumption for SNNs in \citet{sorbaro2020optimizing}.
However, they failed to provide evidence to support their claim, as they mentioned.
As opposed to them, we provide theoretical and experimental proof in the setting of the direct SNN training by the surrogate gradient.


\subsection{Neuron Sparsification}
\label{sec:related-work_others}

Neuron sparsification means increasing the number of permanently zero-valued activations---dead neurons. It is a stronger condition than spike sparsification, which does not force neurons to be permanently inactive.
In ReLU networks, the training with the Adam optimizer and weight decay regularization implicitly induce neuron sparsification~\citep{yaguchi2018adam} because ReLU activations have an inactive state.
However, this claim has yet to be demonstrated in the context of SNNs, where spiking neurons also have an inactive state, even though weight decay is usually adopted in SNNs.
Therefore, to detect the effect of our method correctly, we shall also focus on the weight decay.


\section{Method}
\label{sec:method}

In this section, we propose the spiking synaptic penalty.
First, we describe the spiking neuron model with surrogate gradient mechanism.
Next, we describe the metric for energy consumption, which can be represented by the spiking activity.
Note that we need to optimize both the accuracy and energy efficiency in the training phase.
Finally, we state that the spiking synaptic penalty is the proper penalty term to optimize the energy consumption metric.

\subsection{Neuron Model and Surrogate Gradient}
\label{sec:neuron-model}
In this study, we uses SNNs constructed by single-step spiking neurons, which are superior to the multi-time step SNNs in terms of training and inference costs for static tasks \citep{suetake2023s3nn}.
Note that the single-step spiking neurons are the same setup as that in a previous study of the penalty term~\citep{esser2016cover}.

Let us denote $l \in \{l \in \boldsymbol{Z} \mid 1 \leq l \leq L\}$ as the layer, $d_l$ as the number of neurons in the $l$-th layer, $\boldsymbol{s}_{0} = x \in X \subset \boldsymbol{R}^{d_{0}}$ as the input data, and the subscript $i$ of any vector as its $i$-th component.
Then, the single-step spiking neuron is defined as follows~\citep{suetake2023s3nn}:
\begin{definition}
The forward mode of a single-step spiking neuron consists of two ingredients, the membrane potential $\boldsymbol{u}_{l} \in \boldsymbol{R}^{d_{l}}$ and spikes emitted by neurons $\boldsymbol{s}_{l} \in \{0, 1\}^{d_{l}}$. They are defined using the Heaviside step function $H$ as follows:
\begin{align}
\label{eqn:membrane-potential}
\boldsymbol{u}_{l}
&:= \boldsymbol{W}_{l}\boldsymbol{s}_{l-1}, \\
{}
\label{eqn:spike}
s_{l, i}\left(u_{l, i}\right)
&:= H\left(u_{l, i}-u_{\rm th}\right)
= \left\{
\begin{array}{ll}
1 & \left(u_{l, i}\ge u_{\rm th}\right), \\
0 & \left(u_{l, i}< u_{\rm th}\right),
\end{array}
\right.
\end{align}
where $\boldsymbol{W}_{l} \in \boldsymbol{R}^{d_{l}\times d_{l-1}}$ is the strength of the synapse connections, also called the weight matrix, and $u_{\rm th} \in \boldsymbol{R}$ is the spike firing threshold
(Eq.~\ref{eqn:membrane-potential} for $l=1$ corresponds to the direct encoding~\cite{rueckauer2017conversion}).
\end{definition}

A backward mode of the single-step neuron as it is does not work in the standard backpropagation algorithm because the derivative of Eq.~\ref{eqn:spike} vanishes almost everywhere.
Therefore, we adopt the technique called surrogate gradient, {\textit{i.e.,}} we formally replace the derivative function with some reasonable function, for example, the following one~\citep{suetake2023s3nn}:
\begin{align}
\label{eqn:surrogate-s3nn}
\frac{\partial s_{l, i}}{\partial u_{l, i}}\left(u_{l, i}\right)
:\simeq \left\{
\begin{array}{ll}
\frac{1}{\tau} \frac{1}{u_{l, i}} & \left(u_{l, i} \geq u_{\rm th}\right), \\
\frac{\partial \sigma_\alpha}{\partial u_{l, i}}\left(u_{l, i}\right) & \left(u_{l, i} < u_{\rm th}\right),
\end{array}
\right.
\end{align}
where $\tau$ and $\alpha$ are hyperparameters and $\sigma_\alpha$ is the scaled sigmoid function expressed as follows:
\begin{align}
\sigma_\alpha\left(u_{l, i}\right)
&:= \frac{1}{1 + \exp\left(\left(-u_{l, i}+u_{\rm th}\right)/\alpha\right)}, \label{eqn:scaled-sigmoid}
{} \\
\frac{\partial \sigma_\alpha}{\partial u_{l, i}}\left(u_{l, i}\right)
&= \frac{1}{\alpha} \sigma_\alpha\left(u_{l, i}\right) \left(1 - \sigma_\alpha\left(u_{l, i}\right)\right). \label{eqn:scaled-sigmoid-derivative}
\end{align}
Note that the choice of a function for the surrogate function is irrelevant to our proposal.

\subsection{Metric for Energy Consumption}
\label{sec:metric}

We prepare the symbol $\psi_{l,i}$ for the number of synapses outgoing from the $i$-th neuron in the $l$-th layer, \textit{i.e.}, the number of matrix elements in $(\boldsymbol{W}_{l})_{*, i} \in \boldsymbol{R}^{d_l}$ that is not forced to vanish in terms of network architecture.
Let us denote $W_l, H_l$ as the width and height of the feature map, respectively, $C_l$ as the channel size, and $k_l$ as the kernel size in the $l$-th layer.
We restrict both the kernel width and height to be identical to $k_l$ for the sake of simplicity.
Then, explicit forms of $\psi_{l,i}$ are as follows, \textit{e.g.}, the standard fully connected (fc) and two-dimensional convolutional (conv) layers,
\begin{align}
\psi_{l,i} = \psi_{l}
= \left\{
\begin{array}{ll}
C_{l+1} & {\rm (fc)}, \\
\frac{W_{l+1} H_{l+1}}{W_l H_l} C_{l+1} k_l^2 & {\rm (conv)}, \\
\end{array}
\right.
\end{align}
where the convolutional layer assumes appropriate padding.
Note that some non-trivial padding, stride, and dilation can induce $i$-dependency of $\psi_{l, i}$.

Using $\psi_{l,i}$, we can express the number of floating point operations (FLOPs), which is often used as a metric to measure the computational complexity in ANNs, as follows:
\begin{align}
\label{eqn:flops}
{\text {FLOPs}}(l) &:= \sum_{i=1}^{d_l} \psi_{l, i},
\end{align}
and the layer-wise and total spike firing rates, which are also important metrics to measure the sparsity of spiking activity in SNNs, as follows:
\begin{align}
\label{eqn:spike-firing-rate-layer}
R(l) &:= \mathop{\mathbb{E}}\limits_{x \in X} \left[
\frac{\sum_{i=1}^{d_l}s_{l, i}}{d_l} \right], \\
\label{eqn:spike-firing-rate-balance}
R &:= \sum_{l=1}^L R(l),
\end{align}
where the operation $\mathop{\mathbb{E}}_{x \in X}$ means taking the empirical expectation in the dataset $X$.
Then, the energy consumption metric that we should optimize is defined as follows.
\begin{definition}
Let us denote $T$ as the size of time steps and $E_{\text{AC}}$ as the energy consumption per accumulate operation.
Then, the layer-wise and total energy consumption metrics for the SNN are defined as follows:
\begin{align}
\label{eqn:energy-consumption-layer}
E_{\text {SNN}}(l) &:= T E_{\text {AC}} \mathop{\mathbb{E}}\limits_{x \in X} \left[
\sum_{i=1}^{d_l}\psi_{l,i}s_{l, i} \right], \\
\label{eqn:energy-consumption-total}
E_{\text {SNN}} &:= \sum_{l=1}^{L} E_{\text {SNN}}(l).
\end{align}
\end{definition}
Note that $T$ is equal to one for the single-step neuron model, and we use $E_{\text {AC}} = 0.9$ [pJ]~\citep{horowitz20141}.

If $\psi_{l,i}$ is independent of $i$ ($\exists \psi_l, \forall i, \psi_{l,i}=\psi_{l}$), by combining Eqs.~\ref{eqn:flops} and \ref{eqn:spike-firing-rate-layer}, Eq.~\ref{eqn:energy-consumption-layer} is rewritten as follows:
\begin{align}
E_{\text {SNN}}(l)
&= T E_{\text {AC}} \sum_{i=1}^{d_l} \psi_{l} \mathop{\mathbb{E}}\limits_{x \in X} \left[
\frac{\sum_{i=1}^{d_l}s_{l, i}}{d_l} \right] \nonumber \\
&= T E_{\text {AC}} {\text {FLOPs}}(l) R(l),
\end{align}
which is the same metric as that used in \citet{kim2021revisiting}.
For the sake of simplicity, we consider the case where $\psi_{l, i}$ is independent of $i$ in Sec.~\ref{sec:experiment}.



\subsection{Spiking Synaptic Penalty}
\label{sec:proposal}

To optimize the energy consumption $E_{\text {SNN}}$ (Eq.~\ref{eqn:energy-consumption-total}), we propose the following penalty terms.

\begin{definition}
The layer-wise and total spiking synaptic penalty terms are defined as follows:
\begin{align}
\label{eqn:penalty-synops-layer}
\Omega_{\rm syn}(l) &= \Omega_{\rm syn}(l, \boldsymbol{s}_l)
:= \frac{1}{p} \sum_{i=1}^{d_l} \psi_{l,i} s_{l,i}^p, \\
\label{eqn:penalty-synops}
\Omega_{\rm syn} &= \Omega_{\rm syn}(\boldsymbol{s})
:= \sum_{l=1}^{L} \Omega_{\rm syn}(l, \boldsymbol{s}_l),
\end{align}
where $\boldsymbol{s} := \{\boldsymbol{s}_{l}\}_{l=1}^{L}$ and $p\geq1$.
\end{definition}

The equivalency between the total energy consumption metric and total spiking synaptic penalty immediately follows from their definitions and the equation $s_{l,i}^{p} = s_{l,i}$, which is derived from Eq.~\ref{eqn:spike}.
\begin{theorem}
\label{thm:main-theorem}
The expected value of the layer-wise and total spiking synaptic penalties are precisely proportional to the layer-wise and total energy consumption metrics of SNNs:
\begin{align}
p E_{\text{AC}} \mathop{\mathbb{E}}\limits_{x \in X}\left[ \Omega_{\rm syn}(l) \right] &= E_{\text {SNN}}(l), \\
\label{eqn:main-theorem}
p E_{\text{AC}} \mathop{\mathbb{E}}\limits_{x \in X}\left[ \Omega_{\rm syn} \right] &= E_{\text {SNN}},
\end{align}
for arbitrary $p\geq1$.
\end{theorem}
This fact means that optimizing Eq.~\ref{eqn:penalty-synops} leads to optimizing Eq.~\ref{eqn:energy-consumption-total}.
Hence, we strongly propose to use Eq.~\ref{eqn:penalty-synops} as the penalty term to optimize the energy consumption metric.
In the following,
we indicate the total spiking synaptic penalty when simply referred to as the \textit{spiking synaptic penalty} (including Eqs.~\ref{eqn:spike-firing-rate-balance} and \ref{eqn:energy-consumption-total}).

\begin{remark}
The proposed penalty term can be optimized in the manner of the surrogate gradient as in Sec.~\ref{sec:neuron-model},
and using $p \neq 1$ options controls the backward signal when the spike does not fire as follows:
\begin{align}
\label{eqn:penalty-backward}
\frac{1}{p}\frac{\partial s_{l, i}^p}{\partial u_{l, i}}\left(u_{l, i}\right)
&= s_{l, i}^{p-1} \frac{\partial s_{l, i}}{\partial u_{l, i}} \nonumber \\
&\simeq \left\{
\begin{array}{ll}
\frac{\partial s_{l, i}}{\partial u_{l, i}} & \left(u_{l, i} \geq u_{\rm th}\right), \\
0 & \left(u_{l, i} < u_{\rm th}\right),
\end{array}
\right.
\end{align}
where we used Eq.~\ref{eqn:spike}.
From Eqs.~\ref{eqn:penalty-synops-layer} and \ref{eqn:penalty-backward}, there is no intrinsic difference when $p > 1$; hence, we do not consider $p > 1$ options except for $p=2$, which is commonly used in several studies~\citep{esser2016cover, pellegrini2021low}.
Note that the derivation of Eq.~\ref{eqn:penalty-backward} is clearly independent of the choice of the surrogate gradient.
However, the open problem still remains, {\it i.e.,} it cannot be theoretically decided which choice is better, $p=1$ or $p>1$.
We experimentally examined it for $p=1, 2$ as described in Sec.~\ref{sec:experiment}.
\end{remark}

\subsubsection{Differences from Other Penalty Terms}
\label{sec:difference-from-other-penalty-terms}
The other candidates for the penalty term are as follows:
\begin{align}
\label{eqn:penalty-total}
\Omega_{\rm total} &= \frac{1}{p} \sum_{l=1}^{L} \sum_{i=1}^{d_l} s_{l,i}^p, \\
{}
\label{eqn:penalty-balance}
\Omega_{\rm balance} &= \frac{1}{p} \sum_{l=1}^{L} \sum_{i=1}^{d_l} \frac{1}{d_l} s_{l,i}^p,
\end{align}
where, for $p=2$, $\Omega_{\rm total}$ and $\Omega_{\rm balance}$ are the same as those in \citet{esser2016cover} and \citet{pellegrini2021low}, respectively.
However, we claim that neither can directly optimize the energy consumption because they do not have the proportional nature~(Eq.~\ref{eqn:main-theorem}), although they sparsify the spiking activity to some extent.

Figs.~\ref{fig:comparison-scaling-factor-cnn6-flops} and~\ref{fig:comparison-scaling-factor-vgg11-resnet18-flops} and Table~\ref{tbl:comparison-scaling-factor-metrics} show the discrepancy between the energy consumption metric and penalty terms of the model used in the following experiment.
In these figures and table, we assumed that all the spiking neurons fired, {\textit {i.e.}}, $s_{l,i}=1 (\forall l, i)$, for the sake of simplicity.
In this assumption, the ground truth is proportional to FLOPs without loss of generality.
These figures and table indicate that the spiking synaptic penalty is precisely proportional to the energy consumption metric, but other penalties are not.
In the next section, we will experimentally verify how this claim affects performance.

\begin{figure}[tb]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=\columnwidth]{experiment/comparison-scaling-factor-vgg11-flops.pdf}}
\subcaption{(A) VGG11}
\centerline{\includegraphics[width=\columnwidth]{experiment/comparison-scaling-factor-preresnet18-flops.pdf}}
\subcaption{(B) ResNet18}
\end{center}
\caption{Comparison among layer-wise penalty terms. The network architectures are (A) VGG11 and (B) ResNet18 (App.~\ref{app:detail-of-network-architecture}).
As in Fig.~\ref{fig:comparison-scaling-factor-cnn6-flops}, our penalty term (blue) is precisely proportional to the ground truth (gray).}
\label{fig:comparison-scaling-factor-vgg11-resnet18-flops}
\vskip -0.2in
\end{figure}

\begin{table}[tb]
\caption{Comparison among total penalty terms.
Our penalty term $\Omega_{\rm syn}$ is precisely equal to the ground truth $E_{\rm SNN}/E_{\rm AC}$.}
\label{tbl:comparison-scaling-factor-metrics}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{tabular}{lrrrr}
\toprule
Model & $E_{\rm SNN}/E_{\rm AC}$ & $\Omega_{\rm syn}$ & $\Omega_{\rm total}$ & $\Omega_{\rm balance}$ \\
\midrule
CNN7     &   98895888 &   98895888 &  59688 &  6 \\
VGG11    & 2526060544 & 2526060544 & 249856 & 10 \\
ResNet18 &  553730048 &  553730048 & 671744 & 20 \\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\subsubsection{Normalization of Penalty Terms}

Penalty terms are included in the objective function with their intensity parameter $\lambda$ as the coupling $\lambda \Omega_*$, where the symbol $*$ denotes ``syn", ``total", or ``balance".
For tractable treatment of the intensity parameter between various penalty terms or among models of various scales, we recommend normalizing the penalty terms by $\Omega_{*}(\boldsymbol{1})$, where $\boldsymbol{1}$ indicates $\boldsymbol{s}=\boldsymbol{1}$, {\it{i.e.}}, $s_{l,i}=1 \, \forall l, i$.
Note that replacing $\Omega_{*}$ with $\Omega_{*}/\Omega_{*}(\boldsymbol{1})$ is equivalent to replacing $\lambda$ with
\begin{align}
\label{eqn:normalized-lambda}
\lambda' = \lambda / \Omega_{*}(\boldsymbol{1}).
\end{align}
Hence, we sometimes adopt the normalized notation $\lambda'$ instead of $\lambda$.

\section{Experiment}
\label{sec:experiment}

In this section, we evaluate the effectiveness of the proposed spiking synaptic penalty.
First, we describe the setup for experiments.
Next, we show that the proposed method can decrease energy consumption.
Finally, we show that the proposed method can reduce the energy consumption more than other methods while maintaining the accuracy.
In particular, we show that the proposed method is independent of the choice of a function for the surrogate gradient and better than the conversion approach~\citep{sorbaro2020optimizing}.
Overall, the main objective is to analyze the behavior of our method rather than to achieve state-of-the-art performances.

\subsection{Experimental Setup}
\label{sec:experiment-setup}
As the single-step spiking neuron is developed for static tasks \citep{suetake2023s3nn}, we experimented
for the Fashion-MNIST~\citep{xiao2017fashion}, CIFAR-10, and CIFAR-100~\citep{krizhevsky2009learning} datasets widely used in SNN experiments ~\citep{esser2016cover,zhang2020temporal,chowdhury2021one}
with some network architectures, CNN7, VGG11, and ResNet18 (refer to App.~\ref{app:detail-of-network-architecture} for details).
For these experiments, we implemented the program by the PyTorch framework and used one GPU, an NVIDIA GeForce RTX 3090, with 24 GB (refer to Table~\ref{sec:difference-in-training-time} in the appendix for the difference in the training time).

In these experiments, we used the following objective function:
\begin{align}
\label{eqn:objective-total}
L = &\frac{1}{n}\sum_{n=1}^N \left(
\mathop{CE}(f(x_n), t_n) 
+ \lambda \Omega_{\rm syn}(\boldsymbol{s}(x_n)) \right) \nonumber \\
&+ \lambda_{\rm WD} \left(\| \boldsymbol{W}\|_{L_2}^2 
+ B_{\rm BN} \|\boldsymbol{W}_{\rm BN}\|_{L_2}^2 \right), 
\end{align}
where $(x_n, t_n)$ denotes the pair of input data and its label, $f$ denotes some spiking neural network, $\mathop{CE}$ denotes the cross-entropy function, $\| \boldsymbol{W}\|_{L_2}^2$ denotes the $L_2$ penalty for the weights, {\textit{i.e.}}, weight decay, $\| \boldsymbol{W}_{\rm BN}\|_{L_2}^2$ denotes the $L_2$ penalty for the trainable parameters of batch normalization layers, $\lambda$ and $\lambda_{\rm WD}$ denote the intensity of penalties, and $B_{\rm BN} \in \{0, 1\}$.
Note that we included the weight decay into the objective function to verify its sparsifying effect \citep{yaguchi2018adam} in the context of SNNs.
In addition, we explicitly specified the weight decay for batch normalization layers because it was included in the default setting of the PyTorch framework~\citep{paszke2019pytorch} but not in \citet{yaguchi2018adam}.
The training of $f$ was done by the backpropagation algorithm with the surrogate gradient of Eq.~\ref{eqn:surrogate-s3nn} unless otherwise stated.
The optimizer was selected from the momentum SGD~(mSGD) or Adam to confirm the claims in \citet{yaguchi2018adam}.

Refer to App.~\ref{app:detail-of-experimental-setup} for further details of the experimental setup such as hyperparameters.

\subsection{Energy Reduction by Spiking Synaptic Penalty}
\label{sec:experiment-ablation}

\begin{table*}[tb]
\caption{Performance with respect to varying $\lambda, \lambda_{\rm WD}$, and $B_{\rm BN}$.
The network architecture is VGG11, the dataset is CIFAR-10, and the optimizer is Adam.
$E_{\rm SNN}$ denotes the energy consumption metric for the SNN (Eq.~\ref{eqn:energy-consumption-total}),
$E_{\rm baseline}$ denotes the $E_{\rm SNN}$ for the model with $\lambda = \lambda_{\rm WD} = B_{\rm BN} = 0$,
dead rate denotes the ratio of the number of dead neurons to the number of total neurons,
and $R$ denotes the spike firing rate (Eq.~\ref{eqn:spike-firing-rate-balance}).}
\label{tbl:ablation-vgg11-adam}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{tabular}{ccccccc}
\toprule
$\lambda$ & $\lambda_{\rm WD}$ & $B_{\rm BN}$ & Accuracy [\%] & $E_{\rm SNN}/E_{\rm baseline}$ [\%] & Dead rate [\%] & $R$ [\%] \\
\midrule
    0 &     0 & 0 & 89.09 $\pm$ 0.271 & 100.0 $\pm$ 1.530 & 3.367 $\pm$ 0.221 & 12.11 $\pm$ 0.114 \\
1e-08 &     0 & 0 & 89.23 $\pm$ 0.141 & 73.92 $\pm$ 0.247 & 5.935 $\pm$ 0.822 & 10.34 $\pm$ 0.050 \\
1e-07 &     0 & 0 & 87.81 $\pm$ 0.458 & 34.07 $\pm$ 0.397 & 22.42 $\pm$ 1.005 & 6.046 $\pm$ 0.061 \\
1e-06 &     0 & 0 & 79.75 $\pm$ 0.694 & 11.85 $\pm$ 0.492 & 59.01 $\pm$ 1.174 & 2.176 $\pm$ 0.040 \\
1e-05 &     0 & 0 & 21.81 $\pm$ 14.46 & 1.053 $\pm$ 1.246 & 92.43 $\pm$ 3.149 & 0.174 $\pm$ 0.198 \\
\midrule
    0 & 1e-04 & 0 & 89.15 $\pm$ 0.493 & 81.07 $\pm$ 4.497 & 9.567 $\pm$ 1.005 & 8.474 $\pm$ 0.771 \\
    0 & 1e-03 & 0 & 89.00 $\pm$ 0.832 & 68.64 $\pm$ 6.714 & 18.72 $\pm$ 1.942 & 6.802 $\pm$ 0.492 \\
    0 & 1e-02 & 0 & 84.98 $\pm$ 1.455 & 44.39 $\pm$ 3.540 & 34.13 $\pm$ 0.683 & 4.992 $\pm$ 0.448 \\
    0 & 1e-01 & 0 & 50.51 $\pm$ 36.24 & 17.90 $\pm$ 16.03 & 70.93 $\pm$ 25.66 & 2.988 $\pm$ 2.642 \\
    0 & 1e+00 & 0 & 10.00 $\pm$ 0.000 & 0.000 $\pm$ 0.000 & 100.0 $\pm$ 0.000 & 0.000 $\pm$ 0.000 \\
1e-08 & 1e-03 & 0 & 89.61 $\pm$ 0.026 & 60.28 $\pm$ 1.384 & 20.20 $\pm$ 2.268 & 6.244 $\pm$ 0.155 \\
1e-07 & 1e-03 & 0 & 88.06 $\pm$ 0.653 & 28.72 $\pm$ 0.448 & 35.58 $\pm$ 2.255 & 4.076 $\pm$ 0.064 \\
1e-06 & 1e-03 & 0 & 76.00 $\pm$ 0.962 & 7.829 $\pm$ 0.076 & 68.02 $\pm$ 0.504 & 1.378 $\pm$ 0.016 \\
1e-05 & 1e-03 & 0 & 10.02 $\pm$ 0.029 & 0.007 $\pm$ 0.012 & 99.84 $\pm$ 0.283 & 0.001 $\pm$ 0.002 \\
\midrule
    0 & 1e-05 & 1 & 88.99 $\pm$ 0.626 & 96.29 $\pm$ 2.003 & 6.268 $\pm$ 0.786 & 9.753 $\pm$ 0.246 \\
    0 & 1e-04 & 1 & 89.87 $\pm$ 0.526 & 85.64 $\pm$ 3.412 & 10.24 $\pm$ 1.735 & 7.457 $\pm$ 0.293 \\
    0 & 1e-03 & 1 & 86.91 $\pm$ 0.315 & 51.06 $\pm$ 1.277 & 32.79 $\pm$ 1.257 & 4.376 $\pm$ 0.138 \\
    0 & 1e-02 & 1 & 10.00 $\pm$ 0.000 & 0.000 $\pm$ 0.000 & 100.0 $\pm$ 0.000 & 0.000 $\pm$ 0.000 \\
1e-08 & 1e-04 & 1 & 89.41 $\pm$ 0.147 & 45.87 $\pm$ 0.202 & 18.38 $\pm$ 0.631 & 5.643 $\pm$ 0.019 \\
1e-07 & 1e-04 & 1 & 86.76 $\pm$ 0.957 & 19.21 $\pm$ 0.610 & 41.97 $\pm$ 2.169 & 2.911 $\pm$ 0.043 \\
1e-06 & 1e-04 & 1 & 74.38 $\pm$ 0.465 & 6.582 $\pm$ 0.264 & 71.21 $\pm$ 0.608 & 1.111 $\pm$ 0.028 \\
1e-05 & 1e-04 & 1 & 10.00 $\pm$ 0.000 & 0.003 $\pm$ 0.005 & 99.92 $\pm$ 0.140 & 0.001 $\pm$ 0.001 \\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.1in
\end{table*}

\begin{table*}[tb]
\caption{Performance with respect to varying $\lambda, \lambda_{\rm WD}$, and $B_{\rm BN}$.
The optimizer is mSGD.
The remaining descriptions are the same as those in Table~\ref{tbl:ablation-vgg11-adam}.}
\label{tbl:ablation-vgg11-msgd}
\begin{center}
\begin{small}
\begin{tabular}{ccccccc}
\toprule
$\lambda$ & $\lambda_{\rm WD}$ & $B_{\rm BN}$ & Accuracy [\%] & $E_{\rm SNN}/E_{\rm baseline}$ [\%] & Dead rate [\%] & $R$ [\%] \\
\midrule
    0 &     0 & 0 & 89.14 $\pm$ 0.158 & 100.0 $\pm$ 2.548 & 1.678 $\pm$ 0.167 & 14.19 $\pm$ 0.277 \\
1e-09 &     0 & 0 & 89.36 $\pm$ 0.159 & 78.84 $\pm$ 1.166 & 2.428 $\pm$ 0.188 & 13.35 $\pm$ 0.214 \\
1e-08 &     0 & 0 & 88.81 $\pm$ 0.224 & 40.00 $\pm$ 0.335 & 11.04 $\pm$ 0.139 & 10.07 $\pm$ 0.097 \\
1e-07 &     0 & 0 & 86.80 $\pm$ 0.348 & 18.58 $\pm$ 0.384 & 33.62 $\pm$ 0.875 & 4.476 $\pm$ 0.143 \\
1e-06 &     0 & 0 & 24.32 $\pm$ 7.158 & 73.27 $\pm$ 98.19 & 78.72 $\pm$ 21.06 & 4.207 $\pm$ 4.165 \\
\midrule
    0 & 1e-04 & 0 & 88.59 $\pm$ 0.413 & 84.40 $\pm$ 6.843 & 2.705 $\pm$ 0.754 & 12.73 $\pm$ 0.344 \\
    0 & 1e-03 & 0 & 88.60 $\pm$ 2.134 & 70.96 $\pm$ 8.041 & 3.968 $\pm$ 1.326 & 11.97 $\pm$ 0.582 \\
    0 & 1e-02 & 0 & 70.80 $\pm$ 2.795 & 45.18 $\pm$ 3.008 & 19.57 $\pm$ 2.636 & 10.85 $\pm$ 0.651 \\
    0 & 1e-01 & 0 & 20.23 $\pm$ 17.72 & 18.52 $\pm$ 13.99 & 92.18 $\pm$ 10.94 & 3.322 $\pm$ 4.056 \\
1e-09 & 1e-03 & 0 & 90.96 $\pm$ 0.202 & 64.62 $\pm$ 0.802 & 3.976 $\pm$ 0.196 & 11.82 $\pm$ 0.085 \\
1e-08 & 1e-03 & 0 & 90.06 $\pm$ 0.387 & 35.65 $\pm$ 0.076 & 13.51 $\pm$ 0.195 & 9.043 $\pm$ 0.068 \\
1e-07 & 1e-03 & 0 & 87.74 $\pm$ 0.257 & 16.32 $\pm$ 0.166 & 37.68 $\pm$ 0.598 & 4.163 $\pm$ 0.040 \\
1e-06 & 1e-03 & 0 & 21.42 $\pm$ 9.467 & 26.13 $\pm$ 24.14 & 87.16 $\pm$ 4.910 & 3.255 $\pm$ 1.576 \\
\midrule
    0 & 1e-07 & 1 & 89.30 $\pm$ 0.263 & 99.75 $\pm$ 2.576 & 1.798 $\pm$ 0.260 & 14.17 $\pm$ 0.270 \\
    0 & 1e-06 & 1 & 88.85 $\pm$ 0.277 & 94.59 $\pm$ 2.186 & 1.838 $\pm$ 0.117 & 13.63 $\pm$ 0.108 \\
    0 & 1e-05 & 1 & 88.80 $\pm$ 0.235 & 95.85 $\pm$ 1.806 & 1.834 $\pm$ 0.134 & 13.43 $\pm$ 0.299 \\
    0 & 1e-04 & 1 & 88.51 $\pm$ 0.295 & 87.16 $\pm$ 1.614 & 2.657 $\pm$ 0.105 & 10.12 $\pm$ 0.151 \\
    0 & 1e-03 & 1 & 84.47 $\pm$ 0.757 & 51.21 $\pm$ 2.576 & 16.90 $\pm$ 0.388 & 5.005 $\pm$ 0.310 \\
    0 & 1e-02 & 1 & 19.53 $\pm$ 16.50 & 1.786 $\pm$ 3.093 & 97.64 $\pm$ 4.088 & 0.164 $\pm$ 0.284 \\
1e-09 & 1e-07 & 1 & 89.32 $\pm$ 0.417 & 78.35 $\pm$ 1.694 & 2.606 $\pm$ 0.078 & 13.24 $\pm$ 0.143 \\
1e-08 & 1e-07 & 1 & 88.77 $\pm$ 0.181 & 39.70 $\pm$ 0.844 & 11.29 $\pm$ 0.428 & 10.04 $\pm$ 0.173 \\
1e-07 & 1e-07 & 1 & 86.95 $\pm$ 0.236 & 18.60 $\pm$ 0.285 & 33.49 $\pm$ 0.418 & 4.478 $\pm$ 0.130 \\
1e-06 & 1e-07 & 1 & 28.05 $\pm$ 1.125 & 78.52 $\pm$ 105.4 & 81.75 $\pm$ 13.47 & 4.500 $\pm$ 4.411 \\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\end{table*}

We investigated whether optimizing the spiking synaptic penalty (Eq.~\ref{eqn:penalty-synops}) led to optimizing the energy consumption metric (Eq.~\ref{eqn:energy-consumption-total}) and whether there was any conflict with other terms in the objective function (Eq.~\ref{eqn:objective-total}).
The setting was as follows.
The baseline model was trained for Eq.~\ref{eqn:objective-total} with $\lambda = \lambda_{\rm WD} = 0$.
The other models were trained from scratch with some combinations of the weight decay ($\lambda_{\rm WD} > 0$), $L_2$ penalty for batch normalization layers ($B_{\rm BN} = 1$), and $p=1$ spiking synaptic penalty ($\lambda > 0$).
The results are presented in Tables~\ref{tbl:ablation-vgg11-adam} and \ref{tbl:ablation-vgg11-msgd}.
Note that the values in this table were taken for $\lambda$ and $\lambda_{\rm WD}$ from a point where the accuracy was very low (approximately $20\%$) until the accuracy reached the upper bound and stopped changing.

From the result in Table~\ref{tbl:ablation-vgg11-adam}, we can observe the following.
First, as the intensity of the penalty term increases, the energy consumption metric decreases; the inference accuracy also decreases. Therefore, the intensity parameter $\lambda$ controls the trade-off between them.
Second, the combination of the spiking synaptic penalty and weight decay further reduces the energy consumption metric.
Therefore, we propose to adopt both of them simultaneously.
In addition, we found that the combination of the weight decay and Adam optimizer induces neuron sparsification even without the spiking synaptic penalty, though its contribution to the energy reduction is less than the spiking synaptic penalty.
Furthermore, the neuron sparsification proceeds more strongly for the Adam optimizer than the mSGD optimizer, which is consistent with the claim in~\citet{yaguchi2018adam} (see also Table~\ref{tbl:ablation-vgg11-msgd}).
Note that we cannot observe a remarkable difference between $B_{\rm BN} = 0$ and $1$.
Therefore, we adopt the weight decay with $B_{\rm BN} = 0$ in further experiments to simplify our objective function.
Finally, all the above results hold for not only VGG11 but also CNN7 and ResNet18 (see Tables~ \ref{tbl:ablation-cnn6-adam}--\ref{tbl:ablation-resnet18-msgd} in the appendix).

\subsection{Trade-off between Accuracy and Energy Efficiency}
\label{sec:experiment-trade-off}

\subsubsection{Comparison between Penalties}
\label{sec:comparison-between-penalties}
To see how the difference between various penalty terms affected the training result, we conducted a comparative experiment.
The setting was as follows.
For fair comparison, we used the $\lambda'$ notation for the intensity parameter of penalties rather than the raw $\lambda$ (see Eq.~\ref{eqn:normalized-lambda}).
The baseline model was trained for Eq.~\ref{eqn:objective-total} with $\lambda' = B_{\rm BN} = 0$, and we tuned $\lambda_{\rm WD} > 0$ to obtain the highest accuracy.
Then, the others were trained by varying $\lambda' > 0$ and by replacing $\Omega_{\rm syn}$ in Eq.~\ref{eqn:objective-total} with $\Omega_{\rm total}$ or $\Omega_{\rm balance}$ from scratch.
The results are shown in Fig.~\ref{fig:comparison-scaling-curves} (A) as $\lambda'$-parameterized curves of the energy--accuracy trade-off,
where the energy consumption rate was produced as the energy consumption of each model normalized by that of the baseline model.
Note that it is better for data to be located at the upper left corner in the figure.
Refer to App.~\ref{app:detail-of-hyperparameter} for the sampling of $\lambda'$.
In addition, the quantitative analysis is presented in Table~\ref{tbl:quantitative-comparison-fmnist-cnn6-adam}, where higher scores are better for all the metrics: area under the curve (AUC), Spearman's rank correlation coefficient (Spearman), and the mutual information (MI).
Note that the argument of each metric represents a cutoff parameter, where data with lower accuracy than it are omitted.
We introduced the cutoff parameter because training tended to break as the intensity parameter was increased for all the methods.
Refer to App.~\ref{app:detail-of-metrics} for details of quantitative metrics.

\begin{figure*}[t]
\centering
\begin{minipage}[b]{0.45\linewidth}
    \centering
   
    \includegraphics[width=\linewidth]{experiment/comparison-scaling-factor-cnn6-adam.pdf}\\
   
    \subcaption{(A) Eq.~\ref{eqn:penalty-synops}}
\end{minipage}
\begin{minipage}[b]{0.45\linewidth}
    \centering
    \includegraphics[width=\linewidth]{experiment/comparison-surrogates-triangle-cnn6-adam.pdf}\\
   
   
    \subcaption{(B) Eq.~\ref{eqn:surrogate-triangle}}
\end{minipage}
\\
\begin{minipage}[b]{0.45\linewidth}
    \centering
    \includegraphics[width=\linewidth]{experiment/comparison-surrogates-scaledsigmoid-cnn6-adam.pdf}\\
   
   
    \subcaption{(C) Eq.~\ref{eqn:surrogate-scaled-sigmoid}}
\end{minipage}
\begin{minipage}[b]{0.45\linewidth}
    \centering
    \includegraphics[width=\linewidth]{experiment/comparison-ann2snn-cnn6-adam.pdf}\\
   
   
    \subcaption{(D) Conversion}
\end{minipage}
\caption{Energy--accuracy trade-off curves.
The network architecture is CNN7, the dataset is Fashion-MNIST, and the optimizer is Adam.
The energy consumption rate is the energy consumption of each model normalized by that of the baseline model.
From (A) to (C): each model is trained by the indicated surrogate gradient below each figure.
(D): comparison with the conversion approach~\citep{sorbaro2020optimizing}. $\Omega_{\rm ANN}$ denotes SNNs converted from the QReLU network.
The $x$-axis exceeds one because converted SNNs have higher energy consumption than the SNN baseline for $\Omega_{\rm syn}$.}
\label{fig:comparison-scaling-curves}
\end{figure*}

\begin{table*}[tb]
\caption{Quantitative comparison corresponding to Fig.~\ref{fig:comparison-scaling-curves}~(A). Higher scores are better. The best and the second-best results are highlighted in bold and underlined, respectively.
Refer to App.~\ref{app:detail-of-metrics} for details of the quantitative metrics.
}
\label{tbl:quantitative-comparison-fmnist-cnn6-adam}
\begin{center}
\begin{small}
\scalebox{1.0}{
\begin{tabular}{lcccccc}
\toprule
Method & AUC$(70)$[$\%$] & AUC$(50)$[$\%$] & Spearman$(70)$ & Spearman$(50)$ & MI$(70)$ & MI$(50)$ \\
\midrule
$\Omega_{\rm syn}\,(p=1)$ (Ours)         & \textbf{68.02} & \textbf{79.60} & \textbf{0.9861} & \textbf{0.9865} & \textbf{3.465} & \textbf{3.610} \\
$\Omega_{\rm syn}\,(p=2)$ (Ours)         & 61.62 & 72.69 & 0.9474 & 0.9709 & 3.233 & \underline{3.476} \\
$\Omega_{\rm total}\,(p=1)$              & \underline{63.30} & \underline{76.05} & 0.9766 & \underline{0.9767} & \underline{3.244} & 3.319 \\
$\Omega_{\rm total}\,(p=2)$              & 55.16 & 64.63 & \underline{0.9831} & 0.9701 & 3.218 & 3.295 \\
$\Omega_{\rm balance}\,(p=1)$            & 54.23 & 67.16 & 0.9412 & 0.9412 & 2.978 & 2.978 \\
$\Omega_{\rm balance}\,(p=2)$            & 31.47 & 42.94 & 0.8500 & 0.8946 & 2.708 & 2.833 \\
\bottomrule
\end{tabular}
}
\end{small}
\end{center}
\end{table*}




From the result in Fig.~\ref{fig:comparison-scaling-curves} (A) and Table~\ref{tbl:quantitative-comparison-fmnist-cnn6-adam}, we can observe the following.
First, for each $\Omega_*$, the $p=1$ option is apparently better than the $p=2$ option. Therefore, we propose to adopt the $p=1$ option.
Note that this difference arises from the backward control as Eq.~\ref{eqn:penalty-backward}.
However, the mechanism by which it affects the training result is still unclear; it remains our future work.
Second, for $p = 1$, the trade-off curve of $\Omega_{\rm syn}$ is the best, followed in order by $\Omega_{\rm total}$ and $\Omega_{\rm balance}$.
Therefore, we experimentally clarified the advantage of the coefficient $\psi_{l,i}$ for Eq.~\ref{eqn:penalty-synops-layer}, which had remained an issue in the method proposed by \citet{sorbaro2020optimizing}.
Finally, all the above results hold for not only CNN7 but also VGG11 and ResNet18 (see Fig.~\ref{fig:comparison-scaling-factor-merged} and Tables~\ref{tbl:quantitative-comparison-fmnist-cnn6-msgd}--\ref{tbl:quantitative-cifar100-resnet18-adam} in the appendix).

\subsubsection{Independence of Surrogate Gradients}
\label{sec:experiment-surrogate-gradients}
To see how the choice of a function for the surrogate gradient affected the training result, we conducted the same experiment as that in Sec.~\ref{sec:comparison-between-penalties} except for the choice of a function.
Instead of Eq.~\ref{eqn:surrogate-s3nn}, we adopted the piece-wise linear function \citep{esser2016cover} and scaled sigmoid \citep{pellegrini2021low} function for the surrogate gradient as follows:
\begin{align}
\label{eqn:surrogate-triangle}
\frac{\partial s}{\partial u} &\simeq \max{(1-|u-u_{\rm th}|, 0)},\\
\label{eqn:surrogate-scaled-sigmoid}
\frac{\partial s}{\partial u} &\simeq \frac{\partial \sigma_\alpha}{\partial u},
\end{align}
where $\sigma_\alpha$ is the same as Eq.~\ref{eqn:scaled-sigmoid-derivative}.
The results are presented in Figs.~\ref{fig:comparison-scaling-curves} (B) and (C), and Tables~\ref{tbl:quantitative-comparison-triangle} and \ref{tbl:quantitative-comparison-sigmoid} in the appendix.

From the result in Figs.~\ref{fig:comparison-scaling-curves} (B) and (C), the same observations as those in Sec.~\ref{sec:comparison-between-penalties} hold.
That is, the $p=1$ option is apparently better than the $p=2$ option; the trade-off curve of $\Omega_{\rm syn}$ is the best, followed in order by $\Omega_{\rm total}$ and $\Omega_{\rm balance}$.
Therefore, the spiking synaptic penalty works independent of the choice of a function for the surrogate gradient.

\subsubsection{Superiority to Conversion Approach}
\label{sec:experiment-conversion}

To see how the difference between training methods affected the training result, we also produced the trade-off curve for the conversion approach~\citep{sorbaro2020optimizing}.
We trained a single QReLU network (an ANN with quantized ReLU activations) increasing the intensity of the penalty and evaluated the converted SNNs for each intensity in different time steps: $T=1, 5$, and $10$ (refer to the original paper~\citep{sorbaro2020optimizing} for details).
The results are presented in Fig.~\ref{fig:comparison-scaling-curves} (D) and Table~\ref{tbl:quantitative-comparison-ann2snn} in the appendix, where the energy consumption was normalized by that of the baseline for $\Omega_{\rm syn}$.

From the result in Fig.~\ref{fig:comparison-scaling-curves} (D), we can observe that both the energy consumption and accuracy for the conversion approach are worse than those for the surrogate gradient approach.
This is because the conversion process degrades the accuracy, and the penalty term for the QReLU network cannot directly optimize the energy consumption metric for the SNN.
Hence, we should directly train SNNs by the surrogate gradient and spiking synaptic penalty to avoid such degradation.


\section{Conclusion}
\label{sec:discussion}
We studied the training method to obtain energy-efficient SNNs in terms of the surrogate gradient.
Based on our principle that we should optimize the metric as it is, we derived the spiking synaptic penalty to optimize the energy consumption metric.
Then, we experimentally showed that the spiking synaptic penalty (especially for $p=1$) is superior to the existing penalties and conversion approach.
Furthermore, its effectiveness is indifferent to the type of network architectures and surrogate gradients.
We conclude that our principle has worked well.

An apparent limitation is that the definition of the spiking synaptic penalty depends on that of the energy consumption metric.
However, if the target metric becomes deformed, the penalty should be accordingly deformed in accordance with our principle---even though it is a metric irrelevant to the energy consumption.
Another limitation is that although the target metric is directly included in the objective function, it is just indirectly optimized by the surrogate gradient.
To overcome this issue, we need an alternative training method such as equilibrium training~\citep{xiao2021training}.

We further list some outstanding issues.
First, it is unclear why there was a difference between the training result for $p=1$ and $2$ of the spiking synaptic penalty.
Elucidating the mechanism of this difference could help us understand the surrogate gradient.
Second, we did not focus on the synergy between the spike sparsification and pruning.
A pruning-aware sparsification training will help us obtain more energy-efficient SNNs.
Finally, the high availability of the spiking synaptic penalty should be verified, for example, in the case of real datasets, large networks, and other tasks.
By solving these issues, we can contribute to the realization of genuinely eco-friendly SNNs.

To facilitate further research, we will release the code shortly.




