In <ref>, we plot the time in seconds to construct the generator \(\G\) for a given assessment, using our three different methods.
We consider the average over 7 different assessments.
We see that the first method (x) is always the fastest.
This is because it corresponds to simple subtractions of vectors without additional checks.
At the same time the second method seems to be very close to the first method.

For the third method, it seems to take much longer to construct the generator than for the other two methods, for example for \(C_{\text{imp}}\) for an assessment containing \(15\) pairs of option sets, each containing between two and eight options per pair, it took on average 31 minutes and 36 seconds.
This is because, compared to the second method, for the third method we have to construct \(\G_{\square}\) from \(\D_{\triangle}\) (in memory) while also doing the simplifications that make the difference between \(\G_{\triangle}\) and \(\G_{\square}\).
For the third method, we also plotted error bars for the minimum and maximum time, but not for the others because they are already very close to each other and the error bars were also very small.
From the error bars it can be seen that there is also a large spread in how long it takes to construct the generator \(\G_{\square}\), except for \(C_{\text{lin}}\).
For \(C_{\text{imp}}\) a large jump can be seen in the number of option sets in the average of \(\G_{\square}\), but this is due to a large jump in the maximum value, as the minimum value does not jump up until later.

In <ref>, we plot the time to choose from a single option set; that is, to evaluate the natural extension \(C_{\Ass_{1:\ell}}\) for a single option set with <ref>, using either \(\G_{\times}\), \(\G_{\triangle}\) or \(\G_{\square}\).
The points in the figure are an average over the same 7 assessments \(\Ass_{1:L}\) that we have used in <ref>.

Additionally, for each assessment \(\Ass_{1:L}\), we also generated 7 random options sets to choose from, each containing between 2 and 8 options.
So if the assessments are \(\Ass_{1:L}^1,...,\Ass_{1:L}^7\), we have option sets \(B_1^1,...,B_7^1\) corresponding to \(\Ass_{1:L}^1\) up to \(B_1^7,...,B_7^7\) corresponding to \(\Ass_{1:L}^7\).
Every point on <ref> is the average of the time to calculate \(C_{\Ass_{1:\ell}^r}(B_k^r)\) over \(k\in\{1,...,7\}\) and \(r\in \{1,...,7\}\), for every \(\ell\), using one of our three methods.
The number of options in every option set was chosen uniformly random from \(\{2,...,8\}\), and the options themselves are again generated uniformly random from \([0,1]^4\).
We see in <ref> that our simplifications indeed work, in the sense that more simplifications result in smaller choosing times.
This is what we expected, as a smaller generator will lead to a faster loop over the generator in <ref>.

Next, in <ref>, we plot the sum of the time it takes to construct the generator and the time it takes to use this generator to evaluate the natural extension using <ref>.
For this plot the data points stop when the calculations take more than 50 minutes on average.
This in fact already occurred in <ref>, but is not visible there because we only plot up to 100 seconds there.
For all methods, we see that it is initially faster to use the second method but that the third method eventually becomes faster for larger assessments.
How fast this happens depends on the type of assessment, and we see that for \(C_{\text{max}}\) and \(C_{\text{adm}}\) it happens faster than for \(C_{\text{lin}}\) and \(C_{\text{imp}}\).
For \(C_{\text{imp}}\), it may seem as if the third method is never better, but this is not the case because from \(\ell=13\) onwards, the results of the second method are not displayed because the calculations took too long; so for these higher values of \(\ell\), the third method is again better.

<ref> only considers the time for making one choice though.
If we evaluate the natural extension for several option sets, for the same assessment, then there will come a point where the fact that the time investment \(t_{\text{pre},\square}\) for the preprocessing of the third method is higher than to the time \(t_{\text{pre},\triangle}\) for the preprocessing of the second method is compensated by the fact that the time \(t_{\text{choose},\square}\) to evaluate the natural extension for the third method is lower than the time \(t_{\text{choose},\triangle}\) to evaluate the natural extension for the second method.
The number of evaluations \(n\) at which the total times for both methods break even is the smallest \(n\) such that \(t_{\text{pre},\square}+n t_{\text{choose},\square} \leq t_{\text{pre},\triangle}+n t_{\text{choose},\triangle}\).
If \(t_{\text{choose},\triangle}>t_{\text{choose},\square}\)—which is true in all our experiments—solving for \(n\) gives

\[
n \geq \frac{t_{\text{pre},\square}- t_{\text{pre},\triangle}}{t_{\text{choose},\triangle}- t_{\text{choose},\square}}.
\]

We plot the right-hand side of this inequality in <ref> for the two most interesting cases: \(C_{\text{lin}}\) and \(C_{\text{imp}}\).

The main conclusion is that even for small assessments, the third method is still faster  than the second method, provided that the number \(n\) of option sets that we need to evaluate is high enough.
For \(C_{\text{lin}}\) we found the highest break-even point for \(\ell=9\), for which the third method is faster if we need to evaluate the natural extension for at least 21 option sets.
For \(C_{\text{imp}}\), we found that the highest value was for \(\ell=7\), for which we would need to evaluate at least 70 option sets before the third method becomes faster.

[ ymin=5*10^-3,ymax=10^2,
xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] kiestijden/_1_-stupid.csv;
[mark=triangle,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] kiestijden/_1_-simplify.csv;

[mark=square,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] kiestijden/_1_-xsimplify.csv;

\(C_{\text{lin}}\)

[ ymin=5*10^-3,ymax=10^2,
xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=magenta,solid,magenta,domain=1:4,select coords between index=03] table[x=l,y=d1,col sep=comma] kiestijden/_4_-stupid.csv;
[mark=triangle,mark options=magenta,solid,magenta] table[x=l,y=d1,col sep=comma] kiestijden/_4_-simplify.csv;

[mark=square,mark options=magenta,solid,magenta] table[x=l,y=d1,col sep=comma] kiestijden/_4_-xsimplify.csv;

\(C_{\text{max}}\)

xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] kiestijden/_1,_1,_1_-stupid.csv;
[mark=triangle,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] kiestijden/_1,_1,_1_-simplify.csv;

[mark=square,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] kiestijden/_1,_1,_1_-xsimplify.csv;

\(C_{\text{adm}}\)

xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] kiestijden/_4,_4,_4_-stupid.csv;
[mark=triangle,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] kiestijden/_4,_4,_4_-simplify.csv;

[mark=square,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] kiestijden/_4,_4,_4_-xsimplify.csv;

\(C_{\text{imp}}\)

Time in seconds to evaluate the natural extension for a single option set using \(\G_{\times}\), \(\G_{\triangle}\) or \(\G_{\square}\), as a function of the number \(\ell\) of pairs \((V,W)\) in the assessment. The vertical axis is logarithmic, and the experiments are averaged over 7 separate experiments and 7 option sets to choose from per experiment.

xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1_-stupidcombo.csv;
[mark=triangle,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1_-simplifycombo.csv;

[mark=square,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1_-xsimplifycombo.csv;

\(C_{\text{lin}}\)

[ ymin=5*10^-3,ymax=5*10^3,
xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=magenta,solid,magenta,domain=1:4,select coords between index=03] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4_-stupidcombo.csv;
[mark=triangle,mark options=magenta,solid,magenta] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4_-simplifycombo.csv;

[mark=square,mark options=magenta,solid,magenta] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4_-xsimplifycombo.csv;

\(C_{\text{max}}\)

xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1,_1,_1_-stupidcombo.csv;
[mark=triangle,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1,_1,_1_-simplifycombo.csv;

[mark=square,mark options=blue,solid,blue] table[x=l,y=d1,col sep=comma] koos+kiestijden/_1,_1,_1_-xsimplifycombo.csv;

\(C_{\text{adm}}\)

xmin=1, xmax=15,
xlabel=$\ell$,ylabel=time in seconds,

[mark=x,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4,_4,_4_-stupidcombo.csv;
[mark=triangle,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4,_4,_4_-simplifycombo.csv;

[mark=square,mark options=indigo,solid,indigo] table[x=l,y=d1,col sep=comma] koos+kiestijden/_4,_4,_4_-xsimplifycombo.csv;

\(C_{\text{imp}}\)

Total time in seconds to `construct' a \(\G_{\times}\), \(\G_{\triangle}\) or \(\G_{\square}\) and subsequently use it to evaluate the natural extension for a single option set, as a function of the number \(\ell\) of pairs \((V,W)\) in the assessment. The vertical axis is logarithmic, the experiments are averaged over 7 separate experiments (with for each experiment 7 option sets to choose from) and the lines for each method stop when the average was over 50 minutes.

[xlabel=$\ell$,ylabel=$n$, ]
[mark=pentagon,mark options=LimeGreen,solid,LimeGreen] table[x=l,y=fraction_result,col sep=comma,only marks] ratios/ratios-_1_-simplify.csv;

\(C_{\text{lin}}\)

[xlabel=$\ell$,ylabel=$n$, ]
[mark=pentagon,mark options=indigo,solid,indigo] table[x=l,y=fraction_result,col sep=comma,only marks] ratios/ratios-_4,_4,_4_-simplify.csv;

\(C_{\text{imp}}\)

The vertical axis depicts for how many option sets with (between 2 and 8 options each) you would have to evaluate the natural extension before it is better to use the third method rather than the second method for \(C_{\text{lin}}\) and \(C_{\text{imp}}\).
For \(C_{\text{imp}}\) the plot only goes up to \(\ell=12\) because for higher \(\ell\) it took too long to evaluate the natural extension using the second method.

§ CONCLUSION AND FUTURE WORK

The main conclusion of this work is that choice functions provide a principled viable framework for inferring new decisions from previous ones.
The two key concepts that we introduced to achieve this were consistency and natural extension.
The former allows one to check if there is at least one preference order that is compatible with an assessment, while the latter allows one to infer new choices based on the assessment by considering all compatible orders.
From a practical point of view, our main contributions are algorithms that are able to execute these tasks.
We have also developed effective methods to simplify the information in the assessment, while still retaining all decisive power.
To demonstrate the performance when including these simplifications, we carried out experiments that measure the space and time efficiency under various scenarios.
Our results show that our simplifications allow for larger assessments to be dealt with, and this for various levels of imprecision.
Related to this, we saw that the third method—which implements all the simplifications—outperforms the other two methods in terms of efficiency for larger assessments.
For smaller assessments, the second method is the most efficient, unless one has to evaluate the natural extension for multiple option sets, in which case the third method eventually becomes faster.
A particularly intriguing observation was that for assessments that are generated with maximality or E-admissibility, the generator of the third method even seemed to stabilise, as seen in <ref>.
This is intriguing because it could mean that it is converging towards the `real' model—the model that we used to make the assessment—through simple approximate models.
It would be interesting to investigate further if this trend continues for even larger assessments, and explain why it happens.

Future work could also add onto <ref> by trying to obtain a `simplest' generator for any given assessment, or showing that our methods already achieves this.
Also, how often and in what order one has to apply the simplifications can still be optimised.
Approaches in between method 2 and 3 could for example be considered.
Another direction could be to analyse how efficiency scales with other parameters, both theoretically and experimentally.
This includes the size of the option sets in the assessment, the dimension of the vector space $\V$ and the size of the option set $A$ for which we want to evaluate the natural extension.
One could also consider  alternative forms of assessments, such as bounds on probabilities, bounds on expectations and preference statements, and show how they can be made to fit in our choice functions framework.
Finally, it would be nice to evaluate how our methods perform when applied to real-life decision problems.

§.§.§ Acknowledgements

This work was partially supported by Ghent University, through Jasper De Bock's BOF Starting Grant “Rational decision-making under uncertainty: a new paradigm based on choice functions”, number 01N04819.

[Troffaes, 2007]
authorM. C. M. Troffaes,
titleDecision making under uncertainty using imprecise probabilities,
journalInternational Journal of Approximate Reasoning volume45 (year2007) pages17–29.

[De Bock, Jasper, 2020]
authorDe Bock, Jasper,
titleArchimedean choice functions: an axiomatic foundation for imprecise decision making,
in: booktitleInformation Processing and Management of Uncertainty in Knowledge-Based Systems, volume volume1238, publisherSpringer, year2020, pp. pages195–209.

[De Bock and De Cooman, 2019]
authorJ. De Bock, authorG. De Cooman,
titleInterpreting, axiomatising and representing coherent choice functions in terms of desirability,
in: booktitleProceedings of the eleventh International Symposium on Imprecise Probabilities: Theories and Applications, organizationPMLR, year2019, pp. pages125–134.

[Seidenfeld et al., 2010]
authorT. Seidenfeld, authorM. J. Schervish, authorJ. B. Kadane,
titleCoherent choice functions under uncertainty,
journalSynthese volume172 (year2010) pages157.

[Van Camp et al., 2018]
authorA. Van Camp, authorE. Miranda, authorG. De Cooman,
titleNatural extension of choice functions,
in: booktitleProceedings of the seventeenth International Conference on Information Processing and Management of Uncertainty in Knowledge-Based Systems, organizationSpringer, year2018, pp. pages201–213.

[Sen, 1971]
authorA. K. Sen,
titleChoice functions and revealed preference,
journalThe Review of Economic Studies volume38 (year1971) pages307–317.

[Sen, 1977]
authorA. Sen,
titleSocial choice theory: A re-examination,
journalEconometrica: journal of the Econometric Society  (year1977) pages53–89.

[Birkhoff, 1940]
authorG. Birkhoff, titleLattice theory, editionthird ed., publisherAmerican Mathematical Society, year1940.

[Augustin et al., 2014]
authorT. Augustin, authorF. P. Coolen, authorG. De Cooman, authorM. C. Troffaes, titleIntroduction to Imprecise Probabilities, publisherJohn Wiley & Sons, year2014.

[Van Camp et al., 2018]
authorA. Van Camp, authorG. De Cooman, authorE. Miranda, authorE. Quaeghebeur,
titleCoherent choice functions, desirability and indifference,
journalFuzzy sets and systems volume341 (year2018) pages1–36.

[Quaeghebeur, 2013]
authorE. Quaeghebeur,
titleThe CONEstrip algorithm,
in: booktitleSynergies of Soft Computing and Statistics for Intelligent Data Analysis, publisherSpringer, year2013, pp. pages45–54.

[Matoušek and Gärtner, 2006]
authorJ. Matoušek, authorB. Gärtner, titleUnderstanding and Using Linear Programming, Universitext, publisherSpringer, year2006.

[Vaidya, 1989]
authorP. M. Vaidya,
titleSpeeding-up linear programming using fast matrix multiplication,
in: booktitle30th annual symposium on foundations of computer science, organizationIEEE Computer Society, year1989, pp. pages332–337.

[Barber et al., 1996]
authorC. B. Barber, authorD. P. Dobkin, authorH. Huhdanpaa,
titleThe quickhull algorithm for convex hulls,
journalACM Trans. Math. Softw. volume22 (year1996) pages469–483.

[Fukuda and Prodon, 1996]
authorK. Fukuda, authorA. Prodon,
titleDouble description method revisited,
in: booktitleCombinatorics and Computer Science, publisherSpringer Berlin Heidelberg, addressBerlin, Heidelberg, year1996, pp. pages91–111.

[Fukuda, 2020]
authorK. Fukuda, titlePolyhedral computation, publisherDepartment of Mathematics, Institute of Theoretical Computer Science ETH Zurich, year2020.

[Clarkson, 1994]
authorK. Clarkson,
titleMore output-sensitive geometric algorithms,
in: booktitleProceedings of the 35th Annual Symposium on Foundations of Computer Science, year1994, pp. pages695–702.

[Ottmann et al., 1995]
authorT. Ottmann, authorS. Schuierer, authorS. Soundaralakshmi,
titleEnumerating extreme points in higher dimensions,
in: booktitleSTACS, volume volume95, year1995, pp. pages562–570.

[Legat, 2023]
authorB. Legat,
titlePolyhedral computation,
in: booktitleJuliaCon, year2023. <https://pretalx.com/juliacon2023/talk/JP3SPX/>.

[Legat, 2020]
authorB. Legat, titleSet programming: theory and computation, Ph.D. thesis, UCL-Université Catholique de Louvain, year2020.

[Nakharutai et al., 2018]
authorN. Nakharutai, authorM. C. M. Troffaes, authorC. C. Caiado,
titleImproved linear programming methods for checking avoiding sure loss,
journalInternational Journal of Approximate Reasoning volume101 (year2018) pages293–310.

§ EQUIVALENCE OF AXIOMS

In this appendix we prove the equivalence of our characterisation of coherence and the axiomatic characterisation of De Bock and De Cooman, 2019.
In particular, with \(\R^2_{>0}\coloneqq \{ (\lambda,\nu)\in \R^2 \colon (\lambda,\nu)>0\} \), we consider the following axioms for a rejection function \(R\colon \Q\to \Q\):
for all \(A,B\in \Q\), \(u\in A\) and \((\lambda_1,\lambda_2)\colon A\times B \to \R^2_{>0}\):

${\mathrm{R}}_0$. \(u \in {R}(A)\) if and only if \(0\in {R}(A-u)\);

${\mathrm{R}}_1$. \({R}(A)\neq A\);

${\mathrm{R}}_2$. if \(u>0\) then \(0\in {R}(\{0,u\})\);

${\mathrm{R}}_3$. if \(0\in {R}(A\cup\{0\})\) and  \(0\in {R}(B\cup\{0\})\) then
\[0\in {R}(\{\lambda_1(v,w) v +\lambda_2(v,w) w \colon v\in A,w\in B\}\cup \{0\});\]

${\mathrm{R}}_4$. if \(A\subseteq B\) then \({R}(A)\subseteq {R}(B)\).

Our main result is then the following.

A choice function \(C\colon \Q\to \Q\) is coherent if and only if its corresponding rejection function \(R_C\) satisfies ${\mathrm{R}}_0$–${\mathrm{R}}_4$.

A subtle point is that these are not exactly the axioms used by De Bock and De Cooman, 2019, as they also add \(\emptyset\) to the domain of their choice and rejection functions.
We will however see that this is a mere technicality with no substantial effect.

To be able to use the results of De Bock and De Cooman, 2019 in our proof for <ref>, we too will consider choice and rejection functions whose domain includes \(\emptyset\) and call it them extended choice and rejection functions: a map \(\hat{C}\colon \Qe\to\Qe\) is called an extended choice function if it satisfies \(\hat{C}(A)\subseteq A\) for all \(A\in \Qe\).
With any such extended choice function \(\hat{C}\) we also associate a corresponding extended rejection function \(\hat{R}_{\hat{C}}\colon \Qe\to \Qe\), defined by \(\hat{R}_{\hat{C}}(A)=A\setminus \hat{C}(A)\) for all \(A\in \Qe\).
Since we must have that \(\hat{C}(\emptyset)\subseteq \emptyset\), we always have that \(\hat{C}(\emptyset)=\emptyset\) and \(\hat{R}(\emptyset)=\emptyset\), which makes this part of the domain uninteresting.

With any choice function \(C\colon \Q\to\Q\) we associate the extended choice function

\begin{equation}\label{eq:RextDef}
\overline{C} \colon \Qe\to\Qe\colon A \mapsto \begin{cases} C(A) &\text{if }A\neq \emptyset,\\
\emptyset &\text{if }A=\emptyset.
\end{cases}
\end{equation}

This link between extended choice functions and choice functions is bijective because \(\hat{C}(\emptyset)=\emptyset\) for every extended choice function \(\hat{C}\).

For any extended rejection function \(\hat{R}\colon \Qe\to \Qe\), De Bock and De Cooman, 2019 consider the following axioms:
for all \(A,B\in \Qe\), \(u\in A\) and \((\lambda_1,\lambda_2)\colon A\times B \to \R^2_{>0}\):

$\hat{\mathrm{R}}_0$. \(u \in \hat{R}(A)\) if and only if \(0\in \hat{R}(A-u)\);

$\hat{\mathrm{R}}_1$. \(\hat{R}(\emptyset)=\emptyset\) and if \(A\neq \emptyset\) then \(\hat{R}(A)\neq A\);

$\hat{\mathrm{R}}_2$. if \(u>0\) then \(0\in \hat{R}(\{0,u\})\);

$\hat{\mathrm{R}}_3$. if \(0\in \hat{R}(A\cup\{0\})\) and  \(0\in \hat{R}(B\cup\{0\})\) then
\[0\in \hat{R}(\{\lambda_1(v,w) v +\lambda_2(v,w) w \colon v\in A,w\in B\}\cup \{0\});\]

$\hat{\mathrm{R}}_4$. if \(A\subseteq B\) then \(\hat{R}(A)\subseteq \hat{R}(B)\).

These axioms are very similar to \(\mathrm{R}_0\)–\(\mathrm{R}_4\); the only difference is that they are also imposed for \(A=\emptyset\) and/or \(B=\emptyset\), and that \(\hat{\mathrm{R}}_1\) adds the condition \(\hat{R}(\emptyset)=\emptyset\).

Due to this similarity, it should not come as a surprise that imposing \(\mathrm{R}_0\)–\(\mathrm{R}_4\) on the rejection function \(R_C\) that corresponds to a choice function \(C\) is equivalent to imposing \(\hat{\mathrm{R}}_0\)–\(\hat{\mathrm{R}}_4\) on the extended rejection function \(\hat{R}_{\overline{C}}\) that corresponds to its extension \(\overline{C}\).

For any choice function \(C \colon \Q\to \Q\) and \(k\in \{0,1,...,4\}\) we have that \(R_C\) satisfies  ${\mathrm{R}}_k$ if and only if \(\hat{R}_{\overline{C}}\) satisfies  $\hat{\mathrm{R}}_k$.

It follows immediately from the definitions that for any \(A\in \Q\), \(R_C(A)=\hat{R}_{\overline{C}}(A)\).
Therefore, the implication to the left is immediate.
So now we prove that when at least one of \(A\) and \(B\) are empty, the axioms still hold true.
For \(k=0\), $\hat{\mathrm{R}}_0$ is trivially true for \(A=\emptyset\) as there are no \(u\in A\) then.
For \(k=1\), $\hat{\mathrm{R}}_1$ is true because \(\hat{R}_{\overline{C}}(\emptyset)=\emptyset\) by definition.
\(k=2\) is unrelated to \(A\) and \(B\), so trivially true.
For \(k=3\), $\hat{\mathrm{R}}_3$ is true because if \(A\) and/or \(B\) are empty, then we have \(0\in \hat{R}_{\overline{C}}(\{0\})\) from the antecedent and this is also the implication.
Finally, for \(k=4\), $\hat{\mathrm{R}}_4$ is true because the case where \(A=\emptyset\subseteq B\), for any \(B\in \Qe\), implies that \(\hat{R}_{\overline{C}}(A)=\emptyset \subseteq \hat{R}_{\overline{C}}(B)\) and the case where \(B=\emptyset\) but \(A\neq\emptyset\) is impossible.

Similarly to how we associated a choice function \(C_{\OO}\) with any set of preference orders \(\OO\) in <ref>, we can also associate an extended choice function \(\hat{C}_{\OO}\) with such a set \(\OO\), by letting

\begin{equation}\label{eq:Rdef2}
\hat{C}_{\OO}(A)=\{u\in A\colon (\exists \prec \in \OO)(\forall a \in A) u\not\prec a\}\quad\text{ for all }\quad A\in \Qe.
\end{equation}

We call an extended choice function coherent if and only if it is of this form, with \(\OO\) non-empty.

An extended choice function \(\hat{C}\) is coherent if and only if there is a non-empty set of preference orders \(\OO\subseteq \OOO\) such that \(\hat{C}=\hat{C}_{\OO}\).

A choice function \(C\) is coherent if and only if \(\overline{C}\) is coherent.

The implication to the left is immediate as the definition of coherence is the same for both and \(C\) has a smaller domain than \(\overline{C}\).
For the implication to the right, assume that \(C\) is coherent.
By definition of coherence, there is a set of preference orders \(\OO\) such that \(C=C_{\OO}\).
Then for all \(A\in \Q\) we have \(\overline{C}(A)=C(A)=C_{\OO}(A)=\hat{C}_{\OO}(A)\).
For \(A=\emptyset\) we also have trivially that \(\hat{C}_{\OO}(\emptyset)=\emptyset=\overline{C}(\emptyset)\).

Before we get to the proof of <ref>, we now first use the results of De Bock and De Cooman, 2019 to establish a similar result for extended choice functions.

An extended choice function \(\hat{C}\colon \Qe\to \Qe\) is coherent if and only if its corresponding extended rejection function \(\hat{R}_{\hat{C}}\)  satisfies $\hat{\mathrm{R}}_0$–$\hat{\mathrm{R}}_4$.

Our proof makes use of the following technical lemmata.

Consider any set of preference orders \(\OO\subseteq \OOO\).
Then for all \(A\in \Qe\),
\(
0\in \hat{R}_{\hat{C}_\OO}(A-u) \Leftrightarrow u\in \hat{R}_{\hat{C}_\OO}(A)
\).

\(u\in \hat{R}_{\hat{C}_\OO}(A)\) is by definition equivalent to \((\forall \prec \in \OO)(\exists a \in A) u\prec a\).
This is by <ref> equivalent to \((\forall \prec \in \OO)(\exists v \in A-u) 0\prec v\), which is equivalent to \(0\in \hat{R}_{\hat{C}_\OO}(A-u)\).

Consider an extended rejection function \(\hat{R}\colon \Qe\to\Qe\) that satisfies $\hat{\mathrm{R}}_0$ and let \(K\coloneqq \{ A\in \Qe \colon 0\in \hat{R}(A\cup \{0\})\}\).

\begin{equation}\label{eq:Kdef}
(\forall u \in \V)(\forall A\in \Qe)\; u \in \hat{R}(A\cup \{u\}) \Leftrightarrow A-u\in K.
\end{equation}

\(A-u\in K\) is equivalent to \(0\in \hat{R}((A-u)\cup \{0\})=\hat{R}((A\cup\{u\})-u).\)
By $\hat{\mathrm{R}}_0$ this is equivalent to \(u\in \hat{R}(A\cup\{u\})\).

First we prove the implication to the right.
Consider any non-empty set \(\OO\subseteq \OOO\).
We need to prove that \(\hat{R}_{\hat{C}_\OO}\) satisfies $\hat{\mathrm{R}}_0$ to $\hat{\mathrm{R}}_4$.
Axiom $\hat{\mathrm{R}}_0$ follows immediately from <ref>.

Since \(\hat{R}_{\hat{C}_\OO}(\emptyset)=\emptyset\) follows from Definition (<ref>) and, for any \(A\in \Q\), we have \(\hat{R}_{\hat{C}_\OO}(A)=A\setminus C_{\OO}(A)\), $\hat{\mathrm{R}}_1$ follows immediately from <ref>.
$\hat{\mathrm{R}}_2$ follows immediately from <ref> and Definition (<ref>).
For $\hat{\mathrm{R}}_3$, take any \(\prec\in \OO\).
From Definition (<ref>) and <ref>, we know that there are \(a\in A\) such that \(0\prec a\) and \(b\in B\) such that \(0\prec b\).
Without loss of generality, assume that \( \lambda_1(a,b)>0\).
By <ref>, we have \(0\prec \lambda_1(a,b) a\) and either \(\lambda_2(a,b) b=0\) or \(0\prec \lambda_2(a,b) b\).
In the former case we find that \(0\prec \lambda_1(a,b) a =\lambda_1(a,b) a+\lambda_2(a,b) b\).
In the latter case, it follows from <ref> that  \( \lambda_1(a,b) a\prec\lambda_1(a,b) a+\lambda_2(a,b) b\) and therefore, since \(0\prec \lambda_1(a,b) a\), from <ref> that \(0\prec \lambda_1(a,b) a+\lambda_2(a,b) b\).
Hence, in both cases, this implies that \(0\prec \lambda_1(a,b) a+\lambda_2(a,b) b\).
Since \(\lambda_1(a,b) a+\lambda_2(a,b) b\in \{\lambda_1(v,w) v +\lambda_2(v,w) w \colon v\in A,w\in B\}\cup \{0\}\) and \(\prec \in \OO\) was arbitrary, it therefore follows from Definition (<ref>) that

\[0\in \hat{R}_{\hat{C}_\OO}(\{\lambda_1(v,w) v +\lambda_2(v,w) w \colon v\in A,w\in B\}\cup \{0\}).\]

For $\hat{\mathrm{R}}_4$,
if \(u\in \hat{R}_{\hat{C}_\OO}(A)\), then for all \(\prec \in \OO\) there is some \(a\in A\subseteq B\) such that \(u\prec a\), whence also \(u\in \hat{R}_{\hat{C}_\OO}(B)\).

Next we prove the implication to the left making heavy use of the results in  [De Bock and De Cooman, 2019].
In Ref. [De Bock and De Cooman, 2019], the main tool are so-called `coherent sets of desirable option sets' \(K\subseteq \Qe\), where coherence is characterised through a number of axioms <cit.> that are not important for our purposes.
Since \(\hat{R}_{\hat{C}}\) satisfies $\hat{\mathrm{R}}_0$, we know from <ref> that \(K\coloneqq \{ A\in \Qe \colon 0\in \hat{R}_{\hat{C}}(A\cup \{0\})\}\) satisfies <ref>.
Since \(\hat{R}_{\hat{C}}\) satisfies \(\hat{\mathrm{R}}_0\)–\(\hat{\mathrm{R}}_4\), it therefore  follows from <cit.> that \(K\) is a `coherent set of desirable option sets'.
Therefore, it follows from <cit.> that there is a non-empty set \(\mathscr{D}\subseteq \overline{\mathrm{\mathbf{G}}}\) of coherent sets of desirable options such that

\[K=\bigcap_{G\in \mathscr{D}}\{A\in \Qe\colon A\cap G \neq \emptyset\}.\]

Let \(\OO\coloneqq\{\prec_G\colon G\in \mathscr{D}\}\) and recall from <ref> that this is a set of preference orders, which is furthermore non-empty because \(\mathscr{D}\) is.
Observe also that, for all \(G\in \mathscr{D}\),

\[\{A\in \Qe\colon A\cap G \neq \emptyset\}=\{A\in \Qe \colon (\exists a \in A) 0\prec_G a\}\]
by definition of \(\prec_G\).
We then see that

\begin{align*}
K&=\bigcap_{G\in \mathscr{D}} \{A\in \Qe \colon (\exists a \in A) 0\prec_G a\}\\
&=\{A\in \Qe \colon (\forall G\in \mathscr{D})(\exists a\in A) 0\prec_G a\} \\
&= \{A\in \Qe \colon (\forall \prec \in \OO) (\exists a\in A) 0\prec a\} \\
&= \{A\in \Qe \colon (\forall \prec \in \OO) (\exists a\in A\cup \{0\}) 0\prec a\} \\
&= \{A\in \Qe \colon 0\in \hat{R}_{\hat{C}_\OO}(A\cup \{0\})\},
\end{align*}

where the second to last equality follows from <ref> and the last from <ref>.
By the definition of \(K\), this implies that for any \(A\in \Qe\), \(0\in \hat{R}_{\hat{C}}(A\cup\{0\})\) if and only if \(0\in \hat{R}_{\hat{C}_\OO}(A\cup \{0\}) \).
Now, by $\hat{\mathrm{R}}_0$, the previous observation and <ref> we have for any \(A\in \Qe\) and \(u\in A\) that

\begin{align*}
u\in \hat{R}_{\hat{C}}(A) &\Leftrightarrow 0\in \hat{R}_{\hat{C}}(A-u)\Leftrightarrow 0\in \hat{R}_{\hat{C}}((A-u)\cup\{0\})\\ & \Leftrightarrow 0\in \hat{R}_{\hat{C}_\OO}((A-u)\cup\{0\}) \Leftrightarrow 0\in \hat{R}_{\hat{C}_\OO}(A-u) \Leftrightarrow u\in \hat{R}_{\hat{C}_\OO}(A),
\end{align*}

since \(0\in A-u = (A-u)\cup \{0\}\).
Therefore, \(\hat{R}_{\hat{C}}= \hat{R}_{\hat{C}_\OO}\).

By <ref>, we have that coherence of \(C\) is equivalent to coherence of \(\overline{C}\).
By <ref>, this is equivalent to \(\hat{R}_{\overline{C}}\) satisfying $\hat{\mathrm{R}}_0$–$\hat{\mathrm{R}}_4$.
Finally, by <ref>, this is equivalent to \(R_C\) satisfying ${\mathrm{R}}_0$–${\mathrm{R}}_4$.


# Study of noise in virtual distillation circuits for quantum error mitigation

Pontus Vikstål<EMAIL_ADDRESS>Wallenberg Centre for Quantum Technology,
Department of Microtechnology and Nanoscience, Chalmers University of
Technology, 412 96 Gothenburg, Sweden  Giulia Ferrini  Wallenberg Centre for
Quantum Technology, Department of Microtechnology and Nanoscience, Chalmers
University of Technology, 412 96 Gothenburg, Sweden  Shruti Puri  Department
of Applied Physics, Yale University, New Haven, Connecticut 06511, USA  Yale
Quantum Institute, Yale University, New Haven, Connecticut 06511, USA

###### Abstract

Virtual distillation has been proposed as an error mitigation protocol for
estimating the expectation values of observables in quantum algorithms. It
proceeds by creating a cyclic permutation of $M$ noisy copies of a quantum
state using a sequence of controlled-swap gates. If the noise does not shift
the dominant eigenvector of the density operator away from the ideal state,
then the error in expectation-value estimation can be exponentially reduced
with $M$. In practice, subsequent error-mitigation techniques are required to
suppress the effect of noise in the cyclic permutation circuit itself, leading
to increased experimental complexity. Here, we perform a careful analysis of
noise in the cyclic permutation circuit and find that the estimation of
expectation value of observables diagonal in the computational basis is robust
against dephasing noise. We support the analytical result with numerical
simulations and find that $67\%$ of errors are reduced for $M=2$, with
physical dephasing error probabilities as high as $10\%$. Our results imply
that a broad class of quantum algorithms can be implemented with higher
accuracy in the near-term with qubit platforms where non-dephasing errors are
suppressed, such as superconducting bosonic qubits and Rydberg atoms.

##  I Introduction

Fault-tolerant quantum error correction is necessary for scalable quantum
computation [1], however the associated hardware-performance requirements and
resource overheads are hard to meet with the noisy intermediate-scale quantum
processors available today. Consequently, for near-term applications
alternative techniques to mitigate the effect of noise have been developed.
Some of these techniques are based on scaling noise [2, 3, 4, 5] or learning
about the effect of noise to predict the noise-free behavior of the quantum
protocol [6, 7], while others exploit the symmetry properties of the noise-
free quantum circuit to flag errors [8, 9, 10, 11, 12]. Algorithm- and noise-
specific error-mitigation techniques have also been proposed [13, 14].

Recently an error-mitigation scheme known as virtual distillation, or error
suppression by derangement, has been shown to achieve an exponential
suppression of errors in the estimation of the expectation value of an
observable [15, 16, 17]. The key idea behind this protocol is to compute the
expectation value of an observable by performing measurements on a cyclic-
permutation of $M$ redundant copies of a noisy quantum state. If the effect of
noise is to mix the ideal noise-free state with orthogonal error states, then
symmetries of the cyclic-permutation state suppress the contribution to the
expectation value from the error states exponentially in $M$.

The most straightforward approach to virtual distillation is to prepare the
cyclic-permutation state using an auxiliary qubit and controlled-SWAP (CSWAP)
gates. In practice, this circuit will be prone to errors, limiting the
accuracy of expectation-value estimation without resorting to further noise
mitigation techniques, like zero-noise extrapolation [3, 2]. However, zero-
noise extrapolation not only adds to the sampling cost, but also considerably
increases the circuit complexity as it requires the ability to scale the noise
strength in the quantum circuit either by scaling gate times or by adding more
gates into the circuit [5, 18, 4, 19, 20]. Thus, in this paper we further
investigate the effect of noise in the virtual distillation circuit and
determine conditions under which its faults may be less detrimental, obviating
the need for additional error mitigation.

We consider two commonly studied types of noise: depolarizing and dephasing
noise and find that the mitigated expectation value with virtual distillation
is robust against dephasing noise if the observable is diagonal in the
computational basis, as is typically the case for variational algorithms for
combinatorial optimization [21, 22, 23] or electronic structure calculation of
molecular Hamiltonians [24]. We support our analysis with numerical simulation
of the quantum-approximate optimization algorithm (QAOA) for solving a MaxCut
problem of partitioning the set of vertices in a given graph into two subsets
such that the number of edges shared between the two partitions is maximized.
QAOA is implemented by preparing a variational quantum state using a short-
depth quantum circuit and estimating its energy, i.e. the expectation value of
the Ising-Hamiltonian associated with the MaxCut problem. The parameters in
the circuit are varied until a minimum in the energy landscape is found. In
order to overcome the adverse effects of noise in finding a state that
minimizes the energy, we combine the QAOA protocol with virtual distillation
with $M=2$. We find that the error in estimating the energy with virtual-
distillation with $M=2$ is reduced by $67\%$ when the underlying source of
noise is single-qubit pure dephasing errors at rate of $10\%$, compared to
$20\%$ when the underlying source of noise is single-qubit depolarizing errors
at the same rate. Our findings imply that virtual distillation in a system in
which non-dephasing errors are suppressed compared to dephasing errors is
successful at reducing errors in expectation-value estimation of observables
diagonal in the computational basis without additional error-mitigation
schemes. It is known that such an error channel is relevant for Kerr-cat
qubits in superconducting microwave circuits [25, 26] and Rydberg atomic
qubits [27] not only when the qubits are idle but also during implementation
of Toffoli and controlled-not gates. These two gates can be combined to
implement a CSWAP [28] and thus it is possible to realize robust virtual
distillation in these platforms.

This paper is organized as follows: In Section II we review the virtual
distillation protocol. We analyze the effect of noise in the virtual
distillation circuit on the estimated expectation value in Section III. We
support our analysis with numerical simulations in Section IV and finally, we
give our concluding remarks in Section V.

##  II Virtual Distillation

We begin this section by establishing the notation used throughout this paper,
which is based on Ref. [17]. A boldfaced superscript, for example
$O^{\mathbf{i}}$, will be used to indicate that the operator $O$ acts on the
$i^{\mathrm{th}}$ subsystem. We use superscript with parentheses to indicate
an operator acting on multiple subsystems. For instance, $S^{(M)}$ indicates
that $S$ acts on $M$ subsystems.

Consider the output density operator, $\rho$, of an $N$-qubit noisy quantum
circuit with the spectral decomposition

$\rho=\sum_{k=1}^{d}\lambda_{k}\outerproduct{\psi_{k}}{\psi_{k}}.$ (1)

Here $d=2^{N}$ and $\lambda_{k}$ is the probability that the system is found
in the state $\ket{\psi_{k}}$. We assume, for convenience, that the
probabilities $\lambda_{k}$ are listed in descending order
$\lambda_{1}>\lambda_{2}\ldots>\lambda_{d}$. Exponentiating $\rho$ $M$ times
gives the normalized density operator

$\tilde{\rho}=\frac{\rho^{M}}{\Tr(\rho^{M})}=\frac{\sum_{k=1}^{d}\lambda_{k}^{M}\outerproduct{\psi_{k}}{\psi_{k}}}{\sum_{k=1}^{d}\lambda_{k}^{M}},$
(2)

which approaches the dominant eigenvector $\outerproduct{\psi_{1}}{\psi_{1}}$
exponentially fast with $M$. In virtual distillation, the expectation value of
an observable $O$ is estimated with respect to the exponentiated density
matrix $\tilde{\rho}$,

$\expectationvalue{O}_{\mathrm{mitigated}}:={\Tr(O\tilde{\rho})}=\frac{\Tr(O\rho^{M})}{\Tr(\rho^{M})}.$
(3)

When $\ket{\psi_{1}}$ corresponds to the output of the ideal (noise free)
quantum circuit, then $\expectationvalue{O}_{\mathrm{mitigated}}$ approaches
the ideal expectation value exponentially fast with $M$. This condition is
satisfied when noise in the quantum circuit maps the ideal states to states
that are orthogonal to it, otherwise the dominant eigenvector will drift and
limit the error suppression efficiency [17]. In general, for a multi-qubit
state, both single-qubit depolarizing or dephasing errors can cause drift of
the dominant eigenvector, nevertheless for real world applications this drift
is expected to be small [16, 29, 17], as is also validated by the numerical
simulations in this paper.

Note that, in virtual distillation,
$\expectationvalue{O}_{\mathrm{mitigated}}$ is calculated without explicitly
preparing the state $\tilde{\rho}$, hence the name “virtual”. Instead, virtual
distillation uses $M$ copies of the state $\rho$ together with collective
measurements that only allow symmetric states of the form
$\ket{\psi}\otimes\ket{\psi}\ldots\ket{\psi}$ to contribute to the expectation
value of $O$. More specifically, in Ref. [17] it was shown that Eq. (3) is
equivalent to

$\expectationvalue{O}_{\mathrm{mitigated}}:=\frac{\Tr(O^{(M)}S^{(M)}\rho^{\otimes
M})}{\Tr(S^{(M)}\rho^{\otimes M})},$ (4)

where $O^{(M)}$ is the symmetrized version of the operator $O$,

$O^{(M)}=\frac{1}{M}\sum_{\mathbf{i}=1}^{M}O^{\mathbf{i}},$ (5)

and $S^{(M)}$ is the cyclic shift operator that act on all $M$ subsystems. Its
effect is only to let symmetric states of $\rho^{\otimes M}$ to contribute to
the expectation value of Eq. (4),

$S^{(M)}\ket{\psi_{1}}\otimes\ket{\psi_{2}}\ldots\ket{\psi_{M}}=\ket{\psi_{2}}\otimes\ket{\psi_{3}}\ldots\ket{\psi_{1}}.$
(6)

To measure the observable $S^{(M)}$ in Eq. (4) virtual distillation uses a
procedure similar to the Hadamard test [30]. The procedure begins by preparing
$M$ collective copies of the state $\rho$ together with an auxiliary qubit in
the state $\ket{+}=(\ket{0}+\ket{1})/\sqrt{2}$. Next, a sequence of CSWAP
gates applies $S^{(M)}$ to the $M$ copies of $\rho$ conditioned on the
auxiliary qubit being in state $\ket{1}$. Finally the auxiliary qubit is
measured in the $X$-basis and its expectation value equals to
$\Tr(S^{(M)}\rho^{\otimes M})$, i.e. the denominator of Eq. (4). By
simultaneously also measuring $O^{(M)}$ on the $M$ subsystems $\rho$, the
measurement outcome can be used together with the measurement outcome from the
auxiliary qubit to estimate the numerator $\Tr(O^{(M)}S^{(M)}\rho^{\otimes
M})$ , since $O^{(M)}$ commutes with $S^{(M)}$.

##  III Noise in virtual distillation circuits

We model a noisy gate in the virtual distillation circuit as an ideal gate
followed by independent and identical single-qubit Pauli errors acting on each
qubit participating in the gate. We examine two types of single-qubit noise
channels. The first one is the depolarizing channel which describes a process
where information is completely lost with some probability $\epsilon$, and is
given by [31]

$\Lambda_{\mathrm{dep}}(\rho)=\quantity(1-\epsilon)\rho+\frac{\epsilon}{3}(X\rho
X+Y\rho Y+Z\rho Z),$ (7)

where $\\{X,Y,Z\\}$ are the Pauli operators and $\epsilon$ is the error
probability. The second one is the pure-dephasing channel which describes loss
of phase information with a probability

$\Lambda_{\mathrm{Z}}(\rho)=\quantity(1-\epsilon)\rho+\epsilon Z\rho Z.$ (8)

We use these definitions of the two error channels because their average
channel fidelity is the same. We now evaluate how the above noise channels
affect the mitigation of expectation values as well as their associated
variances.

###  III.1 Noisy mitigated expectation values

In this section we will provide the main results of this paper. For practical
purposes we will limit ourselves to $M=2$. In this case the cyclic shift
operator $S^{(2)}$ factorizes into a tensor product of $N$-SWAPs, and its
controlled version factorizes into a product of $N$ CSWAP-gates. For the
virtual distillation circuit, we assume that a single-qubit noise channel is
applied after every gate to the qubits involved, see FIG. 1.

Figure 1: The virtual distillation circuit for $M=2$ with single-qubit errors
$\Lambda$.

With this, we find that the mitigated expectation value in the presence of
errors in the virtual distillation circuit is given by

$\displaystyle\expectationvalue{O}_{\mathrm{mitigated}}^{\Lambda}=\frac{\Tr(\bar{\Lambda}(O)\rho^{2})}{\Tr(\rho^{2})},$
(9)

where $\bar{\Lambda}=\Lambda\otimes\ldots\otimes\Lambda$ is a tensor product
of $N$ single-qubit error channels and
$\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$. For the full
derivation of Eq. (9) we refer to Appendix A. From Eq. (9) we see that the
influence of errors on the mitigated expectation value will depend on the
observable $O$. Since a general observable on $N$ qubits can be expressed as a
sum of $N$-qubit Pauli strings from the set $\\{I,X,Y,Z\\}^{\otimes N}$, it is
sufficient to consider $O\in\\{I,X,Y,Z\\}^{\otimes N}$. In this case we find
that the mitigated expectation value for the two types of noise are given by

$\displaystyle\expectationvalue{O}_{\mathrm{mitigated}}^{\Lambda_{\mathrm{dep}}}$
$\displaystyle=\quantity(1-\frac{4}{3}\epsilon)^{k}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})},$
(10)
$\displaystyle\expectationvalue{O}_{\mathrm{mitigated}}^{\Lambda_{\mathrm{Z}}}$
$\displaystyle=\quantity(1-2\epsilon)^{k^{\prime}}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})},$
(11)

where $k$ is the number of $\\{X,Y,Z\\}$ Pauli matrices in the tensor product
of $O$, and $k^{\prime}$ is the number of $\\{X,Y\\}$ Pauli matrices in the
tensor product of $O$. Thus we see that errors in the virtual distillation
circuit only attenuate the expectation values. We will study the sampling cost
of this attenuation in the next section. The attenuation can, in principle, be
overcome by error mitigation techniques like polynomial extrapolation.
Importantly, we note that the mitigated expectation value of an observable $O$
that is only a tensor product of Pauli $Z$-operators, so that $k^{\prime}=0$,
is completely immune to pure dephasing in the distillation circuit.

Figure 2: (a) A Graph. (b) The corresponding QAOA circuit for $p=1$ with noise
for the graph in (a). Here
$\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$ is either a
single qubit depolarizing- or dephasing-channel, $\textsc{ZZ}(\alpha)$ on
qubit $j$ and $k$ is $e^{-i\alpha Z_{j}Z_{k}/2}$, and $\textsc{X}(\beta)$ on
qubit $j$ is $e^{-i\beta X_{j}}$. The total gate count in the figure is $6$
single-qubit gates and $8$ two-qubit gates, excluding the error-channels.

###  III.2 Variance of the estimator

We now turn our attention to investigating the sample variance in the
estimation of $\Tr(O^{(2)}S^{(2)}\rho^{\otimes
2})/\allowbreak\Tr(S^{(2)}\rho^{\otimes 2})$. There exists no closed
expression of the sample variance of the quotient between two random
variables, but an approximated one can be obtained by Taylor-expanding the
variance around the mean. If we let $\bar{x}$ denote the sample mean of
$\Tr(O^{(2)}S^{(2)}\rho^{\otimes 2})$, and $\bar{y}$ the sample mean of
$\Tr(S^{(2)}\rho^{\otimes 2})$, then the following unbiased estimator can be
constructed

$\tilde{\theta}:=\frac{\bar{x}}{\bar{y}}.$ (12)

Given sufficiently many samples $R$, the variance of this estimator, for the
noiseless case of FIG. 1, can be approximated as [17]

$\operatorname{Var}\quantity(\mathrm{estim.})\approx\frac{1}{R}\Bigg{(}\frac{\Tr(O^{2}\rho)}{2\Tr(\rho^{2})^{2}}+\frac{\Tr(O\rho)^{2}}{2\Tr(\rho^{2})^{2}}\\\
+\frac{\Tr(O\rho^{2})^{2}}{\Tr(\rho^{2})^{4}}-2\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})^{3}}\Tr(O\rho)\Bigg{)}.$
(13)

Our next step is to find how this variance changes when there is noise in the
virtual distillation circuit. Starting with dephasing errors, we find in
Appendix B that the variance of the estimator is

$\operatorname{Var}_{\Lambda_{\mathrm{Z}}}\quantity(\mathrm{estim.})=\frac{1}{R(1-2\epsilon)^{2N}}\quantity(\frac{\Tr(\bar{\Lambda}(O^{2})\rho)}{2\Tr(\rho^{2})^{2}}+\frac{\Tr(\bar{\Lambda}(O)\rho)^{2}}{2\Tr(\rho^{2})^{2}}+\frac{\Tr(\bar{\Lambda}(O)\rho^{2})^{2}}{\Tr(\rho^{2})^{4}}-2\frac{\Tr(\bar{\Lambda}(O)\rho^{2})}{\Tr(\rho^{2})^{3}}\Tr(\bar{\Lambda}(O)\rho)),$
(14)

which is the same as Eq. (13) but scaled by a factor $(1-2\epsilon)^{-2N}$ and
with $O$ replaced by $\bar{\Lambda}(O)$. In the case of $O$ being a tensor
product of $Z$-operators, $\bar{\Lambda}(O)=O$, the expression reduces to

$\operatorname{Var}_{\Lambda_{\mathrm{Z}}}\quantity(\mathrm{estim.})=\frac{1}{(1-2\epsilon)^{2N}}\operatorname{Var}\quantity(\mathrm{estim.}).$
(15)

This equation shows that $(1-2\epsilon)^{2N}$ extra circuit repetitions are
required for reaching the same level of precision compared to noise free
virtual distillation. However, when $\epsilon\ll 1$ the denominator can be
expanded as $(1-2\epsilon)^{-2N}\approx
1+4N\epsilon+\mathcal{O}(\epsilon^{2})$, which shows that $4N\epsilon$ extra
circuit repetitions are required for small error probabilities. In the
presence of depolarizing errors no simple expression for the variance of the
estimator is found, but we note that for a global depolarizing channel acting
after each CSWAP gate a similar expression to Eq. (15) was found in Ref. [32].

##  IV Numerical Results

To corroborate our analytical results we will now perform numerical
experiments by simulating the quantum approximate optimization algorithm
(QAOA), solving 6-qubit MaxCut problems on 30 randomly generated Erdős–Rényi
graphs with 6-vertices and edge probability of $0.5$.

The MaxCut problem is defined by a graph $G=(V,E)$, where $V$ is the set of
vertices and $E$ is the set of edges. The objective of MaxCut is to partition
the set of vertices into two subsets, such that the number of edges from one
partition to the other is maximum. The problem can be reformulated as finding
the ground state of a Hamiltonian

$H_{\mathrm{MaxCut}}=-\frac{1}{2}\sum_{i,j\in E}(1-Z_{i}Z_{j}),$ (16)

where $Z_{i}$ and $Z_{j}$ are Pauli $Z$ matrices.

We simulate both QAOA and the virtual distillation circuits in the presence of
either single-qubit depolarizing or dephasing errors. Moreover, in order to
also make a comparison that is independent of the state input to the two noisy
virtual distillation circuits, we use the same 30 graphs to create a mixed
state involving the state corresponding to the maximum cut and a thermal state
for each graph. The mixed state is then used as input to the virtual
distillation circuit, and we benchmark its performance in the presence of
either single-qubit depolarizing or dephasing errors.

###  IV.1 Virtual distillation applied to variational states

To approximate the ground state of $H_{\mathrm{MaxCut}}$ with QAOA, the
variational state

$\rho(\boldsymbol{\alpha},\boldsymbol{\beta})=U(\boldsymbol{\alpha},\boldsymbol{\beta})(\outerproduct{+}{+})^{\otimes
N}U^{\dagger}(\boldsymbol{\alpha},\boldsymbol{\beta})$ (17)

is prepared, where

$U(\boldsymbol{\alpha},\boldsymbol{\beta})=\prod_{j=1}^{p}e^{-i\beta_{j}H_{M}}e^{-i\alpha_{j}H_{\mathrm{MaxCut}}}$
(18)

is a unitary operation, $H_{M}=\sum_{i}X_{i}$ is a sum of Pauli $X$ matrices,
$\ket{+}^{\otimes N}$ is a uniform superposition of all computational basis
states, and $\alpha_{j}$, $\beta_{j}$ are $2p$ variational parameters. The
variational parameters are optimized with respect to the expectation value of
the MaxCut Hamiltonian

$C(\boldsymbol{\alpha},\boldsymbol{\beta}):=\Tr(H_{\mathrm{MaxCut}}\rho(\boldsymbol{\alpha},\boldsymbol{\beta})),$
(19)

such that its value is minimized

$(\boldsymbol{\alpha}_{\mathrm{opt}},\boldsymbol{\beta}_{\mathrm{opt}}):=\operatorname*{arg\,min}_{\boldsymbol{\alpha},\boldsymbol{\beta}}C(\boldsymbol{\alpha},\boldsymbol{\beta}).$
(20)

We simulate a noisy QAOA circuit for $p=1$ where $e^{-i\alpha
H_{\mathrm{MaxCut}}}$ is implemented as a product of ZZ–rotations,
$\prod_{j,k\in E}e^{-i\alpha Z_{j}Z_{k}/2}$, followed by single-qubit error-
channels $\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$ acting
on the qubits involved in the gate. The mixer $e^{-i\beta H_{M}}$ is
implemented as single qubit X–rotations, $\prod_{j}e^{-i\beta X_{j}}$,
followed by single-qubit error channel $\Lambda_{\mathrm{dep}}$, to account
for the fact that dephasing errors do not commute with the gate and will thus
depolarize or unbias the noise channel. In practical qubit platforms, two-
qubit gates are generally more noisy than single-qubit gates [33]. Thus, we
reduce the error probability $\epsilon$ by a factor $10$ for the single-qubit
gates. Finally, the same single-qubit error channel $\Lambda$ that is used in
the QAOA circuit is also used in the virtual distillation circuit. We label
the noisy expectation value for error channel
$\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$ obtained using
QAOA without virtual distillation as $C^{\Lambda}(\alpha,\beta)$ and with
virtual distillation as $C^{\Lambda}_{\mathrm{mitigated}}(\alpha,\beta)$.

Figure 3: (a) The distance from the ideal approximation ratio with respect to
the error probability $\epsilon$ averaged over all the instances. (b) The
minimum number of samples $R$ required for the variance of the estimator to be
$\leq 10^{-3}$ averaged over all the instances. Solid (dashed) lines are
obtained with (without) virtual distillation, and blue (red) lines are with
respect to dephasing (depolarizing) errors.

We find the optimal variational parameters for both
$C^{\Lambda}(\alpha,\beta)$ and
$C^{\Lambda}_{\mathrm{mitigated}}(\alpha,\beta)$ for $21$ different error
probabilities $\epsilon$ that are equally spaced between $0$ and $0.1$. To do
this we start by optimizing the energy expectation value for no errors using
brute-force optimization on a $100\times 100$ grid with $\alpha\in[0,\pi]$ and
$\beta\in[0,\pi/2]$ together with an optimization/polishing function that uses
the optimal grid point as an initial guess. The optimal parameters found for
the noiseless QAOA circuit are then used as initial guess to the optimization
function for the first noise iteration, and the initial guess is iteratively
updated for each increasing noise level based on the optimal parameters found
for the previous noise level. This is possible because the optimal parameters
are robust against Pauli channel noise [34, 35].

Figure 4: (a) The approximation ratio with respect to the error probability
$\epsilon$ averaged over all the instances. (b) The minimum number of samples
$R$ required for the variance of the estimator to be $\leq 10^{-3}$ averaged
over all the instances. Solid (dashed) lines are obtained with (without)
virtual distillation, and blue (red) lines are with respect to dephasing
(depolarizing) errors.

After the optimization we compute the approximation ratio by dividing the
expected cuts $C(\alpha_{\mathrm{opt}},\beta_{\mathrm{opt}})$, and
$C^{\Lambda}_{\mathrm{mitigated}}(\alpha_{\mathrm{opt}},\beta_{\mathrm{opt}})$
by the maximum cut $C_{\mathrm{max}}$. FIG. 3a presents the difference between
the ideal (noiseless) approximation ratio and the noisy approximation ratio
averaged over the 30 instances without virtual distillation (unmitigated) for
QAOA with either depolarizing- or dephasing-errors, and with virtual
distillation (mitigated) for either depolarizing- or dephasing-errors. It can
clearly be seen that the best error reduction is obtained for dephasing-errors
in both the QAOA circuit and the virtual distillation circuit. From the figure
the error reduction at $\epsilon=0.1$ with virtual distillation is $67\%$ when
the underlying source of noise is dephasing errors, which shows excellent
error suppression given a high error probability $\epsilon$. For depolarizing
errors the error reduction is $20\%$ with virtual distillation. Interestingly,
we find that dephasing-errors in the QAOA circuit yield a moderately better
mean approximation ratio compared to depolarizing-errors without virtual
distillation, as seen by the blue dashed line being below the red dashed line.
This indicates that QAOA by itself may also be more robust against dephasing
errors.

It is known that single-qubit depolarizing and dephasing error can lead to a
coherent mismatch of the dominant eigenvector with respect to the ideal state
[16, 29, 17]. This means that the fidelity between dominant eigenvector of
$\rho$ and the ideal noiseless state is not unity. In Appendix C we do a
numerical analysis of the coherent mismatch caused by the errors in the QAOA
circuit, and find that it is smaller for dephasing errors in the circuit.

In FIG. 3b we plot the minimum number of samples $R$ required for the variance
of the estimator defined in Eq. (12) to be $\leq 10^{-3}$, averaged over all
the 30 instances. We also do this for the unmitigated variance of the sample
mean, given by

$\operatorname{Var}(\mathrm{estim.})=\frac{1}{2R}\Big{[}\Tr(\rho_{\Lambda}H_{\mathrm{MaxCut}}^{2})\\\
-\Tr(\rho_{\Lambda}H_{\mathrm{MaxCut}})^{2}\Big{]},$ (21)

where $R$ is the total number of samples and $\rho_{\Lambda}$ is the output
from the QAOA circuit with error channel
$\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$. Since the
virtual distillation circuit uses two copies of $\rho_{\Lambda}$ as input but
is only counted as one sample, to make a fair comparison between the
unmitigated and mitigated variances, we scale the variance of the unmitigated
sample mean by $1/2$. Indeed, for a pure state $\rho$, it can be shown that
Eq. (13) is equal to Eq. (21) [17]. From the two solid lines in FIG. 3b, we
see that the mean number of samples grow exponentially with virtual
distillation for both types of errors. However, the difference between the
number of samples needed for dephasing noise is only negligibly ($\ll$ factor
of 10) larger than those needed for depolarizing noise.

###  IV.2 Virtual distillation applied to thermal states

In the previous section, the density matrices that was given as input to the
virtual distillation circuit were different as a result of the error-channel
in the QAOA circuit. In this section, we make a noise agnostic comparison
where the input state is the same to both virtual distillation circuits. The
state that we have chosen is a statistical mixture between the state
corresponding to the maximum-cut and a thermal state. The state corresponding
the maximum-cut is given by the degenerate ground state
$\rho_{\mathrm{ideal}}=\outerproduct{\psi_{\mathrm{GS}}}{\psi_{\mathrm{GS}}}$
of the MaxCut Hamiltonian. We then consider the thermal state

$\rho_{\mathrm{thermal}}=\frac{e^{-\eta H_{\mathrm{MaxCut}}}}{Z(\eta)},$ (22)

where $Z(\eta)=\Tr(e^{-\eta H_{\mathrm{MaxCut}}})$ is the partition function,
and $\eta$ is a constant that is proportional to the inverse temperature. From
this we create an equally mixed state between the two states
$\rho_{\mathrm{ideal}}$ and $\rho_{\mathrm{thermal}}$,

$\rho=\frac{1}{2}\rho_{\mathrm{ideal}}+\frac{1}{2}\rho_{\mathrm{thermal}}.$
(23)

Using $\rho$ as our input to the virtual distillation circuit we perform
virtual distillation on it using the circuit in FIG. 1 with $\Lambda$ as
either $\Lambda_{\mathrm{Z}}$ or $\Lambda_{\mathrm{dep}}$. In the simulations
we choose $\eta=0.1$, and vary the noise-level $\epsilon$ between $0$ and
$0.25$ and calculate the mitigated approximation ratio as well as the number
of repetitions $R$ required for the variance of the estimator to be $\leq
10^{-3}$, averaged over all the 30 instances. FIG. 4a shows the mean
approximation ratio. While the mitigated approximation ratio stays constant
for dephasing errors irrespective of the error probability $\epsilon$, the
mitigated approximation ratio for depolarizing errors quickly decreases and
even becomes lower than the unmitigated approximation ratio as seen by the
solid orange line crossing the black dashed line. This is because, while the
MaxCut Hamiltonian commutes with the error operator $Z$ for the pure dephasing
channel, it does not commute with all the error operators for the depolarizing
one. FIG. 4b shows that the mean number of samples for the mitigated
expectation values grows rapidly with the error probability. However, as
before, the number of samples required when noise is pure dephasing is only
slightly larger ($\lesssim$ factor of 10) than those required when noise is
depolarizing. Note that the two solid lines do not intercept the dashed line
at $\epsilon=0$, which is because the state $\rho$ in Eq. (23) is not a pure
state.

##  V Conclusion

In this work, we have studied the effects of depolarizing- and dephasing-
errors in virtual distillation circuits. We found that depolarizing errors in
the circuit implementation of virtual distillation substantially degrade the
mitigated expectation value. For dephasing errors, we found that the quality
of the mitigated expectation value does not degrade if the observable is a
tensor product of Pauli $Z$ operators. This makes virtual distillation
implemented in a system that strongly favors dephasing errors amenable for
quantum algorithms that use operators diagonal in the computational basis,
such as QAOA. Moreover, trading pure-dephasing noise for depolarizing errors
does not significantly increase the sampling cost.

Our findings suggest that implementing the virtual distillation protocol in a
system that is strongly biased towards dephasing noise thus avoids the need to
rely on other error mitigation techniques in addition to virtual distillation
for mitigating errors. Finally, we mention that bias-preserving CSWAP gates
can be implemented in bosonic cat code systems [36, 25] and Rydberg atoms
[27], which is crucial in order to not unbias or depolarize the noise channel.

###### Acknowledgements.

We thank Oliver Hahn, Timo Hillman, Shahnawaz Ahmed, and Robert Johansson for
fruitful discussions. This work is supported from the Knut and Alice
Wallenberg Foundation through the Wallenberg Center for Quantum Technology
(WACQT). G. F. acknowledges support from the Swedish Research Council
(Vetenskapsrådet) Grant QuACVA. SP was supported by the Air Force Office of
Scientific Research under award number FA9550-21-1-0209.

### CODE AVAILABILITY

The code used for producing the results is made available in Ref. [37]. All
circuit simulations are done using Cirq [38] and Numpy [39]. As optimization
function for the QAOA, minimize implemented in Scipy [40] was used. The random
graphs were generated using NetworkX [41].

## References

  * Lidar and Brun [2013] D. A. Lidar and T. A. Brun,  _Quantum Error Correction_ (Cambridge University Press, 2013).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, “Efficient variational quantum simulator incorporating active error minimization”, Phys. Rev. X 7, 021050 (2017).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, “Error mitigation for short-depth quantum circuits”, Phys. Rev. Lett. 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, “Practical quantum error mitigation for near-future applications”, Phys. Rev. X 8, 031027 (2018).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, “Error mitigation extends the computational reach of a noisy quantum processor”, Nature 567, 491–495 (2019).
  * Strikis  _et al._ [2021] A. Strikis, D. Qin, Y. Chen, S. C. Benjamin, and Y. Li, “Learning-based quantum error mitigation”, PRX Quantum 2, 040330 (2021).
  * Czarnik  _et al._ [2021a] P. Czarnik, A. Arrasmith, P. J. Coles, and L. Cincio, “Error mitigation with clifford quantum-circuit data”, Quantum 5, 592 (2021a).
  * Bonet-Monroig  _et al._ [2018] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. O’Brien, “Low-cost error mitigation by symmetry verification”, Phys. Rev. A 98, 062339 (2018).
  * McArdle  _et al._ [2019] S. McArdle, X. Yuan, and S. Benjamin, “Error-mitigated digital quantum simulation”, Phys. Rev. Lett. 122, 180501 (2019).
  * Sagastizabal  _et al._ [2019] R. Sagastizabal, X. Bonet-Monroig, M. Singh, M. A. Rol, C. Bultink, X. Fu, C. Price, V. Ostroukh, N. Muthusubramanian, A. Bruno, _et al._ , “Experimental error mitigation via symmetry verification in a variational quantum eigensolver”, Phys. Rev. A 100, 010302 (2019).
  * Google AI Quantum and Collaborators  _et al._ [2020] Google AI Quantum and Collaborators, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, _et al._ , “Hartree-fock on a superconducting qubit quantum computer”, Science 369, 1084–1089 (2020).
  * Huggins  _et al._ [2021a] W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, “Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers”, npj Quantum Inf 7, 1–9 (2021a).
  * O’Brien  _et al._ [2021] T. E. O’Brien, S. Polla, N. C. Rubin, W. J. Huggins, S. McArdle, S. Boixo, J. R. McClean, and R. Babbush, “Error mitigation via verified phase estimation”, PRX Quantum 2, 020317 (2021).
  * Maciejewski  _et al._ [2020] F. B. Maciejewski, Z. Zimborás, and M. Oszmaniec, “Mitigation of readout noise in near-term quantum devices by classical post-processing based on detector tomography”, Quantum 4, 257 (2020).
  * Cotler  _et al._ [2019] J. Cotler, S. Choi, A. Lukin, H. Gharibyan, T. Grover, M. E. Tai, M. Rispoli, R. Schittko, P. M. Preiss, A. M. Kaufman, M. Greiner, H. Pichler, and P. Hayden, “Quantum virtual cooling”, Phys. Rev. X 9, 031013 (2019).
  * Koczor [2021a] B. Koczor, “Exponential error suppression for near-term quantum devices”, Phys. Rev. X 11, 031057 (2021a).
  * Huggins  _et al._ [2021b] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, “Virtual distillation for quantum error mitigation”, Phys. Rev. X 11, 041036 (2021b).
  * Dumitrescu  _et al._ [2018] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski, “Cloud quantum computing of an atomic nucleus”, Phys. Rev. Lett. 120, 210501 (2018).
  * Otten and Gray [2019] M. Otten and S. Gray, “Recovering noise-free quantum observables”, Phys. Rev. A 99, 012338 (2019).
  * Giurgica-Tiron  _et al._ [2020] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, “Digital zero noise extrapolation for quantum error mitigation”, 2020 IEEE International Conference on Quantum Computing and Engineering (QCE) , 306–316 (2020).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm”, arXiv:1411.4028 [quant-ph] (2014).
  * Lucas [2014] A. Lucas, “Ising formulations of many NP problems”, Front. Physics 2, 5 (2014).
  * Vikstål  _et al._ [2020] P. Vikstål, M. Grönkvist, M. Svensson, M. Andersson, G. Johansson, and G. Ferrini, “Applying the quantum approximate optimization algorithm to the tail-assignment problem”, Phys. Rev. Applied 14, 034009 (2020).
  * Xia  _et al._ [2018] R. Xia, T. Bian, and S. Kais, “Electronic structure calculations and the ising hamiltonian”, J. Phys. Chem. B 122, 3384-3395 (2018).
  * Puri  _et al._ [2020] S. Puri, L. St-Jean, J. A. Gross, A. Grimm, N. E. Frattini, P. S. Iyer, A. Krishna, S. Touzard, L. Jiang, A. Blais, S. T. Flammia, and S. M. Girvin, “Bias-preserving gates with stabilized cat qubits”, Sci Adv 6, eaay5901 (2020).
  * Grimm  _et al._ [2020] A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, “Stabilization and operation of a kerr-cat qubit”, Nature 584, 205–209 (2020).
  * Cong  _et al._ [2022] I. Cong, H. Levine, A. Keesling, D. Bluvstein, S.-T. Wang, and M. D. Lukin, “Hardware-efficient, fault-tolerant quantum computation with rydberg atoms”, Phys. Rev. X 12, 021049 (2022).
  * Smolin and DiVincenzo [1996] J. A. Smolin and D. P. DiVincenzo, “Five two-bit quantum gates are sufficient to implement the quantum fredkin gate”, Phys. Rev. A 53, 2855 (1996).
  * Koczor [2021b] B. Koczor, “The dominant eigenvector of a noisy quantum state”, New J. Phys. 23, 123047 (2021b).
  * Aharonov  _et al._ [2009] D. Aharonov, V. Jones, and Z. Landau, “A polynomial quantum algorithm for approximating the jones polynomial”, Algorithmica 55, 395–421 (2009).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Czarnik  _et al._ [2021b] P. Czarnik, A. Arrasmith, L. Cincio, and P. J. Coles, “Qubit-efficient exponential suppression of errors”, arXiv:2102.06056 [quant-ph] (2021b).
  * Tannu and Qureshi [2019] S. S. Tannu and M. K. Qureshi, in  _Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems_ (ACM, 2019) pp. 987–999.
  * Xue  _et al._ [2019] C. Xue, Z.-Y. Chen, Y.-C. Wu, and G.-P. Guo, “Effects of quantum noise on quantum approximate optimization algorithm”, arxiv:1909.02196 [quant-ph] (2019).
  * Sharma  _et al._ [2020] K. Sharma, S. Khatri, M. Cerezo, and P. J. Coles, “Noise resilience of variational quantum compiling”, New J. Phys. 22, 043006 (2020).
  * Guillaud and Mirrahimi [2020] J. Guillaud and M. Mirrahimi, “Repetition cat qubits for fault-tolerant quantum computation”, Phys. Rev. X 9, 041053 (2020).
  * Wikståhl [2022] P. Wikståhl, “Wikstahl/error-mitigation: Version 2.0 (v2.0)”, Zenodo, https://doi.org/10.5281/zenodo.7215577 (2022).
  * Cirq Developers [2021] Cirq Developers, “Cirq (v0.12.0)”, Zenodo, https://doi.org/10.5281/zenodo.5182845 (2021).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy”, Nature 585, 357–362 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. aberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”, Nature Methods 17, 261–272 (2020).
  * Hagberg  _et al._ [2008] A. A. Hagberg, D. A. Schult, and P. J. Swart, in  _Proceedings of the 7th Python in Science Conference_ , edited by G. Varoquaux, T. Vaught, and J. Millman (Pasadena, CA USA, 2008) pp. 11–15.
  * van Kempen and van Vliet [2000] G. van Kempen and L. van Vliet, “Mean and variance of ratio estimators used in fluorescence ratio imaging”, Cytometry 39, 300-305 (2000).

##  Appendix A Erroneous mitigated expectation value

In this section we derive an analytical expression of how the mitigated
expectation value is affected by single-qubit depolarizing and dephasing
errors that occur during the virtual distillation circuit. To simplify the
derivation it is useful to use the following convention for the single-qubit
error channels:

$\displaystyle\Lambda_{\mathrm{dep}}(\rho)$
$\displaystyle=\quantity(1-\frac{3}{4}\epsilon)\rho+\frac{3}{4}\epsilon\quantity(X\rho
X+Y\rho Y+Z\rho Z),$ (24) $\displaystyle\Lambda_{\mathrm{Z}}(\rho)$
$\displaystyle=\quantity(1-\frac{\epsilon}{2})\rho+\frac{\epsilon}{2}Z\rho Z,$
(25)

which amounts to substituting $\epsilon\rightarrow 4\epsilon/3$ in the end to
get back to Eq. (7) for the depolarizing channel, and $\epsilon\rightarrow
2\epsilon$ to get back to Eq. (8) for the dephasing channel.

The output from the virtual distillation circuit in FIG. 1 can be written as

$\rho_{\mathrm{out}}=\underbrace{\tilde{\Lambda}_{N}\circ\mathcal{U}_{N}\circ\tilde{\Lambda}_{N-1}\circ\mathcal{U}_{N-1}\ldots\tilde{\Lambda}_{1}\circ\mathcal{U}_{1}}_{\Lambda_{\mathrm{tot}}}\rho_{\mathrm{in}}=\Lambda_{\mathrm{tot}}(\rho_{\mathrm{in}}),$
(26)

where $\rho_{\mathrm{in}}$ is the input state

$\rho_{\mathrm{in}}:=\outerproduct{+}{+}\otimes\rho\otimes\rho,$ (27)

$\mathcal{U}_{i}(\rho_{\mathrm{in}})=\mathrm{CSWAP}_{i}\,\rho_{\mathrm{in}}\,\mathrm{CSWAP}_{i}^{\dagger}$
is the controlled-$\mathrm{SWAP}$ gate that swaps the $i$th qubit of the two
subsystems $\rho^{\otimes 2}$ conditioned on the auxiliary qubit being in the
state $\ket{1}$, and is given by

$\mathrm{CSWAP}_{i}=\outerproduct{0}{0}\otimes I^{\otimes N}\otimes I^{\otimes
N}+\outerproduct{1}{1}\otimes\mathrm{SWAP}_{i},$ (28)

where $I$ is a $2\times 2$ identity matrix; $\tilde{\Lambda}_{i}$ is the
product of the single-qubit error channels that are applied to the qubits
involved in the gate $\mathcal{U}_{i}$,

$\tilde{\Lambda}_{i}=\Lambda\otimes\Lambda_{i}\otimes\Lambda_{i}=\Lambda\otimes\Lambda^{\otimes
2}_{i},$ (29)

with $\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$ being a
single-qubit error channel, and the subscript $i$ of $\Lambda_{i}$ indicates
that $\Lambda$ acts on the $i$th qubit of subsystem $\rho$. The expectation
value of the numerator in Eq. (4) for $M=2$ can be expressed as

$\Tr(X_{\mathrm{aux}}O^{(2)}\Lambda_{\mathrm{tot}}(\rho_{\mathrm{in}}))=\Tr(\Lambda_{\mathrm{tot}}^{\dagger}(X_{\mathrm{aux}}O^{(2)})\rho_{\mathrm{in}}),$
(30)

where

$\Lambda^{\dagger}_{\mathrm{tot}}:=\mathcal{U}^{\dagger}_{1}\circ\tilde{\Lambda}^{\dagger}_{1}\circ\mathcal{U}^{\dagger}_{2}\circ\tilde{\Lambda}^{\dagger}_{2}\ldots\mathcal{U}^{\dagger}_{N}\circ\tilde{\Lambda}^{\dagger}_{N}=\circ_{k=1}^{N}(\mathcal{U}^{\dagger}_{k}\circ\tilde{\Lambda}^{\dagger}_{k})$
(31)

is the adjoint of $\Lambda_{\mathrm{tot}}$, $X_{\mathrm{aux}}$ is the Pauli
$X$ operator on the auxiliary qubit, and $O^{(2)}$ acts on the two subsystems
$\rho$. It can be noted that Eq. (30) is reminiscent of the Heisenberg picture
in which the operators are evolving instead of the quantum states. We can
furthermore drop $\dagger$ from $\mathcal{U}_{k}$ and $\Lambda_{k}$ since they
are both Hermitian quantum maps.

Beginning with the action of $\tilde{\Lambda}_{N}$ in Eq. (31) on
$X_{\mathrm{aux}}O^{(2)}$, the result is

$\displaystyle\tilde{\Lambda}_{N}(X_{\mathrm{aux}}O^{(2)})=\Lambda(X)\otimes\Lambda_{N}^{\otimes
2}(O^{(2)})=(1-\epsilon)X\otimes\Lambda_{N}^{\otimes 2}(O^{(2)}),$ (32)

where we have used the fact that

$\Lambda(X)=\Lambda(\outerproduct{0}{1})+\Lambda(\outerproduct{1}{0})=(1-\epsilon)\outerproduct{0}{1}+(1-\epsilon)\outerproduct{1}{0}=(1-\epsilon)X,$
(33)

for $\Lambda\in\\{\Lambda_{\mathrm{dep}},\Lambda_{\mathrm{Z}}\\}$. Next, Eq.
(32) is followed by the action of $\mathcal{U}_{N}$ which gives

$(\mathcal{U}_{N}\circ\tilde{\Lambda}_{N})(X_{\mathrm{aux}}O^{(2)})=(1-\epsilon)\quantity[\outerproduct{1}{0}\otimes\mathrm{SWAP}_{N}\Lambda_{N}^{\otimes
2}(O^{(2)})+\outerproduct{0}{1}\otimes\Lambda_{N}^{\otimes
2}(O^{(2)})\mathrm{SWAP}_{N}].$ (34)

Repeating this for the remaining terms yields

$\displaystyle\Lambda^{\dagger}_{\mathrm{tot}}(X_{\mathrm{aux}}O^{(2)})$
$\displaystyle=(1-\epsilon)^{N}\quantity[\outerproduct{1}{0}\otimes\prod_{k=1}^{N}\mathrm{SWAP}_{k}\bar{\Lambda}^{\otimes
2}(O^{(2)})+\outerproduct{0}{1}\otimes\bar{\Lambda}^{\otimes
2}(O^{(2)})\prod_{k=1}^{N}\mathrm{SWAP}_{k}]$
$\displaystyle=(1-\epsilon)^{N}X\otimes S^{(2)}\bar{\Lambda}^{\otimes
2}(O^{(2)}),$ (35)

where $\bar{\Lambda}^{\otimes
2}:=\bar{\Lambda}\otimes\bar{\Lambda}=(\Lambda_{1}\circ\ldots\circ\Lambda_{N})\otimes(\Lambda_{1}\circ\ldots\circ\Lambda_{N})=\Lambda\otimes\ldots\otimes\Lambda$
is a tensor product of $2N$ single-qubit error channels $\Lambda$, and where
we have used the fact that $S^{(2)}=\prod_{k=1}^{N}\mathrm{SWAP}_{k}$ and
$[S^{(2)},\bar{\Lambda}^{\otimes 2}(O^{(2)})]=0$ in the last equality. The
expectation value of the numerator is thus calculated to be

$\displaystyle\Tr(\Lambda^{\dagger}_{\mathrm{tot}}(X_{\mathrm{aux}}O^{(2)})\rho_{\mathrm{in}})=(1-\epsilon)^{N}\Tr(X\outerproduct{+}{+})\Tr(\bar{\Lambda}^{\otimes
2}(O^{(2)})S^{(2)}\rho^{\otimes
2})=(1-\epsilon)^{N}\Tr(\bar{\Lambda}(O)\rho^{2}).$ (36)

To calculate the denominator of Eq. (4) for $M=2$, we simply replace $O$ with
the identity $I^{\otimes N}$ for which $\bar{\Lambda}(I^{\otimes
N})=I^{\otimes N}$ because both the depolarizing- and dephasing-channel are
unital maps. Hence the mitigated expectation value with either depolarizing or
dephasing-errors is given by

$\expectationvalue{O}^{\Lambda}_{\mathrm{mitigated}}=\frac{\Tr(\Lambda^{\dagger}_{\mathrm{tot}}(X_{\mathrm{aux}}O^{(2)})\rho_{\mathrm{in}})}{\Tr(\Lambda^{\dagger}_{\mathrm{tot}}(X_{\mathrm{aux}})\rho_{\mathrm{in}})}=\frac{(1-\epsilon)^{N}}{(1-\epsilon)^{N}}\frac{\Tr(\bar{\Lambda}(O)\rho^{2})}{\Tr(\rho^{2})}=\frac{\Tr(\bar{\Lambda}(O)\rho^{2})}{\Tr(\rho^{2})}.$
(37)

We have reached the final expression which shows that single-qubit
depolarizing or dephasing errors that occur during the circuit are equivalent
to a tensor product of $N$-single qubit depolarizing or dephasing channels
$\Lambda$ acting on the observable $O$.

It is useful to consider what happens for the special case when $O$ is a
tensor product of Pauli-operators, e.g. $O=X\otimes Z\otimes I\otimes
Y\ldots$. In this case we have that
$\Lambda_{\mathrm{dep}}(\sigma)=\quantity(1-\epsilon)\sigma$, with
$\sigma\in\\{X,Y,Z\\}$, for the depolarizing channel, and
$\Lambda_{\mathrm{Z}}(\sigma)=\quantity(1-\epsilon)\sigma$, with
$\sigma\in\\{X,Y\\}$, for the dephasing channel. Hence the mitigated
expectation value for the two channels can be written as

$\displaystyle\expectationvalue{O}^{\Lambda_{\mathrm{dep}}}_{\mathrm{mitigated}}$
$\displaystyle=(1-\epsilon)^{k}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})},$ (38)
$\displaystyle\expectationvalue{O}^{\Lambda_{\mathrm{Z}}}_{\mathrm{mitigated}}$
$\displaystyle=(1-\epsilon)^{k^{\prime}}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})},$
(39)

where $k$ is the number of $\\{X,Y,Z\\}$ Pauli-operators in the tensor product
of $O$ and $k^{\prime}$ is the number of $\\{X,Y\\}$ Pauli-operators in the
tensor product of $O$. With the substitution of $\epsilon$ back to the
definitions of Eq. (7) and Eq. (8) the mitigated expectation values becomes

$\displaystyle\expectationvalue{O}^{\Lambda_{\mathrm{dep}}}_{\mathrm{mitigated}}$
$\displaystyle=\quantity(1-\frac{4}{3}\epsilon)^{k}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})},$
(40)
$\displaystyle\expectationvalue{O}^{\Lambda_{\mathrm{Z}}}_{\mathrm{mitigated}}$
$\displaystyle=(1-2\epsilon)^{k^{\prime}}\frac{\Tr(O\rho^{2})}{\Tr(\rho^{2})}.$
(41)

##  Appendix B Variance of the estimator with dephasing errors

In this section, we derive the variance of the estimator in Eq. (12) with
dephasing errors in the virtual distillation circuit. While there exists no
closed-form expression for the variance of the quotient between two random
variables, an approximate expression can be obtained by doing a Taylor
expansion around the mean of the two random variables [42]. By doing so, one
finds that the variance of the estimator is

$\operatorname{Var}(\mathrm{estim.})\approx\frac{1}{R}\quantity(\frac{\operatorname{Var}(X_{\mathrm{aux}}O^{(2)})}{\expectationvalue{X_{\mathrm{aux}}}^{2}}+\frac{\expectationvalue{X_{\mathrm{aux}}O^{(2)}}^{2}}{\expectationvalue{X_{\mathrm{aux}}}^{4}}\operatorname{Var}(X_{\mathrm{aux}})-2\frac{\expectationvalue{O^{(2)}X_{\mathrm{aux}}}}{\expectationvalue{X_{\mathrm{aux}}}^{3}}\operatorname{Cov}(O^{(2)}X_{\mathrm{aux}},X_{\mathrm{aux}})),$
(42)

where $R$ is the number of samples. The terms
$\expectationvalue{X_{\mathrm{aux}}}$ and
$\expectationvalue{O^{(2)}X_{\mathrm{aux}}}$ have been evaluated in Appendix
A, and we are thus left to evaluate the three following terms

$\displaystyle\operatorname{Var}(X_{\mathrm{aux}})$
$\displaystyle=\expectationvalue*{(X_{\mathrm{aux}})^{2}}-\expectationvalue*{X_{\mathrm{aux}}}^{2}=1-\quantity(1-2\epsilon)^{2N}\Tr(\rho^{2})^{2},$
(43) $\displaystyle\operatorname{Var}(O^{(2)}X_{\mathrm{aux}})$
$\displaystyle=\expectationvalue*{(O^{(2)}X_{\mathrm{aux}})^{2}}-\expectationvalue*{O^{(2)}X_{\mathrm{aux}}}^{2}$
$\displaystyle=\expectationvalue*{(O^{(2)})^{2}}-\quantity(1-2\epsilon)^{2N}\Tr(\bar{\Lambda}(O)\rho^{2})^{2}$
$\displaystyle=\frac{1}{4}\expectationvalue{(O^{\mathbf{1}})^{2}}+\frac{1}{4}\expectationvalue{(O^{\mathbf{2}})^{2}}+\frac{1}{2}\expectationvalue{O^{\mathbf{1}}O^{\mathbf{2}}}-\quantity(1-2\epsilon)^{2N}\Tr(\bar{\Lambda}(O)\rho^{2})^{2}$
$\displaystyle=\frac{1}{2}\Tr(\bar{\Lambda}(O^{2})\rho)+\frac{1}{2}\Tr(\bar{\Lambda}(O)\rho)^{2}-\quantity(1-2\epsilon)^{2N}\Tr(\bar{\Lambda}(O)\rho^{2})^{2},$
(44)
$\displaystyle\operatorname{Cov}(O^{(2)}X_{\mathrm{aux}},X_{\mathrm{aux}})$
$\displaystyle=\expectationvalue*{O^{(2)}X_{\mathrm{aux}}^{2}}-\expectationvalue*{O^{(2)}X_{\mathrm{aux}}}\expectationvalue*{X_{\mathrm{aux}}}$
$\displaystyle=\expectationvalue*{O^{(2)}}-(1-2\epsilon)^{2N}\Tr(\bar{\Lambda}(O)\rho^{2})\Tr(\rho^{2})$
$\displaystyle=\Tr(\bar{\Lambda}(O)\rho)-(1-2\epsilon)^{2N}\Tr(\bar{\Lambda}(O)\rho^{2})\Tr(\rho^{2}).$
(45)

Given these calculations we get that the variance of the estimator for
dephasing errors in the virtual distillation circuit is

$\displaystyle\operatorname{Var}(\mathrm{estim.})=\frac{1}{R(1-2\epsilon)^{2N}}\quantity(\frac{\Tr(\bar{\Lambda}(O^{2})\rho)}{2\Tr(\rho^{2})^{2}}+\frac{\Tr(\bar{\Lambda}(O)\rho)^{2}}{2\Tr(\rho^{2})^{2}}+\frac{\Tr(\bar{\Lambda}(O)\rho^{2})^{2}}{\Tr(\rho^{2})^{4}}-2\frac{\Tr(\bar{\Lambda}(O)\rho^{2})}{\Tr(\rho^{2})^{3}}\Tr(\bar{\Lambda}(O)\rho)).$
(46)

Comparing this equation with the expression for the noiseless variance of the
estimator given by Eq. (13), we see that Eq. (46) is scaled by a factor
$1/(1-2\epsilon)^{2N}$ and that $O$ has been replaced by $\bar{\Lambda}(O)$.

##  Appendix C Coherent mismatch

In this section we investigate the coherent mismatch, also known as drift, of
the dominant eigenvector that is caused by the noise in the QAOA circuit.
Since virtual distillation relies on the assumption that the dominant
eigenvector approximates the ideal state, it is of importance to quantify how
good of an approximation this is. The coherent mismatch is defined as [29]

$c:=1-\absolutevalue*{\innerproduct{\psi_{\mathrm{dom}}}{\psi_{\mathrm{ideal}}}}^{2},$
(47)

where $\ket{\psi_{\mathrm{dom}}}$ is the dominant eigenvector of a noisy
quantum state $\rho$ and $\ket{\psi_{\mathrm{ideal}}}$ is with respect to the
ideal (noiseless) state. In FIG. 5a we plot the coherent mismatch, averaged
over all instances, between the ideal noiseless QAOA state, and the dominant
eigenvector of $\rho_{\Lambda}(\alpha_{\mathrm{opt}},\beta_{\mathrm{opt}})$,
where the optimal angles are with respect to
$C_{\mathrm{mitigated}}^{\Lambda}(\alpha,\beta)$. From the figure it can be
seen that the coherent mismatch is smaller for dephasing errors in the QAOA
circuit compared to depolarizing errors. This fact can be explained by
converting an incoherent error to a coherent one. Consider a $Z$ error that
occurs before the $X(\beta)$-gate. By commuting the error through the gate we
get

$X(\beta)Z=e^{-i\beta X}Z=\cos(\beta)Z-i\sin(\beta)XZ=ZX(-\beta),$ (48)

where in the last equality we have used the anti-commutation relation between
$X$ and $Z$, $XZ=-ZX$. Hence, we see that a $Z$ error introduces an incoherent
error that corresponds to an extra $\pi$ rotation of the $X(\beta)$-gate. For
the $ZZ(\alpha)$-gate a $Z$ error acting on any of the two-qubits commute
through the gate, and hence only the $X(\beta)$-gate will be affected by
dephasing errors which causes coherent mismatch. For the depolarizing channel,
however, both $X$ and $Y$ errors will lead to a coherent error of the
$ZZ(\alpha)$-gate which will result in a greater coherent mismatch.
Nevertheless, we see that the coherent mismatch is small for large error
probabilities, less than $2\%$ for depolarizing error.

Figure 5: The coherent mismatch averaged over all instances. The blue line is
the mean coherent mismatch with dephasing errors in the QAOA circuit and the
red line is with depolarizing errors in the QAOA circuit.