Probing Bayesian credible regions intrinsically: a feasible error certification for physical systems

TL;DR

This paper introduces a new method for efficiently certifying the error and credibility of Bayesian estimators in physical systems, especially quantum state tomography, using average log-likelihoods instead of computationally intensive sampling.

Contribution

It reformulates Bayesian credible-region theory to rely on region averages, enabling feasible error certification without rejection sampling or full parameter space knowledge.

Findings

Region-average log-likelihood suffices for size and credibility estimation.

Analytical formulas for $l_2$-norm distance provide asymptotic error certification.

Method applicable to large datasets and quantum systems without Monte Carlo methods.

Abstract

Computing size and credibility of Bayesian credible regions for certifying the reliability of any point estimator of an unknown parameter (such as a quantum state, channel, phase, \emph{etc.}) relies on rejection sampling from the entire parameter space that is practically infeasible for large datasets. We reformulate the Bayesian credible-region theory to show that both properties can be obtained solely from the average of log-likelihood over the region itself, which is computable with direct region sampling. Neither rejection sampling nor any geometrical knowledge about the whole parameter space is necessary, so that general error certification now becomes feasible. We take this region-average theory to the next level by generalizing size to the average $l_p$-norm distance $(p>0)$ between a random region point and the estimator, and present analytical formulas for $p=2$ to estimate distance-induced size and credibility for any physical system and large datasets, thus implying that asymptotic Bayesian error certification is possible without any Monte~Carlo computation. All results are discussed in the context of quantum-state tomography.

Equations

$$S_{\lambda }\equiv$$

$$C_{\lambda }\equiv$$

$$C_{\lambda}=\left(\lambda S_{\lambda}+\int^{1}_{\lambda}\mathrm{d}\,\lambda^{\prime}S_{\lambda^{\prime}}\right)\bigg{/}\int^{1}_{0}\mathrm{d}\,\lambda^{\prime}S_{\lambda^{\prime}}$$

$$\frac{\partial y_{\lambda }}{\partial \lambda }=-\frac{y_{\lambda }}{\lambda u_{\lambda }}$$

$$S_{\mathcal{D},\lambda}\equiv\overline{\mathcal{D}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})}^{\mathcal{R}_{\lambda}}=\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,\mathcal{D}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})\bigg{/}\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\,\boldsymbol{r}^{\prime})$$

$$S_{{\mathcal{D}}_p,\lambda }\to$$

$$\sim$$

$$S_{2\text{A},\lambda }=\mathrm{Tr}\left\{{\mathbf{F}}_{\text{textsc}ml}^{-1}\right\}\frac{(-\log \lambda )}{d/2+1},u_{\text{A},\lambda }=-\frac{2}{d+2}\log \lambda .$$

$$S_{2\text{B},\lambda }=$$

$$u_{\text{B},\lambda }=$$

$$\times \log (\lambda L_{\text{max}})/\log ({\lambda }^{\prime }{L}_{\text{max}}),$$

$$\mathbf{m}=$$

$$\mathbf{M}=$$

Full text

*Introduction.—*Parameter reconstruction from datasets is a preliminary task in the study of natural sciences. In quantum theory, proper reconstruction of quantum states Smithey et al. (1993); Chuang and Nielsen (2000); Řeháček et al. (2007); Teo et al. (2011a); Zhu (2014), quantum channels O’Brien et al. (2004); Teo et al. (2011b); Fiurášek (2015); Varga et al. (2018), interferometric phases Caves (1981); Dorner et al. (2009), etc., is the root to successful executions of all quantum-information protocols Ladd et al. (2010); Demkowicz-Dobrzański et al. (2015); Campbell et al. (2017); Ladd et al. (2010); Lekitsch et al. (2017). A parameter estimator must be accompanied by an appropriate error certification to ascertain its reliability for future physical predictions. Bootstrapping or resampling Efron and Tibshirani (1993); Davison and Hinkley (1997), which generates mock data from collected ones to obtain “error-bars”, can result in highly overoptimistic “error-bar” lengths Suess et al. (2017) that do not accurately characterize the estimator. From the principles of hypothesis testing, one can instead construct Bayesian credible regions Shang et al. (2013); Li et al. (2016) based on the collected data. These credible regions are distinct from the frequentists’ confidence regions Christandl and Renner (2012); Blume-Kohout (2012); Faist and Renner (2016), which are constructed from the complete (often assumed) distribution of estimators that includes all unobserved ones in the experiment.

A credible region $\mathcal{R}$, which is a Bayesian error region constructed from experimentally observed data $\mathbb{D}$, requires the specification of its size and credibility, which is the probability that the true parameter is inside $\mathcal{R}$. It is well-known from Shang et al. (2013) that the latter is readily derived so long as the functional behavior of the former with the shape of $\mathcal{R}$ is known. As the size of $\mathcal{R}$ is defined as the volume fraction of the full parameter space $\mathcal{R}_{0}$, its computation conventionally requires one to first obtain a large sample of points in $\mathcal{R}_{0}$, and later discard (usually very many) points that are outside $\mathcal{R}$. Acquiring a sufficiently large sample of $\mathcal{R}_{0}$ for a subsequently accurate sample filtering is doable with a number of Monte Carlo (MC) methods Shang et al. (2015); Seah et al. (2015), most notably the Hamiltonian Markov-chain MC, provided that $\mathcal{R}$ is not small. In practice, however, when data sample-size $N$ becomes even moderately large, the region $\mathcal{R}$ (of size $\sim N^{-d/2}$ Teo et al. (2018) for a $d$-dimensional parameter) is too tiny for any MC-filtering sampling to be practically feasible. In Teo et al. (2018); Oh et al. (2018), closed-form approximations are given to estimate both region qualities for large $N$ without MC-filtering, with the premise that the volume of $\mathcal{R}_{0}$ is known.

In this Letter, we develop an in-region sampling theory to compute the size and credibility with neither MC-filtering from nor any geometrical knowledge about $\mathcal{R}_{0}$ (such as its volume). We first prove the central lemma which states that both region qualities are computable from the average of log-likelihood over $\mathcal{R}$. We next discuss the hit-and-run MC algorithm Bélisle et al. (1993); Smith (1996); Lovász and Vempala (2006); Kiatsupaibul et al. (2011) as one of the many numerical tools to perform direct region-average computation. As a strategic bonus, we make use of the region-average concept in in-region sampling to define the region capacity of $\mathcal{R}$ induced by an $l_{p}$-norm ($p>0$) between two points in $\mathcal{R}$. This would allow us to derive fully operational asymptotic approximation formulas for $p=2$ (squared-error metric) to carry out rapid error certifications without numerical computations. All results are demonstrated and verified for multi-qubit tomography.

*Error-region size and credibility.—*For a given informationally complete (IC) dataset $\mathbb{D}$, we would like to reconstruct the unknown $d$-dimensional parameter $\boldsymbol{r}$ (vectorial in general) that fully characterizes some physical system. We shall assume that the parameter space $\mathcal{R}_{0}$ (of quantum states, channels, Cartesian-product of independent quantities, etc.) for the physical system of interest is convex, and take the unique estimator $\widehat{\boldsymbol{r}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}$ to be the maximum-likelihood (ML) estimator Aldrich (1997); Řeháček et al. (2007); Teo (2015), that is the estimator that maximizes the likelihood $L=L(\mathbb{D}|\boldsymbol{r}^{\prime})$. It was formally shown in Shang et al. (2013) that the optimal Bayesian credible region (CR) $\mathcal{R}$ for $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ has an isolikelihood boundary $\partial\mathcal{R}$—a boundary of constant likelihood—and every interior point possessing a likelihood $L\geq\lambda\,L_{\text{max}}$ (see Fig. 1). Its size and credibility are

[TABLE]

where the volume measure $(\mathrm{d}\,{\boldsymbol{r}})$ incorporates some prescribed prior distribution $p(\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})$, and $0\leq\lambda\leq 1$ characterizes the shape and size of $\mathcal{R}_{\lambda}$, so that $\mathcal{R}_{\lambda=0}=\mathcal{R}_{0}$ and $\mathcal{R}_{\lambda=1}=\{\widehat{\boldsymbol{r}}_{\textsc{ml}}\}$. Hence, $S_{\lambda}$ measures the total prior content of $\mathcal{R}_{\lambda}$ that monotonically decreases with increasing $\lambda$, and $C_{\lambda}$ its posterior content that expresses the probability that $\boldsymbol{r}\in\mathcal{R}_{\lambda}$. Both $C_{\lambda=\lambda_{0}}$ (pre-chosen to be 0.95 say) and the corresponding $S_{\lambda=\lambda_{0}}$ are reported together with $\widehat{\boldsymbol{r}}_{\textsc{ml}}$. The relation

[TABLE]

means that a single $\boldsymbol{r}^{\prime}$-integration for $S_{\lambda}$ is sufficient to acquire $C_{\lambda}$ Shang et al. (2013). In realistic experiments, where the desired number of data copies $N<\infty$ is usually large (which we assume unless otherwise stated), the likelihood $L$ becomes a Gaussian function owing to the central limit theorem and peaks strongly around $\widehat{\boldsymbol{r}}_{\textsc{ml}}$. In this case, $S_{\lambda}$ becomes very small even for small $\lambda$ or large $C_{\lambda}$ (the desired situation). Therefore, MC-filtering produces almost no yield as such a finite sample would surely miss $\mathcal{R}_{\lambda}$ for a reasonably high $C_{\lambda}$.

We inform that one systematic guide to report error regions is to invoke the elegant notion of evidence, which leads to the so-called plausible region Li et al. (2016); Evans (2016); Al-Labadi et al. (2018); Teo et al. (2018); Oh et al. (2018) for $\mathbb{D}$, in which all points have posterior probabilities larger than or equal to their prior probabilities—a physical measure of statistical significance. Then $C_{\lambda}$ should not exceed the credibility of this plausible region in order for the CR to contain only plausible points (refer to our companion article Oh et al. (2019) for details).

*In-region sampling theory.—*We shall now propose a way to compute both $S_{\lambda}$ and $C_{\lambda}$ without MC-filtering. The physical intuition behind our theory is to realize that if one inspects the average of some quantity $q_{\lambda}$ over the region $\mathcal{R}_{\lambda}$ [formally denoted by $\overline{q_{\lambda}}^{\mathcal{R}_{\lambda}}=\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\boldsymbol{r}^{\prime})q_{\lambda}(\boldsymbol{r}^{\prime})/\int_{\mathcal{R}_{\lambda}}(\mathrm{d}\boldsymbol{r}^{\prime})$], then its rate of change with $\lambda$ actually encodes information about the behavior of $S_{\lambda}$ with $\lambda$. A shrinkage of $\mathcal{R}_{\lambda}$, for example, translates to an exclusion of some $q_{\lambda}$ values from the region-average. More precisely, this leads to the

Region-average computation (RAC) lemma: For any prior $(\mathrm{d}\,\boldsymbol{r}^{\prime})$ and $N$, the prior content $S_{\lambda}$ (up to a multiplicative factor), and hence the credibility $C_{\lambda}$, are all inferable from the $\mathcal{R}$-average quantity $u_{\lambda}=\overline{\log L(\mathbb{D}|\boldsymbol{r}^{\prime})-\log(\lambda\,L_{\text{max}})}^{\mathcal{R}_{\lambda}}$.

We prove this lemma by taking the first-order derivative of $u_{\lambda}S_{\lambda}$ in $\lambda$. Upon noting that $\partial S_{\lambda}/\partial\lambda=\int_{\mathcal{R}_{0}}(\mathrm{d}\,\boldsymbol{r}^{\prime})\,\delta(L-\lambda\,L_{\text{max}})$, we end up with the following first-order differential equation

[TABLE]

that characterizes the full evolution of $y_{\lambda}=u_{\lambda}\,S_{\lambda}$ given the boundary value $S_{\lambda=0}=1$. Equation (3) can be solved easily by iterating $y_{\lambda_{j+1}}=y_{\lambda_{j}}-y_{\lambda_{j}}/(\lambda_{j}u_{\lambda_{j}})$ following Euler’s method Butcher (2003), so that $C_{\lambda}$ can thereafter be computed using Eq. (2). This closes our constructive proof of the RAC lemma.

For any prior distribution $p(\boldsymbol{r})$, there exist many MC Shang et al. (2015); Del Moral et al. (2006) schemes to compute $u_{\lambda}$, many of which use Markov-chain algorithms. Hit-and-run sampling Bélisle et al. (1993); Smith (1996); Lovász and Vempala (2006); Kiatsupaibul et al. (2011) is one such extensively-studied scheme. The mechanism behind hit-and-run starts with the construction of a simple finite convex set $\mathcal{B}\supseteq\mathcal{R}_{\lambda}$. For $N\gg 1$ and some $\lambda$, two general cases exist as shown in Fig. 1. In Case A, we define $\mathcal{B}$ as the hyperellipsoid $\mathcal{E}_{\lambda}$ centered at $\boldsymbol{r}_{\text{c}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}$ that profiles the Gaussian $L$ whenever $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is an interior point. In Case B, where $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is a boundary point on $\partial\mathcal{R}_{0}$, we set $\mathcal{B}$ as the (truncated) hyperellipsoid $\mathcal{E}^{\prime}_{\lambda}$ centered at $\boldsymbol{r}_{\text{c}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}+\boldsymbol{F}_{\textsc{ml}}^{-1}\,\bm{g}_{\textsc{ml}}$, where $\boldsymbol{F}_{\textsc{ml}}$ is the Fisher information evaluated at $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ and $\bm{g}_{\textsc{ml}}=\partial\log L/\partial\boldsymbol{r}^{\prime}|_{\boldsymbol{r}^{\prime}=\widehat{\boldsymbol{r}}_{\textsc{ml}}}$. Next, starting from a reference point in $\mathcal{R}_{\lambda}$, say the ML estimator $\widehat{\boldsymbol{r}}_{\textsc{ml}}$, a finite line segment, with endpoints on $\partial\mathcal{B}$, passing through this point is generated and a random point is picked repeatedly along this line until it lies in $\mathcal{R}_{\lambda}$, thereafter becoming the next reference through which a new finite line segment is generated to find the next point in $\mathcal{R}_{\lambda}$. The final $\mathcal{R}_{\lambda}$ sample is then used to compute any $\mathcal{R}_{\lambda}$-average quantity. The key point is that a hyperellipsoidal $\mathcal{B}$ for hit-and-run is constructed based on the central limit theorem, where the $N\gg 1$ condition guarantees that the physical region is asymptotically contained in $\mathcal{B}$. To play it safe, a good idea would be to choose a hyperellipsoid that is, say, twice the size of the supposed one given by the theorem.

Beginning with $k=1$ and $\boldsymbol{r}_{\text{ref}}=\widehat{\boldsymbol{r}}_{\textsc{ml}}$ of $N\gg 1$, the accelerated version of hit-and-run Bélisle et al. (1993); Smith (1996); Kiatsupaibul et al. (2011) for any given prior distribution $p(\boldsymbol{r})$ runs as follows: 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}$, $\boldsymbol{A}=\boldsymbol{F}_{\textsc{ml}}/(-2\log\lambda^{\prime})$, $2\log(\lambda/\lambda^{\prime})=\bm{g}_{\textsc{ml}}^{\textsc{t}}\boldsymbol{F}_{\textsc{ml}}^{-1}\bm{g}_{\textsc{ml}}$ [$\bm{g}_{\textsc{ml}}=\boldsymbol{0}$ and $\lambda^{\prime}=\lambda$ for Case A]. 2. Define $\beta_{1}\equiv\mu_{\text{min}}=\min\{\mu_{+},\mu_{-}\}$ and $\beta_{2}\equiv\mu_{\text{max}}=\max\{\mu_{+},\mu_{-}\}$. 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. Determine whether $\boldsymbol{r}_{\text{test}}\in\mathcal{R}_{\lambda}$. If so, define $\boldsymbol{r}_{\text{ref}}=\boldsymbol{r}_{\text{test}}$, raise $k$ by 1, and go to 1. If not, set $\beta_{1}=\beta$ if $\beta<0$ or $\beta_{2}=\beta$ if $\beta>0$, and repeat 3 and 4. Sampling terminates when $k>K_{\text{smp}}$ for a prechosen $K_{\text{smp}}$.

We emphasize that the Gaussian approximation serves only as an efficient guide to contain the sampling space. An additional criterion that $\log L>\log(\lambda L_{\text{max}})$ may be used to further ensure that all sampled points truly lie in $\mathcal{R}_{\lambda}$, although this is almost always the case for $N\gg 1$. One main technical issue for Markov-chain schemes is that the convergence rate is strongly dependent on the starting point (finite sample-point correlation). It is well-known, however, that hit-and-run converges fast to $p(\boldsymbol{r})$ (with essentially polynomial complexity) so long as it starts from any interior point. As an example in 4-qubit tomography, such an interior point can be generated in about 10 seconds per $\lambda$ with $N>2\times 10^{6}$ and 4096 measurement outcomes using accelerated projected gradient method Shang et al. (2017) to minimize the function $[1-(\boldsymbol{x}-\boldsymbol{r}_{\text{c}})\bm{\cdot}\boldsymbol{A}\bm{\cdot}(\boldsymbol{x}-\boldsymbol{r}_{\text{c}})]^{2}$ 32 times (see for instance Lovász and Vempala (2006); Lovász (1999) and our companion article Oh et al. (2019) for more technical discussions).

*Region capacity.—*The region-average methodology used to feasibly compute $S_{\lambda}$ (and $C_{\lambda}$) invites more options to gauge the capacity of $\mathcal{R}$. Instead of measuring prior contents, we may check how close is a randomly-chosen point in $\mathcal{R}$ from $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ on average. Formally, the $\mathcal{R}$-average

[TABLE]

for the capacity of $\mathcal{R}_{\lambda}$ now depends additionally on the metric $\mathcal{D}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})$ one chooses to measure this average distance.

One can argue that if the metric is an $l_{p}$-norm of $p>0$, $S_{\mathcal{D},\lambda}$ monotonically decreases with $\lambda$ when $N\gg 1$ for an appropriate $(\mathrm{d}\,\boldsymbol{r}^{\prime})$. To see this we begin with $\mathcal{D}\equiv\mathcal{D}_{p}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})=\left(\sum_{j}|r^{\prime}_{j}-\widehat{r}_{\textsc{ml},j}|^{p}\right)^{1/p}$. According to Fig. 2, after the substitution $\boldsymbol{r}^{\prime\prime}=\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}}$, we have for the more complicated Case B,

[TABLE]

The same conclusion for Case A follows by definition, and remains unchanged also for $\mathcal{D}_{p}(\boldsymbol{r}^{\prime},\widehat{\boldsymbol{r}}_{\textsc{ml}})=\sum_{j}|r^{\prime}_{j}-\widehat{r}_{\textsc{ml},j}|^{p}$ since $S_{\mathcal{D}_{p},\lambda}\sim(-\log\lambda)^{p/2}$ is also monotonic in $\lambda$. These imply that $S_{\mathcal{D}_{p},\lambda}$ induced by any $l_{p}$-norm behaves as a proper capacity measure in the limit $N\gg 1$ under a sufficient class of priors that includes the uniform primitive prior. The new practice for Bayesian CR certification is then to report the three-tuple $\left(\widehat{\boldsymbol{r}}_{\textsc{ml}},C_{\lambda_{0}}(=0.95\text{ say}),S_{\mathcal{D}_{p},\lambda_{0}}\right)$ for some $p>0$.

*Analytical error certification with region capacity.—*It turns out that the approximated extensions of all $\int_{\mathcal{R}_{0}}$ integrals to the whole $\boldsymbol{r}^{\prime}$ space free all $\mathcal{R}$-average quantities from any geometrical dependence on $\mathcal{R}_{0}$, unlike $S_{\lambda}$ that asymptotically depends on $\mathcal{R}_{0}$’s volume Teo et al. (2018). We may then use this observation to acquire asymptotic formulas for $S_{\mathcal{D}_{p},\lambda}$ and $u_{\lambda}$ to perform approximate analytical error certifications. To this end, we regard $S_{2}\equiv S_{\mathcal{D}_{2}}$ induced by the squared $l_{2}$-norm ($p=2$), $\mathcal{D}\equiv|\boldsymbol{r}^{\prime}-\widehat{\boldsymbol{r}}_{\textsc{ml}}|^{2}$, as the prototypical metric-induced capacity measure for $\mathcal{R}_{\lambda}$. Let us first discuss the case in which $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ is an interior point of $\mathcal{R}_{\lambda}$ (Case A). Since $\mathcal{R}_{\lambda}=\mathcal{E}_{\lambda}$, finding $S_{2}$ becomes the business of doing a hyperellipsoidal average of $\mathcal{D}$. This gets us to

[TABLE]

The logarithmic divergences in $\lambda$, a derivation byproduct from Gaussian approximation of $L$ and relaxation of $\partial\mathcal{R}_{0}$, pose no ill consequence so long as $N$ is sufficiently large such that $\mathcal{R}_{\lambda}\subset\mathcal{R}_{0}$ for all $\lambda$ values that give desirably large $C_{\lambda}<1$.

The situation becomes more complicated for Case B, which demands geometrical knowledge about $\partial\mathcal{R}_{0}$ for an exact calculation of $S_{2}$ (see Fig. 3). This tempts us to use a first-order approximation by expanding the likelihood $L$ about $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ to a Gaussian function of hyperellipsoidal-$\mathcal{E}^{\prime}_{\lambda}$ profile centered at $\boldsymbol{r}_{\text{c}}$, and next introducing a hyperplane containing $\widehat{\boldsymbol{r}}_{\textsc{ml}}$ that is tangent to its isoGaussian (constant-Gaussian-value) contour. $S_{2}$ is then a hyperellipsoidal-cap (formed by the hyperplane and the hyperellipsoid from the Gaussian expansion of $L$) average. We refer the Reader to Sec. VII of our companion article for all related technical calculations, and simply state the final formulas:

[TABLE]

involving $V_{d}=\pi^{d/2}/(d/2)!$, $l=\sqrt{\log(\lambda/\lambda^{\prime})/(-\log\lambda^{\prime})}$, $\mathcal{N}_{d,l,x}=V_{d}\,\mathrm{I}_{(1-l)/2}((d+x)/2,(d+x)/2)$ depending on the incomplete beta function $\mathrm{I}_{\cdot}(\cdot,\cdot)$, and

[TABLE]

It is easy to see that Eqs. (7) and (8) include Case A by recognizing that the “effective $\lambda$” ($\lambda^{\prime}$) approaches $\lambda$ ($\bm{g}_{\textsc{ml}}=\boldsymbol{0}$), so that $l\rightarrow 0$ gives $\mathcal{N}_{d,0,x}=V_{d}$ and $\boldsymbol{M}=(-\log\lambda)\,\boldsymbol{F}_{\textsc{ml}}^{-1}/(d+2)$.

*Discussions for quantum-state tomography.—*All results presented thus far apply to arbitrary physical systems. Here, we specifically investigate quantum-state tomography, thereby endowing explicit forms to all important quantities that are pertinent to Bayesian CR error certification.

For an unknown quantum state $\rho$ of Hilbert-space dimension $D$, every data-copy measurement in a tomography experiment is usually mutually independent, so that the log-likelihood $\log L=\sum^{M}_{j=1}n_{j}\log p_{j}$ catalogs the relative frequency data $\sum^{M}_{j=1}n_{j}=N$ of all $M$ measurement outcomes $\Pi_{j}\geq 0$ $(\sum_{j}\Pi_{j}=1)$, each with the Born probability $p_{j}=\mathrm{tr}\!\left\{\rho\Pi_{j}\right\}$. We can express $\rho$ and $\Pi_{j}$ in terms of the Hermitian basis $\{1/\sqrt{D},\Omega_{j}\}^{D^{2}-1}_{j=1}$ such that $\mathrm{tr}\!\left\{\Omega_{j}\right\}=0$ and $\mathrm{tr}\!\left\{\Omega_{j}\Omega_{k}\right\}=\delta_{j,k}$, so that we may denote the ($d=D^{2}-1$)-dimensional $\boldsymbol{r}=\mathrm{tr}\!\left\{\rho\,\boldsymbol{\Omega}\right\}$ and $\boldsymbol{q}_{j}=\mathrm{tr}\!\left\{\Pi_{j}\,\boldsymbol{\Omega}\right\}$. This leads to $\boldsymbol{F}_{\textsc{ml}}=N\sum^{M}_{j=1}\boldsymbol{q}_{j}\boldsymbol{q}^{\textsc{t}}_{j}/p_{\textsc{ml},j}$ ($N\gg 1$) and $\bm{g}_{\textsc{ml}}=\sum^{M}_{j=1}n_{j}\,\boldsymbol{q}_{j}/p_{\textsc{ml},j}$ for the ML state estimator $\widehat{\rho}_{\textsc{ml}}$ of ML probabilities $p_{\textsc{ml},j}=\mathrm{tr}\!\left\{\widehat{\rho}_{\textsc{ml}}\Pi_{j}\right\}$. In concrete terms, for Case A, $\widehat{\rho}_{\textsc{ml}}$ is full rank, such that the CR $\mathcal{R}_{\lambda}\approx\mathcal{E}_{\lambda}$; whereas for Case B, $\widehat{\rho}_{\textsc{ml}}$ is rank-deficient and $\mathcal{R}_{\lambda}\approx\mathcal{R}_{0}\cap\mathcal{E}^{\prime}_{\lambda}$ is therefore approximately a truncated $\mathcal{E}^{\prime}_{\lambda}$ (covariance profile of the Gaussian expansion of $L$ about $\widehat{\boldsymbol{r}}_{\textsc{ml}}$) by the quantum-state space $\mathcal{R}_{0}$—the convex set of unit-trace positive operators. The uniform $(\mathrm{d}\,\boldsymbol{r}^{\prime})$ is assumed.

To compare with the closed-form approximations in Eqs. (6) and (7), we pick the $l_{2}$-norm to measure the region capacity $S_{\textsc{hs}}\equiv S_{2}$ of $\mathcal{R}$, which is equivalent to the Hilbert-Schmidt (HS) distance for quantum states. We emphasize that for sufficiently large $N$, all arguments leading to the monotonicity of $S_{\mathcal{D},\lambda}$ still applies for Case B as $\bm{g}_{\textsc{ml}}\rightarrow\boldsymbol{0}$. Figures 4 and 5 showcase our in-region sampling theory. The matches in both Case A and B between theory and hit-and-run sampling are very good for moderate $D$, but are expected to have some discrepancies for more complex systems due to the more pronounced corners in $\partial\mathcal{R}_{0}$ Bengtsson et al. (2013). Instead, accelerated hit-and-run can be used, the complexity of which are analyzed in our companion article.

*Conclusions.—*In realistic multi-dimensional parameter estimation problems, sufficiently large dataset almost exclusively results in extremely small Bayesian credible regions relative to the entire parameter space. The conventional practice of first doing Monte Carlo to sample the parameter space followed by sample filtering almost always fails to accurately construct such small error regions. Our technique of in-region sampling developed in this Letter is capable of constructing any such small regions efficiently with perfect yield. In-region sampling is equivalent to computing region-averages that is efficient with a wide range of numerical methods. The region-average perspective of in-region sampling allows us to operationally formulate an alternative concept of region capacity through averaging any $l_{p}$ distance norm between two credible-region points, for which, in the special case $p=2$, closed-form approximation formulas to facilitate ultrafast analytical Bayesian error estimations with sufficiently large datasets are readily available. Either way, efficient Bayesian error certifications can now be carried out on physical systems of varying complexity. For exceedingly large quantum systems where Monte Carlo computations start to become visibly taxing, these asymptotic formulas can serve as large-scale approximate certifiers at least for high credibility values.