Introduction

In recent years, great advances have been achieved in the development of technologies based on quantum systems. These quantum technologies aim to provide significant performance improvements compared to their classical counterparts. A contemporary example is quantum computing, where large numbers of two-dimensional quantum systems (qubits) are controlled via canonical sets of operations to solve certain computational tasks more efficiently than classical computers1,2. Several quantum systems, such as trapped ions3,4,5,6, superconducting circuits7,8,9,10, quantum dots11,12,13,14, and photonic platforms15,16,17,18, among others, have been proposed and tested to encode qubits and implement quantum gates, in order to build quantum simulators and quantum computers. Today, private companies provide cloud-based access to their quantum hardware and software.

Currently, the above experimental platforms, among others, have evolved into devices that have been described as noisy intermediate-scale quantum (NISQ) technologies19. These are devices that operate with imperfect gates on a number of qubits between 50 and up to a few hundred qubits. These devices are beyond the simulation capabilities of current classical computing devices but are still far from exhibiting a clear enough improvement in useful computing power. In particular, noisy quantum gates and decoherence limit the depth of the circuits and prevent the realization of complex algorithms20,21,22.

Because of the properties of NISQ computers, the design of methods for the characterization, certification, and benchmarking of these quantum devices has become increasingly relevant and difficult23. For instance, the characterization of states generated by this class of devices requires measurements on a large set of bases. These are physically implemented by applying sets of universal unitary transformations on the qubits followed by measurements on a computational basis. Typically, an effort to reduce the number of bases leads to a new set of bases that require the application of entangling unitary transformations, such as, for instance, the controlled NOT (CNOT) gate, whose implementation is associated with a large error. Consequently, the estimation of states generated by NISQ devices could be affected by a large build-up of inaccuracy in the measurement procedure.

Several methods have been proposed to assess the performance of NISQ devices. Randomized benchmarking24,25,26, direct fidelity estimation27,28, and entanglement certification29,30 are popular tools for this purpose. However, they do not give complete information about the device in question. In contrast, quantum state estimation31 gives plenty of information about the device since it fully determines the state of the system from a set of suitably chosen measurements.

Today, several protocols to perform quantum state estimation of d-dimensional quantum systems are available. These attempt to estimate the d2 − 1 real parameter that characterizes density operators. Standard quantum tomography is based upon the measurement of the d2 − 1 generalized Gell-Mann matrices32,33,34,35. Quantum state estimation based on projections onto the states of d + 1 mutually unbiased bases has been suggested36,37,38,39,40,41,42 to reduce the total number of measurement outcomes. This number can be reduced to a minimum using a symmetric, informationally complete, positive operator-valued measure43,44,45,46,47,48,49, which has exactly d2 measurement outcomes.

Despite the progress made in the reduction of the total number of measurement outcomes used by quantum state estimation methods, the exponential increase of this number prevails in the case of multipartite systems. This is known as the curse of dimensionality. For example, standard quantum tomography for a n-qubit system requires the measurement of 3n local Pauli settings. A solution to overcome this problem is the use of a priori information about the state to be determined. Compressed sensing50 allows a large reduction in the total number of measurements as long as the rank of the state is known. Another alternative is to consider states with special properties, such as, for instance, matrix product states51,52 or permutationally invariant states53, which also leads to a significant reduction in the total number of measurements. The important case of pure states has also been studied. In this case a set of five bases allows one to estimate with high accuracy all pure quantum states54,55,56. Recently, a 3-bases-based tomographic scheme has been introduced to estimate a single qudit57. The generalization of this method to multiple qubits leads to entangled bases. Adaptivity has also been introduced as a means of increasing the estimation accuracy of tomographic methods58,59,60,61,62 and even reaching fundamental accuracy limits63,64.

Here, we present a method to estimate unknown pure quantum states in n-qubit Hilbert spaces, as used by NISQ computers, by means of mn + 1 local bases (with m an integer number), which correspond to the tensor product of n single-qubit bases, or by the computational basis plus m entangled bases, that is, bases that cannot be decomposed as the tensor product of single-qubit bases. Thereby, the number of bases scales linearly as O(mn) and O(m), respectively, with the number n of qubits. This is a great advantage over standard methods, which are of exponential order. The mn + 1 local bases also provide an advantage in the estimation of states generated via NISQ devices, such as current prototypes of quantum processors, because in this case projective measurements can be carried out without the use of typically noisy entangling gates. The computational basis together with m = 2 entangled bases also leads to a clear reduction with respect to the 5-bases tomographic method54,55. The method here proposed also improves over the 3-bases method, since this requires calculating the likelihood function of \({2}^{{2}^{n}-1}\) states to obtain an estimate. For a large number of qubits, this stage exponentially increases the computational cost of the method. Furthermore, in order to obtain distinguishable values of the likelihood, a large ensemble size is required. Our estimation method does not require such a procedure.

The present method estimates pure multi-qubit quantum states using an inductive process. For an arbitrary n-qubit state \(\left|{{\Psi }}\right\rangle\) we define n sets of 2n−j sub-states of dimension 2j, with j = 1, 2, . . . , n and with the sub-state of dimension 2n equal to \(\left|{{\Psi }}\right\rangle\) except a global phase. We first show how to estimate the sub-states of dimension 2 and then show that the knowledge about sub-states of dimension 2j−1 together with some measurement outcomes allows one to estimate the sub-states of dimension 2j. This points to an inductive estimation method: we first estimate the sub-states of dimension 2, and in n iterations we arrive to the sub-state of dimension 2n, completing the estimation procedure. We show that the set of sub-states can be estimated by projecting the state onto mn + 1 local bases or m entangled bases plus the computational basis. Throughout the proposed estimation method, it is necessary to solve several systems of linear equations that might have vanishing determinants if the number of measurements is insufficient, leading to the failure of the method. However, the presence of noise, such as finite statistics, helps to mitigate this problem. We also show that a small increase in the number of local or entangled bases solves this problem and improves the overall fidelity of the estimation. Our method is specifically formulated to estimate pure muti-qubit states. However, the proposed method can also estimate pure states that have been transformed into mixed states by the action of white noise. In particular, this is achieved without modifying the measurement bases or their number. The method provides an estimate of the noise parameter.

We study the present method through Monte Carlo numerical simulations. We randomly generate, according to a Haar-uniform distribution, a set of unknown pure states and calculate the average estimation fidelity as a function of the number of qubits. Projective measurements on each basis are simulated considering a fixed number of repetitions. We show that 2n + 1 separable bases achieve an average estimation fidelity of 0.88 for n = 10. We improve this figure by considering a larger number of bases, that is, 3n + 1 and 4n + 1, which lead to an average estimation fidelity of 0.915 and 0.93 for n = 10, respectively. If the average estimation fidelity is calculated with respect to the set of separable states, then for the above number of bases we obtain the values 0.95, 0.955, and 0.96, correspondingly. The use of entangled bases leads to a less favorable picture, where 2, 3, and 4 entangled bases plus the computational basis lead to an average estimation fidelity of 0.2, 0.6, and 0.8 for n = 6, respectively. These figures increase when estimating separable states, where we achieve an average estimation fidelity of 0.75, 0.95, and 0.95 for 2, 3, and 4 entangled bases plus the computational base, respectively. The large gap between separable and entangled bases can be explained by considering that measurements are simulated with a fixed number of 213 repetitions or, equivalently, with a fixed ensemble size. This is the number of repetitions currently available at IBM’s quantum processors. Thereby, the total ensemble size is much smaller in the case of the entangled bases, which decreases the estimation fidelity. We also conduct simulations with 223 repetitions, which exhibit a large increase in fidelity. For instance, for n = 10 qubits 2n + 1 local bases achieve an average estimation fidelity of 0.999. Our simulations also indicate that the average computation time required to estimate an arbitrary 20-qubit state is of the order of 3 min, when using a portable personal computer.

We test the proposed method by means of experiments carried out using the quantum processor ibmq_montreal provided by IBM. We consider the estimation of fixed quantum states from 2 to 10 qubits using the local bases, and states up to 4 qubits using the entangled bases. We analyze our results using the theoretical fidelity, which is the fidelity between the state to be ideally prepared and its estimate, the experimental fidelity, which is the fidelity between the actually prepared state and its estimate, and which is obtained via direct fidelity estimation, and a purity estimate. We first estimate two completely separable states. One of these states is such that the matrices to be inverted are ill-conditioned. Nevertheless, the estimation of these states via local bases leads to very similar theoretical fidelities above 0.875 for n = 10. The estimation via entangled bases leads to the ill-conditioned state to a theoretical fidelity in the order of 0.15 for n = 4 in the interval [0.7, 0.8] for n = 4. The other state exhibits theoretical fidelities in the interval [0.7, 0.8] for n = 4. We also consider the estimation of an n/2-fold tensor product of a Bell state, where local bases for n = 10 lead to a theoretical fidelity in the interval [0.7, 0.8] and entangled bases lead to a theoretical fidelity above 0.7 for n = 10 and m = 3, 4. Finally, we test the estimation of a Greenberger-Horne-Zeilinger state for 2, 3, and 4 qubits. The proposed method leads to a theoretical fidelity of ~0.87 for n = 4 via local bases while in the case of entangled bases the theoretical fidelity is in the interval [0.6, 0.7]. The large decrease in the theoretical fidelity of entangled states by means of entangled bases can be explained by the use of low-accuracy entangling gates in the preparation of the states as well in the projection of the states of the entangled bases. The experimental fidelity follows the general trend of the theoretical fidelity albeit with lower fidelities and the purity estimate is above 0.8.

Our simulations and experimental results indicate that the method proposed here allows the efficient and scalable estimation of n-qubit pure states by means of mn + 1 separable bases in NISQ computers, while the m entangled bases together with the computational base can offer a large advantage in quantum computers based on high-accuracy entangling gates.

Results

Theory

Let us consider a n-qubit system described by the pure state

$$\left|{{\Psi }}\right\rangle =\mathop{\sum }\limits_{\alpha =0}^{{2}^{n}-1}{c}_{\alpha }{e}^{i{\phi }_{\alpha }}{\left|\alpha \right\rangle }_{n},$$
(1)

where cα ≥ 0 and \({\left|\alpha \right\rangle }_{n}=\left|{\alpha }_{n-1}\right\rangle \otimes \cdots \otimes \left|{\alpha }_{0}\right\rangle\) is the n-qubit computational basis, with \(\alpha =\mathop{\sum }\nolimits_{k = 0}^{n-1}{2}^{k}{\alpha }_{k}\) the integer associated with the n-bit binary string αn−1 ⋯ α0.

Our main aim is to estimate the values of the amplitudes {cα} and the phases {ϕα} with a total number of measurements that do not scale exponentially with the number n of qubits. For this purpose, we employ an iterative algorithm based on estimating sub-states, which are non-normalized j-qubit vectors defined by

$$\left|{{{\Psi }}}_{\beta }^{j}\right\rangle ={\mathop{\sum }\limits_{\alpha = 0}^{{2}^{j}-1}}_{n}\left\langle {2}^{n-j}\beta +\alpha | {{\Psi }}\right\rangle {\left|\alpha \right\rangle }_{j},$$
(2)

with 1 ≤ j ≤ n and 0 ≤ β ≤ 2n−j − 1. Thus, for a given state \(\left|{{\Psi }}\right\rangle\) we have n sets of 2n−j sub-states, one set for each j and one sub-state for each β. Every sub-state has some partial information about the full state. For instance, if j = 1 we have 2n−1 sub-states that are given by

$$\left|{{{\Psi }}}_{\beta }^{1}\right\rangle ={c}_{2\beta }{e}^{i{\phi }_{2\beta }}\left|0\right\rangle +{c}_{2\beta +1}{e}^{i{\phi }_{2\beta +1}}\left|1\right\rangle .$$
(3)

We can see that the amplitudes and the phases are the same enterings the state \(\left|{{\Psi }}\right\rangle\). Thereby, we can reduce the problem of estimating the unknown state \(\left|{{\Psi }}\right\rangle\) to the problem of estimating the set of sub-states. The main obstacle of this approach is that measurements in quantum mechanics do not contain information about global phases, and hence if we try to find every sub-state in an independent way using the results of some measurements, we can only obtain the sub-states \(|{{{\Psi }}}_{\beta }^{j}\rangle\) except a global phase, which in this case, would be a relative phase in \(\left|{{\Psi }}\right\rangle\). Thus, the only sub-state that fully characterizes the system is \(\left|{{{\Psi }}}_{0}^{n}\right\rangle\), which is the same as \(\left|{{\Psi }}\right\rangle\). Nevertheless, the remaining sub-states will also be useful in the task of reconstructing \(\left|{{\Psi }}\right\rangle\), because we can find the sub-states \(|{{{\Psi }}}_{\beta }^{j}\rangle\) (except a global phase) using the sub-states \(|{{{\Psi }}}_{\beta }^{j-1}\rangle\), as we will show.

For a large enough set of identical copies of the state \(\left|{{\Psi }}\right\rangle\), measurements in the computational base \(\{\left|\alpha \right\rangle \}\) gives us a histogram of observations such that we have approximations pα to ∣〈α∣Ψ〉∣2. Then the amplitudes can be estimated as \({c}_{\alpha }=\sqrt{{p}_{\alpha }}\). On the other hand, the phases can be determined from projective measurements on the following states

$$\left|{P}_{a\beta }^{j}\right\rangle ={\left|\beta \right\rangle }_{n-j}\otimes \left|{+}_{a}\right\rangle \otimes {\left|{-}_{a}\right\rangle }^{\otimes j-1},$$
(4)

where \(\left|{+}_{a}\right\rangle ={u}_{a}\left|0\right\rangle +{v}_{a}{e}^{i{\varphi }_{a}}\left|1\right\rangle\) and \(\left|{-}_{a}\right\rangle ={v}_{a}\left|0\right\rangle -{u}_{a}{e}^{i{\varphi }_{a}}\left|1\right\rangle\) are orthonormal single-qubit states, with a = 1, …, m. Here, m is the number of bases \(\{\left|{\pm }_{a}\right\rangle \}\) considered, which must be large enough to carry out the algorithm. In addition, these bases have to be all different from the computational ones.

To start the algorithm, we set j = 1 and the sub-states are given by Eq. (3). The probabilities of projecting \(\left|{{\Psi }}\right\rangle\) onto \(|{P}_{a\beta }^{1}\rangle\) are given by

$$\begin{array}{l}{P}_{a\beta }^{1}={u}_{a}^{2}{c}_{2\beta }^{2}+{v}_{a}^{2}{c}_{2\beta +1}^{2}\\\qquad\;\; +\,2{u}_{a}{v}_{a}{c}_{2\beta }{c}_{2\beta +1}\left[\cos ({\varphi }_{a})\cos ({\phi }_{2\beta +1}-{\phi }_{2\beta })\right.\\\qquad\;\; \,+\,\left.\sin ({\varphi }_{a})\sin ({\phi }_{2\beta +1}-{\phi }_{2\beta })\right].\end{array}$$
(5)

These generate a set of equations, one equation for each value of β, that are linear combinations of the cosine and sine of the relative phase δϕβ = ϕ2β+1 − ϕ2β, with coefficients depending on \(\left|{+}_{a}\right\rangle\). Thereby, we can form a linear system of equations Lβaβ = bβ for trigonometric functions of the relative phase, which is explicitly given by

$$\left[\begin{array}{cc}\cos {\varphi }_{1}&\sin {\varphi }_{1}\\ \vdots &\vdots \\ \cos {\varphi }_{m}&\sin {\varphi }_{m}\end{array}\right]\left[\begin{array}{c}\cos \delta {\phi }_{\beta }\\ \sin \delta {\phi }_{\beta }\end{array}\right]=\left[\begin{array}{c}{\tilde{P}}_{1\beta }^{1}\\ \vdots \\ {\tilde{P}}_{m\beta }^{1}\end{array}\right],$$
(6)

where we have defined

$${\tilde{P}}_{a\beta }^{1}=\frac{1}{2{u}_{a}{v}_{a}{c}_{2\beta }{c}_{2\beta +1}}({P}_{a\beta }^{1}-{u}_{a}^{2}{c}_{2\beta }^{2}-{v}_{a}^{2}{c}_{2\beta +1}^{2}).$$
(7)

Since we know the projectors \(|{P}_{a\beta }^{1}\rangle\) and the coefficients cj, this system can be solved for m ≥ 2 inverting Lβ by the Moore-Penrose pseudo-inverse A+. When the matrix A has linearly independent columns, the Moore-Penrose pseudo-inverse can be calculated as \({A}^{+}={({A}^{T}A)}^{-1}{A}^{T}\). We can guarantee that the pseudo-inversion is carried out with good quality by suitably choosing the coefficients ua, va, and the phases φa entering in the states \(|{P}_{a\beta }^{1}\rangle\) such that Lβ is a well-conditioned matrix. This can be verified by the condition number \({{{\rm{cond}}}}(A)={\sigma }_{\max }(A)/{\sigma }_{\min }(A)\), where \({\sigma }_{\max }(A)\) and \({\sigma }_{\min }(A)\) are the maximum and minimum singular values of A, respectively. If the condition number of A is much larger than 1, the matrix is ill-conditioned, the pseudo-inversion has poor quality, and a low fidelity estimation is achieved. For example, taking the particular case of m = 2 the system has a solution

$$\left[\begin{array}{c}\cos \delta {\phi }_{\beta }\\ \sin \delta {\phi }_{\beta }\end{array}\right]=\frac{1}{\det ({L}_{\beta })}\left[\begin{array}{cc}\sin {\varphi }_{2}&-\sin {\varphi }_{1}\\ -\cos {\varphi }_{2}&\cos {\varphi }_{1}\end{array}\right]\left[\begin{array}{c}{\tilde{P}}_{1\beta }^{1}\\ {\tilde{P}}_{2\beta }^{1}\end{array}\right]$$
(8)

as long as the determinant of Lβ does not vanish, that is, det\(({L}_{\beta })=\cos ({\varphi }_{1})\sin ({\varphi }_{2})-\sin ({\varphi }_{1})\cos ({\varphi }_{2})\,\ne\, 0\). Thereby, taking φ1 = 0 and φ2 = π/2, or equivalently \(\left|{\pm }_{1}\right\rangle ={u}_{1}\left|0\right\rangle \pm {v}_{1}\left|1\right\rangle\), \(\left|{\pm }_{2}\right\rangle ={u}_{2}\left|0\right\rangle \pm i{v}_{2}\left|1\right\rangle\), the system of equations can be always solved since det(Lβ) = cond(Lβ) = 1. The exponential of the relative phase is given by \({e}^{i\delta {\phi }_{\beta }}={\tilde{P}}_{1\beta }^{1}+i{\tilde{P}}_{2\beta }^{1}\). Therefore, the sub-states \(\left|{{{\Psi }}}_{\beta }^{1}\right\rangle\) can be determined up to a global phase,

$$\left|{\tilde{{{\Psi }}}}_{\beta }^{1}\right\rangle ={c}_{2\beta }\left|0\right\rangle +{c}_{2\beta +1}{e}^{i\delta {\phi }_{\beta }}\left|1\right\rangle,$$
(9)

provided we measure the computational basis and at least two projectors \(\left|{P}_{a\beta }^{1}\right\rangle\) for every β in 0, 1, . . . , 2n−1 − 1. Nevertheless, when some of the computational basis probabilities pα are null, we have to consider a particular rule to obtain the sub-state. If any or both of the coefficients \({c}_{2\beta }=\sqrt{{p}_{2\beta +1}}\) or \({c}_{2\beta +1}=\sqrt{{p}_{2\beta +1}}\) are null, the sub-state is determined only by the computational basis, setting \(\left|{\tilde{{{\Psi }}}}_{\beta }^{1}\right\rangle ={c}_{2\beta }\left|0\right\rangle\), \(\left|{\tilde{{{\Psi }}}}_{\beta }^{1}\right\rangle ={c}_{2\beta +1}\left|1\right\rangle\) or \(\left|{\tilde{{{\Psi }}}}_{\beta }^{1}\right\rangle ={{{\bf{0}}}}\) in Eq. (9). This means that there may be null sub-states, which will have an impact on the algorithm in posterior iterations.

In the following iterations, that is, the case j > 1, the sub-states can be expressed as linear combinations of the previous ones,

$$\left|{\tilde{{{\Psi }}}}_{\beta }^{j}\right\rangle =\left|0\right\rangle \otimes \left|{\tilde{{{\Psi }}}}_{2\beta }^{j-1}\right\rangle +{e}^{i\delta {\phi }_{\beta }^{j}}\left|1\right\rangle \otimes \left|{\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\right\rangle,$$
(10)

with \(\left|{\tilde{{{\Psi }}}}_{\beta }^{j}\right\rangle ={e}^{-i{\phi }_{{2}^{j}\beta }}\left|{{{\Psi }}}_{\beta }^{j}\right\rangle\) the corresponding sub-states up to the global phase and \(\delta {\phi }_{\beta }^{j}={\phi }_{{2}^{j}(\beta +1/2)}-{\phi }_{{2}^{j}\beta }\) relative phases. Thus, assuming that we know the sub-states of the previous iteration (except for a global phase), we can determine the sub-state of the current iteration (except for a global phase) simply by determining the relative phase. Analogously to the first iteration, if any or both of the previous sub-states \(\left|{\tilde{{{\Psi }}}}_{2\beta }^{j-1}\right\rangle\) or \(\left|{\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\right\rangle\) are null, there is no relative phase to determine and the next sub-state is simply obtained with the choice δϕβ = 0. Otherwise, we have to determine the relative phase from the projections \(\left|{P}_{\beta }^{j}\right\rangle\). Considering j > 1, the probability of projecting the state \(\left|{{\Psi }}\right\rangle\) onto the state \(\left|{P}_{a\beta }^{j}\right\rangle\) is given by

$$\begin{array}{l}{P}_{a\beta }^{j}={u}_{a}^{2}\left| \left\langle {W}_{a}^{j}\left| {\tilde{{{\Psi }}}}_{2\beta }^{j-1}\right.\right\rangle \right| ^{2}+{v}_{a}^{2}\left| \left\langle {W}_{a}^{j}\left| {\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\right.\right\rangle \right| ^{2}\\\qquad\;\; +\,2{u}_{a}{v}_{a}\Re \left[{e}^{i(\delta {\phi }_{\beta }^{j}-{\varphi }_{a})}\left\langle {\tilde{{{\Psi }}}}_{2\beta }^{j-1}\left| {W}_{a}^{j}\right.\right\rangle \left\langle {W}_{a}^{j}\left| {\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\right.\right\rangle \right],\end{array}$$
(11)

where \(\left|{W}_{a}^{j}\right\rangle ={\left|{-}_{a}\right\rangle }^{\otimes j-1}\). Defining the quantities

$${X}_{a\beta }^{j}={e}^{-i{\varphi }_{a}}\langle {\tilde{{{\Psi }}}}_{2\beta }^{j-1}| {W}_{a}^{j}\rangle \langle {W}_{a}^{j}| {\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\rangle$$
(12)

and

$${\tilde{P}}_{a\beta }^{j}=\frac{1}{2{u}_{a}{v}_{a}}\left({P}_{a\beta }^{j}-{u}_{a}^{2}| \langle {W}_{a}^{j}| {\tilde{{{\Psi }}}}_{2\beta }^{j-1}\rangle {| }^{2}-{v}_{a}^{2}| \langle {W}_{a}^{j}| {\tilde{{{\Psi }}}}_{2\beta +1}^{j-1}\rangle {| }^{2}\right),$$
(13)

we obtain the following system of equations for the relative phases

$$\left[\begin{array}{cc}\Re {X}_{1\beta }^{j}&-\Im {X}_{1\beta }^{j}\\ \vdots &\vdots \\ \Re {X}_{m\beta }^{j}&-\Im {X}_{m\beta }^{j}\end{array}\right]\left[\begin{array}{c}\cos \delta {\phi }_{\beta }^{j}\\ \sin \delta {\phi }_{\beta }^{j}\end{array}\right]=\left[\begin{array}{c}{\tilde{P}}_{1\beta }^{j}\\ \vdots \\ {\tilde{P}}_{m\beta }^{j}\end{array}\right],$$
(14)

where \(\Re\) and \(\Im\) stand for a real and imaginary part, respectively. Again, the phases can be obtained by solving a linear system of equations \({L}_{\beta }^{j}{{{{\bf{a}}}}}_{\beta }^{j}={{{{\bf{b}}}}}_{\beta }^{j}\) given by Eq. (14) through the Moore-Penrose pseudo-inverse. If for a fixed j, we are able to pseudo-invert the matrices \({L}_{\beta }^{j}\), we can find all the \(|{\tilde{{{\Psi }}}}_{\beta }^{j}\rangle\). Then, we continue inductively with j + 1 until reaching j = n. In that case, the protocol ends because the sub-state \(\left|{\tilde{{{\Psi }}}}_{0}^{n}\right\rangle\) is equal to the state \(\left|{{\Psi }}\right\rangle\) up to a global phase.

One might think that, analogously to the case j = 1, the matrix \({L}_{\beta }^{j}\) can be inverted for m≥2. Thereby, to completely estimate the unknown state \(\left|{{\Psi }}\right\rangle\) we would need a minimum of 2(2n − 1) different projections onto states \(|{P}_{a\beta }^{j}\rangle\) in (4), because there are 2n − 1 different sub-states. However, now the matrix \({L}_{\beta }^{j}\) to be pseudo-inverted does not only depend on the projectors \(|{P}_{a\beta }^{j}\rangle\), but it also depends on the sub-states of the previous iteration, so that we cannot guarantee that the pseudo-inversion is feasible. For example, taking m = 2, the solution of the system of equations (14) is

$$\left[\begin{array}{c}\cos \delta {\phi }_{\beta }^{j}\\ \sin \delta {\phi }_{\beta }^{j}\end{array}\right]=\frac{1}{\det ({L}_{\beta }^{j})}\left[\begin{array}{cc}-\Im {X}_{2\beta }^{j}&\Im {X}_{1\beta }^{j}\\ -\Re {X}_{2\beta }^{j}&\Re {X}_{1\beta }^{j}\end{array}\right]\left[\begin{array}{c}{\tilde{P}}_{1\beta }^{j}\\ {\tilde{P}}_{2\beta }^{j}\end{array}\right],$$
(15)

with det\(({L}_{\beta }^{j})=\Im ({X}_{1\beta }^{j}{[{X}_{2\beta }^{j}]}^{* })\) the determinant of \({L}_{\beta }^{j}\). Clearly, this solution is only valid when det\(({L}_{\beta }^{j})\,\ne\, 0\). A non-invertible system of equations for a \(\left|{\tilde{{{\Psi }}}}_{\beta }^{j}\right\rangle\) means that all the equations are linearly dependent and that it is sufficient to consider a single equation \(\Re {X}_{\beta }^{j}\cos (\delta {\phi }_{\beta }^{j})-\Im {X}_{\beta }^{j}\sin (\delta {\phi }_{\beta }^{j})={\tilde{P}}_{\beta }^{j}\), from which the trigonometric functions of the relative phase cannot be solved. However, using the identity \({\cos }^{2}(\delta {\phi }_{\beta }^{j})+{\sin }^{2}(\delta {\phi }_{\beta }^{j})=1\), we can obtain the trigonometric functions except for their corresponding quadrant, that is, we have two possible estimates or ambiguities for the sub-states. Thus, in the worst case, where the matrix \({L}_{\beta }^{j}\) cannot be pseudo-inverted for any sub-state, the algorithm leads to a finite set of 2n−1 possible estimates of \(\left|{{\Psi }}\right\rangle\). This problem can be avoided by increasing the number of m of bases for given values of j and β until the systems can be solved. An alternative solution is to consider an extra adaptive measurement to discriminate between the ambiguities. Furthermore, since the probabilities are experimentally estimated by means of a sample of finite size, with high probability the system of equations can be solved due to the finite sample noise, although not necessarily with high accuracy since \({L}_{\beta }^{j}\) might have a bad condition number.

In order to apply our method based on the inductive estimation of sub-states, we need to find an efficient implementation of projections onto states \(\left|{P}_{a\beta }^{j}\right\rangle\). For this purpose, we consider the measurement of two different sets of bases. The first set of n × m bases are given by

$${{{{\mathcal{L}}}}}_{ab}=\left\{{\left|\beta \right\rangle }_{n-b}\otimes \left|{\pm }_{a}\right\rangle \otimes \cdots \otimes \left|{\pm }_{a}\right\rangle \right\},$$
(16)

with 0 ≤ β ≤ 2n−b − 1. These bases can be implemented via local gates, as shown in Fig. 1. These bases define a larger set of projections than required by the method. However, these extra projections can also be used in the algorithm, since \({\left|\beta \right\rangle }_{n-j}\otimes \left|{\pm }_{a}\right\rangle \otimes \cdots \otimes \left|{\pm }_{a}\right\rangle\) form a system of 2j equations at the jth-iteration. The second set of m bases is given by

$${{{{\mathcal{E}}}}}_{a}=\left\{{\left|\beta \right\rangle }_{n-j}\otimes \left|{+}_{a}\right\rangle \otimes {\left|{-}_{a}\right\rangle }^{\otimes j-1},{\left|{-}_{a}\right\rangle }^{\otimes n}\right\},$$
(17)

with 1 ≤ j ≤ n and 0 ≤ β ≤ 2n−j − 1. Despite projections onto members of the bases are local, these observables are entangled since many multi-controlled U gates are needed to implement them, as shown in Fig. 1. In a quantum computer, these gates can be decomposed in terms of CNOT gates65.

Fig. 1: Quantum circuits to implement the measurement bases.
figure 1

a Local measurements \(\{{{{{\mathcal{L}}}}}_{ab}\}\). b Entangled measurements \(\{{{{{\mathcal{E}}}}}_{a}\}\). The single-qubit gate Ua fulfill \(\left|{+}_{a}\right\rangle ={U}_{a}\left|0\right\rangle\) and \(\left|{-}_{a}\right\rangle ={U}_{a}\left|1\right\rangle\).

Full size image

For example, let us consider a 2-qubit system in the state

$$\left|{{\Psi }}\right\rangle ={c}_{0}{e}^{i{\phi }_{0}}{\left|0\right\rangle }_{2}+{c}_{1}{e}^{i{\phi }_{1}}{\left|1\right\rangle }_{2}+{c}_{2}{e}^{i{\phi }_{2}}{\left|2\right\rangle }_{2}+{c}_{3}{e}^{i{\phi }_{3}}{\left|3\right\rangle }_{2}.$$
(18)

Omitting the Kronecker product for compactness, the corresponding sub-states are

$$\begin{array}{l}\left|{\tilde{{{\Psi }}}}_{0}^{1}\right\rangle ={c}_{0}\left|0\right\rangle +{c}_{1}{e}^{i({\phi }_{1}-{\phi }_{0})}\left|1\right\rangle,\end{array}$$
(19)
$$\begin{array}{l}\left|{\tilde{{{\Psi }}}}_{1}^{1}\right\rangle ={c}_{2}\left|0\right\rangle +{c}_{3}{e}^{i(i{\phi }_{3}-{\phi }_{2})}\left|1\right\rangle,\end{array}$$
(20)
$$\begin{array}{l}\left|{\tilde{{{\Psi }}}}_{0}^{2}\right\rangle =\left|0\right\rangle \left({c}_{0}\left|0\right\rangle +{c}_{1}{e}^{i({\phi }_{1}-{\phi }_{0})}\left|1\right\rangle \right)\\\qquad\;\;+\,{e}^{i({\phi }_{2}-{\phi }_{0})}\left|1\right\rangle \left({c}_{2}\left|0\right\rangle +{c}_{3}{e}^{i(i{\phi }_{3}-{\phi }_{2})}\left|1\right\rangle \right),\end{array}$$
(21)

and, the bases to be measured are

$$\begin{array}{l}{{{{\mathcal{L}}}}}_{a0}=\{\left|0\right\rangle \left|{+}_{a}\right\rangle,\left|0\right\rangle \left|{-}_{a}\right\rangle,\left|1\right\rangle \left|{+}_{a}\right\rangle,\left|1\right\rangle \left|{-}_{a}\right\rangle \},\end{array}$$
(22)
$$\begin{array}{l}{{{{\mathcal{L}}}}}_{a1}=\{\left|{+}_{a}\right\rangle \left|{+}_{a}\right\rangle,\left|{+}_{a}\right\rangle \left|{-}_{a}\right\rangle,\left|{-}_{a}\right\rangle \left|{+}_{a}\right\rangle,\left|{-}_{a}\right\rangle \left|{-}_{a}\right\rangle \},\end{array}$$
(23)

or

$${{{{\mathcal{E}}}}}_{a}=\{\left|0\right\rangle \left|{+}_{a}\right\rangle,\left|1\right\rangle \left|{+}_{a}\right\rangle,\left|{+}_{a}\right\rangle \left|{-}_{a}\right\rangle,\left|{-}_{a}\right\rangle \left|{-}_{a}\right\rangle \}.$$
(24)

Our method has been specifically designed to estimate pure states. However, in practice, the states generated in a quantum processor are affected by noise from various sources. The unintended consequence of this is a decrease in the purity of the generated states, which reduces the feasibility of our estimation scheme. However, as we now show our method can also estimate pure states affected by white noise55. This important property is achieved without modifying the measurement bases or increasing their number, which preserves the main features of the method. In the case of white noise, the unknown state to estimate is given by

$$\rho =(1-\lambda )\left|{{\Psi }}\right\rangle \left\langle {{\Psi }}\right|+\frac{\lambda }{d}{\mathbb{I}},$$
(25)

where \({\mathbb{I}}\) is the 2n-dimensional identity operator and λ ≪ 1 is a non-negative noise parameter. The purity of this state is

$$Tr({\rho }^{2})={(1-\lambda )}^{2}+\frac{\lambda }{d}(2-\lambda ).$$
(26)

In the presence of noise we define the quantities \({P}_{a\beta }^{\rho }=\left\langle {P}_{a\beta }^{1}\right|\rho \left|{P}_{a\beta }^{1}\right\rangle\) and \({p}_{\beta }^{\rho }=\left\langle \beta \right|\rho \left|\beta \right\rangle\). These are the probabilities that are registered when projecting the state ρ on the various bases used by our method. Then, we have that

$$\begin{array}{l}{P}_{a\beta }^{\rho }=(1-\lambda ){c}_{2\beta }^{2}{u}_{a}^{2}+\frac{\lambda }{d}{u}_{a}^{2}+(1-\lambda ){c}_{2\beta +1}^{2}{v}_{a}^{2}+\frac{\lambda }{d}{v}_{a}^{2}\\\qquad\;\; +\,2(1-\lambda ){u}_{a}{v}_{a}{c}_{2\beta }{c}_{2\beta +1}\left[\cos ({\varphi }_{a})\cos ({\phi }_{2\beta +1}-{\phi }_{2\beta })\right.\\ \qquad\;\;+\,\left.\sin ({\varphi }_{a})\sin ({\phi }_{2\beta +1}-{\phi }_{2\beta })\right]\end{array}$$
(27)

and that

$$\begin{array}{l}{p}_{2\beta }^{\rho }=(1-\lambda ){c}_{2\beta }^{2}+\frac{\lambda }{d},\end{array}$$
(28)
$$\begin{array}{l}{p}_{2\beta +1}^{\rho }=(1-\lambda ){c}_{2\beta +1}^{2}+\frac{\lambda }{d}\end{array}$$
(29)

These quantities can be estimated through measurements. We now define the quantity

$$\begin{array}{l}{\tilde{P}}_{a\beta }^{\rho }=\left({P}_{a\beta }^{\rho }-(1-\lambda ){c}_{2\beta }^{2}{u}_{a}^{2}-\frac{\lambda }{d}{u}_{a}^{2}\right.\end{array}$$
(30)
$$\begin{array}{rcl}&&-\left.(1-\lambda ){c}_{2\beta +1}^{2}{v}_{a}^{2}-\frac{\lambda }{d}{v}_{a}^{2}\right)/({u}_{a}{v}_{a}),\end{array}$$
(31)

which can be cast as a linear combination of previously obtained quantities, that is,

$${\tilde{P}}_{a\beta }^{\rho }=\left({P}_{a\beta }^{\rho }-{p}_{2\beta }^{\rho }{u}_{a}^{2}-{p}_{2\beta +1}^{\rho }{v}_{a}^{2}\right)/({u}_{a}{v}_{a}).$$
(32)

This leads to the following system of linear equations

$$\left[\begin{array}{cc}\cos {\varphi }_{1}&\sin {\varphi }_{1}\\ \vdots &\vdots \\ \cos {\varphi }_{m}&\sin {\varphi }_{m}\end{array}\right]\left[\begin{array}{c}2(1-\lambda ){c}_{2\beta }{c}_{2\beta +1}\cos \delta {\phi }_{\beta }\\ 2(1-\lambda ){c}_{2\beta }{c}_{2\beta +1}\sin \delta {\phi }_{\beta }\end{array}\right]=\left[\begin{array}{c}{\tilde{P}}_{1\beta }^{\rho }\\ \vdots \\ {\tilde{P}}_{m\beta }^{\rho }\end{array}\right].$$
(33)

The solution of the above system leads to a value for the quantity

$${{{\Lambda }}}_{\beta }^{(\rho )}=2(1-\lambda ){c}_{2\beta }{c}_{2\beta +1}{e}^{i({\phi }_{2\beta +1}-{\phi }_{2\beta })},$$
(34)

which correspond to a linear combination of the quantities \({\tilde{P}}_{a\beta }^{\rho }\). Since \({\tilde{P}}_{a\beta }^{\rho }\) depends on \({P}_{a\beta }^{\rho }\), \({p}_{2\beta }^{\rho }\) and \({p}_{2\beta +1}^{\rho }\), which are quantities that we can estimate from the projective measurements, we can find an estimate for \({{{\Lambda }}}_{\beta }^{(\rho )}\). Taking the absolute value of Eq. (34) we obtain

$$| {{{\Lambda }}}_{\beta }^{(\rho )}{| }^{2}=4{(1-\lambda )}^{2}{c}_{2\beta }^{2}{c}_{2\beta +1}^{2}.$$
(35)

The equation above allows us to estimate the value of the parameter λ. Replacing \({c}_{2\beta }^{2}\) and \({c}_{2\beta +1}^{2}\) from Eqs. (28) and (29) into Eq. (35) we obtain

$$| {{{\Lambda }}}_{\beta }^{(\rho )}{| }^{2}=4\left({p}_{2\beta }^{\rho }-\frac{\lambda }{d}\right)\left({p}_{2\beta +1}^{\rho }-\frac{\lambda }{d}\right),$$
(36)

which is a quadratic equation for the parameter λ with a solution

$$\lambda =\frac{d}{2}\left[{p}_{2\beta }^{\rho }+{p}_{2\beta +1}^{\rho }-\sqrt{{({p}_{2\beta }^{\rho }-{p}_{2\beta +1}^{\rho })}^{2}+| {{{\Lambda }}}_{\beta }^{(\rho )}{| }^{2}}\right],$$
(37)

where we chose the solution with a negative sign since λ must be a small quantity. Let us note that each pair of probabilities p2β and p2β+1 leads to an estimate of λ which generates a set of estimates {λβ}. These estimates differ due to the finite statistic effects on the values of the probabilities. Averaging the various λβ estimates, we obtain a robust estimate λest of λ.

Numerical simulations

We are interested in how similar an unknown state and its estimate are. To quantify this we use the fidelity defined by

$$F=| \langle {{\Psi }}| {{{\Psi }}}_{{{{\rm{est}}}}}\rangle {| }^{2},$$
(38)

between an arbitrary state \(\left|{{\Psi }}\right\rangle\) and its estimate \(\left|{{{\Psi }}}_{{{{\rm{est}}}}}\right\rangle\). A vanishing fidelity indicates that \(\left|{{\Psi }}\right\rangle\) and \(\left|{{{\Psi }}}_{{{{\rm{est}}}}}\right\rangle\) can be perfectly distinguished. A unitary fidelity indicates that \(\left|{{\Psi }}\right\rangle\) and \(\left|{{{\Psi }}}_{{{{\rm{est}}}}}\right\rangle\) are equal. Thus, a good pure-state estimation method will be characterized by high values of fidelity. We are interested in the behavior of the fidelity as a function of the dimension, or equivalently, the number of qubits, and the impact on the fidelity of shot noise or finite-statistics effects.

To study the performance of the proposed method we carried out numerical experiments. For each number of qubits n = 2, 3, …, 10 we randomly generate a set \(\{\left|{{{\Psi }}}^{(i)}\right\rangle \}\) of 100 Haar-distributed pure states, which play the role of the unknown state to be estimated. Afterward, we simulate projective measurements on each \(\left|{{{\Psi }}}^{(i)}\right\rangle\) using the entangled bases, Eq. (17), and the local bases, Eq. (16), with m = 2, 3, and 4. We use angles φa ∈ [0, π/2, π/4, π/3] and amplitudes \({u}_{a}={v}_{a}=1/\sqrt{2}\). These measurements are simulated considering 213 repetitions for each base. Thereafter, using the set of estimated probabilities obtained from the simulated measurements we applied our estimation method to each state \(\left|{{{\Psi }}}^{(i)}\right\rangle\), which leads to a set \(\{|{{{\Psi }}}_{{{{\rm{est}}}}}^{(i)}\rangle \}\) of estimates. Finally, we calculate the value of the fidelity between each unknown state and its estimate. Thereby, we obtain a set of fidelity values for each number n of qubits, the statistics of which we study in the following figures.

Figure 2a shows the median of the fidelity as a function of the number n of qubits for three different numbers nm + 1 of local bases with m = 2, 3, and 4, from bottom to top. For each number of bases, the infidelity decreases as the number of qubits increases. This is a general feature of quantum state estimation methods. Since for each basis the number of repetitions is fixed, the probabilities entering in the systems of linear equations, such as Eq. (15), are estimated with a decreasing accuracy as n increases, which leads to lower fidelity values. Figure 2a also shows that the median value of the fidelity can be increased, while keeping the number of repetitions fixed, using a larger number of local bases. The proposed method requires the pseudo-inverse for every matrix in Eq. (14). The quality of this inversion can be improved by reducing the condition numbers of the linear equation systems, which can be achieved by increasing the number of bases employed in the estimation. An interesting feature of Fig. 2a is that the interquartile range is very narrow for each basis, which indicates that the estimation method generates very similar fidelity values for all the simulated states.

Fig. 2: Median fidelity between arbitrary states and its estimates as a function of the number of qubits.
figure 2

a, c Arbitrary randomly generated states reconstructed with local bases and entangled bases, respectively. b, d Local randomly generated states reconstructed with with local bases and entangled bases, respectively. The colors label the number of bases: For a and b we have 2n + 1 (blue), 3n + 1 (yellow), and 4n + 1 (purple) local bases. For c and d we have 2 (blue), 3 (yellow), and 4 (purple) entangled bases plus the computational basis. Shaded areas represent the corresponding interquartile range.

Full size image

Figure 2c exhibits the median of the fidelity as a function of the number n of qubits for 2, 3, and 4 entangled bases plus the canonical basis. In this case, the estimation method achieves a poor performance in comparison to the use of local bases. The use of 2 entangled bases leads to two equations for every sub-state, which generates an ill-conditioned equation system. Thereby, a poor estimation of the phases of the probability amplitudes is obtained, which propagates through the iterations of the algorithm leading to a low fidelity. However, as Fig. 2c shows, then fidelity can be greatly improved by increasing the number of entangled bases for the estimation. Also, the simulation was carried out employing a fixed ensemble size of 213 for each base. This was done to allow a comparison between our simulations and the results of experiments on the IBM’s quantum processors. The total ensemble employed with the local and entangled bases scales as mn × 213 and k × 213, respectively. Thereby, the total ensemble used with the local bases is much larger than in the case of the entangled bases, which leads to a better estimation via local bases. An increase in the ensemble size used with the entangled bases could lead to an improvement in the fidelity. Nevertheless, the entangled bases are capable of delivering for n = 6 a fidelity close to 0.9. Figure 2b and d also displays the fidelity as a function of the number of qubits for local and entangled bases, respectively. However, the simulation has been carried out on a set of completely separable states. In both cases we obtain an increase in the performance of the proposed estimation method, achieving fidelities of around 0.95 for 10 qubits when employing local bases. This result indicates that the performance of the proposed estimation method on the set of entangled states should be closer to the depicted in Fig. 2a and c.

Figure 3a shows the effect of increasing the number of repetitions on the median estimation fidelity. The number of repetitions was increased from 213, which is the number used in the previous simulations, to 223. In this case, we have considered up to 20 qubits and 2n + 1 and 4n + 1 local bases. As is apparent from Fig. 3a, for n ≤ 10 the fidelity shows a large increase when compared to the previous simulations shown in Fig. 2a. For instance, for 10-qubit states 2n + 1 bases leads to an average fidelity of 0.999. However, for n ≈ 20, the fidelity is again close to 0.85. The overall impact on the fidelity is similar for both numbers of bases. Figure 3b shows the average estimated time for a n-qubit state, that is, the time required by a portable personal computer to process the experimentally acquired data and output an estimate. According to this figure, for 20 qubits, a time of the order of minutes is required. This short time is a consequence of the fact that the method does not require solving an optimization problem to produce an estimate, but rather requires the solution of a sequence of systems of linear equations.

Fig. 3: Fidelity and execution time as a function of the number of qubits.
figure 3

a Median fidelity between randomly generated states and its estimates as a function of the number of qubits. b Average time used by the algorithm to generate an estimate of a state as a function of the number of qubits. Colors label the number of bases, blue 2n + 1 and purple 4n + 1 local bases.

Full size image

Experimental realization

We carried out an experimental demonstration of our estimation method using the quantum processor ibmq_montreal developed by IBM. This NISQ device has online access and can be programmed with Qiskit, an open-source development framework for working with quantum computers in Python programming language.

To test our estimation method we prepare the following n-qubit states

$$\left|{{{\Phi }}}_{1}^{n}\right\rangle =\frac{1}{\sqrt{5\times {2}^{n}}}{\left(\left|0\right\rangle -2{e}^{i\pi /4}\left|1\right\rangle \right)}^{\otimes n},$$
(39)
$$\left|{{{\Phi }}}_{2}^{n}\right\rangle =\frac{1}{\sqrt{5\times {2}^{n}}}{\left(\left|0\right\rangle +2{e}^{i\pi /4}\left|1\right\rangle \right)}^{\otimes n},$$
(40)

and

$$\left|{{{\Phi }}}_{3}^{n}\right\rangle =\left\{\begin{array}{ll}\frac{1}{{2}^{n/4}}{\left(\left|00\right\rangle +\left|11\right\rangle \right)}^{\otimes n/2},\quad &n\,{{{\rm{even}}}}.\\ \frac{1}{{2}^{(n-1)/4}}{\left(\left|00\right\rangle +\left|11\right\rangle \right)}^{\otimes (n-1)/2}\otimes \left|0\right\rangle,\quad &n\,{{{\rm{odd}}}}.\end{array}\right.$$
(41)

States \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) are strictly separable and thus they can be prepared efficiently with the use of a single-qubit gate acting on each qubit. The state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) has maximal entanglement between pairs of qubits. For n even (odd), this state can be prepared with n/2 ((n − 1)/2) CNOT gates and n/2 ((n − 1)/2) single-qubit gates. Thus, states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\), \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\), and \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) can be generated by circuits with short depth. However, the state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) is generated with a lower preparation fidelity than states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) due to the fact that CNOT gates are implemented with an error much larger than the one achieved in the implementation of single-qubit gates. We also perform the estimation of n-partite GHZ states,

$$\left|{{{\Phi }}}_{4}^{n}\right\rangle =\frac{1}{\sqrt{2}}\left({\left|0\right\rangle }^{\otimes n}+{\left|1\right\rangle }^{\otimes n}\right).$$
(42)

These states are implemented by a long-depth circuit that contains n − 1 CNOT gates. Consequently, the fidelity achieved in its generation can be very low. Thus, we focus our study in small number of qubits n = 2, 3, 4 for both local and entangled bases. Each basis was measured employing 1013 repetitions, that is, the size of the ensemble of equally prepared copies of the unknown quantum state. This is currently the maximal sample size that can be employed in IBM’s quantum processors.

Figures 4, 5, 6, and 7 summarize the results obtained by implementing our estimation method on IBM’s quantum processor for states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\), \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\), \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\), \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\), respectively. In each figure, the upper and lower rows show experimental results obtained through local and entangled bases, correspondingly. The figures display the fidelity between the theoretical states and their estimates, the experimental fidelity evaluated by direct fidelity estimation, and the estimated purity, all of them as a function of the number n of qubits for the cases of mn + 1 local bases (with m = 2, 3 and 4) and m entangled bases plus the computational base. We use angles φa ∈ [0, π/2, π/4, 3π/4] and amplitudes \({u}_{a}={v}_{a}=1/\sqrt{2}\). Shaded areas correspond to the experimental errors (see methods). We use direct fidelity estimation27,28 to compare the estimated state with the actually prepared state. We use a number M = 25 × n of observables for each number of qubit n for the experimental implementation of direct fidelity estimation (see methods).

Fig. 4: Experimental results of the tomographic reconstruction of the local state \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) as a function of the number of qubits.
figure 4

a Theoretical fidelity, b experimental fidelity, and c purity achieved by the tomographic protocol with mn + 1 local bases (upper row, m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). d Theoretical fidelity, e experimental fidelity, and f purity achieved by the tomographic protocol with m entangled bases plus the computational base (lower row, m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). Shaded areas represent the standard deviation obtained by bootstrap method.

Full size image
Fig. 5: Experimental results of the tomographic reconstruction of the local state \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) as a function of the number of qubits.
figure 5

a Theoretical fidelity, b experimental fidelity, and c purity achieved by the tomographic protocol with mn + 1 local bases (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). d Theoretical fidelity, e experimental fidelity, and f purity achieved by the tomographic protocol with m entangled bases plus the computational base (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). Shaded areas represent the standard deviation obtained by bootstrap method.

Full size image
Fig. 6: Experimental results of the tomographic reconstruction of the entangled state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) as a function of the number of qubits.
figure 6

a Theoretical fidelity, b experimental fidelity, and c purity achieved by the tomographic protocol with mn + 1 local bases (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). d Theoretical fidelity, e experimental fidelity, and f purity achieved by the tomographic protocol with m entangled bases plus the computational base (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). Shaded areas represent the standard deviation obtained by bootstrap method.

Full size image
Fig. 7: Experimental results of the tomographic reconstruction of the GHZ state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) as a function of the number of qubits.
figure 7

a Theoretical fidelity, b experimental fidelity, and c purity achieved by the tomographic protocol with mn + 1 local bases (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). d Theoretical fidelity, e experimental fidelity, and f purity achieved by the tomographic protocol with m entangled bases plus the computational base (m = 2 (blue), m = 3 (yellow), and m = 4 (purple)). Shaded areas represent the standard deviation obtained by bootstrap method.

Full size image

According to Figs. 4a and 5a, the estimation of local states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) through local bases lead to values that are comparable with the numerical simulations shown in Fig. 2a and b, which only consider noise due to finite statistics. For the particular case of n = 10 qubits, the theoretical predictions for the estimation fidelity using local bases are within the intervals [0.88, 0.93] for randomly generated states, and [0.95, 0.96] for randomly generated separable states, while the experiment leads to an estimation fidelity in the interval [0.88, 0.91]. Figures 4b and 5b show experimental fidelity values that are below the theoretical fidelity, which for 10 qubits in both cases is still above 0.7. Let us note that the error in the experimental fidelity, which was obtained through direct fidelity estimation, is in the order of 0.02. The purity, displayed in Figs. 4c and 5c, is above 0.92 in both cases.

The estimation of entangled states \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\) by local bases, exhibited in Figs. 6a and 7a, respectively, also leads to good fidelities, albeit lower than in the case of states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\). This is to be expected since states \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\) exhibit entanglement and thus are generated applying several CNOT gates, which increase the preparation error. Nevertheless, for the state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) with n = 10 our method provides an estimation fidelity close to 0.8. A comparison between the estimated fidelity for the bi-local state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and GZH state \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\) for n = 2, 3, and 4 qubits indicates that these entangled states are estimated with similar fidelities. This indicates that the method delivers similar quality estimation for states with a different type of entanglement. Figures 6b and 7b exhibit values of the experimental fidelity that are again consistently below the values of the theoretical fidelity. Nevertheless, these values are above 0.6. The purity of the estimates obtained for state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) are above 0.6 up to 8 qubits. For 9 and 10 qubits there is a sharp decrease in the purity. In the case of the GHZ state \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\), the purity is above 0.8.

The experimental realization of our estimation method shows that the use of local bases allows us to estimate pure states of large numbers of qubits. The fidelity of the estimation is meanly constrained by ensemble size and number of local bases. An increase of any of these quantities leads to an improvement in the quality of the estimation. In particular, states of larger number of qubits can be reliably estimated by increasing the number of local bases. As is shown in Figs. 4a, 5a, 6a, and 7a the use of 4n + 1 local bases leads to higher values of the estimation fidelity than the cases of 2n + 1 and 3n + 1 local bases.

The use of entangled bases shows a different picture, where the physical realization of the method and the particular class of states to be estimated lead to a reduction of the theoretical fidelity when compared to the case of local bases. Figures 4d and 5d exhibit the theoretical fidelity achieved for separable states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\), respectively. A comparison for n = 4 with Figs. 4a and 5a shows a decrease of the theoretical fidelity from ~0.995 to 0.7 in the case of state \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and from approximately 0.995 to 0.2 in the case of state \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\). A similar decrease can also be observed in the case of states \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\) in Figs. 6d and 7d, although in the latter the decrease is less abrupt. The low achieved theoretical fidelity finds its origin in the circuits employed to implement the measurements on the entangled bases. For n = 2, 3, and 4, 2 entangled bases require the use of 2, 7, and 27 CNOT gates, respectively. This number becomes 3027 for n = 10. Since the CNOT gate has a high error rate in NISQ devices, the implemented bases differ significantly from the actual bases \({{{{\mathcal{E}}}}}_{a}\) to be implemented. Besides, this also affects the generated states. Figure 5d shows that for n = 4 the use of 2, 3, and 4 entangled bases plus the computational base lead to a theoretical fidelity value of 0.2. This originates in the state \(\left|{{{\Phi }}}_{2}^{3}\right\rangle\) to be reconstructed, which exhibits ill-conditioned equation systems. In Fig. 6d, the use of two entangled bases plus the computational base leads to a theoretical fidelity of 0.25 while for more entangled bases the value increases up to approximately 0.7. The experimental fidelity exhibited in Figs. 4e, 5e, 6e, and 7e a behavior analogous to theoretical fidelity, albeit reaching similar or slightly lower values. As Figs. 4f and 5f indicate, the purity obtained in the estimation of states \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) is in the interval [0.75, 1], which is lower than the purity obtained using the separable bases. However, a much abrupt decrease in the purity can be observed in Figs. 6f and 7f for entangled states \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\), where the achieved values for n = 4 are in the order of 0.5 and 0.25, respectively. This feature is due to the large number of CNOT gates that are used to generate states and measurement bases. However, for n = 2 and n = 3, the purity is still in the interval [0.7, 1].

Discussion

Estimating states of d-dimensional quantum systems requires a minimal number of measurement outcomes that scales quadratically with the dimension. In the case of composite systems, such as quantum computers, the scaling becomes exponential in the number of subsystems, which makes the estimation by generic methods unfeasible except for few-component systems. NISQ computers are characterized by low-accuracy entangling gates, which are required for implementing measurements of arbitrary observables, and by a fixed number of repetitions, which constrains the size of statistical samples and decreases the estimation accuracy as the number of qubits increases. Therefore, estimating n-qubit states of NISQ computers is a difficult task.

We have proposed a method to estimate pure states of n-qubit systems that is well suited for NISQ computers. The method is based on the reconstruction of the so-called sub-states \(|{\tilde{{{\Psi }}}}_{\beta }^{j}\rangle\), which for an unknown state \(\left|{{\Psi }}\right\rangle\) are n sets of 2n−j non-normalized states, one set for each j = 1, 2, . . . , n. If we know the sub-states \(|{\tilde{{{\Psi }}}}_{\beta }^{j}\rangle\) for a fixed j, then it is easy to reconstruct the sub-states with \(|{\tilde{{{\Psi }}}}_{\beta }^{j+1}\rangle\) from the results of well-defined set of measurements. The method begins by measuring a set of projectors that allow us to reconstructs all the sub-states for j = 1, and iterating, the rest of sub-states until reaching \(|{\tilde{{{\Psi }}}}_{\beta }^{n}\rangle\), the estimate of the unknown pure state. Interestingly, the proposed method is also capable of estimating pure states that are transformed into mixed states by the action of white noise. This does not require changing the measurement bases or increasing their number.

The set of measurements employed by the proposed method corresponds to projective measurements onto a set of bases. We have first shown that mn + 1 bases, with m an integer number equal to or greater than 2, allow us to estimate all pure states. Thereby, a total of (mn + 1)2n measurement outcomes is required. This number compares favorably with other estimation methods. Mutually unbiased bases, SIC-POVM, and compressed sensing require (2n + 1)2n, 22n, and in the order of 22nn2 measurement outcomes, respectively, to reconstruct pure states. Recently, an adaptive version of compressive (ACS) sensing has been introduced66, which does not require a priori information about the state to be reconstructed. Numerical simulations up to 7 qubits have led to the conjecture67 that r-rank states can be estimated via ACS with 4r + 1 bases, that is, only 5 bases in the case of pure states. ACS determines at each iteration whether the measurements performed are informationally complete and adapts the next measurement basis. Each one of these stages requires the solution of optimization problems, such as maximum likelihood estimation, in a space with an exponentially increasing dimension. This severely limits the ability to solve these problems. For a small number of qubits, advanced optimization techniques and programming have allowed quick and reliable implementation of maximum likelihood estimation68. For instance, 10-qubit pure states affected by white noise can be estimated in the order of 104 (s)68. However, as the dimension increases, the exponential behavior prevails. At this point our method provides a solution, since it does not require optimization-based data post-processing. Instead, it requires the solution of small systems of linear equations. For m = 2, the solutions are analytical. Thereby, the trade-of is clear. For a small number n of qubits, where the optimization problems can be solved efficiently, ACS is a very good solution. For a large number of qubits, our method provides a feasible solution by increasing the number of bases to be measured while decreasing the classical computational cost of the estimation process. According to our simulations, estimating an arbitrary 20-qubit state requires a computation time of the order of 3 min on a portable personal computer. This is a direct consequence of the separability of the states belonging to the bases, which makes it possible to use reshaping techniques69. In this way, it is possible to significantly reduce the computation time of matrix products as well as the memory usage.

The mn + 1 bases are local, that is, they can be cast as the tensor product of n single-qubit bases. Projective measurements of the states of these bases are carried out on a NISQ computer by applying local gates onto each qubit followed by a projection onto the computational base. Thus, no entangling gates are necessary. In principle, 2n + 1 bases are enough. We have shown, by means of numerical simulations, that the use of a larger number of local bases increases the estimation accuracy for larger numbers of qubits. A similar effect can be achieved by increasing the ensemble size while using a low number of bases. For instance, in the case of 10-qubit states increasing the ensemble size from 213 to 223 leads to an increase in average fidelity from 0.9 to 0.999. However, if high-accuracy entangling gates are available, then the set of measurements employed by the proposed method can be reduced to projective measurements on m entangled bases plus the computational base, with m an integer number equal to or greater than 2. In this case, an increase in the number of entangled bases also allows us to increase the estimation accuracy.

We tested the proposed estimation method in the IBM’s quantum processing units. The estimation of local pure states via local bases up to 10 qubits provides a theoretical fidelity that agrees with the numerical simulations and is above 90%. The experimental fidelity, which is obtained by direct fidelity estimation, is for 10 qubits above 0.7. The purity is above 0.9. The estimation of entangled states via local bases was also tested. In this case, the estimation of tensor product of two-qubit Bell states up to n = 10 led to theoretical fidelities above 0.7. The experimental fidelity is between 0.5 and 0.7. In the case of the state \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) purity reaches a value close to 0.6 up to 8 qubits after which it exhibits a steep decrease to values below 0.4. In this case, however, the preparation of the state is affected by large errors due to the use of entangling gates. Finally, we tested the estimation of a GHZ state for n = 2, 3, and 4 qubits obtaining theoretical fidelities above 0.86. The experimental fidelity is between 0.7 and 0.8 while purity is above 0.8. Entangled bases were also tested. Nevertheless, due to the massive use of entangling gates in the preparation and measurement stages mixed results were achieved, as expected. States \(\left|{{{\Phi }}}_{1}^{n}\right\rangle\), \(\left|{{{\Phi }}}_{3}^{n}\right\rangle\) and \(\left|{{{\Phi }}}_{4}^{n}\right\rangle\) exhibit theoretical and experimental fidelities above 0.6 up to 4 qubits. The purity shows a steep decrease with respect to the estimation via separable bases. The estimation of state \(\left|{{{\Phi }}}_{2}^{n}\right\rangle\) leads to values of theoretical and experimental fidelities above 0.9 for 2 and 3 qubits. However, for 4 qubits there is a steep decrease in these figures. Purity is above 0.8. The estimation via k entangled bases plus the computational base, where k is high-enough but does not scale with the number n of qubits, can lead to a large fidelity in future fault-tolerating quantum architectures70, where large numbers of qubits and high-accuracy entangling gates are required. Also, it might be possible to increase the estimation fidelity by using alternative implementations of the multi-control NOT gates that employ ancillary qubits. This allows a reduction in the depth of the circuit required to implement measurements on the proposed entangled bases65.

The main characteristic of our protocol is its scalability, which is a consequence of the a priori information about the states to be estimated. This is an approach common among many estimation methods such as, for instance, compressed sensing and matrix product states. The purity assumption is reasonable in systems that are able to prepare high purity states or in systems whose purity can be certified, for instance, through randomized benchmarking71. However, this is not necessarily true on NISQ devices. Decoherence and gate errors can reduce the purity of the generated state in such a way that it becomes a mixed state. It has been shown that experimental results can be improved by reducing the error of the raw data by employing error-mitigation techniques72,73. Combining scalable error-mitigation methods with our estimation method could extend its applicability to high-noise systems. Furthermore, if errors can be approximated by white noise, then our method can still be used without modifying the measurement bases.

Methods

Numerical simulations

In order to compare the estimate \(\left|{{{\Psi }}}_{{{{\rm{est}}}}}\right\rangle\) provided by our method with the unknown state \(\left|{{\Psi }}\right\rangle\), we employ the fidelity \(F(\left|{{\Psi }}\right\rangle,| {{{\Psi }}}_{{{{\rm{est}}}}})=| \langle {{\Psi }}| {{{\Psi }}}_{{{{\rm{est}}}}}\rangle {| }^{2}\). A state-independent figure of merit of the performance is obtained by calculating the mean of the fidelity on the Hilbert space \({{{\mathcal{H}}}}\) of the unknown states, that is,

$$\bar{F}={\int}_{{{{\mathcal{H}}}}}F(\left|{{\Psi }}\right\rangle,| {{{\Psi }}}_{{{{\rm{est}}}}})f(\left|{{\Psi }}\right\rangle )d{{\Psi }},$$
(43)

where \(f(\left|{{\Psi }}\right\rangle )\) is a probability distribution function for the states to be estimated.

We assume that the unknown states are uniformly generated by a source. For a given number n ≥ 2 of qubits, we create the sets \({{{\Omega }}}_{{{{\rm{arb}}}}}^{n}\) and \({{{\Omega }}}_{{{{\rm{sep}}}}}^{n}\) that contain Ns unknown arbitrary and separable states, respectively. The states in each set are generated according to a Haar-uniform distribution and estimated using local bases, Eq. (16), and entangled bases, Eq. (17). This generates four sets of fidelity values, each with Nest elements, whose statistics can be calculated. For example, the mean fidelity \(\bar{F}\) can be estimated by means of

$$\bar{F}=\mathop{\sum }\limits_{i=1}^{{N}_{{{{\rm{est}}}}}}\frac{1}{{N}_{{{{\rm{est}}}}}}F(\left|{{{\Psi }}}_{i}\right\rangle,\left|{{{\Psi }}}_{i,{{{\rm{est}}}}}\right\rangle ).$$
(44)

Other statistical variables, such as median and interquartile range, can be similarly estimated.

The value \(F(\left|{{{\Psi }}}_{i}\right\rangle,\left|{{{\Psi }}}_{i,{{{\rm{est}}}}}\right\rangle )\) is obtained by simulating the projective measurements on the bases \(\{\left|{b}_{j}\right\rangle \}\), local or entangled. A basis \(\{\left|{b}_{j}\right\rangle \}\) and a state \(\left|{{{\Psi }}}_{i}\right\rangle\) define a multinomial probability distribution pi,j = ∣〈Ψi∣bj〉∣2, which is employed to generate a sample of total size 213. The probabilities are estimated as

$${p}_{i,j}=\frac{{n}_{i,j}}{{2}^{13}},$$
(45)

where ni,j is the number of times that the jth outcome is present in the sample. The set of probabilities {pi,j} is then used to obtain the estimate \(\left|{{{\Psi }}}_{i,{{{\rm{est}}}}}\right\rangle\) and thereafter the fidelity value \(F(\left|{{{\Psi }}}_{i}\right\rangle,\left|{{{\Psi }}}_{i,{{{\rm{est}}}}}\right\rangle )\). We set Nest = 102. A sample size of 213 is chosen to agree with the sample size available at IBM’s quantum processors.

The inner products entering in the \({X}_{\alpha \beta }^{j}\) quantities, see Eq. (12), were efficiently calculated by reshaping techniques68,69. These quantities are then used to generate the matrices to be inverted.

Experimental realization

We also tested the performance of the proposed method on the IBM quantum processor ibmq_montreal, which is freely accessible through its cloud service. This quantum processor can implement any unitary transformation acting on its qubits. However, it only allows measurements on a computational basis \(\{\left|j\right\rangle \}\). Therefore, measuring on a non-computational basis \(\{\left|{b}_{j}\right\rangle \}\) requires applying the unitary transformation U, which implements a change of basis, onto the set of qubits and then measuring the resulting state on the computational basis. Using the open-source software development kit Qiskit 62, we created quantum circuits to implement the unitary transformations required to measure on local and entangled bases, Eqs. (16) and Eq. (17), respectively. These quantum circuits are compiled into equivalent circuits involving only universal gates. Similarly, quantum circuits were designed to implement the unitary transformations that generate the unknown states used on the various tests.

Measurements on a given basis are repeated 213 times to generate a sample from which the required statistics can be estimated. This procedure, however, introduces noise in the inference of the probabilities, which affects the performance of the proposed method. This finite statistic noise can be reduced by increasing the number of repetitions. Additionally, systematic errors affect the preparations of the gates, which also reduced the experimental performance of the proposed method. Thereby, transformations that generate states and transformations that change bases are approximately implemented. In particular, due to the highly non-trivial preparation of the CNOT gate, this has an error higher than the one of local gates. This is the main source of error that affects the performance of the proposed method when implemented in the quantum processor. The error in the device can be modeled by white noise \({{{\mathcal{E}}}}\) and the average error per gate r,

$${{{{\mathcal{E}}}}}_{{{{\rm{noise}}}}}({\left|0\right\rangle }^{\otimes n}\left.\right\rangle )=\left(1-\frac{{2}^{nr}}{{2}^{n}-1}\right){{{\mathcal{E}}}}({\left|0\right\rangle }^{\otimes n})+\frac{{2}^{nr}}{{2}^{n}-1}I.$$
(46)

The average error per gate of the quantum device ibmq_montreal is given by r1 = 8 × 10−4 for local gates and r2 = 2 × 10−2 for CNOT gates. IBM provides and updates these values periodically.

The errors in the fidelity F and the noise λ of the experimental estimates are calculated using the bootstrap method. We generate simulated probabilities by resampling the experimentally detected counts. This procedure is repeated 103 times, generating a set of estimators of the state \(\{|{{{\Psi }}}_{{{{\rm{est}}}}}^{j}\rangle \}\), fidelities {Fj} and noise {λj}. The error is given by the standard deviation of these quantities.

The experimental estimation of the fidelity is obtained by direct fidelity estimation27,28. This allows us to estimate the fidelity between an experimental state ρ and a theoretical pure state \(\sigma =\left|\psi \right\rangle \left\langle \psi \right|\) measuring an efficient number of observables. Let orthonormal Pauli operators {Pj} that will be used in the estimation. First, the expected values σj = Tr(Pjσ) are evaluated numerically. After that, we generate a set of M random Pauli operators by the probability distribution \({\sigma }_{m}^{2}\). These operators, labeled by indices {km}, with m = 1, …, M, are measured experimentally on the state ρ in order to obtain the expected values \({\rho }_{{k}_{m}}=Tr(\rho {P}_{{k}_{m}})\). The estimation of the fidelity is given by

$${F}_{{{{\rm{dfe}}}}}=\frac{1}{M}\mathop{\sum }\limits_{m=1}^{M}\frac{{\rho }_{{k}_{m}}}{{\sigma }_{{k}_{m}}}.$$
(47)

The error ΔFdfe in the estimation of the fidelity is calculated using the formula provided by B. P. Lanyon et al.52,

$${{\Delta }}{F}_{{{{\rm{dfe}}}}}^{2}=\frac{1}{{M}^{2}}\mathop{\sum }\limits_{m=1}^{M}\frac{{\rho }_{{k}_{m}}^{2}}{{\sigma }_{{k}_{m}}^{2}}-\frac{1}{{M}^{2}(M-1)}\sum\limits_{\mathop {m,{m^{\prime}} = 1}\limits_{m \ne {m^{\prime}}} }^M \frac{{\rho }_{{k}_{m}}}{{\sigma }_{{k}_{m}}}\frac{{\rho }_{{k}_{m^{\prime} }}}{{\sigma }_{{k}_{m^{\prime} }}}.$$
(48)