1 Introduction↩︎

Standard shadow theory is usually organized around expectation-value prediction, shadow norms, marginal variances, and many-observable prediction guarantees. These quantities control individual reconstructed observables, or collections of such observables, but they do not treat the covariance matrix of the reconstructed shadow-output vector as the primary object. The present paper develops this covariance-matrix viewpoint. For a general shadow-output vector \(\hat{x}\), we consider the covariance matrix \(\Sigma(\rho)=\operatorname{Cov}_\rho(\hat{x})\). In the finite-sample setting, the main object is its fixed selected compression \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\), where \(R_l\) selects a deterministic set of output coordinates before the data are observed. This selected covariance matrix describes the joint fluctuations of reconstructed shadow coordinates, not only their individual means or marginal variances.

This viewpoint is natural because a covariance matrix describes how different coordinates fluctuate together. In multivariate analysis, the mean vector and the covariance matrix are both basic objects: the mean vector gives the first moments of the coordinates, while the covariance matrix describes their second-order joint fluctuations [1], [2]. The covariance matrix also identifies linear combinations of coordinates with large or small variance, as in principal component analysis [3]. For shadow data, the covariance matrix captures the joint fluctuations of reconstructed output coordinates across repeated measurement shots, beyond what is visible from individual expectation estimates or marginal variances alone.

The covariance matrix studied here is not a new physical observable. Its entries are covariances of random reconstructed output coordinates produced by the same measurement-and-reconstruction procedure. Thus off-diagonal entries describe statistical couplings in the shadow post-processing noise. In the Pauli case, for example, \(\operatorname{Cov}_\rho(\hat{x}_P,\hat{x}_Q)\) should not be confused with the physical expectation of a Pauli product such as \(PQ\).

Once the covariance matrix is taken as the object of interest, the finite-sample problem becomes matrix-valued. The empirical covariance matrix must be compared with \(\Sigma_l(\rho)\) in operator norm, and the eigenvalues and spectral projectors of \(\Sigma_l(\rho)\) must be stable under this perturbation. We therefore combine matrix Bernstein concentration, an exact rank-one identity relating true-centered and sample-centered empirical covariances, and standard Weyl and Davis–Kahan perturbation bounds [4]–[9].

The first result is a general finite-sample theorem for fixed selected coordinates. It applies to arbitrary shadow protocols and gives an operator-norm error bound for the selected sample-centered empirical covariance. If the selected number of coordinates and the protocol-dependent constants in this bound are independent of the ambient system size, then the required number of samples is also independent of the ambient dimension. The same operator-norm control gives eigenvalue and spectral-projector guarantees for the selected covariance spectrum.

We then show that the required boundedness condition is satisfied in local measurement settings. For general local product shadow protocols on tensor factors of fixed local dimension, finite-weight product observables depend only on the reconstructed local snapshots on their supports. Consequently, for bounded selected set size, bounded observable weight, and uniformly bounded local reconstruction coefficients, the selected covariance matrix can be estimated with sample complexity independent of the total number of tensor factors. This gives a local-product mechanism for dimension-independent selected covariance estimation, and it is not specific to Pauli measurements.

As an explicit example, we analyze biased local Pauli shadows, introduced for Hamiltonian-estimation tasks by Hadfield–Bravyi–Raymond–Mezzacapo [10]. We use the same measurement mechanism for a different object: the covariance matrix of selected reconstructed Pauli coefficients. In this setting, bounded-weight selected Pauli strings and local basis probabilities bounded away from zero give dimension-independent finite-sample covariance estimation.

For biased local Pauli shadows, we also derive an exact covariance formula. The formula reorganizes the local-Pauli second-moment mechanism already present in the observable-wise variance analysis of Huang–Kueng–Preskill [11] and extends it to biased designs. Pauli compatibility determines which off-diagonal raw second moments can be nonzero, while inverse-probability overlap factors determine their size. Subtracting the corresponding products of means gives the covariance entries. This formula shows how local measurement bias changes both diagonal variances and off-diagonal statistical couplings.

Finally, we contrast these local mechanisms with global Clifford shadows. The general theorem applies to any shadow protocol, but dimension-independent covariance estimation depends on the protocol. Local product protocols can keep the relevant constants in the finite-sample bound independent of the number of tensor factors for finite-weight selected observables. Global Clifford shadows behave differently: the inverse shadow channel amplifies every non-identity Pauli coordinate by the global factor \(d+1\), so the corresponding constants grow with the Hilbert-space dimension even for a nonempty selected Pauli set. Thus the dimension-independent regime obtained for local product protocols relies on protocol-specific structure, not on the abstract finite-sample theorem alone.

The main contributions of the paper are as follows.

• Selected covariance matrices for shadow-output vectors. We formulate the covariance matrix \(\Sigma(\rho)=\operatorname{Cov}_\rho(\hat{x})\) for a general shadow-output vector and study its fixed selected compression \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\). This selected covariance matrix captures joint fluctuations of reconstructed output coordinates, not physical products of observables.

• Finite-sample selected covariance approximation. We prove a finite-sample theorem for selected sample-centered empirical covariances. If the selected dimension and a protocol-dependent one-shot boundedness parameter are independent of the ambient system size, then constant operator-norm accuracy is achieved with sample size independent of the ambient dimension. The proof uses matrix Bernstein concentration, the exact rank-one centering identity, and Weyl and Davis–Kahan perturbation bounds.

• General local product shadows. We show that dimension-independent selected covariance estimation is not specific to Pauli measurements. For local product shadow protocols on fixed-dimensional local systems, finite-weight product observables lead to bounds controlled by support sizes and local reconstruction coefficients, not by the total number of tensor factors.

• Explicit biased local Pauli formulas. For biased local Pauli shadows, we evaluate the quantities entering the general theorem in closed form. For a selected Pauli family \(J_l=\{P_1,\ldots,P_l\}\), the protocol-dependent constants in the finite-sample bound are determined explicitly by the supports of the \(P_a\)’s and the local basis-selection probabilities \(p_{j,r}\). In particular, a bounded number of selected Pauli strings, bounded Pauli weight, and a uniform lower bound on the probabilities \(p_{j,r}\) give a bound independent of the number of qubits. We also derive an exact covariance formula whose entries are governed by Pauli compatibility and inverse-probability overlap factors.

This work is complementary to the standard observable-wise classical-shadow theory, which is organized around shadow norms, marginal variances, and many-observable prediction guarantees [11], [12]. It is also related to Pauli-invariant shadow formalisms, which provide reconstruction maps and sample-complexity results for broad classes of Pauli-invariant ensembles [13]. Our focus is different: we study covariance matrices of reconstructed shadow outputs and the finite-sample recovery of selected covariance spectra.

Covariance quantities estimated from classical-shadow data also appear in the CoVaR approach to variational quantum algorithms [21]. In that setting, covariances between a Hamiltonian and an operator pool are used as nonlinear conditions for finding eigenstates. The present work is different in both object and guarantee: we study the covariance matrix of a shadow-output vector itself and prove finite-sample operator-norm and spectral-perturbation bounds for fixed selected covariance matrices. Table 1 summarizes this positioning.

Other recent directions study dynamics-informed or shallow-shadow regimes, robustness against corruption or heavy tails, and broader group-theoretic or continuous-variable shadow frameworks [15]–[20]. These works broaden the class of shadow protocols or modify the statistical model. The present paper has a different focus: it develops a covariance-matrix analysis for selected shadow-output coordinates, proves finite-sample operator-norm and spectral guarantees, extends the local bounded-output mechanism to general local product shadows, and treats biased local Pauli shadows as the main fully explicit example.

The paper is organized as follows. Section 2 introduces the general shadow-output notation, covariance matrices, empirical covariances, and selected covariance spectra. Section 3 proves the finite-sample selected covariance theorem. Section 4 specializes the framework to biased local Pauli shadows, derives the exact covariance formula, and compares the selected-output scaling with global Clifford shadows. Section 5 extends the local mechanism to general local product shadows and finite-weight product observables. Section 6 discusses the scope, interpretation, and future directions.

2 Preparatory covariance objects and notation↩︎

2.1 General finite-dimensional shadow channel↩︎

This section fixes the general shadow-output notation used throughout the paper. The construction is the standard classical-shadow measurement-channel formalism, written here for a finite randomized POVM. In the usual formulation of Huang–Kueng–Preskill, one samples a unitary, measures in a fixed basis, forms the raw snapshot, and applies the inverse of the associated shadow channel to obtain an unbiased single-shot reconstruction [11]. We use the same measurement-channel and inverse-channel viewpoint, but write it in a form that also covers finite randomized POVMs. This level of generality is consistent with the measurement-frame perspective on shadow tomography, where classical shadows are understood as unbiased estimators associated with a suitable dual frame [14]. For Pauli-invariant unitary ensembles, the corresponding Pauli-diagonal reconstruction-map formalism is developed in Bu–Koh–Garcia–Jaffe [13].

Although the original classical-shadow protocol is often introduced through random unitaries followed by computational-basis measurements [11], the same measurement-channel construction can be written for a finite randomized POVM. We use this slightly more general notation because it makes the later biased local Pauli specialization transparent. Although we use product notation below, the same definitions also apply to nonlocal measurement settings by regarding \(u\) as a global measurement setting and \(M_u\) as an arbitrary POVM on \(\mathcal{H}\).

Let \(\mathcal{H} = \bigotimes_{j=1}^k \mathcal{H}_j\) be a finite-dimensional multipartite Hilbert space. We describe a single round of a general randomized local measurement protocol. For each site \(j=1,\ldots,k\), let \(\mathcal{U}_j\) be a finite set of local measurement settings. For each \(u_j\in\mathcal{U}_j\), let \(M^{(j)}_{u_j} = \{M^{(j)}_{b_j|u_j}\}_{b_j\in\mathcal{B}_{j,u_j}}\) be a POVM on \(\mathcal{H}_j\), so that \[M^{(j)}_{b_j|u_j}\succeq 0, \qquad \sum_{b_j\in\mathcal{B}_{j,u_j}} M^{(j)}_{b_j|u_j} = I_j . \label{eq:general95local95povm95conditions}\tag{1}\]

The notation introduced so far is local to a single tensor factor \(\mathcal{H}_j\). For reference, Table 2 summarizes the local measurement notation before we pass to the global product measurement.

Set \[\mathcal{U} := \mathcal{U}_1\times\cdots\times\mathcal{U}_k . \label{eq:general95setting95space}\tag{2}\] In one round of the protocol, a random measurement setting \[\hat{u} = (\hat{u}_1,\ldots,\hat{u}_k) \in\mathcal{U} . \label{eq:general95random95measurement95setting}\tag{3}\] is drawn according to a prescribed probability distribution \(q\) on \(\mathcal{U}\). We do not assume here that \(q\) is a product distribution. Conditioned on \(\hat{u}=u=(u_1,\ldots,u_k)\), the product POVM \[M_u = \{M_{b|u}\}_{b\in\mathcal{B}_u}. \label{eq:general95product95povm95family}\tag{4}\] is measured, where \[\begin{align} \mathcal{B}_u &:= \mathcal{B}_{1,u_1}\times\cdots\times\mathcal{B}_{k,u_k}, \tag{5} \\ M_{b|u} &:= \bigotimes_{j=1}^k M^{(j)}_{b_j|u_j}, \qquad b=(b_1,\ldots,b_k)\in\mathcal{B}_u . \tag{6} \end{align}\] for \(b=(b_1,\ldots,b_k)\in\mathcal{B}_u\). Thus, conditioned on \(\hat{u}=u\), the measurement outcome \(\hat{b}=(\hat{b}_1,\ldots,\hat{b}_k)\in\mathcal{B}_u\) is obtained with probability \[\mathbb{P}_\rho(\hat{b}=b\mid \hat{u}=u) = \operatorname{Tr}(\rho M_{b|u}). \label{eq:general95conditional95outcome95probability}\tag{7}\]

For notational convenience, we regard the pair \((\hat{u},\hat{b})\) as the single-shot classical outcome of the randomized measurement protocol. The raw operator associated with the outcome \((u,b)\) is the corresponding POVM element \[\sigma_{u,b} := M_{b|u}. \label{eq:general95raw95operator95def}\tag{8}\] For the realized outcome, we write \[\hat{\sigma} := \sigma_{\hat{u},\hat{b}} = M_{\hat{b}|\hat{u}}. \label{eq:general95realized95raw95operator}\tag{9}\]

Let \(\mathsf L(\mathcal{H})\) denote the real vector space of Hermitian operators on \(\mathcal{H}\). Following the standard classical-shadow measurement-channel construction [11], define the associated shadow channel as the linear map \(\mathcal{M}:\mathsf L(\mathcal{H})\to\mathsf L(\mathcal{H})\) by \[\mathcal{M}[\rho] := \mathbb{E}_{\hat{u}} \sum_{b\in\mathcal{B}_{\hat{u}}} M_{b|\hat{u}}\, \operatorname{Tr}(\rho M_{b|\hat{u}}) . \label{eq:general95shadow95channel95expectation95form}\tag{10}\]

Equivalently, \[\mathcal{M}[\rho] = \sum_{u\in\mathcal{U}} q(u) \sum_{b\in\mathcal{B}_u} M_{b|u}\, \operatorname{Tr}(\rho M_{b|u}). \label{eq:general95shadow95channel95povm}\tag{11}\]

Throughout this general formulation, we assume that the randomized measurement protocol is tomographically complete, in the sense that the shadow channel \(\mathcal{M}:\mathsf L(\mathcal{H})\to\mathsf L(\mathcal{H})\) is invertible on the real vector space \(\mathsf L(\mathcal{H})\) of Hermitian operators. Thus \(\mathcal{M}^{-1}\) denotes the inverse linear map on \(\mathsf L(\mathcal{H})\).

The classical-shadow reconstruction applies the inverse channel to the raw snapshot [11]. Thus, for an observed outcome \((\hat{u},\hat{b})\), we set \[\hat{\rho}_{\hat{u},\hat{b}} := \mathcal{M}^{-1}(M_{\hat{b}|\hat{u}}) . \label{eq:general95reconstructed95snapshot}\tag{12}\] Since the mean raw snapshot is \(\mathcal{M}[\rho]\), the reconstructed snapshot is unbiased: \[\mathbb{E}_\rho \bigl[ \hat{\rho}_{\hat{u},\hat{b}} \bigr] = \mathcal{M}^{-1}\mathcal{M}[\rho] = \rho . \label{eq:general95snapshot95unbiasedness}\tag{13}\] Consequently, for every observable \(O_\alpha\) considered below, \[\mathbb{E}_\rho \operatorname{Tr}(\hat{\rho}_{\hat{u},\hat{b}}O_\alpha) = \operatorname{Tr}(\rho O_\alpha). \label{eq:general95snapshot95unbiasedness95observable95form}\tag{14}\] The notation up to this point describes one global measurement shot and its linear reconstruction on \(\mathcal{H}\). Table 3 summarizes the global measurement and reconstruction notation before we introduce the output-coordinate vector.

This formulation includes the usual randomized projective-measurement classical shadows of Huang–Kueng–Preskill as a special case [11]. This is the standard random-unitary shadow channel. The local and global Clifford protocols are obtained by choosing the corresponding unitary ensemble, while the biased local Pauli protocol can be written directly in the finite randomized POVM notation used above. For example, when the setting \(u\) is a unitary \(U\) and the POVM elements are \[M_{b|U} = U^\dagger |b\rangle\!\langle b|U . \label{eq:general95projective95shadow95povm95element}\tag{15}\] the above definition reduces to the standard shadow channel \[\mathcal{M}[\rho] = \mathbb{E}_{\hat{U}} \sum_b \hat{U}^\dagger |b\rangle\!\langle b|\hat{U}\, \operatorname{Tr}\!\left( \rho\,\hat{U}^\dagger |b\rangle\!\langle b|\hat{U} \right). \label{eq:general95projective95shadow95channel}\tag{16}\] The biased local Pauli protocol will be recovered later by taking \(\mathcal{U}_j=\{X,Y,Z\}\) and \(M^{(j)}_{s_j|r_j}=\frac{1}{2}(I+s_j r_j)\).

2.2 Shadow-output vectors and covariance matrices↩︎

We now pass from reconstructed shadow snapshots to the finite-dimensional real random vectors whose covariance matrices will be studied in the rest of the paper. Let \(\mathcal{O} = \{O_1,\ldots,O_p\}\) be a fixed finite collection of Hermitian observables on \(\mathcal{H}\). The identity observable may be included in this collection. These observables are not assumed to be measured directly. Rather, they specify the linear functionals of the state that are evaluated on the reconstructed shadow snapshot.

For a single randomized measurement outcome \((\hat{u},\hat{b})\), define \[\hat{x}_\alpha := \operatorname{Tr}\!\left( \hat{\rho}_{\hat{u},\hat{b}} O_\alpha \right) = \operatorname{Tr}\!\left( \mathcal{M}^{-1}(M_{\hat{b}|\hat{u}}) O_\alpha \right), \label{eq:general95shadow95output95coordinate}\tag{17}\] for \(\alpha=1,\ldots,p\). We write \[\hat{x} := (\hat{x}_1,\ldots,\hat{x}_p)^\top \in\mathbb{R}^p . \label{eq:general95shadow95output95vector}\tag{18}\] for the resulting single-shot shadow-output vector. Its mean is denoted by \[m_\alpha := \operatorname{Tr}(\rho O_\alpha), \qquad m:=(m_1,\ldots,m_p)^\top\in\mathbb{R}^p . \label{eq:general95shadow95output95mean95vector}\tag{19}\] By the unbiasedness of the reconstructed snapshot, \(\mathbb{E}_\rho[\hat{\rho}_{\hat{u},\hat{b}}]=\rho,\) we have \[\mathbb{E}_\rho[\hat{x}_\alpha] = \operatorname{Tr}(\rho O_\alpha) = m_\alpha, \qquad \alpha=1,\ldots,p. \label{eq:general95shadow95output95coordinate95unbiasedness}\tag{20}\] or equivalently \[\mathbb{E}_\rho[\hat{x}]=m. \label{eq:general95shadow95output95unbiasedness}\tag{21}\]

The raw second-moment matrix of the shadow-output vector is \[M^{(2)}(\rho) := \mathbb{E}_\rho[\hat{x}\hat{x}^\top] \in\mathbb{R}^{p\times p}. \label{eq:general95raw95second95moment}\tag{22}\] The associated covariance matrix is \[\Sigma(\rho) := \operatorname{Cov}_\rho(\hat{x}) = M^{(2)}(\rho)-mm^\top \in\mathbb{R}^{p\times p}. \label{eq:general95shadow95output95covariance}\tag{23}\] Equivalently, its entries are \[\Sigma_{\alpha\beta}(\rho) = \operatorname{Cov}_\rho(\hat{x}_\alpha,\hat{x}_\beta) = \mathbb{E}_\rho[\hat{x}_\alpha\hat{x}_\beta] - m_\alpha m_\beta . \label{eq:general95covariance95entries}\tag{24}\]

The off-diagonal entries of \(\Sigma(\rho)\) encode statistical couplings between different shadow-output coordinates. They are properties of the joint distribution of the reconstructed snapshot \(\hat{\rho}_{\hat{u},\hat{b}}\), and should not be confused with directly measuring the products of the observables \(O_\alpha O_\beta\). In particular, the measurement actually performed in one shot is the POVM \(M_{\hat{u}}\), while the quantities \(\hat{x}_\alpha\) are obtained by post-processing the reconstructed shadow snapshot.

We have also introduced the output-coordinate vector obtained by evaluating a fixed observable family on the reconstructed snapshot, together with its moments. Table 4 summarizes this single-shot output and covariance notation.

2.3 Independent samples and empirical covariance matrices↩︎

We now pass from the single-shot shadow-output distribution to a finite sample. Fix a sample size \(N\). For \(i=1,\ldots,N\), let \((\hat{u}_i,\hat{b}_i)\) be independent outcomes generated by applying the same randomized shadow measurement protocol to independent copies of the state \(\rho\). Thus \(\hat{u}_i\) is the randomized measurement setting in the \(i\)-th shot, and \(\hat{b}_i\) is the corresponding measurement outcome.

The reconstructed shadow snapshot in the \(i\)-th shot is \[\hat{\rho}_i := \hat{\rho}_{\hat{u}_i,\hat{b}_i} = \mathcal{M}^{-1}(M_{\hat{b}_i|\hat{u}_i}). \label{eq:general95sample95reconstructed95snapshot}\tag{27}\] For the fixed observable family \(\mathcal{O}=\{O_1,\ldots,O_p\}\), we define the associated shadow-output vector by \[\hat{x}_{i,\alpha} := \operatorname{Tr}(\hat{\rho}_i O_\alpha), \qquad \alpha=1,\ldots,p, \label{eq:general95sample95shadow95output95coordinate}\tag{28}\] and write \[\hat{x}_i := (\hat{x}_{i,1},\ldots,\hat{x}_{i,p})^\top \in\mathbb{R}^p . \label{eq:general95sample95shadow95output95vector}\tag{29}\] Then \(\hat{x}_1,\ldots,\hat{x}_N\) are independent copies of the single-shot shadow-output vector \(\hat{x}\) defined in the previous subsection. In particular, \[\mathbb{E}_\rho[\hat{x}_i]=m, \qquad \operatorname{Cov}_\rho(\hat{x}_i)=\Sigma(\rho), \qquad i=1,\ldots,N. \label{eq:general95sample95mean95covariance}\tag{30}\]

The empirical shadow-output mean is \[\bar{\hat{x}}_N := \frac{1}{N}\sum_{i=1}^N \hat{x}_i . \label{eq:general95empirical95shadow95output95mean}\tag{31}\] We distinguish the true-centered empirical covariance matrix \[\widehat\Sigma_N^{\mathrm{tc}} := \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-m)(\hat{x}_i-m)^\top \in\mathbb{R}^{p\times p} \label{eq:true95centered95empirical95covariance95general}\tag{32}\] from the sample-centered empirical covariance matrix \[\widehat\Sigma_N^{\mathrm{sc}} := \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-\bar{\hat{x}}_N) (\hat{x}_i-\bar{\hat{x}}_N)^\top \in\mathbb{R}^{p\times p}. \label{eq:sample95centered95empirical95covariance95general}\tag{33}\] We use the \(1/N\)-normalized empirical covariance throughout; no unbiased \(1/(N-1)\) correction is used. The true-centered covariance is useful for concentration estimates because its summands are centered with respect to the mean \(m\). The sample-centered covariance is the empirical covariance matrix computed directly from data, since \(m\) is unknown in a tomography problem.

The two empirical covariance matrices are related by the exact identity \[\begin{align} \widehat\Sigma_N^{\mathrm{sc}} - \widehat\Sigma_N^{\mathrm{tc}} = - (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top , \label{eq:exact95rank95one95relation95true95sample95centered95general} \end{align}\tag{34}\] which is shown in Appendix 7.

Thus the sample-centered covariance is a rank-one negative perturbation of the true-centered covariance. This rank-one relation will be used later to transfer the finite-sample operator-norm approximation first proved for the true-centered selected covariance to the sample-centered selected covariance computed from the measurement data.

2.4 Selected covariance matrices and their spectral meaning↩︎

The covariance matrix \(\Sigma(\rho)\) defined above describes the joint fluctuation structure of the full shadow-output vector \(\hat{x}=(\hat{x}_1,\ldots,\hat{x}_p)^\top .\) In the finite-sample analysis below, however, we do not attempt to estimate the full covariance matrix in operator norm. Instead, we restrict attention to a fixed finite set of selected output coordinates.

Let \[J_l = \{\alpha_1,\ldots,\alpha_l\} \subset\{1,\ldots,p\}. \label{eq:general95selected95coordinate95set}\tag{35}\] be a deterministic selected coordinate set, fixed independently of the measurement data. We assume that the indices in \(J_l\) are distinct. Let \[R_l\in\{0,1\}^{l\times p}. \label{eq:general95selection95matrix95def}\tag{36}\] be the corresponding coordinate-selection matrix. Thus the \(a\)-th row of \(R_l\) has a single nonzero entry, equal to \(1\), in the column \(\alpha_a\).

Equivalently, the rows of \(R_l\) are orthonormal coordinate vectors, and hence \[R_lR_l^\top=I_l. \label{eq:general95selection95matrix95property}\tag{37}\] Thus \(R_l\) acts by retaining precisely the coordinates indexed by \(J_l=\{\alpha_1,\ldots,\alpha_l\}\) and discarding all remaining coordinates. More explicitly, for any vector \(x\in\mathbb{R}^p\), \[R_lx = \begin{pmatrix} x_{\alpha_1}\\ \vdots\\ x_{\alpha_l} \end{pmatrix} \in\mathbb{R}^l . \label{eq:general95selection95matrix95action}\tag{38}\] Applying this coordinate projection to the shadow-output vector gives the selected shadow-output vector \[R_l\hat{x} = \begin{pmatrix} \hat{x}_{\alpha_1}\\ \vdots\\ \hat{x}_{\alpha_l} \end{pmatrix}. \label{eq:general95selected95shadow95output95vector}\tag{39}\] Its mean is obtained by applying the same selection matrix to the full mean vector \(m\): \[R_lm = \begin{pmatrix} m_{\alpha_1}\\ \vdots\\ m_{\alpha_l} \end{pmatrix}. \label{eq:general95selected95mean95vector}\tag{40}\]

We define the selected covariance matrix by compressing the full shadow-output covariance matrix with the fixed selection matrix \(R_l\): \[\Sigma_l(\rho) := R_l\Sigma(\rho)R_l^\top \in\mathbb{R}^{l\times l}. \label{eq:selected95covariance95general}\tag{41}\] Since \(R_l\) is deterministic, this compressed matrix is exactly the covariance matrix of the selected shadow-output vector \(R_l\hat{x}\). Using \(\mathbb{E}_\rho[\hat{x}]=m\), we have \[\Sigma_l(\rho) = \operatorname{Cov}_\rho(R_l\hat{x}) = \mathbb{E}_\rho \left[ (R_l\hat{x}-R_lm)(R_l\hat{x}-R_lm)^\top \right]. \label{eq:selected95covariance95cov95form95general}\tag{42}\] Equivalently, if \(J_l=\{\alpha_1,\ldots,\alpha_l\}\), then the \((a,b)\)-entry of \(\Sigma_l(\rho)\) is \[\bigl(\Sigma_l(\rho)\bigr)_{ab} = \operatorname{Cov}_\rho \bigl( \hat{x}_{\alpha_a}, \hat{x}_{\alpha_b} \bigr), \qquad 1\le a,b\le l . \label{eq:selected95covariance95entry95general}\tag{43}\] Thus \(\Sigma_l(\rho)\) records the joint fluctuation structure of the fixed selected coordinates \(\hat{x}_{\alpha_1},\ldots,\hat{x}_{\alpha_l}\), rather than the covariance structure of the full ambient vector \(\hat{x}\in\mathbb{R}^p\).

The spectral meaning of \(\Sigma_l(\rho)\) is immediate from its variational characterization. For any vector \(v\in\mathbb{R}^l\), the scalar random variable \(\langle v,R_l\hat{x}\rangle\) is the shadow-output estimate of the linear combination of selected coordinates specified by \(v\). Its variance is \[\operatorname{Var}_\rho(\langle v,R_l\hat{x}\rangle) = v^\top \Sigma_l(\rho)v . \label{eq:selected95covariance95variance95form95general}\tag{44}\] Consequently, for unit vectors \(v\in\mathbb{R}^l\), we have \[\max_{\|v\|_2=1} \operatorname{Var}_\rho(\langle v,R_l\hat{x}\rangle) = \max_{\|v\|_2=1} v^\top\Sigma_l(\rho)v = \lambda_{\max}(\Sigma_l(\rho)). \label{eq:selected95covariance95max95variance95general}\tag{45}\] Similarly, \[\min_{\|v\|_2=1} \operatorname{Var}_\rho(\langle v,R_l\hat{x}\rangle) = \lambda_{\min}(\Sigma_l(\rho)). \label{eq:selected95covariance95min95variance95general}\tag{46}\] Thus the eigenvalues of \(\Sigma_l(\rho)\) describe the extremal fluctuation scales among all normalized linear combinations of the selected shadow-output coordinates. The corresponding eigenvectors identify the linear combinations that realize these extremal variances.

This spectral interpretation is the reason for estimating selected covariance matrices in operator norm. If an empirical selected covariance matrix \(\widehat\Sigma_{l,N}\) satisfies \[\|\widehat\Sigma_{l,N}-\Sigma_l(\rho)\|_{\mathrm{op}} \le \varepsilon , \label{eq:selected95operator95norm95approx95assumption95general}\tag{47}\] then all selected variance scales are uniformly approximated: \[\left| v^\top\widehat\Sigma_{l,N}v - v^\top\Sigma_l(\rho)v \right| \le \varepsilon \qquad \text{for all } \|v\|_2=1. \label{eq:selected95operator95norm95variance95approx95general}\tag{48}\] In particular, the selected eigenvalues are stable under such an operator-norm perturbation, and isolated selected spectral subspaces can be controlled by standard spectral perturbation theory. These consequences will be made precise in the finite-sample results below.

Finally, we emphasize that the selection matrix \(R_l\) is fixed before the data are observed. The results below are therefore fixed-selection statements. They do not address data-dependent or adaptive choices of selected coordinates, which would require additional post-selection or uniform concentration arguments.

3 Finite-Sample Spectral Approximation of Selected Shadow Covariances↩︎

We now prove the finite-sample selected covariance-estimation theorem in a form that does not rely on the local Pauli structure. The goal is to control the selected covariance matrix \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\) from independent shadow samples, and then to convert this operator-norm control into selected eigenvalue and spectral-projector guarantees.

3.1 Overview of the selected finite-sample argument↩︎

The input is the general shadow-output vector \(\hat{x}=(\hat{x}_1,\ldots,\hat{x}_p)^\top\in\mathbb{R}^p\) introduced in Section 2, together with a fixed coordinate-selection matrix \(R_l\). The covariance matrix is \(\Sigma(\rho)=\operatorname{Cov}_\rho(\hat{x})\), and the selected covariance matrix is \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\).

The point of this section is deliberately selected rather than full-dimensional. We do not attempt to estimate the full \(p\times p\) covariance matrix in operator norm. Instead, all concentration takes place after the fixed compression \(R_l\). Thus the matrix dimension entering the probability bounds is the selected dimension \(l\), not the ambient dimension \(p\).

The probabilistic input is standard. Once the selected centered vectors \(R_l(\hat{x}_i-m)\) are uniformly bounded in Euclidean norm, the true-centered empirical covariance error becomes an average of independent centered self-adjoint matrices. We control this average by the self-adjoint matrix Bernstein inequality in the form of Tropp [4]. This gives an operator-norm perturbation bound for the selected covariance. Standard spectral perturbation results, namely Weyl’s inequality and the Davis–Kahan theorem, then convert this operator-norm control into selected eigenvalue and spectral-projector bounds [5]–[7], [9].

There is one statistical complication: the mean \(m\) is unknown in data analysis. The covariance matrix computed from data is therefore sample-centered, not true-centered. We first prove concentration for the true-centered selected covariance, and then use the exact rank-one centering identity to pass to the sample-centered covariance. The additional empirical-mean term is controlled by a coordinatewise Hoeffding inequality and a union bound over the selected coordinates [9], [22].

The section is organized as follows. We first introduce the selected one-shot radius and state the sample-centered main theorem. We then prove the true-centered matrix concentration bound, derive its spectral consequences, and finally transfer the result to the sample-centered covariance using the rank-one centering identity.

3.2 Selected one-shot radius and main theorem↩︎

This subsection states the main finite-sample result of the paper in the form used in the concrete protocol sections. The statement is intentionally a constant-error theorem under a bounded selected radius. It says that, if the selected dimension and the selected centered one-shot radius remain bounded independently of the ambient system size, then the selected sample-centered covariance matrix can be estimated in operator norm with a sample size independent of the ambient system size.

All protocol dependence enters through the selected one-shot radius. The later local Pauli and general local product sections are devoted to computing or bounding this radius in concrete measurement models. The theorem is not tied to Pauli compatibility, locality, or an explicit covariance formula; those structures are used only later to verify the bounded-radius hypothesis.

The theorem is stated for the sample-centered covariance because this is the matrix computed directly from observed shadow data. Its proof is given in the rest of this section. We first establish a concentration bound for the true-centered selected empirical covariance, and then transfer it to the sample-centered covariance by using the exact rank-one centering identity. Along the way, we also obtain more detailed eigenvalue and spectral-projector perturbation bounds.

Let \(R_l\in\{0,1\}^{l\times p}\) be the fixed coordinate-selection matrix introduced in Section 2.4, so that \[R_lR_l^\top=I_l. \label{eq:fs95general95selection95matrix95orthonormal95rows}\tag{49}\] We define the selected sample-centered and true-centered empirical covariance matrices by \[\begin{align} \widehat\Sigma_{l,N}^{\mathrm{sc}} &:= R_l\widehat\Sigma_N^{\mathrm{sc}}R_l^\top, \tag{50} \\ \widehat\Sigma_{l,N}^{\mathrm{tc}} &:= R_l\widehat\Sigma_N^{\mathrm{tc}}R_l^\top . \tag{51} \end{align}\] Equivalently, we have \[\begin{align} \widehat\Sigma_{l,N}^{\mathrm{sc}} =& \frac{1}{N}\sum_{i=1}^N (R_l\hat{x}_i-R_l\bar{\hat{x}}_N) (R_l\hat{x}_i-R_l\bar{\hat{x}}_N)^\top , \tag{52} \\ \widehat\Sigma_{l,N}^{\mathrm{tc}} =& \frac{1}{N}\sum_{i=1}^N (R_l\hat{x}_i-R_lm) (R_l\hat{x}_i-R_lm)^\top \tag{53}\\ =& \frac{1}{N}\sum_{i=1}^N \hat{z}_i \hat{z}_i^\top, \tag{54} \end{align}\] where the selected centered vector \(\hat{z}_i\) is defined as \[\hat{z}_i:=R_l(\hat{x}_i-m)\in\mathbb{R}^l, \qquad i=1,\ldots,N. \label{eq:selected95centered95vector95general}\tag{55}\]

The next definition isolates the deterministic boundedness assumption needed for the matrix Bernstein argument. It should be read as a selected analogue of the usual bounded-sample assumption in nonasymptotic covariance estimation. The full vector \(\hat{x}\) may live in a very large ambient space, but the theorem uses only the Euclidean radius of the selected vector \(R_l\hat{x}\).

The finite-sample estimates below are built from the one-shot selected output, its covariance, and the corresponding selected-radius parameters. Table 5 collects this notation before we introduce the \(N\)-sample quantities and error levels.

The quantity \(B_l\) bounds the raw selected one-shot output, while \(A_l(\rho)\) bounds the centered selected output. For \(\delta\in(0,1)\), define the true-centered error level by \[\varepsilon_{l,N}^{\mathrm{tc}}(\delta) := A_l(\rho)^2 \left( \sqrt{\frac{8\log(2l/\delta)}{N}} + \frac{4\log(2l/\delta)}{3N} \right). \label{eq:epsilon95tc95general}\tag{60}\] Define also the selected empirical-mean correction by \[\mu_{l,N}(\delta) := \frac{2lA_l(\rho)^2\log(2l/\delta)}{N}. \label{eq:mu95general}\tag{61}\] The sample-centered error level is \[\varepsilon_{l,N}^{\mathrm{sc}}(\delta) := \varepsilon_{l,N}^{\mathrm{tc}}(\delta/2) + \mu_{l,N}(\delta/2). \label{eq:epsilon95sc95def95general}\tag{62}\] Equivalently, substituting 60 and 61 into 62 gives \[\varepsilon_{l,N}^{\mathrm{sc}}(\delta) = A_l(\rho)^2 \left( \sqrt{\frac{8\log(4l/\delta)}{N}} + \left(\frac{4}{3}+2l\right) \frac{\log(4l/\delta)}{N} \right). \label{eq:epsilon95sc95general}\tag{63}\]

We next collect the notation that depends on the \(N\) independent samples and the finite-sample error parameters used in the main theorem. Table 6 separates these empirical objects from the one-shot quantities listed in Table 5.

We now state the main finite-sample theorem for the sample-centered selected covariance. This is the covariance matrix computed directly from data, because the mean \(m\) is unknown. The theorem gives a high-probability operator-norm bound for its distance from the selected covariance \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top .\) The more detailed selected eigenvalue and spectral-projector consequences are stated below as refinements of the same operator-norm approximation.

The error bound depends on four quantities: the selected centered radius \(A_l(\rho)\), the selected dimension \(l\), the sample size \(N\), and the failure probability \(\delta\). The ambient dimension \(p\) does not appear in the logarithmic factor, because the covariance matrix is compressed by the fixed selection matrix \(R_l\) before concentration is applied.

The assumptions \(B_l<\infty\) and \(A_l(\rho)\le A_0\) are stated abstractly in this general section. They should be understood as protocol-dependent boundedness conditions on the selected one-shot output and on its centered version. In applications, one may verify the second condition either by computing the centered radius \(A_l(\rho)\) directly or by using the deterministic upper bound \(A_l(\rho) \le \widetilde{A}_l(\rho) = B_l+\|R_lm\|_2 .\) For the biased local Pauli protocol, the raw selected radius \(B_l\) can be computed explicitly from the local basis-selection probabilities and the supports of the selected Pauli coordinates. Together with the elementary bound \(\|R_lm\|_2\le \sqrt l\) for Pauli expectation values, this gives a state-uniform bound on \(A_l(\rho)\) whenever the selected set size, the Pauli weights, and the inverse local basis probabilities are uniformly bounded.

More generally, for local product shadow protocols with finite-weight selected product observables, the raw selected radius \(B_l\), and hence the centered radius \(A_l(\rho)\), can be bounded in terms of local reconstruction coefficients and support sizes rather than the total number of tensor factors. These protocol-specific estimates are carried out in the later sections and verify the bounded-radius hypothesis used in Theorem 1.

3.3 True-centered selected covariance concentration↩︎

We now begin the proof of Theorem 1. The first step is to analyze the true-centered selected covariance. This is the random-matrix core of the section. After true centering, the covariance error can be written as \[\widehat\Sigma_{l,N}^{\mathrm{tc}}-\Sigma_l(\rho) = \frac{1}{N}\sum_{i=1}^N \left(\hat{z}_i\hat{z}_i^\top-\Sigma_l(\rho)\right),\] where the summands are independent, centered, self-adjoint \(l\times l\) matrices. This is exactly the setting of the self-adjoint matrix Bernstein inequality. We use Tropp’s noncommutative Bernstein bound [4]; the proof below verifies its two inputs: an almost-sure operator-norm bound on each summand and a variance proxy bound for the sum of squared summands.

The role of the selected radius is transparent here. The bound \(\|\hat{z}_i\|_2\le A_l(\rho)\) implies \(\|\hat{z}_i\hat{z}_i^\top\|_{\mathrm{op}}\le A_l(\rho)^2,\) which in turn gives both the summand norm bound and the variance proxy. Thus the concentration rate is governed by \(A_l(\rho)^2\), and the matrix dimension factor is \(l\).

Since 55 implies \[\mathbb{E}_\rho[\hat{z}_i]=0, \label{eq:selected95centered95vector95mean95zero95general}\tag{68}\] the definition of the selected covariance implies \[\Sigma_l(\rho) = \mathbb{E}_\rho[\hat{z}_i\hat{z}_i^\top]. \label{eq:selected95covariance95z95second95moment95general}\tag{69}\] The following elementary bound is the only place where the centered selected radius \(A_l(\rho)\) enters the true-centered concentration argument.

The next proposition is the random-matrix concentration input for the true-centered covariance. We use the self-adjoint matrix Bernstein inequality of Tropp [4], which controls the operator norm of a sum of independent centered self-adjoint random matrices through a uniform summand bound and a matrix variance proxy. In the present setting, the summands are \(\hat{X}_i=\hat{z}_i\hat{z}_i^\top-\Sigma_l(\rho)\), where \(\hat{z}_i=R_l(\hat{x}_i-m)\). The selected centered-radius bound \(\|\hat{z}_i\|_2\le A_l(\rho)\) supplies both Bernstein inputs: it gives the almost-sure summand bound and the variance-proxy bound used in the proof below. No coordinatewise independence of \(R_l\hat{x}_i\) is assumed; only independence across shadow shots is used.

3.4 True-centered selected spectral approximation↩︎

The concentration proposition gives an operator-norm perturbation bound for the selected covariance matrix. The next theorem packages the deterministic spectral consequences of this perturbation. The two-sided matrix inequality and the operator-norm bound are immediate once the deviation matrix is controlled in operator norm. Eigenvalue stability follows from Weyl’s inequality for Hermitian matrices [5]. Stability of isolated spectral subspaces is obtained from the Davis–Kahan theorem, in a standard gap-dependent form [6]–[9].

This theorem is stated separately because it clarifies which part of the argument is probabilistic and which part is deterministic. Probability enters only through the event on which \(\|\widehat\Sigma_{l,N}^{\mathrm{tc}}-\Sigma_l(\rho)\|_{\mathrm{op}}\) is small. Once that event holds, the eigenvalue and projector estimates are ordinary perturbation theory.

The true-centered covariance is an intermediate object. The next subsection transfers the preceding bounds to the sample-centered selected covariance by using the exact rank-one centering relation and a concentration bound for the selected empirical mean.

3.5 Centering correction for selected empirical covariance↩︎

We first record the exact algebraic relation between the true-centered and sample-centered selected empirical covariance matrices. This step is independent of the shadow protocol. It only uses the definitions of empirical centering and the fixed selection matrix \(R_l\).

Recall that \(\hat{z}_i = R_l(\hat{x}_i-m)\) for \(i=1,\ldots,N\), as defined in 55 . The selected empirical mean of the centered outputs is \[\bar{\hat{z}}_N := R_l(\bar{\hat{x}}_N-m) = \frac{1}{N}\sum_{i=1}^N \hat{z}_i . \label{eq:selected95empirical95mean95z95general}\tag{111}\] The point of this subsection is that replacing the unknown mean \(m\) by the empirical mean \(\bar{\hat{x}}_N\) introduces a correction that is not arbitrary: after selection, it is exactly the negative rank-one matrix \(-\bar{\hat{z}}_N\bar{\hat{z}}_N^\top\).

It remains to control the rank-one correction in 113 . This evaluation reduces to bounding the selected empirical mean because \(\left\| \bar{\hat{z}}_N\bar{\hat{z}}_N^\top \right\|_{\mathrm{op}} = \|\bar{\hat{z}}_N\|_2^2\). This remaining issue is done in the next subsection.

3.6 Empirical mean control and sample-centering transfer↩︎

The centering identity in Lemma 2 reduces the passage from true-centered to sample-centered covariance to a bound on \(\|\bar{\hat{z}}_N\|_2^2\). We control this quantity using the selected centered-radius bound 70 , together with coordinatewise Hoeffding inequalities and a union bound over the \(l\) selected coordinates.

We now bound the selected empirical mean. The rank-one correction is controlled by its only nonzero eigenvalue, \(\|\bar{\hat{z}}_N\bar{\hat{z}}_N^\top\|_{\mathrm{op}} = \|\bar{\hat{z}}_N\|_2^2 .\) The coordinates of \(\hat{z}_i\) need not be independent within a single shot. We therefore use only independence across shots: each selected coordinate average is a bounded scalar average. Hoeffding’s inequality is applied coordinatewise, and a union bound over the \(l\) selected coordinates gives the Euclidean-norm bound.

The next lemma records how the true-centered covariance bound is converted into a sample-centered covariance bound. The conversion uses only the rank-one centering identity from Lemma 2. Since sample centering subtracts the positive semidefinite matrix \(\bar{\hat{z}}_N\bar{\hat{z}}_N^\top\), the upper matrix inequality is unchanged, while the lower matrix inequality loses the additional empirical-mean term \(\mu\). This asymmetry is later absorbed into the symmetric operator-norm bound.

3.7 Sample-centered selected spectral approximation↩︎

We now combine the true-centered spectral approximation with the centering correction estimates. The resulting theorem gives the detailed sample-centered operator-norm, eigenvalue, and spectral-projector bounds. This result is stronger than what is needed for the constant-radius main theorem, but it is useful because it makes explicit the spectral consequences of the selected covariance approximation.

The proof combines the three estimates established above: the true-centered spectral approximation, the Hoeffding bound for the selected empirical mean, and the exact rank-one identity relating the true-centered and sample-centered covariances. A union bound allocates half of the failure probability to the true-centered event and half to the empirical-mean event. Once the operator-norm bound for the sample-centered covariance is obtained, the eigenvalue and spectral-projector statements again follow from Weyl’s inequality [5] and the Davis–Kahan perturbation theorem [6]–[9].

4 Local Pauli shadows↩︎

In this section we specialize the general selected-covariance framework to the \(k\)-qubit biased local Pauli shadow protocol. The main purpose of this section is first to evaluate the selected one-shot radius \(B_l\) for a fixed selected Pauli coordinate set and thereby to instantiate the general finite-sample selected covariance theorem in the local Pauli setting.

After this finite-sample specialization, we turn to an analytic calculation of the covariance entries for the same protocol. This latter calculation uses structural features that are specific to qubit Pauli measurements: each non-identity local Pauli observable is one of \(X,Y,Z\), local measurement outcomes take values in \(\{\pm1\}\), and Pauli strings compatible with the same local basis pattern have a simple cancellation structure. Thus the closed covariance formulas below should be understood as additional local-Pauli structure, not as assumptions needed for the abstract finite-sample theorem.

Biased local Pauli measurements have been studied previously as a variance-reduction tool for Hamiltonian expectation estimation [10]. Here we use the same biased basis-selection mechanism to analyze the full covariance structure of selected reconstructed Pauli coordinates.

4.1 Formulation of biased qubit local Pauli shadows↩︎

We now specialize the abstract shadow-output framework of Section 2 to the \(k\)-qubit local Pauli setting. There are two distinct choices to be fixed.

First, we fix the physical quantities whose expectation values are represented by the coordinates of the shadow-output vector. In the general framework, these quantities were specified by an observable family \(\mathcal{O}=\{O_1,\ldots,O_p\}.\) In the present section, we take this observable family to be the set of all non-identity \(k\)-qubit Pauli strings. Let \(\mathcal{P}_k:= \{I,X,Y,Z\}^{\otimes k}\) be the \(k\)-qubit Pauli string set, and let \(\mathcal{P}_k^\times := \mathcal{P}_k\setminus\{I^{\otimes k}\}\) be the set of non-identity Pauli strings. Thus, in the notation of Section 2.2, we specialize \[\mathcal{O} = \{O_1,\ldots,O_p\} = \mathcal{P}_k^\times, \qquad p=4^k-1. \label{eq:local95pauli95observable95family95specialization}\tag{157}\] Each coordinate of the general shadow-output vector \(\hat{x}\in\mathbb{R}^p\) in 18 is therefore indexed by a non-identity Pauli string \(P\in\mathcal{P}_k^\times\). We write this coordinate as \[\hat{x}_P := \operatorname{Tr}(\hat{\rho}\,P), \qquad P\in\mathcal{P}_k^\times, \label{eq:local95pauli95shadow95output95coordinate95specialization}\tag{158}\] where \(\hat{\rho}\) is the single-shot reconstructed shadow snapshot. The corresponding coordinate is the Pauli expectation value \[m_P := \operatorname{Tr}(\rho P), \qquad P\in\mathcal{P}_k^\times. \label{eq:local95pauli95mean95coordinate95specialization}\tag{159}\] At this stage, we have only fixed the target observables whose reconstructed coefficients form the vector \(\hat{x}\); no physical covariance between observables is being assumed or directly measured.

Second, we fix the measurement mechanism that generates these reconstructed coefficients. In the general framework, this was an arbitrary randomized measurement protocol with a corresponding reconstruction map. Here, we use the biased local Pauli shadow protocol. At each site \(j=1,\ldots,k\), a local Pauli basis \(\hat{r}_j\in\{X,Y,Z\}\) is chosen with probabilities \[\begin{align} & \mathbb{P}(\hat{r}_j=r) = p_{j,r}, \qquad r\in\{X,Y,Z\},\notag\\ & p_{j,X},p_{j,Y},p_{j,Z}>0, \quad p_{j,X}+p_{j,Y}+p_{j,Z}=1. \label{eq:local95pauli95basis95probabilities} \end{align}\tag{160}\] We write \[\hat{r} = (\hat{r}_1,\ldots,\hat{r}_k) \in \{X,Y,Z\}^k \label{eq:local95pauli95basis95pattern}\tag{161}\] for the resulting local basis pattern. Conditional on \(\hat{r}\), each site is measured in the corresponding Pauli basis, and the outcome at site \(j\) is \(\hat{s}_j\in\{\pm1\}\). Thus one shot of the measurement protocol produces \[(\hat{r},\hat{s}) = \bigl((\hat{r}_1,\ldots,\hat{r}_k),(\hat{s}_1,\ldots,\hat{s}_k)\bigr). \label{eq:local95pauli95single95shot95data}\tag{162}\]

The full \(k\)-qubit snapshot \(\hat{\rho}\) is given as \[\begin{align} \hat{\rho}=\bigotimes_{j=1}^k \frac{1}{2}\Bigl(I+\frac{\hat{s}_j}{p_{j,\hat{r}_j}}\,\hat{r}_j\Bigr), \label{KJ13} \end{align}\tag{163}\] which is shown in Appendix 8. For \(P\in\mathcal{P}_k^\times\), we write \(\require{physics} \hat{x}_P=\Tr(\hat{\rho} P)\) and \(\require{physics} m_P=\Tr(\rho P)\), identifying each element of \(\mathcal{P}_k^\times\) with its coordinate index. Once these two components have been fixed, the covariance matrix considered in this section is the covariance of the reconstructed Pauli coefficients \((\hat{x}_P)_{P\in\mathcal{P}_k^\times}\). Thus the entries of \(\Sigma(\rho)\) are \[\Sigma_{PQ}(\rho) = \operatorname{Cov}_\rho(\hat{x}_P,\hat{x}_Q) = \mathbb{E}_\rho[\hat{x}_P\hat{x}_Q]-m_Pm_Q \label{eq:local95pauli95covariance95entry95specialization}\tag{164}\] for \(P,Q\in\mathcal{P}_k^\times\). This is a covariance of the random reconstructed coefficients produced by the shadow protocol. It should not be confused with directly measuring a physical product observable such as \(PQ\).

For a Pauli string \(P\in\mathcal{P}_k\), we define its support and its weight as \[\mathop{\mathrm{supp}}(P) := \{j:\,P_j\neq I\}, \quad \mathop{\mathrm{wt}}(P):=| \mathop{\mathrm{supp}}(P)|. \label{eq:local95pauli95support95def}\tag{165}\] We say that a local basis pattern \(r=(r_1,\ldots,r_k)\in\{X,Y,Z\}^k\) is compatible with \(P\), and write \(r\succeq P\), if \(r_j=P_j\) for every \(j\in\mathop{\mathrm{supp}}(P)\). This compatibility condition connects the target observable \(P\) with the measurement mechanism: the reconstructed coefficient \(\hat{x}_P\) can be nonzero only when the chosen basis pattern is compatible with \(P\).

For \(P,Q\in\mathcal{P}_k\), we say that \(P\) and \(Q\) are compatible, and write \(P\sim Q\), if on every site where both are non-identity, they carry the same Pauli. If \(P\sim Q\), define \(P\ominus Q\in\mathcal{P}_k\) by canceling the common non-identity factors and retaining the remaining ones. For later use, define the overlap bias factor \[\begin{align} \beta(P,Q):=\prod_{j\in\mathop{\mathrm{supp}}(P)\cap\mathop{\mathrm{supp}}(Q)} p_{j,P_j}^{-1}, \qquad (P\sim Q), \label{BSHJ3} \end{align}\tag{166}\] which is well defined because compatibility implies \(P_j=Q_j\) on the overlap. The coefficient rule below is the basic local Pauli input for both parts of this section. It is used first to compute the selected one-shot radius, and later to derive the exact covariance formula.

This explicit coefficient rule is the point at which the present subsection becomes specific to qubit local Pauli shadows. The exact covariance formula below is obtained by analyzing products \(\hat{x}_P\hat{x}_Q\) using Lemma 5.

4.2 Selected radius and finite-sample consequences↩︎

We now connect the preceding local Pauli formulation with the general finite-sample selected covariance theorem of Section 3. That theorem applies to arbitrary shadow protocols once the protocol-dependent selected one-shot radius is identified. Thus, in the biased local Pauli setting, the remaining task is to compute the selected radius \(B_l\) appearing in Definition 1 and then substitute it into the sample-centered selected covariance theorem.

Let \(J_l=\{P_1,\ldots,P_l\} \subset \mathcal{P}_k^\times\) be a fixed selected set of distinct non-identity Pauli strings, and let \(R_l\) be the corresponding coordinate-selection matrix. Define \[B_{l,\mathrm{LP}}^2 := \max_{r\in\{X,Y,Z\}^k} \sum_{a=1}^l \left( \prod_{j\in\mathop{\mathrm{supp}}(P_a)} p_{j,P_{a,j}}^{-2} \right) \mathbf{1}_{\{r\succeq P_a\}}. \label{eq:local95pauli95Bl95LP95def}\tag{170}\] This is the deterministic selected one-shot radius predicted by the biased local Pauli coefficient rule.

Define the following quantity: \[\widetilde{A}_{l,\mathrm{LP}}(\rho) := B_{l,\mathrm{LP}}+\|R_lm\|_2 . \label{eq:local95pauli95Al95LP95def}\tag{175}\] By Lemma 6, this is exactly the upper bound \(\widetilde{A}_l(\rho)\) of the radius \(A_l(\rho)\) defined in Definition 1 for the present biased local Pauli protocol.

This selected-radius calculation also identifies the structural features of the local Pauli protocol that make finite-sample selected covariance estimation effective. First, locality ensures that a bounded-weight Pauli coordinate depends only on a bounded number of local measurement choices. Second, the bounded-weight assumption prevents the inverse-probability factors in the reconstructed coefficients from growing with the total number of qubits. Third, local basis compatibility determines which selected Pauli coordinates can be simultaneously active under a single basis pattern, and therefore controls the size of the selected one-shot vector. Finally, the inverse-probability structure shows explicitly how small local basis probabilities amplify both the covariance entries and the selected radius. Thus the dimension-independent selected covariance guarantee is not a black-box consequence of matrix concentration alone. It follows because the local Pauli measurement structure turns the abstract selected-radius condition into a concrete bounded quantity for bounded-size, bounded-weight selected coordinate sets. This point will be useful in the comparison below: the abstract finite-sample theorem applies to any shadow protocol once the selected radius is controlled, but the availability of such control is a protocol-dependent structural feature.

4.3 Comparison with global Clifford shadows↩︎

The preceding local Pauli result illustrates how locality and bounded observable weight can keep the selected one-shot radius independent of the number of qubits. We now contrast this behavior with global Clifford shadows. The purpose of this comparison is not to derive a full covariance formula for global Clifford shadows, but only to show how the selected-radius quantity \(B_l = \operatorname*{ess\,sup}\|R_l\hat{x}\|_2\) behaves when the measurement mechanism is changed from locally chosen Pauli bases to a global Clifford measurement.

The comparison is made within the same abstract framework as before. The target observables are unchanged: we still take the output coordinates to be the non-identity Pauli strings \[\mathcal{O}=\mathcal{P}_k^\times, \qquad p=|\mathcal{P}_k^\times|=d^2-1, \qquad d=2^k.\] Only the measurement mechanism is changed. In the global Clifford protocol, a Clifford unitary \(C\) is drawn, the state is measured in the computational basis after applying \(C\), and a measurement outcome \(b\in\{0,1\}^k\) gives the rank-one stabilizer snapshot \[\hat{\sigma}_{C,b} := C^\dagger |b\rangle\!\langle b| C . \label{eq:global95clifford95raw95snapshot}\tag{188}\] The reconstructed shadow snapshot is \[\hat{\rho}_{C,b} := \mathcal{M}_{\mathrm{Cl}}^{-1}(\hat{\sigma}_{C,b}), \label{eq:global95clifford95reconstructed95snapshot}\tag{189}\] and the Pauli output coordinate is \[\hat{x}_P := \operatorname{Tr}(\hat{\rho}_{C,b}P), \qquad P\in\mathcal{P}_k^\times . \label{eq:global95clifford95pauli95output95coordinate}\tag{190}\]

The key difference from the local Pauli protocol is the action of the shadow channel on the non-identity Pauli sector. For the global Clifford ensemble, the shadow channel acts as a scalar on the traceless Pauli sector. Equivalently, using the standard Clifford-shadow inverse map, or the Pauli-invariant reconstruction formula of Bu–Koh–Garcia–Jaffe [13], one has \(\mathcal{M}_{\mathrm{Cl}}[P] = \frac{1}{d+1}P\) for \(P\in\mathcal{P}_k^\times\), and hence \[\mathcal{M}_{\mathrm{Cl}}^{-1}[P] = (d+1)P. \label{eq:global95clifford95inverse95scalar95action}\tag{191}\] Therefore, for a single global Clifford snapshot \(\hat{\sigma}=\hat{\sigma}_{C,b}\), \[\hat{x}_P = \operatorname{Tr}\!\left( \mathcal{M}_{\mathrm{Cl}}^{-1}(\hat{\sigma})P \right) = (d+1)\operatorname{Tr}(\hat{\sigma} P). \label{eq:global95clifford95coefficient95rule}\tag{192}\]

For a stabilizer snapshot \(\hat{\sigma}\), define its non-identity stabilizer support by \[S(\hat{\sigma}) := \{P\in\mathcal{P}_k^\times: \operatorname{Tr}(\hat{\sigma} P)\neq 0 \}. \label{eq:global95clifford95stabilizer95support95def}\tag{193}\] For a rank-one stabilizer state, the Pauli expectation \(\operatorname{Tr}(\hat{\sigma} P)\) is equal to \(\pm1\) for \(P\in S(\hat{\sigma})\), and is zero otherwise. Thus, we have \[\hat{x}_P^2 = (d+1)^2 \mathbf{1}_{\{P\in S(\hat{\sigma})\}}. \label{eq:global95clifford95coefficient95square}\tag{194}\]

Let \(J_l=\{P_1,\ldots,P_l\}\subset\mathcal{P}_k^\times\) be a fixed selected set of distinct non-identity Pauli strings, and let \(R_l\) be the corresponding selection matrix. Define \[\kappa_{\mathrm{Cl}}(J_l) := \max_{\hat{\sigma}} |J_l\cap S(\hat{\sigma})|, \label{eq:global95clifford95kappa95selected95def}\tag{195}\] where the maximum is over rank-one stabilizer snapshots that can occur in the global Clifford protocol.

4.4 Exact covariance formula for qubit local Pauli shadows↩︎

We now record an exact covariance formula for the same biased local Pauli protocol. This calculation is not needed for the abstract finite-sample theorem or for the selected-radius bound above. Its role is instead structural: it identifies how local Pauli compatibility and local basis-selection probabilities determine the raw second moments of the reconstructed Pauli coefficients and, after subtracting the products of the corresponding means, the entries of the covariance matrix.

For the uniform local Pauli protocol, the underlying second-moment mechanism already appears in the observable-wise variance analysis of Huang–Kueng–Preskill [11]. We first recall the corresponding full covariance-matrix form in the uniform case, and then extend it to biased basis probabilities.

The purpose of this section is to extend this result to the biased local Pauli shadows. Under our setting of the biased local Pauli shadows given in 160 , the above proposition is generalized as follows.

The proof of Theorem 4 proceeds by inserting the coefficient rule for \(\hat{x}_P\) and \(\hat{x}_Q\), separating the incompatible case \(P\not\sim Q\) from the compatible case \(P\sim Q\), identifying the unique common basis event in the latter case, and then subtracting \(m_Pm_Q\) to pass from the raw second moment to the covariance.

5 General Local Product Shadows↩︎

The local Pauli analysis above uses special algebraic features of qubit Pauli measurements. However, the bounded-radius mechanism itself is not restricted to Pauli measurements. In this section we record a more general local product-shadow setting in which the selected observables have finite support. The main point is that, when both the measurement and the reconstruction are local product constructions, a finite-weight observable only sees the local snapshots on its support. Consequently, the selected one-shot radius is controlled by the selected support sizes and local reconstruction coefficients, not by the total number \(k\) of tensor factors.

5.1 Local product measurement model and finite-weight observables↩︎

We set the physical system as \[\mathcal{H} = \bigotimes_{j=1}^k \mathcal{H}_j, \qquad \mathcal{H}_j\simeq \mathbb{C}^t , \label{eq:general95local95product95hilbert95space}\tag{198}\] where the local dimension \(t\) is fixed. At each site \(j\), let \(\mathcal{U}_j\) be a finite set of local measurement settings. A random local setting \(\hat{u}_j\in\mathcal{U}_j\) is drawn with probability \(\mathbb{P}(\hat{u}_j=u_j)=q_{j,u_j}.\) We assume that the local settings are chosen independently, so that we have \(q(u)= \prod_{j=1}^k q_{j,u_j}\) for \(u=(u_1,\ldots,u_k)\).

For each \(u_j\in\mathcal{U}_j\), let \(M^{(j)}_{u_j} = \{M^{(j)}_{b_j|u_j}\}_{b_j\in\mathcal{B}_{j,u_j}}\) be a POVM on \(\mathcal{H}_j\). Conditional on \(u=(u_1,\ldots,u_k)\), the product POVM is defined as \[M_{b|u} := \bigotimes_{j=1}^k M^{(j)}_{b_j|u_j}, \qquad b=(b_1,\ldots,b_k). \label{eq:general95local95product95povm95element}\tag{199}\] Thus one shot produces random data \((\hat{u},\hat{b}) = ((\hat{u}_1,\ldots,\hat{u}_k),(\hat{b}_1,\ldots,\hat{b}_k)).\)

For each site \(j\), define the local shadow channel by \[\mathcal{M}_j[X_j] := \sum_{u_j\in\mathcal{U}_j} q_{j,u_j} \sum_{b_j\in\mathcal{B}_{j,u_j}} M^{(j)}_{b_j|u_j} \operatorname{Tr} \left( X_j M^{(j)}_{b_j|u_j} \right). \label{eq:general95local95product95local95shadow95channel}\tag{200}\] We assume that each \(\mathcal{M}_j\) is invertible on a prescribed local operator subspace \(\mathcal{V}_j\). We also assume that all local POVM elements \(M^{(j)}_{b_j|u_j}\) belong to \(\mathcal{V}_j\), so that \(\mathcal{M}_j^{-1}(M^{(j)}_{b_j|u_j})\) is well defined.

The global shadow channel is defined, as in Section 2.1 has the form \[\mathcal{M}[X] = \sum_{u\in\mathcal{U}} q(u) \sum_{b\in\mathcal{B}_u} M_{b|u} \operatorname{Tr} \left( X M_{b|u} \right). \label{eq:general95local95product95global95shadow95channel}\tag{201}\] In the present local product setting, \(q(u)=\prod_j q_{j,u_j}\) and \(M_{b|u}=\bigotimes_j M^{(j)}_{b_j|u_j}\). Throughout this section, we assume also that each local shadow channel \(\mathcal{M}_j\) is invertible on \(\mathcal{L}(\mathcal{H}_j)\). The next lemma shows that these two product structures imply the product form of the global reconstruction.

The factorization formula 205 in Lemma 7 is the general local product analogue of the local Pauli snapshot formula. We also assume throughout this section that the local snapshots are normalized, \[\operatorname{Tr}(\hat{\rho}_j)=1 \qquad \text{almost surely for every }j. \label{eq:general95local95product95trace95one95snapshot}\tag{206}\] This assumption ensures that identity tensor factors do not contribute to the output coordinate.

We now specify the physical quantities whose expectation values are studied. Let \(O_1,\ldots,O_l\) be selected product observables of the form \[O_a = \bigotimes_{j=1}^k O_{a,j}, \qquad a=1,\ldots,l, \label{eq:general95local95product95selected95observable}\tag{207}\] where \(O_{a,j}\) is a Hermitian operator on \(\mathcal{H}_j\), and \(O_{a,j}=I_j\) for all but finitely many sites. Remember the definitions of the support and the weight as \[\mathop{\mathrm{supp}}(O_a) = \{j:\,O_{a,j}\neq I_j\}, \quad \mathop{\mathrm{wt}}(O_a)=|\mathop{\mathrm{supp}}(O_a)|. \label{eq:general95local95product95support95weight}\tag{208}\] We assume that the selected observables have finite weight, and later we will impose a uniform bound \(\mathop{\mathrm{wt}}(O_a)\le w_0.\)

For each selected observable, define the shadow-output coordinate \[\hat{x}_a := \operatorname{Tr}(\hat{\rho} O_a), \qquad a=1,\ldots,l. \label{eq:general95local95product95shadow95output95coordinate}\tag{209}\] The selected output vector is \[R_l\hat{x} = (\hat{x}_1,\ldots,\hat{x}_l)^\top \in\mathbb{R}^l. \label{eq:general95local95product95selected95output95vector}\tag{210}\] Equivalently, in this section we may regard the selected observables \(O_1,\ldots,O_l\) as the full observable family and take \(R_l=I_l\).

For a local operator \(A\) on \(\mathcal{H}_j\), define the local reconstruction coefficient \[\eta_j(A;\hat{u}_j,\hat{b}_j) := \operatorname{Tr} \left[ \mathcal{M}_j^{-1} \left( M^{(j)}_{\hat{b}_j|\hat{u}_j} \right) A \right] = \operatorname{Tr}(\hat{\rho}_j A). \label{eq:general95local95product95eta95def}\tag{211}\] The following lemma is the general local product analogue of the local Pauli coefficient rule.

In the biased local Pauli protocol, the local coefficient \(\eta_j\) reduces to \(p_{j,P_j}^{-1} \mathbf{1}_{\{\hat{r}_j=P_j\}}\hat{s}_j\) for a non-identity local Pauli \(P_j\). Thus the Pauli compatibility indicator is a special feature of the local Pauli measurement basis. In the general local product setting, the measurement and the local observable need not belong to the same basis, and the corresponding dependence is absorbed into the coefficient \(\eta_j\).

5.2 Selected one-shot radius and centered radius bounds↩︎

We now bound the selected one-shot radius \(B_l = \operatorname*{ess\,sup}\|R_l\hat{x}\|_2\) and the upper bound \(\widetilde{A}_l(\rho)=B_l+\|R_lm\|_2\) of the centered radius \(A_l(\rho)\). The key point is that each coordinate \(\hat{x}_a\) depends only on the local snapshots on \(\mathop{\mathrm{supp}}(O_a)\).

For a local operator \(A\) on \(\mathcal{H}_j\), define the one-site coefficient radius \[\Gamma_j(A) := \operatorname*{ess\,sup} \left| \eta_j(A;\hat{u}_j,\hat{b}_j) \right|. \label{eq:general95local95product95gamma95def}\tag{214}\] For a selected product observable \(O_a\), set \[L_a := \prod_{j\in\mathop{\mathrm{supp}}(O_a)} \Gamma_j(O_{a,j}). \label{eq:general95local95product95La95def}\tag{215}\] Then Lemma 8 immediately gives \[|\hat{x}_a|\le L_a \qquad \text{almost surely}. \label{eq:general95local95product95coordinate95bound}\tag{216}\]

The preceding proposition gives a radius bound in terms of the actual local reconstruction coefficients appearing in the selected observables. The next corollary records a simple uniform version.

A more concrete sufficient condition can be stated in terms of the operator-norm size of the local reconstructed snapshots. Suppose that \[\|\hat{\rho}_j\|_{\mathrm{op}}\le G_0 \qquad \text{almost surely for all }j, \label{eq:general95local95product95local95snapshot95op95bound}\tag{222}\] and that the non-identity local factors of the selected observables satisfy \[\|O_{a,j}\|_{\mathrm{op}}\le 1 \qquad \text{for all }a,\;j\in\mathop{\mathrm{supp}}(O_a). \label{eq:general95local95product95local95observable95op95bound}\tag{223}\] Since \(\dim\mathcal{H}_j=t\), we have \[\|O_{a,j}\|_1\le t\|O_{a,j}\|_{\mathrm{op}}\le t.\] Therefore, \[\begin{align} \Gamma_j(O_{a,j}) &= \operatorname*{ess\,sup} |\operatorname{Tr}(\hat{\rho}_jO_{a,j})| \notag\\ &\le \operatorname*{ess\,sup} \|\hat{\rho}_j\|_{\mathrm{op}}\, \|O_{a,j}\|_1 \le G_0t . \label{eq:general95local95product95G0t95gamma95bound} \end{align}\tag{224}\] Thus, under the finite-weight and bounded-size assumptions of Corollary 3, the choice \(\Gamma_0=G_0t\) implies \[\begin{align} B_l & \le \sqrt{l_0}\,(G_0t)^{w_0}, \tag{225} \\ \widetilde{A}_l(\rho) & \le \sqrt{l_0}\, \bigl( (G_0t)^{w_0}+M_0 \bigr). \tag{226} \end{align}\]

Consequently, for finite-weight selected observables under a local product shadow protocol with uniformly bounded local reconstruction coefficients, the bounded-radius hypothesis in Theorem 1 holds with constants independent of \(k\). The selected sample-centered covariance can therefore be estimated to constant operator-norm accuracy with a sample size independent of the total number of tensor factors.

6 Discussion↩︎

This paper has developed a covariance-matrix viewpoint for classical-shadow outputs. Instead of organizing the analysis only around individual observables, marginal variances, or many-observable prediction bounds, we studied selected covariance matrices of reconstructed shadow-output vectors. The main finite-sample theorem applies to fixed selected coordinates for arbitrary shadow protocols and gives an operator-norm error bound for the selected sample-centered empirical covariance. The bound contains protocol-dependent constants that measure how large the selected reconstructed vector can be in a single measurement round. When these constants and the number of selected coordinates are independent of the ambient system size, the required sample size is also independent of the ambient dimension.

The main structural message is that such protocol-dependent constants can be controlled in local measurement settings. For general local product shadows, finite-weight product observables are controlled by their support sizes and local reconstruction coefficients, rather than by the total number of tensor factors. Biased local Pauli shadows provide a fully explicit instance: the relevant constants are written directly in terms of selected Pauli supports and local basis-selection probabilities. A comparison with global Clifford shadows shows that the same local behavior should not be expected without additional protocol structure. The exact covariance formula for biased local Pauli shadows gives further structural information, but it is separate from the finite-sample estimation argument.

6.1 What the covariance viewpoint adds↩︎

The main conceptual point is that the covariance matrix \(\Sigma(\rho)=\operatorname{Cov}_\rho(\hat{x})\) contains information that is not visible from marginal coordinate variances alone. Its diagonal entries describe the variances of individual reconstructed output coordinates, while its off-diagonal entries describe statistical couplings between different coordinates produced by the same shadow protocol. These off-diagonal entries are covariances of random post-processed shadow outputs, and therefore describe the joint fluctuation structure of the measurement-and-reconstruction procedure. In the Pauli specialization, this statistical covariance should not be confused with physical Pauli product correlations.

This distinction matters whenever one studies more than one shadow-output coordinate at a time. A multi-term observable, or a finite family of observables, does not depend only on the list of marginal variances. Its estimation variance also depends on how the selected output coordinates fluctuate together. The covariance matrix is the object that records this collective statistical structure.

The selected covariance matrix \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\) gives a natural finite-dimensional object associated with a fixed selected coordinate set. Its eigenvalues describe the largest and smallest fluctuation scales among normalized linear combinations of the selected coordinates. Its eigenvectors identify the corresponding principal fluctuation directions inside the selected coordinate set. Thus selected covariance spectra provide a diagnostic that is different from, but complementary to, observable-wise shadow-norm guarantees.

The covariance viewpoint also connects with reconstruction error whenever the chosen output coordinates form a coordinate representation of the reconstructed object. In such settings, traces or selected traces of \(\Sigma(\rho)\) quantify total fluctuation within the chosen coordinate system. More generally, the same covariance matrix organizes single-coordinate variances, off-diagonal statistical couplings, selected spectra, and fluctuations of selected reconstructed coordinates.

6.2 Why the finite-sample theory is selected↩︎

The finite-sample theory in this paper is deliberately formulated for fixed selected coordinate sets. This is not only a technical convenience, but also the natural scale at which one can expect useful operator-norm guarantees. The full shadow-output covariance matrix may have a very large ambient dimension \(p\), and estimating the entire matrix in operator norm would be a much stronger and substantially different task. Instead, we fix a selected coordinate set \(J_l=\{\alpha_1,\ldots,\alpha_l\}\) and study the compressed covariance \(\Sigma_l(\rho)=R_l\Sigma(\rho)R_l^\top\).

This selected formulation has two advantages. First, it matches the statistical question of interest when only a prescribed family of output coordinates is relevant. In that case, the eigenvalues of \(\Sigma_l(\rho)\) describe the fluctuation scales of linear combinations inside that selected family, and the corresponding spectral projectors identify stable or unstable selected directions. Second, the finite-sample concentration is performed after the compression by \(R_l\). Hence the dimension entering the probability bound is the selected dimension \(l\), rather than the ambient output dimension \(p\).

This is why the fixed-selection assumption is important. The selection matrix \(R_l\) is chosen independently of the data. Under this assumption, one may treat the selected covariance as an ordinary \(l\times l\) covariance matrix and apply matrix concentration directly in the selected space. If the coordinate set were chosen after looking at the data, additional arguments would be needed, such as uniform concentration over a larger class of possible selections, sample splitting, or other post-selection controls. Those questions are outside the scope of the present paper.

The selected viewpoint should therefore be read as a finite-sample localization of covariance estimation. It does not claim to reconstruct the entire covariance structure of the full shadow-output space. Rather, it shows that once a coordinate family has been fixed, the covariance spectrum within that family can be estimated with explicit nonasymptotic guarantees. This is the covariance analogue of focusing on a fixed family of observables in ordinary classical-shadow prediction, but with the empirical target changed from expectation values to the selected covariance matrix itself.

6.3 Interpretation of the constant selected-error regime↩︎

The dimension-independent selected regime is the main finite-sample consequence of the selected formulation. The general theorem says that, for any chosen shadow protocol, selected covariance estimation has a sample complexity controlled by the selected dimension and by protocol-dependent constants in the error bound. If these constants remain bounded independently of the ambient system size, then constant operator-norm accuracy for the selected covariance requires a number of samples independent of the ambient dimension.

General local product shadows provide one mechanism for this behavior. When the selected observables have uniformly bounded weight and the relevant local reconstruction coefficients are uniformly bounded, the constants in the finite-sample bound are controlled independently of the number of tensor factors. This explains why dimension-independent selected covariance estimation is not tied to Pauli measurements: it follows from local product structure, finite support, and bounded local reconstruction.

The biased local Pauli protocol gives a fully explicit instance of this mechanism. If the selected set size is bounded by \(l_0\), the selected Pauli weights are bounded by \(w_0\), and the local basis probabilities satisfy \(p_{j,r}\ge p_0>0\), then the finite-sample bound depends only on \(l_0\), \(w_0\), and \(p_0\), not on the total number of qubits. Consequently, the sample size needed to estimate the selected sample-centered covariance to constant operator-norm accuracy can be chosen independently of the total number of qubits.

The comparison with global Clifford shadows separates two issues. The finite-sample theorem is a general statement about selected covariance estimation, but whether its constants remain independent of system size is a protocol-specific question. Local product protocols, including biased local Pauli shadows, give positive examples under finite-support and bounded-coefficient assumptions. Global Clifford shadows show that such behavior should not be expected without additional local or protocol-specific structure.

6.4 Local Pauli covariance formulas as structural information↩︎

The exact local Pauli covariance formula should be read as structural information that complements the finite-sample results. The finite-sample selected covariance theorem does not require an entrywise formula for \(\Sigma(\rho)\). It only requires the protocol-dependent constants in the concentration bound to be controlled. Nevertheless, in the biased local Pauli protocol, the covariance entries can be computed explicitly in terms of Pauli expectations of the state and the local basis-selection probabilities.

The compatibility relation \(P\sim Q\), the cancellation string \(P\ominus Q\), and the overlap factor \(\beta(P,Q)=\prod_{j\in\mathop{\mathrm{supp}}(P)\cap\mathop{\mathrm{supp}}(Q)}p_{j,P_j}^{-1}\) describe how the local Pauli measurement design determines the raw second moments of the reconstructed Pauli coefficients. In particular, Pauli compatibility determines which off-diagonal raw second moments can be nonzero, while inverse-probability overlap factors determine the size of the corresponding raw second-moment contributions. The covariance entries are obtained from these raw second moments by subtracting the products of the corresponding means.

This formula is useful because it shows how local basis bias changes the covariance matrix of the reconstructed Pauli coefficient vector. The bias does not merely rescale marginal variances. It also changes off-diagonal statistical couplings through the same inverse-probability mechanism. Thus the biased local Pauli protocol provides an explicit model in which the covariance matrix can be analyzed entry by entry, while the finite-sample theory explains how selected parts of that covariance matrix can be estimated from data.

This analytic formula is also a reminder that selected covariance estimation and closed-form covariance calculation are different tasks. The selected finite-sample theorem applies even when no explicit formula for \(\Sigma(\rho)\) is available, provided the protocol-dependent constants in the error bound can be controlled. Conversely, an explicit covariance formula can reveal detailed structure, but it does not by itself give finite-sample operator-norm recovery. In the local Pauli setting both ingredients are available: the finite-sample bound gives the estimation guarantee, and the exact covariance formula explains the covariance structure being estimated.

6.5 Scope and outlook↩︎

The results of this paper should be read as selected covariance results. They do not assert operator-norm recovery of the full ambient covariance matrix \(\Sigma(\rho)\), and the selected coordinate set is assumed to be fixed independently of the data. Adaptive or post-selected coordinate choices would require additional tools, such as sample splitting or uniform concentration over candidate selections.

The dimension-independent regime also depends on protocol-specific constants. The finite-sample theorem applies to arbitrary shadow protocols, but it does not guarantee that these constants are small. General local product shadows and biased local Pauli shadows provide positive examples under finite-support and bounded-coefficient assumptions, whereas the global Clifford comparison shows that the corresponding constants can grow with the Hilbert-space dimension.

Finally, the exact covariance formula derived for biased local Pauli shadows is a protocol-specific structural result. It is separate from the finite-sample theorem, which only requires the relevant finite-sample constants to be controlled. Natural future directions include sharper bounds for structured selected sets, adaptive selected covariance estimation, extensions to other shadow protocols, and robust covariance estimation under noise, drift, or corrupted samples.

7 Proof of 34↩︎

Indeed, using \(\hat{x}_i-\bar{\hat{x}}_N = (\hat{x}_i-m)-(\bar{\hat{x}}_N-m),\) we expand the sample-centered covariance as \[\begin{align} & \widehat\Sigma_N^{\mathrm{sc}} \notag\\ =& \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-\bar{\hat{x}}_N)(\hat{x}_i-\bar{\hat{x}}_N)^\top \notag\\ =& \frac{1}{N} \sum_{i=1}^N \bigl((\hat{x}_i-m)-(\bar{\hat{x}}_N-m)\bigr) \bigl((\hat{x}_i-m)-(\bar{\hat{x}}_N-m)\bigr)^\top . \end{align}\] Expanding the product gives \[\begin{align} \widehat\Sigma_N^{\mathrm{sc}} &= \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-m)(\hat{x}_i-m)^\top \notag\\ &\quad - \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-m)(\bar{\hat{x}}_N-m)^\top \notag\\ &\quad - \frac{1}{N} \sum_{i=1}^N (\bar{\hat{x}}_N-m)(\hat{x}_i-m)^\top \notag\\ &\quad + \frac{1}{N} \sum_{i=1}^N (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top . \end{align}\] Since \(\frac{1}{N}\sum_{i=1}^N(\hat{x}_i-m) = \bar{\hat{x}}_N-m,\) the two cross terms become \(- (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top.\) The last term is \[\frac{1}{N} \sum_{i=1}^N (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top = (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top .\] Therefore, we obtain \[\begin{align} \widehat\Sigma_N^{\mathrm{sc}} &= \frac{1}{N} \sum_{i=1}^N (\hat{x}_i-m)(\hat{x}_i-m)^\top - (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top \notag\\ &= \widehat\Sigma_N^{\mathrm{tc}} - (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top , \end{align}\] which implies \(\widehat\Sigma_N^{\mathrm{sc}} - \widehat\Sigma_N^{\mathrm{tc}} = - (\bar{\hat{x}}_N-m)(\bar{\hat{x}}_N-m)^\top .\)

8 Proof of 163↩︎

We briefly derive the reconstruction formula. For a single qubit at site \(j\), the biased Pauli measurement channel acts on an operator \(A\) as \[\begin{align} \mathcal{M}_j[A] & = \sum_{r\in\{X,Y,Z\}} p_{j,r} \sum_{s=\pm1} M^{(j)}_{s|r}\, \operatorname{Tr}\!\left(A M^{(j)}_{s|r}\right),\\ M^{(j)}_{s|r} &= \frac{1}{2}(I+s r). \end{align}\] Since \(\operatorname{Tr}(I M^{(j)}_{s|r})=1\), and \(\operatorname{Tr}(t M^{(j)}_{s|r})=s\,\mathbf{1}_{\{t=r\}}\) for \(t\in\{X,Y,Z\}\), we have \[\mathcal{M}_j[I]=I, \qquad \mathcal{M}_j[t]=p_{j,t}\,t, \quad t\in\{X,Y,Z\}. \label{eq:single95qubit95biased95pauli95channel95diagonal95action}\tag{227}\] Thus \(\mathcal{M}_j\) is diagonal in the single-qubit Pauli basis, and its inverse is \[\mathcal{M}_j^{-1}[I]=I, \qquad \mathcal{M}_j^{-1}[t]=p_{j,t}^{-1}t, \quad t\in\{X,Y,Z\}. \label{eq:single95qubit95biased95pauli95inverse95action}\tag{228}\] Applying this inverse map to the actually observed rank-one projector \(M^{(j)}_{\hat{s}_j|\hat{r}_j}=\frac{1}{2}(I+\hat{s}_j\hat{r}_j)\) gives the single-qubit reconstructed snapshot \[\hat{\rho}_j^{(1)} = \mathcal{M}_j^{-1}\!\left[ M^{(j)}_{\hat{s}_j|\hat{r}_j} \right] = \frac{1}{2} \left( I+\frac{\hat{s}_j}{p_{j,\hat{r}_j}}\hat{r}_j \right). \label{eq:single95qubit95biased95pauli95reconstructed95snapshot}\tag{229}\] Because both the basis choice and the POVM are product over sites, the \(k\)-qubit shadow channel is the tensor product \(\mathcal{M}=\mathcal{M}_1\otimes\cdots\otimes\mathcal{M}_k\), and hence its inverse is \(\mathcal{M}^{-1}=\mathcal{M}_1^{-1}\otimes\cdots\otimes\mathcal{M}_k^{-1}\). Therefore the full reconstructed shadow snapshot is \[\hat{\rho} = \bigotimes_{j=1}^k \hat{\rho}_j^{(1)} = \bigotimes_{j=1}^k \frac{1}{2} \left( I+\frac{\hat{s}_j}{p_{j,\hat{r}_j}}\hat{r}_j \right), \label{eq:biased95local95pauli95reconstructed95snapshot}\tag{230}\] which yields 163 .

9 Proof of Theorem 4↩︎

This appendix extends the exact local-Pauli covariance calculation (Proposition 4) by Huang–Kueng–Preskill [11] to the biased measurement design of Section 4. The proof follows the same overall logic as in the uniform case, but with the deterministic coefficient rule now carrying explicit inverse-probability weights. We first derive the weighted single-shot coefficient formula, then compute the corresponding raw second moments for compatible and incompatible pairs, and finally obtain the exact covariance expression in terms of the overlap bias factor. The result is the closed-form covariance formula for the biased local design stated in Theorem 4.

In this appendix we prove Theorem 4. Throughout, we work with the \(k\)-qubit biased local Pauli classical-shadow protocol described in Section 4.

Acknowledgments↩︎

The author was supported in part by the General R&D Projects of 1+1+1 CUHK-CUHK(SZ)-GDST Joint Collaboration Fund (Grant No. GRDP2025-022), the Guangdong Provincial Quantum Science Strategic Initiative (Grant No. GDZX2505003), and the Shenzhen International Quantum Academy (Grant No. SIQA2025KFKT07). Large language model tools were used as auxiliary aids in preparing this manuscript, including assistance with exposition, literature search, and exploratory discussions of possible approaches to the research problem. The manuscript was written under the author’s direction, and the author is solely responsible for all mathematical content, proofs, references, and conclusions.

References↩︎