\section{}
Artificial neural networks (ANN) have become the mainstream acoustic modeling technique for large vocabulary automatic speech recognition (ASR). A conventional ANN features a multi-layer architecture that requires massive amounts of computation. The brain-inspired spiking neural networks (SNN) closely mimic the biological neural networks and can operate on low-power neuromorphic hardware with spike-based computation. Motivated by their unprecedented energy-efficiency and rapid information processing capability, we explore the use of SNNs for speech recognition. In this work, we use SNNs for acoustic modeling and evaluate their performance on several large vocabulary recognition scenarios. The experimental results demonstrate competitive ASR accuracies to their ANN counterparts, while require significantly reduced computational cost and inference time. Integrating the algorithmic power of deep SNNs with energy-efficient neuromorphic hardware, therefore, offer an attractive solution for ASR applications running locally on mobile and embedded devices.

\tiny
\fontsize{8}{11}\helveticabold {\section{Keywords:} Deep Spiking Neural Networks, Automatic Speech Recognition, Tandem Learning, Neuromorphic Computing, Acoustic Modeling}
\end{abstract}

\section{Introduction}
\label{sec:intro}
Automatic speech recognition (ASR) has enabled the voice interface of mobile devices and smart home appliances in our everyday life. The rapid progress in the integration of voice interfaces has been viable on account of the remarkable performance of the ASR systems using artificial neural networks (ANN) for acoustic modeling~\citep{lippmann1989,lang1990,hinton2012deep,YuDeng15}. Various ANN architectures, either feedforward or recurrent, have been investigated for modeling the acoustic information preserved in speech signals~\citep{dahl2012,graves2013,abdelhamid2014}.

The performance gains come with immense computational requirements often due to the time-synchronous processing of input audio signals. Several techniques have been proposed to reduce the computational load and memory storage of ANNs by reducing the number of parameters that have to be used for inference \citep{sainath2013,xue2013,he2014,povey2018}. Another common solution, for reducing the amount of processed speech, uses a wake word or phrase to initialize the embedded ASR engine and starts listening to input speech~\citep{zehetner2014,Sainath2015ConvolutionalNN,wu2018wakeword}. Moreover, most devices with voice control rely on cloud-based ASR engines rather than local on-device solutions. The necessity of online processing of speech via cloud computing comes with various concerns, such as data security and processing speed, etc. There have been multiple efforts to develop on-device ASR solutions in which the speech signal is processed locally using the computational resources of mobile devices~\citep{Lei2013AccurateAC,mcgraw2016}.

Alternatively, event-driven models such as spiking neural networks (SNNs) inspired by the human brain have attracted ever-growing attention in recent years. The human brain is remarkably efficient and capable of performing complex perceptual and cognitive tasks. Notably, the adult's brain only consumes about 20 watts to solve complex tasks that are equivalent to the power consumption of a dim light bulb \citep{laughlin2003communication}. While brain-inspired ANNs have demonstrated great capabilities in many perceptual \citep{he2016deep, xiong2017toward} and cognitive tasks \citep{silver2017mastering}, these models are computationally intensive and memory inefficient to operate as compared to the biological brains. Unlike ANNs, asynchronous and event-driven information processing of SNNs resembles the computing paradigm that observed in the human brains, whereby the energy consumption matches the activity levels of  sensory stimuli. Given temporally sparse information transmitted in the surrounding environment, the event-driven computation, therefore, exhibits great computational efficiency than the synchronous computation used in ANNs.

Neuromorphic computing (NC), as a non-von Neumann computing paradigm, mimics the event-driven computation of the biological neural systems with SNN in silicon. The emerging neuromorphic computing architectures \citep{davies2018loihi, merolla2014million, furber2012overview} leverage on the massively parallel, low-power computing units to support spike-based information processing. Notably, the design of co-located memory and computing units effectively circumvents the von Neumann bottleneck of low-bandwidth between memory and the processing units \citep{monroe2014neuromorphic}. Therefore, integrating the algorithmic power of deep SNNs with the compelling energy efficiency of NC hardware represents an intriguing solution for pervasive machine learning tasks and always-on applications. Furthermore, growing research efforts are devoted to developing novel non-volatile memory devices for ultra-low-power implementation of biological synapses and neurons \citep{tang2019bridging}.

Some preliminary work on SNN-based phone classification or small-vocabulary speech recognition systems have been explored in \citep{liaw1998,nager2002,loiselle2005,holmberg2005,kroger2009,ref72,ref73,wu2018ijcnn,wu2018spiking,zhang2015digital,zhang2019mpd,bellec2018long}. However, these SNN-based ASR systems are far from the scale and complexity of  modern commercialized ANN-based ASR systems. It is mainly due to lacking effective training algorithms for deep SNNs and efficient software toolbox for SNN-based ASR systems.

Due to the discrete and non-differentiable nature of spike generation, the powerful error back-propagation algorithm is not directly applicable to the training of deep SNNs. Recently, considerable research efforts are devoted to addressing this problem and the resulting learning rules can be broadly categorized into the SNN-to-ANN conversion \citep{cao2015spiking,diehl2015fast}, back-propagation through time with surrogate gradient  \citep{neftci2019surrogate,wu2018direct} and tandem learning \citep{hybridlearning}. Despite several successful attempts on the large-scale image classification tasks with deep SNNs \citep{ethImageNet,hu2018spiking,sengupta2019going,hybridlearning}, their applications to the large-vocabulary continuous ASR (LVCSR) tasks remain unexplored. In this work, we explore an SNN-based acoustic model for LVCSR using a recently proposed tandem learning rule \citep{hybridlearning} that supports an efficient and rapid inference. The experiments on several large vocabulary ASR benchmarks have shown competitive classification accuracies over the baseline ANN models. To the best of our knowledge, this is the first work that successfully applied SNNs for LVCSR tasks. Our preliminary study of an SNN-based acoustic model has also revealed compelling prospects of rapid inference and unprecedented energy efficiency of a neuromorphic approach.

The rest of the paper is organized as follows: In Section \ref{sec:relatedworks}, we first give an overview of spiking neural networks, large vocabulary ASR systems, and existing SNN-based ASR systems. In Section \ref{sec:method}, we introduce the spiking neuron model and the neural coding scheme that converts acoustic features into spike-based representation. We further present a recently introduced tandem learning framework for SNN training and how it is used to train deep SNN-based acoustic models. In Section~\ref{sec:result}, we present experimental results on the learning capability and energy efficiency of SNN-based acoustic models across three different types of recognition tasks including phone recognition, low-resourced and standard large-vocabulary ASR, and compare those to the ANN-based implementations. Finally, a discussion on the experimental findings is given in Section~\ref{sec:dis}.

\section{Fundamentals and Related Work}
\label{sec:relatedworks}
\subsection{Spiking Neural Networks}
The third generation spiking neural networks are originally studied as models to describe the information processing in the biological neural networks, wherein the information is communicated and exchanged via stereotypical action potentials or spikes \citep{gerstner2002spiking}. Neuroscience studies reveal that the temporal structure and frequency of these spike trains are both important information carriers in the biological neural networks. As will be introduced in Section \ref{subsec:neuronmodel}, the spiking neuron operates asynchronously and integrates the synaptic current from its incoming spike trains. An output spike is generated from the spiking neuron whenever its membrane potential crosses the firing threshold, and this output spike will be propagated to the connected neurons via the axon.

Motivated by the same connectionism principle, SNNs share the same network architectures, either feedforward or recurrent, with the conventional ANNs that use analog neurons. As shown in Figure~\ref{fig2}, the early classification decision can be made from the SNN since the generation of the first output spike. However, the quality of the classification decision is typically improved over time with more evidence accumulated. It differs significantly from the synchronous information processing of the conventional ANNs, where the output layer needs to wait until all preceding layers are fully updated. Therefore, despite information is transmitted and processed at a speed that is several orders of magnitude slower in neural substrates than signal processing in modern transistors, biological neural systems can perform complex tasks rapidly. For more overviews about SNNs and their applications, we refer readers to  \citep{pfeiffer2018deep, tavanaei2018deep}.

\subsection{Large Vocabulary Automatic Speech Recognition}
\label{ssec:lvasr}
As shown in Figure~\ref{asrblock}, conventional ASR systems uses acoustic and linguistic information preserved in three distinct components to convert speech signals to the corresponding text: (1) an acoustic model for preserving the statistical representations of different speech units, e.g. phones, from speech features, (2) a language model for assigning probabilities to the co-occurring word sequences and (3) a pronunciation lexicon for mapping the phonetic transcriptions to orthography. These resources are jointly used to determine the most likely hypothesis in the decoding stage.

The acoustic modeling has been achieved using various statistical models such as Gaussian Mixture Models (GMM) for assigning frame-level phone posteriors in conjunction with a Hidden Markov Model (HMM) for duration modeling~\citep{YuDeng15}. More recently, ANN-based approaches have become the standard acoustic models providing state-of-the-art performance across a wide spectrum of ASR tasks~\citep{hinton2012deep}. Together with the numerous ANN architectures explored for acoustic modeling, several end-to-end ANN architectures have been proposed for directly mapping speech features to text with optional use of the other linguistic components~\citep{graves2014,chan2016,watanabe2017}.

The probabilistic definition of acoustic modeling becomes more evident via the Bayesian formulation of the speech recognition task. Given a target speech signal that segmented into T overlapped frames, the resulting frame-wise features can be represented as $\textbf{O} = [\textbf{o}_\textsubscript{1}, \textbf{o}_\textsubscript{2},..., \textbf{o}_\textsubscript{T}]$. An ASR system assigns the probability $P(\textbf{W}|\textbf{O})$ to all possible word sequences $\textbf{W} = [w_\textsubscript{1}, w_\textsubscript{2},...]$, and the word sequence $\hat{\textbf{W}}$ with the highest probability is the recognized output,
\begin{equation}
 \hat{\textbf{W}} = \argmax_{\textbf{W}} P(\textbf{W}|\textbf{O})
\end{equation}
The probability $P(\textbf{W}|\textbf{O})$ can be decomposed into two parts by applying the Bayes' rule as below,
\begin{equation}
 \hat{\textbf{W}} = \argmax_{\textbf{W}} \frac{P(\textbf{O}|\textbf{W})P(\textbf{W})}{P(\textbf{O})}
\end{equation}
$P(\textbf{O})$ can be omitted as it does not depend on $\textbf{W}$. This results in
\begin{equation}
\label{eq:search}
\hat{\textbf{W}} = \argmax_{\textbf{W}} P(\textbf{O}|\textbf{W})P(\textbf{W})
\end{equation}
which formally defines the theoretical foundation that are grounded in conventional ASR systems. $P(\textbf{W})$ is the prior probability of the word sequence $\textbf{W}$ and this probability is provided by the language model which is trained on a large written corpus of the target language. $P(\textbf{O}|\textbf{W})$ is the \textit{likelihood} of the observed feature sequence $\textbf{O}$ given the word sequence $\textbf{W}$, and this probability is associated with the acoustic model. The acoustic model captures the information about the acoustic component of speech signals, aiming to classify different acoustic units accurately. Traditionally, each phone in the phonetic alphabet is modeled using multiple three-state HMM models for different preceding and following phonetic context (triphone)~\citep{lee1990}. The emission probability of these HMM states are shared (tied) among different models to reduce the number of model parameters~\citep{hwang1993}. The output layer of the ANN-based acoustic model is designed accordingly and trained to assign these frame-level tied triphone HMM state (senone) probabilities~\citep{dahl2012}. The output layer uses the softmax function to normalize the output into a probability distribution. These values are scaled with the prior probabilities of each class, obtained from the training data, to determine the likelihood values. These likelihood values are later combined with the probabilities assigned by the language model during the decoding stage so as to find the most likely hypothesis.

Speech features, used as inputs to the acoustic model, describe the spectrotemporal dynamics of the speech signal and discriminate among different phones in the target language. Mel-frequency cepstral coefficients(MFCC)~\citep{davis1980} features are commonly used in conjunction with the GMM-HMM acoustic model. The MFCC features are extracted by (1) performing short-time Fourier transform, (2) applying triangular Mel-scaled filter banks to calculate the power at each Mel frequency in log domain (FBANK) and (3) performing a discrete cosine transform to decorrelate the FBANK features. The third step is often skipped and FBANK features are often used when training ANN-based acoustic models since these models can handle correlation among features. In this work, we incorporate deep SNNs for acoustic modeling instead of the conventional ANNs and compare their ASR performance in different ASR scenarios including phone recognition, low-resourced and standard large vocabulary ASR. The ASR performance obtained using popular speech features have been reported to explore the impact of the feature representation space and its dimensionality for SNN-based acoustic models.

\subsection{Speech Recognition with Spiking Neural Network}
SNNs are well-suited for representing and processing spatial-temporal signals, they hence possess great potentials for speech recognition tasks. Tavanaei et al. \citep{ref72, ref73} proposed SNN-based feature extractors to extract discriminative features from the raw speech signal using unsupervised spiking-timing-dependent plasticity (STDP) rule. While connecting these SNN-based feature extractors with Support Vector Machine (SVM) or Hidden Markov Model (HMM) classifiers, competitive classification accuracies were demonstrated on the isolated spoken digit recognition task. Wu et al. \citep{wu2018ijcnn,wu2018spiking} introduced a SOM-SNN framework for environmental sound and speech recognition. In this framework, the biological-inspired self-organizing map (SOM) is utilized for feature representation, which maps frame-based acoustic features into a spike-based representation that is both sparse and discriminative. The temporal dynamic of the speech signal is further handled by the SNN classifier. Zhang et. al \citep{zhang2019mpd} presented a fully SNN-based speech recognition framework, wherein the spectral information of consecutive frames are encoded with threshold coding and subsequently classified by the SNN that is trained with a novel membrane potential-driven aggregate-labeling learning algorithm.

Recurrent network of spiking neurons (RSNNs) exhibit greater memory capacity than the aforementioned feedforward frameworks. They can capture long temporal information that are useful for speech recognition tasks. In \citep{zhang2015digital}, Zhang et al. presented a spiking liquid-state machine (LSM) speech recognition framework which is attractive for low-power very-large-scale-integration (VLSI) implementation. Bellec et al. recently demonstrated state-of-the-art phone recognition accuracy on the TIMIT dataset by adding neuronal adaptation mechanism to the vanilla RSNNs \citep{bellec2018long}. It is the first time that RSNNs approaching the performance of LSTM networks \citep{greff2016lstm} on the speech recognition task. These preliminary works on the SNN-based ASR systems are however limited to the phone classification or small vocabulary isolated spoken digit recognition tasks. In this work, we apply deep SNNs to LVCSR tasks and demonstrate competitive accuracies over the comparable ANN-based ASR systems.

\section{Methods}
\label{sec:method}
\subsection{Spiking Neuron Model}
\label{subsec:neuronmodel}
As shown in Figure \ref{learningRule}, the frame-based features are extracted and input into the SNN-based acoustic models. Given the short temporal duration of segmented frames and the slow variation of speech signals, these features are typically assumed to be stationary over the short time-period of segmented frames. In this work, we use the integrate-and-fire (IF) neuron model with reset by subtraction scheme \citep{ethImageNet}, which can effectively process these stationary frame-based features with minimal computational costs. At each time step $t$ of a discrete-time simulation, the incoming spikes to neuron $j$ at layer $l$ are transduced into synaptic current as follows
\begin{equation}
z_j^l(t) = \sum\nolimits_i {w_{ji}^{l - 1} \cdot \theta _i^{l - 1}} (t) + b_j^l \slash T_{e}\\
\label{eq1}
\end{equation}
where $\theta _i^{l - 1}(t)$ indicates the occurrence of an input spike from afferent neuron $i$ at time step $t$. In addition, the $w_{ji}^{l - 1}$ denotes the synaptic weight that connects presynaptic neuron $i$ from layer $l-1$. Here, $b_j^l/T_{e}$ can be interpreted as a constant injecting current across the encoding time window of size $T_{e}$, and $b_j$ is determined from the bias term of the coupled analog neurons which will be explained in the tandem learning section. As shown in Figure \ref{ifNeuron}, neuron $j$ integrates the input current $z_j^l(t)$ into its membrane potential $V_j^l(t)$ as per Eq. \ref{eq2}. The $V_j^l(0)$ is reset and initialized to zero for every new frame-based feature input. Without loss of generality, a unitary membrane resistance is assumed here. An output spike is generated whenever $V_j^l(t)$ crosses the firing threshold $\vartheta$ (Eq. \ref{eq4}), which we set to a value of $1$ for all the experiments by assuming that all synaptic weights are normalized with respect to the $\vartheta$.
\begin{equation}
V_j^l(t) = V_j^l(t - 1) + z_j^l(t) - \vartheta \cdot \theta _j^l(t - 1)
\label{eq2}
\end{equation}
\begin{equation}
V_j^l(0) =  0
\label{eq3}
\end{equation}
\begin{equation}
\theta _j^l(t) = \Theta (V_j^l(t) - \vartheta) \;\; with \;\; \Theta (x) = \left\{\begin{array}{l}1,\;\;\; if\;x \ge 0\\
0,\;\;\;otherwise\;
\end{array} \right.
\label{eq4}
\end{equation}
According to Eqs. \ref{eq1} and \ref{eq2}, the free aggregated membrane potential of neuron $j$ (no firing) in layer $l$ can be expressed as
\begin{equation}
V_j^{l,f}= \sum\nolimits_i {w_{ji}^{l - 1} \cdot c_i^{l - 1}} + b_j^l \\
\label{eq5}
\end{equation}
where $c_i^{l - 1}$ is the input spike count from pre-synaptic neuron $i$ at layer $l-1$ as per Eq. \ref{sc}.
\begin{equation}
c_i^{l-1} = \sum\nolimits_{t=1}^{T_e} {\theta _i^{l-1}(t)}. \\
\label{sc}
\end{equation}
The $V_j^{l,f}$ summarizes the aggregate membrane potential contributions of the incoming spikes from pre-synaptic neurons while ignoring their temporal distribution. As will be explained in the tandem learning framework section, this intermediate quantity links the SNN layers to the coupled ANN layers for parameter optimization.

\subsection{Neural Coding Scheme}
SNNs process information transmitted via spike trains, therefore, special mechanisms are required to encode the continuous-valued feature vectors into spike trains and decode the classification results from the activity of output neurons. To this end, we adopt the spiking neural encoding scheme that proposed in the tandem learning framework \citep{hybridlearning}. This encoding scheme first transforms frame-based input feature vector $X^0$ (e.g., MFCC or FBANK features), where ${X^0} = {[x_1^0,x_2^0, \cdot  \cdot  \cdot ,x_n^0]^T}$, through a weighted layer of rectified linear unit (ReLU) neurons as follows
\begin{equation}
V_j^{0,f}(0) \equiv a_j^{0} = f(\sum\nolimits_i {w_{ji}^{0} \cdot x_i^{0}} + b_j^0) \\
\label{encode_eq}
\end{equation}
where $w_{ji}^{0}$ is the strength of the synaptic connection between the input $x_i^{0}$ and ReLU neuron $j$. The $b_j^0$ is the corresponding bias term of the neuron $j$, and $f(\cdot)$ denotes the ReLU activation function. The free aggregate membrane potential $V_j^{0,f}(0)$ is defined to be equal to the activation value $a_j^{0}$ of the ReLU neuron $j$.  We distribute this quantity over the encoding time window $T_e$ and represent it with spike trains as per Eqs. \ref{encode_spike} and \ref{encode_spike2}.
\begin{equation}
\theta _j^0(t) = \Theta (V_j^{0,f}(t-1) - \vartheta)\;\;
\label{encode_spike}
\end{equation}
\begin{equation}
V_j^{0,f}(t) = V_j^{0,f}(t-1) - \vartheta \cdot \theta _j^0(t) \;\;
\label{encode_spike2}
\end{equation}
Altogether, the spike train $s^0$ and spike count $c^0$ that output from the neural encoding layer can be represented as follows
\begin{equation}
{s^0} = \{{\theta ^0}(1),...,{\theta ^0}({T_{e}})\}
\end{equation}
\begin{equation}
{c^0} = \sum\nolimits_{t = 1}^{{T_{e}}} {{\theta ^0}(t)}
\end{equation}
This encoding layer performs weighted transformation inside an end-to-end learning framework. It transforms the original input representation to match the size of the encoding time window $T_e$ and represents the transformed information via spike trains. This encoding scheme is beneficial for rapid inference since the input information can be effectively encoded within a short encoding window. Start from this neural encoding layer, as shown in Figure \ref{learningRule}, we input the spike count ${c^l} $ and ${s^l}$ to subsequent ANN and SNN layers for tandem learning.

To ensure smooth learning with high precision error gradients derived at the output layer, we use the free aggregate membrane potential of output spiking neurons for neural decoding. Considering that the dimensionality of input feature vectors and output classes are much smaller than that of hidden layers, the computation required will be limited when deploying these two layers onto the edge devices.

\subsection{Tandem Learning for Training Deep SNNs}
\label{sec:trainSNN}
Although IF neurons do not emulate rich temporal dynamics of biological neurons, they are however ideal for working with the neural representation that employed in this work, where spike timings play an insignificant role. It is worth noting that connections are commonly drawn between the activation value of ReLU neurons and the steady-state firing rate of IF neurons \citep{ethImageNet}. Here, we present a recently proposed SNN learning rule, under the tandem neural network configuration, that exploits such a connection between the activation value of ANN neurons and the spike count of IF neurons.

By neglecting the temporal dynamic of IF neuron that due to the temporal distribution of incoming spike trains, we may consider the $V_j^{l,f}$ as the main information carrier for SNN layers. The following one-to-one correspondence between the free aggregate membrane potential $V_j^{l,f}$ of spiking neurons and the pseudo `spike count' $a_j^l$ of artificial neurons can be established.
\begin{equation}
\begin{aligned}
a_j^l = min\left(\left\lfloor {V_j^{l,f}}\slash{\vartheta} \right\rfloor \cdot\left( {V_j^{l,f} > 0} \right), T_{e} \right)
\end{aligned}
\label{eq6}
\end{equation}
As shown in Figure \ref{learningRule}, during the activation forward propagation, the SNN layers are used to determine the exact spike representation which then propagate the aggregate spike counts and spike trains to the subsequent ANN and SNN layers, respectively. This interlaced layer structure ensures the information that forward propagated to the coupled ANN and SNN layers are synchronized. Taking Eq. \ref{eq6} as the activation function of ANN layers and using straight-through estimator \citep{bengio2013estimating} to address the discontinuity of the rounding operation, we can use the error gradients derived from ANN layers to approximate those of the coupled SNN layers. It worth noting that the ANN is just an auxiliary structure to facilitates the training of SNN, while only SNN is used during inference.

Notably, competitive classification accuracies are demonstrated with this tandem learning rule for the image classification on the ImageNet dataset \citep{hybridlearning}. By analyzing the relationships between the approximated `spike count' $a_j^l$ and the actual spike count $c_j^{l}$ in a high-dimensional space, Wu et al. have argued that the modified learning dynamic of such a decoupled network can approximate that of an intact ANN. The details of this tandem learning rule are provided in the Algorithm \ref{algo}.

\subsection{SNN-based Acoustic Modeling}
To train the deep SNN-based acoustic models, which is the main contribution of this work, several popular speech features have been extracted from the training recordings as described in Section~\ref{ssec:lvasr}. Before being fed into the SNNs, these input speech features are contextualized by splicing multiple frames so as to exploit more temporal context information. Before training the SNN-based acoustic model, alignments of the speech features with the target senone labels are obtained using a conventional GMM-HMM-based ASR system similar to that described in \citep{dahl2012}. These frame-level alignments enable the training of the deep SNN acoustic model with the tandem learning approach. During the training, the deep SNN learns to map input speech features to posterior probabilities of senones (cf. Section~\ref{ssec:lvasr}) by passing the input speech frames through multiple layers of spiking neurons.

During the inference phase, the acoustic scores provided by the trained SNN model are combined with the information stored in the language model and pronunciation lexicon. It is a common practice to use the weighted finite state transducers (WFST)~\citep{mohri2002} as a unified representation of different ASR resources for creating the search graph containing possible hypotheses. The main motivation for using the WFST-based decoding is: (1) the straightforward composition of different ASR resources for constructing a mapping from HMM states to word sequences and (2) the existence of efficient search algorithms operating on WFST that speed up the decoding process. As a result of the search process, the most likely hypotheses are found and stored in the form of a lattice. The ASR output is chosen based on the weighted sum of the acoustic and language model scores belonging to hypotheses in the lattice. For further details of the WFST-based decoding approach used in this work, we refer the reader to~\citep{povey2012}. In the following sections, we describe the ASR experiments conducted to evaluate the recognition performance of the proposed SNN-based acoustic modeling in several recognition scenarios.

\subsection{Training and Evaluation}
\label{sec:exp}
\subsubsection{Datasets}
The performance of the proposed SNN-based acoustic models is investigated in three different ASR tasks: (1) phone recognition using the TIMIT corpus~\citep{timit}, (2) low-resourced ASR task using the FAME code-switching Frisian-Dutch corpus~\citep{yilmaz2016} and (3) standard large-vocabulary continuous ASR task using the Librispeech corpus~\citep{librispeech}. All speech data used in the experiments has a sampling frequency of 16 kHz.

The train, development and test sets of the standard TIMIT corpus contain 3,696, 400 and 192 utterances from 462, 50 and 24 speakers, respectively. Each utterance is phonetically transcribed using a phonetic alphabet consisting of 48 phones in total. The training data of the FAME corpus comprises of 8.5 hours and 3 hours of broadcast speech from Frisian and Dutch speakers, respectively. The training utterances are spoken by 382 speakers in total. This bilingual dataset contains Frisian-only and Dutch-only utterances as well as mixed utterances with inter-sentential, intra-sentential and intra-word code-switching~\citep{myers1989}. The development and test sets consist of 1 hour of speech from Frisian speakers and 20 minutes of speech from Dutch speakers each. The total number of speakers is 61 in the development set and 54 in the test set.

The Librispeech corpus contains 1,000 hours of reading speech in total collected from audiobooks. This publicly available corpus\footnote{www.openslr.org/resources/12} has been considered as a popular benchmark for ASR algorithms with multiple training and testing settings. In the ASR experiments, we train acoustic models using the 100 (train\_clean\_100) and 360 (train\_clean\_360) hours of speech and apply these models to the clean development (dev\_clean) and test (test\_clean) sets. Further details about this corpus can be found in~\citep{librispeech}.

\subsubsection{Implementation Details}
All ASR experiments are performed using the PyTorch-Kaldi ASR toolkit~\citep{pytorch-kaldi}. This recently introduced toolkit inherits the flexibility of PyTorch toolkit~\citep{paszke2017automatic} for ANN-based acoustic model development and the efficiency of Kaldi ASR toolkit~\citep{povey2011kaldi}. We implement the SNN tandem learning rule in PyTorch and integrate it into the PyTorch-Kaldi toolkit for training the proposed SNN-based acoustic models (cf. Figure~\ref{learningRule}). The PyTorch implementation of the described SNN acoustic models will be made available online soon. For the baseline ANN models, the standard multi-layer perceptron recipes are used. The Kaldi toolkit is used for obtaining the initial alignments, feature extraction, graph creation, and decoding.

For all recognition scenarios, ANNs and SNNs are constructed with 5 hidden layers and 2048 hidden units each using the ReLU activation function. Each fully-connected layer is followed by a batch normalization layer and a dropout layer with a drop probability of 10\% to prevent overfitting. We train these models using various popular speech features including the 13-dimensional Mel-frequency cepstral coefficient (MFCC) feature, 23-dimensional Mel-filterbank (FBANK) feature, and higher resolution 40-dimensional MFCC and FBANK features. We further extract feature space maximum likelihood linear regression (FMLLR)~\citep{gales1998} features to explore the impact of speaker-dependent features. All features include the deltas and delta-deltas; mean and variance normalization are applied before the splicing. The time context size is set to 11 frames by concatenating 5 frames preceding and following. All features are encoded within a short time window of 10-time steps for SNN simulations.

The neural network training is performed by mini-batch Stochastic Gradient Descent (SGD) with an initial learning rate of 0.08 and a minibatch size of 128.~The learning rate is halved if the improvement is less than a preset threshold of 0.001.~The final acoustic models of the TIMIT and FAME corpora are obtained after 24 training epochs, while the models of the Librispeech corpus are trained for 12 epochs.

For the TIMIT and Librispeech ASR tasks, we follow the same language model (LM) and pronunciation lexicon preparation pipeline as provided in the corresponding Kaldi recipes\footnote{https://github.com/kaldi-asr/kaldi/tree/master/egs/\{timit,librispeech\}}. The smallest 3-gram LM (\texttt{tgsmall}) of the Librispeech corpus is used to create the graph for the decoding stage. The details of the LM and lexicon used in the FAME recognition task are given in~\citep{yilmaz2018_1}.

\subsubsection{Evaluation Metrics}
\paragraph{ASR Performance}
The phone recognition on the TIMIT corpus is reported in terms of the phone error rate (PER). The word recognition accuracies on the FAME and Librispeech corpora are reported in terms of word error rate (WER). Both metrics are calculated as the ratio of all recognition errors (insertion, deletion, and substitution) and the total number of phones or words in the reference transcriptions.

\paragraph{Energy Efficiency: Counting Synaptic Operations}
To compare the energy efficiency of ANN and its equivalent SNN implementation, we follow the convention from NC community and compute the total synaptic operations $SynOps$ that required to perform a certain task \citep{merolla2014million,ethImageNet,sengupta2019going}. For ANN, the total synaptic operations (Multiply-and-Accumulate (MAC)) per classification is defined as follows
\begin{equation}
SynOps = \sum\limits_{l = 1}^L {{f_{in}^l}}  \cdot {N_l}
\end{equation}
where $f_{in}^l$ denotes the number of fan-in connections to each neuron in layer $l$, and $N_l$ refers to the number of neurons in layer $l$. In addition, $L$ denotes the total number of network layers. Hence, given a particular network configuration, the total synaptic operations required per classification is a constant number that jointly determined by $f_{in}^l$ and ${N_l}$.

While for SNN, as per Eq. \ref{snnSynOps}, the total synaptic operations (Accumulate (AC)) required per classification are correlated with the spiking neurons' firing rate, the number of fan-out connections $f_{out}$ to neurons in the subsequent layer as well as the simulation time window $T$.
\begin{equation}
SynOps = \sum\limits_{t = 1}^T {\sum\limits_{l = 1}^{L-1} {\sum\limits_{j = 1}^{N_l} {f_{out,j}^l} } }  \cdot s_j^l(t)
\label{snnSynOps}
\end{equation}
where $s_j^l(t)$ indicates whether a spike is generated by neuron $j$ of layer $l$ at time instant $t$.

\section{Results}
\label{sec:result}
\subsection{Phone Recognition on TIMIT Corpus}
We report the PER on the development and test sets of TIMIT corpus in Table~\ref{tab:timit}, with numbers in bold being the best performance given by the speaker-independent features. ASR performances of other state-of-the-art systems using various ANN and SNN architectures are given in the upper panel for reference purposes. As the results shown in Table~\ref{tab:timit}, the proposed SNN-based acoustic models are applicable to different speech features and provide comparable or slightly worse ASR performance than the ANNs with the same network structure. In particular, the ANN system trained with the standard 13-dimensional FBANK feature achieves the best PER of 16.9\% (18.5\%) on the development (test) set. The equivalent SNN system using the same feature achieves slightly worse PER of 17.3\% (18.7\%) on the development (test) set. Although the state-of-the-art ASR systems \citep{ravanelli2018} give approximately 1\% lower PER than the proposed SNN-based phone recognition system, it is largely credit to the longer time context explored by the recurrent Li-GRU model.

It worth mentioning that phone recognition is still a challenging task for spiking neural networks. To the best of our knowledge, only one recent work with recurrent spiking neural networks \citep{Bellec738385} demonstrates some promising test results on this corpus with a PER of 26.4\%. In contrast, our system has achieved significantly lower PER compared to this preliminary study of SNN-based acoustic modeling. However, these results are not directly comparable since the proposed system incorporates both an acoustic and a language model during decoding unlike the system described in~\citep{Bellec738385}.

The experimental results on the TIMIT phone recognition task can be considered as an initial indicator of the compelling prospects of the SNN-based acoustic modeling. Given that the phone recognition task on TIMIT corpus is simplistic compared to the modern LVCSR tasks, we further compare the ANN and SNN performance on newer corpora designed for LVCSR experiments.

\subsection{Low-resourced ASR on FAME Corpus}
In this section, we apply the SNN-based ASR systems to the low-resourced ASR scenario. As summarized in Table~\ref{tab:fame}, the word recognition results on the FAME corpus are reported separately for monolingual Frisian (fy), monolingual Dutch (nl) and code-switched (cs) utterances. The overall performance (all) is also included in the rightmost column. Given that 8.5 hours Frisian and 3 hours of Dutch speech is used during the training phase, we can compare the ASR performance on different subsets, i.e. fy, nl and cs, to identify the variations in the ASR performance for different levels of low-resourcedness. We omit the results on the development set as they follow a similar pattern to the results on the test set.

In this scenario, the SNN acoustic models consistently provide lower WERs than the ANN models for all speech features. Systems with the FBANK features provide lower WERs than those using MFCC features, which is in line with our observations on the TIMIT corpus. The best performance on the test set is obtained using SNN models trained on 40-dimensional FBANK features with an overall WER of 36.9\%. In contrast, the ANN model provides a WER of 39.0\% for the same setting, which is relatively 5.4\% worse than the SNN model. Moreover, the SNN-based acoustic models achieve a relative improvement of 4.7\%, 5.2\% and 8.2\% on the fy, nl and cs subsets of the test set, respectively. These steady improvements in the recognition accuracies highlight the effectiveness of the SNN-based acoustic modeling in scenarios with limited training data compared to the conventional ANN models. The improved ASR performance with SNNs, in the low-resourced setting, may credit to the noisy weight updates derived by the coupled ANN layers of the tandem learning framework \citep{hybridlearning}. It has been recognized that introducing noises into the training stage improves the generalization capability of ANN-based ASR systems~\citep{yin2015noisy}. As a result, the noisy training of the tandem learning is expected to improve the recognition performance in low-resourced scenarios. Further investigation on the impact of this noisy training procedure remains as future work.

\subsection{LVCSR experiments on Librispeech Corpus}
In the final set of ASR experiments, we train acoustic models using the official 100-hour and 360-hour training subsets of the Librispeech corpus to compare the recognition performance of ANN and SNN models in a standard LVCSR scenario. As the results given in the middle panel of Table~\ref{tab:libri}, for 100 hours of training data, the ANN systems perform marginally better than the corresponding SNN systems across all different speech features. The absolute WER differences range from 0.1\% to 0.6\%. These marginal performance degradations of the SNN models is likely due to the reduced representation power of using discrete spike counts. Nevertheless, these results are promising even when comparing to the state-of-the-art ASR systems using more complex ANN architectures as provided in the upper panel of Table~\ref{tab:libri}.

It worth noting that both ANN and SNN systems can take benefit of an increased amount of training data. When increasing the training data from 100 hours to 360 hours, the WERs of the best SNN models reduced from 10.0\% (10.3\%) to 9.2\% (9.4\%) for the development (test) sets, respectively. To the best of our knowledge, it is the first time that SNN-based acoustic models have achieved comparable results over the ANN models for LVCSR tasks. These results suggest that SNNs are potentially good candidates for acoustic modeling.

\subsection{Energy Efficiency of SNN-based ASR Systems}
In addition to the promising modeling capability, the SNN-based ASR systems can achieve unprecedented performance gain when implemented on the low-power neuromorphic chips. In this section, we shed light on this prospect by comparing the energy efficiency of ANN- and SNN-based acoustic models. Given that data movements are the most energy-consuming operations for data-driven AI applications, we calculate the average synaptic operations on 5 randomly chosen utterances from the TIMIT corpus and report the ratio of average synaptic operations required per feature classification (SynOps(SNN)\,/\,SynOps(ANN)). To investigate the effect of different feature representations, we repeat our analysis on the 40-dimensional MFCC, FBANK and FMLLR features as summarized in Table~\ref{synops} and Figure~\ref{spikeRate}.

Taking advantage of the short encoding time window ($T_e=10$), the sparse neuronal activities are observed for all network layers as shown in Figure~\ref{spikeRate}. Among the three features explored in this experiment, it is interesting to note the FMLLR feature achieves the lowest average spike rate. It is likely due to the more discriminative nature of the speaker-dependent feature, while it worth to note that the FMLLR feature is not always available in all ASR scenarios. As provided in Table~\ref{synops}, the SNN implementations taking MFCC, FBANK and FMLLR input features require 1.72, 1.10 and 0.68 times synaptic operations to their ANN counterparts, respectively. Although the average number of synaptic operations required for SNNs that using MFCC and FBANK features are slightly higher than the ANNs, the AC operations performed on SNNs are much cheaper than the MAC operations required for ANNs. One recent study on the Global Foundry 28 nm process has revealed that MAC operations are 14 times more costly than AC operations and requires 21 times more chip area~\citep{ethImageNet}. Therefore, when deploying SNNs onto the emerging neuromorphic chips for inference~\citep{merolla2014million, davies2018loihi}, we expect to receive at least an order of magnitude energy and chip area savings. While the actual energy savings for SNN-based acoustic models are dependent on the chip architectures and materials used, which is outside the scope of this work.

\section{Discussion}
\label{sec:dis}
The remarkable progress in the automatic speech recognition systems has revolutionized the human-computer interface. The rapid growing demands of ASR services have raised concerns on computational efficiency, real-time performance, and data security, etc. It, therefore, motivates novel solutions to address all those concerns. As inspired by the event-driven computation that observed in the biological neural systems, we explore using brain-inspired spiking neural networks for large vocabulary ASR tasks. For this purpose, we proposed a novel SNN-based ASR framework, wherein the SNN is used for acoustic modeling and map the frame-level features into a set of acoustic units. These frame-level outputs will further integrate the word-level information from the corresponding language model to find the most likely word sequence corresponding to the input speech signal.

\subsection{Superior Speech Recognition Performance with SNNs}
The phone and word recognition experiments on the well-known TIMIT and Librispeech benchmarks have demonstrated the promising modeling capacity of SNN acoustic models and their applicability to different input features. These preliminary results have shown that the recognition performance of SNNs is either comparable or slightly worse than the ANNs with the same network architecture on the TIMIT and Librispeech benchmarks. A possible reason for this performance degradation is the reduced representation power of the discrete neural representation (i.e., spike counts) as compared to the continuous floating-point representation of the ANNs \citep{hybridlearning}. This performance gap could potentially be closed by extending the encoding window $T_e$ of SNNs. Moreover, the recognition performance of ANN and SNN models in a low-resourced scenario is also investigated. In this scenario, the SNN acoustic models outperform the conventional ANNs that could be attributed to the noisy training of the tandem learning framework, wherein error gradients of the SNN layers are approximated from the coupled ANN layers.

The neural encoding scheme adopted in this work allows input features to be encoded inside a short encoding time window for rapid processing by SNNs. It is attractive for the time-synchronous ASR tasks that require real-time performance. The preliminary study of the energy efficiency on the TIMIT corpus reveals at least an order of magnitude energy and chip area savings, as compared to the equivalent ANNs, can be achieved when deploying the offline trained SNNs onto neuromorphic chips. The recent study of a keyword spotting task on the Loihi neuromorphic research chip \citep{blouw2019benchmarking} has also demonstrated the compelling energy savings, real-time performance and good scalability of emerging NC architectures over conventional low-power AI chips designed for ANNs.

\subsection{Development of SNN-based ASR Systems}
The active development of open-source software toolkits plays a significant role in the rapid progress of ASR research, instances include the Kaldi \citep{povey2011kaldi} and ESPnet \citep{watanabe2018espnet}. In this work, we demonstrate that state-of-the-art SNN acoustic models can be easily developed in PyTorch and integrated into the PyTorch-Kaldi Speech Recognition Toolkit \citep{pytorch-kaldi}. This software toolkit integrates the efficiency of Kaldi and the flexibility of PyTorch, therefore, it can support the rapid development of SNN-based ASR systems.

\subsection{Future Directions}
The recurrent neural networks have shown great modeling capability for temporal signals by exploring long temporal context information in the input signals \citep{graves2014towards}. As future work, we will explore the recurrent networks of spiking neurons for large-vocabulary ASR tasks to further improve the recognition performance.

The substantial research efforts are devoted to reducing the computational cost and memory footprint of ANNs during inference, instances include network compression \citep{han2015deep}, network quantization \citep{courbariaux2016binarized, zhou2016dorefa} and knowledge distillation \citep{hinton2015distilling}. While the computational paradigm underlying the efficient biological neural networks is fundamentally different from ANNs and hence fosters enormous potentials for neuromorphic computing architectures. Furthermore, grounded on the same connectionism principle, the information of both ANN and SNN are encoded in the network connectivity and connection strength. Therefore, SNN can also take benefits from these early research works on the network compression and quantization of ANNs to further reduce its memory footprint and computation cost \citep{deng2019comprehensive}.

The event-driven silicon cochlea audio sensors \citep{ref45} are designed to mimic the functional mechanism of human cochlea and transform input audio signals into spiking events. Given temporally sparse information is transmitted in the surrounding environment, these sensors have shown greater coding efficiency than conventional microphone sensors \citep{liu2019event}. There are some interesting preliminary ASR studies explore the input spiking events captured by these silicon cochlea sensors \citep{ref69, acharya2018comparison}. However, the scale of the ASR tasks explored in these studies is relatively small comparing to modern ASR benchmarks due to the limited availability of event-based ASR corpora. Pan et al. \citep{pan2019efficient} recently proposed an efficient and perceptually motivated auditory neural encoding scheme to encode the large-scale ASR corpora collected by microphone sensors into spiking events. With this encoding scheme, approximately 50\% spiking events can be reduced with negligible interference to the perceptual quality of inputs audio signals. Taking benefits from these earlier research on the neuromorphic auditory front-end, we are expecting to further improve the energy efficiency of SNN-based ASR systems.

The promising initial results demonstrated by the SNN-based large vocabulary ASR systems in this work is the first step towards a myriad  opportunities for the integration of state-of-the-art ASR engines into mobile and embedded devices with power restrictions. In the long run, the SNN-based ASR systems are expected to take benefits from ever-growing research on novel neuromoprhic auditory front-end, SNN architectures, neuromorphic computing architectures and ultra-low-power non-volatile memory devices to further improve the computing performance.

\section*{Conflict of Interest Statement}
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

\section*{Author Contributions}
JBW and EY designed and conducted all the experiments. All authors contributed to the results interpretation and writing.

\section*{Funding}
This research is supported by Programmatic Grant No. A1687b0033 from the Singapore Government's Research, Innovation and Enterprise 2020 plan (Advanced Manufacturing and Engineering domain). JBW is also partially supported by the Zhejiang Lab’s International Talent Fund for Young Professionals.

\bibliographystyle{frontiersinSCNS_ENG_HUMS}


\section{Introduction}
\label{sec:Introduction}
Quantum computing and quantum communication typically focus on quantum information encoded and processed with quantum bit (qubit) strings,
but replacing qubits by higher-dimensional qudit strings~\cite{GKP01,BDS02} can be advantageous~\cite{EFK+18} for quantum simulation~\cite{NAB+09}, quantum algorithms~\cite{LBA+09,TV16,BRS17}, quantum error correction~\cite{DGP13,MSB+16,MLZ+18}, universal optics-based quantum computation~\cite{NCS18}, quantum communication~\cite{CDB+19,LZE+19} and fault-tolerant quantum computation~\cite{CAB12,CAM14}.
Qudit quantum-information process could reduce space requirements and exploit natural properties such as orbital angular momentum for photons~\cite{BEW+17}, superconductors~\cite{NAB+09} and neutral atoms~\cite{Hec17}.
Specifically, quantum computing on higher-dimensional systems can be more efficient than on qubits~\cite{BRS17,CAB12,GEZ+19,GP13}.
Moreover, there exist genuine entangled states on higher-dimensional systems that cannot be simulated by
the tensor product of pairwise entangled qubit states~\cite{KRB+18}.
Ultimate success of quantum computing, both qubit- and qudit-based, depends on being scalable,
which,
in turn,
requires components meeting fault-tolerance conditions~\cite{Sho96}.

Unitary-gate randomized benchmarking (URB)
is the preferred technique to characterize unitary-gate performance due to its efficiency~\cite{MGE11,MGE12},
which is robust against state-preparation-and-measurement (SPAM) errors
and exponentially superior to the alternative of quantum process tomography (QPT)~\cite{MRL08,JLQ+18}.
URB estimates average fidelity between real and ideal implementation of all~$24$ Clifford gates~in the single-qubit Clifford group, $\mathcal{C}_2$,
which normalizes the Pauli group
\begin{equation}
    \bm\Xi_2:=\langle\text{i}
        \mathds1;X_d;Z_d\rangle
\end{equation}
for angular brackets denoting the generating set~\cite{NC10}.
Average unitary-gate fidelity is obtained
by estimating survival decay rate vs gate-sequence length~\cite{MGE12}.
RB is well developed for qubits but untouched for qudits;
here we introduce qudit URB.

We develop qudit URB by combining
single- and multi-qubit URB theory~\cite{KLR+08,MGE11,MGE12,JS14,ATB16,PRY+17,HFGW18,PCR+18}
with qudit
(Hilbert-space dimension~$d$ with $d=2$ for the qubit)
theory including the generalized Pauli group
$\mathcal{P}_d$~\cite{GKP01,BDS02},
the qudit Clifford group~$\mathcal{C}_d^n$~(\ref{eq:Cdn})~\cite{Got99,App05,Gro06}
for~$n$ qudits
and its connection to the unitary 2-design (U2D)~\cite{GAE07,DCEL09}.
U2Ds are valuable as they enable efficient sampling of a random unitary matrix.
Averaging over a U2D uniform distribution is identical to averaging over the unitary group over the uniform (Haar) measure.

Na\"{i}vely,
qudit URB could be regarded as trivially arising from multi-qubit URB~\cite{MGE12},
which we show is not so,
thereby justifying utilizing U2D properties of the qudit Clifford group.
Our qudit URB scheme reduces to qubit URB as a special case.

Our paper is organized as follows. In \S\ref{sec:Background}, we review background knowledge on generalized Pauli group, qudit Clifford group and U2D for qudits.
Our approach is detailed in \S\ref{sec:Approach}.
We present our results in \S\ref{sec:Results}.
Finally,
\S\ref{sec:Discussion} and \S\ref{sec:Conclusion} provide our discussion and conclusion, repectively.

\section{Background}
\label{sec:Background}
In this section, we provide the required background to address the qudit URB.
As discussing the Clifford group without first an overview of Pauli group is not complete,
we begin the section by explaining generalized Pauli group.
Then we proceed with a discussion of qudit Clifford group.
Finally,
we explain U2D and discuss that qudit Clifford group forms a U2D, which is the most important mathematical concept for qudit URB.
\subsection{Generalized Pauli group}
\label{Generalized Pauli group}
In this subsection, we describe generalization of Pauli group for qudits.
First, we begin with the mathematical concept of qudits and introducing generalized Pauli operations.
Then we explain generalized Pauli group.
Finally, we describe this group for $n$ qudits.

Mathematically, a qudit is a vector
in $d$-dimensional Hilbert space~$\mathscr{H}_d\cong\mathbb{C}^d$
spanned by the orthonormal computational basis~$\{\ket{s}; s\in\mathbb{Z}_d\}$ where
$\mathbb{Z}_d:=\{0,1,\dots,d-1\}$.
Qudit unitary transformations are represented by
unitary matrices $\{U_d\}\in U(d)$,
with~$U(d)$ the $d$-dimensional unitary Lie group.
These unitary transformations include generalized Pauli operations~$\bm\Xi_d$, namely,
\begin{equation}
    X_d\ket{s}=\ket{s\oplus1},\,
    Z_d\ket{s}=\omega^s\ket{s},\,
    \omega:=\exp\left(2\pi\text{i}/d\right)
\end{equation}
being defined by their actions on computational basis states
with $\oplus$
denoting addition modulo~$d$~\cite{GKP01,BDS02}. Both $X_d$ and $Z_d$ have order d, obeying $X^d=Z^d=\mathds{1}_d$.
The usual Pauli operators arise for $d=2$.

The generalized Pauli operators $X_d$ and $Z_d$ satisfy the commutation relation
\begin{equation}
    Z_dX_d=\omega X_dZ_d.
\end{equation}
These operators generate the generalized Pauli group
\begin{equation}
    \mathcal{P}_d:=\langle \tilde{\omega} \mathds{1}_d, X_d,Z_d\rangle,
\end{equation}
for $\tilde{\omega}=\omega$ if $d$ is odd
and $\tilde{\omega}=\omega^{1/2}$ if $d$ is even. If $d$ is odd, the order of $Z_dX_d$ is $d$,
whereas, if~$d$ is even,
then the order of~$Z_dX_d$
is~2$d$, which contributes additional roots of unity.
Therefore, $\tilde{\omega}$ is defined differently for~$d$ even vs odd. The difference between the generalized and usual Pauli operators is that the qudit operators for $d>2$ are only unitary, and not Hermitian.
Hence, any eigenvalue of a qudit Pauli operator can be estimated only via quantum phase estimation algorithm~\cite{NC10}.

The $n$-qudit Hilbert space and the associated state (density operator) space are denoted $\mathscr H_d^{\otimes n}$ and $\mathcal D\left(\mathscr H_d^{\otimes n}\right)$, respectively.
For $n$ qudits,
the generalized Pauli group is
$\mathcal{P}_d^n:=\mathcal{P}_d^{\otimes n}$.
This generalized Pauli group further generalizes to
$\widetilde{\mathcal{P}}_d^n
    :=\mathcal{P}_d^n/\langle\tilde{\omega}\mathds{1}_d\rangle$
without phases for $\langle\tilde{\omega}\mathds{1}_d\rangle$, which indicaties the group generated by $\tilde{\omega}\mathds{1}_d$~\cite{Web16}.
\subsection{Qudit Clifford group}
\label{subsec:qd Clifford group}
In this subsection, we explain $n$-qudit Clifford group. First we begin the section by describing the Clifford group for qudits. Then we explain the simple case of $d=3$.

As the normalizer of the generalized Pauli group,
the qudit Clifford group comprises all unitary operators that map $\mathcal{P}_d^n$ to itself under conjugation.
Hence, the $n$-qudit Clifford group~(\ref{eq:Cdn})
is~\cite{Web16}
\begin{align}
\label{eq:Cdn}
   \mathcal{C}_d^n
    :=&\{c\in U(d^n) ; c \widetilde{\mathcal{P}}_d^nc^{\dagger}\in \mathcal{P}_d^n \}/ \{\text{e}^{\text{i}\theta}\mathds{1}_d^n; ~\theta\in \mathbb{R}\nonumber\}\\
   =&\langle \text{CZ}_d, \text{F}_d, \text{P}_d, Z_d \rangle,
\end{align}
i.e.,
generated by
controlled-$Z$
\begin{equation}
    \mathrm{CZ}|ss'\rangle:=\omega^{ss'}|ss'\rangle,
\end{equation}
quantum Fourier transform
\begin{equation}
    \mathrm{F}|s\rangle:=\frac{1}{\sqrt{d}}\displaystyle\sum_{s'\in \mathbb{Z}_d}\omega^{ss'}|s'\rangle,
\end{equation}
phase gate
\begin{equation}
    \mathrm{P}|s\rangle:=\omega^{\frac{s(s+\varrho_d)}{2}}|s\rangle,\; \varrho_d= \begin{cases}
1, &\text{if $d$ is odd,}\\
0, &\text{otherwise,}
\end{cases}
\end{equation}
and Pauli-Z gates~\cite{Far14,HDD05,Pro19}.
Any $n$-qudit Clifford gate can be decomposed into multiplicative and tensor products of these gates.
The cardinality of the single-qudit Clifford group has been explicitly calculated~\cite{Tol18}.
For $d\in\mathbb{P}$ (prime),
the number of distinct Clifford gates (up to global phase) for the single-qudit case is $d^3\left(d^2-1\right)$~\cite{App05}.

For $d=3$, $\mathcal{C}_3$ comprises 216 Clifford gates~\cite{GRT19},
which are generated by single-qutrit Fourier transform~$F_3$ and phase gate~$S_3$.
Specifically, any single-qutrit Clifford operation can be obtained by a product of three elements in
$\mathscr L$,~$\mathscr M$ and~$\mathscr N$~\cite{GRT19},
where~$\mathscr L$ is the subgroup of~$\mathcal C_3$ generated by~$P_3$ and~$X_3$
and
\begin{equation}
    \mathscr M=\{\mathds 1, F_3^2\},\,
    \mathscr N=\{\mathds 1, F_3, S_3F_3, S_3^2 F_3\}
\end{equation}
for~$F_3^2$ and~$S_3^2$
are the squares of~$F_3$ and~$S_3$, respectively.

 Similar to qubit Clifford circuits, quantum circuits with only prime-dimensional qudit Clifford gates
can be classically efficiently simulated~\cite{ME12}.
To achieve universal quantum computing,
at least one non-Clifford gate must be added to the set of Clifford gates.
An example of single-qudit non-Clifford gates is the generalized T-gate~\cite{Pro19}
\begin{equation}
\label{eq:T}
    T\ket{s}
        :=\omega^{s^3/d^2}\ket{s}.
\end{equation}
For qubits, T-gates can be benchmarked by dihedral-group benchmarking~\cite{CWE15}, instead of benchmarking U2D,
which is our focus.
Whether dihedral benchmarking can be generalized to the qudit case is an open problem.

\subsection{Unitary 2-design for qudits}
\label{subsec:2-design}
In this subsection, we explain U2D, which is an essential concept for a scalable and efficient URB.
Then we discuss U2D for the multi-qudit Clifford group.

A U2D comprises unitary matrices
$\{U_{j}\}_{j=1}^K$
satisfying~\cite{DCEL09}
\begin{equation}\label{eq:2-design}
\frac{1}{K}\displaystyle\sum_{j=1}^K U^{\dagger}_{j}\mathcal E(U_{j}\rho U^{\dagger}_{j})U_{j}=\int_{U(d)}\text{d}U\; U^{\dagger}\mathcal E(U\rho U^{\dagger})U,
\end{equation}
for any quantum channel~$\mathcal E$~\cite{wil11}
and any state~$\rho$
with~$\text{d}U$ denoting the unitarily invariant Haar measure~\cite{BR86}
on the Lie group $U(d)$.
Eq.~(\ref{eq:2-design}) implies that twirling any quantum channel over U2D is equivalent to twirling over a Haar-measure unitary group.
The multi-qubit Clifford group $\mathcal{C}_2^n$ for~$n$ qubits forms a U2D~\cite{DCEL09,DDB02}.

Though not explicitly stated, the multi-qudit Clifford group $\mathcal{C}_d^n$
evidently forms a U2D based on Webb's analysis of Pauli-mixing Clifford ensembles~\cite{Web16}.
A set $\mathcal{S} \subseteq \mathcal{C}_d^n$ is called a \emph{Pauli-mixing Clifford ensemble}
(Def.~3 in~\cite{Web16})
if every pair of non-identity Pauli operators are related,
up to a phase,
by a Clifford conjugation
and the number of Clifford operators for each conjugation is constant.
Thus,
a Clifford operator chosen uniformly at random from $\mathcal{S}$ maps
every non-identity Pauli operator to every other non-identity Pauli operator with equal probability.
Webb's Lemma~3, together with Lemma~2, says that $\mathcal{C}_d^n$ is both Pauli-invariant
(for each $c\in\mathcal S$,
$c \Xi$, up to a phase, is in~$\mathcal S$
as well for all $\Xi\in\mathcal P$)
and Pauli-mixing, and
Lemma~1 says that, if a Pauli-invariant Clifford ensemble is Pauli-mixing,
then this ensemble is a U2D.
Hence, $\mathcal{C}_d^n$ is a U2D.
\section{Approach}
\label{sec:Approach}
In this section, we present our approach to standard URB protocol.
We use the U2D property for the multi-qudit Clifford group to prove twirling over this group is depolarizing.
Then, we discuss averaged sequence fidelity obtained from depolarizng parameter.
The section ends with explaining the fitting model for URB.
\subsection{Twirling over qudit Clifford group}
In this subsection we show how exploiting U2D property of $n$-qudit Clifford group result in depolarizing channel. Twirling a quantum channel with respect to a group of unitary operations is the basic approach utilized by URB.
Twirling~$\mathcal E$ with average fidelity
\begin{equation}
    \bar{F}_{\mathcal E}=\int_{\mathscr H_d} \text{d} \phi \bra{\phi}\mathcal E\left(\ket{\phi}\bra{\phi}\right)\ket{\phi}
\end{equation}
over Haar-random unitary operations
yields a depolarizing channel~\cite{HHH99,Nie02}
\begin{equation}
    \int_{U(d)}\text{d}U\;U^{\dagger}\mathcal E(U\rho U^{\dagger})U=\mathcal E_{\text{dep}}(\rho)
\label{twirled channel},
\end{equation}
with the same average fidelity as for~$\mathcal E$,
where
\begin{equation} \label{relationBetweenpandF}
    \mathcal E_{\text{dep}}(\rho):= p\rho+(1-p)\frac{\mathds1}{d^n},\,
    p=\frac{d^n \bar{F}_{\mathcal E}-1}{d^n-1}.
\end{equation}

\begin{prop}\label{prop:twirling}
Twirling a channel over an $n$-qudit Clifford group yields a depolarizing channel.
\end{prop}
\begin{proof}
Combining Eq.~(\ref{twirled channel})
with the fact
that the $n$-qudit Clifford group forms a U2D
yields
\begin{equation}
\label{twirlCliffordDepolarizing}
\frac{1}{|\mathcal C^n_d|}\sum_{l=1}^{|\mathcal C^n_d|} \mathscr{C}_l^{-1} \circ \mathcal E\circ \mathscr{C}_l=\mathcal E_{\text{dep}}.
\end{equation}
Averaging over the $n$-qudit Clifford group is identical to averaging over the unitary group with respect to the uniform Haar measure
so Eq.~(\ref{twirlCliffordDepolarizing})
follows.
\end{proof}
\subsection{Averaged sequence fidelity}
In this subsection, first, we set some notation that will be used throughout.
Then we discuss how to generate a random sequence of Clifford gates, and obtain averaged sequence fidelity.

A concatenation,
or composition,
of $n$-qudit Clifford gates in a sequence of length $m$ is denoted
$\bigcirc_{j=1}^{m} \mathscr{C}_{i_j}$
for $\bigcirc$ indicating concatenation.
We use $i_j$ as labels for the $j^\text{th}$ element in the $i^\text{th}$ sequence, where $i_j\in\left[\left|\mathcal{C}_d^n\right|\right]$,
with $[k]:=\{1,\dots,k\}$,
and
the gate is denoted~$\mathscr{C}_{i_j}$.
All Clifford gates experience the same noise,
which is represented by a noisy channel~$\Lambda$ following an ideal Clifford gate.

 Averaging over random realizations of a sequence of Clifford gates
\begin{equation}
    S_{\bm{i}_m}=\bigcirc_{j=1}^m
    \Lambda \circ \mathscr{C}_{i_j}, {\bm{i}_m}:=\left(i_1,i_2,\dots,i_m\right)
\end{equation}
with ${\bm{i}_m}$ denoting an $m$-tuple and $\mathscr{C}_{i_m}=\left(\bigcirc_{j=1}^{m-1} \mathscr{C}_{i_{j}}\right)^{-1}$, is equivalent to concatenating
$m-1$ twirled channels~\cite{MGE11,MGE12}
\begin{equation}
      \Lambda_{\textrm{T}}:=\frac{1}{|\mathcal C^n_d|}\sum_{l=1}^{|\mathcal C^n_d|} \mathscr{C}_l^{-1} \circ \Lambda\circ \mathscr{C}_l,
\end{equation}
followed by~$\Lambda$,
i.e., $\Lambda\circ \Lambda_{\textrm{T}}^{\circ m-1}$.
Using Proposition~\ref{prop:twirling}, $\Lambda_{\textrm{T}}^{\circ m}$ can be rewritten as an $m$-fold composition of a depolarizing channel with itself multiple times,
namely,
\begin{equation}\label{mfolddepchannel}
    \Lambda_{\textrm{T},p}^{\circ m-1}(\rho)=p^{m-1}\rho+(1-p^{m-1})\mathds{1}/d^n.
\end{equation}
Hence, for any input state~$\ket{\psi}\in \mathscr{H}_d^{\otimes n}$, \begin{equation}\label{eq:IdealSurvivalProbability}
    \tr\left[\ket{\psi}\bra{\psi}\Lambda\circ \Lambda_{\textrm{T},p}^{\circ {m-1}}(\ket{\psi}\bra{\psi})\right]
\end{equation}
is channel fidelity averaged over random realizations of the sequence.
\subsection{Fitting model}
In this subsection, we present the fitting function for URB, by which we model the behaviour of averaged sequence fidelity.
In practice,
with quantum noise, Eq.~(\ref{eq:IdealSurvivalProbability}) is replaced by
\begin{equation}
\label{eq:rearrangeSurvivalProbability}
    F_\text{seq}(m):=\tr\left[E_\psi\Lambda\circ \Lambda_{\textrm{T},p}^{\circ m-1}(\rho_\psi)\right]
\end{equation}
for $E_\psi$ and $\rho_\psi$ the positive-operator valued measure (POVM)~\cite{NC10} element and quantum state
including SPAM errors, respectively.
Plugging Eq.~(\ref{mfolddepchannel}) into Eq.~(\ref{eq:rearrangeSurvivalProbability}) yields
\begin{equation}
\label{fittingModel}
    F_\text{seq}(m)=A_0 p^{m-1}+B_0,
\end{equation}
which absorbs SPAM errors,
for
\begin{equation*}
A_0:=\tr\left[E_\psi\Lambda\left(\rho_\psi-\frac{\mathds{1}}{d^n}\right)\right],\,
B_0:=\frac{\tr\left[E_\psi\Lambda(\mathds{1})\right]}{d^n}
\end{equation*}
being the coefficients.
\subsection{Summary of approach}
\label{subsec:summary}
In this section, exploiting U2D property of $n$-qudit Clifford group we presented twirling a channel over this group yields a depolarizing channel. Then, we obtained averaged sequence fidelity from depolarizing rate. Furthermore, we introduced fitting model for our URB.
This section relies on the assumption that quantum noise is gate-independent.
As gate-dependent noise decays in the same form
as gate-independent noise plus a perturbation~\cite{W18},
our approach could be naturally extended to the case of gate-dependent noise.

\section{Results}
\label{sec:Results}
In this section we present our main results. The first result of this paper is providing pseudocode for URB procedure. We proceed to explain the URB procedure for qudit Clifford gates. Then, giving a counterexample, we explain why multi-qubit URB does not readily yield qudit URB.
\subsection{Randomized benchmarking as an algorithm}
\label{subsec:algorithm}
In this subsection, we provide pseudocode for our multi-qudit URB procedure,
which is immensely useful to ensure that the procedure flows logically and does not leave out any key steps.
Our algorithm
is designed to estimate average gate fidelity~$\bar{F}_\Lambda$
over a Haar-random set of input states.

Our pseudocode uses the following data types,
expressed conventionally as all capitals.
UNSIGNED INTEGER refers to a positive integer in~$\mathbb{Z}^+$,
and REAL, COMPLEX, BINARY, DARY and INTERVAL refer to real~$\mathbb R$, complex~$\mathbb C$,
binary~$\{0,1\}$,
d-ary~$\{0,1,\dots,d-1\}$
and the unit interval~$[0,1]\subset\mathbb{R}$,
respectively.
Besides classical data types, we introduce quantum data types as well~\cite{knil96}.
QDARY
refers to a qudit of dimension~$d$.
Each of these data types can be an array
with data type followed by brackets~$[~]$;
a sequence of two brackets $[~][~]$
denotes a two-dimensional array,
which is readily generalized to higher-dimensional arrays by adding more brackets.

For pseudocode variables, we use \textsc{camelCase},
so qudit number~$n$
is denoted by \textsc{numQud}
and of type UNSIGNED INTEGER,
and Hilbert-space dimension~$d$
is \textsc{HilbDim}
and also of type UNSIGNED INTEGER.
We use \textsc{numSeq} to denote the maximum number of different gate sequences of a fixed length and of type UNSIGNED INTEGER.
QDARY[~] indicates a multi-qudit state and
QDOP is an operation that maps a QDARY[\textsc{HilbDim}]-typed variable to another QDARY[\textsc{HilbDim}]-typed variable.

In describing our algorithm using pseudocode,
we employ functions from an ideal library explained here.
We use \textsc{rand}(\textsc{maxInt}) to generate a uniformly random integer in [\textsc{maxInt}].
\textsc{prod} maps two Clifford-gate indices to the index corresponding
to the product of these two referenced Clifford gates.
\textsc{inv} maps one Clifford-gate index to the index corresponding to the
inverse of the Clifford gate.
We use data type \textsc{prep}
for preparing a qudit pure state according to a classical description of the state, and
 \textsc{projMeas} denotes qudit-state measurement that yields~$1$ if the qudit state is projected onto a certain pure state and otherwise yields~$0$.
For statistical processes,
we employ \textsc{avg},
which calculates the average of all entries in an array,
and \textsc{fit},
which is a least-squares regression algorithm.
Our randomized benchmarking algorithm comprises input, output and procedure,
which we now describe in plain English.
\subsubsection{Input}
We begin by explaining the input.
The input can be separated into two components,
those that are necessary to specify the benchmarking task and
those that are necessary to specify the benchmarking procedure.
To specify the benchmarking task, the user must specify
the number of qudits $n \in \mathbb{Z}^+$,
the dimension of the Hilbert space $d \in \mathbb{Z}^+$,
and the cardinality $\left|\mathcal C^n_d\right|\in\mathbb{Z}^+$
of the $n$-qudit Clifford group.
Cardinality growing quickly for qudits implies significant experimental challenges.
Circuit depth for testing just the qutrit Clifford group, comprising 216 elements,
would be challenging.
We need a length-$\left|\mathcal C^n_d\right|$ array
of labels for Clifford operators
in order to be able to refer to them individually.

For the benchmarking procedure,
our algorithm caters to an experienced client
who is able to guess good parameters for the number~$k\in\mathbb{Z}^+$ of random Clifford
sequences,
the number $l\in\mathbb{Z}^+$ of repetitions of each
Clifford sequence
and the maximum length $m\in\mathbb{Z}^+$ of a random
Clifford sequence to be executed.
This requirement that the client be able to select good parameters is typical for all qubit-based quantum benchmarking~\cite{JS14}.
Furthermore,
the client is expected to know that the noise model is specified
by an unknown CPTP map, and therefore knows that
the twirled noise model is entirely specified by a single unknown depolarising parameter.
The client aims to estimate this unknown parameter to within a target confidence,
which is not an input to the algorithm.
\subsubsection{Output}
The output of the algorithm is an estimate~$\hat{r}$ of the average gate infidelity, often called average error rate,
\begin{equation}
    r:=1- \bar{F}_\Lambda\in[0,1].
\end{equation}
This estimate is a URB figure of merit
that characterizes average performance of qudit Clifford gates.
\subsubsection{Procedure}
Now we explain the URB procedure
for qudit Clifford gates.
We initialize the $n$-qudit state~$\ket{\psi}$
as the pure state $\ket0^{\otimes n}$.
Then we generate $k$ random sequences of $n$-qudit Clifford gates
$\mathscr{C}_{i_j}$
each of length~$j$,
where $2\le j\le m$, as~$k$ samples of a random sequence. The first~$j-1$ gates in each sequence are uniformly randomly chosen from~$\mathcal C_d^n$,
and the final gate is determined by the first~$j-1$ gates according to
\begin{equation}
    \mathscr{C}_{i_j}=\left(\bigcirc_{j^\prime=1}^{j-1} \mathscr{C}_{i_{j^\prime}}\right)^{-1}.
\end{equation}
As Clifford gates form a group, this final gate is also an element of the group.

We apply each of the $k$ sequences of qudit gates to the initial state.
Then
we apply measurements corresponding to the POVM~$\left\{\Ket{\psi}\bra{\psi}, \mathds{1}-\Ket{\psi}\Bra{\psi}\right\}$
on the output state.
If the measurement outcome corresponds to~$\ket{\psi}\bra{\psi}$,
we assign a value of one,
otherwise, a value of zero.
By averaging over~$k$ different sequences and~$l$ copies of each sequence, we obtain an estimate $\hat{F}(j)$ of the averaged sequence fidelity
\begin{equation}
\label{survivalprobabilityqubit}
        F_\text{seq}(j)=\tr(E_{\psi}\mathcal S_{j}(\rho_\psi)),\,
    \mathcal S_{j}=\frac{1}{k} \sum_{\bm{i}_{\bm{j}-1}\in \left[\left|\mathcal C_d^n\right|\right]^{\otimes j-1} } S_{\bm{i}_j}.
\end{equation}
Now we repeat the above procedure for different values of~$j$, which increases from two to the maximum length~$m$
in succession.
Finally, we fit the estimates $\hat{F}(j)$ to Eq.~(\ref{fittingModel})
with~$p$ the decay parameter
and~$\hat{p}$ its estimate.

Per Eq.~(\ref{relationBetweenpandF}), we see
that fidelity decay parameter $p$ is related to~$r$
    via
\begin{equation}\label{eq:averageErrorRate}
    r=(1-p)\left(1-\frac{1}{d^n}\right).
\end{equation}
Therefore, by estimating $p$ from URB of Clifford gates, we obtain the output average infidelity (algorithm~\ref{alg:quditRB}).
\begin{algorithm}[H]
\caption{Randomized Benchmarking} \label{alg:quditRB}
\begin{algorithmic}[1]
\Require{%
\Statex UNSIGNED INTEGER \textsc{numQud}
\Comment{\# qudits}
\Statex UNSIGNED INTEGER \textsc{HilbDim}
\Comment{Hilbert-space dimension}
\Statex UNSIGNED INTEGER \textsc{cardCliff} \Comment{Cardinality of \textsc{numQud}-qudit Clifford group }
\Statex UNSIGNED INTEGER \textsc{gateIndex}[\textsc{cardCliff}]
\Comment{Array of multi-qudit gates' labels}
\Statex UNSIGNED INTEGER  \textsc{maxLeng} \Comment{Maximal sequence length}
\Statex UNSIGNED INTEGER \textsc{numSeq}
\Comment{\# sequences for each length}
\Statex UNSIGNED INTEGER \textsc{numCop}
\Comment{\# copies for each gate sequence}
}
\Ensure{\Statex  INTERVAL \textsc{avInfid}
\Comment{Average infidelity}
    }
 \Procedure{randBench}{\textsc{numQud},
 \textsc{HilbDim}, \textsc{CliffGate}, \textsc{maxLeng}, \textsc{numSeq}, \textsc{numCop}}

 \State COMPLEX \textsc{initState}$\left[\textsc{HilbDim}\wedge\textsc{numQud}\right] \gets [0,0,\dots,0]$;
  \Comment{The initial state is the \textsc{numQud}-qudit pure state}
\State REAL \textsc{estRegr}[3];
\Comment{Base, slope and intercept for exponential fiting}
\State REAL \textsc{seqFid}[\textsc{maxLeng}];
 \Comment{Fidelity averaged over gate sequences}
 \State UNSIGNED INTEGER \textsc{currIndex};
    \State QDARY \textsc{state}[\textsc{numQud}];
 \For{$k=2:\textsc{maxLeng}$}
  \State BINARY \textsc{outcome}$\left[\textsc{numCop}\right]$;
 \Comment{Array of measurement outcomes}
    \State  REAL \textsc{survProb}[\textsc{numSeq}];
     \Comment{Array of survival probabilities of the initial state}
 \For{$i=1:\textsc{numSeq}$}
 \State $\textsc{currIndex}\gets 1$;
 \Comment{Overwrite $\textsc{currIndex}$ with $1$, the index refered to indentity channel}
    \For{$j=1:k-1$}
     \State $\textsc{gateIndex}[j]
   \gets \textsc{rand}\left(\textsc{cardCliff}\right)$;
    \Comment{Generate a random integer from one to \textsc{cardCliff}.}
     \State $\textsc{currIndex}\gets \textsc{prod}(\textsc{currIndex},\textsc{gateIndex})$;
    \EndFor
         \State $\textsc{gateIndex}\left[k\right]\gets \textsc{inv}(\textsc{currIndex})$;
     \For{$l=1:\textsc{numCop}$}
    \State $\textsc{state} \gets \textsc{prep}(\textsc{initState}$);
 \For{$j=1:k$}
 \State QDOP \textsc{CliffGate}[\textsc{cardCliff}];
\Comment{Array of multi-qudit gates}
    \State $\textsc{state} \gets \textsc{CliffGate}\left[\textsc{gateIndex}[j]\right]*\textsc{state}$;
    \Comment{$j$th Clifford gate maps old state to new state}
    \EndFor
     \algstore{myalg}
\end{algorithmic}
\end{algorithm}
\begin{algorithm}
\begin{algorithmic}
\algrestore{myalg}
    \If{$\textsc{projMeas}(\textsc{state}, \textsc{initState}) =1$ }
    \State $\textsc{outcome}[l]\gets 1$;\Comment{Assign~$1$ if each single-qudit state is projected onto~\textsc{initState}}
    \Else
    \State $\textsc{outcome}[l]\gets 0$;\Comment{Otherwise, assign~$0$}
    \EndIf
\EndFor
\State $\textsc{survProb}[i]\gets \textsc{avg}(\textsc{outcome})$;
\EndFor
\State $\textsc{seqFid}[k]\gets \textsc{avg}(\textsc{survProb})$;
\EndFor
\State \textsc{estRegr}
$\gets \textsc{fit}\left( \textsc{seqFid}\right)$; \Comment{Least-squares regression algorithm with exponential fitting model $y=\textsc{estRegr}[2]*\textsc{estRegr}[1]\wedge x+\textsc{estRegr}[3]$.}
\State \Return $\textsc{avInfid}\gets (1-\textsc{estRegr}[1])
*\left(1-1/\textsc{HilbDim}\wedge\textsc{numQud}\right)$.
  \EndProcedure
\end{algorithmic}
\end{algorithm}

\subsection{Multi-qubit randomized benchmarking versus qudit randomized benchmarking}
Now we explain why previous work on multi-qubit URB does not readily yield qudit URB.
One might expect that solving multi-qubit URB would yield qudit URB trivially.
Such an approach would exploit Schur-Weyl duality~\cite{GW09}.
Schur-Weyl duality,
applied to the symmetry group~$S_n$ and the unitary group $U(d)$,
which have commuting actions on the $n$-fold tensor product of $d$-dimensional Hilbert spaces,
$\mathscr H_d^{\otimes n}$,
states that,
under the joint action of~$S_n$ and~$U(d)$,
the tensor product space decomposes into a direct sum of tensor products of irreducible modules.
The question is whether we can use that to construct Clifford operators for qudits.
We show that this enticing notion is fallacious by falsifying the following proposal.
\begin{proposal}\label{conj:tensor}
A tensor product of any~$n$ single-qubit Clifford operators is a direct sum of a Clifford operator for the $(n+1)$-dimensional symmetric space with any operators for the remaining partially and antisymmetric spaces.
\emph{(FALSE)}
\end{proposal}
\noindent This proposal is enticing because we could simply use existing multi-qubit benchmarking work~\cite{MGE11} instead of producing a new result.

Mathematically,
this proposal can be expressed as follows.
Let~$\left\{C^{(2)}_\imath\right\}_{\imath=1}^n$
be a sequence of~$n$ single-qubit operators,
and let~$C^{(n+1)}$
be any Clifford operator on an $(n+1)$-dimensional Hilbert space~$\mathscr{H}_d$ for $d=n+1$.
The conjecture is then that,
for all $\left\{C^{(2)}_\imath\right\}_{\imath=1}^n$,
a Clifford operator~$C^{(n+1)}$ exists
such that
\begin{equation}
\bigotimes_\imath C^{(2)}_\imath=C^{(n+1)}\oplus \Theta^{(2^n-n-1)},
\end{equation}
for any $(2^n-n-1)$-dimensional unitary operator $\Theta^{(2^n-n-1)}$.
We now demonstrate that this proposal is false by giving a counterexample.

\begin{proof}[\textsl{Counterexample}]
This proposal is falsified with a counterexample,
specifically for $C^{(2)}_\imath=H$
for $\imath\in\{1,2\}$.
The two-qubit Hadamard gate is
\begin{equation}
    H\otimes H
    =R\oplus(-1),\;
R:=\frac{1}{2}
\begin{pmatrix}
 1&\sqrt2&1\\
 \sqrt2&0&-\sqrt2\\
1&-\sqrt2&1
 \end{pmatrix},
\end{equation}
which is block-diagonal on both the symmetric and anti-symmetric subspaces.
However, the block part $R$ on the three-dimensional symmetric subspace is not a qutrit Clifford gate, as $RX_3R\notin \mathcal P_3$.
\end{proof}
\noindent Therefore, the proposal is falsified:
a tensor product of qubit Clifford gates
cannot in general be written as a direct sum involving Clifford gates,
and qudit Clifford gates are not directly
obtained from multi-qubit Clifford gates over the symmetric subspace. This falsification implies a significant difference between quantum computing on qudits vs on multiple qubits,
even for the same total dimension.

\section{Discussion}
\label{sec:Discussion}
We have explained how we can characterize the average performance of the qudit Clifford gates directly through performing qudit URB. The procedure is similar to the qubit URB, just the unitary operators are chosen from qudit Clifford group. We have designed URB for qudit Clifford gates, by synthesizing the U2D property of them with qubit URB. We also devise a pseudocode, which provides the instructions on how to run randomized benchmarking algorithm on a quantum computer.

The U2D property of the qudit Clifford group indicates that twirling a noisy channel over this group yields a depolarizing channel per Proposition~\ref{prop:twirling}. Hence, analogous to qubit URB, we can relate the depolarizing parameter, estimated from URB procedure, to average error rate of qudit Clifford gates.
On the other hand, the natural question that arises is,
given that multi-qubit URB has already been studied, whether qudit URB could be determined from multi-qubit case by considering the symmetric subspace of multi-qubits. We have explained that this symmetrization in Conjecture~\ref{conj:tensor} fails.
 \section{Conclusion}
\label{sec:Conclusion}
Our results extend previous URB results to higher dimensional qudits for estimating average error rate for gate independent errors, and pave the way for experimental characterization of qudit Clifford gates. Recent development on photonic qudit-based quantum computing~\cite{NCS18} provides a good test bed for our qudit URB.
We suggest that quantum optics will provide a good test by exploiting different
photonic degrees of freedom, for example orbital angular-momentum~\cite{BEW+17}, frequency~\cite{LLP+18,IJA+19}, and time~\cite{IJA+19}.

\ack
BCS acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada and from the National Natural Science Foundation of China (NSFC) Grant No.~11675164.
YDW acknowledges support from the Hong Kong Research Grant Council through Grant No.~17300918.
YRS is funded by Australian Research Council Discovery Project DP190102633.