Estimation of pure quantum states in high dimension at the limit of quantum accuracy through complex optimization and statistical inference.

Quantum tomography has become a key tool for the assessment of quantum states, processes, and devices. This drives the search for tomographic methods that achieve greater accuracy. In the case of mixed states of a single 2-dimensional quantum system adaptive methods have been recently introduced that achieve the theoretical accuracy limit deduced by Hayashi and Gill and Massar. However, accurate estimation of higher-dimensional quantum states remains poorly understood. This is mainly due to the existence of incompatible observables, which makes multiparameter estimation difficult. Here we present an adaptive tomographic method and show through numerical simulations that, after a few iterations, it is asymptotically approaching the fundamental Gill-Massar lower bound for the estimation accuracy of pure quantum states in high dimension. The method is based on a combination of stochastic optimization on the field of the complex numbers and statistical inference, exceeds the accuracy of any mixed-state tomographic method, and can be demonstrated with current experimental capabilities. The proposed method may lead to new developments in quantum metrology.

Introduction

The basic principles of quantum mechanics strongly limit the processes that can be performed in nature. The linearity of quantum operations prohibits the possibility of creating in a deterministic way perfect copies of unknown quantum states[1]. Information encoded in non-orthogonal quantum states can be perfectly recovered only at the expense of allowing a probability of error[2-4]. Measurements, a fundamental tool for understanding the natural world, are also limited. In remarkable difference with the classical realm, a measurement in a quantum system alters the quantum state of the system, which prevents further information gain. In this way, an ensemble of many equally prepared systems is necessary to obtain information about a quantum system and, consequently, quantum measurements become statistical in nature[5-7]. Furthermore, quantum mechanics establishes through the quantum Cramér–Rao bound[8-10] a fundamental limit in the best accuracy that a measurement or estimation process can achieve as a function of the ensemble size.

The accuracy limit also bounds the estimation of unknown quantum processes and states. These are central problems in the theory of quantum measurements[11] and play a key role in the control of quantum systems, the benchmarking of quantum technologies[12-14], and quantum metrology[15]. Today there is a large collection of methods[11,16-25] for estimating unknown quantum states, which are collectively known as quantum tomographic methods. These are based on the post-processing of data acquired through the measurement of a set of positive operator-valued measures in an ensemble of N identically, independently prepared copies of the unknown state to be estimated. Recently, there has been a great deal of research activity aimed at designing tomographic methods with increasing accuracy, in particular methods that saturate the quantum accuracy limit. From the theoretical and experimental point of view, it has been shown that a two-stage adaptive tomographic method achieves the quantum accuracy limit in the case of unknown mixed states of a single qubit[26,27]. Unfortunately, this result does not hold in the higher dimensional case[28].

Here, we report an adaptive tomographic method that asymptotically approaches the quantum accuracy limit in the important case of estimating unknown pure quantum states in high dimensions. This is an instance of multi-parameter estimation, where due to the information trade-off among incompatible observables the progress has been slow. In the case of pure states the quantum accuracy is limited by the Gill–Massar lower bound[29]. The method is based on the concatenation of the Complex simultaneous perturbation stochastic approximation (CSPSA), a recently proposed iterative stochastic optimization method on the field of the complex numbers[30], and Maximum likelihood estimation (MLE), a well known statistical inference method[31]. The method here proposed reaches the quantum accuracy limit after a small number of iterations, typically of the order of 8, for all inspected dimensions. Thereby, the method makes an optimal use of the ensemble size and surpasses the estimation accuracy of known methods for pure-state tomography[32-35]. Moreover, the method also surpasses the estimation accuracy of any tomographic method designed to estimate mixed states via separable measurements on the ensemble of equally prepared copies.

Method

The accuracy achieved in the estimation of the parameters of a quantum state can be studied by means of three fundamental inequalities. These are the Cramér–Rao inequality[36] , the quantum Cramér–Rao inequality[36,37] , and the Gill–Massar inequality[29] , where , , and are the dimension of the Hilbert space, the covariance matrix, the classical Fisher information matrix, and the quantum Fisher information matrix, respectively. These inequalities allows one to deduce lower bounds for several metrics of accuracy, such as infidelity or mean square error, considering the impact of the ensemble size. As accuracy metric we employ the infidelity[38,39] between an unknown state and its estimate . The estimation accuracy is given by the expectation or mean value of the infidelity with respect to all possible estimates of a fixed unknown state , that is,where is the probability density function of obtaining the estimate . Depending on the general characteristics of the estimation process, the fundamental inequalities lead to various lower bounds for the estimation accuracy. For the estimation of pure states, which have independent parameters, the ultimate quantum estimation accuracy is given by the inequality , with [29,40,41] the Gill–Massar lower bound. This state-independent bound can be used as a benchmark to assess tomographic methods.

Our main goal is designing a tomographic method for pure quantum states in high dimensions that achieves an estimation accuracy equal to . The design of the proposed tomographic method originates from observing[42-45] that an unknown and fixed pure quantum state can be characterized as the minimizer of the infidelity in the Hilbert space of the estimate , that is, when . Then, we can envision the use of an optimization method to iteratively drive a sequence of estimates toward decreasing values of the infidelity.

The choice of the optimization method requires certain consideration. Traditional optimization methods are based on the evaluation of higher derivatives of the function to be optimized. In the case of the infidelity this is not possible, since the derivatives depend on the state to be estimated, which is unknown. We also require a fast convergence rate from the optimization method. Furthermore, the implementation of the method should be as economical as possible from the point of view of computing and physical resources.

To deal with these situations we resort to the Complex simultaneous perturbation stochastic approximation. This method can optimize real valued functions of complex variables that also depend on unknown complex parameters. The infidelity is a non-holomorphic function, that is, it violates the Cauchy-Schwarz conditions, which are necessary for the existence of a complex derivative. In this case, the usual approach is to optimize with respect to a real parametrization of the complex variables entering in the infidelity, that is, the complex probability amplitudes. This approach has some unwanted side effects. The elements of the real-valued gradient of the infidelity are in general more convoluted than would be those of a complex gradient formed by first order derivatives with respect to the initial complex variables. Additionally, any inherent structures present in the complex derivatives of the infidelity, which could be exploited to enhance the performance of optimization methods, are unused[46]. The CSPSA method has been designed to avoid these unwanted features. The main tool in the formulation of CSPSA is the Wirtinger complex calculus. For a target function the Wirtinger derivatives are defined by[47]where and are the real and imaginary parts of , respectively. These derivatives exist even if is non-holomorphic. Extremal points of are completely characterized by the conditions or, equivalently, [48-50]. Thereby, the complex gradient is defined by with . The CSPSA method is defined by the iterative rule[30]where is a positive gain coefficient and is the estimate of the minimizer of at the k-th iteration. The iteration starts from an initial guess , which is randomly chosen. Instead of employing the complex gradient , CSPSA resorts to an estimator for the Wirtinger gradient of the components of which are defined bywith , a positive gain coefficient, and describes the presence of noise in the values of . The components of the vector are identically and independently distributed random variables in the set . The gain coefficients and control the convergence of CSPSA and are chosen aswhere the values of a, A, s, b and r can be adjusted to increase the rate of convergence. The estimates provided by CSPSA converge asymptotically in mean to the minimizer of of f and is an asymptotically unbiased estimator of the Wirtinger gradient.

The use of an estimator for the Wirtinger gradient allows CSPSA to optimize functions with unknown parameters, provided that the values of can be obtained and the unknown parameters remain constant along the optimization procedure. This is precisely the case of the infidelity. Considering an unknown pure quantum state of a qudit given bywith , and an estimate of given bywith , the infidelity becomesThis function quantifies the deviation of the estimate from the true unknwon state . In this function is a set of complex variables and plays the role of a set of fixed unknown complex parameters. According to Eq. (4) CSPSA evaluates at each iteration k the infidelity at . These values can be obtained by projecting the unknown state onto a d-dimensional orthonormal base (with ) that contains the state . This procedure generates the detection statistics that lead to the probability distributions , where is the total number of copies of the unknown states employed in the projective measurements. These probability distributions are employed to obtain the estimatesfor the infidelity, where we have assumed the convention . These estimates together with Eq. (3) generate the next estimate for the unknown quantum state . The CSPSA method can be understood as a generalization of the Simultaneous perturbation stochastic approach (SPSA)[51,52], a well known stochastic gradient-free optimization method working in the field of the real numbers. SPSA has been proposed[42] and experimentally demonstrated[44] as tomographic method for pure quantum states. In this context, it has been show[30] that CSPSA achieves a higher convergence rate than SPSA.

CSPSA generates 2 probability distributions at each iteration, that is, a total of 2d different probabilities. However, only two of them are employed to estimate the required values of the infidelity. The remaining probabilities are not occupied by the algorithm. Thus, CSPSA generates a large amount of accumulated data that is simply discarded. Here, we show that precisely this information can be used to improve the convergence rate of CSPSA in such a way that the estimation accuracy reaches the Gill–Massar lower bound for the estimation of pure quantum states. This is accomplished by resorting to Maximum likelihood estimation, a well known statistical inference method that is extensively employed as a post-processing stage in quantum tomographic methods. MLE is a method in statistical inference aimed at estimating unknown parameters of a population from observed data. The underlying idea is to choose as estimator the maximizer of the probability of obtaining the observed data[53,54]. MLE was introduced in quantum tomography as a post-processing method[31] to obtain physically acceptable quantum states. A quantum system that undergoes a measurement process described by the set of projectors has a likelihood function given bywhere the detection statistic of each projector is given by the number of counts and the total number of counts is given by . If the quantum state prior to the measurement is , then is the total joint probability of registering data . MLE is defined by the convex optimization problemAt this point we link together CSPSA and MLE. Employing the accumulated data between iterations until we define the accumulated likelihood by the expressionwhich is maximized in the set of pure states employing as starting guess the estimate provided by CSPSA. The refined estimate provided by MLE is then employed as starting guess for the next iteration with CSPSA. We refer to this procedure as the CSPSA-MLE tomographic method. The main steps of the CSPSA-MLE tomographic method have been summarized as pseudocode in Algorithm 1 above.

Results

For a fixed unknown state, the CSPSA-MLE method exhibits three sources of randomness. Since there is no a priori information about the unknown state, the initial guess is chosen according to a Haar-uniform distribution. At each iteration, the vector is also randomly chosen and two measurements are performed on an ensemble of size , which leads to finite statistics effects. Thus, each time the CSPSA-MLE method is employed to estimate a fixed unknown state , a different estimate is generated. In this scenario, the accuracy of the estimation procedure for a fixed unknown state is given by the expectation value of Eq. (1).

To study the performance of the CSPSA-MLE method several Monte Carlo experiments in the regime of a small number of iterations were carried out. A set with pure quantum states of a single d-dimensional quantum system (qudit), uniformly distributed on the unit hypersphere, was generated. Each state in was reconstructed via the CSPSA-MLE- method considering a number G of initial guesses, also uniformly distributed on the unit hypersphere, and R independent simulations for each fixed pair of unknown state and initial guess. At each iteration the values of were estimated considering a multinomial distribution on an ensemble of size . For each state in ; mean, variance, median and interquartile range for the infidelity as functions of the number of iterations k for several values of ensemble size were estimated. Similar numerical experiments were performed via CSPSA without MLE for comparison purposes. Since the optimization of the gains is a computational costly problem, we have resorted to the gains and . These lead to a high rate of convergence for CSPSA in the regime of few iterations. The resting coefficients have been set to , and for , correspondingly.

Figure 1 displays a log-log graphic of the mean infidelity (stars), for two randomly chosen states in , that is, for a single qubit, as a function of the number k of iterations and for several values of the ensemble size . This mean infidelity is estimated aswith and . Within the first 10 iterations the mean infidelity exhibits a fast decrease that is followed by an asymptotic linear behavior. The decrease of becomes more pronounced as increases. After 10 iterations the CSPSA-MLE method leads to a mean infidelity approximately equal to and , for increasing . For the same states CSPSA without MLE yields after 10 iterations a mean infidelity approximately equal to and (see Fig. 1 in[30]), for increasing . Thus, the concatenation of CSPSA to MLE yields a mean infidelity that is to times closer to the true minimum than the one provided by CSPSA alone. Let us note that this comparison considers the same type and amount of physical resources for both methods, that is, ensemble size for the estimation of the infidelity and total number 2dk of measurement outcomes. The estimation via CSPSA can reach similar mean infidelity values to the CSPSA-MLE method but at the expense of more iterations. For instance, with and after 10 iterations the CSPSA-MLE method delivers a mean infidelity of . A similar value can be achieved via estimation with CSPSA with after 100 iterations, which represents an increase of one order of magnitude in the total ensemble size N as well as in the total number of measurements. Thus, the concatenation of MLE to CSPSA provides a significative improvement in the rate of convergence and a large reduction of the required physical resources.

Figure 1

Mean and median infidelity (stars) and (circles), respectively, as functions of the number k of iterations achieved by CSPSA-MLE method for two randomly selected pure states of a single qubit in with with (violet), (green), (yellow), (light blue) and (purple) from top to bottom. Dashed lines indicate the Gill–Massar lower bound for . Shaded areas represent interquartile range. Values of gain coefficients are: , , , , and for , respectively

After 7 iterations and for the mean infidelity asymptotically approximates linear behavior. This resembles the Gill–Massar lower limit for the average infidelity reached in the estimation of pure states, where now is the total ensamble size employed after k iterations. This bound imposes a fundamental precision limit on the achievable mean infidelity: no method for the estimation of pure states attains a mean infidelity lower than . Dashed lines in Fig. 1 correspond to as function of k for several values of . Clearly, the CSPSA-MLE method delivers a mean infidelity of the same order of magnitude than , albeit slightly higher. However, the gap between and tends to close asymptotically as k increases. Thus, the infidelity provided by the CSPSA-MLE method tends to approach the Gill–Massar lower bound with a convergence rate that increases with .

Figure 1 also displays the median infidelity of as a function of k for both randomly chosen states and for several values of . This exhibits a much faster decrease and an earlier onset of the asymptotic linear behavior. The shaded areas in Fig. 1 correspond to the interquartile range, which is divided into two areas above and below the median infidelity . In the lineal regime, the Gill–Massar lower bound lays in the upper half of the interquartile range. This indicates that more than 50% of the reconstruction attempts lead to an infidelity lower than . Nevertheless, the mean infidelity is above the median . This points to the existence of a small fraction of realizations with high values of the infidelity, which increases the value of the mean infidelity above the median infidelity. As the number of iterations increases the impact of these realizations on the value of the mean infidelity decreases and mean and median infidelity tend to reach similar values.

The main features exhibited in Fig. 1 are typical, that is, all randomly chosen states in display a similar behavior. This is shown in Fig. 2 where mean (stars), median (dots) and interquartile range (shaded areas) of the mean infidelity over all states in as a function of k for several values of are depicted and compared to (dashed lines). The mean (stars) correspond to the expectation value of on the Hilbert space, that is,where is the probability density function of randomly and uniformly generating the unknown state . The mean is estimated asFigure 2 exhibits a fast decrease of until approximately iteration 7, which is followed by asymptotic linear behavior. The median shows a similar behavior, although it enters an asymptotic linear behavior faster than the mean in iteration 5 approximately. As the number k of iterations increases mean and median overlap almost perfectly in the linear regime. The interquartile range is very narrow and nearly indistinguishable from these two quantities. This indicates that in the linear regime the mean infidelities of the states in are concentrated in an extremely narrow interval around the mean and the median . Thus, all states in are estimated by the CSPSA-MLE method with an accuracy that is very close to . Furthermore, the estimation accuracy provided by the CSPSA-MLE method tends to converge from above to . Thereby, the CSPSA-MLE method produces a mean infidelity for any unknown pure state that approaches asymptotically the best possible estimation accuracy of pure states allowed by the laws of quantum mechanics. Figure 2 also displays the behavior of (squares) and median (rombos) obtained with the use of CSPSA only, that is, without employing MLE to refine the guesses provided by CSPSA, for the case of . Clearly, the CSPSA-MLE tomographic method outperforms the CSPSA tomographic method by at least 5 orders of magnitude. Therefore, the concatenation of MLE to CSPSA plays a key role in improving the accuracy of the estimation and bringing it closer to the lower bound of Gill–Massar.

Figure 2

Mean and median infidelity (stars, upper row) and (solid circles, lower row) as functions of the number k of iterations obtained via the CSPSA-MLE method for the estimation of single d-dimensional pure quantum states in for , from left to right, with (violet), (green), (yellow), (light blue) and (purple), from top to bottom. Shaded areas describe interquartile range. Dashed lines indicate the value of the Gill–Massar lower bound for . Mean (solid orange squares, upper row) and median (solid orange rombos, lower row) as functions of the number k of iterations obtained via the CSPSA method for the estimation of single d-dimensional pure quantum states in for with . Gain coefficients as in Fig. 1

Once the CSPSA-MLE method enters in a lineal regime, after approximately 10 iterations in the inspected dimensions, delivers an estimation accuracy close to . Thereby, we can writeor since alsoThus, in a given dimension d the CSPSA-MLE tomographic method achieves a predefined estimation accuracy with an ensamble size N that can be divided into 2k ensambles of size each one. Thereby, we can employ a small value of and a large number of iterations or a large value of and a small number of iterations. The last alternative provides a faster convergence rate. In this case, the estimation of the complex gradient is closer to the exact gradient and the convergence of CSPSA becomes similar to a complex formulation of a deterministic first-order iterative optimization algorithm. However, whether the first or second alternative is more appropriate also depends largely on the characteristics of the experimental platform where the estimation is realized.

We can also compare the mean infidelity achieved by the CSPSA-MLE method with the Gill–Massar lower bound for the mean infidelity achieved in the estimation of full-rank mixed states via separable measurements on the ensemble of equally prepared copies, which is given by . Tomographic methods designed to estimate unknown mixed states cannot achieve a better accuracy than as long as they resort to separable measurements on the ensemble of equally prepared copies. This departs quadratically from and as the dimensions increases. Thus, tomographic methods that do not employ the a priori information about the purity of the unknown state cannot estimate pure states with an accuracy better than . Thereby, the CSPSA-MLE method provides an advantage for pure state estimation over standard quantum tomography, two-stage standard quantum tomographic, and methods such as the ones based on mutually unbiased bases, symmetric informationally complete positive-operator-valued measures, and equidistant states.

Discussion

The accurate estimation of quantum states with a limited ensemble size is a difficult task. In the case of 2-dimensional pure quantum states, the best measurement strategy, that is, the one that leads to the (local unbiased) estimator that saturates the quantum Cramér–Rao bound, is generally a function of the parameters of the unknown state itself. Thereby, the use of the optimal estimation strategy is unfeasible[55]. It is possible, however, to employ an adaptive strategy that approaches the quantum Cramér–Rao bound. This strategy consists of a sequence of measurements where each one is optimal for a given guess of the unknown state[56,57]. The multi-parameter estimation of an unknown d-dimensional quantum state, which is defined by independent real numbers, cannot be carried out following a similar strategy since the optimal measurement strategy is unknown.

Instead, we have approached the estimation of unknown d-dimensional quantum states from the optimization of the metric used to characterize the accuracy of the estimation process. The optimization is solved by a method based on the concatenation of CSPSA to MLE. The proposed a method allows for estimating quantum pure states with high accuracy. CSPSA drives a sequence of projective measurements toward the infidelity minimizer. MLE provides at each iteration a refinement of the estimates considering the accumulated data generated by all previous measurements. Monte Carlo experiments for dimensions and 16 indicate that the mean infidelity for a fixed arbitrary unknown state exhibits a fast decrease within few iterations followed by an asymptotic linear trend. In the linear regime the attained mean infidelity closely approaches the fundamental limit on the accuracy established by the Gill–Massar lower bound for the mean infidelity of the estimation of pure states. Hence, the CSPSA-MLE method surpasses the accuracy of known tomographic methods for pure quantum states. The mean infidelity is also lower than the Gill–Massar lower bound for the infidelity of the estimation of mixed quantum states. Consequently, no tomographic method for mixed states can achieve a mean infidelity lower than the one attained by the CSPSA-MLE method. The median infidelity is also below the Gill–Massar lower bound for pure states. Therefore, more than 50% of the estimation attempts leads to lower infidelities than the Gill–Massar lower bound.

The CSPSA-MLE method exhibits a clear trade-off. The concatenation of CSPSA to MLE increases the rate of convergence of the mean infidelity, which leads to a decrease in the number of iterations and in the ensemble size. However, there is an increase in the computational complexity of the algorithm because in each iteration it is now necessary to solve the optimization problem corresponding to MLE. Recently, optimized methods for MLE in high dimension have been proposed[58,59].

An experimental realization of the CSPSA-MLE method applied to state of a single 2-dimensional quantum system can be carried out with current experimental techniques[44,57] for generating and measuring single-photon polarization states. The higher dimensional case can be demonstrated by means of experimental setups based on single photons and concatenated spatial light modulators[32,60,61] or via integrated quantum photonics[62,63]. These two experimental platforms offer the possibility of performing electronically controlled adaptive measurements.

Code availability

The source codes that support the results of this study are available from the corresponding author upon reasonable request.