Fast state tomography with optimal error bounds

We analyze several important and practically relevant measurement schemes: structured POVMs (e.g. symmetrically informationally complete (SIC)–POVMs, mutually unbiased bases (MUBs) and stabilizer states), (global) Pauli observables, Pauli basis measurements, and the uniform/covariant POVM. For details on these measurement classes, and the computation of the LS and PLS estimators we refer to section 4. Here we focus on explaining the novel results and the key techniques used in establishing them. The following theorem shows that the PLS estimator ${\hat{\rho }}_{n}$ converges to the true state ρ in trace distance ||⋅||₁ and provides rank-dependent concentration bounds for each measurement class and arbitrary sample size n.

Theorem 1. (Error bound for ${\hat{\rho }}_{n}$) Let $\rho \in {\mathbb{H}}_{d}$ be a state and fix a number of samples $n\in \mathbb{N}$. Then, for each of the aforementioned measurements, the PLS estimator (table 1) obeys

$$\begin{equation*}\mathrm{Pr}\left[{{\Vert}{\hat{\rho }}_{n}-\rho {\Vert}}_{1}{\geqslant}{\epsilon}\right]{\leqslant}d{\mathrm{e}}^{-\frac{n{{\epsilon}}^{2}}{43g\left(d\right){r}^{2}}}\quad \mathrm{f}\mathrm{o}\mathrm{r}{\epsilon}\in \left[0,1\right],\end{equation*}$$

where $r=\mathrm{min}\left\{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\rho \right),\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\hat{\rho }}_{n}\right)\right\}$ and g(d) specifies dependence on the ambient dimension:

g(d) = 2d for structured POVMs; see equation (8) for ${\hat{L}}_{n}$,

g(d) = d² for Pauli observables; see equation (9) for ${\hat{L}}_{n}$,

g(d) ≃ d^1.6 for Pauli basis measurements; see equation (10) for ${\hat{L}}_{n}$.

The following immediate corollary endows ${\hat{\rho }}_{n}$ with rigorous non-asymptotic error-bars in trace distance.

Corollary 1. The trace-norm ball of size $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\hat{\rho }}_{n}\right)\sqrt{43g\left(d\right){n}^{-1}\mathrm{log}\left(d/\delta \right)}$ around ${\hat{\rho }}_{n}$ (intersected with state space) is a δ-confidence region for the true state ρ.

We emphasize the following aspects of this result:

• (a)  
  (Almost) optimal sampling rate: theorem 1 highlights that

  $$\begin{equation}n{\geqslant}43g\left(d\right)\frac{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}{\left(\rho \right)}^{2}}{{{\epsilon}}^{2}}\mathrm{log}\left(\frac{d}{\delta }\right)\end{equation} \tag{ 5 }$$

  samples suffice to ensure ${{\Vert}{\hat{\rho }}_{n}-\rho {\Vert}}_{1}{\leqslant}{\epsilon}$ with probability at least 1 − δ. For structured POVMs and Pauli observables, this sampling rate is comparable to the best theoretical bounds for alternative tomography algorithms [11, 30]. Moreover, fundamental lower bounds in [13, 20] indicate that this scaling is optimal up to a single log(d)-factor, so it cannot be improved substantially.
• (b)  
  Implicit exploitation of (approximate) low rank: the number of samples required to achieve a good estimator scales quadratically in the rank, rather than the ambient dimension d. This behavior extends to the case where ρ, or ${\hat{\rho }}_{n}$, is well-approximated by a rank-r matrix; see theorem 4 in appendix E. These results are comparable with guarantees for compressed sensing methods [13] that are specifically designed to exploit low-rank. For numerical confirmation we refer to figure 2, and other simulation results in section 5.

Proof sketch for theorem 1. The least-squares estimator ${\hat{L}}_{n}$ can be viewed as a sum of n independent random matrices. To illustrate this, consider a single measurement of the structured POVM type. Then ${\hat{L}}_{1}$ defined in (8) is an instance of the random matrix $X=\left(d+1\right)\vert {v}_{k}\rangle \langle {v}_{k}\vert -\mathbb{I}$, where $k\in \left[m\right]$ occurs with probability $\frac{d}{m}\langle {v}_{k}\vert \rho \vert {v}_{k}\rangle$ (Born's rule). This generalizes to ${\hat{L}}_{n}=\frac{1}{n}{\sum }_{i=1}^{n}{X}_{i}$, where the matrices X_i are statistically independent. Such sums of random matrices concentrate sharply around their expectation value $\mathbb{E}{\hat{L}}_{n}=\rho$, and matrix concentration inequalities [31] quantify this convergence:

$$\begin{equation}\mathrm{Pr}\left[\parallel {\hat{L}}_{n}-\rho {\parallel }_{\infty }{\geqslant}\tau \right]{\leqslant}d{\mathrm{e}}^{-\frac{3n{\tau }^{2}}{8g\left(d\right)}}\quad \tau \in \left[0,1\right].\end{equation} \tag{ 6 }$$

This bound induces a similar operator norm bound for the PLS estimator: the projection is a contraction in operator norm and the shift in eigenvalues is bounded by τ. The resulting bound ${{\Vert}{\hat{\rho }}_{n}-\rho {\Vert}}_{\infty }{\leqslant}2\tau$ may be transformed into a trace-norm bound, as shown in appendix D. This comparison only depends on the (approximate) rank of the states involved. We refer to appendix E for full details of the proof. □

The upper bound in theorem 1 involves a dimension factor d that may be extraneous. This, in turn, introduces an additional log(d)-gap between equation (5) and existing lower bounds [20]. The dimensional factors emerge because we employ matrix-valued concentration inequalities in the proof. Our second main result shows that we can remove this factor for the uniform POVM, which encompasses all rank-one projectors:

Theorem 2. (Convergence of ${\hat{\rho }}_{n}$ for the uniform POVM) For uniform POVM measurements, the PLS estimator obeys

$$\begin{equation*}\mathrm{Pr}\left[\parallel {\hat{\rho }}_{n}-\rho {\parallel }_{1}{\geqslant}{\epsilon}\right]{\leqslant}{\mathrm{e}}^{2.2d-\frac{{{\epsilon}}^{2}n}{480{r}^{2}}}\quad \mathrm{f}\mathrm{o}\mathrm{r}{\epsilon}{ >}0,\end{equation*}$$

where $r=\mathrm{min}\left\{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\rho \right),\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\hat{\rho }}_{n}\right)\right\}$.

This result exactly reproduces the best existing performance guarantees for tomography from independent measurements [22]. The bound follows from standard techniques from high-dimensional probability theory.

Proof sketch of theorem 2. Similarly to the proof of theorm 1, we first establish a (different) concentration bound for the LS estimator ${\hat{L}}_{n}$ with respect to the operator norm, and then use the same arguments as in theorem 1 to convert this into a concentration bound for the PLS with respect to the trace norm. The operator norm has a variational formulation:

$$\begin{equation*}{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }={\mathrm{max}}_{z\in {\mathbb{S}}^{d}}\vert \langle z\vert {\hat{L}}_{n}-\rho \vert z\rangle \vert .\end{equation*}$$

The optimization over the unit sphere may be replaced by a maximization over a finite point set, called a covering net, whose cardinality scales exponentially in d. For any $z\in {\mathbb{S}}^{d}$, $\langle z\vert {\hat{L}}_{n}-\rho \vert z\rangle$ is a sum of n i.i.d. random variables that exhibit subexponential tail decay. (Measuring the uniform POVM allows us to draw this conclusion.) Standard concentration inequalities yield a tail bound that decays exponentially in the number n of samples. Applying a union bound over all points z_i in the net then ensures $\mathrm{Pr}\left[{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\geqslant}\tau \right]{\leqslant}2{\mathrm{e}}^{{c}_{1}d-{c}_{2}n{\tau }^{2}}.$ Subsequently, closeness in operator norm for ${\hat{L}}_{n}$ may be converted into closeness in trace-norm for ${\hat{\rho }}_{n}$ at the cost of an additional (effective) rank factor, cf appendix D. The details of the concentration bound for ${\hat{L}}_{n}$ can be found in appendix F.

□

Tomographically complete measurements can be viewed as injective linear maps $\mathcal{M}:{\mathbb{H}}_{d}\to {\mathbb{R}}^{m}$ with components ${\left[\mathcal{M}\left(X\right)\right]}_{i}=\mathrm{t}\mathrm{r}\left({M}_{i}X\right)$ specified by Born's rule (1). It is well known that the least-squares problem (3) admits the closed-form solution:

$$\begin{equation}{\hat{L}}_{n}={\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left({\mathcal{M}}^{{\dagger}}\;\left({f}_{n}\right)\right).\end{equation} \tag{ 7 }$$

We evaluate this formula for different measurements and content ourselves with sketching key steps and results. The details can be found in appendix B.

4.1.1. Structured POVMs and the uniform POVM

Also known as two-designs, these systems include highly structured, rank-one POVMs ${\left\{\frac{d}{m}\vert {v}_{i}\rangle \langle {v}_{i}\vert \right\}}_{i=1}^{m}$, such as symmetric informationally complete POVMs [32], maximal sets of mutually unbiased bases [33], the set of all stabilizer states [34, 35], as well as the uniform POVM. By definition, for $X\in {\mathbb{H}}_{d}$, all of the above systems obey

$$\begin{equation*}{\mathcal{M}}^{{\dagger}}\mathcal{M}\left(X\right)=\frac{{d}^{2}}{m}\sum _{i=1}^{m}\langle {v}_{i}\vert X\vert {v}_{i}\rangle \vert {v}_{i}\rangle \langle {v}_{i}\vert =\frac{md}{d+1}\left(X+\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}\right).\end{equation*}$$

These equations can readily be inverted, and equation (7) simplifies to

$$\begin{equation}{\hat{L}}_{n}=\left(d+1\right)\sum _{i=1}^{m}{\left[{f}_{n}\right]}_{i}\vert {v}_{i}\rangle \langle {v}_{i}\vert -\mathbb{I}.\end{equation} \tag{ 8 }$$

4.1.2. Pauli observables

Fix d = 2^k (k qubits), and let ${W}_{1},\dots ,\;{W}_{{d}^{2}}\in {\mathbb{H}}_{d}$ be the set of Pauli observables, comprising all possible k-fold tensor products of the elementary 2 × 2 Pauli matrices. We can approximate the expectation value tr(W_i ρ) of each Pauli observable by the empirical mean ${\hat{\mu }}_{i}={\left[{f}_{n/{d}^{2}}^{+}\right]}_{i}-\left[{f}_{n/{d}^{2}}^{-}\right]$ of the two-outcome POVM ${P}_{i}^{{\pm}}=\frac{1}{2}\left(\mathbb{I}{\pm}{W}_{i}\right)$. Pauli matrices form a unitary operator basis, and the evaluation of equation (7) is simple:

$$\begin{equation}{\hat{L}}_{n}=\frac{1}{d}\sum _{i=1}^{{d}^{2}}{\hat{\mu }}_{i}{W}_{i}=\frac{1}{d}\sum _{i=1}^{{d}^{2}}\left({\left[{f}_{n/{d}^{2}}^{+}\right]}_{i}-{\left[{f}_{n/{d}^{2}}^{-}\right]}_{i}\right){W}_{i}.\end{equation} \tag{ 9 }$$

4.1.3. Pauli basis measurements

Rather than approximating (global) expectation values, it is possible to perform different combinations of local Pauli measurements. For d = 2^k , there are 3^k potential combinations in total. Each of the settings $\to {s}\in {\left\{x,y,z\right\}}^{k}$ corresponds to a basis measurement $\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert$, where $\to {o}\in {\left\{{\pm}1\right\}}^{k}$ labels the 2^k potential outcomes. The union $\mathcal{M}$ of all 3^k bases obeys $\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)\left(X\right)={3}^{k}{\mathcal{D}}_{1/3}^{\otimes k}\left(X\right)$, where ${\mathcal{D}}_{1/3}\left(X\right)=\frac{1}{3}\rho +\left(1-\frac{1}{3}\right)\frac{\mathrm{t}\mathrm{r}\left(X\right)}{2}\mathbb{I}$ denotes a single-qubit depolarizing channel. Evaluating equation (7) yields

$$\begin{equation}{\hat{L}}_{n}=\frac{1}{{3}^{k}}\sum _{\to {s},\to {o}}{\left[{f}_{n/{3}^{k}}\right]}_{\to {o}}^{\left(\to {s}\right)}{\left({\mathcal{D}}_{1/3}^{\otimes k}\right)}^{-1}\left(\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \right).\end{equation} \tag{ 10 }$$

Finally, we point out that all of these explicit solutions are guaranteed to have unit trace: $\mathrm{t}\mathrm{r}\left({\hat{L}}_{n}\right)=1$. They are not positive semidefinite in general.

The PLS estimator is defined to be the state closest in Frobenius norm to the least-squares estimator ${\hat{L}}_{n}$. The search (4) admits a simple, analytic solution [27]. Define the all-ones vector $\to {1}\in {\mathbb{R}}^{d}$ and the thresholding function ${\left[\cdot \right]}^{+}$ with components ${\left[\to {y}\right]}_{i}^{+}=\mathrm{max}\left\{{\left[\to {y}\right]}_{i},0\right\}$. Let ${\hat{L}}_{n}=U\mathrm{diag}\left(\to {\lambda }\right){U}^{{\dagger}}$ be an eigenvalue decomposition. Then

$$\begin{equation}{\hat{\rho }}_{n}=U\mathrm{diag}\left({\left[\to {\lambda }-{x}_{0}\to {1}\right]}^{+}\right){U}^{{\dagger}},\end{equation} \tag{ 11 }$$

where ${x}_{0}\in \mathbb{R}$ is chosen so that $\mathrm{t}\mathrm{r}\left({\hat{\rho }}_{n}\right)=1$. The fact that ${\hat{L}}_{n}$ itself has unit trace ensures that this solution to equation (4) is unique. The number x₀ may be determined by applying a root-finding algorithm to the non-increasing function $f\left(x\right)=2+\mathrm{t}\mathrm{r}\left({\hat{L}}_{n}\right)-\;\mathrm{d}x+{\sum }_{i=1}^{d}\left\vert {\lambda }_{i}-x\right\vert$.

The two steps discussed here are inherently scalable: just count frequencies to determine the LS estimators (8)–(10) at a total cost of (at most) $\;\mathrm{min}\left\{m,n\right\}$ matrix additions. The subsequent projection onto state-space is a particular type of soft-thresholding. The associated computational cost is dominated by the eigenvalue decomposition and has runtime (at most) $\mathcal{O}\left({d}^{3}\right)$.

In summary, forming ${\hat{L}}_{n}$ is the dominant cost of a naïve implementation. However, the high degree of structure may allow us to employ techniques from randomized linear algebra [36] to further reduce the cost.

Figure 1 compares ML and PLS for Pauli basis measurements in dimension d = 2⁴. The trace-norm error incurred by PLS is within a factor of two of ML for low-rank states. A more exhaustive simulation study comparing PLS with other 'truncated' estimators for a range of states can be found in [25]. This confirms the statistical efficiency of PLS for the family of low rank states which are the focus of this paper.

Zoom In Zoom Out Reset image size

Computationally, the PLS is significantly faster than ML (implemented using Hradil's algorithm [1, 2]). While an exhaustive comparative analysis of the computational complexity is beyond the scope of the paper, we can get a rough idea by noting that both PLS and ML involve computing a sum of matrices over all outcomes (the LS estimator in the case of PLS and the 'R' operator in the case of ML). For both estimators, this constitutes the most time-consuming subroutine. However, this is done only once in the case of PLS (followed by the projection which is relatively fast), while the ML needs to iterate the step a large number of times to achieve convergence to the maximum. Another advantage of PLS is that the LS part could in principle be computed online while the data is gathered, thus reducing the extra time to that of computing the final projection step.

CS is a natural benchmark for low-rank tomography. The papers [11, 12] apply to Pauli observables, and they show that a random choice of m ⩾ Crdlog⁶(d) Pauli observables is sufficient to reconstruct any rank-r state. The actual reconstruction is performed by solving a convex optimization problem, e.g. the least-squares fit over the set of quantum states [37, 38]. Numerical studies from [13] suggest that m = 256 is appropriate for d = 2⁵ and r = 1. Figure 2 shows that PLS consistently outperforms the CS estimator in this regime. Importantly, PLS was also much faster to evaluate than both, ML and CS.

Zoom In Zoom Out Reset image size

Maximal sets of mutually unbiased bases (MUBs) form a structured POVM (two-design) that lends itself to numerical investigation. Efficient algebraic constructions of MUBs exist for dimensions d = p^k where p is a prime number and k is an arbitrary integer [39–41]. To further underline the implicit advantage of low-rank we fix a prime dimension d and choose a pure state uniformly from the Haar measure on the complex unit sphere in d dimensions. We compute the outcome probabilities for each of the d + 1 different MUB measurements. We then sample outcomes from each distribution a total of $s=\frac{n}{d+1}$ times and compute the estimator ${\hat{\rho }}_{n}$ associated with the total frequency statistics. Figure 3 shows the relation between reconstruction error (in trace distance) and the number of samples per basis on a log–log -scale for a range of different prime dimensions between d = 100 and d = 200. The significant overlap between the graphs for different dimensions suggests that the rate of convergence scales as d/n in terms of total sample size n = s(d + 1) and dimension d; if true, this would mean that the additional log(d)-factor in the norm-one upper bound resulting from theorem 1 may not be necessary, similarly to the upper bound in theorem 2 for uniform measurements. Whether this is the case remains an open question.

Zoom In Zoom Out Reset image size

Corollary 1 is not (yet) optimal. Bootstrapping could be used to obtain tighter confidence regions, and the low computational cost of PLS may speed up this process considerably. It also seems fruitful to combine the ideas presented here with recent insights from [42]. The proof of theorem 1 indicates that PLS is stable with respect to time-dependent state generation (drift). Moreover, the general principle of PLS can be extended to the related problem of process tomography. We intend to address these points in future work.

The uniform/covariant POVM in d-dimensions corresponds to the union of all (properly re-normalized) rank-one projectors: ${\left\{d\vert v\rangle \langle v\vert \mathrm{d}v\right\}}_{v\in \mathbb{S}}^{d}$. Here, dv denotes the unique, unitarily invariant, measure on the complex unit sphere induced by the Haar measure (over the unitary group U(d)). Its high degree of symmetry allows for analyzing this POVM by means of powerful tools from representation theory. This is widely known, see e.g. [43, 44], but we include a short presentation here to be self-contained. Define the frame operator of order k: ${F}_{\left(k\right)}={\int }_{{\mathbb{S}}^{d}}{\left(\vert v\rangle \langle v\vert \right)}^{\otimes k}\mathrm{d}v\in {\mathbb{H}}_{d}^{\otimes k}$. Unitary invariance of dv implies that this frame operator commutes with every k-fold tensor product of a unitary matrix U ∈ U(d):

$$\begin{equation*}{U}^{\otimes k}{F}_{\left(k\right)}={\int }_{{\mathbb{S}}^{d}}{\left(U\vert v\rangle \langle v\vert \right)}^{\otimes k}\mathrm{d}v={\int }_{{\mathbb{S}}^{d}}{\left(\vert \tilde{v}\rangle \langle \tilde{v}\vert U\right)}^{\otimes k}\mathrm{d}\tilde{v}={F}_{\left(k\right)}{U}^{\otimes k}.\end{equation*}$$

Here, we have used a change of variables (v ~ = U v) together with the fact that dv is unitarily invariant (dv ~ = dv). Schur's lemma—one of the most fundamental tools in representation theory—states that any matrix that commutes with every element of a given group representation must be proportional to a sum of the projectors onto the associated irreducible representations (irreps). For the task at hand, the representation of interest is the diagonal representation of the unitary group: U ↦ U^⊗k for all U ∈ U(d). This representation affords, in general, many irreps that may be characterized using Schur–Weyl duality. The symmetric subspace ${\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(k\right)}\subset {\left({\mathbb{C}}^{d}\right)}^{\otimes k}$ is one of them and corresponds to the subspace of all vectors that are invariant under permuting tensor factors. Crucially, F_(k) is an average over rank-one projectors onto vectors |v⟩^⊗k ∈ Sym^(k) and, therefore, its range must be contained entirely within Sym^(k). Combining this with the assertion of Schur's lemma then yields

$$\begin{equation}{F}_{\left(k\right)}={\int }_{{\mathbb{S}}^{d}}{\left(\vert v\rangle \langle v\vert \right)}^{\otimes k}\mathrm{d}v={\left(\genfrac{}{}{0.0pt}{}{d+k-1}{k}\right)}^{-1}{P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(k\right)}}\quad k\in \mathbb{N},\end{equation} \tag{ B.1 }$$

The pre-factor ${\left(\genfrac{}{}{0.0pt}{}{d+k-1}{k}\right)}^{-1}=\mathrm{dim}{\left({\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(k\right)}\right)}^{-1}$ follows from the fact that F_(k) has unit trace.

This closed-form expression is very useful. In particular, it implies that the uniform POVM ${\left\{d\vert v\rangle \langle v\vert \right\}}_{v\in {\mathbb{S}}^{d}}$ is almost an isometry. Fix $X\in {\mathbb{H}}_{d}$ and compute

$$\begin{align}\hfill \left(d+1\right){\int }_{{\mathbb{S}}^{d}}d\langle v\vert X\vert v\rangle \vert v\rangle \langle v\vert \mathrm{d}v=& \left(d+1\right)d{\mathrm{t}\mathrm{r}}_{1}\left({\int }_{{\mathbb{S}}^{d}}{\left(\vert v\rangle \langle v\vert \right)}^{\otimes 2}\mathrm{d}v\left(X\otimes \mathbb{I}\right)\right)\hfill \\ \hfill =& 2{\mathrm{t}\mathrm{r}}_{1}\left({P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(2\right)}}\left(X\otimes \mathbb{I}\right)\right),\hfill \end{align} \tag{ B.2 }$$

where tr₁(A ⊗ B) = tr(A)B denotes the partial trace over the first tensor factor. The projector onto the totally symmetric subspace of two parties has an explicit representation: ${P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(2\right)}}=\frac{1}{2}\left(\mathbb{I}+\mathbb{F}\right)$, where $\mathbb{F}$ denotes the flip operator, i.e. $\mathbb{F}\vert x\rangle \otimes \vert y\rangle =\vert y\rangle \otimes \vert x\rangle$ for all $\vert x\rangle ,\vert y\rangle \in {\mathbb{C}}^{d}$ and extend it linearly to the entire tensor product. Inserting this explicit characterization into equation (B.2) yields

$$\begin{equation}\left(d+1\right){\int }_{{\mathbb{S}}^{d}}d\langle v\vert X\vert v\rangle \vert v\rangle \langle v\vert \mathrm{d}v={\mathrm{t}\mathrm{r}}_{1}\left(\left(\mathbb{I}+\mathbb{F}\right)\left(X\otimes \mathbb{I}\right)\right)=X+\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}.\end{equation} \tag{ B.3 }$$

We emphasize that the full symmetry of the uniform POVM is not required to derive this formula: equation (B.1) for k = 2 is sufficient. This motivates the following definition:

Definition 1 (Two-design).  A (finite) set of m rank-one projectors ${\left\{\vert {v}_{i}\rangle \langle {v}_{i}\vert \right\}}_{i=1}^{m}$ is called a (complex-projective) two-design if

$$\begin{equation*}\frac{1}{m}\sum _{i=1}^{m}{\left(\vert {v}_{i}\rangle \langle {v}_{i}\vert \right)}^{\otimes 2}={\left(\genfrac{}{}{0.0pt}{}{d+1}{2}\right)}^{-1}{P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(2\right)}}.\end{equation*}$$

Each two-design is proportional to a POVM $\mathcal{M}={\left\{\frac{d}{m}\vert {v}_{i}\rangle \langle {v}_{i}\vert \right\}}_{i=1}^{m}$. Clearly, each element is positive semidefinite and $\frac{d}{m}{\sum }_{i=1}^{m}\vert {v}_{i}\rangle \langle {v}_{i}\vert =\mathbb{I}$ follows from taking the partial trace of the defining property. Viewed as maps $\mathcal{M}:{\mathbb{H}}_{d}\to {\mathbb{R}}^{m}$, these POVMs obey

$$\begin{equation*}{\mathcal{M}}^{{\dagger}}\mathcal{M}\left(X\right)=\frac{{d}^{2}}{{m}^{2}}\sum _{i=1}^{m}\langle {v}_{i}\vert X\vert {v}_{i}\rangle \vert {v}_{i}\rangle \langle {v}_{i}\vert =\frac{d\left(X+\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}\right)}{\left(d+1\right)m}\quad \mathrm{f}\mathrm{o}\mathrm{r}X\in {\mathbb{H}}_{d}.\end{equation*}$$

This can be readily inverted

$$\begin{equation}{\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left(X\right)=\frac{m}{d}\left(\left(d+1\right)X-\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}\right).\end{equation} \tag{ B.4 }$$

and the least-squares estimator becomes

$$\begin{align}\hfill {\hat{L}}_{n}=& {\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left({\mathcal{M}}^{{\dagger}}\left(f\right)\right)={\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left(\frac{d}{m}\sum _{i=1}^{m}{f}_{i}\vert {v}_{i}\rangle \langle {v}_{i}\vert \right)\hfill \\ \hfill =& \sum _{i=1}^{m}{f}_{i}\left(\left(d+1\right)\vert {v}_{i}\rangle \langle {v}_{i}\vert -\mathrm{t}\mathrm{r}\left(\vert {v}_{i}\rangle \langle {v}_{i}\vert \right)\mathbb{I}\right)=\left(d+1\right)\sum _{i=1}^{m}{f}_{i}\vert {v}_{i}\rangle \langle {v}_{i}\vert -\mathbb{I}\hfill \end{align} \tag{ B.5 }$$

for any frequency vector ${f}_{n}\in {\mathbb{R}}^{m}$. Mathematically, this is a consequence of the fact that two-design POVMs 'almost' form a tight frame on ${\mathbb{H}}_{d}$. The close connection to well-behaved, tomographically complete, rank-one POVMs has spurred considerable interest in the identification of two-designs. Over the past decades, the following concrete examples have been identified:

• (a)  
  Equiangular lines (SIC POVMs): a family of m unit vectors $\vert {v}_{1}\rangle ,\dots ,\;\vert {v}_{m}\rangle \in {\mathbb{S}}^{d}$ is equiangular, if ${\left\vert \langle {v}_{i},{v}_{j}\rangle \right\vert }^{2}$ is constant for all i ≠ j. The minimal cardinality of such a set is m = d² in which case the angle must be fixed: ${\left\vert \langle {v}_{i},{v}_{j}\rangle \right\vert }^{2}=\frac{1}{d+1}$. Such maximal sets of equiangular lines are known to form two-designs [32] and have been termed symmetric, informationally complete (SIC) POVMs. This nomenclature underlines the importance of equation (B.3) for the original quantum motivation of the study of equiangular lines. While several explicit constructions of SIC POVMs exist, the general question of their existence remains an intriguing open problem.
• (b)  
  Mutually unbiased bases (MUBs): two orthonormal bases ${\left\{\vert {b}_{i}\rangle \right\}}_{i=1}^{d}$ and ${\left\{\vert {c}_{i}\rangle \right\}}_{i=1}^{d}$ of ${\mathbb{C}}^{d}$ are mutually unbiased if ${\left\vert \langle {b}_{i},{c}_{j}\rangle \right\vert }^{2}=\frac{1}{d}$ for all 1 ⩽ i, j ⩽ d. The study of such mutually unbiased bases (MUBs) has a rich history in quantum mechanics that dates back to Schwinger [45]. It is known that at most (d + 1) pairwise mutually unbiased bases can exist in dimension d and explicit algebraic constructions are known for prime power dimensions (d = p^k ). Klappenecker and Roettler [33] showed that maximal sets of MUBs are guaranteed to form two-designs.
• (c)  
  Stabilizer states (STABs): the stabilizer formalism is one of the cornerstones of quantum computation, fault tolerance and error correction, see e.g. [46]. Let ${\mathcal{P}}_{k}$ be the Pauli group on k qubits (d = 2^k ), i.e. the group generated by k-fold tensor products of the elementary Pauli matrices, $\mathbb{I}$ and $i\mathbb{I}$. It is then possible to find maximal abelian subgroups $\mathcal{S}\subset {\mathcal{P}}_{k}$ of size d = 2^k . Since all matrices $W\in \mathcal{S}$ commute, they can be simultaneously diagonalized and determine a single unit vector which is the joint eigenvector with eigenvalue +1 of all the matrices in $\mathcal{S}$ (provided that $-\mathbb{I}\notin \mathcal{S}$). Such vectors are called stabilizer states (STAB) and the group $\mathcal{S}\subset {\mathcal{P}}_{k}$ is its associated stabilizer group. A total of $m={2}^{k}{\prod }_{i=0}^{k}\left({d}^{i}+1\right)={2}^{\frac{1}{2}{k}^{2}+o\left(k\right)}$ different stabilizer states can be generated this way. The union of all of them is actually known to form a three-design [47–49] and, therefore, also a two-design. The latter is also a consequence of earlier results [34, 35]

For d = 2^k , the Pauli matrices ${W}_{1},\dots ,\;{W}_{{d}^{2}}\in {\mathbb{H}}_{d}$ arise from all possible k-fold tensor products of elementary Pauli matrices $\left\{\mathbb{I},{\sigma }_{x},{\sigma }_{y},{\sigma }_{z}\right\}\subset {\mathbb{H}}_{2}$. They are well-known to form a unitary operator basis:

$$\begin{equation}X=\frac{1}{d}\sum _{i=1}^{{d}^{2}}\mathrm{t}\mathrm{r}\left({W}_{i}X\right){W}_{i},\end{equation} \tag{ B.6 }$$

for all $X\in {\mathbb{H}}_{d}$. While they do constitute observables, Pauli matrices by themselves are not POVMs. However, every observable W_i may be associated with a two-outcome POVM ${\mathcal{M}}_{i}=\left\{{P}_{i}^{{\pm}}\right\}=\left\{\frac{1}{2}\left(\mathbb{I}{\pm}{W}_{i}\right)\right\}$. The union ${\bigcup }_{i=1}^{{d}^{2}}{\mathcal{M}}_{i}$ of all these two-outcome POVMs constitutes a linear map $\mathcal{M}:{\mathbb{H}}_{d}\to {\mathbb{R}}^{2\;\mathrm{m}}$ that obeys

$$\begin{equation*}{\mathcal{M}}^{{\dagger}}\mathcal{M}\left(X\right)=\sum _{i=1}^{{d}^{2}}\sum _{o={\pm}}\mathrm{t}\mathrm{r}\left({P}_{i}^{o}X\right){P}_{i}^{o}=\sum _{i=1}^{{d}^{2}}\frac{1}{2}\left(\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}+\mathrm{t}\mathrm{r}\left({W}_{i}X\right){W}_{i}\right)=\frac{1}{2}\left({d}^{2}\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}+X\right),\end{equation*}$$

where the equality follows from equation (B.6). Once more, this expression can be readily inverted:

$$\begin{equation*}{\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left(X\right)=\frac{2}{d}X-\frac{2\mathrm{t}\mathrm{r}\left(X\right)}{{d}^{2}+1}\mathbb{I}.\end{equation*}$$

Before we continue, we note that one Pauli matrix is equal to the identity, say ${W}_{1}=\mathbb{I}$, and the associated POVM is trivial. Hence, we suppose that n copies of ρ are distributed equally among all d² − 1 non-trivial two-Outcome POVMs ${\mathcal{M}}_{i}$. We denote the resulting frequencies by ${\left[f\right]}_{i}^{{\pm}}$ and suppress the dependence on the number of samples. Then, the explicit solution to the least squares problem becomes

$$\begin{align*}\hfill {\hat{L}}_{n}=& {\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left({\mathcal{M}}^{{\dagger}}\left({f}_{n}\right)\right)={\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left(\mathbb{I}\right)+\sum _{i=2}^{{d}^{2}}\sum _{o={\pm}}{\left[f\right]}_{i}^{o}{\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left({P}_{i}^{o}\right)\hfill \\ \hfill =& \frac{2}{d\left({d}^{2}+1\right)}\mathbb{I}+\sum _{i=2}^{{d}^{2}}{\left[f\right]}_{i}^{+}\left(\frac{1}{d}\left(\mathbb{I}+{W}_{i}\right)-\frac{d}{{d}^{2}+1}\mathbb{I}\right)+\sum _{i=2}^{{d}^{2}}{\left[f\right]}_{i}^{-}\left(\frac{1}{d}\left(\mathbb{I}-{W}_{i}\right)-\frac{d}{{d}^{2}+1}\right)\hfill \\ \hfill =& \frac{1}{d}\sum _{i=2}^{{d}^{2}}\left({\left[f\right]}_{i}^{+}-{\left[f\right]}_{i}^{-}\right){W}_{i}+\frac{2+{\sum }_{i=2}^{{d}^{2}}\left({\left[f\right]}_{i}^{+}+{\left[f\right]}_{i}^{-}\right)}{d\left({d}^{2}+1\right)}\mathbb{I}.\hfill \end{align*}$$

We can simplify this expression further by noticing that each two-outcome POVM is dichotomic: either + or − is observed for every run. This implies ${\left[f\right]}_{i}^{+}+{\left[f\right]}_{i}^{-}=1$ and, by extension, ${\sum }_{i=2}^{{d}^{2}}\left({\left[{f}_{i}\right]}_{i}^{+}+{\left[{f}_{i}\right]}_{i}^{-}\right)={d}^{2}-1$. Hence,

$$\begin{equation}{\hat{L}}_{n}=\frac{1}{d}\sum _{i=2}^{{d}^{2}}\left({\left[f\right]}_{i}^{+}-{\left[f\right]}_{i}^{-}\right){W}_{i}+\frac{1}{d}\mathbb{I}=\frac{1}{d}\sum _{i=1}^{{d}^{2}}\left({\left[f\right]}_{i}^{+}-{\left[f\right]}_{i}^{-}\right){W}_{i},\end{equation} \tag{ B.7 }$$

because ${\left[f\right]}_{1}^{+}=1$. This is the formula from the main text and has a compelling interpretation: the difference ${\hat{\mu }}_{i}={\left[{f}_{i}\right]}_{i}^{+}-{\left[{f}_{i}\right]}_{i}^{-}$ is an empirical estimate for the expectation value ${\mu }_{i}=\mathrm{t}\mathrm{r}\left({W}_{i}X\right)$ of the ith Pauli observable. Finally, note that this estimator is again unbiased with respect to random fluctuations in the sample statistics:

$$\begin{equation}\mathbb{E}\left[{\hat{L}}_{n}\right]=\frac{1}{d}\sum _{i=1}^{d}\mathrm{t}\mathrm{r}\left({W}_{i}\rho \right){W}_{i}=\rho .\end{equation} \tag{ B.8 }$$

Before considering the general case, we find it instructive to consider the single qubit case in more detail. For now, fix d = 2 and note that there are three non-trivial Pauli matrices σ_s with $s\in \left\{x,y,z\right\}$. We may associate each σ_s with a two-outcome POVM that is also a basis measurement: $\frac{1}{2}\left(\mathbb{I}{\pm}{\sigma }_{s}\right)=\vert {b}_{{\pm}}^{\left(s\right)}\rangle \langle {b}_{{\pm}}^{\left(s\right)}\vert$. For s ≠ s'

$$\begin{equation*}{\left\vert \langle {b}_{{\pm}}^{\left(s\right)},{b}_{{o}^{\prime }}^{\left({\pm}\right)}\rangle \right\vert }^{2}=\frac{1}{4}\mathrm{t}\mathrm{r}\left(\left(\mathbb{I}{\pm}{\sigma }_{s}\right)\left(\mathbb{I}{\pm}{\sigma }_{{s}^{\prime }}\right)\right)=\frac{1}{2},\end{equation*}$$

because σ_s , σ_s' and σ_s σ_s' have vanishing trace. This implies that the six vectors $\vert {b}_{o}^{\left(s\right)}\rangle \langle {b}_{o}^{\left(s\right)}\vert$ with $o\in \left\{{\pm}1\right\}$ form a maximal set of 3 = (d + 1) mutually unbiased bases. Such vector sets form spherical two-designs and equation (B.3) ensures for any $X\in {\mathbb{H}}_{d}$

$$\begin{equation}\sum _{s,o}\langle {b}_{o}^{\left(s\right)}\vert X\vert {b}_{o}^{\left(s\right)}\rangle \vert {b}_{o}^{\left(s\right)}\rangle \langle {b}_{o}^{\left(s\right)}\vert =X+\mathrm{t}\mathrm{r}\left(X\right)\mathbb{I}=3\left(\frac{1}{3}X+\frac{2}{3}\frac{\mathrm{t}\mathrm{r}\left(X\right)}{2}\mathbb{I}\right)=3{\mathcal{D}}_{1/3}\left(X\right).\end{equation} \tag{ B.9 }$$

Here, ${\mathcal{D}}_{1/3}\left(X\right)=\frac{1}{3}X+\left(1-\frac{1}{3}\right)\frac{\mathrm{t}\mathrm{r}\left(X\right)}{2}\mathbb{I}$ denotes a single-qubit depolarizing channel with loss parameter $p=\frac{1}{3}$.

This behavior extends to multi-qubit systems, i.e. d = 2^k . Suppose that we perform k local (single-qubit) Pauli measurements on a k-qubit state $\rho \in {\mathbb{H}}_{d}$. Then, there are a total of k potential combinations that we label by a string $\to {s}=\left({s}_{1},\dots ,\;{s}_{k}\right)\in {\left\{x,y,z\right\}}^{k}$. Each of them corresponds to a POVM ${\mathcal{M}}^{\left(\to {s}\right)}$ with 2^k = d outcomes that we label by $\to {o}=\left({o}_{1},\dots ,\;{o}_{k}\right)\in {\left\{{\pm}1\right\}}^{k}$. The POVM element associated with index $\to {s}$ and outcome $\to {o}$ has an appealing tensor-product structure: $\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert ={\otimes }_{i=1}^{k}\vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert$. Let $\mathcal{M}={\bigcup }_{\to {s}}{\mathcal{M}}^{\left(\to {s}\right)}:{\mathbb{H}}_{d}\to {\mathbb{R}}^{{3}^{k}}{\times}{\mathbb{R}}^{{2}^{k}}$ denote the union of all such basis measurements. Then, the following formula is true for tensor product matrices $X={\otimes }_{i=1}^{k}{X}_{i}$ and ${X}_{i}\in {\mathbb{H}}_{2}$:

$$\begin{align}\hfill {\mathcal{M}}^{{\dagger}}\mathcal{M}\left(X\right)& =\sum _{\to {s},\to {o}}\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert X\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert =\underset{i=1}{\overset{k}{\otimes }}\left(\sum _{{s}_{i},{o}_{i}}\langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert {X}_{i}\vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert \right)\hfill \\ \hfill & ={3}^{k}\underset{i=1}{\overset{k}{\otimes }}\mathcal{D}\left({X}_{i}\right)={3}^{k}{\mathcal{D}}_{1/3}^{\otimes k}\left(X\right),\hfill \end{align} \tag{ B.10 }$$

where we have used equation (B.9). Linear extension ensures that this formula remains valid for arbitrary matrices $X\in {\mathbb{H}}_{d}$. Since the single qubit depolarizing channel is invertible, the same is true for its k-fold tensor product and we conclude

$$\begin{equation*}{\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left(X\right)=\frac{1}{{3}^{k}}{\left({\mathcal{D}}_{1/3}^{\otimes k}\right)}^{-1}\left(X\right).\end{equation*}$$

Inserting this explicit expression into the closed-form expression for the least squares estimator yields

$$\begin{equation*}{\hat{L}}_{n}={\left({\mathcal{M}}^{{\dagger}}\mathcal{M}\right)}^{-1}\left({\mathcal{M}}^{{\dagger}}\left(\to {f}\right)\right)=\frac{1}{{3}^{k}}\sum _{\to {s},\to {o}}{\left[f\right]}_{\to {o}}^{\left(\to {s}\right)}{\left({\mathcal{D}}_{1/3}^{\otimes k}\right)}^{-1}\left(\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \right),\end{equation*}$$

as advertised in the main text. Here, ${\left[f\right]}_{\to {o}}^{\left(\to {s}\right)}$ is assumed to be a frequency approximation to ${p}_{o}^{\to {\left(s\right)}}=\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \rho \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle$.

We conclude this section with a single-qubit observations that allows for characterizing this expression in a more explicit fashion. Note that one may rewrite ${\mathcal{D}}_{1/3}\left(X\right)$ as $\frac{\mathrm{t}\mathrm{r}\left(X\right)}{2}\mathbb{I}+\frac{1}{6}{\sum }_{s}\mathrm{t}\mathrm{r}\left({\sigma }_{s}X\right){\sigma }_{s}$. This facilitates the computation of the single-qubit inverse:

$$\begin{equation*}{\mathcal{D}}_{1/3}^{-1}\left(X\right)=\frac{\mathrm{t}\mathrm{r}\left(X\right)}{2}\mathbb{I}+\frac{3}{2}\sum _{s}\mathrm{t}\mathrm{r}\left({\sigma }_{s}X\right){\sigma }_{s}\end{equation*}$$

and, in particular

$$\begin{equation*}{\mathcal{D}}_{1/3}^{-1}\left(\vert {b}_{{\pm}}^{\left(s\right)}\rangle \langle {b}_{{\pm}}^{\left(s\right)}\vert \right)=\frac{1}{2}\left({\mathcal{D}}^{-1}\left(\mathbb{I}\right){\pm}{\mathcal{D}}^{-1}\left({\sigma }_{s}\right)\right)=\frac{1}{2}\left(\mathbb{I}{\pm}3{\sigma }_{s}\right)=3\vert {b}_{{\pm}}^{\left(s\right)}\rangle \langle {b}_{{\pm}}^{\left(s\right)}\vert -\mathbb{I}.\end{equation*}$$

This ensures

$$\begin{equation}{\hat{L}}_{n}=\frac{1}{{3}^{k}}\sum _{\to {s},\to {o}}{\left[f\right]}_{\to {o}}^{\left(\to {s}\right)}\underset{i=1}{\overset{k}{\otimes }}\left(3\vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert -\mathbb{I}\right).\end{equation} \tag{ B.11 }$$

which, again, is an unbiased estimator with respect to random fluctuations in the sample statistics.

Finally, we also point out another consequence that will be important later on:

$$\begin{equation}{\left({\mathcal{D}}_{1/3}^{-1}\left(\vert {b}_{o}^{\left(s\right)}\rangle \langle {b}_{o}^{\left(s\right)}\vert \right)\right)}^{2}=5{\mathcal{D}}_{3/5}\left(\vert {b}_{o}^{\left(s\right)}\rangle \langle {b}_{o}^{\left(s\right)}\vert \right),\end{equation} \tag{ B.12 }$$

where ${\mathcal{D}}_{3/5}$ is another single-qubit depolarizing channel.

For structured measurements (two-designs) we may rewrite the (plain) least squares estimator (B.5) as

$$\begin{equation*}{\hat{L}}_{n}=\left(d+1\right)\sum _{i=1}^{m}{\left[{f}_{n}\right]}_{i}\vert {v}_{i}\rangle \langle {v}_{i}\vert -\mathbb{I}=\frac{1}{n}\sum _{i=1}^{n}{X}_{i},\end{equation*}$$

where each X_i is an i.i.d. copy of the random matrix $X\in {\mathbb{H}}_{d}$ that assumes $\left(d+1\right)\vert {v}_{k}\rangle \langle {v}_{k}\vert -\mathbb{I}$ with probability $\frac{d}{m}\langle {v}_{k}\vert \rho \vert {v}_{k}\rangle$ for all $k\in \left[m\right]$. Unbiasedness with respect to the sample statistics ensures $\mathbb{E}\left[{\hat{L}}_{n}\right]=\mathbb{E}\left[X\right]=\rho$. Hence, ${\hat{L}}_{n}-\rho$ is a sum of i.i.d., centered random matrices $\frac{1}{n}\left({X}_{i}-\mathbb{E}\left[{X}_{i}\right]\right)$. These obey

$$\begin{equation*}\frac{1}{n}{{\Vert}{X}_{i}-\mathbb{E}\left[{X}_{i}\right]{\Vert}}_{\infty }=\frac{1}{n}{{\Vert}\left(d+1\right)\vert {v}_{k}\rangle \langle {v}_{k}\vert -\mathbb{I}-\rho {\Vert}}_{\infty }{\leqslant}\frac{d}{n}=:R,\end{equation*}$$

where $k\in \left[m\right]$ is arbitrary. Next, note that the random matrix X obeys

$$\begin{equation*}\mathbb{E}\left[{\left(X-\mathbb{E}\left[X\right]\right)}^{2}\right]=\mathbb{E}\left[{X}^{2}\right]-\mathbb{E}{\left[X\right]}^{2}=\mathbb{E}\left[{X}^{2}\right]-{\rho }^{2}\end{equation*}$$

and also

$$\begin{align*}\hfill \mathbb{E}\left[{X}^{2}\right]=& \sum _{k=1}^{m}\frac{d}{m}\langle {v}_{k}\vert \rho \vert {v}_{k}\rangle {\left(\left(d+1\right)\vert {v}_{k}\rangle \langle {v}_{k}\vert -\mathbb{I}\right)}^{2}=\frac{d\left({d}^{2}-1\right)}{m}\sum _{k=1}^{m}\langle {v}_{k}\vert \rho \vert {v}_{k}\rangle \vert {v}_{k}\rangle \langle {v}_{k}\vert +\mathbb{I}\hfill \\ \hfill =& \left(d-1\right)\left(\rho +\mathbb{I}\right)+\mathbb{I},\hfill \end{align*}$$

according to equation (B.3). This allows us to bound the variance parameter:

$$\begin{equation*}{{\Vert}\sum _{i=1}^{n}\mathbb{E}{\left[\frac{1}{n}\left(X-\mathbb{E}\left[X\right]\right)\right]}^{2}{\Vert}}_{\infty }=\frac{1}{n}{{\Vert}\left(d-1\right)\rho +d\mathbb{I}-{\rho }^{2}{\Vert}}_{\infty }{\leqslant}\frac{2d}{n}=:{\sigma }^{2}.\end{equation*}$$

The ratio $\frac{{\sigma }^{2}}{R}=2$ indicates that any choice of $\tau \in \left[0,2\right]$ will fall into the subgaussian regime of the matrix Bernstein inequality and theorem 3 yields

$$\begin{equation}\mathrm{Pr}\left[{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\geqslant}\tau \right]=\mathrm{Pr}\left[{{\Vert}\frac{1}{n}\sum _{i=1}^{n}\left({X}_{i}-\mathbb{E}\left[{X}_{i}\right]\right){\Vert}}_{\infty }{\geqslant}\tau \right]{\leqslant}d{\mathrm{e}}^{-\frac{3{\tau }^{2}n}{16d}}.\end{equation} \tag{ C.1 }$$

We assume that the total number of samples n is distributed equally among the d² different Pauli measurements. Similar to before, unbiasedness (B.8) and the explicit characterization of the LI estimator (B.7) allow us to write

$$\begin{equation*}{\hat{L}}_{n}-\rho =\sum _{k=1}^{{d}^{2}}\frac{1}{n}\sum _{i=1}^{n/{d}^{2}}\left({X}_{i}^{\left(k\right)}-\mathbb{E}\left[{X}_{i}^{\left(k\right)}\right]\right).\end{equation*}$$

Here, each ${X}_{i}^{\left(k\right)}$ is an independent instance of the random matrix X^(k) = ±dW_k with probability $\frac{1}{2}\left(1{\pm}\mathrm{t}\mathrm{r}\left({W}_{k}\rho \right)\right)$ each. This is a sum of centered random matrices that are independent, but in general not identically distributed. However, independence alone suffices for applying theorem 3. We note in passing that this would not be the case for earlier (weaker) versions of the matrix Bernstein inequality. Bound

$$\begin{equation*}\frac{1}{n}{{\Vert}{X}_{i}^{\left(k\right)}-\mathbb{E}\left[{X}_{i}^{\left(k\right)}\right]{\Vert}}_{\infty }=\frac{d}{n}{{\Vert}\left(1{\pm}\mathrm{t}\mathrm{r}\left({W}_{i}\rho \right)\right)\mathbb{I}{\Vert}}_{\infty }{\leqslant}\frac{2d}{n}=:R\end{equation*}$$

and use $\mathbb{E}\left[{\left(X-\mathbb{E}\left[X\right]\right)}^{2}\right]{\leqslant}\mathbb{E}\left[{X}^{2}\right]$ (in the positive semidefinite order) to considerably simplify the variance computation:

$$\begin{equation*}\frac{{d}^{2}}{n}{{\Vert}\sum _{k=1}^{{d}^{2}}\mathbb{E}\left[{\left({X}^{\left(k\right)}\right)}^{2}\right]{\Vert}}_{\infty }=\frac{1}{n}{{\Vert}\sum _{k=1}^{{d}^{2}}\mathbb{I}{\Vert}}_{\infty }=\frac{{d}^{2}}{n}=:{\sigma }^{2},\end{equation*}$$

because ${\left({X}^{\left(k\right)}\right)}^{2}=\frac{1}{{d}^{2}}\mathbb{I}$. Applying theorem 3 yields

$$\begin{equation*}\mathrm{Pr}\left[{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\geqslant}\tau \right]{\leqslant}d{\mathrm{e}}^{-\frac{3{\tau }^{2}\tilde{n}}{8}}\quad \tau \in \left[0,d/2\right].\end{equation*}$$

Once more, we assume that the total budget of samples n is distributed equally among all 3^k Pauli basis choices. Unbiasedness of the LI estimator together with the explicit description (B.11) allows us to once more interpret ${\hat{L}}_{n}-\rho$ as a sum of independent, centered random matrices:

$$\begin{equation*}{\hat{L}}_{n}-\rho =\sum _{\to {s}}\frac{1}{n}\sum _{i=1}^{n/{3}^{k}}\left({X}_{i}^{\left(\to {s}\right)}-\mathbb{E}\left[{X}_{i}^{\left(\to {s}\right)}\right]\right).\end{equation*}$$

For each $\to {s}\in {\left\{x,y,z\right\}}^{k}$, ${X}_{i}^{\left(\to {s}\right)}$ is an independent copy of the random matrix

$$\begin{equation*}{X}^{\left(\to {s}\right)}=\underset{i=1}{\overset{k}{\otimes }}\left(3\vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert -\mathbb{I}\right)\end{equation*}$$

with probability $\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \rho \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle$ for each $\to {o}\in {\left\{{\pm}1\right\}}^{k}$. Jensen's inequality implies

$$\begin{equation*}\frac{1}{n}{{\Vert}{X}^{\left(\to {s}\right)}-\mathbb{E}\left[{X}^{\left(\to {s}\right)}\right]{\Vert}}_{\infty }{\leqslant}\frac{2}{n}{{\Vert}{X}^{\left(\to {s}\right)}{\Vert}}_{\infty }=\frac{2}{n}\prod _{i=1}^{k}{{\Vert}3\vert {b}_{{o}_{i}}^{\left({s}_{i}\right)}\rangle \langle {b}_{{o}_{i}}^{\left({s}_{i}\right)}\vert -\mathbb{I}{\Vert}}_{\infty }=\frac{{2}^{k+1}}{n}=:R.\end{equation*}$$

For the variance, we once more use $\mathbb{E}\left[{\left(X-\mathbb{E}\left[X\right]\right)}^{2}\right]{\leqslant}\mathbb{E}\left[{X}^{2}\right]$ and compute

$$\begin{align*}\hfill \frac{1}{n}\sum _{\to {s}}\mathbb{E}{\left({X}^{\left(\to {s}\right)}\right)}^{2}=& \sum _{\to {s},\to {o}}\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \rho \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \underset{i=1}{\overset{k}{\otimes }}{\left({\mathcal{D}}^{-1}\left(\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \right)\right)}^{2}\hfill \\ \hfill =& \sum _{\to {s}}\sum _{\to {o}}\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \rho \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \underset{i=1}{\overset{k}{\otimes }}\frac{5}{3}{\mathcal{D}}_{3/5}\left(\vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \right)\hfill \\ \hfill =& {5}^{k}{\mathcal{D}}_{3/5}^{\otimes k}\left(\frac{1}{{3}^{k}}\sum _{\to {s},\to {o}}\langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \rho \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \vert {b}_{\to {o}}^{\left(\to {s}\right)}\rangle \langle {b}_{\to {o}}^{\left(\to {s}\right)}\vert \right)\hfill \\ \hfill =& {5}^{k}{\mathcal{D}}_{3/5}^{\otimes k}\left({\mathcal{D}}_{1/3}^{\otimes k}\left(\rho \right)\right)=\frac{{5}^{k}}{n}{\mathcal{D}}_{1/5}^{\otimes k}\left(\rho \right),\hfill \end{align*}$$

where we have used equation (B.12) and the fact that the combination of two depolarizing channels is again a depolarizing channel. This expression can be evaluated explicitly. For $\alpha \subset \left[k\right]$, let tr_α (ρ) denote the partial trace over all indices contained in α. Then, for $X={\otimes }_{j=1}^{k}{X}_{j}\in {\mathbb{H}}_{d}$

$$\begin{equation*}{5}^{k}{\mathcal{D}}_{1/5}^{\otimes k}\left({\otimes }_{j=1}^{k}{X}_{j}\right)=\underset{j=1}{\overset{k}{\otimes }}\left(\mathrm{t}\mathrm{r}\left({X}_{j}\right)\mathbb{I}+{X}_{j}\right)=\sum _{\alpha \subset \left[k\right]}{2}^{\vert \alpha \vert }{\mathrm{t}\mathrm{r}}_{\alpha }\left(X\right)\otimes {\mathbb{I}}^{\otimes \alpha }\end{equation*}$$

and this extends linearly to all of ${\mathbb{H}}_{d}\simeq {\mathbb{H}}_{2}^{\otimes k}$. Consequently,

$$\begin{align*}\hfill {{\Vert}\frac{1}{n}\sum _{\to {s}}\mathbb{E}\left[{\left({X}^{\left(\to {s}\right)}\right)}^{2}\right]{\Vert}}_{\infty }=& \frac{1}{n}{{\Vert}\sum _{\alpha \subset \left[k\right]}{2}^{\vert \alpha \vert }{\mathrm{t}\mathrm{r}}_{\alpha }\left(\rho \right)\otimes {\mathbb{I}}^{\otimes \alpha }{\Vert}}_{\infty }{\leqslant}\frac{1}{n}\sum _{\alpha \subset \left[k\right]}{2}^{\vert \alpha \vert }{{\Vert}{\mathrm{t}\mathrm{r}}_{\alpha }\left(\rho \right){\Vert}}_{\infty }{{\Vert}{\mathbb{I}}^{\otimes \alpha }{\Vert}}_{\infty }\hfill \\ \hfill {\leqslant}& \frac{1}{n}\sum _{\alpha \subset \left[k\right]}{2}^{\vert \alpha \vert }=\frac{1}{n}\sum _{j=0}^{k}\left(\genfrac{}{}{0.0pt}{}{k}{j}\right){2}^{j}=\frac{{\left(2+1\right)}^{k}}{n}=\frac{{3}^{k}}{n}=:{\sigma }^{2}.\hfill \end{align*}$$

This estimate is actually tight for pure product states of the form ρ = (|ψ⟩⟨ψ|)^⊗k . We may now apply theorem 3 to conclude

$$\begin{equation*}\mathrm{Pr}\left[{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\geqslant}\tau \right]{\leqslant}d\mathrm{exp}\left(-\frac{3n{\tau }^{2}}{8{\times}{3}^{k}}\right)\quad \tau \in \left[0,1\right].\end{equation*}$$

Suppose that we perform n uniform POVM measurements on a fixed quantum state $\rho \in {\mathbb{H}}_{d}$. Then, the least squares estimator is equivalent to a sum of i.i.d. random matrices:

$$\begin{equation*}{\hat{L}}_{n}=\frac{1}{n}\sum _{i=1}^{n}{X}_{i}.\end{equation*}$$

Each X_i is an independent copy of the random matrix X that assumes the value $\left(d+1\right)\vert v\rangle \langle v\vert -\mathbb{I}$ with probability d⟨v|ρ|v⟩dv and v may range over the entire complex unit sphere. Unbiasedness of this estimator in turn implies

$$\begin{equation*}{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }={{\Vert}\frac{1}{n}\sum _{i=1}^{n}\left({X}_{i}-\mathbb{E}\left[{X}_{i}\right]\right){\Vert}}_{\infty }={\mathrm{max}}_{y\in {\mathbb{S}}^{d}}\left\vert \langle y\vert \frac{1}{n}\sum _{i=1}^{n}\left({X}_{i}-\mathbb{E}\left[{X}_{i}\right]\right)\vert y\rangle \right\vert .\end{equation*}$$

Next, we employ a result that is somewhat folklore in random matrix theory, see e.g. [61, lemma 5.3]. It states that the maximum over the entire unit sphere may be replaced by a maximum over certain finite point sets, called covering nets: a covering-net of ${\mathbb{S}}^{d}$ with fineness θ > 0 is a finite set of unit vectors ${\left\{{z}_{j}\right\}}_{j=1}^{N}\subseteq {\mathbb{S}}^{d}$ that covers the entire (complex) unit sphere in the sense that every $y\in {\mathbb{S}}^{d}$ is at least θ-close to a point in the net.

Lemma 3. Let ${\mathcal{N}}_{\theta }={\left\{{z}_{j}\right\}}_{j=1}^{N}$ be a covering net of ${\mathbb{S}}^{d}$ with fineness θ. Then, for any matrix $A\in {\mathbb{H}}_{d}$:

$$\begin{equation*}{\mathrm{max}}_{j\in \left[N\right]}\left\vert \langle {z}_{j}\vert A\vert {z}_{j}\rangle \right\vert {\leqslant}{{\Vert}A{\Vert}}_{\infty }{\leqslant}\frac{1}{1-2\theta }{\mathrm{max}}_{j\in \left[N\right]}\left\vert \langle {z}_{j}\vert A\vert {z}_{j}\rangle \right\vert .\end{equation*}$$

This result highlights that already a rather coarse net suffices to get reasonable approximations to the operator norm. Here, we choose $\theta =\frac{1}{4}$ which, while certainly not optimal, simplifies exposition. In particular,

$$\begin{equation}{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\leqslant}2\;{\mathrm{max}}_{j\in \left[N\right]}\left\vert \frac{1}{n}\sum _{j=1}^{n}\langle {z}_{j}\vert {X}_{i}-\mathbb{E}\left[{X}_{i}\right]\vert {z}_{j}\rangle \right\vert ,\end{equation} \tag{ F.2 }$$

where the maximization is over a covering net of fineness $\theta =\frac{1}{4}$.

Note that the right hand side of equation (F.2) corresponds to a maximum over N different random variables—each of them labeled by a unit vector z_j in the net. Let $z\in {\mathbb{S}}^{d}$ be such a vector. Then, the associated random variable itself corresponds to an empirical average of n i.i.d. variables:

$$\begin{equation*}{s}_{z}=\langle z\vert X-\mathbb{E}\left[X\right]\vert z\rangle .\end{equation*}$$

Clearly, s_z obeys $\mathbb{E}\left[{s}_{z}\right]=0$ and, more importantly, has sub-exponential moment growth. While this follows directly from the fact that s_z is bounded, the following result highlights that this tail-behavior is actually independent of the ambient dimension.

Lemma 4. Fix $z\in {\mathbb{S}}^{d}$. Then for any integer p ⩾ 2, the random variable s_z obeys

$$\begin{equation*}\mathbb{E}\left[\vert {s}_{z}{\vert }^{p}\right]{\leqslant}27{\times}{6}^{p-2}p!\end{equation*}$$

We divert the proof of this statement to the end of this section and content ourselves with emphasizing that the closed form expression of the frame operator (B.1) is essential for bounding all moments simultaneously. More relevant to the task at hand is that such a moment behavior ensures that the tails of the distribution of s_z follow an exponential decay: $\mathrm{Pr}\left[\vert {s}_{z}\vert {\geqslant}t\right]{\leqslant}{\mathrm{e}}^{-ct}$, where c is a constant independent of the dimension d. Strong classical concentration inequalities apply for sums of i.i.d. random variables that exhibit such sub-exponential behavior. We choose to apply a rather general version of the classical Bernstein inequality, see e.g. [58, theorem 7.30].

Theorem 6. Let ${s}_{1},\dots ,\;{s}_{n}\in \mathbb{R}$ i.i.d. copies of a mean-zero random variable s that obeys $\mathbb{E}\left[\vert s{\vert }^{p}\right]{\leqslant}p!{R}^{p-2}{\sigma }^{2}/2$ for all integers p ⩾ 2, where R, σ² > 0 are constants. Then, for all t > 0,

$$\begin{equation*}\mathrm{Pr}\left[\left\vert \sum _{i=1}^{n}{s}_{i}\right\vert {\geqslant}t\right]{\leqslant}2\;\mathrm{exp}\left(-\frac{{t}^{2}/2}{n{\sigma }^{2}+Rt}\right).\end{equation*}$$

Lemma 4 ensures that the random variable s_z meets this requirement with σ² = 54 and R = 6. Hence, the following corollary is an immediate consequence of theorem 6.

Corollary 3. Fix $z\in {\mathbb{S}}^{d}$. Then, for any t ∈ [0, 1]

$$\begin{equation*}\mathrm{Pr}\left[\left\vert \frac{1}{n}\sum _{i=1}^{n}\langle z\vert {X}_{i}-\mathbb{E}\left[{X}_{i}\right]\vert z\rangle \right\vert {\geqslant}t\right]{\leqslant}2{\mathrm{e}}^{-\frac{n{t}^{2}}{120}}.\end{equation*}$$

Recall that equation (F.2) upper-bounds ${{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }$ by a maximum over finitely many random variables, each of which is controlled by the strong exponential tail inequality from corollary 3. To exploit this, we fix $\tau \in \left[0,1\right]$ and apply a union bound (also known as Boole's inequality) over all these different random variables to obtain

$$\begin{align*}\hfill \mathrm{Pr}\left[\parallel {\hat{L}}_{n}-\rho {\parallel }_{\infty }{\geqslant}\tau \right]{\leqslant}& \mathrm{Pr}\left[{\mathrm{max}}_{j\in \left[N\right]}\left\vert \frac{1}{n}\sum _{j=1}^{n}\langle {z}_{j}\vert {X}_{i}-\mathbb{E}\left[{X}_{i}\right]\vert {z}_{j}\rangle \right\vert {\geqslant}\frac{\tau }{2}\right]\hfill \\ \hfill {\leqslant}& N{\mathrm{max}}_{j\in \left[N\right]}\mathrm{Pr}\left[\left\vert \frac{1}{n}\sum _{j=1}^{n}\langle {z}_{j}\vert {X}_{i}-\mathbb{E}\left[{X}_{i}\right]\vert {z}_{j}\rangle \right\vert {\geqslant}\frac{\tau }{2}\right]{\leqslant}2N{\mathrm{e}}^{-\frac{n{\tau }^{2}}{480}},\hfill \end{align*}$$

where the last line is due to corollary 3 Here, $N=\vert {\mathcal{N}}_{\frac{1}{4}}\vert$ denotes the cardinality of a covering net for the complex unit sphere ${\mathbb{S}}^{d}$ with fineness $\theta =\frac{1}{4}$. The complex unit sphere admits an isometric embedding into the real-valued unit sphere in 2d-dimensions: map real- and imaginary parts of each complex vector component onto two distinct real parameters. This map preserves Euclidean lengths and, by extension, also the geometry of the unit sphere. Volumetric upper bounds on the cardinality of covering nets for the 2d-dimensional real-valued unit sphere are widely known, see e.g. [58, proposition C.3] and [61, lemma 5.2]: $\vert {\mathcal{N}}_{\theta }\vert {\leqslant}{\left(1+\frac{2}{\theta }\right)}^{2d}$. Since a fineness of $\theta =\frac{1}{4}$ suffices for our purpose, we can conclude N ⩽ 3^2d and consequently,

$$\begin{equation*}\mathrm{Pr}\left[{{\Vert}{\hat{L}}_{n}-\rho {\Vert}}_{\infty }{\geqslant}\tau \right]{\leqslant}2{\times}{3}^{2d}{\mathrm{e}}^{-\frac{n{\tau }^{2}}{480}}=2{\mathrm{e}}^{2\mathrm{log}\left(3\right)d-\frac{n{\tau }^{2}}{480}}\end{equation*}$$

This concludes the proof of theorem 5.

Recall that, by assumption, the random matrix X assumes the value $X=\left(d+1\right)\vert v\rangle \langle v\vert -\mathbb{I}$ with probability ⟨v|ρ|v⟩dv, where v may range over the entire complex unit sphere ${\mathbb{S}}^{d}$. Moreover, $\mathbb{E}\left[X\right]=\rho$. For fixed $z\in {\mathbb{S}}^{d}$, we may therefore write

$$\begin{equation*}{s}_{z}=\langle z\vert X-\mathbb{E}\left[X\right]\vert z\rangle =\left(d+1\right)\langle v\vert B\vert v\rangle ,\end{equation*}$$

where $B=\vert z\rangle \langle z\vert -\frac{1+\langle z\vert \rho \vert z\rangle }{d+1}\mathbb{I}\in {\mathbb{H}}_{d}$ has bounded trace norm

$$\begin{equation}{{\Vert}B{\Vert}}_{1}{\leqslant}1+\left(1+\langle z\vert \rho \vert z\rangle \right){\leqslant}3.\end{equation} \tag{ F.3 }$$

Next, recall a basic identity from matrix analysis that states

$$\begin{equation*}\vert \langle v\vert B\vert v\rangle \vert =\left\vert \mathrm{t}\mathrm{r}\left(\vert v\rangle \langle v\vert B\right)\right\vert {\leqslant}\mathrm{t}\mathrm{r}\left(\vert v\rangle \langle v\vert \vert B\vert \right),\end{equation*}$$

where $\vert B\vert =\sqrt{{B}^{2}}$ denotes the absolute value of the matrix B. Also, the Schatten-p norms of matrices and their absolute values coincides, in particular ||B||₁ = tr(|B|) = |||B|||₁. We can use this trick to absorb the absolute value in the moment computation. More precisely, fix an integer p ⩾ 2 and note that $\mathbb{E}\left[\vert {s}_{z}{\vert }^{p}\right]$ obeys

$$\begin{equation*}\mathbb{E}\left[{\left\vert \left(d+1\right)\langle v\vert B\vert v\rangle \right\vert }^{p}\right]{\leqslant}{\left(d+1\right)}^{p}\mathbb{E}\left[\mathrm{t}\mathrm{r}{\left(\vert v\rangle \langle v\vert \;\vert B\vert \right)}^{p}\right].\end{equation*}$$

We can now include the distribution of the random matrices X, and (by extension) |v⟩⟨v|, to compute

$$\begin{align*}\hfill \mathbb{E}\left[\vert {s}_{z}{\vert }^{p}\right]{\leqslant}& {\left(d+1\right)}^{p}\mathbb{E}\left[\mathrm{t}\mathrm{r}{\left(\vert v\rangle \langle v\vert \;\vert B\vert \right)}^{p}\right]=d{\left(d+1\right)}^{p}{\int }_{{\mathbb{S}}^{d}}\langle v\vert \rho \vert v\rangle \mathrm{t}\mathrm{r}{\left(\vert v\rangle \langle v\vert \;\vert B\vert \right)}^{p}\;\mathrm{d}v\hfill \\ \hfill =& d{\left(d+1\right)}^{p}\mathrm{t}\mathrm{r}\left({\int }_{{\mathbb{S}}^{d}}{\left.\left(\vert v\rangle \langle v\vert \right)\right)}^{\otimes \left(p+1\right)}\;\rho \otimes \vert B{\vert }^{\otimes p}\right)\hfill \\ \hfill =& d{\left(d+1\right)}^{p}{\left(\genfrac{}{}{0.0pt}{}{d+p}{p+1}\right)}^{-1}\mathrm{t}\mathrm{r}\left({P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(p+1\right)}}\rho \otimes \vert B{\vert }^{\otimes p}\right),\hfill \end{align*}$$

where the last equation is due to equation (B.1). Next, we note that Hoelder's inequality implies

$$\begin{equation*}\mathrm{t}\mathrm{r}\left({P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(p+1\right)}}\rho \otimes \vert B{\vert }^{\otimes p}\right){\leqslant}{{\Vert}{P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{p+1}}{\Vert}}_{\infty }{{\Vert}\rho {\Vert}}_{1}{\Vert}\;\vert B{\vert \;{\Vert}}_{1}^{p}{\leqslant}{3}^{p},\end{equation*}$$

because ${P}_{{\mathrm{S}\mathrm{y}\mathrm{m}}^{\left(p+1\right)}}$ is an orthogonal projector, ρ is a quantum state and B is bounded in trace norm (F.3). For the remaining pre-factor we use the crude bound

$$\begin{equation*}d{\left(d+1\right)}^{p}{\left(\genfrac{}{}{0.0pt}{}{d+p}{p+1}\right)}^{-1}{\leqslant}\left(p+1\right)!{\leqslant}3{\times}{2}^{p-2}p!\end{equation*}$$

to establish the statement.