Error regions in quantum state tomography: computational complexity caused by geometry of quantum states

The outcomes of quantum mechanical experiments are inherently random. It is therefore necessary to develop stringent methods for quantifying the degree of statistical uncertainty about the results of quantum experiments. For the particularly relevant task of quantum state estimation, it has been shown that a significant reduction in uncertainty can be achieved by taking the positivity of quantum states into account. However -- the large number of partial results and heuristics notwithstanding -- no efficient general algorithm is known that produces an optimal uncertainty region from experimental data and the prior constraint of positivity. Here, we make this problem precise and show that the general case is NP-hard. Our result leaves room for the existence of efficient approximate solutions, and therefore does not yet imply that the practical task of quantum uncertainty quantification is intractable. However, it does show that there exists a non-trivial trade-off between optimality and computational efficiency for error regions. We prove two versions of the result: One for frequentist and one for Bayesian statistics.


Introduction
The outcomes of quantum mechanical measurements are subject to intrinsic randomness. As a result, all information we obtain about quantum mechanical systems are subject to statistical uncertainty. It is thus necessary to develop stringent methods for quantifying the degree of uncertainty. These allow one to decide whether an observed feature can be trusted to be real, or whether it may have arisen from mere statistical fluctuations (a "fluke").
In this paper, we concentrate on uncertainty quantification for quantum state tomography (QST) ‡ Here, the task is to infer a density matrix ̺ 0 , associated with a preparation procedure of a finite-dimensional quantum system, from the outcomes of measurements on independent copies of the system. In addition to an estimate̺ for the unknown true state ̺ 0 , a tomography procedure should rigorously quantify the remaining statistical uncertainty.
We note that QST is an established experimental tool -in particular in quantum information-inspired setups. It has been used to characterize quantum states in a large number of different platforms -Refs. [4,5,6,7,8,9,10,11,12] are an incomplete list.
From a technical point of view, uncertainty quantification in QST may seem to be straight-forward. After choosing a measurement to perform (i.e. by specifying a POVM), the probability distribution over the outcomes is a linear function of the unknown state ̺ 0 . Inference and uncertainty quantification in linear models are wellstudied problems of mathematical statistics. What makes the QST problem special is the additional constraint that the density matrix ̺ 0 be positive semi-definite (psd ) and have unit-trace. This shape constraint can lead to a significant reduction in uncertainty -in particular if the true state ̺ 0 is close to the boundary of state space: In this case, it is plausible that a large fraction of possible estimates that seem compatible with the observations can be discarded, as they lie outside of state space.
Indeed, it is known that taking the psd constraint into account can result in a dramatic -even unbounded -reduction in uncertainty. Prime examples are results that employ positivity to show that even informationally incomplete measurements can be used to identify a quantum state with arbitrarily small error [13,14,15,16,17,18]. More precisely, these papers describe ways to rigorously bound the size of a confidence region for the quantum state based only on the observed data and on the knowledge that the data comes from measurements on a valid quantum state. While these uncertainty bounds can always be trusted without further assumptions, only in very particular situations have they been proven to actually become small. These situations include the cases where the true state ̺ 0 is of low rank [16,17], or admits an economical description as a matrix-product state [13].
It stands to reason that there are further cases -not yet identified -for which the size of an error region can be substantially reduced simply by taking into account the quantum shape constraints. This motivates the research program this paper is part of: Understand the general impact of positivity constraints on uncertainty quantification in QST.
The positive results cited above notwithstanding, it is not obvious how to take the a priori information of positive semi-definiteness into account algorithmically. The fact that no practical and optimal general-purpose algorithm for quantum uncertainty quantification has been identified could either reflect a limit in our current understanding -or it could indicate that no efficient algorithm for this problem exists.
In this work, we present first evidence that optimal quantum uncertainty quantification is algorithmically difficult. We give rigorous notions of optimality both from the point of view of Bayesian statistics (where this concept is fairly canonic) and of "orthodox" statistics (where some choices have to be made). We exhibit special cases for which there does exist an efficient algorithm that identifies optimal error regions. However, our main result proves that in general, finding these regions is NP-hard and thus computationally intractable. By working under assumptions that render the unconstrained problem tractable, we show that computational intractability arises solely due to the quantum constraints and not due to the general difficulties of high-dimensional statistics.
The present results do not by themselves imply that the practical problem of uncertainty quantification is unfeasible. For applications, "almost-optimal" regions would be completely satisfactory. And indeed, a number of techniques for tackling this problem in theory and practice have been proposed (e.g. based on sample splitting, resampling, or on approximations for Bayesian posterior distributions -c.f. Sec. 1.4). Each of these methods is known analytically or from numerical experiments to perform well in some regimes. However, this paper does establish that there is a non-trivial trade-off between optimality and computational efficiency in quantum uncertainty quantification. What is more, our work might help guide future efforts that aim to design efficient and optimal estimators: With a very natural construction proven not to possess an efficient algorithm in general, it is now clear that researchers must focus on approximations that circumvent our hardness results. In general, we hope that this work establishes a framework for future positive and negative results, which will eventually allow us to understand which performance can be achieved.
The rest of this paper is structured as follows. In the subsections below, we comment on use-cases of full QST for high-dimensional quantum systems, summarize related works, and clarify the (non-trivial) issue of "optimality" in uncertainty quantification. We then establish the main result for orthodox statistics in Section 2 and follow up with a Bayesian treatment in Section 3.

The need for full tomography
A large number of tomography experiments for quantum systems with hundreds of dimensions has been published, e.g. [19]. However, it is not completely obvious that this approach will continue to make sense as dimensions scale up further.
Indeed, a variety of theoretical tools for quantum hypothesis testing, certification, and scalar quantum parameter estimation [20,21,22,23,24,25] have been developed in the past years, that avoid the costly step of full QST. Examples include entanglement witnesses [22] and direct fidelity estimation [23].
However, there remain use cases that necessitate full-fledged QST. We see a particularly important role in the emergent field of quantum technologies: Any technology requires means of certifying that components function as intended and, should they fail to do so, identify the way in which they deviate from the specification.
As an example, consider the implementation of a quantum gate that is designed to act as a component of a universal quantum computing setup. One could use a certification procedure -direct fidelity estimation, say -to verify that the implementation is sufficiently close to the theoretical target that it meets the stringent demands of the quantum error correction threshold. If it does, the need for QST has been averted. However, should it fail this test, the certification methods give no indication in which way it deviated from the intended behavior. They yield no actionable information that could be used to adjust the preparation procedure. The pertinent question "what went wrong" cannot be cast as a hypothesis test.
Thus, while many estimation and certification schemes can -and should -be formulated without resorting to full tomography, the above example shows that QST remains an important primitive.

Error Regions
As inference based on empirical data is one of the main topics of statistics, it is natural to apply the established notions of uncertainty quantification to QST. These are either confidence regions in orthodox statistics [26] or credible regions in Bayesian statistics [27]. The two approaches give rise to different techniques, but most importantly, have very distinct interpretations [28].
In orthodox (or frequentist) statistics, the task of parameter estimation can be summarized as follows: We assume that the observed data is generated from a parametric model with true parameter Θ, which is unknown. From a finite number of observations X 1 , . . . , X N , we must construct an estimateΘ that should be "close to" the true value Θ in some sense. The function that maps data to such an estimate is called a (point) estimator. A confidence region C with coverage α is a region estimator -that is a function that maps observed data to a subset of the parameter space -such that the true parameter is contained within it with probability greater than α Note that the defining property of a confidence region concerns the behavior of the random function C over the course of many (hypothetical) repetitions of the experiment. No statement is made about a single run. Of course, Eq. (1) does not uniquely determine a confidence region; it does not even guarantee a sensible quantification of uncertainty, as C equal to the whole parameter space fulfills this condition trivially. Therefore, we consider confidence regions that perform well with respect to (w.r.t.) some notion of optimality: In general, smaller regions should be preferred since they convey more confidence in the estimate and exclude more alternatives. But since the size -as measured by volume -of a confidence region may depend on the particular data sample as well as the true value of the parameter, different notions of optimality have been introduced [29].
Bayesian statistics on the other hand treats the parameter Θ itself as a random variable. The distribution over Θ reflects our knowledge about the parameters [27]. Ahead of observing any data, one has to choose a prior distribution, which represents our a priori beliefs. The observed data is then incorporated using Bayes' rule to update the distribution yielding the posterior P(Θ|X 1 , . . . , X N ). A credible region C (we denote both confidence and credibility regions by the same letter) with credibility α is defined as a subset of the parameter space containing at least mass α of the posterior P(Θ ∈ C|X 1 , . . . , X N ) ≥ α.
In contrast to the orthodox setting, here, the data is assumed to be fixed and the probability is assigned w.r.t. Θ.
Since the posterior distribution is uniquely defined by the choice of prior and the data, there is less ambiguity in the choice of a notion of optimality: The most natural choice are minimal-volume credible regions. In case the posterior has the probability density π(θ) w.r.t. the volume measure, these are given by regions of highest posterior density where λ is determined by the saturation of the credibility level condition (2).

Positivity of quantum states
When attempting to construct optimal error regions for QST, we should exploit the physical constraints at hand in order to reduce their size and, therefore, make them more powerful: every valid density matrix ̺ -apart from being Hermitian and normalized -must be positive semidefinite (psd). More formally, in a d-dimensional scenario it is required that Here, S + denotes the set of valid mixed quantum states, which is a proper subset of the real vector space S of Hermitian matrices with unit trace. While the first two properties (hermiticity and normalization) are linear constraints and therefore easy to take into account by virtue of an appropriate parametrization, positivity is far more challenging to employ constructively. A prime example where this structural information is crucial in the construction of optimal error regions is the application of compressed sensing techniques to QST [14,16,30]. Compressed sensing allows to recover a low-rank state from informationally incomplete measurements. Without further assumptions, this can lead to unbounded error regions -c.f. the discussion of Pauli designs in [30] and Sec. 2.1. Nevertheless, the constraints implied by physical states allow for the construction of confidence regions in this setting [30], that are of finite size and that become arbitrarily small as the individual measurement errors tend to zero.
However, as the cited work is specifically tailored to the compressed sensing scenario, it is not clear how to extend it to the general setting of QST. The purpose of this work is to explore the degree to which positivity can be taken into account in general, if one assumes that computational power is bounded.

State of the art
In practice (e.g. [19]), uncertainty quantification for tomography experiments is usually based on general-purpose resampling techniques such as "bootstrapping" [31]. A common procedure is this: For every fixed measurement setting, several repeated experiments are performed. This gives rise to an empirical distribution of outcomes for this particular setting. One then creates a number of simulated data sets by sampling randomly from a multinomial distribution with parameters given by the empirical values. Each simulated data set is mapped to a quantum state using maximum likelihood estimation. The variation between these reconstructions is then reported as the uncertainty region. There is no indication that this procedure grossly misrepresents the actual statistical fluctuations. However, it seems fair to say that its behavior is not well-understood. Indeed, it is simple to come up with pathological cases in which the method would be hopelessly optimistic: E.g. one could estimate the quantum state by performing only one repetition each, but for a large number of randomly chosen settings. The above method would then spuriously find a variance of zero.
On the theoretical side, some techniques to compute rigorously defined error bars for quantum tomographic experiments have been proposed in recent years. The works of Blume-Kohout [32] as well as Christandl, Renner, and Faist [33,34] exhibit methods for constructing confidence regions for QST based on likelihood level sets. While very general, neither paper provides a method that has both a runtime guarantee and also adheres to some notion of non-asymptotic optimality [26,35].
Some authors have proposed a "sample-splitting" approach, where the first part of the data is used to construct an estimate of the true state, whereas the second part serves to construct an error region around it [16] (based on [23]), as well as [30]. These approaches are efficient, but rely on specific measurement ensembles (operator bases with low operator norm), approach optimality only up to poly-logarithmic factors, and -in the case of [16,23] -rely on adaptive measurements.
Regarding Bayesian methods, the Kalman filtering techniques of [21] provide a efficient algorithm for computing credible regions. This is achieved by approximating all Bayesian distributions over density matrices by Gaussians and restricting attention to ellipsoidal credible regions. The authors develop a heuristic method for taking positivity constraints into account -but the degree to which the resulting construction deviates from being optimal remains unknown. A series of recent papers aim to improve this construction by employing the particle filter method for Bayesian estimation and uncertainty quantification [36,37,38]. Here, Bayesian distributions are approximated as superpositions of delta distributions and credible regions constructed using Monte Carlo sampling. These methods lead to fast algorithms and are more flexible than Kalman filters with regard to modelling prior distributions that may not be well-approximated by any Gaussian. However, once more, there seems to be no rigorous estimate for how far the estimated credible regions deviate from optimality. Finally, the work in [39] constructs optimal credible regions w.r.t. a different notion of optimality: Instead of penalizing sets with larger volume, they aim to minimize the prior probability as suggested by [40].

Orthodox Confidence Regions
In this section we are going to present the first major result of this work concerned with orthodox confidence regions in QST. Optimal confidence regions for such highdimensional parameter estimation problems are quite intricate even without any constraints on the allowed parameters. There are only few elementary settings, where optimal error regions are known and easily characterized.
Since the goal of this work is to demonstrate that quantum shape constraints severely complicate even "classically" simple confidence regions, in the further discussion we restrict the discussion to a simplified setting: We focus on confidence ellipsoids for Gaussian distributions, which are one of the few easily characterizable examples. Furthermore, by local asymptotic normality, these arise as a natural approximation in the limit of many measurements. As we show in the following, even characterizing these highly simplifying ellipsoids with the quantum constraints taken into account constitutes a hard computational problem. On the other hand, as indicated in the introduction, these structural assumptions may help to reduce the uncertainty tremendously. Therefore, our work can be interpreted a trade-off between computational efficiency and statistical optimality in QST.

Optimal confidence regions for quantum states
As already indicated in Sec. 1.3, the additional information that the true quantum state ̺ 0 must belong to the set of positive semidefinite matrices S + ⊂ S can be exploited to possibly improve any confidence region for QST. This is especially clear for notions of optimality with a loss function stated in terms of volume § Vol(·), as we will show in this section.
We consider an especially simple procedure to take the positivity constraints into account, namely truncating all non-positive matrices from tractable confidence regions for the unconstrained problem. This approach is mainly motivated by the goal to show that the computational intractability exclusively stems from the quantum constraints and is not caused by difficulties of high-dimensional statistics in general. Furthermore, we prove in Lemma 1 that some notions of optimality, e.g. admissibility, are preserved under truncation. In other words, there are notions of optimality such that truncation of an optimal confidence region for the unconstrained problem gives rise to an optimal region for the constrained one.
First, let us introduce the notion of admissibility as given by [41,Def. 2.2].

Definition 1.
A confidence region C for the parameter estimation of ̺ 0 ∈ S is called (weakly) admissible if there is no other confidence region C ′ that fulfills (i) (equal or smaller volume) Vol(C ′ (y)) ≤ Vol(C(y)) for almost all observations y ∈ R m (ii) (same or better coverage) P(C ′ ∋ ̺ 0 ) ≥ P(C ∋ ̺ 0 ) for all ̺ 0 ∈ S (iii) (strictly better) strict inequality holds for one ̺ 0 ∈ S in (ii) or on a set of positive measure in (i).
In words, C is admissible if there is no other confidence region C ′ that performs at least as good as C and strictly better for some settings. The conditions in Def. 1 are stated only for "almost all" y, since one can always modify the region estimators on sets of measure zero without changing their statistical performance. A different approach is to state condition (i) in terms of the expected volume , which leads to the notion of strong admissibility [41, Def. 7.1].
Def. 1 can also be stated for the parameter estimation with physical constraints, i.e. when ̺ 0 ∈ S + . The question is: Can we obtain admissible confidence regions C + ⊂ S + for the constrained setting from admissible confidence regions C ⊂ S of the unconstrained estimation problem? The following Lemma answers this question with a simple geometric construction: Lemma 1. Let C denote an admissible confidence region for the unconstrained estimation problem for the parameter ̺ 0 ∈ S. Then, C ∩ := C ∩ S + is an admissible confidence region for the constrained problem with ̺ 0 ∈ S + .
Proof. Under the assumption that C ∩ is not admissible, there must exist a "better" confidence region C + for the constrained parameter estimation problem. W.l.o.g. assume that both C + and C ∩ have the same coverage. Therefore, we must have Vol(C + (y)) ≤ Vol(C ∩ (y)) for almost all observations y ∈ R m , and there is a set Y ⊂ R m of non-zero measure such that Vol(C + (y)) < Vol(C ∩ (y)) for y ∈ Y . Define a new confidence region for the unconstrained problem where C c = C \ C ∩ denotes the compliment of C ∩ in C. Then, C ′ has the given coverage level, since C + provides coverage for ̺ 0 ∈ S + , whereas C c provides coverage for the case ̺ 0 ∈ S \ S + . Furthermore, we have for almost all y Finally, strict inequality holds in Eq. (6) for all y ∈ Y due to the assumption on C + . However, this would imply C not being admissible in contradiction to the assumptions of the Lemma.
One criticism raised against the use of the truncated confidence regions is the possibility that they may yield empty realizations and, hence, are considered "unphysical" [42]. However, according to the standard definition in Sec. 1.2, a procedure that reports 95% confidence regions is allowed to give any result 5% of the time.
Furthermore, a different strategy often adopted for point estimator is to use an unconstrained parametrization for the constrained parameter space. A typical example is a coin toss model with bias p ∈ [0, 1]. Instead of p, the problem can also be parameterized in terms of of log-odds log p 1−p , which can take any value in (−∞, ∞). Similar, one could use the following parametrization for quantum states guaranteed to give a positive semidefinite, Hermitian matrix with trace 1 Although this parametrization can certainly be advantageous for point estimation, it is unlikely to be helpful for uncertainty quantification: While X and ρ(X) carry equivalent information, the size of a region measured in "X-space" is hardly related to the size of a region in the physical state space. This is necessarily so, as any map from an unbounded space onto the compact quantum state space must grossly distort the geometry. So, having obtained a "small confidence region" in parameter space does not imply that the state has been well-estimated w.r.t. any physically relevant metric.

Confidence Regions from Linear Inversion
A particularly simple method to transform estimates of measurement data to estimates of quantum states is the method of linear inversion, which we are going to review now: First, assume that the true but unknown quantum state is represented by a d × d density matrix ̺ 0 and the QST performed by measuring m ≥ d 2 − 1 tomographicallycomplete measurement projectors E 1 , . . . , E m . By y k = tr (E k ̺ 0 ), k = 1, . . . , m we denote the (quantum) expectation values of E k for the true state ̺ 0 . Since these relations are linear, we can rewrite them as y = A̺, where ̺ stands for the quantum state interpreted as a vector and A is the measurement (or design) matrix independent of ̺. The desired (pseudo)-inverse of the above relation is and simplifies to Of course, in an experiment, the expectation values y are unknown and can only be approximated by some estimateŷ based on the observed data. The linear inversion estimate for the quantum state̺ is then given by Eq. (7) with the probabilities y replaced by the empirical frequenciesŷ. However, due to statistical fluctuations the estimated state̺ is not necessarily positive semidefinite [43], which led to the development of estimators enforcing the physical constraints such as the maximum likelihood estimator [44]. Although the linear inversion and maximum likelihood estimator solve two distinct problems -namely the unconstrained and constrained one, respectively -in certain cases the two are related. More precisely, if the outcomes approximately follow a Gaussian distribution, a fast projection algorithm computes the maximum likelihood estimate from the linear inversion estimate directly [45].
Here, we take a similar approach. First, the simple geometric interpretation of the linear inversion estimator (see Fig. 1) allows us to map confidence regions for the expectation values to confidence regions for the state without taking into account the positivity constraint: If Cŷ is a confidence region forŷ with confidence level α, then so is its pre-image under the measurement map for̺. Second, the truncation construction from Sec. 2.1 yields an improved confidence region for the problem with quantum constraints taken into account. As shown in Lemma 1, this approach yields admissible confidence region provided the original region was admissible.
The same construction can also be carried out for tomographically incomplete measurements, i.e. for m < d 2 : Since the measurement matrix A is non-invertible in this case, the estimate for the state̺ satisfying A̺ =ŷ is not uniquely defined. However, under additional structural assumptions, one can single out a unique estimate [14,16]. The singularity of the measurement map A also reflects in the confidence region defined by Eq. (8). Even if Cŷ is a bounded region, the confidence region for the state C̺ extends to infinity in the directions "unobserved by A". In both cases, the tomographically complete and incomplete, we can use the intersection with the psd states to reduce the the region's size while not sacrificing coverage. This improvement is especially far-reaching in the latter case, where it turns an unbounded region to a bounded one just by taking into account the physical constraints.
Of course, the question is whether we can somehow characterize the truncated confidence region C ∩ ̺ := A −1 (Cŷ) ∩ S + computationally efficiently. Since our goal is to show that this is an intractable problem exclusively due to the quantum constraints -and not because of the complexity of high-dimensional statistics in general -we are going to make the simplifying assumption that the measured frequencies are approximately Gaussian distributed. Furthermore, we are going to focus on a class of confidence regions that are efficiently characterizable in the unconstrained setting, namely Gaussian confidence ellipsoids or, more precisely, ellipsoidal balls of the form centered at the the empirical frequenciesŷ. The m × m, symmetric, positive semidefinite matrix B completely specifies the ellipsoidal shape of this confidence region. These are the natural generalizations of the well-known 2σ confidence intervals to multivariate Gaussian distributions. However, in the unconstrained setting, the ellipsoidal construction (9) is known to be admissible only for m = {1, 2} [41], while it is not admissible for m ≥ 3 [46] due to Stein's phenomenon [47]. Smaller confidence ellipsoids with the same coverage can be obtained by shifting the center slightly [48,49]. Furthermore, other constructions similar to an egg [50] or the non-convex Pascal limaçon [51] are known to outperform the standard ellipsoids. Nevertheless, non of these constructions is known to be optimal and, to the best of the author's knowledge, no optimal confidence region for multivariate Gaussians in higher dimensions is known.
But, since our discussion is focused on the question how the physical psd constraints can be used to improve confidence regions, we are still going to use the ellipsoids (9) as a tractable example: As we will prove later, it is impossible to characterize the truncated ellipsoids efficiently although they are fully described by a few parameters, namelyŷ and B in the unconstrained case. In other words, we show that there exists a trade-off between computational and statistical efficiency for the problem of determining "good" confidence regions in QST.
In the remainder of this section, we are going to discuss a useful parametrization for the aforementioned ellipsoids (9). To this end we use the fact that any d × d Hermitian matrix can be expanded in a basis formed by the identity 1l and d 2 − 1 traceless Hermitian matrices σ i , i = 1, . . . , d 2 − 1, normalized according to Tr(σ i σ j ) = 2δ ij . With the symbols σ i we associate here the most common choice of the basis elements [52] -explicitly provided in Appendix A -while any other σ ′ i = j O ji σ j , given in terms of an orthogonal d 2 − 1 dimensional matrix O, are valid alternatives. For d = 2 the choice stated in Appendix A is simply the Bloch basis of Pauli matrices: σ 1 ≡ σ x , σ 2 ≡ σ y and σ 3 ≡ σ z . In higher dimensions the matrices σ i maintain the Bloch basis structure: Let then their construction mimics σ Therefore, we are going to refer to the σ i as (generalized) Bloch representation.
We are in position to provide the first result of this paper, falling into the category of geometry of quantum states: Theorem 1. For the tomographically complete case m ≥ d 2 − 1, the pre-image under the measurement matrix of any confidence ellipsoid of the form (9) can be represented as where̺ is the linear inversion estimator, that is a Hermitian matrix with Tr̺ = 1, and the R i > 0 (i = 1, . . . , d 2 − 1) are the ellipsoid's radii in the directions given by 1 furnishes any orientation of the semi-major axes of the ellipsoid.
Proof. Note that whenever the sum has no limits specified (like in Eq. (11)), by default it extends from 1 to d 2 − 1. Let us parameterize both ̺ ∈ C and̺ in the Bloch representation with coordinates w i andŵ i , respectively: Since y = Tr (E̺), andŷ = Tr (E̺) we find where Q is a m × (d 2 − 1) matrix with elements Q ki = Tr (E k σ i ). In other words, the Bloch coordinates satisfy the same ellipsoid equation (9) as the measurement outcomes with B substituted by the d 2 − 1 dimensional square matrix B ′ = Q T BQ. Since B is symmetric and positive definite, the same holds for B ′ . Hence, In the last step of the proof we simply change the orientation of the basis to

Computational Intractability of Truncated Ellipsoids
Guided by the discussion from the previous section we now study the confidence region for the linear inversion QST defined as where C is given by the ellipsoid (11) for the tomographically complete case m = d 2 −1.
In this section, we are going to show that in contrast to the full ellipsoid C, the truncated ellipsoid C ∩ cannot be characterized computationally efficiently. This shows, for example, in the fact that there is no efficient algorithm to answer the following question: How much does taking into account the physical constraints reduce the size of the confidence region on a particular set of observed data? Note that we will not be concerned with properties of the region estimator but with a single instance corresponding to a fixed set of data. By abuse of notation, we are going to refer to these instances as C and C ∩ as well.
More precisely, we are concerned with the question if for a fixed ellipsoid C there is any reduction in size due to constraints in Eq. (15) or if C is fully contained in the set of psd states. For the precise formulation, we use the representation of ellipsoids from Thm. 1.
The main result of this section is the following statement on the computational complexity of the aforementioned problem.
As a consequence of Thm. 2, the problem of "characterizing" the truncated confidence ellipsoids C ∩ ̺ := A −1 (Cŷ)∩S + defined in Sec. 2.2 computationally is hard in general. By "characterizing" we mean computing any property of C ∩ ̺ that is sensitive to whether the truncation influences the original ellipsoid or not, e.g. computing the volume of the truncated ellipsoid or its distance to boundary of the quantum state space with high enough precision. Note, however, that there are also properties such as the diameter that might be unaffected by the truncation in certain special cases and, hence, their computational complexity cannot be classified using Thm. 2. Therefore, the more general problem of computing truncated confidence regions (without the Gaussian approximation) is hard as well since it subsumes Prob. 1.
Another consequence of the theorem concerns confidence regions for the constrained problem, which output "good regions" for the unconstrained problem when the constraints are not active: More precisely, it is extremely natural to use likelihood ratio-based ellipsoidal confidence regions for unconstrained Gaussian data although they cannot be optimal due to Stein's phenomenon. So it is natural to require any quantum region estimator to behave this way in the particular case that the likelihood function is concentrated well away from the boundary of state space. What Thm. 2 shows is that any region estimator subject to this criterion must necessarily solve NPhard problems.
Finally, the remainder of this section is dedicated to give some insight to the proof of the main theorem and to discuss a tractable solvable special case. The proof of Thm. 2 is inspired by a similar result due to Ben-Tal and Nemirovski [53] in robust optimization theory, who showed that the following problem is NP-complete.
Although the two problems are strongly related, the intractability result [53] cannot be applied directly to our tomography related problem due to the following crucial difference: The proof of NP-completeness of Prob. 2 deals with the case k = d(d − 1)/2 + 1 and a set of real symmetric matrices (A k ) k , which are not necessarily pairwise orthogonal to each other [53,Sec. 3.4.1]. However, in Prob. 1, the σ ′ i (i = 1, . . . , d 2 − 1) form an orthogonal basis of the space of complex Hermitian, traceless matrices. Hence, we need to adapt the original proof strategy to deal with the restrictions imposed by our tomography related problem.
Let us start the outline of the proof of Thm. 2 with a simplified example, when the ellipsoid in question is a ball, i.e. when R i = R for all i = 1, . . . , d 2 − 1. With no loss of generality, we can assume σ ′ i = σ i . The following Lemma, which is proven in Appendix C, provides an easily checkable, necessary, and sufficient condition to decide Prob 1 for this special case.
Lemma 2. Let C denote a ball parameterized according to Thm. 1 with with radii R i = R and midpoint̺. C is fully contained in the set of psd density matrices if and only if where mineig̺ denotes the smallest eigenvalue of̺.
The statement is a straightforward but interesting extension of the known result that the largest ball centered at the completely mixed state and fully contained in the set of psd density matrices has radius R max = 1 2d(d−1) . Intuitively, when the center of the ball is moved away from the completely mixed state, the allowed radii become smaller. This correction happens to be quantified by the smallest eigenvalue of the new center. In conclusions, spherical ellipsoids do not constitute hard instances of Problem 1 provided that the minimal eigenvalue of̺ can be computed efficiently with high enough accuracy.
However, it turns out that a slightly more complicated setting is already enough to proof the computational intractability. The ellipsoids under consideration still have their semi-major axes aligned with the generalized Bloch basis, that is we assume σ ′ i = σ i . The only change compared to the previous setting is the choice of radii. We consider the same radius R 1 for all directions generalizing the x-direction to higher dimensions and the distinct radius R 2 for the remaining directions: Recall i d defined in Eq. (10). Now, in order to prove the computational intractability of Problem 1, we use a reduction from the balanced sum problem, which is known to be NP-complete.
In case there is such a vector ψ one says that the instance a allows for a balanced sum partition because the sum of components of a labeled by ψ i = 1 is equal to the sum of components a i labeled by ψ i = −1. The main technical difficulty is now to identify the values of R 1 and R 2 as well asρ depending on an instance of the balanced sum problem a such that the corresponding ellipsoid C given by Thm. 1 contains an element with negative eigenvalues if and only if a has a balanced sum partition. For the details, please see Appendix B.

MVCR for Gaussians
We now turn to the question of minimal volume credible regions (MVCR) in the Bayesian framework: In the unconstrained case, Gaussian posteriors are one of the few examples of multivariate distributions, where the MVCR are simple geometric objects, namely ellipsoids. In practice, Gaussian posteriors arise in the following scenario: Consider a random vector X ∼ N (µ, Σ), where the covariance matrix Σ is known and we wish to estimate its mean µ. If we furthermore assume a Gaussian prior for the mean, the posterior will be Gaussian as well due to the fact that the Gaussian distribution is its own conjugated prior. This is one of the few cases, in which the Bayesian update as well as the computation of an optimal credible region can be carried out analytically. First, computing the parameters for the Gaussian posterior distribution can be done by means of linear Kalman filter update equations, see e.g. [54,Sec. 2.4]. Second, for the credible region, assume that after the Bayesian update, the posterior distribution of µ is parameterized by its mean θ ∈ R N and covariance matrix Σ ∈ R N ×N . Therefore, the posterior of µ has probability density where is the Mahalanobis distance and |Σ| denotes the determinant of Σ. As elaborated in Sec. 1.2, the MVCRs are exactly the highest posterior density sets as defined in Eq. (3). Therefore, the MVCR with credibility α for the density Gaussian (20) is given by This is an ellipsoid centered at θ with radius r α determined by the saturated credibility condition (2): By γ(·, ·) we denote the incomplete Γ-function and P (·, ·) is its normalized version. The above condition fixes r α uniquely since x → P ( N 2 , x) is strictly monotonic for any with credibility α C (yellow) lies completely inside the psd states and is, therefore, equal to the ellipsoid taking into account positivity E(r + α ) with credibility α (blue hatched). Right: Parts of the original ellipsoid E(r α C ) lie outside the psd states (blue). Hence, the ellipsoid that takes into account positivity E(r + α ) has to have a larger radius in order to achieve the sought for credibility.
N > 0. Hence, determining the MVCR for a multivariate Gaussian posterior with known mean and covariances reduces to computing the radius r α , which is formalized in the following problem.
An efficient algorithm for solving Prob. 4 is outlined in the following. To ease notation, we set x = r 2 α /2. (i) W.l.o.g. we can assume that α ≤ 0.9 (or some other arbitrary constant).
Otherwise, the problem can be restated in terms of Q( N 2 , x) = 1 − P ( N 2 , x), which allows for a similar analysis. The condition α ≤ 0.9 restricts the search space for x to some finite interval [0, t max ]. Note that the upper bound t max grows at worst polynomially in N 2 . (ii) The above restriction, the finite precision, and the fact that x → P ( N 2 , x) is strictly monotonic allow for interpreting the problem of finding x given α as a search in an ordered, finite list of size M ∼ tmax δ . (iii) Each entry of this list can be evaluated with exponential precision in polynomial time using a power series expansion of P ( N 2 , x) (for more details see Lemma 10 in Appendix D).
(iv) Since finding x in this list only requires log M evaluations using binary search, the whole problem can be solved in polynomial time.

Bayesian QST
Let us now turn to the application of Bayesian methods to QST, for a more thorough discussion see e.g. [55]. In order to incorporate the prior knowledge of positive semidefiniteness, we chose a prior that is concentrated on S + and vanishes on its complement. As before, we choose a (truncated) Gaussian prior and, therefore, Gaussian posteriors. Hence, the density π + θ,Σ (̺) of a Gaussian posterior on S + with respect to the flat Hilbert-Schmidt measure d̺ on S can be written as Here, π θ,Σ is the multivariate Gaussian from Eq. (20) with θ ∈ S. The other factors in Eq. (24) ensure that π + θ,Σ is a proper probability distribution supported on S + : χ(̺) is the indicator function of S + and C θ,Σ is the normalization constant defined by From now on we will drop the subscripts indicating the mean θ and the covariance matrix Σ if no confusion arises. It is then important to remember that the constant in question is denoted by C, while the credibility region is C. The problem we try to solve is the following: Given the mean θ, covariance matrix Σ, and credibility α, can we find the MVCR for the Gaussian distribution supported on S + ? Since the posterior density (24) is supported on the psd states and MVCRs are highest-density sets due to (3), the MVCR is of the form Similar to Eq. (23), the radius is determined by the credibility condition However, this case involves the normalization constant C from (24) and the integral is restricted to the psd states. Also, there is no closed-form analogue to Eq. (23) due to the psd constraint.

Computational Intractability
Our main result from this section concerns MVCR for Gaussian posteriors that are fully supported on the psd states. We will show that the following problem is computational hard.
In other words, there is no efficient algorithm that outputs smallest volume credibility regions for every Gaussian distribution on S restricted to the positive semidefinite states and every credibility α. Consequently, there cannot be an efficient algorithm to solve the problem of MVCR for QST, since the latter more general problem contains the instances of Prob. 5. To prove Prob. 5, we use a reduction from Problem 1, which has already been shown to be NP-hard. This reduction runs along the following lines: (i) Assume that Prob. 5 can be solved efficiently.
(ii) As we will prove later, every ellipsoid E * in S can be encoded as a minimum volume credible ellipsoid for some Gaussian distribution π with a suitable choice of θ, Σ, and R: Note that only θ is uniquely defined. Σ is defined only up to a multiplicative, positive constant, since every rescaling of Σ can be compensated by an appropriate rescaling of R.
(iii) Using the assumed efficient algorithm for Prob. 5, we can compute the normalization constant C of the truncated distribution (24) for given θ and Σ with sufficient precision in polynomial time.
(iv) Based on this, we can compute a credibility α such that R = r α C and, therefore, (v) The crucial observation is that this ellipsoid is contained in the psd states if and only if the corresponding MVCR for the truncated distribution π + fulfills See Fig. 2 for an illustration. Since we can compute r + α efficiently by assumption, checking Eq. (30) allows us to decide Prob. 1.
In conclusion, the main result from this section is the following lower bound on the computational complexity of Problem 4. Theorem 3. If Problem 5 has a polynomial time algorithm, then we can also decide Problem 1 in polynomial time. Therefore, there is no efficient algorithm for Problem 5 unless P = NP.
The proof runs along the lines outlined above and can be found in Appendix D.
Here, the main technical problem is that we are dealing with finite-precision arithmetic.

Conclusion & Outlook
The goal of this work is to provide an absolute "upper bound" on what we can expect from algorithms computing error regions for QST and to demonstrate that there is a trade-off between optimality and efficiency. This paper should not be understood as providing a no-go theorem for efficient algorithms in practice since the negative result of this work does not rule out efficient algorithms for practically acceptable approximations to optimal regions. Also, there is no indication that the various approaches used in practice give rise to regions that are far from optimal or do not have the advertised coverage. The reason our result leaves room for feasible approaches in practice are twofold: First, like any result showing NP-hardness, we prove that there is no efficient algorithm solving the exact problem deterministically for any instances. Hence, our result neither precludes the existence of efficient approximate or probabilistic algorithms, nor cannot make any statement about average case hardness. Second, although the experimental effort necessary for full-fledged tomography scales polynomially in the dimension of the system -and is, therefore, efficient in the sense of computational complexity -in practice other characterization techniques such as randomized benchmarking or direct fidelity estimation become more important for larger dimensions. It should now be the goal of future work to further close down the gap between existing positive results and the proven no-go theorems from either side.
More specifically, due to the simplifying assumptions made we investigate computational intractability that is solely caused by the quantum constraints and not by the general complications in high-dimensional statistics. In the Bayesian settings we show that minimal volume (w.r.t. the Hilbert-Schmidt measure) credible regions for truncated Gaussian posterior distributions are hard to compute. Therefore, the problem of determining MVCR for QST cannot be solved efficiently as well, since any algorithm solving the latter must also be able to solve instances with the specific prior used in Prob. 5.
The result for frequentist confidence regions is somewhat weaker since optimal confidence regions for high-dimensional Gaussian distributions are not known for most natural notions of optimality. Nevertheless, Gaussian confidence ellipsoids constitute a viable choice due to their simplicity and tractability. However, our results show that the constraints imposed by quantum mechanics render the task of characterizing the confidence regions for the constrained problem computationally intractable -even under the simplifying assumptions made. Of course, any more general setting encompassing the Gaussian approximation will be at least as hard to treat as the one used in this work. Furthermore, it also shows that computing any confidence region estimator yielding ellipsoids when the constraints are not active (and anything possibly better when they are) involves solving NP-hard problems.
Recently, the mathematical statistics community has started to analyze the tradeoffs between computational complexity and optimality in inference problems -see e.g. [56,57,58]. Early papers concentrated on the problem of sparse principal component analysis, which roughly asks whether the covariance matrix of a random vector possess a sparse eigenvector with large eigenvalue [56,57,58]. Later works have addressed the much better-studied problem of sparse inference [58]. The main difference between these papers and the present one is that we always condition on a data set and show that certain operations for quantifying uncertainty given the data are hard. This approach is canonical for a Bayesian analysis, but merely "natural" for orthodox error regions (c.f. Sec. 1.2). In contrast, Refs. [56,57,58] analyze the "global" performance of orthodox estimators -i.e. they do not require looking at worst-case scenarios over the data. References [56,57,58] achieve this by reducing a certain problem ("hidden clique") -that is conjectured to be hard in the average case -to the sparse PCA problem; while [58] employs a more subtle argument involving the non-uniform complexity class P/poly. It would be very interesting to adapt such arguments to the problem of quantum uncertainty quantification.
Of course, from the practical point of view, "positive" results -i.e. new algorithms to solve the problem -would be more beneficial. Here, recent work on sampling distributions restricted to convex bodies [59,60] could be a starting point for further investigations.
Beside quantum state tomography, our results might also be relevant to problems involving psd constraints such as the estimation of covariance matrices. We shall start the current discussion with a word of clarification concerning the dual notation already used in the definition of the Bloch vector. We utilize an alternative representation of the state |Ψ in terms of a complex vector ψ with coordinates specified with respect to the orthonormal basis fixed in Appendix A. Consequently, Ψ|Ψ is the norm of |Ψ , while ψ denotes the norm of ψ. Obviously both norms assume the same value.
In a first step of the proof we write down the positivity condition for the ellipsoid under investigation: The confidence ellipsoid C is fully contained in the set of psd states if and only if for all ρ ∈ C and all |Ψ , Ψ|ρ|Ψ ≥ 0. holds. In the parametrization from Thm. 1, this condition can be rewritten as where we have already restricted our attention to the special case from Eq. (18). Furthermore, we have used the shorthand v i (ψ) = Ψ|σ i |Ψ , which are the rescaled Bloch coordinates of the density matrix |Ψ Ψ|. Condition (B.2) is independent of the norm of |Ψ thus, we can fix Ψ|Ψ = d. Recall that Eq. (B.2) has to hold for all values of u with u T u ≤ 1. Since the left hand side assumes its minimal value for we find that Eq. (B.2) is equivalent to Using the unusual normalization of |Ψ , we find In the following, we restrict our attention to R 1 > R 2 , so that both term inside the square root are manifestly non-negative.
In the second step of the proof we show and utilize the following lemma: If̺ is a symmetric, real matrix w.r.t. |i , then the minimum of g(ψ) is attained by a vector ψ with real coordinates.
Proof. Note that we can decompose any vector |Ψ into its real and imaginary part where the |Ψ i are given by real vectors ψ i . Therefore, for̺ being real and symmetric, we find A similar equality holds with̺ replaced by 1l or σ i for i = 1, . . . , i d , since the latter matrices are symmetric and real as well. To shorten the notation, we now define two i d + 1 dimensional vectors x 1 and x 2 with components (α = 1, 2) (B.9) Since d = ψ 2 = ψ 1 2 + ψ 2 2 , we find where we used triangle inequality in the last step. Therefore g(ψ) ≥ g(ψ 1 ) + g(ψ 2 ) (B.11) so that if g(ψ) is non-negative for all real vectors, it is also non-negative for every complex vector ψ. More intuitively, the above result is true because the construction of g(ψ) utilizes only the generalized σ x Pauli matrices, which by construction pick up certain real parts of ψ * ⊗ ψ (imaginary contribution could appear only due to σ y ).
The next step of the proof, which is crucial for encoding an instance the balanced sum problem, is the choice of the ellipsoid's center. We choosê with q to be specified below and |a = k a k |k denoting a state represented by a real, integral vector a playing the role of the instance of Prob. 3. Since̺ given by Eq. (B.12) is manifestly real and symmetric we can restrict our attention to ψ ∈ R d due to Lemma 3. We find (B.14) Before we will be ready to take an advantage of the above encoding we need to perform a sequence of tedious algebraic manipulations. In short, the function we work with has an algebraic form g(ψ) = κ − √ ∆, with both κ and ∆ being non-negative. Testing if this function is non-negative is thus equivalent to checking the inequality κ 2 − ∆ ≥ 0. If we divide this inequality by 2(R 2 1 − R 2 2 ) and fix q = q + or q = q − with (B.15) we can rearrange it to the convenient form where: Both solutions (B.15) assure that (B.16) is free from additional terms proportional to (a · ψ) 2 , except those already hidden in f . Hence, the original problem of deciding whether the ellipsoid E centered at̺ and with radii (18) is contained in the psd states can be rephrased as deciding whether the maximum of the left hand side of Eq. (B.16) is smaller or equal to some constant: The relation of Problem 1 to the balanced sum problem (Problem 3) is derived in the following Lemma.
Lemma 4. If the instance a of Problem 3 allows for a balanced sum partition, then On the other hand, if there is no such partition, we have where p(x) = 2x 4 is a non-negative polynomial.
For the sake of clarity we relegate the proof of the above lemma to the end of this section. As a consequence of Lemma 4 the choice, implies that an efficient algorithm deciding whether the inequality (B.16) is satisfied or not is also capable of deciding Prob. 3 efficiently. This is exactly the statement of Thm. 2.
The last step we need to make is to find the parameters R 1 and R 2 leading to the choice (B.25). To this end, we set R 2 = ǫR 1 with 0 < ǫ < 1 and introduce two positive parameters (B.26) j the quantity f (ψ j ) is non-negative, so is the right hand side of Eq. (B.23). From (B.24) we can find the bound Rearranging Eq. (B.18), taking the square root and substituting (B.25) we can see that R 1 is implicitly defined by the relation If the left hand side of (B.29) happens to be bigger than 1/2, we need to take the q + solution on the right hand side (and q − in the opposite case). In order for the square roots in Eq. (B.29) to be real-valued, we need to assume The latter condition assures that q ± are real while the former condition, as it does not depend on R 1 , can be immediately solved for ǫ: However, Eq. (B.32) does not yield a universal bound for acceptable values of ǫ since B 1 depends on the particular instance a. To obtain a lower bound independent of a, we use Eq. (B.28), obtaining: Since both sides of (B.29) are non-negative, we can take the square of this relation and turn it it into a quadratic equation for R 1 . Surprisingly, this equation has a trivial solution R 1 = 0 (only relevant while dealing with q − ) and a single non-trivial solution which can be simplified to the form: The condition (B.31) becomes trivially satisfied, while the left hand side of Eq. (B.29) is greater than 1/2 (relevant for q + ) for In the opposite case the inequality is reversed. When (B.35) occurs, we find that while in the opposite case the parameters q + and q − swap. These interrelations between the parameters imply that regardless of the validity of (B.35), the solution (B.34) uniquely determines q initially introduced in (B.12) as given by the formula (B.36). This parameter is manifestly smaller than 1 and due to (B.32) it is also non-negative. With the given choice of parameters (B.34, B.35) and q specified as above, we complete the reduction of the balanced sum problem to Prob. 1. To finalize the proof of Theorem 2, we now state the proof of Lemma 4.
Proof of Lemma 4. The first part of the proof -Eq. (B.22) -follows from a simple calculation utilizing the partition vector ψ defined in (19). Note that as a · ψ = 0, we immediately obtain the first equality in (B.22), which since C 2 is non-negative turns into inequality in (B.23).
To prove (B.24), we define the set of all possible (2 d in total) partition vectors and (for an arbitrary 0 < λ < 1) the set of vectors that are "close" to some element from Z Because a ≥ 1, the set B can be thought of as a disjoint union of 2 d balls centered around the elements of Z. For further convenience we denotez = argmin z∈Z ψ − z , and δ := ψ −z. By constructionz k = sign ψ k so that for all k = 1, . . . , d Using all the above, the fact thatz 2 k = 1 andz 3 k =z k , and the Jensen inequality we can further estimate As a does not allow for a balanced sum partition and both,z and a are integral, we must necessarily have |a ·z| ≥ 1. Thus Taking all the above results together with |a · ψ| ≤ a ψ = a √ d we obtain (B.45) We will now study two cases. For ψ ∈ B, we have 0 ≤ δ ≤ λ/a, so that while for the opposite case (ψ / ∈ B), when δ > λ/a, one finds Therefore, we have for any ψ ∈ R d with ψ so that by setting λ = d −3/4 we obtain the desired result with p(ad) = 2(ad) 4 .

Appendix C. Proof of Lemma 2
To check whether a sphere with radius R centered at̺ is contained in the set of psd states, specialize Eq. (B.4) to the special case R 1 = R 2 : Since for any pure state |Ψ the identity holds (Bloch vectors of pure states live on the hypersphere), the inequality in question becomes Simple minimization with respect to |Ψ leads to the final result stated as Lemma 2.

Appendix D. Proof of Theorem 3
Let us now construct the polynomial time reduction of Prob. 1 to Prob. 4. We will begin with the main observation of this reduction, namely Eq. (30).
Lemma 5. Let π(̺) denote a Gaussian distribution on S and π + (̺) = Cπ(̺)χ(̺) the corresponding restricted Gaussian with the same mean and covariance matrix, as defined in Eq. (24). For any α ∈ [0, 1], the credible ellipsoid E(r α C ) with credibility α C is contained in the psd if and only if the credible ellipsoid for π + , E(r + α ), with credibility α has the same radius, that is Eq. (30) holds.
Proof. The two cases of E(r α C ) being contained and not being contained in the psd states are illustrated in Fig. 2. First, assume that E(r α C ) ⊂ S + , then Note that the right equation is exactly the defining Eq. (27) for the positive radius r + α if r + α = r α C . Now, assume that a part of the ellipsoid O = E(r α C ) \ S + = ∅ lies outside the psd states. Then, as can be seen on the right side of Fig. 2, we need to enlarge r + α to compensate for the lost probability weight of O. The latter cannot be vanishing, since the Gaussian density π(̺) is strictly positive. Therefore, r + α > r α C in this case. Of course, the difference between r α C and r + α may in general become too small to be efficiently detectable. However, we will show that for the instances of the balanced sum problem encoded in Problem 1, this is not the case. A first step toward this is the following Lemma. Lemma 6. Let a ∈ N d describe an instance of the balanced sum problem and the corresponding encoding ellipsoid for Problem 1 defined in Appendix B. There exists a polynomialp such that if E a is not a subset of S + , there is an element ̺ ∈ E a with mineig(̺) ≤ −p( a ) −1 < 0.
By tracing back the steps which lead to this equation, we find for |Ψ : Due to the special choice for ̺ 0 in (B.12) and a · ψ = 0, we have with q defined in (B.15). Therefore, we can rewrite Eq. (D.5) as where we have used Since all the constants on the right hand side of Eq. (D.9) can be expressed as polynomials in the input, it defines the polynomialp( a ) of the lemma. The left hand side of that equation is equal to Ψ|̺|Ψ , where for the special choice of u from (B.2). The claim of the lemma follows for this ̺ using Eq. (D.9).
We will now show how the explicitly parameterized ellipsoid (D.2) can be encoded as a MVCR-ellipsoid of a Gaussian distribution.

Lemma 7. Denote by
an ellipsoid E * ⊂ S, which is axis-aligned with the coordinate axes defined by the generalized Pauli operators.
Then, E * can be encoded as a α C MVCR-ellipsoid for a Gaussian distribution with mean ̺ 0 ∈ S + and covariance matrix Σ. The latter is diagonal in the generalized Bloch basis σ i with entries Σ ij = R 2 i δ ij and for the corresponding radius we have r α C = √ 2. Hence, the credibility is given by which can be calculated efficiently up to exponential precision for given C and N . psd . Same as Fig. 2 (right). Note that the solid blue and hatched blue regions need to have the same volume.
Proof. Since the generalized Pauli operators form an orthogonal system with tr σ i σ j = 2δ ij , we find for ̺ ∈ E * (D.14) Therefore, E * = E( √ 2) with mean ̺ 0 and the stated covariance matrix. The efficient computation of the credibility (D.13) is given later in the proof of Lemma 9.
Based on the gap proven in Lemma 6, we will now turn to the following question: In case Eq. (30) does not hold -that is the corresponding ellipsoid is not fully contained in the psd states -is the corresponding gap always large enough to be efficiently detectable?
Lemma 8. Let a ∈ N d be an instance of the balanced sum problem and denote by E a the corresponding encoding ellipsoid as given by Eq. (D.2). Furthermore, denote by π ̺0,Σ the Gaussian density, which encodes E a = E(r α C ) as an α C credible region as given by Lemma 7. Assume that a has a balanced sum partition and, therefore, E a is not a subset of S + .
Then, there exists a polynomial p such that Here, a 1 = k |a k |. In words, the gap of violation of Eq. (30) can only become polynomially small in the logarithm of the size of the problem specification.
Proof. First, let us lower bound the volume of E(r α C ) that lies outside the psd states (the solid blue region in Fig. D1). From Lemma 6 we know, that there exists a ̺ ∈ E(r α C ) with smallest eigenvalue smaller than −p( a ) −1 for some polynomialp. This also gives us a lower bound on dist(̺, S + ) = inf From [62, Theorem III.2.8] we know that for every ̺ + ∈ S + the following bound holds: This allows us to lower bound the volume of E(r α C ) that lies outside the psd states by an ellipsoid with the same covariance, but radius (2p( a ) maxeig (Σ)) since the solid blue and hatched blue regions in Fig. D1   Finally, note that the following crude inequality holds for x ≥ y, since the integrand is less than 1. Therefore, with Eq. (D.25) which proofs the claim.
We now turn to the problem of computing the normalization constant C for the restricted Gaussian distribution (24). First, we efficiently compute a credibility α ′ ∈ [0, 1] such that the corresponding credible ellipsoid E(r α ′ C ) is guaranteed to be contained in the psd states without knowing the value of C. This allows us to leverage Eq. (30) to compute C. Lemma 9. Let a ∈ N d be an instance of the balanced sum problem and denote by E a the corresponding encoding ellipsoid as defined by Eq. (D.2). Denote by π ̺0,Σ the Gaussian density, which encodes E a as an α credible region according to Lemma 7. Then, the ellipsoid E(r) is fully contained in the psd states provided Proof. We know that for any ̺ ∈ E(r) with r fulfilling (D.28) the following inequalities hold mineig ̺ 0 since mineig Σ −1 = (maxeig Σ) −1 . Therefore, E(r) ⊂ S + due to Lemma 2.
Lemma 10. Using the same notation as Lem. 9 and assuming Prob. 5 can be solved efficiently. Then, for every instance a of the balanced sum problem and the corresponding ̺ 0 , Σ, we can efficiently approximate the normalization constant C of π + ̺0,Σ with exponentially small error. More precisely, we have C =C(1 + ǫ), (D. 29) whereC can be computed in polynomial time making the correction term ǫ exponentially small.
Proof. Due to Lemma 9 and mineig ̺ 0 > 0, we can always find an r > 0 such that E(r) is fully contained in the psd. Indeed, the eigenvalues of ̺ 0 and Σ are readily calculated because of their particular simple form in Eq. (B.12) and Lemma 7: (D.31) Since we can choose r as small as we want, we may assume that x = r 2 2 ≪ 1 < N 2 . In this regime, we can expand the normalized incomplete Γ-function P in a power series [63] P N 2 , x = Since x ≪ 1, the term x k0 ensures that we can make the error in computing α exponentially small using only polynomial time in evaluating P k0 ( N 2 , x).
Assume that we have computedα = α − ǫ for some truncation error ǫ = R k0 ( N 2 , x) > 0. We may now use the (postulated) efficient algorithm for Prob. 5 to compute the radius of the manifestly positive MVCR r + α and, hence, using Eq. (30) the normalization constant: Since C > 1, we have with r α = r Therefore, the ellipsoid with radius r + α is also contained in the psd states. The same holds true if we replace r + α by the actual output r + α ± δ of the postulated efficient algorithm for Prob. 4 Here, δ denotes the (selectable) accuracy. By choosing δ small enough and possibly replacing the original radius r by r − δ, we can ensure that By assumption we can make both ǫ ′ and δ exponentially small using only polynomial time while P k0 ( N 2 , x) ↑ P ( N 2 , x) for k 0 → ∞, the correction tõ We now have all the necessary parts for the proof of the main theorem, which will conclude this section.
Proof of Thm. 3. The proof follows the outline stated at the beginning of this section: First, we encode the ellipsoid of Problem 1 to be checked as a MVCR of a Gaussian with mean ̺ 0 and covariance matrix Σ according to Lemma 7. Using Lemma 10, we compute an estimateC to the normalization constant C. Using the techniques from the proof of the aforementioned Lemma, we may compute an estimate α = C P N 2 , 1 =C(1 + ǫ) P k0 N 2 , 1 + ǫ ′ =α + ǫ ′′ . (D. 46) This can be done for exponential small errors ǫ, ǫ ′ in polynomial time. Here, the computable value is given byα =C P k0 An exponential small difference of α andα also implies an exponential small difference of r + 1−α and r + α : Set x := r + α andx := r + α and assume x >x -the opposite case can be treated along the same lines by choosing a larger constant as a bound forx. Following Eq. (D.25), we have