\section{Overlap estimation via a quasi-probabilistic representation}

A $ n $-qubit factorable, informationally complete positive-operator-valued measure (IC-POVM) is a set of operators $ \{M_{a}^{(n)}=M_{a_1}\otimes\ldots\otimes M_{a_n}\}_a $.
This means the POVM elements are positive semi-definite operators $ M_a^{(n)}\geq0 $ from $ \mathcal{L}(\mathbb{C}^2\otimes\ldots\otimes\mathbb{C}^2) $ that satisfy $ \sum_a M_a^{(n)} = \mathbb{I} $ and span the space of Hermitian operators.
This POVM defines an associated matrix $T = [\text{\normalfont Tr}(M_a^{(n)} M_b^{(n)})]_{ab}$.
Information-completeness allows us to write any $ n $-qubit state $ \rho $ as\cite{Carrasquilla_2019}
\begin{equation}\label{quasi_repr}
\rho = \sum_{a, a'} P_{\rho}(a)\tau_{aa'}M_{a'}^{(n)},
\end{equation}
where $ P_\rho(a) = \text{\normalfont Tr}(\rho M_a^{(n)}) $ and $ \tau$ can be defined as the pseudoinverse of $ T $ (however, as we show below, other choices of $ \tau $ are possible).
Since the POVM elements are not orthogonal, the associated probabilities do not refer to exclusive events, and therefore $ P_\rho $ is called a quasiprobability.

Let us now proceed to characterise all choices of $ \tau $.
Noticing that decomposition (\ref{quasi_repr}) holds for any operator, we can apply it to $ \rho = M_b^{(n)} $ and obtain
\begin{align}\nonumber
[T]_{bb'} =& \text{\normalfont Tr}(M_{b}^{(n)}M_{b'}^{(n)}) \\ =& \sum_{a, a'}\text{\normalfont Tr}(M_b^{(n)}M_a^{(n)})\tau_{aa'}\text{\normalfont Tr}(M_{a'}^{(n)}M_{b'}^{(n)}),
\end{align}
for any $b, b'$.
This yields the relation
\begin{equation}\label{tau}
T = T\tau T.
\end{equation}
Any matrix satisfying (\ref{tau}) is called a generalised inverse for $ T $; in particular, the pseudoinverse fulfils this condition.
Analogously, again from the completeness of the POVM, it is straightforward to check that (\ref{tau}) conversely implies (\ref{quasi_repr}) (see Appendix \ref{app:lemma}).
This proves the following observation.

\begin{lemma}\label{lemma}
	Eq. (\ref{quasi_repr}) holds for any state $ \rho $ if and only if $ T=T\tau T $.
\end{lemma}

Notice that we can use different tensors $ \tau_1,\tau_2 $ for representing different quantum states.
Using this fact to represent states $ \rho $ and $ \sigma $, we find their overlap is given by 
\begin{align}\label{estimator}\nonumber
\text{\normalfont Tr}(\rho \sigma) =& \sum_{a,b}P_\rho(a)[\tau_1T\tau_2^t]_{ab}P_\sigma(b) \\
=& \underset{a\sim P_\rho, b\sim P_\sigma}{\mathbb{E}}[\tau_1T\tau_2^t]_{ab}.
\end{align}
This provides a practical way of estimating the overlap Tr$ (\rho\sigma) $.
First, we repeatedly measure $ M^{(n)} $ on $ \rho $ and $ \sigma $.
This can be done separately and at different times.
Then, according to the obtained outcomes $ a,b $, we approximate the above expectation value by an average of the matrix elements $ \hat\tau=[\tau_1T\tau_2^t]_{ab} $ (see Fig. \ref{fig:estimator}).


\begin{figure}
	\includegraphics[width=\linewidth]{./estimator.pdf}
	\caption{Tensor network representing our state-overlap estimator, based on a quasi-probabilistic representation for quantum states via an IC-POVM $ M $, its associated matrix $ T $ and an auxiliary tensor $ \tau $.}
	\label{fig:estimator}
\end{figure}

This approximation generates a statistical error, which we now consider.
Recall that representation (\ref{quasi_repr}) allows to reconstruct entangled states $ \rho $ as combinations of the separable operators $ M_a $, implying that $ \hat\tau$ forcibly has negative entries.
This negativity can be quantified in terms of the range $ \nu = \max_{a,b}[\hat\tau]_{ab} - \min_{a,b}[\hat\tau]_{ab} $.
Moreover, it controls the number of samples that our method requires in general, as stated in the following.

\begin{thm}
	Let $ \rho, \sigma $ be $ n $-qubit states and $\{\hat\tau_i\}_{i=1}^N $ be $ N $ realisations of $[\hat \tau]_{ab}$, with $ a\sim P_\rho , b\sim P_\sigma$. Then $| \frac1N \sum_i \hat\tau_i - \text{\normalfont Tr}(\rho\sigma)| < \epsilon $ with probability $ 1-\delta $ whenever $ N > \frac{\nu^n}{2\epsilon^2}\log\left(\frac{2}{\delta}\right) $.
\end{thm}

The proof follows directly from Hoeffding's inequality.
Hence, we see the sample complexity scales exponentially on the number of qubits according to $\nu$, which encodes the sample complexity overhead inherent to the quasiprobabilictic representation.

\section{Numerical experiments}

We performed numerical experiments to validate our estimator by taking $ n $ tensor copies of the 6-outcome POVM 
\begin{equation}\label{povm}
\{M_a = \frac13\ket{a}\bra{a}\}_a,
\end{equation}
where $ \ket{a}$ is an eigenvector of either $ X, Y,$ or $ Z $, and $\tau_1=\tau_2=T^+$, where $ T^+ $ denotes the pseudoinverse of $ T $.
Due to the freedom in the choice per equation Eq. (\ref{tau}), we considered different POVMs for our experiments.
We performed a brute-force search on the 4-, 6-, and 8-outcome qubit POVMs by incrementally varying the angles between Bloch-vectors representing the POVM elements.
This revealed the Pauli-6 POVM as the optimal choice.
We also used a Markov chain Monte Carlo\cite{Metropolis1953mh, Sherlock2010rwmh} dynamics to simulate a random walk on the set of solutions to Eq. (\ref{tau}) for minimising $ \nu $.
However, both searches provided only worse negativities.
This is further detailed in Ref. \cite{roeland}.

\begin{figure}[t]
	\includegraphics[width=\linewidth]{./scaling.pdf}
	\caption{Average number of samples required for the absolute mean error of the overlap estimation to be $ <0.05 $, considering 180 different pairs of either product or entangled pure random states, according to the number of qubits in the states.}
	\label{fig:scaling}
\end{figure}

We numerically investigate the number of samples $ N $ required to make the estimation's absolute error, $ |\text{\normalfont Tr}(\rho\sigma) - \frac1N\sum_i\hat\tau_i| $, lower than a fixed threshold.
We estimate the overlap of pure product and entangled pairs of states, constructed as matrix-product states from random local tensors.
In agreement with Ref. \cite{Elben2020}, we find that the latter is a slightly simpler problem.
In Fig. \ref{fig:scaling}, we present the obtained scaling for the fixed error bound of 0.05.
To minimise numerical fluctuations, we check the number of samples required to make the average error of 5 random pairs lower than this bound, for 100 batches of 5 states.
For more details, we refer to our computational appendix \cite{comp_app}.

\section{Connection with previous approaches}

We now investigate other overlap estimators, which we simulate in order to compare them with the quasiprobabilistic method. 
These alternative methods are separated in two groups.
We start from the ones that require coherently-controlled  operations between the two systems, described as quantum circuits.
Next, we consider two direct estimators based on random measurements and the theory of unitary $ t $-designs. 

\textbf{Entanglement-based approaches.---}
There are various quantum circuits designed specifically for calculating the overlap between quantum states.
However, considering realistic implementations, most NISQ devices have a series of restrictions in common.
For instance, the number of implementable quantum gates is limited, as is the connectivity between different qubits.
In addition, although the single-qubit gate noise is usually negligible, the noise present in multiqubit gate implementations cannot be ignored.
To address these points in our simulations of noisy circuits, we used the locally purified density operator (LPDO) paradigm, designed for simulating non-unitary evolutions in tensor networks language (see Appendix \ref{app:lpdo}).

Based on that, we benchmark our method considering quantum computers with the following features: (i) the set of implementable gates comprises noiseless single-qubit gates and noisy CNOTs; (ii) the circuits are one-dimensional; and (iii) the interactions between the qubits are limited to nearest-neighbours, meaning that long-distance CNOTs have to be compiled via first-neighbour CNOTs. 
The optimal compilation then adds $4(d-1)$ extra CNOTs to the original circuit, where $ d $ is the distance between the qubits in the long-distance CNOT.
We assume the noisy CNOTs are well modelled by an ideal CNOT followed by a depolarising channel.
This yields a circuit structure with increasing infidelity.
We emphasise the one-dimensional structure is considered only for simplicity and illustrative purposes.
Our conclusions hold also for 2D structures, but the corresponding classical simulation is more expensive.



Under these conditions, let us recall three circuits for state overlap.
The first is the standard SWAP test\cite{Buhrman2001}, which requires a single-qubit measurement on the ancillary qubit. 
Similarly, the improved SWAP test circuit\cite{Cincio2018} optimises the compilation of the standard circuit in view of constraint (i), demanding a lower number of layers and CNOTs but maintaining the single-qubit measurement on the ancilla.
Finally, on a different line, the Bell-basis circuit\cite{Garcia2013} requires a global measurement on the Bell basis in exchange of reducing the number of gates and layers.
With the limited connectivity assumption (iii), these circuits already involve a high number of CNOT gates for a small number of qubits (see Table \ref{tablee}).

\begin{table}[h]
	\begin{tabular}{c|c|c|c|c|c|}
		\cline{2-6}
		\textbf{}                           & 1Q & 2Q & 3Q  & ... & 8Q   \\ \hline
		\multicolumn{1}{|c|}{Improved SWAP test} & 12 & 60 & 144 & ... & 1104 \\ \hline
		\multicolumn{1}{|c|}{Bell-basis circuit}    & 1  & 10 & 27  & ... & 232  \\ \hline
	\end{tabular}
	\caption{Number of nearest-neighbours CNOTs required from the improved SWAP test and Bell-basis circuits, according to the number of qubits in the states.}
	\label{tablee}
\end{table}



In these approaches, on top of the statistical error of estimating the expectation value of an observable, we have the systemic noise arising from the implementation of two-qubit gates.
In contrast, using the quasi-probabilistic direct approach we obtain a qualitative advantage since the estimation (\ref{estimator}) can be as accurate as needed, as long as the number of samples is large enough.
As shown in Fig. \ref{fig:circuits}, even for a modest number of samples we already see a clear advantage to our method in the cases of 2- and 3-qubit states.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.9\linewidth]{./histograms2.pdf}
	\caption{Histograms of the absolute error in the overlap estimation of 180 pairs of states by the quasi-probabilistic and circuit-based methods. In both circuits we consider a CNOT-noise of 0.005 (single-qubit gates are implemented noiselessly) and measurement noise of 0.01. (a) The improved SWAP-test. The comparison is done based on 500 samples from two-qubit random states (with an averaged overlap of 0.86), for which the near-neighbours circuit has 54 layers. (b) The Bell-basis test. We use 1,500 samples from 3-qubit random states (with an averaged overlap of 0.78), for which the near-neighbours circuit has 28 layers.}
	\label{fig:circuits}
\end{figure}

\textbf{Random unitary-based methods.---}
We now compare our estimator with the two state-of-the-art direct techniques for calculating state overlaps.
The first technique\cite{Elben2020} applies the same random measurements on both states, coordinated by means of classical communication.
Nevertheless, it is also entanglement-free and allows for independent interaction with each system.
Next, it uses the theory of unitary $ t $-designs for  approximating the overlap.
The direct performance comparison with our method (Fig. \ref{fig:inn}) shows that our approach presents a lower mean absolute error for $ n $-qubit states, with $ n<7 $.
However, our method is outperformed for states involving more qubits.
This is confirmed by the sample complexity required to obtain an error bound of $ 0.05 $, $ \sim 2^{1.1n}, $ shown in Fig. \ref{fig:scaling}, while the reported corresponding scaling of the random measurements method is $ \sim 2^{0.8n} $.

\begin{figure}[t]
	\includegraphics[width=\linewidth]{./innsbruck.pdf}
	\caption{Comparison between the quasi-probabilistic and Elben et al. estimator\cite{Elben2020}, for a fixed budget of 10,000 samples, for random pure product and pure entangled states. In the latter, we used 100 samples from each of 100 random measurements. The mean absolute error was computed over 200 instances and the average overlap was $\approx0.6$ in both cases.}
	\label{fig:inn}
\end{figure}


Finally, we discuss the connection of our results with shadow estimation\cite{Huang_2020, kliesch2021}, a recently proposed method to efficiently estimate linear functions of quantum states.
The so-called classical shadow of a state $ \rho $, given by the inverse of the map $ \mathcal{M}(\rho) = \mathbb{E}_{U, b}(U\ket b\bra b U^\dagger) $, is obtained by taking the average over some ensemble of unitaries $ U $ to $ \rho $ and subsequently measuring on the computational basis $ \{\ket b\bra b\}_b $.
Whenever the unitaries are chosen from the local Clifford group, the corresponding shadow fidelity estimator matches Eq. (\ref{estimator}), for $ M_a^{(1)}$ as in Eq. (\ref{povm}) and for a particular choice of $ \tau_1, \tau_2 $ (see Appendix \ref{app:shadow}).
Therefore, these specific instances of each estimator share the same performance\cite{Huang_2020}


\section{Final discussion}

While efficient entanglement-based quantum algorithms for approximating the overlap of quantum states exist, the currently available devices imply unavoidable systematic errors.
Therefore, direct and entanglement-free methods may provide an advantage in the estimation of state overlaps in the near term.
These typically present only statistical errors and treat the overlap estimation as a sample-complexity problem.
Among the entanglement-free direct estimators, we show that our method presents the lowest mean absolute error for n-qubits states, with n < 7.
We also identify that the quasiprobabilistic construction shares some common elements with the shadow fidelity estimation approach.
Moreover, for particular choices on each of the methods, the constructed estimators are the same.

Our work opens up avenues for future explorations.
In here, we have used a quasiprobabilistic representation
specifically for estimating the fidelity between quantum
states.
However, these representations can be used to describe not only quantum states, but also quantum measurements, unitary evolution, and open system quantum dynamics\cite{Pashayan2015, Carrasquilla_21, Luo_21, Reh_21}.
Can they be efficiently used for classical simulations of quantum circuits?
For what circuits the associated sample-complexity overhead becomes too expensive?
It would be interesting to investigate the interplay between these figures in the context of classical-quantum simulations.


\textit{Acknowledgements} ---
The authors thank Ingo Roth for helpful discussions.
LG and LA were financially supported by the Serrapilheira Institute (grant number Serra-1709- 17173) and the Brazilian agencies CNPq (PQ grant No. 311416/2015-2 and INCT-IQ), FAPERJ (PDR10 E- 26/202.802/2016 and JCN E-26/202.701/2018) and S\~ao Paulo Research Foundation (FAPESP) under grants 2016/01343-7 and 2018/04208-9.
LA also acknowledges support from the Technology Innovation Institute of Abu Dhabi, United Arab Emirates.
JC acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC), the Shared Hierarchical Academic Research Computing Network (SHARCNET), Compute Canada, Google Quantum Research Award, and the CIFAR AI chair program.
Resources used in preparing this research were provided, in part, by the Province of Ontario, the Government of Canada through CIFAR, and companies sponsoring the Vector Institute www.vectorinstitute.ai/\#partners.


