Fast and simple quantum state estimation

Quantum state estimation is the process of reconstructing the density matrix from measurements performed over an ensemble of identically prepared quantum systems. In the early days of quantum theory, Pauli posed the question of whether position and momentum probability distributions univocally determine the state of a quantum particle [1], something that holds in classical mechanics. However, quantum states belong to an abstract Hilbert space whose dimension exponentially increases with the system's number of particles. Thus, more information than classically expected is required to determine the state. Since then, there has been increased interest in finding methods to estimate the state of a quantum system from a given set of measurements and several solutions appeared. For instance, standard state tomography [2] reconstructs d-dimensional density matrices from O(d³) rank-one projective valued measures (PVM), whereas mutually unbiased bases (MUBs) [3, 4] and symmetric informationally complete (SIC) positive operator valued measures (POVM) [5] do the same task with O(d²) rank-one measurement projectors. In general, any tight quantum measurement [6], equivalently any complex projective two-design, is informationally complete [7].

Quantum state tomography finds applications in communication systems [8], dissociating molecules [9] and characterization of optical devices [10]. It is a standard tool for verification of quantum devices, e.g. estimating fidelity of two photon CNOT gates [11], and has been used to characterize quantum states of trapped ions [12], cavity fields [13], atomic ensembles [14] and photons [15].

Aside from the experimental procedure of conducting a set of informationally complete measurements on a system, quantum tomography requires an algorithm for reconstructing the state from the measurement statistics. From a variety of techniques proposed, the approaches featuring in the majority of experiments are variants of linear inversion (LI) and maximum-likelihood quantum state estimation (MLE) [16]. As its name suggests, with LI one determines the state of the quantum system under consideration by inverting the measurement map solving a set of linear equations with the measurement data as input. For relevant families of informationally-complete set of measurements, analytical expressions for the inverse maps are known, significantly speeding up the whole reconstruction effort, see e.g. [17]. MLE consists in finding the state that maximizes the probability of obtaining the given experimental data set, among the entire set of density matrices. Within the different implementations of this basic last idea, those currently achieving the best runtimes are variants of a projected-gradient-descent scheme, see [18, 19]. Algorithms based on variants of linear inversion [20, 21] are typically faster than those implementing MLE when the inversion process is taken from already existing expressions [22]. On the other hand, when restrictions on the rank of the state being reconstructed apply, techniques based on the probabilistic method of compressed-sensing have proven to be very satisfactory [23–25]. In particular, the statistics based on five rank-one projective measurements is good enough to have a high–fidelity reconstruction of rank-one quantum states, even under the presence of errors in both state preparation and measurement stages [26]. In this work, we present a general method for quantum state estimation achieving better fidelities than the state-of-the-art implementations of MLE.

This paper is organized as follows. In section 2, we introduce the main ingredient of our algorithm: the physical imposition operator, a linear operator having an intuitive geometrical interpretation. In section 3, we present our iterative algorithm for quantum state estimation based on the physical imposition operator and prove its convergence. In section 3.1 we show that for a wide class of quantum measurements, which include MUBs and tensor product of generalized Pauli observables for N qudit systems, convergence is achieved in a single iteration. In section 4, we numerically study the performance of our algorithm in terms of runtime and fidelity estimation, finding an improvement with respect to the most efficient MLE-based method, as far as we know. Finally, in section 5 we provide conclusions and future lines of research. Proofs of all our results are presented in appendix A.

Consider an experimental procedure $\mathcal{P}$ that prepares a quantum system in some unknown state ρ. Let us assume that, given some prior knowledge about $\mathcal{P}$, our best guess for the state of the system ρ₀, which could be even the maximally mixed state in absence of prior information. Next, we perform a POVM measurement A composed by m_A outcomes, i.e. $A={\left\{{E}_{i}\right\}}_{i{\leqslant}{m}_{A}}$ on an ensemble of systems independently prepared according to $\mathcal{P}$, obtaining the outcome statistics $\overrightarrow {p}={\left\{{p}_{i}\right\}}_{i{\leqslant}{m}_{A}}$. Given this newly acquired information,

how can we update ρ₀ to reflect our new state of knowledge about the system?

To tackle this question, we introduce the physical imposition operator, a linear map that replaces the initial predictions about observable A contained in ρ₀ with an experimentally observed probability p_i .

Definition 2.1. (Physical imposition operator). Let $A={\left\{{E}_{i}\right\}}_{i{\leqslant}{m}_{A}}$ be a POVM acting on a d-dimensional Hilbert space ${\mathcal{H}}_{d}$ and let $\overrightarrow {p}\in {\mathbb{R}}^{{m}_{A}}$ be a probability vector. The physical imposition operator associated to E_i and p_i is the linear map

$$\begin{equation}{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)=\rho +\frac{\left({p}_{i}-\mathrm{Tr}\left[\rho {E}_{i}\right]\right){E}_{i}}{\mathrm{Tr}\left({E}_{i}^{2}\right)},\end{equation} \tag{ 1 }$$

for every i ⩽ m_A .

In order to clarify the meaning of the physical imposition operator (1) let us assume for the moment that A is a projective measurement. In such a case, operator ${T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)$ takes a quantum state ρ, removes the projection along the direction E_i , i.e. it removes the physical information about E_i stored in the state ρ, and imposes a new projection along this direction weighted by the probability p_i . Here, p_i can be either taken from experimental data or simulated by Born's rule with respect to a target state to reconstruct. Note that operator ${\rho }^{\prime }={T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)$ reflects the experimental knowledge about the quantum system. As we will show in section 3, a successive iteration of PIO along an informationally complete set of quantum measurements allows us to reconstruct the quantum state. For POVM in general, operator (1) does not entirely impose the information about the outcome. However, after several imposition of all involved operators PIO the sequence of quantum states successfully converges to a quantum states containing all the physical information, as we demonstrate in theorem 3.1. To simplify notation, along the work we drop the superscript p_i in ${T}_{{E}_{i}}^{{p}_{i}}$ when the considered probability p_i is clear from the context.

Let us now state some important facts about PIOs that easily arise from definition 2.1. From now on, $\mathfrak{D}\left(\rho ,\sigma \right){:=}\mathrm{Tr}\left[{\left(\rho -\sigma \right)}^{2}\right]$ denotes the Hilbert–Schmidt distance between states ρ and σ.

Proposition 2.1. The following properties hold for any POVM ${\left\{{E}_{i}\right\}}_{i{\leqslant}{m}_{A}}$ and any ρ, σ acting on ${\mathcal{H}}_{d}$:

• (a)  
  Imposition of physical information: $\mathrm{Tr}\left[{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right){E}_{i}\right]={p}_{i}$.
• (b)  
  Composition: ${T}_{{E}_{j}}^{{p}_{j}}{\circ}{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)={T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)+{T}_{{E}_{j}}^{{p}_{j}}\left(\rho \right)-\rho -\left({p}_{i}-\mathrm{Tr}\left(\rho {E}_{i}\right)\right)\mathrm{Tr}\left({E}_{i}{E}_{j}\right){E}_{j}/\mathrm{Tr}{\left({E}_{j}\right)}^{2}$.
• (c)  
  Non-expansiveness: $\mathfrak{D}\left({T}_{{E}_{j}}^{{p}_{j}}\left(\rho \right),{T}_{{E}_{j}}^{{p}_{j}}\left(\sigma \right)\right){\leqslant}\mathfrak{D}\left(\rho ,\sigma \right)$.

Some important observations arise from proposition 2.1. First, for j = i in the above item (b) we find that

$$\begin{equation}{T}_{{E}_{i}}^{2}\left(\rho \right)={T}_{{E}_{i}}\left(\rho \right),\end{equation} \tag{ 2 }$$

for any ρ, so operator ${T}_{{E}_{i}}$ is an orthogonal projection, for every i ⩽ m_A and any POVM ${\left\{{E}_{i}\right\}}_{i{\leqslant}{m}_{A}}$. Note that any quantum state $\sigma ={T}_{{E}_{i}}\left(\rho \right)$ is a fixed point of ${T}_{{E}_{i}}$, i.e. ${T}_{{E}_{i}}\left(\sigma \right)=\sigma$, which simply arises from (2). Roughly speaking, quantum states already having the physical information we want to impose are fixed points of the map ${T}_{{E}_{j}}$. This key property allows us to apply dynamical systems theory [27] to study the tomographic problem. We consider the alternating projection method, firstly studied by von Neumann [28] for the case of two alternating projections and generalized by Halperin to any number of projections [29].

In theorem 3.1, we will show that composition of all physical imposition operators associated to an informationally complete set of POVM produces a linear map having a unique attractive fixed point, i.e., the solution to the quantum state tomography problem. The uniqueness of the fixed point guarantees a considerable speed up of the method in practice, as any chosen seed monotonically approaches to the solution of the problem.

To simplify notation, we consider a single physical imposition operation ${\mathcal{T}}_{A}$ for an entire POVM A, defined as follows

$$\begin{equation}{\mathcal{T}}_{A}={T}_{{E}_{{m}_{A}}}{\circ}\dots {\circ}{T}_{{E}_{1}}.\end{equation} \tag{ 3 }$$

Up to a constant factor proportional to identity, that we omit, operator ${\mathcal{T}}_{A}$ reduces to

$$\begin{equation}{\mathcal{T}}_{A}\left(\rho \right)=\sum\limits _{i=1}^{{m}_{A}}{T}_{{E}_{i}}\left(\rho \right),\end{equation} \tag{ 4 }$$

for any PVM A, what follows from considering (3) and proposition 2.1. This additive property holding for PVM measurements plays an important role, as it helps to reduce the runtime of our algorithm. Precisely, this fact allows us to apply Kaczmarz method [30] instead of Halpering alternating projection method, for any informationally complete set of PVM. Kaczmarz method considers projections over the subspace generated by the intersection of all associated hyperplanes, defined by the linear system of equations (Born's rule).

Let us introduce another relevant concept

Definition 2.2. (Generator state). Given a POVM $A={\left\{{E}_{i}\right\}}_{i{\leqslant}{m}_{A}}$ and a probability vector $\overrightarrow {p}\in {\mathbb{R}}^{{m}_{A}}$, a quantum state ρ_gen is called generator state for $\overrightarrow {p}$ if Tr(ρ_gen E_i ) = p_i , for every i ⩽ m_A .

Note that ρ_gen is a fixed point of ${\mathcal{T}}_{{E}_{i}}$, according to (3) and proposition 2.1. State ρ_gen plays an important role to implement numerical simulations, as it guarantees to generate sets of probability distributions compatible with the existence of a positive semidefinite solution to the quantum state tomography problem.

To end this section, note that map ${\mathcal{T}}_{A}$ defined in (4) has a simple interpretation in the Bloch sphere for a qubit system, see figure 1. The image of ${\mathcal{T}}_{A}$, i.e. ${\mathcal{T}}_{A}\left[\mathrm{H}\mathrm{e}\mathrm{r}\mathrm{m}\left({\mathcal{H}}_{2}\right)\right]$, is a plane that contains the disk

$$\begin{equation*}{D}_{A}^{ \overrightarrow {p}}=\left\{z={p}_{2}-{p}_{1}\vert z=\mathrm{Tr}\left(\rho {\sigma }_{z}\right),\enspace {p}_{i}=\mathrm{Tr}\left(\rho {E}_{i}\right),\enspace \rho {\geqslant}0,\mathrm{Tr}\left(\rho \right)=1\right\},\end{equation*}$$

i.e. the disk contains the full set of generator states ρ_gen. Note that ${\mathcal{T}}_{A}$ is not a completely positive trace preserving (CPTP) map, as ${\mathcal{T}}_{A}\left[\mathrm{H}\mathrm{e}\mathrm{r}\mathrm{m}\left({\mathcal{H}}_{2}\right)\right]$ extends beyond the disk ${D}_{A}^{ \overrightarrow {p}}$, i.e. outside the space of states. Indeed, for any state ρ that is not a convex combination of projectors E_i , there exists a probability distribution $\overrightarrow {p}$ such that ${\mathcal{T}}_{A}\left(\rho \right)$ is not positive semi-definite. Roughly speaking, any point inside the Bloch sphere from figure 1 but outside the blue vertical line is projected by ${\mathcal{T}}_{A}$ outside the sphere, for a sufficiently small disk ${D}_{A}^{ \overrightarrow {p}}$.

Zoom In Zoom Out Reset image size

In the practice of quantum state tomography, one collects a set of probability distributions $\overrightarrow {{p}_{1}},\dots , \overrightarrow {{p}_{\ell }}$ from a set of ℓ POVM measurements ${A}_{1}={\left\{{E}_{i}^{\left(1\right)}\right\}}_{i{\leqslant}{m}_{1}},\dots ,{A}_{\ell }={\left\{{E}_{i}^{\left(\ell \right)}\right\}}_{i{\leqslant}{m}_{\ell }}$, implemented over an ensemble of physical systems identically prepared in a quantum state ρ_gen. The statistics collected allows a unique state reconstruction when considering an informationally-complete (IC) sets of observables A₁, ..., A_ℓ . Our algorithm for quantum state estimation, algorithm 1 below, defines a sequence of Hermitian operators ρ_n , not necessarily composed by quantum states, that converges to the unique quantum state that is solution to the tomography problem, i.e. ρ_gen. For the moment, we assume error-free state tomography in our statements. The algorithm applies to any finite dimensional Hilbert space ${\mathcal{H}}_{d}$, and any informationally complete set of quantum observables.

Algorithm 1. Quantum state estimation algorithm.

Input: dimension $\mathrm{d}\in \mathbb{N}$, POVMs A1, ..., Aℓ acting on ${\mathcal{H}}_{d}$,

experimental frequencies ${ \overrightarrow {f}}_{1},\dots ,{ \overrightarrow {f}}_{\ell }\in {\mathbb{R}}^{m}$ and accuracy ∈ [0, 1].

Output: estimate ${\rho }_{\text{est}}\in \mathcal{B}\left({\mathcal{H}}_{d}\right)$.

${\rho }_{0}=\mathbb{I}/d$

$\rho ={\mathcal{T}}_{{A}_{\ell }}{\circ}\dots {\circ}{\mathcal{T}}_{{A}_{1}}\left({\rho }_{0}\right)$

repeat

ρold = ρ

$\rho ={\mathcal{T}}_{{A}_{\ell }}{\circ}\dots {\circ}{\mathcal{T}}_{{A}_{1}}\left({\rho }_{\text{old}}\right)$

Until $\mathfrak{D}\left(\rho ,{\rho }_{\text{old}}\right){\leqslant}{\epsilon}$

Return $\text{arg}\enspace {\text{min}}_{{\rho }_{\text{est}}\in \mathcal{D}\left({\mathcal{H}}_{\mathrm{d}}\right)}\mathfrak{D}\left(\rho ,{\rho }_{\text{est}}\right)$

In algorithm 1, $\mathcal{D}\left({\mathcal{H}}_{\mathrm{d}}\right)$ denotes the set of density operators over ${\mathcal{H}}_{\mathrm{d}}$. Theorem 3.1 below asserts the convergence of algorithm 1 when the input frequencies are exact, i.e. Born-rule, probabilities of an IC set of POVMs.

Theorem 3.1. Let A₁, ..., A_ℓ be a set of informationally complete POVMs acting on a Hilbert space ${\mathcal{H}}_{d}$, associated to a compatible set of probability distributions $\overrightarrow {{p}_{1}},\dots , \overrightarrow {{p}_{\ell }}$. Then, algorithm 1 converges to the unique solution to the quantum state tomography problem.

Here, compatibility refers to the existence of a quantum state associated to exact probability distributions $\overrightarrow {{p}_{1}},\dots , \overrightarrow {{p}_{\ell }}$ what is guaranteed when probabilities come from a generator state ρ_gen. Theorem 3.1 asserts that the composite map ${\mathcal{T}}_{{A}_{\ell }}{\circ}\dots {\circ}{\mathcal{T}}_{{A}_{1}}$ defines a dynamical system having a unique attractive fixed point. The successive iterations of algorithm 1 define a Picard sequence [31]:

$$\begin{align}\hfill {\rho }_{0}& =\mathbb{I}/d,\hfill \\ \hfill {\rho }_{n}& ={\mathcal{T}}_{{A}_{\ell }}{\circ}\dots {\circ}{\mathcal{T}}_{{A}_{1}}\left({\rho }_{n-1}\right),\quad n{\geqslant}1.\hfill \end{align} \tag{ 5 }$$

Note that for arbitrary chosen set of observables, the composition of physical imposition operators depends on its ordering. According to theorem 3.1, this ordering does not affect the success of the convergence in infinitely many steps. However, in practice one is restricted to a finite sequence, where different orderings produce different quantum states as an output. Nonetheless, such difference tends to zero when the state ρ_n is close to the attractive fixed point, i.e. solution to the state tomography problem. According to our experience from numerical simulations, we did not find any advantage from considering a special ordering for composition of operators.

Figure 2 shows the convergence of ρ_n in the Bloch sphere representation for a single qubit system and three PVMs taken at random. For certain families of measurements, e.g. MUBs and tensor product of Pauli matrices, the resulting Picard sequences and, therefore, algorithm 1 converge in a single iteration, see proposition 3.2. That is, ρ_n = ρ₁ for every n ⩾ 1. We numerically observed this same behavior for the 3^N product Pauli eigenbases in the space of N-qubits, with 1 ⩽ N ⩽ 8, conjecturing that it holds for every $N\in \mathbb{N}$, see section 4.2.

Zoom In Zoom Out Reset image size

In a previous work [32], a related algorithm was introduced for quantum state estimation. However, it has several disadvantages with respect to our work, namely: (i) it works for pure states only; (ii) the dynamics is non-linear, requiring a large runtime to converge (iii) convergence to the target state is not guaranteed. The main reason behind this last property is the existence of a large amount of undesired basins of attraction, as the solution to the problem is not the only attractive fixed point; finally, (iv) realistic state reconstruction is not possible due to the impossibility to introduce realistic noise, as it destroys purity. Note that algorithm 1 does not reduce to the one defined in reference [32] when reconstructing pure states, as our imposition operator is linear.

3.1. Ultra-fast convergence

When considering maximal sets of MUBs, the Picard sequences featuring in algorithm 1 converge in a single iteration. This is so because the associated imposition operators commute for MUB. This single-iteration convergence is easy to visualize in the Bloch sphere for a qubit system, as the three disks associated to three MUB are mutually orthogonal, and orthogonal projections acting over orthogonal planes keep the impositions within the intersection of the disks. The same argument also holds in every dimension. Let us formalize this result.

Proposition 3.1. Let T_A and T_B be two physical imposition operators associated to two MUBs A and B. Therefore,

$$\begin{equation}{T}_{B}{\circ}{T}_{A}={T}_{A}{\circ}{T}_{B}={T}_{A}+{T}_{B}-\mathbb{I}.\end{equation} \tag{ 6 }$$

In particular, note that T_A and T_B commute.

Also, it is easy to see from item (b), proposition 2.1 that operators ${T}_{{E}_{i}}$ commute when considering E_i equal to the tensor product local Pauli group. In this case, operators E_i do not form a POVM but given that they define an orthogonal basis in the matrix space, they are an informationaly complete set of observables. Let us now show the main result of this section:

Proposition 3.2. Algorithm 1 converges in a single iteration to the unique solution of the quantum state tomography problem for product of generalized Pauli operators and also for d + 1 MUBs, in any prime power dimension d.

We observe from simulations that the speedup predicted by proposition 3.2 has no consequences in the reconstruction fidelity of our method, which is actually higher than the one provided by MLE.

Theoretical developments from sections 2 and 3 apply to the ideal case of error free probabilities coming from an exact generator state ρ_gen. In practice, probabilities are estimated from frequencies, carrying errors due to finite statistics. Moreover, the states being prepared in each repetition of the experiment are affected by unavoidable systematic errors. These sources of errors imply that the output of algorithm 1 is typically outside the set of quantum states when considering experimental data. We cope with this situation by finding the closest quantum state to the output, called ρ_est in Hilbert–Schmidt (a.k.a. Frobenius) distance, for which there are closed-form expressions [17]. In the following, we provide numerical evidence for robustness of our method in the finite-statistics regime with white noise affecting the generator states, i.e. errors at the preparation stage. That is, we consider noisy states of the form $\tilde {\rho }\left(\lambda \right)=\left(1-\lambda \right)\rho +\lambda \mathbb{I}/d$, where λ quantifies the amount of errors. We understand there are more sophisticated techniques to consider errors, e.g. ill-conditioned measurement matrices [19]. Nonetheless, we believe the consideration of another model to simulate a small amount of errors would not substantially change the exhibited results. We reconstruct the state for N-qubit systems with 1 ⩽ N ⩽ 8, by considering the following sets of measurements: (a) MUBs, (b) tensor product of local Pauli bases and (c) a set of d + 1 informationally complete bases taken at random with Haar distribution. The last case does not have a physical relevance but illustrates performance of our algorithm for a set of measurements defined in an unbiased way. As a benchmark, we compare the performance of our method with the conjugate gradient, accelerated-gradient-descent (CG-AGP) implementation of maximum likelihood estimation (MLE) [18]. Computations were conducted on an Intel core i5-8265U laptop with 8 gb RAM. For the CG-AGP algorithm, we used the implementation provided by authors of reference [18], see reference [33]. We provide an implementation of our algorithm 1 in Python [34], together with the code to run the simulations presented in the current section.

4.1. Mutually unbiased bases

Figure 3 shows performance of algorithm 1 in the reconstruction of N-qubit density matrices from the statistics of a maximal set of 2^N + 1 MUBs. We consider a generator state ρ_gen in dimension d, taken at random according to the Haar measure distribution, with the addition of a 10% level of white noise, i.e. $\tilde {\rho }\left(\lambda \right)=\left(1-\lambda \right)\rho +\lambda \mathbb{I}/{2}^{N}$, with λ = 0.1. Here, it is important to remark that fidelities are compared with respect to the generator state ρ_gen, so that the additional white noise reflects the presence of systematic errors in the state preparation process. Probabilities are estimated from frequencies, i.e. ${f}_{j}={\mathcal{N}}_{j}/\mathcal{N}$ with ${\mathcal{N}}_{j}$ the number of counts for outcome j of some POVM and $\mathcal{N}={\sum }_{j}{\mathcal{N}}_{j}$ the total number of counts. Our simulations consider $\mathcal{N}=100{\times}{2}^{N}$ samples per measurement basis. Our figure of merit is the fidelity $F\left({\rho }_{n},{\rho }_{\text{gen}}\right)=\mathrm{Tr}{\sqrt{\sqrt{{\rho }_{\text{gen}}}{\rho }_{n}\sqrt{{\rho }_{\text{gen}}}}}^{2}$ between the reconstructed state after n iterations ρ_n and the generator state ρ_gen. Runtime of the algorithm is averaged over 50 independent runs, each of them considering a generator state ρ_gen chosen at random according to the Haar measure.

Zoom In Zoom Out Reset image size

4.2.  N-qubit Pauli bases

Here, we consider the reconstruction of N-qubit density matrices from the 3^N PVMs determined by all the products of single qubit Pauli eigenbases, for N = 1, ..., 8. Similarly to the case of MUBs, Picard sequences ${\rho }_{n}={T}_{\text{Pauli}}^{n}\left({\rho }_{0}\right)$ converge in a single iteration when product of Pauli measurements are considered, for any generator state ρ_gen and any initial state ρ₀. Figure 4 shows performance of a single iteration of these Picard sequences, where the generator state ρ_gen is taken at random, according to the Haar measure. Algorithm CG-AGP exploits the product structure of the N-qubit Pauli bases to speedup its most computationally expensive part: the computation of the probabilities given by the successive estimates in the MLE optimization. It does so by working with reduced density matrices which, in turn, imply an efficient use of memory. In order to have a fair comparison with our method, we decided to include the time to compute the N-qubit observables from the single Pauli observables in the total runtime of our algorithm. In practice, however, one would preload them into memory, as they are, of course, not a function of the input, i.e. of the observed probabilities. Nonetheless, figure 4 shows that our algorithm 1 has better fidelities with respect to the algorithm provided in reference [18].

Zoom In Zoom Out Reset image size

4.3. Random measurements for N-qubit systems

The simulations in the preceding subsections correspond to informationally complete sets of measurements for which algorithm 1 converges in a single iteration. To test whether the advantage over [18] hinges critically on this fact, we have numerically tested our algorithm with sets of PVMs selected at random, with respecto to the Haar measure. In figure 5 we show that in this case, the advantage in fidelity increases substantially, compared to figures 3 and 4.

Zoom In Zoom Out Reset image size

Finally, we would like to mention the projective least squares (PLS) quantum state reconstruction [22]. This method outperforms both in runtime and fidelity our algorithm 1. This occurs when the linear inversion procedure required by the method is not solved but taken from analytically existing reconstruction formula. Existing inversion formulas are known for complex projective two-designs, measurement composed by stabilizer states, Pauli observables and uniform/covariant POVM, see [22]. However, when taking into account the cost of solving the linear inversion procedure, our method has a considerable advantage over PLS. For instance, PLS does not have such efficient speed up for a number of physically relevant observables for which there is no explicit inversion known, including the following cases: (a) discrete Wigner functions reconstruction for arbitrary dimensional boson and fermions quantum systems from discrete quadratures, that be treated as observables by considering Ramsey techniques [35], (b) reconstruction of single quantized cavity mode from magnetic dipole measurements with Stern–Gerlach apparatus [36], (c) minimal state reconstruction of d-dimensional quantum systems from POVM consisting on d² elements, inequivalent to SIC-POVM [37], (d) spin s density matrix state reconstruction from Stern–Gerlach measurements [38], (e) quantum state tomography for multiparticle spin 1/2 systems [39], neither reduced to MUBs nor local Pauli measurements.

We introduced an iterative method for quantum state estimation of density matrices from any informationally complete set of quantum measurements in any finite dimensional Hilbert space. We demonstrated convergence to the unique solution for any informationally complete or overcomplete set of POVMs, see theorem 3.1. The method, based on dynamical systems theory, exhibited a simple and intuitive geometrical interpretation in the Bloch sphere for a single qubit system, see figures 1 and 2. Our algorithm revealed an ultra-fast convergence for a wide class of measurements, including MUBs and tensor product of generalized Pauli observables for an arbitrary large number of particles having d internal levels. These results improved fidelities reported by the CG-AGP super-fast MLE estimation [18] for all the studied cases, see section 3.1. Furthermore, numerical simulations revealed strong robustness under the presence of realistic errors in both state preparation and measurement stages, see figures 3–5. We provided an easy to use code developed in Python to implement our algorithm, see [34].

As interesting future lines of research, we pose the following list of open issues: (i) find an upper bound for fidelity reconstruction of algorithm 1 as a function of errors and number of iterations; (ii) characterize the full set of quantum measurements for which algorithm 1 converges in a single iteration; (iii) extend our method to quantum process tomography.

All data that support the findings of this study are included within the article and the computational code provided in Ref. [34].

It is a pleasure to thank Gustavo Cañas Cárdona, Zdenek Hradil, Felix Huber, Santiago Gómez López, Kamil Korzekwa, Andrew Scott, Oliver Reardon-Smith, Stephen Walborn, Andreas Winter and Karol Życzkowski for valuable comments. DG and DU are supported by Grant FONDECYT Iniciación Number 11180474, Chile. DU also acknowledges support from Project ANT1956, Universidad de Antofagasta, Chile. GS acknowledges support from the Government of Spain (FIS2020-TRANQI and Severo Ochoa CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381) and the EU project QRANGE. This work was supported by MINEDUC-UA project, code ANT 1856, which allowed a visit of GS to University of Antofagasta, Chile. As regards to the authorship of the different sections, DUC and DG provided both the theoretical background as well as the new mathematical results, whereas DUC and GS contributed with numerical simulations.

In this section we provide the proofs of all our results.

A.1. Algorithm for quantum state estimation

Proposition 2.1.  The following properties hold for any POVM ${\left\{{E}_{i}\right\}}_{i{\leqslant}m}$ and any ρ acting on ${\mathcal{H}}_{d}$:

• (a)  
  Imposition of physical information: $\mathrm{Tr}\left[{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right){E}_{i}\right]={p}_{i}$.
• (b)  
  Composition: ${T}_{{E}_{j}}^{{p}_{j}}{\circ}{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)={T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)+{T}_{{E}_{j}}^{{p}_{j}}\left(\rho \right)-\rho -\left({p}_{i}-\mathrm{Tr}\left(\rho {E}_{i}\right)\right)\mathrm{Tr}\left({E}_{i}{E}_{j}\right){E}_{j}/\mathrm{Tr}{\left({E}_{j}\right)}^{2}$.
• (c)  
  Non-expansiveness: $\mathfrak{D}\left({T}_{{E}_{j}}^{{p}_{j}}\left(\rho \right),{T}_{{E}_{j}}^{{p}_{j}}\left(\sigma \right)\right){\leqslant}\mathfrak{D}\left(\rho ,\sigma \right)$.

Proof. Items (a) and (b) easily arise from definition 2.1. In order to show the non-expansiveness stated in item (c), let us apply definition 1 to two states ρ and σ, belonging to ${\mathcal{H}}_{d}$, i.e.

$$\begin{equation}{T}_{{E}_{i}}^{{p}_{i}}\left(\rho \right)=\rho +\frac{\left({p}_{i}-\mathrm{Tr}\left[\rho {E}_{i}\right]\right){E}_{i}}{\mathrm{Tr}\left({E}_{i}^{2}\right)},\end{equation} \tag{ A1 }$$

$$\begin{equation}{T}_{{E}_{i}}^{{p}_{i}}\left(\sigma \right)=\sigma +\frac{\left({p}_{i}-\mathrm{Tr}\left[\sigma {E}_{i}\right]\right){E}_{i}}{\mathrm{Tr}\left({E}_{i}^{2}\right)}.\end{equation} \tag{ A2 }$$

Subtracting (A1) from (A2)

$$\begin{equation}{T}_{{E}_{i}}\left(\rho \right)-{T}_{{E}_{i}}\left(\sigma \right)=\left(\rho -\sigma \right)-\frac{\mathrm{Tr}\left[\left(\rho -\sigma \right){E}_{i}\right]{E}_{i}}{\mathrm{Tr}\left({E}_{i}^{2}\right)},\end{equation} \tag{ A3 }$$

where we dropped the upper index p_i from ${T}_{{E}_{i}}^{{p}_{i}}$. Now, let us compute

$$\begin{equation*}\mathfrak{D}{\left({T}_{{E}_{j}}\left(\rho \right),{T}_{{E}_{j}}\left(\sigma \right)\right)}^{2}=\mathrm{Tr}\left[\left({T}_{{E}_{i}}\left(\rho \right)-{T}_{{E}_{i}}\left(\sigma \right)\right){\left({T}_{{E}_{i}}\left(\rho \right)-{T}_{{E}_{i}}\left(\sigma \right)\right)}^{{\dagger}}\right].\end{equation*}$$

Thus,

$$\begin{align}\hfill \mathfrak{D}{\left({T}_{{E}_{j}}\left(\rho \right),{T}_{{E}_{j}}\left(\sigma \right)\right)}^{2}& =\mathfrak{D}{\left(\rho ,\sigma \right)}^{2}-2\frac{\mathrm{Tr}\left[\left(\rho -\sigma \right){E}_{i}\right]\mathrm{Tr}\left[\left(\rho -\sigma \right){E}_{i}\right]}{\mathrm{Tr}\left({E}_{i}^{2}\right)}+\frac{{\left(\mathrm{Tr}\left[\left(\rho -\sigma \right){E}_{i}\right]\right)}^{2}\enspace \mathrm{Tr}\left({E}_{i}^{2}\right)}{{\left(\mathrm{Tr}\left({E}_{i}^{2}\right)\right)}^{2}}\hfill \\ \hfill & =\mathfrak{D}{\left(\rho ,\sigma \right)}^{2}-\frac{{\left(\mathrm{Tr}\left[\left(\rho -\sigma \right){E}_{i}\right]\right)}^{2}}{\mathrm{Tr}\left({E}_{i}^{2}\right)},\hfill \end{align} \tag{ A4 }$$

where $\mathfrak{D}{\left(\rho ,\sigma \right)}^{2}=\mathrm{Tr}\left[\left(\rho -\sigma \right){\left(\rho -\sigma \right)}^{{\dagger}}\right]$. Therefore, $\mathfrak{D}\left({T}_{{E}_{j}}\left(\rho \right),{T}_{{E}_{j}}\left(\sigma \right)\right){\leqslant}\mathfrak{D}\left(\rho ,\sigma \right)$ and item (c) holds. □

Theorem 3.1. Let A₁, ..., A_ℓ be a set of informationally complete POVMs acting on a Hilbert space ${\mathcal{H}}_{d}$, associated to a compatible set of probability distributions $\overrightarrow {{p}_{1}},\dots , \overrightarrow {{p}_{\ell }}$. Therefore, algorithm 1 converges to the unique solution to the quantum state tomography problem.

Proof. First, from item (a) in proposition 2.1 the generator state ρ_gen is a fixed point of each imposition operator ${\mathcal{T}}_{{A}_{i}}$, for every chosen POVM measurement A₁, ..., A_ℓ . Hence, ρ_gen is a fixed point of the composition of all involved operators. Moreover, this fixed point is unique, as there is no other quantum state having the same probability distributions for the considered measurements, as A₁, ..., A_ℓ are informationally complete. Here, we are assuming error-free probability distributions. Finally, convergence of our sequences is guaranteed by the alternating projections method developed by Halperin, which states that successive iterations of non-expansive projections converge to a common fixed point of the involved maps, see theorem 1 in [40]. □

A.2. Single-step convergence

Proposition 3.1. Let ${\mathcal{T}}_{A}$ and ${\mathcal{T}}_{B}$ be physical imposition operators associated to two mutually unbiased bases A and B, for n qudit systems. Therefore

$$\begin{equation}{\mathcal{T}}_{A}{\circ}{\mathcal{T}}_{B}={\mathcal{T}}_{A}+{\mathcal{T}}_{B}-\mathbb{I}.\end{equation} \tag{ A5 }$$

In particular, notice that ${\mathcal{T}}_{A}$ and ${\mathcal{T}}_{B}$ commute.

Proof. First, it is simple to show that ${\mathcal{T}}_{A}\left(\rho \right)={\rho }_{0}+{\sum }_{j=0}^{{m}_{A}-1}{{\Pi}}_{j}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}$ for any PVM A, where Π_j = E_j are the subnormalized rank-one PVM elements. Thus, we have

$$\begin{align*}\hfill {T}_{B}{\circ}{T}_{A}\left({\rho }_{0}\right)& ={\rho }_{0}+\sum\limits _{j=0}^{{m}_{A}-1}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}\hfill \\ \hfill & \quad +\sum\limits _{k=0}^{{m}_{B}-1}{{\Pi}}_{k}^{B}\left[\rho -\left({\rho }_{0}+\sum\limits _{j=0}^{{m}_{A}-1}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}\right)\right]{{\Pi}}_{k}^{B}\hfill \\ \hfill & ={\rho }_{0}+\sum\limits _{j=0}^{{m}_{A}-1}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}+\sum\limits _{k=0}^{{m}_{B}-1}{{\Pi}}_{k}^{B}\left(\rho -{\rho }_{0}\right){{\Pi}}_{k}^{B}+\sum\limits _{j,k}{{\Pi}}_{k}^{B}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}{{\Pi}}_{k}^{B}\hfill \end{align*}$$

On the other hand,

$$\begin{align*}\hfill \sum\limits _{j,k}{{\Pi}}_{k}^{B}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}{{\Pi}}_{k}^{B}& =\sum\limits _{j,k}\enspace \mathrm{Tr}\left({{\Pi}}_{j}^{A}{{\Pi}}_{k}^{B}\right)\mathrm{Tr}\left(\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}\right){{\Pi}}_{k}^{B}\hfill \\ \hfill & =\gamma \left(A,B\right)\sum\limits _{j,k}\enspace \mathrm{Tr}\left(\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}\right){{\Pi}}_{k}^{B}\hfill \\ \hfill & =\gamma \left(A,B\right)\mathrm{Tr}\left(\rho -{\rho }_{0}\right)\hfill \\ \hfill & =0.\hfill \end{align*}$$

Therefore, we have

$$\begin{equation}{T}_{B}{\circ}{T}_{A}\left({\rho }_{0}\right)={\rho }_{0}+\sum\limits _{j=0}^{{m}_{A}-1}{{\Pi}}_{j}^{A}\left(\rho -{\rho }_{0}\right){{\Pi}}_{j}^{A}+\sum\limits _{k=0}^{{m}_{B}-1}{{\Pi}}_{k}^{B}\left(\rho -{\rho }_{0}\right){{\Pi}}_{k}^{B}\end{equation} \tag{ A6 }$$

$$\begin{equation}={T}_{A}\left({\rho }_{0}\right)+{T}_{B}\left({\rho }_{0}\right)-{\rho }_{0},\end{equation} \tag{ A7 }$$

for any initial state ρ₀. So, we have ${T}_{B}{\circ}{T}_{A}={T}_{A}{\circ}{T}_{B}={T}_{A}+{T}_{B}-\mathbb{I}$. □

Proposition 3.2.  Algorithm 1 converges in a single iteration to the unique solution of the quantum state tomography problem for product of generalized Pauli operators and also for d + 1 MUBs, in any prime power dimension d.

Proof. For generalized Pauli operators, commutativity of imposition operators comes from orthogonality condition Tr(E_i E_j ), see item (b) in proposition 2.1. Thus, we have

$$\begin{align}\hfill {\rho }_{n}& ={\left({T}_{{E}_{{d}^{2}}}{\circ}\dots {\circ}{T}_{{E}_{1}}\right)}^{n}\left({\rho }_{0}\right)\hfill \\ \hfill & ={T}_{{E}_{{d}^{2}}}^{n}{\circ}\dots {\circ}{T}_{{E}_{1}}^{n}\left({\rho }_{0}\right)\hfill \\ \hfill & ={T}_{{E}_{{d}^{2}}}{\circ}\dots {\circ}{T}_{{E}_{1}}\left({\rho }_{0}\right),\hfill \end{align} \tag{ A8 }$$

where the second step considers commutativity and the last step the fact that every T_j , j = 1, ..., d + 1 is a projection. On the other hand, from theorem 3.1 we know that ρ_n → ρ_gen when n → ∞, for any generator state ρ_gen. From combining this result with (A8) we have

$$\begin{equation}{T}_{{E}_{{d}^{2}}}{\circ}\dots {\circ}{T}_{{E}_{1}}\left({\rho }_{0}\right)={\rho }_{\text{gen}},\end{equation} \tag{ A9 }$$

for any seed ρ₀ and any generator state ρ_gen, in any prime power dimension d.

For MUB the result holds in the same way, where commutativity between the associated imposition operators associated to every PVM arises from see proposition 3.1. □

Along the work, we implemented simulations considering errors in both state preparation and those arising from finite statistics. In this section, we consider an additional source of errors in the measurement process. Specifically, we consider errors in the measurement apparatus, which is modeled by adding Gaussian perturbations in the direction of spin observables. In figure 6, we show fidelity for quantum state reconstruction for a spin 1/2 particle from three spin observables along orthogonal directions. For the Gaussian noise model, such directions are affected by a Gaussian probability distribution having standard deviation ν, centered in the ideally expected direction. That is, we consider the Gaussian probability distribution $p\left(x\right)\propto {\mathrm{e}}^{-{\left(x-\mu \right)}^{2}/2{\nu }^{2}}$ with μ = 0, for entries of a spin direction n, associated to the observable $S= \overrightarrow {n}\cdot \overrightarrow {\sigma }$, where $\overrightarrow {\sigma }=\left({\sigma }_{x},{\sigma }_{y},{\sigma }_{z}\right)$ is a vector composed by the three Pauli matrices. The amplitude of fluctuations can be controlled by adjusting the standard deviation ν.

Zoom In Zoom Out Reset image size