text:
# Maximum-likelihood estimation of quantum measurement

## I Introduction

Let us imagine that we possess an apparatus which performs some measurement on certain quantum mechanical system such as spin of electron. We do not know which measurement is associated with the device and we would like to find it out.

Obviously, the path to follow here is to perform set of measurements on various known quantum states and then estimate the unknown measurement from the collected data. Such estimation strategy belongs to the broad class of the quantum reconstruction procedures which attracted considerable attention recently. The quantum state reconstruction has been widely studied and now represents a well established technique in many branches of quantum physics (for a review, see, e.g. ). The estimation of the quantum mechanical processes describing input-output transformations of quantum devices has been discussed in  and the problem of complete characterization of arbitrary measurement process has been recently addressed in .

Suppose that the apparatus can respond with $`k`$ different measurement outcomes. As is well known from the theory of quantum measurement , such device is completely characterized by $`k`$ positive operator valued measures (POVM) $`\widehat{\mathrm{\Pi }}_l`$ which govern the measurement statistics,

$$p_{lm}=\mathrm{Tr}[\widehat{\mathrm{\Pi }}_l\widehat{\varrho }_m],$$
(1)
where $`\widehat{\varrho }_m`$ denotes density matrix of the quantum state subject to the measurement, $`p_{lm}`$ denotes the probability that the apparatus would respond with outcome $`\widehat{\mathrm{\Pi }}_l`$ to the quantum state $`\widehat{\varrho }_m`$, and Tr stands for the trace. The POVMs are positive semi-definite Hermitian operators,

$$\widehat{\mathrm{\Pi }}_l0,$$
(2)
which decompose the unit operator,

$$\underset{l=1}{\overset{k}{}}\widehat{\mathrm{\Pi }}_l=\widehat{I}.$$
(3)
The condition (2) ensures that $`p_{lm}0`$ and (3) follows from the requirement that total probability is normalized to unity, $`_{l=1}^kp_{lm}=1`$.

In order to determine the POVMs, one should measure on different known quantum states $`\widehat{\varrho }_m`$ and then estimate the POVMs $`\widehat{\mathrm{\Pi }}_l`$ from the acquired statistics. Let $`F_{lm}`$ denote the total number of detections of $`\widehat{\mathrm{\Pi }}_l`$ for the measurements performed on the quantum state $`\widehat{\varrho }_m`$. Assuming that the theoretical detection probability $`p_{lm}`$ (1) can be replaced with relative frequency, we may write

$$\mathrm{Tr}[\widehat{\mathrm{\Pi }}_l\widehat{\varrho }_m]\underset{i,j=1}{\overset{N}{}}\mathrm{\Pi }_{l,ij}\varrho _{m,ji}=\frac{F_{lm}}{_{l^{}=1}^kF_{l^{}m}},$$
(4)
where $`N`$ is the dimension of the Hilbert space on which the operators $`\widehat{\mathrm{\Pi }}_l`$ act. This establishes a system of linear equations for the unknown elements of the operators $`\mathrm{\Pi }_{l,ij}`$, which may easily be solved if sufficient amount of data is available. This approach is a direct analogue of linear reconstruction algorithms devised for quantum state reconstruction. The linear inversion is simple and straightforward, but it has also one significant disadvantage. The linear procedure cannot guarantee the required properties of $`\widehat{\mathrm{\Pi }}_l`$, namely the conditions (2). Consequently, the linear estimation may lead to unphysical POVMs, predicting negative probabilities $`p_{lm}`$ for certain input quantum states. To avoid such problems, one should resort to more sophisticated nonlinear reconstruction strategy.

In this paper we show that the Maximum-likelihood (ML) estimation is suitable and can be successfully used for the calibration of the measuring apparatus. ML estimation has been recently applied to reconstruction of quantum states  and quantum processes (complete positive maps between density matrices) . Here we employ it to reconstruct an unknown quantum measurement thereby demonstrating again the remarkable versatility and usefulness of ML estimation. The generic formalism is developed in Sec. II and an illustration of our method is provided in Sec. III, where we estimate a measurement performed by Stern-Gerlach apparatus.

## II Maximum-likelihood estimation

The estimated operators $`\widehat{\mathrm{\Pi }}_l`$ maximize the likelihood functional

$$[\{\widehat{\mathrm{\Pi }}_l\}]=\underset{l=1}{\overset{k}{}}\underset{m=1}{\overset{M}{}}\left(\mathrm{Tr}[\widehat{\mathrm{\Pi }}_l\widehat{\varrho }_m]\right)^{f_{lm}},$$
(5)
where

$$f_{lm}=F_{lm}\left[\underset{l^{}=1}{\overset{k}{}}\underset{m^{}=1}{\overset{M}{}}F_{l^{}m^{}}\right]^1$$
(6)
is the relative frequency and $`M`$ is the number of different quantum states $`\widehat{\varrho }_m`$ used for the reconstruction. The maximum of the likelihood functional (5) has to be found in the sub-space of physically allowed operators $`\widehat{\mathrm{\Pi }}_l`$. We can decompose each operator $`\widehat{\mathrm{\Pi }}_l`$ as

$$\widehat{\mathrm{\Pi }}_l=\underset{q=1}{\overset{N}{}}r_{lq}|\varphi _{lq}\varphi _{lq}|,$$
(7)
where $`r_{lq}0`$ are the eigenvalues of $`\widehat{\mathrm{\Pi }}_l`$ and $`|\varphi _{lq}`$ are corresponding orthonormal eigenstates. The maximum of $`[\{\widehat{\mathrm{\Pi }}_l\}]`$ can be found from the extremum conditions. It is convenient to work with the logarithm of the original likelihood functional and the constraint (3) has to be incorporated by introducing a Hermitian matrix of Lagrange multipliers $`\lambda _{ij}=\lambda _{ji}^{}`$. The extremum conditions then read

$`{\displaystyle \frac{}{\varphi _{lq}|}}[{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \underset{m=1}{\overset{M}{}}}f_{lm}\mathrm{ln}\left({\displaystyle \underset{q=1}{\overset{N}{}}}r_{lq}\varphi _{lq}|\widehat{\varrho }_m|\varphi _{lq}\right)`$ (8)
$`{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \underset{q=1}{\overset{N}{}}}r_{lq}\varphi _{lq}|\widehat{\lambda }|\varphi _{lq}]=0.`$ (9)
Thus we immediately find

$$r_{lq}|\varphi _{lq}=\widehat{R}_lr_{lq}|\varphi _{lq},$$
(10)
where

$$\widehat{R}_l=\widehat{\lambda }^1\underset{m=1}{\overset{M}{}}\frac{f_{lm}}{p_{lm}}\widehat{\varrho }_m$$
(11)
and

$$\widehat{\lambda }=\underset{i,j=1}{\overset{N}{}}\lambda _{ij}|ij|.$$
(12)
Let us now multiply (10) by $`\varphi _{lq}|`$ from the right and sum over $`q`$. Thus we obtain

$$\widehat{\mathrm{\Pi }}_l=\widehat{R}_l\widehat{\mathrm{\Pi }}_l.$$
(13)
On averaging (13) and its Hermitean conjugate counterpart, we get

$$\widehat{\mathrm{\Pi }}_l=\frac{1}{2}(\widehat{R}_l\widehat{\mathrm{\Pi }}_l+\widehat{\mathrm{\Pi }}_l\widehat{R}_l^{}).$$
(14)
The matrix of Lagrange multipliers $`\widehat{\lambda }`$ should be determined from the constraint (3). On summing Eq. (13) over $`l`$, we find

$$\widehat{\lambda }=\underset{l=1}{\overset{k}{}}\underset{m=1}{\overset{M}{}}\frac{f_{lm}}{p_{lm}}\widehat{\varrho }_m\widehat{\mathrm{\Pi }}_l.$$
(15)
Eqs. (14) and (15) can be conveniently solved by means of repeated iterations.

If the linear inversion based on Eqs. (4) provides physically sensible result, then the ML estimate agrees with this linear reconstruction. To prove it explicitly, let us assume that the set of POVMs $`\widehat{\mathrm{\Pi }}_l`$ solves the Eqs. (4). Thus we have

$$p_{lm}=\frac{f_{lm}}{_{l^{}=1}^kf_{l^{}m}}.$$
(16)
On inserting this expression into (13), we obtain

$$\widehat{\mathrm{\Pi }}_l=\underset{m}{}^{}\left(\underset{l^{}=1}{\overset{k}{}}f_{l^{}m}\right)\widehat{\lambda }^1\widehat{\varrho }_m\widehat{\mathrm{\Pi }}_l.$$
(17)
Here the prime indicates that we should sum only over those $`m`$ with nonzero $`f_{lm}`$. However, this restriction may be dropped. If $`f_{lm}=p_{lm}=0`$, then $`\widehat{\varrho }_m\widehat{\mathrm{\Pi }}_l=0`$ and the addition of zero to the right-hand side of (17) changes nothing. Thus the set of $`k`$ equations (17) reduces to

$$\widehat{\lambda }=\underset{m=1}{\overset{M}{}}\underset{l=1}{\overset{k}{}}f_{lm}\widehat{\varrho }_m$$
(18)
which is the formula for the operator of Lagrange multipliers valid when the measured data are compatible with some set of POVMs. Notice that $`\widehat{\lambda }`$ is positive definite.

The differences between linear reconstructions and ML estimation occur if the experimental data are not compatible with any physically allowed set of POVMs. The procedure of ML estimation may be interpreted as a synthesis of information from mutually incompatible observations. The ML can correctly handle noisy data and provides reliable estimates in cases when linear algorithms fail.

Notice that the operators $`\widehat{R}_l`$ contain the inversion of the matrix $`\widehat{\lambda }`$. The reconstruction is possible only on such subspace of the total Hilbert space where the inversion $`\widehat{\lambda }^1`$ exists. This restriction can easily be understood if we make use of Eq. (18). The experimental data contain only information on the Hilbert subspace probed by the density matrices $`\widehat{\varrho }_m`$ and the reconstruction of the POVMs must be restricted to this subspace.

One could complain that it is not certain that the positive definiteness of POVMs $`\widehat{\mathrm{\Pi }}_l`$ is preserved during iterations based on Eqs. (14) and (15). We can, however, avoid such complaints by devising an iterative algorithm which exactly satisfies the constraints (2) and (3) at each iteration step. We observe that we can formally decompose the POVMs as

$$\widehat{\mathrm{\Pi }}_l=\widehat{D}_l^{}\widehat{D}_l,$$
(19)
where

$$\widehat{D}_l=\underset{q=1}{\overset{N}{}}\sqrt{r_{ql}}|q\varphi _{lq}|,$$
(20)
and $`|q`$ is some chosen orthonormal basis, $`q|q^{}=\delta _{qq^{}}`$.

From the equations (10) we can derive

$$\widehat{D}_l=\widehat{D}_l\widehat{R}_l^{}.$$
(21)
The constraint (3) provides formula for the operator of Lagrange multipliers,

$$\widehat{\lambda }^1\widehat{G}\widehat{\lambda }^1=\widehat{I},$$
(22)
where $`\widehat{G}`$ is positive operator

$$\widehat{G}=\underset{l=1}{\overset{k}{}}\underset{m,m^{}=1}{\overset{M}{}}\frac{f_{lm}}{p_{lm}}\frac{f_{lm^{}}}{p_{lm^{}}}\widehat{\varrho }_m^{}\widehat{\mathrm{\Pi }}_l\widehat{\varrho }_m.$$
(23)
Upon solving (22) we get $`\widehat{\lambda }=\widehat{G}^{1/2}`$. We fix the branch of the square root of $`\widehat{G}`$ by requiring that $`\widehat{\lambda }`$ should be positive definite operator. We can factorize the matrix $`\widehat{G}`$ as $`\widehat{G}=\widehat{U}^{}\widehat{\mathrm{\Lambda }}\widehat{U}`$ where $`\widehat{U}`$ is unitary matrix and $`\widehat{\mathrm{\Lambda }}`$ is diagonal matrix containing eigenvalues of $`\widehat{G}`$. We define

$$\widehat{\mathrm{\Lambda }}^{1/2}=\mathrm{diag}(\sqrt{\mathrm{\Lambda }_{11}},\mathrm{},\sqrt{\mathrm{\Lambda }_{NN}})$$
(24)
and we can write

$$\widehat{\lambda }=\widehat{U}^{}\widehat{\mathrm{\Lambda }}^{1/2}\widehat{U}.$$
(25)
The advantage of the iterative procedure based on (21) and (25) is that both conditions (2) and (3) are exactly fulfilled at each iteration step. The disadvantage of this approach is the greater numerical complexity in comparison to iterations based on (13) and (15), because we must calculate the eigenvalues of the matrix $`\widehat{G}`$ at each iteration step.

The determination of the quantum measurement can simplify considerably if we have some a-priori information about the apparatus. For example, if we know that we deal with a photodetector, then we have to estimate only a single parameter, the absolute photodetection effeciency $`\eta `$ . Here we briefly consider a broader class of phase-insensitive detectors which are sensitive only to the number of photons in a single mode of electromagnetic field. The POVMs describing phase-insensitive detector are all diagonal in the Fock basis,

$`\widehat{\mathrm{\Pi }}_l={\displaystyle \underset{n}{}}r_{ln}|nn|`$

and the ML estimation reduces to the determination of the eigenvalues $`r_{ln}0`$. The extremum Eqs. (14) and (15) simplify to

$`r_{ln}`$ $`=`$ $`{\displaystyle \frac{r_{ln}}{\lambda _n}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{f_{lm}}{p_{lm}}}\varrho _{m,nn},`$ (26)
$`\lambda _n`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{l}{}}{\displaystyle \frac{f_{lm}}{p_{lm}}}\varrho _{m,nn}r_{ln},`$ (27)
$`p_{lm}`$ $`=`$ $`{\displaystyle \underset{n}{}}\varrho _{m,nn}r_{ln}.`$ (28)
Instead of solving the extremum equations, one may directly search for the maximum of $`[\{\widehat{\mathrm{\Pi }}_l\}]`$ with the help of downhill-simplex algorithm . To implement this algorithm successfully, it is necessary to use a minimal parametrization. If we deal with $`N`$ level system, then each $`\widehat{\mathrm{\Pi }}_l`$ is parametrized by $`N^2`$ real numbers. Since the constraint (3) allows us to determine one POVM in terms of remaining $`k1`$ ones, the number of independent real parameters reads $`N^2(k1)`$. Furthermore we may take advantage of the Cholesky decomposition,

$$\widehat{\mathrm{\Pi }}_l=\widehat{C}_l^{}\widehat{C}_l,$$
(29)
where $`\widehat{C}_l`$ is lower triangular matrix with real elements on its main diagonal. The parametrization (29) is used for the first $`k1`$ operators, and the last one is calculated form (3),

$$\widehat{\mathrm{\Pi }}_k=\widehat{I}\underset{l=1}{\overset{k1}{}}\widehat{C}_l^{}\widehat{C}_l,$$
(30)
thus achieving minimal parametrization. For each parameter set, where $`[\{\widehat{\mathrm{\Pi }}_l\}]`$ is evaluated, one has to check whether (30) is positive semi-definite. If this does not hold, then one sets $`[\{\widehat{\mathrm{\Pi }}_l\}]=0`$, thus restricting the numerical search of the maximum to the domain of physically allowed operators. This domain is a finite volume subspace of a $`N^2(k1)`$ dimensional space.

## III Stern-Gerlach apparatus

In this section we illustrate the developed formalism by means of numerical simulations for Stern-Gerlach apparatus measuring a spin-$`1`$ particle. We choose the three eigenstates of the operator of $`z`$-component of the spin as the basis states, $`\widehat{\sigma }_z|s_z=s_z|s_z`$, $`s_z=1,0,1`$. In this basis, the matrix representation of the spin operators reads

$`\widehat{\sigma }_x={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),\widehat{\sigma }_y={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}0& i& 0\\ i& 0& i\\ 0& i& 0\end{array}\right),`$

$`\widehat{\sigma }_z=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right).`$

In our numerical simulations, we assume $`12`$ different pure quantum states: three eigenstates of $`\widehat{\sigma }_z`$, $`|1_z`$, $`|0_z`$, and $`|1_z`$, and nine superposition states $`2^{1/2}(|1_z+e^{i\psi _j}|0_z),`$ $`2^{1/2}(|0_z+e^{i\psi _j}|1_z),`$ and $`2^{1/2}(|1_z+e^{i\psi _j}|1_z),`$ where $`j=1,2,3`$ and $`\psi _1=0`$, $`\psi _2=\pi /2`$, and $`\psi _3=\pi `$. The measurement on each state is performed $`𝒩`$ times thus the total number of measurements is $`12𝒩`$. We have performed Monte Carlo simulations of the measurement and then we have reconstructed the three POVMs $`\widehat{\mathrm{\Pi }}_l`$ characterizing given apparatus. In the simulations, the device measures the projection of the spin along $`x`$ axis, and the true POVMs are projectors

$$\widehat{\mathrm{\Pi }}_1=|1_x1_x|,\widehat{\mathrm{\Pi }}_2=|0_x0_x|,\widehat{\mathrm{\Pi }}_3=|1_x1_x|,$$
(31)
$`\widehat{\sigma }_x|s_x=s_x|s_x`$. The elements of the POVMs can be conveniently collected in a real vector

$`\stackrel{}{\mathrm{\Pi }}_l`$ $`=`$ $`(\mathrm{\Pi }_{l,11},\mathrm{\Pi }_{l,22},\mathrm{\Pi }_{l,33},\mathrm{Re}\mathrm{\Pi }_{l,12},\mathrm{Im}\mathrm{\Pi }_{l,12},`$ (33)
$`\mathrm{Re}\mathrm{\Pi }_{l,13},\mathrm{Im}\mathrm{\Pi }_{l,13},\mathrm{Re}\mathrm{\Pi }_{l,23},\mathrm{Im}\mathrm{\Pi }_{l,23}).`$
The estimated POVMs as well as the ’true’ POVMs (31) are plotted in Fig. 1. The estimate has been obtained by iteratively solving the extremum equations and both sets (14), (15) and (21), (25) provide identical results. The reconstructed operators are in good agreement with the exact ones. Equally important is the fact that the estimated operators $`\widehat{\mathrm{\Pi }}_l`$ meet the constraints (2) and (3).

In summary, we have shown how to reconstruct a generic quantum measurement with the use of Maximum-likelihood principle. Our method guarantees that the estimated POVMs, which fully describe the measuring apparatus, meet all the required positivity and completness constraints. The numerical feasibility of our technique has been illustrated by means of numerical simulations for Stern-Gerlach apparatus.

###### Acknowledgements.

I would like to thank Z. Hradil, M. Ježek, and J. Řeháček for stimulating discussions. This work was supported by Research Project CEZ: J14/98: 153100009 “Wave and Particle Optics” of the Czech Ministry of Education.