Reliable Optimization of Arbitrary Functions over Quantum Measurements

Abstract

As the connection between classical and quantum worlds, quantum measurements play a unique role in the era of quantum information processing. Given an arbitrary function of quantum measurements, how to obtain its optimal value is often considered as a basic yet important problem in various applications. Typical examples include but are not limited to optimizing the likelihood functions in quantum measurement tomography, searching the Bell parameters in Bell-test experiments, and calculating the capacities of quantum channels. In this work, we propose reliable algorithms for optimizing arbitrary functions over the space of quantum measurements by combining the so-called Gilbert’s algorithm for convex optimization with certain gradient algorithms. With extensive applications, we demonstrate the efficacy of our algorithms with both convex and nonconvex functions.

1. Introduction

In quantum information science, numerous complex mathematical problems remain to be solved. Since the set of quantum states as well as quantum measurements form convex sets, various important tasks in this field, such as the calculation of ground state energy, violation of the Bell inequality, and the detection and quantification of quantum entanglement [1,2], conform to the framework of convex optimization theory. The primary tool in convex optimization is semidefinite programming (SDP) [3,4], which can be used to derive relaxed constraints and provide accurate solutions for a large number of computationally challenging tasks. However, serious drawbacks also exist for SDP including its slow computation speed and low accuracy. For instance, SDP can only compute up to four qubits in quantum state tomography (QST), while improved superfast algorithms [5] can quickly go up to eleven qubits with a higher precision. Consequently, developing more efficient algorithms in convex optimization is becoming more and more crucial as quantum technologies rapidly advance.

Recently, an efficient convex optimization algorithm [6] was proposed by Brierley et al. based on the so-called Gilbert’s algorithm [7]. Concurrently, Ref. [8] used Gilbert’s algorithm to investigate whether nonlocal relationships can be distinguished in polynomial time. In Ref. [9], Gilbert’s algorithm was employed as a tool to satisfy certain constraints, based on which two reliable convex optimization schemes over the quantum state space were proposed. In addition, some nonconvex optimization algorithms were also brought out for QST; for instance, the one in Ref. [10] is faster and more accurate as compared to previous approaches. One notices that all these studies concern only the optimization over quantum state space, with the consideration over quantum measurement space rarely being mentioned.

In fact, various important and meaningful problems related to quantum measurements exist in convex optimization, including, for example, searching the Bell parameters in Bell-test experiments [11], optimizing the correlation of quantum measurements under different measurement settings [12,13,14,15], and maximizing the likelihood functions in quantum measurement tomography. Meanwhile, characterization of quantum measurements forms the basis for quantum state tomography [16,17,18] and quantum process tomography [19,20,21]. Therefore, convex optimization over the quantum measurement space stands as an independent yet important problem in quantum information theory. However, the space of quantum measurements is much more complex as compared to the quantum state space since it is possible to produce an infinite variety of different measurement outcomes as long as the probabilities for these outcomes sum to one. Recently, Ref. [22] proposed a method to optimize over the measurement space based on SDP, but it fails to solve complex tasks due to the intrinsic problem with SDP. Worst of all, nonconvex functions [23] easily appear in the space of quantum measurements. Unlike convex functions, local optima might be found during the process of optimization. Hence, nonconvex optimization is regarded as more difficult than convex optimization. In this work, we propose two reliable algorithms for optimizing arbitrary functions over the space of quantum measurements by combining the so-called Gilbert’s algorithm for convex optimization with the direct-gradient (DG) algorithm as well as the accelerated projected gradient (APG) algorithm. With extensive applications, we demonstrate the efficacy of our algorithms with both convex and nonconvex functions.

This work is organized as follows: In Section 2, we propose two reliable algorithms for optimizing over quantum measurement space by combining Gilbert’s algorithm with the DG and APG algorithms, respectively. The universality of our method is demonstrated by several examples with both convex and nonconvex functions in Section 3. The last Section 4 provides the conclusions.

2. Function Optimization

In the quantum state space $\mathcal{Q}$, an arbitrary state $\rho$ should satisfy the conditions

$$\begin{array}{ccc}\hfill \rho & \ge & 0,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \mathrm{tr}\left(\rho \right)& =& 1.\hfill \end{array}$$

(2)

Given a smaller convex subset $\mathcal{C}\in \mathcal{Q}$, Gilbert’s algorithm can be used to approximately find the closest state ${\rho }^{\mathcal{C}}\in \mathcal{C}$ with respect to $\rho$ [9]. In general, for an arbitrary matrix M in the matrix space $\mathcal{M}$, we employ Gilbert’s algorithm to search for the closest quantum state ${\rho }^{\mathcal{Q}}\in \mathcal{Q}$ with respect to M. Throughout this work, let us denote the operation by using Gilbert’s algorithm as

$${\rho }^{\mathcal{Q}}\equiv \mathcal{S}\left(M\right).$$

(3)

Given experimental data, it is critical to identify the measurement settings that are most compatible with the data. Here, we consider the quantum measurement space $\Omega$ as all the positive operator-valued measures (POVMs). A quantum measurement device is characterized by a set of operators $\left\{{\Pi }_l\right\}$, which have to satisfy two constraints

$$\begin{array}{ccc}\hfill {\Pi }_l& \ge & 0,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \sum _{l=1}^L{\Pi }_l& =& \mathbb{I},\hfill \end{array}$$

(5)

where L is the total number of operators in the set. Denote a function $\mathcal{F}\left[\left\{{\Pi }_l\right\}\right]$ defined over the quantum measurement space $\Omega$. We assume that $\mathcal{F}\left[\left\{{\Pi }_l\right\}\right]$ is differentiable with the gradient $\nabla \mathcal{F}\left[\left\{{\Pi }_l\right\}\right]\equiv G\left[\left\{{\Pi }_l\right\}\right]$. The objective is to optimize $\mathcal{F}\left[\left\{{\Pi }_l\right\}\right]$ over the entire quantum measurement space, and we have

$$\begin{array}{ccc}\hfill \mathrm{optimize}& & \mathcal{F}\left[\left\{{\Pi }_l\right\}\right],\hfill \end{array}$$

(6a)

$$\begin{array}{ccc}\hfill \mathrm{s}.\mathrm{t}.& & \left\{{\Pi }_l\right\}\in \Omega .\hfill \end{array}$$

(6b)

A simple gradient method is very likely to take $\left\{{\Pi }_l\right\}$ outside of the quantum measurement space; for this, we employ Gilbert’s algorithm to guarantee the condition in Equation (4). In addition, we rewrite the POVM as $\left\{{\Pi }_l\right\}=\left\{{\Pi }_1,{\Pi }_2,\cdots ,{\Pi }_{L-1},\mathbb{I}-{\sum }_{l=1}^{L-1}{\Pi }_l\right\}$ to satisfy the condition in Equation (5). Then, the structure of optimization proceeds as follows.

Taking the to-be-minimized objective function as an example, for the $(k+1)$th iteration, first update the $(L-1)$ elements foremost of the measurement operators with the DG scheme to obtain

$$\begin{array}{cc}\hfill {\Pi }_{l,k+1}& ={\Pi }_{l,k}-ϵG\left({\Pi }_{l,k}\right)\hfill \\ \hfill & \equiv DG\left[{\Pi }_{l,k},G\left({\Pi }_{l,k}\right),ϵ\right].\hfill \end{array}$$

(7)

Here, $ϵ$ represents the step size of the update which can be any positive value, and k is the number of iterations. Second, normalize the measurement operators ${\Pi }_{l,k+1}$ as density matrices ${\rho }_{l,k+1}$, such that

$${\rho }_{l,k+1}=\frac{{\Pi }_{l,k+1}}{\mathrm{tr}\left({\Pi }_{l,k+1}\right)},$$

(8)

which could be nonphysical. Third, use Gilbert’s algorithm to project ${\rho }_{l,k+1}$ back to the quantum state space $\mathcal{Q}$, i.e., ${\rho }_{l,k+1}\to {\rho }_{l,k+1}^{\mathcal{Q}}=\mathcal{S}\left({\rho }_{l,k+1}\right)$. Finally, reconstruct the physical measurement operators as

$${\left\{{\Pi }_{l,k+1}^{\Omega }={\rho }_{l,k+1}^{\mathcal{Q}}{t}_{l,k+1}\right\}}_{l=1}^{L-1},$$

(9)

$$\begin{array}{cc}\hfill {\Pi }_{L,k+1}^{\Omega }& =\mathbb{I}-\sum _{l=1}^{L-1}{\Pi }_{l,k+1}^{\Omega },\hfill \end{array}$$

(10)

where the parameter $t_l$ is obtained by fixing the obtained ${\rho }_{l,k+1}^{\mathcal{Q}}$ to obtain ${\left\{t_{l,k+1}\right\}}_{l=1}^{L-1}=\mathrm{argmin}\mathcal{F}\left[{\left\{t_{l,k+1}\right\}}_{l=1}^{L-1}\right]$. Here, to ensure that the first $(L-1)$ measurement operators satisfy condition Equation (4), only $t_{l,k+1}\ge 0$ is required since ${\rho }_{l,k+1}^{\mathcal{Q}}\ge 0$ is guaranteed by using Gilbert’s algorithm. Meanwhile, in order to ensure that the last element of the new POVM satisfies the condition in Equation (4), let

$${\Pi }_{L,k+1}^{\Omega }=\mathbb{I}-\sum _{l=1}^{L-1}\left({\rho }_{l,k+1}^{\mathcal{Q}}{t}_{l,k+1}\right)\ge 0.$$

(11)

Hence, we obtain the new POVM $\left\{{\Pi }_{k+1,l}^{\Omega }\right\}$ that satisfies the condition in Equation (6b) after each iteration. Whenever the difference between the values of the adjacent iterations is less than a certain threshold, the iteration stops, and the optimal POVM is obtained. Otherwise, continue with the iteration and the step size is controlled by a step factor $\beta$. When ${\mathcal{F}}_k<{\mathcal{F}}_{k-1}$, the step size is appropriately selected. When ${\mathcal{F}}_k>{\mathcal{F}}_{k-1}$, it indicates that the step size selection is too large, and the step factor $\beta$ needs to be used to adjust the step size. See the DG algorithm in Algorithm 1.

However, the DG algorithm has some disadvantages, such as slow optimization speed and low accuracy. For faster convergence, one can choose the APG algorithm [5,24]. The APG algorithm adjusts the direction of the gradient at each step, which improves the convergence speed of the algorithm. In simple terms, the APG algorithm has introduced a companion operator $E_{l,k}={\Pi }_{l,k}+\frac{{\theta }_{k-1}-1}{{\theta }_k}\left({\Pi }_{l,k}-{\Pi }_{l,k-1}\right)$, which provides the momentum of the previous step controlled by the parameter $\theta$, in order to update the measurement operators ${\Pi }_{l,k}=E_{l,k-1}-ϵG\left(E_{l,k-1}\right)$. See the specific process shown in Algorithm 2.    

Algorithm 1: DG algorithm

Algorithm 2: APG algorithm

3. Applications

In this section, we demonstrate the efficacy of our algorithms by optimizing arbitrary convex as well as nonconvex functions over the space of quantum measurements.

3.1. Convex Functions

In quantum measurement tomography [25,26,27], a set of known probe states ${\rho }_m$ is measured to provide the information needed to reconstruct an unknown POVM $\left\{{\Pi }_l\right\}$. The probability that the device would respond to the quantum state ${\rho }_m$ by producing the outcome ${\Pi }_l$ is given by

$$p_{lm}=\mathrm{tr}\left({\rho }_m{\Pi }_l\right).$$

(12)

Typically, the linear inversion method [28] can be used to obtain the ideal POVM, but nonphysical results are likely to be obtained. Then, the maximum likelihood estimation (MLE) [29] is proposed to reconstruct the POVM that satisfies all the conditions. However, MLE fails to return any meaningful results when the target POVM is of low rank, which is quite typical, especially in higher-dimensional spaces. These problems can be avoided by using our algorithms.

To estimate the operators $\left\{{\Pi }_l\right\}$, we maximize the likelihood function

$$\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]=\prod _{l=1}^L\prod _{m=1}^M{\left[\mathrm{tr}\left({\rho }_m{\Pi }_l\right)\right]}^{f_{lm}},$$

(13)

where M is the number of different input states ${\rho }_m$, and

$$f_{lm}=\frac{n_{lm}}{n},$$

(14)

with $n_{lm}$ denoting the number of lth outcome when measuring the mth state ${\rho }_m$, and n representing the total number of measured input states. One can see that $\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]$ is not strictly concave, while the log-likelihood $ln\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]$ is. Here, we minimize the negative log-likelihood function $\mathcal{F}\left[\left\{{\Pi }_l\right\}\right]=-ln\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]$ with

$$\begin{array}{c}\hfill ln\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]=\sum _{l=1}^L\sum _{m=1}^M{f}_{lm}ln{p}_{lm}.\end{array}$$

(15)

To satisfy the condition in Equation (5), rewrite the objective function as

$$ln\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]=\sum _{l=1}^{L-1}\sum _{m=1}^M{f}_{lm}ln\left[\mathrm{tr}\left({\rho }_m{\Pi }_l\right)\right]+\sum _{m=1}^M{f}_{Lm}ln\left\{\mathrm{tr}\left[{\rho }_m\left(\mathbb{I}-\sum _{l=1}^{L-1}{\Pi }_l\right)\right]\right\}.$$

(16)

The gradient of $ln\mathcal{L}\left[\left\{{\Pi }_l\right\}\right]$ with respect to ${\Pi }_l$ is

$$\begin{array}{cc}\hfill \nabla ln\mathcal{L}\left({\Pi }_l\right)=& \sum _{m=1}^M\left[\left(\frac{f_{lm}}{p_{lm}}-\frac{f_{lm}}{1-{\sum }_{l=1}^{L-1}{p}_{lm}}\right){\rho }_m\right].\hfill \end{array}$$

(17)

For numerical simulations, we mainly consider Pauli measurements which are the most commonly-used measurements in quantum information processing. Then, the cases of one qubit, one qutrit, two qubits, and two qutrits are used for the experimental setup, respectively. Specifically, the setups of these four scenarios are described below.

3.1.1. One Qubit

For one qubit, we take the eigenstates of ${\sigma }_z$ and the superposition states $-\frac{1}{\sqrt{2}}\left(|0_z\rangle \pm |1_z\rangle \right)$ and $\frac{1}{\sqrt{2}}\left(|0_z\rangle \pm i|1_z\rangle \right)$ as the input states. In the measurement setup, we select the projection of the spin along the x-axis, i.e.,

$${\Pi }_1=|0_x\rangle \langle 0_x|;{\Pi }_2=|1_x\rangle \langle 1_x|.$$

(18)

3.1.2. One Qutrit

For one qutrit, we use 12 different input states: three eigenstates of ${\sigma }_z$, $|-1_z\rangle$, $|0_z\rangle$ and $|1_z\rangle$, and nine superposition states $\frac{1}{\sqrt{2}}\left(|-1_z\rangle +e^{i{\psi }_j}|0_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|0_z\rangle +e^{i{\psi }_j}|1_z\rangle \right)$ and $\frac{1}{\sqrt{2}}\left(|-1_z\rangle +e^{i{\Psi }_j}|1_z\rangle \right)$, where $j=1,2,3$; and ${\psi }_1=0$, ${\psi }_2=\frac{\pi }{2}$, and ${\psi }_3=\pi$. The device measures the projection of the spin along the x-axis, and the POVM are projectors

$${\Pi }_1=|-1_x\rangle \langle -1_x|;{\Pi }_2=|0_x\rangle \langle 0_x|;{\Pi }_3=|1_x\rangle \langle 1_x|.$$

(19)

3.1.3. Two Qubits

In the case of two qubits, we take the tensor products of the four eigenstates of two Pauli-Z operators $|0_z{0}_z\rangle$, $|1_z{1}_z\rangle$, $|0_z{1}_z\rangle$, $|1_z{0}_z\rangle$ and the superposition states $\frac{1}{\sqrt{2}}\left(|0_z{0}_z\rangle +e^{i{\psi }_j}|0_z{1}_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|0_z{0}_z\rangle +e^{i{\psi }_j}|1_z{0}_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|0_z{0}_z\rangle +e^{i{\psi }_j}|1_z{1}_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|0_z{1}_z\rangle +e^{i{\psi }_j}|1_z{0}_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|0_z{1}_z\rangle +e^{i{\psi }_j}|1_z{1}_z\rangle \right)$, $\frac{1}{\sqrt{2}}\left(|1_z{0}_z\rangle +e^{i{\psi }_j}|1_z{1}_z\rangle \right)$ as the probe states, where $j=1,2,3$; ${\psi }_1=0$, ${\psi }_2=\frac{\pi }{2}$, and ${\psi }_3=\pi$. Then, we choose the following POVM for the experimental simulation:

$$\begin{array}{cc}\hfill & {\Pi }_1=|0_x{0}_x\rangle \langle 0_x{0}_x|;{\Pi }_2=|0_x{1}_x\rangle \langle 0_x{1}_x|;\hfill \\ \hfill & {\Pi }_3=|1_x{0}_x\rangle \langle 1_x{0}_x|;{\Pi }_4=|0_x{0}_x\rangle \langle 0_x{0}_x|.\hfill \end{array}$$

(20)

3.1.4. Two Qutrits

Finally, for the case of two qutrits, we perform a numerical simulation of the Stern–Gerlach apparatus measuring two particles with spin-1. We assume 45 different input states: $|1_z-1_z\rangle$, $|-1_z{0}_z\rangle$, $|-1_z{1}_z\rangle$, $|0_z-1_z\rangle$, $|0_z{0}_z\rangle$, $|0_z{1}_z\rangle$, $|1_z{0}_z\rangle$, $|1_z{1}_z\rangle$, $|-1_z-1_z\rangle$ and 36 superposition states. In the simulation, the device measures the projection of the spin along the x -axis, and the POVM are projectors

$$\begin{array}{cc}\hfill {\Pi }_1& =|0_x{1}_x\rangle \langle 0_x{1}_x|;{\Pi }_2=|0_x-1_x\rangle \langle 0_x-1_x|;\hfill \\ \hfill {\Pi }_3& =|1_x{0}_x\rangle \langle 1_x{0}_x|;{\Pi }_4=|1_x-1_x\rangle \langle 1_x-1_x|;\hfill \\ \hfill {\Pi }_5& =|0_x{0}_x\rangle \langle 0_x{0}_x|;{\Pi }_6=|-1_x{0}_x\rangle \langle -1_x{0}_x|;\hfill \\ \hfill {\Pi }_7& =|1_x{1}_x\rangle \langle 1_x{1}_x|;{\Pi }_8=|-1_x{1}_x\rangle \langle -1_x{1}_x|;\hfill \\ \hfill {\Pi }_9& =|-1_x-1_x\rangle \langle -1_x-1_x|.\hfill \end{array}$$

(21)

For each case of simulation, the number of measurements for each probe state is 300, ${10}^5$, ${10}^5$, and $5\times {10}^5$, respectively. Then, according to the frequency obtained by the simulated data, we use our algorithm to reconstruct the POVM. The fidelity between different POVM elements is defined as the fidelity between the two states $\sigma$ and $\rho$, i.e.,

$$\begin{array}{c}\hfill F(\sigma ,\rho ):={\left(\mathrm{tr}\sqrt{\sqrt{\sigma }\rho \sqrt{\sigma }}\right)}^2=F\left(\frac{{\Pi }_l}{\mathrm{tr}\left({\Pi }_l\right)},\frac{{\Pi }_j}{\mathrm{tr}\left({\Pi }_j\right)}\right).\end{array}$$

(22)

In addition, the overall fidelity between two POVMs ${\left\{{\Pi }_l\right\}}_{l=1}^L$ and ${\left\{{\Pi }_j\right\}}_{j=1}^L$ on a d-dimensional Hilbert space is defined by

$$F({\Pi }_l,{\Pi }_j):={\left[\sum _{l=1}^L{w}_l\sqrt{F\left(\frac{{\Pi }_l}{\mathrm{tr}\left({\Pi }_l\right)},\frac{{\Pi }_j}{\mathrm{tr}\left({\Pi }_j\right)}\right)}\right]}^2,$$

(23)

with $w_l=\frac{\sqrt{\mathrm{tr}\left({\Pi }_l\right)\mathrm{tr}\left({\Pi }_j\right)}}{d}$ [30]. The overall fidelities of the reconstructed POVMs are shown in Figure 1. Figure 2 and Figure 3 present the variations of fidelity of the POVM elements reconstructed using the DG algorithm and APG algorithm with respect to the number of iteration steps in different cases. We can see that these two algorithms are almost identical in accuracy, and the fidelities of the measurement operators are close to 1. Generally speaking, the APG algorithm converges faster than the DG algorithm. In addition, one notices that the fidelity of the last element in some of the simulations is not always increasing, which is a result of the constraint that we set in Equation (11).

3.2. Nonconvex Functions

Quantum detector self-characterization (QDSC) tomography is another method for characterizing quantum measurements. Unlike quantum measurement tomography, this method does not require knowing the specific form of the input probe states, but directly optimizes the cost function based on the measurement statistic ${\mathit{f}}_{\mathit{m}}$ to reconstruct the measurements. For POVM with L outcomes detected by m states, a data set of the measurement statistic $f_{lm}$ is obtained. We write the distribution of the data for each state as a vector

$${\mathit{f}}_{\mathit{m}}=\left(\begin{array}{c}{f}_{1m}\\ f_{2m}\\ ⋮\\ f_{Lm}\end{array}\right).$$

(24)

For the one qubit case, define $N_{i,l}={\mathit{b}}_i^T{\mathit{b}}_l$ and write the POVM as

$${\Pi }_l=a_{l}I+{\mathit{b}}_l·\mathit{\sigma }$$

(25)

under the Bloch representation, where i and l represent the number of rows i and columns l of the matrix N, $\mathit{a}={(a_1\cdots a_L)}^T$, ${\mathit{b}}_l=(b_{l,x},b_{l,y},b_{l,z})$, $\mathit{\sigma }=({\sigma }_x,{\sigma }_y,{\sigma }_z)$, $1\le i,l\le L$. The matrix N and vector $\mathit{a}$ can be represented as

$$\begin{array}{cc}\hfill N_{i,l}& ={\mathit{b}}_i^T{\mathit{b}}_l=\frac{1}{2}\mathrm{tr}\left({\Pi }_i{\Pi }_l\right)-\frac{1}{4}\mathrm{tr}\left({\Pi }_i\right)\mathrm{tr}\left({\Pi }_l\right),\hfill \end{array}$$

(26)

$$a_l=\frac{1}{2}\mathrm{tr}\left({\Pi }_l\right).$$

(27)

Then, optimization of the cost function $\mathcal{F}\left(N^{+},\mathit{a}\right)$ is given by [23]

$$\begin{array}{ccc}\hfill \min & & \sum _m{\left[1-{\left({\mathit{f}}_{\mathit{m}}-\mathit{a}\right)}^T{N}^{+}\left({\mathit{f}}_{\mathit{m}}-\mathit{a}\right)\right]}^2,\hfill \end{array}$$

(28a)

$$\begin{array}{ccc}\hfill \mathrm{s}.\mathrm{t}.& & a_l^2-N_{l,l}>=0,\hfill \end{array}$$

(28b)

where $N^{+}$ stands for the Moore–Penrose pseudoinverse of N. One notices that the objective function is nonconvex. Optimization of nonconvex functions is difficult as local minima might be found. Interestingly, we find that our algorithm can also be used to optimize nonconvex functions. Since our algorithm guarantees the conditions for quantum measurements, one only needs to optimize the objective function regardless of the constraint in Equation (28b).

For numerical simulations, we choose 50 probe states:

$$\frac{1}{2}\left(\mathbb{I}+{\sigma }_z\right),\frac{1}{2}\left(\mathbb{I}-{\sigma }_z\right),\frac{1}{2}\left(\mathbb{I}+sin\frac{i\pi }{4}cos\frac{n\pi }{8}{\sigma }_x+sin\frac{i\pi }{4}sin\frac{n\pi }{8}{\sigma }_y+cos\frac{i\pi }{4}{\sigma }_z\right),$$

(29)

where $i=1,2,\cdots ,6$; $n=1,2,\cdots ,8$. In addition, we use the two-dimensional SIC POVM as the measurement device, and each state is measured 200 times. The APG algorithm is used to optimize the objective function. First, select any set of POVM operators in the measurement space, and use Equations (26) and (27) to obtain the initial values $N_k^{+}$ and ${\mathit{a}}_k$, respectively. Similarly, we calculate the gradient of the objective function in Equation (28a). The gradient of the objective function is given by

$$\delta \mathcal{F}\left(\mathit{a}\right)=\sum _{m}2{\left(1-{\mathit{f}}_{\mathit{m}}-\mathit{a}\right)}^T{N}^{+}\left({\mathit{f}}_{\mathit{m}}-\mathit{a}\right)\left\{{\left(N^{+}\right)}^T{\mathit{f}}_{\mathit{m}}+N^{+}{\mathit{f}}_{\mathit{m}}-\left[N^{+}+{\left(N^{+}\right)}^T\right]\mathit{a}\right\},$$

(30)

$$\delta \mathcal{F}\left(N^{+}\right)=\sum _m-2{\left(1-{\mathit{f}}_{\mathit{m}}-\mathit{a}\right)}^T{N}^{+}{\left({\mathit{f}}_{\mathit{m}}-\mathit{a}\right)}^2{\left({\mathit{f}}_{\mathit{m}}-\mathit{a}\right)}^T.$$

(31)

The values of $N_{k+1}$ and ${\mathit{a}}_{\mathit{k}+\mathbf{1}}$ are obtained by iterating over $N_k$ and ${\mathit{a}}_{\mathit{k}}$ using gradient descent; then, ${\mathit{b}}_{l,k+1}$ is obtained by decomposing $N_{k+1}$. In the experiment, we specify that the reference frame, i.e., the vector ${\mathit{b}}_1$ is parallel to the z-direction of the Bloch sphere, and set the $xz$ plane of the Bloch sphere as the plane determined by the vectors ${\mathit{b}}_1$ and ${\mathit{b}}_2$. This is equivalent to $b_{1,x}=b_{1,y}=b_{2,y}=0$. Then, ${\left\{{\Pi }_{l,k+1}\right\}}_{l=1}^{L-1}$ can be obtained by using Equation (25), which is the update for ${\left\{{\Pi }_{l,k}\right\}}_{l=1}^{L-1}$.

The fidelity of each POVM element can approach 1 in a very small number of iteration steps; see Figure 4. Then, the fidelities of the measurements are compared with the ones reported in [23], demonstrating that the performance of our algorithm is slightly better; see Figure 5.

4. Conclusions

We have proposed two reliable algorithms for optimizing arbitrary functions over the quantum measurement space. For a demonstration, we have shown several examples on the convex function of quantum measurement tomography with different dimensions as well as a nonconvex function of one qubit in quantum detector self-characterization tomography. Surprisingly, our method does not encounter the problem of rank deficiency. Compared with SDP, our method can be easily applied to higher-dimensional cases as well as to optimize nonconvex functions. Moreover, our method reports better results as compared to previous approaches. For future work, we will consider the optimization over the joint space of quantum states and quantum measurements, for tasks such as calculating the capacity of quantum channels.