Machine learning for continuous quantum error correction on superconducting qubits

The superconducting qubits are monitored using voltage signals from homodyne measurements of the parity operators that are derived from tones reflected off the resonator (see appendix A). The resonator signal is fed into a Josephson parametric amplifier (JPA) in order to increase the signal strength without adding a significant amount of noise. The amplified radio frequency signals are then demodulated and digitized. After a further digital demodulation, the signals are processed with an exponential anti-aliasing filter with a time constant of 32 ns. This filtered signal, which is averaged in Δt = 32 ns bins, is then streamed from the digitizer card to the computer.

Due to the effects of the amplifier and resonator, we expect that measurements performed on such real superconducting devices will deviate from the idealized behavior predicted by equation (1). In particular, we can anticipate the following three imperfections:

• (a)  
  The noise will possess a high degree of positive auto-correlation at short temporal lags due to the narrow low-pass bandwidth of the JPA and anti-aliasing filter.
• (b)  
  When a bit-flip occurs, the syndrome means will change gradually rather than instantaneously as the resonator reaches its new steady state. These periods are referred to as resonator transients to stress their temporary nature, and arise because of time-dependent changes in the measurement strength ${{\Gamma}}_{m}^{k}$ (see appendix D).
• (c)  
  The values of the syndromes will drift over time due to small changes in experimental conditions (e.g. temperature). Unlike the other imperfections, this effect is only noticeable when comparing across quantum trajectories rather than within them.

These non-ideal behaviors in the measurement signals extracted from our typical physical experiments will be incorporated into our simulated experiments in section 5.

Figure 1 shows experimental dispersive readouts taken from three transmon qubits [26] over the span of 6 μs [22]. The blue and orange lines are a record of the outputs from the two resonators, each measuring a different pair of qubits for their syndromes. The top figure shows the measurement signals from a single experiment, which contain large amounts of auto-correlated noise. During the experiment an X₂ error was injected at 3.0 μs, flipping the system from |100⟩ to |110⟩, but the weak-measurement noise largely obscures its effect on the syndrome values.

Zoom In Zoom Out Reset image size

To reveal these underlying syndromes, the bottom figure of figure 1 shows an average over the measurements from roughly 47 500 experiments, each initialized to |100⟩ and injected with an X₂ error at 3.0 μs. It takes approximately 2 μs after initialization for the syndromes to reach their steady-state values for |100⟩, as the number of photons in each resonator increases from zero gradually. We ignore this effect in our analysis, as it will only occur once at the start of an experiment. After the X₂ error is injected, the syndromes do not instantaneously jump to a new pair of values but instead enter a transitory period which can include significant oscillations. These transients derive from the time-dependent changes in the measurement rate ${{\Gamma}}_{m}^{k}(t)$ analyzed in appendix D. This period lasts for roughly 2 μs, after which the syndromes stabilize at their new steady-state values for |110⟩.

Depending on the underlying hardware, a measurement signal may be generated on a wide variety of different scales, such as the arbitrary voltage scale in figure 1. To denote a signal generically on any scale, we write the measurement samples as

$$\begin{equation}{I}_{k,t}={\bar{S}}_{k,t}+\sqrt{{\tau }_{k}}{\varepsilon }_{t},\end{equation} \tag{ 4 }$$

where ${\bar{S}}_{k,t}$ is the scaled mean of the kth resonator at step t, τ_k is the scaled variance of the kth resonator, and ${\varepsilon }_{t}\sim \mathcal{N}(0,1)$. In this notation, the physical quantities Γ_m and Δt from equation (2) have been absorbed into ${\bar{S}}_{k,t}$ and τ_k .

Unlike the other imperfections, the challenge posed by auto-correlated signal noise can be characterized theoretically. If the Gaussian noise in I_k,t is correlated, then the distribution of noise samples can be parameterized in terms of a covariance matrix Σ whose off-diagonal elements determine the degree of correlation. For simplicity we restrict our analysis to dependencies that are Markovian, such that I_k,t depends only on the preceding measurement I_k,t−1, though our conclusions are not limited to this regime. Using a correlation coefficient of 0 < ρ < 1, the joint Gaussian log-density describing I_k,t and I_k,t−1 is

$$\begin{equation*}\begin{aligned}\hfill \mathrm{log}\,p({I}_{k,t}{I}_{k,t-1}\vert {\bar{S}}_{k,t})=-\frac{1}{2{\tau }_{k}(1-{\rho }^{2})}\left[\begin{matrix}\hfill {\tilde{I}}_{k,t}\hfill & \hfill {\tilde{I}}_{k,t-1}\hfill \end{matrix}\right]\left[\begin{matrix}\hfill 1\hfill & \hfill -\rho \hfill \\ \hfill -\rho \hfill & \hfill 1\hfill \end{matrix}\right]\left[\begin{matrix}\hfill {\tilde{I}}_{k,t}\hfill \\ \hfill {\tilde{I}}_{k,t-1}\hfill \end{matrix}\right]+A,\end{aligned}\end{equation*}$$

where ${\tilde{I}}_{k,j}\equiv {I}_{k,j}-{\bar{S}}_{k,j}$ denotes the centered signal sample at step j and A is the log of the normalization constant. We shall assume hereafter that the signal has been rescaled such that ${\bar{S}}_{k,j}=\pm 1$.

The effect of auto-correlations on error correction is best characterized in terms of how it impacts the usefulness of the syndrome measurements. To be more precise, we know that the purpose of each measurement is to provide some information about whether the underlying syndrome value of the state is 1 or −1. When framed in these terms, we can formalize and quantify a notion of measurement 'usefulness' using Bayesian theory, specifically a ratio called the Bayes factor which we denote as ϕ [27]. This factor can be written in log form as

$$\begin{equation}\mathrm{l}\mathrm{o}\mathrm{g}\;{\phi }_{k,t}=\mathrm{l}\mathrm{o}\mathrm{g}\;\;p({I}_{k,t}\vert {I}_{k,t-1},{\bar{S}}_{k,t}=1)-\mathrm{l}\mathrm{o}\mathrm{g}\;p({I}_{k,t}\vert {I}_{k,t-1},{\bar{S}}_{k,t}=-1),\end{equation} \tag{ 5 }$$

and quantifies how much evidence I_k,t gives about the underlying syndrome value if we have already seen the previous measurement I_k,t−1. The larger the magnitude of log ϕ_k,t the more useful I_k,t is for our task, with its sign simply indicating whether the evidence supports a value of 1 or −1.

Let Q = Σ⁻¹. By making the substitutions σ⁻¹ = Q₂₂ and $\mu ={\bar{S}}_{k,t}-{Q}_{12}/{Q}_{22}({I}_{k,t-1}-{\bar{S}}_{k,t})$ in the unconditional log-densities $-{({I}_{k,t}-\mu )}^{2}/(2\sigma )+A$, each of the conditional log-densities in equation (5) can be written as

$$\begin{equation*}\begin{aligned}\hfill \mathrm{log}\,p({I}_{k,t}\vert {I}_{k,t-1},{\bar{S}}_{k,t})=-\frac{{[{I}_{k,t}-{\bar{S}}_{k,t}-\rho ({I}_{k,t-1}-{\bar{S}}_{k,t})]}^{2}}{2{\tau }_{k}(1-{\rho }^{2})}+A,\end{aligned}\end{equation*}$$

where A is again the normalization constant [28]. Expanding the numerator and keeping only the terms that depend on ${\bar{S}}_{k,t}$ gives

$$\begin{equation*}\mathrm{log}\,p({I}_{k,t}\vert {I}_{k,t-1},{\bar{S}}_{k,t})\to \frac{{S}_{k}^{2}(\rho -1)+2{\bar{S}}_{k,t}({I}_{k,t}-\rho {I}_{k,t-1})}{2{\tau }_{k}(1+\rho )},\end{equation*}$$

where we ignore the other terms since they will cancel when computing log ϕ_k,t . After substituting this representation back into equation (5) we get

$$\begin{equation}\mathrm{log}\,{\phi }_{k,t}=\frac{2({I}_{k,t}-\rho {I}_{k,t-1})}{{\tau }_{k}(1+\rho )},\end{equation} \tag{ 6 }$$

where the value of log ϕ_k,t depends not only on I_k,t and I_k,t−1 but also on the variance and auto-correlation of the measurements.

To see the impact of the auto-covariance more clearly, we compute the expectation value $\mathbb{E}[\mathrm{log}\,{\phi }_{k,t}]$ with respect to a Gaussian distribution centered on the true syndrome value ${S}_{k,t}^{\prime }=\pm 1$. Since equation (6) is linear, we can simply substitute in ${S}_{k,t}^{\prime }$ for I_k,t and I_k,t−1 to get $\mathbb{E}[\mathrm{log}\,{\phi }_{k,t}]$. After taking its magnitude, we have

$$\begin{equation}\vert \mathbb{E}[\mathrm{log}\,{\phi }_{k,t}]\vert =\frac{2(1-\rho )}{{\tau }_{k}(1+\rho )},\end{equation} \tag{ 7 }$$

which decreases as the value of ρ increases. Equation (7) shows that positive auto-correlation (ρ > 0) in the signal makes each of our measurements less useful than if the noise had been uncorrelated (ρ = 0), which means that it will take longer for us to determine the value of ${\bar{S}}_{k,t}$ at a given measurement strength.

This result can be understood by imagining that ${\bar{S}}_{k,t}$ and I_k,t−1 are competing to determine the value of I_k,t , with smaller ρ favoring ${\bar{S}}_{k,t}$. The more that ${\bar{S}}_{k,t}$ affects the measurement, the more that the measurement in turn tells us about ${\bar{S}}_{k,t}$ and thus the more useful it is to us. When ρ is large, the value of I_k,t tends to lie very close to the value of I_k,t−1 regardless of whether ${\bar{S}}_{k,t}$ is 1 or −1, and therefore the measurement does not reveal much new information about the syndrome.

The double threshold protocol from [21] uses two standard signal processing methods, filtering and thresholding, to identify errors. The raw measurement signal is first passed through an exponential filter to smooth out oscillations, and then this averaged value is compared to a pair of adjustable threshold values to determine the state of the system. A slightly different double threshold protocol was proposed in [20], which used boxcar averaging and fixed one of the thresholds at zero.

To estimate the definite error syndromes from the noisy measurements, we first filter the raw signals I_k (t) to obtain corresponding filtered signals ${\mathcal{I}}_{k}(t)$ according to

$$\begin{equation*}{\dot{\mathcal{I}}}_{k}(t)=-\frac{{\mathcal{I}}_{k}(t)}{\tau }+\frac{{I}_{k}(t)}{\tau },\end{equation*}$$

where τ is the averaging time parameter, and whose discretized version is similar. In the regime where t − t₀ ≫ τ where t₀ is at the last filtered signal reset, ${\mathcal{I}}_{k}(t)$ reads as

$$\begin{equation*}{\mathcal{I}}_{k}(t)={\int }_{{t}_{0}}^{t}\mathrm{d}{t}^{\prime }\frac{{\mathrm{e}}^{-\frac{t-{t}^{\prime }}{\tau }}}{\tau }{I}_{k}\left({t}^{\prime }\right).\end{equation*}$$

Thresholds for error correction

After filtering the measurement signals, we then apply a double thresholding protocol to the filtered signals ${\mathcal{I}}_{1}(t)$ and ${\mathcal{I}}_{2}(t)$ that is parameterized by the two thresholds Θ₁ and Θ₂, where Θ₁ is the threshold for the −1 value of the error syndromes and Θ₂ is the threshold for the +1 value of the error syndromes. If at least one of ${\mathcal{I}}_{1}(t)$ or ${\mathcal{I}}_{2}(t)$ is found to lie within the interval (Θ₁, Θ₂), we declare to be uncertain of the error syndromes and do not perform any error correction operation. Otherwise, we apply the following procedure, in accordance with the standard approach for error diagnosis and correction. If both ${\mathcal{I}}_{1}(t) > {{\Theta}}_{2}$ and ${\mathcal{I}}_{2}(t) > {{\Theta}}_{2}$, then we diagnose the error syndromes as (S₁ = +1, S₂ = +1) and accordingly perform no error correction operation. If ${\mathcal{I}}_{1}(t)< {{\Theta}}_{1}$ and ${\mathcal{I}}_{2}(t) > {{\Theta}}_{2}$, then we diagnose the error syndromes as (S₁ = −1, S₂ = +1) and accordingly perform the error correction operation C_op = X₁. If both ${\mathcal{I}}_{1}(t)< {{\Theta}}_{1}$ and ${\mathcal{I}}_{2}(t)< {{\Theta}}_{1}$, then we diagnose the error syndromes as (S₁ = −1, S₂ = −1) and accordingly perform the error correction operation C_op = X₂. If ${\mathcal{I}}_{1}(t) > {{\Theta}}_{2}$ and ${\mathcal{I}}_{2}(t)< {{\Theta}}_{1}$, then we diagnose the error syndromes as (S₁ = +1, S₂ = −1) and accordingly perform the error correction operation C_op = X₃.

In quantum annealing, we note that the error correction operations are applied immediately after the error syndromes are diagnosed to minimize the aforementioned spurious Hamiltonian evolution. The action of an error correction operation C_op, assumed to be instantaneous, changes the three-qubit state ρ(t) according to

$$\begin{equation*}\rho (t)\to {C}_{\text{op}}\rho (t){C}_{\text{op}},\end{equation*}$$

which applies to other models in our work as well. We note that the parameters {τ, Θ₁, Θ₂} constitutes the minimal set of tunable parameters. When the measurement signals I_k have white noise, their optimal values in minimizing the logical error rate can be obtained by equation (43) in [21] together with numerical optimizations.

We further reset the filtered signals ${\mathcal{I}}_{k}(t)$ to the corresponding initial syndrome value, at the same instant to avoid the transient delay in the filtered signals to reflect the application of the error correction operation on the state. Inherent within any error correction protocol, however, is the implicit assumption that the correction properly removes the error, which may not necessarily be the case if the error was misdiagnosed.

We note that the ${\mathcal{I}}_{k}(t)$ used by the double threshold model in CQEC consists of weighted contributions from every raw signal taken prior to t and after the last correction. The discrete Bayesian model and the RNN-based model that we discuss in this work can both be operated on raw signals, using all historical signals taken prior to a given t. This is in contrast to the projective measurement on ancilla superconducting qubits in discrete QEC that applies a matched filter [29] on raw signals taken only within each detection round.

One weakness of the double-threshold scheme is that its predictions are essentially all-or-nothing, since there is no in-built quantity that expresses the model's confidence. This contrasts with probabilistic classifiers, which generate probability values for each prediction class instead of only a single guess. By framing the classification problem in terms of probabilities, we can incorporate our knowledge of the error and noise distributions into our model in a mathematically rigorous manner.

Since each qubit in our system will experience either one or zero net flips after every time step, there are eight different ways that a state can be altered by bit-flips and therefore eight different classes that our classifier must track. We denote each of the possible bit-flip configuration using the state that |000⟩ is taken to by the error, such that |001⟩ denotes a flip on the third qubit, |110⟩ denotes a flip on the first and second qubits, and so on. The goal of a probabilistic error corrector is to accurately determine the probability of all eight 'error states' at time step t given the measurement histories ${\mathcal{M}}_{t}^{k}\equiv {\left\{{I}_{k,{t}^{\prime }}\right\}}_{{t}^{\prime }=0}^{{t}^{\prime }=t}$. We write this posterior probability as

$$\begin{equation}\hat{p}({s}_{t})\equiv p({s}_{t}\vert {\mathcal{M}}_{t}^{1}{\mathcal{M}}_{t}^{2}),\end{equation} \tag{ 8 }$$

where s_t ∈ {0, ..., 7} denotes the digital representation of the error state at step t.

In the remainder of this subsection we consider a probabilistic classifier constructed using Bayes' theorem, which makes prediction based on the posterior probabilities of the different basis states at each time step [30]. Starting with the knowledge of the initial state, this model uses a Markov chain and a set of Gaussian likelihoods to update our beliefs about the system conditioned on the specific measurement values that we observe.

The Bayesian algorithm described in this section is derived by assuming that the mean of a given measurement I_k,t is always determined by the state of the system at the end of the time step. This is equivalent to assuming that errors always happen at the beginning of each time step (see section 2). Since our method for generating quantum trajectories follows this assumption, the Bayesian model is theoretically optimal for the numerical tests carried out in section 5 without mean drift or resonator transients. As the length of the step Δt between measurements goes to zero, this algorithm converges to the Wonham filter [31], which is known to be optimal for continuous quantum filtering of error syndromes [32]. This filter is similar to the discretized, linear Wonham filter derived in [20], except that our filter does not rely on first-order approximations of the Markov evolution or Gaussian functions.

Model structure

Using Bayes' theorem, the posterior probability of equation (8) can be rearranged into the recursive form

$$\begin{equation}\hat{p}({s}_{t})\propto p({I}_{1,t}{I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1}{\mathcal{M}}_{t-1}^{2})\sum\limits _{i=0}^{7}p({s}_{t}\vert {s}_{t-1}=i)\hat{p}({s}_{t-1}=i),\end{equation} \tag{ 9 }$$

where we assume that the occurrence of an error is independent of any previous measurements and that I_k,t depends on the error state at time t along with past signal values due to auto-correlations.

This recursive expression describes a Bayesian filter which takes prior information about the error state of the system and updates it based on the transition probabilities p(s_t |s_t−1) and measurement likelihoods $p({I}_{1,t}{I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1}{\mathcal{M}}_{t-1}^{2})$. The filter can be easily implemented once we have functional forms for these two terms, which we describe next.

Markovian state transitions

The Markovian assumption inherent in p(s_t |s_t−1) is reasonable, given that the net effect of an additional bit-flip error depends only on the error state the system before the error. We assume hereafter that the error rate γ_q is identical for all three qubits, i.e., γ_q = γ. This allows us to model the errors as a Markov chain [33] with an 8 × 8 rate matrix Q given by

$$\begin{equation}{Q}_{ij}=\begin{cases}-3\gamma \quad \text{if}\enspace j=i\hfill \\ \ \ \gamma \quad \text{if}\enspace j\oplus i\in \left\{1,2,4\right\}\hfill \\ \ \ 0\quad \text{otherwise},\hfill \end{cases}\end{equation} \tag{ 10 }$$

where we define our basis such that index i ∈ {0, ..., 7} corresponds to the error state whose classical binary representation is equal to i, e.g. 5 → |101⟩.

Since Q only gives the rate of transition per unit time, we need to compute the transition matrix J in order to get probabilities for a finite step. This matrix can be derived from Q as

$$\begin{equation*}J={\mathrm{e}}^{Q{\Delta}t},\end{equation*}$$

where Δt is the length of the time step. Element J_ij gives the probability of transitioning from error state i to error state j across the time step, so we can relate p(s_t |s_t−1) to J as p(s_t = j|s_t−1 = i) = J_ij . Using J, the sum in equation (9) can be evaluated to give probabilities $\tilde{p}({s}_{t})$

$$\begin{equation}\tilde{p}({s}_{t}=j)\equiv \sum\limits _{i=0}^{7}\hat{p}({s}_{t-1}=i){J}_{ij},\end{equation} \tag{ 11 }$$

which take into account the transitions induced by bit-flip errors during the time step.

Measurement likelihoods

The measurement likelihood $p({I}_{1,t}{I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1}{\mathcal{M}}_{t-1}^{2})$ describes the probability of generating signal values I_1,t and I_2,t given that the system is in error state s_t and that we had previously measured the values in ${\mathcal{M}}_{t-1}^{1}$ and ${\mathcal{M}}_{t-1}^{2}$. Since the noise from each syndrome is independent, we can factor the likelihood as

$$\begin{equation*}p({I}_{1,t}{I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1}{\mathcal{M}}_{t-1}^{2})=p({I}_{1,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1})p({I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{2})\end{equation*}$$

with I_1,t and I_2,t contributing independently to the probability.

If the noise source is assumed to be Gaussian, then the probability density for each I_k,t has the form

$$\begin{equation*}p({I}_{k,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{k})=\frac{1}{2\pi {\sigma }^{2}}\,\mathrm{exp}\left[\frac{-{({I}_{k,t}-{\mu }_{k,t})}^{2}}{2{\sigma }^{2}}\right],\end{equation*}$$

where μ_k,t and σ² are the mean and variance of the signal conditioned on the past measurements ${\mathcal{M}}_{t-1}^{k}$. In practice the auto-correlations rapidly decay, so we only need to condition on a small number of recent measurements. Hence, we let m_k,t−1 be the vector of these measurements, and let c be the vector of their corresponding covariance values. Then

$$\begin{equation}{\mu }_{k,t}={\bar{S}}_{k,t}+{c}^{T}{{\Sigma}}^{-1}({m}_{k,t-1}-{\bar{S}}_{k,t}\vec{1}),\end{equation} \tag{ 12 }$$

$$\begin{equation}{\,\sigma }^{2}=\frac{\tau }{{\Delta}t}-{c}^{T}{{\Sigma}}^{-1}c,\end{equation} \tag{ 13 }$$

where $\vec{1}$ is a vector of ones with the same dimension as m_k,t−1, Σ is the covariance matrix of the variables in m_k,t−1, and ${\bar{S}}_{k,t}$ is the mean corresponding to error state s_t [28]. Since the system always begins in the coding subspace, each error state maps to a definite error subspace and therefore has definite syndrome values regardless of how the logical state was initialized.

After the measurement pair I_k,t is received, the Gaussian likelihood functions are used to convert the probabilities from equation (11) into the next posteriors $\hat{p}({s}_{t})$ as

$$\begin{equation}\hat{p}({s}_{t})\propto \tilde{p}({s}_{t})\cdot p({I}_{1,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{1})p({I}_{2,t}\vert {s}_{t}{\mathcal{M}}_{t-1}^{2}),\end{equation} \tag{ 14 }$$

which will become probabilities after normalization.

Procedure for error correction

The probabilities from equation (14) can be understood as describing how likely it is that the system is in each of the eight error states based on the judgment of the model. Whenever |000⟩ does not have the highest probability, we can infer that at least one error has occurred and take the appropriate action to correct it. This procedure, which effectively takes the argmax of the posteriors, can be altered if certain forms of misclassification are more costly than others, or if the act of making a correction itself carries some cost. The procedure can also be modified so that it is more robust to imperfections in the signal, as we do in section 4.3 by introducing the τ_ignore and τ_streak hyperparameters.

Whenever any correction is made, we must update the model with this information by permuting its probabilities to reflect the applied bit-flip. In our example, a correction on the second qubit would lead us to swap the probabilities between pairs of error states which differ in only the second qubit, e.g., |010⟩ ⇌ |000⟩. Without this update the model will continue to recommend the same correction repeatedly, as it does not realize that the state of system has been changed.

A connection can be made between the Bayesian algorithm described here and the maximum likelihood decoder (MLD) commonly used in discrete error correction [34]. Given a specific noise channel and qubit encoding, the MLD is the protocol with the greatest probability of successfully correcting an error, assuming that we have access to projective measurements of the syndromes. The Bayesian model can be viewed as an extension of the MLD to the continuous measurement regime, where the syndrome measurements provide us with incomplete knowledge of the error subspace. As the variance of the Gaussian measurement noise goes to zero, the Bayesian model reduces to the standard MLD protocol for the three-qubit bit-flip code.

Impact of signal imperfections

Compared to thresholding schemes, the Bayesian classifier described here is far more sensitive to the assumptions we make about the noise and error distributions. Such sensitivity can be an advantage, since it allows for near optimal performance when our knowledge of these distributions is accurate.

Of course, when our assumptions about the distributions are wrong, the accuracy of the model can suffer significantly. Out of the three imperfections described in section 3, only the auto-correlation of neighboring samples is directly accounted for in the model. The resonator transients occur over relatively short time intervals, so they are likely to have only a modest impact on the model's performance. The syndrome drift also has a negative impact, as the mean values of the Gaussian distributions are key parameters in the model. If there is a discrepancy between the actual signal means and our pre-programmed values, then every measurement likelihood calculation will be biased.

We explore the size and significance of these effects for all three of our models in section 5.

Neural networks are a subset of the broader family of machine learning methods based on acquiring a learned representation of the data, which consists of parameterized layers of linear transformations and nonlinear activation functions. RNNs are a class of neural network in which the layers connect temporally, combining the previous time step and a hidden representation into the representation for the current time step. They are thus well suited for representation of the time-dependence of continuously measured error syndromes over discrete time steps. Using a training set of labeled signals, the RNN can learn the properties of the weak measurement signal and the structure of the underlying bit-flip channel, which allows it to accurately detect errors as they occur.

The dynamics of a simple recurrent neural network can be expressed by the following equations:

$$\begin{equation*}\begin{aligned}\hfill & {h}_{t}={\sigma }_{h}\left({W}_{h}{x}_{t}+{U}_{h}{h}_{t-1}+{b}_{h}\right),\hfill \\ \hfill & {\,y}_{t}={\sigma }_{y}\left({W}_{y}{h}_{t}+{b}_{y}\right).\hfill \end{aligned}\end{equation*}$$

For each time step t, the network accepts the input vector x_t and, along with the hidden state vector from the previous time step h_t−1, performs a linear transformation parameterized by the weight matrices W_h and U_h and the bias vector b_h before applying a nonlinear activation function given by σ_h . The result is the hidden state vector for the current time step h_t , which is acted upon by an analogous series of operations defined by W_y , b_y and σ_y to produce the output vector y_t . We note that the hidden state h_t effectively encodes a description of the history of inputs ${\left\{{x}_{{t}^{\prime }}\right\}}_{{t}^{\prime }=0}^{{t}^{\prime }=t}$, which therefore allows the network to extract temporal, non-Markovian features from the data.

In our context, we consider the input at each time step to be the vector of measurement signals plus the initial basis state,

$$\begin{equation}{x}_{t}=\left[\begin{matrix}\hfill {I}_{1,t}\hfill \\ \hfill {I}_{2,t}\hfill \\ \hfill {s}_{0}\hfill \end{matrix}\right].\end{equation} \tag{ 15 }$$

Moreover, instead of the standard recurrent neural network architecture, we use a long short-term memory network (LSTM) [35], which is a particular type of recurrent neural network that involves cell states and various gates to evade the vanishing gradient problem of standard RNN architecture [36]. Nevertheless, the same principle underlying the standard function of RNN applies. The output y_t of the LSTM layer is subsequently passed through a dense layer and a softmax activation to produce the posterior probabilities of the eight basis states $p({s}_{t}\vert {\mathcal{M}}_{t}^{k})$, and we select the basis state with the highest posterior as the prediction ${\hat{s}}_{t}$.

Training

Training samples for the RNN require accurate labeling of the states corresponding to the measurement signals at every time step. However, in reality, decoherence effects such as amplitude damping and thermal excitation prevent us from knowing the correct state of the system at some arbitrary time. As a result, to train the RNN, we have to resort to measurement signals with a well defined underlying quantum state. This can be achieved by simulating the measurement signals on states in the absence of unwanted decoherence effects, which will be described in details in section 5. In the simulations, we provide the measurement strength, the single-qubit bit-flip error rate and the initial quantum state as input parameters, and the simulation produces a large number of quantum trajectories to be the training samples of the RNN. We then train the RNN to diagnose bit-flip errors on the three-qubit system, and the trained RNN can be subsequently used to actively correct for errors that occurred. That said, the same information used to generate the training samples is also provided as prior knowledge to the double threshold and the Bayesian model. The two models both require an explicit estimation of the measurement strength as well as the assumption of a certain error rate.

We maximize the likelihood of the RNN parameters on the training set by minimizing the cross entropy batch total loss function, which is defined as

$$\begin{equation}\mathcal{L}=-\frac{1}{NT}\sum\limits _{n=1}^{N}\sum\limits _{t=1}^{T}\mathrm{log}\,{p}_{n}({s}_{t}),\end{equation} \tag{ 16 }$$

where p_n (s_t ) stands for the posterior probability of the true basis state s_t at time step t in the nth sample, while N denotes the mini-batch size and T denotes the total number of steps in each training sample.

To update the parameters to minimize the loss, we perform an iterative training procedure where for each step and parameter w, one applies a gradient descent update of the form $w{\leftarrow}w-\eta (\partial \mathcal{L}/\partial w)$, where the gradients $\partial \mathcal{L}/\partial w$ are computed via backpropagation through the computation graph of the network.

In our experiments, the gradient descent update is preformed using the ADAM optimizer [37]. We adopt a two-layer stacked LSTM with a hidden state size of 32. This small hidden size limits the largest matrix-vector multiplication in computations, hence the memory required, and also limits the number of parameters, facilitating the implementation of the network in real-time experiments. We further provide a comparison test on the performance of different hidden state sizes in appendix C and show that both smaller LSTM and gated recurrent unit (GRU) architecture [38] offer comparable performance for our purpose. The number of stacked layers of the LSTM/GRU and the hyperparameters, such as the batch size in training, are tuned with the assistance of Ray Tune [39].

Re-calibration method for error correction

When performing active error correction, we once again wish to avoid the delay in the posterior probabilities output by the network to reflect the application of an error correction operation C_op on the system. In the case of the Bayesian classifier, we permute the elements of the vector of posterior probabilities, which encodes the state of the model, in accordance with the error correction operation. For the RNN, however, we cannot apply a particular transformation to the hidden state such that the vector of posterior probabilities outputted by the network is permuted in analogous manner, since the function mapping the hidden state to the output vector of posterior probabilities is highly nontrivial.

Any such delay in the network remaining unaware of the quantum state having been corrected is harmful, because another error X_q occurring during this delay, compounding with the correction C_op on the first error, will induce a logical error at the next error correction operation. To see this clearly, considering that the physical qubits are initially in |000⟩, and the first error X₁ results in the state |100⟩. After detecting the error, the model makes a correction that instantly returns the state back to |000⟩. However, the RNN still has the knowledge of the qubits being in |100⟩ until some time later at t_realize before accepting sufficient number of x_t 's that allows it to predict |000⟩. If a second error X₂ occurs before t_realize, the syndromes become (S₁ = −1, S₂ = −1) because the state becomes |010⟩, whereas the RNN, only knowing the state in |100⟩, will eventually predict |101⟩ that has the same syndromes, which is then equivalent to diagnosing an X₃ error. After applying a second error correction C_op = X₃, the physical qubits are now in |111⟩, constituting a logical error. In other words, since we are not capable of injecting the knowledge of a correction operation into the RNN, a correction operation is equivalent to an error seen by the RNN and active correction effectively increases the bit-flip error rate γ in the eyes of the network. Although the correction is correlated with the detected error, the network is generally trained on quantum trajectories with uncorrelated random bit-flip error instances. As will be explained in section 5.1 that a greater γ will induce more logical errors, we conclude that the naive approach of active correction with the RNN suffers from more logical errors.

Therefore, we propose the following re-calibration protocol to effectively hide the action of any error correction operation from the network, so that there is no longer any delay in the posterior probabilities to begin with.

We specifically keep track of all the error correction operations that has been applied up to the present t,

$$\begin{equation*}{N}_{q,t}=\text{Number}\;\;\text{of}\;\;{X}_{q}\;\;\text{corrections}\;\;\text{applied}.\end{equation*}$$

When the measurement signals I_1,t and I_2,t have symmetric noise around their respective mean values and the possible means of I_k,t are always equal and opposite, each C_op correction changes the mean of I_1,t by a factor of −1 if C_op = X₁, changes the mean of I_2,t by a factor of −1 if C_op = X₃, and changes the mean of both I_k,t by a factor of −1 if C_op = X₂. To hide all the corrections done in the past, the measurement signals that are provided as input to the network for all subsequent time steps are then flipped according to N_q,t ,

$$\begin{equation*}\begin{aligned}\hfill {I}_{1,t}^{\prime }={(-1)}^{{N}_{1,t}+{N}_{2,t}}{I}_{1,t},\\ \hfill {I}_{2,t}^{\prime }={(-1)}^{{N}_{2,t}+{N}_{3,t}}{I}_{2,t},\end{aligned}\end{equation*}$$

which we called the re-calibrated signals. From the perspective of the RNN when taking in ${I}_{k,t}^{\prime }$, it appears as if no error correction operation has been applied to the physical qubits.

Given that at some time step we predict a different state ${\hat{s}}_{t}$, we now perform our error correction operation relative to the previous predicted state ${\hat{s}}_{t-1}$.

Adaption to resonator transients for probabilistic models

When the possible means of I_k,t are not equal and opposite, as occurs in the resonator transients upon applying C_op, the re-calibration method breaks down, because flipping the means of either or both I_k,t does not produce the means as if there was no correction applied. A solution to this is to impose an ignore time period τ_ignore right after the correction is applied at some t. During (t, t + τ_ignore], no input x_t is fed into the RNN. As a result, the hidden state of the network is frozen until the ignore time period ends. The re-calibrated signals are accepted by the network only after t + τ_ignore, which reduces the risk of getting incorrect predictions during the transients, but effectively increases the detection time of any error that occurs during the ignore period.

Imposing τ_ignore should be accompanied by a measure to ensure that the RNN diagnoses any error with sufficiently high confidence so that fewer false alarms of error will be followed by an ignore period τ_ignore upon correction. A feasible measure in practice is to determine an error correction operation only if the RNN predicts the same state ${\left\{{\hat{s}}_{{t}^{\prime }}\right\}}_{{t}^{\prime }=t}^{{t}^{\prime }=t+{\tau }_{\text{streak}}}$ for a streak of time steps τ_streak that is different from the old state ${\hat{s}}_{t-1}$, which is a discrete quantity that is easy to optimize. The {τ_ignore, τ_streak} then constitutes a minimal set of tunable hyperparameters for the task of active correction in the presence of resonator transients, which applies to the Bayesian classifier explained in section 4.2 as well.

In quantum memory, it suffices to track the basis states in response to the bit-flip errors that have occurred and only apply error correction operations when needed. We generated 30 000 trajectories of length T = 20 μs from all four simulation schemes with a pre-defined single-qubit error rate as our testing samples, among which are equal portions of trajectories initialized in one of the eight basis states. While the RNN model employed here is trained on 100 000 quantum trajectories from the corresponding simulation scheme, the error rate, noise variance and auto-correlations input to the Bayesian model are also estimated from those quantum trajectories. The tunable parameters in the double threshold model are numerically optimized in schemes with imperfections; the filtering time τ typically lies in the range 0.3–1.6 μs, with larger τ for smaller γ.

In figure 2, we compare the final fidelity $\mathcal{F}={\vert \langle {\psi }_{T}\vert {\psi }_{0}\rangle \vert }^{2}$ against the initial state of the three models in tracking these quantum trajectories subject to bit-flips. The trend is that the final fidelity decreases as a function of the single-qubit error rate γ. This is because the higher the error rate is, the more chances there will be two different bit-flips before the correction to the first bit-flip is made, resulting in a logical error upon the correction, and therefore a lower final fidelity. For instance, a state starting at |000⟩ is flipped to |001⟩ at t₁ and is later also flipped to |011⟩ at t₂ > t₁, such that t₂ is smaller than t₁ + t_detect where t_detect is the detection time of the first error. Subsequently, the model perceiving syndromes with (S₁ = −1, S₂ = +1) will eventually make a C_op = X₁ correction and change the state to |111⟩, leading to a logical error. From the above argument, it is also evident that a shorter detection time is beneficial.

Zoom In Zoom Out Reset image size

From figure 2, we see that the RNN and the Bayesian classifier outperform the double threshold in all simulation schemes, whereas the RNN approximates the Bayesian classifier in all schemes. As discussed in section 4.2, the Bayesian classifier is the optimal model of the three in schemes A and B where there are only auto-correlations in the signals, which is validated in this task. The fact that their performances in schemes C and D are very similar to that in scheme B indicates that the resonator transient pattern and the drifting of the means do not have a significant effect on all three models.

It is reasonable that the drift has a small negative effect to the two probabilistic models, since the drift is usually on the order of the separation of mean values of the two parities, which is in turn one order of magnitude smaller than the standard deviation of the noise. The large noise variance obscures the drifting means, making the drifted signals appear like more noisy signals with fixed means.

Although the models are motivated by correcting bit-flip errors, they can also be exploited in extending the T₁ time of the logical qubit in |1⟩_L = |111⟩. For this task, actively correcting the state is required as opposed to merely tracking the state. While for practical purpose the RNN model is trained on 30 000 quantum trajectories under bit-flips with a length of T = 120 μs, the Bayesian model, whose parameters are estimated from the same set of trajectories, uses a different transition matrix generated by Q' shown in equation (E1) which takes into account the asymmetric probabilities of transitions between the ground and excited state. The parameters for the double threshold model is numerically optimized on the same set of quantum trajectories.

For the three-qubit system initialized to the fully excited state |111⟩, we inspect the population within a Hamming distance 1 away from the initial state, i.e., the population P_exc of the set of basis states {|111⟩, |110⟩, |101⟩, |011⟩}, since these states can be recovered to the initial state by a majority vote. We compare this P_exc against the population of the excited state |1⟩ of a bare qubit as a function of time in all four simulation schemes, and the results are shown in figure 3. In all schemes, the encoded three-qubit system P_exc decays much slower under active correction by any of the three models than the bare qubit excited state population. In all schemes, both the Bayesian and the RNN-based model outrun the double threshold model.

Zoom In Zoom Out Reset image size

Similar to the task of extending the T₁ time of the state |1⟩_L, here we employ the three models to protect the initial state |1⟩_L from bit-flips. Shown in figure 4, we compare the population P_exc of the three-qubit system against the excited population of the bare qubit in time. For schemes A and B, both the Bayesian and the RNN-based model have an advantage over the double threshold. Furthermore, in figure 4 we extract the initial logical error rate Γ_L as a function of γ by computing the time derivative of P_exc at 9.6 μs at each γ. In either scheme with any of the three models, Γ_L scales approximately quadratic in γ, and we can see a strong suppression of Γ_L relative to a bare qubit or the uncorrected three qubits. We remark that, by introducing feedback based on noisy weak measurements, any correction protocol can underperform a majority vote on the encoded qubits without error correction at sufficiently small γ or runtime.

Zoom In Zoom Out Reset image size

To better understand the performance of the models in this important task, we analyze the detection time spent in true positive detection as well as the number of false alarms when the three-qubit system is in |1⟩_L. The difference between a true positive and a false alarm is illustrated in figure 5, which shows the actual and predicted states of the system when an X₃ error occurs and when the model falsely detects an X₁ error. When a true error occurs, the system remains in the corresponding error subspace for a duration determined by the detection time of the model, after which the error is corrected. By contrast, when the model falsely detects that an error has occurred due to measurement noise, it improperly applies a bit-flip to the system and thus pushes it out of the code subspace. After more measurements are recorded, the model determines that the system is in an error subspace and fixes its mistake by applying another bit-flip.

Zoom In Zoom Out Reset image size

As explained in section 5.1, a shorter detection is favorable and will lead to better error corrections, whereas here we can expect more frequent false alarms arises for models with a shorter detection time as a trade off, since the model is prone to make a correction. This is demonstrated in figure 6, where we can see that the best two models, the Bayesian and the RNN-based, both have a shorter detection time and more frequent false alarms at the same time. Nevertheless, for both of these two models, the overall frequency of all false positive detection remains low and is on the order of 0.1 μs⁻¹.

Zoom In Zoom Out Reset image size

Having demonstrated a clear advantage using the RNN-based protocol for tasks in the quantum memory setting over the double threshold protocol, we now study the performance of our protocol for quantum annealing, using a time-dependent Hamiltonian that does not non-commute with the bit-flip errors. We note that the protocol is also applicable to evolution under quantum gate operations.

In quantum annealing, it is imperative to perform error diagnosis and correction in a manner that is both fast and accurate, in order to avoid accruing these logical errors while single bit-flip errors are being diagnosed and corrected. This is because the action of an error X_q effectively transforms the Hamiltonian from H(t) to X_q H(t)X_q in the Heisenberg picture. Until the error is properly diagnosed and corrected, subsequent coherent evolution of the logical state in the code subspace is due to the modified Hamiltonian X_q H(t)X_q . If the original Hamiltonian does not commute with the error, i.e. X_q H(t)X_q ≠ H(t), then such evolution will be spurious rather than as originally intended, causing logical errors to accrue.

For this simulated experiment (see appendix B), the annealing Hamiltonian with a strength Ω₀ evolving ${\rho }_{0}=\vert {\psi }_{0}\rangle \langle {\psi }_{0}\vert ,\vert {\psi }_{0}\rangle =(\vert {0\rangle }_{\text{L}}+\vert {1\rangle }_{\text{L}})/\sqrt{2}$ is chosen to be

$$\begin{equation}H(t)=-{{\Omega}}_{0}\left[a(t){X}_{1}{X}_{2}{X}_{3}+b(t)\frac{{Z}_{1}+{Z}_{2}+{Z}_{3}}{3}\right],\end{equation} \tag{ 17 }$$

where a(t) = 1 − t/T and b(t) = t/T. In the code subspace, it is equal to

$$\begin{equation}h(t)=-{{\Omega}}_{0}\left[a(t){\sigma }_{x}+b(t){\sigma }_{z}\right],\end{equation} \tag{ 18 }$$

whereas in any error subspace it is equal to the spurious Hamiltonian,

$$\begin{equation*}{h}_{\text{spurious}}(t)=-{{\Omega}}_{0}\left[a(t){\sigma }_{x}+b(t)\frac{{\sigma }_{z}}{3}\right].\end{equation*}$$

We adopt the reduction factor [21] as the metric for evaluating the model performance, which is defined as,

$$\begin{equation}\mathcal{R}=\frac{1-{\mathcal{F}}_{\text{une}}}{1-\mathcal{F}},\end{equation} \tag{ 19 }$$

whose numerator is the final infidelity of an unencoded bare qubit initialized to |0⟩ under the annealing Hamiltonian equation (18), and whose denominator is the final infidelity of the three-qubit encoded state in the code subspace with respect to the target quantum state. As $\dot{a}(t),\dot{b}(t)\to 0$, the target quantum state becomes the ground state of the target Hamiltonian.

As shown in figure 7, at relatively low γ, the Bayesian model achieves the highest reduction factor in scheme A, while both the Bayesian and the RNN-based model outperform the double threshold. However at sufficiently high error rates γ, the encoded qubits under active correction using any of the three models show no improvement over a single unencoded qubit, as expected.

Zoom In Zoom Out Reset image size