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.

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][33][34][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 C ≤ I −1 , the quantum Cramér-Rao inequality 36,37 C ≤ J −1 , and the Gill-Massar inequality 29 Tr(IJ −1 ) ≤ d − 1 , where d, C , I , and J 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 I(|ψ�, |ψ�) = 1 − |�ψ|ψ�| 2 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 p(|ψ�) 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 2(d − 1) independent parameters, the ultimate quantum estimation accuracy is given by the inequality Ī (|ψ�) ≥Ī p , with Ī p = (d − 1)/N 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 Ī p . 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, I = 0 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 f (z, z * ) : C n × C n → R the Wirtinger derivatives are defined by 47 where x i and y i are the real and imaginary parts of z i , respectively. These derivatives exist even if f (z, z * ) is nonholomorphic. Extremal points of f (z, z * ) are completely characterized by the conditions ∂ z * i f = 0 ∀ i = 1, . . . , d or, equivalently, ∂ z i f = 0 ∀ i = 1, . . . , d [48][49][50] . Thereby, the complex gradient is defined by g = ∂ z * f with ∂ z * = (∂ z * 1 , . . . , ∂ z * d ) . The CSPSA method is defined by the iterative rule 30 where a k is a positive gain coefficient and ẑ k is the estimate of the minimizer z of f (z, z * ) at the k-th iteration. The iteration starts from an initial guess ẑ 0 , which is randomly chosen. Instead of employing the complex gradient Scientific RepoRtS | (2020) 10:12781 | https://doi.org/10.1038/s41598-020-69646-z www.nature.com/scientificreports/ g , CSPSA resorts to an estimator ĝ k (ẑ k ,ẑ * k ) for the Wirtinger gradient g of f (z, z * ) the components of which are defined by with ẑ k± =ẑ k ± c k k , c k a positive gain coefficient, and ǫ k,± describes the presence of noise in the values of f (ẑ k± ,ẑ * k± ) . The components of the vector k ∈ C n are identically and independently distributed random variables in the set {±1, ±i} . The gain coefficients a k and c k control the convergence of CSPSA and are chosen as where the values of a, A, s, b and r can be adjusted to increase the rate of convergence. The estimates ẑ k provided by CSPSA converge asymptotically in mean to the minimizer of z of f and ĝ k is an asymptotically unbiased estimator of the Wirtinger gradient.
The use of an estimator ĝ k (ẑ k ,ẑ * k ) for the Wirtinger gradient g allows CSPSA to optimize functions with unknown parameters, provided that the values of f (ẑ k± ,ẑ * k± ) 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 |ψ(z)� of a qudit given by

the infidelity becomes
This function quantifies the deviation of the estimate |ψ(ẑ)� from the true unknwon state |ψ(z)� . In this function ẑ is a set of complex variables and z plays the role of a set of fixed unknown complex parameters. According to Eq. (4) CSPSA evaluates at each iteration k the infidelity I(|ψ(z)�), |ψ(ẑ)� at ẑ =ẑ k± . These values can be obtained by projecting the unknown state |ψ(z)� onto a d-dimensional orthonormal base B k± = {|ψ i,k± �} (with i = 0, . . . , d − 1 ) that contains the state |ψ(ẑ k± )� . This procedure generates the detection statistics n i,k± that lead to the probability distributions p i,k± = n i,k± / j n j,k± , where N est = j n j,k± is the total number of copies of the unknown states employed in the projective measurements. These probability distributions are employed to obtain the estimates for the infidelity, where we have assumed the convention |ψ 0,k± � = |ψ(ẑ k± )� . These estimates together with Eq. (3) generate the next estimate |ψ(ẑ k+1 )� for the unknown quantum state |ψ(z)� . 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. www.nature.com/scientificreports/ 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 2d − 2 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 {|a i ��a i |} has a likelihood function L(ρ) given by where the detection statistic of each projector is given by the number of counts n i and the total number of counts is given by N = i n i . If the quantum state prior to the measurement is ρ , then L(ρ) is the total joint probability of registering data {n i } . MLE is defined by the convex optimization problem At this point we link together CSPSA and MLE. Employing the accumulated data {n i,m± } between iterations m = 1 until m = k we define the accumulated likelihood L k (ρ) by the expression which is maximized in the set of pure states employing as starting guess the estimate |ψ(ẑ k+1 )� 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 k is also randomly chosen and two measurements are performed on an ensemble of size N est , which leads to finite statistics effects. Thus, each time the CSPSA-MLE method is employed to estimate a fixed unknown state |ψ(z)� , a different estimate |ψ(ẑ)� is generated. In this scenario, the accuracy of the estimation procedure for a fixed unknown state |ψ(z)� is given by the expectation value Ī (|ψ(z)�) 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 d with 2 × 10 2 pure quantum states |ψ(z)� of a single  Figure 1 displays a log-log graphic of the mean infidelity Ī (|ψ(z)�) (stars), for two randomly chosen states |ψ(z)� in 2 , that is, for a single qubit, as a function of the number k of iterations and for several values of the ensemble size N est . This mean infidelity is estimated as with G = 500 and R = 20 . Within the first 10 iterations the mean infidelity Ī (|ψ(z)�) exhibits a fast decrease that is followed by an asymptotic linear behavior. The decrease of Ī (|ψ(z)�) becomes more pronounced as N est increases. After 10 iterations the CSPSA-MLE method leads to a mean infidelity Ī (|ψ(z)�) approximately equal to 10 −2 , 7 × 10 −4 , 5 × 10 −5 , 5 × 10 −6 and 5 × 10 −7 , for increasing N est . For the same states CSPSA without MLE yields after 10 iterations a mean infidelity Ī (|ψ(z)�) approximately equal to 2.1 × 10 −1 , 5.2 × 10 −2 , 4.1 × 10 −2 , 4.0 × 10 −2 and 3.9 × 10 −2 (see Fig. 1 in 30 ), for increasing N est . Thus, the concatenation of CSPSA to MLE yields a mean infidelity that is 10 −1 to 10 −5 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 N est 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 N est = 10 4 and after 10 iterations the CSPSA-MLE method delivers a mean infidelity of 5 × 10 −6 . A similar value can be achieved via estimation with CSPSA with N est = 10 4 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.
After 7 iterations and for N est = 10, 10 2 , 10 3 , 10 4 , 10 5 the mean infidelity Ī (|ψ(z)�) asymptotically approximates linear behavior. This resembles the Gill-Massar lower limit Ī p for the average infidelity reached in the estimation of pure states, where now N = 2kN est 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 www.nature.com/scientificreports/ states attains a mean infidelity Ī (|ψ(z)�) lower than Ī p . Dashed lines in Fig. 1 correspond to Ī p as function of k for several values of N est . Clearly, the CSPSA-MLE method delivers a mean infidelity of the same order of magnitude than Ī p , albeit slightly higher. However, the gap between Ī (|ψ(z)�) and Ī p tends to close asymptotically as k increases. Thus, the infidelity provided by the CSPSA-MLE method tends to approach the Gill-Massar lower bound Ī p with a convergence rate that increases with N est . Figure 1 also displays the median infidelity M (|ψ(z)�) of I(|ψ(ẑ)�, |ψ(z)�) as a function of k for both randomly chosen states ψ(ẑ) and for several values of N est . 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 M (|ψ(z)�) . In the lineal regime, the Gill-Massar lower bound Ī p 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 Ī p . Nevertheless, the mean infidelity Ī (|ψ(z)�) is above the median M (|ψ(z)�) . 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 d display a similar behavior. This is shown in Fig. 2 where mean Ī (stars), median M (dots) and interquartile range (shaded areas) of the mean infidelity Ī (|ψ(z)�) over all states in d as a function of k for several values of N est are depicted and compared to Ī p (dashed lines). The mean Ī (stars) correspond to the expectation value of Ī (|ψ(z)�) on the Hilbert space, that is, where f (|ψ�) is the probability density function of randomly and uniformly generating the unknown state |ψ� . The mean Ī is estimated as Figure 2 exhibits a fast decrease of Ī until approximately iteration 7, which is followed by asymptotic linear behavior. The median M 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 M 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 Ī (|ψ(z)�) of the states in d are concentrated in an extremely narrow interval around the mean Ī and the median M . Thus, all states in d  www.nature.com/scientificreports/ 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 Ī p . 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 Ī p allowed by the laws of quantum mechanics. Figure 2 also displays the behavior of Ī (squares) and median M (rombos) obtained with the use of CSPSA only, that is, without employing MLE to refine the guesses provided by CSPSA, for the case of N est = 10 5 . 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.
Once the CSPSA-MLE method enters in a lineal regime, after approximately 10 iterations in the inspected dimensions, delivers an estimation accuracy Ī (|ψ�) close to Ī p . Thereby, we can write or since N = 2N est k also Thus, 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 N est each one. Thereby, we can employ a small value of N est and a large number of iterations or a large value of N est 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 Ī m 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 Ī m = [(d + 1)/2] 2Ī p . Tomographic methods designed to estimate unknown mixed states cannot achieve a better accuracy than Ī m as long as they resort to separable measurements on the ensemble of equally prepared copies. This departs quadratically from Ī p 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 Ī m . 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 2d − 2 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 d = 2, 4, 8 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 www.nature.com/scientificreports/ 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.