Verification of state and entanglement with incomplete tomography

There exists, in general, a convex set of quantum state estimators that maximize the likelihood for informationally incomplete data. We propose an estimation scheme, catered to measurement data of this kind, to search for the exact maximum-likelihood-maximum-entropy estimator using semidefinite programming and a standard multi-dimensional function optimization routine. This scheme can be used to infer the expectation values of a set of entanglement witnesses that can be used to verify the entanglement of the unknown quantum state for composite systems. Next, we establish an alternative numerical scheme that is more computationally robust for the sole purpose of maximizing the likelihood and entropy.


Introduction
Quantum state preparation is the first important step for any protocol that makes use of quantum resources. Examples of such protocols are quantum state teleportation and quantum key distribution which require entangled quantum states. In order to verify the integrity of the true quantum state ρ true prepared by the source, one carries out quantum state tomography to characterize it. Measurements are performed on a collection of identical copies of quantum systems (electrons, photons, etc.) that are emitted from the source. Then, the quantum state of the source is inferred from the measurement data obtained from this collection. The measurements are generically described by a set of positive operators {Π j } that compose a probability operator measurement (POM). Such a procedure of state inference is known as quantum state estimation.
When the measurement outcomes form an informationally complete set, they fully characterize the source and the measurement data obtained will contain maximal information about its state. To infer the unknown state from the data, one can search for a state estimator that maximizes the likelihood functional that yields the probability of obtaining a particular sequence of measurement detections given a quantum state -the maximum-likelihood (ML) estimator [1,2,3,4]. Yet, in tomography experiments performed on complex quantum systems with many degrees of freedom, it is not possible to implement such an informationally complete set of measurement outcomes. Therefore, some information about the source will be missing and its quantum state cannot be unambiguously determined. For instance, if a source produces a mode of light that is described by an infinite-dimensional statistical operator ρ true , then no matter how ingeniously a measurement scheme is designed to probe incoming photons prepared by this source, an infinite amount of information about the mode of light will always remain unknown. The ML estimator obtained from these informationally incomplete data is no longer unique and there will in general be infinitely many other ML estimators that are consistent with the data. These estimators form a convex set under the likelihood plateau.
In order to choose an estimator from the convex set for statistical prediction, we can consider the maximum-entropy principle advocated by E. T. Jaynes [5,6]. In doing so, one obtains a unique estimator that maximizes both the likelihood and the von Neumann entropy functional. Statistically, this estimator is least-biased for the informationally incomplete data. In [7,8,9], we developed and applied an algorithm, based on the steepest-ascent method, to approximately look for the maximum-likelihood-maximum-entropy (MLME) estimator. This algorithm involves a parameter that needs to be chosen just above a minimum threshold to obtain an estimator that is as close to the actual MLME estimator as possible. In general, this threshold depends on the true state ρ true . Therefore, one needs to run the algorithm a few times to estimate this threshold. This paper is organized as follows. After a brief review of this steepest-ascent algorithm in Section 2, we introduce a numerical scheme that is based on completely different principles to directly search for the MLME estimator within the convex set in Section 3. This scheme couples two separate optimization techniquessemidefinite programming (SDP) and a derivative-free optimization method -to maximize the entropy over linear combinations of a maximal set of linearly independent ML estimators that spans the ML convex set. It will be shown that, owing to the mechanisms of SDP, one can make use of this scheme to infer the expectation values of a set of entanglement witnesses to verify the presence of entanglement in the unknown quantum state for composite systems. Finally, in Section 4, we establish a more robust numerical scheme that systematically generates the maximal set of linearly independent ML estimators that defines the convex set without fail. Instead of SDP, this scheme utilizes a nonlinear optimization routine that finds the global maximum of a highly nonlinear functional that is used to generate this maximal set.

Brief review
The likelihood functional log L({n j }; ρ) for a set of measurement data {n j } collected with a POM j Π j = 1 is given by where n j refers to the number of occurrences of the outcome Π j and p j = tr {ρΠ j }.
The corresponding frequencies are given by f j = n j /N . In ML state estimation, the concave log-likelihood functional log L({n j }; ρ) is maximized to obtain the ML estimatorρ ML ‡ for the given set of data. If the number of linearly independent outcomes is D 2 , with D being the dimension of the Hilbert space, the data is informationally complete and the estimatorρ ML is unique. In the case where this number is less than D 2 , the data is informationally incomplete and there exists now a continuous set ofρ ML that yield the same maximal likelihood. This set is convex since the likelihood functional is concave in ρ.
To choose the estimatorρ MLME that has the highest entropy out of this convex set, we can consider the Lagrange functional ‡ The symbol "ˆ" is used to denote all estimators.
where λ is the Lagrange multiplier corresponding to the constraint of maximal entropy S(ρ MLME ) = − tr{ρ MLME logρ MLME } = S max . In doing so, we maximize the two functionals S(ρ) and log(L({n j }; ρ)) simultaneously. We defineρ I,λ to be the estimator that maximizes I(λ; ρ). If λ = 0, I(λ; ρ) = log L({n j }; ρ)/N and maximizing this Lagrange functional is just the procedure of ML. Since the data is informationally incomplete, there exists a convex plateau structure for the log-likelihood functional and maximizing I(λ; ρ) yields a convex set of estimators. For large λ values, the term λS(ρ) dominates, so that the resulting estimatorρ I,λ→∞ = 1/D. When λ takes on a very small positive value [7], the contribution from λS(ρ) becomes relatively much smaller than that of log(L({n j }; ρ))/N , and any variation of the von Neumann entropy functional is only detectable over the state space region that coincides with the likelihood plateau. In other words, maximizing I(0 < λ → 0; ρ) is equivalent to maximizing the entropy over the plateau. Therefore,ρ I,0<λ→0 =ρ MLME . The iterative algorithm that is based on the steepest-ascent method is described in [7]- [9]. In practice, there is a limit to how small λ can be. In particular, when λ is smaller than some numerical threshold λ thres > 0, the gradient of λS(ρ) is no longer visible. In this case, the algorithm treats I(λ < λ thres ; ρ) as I(λ = 0; ρ) and performs ML estimation. Hence, the optimal parameter λ is to be slightly above λ thres so thatρ I,λ λ thres ≈ρ MLME . This inevitably introduces a small bias to the estimator.
Determining the value of λ thres analytically is quite complicated because λ thres is actually a function of the true state and the POM, that is λ thres = λ thres ({Π j }; ρ true ).
Since ρ true is unknown to us, one needs to estimate λ thres through repeated runs of the algorithm. In the next section, we will introduce an alternative scheme to search for the exact MLME estimator.

Algorithm for entanglement and state verification
To obtain the unique MLME estimator that has the highest entropy, a search has to be performed within the ML convex set. In this section, we introduce a numerical procedure to estimateρ MLME in two steps. The first step is to obtain a collection of boundary states of the convex set for the measurement data obtained. In the next step, the operatorρ MLME can be estimated using these boundary states with a standard function optimization routine.
To carry out the first step, we need to identify the boundary of the ML convex set. Especially for large dimensions, the boundary of the convex set has an extremely complicated geometry that is too difficult to be analytically determined. Instead, we investigate its boundary by numerical means. We begin with the fact that the real functional f (H; ρ) = tr{ρH}, where the operator H = H † , is a linear functional of ρ. Therefore, if we try to maximize (or minimize) f (H; ρ) over some subspace of ρ that has a well defined boundary, the maxima (or minima) of this linear functional are always on the boundary of this subspace. We can thus generate boundary states by maximizing or minimizing f (H; ρ), for a given Hermitian operator H, over the ML convex set. This problem is equivalent to the following optimization task: Maximize or minimize f (H; ρ) = tr{ρH}, subjected to the following constraints: This is a standard linear optimization problem, with linear and positivity constraints, that can be solved with the help of SDP [10].
At this point, we need to find out the minimum number of boundary states that is required to search forρ MLME . To do this, we represent a generic ML estimatorρ ML by a set of D 2 linearly independent § trace-orthonormal Hermitian basis operators With these basis operators, we can separateρ into the part in the measurement subspace of dimension D meas , and the rest of the state space that constitutes the unmeasured parameters (the unmeasured subspace) The generation of the basis operators Γ meas j that span the measurement subspace can be done by the Gram-Schmidt orthonormalization procedure on the Π j s, with all conventional inner products for vectors replaced by trace inner products for operators. In other words, the number of linearly independent POM outcomes Π j is D meas . To generate the rest of the basis operators that span the unmeasured subspace, we continue the Gram-Schmidt procedure using randomly generated positive operators instead of the Π j s.
From (3), it can be deduced that the maximum dimension of the ML convex set is D unmeas . To show this, we note that every ML estimator contains the same operator ρ meas since the probabilities of p j = tr{ρ meas Π j } are fixed and tr{ρ unmeas Π j } = 0 for all Π j by definition of the trace-orthonormal basis operators. This implies that § If L Hermitian operators A j are all linearly independent, the rank of the matrix M with elements the only difference between any two ML estimators in the convex set isρ unmeas . As the operatorsρ meas andρ unmeas are linearly independent, it follows that any ML estimator can always be expressed as a linear combination of the uniqueρ meas and the D unmeas linearly independent basis operators that defineρ unmeas . This means that a set of D unmeas + 1 linearly independent boundary states is enough to look for the MLME estimator. As the operatorρ meas is fixed by the measurement operators, the maximal number of free parameters that span the convex set is D unmeas . In some cases, however, the dimension of the ML convex set is lower due to the positivity constraint imposed on the ML estimators. In the extreme case, the convex set is restricted to a single point in state space. In these situations, we do not know its actual dimension and repeated generation of boundary states is necessary to estimate the maximum number of linearly independent boundary states. We remind the reader that with enough linearly independent states, the exact MLME estimator can be obtained up to numerical precision. In Section 4, a different numerical procedure will be introduced to generate the maximal set of linearly independent ML estimators that spans the ML convex set without requiring any knowledge of the convex set.
Apart from serving as a routine to numerically compute the boundary states of the convex set, SDP provides an additional useful function. For composite quantum systems, if one selects the Hermitian operators H to be entanglement witnesses W , one can obtain information about the presence of entanglement in the unknown true state even with informationally incomplete data. These witnesses have the properties that tr{ρ sep W } ≥ 0 for all separable states ρ sep and tr{ρ ent W } < 0 for at least one entangled state ρ ent . The SDP routine, described above, thus looks for the maximum (or minimum) value of f (W ; ρ) for any chosen witness operator W over the space of positive ρs. This way, we can in fact infer a set of maximum (or minimum) witness expectation values over the ML convex set from this optimization procedure. The set of inferred maximum witness expectation values is particularly informative, for if the maximum value of f (W ; ρ) for at least one of the randomly generated operators W is negative, we can immediately conclude that the true state is entangled since this state must lie within the ML convex set that results from the incomplete data. For practical computation, we can choose the witness operators W to be of the decomposable form W = Q tj , where Q is a positive operator with no product kets in its range and the symbol " t j " denotes a partial transposition with respect to the jth subsystem. For bipartite systems, these operators are optimal witnesses [11,12], that is no other witnesses can detect all entangled states that are detected by this witness, as well as other entangled states. For multipartite systems, these operators still serve as entanglement witnesses since tr ρ sep Q t 2 = tr ρ t2 sep Q > 0, although they are no longer optimal in general. To obtain a random set of boundary states of the ML convex set, random statistical operators Q are generated using the relation where X is a random operator which, when expressed in the computational basis, has complex matrix elements that are distributed according to the standard normal distribution of zero mean and unit variance. The second step involves an optimization procedure to maximize the entropy S(ρ) using the generated set of M 0 ≤ D unmeas + 1 linearly independent boundary states . For the purpose of entanglement detection, these boundary states are obtained by maximizing the linear functionals f (W ; ρ) of a set of witness operators W over the ML convex set according to the recipe described above. We start by writing a generic ML estimator as a linear combination of ρ where the t j s are normalized coefficients, such that j t j = 1, that are in general real such thatρ ML ({t j }) ≥ 0. The task now is to look for the values of t j = t max j for which } is maximum over all real normalized coefficients. The unconstrained optimization of this M 0 -dimensional function with respect to t j can be performed with any efficient multi-dimensional unconstrained optimization routine that is included in the standard libraries of commercialized mathematical softwares. In MATLAB, for instance, the function fminsearch does the job using the Nelder-Mead simplex method (NMS). When M 0 is large, it is suggested in [13] that an adaptive version of the Nelder-Mead simplex method (ANMS) may be more advantageous in terms of shorter computation time. We take the resulting operator ρ SDP ≡ρ ML {t max j } as the SDP MLME estimator. There is, however, a caveat to this optimization. Since the positivity ofρ ML ({t j } ) is no longer guaranteed over the entire space of real normalized vectors, the entropy expressed in terms of the eigenvalues ν j of the statistical operatorρ ML ({t j }) = j |ν j ν j ν j | that is represented by its eigenbasis {|ν j }, can take complex values. In order to restrict the optimization to yield only positive SDP MLME estimators, one can replace the entropy function in (6) with the conditional function which effectively restricts the original search region to the admissible state space.
To check if the evaluated operatorρ ML ({t j }) is positive, a highly efficient way is to determine whether or not it admits a Cholesky decomposition.
For full-rank SDP MLME estimatorsρ SDP , there is a simple way to check that ρ SDP is indeed the MLME estimator. We recall that the form of the MLME estimator is given byρ where each λ j is a Lagrange multiplier for the constraint tr{ρ MLME Π j } = tr{ρ ML Π j } for any ML estimatorρ ML in the convex set. Ifρ MLME is full-rank, it follows from (8) that the operator logρ MLME is a linear combination of only the POM outcomes, that is logρ MLME resides in the measurement subspace. This is equivalent to the set of conditions tr{Γ unmeas j logρ SDP } = 0 for all Γ unmeas j s. Defining the variables c j = tr{Γ unmeas j logρ SDP }, the quantity γ ≡ j c 2 j can be used to determine if ρ SDP is close to the actual MLME estimator up to some numerical precision.
In summary, both the MLME estimator and information about the entanglement of ρ true can be obtained with informationally incomplete data using the following algorithm, which we coin as SDP MLME:

SDP MLME
(i) ML -Obtain the ML probabilities for the POM used to collect the measurement data with the ML algorithm.
(ii) SDP -Perform SDP, using (4), to obtain the maximal set of M 0 ≤ D unmeas + 1 linearly independent boundary states that define the ML convex set using witness operators. Inspect the list of maximum expectation values of the witness operators and conclude that the unknown quantum state is entangled if any one of them is negative.
(iii) Entropy maximization -Carry out function maximization with NMS or ANMS on the entropy function S ({t j }) with respect to all real normalized coefficients t j using (5) and (7) to obtainρ SDP . Figure 1 summarizes the results. In Figure (1a), the convergence of γ is consistent with the fact that, in principle, the maximal number of nine linearly independent boundary states is enough to express the MLME estimator. The entanglement detection ratio in Figure 1b is computed as the ratio of the number of detected pure states to the total number of generated states for every number of boundary states used. Ideally, this  Figure 1: Three plots showing various aspects of SDP MLME conducted on 100 random twoqubit entangled pure states as functions of the number of linearly independent ML boundary states used, with plots (a) and (c) averaged over all these pure states. A randomly-generated eight-outcome informationally incomplete POM is used throughout the simulations. In this case, every corresponding ML convex set is specified by nine linearly independent ML estimators. Decomposable optimal witness operators of the form Q t 2 , where the operators Q are entangled pure states, are used to generate the boundary states with SDP.
ratio eventually goes to one if enough witness operators are generated to detect all the random pure states since there are no positive partial-transpose entangled states in this case [11,12]. For benchmarking, Figure 1c is generated to confirm the consistency of SDP MLME with the MLME algorithm described in Section 2. There exists an average bias in 1c that arises from a fixed λ thres = 10 −5 for all the pure states.

Robust algorithm for incomplete state estimation
Despite the usefulness of SDP in entanglement verification and boundary states generation as discussed in Section 3, the speed of the SDP routine strongly depends on the dimension of the Hilbert space and the total number of linear constraints imposed by the measurement data. When the total number of POM outcomes is large, there will generally be a considerable slowdown of the SDP routine as the search accounts for a large set of linear constraints in addition to the positivity constraint. Another feature of the SDP routine is that the sequence of boundary states that are generated from random Hermitian operators are not guaranteed to be linearly independent of one another. This means that typically, one would need to generate a large set of ML estimators that contains the maximal number of linearly independent estimators that span the ML convex set. For convex sets of large dimensions, this approach can be time-consuming. In this section, we propose a different search routine, in place of SDP, to directly look for linearly independent ML estimators within the ML convex set in a deterministic way. With this routine, we establish a feasible algorithm to look for the exact MLME estimator.
To begin, we recall that a given set of M 0 linearly independent ML estimators, as in (5) This hints a straightforward strategy to cumulatively obtain the maximal set of linearly independent ML estimators that spans the entire ML convex set: Starting with a single ML estimator, the next estimator containing the sameρ meas should be chosen such that the smallest eigenvalue σ min (M g ) of the Gram matrix M g for these two estimators is maximized, and so forth, with the maximization performed over positive estimatorsρ (j) ML ≥ 0. In general, σ min (M g ) is a nonlinear functional of M g that has multiple local stationary points. This functional is also not differentiable and has undefined gradients at the boundary of the state space. To search for its global maximum, an appropriate numerical method to use is a nonlinear optimization algorithm that invokes pattern searches [14] and can cope with functionals that have ill-defined gradients. The solver for this algorithm, patternsearch, is readily available in MATLAB. The positivity constraintρ (j) ML ≥ 0 is incorporated into the optimization algorithm using the augmented Lagrangian method with Cholesky decomposition. Another versatile feature of the proposed routine is that it can be set to terminate when the maximal number of linearly independent ML estimators is generated, such that the next ML estimator always yields a zero eigenvalue for M g . We can, therefore, deterministically obtain the maximal set of linearly independent ML estimators that spans the ML convex set in this manner, in contrast with the SDP routine in SDP MLME. In this sense, the routine is operationally robust even without any prior information about the convex set. To ensure that the search is numerically stable, it is favorable to start the pattern search algorithm from a highly-mixed ML estimator. This is obtained by performing SDP as a couple of times and defining the starting estimator as an equal mixture of the resulting ML boundary estimators and the fairly mixed ML estimator obtained from the ML algorithm starting from the maximally-mixed state in computing the ML probabilities. Using this maximal set, the MLME estimator can be directly obtained by maximizing the entropy function in (7) over all linear combinations of the linearly independent ML estimators in the set. These lead to the following pattern search MLME algorithm (PS MLME) for incomplete quantum state estimation:

PS MLME
(i) ML -Obtain the ML probabilities, as well as the corresponding ML estimator that is fairly mixed, for the POM used to collect the measurement data with the ML algorithm starting from the maximally-mixed state.
(ii) Carry out SDP a couple of times to obtain a few ML boundary estimators and check if the ML estimators are close to each other. If the average pairwise norm is smaller than a certain threshold, this means that the ML convex set can be approximated to be a single point and the ML estimator is taken to be the unique estimator. Otherwise, proceed to the next step.
(iii) Definition of the ML convex set -Generate the maximal set of linearly independent ML estimators that defines the ML convex set by maximizing σ min (M g ) via the augmented Lagrangian pattern search algorithm.
(iv) Entropy maximization -Carry out function maximization with NMS or ANMS on the entropy function S({t j }) with respect to all real normalized coefficients t j using (5) and (7) to obtain the MLME estimator.
The performance of PS MLME depends not only on the dimension D unmeas of the unmeasured subspace, but also on the complexity of the functional σ min (M g ). More generally, as the dimension of the Hilbert space, or that of the unmeasured subspace, increases, the number of local maxima of σ min (M g ) increases. This translates to a longer computation time to locate the global maximum in the search for linearly independent ML estimators. Thus, for very large dimensions, PS MLME becomes inefficient and the approximate MLME algorithm, which is based on steepest ascent [7], is a more practical substitute. On the other hand, PS MLME consistently gives more accurate MLME estimators as compared to the steepest-ascent algorithm, which yields biased results for large dimensions because of the finite λ parameter. Hence, there is a tradeoff between computation time and the accuracy of the MLME estimators. Figure 2 compares the performances of the PS MLME and the steepest-ascent MLME algorithms for varying Hilbert space dimensions. Figure 2b shows the consistently more accurate results obtained with PS MLME as compared to the steepest-ascent algorithm, in exchange for its longer computation time, illustrated in Figure 2a, for higher dimensions. The simulations show that for the moderately large dimensions considered in Figure 2, much more accurate MLME estimators can be obtained using PS MLME with longer computation time that is of the same order as that with steepest ascent. To reduce the computation time of PS MLME, one can consider an approximate maximization of σ min (M g ) as long as the resulting value is sufficiently large. Throughout the simulations, the duration of searching for each optimal ML estimator is restricted to five seconds. For even larger dimensions, PS MLME can be computationally demanding even when approximate maximization is carried out, and the steepest-ascent algorithm turns out to be a more realistic option. Efforts to further improve the performance of the PS MLME scheme are in the works. Nevertheless, we hope that the current work can serve as a stepping stone that helps to spur interesting discussions and novel contributions in related fields on numerical optimization over convex sets of positive operators.

Conclusion
We have introduced a scheme to look for the unique estimator that maximizes the likelihood and entropy for informationally incomplete measurement data. This involves two main procedures: a generation of linearly independent boundary maximum-likelihood estimators that spans the convex set with semidefinite programming, and an entropy maximization with these estimators using a standard function optimization routine. Furthermore, for composite quantum systems, one can apply this scheme to infer the expectation values of a set of entanglement witnesses.
This information allows us to verify the entanglement of any quantum state with informationally incomplete data. However, semidefinite programming does not offer a definite control over the generation of maximum-likelihood estimators. This motivated us to develop an alternative scheme that is more operationally robust than the former one to search for the maximum-likelihood-maximum-entropy estimator. This latter scheme makes use of the pattern search optimization algorithm that is suitable for maximizing a nonlinear function required to deterministically generate the maximal set of estimators that defines the convex set. With numerical simulations, we showed that the latter scheme gives much more reliable results than the MLME algorithm discussed in Section 2 at the expense of slightly longer computation time when the dimension of the reconstruction Hilbert space, or the unmeasured subspace, is moderately large.