Adaptive compressive tomography: a numerical study

TL;DR

This paper numerically investigates an adaptive compressive tomography method that efficiently reconstructs low-rank quantum states with fewer measurements, outperforming some existing measurement strategies and proposing a hybrid approach for improved speed.

Contribution

It introduces a faster hybrid compressive tomography scheme combining random and adaptive measurements, with numerical evidence of favorable scaling and performance.

Findings

Adaptive scheme performs comparably to compressed sensing methods.

Outperforms random measurement bases for low-rank states.

Proposes a hybrid measurement approach for efficiency.

Abstract

We perform several numerical studies for our recently published adaptive compressive tomography scheme [D. Ahn et al. Phys. Rev. Lett. 122, 100404 (2019)], which significantly reduces the number of measurement settings to unambiguously reconstruct any rank-deficient state without any a priori knowledge besides its dimension. We show that both entangled and product bases chosen by our adaptive scheme perform comparably well with recently-known compressed-sensing element-probing measurements, and also beat random measurement bases for low-rank quantum states. We also numerically conjecture asymptotic scaling behaviors for this number as a function of the state rank for our adaptive schemes. These scaling formulas appear to be independent of the Hilbert space dimension. As a natural development, we establish a faster hybrid compressive scheme that first chooses random bases, and later adaptive bases as the scheme progresses. As an epilogue, we reiterate important elements of informational completeness for our adaptive scheme.

Full text

I Introduction

The main aim of quantum state tomography Chuang and Nielsen (2000); Paris and Řeháček (2004); Teo (2015) is to reconstruct the unknown true quantum state $\rho_{\text{t}}$ from data ${\mathbb{D}}$ obtained after a set of measurements (collectively represented by the map $\mathcal{M}[\rho_{\text{t}}]\mapsto\mathbb{D}$) are performed. Since quantum mechanics constrains $\rho_{\text{t}}$ to be a positive operator of unit trace, it is unambiguously specified by $d^{2}-1$ free parameters.

To fully characterize $\rho_{\text{t}}$, one performs an informationally complete (IC) set of $O(d^{2})$ measurement outcomes $\Pi_{j}\geq 0$ that form a positive operator-valued measure (POVM) ($\sum_{j}\Pi_{j}=1$). The corresponding physical probabilities $p_{\text{t},j}$ relate the POVM and $\rho_{\text{t}}$ via Born’s Rule, expressed by either $p_{\text{t},j}=\mathrm{tr}\{\rho_{\text{t}}\,\Pi_{j}\}$ or $\mathcal{M}[\rho_{\text{t}}]=\bm{p}_{\text{t}}$. Realistically, $\bm{p}_{\text{t}}$ is always inaccessible. Rather, $\mathbb{D}$ is collected from $N<\infty$ independent sampling copies, so that one acquires the relative frequency $\nu_{j}=n_{j}/N$ for each outcome $\Pi_{j}$ of $n_{j}$ data counts. Working with $\mathbb{D}\equiv\{\nu_{j}\}$, the maximum-likelihood (ML) method Řeháček et al. (2007); Teo et al. (2011); Teo (2015); Shang et al. (2017) may be used to supply the unique state estimator $\widehat{\rho}_{\textsc{ml}}\geq 0$ and its corresponding ML physical probabilities ($\bm{p}_{\textsc{ml}}$). In the limit $N\gg 1$, we have $\widehat{\rho}_{\textsc{ml}}\rightarrow\rho_{\text{t}}$. In both hypothetical noiseless and practical noisy situations, we say that operationally, $\mathbb{D}$ is IC if $\widehat{\rho}_{\textsc{ml}}$ is unique with respect to the whole state space Busch and Lahti (1989); Prugovečki (1977). Incidentally, this notion is synonymous to that of “strictly-complete” in Kalev et al. (2015); Baldwin et al. (2016).

For complex quantum systems in high-dimensional states, measuring $O(d^{2})$ outcomes quickly turns into a resource-intensive practical problem Häffner et al. (2005); Titchener et al. (2018). A prototypical strategy to circumvent this problem is to first assume a priori that $\mathrm{rank}\{\rho_{\text{t}}\}\leq r$ for a given small $r\ll d$, and next invoke the conventional method of compressed sensing (CS) Donoho (2006); Candès and Tao (2006a); Candès and Recht (2009) to find a unique rank-deficient state estimator Gross et al. (2010); Kalev et al. (2015); Baldwin et al. (2016); Steffens et al. (2017); Riofrío et al. (2017). This widely-held standard, however, faces two practical issues. To start off, the a priori assumption about $r$ requires additional verification before it can be applied to CS, since the resulting estimator accuracy hinges on the validity of this assumption. Next, as one has no means of validating the final estimator self-consistently, the usual solution is to compare the estimator with some assumed target state Kalev et al. (2015); Steffens et al. (2017); Riofrío et al. (2017). In the presence of experimental errors, there is simply no guarantee whether such a comparison is actually trustworthy.

In Ref. Ahn et al. (2019), we developed and experimentally carried out an adaptive compressive tomography (ACT) protocol that does not rely on any a priori information about the quantum state apart from its dimension $d$. It consists of a self-consistent informational completeness certification (ICC) stage and an adaptive measurement optimization stage, both of which operate only on accumulated data (we suggest the reference to Appendix A as the list of acronyms lengthens). For practical implementation of ACT, we considered adaptive choices of von Neumann (orthonormal) measurement bases that are generally more easily realized in the laboratory than arbitrary POVMs, as well as the product variant PACT that is compatible with many-body systems by utilizing local product bases instead of entangled ones as in ACT.

In this article, we shall perform a series of numerical studies to explore (P)ACT. In particular (a) we shall demonstrate, up to a reasonably large $d$ and $r$, that in terms of the number of IC bases $k_{\textsc{ic}}$ needed for unique reconstructions, our adaptive schemes outperform known random bases measurements. (b) We next show that the $k_{\textsc{ic}}$ scaling behaviors of both ACT and PACT match the respective ones of the element-probing Baldwin–Flammia POVMs and Baldwin–Goyeneche bases reported in Baldwin et al. (2016), thereby effectively provides conjectured asymptotic scaling behaviors for both adaptive schemes in the limits $r\ll d$ and $d\gg 2$. More specifically, we have $k_{\textsc{ic}}=2r+2$ for ACT and $k_{\textsc{ic}}=4r+1$ for PACT. Finally, (c) we present a “random-adaptive” hybrid compressive tomography scheme that combines the respective speed and compressive advantages of random and adaptive strategies.

The article is organized as follows. Sections II and III respectively provide a preliminary introduction to CS and ACT for subsequent discussions. After providing some details about the simulation specifications in Sec. IV, all numerical results are then presented in Sec. V. We end the article with an additional remarks on the concept of “IC” in ACT in Sec. VI.

II Compressed-sensing tomography

The concept of CS Donoho (2006); Candès and Tao (2006a); Candès and Recht (2009) provides one approach to recover an unknown signal of known degree of sparsity by performing a small set of specialized compressive measurement outcomes, and has since been widely adopted in signal processing Rani et al. (2019). In the context of quantum state tomography Gross et al. (2010); Kalev et al. (2015); Baldwin et al. (2016); Steffens et al. (2017); Riofrío et al. (2017), the low-rank assumption often taken for granted, that is $\mathrm{rank}\{\rho_{\text{t}}\}\leq r$ for a known and sufficiently small integer $r$, permits the utilization of CS to uniquely reconstruct $\rho_{\text{t}}$ with very few outcomes. While the primary focus began with random Pauli observables $M=O(rd(\log\,d)^{2})$ Gross et al. (2010); Flammia et al. (2012), compressive measurements of more variety were later constructed Baldwin et al. (2016).

II.1 Summary of the compressed-sensing procedure

Based on the assumed prior information, for measurement map $\mathcal{M}$ to be compressive and IC with respect to $\rho_{\text{t}}$, it is sufficient to satisfy the so-called restricted isometry property Candès and Tao (2005); Candès et al. (2006). Without spelling out unnecessary mathematical details, we note that there also exist other kinds of compressive measurements without such a property Baldwin et al. (2016). With a compressive $\mathcal{M}$, one seeks the unique estimator $\widehat{\rho}$ from the data convex set $\mathcal{C}$, which contains all possible states that satisfies the constraint $\mathcal{M}[\widehat{\rho}]=\bm{p}_{\text{t}}$ in the noiseless regime.

One possible way is to choose $\widehat{\rho}$ that has the lowest rank out of $\mathcal{C}$. It can be shown that for $\mathcal{M}$ satisfying the restricted isometry property above certain threshold degree of orthonormality, this nonconvex optimization is equivalent to the convex optimization of the trace-norm (nuclear-norm) minimization Recht et al. (2010), where the rank and the trace-norm are respectively the operator version of the $l_{0}$-norm and the $l_{1}$-norm meant for signal recovery Candès and Tao (2006a). This ensures that the set of states consistent with $\mathbb{D}$ intersects the state space at exactly one lowest-rank state that coincides with $\rho_{\text{t}}$ if $N\gg 1$. With real dataset $\mathbb{D}$ and positivity (see also Fig. 1), the CS algorithm (proven to work with measurements possessing restricted isometry property) that uniquely identifies a state estimator may be simplified Kalev et al. (2015) into

**CS procedure with positivity

** For known $r$ and data $\mathbb{D}=\{\nu_{j}\}\equiv\bm{\nu}$:

Minimize $\mathrm{tr}\!\left\{\rho^{\prime}\right\}$, subject to

• •

$\|\mathcal{M}[\rho^{\prime}]-\bm{\nu}\|\leq\epsilon$ for some $\epsilon>0$ that depends on noise,

• •

$\rho^{\prime}\geq 0$. 2. 2.

Trace-renormalize the optimal $\rho^{\prime}$.

It can be shown that the above optimization routine leads to perfect recovery of $\rho_{\text{t}}$ for noiseless data ($\bm{\nu}\rightarrow\bm{p}_{\text{t}}$ with $\epsilon\rightarrow 0$).

II.2 Typical compressive measurements

For subsequent numerical comparisons, we shall consider a few well-known compressive POVMs that were applied to quantum state tomography. The first classic compressive measurement is the set of random Pauli bases (RP) for $n$-qubit systems Gross et al. (2010); Kalev et al. (2015), which maximally comprise $3^{n}$ sets of $n$ tensor products of the standard qubit Pauli operators $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$.

On the other hand, the independent studies of two other types of measurements Flammia et al. (2005); Goyeneche et al. (2015) related to pure-state distinction eventually led to the construction of rank-$r$ generalized compressive measurements Baldwin et al. (2016). These are the element-probing Baldwin–Flammia (BF) POVM and Baldwin–Goyeneche (BG) bases that are constructed using the mathematical concepts of Schur complement and block-diagonalization of the density matrix. More importantly, these compressive measurements have analytical performance scaling behaviors in the regime $r\ll d$. Upon denoting the number of IC outcomes as $M_{\textsc{ic}}$, the BF POVM gives $M_{\textsc{ic}}=(2d-r)r+1$ [two outcomes more than the number of free parameters $(2d-r)r-1$ for a rank-$r$ state], and the BG bases give $k_{\textsc{ic}}=4r+1$, or $M_{\textsc{ic}}=(4r+1)d$.

Random measurements first fueled the progress of CS Candès and Tao (2006b). For state tomography, we look into two classes of random measurements. The first class is the set of random bases generated by Haar unitary operators (RH) that are widely used in quantum information theory Mezzadri (2007); Ćwikliński et al. (2013); Russell et al. (2017); Banchi et al. (2018). The second class is the set of eigenbases of random full-rank states (RS) distributed uniformly according to the Hilbert-Schmidt measure Życzkowski and Sommers (2003). For more details regarding the numerical constructions of random bases, we refer the interested Reader to the Appendix B. Finally, there exists yet another benchmark by Kech and Wolf (KW) for arbitrary von Neumann bases, which states that $k\geq 4r\lceil\frac{d-r}{d-1}\rceil$ is sufficient to uniquely reconstruct any rank-$r$ state Kech and Wolf (2016).

II.3 Issues in conventional compressed sensing

While the CS procedure is a promising candidate for low-rank state recovery. There are two primary concerns that need to be addressed.

First, the proper set of compressive measurements is on the premise that the upper bound $r$ of the rank of $\rho_{\text{t}}$ is known. The assumption of $r$ must thus be experimentally justified. A categorical misclassification of $\rho_{\text{t}}$ would result in either a non-IC (too small an $r$) or an unnecessarily overcomplete dataset (too big an $r$). The former leads to an ambiguous set of estimators, while the latter overuses measurement resources.

Second, there exists no method in CS to systematically validate if the acquired $\mathbb{D}$ is truly IC in real experiments. Instead, the fidelity measure is commonly used as the indicator that the estimator is “close” to some prechosen target state. This approach, at the very least, requires yet another round of certification for these target states, and is evidently not the correct informational completeness indicator as it registers no information about the data convex set $\mathcal{C}$.

III Adaptive compressive tomography

Our recently proposed ACT Ahn et al. (2019) is able to deterministically compress IC datasets for any given rank-$r$ $\rho_{\text{t}}\equiv\rho_{r}$ without relying on any a priori information or ad hoc assumptions about $\rho_{r}$ apart from its dimension $d$. The procedure of ACT is iterative and involves two main stages in every step. The first informational completeness certification stage (ICC) decides, given the measured outcomes and accumulated data $\mathbb{D}$, whether the data convex set $\mathcal{C}$, which is again the set of states $\mathcal{C}=\{\rho^{\prime}\}$ for which $\mathcal{M}[\rho^{\prime}]=\bm{p}_{\textsc{ml}}$ is singleton or not. The second adaptive strategy provides protocols which adaptively chooses the next measurement at each step of ACT according to the accumulated data $\mathbb{D}$. The former would imply a unique state estimator consistent with an IC $\mathbb{D}$, and the procedure terminates. The latter would introduce appropriate protocols which significantly reduces the size of IC $\mathbb{D}$ compared to $O(d^{2})$.

Throughout the article, we shall take the POVM to be sets of von Neumann bases that are each denoted by $\mathcal{B}$ and contains $d$ orthonormal projectors, which are the typical kind of measurement employed in laboratory. A particular iterative adaptation step of ACT is then to search for an optimal basis to measure in the next step (see Fig. 2).

III.1 Informational completeness certification

As ACT progresses iteratively and a sequence of von Neumann bases $\mathcal{B}=\{\mathcal{B}_{1},\mathcal{B}_{2},\ldots,\mathcal{B}_{k}\}$ are measured, the size $\zeta_{k}$ of the data convex set $\mathcal{C}_{k}$ at the $k$th iterative step of ACT indicates whether the corresponding dataset $\mathbb{D}_{k}$ of $k$ measured orthonormal bases is IC or not—$\mathbb{D}_{k}$ is IC if and only if $\zeta_{k}=0$. In this case, we denote the set of IC bases to be $k_{\textsc{ic}}=k$. In other words, the task of ICC is to verify if $\zeta=0$ (or close to some small numerical threshold value in practice).

Fortunately, whilst $\zeta_{k}$ is a complicated state-space integral with respect to some volume measure over $\mathcal{C}_{k}$, there is an alternatively much more feasible way to detect if $\zeta_{k}=0$. Using a randomly chosen full-rank state $Z$ throughout the entire run of ACT, if we define the linear function $f(\rho^{\prime})=\mathrm{tr}\{\rho^{\prime}\,Z\}$ for the variable state $\rho^{\prime}$, then owing to the convexity of $\mathcal{C}_{k}$ regardless of whether data noise is present, it straightforwardly follows that $\zeta_{k}=0$ iff $f_{\text{min},k}=f_{\text{max},k}$, where $f_{\text{min},k}$ and $f_{\text{max},k}$ are respectively the minimum and maximum values of $f$ in $\mathcal{C}_{k}$ (see Fig. 3). Hence, the validity of this if-and-only-if statement is robust against noise. This simple conclusion holds as long as $Z$ is neither $1/d$ nor possesses a kernel that contains $\rho^{\prime}\in\mathcal{C}_{k}$, which are all measure-zero situations. Simply put, the business of ICC is a semidefinite program Vandenberghe and Boyd (1996) that evaluates the quantity $s_{\textsc{cvx},k}=(f_{\text{max},k}-f_{\text{min},k})/(f_{\text{max},1}-f_{\text{min},1})$ at every $k$th iterative step of ACT.

ICC in the $k$th step

Maximize and minimize $f_{Z}(\rho^{\prime})=\mathrm{tr}\!\left\{\rho^{\prime}Z\right\}$ for a fixed, randomly-chosen full-rank state $Z\neq 1/d$ to obtain $f_{\text{max},k}$ and $f_{\text{min},k}$ subject to

• •

$\rho^{\prime}\geq 0$ , $\mathrm{tr}\!\left\{\rho^{\prime}\right\}=1$ ,

• •

$\mathcal{M}[\rho^{\prime}]=\bm{p}_{\textsc{ml},k}$ for $\bm{p}_{\textsc{ml},k}$ obtained from the accumulated $\mathbb{D}_{k}$ thus far . 2. 2.

Compute $0\leq s_{\textsc{cvx},k}\leq 1$ and check if it is smaller than some threshold $\epsilon$. 3. 3.

If $s_{\textsc{cvx},k}<\epsilon$, terminate ACT. Continue otherwise.

More specifically, $s_{\textsc{cvx}}$ is known as the size monotone in the sense that with decreasing $\zeta_{k}$ or increasing $k$, $s_{\textsc{cvx},k}$ never increases in value for near-perfect $\mathbb{D}_{k}$. The reason for this is that with perfect $\mathbb{D}$s, we have the set hierarchy $\mathcal{C}_{1}\supseteq\mathcal{C}_{2}\supseteq\ldots\supseteq\mathcal{C}_{k_{\textsc{ic}}}$ as more linearly-independent bases are measured, so that $f_{\text{max},k}$ monotonically decreases and $f_{\text{min},k}$ monotonically increases for a linear $f$, or $s_{\textsc{cvx},k+1}\leq s_{\textsc{cvx},k}$. Moreover, it is easy to see that if $f_{\text{max},k+1}-f_{\text{min},k+1}<f_{\text{max},k}-f_{\text{min},k}$, then $\mathcal{C}_{k+1}\subset\mathcal{C}_{k}$.

In terms of complexity, unlike the traditional NP-hard problem to determine whether a given set of POVM possesses the restricted isometry property Bandeira et al. (2013), the complexity to decide if this POVM is IC with respect to $\mathbb{D}$ turns out to be, understandably, only as high as that for the semidefinite program employed in ICC.

III.2 Adaptive strategy

Given the collected dataset $\mathbb{D}_{k}$ and corresponding non-singleton data convex set $\mathcal{C}_{k}$ at the $k$th iterative step of ACT, one seeks the optimal basis to be measured in the $(k+1)$th step in order to converge ACT at a reasonably quick pace to the singleton $\mathcal{C}_{k=k_{\textsc{ic}}}$.

If $\rho_{r}$ is a low-rank state ($r\ll d$), the natural motivation for a compressive adaptive strategy, given no other assumption, would be to actively seek the lowest-rank estimator $\widehat{\rho}_{k}$ from the data convex set $\mathcal{C}_{k}$ that ultimately either converges to $\rho_{t}$ for noiseless data, or is very close to it for $N\gg 1$. This is akin to conventional CS where for a properly chosen compressive measurement map $\mathcal{M}$, the unique CS estimator, that converges to $\rho_{r}$ for noiseless data, essentially has the lowest rank out of data convex set. In order to establish a feasible scheme for finding the lowest-rank $\widehat{\rho}_{k}$ in $\mathcal{C}_{k}$, we minimize a concave function over $\mathcal{C}_{k}$ that has rank-minimizing characteristics. For this purpose, two exemplifying concave functions with a global minimum for pure states shall be numerically studied. They are the von Neumann entropy function $S(\rho^{\prime})=-\mathrm{tr}\!\left\{\rho^{\prime}\,\log\rho^{\prime}\right\}$, and the linear entropy function $S_{\textsc{l}}(\rho^{\prime})=1-\mathrm{tr}\!\left\{\rho^{\prime 2}\right\}$, both of which are similar in value near all pure states and take the minimum value 0 for the pure states. We mention that entropy minimization was numerically reported in Tran et al. (2016) to offer a stronger compressive recovery for general low-rank matrices compared to the standard trace-norm minimization in CS.

While in CS, the measurement map is first decided before a unique lowest-rank estimator is obtained, we note here that for ACT, the mechanism that drives the compressive nature is the interplay between positivity and data (ML) constraints on the entropy minimization procedure, which leads to the correct guidance to the unique $\widehat{\rho}_{k_{\textsc{ic}}}\approx\rho_{r}$. It is shown in Ahn et al. (2019) that owing to these two constraints, if one manages to measure the eigenbasis $\mathcal{B}_{\rho_{r}}$ of $\rho_{r}$ that is low-rank ($r\ll d$), then only $k_{0}=\lceil(r^{2}-r)/(d-1)\rceil+1\ll d+1$ bases (including $\mathcal{B}_{\rho_{r}}$) are needed to unambiguously reconstruct $\rho_{r}$. Therefore the eigenbasis of a lowest-(linear-)entropy $\widehat{\rho}_{k}$ (which would essentially have an extremely low rank) in $\mathcal{C}_{k}$ may be assigned to $\mathcal{B}_{k+1}$. In this manner, if $N\gg 1$, it is evident that $\mathcal{B}_{k+1}$ approaches $\mathcal{B}_{\rho_{r}}$, and the action of the optimization procedure under both constraints aid in speeding up this convergence (see Fig. 4).

Generally, the adaptive strategy discussed here yields a sequence of entangled bases. For many-body quantum systems, such bases are often difficult to be realized experimentally. To carry out the adaptive strategy for such systems, we impose, in every iterative step, an additional local structure on the adaptive von Neumann bases by searching for the nearest tensor-product basis that is “closest” to the optimal entangled one. One way to do so is to first minimize the distance between the lowest-(linear-)entropy $\widehat{\rho}_{k}\in\mathcal{C}_{k}$ and another state $\rho^{\prime}$ with product eigenbasis with respect to some metric, and take the eigenbasis of the optimum. This distance reduction helps to find a local von Neumann basis that is close to $\mathcal{B}_{\widehat{\rho}_{k}}$ for the lowest-(linear-)entropy $\widehat{\rho}_{k}$.

Putting everything together, we arrive at the complete (P)ACT procedure:

**(P)ACT procedure

**Beginning with $k=1$ and a random computational basis $\mathcal{B}_{1}$:

Measure $\mathcal{B}_{k}$ and collect the relative frequency data $\sum^{d-1}_{j^{\prime}=0}\nu_{j^{\prime}k}=1$. 2. 2.

From $\left\{\nu_{0k^{\prime}},\ldots,\nu_{d-1\,\,k^{\prime}}\right\}^{k}_{k^{\prime}=1}$, obtain $kd$ physical ML probabilities. 3. 3.

Perform ICC with the ML probabilities and compute $s_{\textsc{cvx},k}$:

• •

If $s_{\textsc{cvx},k}<\epsilon$, terminate ACT and take $\rho_{\text{max}}\approx\rho_{\text{min}}$ as the estimator and report $s_{\textsc{cvx},k}$.

• •

Else Proceed. 4. 4.

Choose a lowest-(linear-)entropy $\widehat{\rho}_{k}\in\mathcal{C}_{k}$ in $\mathcal{C}_{k}$ 5. 5.

Define $\mathcal{B}_{k+1}$ to be the eigenbasis of $\widehat{\rho}_{k}$ for ACT, or a local basis close to it for PACT via some prechosen distance minimization technique. 6. 6.

Set $k=k+1$ and repeat.

IV Simulation Specifications

In this brief section, we clarify all essential technical details involved in our subsequent numerical studies. We begin by summarizing the main goals of our numerical studies:

1. (A)

Compare (P)ACT with the following random measurement schemes (refer to Sec. II.2):

1. (I)

Random Pauli (RP) bases. 2. (II)

Random Haar-uniform (RH) bases. 3. (III)

Eigenbases of random full-rank states (RS). 2. (B)

Benchmark (P)ACT using known analytical scalings of

1. (IV)

Baldwin–Flammia (BF) POVMs, 2. (V)

Baldwin–Goyeneche (BG) bases, 3. (VI)

Kech–Wolf scaling (KW) for arbitrary bases,

and conjecture asymptotic $k_{\textsc{ic}}$ scalings for (P)ACT. 3. (C)

Present and numerically benchmark a new hybrid compressive scheme (HCT).

In what follows, all simulations are performed on ensembles of random quantum states of various ranks $r$ over which important indicators such as $s_{\textsc{cvx}}$ and $k_{\textsc{ic}}$ are averaged. To generate these ensembles, we follow the prescriptions in Życzkowski and Sommers (2003) and distribute the random states according to the Hilbert-Schmidt measure. For a fixed $r$, this is done by first generating $r\times d$ matrices $A_{j}$ with entries i.i.d. standard Gaussian distribution, and next define the ensemble density matrices $\{\rho_{j}\}$ in accordance with $\rho_{j}=A_{j}^{\dagger}A_{j}/\mathrm{tr}\!\left\{A_{j}^{\dagger}A_{j}\right\}$.

Both minimum entropy schemes (with $S$ and $S_{\textsc{l}}$) over convex sets are carried out efficiently using the superfast accelerated projected gradient algorithm Shang et al. (2017). The ICC algorithm is carried out with the CVXOPT package Grant and Boyd (2014, 2008). While the deterministic compressive measurements, namely the BF POVMs and BG bases, possess analytical $k_{\textsc{ic}}$ scaling behaviors that can be readily used for benchmarking purposes, the random measurements, namely RP, RH and RS bases have at most approximated scaling expressions. The generations of both RH and RS bases are described in the appendix. For all local bases schemes, which refer to PACT, (I), (II), and (III), since the quantum systems that we shall be investigating are multi-qubit, these schemes involve basis outcomes defined as projectors onto the tensor-product space of single-qubits.

For the numerical simulation tasks (a)–(c), the results of which are presented in Sec. V, we study all results for noiseless data to understand the underlying theoretical relationships between $\{s_{\textsc{cvx}},k_{\textsc{ic}}\}$ and $r$.

V Numerical results

V.1 Comparisons of (product) adaptive compressive tomography with random-bases measurements

We proceed to systematically compare both ACT and PACT respectively with popular random compressive schemes (I)–(III) in Sec. IV. Figure 5 shows the results of ICC for true quantum states of 4-, 5- and 6-qubit systems. For each type of system, ICC is applied to all schemes on states of ranks $r=2$, 4 and 6. For all the schemes, we confirm that $s_{\textsc{cvx},k}$ monotonically decreases over the number of measured von Neumann bases $k$, and reach $s_{\textsc{cvx},k_{\textsc{ic}}}=0$ (up to some numerical threshold). In terms of the convergence rate of $s_{\textsc{cvx},k}$, it turns out that (P)ACT as well as schemes (II)–(V) are more efficient than RP. Here, ACT scheme shows most efficient performance. Faster convergence of $s_{\textsc{cvx},k}$ directly leads to smaller $k_{\textsc{ic}}$, which implies higher compression in size of IC data.

The compression efficiencies for all schemes are compared via $k_{\textsc{ic}}$ more clearly in Fig. 6, which validates that ACT is the most efficient scheme for all number of qubits, whereas RP turns out to be the most inefficient one. More specifically, it is evident that the $k_{\textsc{ic}}$ gap between ACT and RH, and that between PACT/RH (local) and RP enlarges for larger $r$ and number of qubits, which implies that ACT is relatively more for states of higher rank and Hilbert-space dimension. The apparent coincidence between PACT and RH (local) for all $d$ may be attributed to an unavoidable minimum level of incurred randomness in PACT that arises from the restriction to local-bases measurement during adaptive optimization. Such a level of randomness is sufficient to practically render PACT equivalent to RH.

V.2 Benchmarks for (product) adaptive compressive tomography and asymptotic performances

Benchmarking of the $k_{\textsc{ic}}$ performance for ACT, RH, RS, PACT and local RH schemes is done with the known analytical standards (IV)–(VI) listed in Sec. IV, and the results are presented in Fig. 7. The ACT schemes carried out through minimizing $S$ and $S_{\textsc{l}}$ are compared with the element-probing BF POVMs for arbitrary rank-$r$ states, which contains outcomes of non-unit rank and possess a total number of $M=(2dr-r^{2}+1)$ outcomes. The figure shows that for both ACT protocols, $k_{\textsc{ic}}$ scaling behaviors are in good agreement with the expression $M/d+2$ (a shifted BF scaling). This means that ACT, which yields only rank-1 von Neumann projectors, exceeds only by 2 bases in performance as compared to BF POVMs to perform IC state tomography.

We also note that both ACT schemes defined by the minimization of two different entropy functions give almost indistinguishable $k_{\textsc{ic}}$ curves for all tested $r$ and $d$. This numerical observation confirms that the behaviors of von Neumann and linear entropies are almost the same with respect to the optimization algorithm in the adaptive stage. Thus, the (P)ACT protocols involving the von Neumann entropy minimization is sufficient for the remaining discussions. Indistinguishable $k_{\textsc{ic}}$ curves exist also for random schemes (such as RH and RS), and this amply hints that the effect of informational completeness for the non-adaptive schemes depends weakly on the specific choice of random bases generation algorithm, but rather more strongly on the state rank.

One might, at this point, naively suppose that perhaps ACT should also work sufficiently well if $\widehat{\rho}_{k}$ is some randomly chosen rank-deficient state in $\mathcal{C}_{k}$ for every $k$. Fig. 8 compares the average $k_{\textsc{ic}}$ of min-$S$ adaptive bases with that of eigenbases of random rank-deficient states. It presents clear numerical evidence contrary to the above supposition, which further reinforces the statement that unlike random schemes, an appropriate choice of objective function is necessary for ACT to achieve optimal compression.

In the context of many-body local compressive measurements, we compare the $k_{\textsc{ic}}$s of PACT and local RH scheme with those of the BG bases for arbitrary $r$ and $d$, the latter which employs $4r+1$ specifically constructed orthonormal entangled bases. Figure 7 tells us that both product schemes asymptotically approach the BG scheme in $k_{\textsc{ic}}$ performance as the number of qubits grows.

Hence, we numerically confirm that for all $r$ and $d$, ACT shows stronger compressive behavior than PACT and non-adaptive random schemes. The result of benchmarking with BF POVMs and BG bases reveals $k_{\textsc{ic}}$ behaviors for ACT and PACT that lose the dependency on $d$ in limit of large number of qubits. These conjectured asymptotic scalings are respectively $k_{\textsc{ic}}=2r+2$ and $k_{\textsc{ic}}=4r+1$, and our belief in their validity is further strengthened with 7-qubit systems in Fig. 9.

All the compressive bases schemes presented here are contrasted with the KW benchmark that applies to arbitrary bases measurements, which is effectively a linear function of $r$ that is independent of $d$ due to the ceiling function. We point out that the KW benchmark that estimates the required number of bases needed to perform IC reconstruction using arbitrary von Neumann bases almost always overestimates the $k_{\textsc{ic}}$ performance of ACT. The reason is that the measurement resources for ACT scales nonlinearly with $r$ owing to the low-rank guidance from both positivity and data constraints. The KW benchmark also show signs of overestimation for the random RH and RS schemes for sufficiently large $r$.

V.3 Hybrid compressive tomography

While the complete ACT scheme is highly compressive ($k_{\textsc{ic}}=2r+2$), the computational resources needed to search for the optimal estimators in the data convex set will eventually become noticeably expensive for extremely large physical systems. On the other hand, a fully random scheme simply suggests measurement bases randomly and requires essentially negligible computational resources.

As an attempt to speed up the compression process, we may take advantage of the benefits from both random and adaptive protocols and establish a hybrid compressive tomography (HCT) scheme. This hybrid scheme starts with choosing random von Neumann measurement bases, which may be constructed with either the RH or RS prescription. In fact, any sort of bases would deliver similar performance due to the limited data one always has at the beginning of a compression process. As the accumulated dataset is built up, more information is gained about the unknown $\rho_{r}$ such that it is now justified to spend more effort in searching for good optimal measurements based on more statistically reliable data. Armed with this basic insight, we propose the following procedure:

**HCT procedure

**Beginning with $k=1$, a random computational basis $\mathcal{B}_{1}$ and some positive threshold value $s_{\mathrm{th}}$:

Measure $\mathcal{B}_{k}$ and collect the relative frequency data $\sum^{d-1}_{j^{\prime}=0}\nu_{j^{\prime}k}=1$. 2. 2.

From $\left\{\nu_{0k^{\prime}},\ldots,\nu_{d-1\,\,k^{\prime}}\right\}^{k}_{k^{\prime}=1}$, obtain $kd$ physical ML probabilities. 3. 3.

Perform ICC with the ML probabilities and compute $s_{\textsc{cvx},k}$:

• •

If $s_{\textsc{cvx},k}<\epsilon$, terminate ACT and take $\rho_{\text{max}}\approx\rho_{\text{min}}$ as the estimator and report $s_{\textsc{cvx},k}$.

• •

Else Proceed. 4. 4.

If $s_{\textsc{cvx},k}>s_{\mathrm{th}}$:

• •

Assign a random basis, perhaps from Appendix B, to $\mathcal{B}_{k+1}$.

Else:

• •

Assign an optimal basis obtained using the adaptive strategy in (P)ACT to $\mathcal{B}_{k+1}$. 5. 5.

Set $k=k+1$ and repeat.

Figure 10 reveals that both RH–ACT and RS–ACT type HCT schemes give almost the same averaged $k_{\textsc{ic}}$ behaviors, both of which are essentially identical to ACT as in Fig. 7. Based on these numerical evidences, we may argue that the cause of these equivalent performances comes from almost identical gradual average gains in tomographic information throughout the respective compressive tomography schemes. This supports our observation that, without relying on any a priori knowledge about $\rho_{r}$, adaptive methods only serve as crucial roles in compressive tomography after sufficient amount of tomographic information is acquired for these methods to yield more correct optimized measurement bases. Otherwise, the identical average compression capabilities of both ACT and HCT imply that both schemes effectively acquire the same amount of tomographic information at the early stages of the processes. This justifies the more economical approach of first measuring random bases during the initial stage of HCT, which is far less time-consuming than optimizing for adaptive bases in ACT on highly complex quantum systems.

VI Epilogue: Remarks on informational completeness in adaptive compressive tomography

Besides $s_{\textsc{cvx}}$, we may inspect other aspects of the data convex set in ACT, carried out with $S$ minimization for example, to preempt its correct termination in case one intends to stop the reconstruction procedure with $k<k_{\textsc{ic}}$ measurement bases. Figure 11 showcases the average behaviors of various different minimum-entropy $\widehat{\rho}_{k}$ qualities for ACT under negligible statistical fluctuation ($N\gg 1$). We see that the fidelity, rank and entropy of $\widehat{\rho}_{k}$ reach the correct true values as $k$ increases, and saturate at $k<k_{\textsc{ic}}$. Technically, both the fidelity and rank cannot be used to judge the termination of ACT, since the fidelity requires knowledge of $\rho_{r}$, and it is also impossible to guess whether the rank of $\widehat{\rho}_{k}$ is close to the correct value without $\mathrm{rank}\{\rho_{r}\}$ as reference because of its regular stepped gradient. However, based on our numerical experience, the monotonically increasing entropy $S(\widehat{\rho}_{k})$ (a direct mathematical consequence of the inequality chain $\mathcal{C}_{1}\supseteq\mathcal{C}_{2}\supseteq\ldots\supseteq\mathcal{C}_{k_{\textsc{ic}}}$ for noiseless data) always approaches the final true value with smoothly decreasing gradient $S(\widehat{\rho}_{k})-S({\widehat{\rho}_{k-1}})<S(\widehat{\rho}_{k-1})-S(\widehat{\rho}_{k-2})$. We may then attempt to prematurely stop ACT at $k=k_{0}$ when $S(\widehat{\rho}_{k_{0}})-S(\widehat{\rho}_{k_{0}-1})$ is small enough and use the resulting estimator $\widehat{\rho}_{k_{0}}$ as an approximately IC reconstruction for future statistical predictions, in which case this low-rank $\widehat{\rho}_{k_{0}}$ will be close to $\rho_{r}$ even though the size of $\mathcal{C}_{k_{0}}$ can be appreciably greater than zero. With statistical consistency, this termination method works also for real data of sufficiently large $N$ (or low statistical fluctuation).

The aforementioned remarks eventually lead to a subtle point behind the notion of “IC” adopted by ACT. In general tomography settings, which include the studies in CS, one typically defines “IC” in terms of the class of quantum states of interest, such as the class of rank-$r$ states to name one. In ACT, however, we speak of an IC dataset $\mathbb{D}$ that is statistically related to the particular unknown state we are probing. As a result, the notion of “IC” in ACT, along with the value of $k_{\textsc{ic}}$, strongly depends on both the adaptive bases and $\mathbb{D}$.

A more objective $\mathbb{D}$- or $\rho_{r}$-dependent compressive strategy in exchange for spurious or ad hoc prior assumptions about $\rho_{r}$ also means that an ACT bases POVM that is IC for a given unknown rank-$r$ state is not necessarily going to be IC for another rank-$r$ state. Figure 12 elucidates this fact for a distribution of rank-one true states probed by an ACT bases POVM that is IC for a particular pure state assuming noiseless probabilities.

VII Conclusions

We have performed a comprehensive numerical study of our adaptive compressive tomography schemes, which can uniquely reconstruct any arbitrary rank-deficient quantum state with few measurement bases than the square of its given dimension without any other a priori information about this state. Several numerical results in this article can serve as important guidelines for resource planning, structuring and executing objective compressive tomography for multi-qubit systems, which are widely used in quantum information and quantum optics studies.

After simulating up to reasonably large dimensions $d$ of multi-qubit systems and quantum-state rank $r$, we can summarize these numerical results under low statistical fluctuation:

The average performances (minimal number of informationally complete bases $k_{\textsc{ic}}$) of adaptive compressive tomography schemes almost always beat those of random compressive schemes, which include random Pauli bases, random Haar-unitary bases and eigenbases of random Hilbert-Schmidt-uniform states. 2. 2.

The average $M_{\textsc{ic}}=k_{\textsc{ic}}d$ for the entangled adaptive compressive scheme is greater than that of Baldwin-Flammia measurement outcomes by only $2d$ for all tested range of $r$ and $d$. The asymptotic ($d\gg 2$) average $k_{\textsc{ic}}$ for both entangled and product adaptive schemes are respectively $2r+2$ and $4r+1$. 3. 3.

There is virtually no difference in average $k_{\textsc{ic}}$ performance between a fully adaptive scheme and a hybrid scheme that starts off as a random scheme followed by an adaptive scheme at a later stage for a reasonable transition point. Therefore, the hybrid compressive scheme may be used to speed up the tomography compression process.

acknowledgments

We 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 NRF grant funded by the Korea government (MSIP) (Grant No. 2010-0018295), the European Research Council (Advanced Grant PACART), the Spanish MINECO (Grant FIS2015-67963-P), the Grant Agency of the Czech Republic (Grant No. 18-04291S), and the IGA Project of the Palacký University (Grant No. IGA PrF 2019-007).

Appendix A List of acronyms

[TABLE]

Appendix B Constructions of random bases

We briefly supply two short routines to construct von Neumann bases that are generated by random Haar-distributed unitary operators (RH bases), as well as those that are eigenbases of random full-rank states uniformly distributed according to the Hilbert-Schmidt measure (RS bases). These two unitary sets have generally different operator probability distributions, which may be verified by comparing some of their operator moments.

It is well-known that the QR decomposition generates unitary operators distributed according to the Haar measure Mezzadri (2007), so that the following routine applies:

**Constructing an RH basis

** Starting from a reference basis $\mathcal{B}_{\text{ref}}=\{\left|0\right>,\left|1\right>,\ldots,\left|d-1\right>\}$:

Generate a random $d\times d$ matrix $A$ with entries i.i.d. standard Gaussian distribution. 2. 2.

Compute $Q$ and $R$ from the QR decomposition $A=QR$. 3. 3.

Define $R_{\text{diag}}=\mathrm{diag}\{R\}$. 4. 4.

Define $L=R_{\text{diag}}\oslash|R_{\text{diag}}|$ ($\oslash$ refers to the Hadamard division). 5. 5.

Define $U_{\mathrm{Haar}}=QL$.

Thereafter construct the new basis $\mathcal{B}_{\mathrm{Haar}}=U_{\mathrm{Haar}}\mathcal{B}_{\text{ref}}$.

Steps 3–5 enforces a QR decomposition procedure that produces an effective “$R$” matrix that has positive diagonal entries, which is the correct decomposition procedure we need to generate $U_{\mathrm{Haar}}$.

The second class may be easily generated by the following pseudocode:

**Constructing an RS basis

** Starting from a reference basis $\mathcal{B}_{\text{ref}}=\{\left|0\right>,\left|1\right>,\ldots,\left|d-1\right>\}$:

Generate a random $d\times d$ matrix $A$ with entries i.i.d. standard Gaussian distribution. 2. 2.

Define $\rho_{\textsc{hs}}=\dfrac{A^{\dagger}A}{\mathrm{tr}\!\left\{A^{\dagger}A\right\}}$. 3. 3.

Obtain $U_{\textsc{hs}}$ from the spectral decomposition $\rho_{\textsc{hs}}=U_{\textsc{hs}}D_{\textsc{hs}}U_{\textsc{hs}}^{\dagger}$ of diagonal $D_{\textsc{hs}}$.

Thereafter construct the new basis $\mathcal{B}_{\textsc{hs}}=U_{\textsc{hs}}\mathcal{B}_{\text{ref}}$.