{'fraction_non_alphanumeric': 0.04711674819547794, 'fraction_numerical': 0.015264742905400022, 'mean_word_length': 4.058842495566661, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Experiment Design for Quantum State and Process Tomography and Hamiltonian Parameter Estimation

A number of problems in quantum state and system identification are addressed. Specifically, it is shown that the maximum likelihood estimation (MLE) approach, already known to apply to quantum state tomography, is also applicable to quantum process tomography (estimating the Kraus operator sum representation (OSR)), Hamiltonian parameter estimation, and the related problems of state and process (OSR) distribution estimation. Except for Hamiltonian parameter estimation, the other MLE problems are formally of the same type of convex optimization problem and therefore can be solved very efficiently to within any desired accuracy. Associated with each of these estimation problems, and the focus of the paper, is an optimal experiment design (OED) problem invoked by the Cramer-Rao Inequality: find the number of experiments to be performed in a particular system configuration to maximize estimation accuracy; a configuration being any number of combinations of sample times, hardware settings, prepared initial states, etc. We show that in all of the estimation problems, including Hamiltonian parameter estimation, the optimal experiment design can be obtained by solving a convex optimization problem. Software to solve the MLE and OED convex optimization problems is available upon request from the first author.


Introduction
"In a machine such as this [a quantum computer] there are very many other problems due to imperfections. For example, in the registers for holding the data, there will be problems of cross-talk, interactions between one atom and another in that register, or interaction of the atoms in that register directly with things that are happening along the program line that we didn't exactly bargain for. In other words, there may be small terms in the Hamiltonian besides the ones we've written. Until we propose a complete implementation of this, it is very difficult to analyze. At least some of these problems can be remedied in the usual way by techniques such as error correcting codes and so forth, that have been studied in normal computers. But until we find a specific implementation for this computer, I do not know how to proceed to analyze these effects. However, it appears that they would be very important in practice. This computer seems to be very delicate and these imperfections may produce considerable havoc." -Richard P. Feynman, "Quantum Mechanical Computers," Optics News, February 1985.

Alleviating the "havoc"
The concerns heralded by Feynman remain of concern today in all the implementations envisioned for quantum information systems. In a quantum computer it is highly likely that in order to achieve the desired system objectives, these systems will have to be tuned, or even entirely determined, using estimated quantities obtained from data from the actual system rather than solely relying on an initial design from a theoretical model. The problem addressed here is to design the experiment in order to yield the optimum information for the intended purpose. This goal is not just limited to quantum information systems. It is an essential step in the engineering practice of system identification [22,Ch.14]. That is, the design of the experiment which gives the best performance against a given set of criteria, subject to constraints reflecting the underlying properties of the system dynamics and/or costs associated with the implementation of certain operations or controls.
Clearly each application has a specific threshold of performance. For example, the requirements in quantum chemistry are generally not as severe as in quantum information systems. The objective of a measurement, therefore, depends on the way in which information is encoded into the system to begin with, and this is in turn, depends on the application. In this paper we are concerned with estimating quantum system properties: the state, the process which transforms the state, and parameters in a Hamiltonian model.
The estimation of the state of a quantum system from available measurements is generally referred to as quantum state tomography about which there is extensive literature on both theoretical and experimental aspects, e.g., see [27,Ch.8], [15] and the references therein. The more encompassing procedure of quantum system identification is not so easily categorized as the nomenclature (and methodology) seems to depend on the type and/or intended use of the identified model. For example, quantum process tomography (QPT) refers to determining the Kraus operator-sum-representation (OSR) of the input state to output state (completely positive) map, e.g., [27, §8.4.2], [7]. Hamiltonian parameter estimation refers to determining parameters in a model of the system Hamiltonian, e.g., [24], [6], [12], [37]. Somewhere in between quantum process tomography and Hamiltonian parameter estimation is mechanism identification which seeks an estimate of population transfer between states as the system evolves, e.g., [25].
Maximum likelihood estimation (MLE), a well established method of parameter estimation which is used extensively in current engineering applications, e.g., [22], was proposed in [4,29] and [33] for quantum state tomography of a quantum system with non-continuing measurements, i.e., data is taken from repeated identical experiments. Also, as observed in [29,33], the MLE of the density matrix is a convex optimization problem.
In this paper we address the related problem of optimal experiment design (OED) so as to secure an estimate of the best quality. The approach presented relies on minimizing the Cramer-Rao lower bound [8] where the design parameters are the number of experiments to be performed while the system is in a specified configuration. En route we also show that many related problems in state and process tomography can also be solved using MLE, and moreover, they are all formally the same type of convex optimization problem, namely, a determinant maximization problem, referred to as a maxdet problem [5,36]. Similarly, the OED problem posed here is also of a single general type of convex optimization problem, namely, a semidefinite program (SDP).
Convexity arises in many ways in quantum mechanics and this is briefly discussed in §1.2. The great advantage of convex optimization is a globally optimal solution can be found efficiently and reliably, and perhaps most importantly, can be computed to within any desired accuracy. Achieving these advantages, however, requires the use of specialized numerical solvers. As described in §1.3, the appropriate convex solvers have been embedded in some software tools we have composed which can solve the MLE and OED problems presented here.
In the remainder of the paper we present both MLE and the corresponding OED as applied to: quantum state tomography (MLE in §2.2 and OED in §2.3), estimating the distribution of known input states (MLE in §2.5 and OED in §2.6), quantum process tomography using the Kraus operator sum representation (MLE in §3.1 and OED in §3.2), estimating the distribution of a known OSR set (MLE in §3. 4 and OED in §3. 5), and to Hamiltonian parameter estimation (MLE in §4.1 and OED in §4.2). A summary in table form is presented in §5 followed by a discussion in §6 of the relation of MLE and OED to iterative adaptive control of quantum systems.

Convexity and quantum mechanics
Many quantum operations form convex sets or functions. Consider, for example, the following convex sets which arise from some of the basic aspects of quantum mechanics: probability outcomes {p α ∈ R} α p α = 1, p α ≥ 0 density matrix {ρ ∈ C n×n } Tr ρ = 1, ρ ≥ 0 positive operator valued measure (POVM) operator sum representation (OSR) in fixed basis { B i ∈ C n×n | i = 1, . . . , n 2 } X ∈ C n 2 ×n 2 ij X ij B * i B j = I n , X ≥ 0 An example of a convex function relevant to quantum information is worst-case gate fidelity, a measure of the "distance" between two unitary operations on the same input. As pointed out in [13], there are many ways to define this measure. Consider, for example, where U des ∈ C n×n is the desired unitary and U act ∈ C n×n is the actual unitary. In this case the worst-case fidelity can be interpreted as the minimum probability of obtaining the desired output state U des ψ over all possible pure input states ψ which produce the actual output state U act ψ. If U des and U act differ by a scalar phase then the worst-case fidelity is clearly unity; which is consistent with the fact that a scalar phase cannot be measured. This is not the case for the error norm U des − U act . As shown in Appendix §A.1, obtaining the worst-case fidelity requires solving the following (convex) quadratic programming (QP) problem: with the vectors a, b in R n the real and imaginary parts, respectively, of the eigenvalues of the unitary matrix U * des U act , that is, a = Re eig(U * des U act ), b = Im eig(U * des U act ). In some cases it is possible to compute the worst-case fidelity directly, e.g., in the example in Section §4.3 and in some examples in [27, §9.3]. Although the optimal objective value f wc (U des , U act ) is global, the optimal worst-case state which achieves this value is not unique.
In addition to these examples, convex optimization has been exploited in [3] and [33] in an attempt to realize quantum devices with certain properties. In [9] and [21], convex optimization is used to design optimal state detectors which have the maximum efficiency.
In general, convex optimization problems enjoy many useful properties. From the introduction in [5], and as already stated, the solution to a convex optimization problem can be obtained to within any desired accuracy. In addition, computation time does not explode with problem size, stopping criteria always produce a lower bound on the solution, and if no solution can be found a proof of infeasibility is provided. There is also a complete duality theory which can yield more efficient computation as well as optimality conditions. This is explored briefly in Section §2.3.

Software for tomography & experiment design
We have composed some MATLAB m-files which can be used to solve a subset of the QPT and OED convex optimization problems presented here. The examples shown here were generated using this software. The software, available upon request from the first author, requires the convex solvers YALMIP [23] and SDPT3 [35] which can be downloaded from the internet. These solvers make use of interior-point methods for solving convex optimization problems, e.g., [5,Ch.11], [26]. by the index γ, γ = 1, . . . , n cfg . Configurations can include wave-plate angles for photon counting, sample times at which measurements are made, and settings of any experimental "knobs" such as external control variables, e.g., laser wave shape parameters, magnetic field strengths, and so on. For quantum process tomography ( §3.1) and Hamiltonian parameter estimation ( §4.1), configurations can also include distinctly prepared initial states.
The problem addressed in this section is to determine the minimum number of experiments per configuration in order to obtain a state estimate of a specified quality, i.e., what is the tradeoff between number of experiments per configuration and estimation quality. The method used to solve this problem is based on minimizing the size of the Cramér-Rao lower bound on the estimation error [8].

Data collection
The data is collected using a procedure referred to here as non-continuing measurements. Measurements are recorded from identical experiments in each configuration γ repeated ℓ γ times. The set-up for data collection is shown schematically in Figure 1 for configuration γ.
Outcome counts n αγ , ℓ γ trials α = 1, . . . , n out Here ρ true ∈ C n×n is the true, unknown state to be estimated, σ true γ ∈ C n×n is the reduced density matrix which captures all the statistical behavior of the Q-system under the action of the measurement apparatus, and n αγ is the number of times outcome α is obtained from the ℓ γ experiments. Thus, where ℓ expt is the total number of experiments. The data set consists of all the outcome counts, The design variables used to optimize the experiment are the non-negative integers {ℓ γ } represented by the vector, Let p true αγ denote the true probability of obtaining outcome α when the system is in configuration γ with state input ρ true . Thus, E n αγ = ℓ γ p true αγ (6) where the expectation E(·) taken with repect to the underlying quantum probability distributions.
We pose the following model of the system, where p αγ (ρ) is the outcome probability of measuring α when the system is in configuration γ with input state ρ belonging to the set of density matrices, ρ ∈ C n×n | ρ ≥ 0, Tr ρ = 1 (8) {M αγ } are the POVM elements of the measurement apparatus, and thus, for γ = 1, . . . , n cfg , α M αγ = I n , M αγ ≥ 0, α = 1, . . . , n out (9) and σ γ (ρ) is the reduced density output state of the Q-system model. A general (model) representation of the Q system is the Kraus operator-sum-representation (OSR) which can account for many forms of error sources as well as decoherence [27]. Specifically, in configuration γ, the Q-system model can be parametrized by the set of Kraus matrices, K γ = { K γk ∈ C n×n | k = 1, . . . , κ γ } as follows: with κ γ ≤ n 2 . Implicit in this OSR is the assumption that the Q-system is trace preserving.
Combining this with the measurement model (9) gives the model probability outcomes, In this model, the outcome probabilities are linear in the input state density matrix. Moreover, the set O γ = { O αγ | α = 1, . . . , n out }, satisfies (9), and hence, is a POVM. 2 If the Q-system is modeled as a unitary system, then, The set O γ is still a POVM; in effect the OSR has a single element, namely,

System in the model set
We make the following assumption throughout: the true system is in the model set. This means that, This is always a questionable assumption and in most engineering practice is never true. Relaxing this assumption is an active research topic particularly when identification (state or process) is to be used for control design, e.g., see [19] and [34]. The case when the system is not in the model set will not be explored any further here except for the effect of measurement noise which is discussed next. It is important to emphasize that in order to produce an accurate unbiased estimate of the true density it is necessary to know the noise elements (as described next) which is a consequence of assumption (13).

Noisy measurements
Sensor noise can engender more noisy outcomes than noise-free outcomes. Consider, for example, a photon detection device with two photon-counting detectors. If both are noise-free, meaning, perfect efficiency and no dark count probability, then, provided one photon is always present at the input of the device, there are only two possible outcomes: {10, 01}. If, however, each detector is noisy, then either or both detectors can misfire or fire even with a photon always present at the input. Thus in the noisy case there are four possible outcomes: {10, 01, 11, 00}. Let { M αγ | α = 1, . . . , n out } denote the noisy POVM and let M αγ | α = 1, . . . , n out denote the noise-free POVM with n out ≥ n out where, The {ν αβγ } represents the noise in the measurement, specifically, the conditional probability that α is measured given the noise-free outcome β with the system in configuration γ. Since α ν αβγ = 1, ∀β, γ, it follows that if the noise-free set is a POVM then so is the noisy set.

Maximum likelihood state estimation
The Maximum Likelihood (ML) approach to quantum state estimation presented in this section, as well as observing that the estimation is convex, can be found in [29], [37] and the references therein. Using convex programming methods, such as an interior-point algorithm for computation, was not exploited in these references.
If the experiments are independent, then the probability of obtaining the data (4) is a product of the individual model probabilities (7). Consequently, for an assumed initial state ρ, the model predicts that the probability of obtaining the data set (4) is given by, The data is thus captured in the outcome counts {n αγ } whereas the model terms have a ρdependence. The function Prob {D, ρ} is called the likelihood function and since it is positive, the maximum likelihood estimate (MLE) of ρ is obtained by finding a ρ in the set (8) which maximizes the log-likelihood function, or equivalently, minimizes the negative log-likelihood function, These expressions are obtained by combining (15), (16) and (11). The Maximum Likelihood state estimate, ρ ML , is obtained as the solution to the optimization problem: minimize L(D, ρ) = − α,γ n αγ log Tr O αγ ρ subject to ρ ≥ 0, Tr ρ = 1 L(D, ρ) is a positively weighted sum of log-convex functions of ρ, and hence, is a log-convex function of ρ. The constraint that ρ is a density matrix forms a convex set in ρ. Hence, (17) is in a category of a class of well studied log-convex optimization problems, e.g., [5].

Pure state estimation
Suppose it is known the the true state is pure, that is, ρ true = ψ true ψ * true with ψ true ∈ C n and ψ * true ψ true = 1. In practice we have found that solving (17) when the true state is pure gives solutions which are easily approximated by pure states, that is, the estimated state has one singular value near one and all the rest are very small and positive.
To deal directly with pure state estimation we first need to characterize the set of all density matrices which are pure. This is given by the set { ρ ∈ C n×n | ρ ≥ 0, rank ρ = 1 }, which is equivalent to, ρ ∈ C n×n ρ ≥ 0, Tr ρ = 1, Tr ρ 2 = 1 The corresponding ML estimate is then the solution of, minimize L(ρ) = − α,γ n αγ log Tr O αγ ρ subject to ρ ≥ 0, Tr ρ = 1, Tr This is not a convex optimization problem because the equality constraint, Tr ρ 2 = 1, is not convex. However, relaxing this constraint to the convex inequality constraint, Tr ρ 2 ≤ 1, results in the convex optimization problem: If the solution is on the boundary of the set Tr ρ 2 ≤ 1, then a pure state has been found. There is however, no guaranty that this will occur.

Least-squares (LS) state estimation
In a typical application the number of trials per configuration, ℓ γ , is sufficiently large so that the empirical estimate of the outcome probability, is a good estimate of the true outcome probability p true αγ . The empirical probability estimate also provides the smallest possible value of the negative log-likelihood function, that is, p emp αγ is the solution to, minimize L(p) = − α,γ n αγ log p αγ subject to α p αγ = 1, ∀γ, p αγ ≥ 0, ∀α, γ with optimization variables p αγ , ∀α, γ. Thus, for any value of ρ we have the lower bound, In particular, assuming (6) holds, and the ℓ γ trials are independent, then the variance of the empirical estimate is known to be [28], It therefore follows that for large ℓ γ , p emp αγ ≈ p true αγ , and if as assumed (13), the system is in the model set, then p true αγ = Tr O αγ ρ true . These two conditions lead to taking the state estimate as the solution to the constrained weighted least-squares problem: The weights, w γ , are chosen by the user to emphasize different configurations. A typical choice is the distribution of experiments per configuration, hence, w γ ≥ 0, γ w γ = 1. Because of the semi-definite constraint, this weighted-least-squares problem is a convex optimization in the variable ρ. For large ℓ γ , the solution ought to be a good estimate of the true state. There is, however, little numerical benefit in solving (25) as compared to (17) -they are both convex optimization problems and the numerical complexity is similar provided (17) is solved using an interior-point method [5]. Some advantage is obtained by dropping the semidefinite constraint ρ ≥ 0 in (25) resulting in, This is a standard least-squares problem with a linear equality constraint which can be solved very efficiently using a singular value decomposition to eliminate the equality constraint [14]. For sufficiently large ℓ γ the resulting estimate may satisfy the positivity constraint ρ ≥ 0. If not, it is usually the case that some of the small eigenvalues of the state estimate or estimated outcome probabilities are slightly negative which can be manually set to zero. Solving (26) is numerically faster than solving (17), but not by much. Even with a large amount of data the solution to (26) can produce estimates which are not positive if the data is not sufficiently rich.
In this case the estimates from any procedure which eliminates the positivity constraint can be very misleading. It thus appears that even for large ℓ γ , there is no significant benefit accrued, either because of numerical precision or speed, to using the empirical estimate followed by standard leastsquares. If, however, the ℓ γ are not sufficiently large and/or the data is not sufficiently rich, then it is unlikely that the estimate from (26) will be accurate.
One possible advantage does come about because the solution to (26) can be expressed analytically, and thus it is possible to gain an understanding of how to select the POVM. For example, in [27] special POVM elements are selected to essentially diagonalize the problem, thereby making the least-squares problem (26) simpler, i.e., the elements of the density matrix can be estimated one at a time. However, implementing the requisite POVM set may be very difficult depending on the physical apparatus involved.

Experiment design for state estimation
In this section we describe the experiment design problem for quantum state estimation. The objective is to select the number of experiments per configuration, the elements of the vector ℓ = [ℓ 1 · · · ℓ n cfg ] T ∈ R n cfg , so as to minimize the error between the state estimate, ρ(ℓ), and the true state ρ true . Specifically, we would like to solve for ℓ from: minimize E ρ(ℓ) − ρ true 2 frob subject to γ ℓ γ = ℓ expt integer ℓ γ ≥ 0, γ = 1, . . . , n cfg (27) where ℓ expt is the desired number of total experiments. This is a difficult, if not insoluble problem for several reasons. First, the solution depends on the estimation method which produces ρ(ℓ). Secondly, the problem is integer combinatorial because ℓ is a vector of integers. And finally, the solution depends on ρ true , the very state to be estimated. Fortunately all these issues can be circumvented.
We first eliminate the dependence on the estimation method. The following result can be established using the Cramér-Rao Inequality [8]. The derivation is in Appendix §A.3. 3 Suppose the system generating the data is in the model set used for estimation, i.e., (13) holds. For ℓ = [ℓ 1 · · · ℓ n cfg ] experiments per configuration, suppose ρ(ℓ) is a density matrix and an unbiased estimate of ρ true , i.e., ρ(ℓ) ≥ 0, Tr ρ(ℓ) = 1, and E ρ(ℓ) = ρ true . Under these conditions, the estimation error variance satisfies,

State estimation variance lower bound
and C eq ∈ R n 2 ×n 2 −1 is part of the unitary matrix in the singular value decomposition: vec This theorem states the for any unbiased estimate of ρ true , the variance of the estimate satisfies the inequality (28). The power of the result is that it is independent of how the estimate is obtained, i.e., no estimation algorithm which produces an unbiased estimate can have an estimation error variance smaller than that in (28). There is a generalization for biased estimators but we will not pursue that here.
In general it is difficult to determine if any estimate will achieve the lower bound. However, under the conditions stated in the above result, the ML estimate, ρ ML (ℓ), the solution to (17), approaches ρ true with probability one, asymptotically as ℓ expt increases, and the asymptotic distribution becomes Gaussian with covariance given by the Cramér-Rao bound (see §A.2 for the covariance expression and [22] for a derivation).
The one qualifier to the Cramér-Rao bound as presented is that the indicated inverse exists. This condition, however, is necessary and sufficient to insure that the state is identifiable. More precisely, the state is identifiable if and only if, Under the condition of identifiability, the experiment design problem can be expressed by the following optimization problem in the vector of integers ℓ: where ℓ expt is the desired number of total experiments. The good news is that the objective, Unfortunately, there are still two impediments: (i) restricting ℓ to a vector of integers makes the problem combinatorial; (ii) the lower-bound function V (ℓ, ρ true ) depends on the true value, ρ true . These difficulties can be alleviated to some extent. For (i) we can use the convex relaxation described in [5, §7,5]. For (ii) we can solve the relaxed experiment design problem with either a set of "what-if" estimates as surrogates for ρ true , or use nominal values to start and then "bootstrap" to more precise values by iterating between state estimation and experiment design. We now explain how to perform these steps.

Relaxed experiment design for state estimation
Following the procedure in [5, §7.5], introduce the variables λ γ = ℓ γ /ℓ expt , each of which is the fraction of the total number of experiments performed in configuration γ. Since all the ℓ γ and ℓ expt are non-negative integers, each λ γ is non-negative and rational, specifically an integer multiple of 1/ℓ expt , and in addition, γ λ γ = 1. Let ρ surr denote a surrogate for ρ true , e.g., an estimate or candidate value of ρ true . Using (28)- (29) gives, Using (29), Hence, the objective function V (ℓ, ρ surr ) can be replaced with V (λ, ρ surr ) and the experiment design problem (31) is equivalent to.
The objective is now a convex function of the λ γ , but it is still a combinatorial problem because the λ γ are constrained to each be an integer multiple of 1/ℓ expt . If λ γ is only otherwise constrained to the non-negative reals, then this has the effect of relaxing the constraint that the ℓ γ are integers. As phrased in [5], the relaxed experiment design problem is: The objective is convex, the equality constraint is linear, and the inequality constrains are convex, hence, this is a convex optimization problem in λ ∈ R n cfg . Let λ opt denote the optimal solution to (35). Since the problem no longer depends on ℓ expt , λ opt can be viewed as a distribution of experiments per configuration. 5 Clearly there is no guaranty that ℓ expt λ opt is a vector of integer multiples of 1/ℓ expt . A practical choice for obtaining a vector of integer multiples of 1/ℓ expt is, If ℓ opt is the (unknown) integer vector solution to (31), then we have the relations: The optimal objective is thus bounded above and below by known values obtained from the relaxed optimization. The gap within which falls the optimal solution can be no worse than the difference between V (ℓ round expt , ρ surr ) and V (ℓ expt λ opt , ρ surr ), which can be computed solely from λ opt . If the gap is sufficiently small then for all practical purposes the "optimal" solution is λ opt . From now on we will refer to λ opt as the optimal solution rather than the relaxed optimal.

Performance tradeoff
The optimal distribution λ opt can be used to guide the elimination of small values of λ opt . For example, consider the suboptimal distribution, λ sub , obtained by selecting the largest n sub out of n cfg non-zero values of λ opt . Let ℓ sub expt denote the integer vector of configurations, Using (32), the minimum number of experiments so that V (ℓ sub expt , ρ surr ) ≤ V 0 is given by, As n sub is varied, the graph {n sub , ℓ sub expt } establishes a tradeoff between the number of configurations per experiment versus the total number of experiments such that the lower bound on the estimation variance does not exceed the desired value V 0 . When n sub = n cfg , (39) is identical with (36).
The condition number of the matrix G(λ opt , ρ surr ) gives an indication of the identifiability of the density matrix ρ surr . A very large condition number means that the linear combination of elements of the density matrix associated with the small eigenvalue will be more difficult to obtain then those combinations associated with a large eigenvalue. The condition number of G(λ opt , ρ surr ) is not only affected by the number of experiments per configuration, λ opt , but by the configurations themselves. Examining G(λ opt , ρ surr ) for different ρ surr (surrogates of ρ true ) and different configurations can help establish a good experiment design.

Bootstrapping
A standard approach used to circumvent not knowing the true state needed to optimize the experiment design is to proceed adaptively, or by "bootstrapping." The idea is to use the current estimate of the initial state found from (17), then solve (35), and then repeat. The algorithm at the k-th iteration looks like this: The initial distribution ℓ(0) could be chosen as uniform, e.g., the same for a not too large number of configurations. The algorithm could also start by first solving for a distribution from an initial state surrogate. In each iteration we could also vary ℓ expt . Although each optimization is convex, the joint problem may not be. Conditions for convergence would need to be investigated as well as establishing that this method is efficient, i.e., reduces the number of trials. We will not pursue this any further here.

Dual experiment design problem
Lagrange Duality Theory can provide a lower bound on the objective function in an optimization problem as well as establishing optimality conditions often leading to insights into the optimal solution structure [5,Ch.5]. In many cases the largest lower bound -the solution of the dual problem -is equal to the optimal objective function. The dual problem associated with the experiment design problem (35) is, The optimization variable is W ∈ C n 2 −1×n 2 −1 . The above form of the dual is given in [5, §7.5.2] for a slightly simpler problem (the G γ are dyads) but is essentially the same. A key observation arises from the complementary slackness condition, where λ opt is the solution to the primal problem, (35), and W opt is the solution to the dual problem, (41). Thus, only when the equality constraint holds, Tr W opt G γ (ρ surr ) = 1, is the associated λ opt γ not necessarily equal to zero. It will therefore be usually the case that many of the elements of the optimal distribution will be zero.
Strong duality also holds for this problem, thus the optimal primal and dual objective values are equal, For this problem, a pair (λ, W ) is optimal with respect to ρ surr if and only if:

Example: experiment design for state estimation
A schematic of an apparatus for state tomography of a pair of entangled photons specified by the quantum state (density matrix) ρ is shown in figure 2. The set up has four photon-counting detectors, A, B, C, D. There are four continuous variable settings for the quarter-wave plates and half-wave plates, i.e., q, h, q ′ , h ′ . For any settings of these parameters one of the detectors in each arm will register a photon. The objective is to determine the optimal settings of these parameters and the number of experiments per setting for estimation of the state ρ of the pair using as data the photon counts from the four detectors.
Because the photon sources are not completely efficient, the input quantum state actually consists of either two or zero photons. The detectors register a 0 or 1 depending on whether a photon is incident on them or not. The basis states for the upper arm are therefore: There is a similar set for the lower arm (modes g, h).
The firing patterns for an arbitrary setting of the wave plates, assuming perfect detection efficiency and no dark counts are given in the table: The probabilities for these patterns are given by where {i, j, k, ℓ} ∈ {0, 1}, and M ij AB is the projector for detector A to register count i and simultaneously detector B to register count j. Similarly, M kℓ CD is the projector for detector C to register count j and simultaneously detector D to register count ℓ. The projectors for A and B in the above basis are: A similar set of projectors can be written for C and D with the variables h, q replaced by their primed counterparts h ′ , q ′ . The protocol is to measure the probabilities for enough settings of the variables that the elements of the two-photon density operator can be estimated. The two-photon density operator is the direct product of the one-photon density operator, for which the set of 3 states given above forms a basis. The basis states of the two-photon (9 × 9) density operator, ρ, are: |ijkℓ = |i e |j f |k g |ℓ h with i, j, k, ℓ ∈ {0, 1} Hence, zero photon |0000 one photon |0100 , |1000 , |0001 , |0010 two photon |0101 , |0110 , |1001 , |1010 simulation results: one-arm Consider only one arm of the apparatus in figure 2, say the upper arm with detectors (A,B). Suppose the wave plate settings are, Assume also that the incoming state always is one photon, never none. Hence, ρ ∈ C 2×2 and the projectors are: with ψ 1 , ψ 2 from (46). Assuming each detector has efficiency η, 0 ≤ η ≤ 1 and a non-zero dark count probability, δ, 0 ≤ δ ≤ 1, then there are four possible outcomes at detectors A,B given in the following table: α A B 10 1 0 01 0 1 00 0 0 11 1 1 Following [15,39] the probability of a dark count is denoted by the conditional probability, where 1|0 means the detector has fired "1" given that no photon is present at the detector "0." As shown in [15], it therefore follows that the probability that the detector does not fire "0" although a photon is present at the detector "1" is given by, Here 1 − η is the probability of no detection and 1 − δ is the probability of no dark count. The remaining conditional probabilities are, by definition, constrained to obey: The probabilities for the firing patterns in the above table are thus given by (7) with the following observables M αγ : Numerical computer simulations were performed for two input state cases: pure state: ρ pure = 1 2 For each input state case we computed λ opt with and without "noise:" no noise detector efficiency η = 1 dark count probability δ = 0 yes noise detector efficiency η = 0.75 dark count probability δ = 0.05 For all cases and noise conditions we used the wave plate settings: Both angles are set from 0 to 45 • in 5 • increments. This yields a total of n cfg = 10 2 = 100 configurations corresponding to all the wave plate combinations. Figure 3 shows the optimal distributions λ opt versus configurations γ = 1, . . . , 100 for all four test cases: two input states with and without noise. Observe that the optimal distributions are not uniform but are concentrated near the same particular wave plate settings. These settings are very close to those established in [17].
To check the gap between the relaxed optimum λ opt and the unknown integer optimum we appeal to (36)- (37). The following table shows that these distributions are good approximation to the unknown optimal integer solution for even not so large ℓ expt for the two state cases with no noise. Similar results were obtained for the noisy case.  We explore the two arm case under the assumption that two photons are always present at the input, thereby excluding the zero photon case. 6 The table of detector firing patterns is: α A B C D 0101 0 1 0 1 0110 0 1 1 0 1001 1 0 0 1 1010 1 0 1 0 The following three 4 × 4 input states are considered: with ρ pure , ρ mixd given by (53). We use the angles from the one-arm optimal distribution for the largest This yields a total of n cfg = 4 4 = 256 wave plate settings. Figure 4 shows the optimal distributions λ opt versus configurations γ = 1, . . . , 256 for all three input states. In this case the optimal distributions λ opt are nearly uniform in magnitude. The following table compares the distributions for optimal, suboptimal with 25 angles, uniform at the 25 suboptimal angles, and uniform at all 256 angles by examining the minimum number of experiments required to obtain an RMS estimation error of no more than 0.01. As might be expected it is easier to estimate all pure states than mixed states. The table also shows that the optimal solution can be used effectively to guide the selection of suboptimal distributions. PSfrag replacements γ λ opt pure , no noise λ opt pure , yes noise λ opt mixd , no noise λ opt mixd , yes noise γ λ opt for ρ pure ⊗ ρ pure λ opt for ρ pure ⊗ ρ mixd λ opt for ρ mixd ⊗ ρ mixd Figure 4: optimal distributions -two-arms

Maximum likelihood state distribution estimation
A variation on the state estimation problem is to estimate the distribution of a known set of input states. The set-up for data collection is shown schematically in Figure 5 for configuration γ.
Outcome counts n αγ , ℓ γ trials α = 1, . . . , n out . . , f (n in )] T ∈ R n in is to be estimated. In this case the input state can be represented by the mixed state, The ML estimate of f is then the solution of the optimization problem: As in the MLE for quantum state estimation (17), the objective is log-convex in the state ρ(f ), the state is linear in f ∈ R n in , and the constraints form a convex set in f . Hence, this is a convex optimization problem in the variable f . Combining (59) with (60) gives the more explicit form, Here again as in (25) we could solve for f using the empirical estimate of the outcome probabilities as in (26):

Experiment design for state distribution estimation
Let f surr ∈ R n in be a surrogate for the true state distribution, f true . Following the derivation of (28) in Appendix §A.3, the associated (relaxed) optimal experiment design problem is, with a αγ ∈ R n in from (60) and where C eq ∈ R n in ×n in −1 is part of the unitary matrix W in the singular value decomposition: 1 T n in = USW T , W = [c C eq ] ∈ R n in ×n in .

Quantum Process Tomography: OSR Estimation
We explore two methods of identification for determining the Q-system from data: (i) in this section, estimating the OSR in a fixed basis, and (ii) in Section §4, estimating Hamiltonian parameters. In either case, the set-up is now as shown in Figure 6.
The main difference between state tomography ( Figure 1) and process tomography ( Figure  6) is that in the latter case the input state is prepared at specific values, ρ γ , depending on the configuration, whereas the Q-system does not depend on the configuration. If the process varies with every change in configuration it would be much more difficult to estimate; a model of configuration dependence would need to be established. This situation is perhaps amenable with Hamiltonian parameter estimation but will not be pursued here.

Maximum likelihood OSR estimation
As already stated, the Krause matrices for modeling the (trace preserving) Q-system in this case are not dependent on the configuration γ, specifically, K = { K k | k = 1, . . . , κ } with κ ≤ n 2 . Using (10), the reduced state model in Figure 6 as a function of K is, Combining the above with the measurement model (9) gives the probability outcomes model, The log-likelihood function (16) is, An ML estimate of K is then a solution to, This is not a convex optimization for two reasons: the equality constraint is not linear in K and the objective function is not convex. The problem can be transformed -more accurately, embedded -into a convex optimization problem by expanding the Kraus matrices in a fixed basis. The procedure, described in [27, §8.4.2], is as follows: since any matrix in C n×n can be represented by n 2 complex numbers, let be a basis for matrices in C n×n . The Kraus matrices can thus be expressed as, where the n 2 coefficients {a ki } are complex scalars. Introduce the matrix X ∈ C n 2 ×n 2 , often referred to as the superoperator, with elements, a * ki a kj , i, j = 1, . . . , n 2 As shown in [27], from the requirement to preserve probability, X is restricted to the convex set, The system output state (65) and outcome probabilities (66) now become, where the matrix R αγ ∈ C n 2 ×n 2 has elements, Quantum process tomography is then estimating X ∈ C n 2 ×n 2 from the data set D (4). An ML estimate is obtained by solving for X from: minimize L(D, X) = − α,γ n αγ log Tr XR αγ subject to X ≥ 0, This problem has essentially the same form as (17), and hence is also a convex optimization problem with the optimization variables being the elements of the matrix X. Since X = X * ∈ C n 2 ×n 2 , it can be parametrized by n 4 real variables. Accounting for the n 2 real linear equality constraints, the number of free (real) variables in X is thus n 4 −n 2 . This can be quite large even for a relatively small number of qubits, e.g., for q = [1, 2, 3, 4] qubits, n = 2 q = [2,4,8,16] and n 4 −n 2 = [12, 240, 4032, 65280]. This exponential (in qubit) growth is the main drawback to using this approach. The X (superoperator) matrix can be transformed back to Kraus operators via the singular value decomposition [27, §8.4.2]. Specifically, let X = V SV * with unitary V ∈ C n 2 ×n 2 and S = diag(s 1 · · · s n 2 ) with the singular values ordered so that s 1 ≥ s 2 ≥ · · · ≥ s n 2 ≥ 0. Then the coefficients in the basis representation of the Kraus matrices (70) are, Theoretically there can be fewer then n 2 Kraus operators. For example, if the Q system is unitary, then, In effect, there is one Kraus operator, U, which is unitary and of the same dimension as the input state ρ. The corresponding X matrix is a dyad, hence rank X = 1. A rank constraint is not convex. However, the X matrix is symmetric and positive semidefinite, so the heuristic from [10] applies where the rank constraint is replaced by the trace constraint, From the singular value decomposition of X, Tr X = k s k , and hence, adding the constraint (78) to (75) will force some (or many) of the s k to be small which can be eliminated (postoptimization) thereby reducing the rank. The auxiliary parameter η can be used to find a tradeoff between simpler realizations and performance. The estimation problem is then:

Experiment design for OSR estimation
Let X surr ∈ C n 2 ×n 2 be a surrogate for the true OSR, X true . As derived in Appendix §A.4, the associated (relaxed) optimal experiment design problem is, and C eq ∈ C n 4 ×n 4 −n 2 is part of the unitary matrix W = [C C eq ] ∈ C n 4 ×n 4 in the singular value decomposition of the n 2 × n 4 matrix, [a 1 · · · a n 4 ] = U √ nI n 2 0 n 2 ×n 4 −n 2 W * with with a k = vec(B * i B j ) ∈ C n 2 for k = i + (j − 1)n 2 , i, j = 1, . . . , n 2 . The columns of C eq , i.e., the last n 4 − n 2 columns of W , are a basis for the nullspace of [a 1 · · · a n 4 ].

Example: experiment design for OSR estimation
The 9 combinations of angles together with the 4 combinations of input states gives a total of 36 configurations, γ = 1, . . . , n cfg = 36. Figure 7 shows the optimal distribution of experiments for the 36 configurations using the true OSR corresponding to the Pauli basis set Since the system is simply the identity, Q(ρ) = ρ, with this basis choice, X true = diag(2 0 0 0). (No knowledge of the system being identity is used, hence, all elements of X true are estimated, not just the single element in the "11" location.) The following Approximately 35% fewer experiments are needed using the optimal distribution. Although not dramatic, as in the photon estimation example §2.4, there is a large penalty, in terms of time, for changing wave-plate angles. PSfrag replacements γ λ opt pure , no noise λ opt pure , yes noise λ opt mixd , no noise λ opt mixd , yes noise γ λ opt for ρ pure ⊗ ρ pure λ opt for ρ pure ⊗ ρ mixd λ opt for ρ mixd ⊗ ρ mixd γ λ opt

Maximum likelihood OSR distribution estimation
Suppose the Kraus matrices are known up to a scale factor which is related to its probability of occurrence, that is, One interpretation of this system model is that one of the matrices, say K 1 , is the nominal (unperturbed) system, and the others, K k , k = 2, . . . , κ, are perturbations, each of them occurring with probability q k . Examples of perturbations include the typical errors which can be handled by quantum error correction codes, e.g., depolarization, phase damping, phase and bit flip; see, e.g., [27,Ch.8].
The goal is to use the data to estimate the unknown vector of probabilities, q = [q 1 · · · q κ ] T ∈ R κ . Using the system model (83), the model probability outcomes are, The ML estimate of q ∈ R κ is the solution of the optimization problem, This is a convex optimization problem and is essentially in the same form as problem (60) which seeks the ML estimate of the input state distribution.

Experiment design for OSR distribution estimation
The formulation here is directly analogous to that of experiment design for state distribution estimation §2.6. Let q surr ∈ R n in be a surrogate for the true OSR distribution, q true . Following the lines of the derivation in Appendix §A.4, the associated (relaxed) optimal experiment design problem is, minimize V (λ, q surr ) = Tr with a αγ ∈ R n in from (85) and where C eq ∈ R n in ×n in −1 is part of the unitary matrix W in the singular value decomposition: 1 T n in = USW T , W = [c C eq ] ∈ R n in ×n in .

Example: experiment design for OSR distribution estimation
Consider a quantum process, or channel, where a single qubit state, ρ ∈ C 2×2 , is corrupted by a bit-flip error with occurrence probability q B and a depolarizing error with occurrence probability q D . The process is described by the quantum operation, 7 Q(ρ, q) = q I ρ + q B XρX + q D I/2 is the probability of no error occurring. The probability of observing outcome α with the system in configuration γ is, 8 An interesting aspect of this problem is that not all input states ρ γ lead to identifiability of the occurrence probabilities. And this is independent of the choice of POVM M αγ . To see this consider the single pure input state, The output of the channel (88) is then, Suppose we knew the elements of Q(ψψ * , q) perfectly; call them Q 11 , Q 12 , Q 22 . Then in principal we could solve for the three occurrence probabilities from the linear system of equations,    If det R = 0 then no unique solution exists; the occurrence probabilities are not identifiable. Specifically, det R = 0 for all a, b ∈ C such that, (Re ab * )(|b| 2 − |a| 2 ) = 0, |a| 2 + |b| 2 = 1 (93) 7 X is one of the three 2 × 2 Pauli spin matrices: [27, §8.3], an equivalent set of OSR elements which describe (88) are, Forming the probability outcomes in terms of this expansion results in an overparamtrization.

Hamiltonian Parameter Estimation
The process of modeling a quantum system in this case begins with the construction of a Hamiltonian operator on an infinite dimensional Hilbert space. Eventually, a finite dimensional approximation is invoked in order to calculate anything. (In some cases a finite dimensional model is immediately appropriate, e.g., spin systems, [11,.) The finite dimensional model is the starting point here.

Maximum likelihood Hamiltonian parameter estimation
The quantum system is modeled by a finite dimensional Hamiltonian matrix H(t, θ) ∈ C n×n , having a known dependence on time t, 0 ≤ t ≤ t f , and on an unknown parameter vector θ ∈ R n θ . The model density matrix will depend on θ and the initial (prepared and known) state drawn from the set of states ρ init β ∈ C n×n | β = 1, . . . , n in . Thus, the density matrix associated with initial state ρ init β is ρ β (t, θ) ∈ C n×n which evolves according to, where U(t, θ) ∈ C n×n is the unitary propagator associated with H(t, θ) which satisfies, At each of n sa sample times in a time interval of duration t f , measurements are recorded from identical repeated experiments. Specifically, let { t τ | τ = 1, . . . , n sa } denote the sample times relative to the start of each experiment. Let n αβτ be the number of times the outcome α is recorded at t τ with initial state ρ init β from ℓ βτ experiments. The data set thus consists of all the outcome counts, D = { n αβτ | α = 1, . . . , n out , β = 1, . . . , n in , τ = 1, . . . , n sa } The configurations previously enumerated and labeled by γ = 1, . . . , n cfg are in this case all the combinations of input states ρ init β and sample times τ , thus n cfg = n in n sa . For the POVM M α , the model outcome probability per configuration pair (ρ init β , t τ ) is, The Maximum Likelihood estimate, θ ML ∈ R n θ , is obtained as the solution to the optimization problem: where Θ is a set of constraints on θ. For example, it may be known that θ is restricted to a region near a nominal value, e.g., Θ = { θ | θ − θ nom ≤ δ }. Although this latter set is convex, unfortunately, the likelihood function, L(D, θ), is not guaranteed to be convex in θ. It is possible, however, that it is convex in the restricted region Θ, for example, if δ is sufficiently small.

Experiment design for Hamiltonian parameter estimation
Despite the fact that Hamiltonian parameter estimation is not convex, the (relaxed) experiment design problem is convex. A direct application of the Cramér-Rao bound to the likelihood function in (102) results in the following.
Hamiltonian parameter estimation variance lower bound Suppose the system generating the data is in the model set used for estimation, i.e., (13) holds. For ℓ = [ℓ 1 · · · ℓ n cfg ] experiments per configuration (ρ init β , t τ ), suppose θ(ℓ) ∈ R n θ is an unbiased estimate of θ true ∈ R n θ . Under these conditions, the estimation error variance satisfies, The relaxed experiment design problem with respect to the surrogate θ for θ true is, with optimization variables λ βτ , the distribution of experiments per configuration (ρ init β , t τ ). The difference between this and the previous formulation is that there are no equality constraints on the parameters. The gradient ∇ θ p αβτ (θ) and Jacobian ∇ θθ p αβτ (θ) are dependent on the parametric structure of the Hamiltonian H(t, θ).

Example: experiment design for Hamiltonian parameter estimation
Consider the system Hamiltonian, with constant control ε. The goal is to select the control to make the Hadamard logic gate, U had = (X + Z)/ √ 2. If θ true were known, then the control ε = 1/θ true would produce the Hadamard (to within a scalar phase) at time t = π/2, that is, We assume that only the estimate θ of θ true is available. Using the estimate and knowledge of the Hamiltonian model structure, the control is ε = 1/ θ. This yields the actual gate at t = π/2, Since the parameter estimate, θ, is a random variable, so is the normalized parameter error δ.
Assuming the estimate is unbiased, the expected value of the worst-case gate fidelity (1) is given explicitly by, Consider the case where the system is in the model set, the POVMs are projectors in the computational basis (|0 , |1 ), and the configurations consist of combinations of input states and sample times. Specifically, the example problem is as follows: sample times t k = δ(k − 1), k = 1, . . . , 100, δ = (π/2)/99 with pure input state In this example, with a single parameter and a single input state, the optimal experiment design problem (105) becomes: The trace operation in (105) is eliminated because the matrix G βτ (θ surr ) is now the scalar, The solution can be determined directly: concentrate all the experiments at the recording time t τ where g τ (θ surr ) is a maximum, specifically, The following tables show the minimum number of experiments at the optimal recording time, t opt , in order to achieve 0.01 accuracy (deviation) in the parameter estimate with θ surr = θ true . The two cases shown are for the two input states with the control set to unity.
To make the Hadamard the control update is ε = 1/ θ which follows from the (risky) assumption that the estimate, θ, is perfectly correct. For any of the true values from the above table (114), all the estimates have the same accuracy. Hence, the average value of the worst-case gate fidelity after the update is approximately, The ease of obtaining the estimate by minimizing the negative log-likelihood function can be determined by examining the average likelihood function, which is obtained from (102) and (6). Figure 8 shows plots of normalized 9 E L(θ), 0.8 ≤ θ ≤ 1.2 for the three true parameter values and corresponding optimal recording times from the table in (114) with control ε = 1 and initial state |0 . In all three cases, over the θ range shown, E L(θ) is convex. Figure 9 shows plots of normalized E L(θ), 0.8 ≤ θ ≤ 1.2 for the three true parameter values and corresponding optimal recording times from the table in (114) with control ε = 1 and initial state U had |0 . For initial state |0 , a range of 8876 to 13262 experiments are needed at the optimal recording time to achieve 0.01 deviation in the estimate. With the initial states U had |0 the same accuracy only requires about 2000 experiments. This difference can be inferred partly by comparing the curvature in the plots in Figure 8 with 9; as they are plotted on the same normalized scale. Note the increased curvature of E L(θ) in Figure 9 in the neighborhood of the true value.
If we further increase the control effort, say to ε = 5, the number of experiments required to achieve 0.01 accuracy is significantly reduced as seen in the following table. (117) 9 The plots show E L(θ) divided by its minimum value, thus normalized to have a minimum of unity.

Summarizing Maximum Likelihood Estimation & Optimal Experiment Design
The results presented show that an efficient numerical method based on convex programming is possible for optimizing the experiment for state and process tomography. In addition, the estimation of the state and/or process using data from non-continuing measurements is copacetic with Maximum Likelihood Estimation. Both the experiment design and estimation work naturally together and both can be solved using convex optimization methods.

Maximum likelihood estimation
The general form for estimating the parameter π true is obtained as the solution to the optimization problem, minimize L(π) = − α,γ n αγ log p αγ (π) Under the assumption that the system generating the data is in the model set used for estimation, the estimate, π ML , the solution to (118), is unbiased and has the asymptotic variance, where λ is the vector of fraction of experiments per configuration, and V (λ, π true ) is obtained from the Cramér-Rao inequality.

Optimal experiment design
The general form for estimating the configuration distribution λ is obtained as the solution to the optimization problem, where π is a surrogate for π true . Tables 1 and 2 summarize the class of estimation and experiment design problems, respectively. objective p αγ (π) π ∈ Π comment Hamiltonian parameter estimation p αγ a αγ a T αγ C eq Table 2: Summary of optimal experiment designs. These are convex optimization problems in all cases. The matrix C eq comes from the associated parameter equality constraints, e.g.,

Iterative Adaptive Control
Quantum process tomography can be used for adaptive control design as described in [20]. Adaptive control systems are in general one of two types [2,1]: (1) indirect adaptive controluse the data to first determine parameters in a system model, then based on the model determine the control parameters, and (2) direct adaptive control -use the data to directly select control parameters by comparing actual performance to an ideal. Applications of direct adaptive control (and learning principles) to quantum systems can be found in [30,41,31]. It could be argued that the direct approach is simpler in that a system model is not needed. However, this is a little deceptive because in effect the closed-loop system is being modeled indirectly via the ideal performance. This is made more explicit in the unfalsification/invalidation approach to adaptive control, [34], [19].

Indirect adaptive control
Hamiltonian parameter estimation and associated optimal experiment design can be combined in an iterative indirect adaptive control approach as depicted in Figure 11. control design ε (i) = arg opt ε J(ε, θ (i) ) experiment design ℓ (i) = round(ℓ expt λ (i) ) . . , n out , γ = 1, . . . , n cfg estimate parameters θ (i+1) = arg min L( The "opt" in the control design step could be to maximize worst-case gate fidelity (1), where U(t, ε, θ) is the propagator arising from the parametric Hamiltonian model. The control design step is not necessarily convex. Even if it were, it is not yet known under what conditions the complete iterative procedure will converge to the optimal control, or converge at all [20]. For example, in the simulations to follow convergence to the optimum is dependent on the initial parametrization. The properties of this type of iteration remain an area for further study.
function is obtained. The estimation is followed by a control using the estimated value in the pulse control table. In the two cases shown the algorithm converges to the true value.
Although not shown, the algorithm does not converge from all initial values of ω 1 . Figure 14 shows U(t f ( ω 1 ), ε( ω 1 )) − U des frob vs. estimate ω 1 /ω true 1 with the control from the table. The function is clearly not convex. The region of convergence for n sa = 1 and n sa = 2 sample times are shown in blue and green, respectively. The region of attraction is increased for n sa = 2. These results, of course, are specific to this example and can not be generalized. To re-iterate, conditions for convergence, region of attraction, and so on, are only partially understood, in general, for this type of iteration [16].

Direct adaptive control
In a direct adaptive control system no model is posed for the system; ideally only a performance measure is available and the control parameters are adjusted to improve the performance. The Figure 15 shows the landscape of the system with parameters, ω Qz = 1, ω Qx = 0.01, ω Ez = 1, ω Qx = 0, ω QE = 0.005, ρ init E = I 2 /2 and for constant values of the controls over the ranges 0.96 ≤ ε z ≤ 1.04, 0.1 ≤ ε x ≤ 5.2. where between these two extremes, ranging from having highly constrained knobs to the totally wrong knobs. This comment comes from the observation that physical effects spanning a dynamic range in quantum wavelength of about 10 7 are being controlled by a single type of laser, i.e., the Ti:Sapphire laser, working over a range of about 1-10 on that scale. That is, a domain in laser wavelength of 10 (in some units) is split up into about 100 small pieces as controls, and those very narrow controls manage everything over a huge dynamic range. The fact that the experiments work at all is rather amazing. A guess is that an examination of the landscapes will show considerable detail, much of which is likely false structures arising from having highly constrained controls. From a positive perspective, as more bandwidth becomes available the control landscape will become less complicated and more regular in the sense that more of the local optimum values will provide performance close to the global optimum.

A.2 Cramér-Rao Inequality
The following is the classical form of the Cramér-Rao Inequality.

Cramér-Rao Inequality[8]
Let θ 0 ∈ R p be the true parameter to be estimated from a data set D. Let L(D, θ 0 ) be the true negative log-likelihood function of the system generating the data. Let θ ∈ R p be an unbiased estimate of θ 0 , i.e., E θ = θ 0 . Then, the covariance of the estimate, cov θ = E θ − θ 0 θ − θ 0 T (123) satisfies the matrix inequality, where F (θ 0 ) is the Fisher information matrix, If F (θ 0 ) > 0, then (124) is equivalent to, This famous theorem states the for any unbiased estimator, the covariance of the estimate satisfies the inequality (124), or equivalently (126), provided the Fisher matrix is invertible. Usually only (126) is given as the theorem: the minimum covariance of the estimate is given by the inverse of the Fisher information matrix. The power of the result (124) is that it is independent of how the estimate is obtained. The lower bound only depends on the model structure, the experiment design, and the information in the data.

A.3 Derivation of (28)
Define the vectors r, a αγ ∈ C n 2 as, r = vec ρ, a αγ = vec O αγ Then p αγ = Tr O αγ ρ = a * αγ r and Tr ρ = 1 ⇔ b T r = 1, b = vec I n The next step eliminates the equality constraint b T r = 1 reducing the n 2 unknowns in r to n 2 − 1. Since b T b = n, the SVD of b is, b = W √ n 0 n 2 −n with unitary W ∈ R n 2 ×n 2 , Partition W = [c C eq ] with C eq ∈ R n 2 ×n 2 −1 . Then all r satisfying b T r = 1 are given by r = c/ √ n + C eq z for all z ∈ C n 2 −1 . The likelihood function with the equality constraint (Tr ρ = 1 or b T r = 1) eliminated is then a function only of z, L(D, z) = − α,γ n αγ log(c/ √ n + C eq z) To obtain the Cramér-Rao bound, we first compute, ∇ zz L(D, z) = α,γ n αγ p αγ (z) 2 (C T eq a αγ )(C T eq a αγ ) * Using (6), E n αγ = ℓ γ p αγ (ρ true ), ρ true = c/ √ n + C eq z true , and hence the Fisher information matrix is, with respect to z, F = E L(D, z = z true ) = α,γ ℓ γ p αγ (ρ true ) (C T eq a αγ )(C T eq a αγ ) * = G(ℓ, ρ true ) where the last line comes from the definition of G(ℓ, ρ true ) in (29). Let r = c/ √ n + C eq z be an unbiased estimate of r true . Thus, cov r = C eq cov z C T eq ≥ C eq F −1 C T eq Using var ρ = var r = Tr cov r and the fact that W is unitary, and hence, C T eq C eq = I n 2 −1 , we get, var ρ ≥ Tr C eq F −1 C T eq = Tr F −1 which is the final result (28)- (29).

A.4 Derivation of (80)
Define x, r αγ ∈ C n 4 as, x = vec X, r αγ = vec R αγ The likelihood function in (75) can then be written as, L(D, x) = − α,γ n αγ log r * αγ x and the equality constraint in (75) becomes, Ax = vec I n , A = [a 1 · · · a n 4 ] ∈ C n 2 ×n 4 with a k = vec(B * i B j ) ∈ C n 2 for k = i+(j −1)n 2 , i, j = 1, . . . , n 2 . Perform the singular value decomposition A = U[S 0]W * , W = [C C eq ] with C ∈ C n 2 ×n 2 , C eq ∈ C n 2 ×n 4 −n 2 as given in (82). From the definition of the basis functions (69) it follows that S = √ n I n 2 . Observe also that the columns of C eq are a basis for the nullspace of A. Hence all x satisfying the equality constraint are given by, x = x + C eq z, ∀z with x = (1/ √ n)CU * vec I n The likelihood function with the equality constraint (Ax = vec I n ) eliminated is then a function only of z, L = − α,γ n αγ log(x + C eq z) This is exactly the same form of the likelihood function in §A.3 after the single equality constraint there is eliminated. Hence, to obtain the Cramér-Rao bound (80) repeat, mutatis mutandis, the procedure in §A.3.