Estimation of high-dimensional unitary transformations saturating the Quantum Cramér-Rao bound

We propose an estimation procedure for d -dimensional unitary transformations. For d > 2 , the unitary transformations close to the identity are estimated saturating the quantum Cramér-Rao bound. For d = 2 , the estimation of all unitary transformations is also optimal with some prior information. We show through numerical simulations that, even in the absence of prior information, two-dimensional unitary transformations can be estimated with greater precision than by means of standard quantum process tomography.


Introduction
The continuous advancement in the ability to control quantum systems and its application to the development of quantum technologies has driven the search for high-precision measurements and estimation methods.Quantum metrology aims to develop methods and tools that achieve the ultimate precision in parameter estimation.The enhancements provided by quantum metrology depend on the state of the probe, the quantum measurement, and the landscape of the parameters to be estimated.These are usually related through a multivariate variational problem that generally lacks analytical solutions.Despite this, quantum metrological improvements have already been demonstrated on various experimental platforms in both the single parameter and multiparameter cases [1,2,3,4,5,6,7,8,9,10,11].
Due to its intrinsic difficulty [12], the multiparameter case has remained less explored.In particular, the optimal measurements for different parameters are often incompatible [8] and the J. Escandón-Monardes: jescandon@udec.cloptimal probe states for different parameters can typically be different.Furthermore, in the multiparameter estimation, the quantum Cramér-Rao bound [13,14], which sets a fundamental limit for the covariance matrix, is generally not achievable even asymptotically [15,16].However, there are other Cramér-Rao type bounds, which are defined considering a weighted trace of the covariance matrix and are attainable in some scenarios [17,18,19,20,21].
An instance of multiparamenter estimation is the estimation of d-dimensional unitary transformations.Several methods to accomplish this task have been studied [22,23], particularly standard quantum process tomography (SQPT) [24], which has been successfully implemented for reconstructing quantum gates on ion traps [25], superconducting circuits [26], among many others [27,28].
Here, we propose a novel method for estimating d-dimensional unitary operations.Our approach requires a single target qudit, two control qudits, controlled gates, and Fourier transformations acting on the control qudits.The unknown unitary transformation acts on the target qudit.These resources allow mapping the unitary transformation to a state of both control qudits, which, after performing measurements, lead to an estimate of the coefficients that define the unitary transformation in the Weyl-Heisenberg basis.This is achieved without the need to measure the target qudit.Although our procedure does not provide a quantum metrological advantage, we show that it saturates the quantum Cramér-Rao bound for any initial target qudit state and simply uses measurements in the computational basis.This result holds for any finite dimension d > 2 and for unitary transformations that are close to the identity.In the case d = 2, all unitary transformations can be optimally estimated provided that the octant which the Bloch vector points to is known in advance.In this case, our estimation procedure agrees with previous results [29,30,31].In both cases the estimation scheme is independent of the parameters of the unitary being estimated.
We also estimate 2-dimensional unitary transformations without prior information, at the cost of losing estimation accuracy, and compare with SQPT.We simulate both procedures on Qiskit, IBM's software development platform for quantum processors [32], and study the average gate fidelity [33] as a function of the ensemble size, or number of shots, for a set of randomly generated unitary transformations.In the ideal case, that is, in the absence of error sources, our estimation procedure provides a better mean average gate fidelity than SQPT.Moreover, our estimation procedure shows a narrow standard deviation which means that all unitary transformations are estimated with similar average gate fidelity.In presence of noise affecting state generation, quantum gates, and measurements, our estimation procedure and SQPT exhibit similar mean gate infidelity, although the former with a larger standard deviation.Median gate fidelity of our procedure is larger than that of SQPT, which lays close to the inferior border of the interquartile range of our estimation procedure.Thereby, our estimation procedure provides a better estimation than SQPT in most cases.

Preliminary Material 2.1 Unitary transformations
An arbitrary d-dimensional unitary transformation U can be expanded in the Weyl-Heisenberg basis as where X and Z are the shift and phase operators, respectively.These act onto the canonical basis {|k⟩} with k = 0, . . ., d − 1 as X|k⟩ = |k ⊕ 1⟩ and Z|k⟩ = ω k |k⟩ with ω = exp(2iπ/d).The set {u m,n = r m,n e iϕm,n } of d 2 complex coefficients satisfies the unitarity constraint and characterizes U .With this, a general unitary can be written where we set ϕ 0,0 = 0 without loss of generality.

Quantum estimation theory
The classical Cramér-Rao bound states that the covariance matrix cov( t) of an unbiased estimator t of a parameter vector t is bounded below by the inverse of Fisher information matrix I(t), that is, for n repetitions of the experiment, which leads to limits for the accuracy of the estimate under various figures of merit (for recent reviews on the topic, see Refs.[12] and [34]).
The entries of the Fisher information matrix are defined as where P (y|t) is the probability of observing the value y in an experiment for a given parameter vector t.In the case of quantum mechanics the probability distribution P (y|t) depends on the measurement performed, which leads to different Fisher information matrices.Moreover, it is possible to derive (see Appendix H in Ref. [12]) a new bound known as the quantum Cramér-Rao bound given by where F is the quantum Fisher information matrix.Eq. ( 5) sets a bound to the achievable precision in the estimation of a set of parameters in the context of Quantum Mechanics.In the case of estimating a unitary transformation that acts onto a probe state |ϕ⟩, the quantum Fisher information matrix can be calculated as [12] F where we have

Estimation procedure
To estimate a d-dimensional unitary operation U , we propose a procedure which uses the quantum circuit shown in Fig. 1.This circuit is applied to the following initial quantum state of three qudits where the target qudit is in the arbitrary state |ψ⟩ 0 , and the qudits labeled 1 and 2 are the control states.The superscript i in |Φ i ⟩ 012 indicates the time step in the circuit.Each control qudit is subject to the action of a Fourier transform F |k⟩ = (1/ √ d) ω mk |m⟩, followed by the sequence of controlled shift and phase operators X (0) 02 and Z †(1) 01 defined by The action of the previous transformations leads to the probe state (9) The unitary transformation U to be estimated acts on the target qudit followed by the sequence Z (0) 01 and then X † (1) 02 and inverse Fourier transforms acting on each control qudit.Thereby, the initial state is transformed into the state and Z †(0) 02 disentangles the target qudit from the control qudits, followed by X 2 which correlates the indexes of the coefficients with the computational basis of both control qudits.The final state (see Appendix A for details) is where The coefficients u m,n entering in Eq. ( 1) are now in a one-to-one relation with the states in the canonical basis of the control qudits.Thus, the estimation of state φ 0 12 by any quantum tomographic scheme for pure states leads to the estimation of the unknown unitary transformation U .Moreover, simpler measurement schemes can be used, provided that a pure state of two qudits is defined by a set of 2d 2 − 2 independent real parameters while a unitary transformation acting on a single qudit is characterized only by d 2 − 1 real parameters.Some specific cases are studied below.

Estimation of 2-dimensional unitary transformations
Let us now consider the case of 2-dimensional unitary transformations.These can be written as where σ = (X, Y, Z) T is the Pauli vector, α ∈ [0, π/2], and n ∈ R 3 is a real unitary vector.After carrying out the exponentiation, we obtain the representation where we replaced Y = iXZ, and the coefficients are given by u 0,0 = cos(α), with θ ∈ [0, π] and ϕ ∈ [0, 2π[ being the spherical coordinates for n, and I the 2 × 2 identity matrix.Notice that u 0,0 is always non-negative, whereas the signs of u 1,0 , u 1,1 , u 0,1 depend on the octant in which the vector n points.Estimating an unknown two-dimensional unitary transformation U is thus equivalent to estimating the values of the angles (α, θ, ϕ).Projective measurements of control qudits lead to probabilities P 0,0 = cos 2 (α), ( 16) It follows that These relations allow for estimating the value of α, which is always in the interval [0, π/2].However, parameters θ and ϕ remain ambiguous, since u 0,1 , u 1,0 , u 1,1 are determined up to a sign.This ambiguity is removed when the octant pointed to by n is known beforehand, in which case our estimation procedure characterizes the unknown unitary transformation.Furthermore, it can be shown by direct algebra that our estimation procedure fulfills the equality I = F where and therefore our proposal saturates the quantum Cramér-Rao bound.Moreover, F is diagonal, hence our circuit is optimal in the sense that the three parameters defining U can be estimated simultaneously with the highest possible precision.The quantum Fisher information matrix in Eq. ( 21) was also obtained in other works [29,30,31].
The lack of a priori information does not prevent the use of our estimation procedure.As we show through numerical simulations in section 3.4, our procedure can be complemented with additional measurements and at the same time achieves better estimation accuracy than that obtained by standard process tomography.

Estimation of higher-dimensional unitary transformations
In order to estimate a higher-dimensional unitary transformation we need a clear dependence of U in terms of d 2 − 1 independent real parameters.For example, U can be parametrized in a similar fashion to Eq. (13) as U = exp i d 2 −1 j=1 λ j T j , where T j are the generalized Gell-Mann matrices and λ j are the independent real parameters.For d = 2 we expanded this expression in the Weyl-Heisenberg basis as in Eq. (14), where the complex coefficients u m,n explicitly depend on the corresponding parameters λ j (see Eq. ( 15)).However, it is not possible to obtain equivalent relations for d > 2 analytically, leading to difficulties in both the estimation of λ j and the calculation of the classical and quantum Fisher information matrices.Furthermore, finding a convenient parametrization that facilitates these calculations is an open problem.In this work, we consider the expansion of U from Eq. (2) where the elements in {r m,n , ϕ m,n } are not independent from each other.We restrict the analysis to unitary transformations that are close to the identity, since this approximation uncouples the parameters and allows us to obtain analytically both the classical and quantum Fisher information matrices.
To simplify the notation, we denote the coefficients u px,pz ≡ u p = r p e iϕp , where we have defined a single index p = (p x , p z ).These coefficients are constrained by the conditions and which enforce unitarity (see Appendix B).For unitary transformations close to the identity we have that r m /r 0 ≪ 1 for m ̸ = (0, 0).With this approximation, the only terms contributing in Eq. ( 23) are those where m = (0, 0) and m = ⊖p ≡ (d − p x , d − p z ).Thus, Eq. (23) becomes where we set ϕ 0,0 = 0 without loss of generality.
From the previous equation we obtain and (26) These conditions show that the amplitudes and phases of the coefficients u p are related in pairs.In the case p = ⊖p, Eq. (26) ties the phase to a discrete set, that is, Notice that the last restriction on the phases only occurs when d is even, and only for p = (d/2, 0), p = (0, d/2) and p = (d/2, d/2).Therefore, in the case d = 2, the three coefficients have a fixed phase up to a difference of π.
The constraints in Eqs.(25) and (26) allow us to recognize the d 2 − 1 parameters that characterizes U in the close-to-the-identity approximation.These are, for d = 2, three unpaired amplitudes r p .For d > 2 odd, all the coefficients are paired, implying that the relevant parameters are (d 2 − 1)/2 amplitudes r p and the same number of phases ϕ p .For d > 2 even, we have the three unpaired amplitudes r (d/2,0) , r (0,d/2) and r (d/2,d/2) , and (d 2 − 4)/2 other amplitudes and equal number of phases.
To handle these three cases at once we introduce the following partition of the set Z 2 d of indexes: where S 0 = {(0, 0)}, S u = {(d/2, 0), (0, d/2), (d/2, d/2)}, and S + and S − are any partition such that p ∈ S + if and only if ⊖p ∈ S − .In this way, we can easily identify the set of parameters defining U : Notice that every r a and ϕ a with a ∈ S + is respectively paired with r ⊖a and ϕ ⊖a , with ⊖a ∈ S − , via Eqs.(25) and (26).
Here and in what follows we use indexes f, g ∈ S u to label unpaired amplitudes and a, b ∈ S + to label paired amplitudes and phases.

|n⟩
, for n = 0 and n ∈ S u .
This operation can be understood as a set of Hadamard gates each acting in a subspace labeled with paired indexes, while acting as an identity on the other subspaces.Applying H on φ 0 12 we obtain the state Projective measurements on the computational basis for both qudits leads to the probabilities where ∆ a = ϕ a − ϕ ⊖a is given by the expression These probabilities lead to the estimates for the amplitudes and for the phases Thus, our proposal estimates the amplitudes and phases that characterize any d-dimensional closeto-the-identity unitary gate.In any case, the phases are estimated up to a set of four candidates, as implied by Eq. (36) and in agreement with the 2-dimensional case.The discrimination of the candidates requires prior information or additional experiments.
For dimension d even, the quantum Fisher information matrix characterizing our process is given (see Appendix C for details) by the block matrix where the ordering in the block matrix F even is given by ({r f }, {r a }, {ϕ a }), with f ∈ S u and a ∈ S + , and the explicit form of each block is being δ x,y the Kronecker delta.In particular, for d = 2, the unitary transformation is characterized by three unpaired amplitudes, i.e., S + and S − are empty.Hence, the quantum Fisher information matrix reduces to the upper left block as In the case of dimension d odd all amplitudes and phases are paired, hence S u is empty.Thus, the quantum Fisher information matrix is given by In this way, we have obtained the quantum Fisher information matrix F d for estimating close-to-the-identity unitary transformations in every dimension.Furthermore, we show in Appendix D that the classical Fisher information matrix I d is equal to F d .Thus, our estimation procedure saturates the quantum Cramér-Rao inequality for close-to-the-identity unitary transformations.Let us note that within the approximation r m /r 0 ≪ 1 the non-diagonal terms in F even are O((r m /r 0 ) 2 ), hence they can be neglected and consequently all the Fisher information matrices are nearly diagonal.In this way, in the case of odd dimension, the amplitudes can be estimated independently of each other and with equal precision.For even dimension, the amplitudes are also estimated independently; however, unpaired amplitudes are estimated with half the precision of paired amplitudes.Lastly, the precision in the estimation of the phases is severely restricted since it is proportional to the inverse of the square of the corresponding amplitude.
We must emphasize that our result is only valid within the close-to-the-identity approximation, where the parameters in the Weyl-Heisenberg expansion of d-dimensional unitary matrices, namely the phases {ϕ p } and the amplitudes {r p }, are approximately independent between each other.In the case far from the identity, where the parameters are strongly correlated, the above does not hold: even a slight variation in one of the amplitudes, which produces slight variations in the other ones because of the normalization, leads to violation of the unitarity conditions.Since the relation between the parameters for the general case is not clear, it is not possible to calculate, even numerically, the partial derivatives required to obtain the classical and quantum Fisher information matrices.In Appendix E we use a different parametrization to compare these matrices, showing that the distance between them increases as U is further from the identity.

Numerical simulations
In this section we study the performance of our estimation procedure for the case of qubit gates without prior information.This is achieved by measuring the quantum state in Eq. (12) in three different bases.The resulting statistics completely characterizes the unknown unitary trans-formation, as shown in Appendix F. We simulate our estimation procedure with Qiskit [32], IBM's software development platform for quantum processors, and compare its performance against the built-in function for SQPT.Since the output of SQPT is a quantum channel that is not necessarily unitary, contrary to our procedure which guarantees unitarity, we use as figure of merit the average gate fidelity, which compares a general quantum channel with a unitary quantum transformation [33].We generate a set of 200 single-qubit unitary matrices, which are randomly drawn from a uniform Haar distribution.Each unitary transformation is reconstructed 1000 times using our estimation procedure and additionally SQPT.In both strategies we calculate the average gate fidelity respect to the ideal unitary gate.This is repeated using increasing number of shots (or ensemble sizes) to simulate the measurement results.Finally, we also perform simulations considering various error sources (see Appendix G) and readout error mitigation.The code implementing this method is available in a Github repository [36].
Figure 2 shows the simulation results for our estimation procedure (green solid dots) and SQPT (red solid dots).Figures (a)  In the noiseless case, according to Figs.(a) and (c), both our estimation procedure and SQPT are characterized by almost indistinguishable mean and median average gate fidelity.In addition, in regime of a large number of shots, both estimation procedures exhibit extremely narrow standard deviation and interquartile range.In the small number of shots regime, our estimation procedure has a very rapidly narrowing standard deviation and interquartile range.Therefore, our estimation procedure and SQPT lead to an average gate fidelity that is independent of the unitary transformation.Figures (a) and (c) show that our estimation procedure achieves near-unit gate fidelity for ensemble sizes as small as 2 × 10 3 clearly outperforming SQPT.
Simulations show that the presence of noise affects the estimation of our procedure, decreasing the mean and median average gate fidelity with respect to their noiseless values, as exhibited in Figs.(b) and (d).The mean average gate fidelity of both estimation procedures becomes very similar, while the median gate fidelity of our estimation procedure is above that of SQPT, although standard deviation and interquartile range of our estimation procedure become wider.
In Fig. 3 we study the impact of different error sources on our estimation procedure and compare it to SQPT.Figures (a) and (b) show the mean and median average gate fidelity, respectively, for the noiseless case (solid green dots), noisy control-not gate (solid red pentagons), full noise with ideal control-not gate (solid blue diamonds) and full noise (solid pink squares) for our estimation procedure, and SQPT with full noise (black crosses).All simulations consider readout error mitigation as described in Qiskit documentation [32].The noise models and values used correspond to the ibmq_quito processor, and are provided in Appendix G. Insets depict the small ensemble regime.As expected, Fig. 3 shows that the full noise model leads to the the biggest decrease in the mean and median average gate fidelity.In the small ensemble regime, the estimation considering noisy control-not gates leads to a better mean and median gate fidelity than the estimation considering other error sources.However, as the ensemble size increases, the estimation procedure is clearly more affected by noisy control-not gates.Also, in this regime mean and median average gate fidelity become constant and the increase in the ensemble has no impact in the estimation accuracy.

Summary
In this work, we propose a procedure for estimating d-dimensional unitary gates.Our circuit transcribes the coefficients of the gate in the Weyl-Heisenberg basis into probability amplitudes of two control qudits.For qubit gates whose Bloch  vector points to a known octant, we show that our procedure saturates the quantum Cramér-Rao bound.In this case, the quantum Fisher information matrix is equivalent to the one derived in related works [29,30,31].We extended the analysis to higher dimensions and proved analytically that our procedure is optimal for unitary gates close to the identity in any finite dimension.In Ref. [29] it was shown that the quantum Cramér-Rao bound can be achieved for unitary transformations close to the identity, but no explicit protocol was proposed; our procedure accomplishes this task.Far from the identity our procedure still estimates the amplitudes r p of the coefficients, but the assessment of the precision of this estimation is an open problem.
In addition, we show that our estimation procedure is able to estimate any unitary transformation on qubit systems without requiring a priori information.Unitarity of the estimated operator is guaranteed by construction.Numerical simulations show that our estimation procedure outperforms SQPT in a noiseless scenario for every size of the ensemble and also in noisy scenarios with a small ensemble.Furthermore, considering a noisy scenario with ideal control-not gates, our procedure still outperforms SQPT for any ensemble size.berg uncertainty relations".Science Advances 7, eabd2986 (2021).

A Derivation of the output state
In this appendix we follow the evolution of the state through the circuit illustrated in Fig. 1.Before that, let us start with some preliminary definitions.
Let us consider a d-dimensional Hilbert space H and the d 2 -dimensional space L(H) of linear operators on H.We will denote the d th root of unity as ω ≡ exp(2πi/d), the addition modulo d as ⊕ and the sustraction modulo d as ⊖.Some important unitary operations in L(H) are the following: 1. Identity: 2. Shift operator : 3. Phase operator : 4. Fourier transform: It can be easily shown that the shift and phase operators satisfy the relation Also, we use a short notation for the Weyl-Heisenberg operators, defined as is an orthogonal basis for L(H); indeed, considering the Hilbert-Schmidt inner product for operators ⟨A, B⟩ = T r[A † B] and the relation where δ nm = δ nxmx δ nzmz .Consequently, any operator in L(H) can be written as a linear combination of Weyl-Heisenberg operators.In particular, an unknown unitary U can be expanded as

B Unitarity condition
Let us consider the unitary gate U , written in the Weyl-Heisenberg basis as and its adjoint where D n is defined in Eq. (A.6).We have but since U is unitary we also have U U † = I.Noticing that I = D 0 and considering Eq. (A.7), we have that unitarity of U implies: Let us calculate the left hand side of this expression using Eq.(B.3).We have: but in terms of the shift and phase operators we also have where we have used the commutation relation in Eq. (A.5).Now, summing over k we get In particular, for p = (0, 0) in Eq. (B.4) we obtain a normalization condition, which is equivalent to the total probability rule for the outcomes of our circuit: Besides, for p ̸ = (0, 0), we have

C Quantum Fisher Information Matrix
The entries F ab of the quantum Fisher information matrix are given by with and |Φ⟩ the probe state.In our circuit, the probe state is given by Eq. (A.20) as where |ψ⟩ 0 is the initial state of the target system.Let us consider again the expansion of the unitary gate where n = (n x , n z ) is an ordered pair and D n = D nx,nz = X nx Z nz .The sum can be split using partition (28) as where and we have used r a = r ⊖a in the last sum.Notice that We can now determine the derivatives of U and therefore the operators H a that appears in Eq. (C.1).For f ∈ S u we have: and therefore Similarly, for a ∈ S + we have: and In order to calculate the anti-commutators in Eq. (C.1), let us note that By imposing unitarity on U we have that U U † = I and Eq.(C.13) becomes We can now replace Eqs.(C.6), (C.9) and (C.10) into (C.14) and get: and In Thus, Note that Eq. (C.23) is different to zero only when n x = 0 and n z = 0, in which case Also, We can now find each ⟨Φ| H a |Φ⟩.By using Eq.(C.8), From Eqs.(C.24) and (C.26) we see that the first term in Eq.(C.27) is non-zero only for n = 0, whereas the second term is non-zero only for n = f .Therefore, Similarly, using Eq.(C.11) we have and using Eq.(C.12) we get (C.31) Given the partition (28), we have f, g ∈ S u .It follows that f ̸ = 0 and g ̸ = 0, hence the second and third term vanish.Thus, By replacing Eqs.(C.28) and (C.32) into (C.1),we obtain the entries of the first block of F: In the same way, from Eq. (C.16), we have e iϕa δ a,0 + e −iϕa δ a,0 + e iϕ ⊖a δ ⊖a,0 + e −iϕ ⊖a δ ⊖a,0 Notice that all the Kronecker deltas in this equation are null because their indexes belong to different sets of the partition (28).Thus, and we get the second block of F: From Eq. (C.17), we have Again, the Kronecker deltas are always zero because of the partition, hence and thus From Eq. (C.18) we have r 0 e iϕa δ a,0 + e −iϕa δ a,0 + e iϕ ⊖a δ ⊖a,0 + e −iϕ ⊖a δ ⊖a,0 Now, the only Kronecker deltas that do not vanish are δ a,b .Thus, but some Kronecker deltas are null because of the partition, while the terms in δ a,b cancel out.Thus, and Finally, from Eq. (C.20), we have Summarizing, the quantum Fisher information matrix for estimating a unitary gate close-to-theidentity in dimension even is with A, B, C, and D being the blocks defined in Eqs.(C.33), (C.36), (C.42) and (C.47), respectively.This matrix corresponds to the block matrix in Eq. (37).

D Classical Fisher Information Matrix
The classical Fisher information matrix is given by We can expand the sum using the partition (28): (D.2) Deriving the probablities in Eqs.(33) we get: Considering Eq. (34), Putting all these derivatives in Eq. (D.2), we can calculate I by blocks.We have: and Summarizing, we obtained the following classical Fisher information matrix for our procedure: where each block has the same definition of the corresponding block in the quantum Fisher information matrix in Eq. (C.48).

E Comparison of classical and quantum Fisher information matrices slightly away from the identity E.1 Alternative Parametrization of U : Gell-Mann matrices
In order to compare the classical and quantum Fisher information matrices, we firstly need to guarantee the unitarity of U while mantaining the independence of the parameters.A suitable way to do this is by writing an arbitrary unitary gate as the exponential of a Hermitian operator H: We can always write H as where T j are the generalized Gell-Mann matrices and λ j are d 2 − 1 real parameters (see definition in Ref. [37]).The set {T j } of the d 2 − 1 generalized Gell-Mann matrices plus the operator 2/d I is a set of d 2 linear operators satisfying the orthogonality relation T r[T † i T j ] = 2δ i,j , with δ i,j being the Kronecker delta.Hence, this set is an orthogonal basis for the space of d-dimensional linear operators.As a consequence, any operator U can be spanned as: where u (GM ) k are complex coefficients.In order to relate this expansion with the one in the main text, we will rewrite here Eq. (1) as: where u (W H) k are complex coefficients and D k are the Weyl-Heisenberg operators.Notice that every Gell-Mann matrix can also be spaned in the Weyl-Heisenberg basis.We can write where the j-th component t (k) j of the k-th Gell-Mann matrix can be calculated as Let us highlight that now we have two orthogonal bases in the space of d-dimensional operators.Moreover, we know from Eq. ( 22) that every unitary transformation is associated to a unitary vector when it is expanded in the Weyl-Heisenberg basis.This is not the case if we use the Gell-Mann basis, as can be easily seen by taking the identity, which in the Gell-Mann basis has a unique component different from zero with value d/2.We would like to normalize Gell-Mann matrices in such a way that any unitary operator is still associated to a unitary vector.To this end, we define the normalized Gell-Mann matrices as whose components in the Weyl-Heisenberg basis are These components define the vectors of an alternative orthonormal basis in a d 2 -dimensional vector space.Indeed, let us calculate the inner product of the k-th vector by the j-th vector: In the next subsection, we will use { t(k) } as the measurement basis for the control system in our circuit.

E.2 Estimating the parameteres of U (first order)
Let us consider an arbitrary unitary transformation U close to the identity.If we expand Eq. (E.1) to first order, we get U ≈ I + iH .(E.12) Using Eq. (E.2), Considering the normalized Gell-Mann matrices, we have and we can rewrite Eq. (E.3) as being the components of the unitary vector representing U in the normalized Gell-Mann basis.Thus, we have: Now let us consider again our circuit.By measuring the final state of the control system in the basis induced by the normalized Gell-Mann matrices, we get outcome probabilities from where we get an estimation of the parameters.Assuming λ j ≥ 0, we have Now we can use this estimator to assess the quality of the approximation when we go slightly away from the identity, and also compare numerically the classical and quantum Fisher information matrices.

E.3 Quality of the estimation
We assessed numerically the accuracy of the parameter estimation using Eq.(E.20).For each dimension d = 2, 3, 4 and 5 we created 10000 random unitary transformations close to the identity.Hamiltonian parameters were chosen randomly within a variable range up to [0, 0.1], in such a way that we could scan different distances to the identity.Although it is not properly a distance, we used 1 − r 0 as a measure of how close is U to the identity operator, with r 0 being the component of U in the Weyl-Heisenberg basis as in the main text (we also tried with the trace distance, obtaining equivalent results).For each unitary transformation, we estimated the parameters in the asymptotic limit using Eq.(E.20).The average relative error of the estimated parameters for each unitary is shown in Fig. 4. As expected, the accuracy of the estimation gets worse as U is further from the identity.For dimensions 3, 4 and 5 the approximation leads to an average relative error easily surpassing 1% when 1 − r 0 > 1 × 10 −3 , i.e. for r 0 < 0.999.For dimension d = 2, the average relative error remains below 1% in the range of parameters considered.

E.4 Classical and quantum Fisher Information comparison
For the same unitaries as the previous section, we calculated the quantum Fisher information matrix respect to the Hamiltonian parameters and the classical Fisher information matrix for two different measurement schemes: 1) when the control system of our circuit is measured in the computational basis (corresponding to the Weyl-Heisenberg basis in the space of linear operators) and 2) in the basis induced by the normalized Gell-Mann matrices.We used the trace distance between the quantum an classical Fisher information matrices as figure of merit.Fig 5 shows that for dimensions 3, 4 and 5 the classical Fisher information matrix is closer to the quantum Fisher information matrix when the control system is measured in the basis induced by Gell-Mann matrices.Moreover, this distance approaches to zero as U approaches the identity, while the distance between the classical Fisher information using the computational basis and quantum Fisher information matrices remains constant when U approaches to identity.In the case of dimension d = 2, where both bases are the same, the trace distance achieves lower values than in higher dimensions.

E.5 Far from the identity
The estimation of parameters performed above only considers a first order approximation of U around the identity.It would be interesting to go beyond this approximation to estimate parameters far from the identity.For example, we could approximate U to second order and calculate the Hamiltonian parameters from the outocome probabilities by using equations similar to Eq. (E.20).However, we could not find neither analytical nor numerical solutions for the nonlinear equations system not even in dimension d = 3.An alternative approach consists in having a good guess for the unitary gate (let us  call it Ũ ) and expand U to first order around Ũ .However, this approach is not straightforward, since it involves derivatives of operators which do not necessarily commute.Finally, it would be interesting to explore different parametrizations for unitary operators, what we let for future work.

F Qubit case without prior information
Any qubit unitary gate can be written as (13), which is expanded as where c I = cos(α) and c k = sin(α)n k for k ∈ {x, y, z}, with α ∈ [0, π/2] and n ∈ R 3 a unit vector.Notice that for the case of qubits the quantum state |Φ 9 ⟩ 012 in Eq. (A.37) becomes  For measuring the control qubits in the X basis we need to apply a Hadamard gate on each one (see Fig. 7).We get the new state In order to measure the qubits in the Y eigenbasis, we apply a π/2-phase gate, defined as S = |0⟩⟨0| + i|1⟩⟨1| , (F.5) 1. Get the magnitude of the coefficients r i from P Z jk .
2. If three of the r i are zero, then U ∈ {I, X, Y, Z}.
3. If r z = r x = 0, we only need to find s y using Eq.(F.11).In case it is exactly zero, it means that either r I or r y is close to zero; we choose s y randomly.
4. If r I = r y = 0, we get freedom in the global phase.Hence, we choose s x randomly and calculate s z from s xz and Eq.(F.9).
5. Otherwise, we calculate s xz and s x using Eqs.(F.9) and (F.11), respectively.Then we choose Eq. (F.14) or (F.16) according to which expression has the denominator with larger absolute value, in order to diminish the statistical variation on the estimation of s x .If both denominators are zero, we choose s x as the numerator of Eq. (F.16).
The codes for this algorithm are shown in github [36].

G Noise model
We use the NoiseModel class from Qiskit to incorporate noise in the simulation of our procedure.The simulation with full noise was accomplished by importing an approximate noise model of the ibmq_quito quantum processor of IBM.We consider only qubits 0, 1 and 2 from the five qubits of this processor, assigning qubit 1 for the target qubit in our protocol and qubits 0 and 2 as control register.
Table 1 shows the relevant error parameters for this setting, which are relaxation times T1 and T2, single-qubit gates error SX, readout or measurement errors and control-not gates error.A second noisy simulation was performed considering only control-not gate erros.To this end, we fixed single-qubit gates and readout errors to zero.To avoid relaxation, we kept relaxation times but  set the implementation time of each gate to zero.Finally, a full noise simulation but with ideal controlnot gate was performed.In this case, only control-not error was consider equal to zero.Table 2 shows the error parameters considered in these scenarios.These values are the maximum error parameters in ibmq_quito at the moment of setting the model and we consider the same values for each qubit. Error

Figure 1 :
Figure 1: Quantum circuit implementation of our estimation procedure.F are d-dimensional Fourier transforms acting on control qudits 1 and 2, X (i) tc and Z (i)tc are controlled gates defined in Eq. (8), and U is the unitary transformation to be estimated.State tomography of the control system leads to complete estimation of U .
Figure 2 shows the simulation results for our estimation procedure (green solid dots) and SQPT (red solid dots).Figures (a) and (c) show mean and median gate fidelity, respectively, as functions of the number of shots obtained in absence of error sources, that is, when the operations required by the estimation procedures are carried out perfectly.Figures (b) and (d) show results considering experiments with errors affecting single qubit gates, conditional gates, thermal relaxation, and measurements, an scenario that we call full noise model.This situation corresponds to a unitary gate embedded in a noisy device.Insets illustrate the behavior of estimation procedures in the small number of shots regime.Shaded areas show standard deviation in Figs.(a) and (b) and interquartile range in Figs.(c) and (d).

Figure 2 :
Figure 2: Plots (a) and (b) show the mean and plots (c) and (d) show the median gate fidelities as functions of the number of shots for both our estimation procedure (solid green dots) and SQPT (solid red dots).The left and right plots represent the noiseless and noisy cases, respectively.Shaded areas represent standard deviation or interquartile range.The insets illustrate the low number of shots regime.

Figure 3 :
Figure 3: Mean (a) and median (b) gate fidelities as functions of the number of shots for our estimation procedure with different noise models and SQPT with full noise.Shaded areas represent standard deviation or interquartile range.

Figure 4 :
Figure 4: Average relative error of estimated Hamiltonian parameters using first order approximation.Figures (a), (b), (c) and (d) show unitary transformations close to the identity in dimensions 2, 3, 4 and 5, respectively.

Figure 5 :
Figure 5: Trace distance between classical Fisher information (CFI) and quantum Fisher information (QFI) matrices.For dimension 3, 4 and 5, classical Fisher information is computed from probabilities obtained by measuring the control system in two different bases: the computational basis (WH, blue dots) and the basis induced by normalized Gell-Mann matrices (GM, orange dots).For dimension d = 2 both bases are exactly the same.

Figure 6 :
Figure 6: Measurement in Z eigenbasis.The circuit illustrates the qubit version of the circuit in Fig. 1, without the last X gate in the second control qubit.

Figure 7 :
Figure 7: Measurement in X eigenbasis.The circuit illustrates the qubit version of the circuit in Fig.1, without the last X gate in the second control qubit.Additional rotation is performed in each control qubit in order to measure them in the X eigenbasis.

Figure 8 :
Figure 8: Measurement in Y eigenbasis.The circuit illustrates the qubit version of the circuit in Fig. 1, without the last X gate in the second control qubit.Additional rotation is performed in each control qubit in order to measure them in the Y eigenbasis.
(nx,nz)|n x ⟩ 1 ⊗ |n z ⟩ 2 is the state of the control with exactly the same coefficients of U .For short, we write nz) X nx Z nz , (A.8) order to calculate the expectation values ⟨Φ| H a |Φ⟩ and ⟨Φ| {H a , H b } |Φ⟩ in Eq. (C.1), we are going to need ⟨Φ| D n |Φ⟩ and ⟨Φ| D † a D b |Φ⟩.We can simplify the notation for our probe state and rewrite Eq. (C.2) 1, without the last X gate in the second control qubit.

Table 2 :
For simulations with noiseless control-not gate we set its value to zero.For simulations considering only a noisy control-not gate we set to zero the other noise sources.