Two-stage Estimation for Quantum Detector Tomography: Error Analysis, Numerical and Experimental Results

TL;DR

This paper introduces a two-stage quantum detector tomography method that improves calibration accuracy by combining linear regression with physical constraints, validated through simulations and experiments.

Contribution

A novel two-stage estimation approach for quantum detector tomography that enhances accuracy and guarantees physical validity, with detailed error analysis and optimization.

Findings

Error upper bound established for the method

Simulation results demonstrate improved accuracy

Experimental validation confirms effectiveness

Abstract

Quantum detector tomography is a fundamental technique for calibrating quantum devices and performing quantum engineering tasks. In this paper, a novel quantum detector tomography method is proposed. First, a series of different probe states are used to generate measurement data. Then, using constrained linear regression estimation, a stage-1 estimation of the detector is obtained. Finally, the positive semidefinite requirement is added to guarantee a physical stage-2 estimation. This Two-stage Estimation (TSE) method has computational complexity $O(nd^2M)$, where $n$ is the number of $d$-dimensional detector matrices and $M$ is the number of different probe states. An error upper bound is established, and optimization on the coherent probe states is investigated. We perform simulation and a quantum optical experiment to testify the effectiveness of the TSE method.

Equations

$$\text{i}\hslash \frac{d}{dt}|\psi (t)⟩=\mathcal{H}|\psi (t)⟩,$$

$$p_i=\text{Tr}(P_i\rho ).$$

$$P_i=\sum _{a=1}^{d^2}{\theta }_a^{(i)}{\Omega }_a,$$

$${\rho }_j=\sum _{b=1}^{d^2}{\varphi }_b^{(j)}{\Omega }_b,$$

$$p_{ij}=\text{Tr}(P_i{\rho }_j).$$

$$p_{ij}=\sum _{a=1}^{d^2}{\varphi }_a^{(j)}{\theta }_a^{(i)}={\Phi }_j^T{\Theta }_i,$$

$${\hat{p}}_{ij}={\Phi }_j^T{\Theta }_i+e_{ij}.$$

$$\hat{Y}=X\Theta +e,$$

$$H\Theta =D.$$

$${\hat{\Theta }}_{CLS}={\hat{\Theta }}_{LS}-(X^{T}X{)}^{-1}{H}^T[H(X^{T}X{)}^{-1}{H}^T{]}^{-1}(H{\hat{\Theta }}_{LS}-D),$$

$${\hat{\Theta }}_{LS}=(X^{T}X{)}^{-1}{X}^T\hat{Y}.$$

$$[H(X^{T}X{)}^{-1}{H}^T{]}^{-1}=[H(I_n\otimes Z_0)H^T{]}^{-1}=(n{Z}_0{)}^{-1}=\frac{1}{n}{Z}_0^{-1}.$$

$${\hat{\Theta }}_{LS}=(I_n\otimes Z_0)(I_n\otimes X_0^T)\hat{Y}=(I_n\otimes Z_0{X}_0^T)\hat{Y}.$$

$$\begin{array}[]{rl}&\ \ \ \ \hat{\Theta}_{CLS}\\ &=(I_{n}\otimes Z_{0}X_{0}^{T})\hat{Y}-(I_{n}\otimes Z_{0})\left(\begin{matrix}I_{d^{2}}\\ \vdots\\ I_{d^{2}}\\ \end{matrix}\right)\frac{1}{n}Z_{0}^{-1}\\ &\ \ \ \ \cdot[(I_{d^{2}}\cdots I_{d^{2}})(I_{n}\otimes Z_{0}X_{0}^{T})\hat{Y}-D]\\ &=(I_{n}\otimes Z_{0}X_{0}^{T})\hat{Y}-\frac{1}{n}\left(\begin{matrix}Z_{0}\\ \vdots\\ Z_{0}\\ \end{matrix}\right)Z_{0}^{-1}[(Z_{0}X_{0}^{T}\cdots Z_{0}X_{0}^{T})\hat{Y}-D]\\ &=(I_{n}\otimes Z_{0}X_{0}^{T})\hat{Y}-\frac{1}{n}\left(\begin{matrix}I_{d^{2}}\\ \vdots\\ I_{d^{2}}\\ \end{matrix}\right)[(Z_{0}X_{0}^{T}\cdots Z_{0}X_{0}^{T})\hat{Y}-D]\\ &=(I_{n}\otimes Z_{0}X_{0}^{T})\hat{Y}-\frac{1}{n}\left(\begin{matrix}Z_{0}X_{0}^{T}&\cdots&Z_{0}X_{0}^{T}\\ \vdots&&\vdots\\ Z_{0}X_{0}^{T}&\cdots&Z_{0}X_{0}^{T}\\ \end{matrix}\right)\hat{Y}+\frac{1}{n}\left(\begin{matrix}D\\ \vdots\\ D\\ \end{matrix}\right).\\ \end{array}$$

$$\begin{array}[]{rl}\hat{\Theta}_{CLS}&=\left(\begin{matrix}Z_{0}X_{0}^{T}&&\\ &\ddots&\\ &&Z_{0}X_{0}^{T}\\ \end{matrix}\right)\left(\begin{matrix}\hat{Y}_{1}\\ \vdots\\ \hat{Y}_{n}\\ \end{matrix}\right)\\ &\ \ \ \ -\frac{1}{n}\left(\begin{matrix}Z_{0}X_{0}^{T}&\cdots&Z_{0}X_{0}^{T}\\ \vdots&&\vdots\\ Z_{0}X_{0}^{T}&\cdots&Z_{0}X_{0}^{T}\\ \end{matrix}\right)\left(\begin{matrix}\hat{Y}_{1}\\ \vdots\\ \hat{Y}_{n}\\ \end{matrix}\right)+\frac{1}{n}\left(\begin{matrix}D\\ \vdots\\ D\\ \end{matrix}\right)\\ &=\left(\begin{matrix}Z_{0}X_{0}^{T}\hat{Y}_{1}\\ \vdots\\ Z_{0}X_{0}^{T}\hat{Y}_{n}\\ \end{matrix}\right)-\frac{1}{n}\left(\begin{matrix}Z_{0}X_{0}^{T}\sum_{i}\hat{Y}_{i}\\ \vdots\\ Z_{0}X_{0}^{T}\sum_{i}\hat{Y}_{i}\\ \end{matrix}\right)+\frac{1}{n}\left(\begin{matrix}D\\ \vdots\\ D\\ \end{matrix}\right)\\ &=\left(\begin{matrix}(X_{0}^{T}X_{0})^{-1}X_{0}^{T}(\hat{Y}_{1}-\frac{1}{n}Y_{0})+\frac{1}{n}D\\ \vdots\\ (X_{0}^{T}X_{0})^{-1}X_{0}^{T}(\hat{Y}_{n}-\frac{1}{n}Y_{0})+\frac{1}{n}D\\ \end{matrix}\right).\\ \end{array}$$

$${\hat{K}}_i={\hat{W}}_i^{†}{\hat{F}}_i{\hat{W}}_i-{\hat{W}}_i^{†}{\hat{G}}_i{\hat{W}}_i.$$

$$\min ||{\hat{G}}_i|{|}^2=\min ||{\hat{W}}_i^{†}{\hat{G}}_i^o{\hat{W}}_i|{|}^2=\sum _j\min ({\hat{W}}_i^{†}{\hat{G}}_i^o{\hat{W}}_i{)}_{jj}^2,$$

$${\hat{F}}_i={\hat{W}}_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{W}}_i^{†}$$

$$\begin{array}[]{rl}\hat{G}_{i}&=-\hat{W}_{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{W}_{i}^{\dagger}.\\ \end{array}$$

$$I+\sum _i{\hat{G}}_i=\sum _i{\hat{F}}_i.$$

$$I+\sum _i{\hat{G}}_i=\hat{C}{\hat{C}}^{†}.$$

$$\sum _i{\hat{C}}^{-1}{\hat{F}}_i({\hat{C}}^{†}{)}^{-1}=I.$$

$$\begin{array}[]{rl}||\hat{C}\hat{U}-I||^{2}&=\text{Tr}[(\hat{C}\hat{U}-I)(\hat{U}^{\dagger}\hat{C}^{\dagger}-I)]\\ &=d+\text{Tr}(\hat{C}\hat{C}^{\dagger})-\text{Tr}(\hat{C}\hat{U}+\hat{U}^{\dagger}\hat{C}^{\dagger}).\\ \end{array}$$

$$\frac{\partial L}{\partial {\hat{U}}^{\ast }}=-{\hat{C}}^{†}+(\Lambda +{\Lambda }^{†})\hat{U}=0.$$

$${\hat{C}}^{†}{\hat{U}}^{†}=\Lambda +{\Lambda }^{†}=\hat{U}\hat{C}.$$

$$\sqrt{\hat{C}{\hat{C}}^{†}}{\hat{C}}^{-1}={\hat{U}}_{\beta }{\hat{U}}_{\alpha }^{†}.$$

$$\hat{S}{\hat{U}}_{\gamma }^{†}={\hat{U}}_{\gamma }\hat{S}.$$

$${\hat{U}}_{\gamma }{\hat{S}}^2{\hat{U}}_{\gamma }^{†}={\hat{U}}_{\gamma }\hat{S}{\hat{U}}_{\gamma }\hat{S}=\hat{S}{\hat{U}}_{\gamma }^{†}{\hat{U}}_{\gamma }\hat{S}={\hat{S}}^2.$$

$$\hat{U}={\hat{U}}_{\beta }{\hat{U}}_{\alpha }^{†}=\sqrt{{\hat{C}}^{†}\hat{C}}{\hat{C}}^{-1}.$$

$$\frac{M}{4N}\text{Tr}[(X_0^T{X}_0{)}^{-1}]=O(\frac{d^4}{N}).$$

Full text

I Introduction

Recent decades have witnessed the fast development of quantum science and technology, including quantum computation, quantum communication [1] , quantum sensing [2], etc. Measurement, on a quantum entity or using a quantum object, is the connection between the classical (non-quantum) world and the quantum domain, and plays a fundamental role in investigating and controlling a quantum system [3, 4]. For example, quantum computation can be performed through a series of appropriate measurements in certain schemes [5]. In quantum communication, measurement is a vital part of quantum key distribution [6]. In quantum metrology, adaptive measurement can achieve the Heisenberg limit in phase estimation [7].

Since quantum measurement can also be viewed as a class of quantum resource, its investigation and characterization is fundamentally important [8]. Quantum detector tomography is a technique to characterize quantum measurement devices [9, 10], and thus paves the way for other estimation tasks like quantum state tomography [11]-[16], Hamiltonian identification [17]-[20] and quantum process tomography [22]-[24].

The investigation of protocols for quantum detector tomography dates back to [25], where the Maximum Likelihood Estimation (MLE) method is employed to reconstruct an unknown POVM detector. As one of the most widely recognized methods [11, 26], MLE can preserve the positivity and completeness of the detector, but it is difficult to characterize the error and computational complexity. Phase-insensitive detectors correspond to diagonal matrices in the photon number basis and are thus relatively straightforward to be reconstructed. Ref. [27] modelled this problem as a linear-regression problem and obtained a least squares solution. In [28, 29], phase-insensitive detector tomography was modelled as a convex quadratic optimization problem and an efficient numerical solution was obtained. This method was also experimentally tested in [30, 31], and then was developed in [32] and [33] to model phase-sensitive detector tomography as a recursive constrained convex optimization problem, where the unknown parameters are recursively estimated. For phase-insensitive detectors with a large linear loss, an extension of detector tomography was introduced in [34] and tested on a superconducting multiphoton nanodetector.

In this paper, we propose a novel quantum detector tomography protocol, which is applicable to both phase-insensitive and general phase-sensitive detectors. We first input a series of different states (probe states) to the detector, and then collect all the measurement data. The forthcoming algorithm mainly consists of two stages: in the first stage, we find a constrained least square estimate, which corresponds to a Hermitian estimate satisfying the completeness constraint. However, this estimate can be non-physical; i.e., the estimated detectors may have negative eigenvalues. Hence, in the second stage we further design a series of matrix transformations preserving the Hermitian and completeness constraint to find a physical approximation based on the result in the first stage, and thus obtain the final physical estimate. Our Two-stage Estimation (TSE) method has computational complexity $O(nd^{2}M)$, where $n$ and $d$ are the number and dimension of the detector matrices, respectively, and $M$ is the number of different probe states. This theoretical characterization of the computational complexity is not common in other detector tomography methods. We further prove an error upper bound $O(\frac{d^{5}n^{2}}{N})$ on the condition that the probe states are optimal (if not optimal, the specific form of the bound is also given in Sec. IV), where $N$ is the total copy number of probe states. We then investigated optimization of the types of coherent probe states and the size of their sampling square. We perform numerical simulation to validate the theoretical analysis and compare our algorithm with MLE method. Finally, we slightly modify our method to cater to a practical experiment situation, and we perform quantum optical experiments using two-mode coherent states to testify the effectiveness of our method.

This paper is organized as follows. In Section II, we introduce some preliminary knowledge about quantum physics and formulate our estimation problem. In Section III, we present the procedures of our TSE method and analyze the computational complexity. An upper bound for the estimation error of TSE is established in Section IV. Section V investigates the optimization of the coherent probe states. Section VI presents the numerical simulation results to verify the theoretical analysis in Section IV and Section V, and to compare our method with MLE. Section VII modifies the TSE method according to our practical physical setting and presents the experimental results. Section VIII concludes this paper.

Notation: $A\geq 0$ means $A$ is positive semidefinite. $A^{\dagger}$ is the conjugation ($*$) and transpose ($T$) of $A$. $I$ is the identity matrix. $\mathbb{R}$ and $\mathbb{C}$ are the real and complex domains, respectively. $\otimes$ is the tensor product. $\oplus$ is the matrix direct sum. vec is the column vectorization function. $||\cdot||$ is the Frobenius norm. $\delta$ is the Kronecker delta function. $\text{i}=\sqrt{-1}$. $\text{diag}(X)$ has two effects: it outputs a diagonal matrix with the diagonal elements being the elements in $X$ if $X$ is a vector, or sets all the non-diagonal elements in $X$ to be zero if $X$ is a square matrix. $\hat{X}$ denotes the estimation of variable $X$. For any positive semidefinite $X_{d\times d}$ with spectral decomposition $X=UPU^{\dagger}$, define $\sqrt{X}$ or $X^{\frac{1}{2}}$ as $U\text{diag}(\sqrt{P_{11}},\sqrt{P_{22}},...,\sqrt{P_{dd}})U^{\dagger}$.

II Preliminaries and Problem formulation

II-A Quantum state and measurement

For a $d$-dimensional quantum system, its state is usually described by a $d\times d$ Hermitian matrix $\rho$, which should be positive semidefinite and satisfy $\text{Tr}(\rho)=1$. When $\rho$ is a pure state (satisfying $\text{Tr}(\rho^{2})=1$), we have $\rho=|\psi\rangle\langle\psi|$ where $|\psi\rangle$ is a complex vector on the $d$-dimensional underlying Hilbert space. In this case, we usually identify $\rho$ with $|\psi\rangle$. Otherwise, $\rho$ is called a mixed state, and can be expanded using pure states $\{|\psi_{i}\rangle\}$: $\rho=\sum_{i}c_{i}|\psi_{i}\rangle\langle\psi_{i}|$ where $c_{i}\in\mathbb{R}$ and $\sum_{i}c_{i}=1$. The evolution of a pure state $|\psi\rangle$ is described by the Schrödinger equation

[TABLE]

where $\hbar$ is the reduced Planck constant and $\mathcal{H}$ is the system Hamiltonian.

One of the most common quantum measurement methods is the positive operator valued measure (POVM), and quantum detectors are devices to realize a POVM, especially in the optical domain. A set of POVM elements is a set of operators $\{P_{i}\}$ satisfying the completeness constraint $\sum_{i}P_{i}=I$ and each $P_{i}$ is Hermitian and positive semidefinite. In the case when each operator $P_{i}$ is infinite dimensional, they are usually truncated at a finite dimension $d$ in practice. When the measurements corresponding to operators $\{P_{i}\}$ are performed on $\rho$, the probability of obtaining the $i$-th result is given by the Born Rule

[TABLE]

From the completeness constraint, we thus have $\sum_{i}p_{i}=1$. In practical experiments, suppose that $N$ (also called the resource number) identical copies of $\rho$ are prepared and the $i$-th results occur $N_{i}$ times. Then $N_{i}/N$ is the experimental estimation of the true value $p_{i}$. The measurement apparatus is the physical realization of a quantum detector, and $\{P_{i}\}$ is the mathematical representation. We thus directly call $\{P_{i}\}$ a quantum detector in this paper.

II-B Problem formulation

The technique to deduce an unknown detector from known quantum states and measurement results is called quantum detector tomography. Suppose the true values of a set for a detector are $\{P_{i}\}_{i=1}^{n}$ such that $\sum_{i=1}^{n}P_{i}=I$ with each $P_{i}$ Hermitian and positive semidefinite. We design a series of different quantum states $\rho_{j}$ (called probe states) and record the measurement results $\hat{p}_{ij}$ as the estimate of $p_{ij}=\text{Tr}(P_{i}\rho_{j})$. Assume that $M$ different types of probe states are employed and their total number of copies is $N$. Also assume different probe states use the same number of copies, which is $N/M$. We then aim to solve the following optimization problem:

Problem 1

Given experimental data $\{\hat{p}_{ij}\}$, solve $\min_{\{\hat{P}_{i}\}}\sum_{i=1}^{n}\sum_{j=1}^{M}[\hat{p}_{ij}-\text{Tr}(\hat{P}_{i}\rho_{j})]^{2}$ such that $\sum_{i=1}^{n}\hat{P}_{i}=I$ and $\hat{P}_{i}\geq 0$ for $1\leq i\leq n$.

III Estimation Algorithm

III-A Stage-1 approximation–constrained LRE

We first parameterize the detector and the input (probe) states. Let $\{\Omega_{i}\}_{i=1}^{d^{2}}$ be a complete set of $d$-dimensional traceless Hermitian matrices except $\Omega_{1}=I/\sqrt{d}$, and they satisfy $\text{Tr}(\Omega_{i}^{\dagger}\Omega_{j})=\delta_{ij}$, where $\delta_{ij}$ is the Kronecker function. Denote the detector by $\{P_{i}\}$ which are positive semidefinite and $\sum_{i=1}^{n}P_{i}=I$. Let $\rho_{j}$ be a series of input probe states. Then we can parameterize the detector and probe states as

[TABLE]

[TABLE]

where $\theta_{a}^{(i)}\triangleq\text{Tr}(P_{i}\Omega_{a})$ and $\phi_{b}^{(j)}\triangleq\text{Tr}(\rho_{j}\Omega_{b})$ are real. When $\rho_{j}$ is inputted, the probability to obtain the result corresponding to $P_{i}$ is calculated according to Born’s rule as

[TABLE]

Substituting (1) and (2) into (3), we obtain

[TABLE]

where $\Phi_{j}\triangleq(\phi_{1}^{(j)},\phi_{2}^{(j)},...,\phi_{d^{2}}^{(j)})^{T}$ and $\Theta_{i}\triangleq(\theta_{1}^{(i)},\theta_{2}^{(i)},...,\theta_{d^{2}}^{(i)})^{T}$. Suppose when estimating $\hat{p}_{ij}$, the outcome for $P_{i}$ appears $n_{ij}$ times, then $\hat{p}_{ij}=n_{ij}/(N/M)$. Denote the 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 $(p_{ij}-p_{ij}^{2})/(N/M)$. We thus have the linear regression equation

[TABLE]

Let $\Theta=(\Theta_{1}^{T},\Theta_{2}^{T},...,\Theta_{n}^{T})^{T}$, which is the vector of all the unknown parameters to be estimated. Collect the parametrization of the probe states as $X_{0}=(\Phi_{1},\Phi_{2},...,\Phi_{M})^{T}$. Let $\hat{Y}=(\hat{p}_{11},\hat{p}_{12},...,\hat{p}_{1M},\hat{p}_{21},\hat{p}_{22},...,\hat{p}_{2M},...,\hat{p}_{nM})^{T}$, $X=I_{n}\otimes X_{0}$, $e=(e_{11},e_{12},...,e_{1M},e_{21},e_{22},...,e_{2M},...,e_{nM})^{T}$, $H=(1,1,...,1)_{1\times n}\otimes I_{d^{2}}$, $D_{d^{2}\times 1}=(\sqrt{d},0,...,0)^{T}$. Then the regression equations can be rewritten in a compact form:

[TABLE]

with a linear constraint

[TABLE]

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

Problem 2

Given experimental data $\hat{Y}$, solve $\min_{\{\hat{P}_{i}\}}||\hat{Y}-X\hat{\Theta}||^{2}$ such that $H\hat{\Theta}=D$ and $\hat{P}_{i}\geq 0$ for $1\leq i\leq n$, where $\hat{\Theta}$ is the parametrization of ${\{\hat{P}_{i}\}}$ via (1).

Problem 2 is difficult to solve directly. Hence, we split it into two approximate subproblems:

Problem 2$.1$

Given experimental data $\hat{Y}$, solve $\min_{\{\hat{E}_{i}\}}||\hat{Y}-X\hat{\Theta}||^{2}$ such that $H\hat{\Theta}=D$, where $\hat{\Theta}$ is the parametrization of ${\{\hat{E}_{i}\}}$ via (1).

Problem 2$.2$

Given $\sum_{i=1}^{n}\hat{E}_{i}=I$, solve $\min_{\{\hat{P}_{i}\}}\sum_{i}||\hat{E}_{i}-\hat{P}_{i}||^{2}$ such that $\sum_{i}\hat{P}_{i}=I$ and $\hat{P}_{i}\geq 0$ for $1\leq i\leq n$.

Problem 2$.1$ is a linear regression problem with a linear constraint, and it can be solved analytically via the Constrained Least Squares (CLS) method [35]. Assume the input states have enough diversity such that $X^{T}X$ is nonsingular. This indicates $M\geq d^{2}$ for general complete probe-state sets. The standard CLS solution is [35]

[TABLE]

where $\hat{\Theta}_{LS}$ is unconstrained least square solution

[TABLE]

To further reduce the computational burden, we can simplify the form of (6) and (7). Let $Z_{0}=(X_{0}^{T}X_{0})^{-1}$. Then $(X^{T}X)^{-1}=I_{n}\otimes Z_{0}$, and

[TABLE]

Eq. (7) is in fact

[TABLE]

Also (6) is

[TABLE]

We then partition $\hat{Y}$ as $\hat{Y}^{T}=(\hat{Y}_{1}^{T},\hat{Y}_{2}^{T},...,\hat{Y}_{n}^{T})^{T}$ where $\hat{Y}_{i}=(\hat{p}_{i1},\hat{p}_{i2},...,\hat{p}_{iM})^{T}$ for $1\leq i\leq n$. Denote $Y_{0}=((1,...,1)_{1\times M})^{T}=\sum_{i}\hat{Y}_{i}$. We continue transforming (8) as

[TABLE]

Compared with (6) and (7), Eq. (9) is a faster way to calculate $\hat{\Theta}_{CLS}$.

Let $\hat{\Theta}_{CLS}=(\hat{\Theta}_{1}^{T},...,\hat{\Theta}_{n}^{T})^{T}$ and $\hat{\Theta}_{i}^{T}=(\hat{\theta}_{1}^{(i)},...,\hat{\theta}_{d^{2}}^{(i)})$. From $\hat{\Theta}_{CLS}$, we can obtain a stage-1 estimate $\hat{E}_{i}=\sum_{a=1}^{d^{2}}\hat{\theta}_{a}^{(i)}\Omega_{a}$. The error $||\hat{E}_{i}-E_{i}||$ between the estimate $\hat{E}_{i}$ and its true value $E_{i}$ will be referred to as the CLS error in the rest of this paper. Note that the positive semidefiniteness requirement on $\hat{E}_{i}$ is not considered at this stage, and $\hat{E}_{i}$ may have negative eigenvalues. Hence, we need to further adjust $\hat{E}_{i}$ to obtain a physical estimate.

III-B Difference decomposition

Now we begin to solve Problem 2$.2$. First we decompose each $\hat{E}_{i}$ as the difference of two positive semidefinite matrices $\hat{F}_{i}$ and $\hat{G}_{i}$: $\hat{E}_{i}=\hat{F}_{i}-\hat{G}_{i}$. There are infinitely many such decompositions, because a new decomposition will be obtained once another positive semidefinite matrix is added to $\hat{F}_{i}$ and $\hat{G}_{i}$. We hope to view $\hat{G}_{i}$ as small disturbance, and we thus seek a decomposition method to minimize the norm of $\hat{G}_{i}$.

For each $\hat{E}_{i}$, we perform a spectral decomposition to obtain $\hat{E}_{i}=\hat{W}_{i}\hat{K}_{i}\hat{W}_{i}^{\dagger}$, where $\hat{W}_{i}$ is unitary and $\hat{K}_{i}$ is real diagonal. We have

[TABLE]

Denote the optimal decomposition solution as $\hat{F}_{i}^{o}$ and $\hat{G}_{i}^{o}$. We assert that both $\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i}$ and $\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i}$ must be diagonal. Otherwise, we note that $\text{diag}(\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i})-\text{diag}(\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i})$ still equals to $\hat{K}_{i}$. Since $\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i}$ is positive semidefinite, all of its diagonal elements are thus nonnegative. This indicates that $\text{diag}(\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i})$ is also positive semidefinite. Similarly, $\text{diag}(\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i})$ is also positive semidefinite. Hence, $\text{diag}(\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i})$ and $\text{diag}(\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i})$ are also feasible solutions. Since $||\text{diag}(\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i})||<||\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i}||$, this contradicts the assumption that $\hat{G}_{i}^{o}$ is the optimal solution. Therefore, $\hat{W}_{i}^{\dagger}\hat{F}_{i}^{o}\hat{W}_{i}$ and $\hat{W}_{i}^{\dagger}\hat{G}_{i}^{o}\hat{W}_{i}$ must be diagonal. We then have

[TABLE]

and we can consider its elements: $(\hat{K}_{i})_{jj}=(\hat{W}_{i}^{\dagger}\hat{F}_{i}\hat{W}_{i})_{jj}-(\hat{W}_{i}^{\dagger}\hat{G}_{i}\hat{W}_{i})_{jj}$. If $(\hat{K}_{i})_{jj}>0$, we should take $(\hat{W}_{i}^{\dagger}\hat{F}_{i}\hat{W}_{i})_{jj}=(\hat{K}_{i})_{jj}$ and $(\hat{W}_{i}^{\dagger}\hat{G}_{i}\hat{W}_{i})_{jj}=0$; if $(\hat{K}_{i})_{jj}\leq 0$, we should take $(\hat{W}_{i}^{\dagger}\hat{F}_{i}\hat{W}_{i})_{jj}=0$ and $(\hat{W}_{i}^{\dagger}\hat{G}_{i}\hat{W}_{i})_{jj}=-(\hat{K}_{i})_{jj}$.

The optimal decomposition can be obtained through the following procedure. Assume there are $\hat{n}_{i}$ nonpositive eigenvalues for $\hat{E}_{i}$, and they are in decreasing order in $\text{diag}(\hat{K}_{i})$. Let

[TABLE]

and

[TABLE]

Then we know $\hat{F}_{i}\geq 0$, $\hat{G}_{i}\geq 0$ and $\hat{E}_{i}=\hat{F}_{i}-\hat{G}_{i}$, and this $\hat{G}_{i}$ has the least norm. Also note that for the true values we have $E_{i}=F_{i}=P_{i}$.

III-C Stage-2 approximation

From $I=\sum_{i}\hat{E}_{i}=\sum_{i}\hat{F}_{i}-\sum_{i}\hat{G}_{i}$, we have

[TABLE]

Since each $\hat{G}_{i}$ is positive semidefinite, we can decompose

[TABLE]

Then Eq. (11) is transformed into

[TABLE]

We let $\hat{A}_{i}=\hat{C}^{-1}\hat{F}_{i}(\hat{C}^{\dagger})^{-1}$, and then each $\hat{A}_{i}$ is positive semidefinite and their sum is the identity. Hence, $\{\hat{A}_{i}\}$ is a genuine estimate of the detector and we call $\{\hat{A}_{i}\}$ the stage-2 approximation. A further optimization is needed in order to obtain the final estimation result in the following.

III-D Unitary optimization

When decomposing $I+\sum_{i}\hat{G}_{i}=\hat{C}\hat{C}^{\dagger}$, there is in fact another degree of freedom. For any unitary $\hat{U}$, it holds that $\hat{C}\hat{C}^{\dagger}=\hat{C}\hat{U}\hat{U}^{\dagger}\hat{C}^{\dagger}$. Therefore, $\hat{U}^{\dagger}\hat{A}_{i}\hat{U}$ can also be an estimate of the detector. There might be multiple ways of defining the cost function to choose an appropriate $\hat{U}$, and here we illustrate one effective way. We hope to choose a $\hat{U}$ such that the effect of $\hat{C}$ is (partly) neutralized; i.e., $\hat{U}^{\dagger}\hat{C}^{-1}(\cdot)(\hat{C}^{\dagger})^{-1}\hat{U}$ has minimal effect on $\hat{F}_{i}$. This requires $\hat{C}\hat{U}$ to be close to $I$, and hence, we aim to minimize $||\hat{C}\hat{U}-I||$.

We have

[TABLE]

Let $L=-\text{Tr}(\hat{C}\hat{U}+\hat{U}^{\dagger}\hat{C}^{\dagger})+\text{Tr}[(\Lambda+\Lambda^{\dagger})(\hat{U}\hat{U}^{\dagger}-I)]$ where $\Lambda$ is a Lagrange multiplier matrix and the term $\text{Tr}[(\Lambda+\Lambda^{\dagger})(\hat{U}\hat{U}^{\dagger}-I)]$ amounts to the unitary requirement on $\hat{U}$. By partial differentiation we have

[TABLE]

Therefore,

[TABLE]

We perform a singular value decomposition to obtain $\hat{C}=\hat{U}_{\alpha}\hat{S}\hat{U}_{\beta}^{\dagger}$ where $\hat{U}_{\alpha}$ and $\hat{U}_{\beta}$ are unitary and $\hat{S}$ is diagonal and positive semidefinite. It is straightforward to verify that

[TABLE]

Let $\hat{U}_{\gamma}=\hat{U}_{\beta}^{\dagger}\hat{U}\hat{U}_{\alpha}$. Then (13) is now equivalent to

[TABLE]

Thus we have

[TABLE]

Therefore, $\hat{U}_{\gamma}\hat{S}^{2}\hat{U}_{\gamma}^{\dagger}$ is the spectral decomposition of $\hat{S}^{2}$. Since the probability for $\hat{S}$ to be degenerate is zero, we know $\hat{S}^{2}$ is nondegenerate. Thus we have $\hat{U}_{\gamma}=\text{diag}(\text{e}^{\text{i}\theta_{1}},\text{e}^{\text{i}\theta_{2}},...,\text{e}^{\text{i}\theta_{d}})$ where $\theta_{j}\in[0,2\pi)$ for $1\leq j\leq d$, which indicates that $\hat{U}_{\gamma}$ and $\hat{S}$ commutate. From (14), we then have $\hat{S}=\hat{U}_{\gamma}\hat{S}\hat{U}_{\gamma}=\hat{S}\hat{U}_{\gamma}^{2}$. When the resource number $N$ is large enough, $\hat{C}\hat{C}^{\dagger}$ will be close enough to the identity matrix, and $\hat{C}$ will be close to a unitary matrix (thus nonsingular) and we can view $\hat{S}$ as nonsingular. We thus have $\hat{U}_{\gamma}^{2}=I$, which indicates $\hat{U}_{\gamma}=\text{diag}(\pm 1,\pm 1,...,\pm 1)$. We further have $L=-2\text{Tr}(\hat{U}_{\gamma}\hat{S})$, which indicates $\hat{U}_{\gamma}=I$. Therefore, the optimal solution is

[TABLE]

Hence, the final estimation is $\hat{P}_{i}=\hat{U}^{\dagger}\hat{A}_{i}\hat{U}$ where $\hat{U}$ is determined through (15). The error $||\hat{P}_{i}-P_{i}||$ will be referred to as the final (estimation) error, in contrast to the CLS error $||\hat{E}_{i}-E_{i}||$.

III-E General procedure and computational complexity

We now generalize the procedure of our TSE algorithm and analyze its computational complexity. In this paper, we do not consider the time spent on experiments, since it depends on the experimental realization. In the following, we briefly summarize each step and illustrate their corresponding computational complexity.

Step 1. Stage-1 Approximation. Choose basis sets $\{\Omega_{i}\}$ and probe states $\rho_{j}$ and calculate $\Phi_{j}$. Then perform measurement experiments to collect data $\hat{p}_{ij}$. Obtain the constrained least square solution from (9) and construct the stage-1 approximation $\hat{E}_{i}=\sum_{a}\hat{\theta}_{a}^{(i)}\Omega_{a}$. In (9), both $(X_{0}^{T}X_{0})^{-1}X_{0}^{T}$ and $D$ can be calculated offline prior to the experiments, and the remaining online calculation has computational complexity $O(nd^{2}M)$.

Step 2. Difference Decomposition. Perform spectral decomposition on each $\hat{E}_{i}$ and obtain $\hat{E}_{i}=\hat{F}_{i}-\hat{G}_{i}$. Since the computational complexity of spectral decomposition is cubic in the dimension of a Hermitian matrix [36], this step has total computational complexity $O(nd^{3})$.

Step 3. Stage-2 Approximation. The transformation of $\sum_{i}\hat{F}_{i}$ into $\hat{C}\hat{C}^{\dagger}$ can be accomplished by spectral decomposition. Then, we obtain the stage-2 approximation $\hat{A}_{i}=\hat{C}^{-1}\hat{F}_{i}(\hat{C}^{\dagger})^{-1}$. The complexity of this step is $O(nd^{3})$.

Step 4. Unitary Optimization. Calculate the global unitary matrix $\hat{U}$ according to (15) and obtain the final estimate $\hat{P}_{i}=\hat{U}^{\dagger}\hat{A}_{i}\hat{U}$. This step has computational complexity $O(nd^{3})$.

Since $M\geq d^{2}$ for general complete probe-state sets, we have $nd^{3}\leq nd^{2}M$. Hence, our algorithm has total computational complexity $O(nd^{2}M)$.

IV Error Analysis

In this section, we present a theoretical upper bound for the final estimation error of our TSE algorithm. It is necessary to first characterize the probe states.

Assumption 1

The probe states used are optimal [13, 20]; i.e., they are $d$-dimensional pure states and $X_{0}^{T}X_{0}=c_{0}I$ where $c_{0}\in\mathbb{R}$. From [20], we have the following characterization:

[TABLE]

Let $\text{E}(\cdot)$ denote the expectation w.r.t. all possible measurement results. We present the following theorem to characterize the estimation error:

Theorem 1

Under Assumption 1, the final estimation error of our algorithm $\text{E}(\sum_{i}||\hat{P}_{i}-P_{i}||^{2})$ scales as $O(\frac{d^{5}n^{2}}{{N}})$, where $d$ is the system dimension, $n$ is the number of detector POVM elements and $N$ is the total number of resources.

Proof:

We prove the conclusion through analyzing the error in each step of our algorithm.

IV-A Error in stage-1 approximation

For simplicity, let $Z=(X^{T}X)^{-1}=I\otimes(X_{0}^{T}X_{0})^{-1}$. The estimation error for constrained LRE is

[TABLE]

From [13], we know $E(ee^{T})\leq\frac{M}{4N}I$. Therefore,

[TABLE]

We have

[TABLE]

It is clear that $(HZH^{T})^{-1}=X_{0}^{T}X_{0}/n$ and $HZ^{2}H^{T}=n(X_{0}^{T}X_{0})^{-2}$. Continuing (17), we have

[TABLE]

Hence, we have

[TABLE]

which we refer to as the CLS bound.

Remark 1

In cases when the last POVM element $P_{n}$ is omitted for simplicity, unconstrained LRE can be used for stage-1 approximation, and a corresponding error upper bound can be obtained as in [13]:

[TABLE]

Comparing (18) and (19), we find they are only different by a factor of $\frac{n-1}{n}$. For any given detector, $n$ is fixed and these two bounds behave the same, apart from a constant. We thus omit analysis for unconstrained LRE method.

IV-B Error in $||\hat{F}_{i}-F_{i}||$

We start from the spectral decomposition (10). Since $\hat{W}_{i}^{\dagger}E_{i}\hat{W}_{i}$ is positive semidefinite, its diagonal elements are all nonnegative. Therefore, we have

[TABLE]

IV-C Error in $||\hat{C}\hat{C}^{\dagger}-I||$

We have the following relationship:

[TABLE]

IV-D Error in $||\hat{C}\hat{U}-I||$

Let $\hat{S}^{2}=\text{diag}(1+s_{1},...,1+s_{d})$. We have

[TABLE]

where the last line comes from the fact $||\hat{C}\hat{C}^{\dagger}-I||=||\hat{U}_{\alpha}\hat{S}^{2}\hat{U}_{\alpha}^{\dagger}-I||=||\hat{S}^{2}-I||=\sqrt{\sum_{i}s_{i}^{2}}$.

Using (21), we know

[TABLE]

where we do not incorporate the constant into the $O$ notation before the end of this proof.

IV-E Error in $||(\hat{U}^{\dagger}\hat{C}^{\dagger})^{-1}-I||$

Denote the singular values of $\hat{C}\hat{U}$ as $\hat{\mu}_{i}$ for $1\leq i\leq d$. From (12) we know each $\hat{\mu}_{i}^{2}$ is an eigenvalue of $I+\sum_{i}\hat{G}_{i}$. Hence, we have $\hat{\mu}_{i}\geq 1$ for every $1\leq i\leq d$. Therefore,

[TABLE]

Using (22), we have

[TABLE]

where the last $\sim$ comes from omitting the higher order term.

IV-F Error in $\sum_{i}||\hat{P}_{i}-P_{i}||^{2}$

Since each $F_{i}=E_{i}$ is positive semidefinite, we have

[TABLE]

For each $i$, we have

[TABLE]

Using (20), (22), (23) and (24), we have

[TABLE]

where the first $\sim$ results from omitting the higher order terms, and the second $\sim$ comes from the relationship $||\hat{F}_{i}||\leq||\hat{F}_{i}-F_{i}||+||F_{i}||\sim||F_{i}||$, and the last inequality comes from the Cauchy-Schwarz inequality

[TABLE]

Taking the expectation of (25) and using (18), we have

[TABLE]

Since we have explicitly shown the constants in the $O$ notation, (26) should be interpreted as that the following equation holds asymptotically:

[TABLE]

Using (16), we can further simplify (26) as

[TABLE]

∎

Remark 2

If the probe states are not optimal, (28) might fail and only (27) holds. This proof also indicates that $\text{Tr}[(X_{0}^{T}X_{0})^{-1}]$ is a helpful index to guide the choice of the probe states. If different probe states are highly similar to each other, then they result in a large $\text{Tr}[(X_{0}^{T}X_{0})^{-1}]$ and thus a large estimation error.

V Optimization of the Coherent Probe States

Since the detector to be estimated is usually unknown in practice, the optimization among all the possible probe states should be independent of the specific detector. An advantage of our TSE method is that an explicit error upper bound is presented, which does not involve the specific form of the detector. This can be critical in the optimization of the probe states. Moreover, to adapt to practical applications, we assume the probe states are all coherent states in this section.

V-A On the types of probe states

In quantum optics experiments, the preparation of number states $|k\rangle$ ($k\in\mathbb{N}$) is a difficult task, especially when $k$ is large. Therefore, in practice the input probe states are usually coherent states instead. A coherent state is denoted as $|\alpha\rangle$ where $\alpha\in\mathbb{C}$ and it can be expanded using number states as

[TABLE]

Their inner product relationship is

[TABLE]

We usually identify $\alpha$ with $|\alpha\rangle$ when there is no ambiguity.

Let $|\alpha\rangle_{d}=\text{e}^{-\frac{|\alpha|^{2}}{2}}\sum_{i=0}^{d-1}\frac{\alpha^{i}}{\sqrt{i!}}|i\rangle$. Coherent states are in essence infinite dimensional. To estimate a $d$-dimensional detector, in practice we employ $|\alpha\rangle_{d}$ as the (approximate) mathematical description of $|\alpha\rangle$ in this paper. Throughout this paper we assume the detector gives no signal when saturated, which means the part $|\alpha\rangle-|\alpha\rangle_{d}$ can be distinguished from $|\alpha\rangle_{d}$.

The tail part $|\alpha\rangle-|\alpha\rangle_{d}$ can be viewed as noise, which should be suppressed. This requires the amplitude $|\alpha|$ to be not large. Furthermore, (29) indicates that if $|\alpha|$ and $|\beta|$ are both close to zero, their inner product will also be close to one, which means that coherent states with small amplitudes are very much “alike”. This indicates that we cannot employ probe state sets where all the amplitudes are small. Considering the above two requirements, we design the preparation procedure of the probe states as follows.

Probe States Preparation: Given appropriate $q>0$, generate two random numbers $x$ and $y$ independently with their probability density function uniformly distributed on $[-q,q]$. Then $|x+\text{i}y\rangle$ will be employed as a probe state, with $N/M$ copies. Repeat this process to generate $M$ probe states and employ them to perform detector tomography.

Remark 3

Our sampling procedure is in essence sampling randomly within a given square in the complex plane. Another candidate method is to sample following a certain symmetric fixed pattern within this given square. Since simulation shows little difference in the final estimation error, we stick to our random-sample procedure.

With our probe state preparation procedure, we wonder what is the relationship between $M$ and the final estimation error, when other factors, such as the detector, the total number of copies $N$ for the probe states and the parameter $q$, remain unchanged.

To ensure that the inversion of $X^{T}X$ in (7) exists, it is required that at least $M\geq d^{2}$. We further find that when $M$ is large enough, the final estimation error tends to a constant independent of $M$. We give an explanation as follows.

First, the $j$th probe state $|\alpha\rangle$ is approximately viewed as $|\alpha\rangle_{d}$, which has a corresponding parametrization $\Phi_{j}$. Let $\mathbb{E}(\cdot)$ denote the expectation of functions of $x$ and $y$, in contrast to the expectation $\text{E}(\cdot)$ in Theorem 1. Let $f_{j}=\Phi_{j}-\mathbb{E}(\Phi_{j})$. Then the $f_{j}$s are i.i.d. with respect to the subscript $j$. According to (18), the estimation error upper bound is $\frac{(n-1)M}{4N}\text{Tr}[(X_{0}^{T}X_{0})^{-1}]$. We thus have

[TABLE]

According to the central limit theorem [37], as $M$ tends to infinity, $\mathbb{E}(\frac{1}{M}\sum_{j=1}^{M}f_{j}f_{j}^{T})$ converges to a fixed matrix, and hence the expectation of estimation error tends to a constant.

Two points should be noted: (i) In practice $M$ cannot be arbitrarily large when $N$ is given. (ii) There is usually a gap between this bound (30) and the practical error. However, simulation results imply the effectiveness of the above analysis, which suggests that a modest number of different types of probe states should be enough for practical applications. To investigate the least $M$ that suffices for an estimation task with a given dimension, it only requires us to calculate $\mathbb{E}\{\frac{(n-1)M}{4N}\text{Tr}[(X_{0}^{T}X_{0})^{-1}]\}$ for several candidates of $M$, which is a quantity independent of the specific detector.

V-B Optimization of the size of sampling square for probe states

As analyzed in Sec. V-A, the estimation error would be large if $q$ is too small or too large. Hence, there should be an optimal value for the choice of $q$. This is further verified by the simulation results in Fig. 4.

To locate the optimal value of $q$, we consider the projection of a probe state onto the $d$-dimensional subspace where the detector resides. Theoretically, the optimal value of $q$ should be different for different detectors, even though the dimension is fixed. However, in simulations (for example, Fig. 4), we find that the optimal values for a practical detector and the bound $\mathbb{E}\{\frac{(n-1)M}{4N}\text{Tr}[(X_{0}^{T}X_{0})^{-1}]\}$ usually coincide. Therefore, as an approximation, we can investigate the optimization of this bound w.r.t. $q$. Furthermore, the value of $N$ does not affect this optimization, and from Sec. V-A we know an $M$ not too small will also be irrelevant to the optimization. Hence, we only need to optimize

[TABLE]

which is a quantity uniquely determined by the probe states.

We start from the real function defined on all nonnegative integers $g(k)=\text{e}^{-|\alpha|^{2}}\frac{|\alpha|^{2k}}{k!}$, where $\alpha$ is the corresponding complex number of a probe state $|\alpha\rangle$. From $g(k)-g(k+1)=\text{e}^{-|\alpha|^{2}}|\alpha|^{2k}\frac{k+1-|\alpha|^{2}}{(k+1)!}$, we know $g(k)$ first increases and then decreases after $k\geq|\alpha|^{2}-1$. Hence, $g(k)$ reaches the maximum value around $|\alpha|^{2}-1$.

For any given probe state $|\alpha\rangle$, we consider the amplitude of its projection on each position $|k\rangle\langle j|$, which is

[TABLE]

Fig. 1 shows the grided $h(k,j)$ with $d=8$ and $|\alpha|=2$. Note that $g(k)$ is the restriction of $h(k,j)$ on $j=k$. Using the same technique for analyzing $g(k)$, it is straightforward to prove that grided $h(k,j)$ always has a single peak, with the position of the maximum around $(|\alpha|^{2}-1,|\alpha|^{2}-1)$. Generally, the larger $h(k,j)$ is, the better accuracy one can expect to obtain for estimating the element of a detector at position $|k\rangle\langle j|$. To obtain the least estimation error, a natural idea is to maximize $h(k,j)$ for each position $(k,j)$. However, this is not practical, because from $\sum_{k=0}^{\infty}g(k)=1$ one can see that $\sum_{k,j}h(k,j)$ is bounded. Therefore, to locate the optimal $q$ means to optimally allocate $h(k,j)$ on the $d\times d$ positions.

Generally speaking, when estimating a multivariate target $\{\theta_{i}\}$, the MSE $E(\sum_{i}|\theta_{i}-\hat{\theta}_{i}|^{2})$ is usually dominated by the worst estimated parameter $\max_{i}|\theta_{i}-\hat{\theta}_{i}|$. Hence, the optimal $q$ (denoted as $q_{o}$) should have a good performance for the worst estimation. When $q$ is too small, $\mathbb{E}(|\alpha\rangle\langle\alpha|)$ is overly concentrated near the original point, and the projections on $(k,j)$s far from the original point will be too small, resulting in a large bound in (31); i.e., a bad estimate. Conversely, $q$ should not be too large. If we approximately view $h(k,j)$ as symmetric, it is natural to conclude that the projection of the maximum (or the middle point of the two maxima) of $h(k,j)$ should be at $(\frac{d-1}{2},\frac{d-1}{2})$ for $q_{o}$. More specifically, if $|\alpha|^{2}-1$ is an integer, then $g(|\alpha|^{2}-1)=g(|\alpha|^{2})$ are the two maxima. When $d$ is even, the maximum of $h(k,j)$ should be two contour points $(\frac{d-2}{2},\frac{d-2}{2})$ and $(\frac{d}{2},\frac{d}{2})$, and we should have $\frac{d-2}{2}=|\alpha|^{2}-1$. When $d$ is odd, $h(k,j)$ has one maximum and its projection should be $(\frac{d-1}{2},\frac{d-1}{2})$, which further indicates $\frac{d-1}{2}=\frac{|\alpha|^{2}+|\alpha|^{2}-1}{2}$. Therefore for $q_{o}$, we should always have

[TABLE]

From our sampling scheme for probe states in Sec. V-A, we have

[TABLE]

Combining (32) and (33), we have the following heuristic formula

[TABLE]

Remark 4

If the probe state is the tensor product of single-qubit probe states, then one only needs to optimize each single-qubit probe state, which corresponds to the 2-dimensional edition of (31). This can be straightforwardly achieved by running a numerical simulation, and the result is also covered in Fig. 5.

VI Numerical Results

VI-A Basic performance

We simulate the estimation error under different total resource numbers. We consider a $2$-dimensional system with a detector $P_{1}=\left(\begin{matrix}0&0\\ 0&0.3\\ \end{matrix}\right)$, $P_{2}=\left(\begin{matrix}0.1&-0.02i\\ 0.02i&0.2\\ \end{matrix}\right)$ and $P_{3}=\left(\begin{matrix}0.9&0.02i\\ -0.02i&0.5\\ \end{matrix}\right)$. The sampling parameter for coherent states is $q=0.015$. The number of different types of probe states is $M=40$. We employ our method to estimate the detector using different resource numbers and present the results in Fig. 2. In Fig. 2, the green dashed line is the theoretical CLS error upper bound (18), the black line is the theoretical final error upper bound (26) (or equivalently, (27) without the higher order term), and the blue dots and red diamonds are the CLS error and final error, respectively. The horizontal axis is the logarithm of the total number of copies of probe states $N$ and the vertical axis is the logarithm of the Mean Square Error (MSE) $\text{E}(\sum_{i}||\hat{P}_{i}-P_{i}||^{2})$. Each point in Fig. 2 is the average of 50 simulations.

In Fig. 2, the CLS bound is better than the final bound, which is because more relaxation procedures are used to deduce the final bound and make it looser. When the resource number is large ($N>10^{7}$), the decreasing slope is close to $-1$, which verifies Theorem 1. We also notice that when the resource number is small ($N<10^{6}$), the final estimation error is notably better than the CLS error. This is because the estimation error of arbitrary physical estimation is in essence bounded by a constant, while the CLS estimation can be nonphysical and thus leads to an unbounded error. As a result, when the resource number is not large enough, the CLS estimation is rough and the error exceeds this constant, while the final error is still bounded by this constant. This phenomenon disappears if $q$ is instead set close to the optimal value, because the final error will be too small to be influenced by the constant bound. For example, if $q=1.307$ as predicted by (34), the MSEs in Fig. 2 will decrease by $6$ orders of magnitude.

VI-B On the types of probe states

We simulate the performance of our algorithm with different number of types of probe states. The detector and $q$ are the same as in Sec. VI-A. The total resource number is $1.44\times 10^{9}$. We perform our estimation method with $M$ varying from $4$ to $4000$, and present the results in Fig. 3, where each point is the average of 100 simulations. The legend is the same as Fig. 2, except that the horizontal axis is the number of types of probe states $M$ in logarithm. We can see when $M$ is very small, both the theoretical bound and the practical errors are large, due to the fact that the probe states lack diversity and their linear dependence is high. When $M$ is over $10$, both the bound and the practical errors quickly tend to constants, which validates our analysis in Sec. V-A. Therefore, in practice a moderate number of different probe states should suffice.

VI-C Optimization of the size of sampling square for probe states

We perform simulations to illustrate that the optimal size of the sampling square for coherent probe states coincides with the optimal point of the bound (31). We consider a system with the same detector as that in Sec. VI-A. The total resource number is $N=10^{6}$, and the number of different types of probe states is $M=32$. We perform our estimation method under different $p$, and present the results in Fig. 4. Each point is the average of 200 simulations. We can see that there is indeed an optimal point for the practical estimation error with respect to different sizes $q$, which validates the analysis in Sec. V-B. Also this practical optimal position of $q$ basically coincides with the optimal position of the error bound.

Using the same system we simulate to search for the optimal size $q$ of the sampling square for probe states in different dimensions. The practical optimal positions we search for are the minimum points of the bound (31) under dimensions $d=2,4,8,16$, which are presented as red diamonds in Fig. 5. The blue line is the optimal $q_{o}$ predicted by our formula (34), which are close to the practical optimal values still with improvement space.

VI-D Comparison with MLE using qubit probes

We compare our TSE method with the Maximum Likelihood Estimation (MLE) method, which is one of the most widely used methods. We simulate an $N_{q}$-qubit detector with $P_{1}+P_{2}=I$ where

[TABLE]

and

[TABLE]

The probe states are the tensor product of single-qubit states $\{\frac{I}{2},\frac{I+\sigma_{x}}{2},\frac{I+\sigma_{y}}{2},\frac{I+\sigma_{z}}{2}\}$ where

[TABLE]

For each $N_{q}$, the total resource number of the probe states is $N=10^{3}\times 2^{3N_{q}}$, and they are evenly distributed to each probe state. The MLE algorithm we use is the method in [11]. We compare the estimation results of our TSE method and MLE in Fig. 6, where each point is the average of 10 simulations. The running time ($T$) is the online computational time. Note that the detector tomography method via MLE is in essence a numerical search algorithm and a theoretical characterization of the computational complexity is challenging. For each detector, we first run our algorithm, and then adjust the MLE method such that the averaged estimation error of MLE is within $[95\%,105\%]$ of the error of our algorithm. We see that for $N_{q}\geq 4$ qubits our algorithm can be faster than MLE by over 4 orders of magnitude. In this simulation, $M=d^{2}$, and we thus anticipate the computational complexity is $O(d^{4})$, which indicates a theoretical slope $1.204$ for our running time in the coordinate of Fig. 6. For the simulated running time of our algorithm, the slope of the fitting line of the right three points is $1.060$, which is close to the theoretical value but still with some difference, possibly because the qubit number is not large enough.

It is difficult to rigorously compare the computational complexity of the two algorithms in some averaging sense of all possible detectors. To give a simple illustration, we fix $N_{q}=3$ and $N=10^{3}\times 2^{9}$, and change the detector as $P_{1}^{\prime}+P_{2}^{\prime}=I$ with $P_{1}^{\prime}=kU_{R}P_{1}U_{R}^{\dagger}$, where $k$ is a random variable evenly distributed in $(0,1)$ and $U_{R}$ is a random unitary matrix generated following the algorithm in [38]. We independently generate 10 pairs of $k$ and $U_{R}$ and compare the averaged performance of our algorithm and MLE for these 10 pairs. For each pair, we still first run TSE method (with 10 repetitions) and then adjust MLE so that the averaged estimation error of MLE is within $[95\%,105\%]$ of the error of TSE. The final 10-pair-averaged error of TSE is $0.0581\pm 4.58\times 10^{-3}$, and $0.0584\pm 4.64\times 10^{-3}$ for MLE. The 10-pair-averaged running time of TSE is $1.50\times 10^{-3}\pm 8.14\times 10^{-5}$, and $9.39\pm 1.08$ for MLE, in seconds. This result generally matches the performance in Fig. 6.

We also simulate the case when $n$ increases. We fix $d=4$ and the detector is

[TABLE]

where for $j<n$ we have

[TABLE]

The probe states are the same as those in the above simulation. We choose $n$ to be a power of 2 and run the simulation for different values of $n$. The total resource number of the probe states is fixed as $N=10^{3}\times 2^{3}$, and they are evenly distributed to each probe state. We plot the running time ($T$) versus $n$ in logarithmic coordinates for our TSE method and MLE in Fig. 7, where each point is the average of 10 simulations. We see that TSE can be significantly faster than MLE for large $n$. Theoretically, $T=O(n)$ indicates a slope $0.301$ for our method. For the simulated running time of our algorithm, the slope of the fitting line of the right three points is $0.293$, which is close to the theoretical value. Furthermore, Fig. 6 and Fig. 7 also imply the relationship between the estimation error and $n$ and $d$, which is not close to the prediction of (28). One possible reason is that the bound (28) might not be tight. Also, note that the practical error is dependent on the specific detector and when $n$ and $d$ change the detector necessarily changes. Hence, we leave it an open problem to better characterize the increasing tendency of the error w.r.t. $n$ and $d$.

Remark 5

For practical detectors $n$ is usually smaller than $d$. However, this pattern means a very large $d$ in simulation, which is difficult to perform if we are to simulate the performance of MLE as comparison. Hence, we do not enforce large $d$ when performing simulation in Fig. 7.

VII Experimental Results

VII-A Experimental setup

We first briefly explain the entire experimental setup (as in Fig. 8(a)), which determines the structure of the detector to be estimated. More details about this setup can be found in [39].

In Fig. 8, the purple dashed box corresponds to the emulated quantum detector which works as two-mode inputs - one binary output detector. Two independent quantum modes are encoded within orthogonal polarization modes in one optical beam at the detector input. The two-mode quantum detector consists of two superconducting nanowire detectors (SNSPDs), a polarization beam splitter (PBS), a quarter wave plate (QWP), and a logical OR gate. The polarization of the input beam is first rotated by a QWP0 with the azimuth angle of 45∘ (Fig. 8(b)), or 30∘ (Fig. 8(c)), respectively. Then the beam is split into two spatially separated beams via PBS0, and they are injected into two SNSPDs through optical fibers. The photon counting signals from the two SNSPDs are sent to a logical OR gate, and the final detector output is obtained as on/off signal corresponding to POVMs of $P_{1}$ and $P_{0}$ ($P_{0}+P_{1}=I$). Fig. 8(b) and (c) are different specific settings to generate different emulated detectors.

This experimental setup leads to a special class of detectors. Specifically, we require them to be block diagonal (e.g., see [39]):

[TABLE]

where $m$ is the number of different blocks and $L^{(i)}_{j}\geq 0$ is $d_{j}\times d_{j}$ dimensional, with $\sum_{j=1}^{m}d_{j}=d$. Hence, we need to modify our original TSE method to reconstruct $\{P_{i}\}$.

VII-B Modified TSE protocol

First we choose $\{\Omega_{i}\}_{i=1}^{v}$ to be a complete orthogonal Hermitian basis set for the space of $\{P_{i}\}$ (instead of for $\mathbb{C}_{d\times d}$), where $\Omega_{1}=I_{d}/\sqrt{d}$ and $v$ equals to ${\sum_{j}d_{j}^{2}}$ instead of $d^{2}$. Then we have the parametrization under this basis set as

[TABLE]

[TABLE]

and the theoretical probability is $p_{ij}=\text{Tr}(P_{i}\rho_{j})$, which now becomes

[TABLE]

The linear regression equation is now

[TABLE]

and the error $e_{ij}=\hat{p}_{ij}-p_{ij}$ converges in distribution to a normal distribution with mean 0 and variance $(p_{ij}-p_{ij}^{2})/(N/M)$. Let $\Theta=(\Theta_{1}^{T},\Theta_{2}^{T},...,\Theta_{n}^{T})^{T}$ and $X_{0}=(\Phi_{1},\Phi_{2},...,\Phi_{M})^{T}$. Then $X_{0}$ is $M\times v$ dimensional. Let $Y=(\hat{p}_{11},\hat{p}_{12},...,\hat{p}_{1M},\hat{p}_{21},\hat{p}_{22},...,\hat{p}_{2M},...,\hat{p}_{nM})^{T}$, $X=I_{n}\otimes X_{0}$, $e=(e_{11},e_{12},...,e_{1M},e_{21},e_{22},...,e_{2M},...,e_{nM})^{T}$, $H=(1,1,...,1)_{1\times n}\otimes I_{v}$, $D_{v\times 1}=(\sqrt{d},0,...,0)^{T}$. Then the regression equations can be rewritten in a compact form:

[TABLE]

with a linear constraint

[TABLE]

which is the same form as (4) and (5), but with the dimensions of $\hat{\Theta}_{CLS}$ and $\hat{\Theta}_{LS}$ decreased from $nd^{2}$ to $nv$.

Before proceeding to the CLS solution, we introduce another amendment. In practical experiments, the types of the probe states are not always rich enough, and the resource number can be small. These limitations lead to large CLS error and thus unsatisfactory final errors. More specifically, physical estimations $\{\hat{P}_{i}\}$ always have the eigenvalues of $\hat{P}_{i}$ between 0 and 1, while bad nonphysical estimates usually make some of these eigenvalues far away from the region $[0,1]$, which indicates $||\hat{\Theta}_{CLS}||$ is too large. To avoid a CLS estimate that deviates seriously from the true value, we enforce a further requirement on the cost function of the linear regression process. Note that the original CLS problem is

[TABLE]

We now add an extra penalty item to modify (37) as

[TABLE]

where $\eta>0$. The new cost function is $\hat{Y}^{T}\hat{Y}-2\hat{Y}^{T}X\hat{\Theta}+\hat{\Theta}^{T}(X^{T}X+\eta I)\hat{\Theta}$. Hence, the new CLS solution is obtained by changing all the $X^{T}X$ items in (6) and (7) as $X^{T}X+\eta I$:

[TABLE]

and

[TABLE]

The modification from (37) to (38) is in essence Tikhonov regularization [40], and the optimal parameter $\eta$ is usually difficult to determine by a fixed formula. Note that as the total resource number of all the probe states $N$ increases, $||\hat{Y}-X\hat{\Theta}||$ usually decreases, and $\eta$ should also decrease. We thus choose $\eta=10^{3}/N$ for simplicity. From the CLS solution (40), we obtain the stage-1 estimate $\{\hat{P}_{i}\}$ which might not be positive semidefinite but satisfies all the other requirements.

The block diagonal structure of (36) implies that the detector is decoupled on the subspaces $\mathbb{C}_{d_{1}\times d_{1}},\mathbb{C}_{d_{2}\times d_{2}},...,\mathbb{C}_{d_{m}\times d_{m}}$. We thus can perform the procedures of Sec. III-B, III-C and III-D on these subspaces separately. Specifically, for each $1\leq j\leq m$, $\{\hat{L}^{(i)}_{j}\}_{i=1}^{n}$ is a set of Hermitian estimation on the space $\mathbb{C}_{d_{j}\times d_{j}}$ satisfying $\sum_{i=1}^{n}\hat{L}^{(i)}_{j}=I_{d_{j}}$. We thus employ difference decomposition, stage-2 approximation and unitary optimization in Sec. III-B- Sec. III-D on $\{\hat{L}^{(i)}_{j}\}_{i=1}^{n}$ to obtain a set of physical estimation $\{\hat{Q}^{(i)}_{j}\}_{i=1}^{n}$ for each $j$. The final estimation is thus $\hat{P}_{i}=\hat{Q}^{(i)}_{1}\oplus\hat{Q}^{(i)}_{2}\oplus...\oplus\hat{Q}^{(i)}_{m}$, which is physical and also satisfies the block-diagonal requirement.

Remark 6

An error upper bound similar to Theorem 1 can be given for this modified case. However, the upper bound requires that the form (37) without the penalty item is employed and also that $N$ should be large enough. In practical experiments, $N$ is difficult to be arbitrarily large due to noise and imperfections. Hence, we do not present the similar error bound in this paper.

VII-C Experimental results

We prepare two-mode coherent states for detector tomography by using an adequately attenuated continuous-wave (CW) fiber coupled laser as depicted in the yellow dashed box in Fig. 8(a). We express the general two-mode coherent state without global phase as $|\alpha,\beta\text{e}^{\text{i}\delta}\rangle$ ($\delta\in\mathbb{R}$, $\alpha,\beta\geq 0$), which can be expanded in the photon number basis as

[TABLE]

We can experimentally generate the above two-mode coherent states by attenuating the laser and rotating a QWP1 and a half wave plate (HWP1) after a PBS1. The probe states we used are the following 19 states:

[TABLE]

We performed experiments for two different sets of detectors, denoted as Group I and Group II, respectively. We take $\eta=10^{3}/N$ for both groups. For the true value of Group I (experimental setting as Fig. 8(b)), $P_{1}=L_{1}^{(1)}\oplus L_{2}^{(1)}\oplus L_{3}^{(1)}$, and we have $L_{1}^{(1)}=2.91\times 10^{-4}$,

[TABLE]

and

[TABLE]

For the true value of Group II (experimental setting as Fig. 8(c)), we have $L_{1}^{(1)}=1.27\times 10^{-4}$,

[TABLE]

and

[TABLE]

The error bars are at most 4%, which are derived from the precisions of quantum efficiency measurements for each SNSPD.

We record 100,000 measurement outcomes for each input state, and repeat it 6 times. By truncating the outcome records in the time axis we can obtain data for different resource numbers. We employ our modified algorithm to reconstruct the two sets of detectors, and show the results in Figs. 9 and 10, respectively. We also plot the reconstruction results using simulated measurement data as a comparison. In Fig. 9, the simulation matches the experiment very well. The performance in Fig. 10 is not as good as that for Group I, due to the influence of the nondiagonal elements with amplitudes significantly larger than zero.

VIII Conclusion

In this paper, we have proposed a novel Two-stage Estimation (TSE) quantum detector tomography method. We analysed the computational complexity for our algorithm and established an upper bound for the estimation error. We discussed the optimization of the coherent probe states, and presented simulation results to illustrate the performance of our algorithm. Quantum optical experiments were performed and the results validated the effectiveness of our method.

Acknowledgement

The authors would like to thank Dr. Akira Sone for helpful discussions.