Density matrix reconstruction using non-negative matrix product states

Quantum state tomography is a key technique for quantum information processing but is challenging due to the exponential growth of its complexity with the system size. In this work we propose an algorithm which iteratively finds the best non-negative matrix product state approximation based on a set of measurement outcomes whose size does not necessarily grow exponentially. Compared to the tomography method based on neural network states, our scheme utilizes a so-called tensor train representation that allows straightforward recovery of the unknown density matrix in the matrix product operator form. As applications, the effectiveness of our algorithm is numerically demonstrated to reconstruct the ground state of the XXZ spin chain under depolarizing noise.

Figure 1

A flowchart of our non-negative tensor train state tomography algorithm. Given an unknown $L$-qubit quantum state ${\widehat{\rho }}_0$, an IC-POVM is performed to obtain a sample of bit strings $\left\{{\mathit{a}}^j\right\}$, each of which is then encoded into a one hot MPS $S^j$. Based on these $S^j$, an optimal non-negative MPS $P\left(\mathit{a}\right)$ with a fixed bond dimension is found through optimization, satisfying that it is closest to $P_s\left(\mathit{a}\right)$, namely, the superposition of all states $S^j$, by optimizing each site tensor of the randomly initialized MPS ansatz $P\left(\mathit{a}\right)$. Then an MPO $\widehat{\rho }$ is reconstructed, with the same bond dimension as $P\left(\mathit{a}\right)$, by applying the inverse of the IC-POVM locally on each site of $P\left(\mathit{a}\right)$ using Eq. (8). Simulations in this work are based on synthetic data.

Figure 2

(a) ${\mathcal{I}}_q$ (red dashed line with squares) and ${\mathcal{I}}_c$ (blue solid line with circles) as functions of system size $L$ for fixed $p=0.6$. The inset shows ${\mathcal{I}}_c$ as a function of system size for larger system size (${\mathcal{I}}_q$ for $L>6$ is not shown since it is too expensive to compute). (b) ${\mathcal{I}}_q$ and ${\mathcal{I}}_c$ as functions of depolarized noise strength $p$ for fixed $L=4$. (c) The minimum number of required training data $N$ as a function of $L$ such that ${\mathcal{I}}_c\le 1\%$, with $p=0.6$. (d) The minimum number of required training data $N$ as a function of $p$ such that ${\mathcal{I}}_c\le 1\%$, with $L=4$. The other parameters used in these simulations are $\gamma =2$ and $D=10$.

Figure 3

(a), (c) ${\mathcal{I}}_q$ (a) and ${\mathcal{I}}_c$ (c) as functions of the number of sweeps for different $\gamma$ with $p=0.4$. (b), (d) ${\mathcal{I}}_q$ (b) and ${\mathcal{I}}_c$ (d) as functions of the number of sweeps for different $\gamma$ with $p=0.6$. The inset in (d) shows the tail of the convergence for ${\mathcal{I}}_c$. In (a), (b), (c), (d) the gray lines from darker to lighter correspond to $\gamma =1,1.2,1.6,2$, respectively. We have also chosen the five best results according to their loss values out of 100 trials and plotted the mean values of them. (The standard deviations are shown as error bars.) (e) The left and right axis show the final ${\mathcal{I}}_q$ and ${\mathcal{I}}_c$ as functions of $\gamma$ with $p=0.4$. (f) The left and right axis show the final ${\mathcal{I}}_q$ and ${\mathcal{I}}_c$ as functions of $\gamma$ with $p=0.6$. Here we have chosen $L=6$. We have also used $N=3\times {10}^7$ and $D=10$ in these simulations.

Figure 4

(a) ${\mathcal{I}}_q$ (red dashed line with squares) and ${\mathcal{I}}_c$ (blue solid line with circles) as functions of the number of training data $N$, with $D=10$. (b) ${\mathcal{I}}_q$ (red dashed line with squares) and ${\mathcal{I}}_c$ (blue solid line with circles) as functions of the bond dimension $D$, with $N=3\times {10}^7$. The other parameters used are $L=6, p=0.6$, and $\gamma =2$.

Figure 5

The $x$ axis denotes different labels of the 100 numerical experiments, labeled 1 to 100, sorted by their final loss values from large to small. In our simulation, we have used $P^2\left(\mathit{a}\right)-2P\left(\mathit{a}\right)P_s\left(\mathit{a}\right)$ as the loss value, which simply shifts the original loss value in Eq. (11) by a constant. Here the results are taken from the reconstruction of $P\left(\mathit{a}\right)$ for $L=6, \gamma =2.0$, and $p=0.6$.