Efficient Bayesian credible-region certification for quantum-state tomography

TL;DR

This paper introduces an efficient method for certifying Bayesian credible regions in quantum-state tomography by reformulating the theory, employing accelerated Monte Carlo sampling, and deriving analytical approximations to improve computational efficiency.

Contribution

It presents a reformulated Bayesian credible-region theory, an accelerated hit-and-run Monte Carlo sampling method, and analytical formulas for size and credibility estimation without heavy sampling.

Findings

Region-average quantities effectively characterize credible-region properties.

Accelerated hit-and-run sampling significantly reduces computational complexity.

Analytical approximations accurately estimate size and credibility without Monte Carlo.

Abstract

Standard Bayesian credible-region theory for constructing an error region on the unique estimator of an unknown state in general quantum-state tomography to calculate its size and credibility relies on heavy Monte~Carlo sampling of the state space followed by sample rejection. This conventional method typically gives negligible yield for very small error regions originating from large datasets. We propose an operational reformulated theory to compute both size and credibility from region-average quantities that in principle convey information about behavior of these two properties as the credible-region changes. We next suggest the accelerated hit-and-run Monte~Carlo sampling, customized to the construction of Bayesian error-regions, to efficiently compute region-average quantities, and provide its complexity estimates for quantum states. Finally by understanding size as the region-average distance between two states in the region (measured for instance with either the Hilbert-Schmidt, trace-class or Bures distance), we derive approximation formulas to analytically estimate both distance-induced size and credibility under the pseudo-Bloch parametrization without resorting to any Monte~Carlo computation.

Equations

$$S_{\lambda }\equiv$$

$$C_{\lambda }\equiv$$

$$C_{\lambda }=\frac{\lambda S_{\lambda }+{\int }_{\lambda }^1\mathrm{d}{\lambda }^{\prime }{S}_{{\lambda }^{\prime }}}{{\int }_0^1\mathrm{d}{\lambda }^{\prime }{S}_{{\lambda }^{\prime }}}.$$

$$\log L(\mathbb{D}|{\mathbf{r}}^{\prime })\approx \log L_{\text{max}}-\frac{1}{2}({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})\cdot {\mathbf{F}}_{\text{textsc}ml}\cdot ({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})$$

$$\log L(\mathbb{D}|{\mathbf{r}}^{\prime })\approx$$

$$-\frac{1}{2}({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})\cdot {\mathbf{F}}_{\text{textsc}ml}\cdot ({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})$$

$$=$$

$$u_{\lambda }\equiv {\overline{q_{\lambda }({\mathbf{r}}^{\prime })}}^{{\mathcal{R}}_{\lambda }}=\frac{{\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })q_{\lambda }({\mathbf{r}}^{\prime })}{{\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })}=\frac{1}{K_{\text{smp}}}\sum _{l=1}^{K_{\text{smp}}}{q}_{\lambda ,l}$$

$$\text{updelta}{u}_{\lambda }=\left(\frac{1}{S_{\lambda }}-\frac{\text{updelta}{S}_{\lambda }}{S_{\lambda }^2}\right){\int }_{{\mathcal{R}}_{\text{updelta}\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })q_{\lambda }-\frac{\text{updelta}{S}_{\lambda }}{S_{\lambda }^2}{\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })q_{\lambda }$$

$$\frac{\partial u_{\lambda }{S}_{\lambda }}{\partial \lambda }=$$

$$+$$

$$\frac{\partial }{\partial \lambda }(u_{\lambda }{S}_{\lambda })=-S_{\lambda }\frac{\partial }{\partial \lambda }f(\lambda L_{\text{max}}).$$

$$\frac{\partial y_{\lambda }}{\partial \lambda }=-\frac{y_{\lambda }}{u_{\lambda }}\frac{\partial }{\partial \lambda }f(\lambda L_{\text{max}}).$$

$$y_{{\lambda }_{j+1}}=y_{{\lambda }_j}-\frac{y_{{\lambda }_j}}{u_{{\lambda }_j}}\frac{\partial }{\partial {\lambda }_j}f({\lambda }_j{L}_{\text{max}})$$

$$S_{\mathcal{D},\lambda }\equiv {\overline{\mathcal{D}({\mathbf{r}}^{\prime },{\hat{\mathbf{r}}}_{\text{textsc}ml})}}^{{\mathcal{R}}_{\lambda }}=\frac{{\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })\mathcal{D}({\mathbf{r}}^{\prime },{\hat{\mathbf{r}}}_{\text{textsc}ml})}{{\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })}$$

$${\mathcal{D}}_{\text{textsc}hs}=\mathrm{tr}\left\{({\rho }^{\prime }-{\hat{\rho }}_{\text{textsc}ml}{)}^2\right\}=({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml}{)}^2,$$

$${\mathcal{D}}_{\text{tr}}=\mathrm{tr}\left\{\sqrt{({\rho }^{\prime }-{\hat{\rho }}_{\text{textsc}ml}{)}^2}\right\}$$

$$\mathcal{F}=\mathrm{tr}{\left\{\sqrt{\sqrt{{\hat{\rho }}_{\text{textsc}ml}}{\rho }^{\prime }\sqrt{{\hat{\rho }}_{\text{textsc}ml}}}\right\}}^2$$

$$S_{\mathcal{D},\lambda }\to \frac{\int (\mathrm{d}{\mathbf{r}}^{\prime \prime })\mathcal{D}\eta (1-{{\mathbf{r}}^{\prime \prime }}^{\text{textsc}\mathrm{t}}{\mathbf{F}}_{\text{textsc}ml}{\mathbf{r}}^{\prime \prime }/(-2\log \lambda ))\prod _j\eta ({\mathbf{w}}_j^{\text{textsc}t}{\mathbf{r}}^{\prime \prime })}{\int (\mathrm{d}{\mathbf{r}}^{\prime \prime })\eta (1-{{\mathbf{r}}^{\prime \prime }}^{\text{textsc}\mathrm{t}}{\mathbf{F}}_{\text{textsc}ml}{\mathbf{r}}^{\prime \prime }/(-2\log \lambda ))\prod _j\eta ({\mathbf{w}}_j^{\text{textsc}t}{\mathbf{r}}^{\prime \prime })}$$

$${\int }_{{\mathcal{R}}_0}(\mathrm{d}{\mathbf{r}}^{\prime })\eta (L-\lambda L_{\text{max}})\cdots \to \frac{1}{V_{{\mathcal{R}}_0}}{\int }_{\text{all space}}\prod _j\mathrm{d}{r}_j^{\prime }\cdots$$

$$S_{\text{tr}}\approx \frac{8\sqrt{D{S}_{\text{textsc}hs}}}{3\pi }$$

$${\mathcal{F}}_{\text{A}}\approx 1-\frac{1}{2}{({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})}^{\text{textsc}\mathrm{t}}{\mathbf{Q}}_D({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})$$

$${\mathcal{F}}_{\text{textsc}b}\approx$$

$$+\frac{1}{2}{({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})}^{\text{textsc}\mathrm{t}}\left(\frac{1}{2}\mathrm{tr}\left\{P_r\mathbf{\Omega }\right\}\mathrm{tr}\left\{P_r{\mathbf{\Omega }}^{\text{textsc}\mathrm{t}}\right\}-{\mathbf{Q}}_r\right)({\mathbf{r}}^{\prime }-{\hat{\mathbf{r}}}_{\text{textsc}ml})$$

$${\mathbf{Q}}_r=\sum _{j=1}^r\sum _{k=1}^r\frac{⟨{\lambda }_j|\mathbf{\Omega }|{\lambda }_k⟩⟨{\lambda }_k|{\mathbf{\Omega }}^{\text{textsc}\mathrm{t}}|{\lambda }_j⟩}{{\lambda }_j+{\lambda }_k}.$$

$${\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })=$$

$${\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime })\to$$

$$=$$

$${\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime }){\mathbf{\Delta }}_{\text{textsc}ml}^{\prime }{{\mathbf{\Delta }}^{\prime }}_{\text{textsc}ml}^{\text{textsc}\mathrm{t}}\to$$

$$\mathbf{I}=$$

$$=$$

$${\int }_{{\mathcal{R}}_{\lambda }}(\mathrm{d}{\mathbf{r}}^{\prime }){\mathbf{\Delta }}_{\text{textsc}ml}^{\prime }{{\mathbf{\Delta }}^{\prime }}_{\text{textsc}ml}^{\text{textsc}\mathrm{t}}\to \frac{V_d}{V_{{\mathcal{R}}_0}}(-2\log \lambda {)}^{d/2+1}\det {\left\{{\mathbf{F}}_{\text{textsc}ml}\right\}}^{-1/2}{\mathbf{F}}_{\text{textsc}ml}^{-1}.$$

$$S_{\text{textsc}hs,\lambda }^{\text{ (A)}}$$

Full text

I Introduction

All physical-quantity estimates obtained from collected data should be accompanied by “error-bars” to accurately convey all properties of the physical system of interest. This applies to quantum-state tomography Smithey et al. (1993); Chuang and Nielsen (2000); Řeháček et al. (2007); Teo et al. (2011); Zhu (2014), which is an important preliminary step for implementing all quantum cryptography and computation protocols Ladd et al. (2010); Campbell et al. (2017); Ladd et al. (2010); Lekitsch et al. (2017) reliably.

Bootstrapping procedures Efron and Tibshirani (1993); Davison and Hinkley (1997) are amongst some of the most widely-used techniques for assigning “error-bars” to reconstructed quantum states. Recently, it was pointed out in Suess et al. (2017) that such assignments lack rigorous statistical foundations and may produce “error-bars” that are too small for reliable conclusions. The rather more justified approach falls under the study of hypothesis testing Christensen (2005). Two grand schools of thought exist for this purpose. In the context of quantum-state reconstruction of an unknown state $\rho$, one may treat $\rho$ as “absolute” (the frequentist school) and attempt to extract this knowledge from collected data. This suggests the constructions of confidence regions Christandl and Renner (2012); Blume-Kohout (2012); Faist and Renner (2016), which are error-regions for the state estimator $\widehat{\rho}$ from all plausible datasets, including those unmeasured in the experiment. An accurate $\widehat{\rho}$ for $\rho$ would then entail a collection of typically small confidence regions with high probability that $\rho$ for each plausible dataset lies in the corresponding region.

Given that only one dataset (the measured one) is really available to the observer, we shall focus on the apt Bayesian school of thought that instead regards this dataset as “absolute” and constructs credible regions Shang et al. (2013); Li et al. (2016) as the error regions for $\widehat{\rho}$ beginning with some prior distribution $p(\rho)$ of $\rho$. A fairly accurate estimator $\widehat{\rho}$ for some unknown $\rho$ naturally implies a credible region (generated from the measured dataset) of a small size with a large probability that this $\rho$ is inside the region—a high credibility Teo et al. (2018). In order to obtain a reasonably small error region (be it that of credible- or confidence-type), one may either resort to adaptive strategies Oh et al. (2018) and optimize additional properties of the region, or simply increase the dataset collected in quantum tomography.

The complicated quantum state-space boundary Bengtsson and Życzkowski (2006); Bengtsson et al. (2013) renders any analytical attempt at calculating size and credibility for any credible region futile, leaving numerical Monte Carlo (MC) methods Shang et al. (2015); Seah et al. (2015) as the only viable option. As the size of the credible region is defined as its volume fraction with the quantum state space, one needs an extremely large sample of the state space to finally end up having a reasonable sample for the region. Despite the optimistic advantages that some of these MC schemes may have in generating samples of arbitrary distributions, one deleterious issue for such an MC-filtering strategy becomes apparent when the dataset is large, which is the common situation in any tomography experiment. The resulting credible region eventually becomes too tiny relative to the quantum state space for MC-filtering to produce any effective yield to properly compute the size and credibility.

To solve this problem, we introduced the in-region sampling technique Oh et al. (2019) to feasibly compute all credible-region properties by simply sampling an appropriate quantity over the region itself. This follows from the logic that a change in the region-average quantity encodes the change in both region size and credibility. In our theory, we prove the central lemma stating that the size (and credibility) of any credible region are related to a class of region-average quantities through a first-order differential equation that is solvable numerically. As an example study, we discuss the computation of region-average quantities using the accelerated hit-and-run algorithm, its correlation properties, and estimate its complexity from geometrical considerations of the quantum state space. The region-average formalism encourages the formulation of the region capacity (a different way of stating “how big” a region is) by investigating the average relative distance between two points in the region. This region-average distance may be induced by any of the common measures used in quantum information, and we shall explicitly consider the Hilbert-Schmidt, trace-class, and Bures distance measures as popular examples. It turns out that this perspective, together with the results in Teo et al. (2018), offers closed-form analytical approximation formulas for an alternative rapid approximate Bayesian error certification with no Monte Carlo methods necessary.

This article is organized as follows: A preliminary introduction to the basic notions of standard Bayesian credible-region theory shall ensue in Sec. II, and the stage for discussions with large data is set in Sec. III. Next, we present our in-region sampling theory for size and credibility that works for any kind of data and prior in Sec. IV. Afterwards, we describe how region-average quantities can be numerically computed and estimate the computational complexities in Sec. V. Section VI then proceeds to quantify the region capacity in terms of region-average distances induced by all the three aforementioned distance measures other than the Bayesian size. Finally, for fast analytical Bayesian error estimates, we derive asymptotic formulas for all important region-average quantities in Sec. VII based on the perspective of distance-induced size. All numerical results and computation correlation properties are then presented and discussed in Sec. VIII with explicit examples in quantum tomography.

II Standard Bayesian credible-region theory

Before a quantum-state tomography experiment commences, the observer might have some (justifiable) preconception about the unknown quantum state $\rho\geq 0$ ($\mathrm{tr}\!\left\{\rho\right\}=1$) of Hilbert-space dimension $D$. Such preconception is usually not uniquely privileged, and therefore weighted with some prior probability distribution $p(\rho)$. After the experiment, the observer collects a set of data $\mathbb{D}$ that are informationally complete (IC) such that a unique estimator $\widehat{\rho}$ for $\rho$ is acquired. In quantum theory, the measurements are modeled as a probability-operator measurement (POM) consisting of a set of $M$ positive operators $\Pi_{j}\geq 0$ that sum to the identity. Associated to every such experiment is the likelihood function $L=L(\mathbb{D}|\widehat{\rho})$, with which the observer obtains a posterior probability distribution (knowledge after-the-fact) that is a function of $L$.

It was formerly established in Shang et al. (2013) that for this measured dataset $\mathbb{D}$, if $\widehat{\rho}$ is taken to be the estimator that maximizes $L$—the maximum-likelihood (ML) estimator—, then a Bayesian credible region (CR) $\mathcal{R}$ can be constructed around $\widehat{\rho}_{\textsc{ml}}$, which turns out to have a constant likelihood boundary $\partial{\mathcal{R}}$ within the quantum state space $\mathcal{R}_{0}$. For this CR, which is a subregion of $\mathcal{R}_{0}$, we can specify its size and credibility, the latter which is the probability that $\rho\in\mathcal{R}$. Such a region is optimal in the sense that it gives the largest credibility for a given size, or equivalently possesses the smallest size for a given credibility.

In this article, we shall be interested in reconstructing the $(d=D^{2}-1)$-dimensional real vectorial parameter $\boldsymbol{r}\leftrightarrow\rho$ that characterizes $\rho$. More technically, this equivalent parametrization is achieved with a Hermitian operator basis $\{1/\sqrt{D},\Omega_{j}\}^{d}_{j=1}$ that contains $d$ trace-orthonormal traceless operators $\Omega_{j}$ $\left(\mathrm{tr}\!\left\{\Omega_{j}\Omega_{k}\right\}=\delta_{j,k}\right)$, by which $\boldsymbol{r}=\mathrm{tr}\!\left\{\rho\,\boldsymbol{\Omega}\right\}$ is defined from the column $\boldsymbol{\Omega}$ of $\Omega_{j}$s. Formally, in terms of the multivariate parameter $\boldsymbol{r}$, the size and credibility of $\mathcal{R}=\mathcal{R}_{\lambda}$ for some $0\leq\lambda\leq 1$ are respectively given by Shang et al. (2013)

[TABLE]

where the volume measure $(\mathrm{d}\,{\boldsymbol{r}})$ incorporates the prior distribution $p(\boldsymbol{r})$ for $\boldsymbol{r}$, $\eta$ is the Heaviside function, $L(\mathbb{D})=\int_{\mathcal{R}_{0}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,L(\mathbb{D}|\boldsymbol{r}^{\prime})$. The important variable $0\leq\lambda\leq 1$ specifies the shape and size of $\mathcal{R}_{\lambda}$, from which the limits $\mathcal{R}_{\lambda=0}=\mathcal{R}_{0}$ and $\mathcal{R}_{\lambda=1}=\{\widehat{\boldsymbol{r}}_{\textsc{ml}}\}$ are immediate. Here, we note that the probability parametrization was adopted in Shang et al. (2013). Upon the condition that each datum measurement, corresponding to an outcome $\Pi_{j}$, the inherent statistics of $\mathbb{D}$ is therefore multinomial and the log-likelihood reads $\log L=\sum_{j}n_{j}\log p^{\prime}_{j}$ with the collected relative frequencies $\sum_{j}n_{j}=N$ that make up $N$ measured data copies, and $p^{\prime}_{j}=\mathrm{tr}\!\left\{\rho^{\prime}\Pi_{j}\right\}$ for any state $\rho^{\prime}$.

We can gain a clear physical picture of both size and credibility: they respectively quantify the prior and posterior content $\mathcal{R}$, hence the symbol $S_{\lambda}$ for the former. Owing to the dual nature of size $S_{\lambda}$ and credibility $C_{\lambda}$, it can be shown, in fact, that

[TABLE]

Put differently, $C_{\lambda}$ may be straightforwardly computable through single-parameter integrations in $\lambda$ so long as $S_{\lambda}$ is known up to some arbitrary constant multiple.

Nonetheless, the complicated boundary $\partial\mathcal{R}_{0}$ of the quantum state space makes the computation of $S_{\lambda}$ extremely challenging even numerically. The innate definition of $S_{\lambda}$, namely the volume fraction of $\mathcal{R}_{\lambda}$ to $\mathcal{R}_{0}$, requires, first, the generation of a sufficiently large sample of $\mathcal{R}_{0}$, followed by the filtering of all its sampled points that lie inside $\mathcal{R}_{\lambda}$ for any $\lambda$. There exist various Monte Carlo (MC) methods to sample $\mathcal{R}_{0}$ Shang et al. (2015). Ultimately, this MC-filtering strategy exhibits one major disadvantage: in the limit of large data sample ($N\gg 1$), $\mathcal{R}_{\lambda}$ would become so small relative to $\mathcal{R}_{0}$ that the MC-filtering strategy needs a sufficiently large number of random MC sample points from $\mathcal{R}_{0}$ to produce any useful yield. The scaling of MC sample size needed to maintain a fixed yield, which was estimated to be $O(N^{d/2})$ Teo et al. (2018), thus outgrows the feasible computational yield-rate very quickly. The bottom-line: a much more feasible numerical strategy to perform Bayesian error certification is necessary in this practical data limit.

III The large-data condition

Before presenting an alternative operational theory, unless otherwise stated, we shall consider $N\gg 1$ as the putative limit in pragmatic tomography experiments. We emphasize here that $N$ only has to be sufficiently large for the statistical central limit theorem to dictate a Gaussian form for $L$. In this limit, there can only be one of two cases: either $\mathcal{R}$ is completely inside $\mathcal{R}_{0}$ that contains a full-rank $\widehat{\rho}_{\textsc{ml}}$ (Case A) or partially truncated by the state-space boundary $\partial\mathcal{R}_{0}$ of $\mathcal{R}_{0}$ that houses a rank-deficient $\widehat{\rho}_{\textsc{ml}}$ (Case B) (see Fig. 1).

Case A arises when the unknown state $\rho$ is away from $\partial\mathcal{R}_{0}$, so that a sufficiently large $N$ would produce untruncated regions for $\lambda$ values corresponding to desirably large $C_{\lambda}<1$. This case offers a simple geometrical description for $\mathcal{R}$. Upon invoking the Taylor expansion of

[TABLE]

about the interior $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ up to the second order, we find that the likelihood $L$ is essentially a Gaussian function centered at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ of height $L_{\text{max}}$, with its covariance profile shaped by $\boldsymbol{F}_{\textsc{ml}}$, that is the Fisher information evaluated at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$. The CRs $\mathcal{R}_{\lambda}$ that go with this Gaussian likelihood are, hence, simple hyperellipsoids $\mathcal{E}_{\lambda}$ described by the inequality $(\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}})\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}(\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}})\leq-2\log\lambda$.

If $\boldsymbol{r}$ is located in $\partial\mathcal{R}_{0}$, then as $N$ increases, the ML estimator $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ would eventually approach $\boldsymbol{r}$ and there is a high probability that $\widehat{\boldsymbol{r}}_{\textsc{ml}}\in\partial\mathcal{R}_{0}$ before this happens. For sufficiently large $N$, we have Case B where $\partial\mathcal{R}\cap\partial\mathcal{R}_{0}$ is not disjointed and falls on the side of $\widehat{\boldsymbol{r}}_{\textsc{ml}}$. To asymptotically cope with such a situation, we may again expand

[TABLE]

about the boundary $\widehat{\boldsymbol{r}}_{\textsc{ml}}$, where this time $L$ is a Gaussian function centered at $\boldsymbol{r}_{\text{c}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}+\boldsymbol{F}_{\textsc{ml}}^{-1}\bm{\cdot}\bm{g}_{\textsc{ml}}$ with $\bm{g}_{\textsc{ml}}=\partial\log L(\mathbb{D}|\boldsymbol{r}^{\prime})/\partial\boldsymbol{r}^{\prime}|_{\boldsymbol{r}^{\prime}=\widehat{\boldsymbol{r}}_{\textsc{ml}}}$, and possesses a height $L^{\prime}_{\text{max}}=L_{\text{max}}\,\exp(\bm{g}_{\textsc{ml}}\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}^{-1}\bm{\cdot}\bm{g}_{\textsc{ml}}/2)>L_{\text{max}}$. The covariance profile of this Gaussian function is still governed by $\boldsymbol{F}_{\textsc{ml}}$, which produces hyperellipsoids $\mathcal{E}^{\prime}_{\lambda}$ described according to $(\boldsymbol{r}^{\prime}-\boldsymbol{r}_{\text{c}})\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}(\boldsymbol{r}^{\prime}-\boldsymbol{r}_{\text{c}})\leq-2\log\lambda^{\prime}$ for an “effective $\lambda$” $\lambda^{\prime}$ defined by $2\log(\lambda/\lambda^{\prime})=\bm{g}_{\textsc{ml}}^{\textsc{t}}\boldsymbol{F}_{\textsc{ml}}^{-1}\bm{g}_{\textsc{ml}}$. The CR $\mathcal{R}_{\lambda}$ is then asymptotically $\mathcal{E}^{\prime}_{\lambda}\cap\mathcal{R}_{0}$.

We point out that there is an intermediate case in which $\mathcal{R}_{\lambda}=\mathcal{E}_{\lambda}$, centered at $\widehat{\boldsymbol{r}}_{\textsc{ml}}\notin\partial\mathcal{R}_{0}$, is truncated by $\partial\mathcal{R}_{0}$. Such a situation can happen when $N$ is not sufficiently large, and tends to either Case A or B as $N$ grows. On a separate note, Ref. Teo et al. (2018); Oh et al. (2018) explicitly studies also this intermediate case.

IV In-region sampling theory

Suppose we have a CR $\mathcal{R}_{\lambda}$, with which we define the average quantity

[TABLE]

for some function $q_{\lambda}$, which is approximately equivalent to the discrete-sum average of $q_{\lambda,l}$ values over a sufficiently large number $K_{\text{smp}}$ of region points. If we probe the response of $u_{\lambda}$ with an incremental change $\lambda\rightarrow\lambda+\updelta\lambda$ in $\lambda$ as in Fig. 2, the result is the total change

[TABLE]

after limiting all small changes to the first order, which reveals that a small increment $\updelta u_{\lambda}$ can be explained by a change $\updelta S_{\lambda}$ in size that is accompanied by the (in)exclusion of the annular sum $\int_{\mathcal{R}_{\updelta\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,q_{\lambda}$. Put simply, tracking the change in $u_{\lambda}$ allows us to infer how much $S_{\lambda}$ has changed.

To better utilize this intuition, we first take the derivative of $u_{\lambda}S_{\lambda}$, which gives

[TABLE]

after invoking the derivative identity $\mathrm{d}\eta(x)/\mathrm{d}x=\delta(x)$ between $\eta(x)$ and the Dirac delta function $\delta(x)$. Next, we may impose the following functional form $q_{\lambda}(\boldsymbol{r}^{\prime})\equiv f(L)-f(\lambda L_{\text{max}})$ for $q_{\lambda}$, where $f(L)$ is some arbitrary function of $L$. This simplifies Eq. (7) to

[TABLE]

We now have a first-order differential equation that describes the dynamics of $S_{\lambda}$ according to the parametric region-average $u_{\lambda}$. With the initial condition $S_{\lambda=0}=1$, the entire functional form of $S_{\lambda}$ can then be recovered with Eq. (8). This completes the constructive proof of our so-called

Region-average computation (RAC) lemma: For any prior $(\mathrm{d}\,\boldsymbol{r}^{\prime})$ and measurement data $\mathbb{D}$, the prior content $S_{\lambda}$ (up to a multiplicative factor), and hence the credibility $C_{\lambda}$, are all inferable from $u_{\lambda}$ defined in Eq. (5) with $q_{\lambda}(\boldsymbol{r}^{\prime})\equiv f(L)-f(\lambda L_{\text{max}})$.

To proceed, we first perform the substitution $y_{\lambda}=u_{\lambda}\,S_{\lambda}$ to yield another differential equation

[TABLE]

The solution to $y_{\lambda}$ can then be obtained numerically through Euler’s method Butcher (2003). In practice, we may start from $S_{\lambda\approx 0}\equiv 1$ and iterate

[TABLE]

for a sequence of discretized $\lambda\rightarrow\lambda_{j}$ values ranging from 0 to 1. For feasible computation of $u_{\lambda}$, we shall choose $f(L)=\log L$.

V Region-average numerical computation

V.1 The hit-and-run algorithm

The hit-and-run algorithm is a direct convex-body MC sampling scheme that generates random sample points in the body according to some predefined distribution. This algorithm is thus suited for sampling $\mathcal{R}_{\lambda}$ according to some prior distribution $p(\boldsymbol{r})$ for the unknown $\boldsymbol{r}$.

The sampling principles behind an efficient hit-and-run computation begin with defining the smallest possible convex set $\mathcal{B}\supseteq\mathcal{R}_{\lambda}$ that houses $\mathcal{R}_{\lambda}$ and has an easy-access geometry. Starting from a known point in $\mathcal{R}_{\lambda}$, say the ML estimator $\widehat{\boldsymbol{r}}_{\textsc{ml}}$, a random line segment passing through this point is generated, with its endpoints fixed at $\partial\mathcal{B}$ that are quickly computable because of its simple geometry. Following which, sampling commences by repeatedly picking a random point along this segment until it lies in $\mathcal{R}_{\lambda}$. This point is next taken to be the new reference point through which another line segment is generated to find a new random point in $\mathcal{R}_{\lambda}$, until a set of $K_{\text{smp}}$ points is gathered.

We can make use of the straightforward hyperellipsoidal characteristics inherent from the central limit theorem to construct $\mathcal{B}$. For Case A, where $\mathcal{R}_{\lambda}=\mathcal{E}_{\lambda}$, $\mathcal{B}$ can just be taken to be $\mathcal{E}_{\lambda}$ characterized by $\widetilde{\boldsymbol{F}_{\textsc{ml}}}=\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda)$ from the earlier discussions in Sec. III. We now turn to the more interesting and practically ubiquitous Case B, where the large-$N$ arguments of Sec. III imply that we may fix $\mathcal{B}=\mathcal{E}^{\prime}_{\lambda}$, the profile of which is governed by $\widetilde{\boldsymbol{F}_{\textsc{ml}}^{\prime}}=\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda^{\prime})$. For this case, if $\mathcal{B}$ is much larger than $\mathcal{R}_{\lambda}$, sampling the latter would incur a significant amount of wastage. Fortunately there exists an accelerated version of hit-and-run Kiatsupaibul et al. (2011) that adaptively shrinks the endpoints of the line segment to reduce the search space each step.

To check if the random point chosen from the line segment is in $\mathcal{R}_{\lambda}$, we recall that $\mathcal{R}_{\lambda}$ is equivalently the continuous set of unit-trace operators that are both positive and satisfy the hyperellipsoidal constraint defined by the inequality $L(\mathbb{D}|\boldsymbol{r}^{\prime})>\lambda\,L_{\text{max}}$ under the $N\gg 1$ limit. Since all points in $\mathcal{B}$ essentially fulfill the latter constraint, the primary task is to check if the random point correspond to a legitimate quantum state. It is known that the Cholesky decomposition Horn and Johnson (1985); Higham (2009), a routine that factorizes $\rho^{\prime}=A^{\dagger}A$ for any positive operator $\rho^{\prime}$, is an efficient and numerically stable way for the job, with a computational complexity of $O(D^{3})$. A routine failure implies that the operator corresponding to the selected point is not positive.

The complete pseudo-coded algorithm for the accelerated hit-and-run, tailored for an arbitrary prior distribution $p(\boldsymbol{r})$, is stated below Bélisle et al. (1993); Smith (1996); Kiatsupaibul et al. (2011):

**Accelerated hit-and-run for sampling $\mathcal{R}_{\lambda}$

** Beginning with $k=1$ and $\boldsymbol{r}_{\text{ref}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}$ of $N\gg 1$:

Generate a random line segment characterized by $\boldsymbol{y}=\boldsymbol{r}_{\text{ref}}+\mu\,\boldsymbol{\mathrm{e}_{v}}$, where $\boldsymbol{\mathrm{e}_{v}}=\boldsymbol{v}/|\boldsymbol{v}|$ and $\boldsymbol{v}$ follows the standard Gaussian distribution (mean 0 and variance 1 for each column entry). Its endpoints are parametrized by $\mu_{\pm}=[-b\pm\sqrt{b^{2}-a(c-1)}]/a$, where $\boldsymbol{\Delta}=\boldsymbol{r}_{\text{ref}}-\boldsymbol{r}_{\text{c}}$, $a={\boldsymbol{\mathrm{e}_{v}}}^{\mathrm{\,\textsc{t}}}\,\boldsymbol{A}\,\boldsymbol{\mathrm{e}_{v}}$, $b={\boldsymbol{\Delta}}^{\mathrm{\,\textsc{t}}}\,\boldsymbol{A}\,\boldsymbol{\mathrm{e}_{v}}$, $c={\boldsymbol{\Delta}}^{\mathrm{\,\textsc{t}}}\,\boldsymbol{A}\,\boldsymbol{\Delta}$, and $\boldsymbol{A}=\widetilde{\boldsymbol{F}_{\textsc{ml}}}\,\,\text{or}\,\,\widetilde{\boldsymbol{F}_{\textsc{ml}}}^{\prime}$. 2. 2.

Define $\beta_{1}\equiv\mu_{\text{min}}=\min\{\mu_{+},\mu_{-}\}$ and $\beta_{2}\equiv\mu_{\text{max}}=\max\{\mu_{+},\mu_{-}\}$. 3. 3.

Pick a random number $\beta_{1}\leq\beta\leq\beta_{2}$ according to the marginal probability distribution $p(\boldsymbol{r}_{\text{ref}}+\beta\,\boldsymbol{\mathrm{e}_{v}})/\int\mathrm{d}\beta^{\prime}p(\boldsymbol{r}_{\text{ref}}+\beta^{\prime}\,\boldsymbol{\mathrm{e}_{v}})$ truncated in the interval [$\beta_{1},\beta_{2}$] and obtain $\boldsymbol{r}_{\text{test}}=\boldsymbol{r}_{\text{ref}}+\beta\,\boldsymbol{\mathrm{e}_{v}}$. 4. 4.

Check whether $\rho_{\text{test}}\leftrightarrow\boldsymbol{r}_{\text{test}}$ is positive.

• •

If so, define $\boldsymbol{r}_{\text{ref}}=\boldsymbol{r}_{\text{test}}$, raise $k$ by 1, and go to step 1.

• •

If not, set $\beta_{1}=\beta$ if $\beta<0$, or $\beta_{2}=\beta$ if $\beta>0$, and repeat steps 3 and 4. 5. 5.

End routine if $k>K_{\text{smp}}$, the total number of sample points desired.

To further speed up the algorithm for Case B, one can assign $\mathcal{B}$ to be the hyperellipsoidal cap composed by a hyperplane that is tangent to the isoGaussian level curve of $\mathcal{E}^{\prime}_{\lambda}$ at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ and the part of $\mathcal{E}^{\prime}_{\lambda}$ below it (refer to Sec. VII.3). Numerical experience shows that this speed up is negligible in the presence of the endpoint adaptation mechanism of accelerated hit-and-run.

We end this introduction of hit-and-run by noting that a hyperellipsoidal $\mathcal{B}$ is constructed based on the large-$N$ limit, where the boundary $\partial\mathcal{R}_{\lambda}\cap\mathrm{int}\{\mathcal{R}_{0}\}$ of the physical region is well approximated by this hyperellipsoid. The highly skeptical may insist that, perhaps, for a finite $N$, even if $N\gg 1$, there might still be cases where a part of this boundary protrudes $\mathcal{B}$. To be on the safer side, one may choose a hyperellipsoidal $\mathcal{B}$ of a reasonably larger size (say doubled) than the one given by the central limit theorem. This will almost surely contain the physical error region with a much smaller failure probability. The pertinent question is: “Can we verify that $\mathcal{B}$ contains $\mathcal{R}$ with arbitrary precision?” The answer unfortunately is negative both in theory and practice. This is because a positive answer would entail a complete knowledge about $\partial\mathcal{R}$, obtaining which is either computationally not feasible in general, or an NP-hard problem in some context Suess et al. (2017).

V.2 Numerical complexity estimations

After suppressing dependences on logarithmic factors and error parameters, it was argued that the number of hit-and-run steps needed to gather enough sample points and form an ensemble described by $p(\boldsymbol{r})$ in hit-and-run is $O\!\left(d^{2}R^{2}_{\text{out}}/R^{2}_{\text{in}}\right)=O\!\left(D^{4}R^{2}_{\text{out}}/R^{2}_{\text{in}}\right)$ Cousins and Vempala (2016); Lovász (1999) in the limit $D\gg 2$, where $R_{\text{out}}$ is the radius of the smallest outer sphere that contains $\mathcal{R}_{\lambda}$ and $R_{\text{in}}$ is that of the largest inner sphere that can be inscribed in $\mathcal{R}_{\lambda}$. Together with the floating-point-operations complexity $O(D^{3})$ in a typical Cholesky decomposition algorithm Higham (2009), we have an estimate for the complexity $\mathrm{cmpl}=O(D^{7}R^{2}_{\text{out}}/R^{2}_{\text{in}})$ for the entire hit-and-run scheme.

The treatment of Case A is straightforward as we have the complete information about $\mathcal{R}_{\lambda}\approx\mathcal{E}_{\lambda}$ in the large-$N$ limit. If we denote $\sigma_{>}$ and $\sigma_{<}$ to respectively be the largest and smallest eigenvalue of $\widetilde{\boldsymbol{F}_{\textsc{ml}}}^{-1/2}$, then the corresponding outer and inner radii are $R_{\text{out}}=\sigma_{>}$ and $R_{\text{in}}=\sigma_{<}$ [see Fig. 3(a)], so that $\mathrm{cmpl}_{\text{A}}=O\!\left(D^{7}\mathrm{cond}\!\left(\boldsymbol{F}_{\textsc{ml}}^{-1}\right)\right)$ involves the conditional number $\mathrm{cond}\!\left(\boldsymbol{F}_{\textsc{ml}}^{-1}\right)=\sigma_{>}^{2}/\sigma_{<}^{2}$.

The analysis for Case B requires extra care given the complicated state-space boundary $\partial\mathcal{R}_{0}$. While complete and precise details of $\mathcal{R}_{0}$ are absent so far, from Bengtsson et al. (2013), we know that in the Euclidean space, the largest inner sphere inscribable in $\mathcal{R}_{0}$ has a radius that approaches $1/D$ for $D\gg 2$, and that the smallest outer sphere that contains $\mathcal{R}_{0}$ has a radius going to 1 in the same dimension limit. The overall shape of $\mathcal{R}_{0}$ is therefore a “squashed” convex body for large $D$, such that at least one of its dimensions drops appreciably to zero. To estimate the complexity for Case B, we consider CRs of two tractable types: a Type I CR is located at an extremely sharp corner of $\mathcal{R}_{0}$ that is made from at least one of its rapidly shrinking dimensions, as shown in Fig. 3(b), whereas a Type II CR is situated at an extremely flat boundary of $\mathcal{R}_{0}$ where all of its dimensions remain approximately constant within the CR as in Fig. 3(c). For a conservative estimate of $\mathrm{cmpl}$, we consider an $\mathcal{R}$ such that the longest axis of $\mathcal{E}^{\prime}_{\lambda}$ is aligned with the flat surface. All other types of Case-B CRs may be viewed as intermediate situations of these two and have no analytical complexity estimates known to us. The data-copy number $N\gg 1$ is assumed to be sufficiently large such that $\bm{g}_{\textsc{ml}}\approx\boldsymbol{0}$ and $\boldsymbol{r}_{\text{c}}\approx\widehat{\boldsymbol{r}}_{\textsc{ml}}$.

To estimate $\mathrm{cmpl}$ for a Type I CR, we assume that the corner is extremely sharp in one particular dimension such that the curvature of $\partial\mathcal{R}_{0}$ extending out from $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is almost flat. Then following Fig. 4(a), the concept of similar figures give $R_{\text{out}}/R_{\text{in}}\approx D$, which is independent of $\boldsymbol{F}_{\textsc{ml}}$ for extremely sharp corners, and $\mathrm{cmpl}_{\text{B,I}}=O\!\left(D^{9}\right)$. The complexity for Type II CRs may be estimated with the help of Fig. 4(b), where $R_{\text{out}}/R_{\text{in}}\approx 2\,\mathrm{cond}\!\left(\boldsymbol{F}_{\textsc{ml}}^{-1/2}\right)$ is now independent of $\partial\mathcal{R}_{0}$ due to its extremely mild edge features, leading us to $\mathrm{cmpl}_{\text{B,II}}=O\!\left(D^{7}\mathrm{cond}\!\left(\boldsymbol{F}_{\textsc{ml}}^{-1}\right)\right)=\mathrm{cmpl}_{\text{A}}$.

V.3 Other numerical methods

Other numerical methods apart from hit-and-run may also be used to perform in-region sampling, each of which has its own merits and shortcomings Shang et al. (2015); Del Moral et al. (2006); Metropolis et al. (1953); Hastings (1970). Classical rejection and importance sampling methods are two straightforward ways to acquire samples distributed according to some desired prior distribution. For large $D$ and $N$, these methods rapidly become infeasible due to the decreasing ratio of the CR volume to the full Hermitian sampling volume that includes many more unphysical operators that are not quantum states. To cope with this low-yield difficulty, another Markov-chain method known as the Metropolis-Hastings MC scheme may also be used to do in-region sampling. Such a scheme also suffers from high correlations that are generally dependent on the starting point of a Markov-chain iteration. Hamiltonian MC methods are yet another promising class of algorithms that permit larger sample-point hopping that gives a sample with weak correlations. The scalability of such methods are, however, still work in progress.

VI Distance-induced region capacity

VI.1 The operational definition

The theory of in-region sampling seamlessly paves the way to other creative ways of defining the capacity of a region. Doing so permits us to talk about “how big” a CR is without referring to $\mathcal{R}_{0}$ entirely. To begin, one could measure the region capacity in terms of the average distance between any two points inside $\mathcal{R}$, which is a separate idea from the prior content. Intuitively, the smaller this average distance, the smaller the region and vice versa. Using this simple prescription, we propose the region-average quantity

[TABLE]

to measure the capacity of $\mathcal{R}_{\lambda}$, where $\mathcal{D}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})$ is some pre-chosen distance metric. Notice that the ML estimator $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is selected to be the reference point from which distances are measured without loss of generality.

To concretize all results, we shall look at three distance measures for states that enjoy a good reputation in quantum-information studies. We first mention the Hilbert-Schmidt (HS) distance

[TABLE]

which is equivalent to the squared $l_{2}$-norm of $\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}}$. Closely related to the HS distance is the trace-class distance

[TABLE]

defined by the operator absolute value $|A|=\sqrt{A^{\dagger}A}$. To introduce the third measure, we start by quoting the expression of quantum fidelity Uhlmann (1976)

[TABLE]

between $\rho^{\prime}$ and $\widehat{\rho}_{\textsc{ml}}$, to which we can define the Bures distance Bures (1969); Braunstein and Caves (1994) $\mathcal{D}_{\textsc{b}}=2(1-\sqrt{\mathcal{F}})$. In the limit of large $N$, where the fidelity $\mathcal{F}\approx 1-\epsilon$ differs from 1 by a small amount, $\mathcal{D}_{\textsc{b}}$ is also approximately the infidelity $1-\mathcal{F}$.

VI.2 Monotonic behavior of $S_{\mathcal{D},\lambda}$ for $N\gg 1$

Here, we show that, at least for sufficiently large $N$, $S_{\mathcal{D},\lambda}$, defined by any of these three distance measures, behave correctly as a capacity function in the sense that $S_{\mathcal{D},\lambda}$ should not increase as $\lambda$ increases. We first look at the more complicated Case B, and argue that since the $\mathcal{R}_{0}$ generally has only corners and edges with no other mathematically pathological features, a set of hyperplanes can then be used to model any particular boundary feature on which $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ resides (see Fig. 5). This results in the asymptotic form

[TABLE]

after the substitution $\boldsymbol{r}^{\prime\prime}=\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}}$.

At this stage, we shall consider the asymptotic expressions of the distance measures. The HS distance $\mathcal{D}_{\textsc{hs}}$ takes on the simplest (quadratic) form out of all three, which very straightforwardly gives the asymptotic dependence $S_{\textsc{hs},\lambda}\rightarrow-\log\lambda$ provided the sufficient condition $(\mathrm{d}\,\alpha\,\boldsymbol{r}^{\prime})=g(\alpha)\,(\mathrm{d}\,\boldsymbol{r}^{\prime})$, which includes the uniform primitive prior $(\mathrm{d}\,\boldsymbol{r}^{\prime})=(\mathrm{d}\,\boldsymbol{r}^{\prime})_{\text{unif}}\equiv\prod_{j}\mathrm{d}\,r_{j}$. It is not difficult to see that the same $\lambda$ dependence applies to Case A by taking $\boldsymbol{w}_{j}=\boldsymbol{0}$, so that $S_{\textsc{hs},\lambda}$ is monotonically decreasing with increasing $\lambda$. Next, according to Appendix A, in the limit of large $D$, $S_{\text{tr}}\sim\sqrt{S_{\textsc{hs}}}$, which is also clearly monotonic as well owing to $S_{\textsc{hs}}$’s monotonicity. For $S_{\textsc{b}}$ and $S_{\mathcal{F}}$, one can perform a Taylor expansion on them about $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ (see Sec. VII) and realizes that both functions asymptotically depend on the dyadic $(\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}})(\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}})^{\textsc{t}}$, so that both region-average distances are also asymptotically decreasing with $\lambda$.

VII Approximation formulas for $S_{\mathcal{D},\lambda}$ and $u_{\lambda}$

The prior content $S_{\lambda}$ discussed alongside $C_{\lambda}$ in Secs. II–V quantifies the size of $\mathcal{R}_{\lambda}$ relative to $\mathcal{R}_{0}$. In our earlier article Teo et al. (2018), analytical approximation formulas for $S_{\lambda}$ were proposed in the large-$N$ limit, all of which are scaled with the volume $V_{\mathcal{R}_{0}}$ of $\mathcal{R}_{0}$. As is also shown later in the section, this volume dependence is associated with the extension of every $\mathcal{R}_{0}$ integral

[TABLE]

to the entire $\boldsymbol{r}^{\prime}$ space ascribed with the uniform primitive prior, which is a reasonable step to obtain analytical results under the central limit theorem since $L$ is narrow enough to reside within the confines of $\mathcal{R}_{\lambda}$ under this limit. Therefore, the valid usage of these theoretical expressions hinges on the availability of $V_{\mathcal{R}_{0}}$. In quantum-state tomography where we have no complete theoretical information about $\mathcal{R}_{0}$, $V_{\mathcal{R}_{0}}$ is known only for certain priors and state parametrizations  Życzkowski and Sommers (2003); Andai (2006); Bengtsson and Życzkowski (2006); Bengtsson et al. (2013).

On the other hand, it is obvious that $V_{\mathcal{R}_{0}}$ is canceled out for any region-average quantity after such integral extensions. This allows one to derive operational asymptotic formulas for averages like $S_{\mathcal{D},\lambda}$ and $u_{\lambda}$ regardless of $\mathcal{R}_{0}$ in whichever parametrization. As a calculable standard in this section, we continue to derive expressions in terms of the uniform primitive prior and $\boldsymbol{r}$, although the subsequent instructions may also work for other manageable priors with which $S_{\mathcal{D},\lambda}$ behaves as a proper region-capacity function. We first address the different $\mathcal{D}$ measures in the large-$N$ limit.

VII.1 The various $\mathcal{D}$ measures

VII.1.1 Hilbert-Schmidt and trace-class measures

The HS distance measure $\mathcal{D}_{\textsc{hs}}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})$ takes the very simple quadratic form in (12) under any circumstance, whereas the trace-class distance $\mathcal{D}_{\text{tr}}$ has no easy functional form in terms of $\boldsymbol{r}^{\prime}$ for $D>2$. Nevertheless in the limits $N\gg 1$ and $D\gg 2$, based on the principles of random matrix theory detailed in Appendix A, it is deduced that the asymptotic expression

[TABLE]

relating the final $\mathcal{R}$-averages $S_{\textsc{hs}}$ and $S_{\text{tr}}$ is approximately valid for both Case A and B, which incidentally takes the same form found in Zhu and Englert (2011) that was calculated for statistical-fluctuation studies.

VII.1.2 Bures measure

The Bures distance measure $\mathcal{D}_{\textsc{b}}$ also has no tractable functional form in $\boldsymbol{r}^{\prime}$ for general $D$. To find the asymptotic link with $\boldsymbol{r}^{\prime}$ this time, it is technically more convenient to inspect the behavior of $\mathcal{F}$ around $\widehat{\rho}_{\textsc{ml}}\leftrightarrow\widehat{\boldsymbol{r}}_{\textsc{ml}}$ as $N\gg 1$.

A Taylor expansion about $\widehat{\rho}_{\textsc{ml}}$ as guided in Appendix B, we have

[TABLE]

for Case A and

[TABLE]

for Case B, where $P_{r}$ is the projector onto the support of $\widehat{\rho}_{\textsc{ml}}$ having the rank-deficient spectral decomposition $\widehat{\rho}_{\textsc{ml}}=\sum^{r}_{j=1}\left|{\lambda_{j}}\right>\lambda_{j}\left<{\lambda_{j}}\right|$, and

[TABLE]

VII.2 Case A: hyperellipsoidal theory

The presentation in Sec. VII.1 reduces the necessary ingredients for large-$N$ (or $D$) analytical estimations of $S_{\mathcal{D},\lambda}$ to just the scalar $\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})$, column $\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,\boldsymbol{\Delta}^{\prime}_{\textsc{ml}}$ and dyadic $\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,\boldsymbol{\Delta}^{\prime}_{\textsc{ml}}\,{\boldsymbol{\Delta}^{\prime}}^{\mathrm{\,\textsc{t}}}_{\textsc{ml}}$, where $\boldsymbol{\Delta}^{\prime}_{\textsc{ml}}=\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}}$.

When $\mathcal{R}_{\lambda}\approx\mathcal{E}_{\lambda}$, these three integrals takes on simple analytical forms. We start with

[TABLE]

and transform $\boldsymbol{r}^{\prime}\rightarrow\boldsymbol{r}^{\prime\prime}=\boldsymbol{D}^{1/2}\,{\boldsymbol{O}}^{\mathrm{\,\textsc{t}}}\boldsymbol{\Delta}^{\prime}_{\textsc{ml}}$ to the translated diagonal coordinate variables of $\boldsymbol{F}_{\textsc{ml}}/(-2\,\log\lambda)=\boldsymbol{O}\,\boldsymbol{D}\,{\boldsymbol{O}}^{\mathrm{\,\textsc{t}}}$, so that in the large-$N$ limit and uniform primitive prior, we may relax the boundary of $\mathcal{R}_{0}$ and write

[TABLE]

which is a function of the volume $V_{d}=\pi^{d/2}/(d/2)!$ of the $d$-dimensional unit hyperball, the inverse of $\boldsymbol{F}_{\textsc{ml}}$ that characterizes $\mathcal{E}_{\lambda}$ together with the logarithm of $\lambda$.

In this case, the integral column is zero since the integrand after variable transformation becomes odd in $\boldsymbol{r}^{\prime\prime}$, and we are thus left with

[TABLE]

and

[TABLE]

where the last equality is explained by the orthogonally invariant of the $(d-1)$-dimensional solid-angle measure over the unit columns $\boldsymbol{e}$, and so

[TABLE]

With all these components, the relevant asymptotic formulas concerning all three distance measures

[TABLE]

Here $\mathrm{Tr}$ now addresses the dyadic character, as opposed to $\mathrm{tr}$, and we witness the manifestation of logarithmic divergences from both the relaxation of $\partial\mathcal{R}_{0}$ and Gaussian approximation of $L$.

Next, to analytically calculate $u_{\lambda}$ using $f(L)=\log L$ with which $C_{\lambda}$ can be found, we note that due to the Gaussian form of $L$,

[TABLE]

is a dyadic trace function of $\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,\boldsymbol{\Delta}^{\prime}_{\textsc{ml}}\,{\boldsymbol{\Delta}^{\prime}}^{\mathrm{\,\textsc{t}}}_{\textsc{ml}}$, so that we may use the right-hand side of (25) and put down

[TABLE]

after some basic trace and logarithmic manipulations. It is clear that $d/(d+2)\leq u_{\text{A},\lambda}\leq 1$ is bounded.

VII.3 Case B: hyperellipsoidal-cap theory

In Case B, although the geometry of $\mathcal{R}_{\lambda}\approx\mathcal{E}^{\prime}_{\lambda}\cap\mathcal{R}_{0}$ is now much trickier to deal with, the central limit theorem proposed in Sec. III allows us to approximate $\mathcal{R}_{\lambda}$ by a regular analytical region.

As shown in Fig. 6, one can introduce a hyperplane $\mathcal{P}$, described by $\boldsymbol{n}\bm{\cdot}(\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}})=0$ ($\boldsymbol{n}\propto\bm{g}_{\textsc{ml}}$) that is tangent to the level curve of the Gaussian function in (4) at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$. The hyperspherical cap formed by $\mathcal{P}$ and $\mathcal{E}^{\prime}_{\lambda}$ hence asymptotically contains $\mathcal{R}_{\lambda}$, where we have essentially modeled the highly nontrivial $\partial\mathcal{R}_{\lambda}\cap\partial\mathcal{R}_{0}$ as $\mathcal{P}$. This model implies the estimated assignment

[TABLE]

The change of variable $\boldsymbol{r}^{\prime}\rightarrow\boldsymbol{r}^{\prime\prime}=\boldsymbol{D}^{\prime 1/2}\,{\boldsymbol{O}^{\prime}}^{\mathrm{\,\textsc{t}}}\,(\boldsymbol{r}^{\prime}-\boldsymbol{r}_{\text{c}})$ with respect to the diagonal coordinates of $\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda^{\prime})=\boldsymbol{O}^{\prime}\,\boldsymbol{D}^{\prime}\,{\boldsymbol{O}^{\prime}}^{\mathrm{\,\textsc{t}}}$ leads to

[TABLE]

for any function $q$, which is parametrized by the cap element $(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{cap}}=(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{unif}}\,\eta(1-\boldsymbol{r}^{\prime\prime 2})\,\eta(a-{\boldsymbol{b}}^{\mathrm{\,\textsc{t}}}\,\boldsymbol{r}^{\prime\prime})$, $a={\boldsymbol{n}}^{\mathrm{\,\textsc{t}}}(\widehat{\boldsymbol{r}}_{\textsc{ml}}-\boldsymbol{r}_{\text{c}})$ and $\boldsymbol{b}=\boldsymbol{D}^{\prime-1/2}\,{\boldsymbol{O}^{\prime}}^{\mathrm{\,\textsc{t}}}\,\boldsymbol{n}$. One can check that

[TABLE]

In other words, we have

[TABLE]

and that for any $q_{\lambda}$ belonging to either one of the three distance measures or $\log L-\log(\lambda L_{\text{max}})$, as reasoned in Sec. VII.3, the building blocks of $\overline{q_{\lambda}(\boldsymbol{r}^{\prime})}^{\mathcal{R}_{\lambda}}$ are only $\int(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{cap}}$, $\int(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{cap}}\,\boldsymbol{r}^{\prime\prime}$ and $\int(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{cap}}\,\boldsymbol{r}^{\prime\prime}\,{\boldsymbol{r}^{\prime\prime}}^{\mathrm{\,\textsc{t}}}\vphantom{{M^{M}}^{M}}$. These integrations are all carried out in Appendix C.

In combining all results gathered from Appendices B and C, we denote $\mathcal{N}_{d,l,x}=V_{d}\,\mathrm{I}_{(1-l)/2}((d+x)/2,(d+x)/2)$, which depends on the incomplete Euler’s beta function $\mathrm{I}_{\cdot}(\cdot,\cdot)$, and organize two new auxiliary quantities

[TABLE]

This helps to clean the respective formulas

[TABLE]

for the distance-induced capacity functions and

[TABLE]

We caution the Reader once more regarding the actions of $\mathrm{tr}$ and $\mathrm{Tr}$ at the right-hand side of $S^{\text{ (B)}}_{\textsc{b},\lambda}$ in (34).

For consistency, we end this section by noting that Eqs. (34) and (35) cover Eqs. (26) and (27) because Case A implies that $\lambda^{\prime}=\lambda$ ($\bm{g}_{\textsc{ml}}=\boldsymbol{0}=\boldsymbol{m}$), such that $l=0$ then gives $\mathcal{N}_{d,0,x}=V_{d}$ and $\boldsymbol{M}=(-\log\lambda)\,\boldsymbol{F}_{\textsc{ml}}^{-1}/(d+2)$.

VIII Results and discussions

VIII.1 Region reconstruction

We first present, under the uniform primitive prior $p(\boldsymbol{r})\propto 1$, the computation results of $S_{\lambda}$ and $C_{\lambda}$ from $u_{\lambda}$ in Figs. 7 and 8 for quantum systems of various dimensions $D$. To be more technically precise about our use of Euler’s method described in Sec. IV, we first solve (9) for $y_{\lambda}$ by iterating (10) starting with a numerically small $\lambda$ value to $\lambda\approx 1$ using the function $f(\lambda)=\log L$.

The behavior of $S_{\lambda}$ shows the expected decreasing trend not only in $\lambda$, but also in overall magnitude as $D$ increases. This indicates that the (log-)likelihood is turning into a delta-function peak. For larger $D$ or $N$, the computational accuracy of $S_{\lambda}$ and $C_{\lambda}$ using Euler’s numerical method may be maintained by exploring many more $\lambda$ values near zero, such as via logarithmic scaling of $\lambda$, as all curves possess sharp gradient changes in this $\lambda$ range.

In Fig. 9, both simulated data and theoretical curves of all three capacity functions $S_{\textsc{hs}}$, $S_{\text{tr}}$ and $S_{\textsc{b}}$ are plotted against the credibility $C$ for Case A. In this case, there exists no other factors that could spoil the perfect hyperellipsoidal geometry of $\mathcal{R}_{\lambda}$. As such, the analytical curves fit almost perfectly with the simulated points. We note that even the average trace-class distance $S_{\text{tr}}$, which is approximated with (17) through the theory of random matrices, performs very well relative to the simulated data points.

In Case B, we can start to see discrepancies between theory and simulation from Fig. 10 especially for larger $D$. Such deviations are inevitable as the hyperellipsoidal-cap estimation of the actual CR $\mathcal{R}_{\lambda}$ proposed in Sec. VII.3 introduces additional space outside $\mathcal{R}_{0}$ that is certainly not contained in $\mathcal{R}_{\lambda}$. More specifically, for very large $D$, if the rank-deficient ML estimator $\widehat{\rho}_{\textsc{ml}}$ is located at an extremely sharp state-space corner, which is labeled as the Type-I situation in Sec. V.2, this additional space would be exceedingly large relative to the physical CR, which incurs a proportionately large theory-simulation mismatch. On the other hand, if the rank-deficient $\widehat{\rho}_{\textsc{ml}}$ lies on a relatively flat part of the state-space boundary (the Type-II scenario), then this overestimated space, and hence the mismatch, would be much smaller. ML estimators of ranks 1 and 3, which are considered in Fig. 10, are prime examples of the respective Type-I and II situations. Regardless, the asymptotic formulas in Sec. VII.3 may still be used for an order-of-magnitude estimation of $S_{\mathcal{D}}$ and $C$.

We note that in-region sampling is not restricted to just the uniform distribution, so long as the MC method employed is sufficiently general. Such is the case for hit-and-run. For a calibration check of the general hit-and-run algorithm in Sec. V, we generate and compare both uniform and Gaussian distributions with their respective theoretically derived counterparts for the single-qubit case in Fig. 11.

Next, as a real demonstration, we consider another natural prior that is asymptotically conjugate to the likelihood function, that is the Gaussian form $p(\boldsymbol{r})\propto\exp(-(\boldsymbol{r}-\widehat{\boldsymbol{r}}_{\textsc{ml}})\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}(\boldsymbol{r}-\widehat{\boldsymbol{r}}_{\textsc{ml}})/(2N))$ having a much broader covariance $N\boldsymbol{F}^{-1}_{\textsc{ml}}$ defined by the Fisher information for one copy. This prior is a logical choice given that our knowledge about $\boldsymbol{r}$ is updated to $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ after the measurement, which should be used as the most recent prior information for future Bayesian analyses. The corresponding marginal distribution needed to sample along line segments in hit-and-run is therefore the one-dimensional Gaussian distribution of mean $[\boldsymbol{\mathrm{e}_{v}}\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}(\widehat{\boldsymbol{r}}_{\textsc{ml}}-\boldsymbol{r}_{\text{ref}})]/\boldsymbol{\mathrm{e}_{v}}\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}\boldsymbol{\mathrm{e}_{v}}$ and variance $1/\boldsymbol{\mathrm{e}_{v}}\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}\boldsymbol{\mathrm{e}_{v}}$.

VIII.2 Correlation properties of hit-and-run

We recall that, at least under the uniform primitive prior $(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{unif}}$, hit-and-run converges efficiently to the correct uniform distribution in $O(d^{2}R^{2}_{\text{out}}/R^{2}_{\text{in}})$ as discussed in Sec. V.2. Furthermore, it was argued Lovász and Vempala (2006), as a consequence of the above expression, that given an initial point that has a shortest distance $l$ from the boundary $\partial\mathcal{R}_{\lambda}$, hit-and-run eventually mixes sample points into the uniform distribution after $O(d^{3}R^{2}_{\text{out}}/R^{2}_{\text{in}}\log(R_{\text{out}}/l))$ steps.

This reveals a technical caveat for almost all Markov-chain random-walk algorithms: a random walk starting from a sharp corner point of a convex body requires very many steps to approach the stationary target distribution. It is generally much harder to scout the entire convex region from such a corner than from an interior point, since the Markov chain terminates as soon as a next admissible point is obtained, which is probabilistically near the initial corner point around which the admissible region is tight. Such a situation is essentially status quo for high-dimensional state reconstruction where the state space $\mathcal{R}_{0}$ is filled with plenty of extremely sharp corners.

Doing hit-and-run from an interior point is therefore a primary objective for general CR analysis. Even without the full knowledge about the CR, it is still possible to numerically compute a point that is sufficiently interior for this purpose. The idea is to first find multiple random states on $\partial\mathcal{R}_{\lambda}\cap\mathrm{int}(\mathcal{R}_{0})$, the boundary of $\mathcal{R}_{\lambda}$ in the interior of $\mathcal{R}_{0}$, and next average all these states to obtain an interior state of $\mathcal{R}_{\lambda}$. This is evidently equivalent to the minimization of the convex function $[(\boldsymbol{x}^{\prime}-\boldsymbol{r}_{\text{c}})\bm{\cdot}\boldsymbol{F}_{\textsc{ml}}\bm{\cdot}(\boldsymbol{x}^{\prime}-\boldsymbol{r}_{\text{c}})/(-2\log\lambda^{\prime})-1]^{2}$ with respect to $\rho^{\prime}\geq 0$ for which $\boldsymbol{x}^{\prime}=\mathrm{tr}\!\left\{\rho^{\prime}\boldsymbol{\Omega}\right\}$ multiple times in the large-$N$ limit. Fortunately, this can be carried out extremely quickly by using the superfast accelerated projected gradient routine Shang et al. (2017) (see Appendix E).

Figures 14 and 15 supply graphical visualization of the key sampling activities that goes on for Case B with $D=4$, where a rank-one ML estimator is obtained. It is clear that starting hit-and-run from a corner point (namely the ML estimator, for instance) introduces small average hopping distances for subsequent Markov-chains. This can be interpreted as strong sample correlations that prevent wide coverage of the CR, contrary to performing hit-and-run starting from an interior point.

VIII.3 The constructions of plausible regions

The matter of inspecting $S_{\lambda}$ for a fixed $C_{\lambda}$, say 0.95, is rather subjective, for very often one requires experienced judgment to decide if such a value is sufficient for subsequent prediction tasks. As advocated in Evans (2016); Al-Labadi et al. (2018), there exists a statistically meaningful treatment of the measured dataset $\mathbb{D}$ based on the concept of evidence. It is thus fitting for us to end this article with a short review on how our in-region sampling technique studied here directly supports another interesting kind of Bayesian analysis.

By definition, we say that $\boldsymbol{r}^{\prime}$ is a plausible candidate parameter for the true $\boldsymbol{r}$ if there is indeed evidence in favor of this supposition. That is, its normalized posterior probability $L(\mathbb{D}|\boldsymbol{r}^{\prime})\,p(\boldsymbol{r}^{\prime})/L(\mathbb{D})$ after the measurement is larger than its prior probability $p(\boldsymbol{r}^{\prime})$ before this measurement was performed. Therefore, $\boldsymbol{r}$ is a plausible parameter if the evidence supports the prior knowledge. Under this evidence-belief framework, one can construct another type of Bayesian region—the plausible region (PR)—that contains all plausible choices of $\boldsymbol{r}$. This is really the CR $\mathcal{R}=\mathcal{R}_{\lambda=\lambda_{\rm{crit}}}$ characterized by the critical value Li et al. (2016)

[TABLE]

for which $L(\mathbb{D}|\boldsymbol{r}\in\partial\mathcal{R}_{\lambda=\lambda_{\rm{crit}}})=L(\mathbb{D})$, or the CR that contains all plausible points and nothing else. This follows quickly from the following equality chain:

[TABLE]

so that $L(\mathbb{D}|\boldsymbol{r}\in\partial\mathcal{R}_{\lambda=\lambda_{\rm{crit}}})\equiv\lambda_{\rm{crit}}L_{\text{max}}=L(\mathbb{D})$ gives Eq. (36).

So, constructing a PR is nothing more than one additional step of computing $\lambda_{\text{crit}}$ after a CR construction. In our previous work Teo et al. (2018); Oh et al. (2018), we have supplied MC-less asymptotic approximations to the expression of $\lambda_{\text{crit}}$ for the uniform primitive prior. In the current context, it clearly follows that $\lambda_{\text{crit}}$ is also directly computable by simply doing a Riemann summation of the full $S_{\lambda}$ spectrum obtained through in-region sampling in accordance with Eq. (36).

IX Conclusions

Quantum-state tomography is an important application of multidimensional parameter estimation. The construction of Bayesian credible regions for the reconstructed quantum states after tomography is, unfortunately, a highly nontrivial problem owing to the complex constraints inherited from the state space. Standard numerical recipe of first doing a Monte Carlo sampling of the state space and next discarding points outside the credible region to compute its region qualities (size and credibility) quickly becomes infeasible when the dataset collected in an experiment is relatively big, as the corresponding credible region would be very small with respect to the state space.

In this article, we performed an extensive analysis of our recent in-region sampling technique that can construct credible regions that are usually very small in practice for reasonably high credibility values and large data samples. This technique computes credible-region qualities of a small credible region for any given prior distribution by inspecting how an appropriately chosen region-average quantity changes as the shape of the region varies. This procedure transforms the general credible-region construction into a sequence of direct region sampling followed by simple numerical solution to a single-variable differential equation. This results in no sample wastage since no points are discarded. The method of accelerated hit-and-run is one numerical scheme that can be used to compute region averages rather efficiently provided that sample correlation is properly mitigated with good Monte Carlo starting points. One can also estimate its numerical complexity in the context of tomography despite the complicated state-space boundary.

Furthermore, for highly complex quantum systems of extremely large dimensions, where all numerical methods eventually become practically infeasible, we derive a set of analytical formulas to perform approximate Bayesian error certification through the perspective of distance-induced region capacity measures that alternatively quantifies how large a credible region is. These formulas are now fully operational and further complement those for the conventional size function developed in previous works that require knowledge of the state-space volume.

Acknowledgements.

The authors thank J. Shang for fruitful discussions, and acknowledge financial support from the BK21 Plus Program (21A20131111123) funded by the Ministry of Education (MOE, Korea) and National Research Foundation of Korea (NRF), the framework of international cooperation program managed by the NRF (NRF-2018K2A9A1A06069933), and the Basic Science Research Program through the NRF funded by the Ministry of Education (No. 2018R1D1A1B07048633).

Appendix A The relationships between $S_{\textsc{hs}}$ and $S_{\text{tr}}$ in the large-$N$ (or $D$) limit

Apart from $\mathcal{D}_{\textsc{hs}}$, all other measures have no direct analogs in the $\boldsymbol{r}^{\prime}$ parametrization. However in certain limits, all these measures have approximate relations with $\mathcal{D}_{\textsc{hs}}$.

We start with making an approximate connection between $S_{\text{tr}}$ and $S_{\textsc{hs}}$ by examining the Hermitian operator $\Delta\rho^{\prime}=\rho^{\prime}-\widehat{\rho}_{\textsc{ml}}$ ($\rho^{\prime}\in\mathcal{R}$). In Case A, the distribution of $\Delta\rho^{\prime}$ in $\mathcal{R}$ has zero mean, $\overline{\Delta\rho^{\prime}}^{\mathcal{R}}=0$. This is also approximately true for the Case B situation when $N$ is sufficiently large such that $\mathcal{R}$ is small. Furthermore, the space of $\Delta\rho^{\prime}$ is essentially a bounded set of Hermitian random operators. Here, we shall make the assumption that each matrix entry $\Delta\rho^{\prime}_{jk}$ in the computational basis is an independent random complex number. Under this condition, the $\Delta\rho^{\prime}$s form what is now known as a Wigner ensemble Wigner (1955, 1957); Pizzo et al. (2013); Mehta (2004) with the second moment equal to $\overline{|\Delta\rho^{\prime}_{jk}|^{2}}^{\mathcal{R}}=\mathrm{tr}\!\left\{(\Delta\rho^{\prime})^{2}\right\}=S_{\textsc{hs}}$. Moreover, they are known to have an i.i.d. eigenvalue spectrum that follows the Wigner semicircle law

[TABLE]

in the large-$D$ limit. The trace-class distance $\mathcal{D}_{\text{tr}}$ can thus be calculated with the integral

[TABLE]

so that we end up with (17).

For Case B, that $\overline{\Delta\rho^{\prime}}^{\mathcal{R}}=0$ is obvious, but as we have no means of analytically estimate $\overline{\Delta\rho^{\prime}}^{\mathcal{R}}$, we make a further approximation that as long as $\mathcal{R}$ is sufficiently small, the offset to $\overline{\Delta\rho^{\prime}}^{\mathcal{R}}$ will proportionately be small, so that (17) remains a reasonable asymptotic approximation.

Appendix B Fidelity in the large-$N$ limit

A Taylor expansion of $\mathcal{F}$ about $\widehat{\boldsymbol{r}}_{\textsc{ml}}$, or

[TABLE]

reveals the large-$N$ characteristics that is needed for analysis. The structure of (14), however, demands the operator variation of $\sqrt{A}$ for a positive (semidefinite) $A$. An integral representation of $\sqrt{A}$ exists Bini et al. (2005) and can be written as

[TABLE]

where the limit is understood to be applied at the very end of all calculations so that Eq. (41) is valid even for $A$ with zero eigenvalues.

The first-order variation of $\mathrm{tr}\!\left\{\sqrt{A}\right\}^{2}$ produces

[TABLE]

In terms of $\mathcal{F}$, we substitute $A\equiv\widehat{\rho}_{\textsc{ml}}^{1/2}\,\rho^{\prime}\,\widehat{\rho}_{\textsc{ml}}^{1/2}$, and evaluate the above result with $\rho^{\prime}=\widehat{\rho}_{\textsc{ml}}$, or $A\rightarrow A_{\textsc{ml}}=\widehat{\rho}_{\textsc{ml}}^{2}$, then with $\updelta A_{\textsc{ml}}=\widehat{\rho}_{\textsc{ml}}^{1/2}\,\updelta\widehat{\boldsymbol{r}}_{\textsc{ml}}\bm{\cdot}\boldsymbol{\Omega}\,\widehat{\rho}_{\textsc{ml}}^{1/2}$,

[TABLE]

where we remind the Reader that $\mathrm{tr}$ acts on operators only, not on the vectorial character. For Case B in which $\widehat{\rho}_{\textsc{ml}}=\sum^{r}_{j=1}\left|{\lambda_{j}}\right>\lambda_{j}\left<{\lambda_{j}}\right|$ is rank-deficient, we get, after taking the trace,

[TABLE]

where $P_{r}=\sum^{r}_{j=1}\left|{\lambda_{j}}\right>\left<{\lambda_{j}}\right|$. It is then trivial to realize that this first-order derivative is zero for Case A. Qualitatively, this confirms the fact that when $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is an interior point, $\mathcal{F}$ has a local maximum at this point as it should, while a boundary estimator evaluates to a nonzero $\mathcal{F}$ slope.

Upon denoting $\boldsymbol{W}_{\textsc{ml}}=\widehat{\rho}_{\textsc{ml}}^{1/2}\,\boldsymbol{\Omega}\,\widehat{\rho}_{\textsc{ml}}^{1/2}$, the second-order variation follows from the second line of (42):

[TABLE]

A product-rule dissociation of (45) comprises a $\updelta\,\mathrm{tr}\!\left\{\sqrt{A}\right\}$ and

[TABLE]

After evaluating the variation at $\rho^{\prime}=\widehat{\rho}_{\textsc{ml}}$ and further undoing all integrations with the help of its spectral decomposition, Case B yields

[TABLE]

The counterpart expression for Case A is immediate, of course.

Appendix C Hyperellipsoidal-cap averages

Under the uniform primitive prior, calculations of the hyperellipsoidal-cap integrals

[TABLE]

specified by the uniform cap-volume element $(\mathrm{d}\,\boldsymbol{r}^{\prime\prime})_{\text{cap}}=(\mathrm{d}\,\boldsymbol{r}^{\prime})_{\text{unif}}\,\eta(1-\boldsymbol{r}^{\prime\prime 2})\,\eta(a-{\boldsymbol{b}}^{\mathrm{\,\textsc{t}}}\boldsymbol{r}^{\prime\prime})$ for $0\leq a\leq|\boldsymbol{b}|$ and some column $\boldsymbol{b}$, include systematic manipulations of the double Heaviside functions. One route to take exploits the following integral representation

[TABLE]

with the implicit limit $\epsilon\rightarrow 0^{+}$. We then have, for (48),

[TABLE]

upon noting the well-known $d$-dimensional Gaussian integral result

[TABLE]

for any positive $\boldsymbol{A}$. Let us first do the $t^{\prime}$ integration by invoking the useful transformation

[TABLE]

[TABLE]

As a consequence,

[TABLE]

The above integral in $y$ represents well-known special functions and to see this, we further perform the substitutions $\cos u=(y-a)/|\boldsymbol{b}|$ and $l=a/b$:

[TABLE]

which is a product of the beta function and its normalized incomplete form

[TABLE]

The final answer reads

[TABLE]

For (49),

[TABLE]

where the $\boldsymbol{r}^{\prime\prime}$ integration

[TABLE]

is simplified after a differentiation under the integral sign. Then

[TABLE]

To simplify the $t$ integration, we again recall Eq. (54) to finally get

[TABLE]

We can at least verify the $d=1$ for Eq. (63) after paying attention to the convention $\boldsymbol{b}\rightarrow-b$, for $b\geq 0$. This corresponds to the integral

[TABLE]

By the same token, we may explore the dyadic integral $\boldsymbol{I}_{2}$ in (48) first with (51) to obtain

[TABLE]

where the dyadic $\boldsymbol{r}^{\prime\prime}$ sub-integral

[TABLE]

after an application of Eq. (53) and a dyadic differentiation under the integral sign this time.

This time, we choose to perform the $t$ and $t^{\prime}$ integrals before taking the derivative, inasmuch as

[TABLE]

where the usage of (54) evaluates the $t^{\prime}$ integral

[TABLE]

into another Gaussian integral. Its convenient feature becomes clear when substituted back into (67):

[TABLE]

The end of the tunnel becomes visible after a product-rule differentiation carried out with the basic dyadic identities

[TABLE]

after which we end up with the final answer

[TABLE]

The 1D special case can again be extracted from Eq. (71),

[TABLE]

and compared with the direct calculation

[TABLE]

Appendix D Projection of uniform and Gaussian distributions onto a hyperellipsoid

One may begin with a rotated coordinate system (centered at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$) that diagonalizes the Fisher information $\boldsymbol{F}_{\textsc{ml}}$, so that the projected uniform distribution onto the error region in the large-$N$ approximation is calculated from the integral

[TABLE]

for the eigenvalues $a$, $b$, and $c$ of $\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda)$. The preceding exercises of Appendix C swiftly gives $p_{\text{unif}}\propto\sqrt{1-ax^{2}-by^{2}}$.

The Gaussian distribution, with covariance chosen to be proportional to $\boldsymbol{F}_{\textsc{ml}}$ that possesses the eigenvalues $a^{\prime}$, $b^{\prime}$ and $c^{\prime}$, is given by

[TABLE]

This can be simplified to p_{\text{gauss}}\propto\mathrm{e}^{\mbox{\footnotesize-a^{\prime}x^{2}-b^{\prime}y^{2}}}\,\gamma(1/2,c^{\prime}(1-ax^{2}-by^{2})/c) in terms of the lower incomplete Gamma function $\gamma(\cdot,\cdot)$ using again results from Appendix C.

Appendix E Gradient optimization for obtaining an error-region interior point

To acquire an interior point of a CR for state tomography, which is essentially a Hilbert subspace, it is sufficient to generate very many ($>D$) region boundary points and take the average of these points. In the limit of large $N$, we may approximate the inner boundary of the region as part of the hyperellipsoid described by $(\boldsymbol{x}^{\prime}-\boldsymbol{r}_{\text{c}})\bm{\cdot}\boldsymbol{A}\bm{\cdot}(\boldsymbol{x}^{\prime}-\boldsymbol{r}_{\text{c}})\leq 1$.

If we perform a variation on the relevant function

[TABLE]

where $\boldsymbol{A}=\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda^{\prime})$, we get

[TABLE]

A gradient method, such as the accelerated projected gradient method, requires the definition of the (operator) gradient defined by $\updelta f/\updelta\rho^{\prime}$, which requires the connection $\updelta\,\boldsymbol{x}^{\prime}=\mathrm{tr}\!\left\{\updelta\rho^{\prime}\,\boldsymbol{\Omega}\right\}$. Naturally then, we must acquire the resulting operator

[TABLE]

where the dot products operate only on the vectorial character only.

The mechanisms that drive the accelerated projected gradient search algorithm are beyond the scope of this article. Instead we provide a simple manual to immediately modify and use the open-source MATLAB code file qse_apg.m that is available on https://github.com/qMLE/qMLE. For this purpose, we note the three important variables fval_varrho, fval_new and gradient, which stores the function values of $f$ evaluated with the varrho and rho_new variables, as well as the gradient operator $\updelta f/\updelta\rho^{\prime}$ evaluated with varrho. In order to minimize $f$ with qse_apg.m, one may simply overwrite the existing functional expressions [namely -f’.*log(probs_…) and -qmt(…)] for the three variables with the ones in Eqs. (76) and (78). By our numerical experience with this minimization task, it is advisable to set the parameters defaults.threshold_fval and defaults.imax respectively to eps and $>10^{-8}$ for better accuracies.