# Strictly local Union–Find

Tim Chan<EMAIL_ADDRESS>Department of Materials, University
of Oxford, Parks Road, Oxford OX1 3PH, United Kingdom  Simon C. Benjamin
<EMAIL_ADDRESS>Department of Materials, University of
Oxford, Parks Road, Oxford OX1 3PH, United Kingdom  Quantum Motion, 9 Sterling
Way, London N7 9HJ, United Kingdom

(May 29, 2023)

###### Abstract

Fault-tolerant quantum computing requires classical hardware to perform the
decoding necessary for error correction. The Union–Find decoder is one of the
best candidates for this. It has remarkably organic characteristics, involving
the growth and merger of data structures through nearest-neighbour steps; this
naturally suggests the possibility of realising Union–Find using a lattice of
very simple processors with strictly nearest-neighbour links. In this way the
computational load can be distributed with near-ideal parallelism. Here we
build on earlier work to show for the first time that this strict (rather than
partial) locality is practical, with a worst-case runtime $\mathcal{O}(d^{3})$
and mean runtime subquadratic in $d$ where $d$ is the surface code distance. A
novel parity-calculation scheme is employed, which can also simplify
previously proposed architectures. We compare our strictly local realisation
with one augmented by long-range links; while the latter is of course faster,
we note that local asynchronous logic could largely negate the difference.

##  I Introduction

In the fault-tolerant era of quantum computing, the quantum hardware must be
supported by classical decoders that infer the nature of errors ‘on the fly’
from measurements. It is challenging to find decoder implementations that are
sufficiently fast, compact, and (ideally) with low power requirements. We
focus on decoders for the _surface code_ [1, 2, 3]: one of the most promising
error-correcting codes needed for fault tolerance due to its simplicity. Many
decoders have been developed such as _minimum-weight perfect matching_ (MWPM)
[4, 5], _renormalisation group_ [6, 7], _Markov chain Monte Carlo_ [8], and
various using _belief propagation_ [9, 10, 11, 12] and _neural networks_ [13,
14]. A recent development is the _hierarchical decoder_ [15, 16, 17, 18] in
which a simple, partially capable decoder handles most of the decoding,
sending the trickier problems to a more complex decoder.

One simple and fast but fully capable decoder is the Union–Find decoder (UF)
designed by Delfosse et al. [19, 20]. It has a relatively high accuracy and a
mean runtime slightly higher than cubic in $d$ [21, §2.3]. Liyanage et al.
[21] recently proposed _Helios_ , the FPGA implementation of an almost-local
version of UF whose output is identical to that of original UF. By _local_ we
mean runnable on a grid of identical nodes (classical processors), each
communicating only with their nearest neighbours. Helios is local save for an
additional controller which communicates with each node directly. They report
an improved mean runtime: sublinear in $d$.

Extending this literature, and especially the ideas proposed for Helios, we
introduce a paradigm which reduces the memory requirements of each node, the
number and size of messages passed around, the complexity of the grid
architecture, and the number of stages in the algorithm as a whole. This
improvement is a simple scheme inspired by anyon annihilation to calculate
parities. Moreover we also design a _strictly local_ version of UF (SLUF), the
first of its kind, whose mean runtime is subquadratic in $d$ for error
probabilities below threshold. Having compared these strictly local and
almost-local models we note that the use of a fast communication relay,
potentially implemented with asynchronous logic, would constitute a third
version of UF which rivals the speed of almost-local UF whilst maintaining
strict locality.

The paper proceeds as follows. In Section II we cover the prerequisite
background theory. In Section III we discuss almost-local UF as introduced in
the Helios scheme, and describe our variation of this idea in the almost-local
paradigm, before describing SLUF in detail in Section IV. We evaluate runtime
for both versions of UF and discuss how SLUF can be sped up using asynchronous
logic in Section V before concluding in Section VI.

##  II Background

Section II.1 covers the surface code with a focus on decoding. We approach
this standard material using graph theory as it facilitates a natural
explanation of UF presented in Section II.2.

###  II.1 Surface Code and Decoding

As its name suggests, the surface code arranges the physical qubits on a 2D
grid. It encodes one logical qubit and corrects for both bit- and phaseflip
errors (the ability to correct these two discrete errors allows us to correct
any qubit error [22, §10.3.1]). Here, we describe how to correct bitfip
errors; phaseflip errors are corrected analogously. We explain the decoding
cycle assuming measurements to be perfect in Section II.1.1 before
generalising to faulty measurements in Section II.1.2.

####  II.1.1 Perfect Measurements

Figure 1:  (a) The graph representing a distance-5 surface code. Bulk nodes
are green; boundary nodes, blue. (b) An example of an error $\mathbb{E}$ is
the set of red edges whose corresponding syndrome $\mathbb{S}$ is the set of
red nodes. Note a path in $\mathbb{E}$ (like that from node 26 to 29) will
make a defect only at each endpoint; a cycle will make no defects. The decoder
sees only $\mathbb{S}$, not $\mathbb{E}$. (c) A correction $\mathbb{C}$ for
$\mathbb{E}$ is the set of red edges. The syndrome made by $\mathbb{C}$ is
always the same as that of $\mathbb{E}$. (d) The leftover $\mathbb{L}$ after
combining $\mathbb{E}$ and $\mathbb{C}$. This example comprises a cycle and
two paths: one between nodes on the same boundary and one between opposite
boundaries. The latter represents a logical bitflip.

The surface code is based on the graph $G=(V,E)$ in Fig. 1(a) which comprises
a set $V=V_{\text{bulk}}\cup V_{\text{boundary}}$ of nodes arranged on a
$d\times(d+1)$ grid and a set $E$ of edges. Members of $V_{\text{bulk}}$ are
called bulk nodes; $V_{\text{boundary}}$, boundary nodes. Each edge
corresponds to a set of two nodes from $V$. Physically, each bulk node
represents an ancilla qubit used for measurement (sometimes called a _measure
qubit_) while each edge represents a data qubit. Boundary nodes represent
nothing physical but exist just to allow every data qubit to be represented by
an edge. We will treat the decoding problem abstractly in terms of nodes and
edges. In a decoding cycle:

  1. 1.

Noise corrupts our system: Each edge bitflips with some probability $p$ called
the _physical error probability_. The _error_ $\mathbb{E}\subseteq E$ is the
set of bitflipped edges.

  2. 2.

Ancilla qubits flag the noise: Every bulk node records the parity of the edges
in $\mathbb{E}$ which are incident to it. If it records odd parity, it is a
_defect_. The _syndrome_ $\mathbb{S}\subseteq V_{\text{bulk}}$ is the set of
defects. Figure 1(b) shows an error and its syndrome. Mathematically,
$\mathbb{S}=\sigma(\mathbb{E})$ where

$\sigma(\mathbb{E}):=V_{\text{bulk}}\cap\operatorname*{{\triangle}}_{e\in\mathbb{E}}e$
(1)

and
$\operatorname*{{\triangle}}_{e\in\mathbb{E}}e:=f\mathbin{\vartriangle}\dots\mathbin{\vartriangle}g$
for $\mathbb{E}=\\{{f,\dots,g}\\}$ where $\mathbin{\vartriangle}$ denotes
symmetric difference.

  3. 3.

The classical decoder acts: Given $\mathbb{S}$, the decoder must output a
_correction_ $\mathbb{C}\subseteq E$ which, if treated as an error, would make
the same syndrome $\mathbb{S}$. In other words, $\mathbb{C}$ must satisfy
$\sigma(\mathbb{C})=\mathbb{S}$. Figure 1(c) shows a correction.

  4. 4.

Success or failure: The _leftover_ $\mathbb{L:=E\mathbin{\vartriangle}C}$ is
then guaranteed to satisfy the following Lemma (for completeness, a proof is
given in Appendix A).

######  Lemma 1.

$\mathbb{L}$ comprises only cycles or paths between distinct boundary nodes.

The question of whether the decoder has been successful depends on whether a
logical error has resulted from the correction. Any path between opposite
boundaries represents a logical bitflip. Hence, if $\mathbb{L}$ has an odd
number of such paths the logical qubit has picked up a logical error. Figure
1(d) shows a leftover with one such path.

The aim of the decoder is to output a correction which is least likely to lead
to a logical error. Each decoder has a different _threshold_ i.e. physical
error probability $p_{\textnormal{th}}$ below (above) which increasing $d$
decreases (increases) the logical error probability. Clearly a decoder is
effective only below its threshold i.e. $p<p_{\textnormal{th}}$. Moreover, the
further below threshold a decoder operates, the more effective it is.
Practically, it is useless to operate a decoder at its threshold; significant
benefits of error correction are reaped only when $p$ is several times below
threshold.

####  II.1.2 Faulty Measurements

Figure 2: The 3D generalisation of Fig. 1 to account for faulty measurements.
(a) Instead of distance-5, we show a distance-3 code for visual clarity. The
time axis points upward. (b) The error $\mathbb{E}$ can now include vertical
‘timelike’ edges like that from node 18 to 19. The syndrome $\mathbb{S}$ is
the set of red nodes throughout the whole lattice. (c) A correction
$\mathbb{C}$ for $\mathbb{E}$. (d) The leftover $\mathbb{L}$ comprises one
path between the same boundary and one between opposite boundaries. The latter
represents a logical bitflip.

In reality, ancilla qubit measurement is faulty i.e. in Item 2 each bulk node
records the _wrong_ parity with some probability $q$. The result is that the
syndrome is subject to some noise:
$\mathbb{S}=\sigma(\mathbb{E})\mathbin{\vartriangle}\mathbb{F}$, where
$\mathbb{F}\subseteq V_{\text{bulk}}$ is the set of nodes which record the
wrong parity for that measurement round. Note boundary nodes cannot be in
$\mathbb{F}$ as they do not represent ancilla qubits which are measured.

The usual way to account for faulty measurements is to measure each ancilla
qubit not once but $\tau+1$ times per decoding cycle, where $\tau>0$
(intuitively, the higher $q$ is, the higher we should make $\tau$). We then
have not one 2D sheet of syndrome data as in Fig. 1(b), but $\tau+1$ of them,
which we label $\\{{\mathbb{S}_{0},\dots,\mathbb{S}_{\tau}}\\}$. In general
there is a new $\mathbb{E}$ and $\mathbb{F}$ for each sheet, but the data is
cumulative in $\mathbb{E}$ i.e.

$\mathbb{S}_{t}=\big{[}{\operatorname*{{\triangle}}_{u=0}^{t}\sigma(\mathbb{E}_{u})}\big{]}\mathbin{\vartriangle}\mathbb{F}_{t}\qquad\forall
t\in 0..\tau$ (2)

so to localise each $\mathbb{E}_{t}$ we take the difference between
consecutive sheets:

$\displaystyle\textnormal{}\mathbb{S}_{t}$
$\displaystyle:=\mathbb{S}_{t}\mathbin{\vartriangle}\mathbb{S}_{t-1}\qquad\forall
t\in 1..\tau$
$\displaystyle\overset{\eqref{eq:2D_sheet}}{=}\sigma(\mathbb{E}_{t})\mathbin{\vartriangle}\mathbb{F}_{t}\mathbin{\vartriangle}\mathbb{F}_{t-1}.$
(3)

We can then stack these 2D ‘difference sheets’
$\\{{\textnormal{}\mathbb{S}_{1},\dots,\textnormal{}\mathbb{S}_{\tau}}\\}$ to
make a 3D simple cubic lattice as in Fig. 2(a). Sheets are stacked such that
$\textnormal{}\mathbb{S}_{1}$ is on the bottom, $\textnormal{}\mathbb{S}_{2}$
is the next one above, etc. This way, sheets which derive from later
measurements are higher up the stack so we can interpret the upward direction
as going forward in time (hence the subscript $t$).

Each bulk node no longer represents a single ancilla qubit; rather, the
difference between two consecutive measurements of that ancilla qubit at a
given point in time. Each horizontal edge no longer represents a single data
qubit but a possible time at which data qubit could bitflip. Vertical edges
are a new addtion: each one represents a single possible faulty measurement.
Specifically, each edge between $\textnormal{}\mathbb{S}_{t}$ and
$\textnormal{}\mathbb{S}_{t+1}$ corresponds to a possible member of
$\mathbb{F}_{t}$. This makes sense e.g. a member of $\mathbb{F}_{1}$ effects a
change in both $\textnormal{}\mathbb{S}_{1}$ and $\textnormal{}\mathbb{S}_{2}$
as per Eq. 3.

If we redefine $G$ as this 3D simple cubic lattice, all steps of the decoding
cycle are the same as in Section II.1 with the modification in Item 1 that
each _horizontal_ edge bitflips with probability $p$; each _vertical_ edge,
$q$. Figure 2(b–d) shows a decoding cycle example.

In our numerical emulations we make two simplifications:

  1. 1.

We assume the _global decoding_ scheme i.e. we set
$\mathbb{F}_{0}=\mathbb{F}_{\tau}=\varnothing$, which explains the absence of
vertical edges below $\textnormal{}\mathbb{S}_{1}$ and above
$\textnormal{}\mathbb{S}_{\tau}$ in Fig. 2. This allows us to determine the
success of the decoder for a batch of $\tau+1$ measurement rounds in
isolation. Considering how consecutive batches join at their time boundaries
and decoding them is the topic of _sliding window decoding_ which we do not
treat in this paper. Global decoding, though less realistic, captures the
relevant behaviour of decoders like MWPM and UF, and leads to a negligible
difference in decoding accuracy [23, Appendix C].

  2. 2.

We set $p=q$ which suggests $\tau$ should be on the same order as $d$, hence
we set $\tau=d$. Our versions of UF could account for $p\neq q$ by weighting
the edges of $G$ in the syndrome validation stage.

MWPM is a popular decoder which uses the following idea: for small $p$ the
number of edges in $\mathbb{E}$ is likely to be small, so it seems sensible to
minimise the number of edges in the correction. MWPM finds a correction whose
size is minimal. While accurate, it has a worst-case runtime
$\mathcal{O}(N^{3}\lg N)$ [24, §3] where $N:=|V|=\mathcal{O}(d^{2})$ for
perfect measurements and $\mathcal{O}(d^{3})$ for faulty measurements. UF can
be seen as an approximation of MWPM [25] with an improved worst-case runtime
$\mathcal{O}\big{(}{N\alpha(N)}\big{)}$ where $\alpha$, the inverse of
Ackermann’s function, is less than three for all practical surface code
distances.

For visual clarity, Fig. 3 and onwards use the graph $G$ for perfect
measurements. However, the reader should keep in mind all the material which
follows (including our versions of UF) also applies to $G$ for faulty
measurements.

###  II.2 UF

Conceptually, UF groups defects into spatial clusters then finds a correction
only from edges within a cluster. UF keeps the total size of all clusters
small to ensure the correction is small. The algorithm comprises the following
stages:

  1. 1.

Syndrome Validation An active cluster nucleates at each defect which grows
outward in all directions at the same speed.

  2. 2.

Spanning Forest Growth A spanning tree is grown within each cluster. In our
modification of UF, this stage comes for free.

  3. 3.

Peeling Each edge of each tree is peeled i.e. removed from the tree and
included or excluded from the correction depending on the syndrome.

####  II.2.1 Syndrome Validation

Figure 3: A syndrome validation example which needs 4 growth rounds. Ungrown
edges are invisible; half-grown, dotted; fully grown, solid. Edges which have
grown since the last round are black; else, grey. (a) Each defect grows its
incident edges by half. Note some edges have grown by half twice, so are fully
grown. (b) The two active clusters grow again. The west one touches a boundary
hence inactivates. (c) The remaining active cluster grows and touches hence
merges with the inactive cluster north of it. (d) The resultant active cluster
grows and touches the boundary hence inactivates.

Figure 3 shows a syndrome validation example. Mathematically, a cluster
$C=(V_{C},E_{C})$ is a connected subgraph of $G$. A cluster is _active_ iff
there exists no correction for its defects using only edges in the cluster
i.e. $\nexists\mathbb{C}\subseteq E_{C}:\sigma(\mathbb{C})=\mathbb{S}\cap
V_{C}$. Else it is _inactive_. Intuitively, an inactive cluster is one where
we have accounted for all defects with a correction local to that cluster; an
active cluster still seeks such a correction. In practice we determine
activity by the following lemma adapted from [20, Lemma 1] and proved in
Appendix B.

######  Lemma 2.

Cluster $C$ is active iff it has an odd number of defects and touches no
boundary i.e.

$\big{(}{\textnormal{{\hbox{|\mathbb{S} \cap V_C|}}
odd}}\big{)}\wedge\big{(}{V_{\textnormal{boundary}}\cap
V_{C}=\varnothing}\big{)}.$ (4)

Before syndrome validation, a cluster $(\\{{v}\\},\varnothing)$ is initialised
for every $v\in V$ (so there is an active cluster for each defect and an
inactive one for each nondefect). Each edge in $E$ has a _growth value_ in
$\\{{0,\frac{1}{2},1}\\}$ which is initialised at 0. Syndrome validation
proceeds via some number of growth rounds and stops when all clusters are
inactive. In a growth round:

  1. 1.

Each active cluster grows all edges around it by a half.

  2. 2.

Each new fully grown edge merges the two clusters at its endpoints into one
bigger cluster, if the endpoints did not already belong to the same cluster.
Mathematically, if clusters $A$ and $B$ merge along edge $e$ the resultant
cluster is

$C=\big{(}{V_{A}\cup V_{B},E_{A}\cup E_{B}\cup\\{{e}\\}}\big{)}.$ (5)

One may wonder why edges grow by only half per growth round. If they grew by 1
per round, a pair of defects separated by two edges would merge just as fast
as a pair separated by one edge. Growing edges by half per round distinguishes
these cases, reducing the number of fully grown edges hence the size of all
clusters.

Syndrome validation has a runtime $\mathcal{O}\big{(}{N\alpha(N)}\big{)}$ if
the Union–Find data structure (after which the Union–Find decoder is named) is
used [26]. It ensures the _existence_ of a correction using only cluster edges
while the next two stages, shown in Fig. 4, find an _instance_ of one.

Figure 4: Continuing the example in Fig. 3: (a) A spanning forest is grown in
the direction of the arrows. Each arrow is an edge in the forest and points
from parent to child node. (b) The forest is peeled in the direction of the
arrows. The set of red arrows is the correction returned by the decoder.

####  II.2.2 Spanning Forest Growth

For each cluster $C$ a root node is chosen as follows: if $C$ touches a
boundary, one of its boundary nodes is chosen; else, any node is chosen. A
spanning tree of $C$ (a subgraph of $C$ that is connected, acyclic, and spans
$V_{C}$) is then grown from this root. The union of all the trees is the
forest. This stage is done by breadth- or depth-first search; either way, the
runtime is linear in the number of edges in the forest i.e. $\mathcal{O}(N)$.

####  II.2.3 Peeling

A correction $\mathbb{C}=\varnothing$ is initialised then the forest is
traversed edge by edge, starting from the leaves. Each time a defect is
encountered all edges traversed thereafter, until another defect is
encountered, are added to $\mathbb{C}$. Once the whole forest has been
traversed $\mathbb{C}$ will satisfy $\sigma(\mathbb{C})=\mathbb{S}$ [19,
Theorem 1]. Since each edge in the forest is traversed exactly once the
runtime is $\mathcal{O}(N)$.

This concludes our review of the well-established surface code and UF
paradigms, using the notation that we require for the remainder of the paper.
We now briefly discuss recent literature working towards local architectures
for implementing UF, especially Helios.

##  III Almost-Local UF

As we can see the nature of UF is quite local: all edges involved connect
neighbouring nodes, and clusters communicate with each other only when they
touch; so it is natural to try to implement UF as locally as possible. This
may have advantages over the standard picture in which the decoder protocol is
implemented using a single, fully fledged classical processor, presumably
located remotely from the actual qubits. Instead, one can imagine a more
distributed decoder comprising a very simple classical processor corresponding
to each node and an interprocessor communication link corresponding to each
edge. We can visualise this local architecture as directly corresponding to
Fig. 2(a) if nodes represent processors and edges represent the links.

This alternative paradigm has various benefits. Firstly, its computation is
parallelised: non-neighbouring computations can be run independently and
simultaneously; we will see later this leads to superior runtime scaling.
Secondly, it is faster as each processor is simpler and sits closer to its
data source (the ancilla qubit). Thirdly, there may be implications in terms
of robustness: if one part breaks down then only that part, rather than the
whole decoder, must be replaced.

In Section III.1 we give a conceptual overview of Helios developed by Liyanage
et al. [21] and in Section III.2 we describe almost-local UF in detail with
our modification _anyon annihilation_ which simplifies the physical
architecture and allows a spanning forest to be grown in 1 timestep regardless
of $d$.

###  III.1 Helios

The authors consider only syndrome validation as this is the only stage that
is nontrivial to implement locally. Their implementation uses the local
architecture just described i.e. there is a processor for each node in $G$ and
a communication link for each edge. For the rest of this section we refer to
processors as nodes and communication links as edges.

Clusters are identified by a unique integer stored by all its nodes. When two
clusters merge the higher-integer cluster becomes part of the other one. This
happens via a flooding algorithm: each node looks at its neighbours along
fully grown edges; if it sees a lower integer than its own, it stores that
lower integer instead. This results in a wave of change propagating across the
higher-integer cluster, starting from where the two clusters touched. After
the flood, all the nodes of the higher-integer cluster will store the lower
integer. Generalising to three or more clusters merging at once is simple: a
flood propagates from each touching point but eventually only the lowest-
integer cluster will remain.

As well as merging, the other task in syndrome validation is determining the
activity of each cluster. The authors devised the following messaging system
to tackle this. Each cluster always has a well-defined root node. Whenever a
defect changes its stored integer, it sends a message to the root of the
cluster it has just joined. Every root stores a bit and flips it whenever it
receives a message. Once all flooding has finished and all messages have
reached their destination, the bit of each root will equal the parity of the
number of defects in its cluster. This, by 2, determines the activity of
clusters which do not touch a boundary. The authors do not mention treatment
of clusters which touch a boundary but we do so in Section III.2.

These messages are sent via a network of communication links comprising the
edges of $G$ and additional edges between non-neighbouring nodes [21, Figure
3] but our modification removes the need for them.

In addition to the nodes, there is a single controller which manages the
global variables of the UF algorithm. It connects to each node directly (in
general it connects to each node indirectly via a hierarchical tree structure
but for low enough code distances this tree is not used) but the type of
information communicable is restricted: the controller can broadcast the same
message to all nodes and receive Boolean messages from any node (but does not
know which node sent the message). This controller–node connectivity is the
reason we say Helios is almost local.

The worst-case runtime of Helios is $\mathcal{O}(N)$ because the maximum
possible cluster length is $\mathcal{O}(N)$ (see example cluster in Fig. 5)
and the duration to flood such a cluster scales linearly with its length. This
is formally better than the $\mathcal{O}\big{(}{N\alpha(N)}\big{)}$ of
original UF owing to the parallelism of Helios. More impressive is the _mean_
runtime which, for faulty measurements i.e. $N=\mathcal{O}(d^{3})$, the
authors report as sublinear in $d$ from numerics.

Figure 5: A maximum-length cluster is a winding snake-like path of length
$\mathcal{O}(N)$ where $N$ is the number of nodes in the graph $G$.

###  III.2 Our Version of Almost-Local UF

Figure 6: Flowchart for the stage variable which alternates between fixed- and
arbitrary-duration stages.

We assume decoding to be cellular automata-like i.e. via discrete timesteps
with each component performing one operation, and messages travelling one
edge, per timestep.

The controller has a variable attribute `stage` indicating which stage it is
in. It is initialised to `growing` and follows the flowchart in Fig. 6. There
are seven stages: `growing`, `merging`, `presyncing`, `syncing`, `burning`,
`peeling`, `done`. The odd stages last exactly 1 timestep; the even stages, an
arbitrary number of timesteps. During syndrome validation the controller
cycles through the first four stages once per growth round. `burning` is for
spanning forest growth, `peeling` is for peeling, and `done` indicates the end
of the decoding cycle.

Each node is assigned a unique integer `ID` as in Fig. 2(a) and has the
following variable attributes:

  * •

`defect` is a boolean indicating whether the node is a defect.

  * •

`active` is a boolean indicating whether the cluster the node is in is active.
Since each node starts in its own cluster, $v.\verb|active|$ is initialised as
$v.\verb|defect|~{}\forall v\in V$. If $v.\verb|active|=\verb|true|$ we say
$v$ is active.

  * •

`CID` is an integer indicating the cluster the node belongs to. Clusters are
identified by the lowest `ID` of all its nodes. The node with this `ID` is the
_root_ of the cluster. $v.\verb|CID|$ is initialised as $v.\verb|ID|~{}\forall
v\in V$.

  * •

`anyon` is a boolean indicating whether the node has an anyon. This anyon can
be thought of as a particle which is passed between nodes during syndrome
validation. When two anyons meet i.e. a node receives two anyons in a
timestep, they annihilate. $v.\verb|anyon|$ is initialised as
$v.\verb|defect|~{}\forall v\in V$. The idea is to accumulate at the root all
the anyons in a cluster so they annihilate on the way or on arrival.
Eventually zero or one anyon will remain at the root, indicating the parity of
the number of defects in the cluster.

  * •

`pointer` is the direction the node should relay the anyon. Possible values
(assuming faulty measurements) are `C`, `N`, `W`, `E`, `S`, `D`, `U`,
representing centre, north, west, east, south, down, up, respectively.
Following the pointers starting from any node should eventually lead to the
root of the cluster it is in. Only roots have a `pointer` of `C` as they do
not relay anyons but accumulate them. $v.\verb|pointer|$ is initialised to
$\verb|C|~{}\forall v\in V$. Together, `anyon` and `pointer` replace Helios’
messaging system.

  * •

`busy` is a boolean used to indicate whether the node is busy. $v.\verb|busy|$
is initialised to $\verb|false|~{}\forall v\in V$.

Each edge has a variable attribute `support` equal to its growth value, so
$e.\verb|support|$ is initialised as $0~{}\forall e\in E$. Since processing is
done at nodes rather than edges we can for each edge pick one of its
endpoints, say the one with the lower `ID`, to be responsible for storing its
`support`. The set of edges a node is responsible for, is then the _owned
edges_ :

$E_{\text{owned}}(v):=\\{{uv:uv\in
E~{}\text{and}~{}v.\texttt{ID}<u.\texttt{ID}}\\}.$ (6)

In each timestep the controller runs Algorithm 1, which simply moves `stage`
on if no node is busy. This is where the controller–node connectivity is used:
the controller need not know _which_ node is busy; rather, if _a_ node is
busy. After every `syncing` stage the controller also checks for any active
clusters. If there are, another growth round is needed so `stage` is reset to
`growing`. If none, syndrome validation is done so `stage` is set to
`burning`. Once `stage` reaches `done` the controller stops running Algorithm
1 and instead stores the correction $\mathbb{C}$ made during peeling.

Algorithm 1 Run by controller every timestep.

AdvanceAlmostLocalController $(\text{not}~{}v.\texttt{busy})~{}\forall v\in
V$$\texttt{stage}=\texttt{syncing}$$(\text{not}~{}v.\texttt{active})~{}\forall
v\in V$

$\texttt{stage}\leftarrow\texttt{burning}$

$\texttt{stage}\leftarrow\texttt{growing}$

$\texttt{stage}\leftarrow$ next stage $\triangleright$ Stage order is as
introduced in Section III.2.

####  III.2.1 Syndrome Validation

Figure 7: An example of a merging stage which needs 3 timesteps. For each
node, ID is shown in the top-left if $\texttt{ID}=\texttt{CID}$; else, CID is
shown in the top-right. If the node has an anyon, a 1 is shown in the bottom-
left. Pointers are shown in the bottom-right if they differ from C. (a) The
first growth stage has just finished. (b) Nodes 7 and 8 change their CID to 6
and point to node 6 which is now the root of a 3-node cluster. (c) Nodes 7 and
8 relay their anyons to the root; the two anyons annihilate and one remains,
indicating the cluster has an odd number of defects. This cluster touches no
boundary so by 2 is active.

It is helpful to define the _accessibles_ of a node as the set of neighbours
along fully grown edges:

$\operatorname{accessibles}(v):=\\{{u\in V:uv\in E\wedge
uv.\texttt{support}=1}\\}.$ (7)

In a timestep each node runs whichever procedure matches `stage` out of the
four in Algorithm 2.

Algorithm 2 The four procedures each node can run during syndrome validation.

Growing$v$$v.\texttt{active}$unfully grown edges $e$incident to $v$

$e.\texttt{support}\leftarrow e.\texttt{support}+\frac{1}{2}$

Merging$v$

$v.\texttt{busy}\leftarrow\texttt{false}$ ($v\in V_{\text{bulk}}$) and ($v$
not a root) and ($v.\texttt{anyon}$)

$v.\texttt{busy}\leftarrow\texttt{true}$

relay anyon in direction of $v.\texttt{pointer}$
$u\in\operatorname{accessibles}(v)$ $u.\texttt{CID}<v.\texttt{CID}$

$v.\texttt{busy}\leftarrow\texttt{true}$

$v.\texttt{pointer}\leftarrow$ toward $u$

$v.\texttt{CID}\leftarrow u.\texttt{CID}$

Presyncing$v$

$v.\texttt{active}\leftarrow$ ($v\in V_{\text{bulk}}$) and ($v$ a root) and
($v.\texttt{anyon}$)

Syncing$v$

$v.\texttt{busy}\leftarrow$ (not $v.\texttt{active}$) and ($\exists
u\in\operatorname{accessibles}(v):u.\texttt{active}$)

$v.\texttt{active}\leftarrow v.\texttt{active}$ or $v.\texttt{busy}$

Stage `growing` is simply a change of `support`s: each node in an active
cluster increments the `support` of edges around it, if not already fully
grown.

Figure 7 shows a `merging` example. In `merging` two processes happen:
flooding and anyon annihilation. Flooding is as described in Section III.1:
each node looks at its accessibles; if it sees a lower `CID` than its own, it
updates its own to match it. Additionally, the node updates its `pointer`
toward that accessible. This ensures following the pointers always leads to
the correct root promised by `CID`. In anyon annihilation, any node with an
anyon relays it in the direction of its pointer. At the end of this stage, a
cluster has an odd number of defects iff its root has an anyon.

Stage `presyncing` determines which clusters are active. Note in Fig. 2(a)
boundary nodes have lower `ID`s than bulk nodes. This ensures that if a
cluster touches a boundary, its root is a boundary node. Each root can thus
determine if its cluster is active by checking if itself is a boundary node
and if it has an anyon, by 2. In `presyncing`, roots of active clusters set
their `active` attribute to `true`; every other node sets theirs to `false`.

Figure 8: An example of a syncing stage which needs 2 timesteps. Shown
attributes are same as in Fig. 7 except the bottom-right which, instead of the
pointer, shows ‘a’ if the node is active. (a) A presyncing stage has just
finished. Note how the only active node is the root of the 3-node active
cluster. (b) Nodes 7 and 8 see this active node along a fully grown edge so
become active themselves. This is the end of the flood as all the nodes in the
active cluster are active.

Figure 8 shows a `syncing` example. In `syncing` the activity of each cluster
is propagated from root to all its other nodes via a flood: each node looks at
its accessibles; if it sees an active node, it becomes active, if not already.
At the end of this stage, a node is in an active cluster iff it is active.

In Helios’ messaging system, messages contain the `ID` of the destination
node. The memory needed to store this scales logarithmically in $d$. A new
message is emitted whenever a defect changes its `CID` which leads to
$\mathcal{O}(N^{2})$ messages throughout syndrome validation in the worst
case. Messages must be kept separate so buffers are needed but these can stall
if full. Anyon annihilation improves upon this system as `anyon` is a single
bit regardless of $d$. Throughout syndrome validation, the number of anyons
never exceeds the number of defects in the graph $G$ and only ever decreases,
so is $\mathcal{O}(N)$ in the worst case. There is no need for buffers (as
anyons need not be kept separate) nor additional links between non-
neighbouring nodes.

####  III.2.2 Burning

Figure 9: An example of burning and peeling. Pointers are shown as in Fig. 7
but in the middle of each node. (a) Before burning, a cluster contains a
cycle, within which there is one edge no pointer points along. (b) This edge
is burned i.e. removed to leave behind a spanning forest. (c) The four leaves
of the forest are peeled resulting in a smaller forest. (d) The two leaves of
the smaller forest are peeled resulting in an even smaller forest. These two
edges are added to the correction as the leaf nodes were defects. Edges so far
added to the correction are in red. (e) The final leaf is peeled and added to
the correction.

In original UF a spanning forest is grown in $\mathcal{O}(N)$ timesteps via
breadth- or depth-first search. Our version does this in 1 timestep by
retrieving the forest directly from the `pointer`s in a stage we call
_burning_ , described in Algorithm 3. An example is shown in Fig. 9(a, b). In
burning, each fully grown edge looks at its two endpoints; if neither of their
`pointer`s point along the edge, that edge is ‘burned’ i.e. it sets its
`support` to zero. The spanning forest is the set of fully grown edges which
remain. This is because the pointers already define a spanning tree within
each cluster; burning simply removes the edges ‘unused’ by these pointers.

Algorithm 3 Run by each node during burning.

Burning$v$$uv\in E_{\text{owned}}(v)$$uv.\texttt{support}=1$and not
($u.\texttt{pointer}$toward $v$or $v.\texttt{pointer}$toward $u$)

$uv.\texttt{support}\leftarrow 0$

####  III.2.3 Peeling

Figure 9(b–e) shows our version of peeling which is simply a local version of
that described in Section II.2.3 and takes $\mathcal{O}(N)$ timesteps. In a
timestep each node runs the procedure in Algorithm 4: each node checks if it
is a leaf node of the spanning forest; if so it will peel its corresponding
leaf edge.

Algorithm 4 Run by each node during peeling.

Peeling$v$

$v.\texttt{busy}\leftarrow\texttt{false}$ ($v$ not a root) and
$\big{(}{|{\operatorname{accessibles}(v)}|=1}\big{)}$

2:$v.\texttt{busy}\leftarrow\texttt{true}$

$\\{{u}\\}:=\operatorname{accessibles}(v)$

4:$uv.\texttt{support}\leftarrow 0$ $v.\texttt{defect}$

$v.\texttt{defect}\leftarrow\texttt{false}$

6:add $uv$ to $\mathbb{C}$

flip $u.\texttt{defect}$

####  III.2.4 Worst-Case Runtime Scaling

Our version’s worst-case decoding runtime is $\mathcal{O}(N)$ due to the same
reason as Helios’: flooding a cluster of length $\mathcal{O}(N)$ is the
process with the worst scaling in the whole decoding cycle. Specifically, the
number of timesteps to flood a cluster equals its length in edges. We expect
mean runtime to scale almost identically to Helios’ which we indeed observe in
Section V.

##  IV Strictly Local UF

We now present a version of UF without a direct communication link between the
controller and each node. Instead, the controller communicates only with the
node whose `ID` is 0. Broadcasting of `stage` from controller to nodes happens
via a flood through the whole graph $G$, and `busy` and `active` propagate as
signals from node to controller by relaying. We call these two new processes
_staging_ and _signalling_ , respectively. Signalling is similar to how anyons
are relayed but instead of being relayed along pointers these signals are
simply relayed toward the controller, one edge per timestep. This requires the
controller and nodes to store additional attributes which are described in
Section IV.1.

The change clearly increases runtime but has the advantage of being strictly
local, hence much easier to physically implement. The worst-case runtime
remains $\mathcal{O}(N)$ and is discussed in Section IV.3. We will see in
Section V that the mean runtime stays below $\mathcal{O}(d^{2})$ for physical
error probabilities below threshold (recall from Section II.1.1 this is the
decoder’s operable region). Moreover we presently note that a rapid-
transmission variant can substantially recover the speed of the almost-local
architecture within the strictly local paradigm.

###  IV.1 Additional Attributes

The controller stores a constant integer `span` which equals the number of
edges to the furthest node. This node is always the boundary node with the
highest `ID` e.g. node 17 in Fig. 2(a). The controller also has the following
variable attributes:

  * •

`countdown` is an integer which tracks how long the controller must wait until
it is sure staging or signalling is done. It is initialised to 0.

  * •

`busy_signal` is a boolean indicating whether the controller receives a
`busy_signal` in the current timestep.

  * •

`active_signal` is a boolean indicating whether the controller _has_ received
an `active_signal` during the current `syncing` _stage_. It is initialised to
`false`.

Each node stores the following constants:

  * •

`span` is an integer which equals the number of edges from itself to the node
furthest from the controller.

  * •

`signalee` is the neighbour in the direction toward the controller. This is
the neighbour the node looks to for staging, and sends `busy_signal`s and
`active_signal`s to. For faulty measurements: it equals the neighbour below;
if the node is already in the lowest sheet $\textnormal{}\mathbb{S}_{1}$, it
is instead the east neighbour; if the node is already eastmost, it is instead
the north neighbour. For node 0, `signalee` is the controller.

Each node also has the following variable attributes:

  * •

`stage` is same as that of controller. This is used for staging and is
initialised to $\verb|growing|~{}\forall v\in V$.

  * •

`countdown` is same as that of controller but used only to track how long the
node must wait until staging is done. It is initialised to $0~{}\forall v\in
V$.

  * •

`busy_signal` is a boolean indicating if in the current timestep the node
either is busy or receives a `busy_signal` from another node.

  * •

`active_signal` is a boolean indicating if in the current timestep the node
either has `stage` as `presyncing` and becomes active or receives an
`active_signal` from another node.

To manage these additional attributes the procedures have additional steps.
These are explained in the next subsection.

###  IV.2 Modified Procedures

For the controller, Algorithm 1 is replaced by Algorithm 5. Main differences
are: checking for any busy nodes is replaced by checking if `countdown` has
reached zero; checking for any active nodes, by checking if an `active_signal`
has been received.

Algorithm 5 Run by controller every timestep in strictly local UF.

AdvanceLocalController $\texttt{countdown}=0$stagechange block.
$\texttt{stage}=\texttt{syncing}$not active_signal

$\texttt{stage}\leftarrow\texttt{burning}$

2:$\texttt{stage}\leftarrow\texttt{growing}$

$\texttt{stage}\leftarrow$ next stage countdown reset block.
$\texttt{stage}\in\\{{\texttt{merging},\texttt{peeling}}\\}$

4:$\texttt{countdown}\leftarrow\texttt{span}+1$

$\texttt{countdown}\leftarrow\texttt{span}$ busy_signal and
$\texttt{countdown}\in\\{{1,2}\\}$

6:$\texttt{countdown}\leftarrow 2$ $\triangleright$ Only happens during
arbitrary-duration stages.

$\texttt{countdown}\leftarrow\texttt{countdown}-1$

In Section III.2 the four stages that lasted 1 timestep each follow an
arbitrary-duration stage, so in SLUF, staging must be used to broadcast these
four stages to the nodes. As indicated in 5 it takes `span` timesteps for
staging to complete i.e. for a `stage` flood to propagate from controller to
the furthest node. Hence these four stages are still of fixed duration.

The other three (arbitrary-duration) stages need not be broadcast as they each
follow a fixed-duration stage so the nodes know exactly when to start them. In
`merging` and `peeling` the controller must wait at least $\verb|span|+1$
timesteps as it takes `span` timesteps for a signal to go from the furthest
node to the controller. 4 indicates this. In `syncing` the furthest node (a
boundary node) is never busy nor active so the controller need only wait at
least `span` timesteps (5). In any of these three stages, after this minimum
duration, `busy_signal`s are no more than 2 edges apart due to 1 so
`countdown` resets to 2 in 6.

######  Proposition 1 (Doppler effect).

For any `busy_signal` received by the controller after $\verb|span|+1$
timesteps in `merging`, another `busy_signal` must have been received $\leq 2$
timesteps earlier.

###### Proof.

All `busy_signal`s emitted in the first `merging` timestep reach the
controller in $\verb|span|+1$ timesteps or less. Hence any `busy_signal` $s$
received after this time must have been emitted by some node $v$ at some time
$t_{s}$ _after_ the first `merging` timestep. This implies at $t_{s}-1$ some
node $u\in\\{{v}\\}\cup\operatorname{accessibles}(v)$ was busy, so emitted a
`busy_signal` $s^{\prime}$. The duration between the controller receiving
$s^{\prime}$ and $s$ is maximised when $u$ is one edge closer to the
controller than $v$. In this case $s^{\prime}$ arrives 2 timesteps before $s$.
∎

The proof of 1 can be understood as a Doppler effect: waves of busyness travel
through clusters at a speed of 1 edge per timestep. Wavefronts emit
`busy_signal`s at a frequency of 1 per timestep. If a wave travels away from
the controller, the controller will receive these signals at a frequency of 1
per 2 timesteps.

In a timestep each node runs AdvanceLocalNode in Algorithm 6. This reuses
Algorithms 2, 3 and 4 but wraps each procedure in either AdvanceFixed or
AdvanceArbitrary depending on whether the stage is of fixed or arbitrary
duration.

Algorithm 6 Each node in strictly local UF runs AdvanceLocalNode every
timestep.

AdvanceLocalNode$v$

$\Call{Proc}{}\leftarrow$ procedure matching $v.\texttt{stage}$ in Algorithms
2, 3 and 4
$v.\texttt{stage}\in\\{{\texttt{merging},\texttt{syncing},\texttt{peeling}}\\}$

2:AdvanceArbitrary$v$, Proc

AdvanceFixed$v$, Proc

4:AdvanceFixed$v$, Proc $\triangleright$ For fixed-duration stages.
$v.\texttt{countdown}=0$

Proc$v$ $v.\texttt{stage}=\texttt{presyncing}$

6:$v.\texttt{active\\_signal}\leftarrow v.\texttt{active}$

$v.\texttt{stage}\leftarrow$ next stage

8:$\texttt{countdown}\leftarrow\texttt{countdown}-1$

AdvanceArbitrary$v$, Proc $\triangleright$ For arbitrary-duration stages.
$v.\texttt{stage}=v.\texttt{signalee}.\texttt{stage}$

10:Proc$v$

$v.\texttt{busy\\_signal}\leftarrow v.\texttt{busy\\_signal}$ or
$v.\texttt{busy}$

12:relay $v.\texttt{busy\\_signal}$ and $v.\texttt{active\\_signal}$ to
$v.\texttt{signalee}$

$v.\texttt{stage}\leftarrow v.\texttt{signalee}.\texttt{stage}$
$\triangleright$ Staging.

14:$v.\texttt{countdown}\leftarrow v.\texttt{span}$

A node’s `countdown` is used only when its `stage` is of fixed duration,
during which it should equal the controller’s `countdown`. It tells the node
exactly when to change `stage`: as soon as it reaches zero in 4. At other
times i.e. after arbitrary-duration stages, the update of `stage` happens via
staging (13).

###  IV.3 Worst-Case Runtime Scaling

The new processes in SLUF, staging and signalling, take $\mathcal{O}(d)$
timesteps and are used $\mathcal{O}(1)$ times per growth round. A decoding
cycle needs a maximum of $\mathcal{O}(d)$ growth rounds so these processes
contribute a maximum $\mathcal{O}(d^{2})$ timesteps to the overall runtime.
This scaling does not exceed $\mathcal{O}(N)$ regardless of whether
$N=\mathcal{O}(d^{2})$ or $\mathcal{O}(d^{3})$ so the worst-case runtime
scaling of SLUF remains $\mathcal{O}(N)$.

##  V Runtime Data

Figure 10: Mean runtime to complete syndrome validation against code distance
for four different physical error probabilities $p$ (we set measurement error
probability $q=p$). Each datapoint is the mean of $10^{3}$ samples of decoding
cycle, with errorbars showing standard error. (a) Almost-local UF. (b) Same as
(a) but on a log–log scale. (c) Strictly local UF. (d) Same as (c) but on a
log–log scale.

We test the runtime of our version of almost-local UF (ALUF) and SLUF by
emulating many decoding cycles and counting the number of timesteps to
complete syndrome validation. This is a good proxy for decoding runtime as we
already know the exact runtime of burning, and the runtime of peeling equals
the length of the longest cluster which is at most the runtime of the last
round of merging.

In fact, `burning` and `peeling` need not be standalone stages. We could have
incorporated peeling into merging if, instead of relaying an anyon along an
edge, we peeled that edge, with the modifications in Algorithm 4 that 4 is not
done and 6 becomes
$\mathbb{C}\leftarrow\mathbb{C}\mathbin{\vartriangle}\\{{uv}\\}$. Defects
would then take over the role of anyons and the decoding cycle would finish
after the last `syncing` stage, without need for burning. Although small, this
modification may be useful in future implementations for which we want to
minimise the number of stages e.g. a sliding window local UF.

Figure 10 shows the mean runtime for both versions of UF. The gradients of the
best-fit lines for the log–log plots (b, d) estimate the mean runtime scaling
with code distance. These are tabulated in Table 1. Note $p=2.6\cdot 10^{-2}$
is the threshold of UF so $5.0\cdot 10^{-3}$ is roughly the highest practical
$p$ at which UF would operate. For this value of $p$ and below: ALUF shows
sublinear scaling; SLUF, subquadratic scaling. We expect the former as ALUF is
closely related to Helios which Liyanage et al. [21] observe to have sublinear
scaling.

Table 1: The best estimate of $m$ such that the mean runtime $\propto d^{m}$ where $d$ is code distance. $p$ | almost-local UF | strictly local UF
---|---|---
$5.0\cdot 10^{-4}$ | 0.35(2) | 1.09(2)
$1.0\cdot 10^{-3}$ | 0.48(3) | 1.21(4)
$5.0\cdot 10^{-3}$ | 0.83(1) | 1.39(3)
$2.6\cdot 10^{-2}$ | 1.61(1) | 1.63(3)

As well as the mean, we also evaluate the distribution of runtimes which are
shown in Figs. 11 and 12. SLUF has equally spaced sharp peaks in runtime; each
one corresponds to a certain number of growth rounds. The wide spacing implies
most of the runtime is due to staging and signalling i.e. the controller
waiting until countdown reaches zero before moving on [this can also be
inferred by the difference in absolute runtime between ALUF (which does not
need a countdown) and SLUF: at $d=25$ the former is more than 20 times faster
than the latter]. The countdown start, hence the peak spacing, increases
linearly with $d$.

The peak locations can be derived as follows. If syndrome validation comprises
0 or 1 growth rounds, `countdown` resets to `span` or $\verb|span|+1$ three
times (see Fig. 6). For 2 growth rounds, `countdown` resets seven times. For
$r\geq 0$ growth rounds, `countdown` resets $\max\\{{3,4r-1}\\}$ times. Hence
the peak corresponding to $r$ growth rounds is roughly located at
$\verb|span|\cdot\max\\{{3,4r-1}\\}$ timesteps.

Figure 11: Distribution of runtime to complete syndrome validation using
almost-local UF, for physical error probability $p=5.0\cdot 10^{-3}$. Within
each histogram: the horizontal magenta line shows the mean; red line, the
maximum, of $10^{3}$ samples of decoding cycle. Each bar has width 1. Figure
12: Same as Fig. 11 but for strictly local UF. Each bar has width 21. Note the
$x$-axes show counts on a log scale.

###  V.1 Asynchronous Logic

As just mentioned, most of the runtime of SLUF is due to staging and
signalling (defined in Section IV). The fraction of runtime spent on these
operations for the realistic scenario of $p=5.0\cdot 10^{-3}$ and $d=25$, is
0.955. This motivates the question: can we recover the speed of ALUF whilst
maintaining strict locality? So far we have only considered versions of UF
using _synchronous_ logic i.e. where each node changes its state once per
clock cycle, based on states in previous clock cycles. Signals between nodes
therefore are limited to travel at one edge per clock cycle. _Asynchronous_
logic on the other hand, has no clock cycle. Instead, each node changes its
state as soon as it receives the information to do so without having to wait
until the next clock cycle.

This has two advantages. Firstly, less heat is generated as no clock pulse
need be distributed to every flip-flop on the circuit. Secondly, it is faster:
signal speed is limited only by how fast each node responds to its inputs.
This speed advantage is greatest if used for simple, indiscriminate operations
like staging and signalling.

Therefore, it makes sense to consider a strictly local version of UF which
uses asynchronous logic for staging and signalling, and synchronous logic for
all other operations. In the case that it takes less than one synchronous
clock cycle for a signal to travel across the whole array via asynchronous
logic, the runtime will be indentical to ALUF whilst maintaining strict
locality.

We stress that within our SLUF paradigm, staging and signalling are ideally
suited to implementation with asynchronous logic since there is no risk of
race conditions leading to indeterminate behaviour. The signals travel through
the array in such a fashion that it does not matter if propagation along one
path is faster or slower than an alternate; we need only understand the
timescale for the process to be complete.

##  VI Conclusion

We explore UF realised by local architectures, following recent developments
including the Helios implementation [21]. We present an improved ‘almost-
local’ model, and moreover a new implementation that is strictly local; the
first of its kind. We numerically reproduce the sublinear in $d$ mean runtime
scaling of Helios then observe a subquadratic scaling for our strictly local
variant of UF. Both scalings are better than the slightly higher-than-cubic
scaling of original UF. Further, strictly local UF offers architecture design
advantages due to its local nature. Lastly, we note our paradigm is compatible
with the use of asynchronous logic to massively speed up runtime without
sacrificing locality, making it an attractive design option.

##  VII Acknowledgements

We thank Armands Strikis and David Garner for useful discussions. The authors
would like to acknowledge the use of the University of Oxford Advanced
Research Computing (ARC) facility [27] in carrying out this work and
specifically the facilities made available from the EPSRC QCS Hub grant
(agreement No. EP/T001062/1). The authors also acknowledge support from
EPSRC’s Robust and Reliable Quantum Computing (RoaRQ) project (EP/W032635/1).
TC acknowledges an EPSRC DTP studentship.

##  VIII Author Contributions

SCB contributed to the core concepts for ALUF/SLUF, and authored elements of
the paper. TC established the core concepts ALUF/SLUF, implemented, tested and
validated these concepts, wrote the simulation code, collected and processed
the numerical data, and was the primary author of the paper.

Note added: At the time of submission of this preprint, we noted that Ref.
[21] had recently been updated to incorporate a concept comparable to our
pointer mechanism. The revised version of [21] notes the advantages versus its
earlier messaging system, which we allude to in our discussions of this work.

## References

  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Topological quantum memory, Journal of Mathematical Physics 43, 4452 (2002).
  * Fowler  _et al._ [2012a] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Surface codes: Towards practical large-scale quantum computation, Physical Review A 86, 032324 (2012a).
  * Litinski [2019] D. Litinski, A game of surface codes: Large-scale quantum computing with lattice surgery, Quantum 3, 128 (2019).
  * Edmonds [1965] J. Edmonds, Paths, trees, and flowers, Canadian Journal of Mathematics 17, 449 (1965).
  * Fowler  _et al._ [2012b] A. G. Fowler, A. C. Whiteside, and L. C. L. Hollenberg, Towards practical classical processing for the surface code, Physical Review Letters 108, 180501 (2012b).
  * Duclos-Cianci and Poulin [2010a] G. Duclos-Cianci and D. Poulin, Fast decoders for topological quantum codes, Physical Reviw Letters 104, 050504 (2010a).
  * Duclos-Cianci and Poulin [2010b] G. Duclos-Cianci and D. Poulin, A renormalization group decoding algorithm for topological quantum codes, in  _2010 IEEE Information Theory Workshop_ (2010) pp. 1–5.
  * Wootton and Loss [2012] J. R. Wootton and D. Loss, High threshold error correction for the surface code, Physical Review Letters 109, 160503 (2012).
  * Criger and Ashraf [2018] B. Criger and I. Ashraf, Multi-path summation for decoding 2D topological codes, Quantum 2, 102 (2018).
  * Higgott  _et al._ [2022] O. Higgott, T. C. Bohdanowicz, A. Kubica, S. T. Flammia, and E. T. Campbell, Fragile boundaries of tailored surface codes and improved decoding of circuit-level noise (2022), arXiv:2203.04948 [quant-ph] .
  * Higgott and Breuckmann [2022] O. Higgott and N. P. Breuckmann, Improved single-shot decoding of higher dimensional hypergraph product codes (2022).
  * Kuo and Lai [2022] K.-Y. Kuo and C.-Y. Lai, Exploiting degeneracy in belief propagation decoding of quantum codes, npj Quantum Information 8, 111 (2022).
  * Sheth  _et al._ [2020] M. Sheth, S. Z. Jafarzadeh, and V. Gheorghiu, Neural ensemble decoding for topological quantum error-correcting codes, Physical Review A 101, 032338 (2020).
  * Overwater  _et al._ [2022] R. W. J. Overwater, M. Babaie, and F. Sebastiano, Neural-network decoders for quantum error correction using surface codes: A space exploration of the hardware cost-performance tradeoffs, IEEE Transactions on Quantum Engineering 3, 1 (2022).
  * Delfosse [2020] N. Delfosse, Hierarchical decoding to reduce hardware requirements for quantum computing (2020), arXiv:2001.11427 [quant-ph] .
  * Meinerz  _et al._ [2022] K. Meinerz, C.-Y. Park, and S. Trebst, Scalable neural decoder for topological surface codes, Physical Review Letters 128, 080505 (2022).
  * Ravi  _et al._ [2022] G. S. Ravi, J. M. Baker, A. Fayyazi, S. F. Lin, A. Javadi-Abhari, M. Pedram, and F. T. Chong, Better than worst-case decoding for quantum error correction (2022), arXiv:2208.08547 [quant-ph] .
  * Smith  _et al._ [2023] S. C. Smith, B. J. Brown, and S. D. Bartlett, Local predecoder to reduce the bandwidth and latency of quantum error correction, Physical Review Applied 19, 034050 (2023).
  * Delfosse and Zémor [2020] N. Delfosse and G. Zémor, Linear-time maximum likelihood decoding of surface codes over the quantum erasure channel, Physical Review Research 2, 033042 (2020).
  * Delfosse and Nickerson [2021] N. Delfosse and N. H. Nickerson, Almost-linear time decoding algorithm for topological codes, Quantum 5, 595 (2021).
  * Liyanage  _et al._ [2023] N. Liyanage, Y. Wu, A. Deters, and L. Zhong, Scalable quantum error correction for surface codes using FPGA (2023), arXiv:2301.08419 [quant-ph] .
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Skoric  _et al._ [2023] L. Skoric, D. E. Browne, K. M. Barnes, N. I. Gillespie, and E. T. Campbell, Parallel window decoding enables scalable fault tolerant quantum computation (2023), arXiv:2209.08552 [quant-ph] .
  * Higgott [2022] O. Higgott, PyMatching: A Python package for decoding quantum codes with minimum-weight perfect matching, ACM Transactions on Quantum Computing 3, 1 (2022).
  * Wu  _et al._ [2022] Y. Wu, N. Liyanage, and L. Zhong, An interpretation of union-find decoder on weighted graphs (2022), arXiv:2211.03288 [quant-ph] .
  * Tarjan [1975] R. E. Tarjan, Efficiency of a good but not linear set union algorithm, Journal of the ACM 22, 215 (1975).
  * Richards [2015] A. Richards,  _University of Oxford Advanced Research Computing_ (2015).

##  Appendix A Proof of 1

First consider the following.

######  Lemma 3 (Distributivity of $\sigma$).

The syndrome of the symmetric difference of errors is the symmetric difference
of their syndromes i.e. for $n\geq 2$:

$\sigma\bigg{(}{\operatorname*{{\triangle}}_{i=1}^{n}\mathbb{E}_{i}}\bigg{)}=\operatorname*{{\triangle}}_{i=1}^{n}\sigma(\mathbb{E}_{i}).$
(8)

###### Proof.

By induction. Consider base case $n=2$:

$\sigma(\mathbb{E}_{1}\mathbin{\vartriangle}\mathbb{E}_{2})\overset{\eqref{eq:syndrome_definition}}{=}V_{\text{bulk}}\cap\operatorname*{{\triangle}}_{e\in\mathbb{E}_{1}\mathbin{\vartriangle}\mathbb{E}_{2}}e.$
(9)

Since $e\mathbin{\vartriangle}e\equiv\varnothing$ we have
$\operatorname*{{\triangle}}_{e\in\mathbb{E}_{1}\mathbin{\vartriangle}\mathbb{E}_{2}}e\equiv(\operatorname*{{\triangle}}_{e\in\mathbb{E}_{1}}e)\mathbin{\vartriangle}(\operatorname*{{\triangle}}_{e\in\mathbb{E}_{2}}e)$
so

$\sigma(\mathbb{E}_{1}\mathbin{\vartriangle}\mathbb{E}_{2})\overset{\eqref{eq:syndrome_definition}}{=}\sigma(\mathbb{E}_{1})\mathbin{\vartriangle}\sigma(\mathbb{E}_{2})$
(10)

hence true for $n=2$. Now for the inductive step; assume true for $n=k$ and
consider for $n=k+1$:

$\displaystyle\operatorname*{{\triangle}}_{i=1}^{k+1}\mathbb{E}_{i}$
$\displaystyle=\bigg{(}{\operatorname*{{\triangle}}_{i=1}^{k}\mathbb{E}_{i}}\bigg{)}\mathbin{\vartriangle}\mathbb{E}_{k+1}$
$\displaystyle\bigg{(}{\operatorname*{{\triangle}}_{i=1}^{k+1}\mathbb{E}_{i}}\bigg{)}$
$\displaystyle\overset{\eqref{eq:syndrome_of_2_errors}}{=}\sigma\bigg{(}{\operatorname*{{\triangle}}_{i=1}^{k}\mathbb{E}_{i}}\bigg{)}\mathbin{\vartriangle}\sigma(\mathbb{E}_{k+1})$
$\displaystyle\overset{\eqref{eq:syndrome_of_n_errors}}{=}\operatorname*{{\triangle}}_{i=1}^{k+1}\sigma(\mathbb{E}_{i}).$
(11)

Hence, true for $n=k$ implies true for $n=k+1$. Together, the base case and
inductive step imply true for $n\geq 2$. ∎

Now we can prove 1.

###### Proof of 1.

Consider the syndrome produced by the leftover:

$\displaystyle(\mathbb{L})$
$\displaystyle\overset{\eqref{eq:syndrome_of_2_errors}}{=}\sigma(\mathbb{E})\mathbin{\vartriangle}\sigma(\mathbb{C})$
$\displaystyle=\mathbb{S\mathbin{\vartriangle}S}$ $\displaystyle=\varnothing.$
(12)

Hence, every bulk node is incident to an even number of edges in $\mathbb{L}$.
This can only happen if $\mathbb{L}$ comprises cycles or paths whose endpoints
are both boundary nodes. Since each boundary node is incident to exactly one
edge and $\mathbb{L}$ comprises distinct edges, these endpoints must be
distinct. ∎

##  Appendix B Proof of 2

We prove the following equivalent statement.

######  Lemma 4.

Cluster $C$ is inactive iff it has an even number of defects or touches a
boundary i.e.

$\big{(}{\textnormal{{\hbox{|\mathbb{S} \cap V_C|}}
even}}\big{)}\vee\big{(}{V_{\textnormal{boundary}}\cap
V_{C}\neq\varnothing}\big{)}.$ (13)

###### Proof.

  * $\Rightarrow$

We know $\exists\mathbb{C}\subseteq E_{C}:\sigma(\mathbb{C})=\mathbb{S}\cap
V_{C}$. If $C$ touches…

    * –

no boundary, $|\sigma(\mathbb{C})|$ even $\forall\mathbb{C}\subseteq E_{C}$ as
adding any edge in $E_{C}$ to $\mathbb{C}$ changes $|\sigma(\mathbb{C})|$ by
either 0 or 2.

    * –

a boundary, $|\sigma(\mathbb{C})|$ can be even or odd as adding an edge
containing a boundary node to $\mathbb{C}$ changes $|\sigma(\mathbb{C})|$ by
1.

  * $\Leftarrow$

Let $p:=|\mathbb{S}\cap V_{C}|$. If $p$ even, partition $\mathbb{S}\cap V_{C}$
into $p/2$ pairs. $C$ connected so there exists a path $P_{i}\subseteq E_{C}$
between the defects in each pair. Be
$\mathbb{C}:=\operatorname*{{\triangle}}_{i=1}^{p/2}P_{i}$ the symmetric
difference of these paths then

$\sigma(\mathbb{C})\overset{\eqref{eq:syndrome_of_n_errors}}{=}\operatorname*{{\triangle}}_{i=1}^{p/2}\sigma(P_{i})$
(14)

where $\sigma(P_{i})$ is precisely the $i^{\textnormal{th}}$ pair so

$\displaystyle(\mathbb{C})$ $\displaystyle=\bigcup_{i=1}^{p/2}\sigma(P_{i})$
$\displaystyle=\mathbb{S}\cap V_{C}.$ (15)

The existence of a correction satisfying Eq. 15 is the definition for $C$ to
be inactive. ∎


# Mu-suppression detection in motor imagery electroencephalographic signals
using the generalized extreme value distribution

Antonio Quintero-Rincón, Carlos D’Giano
Centro Integral de Epilepsia y Telemetría
Fundación Lucha contra las Enfermedades Neurológicas Infantiles (FLENI)
Buenos Aires, Argentina
<EMAIL_ADDRESS><EMAIL_ADDRESS>
&Hadj Batatia
IRIT - INPT
University of Toulouse
Toulouse, France
<EMAIL_ADDRESS>

###### Abstract

This paper deals with the detection of mu-suppression from
electroencephalographic (EEG) signals in brain-computer interface (BCI). For
this purpose, an efficient algorithm is proposed based on a statistical model
and a linear classifier. Precisely, the generalized extreme value distribution
(GEV) is proposed to represent the power spectrum density of the EEG signal in
the central motor cortex. The associated three parameters are estimated using
the maximum likelihood method. Based on these parameters, a simple and
efficient linear classifier was designed to classify three types of events:
imagery, movement, and resting. Preliminary results show that the proposed
statistical model can be used in order to detect precisely the mu-suppression
and distinguish different EEG events, with very good classification accuracy.

_K_ eywords Motor imagery $\cdot$ Mu-suppression $\cdot$ Generalized extreme
value $\cdot$ Electroencephalography $\cdot$ Brain-computer interface

##  1 Introduction

Electroencephalograms (EEG) are a non-invasive longstanding medical modality
that measures the brain’s activity by recording the electromagnetic field at
the scalp. Since its creation, EEG has played a fundamental role in
understanding several major neurological disorders, by analyzing their
manifestation into brain rhythms. For example, the study of deceases such as
depression, age-related cognitive deterioration, epilepsy, anxiety disorders
and subnormal brain development in children have benefited from this
technology. The typical brain rhythms are distinguished by their different
frequency ranges, called delta ($\delta$) within the range 0.5 to 4Hz, theta
($\theta$) within the range 4 to 7.5Hz, alpha ($\alpha$) within the range 8 to
13Hz, beta ($\beta$) within the range 14 to 30Hz, and gamma ($\gamma$) within
the range 30 to 64Hz. In this study, we focus on the brain rhythm called mu
($\mu$) within the range 7.5 to 11.5Hz. Mu-waves are considered to emerge
naturally and may convey information about what the functioning of brain
hierarchies [1]. According to [2], there exist three historical theoretical
hypotheses to explaining the mu-brain rhythm: i) the neuronal
hyperexcitability related to the rolandic cortex; ii) the superficial cortical
inhibition explaining its suppression with motor activity; and iii) the
somatosensory cortical idling, related to the afference-dependent phenomenon.
This study considers the second hypothesis, as the mu-rhythm relates strongly
to the sensorimotor cortex and associated areas, in particular, the changes in
the bilateral brain activities subject to physical and imaginary movements
[3]. Based on the same consideration, this rhythm has been studied in brain-
computer interface (BCI) [4, 5]. The underlying idea of BCI is to supply
communication and control of devices through the monitoring of brain activity,
by using EEG channels.

The generalized extreme value (GEV) distribution is a family that includes
continuous probability distributions obtained as the limit of maxima of a
sequence of independent and identically distributed random variables [6]. This
distribution has been used in [7] to detect interictal spikes in epileptic EEG
signals, using time-frequency properties. The underlying idea was to identify
strong outliers using GEV to model normalized EEG data. In [8], showed that
both EEG and MEG signals can be correctly modeled using GEV distribution. Luca
et al. [9] used an unsupervised method to detect hyper-motor epileptic
seizures, where GEV was applied for extracting maxima in EEG signals, using
multivariate kernel density. In [10], GEV was used to assess characteristics
of Alzheimer’s disease using EEG signals, where the variance of the power of
each frequency were used to derive an index of neuronal abnormality. For
applications in other biomedical signals see [11].

The purpose of this paper is to present a novel and rapid algorithm for
detecting mu-suppression in EEG signals by using the generalized extreme value
distribution. The underlying idea is to estimate the maximum and minimum
values of the signal in the central motor cortex using statistical modeling,
for motor imagery events, corresponding to mu-suppression. To the best of our
knowledge, this statistical model has not been investigated yet for detecting
the mu-suppression in EEG signals, despite the extensive study of this
phenomenon [12, 13, 14]. Several other methods have been proposed in the
literature to estimate mu-suppression in motor imagery, see [15, 16] for a
comprehensive state-of-the-art. Table 1 summarizes the most common methods,
such as independent components analysis, common spatial patterns, linear
discriminant analysis, and the support vector machines.

Method |  Features |  Electrodes |  Classf. |  Perf. |  Ref.
---|---|---|---|---|---
GEV from periodogram |  GEV parameters |  C5, C3, C1, Cz, C2, C4, C6 |  LDA |  100% |  our
AAR + BP |  Coefficients from AAR and BP |  C3, Cz |  LDA |  80% |  [17]
RFE |  Welch method from RFE |  C1, C2, C3, C4, FC3, FC4, CP3, CP4 |  SVM |  87% |  [18]
ICA vs. CSP |  log-transformed normalized variances |  C3, C4, Cz, CP1, CP2, CPz, Fz, Pz |  LDA |  84% |  [19]
CSP |  Eigenvalue spectrum |  multi-channels |  LDA |  84% |  [20]
EMD |  Periodograms fluctuations from intrinsic mode functions and the Hilbert evelope |  C3, C4 |  NR |  NR |  [21]
SCSP |  Variance of the spatially filtered signals |  C3, C4, Cz |  SVM |  94% |  [22]
ICA |  time-frequency average |  multi-channels |  LAT |  90% |  [16]
Table 1: State-of-the-art methods to the mu-suppression estimation in EEG
signals, compared in terms classification techniques, features, and reported
performance. Abbreviations are as follows: Adaptive autoregressive model
(AAR), Band power estimation (BP), Recursive feature elimination (RFE),
Independent components analysis (ICA), Common spatial patterns (CSP),
Empirical mode decomposition (EMD), Sparse common spatial pattern (SCSP),
Linear discriminant analysis (LDA), Support Vector Machine (SVM), According to
the latency (LAT), Not reported (NR). Figure 1(a) shows channel locations.
Note that our method (highlighted in red) had excellent accuracy results.

The remainder of the paper is organized as follows. In Section 2.2, the
generalized extreme value distribution is presented. Section 2 describes the
database used and the proposed methodology. Experimentation using the EEG
dataset of motor imagery from [23] is presented in Section 3 and the results
discussed. Finally, Section 4 draws conclusions and future works.

##  2 Materials and Methods

###  2.1 Motor Imagery database

A motor imagery EEG database of the left and right hands taken from [23] was
used. EEG data of 52 healthy subjects were recorded using a 64 channel montage
based on the international 10-10 system, with a 512 Hz sampling rate. In this
database, 3 events were taken into account. First, non-task, with the motor
cortex, data was recorded in order to understand each subject’s performance
related to their eye blinking, eyeball up/down movement, eyeball left/right
movement, head movement, jaw clenching, and resting state. Second, real right
and left-hand movement was measured, where subjects had to follow certain
established patterns. Finally, motor imagery data was computed when asking
subjects to imagine hand movement. The signals were preprocessed with a high
pass and low pass Butterworth filter of order 4th, having an 8-30 Hz bandpass.
Data was validated by using the percentage of bad trials, event-related
desynchronization/synchronization (ERD/ERS) analysis, and classification
analysis. The central motor cortex was selected for this study, see Fig. 1(a).

(a) Central motor cortex

(b) Channel Cz

Figure 1: Figure (a) shows the central motor cortex channels
(C5,C3,C1,Cz,C2,C4,C6) following the international 10-10 system. Figure (b)
shows an example of the central channel for each type of event. Imagery event
(grey) having the highest amplitude, movement event (yellow) with a medium
amplitude, and resting event (blue) corresponds to the lower amplitude. Note
that, imagery signals and movement signals may have very close amplitudes.

###  2.2 Statistical modeling

We propose to model the periodogram power spectral density of the EEG data
using the generalized extreme value distribution (GEV). The choice of this
distribution is motivated by its flexibility. It is a family that combines the
Gumbel (type I), Fréchet (type II) and Weibull (type III) maximum extreme
value distributions through three parameters, namely location
($\mu\in\mathcal{R}$), scale ($\sigma>0$), and shape ($\xi\in\mathcal{R}$ ).
Its probability density function (PDF) is given by

$\displaystyle f(x,\mu,\sigma,\xi)$
$\displaystyle=\frac{1}{\sigma}\;t(x)^{\xi+1}\exp^{-t(x)}$ (1)

where

$\displaystyle t(x)$
$\displaystyle=\begin{cases}(1+\xi(\frac{x-\mu}{\sigma}))^{-1/\xi}&\text{,
if}\ \xi\neq 0\\\ \exp^{\frac{x-\mu}{\sigma}}&\text{, if}\ \xi=0\end{cases}$
(2)

$\xi$ determines the domain of the PDF as follows:

$\displaystyle x\in\left[\frac{\mu-\sigma}{\xi},+\infty\right)$
$\displaystyle\text{, when }\xi>0\text{, type II}$ (3) $\displaystyle
x\in\left(-\infty,+\infty\right)$ $\displaystyle\text{, when }\xi=0\text{,
type I}$ $\displaystyle x\in\left(-\infty,\frac{\mu-\sigma}{\xi}\right]$
$\displaystyle\text{, when }\xi<0\text{, type III}$

In order to fit the GEV distribution to the data, we calculated the parameters
$\mu$ and $\sigma$ using the maximum likelihood estimators [24]

$\displaystyle\hat{\mu}$
$\displaystyle=-\hat{\sigma}\log\left(\frac{1}{n}\sum_{i=1}^{n}\exp^{\frac{X_{i}}{\hat{\sigma}}}\right)$
(4) $\displaystyle\hat{\sigma}$
$\displaystyle=\bar{X}-\frac{\sum_{i=1}^{n}X_{i}\exp^{-\frac{X_{i}}{\hat{\sigma}}}}{\sum_{i=1}^{n}\exp^{-\frac{X_{i}}{\hat{\sigma}}}}$
(5)

The cumulative density function (CDF) is given by

$\displaystyle f(x,\mu,\sigma,\xi)$ $\displaystyle=\exp^{-t(x)}$ (6)

for $x\in$ type II and type III.

We refer the reader to [25] for a comprehensive treatment of the mathematical
properties of the extreme value distribution.

###  2.3 Proposed methodology

Let $\boldsymbol{X}\in\mathbb{R}^{N\times M}$ denote the matrix of $M$ EEG
signals $\boldsymbol{x}_{m}\in\mathbb{R}^{N\times 1}$ measured simultaneously
on different channels related to the central motor cortex and at $N$ discrete-
time instants. The proposed methodology is composed of three stages. In the
first-stage, the periodogram power spectral density (PSD) is estimated from
the EEG signal on the mu-brain rhythm, in steps of half the sampling frequency
without overlapping using a symmetric Hamming window

$\displaystyle P(f)=\frac{\Delta
t}{N}\left|\sum_{n=0}^{N-1}x_{m,n}h(n)\exp^{-2j\pi f\Delta
t\;n}\right|^{2},\\-\frac{1}{2}\Delta t<f\leq\frac{1}{2}\Delta t$ (7)

where $\Delta t$ is the sampling interval, $1\leq m\leq 7$ is related to the
seven channels in the central motor cortex (C5, C3, C1, Cz, C2, C4, C6) and
the Hamming window defined as
$h(n)=0.54-0.46\;cos\left(2\pi\frac{n}{N}\right)$, $0\leq n\leq N$.

In the second-stage, the generalized extreme value distribution is fitted to
the PSD of each event (imagery, movement and resting), yielding the
corresponding parameters ($\mu$, $\sigma$ and $\xi$):

$\displaystyle I_{i}$ $\displaystyle=[\mu_{i},\sigma_{i},\xi_{i}]$ (8)
$\displaystyle M_{i}$ $\displaystyle=[\mu_{i},\sigma_{i},\xi_{i}]$ (9)
$\displaystyle R_{i}$ $\displaystyle=[\mu_{i},\sigma_{i},\xi_{i}]$ (10)

Finally in the third stage, the vectors (8),(9) and (10) are classified using
a linear classifier derived from a linear discriminant analysis (LDA)[26, 27].
Three classes $\\{0,1,-1\\}$ are possible, respectively for imagery, movement
and resting.

##  3 Results and discussion

Figures 2, 3, and 4 show the Goodness-of-fit measures for the probability
density function (PDF) and the cumulative distribution density (CDF), and the
CDF error for the central motor cortex data. These measures are clearly
significant. This is corroborated by the CDFs errors ($-0.05\leq\text{CDF
error}\leq 0.05$), and the $p$-values $\leq 0.05$ using the Kolmogorov-Smirnov
(KS) score, where the imagery events exhibit a $p=0.0021$, the movement events
a $p=0.0094$, and the resting events a $p=1.9847e-04$.

Tables 2 to 4 show the mean values, with their $95\%$ confidence intervals,
for the three GEV parameters corresponding to the different events. It appears
clearly that the proposed statistical model allows for a threshold approach to
detect the mu-suppression by distinguishing imagery, from movement and resting
events. For illustration, in EEG data shown in Figure 1(b), the imagery event
has the highest amplitude, the movement event has the medium amplitude, and
the resting event has a lower amplitude. While in Figure 5, showing averages
from Table 2, the scale ($\sigma$) and the location ($\mu$) parameters have
the highest amplitude for resting state, the medium amplitude for imagery
event, and the lower amplitude for movement event. Whereas, the shape ($\xi$)
parameter shows the movement event with the highest amplitude, the resting
event with the medium amplitude and the imagery state with the lower
amplitude. The combination of these different thresholds makes it possible to
use the three parameters, estimated by fitting the GEV distribution to the
periodogram, for detecting the mu-suppression in motor imagery EEG signals.

Figure 2: Mu-band in imagery event

Figure 3: Mu-band in movement event

(a) Resting (b) Resting

Figure 4: Mu-band in resting event | shape ($\xi$)
---|---
Ima | 0.09 [-0.21,0.38]
Mov | 0.26 [-0.06,0.58]
Res | 0.15 [-0.26,0.56]
Table 2: Mean values of the shape parameter of the GEV adjusted to the
periodogram of mu-band, with the 95% confidence intervals.

| scale ($\sigma$)
---|---
Ima | 143275.14 | [110712.36, | 185415.39]
Mov | 31104.66 | [23545.71, | 41090.27]
Res | 1273712.63 | [945447.56, | 1715953.13]

Table 3: Mean values of the scale parameter of the GEV adjusted to the
periodogram of mu-band, with the 95% confidence intervals

.

| location ($\mu$)
---|---
Ima | 196975.08 | [149574.66, | 244375.50]
Mov | 33551.02 | [23148.30, | 43953.76]
Res | 1241060.75 | [778886.63, | 1703234.88]

Table 4: Mean values of the location parameter of the GEV adjusted to the
periodogram of mu-band, with the 95% confidence intervals

.

Figure 5: Plot of the mean values from Tables 2 to 4. Note that, these
parameters allow the identification of imagery events (grey), the movement
events (yellow), and the resting events (blue) using a threshold approach.

To assess the quality of the proposed method for mu-band suppression, we
adopted a supervised classification approach with 30 thousand randomized
signals for training and testing from the motor imagery database with the 3
GEV features for each vector as follows:

  * •

$\theta_{1}=[\mu_{1},\sigma_{1},\xi_{1}]$ 10 thousand randomized imagery
events.

  * •

$\theta_{0}=[\mu_{0},\sigma_{0},\xi_{0}]$ 10 thousand randomized movement
events.

  * •

$\theta_{-1}=[\mu_{-1},\sigma_{-1},\xi_{-1}]$ 10 thousand randomized resting
events.

A linear classifier derived from the linear discriminant analysis (LDA) was
designed to distinguish between the different events
$[\theta_{1},\theta_{0},\theta_{-1}]$. Note that, the LDA computational
complexity is lower with respect to support vector machines (SVM), which is
one of the most used classifiers. The LDA computational complexity is
$O(mn+mt+nt)$, where $m$ is the number of samples, $n$ is the number of
features, and $t=min(m,n)$. While SVM complexity is $O(n^{2})$ or $O(n^{3})$
depending on the parametrization, with $n$ being the size of the dataset,
which is generally large when EEG signals are analyzing [27].

Due to the big data volume generated, a 20-fold cross-validation technique was
used for training and testing, with 10 times repetition to ensure no bias in
the partitioning and to avoid imbalanced data [28]. Figure 6, shows the
scatter distribution of the different randomized events, namely imagery,
movement, and resting. It is interesting to note the separability of the
classes. This justifies evidently the use of linear classification to
distinguish events. The classifier achieves a 100% sensitivity (True positive
rate) and specificity (True false rate) for mu-suppression detection.

Figure 6: Scatter plot example with randomized data. It is clearly possible to
use a linear classifier as imagery events (grey), the movement events
(yellow), and the resting events (blue) are well separated in the GEV
parameter space.

##  4 Conclusions

This work presented a study of detecting mu-suppression in EEG signals, based
on the generalized extreme value (GEV) distribution. The three parameters of
this distribution were estimated from its periodogram power spectral density
for three events: imagery, movement, and resting. A linear classifier derived
from the linear discriminant analysis (LDA) was designed to distinguish
between event types using these parameters. The performance of the proposed
method was evaluated on a real dataset from 52 subjects achieving 100%
accuracy. In addition to its performance, an advantage of this method is its
low computational cost compared to existing methods. These good results have
the potential to shed new light on mu-suppression detection in motor imagery
EEG signals in the central motor cortex.

The noise and artifacts were not taken into consideration in this work, which
constitutes its main limitation. Future work will focus on the study of the
noise and artifacts in order to make the method applicable in real-time. A
large scale testing campaign will also be undertaken with other databases and
other possible channel locations with the idea of using the least amount of
channels possible.

##  5 Acknowledgments

The authors would like to thank Joaquin Ems, Lourdes Hirschson, and Catalina
Carenzo for useful comments on an earlier version of the manuscript. We are
grateful to the professor emeritus Jaime A. Pineda from the University of
California, San Diego, La Jolla, CA, USA, for all the useful discussions
during the process of this study.

## References

  * [1] K. Friston, “Waves of prediction,” _PLOS Biology_ , vol. 10, no. 17, p. e3000426, 2019.
  * [2] D. L. Schomer and F. L. da Silva, _Niedermeyer’s Electroencephalography: Basic Principles, Clinical Applications, and Related Fields_. Lippincott Williams and Wilkins, 2010.
  * [3] S. Sanei, _Adaptive Processing of Brain Signals_. Wiley-Blackwell, 2013.
  * [4] S. Arroyo, R. P. Lesser, B. Gordon, S. Uematsu, D. Jackson, and R. Webber, “Functional significance of the mu rhythm of human cortex: an electrophysiologic study with subdural electrodes,” _Electroencephalography and Clinical Neurophysiology_ , vol. 3, no. 87, pp. 76–87, 1993.
  * [5] G. Pfurtscheller and C. Neuper, “Event-related synchronization of mu rhythm in the EEG over the cortical hand area in man,” _Neuroscience Letters_ , vol. 1, no. 174, pp. 93–96, 1994.
  * [6] E. Castillo, A. S. Hadi, N. Balakrishnan, and J. M. Sarabia, _Extreme Value and Related Models with Applications in Engineering and Science_. Wiley-Interscience, 2004.
  * [7] N. Gaspard, R. Alkawadri, P. Farooque, I. I. Goncharova, and H. P. Zaveri, “Automatic detection of prominent interictal spikes in intracranial EEG: Validation of an algorithm and relationsip to the seizure onset zone,” _Clinical Neurophysiology_ , vol. 6, no. 125, pp. 1095–1103, 2014.
  * [8] N. Nazmia, S. A. Mazlana, H. Zamzuria, M. Azizi, and A. Rahman, “Fitting distribution for electromyography and electroencephalography signals based on goodness-of-fit tests,” _IEEE International Symposium on Robotics and Intelligent Sensors_ , no. 76, pp. 468–473, 2015.
  * [9] S. Luca, P. Karsmakers, K. Cuppens, T. Croonenborghs, A. V. de Vel, B. Ceulemans, L. Lagaee, and B. V. Sabine Van Huffel a, b, “Detecting rare events using extreme value statistics applied to epileptic convulsions in children,” _Artificial Intelligence in Medicine_ , no. 60, pp. 89–96, 2014\.
  * [10] T. Ueda, T. Musha, and T. Yagi, “Research of the characteristics of alzheimer’s disease using EEG,” _Annual International Conference of the IEEE EMBS_ , no. 2009, pp. 4998–5001, 2009.
  * [11] S. J. Roberts, “Extreme value statistics for novelty detection in biomedical data processing,” _IEE Proceedings - Science, Measurement and Technology_ , no. 147, pp. 363–367, 2000.
  * [12] J. Pineda, “The functional significance of mu rhythms: Translating "seeing" and "hearing" into "doing",” _Brain Research Reviews_ , vol. 68, no. 57, pp. 57–68, 2005.
  * [13] H. M. Hobson and D. V. M. Bishop, “The interpretation of mu suppression as an index of mirror neuron activity: past, present and future,” _Royal Society Open Science_ , vol. 3, no. 4, p. 160662, 2017.
  * [14] Y. Nishimura1, Y. Ikeda1, A. Suematsu1, and S. Higuchi, “Effect of visual orientation on mu suppression in children: a comparative EEG study with adults,” _Journal of Physiological Anthropology_ , vol. 16, no. 37, pp. 1–9, 2018.
  * [15] M. Hamedi, S.-H. Salleh, and A. M. Noor, “Electroencephalographic motor imagery brain connectivity analysis for BCI: A review,” _Neural Computation_ , no. 6, pp. 999–1041, 2016.
  * [16] K. Xu, Y.-Y. Huang, and J.-R. Duann, “The sensitivity of single-trial mu-suppression detection for motor imagery performance as compared to motor execution and motor observation performance,” _Frontiers in Human Neuroscience_ , vol. 13, no. 302, pp. 1–12, 2019.
  * [17] C. Guger, G. Edlinger, W. Harkam, I. Niedermayer, and G. Pfurtscheller, “How many people are able to operate an EEG-based brain-computer interface (BCI)?” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 2, no. 11, pp. 145–147, 2003.
  * [18] M. Schroder, T. N. Lal, T. Hinterberger, M. Bogdan, N. Hill, N. Birbaumer, W. Rosenstiel, and B. Scholkopf, “Robust EEG channel selection across subjects for brain-computer interfaces,” _EURASIP Journal on Applied Signal Processing_ , no. 19, pp. 3103–3112, 2005.
  * [19] M. Naeem, C. Brunner, R. Leeb, B. Graimann, and G. Pfurtscheller, “Seperability of four-class motor imagery data using independent components analysis,” _Journal of Neural Engineering_ , no. 3, pp. 208–216, 2006.
  * [20] B. Blankertz, F. Losch, M. Krauledat, G. Dornhege, G. Curio, and K.-R. Muller, “The Berlin brain-computer interface: Accurate performance from first-session in BCI-naive subjects,” _IEEE Transactions on Biomedical Engineering_ , vol. 10, no. 55, pp. 2452–2462, 2008.
  * [21] B. Wan, Z. Zhou, L. Xu, D. Ming, H. Qi, and L. Cheng, “Mu rhythm desynchronization detection based on empirical mode decomposition,” _31st Annual International Conference of the IEEE EMBS_ , pp. 2232–2235, 2009\.
  * [22] M. Arvaneh, C. Guan, K. K. Ang, and C. Quek, “Optimizing the channel selection and classification accuracy in EEG-based BCI,” _IEEE Transactions on Biomedical Engineering_ , vol. 6, no. 58, pp. 1865–1873, 2011.
  * [23] H. Cho, M. Ahn, S. Ahn, M. Kwon, and S. C. Jun, “EEG datasets for motor imagery brain computer interface,” _Gigascience_ , vol. 7, no. 6, pp. 1–8, 2017.
  * [24] R. B. D’Agostino and M. A. Stephens, _Goodness-of-Fit Techniques_. Dekker, 1986.
  * [25] S. Kotz and S. Nadarajah, _Extreme Value Distributions: Theory and Applications_. ICP, 2001.
  * [26] A. Quintero-Rincón, M. Pereyra, C. D’Giano, M. Risk, and H. Batatia, “Fast statistical model-based classification of epileptic EEG signals,” _Biocybernetics and Biomedical Engineering_ , vol. 38, no. 4, pp. 877–889, 2018.
  * [27] A. Quintero-Rincón, C. D’Giano, and M. Risk, “Epileptic seizure prediction using Pearson’s product-moment correlation coefficient of a linear classifier from generalized Gaussian modeling,” _Neurología Argentina_ , vol. 10, no. 4, pp. 201–217, 2018.
  * [28] A. Quintero-Rincón, M. Flugelman, J. Prendes, and C. D’Giano, “Study on epileptic seizure detection in EEG signals using largest Lyapunov exponents and logistic regression,” _Revista Argentina de Bioingeniería, Bioengineering Argentinian Society_ , vol. 23, no. 2, pp. 17–24, 2019.