5 Literature review↩︎

This section provides a comprehensive review of related literature, which consists of the following three parts.

Classical shadow in quantum observable estimation Classical shadow [7] is a protocol that allows experimenters to predict \(M\) observables of a quantum state \(\rho\) with only \(\log(M)\) copies of \(\rho\). The procedure only requires repeated single-copy measurements of \(\rho\). The terminology is inspired by shadow tomography, which Aaronson coined in [14]. One can also view the classical shadow estimator as a projected least-squares predictor [15]. Classical shadow has different procedures to perform random uniform scrambling on the state \(\rho\). In [7], the authors propose random Pauli measurement and random Clifford scrambling, the former of which can be implemented with a circuit of depth one, whereas the latter is significantly more involved in terms of circuit complexity. One can also consider other random scrambling protocols of intermediate complexity [16]–[18]. While our work only considers the random Pauli measurement setting, future work might also consider other random scrambling protocols.

Quantum state tomography Quantum state tomography (QST) is a fundamental task related to learning and engineering high-dimensional quantum systems in many-body physics. Specifically, it aims to reconstruct an unknown quantum state from a collection of experimental measurements. An extensive body of work has explored QST from many different aspects. For instance, one prominent line of work formulates QST as a matrix recovery problem [15], [19]–[34], where the density matrix is estimated via statistical techniques such as projected least squares, semidefinite programming, compressed sensing [35]–[38] and non-convex programming. Moreover, to alleviate the high computational cost of general QST, researchers have adopted simplifying assumptions by considering quantum states that can be well approximated by matrix product states with low bond dimension. This leads to the development of matrix product state tomography [3], [39]–[47] and variants such as multiscale entangled states tomography [48]. Furthermore, an alternative approach to tackling the high dimensionality of quantum state tomography is to employ machine learning-based methods, such as neural networks [4], [49]–[58] and generative models [12], [59]. Building upon the aforementioned methodologies and measurement protocols such as classical shadows [7], QST and its variants have inspired numerous related studies and applications, such as theoretical analyses of QST [60]–[62], quantum measurement tomography [63]–[68] and quantum process tomography [69]–[72].

We discuss a few related works that our work bears the most resemblance to. Our work considers performing QST procedures based on the classical shadow approximation, and our work uses a sketching algorithm to reconstruct the true state. In [8], [45], the authors consider a sketching-based procedure for MPS, but they do not use classical shadow, and the procedure does not have a convergence guarantee. In [46], [73], [74], the authors consider a variational approach to train a matrix product operator (MPO) ansatz based on minimizing measurement error, where the proposed algorithm does not have a formal convergence guarantee, and the system size in the experiments remains relatively small. In [75], the authors consider directly fitting an MPO ansatz to fit the classical shadow approximation \(\hat{\rho}\), but the procedure uses classical shadow protocol under Haar-random projective measurements, whereas our work focuses on the Pauli measurement setting.

Tensor trains (matrix product states) in quantum physics The tensor train (TT) or matrix product state (MPS) ansatz, which originates from the density matrix renormalization group (DMRG) algorithm [5], provides an efficient parametrization of entangled quantum states in high-dimensional Hilbert spaces. For a complete discussion of TT/MPS and its variants, we refer the readers to the following review articles [76]–[79]. As one of the most powerful tools in quantum many-body physics, the TT/MPS ansatz and its variants have been applied to a wide range of related tasks. Examples of such tasks include but are not limited to computing ground states via quantum Monte Carlo methods [80]–[83], studies of electronic structure theory and quantum chemistry [84], [85], simulating open quantum systems and quantum dynamics in general [86]–[113], solving quantum impurity models [114]–[133], quantum simulation and quantum computing [134]–[140], studies of quantum field theories [141]–[146], representing and learning Feynman diagrams [147], [148], etc.

6 Notations↩︎

In addition to the notations of tensors and tensor networks covered above, we provide a brief review of other mathematical notations in this section. We use \(\otimes\) to denote the tensor product. For measuring distances between probability distributions, we use \(D_{\mathrm{KL}}(\cdot,\cdot)\) to denote the Kullback-Leibler (KL) divergence between any two probability distributions. For two probability distributions \(P, Q\), we use \(P \ll Q\) to mean that \(P\) is absolutely continuous with respect to \(Q\). With fixed \(d \in {\mathbb{N}}\), we use \({\boldsymbol{I}}_d\) to denote the identity matrix of size \(d \times d\). For any complex vector \({\boldsymbol{a}}\in \mathbb{C}^d\), we use \({\boldsymbol{a}}^\intercal\) and \({\boldsymbol{a}}^\ast\) to denote the transpose and conjugate transpose of \({\boldsymbol{a}}\), respectively. Moreover, the dot product between any two square matrices \(A,B \in \mathbb{C}^{d \times d}\) is denoted by \(\langle A,B \rangle_F = \mathrm{tr}\left[A^* B\right]\), where \(A^\ast\) denotes the conjugate transpose of \(A\). Regarding the notation of norms, we use \(\left\|\cdot\right\|_1\) and \(\left\|\cdot\right\|_{F}\) to denote the \(l_1\) norm and the Frobenius norm, respectively. In particular, for any vector \(v \in \mathbb{C}^n\) and matrix \(M \in \mathbb{C}^{m \times n}\), we use \(\|v\|\) and \(\|M\| = \sup_{\|u\|=1}\|Mu\|\) to denote the \(l_2\) norm of \(v\) and the induced operator \(2\)-norm of \(M\), respectively. The singular values of \(M\) are denoted by \(s_1(M) \geq s_2(M) \geq \cdots \geq s_l(M)\), where \(l = \min\{m,n\}\). For any \(p \in (0,1)\), the Bernoulli distribution with mean \(p\) is denoted by \(\text{Ber}(p)\). For any two given quantities \(f\) and \(g\), we write \(f \gtrsim g\) when the inequality \(f \geq Cg\) holds for some fixed constant \(C > 0\). Finally, we use the standard symbols \(X,Y,Z\) to denote the Pauli matrices in \(\mathbb{C}^{2 \times 2}\), which satisfy \[X = \begin{bmatrix} 0 &1\\ 1 &0 \end{bmatrix}, \; Y = iXZ = \begin{bmatrix} 0 &-i\\ i &0 \end{bmatrix}, \; Z = \begin{bmatrix} 1 &0\\ 0 &-1 \end{bmatrix}.\]

7 Details for the sketch tomography procedure↩︎

We go through the detailed procedure and derivation for sketch tomography. We assume basic familiarity with the tensor network ansatz and tensor diagrams.

7.1 TT representation of \(\rho\)↩︎

We first verify the claim in the main text that \(\rho\) has a TT format when the target state \(\ket{\psi}\) is a matrix product state. In terms of a tensor diagram, one can represent \(\ket{\psi}\) as follows: \[\label{eqn:32State32representation} \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor,minimum width=4cm] (psi) {\Ket{\psi}}; \foreach \x/\lab in {0.5/{j_1}, 1.5/{j_2}, 3.5/{j_n}}{ \draw[leg] ((psi.south west)+(\x,0)) -- ++(0,-\TNVerticalLeg); \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } \foreach \x/\lab in {2.5/{\cdots}}{ \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } } } \; = \; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[right=\TNHorizontalLeg of A2] (dots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of dots] (AN) {F_n}; \draw[leg] (A1.east)--(A2.west); \draw[leg] (A2.east)--(dots.west); \draw[leg] (dots.east)--(AN.west); \foreach \T/\lab in {A1/{j_1}, A2/{j_2}, AN/{j_n}}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); \node[below] at ((\T.south)+(0,-\TNVerticalLeg)) {\lab}; } } }\,.\tag{9}\]

As is shown in 9 , there exists a collection of \(n\) tensor components \((F_k)_{k = 1}^{n}\), where \(F_1 \in {\mathbb{C}}^{2 \times a_1}\), \(F_k \in {\mathbb{C}}^{a_{k - 1} \times 2 \times a_k}\) for \(k = 2, \ldots, n - 1\), and \(F_n \in {\mathbb{C}}^{a_{n-1} \times 2}\). The evaluation of any entry in \(\ket{\psi}\) is shown in 9 , and one can write it equivalently as follows: \[\label{eqn:32State32representation32ver322} \ket{\psi}(j_1, \ldots, j_n) = \sum_{\gamma_{1}, \ldots, \gamma_{n-1}}F_1(j_1, \gamma_1)F_2(\gamma_1, j_2, \gamma_2)\cdots F_n(\gamma_{n-1}, j_n).\tag{10}\]

The associated density matrix \(\rho = \ket{\psi}\bra{\psi}\) is the target object that sketch tomography approximates. One can first use 9 to directly write \(\rho\) in a tensor diagram as follows: \[\label{eqn:32density32matrix32representation} \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor,minimum width=4cm] (psi) {\rho}; \foreach \x/\lab in {0.5/{j_1}, 1.5/{j_2}, 3.5/{j_n}}{ \draw[leg] ((psi.south west)+(\x,0)) -- ++(0,-\TNVerticalLeg); \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; \draw[leg] ((psi.north west)+(\x,0)) -- ++(0,\TNVerticalLeg); \node[above] at ((psi.north west)+(\x,\TNVerticalLeg)) {\lab'}; } \foreach \x/\lab in {2.5/{\cdots}}{ \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; \node[above] at ((psi.north west)+(\x,\TNVerticalLeg)) {\lab}; } } } \; = \; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\Bar{F_1}}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[tensor, above = \TNVerticalLeg of A2] (B2) {\Bar{F_2}}; \node[right=\TNHorizontalLeg of A2] (Adots) {\ldots}; \node[right=\TNHorizontalLeg of B2] (Bdots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of Adots] (AN) {F_n}; \node[tensor, above = \TNVerticalLeg of AN] (BN) {\Bar{F_n}}; \foreach \a/\b in {A1/A2, A2/Adots, Adots/AN, B1/B2, B2/Bdots, Bdots/BN}{ \draw[leg] (\a.east)--(\b.west); } \foreach \T/\lab in {A1/{j_1}, A2/{j_2}, AN/{j_n}}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); \node[below] at ((\T.south)+(0,-\TNVerticalLeg)) {\lab}; } \foreach \T/\lab in {B1/{j_1'}, B2/{j_2'}, BN/{j_n'}}{ \draw[leg] (\T.north)--++(0,\TNVerticalLeg); \node[above] at ((\T.north)+(0,\TNVerticalLeg)) {\lab}; } } }\,.\tag{11}\]

We derive the TT format of \(\rho\) by writing 11 with the Pauli matrices as the basis. For \(k \in [n]\), we let \((\sigma^{X}_k, \sigma^{Y}_k, \sigma^{Z}_k)\) denote the Pauli matrices on site \(k\), and we write \((\sigma^{1}_{k}, \sigma^{2}_{k}, \sigma^{3}_{k}, \sigma^{4}_{k}) = (\frac{1}{\sqrt{2}}I_{2}, \frac{1}{\sqrt{2}}\sigma^{X}_{k}, \frac{1}{\sqrt{2}}\sigma^{Y}_{k}, \frac{1}{\sqrt{2}}\sigma^{Z}_{k})\). The Pauli matrices \((\sigma^{i}_{k})_{i \in [4], k \in [n]}\) allow one to write 11 in a TT format. We construct tensor components \((G_k)_{k = 1}^{n}\) directly from the MPS ansatz of \(\ket{\psi}\) by the following tensor diagrams: \[\label{eqn:32def32of32Gk32ver322} \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {G_k}; \node[left=\TNHorizontalLeg of A1] (A0) {}; \node[right=\TNHorizontalLeg of A1] (A2) {}; \foreach \x/\lab in {-0.4/{i_1}, -0.6/{i_{k-1}}}{ \draw[leg] ((A1.north east)+(0,\x)) -- ++(\TNHorizontalLeg,0); \draw[leg] ((A1.north west)+(0,\x)) -- ++(-\TNHorizontalLeg,0); } \draw[leg] (A1.north) -- ++(0,\TNVerticalLeg); \node[above] at ((A1.north)+(0,\TNVerticalLeg)) {i}; \node[below] at ((A1.south)+(0,-\TNVerticalLeg)) {}; } } \;=\; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_k}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\sigma_{k}^{i}}; \node[tensor, above = \TNVerticalLeg of B1] (C1) {\Bar{F_k}}; \draw[leg] (A1.north)--(B1.south); \draw[leg] (B1.north)--(C1.south); \node[left=\TNHorizontalLeg of A1] (A0) {}; \node[left=\TNVerticalLeg of C1] (C0) {}; \node[right=\TNHorizontalLeg of A1] (A2) {}; \node[right=\TNVerticalLeg of C1] (C2) {}; \foreach \a/\b in {A0/A1, C0/C1, A1/A2, C1/C2}{ \draw[leg] (\a.east)--(\b.west); } } }, \, \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {G_1}; \node[right=\TNHorizontalLeg of A1] (A2) {}; \foreach \x/\lab in {-0.4/{i_1}, -0.6/{i_{k-1}}}{ \draw[leg] ((A1.north east)+(0,\x)) -- ++(\TNHorizontalLeg,0); } \draw[leg] (A1.north) -- ++(0,\TNVerticalLeg); \node[above] at ((A1.north)+(0,\TNVerticalLeg)) {i}; \node[below] at ((A1.south)+(0,-\TNVerticalLeg)) {}; } } \;=\; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\sigma_{1}^{i}}; \node[tensor, above = \TNVerticalLeg of B1] (C1) {\Bar{F_1}}; \draw[leg] (A1.north)--(B1.south); \draw[leg] (B1.north)--(C1.south); \node[right=\TNHorizontalLeg of A1] (A2) {}; \node[right=\TNVerticalLeg of C1] (C2) {}; \foreach \a/\b in {A1/A2, C1/C2}{ \draw[leg] (\a.east)--(\b.west); } } }, \, \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {G_n}; \node[left=\TNHorizontalLeg of A1] (A0) {}; \foreach \x/\lab in {-0.4/{i_1}, -0.6/{i_{k-1}}}{ \draw[leg] ((A1.north west)+(0,\x)) -- ++(-\TNHorizontalLeg,0); } \draw[leg] (A1.north) -- ++(0,\TNVerticalLeg); \node[above] at ((A1.north)+(0,\TNVerticalLeg)) {i}; \node[below] at ((A1.south)+(0,-\TNVerticalLeg)) {}; } } \;=\; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_n}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\sigma_{n}^{i}}; \node[tensor, above = \TNVerticalLeg of B1] (C1) {\Bar{F_n}}; \draw[leg] (A1.north)--(B1.south); \draw[leg] (B1.north)--(C1.south); \node[left=\TNHorizontalLeg of A1] (A0) {}; \node[left=\TNVerticalLeg of C1] (C0) {}; \foreach \a/\b in {A0/A1, C0/C1}{ \draw[leg] (\a.east)--(\b.west); } } }\,.\tag{12}\] In other words, for \(k \not \in \{1, n\}\), we construct \(G_k \in \mathbb{R}^{a_{k-1}^2 \times 4 \times a_{k}^2}\) by \[\label{eqn:32construction32of32G95k} G_{k}((\gamma_{k-1}, \gamma_{k-1}'), i, (\gamma_{k}, \gamma_{k}')) = \sum_{j, j'}\sigma^{i}_{k}(j, j')F_{k}(\gamma_{k-1}, j, \gamma_{k})\Bar{F_{k}}(\gamma'_{k-1}, j', \gamma'_{k}),\tag{13}\] and the construction for \(G_1\) and \(G_n\) can be likewise derived from 13 by respectively omitting the \((\gamma_{k-1},\gamma'_{k-1})\) variable and the \((\gamma_{k},\gamma'_{k})\) variable. Lastly, 13 shows that the entries of \(G_k\) are defined by an inner product between two Hermitian matrices, and so all entries of \(G_k\) are real numbers.

We now verify that \(\rho\) admits a TT format given by \((G_k)_{k=1}^{n}\). Because multi-linear products of Pauli matrices form a basis in the space of Hermitian matrices in \({\mathbb{C}}^{2^n \times 2^n}\), there exists a tensor \(C \colon [4]^n \to \mathbb{R}\) whereby the term \(\rho\) is given by \[\rho = \sum_{i_1, \ldots, i_n =1}^{4}C(i_1, \ldots, i_n)\prod_{l = 1}^{n}\sigma_{l}^{i_l}.\] Moreover, the orthonormality of the Pauli matrices implies the following equation \[C(i_1, \ldots, i_n) = \sum_{j_1, j_1', \ldots, j_n, j_n'}\rho((j_1, j_1'), \ldots, (j_n, j_n'))\prod_{l = 1}^{n}\sigma_{l}^{i_l}(j_l, j_l').\] Therefore, one can plug in the definition of \(\rho\) in 11 to directly derive \(C\). One has \[\vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor,minimum width=4cm] (psi) {C}; \foreach \x/\lab in {0.5/{i_1}, 1.5/{i_2}, 3.5/{i_n}}{ \draw[leg] ((psi.south west)+(\x,0)) -- ++(0,-\TNVerticalLeg); \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } \foreach \x/\lab in {2.5/{\cdots}}{ \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; \node[above] at ((psi.north west)+(\x,2*\TNVerticalLeg)) {}; } } } \; = \; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\sigma_{1}^{i_1}}; \node[tensor, above = \TNVerticalLeg of B1] (C1) {\Bar{F_1}}; \node[tensor, right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[tensor, above = \TNVerticalLeg of A2] (B2) {\sigma_{2}^{i_2}}; \node[tensor, above = \TNVerticalLeg of B2] (C2) {\Bar{F_2}}; \node[right=\TNHorizontalLeg of A2] (Adots) {\ldots}; \node[right=\TNHorizontalLeg of C2] (Cdots) {\ldots}; \node[tensor, right=\TNHorizontalLeg of Adots] (AN) {F_n}; \node[tensor, above = \TNVerticalLeg of AN] (BN) {\sigma_{n}^{i_n}}; \node[tensor, above = \TNVerticalLeg of BN] (CN) {\Bar{F_n}}; \foreach \a/\b in {A1/A2, A2/Adots, Adots/AN, C1/C2, C2/Cdots, Cdots/CN}{ \draw[leg] (\a.east)--(\b.west); } \foreach \T in {A1,A2,AN}{ \draw[leg] (\T.north)--++(0,\TNVerticalLeg); } \foreach \T in {C1,C2,CN}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); } } } = \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {G_1}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {G_2}; \node[right=\TNHorizontalLeg of A2] (dots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of dots] (AN) {G_n}; \draw[leg] (A1.east)--(A2.west); \draw[leg] (A2.east)--(dots.west); \draw[leg] (dots.east)--(AN.west); \foreach \T/\lab in {A1/{i_1}, A2/{i_2}, AN/{i_n}}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); \node[below] at ((\T.south)+(0,-\TNVerticalLeg)) {\lab}; \node[above] at ((\T.north)+(0,2*\TNVerticalLeg)) {}; } } }\, ,\] where the second equality is from 12 .

7.2 Sketch tomography algorithm for quantum state tomography↩︎

We cover how to approximate \(\rho\) with the classical shadow protocol. The main idea for the procedure is covered in the main text, and here we go through the full implementation details.

We explain how one applies the classical shadow protocol from Pauli measurements. Let \(B\) be the number of recorded measurements. The classical shadow protocol uses the measurement outcome to produce a collection of \(2 \times 2\) matrices \(\{M_{l}^{(j)} \in {\mathbb{C}}^{2 \times 2}\}_{l \in [n], j \in [B]}\). The resultant classical shadow approximation \(\hat{\rho}\) is written as a sum of separable matrices as follows: \[\label{eqn:32formula32of32rho32hat} \hat{\rho} = \frac{1}{B}\sum_{j = 1}^{B} \bigotimes_{l =1}^{n} M_{l}^{(j)}.\tag{14}\] 14 allows for an efficient calculation of the partial trace of \(\hat{\rho}\). Therefore, performing observable approximation over 14 is efficient for any local observable \(O\). When \(O\) is \(k\)-local, the estimator \(\mathrm{tr}\left[O\hat{\rho}\right]\) is unbiased with a variance bounded from above by \({\mathcal{O}}(4^k)\), which is moderate for small \(k\).

In practice, instead of a single classical shadow approximation \(\hat{\rho}\) obtained from \(B\) experiments on \(\rho\), one can repeat the same protocol \(W\) times to form classical shadow approximators \(\hat{\rho}_{1}, \ldots, \hat{\rho}_{W}\). Then, the observable approximation is done by \[\langle O \rangle = \mathrm{tr}\left[O\rho\right] \approx \mathrm{median}\left\{\mathrm{tr}\left[O\hat{\rho}_1\right], \ldots , \mathrm{tr}\left[O\hat{\rho}_W\right]\right\},\] where \(\mathrm{median}\{a_1, \ldots, a_W\}\) is the median of the \(W\) scalars \(a_1, \ldots, a_W\). The median-of-mean estimator is unbiased and is more robust with respect to outliers.

We go through our QST procedure for obtaining \(\rho\). For notational compactness, for an index set \(S = \{s_1, \ldots, s_l\} \subset [d]\), we use a multi-index notation by letting \(i_{S} := (i_{s_1}, \ldots, i_{s_l})\), i.e. \(i_{S}\) stands for the subvector of variables with entries from the index set \(S\). Using the multi-index notation, any density matrix \(\rho\) can be uniquely represented by a tensor \(C \colon [4]^n \to \mathbb{R}\) as follows: \[\rho = \sum_{i_{[n]} \in [4]^n}C(i_{[n]})\prod_{l = 1}^{n}\sigma_{l}^{i_l}.\]

Given the MPS assumption on \(\ket{\psi}\), the tensor \(C\) can be written in terms of tensor components \((G_{k})_{k = 1}^{n}\), where \(G_{1} \in \mathbb{R}^{ 4 \times r_{1}}\), \(G_{k} \in \mathbb{R}^{r_{k-1} \times 4 \times r_{k}}\) for \(k = 2,\ldots, n-1\), and \(G_{n} \in \mathbb{R}^{r_{n-1} \times 4 }\). The equation for \(C\) is written as follows: \[\label{eq:32TT32supplemental} C(i_{[n]}) =\sum_{\alpha_{[n-1]}}G_{1}(i_{1}, \alpha_{1})G_{2}(\alpha_{1}, i_{2}, \alpha_{2})\cdots G_{n}(\alpha_{n-1}, i_{n}).\tag{15}\]

To obtain \(\rho\), we obtain each \(G_{k}\) for \(k = 1,\ldots,n\). We first cover the case where \(k \not \in \{1, n\}\). Due to the structural equation in 15 , there exist tensors \(C_{<k} \colon [4]^{k-1} \times [r_{k-1}] \to {\mathbb{R}}\) and \(C_{>k} \colon [r_{k}] \times [4]^{n-k} \to {\mathbb{R}}\) so that the following equation holds: \[\label{eq:32full32equation32for32Gk} \sum_{\alpha_{k-1}, \alpha_{k}}C_{<k}(i_{[k-1]}, \alpha_{k-1})G_{k}(\alpha_{k-1}, i_k, \alpha_{k})C_{>k}(\alpha_{k}, i_{[n] - [k]}) = C(i_{[n]}).\tag{16}\]

One then applies sketching to 16 to form a tractable linear system for \(G_k\). We let \(\tilde{r}_{k-1}, \tilde{r}_{k}\) be two integers such that \(\tilde{r}_{k-1} \geq r_{k-1}\) and \(\tilde{r}_{k} \geq r_{k}\). We introduce two sketch tensors \(S_{<k} \colon [\tilde{r}_{k-1}] \times [4]^{k-1} \to \mathbb{R}, S_{>k} \colon [4]^{n-k} \times [\tilde{r}_{k}] \to \mathbb{R}\). By contracting 16 with \(S_{>k}\) and \(S_{<k}\), one obtains \[\label{eqn:32sketched32linear32system32for32Gk} \sum_{\alpha_{k-1}, \alpha_{k}}A_{<k}(\zeta, \alpha_{k-1})G_{k}(\alpha_{k-1}, i_{k}, \alpha_{k})A_{>k}(\alpha_{k}, \mu) = B_{k}(\zeta, i_k, \mu),\tag{17}\] where \(A_{>k}, A_{<k}\) are respectively the contraction of \(C_{>k}, C_{<k}\) by \(S_{>k}, S_{<k}\) over \(i_{[n] - [k]}, i_{[k-1]}\) and \(B_{k}\) is the contraction of \(C\) by \(S_{>k} \otimes S_{<k}\) over \(i_{[n] - \left\{ k \right\}}\).

In practice, the contractions are done implicitly through taking observables. We let \(\{L_{k}^{\zeta}\}_{\zeta \in [\tilde{r}_{k-1}]}\) and \(\{R_{k}^{\mu}\}_{\mu \in [\tilde{r}_{k}]}\) be observables associated with \(S_{<k}, S_{>k}\) through the following equation: \[\label{eqn:32implicit32construction32of32sketch32tensor} L_{k}^{\zeta} = \sum_{i_{[k-1]}}S_{<k}(\zeta, i_{[k-1]})\prod_{l =1}^{k-1}\sigma_{l}^{i_l}, \quad R_{k}^{\mu} = \sum_{i_{[n] - [k]}}S_{>k}(i_{[n] - [k]}, \mu)\prod_{l = k+1}^{n}\sigma_{l}^{i_l}.\tag{18}\] Then, to calculate \(B_k\), we utilize the fact that \(\{\prod_{l = 1}^{n}\sigma_{l}^{i_l}\}_{i_{[n]} \in [4]^n}\) is an orthonormal basis in the space of Hermitian matrices in \({\mathbb{C}}^{2^n \times 2^n}\), and so the following holds for \(B_k\): \[B_{k}(\zeta, i_k, \mu) = \sum_{i_{[n] - \left\{ k \right\}}}C(i_{[n]})S_{<k}(\zeta, i_{[k-1]})S_{>k}(i_{[n] - [k]}, \mu)=\mathrm{tr}\left[\rho L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right],\] where the second equality holds by the definition of \(\{L_{k}^{\zeta}\}_{\zeta \in [\tilde{r}_{k-1}]}\) and \(\{R_{k}^{\mu}\}_{\mu \in [\tilde{r}_{k}]}\). With the classical shadow protocol, we approximate \(B_k\) by \[\label{eqn:32definition32of32B95k} \begin{align} B_{k}(\zeta, i_k, \mu) \approx \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right]\right\}. \end{align}\tag{19}\]

There exists a gauge degree of freedom for the TT representation of \(C\) at each internal bond of the tensor. More explicitly, the pair of tensors \((C_{>k}, C_{<(k+1)})\) determines the value of \((A_{>k}, A_{<(k+1)})\), and so obtaining \(A_{>k}\) is dependent on fixing the gauge on \((C_{>k}, C_{<(k+1)})\). To fix the gauge and obtain \(A_{<k}, A_{>k}\), we assume that all sketch tensors \(\{S_{>1}\} \cup \{S_{>j}, S_{<j}\}_{j=2}^{n-1} \cup \{S_{<n}\}\) have been defined via 18 . We use sketching of \(C\) to uniquely determine a gauge on \((C_{>k}, C_{<(k+1)})\).

To obtain \(A_{>k}\), we construct a tensor \(Z_{k}\) from contracting \(C\) with \(S_{<(k+1)} \otimes S_{>k}\). By associating \(S_{<(k+1)}\) with observables \(\{L_{k+1}^{\zeta'}\}_{\zeta' \in [\tilde{r}_{k}]}\), we have \[\label{eqn:32Z95k32formula} Z_{k}(\zeta', \mu) = \mathrm{tr}\left[\rho L_{k+1}^{\zeta'}R_{k}^{\mu}\right] \approx \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{k+1}^{\zeta'}R_{k}^{\mu}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{k+1}^{\zeta'}R_{k}^{\mu}\right]\right\},\tag{20}\] where the equality holds by the definition of \(\{L_{k+1}^{\zeta'}\}_{\zeta' \in [\tilde{r}_{k}]}\) and \(\{R_{k}^{\mu}\}_{\mu \in [\tilde{r}_{k}]}\). The SVD of \(Z_{k}\) uniquely determines the best rank-\(r_k\) factorization \(Z_{k} = USV^{\top}\). For example, the choice of setting \(A_{<(k+1)} = U\) and \(A_{>k} = SV^{\top}\) uniquely fixes the gauge of \((C_{>k}, C_{<(k+1)})\). One way to understand the construction is that \(C_{>k}\) is chosen to be the unique tensor for which the contraction of \(C_{>k}\) with \(S_{>k}\) is \(U\). Finally, obtaining \(A_{<k}\) is done through the best rank-\(r_{k-1}\) factorization of \(Z_{k-1}\).

We now detail the procedure for \(k = 1\) and \(k = n\). For \(G_1\) and \(G_n\), one can likewise use sketch tensors \(S_{>1}\) and \(S_{<n}\) to form the following linear system \[\label{eqn:32sketched32linear32system32boundary32case} \begin{align} &\sum_{\alpha_{1}}G_{1}(i_{1}, \alpha_{1})A_{>1}(\alpha_{1}, \mu) = B_{1}(i_1, \mu),\\ &\sum_{\alpha_{n-1}}A_{<n}(\zeta, \alpha_{n-1})G_{n}(\alpha_{n-1}, i_{n}) = B_{n}(\zeta, i_n). \end{align}\tag{21}\]

The term \(A_{>1}\) and \(A_{<n}\) are respectively obtained using \(Z_{1}\) and \(Z_{n-1}\) obtained via 20 . Lastly, the sketch tensors \(S_{>1}, S_{<n}\) are respectively defined implicitly with observables \(\{R_{1}^{\mu}\}_{\mu \in [\tilde{r}_{1}]}\) and \(\{L_{n}^{\zeta}\}_{\zeta \in [\tilde{r}_{n-1}]}\). One obtains \(B_1\) and \(B_n\) by \[\label{eqn:32definition32of32B95k32boundary32case} \begin{align} &B_{1}(i_1, \mu) \approx \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 \sigma_{1}^{i_1}R_{1}^{\mu}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W \sigma_{1}^{i_1}R_{1}^{\mu}\right]\right\}, \\ &B_{n}(\zeta, i_n) \approx \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{n}^{\zeta}\sigma_{n}^{i_n}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{n}^{\zeta}\sigma_{n}^{i_n}\right]\right\}. \end{align}\tag{22}\]

We summarize our approach in 4. We use the classical shadow to obtain each \(Z_{k}, B_{k}\) and we use SVD on each \(Z_{k}\) to obtain all \(A_{>k}, A_{<k}\). Then, solving the corresponding sketched linear equation for all \(G_k\) allows one to obtain a TT approximation of \(\rho\).

8 Proof of the upper bound↩︎

This section goes through the convergence guarantee for the sketch tomography procedure. The ground truth is a density matrix \(\rho\) whose coefficient tensor is a tensor train. Sketch tomography takes in the classical shadow approximation \(\hat{\rho}\) and outputs an approximate density matrix \(\tilde{\rho}\) with tensor component \((\hat{G}_{k})_{k = 1}^{n}\). The organization of this section is as follows. In 8.1, we state the necessary assumptions for the results of this section. In 8.2, we show that each observable estimation task required by 4 is accurate. In 8.3, we prove that each tensor component of \(\tilde{\rho}\) is close to the tensor component of \(\rho\). In 8.4, we prove that \(\tilde{\rho}\) is close to \(\rho\) in the Frobenius norm.

Notations We summarize the necessary notations for the analysis. For a vector \(v\), the term \(\lVert v \rVert\) is its \(l_2\) norm. For a matrix \(M\), the term \(\lVert M \rVert\) is its operator norm, and the term \(s_j(M)\) is the \(j\)-th singular value of \(M\). For a tensor \(G\), the term \(\lVert G\rVert_{F}\) is the Frobenius norm of \(G\). We also assume the multi-index notation used in Appendix 7.

8.1 Assumptions↩︎

We list out the assumptions for the convergence guarantee of sketch tomography. As in Appendix 7, we use \(\rho\) to represent the target density matrix on \(n\) qubits. For easy reference, we summarize the construction of the Pauli-basis representation of \(\rho\):

Following 1, we formalize the tensor train (TT) property on \(\rho\) as an assumption. In particular, Appendix 7 shows that the following TT assumption on \(\rho\) holds when \(\rho\) is the density matrix of a matrix product state \(\ket{\psi}\).

Secondly, we place assumptions to constrain the design space for sketch tensors. As shown in Appendix 7, the sketch tensors are defined implicitly through observables. The integers \(\{\tilde{r}_k\}_{k = 1}^{n-1}\) represent the number of sketch functions used during sketch tomography. The sketching procedure uses sketch tensors \(\{S_{<k}, S_{>k}\}_{k = 1}^{n}\) which are defined implicitly through observables. We summarize the construction here:

The second assumption requires that the associated observables of the sketch tensors should come from local observables. Essentially, this assumption is to ensure that the classical shadow protocol under random Pauli measurement can reliably perform observable estimations required by 4.

2 is a mild assumption on the design of sketch tensors. One way to satisfy 2 is to set \(\{L_{k}^{\zeta}\}_{\zeta}\) and \(\{R_{k}^{\mu}\}_{\mu}\) to be a linear combination of 1-local and 2-local Pauli observables with bounded total weights. When running 4 with sketch tensors satisfying 2, one can ensure that all observable estimation done by \(\hat{\rho}\) can achieve a high accuracy with a sufficient number of state copies.

The third assumption chooses how one obtains \(A_{<(k+1)}\) and \(A_{>k}\) from \(Z_{k}\). While this assumption could be precluded, its inclusion significantly simplifies the error analysis. We state it here:

This particular gauge choice is made only for the sake of analysis. In practice, one may instead take \((A_{<(k+1)}, A_{>k})\) from the SVD of \(Z_k\); this corresponds to an orthogonal change of basis on the internal legs and does not affect the asymptotic scaling of our bounds. However, using SVD complicates the notation for error analysis.

8.2 Error bound on observable estimation↩︎

1 and 2 allow us to provide an error bound on all of the observable estimation tasks that are performed when running 4. To obtain \(G_{k}\), we have seen in 2 that one needs to form three sketches \(B_{k} \in \mathbb{R}^{\tilde{r}_{k-1} \times 4 \times \tilde{r}_k}, Z_{k-1} \in \mathbb{R}^{\tilde{r}_{k-1} \times \tilde{r}_{k-1}}, Z_{k} \in \mathbb{R}^{\tilde{r}_{k} \times \tilde{r}_{k}}\) by performing observable estimation on \(\rho\). This section bounds the error rate on the sketch estimation.

For the error estimation, we assume that the reader is familiar with the definition of the shadow norm \(\lVert \cdot \rVert_{\mathrm{shadow}}\) in [7]. For readers unfamiliar with the concept, one can still follow the derivation by assuming the results on \(\lVert \cdot \rVert_{\mathrm{shadow}}\) from [7]. For a brief explanation, the shadow norm essentially measures the complexity of the observable in a given classical shadow protocol. For random Pauli measurement, estimating \(\mathrm{tr}\left[\rho O\right]\) is only accurate when \(\lVert O \rVert_{\mathrm{shadow}}\) is small. Moreover, \(\lVert O \rVert_{\mathrm{shadow}}\) can be reliably bounded from above when \(O\) is a local observable or a sum of local observables. Therefore, 2 essentially forces the sketch tensors to correspond to observables that have a small shadow norm.

From the formula in 2, the sketches \(B_{k}, Z_{k-1}, Z_{k}\) satisfy the following equations: \[\label{eqn:32true32sketch} \begin{align} &B_{k}(\zeta, i_k, \mu) = \mathrm{tr}\left[\rho L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right],\; \\ &Z_{k-1}(\zeta, \mu') = \mathrm{tr}\left[\rho L_{k}^{\zeta}R_{k-1}^{\mu'}\right],\;\\ &Z_{k}(\zeta', \mu) = \mathrm{tr}\left[\rho L_{k+1}^{\zeta'}R_{k}^{\mu}\right]. \end{align}\tag{24}\]

During the procedure of 4, we use \(W\) copies of the classical shadow estimator to form the approximated sketches \(\hat{B}_{k}, \hat{Z}_{k-1}, \hat{Z}_{k}\) with the following formula: \[\label{eqn:32approximated32sketch} \begin{align} &\hat{B}_{k}(\zeta, i_k, \mu) = \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{k}^{\zeta}\sigma_{k}^{i_k}R_{k}^{\mu}\right]\right\}, \\ &\hat{Z}_{k-1}(\zeta, \mu') = \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{k}^{\zeta}R_{k-1}^{\mu'}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{k}^{\zeta}R_{k-1}^{\mu'}\right]\right\}, \\ &\hat{Z}_{k}(\zeta', \mu) = \mathrm{median}\left\{\mathrm{tr}\left[\hat{\rho}_1 L_{k+1}^{\zeta'}R_{k}^{\mu}\right], \ldots, \mathrm{tr}\left[\hat{\rho}_W L_{k+1}^{\zeta'}R_{k}^{\mu}\right] \right\}. \end{align}\tag{25}\]

2 allows us to bound the error in approximating the formed sketches. We show that one has the following result:

For the constant \(c_S\), one typically takes \(D_1 = 2\) and \(D_2 = 1\) to ensure the classical shadow procedure is successful. In practice, instead of enforcing 2, one can also directly design sketch tensors so that they correspond to observables with a small shadow norm.

8.3 Error bound on tensor component \(G_k\)↩︎

This section proves the component-wise convergence result. 4 outputs a series of tensor components \((\hat{G}_{k})_{k = 1}^{n}\) to form an approximation \(\tilde{\rho}\). We prove that the output \(\hat{G}_{k}\) accurately approximates \(G_{k}\).

For the error analysis, we need to invoke a result on the stability of the least-squares solution under perturbations of the linear system.

We now present the proof of 4.

8.4 Global error bound on estimating \(\rho\)↩︎

We have established that the tensor components \((\hat{G}_k)_{k = 1}^{n}\) converge to \((G_k)_{k = 1}^{n}\) at an \({\mathcal{O}}(B^{-1/2})\) rate when given \({\mathcal{O}}(B)\) copies of \(\rho\). This subsection shows how the local component-wise convergence ensures that the final output \(\tilde{\rho}\) converges to \(\rho\). The analysis in this section uses a novel perturbational result for tensor train. This part of the error analysis heavily uses tensor reshaping, and we frequently reshape tensors into matrices. Therefore, we introduce an additional matricization notation that simplifies the analysis.

We present a result on how local component-wise error propagates to the global tensor train error in the Frobenius norm. To the best of our knowledge, this result is a novel contribution first derived in this work. We anticipate that this bound is of independent interest for the tensor network community. Minimal modifications would generalize the result to tree-based tensor networks as defined in [150].

One can see that 2 can not be directly applied to our setting unless the coefficient tensor \(C\) has maximal rank \(r \leq 2\). However, the result in 2 becomes applicable after one merges several tensor components into a bigger tensor component. To use 2, we shall merge connected tensor components. For example, if one has \(L\) tensor components \(G_{k}, \ldots, G_{k+L-1}\), then merging these components leads to \[\vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor,minimum width=3cm] (psi) {F}; \node[right=\TNHorizontalLeg of psi] (empty) {}; \draw[leg] (empty.west)--(psi.east); \foreach \x/\lab in {0.4/{i_{k}}, 2.6/{i_{k+L-1}}}{ \draw[leg] ((psi.south west)+(\x,0)) -- ++(0,-\TNVerticalLeg); \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } \foreach \x/\lab in {1.5/{\cdots}}{ \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } \draw[leg] (psi.west)--++(-\TNHorizontalLeg,0); } } \;=\; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {G_{k}}; \node[left=\TNHorizontalLeg of A1] (empty) {}; \node[right=\TNHorizontalLeg of A1] (dots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of dots] (AN) {G_{k+L-1}}; \draw[leg] (empty.east)--(A1.west); \draw[leg] (A1.east)--(dots.west); \draw[leg] (dots.east)--(AN.west); \foreach \T/\lab in {A1/{i_{k}}, AN/{i_{k+L-1}}}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); \node[below] at ((\T.south)+(0,-\TNVerticalLeg)) {\lab}; } \draw[leg] (AN.east)--++(\TNHorizontalLeg,0); } },\] which one views as a 3-tensor of size \(r_{k-1} \times 4^L \times r_{k+L-1}\). One sees that it suffices to take merge \(L = \lceil \log_2{(r_{\mathrm{max}}}) \rceil\) blocks to ensure that the condition in 2 can hold.

We are now ready to prove our main result.

We give some remarks on the bound that we have obtained in 5. The \(c_G^{L}\) dependency is mild because one has \(L = \lceil \log_2(r_{\mathrm{max}}) \rceil\). However, the \(c_G^{L}\) dependency is due to the loose bound in 36 . We conjecture that a more refined analysis can improve the \(c_G^{L}\) dependence to an \({\mathcal{O}}(1)\) dependence that is independent of \(L\). We leave a more refined analysis for future work.

Moreover, the assumption of \(n \geq L^2\) is not necessary. For the regime of \(L \leq n < L^2\), one is not guaranteed to write an integer less than \(L^2 - L - 1\) by a weighted sum of \(L\) and \(L+1\). Therefore, one would require a more complicated construction for the merging operation to construct \((F_j)_{j = 1}^{d}\). That said, for \(n \leq L^2\), one can obtain a similar result as obtained in 5, but the worst case would involve an \({\mathcal{O}}(c_G^{4L})\) dependence.

9 Detail on maximum likelihood estimation training using MPS↩︎

9.1 Methodology↩︎

We follow the procedure to utilize the Pauli measurement data for MLE training. One has access to \(B\) unitary transformations \(\{U^{(j)} \in U(2^n)\}_{j = 1}^{B}\) and the corresponding binary measurement outcomes \(b^{(j)} = (b^{(j)}_1, \ldots, b^{(j)}_n) \in \{0, 1\}^{n}\). The \(j\)-th measurement outcome is also a state \(\ket{b^{(j)}} \in {\mathbb{C}}^{2^n}\) by the computational basis encoding. In the MLE framework, one minimizes a negative-log-likelihood (NLL) loss defined as follows \[\label{eqn:32NLL} \mathcal{L}(\ket{\phi}) = -\frac{1}{B}\sum_{j = 1}^{B}{\log\left(\lvert\bra{b^{(j)}} U^{(j)}\ket{\phi}\rvert^2\right)} + \log(\braket{\phi | \phi}),\tag{38}\] and \(\ket{\phi}\) is optimized over an MPS class. We illustrate the key concepts in a tensor diagram. The ansatz \(\ket{\phi}\) is parameterized by a collection of tensor component \((F_k)_{k = 1}^{n}\) and the representation of \(\ket{\phi}\) is as follows: \[\vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor,minimum width=4cm] (psi) {\Ket{\phi}}; \foreach \x/\lab in {0.5/{j_1}, 1.5/{j_2}, 3.5/{j_n}}{ \draw[leg] ((psi.south west)+(\x,0)) -- ++(0,-\TNVerticalLeg); \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } \foreach \x/\lab in {2.5/{\cdots}}{ \node[below] at ((psi.south west)+(\x,-\TNVerticalLeg)) {\lab}; } } } \; = \; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[right=\TNHorizontalLeg of A2] (dots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of dots] (AN) {F_n}; \draw[leg] (A1.east)--(A2.west); \draw[leg] (A2.east)--(dots.west); \draw[leg] (dots.east)--(AN.west); \foreach \T/\lab in {A1/{j_1}, A2/{j_2}, AN/{j_n}}{ \draw[leg] (\T.south)--++(0,-\TNVerticalLeg); \node[below] at ((\T.south)+(0,-\TNVerticalLeg)) {\lab}; } } }\,.\]

The normalization constant \(Z = \braket{\phi | \phi}\) is evaluated in a tensor diagram as follows: \[\label{eqn:32def32of32Z32in32phi} Z \; =\; \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {\Bar{F_1}}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[tensor, above = \TNVerticalLeg of A2] (B2) {\Bar{F_2}}; \node[right=\TNHorizontalLeg of A2] (Adots) {\ldots}; \node[right=\TNHorizontalLeg of B2] (Bdots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of Adots] (AN) {F_n}; \node[tensor, above = \TNVerticalLeg of AN] (BN) {\Bar{F_n}}; \foreach \a/\b in {A1/A2, A2/Adots, Adots/AN, B1/B2, B2/Bdots, Bdots/BN}{ \draw[leg] (\a.east)--(\b.west); } \foreach \T/\lab in {A1/B1, A2/B2, AN/BN}{ \draw[leg] (\T.north)--++(0,\TNVerticalLeg); } } }\, .\tag{39}\]

Similarly, the evaluation of \(\bra{b^{(j)}} U^{(j)}\ket{\phi}\) is done in a tensor diagram. Let \(U^{(j)} = \otimes_{k = 1}^{n}U_{k}^{(j)}\) be the unitary transformation for the \(j\)-th measurement. One can write the diagram for the evaluation as follows: \[\label{eqn:32def32of32eval32in32phi} \bra{b^{(j)}} U^{(j)}\ket{\phi} = \vcenter{ \tikz[baseline=-0.6ex,scale=0.5,transform shape, every node/.style={font=\Large}]{ \node[tensor] (A1) {F_1}; \node[tensor, above = \TNVerticalLeg of A1] (B1) {U_1^{(j)}}; \node[tensor,right=\TNHorizontalLeg of A1] (A2) {F_2}; \node[tensor, above = \TNVerticalLeg of A2] (B2) {U_2^{(j)}}; \node[right=\TNHorizontalLeg of A2] (Adots) {\ldots}; \node[right=\TNHorizontalLeg of B2] (Bdots) {\ldots}; \node[tensor,right=\TNHorizontalLeg of Adots] (AN) {F_n}; \node[tensor, above = \TNVerticalLeg of AN] (BN) {U_n^{(j)}}; \foreach \a/\b in {A1/A2, A2/Adots, Adots/AN, B1/B2, B2/Bdots, Bdots/BN}{ \draw[leg] (\a.east)--(\b.west); } \foreach \T/\lab in {A1/B1, A2/B2, AN/BN}{ \draw[leg] (\T.north)--++(0,\TNVerticalLeg); } \foreach \T/\lab in {B1/{b_1^{(j)}}, B2/{b_2^{(j)}}, BN/{b_n^{(j)}}}{ \draw[leg] (\T.north)--++(0,\TNVerticalLeg); \node[above] at ((\T.north)+(0,\TNVerticalLeg)) {\lab}; } } }\,.\tag{40}\]

The diagrams in 39 and 40 show that \(\mathcal{L}\) in 38 is efficient to evaluate and differentiate against in an MPS class. For the optimization algorithm, we use a one-site DMRG approach summarized in 5, whereby DMRG refers to the fact that the mixed canonical form has been used and that the parameter update is done sequentially.

9.2 Implementation detail↩︎

We report the implementation details for our MLE benchmark. For the 1D Heisenberg experiment and the 1D TFIM experiment in the main text, we choose a learning rate of \(\eta = 0.1\) and we obtain \(\ket{\phi}\) by running many iterations of 5 until \(\mathcal{L}(\ket{\phi}) \leq \mathcal{L}(\ket{\psi})\). By the definition of the loss function in the MLE procedure, this means that the obtained state \(\ket{\phi}\) is better at explaining the generated data than the true state \(\ket{\psi}\). The ansatz for \(\ket{\phi}\) is taken to be in the same class as \(\ket{\psi}\) in the sense that \(\ket{\phi}\) is optimized over the complex MPS and all tensor components of \(\ket{\phi}\) are of the same shape as \(\ket{\psi}\).

The main text shows that the obtained MLE solution is not accurate in certain observable estimation tasks. One can conclude from the numerical performance that, for practical settings, an MPS model that is competitive in the log likelihood metric can have poor observable estimation accuracy. MLE is trained on the log likelihood metric, and it does not include a penalization term for how well it predicts individual observables. Therefore, when the sample size is not sufficiently large, the fact that \(\ket{\psi}\) is successful in the log likelihood metric is insufficient to guarantee an accurate observable prediction. On the other hand, MLE training is more computationally challenging when the sample size is large. We remark that the conducted experiment is already more ideal than practical settings by using the log likelihood level of the true model \(\ket{\psi}\) as a termination criterion. In practice, one might even prematurely terminate training when in a local minimum, which leads to a worse quality for the MLE solution. Therefore, while MLE is a good practical proposal for quantum state tomography, our proposed sketch tomography procedure is, in some cases, more advantageous due to its performance guarantee in 1 and due to practical performances shown in the main text.

We explain why we do not include a benchmark in the 2D Heisenberg experiment. To calculate the evaluation diagram used in 40 , one needs to at least store a complex tensor of size \(a \times 2 \times a\). Therefore, the memory complexity of calculating \(B\) samples in parallel is on the order of \(4a^2B\) floating points. For \(B = 10^4\) samples, the memory requirement is 6.4 GB, which is already quite substantial. Therefore, training the MLE model has a substantial memory complexity even for a small sample size. Moreover, as mentioned in the main text, the Pauli measurement samples for \(\ket{\psi}\) are too small and will likely lead to overfitting when one chooses \(\ket{\phi}\) to be in the same MPS class as \(\ket{\psi}\). Therefore, to have a fair assessment of MLE, one would also need a much larger sample size, which would drastically increase the computational complexity to perform the training. Therefore, we choose to keep the sample size at \(B = 8 \times 10^5\) and only present the result for sketch tomography. For researchers who wish to train an MLE benchmark in this 2D experiment, we suggest that training the MPS model with a mini-batch stochastic gradient descent strategy might allow one to compare MLE with shadow tomography.

10 Proof of the lower bound↩︎

In this section, we provide both a formal version and a complete proof of the information-theoretic lower bound presented in 2 of 2.2 above. We begin by presenting a mathematical definition of tensor trains as follows, which serves as a complement to the definition based on tensor diagrams presented in (2 ) of 2 above.

With the formal definition of TTs given above, we then present a mathematically rigorous version of the information-theoretic lower bound as below.

Before discussing the proof of 6 above, we need to present a few lemmas. In particular, 3 below discusses properties of the TT class specified in 4 above under tensor products, while 4 provides a universal upper bound on the rank of any TT in \(\mathbb{C}^{2^n \times 2^n}\).

In addition to the results on TTs presented above, we also need a few results from nonparametric statistics in our proof. Specifically, we need Theorem 20 from [152], which is presented as below:

One other result we need here is the Varshamov-Gilbert (VG) Lemma, which is presented as Theorem 2.9 in [153].

We note that the VG lemma is also a key component in the proof of the information-theoretic lower bound established in [154], where it is employed to construct a family of pairwise-separated density matrices. Such a family, which is often referred to as a covering, plays a central role in characterizing fundamental information-theoretic limits of estimation problems from not only statistics [153] but also quantum physics [155], [156]. With all essential tools listed above, we now present a complete proof of the lower bound provided in 6 as follows.

Proof of 6 Without loss of generality, here we assume that \(n\) is some integer divisible by \(8\), as this will only impact the information-theoretic lower bound up to some constant. We begin by taking \(\gamma = \frac{1}{1600B}\), which satisfies \(\epsilon \in \left(0,\frac{1}{1600n}\right)\) from the given assumption \(B \geq n\). Then we consider the following two density matrices \(\sigma,\omega\) in \(\mathbb{C}^{2 \times 2}\) defined as follows: \[\label{eqn:32defn32of32basis32density32matrices} \sigma = \frac{1}{2}{\boldsymbol{I}}_2 -\frac{1}{2}\sqrt{\gamma}X + \frac{1}{2}\sqrt{1-\gamma}Y, \;\omega := \frac{1}{2}{\boldsymbol{I}}_2 + \frac{1}{2}\sqrt{\gamma}X + \frac{1}{2}\sqrt{1-\gamma}Y\tag{41}\] Moreover, applying the VG Lemma listed in 8 above yields a collection \(\mathcal{V}_{n} := \left\{\tau^{(0)},\cdots,\tau^{(2^{\frac{n}{8}})}\right\}\) of \(\left(2^{\frac{n}{8}} + 1\right)\) binary strings in \(\{0,1\}^{n}\) such that \[\label{eqn:32VG32lemma32l132norm32lower32bound} \left\|\tau^{(j)}-\tau^{(k)}\right\|_{1} = \sum_{l=1}^{n}\left|\tau^{(j)}_l-\tau^{(k)}_l\right| \geq \frac{n}{8},\tag{42}\] holds for any \(0 \leq j \neq k \leq 2^{\frac{n}{8}}\). Based on the two density matrices \(\sigma, \omega\) and the collection \(\mathcal{V}_{n}\) of binary strings constructed above, we further define a density matrix \(\rho^{(j)}\) in \({\mathbb{C}}^{n \times n}\) associated with the string \(\tau^{(j)}\) for any \(j \in \left\{0,1,\cdots,2^{\frac{n}{8}}\right\}\), which takes the form of a tensor product as follows: \[\label{eqn:32defn32of32density32matrices32via32tensor32product} \begin{align} \rho^{(j)} &:= \left(\tau^{(j)}_1 \sigma + \left(1-\tau^{(j)}_1\right)\omega\right) \otimes \cdots \otimes \left(\tau^{(j)}_n \sigma + \left(1-\tau^{(j)}_n\right)\omega\right). \end{align}\tag{43}\] Let \(\mathcal{C}_{n} := \left\{\rho^{(0)}, \rho^{(1)}, \cdots,\rho^{(2^{\frac{n}{8}})}\right\}\) be the corresponding collection of density matrices defined in (43 ) above. We begin by first verifying that \({\mathcal{C}}_n \subseteq {\mathcal{M}}_n(r,1)\) for constant \(r \geq 2\). On the one hand, since \(\left\|\sigma\right\|_{F} = \left\|\omega\right\|_{F} = \frac{1}{2}\left(1+\epsilon +1-\epsilon\right) = 1\), a direct computation of the Frobenius norm \(\left\|\rho^{(j)}\right\|_{F}\) yields that \[\label{eqn:32bound32on32frobenius32norm32of32each32rho95j} \begin{align} \left\|\rho^{(j)}\right\|_{F} = \prod_{i=1}^{n}\left\|\tau^{(j)}_i \sigma + \left(1-\tau^{(j)}_i\right)\omega\right\|_{F} = 1, \end{align}\tag{44}\] for any \(0 \leq j \leq 2^{\frac{L}{8}}\). On the other hand, we apply 4 to the \(i\)-th component \(\tau^{(j)}_i \sigma + \left(1-\tau^{(j)}_i\right)\omega\) in the tensor product \(\rho^{(j)}\) to deduce that \[\tau^{(j)}_i \sigma + \left(1-\tau^{(j)}_i\right)\omega \in {\mathcal{M}}_{1}(2,1) \subset {\mathcal{M}}_{1}(r,1),\] for any \(0 \leq j \leq 2^{\frac{L}{8}}\) and \(1 \leq i \leq L\). Then we may further apply 3 to the tensor product \(\rho^{(j)} = \left(\tau^{(j)}_1 \sigma + \left(1-\tau^{(j)}_1\right)\omega\right) \otimes \cdots \otimes \left(\tau^{(j)}_L \sigma + \left(1-\tau^{(j)}_L\right)\omega\right)\) to deduce that \(\rho^{(j)} \in {\mathcal{M}}_{1 \times n}(r,1) = {\mathcal{M}}_{n}(r,1)\) for any \(0 \leq j \leq 2^{\frac{L}{8}}\), as desired.

Now we proceed to verify that density matrices in the class \({\mathcal{C}}_n\) are pairwisely separated before applying Fano’s method. Specifically, for any \(0 \leq j \neq k \leq 2^{\frac{L}{8}}\), we may compute the Frobenius norm of the difference \(\rho^{(j)} - \rho^{(k)}\) as below: \[\label{eqn:32lower32bound32between32difference32of32density32matrices} \begin{align} \left\|\rho^{(j)} - \rho^{(k)}\right\|^2_{F} &= \left\langle \rho^{(j)} - \rho^{(k)}, \rho^{(j)} - \rho^{(k)} \right\rangle_F \\ &= \left\|\rho^{(j)}\right\|^2_{F} + \left\|\rho^{(k)}\right\|^2_{F} - 2\left\langle \rho^{(j)}, \rho^{(k)}\right\rangle_F = 2- 2\left\langle \rho^{(j)}, \rho^{(k)} \right\rangle_F\\ &= 2-2\prod_{i=1}^{n}\left\langle \tau^{(j)}_i \sigma + \left(1-\tau^{(j)}_i\right)\omega, \tau^{(k)}_i \sigma + \left(1-\tau^{(k)}_i\right)\omega \right\rangle_F\\ &= 2-2\left\langle\sigma, \omega\right\rangle_{F}^{\left\|\tau^{(j)}-\tau^{(k)}\right\|_{1}}. \end{align}\tag{45}\] Plugging in the expressions of \(\sigma\) and \(\omega\) in (41 ) above gives us that \[\label{eqn:32computation32of32basis32dot32product} \left\langle\sigma, \omega\right\rangle_{F} = \frac{1}{2} -\frac{1}{2}\gamma + \frac{1}{2}(1-\gamma) = 1-\gamma\tag{46}\] Then we may substitute (46 ) and the lower bound in (42 ) into (45 ), which yields the following lower bound \[\label{eqn:32lower32bound32on32separation32distance} \begin{align} \left\|\rho^{(j)} - \rho^{(k)}\right\|^2_{F} &\geq 2-2\left(1-\gamma\right)^{\frac{n}{8}} \geq 2 \times \frac{1}{2}\times \gamma \times \frac{n}{8} \geq \frac{\gamma n}{8}, \end{align}\tag{47}\] where the second last inequality follows from the assumption \(B \geq n\) and the fact that \[1-(1-t)^n \geq nt - \frac{n(n-1)}{2}t^2 \geq \frac{1}{2}nt\] holds for any \(n \in {\mathbb{N}}\) and \(t \in (0, \frac{1}{n}]\). Then we can define \(\delta:= \sqrt{\frac{\gamma n}{8}}\) to be the separation distance between any two density matrices.

Furthermore, we now consider the probability distributions associated with the samples obtained by measuring the collection \({\mathcal{C}}_n\) of pairwisely separated density matrices via the classical shadow protocol based on \(\mathcal{U}\). Specifically, we use the binary string \(b^{(j)}_k = \left(b^{(j)}_k(l)\right)_{l=1}^{n} \in \{0,1\}^{n}\) to denote the outcome obtained from the \(k\)-th measurement of the density matrix \(\rho^{(j)}\) under the classical shadow protocol based on \(\mathcal{U}\) for any \(1 \leq k \leq B\). From the tensor product structure of each \(\rho^{(j)}\) in \(\mathcal{C}_n\), we can deduce that the joint distribution formed by \(\left\{b^{(j)}_k\right\}_{k=1}^{B}\) satisfies \[P_{j} := \prod_{k=1}^{B}P_{b^{(j)}_k} = \prod_{k=1}^{B}\prod_{l=1}^{n}P_{b^{(j)}_k(l)}\] for any \(0 \leq j \leq 2^{\frac{L}{8}}\). Then a direct computation further implies \[\label{eqn:32upper32bound32on32KL32between32product32of32dists} \begin{align} D_{\mathrm{KL}}(P_j, P_0) &= D_{\mathrm{KL}}\left(\prod_{k=1}^{B}\prod_{l=1}^{n}P_{b^{(j)}_k(l)}, \prod_{k=1}^{B}\prod_{l=1}^{n}P_{b^{(0)}_k(l)}\right)\\ &= \sum_{k=1}^{B}\sum_{l=1}^{n}D_{\mathrm{KL}}\left(P_{b^{(j)}_k(l)}, P_{b^{(0)}_k(l)}\right). \end{align}\tag{48}\] From our construction of the density matrices \(\rho^{(j)} \;(1 \leq j \leq 2^{\frac{L}{8}})\) given above, we have that each \(P_{b^{(j)}_k(l)}\) must be a Bernoulli distribution. Moreover, we note that the difference between \(\rho^{(j)}\) and \(\rho^{(0)}\) is determined by the difference between \(\sigma\) and \(\omega\) for any \(1 \leq j \leq 2^{\frac{n}{8}}\). By using \(\text{Ber}(p)\) to denote the Bernoulli distribution with mean \(p\) for any \(p \in (0,1)\), we can use the expressions of \(\sigma\) and \(\omega\) given in (41 ) above to deduce that \(D_{\mathrm{KL}}\left(P_{b^{(j)}_k(l)}, P_{b^{(0)}_k(l)}\right)\) can be upper bounded as follows \[\begin{align} D_{\mathrm{KL}}\left(P_{b^{(j)}_k(l)}, P_{b^{(0)}_k(l)}\right) &\leq D_{\mathrm{KL}}\left(\text{Ber}\left(\frac{1+\sqrt{\gamma}}{2}\right), \text{Ber}\left(\frac{1-\sqrt{\gamma}}{2}\right)\right)\\ &= \frac{1+\sqrt{\gamma}}{2}\log\left(\frac{\frac{1+\sqrt{\gamma}}{2}}{\frac{1-\sqrt{\gamma}}{2}}\right) + \frac{1-\sqrt{\gamma}}{2}\log\left(\frac{\frac{1-\sqrt{\gamma}}{2}}{\frac{1+\sqrt{\gamma}}{2}}\right)\\ & = \sqrt{\gamma}\log\left(\frac{1+\sqrt{\gamma}}{1-\sqrt{\gamma}}\right) \leq \frac{2\gamma}{1-\gamma} \leq 4\gamma, \end{align}\] where the second last inequality above follows from the fact that \(\log\left(\frac{1+t}{1-t}\right) \leq \frac{2t}{1-t^2}\) for \(t \in (-1,1)\) and the last inequality above follows from our choice of \(\gamma\). Substituting the upper bound derived above into (48 ) then yields that \[\begin{align} \frac{1}{2^{\frac{n}{8}}}\sum_{j=1}^{2^{\frac{n}{8}}}D_{\mathrm{KL}}(P_j, P_0) \leq 4nB\gamma = \frac{n}{400}. \end{align}\] By letting \(\mathcal{B}:= \{\widehat{\rho}_{(i)}\}_{i=1}^{B}\) be the collection of \(B\) outcomes obtained via the \(B\) measurements under the classical shadow protocol and taking \(\alpha^\ast:= \frac{n}{400}\), \(M:= 2^{\frac{n}{8}}\) here, we can further apply Fano’s method listed in 7 above to deduce the following lower bound on the probability of misclassification for any classifier \(\Theta: \mathcal{B} \rightarrow \{0,1,\cdots,M\}\): \[\label{eqn:32lower32bound32on32misclassification32probability} \begin{align} \max_{0 \leq j \leq M}P_j\left(\Theta(\mathcal{B}) \neq j\right) &\geq \frac{2^{\frac{n}{16}}}{1+2^{\frac{n}{16}}}\left(1 - \frac{3n}{400 \times \frac{n}{8}\log(2)} -\frac{1}{\frac{n}{2}\log(2)}\right) \\ &\geq \frac{1}{2}\left(1- \frac{3}{50\log(2)} - \frac{1}{2\log(2)}\right) \geq \frac{1}{20} \end{align}\tag{49}\] where the second last inequality above follows from the assumption \(n \geq 8\).

Finally, we will prove the desired minimax lower bound by combining the lower bound on the misclassification probability proved in (49 ) above with Markov’s inequality and inequality (47 ) proved earlier. Specifically, for any estimator \(\widehat{\Theta}:\{\widehat{\rho}_{(i)}\}_{i=1}^{B} \rightarrow {\mathbb{C}}^{2^n \times 2^n}\) of the target density matrix, we define a corresponding classifier \(\widetilde{\Theta}:\mathcal{B} \rightarrow \{0,1,\cdots,M\}\) as follows: \[\widetilde{\Theta} \in \mathop{\mathrm{arg\,min}}_{0 \leq j \leq M} \left\|\widehat{\Theta} -\rho^{(j)}\right\|_{F}\] Given that the collection \({\mathcal{C}}_n = \{\rho^{(j)}\}_{j=0}^{M}\) has been proved to be \(\delta\)-seperated in (45 ), we then have the following inequality for any \(0 \leq j \leq M\): \[P_j\left(\widetilde{\Theta}(\mathcal{B}) \neq j\right) \leq P_j\left(\left\|\widehat{\Theta} -\rho^{(j)}\right\|_{F} \geq \frac{\delta}{2}\right)\] from the triangle inequality. Applying Markov’s inequality to the RHS above and the lower bound in (49 ) to the LHS above further implies \[\label{eqn:32application32of32Markov39s32ineq} \begin{align} \max_{0 \leq j \leq M}\frac{\mathbb{E}_{P_j}\left[\left\|\widehat{\Theta} -\rho^{(j)}\right\|_{F}\right]}{\frac{\delta}{2}} &\geq \max_{0 \leq j \leq M}P_j\left(\left\|\widehat{\Theta} -\rho^{(j)}\right\|_{F} \geq \frac{\delta}{2}\right)\\ &\geq \max_{0 \leq j \leq M}P_j\left(\widetilde{\Theta}(\mathcal{B}) \neq j\right) \geq \frac{1}{20}. \end{align}\tag{50}\] Finally, plugging in (50 ) indicates that \[\label{eqn:32final32minimax32lower32bound} \begin{align} \inf_{\widehat{\Theta}}\sup_{\rho \in {\mathcal{M}}_n(r,1)}\mathbb{E}_{\{\widehat{\rho}_{(i)}\}_{i=1}^{B}}\left[\left\|\widehat{\Theta} - \rho\right\|_{F}\right] &\geq \inf_{\widehat{\Theta}}\max_{0 \leq j \leq M}\mathbb{E}_{P_j}\left[\left\|\widehat{\Theta} - \rho^{(j)}\right\|_{F}\right]\\ &\geq \frac{\delta}{40} = \frac{1}{40}\sqrt{\frac{\gamma n}{8}} \gtrsim \sqrt{n\gamma} \gtrsim \sqrt{\frac{n}{B}}, \end{align}\tag{51}\] as desired. This concludes our proof of the information-theoretic lower bound.