# Virtual linear map algorithm
for classical boost in near-term quantum computing

Guillermo García-Pérez<EMAIL_ADDRESS>Algorithmiq Ltd, Kanavakatu 3 C,
FI-00160 Helsinki, Finland  Elsi-Mari Borrelli  Algorithmiq Ltd, Kanavakatu 3
C, FI-00160 Helsinki, Finland  Matea Leahy  Algorithmiq Ltd, Kanavakatu 3 C,
FI-00160 Helsinki, Finland  Joonas Malmi  Algorithmiq Ltd, Kanavakatu 3 C,
FI-00160 Helsinki, Finland  Sabrina Maniscalco  Algorithmiq Ltd, Kanavakatu 3
C, FI-00160 Helsinki, Finland  Matteo A. C. Rossi  Algorithmiq Ltd, Kanavakatu
3 C, FI-00160 Helsinki, Finland  Boris Sokolov  Algorithmiq Ltd, Kanavakatu 3
C, FI-00160 Helsinki, Finland  Daniel Cavalcanti  Algorithmiq Ltd, Kanavakatu
3 C, FI-00160 Helsinki, Finland

###### Abstract

The rapid progress in quantum computing witnessed in recent years has sparked
widespread interest in developing scalable quantum information theoretic
methods to work with large quantum systems. For instance, several approaches
have been proposed to bypass tomographic state reconstruction, and yet retain
to a certain extent the capability to estimate multiple physical properties of
a given state previously measured. In this paper, we introduce the Virtual
Linear Map Algorithm (VILMA), a new method that enables not only to estimate
multiple operator averages using classical post-processing of informationally
complete measurement outcomes, but also to do so for the image of the measured
reference state under low-depth circuits of arbitrary, not necessarily
physical, $k$-local maps. We also show that VILMA allows for the variational
optimisation of the virtual circuit through sequences of efficient linear
programs. Finally, we explore the purely classical version of the algorithm,
in which the input state is a state with a classically efficient
representation, and show that the method can prepare ground states of many-
body Hamiltonians.

##  I Introduction

The last decade has witnessed a tremendous progress in the field of quantum
computing. Today, we have access to several devices composed by tens to
hundreds of physical qubits that can be used as test-beds for small quantum
simulations and even to demonstrate the superiority of quantum information
processing over its classical counterpart [1, 2, 3]. However, the small number
of qubits, limited connectivity, and presence of noise remain strong barriers
for the implementation of quantum protocols offering true practical advantage
over classical computing. Due to the limitations in the quantum hardware,
hybrid algorithms incorporating classical pre- or post-processing
methodologies are often utilised [4, 5]. Despite the promise that such
techniques hold in demonstrating the first useful applications for quantum
computing in the near future, so far such showcases are still missing. Thus,
the development of further classical pre- and post-processing techniques
remains of central importance for the future of quantum computing and
simulation [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

A key limitation that we encounter when post-processing the result of a
quantum computation is the poor statistics obtained in such calculations. Due
to the exponential growth of the Hilbert space, the number of shots one can
access in a typical experiment is not enough to perform quantum state
tomography. This prevents us from studying many properties of the quantum
state produced by the quantum processor. Even determining the expected value
of relevant observables within a satisfactory accuracy may be difficult in
this regime [18, 19, 20, 21]. Because of this, alternative estimation
techniques must be developed [22, 23, 24, 25, 26, 27, 13, 28, 29, 30, 31]

In this article we introduce a classical algorithm that can be used as a post-
processing method to assist computation on quantum devices. The main goal of
the algorithm is to estimate physical properties of a modified state
corresponding to the result of a virtual transformation applied on the state
produced by the quantum processor. By virtual, we mean that these operations
are not physically implemented in the quantum system, but rather implemented
in a classical device. As such, these virtual maps can generally represent
non-physical operations, that is, they do not need to be described by
completely positive maps. The key feature of the algorithm is the fact that
its input is the statistics obtained from an informationally-complete set of
measurements applied on a quantum state, but does not require full state
tomography. We call this algorithm VILMA, from _Virtual Linear Map Algorithm_.

The paper is structured as follows: In Section II we give a general overview
of the algorithm, formalise the VILMA method and provide some mathematical
details. In Section III we show some numerical results on the reconstruction
of multiple expectation values on transformed states. Section IV presents a
methodology to optimise the VILMA maps. In Section V we discuss VILMA as a
purely classical method and, finally, in Section VI we present our conclusion.

##  II The VILMA method

Figure 1: A) A quantum processor produces an $N$-qubit state $\varrho$ that is
then measured by a tomographically complete set of measurements. The
measurement statistics is then post-processed classically through VILMA to
provide an estimate of $\mathrm{Tr}[\Lambda(\varrho)\mathcal{O}]$, where
$\Lambda$ is a linear map and $\mathcal{O}$ an efficiently represented
observable. B-D) Examples of lower depth circuits composed by a sequence of
linear maps acting on few qubits that can be used to define VILMA.

Suppose we have at our disposal a quantum processor that can produce an
$N$-qubit state $\varrho$. To characterise the system fully, estimating the
state $\varrho$ is necessary, but becomes practically impossible even for tens
of qubits. However, even if full state reconstruction is not possible one can
still estimate the expected value $\mathrm{Tr}[\varrho\mathcal{O}_{i}]$ of a
set of observables $\\{\mathcal{O}_{i}\\}$. This can be done directly, i.e. by
measuring the observables, if these observables can be efficiently implemented
in the physical apparatus. Or they can be indirectly estimated from the
results of other measurements [27, 28, 13, 32, 30], in the case that these
observables have an efficient classical representation (for instance,
$\mathcal{O}$ can be written as a linear combination of few Pauli strings).

Here we go a step beyond and provide a method to estimate the expected value
of measurements applied to a modified state
$\varrho^{\prime}=\Lambda(\varrho)$, where $\Lambda$ is a linear map. More
specifically, VILMA allows us to estimate
$\mathrm{Tr}[\Lambda(\varrho)\mathcal{O}_{i}]$ given that we have the
statistics of a tomographically complete set of measurements applied on
$\varrho$. We insist that due to the system’s size we can only implement a low
number of measurement rounds so we do not have enough statistics nor classical
memory to reconstruct the state $\varrho$. We also stress that in order to
perform the computation efficiently on a classical computer, we need to
restrict ourselves to observables that have an efficient classical
representation, for instance, as a linear combination of few Pauli strings.
This is the case for many problems of interest, such as estimating the ground
state of local Hamiltonians.

We will also consider $\Lambda$ to be the composition of $k$-local maps
forming a low-depth circuit. Examples of such maps can be seen in Fig. 1B-D).
As we will see, this assumption will allow us to not only carry on the
calculations efficiently, but also to consider optimisation problems that
variationally adjust each of these $k$-local maps.

In what follows, we describe the methodology used in order to apply virtual
post-measurement maps to the measured state. It should be stated that, while
the main goal of the method is to calculate quantities of the form
$\mathrm{Tr}[\Lambda(\varrho)\mathcal{O}]$ efficiently, the techniques can be
easily generalised for other purposes. We introduce some of these additional
features as well.

###  II.1 Operator averages on transformed states

####  II.1.1 Statistical estimation with IC measurement data

Suppose the quantum processing unit (QPU) is in some $N$-qubit state
$\varrho$, which is measured by a $k$-local informationally complete (IC)
measurement, that is given by an IC-POVM whose statistics singles out the
state on which the measurement is applied [33]. By $k$-local here we mean that
it is composed of POVM elements that act non-trivially on $k$ qubits at most.
In principle, the IC-POVM need not be qubit-local nor minimal (i.e. contain
the minimum number of POVM elements), so that they can in principle act on
$k>1$ qubits and be overcomplete. However, the number of qubits $k$ they act
upon must be small enough so that operators of dimension $2^{k}\times 2^{k}$
can be dealt with classically. Moreover, the number of POVM elements (which
correspond to the number of outcomes of the measurement) should also be small
enough so that one can store the measurement statistics in classical memory.
In what follows we will assume $k=1$ and that the POVMs have the minimal
number of elements to be considered IC, although all the results included here
can be easily generalised to larger $k$ and non-minimal measurements. In
particular, we will consider POVMs of the form
$\\{\Pi_{\mathbf{m}}=\bigotimes_{i=1}^{N}\Pi_{m_{i}}^{(i)}\\}$, where each
$\\{\Pi_{m_{i}}^{(i)}\\}$ is a four-outcome IC-POVM on the Hilbert space of
qubit $i$.

Given such a local IC-POVM, one can find a set of dual effects
$\\{D_{n_{i}}^{(i)}\\}$ satisfying
$\mathcal{O}=\sum_{\mathbf{m}}\mathrm{Tr}[\mathcal{O}\Pi_{\mathbf{m}}]D_{\mathbf{m}}$
for all $\mathcal{O}$. A general method to compute the set of duals for a
given POVM can be found in Ref. [34]. The decomposition of an $N$-qubit state
$\varrho$ in terms of these dual effects therefore reads
$\varrho=\sum_{\mathbf{m}}p_{\mathbf{m}}D_{\mathbf{m}}$, where
$D_{\mathbf{m}}=\bigotimes_{i=1}^{N}D_{m_{i}}^{(i)}$ and
$p_{\mathbf{m}}=\mathrm{Tr}[\varrho\Pi_{\mathbf{m}}]$ is the probability of
obtaining outcome $\mathbf{m}$ once the measurement is performed.

We also assume that the observables of interest admit an efficient classical
representations known to us. In particular, we assume that they are of the
form $\mathcal{O}=\sum_{\mathbf{k}}c_{\mathbf{k}}P_{\mathbf{k}}$, where
$c_{\mathbf{k}}\in\mathbb{C}$ and
$P_{\mathbf{k}}=\bigotimes_{i=1}^{N}P_{k_{i}}^{(i)}$. Each $P_{k_{i}}^{(i)}$
is a single-qubit operator, such as a Pauli operator. In fact, in most
applications, $\mathcal{O}$ is given in terms of such a linear combination of
Pauli strings. Importantly, for a wide class of relevant problems in many-body
physics, the observables of interest can be mapped into low-weight $N$-qubit
operators, meaning that only a small fraction of the Pauli operators in each
string is different from identity. Their corresponding observable averages can
be estimated using local IC-POVM with polynomially scaling number of
measurements [13, 28, 30].

In experiments, only a finite number of measurement rounds $S$ can be
performed. In this case such a procedure will result in a random sequence of
outcomes $(\mathbf{m}_{1},\ldots,\mathbf{m}_{S})$, according to which we can
write a crude approximation to the state,
$\varrho_{S}=\sum_{i=1}^{S}D_{\mathbf{m}_{i}}/S$. Notice that
$\lim_{S\rightarrow\infty}\varrho_{S}=\varrho$.

Assuming a linear map $\Lambda$ that is independent of $\varrho$ and
$\mathcal{O}$, we can write

$\displaystyle\mathrm{Tr}[\Lambda(\varrho)\mathcal{O}]$ $\displaystyle=$
$\displaystyle\lim_{S\rightarrow\infty}\mathrm{Tr}[\Lambda(\varrho_{S})\mathcal{O}]$
(1) $\displaystyle=$
$\displaystyle\lim_{S\rightarrow\infty}\sum_{i=1}^{S}\frac{1}{S}\sum_{\mathbf{k}}c_{\mathbf{k}}\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}].$
(2)

This implies that the quantity
$\bar{\mathcal{O}}_{\Lambda}=\sum_{i=1}^{S}\omega_{\mathbf{m}_{i}}/S$, with
$\omega_{\mathbf{m}_{i}}=\sum_{\mathbf{k}}c_{\mathbf{k}}\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$,
is a consistent estimator of the mean value of $\mathcal{O}$ for state
$\Lambda(\varrho)$, that is,
$\lim_{S\rightarrow\infty}\bar{\mathcal{O}}_{\Lambda}=\mathcal{O}_{\Lambda}\equiv\mathrm{Tr}[\Lambda(\varrho)\mathcal{O}]$.
In addition, it is easy to see from the linearity of the expression that
$\bar{\mathcal{O}}_{\Lambda}$ is also unbiased, meaning that its average over
$S$-measurement sampling experiments is also equal to the observable average,
$\mathbb{E}_{S}[\bar{\mathcal{O}}_{\Lambda}]=\mathcal{O}_{\Lambda}$.
Accordingly, the mean squared error of the estimation is given by
$\mathbb{E}_{S}[(\bar{\mathcal{O}}_{\Lambda}-\mathcal{O}_{\Lambda})^{2}]=\mathrm{Var}(\omega_{\mathbf{m}})/S$,
where $\mathrm{Var}(\omega_{\mathbf{m}})$ is the variance of
$\omega_{\mathbf{m}}$ over the probability distribution of the outcomes. In
realistic scenarios, the variance $\mathrm{Var}(\omega_{\mathbf{m}})$ is not
known, but it can be estimated from the measurement outputs using the unbiased
estimator
$\bar{V}(\mathcal{O}_{\Lambda})\equiv[\sum_{i=1}^{S}\omega_{\mathbf{m}_{i}}^{2}/S-(\bar{\mathcal{O}}_{\Lambda})^{2}]S/(S-1)$.
In short, given $S$ IC-POVM outcomes from state $\varrho$, computing
$\bar{\mathcal{O}}_{\Lambda}$ and
$\sigma\equiv\sqrt{\bar{V}(\mathcal{O}_{\Lambda})/S}$ enables the estimation
of the observable average $\mathcal{O}_{\Lambda}$ and the corresponding
statistical error, respectively. Regarding the latter, notice that while it is
possible to derive analytical bounds for the statistical error based on the
weight of the Pauli strings in the operator $\mathcal{O}$ when using IC POVM-
based estimators [13, 28, 32, 30], the calculations cannot be
straightforwardly generalised to VILMA (in fact, the statistical error must
depend on the details of the map $\Lambda$ as well). It is nevertheless
possible to guarantee the polynomial scaling of the variance for $k$-local
observables and certain maps.

In practice, a major difficulty in dealing with the above terms lies in
computing the terms $\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$.
In some cases, this may be easy to do, for instance if the map
$\Lambda(\cdot)=\sum_{ij}\lambda_{ij}B_{i}\cdot B_{j}^{\dagger}$ only involves
a moderate number of terms $\lambda_{ij}$ and the operators $B_{i}$ have
bounded locality; one can then simply compute
$\mathrm{Tr}[D_{\mathbf{m}_{i}}\Lambda^{\dagger}(P_{\mathbf{k}})]$ in
polynomial time, and the variance will be polynomially bounded too. However,
if the map has a complex structure, the computation may be much more
challenging. In particular, notice that the explicit representation of
$\Lambda(D_{\mathbf{m}_{i}})$ can generally be classically prohibitive, even
for modest $N$. In the next section, we provide an algorithm to compute these
traces for a class of maps of particular relevance for quantum computing,
namely circuits of $K$-local maps. We show that, by restricting the circuit
structure, we can make such a computation classically amenable.

####  II.1.2 Traces involving mapped dual effects

The main idea of VILMA is to restrict the structure and complexity of the map
$\Lambda$ so that the calculation of each term
$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$ only involves dealing
with operators of bounded and efficient dimension during all the intermediate
steps. More precisely, we ensure that the trace in
$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$ can be computed
through intermediate partial traces that keep the dimension of non-trivial
operators (meaning, those that are not tensor products of single-qubit ones)
under control. The details will become clear in what follows.

Figure 2: A circuit decomposition of a VILMA map composed by two layers of
sequential 2-qubit gates. $q_{k}$ refers to qubit $k$. A) In the first step of
the algorithm we apply the maps $\Lambda_{12}^{1}$, $\Lambda_{23}^{2}$ and
$\Lambda_{12}^{2}$ (green boxes), which correspond to the causal cone of qubit
1 (see main text), to $D_{m_{1}}\otimes D_{m_{2}}\otimes D_{m_{3}}$. After
multiplying with $P_{k_{1}}\otimes\mathbb{1}\otimes\mathbb{1}$ and tracing out
qubit 1, we are left with a residual operator $R_{23}$ on qubits $2$ and $3$,
and move to the next step. B) In the second step we apply maps
$\Lambda_{34}^{1}$ and $\Lambda_{23}^{2}$, which correspond to the causal cone
of qubit 2, to $R_{23}\otimes D_{m_{4}}$. After multiplying the resulting
operator with $P_{k_{2}}\otimes\mathbb{1}\otimes\mathbb{1}$, we trace out
qubit 2 to obtain a residual operator $R_{34}$ on qubits $3$ and $4$, and move
to the next step. The algorithm then proceeds similarly, by applying two-qubit
maps on three-qubits and performing partial traces until all the maps are
applied and the traces are performed.

We proceed by considering maps $\Lambda$ that can be decomposed in terms of a
circuit of $k$-qubit maps with $k\leq K$ for some $K$, and with an adequate
causal structure. For the sake of clarity, in what follows we will present
this idea using a concrete example of maps composed of two layers of 2-qubit
gates, shown in Figure 2:

$\Lambda=\cdots\circ\Lambda_{34}^{2}\circ\Lambda_{23}^{2}\circ\Lambda_{12}^{2}\circ\cdots\circ\Lambda_{34}^{1}\circ\Lambda_{23}^{1}\circ\Lambda_{12}^{1}.$
(3)

Indices $ij$ in $\Lambda_{ij}^{l}$ indicate that the map is applied to qubits
$ij$, and index $l$ indicates that it is part of the $l^{th}$ layer. For such
VILMA map, $\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$ can be
computed with no need for calculating any $m$-qubit operator with $m>3$ at any
point (even if the same structure is extended to $N>4$ qubits). To see how
this can be done, first notice that the 2-qubit maps $\Lambda_{ij}^{l}$ acting
on different qubits commute, so that we can write, for instance,

$\Lambda=\cdots\circ\Lambda_{23}^{2}\circ\Lambda_{34}^{1}\circ\Lambda_{12}^{2}\circ\Lambda_{23}^{1}\circ\Lambda_{12}^{1}.$
(4)

This ordering presents the following advantage: the first three maps act non-
trivially only on qubits 1, 2, and 3, and the rest of them act trivially on
qubit 1. Thus, we can first compute
$\Lambda_{12}^{2}\circ\Lambda_{23}^{1}\circ\Lambda_{12}^{1}(D_{m_{1}}\otimes
D_{m_{2}}\otimes D_{m_{3}})$, which requires the explicit computation of a
three-qubit operator. Now, the only non-trivial operation on qubit 1 is the
multiplication with $P_{k_{1}}$, after which no other operation takes place in
such Hilbert space when computing $\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}$.
Therefore, we can simply trace out qubit 1, and keep track of the resulting
residual operator

$R_{2,3}=\mathrm{Tr}_{1}[\Lambda_{12}^{2}\circ\Lambda_{23}^{1}\circ\Lambda_{12}^{1}(D_{m_{1}}\otimes
D_{m_{2}}\otimes
D_{m_{3}})P_{k_{1}}\otimes\mathbb{I}_{2}\otimes\mathbb{I}_{3}]$ (5)

In Appendix A, we present a more formal explanation of this step.

We can then proceed with maps $\Lambda_{34}^{1}$ and $\Lambda_{23}^{2}$, which
must be applied to $R_{2,3}\otimes D_{m_{4}}$, yielding another three-qubit
operator. Once again, we can trace out qubit 2 after multiplying the
corresponding $P_{k_{2}}$ and keep track of the residual operator
$R_{3,4}=\mathrm{Tr}_{2}[\Lambda_{23}^{2}\circ\Lambda_{34}^{1}(R_{2,3}\otimes
D_{m_{4}})P_{k_{2}}\otimes\mathbb{I}_{3}\otimes\mathbb{I}_{4}]$.

This method can then be iterated until all the traces are performed. Notice
that the algorithm can be applied backwards if we swap the roles of
$D_{\mathbf{m}_{i}}$ and $P_{\mathbf{k}}$ and use the adjoint map
$\Lambda^{\dagger}$, since
$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]=\mathrm{Tr}[\Lambda^{\dagger}(P_{\mathbf{k}})D_{\mathbf{m}_{i}}]$.
The adjoint map $\Lambda^{\dagger}$ can be easily obtained by mirroring the
circuit (i.e. reversing the order of the maps) and substituting each map with
its adjoint. Depending on the map, the forwards or backwards algorithms can
have different computational complexity. It is also possible to combine them
in the same calculation, as will be discussed in Section IV.

The key feature that makes it possible to split the calculation into a
sequence of few-qubit calculations is the causal cone structure of the VILMA
maps represented in Fig. 2. More precisely, the past causal cone of qubit 1
(that is, the maps on which the reduced operator
$\mathrm{Tr}_{2,\ldots,6}[\Lambda(D_{\mathbf{m}_{i}})]$ depends) only involves
maps $\Lambda_{12}^{1}$, $\Lambda_{23}^{1}$, and $\Lambda_{12}^{2}$, and hence
qubits 1, 2, and 3. Likewise, the past causal cone of qubit 2 additionally
involves maps $\Lambda_{34}^{1}$ and $\Lambda_{23}^{2}$, that is, it requires
the inclusion of qubit 4. However, since we can previously trace out qubit 1
in the calculation, we only need to take into account non-trivial operators in
the joint Hilbert space of qubits 2, 3, and 4 in this step.

In general, other circuit structures may be used for the maps in the
algorithm, but the efficiency of the algorithm highly depends on the causal
structure of the circuit. For the structure considered above, the scaling of
the method is linear in the number of qubits $N$. However, as we add more
layers to the map circuit (that is, sequences of maps $\Lambda_{ij}^{3}$,
$\Lambda_{ij}^{4}$, etc), the past causal cone of each qubit involves more
qubits (e.g., four qubits for three layers, and so on), which results in
higher-dimensional intermediate operators and thus higher computational cost.
Moreover, recall that the sequential computation presented above must be
performed for all POVM outcomes $\mathbf{m}_{i}$ and terms $\mathbf{k}$ in
$\mathcal{O}$.

##  III Numerical simulations

Let us illustrate the method with the following numerical experiment. As a
reference state, that is, the one in the QPU, we consider a perturbed version
of the ground state of the Hamiltonian $H$ of the hydrogen molecule
$\rm{H}_{2}$ in the 6-31G basis with a stretched geometry (bond distance
$d=2$Å) mapped to qubit space using the Jordan-Wigner transformation. The
resulting Hamiltonian can be encoded with eight qubits. As a perturbation, we
consider the image $\varrho$ of the ground state under a sequence of channels
$\mathcal{N}_{ij}^{1}$ forming a circuit like that in Fig. 1A). In this way,
by construction, there exists a single-layer VILMA circuit capable of
cancelling the perturbation. To control the magnitude of the perturbation, the
$\mathcal{N}_{ij}^{1}$ are chosen as
$\mathcal{N}_{ij}^{1}=(1-p)\mathbb{1}+p\mathcal{E}_{ij}^{1}$, where
$\mathcal{E}_{ij}^{1}$ is randomly sampled from the space of two-qubit CPTP
maps, and $p$ is a parameter that we fix to $p=0.05$.

Next, we sample from $\varrho$ using a POVM composed of single-qubit symmetric
IC-POVM. We sample 50 realisations of $S=10^{4}$ samples each. Our aim is to
show that, for each such batch of $S=10^{4}$ samples, we can compute multiple
expectation values, for different operators and maps. In particular, we
consider three maps: identity, as a reference with the original perturbed
state, the inverse of the perturbation, to illustrate that non-physical
operations can be used and, finally, a layer of randomly chosen unitaries,
which could represent a sequence of gates that one may consider applying on
the quantum computer. As operators, we consider the Hamiltonian, along with
the real parts of 1-body and 2-body reduced density matrices,
$\mathrm{Re}\langle a_{i}^{\dagger}a_{j}\rangle$ and $\mathrm{Re}\langle
a_{i}^{\dagger}a_{j}a_{k}^{\dagger}a_{l}\rangle$.

Figure 3: Statistical estimations of three observables, the Hamiltonian $H$, a
1-RDM element $\mathrm{Re}\langle a_{2}^{\dagger}a_{7}\rangle$, and a 2-RDM
element $\mathrm{Re}\langle a_{0}^{\dagger}a_{1}a_{4}^{\dagger}a_{5}\rangle$,
using three maps, the identity map $\Lambda_{0}$, the inverse of the
perturbation map $\Lambda_{1}$, and a layer of random unitaries $\Lambda_{2}$,
over 50 realisations of $S=10^{4}$ measurement rounds each. The histograms
show the distribution of estimates $\bar{\mathcal{O}}_{\Lambda}$ for two
different pairs of observable and map, along with the exact value (red
vertical line). The blue vertical lines indicate the mean of the estimations
over the 50 realisations.  Figure 4: Statistical estimations of three
observables, the Hamiltonian $H$, a 1-RDM element $\mathrm{Re}\langle
a_{2}^{\dagger}a_{7}\rangle$, and a 2-RDM element $\mathrm{Re}\langle
a_{0}^{\dagger}a_{1}a_{4}^{\dagger}a_{5}\rangle$, using three maps, the
identity map $\Lambda_{0}$, the inverse of the perturbation map $\Lambda_{1}$,
and a layer of random unitaries $\Lambda_{2}$, over 50 realisations of
$S=10^{4}$ measurement rounds each. A–C) With the sampling data obtained in
every realisation, we compute $\bar{\mathcal{O}}_{\Lambda}$ for every
combination of map and observable and, with it, the error in the estimation,
$|\bar{\mathcal{O}}_{\Lambda}-\mathcal{O}_{\Lambda}|$. We also compute each
estimate of the error $\sigma$. With high probability, the error
$|\bar{\mathcal{O}}_{\Lambda}-\mathcal{O}_{\Lambda}|$ should be smaller than
$3\sigma$. To show that this is indeed the case, for every realisation, and
for every pair $(\mathcal{O},\Lambda)$, we draw a point
$(3\sigma,|\bar{\mathcal{O}}_{\Lambda}-\mathcal{O}_{\Lambda}|)$. It can be
seen that in nearly all cases, points lie below the diagonal (dashed line).

In Fig. 3, we show the resulting distribution of $\bar{\mathcal{O}}_{\Lambda}$
for the different maps and observables. The red and blue vertical lines
indicate the exact values $\mathcal{O}_{\Lambda}$ and the average of
$\bar{\mathcal{O}}_{\Lambda}$ over realisations, respectively. The obtained
values fluctuate around the correct value. Importantly, as discussed in
Section II, the method also produces, along with the estimate
$\bar{\mathcal{O}}_{\Lambda}$, an estimation of the error incurred,
$\sigma\equiv\sqrt{\bar{V}(\mathcal{O}_{\Lambda})/S}$. To assess whether the
error estimations $\sigma$ are meaningful, we compare, for each estimation,
the actual error and the estimated one. In particular, the error
$|\bar{\mathcal{O}}_{\Lambda}-\mathcal{O}_{\Lambda}|$ should be smaller than
$3\sigma$ with high probability. In Fig. 4 we show that, indeed, for nearly
all estimations, this is the case.

Importantly, while VILMA allows us to reuse the IC measurement data to
estimate many expectation values on many different states, these estimations
are statistically correlated. Thus, one should keep in mind that
$\sqrt{\bar{V}(\mathcal{O}_{\Lambda})/S}$ yields an estimate of the error of
$\bar{\mathcal{O}}_{\Lambda}$, but two different estimations
$\bar{\mathcal{O}}_{\Lambda}$ and
$\bar{\mathcal{O^{\prime}}}_{\Lambda^{\prime}}$ will generally have non-zero
covariance.

##  IV Variational optimisation with VILMA

For some applications, VILMA may be used as a classical boost to a variational
calculation. In such situations, we are interested in finding a map $\Lambda$
that minimises or maximises an observable average $\mathcal{O}_{\Lambda}$ for
the given input state. In order to carry out the optimisation efficiently, we
can also use the specific structure of the VILMA circuit to evaluate the
observable average as a function of a single map component $\Lambda_{ij}^{l}$,
while keeping all other components fixed. It is therefore unnecessary to
repeat the whole algorithm presented above for all the other maps at each
optimisation step if only $\Lambda_{ij}^{l}$ is modified. In fact, it is
possible to write down an expression, linear in $\Lambda_{ij}^{l}$, for
$\bar{\mathcal{O}}_{\Lambda}$ that only depends on operators defined in the
local Hilbert space of qubits $i$ and $j$. In what follows, we outline how
this can be done.

Suppose we want to optimise an observable average with respect to a specific
map component $\Lambda_{ij}^{l}$. As discussed above, VILMA can be applied
both in the forwards and backwards directions. If both directions can be
applied efficiently, as for the symmetric structure discussed in the previous
section, the strategy is to apply the algorithm forwards until the stage in
which $\Lambda_{ij}^{l}$ would be applied, and then backwards for the rest of
the maps, so that eventually we arrive at an expression of the form
$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]=\sum_{a}\mathrm{Tr}[\Lambda_{ij}^{l}(R_{a}^{(\mathbf{m}_{i},\mathbf{k})})\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}]$.
The estimator as a function of $\Lambda_{ij}^{l}$ then reads

$\bar{\mathcal{O}}_{\Lambda}=\sum_{i=1}^{S}\frac{1}{S}\sum_{\mathbf{k}}c_{\mathbf{k}}\sum_{a}\mathrm{Tr}[\Lambda_{ij}^{l}(R_{a}^{(\mathbf{m}_{i},\mathbf{k})})\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}]$
(6)

The linearity of this expression enables using linear programming, including
semidefinite programming (SDP) if, for instance, one is interested in imposing
the complete positivity of the map $\Lambda_{ij}^{l}$. We illustrate the
derivation of the above expression with an example in what follows, but the
generalisation to other situations should be clear.

Consider again the map in Fig. 2 and let us single out the map
$\Lambda_{2,3}^{2}$. In the calculation of
$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]$, the forward
algorithm starts by computing $R_{2,3}$ (Eq. (5)). Next, we would need to
apply $\Lambda_{34}^{1}$ and $\Lambda_{23}^{2}$ to $R_{2,3}\otimes D_{m_{4}}$.
Instead, however, we apply only $\Lambda_{34}^{1}$ and calculate explicitly

$R_{2,3,4}=\mathbb{1}_{q_{2}}\otimes\Lambda_{34}^{1}(R_{2,3}\otimes
D_{m_{4}}).$ (7)

The backwards algorithm is then applied until qubit 5 is traced out, leaving
the residual operator $\bar{R}_{3,4}$. Let us define

$\bar{R}_{2,3,4}=P_{k_{2}}\otimes\bar{R}_{3,4}.$ (8)

Notice that we can now write

$\mathrm{Tr}[\Lambda(D_{\mathbf{m}_{i}})P_{\mathbf{k}}]=\mathrm{Tr}[\Lambda_{23}^{2}\otimes\mathbb{1}_{q_{4}}(R_{2,3,4})\bar{R}_{2,3,4}].$
(9)

Finally, if we write
$R_{2,3,4}=\sum_{a}R_{a}^{(\mathbf{m}_{i},\mathbf{k})}\otimes B_{a}$ and
$\bar{R}_{2,3,4}=\sum_{a}\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}\otimes
B_{a}$, where $\\{B_{a}\\}$ is the normalised Pauli basis in the Hilbert space
of qubit 4, on which $\Lambda_{23}^{2}$ acts trivially, we see that
$\mathrm{Tr}[\Lambda_{23}^{2}\otimes\mathbb{1}_{q_{4}}(R_{2,3,4})\bar{R}_{2,3,4}]=\sum_{a,a^{\prime}}\mathrm{Tr}[\Lambda_{23}^{2}(R_{a}^{(\mathbf{m}_{i},\mathbf{k})})\bar{R}_{a^{\prime}}^{(\mathbf{m}_{i},\mathbf{k})}\otimes
B_{a}B_{a^{\prime}}]=\sum_{a}\mathrm{Tr}[\Lambda_{23}^{2}(R_{a}^{(\mathbf{m}_{i},\mathbf{k})})\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}]$.
This procedure allows us to pre-compute a set of low-dimensional operators,
$\\{R_{a}^{(\mathbf{m}_{i},\mathbf{k})},\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}\\}_{a,(\mathbf{m}_{i},\mathbf{k})}$,
that capture all the dependence of $\bar{\mathcal{O}}_{\Lambda}$ on a specific
map component $\Lambda_{ij}^{l}$.

The algorithm above can be summarised in the following general steps: 1) run
the forward algorithm, stopping right before the singled-out map $\Lambda_{s}$
is applied. Decompose the resulting operator $R$ as
$R=\sum_{a}R_{a}^{(\mathbf{m}_{i},\mathbf{k})}\otimes B_{a}$, where
$\\{B_{a}\\}$ is the normalised Pauli basis in the Hilbert space of the qubits
on which $\Lambda_{s}$ acts trivially, and store
$\\{R_{a}^{(\mathbf{m}_{i},\mathbf{k})}\\}$. 2) run the backward algorithm,
stopping right before the singled-out map $\Lambda_{s}$ is applied. Decompose
similarly the resulting operator $\bar{R}$ as
$\bar{R}=\sum_{a}\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}\otimes B_{a}$ and
store $\\{\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}\\}$. Repeating the process
for every measurement outcome and Pauli string pair
$(\mathbf{m}_{i},\mathbf{k})$ allows us to use Eq. (6) for optimisation.

Figure 5: Relative error between the energy of the image of the input state
through the VILMA map, $\bar{\mathcal{O}}_{\Lambda}$, and the ground state
energy, $E_{0}$, along the optimisation process for two input states, a noisy
VQE simulation and $\ket{0}^{\otimes N}$, and four XX model Hamiltonians, for
two system sizes and two values of the magnetic field. A) Right below the
phase transition of the model, the energy converges to within the numerical
precision of the optimiser, regardless of the input state. The number of
iterations required changes notably from one example to another. For $N=6$,
convergence is reached faster with the input vacuum state than with the noisy
VQE input, while the opposite is true for $N=7$. B) For the zero-field XX
model, the optimisation does not reach the ground state energy. However, with
the noisy circuit input state, VILMA decreases the energy with respect to the
initial reference state and, in both cases, below the energy that can be
reached using the classical input state $\ket{0}^{\otimes N}$.

We apply this procedure on a single-layer VILMA circuit for four different XX
model Hamiltonians,
$H=-J[\sum_{i}(\sigma_{x}^{(i)}\sigma_{x}^{(i+1)}+\sigma_{y}^{(i)}\sigma_{y}^{(i+1)})/2+B\sigma_{z}^{(i)}]$,
where $J$ is the coupling constant, that we set to one, $B$ is the magnetic
field, and $\sigma_{k}^{(i)}$ with $k=x,y,z$ are Pauli matrices. Periodic
boundary conditions, $\sigma_{k}^{(N+1)}=\sigma_{k}^{(1)}$, are assumed. The
cases that we consider correspond to $B=0$ and $B\approx 1$ (right below the
phase transition exhibited by the model in the thermodynamic limit [35]). At
finite size, the ground state exhibits a series of level crossings as a
function of $B$, as a result of which the entanglement structure of the state
changes notoriously [35, 36]. We also use two different system sizes, $N=6$
and $N=7$ spins.

As reference states we take the outputs of noisy circuits. The circuits are
pre-trained VQE ansätze that, in the absence of noise, would approximate the
ground states [5]. However, we apply noisy gates. Every CNOT is composed with
a channel $\mathcal{N}_{\mathrm{CNOT}}$, which itself is a composition of
depolarising noise (with $p=10^{-3}$) and coherent noise (the two qubits are
rotated by applying two single-qubit rotations $R_{x}(\theta)R_{z}(\theta)$
with an angle $\theta=0.05$ to each of them.

The single-layer VILMA map is optimised in the following way. Every map
$\Lambda_{ij}^{1}$ is chosen sequentially and, for each of them, we compute
the operators
$\\{R_{a}^{(\mathbf{m}_{i},\mathbf{k})},\bar{R}_{a}^{(\mathbf{m}_{i},\mathbf{k})}\\}_{a,(\mathbf{m}_{i},\mathbf{k})}$
to then minimise $\bar{\mathcal{O}}_{\Lambda}$ over $\Lambda_{ij}^{1}$. In the
optimisation, we impose each map to be CPTP to ensure that the output operator
is a valid state. The resulting constrained optimisation of Eq. (6) is an SDP,
which can be solved efficiently. Since we are now only interested in
illustrating the local optimisation strategy described earlier, in this
artificial numerical experiment we choose to bypass the additional
complications stemming from the optimisation with finite statistics (e.g.
over-fitting) by considering the exact noisy state as input, that is,
$\varrho$ instead of $\varrho_{S}$. In practice, however, finite statistics-
related issues need to be handled.

In Fig. 5, we plot the relative error of the resulting step at each iteration
(solid lines). For $B\approx 1$, the optimisation finds the ground state
within the numerical precision of the optimiser for both system sizes. For
$B=0$, on the other hand, it finds lower energies, but does not reach the
ground state. Importantly, recall that the optimisation is constrained to the
complete positivity of $\Lambda$ because we are neglecting all information
about the origin of the input state and still enforcing the positivity of the
output one. As a consequence, even if a single-layer map mapping the input
state to the ground state exists, it may not be CP, and hence remain
inaccessible to the optimisation procedure. In practice, one should use
strategies in which information about the noise in the device is taken into
account, which is in principle possible given that VILMA per se does not
require the positivity of the maps. In any case, notice that even with one
layer of maps and modest classical compute (indeed, we never reconstruct
operators of dimension larger than $4\times 4$), the algorithm could boost the
result obtained with a quantum computer. Once again, we stress that, in real
experiments, this will require managing additional finite statistics-related
matters.

##  V VILMA as a classical ansatz

So far, we considered VILMA as a classical algorithm that takes the
measurement statistics from a quantum processor as an input. Another
possibility would be to input the classical description of a quantum state
(for instance $\ket{0}^{\otimes N}$) and look for the optimal VILMA maps that,
say, produce a state achieving the minimum expected value of a given
observable (e.g. a Hamiltonian). In this way, VILMA can be used as an
algorithm to produce a classically efficient ansatz.

In Fig. 5, we compare this technique with the results obtained in the previous
section. On the one hand, we can see that, for $B$ nearing the transition,
VILMA can find the ground state even with the $\ket{0}^{\otimes N}$ input
state. However, if the maps are initially set to identity, the local
optimisation sweeps are futile, as the system seems to be at a local minimum.
By initialising the VILMA map with randomly chosen unitary gates, the problem
is overcome. This indicates that there is potential room for improvement in
the optimisation strategy. In fact, finding the ground state with a single-
layer of VILMA should be expected, given that the ground state of the
Hamiltonian is a W state, which can be prepared by a single layer of unitary
gates [37].

We also notice that in Ref. [38] it is shown that, in the case of a layered
circuit of two-qubit maps like the one employed here set to be unitary
operations, the produced state can encode matrix product states [39] whose
bond dimension depends on the number of layers in the circuit. The purely
classical version of VILMA can go beyond this in two ways: (i) it is not
constrained to unitary maps, and (ii) different topologies can be considered.
We leave as an open question the precise characterisation of the class of
states reachable with classical VILMA and the connection with MPS and other
tensor network methods.

In any case, it is worth discussing the connection between VILMA with
classical and quantum inputs. In particular, we make the following
observation: if a state $\varrho^{\prime}$ is classically reachable with
VILMA, that is, if we can access it by applying an affordable VILMA map
$\Lambda$ to $\ket{0}^{\otimes N}$, then $\varrho^{\prime}$ is reachable with
the same VILMA map structure if the input is the state of the quantum
processor $\varrho$. This is so because it is possible to compose the maps in
the first layer of VILMA with single-qubit CPTP channels that map any input
state to $\ket{0}$ in such a way that the $N$-qubit input state is initially
mapped to $\ket{0}^{\otimes N}$ 111For instance, for a layered structure like
the one in Fig. 2, this can be done in the following way. Let $\Lambda$ be a
map such that $\varrho^{\prime}=\Lambda(\outerproduct{0}{0}^{\otimes N})$, and
$\mathcal{Z}$ be the single-qubit CPTP map defined by
$\mathcal{Z}(\sigma)=\outerproduct{0}{0}$. Now, consider the map
$\Lambda^{\prime}$ with the same structure as $\Lambda$ and in which
$\Lambda_{ij}^{{}^{\prime}l}=\Lambda_{ij}^{l}$ for $l>1$. For the first layer
$l=1$, let
$\Lambda_{12}^{{}^{\prime}1}=\Lambda_{12}^{1}\circ(\mathcal{Z}_{1}\otimes\mathcal{Z}_{2})$
and
$\Lambda_{ij}^{{}^{\prime}1}=\Lambda_{ij}^{1}\circ(\mathbb{1}_{i}\otimes\mathcal{Z}_{j})$
for all other maps in the layer. The map $\Lambda^{\prime}$ maps any input
state to $\varrho^{\prime}$.. The converse is clearly not true, since the
input data itself can correspond to a state not classically reachable by
VILMA. Therefore, VILMA may enable computing observable averages that may be
inaccessible using only classical methods or only limited quantum resources.
The results in Fig. 5 B) can be interpreted as an example of this: the
algorithm can benefit from the correlations in the input state to reach states
(in this case, lower-energy ones) that would not be accessible with the same
structure using a product input state. In particular, in the simulation for
$N=7$, the state resulting from the noisy circuit execution has a higher
energy than what can be achieved with a single-layer VILMA map. While one
would deem such quantum computer output useless, given that a better result is
achievable classically, the algorithm exploits the data to achieve a lower
energy than with the separable input.

##  VI Conclusions

In this paper we have introduced VILMA, a scalable algorithm to efficiently
apply a sequence of (not necessarily physical) linear maps on quantum states
and estimate the result of measurements applied on this modified state. Most
importantly, VILMA can be used in the case where the states under scrutiny are
not known, but we have access to just a finite sample of outcomes of a
tomographically complete set of measurements applied to these states. Yet, the
resulting estimates of expectation values are unbiased, and the algorithm
provides a meaningful estimate of the error incurred due to finite statistics.
This allows to use VILMA to post-process the results of a quantum processor in
a regime where performing full state quantum tomography is out of reach.

Besides estimating mean values of post-processed states, VILMA also enables
optimising over the sequence of linear maps, opening the possibility for
applying it to several ends, such as noise mitigation and variational search
of Hamiltonian ground states. Our results suggest that, on the one hand, this
method has a very natural classical counterpart, in which the input state is a
product state, that may have interesting connections with matrix product
states or other tensor network methods. We have shown that purely classical
VILMA can indeed find ground states of many-body Hamiltonians efficiently in
some cases. On the other hand, the algorithm can take as input the IC
measurement data from a quantum computer and boost the results, and we have
discussed that, setting optimisation-related issues aside, the class of states
reachable with such inputs is strictly larger than using product input states.
We have provided proof-of-concept simulation results in which noisy circuit
states are used to achieve lower energies than can be found with product input
states.

There are several open questions related to VILMA that we leave for future
investigations. First of all, one must be careful when using VILMA
variationally on applications involving real experimental data to not over-fit
the results due to statistical fluctuations. One possibility to avoid this is
to perform the optimisation of the maps on a different set of experimental
data than the ones used for the final estimation of the expected values of
interest. Another interesting line of research is to study to what extent
VILMA can be used as a noise mitigation scheme for realistic types and levels
of noise in current hardware. Finally, the connection between VILMA, MPS and
other ansätze is certainly worth further investigation.

Acknowledgements: We acknowledge Zoltán Zimborás, Adam Glos, and Stefan Knecht
for interesting discussions. Competing interests: Elements of this work are
included in a patent filed by Algorithm Ltd with the European Patent Office.
Authors contributions: GGP conceived the algorithm. EMB, DC, GGP designed and
directed the research. GGP, ML, JM, MACR, BS implemented the algorithm and ran
simulations. EMB, DC, GGP wrote the first version of the manuscript. All
authors contributed to scientific discussions and to the writing of the
manuscript.

##  Appendix A Further mathematical details

The main text intended to introduce the method to calculate the traces in an
intuitive manner with an explicit example. In this section, we formalise the
ideas in more detail.

Let $\Lambda$ be a linear map acting in the space of linear operators
$L(\bigotimes_{i=1}^{N}\mathcal{H}_{i})$ of a set of $N$ quantum systems
$\\{q_{1},\ldots,q_{N}\\}$. The map is composed of $L$ $k$-local maps,
$\Lambda=\bigcirc_{r=1}^{L}\Lambda_{r}\equiv\Lambda_{L}\circ\cdots\circ\Lambda_{2}\circ\Lambda_{1}$,
that is, $\Lambda_{r}$ acts non-trivially on at most $k$ subsystems. Since
some of the map components may commute with one another, the map $\Lambda$ may
admit multiple such decompositions. In particular, let $\mathcal{T}=\\{T\in
S_{L}|\bigcirc_{r=1}^{L}\Lambda_{T(r)}=\Lambda\\}$, where $S_{L}$ is the
permutation group of $L$ elements, be the set of all possible permutations of
map components that preserve $\Lambda$.

Consider an element $T\in\mathcal{T}$, and an integer $1\leq l<L$ such that
$\bigcirc_{r=l+1}^{L}\Lambda_{T(r)}$ acts trivially on some subsystem $q_{i}$.
Denote $\Lambda^{\bar{q}_{i}}\equiv\bigcirc_{r=l+1}^{L}\Lambda_{T(r)}$ and
$\Lambda^{q_{i}}\equiv\bigcirc_{r=1}^{l}\Lambda_{T(r)}$. Obviously, we have
$\Lambda=\Lambda^{\bar{q}_{i}}\circ\Lambda^{q_{i}}$. Furthermore, let
$\mathcal{F}(\Lambda^{q_{i}})\subset\\{1,\ldots,N\\}$ indicate the set of
subsystems on which $\Lambda^{q_{i}}$ acts non-trivially and
$\bar{\mathcal{F}}(\Lambda^{q_{i}})=\\{1,\ldots,N\\}\setminus\mathcal{F}(\Lambda^{q_{i}})$
its complement.

To simplify the notation, we now denote the dual effects by
$D=\bigotimes_{i=1}^{N}D_{i}$ and the Pauli strings by
$P=\bigotimes_{i=1}^{N}P_{i}$. We can write

$\displaystyle\mathrm{Tr}[\Lambda(D)P]$
$\displaystyle=\mathrm{Tr}\left[\Lambda^{\bar{q}_{i}}\left(\Lambda^{q_{i}}\left(\bigotimes\limits_{k\in\mathcal{F}(\Lambda^{q_{i}})}D_{k}\right)\otimes\bigotimes\limits_{k^{\prime}\in\bar{\mathcal{F}}(\Lambda^{q_{i}})}D_{k^{\prime}}\right)\bigotimes\limits_{k^{\prime\prime}=1}^{N}P_{k^{\prime\prime}}\right]$
(10)

In order to make the next steps clear and explicit, it is useful to write,
using the fact that $\Lambda^{\bar{q}_{i}}$ acts trivially on subsystem
$q_{i}$,
$\Lambda^{\bar{q}_{i}}(\cdot)=\sum_{a,b}\lambda_{ab}\mathbb{I}_{i}\otimes
B_{a}\cdot\mathbb{I}_{i}\otimes B_{b}^{\dagger}$, where $\\{B_{a}\\}$ is a
basis of $L(\bigotimes_{k\neq i}\mathcal{H}_{k})$. Therefore, we have

$\displaystyle\mathrm{Tr}[\Lambda(D)P]$
$\displaystyle=\sum_{a,b}\lambda_{ab}\mathrm{Tr}\Bigg{[}\left(\mathbb{I}_{i}\otimes
B_{a}\right)\left(\Lambda^{q_{i}}\left(\bigotimes\limits_{k\in\mathcal{F}(\Lambda^{q_{i}})}D_{k}\right)\otimes\bigotimes\limits_{k^{\prime}\in\bar{\mathcal{F}}(\Lambda^{q_{i}})}D_{k^{\prime}}\right)\left(\mathbb{I}_{i}\otimes
B_{b}^{\dagger}\right)\left(P_{i}\otimes\mathbb{I}\right)\left(\mathbb{I}_{i}\otimes\bigotimes\limits_{k^{\prime\prime}\neq
i}P_{k^{\prime\prime}}\right)\Bigg{]}$ (11)
$\displaystyle=\sum_{a,b}\lambda_{ab}\mathrm{Tr}\left[B_{a}\left(\mathrm{Tr}_{i}\left[\left(P_{i}\otimes\mathbb{I}\right)\Lambda^{q_{i}}\left(\bigotimes\limits_{k\in\mathcal{F}(\Lambda^{q_{i}})}D_{k}\right)\right]\otimes\bigotimes\limits_{k^{\prime}\in\bar{\mathcal{F}}(\Lambda^{q_{i}})}D_{k^{\prime}}\right)B_{b}^{\dagger}\bigotimes\limits_{k^{\prime\prime}\neq
i}P_{k^{\prime\prime}}\right]$

Let us now redefine
$\Lambda^{\bar{q}_{i}}(\cdot)=\sum_{a,b}\lambda_{ab}B_{a}\cdot
B_{b}^{\dagger}$ and denote
$R_{\mathcal{F}(\Lambda^{q_{i}})\setminus\\{i\\}}=\mathrm{Tr}_{i}\left[\left(P_{i}\otimes\mathbb{I}\right)\Lambda^{q_{i}}\left(\bigotimes\limits_{k\in\mathcal{F}(\Lambda^{q_{i}})}D_{k}\right)\right]$.
The above expression now reads

$\mathrm{Tr}[\Lambda(D)P]=\mathrm{Tr}\left[\Lambda^{\bar{q}_{i}}\left(R_{\mathcal{F}(\Lambda^{q_{i}})\setminus\\{i\\}}\otimes\bigotimes\limits_{k\in\bar{\mathcal{F}}(\Lambda^{q_{i}})}D_{k}\right)\bigotimes\limits_{k^{\prime}\neq
i}P_{k^{\prime}}\right].$ (12)

Writing $\Lambda^{\bar{q}_{i}}=\Lambda^{\bar{q}_{j}}\circ\Lambda^{q_{j}}$
where, this time, $\Lambda^{\bar{q}_{j}}$ acts trivially on subsystem $q_{j}$,
we can iterate the process by applying $\Lambda^{q_{j}}$ to
$R_{\mathcal{F}(\Lambda^{q_{i}})\setminus\\{i\\}}\otimes\bigotimes_{k\in\bar{\mathcal{F}}(\Lambda^{q_{i}})}D_{k}$
(or, rather, to the relevant subsystems
$(\mathcal{F}(\Lambda^{q_{i}})\setminus\\{i\\})\cup\mathcal{F}(\Lambda^{q_{j}})$
if
$(\mathcal{F}(\Lambda^{q_{i}})\setminus\\{i\\})\cap\mathcal{F}(\Lambda^{q_{j}})\neq\emptyset$
or to $\mathcal{F}(\Lambda^{q_{j}})$ otherwise) and taking the corresponding
partial trace over $q_{j}$.

Importantly, the highest-dimensional operators that must be computed in the
above calculation are the terms such as
$\left(P_{i}\otimes\mathbb{I}\right)\Lambda^{q_{i}}\left(\bigotimes_{k\in\mathcal{F}(\Lambda^{q_{i}})}D_{k}\right)$.
Thus, the efficiency of the algorithm relies on the appropriate choice of
integer $l$ at each step (typically, a value that minimises the number of
subsystems on which the corresponding $\Lambda^{q_{i}}$ acts, that is,
$|\mathcal{F}(\Lambda^{q_{i}})|$, is desirable). Similarly, in order to make
sure that the procedure only involves the explicit calculation of operators of
small enough dimension, a proper choice of map structure and permutation
$T\in\mathcal{T}$ are necessary as well. Indeed, notice that the residual
operators $R$ live in the Hilbert space of the subsystems in the union of
intermediate relevant spaces (minus the traced-out parties). Hence, it is in
general a good idea to design the map structure so that only one party needs
to be added in every iteration, so the partial trace balances the total count.
The text presents an example in which all the elements have been chosen in
such a way that the computation is efficient for any number of qubits $N$.

## References

  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Science 370, 1460 (2020).
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, A. E. Lita, T. Gerrits, S. W. Nam, V. D. Vaidya, M. Menotti, I. Dhand, Z. Vernon, N. Quesada, and J. Lavoie, Nature 606, 75 (2022).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Journal of the Physical Society of Japan 90, 032001 (2021), publisher: The Physical Society of Japan.
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Reviews of Modern Physics 94, 015004 (2022), publisher: American Physical Society.
  * Urbanek  _et al._ [2021] M. Urbanek, B. Nachman, V. R. Pascuzzi, A. He, C. W. Bauer, and W. A. de Jong, Physical Review Letters 127, 270502 (2021).
  * Wallman and Emerson [2016] J. J. Wallman and J. Emerson, Physical Review A 94, 052325 (2016).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Physical Review X 7, 021050 (2017), publisher: American Physical Society.
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Physical Review Letters 119, 180509 (2017), publisher: American Physical Society.
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Physical Review X 8, 031027 (2018), publisher: American Physical Society.
  * Smart and Mazziotti [2019] S. E. Smart and D. A. Mazziotti, Physical Review A 100, 022517 (2019), publisher: American Physical Society.
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019), number: 7749 Publisher: Nature Publishing Group.
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Nature Physics 16, 1050 (2020).
  * Suchsland  _et al._ [2021] P. Suchsland, F. Tacchino, M. H. Fischer, T. Neupert, P. K. Barkoutsos, and I. Tavernelli, Quantum 5, 492 (2021), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Wiersema  _et al._ [2022] R. Wiersema, L. Guerini, J. F. Carrasquilla, and L. Aolita, “Circuit connectivity boosts by quantum-classical-quantum interfaces,”  (2022).
  * Ravi  _et al._ [2022] G. S. Ravi, P. Gokhale, Y. Ding, W. M. Kirby, K. N. Smith, J. M. Baker, P. J. Love, H. Hoffmann, K. R. Brown, and F. T. Chong,  _CAFQA: Clifford Ansatz For Quantum Accuracy_, Tech. Rep. arXiv:2202.12924 (arXiv, 2022) arXiv:2202.12924 [quant-ph] type: article.
  * Rakyta and Zimborás [2022] P. Rakyta and Z. Zimborás,  _Efficient quantum gate decomposition via adaptive circuit compression_, Tech. Rep. arXiv:2203.04426 (arXiv, 2022) arXiv:2203.04426 [quant-ph] type: article.
  * McClean  _et al._ [2014] J. R. McClean, R. Babbush, P. J. Love, and A. Aspuru-Guzik, The Journal of Physical Chemistry Letters 5, 4368 (2014), publisher: American Chemical Society.
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Physical Review A 92, 042303 (2015), publisher: American Physical Society.
  * Babbush  _et al._ [2018] R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven, and G. K.-L. Chan, Physical Review X 8, 011044 (2018), publisher: American Physical Society.
  * Cai [2020] Z. Cai, Physical Review Applied 14, 014059 (2020), publisher: American Physical Society.
  * Cramer  _et al._ [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Nature Communications 1 (2010), 10.1038/ncomms1147.
  * da Silva  _et al._ [2011] M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Physical Review Letters 107 (2011), 10.1103/physrevlett.107.210404.
  * Torlai  _et al._ [2018] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Nature Physics 14, 447 (2018).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Nature Machine Intelligence 1, 155 (2019).
  * Paini and Kalev [2019] M. Paini and A. Kalev, “An approximate description of quantum states,”  (2019).
  * Morris and Dakić [2020] J. Morris and B. Dakić, _Selective Quantum State Tomography_, Tech. Rep. arXiv:1909.05880 (arXiv, 2020) arXiv:1909.05880 [quant-ph] type: article.
  * Jiang  _et al._ [2020] Z. Jiang, A. Kalev, W. Mruczkiewicz, and H. Neven, Quantum 4, 276 (2020).
  * Huggins  _et al._ [2021] W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, npj Quantum Information 7, 1 (2021), number: 1 Publisher: Nature Publishing Group.
  * García-Pérez  _et al._ [2021] G. García-Pérez, M. A. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, PRX Quantum 2 (2021), 10.1103/prxquantum.2.040342.
  * Morris  _et al._ [2022] J. Morris, V. Saggio, A. Gocanin, and B. Dakić, Advanced Quantum Technologies 5, 2100118 (2022).
  * Acharya  _et al._ [2021] A. Acharya, S. Saha, and A. M. Sengupta,  _Informationally complete POVM-based shadow tomography_, Tech. Rep. arXiv:2105.05992 (arXiv, 2021) arXiv:2105.05992 [quant-ph] type: article.
  * Ariano  _et al._ [2004] G. M. D. Ariano, P. Perinotti, and M. F. Sacchi, Journal of Optics B: Quantum and Semiclassical Optics 6, S487 (2004), publisher: IOP Publishing.
  * Guerini  _et al._ [2021] L. Guerini, R. Wiersema, J. F. Carrasquilla, and L. Aolita,  _Quasiprobabilistic state-overlap estimator for NISQ devices_, Tech. Rep. arXiv:2112.11618 (arXiv, 2021) arXiv:2112.11618 [quant-ph] type: article.
  * Son and Vedral [2009] W. Son and V. Vedral, Open Systems & Information Dynamics 16, 281 (2009), publisher: World Scientific Publishing Co.
  * Sokolov  _et al._ [2022] B. Sokolov, M. A. C. Rossi, G. García-Pérez, and S. Maniscalco, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20200421 (2022), publisher: Royal Society.
  * Cruz  _et al._ [2019] D. Cruz, R. Fournier, F. Gremion, A. Jeannerot, K. Komagata, T. Tosic, J. Thiesbrummel, C. L. Chan, N. Macris, M.-A. Dupertuis, and C. Javerzac-Galy, Advanced Quantum Technologies 2, 1900015 (2019), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/qute.201900015.
  * Ran [2020] S.-J. Ran, Physical Review A 101 (2020), 10.1103/physreva.101.032310.
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Reviews of Modern Physics 93 (2021), 10.1103/revmodphys.93.045003.
  * Note [1] For instance, for a layered structure like the one in Fig. 2, this can be done in the following way. Let $\Lambda$ be a map such that $\varrho^{\prime}=\Lambda(\outerproduct{0}{0}^{\otimes N})$, and $\mathcal{Z}$ be the single-qubit CPTP map defined by $\mathcal{Z}(\sigma)=\outerproduct{0}{0}$. Now, consider the map $\Lambda^{\prime}$ with the same structure as $\Lambda$ and in which $\Lambda_{ij}^{{}^{\prime}l}=\Lambda_{ij}^{l}$ for $l>1$. For the first layer $l=1$, let $\Lambda_{12}^{{}^{\prime}1}=\Lambda_{12}^{1}\circ(\mathcal{Z}_{1}\otimes\mathcal{Z}_{2})$ and $\Lambda_{ij}^{{}^{\prime}1}=\Lambda_{ij}^{1}\circ(\PazoBB{1}_{i}\otimes\mathcal{Z}_{j})$ for all other maps in the layer. The map $\Lambda^{\prime}$ maps any input state to $\varrho^{\prime}$.


Spanning Tree Matching Decoder
for Quantum Surface Codes

Diego Forlivesi, Student Member, IEEE,
Lorenzo Valentini, Member, IEEE,
and Marco Chiani, Fellow, IEEE
The authors are with the Department of Electrical, Electronic, and Information Engineering “Guglielmo Marconi” and CNIT/WiLab, University of Bologna, V.le Risorgimento 2, 40136 Bologna, Italy. E-mail: {diego.forlivesi2,lorenzo.valentini13<EMAIL_ADDRESS>Work funded in part by the European Union - Next Generation EU, PNRR project PRIN n. 2022JES5S2.

Forlivesi, Valentini, Chiani: Spanning Tree Matching Decoder for Quantum Surface Codes

We introduce the STM decoder for surface codes, which guarantees the error correction capability up to the code's designed distance by first employing an instance of the MST on a subset of ancilla qubits within the lattice. Then, a perfect matching graph is simply obtained, by selecting the edges more likely to be faulty. A comparative analysis reveals that the STM decoder, at the cost of a slight performance degradation, provides a substantial advantage in decoding time compared to the MWPM decoder. Finally, we propose an even more simplified and faster algorithm, the RFire decoder, designed for scenarios where decoding speed is a critical requirement.

Quantum Error Correcting Codes, Quantum Communications, Quantum Computing, Surface Codes

§ INTRODUCTION

Leveraging the unique characteristics of quantum mechanics has considerably expanded the possibilities within the realm of information management, encompassing various domains such as sensing, processing, and communication <cit.>. The key hurdle in constructing a quantum computer is the inevitable presence of errors that, if not addressed, rapidly degrade quantum information. Therefore, error correction is crucial for meaningful quantum computation <cit.>. Surface codes are considered central to the architecture for the first generation of quantum computers, thanks to their high error thresholds, planar structure, and locality <cit.>. The MWPM decoder is presently the most widely used decoder for surface codes <cit.>. This decoder results in large threshold error rates, but in practice, its polynomial time complexity introduces a latency that can be a bottleneck for fault-tolerant quantum computing architectures <cit.>. Considerable efforts have been invested in optimizing the performance of this decoder <cit.>. Sub-optimal algorithms based on the UF decoder have also been proposed, achieving almost linear time in code length <cit.>. Additionally, there is ongoing research into neural network-based solutions <cit.>.

In this letter, we present a fast decoding technique tailored for surface codes, called the STM decoder. The STM algorithm involves implementing an instance of the MST on a subset of the ancilla qubits within the lattice. This is followed by a simple and fast construction of a perfect matching graph, resulting in the estimated error pattern. Finally, we propose an even more simplified and faster algorithm, the RFire decoder, tailored for situations where decoding speed is of paramount importance. We conduct a comparative analysis of the performance, considering logical error rates and execution times, between the STM, RFire, and MWPM decoders. We show that, at the price of some performance degradation, the proposed algorithms offer significant advantages in terms of decoding time.

§ PRELIMINARIES AND BACKGROUND

§.§ Quantum Stabilizer Error-Correcting Codes

The Pauli operators are denoted as $\PauliX, \PauliY$, and $\PauliZ$.
We indicate with $[[n,k,d]]$ a  QECC encoding $k$ logical qubits $\ket{\varphi}$ into a codeword of $n$ data qubits $\ket{\psi}$, with minimum distance $d$. The code allows the correction of all patterns with up to $t = \lfloor(d-1)/2 \rfloor$ data qubit errors.

Employing the stabilizer formalism, each code is characterized by $n-k$ independent and commuting operators $\M{G}_i \in \mathcal{G}_n$, termed stabilizer generators or simply generators, with $\mathcal{G}_n$ being the Pauli group on $n$ qubits <cit.>.

The subgroup of $\mathcal{G}_n$ generated by all combinations of the $\M{G}_i \in \mathcal{G}_n$ is called stabilizer and indicated as $\mathcal{S}$.
The code $\mathcal{C}$ is the set of quantum states $\ket{\psi}$ stabilized by $\mathcal{S}$, i.e., satisfying
$\M{S}\ket{\psi}=\ket{\psi} \, \forall \M{S} \in \mathcal{S}$, or, equivalently, $\M{G}_i \ket{\psi}=\ket{\psi},\, i=1, 2, \ldots, n-k$.
The operators that commute with the stabilizer group but are not part of it are called logical operators.

The generators specify measurements on quantum codewords that do not disturb the original quantum state.
These measurements are carried out using extra qubits, named ancillas.
In fact, assume an error $\M{E}\in \mathcal{G}_n$ affects a codeword, so that the state becomes  $\M{E}\ket{\psi}$. It is possible to extract a binary sequence $\V{s}$ (also referred to as error syndrome) where the $i$-th entry $s_i$ is zero if $\M{G}_i$ commutes with $\M{E}$, while $s_i=1$ if $\M{G}_i$ anticommutes with it.
This enable the possibility to perform quantum error correction through error syndrome decoding using as input the binary sequence $\V{s}$.
We will refer to ancillas measuring $s_i = 1$ as defects.

Among stabilizer codes we find the surface codes. These codes arrange qubits on a planar sheet <cit.>.
One of the advantages of such an arrangement, is that it requires only nearest-neighbor interactions between qubits.
Also, this structure enables the possibility to perform a single round of stabilizer measurements with parallel operations <cit.>. Logical operators can be easily identified on the surface codes: $\M{Z}_L$ ($\M{X}_L$) operator consists of a tensor product of $\M{Z}$'s ($\M{X}$'s) crossing horizontally (vertically) the lattice <cit.>.

In the case of a QECC over the depolarizing channel with data qubit error probability $p \ll 1$, we can approximate the logical error rate
as <cit.>

\begin{align}
\label{eq:error_probWithBetaApprox}
&\approx \left(1-\beta_{t+1}\right) \binom{n}{t+1}p^{t+1} %+ (1-\beta_{t+2}) \binom{n}{t+2}p^{t+2}
\end{align}

where $\beta_{t+1}$

is the fraction of errors of weight $t+1$ that the decoder is able to correct.

§.§ Minimum Weight Perfect Matching

A key characteristic of surface codes is the availability of a minimum weight decoder, known as the MWPM. This syndrome decoder constructs a graph in which vertices represent defects, and edges are assigned weights depending on the error probability of the qubits separating the pair. For instance, in case of i.i.d. data qubit errors, these weights correspond to the number of qubits between the pair.
In general, such a weighting procedure is performed using the Dijkstra algorithm to connect the defects on the graph representing the full lattice <cit.>.
By suitably assigning distances on the lattice, Dijkstra's algorithm allows for the consideration of different systems and qubit error statistics.
For surface and rotated surface structures <cit.>, considering i.i.d. data qubit errors, this step can be greatly simplified adopting the Manhattan distance between defects.
Then, the MWPM proceeds by matching the graph to estimate the data qubit errors. A matching of a graph is a set of edges such that no two edges in the matching share a common vertex. A perfect matching is a matching that includes all vertices in the graph <cit.>.

For topological quantum codes without boundaries (toric codes),
we can directly search for perfect matches using the MWPM algorithm on the resulting graph <cit.>.
For surface codes, which have boundaries, it is necessary to add ghost defects (also known as ghost ancillas <cit.>) before applying Dijkstra's algorithm.
This is because errors occurring along a boundary excite only one ancilla, leading to an odd number of defects which impedes running the MWPM.
Therefore, for each defect, a corresponding ghost defect is considered by the decoder.
In the final graph, these defects are also connected between themselves with zero distance in order to preserve the error correction capability of the code <cit.>.

The decoder based on the MWPM is able to guarantee the error correction capability $t$.
In addition, it is also able to correct a significant amount of error patterns of weight greater than $t$.
As an example, for the $[[9, 1, 3]]$ rotated surface code over the depolarizing channel, we have $\beta_2 = 0.5$, i.e., the decoder is able to correct $50\%$ of the error patterns of weight two <cit.>.
As implemented in the LEMON, the MWPM algorithm complexity is $O(N M \, \log N)$, where $N$ and $M$ represent the number of vertices and edges in the graph, respectively <cit.>.

Specifically, for a toric code we have $N = n_\mathrm{d}$ and $M = \binom{n_\mathrm{d}}{2}$, while for a surface code $N = 2\,n_\mathrm{d}$ and $M = \binom{2\,n_\mathrm{d}}{2}$ due to ghost ancillas, where $n_\mathrm{d}$ denotes the number of defects.

Since planar architectures are favorable from an implementation point of view, we will focus on surface codes in the following.
The STM decoder for toric codes can be easily constructed by straightforward variations on the algorithm described below.

§ SPANNING TREE MATCHING DECODER

In this Section, we introduce the STM decoder for surface codes.
First, from (<ref>) we observe that, for a code with given $n$ and $t$, the error correction capability is determined by the fraction of errors of weight $ w = t + 1$ that the decoder is able to correct.

Since the surface codes belong to the class of CSS codes, for the sake of simplicity we will refer to the lattice where the sites constitute the $\M{X}$ generators.
The same reasoning can be applied to the dual lattice.

The STM decoder consists of three phases: MST evaluation, tree matching procedure, and error correction.

§.§ Minimum Spanning Tree Phase

We first construct the complete graph $\mathcal{G} = (n_\mathrm{d},\binom{n_\mathrm{d}}{2})$ connecting all defects, by using Dijkstra or the Manhattan distances.

Note that, since there are no ghost ancillas at this stage, this graph has fewer edges and vertices than the one used by the MWPM.

Then, we execute the MST algorithm on the graph.
This can be achieved with a complexity of $O(M\log N)$, where $M= \binom{n_\mathrm{d}}{2}$ and $N = n_\mathrm{d}$ stand for the edges and the vertices of the graph, respectively <cit.>.
At this stage, we construct two MST starting from the obtained one.
If the number of defects is already even, one of the output trees is the original one, while the other is derived by introducing a ghost defect on both the left and right sides.
If the number of defects is odd, an output tree is obtained by adding a ghost defect to the left, while the alternative output tree is built by adding a ghost defect to the right. Then, the added ghost defects are connected to the nearest defect in the lattice. An example of this phase is reported in Fig. <ref>(b).
In case of multiple defects with a ghost defect at the same distance, we select the one with the highest minimum distance to other non-ghost defects.

Spanning tree matching decoder with a $[[85,1,7]]$ surface code. a) Three $\M{Z}$ channel errors occur on the lattice. Exited ancillas are depicted in red. b) Two alternative MST obtained with the nearest ghost ancilla to the left (above) and to the right (below) boundary, respectively. c) Resulting $\mathcal{E}$ from the tree matching procedure.

§.§ Tree Matching Phase

During this step, we match the trees to obtain the estimated error patterns. For a tree $\mathcal{T}$ let us define as $\mathcal{E}$ its perfect matching graph.

We use $\text{adj}(v, \mathcal{G})$ to indicate the set of vertices adjacent to the vertex $v$ in the graph $\mathcal{G}$, and $\text{deg}(v, \mathcal{G})$ for the degree of $v$ in the graph $\mathcal{G}$. $\mathcal{B}(\mathcal{T})=\{v \in \mathcal{T} \, | \,\text{deg}(v, \mathcal{T})=1\}$ is the set of boundary vertices in  $\mathcal{T}$.

A simple algorithm for obtaining the perfect matching of an MST consists of the following steps.

1. For each $b \in \mathcal{B}(\mathcal{T})$

Let $a=\text{adj}(b, \mathcal{T})$;

If $\text{deg}(a, \mathcal{T}) = 2$: call the edges $e_1 = (a, b)$ and $e_2 = (a, \, \text{adj}(a, \mathcal{T}) \setminus \{b\})$. Add $e_1$ to $\mathcal{E}$, and remove the subgraph $(\{ a, b\}, \{e_1, e_2\})$ from $\mathcal{T}$.

2. For each $b \in \mathcal{B}(\mathcal{T})$

Let $a=\text{adj}(b, \mathcal{T})$;

If $\text{deg}(a, \mathcal{T}) = 3$: call $\{v_1,v_2\}= \text{adj}(a, \mathcal{T}) \setminus \{b\}$, and the edges $e_1 = (a, b)$, $e_2 = (a, v_1)$, $e_3 = (a, v_2)$. Add the edge $e_1$ to $\mathcal{E}$ and remove from $\mathcal{T}$ the subgraph $(\{ a, b\}, \{e_1, e_2, e_3\})$. Insert the edge $e = (v_1, v_2)$ with weight $w(e) = w(e_2) + w(e_3)$ in $\mathcal{T}$. Return to step 1.

3. For each $b \in \mathcal{B}(\mathcal{T})$

Let $a=\text{adj}(b, \mathcal{T})$;

If $\text{deg}(a, \mathcal{T}) = 4$: remove from $\mathcal{T}$ the edge $e = (a, v)$ where $v$ is the sole vertex in $\text{adj}(a, \mathcal{T})$ with $\text{deg}(v, \mathcal{T}) > 1$. Return to step 2.

Note that, considering i.i.d. data qubit errors, a vertex $v$ in the MST cannot have $\text{deg}(v, \mathcal{T}) > 4$.

The algorithm terminates when the graph $\mathcal{T}$ becomes empty, and gives a set of edges $\mathcal{E}$ representing a valid solution to the error correction problem. This procedure is carried out for both the MST obtained in Section <ref>, and the two solutions are indicated as $\mathcal{E}_1$ and $\mathcal{E}_2$.
Examples of the tree matching procedure described above are depicted in Fig. <ref>(c) and in Fig. <ref>(c).

§.§ Error Correction

After the previous phase we have two possible sets of faulty qubits, $\mathcal{E}_1$, $\mathcal{E}_2$, with a total number of data qubit errors $w_1$ and $w_2$, respectively. An instance of STM decoding, involving both MST, is illustrated in Fig. <ref>.
If $w_1 \le t+1$ or $w_2 \le t+1$, then we have found the correction operator[Note that, if $w_1 \leq t$, the matching of the second tree can be avoided to save time.].

If, instead, both $w_1$ and $w_2$ have weight $> t+1$, a simple processing considering both $\mathcal{E}_1$ and $\mathcal{E}_2$ can be applied, which guarantees the correction if the channel errors has weight $\le t$ (see  Section <ref>).
The basic idea is choosing the solution with the smallest number of horizontally traversed edges.

This is motivated by the fact that logical operators traverse the lattice from one side to the other one.
Hence, the correction with more operators along the horizontal dimension will more likely cause a logical operator (i.e., an undetected error).

A column is the set of horizontal edges aligned in the vertical direction of the lattice, as shown in Fig. <ref>(a).

In this way, we can enumerate the columns from left to right ranging from $1$ to $d$.

Then, we can define two vectors $\V{c}_{i}$ with entries $c_{i,j}$, where $i=1, 2$ and $j=1, \dots, d$, representing the cardinality of the intersection between $\mathcal{E}_i$ and the $j$-th column.
Since a solution $\mathcal{E}_i$ with two edges in the $j$-th column is equivalent, by adding a stabilizer, to another solution without any edge in the $j$-th column,

we also define as $\V{u}_{i}$ a vector with entries $u_{i,j} = c_{i,j} \mod 2$.
As a consequence, the function

\begin{align}
\label{eq:metricCol}
f(\mathcal{E}_i) = \sum_{j=1}^{d} u_{i,j}
\end{align}

can be used as a metric to quantify how much correcting $\mathcal{E}_i$ is likely to cause a logical operator.
Thus, the final solution is that minimizing (<ref>).
For instance, in Fig. <ref>(c) there are three edges in solution $\mathcal{E}_1$, with weights $w = 1$, $w = 2$, and $w = 5$. They consist of one, two, and three horizontal qubits, respectively. However, the horizontal qubit of the edge with weight $w = 1$ belongs to the same column as the edge with weight $w = 5$. Hence, $f(\mathcal{E}_i) = 4$.

Referring to the lattice in Fig. <ref>a, let us designate the horizontal qubits with numbers ranging from $0$ to $d^2 - 1$, progressing from left to right and top to bottom. Starting from $\mathcal{E}$, we can derive a vector $\mathbf{e}$ of length $d^2 - 1$ where a value of $1$ indicates the presence of the corresponding qubit in the solution, and $0$ signifies its absence. Subsequently, we define a vector $\mathbf{v}$ of length $d$ as

\begin{align}
v_i = \sum_{l = 0}^{d^2 - 1} e_l \, [\,i &\equiv l \!\!\! \mod (d^2 - 1)\,]
\end{align}

where $[\cdot]$ equals 1 when the condition within the brackets is true and 0 otherwise. Hence, we introduce a vector $\mathbf{u}$ of length $d$ where

\begin{align}
u_i =  v_i \!\!\! \mod 2.
\end{align}

Finally, we define a function $f(\cdot)$, which takes a vector of length $d$ as input and yields the sum of its elements as output. The STM decoder, in the event that no solution satisfies $ w(\mathcal{E}) \leq t + 1$, selects as

\begin{align}
\label{eq:min}
\mathcal{E} =  \argmin_{\mathcal{E}_i; \,\, i = 1,2} f(\mathbf{u}_i)
\end{align}

where $\mathcal{E}_1$ and $\mathcal{E}_2$ refer to the first and to the second solution obtained through the tree pruning.

Spanning tree matching decoder with a $[[85,1,7]]$ surface code. Both matched spanning trees have $w > t + 1$. Hence, the error correction is performed according to (<ref>).

§.§ Distance-preserving decoding

In this section we show that (<ref>), in the case of a $[[n,k,d=2t +~1]]$ surface code, assures the correction of channel errors with weight $w \leq t$.

Given a surface code, let us call $\mathcal{C}$ an arbitrary Pauli $\M{Z}$ error chain connecting two sites defects $v_1$ and $v_2$. Then, the intersection of $\mathcal{C}$ with any column between $v_1$ and $v_2$ has cardinality $1 \pmod 2$. The intersection with any of the other columns has cardinality $0 \pmod 2$.

The minimal-weight chain linking $v_1$ and $v_2$ clearly satisfies the theorem. Every other chain can be obtained from this one, through the application of plaquette generators. Since each plaquette has two edges in the same column, it does not affect the modulo 2 counting.

Given any error pattern over an $[[n,k,d]]$ surface code, resulting in an even number of site defects (including also ghost defects), all possible perfect matchings $\mathcal{M}$ have the same metric $f(\mathcal{M})$.

We label the horizontal position of the sites in the lattice, from left to right, as $0, 1, \ldots, d+1$, with $0$ for the left ghost defect and $d+1$ for the right ghost defect.
Vertically aligned sites share the same label.
Let us firstly examine the case where the lattice has four site defects.
Reordering them we obtain four indexes $0 \le j_1 \le j_2 \le j_3 \le j_4 \le d+1$, representing the site positions.
In this setup we can have three possible perfect matchings $\mathcal{M}_i$, with $i=1,2,3$.
Considering  $\mathcal{M}_1$ as the one connecting site $1$ with $2$, and site $3$ with $4$, by application of Lemma <ref>, we have that $f(\mathcal{M}_1) = (j_2 - j_1) + (j_4 - j_3)$. It is easy to check that the other two matchings have the same metric.

Finally, we observe that each perfect matching $\mathcal{M}$ in any graph $\mathcal{G}$ with an even number of defects can be obtained by iteratevely removing two edges and reconnecting them as needed.
Hence, the claim follows since this operation does not affect the metric $f(\mathcal{M})$.

Each perfect matching has the same metric (<ref>) as the MWPM.

It follows that (<ref>) can be interpreted as the minimum number of traversed columns for all perfect matchings of the same graph.

We have seen in Section <ref> that, starting from a complete graph $\mathcal{G}$ formed by connecting all defects, we obtain two distinct graphs $\mathcal{G}_1$ and $\mathcal{G}_2$ by adding ghost ancillas. Due to the even number of defects in these graphs, we can compute two perfect matchings, denoted as $\mathcal{E}_1$ and $\mathcal{E}_2$. It is noteworthy that applying a logical operator to one matching we get the other. Consequently, one of the solutions will correct the error (call it $\mathcal{E}_c$), while the alternative, $\mathcal{E}_a$, will produce a logical error.

Let us consider an $[[n,k,d]]$ surface code and an error pattern of weight $w \le t$. A perfect matching leading to the correct solution, $\mathcal{E}_c$,
has $f(\mathcal{E}_c) < f(\mathcal{E}_a)$, with $\mathcal{E}_a$ being any perfect matching on the alternative graph.

The number of columns in the lattice is $d$.
Since we are considering error patterns with weight $w \le t$, we have that any perfect matching $\mathcal{E}_\mathrm{c}$ representing the correct solution, due to Lemma <ref>, satisfies $f(\mathcal{E}_\mathrm{c}) \le t$.
Regarding the alternative solution $\mathcal{E}_\mathrm{a}$, we have that $f(\mathcal{E}_\mathrm{a}) = d - f(\mathcal{E}_\mathrm{c})$, since they differ by a logical operator.
This leads to $f(\mathcal{E}_\mathrm{a}) \ge t+1$ for $d$ odd and $f(\mathcal{E}_\mathrm{a}) \ge t+2$ for $d$ even,  proving the statement.

The previous theorem states that a surface decoder chosing a matching with the minimum traversed column metric (<ref>) preserves the error correction capability of the code.

The same theorem applies to rotated surface codes as well, owing to their similar lattice structure.

§ RAPID-FIRE DECODER

We show now how to design an even faster decoder that ensures the correction for all errors with weights $w \leq t$.
Due to Corollary <ref>, each matching is equivalent to the minimum one in terms of (<ref>).
Hence, it is possible to compute $\mathcal{E}_i$ with $i=1,2$ without evaluating and matching the MST.
Firstly, we compute the complete graph $\mathcal{G}$ on the defects.
Then, we construct two graphs $\mathcal{S}_i$ by adding to $\mathcal{G}$ ghost ancillas on boundaries with the same strategy described in Sections <ref>.

For each $\mathcal{S}_i$, we iteratively pair each defect with its closest one in a greedy fashion, we update the solution $\mathcal{E}_i$, and we remove both the defects from $\mathcal{S}_i$.
This process results in two potential sets of faulty qubits, $\mathcal{E}_1$ and $\mathcal{E}_2$.
Finally, we determine the error correction operator to adopt, following the procedure outlined in Section <ref>.
In this way, we guarantee the correction of all errors of weight up to $w = t$ by Theorem <ref>.
We call this the RFire decoder.

§ NUMERICAL RESULTS

name = plot,
scale only axis,
xlabel style=font=, font = ,
xlabel= $p$,
tick label style=black, semithick, font=,
ylabel style=font=, font = ,
ylabel= $p_L$,
axis background/.style=fill=white,
legend style=legend cell align=left, align=left, draw=white!15!black
[color=brightRed, line width=1.5pt, mark=o, mark size = 2.5pt, mark repeat = 2, mark phase = 2 mark options=solid, brightRed]
table[row sep=crcr]
0.001	2e-05

0.005	0.00049

0.008	0.00118

0.01	0.00184

0.05	0.03946

[color=brightRed, line width=1.5pt, dashed]

table[row sep=crcr]
0.001	1.64e-05

0.005	0.00039113

0.008	0.000997308

0.01	0.001562046

0.05	0.0352298

[color=brightRed, line width=1.5pt, mark=square, mark size = 2.5pt, mark options=solid, brightRed]
table[row sep=crcr]
0.003	4.94286e-06

0.005	2.15675e-05

0.008	9.10919e-05

0.01	0.000172973

0.05	0.026

[color=brightRed, line width=1.5pt, mark=diamond, mark size = 4pt, mark options=solid, brightRed]
table[row sep=crcr]
0.005   2.1e-6

0.008	2e-05

0.02	0.001

0.05	0.04

[color=brightBlue, line width=1.5pt, mark=o, mark size = 2.5pt, mark repeat = 2, mark options=solid, brightBlue]
table[row sep=crcr]
0.001	2e-05

0.005	0.00046

0.008	0.00114

0.01	0.0017

0.05	0.03819

[color=brightBlue, line width=1.5pt, dashed]
table[row sep=crcr]

0.001	1.54e-05

0.005	0.0003926

0.008	0.0010174

0.01	0.0015618

0.05	0.0348174

[color= graphYellow, line width=1.5pt, dashed]
table[row sep=crcr]
0.001	1.7e-05

0.005	0.000427084

0.008	0.001063158

0.01	0.001675

0.05	0.0357

[color=brightBlue, line width=1.5pt, mark=square, mark size = 2.5pt, mark options=solid, brightBlue]
table[row sep=crcr]
0.003	3.8e-06

0.005	1.8e-05

0.008	7.1e-05

0.01	0.000124

0.05	0.01565

[color=brightBlue, line width=1.5pt, mark=diamond, mark size = 4pt, mark options=solid, brightBlue]
table[row sep=crcr]
0.008	3e-06

0.02	0.00015

0.05	0.0063

[color=graphYellow, line width=1.5pt, mark=o, mark size = 2.5pt, mark options=solid, graphYellow]
table[row sep=crcr]
0.001	2.2e-05

0.005	0.0005025

0.008	0.00127778

0.01	0.00190909

0.05	0.0406

[color=graphYellow, line width=1.5pt, mark=square, mark size = 2.5 pt,  mark options=solid, graphYellow]
table[row sep=crcr]
0.003	0.62e-05

0.005	2.7e-05

0.008	0.000108

0.01	0.000205

0.05	0.0312

[color=graphYellow, line width=1.5pt, mark=diamond, mark size = 4 pt, mark options=solid, graphYellow]
table[row sep=crcr]
0.005   5.73298e-06

0.008	3.41702e-05

0.02	0.00138

0.05	0.0462

[draw,fill=white,inner sep=1.5pt, below right=0.5em] at (plot.north west)

RFire      STM     MWPM

$[[9,1,3]]$      <ref>     <ref>     <ref>

$[[13,1,3]]$      <ref>     <ref>     <ref>

$[[41,1,5]]$     <ref>     <ref>     <ref>

$[[85,1,7]]$       <ref>     <ref>     <ref>

Logical error probability vs. physical error probability of the channel. Rotated $[[9,1,3]]$, and standard $[[13,1,3]]$, $[[41,1,5]]$, and $[[85,1,7]]$ surface codes over depolarizing channel.

AVERAGE EXECUTION TIMES $[\mu s]$

$n_\mathrm{d}=1$     $n_\mathrm{d}=2$     $n_\mathrm{d}=3$     $n_\mathrm{d}=4$     $n_\mathrm{d}=5$     $n_\mathrm{d}=6$     $n_\mathrm{d}=7$     $n_\mathrm{d}=8$

RFire          $0.0008$    $0.0015$    $0.0021$     $0.0025$

STM     $[[13,1,3]]$     $0.0025$    $0.0031$    $0.0051$     $0.0123$

MWPM          $17.329$    $30.618$    $62.254$     $63.653$

RFire          $0.0019$     $0.0044$     $0.0051$    $0.0081$     $0.0085$     $0.0108$     $0.0110$    $0.0341$

STM     $[[85,1,7]]$     $0.0031$     $0.0047$     $0.0091$    $0.0121$     $0.0695$     $0.2071$     $1.2839$    $1.9302$

MWPM          $19.156$     $32.247$     $55.864$    $67.591$     $100.55$     $115.04$     $157.67$    $171.78$

In this section we compare the performance of surface and rotated surface codes using MWPM, STM, and RFire decoding via Monte Carlo simulations.
These decoders are implemented in C++ and run with an Apple Silicon M2 processor.
In the implementation, we exploit the LEMON C++ library for an efficient MWPM algorithm <cit.>.
To evaluate the complexity of the decoders we measure the average execution time of the processing that starts with the initial complete graph $\mathcal{G}$ (generation of $\mathcal{G}$ not included) and ends returning the solution $\mathcal{E}$.
The evaluation, carried out when varying the lattice size and the number of defects, is reported in Table <ref>.
As expected, the computational time saving is remarkable when STM and RFire are adopted.

Furthermore, in Fig. <ref>, the logical error rate of some surface codes over a depolarizing channel is depicted.

In particular, the performance gap between MWPM and the proposed fast decoders widens progressively as the distance of the code increases. This is due to the fact that, although the STM and the RFire decoders effectively correct all the errors with weights up to $t$, the MWPM decoder corrects a larger fraction of those with weight $t+1$.
In terms of speed, the RFire decoder is the fastest, but it looses in error correction capability.
Note that although we showcased our decoder using standard surface codes for the sake of clarity, its performance remains consistent when applied to rotated surface types, as illustrated in Fig. <ref> for the $[[9,1,3]]$ code.
Among the decoders available in the literature, we concentrate on comparing with the only current fast decoder, which is the UF <cit.>. To accomplish this, we utilized the Qsurface library <cit.>. Specifically, for the $[[85,1,7]]$ surface code and a physical error rate $p_Z = 0.03$, the UF shows a speedup by a factor of approximately $\cross 2$ with respect to MWPM.
In the same settings, the speedup of RFire and MST decoders with respect to MWPM is $\cross 3625$, and $\cross 164$, resulting in the quickest execution times.

§ CONCLUSIONS

We have introduced two fast decoders designed for surface codes and compared their performance with the widely used MWPM decoder. While observing a slight reduction in error correction capability for errors of weight $j \geq t + 1$, our findings highlight a substantial advantage in terms of execution times, with a speedup of a factor $\cross 10.000$, for both the STM and the RFire decoders.