Quantum measurement optimization by decomposition of measurements into extremals

Using the convex structure of positive operator value measurements and several quantities used in quantum metrology, such as quantum Fisher information or the quantum Van Trees information, we present an efficient numerical method to find the best strategy allowed by quantum mechanics to estimate a parameter. This method explores extremal measurements thus providing a significant advantage over previously used methods. We exemplify the method for different cost functions in a qubit and in a harmonic oscillator and find a strong numerical advantage when the desired target error is sufficiently small.


Quantum measurement optimization by decomposition of measurements into extremals esteban Martínez-Vargas 1 , carlos pineda 2 & pablo Barberis-Blostein 3 ✉
Using the convex structure of positive operator value measurements and several quantities used in quantum metrology, such as quantum fisher information or the quantum Van trees information, we present an efficient numerical method to find the best strategy allowed by quantum mechanics to estimate a parameter. This method explores extremal measurements thus providing a significant advantage over previously used methods. We exemplify the method for different cost functions in a qubit and in a harmonic oscillator and find a strong numerical advantage when the desired target error is sufficiently small.
One of the goals of metrology is to provide an optimal strategy to measure the value of a parameter under certain fixed conditions. For completeness, both an approximate value of the parameter and an estimation of the error must be given. If the physical system from which the parameter is to be estimated is analyzed within the framework of quantum mechanics, we shall speak about quantum metrology [1][2][3][4][5] . Motivated by the fact that quantum systems can offer an important advantage over classical systems in precision when estimating a parameter 6 , there has been intense theoretical and experimental advances in the area in the last years 7,8 . Moreover, precision in the estimation of parameters has several applications in the development of quantum technologies 9,10 and quantum state manipulation 11 .
To estimate a parameter of a physical system one acquires data through measurements; the estimation of the parameter is obtained by applying a function, known as the estimator, to the data. The probability distribution of measurement outcomes can be modelled using a statistical model of the experiment: the probability distribution of outcomes conditioned to the value of the parameter. This statistical model might describe, say, a noisy measurement apparatus. In this article we specialize to the case of quantum mechanics, where the probability distribution is given by the Born rule and the statistical model is obtained once it is decided which operator is going to be measured. In addition to the random component of measurement, we also consider classical noise, which we include through the density matrix formalism.
The cost function that quantifies the error of the parameter estimation of a quantum system (for example the mean square error) depends on both the measurement to be performed and the estimator; the optimal measurement strategy consists of the quantum measurement and the estimator that extremizes it. However, finding the extreme of a cost function over all possible quantum measurements and estimators is not simple. Alternatively, when a Cramér-Rao type inequality exists, the problem can be reduced to find the extreme of another cost function (for example the Fisher information) over all the quantum measurements. This simplifies the problem because it is not longer necessary to maximize over the space of estimators. One still has to deal with extremizing a cost function over all quantum measurements, which in general is difficult. However, under special circumstances, such as symmetries, the problem can be simplified 12 .
It is possible, though very costly, to numerically find the maximum over all quantum measurements of cost functions. The straightforward way to solve the problem is to randomly sample the space of all positive operator value measures (POVMs), evaluate the cost function on this sample, and keep the maximum value obtained. This method, that we call the random sampling method (RSM), is very inefficient since the POVM space is large.
In this paper, we show how to numerically find the maximum over all quantum measurements of the Fisher information and the Bayesian version of this bound, the Van Trees bound 13 , in a way that is orders of magnitude faster than using the RSM. We rely on the following: (i) the cost functions of interest in quantum metrology are convex with respect to the POVMs and (ii) the set of POVMs is convex 14 . Therefore, following the maximum principle 15,16 , the maximum over the POVMs must be on an extremal point of the POVM set. Our approach is simply to sample randomly extremal POVMs and get the highest value of the appropriate cost function. It is easy to produce efficiently a general POVM from a random unitary matrix, however, it is not trivial to produce randomly extremal POVMs. For this we used the algorithm proposed by Sentís et al. 17 .
We call our method the random extreme sampling method (RESM). The techniques presented here can be used, for example, to find numerically the maximum value of the quantum Van Trees information 18 , the quantum Fisher information (even if the input state is not pure), or any convex cost function. Together with the maximum, the corresponding quantum measurement is obtained. The method can also be applied to the case of multivariable convex cost functions, and thus used to find the optimal measurement in the estimation of several parameters. For example: Given a multivariable statistical model we can construct the Fisher matrix, its diagonal elements gives a bound on the variance of each parameter. We can use our method to find the quantum measurement that maximizes any convex cost function of the diagonal elements.
Our approach can also be used in other areas different from quantum metrology, where a maximization over POVMs is required. An example is statistical decision theory 19 . The problem is the following: we have a set of possible decisions, from which one is chosen depending on the outcome of a quantum system measurement. The distribution probability for the decisions is given by a quantum measurement. We look for the best decision. How good is the decision is rated by a cost function defined in the space of quantum measurements. The problem of finding the best decision is mapped to the problem of maximizing this cost function over all the quantum measurements. Note that quantum parameter estimation can be cast as a problem of statistical decision theory, where the decision to be chosen is an estimation of the parameter.

Results
In this section we first introduce some convex cost functions that give bounds to the error of parameter estimation. Then we present our main result: a numerical algorithm that allow us to find the maximum over all quantum measurement of any convex cost. We finish with examples of interest where we apply the algorithm.
Bounds in the error of parameter estimation. We discuss how to get bounds for the error made in an estimation process. We use two error measures: the mean squared error and the Bayesian mean squared error, which is used when some information is known about the parameter. This discussion is general and just assumes that a statistical model is given. We also discuss how to apply these ideas for a quantum system.

Crámer-Rao inequality.
In this section, we introduce some basic quantities needed to develop further the discussion. Let be the distribution probability of the outcomes y, of the random variable Y, conditioned to a fixed value of the real parameter θ. We assume that each y is a set of real numbers of fixed finite size. The function θ | p y ( ) is the statistical model; y represents what is measured in an experiment and its probability distribution, p, depends on the parameter θ. Using the outcomes y, the parameter is estimated through the real-valued function θ y ( ), known as the estimator. The estimator is unbiased if it is on average correct, meaning that its expected value is equal to the actual value of the parameter,ˆ∫ The uncertainty of the estimator is given by the mean squared error, defined as 2 2 We say that the measurement strategy is optimal if the estimator minimizes the mean squared error. Finally, let us define the Fisher information If the estimator is unbiased (i.e. if Eq. (2) holds), using the Cauchy-Schwarz inequality, one arrives to the Cramér-Rao inequality 20,21 : Note that ς 2 depends on the choice of the specific estimator θˆy ( ), whereas the Fisher information depends only on the distribution probability of the random variable. From the Cramér-Rao inequality, we see that the inverse of the Fisher information bounds from below the mean squared error independently of the estimator we use: the larger the Fisher information the smaller the error bound. Fisher showed that in the limit where the number of measurements goes to infinity, the maximum likelihood estimator saturates this inequality 22 .
Quantum Cramér-Rao inequality. We want to find the best measurement strategy to estimate a parameter, θ, that appears in the Hamiltonian of a quantum system. In order to estimate the parameter we proceed as follows: we start with an initial state and let the system evolve some fixed time. The dynamics of the system depends on the parameters of the Hamiltonian; after the evolution, the state of the system depends on the parameter we want to estimate: ρ ρ θ = ( ). We then measure some observable of the system and use the result to estimate θ. Fixing the Hamiltonian, time and initial state, we want to know if the strategy we are using minimizes the error in the parameter estimation.
General measurements in quantum mechanics are described by the positive operator valued measure (POVM) formalism, which we briefly recall in order to fix the notation. If E { ( )} ξ is a POVM parametrized by the real parameter ξ, for each value of ξ, ξ E( ) is a self-adjoint operator on the system Hilbert space. They satisfy the completeness relationÊ ∫ ξ ξ = and the probability of measuring the result ξ is In order for Eq. (7) to be a probability distribution we require that the elements E( ) ξ to be positive semidefinite, Notice that ξ can also belong to a finite set (or a combination of several discrete and continuous indices) if the number of possible outcomes is finite. The expressions throughout this article generalize replacing ∫ ξ d by ∑ ξ . Fixing the POVM and thinking of (7) as the distribution probability of the outcomes [as in Eq. (1)], one can use the tools introduced in Section (2.1.1). In particular, we can calculate the Fisher information and use the Cramér-Rao inequality to know if a given estimator is optimal. However, note that there is a dependence of the Fisher information on the POVM we choose. In order to have the lowest bound for the error, we maximize the Fisher information over all the possible measurements 23 is known as the quantum Fisher information, and through the Cramér-Rao inequality, Q 2 tells us the minimal possible error for the best measurement strategy for estimating a parameter appearing in the Hamiltonian of a quantum system. Equation (10) is a direct result from Eq. (5). Since Eq. (5) is valid for every POVM, it is valid for the one in which the maximum Fisher information is attained. The POVM that maximizes F Q is the one that should be used to get the smallest error in the parameter estimation 1 ; we call this POVM the optimal POVM. If the quantum state is pure there are analytical formulas for finding F Q . Otherwise no general formula is known and one must rely in numerical methods.

Bayesian Cramér-Rao inequality.
We consider the case where we have some partial knowledge of the parameter to be estimated. An example: We want to estimate the velocity of one particle in a dilute gas at temperature T . Without doing any measurement, we know that the particle velocity can be interpreted as a random variable that satisfies the Maxwell-Boltzmann distribution. Note that with the information we already have, we can estimate the velocity as the mean of the Maxwell- Here K b is the Boltzmann constant and m the mass of the particle. It is reasonable to design the experiment to measure velocities around v. The result of the measurement should improve the estimation, giving an error smaller than µ. Note that, in general, experiments designed to measure a parameter usually work for some expected range of its value, which implies assumptions were taken about the value of the parameter.
The Bayesian Cramér-Rao inequality can be used to decide what is the best estimator in the situation where we have partial knowledge of the parameter. We model the parameter as the random variable Θ, with outcomes θ, and probability distribution λ θ ( ). The outcomes of the experiment are modelled as the random variable Y, with outcomes y, and probability distribution p y ( ) θ | . The cost function we want to minimize is the mean square error It can be shown that the error is bound from below by the Cramér-Rao type inequality 13 , The first term of the sum is the expectation value of the Fisher information; the second term is the Fisher information of the probability distribution of the possible values of the parameter. The last term codifies what we already know about the parameter. As can be seen from the previous equation, the generalized Fisher information is larger than the Fisher information due to the knowledge we already have of the parameter. This has a simple interpretation: we can use λ θ ( ) to estimate the parameter, and measuring the system necessarily diminishes the error in the estimation of the parameter. The best strategy for measuring a parameter, with known information codified in a probability distribution, is given by the estimator, θˆ, that saturates inequality (12).
If we want to estimate the outcomes of a random variable, the problem of finding the best strategy is exactly the same as discussed in this section. In this case the experiment is repeated several times with the values of the parameter satisfying the probability distribution of the random variable. An example: measure the velocity of several particles in thermal equilibrium in a dilute gas. The velocities for the classical particles obey a Maxwell-Boltzmann distribution.
In this context, we found useful 24 , a review of bayesian inference in physics.
Bayesian quantum Crámer-Rao inequality. We want to estimate the parameter θ of a Hamiltonian in a quantum system where the known information about θ is codified in the probability distribution λ θ ( ). Given a POVM, ξ , a statistical model can be built through Eq. (7), and the problem is reduced to the classical one.
, that maximizes the generalized Fisher information (13), together with the appropriate estimator, saturates the Cramér-Rao type inequality   . For a combination with other weights, continuity and a recursive procedure imply convexity of F. From these two observations, we infer that the Van Trees information is also convex, as the integral of convex functions is also convex.
Since the maximum of the convex cost functions lies on the extremal points of all POVMs, we only need to search in this subset simplifying greatly the optimization task. A way to sample randomly such a set is presented in the following sub-section.
The algorithm. The outline of the algorithm is as follows: We produce a random POVM, and decompose it in extremals. We then evaluate the cost function using extremal POVMs and choose the one which yields the highest value. We repeat the procedure several times and keep the optimal POVM. We provide an implementation in an online repository 26 .
Random Sampling Method. To produce a random POVM, we use the purification algorithm 27 backwards, which transforms a general POVM into a usual projective measurement in an enlarged space.
The first step is to produce a random unitary matrix that acts on both the original space and an ancilla Hilbert space. The dimension of this ancilla space is the number of outputs of the initial POVM. Because the extremal POVMs have at most d 2 elements, where d is the Hilbert space dimension 14 , nothing is gained if the dimension of the ancilla space is larger than d 2 . For example, if we are trying to estimate a parameter of a qubit, we only need to use a dimension 2 Hilbert space and a dimension 4 ancilla space. If the Hilbert space is large or has an infinite dimension, we run the algorithm for some arbitrary dimension of the ancilla space and then increase its dimension until stable results are obtained.
The aforementioned unitary matrix is chosen with a measure invariant with respect to multiplication by unitary operators, i.e., with the Haar measure. The ensemble induced by the aforementioned measure is called the circular unitary ensemble (CUE) 28 . The easiest way to construct a representative member of the CUE is to construct a member of the Gaussian unitary ensemble (GUE) 28 , the ensemble of hermitian matrices invariant under unitary conjugation subject to the condition that the ensemble average of the trace of the square of the matrices is fixed. Generating a member H of such an ensemble is simple: we build a matrix A with Gaussian complex numbers, all with equal standard deviation and zero mean. Let † = + H A A, we can then calculate the matrix elements of † UHU , for any unitary U, and verify that the distribution of all matrix elements of the rotated matrix remains invariant. Thus, the eigenbasis of H has the Haar measure, and then the matrix U that diagonalizes H is a member of the CUE with the appropriate measure. The routine CUEMember calculates such an operator and can be found in 26 .
It is well known that we can interpret an arbitrary n output POVM as a projective measurement in a space built from the original one and a coupled n-dimensional ancilla space. If µ | 〉 µ= … { } n 1, , is an orthonormal base in the ancilla space, |Ψ〉 a state in the original Hilbert space, and Q { } m the operators characterizing the POVM, we can define a unitary operator U such that

⊗ | 〉〈 〉 = … is equivalent to the measurement defined by Q
{ } m , in the sense that the probabilities are the same, and the states, after discarding the ancilla space, are also identical, see 27 . Equation (17) can also be interpreted as a way to induce a POVM measurement in the original Hilbert space, starting from a unitary operator in an extended Hilbert space. In fact, if we replace |Ψ〉 by the basis state | 〉 j and premultiply by the bra 〈 |〈 | i m , we obtain Thus, from the random unitary U we can get all matrix elements of each of the Q m , according to Eq. (18). Since for all POVMs one can build a unitary transformation in the extended space such that Eq. (18) holds 27 , sampling all unitaries in the extended space guaranties sampling all POVMs with the corresponding number of outcomes. The routine POVM calculates the POVM in this way and can be found in 26 . The aforementioned method to sample POVMs will be called Random sampling method, or RSM for short. Naimark dilationAt this point, we would like to mention another POVM sampling method, inspired in Naimark's dilation theorem 29 . We start with a unitary matrix, acting on the original space tensored with an ancilla space  ancilla . The columns | 〉 v i of this matrix can be interpreted as an orthonormal basis on the extended space of dimension d′. Notice that d′ is a multiple of d, the dimension of the original Hilbert space. Let us define the d′ operators Conversion to a rank-1 POVM. To proceed further, we need a rank-1 POVM, so we must transform the aforementioned POVM accordingly. Recall that a rank-1 POVM is one whose elements are all rank-1 operators. However, typically, the Q { } m are not rank-1 operators.

For any given
that is not a rank-1 operator, we perform the spectral decomposition using a standard algorithm. Notice that the Q i are projectors, and 0 i λ > . We then replace Q with the l operators λQ i i , and the completeness relation remains valid since Notice that the number of elements at the end of the algorithm will be larger than the number of outputs of the initial POVM.
This Even though x a = is a solution to the problem, the usual numerical algorithms provide an extremal point. Notice that if an element x i of the solution is 0, it means that we do not include the operator in the POVM. The construction of the matrix A and vector b is done with the routine AConstruction. Let this extremal solution be x ext . The Mathematica LinearProgramming[c,A,b] implementation solves the following linear program: T where c is a constant vector. As we are interested only in finding a solution and not in optimizing a given vector, we give a random vector c to the usual Mathematica LinearProgramming routine. This vector c has random entries between 0 and 1 given by the uniform distribution with the Mathematica routine RandomReal. This process is done by our routine LinearProg.
To obtain the extremal POVM we start by defining x′ via ext with p a scalar. Requiring that x 0 ′ ≥ can be enforced letting i i i ext which in turn implies that p is a probability and that for some i, Indeed, Q ext is an extremal POVM 17 , and since one of the elements of x′ is null, ′ Q is a n 1 − output POVM (given that Q is an n output POVM) for which we can iterate the algorithm until a single output POVM is obtained. Notice that with this algorithm, all POVMs with a given number of outputs can in principle be sampled.
The routine BuildExtremal constructs an extremal POVM from the output of the Linear Program. The routine CalculateP calculates the p probability from Eqs. (21) and (22). Finally, the auxiliary solution ′ Q is built with the routine AuxiliarSol. Examples. In this section we apply the method described in the section Numerical algorithm to estimate the quantum Fisher information and the quantum Van Trees information. In the examples we observe a computational speedup when using the algorithm presented here compared with methods that randomly sample the whole POVMs space.
Qubit. We use the algorithm to calculate the quantum Fisher information for estimating a parameter encoded in a pure qubit state. We use known analytical results to benchmark our numerical method.
We consider a spin 1/2 particle in the state, where θ π ∈ [0, 2 ) is the phase between the two basis states and is a known parameter that characterizes the weight of each element of the superposition. The problem is the following: we want to find the best strategy to estimate the phase θ given that η is known 30 .
Given a set of states parametrized by ξ, i.e. ρ ξ ( ), we consider the following family of POVMs: Each of these POVMs has two elements that corresponds to the outcomes 1 and 2. In this subsection we shall consider the particular case Using Eq. (7) we find that the probability of measuring outcome 1 or 2 for the POVM ξ is given by  which is a function of the parameter we want to estimate, i.e. θ. Notice that there is a dependence on the initial state, via η, and on the POVM used, via ξ; we make these dependencies on the POVM explicit via a superscript. When the state is pure, the maximum quantum Fisher information can be analytically calculated 23 . In this example the quantum Fisher information is the maximum of F ( ) ( , ) θ ξ η with respect to ξ: In order to evaluate the performance of the proposed algorithm, we apply the RSM and RESM methods and compare their performance with the exact result (28). We define the errors where F RSM and F RESM are the Fisher information numerically calculated using the RSM and RESM respectively. In Fig. 1(b), we plot running time vs error for the two methods. It is clear from the plot that RESM is better, and the longer the program runs, the better the results using RESM compared with RSM. For this example, we obtain an error two orders of magnitude smaller running the program for the same length of time. Now we consider that we have some information about θ codified in the probability distribution θ p( ); limits to the error in the estimation are given by the Cramér-Rao type inequality (14). We assume that the angle θ has the uniform distribution p( ) 1/2 θ π = in Eq. (29). First we consider the maximization of the generalized Fisher information over the family of POVMs, ξ P( ), given by Eqs. (24) and (25) Because we are using a subset of all the POVMS ˜Z Z Q Q P ≤ , this approach allows us to get an analytic approximation for the quantum Van Trees information. For a uniform superposition (η π = /2), the Fisher information becomes independent of ξ; in fact ˜Z Z 1 Q Q P = = , see Eq. (27). This implies that any POVM from the family ξ P( ) maximizes the Fisher information. In general we obtain Q P so we can assert that if only POVMs of the family ξ P( ) are allowed, the best estimation is in the case where /2 η π = . Now we apply RESM to calculate Z Q and compare with Z ( ) Q P η , see Fig. 1(a). The maximum of Z Q is again obtained when /2 η π = . That means that the lowest error in the phase estimation is obtained when the weights of the superposition are the same. The figure suggests that Z Z Q Q P= . We observed that surprisingly, almost any extremal POVM is useful for finding the maximum Z Q . We require very few samplings (≈10) to observe good agreement between Eq. (30) and the numerical calculations. We also arrive at a reasonable estimate with 1 sample for most ηs.
Phase estimation. We calculate the Van Trees information for the phase estimation problem, one of the workhorses of quantum metrology. Since for this case no analytical solutions are known, this is an interesting testbed for our method.
Initial coherent state. We want to estimate the phase difference θ between two paths that light can follow, see 18 for a similar calculation. We probe the system with a coherent state, such that one path yields the state α | 〉 (with α a complex number) and the other where n is the number operator. Assume that we know, with some error, the size and the refractive index of the object that creates the phase difference. We can make an initial estimation of the phase difference between the two paths because it is proportional to the length travelled inside the object. We can model this situation assuming that θ is a random variable. As an example, we consider a Gaussian distribution centered at π, with standard deviation π/4 and trimmed at the edges (0 and 2π). Using RESM, we calculate the quantum Van Trees information for different values of α | |, see Fig. 2(a). The line is obtained using Eq. (24) with ( ) ( ) ( ) ρ ξ φ ξ φ ξ = | 〉〈 | as an ansatz. The figure suggests that the family of POVMS proposed is a good ansatz. When α | | decreases, the Van Trees information decreases, but when α | | = 0 still Z 0 Q > , since an estimation of the phase difference can be done with the information we already have about the parameter prior to any measurements. In order to do the calculation we approximate the coherent states with a truncated Hilbert space and limited the number of outcomes of the POVMs. We observe the results as a www.nature.com/scientificreports www.nature.com/scientificreports/ function of Hilbert space dimension and the number of outcomes of the POVMs. We find that a Hilbert space of dimension 7 and a POVM of 10 outcomes give us stable results. Initial displaced thermal state. As a final example, we consider estimating a parameter, chosen from a given distribution, encoded in a non-pure state. Again, there are no analytical expressions for the quantum Fisher information for this case. We calculate Z Q in order to bound the error in estimating the parameter. We build upon the last example, considering a thermal state displaced by the same operator that would give the state (31). Let be the state in which the parameter (θ) is encoded. For the Fig. 2(b), we used again a Gaussian distribution with mean π and standard deviation π/4. In Fig. 2(b) we show the numerical calculations of Z Q using the algorithm RESM. We compare it with the the ansatz composed of the two outcome POVM (24) with ρ ξ φ ξ φ ξ = | 〉〈 | ( ) ( ) ( ) , see Eq. (32). We see that the points calculated with the RESM algorithm which are a lower bound of Z Q beat the ansatz case for most points. We expect such behavior as the state in consideration (32) is a mixed state.

Discussion
The random extreme sampling method (RESM) can be used to find efficiently the maximum of a cost function over all possible quantum measurements. Particularly it is useful for finding limits in the precision of parameter estimation, through the cost function known as the quantum Fisher information, when the state to be measured is a mixed state. It can also be used to find the optimal measurement strategy by a given convex cost function by finding the POVM that maximizes it, at a considerable lower computational cost.  (29), resulting in Eq. (30), is shown as a blue line. We consider different families of states parametrized by η, see Eq. (23). The number of outcomes of the POVM is fixed to 4, as this is the maximum number of outcomes for an extremal POVM in a 2-dimensional Hilbert space 25 . Note that the RESM method gives much better results with one sample than the RSM method with 1000. (b) Error in the numerical estimation of F Q , see Eq. (28), sampling directly the whole space of POVMs (RMS) or only its extremal points (RESM) for η π = /2, with respect to the computational time invested. The slope for RSM case is m 0 63 RSM = − . and for RESM is m 1 84 RESM = − . . The error with the method proposed decreases much faster using RESM rather than RMS.

Methods
The code implementation can be obtained in the repository 26   (we varied r) with a random phase chosen from Gaussian distribution. It can be seen how the ansatz does not work. For the figure (a) we sampled 150 times, we considered a POVM with 10 outcomes and we truncated the photon space to the lowest seven states of the harmonic oscillator. For the figure (b) we sampled 200 times and the approximation of the dimension of the Hilbert space and the number of outcomes was the same, the dimensionless temperature is = − k T 10 .