1

## Bayesian mitigation of measurement errors in multi-qubit experiments

F. Cosco, 1, ∗ F. Plastina, 2, 3 and N. Lo Gullo 2, 3, 4, †

Quantum algorithms and software, VTT Technical Research Centre of Finland Ltd, Tietotie 3, 02150 Espoo, Finland 2 Dipartimento di Fisica, Universit` della Calabria, 87036 Arcavacata di Rende (CS), Italy a 3 INFN, Sezione LNF, gruppo collegato di Cosenza 4

ICAR, Consiglio Nazionale delle Ricerche, Italy

We introduce an implementation of Bayesian measurement error mitigation tailored for multiqubit experiments on near-term quantum devices. Our approach leverages complete information from the readout signal, which is available before any binary state assignment of the qubits. We provide a detailed algorithm workflow, from the calibration of detector response functions to the postprocessing of measurement outcomes, offering a computationally efficient solution suitable for the output size typical of current quantum computing devices and quantum algorithms. We benchmark our protocol on actual quantum computers with superconducting qubits, where the readout signal encodes the measurement information in the IQ clouds before qubit state assignment. Finally, we compare the performance of our algorithm against other measurement error mitigation methods, such as iterative Bayesian unfolding and noise matrix inversion.

## I. INTRODUCTION

Despite the rapid advancements in quantum computing technologies, practical applications of noisy intermediate-scale quantum (NISQ) processing units (QPUs) are still out of reach. While QPU manufacturers keep improving quality metrics [1-3], be it decoherence times, gate and readout fidelities, a lot of effort is going into developing software tools to shrink the gap between acceptable noise levels and utility-scale algorithms [4, 5]. Quantum error mitigation (QEM) embodies such a goal, even though a comprehensive and universally accepted definition for the field is lacking [6]. In general, QEM encompasses all those algorithmic schemes that reduce noise and errors in the expectation value of some figure of merit, such as the average of observables or quantum state probabilities. This is typically achieved by postprocessing the outputs from an ensemble of circuit runs, either the same circuit or a family of properly designed ones [7]. This statistical approach starkly contrasts with quantum error correction, which aims to reduce noise after each shot. Among the most advanced techniques, zero-noise extrapolation [8] with probabilistic error amplification and probabilistic error cancellation [9] have recently been found successful in the largest available QPUs [10-12]. Some specialised techniques aim to target specific sources of errors, such as measurement error mitigation strategies. These strategies are designed to improve the readout step of quantum algorithms execution by modeling the measurement error, which is a result of the faulty physical measurement process. Most readout mitigation techniques rely on detector or measurement tomography [13-16] to obtain the assignment, noise, or confusion matrix. The readout error is then addressed in various ways. One straightforward strategy involves

∗ francesco.cosco@vtt.fi

† nicolino.logullo@unical.it

Figure 1. Pictorial sketch of the heuristic algorithm to perform Bayesian measurement error mitigation in a multiqubit experiment. In the first step (I), we collect the noisy bitstrings counts and reduce the 2 N q outcome subspace to mitigate to the M measured noisy counts for which we have at least one measurement. In step (II), we select two populations to use as variables while keeping all others constant, and apply Bayesian measurement error mitigation to obtain two new mitigated values. The Bayesian update is performed for all possible pairs. Step (II) is repeated until a specified convergence criterion or exit condition is met, resulting in the full set of mitigated values (III).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-6aab06e1-78e6-42f3-bc87-2e3ed5c65c00.jpg)

manipulating the outcome count statistics by multiplying the inverse of the confusion matrix with the measured quantum state probability vector [17]. However, this approach can yield non-physical values, such as negative probabilities, an issue that can be accounted for by incorporating a convex optimizer into the mitigation pipeline [18]. Alternatively, some proposals have developed machine learning-based techniques to perform the assignment process [19-22].

In a recent work [23], we have introduced an efficient and accurate readout measurement scheme for single and multiqubit states based on Bayesian inference, inspired by successful application of Bayesian inference in quantum metrology and sensing [24, 25] or parameter estimation for quantum circuits and states [14, 26-32]. Our method, bypassing the need to invert the noise matrix, leverages the knowledge of the detector faulty response and Bayesian inference to obtain a probability distribution for each qubit state. Moreover, such an approach is not limited to use binary outcomes from the detector, but it can exploit the detector unprocessed analog readout data, before any binary assignment is made. The idea of using analog data has found traction also in quantum error correction, where it has been shown to lead to improved decoders and logical qubits error rates [33, 34]. Still, all measurement error mitigation schemes suffer some form of scalability issue when moving to many-qubit systems, as intuitively the noise matrix size scales exponentially with the number of qubits. Often this issue is tackled by not using the full assignment matrix, or its inverse, and instead working in a subspace defined by the noisy outcomes [35]. In this work, we aim to tackle thoroughly the scaling problem when aiming to mitigate measurement errors with our formerly introduced Bayesian readout and further show how the use of the readout analog data leads to increased performances. Beyond showing a solid proof of principle of the method, we make a full scaling analysis of the algorithm and provide an estimate of the classical computation time needed to post-process the measurement results obtained from currently available superconducting QPUs.

The paper is organised as follows. In section II, we introduce all the necessary theoretical concepts to model measurement errors in quantum computers and how a Bayesian approach can be used to infer the correct outcomes from the noisy ones, provided the knowledge of the detector noise matrix. Furthermore, we introduce a novel algorithm to implement the Bayesian measurement error mitigation scheme when dealing with multi-qubit measurements. We provide a complete step-by-step procedure which is efficient and feasible for current quantum computing devices, and show how to adapt it for different type of available measurement data, being it binary outcome, i. e. '0's or '1's, or raw analog readout detector signal, e.g. IQ-clouds for superconducting QPUs.

In section III, we benchmark our implementation of Bayesian measurement error mitigation running quantum algorithms on an actual quantum computer. We show how the algorithm can be efficiently applied to improve readout accuracy for multi-qubit experiments and how the performances of the mitigation scheme scale with the number of qubit and the number of measurements.

In section IV, we focus on the Bayesian measurement error mitigation scheme applied to the detector analog data, showing how using the IQ-readout signal, instead of a faulty binary outcome for each qubit, can improve dramatically the performances of the mitigation algorithm in a multi-qubit scenario. In order to highlight the merit and qualities of our implementation of Bayesian measurement error mitigation, in section V, we compare its performances against two other measurement error mitigation techniques, i.e. the iterative Bayesian unfolding and the noise matrix inversion with positivity constraints. In section VI, we discuss the theoretical computational cost of our implementation of the algorithm for multi-qubit Bayesian measurement error mitigation and give, as an example, the actual classical CPU time spent to perform the mitigation routine for some of the results shown the previous sections. Throughout this work, we run our quantum algorithms on two superconducting qubit devices, the superconducting qubits QPUs VTTQ20 (Leena from IQM) and IBM127 (Nazca). Finally, we present the conclusions and outline some possible future research directions and applications in Section VII.

## II. THEORETICAL BACKGROUND

Let's start with a mathematical description of a noisy measurement process for a single qubit. Let s ∈ { 0 1 , } be the measurement outcome observed in a single experiment. In the absence of measurement errors, the probability to register either one is fully determined by the quantum state of the system. Formally, a conventional single qubit measurement is modelled with a Positive-Operator Valued Measure (POVM) with two effects { M 0 = | 0 ⟩⟨ 0 | , M 1 = | 1 ⟩⟨ 1 |} . Consequently, the probability of observing a measurement outcome s , when the quantum system is in a state ˆ, can be calculated as ρ Pr( ) = Tr(ˆ s ρM s ) ≡ ρ s .

However, in real a setup, the readout step is noisy and errors emerge in the process of collecting the signal from the detector and in the assignment of the correct qubit state. In practice, we can model the experimental noisy measurement trough a set of noisy effects M noisy i = ∑ j ⟨ | i Λ | j ⟩ M j = ∑ j Λ ij M j , where Λ is the noise matrix and Λ ij is the probability of observing/assigning outcome i when the true outcome is j . Therefore, after repeated measurements, we have direct access to the noisy probabilities

$$\rho _ { s } ^ { \text{noisy} } = \sum _ { i j } \Lambda _ { i j } \rho _ { j } \delta ( i - s ), \text{ \quad \ \ } ( 1 )$$

which are connected to the true qubit state probabilities trough the noise matrix Λ. Thus, the goal of measurement error mitigation is to reconstruct the true qubit state probabilities after measuring ρ noisy s . Assuming that we have access to the noise matrix Λ, from Eq. (1), it is straightforward to attempt to invert the noise matrix and use it to reconstruct the correct qubit state probability vector as ρ = Λ -1 ρ noise . However, this simple measurement error mitigation strategy might result in a probability vector containing negative values due to statistical errors. This issue is often accounted for by employing a convex optimizer which project the corrected vector to a

physical one, even though the issue and the extra step required become more concerning and demanding when dealing with multi-qubit states.

To avoid these issues, we propose here a strategy to bypass the need to invert the noise matrix by noting how the measured outcome ρ noisy s , as introduced in Eq. (1), can be read as the conditional probability of measuring s when the correct qubit state probabilities are ρ 0 and ρ 1 , i.e. ρ noisy s = Pr( s ρ | 0 , ρ 1 ).

Thus, we can reconstruct a probability distribution for the correct probabilities ρ 0 and ρ 1 trough Bayes' theorem as Pr( ρ , ρ 0 1 | s ) ∝ [Λ s 0 ρ 0 +Λ s 1 ρ 1 ] Pr( ρ , ρ 0 1 ). It is convenient to follow an iterative approach and consider individually the N shots values collected after each measurement. In this way, we can update the conditional probability for the quantum state probabilities ρ 0 and ρ 1 shot after shot as

$$\Pr ( \rho _ { 0 }, \rho _ { 1 } | s ^ { ( n ) } ) \in [ \Lambda _ { s 0 } \rho _ { 0 } + \Lambda _ { s 1 } \rho _ { 1 } ] \Pr ( \rho _ { 0 }, \rho _ { 1 } | s ^ { ( n - 1 ) } ), \ ( 2 )$$

where s ( n ) ∈ { 0 1 , } is the bit registered in the n-th experiment, and Pr( ρ , ρ 0 1 | s ( n -1) ) is the prior probability distribution which corresponds to the posterior distribution calculated at the previous step. The idea is to update the prior probability distribution after each iteration step with the newly calculated posterior probability distribution. When assuming no prior knowledge about the system state, we start the iteration with a uniform prior probability distribution Pr( ρ , ρ 0 1 | s (0) ) = const.

The end result of the Bayesian measurement error mitigation scheme is the probability distribution Pr( ρ , ρ 0 1 | s (N shots ) ) which represents the posterior probability distribution for the qubit state probabilities coherent with all the collected measurement data. Finally, this distribution can be used to estimate the qubit state probabilities P 0 / 1 . In our previous work [23] we proposed to use the first moment of the distribution Pr( ρ , ρ 0 1 | s (N shots ) ), namely P 0 = ∫ dρ dρ 0 1 ρ 0 Pr( ρ , ρ 0 1 | s (N shots ) ) and similarly for P 1 . Nevertheless, we found that results are highly improved by choosing the values ρ 0 / 1 for which the probability Pr( ρ , ρ 0 1 | s (N shots ) ) has its maximum , i.e. ( P , P 0 1 ) = arg max Pr( ρ , ρ 0 1 | s (N shots ) ).

ρ ,ρ 0 1

With these ingredients, we can extend the Bayesian measurement error mitigation scheme to the multi-qubit case with relative ease. For a system of N q qubits, the measurement outcome obtained from a single experiment is a N q digits string, s ∈ { 0 1 , } N q . Equivalent to Eq. (1), the conditional probability of observing the multi-qubit string s is related to the correct outcome trough a multiqubit noise matrix as

$$\Pr ( s | \boldsymbol \rho \equiv \boldsymbol \rho _ { 0 }, \dots, \boldsymbol \rho _ { 2 ^ { N } q - 1 } } ) = \sum _ { \text{ij} } \Lambda _ { \text{ij} } \boldsymbol \rho _ { \text{j} } \delta ( \text{i} - s ). \quad ( 3 ) \quad \text{thou} _ { \text{num} }.$$

From the formal expression of Eq. (3), we can see how the number of terms appearing in the noisy conditional

Figure 2. (a): (Top) Measurement clouds in the IQ plane, where each point corresponds to a single shot. The blue cloud represents measurements of the '0' state, i.e., when no operation is performed on the qubit besides the measurement. The orange cloud represents measurements of the '1' state, where a single X-gate is applied to the qubit before measurement. (Bottom) Projected histograms from the IQ plane onto the Q-axis.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-472b7305-ed20-47e6-98bd-5879d363f2e9.jpg)

probability increases exponentially with the number of qubits composing the measured string s . Nonetheless, from a purely mathematical point of view, the Bayesian inference of the correct multi-qubit state probabilities can be achieved in a similar fashion to the single qubit case by employing the iterative Bayesian update in the higher dimension, i.e.

$$\frac { \psi } { 1 } \sum _ { \psi } \Lambda _ { i j } \rho _ { j } \delta ( \mathbf i - s ^ { ( n ) } ) \Pr ( \rho | s ^ { ( n - 1 ) } ), \quad ( 4 )$$

within the parameters space which satisfies the normalisation condition ∑ i ρ i = 1. However, since the computational cost increases exponentially with the number of qubits such a scheme becomes of limited utility beyond a few qubits.

For example, for a system with N q qubits, Eq. (4) requires handling a multi-parameter probability distribution Pr( ρ | s ( n ) ) at each iteration step of the Bayesian update. This distribution involves 2 N q variables ρ i , each sampled within the interval [0 , 1]. If each variable is sampled with n p points, the total probability distribution would necessitate a mesh-grid with n 2 Nq p points. Even though the normalization condition reduces the effective number of independent coordinates the computational cost remains high and out of reach.

Additional complications arise from the noise matrix resulting of the multi-qubit measurement process. Ide-

ally, the reconstruction of the multi-qubit noise matrix appearing in Eq. (4) would require a full detector tomography, which is impractical and unfeasible as the number of necessary measurements scales exponentially with the number of qubits. Consequently, it is common practice to assume uncorrelated readout noise, allowing the multiqubit noise matrix to be constructed as a tensor product of the single-qubit noise matrices as

$$\Lambda = \Lambda ^ { ( q _ { 1 } ) } \Lambda ^ { ( q _ { 2 } ) } \dots \Lambda ^ { ( q _ { N _ { q } } ) }. \quad \ \ ( 5 ) \quad \ e r t { 1 } { 2 }$$

̸

While making the tomographic step easier, even the multi-qubit noise matrix of Eq. (5) does not reduce the computational cost of implementing Eq. (4) and the Bayesian measurement error mitigation remains resource intensive and cannot be employed as it is. For all these reasons, we have developed a heuristic strategy that enables the use of the scheme in the multi-qubit case, as depicted in Fig. 1. The first key component involves reducing the number of populations ρ i included in the iterative procedure. We achieve this by running the mitigation algorithm only on the subspace of multi-qubit state probabilities where noisy outcomes have been measured. In other words, we exclude from the Bayesian cycle any states for which no strings are obtained in the experiment, limiting it to the strings i for which ρ noisy i = 0. Subspace reduction is a common strategy for measurement error mitigation techniques which suffer from the same exponential scaling problem, as discussed in [35, 36].

While parameter reduction is crucial to correct measurement errors in multi-qubit outcomes, it does not make handling the multi-variable probability distribution of Eq. (4) much easier. Therefore, the next ingredient is to effectively reduce the number of variables to be used when handling the probability distribution which undergoes the Bayesian update. We achieve this by demoting all but two probabilities, ρ i and ρ j , to constant estimates R k , and redefining the conditional probability function for measuring a specific multi-qubit string as

̸

$$\Pr ( s | \rho _ { i }, \rho _ { j } ) = \Lambda _ { s i } \rho _ { i } + \Lambda _ { s j } \rho _ { j } + \sum _ { k \neq i, j } \Lambda _ { s k } R _ { k }. \quad ( 6 ) \quad \text{the I} \\ \text{prob}.$$

̸

For a single measurement, Eq. (6) is still the conditional probability of measuring the multi-qubit string s when the correct probabilities for strings i , j , and the remaining k (with i = j = k ) are ρ i , ρ j , and R k respectively. However, only the first two are treated as variables of the probability distribution, while the remaining ones are kept constant. Using this form for the conditional probability, we can apply the Bayesian iterative cycle in the same fashion trough

̸

̸

$$\Pr ( \rho _ { i }, \rho _ { j } | \mathbf s ^ { ( n ) } ) \in & \left [ \Lambda _ { s i } \rho _ { i } + \Lambda _ { s j } \rho _ { j } + \sum _ { k \neq i, j } \Lambda _ { s k } R _ { k } \right ] \quad ( 7 ) \quad \text{mean} \\ & \times \Pr ( \rho _ { i }, \rho _ { j } | \mathbf s ^ { ( n - 1 ) } ),$$

and obtain a final probability distribution for the two variables. The posterior probability distribution Pr( ρ , ρ i j | s (N shots ) ) can then be used to estimate the correct values for the probabilities of measuring i and j , thus obtaining some updated values for R i and R j . After this step, the strategy is to repeat the procedure for a new pair of probabilities until all possible pair combinations have been updated once. This sequential heuristic binary optimization is then repeated until a chosen convergence criterion is achieved. A simple rule of thumb is to stop the algorithm repetitions if the change in the obtained populations does not change with respect to the previous iteration within a chosen accuracy. All the steps of the algorithm are outlined in Algorithm 1.

We have outlined the general theory behind Bayesian measurement error mitigation using a given bit string outcome s . However, this framework can also be employed when we have access to analog detector data before any assignment to 0s and 1s is made. For example, in the case of superconducting qubits, the readout step exploits the dispersive interaction between the superconducting qubit and a readout resonator, which is the actual component of the QPU targeted by the readout pulse [37, 38]. Each measurement registers a point in the IQ-plane for each shot, which is then converted to a binary outcome, often according to its relative position with respect to some threshold value, as shown in Fig. 2. Within the Bayesian mitigation framework, we can leverage the full continuous range of the physical measurement, albeit with the additional cost of calibrating the detector response, i.e. performing detector tomography, accordingly.

This means that if the detector measures the physical value Q , the detector tomography aims to reconstruct a noise matrix with a continuous index Λ Qj ≡ Λ ( j Q ), which can be interpreted as the probability of measuring Q if the qubit was in | j ⟩ . The Bayesian measurement error mitigation algorithm can be applied in the same fashion as for bitstring measurements, noting that in the multi-qubit case we collect the set Q = Q (1) , ..., Q (N ) q instead of the strings s = s (1) , ..., s (N ) q for each shot. Thus, the Bayesian update rule makes use of the conditional probability

̸

$$\underset { \underset { \underset { \text{$1,1} } { \text{ng s} } } { \text{ndi-} } } \, \Pr ( Q | \rho _ { i }, \rho _ { j } ) = \Lambda _ { i } ( Q ) \rho _ { i } + \Lambda _ { j } ( Q ) \rho _ { j } + \sum _ { k \neq i, j } \Lambda _ { k } ( Q ) R _ { k }, \ ( 8 )$$

which is the equivalent of Eq. (6).

In the following sections, we will demonstrate that the detector analog data, when available, contains more information than binary outcomes and its use enhances the effectiveness of the mitigation algorithm with negligible increase in computational time.

Finally, it is worth mentioning that our approach to measurement error mitigation through Bayesian inference shares some similarities with the Iterative Bayesian Unfolding (IBU). IBU is a form of regularized matrix inversion applied to the noise matrix [39, 40]. The IBU mitigation scheme starts with the noise matrix Λ and

̸

̸

## Algorithm1 Algorithm for multi-qubit Bayesian measurement error mitigation

Input: Measurement outcomes for each shot, convergence tolerance ϵ

Output: Mitigated Multi-qubit state populations

1: Compute initial noisy multi-qubit state populations R k = 0, with k = 1 , ..., M ≤ 2 N q

2: Reduce subspace to mitigate to M states

3: while ∆( R ( n +1) , R ( n ) ) ≥ ϵ do

4: for all i = , with j i, j = 1 , ..., M do

5: Promote R i and R j to variables ρ i and ρ j

6: Perform Bayesian cycle to obtain Pr( ρ , ρ i j )

7: Update R i and R j

8: end for

9: Update M if some R i becomes zero

10: Compute ∆( R ( n +1) , R ( n ) )

11: end while

the noisy empirical distribution ρ noisy , and by applying Bayes' rule to an initial guess ρ 0 j it calculates the mitigated probability vectors ρ n +1 j via iteration, i.e.,

$$\rho _ { j } ^ { n + 1 } = \sum \rho _ { i } ^ { \text{noisy} } \frac { \Lambda _ { i j } \rho _ { j } ^ { n } } { \sum _ { m } \Lambda _ { i m } \rho _ { m } ^ { n } }. \quad \quad ( 9 ) \quad \quad,$$

However, while the IBU method uses Bayes' rule and starts with a prior guess, it does not result in a proper posterior distribution in the Bayesian sense. For this reason, it has been more appropriately referred to as iterative expectation maximization unfolding, and it is worth mentioning that more recently the IBU update rule has been vectorized to enable fast parallel computation on GPUs [36].

## III. BAYESIAN MITIGATION WITH BINARY OUTCOMES

In this section, we apply our Bayesian measurement error mitigation strategy to binary multiqubit outcomes obtained by executing quantum circuits on the superconducting qubits quantum computer VTTQ20 (Leena, QPU from IQM). The protocol begins with single-qubit measurement tomography to determine the uncorrelated error rates. While the service provider typically provides these error rates, we prefer extracting them explicitly in our work. We characterize the detector by performing two simple experiments: qubit reset followed by measurement and qubit flip to the excited state followed by measurement. Using N shots = 10 5 shots for each experiment, we obtain a two-by-two noise matrix Λ ( q ) for each qubit in the QPU. With these, we have access to the multiqubit noise matrix and we are now able to apply the Bayesian measurement error mitigation algorithm to the outcome of any quantum circuit.

To showcase the effectiveness of the measurement error mitigation scheme, we select a benchmark algorithm

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-33f32d78-9bdd-4103-a06d-903deb052c44.jpg)

Figure 3. (a): Quantum circuit for the preparation of a random multiqubit string: single qubit X gates are applied on random qubits. (b): Comparison between the mitigated (solid orange) and unmitigated (dashed gray) success probability of measuring the prepared bitstring as a function of the string length. (c): Mitigated (solid orange) success probability of measuring a prepared 16 qubits bitstring as a function of the number of shots. The results shown in (b) and (c) are averaged over 20 realizations of random bitstrings.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-c9bc4a0e-bfc5-4f55-a245-69863b60c09f.jpg)

designed to minimize errors from faulty gates and decoherence. Specifically, we use a multiqubit string preparation algorithm involving a single layer of single-qubit Xgates. As illustrated in Fig. 3(a), for a chosen bit string (e.g., '10..101'), we apply X-gates to specific qubits to flip their state from '0' to '1'. In these circuits, noisy readout is the primary source of errors, with potential imperfections in the rotation gate calibration accounted

Figure 4. (a): Success probability of measuring a single qubit as a function of N bin , representing the number of partitions used to construct the detector response function defined on the Q-axis. (b): Success probability of measuring a 20-qubit string as a function of N bin . (c): Comparison of mitigated success probabilities for measuring a prepared bitstring as a function of string length, for different number of partitions of the Q-axis established when reconstructing the detector response functions. The results shown in (a), (b), and (c) are obtained with 10 3 shots and averaged over 50 realizations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-e4f31cc6-42e1-4cdf-aaf3-5715bdb0dfae.jpg)

for within the noise matrix.

Thus, we prepare 20 random bit strings of different length, ranging from N q = 1 to N q = 19. After executing the quantum circuits and measuring the results, we obtain the unmitigated initial estimates. We then apply the Bayesian mitigation procedure outlined in the previous section, as described by Eqs. (4)-(7).

In Fig. 3(b), we display the success probability averaged over the 20 realizations as a function of the multiqubit string length, comparing it with the unmitigated probability. By comparing the mitigated results with the original unmitigated estimates, we can quantitatively assess the effectiveness of our error mitigation approach across different qubit lengths. We observe that the unmitigated success probability decreases rapidly as the string length increases, highlighting the challenge posed by measurement errors in longer sequences. On the other hand, applying the Bayesian measurement error mitigation results in significantly higher success probabilities across all lengths. Notably, for the longest string considered (19 qubits), the success probability improves from approximately 26% without mitigation to 92% after mitigation with just 10 3 shots.

length of N q = 16 and plot the mitigated and unmitigated success probabilities as a function of the number of shots. We observe that, while the unmitigated success probability shows only slight improvement with an increasing number of shots, the mitigated probability experiences a substantial enhancement. For instance, with 10 4 shots, the mitigated success probability reaches approximately 93%. This demonstrates the significant effectiveness of the Bayesian error mitigation technique in improving the accuracy of quantum measurements, even as the number of shots increases.

In Fig. 3 (c), we focus on a fixed multiqubit string

Furthermore, in the context of this string preparation algorithm, the unmitigated result, on average, tends to plateau near the product of the individual readout fidelities, which for the 16 qubits used in this experiment corresponds to approximately ∼ 1 / 3. In contrast, the Bayesian measurement error mitigation algorithm continuously refines with each shot, offering further advantages in its application. Although the rate of improvement diminishes over time, this behaviour is expected as Bayesian inference benefits from a larger number of measurements, allowing for a more precise estimation of the multiqubit state populations.

## IV. BAYESIAN MITIGATION WITH ANALOG DATA

In this section, we demonstrate the enhancements to the Bayesian measurement error mitigation scheme that can be achieved by leveraging readout analog data from the detector. Here, to validate these improvements, we run the algorithms on the superconducting qubits quantum computer 'Nazca' - IBM127.

The initial step involves performing measurement tomography to reconstruct the detector response functions. To do this, similar to measuring error rates, we execute two algorithms: one where we measure the qubit after reset, and one where we measure the qubit after applying an X-gate flip. In a superconducting qubit device, the typical measurement outcome is illustrated in Fig. 2 (top panel), where each measurement is depicted as a dot in the I-Q plane. Based on prior calibration, the standard method for assigning a measurement outcome as either '0' or '1' involves defining some threshold values or curves in this two-dimensional space. This assignment method is clearly not perfect. Instead of producing two distinct and sharp signals corresponding to the qubit's two possible states, the detector's response functions often exhibit significant overlap. This overlap is evident in the projected distributions shown in Fig. 2 (bottom panel), where the signals for the '0' and '1' states are not completely separated, causing then the assignment errors.

When dealing with readout analog data, the reconstruction of these histograms is analogous to the process of reconstructing the noise matrix for binary outcomes. Formally, the histogram shown in Fig. 2 can be described by a function such as

$$\Lambda _ { j } ( Q ) = \sum _ { i = 1 } ^ { N _ { \text{bin} } } \lambda _ { j i } \mathbb { I } [ Q _ { i }, Q _ { i + 1 } ] ( Q ), \quad \ \ ( 1 0 ) \, \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where Λ ( j Q ) represents the probability of the state being j when the detector measures the physical value Q , and I [ Q , Q i i +1 ]( Q ) is the indicator function defined as I [ Q , Q i i +1 ]( Q ) = 1 if Q ∈ [ Q , Q i i +1 ] and 0 otherwise. λ ji is then the probability of the state being j when the measured Q is in the i-th interval [ Q , Q i i +1 ]. Here, N bin denotes the number of intervals used to discretize the Q-axis of the IQ-plane and serves as a parameter representing the discretization applied in modelling the detector response functions. While an analytical expression for these distributions would be ideal, and although a bimodal Gaussian can often adequately capture the behaviour of the readout functions, it is typically more straightforward and more efficient to use directly the histograms. This approach is somewhat similar to defining Q-dependent error rates along the Q-axis partition.

Therefore, once we gather data from N shots = 10 5 shots, we establish a number of partitions, denoted as N bin , and reconstruct response functions for each individual qubit. These response functions are then utilized in the multiqubit conditional probability of Eq. (8) to be used in the the update cycle of the Bayesian measurement error mitigation algorithm. In this framework, when we consider N bin = 2, this approach is equivalent to working with binary outcomes.

In Fig. 4, we investigate the influence of varying the number of partitions, N bin , on Bayesian measurement error mitigation when using analog data. Therefore, we prepare 50 random bit strings across various qubit lengths and record the measurement outcomes in analog form. After obtaining the unmitigated initial estimates, we apply the Bayesian measurement error mitigation procedure, leveraging the previously characterized Q-dependent detector response functions. By examining how different values of N bin impact the mitigation process, we aim to optimize the partitioning strategy for improved readout error mitigation. This analysis allows us to understand the trade-offs between the granularity in the analog data and the effectiveness of the error mitigation, providing insights into how to best utilize the analog information obtained in quantum measurement. In Fig. 4(a), we present the mitigated single-qubit success probability averaged over the 50 realisation as a function of N bin . We observe a noticeable trend where the success probability, which in the single qubit case is akin to the readout fidelity, improves as N bin increases. This trend becomes more pronounced and stable when analyzing longer bit strings. In Fig. 4(b), we display the mitigated success probability averaged over the 50 realisations for a 20-qubit string. Also in this case, we observe that the mitigated success probability increases as the response functions are discretized with a greater number of partitions N bin . However, these improvements eventually plateau beyond a certain threshold value of N bin .

In Fig. 4(c), we display the averaged mitigated success probabilities as a function of the string length for different values of N bin . We see that increasing the granularity of the detector response function constantly result in better mitigated values for all string lengths. For the case of N q = 20, the success probability significantly improves from approximately ∼ 89% to about ∼ 94% when the discretization of the response functions is increased from N bin = 2 to N bin = 10. Both cases represent a substantial improvement over the unmitigated result, which even though is not shown, it is around ∼ 25%. These findings underscore how leveraging the additional information encoded in the analog data can enhance the effectiveness of the Bayesian measurement error mitigation algorithm with minimal to none additional computational cost. Therefore, by optimizing the partitioning strategy for the detector analog data, we can achieve higher accuracy in the quantum measurements, thereby improving the overall performance of quantum algorithms on NISQ devices.

Figure 5. Comparison between success probabilities of measuring the prepared bit-strings as a function of the string length mitigated according to different readout error mitigation strategies. In solid blue and solid green we report the success probability mitigated trough Iterative Bayesian Unfolding and noise Matrix Invervesion Method respectively. In dashed orange and dashed yellow, we display the success probabilities mitigated trough our implementation of Bayesian measurement error mitigation using two partitions of the Qaxis for the response functions, N bin = 2 and N bin = 20 respectively. All the results displayed are obtained by averaging over 50 realisations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-74dc8717-e956-4526-a3e8-bc0f94bd7ba5.jpg)

## V. COMPARISON WITH OTHER MITIGATION TECHNIQUES

The results presented in the previous section provide compelling evidence of the efficacy of our implementation of Bayesian measurement error mitigation in reducing readout errors in multiqubit algorithms. However, it is relevant to compare the performances of our Bayesian approach with other established measurement error mitigation techniques. For benchmark purposes, we selected two techniques. The first benchmark measurement error mitigation scheme is the Iterative Bayesian Unfolding method (IBU) introduced at the end of Sec. II, which shares a probabilistic formulation linking correct qubit populations with the measured outcomes. However, IBU employs iterative expectation maximization rather than a formal prior-posterior update iteration in the Bayesian sense. The second measurement error mitigation technique is based on the Noise Matrix Inversion Method (MIM). This method applies the inverse of the noise matrix to the measured outcomes and, if the resulting mitigated probability vector is non-physical, minimize the squared distance squared | Λ -1 ρ noise - | ρ 2 , enforcing the constraints of positivity and normalization.

In Fig. 5, we compare the success probabilities of the previously prepared random bitstrings mitigated using IBU, MIM, and our Bayesian measurement error mitigation method as a function of the multiqubit string length, averaged over 50 realizations. All methods effectively mitigate measurement errors, demonstrating their viability. For our Bayesian method we employ two extreme values for the discretization parameter N bin of the detector response function.

When applying our Bayesian mitigation algorithm to binary outcomes, i.e fixing N bin = 2, we achieve results very similar to IBU, with marginal improvements observed for the longer strings considered in our study. This similarity is expected because both methods start from similar probabilistic foundations but employ different iterative approaches to reach the optimization/maximization. On the other hand, by increasing the discretization parameter of the response function, our Bayesian method outperforms IBU. This enhancement suggests that our approach effectively leverages additional information encoded in the readout signal, which contributes to improved mitigated outcomes. Comparing our Bayesian method with MIM reveals a peculiar trend. Initially, MIM consistently outperforms our Bayesian method when applied to binary outcomes. However, with N bin = 20, we observe a transition across different qubit numbers: MIM excels for shorter strings, whereas our method begins to surpass MIM for longer strings. Specifically, this transition appears for qubit numbers greater than 5 in our experimental results.

These comparisons highlight the strengths of each method under different conditions, showcasing how our Bayesian measurement error mitigation approach applied to the detector analog data can adapt and potentially outperform more established traditional methods like IBU and MIM when mitigating the outcome of large multiqubit experiments.

## VI. SCALING ANALYSIS AND OTHER CONSIDERATIONS

̸

To conclude, we estimate the computational complexity for the Bayesian multiqubit measurement error mitigation algorithm 1. A single Bayesian iteration, as described in Eq. (2), scales linearly with the number of shots N shots . This process must be repeated for each pair of states in the truncated subspace. If initially we have M strings with counts such that ρ noisy i = 0, we have M M ( -1) 2 populations pairs to consider and update. Thus, if all pairs are updated N r times, the overall computational complexity to post-process the results of a multiqubit experiment scales as ∝ N r N shots M M ( -1) 2 . This indicates a linear relationship with the number of shots and iterations and a quadratic relationship with the number of states after truncation based on the initial noisy counts. Therefore, the primary computational cost for large systems arises from the number of initial noisy populations, which is intricately linked to the number of shots and the extent to which the quantum state is spread over the multiqubit Hilbert space. Hence, the algorithm is more efficient when the outcome of the quantum algorithm is concentrated on a limited number of populations.

To determine the necessary number of updates for the multiqubit state populations, Fig. 6 (a) shows the total variation distance between the probability vectors be-

Figure 6. (a): Total variation distance between the mitigated probabilities after each Bayesian iteration as a function of the iteration number for different string lengths. (b): Time required by our Bayesian mitigation scheme to mitigate noisy outcomes as a function of the number of qubits measured. (c): Time required by our Bayesian mitigation scheme to mitigate the outcomes of the string preparation experiment run on 20 qubits as a function of the number of shots. The results are averaged over 50 realizations and correspond to the algorithm using N bin = 20 for the discretization of the detector response function.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00869v1-figure-61099115-0935-499a-8738-d57aeb08f1cd.jpg)

fore and after each Bayesian update from Eq. (7) for each population pair. The results for qubit counts of N q = 5 10 and 20 indicate that convergence is achieved , after a few iterations as the figure of merit reaches a plateau. Consequently, we set a threshold to terminate the Bayesian measurement error mitigation algorithm when the total variation distance between probability vectors before and after an iteration falls below ϵ = 10 -3 . To ensure the algorithm exits the iteration loop, we imposed a maximum of 20 iterations, which was never reached in our simulations.

## VII. CONCLUSIONS

With this exit condition, we recorded the time required to run the Bayesian mitigation algorithm. In Fig. 6 (b) and (c), we display the average time to mitigate the outcome of the random string state preparation algorithm as a function of the number of qubits and the number of shots for a circuit with 20 qubits, respectively. The timing includes all the steps necessary to perform the mitigation algorithm as detailed in 1. Using a conventional quad-core laptop, the post-processing of the experiment data was within 1 second in the worst case considered here, i.e. N = 20 and N q shots = 2000.

In this work, we have introduced and conducted a comprehensive analysis of an efficient measurement error mitigation scheme that significantly enhances the accuracy of the readout step in multiqubit experiments on quantum computers. The algorithm combines Bayesian inference with a heuristic binary optimization cycle to construct a probability distribution for the multiqubit state populations, based on the known noise matrix affecting the measurement and state assignment. This approach effectively incorporates noise and imperfections in the post-processing of measurement outcomes. Crucially, our Bayesian framework is designed to avoid assigning unphysical values, ensuring that the mitigated outcomes remain consistent with the experimental data at each step of the mitigation algorithm. We tested the scalability and performance of the algorithm on two real superconducting qubit quantum computers, demonstrating a substantial improvement in the precision of multiqubit bitstrings with minimal computational overhead. Additionally, we have shown that it is more advantageous to apply readout mitigation techniques at the level of the detector's analog data. For superconducting QPUs, this involves working with data in the IQ-plane rather than aggregated binary counts. We compared our method against the state-

of-the-art Iterative Bayesian Unfolding (IBU) mitigation technique and found that our approach can outperform IBU, especially when incorporating analog readout data into the mitigation scheme. Although we focused on superconducting qubit devices, our mitigation approach is versatile and could be applied to other quantum computing platforms, either directly applying the mitigation procedure on the aggregated bitstring counts or using any detector analog data if available. This flexibility underscores the broad applicability and potential impact of our Bayesian measurement error mitigation framework in advancing the field of quantum computing.

- [1] E. Hyypp¨ a, S. Kundu, C. F. Chan, A. Gunyh´, J. Hotari, o D. Janzso, K. Juliusson, O. Kiuru, J. Kotilahti, A. Landra, W. Liu, F. Marxer, A. M¨ akinen, J.-L. Orgiazzi, M. Palma, M. Savytskyi, F. Tosto, J. Tuorila, V. Vadimov, T. Li, C. Ockeloen-Korppi, J. Heinsoo, K. Y. Tan, J. Hassel, and M. M¨tt¨nen, Nature Communications o o 13 , 6895 (2022).
- [2] F. Marxer, A. Veps¨ al¨inen, a S. W. Jolin, J. Tuorila, A. Landra, C. Ockeloen-Korppi, W. Liu, O. Ahonen, A. Auer, L. Belzane, V. Bergholm, C. F. Chan, K. W. Chan, T. Hiltunen, J. Hotari, E. Hyypp¨ a, J. Ikonen, D. Janzso, M. Koistinen, J. Kotilahti, T. Li, J. Luus, M. Papic, M. Partanen, J. R¨bin¨, J. Rosti, M. Savyta a skyi, M. Sepp¨ al¨, a V. Sevriuk, E. Takala, B. Tarasinski, M. J. Thapa, F. Tosto, N. Vorobeva, L. Yu, K. Y. Tan, J. Hassel, M. M¨tt¨nen, and J. Heinsoo, PRX Quantum o o 4 , 010314 (2023).
- [3] U. R´ eglade, A. Bocquet, R. Gautier, J. Cohen, A. Marquet, E. Albertinale, N. Pankratova, M. Hall´n, e F. Rautschke, L.-A. Sellem, et al. , Nature , 1 (2024).
- [4] P. S. Mundada, A. Barbosa, S. Maity, Y. Wang, T. Merkh, T. Stace, F. Nielson, A. R. Carvalho, M. Hush, M. J. Biercuk, and Y. Baum, Phys. Rev. Appl. 20 , 024034 (2023).
- [5] L. Skoric, D. E. Browne, K. M. Barnes, N. I. Gillespie, and E. T. Campbell, Nature Communications 14 , 7040 (2023).
- [6] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O'Brien, Rev. Mod. Phys. 95 , 045005 (2023).
- [7] S. N. Filippov, S. Maniscalco, and G. Garc´ ıa-P´rez, arXiv e preprint arXiv:2403.13542 (2024).
- [8] Y. Li and S. C. Benjamin, Phys. Rev. X 7 , 021050 (2017).
- [9] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119 , 180509 (2017).
- [10] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, et al. , Nature 618 , 500 (2023).
- [11] E. Van Den Berg, Z. K. Minev, A. Kandala, and K. Temme, Nature Physics 19 , 1116 (2023).
- [12] D. Bluvstein, S. J. Evered, A. A. Geim, S. H. Li, H. Zhou, T. Manovitz, S. Ebadi, M. Cain, M. Kalinowski, D. Hangleiter, et al. , Nature 626 , 58 (2024).
- [13] A. Feito, J. S. Lundeen, H. Coldenstrodt-Ronge, J. Eisert, M. B. Plenio, and I. A. Walmsley, New Journal of Physics 11 , 093038 (2009).

## ACKNOWLEDGMENTS

We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team. This work has been partly funded by the Business Finland Co-Innovation Project 40561/31/2020 Quantum Technologies Industrial (QuTI).

- [14] J. S. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. L. Pregnell, C. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, Nature Physics 5 , 27 (2009).
- [15] F. B. Maciejewski, Z. Zimbor´ as, and M. Oszmaniec, Quantum 4 , 257 (2020).
- [16] M. Cattaneo, M. A. C. Rossi, K. Korhonen, E.-M. Borrelli, G. Garc´ ıa-P´rez, e Z. Zimbor´s, and D. Cavalcanti, a Phys. Rev. Res. 5 , 033154 (2023).
- [17] S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, Phys. Rev. A 103 , 042605 (2021).
- [18] Y. Chen, M. Farahzad, S. Yoo, and T.-C. Wei, Phys. Rev. A 100 , 052315 (2019).
- [19] R. Navarathna, T. Jones, T. Moghaddam, A. Kulikov, R. Beriwal, M. Jerger, P. Pakkiam, and A. Fedorov, Applied Physics Letters 119 , 114003 (2021).
- [20] B. Lienhard, A. Veps¨ al¨inen, a L. C. Govia, C. R. Hoffer, J. Y. Qiu, D. Rist` e, M. Ware, D. Kim, R. Winik, A. Melville, B. Niedzielski, J. Yoder, G. J. Ribeill, T. A. Ohki, H. K. Krovi, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Phys. Rev. Appl. 17 , 014024 (2022).
- [21] U. Azad and H. Zhang, 2022 IEEE/ACM 7th Symposium on Edge Computing (SEC) , 362 (2022).
- [22] H. G. Ahmad, R. Schiattarella, P. Mastrovito, A. Chiatto, A. Levochkina, M. Esposito, D. Montemurro, G. P. Pepe, A. Bruno, F. Tafuri, et al. , Advanced Quantum Technologies , 2300400 (2024).
- [23] F. Cosco and N. Lo Gullo, Phys. Rev. A 108 , L060402 (2023).
- [24] H. T. Dinani, D. W. Berry, R. Gonzalez, J. R. Maze, and C. Bonato, Phys. Rev. B 99 , 125413 (2019).
- [25] R. Puebla, Y. Ban, J. Haase, M. Plenio, M. Paternostro, and J. Casanova, Phys. Rev. Appl. 16 , 024044 (2021).
- [26] D. Leibfried, E. Knill, S. Seidelin, J. Britton, R. B. Blakestad, J. Chiaverini, D. B. Hume, W. M. Itano, J. D. Jost, C. Langer, et al. , Nature 438 , 639 (2005).
- [27] B. Teklu, S. Olivares, and M. G. A. Paris, Journal of Physics B: Atomic, Molecular and Optical Physics 42 , 035502 (2009).
- [28] K. N. Cassemiro, K. Laiho, and C. Silberhorn, New Journal of Physics 12 , 113052 (2010).
- [29] S. Paesani, A. A. Gentile, R. Santagati, J. Wang, N. Wiebe, D. P. Tew, J. L. O'Brien, and M. G. Thompson, Phys. Rev. Lett. 118 , 100503 (2017).
- [30] K. T. Laverick, I. Guevara, and H. M. Wiseman, Phys. Rev. A 104 , 032213 (2021).

- [31] V. Gebhart, A. Smerzi, and L. Pezz`, Phys. Rev. Appl. e 16 , 014035 (2021).
- [32] S. Duffield, M. Benedetti, and M. Rosenkranz, Machine Learning: Science and Technology 4 , 025007 (2023).
- [33] C. A. Pattison, M. E. Beverland, M. P. da Silva, and N. Delfosse, arXiv preprint arXiv:2107.13589 (2021).
- [34] H. Ali, J. Marques, O. Crawford, J. Majaniemi, M. Serra-Peralta, D. Byfield, B. Varbanov, B. M. Terhal, L. DiCarlo, and E. T. Campbell, arXiv preprint arXiv:2403.00706 (2024).
- [35] P. D. Nation, H. Kang, N. Sundaresan, and J. M. Gambetta, PRX Quantum 2 , 040326 (2021).
- [36] B. Pokharel, S. Srinivasan, G. Quiroz, and B. Boots, Phys. Rev. Res. 6 , 013187 (2024).
- [37] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 69 , 062320 (2004).
- [38] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93 , 025005 (2021).
- [39] B. Nachman, M. Urbanek, W. A. de Jong, and C. W. Bauer, npj Quantum Information 6 , 84 (2020).
- [40] R. Hicks, C. W. Bauer, and B. Nachman, Phys. Rev. A 103 , 022407 (2021).
F. Cosco F. Plastina N. Lo Gullo
Bayesian mitigation of measurement errors in multi-qubit experiments
We introduce an implementation of Bayesian measurement error mitigation
tailored for multiqubit experiments on near-term quantum devices. Our approach
leverages complete information from the readout signal, which is available
before any binary state assignment of the qubits. We provide a detailed
algorithm workflow, from the calibration of detector response functions to the
post-processing of measurement outcomes, offering a computationally efficient
solution suitable for the output size typical of current quantum computing
devices and quantum algorithms. We benchmark our protocol on actual quantum
computers with superconducting qubits, where the readout signal encodes the
measurement information in the IQ clouds before qubit state assignment.
Finally, we compare the performance of our algorithm against other measurement
error mitigation methods, such as iterative Bayesian unfolding and noise matrix
inversion.