\section{Introduction}
Large-scale online experiments allow researchers to recruit diverse populations from around the world, thus reducing the ``WEIRD" (Western, Educated, Industrial, Rich, Democratic) recruiting bias in cognitive science \cite{barrett2020crosscultural, henrich2010weird}. However, in online experiments, there is often less control on recruiting compared with a lab study. Thus it is critical to quickly and accurately assess the language proficiency of participants and not, which is common research practice, blindly trust self-report questionnaires \cite{lemhoefer2011LexTALE}.

Language proficiency is usually measured using vocabulary tests, which test the participants' knowledge of relatively uncommon words.
In contrast to the majority of vocabulary tests, which are designed in the context of second-language education \cite{alderson2005DIALANG}, the LexTALE test \cite{lemhoefer2011LexTALE} has been developed to meet the needs of behavioral lab experiments. In this short pen-and-paper test, participants are asked to go through the list and determine if each test item is a real word or a pseudoword (i.e. not a real word). Real words are infrequent words (occurring approximately 1--20$\times$ per million) which are assumed to be only recognized by highly proficient speakers. Pseudowords are created by changing letters in existing words or by recombining existing morphemes. In this test, participants should obtain a higher accuracy if they are more proficient in the language. This procedure makes the test faster than educational vocabulary tests, which can take up to 45 minutes to complete \cite{university2001qpt}.

 The test was originally developed for English, Dutch, and German, and was later adapted to Chinese \cite{chan2018LEXTALECH}, Italian \cite{amenta2020LexITA}, Finnish \cite{salmela2021Lexize}, French \cite{brysbaert2013LextaleFR}, Spanish \cite{izura2014LextaleEsp}, and Portuguese \cite{zhou2021LextPT}. Extending the test to another language however, requires extensive manual work and experimental validation.

While LexTALE was developed for laboratory experiments, it does not work well in online experiments. For example, LexTALE uses a fixed presentation order (participants who participate in multiple online experiments can memorize the order); it presents twice as many words as pseudowords (participants can exploit this imbalance); and items are not presented with a fixed pace (without time limitation participants can use an online dictionary search to look up the word). Here, we address these issues and additionally overcome two fundamental limitations in the LexTALE design. LexTALE relies on (1) the existence of curated word frequency databases and (2) the domain knowledge of the creator, which can be up to a certain degree of subjectivity (e.g., which linguistic rules were applied for the creation of pseudowords). These limitations have likely contributed to the fact that LexTALE is only available in a small set of languages and consists of a small number of words, which makes repeated participation impossible due to memory effects.

\begin{figure*}[ht!]
    \begin{center}
        \includegraphics{figures/approach.pdf}
    \end{center}
    \caption{Stimulus creation and study design. \textbf{A} Real words were extracted from text found on Wikipedia. Each token was lemmatized and spellchecked. \textbf{B} From all accepted lemmas in a language, we compute the n-gram transitional probabilities of letters within words. Pseudowords are created by drawing samples from the transitional probability. Here for demonstration purposes, we use 3-grams without * padding. Character-based languages, such as Chinese, are first converted into letter-based representations (e.g., Pinyin) to compute the transitional probabilities. The letters are later converted back to character sequences.  \textbf{C} About 40 monolingually-raised participants from the United Kingdom, Germany, the Netherlands, Spain, France, and Italy participate in a language test in their native (L1) and in a foreign language (L2). \textbf{D} Interface for an Italian participant. Participants use their keyboard for quick responses.}
    \label{fig:approach}
\end{figure*}

To alleviate these issues, we developed WikiVocab -- an automatic pipeline to generate vocabulary tests based on Wikipedia. Wikipedia is one of the largest collaborative internet projects and available in 329 languages\footnote{https://en.wikipedia.org/wiki/Wikipedia}. Because of this it, can be used to estimate vocabulary frequency in multiple languages. In our approach, we obtain a word frequency distribution from Wikipedia in a given language using automated heuristics. We sample rare words from this distribution, compute the n-gram transitional probabilities of letters within words, and create pseudowords that exhibit comparable n-gram probabilities. Since we obtain the word frequency distribution directly from a text corpus and implicitly capture orthographic regularities via n-gram probabilities, our approach does not rely on domain expertise nor on curated databases. The automatic pipeline allows us to create a very large pool of real words and pseudowords.

In what follows, we will describe how we generate the vocabulary test and benchmark it against available LexTALE tests in eight languages using six experiments with a total of 236 participants recruited from six countries. We show that our test can separate native speakers of closely related languages, that participants' performance in our test strongly correlates with LexTALE and self-reports, and finally, that participants' scores in different foreign languages reflect linguistic distances from these languages to their native language.

\section{Methods}
\subsection{Extract and collect texts}
\paragraph{Text processing} For each language, we process the longest 10,000 articles on Wikipedia (see Figure \ref{fig:approach}A), as the same articles in different languages fluctuate in their length. We use UDPipe 2.0 \cite{straka2018udpipe} to lemmatize all words and obtain a Part-Of-Speech tag (POS) to exclude proper nouns such as ``Eli Abbott" from the vocabulary test (see Figure \ref{fig:approach}A for the example). Since each Wikipedia page from a given language may contain words from multiple languages, we automatically check the original language of each word, using two approaches: (1) by using \texttt{fasttext} \cite{bojanowski2016fasttext} and (2) by using open-source spellcheckers from LibreOffice dictionaries\footnote{\url{https://github.com/LibreOffice/dictionaries}} using \texttt{guess\_language-spirit}\footnote{\url{https://pypi.org/project/guess\_language-spirit/}} and \texttt{pyenchant}\footnote{\url{https://github.com/pyenchant/pyenchant}}. A word is considered to be part of the associated language if either the spellchecker or \texttt{fasttext} correctly predicts that language. In addition to those filters, we add some language-agnostic preprocessing steps. For example, we reject words that are written in exclusively capital letters, exclude one-letter words, and remove words containing numbers or punctuation. From all words that fulfill these criteria, we only select the nouns to make the test more comparable across languages. In addition, we only use the lemmatized words to avoid multiple word forms of the same word.

\paragraph{Word filtering}
After we collect all the words per language, we convert all words to lowercase. Next, we remove words that exclusively occur in a small set of articles as these words tend to be jargon (e.g. `hippocampus') and are therefore not a good fit for a general-purpose vocabulary test. Thus, we calculate the ratio of word counts to the number of different articles in which the word occurs, and keep words in the 95th percentile, thus cutting off the top 5\%.

\paragraph{Compound word removal}
We train \texttt{charsplit} \cite{tuggener2016charsplit} on all the accepted lemmas in a language to detect potential compound words (for example ``snowball'' in English). We consider a word to be a compound word if there is a plausible word boundary between two words (i.e., $>$ 0; threshold proposed by \citeA{tuggener2016charsplit}) and if the last segment of the detected compound is a valid word in that language. The exclusion of compound words is relevant for many languages, as sampling from n-grams of compound words tends to generate even longer compound words, which are never used in the language, but do have a lexical meaning (for example ``Krankenversicherungsänderungsgesetzesentwurf" in German, which would mean: a draft of a law to change a health insurance provider). Another issue with compound words is that they pollute the word frequency distribution: compound words typically occur less often than each of their compounds, but this does not make the compound word harder to recognize than its elements (e.g., ``snowball" is not much harder than ``snow" and ``ball"). We therefore automatically remove all detected compound words. We could not apply compound word elimination to Chinese, as most multi-character words in Chinese would be considered compound words according to the heuristics above.

\subsection{Create pseudowords}
\paragraph{Convert to letters}
For Chinese, we use Pinyin\footnote{\url{https://github.com/mozillazg/python-pinyin}} to convert characters to sequences of letters from which we can obtain letters' n-grams (see Figure \ref{fig:approach}B). We store the mapping between characters and letters to convert the letters back to characters after sampling the n-grams.

\paragraph{Compute n-grams}
We compute n-grams of letter-based sequences from all words that are valid nouns and fulfill all filtering criteria. Existing linguistic work on spoken lexicons of multiple languages has shown that pseudowords generated from a 5-phone model capture most phonotactic regularities across the real words of the language \cite{dautriche2017_phonotactic_regularities, trott2020_phonotactic_regularities}. While, n-phone statistics on the phonetic transcriptions of words are not identical to n-gram statistics on written words\footnote{For example, a multi-letter grapheme $<$ph$>$ can refer to a single phoneme /f/, or certain graphemes are not pronounced such as the final $<$e$>$ in ``lak\textbf{e}".}, n-phones are similar to n-gram statistics. We, therefore, use 5-gram transitional probabilities to create pseudowords, since they are the closest equivalents to 5-phone transitional probabilities for the written language. To be able to create pseudowords that contain less than 5 letters, we pad each word with four asterisks on each side.

\paragraph{Sample from n-grams}
We start by sampling a 5-gram sequence which starts with four asterisks (`****') as the start symbols of a word. When sampling the next 5-gram, we make sure that the first four letters are identical to the last four letters of the previous 5-gram. We repeat this procedure until a 5-gram is drawn, which ends with an asterisk that indicates the end of the word. We remove the asterisks and validate the created string of letters. We keep on sampling until we create 1,000 unique pseudowords.

\paragraph{Validate pseudowords}
We reject generated pseudowords that correspond to real words in the language. We also reject pseudowords that contain too few or too many letters based on the range of word length of all real words of each language. For Chinese, we convert each pseudoword from the Pinyin representation back to the character representation. We do this by replacing all characters. To replace the longest letter sequences first, we sort the letter-character mapping by the length of the letter string. If there are left-over letters in the pseudowords, the word is rejected. For all other languages, we check if the created pseudowords are likely to be compound words. We reject the word if this is the case. To avoid the creation of pseudowords that look similar to existing words and are potential typos, we compute a fuzzy search using \texttt{thefuzz}\footnote{\url{https://github.com/seatgeek/thefuzz}}. Since the total number of lemmas and tokens is extremely large, we limit the search to words that start with the first and last three letters and are of a similar length (10\% difference in length allowed). We store the maximum match between the pseudowords and any word in that language.

\paragraph{Identify difficult words}
To match the task difficulty of LexTALE, we first look up the real LexTALE items in our word frequency distribution (log10-scale). We then  compute the mean and standard deviation of the LexTALE items per language. We select real words by sampling from a normal distribution centered at the average LexTALE word frequency with the standard deviation being computed over all languages with LexTALE.

\paragraph{Pair words}
Out of the 1,000 created pseudowords, we select 500 that best match the words in that language. To do so, we obtain the logarithm of the transitional probabilities of the letters' 5-grams for both the words and the pseudowords. We then computed the average absolute difference between the words and pseudowords that have the same number of 5-grams. On the resulting distances, we keep matching the word and pseudowords with the smallest distance. So, we match words and pseudowords that have a similar `rarity' at the same positions in the word. We repeat this procedure until we match all pseudowords. From the matched list, we only include 500 matched pseudowords that have the smallest fuzzy match ratio to any of the words in that language. By doing so, we select the 500 pseudowords which are least likely to be typos. For compatibility with LexTALE, we only include the first 60 items for the behavioral experiment.

\subsection{Experimental Design}
We selected eight languages to benchmark our test on (see Figure \ref{fig:approach}C). Six languages are Indo-European languages and can further be grouped into Germanic (English, \texttt{en}; German, \texttt{de}; Dutch, \texttt{nl}) and Italic (Spanish, \texttt{es}; French, \texttt{fr}; Italian, \texttt{it}) subfamilies. Chinese (\texttt{zh}) is a Sino-Tibetan language and Finnish (\texttt{fi}) is a Uralic language. We use a semi-balanced design, in which native speakers of the six Indo-European languages participate in the experiment. Each participant perform LexTALE and WikiVocab in their native language and in one of the other seven foreign languages selected at random. The order of the tests is also randomized. One real word and pseudoword for each language are presented in Table \ref{tab:example_items} as examples.\footnote{The list of all the stimuli can be downloaded here: \url{https://cogsci-23.s3.amazonaws.com/wikivocab.zip}.}

\begin{table}[ht]
\centering
\begin{tabular}{rll}
  \hline
 & real & pseudo \\
  \hline
  de &  fasanerie & fälschuld  \\
  nl & verrechtsing & rugzeeschaf \\
  en & carrion & washioneer  \\
  fr & envasement & paléonton \\
  it & sciocchezza & giustibola  \\
  es & demonio & tonomía \\
  fi & figuuri & antalaisu  \\
  zh & \begin{CJK*}{UTF8}{gbsn}经销\end{CJK*}  &\begin{CJK*}{UTF8}{gbsn}冰浑\end{CJK*}  \\
   \hline
\end{tabular}
\caption{Example items}
\label{tab:example_items}
\end{table}

\subsection{Participants}
To collect the data in the six countries, we ran online experiments on Prolific\footnote{\url{https://www.prolific.co}}. All texts in the interface of the experiment (e.g., buttons, instructions, etc) were presented in the native language of the participant. Non-English texts were first automatically translated using DeepL, then manually checked and corrected by a native speaker of each language. The recruitment and experimental pipelines were automated using PsyNet~\cite{harrison2020psynet}\footnote{\url{https://www.psynet.dev}}, a modern framework for experiment design and deployment which builds on the Dallinger\footnote{\url{https://dallinger.readthedocs.io/en/latest/}} platform for recruitment automation. Participants have to be born and currently reside in the targeted country, raised monolingually, hold nationality from that country, and speak the target language as their mother tongue. Overall, 236 participants with a mean age of 33 (SD = 12) complete the study. 61 \% of the participants have at least a Bachelor's degree. Participants were paid 9 GBP per hour and provided informed consent according to an approved protocol.

\subsection{Interface}
Words and pseudowords were presented in random order. To reduce the chance that a participant will search the word on the internet, we disabled the ability to copy the target word and limited the display time to two seconds. Participants were asked to respond as fast as possible by pressing two dedicated keys on their keyboard (see Figure \ref{fig:approach}D). To estimate the reliability of the test, participants did two batches of trials per language and test. Each batch contained 30 trials.

\begin{figure*}[ht!]
    \begin{center}
        \includegraphics{figures/figure_2.pdf}
    \end{center}
    \caption{Results. \textbf{A} Scatter plots of WikiVocab accuracy compared with LexTALE accuracy (left) and language self-report (right). \textbf{B} Consistency in the accuracy of the same participant for WikiVocab (left) and LexTALE (right). \textbf{C} WikiVocab accuracy on native language (i.e., L1), same language sub-family (e.g., L1 = Italian, L2 = Spanish), different sub-family (e.g., L1 = English, L2 = Spanish), and not Indo-European (e.g., L1 = German, L2 = Finnish). The dots represent measurements from all languages. All error bars are standard errors from the mean.}
    \label{fig:results}
\end{figure*}

\section{Results}
We computed the correlation between the average performance in the two blocks for WikiVocab and LexTALE for each language per participant. As shown in Figure \ref{fig:results}A the performance in LexTALE strongly correlates ($r = 0.86$, $p < .001$) with the performance in WikiVocab. The score in WikiVocab also strongly correlates with the language self-report of the participant ($r = 0.86$, $p < .001$). To measure the consistency on the test, we measured the correlation between performances on both blocks: $r = 0.82$ ($p < .001$) for WikiVocab (Figure \ref{fig:results}B) which is similarly high as the consistency measured for LexTALE ($r = 0.88$,  $p < .001$). In Figure \ref{fig:results}D, we show that the accuracy of the native language is much higher than for the other foreign languages.

To further benchmark the human performances on WikiVocab, we also ran the same test on GPT-3 (see Table \ref{tab:accuracies}), a large language model pre-trained on a massive corpus of language data~\cite{brown2020language}. To elicit GPT-3 responses we used the following prompt ``\textit{You will be shown a series of words and must determine if each one is a real or fake word in \{language\}. If it is a real word in \{language\} respond with 1. Otherwise, respond with 0.
Word: \{word\}
Response:}", where \{word\} was replaced by the target word and \{language\} by the target language.
GPT-3 outperformed native speakers in all languages except for Dutch and German. Both languages unite a strong tendency to create long compound words. In Dutch, 9/14 mistakes were made because, by misclassifying, the residual compound words were not detected as real words. Importantly, GPT3 performed nearly perfectly on what can be seen as its ``L1" (accuracy in English was 0.98) compared with lower performances on its ``L2s"  (performance in other languages ranges between 0.67-0.95). Moreover, languages that are likely to be well-represented in the training set of GPT-3 (es, fr and it) have relatively high scores (larger than 0.88) suggesting that ``language familiarity" also affects task performance for GPT-3.

\begin{table}[ht]
\centering
    \begin{tabular}{rlrrrrrr}
      \hline
    & DE & NL & GB & FR & IT & ES & GPT-3 \\
      \hline
      de & \textbf{0.90} & 0.69 & 0.52 & 0.53 & 0.60 & 0.53 & 0.85 \\
      nl & 0.55 & \textbf{0.84} & 0.65 & 0.54 & 0.54 & 0.56 & 0.77 \\
      en & 0.82 & 0.84 & \textbf{0.86} & 0.84 & 0.84 & 0.79 & 0.98 \\
      fr & 0.59 & 0.65 & 0.58 & \textbf{0.83} & 0.66 & 0.60 & 0.88 \\
      it & 0.61 & 0.51 & 0.60 & 0.66 & \textbf{0.92} & 0.70 & 0.92 \\
      es & 0.72 & 0.65 & 0.74 & 0.69 & 0.73 & \textbf{0.93} & 0.95 \\
      fi & NA & 0.51 & 0.51 & 0.53 & 0.50 & 0.45 & 0.73 \\
      zh & 0.57 & 0.54 & 0.55 & 0.57 & 0.56 & 0.52 & 0.67 \\
       \hline
    \end{tabular}
    \caption{Human vs. GPT-3 performance on WikiVocab. The rows are the tested languages, the columns are the country of the native speakers. Note: No German participant completed the Finnish WikiVocab task, so there is a missing data point.}
    \label{tab:accuracies}
\end{table}

In Figure \ref{fig:results_across}A, we show the average self-report proficiency and the performance on LexTALE and WikiVocab in the foreign and native languages. Interestingly, most participants score well on the English test and also indicate in the self-report they are proficient in English. Other than that, in all three cases, the diagonals (reflecting L1) are higher than non-diagonal items.

To assess whether participants perform better in their native language compared with a foreign language, we computed the difference between their performances (Figure \ref{fig:results_across}B). We found  that the gap between one's native language and the foreign language is largest for Dutch, whereas the difference in performance between native languages and English was lowest. The strong horizontal line in Figure \ref{fig:results_across}A for English shows that most participants in the test were fluent in English and the gap to L1 (the diagonal) was smallest. Nonetheless, for all languages, the performance in the native language was significantly higher than in the foreign language ($p<0.05$).

Finally, we correlated the mean accuracy on the test with linguistic distance between the tested and native language from \citeA{beaufils2020lexicaldistance}. Test score and linguistic distance are negatively correlated $r = -0.71$ ($p < 0.001$), indicating that the performance on the test is lower if the tested language is further away from one's native language (so German and Dutch are ``closer'' in linguistic distance compared with Dutch and Finnish).

\begin{figure*}[ht!]
    \begin{center}
        \includegraphics{figures/figure_3.pdf}
    \end{center}
    \caption{Results across languages. \textbf{A} Heatmaps of the mean self-report (left), LexTALE accuracy (center), and WikiVocab accuracy (right). \textbf{B} For each of the native languages (x-axis), we computed the difference in performance compared to the tested foreign language. The asterisks indicate the level of significance of a paired t-test. Significance codes: * = $p < 0.05$, *** = $p < 0.001$. Error bars represent the standard error of the mean. \textbf{C} Scatter plot of linguistic distance between the tested and native language and mean WikiVocab accuracy.}
    \label{fig:results_across}
\end{figure*}
\section{Discussion}
In this paper, we propose a new approach, ``WikiVocab'', to create automated vocabulary tests for online research. We generated the test for eight languages and showed that our test allows us to dissect between native speakers of closely related languages and that the performance strongly correlates with existing tests and self-reports. When looking into the differences between the languages, we showed that the test performances correlate with linguistic distance and that the native language performance is significantly higher than the same participants' performance in foreign languages.

Our results show that most participants are highly proficient in English, in their self-reports and on both tests. One possible explanation for this behavior might be the relatively high degree of formal education of the participants ($>$ 60\% have at least a Bachelor's degree) and the fact that most experiments on Prolific are in English. To increase the contrast between native and second language learners of English, one could further decrease the frequency of occurrence of real words which makes them more difficult to recognize.

The pairwise comparisons between the countries also reveal that the findings are strongly directional. Whereas the Dutch participants almost reach 70\% accuracy in the German test, German participants are at chance-level for the Dutch tests. These differences might be attributed to which languages are taught in school in these two countries. Generally, we can expect an increased accuracy for languages that are part of the curriculum in the respective countries.

GPT-3 generally performs well on most languages, indicating that large language models can outperform humans on the task and bots can potentially mimic human behavior. One way to avoid this would be to present text as images, continue preventing copying, limit viewing and response times, and possibly add modifications similar to CAPTCHA \cite{bursztein2010captcha}.

One important limitation of our method is it relied on linguistic resources (such as parsers, lemmatizers, and spell checkers). We found that these tools are available in at least 60 languages, suggesting that our method may apply to a large number of languages. For low resources languages where these are not available, in future work we plan to test if removing some of the processing steps would still allow good task performance. Alternatively, it is possible to recruit native speakers to test performance on candidate words and select only words that can be reliably detected by native speakers.

To conclude, WikiVocab offers a systematic and data-driven way to automatically produce language tests, which can be used to screen participants, particularly in the online environment. The generative nature of the test allows researchers to increase the diversity of  online research and we believe would be an important tool to assess language proficiency in human behavioral experiments. More broadly, our work shows how to scale experimental methods from the lab to the online world, thus contributing to the scale and diversity of cognitive science.

\bibliographystyle{apacite}

\setlength{\bibleftmargin}{.125in}
\setlength{\bibindent}{-\bibleftmargin}


\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.


\section{Introduction}
Medical imaging, like Computed Tomography\,(CT), requires lossless compression due to the essential content. However, lossless compression leads to high bit rates, which is challenging for fast transmitting and archiving the data. Especially modern multi-slice CT scanners, which allow for recording volumes of the human body over a specific time, can generate up to $200$\,MB/s~\cite{Ohnesorge2007}. The structure of such a dynamic 3-D+$t$ volume can be seen on the left side of Fig.\,\ref{fig:data}. Here, it consists of $T{=}4$ temporally equidistant spatial \mbox{3-D} volumes of size $X{\times}Y{\times}Z$ and can be splitted into $Z{=}3$ temporal sequences. Then, every temporal sequence can be encoded using state-of-the-art compression techniques.

Established video compression standards, like the High Efficiency Video Coding\,(HEVC) standard~\cite{2016}, have mainly been developed for consumer applications and are specialized for lossy compression. In~\cite{Lanz2018}, it has been shown that for lossless compression of dynamic CT data, 3-D subband coding\,(SBC)~\cite{Karlsson1988} performs better than HEVC. Since 3-D SBC is based on discrete wavelet transforms\,(WT), it inherently provides scalability features. This allows for transmitting a base layer\,(BL) with coarser quality and an enhancement layer\,(EL) on demand, comprising the residual data to reconstruct the original sequence. In Fig.\,\ref{fig:data}, a WT in temporal direction is performed, decomposing every single temporal sequence into two layers, both comprising only half the frame rate.
This is very advantageous for telemedicine applications, like remote surgery systems or transmission of medical content for remote assessment over wireless networks with limited channel capacity. Incorporating motion compensation\,(MC) methods into the lifting structure of a WT is called motion compensated temporal filtering\,(MCTF)~\cite{334985} and increases the quality of the BL significantly, but results also in a lower compression efficiency. To overcome this drawback, wavelet lifting with denoised update\,(WLDU)~\cite{Lanz2018} was recently introduced.

To further improve the compression efficiency of dynamic CT data, in this paper we propose to change the processing order of MC and denoising in the update step of WLDU. Besides this, we introduce a second denoising operation in the prediction step, which leads to additional coding gains.
After a short recap of compensated wavelet lifting in Section~\ref{sec:sbc}, the two proposed modifications will be explained in detail in Section~\ref{sec:prop}. Simulation results are given in Section~\ref{sec:results}, while Section~\ref{sec:conclusion} concludes this work.

\begin{figure}
	\includegraphics[width=\columnwidth]{dynamicVolume}
	\caption{Considered scenario to achieve temporal scalability of a 3-D+$t$ volume using a wavelet transform in temporal direction.}
	\pgfkeysvalueof{/tikz/dsp/label}}{fig:data}
	\vspace{-3ex}
\end{figure}

\section{Compensated Wavelet Lifting}
\pgfkeysvalueof{/tikz/dsp/label}}{sec:sbc}
Discrete wavelet transforms can be implemented very efficiently by the so-called lifting structure~\cite{Sweldens1995}, which consists of three steps. In Fig.\,\ref{fig:liftingMod}, a block diagram of the lifting structure for a WT in temporal direction is shown. In the first step, named splitting, the input signal is decomposed into even- and odd-indexed frames $f_{2t}$ and $f_{2t-1}$. In the next step, called prediction step, the even frames are predicted from the odd frames by a prediction operator $\mathcal{P}$. After subtracting the predicted values from $f_{2t}$, the highpass\,(HP) coefficients $\text{HP}_t$ are achieved. The last step is called update step and filters the HP coefficients by an update operator $\mathcal{U}$. Then, the HP coefficients are added to the odd frames, resulting in the lowpass\,(LP) coefficients $\text{LP}_t$. Since the resulting LP subband is very similar to the original signal, it directly serves as the BL. In contrast, the HP subband comprises the residual data and therefore corresponds to the EL.

\begin{figure}
	\centering
	\includegraphics[width=0.73\columnwidth]{liftingStructureMod}
	\caption{Block diagram of the motion compensated lifting structure containing the denoising filter (DN) of WLDU~\cite{Lanz2018}.}
	\pgfkeysvalueof{/tikz/dsp/label}}{fig:liftingMod}
	\vspace{-2ex}
\end{figure}

The lifting structure offers a very flexible framework and can be modified in multiple ways. To guarantee lossless reconstruction, rounding operators are applied to achieve integer-to-integer transforms~\cite{647983}. Additionally, MC methods can directly be incorporated into the lifting structure to avoid blurring and ghosting artifacts due to temporal displacements in the sequence. This is called MCTF and is achieved by realizing the prediction operator $\mathcal{P}$ by the warping operator $\mathcal{W}_{2t-1\rightarrow 2t}$. Accordingly, for the Haar wavelet filters, the HP coefficients are calculated by
\begin{equation}
	\text{HP}_t = f_{2t} - \lfloor \mathcal{W}_{2t-1\rightarrow 2t}(f_{2t-1})\rfloor.
\end{equation}
To achieve an equivalent WT, in the update step, the compensation has to be inverted. Consequently, the LP coefficients of the Haar transform are calculated by reversing the index of $\mathcal{W}$
\begin{equation}
	\text{LP}_t = f_{2t-1} + \lfloor \frac{1}{2} \mathcal{W}_{2t\rightarrow 2t-1}(\text{HP}_t)\rfloor.
\end{equation}
However, for noisy CT data, the required bit rate to encode the subbands will rise by applying MCTF. This is due to the fact that not only structural information is warped from one frame to the other, but also noise is warped back and forth.

To overcome the noise problem, prefiltering of the even and odd frames can be applied. However, lossless reconstruction is not provided anymore and therefore, such schemes are not suitable for professional applications, in which lossless compression is claimed.
One way to reduce this effect without compromising lossless reconstruction is to skip the update step entirely, which is called truncated WT. However, the resulting LP subband suffers not only from temporal aliasing, but loses also the property of being a good representative for both even- and odd-indexed frames. This makes the truncated WT useless for professional applications.
Another way to reduce the increase of the noise variance in the LP subband and to preserve lossless reconstruction was recently investigated in~\cite{Lanz2018} and is called wavelet lifting with denoised update\,(WLDU). The basic scheme can be seen in Fig.\,\ref{fig:liftingMod} highlighted in red. It has been shown that the application of suitable filters for denoising\,(DN) in the update step leads to a rising compression ratio of the LP subband, while the visual quality of the LP subband can be preserved for being used in telemedicine applications.

\section{Proposed Modifications}
\pgfkeysvalueof{/tikz/dsp/label}}{sec:prop}

\begin{figure}
	\centering
	\subfloat{\includegraphics[width=0.6\columnwidth]{liftingStructureWLDUr}}\\ \vspace{-1.5ex}
	\subfloat{\includegraphics[width=0.6\columnwidth]{liftingStructureWLDP}}\\ \vspace{-1.5ex}
	\subfloat{\includegraphics[width=0.6\columnwidth]{liftingStructureWLDPU}}
	\caption{Block diagrams of the different proposed setups. From top to bottom: WLDUr, WLDP, and WLDPU.}
	\pgfkeysvalueof{/tikz/dsp/label}}{fig:liftingMods}
	\vspace{-2ex}
\end{figure}

\begin{figure*}[t]
	\centering
	\resizebox{\textwidth}{!}{%
		\input{Figures/groupplot_BM3D_all.tikz}}
	\caption{$\text{PSNR}_{\text{LP}_t}$ (left) and $\text{SSIM}_{\text{LP}_t}$ (right) results, respectively, against the overall file size comparing WLDU, WLDUr, WLDP, and WLDPU. The arrow above the diagrams indicates the direction for increasing values of the filter strength $h{=}\xi{\cdot}\sigma_n^2$. For better presentation, only every eighth value for $h$ is plotted. Results are averaged over all 327 sequences.}
	\pgfkeysvalueof{/tikz/dsp/label}}{fig:val}
	\vspace{-2ex}
\end{figure*}

To further increase the compression ratio, we introduce two novel modifications in this paper, resulting in three different setups. First, we propose to change the processing order in the update step of WLDU as it is shown in Fig.\,\ref{fig:liftingMods} on the top. In the following, this setup will be called WLDU with reversed order\,(WLDUr).
However, by applying DN in the update step, only the compression ratio of the LP subband is affected. The HP subband is not touched so far. Therefore, the second setup targets the prediction step, as can be seen in Fig.\,\ref{fig:liftingMods} in the middle, called wavelet lifting with denoised prediction\,(WLDP).
The third setup is given by the combination of WLDUr and WLDP, called wavelet lifting with denoised prediction and update\,(WLDPU) and can be seen in Fig.\,\ref{fig:liftingMods} on the bottom.

\subsection{WLDU with reversed order}
\pgfkeysvalueof{/tikz/dsp/label}}{subsec:order}
As already described above, MCTF reduces the appearance of ghosting artifacts. However, it also increases the noise variance $\sigma_n^2$ of the single subbands, leading to a higher overall entropy~\cite{lnt2012-10} and therefore to a higher number of bits required for coding the subbands.
To compensate this drawback, WLDU was introduced. By denoising the HP frames before warping them back to the time step of the odd frames, the undesired increase of $\sigma_n^2$ in the HP subband caused by the MC in the prediction step was reduced.
However, by inversely compensating the HP frames and adding them to the odd frames in the update step, the noise variance of the LP subband can be increased further due to the inverse compensation process. Therefore, we propose to warp the HP frames at first to the time step of the odd frames and to denoise them by a suitable filter afterwards. By applying WLDUr, erroneously warped objects and noisy structures from both warping operations can be suppressed more efficiently, leading to a better compression ratio of the LP frames.

\subsection{Wavelet Lifting with Denoised Prediction}
\pgfkeysvalueof{/tikz/dsp/label}}{subsec:wldp}
The second proposed modification concerns the prediction step. So far, the DN filter of WLDU and WLDUr has only influenced the compression ratio of the resulting LP subband. However, for lossless reconstruction of the original sequence, the file size of the corresponding HP subband is also of interest. Therefore, we introduce a second DN filter in the prediction step similar to~\cite{faucris.122127984}, where in-loop denoising is performed on reference frames in the context of HEVC.
According to them, in-loop prediction filters can be applied on the reference frames themselves or on the motion compensated frames. However, by performing DN on the reference frames, for every DN setup, including the specific filter type and the chosen filter strength, new motion vectors have to be calculated.
This makes the whole framework very inflexible. Therefore, we decided to perform DN on the motion compensated frames.

\subsection{Wavelet Lifting with Denoised Prediction and Update}
\pgfkeysvalueof{/tikz/dsp/label}}{subsec:wldpu}
Since WLDP influences also the performance of the denoising in the update step, the combination of both modifications is also considered. By forcing noise reduction in both lifting steps, we expect a further increase of the compression ratio.

\section{Simulation Results}
\pgfkeysvalueof{/tikz/dsp/label}}{sec:results}

\begin{table*}[t]
	\begin{minipage}{0.74\textwidth}
		\caption*{TABLE I: Overall file size and $\text{PSNR}_{\text{LP}_t}$ for certain values of Fig.~\ref{fig:val} using BM3D filtering. Absolute and relative differences against MCTF are given. Negative numbers denote a smaller file size and a lower $\text{PSNR}_{\text{LP}_t}$ compared to MCTF. Results are averaged over all 327 sequences.}
		\setlength{\tabcolsep}{2pt}
		\resizebox{\textwidth}{!}{%
			\begin{tabular}{l|c|c|c|c||c|c}
				Method & Noise & File size & \multicolumn{2}{c||}{$\Delta$ to MCTF}   & $\text{PSNR}_{\text{LP}_t}$ {[}dB{]} & $\Delta$ to MCTF {[}dB{]} \\ \hhline{~|~|~|--||~|~}
				& parameter & LP + HP + MV {[}MB{]} & Absolute {[}MB{]} & Relative {[}$\%${]} &   &     \\ \hline
				MCTF                       &             & 1.2821 & -       & -       & 47.7322 & -     \\
				$\text{WLDU}$  & $\xi{=}16$  & 1.2760 & -0.0061 & -0.4758 & 47.5859 & -0.1463 \\
				\textcolor{mycolorblue}{$\text{WLDUr}$} & $\xi{=}\hphantom{1}8$  & 1.2573 & -0.0248 & -1.9343 & 47.5861 & -0.1461 \\
				\textcolor{mycolorgreen}{$\text{WLDP}$}  & $\xi{=}24$ & 1.2429 & -0.0392 & -3.0575 & 47.6041 & -0.1281 \\
				\textcolor{mycolorpurple}{$\text{WLDPU}$} & $\xi{=}\hphantom{1}8$  & 1.2258 & -0.0563 & -4.3912 & 47.5885 & -0.1437\\
				Truncated WT               &            & 1.2567 & -0.0254 & -1.9811 & 44.8806 & -2.8516
			\end{tabular}}
			\pgfkeysvalueof{/tikz/dsp/label}}{tab:results}
		\end{minipage}\hspace{0.5cm}
		\begin{minipage}{0.22\textwidth}
			\centering
			\input{Figures/groupplot_variances.tikz
			\vspace{-1ex}
			\caption*{{Fig.\,5: Range of relative bit rate savings for all proposed setups compared to MCTF.}}
			\pgfkeysvalueof{/tikz/dsp/label}}{fig:variances}
		\end{minipage}
		\vspace{-2ex}
	\end{table*}
	\setcounter{figure}{5}

	All experiments are evaluated on two 3-D+t medical CT data sets\footnote{The CT volume data sets were kindly provided by Siemens Healthineers.}, both describing a beating heart and comprising $10$ time steps and a spatial resolution of $512\times512$ pixels at a bit depth of $12$ bits per sample. The data sets consist of $127$ and $200$ slices per time step, respectively, summing up to $327$ single temporal sequences.
	The MC is realized by a mesh-based approach~\cite{8066388} using a grid size of $8$ pixels. The resulting LP and HP subbands are encoded without any loss using the wavelet-based volume coder JPEG\,2000~\cite{ITU-T-800}. Here, the \mbox{OpenJPEG}~\cite{openjpeg} implementation with four spatial decomposition steps in $xy$-direction is applied. The arising motion vectors\,(MV) are encoded using the QccPack library~\cite{fowler2000qccpack}. To evaluate the quality of the LP subband, we use the $\text{PSNR}_{\text{LP}_t}$ metric as well as the modified structural similarity index $\text{SSIM}_{\text{LP}_t}$, which were both introduced in~\cite{Lanz2018}. Both consider not only the similarity of the LP frames to the odd-indexed frames, but also the similarity of the LP frames to the even-indexed frames.
	All results are averaged over all $327$ sequences.

	To compare our novel modifications to WLDU, we apply Block Matching and 3-D Filtering\,(BM3D)~\cite{Dabov2007} for denoising, which was already used in~\cite{Lanz2018} and has given the highest compression gains.
	The corresponding filter parameters are also chosen identically. The denoising is evaluated for different values of the filter strength $h$. The filter strength $h$ is an input parameter of BM3D and is calculated by
	\begin{equation}
		h = \xi \cdot \sigma_n^2,
	\end{equation}
	where $\sigma_n^2$ describes the noise variance of the input image and $\xi$ denotes an arbitrary noise parameter which controls the strength of denoising in order to improve the compression efficiency. A higher value of $\xi$ leads to a lower noise variance of the output image and according to~\cite{lnt2012-10} to a better compression ratio. We examine the influence of $h$ by varying $\xi$ in a range of integer values from $1$ to $100$. $\sigma_n^2$ is estimated by a low-complexity algorithm proposed in~\cite{IMMERKAER1996300} at both encoder and decoder side.

	\subsection{Comparing WLDUr, WLDP, and WLDPU}

	In Fig.~\ref{fig:val}, all considered setups are evaluated for the BM3D filter. For better displaying, only every eighth value of $h$ is plotted. Additionally to the $\text{PSNR}_{\text{LP}_t}$ and rate results on the left side, the $\text{SSIM}_{\text{LP}_t}$ results are provided on the right side. The single curves for both metrices show the same behavior. Therefore, all further analyses are evaluated only by $\text{PSNR}_{\text{LP}_t}$.

	As mentioned in Section~\ref{sec:sbc}, the diagram shows that the truncated WT results in a very low $\text{PSNR}_{\text{LP}_t}$ value. However, the compression efficiency of the truncated WT is higher than for WLDU with a very strong filter strength of $h{=}100{\cdot}\sigma_n^2$.
	Further, it can be observed that all proposed setups outperform the compression efficiency of WLDU. By applying WLDUr, the quality of the LP subband is kept at a very high level until a filter strength of $h{=}8{\cdot}\sigma_n^2$, afterwards $\text{PSNR}_{\text{LP}_t}$ drops very fast. However, the achieved compression of BM3D for this filter strength gives similar results as the truncated WT.
	By using WLDP, further bit rate savings compared to WLDUr can be observed. Until a filter strength of $h{=}24\cdot\sigma_n^2$ the rate decreases significantly, while the quality in terms of $\text{PSNR}_{\text{LP}_t}$ is nearly kept constant. Afterwards the rate starts increasing and the $\text{PSNR}_{\text{LP}_t}$ value decreases simultaneously.
	The combination of both proposed modifications, shown by the violet curve for WLDPU, merges the properties of WLDUr and WLDP. This means that additional coding gains are achieved due to the influence of the DN operations to both subbands.

	For a fair comparison of all proposed setups, we draw a horizontal line at a constant level of $\text{PSNR}_{\text{LP}_t}$ in Fig.~\ref{fig:val} and choose those filter strengths which operate at this level. These values are marked with black circles in Fig.~\ref{fig:val} and are listed in Tab.~I. The differences against MCTF are also provided.
	Accordingly, with $\text{WLDU}$ and a filter strength of $h{=}16{\cdot}\sigma_n^2$ bit rate savings of only $0.4758\%$ are achievable. For $\text{WLDUr}$ with a filter strength of $h{=}8{\cdot}\sigma_n^2$, coding gains of $1.9343\%$ compared to MCTF are possible. Considering $\text{WLDP}$ with a filter strength of $h{=}24{\cdot}\sigma_n^2$, bit rate savings of $3.0575\%$ compared to MCTF can be achieved. Finally, $\text{WLDPU}$ achieves coding gains of $4.3912\%$ compared to MCTF for a filter strength of $h{=}8{\cdot}\sigma_n^2$.
	To confirm that these averaged results are representative for all 327 sequences, we additionally determine the minimum and maximum achieved compression results for the specific noise parameters from Tab.~I. for all proposed setups. The results can be seen in Fig.~5, where the whole range of achieved bit rate savings per sequence is depicted.

	In Fig.~\ref{fig:visual_results}, some visual results are provided. From left to right, the reference frame $f_1$, the current frame $f_2$, $\text{LP}_1$ for MCTF, and $\text{LP}_1$ for WLDPU from data set $1$ for slice position $z{=}1$ can be seen.
	\begin{figure*}
		\centering
		\begin{tikzpicture}
		\node[align=center] (ref) at (0,0) {\includegraphics[width=0.25\textwidth]{Figures/ref}};
		\node[align=center](ref_label) at (0,-1.8) {Reference frame $f_1$};
		\draw[draw=red]  (0.1,0.75) rectangle (0.1+0.4,0.75+0.3);
		\draw[draw=red]  (-0.7,-0.55) rectangle (-0.7+0.4,-0.55+0.3);
		\draw[draw=red]  (-0.6,0.55) rectangle (-0.6+0.4,0.55+0.4);
		\node[align=center] (cur) at (4.5,0) {\includegraphics[width=0.25\textwidth]{Figures/cur}};
		\node[align=center](cur_label) at (4.5,-1.8) {Current frame $f_2$};
		\draw[draw=red]  (4.6,0.75) rectangle (4.6+0.4,0.75+0.3);
		\draw[draw=red]  (3.8,-0.55) rectangle (3.8+0.4,-0.55+0.3);
		\draw[draw=red]  (3.9,0.55) rectangle (3.9+0.4,0.55+0.4);
		\node[align=center] (mctf) at (9,0) {\includegraphics[width=0.25\textwidth]{Figures/LP_mctf}};
		\node[align=center](mctf_label) at (9,-1.8) {$\text{LP}_1$ from MCTF};
		\node[fill=white,position=-54:{-12mm} from mctf] (n1) {\small{$49.20$\,dB}};
		\node[fill=white,position=234:{-12mm} from mctf] (n2) {\small{$303.21$\,kB}};
		\node[align=center] (wldpu) at (13.5,0) {\includegraphics[width=0.25\textwidth]{Figures/LP_wldpu}};
		\node[align=center](wldpu_label) at (13.5,-1.8) {$\text{LP}_1$ from WLDPU};
		\node[fill=white,position=-54:{-12mm} from wldpu] (n3) {\small{$49.10$\,dB}};
		\node[fill=white,position=234:{-12mm} from wldpu] (n4) {\small{$292.84$\,kB}};
		\end{tikzpicture}
		\caption{Visual results comparing LP frames from WLDPU and MCTF against original frames. Dominant moving areas are highlighted. $\text{PSNR}_{\text{LP}_t}$ and file size to reconstruct $f_1$ and $f_2$, comprising $\text{LP}_1$, $\text{HP}_1$, and the corresponding MVs, are provided.}
		\pgfkeysvalueof{/tikz/dsp/label}}{fig:visual_results}
		\vspace{-3ex}
	\end{figure*}
	In the original frames, some dominant moving areas are marked by red boxes. As can be seen, both LP frames provide very similar visual results, which is also confirmed by the $\text{PSNR}_{\text{LP}_t}$ values. Additionally, the file size to reconstruct $f_1$ and $f_2$ at decoder side, including the file size of $\text{LP}_1$, $\text{HP}_1$, and the corresponding MVs, is given. This rate is reduced by $3.42\%$ from MCTF to WLDPU.

	\subsection{Comparing different filters}

	\begin{figure}[t]
		\centering
		\resizebox{0.47\textwidth}{!}{%
			\input{Figures/groupplot_WLDUr.tikz}}
		\caption{Comparing different filters for WLDUr, WLDP, and WLDPU. Results are averaged over all 327 sequences.
		}
		\pgfkeysvalueof{/tikz/dsp/label}}{fig:WLDUr}
		\vspace{-3ex}
	\end{figure}

	As a second experiment, different filters for each setup are evaluated. In addition to BM3D filtering, the Adaptive Wiener Filter\,(AWF)~\cite{1990ph...book.....L}, the Non-Local Means algorithm\,(NLM)~\cite{1467423}, and Guided Image Filtering\,(GIF)~\cite{He2010} is applied.
	These filters were already applied in~\cite{Lanz2018} to evaluate WLDU.
	The upper diagram of Fig.~\ref{fig:WLDUr} shows the performance of all filters and filter strengths for WLDUr. For AWF and GIF, the quality of the LP subband is kept at a very high level for the whole range of tested filter strengths, while the overall file size is reduced significantly. In contrast, the NLM and BM3D filters show this behavior only until a filter strength of $h{=}8{\cdot}\sigma_n^2$, afterwards the quality drops very fast and the rate increases. However, the achieved compression ratio of BM3D for this filter strength is outstanding compared to all other filter setups.
	In the middle diagram of Fig.~\ref{fig:WLDUr}, the performance of all filters and filter strengths for WLDP is given. The same behavior for all filters like for WLDUr can be observed, but in a larger stretch. This means, further bit rate savings compared to WLDUr are achieved at approximately the same level of $\text{PSNR}_{\text{LP}_t}$. Again, the BM3D filter performs best.
	The performance of all filters and filter strengths for WLDPU can be observed in the lower diagram of Fig.~\ref{fig:WLDUr}. It shows that all filters profit by the combination of WLDP and WLDUr. Consistently, BM3D filtering achieves the best results.

	Obviously, the right choice of the filter strength $h$ is very important for achieving a high compression ratio and keeping the quality of the LP subband in terms of $\text{PSNR}_{\text{LP}_t}$ at a high level.
	Therefore, the next step of our work targets the generation of a suitable signal model, which allows for a theoretical analysis of the noise impact to the entire framework as well as an optimum choice of the filter strength by minimizing $\sigma_n^2$ in both subbands.
	Nevertheless, we showed that it is possible to outperform state-of-the-art methods like WLDU and the truncated WT with our proposed modifications. This is very advantageous for professional applications, where a scalable representation always requires a suitable trade-off regarding the quality of the LP subband and the corresponding compression ratio.

	\section{Conclusion}
	\pgfkeysvalueof{/tikz/dsp/label}}{sec:conclusion}
	A scalable representation of video data is often required in professional applications. This can efficiently be reached by wavelet-based approaches. By incorporating MC methods into the lifting structure, the quality of the LP subband, which serves as a temporally downscaled representative of the original sequence, can be increased significantly. However, the compression efficiency of dynamic CT data suffers enormously from MCTF. By integrating denoising filters into the lifting structure, this drawback can be lessened. In this work, a novel processing order for an already existing approach of denoising in the update step was introduced. Further, a second denoising filter in the prediction step was proposed. Both modifications offer a significant enhancement of the compression efficiency for dynamic CT data, while the quality of the LP subband is kept nearly constant for a certain range of filter strengths. Additionally, the property of lossless reconstruction is still given, which makes our novel approach very attractive for high-fidelity application scenarios, where the channel capacity is limited.

	\bibliographystyle{IEEEbib}