\floatsetup

[table]capposition=top \floatsetup[figure]capposition=bottom \newfloatcommandcapbtabboxtable[][\FBwidth]

A two-stage solution to quantum process tomography: error analysis and optimal design

Shuixin Xiao, Yuanlong Wang, Jun Zhang, Daoyi Dong, Gary J. Mooney, Ian R. Petersen, and Hidehiro Yonezawa This research was supported by the Australian Research Council Future Fellowship Funding Scheme under Project FT220100656, Discovery Project Funding Scheme under Project DP210101938, the Centres of Excellence under Grant CE170100012, and the National Natural Science Foundation of China (62173229, 12288201). Shuixin Xiao is with School of Engineering and Technology, University of New South Wales, Canberra ACT 2600, Australia, and University of Michigan–Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai 200240, China, and also with School of Engineering, The Australian National University, Canberra, ACT 2601, Australia (e-mail: [email protected]). Yuanlong Wang is with Key Laboratory of Systems and Control, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China (e-mail: [email protected]).Jun Zhang is with University of Michigan–Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai 200240, China (e-mail: [email protected]).Daoyi Dong is with CIICADA Lab, School of Engineering, Australian National University, ACT 2601, Australia, and School of Physics, The University of Melbourne, Parkville, Victoria 3010, Australia (e-mail: [email protected]).Gary J. Mooney is with School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia (e-mail: [email protected]).Ian R. Petersen is with School of Engineering, Australian National University, ACT 2601, Australia (e-mail: [email protected]). Hidehiro Yonezawa is with RIKEN Center for Quantum Computing, Japan, and School of Engineering and Technology, University of New South Wales, Canberra ACT 2600, Australia, and also with Centre for Quantum Computation and Communication Technology, Australian Research Council, Canberra, ACT 2600, Australia (e-mail: [email protected]).

Abstract

Quantum process tomography is a critical task for characterizing the dynamics of quantum systems and achieving precise quantum control. In this paper, we propose a two-stage solution for both trace-preserving and non-trace-preserving quantum process tomography. Utilizing a tensor structure, our algorithm exhibits a computational complexity of $O(MLd^{2})$ where $d$ is the dimension of the quantum system and $M$ , $L$ represent the numbers of different input states and measurement operators, respectively. We establish an analytical error upper bound and then design the optimal input states and the optimal measurement operators, which are both based on minimizing the error upper bound and maximizing the robustness characterized by the condition number. Numerical examples and testing on IBM quantum devices are presented to demonstrate the performance and efficiency of our algorithm.

Index Terms:

Quantum process tomography, quantum system identification, error analysis, quantum system.

I Introduction

In the past decades, quantum science and technologies have significantly advanced in many fields such as quantum computation [1, 2], quantum communication [3, 4], quantum sensing [5] and quantum control [6, 7, 8]. To fully develop and realize these technologies, a fundamental challenge to overcome is to characterize unknown quantum dynamics. A typical framework to formulate this problem is quantum process tomography (QPT), where the parameters determining the map from input quantum states to output states need to be estimated.

For a closed quantum system, QPT is reduced to Hamiltonian identification and many results have been obtained in this scenario [9, 10, 11, 12, 13, 14]. For an open quantum system characterized by time-independent parameters, the process is completely positive (CP) and attention is usually restricted to the trace-preserving (TP) case. Using direct and convex optimization methods, Zorzi et al. [15, 16] gave the minimum positive operator-valued measure (POVM) resources required to estimate a CPTP quantum process, and Knee et al. [17] and Surawy-Stepney et al. [18] proposed iterative projection algorithms to identify such a process. The efficiency of parameter estimation for quantum processes was studied in [19]. Zhang and Sarovar [20] utilized the time traces of observable measurements to identify the parameters in a master equation. The identifiability and identification for the passive quantum systems were investigated in [21]. When the quantum process is non-trace-preserving (non-TP), Huang et al. [22] proposed a convex optimization method and Bongioanni et al. [23] demonstrated non-TP QPT based on a maximum likelihood approach in an experiment. In practice, a quantum process may be affected by non-Markovian noise. White et al. [24] developed a framework to characterize non-Markovian dynamics and the work in [25] proposed process tensor tomography for non-Markovian QPT. Moreover, adaptive strategies for QPT were proposed in [26, 27, 28].

In this paper, we focus on Standard Quantum Process Tomography (SQPT) [2, 14, 29, 30, 31] where states with the same dimension as the process are inputted, and the corresponding outputs are measured to reconstruct the process. We extend the framework in [2, 13] to allow for general input states and non-TP processes, and also directly connect the process with measurement data in a unified system equation. Using a mathematical representation, we deduce the special structural properties of the system equation, allowing for more efficient analysis, solution and optimal design. We give an analytical two-stage solution (TSS) to a general QPT problem, as a closed-form estimation which is not common among the existing QPT solutions. For $d$ -dimensional quantum systems, our algorithm has computational complexity $O(MLd^{2})$ and storage requirement $O(ML)$ , where $M\geq d^{2}$ is the number of different informationally complete or overcomplete input states and $L\geq d^{2}$ is the number of different measurement operators. We also establish an analytical error upper bound for our algorithm. Then we study the optimization of input states and measurement operators, based on this error upper bound and the condition number of the system equation, which reflects the robustness of our algorithm. We give theorems to characterize lower bounds on the error and condition number for both input states and measurement operators. We prove that SIC (symmetric informationally complete) states and MUB (mutually unbiased bases) states [32] are two examples achieving the lower bounds for our algorithm. We also prove MUB measurement achieves the lower bounds among measurement operators. Numerical examples and testing on IBM Quantum devices demonstrate the effectiveness of our algorithm and validate the theoretical error analysis. We compare our TSS algorithm with the convex optimization method in [14] and the results show that our algorithm is more efficient in both time and space costs during calculation. The main contributions of this paper are summarized as follows.

(i)

We propose an analytical two-stage solution (TSS) to a general QPT problem. Our TSS algorithm does not require specific prior knowledge about the unknown process and can be applied in both TP and non-TP quantum processes.
(ii)

Utilizing a tensor structure, our TSS algorithm has the computational complexity $O(MLd^{2})$ and a clear storage requirement $O(ML)$ , resulting from our closed-form estimation formula.
(iii)

Our TSS algorithm can give an analytical error upper bound which can be further utilized to optimize the input states and the measurement operators.
(iv)

Numerical examples and testing on IBM Quantum devices demonstrate the theoretical results and the effectiveness of our method. Compared to the convex optimization method in [14], our algorithm is more efficient in both time and space costs.

The organization of this paper is as follows. Section II introduces some preliminary knowledge and establishes our framework for QPT. Section III studies the tensor structure and presents our TSS algorithm for QPT. In Section IV, we analyze the computational complexity and the storage requirement. In Section V, we show the error analysis and in Section VI, we study the optimization of the input states and the measurement operators. Numerical examples are presented in Section VII and Section VIII concludes this paper.

Notation: The $i$ -th row and $j$ -th column of a matrix $X$ is $X_{ij}$ . The $j$ -th column of $X$ is $\operatorname{col}_{j}(X)$ . The transpose of $X$ is $X^{T}$ . The conjugate $(*)$ and transpose of $X$ is $X^{\dagger}$ . The sets of real and complex numbers are $\mathbb{R}$ and $\mathbb{C}$ , respectively. The sets of $d$ -dimension complex vectors and $d\times d$ complex matrices are $\mathbb{C}^{d}$ and $\mathbb{C}^{d\times d}$ , respectively. The identity matrix is $I$ . $\rm i=\sqrt{-1}$ . The trace of $X$ is $\text{Tr}(X)$ . The Frobenius norm of a matrix $X$ is denoted as $||X||$ and the 2-norm of a vector $x$ is $||x||$ . We use density matrix $\rho$ to represent a quantum state where $\rho=\rho^{\dagger}$ , $\rho\geq 0$ and $\operatorname{Tr}(\rho)=1$ and use a unit complex vector $|\psi\rangle$ to represent a pure state. The estimate of $X$ is $\hat{X}$ . The inner product of two matrices $X$ and $Y$ is defined as $\langle X,Y\rangle\triangleq\text{Tr}(X^{\dagger}Y)$ . The inner product of two vectors $x$ and $y$ is defined as $\langle x,y\rangle\triangleq x^{\dagger}y$ . The tensor product of $A$ and $B$ is denoted $A\otimes B$ . Hilbert space is $\mathbb{H}$ . The Kronecker delta function is $\delta$ . The $\operatorname{diag}(a)$ denotes a diagonal matrix with the $i$ -th diagonal element being the $i$ -th element of the vector $a$ . For any positive semidefinite $X_{d\times d}$ with spectral decomposition $X=UPU^{\dagger}$ , we define $\sqrt{X}$ or $X^{\frac{1}{2}}$ as $U\operatorname{diag}\left(\sqrt{P_{11}},\sqrt{P_{22}},\cdots,\sqrt{P_{dd}}\right)U^{\dagger}$ .

II Preliminaries and quantum process tomography

II-A Vectorization function

We introduce the (column-)vectorization function $\text{vec}:\mathbb{C}^{m\times n}\mapsto\mathbb{C}^{mn}$ . For a matrix $X_{m\times n}$ ,

	$\displaystyle\text{vec}(X_{m\times n})\triangleq$	$\displaystyle[X_{11},X_{21},\cdots,X_{m1},X_{12},\cdots,X_{m2},$		(1)
		$\displaystyle\cdots,X_{1n},\cdots,X_{mn}]^{T}.$		(1)

Thus, $\text{vec}(\cdot)$ is a linear and basis dependent map. We also define the inverse $\text{vec}^{-1}(\cdot)$ which maps a $d^{2}\times 1$ vector into a $d\times d$ square matrix. The common properties of $\text{vec}(\cdot)$ are listed as follows [33, 34]:

\langle X,Y\rangle=\langle\text{vec}(X),\text{vec}(Y)\rangle,

(2)

\text{vec}(XYZ)=(Z^{T}\otimes X)\text{vec}(Y),

(3)

\text{Tr}_{1}(\text{vec}(X)\text{vec}(Y)^{\dagger})=XY^{\dagger},

(4)

where $\text{Tr}_{1}(X)$ denotes the partial trace on the space $\mathbb{H}_{1}$ with $X$ belonging to the space $\mathbb{H}_{1}\otimes\mathbb{H}_{2}$ .

II-B Quantum process tomography

According to the system architecture, QPT generally can be divided into three classes: Standard Quantum Process Tomography (SQPT) [2, 14, 29, 30, 31], Ancilla-Assisted Process Tomography (AAPT) [35, 36, 37, 38] and Direct Characterization of Quantum Dynamics (DCQD) [39, 40, 41]. In SQPT, one inputs quantum states with the same dimension as the process, and then we reconstruct the process from quantum state tomography (QST) on the output states. In AAPT and DCQD, an auxiliary system (ancilla) is attached to the principal system, and the input states and the measurement of the outputs are both carried out on the extended Hilbert space. AAPT transforms QPT to QST in a larger Hilbert space and the input state can be separable [42]. DCQD requires entangled input states and measurements such that the measured probability distributions are directly related to the elements of the process representation [42]. Since AAPT needs high dimensional states in the extended Hilbert space and DCQD usually requires entanglement which are both more technologically demanding, we focus on the SQPT scenario in this paper.

Here, we briefly introduce the SQPT framework as given in [2, 13] and extend it to a more general framework based on the following two aspects:

(i)

Input Fock states are extended to be arbitrary states.
(ii)

The unknown process is extended from the TP case to both the TP and non-TP cases.

For a $d$ -dimensional quantum system, its dynamics can be described by a completely-positive (CP) linear map $\mathcal{E}$ . If we input a quantum state $\rho^{\text{in}}$ (input state), using the Kraus operator-sum representation [2], the output state $\rho^{\text{out}}$ is given by

\rho^{\text{out }}=\mathcal{E}\left(\rho^{\text{in }}\right)=\sum_{i=1}^{d^{2}}\mathcal{A}_{i}\rho^{\text{in }}\mathcal{A}_{i}^{\dagger},

(5)

where $\mathcal{A}_{i}\in\mathbb{C}^{d\times d}$ and satisfy

\sum_{i=1}^{d^{2}}\mathcal{A}_{i}^{\dagger}\mathcal{A}_{i}\leq I_{d}.

(6)

When the equality in (6) holds, the map $\mathcal{E}$ is said to be trace-preserving (TP). Otherwise, it is non-trace-preserving (non-TP). Let $\{E_{i}\}_{i=1}^{d^{2}}$ be a fixed basis of $\mathbb{C}^{d\times d}$ . From $\{\mathcal{A}_{i}\}_{i=1}^{d^{2}}$ , we can obtain process matrix $X$ [2, 13] which is a $d^{2}\times d^{2}$ Hermitian and positive semidefinite matrix such that

\mathcal{E}(\rho^{\text{in}})=\sum_{j,k=1}^{d^{2}}E_{j}\rho^{\text{in}}E_{k}^{\dagger}X_{jk}.

(7)

In this paper, we consider both TP and non-TP quantum processes and (6) becomes

\sum_{j,k=1}^{d^{2}}X_{jk}E_{k}^{\dagger}E_{j}\leq I_{d}.

(8)

Since the relationship between $X$ and $\mathcal{E}$ is one-to-one [2, 13], the identification of $\mathcal{E}$ is equivalent to identifying $X$ . The number of independent elements in $X$ is $d^{4}-d^{2}$ for the TP case and $d^{4}$ for the non-TP case.

Let $\{\sigma_{n}\}_{n=1}^{d^{2}}$ be a complete basis set of $\mathbb{C}^{d\times d}$ . Denote the set of all the input states as $\{\rho_{m}^{\text{in}}\}_{m=1}^{M}$ . We can expand each output state uniquely in $\{\sigma_{n}\}_{n=1}^{d^{2}}$ as

\rho_{m}^{\text{out}}=\mathcal{E}(\rho_{m}^{\text{in}})=\sum_{n=1}^{d^{2}}\alpha_{mn}\sigma_{n}.

(9)

We also define $\beta_{mn}^{jk}$ such that

E_{j}\rho_{m}^{\text{in}}E_{k}^{\dagger}=\sum_{n=1}^{d^{2}}\beta_{mn}^{jk}\sigma_{n}.

(10)

From the linear independence of $\{\sigma_{n}\}_{n=1}^{d^{2}}$ , the relationship between $X$ and $\alpha$ is [2, 13]

\sum_{j,k=1}^{d^{2}}\beta_{mn}^{jk}X_{jk}=\alpha_{mn}.

(11)

To guarantee that (11) has a unique solution, the maximal linear independent subset of $\{\rho_{m}^{\text{in}}\}_{m=1}^{M}$ must have $d^{2}$ elements. If this is achieved with $M=d^{2}$ , we call the input state set informationally complete. If this is achieved with $M>d^{2}$ , we call it informationally overcomplete. The work in [2, 13] is restricted to a special complete case, where $\{\rho_{m}^{\text{in}}\}_{m=1}^{M}=\{|j\rangle\langle k|\}_{j,k=1}^{d}$ (after a linear combination of the practical experiment results). Here we extend the framework of [2, 13] to the general informationally complete/overcomplete case with $M\geq d^{2}$ . We define the matrix $A$ where $A_{mn}\triangleq\alpha_{mn}$ and arrange the elements $\beta_{mn}^{jk}$ into a matrix $B$ :

B=\left[\begin{array}[]{*{7}{c}}\beta_{11}^{11}&\beta_{11}^{21}&\cdots&\beta_{11}^{12}&\beta_{11}^{22}&\cdots&\beta_{11}^{d^{2}d^{2}}\\ \beta_{21}^{11}&\beta_{21}^{21}&\cdots&\beta_{21}^{12}&\beta_{21}^{22}&\cdots&\beta_{21}^{d^{2}d^{2}}\\ \lx@intercol\dotfill\hfil\hfil\lx@intercol\\ \beta_{12}^{11}&\beta_{12}^{21}&\cdots&\beta_{12}^{12}&\beta_{12}^{22}&\cdots&\beta_{12}^{d^{2}d^{2}}\\ \beta_{22}^{11}&\beta_{22}^{21}&\cdots&\beta_{22}^{12}&\beta_{22}^{22}&\cdots&\beta_{22}^{d^{2}d^{2}}\\ \lx@intercol\dotfill\hfil\hfil\lx@intercol\\ \beta_{Md^{2}}^{11}&\beta_{Md^{2}}^{21}&\cdots&\beta_{Md^{2}}^{12}&\beta_{Md^{2}}^{22}&\cdots&\beta_{Md^{2}}^{d^{2}d^{2}}\\ \end{array}\right]

(12)

which is an $Md^{2}\times d^{4}$ matrix. We rewrite (11) into a compact form as

B\text{vec}(X)=\text{vec}(A),

(13)

where $B$ is determined once the bases $\{E_{i}\}_{i=1}^{d^{2}}$ and $\{\sigma_{n}\}_{n=1}^{d^{2}}$ are chosen.

In an experiment, assume that the number of copies for each input state is $N$ . Assume $P_{i,j}$ is the $i$ -th POVM element (i.e., measurement operator) in the $j$ -th POVM set where $1\leq j\leq J$ and $1\leq i\leq n_{j}$ . Therefore, the total number of different measurement operators is $L=\sum_{j=1}^{J}n_{j}$ and we have

P_{1,j}+P_{2,j}+\cdots+P_{n_{j},j}=I_{d},

(14)

for $1\leq j\leq J$ , which is the completeness constraint for each POVM set. We then vectorize these POVM elements as $\{P_{l}\}_{l=1}^{L}$ where, e.g., $P_{i+\sum_{q=1}^{j-1}n_{q}}\leftarrow P_{i,j}$ . These are the measurement operators on the output states. To extract all of the information in the output state, these measurement operators should be informationally complete or overcomplete, and thus $L\geq d^{2}$ . Some efficient methods for QST with fewer measurements are also discussed in [43, 44, 45]. Then, we can also uniquely expand $P_{l}$ in $\{\sigma_{n}\}_{n=1}^{d^{2}}$ as

P_{l}=\sum_{j=1}^{d^{2}}\mu_{lj}\sigma_{j},

(15)

and the ideal measurement probability $p_{ml}$ of the $m$ -th output state and the $l$ -th measurement operator is

$\displaystyle p_{ml}$	$\displaystyle=\operatorname{Tr}\left(\rho^{\text{out}}_{m}P_{l}\right)$
	$\displaystyle=\operatorname{Tr}\left(\sum_{n=1}^{d^{2}}\alpha_{mn}\sigma_{n}\sum_{j=1}^{d^{2}}\mu_{lj}\sigma_{j}\right)$
	$\displaystyle=\sum_{n,j=1}^{d^{2}}\alpha_{mn}\mu_{lj}\operatorname{Tr}\left(\sigma_{n}\sigma_{j}\right).$	(16)

Define the matrix $\mathcal{P}$ with $\mathcal{P}_{ml}\triangleq p_{ml}$ and let

C\!\triangleq\!\!\left[\begin{array}[]{cccc}\operatorname{Tr}\left(\sigma_{1}\sum_{j=1}^{d^{2}}\mu_{1j}\sigma_{j}\right)&\cdots&\operatorname{Tr}\left(\sigma_{d^{2}}\sum_{j=1}^{d^{2}}\mu_{1j}\sigma_{j}\right)\\ \operatorname{Tr}\left(\sigma_{1}\sum_{j=1}^{d^{2}}\mu_{2j}\sigma_{j}\right)&\cdots&\operatorname{Tr}\left(\sigma_{d^{2}}\sum_{j=1}^{d^{2}}\mu_{2j}\sigma_{j}\right)\\ \vdots&\vdots&\vdots\\ \operatorname{Tr}\left(\sigma_{1}\sum_{j=1}^{d^{2}}\mu_{Lj}\sigma_{j}\right)&\cdots&\operatorname{Tr}\left(\sigma_{d^{2}}\sum_{j=1}^{d^{2}}\mu_{Lj}\sigma_{j}\right)\end{array}\right],

(17)

which is determined by the measurement operators $\{P_{l}\}_{l=1}^{L}$ . Thus we have

\left(I_{M}\otimes C\right)\operatorname{vec}(A^{T})=\operatorname{vec}(\mathcal{P}^{T}).

(18)

Define $K$ such that $K\operatorname{vec}(A)=\operatorname{vec}(A^{T})$ and thus $K$ is a $Md^{2}\times Md^{2}$ commutation matrix. Using (13) and (18), we have

\left(I_{M}\otimes C\right)KB\operatorname{vec}(X)=\operatorname{vec}(\mathcal{P}^{T}),

(19)

which directly connects the unknown quantum process $X$ with measurement data $\mathcal{P}$ in the experiment.

Assume that the practical measurement result is $\hat{p}_{ml}$ and the measurement error is $e_{ml}=\hat{p}_{ml}-p_{ml}$ . According to the central limit theorem, $e_{ml}$ converges in distribution to a normal distribution with mean zero and variance $\left(p_{ml}-p_{ml}^{2}\right)/(N/J)$ [46, 47]. We consider the general scenario where there is no prior knowledge about whether the process is TP or non-TP, and the QPT can be formulated as the following optimization problem:

Problem 1

Given the parameterization matrix $B$ for the input states, the parameterization matrix $C$ for the measurement operators, the permutation matrix $K$ and the experimental data $\hat{\mathcal{P}}$ , find a Hermitian and positive semidefinite estimate $\hat{X}$ minimizing $\left\|\left(I_{M}\otimes C\right)KB\operatorname{vec}(\hat{X})-\operatorname{vec}(\widehat{\mathcal{P}}^{T})\right\|$ and satisfying the constraint (8).

III Tensor structure and two-stage quantum process tomography solution

In this section, we discuss how to simplify the structure of $B$ under a suitably chosen representation, which can reduce the computational complexity and storage requirement. Then we propose a two-stage solution (TSS) to obtain analytical tomography formulas for Problem 1. When there is prior knowledge that the process is TP, most of the following analysis still applies, unless otherwise specified.

III-A Tensor structure of QPT

Firstly, we consider the basis sets $\{E_{i}\}_{i=1}^{d^{2}}$ and $\{\sigma_{n}\}_{n=1}^{d^{2}}$ , which, if properly chosen, can simplify the QPT problem and greatly benefit the time and space costs of our algorithm. We choose these basis sets as the natural basis $\{|j\rangle\langle k|\}_{1\leq j,k\leq d}$ where $i=(j-1)d+k,n=(k-1)d+j$ . The advantages of this choice can be demonstrated as follows.

Lemma 1

[13, 15] If $\{E_{i}\}_{i=1}^{d^{2}}$ is chosen as the natural basis $\{|j\rangle\langle k|\}_{1\leq j,k\leq d}$ where $i=(j-1)d+k$ , then the constraint (8) becomes $\operatorname{Tr}_{1}(X)\leq I_{d}$ .

Lemma 1 fellows from the Choi-Jamiołkowski isomorphism and the proof for $\operatorname{Tr}_{1}(X)=I_{d}$ in the TP case can be found in [13, 15], which can be straightforwardly extended to non-TP cases. Define $F\triangleq\operatorname{Tr}_{1}(X)\leq I_{d}$ , and $F$ represents the success probability of the quantum process [23]. Furthermore, we have the following benefit.

Proposition 1

Define the collection of all the vectorized input states as

V\triangleq\left[\operatorname{vec}\left(\rho_{1}^{\text{in}}\right),\operatorname{vec}\left(\rho_{2}^{\text{in}}\right),\cdots,\operatorname{vec}\left(\rho_{M}^{\text{in}}\right)\right].

(20)

If $\{E_{i}\}_{i=1}^{d^{2}}$ and $\{\sigma_{n}\}_{n=1}^{d^{2}}$ are chosen as the natural basis $\{|j\rangle\langle k|\}_{1\leq j,k\leq d}$ where $i=(j-1)d+k,n=(k-1)d+j$ , then we have

\left(I_{d^{2}}\otimes V^{T}\right)R=B,

(21)

where $R$ is a $d^{4}\times d^{4}$ permutation matrix.

Proof:

Using (3), (10) becomes

\left(E_{k}^{*}\otimes E_{j}\right)\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)=\sum_{n=1}^{d^{2}}\beta_{{mn}}^{{jk}}\operatorname{vec}\left(\sigma_{n}\right).

(22)

Since $\{\sigma_{n}\}_{n=1}^{d^{2}}$ the is natural basis, using (20), we have

\left(E_{k}^{*}\otimes E_{j}\right)V=\left[\begin{array}[]{ccc}\beta_{1{1}}^{{jk}}&\cdots&\beta_{{M1}}^{{jk}}\\ \vdots&\vdots&\vdots\\ \beta_{1{d^{2}}}^{{jk}}&\cdots&\beta_{{Md^{2}}}^{{jk}}\end{array}\right],

(23)

and $\operatorname{vec}\left(\left[\left(E_{k}^{*}\otimes E_{j}\right)V\right]^{T}\right)$ is the $(j+(k-1)d^{2})$ -th column of $B$ . Let $k=(u-1)d+v$ , $j=(x-1)d+y$ . Then

E_{k}^{*}\otimes E_{j}=|u\rangle\langle v|\otimes|x\rangle\langle y|,

(24)

where only the element at $\left((u-1)d+x,(v-1)d+y\right)$ is non-zero. Then

\displaystyle\left(E_{k}^{*}\otimes E_{j}\right)\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)

\displaystyle=[0,\cdots,0,\underbrace{\left(\rho_{m}^{\text{in}}\right)_{yv}}_{((u-1)d+x)\text{th}},0,\cdots,0]^{T},

(25)

and therefore,

\left[\left(E_{k}^{*}\otimes E_{j}\right)V\right]^{T}=\left[\begin{array}[]{ccccccc}0&\cdots&0&\left(\rho_{1}^{\text{in}}\right)_{yv}&0&\cdots&0\\ \vdots&&\vdots&\vdots&\vdots&&\vdots\\ 0&\cdots&0&\left(\rho_{M}^{\text{in}}\right)_{yv}&0&\cdots&0\end{array}\right],

(26)

where only the $(u-1)d+x$ -th column is non-zero. Thus the $(j+(k-1)d^{2})$ -th column of $B$ is

\left[0,\cdots,0,\left(\rho_{1}^{\text{in}}\right)_{yv},\cdots,\left(\rho_{M}^{\text{in}}\right)_{yv},0,\cdots,0\right]^{T}.

(27)

The matrix $I_{d^{2}}\otimes V^{T}=\operatorname{diag}\left(\left[V^{T},\cdots,V^{T}\right]\right)$ is block diagonal. Thus the $(d^{2}\left[(u-1)d+x-1\right]+(v-1)d+y)$ -th column of $I_{d^{2}}\otimes V^{T}$ is the same as the $(j+(k-1)d^{2})$ -th column of $B$ . Therefore,

\left(I_{d^{2}}\otimes V^{T}\right)R=B,

(28)

where $R$ is a $d^{4}\times d^{4}$ permutation matrix. ∎

Since we choose $\{\sigma_{n}\}_{n=1}^{d^{2}}$ as the natural basis, we have

C=\left[\operatorname{vec}\left(P_{1}\right),\operatorname{vec}\left(P_{2}\right),\cdots,\operatorname{vec}\left(P_{L}\right)\right]^{T}.

(29)

Using Proposition 1 and (19), we have

\left(I_{M}\otimes C\right)K\left(I_{d^{2}}\otimes V^{T}\right)R\operatorname{vec}(X)=\operatorname{vec}({\mathcal{P}^{T}}).

(30)

Define $Y\triangleq\left(I_{M}\otimes C\right)K\left(I_{d^{2}}\otimes V^{T}\right)R$ . Since we assume the input states and measurement operators are both informationally complete or overcomplete, the unique least squares solution to $\min||Y\text{vec}(\hat{X})-\text{vec}(\widehat{\mathcal{P}}^{T})||$ is

\operatorname{vec}\left(\hat{X}_{\text{LS}}^{(1)}\right)=\left(Y^{\dagger}Y\right)^{-1}Y^{\dagger}\operatorname{vec}(\widehat{\mathcal{P}}^{T}).

(31)

However, the computational complexity for this least squares solution is usually $O(d^{12})$ [17] which is quite high. To reduce the computational complexity, we do not actually utilize (31). Noting the special structure property of (30), we firstly reconstruct $\hat{A}$ and then obtain $\hat{X}$ , which is more efficient in computation. To reconstruct $\hat{A}$ from (18), we have

\operatorname{vec}(\hat{A})=K^{T}\left(I_{M}\otimes\left(C^{\dagger}C\right)^{-1}C^{\dagger}\right)\operatorname{vec}(\widehat{\mathcal{P}}^{T}),

(32)

and then the least squares solution $\operatorname{vec}\left(\hat{X}_{\text{LS}}^{(2)}\right)$ is given by

	$\displaystyle\operatorname{vec}\left(\hat{X}_{\text{LS}}^{(2)}\right)$	$\displaystyle=\left(B^{\dagger}B\right)^{-1}B^{\dagger}\operatorname{vec}(\hat{A})$		(33)
		$\displaystyle=R^{T}\left(I_{d^{2}}\otimes\left(V^{}V^{T}\right)^{-1}V^{}\right)\operatorname{vec}(\hat{A}).$		(33)

Therefore, we can directly obtain the least squares solution $\operatorname{vec}\left(\hat{X}_{\text{LS}}^{(2)}\right)$ as

\displaystyle\operatorname{vec}\left(\hat{X}_{\text{LS}}^{(2)}\right)=

\displaystyle R^{T}\left(I_{d^{2}}\otimes\left(V^{*}V^{T}\right)^{-1}V^{*}\right)K^{T}\left(I_{M}\otimes\left(C^{\dagger}C\right)^{-1}C^{\dagger}\right)\operatorname{vec}(\widehat{\mathcal{P}}^{T}).

(34)

Remark 1

The procedure to obtain $\hat{A}$ in (32) is usually known as QST for output states [2, 13] but our method is different from such a procedure. Here we do not consider any constraints, i.e., Hermitian, positive semidefinite and unit trace, on the estimate $\hat{\rho}_{m}^{\text{out}}$ and we directly reconstruct the process from the data vector by linear inversion. Of course, it is possible to obtain an alternative physical estimate from $\hat{\rho}_{m}^{\text{out}}$ as in [48] or via Maximum Likelihood Estimation (MLE). However, here $\hat{A}$ is unnecessary to be a physical estimate because it is only an intermediate product, and our direct method can already guarantee the ultimate process estimation to be physical by later solving Problem 2 $.1$ and Problem 2 $.2$ in Sec. III-B.

Considering the positive semidefinite requirement $X\geq 0$ and the constraint $\operatorname{Tr}_{1}(X)\leq I_{d}$ , we convert QPT into an optimization problem as follows.

Problem 2

Given the parameterization matrix $V$ of all the input states, the permutation matrix $R$ and reconstructed $\hat{A}$ , find a Hermitian and positive semidefinite estimate $\hat{X}$ minimizing

\left\|\hat{X}-\operatorname{vec}^{-1}\left(R^{T}\left(I_{d^{2}}\otimes\left(V^{*}V^{T}\right)^{-1}V^{*}\right)\operatorname{vec}(\hat{A})\right)\right\|,

such that $\operatorname{Tr}_{1}(\hat{X})\leq I_{d}$ .

III-B Two-stage solution of QPT

The direct solution to Problem 2 is usually difficult and we split it into two optimization sub-problems (i.e., a two-stage solution):

Problem 2 $.1$

Let

\hat{D}\triangleq\hat{X}_{\text{LS}}^{(2)}=\operatorname{vec}^{-1}\left(R^{T}\left(I_{d^{2}}\otimes\left(V^{*}V^{T}\right)^{-1}V^{*}\right)\operatorname{vec}(\hat{A})\right)

be a given matrix. Find a Hermitian and positive semidefinite $d^{2}\times d^{2}$ matrix $\hat{G}$ minimizing $||\hat{G}-\hat{D}||$ .

Problem 2 $.2$

Let $\hat{G}\geq 0$ be given. Find a Hermitian and positive semidefinite $d^{2}\times d^{2}$ matrix $\hat{X}$ minimizing $||\hat{X}-\hat{G}||$ , such that $\operatorname{Tr}_{1}(\hat{X})\leq I_{d}$ .

For Problem 2 $.1$ , note that if $S=S^{\dagger},W=-W^{\dagger}$ , we have $\|S+W\|^{2}=\|S\|^{2}+\|W\|^{2}$ . Therefore,

\|\hat{G}-\hat{D}\|^{2}=\left\|\hat{G}-\frac{\hat{D}+\hat{D}^{\dagger}}{2}\right\|^{2}+\left\|\frac{\hat{D}-\hat{D}^{\dagger}}{2}\right\|^{2}.

(35)

We perform the spectral decomposition as $\frac{\hat{D}+\hat{D}^{\dagger}}{2}=U\hat{K}U^{\dagger}$ where $\hat{K}=\operatorname{diag}\left(k_{1},\cdots,k_{d^{2}}\right)$ is a diagonal matrix. We define $\hat{Z}\triangleq U^{\dagger}\hat{G}U$ . Since $\hat{G}$ is positive semidefinite, we have $\hat{Z}_{ii}\geq 0$ . Therefore,

$\displaystyle\left\\|\hat{G}-\frac{\hat{D}+\hat{D}^{\dagger}}{2}\right\\|^{2}$	$\displaystyle=\\|\hat{K}-\hat{Z}\\|^{2}$
	$\displaystyle=\sum_{i\neq j}\left\|\hat{Z}_{ij}\right\|^{2}+\sum_{i}\left(k_{i}-\hat{Z}_{ii}\right)^{2}$
	$\displaystyle\geq\sum_{k_{i}<0}\left(k_{i}-\hat{Z}_{ii}\right)^{2}$
	$\displaystyle\geq\sum_{k_{i}<0}k_{i}^{2}.$	(36)

Therefore, the unique optimal solution is $\hat{Z}=\operatorname{diag}(b)$ where

b_{i}=\begin{cases}k_{i},&k_{i}\geq 0,\\ 0,&k_{i}<0,\end{cases}

(37)

and $\hat{G}=U\operatorname{diag}(b)U^{\dagger}$ .

After solving Problem 2 $.1$ , we rewrite the spectral decomposition of $\hat{G}$ as

\displaystyle\hat{G}=\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger},

(38)

where each $\hat{S}_{i}\in\mathbb{C}^{d\times d}$ . Define $\hat{F}\triangleq\sum_{i=1}^{d^{2}}\hat{S}_{i}\hat{S}_{i}^{\dagger}$ . Using (4), we have

	$\displaystyle\operatorname{Tr}_{1}(\hat{G})$	$\displaystyle=\operatorname{Tr}_{1}\left(\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\right)$		(39)
		$\displaystyle=\sum_{i=1}^{d^{2}}\hat{S}_{i}\hat{S}_{i}^{\dagger}=\hat{F}.$		(39)

Then, we consider the partial trace constraint by correcting $\hat{G}$ . Assume that the spectral decomposition of $\operatorname{Tr}_{1}(X)=F$ is

F=U_{F}\operatorname{diag}\left(f_{1},\cdots,f_{d}\right)U_{F}^{\dagger},

(40)

where $1\geq f_{1}\geq\cdots f_{d}\geq 0$ and the spectral decomposition of $\hat{F}$ is

\hat{F}=U_{\hat{F}}\operatorname{diag}\left(\hat{f}_{1},\cdots,\hat{f}_{d}\right)U_{\hat{F}}^{\dagger},

(41)

where $\hat{f}_{1}\geq\cdots\geq\hat{f}_{c}>0$ and $\hat{f}_{c+1}=\cdots\hat{f}_{d}=0$ , i.e., $\operatorname{Rank}(\hat{F})=c$ . Then, we define

\bar{F}\triangleq U_{\hat{F}}\operatorname{diag}\left(\bar{f}_{1},\cdots,\bar{f}_{d}\right)U_{\hat{F}}^{\dagger},

(42)

where $\bar{f}_{i}=\hat{f}_{i}$ for $1\leq i\leq c$ , $\bar{f}_{i}=\frac{\hat{f}_{c}}{N}$ for $c+1\leq i\leq d$ , and $N$ is the number of copies for each input states. Since $\bar{F}$ is invertible, $\bar{F}^{-1/2}$ is well defined and we also define

\tilde{F}\triangleq U_{\hat{F}}\operatorname{diag}\left(\tilde{f}_{1},\cdots,\tilde{f}_{d}\right)U_{\hat{F}}^{\dagger},

(43)

where $\tilde{f}_{i}=\min\left({\bar{f}}_{i},1\right)$ for $1\leq i\leq d$ . Thus, $\tilde{f}_{i}\leq 1$ for $1\leq i\leq d$ . To solve Problem 2 $.2$ and satisfy the partial trace constraint, let $\hat{Q}_{i}\triangleq\tilde{F}^{1/2}\bar{F}^{-1/2}\hat{S}_{i}$ . We then obtain $\hat{X}$ as

\displaystyle\hat{X}=\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger},

(44)

where $\left\{\operatorname{vec}\left(\hat{Q}_{i}\right)\right\}_{i=1}^{d^{2}}$ might not be an orthogonal basis. From (44), $\hat{X}$ is Hermitian, positive semidefinite and also satisfies

$\displaystyle\operatorname{Tr}_{1}(\hat{X})$	$\displaystyle=\operatorname{Tr}_{1}\left(\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}\right)$
	$\displaystyle=\sum_{i=1}^{d^{2}}\hat{Q}_{i}\hat{Q}_{i}^{\dagger}=\tilde{F}^{1/2}\bar{F}^{-1/2}\sum_{i=1}^{d^{2}}\hat{S}_{i}\hat{S}_{i}^{\dagger}\bar{F}^{-1/2}\tilde{F}^{1/2}$
	$\displaystyle=\tilde{F}^{1/2}\bar{F}^{-1/2}\hat{F}\bar{F}^{-1/2}\tilde{F}^{1/2}$
	$\displaystyle=\tilde{F}^{1/2}U_{\hat{F}}\operatorname{diag}\left(1,\cdots,1,0,\cdots,0\right)U_{\hat{F}}^{\dagger}\tilde{F}^{1/2}$
	$\displaystyle=U_{\hat{F}}\operatorname{diag}\left(\tilde{f}_{1},\cdots,\tilde{f}_{c},0,\cdots,0\right)U_{\hat{F}}^{\dagger}$
	$\displaystyle\leq I_{d}.$	(45)

The Stage 1 solution $\hat{G}$ is a projection under the Frobenius norm satisfying the CP constraint, which was also discussed in [17, 18]. The TP constraint makes QPT a nontrivial extension of state tomography [17]. The work in [17, 18] introduced two iterative algorithms designed to ensure that the estimate satisfies the CPTP constraint. While these algorithms may offer good accuracy, our two-stage algorithm can achieve a closed-form solution. We can also efficiently optimize the input states and the measurement operators using the TSS algorithm. It is worth highlighting that our algorithm is versatile and can be applied to non-TP cases. Additionally, considering Hamiltonian identification in [13] where $\operatorname{Rank}(\hat{X})=1$ and $\operatorname{Tr}_{1}(\hat{X})=I_{d}$ , our Stage 2 solution is optimal for Problem 2 $.2$ which has been proved in [13]. For general QPT, the estimated process matrix $\hat{X}$ is a sub-optimal solution.

Remark 2

If we have the prior knowledge that the quantum process is TP, we can assume that $\hat{F}$ is non-singular (i.e., $\hat{F}>0$ ) and it is true when the number of copies is sufficiently large, because $\hat{F}$ converges to $I_{d}$ as $N$ tends to infinity. Thus, we correct $\hat{Q}_{i}$ as $\hat{Q}_{i}=\hat{F}^{-1/2}\hat{S}_{i}$ and $\hat{X}=\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}$ where $\operatorname{Tr}_{1}(\hat{X})=I_{d}$ .

IV Computational complexity and storage requirements

Refer to caption — Figure 1: Procedure for our TSS algorithm for QPT which has four steps. In Step 1, using (32), we reconstruct the parameterization matrix of all of the output states $\hat{A}$ directly from measurement data which is different from the QST in [2, 13]. Here $\hat{A}$ is not needed be a physical estimate because it is only an intermediate product. In Step 2, we utilize (33). Step 3 and Step 4 address the positive semidefinite constraint and partial trace constraint, respectively, by solving Problems 2 $.1$ and 2 $.2$ . These two steps together constitute our two-stage solution.

We summarize our procedure for the QPT framework with the TSS algorithm in Fig. 1. In this paper, we do not consider the time spent on experiments. Here we discuss the computational complexity and storage requirement in each step in the red box of Fig. 1 and give the total time and space complexity in Table I.

Step 1. We reconstruct $\hat{A}$ using (32). The computational complexity is $O(Ld^{4})$ for $\left(C^{\dagger}C\right)^{-1}C^{\dagger}$ , $O(MLd^{2})$ for $\left(I_{M}\otimes\left(C^{\dagger}C\right)^{-1}C^{\dagger}\right)\operatorname{vec}(\widehat{\mathcal{P}}^{T})$ and $O(Md^{2})$ for $K^{T}\left(I_{M}\otimes\left(C^{\dagger}C\right)^{-1}C^{\dagger}\right)\operatorname{vec}(\widehat{\mathcal{P}}^{T})$ where $M\geq d^{2}$ is the number of different input states and $L\geq d^{2}$ is the number of different measurement operators. Thus, the total computational complexity is $O(MLd^{2})$ . For storage requirement, the storage is $O(Ld^{2})$ for the parameterization matrix $C$ of measurement operators, $O(ML)$ for measurement data $\widehat{\mathcal{P}}$ and $O(Md^{2})$ for $\hat{A}$ . Since $K$ is a permutation matrix, it can be stored with an expense $O(Md^{2})$ .

Step 2. Using (33), the computational complexity is $O(Md^{4})$ for $\left(V^{*}V^{T}\right)^{-1}V^{*}$ , $O(Md^{4})$ for
$\left(I_{d^{2}}\otimes\left(V^{*}V^{T}\right)^{-1}V^{*}\right)\operatorname{vec}(\hat{A})$ and $O(d^{4})$ for $R^{T}\left(I_{d^{2}}\otimes\left(V^{*}V^{T}\right)^{-1}V^{*}\right)\operatorname{vec}(\hat{A})$ . Thus, the total computational complexity is $O(Md^{4})$ . We need to store $\left(V^{*}V^{T}\right)^{-1}V^{*}$ and $R^{T}$ . Since $\left(V^{*}V^{T}\right)^{-1}V^{*}$ is a $d^{2}\times M$ matrix and $R$ is a permutation matrix, the storage requirements are $O(Md^{2})$ and $O(d^{2})$ , respectively.

Step 3. To solve Problem 2 $.1$ , the computational complexity is determined by the spectral decomposition of $\frac{\hat{D}+\hat{D}^{\dagger}}{2}$ , which is $O(d^{6})$ . The storage requirements are $O(d^{4})$ for $\hat{D}$ , $O(d^{4})$ for the spectral decomposition of $\frac{\hat{D}+\hat{D}^{\dagger}}{2}$ , and $O(d^{4})$ for $\hat{G}$ .

Step 4. To solve Problem 2 $.2$ , the computational complexity for calculating $\hat{F}$ , $\bar{F}^{-1/2}$ and $\tilde{F}^{1/2}$ are $O(d^{4})$ , $O(d^{3})$ and $O(d^{3})$ , respectively. Then the computational complexity for calculating $\hat{X}$ is $O(d^{6})$ . The storage requirements for $\{\hat{S}_{i}\}_{i=1}^{d^{2}}$ and $\{\hat{Q}_{i}\}_{i=1}^{d^{2}}$ are both $O(d^{4})$ . Then the storage for $\hat{F}$ , $\bar{F}$ and $\tilde{F}$ are all $O(d^{2})$ , and for $\hat{X}$ is $O(d^{4})$ .

TABLE I: Computational complexity and storage requirement

	Computational complexity	Storage requirement
Step 1	$O(MLd^{2})$	$O(ML)$
Step 2	$O(Md^{4})$	$O(Md^{2})$
Step 3	$O(d^{6})$	$O(d^{4})$
Step 4	$O(d^{6})$	$O(d^{4})$
Total	$O(MLd^{2})$	$O(ML)$

The total computational complexity and storage requirements for our TSS algorithm are presented in Table I. Since $M\geq d^{2}$ and $L\geq d^{2}$ , the total computational complexity is $O(MLd^{2})$ . From [2, 42, 17], the computational complexity for a direct least squares solution with $d^{2}$ different input states is $O(d^{12})$ , which is significantly higher than $O(MLd^{2})$ . The reason we obtain a lower computational complexity here is that we utilize the special structure of $B$ described in Proposition 1. As for the total storage requirement, our TSS algorithm is $O(ML)$ . Without employing the tensor structure as (33), the storage space will be $O(Md^{6})$ because the storage requirement of $B$ will be $O(Md^{6})$ .

V Error analysis

From Fig. 1, our TSS algorithm has four steps. In this section, we analyze the error upper bound and give the following theorem to characterize the error upper bound analytically.

Theorem 1

If $\{E_{i}\}_{i=1}^{d^{2}}$ and $\{\sigma_{n}\}_{n=1}^{d^{2}}$ are chosen as the natural basis and the input quantum states are $\{\rho_{m}^{\text{in}}\}_{m=1}^{M}$ , then the estimation error of the TSS algorithm $\mathbb{E}||\hat{X}-X||$ scales as

O\!\!\left(\!\!\frac{\sqrt{d}\operatorname{Tr}({F})\sqrt{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)}\!\sqrt{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)}}{\sqrt{N}}\!\right),

where $N$ is the number of copies for each output state, $J$ is the number of POVM sets, $C$ is defined as (29), $V$ is defined as (20), $F=\operatorname{Tr}_{1}(X)$ and $\mathbb{E}(\cdot)$ denotes expectation with respect to all possible measurement results.

Proof:

We will first calculate the error for each of the four steps as shown in Fig. 1 and then present the final error bound. The main tools in this proof include the error analysis of QST in [46] and some matrix inequalities in Appendix A. For any quantity $S$ , we denote its estimation error as $\Delta_{S}\triangleq||\hat{S}-S||$ .

V-A Error in Step 1

The MSE of the $m$ -th estimated output state, i.e., $\operatorname{col}_{m}(\hat{A}^{T})$ is asymptotically (see [46])

		$\displaystyle\mathbb{E}\left\\|\operatorname{col}_{m}(\hat{A}^{T})-\operatorname{col}_{m}(A^{T})\right\\|^{2}$		(46)
	$\displaystyle=$	$\displaystyle\frac{J}{N}\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}C^{\dagger}\mathrm{P}C\left(C^{\dagger}C\right)^{-1}\right),$		(46)

where $N$ is the number of copies for each input state and

\mathrm{P}\triangleq\operatorname{diag}\left(p_{m1}-p_{m1}^{2},\cdots,p_{mL}-p_{mL}^{2}\right).

Therefore, the MSE is bounded by [46]

\mathbb{E}\left\|\operatorname{col}_{m}(\hat{A}^{T})-\operatorname{col}_{m}(A^{T})\right\|^{2}\leq\frac{J}{4N}\operatorname{Tr}\left(\!\left(C^{\dagger}C\right)^{-1}\right).

(47)

Then, for $\hat{A}$ , we have

\begin{array}[]{rl}\mathbb{E}\Delta_{A}^{2}&=\sum_{m=1}^{M}\mathbb{E}\left\|\operatorname{col}_{m}(\hat{A}^{T})-\operatorname{col}_{m}(A^{T})\right\|^{2}\\ &\leq\dfrac{MJ}{4N}\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right).\end{array}

(48)

V-B Error in Step 2

Since

	$\displaystyle\left\\|\left(V^{}V^{T}\right)^{-1}V^{}\right\\|^{2}$	$\displaystyle=\operatorname{Tr}\left(V^{T}\left(V^{}V^{T}\right)^{-1}\left(V^{}V^{T}\right)^{-1}V^{*}\right)$		(49)
		$\displaystyle=\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right),$		(49)

we have

	$\displaystyle\quad\Delta_{D}^{2}=\\|\hat{D}-D\\|^{2}$
	$\displaystyle=\left\\|\operatorname{vec}\left(\hat{D}\right)-\operatorname{vec}\left(D\right)\right\\|^{2}$
	$\displaystyle=\left\\|R^{T}\left(I\otimes\left(V^{}V^{T}\right)^{-1}V^{}\right)(\operatorname{vec}(\hat{A})-\operatorname{vec}(A))\right\\|^{2}$
	$\displaystyle\leq\left\\|\left(V^{}V^{T}\right)^{-1}V^{}\right\\|^{2}\\|\operatorname{vec}(\hat{A})-\operatorname{vec}(A)\\|^{2}$
	$\displaystyle=\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\Delta_{A}^{2},$		(50)

where we use Proposition 3 for the inequality.

V-C Error in Step 3

Since $\hat{G}$ minimizes $||\hat{G}-\hat{D}||$ , we have $\|\hat{G}-\hat{D}\|\leq\|\hat{D}-D\|$ and thus

\displaystyle\Delta_{G}=

\displaystyle\|\hat{G}-G\|=\|\hat{G}-D\|\leq\|\hat{G}-\hat{D}\|+\|\hat{D}-D\|\leq 2\Delta_{D},

(51)

because the true value $G=D$ .

V-D Error in Step 4

For $\hat{F}$ , if $1\geq\hat{f}_{1}\geq\cdots\geq\hat{f}_{d}\geq 0$ , i.e., $\hat{F}\leq I$ , we have $\tilde{F}^{1/2}\bar{F}^{-1/2}=I$ and thus $\hat{X}=\hat{G}$ . Therefore, using (48), (50) and (51), we obtain the final error bound as

\displaystyle\mathbb{E}\|\hat{X}-X\|=\mathbb{E}\|\hat{G}-G\|\sim O\left(\frac{\sqrt{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)}\sqrt{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)}}{\sqrt{N}}\right).

(52)

Then we consider the case where at least one eigenvalue of $\hat{F}$ is larger than $1$ . Assume $\hat{f}_{i}>1$ for $1\leq i\leq a$ and $0\leq\hat{f}_{i}\leq 1$ for $a+1\leq i\leq d$ in (41). From Lemma 2 and (39), we have

$\displaystyle\\|\hat{F}-F\\|$	$\displaystyle=\left\\|\operatorname{Tr}_{1}(\hat{G})-\operatorname{Tr}_{1}(G)\right\\|$	(53)
	$\displaystyle\leq\sqrt{d}\\|\hat{G}-G\\|$
	$\displaystyle=\sqrt{d}\Delta_{G}.$

Using (40), (41), and (131) in Appendix A, we can obtain

\max_{i}\left|\hat{f}_{i}-f_{i}\right|\leq\|\hat{F}-F\|\leq\sqrt{d}\Delta_{G},

(54)

and therefore we have

\operatorname{Tr}({F})-d\sqrt{d}\Delta_{G}\leq\operatorname{Tr}(\hat{F})\leq\operatorname{Tr}({F})+d\sqrt{d}\Delta_{G}.

(55)

Using (42) and (132) in Appendix A, we have

$\displaystyle\\|\tilde{F}-\hat{F}\\|^{2}$	$\displaystyle=\sum_{i=1}^{d}\left(\tilde{f}_{i}-\hat{f}_{i}\right)^{2}$
	$\displaystyle=\sum_{i=1}^{a}\left(1-\hat{f_{i}}\right)^{2}+\sum_{i=c+1}^{d}\left(\frac{\hat{f}_{c}}{N}\right)^{2}$
	$\displaystyle\leq\sum_{i=1}^{d}\left(f_{i}-\hat{f_{i}}\right)^{2}+(d-c)\left(\frac{\hat{f}_{c}}{N}\right)^{2}$
	$\displaystyle\leq\\|\hat{F}-{F}\\|^{2}+O\left(\frac{1}{N^{2}}\right).$	(56)

Since $\lim_{N\rightarrow\infty}\Delta_{G}=0$ , we have $\operatorname{Tr}(\hat{G})=\operatorname{Tr}(\hat{F})\sim\operatorname{Tr}({F})$ , and thus (56) indicates

\operatorname{Tr}(\hat{F})\sim\operatorname{Tr}(\tilde{F})\sim\operatorname{Tr}({F}).

(57)

Then we consider $\tilde{F}^{1/2}\bar{F}^{-1/2}$ as

\tilde{F}^{1/2}\bar{F}^{-1/2}=U_{\hat{F}}\operatorname{diag}\left(\sqrt{\frac{\tilde{f}_{1}}{\bar{f}_{1}}},\cdots\sqrt{\frac{\tilde{f}_{d}}{\bar{f}_{d}}}\right)U_{\hat{F}}^{\dagger}.

(58)

For $a+1\leq i\leq d$ , $\tilde{f}_{i}=\bar{f}_{i}$ . For $1\leq i\leq a$ , $\bar{f}_{i}=\hat{f}_{i}$ , and we have

\sqrt{\frac{\tilde{f}_{i}}{\bar{f}_{i}}}-1=\sqrt{\frac{1}{1+\left(\hat{f_{i}}-1\right)}}-1=\frac{1-\hat{f}_{i}}{2}+o\left(\frac{\hat{f_{i}}-1}{2}\right).

(59)

Thus, using (59) and (132) in Appendix A, we have

		$\displaystyle\quad\left\\|\tilde{F}^{1/2}\bar{F}^{-1/2}-I_{d}\right\\|^{2}=\sum_{i=1}^{a}\left(\frac{1-\hat{f}_{i}}{2}+o\left(\frac{1-\hat{f}_{i}}{2}\right)\right)^{2}$		(60)
		$\displaystyle\sim\frac{\sum_{i=1}^{a}\left(\hat{f}_{i}-1\right)^{2}}{4}\leq\frac{\sum_{i=1}^{d}\left(\hat{f_{i}}-f_{i}\right)^{2}}{4}\leq\frac{\\|\hat{F}-{F}\\|^{2}}{4}.$		(60)

Then using (45) and (57), we have

	$\displaystyle\operatorname{Tr}(\hat{X})$	$\displaystyle=\operatorname{Tr}\left(\operatorname{Tr}_{1}(\hat{X})\right)=\sum_{i=1}^{c}\tilde{f}_{i}$		(61)
		$\displaystyle=\operatorname{Tr}(\tilde{F})-\frac{(d-c)\hat{f}_{c}}{N}\sim\operatorname{Tr}({F}).$		(61)

Using (44), we thus have

	$\displaystyle\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\right\\|^{2}$	$\displaystyle=\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}\operatorname{vec}\left(\hat{Q}_{i}\right)$		(62)
		$\displaystyle=\operatorname{Tr}(\hat{X})\sim\operatorname{Tr}({F}).$		(62)

Similarly, using (38) and (39), we have

$\displaystyle\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$	$\displaystyle=\sum_{i=1}^{d^{2}}\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\operatorname{vec}\left(\hat{S}_{i}\right)$	(63)
	$\displaystyle=\operatorname{Tr}(\hat{G})=\operatorname{Tr}\left(\operatorname{Tr}_{1}(\hat{G})\right)$
	$\displaystyle=\operatorname{Tr}(\hat{F})\sim\operatorname{Tr}(F).$

Thus, from (3) and the definition of $\hat{Q}_{i}$ , we know $\operatorname{vec}\left(\hat{Q}_{i}\right)=\left(I_{d}\otimes\tilde{F}^{1/2}\bar{F}^{-1/2}\right)\operatorname{vec}\left(\hat{S}_{i}\right)$ , and we have

	$\displaystyle\quad\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle=\sum_{i=1}^{d^{2}}\left\\|\left(I_{d}\otimes\tilde{F}^{1/2}\bar{F}^{-1/2}-I_{d}\otimes I_{d}\right)\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle\leq\left\\|\tilde{F}^{1/2}\bar{F}^{-1/2}-I_{d}\right\\|^{2}\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle\sim\frac{\operatorname{Tr}({F})}{4}\\|\hat{F}-F\\|^{2}\sim\frac{d\operatorname{Tr}({F})}{4}\\|G-\hat{G}\\|^{2},$		(64)

where we use Proposition 3 in the inequality. Therefore, using (62)–(64), we have

	$\displaystyle\quad\\|\hat{X}-\hat{G}\\|^{2}$
	$\displaystyle=\left\\|\sum_{i=1}^{d^{2}}\left(\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\right)\right\\|^{2}$
	$\displaystyle\leq\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\right\\|\right)^{2}$
	$\displaystyle=\Bigg{(}\sum_{i=1}^{d^{2}}\Big{\\|}\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}+\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\Big{\\|}\Bigg{)}^{2}$
	$\displaystyle\leq\Bigg{(}\sum_{i=1}^{d^{2}}\left(\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}\right\\|+\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\right)\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\Bigg{)}^{2}$
	$\displaystyle\leq\left(\sum_{i=1}^{d^{2}}\left(\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\right\\|+\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\right)^{2}\right)\cdot\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)$
	$\displaystyle\leq 2\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\right\\|^{2}+\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)\cdot\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)$
	$\displaystyle\sim d\left(\operatorname{Tr}({F})\right)^{2}\\|G-\hat{G}\\|^{2},$		(65)

where we use the Cauchy-Schwarz inequality in the third inequality and AM-QM inequality [49] in the fourth inequality. Thus, using (48), (50), (51) and (65), the final error bound is

		$\displaystyle\quad\mathbb{E}\\|\hat{X}-X\\|\leq\mathbb{E}\\|\hat{X}-\hat{G}\\|+\mathbb{E}\\|\hat{G}-G\\|$		(66)
		$\displaystyle\leq(\sqrt{d}\operatorname{Tr}({F})+1)\mathbb{E}\\|\hat{G}-G\\|$
		$\displaystyle\leq O\left(\!\!\frac{\sqrt{d}\operatorname{Tr}({F})\sqrt{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)}\sqrt{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)}}{\sqrt{N}}\right).$

Comparing (52) and (66), we obtain the final error bound for any quantum process in the approximation sense as

\displaystyle\mathbb{E}\|\hat{X}-X\|\leq O\left(\frac{\sqrt{d}\operatorname{Tr}({F})\sqrt{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)}\sqrt{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)}}{\sqrt{N}}\right),

(67)

where $F=\operatorname{Tr}_{1}(X)$ represents the success probability of the quantum process [23]. ∎

Note that when the quantum process is TP, $\operatorname{Tr}(F)=d$ and the error analysis is similar to non-TP case. Then the final error bound is

\displaystyle\mathbb{E}\|\hat{X}-X\|\leq O\left(\frac{d^{3/2}\sqrt{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)}\sqrt{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)}}{\sqrt{N}}\right).

(68)

VI Optimization of input states and measurement operators

To minimize the tomography error, the optimal input states $\{\rho_{m}^{\text{in}}\}_{m=1}^{M}$ and optimal measurement operators $\{P_{l}\}_{l=1}^{L}$ are usually dependent on the specific process. One advantage of our TSS algorithm is that we can give an analytical error upper bound as (67) or (68) which depends on the input states and measurement operators instead of the unknown process. Therefore, we can optimize the input states and measurement operators based on the upper bound. In addition, we also consider how robust the estimated process $\hat{X}$ is w.r.t. the measurement error, which is determined by (32) and (33). Therefore, we use the condition numbers of $I_{M}\otimes C$ and $B$ to describe the robustness, which reflects the sensitivity of the solution to the perturbations in data. Moreover, we consider the relationship between $M$ and the final estimation error with randomly generated input states.

VI-A Optimal input quantum states

Define the set with $M$ different types of input states as

\displaystyle\mathcal{D}(d,M)\triangleq

\displaystyle\Big{\{}\{\rho_{m}^{\text{in}}\}_{m=1}^{M}\mid\forall 1\leq m\leq M,\rho_{m}^{\text{in}}\in\mathbb{C}^{d\times d},\rho_{m}^{\text{in}}=\left(\rho_{m}^{\text{in}}\right)^{\dagger},\rho_{m}^{\text{in}}\geq 0,\operatorname{Tr}(\rho_{m}^{\text{in}})=1\Big{\}}.

(69)

Similarly to the optimal probe states in quantum detector tomography [32, 50], we consider the optimal input states for QPT. The part of (67) or (68) dependent on the input states is $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ , which will be our first cost function. Therefore, we define the set of all the optimal input states to minimize this bound as

\mathrm{OIS}_{1}(d,M)\triangleq\underset{\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}\in\mathcal{D}(d,M)}{\arg\min}M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right).

(70)

For general QPT, we have no prior knowledge of the unknown process. From (33) with the reconstructed $\hat{A}$ , the following conditions are equivalent: (i) the process cannot be uniquely identified; (ii) the input states do not span the space of all $d$ -dimensional Hermitian matrices; (iii) $V^{*}V^{T}$ is singular; (iv) $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ is infinite. From (i) and (iv), it is thus reasonable to take $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ as a cost function.

In addition, maximizing the robustness of the reconstructed process with respect to measurement noise amounts to minimizing the condition number. Therefore, from (33), the second cost function is the condition number of $B$ ; i.e., $\operatorname{cond}(B)$ , defined as

\operatorname{cond}(B)\triangleq\frac{\sigma_{\max}(B)}{\sigma_{\min}(B)}

and $\sigma_{\max/\min}(B)$ is the maximum/minimum singular value of $B$ . Using (21), we have $\operatorname{cond}(B)=\operatorname{cond}(I)\operatorname{cond}(V^{T})=\operatorname{cond}(V)$ which also only depends on the input states. Similarly, we also define the set of all the optimal input states to minimize the condition number as

\mathrm{OIS}_{2}(d,M)\triangleq\underset{\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}\in\mathcal{D}(d,M)}{\arg\min}\operatorname{cond}(V).

(71)

We then define two sets $\underline{\mathrm{OIS}_{1}}(d,M)$ and $\underline{\mathrm{OIS}_{2}}(d,M)$ to characterize two lower bounds of these two cost functions as

\displaystyle\underline{\mathrm{OIS}_{1}}(d,M)\triangleq

\displaystyle\Big{\{}{\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}}\in\mathcal{D}(d,M)\mid M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)=d^{4}+d^{3}-d^{2}\Big{\}},

(72)

\underline{\mathrm{OIS}_{2}}(d,M)\triangleq\left\{{\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}}\in\mathcal{D}(d,M)\mid\operatorname{cond}(V)=\sqrt{d+1}\right\},

(73)

and the proof is as given in Theorem 2. In addition, we define the set of the optimal input states ( $\mathrm{OIS}(d,M)$ ) and the set of the optimal input states achieving both of the lower bounds ( $\underline{\mathrm{OIS}}(d,M)$ ) as

\displaystyle\mathrm{OIS}(d,M)\triangleq\mathrm{OIS}_{1}(d,M)\cap\mathrm{OIS}_{2}(d,M),

(74)

and

\displaystyle\underline{\mathrm{OIS}}(d,M)\triangleq\underline{\mathrm{OIS}_{1}}(d,M)\cap\underline{\mathrm{OIS}_{2}}(d,M).

(75)

Thus when $\mathrm{OIS}(d,M)$ is non-empty, it should minimize $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ and minimize $\operatorname{cond}(V)$ simultaneously. Additionally, when $\underline{\mathrm{OIS}}(d,M)$ is non-empty, it can achieve the lower bounds of $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ and $\operatorname{cond}(V)$ simultaneously. Then, we present the following theorem for $\underline{\mathrm{OIS}}(d,M)$ .

Theorem 2

For a $d$ -dimensional quantum process, with $M$ different types of input states and the parameterization matrix $V$ as in (20), we have $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}\geq d^{4}+d^{3}-d^{2}$ and $\operatorname{cond}(V)\geq\sqrt{d+1}$ . These lower bounds are achieved simultaneously; i.e., ${\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}}\in\underline{\mathrm{OIS}}(d,M)$ if and only if the eigenvalues of $V^{*}V^{T}$ are $\tau_{1}=\frac{M}{d}$ and $\tau_{2}=\dots=\tau_{d^{2}}=\frac{M}{d(d+1)}$ .

Proof:

Denote the eigenvalues of $V^{*}V^{T}$ as $\tau_{1}\geq\tau_{2}\geq\cdots\geq\tau_{d^{2}}>0$ . Since $\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}\in\mathcal{D}(d,M)$ , we can obtain two constraints $\sum_{i=1}^{d^{2}}\tau_{i}\leq M$ and $\tau_{1}\geq\frac{M}{d}$ coming from the purity requirement and unit trace of quantum states, respectively, which have been proved in [32]. To minimize $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ and $\operatorname{cond}(V)$ , we relax the optimization problems as follows:

	min	$\displaystyle\sum_{i=1}^{d^{2}}\frac{M}{\tau_{i}}$		(76)
	s.t.	$\displaystyle\sum_{i=1}^{d^{2}}\tau_{i}\leq M,\tau_{1}\geq\frac{M}{d},$		(76)

and

	min	$\displaystyle\sqrt{\frac{\tau_{1}}{\tau_{d^{2}}}}$		(77)
	s.t.	$\displaystyle\sum_{i=1}^{d^{2}}\tau_{i}\leq M,\tau_{1}\geq\frac{M}{d}.$		(77)

Using Lagrange multiplier method, we have

M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\geq d^{4}+d^{3}-d^{2},

(78)

and

\operatorname{cond}(V)\geq\sqrt{\frac{M/d}{M/\left(d(d+1)\right)}}=\sqrt{d+1}.

(79)

These lower bounds hold if and only if $\tau_{1}=\frac{M}{d}$ and $\tau_{2}=\dots=\tau_{d^{2}}=\frac{M}{d(d+1)}$ . ∎

Since $\mathcal{D}(d,M)$ is a closed space, $\text{OIS}_{1}(d,M)$ and $\text{OIS}_{2}(d,M)$ are always non-empty for any $M\geq d^{2}$ . However, it is not clear whether this also holds for $\text{OIS}(d,M)$ . For the lower bounds of these two cost functions, if one of them is achieved, we have $\tau_{1}=\frac{M}{d},\tau_{2}=\dots=\tau_{d^{2}}=\frac{M}{d(d+1)}$ and thus

\displaystyle\underline{\mathrm{OIS}}(d,M)=\underline{\mathrm{OIS}_{1}}(d,M)=\underline{\mathrm{OIS}_{2}}(d,M).

(80)

If the two lower bounds cannot be achieved, the above three sets are empty. Therefore, (80) always holds for any $M\geq d^{2}$ . However, we still have not fully determined what values of $M$ make (80) non-empty. In addition, when $\underline{\mathrm{OIS}}(d,M)$ is not empty, we have $\underline{\mathrm{OIS}}(d,M)={\mathrm{OIS}}(d,M)$ . If ${\left\{\rho_{m}^{\text{in}}\right\}_{m=1}^{M}}\in\underline{\mathrm{OIS}}(d,M)$ , all of them must be pure states from the purity requirement. This indicates that pure states may be better than mixed states in QPT for reducing the estimation error. We provide the following two examples in Appendix B which belong to $\underline{\mathrm{OIS}}(d,M)$ and achieve the lower bounds: SIC (symmetric informationally complete) states with the smallest $M=d^{2}$ and MUB (mutually unbiased bases) states for $M=d(d+1)$ . Moreover, we consider that the system is restricted to a multi-qubit system where $d=2^{m}$ and each input 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_{i}\leq M_{i}$ and there are $M_{i}$ different types of single qubits for the $i$ -th qubit of the product states. Thus, the total number of these $m$ -qubit tensor product states is thus $M=\prod_{i=1}^{m}{M_{i}}$ . We have the following theorem.

Theorem 3

For $m$ -qubit product input states, we have $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}\geq 20^{m}$ and $\operatorname{cond}(V)\geq\sqrt{3^{m}}$ . These lower bounds are achieved simultaneously if and only if for each $1\leq i\leq m$ , $\left\{\rho^{(i)}_{j_{i}}\right\}_{i=1}^{M_{i}}\in\underline{\mathrm{OIS}}(2,M)$ .

Proof:

Assuming that the parameterization matrix in (20) of all the $i$ -th single-qubit states $\left\{\rho_{j_{i}}^{(i)}\right\}_{j_{i}=1}^{M_{i}}$ is $V_{i}$ . We thus have $V=V_{1}\otimes V_{2}\cdots\otimes V_{m}$ . Therefore,

\displaystyle M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)

\displaystyle=M\prod_{i=1}^{m}\operatorname{Tr}\left(\left(V_{i}^{*}V_{i}^{T}\right)^{-1}\right)\geq M\prod_{i=1}^{m}\frac{20}{M_{i}}=20^{m},

(81)

and

\operatorname{cond}(V)=\prod_{i=1}^{m}\operatorname{cond}\left(V_{i}\right)\geq\prod_{i=1}^{m}\sqrt{3}=\sqrt{3^{m}}.

(82)

The above two equalities hold if and only if for each $1\leq i\leq m$ , $\left\{\rho^{(i)}_{j_{i}}\right\}_{i=1}^{M_{i}}\in\underline{\mathrm{OIS}}(2,M)$ . ∎

A similar proof for quantum detector tomography can also be found in [32]. For one-qubit input states, they are in the Bloch sphere and [32] has proved that platonic solids–tetrahedron, cube, octahedron, dodecahedron and icosahedron construct five examples, which belong to $\underline{\mathrm{OIS}}(2,M)$ with $M=4,6,8,12,20$ and achieve one-qubit lower bounds. $\underline{\mathrm{OIS}}(2,4)$ and $\underline{\mathrm{OIS}}(2,6)$ are one-qubit SIC states and one-qubit MUB states, respectively. For two-qubit product input states, Cube states (product states of two one-qubit MUB states) can achieve the lower bounds $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}=400$ and $\operatorname{cond}(V)=3$ .

VI-B Optimal measurement operators

For our error bound (67) or (68), the part depending on the measurement operators is $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ and we need to minimize it. Since $\operatorname{cond}(I_{M}\otimes C)=\operatorname{cond}(I_{M})\operatorname{cond}(C)=\operatorname{cond}(C)$ , we also consider the minimum of $\operatorname{cond}(C)$ .

Assume $P_{i,j}$ is the $i$ -th POVM element (i.e., a measurement operator) in the $j$ -th POVM set where $1\leq j\leq J$ and $1\leq i\leq n_{j}$ . Define the set size $\mathcal{n}\triangleq[n_{1},\cdots,n_{J}]$ . Therefore, the total number of measurement operators $L=\sum_{j=1}^{J}n_{j}$ . Let $\left\{\Omega_{i}\right\}_{i=1}^{d^{2}}$ be a complete set of $d$ -dimensional traceless Hermitian matrices except $\Omega_{1}=I/\sqrt{d}$ , and they satisfy $\operatorname{Tr}\left(\Omega_{i}^{\dagger}\Omega_{j}\right)=\delta_{ij}$ where $\delta_{ij}$ is the Kronecker function. Then we can parameterize the POVM element as

\displaystyle P_{i,j}

\displaystyle=\sum_{a=1}^{d^{2}}\phi_{i,j}^{a}\Omega_{a},

(83)

where $\phi_{i,j}^{a}\triangleq\operatorname{Tr}\left(P_{i,j}\Omega_{a}\right)$ is real. Thus we use $\phi_{i,j}\triangleq\left[\phi_{i,j}^{1},\cdots,\phi_{i,j}^{d^{2}}\right]^{T}$ as the parameterization of $P_{i,j}$ and the completeness constraint (14) becomes

\sum_{i=1}^{n_{j}}\phi_{i,j}=[\sqrt{d},0,\cdots,0]^{T},

(84)

and thus

\sum_{i=1}^{n_{j}}\phi_{i,j}^{1}=\sqrt{d},\quad\left\|\sum_{j=1}^{n_{j}}\phi_{i,j}\right\|=\sqrt{d}.

(85)

Conversely, using (85), we can also obtain (84). Define

\tilde{C}=\left[\phi_{1,1},\cdots\phi_{n_{1},1},\phi_{1,2},\cdots,\phi_{n_{J},J}\right]^{T},

(86)

which is the parameterization matrix of all the measurement operators in the basis of $\left\{\Omega_{i}\right\}_{i=1}^{d^{2}}$ and is a real matrix. Since $\left\{\Omega_{i}\right\}_{i=1}^{d^{2}}$ and natural basis are both orthonormal basis, the relationship between $C$ and $\tilde{C}$ is $\tilde{C}=CU$ where $U$ is a unitary matrix. Therefore, the eigenvalues of $C^{\dagger}{C}$ and $\tilde{C}^{T}\tilde{C}$ are the same, and $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)=J\operatorname{Tr}\left(\left(\tilde{C}^{T}\tilde{C}\right)^{-1}\right)$ , $\operatorname{cond}(\tilde{C})=\operatorname{cond}(C)$ .

Similarly, we define the set of all $d$ -dimensional measurement operators with set number $J$ and set size $\mathcal{n}$ as

\displaystyle\mathcal{M}(d,J,\mathcal{n})\triangleq\Big{\{}\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\mid\forall 1\leq j\leq J,\forall 1\leq i\leq n_{j},P_{i,j}\in\mathbb{C}^{d\times d},P_{i,j}=P_{i,j}^{\dagger},P_{i,j}\geq 0,\sum_{i=1}^{n_{j}}P_{i,j}=I_{d},\Big{\}}.

(87)

For general QPT, we have no prior knowledge of the output states. From (32), the following conditions are equivalent: (i) the output states cannot be uniquely identified; (ii) the measurement operators do not span the space of all $d$ -dimensional Hermitian matrices; (iii) $C^{\dagger}C$ is singular; (iv) $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ is infinite. From (i) and (iv), it is also reasonable to choose $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ as the first cost function. Thus we define the set of all the optimal measurement operators to minimize $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ ( $\mathrm{OMO}_{1}(d,J,\mathcal{n})$ ) as

\mathrm{OMO}_{1}(d,J,\mathcal{n})\triangleq\underset{\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\mathcal{M}(d,J,\mathcal{n})}{\arg\min}J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right).

(88)

Since maximizing the robustness of the reconstructed output states with respect to measurement noise is equivalent to minimizing the condition number, from (32), we also choose $\operatorname{cond}(C)$ as the second cost function. Thus, we define the set of all the optimal measurement operators to minimize $\operatorname{cond}(C)$ ( $\mathrm{OMO}_{2}(d,J,\mathcal{n})$ ) as

\mathrm{OMO}_{2}(d,J,\mathcal{n})\triangleq\underset{\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\mathcal{M}(d,J,\mathcal{n})}{\arg\min}\operatorname{cond}(C).

(89)

Then we define $s\triangleq\sum_{j=1}^{J}\frac{d}{n_{j}}$ and two sets $\underline{\mathrm{OMO}_{1}}(d,J,\mathcal{n})$ and $\underline{\mathrm{OMO}_{2}}(d,J,\mathcal{n})$ to characterize the lower bounds of these two cost functions as

\displaystyle\underline{\mathrm{OMO}_{1}}(d,J,\mathcal{n})\triangleq\Big{\{}\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\mathcal{M}(d,J,\mathcal{n})\mid J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)=J\left(\frac{1}{s}+\frac{\left(d^{2}-1\right)^{2}}{Jd-s}\right)\Big{\}},

(90)

\displaystyle\underline{\mathrm{OMO}_{2}}(d,J,\mathcal{n})\triangleq

\displaystyle\Big{\{}\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\mathcal{M}(d,J,\mathcal{n})\mid\operatorname{cond}({C})=\sqrt{\frac{\left(d^{2}-1\right)s}{Jd-s}}\Big{\}},

(91)

where the proof is given as in Theorem 4. We also define the set of all the optimal measurement operators ( $\mathrm{OMO}(d,J,\mathcal{n})$ ) as

\displaystyle\mathrm{OMO}(d,J,\mathcal{n})\triangleq\mathrm{OMO}_{1}(d,J,\mathcal{n})\cap\mathrm{OMO}_{2}(d,J,\mathcal{n}),

(92)

and the set of all the optimal measurement operators achieving both of the lower bounds ( $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ ) as

\displaystyle\underline{\mathrm{OMO}}(d,J,\mathcal{n})\triangleq\underline{\mathrm{OMO}_{1}}(d,J,\mathcal{n})\cap\underline{\mathrm{OMO}_{2}}(d,J,\mathcal{n}).

(93)

Thus, $\mathrm{OMO}(d,J,\mathcal{n})$ , when it is non-empty, should minimize $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ and minimize $\operatorname{cond}(C)$ simultaneously. When $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ is non-empty, it can achieve the lower bounds of $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ and $\operatorname{cond}(C)$ simultaneously. Then we present the following theorem for $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ .

Theorem 4

For a $d$ -dimensional quantum process, let the total number of POVM sets be $J$ and the number of POVM elements in the $j$ -th POVM set be $n_{j}$ . Then

J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)\geq J\left(\frac{1}{s}+\frac{\left(d^{2}-1\right)^{2}}{Jd-s}\right)

and

\operatorname{cond}(C)\geq\sqrt{\frac{\left(d^{2}-1\right)s}{Jd-s}}.

These lower bounds are achieved simultaneously; i.e., $\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ , if and only if the eigenvalues of $C^{\dagger}C$ are $v_{1}=s$ and $v_{2}=\dots=v_{d^{2}}=\frac{Jd-s}{d^{2}-1}$ .

Proof:

Denote the eigenvalues of $C^{\dagger}{C}$ or $\tilde{C}^{T}\tilde{C}$ as $v_{1}\geq v_{2}\geq\cdots\geq v_{d^{2}}>0$ . Since all the measurement operators are positive semidefinite, we have

\operatorname{Tr}\left((P_{m,j})^{\dagger}P_{n,j}\right)=\left(\phi_{m,j}\right)^{T}\phi_{n,j}\geq 0.

(94)

Using (84), we have

	$\displaystyle d$	$\displaystyle=\left\\|\sum_{i=1}^{n_{j}}\phi_{i,j}\right\\|^{2}=\sum_{i=1}^{n_{j}}\left\\|\phi_{i,j}\right\\|^{2}+\sum_{m\neq n}\left(\phi_{m,j}\right)^{T}\phi_{n,j}$		(95)
		$\displaystyle\geq\sum_{i=1}^{n_{j}}\left\\|\phi_{i,j}\right\\|^{2}.$		(95)

Therefore,

\sum_{i=1}^{d^{2}}v_{i}=\operatorname{Tr}\left(\tilde{C}^{T}\tilde{C}\right)=\sum_{j=1}^{J}\left(\sum_{i=1}^{n_{j}}\left\|\phi_{i,j}\right\|^{2}\right)\leq Jd.

(96)

Since

\sum_{i=1}^{n_{j}}\left(\phi_{i,j}^{1}\right)^{2}\geq\frac{\left(\sum_{i=1}^{n_{j}}\phi_{i,j}^{1}\right)^{2}}{n_{j}}=\frac{d}{n_{j}},

(97)

the first diagonal element of $\tilde{C}^{T}\tilde{C}$ is

\displaystyle\left(\tilde{C}^{T}\tilde{C}\right)_{11}

\displaystyle=\sum_{j=1}^{J}\left(\sum_{i=1}^{n_{j}}\left(\phi_{i,j}^{1}\right)^{2}\right)\geq\sum_{j=1}^{J}\frac{d}{n_{j}},

(98)

and thus

v_{1}\geq\left(\tilde{C}^{T}\tilde{C}\right)_{11}\geq\sum_{j=1}^{J}\frac{d}{n_{j}}.

(99)

To minimize $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)$ and $\operatorname{cond}(C)$ , the optimization problems can be relaxed as follows

	min	$\displaystyle\sum_{i=1}^{d^{2}}\frac{J}{v_{i}}$		(100)
	s.t.	$\displaystyle\sum_{i=1}^{d^{2}}v_{i}\leq Jd,v_{1}\geq\sum_{j=1}^{J}\frac{d}{n_{j}},$		(100)

and

	min	$\displaystyle\sqrt{\frac{v_{1}}{v_{d^{2}}}}$		(101)
	s.t.	$\displaystyle\sum_{i=1}^{d^{2}}v_{i}\leq Jd,v_{1}\geq\sum_{j=1}^{J}\frac{d}{n_{j}}.$		(101)

Using the Lagrange multiplier method, we have

\displaystyle J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)=J\operatorname{Tr}\left(\left(\tilde{C}^{T}\tilde{C}\right)^{-1}\right)\geq J\left(\frac{1}{s}+\frac{\left(d^{2}-1\right)^{2}}{Jd-s}\right),

(102)

and

\operatorname{cond}({C})\geq\sqrt{\frac{s}{(Jd-s)/(d^{2}-1)}}=\sqrt{\frac{\left(d^{2}-1\right)s}{Jd-s}}.

(103)

These lower bounds are achieved simultaneously if and only if $v_{1}=\sum_{j=1}^{J}\frac{d}{n_{j}}$ and $v_{2}=\dots=v_{d^{2}}=\frac{Jd-v_{1}}{d^{2}-1}$ . ∎

Since $\mathcal{M}(d,J,\mathcal{n})$ is a closed space, there always exist non-empty $\mathrm{OMO}_{1}(d,J,\mathcal{n})$ and $\mathrm{OMO}_{2}(d,J,\mathcal{n})$ for arbitrary set number $J$ and set size $\mathcal{n}$ . However, it is not clear whether this also holds for $\mathrm{OMO}(d,J,\mathcal{n})$ . For the lower bounds of these two cost functions, if one of them is achieved, we have $v_{1}=\sum_{j=1}^{J}\frac{d}{n_{j}},v_{2}=\dots=v_{d^{2}}=\frac{Jd-v_{1}}{d^{2}-1}$ and thus

\displaystyle\underline{\mathrm{OMO}}(d,J,\mathcal{n})=\underline{\mathrm{OMO}_{1}}(d,J,\mathcal{n})=\underline{\mathrm{OMO}_{2}}(d,J,\mathcal{n}).

(104)

If the two lower bounds can not be achieved, the above three sets are empty. Therefore, (104) always holds for any $J$ and $\mathcal{n}$ . Our results show that if $\left\{P_{i,j}\right\}_{i,j=1}^{n_{j},J}\in\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ , then in each POVM set, these POVM elements are orthogonal; i.e., $\operatorname{Tr}\left((P_{m,j})^{\dagger}P_{n,j}\right)=0$ $(m\neq n)$ from the first constraint. To achieve these lower bounds, we must have $J\geq d$ because the maximum number of orthogonal POVM elements in one set is $d$ . It is an open problem as to whether $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ is non-empty for arbitrary set number $J$ and set sizes $n_{j}$ ( $1\leq j\leq J$ ). When it is non-empty, we have $\underline{\mathrm{OMO}}(d,J,\mathcal{n})={\mathrm{OMO}}(d,J,\mathcal{n})$ .

Here, we give one example of $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ : MUB (mutually unbiased bases) measurement with $J=d+1$ and $n_{j}=d$ for $1\leq j\leq d+1$ , which can achieve these lower bounds. Two sets of orthogonal bases $\mathcal{H}^{m}=\{\left|\psi_{i}^{m}\right\rangle\}_{i=1}^{d}$ and $\mathcal{H}^{n}=\{\left|\psi_{j}^{n}\right\rangle\}_{j=1}^{d}$ are called mutually unbiased if and only if [51]

\left|\left\langle\psi_{i}^{m}|\psi_{j}^{n}\right\rangle\right|^{2}=\left\{\begin{array}[]{ll}\frac{1}{d}&\text{ for }m\neq n,\\ \delta_{ij}&\text{ for }m=n.\end{array}\right.

(105)

For a prime $p$ and a positive integer $k$ , there exist maximally $d+1$ sets of mutually unbiased bases in Hilbert spaces of prime-power dimension $d=p^{k}$ [52]. For other values, it is still an open problem as to whether there exist $d+1$ sets of mutually unbiased bases. The corresponding MUB measurements are $\left\{\left|\psi_{i}^{m}\right\rangle\left\langle\psi_{i}^{m}\right|\right\}_{m,i=1}^{d+1,d}$ . Then we give the following proposition with a similar proof to that of Proposition 2 in [32].

Proposition 2

$d$ -dimensional MUB measurements (when they exist) belong to $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ with $J=d+1$ and $n_{j}=d$ for $1\leq j\leq d+1$ .

Proof:

We assume that the standard singular value decomposition (SVD) of $C$ is $C=U_{c}\Sigma V^{\dagger}_{c}$ . Thus, using (29) and (105), we have

CC^{\dagger}=\left[\begin{array}[]{cccc}I_{d}&\left(\frac{1}{d}\right)_{d}&\cdots&\left(\frac{1}{d}\right)_{d}\\ \left(\frac{1}{d}\right)_{d}&I_{d}&\cdots&\left(\frac{1}{d}\right)_{d}\\ \vdots&\vdots&\vdots&\vdots\\ \left(\frac{1}{d}\right)_{d}&\left(\frac{1}{d}\right)_{d}&\cdots&I_{d}\end{array}\right]=U_{c}\Sigma\Sigma^{T}U^{\dagger}_{c}

(106)

where $\left(\frac{1}{d}\right)_{d}$ denotes the $d\times d$ matrix where all the elements are $\frac{1}{d}$ . Therefore,

\displaystyle C^{\dagger}C

\displaystyle=V_{c}\Sigma^{T}U^{T}_{c}U_{c}\Sigma V^{\dagger}_{c}=V_{c}\operatorname{diag}\left(d+1,I_{d^{2}-1}\right)V^{\dagger}_{c}.

(107)

Thus the eigenvalues of $C^{\dagger}C$ are $v_{1}=d+1,v_{2}=\cdots=v_{d^{2}}=1$ for MUB measurements and we can then calculate $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)=d^{3}+d^{2}-d$ and $\operatorname{cond}(C)=\sqrt{d+1}$ , which achieve the lower bounds. ∎

Remark 3

SIC-POVM are usually thought as optimal (in certain senses) measurements in quantum physics. The simplest mathematical definition of SIC-POVM is a set of $d^{2}$ normalized vectors $\left|\psi_{k}\right\rangle$ in $\mathbb{C}^{d}$ satisfying [53]

\left|\left\langle\psi_{j}|\psi_{k}\right\rangle\right|^{2}=\frac{1}{d+1},j\neq k.

(108)

The corresponding POVM elements are $P_{i}=\frac{1}{d}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$ for $1\leq i\leq d^{2}$ and $\sum_{i=1}^{d^{2}}P_{i}=I$ with $J=1$ , $n_{1}=d^{2}$ . However, here SIC-POVM does not belong to $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ because these POVM elements are not orthogonal with each other. Therefore, the first constraint is not satisfied. We can calculate $J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)=d^{3}+d^{2}-d$ and $\operatorname{cond}(C)=\sqrt{d+1}$ for SIC-POVM. These values are the same as MUB measurement which has $J=d+1$ and $n_{j}=d$ for $1\leq j\leq d+1$ . Thus we conjecture that SIC-POVM belongs to $\mathrm{OMO}(d,J,\mathcal{n})$ for $J=1$ , $n_{1}=d^{2}$ .

Remark 4

Our analysis for optimal measurement operators can also be applied to QST. A similar problem for QST has been discussed in [46]. However, their work is restricted to optimal projective measurement, while our result considers general POVM. Miranowicz et al. [54] also considered the minimum condition number for the optimal measurement operators in QST and proposed generalized Pauli operators whose condition number is $1$ . However, generalized Pauli operators are realized from projective measurements and their condition number does not directly link to the experimental data.

With input states like SIC states or MUB states belonging to $\underline{\mathrm{OIS}}(d,M)$ and measurement operators like MUB measurements belonging to $\underline{\mathrm{OMO}}(d,J,\mathcal{n})$ , the error upper bound of (67) becomes $O\left(\frac{d^{4}\operatorname{Tr}(F)}{\sqrt{N}}\right)$ . The limitation of our error upper bound is that it might be loose with respect to the system dimension $d$ and we leave it an open problem to characterize $d$ more accurately. We can obtain a tighter upper bound if we use the natural basis states $\{|j\rangle\langle k|\}_{1\leq j,k\leq d}$ and $M=d^{2}$ with MUB measurements. Note that in experiments, all of the quantum states are Hermitian. Hence, we cannot actually use $|j\rangle\langle k|$ ( $j\neq k$ ) as input quantum states. According to [2], since we have

\displaystyle|j\rangle\langle k|=

\displaystyle|+\rangle\langle+|+\mathrm{i}|-\rangle\langle-|-\frac{1+\mathrm{i}}{2}|j\rangle\langle j|-\frac{1+\mathrm{i}}{2}|k\rangle\langle k|,

(109)

where $|\pm\rangle=(|j\rangle\pm|k\rangle)/\sqrt{2}$ , $\mathcal{E}(|j\rangle\langle k|)$ can be obtained as

\begin{array}[]{rl}\mathcal{E}\left(|j\rangle\langle k|\right)=&\mathcal{E}(|+\rangle\langle+|)+\mathrm{i}\mathcal{E}(|-\rangle\langle-|)-\dfrac{1+\mathrm{i}}{2}\mathcal{E}(|j\rangle\langle j|)-\dfrac{1+\mathrm{i}}{2}\mathcal{E}(|k\rangle\langle k|).\\ \end{array}

(110)

Using natural basis states which we define as all the states needed by the RHS of (109), Wang et al. [13] have proved that $B$ is a permutation matrix. Thus, the error analysis in Step 2 becomes

\displaystyle\|\hat{D}-D\|

\displaystyle=\left\|\operatorname{vec}^{-1}\left(B^{\dagger}\operatorname{vec}(\hat{A})\right)-\operatorname{vec}^{-1}\left(B^{\dagger}\operatorname{vec}(A)\right)\right\|=\left\|B^{\dagger}(\operatorname{vec}(\hat{A})-\operatorname{vec}(A))\right\|=\|\hat{A}-A\|,

(111)

and the error analysis in other steps are the same as our analysis. Therefore, we obtain a tighter upper bound as $O\left(\frac{d^{3}\operatorname{Tr}(F)}{\sqrt{N}}\right)$ .

VI-C On the different types of input quantum states and measurement operators

In practice, the different types of input quantum states $M$ and measurement operators $J$ may vary and here we discuss the MSE scaling with respect to $M$ and $J$ . We first allow $M$ to change, with the dimension $d$ , the measurement operators and the copy number of each input state $N$ given. The different input states are randomly generated as [55, 56] or as [57] for truncated coherent states. Coherent states are more straightforward to be prepared than Fock states in quantum optics and have been used in QPT in [58, 59]. Each scenario with $N$ copies generates the input states $\rho_{m}^{\text{in}}$ according to a certain probability distribution and the total copies are $N_{t}=NM$ , based on which $\mathrm{E}(\cdot)$ denotes the expectation, different from $\mathbb{E}(\cdot)$ in Theorem 1. Define $\mathrm{f}_{m}\triangleq\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)-\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)\right)$ , which is i.i.d. with respect to the subscript $i$ . Thus we have

	$\displaystyle\quad\mathrm{E}\left\{\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\right\}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}M\mathrm{E}\Big{[}\Big{(}\sum_{m=1}^{M}\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)+\mathrm{f}_{m}^{}\Big{)}\cdot\Big{(}\mathrm{E}\Big{(}\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\Big{)}+\mathrm{f}_{m}^{T}\Big{)}\Big{)}^{-1}\Big{]}\Big{\}}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}M\mathrm{E}\Big{[}\Big{(}M\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\right)+\sum_{m=1}^{M}\mathrm{f}_{m}^{}\mathrm{f}_{m}^{T}\Big{)}^{-1}\Big{]}\Big{\}}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}\Big{[}\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)\mathrm{E}\Big{(}\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\Big{)}+\mathrm{E}\Big{(}\frac{1}{M}\sum_{m=1}^{M}\mathrm{f}_{m}^{}\mathrm{f}_{m}^{T}\Big{)}\Big{]}^{-1}\Big{\}}.$		(112)

According to the central limit theorem, as $M$ tends to infinity, $\mathrm{E}(\frac{1}{M}\sum_{m=1}^{M}\mathrm{f}_{m}^{*}\mathrm{f}_{m}^{T})$ converges to a fixed matrix [57]. Therefore the expectation of $\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)$ is

\displaystyle\mathrm{E}\left\{\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\right\}=O\left(\frac{1}{{M}}\right).

(113)

In Step 1, since $N$ is given, using (46), we have

\displaystyle\mathrm{E}\left\|\operatorname{col}_{m}(\hat{A}^{T})-\operatorname{col}_{m}(A^{T})\right\|^{2}=O\left(1\right).

(114)

In Step 2, assume that the SVD of $\left(V^{*}V^{T}\right)^{-1}V^{*}$ is

\left(V^{*}V^{T}\right)^{-1}V^{*}=U\left[\Sigma_{1},0\right]W^{\dagger},

(115)

where $\Sigma_{1}$ is a $d^{2}\times d^{2}$ diagonal matrix, and $U$ and $W$ are unitary. Using (49) and (113), $\mathrm{E}\left\|\Sigma_{1}\right\|^{2}=O\left(\frac{1}{{M}}\right)$ . Let

h_{m}=W^{\dagger}\left(\operatorname{col}_{m}(\hat{A})-\operatorname{col}_{m}(A)\right)=\left[\begin{array}[]{l}h_{m}^{(1)}\\ h_{m}^{(2)}\end{array}\right],

(116)

where $h_{m}^{(1)}$ is a $d^{2}\times 1$ vector, $h_{m}^{(2)}$ is a $(M-d^{2})\times 1$ vector and each element in $h_{m}$ scales as $O(1)$ from (114). Then for the error in Step 2, we have

	$\displaystyle\quad\mathrm{E}\left\\|R^{T}\left(I_{d^{2}}\otimes\left(V^{}V^{T}\right)^{-1}V^{}\right)(\operatorname{vec}(\hat{A})-\operatorname{vec}(A))\right\\|^{2}$
	$\displaystyle=\mathrm{E}\sum_{m=1}^{d^{2}}\left\\|\left(V^{}V^{T}\right)^{-1}V^{}\left(\operatorname{col}_{m}(\hat{A})-\operatorname{col}_{m}(A)\right)\right\\|^{2}$
	$\displaystyle=\mathrm{E}\sum_{m=1}^{d^{2}}\left\\|U\left[\Sigma_{1},0\right]\left[\begin{array}[]{l}h_{m}^{(1)}\\ h_{m}^{(2)}\end{array}\right]\right\\|^{2}$
	$\displaystyle=\mathrm{E}\sum_{m=1}^{d^{2}}\left\\|\Sigma_{1}h_{m}^{(1)}\right\\|^{2}=O\left(\frac{1}{{M}}\right).$		(117)

Since Steps 3–4 are not related to the scaling on $M$ , the final error scaling on $M$ is

\mathrm{E}\|\hat{X}-X\|=O\left(\frac{1}{\sqrt{M}}\right).

(118)

From our error upper bound in (67), using (113), we have

\displaystyle\mathrm{E}\left\{M\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\right\}=O\left(1\right),

(119)

and thus the scaling of (67) is $O\left(1\right)$ with given $N$ and $d$ . But the accurate scaling on $M$ is $\mathrm{E}\|\hat{X}-X\|=O\left(\frac{1}{\sqrt{M}}\right)$ . Therefore, our error upper bound (67) is not always tight with respect to $M$ .

Then we consider the MSE scaling of our TSS on $J$ . With given input states, $N$ and $d$ , we can also generate different measurement operators according to certain probability distributions. Similarly, we can prove

\mathrm{E}\left\|\operatorname{col}_{m}(\hat{A}^{T})-\operatorname{col}_{m}(A^{T})\right\|^{2}=O(1),

(120)

and

\displaystyle\mathrm{E}\left\{J\operatorname{Tr}\left(\left(C^{\dagger}C\right)^{-1}\right)\right\}=O\left(1\right).

(121)

Therefore, our error upper bound (67) is tight with respect to $J$ .

VII Numerical examples

To perform measurement on the output states, there are different measurement bases such as SIC-POVM [53], MUB measurement [52], and Cube bases [60]. In this section, we use Cube bases, as it is relatively easy to be implemented in experiments. For one-qubit systems, the Cube bases are $\left\{\frac{I\pm\sigma_{x}}{2},\frac{I\pm\sigma_{y}}{2},\frac{I\pm\sigma_{z}}{2}\right\}$ where $\sigma_{x}$ , $\sigma_{y}$ , $\sigma_{z}$ are Pauli matrices. For multi-qubit systems, the Cube bases are the tensor products of one-qubit Cube bases.

VII-A Performance illustration

We consider a $4$ -dimensional quantum system and the TP process matrix $X$ is determined by $\{\mathcal{A}_{i}\}_{i=1}^{3}$ where

		$\displaystyle\mathcal{A}_{1}=U_{1}\operatorname{diag}(0.5,0.4,0,0),$		(122)
		$\displaystyle\mathcal{A}_{2}=U_{2}\operatorname{diag}(0.1,0.2,0,0),$
		$\displaystyle\mathcal{A}_{3}=U_{3}\sqrt{I-\mathcal{A}_{1}^{\dagger}\mathcal{A}_{1}-\mathcal{A}_{2}^{\dagger}\mathcal{A}_{2}},$

where $U_{1}$ , $U_{2}$ and $U_{3}$ are random unitary matrices [55, 56].

We use four different sets of probe states. The first set is SIC state which is given in (133) in Appendix B and $M=16$ . The second set is MUB states which is given in (135) in Appendix B and $M=20$ . The third one is random states where we randomly generate $20$ input states using the algorithm in [55, 56]. The fourth set is natural basis states in the RHS of (109). The result of MSE (mean squared error) versus the total number of copies $N_{t}=MN$ is shown in Fig. 2(a), where $\operatorname{MSE}=\mathbb{E}||\hat{X}-X||^{2}$ and we make all of the four sets have the same $N_{t}$ . For each number of copies $N_{t}$ , we repeat our algorithm $20$ times and obtain the average MSE and error bars. For all the four sets of input states, the scaling of MSE is basically $O(1/N_{t})$ which satisfies Theorem 1. Moreover, the MSEs of SIC states, MUB states, natural basis states are all smaller than random states. The MSEs of SIC states and MUB states are quite similar and both smaller than natural basis states. For non-TP processes, we consider the same $\mathcal{A}_{1},\mathcal{A}_{2}$ and different $\mathcal{A}_{3}$ as

\mathcal{A}_{3}=U_{3}\sqrt{U_{4}\operatorname{diag}(1,0.8,0.7,0.5)U_{4}^{\dagger}-\mathcal{A}_{1}^{\dagger}\mathcal{A}_{1}-\mathcal{A}_{2}^{\dagger}\mathcal{A}_{2}},

(123)

where $U_{4}$ is a random unitary matrix. The result for non-TP processes is shown in Fig. 2(b) which is similar to Fig. 2(a). The MSE of non-TP processes is a little smaller than that of TP processes because the constraint for TP processes is stricter than that for non-TP processes.

We also test our algorithm on the IBM Quantum device [61]. Here we perform QPT on the CNOT gate defined as

\text{CNOT}=\left[\begin{array}[]{llll}1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ 0&1&0&0\end{array}\right],

(124)

in the IBM Quantum Composer. For the input states, we use states $\left\{\frac{I+\sigma_{x}}{2},\frac{I+\sigma_{y}}{2},\frac{I\pm\sigma_{z}}{2}\right\}^{\otimes 2}$ which can be generated by applying Hadamard and $S^{\dagger}$ gates to qubits initialized in the $|0\rangle$ state [62]. In addition, we use Cube measurements. Fig. 3 is an example of IBM Quantum Composer where the input state is $\frac{I-\sigma_{z}}{2}\otimes\frac{I+\sigma_{x}}{2}$ and the measurement operators are $\frac{I\pm\sigma_{y}}{2}\otimes\frac{I\pm\sigma_{x}}{2}$ .

If we have prior knowledge that the process is unitary, some references have proposed effective identification algorithms [63, 64]. Here we assume that we do not have prior knowledge that the unknown process is in essence unitary. We perform these numerical examples using the MATLAB on a classical computer and the $ibmq\_qasm\_simulator$ . The total numbers of copies are $N_{t}=1.08\times 10^{5},5.40\times 10^{5},2.70\times 10^{6},1.35\times 10^{7}$ and the experiments are repeated 5 times. Then we apply our TSS algorithm and the results are shown in Fig. 4. Besides MSE, we also consider another common fidelity metric, defined as [23]

F(\hat{X},X)\triangleq\left[\operatorname{Tr}\sqrt{\sqrt{\hat{X}}X\sqrt{\hat{X}}}\right]^{2}/\left[\operatorname{Tr}\left(X\right)\operatorname{Tr}\left(\hat{X}\right)\right],

(125)

and the correspond infidelity is defined as $1-F(\hat{X},X)$ . From Fig. 4, the scalings of the MSEs from $ibmq\_qasm\_simulator$ and simulation are also both $O(1/N_{t})$ which satisfies Theorem 1. The scalings of the infidelities are both $O(1/\sqrt{N_{t}})$ because the CNOT process is rank-deficient. A similar infidelity scaling has also been studied in QST [51, 65]. In addition, the errors between simulation and $ibmq\_qasm\_simulator$ are close. Then we apply our algorithm on the $ibmq\_quito$ 5-qubit system. The MSE and fidelity are presented in Fig. 5 where the total number of copies are $N_{t}=4.32\times 10^{5},5.76\times 10^{5},7.20\times 10^{5},8.64\times 10^{5},10.08\times 10^{5},11.52\times 10^{5}$ , respectively, and we repeat $3$ times. These results are both worse than $ibmq\_qasm\_simulator$ and simulation because there exist state preparation error and readout error in $ibmq\_quito$ [61]. The average CNOT error is $1.253\times 10^{-2}$ and average readout error is $4.430\times 10^{-2}$ for $ibmq\_quito$ . Thus, the approximate fidelity is $(1-1.253\times 10^{-2})(1-4.430\times 10^{-2})=0.9437$ .

VII-B Performance comparison

Here, we compare the performance between our TSS algorithm and the convex optimization method in [14]. Baldwin et al. [14] formulated the TP QPT problem as a convex optimization problem

$\displaystyle\min$	$\displaystyle\sum_{m,l}\left\|\mathcal{g}_{ml}-\operatorname{Tr}\left(\mathcal{O}_{ml}^{\dagger}\hat{X}\right)\right\|^{2}$	(126)
s.t.	$\displaystyle\sum_{j,k}\hat{X}_{jk}E_{k}^{\dagger}E_{j}=I,$
	$\displaystyle\hat{X}=\hat{X}^{\dagger},\quad\hat{X}\geq 0,$

where $\mathcal{g}_{ml}$ is measurement data, $\mathcal{O}_{ml}$ is a constant matrix which is the tensor product between the transpose of the $m$ -th input states and $l$ -th measurement bases where $1\leq m\leq M$ , $1\leq l\leq L$ , and $\hat{X}$ is the estimated process matrix. The detailed description can be found in equation (13) in [14]. This problem can be solved as a convex semidefinite program, which is realized by CVX [66, 67].

We randomly generate input quantum states as [55, 56] and $M=d(d+1)$ . With given number of copies $N=3\times 10^{4}$ for each output state and randomly generated process as in (122), a comparison of results is shown in Fig. 6 where the horizontal axis is the number of qubits $N_{q}$ , the left vertical axis is the MSE on a logarithmic scale and the right vertical axis is the running time $T$ (in seconds) on a logarithmic scale. We repeat the simulation $10$ times and obtain the average MSE, the average running time and error bars.

When the two algorithms have a similar MSE (we set the error of convex optimization a little larger than TSS by changing CVX precision), the running time of TSS algorithm is smaller than convex optimization by approximately three orders of magnitude for three-qubit systems, indicating that our TSS algorithm is more efficient than the convex optimization method. Since $M=d(d+1)$ and we choose Cube bases where $L=6^{N_{q}}$ , the computational complexity for our algorithm is $O(96^{N_{q}})$ . For the simulated running time of our algorithm, since our computational complexity analysis is asymptotical on $d$ , we fit the rightmost three points and the slope is $2.0843$ , which is close to the theoretical value $\log_{10}(96)=1.9823$ . The slight difference might be attributed to fact that the qubit number is not large enough. The storage requirement is also labeled by the black numbers in Fig. 6. The storage requirement for our algorithm is $O(ML)$ . Since we use Cube bases of measurement, $L=6^{N_{q}},M=4^{N_{q}}+2^{N_{q}}$ and thus $ML=24^{N_{q}}+12^{N_{q}}$ . As $N_{q}$ increases from $1$ to $5$ , $ML$ are $36$ , $720$ , $15552$ , $352512$ and $8211456$ , respectively. For the convex optimization method, the main storage requirement is the cost of storing $\mathcal{O}_{ml}$ , which is a $d^{2}\times d^{2}$ matrix and the total number of these matrices is $ML$ . Thus, the storage requirement is $O(MLd^{4})$ . As $N_{q}$ increases from $1$ to $3$ , $MLd^{4}=384^{N_{q}}+192^{Nq}$ are $576$ , $184320$ and $63700992$ , respectively. Thus, our algorithm needs less storage than the convex optimization method.

We then compare the two algorithms by changing the types of input states $M$ with $d=4$ , where we randomly generate input states as in [55, 56] and the number of copies for each output state is $N=3^{2}\times 10^{4}$ . The comparison results are shown in Fig. 7 where the horizontal axis is the different types of input states $M$ , the left vertical axis is the MSE and the right vertical axis is the running time $T$ , all on a logarithmic scale. From Fig. 7, the scaling between MSE and $M$ is $-1.1474$ , which is close to $-1$ and consistent with (118). This scaling is different from the scaling on $N_{t}$ in Fig. 2 because here $N$ is given and $M$ increases, while in Fig. 2, $M$ is given and $N$ increases. For the simulated running time of our algorithm, we fit the rightmost three points to obtain the slope $1.0090$ , which is close to the theoretical value $1$ . For the convex optimization method, the slope of the fitted line of the rightmost three points is $1.9186$ , which is larger than that of our TSS algorithm, indicating that our TSS algorithm is more efficient than the convex optimization method.

VIII Conclusion

In this paper, we have introduced an analytical two-stage solution (TSS) applicable to both trace-preserving and non-trace-preserving QPT. Leveraging the natural basis, we have utilized the tensor structure of the coefficient matrix and provided insights into the computational complexity and storage requirements for our algorithm. Our contributions include the establishment of an analytical upper bound for errors and the optimization of input quantum states and measurement operators based on this error upper bound and condition number. The effectiveness and theoretical soundness of our algorithm have been demonstrated through numerical examples and testing on an IBM Quantum device. Furthermore, we have benchmarked our TSS algorithm against the convex optimization method outlined in [14], and the results show that our algorithm is more efficient in terms of both time and space costs. Further work will focus on extending our algorithm to AAPT.

Acknowledgments

The authors would like to thank Prof. Lloyd C. L. Hollenberg for the helpful discussion. This research was supported by the University of Melbourne through the establishment of the IBM Quantum Network Hub at the University. Shuixin Xiao would like to thank the support from the IEEE Control Systems Society Graduate Collaboration Fellowship.

Appendix A Several propositions and lemmas

Here, we give some propositions and lemmas which are utilized in proving Theorem 1.

Proposition 3

For any $A\in\mathbb{C}^{n\times k}$ and $b\in\mathbb{C}^{mk}$ , we have

\left\|\left(I_{m}\otimes A\right)b\right\|\leq\|A\|\|b\|

(127)

Proof:

Let

b=\left[b_{1}^{T},b_{2}^{T},\cdots,b_{m}^{T}\right]^{T},

(128)

where $b_{i}$ is a $k\times 1$ vector. Therefore,

\displaystyle\left\|\!\left(I_{m}\otimes A\right)\!b\right\|^{2}\!=\!\sum_{i=1}^{m}\left\|Ab_{i}\right\|^{2}\!\leq\!\|A\|^{2}\sum_{i=1}^{m}\left\|b_{i}\right\|^{2}\!=\!\|A\|^{2}\|b\|^{2}.

(129)

∎

Lemma 2

[68] Let $\mathbb{H}_{A}$ and $\mathbb{H}_{B}$ be finite-dimensional Hilbert spaces of dimensions $d_{A}$ and $d_{B}$ , respectively, and let $X\in\mathbb{H}_{A}\otimes\mathbb{H}_{B}$ . Then for any unitarily invariant norm that is multiplicative over tensor products, the partial trace satisfies the norm inequality

\left\|\operatorname{Tr}_{A}(X)\right\|\leq\frac{d_{A}}{\left\|I_{A}\right\|}\|X\|,

(130)

where $I_{A}$ is the identity operator.

Lemma 3

([69] Theorem 8.1 and Theorem 28.3) Let $X$ , $Y$ be Hermitian matrices with eigenvalues $\lambda_{1}(X)\geq\cdots\geq\lambda_{n}(X)$ and $\lambda_{1}(Y)\geq\cdots\geq\lambda_{n}(Y)$ , respectively. Then

\max_{j}|\lambda_{j}(X)-\lambda_{j}(Y)|\leq||X-Y||,

(131)

and

\sum_{j=1}^{n}\left(\lambda_{j}(X)-\lambda_{j}(Y)\right)^{2}\leq||X-Y||^{2}.

(132)

Appendix B SIC states and MUB states

For input states achieving the lower bounds of $M\operatorname{Tr}\big{(}(V^{*}V^{T})^{-1}\big{)}$ and $\operatorname{cond}(V)$ , we give two examples, SIC states and MUB states, as the following propositions.

Proposition 4

$d$ -dimensional SIC states (when they exist) belong to $\underline{\mathrm{OIS}}(d,M)$ with the smallest $M$ as $M=d^{2}$ .

Proposition 5

$d$ -dimensional MUB states (when they exist) belong to $\underline{\mathrm{OIS}}(d,M)$ for $M=d(d+1)$ .

The proofs are similar to the optimal probe states in quantum detector tomography and can be found in [32].

For SIC-POVM in $d=4$ , Bengtsson [70] gave one expression ignoring overall phases and normalization in the natural basis as

\left[\begin{array}[]{cccccccccccccccc}x&x&x&x&\mathrm{i}&\mathrm{i}&-\mathrm{i}&-\mathrm{i}&\mathrm{i}&\mathrm{i}&-\mathrm{i}&-\mathrm{i}&\mathrm{i}&\mathrm{i}&-\mathrm{i}&-\mathrm{i}\\ 1&1&-1&-1&x&x&x&x&\mathrm{i}&-\mathrm{i}&\mathrm{i}&-\mathrm{i}&1&-1&1&-1\\ 1&-1&1&-1&1&-1&1&-1&x&x&x&x&-\mathrm{i}&\mathrm{i}&\mathrm{i}&-\mathrm{i}\\ 1&-1&-1&1&-\mathrm{i}&\mathrm{i}&\mathrm{i}&-\mathrm{i}&-1&1&1&-1&x&x&x&x\end{array}\right],

(133)

where $x=\sqrt{2+\sqrt{5}}$ . Let $\left|\psi_{n}^{(\text{SIC})}\right\rangle$ be the $n$ -th column of (133). In this paper, we call the set of

\rho_{n}=\frac{\left|\psi_{n}^{(\text{SIC})}\right\rangle\left\langle\psi_{n}^{(\text{SIC})}\right|}{\operatorname{Tr}\left(\left|\psi_{n}^{(\text{SIC})}\right\rangle\left\langle\psi_{n}^{(\text{SIC})}\right|\right)}

as SIC states and $\left\{P_{n}=\frac{1}{d}\rho_{n}\right\}_{n=1}^{d^{2}}$ as SIC-POVM.

For $d=4$ , five MUB measurement sets are

\left\{\!|\psi_{n}^{(\text{MUB})}\rangle\!\right\}\!=\!\left\{\!\left\{\left|\psi_{n}^{A}\right\rangle\!\right\}\!,\!\left\{\left|\psi_{n}^{B}\right\rangle\!\right\}\!,\!\left\{\left|\psi_{n}^{C}\right\rangle\!\right\}\!,\!\left\{\left|\psi_{n}^{D}\right\rangle\!\right\}\!,\!\left\{\left|\psi_{n}^{E}\right\rangle\!\right\}\!\right\},

(134)

and

	$\displaystyle\left\{\left\|\psi_{n}^{A}\right\rangle\right\}=\{\|00\rangle,\|01\rangle,\|10\rangle,\|11\rangle\},$
	$\displaystyle\left\{\left\|\psi_{n}^{B}\right\rangle\right\}=\{\|R\pm\rangle,\|L\pm\rangle\},$
	$\displaystyle\left\{\left\|\psi_{n}^{C}\right\rangle\right\}=\{\|\pm R\rangle,\|\pm L\rangle\},$
	$\displaystyle\left\{\left\|\psi_{n}^{D}\right\rangle\right\}=\left\{\frac{1}{\sqrt{2}}(\|R0\rangle\pm\mathrm{i}\|L1\rangle),\frac{1}{\sqrt{2}}(\|R1\rangle\pm\mathrm{i}\|L0\rangle)\right\},$
	$\displaystyle\left\{\left\|\psi_{n}^{E}\right\rangle\right\}=\left\{\frac{1}{\sqrt{2}}(\|RR\rangle\pm\mathrm{i}\|LL\rangle),\frac{1}{\sqrt{2}}(\|RL\rangle\pm\mathrm{i}\|LR\rangle)\right\},$		(135)

where $|\pm\rangle=(|0\rangle\pm|1\rangle)/\sqrt{2}$ , $|R\rangle=(|0\rangle-{\mathrm{i}}|1\rangle)/\sqrt{2}$ , and $|L\rangle=(|0\rangle+{\mathrm{i}}|1\rangle)/\sqrt{2}$ in the natural basis. Adamson et al. [51] utilized these Mutually Unbiased Bases in the QST. In this paper, we call $\rho_{n}=|\psi_{n}^{(\text{MUB})}\rangle\langle\psi_{n}^{(\text{MUB})}|$ a MUB state and $P_{n}=|\psi_{n}^{(\text{MUB})}\rangle\langle\psi_{n}^{(\text{MUB})}|$ a MUB measurement operator.

References

[1] D. P. DiVincenzo, “Quantum computation,” Science, vol. 270, no. 5234, pp. 255–261, 1995.
[2] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information. Cambridge University Press, 2010.
[3] N. Gisin and R. Thew, “Quantum communication,” Nature Photonics, vol. 1, no. 3, p. 165, 2007.
[4] M.-H. Hsieh and M. M. Wilde, “Trading classical communication, quantum communication, and entanglement in quantum Shannon theory,” IEEE Transactions on Information Theory, vol. 56, no. 9, pp. 4705–4730, 2010.
[5] C. L. Degen, F. Reinhard, and P. Cappellaro, “Quantum sensing,” Reviews of Modern Physics, vol. 89, p. 035002, 2017.
[6] D. Dong and I. R. Petersen, “Quantum control theory and applications: a survey,” IET Control Theory & Applications, vol. 4, no. 12, pp. 2651–2671, 2010.
[7] D. Dong and I. R. Petersen, “Quantum estimation, control and learning: Opportunities and challenges,” Annual Reviews in Control, vol. 54, pp. 243–251, 2022.
[8] D. Dong and I. R. Petersen, Learning and Robust Control in Quantum Technology. Springer Nature Switzerland AG, 2023.
[9] D. Burgarth and K. Yuasa, “Quantum system identification,” Physical Review Letters, vol. 108, p. 080502, 2012.
[10] A. Sone and P. Cappellaro, “Hamiltonian identifiability assisted by a single-probe measurement,” Physical Review A, vol. 95, p. 022335, 2017.
[11] Y. Wang, D. Dong, A. Sone, I. R. Petersen, H. Yonezawa, and P. Cappellaro, “Quantum Hamiltonian identifiability via a similarity transformation approach and beyond,” IEEE Transactions on Automatic Control, vol. 65, no. 11, pp. 4632–4647, 2020.
[12] J. Zhang and M. Sarovar, “Quantum Hamiltonian identification from measurement time traces,” Physical Review Letters, vol. 113, p. 080401, 2014.
[13] Y. Wang, D. Dong, B. Qi, J. Zhang, I. R. Petersen, and H. Yonezawa, “A quantum Hamiltonian identification algorithm: Computational complexity and error analysis,” IEEE Transactions on Automatic Control, vol. 63, no. 5, pp. 1388–1403, 2018.
[14] C. H. Baldwin, A. Kalev, and I. H. Deutsch, “Quantum process tomography of unitary and near-unitary maps,” Physical Review A, vol. 90, p. 012110, 2014.
[15] M. Zorzi, F. Ticozzi, and A. Ferrante, “Minimal resources identifiability and estimation of quantum channels,” Quantum Information Processing, vol. 13, no. 3, pp. 683–707, 2014.
[16] M. Zorzi, F. Ticozzi, and A. Ferrante, “Estimation of quantum channels: Identifiability and ML methods,” in 2012 51st IEEE Conference on Decision and Control (CDC), pp. 1674–1679, 2012.
[17] G. C. Knee, E. Bolduc, J. Leach, and E. M. Gauger, “Quantum process tomography via completely positive and trace-preserving projection,” Physical Review A, vol. 98, p. 062336, 2018.
[18] T. Surawy-Stepney, J. Kahn, R. Kueng, and M. Guţă, “Projected least-squares quantum process tomography,” Quantum, vol. 6, p. 844, 2022.
[19] Z. Ji, G. Wang, R. Duan, Y. Feng, and M. Ying, “Parameter estimation of quantum channels,” IEEE Transactions on Information Theory, vol. 54, no. 11, pp. 5172–5185, 2008.
[20] J. Zhang and M. Sarovar, “Identification of open quantum systems from observable time traces,” Physical Review A, vol. 91, no. 5, p. 052121, 2015.
[21] M. Guţă and N. Yamamoto, “System identification for passive linear quantum systems,” IEEE Transactions on Automatic Control, vol. 61, no. 4, pp. 921–936, 2016.
[22] X.-L. Huang, J. Gao, Z.-Q. Jiao, Z.-Q. Yan, Z.-Y. Zhang, D.-Y. Chen, X. Zhang, L. Ji, and X.-M. Jin, “Reconstruction of quantum channel via convex optimization,” Science Bulletin, vol. 65, no. 4, pp. 286–292, 2020.
[23] I. Bongioanni, L. Sansoni, F. Sciarrino, G. Vallone, and P. Mataloni, “Experimental quantum process tomography of non-trace-preserving maps,” Physical Review A, vol. 82, p. 042307, 2010.
[24] G. A. White, C. D. Hill, F. A. Pollock, L. C. Hollenberg, and K. Modi, “Demonstration of non-Markovian process characterisation and control on a quantum processor,” Nature Communications, vol. 11, p. 6301, 2020.
[25] G. A. White, F. A. Pollock, L. C. Hollenberg, K. Modi, and C. D. Hill, “Non-markovian quantum process tomography,” PRX Quantum, vol. 3, p. 020344, 2022.
[26] H. Wang, W. Zheng, N. Yu, K. Li, D. Lu, T. Xin, C. Li, Z. Ji, D. Kribs, B. Zeng, X. Peng, and J. Du, “Quantum state and process tomography via adaptive measurements,” Science China Physics, Mechanics & Astronomy, vol. 59, no. 10, p. 100313, 2016.
[27] Y. S. Teo, B.-G. Englert, J. Řeháček, and Z. Hradil, “Adaptive schemes for incomplete quantum process tomography,” Physical Review A, vol. 84, p. 062125, 2011.
[28] I. A. Pogorelov, G. I. Struchalin, S. S. Straupe, I. V. Radchenko, K. S. Kravtsov, and S. P. Kulik, “Experimental adaptive process tomography,” Physical Review A, vol. 95, p. 012302, 2017.
[29] I. L. Chuang and M. A. Nielsen, “Prescription for experimental determination of the dynamics of a quantum black box,” Journal of Modern Optics, vol. 44, no. 11-12, pp. 2455–2467, 1997.
[30] J. F. Poyatos, J. I. Cirac, and P. Zoller, “Complete characterization of a quantum process: The two-bit quantum gate,” Physical Review Letters, vol. 78, pp. 390–393, 1997.
[31] J. Fiurášek and Z. Hradil, “Maximum-likelihood estimation of quantum processes,” Physical Review A, vol. 63, p. 020101, 2001.
[32] S. Xiao, Y. Wang, D. Dong, and J. Zhang, “Optimal and two-step adaptive quantum detector tomography,” Automatica, vol. 141, p. 110296, 2022.
[33] R. A. Horn and C. R. Johnson, Matrix Analysis. Cambridge University Press, 2012.
[34] J. Watrous, The Theory of Quantum Information. Cambridge University Press, 2018.
[35] D. W. Leung, “Choi’s proof as a recipe for quantum process tomography,” Journal of Mathematical Physics, vol. 44, no. 2, pp. 528–533, 2003.
[36] G. M. D’Ariano and P. Lo Presti, “Quantum tomography for measuring experimentally the matrix elements of an arbitrary quantum operation,” Physical Review Letters, vol. 86, pp. 4195–4198, 2001.
[37] J. B. Altepeter, D. Branning, E. Jeffrey, T. C. Wei, P. G. Kwiat, R. T. Thew, J. L. O’Brien, M. A. Nielsen, and A. G. White, “Ancilla-assisted quantum process tomography,” Physical Review Letters, vol. 90, p. 193601, 2003.
[38] G. M. D’Ariano and P. Lo Presti, “Imprinting complete information about a quantum channel on its output state,” Physical Review Letters, vol. 91, p. 047902, 2003.
[39] M. Mohseni and D. A. Lidar, “Direct characterization of quantum dynamics,” Physical Review Letters, vol. 97, p. 170501, 2006.
[40] M. Mohseni and D. A. Lidar, “Direct characterization of quantum dynamics: General theory,” Physical Review A, vol. 75, p. 062331, 2007.
[41] Z.-W. Wang, Y.-S. Zhang, Y.-F. Huang, X.-F. Ren, and G.-C. Guo, “Experimental realization of direct characterization of quantum dynamics,” Physical Review A, vol. 75, p. 044304, 2007.
[42] M. Mohseni, A. T. Rezakhani, and D. A. Lidar, “Quantum-process tomography: Resource analysis of different strategies,” Physical Review A, vol. 77, p. 032322, 2008.
[43] M. Zorzi, F. Ticozzi, and A. Ferrante, “Minimum relative entropy for quantum estimation: Feasibility and general solution,” IEEE Transactions on Information Theory, vol. 60, no. 1, pp. 357–367, 2014.
[44] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, “Sample-optimal tomography of quantum states,” IEEE Transactions on Information Theory, vol. 63, no. 9, pp. 5628–5641, 2017.
[45] M. Berta, J. M. Renes, and M. M. Wilde, “Identifying the information gain of a quantum measurement,” IEEE Transactions on Information Theory, vol. 60, no. 12, pp. 7987–8006, 2014.
[46] B. Qi, Z. Hou, L. Li, D. Dong, G.-Y. Xiang, and G.-C. Guo, “Quantum state tomography via linear regression estimation,” Scientific Reports, vol. 3, p. 3496, 2013.
[47] B. Mu, H. Qi, I. R. Petersen, and G. Shi, “Quantum tomography by regularized linear regressions,” Automatica, vol. 114, p. 108837, 2020.
[48] J. A. Smolin, J. M. Gambetta, and G. Smith, “Efficient method for computing the maximum-likelihood quantum state from measurements with additive Gaussian noise,” Physical Review Letters, vol. 108, p. 070502, 2012.
[49] H. Sedrakyan and N. Sedrakyan, Algebraic Inequalities. Springer, 2018.
[50] S. Xiao, Y. Wang, D. Dong, and J. Zhang, “Optimal quantum detector tomography via linear regression estimation,” in 2021 60th IEEE Conference on Decision and Control (CDC), pp. 4140–4145, 2021.
[51] R. B. A. Adamson and A. M. Steinberg, “Improving quantum state estimation with mutually unbiased bases,” Physical Review Letters, vol. 105, p. 030406, 2010.
[52] T. Durt, B.-G. Englert, I. Bengtsson, and K. Życzkowski, “On mutually unbiased bases,” International Journal of Quantum Information, vol. 08, no. 04, pp. 535–640, 2010.
[53] J. M. Renes, R. Blume-Kohout, A. J. Scott, and C. M. Caves, “Symmetric informationally complete quantum measurements,” Journal of Mathematical Physics, vol. 45, no. 6, pp. 2171–2180, 2004.
[54] A. Miranowicz, K. Bartkiewicz, J. Peřina, M. Koashi, N. Imoto, and F. Nori, “Optimal two-qubit tomography based on local and global measurements: Maximal robustness against errors as described by condition numbers,” Physical Review A, vol. 90, p. 062123, 2014.
[55] J. A. Miszczak, “Generating and using truly random quantum states in Mathematica,” Computer Physics Communications, vol. 183, no. 1, pp. 118–124, 2012.
[56] N. Johnston, “QETLAB: A MATLAB toolbox for quantum entanglement, version 0.9,” Jan. 2016.
[57] Y. Wang, S. Yokoyama, D. Dong, I. R. Petersen, E. H. Huntington, and H. Yonezawa, “Two-stage estimation for quantum detector tomography: Error analysis, numerical and experimental results,” IEEE Transactions on Information Theory, vol. 67, no. 4, pp. 2293–2307, 2021.
[58] M. Lobino, D. Korystov, C. Kupchak, E. Figueroa, B. C. Sanders, and A. I. Lvovsky, “Complete characterization of quantum-optical processes,” Science, vol. 322, no. 5901, pp. 563–566, 2008.
[59] S. Rahimi-Keshari, A. Scherer, A. Mann, A. T. Rezakhani, A. I. Lvovsky, and B. C. Sanders, “Quantum process tomography with coherent states,” New Journal of Physics, vol. 13, no. 1, p. 013006, 2011.
[60] M. D. de Burgh, N. K. Langford, A. C. Doherty, and A. Gilchrist, “Choice of measurement sets in qubit tomography,” Physical Review A, vol. 78, p. 052122, 2008.
[61] “IBM Quantum.” https://quantum-computing.ibm.com, 2021.
[62] A. Dang, G. A. White, L. C. Hollenberg, and C. D. Hill, “Process tomography on a 7-qubit quantum processor via tensor network contraction path finding,” arXiv preprint arXiv:2112.06364, 2021.
[63] G. Gutoski and N. Johnston, “Process tomography for unitary quantum channels,” Journal of Mathematical Physics, vol. 55, no. 3, p. 032201, 2014.
[64] Y. Wang, Q. Yin, D. Dong, B. Qi, I. R. Petersen, Z. Hou, H. Yonezawa, and G.-Y. Xiang, “Quantum gate identification: Error analysis, numerical results and optical experiment,” Automatica, vol. 101, pp. 269–279, 2019.
[65] B. Qi, Z. Hou, Y. Wang, D. Dong, H.-S. Zhong, L. Li, G.-Y. Xiang, H. M. Wiseman, C.-F. Li, and G.-C. Guo, “Adaptive quantum state tomography via linear regression estimation: Theory and two-qubit experiment,” npj Quantum Information, vol. 3, p. 19, 2017.
[66] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1.” http://cvxr.com/cvx, Mar. 2014.
[67] M. Grant and S. Boyd, “Graph implementations for nonsmooth convex programs,” in Recent Advances in Learning and Control (V. Blondel, S. Boyd, and H. Kimura, eds.), Lecture Notes in Control and Information Sciences, pp. 95–110, Springer-Verlag Limited, 2008. http://stanford.edu/~boyd/graph_dcp.html.
[68] D. A. Lidar, P. Zanardi, and K. Khodjasteh, “Distance bounds on quantum dynamics,” Physical Review A, vol. 78, p. 012308, 2008.
[69] R. Bhatia, Perturbation Bounds for Matrix Eigenvalues. Society for Industrial and Applied Mathematics, 2007.
[70] I. Bengtsson, “From SICs and MUBs to Eddington,” Journal of Physics: Conference Series, vol. 254, p. 012007, 2010.

	$\displaystyle\quad\Delta_{D}^{2}=\\|\hat{D}-D\\|^{2}$
	$\displaystyle=\left\\|\operatorname{vec}\left(\hat{D}\right)-\operatorname{vec}\left(D\right)\right\\|^{2}$
	$\displaystyle=\left\\|R^{T}\left(I\otimes\left(V^{}V^{T}\right)^{-1}V^{}\right)(\operatorname{vec}(\hat{A})-\operatorname{vec}(A))\right\\|^{2}$
	$\displaystyle\leq\left\\|\left(V^{}V^{T}\right)^{-1}V^{}\right\\|^{2}\\|\operatorname{vec}(\hat{A})-\operatorname{vec}(A)\\|^{2}$
	$\displaystyle=\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\Delta_{A}^{2},$		(50)

$\displaystyle\\|\hat{F}-F\\|$	$\displaystyle=\left\\|\operatorname{Tr}_{1}(\hat{G})-\operatorname{Tr}_{1}(G)\right\\|$	(53)
	$\displaystyle\leq\sqrt{d}\\|\hat{G}-G\\|$
	$\displaystyle=\sqrt{d}\Delta_{G}.$

	$\displaystyle\quad\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle=\sum_{i=1}^{d^{2}}\left\\|\left(I_{d}\otimes\tilde{F}^{1/2}\bar{F}^{-1/2}-I_{d}\otimes I_{d}\right)\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle\leq\left\\|\tilde{F}^{1/2}\bar{F}^{-1/2}-I_{d}\right\\|^{2}\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}$
	$\displaystyle\sim\frac{\operatorname{Tr}({F})}{4}\\|\hat{F}-F\\|^{2}\sim\frac{d\operatorname{Tr}({F})}{4}\\|G-\hat{G}\\|^{2},$		(64)

	$\displaystyle\quad\\|\hat{X}-\hat{G}\\|^{2}$
	$\displaystyle=\left\\|\sum_{i=1}^{d^{2}}\left(\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\right)\right\\|^{2}$
	$\displaystyle\leq\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\right\\|\right)^{2}$
	$\displaystyle=\Bigg{(}\sum_{i=1}^{d^{2}}\Big{\\|}\operatorname{vec}\left(\hat{Q}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}+\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}-\operatorname{vec}\left(\hat{S}_{i}\right)\operatorname{vec}\left(\hat{S}_{i}\right)^{\dagger}\Big{\\|}\Bigg{)}^{2}$
	$\displaystyle\leq\Bigg{(}\sum_{i=1}^{d^{2}}\left(\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)^{\dagger}\right\\|+\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\right)\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\Bigg{)}^{2}$
	$\displaystyle\leq\left(\sum_{i=1}^{d^{2}}\left(\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\right\\|+\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|\right)^{2}\right)\cdot\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)$
	$\displaystyle\leq 2\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)\right\\|^{2}+\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)\cdot\left(\sum_{i=1}^{d^{2}}\left\\|\operatorname{vec}\left(\hat{Q}_{i}\right)-\operatorname{vec}\left(\hat{S}_{i}\right)\right\\|^{2}\right)$
	$\displaystyle\sim d\left(\operatorname{Tr}({F})\right)^{2}\\|G-\hat{G}\\|^{2},$		(65)

	$\displaystyle\quad\mathrm{E}\left\{\operatorname{Tr}\left(\left(V^{*}V^{T}\right)^{-1}\right)\right\}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}M\mathrm{E}\Big{[}\Big{(}\sum_{m=1}^{M}\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)+\mathrm{f}_{m}^{}\Big{)}\cdot\Big{(}\mathrm{E}\Big{(}\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\Big{)}+\mathrm{f}_{m}^{T}\Big{)}\Big{)}^{-1}\Big{]}\Big{\}}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}M\mathrm{E}\Big{[}\Big{(}M\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\right)+\sum_{m=1}^{M}\mathrm{f}_{m}^{}\mathrm{f}_{m}^{T}\Big{)}^{-1}\Big{]}\Big{\}}$
	$\displaystyle=\frac{1}{M}\operatorname{Tr}\Big{\{}\Big{[}\mathrm{E}\left(\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{}\right)\mathrm{E}\Big{(}\operatorname{vec}\left(\rho_{m}^{\text{in}}\right)^{T}\Big{)}+\mathrm{E}\Big{(}\frac{1}{M}\sum_{m=1}^{M}\mathrm{f}_{m}^{}\mathrm{f}_{m}^{T}\Big{)}\Big{]}^{-1}\Big{\}}.$		(112)

A two-stage solution to quantum process tomography: error analysis and optimal design

Abstract

Index Terms:

I Introduction

II Preliminaries and quantum process tomography

II-A Vectorization function

II-B Quantum process tomography

Problem 1

III Tensor structure and two-stage quantum process tomography solution

III-A Tensor structure of QPT

Lemma 1

Proposition 1

Proof:

Remark 1

Problem 2

III-B Two-stage solution of QPT

Problem 2.1.1

Problem 2.2.2

Remark 2

IV Computational complexity and storage requirements

V Error analysis

Theorem 1

Proof:

V-A Error in Step 1

V-B Error in Step 2

V-C Error in Step 3

V-D Error in Step 4

VI Optimization of input states and measurement operators

VI-A Optimal input quantum states

Theorem 2

Proof:

Theorem 3

Proof:

VI-B Optimal measurement operators

Theorem 4

Proof:

Proposition 2

Proof:

Remark 3

Remark 4

VI-C On the different types of input quantum states and measurement operators

VII Numerical examples

VII-A Performance illustration

VII-B Performance comparison

VIII Conclusion

Acknowledgments

Appendix A Several propositions and lemmas

Proposition 3

Proof:

Lemma 2

Lemma 3

Appendix B SIC states and MUB states

Proposition 4

Proposition 5

References

Problem 2 $.1$

Problem 2 $.2$