#  Sup3r: A Semi-Supervised Algorithm for increasing Sparsity, Stability, and
Separability in Hierarchy Of Time-Surfaces architectures.

Marco Rasetto
Department of Bioengineering
University of Pittsburgh, Center of Neural Basis of Cognition
Pittsburgh, USA
<EMAIL_ADDRESS>
&Himanshu Akolkar
Department of Neurobiology
University of Pittsburgh,
Pittsburgh, USA
<EMAIL_ADDRESS>
&Ryad Benosman
Meta Reality Labs,
New York, USA

###### Abstract

The Hierarchy Of Time-Surfaces (HOTS) algorithm, a neuromorphic approach for
feature extraction from event data, presents promising capabilities but faces
challenges in accuracy and compatibility with neuromorphic hardware. In this
paper, we introduce Sup3r, a Semi-Supervised algorithm aimed at addressing
these challenges. Sup3r enhances sparsity, stability, and separability in the
HOTS networks. It enables end-to-end online training of HOTS networks
replacing external classifiers, by leveraging semi-supervised learning. Sup3r
learns class-informative patterns, mitigates confounding features, and reduces
the number of processed events. Moreover, Sup3r facilitates continual and
incremental learning, allowing adaptation to data distribution shifts and
learning new tasks without forgetting. Preliminary results on N-MNIST
demonstrate that Sup3r achieves comparable accuracy to similarly sized
Artificial Neural Networks trained with back-propagation. This work showcases
the potential of Sup3r to advance the capabilities of HOTS networks, offering
a promising avenue for neuromorphic algorithms in real-world applications.

_K_ eywords Neuromorphic $\cdot$ Learning $\cdot$ Rule $\cdot$ Supervised
$\cdot$ Semi Supervised $\cdot$ Continual Learning $\cdot$ Incremental
Learning $\cdot$ HOTS $\cdot$ Time Vectors $\cdot$ Time Surfaces

## Introduction

Hierarchy Of Time-Surfaces is a neuromorphic algorithm used to extract
features from patterns of events [1]. This is possible thanks to a type of
representation called time-surface or time vector, where events are
interpolated by exponential decay kernels and collected to represent relative
time differences between the activation of units in the network. Time surfaces
are one of the most common representations in the neuromorphic field since
they allow to interface event data with traditional machine learning and
computer vision algorithms [2, 3]. In HOTS, time surfaces are clustered
together using algorithms like $k$-means to extract common activity patterns,
and layers of units are built by considering each centroid as a neuron that
can emit a new event when an input time surface is assigned to it. For this
reason, HOTS shares many points in common with bag-of-words or bag-of-features
algorithms[4]. For instance, HOTS requires an external classifier on
histograms of features to classify information. Similarly to bag-of-words
algorithms, HOTS classifiers are histograms that accumulate features over a
given temporal window to produce an input vector to traditional machine
learning algorithms like Support Vector Machines and Multi-Layer
Perceptrons[5, 6, 1, 7]. This approach limits compatibility with neuromorphic
hardware and can nullify latency and energy efficiency advantages that are
found in neuromorphic systems. Compared to Spiking Neural Networks (SNNs)
trained with backpropagation through time, HOTS lags in accuracy [8, 9, 10,
1]. However, feature engineering and model tweaking can often bridge the gap
on selected tasks [9, 8]. This is not surprising as it is known that
clustering algorithms can reach comparable accuracy to small Convolutional
Neural Networks (CNNs) when carefully tuned [11]. Moreover, clustering-based
algorithms like HOTS are still worth investigating because of properties that
might be inherited from clustering, like better explainability [12, 13]
compared to more complex SNNs and continual and incremental learning[14, 15,
16]. This work presents Sup3r, a Semi-Supervised algorithm for increasing
Sparsity Stability and Separability in the Hierarchy of Time-Surfaces. This
algorithm can train end-to-end HOTS networks online, ditching the classifier
for a Deep Neural Network architecture where every last unit encodes for a
different label. We show that Sup3r can learn class-informative patterns and
rejects events from confounding features, reducing the number of processed
events. Moreover, Sup3r allows the HOTS network to adapt to distribution
shifts in the data (continual learning) and learn new tasks without forgetting
(incremental learning). Finally, preliminary results on N-MNIST show that
Sup3r can reach comparable accuracy to a similarly sized Artificial Neural
Network (ANN) trained with back-propagation.

##  1 Time Surfaces and Time Vectors

HOTS originally introduced the concept of time surfaces [1]. These descriptors
encode relative timings between events of units within a defined neighborhood
and were historically used on events generated by event-based cameras [17].
This type of camera does not output frames but an asynchronous stream of
events generated by individual pixels detecting brightness changes. Thus, in
the context of vision, an event $i$ can be described as the tuple:

$ev_{i}=(x_{i},y_{i},p_{i},t_{i})$ (1)

where $t_{i}$ is the timestamp of the event, $x_{i}$ and $y_{i}$ are the
horizontal and vertical coordinates of the pixels generating the event, and
$p_{i}\in[0,1]$ is the polarity of the event, which is 0 when the event
signifies a decrease in brightness for the given pixel and 1 when it signifies
an increase. To generalize to different modalities and uses, we can redefine
events as:

$ev_{i}=(\mathbf{x_{i}},p_{i},t_{i})$ (2)

where $\mathbf{x}$ is now a vector (in bold) representing any possible space
coordinates depending on the event-based sensor type. In the case of events
generated by event-based cameras $\mathbf{x}=(x,y)$ [1], whereas in the case
of neuromorphic cochleas $\mathbf{x}=ch$ representing the channel index [6].
To generalize, we consider $p_{i}\in\mathbb{N}$ as we now use it to index
events produced by multiple units at the same spatial location [1, 6]. Time
surfaces are built by centering a squared window $R_{i}$ of lateral size $l$
around every $\mathbf{x_{i}}$ and collecting the last event emitted by
neighboring pixels in a time context $T_{i}$:

$\mathbf{T_{i}}=\max_{j\leq i}\\{t_{j}|\mathbf{x_{j}}\in R_{i}\\}$ (3)

Finally, time surfaces are built by applying an exponential decay kernel on
the time context:

$\mathbf{ts_{i}}=e^{-\frac{(t_{i}-\mathbf{T_{i}})}{\tau}}$ (4)

Time vectors $\mathbf{tv_{i}}$ are a more general definition of time surface
where the number of dimensions is one or more than two [6]. In HOTS, time
surfaces are clustered with clustering algorithms like $k$-means. Whenever a
time surface is assigned to a centroid, it produces an event. These events can
be used to generate new time surfaces that can be clustered with longer taus
to integrate and extract more complex patterns of activities [1]. Events
outputted by a layer $k$ are so defined:

$ev_{i}^{k+1}=(\mathbf{x_{i}}^{k+1},p_{i}^{k+1},t_{i}^{k+1})$ (5)

where $\mathbf{x_{i}}^{k+1}=\mathbf{x_{i}}^{k}$, $t_{i}^{k+1}=t_{i}^{k}$, and
$p_{i}^{k+1}=f_{i}^{k}$ with $f$ being the "firing" feature/centroid assigned
to the input time surface $ts_{i}^{k}$.

##  2 A hierarchy of time surfaces and time vectors

A fundamental principle behind HOTS networks is the increase of spatial
neighborhood $R^{k}$ and $\tau^{k}$ at every layer to extract longer and more
complex spatial and temporal patterns. In this work, we use the same
principle, with the difference of the last layer $K$ (where $K$ denotes the
last layer of a network with layer index $k$) where $R^{K}$ is always
including all coordinates $\mathbf{x}$. This is somewhat different from
traditional HOTS implementation, where it is often the external classifier
that integrates all the possible spatial coordinates.

##  3 Sparsity, Separability, and Stability in spiking networks

Sup3r aims to remove external classifiers from HOTS by using the network’s
last layer as a classifier for a pattern recognition problem. In that case, we
expect the number of clusters to equal the number of classes in the problem
and for each cluster to generate events only when the correct class is
presented to the network. In this context, the classifier will have to limit
the firing to only one unit (high Sparsity) for the whole duration of the
stimuli (high Stability). In the context of clustering, this also means that
classes assigned to different centroids must be well Separated. Conveniently,
we can represent these three concepts using a time vector called feedback-
time-vector $\mathbf{ftv}^{K}$, built on the output events $ev_{i}^{K+1}$ of
our classifier. Rather than pulling the last event in a given neighborhood $R$
(eq. 3), we pull the last event of every polarity in the last layer:

$\mathbf{FT_{i}^{K}}=\max_{j\leq i}\\{t_{j}^{K+1}|\mathbf{p_{j}^{K+1}}\\}$ (6)

We use $\mathbf{ftv}^{k}$ to differentiate feedback time vectors from standard
time vectors $\mathbf{tv}^{K}$ or time surfaces $\mathbf{ts^{K}}$ built on the
input events $ev_{i}^{K}$. We then use it to build the descriptor $S$:

$S^{K}_{i}=G(\mathbf{ftv}^{K}_{i}[f]-\sum\limits_{l\neq
f}^{N^{K}}\frac{\mathbf{ftv}^{K}_{i}[l]}{K-1})$ (7)

Where $f$ is the index of the feedback-time-vector element corresponding to
the ”firing” centroid (last event) that triggered its generation, and $N^{K}$
is the number of clusters for layer $K$. $G$ is a sign function that is
positive if the ”firing” centroid corresponds to the correct class label,
negative if otherwise:

$G=\begin{cases}+&\text{for }f=cc\\\ -&\text{for }f\neq cc\end{cases}$ (8)

where $cc$ is the index of the correct class. To see how this descriptor
relates to Sparsity, Stability, and Separability, let us look at the limits
for $S$. If the classifier has a stable, sparse, and correct response (only
one unit active, corresponding to the correct class index $cl$), $S$ will tend
to 1:

$S^{K}_{i}=+(1-\sum\limits_{l\neq f}^{N^{K}}\frac{0}{K-1})=1$ (9)

Conversely, if only one unit is active, but it corresponds to the wrong class
label, G will be negative:

$S^{K}_{i}=-(1-\sum\limits_{l\neq f}^{N^{K}}\frac{0}{K-1})=-1$ (10)

indicating a fully sparse and stable but wrong response, meaning the
classifier does not correctly separate the classes.

In case the activity of the classifier is dense, meaning that all clusters are
firing events at approximately the same time, the feedback time vectors will
become fully populated by ones, making $S=0$:

$S^{K}_{i}=G(1-\sum\limits_{l\neq
f}^{N^{K}}\frac{1}{K-1})=G(1-\frac{K-1}{K-1})=G(0)=0$ (11)

Thus, $S$ amplitude $\in[0,1]$ represents the degree of sparsity and stability
of the last layer activation, where 1 denotes a fully sparse and stable
response, and 0 indicates a dense and unstable activity. We can extend this
descriptor to every layer of the network $k$, to characterize the degree of
Sparsity, Stability, and Separability of all the centroids in multi-layer
networks. The only difference with layer K is that intermediate layers will
output events with spatial coordinates since they are not integrating over all
coordinates $\mathbf{x}$ (sec. 2). Therefore, the feedback time context
$FT_{i}^{k}$ will be defined as such:

$\mathbf{FT_{i}^{k}}=\max_{j\leq
i}\\{t_{j}^{k+1}|\mathbf{x_{i}^{k+1}},\mathbf{p_{j}^{k+1}}\\}$ (12)

meaning that $\mathbf{FT_{i}^{k}}$ pulls events with different polarities
$\mathbf{p_{j}^{k+1}}$ and the same spatial coordinates
$\mathbf{x_{i}^{k+1}}$.

##  4 Learning class-relevant features

Since we want to develop an online and local learning rule, we cannot resort
to backpropagation to maximize $S$, especially if we consider that HOTS
neurons are not differentiable. Therefore, we resort to a different approach,
focusing on using a modified online $k$-means learning rule. The learning rule
we propose is described by the following equation:

$\Delta\mathbf{c}^{k}_{f}=\alpha\Delta S^{k+1}_{i}\mathbf{q}^{k}_{i,f}+\beta
S^{k+1}_{i}\mathbf{q}^{k}_{i,f}$ (13)

where $\alpha$ and$\beta$ are learning rates, $\mathbf{c}^{k}_{f}$ represents
the position of the ”firing” centroid, assigned to a novel time surface
$\mathbf{ts}^{k}_{i}$ and
$\mathbf{q}^{k}_{i,f}=\mathbf{ts}^{k}_{i}-\mathbf{c}^{k}_{f}$. If we assume
$S^{k}_{i},\Delta S^{k}_{i}=1$ and $\eta=\alpha+\beta$ the equation takes the
form of online $k$-means [14]:

$\Delta\mathbf{c}^{k}_{f}=\eta\,\mathbf{q}^{k}_{i,f}=\eta\,(\mathbf{ts}^{k}_{i}-\mathbf{c}^{k}_{f})$
(14)

As shown (eq: 13), the proposed learning rule comprises two terms that weight
the online $k$-means equation (eq: 14) by $\Delta S$ and $S$. The first term
allows centroids to only move in the direction of a set of features when they
are associated with an increase in $S$. Conversely, centroids will move away
from features that decrease $S$ (such as features of the wrong class or class-
unrelated features), ensuring that only class-relevant features are learned.
As $S\rightarrow 1$, $\Delta S\rightarrow 0$ that makes the first term of
equation 13 tend to 0. This ensures the learning rate will decrease as we
approach a maximum, but it can cause the centroids never to reach the center
of mass of clusters. For this reason, the second term of the learning rule is
weighted by $S$ but with a different learning rate $\beta<\alpha$ to ensure it
is not adding instability to local maxima.

##  5 Individual thresholds

One problem arising from the presented learning rule is the choice of the
decision boundaries for the cluster assignments. In $k$-means, the decision
boundaries are defined as the equidistant points between two adjacent
clusters, similar to the borders of a Voronoi diagram. Inconveniently, the
learning rule presented with equation 13 is not compatible with this type of
decision boundary. The reason is that this type of segmentation will always
assign every sample to a centroid, making it impossible to learn only class-
relevant features. This means that when some centroids are pushed to more
relevant features, the less relevant features will be assigned to other nearby
clusters, reducing their specificity and ultimately lowering classification
accuracy. To avoid this scenario, we introduce individual thresholds
$th^{k}_{n}$, where n is the cluster index for layer $k$. We initially assign
every sample to the closest cluster as done in $k$-means, but now, every
cluster will generate a new event only when the sample distance to the cluster
is less than $th^{k}_{n}$. When this happens, we tune individual thresholds
with the equation:

$\Delta th^{k}_{f}=\gamma\Delta
S^{k+1}_{i}e^{-\frac{\lVert\mathbf{q}^{k}_{i,f}\rVert_{2}}{d^{k}}}+\delta
S^{k+1}_{i}e^{-\frac{\lVert\mathbf{q}^{k}_{i,f}\rVert_{2}}{d^{k}}}$ (15)

where $d$ is a parameter used to control how close to the border samples need
to be to have an effect on the threshold update. We also update nearby
thresholds with a competitive rule to reduce the number of events produced by
the network and threshold overlap between adjacent clusters. If a sample falls
within the decision boundaries of other ”non-firing” clusters, represented by
the index $nf$, and $\Delta S$ and $S$ are positive, we update their
thresholds with the following equation:

$\Delta th^{k}_{nf}=-\gamma\Delta
S^{k+1}_{i}e^{-\frac{\lVert\mathbf{q}^{k}_{i,nf}\rVert_{2}}{d^{k}}}-\delta
S^{k+1}_{i}e^{-\frac{\lVert\mathbf{q}^{k}_{i,nf}\rVert_{2}}{d^{k}}}$ (16)

Similarly to equation 13, both rule presents two terms: one weighted by
$\Delta S^{k+1}_{i}$ and one weighted by $S^{k+1}_{i}$. The first term changes
the decision boundaries to include samples associated with an increase of $S$
and reject samples associated with its decrease. In the situation where the
network is consistently wrong or consistently right, this first term will tend
to 0, requiring a second term weighted by $S^{k+1}_{i}$ to increase network
firing rate and stability for correct classes and reduce them for incorrect
ones. As noted for equation 13, $\delta<\gamma$ to ensure the stability of the
learning rule.

##  6 Network Initialization

Similarly to $k$-means, centroid initialization is an important step to ensure
good results [18]. Random initialization might put the centroids too far away
from the feature spaces or regions where the sum of updates from the learning
rule might push centroids even further from any meaningful features. However,
we are not interested in state-of-the-art initialization techniques since they
might obfuscate the actual performance results of Sup3r. Therefore, we
initialize centroids on the average time surfaces from a small batch of data
and add noise drawn from a uniform distribution:

$\mathbf{c_{j}^{k}}=(1-\zeta)\,\overline{\mathbf{ts^{k}}}+\zeta\,\mu_{ts}\,\mathbf{U_{j}^{k}}$
(17)

Were $\mathbf{U}$ represents a matrix with the same shape as
$\mathbf{c_{j}^{k}}$ and composed of random values pulled from a uniform
distribution. Since $\mathbf{ts_{i}}$ with small $\tau$ are mostly 0 while
$\mathbf{ts_{i}}$ with high $\tau$ are mostly 1, we scale $\mathbf{U}$ by the
average of the elements in $\overline{\mathbf{ts^{k}}}$, defined as
$\mu_{ts}$. This approach is extremely simple. In the future, more complex
approaches could be explored to ensure the convergence of most features and
boost learning speed and accuracy.

## Results

###  6.1 Synthetic benchmark

Hierarchy of Time-Surfaces uses unsupervised clustering algorithms such as
online k-means to extract recurring features in the data. While this approach
can be effective for pre-training on large quantities of unlabeled data, it
presents several drawbacks compared to state-of-the-art neuromorphic networks
that use approximate gradient methods based on Back Propagation Through Time.
One problem is the need for an external non-neuromorphic classifier, which we
outlined in the introduction. A second problem is the loss of accuracy given
by learning non-class-relevant features. To best represent this problem, we
design a tailored neuromorphic synthetic benchmark. We propose a visual
pattern classification task requiring a three-layered network (fig. 1). The
proposed task is to classify two ”sentences” (”v/v yty” and ”vxv yty”)
composed of two ”words” each made of three ”characters” drawn by lattice of
5x5 Poisson spiking neurons, for a total of 5x30 neurons (as can be seen in
fig. 1a). Yellow pixels correspond to ”high” firing rate while blue pixels
correspond to ”low” firing rate. We use Poisson neurons and background firing
rate to introduce variability per every sample. The training set consists of
1000 different recordings (500 for each sentence) of 10ms. This visual pattern
is then classified by a three-layered convolutional HOTS neural network (fig.
1b) where the first layer ($N^{1}=3$) responds to individual characters, the
second layer ($N^{2}=6$) responds to single words, and the last layer
($N^{3}=2$) classifies the entire sentence. It should be noted that the number
of HOTS centroids for each layer ($N_{i}$) does not correspond to the actual
number of clusters in the data. In the first layer, the number of clusters is
five (all the possible characters: ”v, /, x, y, t”), while in the second
layer, this number is three (all the possible words in the dataset: ”v/v, vxv,
yty”). This choice is motivated by the need to represent real-world scenarios
in which the number of clusters is unknown. This can be especially problematic
when the number of centroids is less than the actual number of clusters in the
data, as similar clusters representing two features of distinct classes might
be assigned a single centroid. This is the case of layer 1, making it likely
for the same centroid to represent similar characters like ”x” and ”/”. These
two characters are the only difference in the two sentences in the dataset
(fig 1a), making them class-relevant features. Failure to assign separate
centroids to ”x” and ”/” means that layer 1 will respond similarly to the two
sentences, causing the other network layers to fail to separate the two
sentences.

Figure 1: The visual classification benchmark to test Sup3r on a small
3-layered network. a) The test is composed of a ”sentence” (blue)
classification task, where two different sentences are each composed of two
among three different ”words” (green), which are composed of three of five
different ”characters” (orange). The only difference between the two sentences
consists of two characters, ”x” and ”/”, which are contained in the top word.
To generate spiking data, sentences are written by a lattice of 5x30 Poisson
neurons firing at a high firing rate for the yellow pixels and a low firing
rate for the blue pixels. b) The convolutional HOTS network used to solve the
task. The gray rectangles represent the output dimensionality of each layer.
$N$ is the number of clusters responding to the six different characters in a
sentence (Layer 1), to the two words in a sentence (Layer 2), and finally to
the entire sentence (Layer 3).

####  6.1.1 Classification Results

We show the classification results of a separate test set of 1000 recordings
(500 recordings for each sentence). We test Sup3r against a comparable network
with k-means features. Figure 2 shows the features from layer 1 for both
architectures. As expected, Sup3r can leverage Separability Sparsity and
Stability from higher layers to extract ”x” and ”/” (fig. 2a) without
utilizing the third centroid. Oppositely, k-means must assign every centroid
to a group of samples and minimize their variance, an optimization technique
that forces similar features to share a centroid when the number of centroids
is below the number of actual clusters in the dataset. In figure 2b, we define
the network accuracy by calculating the percentage of events assigned to the
right sentence in a single recording. The k-means algorithm causes HOTS to
fail to recognize the two sentences, resulting in $\approx$55.36% accuracy
(close to the chance level). Conversely, Sup3r reaches $\approx$99.92%
accuracy while using only two centroids. Since Sup3r can reject any event
produced by the non-class-relevant features, it outputs only $\approx$14.33%
of the input events (fig. 2c), potentially reducing energy consumption in
hardware implementations. To train the network we used the following
parameters: $\alpha=1\cdot 10^{-4}$, $\beta=1\cdot 10^{-5}$, $\gamma=1\cdot
10^{-4}$, $\delta=5\cdot 10^{-6}$, $\tau_{1}=1s$,$\tau_{2}=1ms$,
$\tau_{3}=1ms$ and ftv taus $f\tau_{2}=100ms$ and $f\tau_{3}=10ms$. We also
report the results in table 1.

Figure 2: Neuromorphic benchmark results, comparing Sup3r against $k$-means for feature extraction. a) Layer 1 features extracted with $k$-means (light blue) and Sup3r (orange). $k$-means fails to extract class-relevant features, while Sup3r can correctly extract ”x” ”/”, leaving one centroid unused. b) Consequently, $k$-means fails the classification task, performing close to the chance level ($\approx$55.36% accuracy, calculated as the percentage of events assigned to the correct class), while Sup3r can solve the task with $\approx$99.92% accuracy. Since Sup3r rejects events from non-class-relevant characters, only $\approx$14.33% of events are propagated in the network (calculated as the ratio between the number of the last layer output events divided by the number of the first layer input events). Conversely, the k-means network outputs the same number of input events. Table 1: Classification results (Mean on 10 networks trained and tested on independently generated synthetic datasets) | Accuracy | Processed events
---|---|---
Sup3r Network | 99.92% | 14.33%
$k$-means | 55.53% | 100%

####  6.1.2 Continual Learning

Sup3r learning signal $S_{i}^{K}$ is composed of a local signal (coming only
from the next layer ) and a global signal $G$. It is important to note that
the local signal is entirely unsupervised, as there is no need to know the
correct label to calculate it (eq. 7). This opens the possibility of keeping
the unsupervised learning running while inferring to adapt to potential shifts
in the class distributions. To test this feature, we design a test set of 1750
recordings by shifting ”x” and ”/” characters in the test sentences. To shift
the characters, we subtract pixels from “x” to add them to ”/” incrementally,
with the final shift causing the two letters to swap entirely (fig. 3a). Every
shift lasts for 250 recordings, with the last shift lasting for 1000 to ensure
convergence and stability. Looking at the Sup3r centroids, we can see them
adapting to the shift so that ”x” and ”/” become swapped (fig. 3b). To
evaluate Sup3r accuracy, in figure 3c, we test the network against an ablated
network with no unsupervised learning. As expected, Sup3r adapts to
distribution shifts while maintaining the original accuracy (>99%), while the
ablated network fails the task, ending with 0% accuracy. This final test was
performed with 10 independently trained networks with 10 distinct test sets.
To train the network, we used the same parameters we used for the
classification results.

Figure 3: Continual learning test with Sup3r. In this test, we perform a class
shift in the two sentences by subtracting pixels from ”x” and adding them to
”/” every 250 sentences (a), with the final shift lasting 1000 sentences.
Sup3r centroids slowly adapt to the shift (b). In (c), Sup3r network accuracy
remains high and returns to >99%, while the same network without unsupervised
learning fails the task (Ablation). This test is performed with ten
independently trained networks and randomized test sets. Accuracy is reported
as the average across runs in dark colors with min and max values in light
colors.

####  6.1.3 Incremental Learning

Our hypothesis is that Sup3r ability to extract sparse class-relevant features
makes it a natural candidate for Incremental Learning. Since we are extracting
a set of sparse class-relevant features, Sup3r might learn new tasks simply by
adding new centroids and fixing the previous ones to avoid catastrophic
forgetting.

To test this hypothesis, we design a new task (task 2) in which a Sup3r
network has to learn two novel sentences, "vpv yty" and "vov yty", after being
trained on the previous classification task (task 1) in section 6.1. We then
tested the same network on task 1 after retraining on task 2 to demonstrate
that it did not suffer catastrophic forgetting. We show the results in Table
2.

Table 2: Incremental Learning results (Mean on 10 networks trained and tested on independently generated synthetic datasets) | Task 1 | Task 2 | Task 1 after Task 2
---|---|---|---
Sup3r Network | 99.92% | 99.99% | 99.96%

The results show that the network can train on Task 2 and preserve the
accuracy previously reached on Task 1. To train the network, we used the same
parameters we used for the classification results.

##  7 Sup3r against backpropagation

We present preliminary results to show how Sup3r compares against
backpropagation. We use the N-MNIST dataset [19] a neuromorphic version of the
popular MNIST dataset [20]. This allows us to compare a Sup3r-HOTS network
with a traditional ANN. The Sup3r HOTS network has a single hidden layer with
32 centroids and spatial neighborhood $R^{1}=9$. The ANN has the same
architecture with a single convolutional hidden layer with 32 units with
sigmoid activation and filter size 9x9. We train the ANN with stochastic
gradient descent and a learning rate of 0.05. Other Sup3r parameters where
$\tau_{1}=100ms$, $\tau_{2}=1ms$ , $\alpha=5\cdot 10^{-3}$, $\beta=5\cdot
10^{-6}$, $\gamma=1\cdot 10^{-4}$, $\delta=1\cdot 10^{-7}$. The reason for
such low learning rates was the batch-based training implemented on GPU via
OpenCL to achieve higher parallelism and reduce computation time. In this
implementation, $\Delta c_{f}^{k}$ are summed for every event, causing large
centroid and threshold updates, which are compensated by lowering the learning
rate. Table 3 shows the result of this comparison.

Table 3: Sup3r NMNIST comparison. *HOTS implementation consisting of two hidden layers (32-96 units) of double exponential time surfaces and a SVM polynomial classifier [7]. Sup3r and ANN consisting of a single hidden layer with 32 units. | Sup3r on N-MNIST | ANN on MNIST | HOTS on N-MNIST[7],*
---|---|---|---
Test Accuracy | 94.98% | 95.61% | 91.27%
% Processed Events | 72.61% | N/A | 100%

Sup3r accuracy is comparable to an ANN on MNIST and higher than HOTS best
results [7] on a two-layered network (32-96 centroids) with double exponential
decay time surfaces and a polynomial SVM classifier, while cutting off 28$\%$
of events. While these results are promising, it is important to underline
that they are preliminary. The reason is that this test was implemented on
OpenCL for GPU acceleration prioritizing code readability and debugging rather
than hardware optimization. This limited the number of experiments we run on
N-MNIST. A single epoch of the N-MNIST dataset runs between 2 hours and 45
minutes on an NVIDIA 2060 Super, depending on the number of processed events
in the network, as computation is done event by event. This is due to the high
CPU usage to queue multiple kernels per event, the much bigger N-MNIST dataset
compared to his frame-based counterpart (several thousands of sparse events
compared to a single 28x28 frame), and to HOTS, which requires several small
matrix operations like multiplications and accumulation to process a single
event. Exploiting vendor-specific optimized kernels, sparse libraries, kernel
launching from GPU, and matrix operations should significantly reduce
computation time and allow for an exhaustive hyperparameter analysis and
testing on more complex architectures and datasets.

## Conclusion

In this paper, we presented Sup3r, a novel learning algorithm for training
Hierarchy of Time Surfaces algorithms. Our results show that this algorithm
can extract class-relevant features, reduce the number of processed events,
and enable incremental and continual learning without the need for external
classifiers. Moreover, preliminary results show that this algorithm could
reach comparable accuracy to backpropagation on a small network. All the
results together make Sup3r the best way to train HOTS networks end-to-end.
Future works will have to focus on accelerating Sup3r on GPUs, testing it on
more complex datasets and deep-layered architectures, and reducing the
excessive number of hyperparameters. Another point of interest could be to
adapt Sup3r to work on much more popular Spiking Neural Network models like
Integrate and Fire neurons.

## References

  * [1] Xavier Lagorce, Garrick Orchard, Francesco Galluppi, Bertram E Shi, and Ryad B Benosman.  Hots: a hierarchy of event-based time-surfaces for pattern recognition.  IEEE transactions on pattern analysis and machine intelligence, 39(7):1346–1359, 2016.
  * [2] Guillermo Gallego, Tobi Delbrück, Garrick Orchard, Chiara Bartolozzi, Brian Taba, Andrea Censi, Stefan Leutenegger, Andrew J Davison, Jörg Conradt, Kostas Daniilidis, et al.  Event-based vision: A survey.  IEEE transactions on pattern analysis and machine intelligence, 44(1):154–180, 2020.
  * [3] Daniel Gehrig, Antonio Loquercio, Konstantinos G Derpanis, and Davide Scaramuzza.  End-to-end learning of representations for asynchronous event-based data.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5633–5643, 2019.
  * [4] Gabriella Csurka, Christopher Dance, Lixin Fan, Jutta Willamowski, and Cédric Bray.  Visual categorization with bags of keypoints.  In Workshop on statistical learning in computer vision, ECCV, volume 1, pages 1–2. Prague, 2004.
  * [5] Ricardo Tapiador-Morales, Jean-Matthieu Maro, Angel Jimenez-Fernandez, Gabriel Jimenez-Moreno, Ryad Benosman, and Alejandro Linares-Barranco.  Event-based gesture recognition through a hierarchy of time-surfaces for fpga.  Sensors, 20(12):3404, 2020.
  * [6] Marco Rasetto, Juan P Dominguez-Morales, Angel Jimenez-Fernandez, and Ryad Benosman.  Event based time-vectors for auditory features extraction: a neuromorphic approach for low power audio recognition.  arXiv preprint arXiv:2112.07011, 2021.
  * [7] Marco Rasetto, Qingzhou Wan, Himanshu Akolkar, Feng Xiong, Bertram Shi, and Ryad Benosman.  Building time-surfaces by exploiting the complex volatility of an ecram memristor.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 2023.
  * [8] Amos Sironi, Manuele Brambilla, Nicolas Bourdis, Xavier Lagorce, and Ryad Benosman.  Hats: Histograms of averaged time surfaces for robust event-based object classification.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1731–1740, 2018.
  * [9] Antoine Grimaldi, Victor Boutin, Sio-Hoi Ieng, Ryad Benosman, and Laurent Perrinet.  A robust event-driven approach to always-on object recognition.  Authorea Preprints, 2023.
  * [10] Sumit B Shrestha and Garrick Orchard.  Slayer: Spike layer error reassignment in time.  Advances in neural information processing systems, 31, 2018.
  * [11] Adam Coates and Andrew Y. Ng.  Learning Feature Representations with K-Means, pages 561–580.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.
  * [12] Jelili Oyelade, Itunuoluwa Isewon, Funke Oladipupo, Olufemi Aromolaran, Efosa Uwoghiren, Faridah Ameh, Moses Achas, and Ezekiel Adebiyi.  Clustering algorithms: their application to gene expression data.  Bioinformatics and Biology insights, 10:BBI–S38316, 2016.
  * [13] Adam Jaeger and David Banks.  Cluster analysis: A modern statistical review.  Wiley Interdisciplinary Reviews: Computational Statistics, 15(3):e1597, 2023.
  * [14] Leon Bottou and Yoshua Bengio.  Convergence properties of the k-means algorithms.  Advances in neural information processing systems, 7, 1994.
  * [15] Mahardhika Pratama, Andri Ashfahani, and Edwin Lughofer.  Unsupervised continual learning via self-adaptive deep clustering approach.  In International Workshop on Continual Semi-Supervised Learning, pages 48–61. Springer, 2021.
  * [16] Christiaan Lamers, René Vidal, Nabil Belbachir, Niki van Stein, Thomas Bäeck, and Paris Giampouras.  Clustering-based domain-incremental learning.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 3384–3392, 2023.
  * [17] Teresa Serrano-Gotarredona and Bernabé Linares-Barranco.  A 128$\,\times$ 128 1.5% contrast sensitivity 0.9% fpn 3 µs latency 4 mw asynchronous frame-free dynamic vision sensor using transimpedance preamplifiers.  IEEE Journal of Solid-State Circuits, 48(3):827–838, 2013.
  * [18] Marina Gul and M Abdul Rehman.  Big data: an optimized approach for cluster initialization.  Journal of Big Data, 10(1):120, 2023.
  * [19] Garrick Orchard, Ajinkya Jayawant, Gregory K Cohen, and Nitish Thakor.  Converting static image datasets to spiking neuromorphic datasets using saccades.  Frontiers in neuroscience, 9:437, 2015.
  * [20] Li Deng.  The mnist database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine, 29(6):141–142, 2012.


# Information Reconciliation for High-Dimensional Quantum Key Distribution
using Nonbinary LDPC codes
††thanks: We acknowledge funding from: The Center of Excellence SPOC (ref
DNRF123)

Ronny Müller2, Davide Bacco1, Leif Katsou Oxenløwe2, and Søren Forchhammer2
2Department of Electrical and Photonics Engineering
Technical University of Denmark, Lyngby, Denmark
Email<EMAIL_ADDRESS>1Department of Physics and Astronomy
University of Florence, Florence, Italy

###### Abstract

Information Reconciliation is an essential part of Quantum Key distribution
protocols that closely resembles Slepian-Wolf coding. The application of
nonbinary LDPC codes in the Information Reconciliation stage of a high-
dimensional discrete-variable Quantum Key Distribution setup is proposed. We
model the quantum channel using a $q$-ary symmetric channel over which qudits
are sent. Node degree distributions optimized via density evolution for the
Quantum Key Distribution setting are presented, and we show that codes
constructed using these distributions allow for efficient reconciliation of
large-alphabet keys.

###### Index Terms:

Information reconciliation, High-dimensional QKD, nonbinary LDPC codes

##  I Introduction

Quantum Key Distribution (QKD) aims to share a symmetric secret key usable for
encryption between two remote parties. It consists of a quantum phase in which
the two parties (Alice and Bob) share some quantum information sent over a
quantum channel, followed by a classical post-processing phase in which a
shared secret key is extracted [1]. While QKD using quantum bits (qubits) is
well established, its High-dimensional version (HD-QKD) (using qudits) has
seen significantly less research effort. HD-QKD offers many potential
advantages, most importantly a significantly increased resilience to noise
[2]. Physical implementations of HD-QKD systems have been shown using a
variety of different technologies, including but not limited to time-bin
encoding [3], time-energy encoding [4], and orbital angular momentum [5].
An essential part of the post-processing is the information reconciliation
phase during which a common string is extracted from the results of the
quantum state measurements. This is followed by privacy amplification [6],
during which a secret key is extracted from that common string, e.g. using
universal hash functions. While a lot of research has analyzed and optimized
the information reconciliation phase for the two-dimensional case, little
regard has yet been paid to the high-dimensional case, apart from the
introduction of the layered scheme in 2013 [7]. In contrast to most channel
coding applications, the (HD)-QKD setting has lower requirements on latency
and throughput but puts a lot of value on a minimized information leakage.
Motivated by this setting, the excellent decoding performance of nonbinary
LDPC codes [8], and their native handling of high dimensions, we explore the
design and application of nonbinary LDPC codes for the post-processing in HD-
QKD protocols.

##  II Background

###  II-A Information Reconciliation

The information reconciliation stage in QKD aims to correct any mismatch
between the keys of the two parties while leaking as little information as
possible to a potential eavesdropper. In general, Alice sends a random string
$\mathbf{x}=(x_{0},...,x_{n-1})$, $x_{i}=\\{0,...,q-1\\}$ of $n$ qudits of dim
$q$ to Bob who proceeds to measure them, resulting in his version of the
string $\mathbf{y}=(y_{0},...,y_{n-1})$, $y_{i}=\\{0,...,q-1\\}$. We assume
that the quantum channel can be accurately modeled by a substitute channel
where $\mathbf{x}$ and $\mathbf{y}$ are correlated as a $q$-ary symmetric
channel, as errors are usually uncorrelated and symmetric. The transition
probabilities of such a channel is given by:

$\text{P}(y_{i}|x_{i})=\begin{cases}1-p&y_{i}=x_{i},\\\
\frac{p}{1-q}&\text{else}.\end{cases}$ (1)

Here, the parameter $p$ denotes the channel transition probability. We refer
to the symbol error rate between $\mathbf{x}$ and $\mathbf{y}$ as the quantum
bit error rate (QBER) in a slight abuse of notation but in accordance with
experimental works on HD-QKD. We assume the QBER to be an underlying channel
property in our simulations, making it equivalent to the channel parameter
$p$. Besides the qudits, Alice additionally sends a message $\mathbf{s}$
(usually the syndrome) of length $m$ over a classical communication channel
which is assumed to be error-free. From a coding point of view, this is equal
to asymmetric Slepian-Wolf coding with side information at the receiver, where
$\mathbf{s}$ represents the compressed version of $\mathbf{x}$ and
$\mathbf{y}$ is the side information. A more detailed description of this
equivalence can be found in [9]. Any information that is leaked to a potential
eavesdropper at any point of the quantum key distribution has to be subtracted
from the final secret key [10] during privacy amplification. The information
leaked during the information reconciliation stage will be denoted by
leak${}_{\text{IR}}$. Given no rate adaption, it can be upper-bounded by the
syndrome length in symbols, $\text{leak}_{\text{IR}}\leq m$. Using the
Slepian-Wolf bound [11], the minimum amount of leaked information to
successfully reconcile with arbitrary low failure chance is therefore given by
the conditional entropy:

$m\geq n\text{H}(X|Y).$ (2)

The conditional entropy of the $q$-ary symmetric channel, assuming independent
and identically distributed input $X$, can be written as

$\text{H}(X|Y)=-((1-p)\text{log}_{q}(1-p)-p\cdot\text{log}_{q}(\frac{p}{q-1})).$
(3)

The performance of a code with respect to information leakage can be measured
by its efficiency $f$, given by

$f=\frac{m}{n\text{H}(X|Y)}.$ (4)

Note that an efficiency of $f>1$ corresponds to leaking more bits than
required by the theoretical minimum of $f=1$, which represents the best
possible performance according to the Slepian-Wolf bound. Practical systems
have $f>1$.

TABLE I: The code rate, the best possible theoretical threshold (TT) according to the Slepian-Wolf bound, the ensemble threshold (ET) using density evolution, the ensemble efficiency (EE), the evaluated code efficiency (CE), and the respective degree polynomial for an 8-ary symmetric channel. Rate | TT | ET | EE | CE | Ensemble (edge view)
---|---|---|---|---|---
0.90 | 0.033 | 0.031 | 1.060 | 1.17 | $\lambda(x)=0.112x+0.103x^{2}+0.194x^{3}+0.146x^{9}+0.163x^{10}+0.003x^{17}$
|  |  |  |  |  $+0.173x^{19}+0.049x^{26}+0.006x^{28}+0.052x^{29}$
0.85 | 0.053 | 0.052 | 1.080 | 1.12 | $\lambda(x)=0.125x+0.165x^{2}+0.163x^{5}+0.11x^{7}+0.073x^{12}+0.089x^{18}$
|  |  |  |  |  $+0.122x^{27}+0.154x^{32}$
0.80 | 0.0758 | 0.0724 | 1.038 | 1.12 | $\lambda(x)=0.146x+0.177x^{2}+0.130x^{4}+0.084x^{7}+0.149x^{10}+0.035x^{19}$
|  |  |  |  |  $+0.029x^{22}+0.087x^{25}+0.163x^{26}$
0.75 | 0.100 | 0.096 | 1.030 | 1.10 | $\lambda(x)=0.165x+0.192x^{2}+0.092x^{5}+0.176x^{7}+0.019x^{10}+0.086x^{17}$
|  |  |  |  |  $+0.129x^{19}+0.103x^{30}+0.038x^{31}$
0.70 | 0.126 | 0.121 | 1.030 | 1.10 | $\lambda(x)=0.160x+0.208x^{2}+0.140x^{5}+0.096x^{8}+0.028x^{10}+0.013x^{11}$
|  |  |  |  |  $+0.113x^{18}+0.032x^{21}+0.211x^{27}$
0.65 | 0.154 | 0.147 | 1.032 | 1.10 | $\lambda(x)=0.173x+0.228x^{2}+0.092x^{4}+0.169x^{8}+0.112x^{14}+0.019x^{23}$
|  |  |  |  |  $+0.012x^{24}+0.195x^{28}$
0.60 | 0.183 | 0.177 | 1.026 | 1.10 | $\lambda(x)=0.192x+0.196x^{2}+0.222x^{5}+0.104x^{13}+0.114x^{23}+0.055x^{25}$
|  |  |  |  |  $+0.117x^{27}$
0.55 | 0.214 | 0.207 | 1.024 | 1.10 | $\lambda(x)=0.183x+0.269x^{2}+0.124x^{6}+0.036x^{8}+0.097x^{10}+0.004^{21}$
|  |  |  |  |  $+0.116x^{25}+0.171x^{26}$
0.50 | 0.247 | 0.239 | 1.024 | 1.11 | $\lambda(x)=0.215x+0.256x^{2}+0.030x^{4}+0.154x^{7}+0.065x^{11}+0.050^{13}$
|  |  |  |  |  $+0.072x^{21}+0.128x^{27}$
Figure 1: The frame error rate (FER) of 9 codes of coding rate [0.50,
0.55,…,0.90] (left to right) for different values of the QBER on a
$q=8$-dimensional $q$-ary symmetric channel. Each data point has been
evaluated with 5000 samples or until 100 errors have been collected.

###  II-B Nonbinary LDPC codes

We assume some familiarity with belief propagation decoding of binary Low-
density parity-check (LDPC) codes and Tanner graphs. For an overview see [12].
Nonbinary LDPC codes can be characterized by their parity check Matrix
$\mathbf{H}$, with $m$ rows and $n$ columns, which contain values in a Galois
Field (GF) of order $q$ as its entries. To increase clarity in this section,
we will mark all variables that represent a Galois field element with a hat,
e.g. $\hat{a}$. Furthermore, let $\oplus,\ominus,\otimes,$ and $\oslash$ refer
to the common operations on Galois field elements. An LDPC code can be
represented as a bipartite graph, the Tanner graph. There, the parity-check
equations constitute one side of the graph, called the check nodes, while the
symbols of the codeword represent the other side of the graph, the variable
nodes. The Tanner graph of a nonbinary LDPC code additionally has weighted
edges between the check and variable nodes where the weight corresponds to the
respective entry of $\mathbf{H}$. The syndrome $\mathbf{s}$ of the $q$-ary
string $\mathbf{x}$ is calculated as $\mathbf{s}=\mathbf{H}\mathbf{x}$.

In this work, we deploy a log-domain FFT-SPA [13, 14] for decoding. Detailed
descriptions of this algorithm can be found in [15, 16] but will be reiterated
here for completeness of the work. Let $Z$ denote a random variable that takes
values in GF($q$) such that $\text{P}(Z_{i}=k)$ denotes the probability that
the qudit $i$ has value $k=0,...,q-1$. The probability vector
$\mathbf{p}=(p_{0},...p_{q-1})$, $p_{j}=\text{P}(Z=j)$ can be transformed into
the log-domain using the generalized equivalent of the log-likelihood-ratio
(LLR) in the binary case, $\mathbf{m}=(m_{0},...,m_{q-1})$,
$m_{j}=\text{log}\frac{\text{P}(Z=0)}{\text{P}(Z=j)}=\log(\frac{p_{0}}{p_{j}})$.
Given the LLR presentation, one can go back to probabilities via
$p_{j}=\exp(-m_{j})/\sum_{k=0}^{q-1}\exp(-m_{k})$. We denote these transforms
with $p(\cdot)$ and $m(\cdot)$. To further simplify notation, we define the
multiplication and division of an element $\hat{a}$ in GF$(q)$ and an LLR
message as a permutation of the indices of the vector:

$\displaystyle\hat{a}\cdot\mathbf{m}$
$\displaystyle:=(m_{\hat{0}\oslash\hat{a}},...,m_{\hat{q-1}\oslash\hat{a}})$
(5) $\displaystyle\mathbf{m}/\hat{a}$
$\displaystyle:=(m_{\hat{0}\otimes\hat{a}},...,m_{\hat{q-1}\otimes\hat{a}}),$
(6)

where the multiplication and division of the indices take place in the Galois
Field. These permutations are required as we have to weigh messages according
to their edge weight during the decoding. We further define two
transformations involved with the decoding,

$\displaystyle\bar{\mathcal{F}}(\mathbf{m},\hat{H}_{ij})$
$\displaystyle=\mathcal{F}(p(\hat{H}_{ij}\cdot\mathbf{m}))$ (7)
$\displaystyle\bar{\mathcal{F}}(\mathbf{m},\hat{H}_{ij})^{-1}$
$\displaystyle=m(\mathcal{F}^{-1}(\mathbf{m}))/\hat{H}_{ij},$ (8)

where $\mathcal{F}$ is the discrete Fourier transform. Note that for $q$ being
a power of 2, the Fast Walsh Hadamard Transform can be used. The decoding
process then consists of two iterating message-passing phases, from check
nodes to variable nodes and vice versa. The message update rule at iteration
$l$ for the check node corresponding to the parity check matrix entry at
$(i,j)$ can then be written as

$\mathbf{m}^{(l)}_{ij,\text{CV}}=\mathcal{A}(\hat{s}^{\prime}_{i})\bar{\mathcal{F}}^{-1}(\underset{j\prime\in\mathcal{M}(i)/j}{\Pi}\bar{\mathcal{F}}(\mathbf{m}^{(l-1)}_{ij\prime},\hat{H}_{ij\prime}),\hat{H}_{ij}),$
(9)

where $\mathcal{M}(i)$ is the set of all check nodes in row $i$ of
$\mathbf{H}$. $\mathcal{A}$, defined as
$\mathcal{A}_{kj}(\hat{a})=\delta(\hat{a}\oplus k\ominus j)-\delta(a\ominus
j)$, accounts for the nonzero syndrome [16]. The weighted syndrome value is
calculated as $\hat{s}^{\prime}_{i}=\hat{s}_{i}\oslash\hat{H}_{ij}$. The a
posteriori message of column $j$ can be written as

$\tilde{\mathbf{m}}^{(l)}_{j}=\mathbf{m}^{(0)}(j)+\underset{i^{\prime}\in\mathcal{N}(j)}{\sum}\mathbf{m}^{(l)}_{i^{\prime}j,\text{CV}},$
(10)

where $\mathcal{N}(j)$ is the set of all check nodes in column $j$ of
$\mathbf{H}$. The best guess $\tilde{\mathbf{x}}$ at each iteration $l$ can be
calculated as the minimum value of the a posteriori,
$\tilde{x}_{j}^{(l)}=\text{argmin}(\tilde{m}_{j})$. The second message
passings, from variable to check nodes, are given by

$\mathbf{m}_{ij,\text{VC}}^{(l)}=\tilde{\mathbf{m}}^{(l)}_{j}-\mathbf{m}^{(l)}_{ij,\text{CV}}.$
(11)

The message passing is repeated until either
$\mathbf{H}\tilde{\mathbf{x}}=\mathbf{s}$ or the maximum number of iterations
is reached.

###  II-C Density Evolution

For a uniform edge weight distribution, the asymptotic decoding performance of
LDPC codes in the regime of infinite code length is completely characterized
by two polynomials [17, 18]:

$\lambda(x)=\sum_{i=0}^{d_{\text{v,
max}}}\lambda_{i}x^{i-1}\quad\rho(x)=\sum_{i=0}^{d_{\text{c,
max}}}\rho_{i}x^{i-1}.$ (12)

Here, $\lambda_{i}$ ($\rho_{i}$) are the fraction of edges connected to
variable (check) nodes of degree $i$ while $d_{\text{v, max}}$ ($d_{\text{c,
max}}$) is the highest degree of the variable (check) nodes. One can then
define the code ensemble $\mathcal{E}(\lambda,\rho)$, representing all codes
of infinite length with degree distributions given by $\lambda$ and $\rho$.
The threshold $p_{t}(\lambda,\rho)$ of the code ensemble
$\mathcal{E}(\lambda,\rho)$ is then defined as the worst channel parameter
(QBER) for which decoding is still possible with an arbitrary small failure
probability. This threshold can be approximated using Monte-Carlo Density
Evolution (MC-DE), which is described in detail in [19]. This method
repeatedly samples the node degrees according to $\lambda$ and $\rho$, and
draws random connections between nodes for each iteration. Given a high enough
sample size, this simulates the performance of a cycle-free code. Note that
MC-DE is especially well suited for nonbinary LDPC codes as the different edge
weights help in decorrelating messages [19]. During the simulation, we monitor
the mean entropy of all messages. Once it drops below a certain value, the
decoding is declared a success. If this does not happen after a maximum number
of iterations, the evaluated channel parameter is above the threshold of
$\mathcal{E}(\lambda,\rho)$. Using a concentrated check node distribution
(which is favorable according to [20]) and a fixed code rate, we can further
simplify to $\mathcal{E}(\lambda)$. One can then use the threshold as an
objective function to optimize the code design, which is commonly done using
the Differential Evolution algorithm [21].

##  III Results

While the described approach is feasible for any dimension $q$, we chose to
analyze for an 8-dimensional system as this is in the range of reported
implementations [22]. We designed 9 codes with code rates between 0.50 and
0.90, corresponding to a QBER range between 0 and 24.7%. We used 100000 nodes
with a maximum of 150 iterations for the MC-DE and swept the QBER in 20 steps
in a short range below the best possible threshold. For the Differential
Evolution, we used population sizes between 15 and 50, a differential weight
of 0.85, and a crossover probability of 0.7. We enforced a sparsity of at most
10 different degrees in the polynomial and choose $d_{\text{v, max}}=40$. The
sparsity allowed for reasonable optimization complexity. We observed that
higher $d_{\text{v, max}}$ can lead to numerical instability during decoding.
The results of the optimization can be found in Tab. I. For the optimization
and evaluation, we used the all-zero codeword assumption, which holds for the
given scenario of a symmetric channel [15]. For all rates, the designed
thresholds are close to the theoretical bound. We constructed LDPC codes with
a length of $n=30000$ symbols using Progressive Edge Growth [23] to confirm
the results of the density evolution experimentally. A log-FFT-SPA decoder is
used to reconcile the messages. The performance of the finite-size codes can
be seen in Fig. 1. As expected for finite-size codes, they do not reach the
ensemble threshold but show sub-optimal performance. To estimate the impact on
the secret key rate, we also evaluate the efficiency $f$ of the ensembles and
the constructed codes at their working points. We measured the efficiency at a
frame error rate (FER) of 1% which is a common heuristic. Note that there is a
direct trade-off between efficiency and FER which can be optimized to maximize
the overall secret key generation. After the reconciliation, the correctness
of the frames is asserted by exchanging short hashes. Incorrect frames are
then simply discarded, only lowering the secret key rate. QKD systems can
therefore efficiently operate at significantly higher FERs compared to what
would be acceptable in a classical communication setting. The resulting
efficiencies can be seen in Tab. I. The decoder runs a maximum of 100 decoding
iterations.

##  IV Discussion

Nonbinary LDPC codes pose as a natural candidate for the information
reconciliation stage of HD-QKD as their order can be matched to the dimension
of the used qudits, and nonbinary codes are known to have good decoding
performance [8]. Their main disadvantage, high decoding complexity, is less of
an obstacle in this setting as the keys can be processed and stored before
their usage in real-time applications, reducing the importance of decoding
latency. We expect HD-QKD systems to be used in a long-haul setting with a
high error rate and small count rate resulting in low throughput, further
reducing the requirements on throughput volume. Nevertheless, less complex
decoder algorithms like EMS [24] or TEMS [25] can be considered to allow the
usage of longer codes. This is regarded beneficial as they reduce finite-size
penalties of the secret key rate. Further, our codes allow for straightforward
rate adaption and integration into existing schemes developed for binary
channels, like Blind Reconciliation [26].

While the efficiencies we report are between 1.10 and 1.17, and therefore in a
similar range to the efficiencies of binary LDPC codes for conventional two-
dimensional QKD [27], it is difficult to make statements about their
performance in comparison to alternative schemes. To the best of our
knowledge, only the layered scheme [7], which handles the decoding using
$\lceil\log_{2}(q)\rceil$ binary decoders on bit planes, has been introduced
for information reconciliation with regard to HD-QKD. It is similar in concept
to the multilevel coding and multistage decoding methods used in slice
reconciliation for continuous-variable (CV) QKD [28]. While the layered scheme
allows for reconciliation using binary LDPC codes only, it brings its own
drawbacks, like error propagation, bit mapping, and interactive communication.
Unfortunately, little to no results that would allow for a fair comparison
have been published yet. Simulation results for $q=32$, a code length per
layer of 800 symbols, and a QBER of 20%, show an efficiency of the layered
scheme of $f\approx 1.4$ (Fig. 5 in [7]). Later experimental implementations
report efficiencies of 1.25 [29] ($q=3$, $n=1944$, $p=8\%$) and 1.17 [30]
($q=1024$, $n=4000$, $p=39.6\%)$. These papers report their efficiencies in
the $\beta$-notation. $\beta$ is commonly used in the continuous-variable QKD
community, whereas $f$ is more widespread with respect to discrete-variable
QKD. They can be related via
$\beta(\text{H}(X)-\text{H}(X|Y))=\text{H}(X)-f\text{H}(X|Y)$. Note that this
conversion is dependent on $p$ and $q$. We want to emphasize that a direct
comparison to the efficiencies reported in our work cannot be considered fair,
as the settings are not identical. A fair comparison between the two
approaches that additionally considers practicality remains the object of
future work.

While the found ensembles show thresholds close to the Slepian-Wolf bound, we
believe that even better results could be achieved by an extended search of
the hyperparameters involved in the optimization, e.g. the enforced sparsity
and highest degree of $\lambda$, and a finer sweep of the QBER during density
evolution. Similarly, we believe that the performance of the constructed codes
can be increased by using longer codes and improved versions of the PEG
algorithm [31, 32].

##  V Conclusion

In this paper, information reconciliation for High-Dimensional Quantum Key
Distribution using nonbinary LDPC codes over the $q$-ary symmetric channel is
introduced. We present optimized node degree distributions for the case of
$q=8$ that reach close to the Slepian-Wolf bound. Codes constructed according
to these distributions with a length of 30000 symbols reach efficiencies of
$f\approx 1.10$.

## References

  * [1] C. H. Bennett and G. Brassard, “Quantum cryptography: Public key distribution and coin tossing,” _Theoretical Computer Science_ , vol. 560, pp. 7–11, dec 2014. [Online]. Available: https://doi.org/10.1016%2Fj.tcs.2014.05.025
  * [2] H. Bechmann-Pasquinucci and W. Tittel, “Quantum cryptography using larger alphabets,” _Phys. Rev. A_ , vol. 61, p. 062308, May 2000. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevA.61.062308
  * [3] I. Vagniluca, B. Da Lio, D. Rusca, D. Cozzolino, Y. Ding, H. Zbinden, A. Zavatta, L. K. Oxenløwe, and D. Bacco, “Efficient time-bin encoding for practical high-dimensional quantum key distribution,” _Phys. Rev. Appl._ , vol. 14, p. 014051, Jul 2020. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevApplied.14.014051
  * [4] T. Zhong, H. Zhou, R. Horansky, C. Lee, V. Verma, A. Lita, A. Restelli, J. Bienfang, R. Mirin, T. Gerrits, S. Nam, F. Marsili, M. Shaw, Z. Zhang, L. Wang, D. Englund, G. Wornell, J. Shapiro, and F. Wong, “Photon-efficient quantum key distribution using time–energy entanglement with high-dimensional encoding,” _New Journal of Physics_ , vol. 17, p. 022002, 02 2015.
  * [5] D. Cozzolino, D. Bacco, B. Da Lio, K. Ingerslev, Y. Ding, K. Dalgaard, P. Kristensen, M. Galili, K. Rottwitt, S. Ramachandran, and L. K. Oxenløwe, “Orbital angular momentum states enabling fiber-based high-dimensional quantum communication,” _Phys. Rev. Appl._ , vol. 11, p. 064058, Jun 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevApplied.11.064058
  * [6] C. H. Bennett, G. Brassard, and J.-M. Robert, “Privacy amplification by public discussion,” _SIAM Journal on Computing_ , vol. 17, no. 2, pp. 210–229, 1988\. [Online]. Available: https://doi.org/10.1137/0217014
  * [7] H. Zhou, L. Wang, and G. Wornell, “Layered schemes for large-alphabet secret key distribution,” in _2013 Information Theory and Applications Workshop (ITA)_ , 2013, pp. 1–10.
  * [8] M. Davey and D. MacKay, “Low density parity check codes over gf(q),” in _1998 Information Theory Workshop (Cat. No.98EX131)_ , 1998, pp. 70–71.
  * [9] A. Liveris, Z. Xiong, and C. Georghiades, “Compression of binary sources with side information at the decoder using ldpc codes,” _IEEE Communications Letters_ , vol. 6, no. 10, pp. 440–442, 2002.
  * [10] K. Brádler, M. Mirhosseini, R. Fickler, A. Broadbent, and R. Boyd, “Finite-key security analysis for multilevel quantum key distribution,” _New Journal of Physics_ , vol. 18, no. 7, p. 073030, jul 2016. [Online]. Available: https://doi.org/10.1088%2F1367-2630%2F18%2F7%2F073030
  * [11] D. Slepian and J. Wolf, “Noiseless coding of correlated information sources,” _IEEE Transactions on Information Theory_ , vol. 19, no. 4, pp. 471–480, 1973\.
  * [12] S. Lin and J. Li, _Fundamentals of Classical and Modern Error-Correcting Codes_. Cambridge University Press, 2022\.
  * [13] S. Aruna and M. Anbuselvi, “Fft-spa based non-binary ldpc decoder for ieee 802.11n standard,” in _2013 International Conference on Communication and Signal Processing_ , 2013, pp. 566–569.
  * [14] H. Wymeersch, H. Steendam, and M. Moeneclaey, “Log-domain decoding of ldpc codes over gf(q),” in _2004 IEEE International Conference on Communications (IEEE Cat. No.04CH37577)_ , vol. 2, 2004, pp. 772–776 Vol.2.
  * [15] G. Li, I. J. Fair, and W. A. Krzymien, “Density evolution for nonbinary ldpc codes under gaussian approximation,” _IEEE Transactions on Information Theory_ , vol. 55, no. 3, pp. 997–1015, 2009.
  * [16] E. Dupraz, V. Savin, and M. Kieffer, “Density evolution for the design of non-binary low density parity check codes for slepian-wolf coding,” _IEEE Transactions on Communications_ , vol. 63, no. 1, pp. 25–36, 2015.
  * [17] V. Savin, “Non binary ldpc codes over the binary erasure channel: Density evolution analysis,” in _2008 First International Symposium on Applied Sciences on Biomedical and Communication Technologies_ , 2008, pp. 1–5.
  * [18] T. Richardson and R. Urbanke, “The capacity of low-density parity-check codes under message-passing decoding,” _IEEE Transactions on Information Theory_ , vol. 47, no. 2, pp. 599–618, 2001.
  * [19] M. Gorgoglione, V. Savin, and D. Declercq, “Optimized puncturing distributions for irregular non-binary ldpc codes,” in _2010 International Symposium On Information Theory & Its Applications_, 11 2010, pp. 400 – 405.
  * [20] S.-Y. Chung, T. Richardson, and R. Urbanke, “Analysis of sum-product decoding of low-density parity-check codes using a gaussian approximation,” _IEEE Transactions on Information Theory_ , vol. 47, no. 2, pp. 657–670, 2001\.
  * [21] R. Storn and K. Price, “Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces,” _Journal of Global Optimization_ , vol. 11, pp. 341–359, 01 1997.
  * [22] X.-M. Hu, C. Zhang, Y. Guo, F.-X. Wang, W.-B. Xing, C.-X. Huang, B.-H. Liu, Y.-F. Huang, C.-F. Li, G.-C. Guo, X. Gao, M. Pivoluska, and M. Huber, “Pathways for entanglement-based quantum communication in the face of high noise,” _Phys. Rev. Lett._ , vol. 127, p. 110505, Sep 2021. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.127.110505
  * [23] X.-Y. Hu, E. Eleftheriou, and D. Arnold, “Regular and irregular progressive edge-growth tanner graphs,” _IEEE Transactions on Information Theory_ , vol. 51, no. 1, pp. 386–398, 2005.
  * [24] D. Declercq and M. Fossorier, “Decoding algorithms for nonbinary ldpc codes over gf$(q)$,” _IEEE Transactions on Communications_ , vol. 55, no. 4, pp. 633–643, 2007.
  * [25] E. Li, K. Gunnam, and D. Declercq, “Trellis based extended min-sum for decoding nonbinary ldpc codes,” in _2011 8th International Symposium on Wireless Communication Systems_ , 2011, pp. 46–50.
  * [26] J. Martinez-Mateo, D. Elkouss, and V. Martin, “Blind reconciliation,” _Quantum Info. Comput._ , vol. 12, no. 9–10, p. 791–812, sep 2012.
  * [27] E. Kiktenko, A. Trushechkin, C. Lim, Y. Kurochkin, and A. Fedorov, “Symmetric blind information reconciliation for quantum key distribution,” _Physical Review Applied_ , vol. 8, no. 4, oct 2017. [Online]. Available: https://doi.org/10.1103%2Fphysrevapplied.8.044017
  * [28] G. Van Assche, J. Cardinal, and N. Cerf, “Reconciliation of a quantum-distributed gaussian key,” _IEEE Transactions on Information Theory_ , vol. 50, no. 2, pp. 394–400, 2004.
  * [29] J. Liu, Z. Lin, D. Liu, X. Feng, F. Liu, K. Cui, Y. Huang, and W. Zhang, “High-dimensional quantum key distribution using energy-time entanglement over 242 km partially deployed fiber,” 2022. [Online]. Available: https://arxiv.org/abs/2212.02695
  * [30] T. Zhong, H. Zhou, R. D. Horansky, C. Lee, V. B. Verma, A. E. Lita, A. Restelli, J. C. Bienfang, R. P. Mirin, T. Gerrits, S. W. Nam, F. Marsili, M. D. Shaw, Z. Zhang, L. Wang, D. Englund, G. W. Wornell, J. H. Shapiro, and F. N. C. Wong, “Photon-efficient quantum key distribution using time-energy entanglement with high-dimensional encoding,” _New Journal of Physics_ , vol. 17, no. 2, p. 022002, Feb. 2015\.
  * [31] H. Xiao and A. Banihashemi, “Improved progressive-edge-growth (peg) construction of irregular ldpc codes,” _IEEE Communications Letters_ , vol. 8, no. 12, pp. 715–717, 2004.
  * [32] A. Venkiah, D. Declercq, and C. Poulliat, “Randomized progressive edge-growth (randpeg),” in _2008 5th International Symposium on Turbo Codes and Related Topics_ , 2008, pp. 283–287.