A variational toolbox for quantum multi-parameter estimation

With an ever-expanding ecosystem of noisy and intermediate-scale quantum devices, exploring their possible applications is a rapidly growing field of quantum information science. In this work, we demonstrate that variational quantum algorithms feasible on such devices address a challenge central to the field of quantum metrology: The identification of near-optimal probes and measurement operators for noisy multi-parameter estimation problems. We first introduce a general framework that allows for sequential updates of variational parameters to improve probe states and measurements and is widely applicable to both discrete and continuous-variable settings. We then demonstrate the practical functioning of the approach through numerical simulations, showcasing how tailored probes and measurements improve over standard methods in the noisy regime. Along the way, we prove the validity of a general parameter-shift rule for noisy evolutions, expected to be of general interest in variational quantum algorithms. In our approach, we advocate the mindset of quantum-aided design, exploiting quantum technology to learn close to optimal, experimentally feasible quantum metrology protocols.

Quantum metrology exploits non-classical effects to extend the sensitivity of sensing and parameter estimation methods beyond classical limits.The achievable precision in quantum metrology depends on the interplay of quantum correlations of the probe states, the unavoidable quantum noise present in the scheme, the geometry of the parameters to be estimated and the information that is gained by possibly intricate measurements [1,2].Carefully designing suitable quantum probes and measurements is therefore a frontier of quantum metrology [3] and can lead to significantly improved sensitivity as has been demonstrated in numerous works [4][5][6][7][8].Such methods have recently been generalized to a multi-parameter setting where a set of spatially distributed sensors is used to simultaneously estimate a number of distinct parameters or a function thereof [9][10][11][12][13], setting which is relevant for a broad range of applications ranging from nanoscale NMR [14] to networks of atomic clocks [15].
Designing close to optimal quantum protocols for a specific metrological task is, however, highly challenging.Identifying suitable probes and measurement schemes can be a classically intractable task even in the absence of errors.It requires optimizing over quantum states of high dimension, which soon becomes infeasible in practice due to the exponential growth of dimension with the system size.To add insult to injury, the solution of this task is even less obvious when realistic constraints such as experimental imperfections and decoherence are taken into account: Such decoherence processes, however, are at the heart of the matter when designing realistic quantum metrological protocols in the first place.This observation gives rise to the insight that in many relevant and meaningful scenarios, one cannot help but resort to quantum tools to find such protocols.
Fortunately, the increasing body of tools originating from the study of near-term noisy intermediate scale quantum (NISQ) computers [16] may come to the rescue here: This work provides a variational quantum algorithm for the optimization of quantum sensing protocols, applicable to a wide range of realistic metrology problems in the general multi- I f C(θ (t) , φ, µ (t) ) ρ(θ, φ) Figure 1.Illustration of the hybrid approach presented in this work to variationally optimize the probe state and measurement for a noisy unitary multi-parameter estimation problem.The various parameters are defined in the main text.
parameter setting.It can be implemented directly on controllable quantum devices, thereby circumventing the need to simulate large quantum systems classically.Variational quantum algorithms combine hybrid quantum-classical algorithms [17] with variational approaches that have been used in the context of quantum metrology to identify good probes in the single parameter setting [18] and for generating optimal spin squeezed states in optical tweezer arrays [19].We consider the generic task of estimating a multivariate function of the sensing parameters and use the multiparameter Cramér-Rao bound to quantify the quality of probe state and measurement as it provides for a saturable lower bound on the quality of any unbiased estimator.Importantly, we describe how this is extracted experimentally based on a novel insight regarding the parameter-shift rule [20] in a noisy setting, an idea that has arisen in the context of quantum op-timal control [21] and learning applications of quantum circuits [22] and which we develop further, expected to be interesting in its own right.Due to the variational nature of our quantum algorithm -involving a quantum circuit in its core instead of a possibly inefficient classical prescriptionwe can tailor estimation strategies to the specific experimental capabilities and dominating noise sources.Our techniques are very general and can be applied to both discrete variable architectures based on atomic systems [23,24] or superconducting qubits [25] and continuous variable architectures based on Gaussian light sources [26,27].As such, we add quantum metrology to the possible applications of NISQ devices while at the same time providing a solution to the challenging problem of optimizing quantum metrological prescriptions.
The structure of the algorithm.A high-level description of our algorithm is sketched in Fig. 1.In each step of the algorithm, a probe quantum state ρ 0 (θ) is generated by a quantum circuit parametrized by θ.It undergoes a noisy unitary transformation, consisting of a unitary evolution U (φ) that encodes the parameters, which are the arguments of f , the multi-variate function to be estimated, and an arbitrary, possibly non-local, quantum channel N accounting for the system noise.The state of the system before the measurement is The subsequent measurement is most generally described by a parametrized positive-operator valued measure (POVM) M = {Π l (µ)} resulting in measurement output probabilities with µ being the POVM parametrization.This concludes the part of the protocol that runs on the quantum device.This step is repeated a number of times to get accurate estimates of the outcome probabilities.These are then used to classically compute a cost function quantifying the estimation quality of the probe state.From this, both the state preparation circuit and the measurement procedure are updated to further increase the estimation quality based on gradient-descent techniques [28].The entire procedure is iterated until a minimum is reached, yielding a close to optimal sensing protocol within the variational manifold of probe state preparation and measurement procedure.
Cost function.The central challenge when constructing a variational quantum algorithm is to identify a cost function that captures the nature of the problem and that can be effectively evaluated and differentiated on actual quantum hardware.We want to quantify the performance of an estimator for a general multi-variate function f (φ) of the parameters φ.We therefore consider the classical Fisher information matrix (CFIM) [29] I f = J T I φ J where J is the Jacobian of f with entries J j,k = ∂f j /∂φ k and where we have used ∂ j as a shorthand notation for ∂/∂φ j .The CFIM gives a fundamental lower bound to the covariance matrix of any unbiased estimator f of f (φ) according to the Cramér-Rao bound where n is the number of samples.We note that the Cramér-Rao bound can, in principle, always be saturated in the limit n → ∞ by maximum-likelihood estimation (MLE) [30].
For single-parameter estimation, the CFIM reduces to a scalar quantity and the Cramér-Rao bound can directly be used as a cost function.For the multi-variate setting, we will follow the approach of Ref. [9] and apply a positive semidefinite weighting matrix W to both sides of the Cramér-Rao bound (4) and perform a trace to obtain the scalar inequality The right-hand side is the natural choice for the cost function While the outcome probabilities {p l } can be readily estimated through repeated measurements, the CFIM also contains their derivatives.These derivatives, and the derivatives of the cost function necessary for gradient-based optimization, are computed using the parameter-shift rule, which enables the calculation of analytic derivatives of expectation values on quantum hardware [20] for a range of fundamental quantum gates in discrete variable and Gaussian bosonic circuits.Each derivative is computed from the expectation values evaluated at two different points in parameter space.Recent work [31] generalized this approach and showed that analytic gradients for more intricate evolutions can be obtained via stochastic methods, even reaching beyond the noisy unitary model considered in this work.While initially derived for unitary quantum circuits, we prove in Appendix A that the parameter-shift rule and its stochastic generalization extend to noisy unitaries.As the most common elementary operations considered both in the discrete variable and the Gaussian continuous variable case admit a parameter-shift rule, we will in the following assume that the state preparation and the POVM in our scheme can be trained using the parameter-shift rule.
Realizing the optimization.The difficulty of running the variational algorithm naturally depends on the capabilities of the quantum device available.If we have the opportunity to perform the parametrized state preparation and measurement on the sensing platform, we can use it directly to run the algorithm.This is of course desirable since it circumvents the need to simulate the encoding and noise channel on another device.We can even do without a tomography of the noise channel as only the unitary components as such need to be known to good precision and have to be certified and benchmarked [32].The procedure is illustrated in Fig. 2(a).The output probabilities {p l } are directly accessible from the outputs of the experiment.To obtain estimates of the derivatives {∂ j p l }, it is, however, necessary to be able to re-run the same evolution with the parameter in question shifted by a fixed amount in order to use the parameter-shift rule.Whether this is possible depends on the exact nature of the unitary in the encoding evolution (see Appendix B for details).
If the necessary level of control is not achievable on the sensing platform directly, another quantum system with more control can be used to run the algorithm and optimize a sensing strategy.In particular, NISQ devices with high levels of control could be employed.We stress, however, that the inherent errors of the NISQ device need to be sufficiently low to provide for a faithful emulation of the sensing platform.The noise of the sensing platform can be recreated either through the use of ancillary qubits or repeated sampling of unitaries (see Appendix C).
If the noise process in question is hard to emulate on the NISQ device, or if the parameter-shift rule cannot be implemented for the encoding process, part of the algorithm can be moved from the quantum device to the classical hardware.This can alleviate the requirements for the NISQ device for a, in many realistic cases, modest overhead on the classical computation.The procedure in this situation is shown in 2(b).
We rewrite the output probabilities of the quantum circuit as (7) where N † is the adjoint channel of N and we defined the noisy POVM operators To compute the derivatives of the output probabilities, we use results on unitary estimation problems [33].The derivative of the probe after the noiseless part of the evolution with the Hermitian generators For a time evolution under commuting Hamiltonians H j , the generators are identical to the Hamiltonians H j = H j .
Using Eq. ( 8), the derivatives of the output probabilities become Defining the derivative operators allows us to recast the CFIM as a function of expectation values, similar in form to Eq. ( 3): where all expectation values are evaluated on the noise-free encoded state ρ(φ).
We note that the operators {Π l } and {∆ j l } have to be precomputed and stored classically in this hybrid algorithm.This requires an efficient classical representation of both the error channel and the POVM.This is feasible in practice, as the POVMs usually considered have unit rank.The rank of the noisy POVM operators is then maximally the Kraus rank of N .Locality of measurements and noise further simplifies the representation of the POVM operators.If necessary, they can also be stored in a format of tensor networks [7].
Because the CFIM needs to be evaluated on the noise-free encoded state, we need to parametrize the POVM {Π l (µ) } classically.If the POVM is parametrized in terms of quantum operations, we can again use the parameter-shift rule to compute the derivatives with respect to the parameters {µ k }.
Numerical experiments.To showcase possible applications of our variational approach, we numerically investigate two exemplary noisy estimation problems.The experiments have been implemented using the PennyLane library for quantum machine learning [34] and the QuEST mixed state simulator [35].
We first employ the setting of Ramsay spectroscopy, a widely used technique for quantum metrology with atoms and ions.The metrological parameters are phase shifts φ arising from the interaction of probe ions modeled as two-level systems with an external driving force.We follow Ref. [36] and model the noise in the parameter encoding as local dephasing with dephasing probability p.We consider a pure probe state and a projective measurement, where the computational basis is parametrized by local unitaries.We optimize a sensing strategy for a system of three ions and instead of estimating the phases independently, we analyze the task of estimating their average 1 N j φ j .We first reproduce known results to validate the performance of our approach.Ref. [9] has shown that generalized GHZ states are optimal when the encoding is noise-free.We have emulated the sensing task with φ j = π 6 , the sweet spot for GHZ sensing and performed co-optimization of state preparation and measurement.The parametrization was chosen so that it can produce any three qubit quantum state and any local measurement.Fig. 3 shows that we recover the optimality of the GHZ sensing procedure with a measurement in the Hadamard basis for p = 0.At .Semi-logarithmic plot of the Cramér-Rao bound for average phase sensing over the dephasing probability of (blue) the standard Ramsey probe and (orange) a GHZ probe.Each green dot marks the bound for a probe optimized for the particular dephasing probability.The phase shifts were fixed at the optimal spot for GHZ sensing at φj = π 6 .We recover the optimality of the GHZ probe in the noiseless case and find improved sensing protocols in the noisy case.increasing noise levels, the advantage of GHZ sensing disappears as expected [36] and we find sensing procedures that outperform both GHZ and standard Ramsay spectroscopy.
As a second task, we will apply our algorithm to the setting of spin imaging.This is at the heart of nanoscale NMR, which has a wide range of applications within chemistry and biological imaging.In particular, nanoscale sensors based on nitrogen-vacancy (NV) centers have shown great potential [14,[37][38][39][40] and we will consider the task of determining the position of a spin by triangulation with three NV center probes.In short, the dipole-dipole interaction between the NV centers and the spin shifts the energy levels of the NV centers resulting in a position dependent phase shift.From measurements of this phase, the position of the spin can be determined.We neglect the dipole-dipole interaction between the three NV centers assuming that they can be decoupled be appropriate pulse sequences [41].Furthermore, we perform a secular approximation, simplifying the interaction to a Pauli rotation about the symmetry axis of the NV centers [42].We model the encoding noise as local dephasing with a fixed probability of p e = 0.1.We refer to Appendix E for a more detailed derivation of our model.We use our algorithm to study the influence of gate noise, which we model as depolarizing noise with depolarization probability p g , on the quality of the sensing protocol.We compare a local probe to a probe that has two additional entangling operations while the measurement is always local.The exact structure of the sensing protocols is detailed in Appendix E. Fig. 4 shows that the shallow entangled probe yields an advantage in the absence of gate noise.But even at very small gate noise levels, the additional noise from the entangling operations will make the entangled probe inferior to the local probe.Although we only consider one specific ansatz for an entangled probe, this sheds light on the importance of gate noise in realistic probe preparations.
Summary and outlook.In this work, we have introduced a variational approach that can be applied to a wide variety of metrological problems to optimize multi-parameter estimation schemes both on sensing platforms and on quantum devices.These approaches can be tailored to specifically match the available hardware and the sensing problem in question.Instead of addressing the challenges of noisy quantum metrology by classically computing optimal probes and measurements, we have shown how to learn variational parameters of a quantum circuit for this purpose.The machinery used can also be used to improve upon an already good guess.While the approach outlined in the main text is widely applicable, covering arbitrary noise models and post-processing functions, it still does not cover all possible metrological tasks.In Appendix D, we therefore present extensions of our algorithm that take into account prior knowledge about the underlying parameters, handle mutual time-dependence of unitary evolution and noise, and outline how we can benchmark our results against the ultimately attainable precision given by the quantum Cramér-Rao bound.The results of Appendix A furthermore allow the extension of our method to the metrology of error channel parameters.We hope that the present work contributes to a shift in mindset -somewhat reminiscent of the idea of quantum computer-aided design of physical platforms [43] -so that not all sensing protocols have to be designed classically, but that one can make use of quantum technologies themselves to actually learn improved quantum protocols.
If the variational approach is directly applied to the sensing platform, the Bayesian nature of the underlying parameters is automatically taken into account.If the sensing is simulated on a quantum device, we can include the prior distribution by performing multiple runs of the experiment where φ is freshly pulled from P (φ) during each run.The CFIM can then be estimated from the results of the different runs.

Including time dependence
While the noisy unitary model is quite general, in many relevant cases both the encoded parameters and the noise depend on the exposure time t to the physical encoding process.Due to the interdependence of noise and unitary evolution, we can no longer use the parameter-shift rule.But we can still resort to finite differences to obtain the derivatives, both on the sensing platform and in the emulated context.In the simulation part, we have the additional possibility to compute the derivative directly.We will now assume that both the POVM (via the error channel) and the noise-free state depend on t.The time-derivative of the probabilities is then To calculate the channel-dependent part, we need the time-derivative of the error channel on the classical side of the optimization For the calculation of the state-dependent part, we can re-use the output probabilities we already calculated If the parameter-encoding is parameter-shift differentiable, this is already sufficient.If not, we also have to take care of the time-derivatives of the partial derivatives.The desired quantity is The first term encodes the error dependence, we analyze it using the definition of the derivative operators ( 11) where we have defined The second term encodes the dependence of the encoded state.We can analyze it by using Eq.(D4) and the differentiation rule Eq. ( 8): where we have defined With of the above definitions we arrive at We now have again reformulated all derivatives as expectation values at the cost of more complicated expressions.The derivatives can then be used in conjunction with Eq. ( 3) to obtain the CFIM and subsequently the cost function.

Evaluation of the quantum Cramér-Rao bound
Helstrom and Holevo [49,50] have shown that there exists a more fundamental bound than the Cramér-Rao bound that is independent of the chosen POVM, the Quantum Cramér-Rao bound (QCRB) where F φ is the quantum Fisher information matrix (QFIM) defined as where the symmetric logarithmic derivative (SLD) operators L j are implicitly defined by the equation The properties of the QFIM and its various applications in quantum metrology are discussed in detail in Ref. [51].
It would be satisfying to also be able to evaluate the fundamental precision limit for the given sensing problem that is given by the QCRB.Nonetheless, the calculation of the Quantum Cramér-Rao bound is much more intricate than the calculation of the classical Cramér-Rao bound.Formulas exist for unitary encoding processes, but calculating the QFIM for an arbitrary mixed probe state requires either full tomography or solving an infeasibly large semidefinite program [52].
For a pure probe state on the other hand, the calculation of the QFIM is straightforward, as it reduces to the fourfold covariance matrix of the generators [51] evaluated on the final state vector |ψ(θ, φ) .It can thus be estimated from measurements of the observables {H j } and their products and therefore also be differentiated via the parameter-shift rule.We can thus construct a similar cost function to Eq. ( 6), replacing the CFIM with the QFIM, and carry out an optimization.Current quantum devices usually do not reach the purity required to perform this optimization, so we have to resort to classical simulation instead.
Convexity of the quantum Fisher information ensures that the maximum Fisher information is attained on pure states.We can thus expect the pure picture to capture the optimal sensing performance when all noise is suppressed.This can serve as a benchmark to judge the performance of the probes obtained by the methods outlined above.Circuit.The following quantum circuit was used to perform the numerical simulations: The preparation and measurement circuits considered are detailed below.The encoding models the phase shift picked up during the interaction with additional dephasing noise D with dephasing probability p.We consider the following preparation circuit that parametrizes any pure three qubit state with the minimal number of parameters: Here, R P (θ) is a rotation generated by the Pauli word P .The measurement is parametrized by a local unitary transformation: Optimization.The encoded phases were fixed at φ j = π/6, the optimal parameters for a GHZ probe [36].The dephasing probability of the encoding was varied in the interval p ∈ [0, 0.5].The weights of the variational circuit was initialized at random, with each weight θ j , µ k randomly drawn from a uniform distribution over [0, 2π].The optimization was executed using gradient descent for 1000 steps with an initial step size of 0.05.A regularization term of 1 × 10 −10 × I was added to the computed CFIM before inversion during the calculation of the cost function to avoid inverting a singular matrix.Multiple runs of the above optimization were repeated and the best results selected.

NV trilateration
We consider a setup with three negatively charged NV centers (NV1, NV2, NV3) located in a diamond sample at positions r j = (x j , y j , z j ).We assume that these positions are known due to prior characterization of the system.The sensing task is to locate the position of a target spin J at an unknown position r = (x, y, z).As in Ref. [42], we are using the m s = 0 and m s = −1 ground states of the NV centers for the sensing task.
We consider the interactions between the NV centers and the target spin to be governed by the following dipole interaction Hamiltonian under the secular approximation [42] where S f n is the spin operator of the j'th sensing spin projected on its symmetry axis, which has unit vector n j .The spin to be sensed (target spin) has spin operator J b projected on the magnetic field axis (B with unit vector b).The distance between the j'th sensing spin and the target spin is δr j and has corresponding unit vector e j .The parameters ∆, γ e , and µ e are the zero field splitting of the NV ground states, the electron's gyromagnetic ratio, and the magnetic constant, respectively.Note that we have neglected interaction between the sensing spins, assuming that they can be decoupled by suitable pulse sequences [41].Performing the DEER sequence, the Hamiltonian simplifies to where we have neglected the local evolution of the NV centers ( j ∆(S n j ) 2 ) since this only adds a constant phase that can be mitigated.
From the Hamiltonian, we can identify the independent phases where and we have absorbed all constants into C.With this expression at hand, we use the autodifferentiation features of sympy [53] to calculate the Jacobian J j,k = ∂φ j /∂r k .The Jacobian necessary for our optimization can then be calculated by inverting J.
Circuit.The following quantum circuit is used to perform the numerical simulations: The preparation and measurement circuits considered are detailed below.The encoding models the evolution under the aforementioned Hamiltonian with additional dephasing noise D with dephasing probability p e .We seek to study the trade-off between a shallow entangled probe and a local probe at different gate noise levels.The gate noise will be modeled as depolarizing noise P with depolarization probability p g that is equal for single-and two-qubit gates.We consider two different state preparations which use the following primitives: First, an arbitrary single qubit operation Second, an entangling two-qubit operation The first preparation we consider is a simple local state preparation The comparison is made with a shallow entangling state preparation that extends the local state preparation by executing two entangling operations We note that this state preparation can create a GHZ state.For both preparations, we use a local unitary to parametrize the measurement Table I.Ground truth used for the NV trilateration experiment.The location is given in unit-free coordinates, the symmetry axis is given in spherical coordinates.Note that the symmetry axis in the case of the target spin is the precession axis identical to the magnetic field direction.
Optimization.The ground truth used for the numerical experiment is shown in Tab.I.The evolution time was fixed so that the total accumulated phase φ 1 + φ 2 + φ 3 = π/2.The dephasing probability of the encoding was fixed at p e = 0.1 while the gate depolarization probability was varied in the interval p g ∈ [0, 0.1].
The weights of the variational circuits have been initialized at random, with each weight θ j , µ k randomly drawn from a uniform distribution over [0, 2π].The optimization was executed using gradient descent for 1000 steps with an initial step size of 0.01.To avoid issues with overly large step sizes we used the following strategy: every time an increase of the cost function after a step was detected, the step size was multiplied with a factor of 0.7 and the step repeated.A regularization term of 1 × 10 −10 × I was added to the computed CFIM before inversion during the calculation of the cost function to avoid inverting a singular matrix.Multiple runs of the above optimization were repeated and the best results selected.

Figure 2 .
Figure 2. Realization of the approach proposed in Fig. 1. a) On a sensing platform or during its emulation, the output probabilities are directly computed from the POVM measurements.The derivatives are computed using the parameter-shift rule.b) When simulating on a quantum computer, the action of the noise channel is simulated by altering the POVM.The derivatives are computed by measuring the expectation values of a set of operators, used in training the network.

Figure 3
Figure3.Semi-logarithmic plot of the Cramér-Rao bound for average phase sensing over the dephasing probability of (blue) the standard Ramsey probe and (orange) a GHZ probe.Each green dot marks the bound for a probe optimized for the particular dephasing probability.The phase shifts were fixed at the optimal spot for GHZ sensing at φj = π 6 .We recover the optimality of the GHZ probe in the noiseless case and find improved sensing protocols in the noisy case.

Figure 4 .
Figure 4. Semi-logarithmic plot of the weighted Cramér-Rao bound for NV trilateration over the gate depolarization probability of (blue) local probes and (orange) shallow entangled probes.Each dot marks the bound for a probe optimized for the particular gate depolarization probability.The measurement is always local.The inset shows the positions of the NV centers and the target spin used in the experiment.The entangling operations only yield an advantage for vanishing gate noise.