Optimal and two-step adaptive quantum detector tomography

For a $d$-dimensional quantum system, its state can be described by a $d\times d$ Hermitian matrix $\rho$, satisfying $\rho\geq 0,\operatorname{Tr}(\rho)=1$. When $\rho=|\psi\rangle\langle\psi|$ for some $|\psi\rangle\in\mathbb{C}^{d}$, we call $\rho$ a pure state, and its purity $\operatorname{Tr}(\rho^{2})$ reaches the maximum value $1$. Otherwise, $\rho$ is called a mixed state, and can be represented using pure states $\left\{\left|\psi_{i}\right\rangle\right\}:\rho=\sum_{i}c_{i}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$ where $c_{i}\in\mathbb{R}$ and $\sum_{i}c_{i}=1$.

A quantum detector connects the classical and quantum world through a set of operators known as positive-operator-valued measure (POVM). A set of POVM elements is a set of Hermitian and positive semidefinite operators $\left\{P_{i}\right\}$, which is the mathematical representation of quantum detectors, satisfying the completeness constraint $\sum_{i}P_{i}=I$. We directly call $\left\{P_{i}\right\}$ a POVM element in this paper. The operators $P_{i}$ may be finite or infinite dimensional in theory. When infinite dimensional, we usually truncate them at a finite dimension $d$ in practice. When $\rho$ is measured using $\left\{P_{i}\right\}$, the probability of obtaining the $i$-th result is given by the Born’s rule

Because of the completeness constraint, we have $\sum_{i}p_{i}=1$. In practical experiments, suppose $N$ (called resource number) identical copies of $\rho$ are prepared and the $i$-th result occurs $N_{i}$ times. Then $\hat{p}_{i}=N_{i}/N$ is the experimental estimation of the true value $p_{i}$.

By applying known quantum states to an unknown quantum detector and obtaining the measurement results, one can estimate the detector, which is called quantum detector tomography (QDT). Denote the true value of the detector as $\left\{P_{i}\right\}_{i=1}^{n}$. We design $M$ different types of quantum states (called probe states), where each state $\rho_{j}$ use the same resource number $N/M$ thus their total number of copies is $N$. One can formulate the problem of QDT as [28]:

There are various methods to formulate and solve the problem of QDT. For maximum likelihood estimation (MLE) [18], if large data are given, the MLE can asymptotically reach the Cram$\acute{e}$r-Rao bound for parameter estimation while the computational complexity is usually high. There also exist several convex optimization approaches [21, 22, 23] which might be more efficient. However, an analytical error upper bound is missing for MLE and convex optimization approaches, making it difficult to optimize the input probe states in the general case without available prior information about the detector. Here we choose the linear regression method due to its simplicity and computational efficiency. Moreover, Ref. [28] presents the two-stage QDT solution to Problem 1, giving an analytical error upper bound in favor of optimizing general input states. We thus review this solving procedure as follows.

Let $\left\{\Omega_{i}\right\}_{i=1}^{d^{2}}$ be a complete basis set of orthonormal operators with $d$-dimension. Without loss of generality, let $\operatorname{Tr}\left(\Omega_{i}^{\dagger}\Omega_{j}\right)=\delta_{ij}$, $\Omega_{i}=\Omega_{i}^{\dagger}$ where $\operatorname{Tr}\left(\Omega_{i}\right)=0$ except $\Omega_{1}=I/\sqrt{d}$. Then we parameterize the detector and probe states as

where $\theta_{i}^{a}\triangleq\operatorname{Tr}\left(P_{i}\Omega_{a}\right)$ and $\phi_{j}^{b}\triangleq\operatorname{Tr}\left(\rho_{j}\Omega_{b}\right)$ are real. Using Born’s rule, we can obtain

where $\phi_{j}$, $\theta_{i}$ is the parameterization of $\rho_{j}$ and $P_{i}$, respectively. Suppose the outcome for $P_{i}$ appears $n_{ij}$ times, then $\hat{p}_{ij}=n_{ij}/(N/M)$. Denote the estimation error as $e_{ij}=\hat{p}_{ij}-p_{ij}$. According to the central limit theorem, $e_{ij}$ converges in distribution to a normal distribution with mean zero and variance $\left(p_{ij}-p_{ij}^{2}\right)/(N/M)$. We thus have the linear regression equation

Let $\Theta=\left(\theta_{1}^{T},\theta_{2}^{T},\ldots,\theta_{n}^{T}\right)^{T}$ be the vector of all the unknown parameters to be estimated. Collect the parameterization of the probe states as $X=\left(\phi_{1},\phi_{2},\ldots,\phi_{M}\right)^{T}$. Let $\hat{\mathcal{y}}=$ $(\hat{p}_{11},\hat{p}_{12},\ldots,\hat{p}_{1M},\hat{p}_{21},\hat{p}_{22}$, $\ldots,\hat{p}_{2M},\ldots,\hat{p}_{nM})^{T}$, $\mathcal{X}=I_{n}\otimes X$, $\mathcal{e}=(e_{11},e_{12},\ldots,e_{1M}$, $e_{21},e_{22},\ldots,e_{2M},\ldots,e_{nM})^{T}$, $\mathcal{H}=(1,1,\ldots,1)_{1\times n}\otimes I_{d^{2}},\mathcal{d}_{d^{2}\times 1}=(\sqrt{d},0,\ldots,0)^{T}$. Then the regression equations can be rewritten in a compact form [28]:

with a linear constraint

Now Problem 1 can be transformed into the following equivalent form:

In [28], Problem 2 is split into two approximate subproblems:

The reconstruction algorithm has two stages, as shown in Fig. 1. In Stage 1, Problem 2.1 is solved by Constrained Least Squares (CLS) method and the form of the solution is further simplified. In Stage 2, Problem 2.2 is solved by matrix decomposition. We briefly review the two-stage QDT reconstruction algorithm in [28] in the next subsection.

In Stage 1, we directly give the simplified form of the CLS solution to Problem 2.1 as

where $\hat{y}_{i}=\left(\hat{p}_{i1},\hat{p}_{i2},\ldots,\hat{p}_{iM}\right)^{T}$ for $1\leq i\leq n$ and $y_{0}=\left((1,\ldots,1)_{1\times M}\right)^{T}=\sum_{i}\hat{y}_{i}$. Let the $i$-th block be $\hat{\Theta}_{i,\text{CLS}}=\left(\hat{\theta}_{i}^{1},\ldots,\hat{\theta}_{i}^{d^{2}}\right)^{T}$ and the Stage 1 estimate is $\hat{E}_{i}=\sum_{a=1}^{d^{2}}\hat{\theta}_{i}^{a}\Omega_{a}$ which may have negative eigenvalues. Thus, Stage 2 is designed to obtain a physical estimate $\hat{P_{i}}$ by eigenvalue correction with three substages. Firstly, we decompose $\hat{E}_{i}=\hat{F}_{i}-\hat{G}_{i}$ with $\hat{F}_{i}\geq 0,\hat{G}_{i}\geq 0$. We then perform a spectral decomposition to obtain $\hat{E}_{i}=\hat{R}_{i}\hat{K}_{i}\hat{R}_{i}^{\dagger}$. Assume there are $\hat{n}_{i}$ nonpositive eigenvalues for $\hat{E}_{i}$, and they are in decreasing order. Thus, the best decomposition in the sense of minimizing $\|\hat{G}_{i}\|$ is $\hat{F}_{i}=\hat{R}_{i}\text{diag}((\hat{K}_{i})_{11},(\hat{K}_{i})_{22},...,(\hat{K}_{i})_{(d-\hat{n}_{i})(d-\hat{n}_{i})},0,...,0)\hat{R}_{i}^{\dagger}$, and $\hat{G}_{i}=-\hat{R}_{i}\text{diag}(0,...,0,(\hat{K}_{i})_{(d-\hat{n}_{i}+1)(d-\hat{n}_{i}+1)}$,
$(\hat{K}_{i})_{(d-\hat{n}_{i}+2)(d-\hat{n}_{i}+2)},...,(\hat{K}_{i})_{dd})\hat{R}_{i}^{\dagger}$. Secondly, we apply decomposition $I+\sum_{i}\hat{G}_{i}=\sum_{i}\hat{F}_{i}=\hat{C}\hat{C}^{\dagger}$, and hence $\sum_{i}\hat{C}^{-1}\hat{F}_{i}\hat{C}^{-\dagger}=I$. Let $\hat{B}_{i}=\hat{C}^{-1}\hat{F}_{i}\hat{C}^{-\dagger}$ which is positive semidefinite and $\sum_{i}\hat{B}_{i}=I$. For any unitary $\hat{U}$, $\hat{C}\hat{C}^{\dagger}=\hat{C}\hat{U}\hat{U}^{\dagger}\hat{C}^{\dagger}$ holds. Therefore, $\hat{U}^{\dagger}\hat{B}_{i}\hat{U}$ can also be an estimate of the detector. To neutralize the effect of $\hat{U}^{\dagger}\hat{C}^{-1}(\cdot)(\hat{C}^{\dagger})^{-1}\hat{U}$ on $\hat{F}_{i}$, we define an optimal unitary matrix $\hat{U}_{\text{op}}$ to minimize $\|\hat{C}\hat{U}_{\text{op}}-I\|$ [28] and its analytical expression is $\hat{U}_{\text{op}}=\sqrt{\hat{C}^{\dagger}\hat{C}}\hat{C}^{-1}$. Therefore, the final estimate is $\hat{P}_{i}=\hat{U}_{\text{op}}^{\dagger}\hat{B}_{i}\hat{U}_{\text{op}}$. This general algorithm can be used for arbitrary $n\geq 2$.

One index to evaluate probe states is to score their worst performance, i.e., their upper bounds of estimation errors. Let $\mathbb{E}$ denote the expectation w.r.t. all possible measurement results. The Stage 1 error $\left\|\hat{E}_{i}-P_{i}\right\|$ between the estimate $\hat{E}_{i}$ and its true value $P_{i}$ is bounded by UMSE (upper bound of the mean squared error ) [28]

and the final estimation error for all detectors
$\mathbb{E}\left(\Sigma_{i}\left\|\hat{P}_{i}-P_{i}\right\|^{2}\right)$ is bounded by [28]

Since in (8) and (9) the parts dependent on probe states are both $M\operatorname{Tr}\big{[}(X^{T}X)^{-1}\big{]}$, we take it as our first criterion. For general QDT, we have no special prior knowledge on the detectors, and the following conditions are equivalent: (i) the detector can not be uniquely identified; (ii) the probe states do not span the space of all $d$-dimensional states; (iii) $X^{T}X$ is singular; (iv) the UMSE is infinite. From (i) and (iv), it is thus reasonable to take $M\operatorname{Tr}(X^{T}X)^{-1}$ as an evaluation index. We require the optimal probe states to minimize $M\operatorname{Tr}(X^{T}X)^{-1}$. This criterion is conservative and its limitation is that the UMSE might not be tight. The relation (8) is only tight for the case $P_{1}=P_{2}=I/2$ and (9) is always loose. Therefore, minimizing UMSE is not equivalent to minimizing MSE. In addition, UMSE depends on the assumption that there only exists statistics error in the measurement data from finite state copies (discussed in Remark 4). Despite these limitations, numerical results in [28] (e.g., Fig. 4 therein) indicate a very similar behavior between UMSE and MSE when changing the probe states while maintaining the other parameters $n$, $M$, $N$ fixed. Hence, UMSE is also a useful index evaluating the probe state set.

The second index we consider is, for a set of probe states, how robust the generated estimation result is w.r.t. measurement noise. Note that our CLS estimation (7) is in fact equivalent to the least squares estimation of the linear regression problem

Typically, the sensitivity of a linear system solution to perturbations in the data is evaluated by the condition number of the coefficient matrix, defined (among several possible choices) in this paper as $\operatorname{cond}(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}$, where $\sigma_{\max/\min}(A)$ is the maximum/minimum singular value of $A$. Hence, to maximize the estimation robustness w.r.t. measurement noise amounts to minimizing the condition number in (10) $\operatorname{cond}(I_{n}\otimes X)=\operatorname{cond}(I_{n})\operatorname{cond}(X)=\operatorname{cond}(X)$. Therefore, the second evaluation index is chosen as $\operatorname{cond}(X)$. To sum up, the optimal probe states should minimize $M\operatorname{Tr}(X^{T}X)^{-1}$ and minimize $\operatorname{cond}(X)$ simultaneously. In the following we give specific characterization and show that the two optimal indices can be achieved simultaneously for several examples.

We now give the condition on optimal probe states (OPS).

We assume that $M\geq d^{2}$ different types of probe states are used, since clearly $M<d^{2}$ is not optimal. Denote the eigenvalues of $X^{T}X$ as $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d^{2}}$. For minimizing UMSE, it can be formulated as

For maximizing robustness, it can be formulated as

The constraint $\sum_{i=1}^{d^{2}}\lambda_{i}\leq M$ comes from the purity requirement. For the parameterization of each probe state, we have

Thus for the parameterization matrix of probe states, we have

The second constraint $\lambda_{1}\geq\frac{M}{d}$ is from the matrix element satisfying

in our chosen basis $\left\{\Omega_{i}\right\}_{i=1}^{d^{2}}$.

For (11) and (12), using Lagrange multiplier method as in [4], the minimum of $\sum_{i=1}^{d^{2}}\frac{M}{\lambda_{i}}$ is $d^{4}+d^{3}-d^{2}$, achieved when $\lambda_{1}=\frac{M}{d}$ and $\lambda_{2}=\dots=\lambda_{d^{2}}=\frac{M}{d(d+1)}$. Hence, the minimum of UMSE is

which does not rely on the type number of probe states $M$. The minimum of $\sqrt{\frac{\lambda_{1}}{\lambda_{d^{2}}}}$ is $\sqrt{d+1}$ because $\lambda_{1}\geq\frac{M}{d}$ and $\lambda_{d^{2}}\leq\frac{M}{d(d+1)}$. These minima can be achieved at the same time if and only if $\lambda_{1}=\frac{M}{d}$ and $\lambda_{2}=\dots=\lambda_{d^{2}}=\frac{M}{d(d+1)}$.

Now, we show that $X^{T}X$ is a diagonal matrix. Assume $X^{T}X=V\operatorname{diag}\left(\frac{M}{d},\frac{M}{d(d+1)}I_{d-1}\right)V^{T}$. For $\left(X^{T}X\right)_{11}$, we have

where $v=\left[1,0,\cdots,0\right]^{T}$. Thus, $V_{11}=1$ and $V=\operatorname{diag}(1,\tilde{V}_{d-1})$ where $\tilde{V}_{d-1}$ is a $d-1$ dimensional orthogonal matrix. Therefore, we have

which is a diagonal matrix. $\Box$

Our results also show that the description of OPS in Assumption $1$ in [28] needs to be more precise. For certain $M$, if there exist OPS, all of them must be pure states because $\sum_{j=1}^{M}\operatorname{Tr}(\rho_{j}^{2})=\sum_{j=1}^{d^{2}}\lambda_{i}=M$ from (14). This indicates that pure states are better than mixed states. However, it is not clear whether these OPS exist for arbitrary $M\geq d^{2}$. Thus, whether we can always find a pure state set which is better than arbitrary mixed state set for given $M$ is still an open problem.

We then give two examples of OPS for $M=d^{2}$ and $M=d(d+1)$ which are motivated by projection measurements.

The first example is motivated from SIC-POVM. To reconstruct an unknown quantum state $\rho$, a generalized measurement must have at least $d^{2}$ linear independent elements, which is called informationally complete. Furthermore, if the measurement results are maximally independent, the POVM is called symmetric informationally complete POVM (SIC-POVM). The simplest mathematical definition of an SIC-POVM is a set of $d^{2}$ normalized vectors $\left|\phi_{k}\right\rangle$ in $\mathbb{C}^{d}$ satisfying [40]

It has been conjectured that SIC-POVMs exist for all dimensions [40] and their existence has been given for $d\leq 121$, and some other sporadic values [41]. When SIC-POVMs exist in $d$ dimension, we define the $d^{2}$ pure states $|\phi_{k}\rangle\in\mathbb{C}^{d}(1\leq k\leq d^{2})$ to be SIC states if they satisfy (19).

Note that choosing different orthogonal bases does not change the inner product between vectors. Therefore, for parameterization of the SIC states, we have

and $M=d^{2}$. We assume the standard singular value decomposition of $X$ is $X=U\Sigma V^{T}$ where $\Sigma$ is a diagonal matrix, and $U$ and $V$ are unitary matrices. According to (20), we have

Thus, $\Sigma=\operatorname{diag}(d,{\frac{d}{d+1}I_{d^{2}-1}})$ and

The largest eigenvalue of $X^{T}X$ is $d$ and the other eigenvalues are $\frac{d}{d+1}$, which satisfies Theorem 1. Therefore, SIC states are OPS. Also, among all $d$-dimensional OPS, SIC states have the smallest $M$ as $M=d^{2}$. If $M<d^{2}$, it is not informationally complete, and the detector cannot be uniquely determined without other prior knowledge. $\Box$ The second example is motivated from MUB measurement. Two sets of orthogonal bases $\mathcal{B}^{k}=\{\left|\psi_{i}^{k}\right\rangle:i=1,\ldots$$,d\}$ and $\mathcal{B}^{\ell}=\left\{\left|\psi_{j}^{\ell}\right\rangle:j=1,\ldots,d\right\}$ are called mutually unbiased if and only if [42]

In particular, one can find maximally $d+1$ sets of mutually unbiased bases in Hilbert spaces of prime-power dimension $d=p^{k},$ with $p$ a prime and ${k}$ a positive integer [43]. When $d+1$ sets of MUB measurement exist in $\mathbb{C}^{d}$, we view each projective MUB measurement operator as a pure state and we call them MUB states. Thus, MUB states always exist for $m$-qubit systems $(d=2^{m})$.

For MUB states in (23), we have $M=d(d+1)$. For their parameterization, we have

We assume the standard singular value decomposition of $X$ is $X=U\Sigma V^{T}$. Thus, we have

where $\left(\frac{1}{d}\right)_{d}$ denotes a $d\times d$ matrix and all the elements are $\frac{1}{d}$. We have

where the largest eigenvalue is $d+1$ and the remaining $d^{2}-1$ eigenvalues are $1$. They are OPS because their eigenvalues satisfy Theorem 1. $\Box$ For $M\neq d^{2}$ and $M\neq d(d+1)$, we leave it an open problem when there always exist OPS satisfying the two indices simultaneously. For two-qubit detectors, from the above results we know the optimal probe states can be constructed using $4$-dimensional MUB states and SIC states as shown in Appendix A. A similar two-qubit problem for QST was discussed in [4] to determine the optimal measurement based on UMSE and in [39] based on condition number.

For one-qubit probe states, they have simple geometric property; i.e., they are in the Bloch sphere. The pure states are on the surface and the mixed states are inside the sphere. Hence, an alternative method to search for one-qubit OPS is based on geometric symmetry. For $M=4$, when the probe states are on the surface of Bloch sphere and become the four vertices of a tetrahedron concentric with Bloch sphere, they are OPS. In fact, they are also SIC states. For $M=6$, one can construct OPS similarly using octahedron, and they are also MUB states. We also have cube, icosahedron, dodecahedron for $M=8,12,20$, respectively, which are OPS. We conjecture all one-qubit OPS are constructed from the five platonic solids on the Bloch sphere in this way, and we show there do not exist OPS for $M=5$ in Appendix B. For multi-qubit states, we still do not fully know the geometric property and it is thus an open problem to find other optimal probe states.

In experiment, product states are among the ones most straightforwardly to be implemented. In this subsection we consider the case $d=2^{m}$ for some integer $m$ and each probe state is an $m$-qubit tensor product state as $\rho_{j}=\rho^{(1)}_{j_{1}}\otimes\cdots\otimes\rho^{(m)}_{j_{m}}$ where ${1\leq j_{1}\leq M_{1},\cdots,1\leq j_{m}\leq M_{m}}$ and there are $M_{i}$ different types of one-qubits for the $i$-th qubit of the product states. The total number of these $m$-qubit tensor product states is thus $\prod_{i=1}^{m}{M_{i}}$. We then give their optimal UMSE and condition number.

Assume the parameterization of all the $i$-th single-qubit states $\left\{\rho_{j_{i}}^{(i)}\right\}_{j_{i}=1}^{M_{i}}$ is $X_{i}$. We thus have

and therefore,

where $20$ comes from $d^{4}+d^{3}-d^{2}$ for $d=2$. Thus, the UMSE is $\frac{(n-1)M}{4N}\operatorname{Tr}\left[\left(X^{T}X\right)^{-1}\right]\geq\frac{20^{m}(n-1)}{4N}$. For condition number,

where $3$ comes from $d+1$ for $d=2$. We can reach these minima if and only if for each $i$, $\left\{\rho_{j_{i}}^{(i)}\right\}_{j_{i}=1}^{M_{i}}$ is an OPS set. $\Box$

Here we compare the two criteria–UMSE and condition number between OPS and optimal product states for an $m$-qubit detector with the dimension $d=2^{m}$. The UMSEs are $\frac{(n-1)(16^{m}+8^{m}-4^{m})}{4N}$ for OPS and $\frac{20^{m}(n-1)}{4N}$ for optimal product states. UMSE of optimal probe states is always smaller than that of optimal product states for $m\geq 2$. For $m=1$, they are both $\frac{20(n-1)}{4N}$. For condition number, it is $\sqrt{2^{m}+1}$ for OPS and $\sqrt{3^{m}}$ for optimal product states. Thus, for both UMSE and condition number, OPS play better than optimal product states in multi-qubit systems. Because most of the OPS are entangled states, this is an example showcasing the advantage of entanglement in QDT.

We give an example of two-qubit product probe states. Let the parameterization of the $j$-th two-qubit product probe state $\rho_{j}$ be

where $\phi_{j_{i}}^{(i)}=\left[\frac{1}{\sqrt{2}},\psi_{j_{i}}^{(i)}\right]^{T}$ is the parameterization of the $i$-th qubit of $\rho_{j}$. Since optimal one-qubit probe states are pure states, $\left\|\psi_{j_{1}}^{(1)}\right\|=\left\|\psi_{j_{2}}^{(2)}\right\|=\frac{1}{\sqrt{2}}$. Denote the eigenvalues of $X^{T}X$ as $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d^{2}}$. We can obtain three constraints as

Thus, the optimization problem is

It can be proven that $\sum_{i=1}^{16}\frac{M}{\lambda_{i}}$ reaches its minimum $400$ and the minimum UMSE is $\frac{100(n-1)}{N}$ when $\lambda_{1}=\frac{M}{4}$$,\lambda_{2}=\cdots=\lambda_{7}=\frac{M}{12}$, $\lambda_{8}=\cdots=\lambda_{16}=\frac{M}{36}$. It can be verified that this minimum UMSE can be reached by using the tensor of platonic solid states such as $M=36$ (Cube). The minimum condition number is $\sqrt{\frac{\lambda_{1}}{\lambda_{16}}}=\sqrt{\frac{M/4}{M/36}}=3$.