Improved Quantum Algorithms for Fidelity Estimation

András Gilyén Alexander Poremba Alfréd Rényi Institute of Mathematics, Hungary. Supported in part by the AWS Center for Quantum Computing. [email protected]Caltech, Pasadena, CA, USA. Funding provided by Samsung Electronics Co., Ltd., for the project “The Computational Power of Sampling on Quantum Computers”. Additional support was provided by the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1733907).Caltech, Pasadena, CA, USA. Partially supported by AFOSR YIP award number number FA9550-16-1-0495 and the Institute for Quantum Information and Matter (an NSF Physics Frontiers Center; NSF Grant PHY-1733907). [email protected]

Abstract

Fidelity is a fundamental measure for the closeness of two quantum states, which is important both from a theoretical and a practical point of view. Yet, in general, it is difficult to give good estimates of fidelity, especially when one works with mixed states over Hilbert spaces of very high dimension. Although, there has been some progress on fidelity estimation, all prior work either requires a large number of identical copies of the relevant states, or relies on unproven heuristics. In this work, we improve on both of these aspects by developing new and efficient quantum algorithms for fidelity estimation with provable performance guarantees in case at least one of the states is approximately low-rank. Our algorithms use advanced quantum linear algebra techniques, such as the quantum singular value transformation, as well as density matrix exponentiation and quantum spectral sampling. As a complementary result, we prove that fidelity estimation to any non-trivial constant additive accuracy is hard in general, by giving a sample complexity lower bound that depends polynomially on the dimension. Moreover, if circuit descriptions for the relevant states are provided, we show that the task is hard for the complexity class called (honest verifier) quantum statistical zero knowledge via a reduction to a closely related result by Watrous.

1 Introduction

Today’s quantum computers suffer from various kinds of incoherent noise (for example, as a result of $T_{1}$ and $T_{2}$ processes), which makes it difficult to use current quantum technologies to their best advantage [NC00]. Characterizing noise in quantum systems is therefore a fundamental problem for quantum computation and quantum information. Since quantum states have a much finer structure than their classical counterparts comprised of probability distributions, this calls for sophisticated distance measures between quantum states, such as trace distance, the Bures metric and fidelity.

Each distance measure captures slightly different aspects of how two quantum states differ. While fidelity is not a metric on the space of density matrices, it stands out by its versatility and applicability, and naturally appears in many practical scenarios. For example, it captures the geometric distance between thermal states of condensed matter systems nearing phase transitions, and can thus provide useful information about the zero temperature phase diagram [ZQWS07, QC09]. In other contexts, the fidelity of quantum states allows one to infer chaotic behavior of thermofield dynamics of many-body quantum systems [XCPdC21].

The fidelity of two positive semi-definite operators $\rho$ and $\sigma$ on a Hilbert space¹¹1In the infinite dimensional setting one should also assume that $\rho$ and $\sigma$ are trace-class operators. To avoid similar difficulties in this paper we restrict our attention to the finite-dimensional case. $\mathcal{H}$ is defined as

\displaystyle F(\rho,\sigma)=\mbox{Tr}\left[\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right].

(1)

The fidelity is symmetric in $\rho$ and $\sigma$ , and for quantum states its value lies between $0$ and $1$ , equalling $1$ if and only if the states are identical. In this work, we are concerned with the problem of estimating the fidelity up to some additive error. In other words, given two density operators $\rho,\sigma\in\mathbb{C}^{d\times d}$ and $\varepsilon\in(0,1)$ , the problem is to output an additive $\varepsilon$ -approximation $\hat{F}(\rho,\sigma)$ such that

\displaystyle F(\rho,\sigma)-\varepsilon\,\leq\,\hat{F}(\rho,\sigma)\,\leq\,F(\rho,\sigma)+\varepsilon.

(2)

We study two input models. In the weaker input model called sampling access, we only assume access to identical independent copies of the states, whereas in the stronger model called purified access, we assume access to quantum circuits $U_{\rho}$ and $U_{\sigma}$ that allow one to prepare a purification of the quantum states. In the latter model, we denote by $T_{\rho}$ and $T_{\sigma}$ the time complexity of preparing the purifications of $\rho$ and $\sigma$ , respectively. Let us also denote by $r\in[d]$ the smallest rank of the two states. Without loss of generality, we can always assume that $r=\mathrm{rk}\left(\rho\right)$ . In general it is computationally difficult to give good estimates of the fidelity $F(\rho,\sigma)$ , especially when one works with mixed states over Hilbert spaces of very high dimension.

In this work, we present new and efficient approximation algorithms for fidelity estimation which have $\operatorname{poly}(r,1/\varepsilon)$ time and sample/query complexity, and far outperform previous algorithms in the literature. As a complementary result, we prove new hardness results and show that the task of approximating fidelity to any non-trivial constant error is hard for the complexity class called (honest verifier) quantum statistical zero knowledge.

1.1 Related work

We now give an overview of approximation algorithms for fidelity estimation. Let us first discuss the setting in which one of the states is pure, which is fairly well-understood as the fidelity reduces to the the simple quantity $F(\lvert\psi\rangle\langle\psi\rvert,\sigma)=\sqrt{\langle\psi\rvert\sigma\lvert\psi\rangle}$ . Buhrmann et al. [BCWdW01] gave an efficient quantum algorithm known as the Swap Test which allows one to obtain an additive $\varepsilon$ -approximation of $\langle\psi\rvert\sigma\lvert\psi\rangle$ given $\tilde{O}(1/\varepsilon^{2})$ indentical copies of $\lvert\psi\rangle$ and $\sigma$ (see also [CSSC18]). Flammia and Liu [FL11] subsequently gave a randomized $\varepsilon$ -approximation algorithm known as direct fidelity estimation which only involves Pauli measurements on $\tilde{O}(d/\varepsilon^{2})$ samples of $\lvert\psi\rangle$ and $\sigma$ . The general task of fidelity estimation in which both density operators are mixed states is far less understood and requires a much more careful approach.

A simple and direct method for estimating mixed state fidelity is through the use of quantum state tomography. O’Donnell and Wright [OW16] showed that, given $O(r\cdot d/\varepsilon^{2})$ copies of $\rho$ , one can obtain an estimate $\hat{\rho}$ such that $\|\hat{\rho}-\rho\|_{2}\leq\varepsilon$ . Using quantum state tomography, one can therefore construct a simple $\operatorname{poly}(d,1/\varepsilon)$ -time approximation algorithm which evaluates the fidelity $F(\hat{\rho},\hat{\sigma})$ directly given in the order of $\operatorname{poly}(d,1/\varepsilon)$ many copies of $\rho$ and $\sigma$ .

Cerezo et al. [CPCC20] later studied the problem of low-rank fidelity estimation on near-term quantum computers via heuristic variational quantum algorithms that require many identical copies of $\rho$ and $\sigma$ . In the same work, the authors also showed that low-rank fidelity estimation is hard for the complexity class called $\mathsf{DQC1}$ , which consists of all problems that can be efficiently solved with bounded error in the one clean-qubit model of quantum computation. Agarwal et al. [ARSW21] recently considered variational algorithms for fidelity estimation in other special cases. As a complementary result, the authors also showed that fidelity estimation in the case where one state is pure and the other is mixed is $\mathsf{BQP}$ -complete. In the meantime Wang et al. [WZC⁺21] proposed a quantum algorithm in the purified access model that utilizes block-encoding techniques and computes an $\varepsilon$ -approximation to $F(\rho,\sigma)$ in time $O\left(\frac{r^{21.5}}{{\varepsilon^{23.5}}}(T_{\rho}+T_{\sigma})\right)$ , where $T_{\rho}$ and $T_{\sigma}$ correspond to the time complexity of purified access, i.e., the complexity of the circuits $U_{\rho}$ and $U_{\sigma}$ preparing purifications of $\rho$ and $\sigma$ respectively. Crucially, the work of Wang et al. [WZC⁺21] does not take the special case into consideration where one of the states is approximately low-rank.

We give a summary of the most relevant results for fidelity estimation in the table below.

Approximation method	Time/query/sample complexity	Assumptions
Quantum state tomography [OW16]	$\operatorname{poly}(d,1/\varepsilon)$	identical copies
Variational fidelity estimation [CPCC20]	N/A (heuristic)	identical copies
Block-encoding algorithm [WZC⁺21]	$\tilde{O}\left(\frac{r^{21.5}}{{\varepsilon^{23.5}}}(T_{\rho}+T_{\sigma})\right)$	purified access
Our block-encoding algorithm (Section 3)	$\tilde{O}\left(\frac{r^{2.5}}{\varepsilon^{5}}(T_{\rho}+T_{\sigma})\right)$	purified access
Our block-encoding algorithm (Section 3)	$\tilde{O}\left(\frac{r^{5.5}}{\varepsilon^{12}}\right)$	identical copies
Our spectral sampling algorithm (Section 4)	$\tilde{O}\left(\frac{r^{10.5}(T_{\rho}+T_{\sigma})}{\varepsilon^{25}\Delta}+\frac{r^{3}T_{\rho}}{\min\{\frac{\varepsilon^{7}}{r^{3}},\Delta\}^{3}}\right)$	purified access, spectrum*
Any (i.e., lower bound) [BOW19, OW15]	$\Omega\left(\frac{r}{\varepsilon}\right)$	identical copies

^∗We remark that our spectral sampling-based algorithm assumes that $\rho$ has a $\Delta$ -gapped spectrum.

Our improved quantum algorithms for fidelity estimation in Section 3 and Section 4 achieve $\operatorname{poly}(r,1/\varepsilon)$ time and sample/query complexity in order to output an $\varepsilon$ -estimate for the fidelity $F(\rho,\sigma)$ . We remark that our algorithms give further significant improvements in the case in which at least one of the states is approximately low-rank through the use of truncation.²²2Just before submitting this manuscript we noticed the concurrent work of Wang et al. [WGL⁺22]. While our block-encoding based algorithm still has a far better complexity, our spectral sampling algorithm is less favourable in the worst case. However, it offers significant improvements in the approximate low-rank regime through the use of truncation.

1.2 Our results

Fidelity estimation with block-encoding based algorithms.

Our first algorithm is based on advanced quantum linear algebra techniques, such as block-encodings and the quantum singular value transformation (QSVT) [GSLW19]. Our algorithm obtains an estimate $F(\rho_{\theta},\sigma)$ , where $\rho_{\theta}$ is a so-called “soft-thresholded” version of $\rho$ in which eigenvalues of $\rho$ below $(1-\delta)\theta$ are completely removed and eigenvalues above $\theta$ are kept intact, while eigenvalues in the interval $[(1-\delta)\theta,\theta]$ are potentially missing or decreased by some amount. In the purified access model our algorithm has the complexity

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}_{\theta}}{\varepsilon^{2}\delta\theta^{\frac{3}{2}}}T_{\rho}+\frac{\mathrm{rk}_{\theta}^{2}}{\varepsilon^{4}\theta^{\frac{1}{2}}}T_{\sigma}\right),

where $\mathrm{rk}_{\theta}$ is any upper bound on the number of eigenvalues of $\rho$ in the interval $[(1-\delta)\theta,1]$ . Since $\mathrm{rk}_{\theta}\leq\frac{1}{(1-\delta)\theta}$ , for any $\delta\in[0,\frac{1}{2}]$ , the above can be always be bounded above by

\widetilde{\mathcal{O}}\left(\frac{T_{\rho}}{\varepsilon^{2}\delta\theta^{\frac{5}{2}}}+\frac{T_{\sigma}}{\varepsilon^{4}\theta^{\frac{5}{2}}}\right).

On the other hand, if we know that $\mathrm{rk}\left(\rho\right)\leq r$ , then choosing $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ and $\delta=\frac{1}{2}$ our algorithm obtains an $\varepsilon$ -precise estimate of $F(\rho,\sigma)$ in complexity

\widetilde{\mathcal{O}}\left(\frac{r^{\frac{5}{2}}}{\varepsilon^{5}}(T_{\rho}+T_{\sigma})\right),

as we show that $|F(\rho,\sigma)-F(\rho_{\theta},\sigma)|\leq\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[0,\theta)}\rho\Pi^{\phantom{\rho}}_{[0,\theta)}\right]}$ in Section 2.7³³3Here $\Pi_{[0,\theta)}$ projects out the eigenvalues of $\rho$ below $\theta$ as defined in Corollary 24. analogously to [CPCC20].

At its core, our algorithm builds on the Hadamard Test which, given a quantum state $\rho$ and a block-encoding of a matrix $A$ , outputs $0$ with probability $\mathrm{Re}(\mbox{Tr}\left[\rho(I+A)/2\right])$ . Denoting by $U\Sigma V^{\dagger}$ a singular value decomposition of $\sqrt{\rho}\sqrt{\sigma}$ , we can then write the fidelity as follows:

\displaystyle F(\rho,\sigma)=\left\lVert\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}=\mbox{Tr}\left[U^{\dagger}\sqrt{\rho}\sqrt{\sigma}V\right]=\mbox{Tr}\left[\sqrt{\rho}\sqrt{\sigma}VU^{\dagger}\right]=\mbox{Tr}\left[\rho(\rho^{-\frac{1}{2}}\sqrt{\sigma}VU^{\dagger})\right].

Therefore, it suffices to construct a (subnormalized) block-encoding of $\rho^{-\frac{1}{2}}\sqrt{\sigma}VU^{\dagger}$ in order to compute an estimate of $F(\rho,\sigma)$ . When working with $\rho_{\theta}$ instead of $\rho$ we can effectively bound $\left\lVert\rho^{-\frac{1}{2}}\right\rVert$ by $\mathcal{O}\left(\theta^{-\frac{1}{2}}\right)$ , so the block-encoding will have subnormlaization $\sim\theta^{-\frac{1}{2}}$ . This means that we can apply approximately $\frac{1}{\sqrt{\theta}\varepsilon}$ -rounds of amplitude amplification to obtain an $\varepsilon$ -precise estimate of $F(\rho_{\theta},\sigma)$ with high probability.

In order to implement a block-encoding of $\rho^{-\frac{1}{2}}\sqrt{\sigma}VU^{\dagger}$ we implement each of $\rho^{-\frac{1}{2}}$ , $\sqrt{\sigma}$ , and $VU^{\dagger}$ as individual block-encodings and simply take the products of the block-encodings, which is a native quantum operation [GSLW19]. Here, a key observation is that the purified access model implies that we also have a unitary block-encoding of $\rho$ and $\sigma$ [GSLW19], so we can obtain block-encodings of $\rho^{-\frac{1}{2}}$ and $\sqrt{\sigma}$ by applying the QSVT on the block-encodings of $\rho$ and $\sigma$ . We can obtain an approximate implementation of $VU^{\dagger}$ by implementing an (approximate) block-encoding of $\sqrt{\sigma}\sqrt{\rho}$ then applying “singular vector transformation” [GSLW19]. Again in order to implement a block-encoding of $\sqrt{\sigma}\sqrt{\rho}$ we can simply implement block-encodings of both $\sqrt{\sigma}$ and $\sqrt{\rho}$ via the QSVT. On a high level, the above describes the essence of our block-encoding-based algorithm. Since the QSVT only allows for polynomial transformations, we need to give appropriate polynomial approximations of $x^{\pm\frac{1}{2}}$ – the details of the approximation error and complexity analysis can be found in Section 3.

In case we only have access to samples of $\rho$ and $\sigma$ , we can still use density matrix exponentiation [LMR14, KLL⁺17] in combination with the QSVT to implement approximate block-encodings of $\rho$ and $\sigma$ , and use essentially the same strategy as described above. However, in order to maintain the required accuracy throughout the circuit our algorithm requires the use of a large number of samples.

Fidelity estimation via quantum spectral sampling.

Our second approximation algorithm for fidelity estimation exploits the fact that it is possible to “sample” from the spectrum of density operators. The main idea is the following. Suppose we wish to estimate the fidelity $F(\rho,\sigma)$ , for density matrices $\rho,\sigma\in\mathbb{C}^{d\times d}$ . Let $r=\mathrm{rk}\left(\rho\right)$ be the rank of $\rho$ , and let $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{r})$ be the spectrum of $\rho$ (with multiplicity). Expanding $\sigma$ in the eigenbasis of $\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert$ , we find that

\displaystyle F(\rho,\sigma)=\mbox{Tr}{\left[\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right]}=\mbox{Tr}{\left[\sqrt{\sum_{i,j\in[r]}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle\,\,\lvert\psi_{i}\rangle\langle\psi_{j}\rvert}\right]}.

(3)

In other words, we can write the fidelity between $\rho$ and $\sigma$ as the quantity $F(\rho,\sigma)=\mbox{Tr}{[\sqrt{\Lambda}]}$ , where $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}\in\mathbb{C}^{d\times d}$ has the following non-trivial entries in the eigenbasis of $\rho$ :

\displaystyle\Lambda_{ij}=\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle,\quad\quad\forall i,j\in\{1,\dots,r\}.

Hence, it suffices to directly compute the matrix elements of $\Lambda$ in order to estimate the fidelity $F(\rho,\sigma)$ . Let us first consider the case of exact fidelity estimation in order to illustrate how our spectral sampling algorithm works. We remark that our spectral sampling algorithm can handle the case when $\rho$ is approximately low-rank using a soft-thresholding approach similar to our block-encoding algorithm.

To estimate the eigenvalues of $\rho$ , we use the idea of quantum spectral sampling first introduced by Lloyd, Mohseni and Rebentrost [LMR14] in the context of quantum principal component analysis, and later extended by Prakash [Pra14]. This subroutine allows us to approximately perform quantum phase estimation on $\rho$ with respect to a unitary $e^{-2\pi i\rho}$ , resulting in the operation

\displaystyle\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert\quad\mapsto\quad\sum_{i=1}^{r}\lambda_{i}\,\lvert\psi_{i}\rangle\langle\psi_{i}\rvert\otimes\lvert\tilde{\lambda}_{i}\rangle\langle\tilde{\lambda}_{i}\rvert.

(4)

By repeatedly performing the operation in (4), we can sample random pairs of eigenstates and eigenvalues $(\lvert\psi_{j}\rangle,\tilde{\lambda}_{j})$ , where $\tilde{\lambda}_{j}\approx\lambda_{j}$ . In order to obtain a full collection of all eigenvalues of $\operatorname{spec}(\lambda)$ , we have to repeat this procedure multiple times, which raises the question: How many repetitions of the quantum spectral sampling procedure are necessary to find a full collection of $r$ distinct eigenvalues? To distinguish between the different eigenvalues, we must assume that $\rho$ has a non-degenerate spectrum $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{r})$ , where each eigenvalue is separated by a gap $\Delta>0$ with

\displaystyle\Delta=\min_{i\in[r-1]}\big{|}\lambda_{i+1}-\lambda_{i}\big{|}.

(5)

In Section 4.6, we give concrete upper bounds on the number of repetitions needed to complete a full collection. In particular, we analyze the non-uniform coupon collector problem which asks how many draws are needed to collect all $r$ eigenvalues, where the $i$ -th eigenvalue is drawn with probability $\lambda_{i}\in(0,1]$ . Denoting by $T_{\operatorname{spec}(\rho)}$ the random variable for the number of draws needed to complete the collection, we have by an identity due to Flajolet et al. [FGT92],

\displaystyle\mathbb{E}[T_{\operatorname{spec}(\rho)}]=\int_{0}^{\infty}\Big{(}1-\prod_{i=1}^{r}(1-e^{-\lambda_{i}t})\Big{)}\,\mathrm{d}t.

(6)

Our first result is a non-trivial upper bound on the average number of draws in the non-uniform coupon collector problem. In Lemma 36, we show that

\displaystyle\mathbb{E}[T_{\operatorname{spec}(\rho)}]\,\leq\,r\cdot H(\operatorname{spec}(\rho))^{-1},

(7)

where $H(x)=r/\sum_{i=1}^{r}x_{i}^{-1}$ is the harmonic mean of $x=(x_{1},\dots,x_{r})$ . This allows us to directly relate the average number of draws necessary to complete the collection to spectral properties of $\rho$ . Unfortunately, our initial bound in (7) is not tight. In particular, for the uniform spectrum $(\frac{1}{r},\dots,\frac{1}{r})$ , our bound tells us that $\mathbb{E}[T_{(\frac{1}{r},\dots,\frac{1}{r})}]\leq r^{2}$ , whereas a well known result on the (standard) uniform coupon collector problem states that the average number of draws is in the order of $\Theta(r\log r)$ .

In order to further improve on the bound in (7), we use a coupling argument which allows us to relate instances of the non-uniform coupon collector problem to worst-case instances of the uniform coupon collector problem. For example, we show in Lemma 38 that, if $\kappa\in(0,1)$ is a lower bound on the smallest eigenvalue of $\rho$ , then it holds that

\mathbb{E}[T_{\operatorname{spec}(\rho)}]\,\leq\,\mathbb{E}[T_{(\frac{1}{m},\dots,\frac{1}{m})}]=\Theta\left(\log(1/\kappa)/\kappa\right),

where we choose $m=\lceil\frac{1}{2\kappa}\rceil$ . In Corollary 39, we generalize the former by introducing a threshold parameter $\theta\in(0,1)$ and only considering eigenvalues of $\rho$ which lie above $\theta$ . This allows us to obtain upper bounds that asymptotically match the bounds for the uniform coupon collector problem.

Going back to fidelity estimation, let us now describe how we can approximate $\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ , which is an additional quantity required to estimate the matrix elements $\Lambda_{ij}$ , for all $i,j\in[r]$ . Our first observation is that the diagonal entries of $\Lambda$ can easily be estimated via the Swap Test introduced by Buhrmann et al. [BCWdW01]. In particular, once we have obtained a complete collection of pairs of eigenstates and eigenvalues $(\lvert\psi_{i}\rangle,\tilde{\lambda}_{i})$ , we can use the Swap Test on input $|\psi_{i}\rangle\!\langle\psi_{i}|$ and $\sigma$ to estimate $\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle$ up to inverse-polynomial (in $\log d$ ) additive error. Unfortunately, estimating the off-diagonal entries of the matrix $\Lambda$ is a lot more involved, since $\sigma_{ij}=\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ is, in general, a complex number which contains both a real and an imaginary part.

One possible solution for estimating $\sigma_{ij}$ is to use density matrix exponentiation introduced by Lloyd, Mohseni and Rebentrost [LMR14] which allows us to approximately implement a unitary $e^{-i\sigma t}$ , for small $t\in(0,1)$ . Let $j\in[r]$ be an index. A second-order Taylor expansion of $e^{-i\sigma t}$ reveals that

\displaystyle e^{-i\sigma t}\lvert\psi_{j}\rangle=\lvert\psi_{j}\rangle-it\,\sigma\lvert\psi_{j}\rangle+O(t^{2}).

(8)

Therefore, for any index $i\in[r]$ , we obtain the following identity,

\displaystyle\langle\psi_{i}\rvert e^{-i\sigma t}\lvert\psi_{j}\rangle=\langle\psi_{i}|\psi_{j}\rangle-it\,\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle+O(t^{2}).

(9)

Re-arranging the quantity in (9), we find that

\displaystyle\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle=\frac{-i}{t}\cdot(\langle\psi_{i}|\psi_{j}\rangle-\langle\psi_{i}\rvert e^{-i\sigma t}\lvert\psi_{j}\rangle)+O(t),

(10)

which yields an approximate formula for $\sigma_{ij}=\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ (up to first order in $t$ ). Notice that the right-hand side of Eq. (10) consists of simple overlaps between pure states. Unfortunately, we cannot apply the Swap Test to estimate the above quantities, since we are dealing with complex-valued inner products. Hence, we must rely on the so-called Hadamard Test due to Aharanov, Jones and Landau [AJL06], which allows one to estimate the real and imaginary parts of $\langle\psi\rvert U\lvert\psi\rangle$ , for a state $\lvert\psi\rangle$ and unitary $U$ . However, in order to estimate the required quantities in Eq. (10), we have to make use of an additional technique. Namely, we use quantum eigenstate filtering due to Lin and Tong [LT20] in order to approximately obtain circuits $U_{i}$ (and $U_{i}^{\dagger}$ ) that prepare (and uncompute) eigenstates $\lvert\psi_{i}\rangle$ of the state $\rho$ via purified access to $U_{\rho}$ , for every index $i\in[r]$ . This allows us to estimate $\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ by instead approximating the following simple quantities up to inverse-polynomial (in $\log d$ ) precision via the Hadamard Test:

	$\displaystyle\langle\psi_{i}\|\psi_{j}\rangle$	$\displaystyle=\mathrm{Re}{\langle\psi_{i}\rvert U_{j}U_{i}^{\dagger}\lvert\psi_{i}\rangle}+\mathrm{Im}{\langle\psi_{i}\rvert U_{j}U_{i}^{\dagger}\lvert\psi_{i}\rangle}\vspace{4mm}$
	$\displaystyle\langle\psi_{i}\rvert e^{-i\sigma t}\lvert\psi_{j}\rangle$	$\displaystyle=\mathrm{Re}{\langle\psi_{i}\rvert e^{-i\sigma t}U_{j}U_{i}^{\dagger}\lvert\psi_{i}\rangle}+\mathrm{Im}{\langle\psi_{i}\rvert e^{-i\sigma t}U_{j}U_{i}^{\dagger}\lvert\psi_{i}\rangle}.$

Another possible – and much more efficient – solution for estimating the quantity $\sigma_{ij}=\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ is the following. Rather than using density matrix exponentiation, our spectral sampling algorithm implements a block-encoding of $\sigma$ which we can easily construct via purified access to the state $\sigma$ . Letting $U$ denote the associated block-encoding unitary, we perform a Hadamard test with respect to $\lvert 0\rangle\langle 0\rvert$ and $U_{j}^{\dagger}UU_{i}$ to directly estimate the real and imaginary parts of

\mathrm{tr}[\lvert 0\rangle\langle 0\rvert U_{i}^{\dagger}UU_{j}]=\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle=\sigma_{ij},\quad\forall i,j\in[r].

Therefore, we can estimate the matrix entries $\Lambda_{ij}=\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle$ , for every pair of indices $i,j\in[r]$ . Denoting our estimate by $\hat{\Lambda}$ , we can then obtain an approximate fidelity estimate by computing $\hat{F}(\rho,\sigma)=\mbox{Tr}{[\sqrt{\hat{\Lambda}_{+}}}]$ , where $\hat{\Lambda}_{+}$ is the projection of $\hat{\Lambda}$ onto the positive semidefinite cone. We show in in Theorem 40 that our spectral sampling algorithm obtains an $\varepsilon$ -estimate $F(\rho_{\theta},\sigma)$ with high probability, where $\rho_{\theta}$ is a “soft-thresholded” version of $\rho$ with $\theta\in(0,1)$ , in time

\tilde{O}\left(\frac{T_{\rho}+T_{\sigma}}{\theta^{10.5}\varepsilon^{4}\Delta}+\frac{T_{\rho}}{\theta^{3}\min\{\theta^{3}\varepsilon,\Delta\}^{3}}\right).

Finally, if we know that the rank of $\rho$ is at most $\mathrm{rk}\left(\rho\right)\leq r$ , then choosing $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ we obtain an $\varepsilon$ -precise estimate of $F(\rho,\sigma)$ with high probability in time

\tilde{O}\left(\frac{r^{10.5}(T_{\rho}+T_{\sigma}\big{)}}{\varepsilon^{25}\Delta}+\frac{r^{3}T_{\rho}}{\min\{\frac{\varepsilon^{7}}{r^{3}},\Delta\}^{3}}\right).

While our spectral-sampling based algorithm for fidelity estimation performs significantly worse than our block-encoding algorithm, it may be easier to implement in certain settings; for example, when it is easy to obtain circuits that prepare the eigenstates of one of the density operators.

1.3 $\mathsf{QSZK_{HV}}$ -hardness of fidelity estimation to any non-trivial accuracy

Now we show that fidelity estimation to any non-trivial fixed precision is $\mathsf{QSZK_{HV}}$ -hard. This provides evidence for the intractability of the problem in general without further assumptions on the states.

Theorem 1 ( $\mathsf{QSZK_{HV}}$ -hardness of non-trivial fidelity estimation).

Consider the following problem: one is given (the description) of two quantum circuits $U,V$ preparing purifications of quantum states $\rho$ and $\sigma$ respectively, and the task is to output a number $\hat{F}(\rho,\sigma)$ such that $|\hat{F}(\rho,\sigma)-F(\rho,\sigma)|\leq\frac{1}{2}-\delta$ . This problem is $\mathsf{QSZK_{HV}}$ -hard for every $\delta\in(0,\frac{1}{2}]$ .

Proof.

By the Fuchs–van de Graaf inequalities, we have

1-F(\rho,\sigma)\leq\frac{1}{2}\left\lVert\rho-\sigma\right\rVert_{1}\leq\sqrt{1-F(\rho,\sigma)^{2}}.

(11)

Suppose we are given quantum circuits preparing purifications of $\rho$ and $\sigma$ and we are promised that either $\frac{1}{2}\left\lVert\rho-\sigma\right\rVert_{1}\leq\varepsilon$ or $\frac{1}{2}\left\lVert\rho-\sigma\right\rVert_{1}\geq 1-\varepsilon$ for some constant $\varepsilon\in(0,\frac{1}{8}]$ . Watrous proved that this problem is $\mathsf{QSZK_{HV}}$ -complete [Wat02]. By Equation 11 $\frac{1}{2}\left\lVert\rho-\sigma\right\rVert_{1}\leq\varepsilon$ implies $1-\varepsilon\leq F(\rho,\sigma)$ , and $1-\varepsilon\leq\frac{1}{2}\left\lVert\rho-\sigma\right\rVert_{1}$ implies $F(\rho,\sigma)\leq\sqrt{2\varepsilon}$ . In particular estimating the fidelity $F(\rho,\sigma)$ to precision $\frac{1}{2}-\sqrt{2\varepsilon}$ solves the distinguishing problem. Substituting $\delta:=\sqrt{2\varepsilon}$ this implies that for every $\delta\in(0,\frac{1}{2}]$ fidelity estimation to precision $\frac{1}{2}-\delta$ is $\mathsf{QSZK_{HV}}$ -hard. Since estimating the fidelity to precision $\frac{1}{2}$ is trivial (taking estimate $\frac{1}{2}$ ), this means that fidelity estimation to any fixed non-trivial accuracy is $\mathsf{QSZK_{HV}}$ -hard in general. ∎

1.4 A sample complexity lower bound for constant precision fidelity estimation

Now we prove that any non-trivial fidelity estimation algorithm must use at least a polynomially large number of copies even if one of the states is known in advance.

As Bădescu, O’Donnell, and Wright [BOW19] pointed out testing closeness with respect to fidelity requires a number of copies of the states proportional to the dimension of the states even if one of the states is a fixed known state, namely the completely mixed state. Their observation follows from a reduction to the earlier results of O’Donnell, and Wright [OW15].

Corollary 2.

Let $\delta\in(0,1/7]$ , and consider the following problem: Given a known quantum state $\sigma$ of rank $r$ and copies of a state $\rho$ with the promise that $\mathrm{rk}\left(\rho\right)\leq r$ , then computing an estimate $\hat{F}(\rho,\sigma)$ such that $|\hat{F}(\rho,\sigma)-F(\rho,\sigma)|\leq\delta$ requires using $\Omega(r/\delta)$ copies in general.

Proof.

We proceed by a reduction to [BOW19, Theorem 1.7], which considers $\sigma$ to be the completely mixed state in dimension $r$ and $\rho$ to be an arbitrary state having half of its eigenvalues $\frac{1-\varepsilon}{r}$ and the other half $\frac{1-\varepsilon}{r}$ . As [BOW19] notes it follows from the results of [OW15] that distinguishing $\sigma$ from such states $\rho$ requires using $\Omega(\frac{r}{\varepsilon^{2}})$ samples for every $\varepsilon\in[0,1]$ . Although they state the result in term of the dimensionality of the Hilbert space, adding extra dimension to the Hilbert space will not reduce the sample complexity, so this result can also be stated in terms of rank.

On the other hand a $\delta:=\frac{F(\sigma,\sigma)-F(\rho,\sigma)}{2}=\frac{1-F(\rho,\sigma)}{2}$ -precise fidelity estimation algorithm can in particular distinguish $\sigma$ and $\rho$ , and thereby any such algorithm must use at least $\Omega(\frac{r}{\varepsilon^{2}})=\Omega(\frac{r}{\delta})$ samples, since $F(\rho,\sigma)=\frac{1}{2}(\sqrt{1+\varepsilon}+\sqrt{1-\varepsilon})\leq 1-\frac{\varepsilon^{2}}{8}$ for every $\varepsilon\in[0,1]\Rightarrow\delta\geq\frac{\varepsilon^{2}}{16}$ .⁴⁴4The Taylor series of $f(x):=\sqrt{1+x}$ is $1+\frac{x}{2}-\frac{x^{2}}{8}+\frac{x^{3}}{16}-\mathcal{O}\left(x^{4}\right)$ . By Lagrange’s remainder theorem we get that for every $x\in(-1,1)$ we have $f(x)=1+\frac{x}{2}-\frac{x^{2}}{8}+\frac{x^{3}}{16}+\frac{f^{(4)}(\eta)}{4!}\eta^{4}$ for some $\eta\in[-|x|,|x|]$ . Since the fourth derivative of $f(x)=\sqrt{1+x}$ is $f^{(4)}(x)=\frac{15}{16(x+1)^{7/2}}$ which is non-negative on $(-1,1)$ we get the inequality $\sqrt{1+x}\leq 1+\frac{x}{2}-\frac{x^{2}}{8}+\frac{x^{3}}{16}$ for every $x\in[-1,1]$ . ∎

Proposition 3.

Let $\delta\in(0,\frac{1}{2})$ , and consider the following problem: Given a known quantum state $\sigma$ of rank $r$ and copies of a state $\rho$ with the promise that $\mathrm{rk}\left(\rho\right)\leq r$ , then computing an estimate $\hat{F}(\rho,\sigma)$ such that $|\hat{F}(\rho,\sigma)-F(\rho,\sigma)|<\frac{1}{2}-\delta$ requires using $\Omega(\delta\sqrt{r})$ copies in general.

Proof.

Consider $\sigma$ to be the uniform distribution over $[r]$ , and the set of states $\rho$ that are uniform over a $\delta^{2}r$ -sized subset of $[r]$ . Then $F(\sigma,\rho)=\delta$ . Suppose we either get copies from $\sigma$ or a uniformly random $\rho$ . Until an element is repeated all that we see are distinct uniformly random elements of $[0,r]$ and in order to find a repetition with non-negligible probability we need to obtain at least $\Omega(\delta r)$ samples. On the other hand estimating fidelity to precision better than $\frac{1}{2}-\frac{\delta}{2}$ can distinguish the two cases. ∎

We think that the above bound can be improved to $\Omega(r)$ using the techniques of [OW15], when one considers the set of states $\rho$ that are uniform over a $\delta^{2}r$ -dimensional subspace of the support of $\sigma$ , however this result does not seem to directly follow from the results of [OW15].

2 Preliminaries

For a matrix $A\in\mathbb{C}^{m\times n}$ and $p\in[1,\infty]$ , we denote by $\left\lVert A\right\rVert_{p}$ the Schatten $p$ -norm, which is the $\ell^{p}$ -norm of the singular values $(\sum_{i}\varsigma_{i}^{p}(A))^{1/p}$ . In particular, we use the notation $\|A\|=\|A\|_{\infty}$ . We recall some useful inequalities [Bha97, Section IV.2]. Hölder’s inequality states that for all $B\in\mathbb{C}^{n\times k}$ and $r,p,q\in[1,\infty]$ such that $\frac{1}{p}+\frac{1}{q}=\frac{1}{r}$ , we have $\left\lVert AB\right\rVert_{r}\leq\left\lVert A\right\rVert_{p}\left\lVert B\right\rVert_{q}$ .⁵⁵5Note that the expression $(\sum_{i}\varsigma_{i}^{p}(A))^{1/p}$ makes sense for every $p>0$ , but will not give a norm for $p\in(0,1)$ (due to violating the triangle inequality). Nevertheless, Hölder’s inequality holds for these quantities as well, which can be formulated as follows [Bha97, Exercise IV.2.7]: $\left\lVert|AB|^{r}\right\rVert_{1}^{\frac{1}{r}}\leq\left\lVert|A|^{p}\right\rVert_{1}^{\frac{1}{p}}\left\lVert|B|^{q}\right\rVert_{1}^{\frac{1}{q}}$ for all $r,p,q\in(0,\infty]$ such that $\frac{1}{p}+\frac{1}{q}=\frac{1}{r}$ , where $|X|=\sqrt{X^{\dagger}X}$ . The trace-norm inequality states that if $n=m$ , then $|\mbox{Tr}(A)|\leq\left\lVert A\right\rVert_{1}$ . For a hermitian matrix $A\in\mathbb{C}^{n\times n}$ with spectral decomposition $A=UDU^{\dagger}$ and $D=\mathrm{diag}(\lambda_{1},\dots,\lambda_{n})$ , we denote by $A_{+}=UD_{+}U^{\dagger}$ and $A_{-}=UD_{-}U^{\dagger}$ the projections onto the positive semidefinite and negative semidefinite cone, respectively, where we let $D_{+}=\mathrm{diag}(\max\{0,\lambda_{1}\},\dots,\max\{0,\lambda_{n}\})$ and $D_{-}=\mathrm{diag}(\min\{0,\lambda_{1}\},\dots,\min\{0,\lambda_{n}\})$ . For an integer $m\in\mathbb{N}$ , we denote by $H_{m}$ the $m$ -th harmonic number given by $H_{m}=1+\frac{1}{2}+\frac{1}{3}+\dots+\frac{1}{m}$ .

For a function $f\colon\mathbb{R}\mapsto\mathbb{C}$ and set $S\subseteq\mathbb{R}$ we use the notation $\left\lVert f\right\rVert_{S}:=\sup_{x\in S}|f(x)|$ .

Definition 4 (Purified access).

Let $\rho\in\mathbb{C}^{d\times d}$ be a density operator. We say that we have purified access to the state $\rho$ if we have access to a unitary $U_{\rho}$ (and its inverse) acting as follows:

U_{\rho}\lvert 0\rangle_{A}\lvert 0\rangle_{B}=\lvert\psi_{\rho}\rangle_{AB}=\sum_{i=1}^{d}\lambda_{i}\lvert\phi_{i}\rangle_{A}\lvert\psi_{i}\rangle_{B},

where $\mbox{Tr}_{A}\left[\lvert\psi_{\rho}\rangle\langle\psi_{\rho}\rvert_{AB}\right]=\rho$ , and where it holds that $\langle\phi_{i}|\phi_{j}\rangle=\langle\psi_{i}|\psi_{j}\rangle=\delta_{ij}$ , for all $i,j\in[d]$ . In this context, we denote by $T_{\rho}$ the time it takes to implement the unitary $U_{\rho}$ .

We use the following result which is a slight adaptation of [HJ90, Fact 7.4.9.2].

Lemma 5 (Projection onto the positive semidefinite cone).

Let $A\in\mathbb{C}^{n\times n}$ be a hermitian matrix with spectral decomposition $A=UDU^{\dagger}$ , where $D=\mathrm{diag}(\lambda_{1},\dots,\lambda_{n})$ , and let $A_{+}=UD_{+}U^{\dagger}$ be the projection onto the positive semidefinite cone with spectrum $D_{+}=\mathrm{diag}(\max\{0,\lambda_{1}\},\dots,\max\{0,\lambda_{n}\})$ . Let ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ be any unitarily invariant norm over $\mathbb{C}^{n\times n}$ . Then, it holds that

A_{+}=\underset{X\succeq 0}{\mathrm{argmin}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A-X\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.

In other words, $A_{+}$ is the closest positive semidefinite matrix to $A$ with respect to the norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ .

2.1 Matrix Arithmetics using blocks of unitaries

In this section we recall some basic results from the generic matrix arithmetic toolbox described in [GSLW18], which is a distilled version of the results of a series of works on quantum algorithms [HHL09, BCC⁺15, CKS17, LC19, vAGGdW20, CGJ19].

First we introduce the definition of block-encoding which the main idea of which is to represents a subnormalized matrix as the upper-left block of a unitary.

U=\left[\begin{array}[]{cc}A/\alpha&.\\ .&.\end{array}\right]\kern 28.45274pt\Longrightarrow\kern 28.45274ptA=\alpha(\langle 0\rvert\otimes I)U(\lvert 0\rangle\otimes I)

Definition 6 (Block-encoding).

Suppose that $A$ is an $s$ -qubit operator, $\alpha,\varepsilon\in\mathbb{R}_{+}$ and $a\in\mathbb{N}$ , then we say that the $(s+a)$ -qubit unitary $U$ is an $(\alpha,a,\varepsilon)$ -block-encoding of $A$ , if

\left\lVert A-\alpha(\langle 0\rvert^{\otimes a}\otimes I)U(\lvert 0\rangle^{\otimes a}\otimes I)\right\rVert\leq\varepsilon.

In case $\alpha=1$ and $\varepsilon=0$ we simply call an $(\alpha,a,\varepsilon)$ -block-encoding a block-encoding (with $a$ ancillas).

There are several ways to construct block-encodings, for a summary of the techniques we refer to [GSLW19]. For our work the most important is the following result due to Low and Chuang [LC19]:

Lemma 7 (Block-encoding of density operators with purified access [GSLW18, Lemma 45]).

Suppose that $\rho$ is an $s$ -qubit density operator and $G$ is an $(a+s)$ -qubit unitary that on the $\lvert 0\rangle\lvert 0\rangle$ input state prepares a purification $\lvert 0\rangle\lvert 0\rangle\rightarrow\lvert\rho\rangle$ , s.t. $\mathrm{Tr}_{a}{|\rho\rangle\!\langle\rho|}=\rho$ . Then $(G^{\dagger}\otimes I_{s})(I_{a}\otimes\mathrm{SWAP}_{s})(G\otimes I_{s})$ is a $(1,a+s,0)$ -block-encoding of $\rho$ .

Block-encodings are convenient to work with, in particular one can efficiently construct linear combinations of block-encodings via the the so-called linear combination of unitaries technique. Moreover, one can also easily form products of block-encodings as follows:

Lemma 8 (Product of block-encoded matrices [GSLW18, Lemma 53]).

If $U$ is an $(\alpha,a,\delta)$ -block-encoding of an $s$ -qubit operator $A$ , and $V$ is an $(\beta,b,\varepsilon)$ -block-encoding of an $s$ -qubit operator $B$ then⁶⁶6The identity operators act on each others ancilla qubits, which is hard to express properly using simple tensor notation, but the reader should read this tensor product this way. $(I_{b}\otimes U)(I_{a}\otimes V)$ is an $(\alpha\beta,a+b,\alpha\varepsilon+\beta\delta)$ -block-encoding of $AB$ .

2.2 The Swap Test

Let us now recall the Swap Test introduced by Buhrmann et al. [BCWdW01]. We remark that detailed circuits for the general case can also be found in the work of Cincio et al. [CSSC18]. Given as input identical copies of density operators $\lvert\psi\rangle\langle\psi\rvert$ and $\sigma$ , we can repeat the following quantum circuit

\Qcircuit@C=2.47mm@R=1.32mm{\lstick{\lvert 0\rangle}&\gate{H}\ctrl{1}\gate{H}\measureD{0/1}\\ \lstick{\lvert\psi\rangle}\qw\multigate{1}{\mathsf{SWAP}}\qw\qw\\ \lstick{\sigma}\qw\ghost{\mathsf{SWAP}}\qw\qw}

$\log\big{(}\frac{2}{\nu}\big{)}/{\eta^{2}}$ many times with parameters $\eta\in(0,1)$ and $\nu\in(0,1)$ to obtain an additive $\eta$ -approximation to the quantity $\langle\psi\rvert\sigma\lvert\psi\rangle$ with probability at least $1-\nu$ .

2.3 The Hadamard test and block-measurements

In this section, we recall the Hadamard Test due to Aharanov, Jones and Landau [AJL06]. Given as input identical copies of a state $\lvert\psi\rangle$ and a unitary $U$ , we can repeat the following circuit

\Qcircuit@C=2.47mm@R=1.32mm{\lstick{\lvert 0\rangle}&\gate{H}\ctrl{1}\gate{H}\measureD{0/1}\\ \lstick{\lvert\psi\rangle}\qw\gate{U}\qw\qw}

$4\log\big{(}\frac{4}{\delta}\big{)}/\xi^{2}$ times with parameters $\xi\in(0,1)$ and $\delta\in(0,1)$ to approximate $z=\langle\psi\rvert U\lvert\psi\rangle$ as $\mathrm{Re}{\langle\psi\rvert U\lvert\psi\rangle}$ and $\mathrm{Im}{\langle\psi\rvert U\lvert\psi\rangle}$ , each within an additive error of $\xi/\sqrt{2}$ with probability $1-\delta/2$ . Note that, to obtain the imaginary part, we have to replace the final $X$ basis measurement of the circuit with a $Y$ basis measurement. By the union bound, we obtain an estimate $\tilde{z}$ with $|\tilde{z}-z|\leq\xi$ with probability $1-\delta$ .

We now generalize the Hadamard Test to expectation values of block-encoded matrices as follows.

Lemma 9 (Hadamard test for estimating the expectation value of block-encoded matrices).

Suppose that $U$ is a block-encoding of $A\in\mathbb{C}^{d\times d}$ . Then, performing the Hadamard test on a quantum state $\rho\in\mathbb{C}^{d\times d}$ with a controlled- $U$ and a subsequent $X$ or $Y$ basis measurement yields outcome $0$ with probability $\mathrm{Re}(\mathrm{Tr}[{\rho(I+A)/2}])$ and $\mathrm{Im}(\mathrm{Tr}[{\rho(I+A)/2}])$ , respectively.

Proof.

The probability of getting outcome $0$ with respect to an $X$ basis measurement is

\left\lVert(\langle+\rvert\otimes I)(\lvert 0\rangle\lvert\psi\rangle+\lvert 1\rangle U\lvert\psi\rangle)/\sqrt{2}\right\rVert^{2}=\frac{1}{4}\left\lVert\lvert\psi\rangle+U\lvert\psi\rangle\right\rVert^{2}=\frac{1}{2}+\frac{1}{2}\mathrm{Re}(\langle\psi|U|\psi\rangle)=\frac{1}{2}+\frac{1}{2}\mathrm{Re}(\mbox{Tr}\left[|\psi\rangle\!\langle\psi|U\right]).

By linearity we see that for a mixed input state $\rho$ the probability of getting outcome $0$ is

\frac{1}{2}+\frac{1}{2}\mathrm{Re}(\mbox{Tr}\left[\rho U\right])=\mathrm{Re}(\mbox{Tr}\left[\rho(I+U)/2\right]).

Now suppose that $A=(\langle 0\rvert\otimes I)U(\lvert 0\rangle\otimes I)$ , and the input state is $|0\rangle\!\langle 0|\otimes\rho$ , then the probability of getting outcome $0$ is

\mathrm{Re}(\mbox{Tr}\left[(|0\rangle\!\langle 0|\otimes\rho)(I+U)/2\right])=\mathrm{Re}(\mbox{Tr}\left[(|0\rangle\!\langle 0|\!\otimes\!I)^{\!2}(|0\rangle\!\langle 0|\otimes\rho)(I+U)/2\right])=\mathrm{Re}(\mbox{Tr}\left[\rho(I+A)/2\right]).

We remark that the imaginary part can be analogously obtained via a final $Y$ basis measurement. ∎

Let us also describe and alternative method for estimating $\mbox{Tr}\left[\rho\cdot A^{\dagger}A\right]$ using a block-encoding of $A$ . Suppose that $U$ is a block-encoding of $A$ , and we apply $U$ to a quantum state $\rho$ and then measure the defining auxiliary qubits of the block-encoding. The probability of finding them in the $\lvert\bar{0}\rangle$ state is

\displaystyle\mbox{Tr}\left[U(|\tilde{0}\rangle\!\langle\tilde{0}|\otimes\rho)U^{\dagger}(|\bar{0}\rangle\!\langle\bar{0}|\otimes I)\right]=\mbox{Tr}\Big{[}\underset{A}{\underbrace{(\langle\bar{0}\rvert\otimes I)U(\lvert\tilde{0}\rangle\otimes I)}}\cdot\rho\cdot\underset{A^{\dagger}}{\underbrace{(\langle\tilde{0}\rvert\otimes I)U^{\dagger}(\lvert\bar{0}\rangle\otimes I)}}\Big{]}=\mbox{Tr}\left[\rho\cdot A^{\dagger}A\right].

Note that $A$ can be a rectangular matrix, which is the case if $\bar{0}$ and $\tilde{0}$ have an unequal number of qubits.

2.4 Quantum Singular Value Transformation

Let $f\colon\mathbb{R}\mapsto\mathbb{C}$ be an odd function, i.e., $f(-x)=-f(x)$ . For a matrix $A$ with singular value decomposition $A=U\Sigma V^{\dagger}$ let us denote by $\mathrm{SV}^{(f)}\!(A)$ the singular value transform of $A$ defined as $\mathrm{SV}^{(f)}\!(A):=Uf(\Sigma)V^{\dagger}$ . We will heavily use the following fundamental result about quantum implementations of the singular value transform:

Theorem 10 (Quantum Singular Value Transformation [GSLW18, Corollary 18 & Lemma 19]).

If $p\in\mathbb{R}[x]$ is an odd polynomial such that $\left\lVert p(x)\right\rVert_{[-1,1]}\leq 1$ and $W$ is a block-encoding of $A$ with $a$ ancillas, then we can implement a block-encoding of $\mathrm{SV}^{(f)}\!(A)$ with $(\deg(p)+1)/2$ uses of $U$ , $(\deg(p)-1)/2$ uses of $U^{\dagger}$ , and $\mathcal{O}\left(a\deg(p)\right)$ other two-qubit gates.

In case we need to approximate the singular value transform by using an approximate block-encoding we have the following robustness result from [GSLW19]:

Lemma 11 (Robustness of singular value transformation [GSLW18, Lemma 23]).

If $p\in\mathbb{C}[x]$ is an odd degree- $d$ polynomial such that $\left\lVert p\right\rVert_{[-1,1]}\leq 1$ , moreover $A,\tilde{A}\in\mathbb{C}^{\tilde{n}\times n}$ are matrices of operator norm at most $1$ , such that

\left\lVert A-\tilde{A}\right\rVert+\left\lVert\frac{A+\tilde{A}}{2}\right\rVert^{2}\leq 1,

(12)

then we have that

\left\lVert\mathrm{SV}^{(p)}\!(A)-\mathrm{SV}^{(p)}\!(\tilde{A})\right\rVert\leq d\sqrt{\frac{2}{1-\left\lVert\frac{A+\tilde{A}}{2}\right\rVert^{2}}}\left\lVert A-\tilde{A}\right\rVert.

The above lemma is proven in [GSLW18] by showing that there are block-encodings $U,\tilde{U}\in\mathbb{C}^{4(\tilde{n}+n)\times 4(\tilde{n}+n)}$ of $A$ and $\tilde{A}$ respectively so that

\left\lVert U-\tilde{U}\right\rVert\leq\sqrt{\frac{2}{1-\left\lVert\frac{A+\tilde{A}}{2}\right\rVert^{2}}}\left\lVert A-\tilde{A}\right\rVert.

(13)

Using singular value decomposition one can show that if $U,W\in\mathbb{C}^{d\times d}$ are block-encodings of $A\in\mathbb{C}^{\tilde{n}\times n}$ , then there are unitaries $V\in\mathbb{C}^{d-n\times d-n}$ , $\tilde{V}\in\mathbb{C}^{d-\tilde{n}\times d-\tilde{n}}$ such that $(I_{\tilde{n}}\oplus\tilde{V})U(I_{n}\oplus V^{\dagger})=W$ .⁷⁷7In fact this is implicit in the proof of the main results of [GSLW18].

Corollary 12.

For every block-encoding $U\in\mathbb{C}^{d\times d}$ of $A\in\mathbb{C}^{\tilde{n}\times n}$ with $d\geq 4(\tilde{n}+n)$ and for any $\tilde{A}\in\mathbb{C}^{\tilde{n}\times n}$ satisfying Equation 12 we have a block-encoding $\tilde{U}\in\mathbb{C}^{d\times d}$ of $\tilde{A}$ satisfying Equation 13.

Corollary 13 (Error propagation in block-encodings).

Let $d\geq 4(\tilde{n}+n)$ ; let $W(U)$ be a quantum circuit that uses $U$ and $U^{\dagger}$ a total of $T$ -times. If $W(U)$ implements a block-encoding of $B$ for every block-encoding $U^{d\times d}$ of $A\in\mathbb{C}^{\tilde{n}\times n}$ , then $W(\tilde{U})$ implements a block-encoding of $\tilde{B}$ such that if $\tilde{U}$ is a block-encoding of $\tilde{A}\in\mathbb{C}^{\tilde{n}\times n}$ satisfying Equation 12 then

\left\lVert B-\tilde{B}\right\rVert\leq T\sqrt{\frac{2}{1-\left\lVert\frac{A+\tilde{A}}{2}\right\rVert^{2}}}\left\lVert A-\tilde{A}\right\rVert.

(14)

Proof.

Take some $A\in\mathbb{C}^{\tilde{n}\times n}$ , $d\geq 4(\tilde{n}+n)$ , and a block-encoding $\tilde{U}$ of some $\tilde{A}\in\mathbb{C}^{\tilde{n}\times n}$ satisfying Equation 12. By Corollary 12 there exists some $U$ block-encoding of $A$ such that Equation 13 holds. Then $W(U)$ is a block-encoding of $B$ , and $\left\lVert W(U)-W(\tilde{U})\right\rVert\leq T\left\lVert U-\tilde{U}\right\rVert\leq T\sqrt{\frac{2}{1-\left\lVert\frac{A+\tilde{A}}{2}\right\rVert^{2}}}\left\lVert A-\tilde{A}\right\rVert$ , and therefore Equation 14 holds as well. ∎

2.5 Low-degree polynomial approximations

As Theorem 10 suggests in order to optimize our algorithm we will need low-degree polynomial approximations of various functions. We list some polynomial approximation results that we need in our complexity analysis.

Later we will use a low-degree polynomial $p$ approximating the sign function $\mathrm{sgn}\mapsto\{-1,0,1\}$ . To construct that we invoke a result of Low and Chuang [LC17, Corollary 6] about constructive polynomial approximations of the sign function – the error of the optimal approximation, studied by Eremenko and Yuditskii [EY07], achieves similar scaling but is non-constructive.

Lemma 14 (Polynomial approximations of the sign function).

For all $\delta>0$ , $\varepsilon\in(0,1/2)$ there exists an efficiently computable odd polynomial $p\in\mathbb{R}[x]$ of degree $n=\mathcal{O}\left(\frac{\log(1/\varepsilon)}{\delta}\right)$ , such that

•

for all $x\in[-2,2]\colon|p(x)|\leq 1$ , and
•

for all $x\in[-2,2]\setminus(-\delta,\delta)\colon|p(x)-\mathrm{sgn}(x)|\leq\varepsilon$ .

Corollary 15 (Polynomial approximations of rectangle functions [GSLW19, Corollary 16]).

Let $\delta,\varepsilon$ $\in(0,\frac{1}{2})$ and $t\in[-1,1]$ . There exists an even polynomial $P^{\prime}\in\mathbb{R}[x]$ of degree $\mathcal{O}\left(\log(\frac{1}{\varepsilon})/\delta\right)$ , such that $|P^{\prime}(x)|\leq 1$ for all $x\in[-1,1]$ , and

\left\{\begin{array}[]{rcl}P^{\prime}(x)\in&[0,\varepsilon]&\text{for all }x\in[-1,-t-\delta]\cup[t+\delta,1]\text{, and}\\ P^{\prime}(x)\in&[1-\varepsilon,1]&\text{for all }x\in[-t+\delta,t-\delta].\end{array}\right.

(15)

Lemma 16 (Polynomial approximations of negative power functions [Gil19, Cor. 3.4.13]).

Let $\delta,\varepsilon\in(0,\frac{1}{2}]$ , $c>0$ and let $f(x):=\frac{\delta^{c}}{2}x^{-c}$ , then there exist even/odd polynomials $P,P^{\prime}\in\mathbb{R}[x]$ , such that $\left\lVert P-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P\right\rVert_{[-1,1]}\leq 1$ and similarly $\left\lVert P^{\prime}-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P^{\prime}\right\rVert_{[-1,1]}\leq 1$ , moreover the degrees of the polynomials are $\mathcal{O}\left(\frac{\max[1,c]}{\delta}\log\left(\frac{1}{\varepsilon}\right)\right)$ .

Lemma 17 (Polynomial approximations of positive power functions [Gil19, Cor. 3.4.14]).

Let $\delta,\varepsilon\in(0,\frac{1}{2}]$ , $c\in(0,1]$ and let $f(x):=\frac{1}{2}x^{c}$ , then there exist even/odd polynomials $P,P^{\prime}\in\mathbb{R}[x]$ , such that $\left\lVert P-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P\right\rVert_{[-1,1]}\leq 1$ and similarly $\left\lVert P^{\prime}-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P^{\prime}\right\rVert_{[-1,1]}\leq 1$ , moreover the degree of the polynomials are $\mathcal{O}\left(\frac{1}{\delta}\log\left(\frac{1}{\varepsilon}\right)\right)$ .

Corollary 18 (Polynomial approximations of positive power functions).

Let $\delta,\varepsilon\in(0,\frac{1}{2}]$ , $c\in(0,1]$ and let $f(x):=\frac{1}{2}x^{c}$ , then there exist even/odd polynomials $P,P^{\prime}\in\mathbb{R}[x]$ , such that $\left\lVert P-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P\right\rVert_{[-1,1]}\leq 1$ , $\forall x\in[-1,1]\colon|P(x)|\leq f(|x|)+\varepsilon$ , and similarly $\left\lVert P^{\prime}-f\right\rVert_{[\delta,1]}\leq\varepsilon$ , $\left\lVert P^{\prime}\right\rVert_{[-1,1]}\leq 1$ , $\forall x\in[-1,1]\colon|P^{\prime}(x)|\leq f(|x|)+\varepsilon$ , moreover the degrees of the polynomials are $\mathcal{O}\left(\frac{\max[1,c]}{\delta}\log\left(\frac{1}{\varepsilon}\right)\right)$ .

Proof.

Choose $\delta^{\prime}:=\delta/2$ and $\varepsilon^{\prime}:=\varepsilon/2$ in Lemma 17. Then take a polynomial given by Corollary 15 for parameters $\delta^{\prime}:=\delta/4$ , $t:=\delta/4$ and $\varepsilon^{\prime}:=\varepsilon/2$ . The product of the corresponding polynomial satisfies the desired properties. ∎

2.6 Density matrix exponentiation and block-encodings

In case we do not have purified access to the density operators, but can only get independent copies we do not have a direct analog of Lemma 7. Instead we will rely on the technique of density matrix exponentiation [LMR14]. For this we invoke the following form of the result:

Theorem 19 (Density matrix exponentiation [KLL⁺17, Theorem 5 & Theorem 20]).

For an unkown quantum state $\rho\in\mathbb{C}^{2^{q}\times 2^{q}}$ the sample complexity of implementing the controlled- $e^{it\rho}$ unitary to diamond-norm error $\delta$ is $\Theta(\frac{t^{2}}{\delta})$ (for the lower-bound one needs to assume $\delta\leq\frac{1}{6}\min(1,\frac{|t|}{\pi})$ ). Moreover, the implementation uses $\mathcal{O}\left(q\cdot\frac{t^{2}}{\delta}\right)$ two-qubit quantum gates.

As pointed out in [GLM⁺20], using the results of [GSLW19], in particular the result about taking the logarithm of a unitary [GSLW18, Corollary 71] one can implement a block-encoding of $\rho$ also in the sampling access model. To see this we first recall the following result:

Lemma 20 (Implementing the logarithm of unitaries [GSLW18, Corollary 71]).

Suppose that $U=e^{iH}$ , where $H$ is a Hamiltonian of norm at most $\frac{1}{2}$ . Let $\varepsilon\in(0,\frac{1}{2}]$ , then we can implement a $(\frac{2}{\pi},2,\varepsilon)$ -block-encoding of $H$ with $\mathcal{O}\left(\log\left(\frac{1}{\varepsilon}\right)\right)$ uses of controlled- $U$ and its inverse, using $\mathcal{O}\left(\log\left(\frac{1}{\varepsilon}\right)\right)$ two-qubit gates and using a single ancilla qubit.

Applying this result to the controlled- $e^{\pm i\frac{\rho}{2}}$ evolution given by Theorem 19 we get the following corollary of Theorem 19:

Corollary 21 (Sampling to block-encoding).

For an unknown quantum state $\rho\in\mathbb{C}^{2^{q}\times 2^{q}}$ we can implement a quantum operation (quantum channel) that is $\delta$ -close in the diamond-norm to an $(\frac{4}{\pi},2,\varepsilon)$ -block-encoding of $\rho$ by using $\mathcal{O}\left(\frac{\log^{2}(1/\varepsilon)}{\delta}\right)$ samples of $\rho$ and $\mathcal{O}\left(q\cdot\frac{\log^{2}(1/\varepsilon)}{\delta}\right)$ two-qubit quantum gates. In particular we can implement a quantum operation (quantum channel) that is $\delta$ -close in the diamond-norm to an $(\frac{4}{\pi},3,0)$ -block-encoding of $\rho$ by using $\mathcal{O}\left(\frac{\log^{2}(1/\delta)}{\delta}\right)$ samples of $\rho$ and $\mathcal{O}\left(q\cdot\frac{\log^{2}(1/\delta)}{\delta}\right)$ two-qubit quantum gates.

Proof.

Due to Lemma 20 we can implement an $(\frac{4}{\pi},2,\varepsilon)$ -block-encoding of $\rho$ with $\mathcal{O}\left(\log(1/\varepsilon)\right)$ uses of controlled- $e^{\pm i\frac{\rho}{2}}$ unitary. Due to Theorem 19 we can approximate each application of the controlled- $e^{\pm i\frac{\rho}{2}}$ unitary to diamond error $\mathcal{O}\left(\delta/\log(1/\varepsilon)\right)$ using $\mathcal{O}\left(\frac{\log(1/\varepsilon)}{\delta}\right)$ samples. This amounts to an overall number of $\mathcal{O}\left(\frac{\log^{2}(1/\varepsilon)}{\delta}\right)$ samples of $\rho$ . By the Lemma 20 and Theorem 19 the gate complexity of this implementation is $\mathcal{O}\left(q\cdot\frac{\log^{2}(1/\varepsilon)}{\delta}\right)$ .

On the other hand a $(\frac{4}{\pi},2,\varepsilon)$ -block-encoding also trivially gives rise to a $(\frac{4}{\pi},3,\varepsilon)$ -block-encoding by simply adding an extra qubit with the identity operation on it. Furthermore, due to Corollary 13 a unitary $(\frac{4}{\pi},3,\varepsilon)$ -block-encoding of $\rho$ is $\mathcal{O}\left(\varepsilon\right)$ -close to a unitary $(\frac{4}{\pi},3,0)$ -block-encoding of $\rho$ , and by definition these unitaries are also $\mathcal{O}\left(\varepsilon\right)$ -close in the diamond norm. So choosing $\delta\leftarrow\frac{\delta^{\prime}}{2}$ and $\varepsilon\leftarrow\mathcal{O}\left(\delta^{\prime}\right)$ the above construction gives a $\delta^{\prime}$ -precise implementation of a $(\frac{4}{\pi},3,0)$ -block-encoding of $\rho$ in the diamond norm. ∎

2.7 Truncated fidelity bounds

In general it is difficult to work with density operators that have a large number of tiny eigenvalues that all together represent a significant contribution to the trace. On the other hand, if we filter out small eigenvalues then the problem becomes tractable. Since in general we can only apply soft versions of filtering we need to understand how big is the inaccuracy introduced by such soft truncation. Therefore we devise some slight generalizations of the truncation bounds from [CPCC20, Section 2].

Lemma 22 (Monotonicity of fidelity).

Let $0\preceq A\preceq B$ such that $A,B$ , and $\rho$ commute with each other, then

\displaystyle F(A\rho A,\sigma)\leq F(B\rho B,\sigma).

Proof.

Since $A\rho A\preceq B\rho B$ , and the $\sqrt{\cdot}$ function is operator monotone [HP14, Chapter 4.1] we have

F(A\rho A,\sigma)=\mbox{Tr}\left[\sqrt{\sqrt{\sigma}A\rho A\sqrt{\sigma}}\right]\leq\mbox{Tr}\left[\sqrt{\sqrt{\sigma}B\rho B\sqrt{\sigma}}\right]=F(B\rho B,\sigma).\qed

Lemma 23 (Hard truncation bounds).

Let $\rho,\sigma\succeq 0$ and $\Pi$ be an orthogonal projector that commutes with $\rho$ , then

\displaystyle F(\rho,\sigma)\leq F(\Pi\rho\Pi,\sigma)+\sqrt{\mbox{Tr}\left[(I-\Pi)\rho(I-\Pi)\right]}\sqrt{\mbox{Tr}\left[(I-\Pi)\sigma(I-\Pi)\right]}.

(16)

Proof.

$\displaystyle F(\rho,\sigma)=\left\lVert\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}$	$\displaystyle\leq\left\lVert\sqrt{\rho}\Pi\sqrt{\sigma}\right\rVert_{1}+\left\lVert\sqrt{\rho}(I-\Pi)\sqrt{\sigma}\right\rVert_{1}$	(by the triangle inequality)
	$\displaystyle=\left\lVert\sqrt{\Pi\rho\Pi}\sqrt{\sigma}\right\rVert_{1}+\left\lVert\sqrt{\rho}(I-\Pi)^{2}\sqrt{\sigma}\right\rVert_{1}$	( $\Pi\rho=\rho\Pi$ and $(I-\Pi)=(I-\Pi)^{2}$ )
	$\displaystyle\leq\left\lVert\sqrt{\Pi\rho\Pi}\sqrt{\sigma}\right\rVert_{1}+\\|\sqrt{\rho}(I-\Pi)\\|_{2}\cdot\\|\sqrt{\sigma}(I-\Pi)\\|_{2}$	(by Hölder’s inequaltiy)
	$\displaystyle=F(\Pi\rho\Pi,\sigma)+\sqrt{\mbox{Tr}\left[(I-\Pi)\rho(I-\Pi)\right]}\cdot\sqrt{\mbox{Tr}\left[(I-\Pi)\sigma(I-\Pi)\right]}.\qed$

Corollary 24 (Soft truncation bounds).

Let $\rho,\sigma$ be quantum states and $\Pi^{\phantom{\rho}}_{I}$ be the orthogonal projector to the subspace spanned by eigenvectors of $\rho$ with eigenvalues in $I$ . Let $0\leq\alpha\leq\beta$ and $f\colon\mathbb{R}\mapsto[0,1]$ be such that for all $x<\alpha$ we have $f(x)=0$ and for all $x\geq\beta$ we have $f(x)=1$ , then

F(\Pi^{\phantom{\rho}}_{[\beta,1)}\rho\Pi^{\phantom{\rho}}_{[\beta,1)},\sigma)\leq F(f(\rho)\cdot\rho,\sigma)\leq F(\Pi^{\phantom{\rho}}_{[\alpha,1)}\rho\Pi^{\phantom{\rho}}_{[\alpha,1)},\sigma)\leq F(\rho,\sigma),

(17)

and

F(\rho,\sigma)-F(f(\rho)\cdot\rho,\sigma)\leq\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[0,\beta)}\rho\Pi^{\phantom{\rho}}_{[0,\beta)}\right]}\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[0,\beta)}\sigma\Pi^{\phantom{\rho}}_{[0,\beta)}\right]}.

(18)

Proof.

Let $P:=\sqrt{f(\rho)}$ , then $\Pi^{\phantom{\rho}}_{[\beta,1)}\preceq P\preceq\Pi^{\phantom{\rho}}_{[\alpha,1)}\preceq I$ so Equation 17 follows from Lemma 22 and Equation 18 follows from Equation 17 and Lemma 23. ∎

Replacing $\rho$ in Lemma 23 by $\Pi_{\alpha}^{\phantom{\rho}}\rho\Pi_{\alpha}^{\phantom{\rho}}$ (letting $\Pi_{\theta}^{\phantom{\rho}}:=\Pi_{[\theta,1)}^{\phantom{\rho}}$ , $\alpha\leq\beta$ ) we also get the following bound:

	$\displaystyle F(\Pi_{\alpha}^{\phantom{\rho}}\rho\Pi_{\alpha}^{\phantom{\rho}},\sigma)-F(\Pi_{\beta}^{\phantom{\rho}}\rho\Pi_{\beta}^{\phantom{\rho}},\sigma)$	$\displaystyle=F(\Pi_{\alpha}^{\phantom{\rho}}\rho\Pi_{\alpha}^{\phantom{\rho}},\Pi_{\alpha}^{\phantom{\rho}}\sigma\Pi_{\alpha}^{\phantom{\rho}})-F(\Pi_{\alpha}^{\phantom{\rho}}\Pi_{\beta}^{\phantom{\rho}}\rho\Pi_{\beta}^{\phantom{\rho}}\Pi_{\alpha}^{\phantom{\rho}},\Pi_{\alpha}^{\phantom{\rho}}\sigma\Pi_{\alpha}^{\phantom{\rho}})$
		$\displaystyle\leq\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[\alpha,\beta)}\rho\Pi^{\phantom{\rho}}_{[\alpha,\beta)}\right]}\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[\alpha,\beta)}\sigma\Pi^{\phantom{\rho}}_{[\alpha,\beta)}\right]},$		(19)

where in the equality we used that for any $\rho\succeq 0$ and $\Pi$ orthogonal projector, we have

F(\Pi\rho\Pi,\sigma)=\mbox{Tr}\left[\sqrt{\sqrt{\Pi\rho\Pi}\sigma\sqrt{\Pi\rho\Pi}}\right]=\mbox{Tr}\left[\sqrt{\sqrt{\Pi\rho\Pi}\Pi\sigma\Pi\sqrt{\Pi\rho\Pi}}\right]=F(\Pi\rho\Pi,\Pi\sigma\Pi).

3 Fidelity estimation with block-encoding based algorithms

Let us now sketch the main idea behind our block-encoding algorithm that builds on the Hadamard test. Suppose that $\sqrt{\rho}\sqrt{\sigma}$ has singular value decomposition $U\Sigma V^{\dagger}$ , then

\displaystyle F(\rho,\sigma)=\left\lVert\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}=\mbox{Tr}\left[U^{\dagger}\sqrt{\rho}\sqrt{\sigma}V\right]=\mbox{Tr}\left[\sqrt{\rho}\sqrt{\sigma}VU^{\dagger}\right]=\mbox{Tr}\left[\rho(\rho^{-\frac{1}{2}}\sqrt{\sigma}VU^{\dagger})\right].

(20)

By Lemma 9 it suffices to implement a (subnormalized) block-encoding of $(\rho^{-\frac{1}{2}}\sqrt{\sigma}VU^{\dagger})$ in order to use the Hadamard test for computing an estimate of the fidelity. The main issue with this approach is that $\rho^{-\frac{1}{2}}$ can, in general, have arbitrarily large eigenvalues.

In order to deal with singularities arising from small eigenvalues we modify the task as follows. Let $\rho_{I}:=\Pi^{\phantom{\rho}}_{I}\rho\Pi^{\phantom{\rho}}_{I}$ denote the subnormalized density matrix we get after throwing away eigenvalues outside the interval $I$ . For some $\delta,\varepsilon\in[0,\frac{1}{2}]$ we wish to provide an estimate $f$ of $F(\rho_{[\theta,1]},\sigma)$ to precision

\sqrt{\mbox{Tr}\left[\rho_{[(1-\delta)\theta,\theta)}\right]}\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[(1-\delta)\theta,\theta)}\sigma\Pi^{\phantom{\rho}}_{[(1-\delta)\theta,\theta)}\right]}+\varepsilon,

(21)

in turn providing an estimate of $F(\rho,\sigma)$ with precision

\sqrt{\mbox{Tr}\left[\rho_{[0,\theta)}\right]}\sqrt{\mbox{Tr}\left[\Pi^{\phantom{\rho}}_{[0,\theta)}\sigma\Pi^{\phantom{\rho}}_{[0,\theta)}\right]}+\varepsilon.

(22)

For this we shall use some soft threshold function $t\colon\mathbb{R}\mapsto[0,1]$ such that for all $x<(1-\delta)\theta$ we have $t(x)=0$ and for all $x\geq\theta$ we have $t(x)=1$ . By Equation 19 and Corollary 24 we have that $f:=F(t^{2}(\rho)\rho,\sigma)$ satisfies both the above requirements with $\varepsilon=0$ , so it suffices to compute $F(t^{2}(\rho)\rho,\sigma)$ with $\varepsilon$ -precision.

In the following we analyze the propagation of errors and the complexity of the implementation.

3.1 Polynomial approximations and error bounds

In order to make the procedure more efficient we can approximate $\sqrt{\sigma}$ . Let $\frac{1}{2}s\colon\mathbb{R}\mapsto[0,1]$ be a polynomial function provided by Corollary 18 with parameters $\delta^{\prime}\leftarrow\frac{\varepsilon^{2}}{160\mathrm{rk}\left(t(\rho)\right)}$ and $\varepsilon^{\prime}\leftarrow\frac{\varepsilon}{20\sqrt{\mathrm{rk}\left(t(\rho)\right)}}$ , then $\left\lVert s\right\rVert_{[-1,1]}\leq 1+\varepsilon^{\prime}\leq 1+\varepsilon\leq 2$ and

$\displaystyle\left\|\left\lVert t(\rho)\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}-\left\lVert t(\rho)\sqrt{\rho}s(\sigma)\right\rVert_{1}\right\|$	$\displaystyle\leq\left\lVert t(\rho)\sqrt{\rho}(\sqrt{\sigma}-s(\sigma))\right\rVert_{1}$	(by the triangle inequality)
	$\displaystyle\leq\left\lVert t(\rho)\right\rVert_{2}\left\lVert\sqrt{\rho}\right\rVert_{2}\left\lVert\sqrt{\sigma}-s(\sigma)\right\rVert$	(by Hölder’s inequality)
	$\displaystyle\leq\sqrt{\mathrm{rk}\left(t(\rho)\right)}\cdot 1\cdot\frac{\varepsilon}{5\sqrt{\mathrm{rk}\left(t(\rho)\right)}}=\frac{\varepsilon}{5}.$	(23)

Let us also approximate $\sqrt{\rho}$ by a polynomial. We construct a polynomial $\frac{\sqrt{\theta}}{2\sqrt{2}}q$ by Lemma 16, with parameters $\delta^{\prime}\leftarrow\frac{\theta}{2}$ and $\varepsilon^{\prime}\leftarrow\frac{\varepsilon\sqrt{\theta}}{20\sqrt{2}}$ , then $\left\lVert xq(x)\right\rVert_{[-1,1]}\leq 2$ and

$\displaystyle\left\|\left\lVert t(\rho)\sqrt{\rho}s(\sigma)\right\rVert_{1}-\left\lVert t(\rho)\rho q(\rho)s(\sigma)\right\rVert_{1}\right\|$	$\displaystyle\leq\left\lVert t(\rho)\rho\left(\frac{1}{\sqrt{\rho}}-q(\rho)\right)s(\sigma)\right\rVert_{1}$	(by the triangle inequality)
	$\displaystyle\leq\left\lVert\rho\right\rVert_{1}\left\lVert\frac{t(\rho)}{\sqrt{\rho}}-t(\rho)q(\rho)\right\rVert\left\lVert s(\sigma)\right\rVert$	(by Hölder’s inequality)
	$\displaystyle\leq 1\cdot\frac{\varepsilon}{10}\cdot 2=\frac{\varepsilon}{5}.$	(24)

Let $Q:=t(\rho)\rho q(\rho)s(\sigma)$ , by definition we have that

\left\lVert Q\right\rVert_{1}=\mbox{Tr}\left[Q\cdot SV^{(\mathrm{sgn})}\!(Q^{\dagger})\right].

(25)

Observe that $Q$ has at most $\mathrm{rk}\left(t(\rho)\right)$ non-zero singular values. Also by Hölder’s inequality we know that $\left\lVert t(\rho)\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}\leq\left\lVert t(\rho)\right\rVert\left\lVert\sqrt{\rho}\right\rVert_{2}\left\lVert\sqrt{\sigma}\right\rVert_{2}\leq 1$ , so by Equations 23 and 24 we get $\left\lVert Q\right\rVert_{1}\leq 1+\frac{2\varepsilon}{5}\leq\frac{6}{5}$ . Let $\varsigma_{i}$ denote the singular values of $Q$ . Let $p\colon[-1,1]\mapsto[-1,1]$ be an odd function such that for all $x\in[-1,1]$ with $|x|\geq\frac{\varepsilon}{160\mathrm{rk}\left(t(\rho)\right)}$ we have $\mathrm{sgn}(x)-p(x)\leq 5\varepsilon/60$ , then

$\displaystyle\left\|\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(\mathrm{sgn})}\!(Q^{\dagger})\right]-\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)\right]\right\|$	$\displaystyle=\left\|\sum_{i=1}\varsigma_{i}-\varsigma_{i}p(\varsigma_{i}/8)\right\|$
	$\displaystyle\leq\sum_{i=1}\varsigma_{i}\|1-p(\varsigma_{i}/8)\|$
	$\displaystyle=\kern-2.84526pt\sum_{\varsigma_{i}\geq\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\kern-8.53581pt\varsigma_{i}\|1-p(\varsigma_{i}/8)\|+\kern-2.84526pt\sum_{\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\kern-8.53581pt\varsigma_{i}\|1-p(\varsigma_{i}/8)\|$
	$\displaystyle\leq\sum_{\varsigma_{i}\geq\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\varsigma_{i}\frac{5\varepsilon}{60}+\sum_{\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}2\varsigma_{i}$
	$\displaystyle\leq\mbox{Tr}\left[\Sigma\right]\frac{\varepsilon}{10}+\sum_{0<\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}2\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}\leq\frac{\varepsilon}{5}.$	(26)

So performing singular value transformation according to $p$ as opposed to $\mathrm{sgn}$ introduces an error in the estimation that is bounded by $\varepsilon/5$ . Moreover, due to Lemma 14 we can find such a polynomial $p$ with degree $\mathcal{O}\left(\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon}\log(1/\varepsilon)\right)$ .

Now we find a polynomial approximation of $t(x)$ . We choose a polynomial $\tilde{t}(x):=1-P^{\prime}(x)$ for a rectangle function from Corollary 15 with parameters $t^{\prime}=(1-\frac{\delta}{2})\theta$ , $\delta^{\prime}\leftarrow\frac{\delta}{2}\theta$ , and $\varepsilon^{\prime}\leftarrow\min\{\frac{\sqrt{\theta}\varepsilon}{20\sqrt{2}},\frac{5\varepsilon}{60\deg(p)}\}$ . The degree of $\tilde{t}$ is $\mathcal{O}\left(\frac{1}{\delta\theta}\log\frac{\mathrm{rk}\left(t(\rho)\right)}{\theta\varepsilon}\right)$ . Note that we did not yet even specify the shape of $t(x)$ for $x\in[(1-\frac{\delta}{2})\theta,\theta]$ , so we simply define it to be $\tilde{t}(x)$ – resulting in $|t(x)-\tilde{t}(x)|\leq\varepsilon^{\prime}$ for $x\in[0,1]$ . Let us define $\tilde{Q}:=\tilde{t}(\rho)\rho q(\rho)s(\sigma)$ , then we wish to bound

$\displaystyle\left\|\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)\right]-\mbox{Tr}\left[\tilde{Q}\cdot\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right]\right\|$	$\displaystyle\leq\left\|\mbox{Tr}\left[Q\cdot\left(\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)-\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right)\right]\right\|$
	$\displaystyle\phantom{=}+\left\|\mbox{Tr}\left[\left(Q-\tilde{Q}\right)\cdot\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right]\right\|$	(triangle inequality)
	$\displaystyle\leq\left\lVert Q\cdot\left(\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)-\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right)\right\rVert_{1}$
	$\displaystyle\phantom{=}+\left\lVert\left(Q-\tilde{Q}\right)\cdot\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right\rVert_{1}$	(trace-norm inequality)
	$\displaystyle\leq\underset{\leq\frac{6}{5}}{\underbrace{\left\lVert Q\right\rVert_{1}}}\cdot\left\lVert\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)-\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right\rVert$	(Hölder’s ineq.)
	$\displaystyle\phantom{=}+\left\lVert Q-\tilde{Q}\right\rVert_{1}\cdot\underset{\leq 1}{\underbrace{\left\lVert\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right\rVert}}.$

Now we bound both terms individually as follows

$\displaystyle\left\lVert Q-\tilde{Q}\right\rVert_{1}$	$\displaystyle=\left\lVert t(\rho)\rho q(\rho)s(\sigma)\right\rVert_{1}-\left\lVert\tilde{t}(\rho)\rho q(\rho)s(\sigma)\right\rVert_{1}$	(by definition)
	$\displaystyle\leq\left\lVert(t(\rho)-\tilde{t}(\rho))\rho q(\rho)s(\sigma)\right\rVert_{1}$	(by the triangle inequality)
	$\displaystyle\leq\left\lVert(t(\rho)-\tilde{t}(\rho))\right\rVert\left\lVert\rho\right\rVert_{1}\left\lVert q(\rho)\right\rVert\left\lVert s(\sigma)\right\rVert$	(by Hölder’s inequality)
	$\displaystyle\leq\varepsilon^{\prime}\cdot 1\cdot\frac{2\sqrt{2}}{\sqrt{\theta}}\cdot 2=\frac{\varepsilon}{10},$

and by Lemma 11 (choosing $\bar{A}:=\bar{Q}/8$ , and observing $\left\lVert\bar{Q}\right\rVert\leq 4$ for $\bar{Q}\in\{Q,\tilde{Q}\}$ ) we have that

	$\displaystyle\left\lVert\mathrm{SV}^{(p)}\!(Q/8)-\mathrm{SV}^{(p)}\!(\tilde{Q}/8)\right\rVert$	$\displaystyle\leq\deg(p)\cdot\frac{1}{4}\cdot\left\lVert\tilde{Q}-Q\right\rVert$
		$\displaystyle=\frac{\deg(p)}{4}\left\lVert(t(\rho)-\tilde{t}(\rho))\rho q(\rho)s(\sigma)\right\rVert$
		$\displaystyle\leq\frac{\deg(p)}{4}\left\lVert t(\rho)-\tilde{t}(\rho)\right\rVert\left\lVert\rho q(\rho)\right\rVert\left\lVert s(\sigma)\right\rVert$
		$\displaystyle\leq\frac{\deg(p)}{4}\cdot\varepsilon^{\prime}\cdot 4\leq\frac{5\varepsilon}{60},$

ultimately resulting in

\left|\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)\right]-\mbox{Tr}\left[\tilde{Q}\cdot\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right]\right|\leq\frac{\varepsilon}{5}.

(27)

Therefore, combining Equations 23, 24, 25, 26 and 27 we ultimately get

\left|\left\lVert t(\rho)\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}-\mbox{Tr}\left[\tilde{Q}\cdot\mathrm{SV}^{(p)}\!(\tilde{Q}^{\dagger}/8)\right]\right|\leq\frac{4\varepsilon}{5}.

(28)

3.2 Complexity analysis – with purified access

Assuming the purified access model, we can show the following result for our block-encoding algorithm.

Theorem 25.

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be arbitrary density matrices. Suppose also that $\rho$ has the smallest rank of the two states. Let $\varepsilon\in(0,1)$ , $\delta\in(0,1)$ and $\theta\in(0,1)$ be parameters. Then, given purified access to $\rho$ and $\sigma$ , our block-encoding algorithm in Section 3 runs in time

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(\Pi^{\rho}_{[(1-\delta)\theta,1]}\right)}{\varepsilon^{2}\delta\theta^{\frac{3}{2}}}T_{\rho}+\frac{\mathrm{rk}\left(\Pi^{\rho}_{[(1-\delta)\theta,1]}\right)^{2}}{\varepsilon^{4}\theta^{\frac{1}{2}}}T_{\sigma}\right)

and outputs (with high probability) an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that

|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon,

where $\rho_{\theta}$ is a “soft-thresholded” version of $\rho$ in which eigenvalues of $\rho$ below $(1-\delta)\theta$ are removed and those above $\theta$ are kept intact, while eigenvalues in $[(1-\delta)\theta,\theta]$ are decreased by some amount.

Proof.

By Lemma 9 it suffices to implement a ( $\frac{\sqrt{\theta}}{4\sqrt{2}}$ subnormalized) block-encoding of

\tilde{t}(\rho)q(\rho)s(\sigma)\mathrm{SV}^{(p)}\!(\tilde{t}(\rho)\rho q(\rho)s(\sigma)/8)

in order to use the Hadamard test for computing an estimate of $\mbox{Tr}\left[\rho\tilde{t}(\rho)q(\rho)s(\sigma)\mathrm{SV}^{(p)}\!(\tilde{t}(\rho)\rho q(\rho)s(\sigma)/8)\right]$ – which by Equation 28 is $\frac{4\varepsilon}{5}$ -close to the fidelity. Using amplitude estimation we can estimate the success probability of the Hadamard test to precision $\mathcal{O}\left(\frac{\varepsilon\sqrt{\theta}}{40\sqrt{2}}\right)$ , which then results in an $\varepsilon$ precise estimate of

F(t^{2}(\rho)\rho,\sigma)=\left\lVert t(\rho)\sqrt{\rho}\sqrt{\sigma}\right\rVert_{1}.

(29)

In order to implement a block-encoding of $\tilde{Q}/8$ we implement both $\tilde{t}(\rho)\rho q(\rho)/4$ and $s(\sigma)/2$ using QSVT and take their product [GSLW19]. Sine $\deg(s)=\mathcal{O}\left(\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon^{2}}\log\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon}\right)$ and $\deg(\tilde{t}(x)\cdot x\cdot q(x))$ is $\mathcal{O}\left(\frac{1}{\delta\theta}\log\frac{\mathrm{rk}\left(t(\rho)\right)}{\theta\varepsilon}\right)$ , the complexity of implementing $\frac{\tilde{t}(\rho)\rho q(\rho)s(\sigma)}{8}$ is $\mathcal{O}\left(\frac{T_{\rho}}{\delta\theta}\log\frac{\mathrm{rk}\left(t(\rho)\right)}{\theta\varepsilon}+\frac{T_{\sigma}\mathrm{rk}\left(t(\rho)\right)}{\varepsilon^{2}}\log\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon}\right)$ . Applying QSVT by the polynomial $p$ uses the above block-encoding a total of $\mathcal{O}\left(\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon}\log(1/\varepsilon)\right)$ times resulting in complexity $\widetilde{\mathcal{O}}\left(\frac{T_{\rho}\mathrm{rk}\left(t(\rho)\right)}{\varepsilon\delta\theta}+\frac{T_{\sigma}\mathrm{rk}\left(t(\rho)\right)^{2}}{\varepsilon^{3}}\right)$ for implementing $\mathrm{SV}^{(p)}\!(\tilde{t}(\rho)\rho q(\rho)s(\sigma)/8)$ . Similarly we can implement a block-encoding pf $\tilde{t}(\rho)q(\rho)\cdot\frac{\sqrt{\theta}}{2\sqrt{2}}$ by QSVT and take its product with $s(\sigma)/2$ , then take a product with $\mathrm{SV}^{(p)}\!(\tilde{t}(\rho)\rho q(\rho)s(\sigma)/8)$ , which only gives a lower order contribution to the running time. Before providing the final runtime bound let us note that $\mathrm{rk}\left(t(\rho)\right)\leq\mathrm{rk}\left(\Pi^{\rho}_{[(1-\delta)\theta,1]}\right)$ . This all together gives the runtime bound

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(\Pi^{\rho}_{[(1-\delta)\theta,1]}\right)}{\varepsilon^{2}\delta\theta^{\frac{3}{2}}}T_{\rho}+\frac{\mathrm{rk}\left(\Pi^{\rho}_{[(1-\delta)\theta,1]}\right)^{2}}{\varepsilon^{4}\theta^{\frac{1}{2}}}T_{\sigma}\right).

This proves the claim. ∎

If we have an upper bound on the smallest rank of the two states, we obtain the following result.

Corollary 26.

\widetilde{\mathcal{O}}\left(\frac{r^{\frac{5}{2}}}{\varepsilon^{5}}(T_{\rho}+T_{\sigma})\right).

and outputs (with high probability) an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that

|F(\rho,\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon,

where $\rho_{\theta}$ is a “soft-thresholded” version of $\rho$ with parameters $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ and $\delta=\frac{1}{2}$ .

Proof.

Assuming that $\mathrm{rk}\left(\rho\right)\leq r$ , we can use Equation 22 and set $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ and $\delta=\frac{1}{2}$ to obtain an $\varepsilon$ -precise estimate $\hat{F}(\rho_{\theta},\sigma)$ to the fidelity $F(\rho,\sigma)$ in complexity $\widetilde{\mathcal{O}}\left(\frac{r^{\frac{5}{2}}}{\varepsilon^{5}}(T_{\rho}+T_{\sigma})\right)$ . ∎

3.3 Complexity analysis – with sampling access

Assuming the purified access model, we can show the following result for our block-encoding algorithm.

Theorem 27.

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be arbitrary density matrices. Suppose also that $\rho$ has the smallest rank of the two states. Let $\varepsilon\in(0,1)$ , $\delta\in(0,1)$ and $\theta\in(0,1)$ be parameters. Then, given sampling access to $\rho$ and $\sigma$ , our block-encoding algorithm in Section 3 uses

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)^{2}}{\varepsilon^{5}\delta^{2}\theta^{\frac{7}{2}}}\right)\text{ copies of }\rho\text{ and}\,\,\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)^{4}}{\varepsilon^{9}\theta^{\frac{3}{2}}}\right)\text{ copies of }\sigma,

where $t$ is the threshold function in Section 3, and outputs (with high probability) $\hat{F}(\rho_{\theta},\sigma)$ such that

|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon,

Proof.

The overall approach is analogous to Section 3.2. We implement a $\Theta(\sqrt{\theta})$ -subnormalized ( $\mathcal{O}(\varepsilon\sqrt{\theta})$ -approximate) block-encoding of

\tilde{t}(\rho)q(\rho)s(\sigma)\mathrm{SV}^{(p)}\!(\tilde{t}(\rho)\rho q(\rho)s(\sigma)/8),

(30)

and apply the block-Hadamard test on a copy of $\rho$ as described in Lemma 9. Then it suffices to estimate the probability of outcome $0$ to precision $\mathcal{O}\left(\varepsilon\sqrt{\theta}\right)$ . Since in this scenario we only get copies of $\rho$ we cannot implement amplitude estimation (which would require the ability to prepare $\rho$ ), but need to repeat the Hadamard test a total of $\frac{1}{\varepsilon^{2}\theta}$ times to get an estimate with such precision.

Another difference from Section 3.2 is that we do not natively have a perfect block-encoding of $\rho$ and $\sigma$ . However, using density matrix exponentiation by Corollary 21 we can get an approximate block-encoding of a $\frac{\pi}{4}$ -subnormalized block-encoding of $\rho$ and $\sigma$ . The $\frac{\pi}{4}$ -subnormalization constant only induces constant factor changes in our analysis in Section 3.1. In particular our earlier analysis in Section 3.2 still shows that a $\Theta(\sqrt{\theta})$ -subnormalized block-encoding of (30) can be implemented with $b_{\rho}:=\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)}{\varepsilon\delta\theta}\right)$ uses of $U_{\rho}$ a block-encoding of $\frac{\pi}{4}\rho$ and $b_{\sigma}:=\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)^{2}}{\varepsilon^{3}}\right)$ uses of $U_{\sigma}$ a block-encoding of $\frac{\pi}{4}\sigma$ . We can implement $U_{\rho}$ to $\mathcal{O}\left(\frac{\varepsilon\sqrt{\theta}}{b_{\rho}}\right)$ -error in the diamond norm by using $\widetilde{\mathcal{O}}\left(\frac{b_{\rho}}{\varepsilon\sqrt{\theta}}\right)$ copies of $\rho$ and similarly implement $U_{\sigma}$ to $\mathcal{O}\left(\frac{\varepsilon\sqrt{\theta}}{b_{\sigma}}\right)$ -error in the diamond norm by using $\widetilde{\mathcal{O}}\left(\frac{b_{\sigma}}{\varepsilon\sqrt{\theta}}\right)$ copies of $\sigma$ due to Corollary 21. This then results in an implementation of a block-encoding of a $\Theta(\sqrt{\theta})$ -subnormalized version of (30) up to diamond-norm error $\mathcal{O}\left(\varepsilon\sqrt{\theta}\right)$ using $\widetilde{\mathcal{O}}\left(\frac{b_{\rho}^{2}}{\varepsilon\sqrt{\theta}}\right)$ copies of $\rho$ and $\widetilde{\mathcal{O}}\left(\frac{b_{\sigma}^{2}}{\varepsilon\sqrt{\theta}}\right)$ copies of $\sigma$ . This construction ensures that the probability of getting outcome $0$ in the Hadamard-test is $\mathcal{O}\left(\varepsilon\sqrt{\theta}\right)$ -close to the outcome one would get by using an exact block-encoding of (30). By appropriately choosing the constants this ensures that repeating the Hadamard-test $\frac{1}{\varepsilon^{2}\theta}$ times with high-probability we get an $\varepsilon$ -precisie estimate of (29).

This amounts to an algorithm that uses $\widetilde{\mathcal{O}}\left(\frac{b_{\rho}^{2}}{\varepsilon^{3}\theta^{\frac{3}{2}}}\right)$ and $\widetilde{\mathcal{O}}\left(\frac{b_{\sigma}^{2}}{\varepsilon^{3}\theta^{\frac{3}{2}}}\right)$ copies of $\rho$ and $\sigma$ respectively, i.e., the ultimate algorithm uses

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)^{2}}{\varepsilon^{5}\delta^{2}\theta^{\frac{7}{2}}}\right)\text{ copies of }\rho

and

\widetilde{\mathcal{O}}\left(\frac{\mathrm{rk}\left(t(\rho)\right)^{4}}{\varepsilon^{9}\theta^{\frac{3}{2}}}\right)\text{ copies of }\sigma.

The implementation is also gate efficient, the gate complexity overhead of the implementation is $\mathcal{O}\left(\log(\dim(\rho))\right)$ as follows from Corollary 21. This proves the claim. ∎

If we have an upper bound on the smallest rank of the two states, we get the following:

Corollary 28.

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be arbitrary density matrices. Suppose also that $\rho$ has the smallest rank of the two states, where $\mathrm{rk}\left(\rho\right)\leq r$ . Let $\varepsilon\in(0,1)$ be a parameter. Then, given sampling access to $\rho$ and $\sigma$ , our block-encoding algorithm in Section 3 uses $\widetilde{\mathcal{O}}\left(\frac{r^{5.5}}{\varepsilon^{12}}\right)$ copies of $\rho$ and $\sigma$ and outputs (with high probability) an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that

|F(\rho,\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon,

where $\rho_{\theta}$ is a “soft-thresholded” version of $\rho$ with parameters $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ and $\delta=\frac{1}{2}$ .

Proof.

In case we know that $\mathrm{rk}\left(\rho\right)\leq r$ , by Equation 22 we can see that by setting $\theta=\Theta(\frac{\varepsilon^{2}}{r})$ and $\delta=\frac{1}{2}$ we can obtain an $\varepsilon$ -precise estimate with $\widetilde{\mathcal{O}}\left(\frac{r^{5.5}}{\varepsilon^{12}}\right)$ copies of $\rho$ and $\sigma$ . ∎

4 Fidelity estimation via spectral sampling

In this section, we present our second approximation algorithm for the problem of fidelity estimation. Recall that we can write the fidelity between two states $\rho,\sigma\in\mathbb{C}^{d\times d}$ as the quantity $F(\rho,\sigma)=\mbox{Tr}{[\sqrt{\Lambda}]}$ , where $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}\in\mathbb{C}^{d\times d}$ has the following non-trivial entries in the eigenbasis of $\rho$ :

\displaystyle\Lambda_{ij}=\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle,\quad\quad\forall i,j\in\{1,\dots,\mathrm{rk}\left(\rho\right)\}.

Hence, it suffices to directly compute the matrix elements of $\Lambda$ in order to estimate the fidelity $F(\rho,\sigma)$ . In Section 4.1, we show a continuity bound for fidelity estimation. This allows us to quantify the approximation error of our estimate $\hat{F}(\rho,\sigma)=\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}$ for $F(\rho,\sigma)=\mbox{Tr}{[\sqrt{\Lambda}]}$ in terms of the approximation error between $\Lambda$ and $\hat{\Lambda}$ , where $\hat{\Lambda}$ is our initial estimate and $\hat{\Lambda}_{+}$ is the projection onto the positive semidefinite cone. In particular, we show in Theorem 29 that

|F(\rho,\sigma)-\hat{F}(\rho,\sigma)|\,\leq\,\sqrt{2}r\cdot\sqrt{\left\|\Lambda-\hat{\Lambda}\right\|_{1}},

where $r=\min\{\mathrm{rk}\left(\rho\right),\mathrm{rk}\left(\sigma\right)\}$ denotes the smallest rank of the states $\rho$ and $\sigma$ . To estimate the eigenvalues of $\rho$ , we rely on the technique called quantum spectral sampling first introduced by Lloyd, Mohseni and Rebentrost [LMR14] in the context of quantum principal component analysis.

4.1 Continuity bound for fidelity estimation

In this section, we prove an important technical result which allows us to relate the approximation error of a fidelity estimate for $F(\rho,\sigma)=\mbox{Tr}{[\sqrt{\Lambda}]}$ in terms of the approximation error for the matrix $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}$ . In other words, if $\Lambda$ and $\hat{\Lambda}$ are close, then the fidelity estimate $\hat{F}(\rho,\sigma)=\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}$ is close to $F(\rho,\sigma)=\mbox{Tr}{[\sqrt{\Lambda}]}$ , where $\hat{\Lambda}_{+}$ is the projection onto the positive semidefinite cone.

Theorem 29 (Continuity bound).

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be density matrices and $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}$ . Let $\hat{\Lambda}\in\mathbb{C}^{d\times d}$ be an arbitrary hermitian matrix with $\mathrm{rk}(\hat{\Lambda})\leq\mathrm{rk}\left(\Lambda\right)$ and suppose that $\hat{F}(\rho,\sigma)=\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}$ , where $\hat{\Lambda}_{+}$ is the projection of $\hat{\Lambda}$ onto the positive semidefinite cone. Then,

|F(\rho,\sigma)-\hat{F}(\rho,\sigma)|=\big{|}\mathrm{Tr}\big{[}{\sqrt{\Lambda}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}\big{|}\,\leq\,\sqrt{2}r\cdot\sqrt{\left\|\Lambda-\hat{\Lambda}\right\|_{1}},

where we let $r=\min\{\mathrm{rk}\left(\rho\right),\mathrm{rk}\left(\sigma\right)\}$ denote the smallest rank of the states $\rho$ and $\sigma$ .

Proof.

Note that $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}$ is positive semidefinite and hermitian with $\mathrm{rk}\left(\Lambda\right)\leq r$ . Let $\operatorname{spec}(\Lambda)=(\lambda_{1},\dots,\lambda_{r})$ and $\operatorname{spec}(\hat{\Lambda}_{+})=(\hat{\lambda}^{+}_{1},\dots,\hat{\lambda}^{+}_{r})$ be the spectra of $\Lambda$ and $\hat{\Lambda}_{+}$ , respectively. Then,

$\displaystyle\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}^{+}}}\big{]}\big{\|}$	$\displaystyle=\big{\|}\sum_{i=1}^{r}\sqrt{\lambda_{i}}-\sum_{i=1}^{r}\sqrt{\hat{\lambda}^{+}_{i}}\big{\|}$	$\displaystyle(\text{by definition})$
	$\displaystyle\leq\sum_{i=1}^{r}\big{\|}\sqrt{\lambda_{i}}-\sqrt{\hat{\lambda}_{i}^{+}}\big{\|}$	$\displaystyle(\text{triangle inequality})$
	$\displaystyle\leq r\cdot\max_{1\leq i\leq r}\big{\|}\sqrt{\lambda_{i}}-\sqrt{\hat{\lambda}_{i}^{+}}\big{\|}$
	$\displaystyle\leq r\cdot\left\\|\sqrt{\Lambda}-\sqrt{\hat{\Lambda}_{+}}\right\\|_{2}$	$\displaystyle(\text{\cite[cite]{[\@@bibref{}{bhatia1997MatrixAnalysis}{}{}, Problem III.6.13]}})$
	$\displaystyle\leq r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}_{+}\right\\|_{1}}$	$\displaystyle(\text{\cite[cite]{[\@@bibref{}{cmp/1103842028}{}{}, Lemma 4.1]}})$

Note that $\hat{\Lambda}=\hat{\Lambda}_{+}+\hat{\Lambda}_{-}$ . Thus, by the triangle inequality, we obtain

\displaystyle\left\|\Lambda-\hat{\Lambda}_{+}\right\|_{1}=\left\|\Lambda-\hat{\Lambda}+\hat{\Lambda}_{-}\right\|_{1}\leq\left\|\Lambda-\hat{\Lambda}\right\|_{1}+\left\|\hat{\Lambda}_{-}\right\|_{1}.

(31)

Let us now bound the quantity $\|\hat{\Lambda}_{-}\|_{1}$ . Recall from Lemma 5 that the projection $\hat{\Lambda}_{+}$ satisfies

\displaystyle\hat{\Lambda}_{+}=\underset{X\succeq 0}{\mathrm{argmin}}\left\|X-\hat{\Lambda}\right\|_{1}.

(32)

In other words, Eq. (32) yields the following inequality

\displaystyle\left\|X-\hat{\Lambda}\right\|_{1}\geq\left\|\hat{\Lambda}_{+}-\hat{\Lambda}\right\|_{1}=\left\|\hat{\Lambda}_{-}\right\|_{1},\quad\quad\forall X\succeq 0.

(33)

Using the fact that $\Lambda\succeq 0$ , Eq. (33) implies the following upper bound given by

\displaystyle\left\|\hat{\Lambda}_{-}\right\|_{1}\leq\left\|\Lambda-\hat{\Lambda}\right\|_{1}.

(34)

Putting everything together, we can use (31) and (34) to conclude that

	$\displaystyle\|F(\rho,\sigma)-\hat{F}(\rho,\sigma)\|$	$\displaystyle=\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}\big{\|}$
		$\displaystyle\leq r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}_{+}\right\\|_{1}}$
		$\displaystyle=\sqrt{2}r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}\right\\|_{1}}.\quad\quad\quad$

This proves the claim. ∎

4.2 Algorithm

Let us now state our second approximation algorithm for fidelity estimation via quantum spectral sampling. Recall that, in Section 2.2 and Section 2.3, we reviewed the Swap Test and Hadamard Test, respectively. We will review the quantum spectral sampling algorithm in Section 4.3 and the quantum eigenstate filtering algorithm in Section 4.4.

Finally, we also remark that the analysis of the algorithm is given in Section 4.7.

Input: Purified access to density operators

\rho,\sigma\in\mathbb{C}^{d\times d}

via unitaries

U_{\rho}

and

U_{\sigma}

Promise:

\rho=\sum_{i=1}^{\mathrm{rk}\left(\rho\right)}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert

has well-separated spectrum with gap

\Delta>0

Parameters:

\varepsilon\in(0,1)

\delta\in(0,1)

\theta\in(0,1)

1 Set

m\leftarrow\lceil\frac{1}{2\theta}\rceil

2 repeat $\left\lceil\frac{16\cdot mH_{m}}{\delta\cdot\theta^{2}}\right\rceil$ times

3 Run the subroutine Quantum_Spectral_Sampling

(U_{\rho},\gamma,\ell)

with

\gamma=\min\{\frac{\theta^{3}\varepsilon}{\sqrt{2}},\frac{\Delta}{2}\}

and

\ell=\left\lceil\log\big{(}\frac{16}{\delta\cdot\theta^{2}}\big{)}+\log\lceil\frac{16mH_{m}}{\delta\cdot\theta^{2}}\rceil\right\rceil

to generate the bipartite state

\sum_{i=1}^{\mathrm{rk}\left(\rho\right)}\lambda_{i}\,\lvert\psi_{i}\rangle\langle\psi_{i}\rvert\otimes\lvert\tilde{\lambda}_{i}\rangle\langle\tilde{\lambda}_{i}\rvert

4 Measure the second register in the computational basis to obtain a sample

\tilde{\lambda}_{j}

5 Discard

\tilde{\lambda}_{j}

if it is smaller than

\frac{3\theta}{4}

or within distance

\frac{\Delta}{2}+\gamma

of any previously seen sample.

6 end

7Let

(\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{r_{\theta}})

be the collected eigenvalues in decreasing order.

8 for $i=1$ to $r_{\theta}$ do

9 for $j=i$ to $r_{\theta}$ do

10 if $i=j$ then

11 Run Swap_Test

(U_{i}\lvert 0^{\log d}\rangle,\sigma,\xi,\nu)

with

\xi=\frac{\varepsilon^{2}}{8r_{\theta}^{4}}

and

\nu=\frac{\delta}{2r_{\theta}^{2}}

to obtain

\tilde{\sigma}_{ii}\approx\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle

, for the unitary

U_{i}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{i},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})

in Theorem 32

12 Set

\hat{\Lambda}_{ii}^{\theta}\leftarrow\tilde{\lambda}_{i}\tilde{\sigma}_{ii}

14 else

15 Run Hadamard_Test

(\lvert 0^{\log d}\rangle,U_{i}^{\dagger}UU_{j},\xi,\nu)

with parameters

\xi=\frac{\varepsilon^{2}}{8r_{\theta}^{4}}

and

\nu=\frac{\delta}{2r_{\theta}^{2}}

, for the unitaries

U_{j}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{j},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})

and

U_{i}^{\dagger}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{i},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})^{\dagger}

in Theorem 32, and where

U

is a block-encoding of

\sigma

as in Lemma 7, to obtain

\tilde{\sigma}_{ij}\approx\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle

16 Set

\hat{\Lambda}^{\theta}_{ij}\leftarrow\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}\,\tilde{\sigma}_{ij}

17 Set

\hat{\Lambda}^{\theta}_{ji}=\hat{\Lambda}^{\theta^{*}}_{ij}

18 end if

20 end for

22 end for

Output the estimate

\hat{F}(\rho_{\theta},\sigma)=\mathrm{Tr}\left[\sqrt{\hat{\Lambda}_{+}^{\theta}}\right]

Algorithm 1 Fidelity estimation via spectral sampling

4.3 Quantum spectral sampling

The following definition of the infinitesimal SWAP operation allows us to approximately implement the density matrix exponential $U=e^{-2\pi i\rho}$ , as shown by Lloyd, Mohseni and Rebentrost [LMR14] in the context of quantum principal component analysis, and later extended by Prakash [Pra14].

Definition 30 (Infinitesimal SWAP operation).

Let $k\in\mathbb{N}$ be a parameter and let $\rho$ be a mixed state. We define the action $U^{(k)}\sigma U^{(k)}{}^{\dagger}:=\sigma^{(k)}$ on an input state $\sigma$ implicitly via the following iteration:

	$\displaystyle\sigma^{(0)}=\sigma$
	$\displaystyle\sigma^{(n+1)}=\mathrm{Tr}_{B}\big{[}e^{-2\pi iS/k}\big{(}\sigma^{(n)}_{A}\otimes\rho_{B}\big{)}e^{2\pi iS/k}\big{]},\quad 0\leq n\leq k-1,$

where $S$ is the SWAP operator.

Prakash [Pra14] showed that the procedure $U^{(k)}\approx e^{-2\pi i\rho}$ can be implemented in time $O(kT_{\rho})$ , where $T_{\rho}$ is the time it takes to prepare the state $\rho$ given purified access $U_{\rho}$ . In particular, the following theorem due Prakash [Pra14, Theorem 3.2.1] states that it is possible to sample from the eigenvalues of a density operator $\rho\in\mathbb{C}^{d\times d}$ in time $\tilde{O}(T_{\rho}/\gamma^{3})$ with high probability.

Theorem 31 ([Pra14]).

Let $\rho\in\mathbb{C}^{d\times d}$ be a density matrix, and let $T_{\rho}$ denote the time it takes to prepare $\rho$ via purified access to the unitary $U_{\rho}$ . Let $\gamma\in(0,1)$ and $\ell\in\mathbb{N}$ be parameters. Then, the procedure Quantum_Spectral_Sampling $(U_{\rho},\gamma,\ell)$ in Algorithm 2 runs in time $\tilde{O}(\ell\cdot T_{\rho}/\gamma^{3})$ and produces eigenvalue estimates $\tilde{\lambda}_{j}$ such that $|\tilde{\lambda}_{j}-\lambda_{j}|\leq\gamma$ with probability at least $1-2^{-\ell}$ .

Input: Purified access to

\rho\in\mathbb{C}^{d\times d}

with spectral decomposition

\rho=\sum_{i=1}^{\mathrm{rk}\left(\rho\right)}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert

via

U_{\rho}

. Parameters

\gamma\in(0,1)

and

\ell\in\mathbb{N}

2repeat $\ell$ times

3 Run Quantum_Phase_Estimation on

\rho

and the simulated unitary

U_{k}

in Definition 30 with precision

\gamma

and parameter

k=\lceil\frac{200\log(d\log(1/\gamma)/\gamma))}{\gamma^{2}}\rceil

5 end

Post-select the most frequently observed estimates

\tilde{\lambda}_{i}

to obtain the state

\sum_{i=1}^{\mathrm{rk}\left(\rho\right)}\lambda_{i}\,\lvert\psi_{i}\rangle\langle\psi_{i}\rvert\otimes\lvert\tilde{\lambda}_{i}\rangle\langle\tilde{\lambda}_{i}\rvert

Algorithm 2 Quantum_Spectral_Sampling

4.4 Quantum eigenstate filtering

The following theorem is implicit in the work of Lin and Tong [LT20, Theorem 3] and states that we can approximately prepare a given eigenstate $\lvert\psi_{i}\rangle$ of a density operator $\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert$ up to precision $\varepsilon\in(0,1)$ in time $\operatorname{poly}\big{(}\Delta^{-1},1/\sqrt{\lambda_{i}},\log\left(1/\varepsilon\right),T_{\rho}\big{)}$ , where $r=\mathrm{rk}\left(\rho\right)$ is the rank of $\rho$ and where $T_{\rho}$ is the time it takes to prepare a purification if $\rho$ .

Theorem 32 ([LT20]).

Let $\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert\in\mathbb{C}^{d\times d}$ be a density operator with rank $r=\mathrm{rk}\left(\rho\right)$ and $\Delta$ -gapped spectrum $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{r})$ , for some $\Delta>0$ . Then, given purified access $U_{\rho}$ to $\rho$ , one can construct a unitary quantum circuit $U_{i}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\lambda_{i},\varepsilon)$ (and its inverse $U_{i}^{\dagger}$ ) with respect to an eigenvalue $\lambda_{i}$ and parameter $\varepsilon\in(0,1)$ which takes as input $\lvert 0^{\log d}\rangle$ and generates an approximate eigenstate $\lvert\tilde{\psi}_{i}\rangle$ in time $O\left(\frac{\log\left(\frac{1}{\varepsilon}\right)T_{\rho}}{\Delta\cdot\sqrt{\lambda_{i}}}\right)$ such that:

\|\lvert\tilde{\psi}_{i}\rangle-\lvert\psi_{i}\rangle\|\leq\varepsilon.

4.5 Technical lemmas

Lemma 33.

Let $\lvert\psi\rangle,\lvert\phi\rangle\in\mathbb{C}^{d}$ be pure states and $\sigma\in\mathbb{C}^{d\times d}$ a density matrix, and suppose that there exist pure states $\lvert\tilde{\psi}\rangle,\lvert\tilde{\phi}\rangle\in\mathbb{C}^{d}$ such that $\|\lvert\tilde{\psi}\rangle-\lvert\psi\rangle\|\leq\varepsilon$ and $\|\lvert\tilde{\phi}\rangle-\lvert\phi\rangle\|\leq\varepsilon$ , for $\varepsilon\in(0,1)$ . Then,

\big{|}\langle\tilde{\psi}\rvert\sigma\lvert\tilde{\phi}\rangle-\langle\psi\rvert\sigma\lvert\phi\rangle\big{|}\leq 2\varepsilon.

Proof.

Using that the states are normalized, we find that

$\displaystyle\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\tilde{\phi}\rangle-\langle\psi\rvert\sigma\lvert\phi\rangle\big{\|}$	$\displaystyle\leq\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\tilde{\phi}\rangle-\langle\tilde{\psi}\rvert\sigma\lvert\phi\rangle\big{\|}+\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\phi\rangle-\langle\psi\rvert\sigma\lvert\phi\rangle\big{\|}$
	$\displaystyle=\big{\|}\langle\tilde{\psi}\rvert\sigma(\lvert\tilde{\phi}\rangle-\lvert\phi\rangle)\big{\|}+\big{\|}(\langle\tilde{\psi}\rvert-\langle\psi\rvert)\sigma\lvert\phi\rangle\big{\|}$	$\displaystyle(\text{linearity})$
	$\displaystyle=\big{\|}\langle\tilde{\psi}\rvert\sigma(\lvert\tilde{\phi}\rangle-\lvert\phi\rangle)\big{\|}+\big{\|}\langle\phi\rvert\sigma(\lvert\tilde{\psi}\rangle-\lvert\psi\rangle)^{*}\big{\|}$	$\displaystyle(\text{skew symmetry})$
	$\displaystyle\leq\\|\tilde{\psi}\\|\cdot\\|\sigma\\|\cdot\\|\lvert\tilde{\phi}\rangle-\lvert\phi\rangle\\|+\\|\phi\\|\cdot\\|\sigma\\|\cdot\\|\lvert\tilde{\psi}\rangle-\lvert\psi\rangle\\|$	$\displaystyle(\text{Cauchy-Schwarz})$
	$\displaystyle\leq 2\varepsilon.$

∎

Lemma 34 (Estimation errors for diagonal terms).

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be density matrices, where $\rho$ has the spectral decomposition $\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert$ . Suppose that, for every $i\in[r]$ , there exist estimates $\tilde{\lambda}_{i}\in(0,1]$ and $\tilde{\sigma}_{ii}\in(0,1]$ as well as parameters $\gamma,\xi,\tau,\nu\in(0,1)$ such that

•

$|\tilde{\lambda}_{i}-\lambda_{i}|\leq\gamma$ with probability at least $1-\tau$ ,
•

$|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle|\leq\xi$ with probability at least $1-\nu$ .

Then, for $i\in[r]$ , it holds with probability at least $1-\tau-\nu$ :

\big{|}\tilde{\lambda}_{i}\tilde{\sigma}_{ii}-\lambda_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\big{|}\leq\gamma+\xi.

Proof.

Fix $i\in[r]$ . By the union bound we get that with probability at least $1-\tau-\nu$ :

$\displaystyle\big{\|}\tilde{\lambda}_{i}\tilde{\sigma}_{ii}-\lambda_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\big{\|}$	$\displaystyle=\|\tilde{\lambda}_{i}\tilde{\sigma}_{ii}-\tilde{\lambda}_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle+\tilde{\lambda}_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\lambda_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\|\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{triangle inequality})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\psi_{i}\\|\cdot\\|\sigma\lvert\psi_{i}\rangle\\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{Cauchy-Schwarz})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\psi_{i}\\|\cdot\\|\sigma\\|\cdot\\|\psi_{i}\\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{consistency of norm})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\sigma\\|_{1}\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{since }\\|\sigma\\|\leq\\|\sigma\\|_{1})$
	$\displaystyle\leq\varepsilon+\xi,$

where in the last last we used that $0\leq\tilde{\lambda}_{i}\leq 1$ , for all $i\in[r]$ , and that $\|\sigma\|_{1}=1$ . ∎

Lemma 35 (Estimation errors for off-diagonal terms).

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be density matrices, where $\rho$ has the spectral decomposition $\rho=\sum_{i=1}^{r}\lambda_{i}\lvert\psi_{i}\rangle\langle\psi_{i}\rvert$ . Let $t\in(0,1)$ and suppose that, for every $i,j\in[r]$ with $i<j$ , there exist $\tilde{\lambda}_{i},\tilde{\sigma}_{ii}\in(0,1]$ as well as parameters $\gamma,\tau,\zeta,\nu\in(0,1)$ such that

•

$|\tilde{\lambda}_{i}-\lambda_{i}|\leq\gamma$ with probability at least $1-\tau$ , and
•

$|\tilde{\sigma}_{ij}-\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle|\leq\zeta$ with probability at least $1-\nu$ .

Then, for any fixed pair $i<j$ , it holds with probability at least $1-\nu-2\tau$ :

\Big{|}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}\,\tilde{\sigma}_{ij}\Big{|}\leq\zeta+\frac{\gamma^{2}}{2\kappa^{2}},

where $\kappa$ is a lower bound on the smallest eigenvalue of $\rho$ .

Proof.

Let us first bound the additive error of the estimate $\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}$ . Recall that a simple first-order Taylor expansion for $f(x)=\sqrt{x}$ with additive error $\Delta x$ with respect to $x$ reveals that

\sqrt{x+\Delta x}=\sqrt{x}\cdot\left(1+\frac{1}{2}\frac{\Delta x}{x}\right)+O(\Delta x^{2}),\quad\quad\forall x\geq 0.

Using Lagrange’s remainder theorem, we can bound the remainder as $\frac{1}{8}\frac{\Delta x^{2}}{x^{2}}$ . Assuming that $\kappa\in(0,1)$ is a lower bound on the smallest eigenvalue of $\rho$ , we get with probability $1-\tau$ ,

\displaystyle|\sqrt{\tilde{\lambda}_{i}}-\sqrt{\lambda_{i}}|\leq\frac{\gamma^{2}}{8\kappa^{2}}.

(35)

Consequently, with probability at least $1-2\tau$ , we have

\displaystyle|\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}-\sqrt{{\lambda}_{i}}\sqrt{{\lambda}_{j}}|\leq\frac{\gamma^{2}}{4\kappa^{2}}+\frac{\gamma^{4}}{64\kappa^{4}}\leq\frac{\gamma^{2}}{2\kappa^{2}}.

(36)

where we used the fact that $\gamma\in(0,1)$ . Therefore, by the triangle inequality and (36), we obtain the following upper bound with probability at least $1-\nu-2\tau$ :

	$\displaystyle\Big{\|}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}\,\tilde{\sigma}_{ij}\Big{\|}$
	$\displaystyle\leq\Big{\|}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\,\tilde{\sigma}_{ij}\Big{\|}+\Big{\|}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\,\tilde{\sigma}_{ij}-\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}\,\tilde{\sigma}_{ij}\Big{\|}$
	$\displaystyle\leq\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\cdot\big{\|}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\tilde{\sigma}_{ij}\big{\|}+\tilde{\sigma}_{ij}\cdot\Big{\|}\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}-\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\Big{\|}$
	$\displaystyle\leq\zeta+\frac{\gamma^{2}}{2\kappa^{2}}$

This proves the claim. ∎

4.6 Bounds for the non-uniform coupon collector problem

The non-uniform coupon collector problem was first analyzed by Flajolet et al. [FGT92] and asks the following: Given a (possibly non-uniform) probability distribution $p=(p_{1},\dots,p_{n})$ with $p_{i}>0$ over a set of $n$ coupons, where the $i$ -th coupon is sampled independently with probability $p_{i}$ , how many independent samples from $p$ are necessary to obtain a full collection? If $p$ is the uniform distribution, the problem is equivalent to the well-known (standard) coupon collector problem.

Let $T_{(p_{1},\dots,p_{n})}$ be a random variable for the number of samples that are necessary to obtain a complete collection of all $n$ coupons. Our first result is the following non-trivial upper bound on the average waiting time $\mathbb{E}[T_{(p_{1},\dots,p_{n})}]$ with respect to the random variable $T_{(p_{1},\dots,p_{n})}$ .

Lemma 36.

Let $p=(p_{1},\dots,p_{n})$ be a distribution such that $p_{i}>0$ , for all $i\in[n]$ . Then,

\mathbb{E}[T_{(p_{1},\dots,p_{n})}]=\int_{0}^{\infty}\Big{(}1-\prod_{i=1}^{n}(1-e^{-p_{i}t})\Big{)}\,\mathrm{d}t\,\leq\,n\cdot H(p_{1},\dots,p_{n})^{-1},

where $H(x_{1},\dots,x_{n})=n/(\sum_{i=1}^{n}x_{i}^{-1})$ is the harmonic mean.

Proof.

Using an identity due to Flajolet et al. [FGT92], we obtain

$\displaystyle\mathbb{E}[T_{(p_{1},\dots,p_{n})}]$	$\displaystyle=\int_{0}^{\infty}\Big{(}1-\prod_{i=1}^{n}(1-e^{-p_{i}t})\Big{)}\,\mathrm{d}t$
	$\displaystyle\leq\int_{0}^{\infty}\Big{(}\sum_{i=1}^{n}e^{-p_{i}t}\Big{)}\,\mathrm{d}t$	$\displaystyle(\text{Weierstrass product inequality})$
	$\displaystyle=\sum_{i=1}^{n}\int_{0}^{\infty}e^{-p_{i}t}\,\mathrm{d}t$	$\displaystyle(\text{switching order of summation and integration})$
	$\displaystyle=\sum_{i=1}^{n}\left[-\frac{e^{-p_{i}t}}{p_{i}}\right]_{0}^{\infty}$
	$\displaystyle=\sum_{i=1}^{n}p_{i}^{-1}\quad=\quad n\cdot H(p_{1},\dots,p_{n})^{-1}.$	$\displaystyle(\text{by definition})$

∎

We can now apply the previous lemma in the context of spectral sampling as follows. Let $\rho\in\mathbb{C}^{d\times d}$ be a density matrix with spectrum $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{r})$ . Suppose we have access to a subroutine as in (4) that allows us to sample random pairs of eigenstates and eigenvalues $(\lvert\psi_{i}\rangle,\tilde{\lambda}_{i})$ with probability $\lambda_{i}\in(0,1]$ , with the promise that each approximate eigenvalue $\tilde{\lambda}_{i}$ is sufficiently close to $\lambda_{i}$ and separated from the remaining spectrum of $\rho$ . Let $\mathbb{E}[T_{\operatorname{spec}(\rho)}]$ denote the average number of repetitions needed to obtain a full collection of $r$ distinct eigenvalues.

Applying the inequality in Lemma 36 in the context of spectral sampling of a density operator $\rho\in\mathbb{C}^{d\times d}$ , we obtain the following upper bound which directly relates the average time to complete the collection to the spectral properties of $\rho$ .

Corollary 37.

Let $\rho\in\mathbb{C}^{d\times d}$ be a density matrix and let $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{r})$ with $r=\mathrm{rk}\left(\rho\right)$ . Then,

\mathbb{E}[T_{\operatorname{spec}(\rho)}]=\int_{0}^{\infty}\Big{(}1-\prod_{i=1}^{r}(1-e^{-\lambda_{i}t})\Big{)}\,\mathrm{d}t\,\leq\,r\cdot H(\operatorname{spec}(\rho))^{-1}.

Unfortunately, the above result is not tight. In particular, for the uniform spectrum $(\frac{1}{r},\dots,\frac{1}{r})$ , our bound tells us that $\mathbb{E}[T_{(\frac{1}{r},\dots,\frac{1}{r})}]\leq r^{2}$ , whereas a well known result on the (standard) uniform coupon collector problem states that the average number of draws is in the order of $\Theta(r\log r)$ .

In order to find an improved bound which asymptotically matches the standard coupon collector result, we use a coupling argument which allows us to relate an instance of the non-uniform coupon collector problem to a worst-case instance of the uniform coupon collector problem.

To this end, it is convenient to formalize the coupon collector problem with respect to the uniform distribution $\mathcal{U}([0,1])$ on the interval $[0,1]$ as follows. Let $p=(p_{1},\dots,p_{n})$ be a (possibly non-uniform) distribution with $p_{i}>0$ , for all $i\in[n]$ , and consider the following experiment:

1.

Sample $x\sim\mathcal{U}([0,1])$ according to the uniform distribution over $[0,1]$ .
2.

Assign $x$ to a coupon by bucketing it according to the probability distribution $p$ .

The problem then becomes: How many samples from $\mathcal{U}([0,1])$ are needed for a full collection? It is easy to see that the experiment is equivalent to the standard (non-uniform) coupon collector problem.

As an example, consider the uniform distribution $p=(\frac{1}{n},\dots,\frac{1}{n})$ . In this case, we bucket using the function $f(x)=\lceil x\cdot n\rceil$ , which means that we obtain the $i$ -th coupon with probability

\underset{{x\sim\mathcal{U}([0,1])}}{\mbox{Pr}}\left[f(x)=i\right]=\frac{1}{n}.

In the non-uniform case, $p$ implicitly defines a partition of $[0,1]$ since $\sum_{i=1}^{n}p_{i}=1$ . We let $g$ define the function that maps $x$ to a unique coupon by bucketing it into the respective interval in $[0,1]$ .

We prove the following result.

Lemma 38 (Worst-case bound for the non-uniform coupon collector’s problem).

Let $p=(p_{1},\dots,p_{n})$ be a probability distribution with $p_{1}\geq p_{2}\geq\dots\geq p_{n}>0$ and let $m$ be the smallest integer such that $p_{n}\geq\frac{1}{2m}$ . Then, we obtain the following upper bound to sample a full collection

\mathbb{E}\big{[}T_{(p_{1},\dots,p_{n})}\big{]}\leq\mathbb{E}\big{[}T_{(\frac{1}{m},\dots,\frac{1}{m})}\big{]}.

Proof.

We use a coupling argument by analyzing the two instances of the coupon collector over the probability measure $\mathcal{U}([0,1])$ . By using $f$ and $g$ from before, we can see that the marginal distributions of our coupling match the original instances. In fact, the $i$ -th coupon is selected with probability

	$\displaystyle\underset{{x\sim\mathcal{U}([0,1])}}{\mbox{Pr}}\left[g(x)=i\right]=p_{i}$	$\displaystyle(\text{non-uniform distribution})$
	$\displaystyle\underset{{x\sim\mathcal{U}([0,1])}}{\mbox{Pr}}\left[f(x)=i\right]=\frac{1}{m}$	$\displaystyle(\text{uniform distribution})$

Therefore, our aforementioned coupling over the probability measure $\mathcal{U}([0,1])$ is well-defined. Let $X_{t}$ be the event that the collection is incomplete at step $t-1$ for the distribution $p$ . Similarly, we define $Y_{t}$ to be the event that the collection is incomplete at step $t-1$ for the distribution $q=(\frac{1}{m},\dots,\frac{1}{m})$ .

We now claim that $X_{t}\subseteq Y_{t}$ , for every $t\geq 1$ . In other words, if the collection is incomplete with respect to $p$ , then it must be also be incomplete with respect to $q$ . Suppose that event $X_{t}$ occurs, hence there exists a coupon $j$ which has not been collected. Because of the coupling this means that there exists an interval $I_{j}\subset[0,1]$ for which no $x$ has appeared such that $x\in I_{j}$ . Now, because $\min_{i\in[n]}p_{i}\geq\frac{1}{2m}$ , it follows that $I_{j}$ contains at least one interval of size $\frac{1}{m}$ in the range of $f$ . Hence, there exists a coupon for the distribution $q$ which has not been collected. This proves the claim.

We can now show the following upper bound:

$\displaystyle\mathbb{E}\big{[}T_{(p_{1},\dots,p_{n})}\big{]}$	$\displaystyle=\mathbb{E}\left[\sum_{t=1}^{\infty}\mathbb{I}_{X_{t}}\right]$	$\displaystyle(\text{by definition})$
	$\displaystyle=\sum_{t=1}^{\infty}\mathbb{E}[\mathbb{I}_{X_{t}}]$	$\displaystyle(\text{linearity of expectation})$
	$\displaystyle=\sum_{t=1}^{\infty}\mbox{Pr}[X_{t}]$
	$\displaystyle\leq\sum_{t=1}^{\infty}\mbox{Pr}[Y_{t}]$	$\displaystyle(\text{using that }X_{t}\subseteq Y_{t})$
	$\displaystyle=\mathbb{E}\big{[}T_{(\frac{1}{m},\dots,\frac{1}{m})}\big{]}.$

∎

The following corollary is a simple consequence of Lemma 38.

Corollary 39 (Worst-case bound for the thresholded non-uniform coupon collector’s problem).

Let $p=(p_{1},\dots,p_{n})$ be an arbitrary probability distribution with $p_{1}\geq p_{2}\geq\dots\geq p_{n}>0$ and let $\theta\in(0,1)$ be a threshold value and let $m=\lceil\frac{1}{2\theta}\rceil$ . Let $T_{(p_{1},\dots,p_{n})}^{\theta}$ be the the random variable for the number of repetitions it takes to sample all coupons which occur with probability at least $\theta$ . Then,

\mathbb{E}\big{[}T_{(p_{1},\dots,p_{n})}^{\theta}\big{]}\leq\mathbb{E}\big{[}T_{(\frac{1}{m},\dots,\frac{1}{m})}\big{]}.

4.7 Analysis of the Algorithm

Let us now analyze our spectral sampling-based algorithm for (truncated) fidelity estimation.

Theorem 40.

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be arbitrary density matrices, and suppose that $\rho$ has a well-separated spectrum with a gap $\Delta>0$ . Suppose also that $\rho$ has the smallest rank of the two states. Let $\varepsilon\in(0,1)$ , $\delta\in(0,1)$ , $\theta\in(0,1)$ be parameters. Then, Algorithm 1 runs in time

\tilde{O}\left(\frac{T_{\rho}+T_{\sigma}}{\theta^{10.5}\varepsilon^{4}\Delta}+\frac{T_{\rho}}{\delta\theta^{3}\gamma^{3}}\right)

where $\gamma=\min\{\frac{\theta^{3}\varepsilon}{\sqrt{2}},\frac{\Delta}{2}\}$ and outputs an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that with probability $1-\delta$ :

|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon,

where $\rho_{\theta}$ is a “soft-thresholded” version of $\rho$ in which eigenvalues of $\rho$ below $\theta/2$ are completely removed and those above $\theta$ are kept intact, while eigenvalues in $[\theta/2,\theta]$ are decreased by some amount.

Proof.

Let $\varepsilon\in(0,1)$ and $\delta\in(0,1)$ be parameters, and let $\theta\in(0,1)$ be the truncation parameter. Let $\operatorname{spec}(\rho)=(\lambda_{1},\dots,\lambda_{\mathrm{rk}\left(\rho\right)})$ be the eigenvalue spectrum of $\rho\in\mathbb{C}^{d\times d}$ and let $\mathrm{rk}_{\theta}$ be the number of eigenvalues in the interval $(\theta,1]$ . Let $m=\lceil\frac{1}{2\theta}\rceil$ be a parameter and note that $\mathrm{rk}_{\theta}\leq\frac{1}{\theta}$ .

We show that it suffices to run Quantum_Spectral_Sampling $(U_{\rho},\gamma,\ell)$ at most $M=\lceil\frac{16mH_{m}}{\delta\cdot\theta^{2}}\rceil$ times with parameters $\gamma=\min\{\frac{\theta^{3}\varepsilon}{\sqrt{2}},\frac{\Delta}{2}\}$ and $\ell=\lceil\log\big{(}\frac{16}{\delta\cdot\theta^{2}}\big{)}+\log M\rceil$ to find $r_{\theta}\geq\mathrm{rk}_{\theta}$ accurate estimates $\tilde{\lambda}_{i}\geq\frac{3\theta}{4}$ of distinct eigenvalues $\lambda_{i}\in[\frac{\theta}{2},1]$ , including the complete set of eigenvalues in the interval $(\theta,1]$ , with probability at least $1-\frac{\delta\cdot\theta^{2}}{8}$ . Let us now analyze the procedure in more detail.

In each iteration of Quantum_Spectral_Sampling $(U_{\rho},\gamma,\ell)$ , we obtain an eigenvalue estimate $\tilde{\lambda}_{i}$ such that $|\tilde{\lambda}_{i}-\lambda_{i}|\leq\gamma$ with probability $1-2^{-l}$ . Hence, by the union bound, all $M$ trials will produce accurate estimates with probability at least $1-\frac{\delta\cdot\theta^{2}}{16}$ . Moreover, by our choice of parameters, each eigenvalue estimate $\tilde{\lambda}_{i}$ falls within distance $\gamma<\frac{\Delta}{2}$ of the eigenvalue $\lambda_{i}$ , enabling Algorithm 1 to uniquely identify each sample. Note also that, since $\gamma<\frac{\theta}{4}$ , Algorithm 1 never accepts any eigenvalue estimates $\tilde{\lambda}_{i}\geq\frac{3\theta}{4}$ that correspond to eigenvalues $\lambda_{i}$ in the interval $(0,\frac{\theta}{2})$ , thus implying that $r_{\theta}\leq\frac{2}{\theta}$ .

Let us now bound the probability of error, i.e. the probability of not finding a complete set of at least $\mathrm{rk}_{\theta}$ distinct samples, including all eigenvalues in $(\theta,1]$ , after $M$ trials. Let $T_{(\lambda_{1},\dots,\lambda_{\mathrm{rk}(\rho)})}^{\theta}$ be a random variable for the number of samples needed to obtain such a collection. By Markov’s inequality,

\displaystyle\mbox{Pr}\left[T_{(\lambda_{1},\dots,\lambda_{\mathrm{rk}(\rho)})}^{\theta}\geq\frac{16mH_{m}}{\delta\cdot\theta^{2}}\right]\leq\frac{\delta\cdot\theta^{2}}{16mH_{m}}\mathbb{E}[T_{(\lambda_{1},\dots,\lambda_{\mathrm{rk}(\rho)})}^{\theta}].

(37)

Let $T_{(\frac{1}{m},\dots,\frac{1}{m})}$ denote the number of repetitions needed to draw $m$ distinct coupons in the uniform coupon collector problem, where $m=\lceil\frac{1}{2\theta}\rceil$ . From Corollary 39 it follows that

\mathbb{E}[T_{(\lambda_{1},\dots,\lambda_{\mathrm{rk}(\rho)})}^{\theta}]\,\leq\,\mathbb{E}[T_{(\frac{1}{m},\dots,\frac{1}{m})}]=m\cdot H_{m},

where $H_{m}$ is the $m$ -th harmonic number and $m\cdot H_{m}=\Theta(m\log m)$ . Hence, the probability in (37) is at most $\frac{\delta\cdot\theta^{2}}{16}$ , as required. Therefore, by the union bound, the quantum spectral sampling procedure of Algorithm 1 succeeds at collecting $r_{\theta}\geq\mathrm{rk}_{\theta}$ both distinct and accurate eigenvalue estimates after

M=\left\lceil\frac{16\cdot mH_{m}}{\delta\cdot\theta^{2}}\right\rceil=\Theta\left(\frac{\log\big{(}\theta^{-1}\big{)}}{\delta\cdot\theta^{3}}\right)

repetitions of Quantum_Spectral_Sampling $(U_{\rho},\gamma,\ell)$ with probability at least $1-\frac{\delta\cdot\theta^{2}}{8}$ . Because $r_{0}\leq\frac{2}{\theta}$ , this in turn implies that the spectral sampling procedure succeeds with probability at least $1-\frac{\delta}{2r_{0}^{2}}$ .

Recall that $\Lambda(\rho,\sigma)=\sqrt{\rho}\sigma\sqrt{\rho}\in\mathbb{C}^{d\times d}$ has the following non-trivial entries in the eigenbasis of $\rho$ :

\displaystyle\Lambda_{ij}=\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}\langle\psi_{i}\rvert\sigma\lvert\psi_{j}\rangle,\quad\quad\forall i,j\in\{1,\dots,\mathrm{rk}\left(\rho\right)\}.

We show that Algorithm 1 obtains an estimate $\hat{\Lambda}^{\theta}$ of a matrix $\Lambda^{\theta}=\Lambda(\rho_{\theta},\sigma)$ , where $\rho_{\theta}$ is a soft truncation of $\rho$ in which eigenvalues of $\rho$ below $\theta/2$ are completely removed and those above $\theta$ are kept intact, while eigenvalues in $[\theta/2,\theta]$ are potentially incomplete or decreased by some amount.

Let us now consider the estimates for the diagonal entries of the matrix $\Lambda^{\theta}$ (lines $11$ and $12$ ). Let $i\in[r_{\theta}]$ be an index. To estimate $\Lambda_{ii}^{\theta}$ , we first run Swap_Test $(U_{i}\lvert 0^{\log d}\rangle,\sigma,\xi,\nu)$ with parameters $\xi=\frac{\varepsilon^{2}}{8r_{\theta}^{4}}$ and $\nu=\frac{\delta}{2r_{\theta}^{2}}$ to obtain an estimate $\tilde{\sigma}_{ii}$ such that $|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle|\leq\xi$ with probability at least $1-\nu$ (by Lemma 33), where $U_{i}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{i},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})$ is the unitary in Theorem 32. Letting $\hat{\Lambda}_{ii}^{\theta}\leftarrow\tilde{\lambda}_{i}\tilde{\sigma}_{ii}$ , it then follows from Lemma 34 that for every index $i\in[r_{\theta}]$ ,

\displaystyle\mbox{Pr}\left[\Big{|}{\Lambda}_{ii}^{\theta}-\hat{\Lambda}_{ii}^{\theta}\Big{|}\leq\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]\geq 1-\frac{\delta}{r_{\theta}^{2}}.

(38)

Let us now consider the estimates for the off-diagonal entries of the matrix $\Lambda^{\theta}$ (lines $14$ to $16$ ). Let $i,j\in[r_{\theta}]$ be a pair of indices with $i<j$ . We run Hadamard_Test $(\lvert 0^{\log d}\rangle,U_{i}^{\dagger}UU_{j},\xi,\nu)$ with parameters $\xi=\frac{\varepsilon^{2}}{8r_{\theta}^{4}}$ and $\nu=\frac{\delta}{2r_{\theta}^{2}}$ to obtain an estimate $\tilde{\sigma}_{ij}$ , where $U_{j}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{j},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})$ and $U_{i}^{\dagger}=\texttt{Quantum\_Eigenstate\_Filtering}(U_{\rho},\tilde{\lambda}_{i},\frac{\varepsilon^{2}}{16r_{\theta}^{4}})^{\dagger}$ are unitary as in Theorem 32, and where $U$ is a block-encoding of $\sigma$ as in Lemma 7. Letting $\hat{\Lambda}_{ij}^{\theta}\leftarrow\sqrt{\tilde{\lambda}_{i}}\sqrt{\tilde{\lambda}_{j}}\,\tilde{\sigma}_{ij}$ and $\hat{\Lambda}^{\theta}_{ji}=\hat{\Lambda}^{\theta^{*}}_{ij}$ for $i\neq j$ , we then have by Lemma 35 that

\mbox{Pr}\left[\Big{|}{\Lambda}_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\Big{|}\leq\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]\geq 1-\frac{\delta}{r_{\theta}^{2}}.

Given our choice of parameters $\varepsilon$ and $\delta$ , we have that for every fixed $i,j\in[r_{\theta}]$ :

\mbox{Pr}\left[\Big{|}{\Lambda}_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\Big{|}>\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]\leq\frac{\delta}{r_{\theta}^{2}}.

Using the continuity bound for fidelity estimation (Theorem 29), we obtain

$\displaystyle\|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)\|$	$\displaystyle=\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda^{\theta}}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}^{\theta}}}\big{]}\big{\|}$	(39)
	$\displaystyle\leq\sqrt{2}r_{\theta}\cdot\sqrt{\\|\Lambda^{\theta}-\hat{\Lambda}^{\theta}\\|_{1}}$
	$\displaystyle\leq\sqrt{2}r_{\theta}^{2}\cdot\sqrt{\\|\Lambda^{\theta}-\hat{\Lambda}^{\theta}\\|_{\max}}$
	$\displaystyle=\sqrt{2}r_{\theta}^{2}\cdot\sqrt{\max_{1\leq i,j\leq r_{\theta}}\|\Lambda_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\|}.$	(40)

Putting everything together and using the inequality in (40), we find that

$\displaystyle\mbox{Pr}\Big{[}\|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)\|\leq\varepsilon\Big{]}$	$\displaystyle=\mbox{Pr}\Big{[}\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda^{\theta}}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}^{\theta}}}\big{]}\big{\|}\leq\varepsilon\Big{]}$	$\displaystyle(\text{by definition})$
	$\displaystyle\geq\mbox{Pr}\left[\max_{1\leq i,j\leq r_{\theta}}\,\Big{\|}\Lambda_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\Big{\|}\leq\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]$
	$\displaystyle\geq\mbox{Pr}\left[\forall i,\forall j\,:\,\Big{\|}\Lambda_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\Big{\|}\leq\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]$
	$\displaystyle\geq 1-\sum_{i=1}^{r_{\theta}}\sum_{j=1}^{r_{\theta}}\mbox{Pr}\left[\Big{\|}\Lambda_{ij}^{\theta}-\hat{\Lambda}_{ij}^{\theta}\Big{\|}>\frac{\varepsilon^{2}}{2r_{\theta}^{4}}\right]$	$\displaystyle(\text{union bound})$
	$\displaystyle\geq 1-r_{\theta}^{2}\cdot\delta/r_{\theta}^{2}\quad=\quad 1-\delta.$

This proves the claim. ∎

Finally, we obtain the following result.

Theorem 41.

Let $\rho,\sigma\in\mathbb{C}^{d\times d}$ be arbitrary density matrices, and suppose that $\rho$ has a well-separated spectrum with a gap $\Delta>0$ and that $\rho$ has the smallest rank with $\mathrm{rk}\left(\rho\right)\leq r$ . Let $\varepsilon\in(0,1)$ and $\delta\in(0,1)$ , and fix $\theta=\frac{\varepsilon^{2}}{4r}$ . Then, Algorithm 1 with parameters $\varepsilon^{\prime}=\frac{\varepsilon}{2}$ , $\delta$ and $\theta$ runs in time

\tilde{O}\left(\frac{r^{10.5}}{\varepsilon^{25}\Delta}\cdot\big{(}T_{\rho}+T_{\sigma}\big{)}+\frac{r^{3}T_{\rho}}{\delta\min\{\frac{\varepsilon^{7}}{r^{3}},\Delta\}^{3}}\right)

and outputs an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that with probability $1-\delta$ :

|F(\rho,\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\varepsilon.

Proof.

Given the set of parameters $\varepsilon^{\prime}=\frac{\varepsilon}{2}\in(0,1)$ , $\delta\in(0,1)$ and $\theta=\frac{\varepsilon^{2}}{4r}$ , it follows from Theorem 40 that Algorithm 1 outputs an estimate $\hat{F}(\rho_{\theta},\sigma)$ such that with probability $1-\delta$ :

|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)|\leq\frac{\varepsilon}{2},

	$\displaystyle\|F(\rho,\sigma)-\hat{F}(\rho_{\theta},\sigma)\|$	$\displaystyle\leq\|F(\rho,\sigma)-F(\rho_{\theta},\sigma)\|+\|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)\|$
		$\displaystyle\leq\sqrt{\mbox{Tr}\left[\Pi_{[0,\theta)}\rho\Pi_{[0,\theta)}\right]}+\|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)\|$
		$\displaystyle\leq\sqrt{r\cdot\theta}+\|F(\rho_{\theta},\sigma)-\hat{F}(\rho_{\theta},\sigma)\|$
		$\displaystyle\leq\frac{\varepsilon}{2}+\frac{\varepsilon}{2}$
		$\displaystyle=\varepsilon.$

This proves the claim. ∎

5 Discussion

We give an efficient and versatile algorithm for (truncated) fidelity estimation. Our algorithm demonstrates the potential of block-encoding and quantum singular value transformation techniques for quantum information processing tasks. We also demonstrate and work out the specifics of a generic method suggested by [GLM⁺20] for utilizing these techniques in the scenario when one only has access to copies of the states. This method might be of independent interest.

Our alternative spectral-sampling-based algorithm for fidelity estimation performs significantly worse in general compared to our block-encoding algorithm, however it may be easier to implement in certain settings, e.g., when it is easy to obtain circuits that prepare the eigenstates of one of the density operators. For example, consider the problem of exactly simulating the one-dimensional Ising chain [CL18]. There, it is possible to efficiently prepare all eigenstates of the Ising Hamiltonian without relying on quantum eigenstate filtering.

Acknowledgements

A.P. would like to thank Ronald de Wolf and Thomas Vidick for many useful suggestions that helped improve the spectral sampling algorithm.

Part of this work was done while the authors visited the Simons Institute for the Theory of Computing; we gratefully acknowledge its hospitality.

References

[AJL06] Dorit Aharonov, Vaughan Jones, and Zeph Landau. A polynomial quantum algorithm for approximating the Jones polynomial. In Proceedings of the 38th ACM Symposium on the Theory of Computing (STOC), pages 427–436, 2006. arXiv: quant-ph/0511096
[vAGGdW20] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf. Quantum SDP-solvers: Better upper and lower bounds. Quantum, 4:230, 2020. Earlier version in FOCS’17. arXiv: 1705.01843
[ARSW21] Rochisha Agarwal, Soorya Rethinasamy, Kunal Sharma, and Mark M. Wilde. Estimating distinguishability measures on quantum computers, 2021. arXiv: 2108.08406
[BCC⁺15] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma. Simulating Hamiltonian dynamics with a truncated Taylor series. Physical Review Letters, 114(9):090502, 2015. arXiv: 1412.4687
[BCWdW01] Harry Buhrman, Richard Cleve, John Watrous, and Ronald de Wolf. Quantum fingerprinting. Physical Review Letters, 87(16):167902, 2001. arXiv: quant-ph/0102001
[Bha97] Rajendra Bhatia. Matrix Analysis, volume 169 of Graduate Texts in Mathematics. Springer, 1997.
[BOW19] Costin Bădescu, Ryan O’Donnell, and John Wright. Quantum state certification. In Proceedings of the 51st ACM Symposium on the Theory of Computing (STOC), pages 503--514, 2019. arXiv: 1708.06002
[CGJ19] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery. The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), pages 33:1--33:14, 2019. arXiv: 1804.01973
[CKS17] Andrew M. Childs, Robin Kothari, and Rolando D. Somma. Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM Journal on Computing, 46(6):1920--1950, 2017. arXiv: 1511.02306
[CL18] Alba Cervera-Lierta. Exact Ising model simulation on a quantum computer. Quantum, 2:114, 2018. arXiv: 2103.09076
[CPCC20] Marco Cerezo, Alexander Poremba, Lukasz Cincio, and Patrick J. Coles. Variational quantum fidelity estimation. arXiv: 1906.09253, 2020.
[CSSC18] Lukasz Cincio, Yiğit Subaşı, Andrew T Sornborger, and Patrick J Coles. Learning the quantum algorithm for state overlap. New Journal of Physics, 20(11):113022, 2018. arXiv: 1803.04114
[EY07] Alexandre Eremenko and Peter Yuditskii. Uniform approximation of sgn(x) by polynomials and entire functions. Journal d’Analyse Mathématique, 101(1):313--324, 2007. arXiv: math/0604324
[FGT92] Philippe Flajolet, Danièle Gardy, and Loÿs Thimonier. Birthday paradox, coupon collectors, caching algorithms and self-organizing search. Discrete Applied Mathematics, 39(3):207--229, 1992.
[FL11] Steven T. Flammia and Yi-Kai Liu. Direct fidelity estimation from few pauli measurements. Phys. Rev. Lett., 106:230501, Jun 2011.
[Gil19] András Gilyén. Quantum Singular Value Transformation & Its Algorithmic Applications. PhD thesis, University of Amsterdam, 2019.
[GLM⁺20] András Gilyén, Seth Lloyd, Iman Marvian, Yihui Quek, and Mark M. Wilde. Quantum algorithm for Petz recovery channels and pretty good measurements. arXiv: 2006.16924, 2020.
[GSLW18] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics [full version], 2018. arXiv: 1806.01838
[GSLW19] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st ACM Symposium on the Theory of Computing (STOC), pages 193--204, 2019. arXiv: 1806.01838
[HHL09] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd. Quantum algorithm for linear systems of equations. Physical Review Letters, 103(15):150502, 2009. arXiv: 0811.3171
[HJ90] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 1990.
[HP14] Fumio Hiai and Dénes Petz. Introduction to Matrix Analysis and Applications. Universitext. Springer, 2014.
[KLL⁺17] Shelby Kimmel, Cedric Yen-Yu Lin, Guang Hao Low, Maris Ozols, and Theodore J. Yoder. Hamiltonian simulation with optimal sample complexity. npj Quantum Information, 3(1):13, 2017. arXiv: 1608.00281
[LC17] Guang Hao Low and Isaac L. Chuang. Hamiltonian simulation by uniform spectral amplification. arXiv: 1707.05391, 2017.
[LC19] Guang Hao Low and Isaac L. Chuang. Hamiltonian simulation by qubitization. Quantum, 3:163, 2019. arXiv: 1610.06546
[LMR14] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost. Quantum principal component analysis. Nature Physics, 10:631--633, 2014. arXiv: 1307.0401
[LT20] Lin Lin and Yu Tong. Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems. Quantum, 4:361, 2020. arXiv: 1910.14596
[NC00] Michael A. Nielsen and Isaac L. Chuang. Quantum computation and quantum information. Cambridge University Press, 2000.
[OW15] Ryan O’Donnell and John Wright. Quantum spectrum testing. In Proceedings of the 47th ACM Symposium on the Theory of Computing (STOC), pages 529--538, 2015. arXiv: 1501.05028
[OW16] Ryan O’Donnell and John Wright. Efficient quantum tomography. In Proceedings of the 48th ACM Symposium on the Theory of Computing (STOC), pages 899--912, 2016. arXiv: 1508.01907
[Pra14] Anupam Prakash. Quantum Algorithms for Linear Algebra and Machine Learning. PhD thesis, University of California at Berkeley, 2014.
[PS70] Robert T. Powers and Erling Størmer. Free states of the canonical anticommutation relations. Communications in Mathematical Physics, 16(1):1 -- 33, 1970.
[QC09] H. T. Quan and F. M. Cucchietti. Quantum fidelity and thermal phase transitions. Phys. Rev. E, 79:031101, Mar 2009.
[Wat02] John Watrous. Limits on the power of quantum statistical zero-knowledge. In Proceedings of the 43rd IEEE Symposium on Foundations of Computer Science (FOCS), pages 459--468, 2002. arXiv: quant-ph/0202111
[WGL⁺22] Qisheng Wang, Ji Guan, Junyi Liu, Zhicheng Zhang, and Mingsheng Ying. New quantum algorithms for computing quantum entropies and distances, 2022. arXiv: 2203.13522
[WZC⁺21] Qisheng Wang, Zhicheng Zhang, Kean Chen, Ji Guan, Wang Fang, and Mingsheng Ying. Quantum algorithm for fidelity estimation, 2021. arXiv: 2103.09076
[XCPdC21] Zhenyu Xu, Aurelia Chenu, Toma ž Prosen, and Adolfo del Campo. Thermofield dynamics: Quantum chaos versus decoherence. Phys. Rev. B, 103:064309, Feb 2021.
[ZQWS07] Paolo Zanardi, H. T. Quan, Xiaoguang Wang, and C. P. Sun. Mixed-state fidelity and quantum criticality at finite temperature. Physical Review A, 75:032109, Mar 2007.

$\displaystyle\left\|\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(\mathrm{sgn})}\!(Q^{\dagger})\right]-\mbox{Tr}\left[Q\cdot\mathrm{SV}^{(p)}\!(Q^{\dagger}/8)\right]\right\|$	$\displaystyle=\left\|\sum_{i=1}\varsigma_{i}-\varsigma_{i}p(\varsigma_{i}/8)\right\|$
	$\displaystyle\leq\sum_{i=1}\varsigma_{i}\|1-p(\varsigma_{i}/8)\|$
	$\displaystyle=\kern-2.84526pt\sum_{\varsigma_{i}\geq\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\kern-8.53581pt\varsigma_{i}\|1-p(\varsigma_{i}/8)\|+\kern-2.84526pt\sum_{\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\kern-8.53581pt\varsigma_{i}\|1-p(\varsigma_{i}/8)\|$
	$\displaystyle\leq\sum_{\varsigma_{i}\geq\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}\varsigma_{i}\frac{5\varepsilon}{60}+\sum_{\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}2\varsigma_{i}$
	$\displaystyle\leq\mbox{Tr}\left[\Sigma\right]\frac{\varepsilon}{10}+\sum_{0<\varsigma_{i}<\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}}2\frac{\varepsilon}{20\mathrm{rk}\left(t(\rho)\right)}\leq\frac{\varepsilon}{5}.$	(26)

$\displaystyle\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}^{+}}}\big{]}\big{\|}$	$\displaystyle=\big{\|}\sum_{i=1}^{r}\sqrt{\lambda_{i}}-\sum_{i=1}^{r}\sqrt{\hat{\lambda}^{+}_{i}}\big{\|}$	$\displaystyle(\text{by definition})$
	$\displaystyle\leq\sum_{i=1}^{r}\big{\|}\sqrt{\lambda_{i}}-\sqrt{\hat{\lambda}_{i}^{+}}\big{\|}$	$\displaystyle(\text{triangle inequality})$
	$\displaystyle\leq r\cdot\max_{1\leq i\leq r}\big{\|}\sqrt{\lambda_{i}}-\sqrt{\hat{\lambda}_{i}^{+}}\big{\|}$
	$\displaystyle\leq r\cdot\left\\|\sqrt{\Lambda}-\sqrt{\hat{\Lambda}_{+}}\right\\|_{2}$	$\displaystyle(\text{\cite[cite]{[\@@bibref{}{bhatia1997MatrixAnalysis}{}{}, Problem III.6.13]}})$
	$\displaystyle\leq r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}_{+}\right\\|_{1}}$	$\displaystyle(\text{\cite[cite]{[\@@bibref{}{cmp/1103842028}{}{}, Lemma 4.1]}})$

	$\displaystyle\|F(\rho,\sigma)-\hat{F}(\rho,\sigma)\|$	$\displaystyle=\big{\|}\mathrm{Tr}\big{[}{\sqrt{\Lambda}}\big{]}-\mathrm{Tr}\big{[}{\sqrt{\hat{\Lambda}_{+}}}\big{]}\big{\|}$
		$\displaystyle\leq r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}_{+}\right\\|_{1}}$
		$\displaystyle=\sqrt{2}r\cdot\sqrt{\left\\|\Lambda-\hat{\Lambda}\right\\|_{1}}.\quad\quad\quad$

$\displaystyle\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\tilde{\phi}\rangle-\langle\psi\rvert\sigma\lvert\phi\rangle\big{\|}$	$\displaystyle\leq\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\tilde{\phi}\rangle-\langle\tilde{\psi}\rvert\sigma\lvert\phi\rangle\big{\|}+\big{\|}\langle\tilde{\psi}\rvert\sigma\lvert\phi\rangle-\langle\psi\rvert\sigma\lvert\phi\rangle\big{\|}$
	$\displaystyle=\big{\|}\langle\tilde{\psi}\rvert\sigma(\lvert\tilde{\phi}\rangle-\lvert\phi\rangle)\big{\|}+\big{\|}(\langle\tilde{\psi}\rvert-\langle\psi\rvert)\sigma\lvert\phi\rangle\big{\|}$	$\displaystyle(\text{linearity})$
	$\displaystyle=\big{\|}\langle\tilde{\psi}\rvert\sigma(\lvert\tilde{\phi}\rangle-\lvert\phi\rangle)\big{\|}+\big{\|}\langle\phi\rvert\sigma(\lvert\tilde{\psi}\rangle-\lvert\psi\rangle)^{*}\big{\|}$	$\displaystyle(\text{skew symmetry})$
	$\displaystyle\leq\\|\tilde{\psi}\\|\cdot\\|\sigma\\|\cdot\\|\lvert\tilde{\phi}\rangle-\lvert\phi\rangle\\|+\\|\phi\\|\cdot\\|\sigma\\|\cdot\\|\lvert\tilde{\psi}\rangle-\lvert\psi\rangle\\|$	$\displaystyle(\text{Cauchy-Schwarz})$
	$\displaystyle\leq 2\varepsilon.$

$\displaystyle\big{\|}\tilde{\lambda}_{i}\tilde{\sigma}_{ii}-\lambda_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\big{\|}$	$\displaystyle=\|\tilde{\lambda}_{i}\tilde{\sigma}_{ii}-\tilde{\lambda}_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle+\tilde{\lambda}_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle-\lambda_{i}\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\|\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{triangle inequality})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\psi_{i}\\|\cdot\\|\sigma\lvert\psi_{i}\rangle\\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{Cauchy-Schwarz})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\psi_{i}\\|\cdot\\|\sigma\\|\cdot\\|\psi_{i}\\|\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{consistency of norm})$
	$\displaystyle\leq\tilde{\lambda}_{i}\cdot\|\tilde{\sigma}_{ii}-\langle\psi_{i}\rvert\sigma\lvert\psi_{i}\rangle\|+\\|\sigma\\|_{1}\cdot\|\tilde{\lambda}_{i}-{\lambda}_{i}\|$	$\displaystyle(\text{since }\\|\sigma\\|\leq\\|\sigma\\|_{1})$
	$\displaystyle\leq\varepsilon+\xi,$

Improved Quantum Algorithms for Fidelity Estimation

Abstract

1 Introduction

1.1 Related work

1.2 Our results

Fidelity estimation with block-encoding based algorithms.

Fidelity estimation via quantum spectral sampling.

1.3 𝖰𝖲𝖹𝖪𝖧𝖵\mathsf{QSZK_{HV}}-hardness of fidelity estimation to any non-trivial accuracy

Theorem 1 (𝖰𝖲𝖹𝖪𝖧𝖵\mathsf{QSZK_{HV}}-hardness of non-trivial fidelity estimation).

Proof.

1.4 A sample complexity lower bound for constant precision fidelity estimation

Corollary 2.

Proof.

Proposition 3.

Proof.

2 Preliminaries

Definition 4 (Purified access).

Lemma 5 (Projection onto the positive semidefinite cone).

2.1 Matrix Arithmetics using blocks of unitaries

Definition 6 (Block-encoding).

Lemma 7 (Block-encoding of density operators with purified access [GSLW18, Lemma 45]).

Lemma 8 (Product of block-encoded matrices [GSLW18, Lemma 53]).

2.2 The Swap Test

2.3 The Hadamard test and block-measurements

Lemma 9 (Hadamard test for estimating the expectation value of block-encoded matrices).

Proof.

2.4 Quantum Singular Value Transformation

Theorem 10 (Quantum Singular Value Transformation [GSLW18, Corollary 18 & Lemma 19]).

Lemma 11 (Robustness of singular value transformation [GSLW18, Lemma 23]).

Corollary 12.

Corollary 13 (Error propagation in block-encodings).

Proof.

2.5 Low-degree polynomial approximations

Lemma 14 (Polynomial approximations of the sign function).

Corollary 15 (Polynomial approximations of rectangle functions [GSLW19, Corollary 16]).

Lemma 16 (Polynomial approximations of negative power functions [Gil19, Cor. 3.4.13]).

Lemma 17 (Polynomial approximations of positive power functions [Gil19, Cor. 3.4.14]).

Corollary 18 (Polynomial approximations of positive power functions).

Proof.

2.6 Density matrix exponentiation and block-encodings

Theorem 19 (Density matrix exponentiation [KLL+17, Theorem 5 & Theorem 20]).

Lemma 20 (Implementing the logarithm of unitaries [GSLW18, Corollary 71]).

Corollary 21 (Sampling to block-encoding).

Proof.

2.7 Truncated fidelity bounds

Lemma 22 (Monotonicity of fidelity).

Proof.

Lemma 23 (Hard truncation bounds).

Proof.

Corollary 24 (Soft truncation bounds).

Proof.

3 Fidelity estimation with block-encoding based algorithms

3.1 Polynomial approximations and error bounds

3.2 Complexity analysis – with purified access

Theorem 25.

Proof.

Corollary 26.

Proof.

3.3 Complexity analysis – with sampling access

Theorem 27.

Proof.

Corollary 28.

Proof.

4 Fidelity estimation via spectral sampling

4.1 Continuity bound for fidelity estimation

Theorem 29 (Continuity bound).

Proof.

4.2 Algorithm

4.3 Quantum spectral sampling

Definition 30 (Infinitesimal SWAP operation).

Theorem 31 ([Pra14]).

4.4 Quantum eigenstate filtering

Theorem 32 ([LT20]).

4.5 Technical lemmas

Lemma 33.

Proof.

Lemma 34 (Estimation errors for diagonal terms).

Proof.

Lemma 35 (Estimation errors for off-diagonal terms).

Proof.

1.3 $\mathsf{QSZK_{HV}}$ -hardness of fidelity estimation to any non-trivial accuracy

Theorem 1 ( $\mathsf{QSZK_{HV}}$ -hardness of non-trivial fidelity estimation).

Theorem 19 (Density matrix exponentiation [KLL⁺17, Theorem 5 & Theorem 20]).