# A Logarithmic Bayesian Approach to Quantum Error Detection

Ian Convy<EMAIL_ADDRESS>Department of Chemistry, University of
California, Berkeley, CA 94720, USA  Berkeley Quantum Information and
Computation Center, University of California, Berkeley, CA 94720, USA  K.
Birgitta Whaley  Department of Chemistry, University of California, Berkeley,
CA 94720, USA  Berkeley Quantum Information and Computation Center, University
of California, Berkeley, CA 94720, USA

###### Abstract

We consider the problem of continuous quantum error correction from a Bayesian
perspective, proposing a pair of digital filters using logarithmic
probabilities that are able to achieve near-optimal performance on a three-
qubit bit-flip code while still being reasonable to implement on low-latency
hardware. These practical filters are approximations of an optimal filter that
we derive explicitly for finite time steps, in contrast with previous work
that has relied on stochastic differential equations such as the Wonham
filter. By utilizing logarithmic probabilities, we are able to eliminate the
need for explicit normalization and reduce the Gaussian noise distribution to
a simple quadratic expression. The state transitions induced by the bit-flip
errors are modeled using a Markov chain, which for log-probabilties must be
evaluated using a LogSumExp function. We develop the two versions of our
filter by constraining this LogSumExp to have either one or two inputs, which
favors either simplicity or accuracy, respectively. Using simulated data, we
demonstrate that the one-term and two-term filters are able to significantly
outperform a double threshold scheme and a linearized version of the Wonham
filter in tests of error detection under a wide variety of error rates and
time steps.

##  1 Introduction

One of the major obstacles to large-scale quantum computation is the high
frequency of qubit errors. Small interactions with the environment or
imperfections in gate implementation can perturb the underlying quantum state
throughout a computation and ultimately render the output useless. These
issues only magnify as the size and complexity of the quantum computer
increases, jeopardizing any attempt to demonstrate quantum advantage in key
tasks such as the factorization of large integers [1]. The obvious challenge
posed by quantum errors has given rise to the field of quantum error
correction [2] and spurred the development of numerous techniques to identify
and correct errors when they occur. These schemes, referred to as error
correction codes, typically operate by encoding the quantum state into a
larger Hilbert space and then monitoring the location of the state within this
space using a specific set of observables.

The present work is focused on continuous or “always-on” error correction,
where the state of the system is continuously monitored through a noisy signal
channel in order to diagnose errors and to make corrections if necessary.
Since the signal is noisy, we cannot know with certainty whether an error has
occurred, and must instead make subjective judgements about the condition of
the system based on past measurements and any other information we might have.
In order to quantify the uncertainty in this task it is natural to adopt a
Bayesian framework, which allows for new information to be easily combined
with prior knowledge of the system using Bayes’ theorem.

Assuming that the errors experienced by a quantum computer are uncorrelated,
we can model the system as a Markov chain whose state is imperfectly observed
by a measurement apparatus. Optimal treatments of this problem originate with
work done by Wonham and his derivation of what is today known today as the
Wonham filter, which describes the evolution of the Markov state probabilities
conditioned on a particular set of noisy observations [3]. Despite its
optimality, this filter has not found widespread use for error correction, in
part because it assumes that the value of the signal is characterized as a
continuous function of time, while in practice the signal is typically only
observed as a discrete set of measurement samples. As a non-linear stochastic
differential equation, the Wonham filter can also be challenging to evaluate
numerically even if the signal function is known. Due to these limitations,
most continuous error correction schemes avoid using Bayesian theory and
instead rely on various thresholding procedures which are easier to implement
but known to provide suboptimal performance [4][5][6].

The purpose of this work is to introduce a discrete-time filter for continuous
quantum error correction that significantly outperforms common thresholding
schemes, while still being practical to implement on real hardware. We
construct our filter using logarithmic probabilities, which are numerically
stable, easy to normalize, and allow for straightforward evaluation of
Gaussian noise distributions using only arithmetic operations. The filter has
two different forms based on how it updates the posterior, with one method
emphasizing accuracy while the other emphasizes computational simplicity. Both
versions of the filter are able to achieve near optimal performance under a
wide variety of error rates and time steps.

This paper has the following structure. Sec. 2 reviews the three-qubit code
that we will use as a simple error correction primitive, and describes the
principles of continuous measurement. Sec. 3 derives an optimal Bayesian
filter for identifying bit-flip errors, which cannot be practically
implemented in its exact form. Sec. 4 discusses prior work in the field of
continuous quantum error correction and outlines several obstacles to
implementing a Bayesian filter. Sec. 5 proposes an approximate version of the
optimal filter which overcomes those obstacles, and provides two different
implementations which are both experimentally feasible using today’s
technology. Sec. 6 tests our filter on simulated data and compares the
performances of the two implementations against other schemes found in the
literature. Sec. 7 summarizes our findings and identifies areas for future
work. An Appendix is also provided, which contains technical details that
supplement the body of the paper.

##  2 Quantum Error Correction

###  2.1 Three-qubit coding scheme

To insulate a quantum system from errors, we must add some level of redundancy
to its Hilbert space. This is achieved by assembling a set of physical qubits
which is larger than necessary to perform the desired task, with the idea that
these extra degrees of freedom will be used to implement an error correction
scheme. The computation is then understood in terms of its effect on the
coding subspace of the system, which describes the portion of the expanded
Hilbert space where the intended quantum evolution will take place.

For simplicity we focus on a three-qubit repetition code, which is designed to
protect the state of a singe qubit from the effects of bit-flip errors. This
encoded qubit is the true object of interest, so the computational task is
conceptualized in terms of single-qubit operations which then need to be
mapped to the full three-qubit system. We can perform this mapping by
considering the so-called logical qubit, defined as

$\begin{split}\ket{0_{L}}&\equiv\frac{1}{\sqrt{2}}(\ket{000}+\ket{111})\\\
\ket{1_{L}}&\equiv\frac{1}{\sqrt{2}}(\ket{000}-\ket{111}),\end{split}$ (1)

which occupies the two-dimensional subspace spanned by $\ket{000}$ and
$\ket{111}$ that we designate as the coding subspace.

The logical states in Eq. (1) are used to encode the desired behavior of the
system in the absence of errors, while the other parts of the Hilbert space
serve to catch any bit-flips that occur. For example, a bit-flip error on the
first qubit will transform $\ket{000}$ into $\ket{100}$ and $\ket{111}$ into
$\ket{011}$, so the system will be shifted into a different two-dimensional
subspace. As long as we can detect this change in the subspace, it will be
possible for us to track and then correct the errors that have occurred on the
system.

To perform this detection, we require a set of observables that uniquely
identify the coding subspace and each of the three error subspaces generated
from bit-flips on the first, second, and third qubits. This is achieved using
the parity operators $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$, whose eigenstates and
eigenvalues are

$(Z_{1}Z_{2}\ket{\psi},Z_{2}Z_{3}\ket{\psi})=\begin{cases}(\ket{\psi},\ket{\psi})&\ket{\psi}=a\ket{000}+b\ket{111}\\\
(-\ket{\psi},\ket{\psi})&\ket{\psi}=a\ket{100}+b\ket{011}\\\
(\ket{\psi},-\ket{\psi})&\ket{\psi}=a\ket{001}+b\ket{110}\\\
(-\ket{\psi},-\ket{\psi})&\ket{\psi}=a\ket{010}+b\ket{101},\end{cases}$ (2)

for any complex numbers $a$ and $b$. Eq. (2) shows that the combined action of
the parity operators splits the three-qubit Hilbert space into four distinct
2-D subspaces, each characterized by a different pair of eigenvalues in
$\\{-1,1\\}$ which are referred to as syndromes. The $(1,1)$ syndrome
identifies the coding subspace, while the other three combinations correspond
to states which result from a bit-flip on one of the qubits. We refer to the
latter as the error subspaces, and note that the syndrome value identifies on
which qubit the error has occurred.

The observables $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$ offer us a simple method for
bit-flip error detection, assuming that they can in fact be measured. After
initializing the system into the coding subspace, we simply perform parity
measurements of our system as needed and monitor the response of the system.
If our results indicate that the system has moved out of the coding subspace,
then we can record this deviation and apply a bit-flip correction based on
which error subspace the system has been moved to. For example, if one error
occurs on the first qubit then we have

$a\ket{000}+b\ket{111}\xrightarrow[\text{qubit
1}]{\text{error}}a\ket{100}+b\ket{011}\xrightarrow[\text{qubit
1}]{\text{correction}}a\ket{000}+b\ket{111},$ (3)

and the system has been returned to its proper state. Each time a bit-flip is
diagnosed to have occurred on a given qubit, we can simply apply another bit-
flip to that same qubit and the error will be undone.

It is important to note that if two bit-flip errors occur simultaneously (or
rapidly enough to appear indistinguishable) then the above detection scheme
will mistakenly assume that only one error has occurred and perform the wrong
correction. For example, if errors occur on the second and third qubits then

$a\ket{000}+b\ket{111}\xrightarrow[\text{qubits
2,3}]{\text{errors}}a\ket{011}+b\ket{100}\xrightarrow[\text{qubit
1}]{\text{correction}}a\ket{111}+b\ket{000},$ (4)

and the system will experience a logical error since the $a$ and $b$
coefficients have been exchanged. This type of misdiagnosis occurs because the
errors in Eq. (3) and Eq. (4) both move the system into the same parity
subspace and thus return the same syndrome values, so we must respond with the
same correction in both cases. Given that single errors are more common then
double errors under realistic conditions, the most reasonable response is to
always assume that one bit-flip has occurred whenever the system enters an
error subspace. This compromise is an inescapable consequence of using a short
code, though we can reduce its impact by detecting errors more quickly and
thus better resolving the difference between one bit-flip and two sequential
bit-flips.

Throughout this paper we will focus our attention on quantum memory, where it
is assumed that the system is initialized into $\ket{000}$ and only evolves
due to random bit-flip errors. However, Eq. (2) shows that the parity
operators act uniformly on states within a given coding or error subspace,
meaning that the measurement signal will be identical regardless of whether
the system is in a basis state or an arbitrary superposition. Even if a
Hamiltonian is applied to the system, as would occur during a process like
quantum annealing [7], the measurement signal will not be affected so long as
the Hamiltonian does not mix together the four parity subspaces. As a result,
any filter designed to handle quantum memory can be seamlessly applied to more
complicated processes by just changing the metric used to evaluate its
performance.

###  2.2 Error correction with continuous measurements

The three-qubit code described in Sec. 2.1 is agnostic to the actual method
used to probe $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$, since the separation of the
Hilbert space into coding and error subspaces is a property of the observables
themselves rather than of any particular measurement scheme. Most work on
quantum errror detection has focused on discrete error correction, which
utilizes projective measurements to resolve the exact syndrome of the state at
periodic intervals [2]. However, recent advances in superconducting qubit
architecture has led to increased interest in continous error correction,
where the parity operators are monitored at all times and a noisy readout of
the underlying syndrome values is generated [8][9]. Under these noisy
conditions a filter will be needed in order to extract relevant information
from the signal.

In an idealized setting, we can formulate a simple mathematical description of
the continuous measurement process using Gaussian POVMs. At time $t$, the
probability of observing signal readout $\alpha(t)$ from a parity operator
with syndrome $S(t)$ is

$\text{Prob}(\alpha(t))\propto\exp[-\frac{(S(t)-\alpha(t))^{2}}{2\sigma^{2}}],$
(5)

which implies that the signal values are distributed as
$\mathcal{N}(S(t),\sigma^{2})$. The syndrome value $S(t)\in\\{-1,+1\\}$ is
determined by the state of the system at time $t$, and will therefore change
as errors occur. We assume that the variance $\sigma^{2}$ of the noisy signal
is identical across all states and is thus independent of time.

The weak measurement described by Eq. (5) must occur over some finite time
$\Delta t$, whose duration will determine the variance of the signal since
$\sigma^{2}=\frac{k}{\Delta t}$ for some positive $k$ (the strength of the
measurement is thus proportional to $\frac{1}{k}$). A longer measurement
$\tilde{M}$ of time $T$ can be constructed by averaging over a set of these
weaker measurements

$\tilde{M}(t,T,\Delta t)=(\frac{T}{\Delta t})^{-1}\sum^{\frac{T}{\Delta
t}}_{n=1}\alpha(t+n\Delta t)=\frac{1}{T}\sum^{\frac{T}{\Delta
t}}_{n=1}[S(t+n\Delta t)+\epsilon]\Delta t,$ (6)

where the signal at time increment $t+n\Delta t$ is the sum of the syndrome
$S(t+n\Delta t)$ and a Gaussian noise term
$\epsilon\sim\mathcal{N}(0,\frac{k}{\Delta t})$. The value of $t$ represents
the time at which the averaging process starts, while the integer $n$ indexes
the set of sequential weak measurements being averaged over, each of duration
$\Delta t$. Since the average is taken over a time interval of length $T$, we
can perform $\frac{T}{\Delta t}$ weak measurements within the interval.

Given that the Gaussian noise $\epsilon$ in Eq. (6) is additive, its average
can be calculated separately from the syndrome average as

$\xi(T)\equiv\frac{\Delta t}{T}\sum^{\frac{T}{\Delta
t}}_{n=1}\epsilon=\frac{\Delta t}{T}\tilde{\epsilon},$ (7)

where $\tilde{\epsilon}\sim\mathcal{N}(0,\frac{kT}{(\Delta t)^{2}})$ due to
the additive property of the variance and thus
$\xi(T)\sim\mathcal{N}(0,\frac{k}{T})$ after scaling by $\frac{\Delta t}{T}$.
The value of $\tilde{M}$ is therefore equal to the averaged syndrome plus a
new noise term $\xi(T)$, whose variance is inversely proportional to $T$:

$\tilde{M}(t,T,\Delta t)=\frac{1}{T}\sum^{\frac{T}{\Delta
t}}_{n=1}[S(t+n\Delta t)\Delta t]+\xi(T).$ (8)

From Eq. (8), we define the continuous measurement $M$ as the limit of
$\tilde{M}$ for $\Delta t\rightarrow 0$, which turns the discrete sum into an
integral,

$M(t,T)=\frac{1}{T}\int^{T}_{0}S(t+\tau)d\tau+\xi(T),$ (9)

that is distributed as a Gaussian with variance $\frac{k}{T}$ and mean value
given by the integrated average of the syndrome over the interval.

Figure 1: Plots A and B show examples of 30 $\mu$s measurement sequences for
$T=1\ \mu\text{s}$ and $T=0.1\ \mu\text{s}$ respectively, with $k=0.5$ for
both plots. Each point represents the unitless value of $M(t,T)$ taken at the
time shown on the x-axis for one of the two parity operators, with the blue
lines tracking the noisy measurement of $Z_{1}Z_{2}$ and the orange lines
tracking $Z_{2}Z_{3}$. The underlying state evolution in both plots is
identical, with the vertical black lines indicating the times at which bit-
flip errors occurred (10.7 $\mu$s and 20 $\mu$s). The state of the system
before and after each error is given by the three-bit number above plot A.
Plot C shows the true syndrome value from $Z_{1}Z_{2}$ without noise, zoomed
in on times near the first bit-flip. The vertical gray lines indicate the
integration intervals used to generate the averaged syndromes found in Eq. (9)
when $T=1\ \mu\text{s}$, with the numerical values of these averages given at
the bottom of the plot. The shaded regions indicate the “area under the curve”
contribution at a given time, with green indicating a positive contribution
and red a negative contribution. The bit-flip error occurs within the second
interval, causing its syndrome average to deviate from the expected values of
$\pm 1$ and instead lie somewhere in between.

In Figure 1 we provide examples of measurement sequences that could be
observed as a three qubit system evolves under the influence of random bit-
flip errors for 30 $\mu$s. For a system starting in $\ket{000}$, the
measurements $M(t,T)$ for each syndrome will be centered at +1 and oscillate
randomly due to the additive noise term. Comparing Figures 1A and 1B, it is
clear that increasing the integration time $T$ dramatically reduces the
fluctuations of the signal, making it much easier to see visually the effects
of the bit-flip errors on the measurements. Perhaps counterintuitively, this
reduction in the variance does not actually make the task of identifying a
syndrome value from its noisy measurements any easier, since the standard
error (SE) of mean value estimation scales with both the noise strength
$\frac{k}{T}$ and the number of samples to be averaged. Given that the number
of measurements $n$ in a time interval of length $\beta$ is inversely
proportional to $T$, i.e. $n=\frac{\beta}{T}$, we have

$\text{SE}^{2}=\frac{\sigma^{2}}{n}=\frac{k}{T}\frac{T}{\beta}=\frac{k}{\beta}$
(10)

and thus the uncertainty in our estimate of the syndrome will not depend on
$T$.

This is not to say, however, that the length of the integration period has no
effect on the state tracking problem. Indeed, Eq. (10) simply states that
there is a proportionate trade-off between the informativeness of a
measurement and the number of such measurements if the syndrome value is
fixed. When we introduce bit-flip errors into the system, and thus consider
situations where the syndrome changes with time, the analysis is significantly
more complicated. Figure 1C shows that when an error occurs within an
integration period, a portion of the syndrome average will consist of +1
contributions and the other portion will consist of -1 contributions. The
resulting syndrome mean can therefore take on any value in the interval
$[-1,+1]$, which suggests that the measurement will be inherently less
informative than a measurement that occurs when the state of the system is
constant. As the value of $T$ grows larger, more measurements will fall into
this category. The incorporation of these intermediate syndrome means into an
optimal state tracking model is explored in Sec. 3.

###  2.3 Bayesian treatment of measurement

Determining the state of a quantum system from noisy signals is not a new
problem, as the field of quantum mechanics has long grappled with how to
properly describe the mechanism and effects of measurement [10]. One of the
key challenge lies in knowing how to properly update the density matrix
$\rho_{s}$ of a system to reflect the fact that we have performed a
measurement on it. A common procedure involves coupling the system to a
detector and then tracing out the detector’s degrees of freedom,

$\rho_{s}=\text{Tr}_{d}(\rho_{sd})$ (11)

where $\rho_{sd}$ is the joint density matrix of the system and detector. By
its very nature this approach cannot offer any information about the state of
the system after a specific measurement, since we are averaging over all of
the possible detector configurations without reference to the outcome that
really occurred. Instead, what Eq. (11) describes is the behavior of an
ensemble of systems after they are all measured, since the frequency of each
measurement outcome can be specified without knowing the result for any given
system. This approach is therefore insufficient for an error correction
scheme, as the very purpose of error correction is to monitor a specific
system and respond to the errors that actually occur rather than to the
average over all possible errors.

The issue of specific outcomes versus ensemble behavior is reminiscent of the
disagreement between Bayesian and frequentist statistics [11], and it is
therefore unsurprising that a Bayesian formalism for quantum measurement was
developed. This formalism was pioneered by Korotkov, who proposed an update
rule for the diagonal elements of the density matrix after a measurement $M$
using Bayes’ theorem, which has the form

$\tilde{\rho}_{ii}=\frac{P(M|i)\rho_{ii}}{\sum_{j}P(M|j)\rho_{jj}},$ (12)

where $\rho_{ii}$ is the original diagonal element acting as the prior and
$\tilde{\rho}_{ii}$ is the updated (posterior) element [12]. The term $P(M|i)$
describes the probability of observing measurement result $M$ given that the
system is in state $i$.

##  3 An Optimal Bayesian Filter

For our purposes, the significance of Korotkov’s measurement formalism from
Eq. (12) lies in its treatment of the density matrix elements as classical
Bayesian probability distributions, which allows us to leverage the existing
toolkit of Bayesian statistics [13][14] to develop an error correction
algorithm. Fundamentally, the goal of this algorithm is to determine, using
all available information, the most probable state of the system at each time
step in order to detect errors. As in Eq. (12), we seek to to derive a
posterior probability distribution for the state of our system after each
measurement, using knowledge of the underlying bit-flip dynamics and of the
Gaussian noise corrupting the measurement signal.

###  3.1 Recursive form of the posterior probability

For notational simplicity, we introduce the vector-valued quantity

$\vec{M}_{i}\equiv\begin{bmatrix}M^{(1)}(iT,T)\\\ M^{(2)}(iT,T)\end{bmatrix},$
(13)

where $M^{(1)}$ and $M^{(2)}$ are measurements of the parity operators
$Z_{1}Z_{2}$ and $Z_{2}Z_{3}$ respectively and $i$ is a non-negative integer.
The time argument $iT$ appears because the sets of sequential measurements
must necessarily be spaced out in increments of $T$. To denote the state of
the system after the $i$th measurement, i.e. at time $t=(i+1)T$, we employ the
compound index $0\leq\ell_{i}\leq 7$. Using this notation, we write the
posterior probability of the $\ell_{i}$th state as

$\hat{P}(\ell_{i})\equiv P(\ell_{i}|\vec{M}_{0}...\vec{M}_{i}),$ (14)

where $P(C|AB)$ is generically the probability distribution of variable $C$
given knowledge of the state of variables $A$ and $B$. Using the chain rule of
probability, Eq. (14) can be put into a recursive form as

$\begin{split}\hat{P}(\ell_{i})&=\frac{P(\vec{M}_{i}\ell_{i}|\vec{M_{0}}...\vec{M}_{i-1})}{P(\vec{M}_{i}|\vec{M}_{0}...\vec{M}_{i-1})}\\\
&=\frac{1}{P(\vec{M}_{i}|\vec{M}_{0}...\vec{M}_{i-1})}\sum^{7}_{\ell_{i-1}=0}P(\vec{M}_{i}\ell_{i}|\ell_{i-1}\vec{M_{0}}...\vec{M}_{i-1})\hat{P}(\ell_{i-1})\\\
&=\frac{1}{P(\vec{M}_{i}|\vec{M}_{0}...\vec{M}_{i-1})}\sum^{7}_{\ell_{i-1}=0}P(\vec{M}_{i}|\ell_{i}\ell_{i-1}\vec{M_{0}}...\vec{M}_{i-1})P(\ell_{i}|\ell_{i-1}\vec{M_{0}}...\vec{M}_{i-1})\hat{P}(\ell_{i-1})\\\
&=\frac{1}{P(\vec{M}_{i}|\vec{M}_{0}...\vec{M}_{i-1})}\sum^{7}_{\ell_{i-1}=0}P(\vec{M}_{i}|\ell_{i}\ell_{i-1})P(\ell_{i}|\ell_{i-1})\hat{P}(\ell_{i-1}),\end{split}$
(15)

where in the last line we explicitly assume that $\vec{M}_{i}$ and $\ell_{i}$
are conditionally independent of $\\{\vec{M}_{0},...,\vec{M}_{i-1}\\}$ given
knowledge of $\ell_{i-1}$. This is equivalent to assuming that the state
transitions depend only on the prior state (Markovian assumption), and that
the additive measurement noise is uncorrelated across time. Since
$P(\vec{M}_{i}|\vec{M}_{0}...\vec{M}_{i-1})$ has no explicit $\ell_{i}$
dependence and is not otherwise coupled to any of the terms in the sum, it
will be ignored in our subsequent analysis and treated as simply a
normalization factor.

If we assume that the system begins in a known state, which is reasonable at
the start of a quantum experiment, then Eq. (15) provides a recipe for
iteratively computing the probabilities of future states in light of new
measurement results. The two quantities that must be solved for as part of
this procedure are

  1. 1.

$P(\ell_{i}|\ell_{i-1})$, which is the probability of jumping from state
$\ell_{i-1}$ to state ${\ell_{i}}$ in time $T$ due to bit-flip errors, and

  2. 2.

$P(\vec{M}_{i}|\ell_{i}\ell_{i-1})$, which is the probability of measuring
values $\vec{M}_{i}$ from the parity operators given knowledge of the system
states at the beginning and end of the integration period.

The remainder of this section will be dedicated to deriving explicit
functional forms for these two probability terms.

###  3.2 Analysis of the transition probability

Our treatment of bit-flip errors will assume that they act on the system
independently of one another and occur at a fixed rate $\mu$ that is known to
us. In saying that the errors are independent, we mean that an error occurring
in one time interval gives us no information about whether an error will occur
in any other non-overlapping interval. The number $e_{k}$ of such errors
acting on the $k$th qubit in a time interval of length $T$ will obey the
following Poisson distribution

$P(e_{k})=\frac{(\mu T)^{e_{k}}}{e_{k}!}e^{-\mu T},$ (16)

where we emphasize that the error rate is identical across all three qubits.

In order to use Eq. (16) to derive a functional form for
$P(\ell_{i}|\ell_{i-1})$, we must specify the number of errors needed to
connect state $\ell_{i-1}$ to state $\ell_{i}$. Taking for example a system
that starts as $\ket{000}$ at step $i-1$ and ends up as $\ket{100}$ at step
$i$, corresponding to $\ell_{i-1}=0$ and $\ell_{i}=4$, it is clear that the
first qubit must have been flipped at some point in the interval between the
two steps. However, this does not mean that only a single error occurred on
the first qubit, or that none of the other qubits experienced any errors.
Rather, it means that the first qubit experienced an odd number of errors
while the rest of the qubits experienced an even number of errors (including
no error at all), such that the net number of bit-flips works out to be one
for the first qubit and zero for the others. We can therefore express
$P(\ell_{i}|\ell_{i-1})$ as a sum over the Poisson probabilities of every
possible error combination consistent with the net number of flips between
$\ell_{i}$ and $\ell_{i-1}$.

Rather than evaluate this sum directly, it is easier to look at each qubit
separately and calculate the probability that it will experience either an
even or odd number of errors. Since the errors on each qubit are independent,
we can then take appropriate products of these probabilities to calculate
$P(\ell_{i}|\ell_{i-1})$. The probability that the $k$th qubit will experience
an even number of errors is

$P(e_{k}\text{ is even})=\sum^{\infty}_{j=0}P(e_{k}=2j)=e^{-\mu
T}\sum^{\infty}_{j=0}\frac{(\mu T)^{2j}}{(2j)!}=e^{-\mu T}\cosh(\mu T),$ (17)

and therefore the probability of an odd number of errors is

$P(e_{k}\text{ is odd})=1-P(e_{k}\text{ is even})=e^{-\mu T}\sinh(\mu T).$
(18)

Using Eq. (17) and Eq. (18), the transition probability
$P(\ell_{i}|\ell_{i-1})$ is given by

$P(\ell_{i}|\ell_{i-1})=[\sinh(\mu T)]^{d(\ell_{i},\ell_{i-1})}[\cosh(\mu
T)]^{3-d(\ell_{i},\ell_{i-1})}\exp[-3\mu T],$ (19)

where $d(\ell_{i},\ell_{i-1})$ is the Hamming distance between the bit
representation of $\ell_{i}$ and the bit representation of $\ell_{i-1}$. In
words, the probability consist of an exponential multiplied by a sinh term for
every qubit that experiences a net flip and a cosh term for every qubit whose
state is ultimately left unchanged.

Since the error rate $\mu$ is independent of time, the value of index $i$ is
irrelevent to the value of $P(\ell_{i}|\ell_{i-1})$. This allows us to define
a single $8\times 8$ time-invariant transition matrix $J$ with elements
$J_{\ell_{i-1},\ell_{i}}\equiv P(\ell_{i}|\ell_{i-1})$. This matrix can be
understood as the parameterization of a discrete Markov chain [15] which
describes how bit-flip errors can alter the state of our system over time
intervals of length $T$. For the remainder of this paper we will use the
elements of $J$ to denote the transition probabilities
$P(\ell_{i}|\ell_{i-1})$.

###  3.3 Analysis of the measurement density

From our discussion of continuous measurement in Sec. 2.2, we know that the
value of the $i$th measurement will depend on the average value of the
syndrome from time $iT$ to $(i+1)T$, plus an additive Gaussian noise term.
Since it will be necessary to discuss measurements of the two parity operators
separately, we unpack the vector quantity $\vec{M}_{i}$ as per Eq. (13) into
its components $M^{(1)}_{i}$ and $M^{(2)}_{i}$, which describe the results of
measuring $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$ respectively. The measurement
distribution $P(\vec{M}_{i}|\ell_{i}\ell_{i-1})$ will be hereafter written as
$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ to make this distinction
explicit. For convenience we denote the averaged syndrome values as
$\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$, such that

$\bar{S}^{(j)}_{i}\equiv\frac{1}{T}\int^{T}_{0}S^{(j)}(iT+\tau)d\tau$ (20)

where $S^{(j)}(t)$ is the value of the $j$th syndrome at time $t$.

Using this labeling, we have from Eq. (9) that
$M^{(j)}_{i}\sim\mathcal{N}(\bar{S}^{(j)}_{i},\frac{k}{T})$, so
$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ will broadly resemble a
bivariate Gaussian distribution. However, since the means $\bar{S}^{(1)}_{i}$
and $\bar{S}^{(2)}_{i}$ are themselves random variables due to errors
occurring within the integration period, the measurement distribution will
consist of a Gaussian function integrated over mean values in the continuous
interval [-1, +1]. This is given by

$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})=\int^{1}_{-1}\int^{1}_{-1}\prod^{2}_{j=1}\frac{\exp[\frac{-T}{2k}(M^{(j)}_{i}-\bar{S}^{(j)}_{i})^{2}]}{\sqrt{2\pi\frac{k}{T}}}P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})d\bar{S}^{(1)}_{i}d\bar{S}^{(2)}_{i},$
(21)

where we assume that the measurement noise of $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$ is
uncorrelated and of equal variance. Parsing Eq. (21), we see that the
measurement distribution has the form of a continuous Gaussian mixture, with
each component centered at a different point in the 2D interval
$[-1,+1]\times[-1,+1]$ and weighted by the probability
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$ of observing those
mean values, assuming that $\ell_{i}$ and $\ell_{i-1}$ are known. Figure 2
illustrates the Gaussian-like behavior of
$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$, which is especially clear
under strong noise. We exploit this similarity in our Bayesian algorithm to
efficiently approximate the measurement log-likelihood as described in Sec.
5.1.

Figure 2: Contour plots of $P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$
derived from Eq. (21) using a uniform distribution over the syndrome means
$\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$. Plot A was generated using a
moderate amount of noise ($\frac{k}{T}=1)$, and possesses a Gaussian-like
shape. Plot B contains far less noise ($\frac{k}{T}=0.05)$, so the square
shape of the underlying uniform distribution starts to show through.

The fact that Eq. (21) is not simply a bivariate Gaussian distribution with a
fixed mean reflects the uncertainty that is inherent in our problem. Since
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$ is conditioned only
on the state of the system at $iT$ and $(i+1)T$, we must predict the values of
the syndrome averages using this information alone. Figure 1C demonstrates the
difficulty of such an inference, as the value of $\bar{S}^{(j)}_{i}$ is
dependent not only on the precise number of errors in the integration interval
but more importantly on the locations of these errors within the interval, all
of which are unknown to us. Altogether there are an infinite number of
possible values that $\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$ can take for
any given $(\ell_{i-1},\ell_{i})$ pair, so we must marginalize over these
degrees of freedom in order to derive the measurement distribution.

To better understand where the probability of $\bar{S}^{(1)}_{i}$ and
$\bar{S}^{(2)}_{i}$ is concentrated, we can condition
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$ on the number of
errors occurring in the $i$th interval, which gives

$\begin{split}P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})&=\sum^{\infty}_{e_{1},e_{2},e_{3}=0}P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1};e_{1}e_{2}e_{3})P(e_{1}e_{2}e_{3}|\ell_{i}\ell_{i-1})\\\
&=\sum^{\infty}_{e_{1},e_{2},e_{3}=0}P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})P(e_{1}|\ell_{i}\ell_{i-1})P(e_{2}|\ell_{i}\ell_{i-1})P(e_{3}|\ell_{i}\ell_{i-1}),\end{split}$
(22)

where $e_{k}$ is the number of errors occurring on the $k$th qubit during the
$i$th interval (this time index is suppressed for notational convenience). In
the final line we exploit the fact that errors on different qubits are
independent, and that $\ell_{i}$ is a deterministic function of $e_{1}$,
$e_{2}$, $e_{3}$, $\ell_{i-i}$ and therefore does not need to be explicitly
conditioned on.

From our discussion in Sec. 3.2, we know that $\ell_{i}$ and $\ell_{i-1}$ only
specify whether an even or odd number of errors occurred on each qubit, so
$P(e_{k}|\ell_{i}\ell_{i-1})$ can be rewritten as $P(e_{k}|\text{odd})$ or
$P(e_{k}|\text{even})$ depending on whether $\ell_{i}$ and $\ell_{i-1}$ differ
in the $k$th qubit. These distributions are identical to the Poisson
distribution from Eq. (16) except that the probability of either even or odd
$e_{k}$ is set to zero. After re-normalizing, we get

$P(e_{k}|\text{odd})=\begin{cases}\frac{(\mu T)^{e_{k}}}{e_{k}!\sinh(\mu
T)}&e_{k}\text{ odd}\\\ 0&e_{k}\text{ even}\end{cases},\quad\quad\quad
P(e_{k}|\text{even})=\begin{cases}\frac{(\mu T)^{e_{k}}}{e_{k}!\cosh(\mu
T)}&e_{k}\text{ even}\\\ 0&e_{k}\text{ odd}\end{cases}.$ (23)

As an example of how these distributions are used, if $\ell_{i}=4$ and
$\ell_{i-1}=0$ then the error distribution factorizes as
$P(e_{1}e_{2}e_{3}|4,0)=P(e_{1}|\text{odd})P(e_{2}|\text{even})P(e_{3}|\text{even})$.

Referring back to Eq. (22), the final term left to characterize is
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})$, which
describes how the syndrome averages are distributed given a starting state and
the number of errors that occurred. This term is challenging to evaluate,
since errors on the second qubit will flip both syndromes simultaneously and
thus prevent the joint distribution from factorizing. If we consider errors on
the second qubit separately from errors on the first and third qubits, i.e.
focus on cases where $e_{1},e_{3}=0$ and where $e_{2}=0$, then we can derive
(see Appendix A.1) the following expressions

$\displaystyle P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}0e_{3})$
$\displaystyle=P(\bar{S}^{(1)}_{i}|\ell_{i-1};e_{1}00)P(\bar{S}^{(2)}_{i}|\ell_{i-1};00e_{3})$
(24) $\displaystyle P(\bar{S}^{(1)}_{i}|\ell_{i-1};e_{1}00)$
$\displaystyle=\begin{cases}\delta(1\mp_{1}\bar{S}^{(1)}_{i})&e_{1}=0\\\\[2.84544pt]
\dfrac{e_{1}!}{2a_{1}!b_{1}!}\left(\dfrac{1\pm_{1}\bar{S}^{(1)}_{i}}{2}\right)^{a_{1}}\left(\dfrac{1\mp_{1}\bar{S}^{(1)}_{i}}{2}\right)^{b_{1}}&e_{1}>0\end{cases}$
(25) $\displaystyle P(\bar{S}^{(2)}_{i}|\ell_{i-1};00e_{3})$
$\displaystyle=\begin{cases}\delta(1\mp_{2}\bar{S}^{(2)}_{i})&e_{3}=0\\\\[2.84544pt]
\dfrac{e_{3}!}{2a_{3}!b_{3}!}\left(\dfrac{1\pm_{1}\bar{S}^{(1)}_{i}}{2}\right)^{a_{3}}\left(\dfrac{1\mp_{1}\bar{S}^{(1)}_{i}}{2}\right)^{b_{3}}&e_{3}>0\end{cases}$
(26) $\displaystyle P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};0e_{2}0)$
$\displaystyle=\begin{cases}\delta(1\mp_{1}\bar{S}^{(1)}_{i})\delta(\pm_{1}\bar{S}^{(1)}_{i}\mp_{2}\bar{S}^{(2)}_{i})&e_{2}=0\\\\[2.84544pt]
\dfrac{e_{2}!}{2a_{2}!b_{2}!}\left(\dfrac{1\pm_{2}\bar{S}^{(1)}_{i}}{2}\right)^{a_{2}}\left(\dfrac{1\mp_{2}\bar{S}^{(1)}_{i}}{2}\right)^{b_{2}}\delta(\pm_{1}\bar{S}^{(1)}_{i}\mp_{2}\bar{S}^{(2)}_{i})&e_{2}>0\end{cases}$
(27)

where $a_{k}$ and $b_{k}$ are the integer parts of $\frac{e_{k}}{2}$ and
$\frac{e_{k}-1}{2}$ respectively, and $\delta(x)$ is the Dirac delta function.
The “$\pm_{j}$” term indicates whether state $\ell_{i-1}$ has even “$+$” or
odd “$-$” parity with respect to the $j$th measurement operator, and therefore
whether the syndrome should be added or subtracted. We note that the syndrome
density is uniform (on its support) when only a single error occurs (i.e.
$a_{k}=b_{k}=0$). This simplification will be utilized in Sec. 5.1 as we
construct our logarithmic filter.

Unfortunately, when errors occur on the second qubit and on at least one of
the other qubits, the joint distributions are non-smooth and appear to lack
convenient analytic forms (see Appendix A.1.3). As an example, if the system
starts in $\ket{000}$ and experiences one error on each qubit, then the
syndrome distribution is given by

$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|0;111)=\frac{1}{4}\left[\min(\bar{S}^{(1)}_{i},\bar{S}^{(2)}_{i})+\max(0,\bar{S}^{(1)}_{i}+\bar{S}^{(2)}_{i})\right],$
(28)

which is a piecewise function with a discontinuous first derivative. Figure 3
shows a plot of $P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|0;111)$, together with
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|0;011)$ which is even less smooth.

Due to the poor behavior of the syndrome densities, we are unable to make
further progress toward an exact, analytic expression for the measurement
probability $P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$. Most of this
difficulty stems from the fact that we must, in general, account for an
infinite number of possible errors that could occur across all three qubits
within a given integration interval. If, however, we restrict ourselves to
situations where only a single error occurs, which is reasonable when $\mu$
and $T$ are small, then the expressions in Eqs. (25 – 27) can be used directly
and we will be able to evaluate the integrals in Eq. (21). This is the
approach we take in Sec. 5.1, although it is no longer an exact treatment of
the problem.

Figure 3: Plot A shows the probability density
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|0;011)$, which describes the syndrome
means when the system starts in $\ket{000}$ and then experiences an error on
the second and third qubits. Plot B shows
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|0;111)$, whose form is given in Eq.
(28). As a general rule, the more errors that occur in the interval the more
“smooth” the syndrome density will appear, with plot A being outright
discontinuous while plot B has a discontinuous first derivative. These plots
were generated from $10^{8}$ samples using histograms with 50 equally spaced
bins along each axis. The color gradient indicates the relative magnitude of
each bin, with different scales for the two plots.

###  3.4 Constructing the Bayesian Filter

Although Sec. 3.3 suggests that the state tracking problem does not possess a
convenient functional form, it is still useful to synthesize our results from
this section into an algebraic description of how the optimal model would
operate. Referring back to Eq. (15), the posterior state probabilities
$\hat{P}(\ell_{i})$ are given recursively in terms of the measurement
probabilities $P({M}^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ and transition
probabilities $J_{\ell_{i-1}\ell_{i}}\equiv P(\ell_{i}|\ell_{i-1})$. Noting
the dependence of both terms on $\ell_{i}$ and $\ell_{i-1}$, the update rule
is given by

$\begin{split}\hat{P}(\ell_{i})&\propto\sum^{7}_{\ell_{i-1}=0}P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})J_{\ell_{i-1}\ell_{i}}\hat{P}(\ell_{i-1})\\\
&\propto\sum^{7}_{\ell_{i-1}=0}\hat{P}(\ell_{i-1})\tilde{J}(M^{(1)}_{i}M^{(2)}_{i})_{\ell_{i-1}\ell_{i}},\end{split}$
(29)

where $\tilde{J}(M^{(1)}_{i}M^{(2)}_{i})_{\ell_{i-1}\ell_{i}}\equiv
P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})J_{\ell_{i-1},\ell_{i}}$ can be
viewed as a weighted transition matrix which combines the Markov dynamics of
the error process with the likelihood of the observed measurement outcomes.

Due to the recursive nature of Eq. (29), the state tracking model operates
naturally as a digital filter, taking in syndrome measurements at each time-
step and outputting a set of posterior probabilities that incorporate all of
the information available to us. Setting aside the challenges of computing
$P({M}^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$, this filter is optimal in
the sense that it was constructed from the formal probability manipulations in
Eq. (15), which assumed only that the measurements and errors were Markovian.
In Sec. 3.3 we made further assumptions about the structure of the measurement
signal, and insofar as these assumptions are valid there is nowhere for the
algorithm to be improved. Of course, should the physical system not obey the
idealized model outlined in Sec. 2.2, then the filter will need to be modified
in order to remain optimal.

From a computational standpoint, Eq. (29) represents the filter as a simple
matrix-vector product, with the difficulty centered on calculating
$\tilde{J}$. More specifically, Sec. 3.3 made it clear that computing
$P({M}^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ requires us to evaluate the
integral in Eq. (21), which appears to be analytically intractable when the
error combinations are not restricted. This limitation, together with other
issues that will be discussed in Sec. 4, rules out using the filter for real-
time error tracking. That said, accurate numerical evaluation of the integral
is still possible via sampling techniques, which are discussed further in
Appendix A.2, so we will use this optimal filter as a benchmark for other
state tracking methods.

##  4 Obstacles and Prior Work

In the context of real-time quantum error tracking, where latency is measured
in nanoseconds and hardware resources can be tight, the Bayesian filter
outlined in Sec. 3 has four significant challenges:

  1. 1.

The Gaussian integral in Eq. (21) cannot be evaluated analytically, and
sampling methods are too resource-intensive for real-time filtering.

  2. 2.

Computing even a single Gaussian requires the implementation of an exponential
function, which can be challenging on low-latency, high-throughput devices
such as FPGAs [16].

  3. 3.

The outputs of the Gaussians will likely need to be represented using
floating-point numbers to capture the necessary range and precision, which
adds further computational burden.

  4. 4.

The probabilities generated by Eq. (29) will need to be periodically
normalized in order to prevent overflow or underflow, which requires a
dedicated division routine that will take up resources and be slow to run.

Overcoming all of these obstacles while preserving the underlying Bayesian
framework is not trivial.

The first derivation of an optimal continuous-time filter for the three qubit
bit-flip code was published by van Handel and Mabuchi, who recovered the well-
known Wonham filter after solving for the least-squares estimator of the
density matrix [17]. This classical filter was designed to output
probabilities for the states of a Markov chain observed continuously by a
signal with additive Gaussian noise. Using a modified version of our notation,
the canonical Wonham filter has the form

$d\hat{P}_{\ell}(t)=\sum^{7}_{\ell^{\prime}=0}\hat{P}_{\ell^{\prime}}(t)Q_{\ell^{\prime}\ell}dt+\frac{1}{k}\sum^{2}_{j=1}[S^{(j)}_{\ell}-\bar{\bar{S}}^{(j)}(t)]\hat{P}_{\ell}(t)[dM^{(j)}(t)-\bar{\bar{S}}^{(j)}(t)dt],$
(30)

where $\hat{P}_{\ell}(t)$ is the probability of the $\ell$th state at time
$t$, $S^{(j)}_{\ell}\in\\{-1,+1\\}$ is the syndrome of the $\ell$th state from
the $j$th parity operator,
$\bar{\bar{S}}^{(j)}(t)\equiv\sum_{\ell}\hat{P}_{\ell}(t)S^{(j)}_{\ell}$ is
the average of the $j$th syndrome at time $t$, and $Q=\frac{1}{T}\ln J$ is the
rate matrix of the continuous Markov chain which describes the state
transitions.

It is important to emphasize that Eq. (29) and Eq. (30) are optimal with
respect to different measurement regimes. The Wonham filter is optimal under
the assumption that we have access to instantaneous signal readouts at
arbitrary $t$, while our filter is optimal when the measurements are
restricted to being integrated averages of the signal over some finite period
$T$ as in Eq. (9). As $T\rightarrow 0$ these two types of measurement
converge, and thus Eq. (29) converges to Eq. (30) after normalization (see
Appendix A.3). The general problem of filtering Markov jump processes using
discrete observations is explored by Borisov [18], who derives results that
are in agreement with our work in Sec. 3.

As a non-linear stochastic differential equation, the Wonham filter is not
practical to use in its exact form, though it is still possible to create
discretized approximations of Eq. (30) which are effective. For example, Gange
George et al. applied the Euler-Maruyama method to a logarithmic
transformation of the Wonham filter in order to derive a first-order update
rule that was numerically stable [19]. In the specific context of continous
error correction, Mohseninia et al. [5] proposed a linearized version of Eq.
(30) equivalent to

$\hat{P}(\ell_{i})\propto\hat{P}(\ell_{i-1}=\ell_{i})+T\sum^{7}_{\ell_{i-1}=0}\left[Q_{\ell_{i-1}\ell_{i}}+\frac{M^{(1)}_{i}S^{(1)}_{\ell_{i-1}}+M^{(2)}_{i}S^{(2)}_{\ell_{i-1}}}{k}\delta_{\ell_{i-1}\ell_{i}}\right]\hat{P}({\ell_{i-1}}),$
(31)

where $S^{(j)}_{\ell_{i-1}}$ and $Q$ are defined as in Eq. (30). The filter
described by Eq. (31) does not involve any Gaussian functions and therefore
avoids the first three obstacles in our list, but will still require periodic
normalization. Additionally, as $T$ grows larger the accuracy of this first-
order approximation will decline.

The filters discussed so far have all been motivated at some level by Bayesian
probability analysis, but there exists another class of algorithm, referred to
as threshold or boxcar filters, which eschews connections to probability
theory in favor of computational simplicity [20]. At a basic level, these
models take a pair of integrated measurement values
$(M^{(1)}_{i},M^{(2)}_{i})$ and compares them to a set of predefined
thresholds, after which an appropriate action is taken. Atalaya et al.
developed a double-threshold algorithm for the Bacond-Shor code using non-
commuting observables [4], while Mosheninia et al. proposed a non-Markovian
boxcar filter and a double threshold scheme for the three-qubit repetition
code which were easy to implement and attained good performance [5]. Atalaya
and Zhang et al. applied a flexible double thresholding scheme to a system
undergoing quantum annealing, demonstrating that these algorithms were
effective at error correction even in the presence of Hamiltonian evolution
[6].

##  5 A Practical Bayesian Filter

To address the obstacles identified in Sec. 4, we propose an algorithm that
avoids any exponentiation or division operations, while also approximating the
analytically intractable integral in Eq. (21). In Sec. 5.1 we simplify Eq.
(21) by assuming that only a single error occurs during each integration
period, which yields simple expressions for
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$. In Sec. 5.2 we
eliminate the need for exponentiation and division by shifting our analysis
into log-probability space, where the multiplication of probabilities becomes
addition, normalization becomes subtraction, and Gaussian distributions
transform into simple quadratic functions. Finally, Sec. 5.4 introduces what
we call “single-term” and “two-term” strategies for approximating the
LogSumExp functions needed for Markov chain evolution using log-probabilities.
The performance of our algorithm is tested numerically in Sec. 6.

###  5.1 Single-error approximation

The largest obstacle to using Eq. (29) as a filter is that
$P(M^{(1)}M^{(2)}|\ell_{i}\ell_{i-1})$ cannot be easily evaluated. As explored
in Sec. 3.3, this difficulty arises because the syndrome density terms
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})$ do not in
general have convenient analytic forms, especially when the number of errors
is large. This in turn prevents us from solving the integral in Eq. (21).
However, if we were to assume that only a single error can occur within a
given measurement interval, then Eqs. (25 – 27) would reduce to

$\displaystyle P(\bar{S}^{(1)}_{i}|\ell_{i-1};e_{1}00)$
$\displaystyle=\begin{cases}\delta(1\mp_{1}\bar{S}^{(1)})&e_{1}=0\\\\[2.84544pt]
\frac{1}{2}&e_{1}=1\end{cases}$ (32) $\displaystyle
P(\bar{S}^{(1)}_{i}|\ell_{i-1};00e_{3})$
$\displaystyle=\begin{cases}\delta(1\mp_{2}\bar{S}^{(2)}_{i})&e_{3}=0\\\\[2.84544pt]
\frac{1}{2}&e_{3}=1\end{cases}$ (33) $\displaystyle
P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};0e_{2}0)$
$\displaystyle=\begin{cases}\delta(1\mp_{1}\bar{S}^{(1)}_{i})\delta(1\mp_{2}\bar{S}^{(2)}_{i})&e_{2}=0\\\\[2.84544pt]
\delta(\pm_{1}\bar{S}^{(1)}_{i}\mp_{2}\bar{S}^{(2)}_{i})&e_{2}=1,\end{cases}$
(34)

which all lead to tractable integrals when substituted into Eq. (21). This
single-error approximation, which is reasonable when the average number of
error per interval $\mu T$ is small, leads to Gaussian-like measurement
distributions which can be easily incorporated into our log-probability
filter.

For simplicity, we first consider an error that occurs on either the first
($e_{1}=1)$ or third ($e_{3}=1)$ qubit, since it will affect only one of the
syndromes and thus allow the joint distribution to separate as in Eq. (24). By
substituting $P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};100)$ into Eq.
(21) and requiring that $\ell_{i}$ and $\ell_{i-1}$ differ only in the first
qubit, we get

$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};100)=\frac{1}{4}\left[\text{erf}\left(\frac{M^{(1)}_{i}+1}{\frac{2k}{T}}\right)-\text{erf}\left(\frac{M^{(1)}_{i}-1}{\frac{2k}{T}}\right)\right]\cdot\frac{\exp[\frac{-T}{2k}(M^{(2)}_{i}\mp_{2}1)^{2}]}{\sqrt{2\pi\frac{k}{T}}},$
(35)

where $\text{erf}(x)\equiv\frac{2}{\sqrt{\pi}}\int^{x}_{0}e^{-y^{2}}dy$ is the
error function. The last line of Eq. (35) consists of two distinct terms, one
of which is a function of $M^{(1)}_{i}$ and the other of $M^{(2)}_{i}$. The
$M^{(2)}_{i}$ term is simply a Gaussian centered at $\pm_{2}1$, since the
second syndrome is not sensitive to an error on the first qubit. The
$M^{(1)}_{i}$ term, by contrast, arises because the mean of the first syndrome
is distributed uniformly over the range $[-1,+1]$, which yields a pair of
error functions. Notably, this second term does not depend on $\pm_{1}$, so it
is independent of the state of the first two qubits.

While the second term in Eq. (35) is an exact result from of our single-error
integration, it is not very convenient to evaluate numerically. As described
in Appendix A.4, the measurement distribution can be approximated as

$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};100)=\frac{\exp[\frac{-3}{2(1+3\frac{k}{T})}(M^{(1)}_{i})^{2}]}{\sqrt{2\pi(\frac{1}{3}+\frac{k}{T}})}\cdot\frac{\exp[\frac{-T}{2k}(M^{(2)}_{i}\mp_{2}1)^{2}]}{\sqrt{2\pi\frac{k}{T}}},$
(36)

which will be easy to evaluate numerically after we shift to log-probabilities
in Sec. 5.2. The measurement density for an error on the third qubit can be
attained from Eq. (36) by simply swapping $M^{(1)}_{i}$ and $M^{(2)}_{i}$ and
replacing $\mp_{2}$ with $\mp_{1}$. Carrying out a similar procedure on the
second qubit gives

$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};010)=\frac{1}{2}\frac{\exp[\frac{-T}{k}(\frac{M^{(1)}_{i}\mp_{c}M^{(2)}_{i}}{2})^{2}]}{\sqrt{\pi\frac{k}{T}}}\cdot\frac{\exp[\frac{-3}{2(1+3\frac{k}{2T})}(\frac{M^{(1)}_{i}\pm_{c}M^{(2)}_{i}}{2})^{2}]}{\sqrt{2\pi(\frac{1}{3}+\frac{k}{2T}})},$
(37)

where $\pm_{c}$ is positive when $\ell_{i-1}$ has the same parity with respect
to $Z_{1}Z_{2}$ and $Z_{2}Z_{3}$ but negative when the parities are opposite
(roughly, “$\pm_{c}\equiv\pm_{1}\cdot\pm_{2}$”).

The final case to consider is when no errors occur in an interval at all. With
reasonable values of $\mu$ and $T$ this is the most likely outcome for any
given interval, and it’s measurement distribution can be solved for easily by
substituting
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})=P(\bar{S}^{(1)}_{i}|\ell_{i-1};000)P(\bar{S}^{(2)}_{i}|\ell_{i-1};000)$
into Eq. (21). This gives

$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};000)=\frac{\exp[\frac{-T}{2k}(M^{(1)}_{i}\mp_{1}1)^{2}]}{\sqrt{2\pi\frac{k}{T}}}\cdot\frac{\exp[\frac{-T}{2k}(M^{(2)}_{i}\mp_{2}1)^{2}]}{\sqrt{2\pi\frac{k}{T}}},$
(38)

which is simply the product of Gaussian distributions with means corresponding
to the parities of $\ell_{i-1}$.

Taken together, Eqs. (36, 37, 38) constitute the approximate description of
$P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ that we will use to construct
our log-probability filter. The Gaussian form of each equation is especially
important on a practical level, as the log-densities will all reduce to sums
of simple quadratic equations which can be easily computed.

###  5.2 Moving to logarithmic probability

To begin the transformation from probabilities to log-probabilities, we take
the logarithm of the optimal Bayesian filter from Eq. (29)

$\begin{split}\log\hat{P}(\ell_{i})&=\log\sum^{7}_{\ell_{i-1}=0}P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})J_{\ell_{i-1}\ell_{i}}\hat{P}(\ell_{i-1})\\\
&=\log\sum^{7}_{\ell_{i-1}=0}\exp[\log\hat{P}(\ell_{i-1})+\log\tilde{J}(M^{(1)}_{i}M^{(2)}_{i})_{\ell_{i-1}\ell_{i}}],\end{split}$
(39)

where in the last line we exponentiate the logarithm of each term in the sum
in order to preserve the recursive structure of the filter. Since
$\log\tilde{J}\equiv\log P(\ell_{i}|\ell_{i-1})+\log
P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$, the terms in Eq. (39) can be
evaluated using our results from Sec. 3.2 and Sec. 5.1.

The transition log-probability is given by the logarithm of Eq. (19), which
yields

$\log J_{\ell_{i-1}\ell_{i}}=d(\ell_{i},\ell_{i-1})\log\sinh(\mu
T)+[3-d(\ell_{i},\ell_{i-1})]\log\cosh(\mu T)-3\mu T$ (40)

where $d(\ell_{i},\ell_{i-1})$ is the Hamming distance between $\ell_{i}$ and
$\ell_{i-1}$. While the form of Eq. (40) is not especially illuminating, the
value of $\mu T$ will be known to us and assumed to be constant across an
experiment. We are therefore able to compute the values of $\log
J_{\ell_{i-1}\ell_{i}}$ in advance and use them in the filter without needing
to evaluate the $\log\sinh$ or $\log\cosh$ functions in real time.

The measurement log-probabilities, by contrast, must be calculated anew each
time a measurement is recorded. Such rapid computation is feasible because
Eqs. (36, 37, 38) are all in Gaussian form, and therefore have the following
quadratic log-probabilities

$\displaystyle\log P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};000)$
$\displaystyle=-\frac{T}{2k}[(M^{(1)}_{i}\mp_{1}1)^{2}-(M^{(2)}_{i}\mp_{2}1)^{2}]+\log
N$ (41) $\displaystyle\log P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};100)$
$\displaystyle=-\frac{3}{2(1+3\frac{k}{T})}(M^{(1)}_{i})^{2}-\frac{T}{2k}(M^{(2)}_{i}\mp_{2}1)^{2}+\log
N^{\prime}$ (42) $\displaystyle\log P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};010)$
$\displaystyle=-\frac{T}{k}(\frac{M^{(1)}_{i}\mp_{c}M^{(2)}_{i}}{2})^{2}-\frac{3}{2(1+3\frac{k}{2T})}(\frac{M^{(1)}_{i}\pm_{c}M^{(2)}_{i}}{2})^{2}+\log
N^{\prime\prime}$ (43) $\displaystyle\log
P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i-1};001)$
$\displaystyle=-\frac{T}{2k}(M^{(1)}_{i}\mp_{1}1)^{2}-\frac{3}{2(1+3\frac{k}{T})}(M^{(2)}_{i})^{2}+\log
N^{\prime\prime\prime},$ (44)

where the normalization constants have been left unspecified for convenience.
The different variances, all functions of $k$ and $T$, are constant during an
experiment and therefore do not need to be computed at each step.

Using Eqs. (40 – 44), we can calculate $\log\tilde{J}$ and therefore evaluate
the exponent of each term in the sum of Eq. (39). The last line of Eq. (39) is
known as a LogSumExp [21], and its evaluation is the final obstacle to
implementing our filter. For convenience we adopt the more compact notation
$L^{\ell_{i}}_{\ell_{i-1}}\equiv\log\tilde{J}_{\ell_{i-1}\ell_{i}}+\log\hat{P}(\ell_{i-1})$,
which serve as the inputs of the LogSumExp function,

$\log\hat{P}(\ell_{i})=\log\sum^{7}_{\ell_{i-1}=0}\exp[L^{\ell_{i}}_{\ell_{i-1}}].$
(45)

Using this form, we exploit a well-known trick for evaluating LogSumExp
functions by pulling the largest exponential outside of the sum

$\log\sum^{7}_{\ell_{i-1}=0}\exp[L^{\ell_{i}}_{\ell_{i-1}}]=L^{\ell_{i}}_{*}+\log(1+\sum_{\ell_{i-1}\neq*}\exp[L^{\ell_{i}}_{\ell_{i-1}}-L^{\ell_{i}}_{*}]),$
(46)

where $L^{\ell_{i}}_{*}$ is the largest input to the LogSumExp function. The
second term in Eq. (46) is now the logarithm of one plus a set of terms whose
exponents are all non-positive. Figure 4 shows that as the difference between
$L^{\ell_{i}}_{*}$ and the other inputs grows, the output of Eq. (46) will
converge to $L^{\ell_{i}}_{*}$. If the difference is large for some subset of
the inputs, the expression can be simplified by removing that subset from the
sum.

Figure 4: Plot of $\log(1+e^{\Delta L})$, which goes to zero as the magnitude
of $\Delta L\equiv L^{\ell_{i}}_{\ell_{i-1}}-L^{\ell_{i}}_{*}$ increases. This
implies that $\text{LogSumExp}(x)\equiv L^{\ell_{i}}_{*}+\log(1+e^{\Delta L})$
rapidly approaches $L^{\ell_{i}}_{*}$ as the gap between $L^{\ell_{i}}_{*}$
and the other terms increases. Note that $\Delta L$ is never positive, since
$L^{\ell_{i}}_{*}$ is the largest of all $L^{\ell_{i}}_{\ell_{i-1}}$ by
definition.

###  5.3 Evaluating the LogSumExp

To take advantage of Eq. (46), we must know the relative magnitudes of the
various $L^{\ell_{i}}_{\ell_{i-1}}$ terms in a typical run. For convenience,
we introduce the non-positive quantity $\Delta
L^{\ell_{i}}_{\ell_{i-1}}=L^{\ell_{i}}_{\ell_{i-1}}-L^{\ell_{i}}_{*}$, which
is simplified to just $\Delta L$ when the specific indices are not relevant.
Using Figure 4 as a guide, it is clear that when

$\Delta L\leq-4$ (47)

the magnitude of $\log(1+e^{\Delta L})$ is negligible. This provides a
threshold to determine whether a given $L^{\ell_{i}}_{\ell_{i-1}}$ contributes
significantly to Eq. (46) or if it can instead be reasonably ignored.

Figure 5: Plots of average $L^{\ell_{i}}_{\ell_{i-1}}$ contributions to the
posteriors of $\ket{000}$ (plot A) and of $\ket{001}$ (plot B) taken from the
optimal filter during a 15 $\mu s$ run with $\frac{k}{T}=4$ that was
initialized into $\ket{000}$. The dashed black lines mark the occurrence of a
bit-flip error at $7.5\ \mu$s on the third qubit, which takes the system to
$\ket{001}$. The plotted data represent averages over the measurement noise
taken from $10^{5}$ different runs. When the system is $\ket{000}$, plot A
shows clearly that only $L^{0}_{0}$ is significant, while plot B shows that
both $L^{1}_{0}$ and $L^{1}_{1}$ are significant. After the bit-flip occurs
the situation is reversed, with $L^{1}_{1}$, $L^{0}_{0}$, and $L^{0}_{1}$
contributing significantly.

In Figure 5, we plot the average values of $L^{\ell_{i}}_{\ell_{i-1}}$ for
$\ell_{i}=0$ and $\ell_{i}=1$ before and after a bit-flip error occurs. These
values are taken from the optimal filter in Eq. (29), which we simulate
numerically using methods described in Appendix A.2. The plots of Figure 5
demonstrate clearly that only one $L^{\ell_{i}}_{\ell_{i-1}}$ term is
significant when $\ell_{i}$ is equal to the true system state and two of the
terms are significant when $\ell_{i}$ is separated from the true state by one
bit-flip. Plots of the other $\ell_{i}$ indicate that states separated from
the true state by two bit-flips have four significant terms and the state
separated by three bit-flips has eight significant terms.

To understand the behavior shown in Figure 5, it is necessary to determine the
typical magnitudes of Eq. (40) and Eqs.(41 – 44) under realistic conditions.
Starting first with the measurement log-densities, we are interested in their
expected values when the system is in some state $\ell^{\prime}$, so we take
the average of Eqs.(41 – 44) with respect to measurements distributed as
$P(M^{(1)}_{i}M^{(2)}_{i}|\ell^{\prime}\ell^{\prime})$. For convenience we
label this quantity as

$G^{\ell_{i}}_{\ell_{i-1}}\equiv E[\log
P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})],$ (48)

and after some algebra the average log-densities are shown to scale as

$G^{\ell_{i}}_{\ell_{i-1}}=\begin{cases}\mathcal{O}(-\log\frac{T}{k})&\frac{T}{k}\ll
1,\ \text{all }\ell_{i},\ell_{i-1}\\\
\mathcal{O}(-\log\frac{T}{k})&\frac{T}{k}\gg 1,\ \ell_{i}\text{ or
}\ell_{i-1}\in\\{\ell^{\prime},\ell^{\prime}\oplus 7\\}\\\
\mathcal{O}(-\frac{T}{k})&\frac{T}{k}\gg 1,\ \ell_{i}\text{ or
}\ell_{i-1}\notin\\{\ell^{\prime},\ell^{\prime}\oplus 7\\},\end{cases}$ (49)

where $\oplus$ is the bitwise exclusive-or and $\frac{T}{k}$ is the precision
of the noise (reciprocal of its variance). In words, Eq. (49) states that the
log-densities all have a similar scale when the noise is large
($\frac{T}{k}\ll 1$) but begin to diverge when the noise is small
($\frac{T}{k}\gg 1$), with transitions to/from $\ell^{\prime}$ and its
complement having a disproportionately greater likelihood. This is expected,
as the measurement densities differ primarily in the locations of their mean
values, which become more obvious as the widths of the distributions shrink.
With respect to what is shown in Figure 5, Eq. (49) indicates that unless the
noise magnitude is very small, the $L^{\ell_{i}}_{\ell_{i-1}}$ will not differ
greatly from one another based on their $G^{\ell_{i}}_{\ell_{i-1}}$ values.

Turning our attention to the transition probabilities, we know from the
numerical tests in Sec. 6.2 that when $\mu T$ is on the order of $10^{-2}$ it
becomes impossible to accurately track the state of the system at moderate run
times, even for the optimal filter. We can therefore expand Eq. (40) for small
$\mu T$ as

$\log J_{\ell_{i-1}\ell_{i}}=d(\ell_{i},\ell_{i-1})\log(\mu T)+\mathcal{O}(\mu
T),$ (50)

and keep only the leading $\mu T$ term. For fixed $\mu T$, the values of $\log
J_{\ell_{i-1}\ell_{i}}$ differ only in how many factors of $\log(\mu T)$ are
included based on the Hamming distance $d(\ell_{i},\ell_{i-1})$. Since
$\log(10^{-2})\approx-4.6$, the magnitude of a given
$L^{\ell_{i}}_{\ell_{i-1}}$ decreases rapidly with every bit-flip that
separates $\ell_{i}$ from $\ell_{i-1}$. Indeed, $\log(10^{-2})$ lies below the
threshold identified in Eq. (47), so if two $L^{\ell_{i}}_{\ell_{i-1}}$ terms
differ by this amount then only the larger of the two will be significant.

Using Eq. (50) and the fact that $L^{\ell_{i}}_{\ell_{i-1}}\equiv
G^{\ell_{i}}_{\ell_{i-1}}+\log
J_{\ell_{i-1}\ell_{i}}+\log\hat{P}(\ell_{i-1})$, it is possible to explain the
broad patterns observed in Figure 5. Beginning with the contributions to
$\ket{000}$, the magnitude of $L^{0}_{0}$ (blue line in plot A) will dominate,
since the system begins in $\ket{000}$ (large prior) and has a high
probability of remaining in $\ket{000}$ (large transition element). All other
$L^{0}_{i-1}$ will involve at least one bit-flip and therefore will not be
significant relative to $L^{0}_{0}$. By contrast, among the contributions to
$\ket{001}$ (plot B), both $L^{1}_{0}$ and $L^{1}_{1}$ (blue and orange lines
respectively) are significant, as each term has one large component (prior for
$L^{1}_{0}$ and transition element for $L^{1}_{1}$) and one small component
(transition element for $L^{1}_{0}$ and prior for $L^{1}_{1}$). After the bit-
flip occurs the measurement likelihoods begin to favor $\ket{001}$ over
$\ket{000}$, so the prior terms all shift toward $\ket{001}$ as well.

###  5.4 Single-term and two-term filters

The combined results of Sec. 5.2 and Sec. 5.3 offer a straightforward recipe
for implementing the log-probability filter of Eq. (39). First, the
measurement values are passed through the quadratic expressions in Eqs.
(41–44) to compute the measurement log-densities $\log
P(M^{(1)}_{1}M^{(2)}|\ell_{i}\ell_{i-1})$, which are then added to the
transition log-probabilities $\log J_{\ell_{i-1}\ell_{i}}$ and log-priors
$\log\hat{P}(\ell_{i-1})$ to generate the $L^{\ell_{i}}_{\ell_{i-1}}$ terms.
These terms must then be fed into Eq. (46), but we are free to choose how many
of them to include based on their relative magnitudes. It is here that the
algorithm splits into two different paths depending on whether we favor
accuracy or simplicity.

If we desire a filter that is highly accurate, then we can take the two
largest values of $L^{\ell_{i}}_{\ell_{i-1}}$ for each $\ell_{i}$ and
substitute them into Eq. (46). In this scenario the LogSumExp reduces to

$\log\hat{P}(\ell_{i})=\log\sum_{\text{pair}}\exp[L^{\ell_{i}}_{\ell_{i-1}}]=L^{\ell_{i}}_{*}+\log(1+\exp[\Delta
L]),$ (51)

where $\Delta L$ is the (negative) difference between the largest and second-
largest $L^{\ell_{i}}_{\ell_{i-1}}$ terms. While Eq. (51) can be inconvenient
to evaluate exactly, methods have been developed in the context of logarithmic
number systems that utilize lookup tables and various interpolation schemes to
yield efficient and accurate estimates [22]. By keeping the two largest values
of $L^{\ell_{i}}_{\ell_{i-1}}$ we guarantee that the log-probabilities of the
true state and its three nearest states are all correctly propagated into the
next time step. The log-probabilities of the other four states will not be
computed as accurately, although this is unimportant since these states are
not involved in the single-flip transitions which the three-qubit code is
designed to detect.

Although the logarithm term in Eq. (51) can be evaluated in a reasonable
manner, it still represents an extra computational step that increases the
overall complexity of the algorithm. If simplicity is valued over accuracy,
then we can alternatively choose to keep only the largest of the
$L^{\ell_{i}}_{\ell_{i-1}}$ after each time step. In this extreme case the
LogSumExp collapses to

$\log\hat{P}(\ell_{i})=\log\exp[L^{\ell_{i}}_{*}]=L^{\ell_{i}}_{*},$ (52)

and $L^{\ell_{i}}_{*}$ simply becomes the new posterior. Figure 5 shows that
when $\ell_{i}$ equals the true state only one $L^{\ell_{i}}_{\ell_{i-1}}$
will be significant, so nothing important is being ignored there. For states
adjacent to the true state, however, we will necessarily be ignoring one of
the two significant contributions. While this single-term simplification
appears quite restrictive, we show in Sec. 6 that Eq. (52) achieves high
accuracy, though as expected it performs somewhat worse than the two-term
approach.

Figure 6: Plots of average posterior log-probabilities (normalized) for the
optimal filter (plot A), two-term filter (plot B), and single-term filter
(plot C) generated from the same parameters and measurement set used in Figure
5. The solid gray line indicates the time at which the filter correctly
identifies that an error has occurred on the third qubit (dashed black line
marks the time of the error). The optimal and two-term filter show almost
identical behavior, while the single-term filter takes longer to detect the
error.

In Figure 6 we plot the average evolution of the posterior log-probabilities
in the presence of an error using outputs from the optimal filter versus those
of the single-term and two-term filters. The plots were generated from the
same data used to create Figure 5, and clearly show that all three filters are
able to detect the error on the third qubit. The single-term filter performs
worst as expected, taking longer to identify the error. Impressively, the two-
term filter behaves almost identically to the optimal filter, with small
discrepancies emerging for $\ell_{i}\in\\{3,5,6,7\\}$ since these states all
have more than two significant LogSumExp terms.

###  5.5 Normalization

As discussed in Sec. 4, one of the practical challenges of using Bayesian
methods is the need for periodic normalization of the probability outputs,
without which the values will grow or shrink rapidly. This behavior is shown
in Figure 7, where the outputs of the linearized Wonham filter from Eq. (31)
grow exponentially with time. Since exponential changes become linear on a
logarithmic scale, the outputs of our logarithmic filters instead scale
linearly with time. In situations where error tracking needs to be done for
only a short duration, the linear growth of the outputs is likely tolerable
and can simply be ignored. This represents a significant improvement over the
linearized Wonham filter, which must incorporate a costly normalization
routine to be viable over virtually any time scale.

Figure 7: This plot tracks the magnitude of the largest probability output as
a function of time for the linear Wonham filter and two different versions of
our two-term filter (the single-term version behaves identically). The solid
lines are averages across $10^{5}$ trials, while the shaded regions enclose
one standard deviation (omitted for the Wonham filter). The blue line records
the magnitude of the two-term filter without any normalization, demonstrating
a clear linear pattern. The orange line shows the same algorithm with the
analytic correction factor from Eq. (53) included, which dramatically reduces
the largest magnitude (mean and standard deviation of about 30 after 100 $\mu
s$).

For longer runs, or in cases where the range of output values is limited due
to hardware restrictions, we may wish to slow the growth of the output
magnitudes even further. This can be achieved by calculating the average rate
of change per time step and simply subtracting this quantity at the end of
each step. The average change $\Delta$ is given by

$\Delta=-[1+\log(2\pi\frac{k}{T})-3\log\cosh(\mu T)+3\mu T],$ (53)

which is the sum of $\log J_{\ell\ell}$ and $G^{\ell}_{\ell}$ when the true
state of the system is $\ell$ (since this is the only significant
contribution). The corrected outputs are represented by the orange line in
Figure 7. Applying the correction greatly slows the growth of the unnormalized
log-probabilities, which allows the filters to be run within a very narrow
numerical range.

##  6 Performance

To evaluate the performance of our single-term and two-term filters in an
error-correction setting, we simulate a large number of trajectories (see
Appendix A.5) and then record how accurately the filters are able to correctly
identify the final state. As our definition of “accuracy” we adopt a binary
measure which is 1 when the filter predicts either the true state or a state
that differs from the true state by a single bit-flip, and is 0 otherwise. At
any given time there will be four states considered to be correct and another
four states considered to be incorrect, which gives an expected accuracy of
50% when guessing randomly. We choose this measure of accuracy because errors
on a single qubit can be corrected via simple majority vote, while errors on
two or more qubits will signify a logical error which the filter was supposed
to have prevented.

Secs. 6.1, 6.2, and 6.3 present and discuss the performance of the filters
under different run durations, error rates $\mu$, and step sizes $T$
respectively. For reference, the performances of our single-term and two-term
filters are plotted alongside those of the optimal filter from Eq. (29) to
assess the practical impact of the simplifying assumptions made in Sec. 5.4.
We also include plots for the linearized Wonham filter from Eq. (31) and the
double-threshold scheme from Atalaya and Zhang et al. to see how well our
filters perform relative to existing algorithms.

###  6.1 Experiment duration

Figure 8: Plot of average inaccuracy (1 - accuracy) of the five different
filters as a function of run time. The error rate was fixed at $2.5\times
10^{-3}\mu\text{s}$ and the time step at $10^{-1}\mu\text{s}$. The accuracy of
the two-term filter is effectively equal to that of the optimal filter, such
that the two curves completely overlap, with the single-term version
performing slightly worse. The linear Wonham filter and double threshold
performed comparably, though both did significantly worse than the logarithmic
filters.

Figure 8 shows the performance of our single-term and two-term filters
alongside the reference algorithms as a function of run duration, with the
longest runs lasting a full millisecond. The accuracy is a decreasing function
of time, since the probability of a logical error will increase as the run
grows longer. Mohseninia et al. showed that the accuracy will decrease
linearly at small error rates, a pattern that is clearly visible in our data
across all filters [5].

As expected, the optimal filter achieves the highest accuracy, though it is
almost exactly matched by the two-term filter. Despite its greater simplicity,
the single-term filter performs within a percentage point of the optimal
model. By contrast, the linear Wonham filter and double threshold algorithm
both perform significantly worse, though for different reasons. The linearized
Wonham filter is constructed explicitly from a first-order approximation,
which is exact as $T\rightarrow 0$ but suboptimal for any finite step size.
The double threshold scheme by contrast is simply suboptimal by design, as it
does not directly incorporate the error probabilities or the Gaussian
distribution of the noise into its filter.

###  6.2 Error rate

Figure 9: Plot of average inaccuracy (1 - accuracy) of the five different
filters as a function of error rate, with both axes on a log scale. The
duration of the runs was fixed at 100 $\mu\text{s}$ and the time step at
$10^{-1}\ \mu\text{s}$. The accuracy of the two-term filter is effectively
equal to that of the optimal filter (curves overlap completely), with the
single-term version performing slightly worse. The linear Wonham filter and
double threshold performed comparably, though both did significantly worse
than the logarithmic filters.

Figure 9 shows the accuracy of the five filter as a function of the per-qubit
bit-flip error rate $\mu$. On a log-log plot the dependence between the
accuracy and error rate is linear for small error rates, but levels off at 0.5
as the error rate approaches $0.1\ (\mu s)^{-1}$. In the latter regime the
system experiences enough logical errors that all memory of the initial state
is lost, and the system is therefore unable to distinguish the true state from
its complement. As before, the optimal and two-term filters perform almost
identically, while single-term term algorithm is slightly less accurate. The
negative effect of non-zero $T$ on the linear Wonham filter is especially
obvious at small error rates, where the accuracies of all the other filters
converge together at roughly double the accuracy of the Wonham filter.

###  6.3 Time step

Figure 10: Plot of average inaccuracy (1 - accuracy) of the five different
filters as a function of time step $T$. The duration of the run was fixed at
100 $\mu\text{s}$ and the error rate at $2.5\times 10^{-3}$. The accuracy of
the two-term filter matches that of the optimal curve at most step sizes, but
performs slightly worse at very small and very large values of $T$. The
single-term filter shows a significant improvement as the time step increases,
while the accuracy of the linear Wonham filter decreases rapidly. The double
threshold appears largely unaffected by the step size.

Lastly, Figure 10 shows the accuracy of the filters as a function of the
measurement integration time $T$, and reveals a large number of different
trends. The time steps range from $1$ ns to 1 $\mu s$, mapping to noise
variances on the order of $100$ to $0.1$ respectively, which allows us to
explore the behavior of these filters under both high-noise and low-noise
conditions.

Focusing first on the double-threshold scheme, it is clear that the step size
does not have a significant effect on its performance. This is expected, as
the thresholding procedure already includes an exponential moving average
which smooths over the Gaussian noise regardless of step size. Indeed, since
smaller values of $T$ result in more measurements per time interval, the
stronger noise is balanced out by the greater sample size of the average. This
is the same trade-off described in Eq. (10) in the context of mean value
estimation.

Looking next at the linear Wonham filter, Figure 10 provides a very clear
demonstration of its dependence on $T$. At a time step of $1$ ns the filter is
effectively optimal, which is expected given that the differential Wonham
filter of Eq. (30) is known to be exact. However, under more realistic
conditions in which hardware latency and other factors generate longer time
steps, the first-order nature of the filter starts to negatively impact its
performance.

For the single-term and two-term filters, several different trends emerge. The
single-term filter has a much poorer performance at small time steps than the
two-term filter, while at larger time steps their performances converge. This
occurs because the transition probabilities are suppressed as $T\rightarrow
0$, which means that almost all of the contributions to the LogSumExp are
small. For the two-term filter this is not a significant problem, since taking
both the largest and second-largest terms allows for the log-probability to
accumulate over time. The outputs of the one-term filter cannot build up in
this manner, so the log-probabilities are reduced far below their true values.
As $T$ increases this effect is diminished, so the performance of the two
filters will begin to converge. Figure 10 also shows that the two-term filter
diverges from the optimal filter at large $T$, since the single-error Gaussian
approximation made in Sec. 5.1 starts to break down when $\mu T$ is large and
$\frac{k}{T}$ is small.

##  7 Discussion

The core objective of our work here was to devise a quantum error detection
filter that is grounded in formal Bayesian analysis, yet also practical to
implement on real-world quantum hardware. We have assumed that the continuous
measurements are represented as a set of discrete signal averages taken over a
time $T$, which then serve as inputs to the filter. This time step $T$ cannot
be avoided when working with actual hardware, since even if weak measurements
could be made on arbitrarily short time scales there would still be latency
associated with transferring and processing the data. As discussed in Sec.
2.2, the fact that $T$ must be non-zero has important effects on the
theoretical distribution of the measurement values, which become more
significant as $T$ and $\mu$ increase.

The filter derived through our Bayesian analysis in Sec. 3 can be understood
as as a discretized version of the Wonham filter for finite $T$, though our
derivation did not originate from its stochastic differential equation. Much
of the difficulty in that analysis stemmed from the fact that errors can occur
within the integration period of the measurement, which causes the underlying
syndrome means to be distributed across the entire $[-1,1]$ interval instead
of being constrained to only $\pm 1$. Even with this additional degree of
freedom, a closed-form solution for the marginal distribution of either
syndrome can be derived by simply summing Eq. (32) over $0\leq e_{k}<\infty$.
The truly challenging part of the analysis comes when the two syndrome
distributions are coupled together by errors on the second qubit, which
induces complex inter-dependencies between $\bar{S}^{(1)}_{i}$ and
$\bar{S}^{(2)}_{i}$. One avenue for future work could lie in analyzing the
form of these joint distributions, with the goal of determining whether a
convenient or illuminating analytic form exists. In the absence of such a
form, some type of simplification will always be necessary in order to utilize
the filter.

Indeed, the single-error assumption that we made when deriving our log-
probability filters in Sec. 5 was designed explicitly to avoid the
mathematical challenges associated with the joint syndrome distributions. This
assumption was the most restrictive that we made when developing the filters,
and it becomes increasingly inaccurate as the average number of errors per
step ($\mu T)$ increases. While it’s true that the error tracking problem as a
whole becomes quite challenging when $\mu T$ is large, the question of how a
filter could best be designed for such a situation warrants further
exploration.

However, in the preferred regime where $\mu T$ is small, the numerical results
in Sec. 6 show that our single-term and two-term log-probability filters are
highly effective. In every test that we conducted these two filters
outperformed both the first-order Wonham filter and double threshold algorithm
by a significant margin, with the two-term filter being virtually optimal over
a wide range of $\mu$ and $T$ values. The single-term filter performs slightly
worse as expected, but is surprisingly effective given its simplicity.
Further, Figure 10 shows that the performance of the single-term filter
increases as $T$ grows larger, up to a value of about $10^{-1}\ \mu\text{s}$.
This, combined with the fact that the optimal filter experiences only a
negligible improvement in performance when moving to smaller $T$, suggests
that effort spent on latency reduction is likely to provide diminishing
returns for filter accuracy.

Given the rapidly growing size of modern quantum hardware, it is safe to
assume that error correction will remain an integral component of quantum
computation for the foreseeable future. The three-qubit toy model analyzed
here is insufficient to protect against arbitrary errors, so an obvious
extension of our work would focus on developing continuous error correction
filters for larger codes e.g. the Shor, Steane, or subsystem codes which
provide more comprehensive protection [23][24]. Our results indicate that
Bayesian filters similar to those presented here would find great success on
these more complex systems.

## Acknowledgements

I.C. was supported by the National Aeronautics and Space Administration under
Grant/Contract/ Agreement No. 80NSSC19K1123 issued through the Aeronautics
Research Mission Directorate. This material is also based upon work supported
by the U.S. Department of Energy, Office of Science, National Quantum
Information Science Research Centers, Quantum Systems Accelerator.

## References

  * [1] Peter W. Shor  “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer” arXiv: quant-ph/9508027  In _SIAM Journal on Computing_ 26.5, 1997, pp. 1484–1509  DOI: 10.1137/S0097539795293172
  * [2] “Quantum Error Correction”  Cambridge University Press, 2013  DOI: 10.1017/CBO9781139034807
  * [3] W.. Wonham  “Some Applications of Stochastic Differential Equations to Optimal Nonlinear Filtering”  In _Journal of the Society for Industrial and Applied Mathematics Series A Control_ 2.3 Society for IndustrialApplied Mathematics, 1964, pp. 347–369  DOI: 10.1137/0302028
  * [4] Juan Atalaya, Mohammad Bahrami, Leonid P. Pryadko and Alexander N. Korotkov  “Bacon-Shor code with continuous measurement of noncommuting operators”  In _Physical Review A_ 95.3 American Physical Society, 2017, pp. 032317  DOI: 10.1103/PhysRevA.95.032317
  * [5] Razieh Mohseninia et al.  “Always-On Quantum Error Tracking with Continuous Parity Measurements”  In _Quantum_ 4 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2020, pp. 358  DOI: 10.22331/q-2020-11-04-358
  * [6] J. Atalaya et al.  “Continuous quantum error correction for evolution under time-dependent Hamiltonians”  In _Physical Review A_ 103.4 American Physical Society, 2021, pp. 042406  DOI: 10.1103/PhysRevA.103.042406
  * [7] Satoshi Morita and Hidetoshi Nishimori  “Mathematical foundation of quantum annealing”  In _Journal of Mathematical Physics_ 49.12 American Institute of Physics, 2008, pp. 125210  DOI: 10.1063/1.2995837
  * [8] Kurt Jacobs and Daniel A. Steck  “A straightforward introduction to continuous quantum measurement”  In _Contemporary Physics_ 47.5, 2006, pp. 279–303  DOI: 10.1080/00107510601101934
  * [9] Alexander N. Korotkov  “Quantum Bayesian approach to circuit QED measurement with moderate bandwidth”  In _Physical Review A_ 94.4, 2016, pp. 042326  DOI: 10.1103/PhysRevA.94.042326
  * [10] Howard M. Wiseman and Gerard J. Milburn  “Quantum Measurement and Control”  Cambridge University Press, 2009  DOI: 10.1017/CBO9780511813948
  * [11] Vic Barnett  “Comparative statistical inference”, Wiley series in probability and mathematical statistics  Wiley, 1999
  * [12] Alexander N. Korotkov  “Continuous quantum measurement of a double dot”  In _Physical Review B_ 60.8 American Physical Society, 1999, pp. 5737–5742  DOI: 10.1103/PhysRevB.60.5737
  * [13] Bruno De Finetti, Antonio Machì, Adrian F.. Smith and Bruno De Finetti  “Theory of probability: a critical introductory treatment”  John Wiley & Sons, 2017
  * [14] D.. Sivia and J. Skilling  “Data analysis: a Bayesian tutorial”, Oxford science publications  Oxford University Press, 2006
  * [15] J.. Norris  “Markov Chains”  Cambridge University Press, 1997  DOI: 10.1017/CBO9780511810633
  * [16] Aiken Pang and Peter Membrey  “Beginning FPGA: Programming Metal”  Apress, 2017  DOI: 10.1007/978-1-4302-6248-0
  * [17] Hideo Mabuchi  “Continuous quantum error correction as classical hybrid control”  In _New Journal of Physics_ 11.10, 2009, pp. 105044  DOI: 10.1088/1367-2630/11/10/105044
  * [18] A.. Borisov  “L1-Optimal Filtering of Markov Jump Processes. I. Exact Solution and Numerical Implementation Schemes”  In _Automation and Remote Control_ 81.11, 2020, pp. 1945–1962  DOI: 10.1134/S0005117920110016
  * [19] Yin Gang George, Qing Zhang and Yuanjin Liu  “Discrete-time approximation of Wonham filters”  In _Journal of Control Theory and Applications_ 2.1, 2004, pp. 1–10  DOI: 10.1007/s11768-004-0017-7
  * [20] G.. Hieftje  “Signal-to-noise enhancement through instrumental techniques. II. Signal averaging, boxcar integration, and correlation techniques”  In _Analytical Chemistry_ 44.7, 1972, pp. 69A–78a  DOI: 10.1021/ac60315a036
  * [21] Pierre Blanchard, Desmond J Higham and Nicholas J Higham  “Accurately computing the log-sum-exp and softmax functions”  In _IMA Journal of Numerical Analysis_ , 2020, pp. draa038  DOI: 10.1093/imanum/draa038
  * [22] M. Haselman et al.  “A Comparison of Floating Point and Logarithmic Number Systems for FPGAs”  In _13th Annual IEEE Symposium on Field-Programmable Custom Computing Machines (FCCM’05)_ IEEE, 2005, pp. 181–190  DOI: 10.1109/FCCM.2005.6
  * [23] Peter W. Shor  “Scheme for reducing decoherence in quantum computer memory”  In _Physical Review A_ 52.4, 1995, pp. R2493–R2496  DOI: 10.1103/PhysRevA.52.R2493
  * [24] A.. Steane  “Error Correcting Codes in Quantum Theory”  In _Physical Review Letters_ 77.5 American Physical Society, 1996, pp. 793–797  DOI: 10.1103/PhysRevLett.77.793

##  Appendix A Appendix

###  A.1 Syndrome distributions

####  A.1.1 General form

In Eq. (22), we expand the syndrome distribution
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$ in terms of the
error probabilities $p(e_{k}|\ell_{i}\ell_{i-1})$ and the conditional
distribution
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})$ of
observing a particular pair of syndrome means given a number of errors $e_{k}$
on the $k$th qubit. Figure 1C shows that the values of $\bar{S}^{(1)}_{i}$ and
$\bar{S}^{(2)}_{i}$ depend explicitly on the number and location of the bit-
flip errors within the integration interval, so we introduce a new set of
variables $\\{x^{k_{j}}_{j}\\}$ which denote the positions of the $N\equiv
e_{1}+e_{2}+e_{3}$ errors. The index $j$ indicates the order that the errors
occur, such that $x^{k_{j-1}}_{j-1}\leq x^{k_{j}}_{j}$, while
$k_{j}\in\\{1,2,3\\}$ denotes on which qubit the error occurs. Note that these
“positions” are in fact times $t$ during the run, but for convenience we
define the beginning of the $i$th integration period to be at $t=0$ and then
divide the times by $T$ so that they all lie on the unit-less interval
$[0,1]$. The value $0\leq x^{k_{j}}_{j}\leq 1$ then denotes a location within
this interval. Using this new set of variables, the syndrome density is given
by the marginalization over all error locations and all possible qubit
assignments:

$\begin{split}P(\bar{S}^{(1)}_{i}&\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})=\\\
&\sum^{k_{1}=3}_{k_{1}=1}...\sum^{k_{N}=3}_{k_{N}=1}\int^{1}_{0}dx^{k_{1}}_{1}\int^{1}_{x^{k_{1}}_{1}}dx^{k_{2}}_{2}...\int^{1}_{x^{k_{N-1}}_{N-1}}dx^{k_{N}}_{N}P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1}\\{x^{k_{j}}_{j}\\})P(\\{x^{k_{j}}_{j}\\}|e_{1}e_{2}e_{3}),\end{split}$
(54)

where we have now further factorized the syndrome density into the product of
two new distributions.

Starting first with $P(\\{x^{k_{j}}_{j}\\}|e_{1}e_{2}e_{3})$, we know that the
location of an error and the qubit that it is assigned to are independent, so
we have

$P(\\{x^{k_{j}}_{j}\\}|e_{1}e_{2}e_{3})=P(\\{x_{j}\\}|e_{1}e_{2}e_{3})P(\\{k_{j}\\}|e_{1}e_{2}e_{3}).$
(55)

Given that each qubit is equally likely to experience an error,
$P(\\{k_{j}\\}|e_{1}e_{2}e_{3})$ must be uniform over all valid assignments of
$\\{k_{j}\\}$, where an assigment is valid if exactly $e_{k}$ errors are
assigned to the $k$th qubit. Since the distribution must be normalized, each
configuration will have a probability equal to one over the total number of
valid configurations:

$P(\\{k_{j}\\}|e_{1}e_{2}e_{3})=\begin{cases}\dfrac{e_{1}!e_{2}!e_{3}!}{N!}&\text{when
valid}\\\ 0&\text{otherwise}.\end{cases}$ (56)

Similarly, the errors themselves will be uniformly distributed within the
interval, so the joint distribution of all error locations $\\{x_{j}\\}$ is
one over the integral across all positions

$P(\\{x_{j}\\}|e_{1}e_{2}e_{3})=\left[\int^{1}_{0}dx^{k_{1}}_{1}\int^{1}_{x^{k_{1}}_{1}}dx^{k_{2}}_{2}...\int^{1}_{x^{k_{N-1}}_{N-1}}dx^{k_{N}}_{N}\right]^{-1}=N!,$
(57)

which depends only on the total number of errors.

Turning now to
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1}\\{x^{k_{j}}_{j}\\})$, the
location of the errors and their qubit assignments will completely determine
the value of the syndromes, so we will have a product of Dirac delta
distributions. Figure 1C shows that the syndrome values are given by the sum
of the gaps between errors, with the sign of the gap value alternating after
every flip. Therefore

$\begin{split}P(\bar{S}^{(m)}_{i}|\ell_{i-1}\\{x^{(m)}_{j}\\})&=\delta\left(\pm_{m}\bar{S}^{(m)}_{i}-x^{(m)}_{1}-\sum^{N^{(m)}}_{j=2}(-1)^{j-1}(x^{(m)}_{j}-x^{(m)}_{j-1})-(-1)^{N^{(m)}}(1-x^{(k)}_{N^{(m)}})\right)\\\
&=\delta\left(\pm_{m}\bar{S}^{(m)}_{i}-(-1)^{N^{(m)}}-2\sum^{N^{(m)}}_{j=1}(-1)^{j-1}x^{(m)}_{j}\right),\end{split}$
(58)

where $\\{x^{(m)}_{j}\\}$ is the set of errors which affect the $m$th parity
operator and $N^{(m)}$ is the number of errors in this set.

Using Eqs. (55–58), we can rewrite Eq. (54) as

$\begin{split}P(\bar{S}^{(1)}_{i}&\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})=\\\
&e_{1}!e_{2}!e_{3}!\sum_{\\{k_{j}\\}\text{ valid
}}\int^{1}_{0}dx^{k_{1}}_{1}\int^{1}_{x^{k_{1}}_{1}}dx^{k_{2}}_{2}...\int^{1}_{x^{k_{N-1}}_{N-1}}dx^{k_{N}}_{N}P(\bar{S}^{(1)}_{i}|\ell_{i-1}\\{x^{(1)}_{j}\\})P(\bar{S}^{(2)}_{i}|\ell_{i-1}\\{x^{(2)}_{j}\\}),\end{split}$
(59)

which can be understood as the sum of volume integrals over the values of
$\\{x^{(k_{j})}_{j}\\}$ which satisfy the Dirac distributions from Eq. (58).
In principle, all that remains is to solve for the appropriate integral
boundaries and then evaluate the integrals, though this is difficult to do for
an arbitrary distribution of errors. In Sec. A.1.2 we describe how Eq. (59)
can be evaluated when errors occur on only a single qubit, while in Sec. A.1.3
we discuss the challenges of generalizing this to errors on multiple qubits,
specifically when the second qubit is involved.

####  A.1.2 Distributions for single-qubit errors

If errors are constrained to occur on only a single qubit, then one of the
syndrome means will behave in a trivial manner. For errors on the first qubit
we will have $\bar{S}^{(2)}_{i}=\pm_{2}1$, and for errors on the the third
qubit we will have $\bar{S}^{(1)}_{i}=\pm_{1}1$. With errors on the second
qubit both syndromes will change in an identical manner, so we can choose to
model explicitly the behavior of $\bar{S}^{(1)}_{i}$ and then just set
$\bar{S}^{(2)}_{i}=\bar{S}^{(1)}_{i}$.

We begin by considering a scenario where the first qubit experiences $N$
errors, while $e_{2}=e_{3}=0$. Under these conditions,
$P(\bar{S}^{(2)}_{i}|\ell_{i-1}\\{x^{(2)}_{j}\\})=\delta(\bar{S}^{(2)}_{i}\mp_{2}1)$,
as $\\{x^{(2)}_{j}\\}$ is empty and thus there are no errors on the second or
third qubits to flip the initial syndrome of the second parity operator.
Indeed, the syndrome distribution of Eq. (59) factorizes completely, with the
marginal distribution of $\bar{S}^{(1)}_{i}$ given by

$P(\bar{S}^{(1)}_{i}|\ell_{i-1};N00)=N!\int^{1}_{0}dx_{1}\int^{1}_{x_{1}}dx_{2}...\int^{1}_{x_{N-1}}dx_{N}\delta\left(\pm_{1}\bar{S}^{(1)}_{i}-(-1)^{N}-2\sum^{N}_{j=1}(-1)^{j-1}x_{j}\right).$
(60)

Note that the sum of Eq. (59) has collapsed into a single term, since there is
only one way to assign $N$ errors to a single qubit. Given that $k_{j}=1$ for
all $j$, we have removed the superscripts from the error locations.

As mentioned previously, we can view Eq. (60) as a volume integral over the
regions of the variable space which give the desired value of
$\bar{S}^{(1)}_{i}$. To derive the integration bounds, we consider first the
constraints placed on $x_{1}$. The first error will contribute
$\pm_{1}(x_{1}-0)=\pm_{1}x_{1}$ to the value of the syndrome, which cannot be
altered by any subsequent errors. Therefore, the range of $\bar{S}^{(1)}_{i}$
after the first error is

$\pm_{1}x_{1}-(1-x_{1})\leq\bar{S}^{(1)}_{i}\leq\pm_{1}x_{1}+(1-x_{1}),$ (61)

since the rest of interval can be entirely negative or entirely positive at
the two extremes. Writing the limits with respect to $x_{1}$ gives

$x_{1}\leq\frac{1\pm_{1}\bar{S}^{(1)}_{i}}{2},$ (62)

which becomes the upper integration bound of the first integral. We can repeat
this same process for each $x_{j}$, taking the sum of the contributions from
errors up to $x_{j}$ and either adding or subtracting the remaining interval
length to get the upper and lower bounds respectively on $\bar{S}^{(1)}_{i}$.
These bounds are then rearranged to derive an upper bound on $x_{j}$, which
has the generic form

$B_{j}\equiv\frac{1}{2}+(-1)^{j-1}\left(\frac{\pm_{1}\bar{S}^{(1)}_{i}}{2}+\sum^{j-1}_{n=1}(-1)^{n}x_{n}\right),$
(63)

such that the integration limits are $x_{j-1}\leq x_{j}\leq B_{j}$. These
bounds have the recursion relation

$B_{j}=(1+x_{j-1})-B_{j-1}$ (64)

which is key to evaluating the integrals.

Using Eq. (63) and Eq. (64), the volume integral in Eq. (60) can be easily
solved. The integral over $x_{N}$ will collapse to $\frac{1}{2}$ due to the
Dirac delta function, and contribute nothing further to the volume.
Considering the next two integrals we have

$\begin{split}\int^{B_{N-2}}_{x_{N-3}}dx_{N-2}\int^{B_{N-1}}_{x_{N-2}}dx_{N-1}&=\int^{B_{N-2}}_{x_{N-3}}dx_{N-2}(B_{N_{1}}-x_{N-2})\\\
&=(1-B_{N-2})\int^{B_{N-2}}_{x_{N-3}}dx_{N-2}\\\
&=(1-B_{N-2})(1-B_{N-3}),\end{split}$ (65)

where we used the recursive expression $B_{j}-x_{j-1}=1-B_{j-1}$. The result
of the first integral does not have any dependence on $x_{N-2}$ and therefore
is unaffected by the second integral. Indeed, the generic action of the
integral for $x_{j-1}$ on $(1-B_{j})^{a}$ is

$\int^{B_{j-1}}_{x_{j-2}}dx_{j-1}(1-B_{j})^{a}=\frac{1}{a+1}(1-B_{j-2})^{a+1},$
(66)

which maps $B_{j}\rightarrow B_{j-2}$. By repeated application of Eq. (66),
the probability density in Eq. (60) reduces to

$P(\bar{S}^{(1)}_{i}|\ell_{i-1};N00)\propto(1-B_{0})^{N//2}(1-B_{1})^{(N-1)//2},$
(67)

where $x//y$ is $x$ floor-divided by $y$, and where $B_{0}$ is defined from
the recursion relation in Eq. (64) with $j=1$ and $x_{0}\equiv 0$.
Substituting the expressions for $B_{1}$ and $B_{0}$ into Eq. (67) gives

$P(\bar{S}^{(1)}_{i}|\ell_{i-1};N00)=\frac{N!}{2(N//2)!((N-1)//2)!}(\frac{1\pm_{1}\bar{S}^{(1)}_{i}}{2})^{N//2}(\frac{1\mp_{1}\bar{S}^{(1)}_{i}}{2})^{(N-1)//2}$
(68)

after normalization, which is simply a Beta distribution with respect to
$\frac{1\pm\bar{S}^{(1)}_{i}}{2}$ on the interval $[-1,1]$. This result can be
easily generalized to $N$ errors on the third qubit by replacing
$\bar{S}^{(1)}_{i}$ with $\bar{S}^{(2)}_{i}$ and $\pm_{1}$ with $\pm_{2}$.

####  A.1.3 Challenges of multi-qubit errors

When errors occur on the second qubit and either the first or third qubits (or
both), the evaluation of Eq. (54) is far more challenging. First, the sum over
$\\{k_{j}\\}$ will not collapse to a single term, so we must evaluate multiple
volume integrals. More importantly, the error locations assigned to the second
qubit in each of those volume integrals will have to simultaneously satisfy
constraints imposed by both $\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$, which
means that the bounds $B_{j}$ on those errors will be given by

$B_{j}=\min\left[\frac{1}{2}+(-1)^{j-1}\left(\frac{\pm_{1}\bar{S}^{(1)}_{i}}{2}+\sum^{N^{(1)}_{j}}_{n=1}(-1)^{n}x^{(1)}_{n}\right),\frac{1}{2}+(-1)^{j-1}\left(\frac{\pm_{1}\bar{S}^{(2)}_{i}}{2}+\sum^{N^{(2)}_{j}}_{n=1}(-1)^{n}x^{(2)}_{n}\right)\right],$
(69)

where $N^{(m)}_{j}$ is the number of errors occurring earlier than
$x^{k_{j}}_{j}$ that affect the $m$th parity operator. These bounds create
dependencies between the variables in $\\{x^{(1)}_{j}\\}$ and
$\\{x^{(2)}_{j}\\}$ which prevent us from using the simple recursion strategy
employed in Eq. (66). That said, an expression for
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i-1};e_{1}e_{2}e_{3})$ can still
be calculated from Eq. (54) in a brute-force manner using symbolic computation
tools, though this becomes increasingly impractical as the values of $e_{1}$,
$e_{2}$, and $e_{3}$ grow.

###  A.2 Implementing the optimal filter

From Eq. (29), the optimal Bayesian filter is a function of the measurement
density $P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ and transition elements
$J_{\ell_{i-1}\ell_{i}}$. The matrix $J$ is easy to compute, but in order to
calculate $P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})$ we must evaluate the
integral in Eq. (21). Although the analytic solution is not available to us,
we can approximate the value of the integral using a Riemann sum

$\begin{split}&P(M^{(1)}_{i}M^{(2)}_{i}|\ell_{i}\ell_{i-1})\approx\\\
&\frac{1}{n^{2}}\sum^{n-1}_{m^{(1)}=-n}\sum^{n-1}_{m^{(2)}=-n}\prod^{2}_{j=1}\frac{\exp[\frac{-T}{2k}(M^{(j)}_{i}-(\frac{m^{(j)}+\frac{1}{2}}{n}))^{2}]}{\sqrt{2\pi\frac{k}{T}}}P(\bar{S}^{(1)}_{i}=\frac{m^{(1)}+\frac{1}{2}}{n},\bar{S}^{(2)}_{i}=\frac{m^{(2)}+\frac{1}{2}}{n}|\ell_{i}\ell_{i-1}),\end{split}$
(70)

where the interval $[-1,1]$ has been discretized into $2n$ evenly-spaced
segments.

The syndrome density is estimated by using Eq. (22) to factorize it into
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|e_{1}e_{2}e_{3})$ and
$P(e_{1}e_{2}e_{3}|\ell_{i}\ell_{i-1})$, the latter of which can be evaluated
analytically. To approximate the value of
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|e_{1}e_{2}e_{3})$, we construct
histograms of $\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$ with $4n^{2}$ bins
for values of $e_{1}$, $e_{2}$, and $e_{3}$ such that $e_{1}+e_{2}+e_{3}\leq
N$. The value of $N$ is chosen so that the probability of experiencing more
than $N$ errors in a given interval is below some threshold value. Since $\mu$
and $T$ only impact $P(e_{1}e_{2}e_{e}|\ell_{i}\ell_{i-1})$, these histograms
can be reused for multiple error rates and integration times. To reconstruct
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$, we perform the sum
in Eq. (22) using the histograms and the analytic error probabilities to
generate a single combined histogram.

The sum in Eq. (70) must be evaluated at every time step, which becomes a
significant computational bottleneck for large $n$. To speed up the filter, we
constructed a lookup table for the Gaussian distribution at discrete values
across $M^{(1)}_{i}$ and $M^{(2)}_{i}$ and then used bi-linear interpolation
to compute values for continuous inputs. Given the smoothness of the Gaussian
function, these interpolations were highly accurate and very fast to compute.

###  A.3 Convergence to the Wonham filter

The (unnormalized) Bayesian filter given in Eq. (29) can be shown to converge
to the linear Wonham filter of Mohseninia et al. as $T\rightarrow 0$, and thus
converge to the Wonham filter after normalization. Up to the smallest order in
$T$, the transition matrix elements are given by

$J_{\ell_{i-1}\ell_{i}}\approx(\mu T)^{d(\ell_{i-1},\ell_{i})}(1-3\mu T),$
(71)

where $d(\ell_{i-1},\ell_{i})$ is the Hamming distance between the three-bit
representations of $\ell_{i}$ and $\ell_{i-1}$. Ignoring the normalization
factor, the measurement density can be expanded up to its smallest order in
$T$ as

$\begin{split}P(M^{(1)}_{i}&M^{(2)}_{i}|\ell_{i}\ell_{i-1})\propto\\\
&\int^{1}_{-1}\int^{1}_{-1}\left(1-\frac{T}{2k}\left[(M^{(1)}_{i}-\bar{S}^{(1)}_{i})^{2}+(M^{(2)}_{i}-\bar{S}^{(2)}_{i})^{2}\right]\right)\tilde{P}(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})d\bar{S}^{(1)}_{i}d\bar{S}^{(2)}_{i},\end{split}$
(72)

where $\tilde{P}(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})$ is
equal to $P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i})|e_{1}e_{2}e_{3})$ such that
$e_{k}=1$ if there is a net flip on the $k$th qubit when moving from
$\ell_{i-1}$ to $\ell_{i}$ and $e_{k}=0$ otherwise. This is the syndrome
density with the fewest errors that is still consistent with the state
transition.

Substituting Eq. (71) and Eq. (72) into Eq. (29) gives

$\hat{P}(\ell_{i})\propto\sum_{\ell_{i-1}}(\mu
T)^{d(\ell_{i-1},\ell_{i})}(1-3\mu
T)(1-\frac{T}{2k}I_{\ell_{i-1}\ell_{i}})P(\ell_{i-1}),$ (73)

where $I$ is a matrix with elements given by

$I_{\ell_{i-1}\ell_{i}}=\int^{1}_{-1}\int^{1}_{-1}\left[(M^{(1)}_{i}-\bar{S}^{(1)}_{i})^{2}+(M^{(2)}_{i}-\bar{S}^{(2)}_{i})^{2}\right]\tilde{P}(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|\ell_{i}\ell_{i-1})d\bar{S}^{(1)}_{i}d\bar{S}^{(2)}_{i}.$
(74)

The sum in Eq. (73) can be simplified by keeping only those terms which are at
most first-order in $T$. With $\mathcal{I}\equiv\\{1,2,4\\}$ containing the
set of indices corresponding to bit-flips on the first, second, and third
qubits respectively, we have

$\hat{P}(\ell_{i})\propto
P(\ell_{i-1}=\ell_{i})+T\left[\mu\sum_{i\in\mathcal{I}}P(\ell_{i-1}=\ell_{i}\oplus
i)-(3\mu+\frac{1}{2k}I_{\ell_{i}\ell_{i}})P(\ell_{i-1}=\ell_{i})\right],$ (75)

which is equivalent to Eq. (31) after introducing $Q$ and evaluating
$I_{\ell_{i}\ell_{i}}$.

###  A.4 Gaussian fit of single-error distribution

In Eq. (35) we considered the distribution of $M^{(1)}_{i}$ when a single
error occurs on the first qubit, deriving the expression

$P(M^{(1)}_{i}|\ell_{i-1};100)=\frac{1}{4}\left[\text{erf}(\frac{M^{(1)}_{i}+1}{\frac{2k}{T}})-\text{erf}(\frac{M^{(1)}_{i}-1}{\frac{2k}{T}})\right].$
(76)

From a practical perspective it would be convenient if the distribution were
Gaussian, so we compute the mean value $s$ and variance $\sigma^{2}$ of Eq.
(76) and then construct a Gaussian distribution off of these parameters. This
can be done by going back to the integral which generated Eq. (76) and simply
switching the order of integration

$\displaystyle s$
$\displaystyle=\frac{1}{2}\int^{1}_{-1}d\bar{S}^{(1)}_{i}\int^{\infty}_{-\infty}dM^{(1)}_{i}M^{(1)}_{i}\frac{\exp[\frac{-T}{2k}(M^{(1)}_{i}-\bar{S}^{(1)}_{i})^{2}]}{\sqrt{2\pi\frac{k}{T}}}=\frac{1}{2}\int^{1}_{-1}\bar{S}^{(1)}_{i}d\bar{S}^{(1)}_{i}=0$
(77) $\displaystyle\sigma^{2}$
$\displaystyle=\frac{1}{2}\int^{1}_{-1}d\bar{S}^{(1)}_{i}\int^{\infty}_{-\infty}dM^{(1)}_{i}(M^{(1)}_{i})^{2}\frac{\exp[\frac{-T}{2k}(M^{(1)}_{i}-\bar{S}^{(1)}_{i})^{2}]}{\sqrt{2\pi\frac{k}{T}}}=\frac{1}{2}\int^{1}_{-1}[(\bar{S}^{(1)}_{i})^{2}+\frac{k}{T}]d\bar{S}^{(1)}_{i}=\frac{1}{3}+\frac{k}{T}.$
(78)

###  A.5 Simulating trajectories

To generate synthetic data of duration $nT$, we must correctly evolve the
state of the system across each of the $n$ integration intervals (of length
$T$), and then generate an appropriate measurement record. To perform this
evolution we utilize the so-called “jump, no-jump” approach, where individual
bit-flips are sampled at each time step and then applied to the system such
that the state remains pure across the run. For particular values of $\mu$ and
$T$ the probability of experiencing $e_{k}$ bit-flip errors on the $k$th qubit
is given by the Poisson distribution in Eq. (16), which can be easily sampled
from. After $e_{1}$, $e_{2}$, and $e_{3}$ have been determined for each of the
$n$ time steps, we apply the errors to our initial state $\ket{000}$ in the
proper order to evolve the system.

To generate the corresponding measurements, we first sample the syndrome means
$\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$ from the measurement distribution
$P(\bar{S}^{(1)}_{i}\bar{S}^{(2)}_{i}|e_{1}e_{2}e_{3})$ conditioned on the
number of errors that occurred in the step. This can be done in
straightforward manner by sampling $e_{k}$ uniform values on the interval
$[0,T]$ for each qubit and then sorting them by value in ascending order. To
calculate $\bar{S}^{(1)}_{i}$, the gaps between errors on the first and second
qubit are added and subtracted in an alternating order, with the first gap
given by the value of the first error and the last gap given by one minus the
value of the last error (see Eq. (58) for more detail). After the gaps have
been summed the syndrome value is multiplied by $\pm_{1}$ to take into account
the parity of the state at the start of interval. An analogous procedure can
be be carried out for $\bar{S}^{(2)}_{i}$ using $\pm_{2}$ and the gaps between
errors on the second and third qubits. Once the syndrome values have been
sampled, the values of $M^{(1)}_{i}$ and $M^{(2)}_{i}$ are generated by simply
adding Gaussian noise with mean zero and variance $\frac{k}{T}$ to
$\bar{S}^{(1)}_{i}$ and $\bar{S}^{(2)}_{i}$.