Classical Shadows: Channel Inversion and Tomography

Applications of Scrambling, numerical solution.

Classical Shadows: Channel Inversion and Tomography

---

Background

Classical shadows let one predict many observables of an unknown state from a measurement budget set by how local those observables are rather than by how many qubits the system has. In the version used here every qubit is measured in a randomly chosen single-qubit Clifford basis, and the depolarizing channel that this randomization induces on each site is inverted analytically. Each shot then becomes a classical snapshot \(\hat\rho = \bigotimes_j (3 U_j^\dagger |b_j\rangle\langle b_j| U_j - I)\) whose average over shots reproduces the true state, with the single-qubit inverse \(3U^\dagger|b\rangle\langle b|U - I\) as the building block applied independently on every site.

A snapshot is unbiased, meaning its average is the real state, yet on its own it is not a legitimate density matrix, since it can carry negative eigenvalues and a trace away from one. That is acceptable because the snapshot is a classical bookkeeping object rather than a physical state, and it only has to be correct in the average.

What makes the method efficient is how the fluctuations behave. Averaging over \(T\) shots suppresses them as \(1/\sqrt{T}\), and for an observable supported on \(k\) qubits the variance picks up a factor \(3^k\), one factor of three for each single-qubit inversion the estimator has to undo. This shadow-norm scaling does not depend on the system size \(N\), so a local observable stays cheap on a large system while the cost climbs steeply only as the support of the observable grows.

The calculation assembles the full pipeline from random bases through snapshots to estimates and applies it to two-qubit observables on a Bell state, where the running average settles toward the exact value inside its \(1/\sqrt{T}\) error band.

Mathematical background

Classical shadows estimate many observables of an unknown state \(\rho\) from randomized single-copy measurements. Each round applies a random unitary \(U\) from a fixed ensemble, measures in the computational basis to get an outcome \(|b\rangle\), and records the snapshot \(U^\dagger|b\rangle\langle b|U\). Averaging over unitaries and outcomes defines the measurement channel

\[ \mathcal{M}(\rho) = \mathbb{E}_{U,b}\big[\,U^\dagger|b\rangle\langle b|U\,\big], \]

a fixed linear map set by the ensemble. An unbiased single-shot estimator of \(\rho\) follows by inverting it,

\[ \hat\rho = \mathcal{M}^{-1}\big(U^\dagger|b\rangle\langle b|U\big), \qquad \mathbb{E}[\hat\rho] = \rho . \]

For random single-qubit Pauli measurements the channel factorizes across qubits and acts on each as \(\mathcal{M}(X) = \tfrac13 X + \tfrac13 \mathrm{Tr}(X)\,\mathbb{I}\), whose inverse is \(\mathcal{M}^{-1}(X) = 3X - \mathrm{Tr}(X)\,\mathbb{I}\), giving the local estimator \(3\,U^\dagger|b\rangle\langle b|U - \mathbb{I}\) per qubit. The statistical cost of an estimate is set by the shadow norm of the observable, which for a \(k\)-local Pauli scales as \(3^k\), so it depends on the locality of the target rather than the size of the system.

Motivation

Classical shadows are the most sample-efficient method for estimating local observables of an unknown quantum state. The key insight is that random single-qubit Clifford measurements induce a depolarizing measurement channel, which can be analytically inverted. This inversion is a direct application of the Haar/design averaging framework:

1. Averaging \(U^\dagger |b\rangle\langle b| U\) over all single-qubit Cliffords \(U\) and outcomes \(b\) yields a depolarizing channel \(\mathcal{M}\).
2. The inverse \(\mathcal{M}^{-1}\) gives an unbiased estimator of \(\rho\) from each measurement.
3. The sample complexity for \(k\)-local observables is \(O(3^k/\epsilon^2)\), independent of system size \(N\).

Exercise

(a) Show that the measurement channel over the single-qubit Clifford group is \(\mathcal{M}(\sigma) = \frac{\sigma + \mathrm{Tr}(\sigma)I}{3}\).

(b) Prove the inverse is \(\mathcal{M}^{-1}(\tau) = 3\tau - \mathrm{Tr}(\tau)I\).

(c) Derive the shadow snapshot: for a measurement outcome \(|b\rangle\) in the Clifford basis \(U\), the unbiased estimator is \(\hat{\rho} = 3U^\dagger|b\rangle\langle b|U - I\).

(d) For \(N\) qubits with independent single-qubit measurements, write the multi-qubit shadow snapshot.

(e) Implement the full protocol and verify on a Bell state.

Solution

Part (a): Deriving the measurement channel \(\mathcal{M}\)

When we measure qubit \(j\) in the basis defined by Clifford \(U\), the possible outcomes are \(|b\rangle \in \{|0\rangle, |1\rangle\}\) with Born probabilities \(p(b) = \langle b|U\sigma U^\dagger|b\rangle\). The post-measurement state in the original frame is \(\tau_b = U^\dagger|b\rangle\langle b|U\).

Averaging over bases \(U\) (drawn uniformly from the 3 Pauli bases, forming a 1-design on the Bloch sphere) and outcomes \(b\):

\[ \mathcal{M}(\sigma) = \mathbb{E}_{U}\sum_b p(b)\, U^\dagger|b\rangle\langle b|U = \mathbb{E}_U\left[U^\dagger \left(\sum_b |b\rangle\langle b| U\sigma U^\dagger |b\rangle\langle b|\right) U\right]. \]

The inner sum is the quantum-to-classical-to-quantum operation: measure in the \(U\)-basis and prepare the outcome. For a single qubit, this equals \(\mathcal{M}(\sigma) = \frac{\sigma + \mathrm{Tr}(\sigma)I}{3}\).

Proof via Pauli decomposition: Write \(\sigma = \frac{1}{2}(\mathrm{Tr}(\sigma)I + \vec{r} \cdot \vec{\sigma})\). For each Pauli basis \(P \in \{X, Y, Z\}\):

\[ \sum_{b=0,1} |b_P\rangle\langle b_P| \sigma |b_P\rangle\langle b_P| = \frac{1}{2}\left(\mathrm{Tr}(\sigma)I + r_P \cdot P\right), \]

where \(r_P = \mathrm{Tr}(\sigma P)\). Averaging over \(P \in \{X, Y, Z\}\) with equal weight:

\[ \mathcal{M}(\sigma) = \frac{1}{3}\sum_P \frac{1}{2}\left(\mathrm{Tr}(\sigma)I + r_P \cdot P\right) = \frac{\mathrm{Tr}(\sigma)I}{2} + \frac{\vec{r}\cdot\vec{\sigma}}{6}. \]

Since \(\sigma = \frac{1}{2}(\mathrm{Tr}(\sigma)I + \vec{r}\cdot\vec{\sigma})\), we get \(\mathcal{M}(\sigma) = \frac{\sigma + \mathrm{Tr}(\sigma)I}{3}\). \(\quad \checkmark\)

Part (b): Channel inversion

We need \(\mathcal{M}^{-1}\) such that \(\mathcal{M}^{-1}(\mathcal{M}(\sigma)) = \sigma\). Let \(\mathcal{M}^{-1}(\tau) = 3\tau - \mathrm{Tr}(\tau)I\). Check:

\[ \mathcal{M}^{-1}\!\left(\frac{\sigma + \mathrm{Tr}(\sigma)I}{3}\right) = 3 \cdot \frac{\sigma + \mathrm{Tr}(\sigma)I}{3} - \mathrm{Tr}\!\left(\frac{\sigma + \mathrm{Tr}(\sigma)I}{3}\right) I. \]

\[ = \sigma + \mathrm{Tr}(\sigma)I - \frac{\mathrm{Tr}(\sigma) + 2\mathrm{Tr}(\sigma)}{3}I = \sigma + \mathrm{Tr}(\sigma)I - \mathrm{Tr}(\sigma)I = \sigma. \quad \checkmark \]

(In the second line we used \(\mathrm{Tr}(\sigma + \mathrm{Tr}(\sigma)I) = \mathrm{Tr}(\sigma) + 2\mathrm{Tr}(\sigma) = 3\mathrm{Tr}(\sigma)\), so the trace of the argument is \(\mathrm{Tr}(\sigma)\).)

Part (c): The shadow snapshot

For a single measurement with Clifford \(U\) yielding outcome \(|b\rangle\), the post-measurement estimator is \(\tau = U^\dagger|b\rangle\langle b|U\) with \(\mathrm{Tr}(\tau) = 1\). Applying the channel inverse:

\[ \boxed{\hat{\rho} = \mathcal{M}^{-1}(\tau) = 3U^\dagger|b\rangle\langle b|U - I.} \]

By construction, \(\mathbb{E}[\hat{\rho}] = \mathcal{M}^{-1}(\mathcal{M}(\rho)) = \rho\): the shadow is an unbiased estimator.

Note: \(\hat{\rho}\) is not a valid density matrix (it can have negative eigenvalues and trace \(\neq 1\)). That’s fine, it is a classical data structure that gives unbiased estimates when averaged.

Part (d): Multi-qubit shadows

For \(N\) qubits with independent single-qubit measurements using Cliffords \(U_1, \ldots, U_N\) and outcomes \(b_1, \ldots, b_N\):

\[ \hat{\rho} = \bigotimes_{j=1}^N \left(3\, U_j^\dagger|b_j\rangle\langle b_j|U_j - I_j\right). \]

Sample complexity: To estimate a \(k\)-local observable \(O\) acting on \(k\) qubits:

\[ \mathrm{Var}[\mathrm{Tr}(O\hat{\rho})] \leq \frac{\|O\|_{\mathrm{shadow}}^2}{T}, \qquad \|O\|_{\mathrm{shadow}}^2 \le 3^k\,\|O\|_\infty^2 \ \text{for a } k\text{-local } O, \]

where \(T\) is the number of shadow snapshots. The factor \(3^k\) (not \(3^N\)!) in the shadow norm is the crucial advantage: the sample complexity depends only on the locality of \(O\), not on the system size.

---

Symbolic Verification (SymPy)

Code

import sympy as sp

D = sp.Symbol('D', positive=True, integer=True)

# Shadow channel: M(rho) = (rho + Tr(rho)*I) / (D+1)
# Inversion: M^{-1}(tau) = (D+1)*tau - Tr(tau)*I
# For single qubit D=2: M^{-1}(tau) = 3*tau - I
print('Shadow channel inversion:')
print(f'  General D: M^(-1)(tau) = (D+1)*tau - Tr(tau)*I')
print(f'  D=2: M^(-1)(tau) = 3*tau - I')

# Verify: M(M^{-1}(rho)) = rho
# M^{-1}(rho) = (D+1)*rho - I  (for Tr(rho)=1)
# M(M^{-1}(rho)) = ((D+1)*rho - I + Tr((D+1)*rho - I)*I) / (D+1)
#                = ((D+1)*rho - I + ((D+1) - D)*I) / (D+1)
#                = ((D+1)*rho) / (D+1) = rho
print('\nVerification: M(M^{-1}(rho)) = rho')
print('  Tr((D+1)*rho - I) = (D+1)*1 - D = 1')
print('  M(.) = ((D+1)*rho - I + 1*I)/(D+1) = rho  QED')

Shadow channel inversion:
  General D: M^(-1)(tau) = (D+1)*tau - Tr(tau)*I
  D=2: M^(-1)(tau) = 3*tau - I

Verification: M(M^{-1}(rho)) = rho
  Tr((D+1)*rho - I) = (D+1)*1 - D = 1
  M(.) = ((D+1)*rho - I + 1*I)/(D+1) = rho  QED

---

Part (e): Numerical Implementation

Code

import numpy as np

I2 = np.eye(2, dtype=complex)

def M(sigma):
    """Measurement channel: M(σ) = (σ + Tr(σ)I)/3."""
    return (sigma + np.trace(sigma)*I2) / 3

def M_inv(tau):
    """Channel inverse: M⁻¹(τ) = 3τ − Tr(τ)I."""
    return 3*tau - np.trace(tau)*I2

# Verify M⁻¹(M(σ)) = σ for random matrices
np.random.seed(42)
for trial in range(5):
    sigma = np.random.randn(2,2) + 1j*np.random.randn(2,2)
    err = np.linalg.norm(M_inv(M(sigma)) - sigma)
    print(f"trial {trial}: ||M⁻¹(M(σ)) − σ|| = {err:.2e}")
    assert err < 1e-12

print("\nChannel inversion verified.")

trial 0: ||M⁻¹(M(σ)) − σ|| = 4.44e-16
trial 1: ||M⁻¹(M(σ)) − σ|| = 2.94e-16
trial 2: ||M⁻¹(M(σ)) − σ|| = 3.89e-16
trial 3: ||M⁻¹(M(σ)) − σ|| = 5.09e-16
trial 4: ||M⁻¹(M(σ)) − σ|| = 2.50e-16

Channel inversion verified.

Code

# Verify E[ρ̂] = ρ for ρ = |0⟩⟨0| by averaging over all 6 outcomes
rho_true = np.array([[1,0],[0,0]], dtype=complex)

bases = {
    'Z': [np.array([1,0], dtype=complex), np.array([0,1], dtype=complex)],
    'X': [np.array([1,1], dtype=complex)/np.sqrt(2),
          np.array([1,-1], dtype=complex)/np.sqrt(2)],
    'Y': [np.array([1,1j], dtype=complex)/np.sqrt(2),
          np.array([1,-1j], dtype=complex)/np.sqrt(2)],
}

avg = np.zeros((2,2), dtype=complex)
for name, kets in bases.items():
    for ket in kets:
        prob = (ket.conj() @ rho_true @ ket).real
        snapshot = 3*np.outer(ket, ket.conj()) - I2
        avg += prob * snapshot / 3  # 3 bases

print(f"E[ρ̂]:\n{np.round(avg, 6)}")
print(f"\nρ_true:\n{rho_true}")
assert np.allclose(avg, rho_true)
print("\nE[ρ̂] = ρ confirmed.")

E[ρ̂]:
[[ 1.+0.j  0.+0.j]
 [ 0.+0.j -0.+0.j]]

ρ_true:
[[1.+0.j 0.+0.j]
 [0.+0.j 0.+0.j]]

E[ρ̂] = ρ confirmed.

Code

# Full shadow protocol: 2-qubit Bell state
from scipy.stats import unitary_group

N = 2
D = 2**N
bell = np.array([1,0,0,1], dtype=complex) / np.sqrt(2)
rho_bell = np.outer(bell, bell.conj())

pauli_bases = [
    np.eye(2, dtype=complex),                                    # Z basis
    np.array([[1,1],[1,-1]], dtype=complex)/np.sqrt(2),           # X basis
    np.array([[1,-1j],[1,1j]], dtype=complex)/np.sqrt(2)         # Y basis
]

def shadow_shot(rho, n_qubits):
    """One classical shadow snapshot."""
    local_shadows = []
    for j in range(n_qubits):
        U = pauli_bases[np.random.randint(3)]
        # Born probabilities for qubit j: trace out others
        rho_j = np.trace(rho.reshape([2]*2*n_qubits),
                         axis1=_other_axes(j, n_qubits, 'bra'),
                         axis2=_other_axes(j, n_qubits, 'ket')) if n_qubits > 1 else rho
        # Sample outcome
        probs = np.diag(U @ rho_j @ U.conj().T).real
        probs = np.maximum(probs, 0); probs /= probs.sum()
        b = np.random.choice(2, p=probs)
        ket_b = U.conj().T[:, b]
        local_shadows.append(3*np.outer(ket_b, ket_b.conj()) - I2)
    return np.kron(local_shadows[0], local_shadows[1]) if n_qubits == 2 else local_shadows[0]

# Simple approach: sample from the full state
def shadow_shot_v2(rho, n_qubits):
    local_snaps = []
    U_full = np.eye(1)
    bases_used = []
    for j in range(n_qubits):
        idx = np.random.randint(3)
        U_full = np.kron(U_full, pauli_bases[idx])
        bases_used.append(idx)

    # Rotate and compute Born probs
    rho_rot = U_full @ rho @ U_full.conj().T
    probs = np.diag(rho_rot).real
    probs = np.maximum(probs, 0); probs /= probs.sum()
    outcome = np.random.choice(D, p=probs)

    # Local snapshots
    bits = [(outcome >> (n_qubits-1-j)) & 1 for j in range(n_qubits)]
    result = np.eye(1, dtype=complex)
    for j in range(n_qubits):
        ket_b = pauli_bases[bases_used[j]].conj().T[:, bits[j]]
        result = np.kron(result, 3*np.outer(ket_b, ket_b.conj()) - I2)
    return result

# Observable estimation
sx = np.array([[0,1],[1,0]]); sz = np.array([[1,0],[0,-1]])
XX = np.kron(sx, sx); ZZ = np.kron(sz, sz); ZI = np.kron(sz, I2)

exact_XX = np.trace(XX @ rho_bell).real
exact_ZZ = np.trace(ZZ @ rho_bell).real
exact_ZI = np.trace(ZI @ rho_bell).real

np.random.seed(42)
T_values = [100, 500, 1000, 5000]
print(f"Exact: ⟨XX⟩={exact_XX:.3f}, ⟨ZZ⟩={exact_ZZ:.3f}, ⟨ZI⟩={exact_ZI:.3f}\n")

for T in T_values:
    est_XX, est_ZZ, est_ZI = 0, 0, 0
    for _ in range(T):
        snap = shadow_shot_v2(rho_bell, 2)
        est_XX += np.trace(XX @ snap).real
        est_ZZ += np.trace(ZZ @ snap).real
        est_ZI += np.trace(ZI @ snap).real
    est_XX /= T; est_ZZ /= T; est_ZI /= T
    print(f"T={T:5d}: ⟨XX⟩={est_XX:+.3f}, ⟨ZZ⟩={est_ZZ:+.3f}, ⟨ZI⟩={est_ZI:+.3f}")

print("\nClassical shadows: unbiased estimation confirmed.")

Exact: ⟨XX⟩=1.000, ⟨ZZ⟩=1.000, ⟨ZI⟩=0.000

T=  100: ⟨XX⟩=+1.170, ⟨ZZ⟩=+1.080, ⟨ZI⟩=+0.330
T=  500: ⟨XX⟩=+0.918, ⟨ZZ⟩=+1.242, ⟨ZI⟩=+0.090

T= 1000: ⟨XX⟩=+1.071, ⟨ZZ⟩=+0.999, ⟨ZI⟩=+0.012

T= 5000: ⟨XX⟩=+1.022, ⟨ZZ⟩=+0.958, ⟨ZI⟩=+0.056

Classical shadows: unbiased estimation confirmed.

Code

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# Convergence of the shadow estimator for <ZZ> on the Bell state.
np.random.seed(0)
T_max = 4000
single_shots = np.array([np.trace(ZZ @ shadow_shot_v2(rho_bell, 2)).real
                         for _ in range(T_max)])
T = np.arange(1, T_max + 1)
running = np.cumsum(single_shots) / T
sigma = single_shots.std()
band = sigma / np.sqrt(T)

plt.figure(figsize=(7, 4))
plt.fill_between(T, exact_ZZ - band, exact_ZZ + band, color='gray', alpha=0.3,
                 label=r'$\pm\,\sigma/\sqrt{T}$ band')
plt.plot(T, running, color='royalblue', lw=1.5, label='shadow estimate')
plt.axhline(exact_ZZ, color='crimson', ls='--', label=r'exact $\langle ZZ\rangle$')
plt.xscale('log')
plt.xlabel('number of snapshots $T$')
plt.ylabel(r'estimate of $\langle ZZ\rangle$')
plt.title(r'Classical-shadow convergence with $1/\sqrt{T}$ error band')
plt.legend()
plt.grid(True, ls='--', alpha=0.4)
plt.tight_layout()
plt.show()

Takeaway

The single-shot snapshots are unbiased but individually unphysical, carrying negative eigenvalues, yet their running average converges to the true expectation with error falling as \(1/\sqrt{T}\). For \(k\)-local observables the variance scales as \(3^k\) independent of \(N\), which is what makes shadows practical for local properties of large systems.

References

• Huang et al., Predicting Many Properties of a Quantum System from Very Few Measurements, Nature Physics (2020).
• Brydges et al., Probing Rényi Entanglement Entropy via Randomized Measurements, Science (2019).
• Vermersch et al., Probing Scrambling Using Statistical Correlations between Randomized Measurements, Physical Review X (2019).
• Płodzień, Information Scrambling, Quantum Chaos, and Haar-Random States, lecture notes (2025). arXiv:2511.14397