Why is the complexity of $n$-qubit state tomography not upper bounded as $O(3^n)$?

Question

Consider the task of fully determining an $n$-qubit state $\rho$ which can be written as

\begin{equation}\tag{1} \rho = \sum_{p \in \{I, X, Y, Z\}^n} \text{Tr}(\rho P_{p}) P_{p} \end{equation}

and each $P_{p} = P_{p_1} \otimes \dots \otimes P_{p_n}$ is a tensor product of Pauli matrices. This suggests that I could perform state tomography evaluate each expectation value $\langle P_p \rangle = \text{Tr}(\rho P_p)$. I would plan on having $3^n$ distinct experimental configurations, one for each combination of local measurement bases $\{X, Y, Z\}^n$.

I thought that the discrepancy between $3^n$ measurement configurations and $4^n-1$ coefficients needed to specify $(1)$ would be resolved because an expectation value of the form $\langle X \otimes I \otimes X \rangle$ could be computed using a marginal distribution over the bitstrings results from the experiment that measures $\langle X \otimes X \otimes X\rangle$ (or the experiments used to compute $\langle X \otimes Y \otimes X\rangle$ or $\langle X \otimes Z \otimes X\rangle$). So any experiment to determine a term $\text{Tr}(\rho P_p)$ Equation $(1)$ where $P_p$ contained an $I$ would be redundant with some other experiment. This is one of the features motivating the method of (Huang, 2021): If you instantiate Theorem 1 therein with $L=4^n$ and $w=n$, it asserts that $4^n$ many $\epsilon$-accurate estimators for Pauli expectation values can be computed in $M = O(n3^n / \epsilon^2)$ total experiments.

But when I look elsewhere in the literature (e.g. Franca, 2017) it suggests that for an arbitrary full-rank, $2^n$-dimensional state $\rho$ you do indeed need $\Omega(4^n)$ measurement configurations for quantum state tomography.

How do I resolve the discrepancy between these two scaling behaviors?

Answer 1

Let's say you have a magical machine that gives you $\langle P_{p} \rangle$ (which are expectation values and therefore, well, numbers) and only the $\langle P_{p} \rangle$. It does this for all the $3$- (n)-fold tensor products of the traceless Paulis. That is: $p \in \{X,Y,Z\}^{\otimes 3 (n)}$, for a total of $3^{3 (n)}$ Paulis, where $n$ is the general case and $3$ is specific to your examples.

The process that you describe to obtain e.g. $\langle X \otimes I \otimes X \rangle$ from, say, $\langle X \otimes X \otimes X \rangle$ does not work here. For that, you need more. These $3^{n}$ values are not enough.

The trick is, the usual method of obtaining these Pauli expectation values is by measuring every qubit in a Pauli basis separately. A single-qubit Pauli measurement is a projection upon either of its eigenspaces. For instance, for the $Z$ operator has projectors $|0\rangle\langle0|$ and $|1\rangle\langle1|$, and the $X$ operator has projectors $|+\rangle\langle+|$ and $|-\rangle\langle-|$. You thus gather (and use) the statistics of six operators per qubit to perform full QST! This results in a total of $6^{n}$ (positive) operators that one uses in standard QST - in the case of the Paulis these are projectors, but the more general POVM can also work.

With all these operators one can reconstruct all Paulis (including the containing-$I$-terms) and therefore reconstruct the density matrix. This is exactly possible because these $6^{n}$ operators form a spanning set of the space of density matrices. The Paulis do too$^{1}$ and therefor we can use either of them reconstruct the density matrix.

But wait - now we have $6^{n}$ terms! We only needed $4^{n}$ terms, right? Well, yes! The space of density matrices dimension' scales as $4^{n}$, so if you can find a set of (positive) operators of size $4^{n}$ that are independent, this should work. A particular nice example is a SIC-POVM; a Symmetric and Informationally-Complete POVM. If you have some magical machine that can perform this SIC-POVM $\{A_{k}\}$ measurement on all your qubits, it would have to perform $4^{n}$ measurements and return single values $\langle A_{k} \rangle = \mathrm{tr}\big[\rho A_{k}\big]$, which would be enough to reconstruct the entire density matrix.

A more straightforward but also easier example of a $4$-partite POVM (actually, PVM) is the set consisting of the $+1$ and $-1$ eigenstates of the $Z$ operator and the $+1$ eigenstates of the $X$ and $Y$ operators: $\{|0\rangle\langle0|,|1\rangle\langle1|,|+\rangle\langle+|,|+i\rangle\langle+i|\}$. You can check that these form a basis for the space of density matrices. However, please realize that to experimentally implement this PVM would be a daunting task (if not completely impossible) without just performing all $3^{n}$ Pauli measurements. But yeah - you could perform Pauli measurements in all $3^{n}$ basis and then only use the statistics for the above projectors, and throw the $-1$ eigenstates of $X$ and $Y$ away!

$^{1}$Actually, they don't. They form a spanning set (and a basis) for the space of Hermitian matrices, of which density matrices are a subset; they're not even density matrices themselves because they're not positive semidefinite nor have trace $1$. But this is not important for the current discussion.

Answer 2

Let $\mathcal X$ be an $N$-dimensional (complex) vector space space. The (real) vector space of Hermitian operators defined on it, $\mathrm{Herm}(\mathcal X)$, has dimension $N^2$. An easy way to see this is to realise that generic complex matrices are characterised by $2N^2$ real parameters, and $A^\dagger=A$ imposes $N^2$ independent constraints. Equivalently, you can notice that the set of complex orthonormal bases also has dimension $N^2$ (as a manifold, not as a vector space), and exploit the close relation between (skew-)Hermitian matrices and unitaries.

The set of density matrices on $\mathcal X$ is a convex subset of the affine subspace of $\mathrm{Herm}(\mathcal X)$ defined by the normalisation condition $\operatorname{Tr}(\rho)=1$. This makes it a subset of an $(N^2-1)$-affine subspace. More specifically, it is the (convex) subset defined by the non-negativity condition $\rho\ge0$.

Applying this to an $n$-qubit system, we have $N=2^n$ and thus the set of states is a convex subset of an affine $((2^n)^2-1)=(4^n-1)$-dimensional space.

The reason you got a wrong result is that you overestimated the contribution of the normalisation to the counting. Expectation values such as $X\otimes I$ are independent from the others. To see it, notice that you can take states of the form $\sum_k c_p P_p$, and arbitrarily change the $X\otimes I$ without touching the other coefficients (assuming all coefficients are small enough that you still get a valid quantum state doing such perturbation).

With regards to the number of needed measurement configurations, one should keep into account that the number of coefficients characterising a state needs not be strictly related to the number of measurement configurations used to do it. As an extreme example, one could use a single measurement configuration to estimate all of the expectation values. To see it, evolve the state through a random isometry that sufficiently enlarges the space, and measure in the outcome computational basis. This is a single measurement setting, but in principle allows you to estimate all the expectation values at once. On the other extreme, if you consider measurement settings with binary outcomes, you will need $\mathcal O(\log_2 N)$ of them.