In shadow tomography, how is the state reconstructed from its shadows?

Question

I'm reading Huang et al. (2020) (nature physics), where the authors present a version of Aaronson's shadow tomography scheme as follows (see page 11 in the arXiv version):

We want to estimate a state $\rho$. We apply a number of random unitary evolutions, $\rho\mapsto U\rho U^\dagger$, picking $U$ from an ensemble $\mathcal U$. For each choice of unitary $U$, we perform a measurement, observing a state $|b\rangle$. We then apply the inverse evolution to this state, obtaining $U^\dagger|b\rangle$. On average, this procedure leaves us with the state $$\mathcal M(\rho) =\mathbb E_{U\sim\mathcal U}\Bigg[ \sum_b \underbrace{\langle b|U\rho U^\dagger|b\rangle}_{\text{prob of observing $|b\rangle$}} \!\!\!\!\!\!\! \overbrace{(U^\dagger |b\rangle\!\langle b|U)}^{\text{post-measurement state}} \Bigg].$$ The claim is then that $\mathcal M$ can be inverted to obtain, on average, the original state: $$ \hat\rho=\mathcal M^{-1}(U^\dagger |\hat b\rangle\!\langle \hat b| U) \,\,\text{ is such that } \,\, \mathbb E[\hat \rho]=\rho. $$

Is there an easy way to see how one would go in performing such inversion? The authors mention that we are thinking here of $\mathcal M$ as a linear map, so I suppose we represent $\rho$ as a vector in some operatorial basis of Hermitians, which is fine, but to then perform the inversion of this linear map we would need to characterise it, which in this case I think would mean to know the values of $\mathcal M(\sigma_i)$ for some complete basis of Hermitians $\{\sigma_i\}_i$.

Answer 1

Under the assumption that the ensemble $\mathcal{U}$ faithfully produces the Haar expectations at least to the second moment, the inversion can be performed as suggested in the last paragraph of the question:

Define:

$$\theta_b = U^{\dagger}|b\rangle\langle b| U$$

and replace the averaging over the ensemble by Haar averaging. (This step is is done only to compute the resulting quantity. In practice we average over the ensemble):

$$\mathcal{M}(\rho) = \int d\mu_{\small{HARR}} \sum_b \mathrm{Tr}(\rho \theta_b) \theta_b$$

($d\mu_{\small{HARR}} $ is normalized to a unit volume). Suppose that we need to measure the expectation of an operator $O$ in the state $\rho$. We measure this operator in the final state average the expectations over the ensemble and the computational basis projectors, we obtain the following result:

$$\mathrm{Tr}(\mathcal{M}(\rho) O) = \int d\mu_{\small{HARR}} \sum_b \mathrm{Tr}(\rho \theta_b) \mathrm{Tr}(O \theta_b)$$

Using the following identity:

$$\int d\mu_{\small{HARR}} \mathrm{Tr}(\rho \theta_b) (O \theta_b) = \frac{ \mathrm{Tr}(\rho) \mathrm{Tr}(O)+ \mathrm{Tr}(\rho O) }{N(N+1)}$$

($N$ is the dimension of the Hilbert space), thus we obtain:

$$\mathrm{Tr}(\mathcal{M}(\rho) O) = \frac{ \mathrm{Tr}(O)+ \mathrm{Tr}(\rho O) }{(N+1)}$$

From which the required expectation $\mathrm{Tr}(\rho O)$ can be computed.

Of course, we can repeat the procedure for a complete basis of the Lie algebra $\mathfrak{u}(N)$ for a full tomography of $\rho$.

Explanation of the integration formula

$\theta_b = U^{\dagger}|b\rangle\langle b| U$ is a one-dimensional projector can be written as:

$$\theta_b = |\Psi \rangle\langle\Psi|$$

where $ |\Psi\rangle$ is a unit vector

When $U$ uniformly scans the unitary group manifold, $ |\Psi\rangle$ will uniformly scan a $2N-1$ dimensional round sphere. This is because the metric of a round sphere is Haar. Thus, we need to evaluate an integral of the form:

$$\int_{S^{2N-1} }\langle\Psi| A |\Psi \rangle \langle\Psi| B |\Psi\rangle d\mu_{\small{ROUND}}(S^{2N-1}) $$

The integrand is a homogeneous polynomial of the components of the unit vector $ |\Psi\rangle = \sum_{a=1}^N \psi_a |a\rangle$:

$$\langle\Psi| A |\Psi \rangle \langle\Psi| B |\Psi\rangle = \sum_{a,b,c,d} A_{ab} B_{cd} \bar{\psi_a}\bar{\psi_c}\psi_b \psi_d$$

Integrals of homogeneous polynomials on the sphere can be evaluated as Gaussian integrals up to an overall normalization, which can be computed from a trivial particular example. In this case we can use the Wick's theorem:

The constant can be evaluated by taking $A=B=I_{N\times N}$, we get:

$$\mathrm{Const.} = \frac{1}{N(N+1)}$$

Answer 2

This is not a complete answer but a couple of comments that should help clarify these details:

(i) $\mathcal{M}$ is not invertible as a quantum channel but as a linear map; namely, even though the inverse exists, $\mathcal{M}^{-1}$ is not CP, even though it is a linear map. Assuming the input-output dimensions are the same, the only CP maps that are invertible are unitaries. This is the same problem that arises when studying dynamical semigroups for Markovian evolutions, they are described as a family $\{ e^{t \mathcal{L}} \}_{t \geq 0}$ where $\mathcal{L}$ is the Lindbladian; clearly $e^{-t \mathcal{L}}$ is the inverse, but it is not CP.

(ii) The authors assume ``tomographic completeness,'' that is, for each $\sigma \neq \rho, \exists U \in \mathcal{U} \text{ s.t. } \langle b| U \sigma U^{\dagger} | b \rangle \neq \langle b| U \rho U^{\dagger} | b \rangle$ which ensures that $\mathcal{M}^{-1}$ exists and is unique.

(iii) The classical shadow, $\hat{\rho}$ is not necessarily positive semidefinite since it originates from the application of a not CP map on some state.

(iv) $\mathcal{M}^{-1}$ is applied to the classically stored measurement outcomes, $U^{\dagger} | \hat{b} \rangle \langle \hat{b} | U$ and not the quantum states themselves.