Efficient computation of the Nagaoka--Hayashi bound for multi-parameter estimation with separable measurements

Finding the optimal attainable precisions in quantum multiparameter metrology is a non trivial problem. One approach to tackling this problem involves the computation of bounds which impose limits on how accurately we can estimate certain physical quantities. One such bound is the Holevo Cramer Rao bound on the trace of the mean squared error matrix. The Holevo bound is an asymptotically achievable bound when one allows for any measurement strategy, including collective measurements on many copies of the probe. In this work we introduce a tighter bound for estimating multiple parameters simultaneously when performing separable measurements on finite copies of the probe. This makes it more relevant in terms of experimental accessibility. We show that this bound can be efficiently computed by casting it as a semidefinite program. We illustrate our bound with several examples of collective measurements on finite copies of the probe. These results have implications for the necessary requirements to saturate the Holevo bound.

Except for special cases involving qubits [46] or estimating Gaussian amplitudes [47,48], in general the problem of finding the optimal measurement that minimises the sum of the mean squared error in multi-parameter estimation is a non-trivial problem. Instead, one resorts to finding bounds on these errors [49]. Some of these bounds are the bounds based on the symmetric logarithmic derivatives (SLD) [15,16] and the right logarithmic derivatives [50] as well as the Gill-Massar [51] bound. While these bounds are easy to compute, they are in general not tight. A tighter bound for the sum of the mean squared error which can be achieved in the asymptotic limit is given by the Holevo Cramér-Rao bound [52]. The computation of the Holevo bound was recently cast as a semidefinite program which has made it easy to compute. This was first performed for the Gaussian amplitude estimation problem [47] and was later generalised to an arbitrary model [53]. Furthermore analytic expressions which upper and lower bound the Holevo bound have recently been found [54]. In some special cases the measurement strategy required to reach the Holevo bound is known, for example with pure state probes [55] or for estimating a single parameter.
In general the Holevo bound is only asymptotically achievable [56][57][58], requiring a collective measurement over infinitely many copies of the probe state. A collective measurement here means that all copies of the probe state are measured simultaneously. In contrast a separable measurement restricts the probe states to be measured individually.
In practice collective measurements are extremely challenging to perform and are not accessible to most experimental teams. Thus it would be useful to have a tighter bound on the minimum achievable error when restricted to separable, single-copy measurements. One such bound for simultaneously estimating two-parameters was introduced by Nagaoka [59]. This bound is at least as tight as the Holevo bound and it can be saturated for probes in a twodimensional Hilbert space [60]. However, just like the Holevo bound, Nagaoka's bound is not an explicit bound-it requires a further non-trivial minimisation.
In this work we generalise the Nagaoka bound to estimating more than two parameters, and we call this generalised bound the Nagaoka-Hayashi bound. This bound applies to separable measurements on a finite number of copies of the probe state, unlike the Holevo bound which, as mentioned above, is only asymptotically attainable in general.
We further show that the minimisation required in the Nagaoka-Hayashi bound can computed using a semidefinite program. This makes its computation accessible. We illustrate our results with two examples which highlight some of the interesting features of finite copy metrology which are inaccessible with conventional techniques. In both of these examples we are able to find the positive operator valued measure (POVM) which saturates the bound, however whether this is always possible remains an open question.

II. RESULTS
Consider an n-parameter family of states {S θ |θ ∈ Θ ⊆ R n } in a finite d dimensional Hilbert space H q with θ = (θ 1 , . . . , θ n ) denoting the n independent true values that we wish to estimate. Let Π = (Π 1 , . . . , Π M ) be a column vector of M POVM elements, where (·) denotes partial transpose with respect to the classical subsystem. The quantum operators Π m are not transposed. This means Π m ≥ 0 and m Π m = 1. Each outcome m assigns an estimated value for θ j through the classical estimator functionθ jm . The standard measure of estimation error when restricted to separable measurements is through the n-by-n mean squared error (MSE) matrix V θ (Π,θ) with entries for j, k = 1, . . . , n .
The notation Tr [·] in serif font is used to represent the trace of an operator in H q , the Hilbert space of the quantum system. For brevity of notation, hereafter we drop the argument and write the MSE matrix as V θ . We aim to minimise the trace of the MSE matrix under the condition that our estimates are locally unbiased m Tr[S θ Π m ]θ jm = θ j and The Nagaoka bound for two-parameter estimation gives a lower bound on the trace of the MSE matrix as [59] Tr where the sans-serif font Tr[·] denotes the trace of a classical matrix in H c , TrAbs A is the sum of the absolute values of the eigenvalues of the operator A, and X = (X 1 , X 2 , ..., X n ) is a vector of Hermitian estimator observables X j that satisfy the locally unbiased condition at θ The Nagaoka bound was conjectured to be a tight bound for Tr[V θ ] [60].

A. Computable multi-parameter bound
As we shall shortly prove, the Nagaoka bound can be generalised to more than two parameters. This result is stated as the following theorem.
Theorem 1 (Nagaoka-Hayashi bound). Let V θ be the MSE matrix of an unbiased estimate of θ for a separable measurement on a model S θ . Then the trace of V θ is bounded by where S θ = 1 n ⊗ S θ and L is an n-by-n matrix of Hermitian operators L jk .
We use the symbol Tr[·] to denote trace over both classical and quantum systems, i.e. over both H q and H c . We call this bound the Nagaoka-Hayashi bound. However the Nagaoka-Hayashi bound is not an explicit bound. Our second main result is that this bound, c NH can be computed as a semidefinite program.

B. Related bounds
Before proceeding on the proof and computation of the Nagaoka-Hayashi bound, we digress briefly to mention two related bounds. The first is the Holevo bound which can be written as [52] Tr where we require the classical estimatorx and POVM Π to satisfy mx jm Π m = X j .
The derivatives of the state S with respect to θ do not play any role here. A bound on the trace ofŨ θ is given by Hayashi's bound [61] Tr As Hayashi's work is only available in Japanese, we summarise its main results in appendix A. If the given matrices X happen to satisfy the locally unbiasedness condition (4) for θ, thenŨ θ also forms a valid parameter-MSE matrix for those θ. In this case, because of the additional restriction (8), it is clear that c NH−U ≥ c NH . Also in this setting, Watanabe et al. [62] derived bounds for estimating two observables when restricted to certain classes of random and noisy measurements. In the case when both the observables and state S are two-dimensional, these bounds are achievable. In fact, when the number of observables n = 2, the minimisation over L can be performed analytically and c NH−U takes the explicit form

C. Proof of main results
In this section, we shall prove Theorem 1. To that end, we need to introduce some definitions. We rewrite the elements of the MSE matrix as where the MSE-matrix operator L θ (Π,θ) is an n-by-n matrix with operator elements. We introduce a classical matrix ξ with elements ξ jm :=θ jm − θ j so that where we have set n = 3 to simplify the presentation. The generalisation to arbitrary n is straight-forward. With this notation, it is clear that L θ is an operator on the extended Hilbert space H c ⊗ H q . To anticipate the proof, it is useful to write L θ in the following form where M is the number of POVM outcomes and Ξ ij = ξ ij 1. We can also introduce the following extension to S θ , S θ = 1 ⊗ S θ so that the expression for the MSE matrix can be written as We are now ready to prove Theorem 1.
Proof. Suppose the optimal POVM and unbiased estimator have been found and are given by Π andθ which leads We use asterisk to denote the optimal values and optimal operators. From Π andθ, we can construct the estimator matrices so that Comparing the above with (14) and using the result which holds because Π j are positive operators that sums up to 1 (see Proposition II.9.1 of Holevo [52]), we arrive at ≥ min In the two parameter case, we show in appendix B that c NH reduces to the original Nagaoka bound c N in (3). More where L jk = L kj Hermitian and X j Hermitian satisfying the conditions (4) for local unbiasedness. The conversion to a standard semidefinite program is performed in appendix D. We also show in the same appendix that the worst case computational complexity for solving the SDP to an accuracy is O (nd) 3/2 log(1/ ) .
The computation of the Holevo bound c H was shown to be a semidefinite program by Albarelli et al. [53]. The difference between the Holevo bound and the Nagaoka-Hayashi bound is that in the former, the optimisation is performed directly on the covariance matrix V = Tr[S θ L] while in the latter the optimisation is performed on the operators L.
We note that both programs can also be applied to compute the bound on the operator-MSE c NH−U (9) with little modification-the only changes needed are to replace the minimisation variables X with the given observables and ignore the conditions (4).

III. EXAMPLES
In the following, we demonstrate our results by computing the Holevo and Nagaoka-Hayashi bounds for two illustrative examples-the estimation of orthogonal qubit rotations on the Bloch sphere in a phase damping channel and the simultaneous estimation of phase and loss in an interferometer. In the former we find that the Holevo bound is always smaller than the Nagaoka-Hayashi bound, and in the latter we find that the two bounds are always equal.
The minimisation problem was solved with the Yalmip toolbox [63] for Matlab using the Mosek solver [64].
Even though the semidefinite program only returns numerical values for X and L, in some of these examples, the analytical forms for them can be inferred from the numerical solutions. Furthermore, every semidefinite program (24) has a dual program that involves performing a maximisation over the Lagrange multipliers associated with the primal program [65]. That the inferred solutions are indeed optimal can then be verified by checking that the values for the primal and dual programs coincide. For both of the examples considered we present the dual solutions in appendix G.
A. Example 1: Estimation of qubit rotations with a two-qubit probe Our first example concerns estimating the rotation experienced by qubit probes subject to the phase damping channel. This channel has particular relevance for modelling decoherence in trapped ions [66][67][68]. We consider the maximally entangled two-qubit state (|01 + |10 ) / √ 2 as a probe. The first qubit acts as a signal-probe which passes through a channel imparting three small rotations: θ x , θ y and θ z about the x, y and z axis of the Bloch sphere. The rotated probe is then subject to the phase damping channel E with a known damping strength The second idler-qubit is stored in a perfect quantum memory and remains unaffected by the rotation or phase damping. The resulting two-qubit state then has an approximate matrix representation in the computational basis as which is valid to the first order in θ. The partial derivatives of S θ with respect to θ evaluated at θ = 0 are

Single parameter estimation
Let's start with the simple case when θ y = θ z = 0 and we are only estimating the single parameter θ x . In a single parameter estimation problem, the Holevo bound coincides with the Nagaoka-Hayashi bound and can always be saturated by a separable measurement. In this case, the two bounds can be achieved by the estimator operator which gives c H,1 = c NH,1 = 1, independent of . The optimal measurement that saturates this bound is a projective measurement on the four orthogonal eigenvectors of X x This together with the estimation coefficients ξ = (1, −1, −1, 1) gives an estimation variance of v x = 1. The phase damping channel has no effect on the estimation precision.

Two parameter estimation
Next, for estimating the two parameters θ x and θ y when θ z = 0, the Holevo and Nagaoka bounds no longer coincide.
The optimal matrices that achieve the minimum in the Holevo bound are found to be which gives c H,2 = 2. This means that there exists a sequence of collective measurements that can saturate a variance of v x = v y = 1 in the asymptotic limit.
Unlike the single parameter case, the optimal X x and X y operators for the Nagaoka bound are different from those which optimise the Holevo bound. For the Nagaoka bound the optimal matrices are which gives c NH,2 = 4/(2 − ). Since there is a gap between the Holevo and Nagaoka bounds, a separable measurement cannot saturate the Holevo bound-a collective measurement is required. We show in appendix E that the Nagaoka bound is saturated by a family of five-outcome POVMs which gives v x = v y = 2/(2 − ). This means that when restricted to separable measurements, this is the smallest pair of variances possible.
FIG. 1. Holevo bounds (solid lines) and Nagaoka-Hayashi bounds (dashed lines) in terms of average preciseness for estimating two (blue) and three (red) orthogonal rotation parameters simultaneously using a maximally entangled two-qubit probe under the action of the phase damping channel. The Nagaoka-Hayashi bounds can be achieved by a separable measurement on a single probe, while the Holevo bound require a collective measurement on possibly infinite number of copies. The shaded area shows the gap between the two bounds. For estimating a single parameter, the Holevo and Nagaoka-Hayashi bounds coincide and are equal to the two parameter Holevo bound.

Three parameter estimation
Finally for estimating all three angles θ x , θ y and θ z simultaneously we find the Holevo and Nagaoka-Hayashi bounds are c H,3 = 2 + 1 (1 − ) 2 and Just like the two parameter case, the gap between the two bounds implies that a collective measurement is required to saturate the Holevo bound. These bounds are achieved by the same estimator operators (30) for the Holevo bound and (31) for the Nagaoka-Hayashi bound with the additional We write down an explicit POVM that can approach c NH, 3 with showing that this bound is tight.
In order to quantify the estimation accuracy, we define the average preciseness for simultaneous estimation of n parameters with n/(v 1 + · · · + v n ) as a figure of merit on how good the estimators perform. By construction, a large average preciseness implies that all n parameters can be determined accurately. We plot this quantity in Fig. 1 for all three estimation cases. We also note that in the two and three parameter examples, it is easy to check that the SLD Fisher information matrix is diagonal. Furthermore the model is asymptotically classical and the Holevo bound coincides with the SLD bound [21,69].

The Nagaoka-Hayashi bound for multiple copies of the probe state
We now demonstrate the usefulness of the Nagaoka-Hayashi bound and the associated SDP by examining the precision limits when we perform collective measurements on finite copies of the probe state. We denote the Nagaoka-  Fig. 2 shows how the gap between the two bounds shrinks for an increasing number of copies of the probe state. We consider up to three copies of the probe state. Without the Nagaoka-Hayashi bound a brute force search for the optimal measurement strategy for three copies would require optimising an M outcome POVM, where each outcome is a 64-by-64 matrix. Thus the Nagaoka-Hayashi bound and the associated SDP offer an efficient way to investigate the asymptotic attainability of the Holevo bound. It provides a tool to address how fast optimal estimators on finite copies converge to the asymptotic bound.

Discussion of qubit rotation example
This example demonstrates several interesting features of finite copy metrology. First we are able to definitively show that there exists a gap between the attainable precision with collective and separable measurements. Without a separable measurement bound such a claim is not possible as any gap between a numerically optimal POVM and the Holevo bound may be a result of a deficiency in the numerical search as opposed to a physically meaningful gap.
Secondly as we are able to find a POVM which coincides with the Nagaoka-Hayashi bound, we are able to say with certainty that this POVM is optimal. Finally we are able to investigate the attainability of the Holevo bound. While it is known that the Holevo bound is asymptotically attainable, it is not known how many copies of the probe state are required to get close to the Holevo bound. As mentioned above to investigate this numerically with a POVM search is computationally very expensive. The SDP presented circumvents this and allows us to investigate the attainability of the Holevo bound in a numerically efficient manner.

B. Example 2: Phase and transmissivity estimation in interferometry
In our next example, we consider the problem of estimation of phase change φ and transmissivity η in one arm of an interferometer as shown in Fig. 3. Following Crowley et al. [33], we consider initial pure states with a definite photon number N across the two modes |ψ in = N k=0 |k, N − k a k , where |N, M represents a state with N photons in the first mode and M photons in the second mode. One family of states with a fixed photon number is the Holland-Burnett states which are obtained by interfering two Fock states with an equal number of photons on a balanced beam splitter. These states lead to a phase estimation precision better than an interferometer driven by a coherent light source with the same number of photons [70]. The Holevo bound for the Holland-Burnett state was computed by Albarelli et al. [53] for up to N = 14. In general, the Holevo bound requires a collective measurement on several probes to be saturated. But for some values of N and η, the Holevo bound can be saturated by a separable measurement, Π (φ) that optimally measures the phase [53].
We compute the Nagaoka bound for these states for different values of η with φ = 0 for N up to 14 using our SDP.
We find that the Nagaoka and Holevo bounds always coincide (up to numerical noise). This is to be expected when Π (φ) saturates the Holevo bound, but is not so obvious when it does not. The fact that there is no gap between the Holevo and Nagaoka bound implies one of two possibilities: either (i) the Nagaoka bound is not tight or (ii) separable measurements are always optimal for simultaneous estimation of φ and η, in other words, collective measurements cannot do better. In the following, we show that the second statement is true.

Measurement saturating the Nagaoka bound
The initial pure state |ψ in = N k=0 |k, N − k a k transforms in the lossy interferometer channel to the following state where each term in the direct sum represents a state with l lost photons. The state S φ,η is a mixed state with rank N + 1. Here b kl = k l η k−l (1 − η) l are the beam-splitter coefficients and p l represents the probability of losing l photons. The partial derivatives of S φ,η share the same direct sum structure with each block having at most rank 2. Thus what we have is a direct sum of pure state models, and for such a model, we have a separable measurement with a direct sum structure that can achieve the Holevo bound [55]. Each block can be measured separately but we cannot minimise v η + v φ separately in each block. This is because how much weight we attach to η or φ in one block will depend on how much information about them that we can get from the other blocks. But regardless of the weights, each l = N block requires at most a 3 outcome POVM to saturate the Holevo bound, so the total number of POVM outcomes needed is at most 3N + 1. The extra 1 comes from the l = N block where all photons are lost. An analytic POVM that saturates the Holevo bound for the N = 1 case is given in appendix F. The dual solution to the Nagaoka-Hayashi bound is presented in appendix G 2.

Discussion of optical interferometry example
This problem demonstrates a very different but equally insightful feature of finite copy metrology compared to the qubit rotation problem. The simultaneous estimation of phase and loss has been very well studied in the literature [6,8,33,53], however until now the fact that separable measurements are sufficient to reach the ultimate attainable precision had remained unknown. This insight was only possible with our SDP, which allowed the Nagaoka-Hayashi and Holevo bounds to be compared for large N . We plot the numerically calculated Nagaoka-Hayashi and Holevo bounds for different N and η in Fig. 4 is essential and when it is not. Our results can be applied to many other problems in multi-parameter quantum metrology and will help quantify the maximal advantage collective measurements have to offer. In some cases, a separable measurement is already optimal, simplifying any experimental realisation.
In the first example, we have assumed that the damping strength is known. However in a practical setting, it would be more realistic to consider as a nuisance parameter, an unknown parameter that we are not interested in which nevertheless may hinder our measurement precision [25,42,71]. The quantum Cramér-Rao bound in the presence of nuisance parameters can be computed utilising a low-rank weight matrix [25,71]. As we show in appendix C, our SDP formalism can be immediately applied to such cases. An interesting extension to this work would be to investigate examples which incorporate nuisance parameters.

DATA AVAILABILITY
The data that supports the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code that supports the findings of this study are available from the corresponding author upon reasonable request.

Computation and Communication Technology (Grant No. CE170100012). JS is supported by the UEC Research
Support Program, the University of Electro-Communications. We are grateful to Professor Nagaoka and Professor Hayashi for helpful discussions.

COMPETING INTERESTS
The authors declare that there are no competing interests. We summarise Hayashi's result [61] which was published in the proceedings of a domestic workshop in the Research Institute for Mathematical Sciences (RIMS) at Kyoto University in Japanese for the reader's convenience. Let H q be a finite d-dimensional Hilbert space and consider a set of observables (Hermitian matrices) X = (X 1 , X 2 , ..., X n ) on it. We say that a POVM Π = {Π m } is a simultaneous measurement of the given observables X, if holds for all j. In general, a projection measurement does not exist unless the X j commute with each other, but a POVM Π exists. Given a state S on H q , we define the expectation value of X j by We define the covariance matrix by We are interested in minimizing the sum of the diagonal elements ofŨ(Π,x). As the second term is constant this is equivalent to minimising Tr[U]. Indeed the second term can be ignored for all practical purposes. We define the precision limit as Note here that C depends on the given state S and the set of observables X. Hayashi derived the following two bounds for C.
Theorem A.1 (Hayashi). The following are lower bounds for C and further that C ≥ C 1 ≥ C 2 holds.
where S = 1 ⊗ S and U are complex matrices on the extended Hilbert space H c ⊗ H q .
Hayashi's first bound C 1 is considered as the generalisation of the Nagaoka bound for simultaneous measurement of non-commuting observables [59]. Using the linear programming approach, Hayashi further derived the following alternative forms for C 1 and C 2 where Sym ± (A) = 1 2 (A ± A ) is the symmetrised (anti-symmetrized) matrix of A on H c ⊗ H q with respect to the classical index.
Finding the fundamental limit C is still an open problem. For two observables, Nagaoka conjectured that the bound C 1 is tight [60]. In other words, C = C 1 .

Appendix B: Nagaoka bound for two parameter estimation
The Nagaoka bound for the two parameter estimation case is [59] c N = min with X j Hermitian satisfying (4) in the main text. In this appendix we show that in the two-parameter case, the Nagaoka-Hayashi bound, (5) in the main text, coincides with the original Nagaoka bound. When n = 2, the Nagaoka-

Hayashi bound is
with L jk Hermitian and X j Hermitian satisfying (4) in the main text. We can write the condition in (B2) as Recognising that [X 1 , X 2 ]/2 is an antihermitian matrix which we label as iH, we can rewrite the condition as where L denotes the matrix on the left hand side of (B3). In order for this matrix to be positive we require [72] for any unitarily invariant norm. This inequality can be saturated by the choice where |H| = √ H 2 . The following corollary ensures (B4) is satisfied.
Corollary B.1. Let A be any matrix. Then the matrix The matrix L can be chosen in this way by optimising over the matrix L so that   L 11 L 12 We let this norm be TrAbs, which is equal to the trace for the left hand side of this equation and so the condition Rearranging and including S θ we arrive at We now present an alternative method of proving the equivalence between the two bounds for the benefit of the interested reader. We define a complex weight matrix on the extended Hilbert space, whose real part is the identity for simplicity.
From positivity of W, the parameter w takes values in [−1, 1]. By partial tracing over the parameter space after multiplying √ W from both sides, Eq. B3 implies Since w is arbitrary real in [−1, 1], we obtain We then use the following well-known Lemma (see for example Lemma 6.6.1 of Holevo [52]).
Lemma B.2. For a given Hermitian matrix Z, suppose Y obeys inequalities Y ≥ ±Z. Then, the minimum of the trace of Y is given by This gives the bound As before the left hand side of Eq. B14 can be chosen so that the two sides are equal. Thus, min Tr[S θ L] under the constraints is given by Note that this method works for more general weight matrices.

Appendix C: Generalisation to arbitrary weight matrix
We present a generalisation of our main results to an arbitrary weight matrix W ≥ 0. In the case where the weight matrix W > 0 is full rank, it can be set to the identity after a suitable reparametrisation for the model (see for example, Sec. V of Fujiwara and Nagaoka [73]). Since we are only interested in local bound, this reparametrisation does not matter. Specifically, we can reparametrise the model as ϕ j = k H jk θ k where H = √ W is a real and regular matrix. Estimating the new parameters ϕ is equivalent to estimating the original parameters θ with a weight matrix W.
When W is not full rank, a bit more care is required in reparametrising the model because it might be possible that some of the new parameters ϕ j are exactly zero or that two of the ϕ j 's might be identical. This situation is common when studying parameter estimation in the presence of nuisance parameters [25,42,71]. Nonetheless, it is still easy to incorporate the weight matrix W into our original framework. We now wish to minimise Recalling that the MSE matrix can be written as V θ = Tr[S θ L θ ], this is handled by noting the following where S θ = (W ⊗ 1)S θ = W ⊗ S θ is a positive semidefinite matrix. Thus, by changing from S θ to S θ , nothing about the problem changes and it can be solved using the same SDP as in the main text.

Appendix D: Conversion to standard SDP and complexity discussions
Here we show that the program with L jk = L kj Hermitian and X j Hermitian satisfying (4) in the main text can be converted to the standard SDP program where Y is a positive-semidefinite Hermitian matrix of size nd + d having the form Y = of H q and m is the total number of constraints on Y . The objective function to be minimised is handled with There are five groups of constraints on Y that have to be implemented through F k and c k . Denoting S j = ∂S θ ∂θj , the constraints are: 3. X j Hermitian. 4. L jk = L kj Hermitian.
5. The lower n-by-n block of Y equals the identity operator.
In the following, we set n = 3 to simplify the notations. The group 1 constraints are achieved with the n matrices and constants The group 2 constraints are achieved with the n×n matrices and constants for j = 1, . . . , n. To implement the rest of the constraints, we introduce d 2 Hermitian basis-operators B j for L(H q ) where L(H q ) denote the space of Hermitian operators in H q , Tr[B j B k ] = δ jk and B 1 proportional to the identity [74][75][76]. If S θ is not full rank, the number of basis operators can be reduced by (d − r) 2 where r is the rank of S θ by restricting B j to the quotient space L(H q )/L(ker(S θ )). See for example the discussions in [52,Sec. 2.10] or [53]. The group 3 constraints are then implemented by n×d 2 matrices and constants for j = 1, . . . , d 2 . The group 4 constraints are implemented with n 2 − n 2 ×d 2 matrices and constants 1,2,j = 0 , 1,3,j = 0 , 2,3,j = 0 . (D7) for j = 1, . . . , d 2 . Finally, the group 5 constraints are implemented with d 2 matrices and constants for j = 2, 3, . . . , d 2 .
The worst-case time complexity for solving the SDP (D1) or (D2) to a desired accuracy is O( √ N log(1/ )) where N = (n + 1)d is the size of the matrix F 0 [65,77]. However in our simulations, we observed that the time complexity is independent of N . This is consistent with reports in the literature that in practice, the SDP algorithms perform much better than its worst-case bound [77]. Each time step requires solving a system of linear equations with a computational complexity of O(N 3 ). Therefore, the overall worst-case computational complexity is O N 3/2 log(1/ ) .
Appendix E: Estimation of qubit rotations under phase damping channel with a two-qubit probe-analytic POVM saturating the Nagaoka-Hayashi bound We now present an analytic measurement strategy that saturates the Nagaoka-Hayashi bound for the qubit rotation estimation problem. We first define the four sub-normalised projectors where a and b are two non-zero real parameters satisfying a 2 +b 2 ≤ 1. An optimal strategy that saturates the Nagaoka bound for estimating θ x and θ y consists of measuring the five-outcome POVM with Π j = |φ j φ j | for j = 1, 2, 3, 4 and The probability for each POVM outcome is We can use this to construct unbiased estimators for θ x and θ y with In this construction, the fifth outcome Π 5 does not give any additional information about θ x or θ y . Nonetheless, it is still necessary to be included so that the POVM outcomes sum up to 1. For a finite sample, to have a better estimate of θ x and θ y , it is thus beneficial to have both a and b large so the outcomes Π 1 to Π 4 occur more often. However, in the asymptotic limit, the variances in our estimate of θ x and θ y are which do not depend on a or b. The sum v x + v y = 4/(2 − ) saturates the Nagaoka bound as claimed.
For estimating all three parameters θ x , θ y and θ z , one measurement strategy is to use the same POVM outcomes for estimating θ x and θ y but splitting Π 5 to get some information on θ z . Ideally, we would like to use these four projectors we get when setting a = b = 0, to obtain the most information on θ z without affecting the estimate of θ x and θ y . But the problem is that at this singular point, the first four outcomes Π 1 , Π 2 , Π 3 and Π 4 do not give any information on θ x and θ y . To fix this, we need both a and b to be close to but not exactly zero. Writing δ = (a 2 + b 2 )/2, we can split Π 5 as which has outcome probabilities This together with ξ z,1 = ξ z,2 = ξ z,3 = ξ z,4 = ξ z,5 = 0 , and ξ z,6 = −ξ z,7 = 1 give a variance for estimating θ z as v z = 1 which approaches v z = 1 (1 − ) 2 as δ tends to zero.
whose derivatives evaluated at φ = 0 are where the matrix basis is {|00 , |01 , |10 }. The Holevo bound for this model was computed by Albarelli et al. [53] to be otherwise. (F3) In the following, we show that this bound can be saturated by a separable measurement. There exist a family of measurements that can saturate the Holevo bound. One of them is the four-outcome POVM together with the estimation coefficients One can verify that when η < (a 2 0 −a 2 1 )/2, these outcomes are non-negative operators that satisfy Π 1 +Π 2 +Π 3 +Π 4 = 1. The estimator matrices X η = ξ η,1 Π 1 + ξ η,2 Π 2 + ξ η,3 Π 3 + ξ η,4 Π 4 and X φ = ξ φ,3 Π 3 + ξ φ,4 Π 4 satisfy the unbiased conditions, (4) in the main text. The probability for each outcome to occur is The variances of these two estimators are which together gives v η + v φ = (1 + 3η − 4η 3 )/4ηa 2 1 saturating the Holevo bound (F3) as claimed. At the boundary η = (a 2 0 − a 2 1 )/2a 2 0 , the POVM outcome Π 2 = 0 while the remaining three reduce to a projective measurement on the eigenstate of the SLD operator [53] In this case, the estimator coefficients are This measurement scheme remains optimal even when η > (a 2 0 − a 2 1 )/2a 2 0 . Comparing the 4-outcome POVM (F4) to the 3-outcome POVM (F8), we see that the role played by Π 2 is to obtain a better estimate of η, but at the expense of a worse estimate of φ. Whether this trade-off improves the overall sum of the MSE depends on the exact form of the probe and the value of η. We note that the estimators presented here depend on the unknown parameter η.
Although this would be an issue if we were interested in global parameter estimation, for local estimation this is not an issue, as we are only interested in estimating η in the local neighbourhood of some a priori known value, η 0 .

Appendix G: Dual solutions for the semidefinite program
From the constructed POVM, we can arrive at a candidate for the optimal X and L matrices using (17) and (14) from the main text which gives an upper bound to the primal solution. In this appendix, we write down the dual problem and provide its solution which gives a lower bound to the primal solution. One can easily check that the lower and upper bounds coincide which implies that the candidate solution is indeed an optimal solution for the Nagaoka-Hayashi bound.
The dual problem isc where the matrices F k and constants c k implements the constraints on the primal SDP as defined in appendix D.

Phase and transmissivity estimation in interferometer-dual solutions
We now write down the dual solution to the Nagaoka-Hayashi bound for the second example, phase and transmissivity estimation in an interferometer, when N = 1. To do this, we use the following 9 matrices as basis matrices: and all other y k zero. When the condition a 1 < 1/ √ 2 and η < (a 2 0 − a 2 1 )/2a 2 0 is not satisfied, one solution to the dual problem is given by: and all other y k zero. One can check that these solutions coincide with the primal solution in appendix F.