# Simulation of the five-qubit quantum error correction code on
superconducting qubits

I. A. Simakov  Moscow Institute of Physics and Technology, 141701 Dolgoprudny,
Russia  Russian Quantum Center, 143025 Skolkovo, Moscow, Russia  National
University of Science and Technology ‘‘MISIS’’, 119049 Moscow, Russia  I. S.
Besedin  Russian Quantum Center, 143025 Skolkovo, Moscow, Russia  National
University of Science and Technology ‘‘MISIS’’, 119049 Moscow, Russia  A. V.
Ustinov  Russian Quantum Center, 143025 Skolkovo, Moscow, Russia  National
University of Science and Technology ‘‘MISIS’’, 119049 Moscow, Russia
Physikalisches Institut, Karlsruhe Institute of Technology, 76131 Karlsruhe,
Germany

(August 27, 2024)

###### Abstract

Experimental realization of stabilizer-based quantum error correction (QEC)
codes that would yield superior logical qubit performance is one of the
formidable task for state-of-the-art quantum processors. A major obstacle
towards realizing this goal is the large footprint of QEC codes, even those
with a small distance. We propose a circuit based on the minimal distance-3
QEC code, which requires only 5 data qubits and 5 ancilla qubits, connected in
a ring with iSWAP gates implemented between neighboring qubits. Using a
density-matrix simulation, we show that, thanks to its smaller footprint, the
proposed code has a lower logical error rate than Surface-17 for similar
physical error rates. We also estimate the performance of a neural network-
based error decoder, which can be trained to accommodate the error statistics
of a specific quantum processor by training on experimental data.

##  I Introduction

The realization of quantum error correction codes is an essential step towards
scalable quantum computation. Stabilizer codes Shor (1995); DiVincenzo and
Shor (1996); Gottesman (1998); Steane (1996a, b) encode ‘‘logical’’ qubits in
‘‘physical’’ qubits. If the error rate of a physical qubit is less than a
certain threshold level, then the composite logical qubit will have a lower
error rate Knill et al. . Among the plethora of different QEC codes that have
been proposed, topological QEC codes, such as surface codes Bravyi and Kitaev
; Fowler et al. (2012a); Bonilla Ataides et al. (2021) and color codes Bombin
and Martin-Delgado (2006, 2007) particularly stand out. The topological
property means that there exists an arrangement of the physical qubits on a
finite-dimensional lattice such that only local operators need to be measured
each error correction cycle. This property is important for practical
realization on a number of physical platforms, such as superconducting qubits,
where high-fidelity two-qubit gates are often only available between
neighboring qubits. As a consequence, many recent experiments and near-term
proposals for stabilizer codes on superconducting qubits involve small-scale
surface codes Córcoles A. (2015); Versluis et al. (2017); Varbanov et al.
(2020); AI (2021); Andersen et al. (2020); Marques et al. .

Performance of error correction codes depends not only on the frequency of
errors, but also on the type of errors. Simulation of the performance of
different QEC codes relies on an accurate and physically motivated error
model, which can be computationally expensive, but still viable for small-
scale circuits, such as Surface-17 O’Brien et al. (2017); Tomita and Svore
(2014). For larger codes, approximate error models have been proposed Bravyi
et al. (2017); Ghosh et al. (2012).

Furthermore, the problem of an optimal quantum error decoder remains relevant
Fowler et al. (2012b); Varsamopoulos et al. (2020). One of the most reliable
conventional decoding algorithms is the minimum weight-perfect matching (MWPM)
algorithm, which shows good performance for different codes and in particular
for surface codes O’Brien et al. (2017); Tomita and Svore (2014); Das et al. .
Accuracy of the MWPM decoder relies on an approximate error model, which may
or may not be justified for real devices. The key advantage of MWPM decoders
is speed, as the runtime of the decoder scales polynomially with the size of
the circuit. This makes MWPM decoders good candidates for large-scale error
correction, but questionable for early small-scale codes with qubits and gate
fidelities barely over the error correction threshold. Artificial neural
network-based decoders are another promising type of decoders Sheth et al.
(2020); Baireuther et al. (2019). Their runtime and accuracy varies, and
thanks to their architecture they can accommodate error statistics directly
from measured data Baireuther et al. (2018), demonstrating better logical
error rates compared to MWPM.

Even though surface codes and color codes have numerous important properties
that make them good candidates for scalable universal quantum computation,
they are by far not the most efficient codes in terms of code distance per
physical qubit. Unfortunately, the stabilizers in these more efficient codes
have a complex structure, and designing a circuit to measure them with high
fidelity has proven to be a challenging task. The smallest possible distance-3
code that corrects all possible single-qubit errors is the 5-qubit ZXXZ code.
Schuch and Siewert made a proposal of a quantum circuit that utilizes 9 qubits
connected in a ring topology with iSWAP gates available on neighboring qubits
Schuch and Siewert (2003) implementing this code.

In this paper, we propose a quantum circuit for the 5-qubit distance-3 QEC
code appropriate for experimental realization on a ring of 10 superconducting
qubits. We estimate the upper bound performance of the code using a time-
domain density-matrix simulation of the circuit under a realistic noise model
and compare its performance to the surface code of the same distance. For
forthcoming experiments, we present a neural network decoder based on long
short-term memory architecture, to be trained on only experimentally available
data.

Figure 1: a) The quantum circuit for the simulated state preservation
experiment for the five-qubit code. The scheme is divided into three stages:
initialization of the data qubits, repetitive error correction cycles, and
final measurement of the data qubits. During the correction cycle there are
four series of simultaneous iSWAP gates between the neighboring qubits,
highlighted with the dashed purple lines. b) Schematic layout of computational
qubit positions during the correction cycle. Qubits that are paired in the
same iSWAP gate at each step are outlined with a dashed rectangle. The
numbered dashed purple lines correspond to the ones in panel (a).

##  II Methods

In this section we describe the quantum circuit of the five-qubit code and
give the details of the simulated experiment.

###  II.1 Quantum circuit

The five-qubit QEC code encodes 1 logical qubit in 5 data qubits, protecting
it against arbitrary single-qubit errors. Its four syndrome generators are

$\displaystyle g_{0}$ $\displaystyle=Z_{0}X_{1}X_{2}Z_{3}I_{4}$ (1)
$\displaystyle g_{1}$ $\displaystyle=I_{0}Z_{1}X_{2}X_{3}Z_{4}$ $\displaystyle
g_{2}$ $\displaystyle=Z_{0}I_{1}Z_{2}X_{3}X_{4}$ $\displaystyle g_{3}$
$\displaystyle=X_{0}Z_{1}I_{2}Z_{3}X_{4}$

and Pauli operators on a logical qubit are expressed by the following single-
qubit operators:

$\displaystyle X_{L}$ $\displaystyle=X_{0}X_{1}X_{2}X_{3}X_{4}$ (2)
$\displaystyle Y_{L}$ $\displaystyle=Y_{0}Y_{1}Y_{2}Y_{3}Y_{4}$ $\displaystyle
Z_{L}$ $\displaystyle=Z_{0}Z_{1}Z_{2}Z_{3}Z_{4},$

where $X_{i}$, $Y_{i}$, $Z_{i}$, $I_{i}$ are Pauli operators applied to the
$i$-th qubit.

The quantum circuit for this experimental procedure is presented in Fig. 1a.
It consists of ten physical qubits, five of which are data qubits and the rest
are ancillary qubits used to extract the error syndromes. The circuit
operation can be divided into three stages: data qubit initialization, $k$
repetitions of the error correction cycle, and measurement of the data qubits.

In the initialization stage, we assume that each qubit can be prepared in the
$|0\rangle$ state. To initialize other states, a single-qubit gate denoted as
$P=X,H,ZH,SH\text{ or }ZSH$ can transform the $|0\rangle$ state into
$|1\rangle,|+\rangle,|-\rangle,|i+\rangle,\text{ or }|i-\rangle$,
respectively.

The syndrome measurement stage consists of $k$ cycles. Each cycle starts with
the initialization of the ancillary qubits. The most important parts of the
cycle are four series of iSWAP gates that entangle neighboring qubits. The
action of the iSWAP gate on two qubits can be decomposed into a swap of their
computational states, a CZ gate, and single-qubit $S$ gates on both qubits.
Thereby, it appears that data qubits move clockwise in the ring and ancillas
move counterclockwise, as it is illustrated in Fig 1b. In this cartoon we
schematically show the position of computational qubits before and after each
series of iSWAP gates. Qubits that are paired in the same iSWAP gate at each
step are outlined with a dashed rectangle. At the end of the correction cycle
we measure all ancillas and reset them for the next error correction cycle. It
should be noted that one ancilla measurement is redundant. The extra syndrome
$g_{4}=X_{0}X_{1}Z_{2}I_{3}Z_{4}$ can be expressed as the product of the other
syndromes.

Finally, after $k$ QEC cycles we act on the data qubits with the same gate $P$
as in the initializing procedure and projectively measure them in the $Z$
basis.

In theory the five-qubit code, as the majority of other QEC codes, protects
the logical state against one arbitrary single-qubit error that happens to
data qubit between the correction cycles. It is vulnerable for measurement
errors, initialization errors, multi-qubit errors and imperfection of gates.
But, fortunately, the impact of these unwanted errors can be significantly
reduced by using a smart decoding algorithm that will have access to all
measured data over several sequential correction cycles.

###  II.2 Simulated experiment

To quantify the performance of the five-qubit code we focus on simulation of
the state preservation experiment. The experimental procedure is the
following: we initialize data qubits in one of the logical states
$|0\rangle,|1\rangle,|+\rangle,|-\rangle,|i+\rangle,|i-\rangle$, execute the
syndrome measurement gates for $k$ cycles, and then measure the final state.

For the simulation we calculate time-domain density-matrix evolution and use
operator-sum representation for gates and error channels. In our model we take
into account amplitude and phase damping during qubit idling and gates, and
consider the phenomenological depolarizing effect and stochastic errors of the
rotation gate phase in order to resemble gates implemented in Rol et al.
(2017); Quintana et al. (2017); Yan et al. (2016). Also, we take into account
measurement errors that are similar to the ones observed in a real experiment
Bultink et al. (2016). The main time parameters used in the simulation are
given in Table 3. We calculate the gate fidelities with the formula Nielsen
(2002)

$F_{\text{gate}}=\frac{\text{Tr}(R_{\text{ideal}}^{\dagger}R)+d}{d(d+1)},$ (3)

where $d=2$ for single-qubit gates, $d=4$ for two-qubit gates, and
$R_{\text{ideal}}$ and $R$ are Pauli Transfer matrices corresponding to the
actions of the ideal and simulated gates. The resulting fidelities in our
simulation are 0.9995 for single-qubit rotation gates and 0.998 for two-qubit
iSWAP gates. More detailed information about the error models is given in the
Appendix B.

Figure 2: Comparison of the ideal (upper bound) decoder performance of (a) the
five-qubit code and (b) Surface-17 for different initial states of the data
qubits. The dashed curve is the least-squares fit of the corrected state
fidelity with Eq. (5). The inset plot in figure (a) shows the decay of the
exact logical states $|0\rangle$ and $|1\rangle$ of the five-qubit code.

##  III Results

In the current section we estimate the best possible performance of the five-
qubit code, compare it to surface-17 code and present the logical error rate
dependence on the qubit coherence. Further, we suggest a budding decoder based
on an artificial neural network and test it on a problem of an arbitrary
initial state correction.

###  III.1 Upper bound performance

The decoder is a function that gets the ancilla qubit measurement results as
an input, and returns a unitary operator $U$ that should be applied to the
final logical state to reconstruct the initial state as its output. The ideal
decoder returns the unitary operator that yields the best possible fidelity of
the final state with respect to the pure initial state. This fidelity cannot
exceed the largest eigenvalue of the final density matrix, that gives a
decoder-independent upper bound for the error correction code performance.

Using the density matrix of the final state of the 5 data qubits $\rho$ we
construct the density matrix of the logical qubit $\rho_{L}$:

$\rho_{L}=\cfrac{I+\text{Tr}(X_{L}\rho)X+\text{Tr}(Y_{L}\rho)Y+\text{Tr}(Z_{L}\rho)Z}{2}.$
(4)

The time dependence of the recovered state fidelity is presented in Fig. 2a.
The data points are averages of upper bound fidelities in independent runs of
the density matrix simulation with different ancilla measurement results. For
each of the six initial states $|0\rangle,|1\rangle,|+\rangle$,
$|-\rangle,|i+\rangle$ and $|i-\rangle$ 10000 simulations have been run. The
simulated points are approximated with the function O’Brien et al. (2017)

$F(t)=\frac{1}{2}+\frac{1}{2}(1-2\varepsilon)^{t-t_{0}},$ (5)

where $\varepsilon$ is the error rate per $1\,\mu$s, $t$ is time in $\mu$s.
The parameters $\varepsilon$ and $t_{0}$ are obtained from a the least squares
fit and their values are presented in Table 1. As expected for the 5-qubit
code, the state preservation fidelities of the six select states do not differ
from each other significantly.

| Error rate $\varepsilon$ in % per $1\,\mu$s
---|---
Initial state | 5-qubit | Surface-17
$|0\rangle$ | $0.644\pm 0.006$ | $0.604\pm 0.005$
$|1\rangle$ | $0.629\pm 0.006$ | $0.621\pm 0.005$
$|+\rangle$ | $0.675\pm 0.002$ | $0.466\pm 0.007$
$|-\rangle$ | $0.665\pm 0.004$ | $0.442\pm 0.004$
$|i+\rangle$ | $0.638\pm 0.003$ | $1.143\pm 0.020$
$|i-\rangle$ | $0.623\pm 0.005$ | $1.123\pm 0.015$
Mean | $0.646\pm 0.005$ | $0.703\pm 0.009$
Table 1: Upper bound error rates $\varepsilon$ as obtained from least-squares
fitting of fidelity upper bounds shown in Fig. 2 with Eq. 5 for the six
different initial states in the five-qubit QEC code and Surface-17. Errors
estimates for $\varepsilon$ are computed from the variance matrix of the least
squares method. ‘‘Mean’’ error rate corresponds to the average over the six
initial states.

The first point of each time dependence is excluded from the fit, as it
clearly does not lie on the same curve as the others. The reason for the
abrupt fidelity drop after the first cycle is that our proposed logical qubit
initialization procedure is flawed. While it does prepare an eigenstate of the
logical qubit’s operators with only single-qubit gates, preparing a state that
is also an eigenstate of the stabilizers requires a complex circuit, which is
challenging to implement with high fidelity in the proposed device topology.
As a result, during the first cycle the ancilla measurements randomly project
the system into one of the stabilizer eigenstates, and actual error correction
starts only from the second cycle.

If we numerically prepare an eigenstate of the error correction code instead
of the described initialization sequence, the curve starts from the first
cycle, as shown in the inset plot in Fig. 2a for the $|0\rangle_{L}$ and
$|1\rangle_{L}$ states. The fitted error rates $\varepsilon$ are $0.598\pm
0.006$ and $0.588\pm 0.009$% per $1\,\mu$s.

To evaluate how well the five-qubit code fares in compare to other stabilizer
codes, we simulate Surface-17 under the same noise conditions (Appendix C).
The time dependence of the recovered state fidelity for Surface-17 averaged
over 5000 independent runs per each of the six initial states is presented in
Fig. 2b and in Table 1. In contrast to the five-qubit code it appears that
different states are preserved with considerably different fidelity, which
leads to worse preservation of arbitrary states on average.

Figure 3: Upper bound fidelities of the five-qubit (blue) and Surface-17
(green) error correction codes, averaged over the six initial states
$|0\rangle,|1\rangle,|+\rangle$, $|-\rangle,|i+\rangle$ and $|i-\rangle$,
calculated via the density-matrix simulation. The grey curve corresponds to
the average state preservation fidelity in a single physical qubit.

Also we compare the both correction codes to the single qubit relaxation with
times $T_{1}$ and $T_{2}$. In this case the fidelity averaged over the six
initial states is given by the formula

$F_{\text{single}}(t)=\frac{1}{6}\left(1+e^{-t/T_{1}}\right)+\frac{1}{3}\left(1+e^{-t/T_{2}}\right).$
(6)

In Fig. 3 we plot the mean state preservation fidelity of the considered QEC
codes over the six cardinal states.

Figure 4:  Upper bound logical error rate $\varepsilon$ of the five-qubit code
in presence of only relaxation and dephasing ($T_{1}=T_{2})$, averaged over
the three initial states $|0\rangle,|+\rangle$ and $|i+\rangle$. The blue
points are data calculated via the density-matrix simulation; blue curve is
their fit with an Eq. 7; orange line corresponds to a physical error rate of a
single qubit.

A noteworthy question is the dependency of the logical error rate on the
quality of the qubits. To study this issue, we simplify the error model and
take into account only the decoherence of qubits described by the relaxation
times $T_{1}$ and $T_{2}$. It is essential to note that qubits relax not only
between correction cycles, when the ancillas are being measured but also
during the whole correction procedure. We simulate the state preservation
experiment for several times $T_{1}$ (the transverse relaxation time $T_{2}$
is chosen equal to $T_{1}$) and obtain the corresponding upper-bound logical
error rate $\varepsilon$) shown in Fig. 4.

To approximate the data we suggest following simple model. Let $N$ qubits
participate in the experiment and $p$ is an independent probability of a
single-qubit error per unit of time. Then the probability that more than two
single-qubit errors happen during a correction cycle can be estimated by the
following expression:

$P=1-\left[(1-p)^{N}+Np(1-p)^{N-1}\right],$ (7)

that for small $p$ it can be expressed in the series:

$P=\frac{(N-1)N}{2}p^{2}+\mathcal{O}(p^{3}).$ (8)

Since the five-qubit code protects the logical state against one single-qubit
error, then it is reasonable to approximate the dependence of the logical
error rate on time $T_{1}$ by formula 7, where $p$ is intuitively proportional
to $p=\beta\exp(-t/T_{1})$ with some proportionality coefficient $\beta$.
Optimizing this parameter $\beta$ ($\beta_{\mathrm{fit}}=0.44$) it turns out
that the theoretical curve fits well the simulated values at $T_{1}\geqslant
20\mu$s. For smaller $T_{1}$ some divergence is observed, which is explained
by the fact that the model described above can be valid only for small $p$.

In addition, comparing the formulas 5 and 6, one can see that in case
$T_{1}=T_{2}$ logical error rate per $1\mu$s for a single qubit is
$\varepsilon=0.5\exp(-1/T_{1})$ [$T_{1}$ in $\mu$s] and, therefore, consider
it as a physical error rate.

Figure 5: Architecture of the neural network decoder. The network takes
measurements of the five ancillary qubits during $k$ cycles as an input,
processes them with 5 LSTM layers with hidden size 256, then passes them
through a perceptron with 2 hidden layers of size 64, and output size of 4
neurons, and ReLU activation functions. After a final softmax activation
function, the outputs are interpreted as probability predictions for each
Pauli operator error $p(X),p(Y),p(Z)$ or $p(I)$.

###  III.2 Neural network decoder

The upper bound is the theoretical limit of any decoder performance for a
given error model. Realizing this decoder amounts to simulating the system’s
evolution with this error model, which can be challenging for large systems
and systems with complex error statistics.

Artificial neural network decoders are an alternative approach altogether. The
problem of any decoder is to find and analyze patterns in the sequences of
ancillary qubit measurements. Instead of programming an algorithm that would
detect patterns that physical errors produce in the ancilla measurements, one
can train a neural network to detect these patterns from a training dataset.
Due to the nature of the errors, data about the ancilla measurements from all
previous steps may contain information about which logical error has which
probability on each step. This time series-like type of feature data can be
effectively accommodated by recursive neural network architectures, such as
the Long Short-Term Memory (LSTM) Hochreiter and Schmidhuber (1997). The
performance of LSTM-based neural network decoders has been studied for color
codes Baireuther et al. (2019) and surface codes Baireuther et al. (2018); the
state preservation fidelity obtained with these decoders is on par or better
than with more traditional decoders.

These inputs are fed into a five-layer LSTM network with a hidden layer size
of 256 neurons Fig. 5. The LSTM’s outputs are passed on to a perceptron with
two hidden layers with 64 neurons, an output layer of 4 neurons, and a ReLU
activation function. After a final softmax activation layer, the neural
network outputs are interpreted as probability predictions for each Pauli
operator error $p(X),p(Y),p(Z)$ or $p(I)$ (no error). Our neural network
implementation is based on the PyTorch framework.

A training dataset contains the simulated ancilla measurement results and
readout results of the final data qubit states. The neural network is trained
simultaneously with three training datasets with different initialization and
projection operators $P$ corresponding to three different cardinal points of
the Bloch sphere $\\{|+\rangle,|i+\rangle,|0\rangle\\}$, with $4\times 10^{6}$
readout results in each dataset.

From formula (2) it follows that the readout result of the logical qubit’s
Pauli operators $\Pi_{L}\in\\{X_{L},Y_{L},Z_{L}\\}$ can be expressed as

$\Pi_{L}=\prod\limits_{i=0}^{4}\Pi_{i}.$ (9)

If $\Pi_{L}=+1$, then effectively the qubit has not flipped, and the logical
error operator is either $I$ or $\Pi$. Otherwise, it is one of the other two
Pauli operators. We define a loss function based on the cross-entropy of the
experimentally measured logical error operator and the neural network outputs
over the set of logical error operators as

$L=\frac{1}{2}\begin{cases}\log p(I)+\log p(\Pi),&\text{ if }\Pi_{L}=+1,\\\
\sum\limits_{\mathcal{P}\in\\{X,Y,Z\\}/\\{\Pi\\}}\log p(\mathcal{P}),&\text{
otherwise}.\end{cases}$ (10)

Figure 6: Comparison of recovered logical qubit state fidelity in the five-
qubit QEC code for a perfect decoder (upper bound) and the neural network
decoder applied to $|0\rangle,|+\rangle,|i+\rangle$ initial states and to
arbitrary uniformly distributed pure states.

The decoder performance is tested on 150000 runs from the same initial states,
but as opposed to the training procedure we do not measure the final logical
state but just calculate it from the density matrix of data qubits via the
formula 4 for better comparison with the upper bound decoder. The fidelity of
the neural network decoder-recovered state is shown in Fig. 6 (violet curve)
and the fitted error rate is given in Table 2.

However, the goal is to correct not only the six cardinal points on the Bloch
sphere but an arbitrary state. To check how well the decoder fares in this
case we generate a random pure single-qubit state
$\alpha|0\rangle+\beta|1\rangle$, uniformly distributed on the Bloch sphere,
and encode it into the state $\alpha|00000\rangle+\beta|11111\rangle$ by the
circuit shown in Fig. 7. We assume that the encoding circuit is error-free.
After initialization, we simulate the state preservation experiment and apply
the neural network decoder to the final density matrix. We repeat the whole
procedure for 17000 different initial states and plot the average fidelity in
Fig. 6 (pink curve). Taking into account the error, we observe approximately
the nearly same $\varepsilon$ as for the cardinal points (see Table 2).

|  $\varepsilon$, % per $1\,\mu$s
---|---
Initial state | Neural Network Decoder
$|0\rangle$ | $0.88\pm 0.01$
$|+\rangle$ | $1.03\pm 0.01$
$|i+\rangle$ | $1.05\pm 0.01$
Mean over |
$|0\rangle,|+\rangle,|i+\rangle$ states | $0.99\pm 0.01$
Arbitrary states | $0.993\pm 0.003$
Table 2: Error rates $\varepsilon$ obtained with least-squares fitting of the fidelity time dependence for the neural network decoder shown in Fig. 6 for different initial states. The error is computed via the k-fold cross-validation method with $k=5$. Figure 7: The quantum circuit for arbitrary initial state preparation consisting of four CNOT gates that encodes $\alpha|0\rangle+\beta|1\rangle$ state into $\alpha|00000\rangle+\beta|11111\rangle$. Parameter | Symbol | Value
---|---|---
Qubit relaxation time | $T_{1}$ |  $30$ $\mu\text{s}$
Qubit dephasing time | $T_{2}$ |  $40$ $\mu\text{s}$
Hadamard gate duration | $\tau_{1}$ |  $20$ ns
iSWAP gate duration | $\tau_{2}$ |  $40$ ns
Measurement time | $t_{m}$ |  $300$ ns
Reset time | $t_{r}$ |  $300$ ns
5-qubit cycle time | $t_{c}$ |  $840$ ns
Table 3: Time parameters used in the density-matrix simulation Bultink et al.
(2016); Rol et al. (2017).

##  IV Conclusion

In conclusion, we propose a circuit that implements the five-qubit quantum
error correction code on a realistic topology of superconducting qubits. The
key advantage of the scheme is that it requires only 10 qubits connected in a
ring with two-qubit iSWAP gates available between any two neighboring qubits.
The small footprint of our circuit makes it one of the least hardware-
demanding proposals for stabilizer-based error correction. Using a density-
matrix simulation of the state preservation experiment under a realistic noise
model, we obtain estimates for logical error rates. We propose a decoder based
on an artificial neural network, that can be trained on experimentally
available data, which shows satisfactory performance close to the upper bound
for preservation of arbitrary logical states.

## Acknowledgements

This work was performed with the financial support of the Russian Science
Foundation, Project No. 21-72-30026.

## Code availability

The computer code developed for simulations is available from the
corresponding author<EMAIL_ADDRESS>upon request.

## Author contributions

I.A.S. performed simulations and calculations. I.S.B. proposed the stabilizer
circuit and the neural network decoder architecture. A.V.U. supervised the
project. All authors contributed to the manuscript.

## References

  * Shor [1995] Peter W. Shor.  Scheme for reducing decoherence in quantum computer memory.  _Phys. Rev. A_ , 52:R2493–R2496, Oct 1995.  doi: 10.1103/PhysRevA.52.R2493.
  * DiVincenzo and Shor [1996] David P. DiVincenzo and Peter W. Shor.  Fault-tolerant error correction with efficient quantum codes.  _Phys. Rev. Lett._ , 77:3260–3263, Oct 1996.  doi: 10.1103/PhysRevLett.77.3260.
  * Gottesman [1998] Daniel Gottesman.  Theory of fault-tolerant quantum computation.  _Phys. Rev. A_ , 57:127–137, Jan 1998.  doi: 10.1103/PhysRevA.57.127.
  * Steane [1996a] A. M. Steane.  Error correcting codes in quantum theory.  _Phys. Rev. Lett._ , 77:793–797, Jul 1996a.  doi: 10.1103/PhysRevLett.77.793.
  * Steane [1996b] Andrew Steane.  Multiple particle interference and quantum error correction.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 452, 01 1996b.  doi: 10.1098/rspa.1996.0136.
  * [6] Emanuel Knill, Raymond Laflamme, and Wojciech Zurek.  Threshold accuracy for quantum computation.  _arXiv:9610011 (1996)_.
  * [7] S. B. Bravyi and A. Y. Kitaev.  Quantum codes on a lattice with boundary.  _arXiv:quant-ph/9811052 (1998)_.  doi: 10.1098/rspa.1996.0136.
  * Fowler et al. [2012a] Austin G. Fowler, Matteo Mariantoni, John M. Martinis, and Andrew N. Cleland.  Surface codes: Towards practical large-scale quantum computation.  _Physical Review A_ , 86(3), Sep 2012a.  ISSN 1094-1622.  doi: 10.1103/physreva.86.032324.
  * Bonilla Ataides et al. [2021] Juan Pablo Bonilla Ataides, David Tuckett, Stephen Bartlett, Steven Flammia, and Benjamin Brown.  The xzzx surface code.  _Nature Communications_ , 12, 04 2021.  doi: 10.1038/s41467-021-22274-1.
  * Bombin and Martin-Delgado [2006] H. Bombin and M. A. Martin-Delgado.  Topological quantum distillation.  _Phys. Rev. Lett._ , 97:180501, Oct 2006.  doi: 10.1103/PhysRevLett.97.180501.
  * Bombin and Martin-Delgado [2007] H. Bombin and M. A. Martin-Delgado.  Topological computation without braiding.  _Phys. Rev. Lett._ , 98:160502, Apr 2007.  doi: 10.1103/PhysRevLett.98.160502.
  * Córcoles A. [2015] Srinivasan S. et al. Córcoles A., Magesan E.  Demonstration of a quantum error detection code using a square lattice of four superconducting qubits.  _Nature Communications_ , 6:6979, 2015.  doi: 10.1038/ncomms7979.
  * Versluis et al. [2017] R. Versluis, S. Poletto, N. Khammassi, B. Tarasinski, N. Haider, D. J. Michalak, A. Bruno, K. Bertels, and L. DiCarlo.  Scalable quantum circuit and control for a superconducting surface code.  _Phys. Rev. Applied_ , 8:034021, Sep 2017.  doi: 10.1103/PhysRevApplied.8.034021.
  * Varbanov et al. [2020] Boris Varbanov, Francesco Battistel, Brian Tarasinski, Viacheslav Ostroukh, Thomas O’Brien, Leonardo DiCarlo, and Barbara Terhal.  Leakage detection for a transmon-based surface code.  _npj Quantum Information_ , 6, 12 2020.  doi: 10.1038/s41534-020-00330-w.
  * AI [2021] Google Quantum AI.  Exponential suppression of bit or phase flip errors with repetitive error correction.  _Nature_ , 595:383–387, 7 2021.  doi: 10.1038/s41586-021-03588-y.
  * Andersen et al. [2020] Christian Kraglund Andersen, Ants Remm, Stefania Lazar, Sebastian Krinner, Nathan Lacroix, Graham J Norris, Mihai Gabureac, Christopher Eichler, and Andreas Wallraff.  Repeated quantum error detection in a surface code.  _Nature Physics_ , 16(8):875–880, 2020.  doi: 10.1038/s41567-020-0920-y.
  * [17] JF Marques, BM Varbanov, MS Moreira, Hany Ali, Nandini Muthusubramanian, Christos Zachariadis, Francesco Battistel, Marc Beekman, Nadia Haider, Wouter Vlothuizen, et al.  Logical-qubit operations in an error-detecting surface code.  _arXiv preprint arXiv:2102.13071 (2021)_.  doi: 10.21203/rs.3.rs-380732/v1.
  * O’Brien et al. [2017] T. E. O’Brien, B. Tarasinski, and L. DiCarlo.  Density-matrix simulation of small surface codes under current and projected experimental noise.  _npj Quantum Information_ , 3(1), Sep 2017.  ISSN 2056-6387.  doi: 10.1038/s41534-017-0039-x.
  * Tomita and Svore [2014] Yu Tomita and Krysta M. Svore.  Low-distance surface codes under realistic quantum noise.  _Phys. Rev. A_ , 90:062320, Dec 2014.  doi: 10.1103/PhysRevA.90.062320.
  * Bravyi et al. [2017] Sergey Bravyi, Matthias Englbrecht, Robert Koenig, and Nolan Peard.  Correcting coherent errors with surface codes.  _npj Quantum Information_ , 4, 10 2017.  doi: 10.1038/s41534-018-0106-y.
  * Ghosh et al. [2012] Joydip Ghosh, Austin G. Fowler, and Michael R. Geller.  Surface code with decoherence: An analysis of three superconducting architectures.  _Phys. Rev. A_ , 86:062318, Dec 2012.  doi: 10.1103/PhysRevA.86.062318.
  * Fowler et al. [2012b] Austin G. Fowler, Adam C. Whiteside, and Lloyd C. L. Hollenberg.  Towards practical classical processing for the surface code.  _Phys. Rev. Lett._ , 108:180501, May 2012b.  doi: 10.1103/PhysRevLett.108.180501.
  * Varsamopoulos et al. [2020] Savvas Varsamopoulos, Koen Bertels, and Carmen G. Almudever.  Decoding surface code with a distributed neural network based decoder.  2, 2020.  doi: 10.1007/s42484-020-00015-9.
  * [24] Poulami Das, Christopher A Pattison, Srilatha Manne, Douglas Carmean, Krysta Svore, Moinuddin Qureshi, and Nicolas Delfosse.  A scalable decoder micro-architecture for fault-tolerant quantum computing.  _arXiv preprint arXiv:2001.06598 (2020)_.
  * Sheth et al. [2020] Milap Sheth, Sara Zafar Jafarzadeh, and Vlad Gheorghiu.  Neural ensemble decoding for topological quantum error-correcting codes.  _Physical Review A_ , 101(3):032338, 2020.  doi: 10.1103/PhysRevA.101.032338.
  * Baireuther et al. [2019] Paul Baireuther, Marcello D Caio, Ben Criger, Carlo WJ Beenakker, and Thomas E O’Brien.  Neural network decoder for topological color codes with circuit level noise.  _New Journal of Physics_ , 21(1):013003, 2019\.  doi: 10.1088/1367-2630/aaf29e.
  * Baireuther et al. [2018] Paul Baireuther, Thomas E O’Brien, Brian Tarasinski, and Carlo WJ Beenakker.  Machine-learning-assisted correction of correlated qubit errors in a topological code.  _Quantum_ , 2:48, 2018.  doi: 10.22331/q-2018-01-29-48.
  * Schuch and Siewert [2003] Norbert Schuch and Jens Siewert.  Natural two-qubit gate for quantum computation using the xy interaction.  _Physical Review A_ , 67(3):032301, 2003.  doi: 10.1103/PhysRevA.67.032301.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9:1735–80, 12 1997.  doi: 10.1162/neco.1997.9.8.1735.
  * Rol et al. [2017] M. A. Rol, C. C. Bultink, T. E. O’Brien, S. R. de Jong, L. S. Theis, X. Fu, F. Luthi, R. F. L. Vermeulen, J. C. de Sterke, A. Bruno, D. Deurloo, R. N. Schouten, F. K. Wilhelm, and L. DiCarlo.  Restless tuneup of high-fidelity qubit gates.  _Phys. Rev. Applied_ , 7:041001, Apr 2017.  doi: 10.1103/PhysRevApplied.7.041001.
  * Quintana et al. [2017] C. M. Quintana, Yu Chen, D. Sank, A. G. Petukhov, T. C. White, Dvir Kafri, B. Chiaro, A. Megrant, R. Barends, B. Campbell, Z. Chen, A. Dunsworth, A. G. Fowler, R. Graff, E. Jeffrey, J. Kelly, E. Lucero, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, P. Roushan, A. Shabani, V. N. Smelyanskiy, A. Vainsencher, J. Wenner, H. Neven, and John M. Martinis.  Observation of classical-quantum crossover of $1/f$ flux noise and its paramagnetic temperature dependence.  _Phys. Rev. Lett._ , 118:057702, Jan 2017.  doi: 10.1103/PhysRevLett.118.057702.
  * Yan et al. [2016] Fei Yan, Simon Gustavsson, Archana Kamal, Jeffrey Birenbaum, Adam P Sears, David Hover, Ted J Gudmundsen, Danna Rosenberg, Gabriel Samach, Steven Weber, et al.  The flux qubit revisited to enhance coherence and reproducibility.  _Nature communications_ , 7(1):1–9, 2016.  doi: 10.1038/ncomms12964.
  * Bultink et al. [2016] C. C. Bultink, M. A. Rol, T. E. O’Brien, X. Fu, B. C. S. Dikken, C. Dickel, R. F. L. Vermeulen, J. C. de Sterke, A. Bruno, R. N. Schouten, and L. DiCarlo.  Active resonator reset in the nonlinear dispersive regime of circuit qed.  _Phys. Rev. Applied_ , 6:034008, Sep 2016.  doi: 10.1103/PhysRevApplied.6.034008.
  * Nielsen [2002] Michael A Nielsen.  A simple formula for the average gate fidelity of a quantum dynamical operation.  _Physics Letters A_ , 303(4):249–252, 2002.  doi: 10.1016/S0375-9601(02)01272-0.
  * Nielsen and Chuang [2011] Michael A. Nielsen and Isaac L. Chuang.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_.  Cambridge University Press, USA, 10th edition, 2011.  ISBN 1107002176.
  * Ganzhorn et al. [2020] M. Ganzhorn, G. Salis, D. J. Egger, A. Fuhrer, M. Mergenthaler, C. Müller, P. Müller, S. Paredes, M. Pechal, M. Werninghaus, and S. Filipp.  Benchmarking the noise sensitivity of different parametric two-qubit gates in a single superconducting quantum computing platform.  _Phys. Rev. Research_ , 2:033447, Sep 2020.  doi: 10.1103/PhysRevResearch.2.033447.

##  Appendix A Density-matrix simulation details

In order to estimate the limit of the QEC code performance we need to simulate
the dynamics of the investigated open quantum system, described with the
density matrix $\rho$. The time evolution $\rho\rightarrow\mathcal{E}(\rho)$
can be modeled with the operator-sum representation Nielsen and Chuang [2011]:

$\mathcal{E}(\rho)=\sum_{i}E_{i}\rho E_{i}^{\dagger},$ (11)

where $E_{i}$ are Kraus operators:

$\sum_{i=0}^{n}E_{i}^{\dagger}E_{i}=I.$ (12)

For the convenience of calculations we use the Pauli transfer matrix
representation of the superoperator $\mathcal{E}$:

$(R_{\mathcal{E}})_{ij}=\frac{1}{2}\text{Tr}\left(\sigma_{i}\mathcal{E}\left(\sigma_{j}\right)\right),$
(13)

where $\sigma_{0}$ is the identity matrix $I$ and
$\sigma_{1},\sigma_{2},\sigma_{3}$ are Pauli matrices $X,Y,Z$. In a such form
the time evolution, including quantum gates and noise channels, can be
expressed as simple multiplication of the transfer matrix
$(R_{\mathcal{E}})_{ij}$ and vector representation of the density matrix in
the Pauli basis.

The Hilbert space of the quantum system implementing the QEC code (Fig. 1)
consists of 10 qubits. The density matrix of the system is represented as a
vector of dimension $4^{10}$ and transfer matrices, acting on this vector,
have the size $4^{10}\times 4^{10}$. In order to accelerate calculation we use
sparse matrices, batches and GPUs.

##  Appendix B Error models

This appendix contains the mathematical description of the noise channels that
we take into account in the density-matrix simulation. The basic motivation of
all the error models introduced in the paper is to get the simulation closer
to real experiments Bultink et al. [2016], Rol et al. [2017], Quintana et al.
[2017], Yan et al. [2016]. The values of the error model parameters are given
in Table 4, excluding measurement errors which are listed separately.

Parameter | Symbol | Value
---|---|---
Qubit relaxation time | $T_{1}$ |  $30$ $\mu\text{s}$
Qubit dephasing time | $T_{2}$ |  $40$ $\mu\text{s}$
single-qubit gate duration | $\tau_{1}$ |  $20$ ns
two-qubit gate duration | $\tau_{2}$ |  $40$ ns
Measurement time | $t_{m}$ |  $300$ ns
Reset time | $t_{r}$ |  $300$ ns
In-axis rotation error | $p_{\text{axis}}$ | $10^{-4}$
In-plane rotation error | $p_{\text{plane}}$ | $5\cdot 10^{-4}$
Root-mean-square error | $\varphi_{\text{rms}}$ | $0.01$
Table 4: Parameters of error models Bultink et al. [2016], Rol et al. [2017],
Quintana et al. [2017], Yan et al. [2016].

###  B.1 Idling

In the simulation we use the standard model of amplitude and phase damping
channels. The operation elements of these channels are

$E^{\text{AD}}_{0}=\left(\begin{array}[]{cc}1&0\\\
0&\sqrt{1-\gamma_{1}}\end{array}\right),\quad
E^{\text{AD}}_{1}=\left(\begin{array}[]{cc}0&\sqrt{\gamma_{1}}\\\
0&0\end{array}\right)$ (14)

and

$E^{\text{PD}}_{0}=\left(\begin{array}[]{cc}1&0\\\
0&\sqrt{1-\gamma_{\varphi}}\end{array}\right),\quad
E^{\text{PD}}_{1}=\left(\begin{array}[]{cc}0&0\\\
0&\sqrt{\gamma_{\varphi}}\end{array}\right),$ (15)

where $\gamma_{1}=1-e^{-t/T_{1}}$, $\gamma_{\varphi}=1-e^{-2t/T_{\varphi}}$,
$\frac{1}{T_{2}}=\frac{1}{2T_{1}}+\frac{1}{T_{\varphi}}$. Corresponding Pauli
transfer matrices are

$R_{\mathcal{E}_{\text{AD}}}=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\sqrt{1-\gamma_{1}}&0&0\\\ 0&0&\sqrt{1-\gamma_{1}}&0\\\
\gamma_{1}&0&0&1-\gamma_{1}\end{array}\right),$ (16)
$R_{\mathcal{E}_{\text{PD}}}=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\sqrt{1-\gamma_{\varphi}}&0&0\\\ 0&0&\sqrt{1-\gamma_{\varphi}}&0\\\
0&0&0&1\end{array}\right).$ (17)

Thus, the evolution of an idling qubit for the duration $t$ is described by
transfer matrix

$R_{\text{idle}}(t)=R_{\mathcal{E}_{\text{AD}}}R_{\mathcal{E}_{\text{PD}}}.$
(18)

###  B.2 Single-qubit gates

There are three single-qubit gates in the circuit Fig. 1: Hadamard ($H$),
Phase ($S$) and Pauli-Z ($Z$). The last two gates are considered to be
virtual, instantaneous and thereby ideal. The Hadamard gate can be decomposed
into rotations around the $z$-axis and the $y$-axis: $H=ZR_{y}(-\pi/2)$.
Thereby, it is necessary to build an error model for the $y$-axis rotation.
This model should include amplitude and phase damping, depolarizing and
statistical error of the angle of rotation.

Processes of amplitude and phase decoherence are taken into account with the
sandwich model, where the gate transfer matrix is expressed as the product of
the two idling process matrices interleaved with the ideal gate process
matrix.

In order to have better conformity with tomography experiments, we
additionally introduce a depolarization channel, that geometrically
corresponds to squeezing along the $x$ and $z$ axes with the factor
$1-p_{\text{plane}}$ and the $y$ axis with the factor $1-p_{\text{axis}}$. The
depolarization Pauli transfer matrix is

$R_{\mathrm{dep}}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&1-p_{\text{plane
}}&0&0\\\ 0&0&1-p_{\text{axis }}&0\\\ 0&0&0&1-p_{\text{plane
}}\end{array}\right).$ (19)

Finally, we include the statistical error for the angle of rotation $\varphi$.
We consider $\varphi$ to be normally distributed with mean $\varphi_{0}$ and
variance $\sigma^{2}$. The expectation value of the exponential function is

$\langle\mathrm{e}^{i\varphi}\rangle=\int\limits_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}\sigma}\mathrm{e}^{-\frac{(\varphi-\varphi_{0})^{2}}{2\sigma^{2}}}\mathrm{e}^{i\varphi}\mathrm{d}\varphi=\mathrm{e}^{i\varphi}\mathrm{e}^{-\frac{\sigma^{2}}{2}}.$
(20)

Using this result we add stochastic error to the matrix elements of the
rotation gates. These stochastic errors can be attributed to flux noise and
long-term drift of microwave signal source phase and amplitude.
Phenomenological estimates yield an RMS value in the order of
$\sigma=\varphi_{\text{rms}}\approx 0.01$ rad.

According to the model above, the Hadamard gate can be expressed as

$R_{H}=R_{\text{idle}}\left(\frac{\tau_{1}}{2}\right)R_{\mathrm{dep}}R_{R_{z}(\pi)}R^{\prime}_{R_{y}\left(-\frac{\pi}{2}\right)}R_{\text{idle}}\left(\frac{\tau_{1}}{2}\right),$
(21)

where $\tau_{1}$ is the single-qubit gate duration and
$R^{\prime}_{R_{y}\left(-\frac{\pi}{2}\right)}$ is the PTM of the rotation
gate with the stochastic error.

###  B.3 Two-qubit gates

The only two-qubit gate used in the present five-qubit code is an iSWAP gate.
This gate is generated by XY-type interaction, which is natural for many
platforms for quantum computing, including superconducting qubits Ganzhorn et
al. [2020], Schuch and Siewert [2003], and the unitary matrix of such
interaction is

$U(\theta,\eta,\zeta)=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\cos\theta/2&ie^{i\eta}\sin\frac{\theta}{2}&0\\\
0&ie^{-i\eta}\sin\theta/2&\cos\theta/2&0\\\
0&0&0&e^{i\zeta}\end{array}\right).$ (22)

For iSWAP the parameters are $U_{\mathrm{iSWAP}}=U(\pi,0,0)$. Here we assume
the same stochastic error independent and equal for both parameters $\theta$
and $\eta$ as it is described in the section above for the single-qubit
rotation gate. The parameter $\zeta$ takes into account the conditional phase
arising in the dispersive limit due to the
$|11\rangle\leftrightarrow|20\rangle$ interaction in transmon qubits. The
error for $\zeta$ is correlated with the single-qubit stochastic RMS deviation
due to flux noise, such that $\sigma_{\zeta}=\varphi_{\text{rms}}/2=0.005$
rad.

The Surface-17 code uses the controlled-Z (CZ) gate with the following unitary
matrix: $U_{\text{CZ}}=U(0,0,\pi)$. For CZ we take into account stochastic
error only for parameter $\zeta$.

Decoherence is modeled by a sandwich structure, as it was described above for
single-qubit rotation gates, and the final PTM of the two-qubit gate $G$ is

$R_{G}=R_{\text{idle}}\left(\frac{\tau_{2}}{2}\right)R^{\prime}_{G}R_{\text{idle}}\left(\frac{\tau_{2}}{2}\right).$
(23)

Here $\tau_{2}$ is the duration of the gate, and $R^{\prime}_{G}$ is the PTM
of the two-qubit gate with stochastic errors.

###  B.4 Measurement

The measurement procedure is realized by the method, described below, in order
to approximated the experimental results obtained in Bultink et al. [2016]. We
use the model of back-to-back measurements. Firstly, we calculate from the
density matrix the probabilities of a qubit to be projected in $|0\rangle$ and
$|1\rangle$ states and then randomly with respect to these probabilities
choose the initial state $i\in\\{|0\rangle,|1\rangle\\}$. The second step is
to determine the measurement result $m\in\\{+1,-1\\}$ and outcome state
$o\in\\{|0\rangle,|1\rangle\\}$ of the qubit after the measurement procedure.
We match with initial and outcome states according to the experimentally
obtained probabilities given in the Table 5.

Figure 8: Quantum circuit of the QEC cycle of Surface-17 code, consisting of single-qubit rotations about the $y$-axis and eight series of simultaneous CZ gates. Probability | Value | Probability | Value
---|---|---|---
$\varepsilon_{0}^{+1,0}$ | $0.9985$ | $\varepsilon_{1}^{+1,0}$ | $0.0050$
$\varepsilon_{0}^{+1,1}$ | $0.0000$ | $\varepsilon_{1}^{+1,1}$ | $0.0015$
$\varepsilon_{0}^{-1,0}$ | $0.0015$ | $\varepsilon_{1}^{-1,0}$ | $0.0149$
$\varepsilon_{0}^{-1,1}$ | $0.0000$ | $\varepsilon_{1}^{-1,1}$ | $0.9786$
Table 5: Probabilities of different errors in measurement procedure O’Brien et
al. [2017], Bultink et al. [2016].

##  Appendix C Surface-17

The simulated experimental procedure for the surface code simulation is
similar to the five-qubit code simulation procedure. It can be divided into
three stages: data qubit preparation, correction cycles and data qubit
measurement. The first and the third stages are the same as in the five-qubit
code with the only difference being the amount data qubits. The quantum
circuit for the QEC cycle of Surface-17 is presented in Fig. 8 O’Brien et al.
[2017]. It consists of single-qubit rotations about $y$-axis and four series
of simultaneous CZ gates (three gates per series). The syndromes for
Surface-17 are

$\displaystyle S_{A0}$ $\displaystyle=X_{1}X_{2}$ (24) $\displaystyle S_{A1}$
$\displaystyle=Z_{0}Z_{3}$ $\displaystyle S_{A2}$
$\displaystyle=X_{0}X_{1}X_{3}X_{4}$ $\displaystyle S_{A3}$
$\displaystyle=Z_{1}Z_{2}Z_{4}Z_{5}$ $\displaystyle S_{A4}$
$\displaystyle=Z_{3}Z_{4}Z_{6}Z_{7}$ $\displaystyle S_{A5}$
$\displaystyle=X_{4}X_{5}X_{7}X_{8}$ $\displaystyle S_{A6}$
$\displaystyle=Z_{5}Z_{8}$ $\displaystyle S_{A7}$ $\displaystyle=X_{6}X_{7},$

where the indices of Pauli $X$ and $Z$ operations relate to the number of data
qubit.

It is remarkable that $X$-ancilla measurement starts right after the
$R_{y}^{\pi/2}$ gate on qubits {A0, A2, A5, A7}, and $Z$-ancilla measurement
continues during the first four series of CZ gates as it is shown in the
timeline in Fig 8. Thereby, the full correction cycle time consists of 4 two-
qubit gates, 2 single-qubit rotations, 1 measurement and reset procedure,
totaling at 800 ns.

The circuit is simulated under the same noise models as the five-qubit code.


# Conformal Thresholded Intervals for Efficient Regression

###### Abstract.

This paper introduces Conformal Thresholded Intervals (CTI), a novel conformal
regression method that aims to produce the smallest possible prediction set
with guaranteed coverage. Unlike existing methods that rely on nested
conformal framework and full conditional distribution estimation, CTI
estimates the conditional probability density for a new response to fall into
each interquantile interval using off-the-shelf multi-output quantile
regression. CTI constructs prediction sets by thresholding the estimated
conditional interquantile intervals based on their length, which is inversely
proportional to the estimated probability density. The threshold is determined
using a calibration set to ensure marginal coverage. Experimental results
demonstrate that CTI achieves optimal performance across various datasets.

Rui Luo and Zhixin Zhou

##  1\. Introduction

Conformal prediction is a powerful framework for constructing prediction
intervals with finite-sample validity guarantees. By leveraging
exchangeability of the data, conformal methods can convert the output of any
machine learning algorithm into a set-valued prediction satisfying the
required coverage level, without assumptions on the data distribution. This
paper develops a novel conformal prediction method for regression that aims to
produce the smallest possible prediction set with guaranteed coverage.

Most existing conformal methods for regression either directly predict the
lower and upper endpoints of the interval using quantile regression models
[30, 17, 32, 11] or first estimate the full conditional distribution of the
response and then invert it to obtain prediction sets [13, 4]. While these
approaches perform well in many situations, they may produce sub-optimal
intervals if the conditional distribution is skewed. Conformalized quantile
regression typically yields equal-tailed intervals, but the shortest valid
interval may be unbalanced. On the other hand, density-based methods can adapt
to skewness but typically involve many tuning parameters and more difficult
interpretation, which can be complex for practitioners.

To address these limitations, we propose conformal thresholded intervals
(CTI), a conformal inference method that seeks the smallest possible
prediction set. Instead of relying on an estimate of the full conditional
distribution, we use off-the-shelf multi-output quantile regression and
construct prediction set by thresholding the estimated conditional
interquantile intervals. Compared with conformal histogram regression (CHR)
[33], which first partitions the response space into bins, CTI directly trains
a multi-output quantile regression model that uses equiprobable quantiles.
This allows us to estimate the conditional probability density for a new
response to fall into each interquantile interval, without the need for
explicitly binning the response space.

For each sample in the calibration set, we obtain the interquantile interval
that its response falls into and find the corresponding probability density
estimate. We compute the non-conformity scores based on these estimates.
Intuitively, the non-conformity score is higher for a sample that falls into a
long interquantile interval and lower for a sample that falls into a short
interquantile interval. By adopting a similar thresholding idea as in
conformal classification [31, 23], we threshold the intervals according to
their length, the inverse of which corresponds to the probability density
estimate. At test time, the threshold, i.e., the quantile for non-conformity
scores, is used in constructing prediction sets for test samples.
Specifically, the interquantile intervals are sorted in ascending order of
length, and the first ones shorter than or equal to the threshold are kept. We
show that the prediction sets generated from thresholding interquantile
intervals guarantee marginal coverage and can achieve desired conditional
coverage as well as the smallest expected prediction interval length if the
multi-output quantile regression model produces true conditional probability
density estimates.

The rest of this paper is organized as follows. We discuss related work in
Section 2. Section 3 describes the proposed CTI method in detail. Section 4
presents a theoretical analysis of CTI. Section 5 provides numerical
experiments comparing CTI to existing conformal regression methods on both
simulated and real data. Finally, Section 6 concludes with a discussion of the
main results and future directions.

##  2\. Related Work

Quantile regression, introduced by [18], estimates the $\tau$-th conditional
quantile function by minimizing the check function loss:

(1) $\min_{f_{\tau}}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-q_{\tau}(x_{i})),$

where

(2) $\rho_{\tau}(r)=\begin{cases}\tau r&\textrm{if $r>0$}\\\
-(1-\tau)r&\textrm{otherwise}\end{cases}$

is the check function representing the absolute loss.

Quantile regression has been widely used to construct prediction intervals by
estimating conditional quantile functions at specific levels, such as the 5%
and 95% levels for 90% nominal coverage [12, 36, 25, 35, 34]. This approach
adapts to local variability, even for highly heteroscedastic data.

[5] showed that simultaneous estimation of multiple quantiles is
asymptotically more efficient than separate estimation of individual
regression quantiles or ignoring within-subject dependency. However, this
approach does not guarantee non-crossing quantiles, which can affect the
validity of the predictions and introduce critical issues in certain
scenarios. To address this limitation, research on non-crossing multiple
quantile regression has gained attention in recent years, with several methods
proposed to ensure non-crossing quantile estimates, including stepwise
approaches [21], non-parametric techniques [3], and deep learning-based models
[27, 2].

However, the validity of the produced intervals is only guaranteed for
specific models under certain regularity and asymptotic conditions [34, 35,
25]. Many related methods for constructing valid prediction intervals can be
encompassed within the nested conformal prediction framework, where a nested
sequence of prediction sets is generated by thresholding nonconformity scores
derived from various approaches, such as residual-based methods [29, 1, 19],
quantile regression [30, 17, 32, 4], density estimation [14, 33, 15], and
their combinations with ensemble methods [11] and localized methods [28, 6,
22]. However, as noted by [20], the optimal conditionally-valid prediction
regions are level sets of conditional densities, which need not be intervals,
suggesting that constructing possibly non-convex prediction sets might lead to
more efficient conformal predictors.

Our proposed method for constructing non-convex prediction sets is related to
the work of [15], who introduce a profile distance to measure the similarity
between features and construct prediction sets based on neighboring samples.
In contrast, our method directly estimates the conditional probability density
for a new response to fall into each interquantile interval based on a multi-
output quantile regression model. By thresholding the interquantile intervals
based on their length, which is inversely proportional to the estimated
probability density, we can construct efficient prediction sets that adapt to
the local density of the data. This approach allows us to generate prediction
sets that are not restricted to intervals and can potentially achieve better
coverage and efficiency compared to interval-based methods.

Another related approach [10] converts regression to a classification problem
and employs a conditional distribution with a smoothness-enforcing penalty.
This method is orthogonal to our approach and can be potentially combined with
our multi-output quantile regression framework to further improve the
efficiency of the constructed prediction sets.

##  3\. Proposed Method

###  3.1. Problem Setup

We consider a general regression problem with a dataset
$\\{(x_{i},y_{i})\\}_{i=1}^{n}$, where
$x_{i}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ is the input feature vector and
$y_{i}\in\mathcal{Y}\subseteq\mathbb{R}$ is the corresponding continuous
response variable. The dataset is split into three parts: a training set
$\mathcal{D}_{\text{train}}$, a calibration set $\mathcal{D}_{\text{cal}}$,
and a test set $\mathcal{D}_{\text{test}}$. The corresponding indices set are
denoted by $\mathcal{I}_{\textrm{train}},\mathcal{I}_{\textrm{cal}}$ and
$\mathcal{I}_{\textrm{test}}$ respectively. We assume that the examples in
these sets are exchangeable.

Our goal is to construct a conformal predictor that outputs a prediction set
$\mathcal{C}(X)\subseteq\mathcal{Y}$ for each test input $X$ such that the
true response value $Y$ is included in $\mathcal{C}(X)$ with a probability of
at least $1-\alpha$, where $\alpha\in(0,1)$ is a user-specified significance
level. Formally, we aim to achieve the following marginal coverage guarantee:

$\mathbb{P}(Y\in\mathcal{C}(X))\geq 1-\alpha$

for joint distribution for $X$ and $Y$. The probability in this statement is
marginal, being taken over all the samples in $\mathcal{D}_{\text{cal}}$ and
$\mathcal{D}_{\text{test}}$.

While achieving valid coverage, we aim to construct prediction sets that are
as informative as possible. Specifically, we seek to minimize the expected
length of the prediction sets:

$\mathbb{E}[\mu(\mathcal{C}(X))]=\int_{\mathcal{X}}\mu(\mathcal{C}(x))dP(x),$

where $\mu$ denotes the Lebesgue measure on $\mathbb{R}$ and $P(x)$ is the
marginal distribution of the input features.

###  3.2. Our Method

First, we apply quantile regression on the training set
$\mathcal{D}_{\text{train}}$ to predict the $\tau$-th quantile of the
conditional distribution $Y|X=x$ for every $x\in\mathcal{X}$, where $\tau$
takes values from 0 to 1 in increments of $1/K$. The estimated quantile for
$\tau=k/K$ is denoted by

(3) $\displaystyle\widehat{q}_{k}(x)\quad\text{for}\quad k=0,1,\dots,K.$

We then define the interquantile intervals as

(4) $\displaystyle I_{k}(x)=(\widehat{q}_{k-1}(x),\
\widehat{q}_{k}(x)]\quad\text{for}\quad k=\ 1,\dots,K.$

Assuming the quantile regression provides sufficiently accurate estimations,
each interval should have approximately the same probability, $1/K$, of
covering the true label $Y$. To minimize the size of the prediction set, it is
more efficient to include intervals with smaller sizes. This strategy leads us
to define the confidence set as:

(5) $\displaystyle\mathcal{C}(x)=\bigcup\\{I_{k}(x):\mu(I_{k}(x))\leq
t,k=1,\dots,K\\},$

where $t$ is a threshold determined in a marginal sense, meaning it is
independent of $x$. To determine $t$, we utilize the calibration set. We want
$t$ to satisfy the condition that $y_{i}\in\mathcal{C}(x_{i})$ for at least
$\lceil(1+|\mathcal{I}_{\text{cal}}|)(1-\alpha)\rceil$ instances in the
calibration set, where $i\in\mathcal{I}_{\text{cal}}$. We define $t$ as the
smallest value that satisfies this condition:

(6) $\displaystyle
t=\lceil(1+|\mathcal{I}_{\text{cal}}|)(1-\alpha)\rceil\text{-th smallest value
of }\mu(I_{k(y_{i})}(x_{i}))\text{ for }i\in\mathcal{I}_{\textrm{cal}},$

where $k(y)$ is the index that of the interval that $y$ belongs, i.e., $y\in
I_{k(y)}(x)$. By plugging $t$ back into (5), we obtain the prediction set for
every $x\in\mathcal{X}$.

Algorithm 1 Conformal thresholding intervals

1:Input: labeled data $\\{(x_{i},y_{i})\\}_{i\in\mathcal{I}}$, unlabel test
data $\\{x_{i}\\}_{\mathcal{I}_{\textrm{test}}}$, a data split ratio, black-
box learning algorithm $\mathcal{B}$, level $\alpha\in(0,1)$, number of
interquantile intervals $K$

2:Randomly split the indices $\mathcal{I}$ into $\mathcal{I}_{\textrm{train}}$
and $\mathcal{I}_{\textrm{cal}}$.

3:Train $\mathcal{B}$ on samples in $\mathcal{I}_{\textrm{train}}$, and obtain
quantile estimation functions $\widehat{q}_{k}$ for $k=0,1,\dots,K$.

4:For every $i\in\mathcal{I}_{\textrm{cal}}\cup\mathcal{I}_{\textrm{test}}$,
evaluate $\widehat{q}_{k}(x_{i})$ for $k=0,1,\dots,K$.

5:For every $i\in\mathcal{I}_{\textrm{cal}}\cup\mathcal{I}_{\textrm{test}}$,
define the interquantile intervals
$I_{k}(x_{i})=(\widehat{q}_{k-1}(x_{i}),\widehat{q}_{k}(x_{i})]$ for
$k=1,\dots,K$.

6:$t\leftarrow\lceil(1+|\mathcal{I}_{\text{cal}}|)(1-\alpha)\rceil\text{-th
smallest value of }\mu(I_{k(y_{i})}(x_{i}))\text{ for
}i\in\mathcal{I}_{\textrm{cal}}$.

7:For
$i\in\mathcal{I}_{\textrm{test}},\mathcal{C}(x_{i})=\bigcup\\{I_{k}(x_{i}):\mu(I_{k}(x_{i}))\leq
t,k=1,\dots,K\\}.$

8:Output: $\mathcal{C}(x_{i})$ for $i\in\mathcal{I}_{\textrm{test}}$.

######  Remark 3.3.

The procedure we propose is generic and can be applied to any multi-output
quantile regression method. In Section 5, we present results based on quantile
regression methods using both a neural network and a random forest. In some
cases, quantile regression algorithms may not support $\tau$ values of exactly
0 or 1. To overcome this limitation, we can replace these extreme values with
values very close to 0 or 1, respectively. For example, we might use
$\tau=0.001$ instead of $\tau=0$, and $\tau=0.999$ instead of $\tau=1$. For
the sake of notational simplicity, we will continue to refer to these values
as 0 and 1 throughout the discussion, unless otherwise specified.

######  Remark 3.4.

Our approach can also be considered in terms of conformity scores. Using the
definition of the prediction set in equation (5), the value of label $y$ is
contained within a small interval. More formally, let $k(y)$ be the index such
that $y\in I_{k(y)}(x)$. We can then define the conformity score function for
our proposed method as:

$\displaystyle s(x,y)=\mu(I_{k(y)}(x)).$

This conformity score function assigns a score to each label $y$ based on the
size of the interval $I_{k(y)}(x)$ in which it falls. A smaller score
indicates that the label $y$ is more likely to be the true label for the input
$x$. In the context of conformal prediction, labels with smaller conformity
scores are given priority for inclusion in the prediction set.

##  4\. Theoretical Analysis

In the context of the entire population, CTI shares a very similar formulation
with the Least Ambiguous Set method used for classification, as described in
[31]. If we assume that our quantile regression model is sufficiently
accurate, CTI has the potential to achieve the optimal size for prediction
sets when considering the marginal distribution. To understand this better,
let’s first take a look at the Neyman-Pearson Lemma:

######  Lemma 4.1 (Neyman-Pearson).

Let $f$ and $g$ be two nonnegative measurable functions. Then the optimizer of
the problem

$\displaystyle\min_{C}\int_{C}g\quad\text{subject to}\quad\int_{C}f\geq
1-\alpha,$

is given by $C=\\{x:f(x)/g(x)\geq t^{\prime}\\}$ if there exists $t$ such that
$\int_{f/g\geq t^{\prime}}f=1-\alpha$.

To formalize the problem of minimizing the expected length of the prediction
set subject to $1-\alpha$ coverage, we can write the problem as:

$\displaystyle\min_{\mathcal{C}(x)}\int_{\mathcal{X}}\int_{\mathcal{C}(x)}1d\mu(y)dP(x)\quad\text{subject
to}\quad\int_{\mathcal{X}}\int_{\mathcal{C}(x)}f(y|x)d\mu(y)dP(x)\geq
1-\alpha.$

The Neyman-Pearson Lemma implies that the optimal solution for
$\mathcal{C}(x)$ has the form:

(7) $\displaystyle\mathcal{C}(x)=\\{y:f(y|x)\geq t^{\prime}\\}$

for some suitable threshold $t^{\prime}$. Indeed, this threshold can be
defined as

(8) $\displaystyle t^{\prime}=\inf\\{t\in\mathbb{R}:\mathbb{P}(f(Y|X)\geq
t)\geq 1-\alpha\\}.$

which will be shown in Lemma A.2. Our algorithm is an empirical construction
of such an interval. Suppose the quantile regression approximates
$\widehat{q}_{\tau}$ well. In that case, we have:

$\displaystyle\int_{y\in I_{k}(x)}f(y|x)d\mu(y)=\mathbb{P}(Y\in
I_{k}(X))\approx 1/K.$

As $K$ approaches infinity, $\mu(I_{k}(x))$ tends to 0. If $f(y|x)$ is
sufficiently smooth, then $f(y|x)\approx 1/(K\mu(I_{k}(x)))$. The threshold on
the length of intervals $\mu(I_{k}(x))\leq t$ in equation (5) approximately
implies $f(y|x)\geq 1/(Kt)$, which is optimal in the sense of the Neyman-
Pearson Lemma. This means that our algorithm, which constructs prediction sets
based on the threshold on interval lengths, is an empirical approximation of
the optimal solution prescribed by the Neyman-Pearson Lemma.

The demonstration of the coverage probability for CTI follows the same
reasoning as the traditional proof used in the general conformal prediction
framework.

######  Theorem 4.2 (Coverage Probability).

Suppose the samples in
$\\{(X_{i},Y_{i})\\}_{i\in\mathcal{I}_{\textrm{cal}}\cup\mathcal{I}_{\textrm{test}}}\\}$
are exchangeable, then for $(X,Y)$ in the test set, the coverage probability

$\displaystyle\mathbb{P}\left(Y\in\mathcal{C}(X)\right)\geq 1-\alpha.$

###### Proof.

The score function $s(X_{i},Y_{i})$ for
$i\in\mathcal{I}_{\textrm{cal}}\cup\mathcal{I}_{\textrm{test}}$ are also
exchangeable. For any $(X,Y)$ in the test set, the rank of $s(X,Y)$ is smaller
than $t$ defined in (6) with probability
$\frac{\lceil(1+|\mathcal{I}_{\textrm{cal}}|)(1-\alpha)\rceil}{1+|\mathcal{I}_{\textrm{cal}}|}\geq
1-\alpha$. ∎

The upcoming proposition will demonstrate that the threshold $t$, which is
defined in equation (6) for the length of interquantile intervals, results in
a suitable threshold for the distribution of the conditional density $f(Y|X)$,
as shown in equation (8).

######  Proposition 1 (Threshold Consistency).

Suppose the interquantile intervals satisfy

(9)
$\displaystyle\sup_{t^{\prime}}\left|F_{f(Y|X)}(t^{\prime})-\frac{1}{\mathcal{I}_{\textrm{cal}}}\sum_{i\in\mathcal{I}_{\textrm{cal}}}\mathbf{1}\big{\\{}\mu(I_{k(y_{i})}(x_{i}))>\frac{1}{Kt^{\prime}}\big{\\}}\right|\leq\epsilon$

Let $t$ be as defined in (6) and let
$1-\alpha^{\prime}=\lceil(1+|\mathcal{I}_{\text{cal}}|)(1-\alpha)\rceil/|\mathcal{I}_{\textrm{cal}}|$,
then $F_{f(Y|X)}(1/(Kt))\geq\alpha^{\prime}-\epsilon$.

###### Proof.

Let $t^{\prime}=1/(Kt)$ in the assumption (9), then by the definition of $t$,

$\displaystyle\sum_{i\in\mathcal{I}_{\textrm{cal}}}\mathbf{1}\big{\\{}\mu(I_{k(y_{i})}(x_{i}))>\frac{1}{Kt^{\prime}}\big{\\}}=|\mathcal{I}_{\textrm{cal}}|-\lceil(1+|\mathcal{I}_{\text{cal}}|)(1-\alpha)\rceil=:|\mathcal{I}_{\text{cal}}|\alpha^{\prime}.$

Under the assumption, $F_{f(Y|X)}(t^{\prime})\geq\alpha^{\prime}-\epsilon$. ∎

This theorem establishes that if the value of $\epsilon$ is sufficiently small
and the size of the calibration set is large enough, then the value
$t^{\prime}=1/(Kt)$, where $t$ is the threshold defined in equation (6), will
be close to the $\alpha$-th quantile of the distribution $f(Y|X)$.

Lastly, the theorem presented below demonstrates that the size of the
prediction set obtained using CTI will not significantly exceed the size of
the theoretically optimal prediction set, which is defined in equation (7).

######  Theorem 4.3 (Prediction Set).

Suppose for $x\in\mathcal{X}$,

(10) $\displaystyle P(Y\in I_{k}(x)|X=x)=\int_{y\in
I_{k}(x)}f(y|x)dy\geq\frac{1-\delta(x)}{K},$

and suppose $f(y|x)$ has a Lipschitz constant $L(x)$, then

$\displaystyle\mathcal{C}_{1-\alpha}(x)\subseteq\left\\{y:f(y|x)\geq\frac{1-\delta_{k}(x)}{Kt}-\frac{L(x)t}{2}\right\\}.$

###### Proof.

In the construction of $\mathcal{C}_{1-\alpha}(x)$, we consider the intervals
with length$\mu(I_{k}(x))\leq t$. Given the assumption (10) in the theorem,
and combine it with Lemma A.1, we have

$\displaystyle\min_{y\in
I_{k}}f(y|x)\geq\frac{1-\delta_{k}(x)}{Kt}-\frac{L(x)t}{2}=:t^{\prime}.$

By the construction of $\mathcal{C}_{1-\alpha}(x)$, we have the conclusion of
the theorem. ∎

The previous proposition demonstrates that the value $1/(Kt)$, where $t$ is
the threshold defined in equation (6), is slightly smaller than the
$\alpha$-th quantile of $f(Y|X)$, which represents the theoretically optimal
threshold as shown in equation (8). This theorem further illustrates that the
prediction set obtained using CTI will include values that are even more
conservative. However, we understand that $1/Kt$, where $t$ is the threshold
defined in equation (6), serves as a relatively stable estimate of the
quantile, which is asymptotically equivalent to a constant value. As the
number of interquantile intervals $K$ approaches infinity, the threshold $t$
converges to 0. If we additionally assume that our quantile regression model
is accurate, meaning that the error term $\delta_{k}(x)$ is small, then the
prediction set $\mathcal{C}_{1-\alpha}(x)$ obtained using CTI will be close to
the theoretically optimal prediction set defined in equation (5).

##  5\. Experiments

This section presents experiments that evaluate the performance of prediction
sets generated by Split Conformal [19], CQR [30], CHR [33], and our proposed
CTI, on various simulated and real-world datasets as in [30, 17, 32]. For a
detailed introduction of the datasets, please refer to [32]. We have released
our code on GitHub, which is built upon TorchCP [37], CHR [33], and QRF [16].

Following similar procedures, we standardize the response variables $Y$ for
all datasets. We randomly hold out 20% of the samples for testing and divide
the remaining observations into two disjoint sets, $\mathcal{I}_{1}$ and
$\mathcal{I}_{2}$, for training the multi-output quantile regression model and
computing the non-conformity scores, respectively. We use 70% of the samples
for training, which achieves a reasonable compromise for all datasets [32]. We
repeat all experiments 10 times, starting from the initial data splitting. We
evaluate the performance of the generated prediction intervals in terms of
coverage and efficiency.

As shown in Table 1, CTI achieves optimal performance across all datasets
except for the star dataset, which has a relatively small sample size
($n=2161,d=39$). The limited number of samples in this dataset may hinder the
performance of the multi-output quantile regression model, as it requires
sufficient data to accurately capture the underlying relationships between the
features and the response variable. We also notice a similar trend in the
relative performance comparison of CTI based on random forest and CTI based on
neural network, as well as CHR based on random forest and CHR based on neural
network. This suggests that the efficiency of the conformal prediction sets
depends on the quality of the multi-output quantile regression. The choice of
the underlying model plays a crucial role in the performance of the conformal
prediction methods.

To further analyze the superior performance of CTI, we compare the interval
lengths between the response intervals (i.e., intervals containing the actual
responses) and all intervals generated by the multi-output quantile regression
model for the MEPS-19 dataset. Figure 1 presents the distribution of interval
lengths for both the response intervals (blue histogram) and all intervals
(red histogram) on the test set. The difference in means between the two
distributions is -0.0415, indicating that the response intervals have a
smaller average length compared to all intervals.

This finding suggests that the multi-output quantile regression model
effectively captures the underlying structure of the data, producing tighter
intervals around the true responses. By estimating the conditional quantile
functions simultaneously, the model can leverage the dependencies between
different quantile levels and generate more precise intervals. The multi-
output approach allows the model to share information across different
quantile levels, leading to improved interval estimates.

Consequently, CTI can construct more efficient prediction sets by thresholding
the intervals based on their estimated probability density. By selecting
intervals with higher probability density, CTI focuses on regions of the
feature space where the true response is more likely to occur. This adaptive
thresholding approach leads to improved coverage and shorter interval lengths
compared to other methods that rely on fixed thresholds or do not consider the
local density of the data.

Dataset Metric CTI(RF) CTI(NN) CHR(RF) CHR(NN) CQR Split bike Coverage 0.898
(0.007) 0.899 (0.007) 0.898 (0.010) 0.900 (0.006) 0.906 (0.009) 0.899 (0.008)
Size 1.032 (0.029) 0.720 (0.028) 1.124 (0.028) 0.758 (0.047) 1.599 (0.054)
1.345 (0.053) blog Coverage 0.910 (0.002) 0.900 (0.004) 0.902 (0.004) 0.902
(0.003) 0.940 (0.009) 0.910 (0.006) Size 0.709 (0.031) 1.003 (0.024) 1.567
(0.074) 1.737 (0.154) 3.259 (0.327) 1.453 (0.113) community Coverage 0.909
(0.018) 0.908 (0.021) 0.903 (0.015) 0.905 (0.021) 0.889 (0.024) 0.902 (0.024)
Size 1.611 (0.088) 1.275 (0.095) 1.637 (0.096) 1.588 (0.100) 1.680 (0.078)
2.132 (0.188) star Coverage 0.903 (0.018) 0.910 (0.017) 0.907 (0.021) 0.897
(0.018) 0.901 (0.016) 0.910 (0.024) Size 0.186 (0.006) 0.197 (0.009) 0.182
(0.005) 0.204 (0.009) 0.181 (0.005) 0.181 (0.008) homes Coverage 0.900 (0.005)
0.900 (0.006) 0.899 (0.005) 0.895 (0.007) 0.898 (0.006) 0.897 (0.005) Size
0.640 (0.011) 0.515 (0.008) 0.682 (0.012) 0.535 (0.010) 0.851 (0.052) 0.825
(0.072) facebook1 Coverage 0.909 (0.003) 0.899 (0.003) 0.901 (0.004) 0.900
(0.004) 0.945 (0.009) 0.903 (0.002) Size 0.766 (0.033) 0.780 (0.023) 1.595
(0.088) 1.379 (0.086) 2.627 (0.329) 2.252 (0.208) facebook2 Coverage 0.911
(0.002) 0.900 (0.001) 0.899 (0.002) 0.899 (0.002) 0.943 (0.006) 0.904 (0.002)
Size 0.735 (0.017) 0.773 (0.023) 1.533 (0.053) 1.382 (0.057) 2.661 (0.272)
2.100 (0.108) bio Coverage 0.900 (0.004) 0.902 (0.004) 0.899 (0.005) 0.900
(0.004) 0.900 (0.003) 0.901 (0.004) Size 1.295 (0.018) 1.474 (0.030) 1.450
(0.023) 1.576 (0.012) 2.005 (0.016) 1.961 (0.039) concrete Coverage 0.908
(0.024) 0.900 (0.031) 0.899 (0.022) 0.900 (0.023) 0.901 (0.024) 0.896 (0.021)
Size 0.967 (0.035) 0.473 (0.050) 0.933 (0.041) 0.505 (0.144) 0.692 (0.051)
0.619 (0.029) meps19 Coverage 0.907 (0.008) 0.902 (0.007) 0.901 (0.007) 0.902
(0.004) 0.932 (0.007) 0.902 (0.010) Size 1.760 (0.087) 1.795 (0.061) 2.388
(0.195) 2.602 (0.128) 2.923 (0.170) 3.092 (0.377) meps20 Coverage 0.904
(0.004) 0.901 (0.007) 0.901 (0.007) 0.901 (0.006) 0.927 (0.009) 0.902 (0.005)
Size 1.883 (0.067) 1.921 (0.091) 2.376 (0.105) 2.594 (0.140) 2.925 (0.193)
3.154 (0.217) meps21 Coverage 0.906 (0.005) 0.900 (0.008) 0.900 (0.006) 0.898
(0.004) 0.928 (0.007) 0.905 (0.004) Size 1.832 (0.089) 1.866 (0.076) 2.510
(0.167) 2.609 (0.145) 2.971 (0.179) 3.046 (0.199)

Table 1. Coverage and Size of Different Methods Figure 1. Comparison of
interval lengths between response intervals and all intervals for the MEPS-19
dataset. The blue histogram represents the distribution of interval lengths
for the intervals that contain the actual responses, while the red histogram
shows the distribution of all intervals generated by the multi-output quantile
regression model. Both all results on the test set. The difference in means
between the two distributions is -0.0415, indicating that the response
intervals have a smaller average length compared to all intervals.

##  6\. Discussion and Future Work

Our proposed conformal prediction method for regression tasks has demonstrated
promising results in achieving high performance. However, there are several
areas where further research and improvements can be made.

Computing prediction intervals as opposed to prediction sets. There might be
concerns regarding constructing disjoint prediction sets instead of prediction
intervals, as it may be more desirable to have a single contiguous interval
for interpretability and practicality. To address this issue, we can adopt
ideas similar to those in [7], which deals with the problem of non-contiguous
prediction sets for ordinal classification tasks. The authors propose training
the predictor with a unimodal posterior over classes, allowing the
construction of contiguous prediction sets. In the context of regression, we
can impose unimodality constraints on the estimated conditional density or
quantile functions to ensure that the resulting prediction sets are contiguous
intervals. This approach would provide a more interpretable and user-friendly
output while maintaining the advantages of our proposed method.

Weighted aggregation of score functions. In this work, we implement the
concept of weighted aggregation of score functions for conformal
classification [24] by applying it to multi-output quantile regression models.
Specifically, we combine the results from a neural network model and a random
forest model. To ensure the combined results form a valid non-conformity score
and accurately reflect predictions from both models, we use a harmonic mean
for aggregation. This approach allows us to leverage the strengths of
different models and improve the overall performance of the conformal
prediction sets. Future research can explore alternative aggregation methods
and investigate the optimal combination of models for different datasets and
problem settings.

Extension to other regression settings. Our current work focuses on the
standard regression setting with continuous response variables. However, there
are various other regression settings, such as multivariate regression or
functional regression [8, 26, 9], where conformal prediction methods could be
beneficial. Adapting our multi-output quantile regression approach to these
settings would require addressing the specific challenges and characteristics
of each problem. For example, in multivariate regression, the construction of
prediction sets would involve estimating joint quantile functions and handling
the dependencies between multiple response variables. In functional
regression, the response variable is a function rather than a scalar,
requiring the development of appropriate non-conformity measures and
prediction set construction techniques.

##  7\. Conclusion

Conformal Thresholded Intervals (CTI) is a novel conformal prediction method
for regression that constructs the smallest possible prediction set with
guaranteed coverage. By leveraging multi-output quantile regression and
thresholding the estimated conditional interquantile intervals based on their
length, CTI adapts to the local density of the data and generates non-convex
prediction sets. Experimental results demonstrate the effectiveness of CTI in
achieving optimal performance across various datasets. However, the
performance of CTI depends on the quality of the underlying multi-output
quantile regression model, emphasizing the importance of model selection in
conformal prediction methods. Future research directions include optimizing
model selection as well as model aggregation and extending CTI to other types
of regression problems.

## References

  * [1] Balasubramanian, V., Ho, S.-S., and Vovk, V. Conformal prediction for reliable machine learning: theory, adaptations and applications.  Newnes, 2014.
  * [2] Brando, A., Center, B. S., Rodríguez-Serrano, J., Vitrià, J., et al. Deep non-crossing quantiles through the partial derivative.  In International Conference on Artificial Intelligence and Statistics (2022), PMLR, pp. 7902–7914.
  * [3] Cannon, A. J. Non-crossing nonlinear regression quantiles by monotone composite quantile regression neural network, with application to rainfall extremes.  Stochastic environmental research and risk assessment 32, 11 (2018), 3207–3225.
  * [4] Chernozhukov, V., Wüthrich, K., and Zhu, Y. Distributional conformal prediction.  Proceedings of the National Academy of Sciences 118, 48 (2021), e2107794118.
  * [5] Cho, H., Kim, S., and Kim, M.-O. Multiple quantile regression analysis of longitudinal data: Heteroscedasticity and efficient estimation.  Journal of Multivariate Analysis 155 (2017), 334–343.
  * [6] Colombo, N. On training locally adaptive cp.  In Conformal and Probabilistic Prediction with Applications (2023), PMLR, pp. 384–398.
  * [7] DEY, P., Merugu, S., and Kaveri, S. R. Conformal prediction sets for ordinal classification.  In Thirty-seventh Conference on Neural Information Processing Systems (2023).
  * [8] Diquigiovanni, J., Fontana, M., and Vantini, S. Conformal prediction bands for multivariate functional data.  Journal of Multivariate Analysis 189 (2022), 104879.
  * [9] Feldman, S., Bates, S., and Romano, Y. Calibrated multiple-output quantile regression with representation learning.  Journal of Machine Learning Research 24, 24 (2023), 1–48.
  * [10] Guha, E. K., Natarajan, S., Möllenhoff, T., Khan, M. E., and Ndiaye, E. Conformal prediction via regression-as-classification.  In The Twelfth International Conference on Learning Representations (2024).
  * [11] Gupta, C., Kuchibhotla, A. K., and Ramdas, A. Nested conformal prediction and quantile out-of-bag ensemble methods.  Pattern Recognition 127 (2022), 108496.
  * [12] Hunter, D. R., and Lange, K. Quantile regression via an mm algorithm.  Journal of Computational and Graphical Statistics 9, 1 (2000), 60–77.
  * [13] Izbicki, R., Shimizu, G., and Stern, R. Flexible distribution-free conditional predictive bands using density estimators.  In International Conference on Artificial Intelligence and Statistics (2020), PMLR, pp. 3068–3077.
  * [14] Izbicki, R., Shimizu, G., and Stern, R. Flexible distribution-free conditional predictive bands using density estimators.  In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics (26–28 Aug 2020), S. Chiappa and R. Calandra, Eds., vol. 108 of Proceedings of Machine Learning Research, PMLR, pp. 3068–3077.
  * [15] Izbicki, R., Shimizu, G., and Stern, R. B. Cd-split and hpd-split: Efficient conformal regions in high dimensions.  Journal of Machine Learning Research 23, 87 (2022), 1–32.
  * [16] Johnson, R. A. quantile-forest: A python package for quantile regression forests.  Journal of Open Source Software 9, 93 (2024), 5976.
  * [17] Kivaranovic, D., Johnson, K. D., and Leeb, H. Adaptive, distribution-free prediction intervals for deep networks.  In International Conference on Artificial Intelligence and Statistics (2020), PMLR, pp. 4346–4356.
  * [18] Koenker, R. Quantile regression, vol. 38.  Cambridge university press, 2005.
  * [19] Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. Distribution-free predictive inference for regression.  Journal of the American Statistical Association 113, 523 (2018), 1094–1111.
  * [20] Lei, J., Robins, J., and Wasserman, L. Distribution-free prediction sets.  Journal of the American Statistical Association 108, 501 (2013), 278–287.
  * [21] Liu, Y., and Wu, Y. Stepwise multiple quantile regression estimation using non-crossing constraints.  Statistics and its Interface 2, 3 (2009), 299–310.
  * [22] Luo, R., and Colombo, N. Conformal load prediction with transductive graph autoencoders.  arXiv preprint arXiv:2406.08281 (2024).
  * [23] Luo, R., and Zhou, Z. Trustworthy classification through rank-based conformal prediction sets.  arXiv preprint arXiv:2407.04407 (2024).
  * [24] Luo, R., and Zhou, Z. Weighted aggregation of conformity scores for classification, 2024.
  * [25] Meinshausen, N. Quantile regression forests.  Journal of Machine Learning Research 7 (2006), 983–999.
  * [26] Messoudi, S., Destercke, S., and Rousseau, S. Ellipsoidal conformal inference for multi-target regression.  In Conformal and Probabilistic Prediction with Applications (2022), PMLR, pp. 294–306.
  * [27] Moon, S. J., Jeon, J.-J., Lee, J. S. H., and Kim, Y. Learning multiple quantiles with neural networks.  Journal of Computational and Graphical Statistics 30, 4 (2021), 1238–1248.
  * [28] Papadopoulos, H., Gammerman, A., and Vovk, V. Normalized nonconformity measures for regression conformal prediction.  In Proceedings of the IASTED International Conference on Artificial Intelligence and Applications (AIA 2008) (2008), pp. 64–69.
  * [29] Papadopoulos, H., Proedrou, K., Vovk, V., and Gammerman, A. Inductive confidence machines for regression.  In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19–23, 2002 Proceedings 13 (2002), Springer, pp. 345–356.
  * [30] Romano, Y., Patterson, E., and Candès, E. Conformalized quantile regression.  In Advances in Neural Information Processing Systems (2019), pp. 3538–3548.
  * [31] Sadinle, M., Lei, J., and Wasserman, L. Least ambiguous set-valued classifiers with bounded error levels.  Journal of the American Statistical Association 114, 525 (2019), 223–234.
  * [32] Sesia, M., and Candès, E. A comparison of some conformal quantile regression methods.  Stat 9, 1 (2020).
  * [33] Sesia, M., and Romano, Y. Conformal prediction using conditional histograms.  Advances in Neural Information Processing Systems 34 (2021), 6304–6315.
  * [34] Steinwart, I., and Christmann, A. Estimating conditional quantiles with the help of the pinball loss.
  * [35] Takeuchi, I., Le, Q. V., Sears, T. D., Smola, A. J., and Williams, C. Nonparametric quantile estimation.  Journal of machine learning research 7, 7 (2006).
  * [36] Taylor, J. W. A quantile regression neural network approach to estimating the conditional density of multiperiod returns.  Journal of Forecasting 19, 4 (2000), 299–311.
  * [37] Wei, H., and Huang, J. Torchcp: A library for conformal prediction based on pytorch.  arXiv preprint arXiv:2402.12683 (2024).

##  Appendix A Ancillary Lemmas

######  Lemma A.1.

For a Lipschitz function $f:\mathbb{R}\to\mathbb{R}$ with Lipschitz constant
$L$, if $\int_{a}^{b}f(x)d\mu(x)=c$, then
$f(x)\in[\frac{c}{b-a}-\frac{L(b-a)}{2},\frac{c}{b-a}+\frac{L(b-a)}{2}]$ for
all $x\in[a,b]$.

###### Proof.

Suppose $f(x^{\prime})=d$ for some $x^{\prime}\in[a,b]$. Then, by the
Lipschitz condition, we have:

$\displaystyle|f(x)-f(x^{\prime})|\leq L|x-x^{\prime}|\implies f(x)\leq
d+L|x-x^{\prime}|$

for all $x\in[a,b]$. Integrating both sides over $[a,b]$, we get:

$\displaystyle\int_{a}^{b}f(x)d\mu(x)$
$\displaystyle\leq\int_{a}^{b}(d+L|x-x^{\prime}|)d\mu(x)$
$\displaystyle=d(b-a)+L\int_{a}^{b}|x-x^{\prime}|d\mu(x)$ $\displaystyle\leq
d(b-a)+\frac{1}{2}L(b-a)^{2},$

where the last inequality follows from the fact that
$\int_{a}^{b}|x-x^{\prime}|d\mu(x)\leq\frac{1}{2}(b-a)^{2}$. Since
$\int_{a}^{b}f(x)d\mu(x)=c$, we have:

$\displaystyle c$ $\displaystyle\leq d(b-a)+\frac{1}{2}L(b-a)^{2}$
$\displaystyle\implies d$
$\displaystyle\geq\frac{1}{b-a}\left(c-\frac{1}{2}L(b-a)^{2}\right)$
$\displaystyle=\frac{c}{b-a}-\frac{L(b-a)}{2}.$

Since $x^{\prime}$ was arbitrary, this lower bound holds for all $x\in[a,b]$.
The upper bound can be proved analogously. ∎

######  Lemma A.2.

Let $\mathcal{C}(x)=\\{y:f(y|x)\geq t^{\prime}\\}$. Then the smallest
$t^{\prime}$ satisfying

$\displaystyle\int_{\mathcal{X}}\int_{\mathcal{C}(x)}f(y|x)d\mu(y)dP(x)\geq
1-\alpha$

is given by

$\displaystyle t^{\prime}=\inf\\{t\in\mathbb{R}:\mathbb{P}(f(Y|X)\geq t)\geq
1-\alpha\\}.$

###### Proof.

By direct calculation.

$\displaystyle\mathbb{P}(f(Y|X)\geq
t)=\int_{\mathcal{X}}\int_{\mathcal{Y}}\mathbf{1}\\{f(y|x)\geq
t\\}dP(y|x)dP(x)=\int_{\mathcal{X}}\int_{\\{f(y|x)\geq
t\\}}f(y|x)d\mu(y)dP(x)$

The proof is complete. ∎