Spectral thresholding quantum tomography for low rank states

We introduce here the rank-penalized nonlinear estimator, which can be computed from the LSE by truncation to an appropriately chosen rank.

As noted earlier, while the LSE is unbiased, it has a large variance due to the fact that it does not take into account the physical constraints encoded in the unknown parameter ρ. A possible remedy is to 'project' the LSE onto the space of physical states, i.e. positive, trace-one matrices. This method will be discussed in the following subsection. Another improvement can be obtained by taking into account the 'sparsity' properties of the unknown state. For instance, in many experimental situations the goal is to create a particular low rank, or even pure state. The fact that such states can be characterized with a smaller number of parameters than a general density matrix, has two important consequences. Firstly, they can be estimated by measuring an 'informationally incomplete' set of observables, as demonstrated in [25, 26]. Secondly, the prior information can be used to design estimators with reduced estimation error compared with generic methods which do not take into account the structure of the state. Roughly speaking, this is because each unknown parameter brings its own contribution to the overall error of the estimator.

However, the downside of working with a lower dimensional model is that it contains built-in assumptions which may not be satisfied by the true (unknown) physical state. Preparing a pure state is strictly speaking rarely achievable due to various experimental imperfections, so using a pure state statistical model is in fact an oversimplification and can lead to erroneous conclusions about the true state. On the other hand, one can argue that when the (small) experimental noises are taken into account, the actual state is 'effectively' low rank, i.e. it has a small number of significant eigenvalues and a large number of eigenvalues which are so close to zero that they cannot be distinguished from it. Then, the interesting question is how to decide on where to make the cut-off between statistically relevant eigenvalues and pure statistical noise. This is a common problem in statistics which is closely related to that of model selection [30]. Below we describe the rank-penalized estimator addressing this problem, and show that its theoretical and practical performance is superior to the LSE, and is close to what one would expect from an optimal estimator. In addition, its computation requires only the diagonalization of the LSE.

Before presenting a simple algorithm for computing the estimator, we briefly discuss the idea behind its definition. Let

$$\begin{eqnarray}&&{\widehat{\rho }}_{n}^{(\mathrm{ls})}=\displaystyle \sum _{i=1}^{d}{\widehat{\lambda }}_{i}| {\hat{\psi }}_{i}\rangle \langle {\hat{\psi }}_{i}| .\end{eqnarray} \tag{ 4.1 }$$

Be the spectral decomposition of the LSE, with eigenvalues ordered such that $| {\widehat{\lambda }}_{1}| \geqslant ...\geqslant | {\widehat{\lambda }}_{d}| .$ For each given rank $\kappa \in \{1,...,d\}$ we can project ${\widehat{\rho }}_{n}^{(\mathrm{ls})}$ onto the space of matrices of rank κ by computing the matrix which is the closest to ${\widehat{\rho }}_{n}^{(\mathrm{ls})}$ with respect to the norm-two distance

$$\begin{eqnarray*}&&{\widehat{\rho }}_{n}(\kappa ):= \underset{\tau :{\rm{rank}}(\tau )=\kappa }{\mathrm{arg}\;\mathrm{min}}\parallel \tau -{\widehat{\rho }}_{n}^{(\mathrm{ls})}{\parallel }_{2}^{2}.\end{eqnarray*}$$

Although the projection is not a linear operator, ${\widehat{\rho }}_{n}(\kappa )$ is easy to compute, and is obtained by truncating the spectral decomposition (4.1) to the most significant κ eigenvalues

$$\begin{eqnarray*}&&{\widehat{\rho }}_{n}(\kappa )=\displaystyle \sum _{i=1}^{\kappa }{\widehat{\lambda }}_{i}| {\hat{\psi }}_{i}\rangle \langle {\hat{\psi }}_{i}| .\end{eqnarray*}$$

The question is now how to choose the rank κ in order to obtain a good estimator. In figure 2 we illustrate the dependence of the norm-two square error $e(\kappa ):= \parallel {\widehat{\rho }}_{n}(\kappa )-\rho {\parallel }_{2}^{2}$ on the rank, for a particular dataset generated with a rank 6, 4-ions state. As the rank is increased starting with $\kappa =1$ (pure states), the error decreases steeply as ${\widehat{\rho }}_{n}(\kappa )$ becomes less biased, it reaches a minimum close to the true rank, and increases slowly as added parameters increase the variance of the estimator. However, since the state ρ is unknown, the norm-two error and optimal rank for which the minimum is achieved, are unknown. To go around this, we can estimate the error $e(\kappa )$ from the data by means of e.g. cross-validation, as it will be described in section 6. However, in this section we follow a different path, and we define the rank-penalized estimator [8, 35] as the minimizer over κ of the following expression:

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}(\kappa )-{\widehat{\rho }}_{n}^{(\mathrm{ls})}{\parallel }_{2}^{2}+{\nu }^{2}\cdot \kappa =\displaystyle \sum _{i=\kappa +1}^{d}{\widehat{\lambda }}_{i}^{2}+{\nu }^{2}\cdot \kappa ,\end{eqnarray*}$$

where ν is a constant which will be tuned appropriately. The first term quantifies the fit of the truncated estimator with respect to the LSE, while the second term is a penalty which increases with the complexity of the model, i.e. the rank. The rank penalized estimator ${\widehat{\rho }}_{n}^{(\mathrm{pen})}$ is thus the solution of the simple optimization problem

$$\begin{eqnarray}{\widehat{\rho }}_{n}^{(\mathrm{pen})} & := & {\widehat{\rho }}_{n}(\hat{\kappa }),\qquad {\rm{where}}\qquad \hat{\kappa }:= \underset{\kappa =1,\ldots ,d}{\mathrm{arg}\;\mathrm{min}}\left\{\displaystyle \sum _{i=\kappa +1}^{d}{\widehat{\lambda }}_{i}^{2}+{\nu }^{2}\cdot \kappa \right\}=\mathrm{max}\{\kappa \ :{\widehat{\lambda }}_{\kappa }^{2}\geqslant {\nu }^{2}\}.\end{eqnarray} \tag{ 4.2 }$$

This means that the eigenvalues below a certain noise threshold are set to zero while those above the threshold remain unchanged. The following theorem is our first main result, and shows that the appropriate threshold is given by the upper bound on the operator norm error of the LSE, as established in proposition 1.

Zoom In Zoom Out Reset image size

Theorem 1. Let $\theta \gt 0$ be an arbitrary constant, let $c(\theta ):= 1+2/\theta ,$ and let $\varepsilon \gt 0$ be a small parameter. Then with probability larger than $1-\varepsilon ,$ we have

$$\begin{eqnarray}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{2}^{2}\leqslant \underset{\kappa =1,\ldots ,d}{\mathrm{min}}\left\{{c}^{2}(\theta )\displaystyle \sum _{j\gt \kappa }{\lambda }_{j}^{2}(\rho )+2c(\theta )\nu {(\varepsilon )}^{2}\kappa \right\},\end{eqnarray} \tag{ 4.3 }$$

where ${\widehat{\rho }}_{n}^{(\mathrm{pen})}$ is the penalized estimator defined in (4.2) with threshold $\nu {(\varepsilon )}^{2}$ given by

$$\begin{eqnarray*}&&\nu {(\varepsilon )}^{2}:= \displaystyle \frac{2}{n}{\left(\displaystyle \frac{2}{3}\right)}^{k}\mathrm{log}\left(\displaystyle \frac{{2}^{k+1}}{\varepsilon }\right)=\displaystyle \frac{2d}{N}\mathrm{log}\displaystyle \frac{2d}{\varepsilon },\end{eqnarray*}$$

which is assumed to be o(1) with increasing n and k.

Proof. The upper bound follows directly from proposition 1 combined with the following oracle inequality established in [8],

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{2}^{2}\leqslant \underset{\kappa =1,\ldots ,d}{\mathrm{min}}\left\{{c}^{2}(\theta )\underset{R:{\rm{rank}}(R)=\kappa }{\mathrm{min}}\parallel R-\rho {\parallel }_{2}^{2}+2c(\theta )\cdot \kappa \cdot \nu {(\varepsilon )}^{2},\right\}\end{eqnarray*}$$

which holds true provided that $(1+\theta )\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho {\parallel }^{2}\leqslant \nu {(\varepsilon )}^{2}.$ This event occurs with probability larger than $1-\varepsilon .$

Let us make some explanatory remarks on the above result. Firstly, the bound (4.3) applies to all states ρ, not only 'small' rank ones. Recall that ${{\mathcal{S}}}_{d,r}$ denotes the set of states of rank-r states on ${{\mathbb{C}}}^{d}.$ In the special case when ${\rm{rank}}(\rho )=r\leqslant d,$ the theorem implies that, with probability larger than $1-\varepsilon ,$

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{2}^{2}\leqslant 2c(\theta )\cdot r\cdot \nu {(\varepsilon )}^{2}.\end{eqnarray*}$$

If the rank r is much smaller than d, this bound is a significant improvement to the corresponding upper bound $d\cdot \nu {(\varepsilon )}^{2}$ for the LSE, which can be derived by combining the operator norm bound of proposition 1 and the matrix norms inequality $\parallel \tau {\parallel }_{2}^{2}\leqslant d\parallel \tau {\parallel }^{2},$ for $\tau \in M({{\mathbb{C}}}^{d}).$ Moreover, up to a constant factor the rate $\nu {(\varepsilon )}^{2}$ is equal to $d\cdot r\mathrm{log}(d)/N$ which is essentially the ratio of the number of parameters and total number of measurements. In section 5 we will show that apart from the log factor this rate is also optimal, and cannot be improved even if the rank of the state is known, which indicates that the estimator adapts to the complexity of the true parameter. Furthermore, we stress the fact that the bound (4.3) holds true for growing dimension $d={2}^{k}$ as well as the number of measurements n; the bound remains meaningful as far as $d\mathrm{log}d/N\to 0.$

The second observation is that our procedure selects the true rank consistently. Denote by $\hat{\kappa }$ the rank of the resulting estimator ${\widehat{\rho }}_{n}^{(\mathrm{pen})}.$ Following [8] we can prove that, if there exists some κ such that ${\lambda }_{\kappa }(\rho )\gt (1+\delta )\sqrt{\nu (\varepsilon )}$ and ${\lambda }_{\kappa +1}(\rho )\lt (1-\delta )\sqrt{\nu (\varepsilon )}$ for some $\delta \in (0,1),$ then

$$\begin{eqnarray*}&&{\mathbb{P}}(\hat{\kappa }=\kappa )=1-{\mathbb{P}}(\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel \geqslant \delta \sqrt{\nu (\varepsilon )}).\end{eqnarray*}$$

This stresses the fact that the procedure detects the eigenvalues above a threshold related to the error of the LSE. If the true rank of ρ is r and if k and n are such that ν tends to 0 (which always occurs for fixed number of ions k), then ${\lambda }_{r}\gt \nu$ asymptotically and the probability that $\hat{\kappa }=r$ tends to 1.

We can also project ${\widehat{\rho }}_{n}^{(\mathrm{pen})}$ on the matrices with trace 1, to get

$$\begin{eqnarray}&&{\widehat{\rho }}_{n}^{(\mathrm{pen},n)}=\mathrm{arg}\;\underset{R\in {\mathcal{S}}}{\mathrm{min}}\parallel R-{\widehat{\rho }}_{n}^{(\mathrm{pen})}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 4.4 }$$

where ${\mathcal{S}}$ is the set of all density matrices on ${{\mathbb{C}}}^{d}.$ The following corollary shows that the key properties of the estimator are preserved if we additionally normalize it to trace-one after thresholding.

Corollary 1. Under the notation and assumptions of theorem 1 if ρ is an arbitrary state in ${{\mathcal{S}}}_{d,r}$ and if k and n are such that ${\lambda }_{r}\gt \nu (\varepsilon )$ for some $\varepsilon \in (0,1),$ then, with probability larger than $1-\varepsilon ,$

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2}\leqslant 8c(\theta )\cdot r\cdot \nu {(\varepsilon )}^{2}.\end{eqnarray*}$$

Moreover, there exists an absolute constant $C\gt 0$ such that

$$\begin{eqnarray*}&&\underset{\rho \in {{\mathcal{S}}}_{d,r}}{\mathrm{sup}}{E}_{\rho }\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2}\leqslant C\displaystyle \frac{{rd}}{N}\mathrm{log}\left(\displaystyle \frac{2d}{\varepsilon }\right).\end{eqnarray*}$$

Proof. See appendix A.4.

Although the rank-penalized estimator performs well in terms of its risk, it is not necessarily positive and trace-one and therefore it may not represent a physical state. In this section we propose and analyse the following 'physical estimator'

$$\begin{eqnarray}&&{\widehat{\rho }}_{n}^{(\mathrm{phys})}=\mathrm{arg}\;\underset{\sigma \in {\mathcal{S}}(\nu )}{\mathrm{min}}\parallel \sigma -{\tilde{\rho }}_{n}^{(\mathrm{ls})}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 4.5 }$$

where ${\tilde{\rho }}_{n}^{(\mathrm{ls})}$ is the 'normalized LSE' defined in (3.2), and ${\mathcal{S}}(\nu )$ denotes the set of states at noise level ν

$$\begin{eqnarray*}&&{\mathcal{S}}(\nu )=\left\{\sigma {\rm{}}\;{\rm{density}}\;{\rm{matrix}}\;{\rm{with}}\;{\rm{eigenvalues}}\;{\rm{}}{\lambda }_{j}\in \{0\}\cup (4\nu ,1],j=1,\ldots ,d\right\}.\end{eqnarray*}$$

In particular, the space of all density matrices correspond to $\nu =0$ and is denoted ${\mathcal{S}}.$ The estimator ${\widehat{\rho }}_{n}^{(\mathrm{phys})}$ is therefore the physical state which is closest to the (normalized) LSE, and whose non-zero eigenvalues are above the threshold$4\nu .$

Before analysing the performance of the estimator, we describe its numerical implementation through the following simple iterative algorithm.

Let ${\tilde{\lambda }}_{1}\geqslant ...\geqslant {\tilde{\lambda }}_{d}$ denote the eigenvalues of ${\tilde{\rho }}_{n}^{(\mathrm{ls})}.$ Let ${\ell }=0,$ and define ${\widehat{\lambda }}_{j}^{(0)}={\tilde{\lambda }}_{j}$ for $j=1,...,d.$

For ${\ell }=1,\ldots ,d,$ do

if ${\widehat{\lambda }}_{d-{\ell }+1}^{({\ell }-1)}\gt 4\nu ,$ STOP;

else, put ${\widehat{\lambda }}_{d-{\ell }+1}^{({\ell })}=0$ and

$$\begin{eqnarray*}&&{\widehat{\lambda }}_{j}^{({\ell })}={\widehat{\lambda }}_{j}^{({\ell }-1)}+\displaystyle \frac{1}{d-{\ell }}\left(1-\displaystyle \sum _{k=1}^{d-{\ell }}{\widehat{\lambda }}_{k}^{({\ell }-1)}\right),{\rm{}}\;{\rm{for}}\;{\rm{}}j=1,\ldots ,d-{\ell };\end{eqnarray*}$$

${\ell }={\ell }+1.$

The algorithm checks whether the smallest eigenvalue is larger than the noise level $4\nu$ and if it is not, then sets its value to 0 and distributes the mass of the erased eigenvalue in such a way that they sum to 1. This algorithm is similar to that proposed by Smolin et al  [4], with the important difference that we do not keep all positive eigenvalues but only significantly positive eigenvalues. Here, significant means larger than the noise threshold of the order of the operator-norm error of the LSE. The reason for keeping only eigenvalues above the threshold is similar to the case of the penalized estimator. By proposition 1, eigenvalues below the threshold are of the order of the estimation error $\nu (\epsilon ).$ Therefore, if the state is low rank, such eigenvalues are likely to correspond to zero eigenvalues of the true state and keeping them unchanged would induce significant statistical noise.

If the total number of iterations is $\widehat{{\ell }}=d-\widehat{r}$ then the estimator ${\widehat{\rho }}_{n}^{(\mathrm{phys})}$ has rank $\widehat{r}.$ Its eigenvalues are equal to 0 for $j\gt \widehat{r},$ while, for $j\leqslant \widehat{r}$ they are given by

$$\begin{eqnarray*}&&{\widehat{\lambda }}_{j}^{(\mathrm{phys})}={\tilde{\lambda }}_{j}+\displaystyle \frac{L}{2},\qquad {\rm{}}\;{\rm{where}}\;{\rm{}}\qquad \displaystyle \frac{L}{2}\widehat{r}=\displaystyle \sum _{k\gt \widehat{r}}{\tilde{\lambda }}_{k}=1-\displaystyle \sum _{k\leqslant \widehat{r}}{\tilde{\lambda }}_{k}.\end{eqnarray*}$$

This implies that ${\widehat{\rho }}_{n}^{(\mathrm{phys})}$ has decreasing eigenvalues and ${\widehat{\lambda }}_{\widehat{r}}^{(\mathrm{phys})}\geqslant 4\nu .$ The following theorem shows that ${\widehat{\rho }}_{n}^{(\mathrm{phys})}$ is rank-consistent and its MSE has the same scaling as that of penalized estimator ${\widehat{\rho }}_{n}^{(\mathrm{pen})}.$

Theorem 2. Assume that the state ρ has rank r, i.e. belongs to ${{\mathcal{S}}}_{d,r}.$ For small $\varepsilon \gt 0,$ let $\nu =\nu (\varepsilon )$ be defined as in theorem 1, and assume that ${\lambda }_{r}\gt 8\nu (\varepsilon ).$ Then, with ${{\mathbb{P}}}_{\rho }$ probability larger than $1-\varepsilon$ we have $\widehat{r}=r$ and

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\leqslant 48\cdot r\cdot \nu {(\varepsilon )}^{2}.\end{eqnarray*}$$

Moreover, there exists an absolute constant $C\gt 0$ such that

$$\begin{eqnarray*}&&\underset{\rho \in {{\mathcal{S}}}_{d,r}}{\mathrm{sup}}{E}_{\rho }\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\leqslant C\displaystyle \frac{{rd}}{N}\mathrm{log}\left(\displaystyle \frac{2d}{\varepsilon }\right).\end{eqnarray*}$$

Proof. See appendix A.5.

In order to generate a density matrix of rank r, we first create a rank r upper triangular matrix T in which

• (i)  
  the off-diagonal elements of the first r rows are random complex numbers,
• (ii)  
  the diagonal elements ${T}_{22},...,{T}_{{rr}}$ are real, positive random numbers,
• (iii)  
  all elements of the rows $r+1,...,d$ are zero.

The matrix T is completed by setting T₁₁ such that ${T}_{11}\geqslant 0,$ and $\parallel T{\parallel }_{2}^{2}=1.$ If these conditions cannot be satisfied we repeat the procedure by generating a new set of matrix elements for T. When successful, we set $\rho := {T}^{*}T$ which by construction is a density matrix of rank r. We note that it is not our purpose to generate matrices from a particular 'uniform' ensemble, but merely to have a state with reasonably random eigenvectors, and whose r eigenvalues are not significantly smaller than $1/r.$ Following this procedure we have generated 4 states of 4 ions ($d={2}^{4}$) with ranks $1,2,6,10.$ The rank 6 state for instance, has non-zero eigenvalues $(0.47,0.19,0.12,0.11,0.07,0.04).$

For each state, we have then simulated a number of 100 independent datasets with a given number of repetitions chosen from the range $20,100,500,$ and 2500. In this way we can study the dependence of the MSE of each estimator on state (or rank) and number of repetitions.

We conducted the following simulation study for all the possible combinations between the states and the total number of cycles (i.e. $4\times 4=16$ different scenarios). Below, we denote by r the rank of the 'true' state ρ, from which the data has been generated. The procedure has the following steps:

• (1)  
  For a given number of repetitions n, we simulate 5 independent datasets ${D}_{1},...,{D}_{5},$ each with $n/5$ repetitions. By simply adding the number of counts for each setting and outcome, we obtain a dataset D of n repetitions. However, as we will see below, having 5 separate 'smaller' datasets is important for the purpose of applying cross-validation. Note that such a procedure can be easily implemented in an experimental setting.
• (2)  
  We compute the LSE ${\widehat{\rho }}_{n}^{(\mathrm{ls})}$ based on the full dataset D with total number of cycles n.
• (3)  
  We compute the oracle 'estimator' as follows:

  • (a)  
    We compute the spectral decomposition (4.1) of ${\widehat{\rho }}_{n}^{(\mathrm{ls})},$ with the eigenvalues ${\widehat{\lambda }}_{i}$ arranged in decreasing order of their absolute values. For each rank $1\leqslant \kappa \leqslant d$ we define the truncated (least squares) matrix

    $$\begin{eqnarray*}&&{\widehat{\rho }}_{n}(\kappa )=\displaystyle \sum _{i=1}^{\kappa }{\widehat{\lambda }}_{i}| {\hat{\psi }}_{i}\rangle \langle {\hat{\psi }}_{i}| .\end{eqnarray*}$$
  • (b)  
    We then evaluate the norm-two distance squared $e(\kappa ):= \parallel \rho -{\widehat{\rho }}_{n}(\kappa ){\parallel }_{2}^{2}$ and define the oracle estimator as the truncated estimator with minimal norm two error

    $$\begin{eqnarray*}&&{\widehat{\rho }}_{n}^{(\mathrm{oracle})}={\widehat{\rho }}_{n}({\kappa }_{0}),\qquad {\kappa }_{0}=\mathrm{arg}\;\underset{\kappa }{\mathrm{min}}e(\kappa ).\end{eqnarray*}$$

    Note that the oracle estimator relies on the knowledge of the true state ρ which is not available in a real data set-up. It is nevertheless useful as a benchmark for judging the performance of other estimators in simulation studies. At the next point we define the cross-validation estimator which tries to find the 'optimal' rank ${\kappa }_{0}$ by replacing the unknown state ρ with the LSE computed on a separate batch of data.
• (4)  
  We compute the cross validation estimator as follows.

  • (a)  
    For each $j\in \{1,...,5\}$ we compute the following estimators. While holding the batch D_j out, we compute the LSE ${\widehat{\rho }}_{n;-j}^{(\mathrm{ls})}$ for the dataset consisting of joining the remaining four batches together. Similarly to the point above, we define the rank κ truncation of this estimator by ${\widehat{\rho }}_{n;-j}(\kappa ).$ We also compute the LSE for the remaining batch j, denoted by ${\widehat{\rho }}_{n;j}^{(\mathrm{ls})}.$
  • (b)  
    For each rank κ we evaluate the 'empirical discrepancy'

    $$\begin{eqnarray*}&&{\rm{CV}}(\kappa )=\displaystyle \frac{1}{5}\displaystyle \sum _{i=1}^{5}{\| {\widehat{\rho }}_{n;-j}(\kappa )-{\widehat{\rho }}_{n;j}^{(\mathrm{ls})}\| }_{2}^{2}.\end{eqnarray*}$$

    Since ${\widehat{\rho }}_{n;-j}(\kappa )$ and ${\widehat{\rho }}_{n;j}^{(\mathrm{ls})}$ are independent, and the LSE is unbiased ${\mathbb{E}}({\widehat{\rho }}_{n;j}^{(\mathrm{ls})})=\rho ,$ the expected value of ${\rm{CV}}(k)$ is

    $$\begin{eqnarray*}{\mathbb{E}}\left[{\rm{CV}}(\kappa )\right] & = & {{\mathbb{E}}}_{-1}\left[{\rm{Tr}}\left({\widehat{\rho }}_{n;-1}{(\kappa )}^{2}\right)\right]-2{\rm{Tr}}\left({{\mathbb{E}}}_{-1}\left[{\widehat{\rho }}_{n;-1}(\kappa )\right]\cdot {{\mathbb{E}}}_{1}\left[{\widehat{\rho }}_{n;1}^{(\mathrm{ls})}\right]\right)+{{\mathbb{E}}}_{1}\left[{\rm{Tr}}\left({\left({\widehat{\rho }}_{n;1}^{(\mathrm{ls})}\right)}^{2}\right)\right]\\ & = & {{\mathbb{E}}}_{-1}\left[{\rm{Tr}}\left({\widehat{\rho }}_{n;-1}{(\kappa )}^{2}\right)\right]-2{\rm{Tr}}\left({{\mathbb{E}}}_{-1}\left[{\widehat{\rho }}_{n;-1}(\kappa )\right]\rho \right)+{\rm{Tr}}({\rho }^{2})+C\\ & = & {{\mathbb{E}}}_{-1}\left[\parallel {\widehat{\rho }}_{n;-1}(\kappa )-\rho {\parallel }^{2}\right]+C,\end{eqnarray*}$$

    where we denoted ${{\mathbb{E}}}_{-1}$ and ${{\mathbb{E}}}_{1}$ the expectation over all batches except the first, and respectively over the first batch. Therefore, the average of ${\rm{CV}}(\kappa )$ is equal to the MSE of the truncated estimator ${\widehat{\rho }}_{n;-1}(\kappa ),$ up to a constant C which is independent of κ.
  • (c)  
    Based on the above observation we use the the cross-validation method as a proxy for the oracle estimator. Concretely, we minimize CV $(\kappa )$ with respect to κ

    $$\begin{eqnarray*}&&{\hat{\kappa }}_{\mathrm{cv}}=\underset{k}{\mathrm{arg}\;\mathrm{min}}{\rm{}}\;{\rm{CV}}(\kappa ),\end{eqnarray*}$$

    and define the cross-validation estimator as the truncation to rank ${\hat{\kappa }}_{\mathrm{cv}}$ of the full data LSE ${\widehat{\rho }}_{n}^{(\mathrm{cv})}:= {\widehat{\rho }}_{n}({\hat{\kappa }}_{\mathrm{cv}}).$
• (5)  
  We compute the cross-validated rank-penalized estimator as follows.

  • (a)  
    Let c be a penalization constant chosen from a suitable set of discrete values in the interval $[0,3].$ Similarly to the cross-validation procedure, we hold out batch j, and we compute the rank-penalized estimator (4.2), with penalty constant $c{\nu }^{2}$ for $j=1,\ldots ,5.$ We denote these estimators by ${\widehat{\rho }}_{n;-j}^{(\mathrm{pen})}(c).$ We will also need the LSE ${\rho }_{n;j}^{(\mathrm{ls})}$ for batch j computed above.
  • (b)  
    For each value of c we evaluate the empirical discrepancy

    $$\begin{eqnarray*}&&{\rm{CV}}(c)=\displaystyle \frac{1}{5}\displaystyle \sum _{i=1}^{5}{\| {\widehat{\rho }}_{n;-j}^{(\mathrm{pen})}(c)-{\rho }_{n;j}^{(\mathrm{ls})}\| }_{2}^{2}\end{eqnarray*}$$

    and minimize ${\rm{CV}}(c)$ with respect to the constant c

    $$\begin{eqnarray*}&&\widehat{c}=\underset{c}{\mathrm{arg}\;\mathrm{min}}{\rm{}}\;{\rm{CV}}(c).\end{eqnarray*}$$

    Finally we compute the cross-validated rank penalized estimator ${\widehat{\rho }}_{n}^{(\mathrm{rk}-\mathrm{cv})}$ which is defined as in (4.2), with constant $\widehat{c}{\nu }^{2},$ on the whole dataset D.
• (6)  
  We compute the physical estimator as follows.

  • (a)  
    As above we choose a constant c from a grid over the interval $[0,3].$ We hold out batch j, and we compute the physical threshold estimator (4.5), using the algorithm below this equation, with threshold $c\cdot 4\nu .$ We denote the resulting estimators by ${\widehat{\rho }}_{n;-j}^{(\mathrm{phys})}(c),$ for $j=1,\ldots ,5.$
  • (b)  
    For each value of c we evaluate the empirical discrepancy

    $$\begin{eqnarray*}&&{\rm{CV}}(c)=\displaystyle \frac{1}{5}\displaystyle \sum _{i=1}^{5}{\| {\widehat{\rho }}_{n;-j}^{(\mathrm{phys})}(c)-{\widehat{\rho }}_{n;j}^{(\mathrm{ls})}\| }_{2}^{2}.\end{eqnarray*}$$

    We then minimize ${\rm{CV}}(c)$ with respect to the constant c.

    $$\begin{eqnarray*}&&\hat{c}=\underset{c}{\mathrm{arg}\;\mathrm{min}}{\rm{}}\;{\rm{CV}}(c)\end{eqnarray*}$$

    Finally we compute the cross-validated physical estimator ${\widehat{\rho }}_{n}^{(\mathrm{phys}-\mathrm{cv})}$ which is defined as in (4.5), with constant $\hat{c}\cdot 4\nu .$

We collect here the results of the simulation study described in the previous section. As a figure of merit we focus on the MSE ${\mathbb{E}}\left((\parallel {\widehat{\rho }}_{n}-\rho {\parallel }_{2}^{2}\right)$ of each estimator, which is estimated by averaging the square errors over the 100 independent repetitions of the procedure. We are also interested in how the different methods perform relative to each other, and whether the selected rank is consistent, i.e. it concentrates on the rank of the true state for large number of repetitions.

The four panels in figure 3 represent the boxplots of the square errors $\parallel {\widehat{\rho }}_{n}-\rho {\parallel }_{2}^{2}$ for the different estimators, and different states, when the number of repetitions is n = 20. Similarly, figure 4 shows the same boxplots at n = 100. As expected, in both cases the least squares performs significantly worse than the other estimators, and the discrepancy is larger for small rank states. The remaining 4 estimators have similar MSE's with the physical one performing slightly better than the rest, followed by the oracle. Note also that the estimators' variances (indicated by the size of the boxes) are larger for the least squares than the other estimators. A similar behaviour has been observed for $n=500,2500$ repetitions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 illustrates the dependence of the MSE of a given estimator, as a function of n, for the four different states which have been analysed. Since the MSE decreases as ${n}^{-1}$ we have chosen to plot the 'renormalized' MSE given by $n\cdot {\mathbb{E}}\parallel {\widehat{\rho }}_{n}-\rho {\parallel }_{2}^{2},$ which converges to a constant value for large n. As expected the limiting value increases with the rank of the state, as a proxy for the number of parameters to be estimated.

Zoom In Zoom Out Reset image size

The histograms in figure 6 show the probability distributions for the chosen rank of each given estimator, as a function of the number of measurement repetitions n, for the state of rank 6. We note that in all cases the proportion of times that the chosen rank is equal to the true rank of the state increases as with the number of repetitions. However, this convergence towards a 'rank-consistent' estimator is rather slow, as the proportion surpasses 80% only when n = 2500. Another observation is that the penalty and threshold estimators appear to have different behaviours: the former tends to underestimate the true rank, while the latter tends to overestimate it. As expected, the oracle estimator is more likely to choose the correct rank for large number of repetitions. Perhaps slightly more surprising, for small number of repetitions (n = 20), the oracle choose a pure state in most cases.

Zoom In Zoom Out Reset image size

Proof. Since the matrix elements of ${\bf{A}}$ are ${{\bf{A}}}_{{\bf{b}},({\bf{o}},{\bf{s}})}={A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}),$ the product ${{\bf{A}}}^{*}\cdot {\bf{A}}$ has elements

$$\begin{eqnarray*}&&{[{{\bf{A}}}^{*}\cdot {\bf{A}}]}_{{\bf{b}},{{\bf{b}}}^{\prime }}=\displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{{\bf{o}}}{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}){A}_{{{\bf{b}}}^{\prime }}({\bf{o}}| {\bf{s}}).\end{eqnarray*}$$

If ${\bf{b}}={{\bf{b}}}^{\prime }$ it is easy to see that

$$\begin{eqnarray*}&&{[{{\bf{A}}}^{*}\cdot {\bf{A}}]}_{{\bf{b}},{\bf{b}}}=\displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{{\bf{o}}}\displaystyle \prod _{j\not\hspace{-2pt}{\in }{E}_{{\bf{b}}}}{\delta }_{{s}_{j},{b}_{j}}={2}^{k}{3}^{d({\bf{b}})}.\end{eqnarray*}$$

If ${\bf{b}}\ne {{\bf{b}}}^{\prime },$ we have either ${E}_{{\bf{b}}}={E}_{{{\bf{b}}}^{\prime }}$ or ${E}_{{\bf{b}}}\;\not\hspace{-2pt}{=}\;{E}_{{{\bf{b}}}^{\prime }}.$ On the one hand, in case the sets ${E}_{{\bf{b}}}$ and ${E}_{{{\bf{b}}}^{\prime }}$ are equal, we have

$$\begin{eqnarray*}&&{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}){A}_{{{\bf{b}}}^{\prime }}({\bf{o}}| {\bf{s}})=\displaystyle \prod _{j\not\hspace{-2pt}{\in }{E}_{{\bf{b}}}}{\delta }_{{s}_{j},{b}_{j}}\cdot {\delta }_{{s}_{j},{b}_{j}^{\prime }}.\end{eqnarray*}$$

For each ${\bf{s}},$ the previous product is 0. Indeed, if different from 0 then ${b}_{j}={b}_{j}^{\prime }$ for all j not in ${E}_{{\bf{b}}}.$ As ${b}_{j}={b}_{{j}^{\prime }}=I$ for j in ${E}_{{\bf{b}}},$ it implies that ${\bf{b}}={{\bf{b}}}^{\prime }$ which contradicts the assumption here.

On the other hand, if the sets ${E}_{{\bf{b}}}$ and ${E}_{{{\bf{b}}}^{\prime }}$ are different, there exists at least one coordinate j₀ in the symmetric difference ${E}_{{\bf{b}}}{\rm{\Delta }}{E}_{{b}^{\prime }}$ and the sum over outcomes ${\bf{o}}$ will split over values of ${\bf{o}}$ where ${o}_{{j}_{0}}=1$ and values where ${o}_{{j}_{0}}=-1:$

$$\begin{eqnarray*} & & \displaystyle \sum _{{\bf{o}}}{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}){A}_{{{\bf{b}}}^{\prime }}({\bf{o}}| {\bf{s}})\\ & & =\displaystyle \sum _{{\bf{o}}:{o}_{{j}_{0}}=1}{\delta }_{{s}_{{j}_{0}},{b}_{{j}_{0}}}\displaystyle \prod _{j\not\hspace{-2pt}{\in }{E}_{{\bf{b}}}}{o}_{j}{\delta }_{{s}_{j},{b}_{j}}\cdot \displaystyle \prod _{l\not\hspace{-2pt}{\in }{E}_{{{\bf{b}}}^{\prime }}}{o}_{l}{\delta }_{{s}_{l},{b}_{l}^{\prime }}\\ & & -\displaystyle \sum _{{\bf{o}}:{o}_{{j}_{0}}=-1}{\delta }_{{s}_{{j}_{0}},{b}_{{j}_{0}}}\displaystyle \prod _{j\not\hspace{-2pt}{\in }{E}_{{\bf{b}}}}{o}_{j}{\delta }_{{s}_{j}={b}_{j}}\cdot \displaystyle \prod _{l\not\hspace{-2pt}{\in }{E}_{{{\bf{b}}}^{\prime }}}{o}_{l}{\delta }_{{s}_{l},{b}_{l}^{\prime }}=0.\end{eqnarray*}$$

We assumed here that j₀ belongs to ${E}_{{\bf{b}}}\setminus {E}_{{{\bf{b}}}^{\prime }}$ and the same holds for j₀ in ${E}_{{{\bf{b}}}^{\prime }}\setminus {E}_{{\bf{b}}}.$

Proof of proposition 1. For more clarity, in this proof we will use the notation $I(a=b)={\delta }_{a,b}.$ Note that $f({\bf{o}}| {\bf{s}})=\frac{1}{n}{\displaystyle \sum }_{i}I({X}_{{\bf{s}},i}={\bf{o}}),$ where the random variables ${X}_{{\bf{s}},i}$ are independent for all settings ${\bf{s}}$ and all i from 1 to n. To estimate the risk of the linear estimator we write

$$\begin{eqnarray*}{\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho & = & \displaystyle \sum _{{\bf{b}}}\displaystyle \sum _{{\bf{o}}}\displaystyle \sum _{{\bf{s}}}(f({\bf{o}}| {\bf{s}}))-p({\bf{o}}| {\bf{s}}))\displaystyle \frac{{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}})}{{2}^{k}{3}^{d({\bf{b}})}}{\sigma }_{{\bf{b}}}\\ & = & \displaystyle \sum _{{\bf{b}}}\displaystyle \sum _{{\bf{o}}}\displaystyle \sum _{{\bf{s}}}\displaystyle \frac{1}{n}\displaystyle \sum _{i}(I({X}_{{\bf{s}},i}={\bf{o}})-p({\bf{o}}| {\bf{s}}))\displaystyle \frac{{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}})}{{2}^{k}{3}^{d({\bf{b}})}}{\sigma }_{{\bf{b}}}\\ & := & \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{i}{W}_{{\bf{s}},i},\end{eqnarray*}$$

where ${W}_{{\bf{s}},i}$ are independent and centred Hermitian random matrices. We will apply the following extension of the Bernstein matrix inequality [39] due to [40], see also [32, 33].

Proposition 3. (Bernstein inequality, [40].) Let ${Y}_{1},\ldots ,{Y}_{n}$ be independent, centred, m × m Hermitian random matrices. Suppose that, for some constants $V,W\gt 0$ we have $\parallel {Y}_{j}\parallel \leqslant V,$ for all j from 1 to n, and that $\parallel {\displaystyle \sum }_{j}{\mathbb{E}}({Y}_{j}^{2})\parallel \leqslant W.$ Then, for all $t\geqslant 0,$

$$\begin{eqnarray*}&&{\mathbb{P}}(\parallel {Y}_{1}+...+{Y}_{n}\parallel \geqslant t)\leqslant 2m\mathrm{exp}\left(-\displaystyle \frac{{t}^{2}/2}{W+{tV}/3}\right).\end{eqnarray*}$$

In our setup we bound $\parallel {W}_{{\bf{s}},i}\parallel \leqslant V$ for all ${\bf{s}}$ and i and $\parallel {\displaystyle \sum }_{{\bf{s}}}{\mathbb{E}}({W}_{{\bf{s}},i}^{*}{W}_{{\bf{s}},i})\parallel \leqslant W,$ where $V,W$ are evaluated below. We have

$$\begin{eqnarray*}\parallel {W}_{{\bf{s}},i}\parallel & \leqslant & \displaystyle \frac{1}{n}\displaystyle \sum _{{\bf{b}}}\displaystyle \sum _{{\bf{o}}}\left|\displaystyle \frac{{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}})}{{2}^{k}{3}^{d({\bf{b}})}}\right|\cdot | I({X}_{{\bf{s}},i}={\bf{o}})-p({\bf{o}}| {\bf{s}})| \cdot \parallel {\sigma }_{{\bf{b}}}\parallel \\ & \leqslant & \displaystyle \frac{1}{n}\displaystyle \sum _{{\bf{b}}}\displaystyle \frac{1}{{2}^{k}{3}^{d(b)}}\displaystyle \prod _{j\not\hspace{-2pt}{\in }{E}_{{\bf{b}}}}I({b}_{j}={s}_{j})\displaystyle \sum _{{\bf{o}}}| I({X}_{{\bf{s}},i}={\bf{o}})-p({\bf{o}}| {\bf{s}})| \\ & \leqslant & \displaystyle \frac{2}{n{2}^{k}}\displaystyle \sum _{{\ell }=0}^{k}\;\displaystyle \sum _{b:d(b)={\ell }}\displaystyle \frac{1}{{3}^{{\ell }}}=\displaystyle \frac{2}{n{2}^{k}}\displaystyle \sum _{{\ell }=0}^{k}\displaystyle \left(\genfrac{}{}{0em}{}{k}{{\ell }}\right)\displaystyle \frac{1}{{3}^{{\ell }}}=\displaystyle \frac{2}{n{2}^{k}}{\left(1+\displaystyle \frac{1}{3}\right)}^{k}=\displaystyle \frac{2}{n}{\left(\displaystyle \frac{2}{3}\right)}^{k}:= V.\end{eqnarray*}$$

Let us denote $B({\bf{o}}| {\bf{s}}):= {\displaystyle \sum }_{{\bf{b}}}{2}^{-k}{3}^{-d({\bf{b}})}{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}){\sigma }_{{\bf{b}}}.$ Then

$$\begin{eqnarray} & & \parallel \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{i}{\mathbb{E}}({W}_{{\bf{s}},i}^{*}{W}_{{\bf{s}},i})\parallel \\ & = & \displaystyle \frac{1}{{n}^{2}}\parallel \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{i}\displaystyle \sum _{{\bf{o}},{{\bf{o}}}^{\prime }}{B}^{*}({\bf{o}}| {\bf{s}})\cdot {\rm{Cov}}(I({X}_{{\bf{s}},i}={\bf{o}}),I({X}_{{\bf{s}},i}={{\bf{o}}}^{\prime }))\cdot B({{\bf{o}}}^{\prime }| {\bf{s}})\parallel \\ & \leqslant & \displaystyle \frac{1}{{n}^{2}}\parallel \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{i}\displaystyle \sum _{{\bf{o}}}{B}^{*}({\bf{o}}| {\bf{s}})\cdot B({\bf{o}}| {\bf{s}})\parallel .\end{eqnarray} \tag{ A.1 }$$

The last inequality follows from the fact that the covariance matrix can be written as the difference of two positive matrices ${{\rm{Cov}}}_{{\bf{o}},{{\bf{o}}}^{\prime }}:= {\rm{Cov}}(I({X}_{{\bf{s}},i}={\bf{o}}),I({X}_{{\bf{s}},i}={{\bf{o}}}^{\prime }))=p({\bf{o}}| {\bf{s}}){\delta }_{{\bf{o}},{{\bf{o}}}^{\prime }}-p({\bf{o}}| {\bf{s}})p({{\bf{o}}}^{\prime }| {\bf{s}})$ and since $p({\bf{o}}| {\bf{s}})\leqslant 1,$ we get ${\rm{Cov}}\leqslant {\bf{1}}.$

We replace $B({\bf{o}}| {\bf{s}})$ in (A.1) and use lemma 1 to get

$$\begin{eqnarray*}&&\parallel \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{i}{\mathbb{E}}({W}_{{\bf{s}},i}^{*}{W}_{{\bf{s}},i})\parallel \\ &&\quad \leqslant \displaystyle \frac{1}{n}\| \displaystyle \sum _{{\bf{s}}}\displaystyle \sum _{{\bf{o}}}\displaystyle \sum _{{\bf{b}},{{\bf{b}}}^{\prime }}\displaystyle \frac{1}{{2}^{2k}{3}^{d({\bf{b}})+d({{\bf{b}}}^{\prime })}}{A}_{{\bf{b}}}({\bf{o}}| {\bf{s}}){A}_{{{\bf{b}}}^{\prime }}({\bf{o}}| {\bf{s}})\cdot {\sigma }_{{\bf{b}}}{\sigma }_{{{\bf{b}}}^{\prime }}\| \\ &&\quad \leqslant \displaystyle \frac{1}{n{2}^{k}}\displaystyle \sum _{{\ell }=0}^{k}\;\displaystyle \sum _{{\bf{b}}:d({\bf{b}})={\ell }}\displaystyle \frac{1}{{3}^{{\ell }}}=\displaystyle \frac{1}{n{2}^{k}}\displaystyle \sum _{{\ell }=0}^{k}\displaystyle \left(\genfrac{}{}{0em}{}{k}{{\ell }}\right)\displaystyle \frac{1}{{3}^{{\ell }}}=\displaystyle \frac{1}{n}{\left(\displaystyle \frac{2}{3}\right)}^{k}:= W.\end{eqnarray*}$$

Apply the matrix Bernstein inequality in proposition 3 to get, for any $t\gt 0:$

$$\begin{eqnarray*}&&{{\mathbb{P}}}_{\rho }\left(\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel \geqslant t\right)\leqslant {2}^{k+1}\mathrm{exp}\left(-\displaystyle \frac{{t}^{2}}{1+2t/3}\displaystyle \frac{n}{2}{\left(\displaystyle \frac{3}{2}\right)}^{k}\right).\end{eqnarray*}$$

We choose $t\gt 0$ such that

$$\begin{eqnarray*}&&\varepsilon ={2}^{k+1}\mathrm{exp}\left(-\displaystyle \frac{{t}^{2}}{1+2t/3}\cdot \displaystyle \frac{n}{2}{\left(\displaystyle \frac{3}{2}\right)}^{k}\right),\end{eqnarray*}$$

which leads to t such that

$$\begin{eqnarray*}&&\displaystyle \frac{{t}^{2}}{1+2t/3}=\displaystyle \frac{2}{n}{\left(\displaystyle \frac{2}{3}\right)}^{k}\mathrm{log}\left(\displaystyle \frac{{2}^{k+1}}{\varepsilon }\right)=\displaystyle \frac{3}{2}v(\varepsilon ).\end{eqnarray*}$$

Then the convenient choice of t is $\nu (\varepsilon )=v(\varepsilon )/2+\sqrt{{v}^{2}(\varepsilon )/2+3/2\cdot v(\varepsilon )},$ which is equivalent to $\sqrt{3/2\cdot v(\varepsilon )}$ when this last term tends to 0.□

Proof of proposition 1. Let us denote by $({\widehat{\lambda }}_{1}\geqslant ...\geqslant {\widehat{\lambda }}_{d})$ the eigenvalues of $| {\widehat{\rho }}_{n}^{(\mathrm{ls})}|$ and by ${\tilde{\lambda }}_{1},\ldots ,{\tilde{\lambda }}_{d}$ the eigenvalues of the resulting estimator ${\widehat{\rho }}_{n}^{(\mathrm{ls},n)}.$ It is easy to see that the latter has the same eigenvectors as ${\widehat{\rho }}_{n}^{(\mathrm{ls})}.$

Therefore

$$\begin{eqnarray*}&&({\tilde{\lambda }}_{1},\ldots ,{\tilde{\lambda }}_{d})=\underset{\lambda :{\displaystyle \sum }_{j=1}^{d}{\lambda }_{j}=1}{\mathrm{argmin}}\;{\displaystyle \sum }_{j=1}^{d}{({\lambda }_{j}-{\widehat{\lambda }}_{j})}^{2}.\end{eqnarray*}$$

This optimization has the explicit solution

$$\begin{eqnarray*}&&{\tilde{\lambda }}_{j}={\widehat{\lambda }}_{j}+\displaystyle \frac{L}{2},{\rm{where}}\;{\rm{}}\displaystyle \frac{L}{2}d=1-\displaystyle \sum _{j=1}^{d}{\widehat{\lambda }}_{j}.\end{eqnarray*}$$

Note that

$$\begin{eqnarray*}&&\displaystyle \frac{L}{2}=\displaystyle \frac{1}{\;}d{\rm{tr}}(\rho -{\widehat{\rho }}_{n}^{(\mathrm{ls})})\leqslant \parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel .\end{eqnarray*}$$

Therefore, $\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls},n)}-{\widehat{\rho }}_{n}^{(\mathrm{ls})}\parallel =L/2\leqslant \parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel$ and thus

$$\begin{eqnarray*}\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls},n)}-\rho \parallel & \leqslant & \parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel +\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls},n)}-{\widehat{\rho }}_{n}^{(\mathrm{ls})}\parallel \leqslant 2\nu (\varepsilon ).\end{eqnarray*}$$

□

Proof of corollary 1. We order the eigenvalues of ${\widehat{\rho }}_{n}^{(\mathrm{pen})}$ in decreasing order of absolute values $(| {\widehat{\lambda }}_{1}| \geqslant ...\geqslant | {\widehat{\lambda }}_{d}| ),$ and denote by ${\tilde{\lambda }}_{1},\ldots ,{\tilde{\lambda }}_{d}$ the eigenvalues of ${\widehat{\rho }}_{n}^{(\mathrm{pen},n)}.$ As in the previous proof, we can see that both matrices will share the same eigenvectors and that the relation between the eigenvalues is

$$\begin{eqnarray*}&&{\tilde{\lambda }}_{j}={\widehat{\lambda }}_{j}+\displaystyle \frac{L}{2},{\rm{}}\;{\rm{where}}\;{\rm{}}\displaystyle \frac{L}{2}d=1-\displaystyle \sum _{j=1}^{d}{\widehat{\lambda }}_{j}.\end{eqnarray*}$$

Thus

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-{\widehat{\rho }}_{n}^{(\mathrm{pen},n)}{\parallel }_{F}^{2}=d{(L/2)}^{2}\leqslant {{\rm{Tr}}}^{2}({\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho )/d\leqslant \parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{F}^{2}.\end{eqnarray*}$$

We deduce that $\parallel \rho -{\widehat{\rho }}_{n}^{(\mathrm{pen},n)}{\parallel }_{F}\leqslant \parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-{\widehat{\rho }}_{n}^{(\mathrm{pen},n)}{\parallel }_{F}+\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{F}\leqslant 2\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen})}-\rho {\parallel }_{F}.$

Therefore

$$\begin{eqnarray*}&&{\mathbb{P}}\left(\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{F}\gt 8c(\theta )r\nu {(\varepsilon )}^{2}\right)\lt \varepsilon \end{eqnarray*}$$

and the previous inequality remains true for all $0\lt e\leqslant \varepsilon \lt 1.$ Let us denote by

$$\begin{eqnarray*}&&x(e)=8r{\nu }^{2}(e)=16\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{e}),\end{eqnarray*}$$

which is a decreasing function of e. This implies that $e=2d\mathrm{exp}\left(-\displaystyle \frac{{Nx}(e)}{16{rd}}\right).$ Thus

$$\begin{eqnarray*}{{\mathbb{E}}}_{\rho }\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2} & = & {\displaystyle \int }_{0}^{\infty }{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x\\ & = & {\displaystyle \int }_{0}^{x(\varepsilon )}{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x+{\displaystyle \int }_{x(\varepsilon )}^{\infty }{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{pen},n)}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x\\ & \leqslant & x(\varepsilon )+{\displaystyle \int }_{x(\varepsilon )}^{\infty }2d\mathrm{exp}(-\displaystyle \frac{{Nx}}{16{rd}}){\rm{d}}x\\ & \leqslant & 16\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{\varepsilon })+2d\cdot 16\displaystyle \frac{{rd}}{N}\mathrm{exp}(-\displaystyle \frac{N}{16{rd}}x(\varepsilon ))\\ & \leqslant & 16\displaystyle \frac{{rd}}{N}\left(\mathrm{log}(\displaystyle \frac{2d}{\varepsilon })+\varepsilon \right)\leqslant C\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{\varepsilon }).\end{eqnarray*}$$

□

Proof of theorem 2. We recall from proposition 2 that with probability larger than $1-\varepsilon ,$ we have $\parallel {\tilde{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel \leqslant 2\nu (\varepsilon ).$ In particular, by using the Weyl inequality [41], this implies that $| {\tilde{\lambda }}_{k}-{\lambda }_{k}| \leqslant 2\nu (\varepsilon )$ for all k from 1 to d, where ${\lambda }_{1},\ldots ,{\lambda }_{d}$ are the eigenvalues of ρ arranged in decreasing order. After a total of $\hat{{\ell }}=d-\hat{r}$ iterations the algorithm stops and we have [42]

$$\begin{eqnarray*}&&{\widehat{\lambda }}_{\hat{r}+1}^{(\mathrm{phys})}=...={\widehat{\lambda }}_{d}^{(\mathrm{phys})}=0\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\widehat{\lambda }}_{j}^{(\mathrm{phys})}:= {\widehat{\lambda }}_{j}^{(\hat{{\ell }})}={\tilde{\lambda }}_{j}+\displaystyle \frac{1}{\hat{r}}\displaystyle \sum _{k\gt \hat{r}}{\tilde{\lambda }}_{k}\qquad j=1,...,\hat{r}.\end{eqnarray*}$$

From now on we assume that $\parallel {\tilde{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel \leqslant 2\nu (\varepsilon )$ which holds with probability larger than $1-\varepsilon ,$ see corollary 3.2. We will show that under this assumption $\hat{r}=r.$ For this we consider the two cases $\hat{r}\gt r$ and $\hat{r}\lt r$ separately.

Suppose $\hat{r}\gt r.$ For $j\leqslant \hat{r}$ we have

$$\begin{eqnarray*}| {\widehat{\lambda }}_{j}^{(\mathrm{phys})}-{\lambda }_{j}| & = & \left|{\tilde{\lambda }}_{j}-{\lambda }_{j}+\displaystyle \frac{1}{\hat{r}}\left(1-\displaystyle \sum _{k\leqslant \hat{r}}{\tilde{\lambda }}_{k}\right)\right|\\ & = & \left|{\tilde{\lambda }}_{j}-{\lambda }_{j}+\displaystyle \frac{1}{\hat{r}}\displaystyle \sum _{k\leqslant \hat{r}}({\lambda }_{k}-{\tilde{\lambda }}_{k})\right|\\ & \leqslant & | {\tilde{\lambda }}_{j}-{\lambda }_{j}| +\displaystyle \frac{1}{\hat{r}}\displaystyle \sum _{k\leqslant \hat{r}}| {\tilde{\lambda }}_{k}-{\lambda }_{k}| \leqslant 2\nu +2\nu =4\nu .\end{eqnarray*}$$

This implies that ${\widehat{\lambda }}_{\hat{r}}^{(\mathrm{phys})}\leqslant {\lambda }_{\hat{r}}+4\nu =4\nu$ since $\hat{r}\gt r.$ However, as discussed above, by the definition of the threshold estimator, ${\widehat{\lambda }}_{\hat{r}}^{(\mathrm{phys})}\geqslant 4\nu .$ Therefore $\hat{r}$ cannot be larger than r.

Suppose $\hat{r}\lt r.$ Then ${\widehat{\lambda }}_{r}^{(\mathrm{phys})}\geqslant {\lambda }_{r}-| {\widehat{\lambda }}_{r}^{(\mathrm{phys})}-{\lambda }_{r}| \gt 8\nu -4\nu =4\nu .$ However this is not possible since ${\widehat{\lambda }}_{r}^{(\mathrm{phys})}=0$ under the assumption $\hat{r}\lt r.$ Thus, $\widehat{r}=r$ with probability larger than $1-\epsilon .$

We consider now the MSE of the estimator. As shown above, we have $| {\widehat{\lambda }}_{j}^{(\mathrm{phys})}-{\lambda }_{j}| \leqslant 4\nu .$ Let $\rho ={{UDU}}^{*}$ be the diagonalization of ρ, where U is unitary and D is the diagonal matrix of eigenvalues. Similarly, let ${\widehat{\rho }}_{n}^{(\mathrm{ls})}={\widehat{U}}_{n}{D}_{n}^{(\mathrm{ls})}{\widehat{U}}_{n}^{*}$ and ${\widehat{\rho }}_{n}^{(\mathrm{phys})}={\widehat{U}}_{n}{D}_{n}^{(\mathrm{phys})}{\widehat{U}}_{n}^{*}$ be the decompositions of the least squares and the physical threshold estimator, where we take into account that the latter two have the same eigenbasis. Then

$$\begin{eqnarray}\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2} & = & \parallel {\widehat{U}}_{n}{D}_{n}^{(\mathrm{phys})}{\widehat{U}}_{n}^{*}-{{UDU}}^{*}{\parallel }_{2}\\ & \leqslant & \parallel {\widehat{U}}_{n}{D}_{n}^{(\mathrm{phys})}{\widehat{U}}_{n}^{*}-{\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}{\parallel }_{2}+\parallel {\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}-{{UDU}}^{*}{\parallel }_{2}\\ & \leqslant & \parallel {D}_{n}^{(\mathrm{phys})}-D{\parallel }_{2}+\parallel {\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}-{{UDU}}^{*}{\parallel }_{2}\end{eqnarray} \tag{ A.2 }$$

The first norm on the right-hand side of the previous inequality is bounded as

$$\begin{eqnarray}&&\parallel {D}_{n}^{(\mathrm{phys})}-D{\parallel }_{2}^{2}=\displaystyle \sum _{j\leqslant \widehat{r}}{({\widehat{\lambda }}_{j}^{(\mathrm{phys})}-{\lambda }_{j})}^{2}\leqslant 16r{\nu }^{2}.\end{eqnarray} \tag{ A.3 }$$

For the second norm we use a similar triangle inequality for the operator norm

$$\begin{eqnarray*}\parallel {\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}-{{UDU}}^{*}\parallel & \leqslant & \parallel {\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}-{\widehat{U}}_{n}{D}_{n}^{(\mathrm{ls})}{\widehat{U}}_{n}^{*}\parallel +\parallel {\widehat{U}}_{n}{D}_{n}^{(\mathrm{ls})}{\widehat{U}}_{n}^{*}-{{UDU}}^{*}\parallel \\ & = & \parallel D-{D}_{n}^{(\mathrm{ls})}\parallel +\parallel {\widehat{\rho }}_{n}^{(\mathrm{ls})}-\rho \parallel \leqslant \nu +\nu =2\nu .\end{eqnarray*}$$

The first term is smaller than ν since $| {\lambda }_{i}-{\widehat{\lambda }}_{i}^{(\mathrm{ls})}| \leqslant \nu$ as it follows from proposition 1, and the Weyl inequality [41]. The second term is also bounded by ν, by proposition 1. Therefore, since ${\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}$ and ${{UDU}}^{*}$ are rank r matrices, the difference is at most rank $2r$ and $\parallel {\widehat{U}}_{n}D{\widehat{U}}_{n}^{*}-{{UDU}}^{*}{\parallel }_{2}\leqslant 2\nu \sqrt{2r}.$ By plugging this together with (A.3) into the right side of (A.2) we get

$$\begin{eqnarray*}&&\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\leqslant {(2\sqrt{2}+4)}^{2}r{\nu }^{2},\end{eqnarray*}$$

with probability larger than $1-\epsilon .$

Moreover, the previous inequality remains true for all $0\lt e\leqslant \varepsilon \lt 1.$ Let us denote by

$$\begin{eqnarray*}&&x(e)=48r{\nu }^{2}(e)=96\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{e}),\end{eqnarray*}$$

which is a decreasing function of e. This implies that $e=2d\mathrm{exp}\left(-\displaystyle \frac{{Nx}(e)}{96{rd}}\right).$ Thus

$$\begin{eqnarray*}{{\mathbb{E}}}_{\rho }\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2} & = & {\displaystyle \int }_{0}^{\infty }{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x\\ & = & {\displaystyle \int }_{0}^{x(\varepsilon )}{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x+{\displaystyle \int }_{x(\varepsilon )}^{\infty }{{\mathbb{P}}}_{\rho }(\parallel {\widehat{\rho }}_{n}^{(\mathrm{phys})}-\rho {\parallel }_{2}^{2}\gt x){\rm{d}}x\\ & \leqslant & x(\varepsilon )+{\displaystyle \int }_{x(\varepsilon )}^{\infty }2d\mathrm{exp}(-\displaystyle \frac{{Nx}}{96{rd}}){\rm{d}}x\\ & \leqslant & 96\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{\varepsilon })+2d\cdot 96\displaystyle \frac{{rd}}{N}\mathrm{exp}(-\displaystyle \frac{N}{96{rd}}x(\varepsilon ))\\ & \leqslant & 96\displaystyle \frac{{rd}}{N}\left(\mathrm{log}(\displaystyle \frac{2d}{\varepsilon })+\varepsilon )\right)\leqslant C\displaystyle \frac{{rd}}{N}\mathrm{log}(\displaystyle \frac{2d}{\varepsilon }).\end{eqnarray*}$$

In this section we present a detailed calculation of the average quantum Fisher information at a the rank r state ${\rho }_{0}$ defined in (5.2), for random measurements with uniform distribution over the measurement basis. We consider the full parametrization by $\theta \in {{\mathbb{R}}}^{{d}^{2}}$ given by equation (5.3). As explained in the proof, the Fisher information matrix for the parametrization $\tilde{\theta }$ of the rotation model ${{\mathcal{R}}}_{d,r}$ is a particular block of the larger Fisher matrix computed here. We will come back to this at the end of the computation.

Let

$$\begin{eqnarray*}&&{{\bf{B}}}_{U}:= \{| {\bf{o}};U\rangle := U| {\bf{o}}\rangle :{\bf{o}}=1,...,{2}^{k}\}\end{eqnarray*}$$

denote the ONB obtained by rotating the standard basis by the unitary U.

The Fisher info $I(\rho | {{\bf{B}}}_{U})$ associated to the measurement with ONB ${{\bf{B}}}_{U}$ is given by the matrix

$$\begin{eqnarray*}&&I{(\rho | {{\bf{B}}}_{U})}_{a,b}=\displaystyle \sum _{{\bf{o}}:{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})\gt 0}\displaystyle \frac{1}{{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}\displaystyle \frac{\partial {p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\partial {\theta }_{a}}\displaystyle \frac{\partial {p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\partial {\theta }_{b}},\end{eqnarray*}$$

where ${p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})$ is the probability of the outcome ${\bf{o}}$

$$\begin{eqnarray*}{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U}) & = & \langle {\bf{o}},U| \rho | {\bf{o}},U\rangle =\displaystyle \sum _{i}{\rho }_{{ii}}| \langle {\bf{o}},U| i\rangle {| }^{2}\\ & + & 2\displaystyle \sum _{i\lt j}{\rm{Re}}({\rho }_{i,j}){\rm{Re}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle )+2\displaystyle \sum _{i\lt j}{\rm{Im}}({\rho }_{i,j}){\rm{Im}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle )\end{eqnarray*}$$

and the partial derivatives are given by

$$\begin{eqnarray*}\displaystyle \frac{\partial {p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\partial {\rho }_{i,i}} & = & | \langle {\bf{o}},U| i\rangle {| }^{2}\\ \displaystyle \frac{\partial {p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\partial {\rm{Re}}{\rho }_{i,j}} & = & 2{\rm{Re}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle )\\ \displaystyle \frac{\partial {p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\partial {\rm{Im}}{\rho }_{i,j}} & = & 2{\rm{Im}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle ).\end{eqnarray*}$$

The Fisher information matrix $I(\rho | {{\bf{B}}}_{U})$ has the following block structure

$$\begin{eqnarray*}I(\rho | {{\bf{B}}}_{U})=\left(\begin{array}{ccc}{I}^{{dd}}(\rho | {{\bf{B}}}_{U}) & {I}^{{dr}}(\rho | {{\bf{B}}}_{U}) & {I}^{{di}}(\rho | {{\bf{B}}}_{U})\\ {I}^{{rd}}(\rho | {{\bf{B}}}_{U}) & {I}^{{rr}}(\rho | {{\bf{B}}}_{U}) & {I}^{{ri}}(\rho | {{\bf{B}}}_{U})\\ {I}^{{id}}(\rho | {{\bf{B}}}_{U}) & {I}^{{ir}}(\rho | {{\bf{B}}}_{U}) & {I}^{{ii}}(\rho | {{\bf{B}}}_{U})\end{array}\right)\end{eqnarray*}$$

with superscripts indicating the type of parameter considered: diagonal, real or imaginary part of off-diagonal element. The average Fisher information for a randomly chosen basis is

$$\begin{eqnarray*}&&\bar{I}(\rho ):= \int \mu ({\rm{d}}U)I(\rho | {{\bf{B}}}_{U}),\end{eqnarray*}$$

where $\mu ({\rm{d}}U)$ is the Haar measure over unitaries used for choosing the random basis. Note that by symmetry $\bar{I}(\rho )$ only depends on the spectrum of ρ. We will not compute $\bar{I}(\rho )$ for an arbitrary state ρ but only at $\rho ={\rho }_{0}.$ The corresponding Fisher information will be denoted $\bar{I}=\bar{I}({\rho }_{0})$ and is a function of d and r. Below we compute the different blocks of $\bar{I}.$ For the matrix elements we will use a suggestive notation, e.g. ${\bar{I}}_{{ii};{jk}}^{d,r}$ denotes the element corresponding to the diagonal parameter ${\theta }_{{ii}}^{d}={\rho }_{{ii}}$ and the real part of the off-diagonal element ${\rho }_{{jk}},$ etc.

(A) Diagonal–diagonal block.

$$\begin{eqnarray*}&&I{(\rho | {{\bf{B}}}_{U})}_{{ii};{jj}}^{{dd}}=\displaystyle \sum _{{\bf{o}}:{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})\gt 0}\displaystyle \frac{1}{{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}| \langle {\bf{o}},U| i\rangle {| }^{2}| \langle {\bf{o}},U| j\rangle {| }^{2}.\end{eqnarray*}$$

By integrating over unitaries we obtain the corresponding matrix element of the average Fisher matrix $\bar{I}.$ Since ${p}_{{\rho }_{0}}({\bf{o}}| {{\bf{B}}}_{U})\gt 0$ is true for all ${\bf{o}},$ with probability one with respect to ${\mu }^{d}({\rm{d}}U),$ we drop the condition from the sum. At the state $\rho ={\rho }_{0}$ defined in (5.2), we have

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\displaystyle \sum _{{\bf{o}}}\displaystyle \int \displaystyle \frac{{| \langle {\bf{o}},U| i\rangle | }^{2}{| \langle {\bf{o}},U| j\rangle | }^{2}}{{\displaystyle \sum }_{k=1}^{r}{| \langle {\bf{o}},U| k\rangle | }^{2}}{\mu }^{d}({\rm{d}}U)\\ & = & r\cdot d\displaystyle \int \displaystyle \frac{{| \langle {\bf{1}},U| i\rangle | }^{2}{| \langle {\bf{1}},U| j\rangle | }^{2}}{{\displaystyle \sum }_{k=1}^{r}{| \langle {\bf{1}},U| k\rangle | }^{2}}{\mu }^{d}({\rm{d}}U)\\ & = & r\cdot d\displaystyle \int \displaystyle \frac{{| \langle \psi | i\rangle | }^{2}{| \langle \psi | j\rangle | }^{2}}{{\displaystyle \sum }_{k=1}^{r}{| \langle \psi | k\rangle | }^{2}}{\nu }^{d}({\rm{d}}\psi ),\end{eqnarray*}$$

where ${\nu }^{d}({\rm{d}}\psi )$ is the uniform measure over the projective space on ${{\mathbb{C}}}^{d}.$ To compute the integral we decompose $| \psi \rangle$ as

$$\begin{eqnarray*}&&| \psi \rangle =q| {\psi }_{1}\rangle +\sqrt{1-{q}^{2}}| {\psi }_{2}\rangle \end{eqnarray*}$$

with $| {\psi }_{1}\rangle \in {{\mathcal{H}}}_{r}:= {\rm{Span}}\{| 1\rangle ,...,| r\rangle \}$ and $| {\psi }_{2}\rangle \in {{\mathcal{H}}}_{r}^{\perp }$ normalized (orthogonal) vectors, and $0\leqslant q\leqslant 1.$ The uniform measure can be expressed as

$$\begin{eqnarray*}&&{\nu }^{d}({\rm{d}}\psi )={m}^{d,r}({\rm{d}}q)\times {\nu }^{r}({\rm{d}}{\psi }_{1})\times {\nu }^{d-r}({\rm{d}}{\psi }_{2}),\end{eqnarray*}$$

where ${m}^{d,r}({\rm{d}}q)$ is the distribution of the length of the projection of a random vector in ${{\mathbb{C}}}^{d}$ onto an r-dimensional subspace. With this notation we have

$$\begin{eqnarray}&&{\bar{I}}_{{ii};{jj}}^{{dd}}=r\cdot d\int \int \int \displaystyle \frac{1}{{q}^{2}}| \langle \psi | i\rangle {| }^{2}| \langle \psi | j\rangle {| }^{2}{m}^{d,r}({\rm{d}}q){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{d-r}({\rm{d}}{\psi }_{2}).\end{eqnarray} \tag{ A.4 }$$

We distinguish four sub-cases depending on whether each of the indices $i,j$ belongs to $\{1,...,r\}$ or $\{r+1,...,d\}.$

Sub-case 1: $i,j\leqslant r.$ In this case (A.4) becomes

$$\begin{eqnarray}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\cdot d\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{q}^{4}{m}^{d,r}({\rm{d}}q)\times \displaystyle \int | \langle {\psi }_{1}| i\rangle {| }^{2}| \langle {\psi }_{1}| j\rangle {| }^{2}{\nu }^{r}({\rm{d}}{\psi }_{1})\\ & = & r\cdot d\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{q}^{4}{m}^{d,r}({\rm{d}}q)\times \displaystyle \int {U}_{1,i}^{*}{U}_{i,1}{U}_{1,j}^{*}{U}_{j,1}{\mu }^{r}({\rm{d}}U).\end{eqnarray} \tag{ A.5 }$$

In the last line we have re-written $| {\psi }_{1}\rangle$ as $U| 1\rangle$ in order to use the existing formulas for integrals of monomials over the unitary group [43]. Since we will use these formulas repeatedly, we recall that

$$\begin{eqnarray}&&\int {U}_{{i}_{1},{j}_{1}}{U}_{{i}_{2},{j}_{2}}{U}_{{j}_{1}^{\prime },{i}_{1}^{\prime }}^{*}{U}_{{j}_{2}^{\prime },{i}_{2}^{\prime }}^{*}{\mu }^{l}({\rm{d}}U)=\displaystyle \sum _{\sigma ,\tau \in {S}_{2}}{\delta }_{{i}_{1},{i}_{\sigma 1}^{\prime }}{\delta }_{{i}_{2},{i}_{\sigma 2}^{\prime }}{\delta }_{{j}_{1},{j}_{\tau 1}^{\prime }}{\delta }_{{j}_{2},{j}_{\tau 2}^{\prime }}{{\rm{Wg}}}^{l}(\sigma {\tau }^{-1}),\end{eqnarray} \tag{ A.6 }$$

where σ and τ are permutations in ${S}_{2}=\{(1,1),(2)\}$ and ${{\rm{Wg}}}^{l}(\cdot )$ is the Weingarten function over S₂:

$$\begin{eqnarray*}&&\ {{\rm{Wg}}}^{l}(1,1)=\displaystyle \frac{1}{{l}^{2}-1},\qquad {{\rm{Wg}}}^{l}(2)=\displaystyle \frac{-1}{l({l}^{2}-1)}.\end{eqnarray*}$$

The integral over q can be easily evaluated as

$$\begin{eqnarray}&&\int {q}^{2}{m}^{d,r}({\rm{d}}q)=\displaystyle \sum _{k=1}^{r}\int | \langle k| \psi \rangle {| }^{2}{\mu }^{d}({\rm{d}}\psi )=\displaystyle \frac{r}{d}.\end{eqnarray} \tag{ A.7 }$$

By inserting (A.6) and (A.7) into (A.5) we get

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\cdot d\cdot \displaystyle \frac{r}{d}\left(\displaystyle \frac{1}{{r}^{2}-1}-\displaystyle \frac{1}{r({r}^{2}-1)}\right)=\displaystyle \frac{r}{r+1},\quad i\ne j\\ {\bar{I}}_{{ii};{ii}}^{{dd}} & = & r\cdot d\cdot \displaystyle \frac{r}{d}\cdot 2\left(\displaystyle \frac{1}{{r}^{2}-1}-\displaystyle \frac{1}{r({r}^{2}-1)}\right)=\displaystyle \frac{2r}{r+1}.\end{eqnarray*}$$

Sub-case 2: $i\leqslant r$and$j\gt r.$ From (A.4) we get

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\cdot d\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{q}^{2}(1-{q}^{2}){m}^{d,r}({\rm{d}}q)\times \displaystyle \int | \langle {\psi }_{1}| i\rangle {| }^{2}{\nu }^{r}({\rm{d}}{\psi }_{1})\times \displaystyle \int | \langle {\psi }_{2}| j\rangle {| }^{2}{\nu }^{d-r}({\rm{d}}{\psi }_{2})\\ & = & r\cdot d\left(1-\displaystyle \frac{r}{d}\right)\cdot \displaystyle \frac{1}{r}\cdot \displaystyle \frac{1}{d-r}=1.\end{eqnarray*}$$

Sub-case 3: $i,j\gt r.$ Similarly, in this case

$$\begin{eqnarray}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\cdot d\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{(1-{q}^{2})}^{2}{m}^{d,r}({\rm{d}}q)\times \displaystyle \int | \langle {\psi }_{2}| i\rangle {| }^{2}| \langle {\psi }_{2}| j\rangle {| }^{2}{\nu }^{d-r}({\rm{d}}{\psi }_{2})\\ & = & r\cdot d\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{(1-{q}^{2})}^{2}{m}^{d,r}({\rm{d}}q)\times \displaystyle \int {U}_{1,i}^{*}{U}_{i,1}{U}_{1,j}^{*}{U}_{j,1}{\mu }^{d-r}({\rm{d}}U).\end{eqnarray} \tag{ A.8 }$$

We first simplify the integral on the right side

$$\begin{eqnarray}&&\int \displaystyle \frac{1}{{q}^{2}}{(1-{q}^{2})}^{2}{m}^{d,r}({\rm{d}}q)=\int \displaystyle \frac{1}{{q}^{2}}{m}^{d,r}({\rm{d}}q)-2+\displaystyle \frac{r}{d}.\end{eqnarray} \tag{ A.9 }$$

To evaluate the remaining integral, we consider a multivariate Gaussian random variable $c=({c}_{1},...,{c}_{2d})\sim N(0,{I}_{2d}),$ and denote by ${g}^{2d}({\rm{d}}c)$ its probability distribution. From this we construct the complex vector with uniform distribution over the unit ball in ${{\mathbb{C}}}^{d}$

$$\begin{eqnarray*}&&| \psi \rangle := \displaystyle \frac{1}{\parallel c\parallel }\displaystyle \sum _{i=1}^{d}({c}_{i}+{\rm{i}}{c}_{d+1})| i\rangle ,\qquad \parallel c{\parallel }^{2}=\displaystyle \sum _{l=1}^{2d}{c}_{i}^{2}.\end{eqnarray*}$$

We can now write

$$\begin{eqnarray}\displaystyle \int \displaystyle \frac{1}{{q}^{2}}{m}^{d,r}({\rm{d}}q) & = & \displaystyle \int {g}^{2d}({\rm{d}}c)\displaystyle \frac{{\displaystyle \sum }_{l=1}^{2d}{c}_{l}^{2}}{{\displaystyle \sum }_{l=1}^{2r}{c}_{l}^{2}}=1+\displaystyle \int {g}^{2d}({\rm{d}}c)\displaystyle \frac{{\displaystyle \sum }_{l=2r+1}^{2d}{c}_{l}^{2}}{{\displaystyle \sum }_{l=1}^{2r}{c}_{l}^{2}}\\ & = & 1+\displaystyle \int \parallel {c}^{1}{\parallel }^{2}\;{g}^{2(d-r)}({\rm{d}}{c}^{1})\cdot \displaystyle \int \displaystyle \frac{1}{\parallel {c}^{2}{\parallel }^{2}}{g}^{2r}({\rm{d}}{c}^{2})\\ & = & 1+2(d-r)\cdot \displaystyle \frac{1}{2r-2}=\displaystyle \frac{d-1}{r-1}.\end{eqnarray} \tag{ A.10 }$$

Above we used the fact that $\parallel {c}^{2}{\parallel }^{2}$ is a ${\chi }^{2}$ variable with $2r$ degrees of freedom, so the second integral is the mean of its inverse. By inserting (A.10) into (A.9) we obtain

$$\begin{eqnarray}&&\int \displaystyle \frac{1}{{q}^{2}}{(1-{q}^{2})}^{2}{m}^{d,r}({\rm{d}}q)=-2+\displaystyle \frac{r}{d}+\displaystyle \frac{d-1}{r-1}=\displaystyle \frac{(d-r)(d-r+1)}{d(r-1)}.\end{eqnarray} \tag{ A.11 }$$

Finally, by inserting (A.11) into (A.8) and applying (A.6) we obtain

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jj}}^{{dd}} & = & r\dot{d}\cdot \displaystyle \frac{(d-r)(d-r+1)}{d(r-1)}\cdot \left(\displaystyle \frac{1}{{(d-r)}^{2}-1}-\displaystyle \frac{1}{(d-r)({(d-r)}^{2}-1)}\right)\\ & = & \displaystyle \frac{r}{r-1},\quad i\ne j\\ {\bar{I}}_{{ii};{ii}}^{{dd}} & = & \displaystyle \frac{2r}{r-1}.\end{eqnarray*}$$

Note that for r = 1 these matrix elements are infinite. This is due to large contributions from measurements which have one basis vector close to being orthogonal to the one-dimensional vector state.. This is somewhat akin to what happens in the case of a Bernoulli variable (coin toss), in the case when the probability is close to zero or one. While this phenomenon is interesting, it does not play any role in our analysis for which the diagonal matrix elements are not relevant parameters.

(B) Diagonal-real block.

$$\begin{eqnarray*}&&I{(\rho | {{\bf{B}}}_{U})}_{{ii};{jk}}^{{dr}}=2\displaystyle \sum _{{\bf{o}}:{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})\gt 0}\displaystyle \frac{1}{{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}| \langle {\bf{o}},U| i\rangle {| }^{2}\cdot {\rm{Re}}(\langle j| {\bf{o}},U\rangle \langle {\bf{o}},U| k\rangle ).\end{eqnarray*}$$

As before, the corresponding entry of the average Fisher information matrix is

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jk}}^{{dr}} & = & 2\cdot d\cdot r\displaystyle \int \displaystyle \frac{{| \langle \psi | i\rangle | }^{2}\cdot {\rm{Re}}(\langle j| \psi \rangle \langle \psi | k\rangle )}{{\displaystyle \sum }_{k=1}^{r}{| \langle \psi | k\rangle | }^{2}}{\nu }^{d}({\rm{d}}\psi )\\ & = & 2\cdot d\cdot r\displaystyle \int \displaystyle \int \displaystyle \int \displaystyle \frac{1}{{q}^{2}}| \langle \psi | i\rangle {| }^{2}\cdot {\rm{Re}}(\langle j| \psi \rangle \langle \psi | k\rangle ){m}^{d,r}({\rm{d}}q){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{d-r}({\rm{d}}{\psi }_{2}).\end{eqnarray*}$$

Sub-case 1: $i,j,k\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jk}}^{{dr}} & = & 2\cdot d\cdot r\displaystyle \int {q}^{2}{m}^{d,r}({\rm{d}}q)\displaystyle \int {| \langle {\psi }_{1}| i\rangle | }^{2}{\rm{Re}}(\langle j| {\psi }_{1}\rangle \langle {\psi }_{1}| k\rangle ){\nu }^{r}({\rm{d}}{\psi }_{1})\\ & = & d\cdot r\cdot \displaystyle \frac{r}{d}\displaystyle \int {\mu }^{r}({\rm{d}}U){U}_{i1}{U}_{1i}^{*}({U}_{j1}{U}_{1k}^{*}+{U}_{k1}{U}_{1j}^{*})=0,\end{eqnarray*}$$

where we applied formula (A.6) for the integrals over the unitaries.

Sub-case 2: $i,j\leqslant r$ and $k\gt r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int q\sqrt{1-{q}^{2}}{m}^{d,r}({\rm{d}}q)\\ &&\times {\rm{Re}}\displaystyle \int | \langle {\psi }_{1}| i\rangle {| }^{2}\langle j| {\psi }_{1}\rangle {\nu }^{r}({\rm{d}}{\psi }_{1})\displaystyle \int \langle {\psi }_{2}| k\rangle {\nu }^{d-r}({\rm{d}}{\psi }_{2})=0.\end{eqnarray*}$$

Sub-case 3: $i\leqslant r$ and $j,k\gt r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int (1-{q}^{2}){m}^{d,r}({\rm{d}}q)\\ &&\times \displaystyle \int | \langle {\psi }_{1}| i\rangle {| }^{2}{\nu }^{r}({\rm{d}}{\psi }_{1})\cdot {\rm{Re}}\displaystyle \int \langle j| {\psi }_{2}\rangle \langle {\psi }_{2}| k\rangle {\nu }^{d-r}({\rm{d}}{\psi }_{2})=0.\end{eqnarray*}$$

Sub-case 4: $i\gt r$ and $j,k\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int (1-{q}^{2}){m}^{d,r}({\rm{d}}q)\\ &&\times \displaystyle \int | \langle {\psi }_{2}| i\rangle {| }^{2}{\nu }^{d-r}({\rm{d}}{\psi }_{2})\cdot {\rm{Re}}\displaystyle \int \langle j| {\psi }_{1}\rangle \langle {\psi }_{1}| k\rangle {\nu }^{r}({\rm{d}}{\psi }_{1})=0.\end{eqnarray*}$$

Sub-case 5: $i\gt r$ and $j,k\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int (1-{q}^{2}){m}^{d,r}({\rm{d}}q)\\ &&\times \displaystyle \int | \langle {\psi }_{2}| i\rangle {| }^{2}{\nu }^{d-r}({\rm{d}}{\psi }_{2})\cdot {\rm{Re}}\displaystyle \int \langle j| {\psi }_{1}\rangle \langle {\psi }_{1}| k\rangle {\nu }^{r}({\rm{d}}{\psi }_{1})=0.\end{eqnarray*}$$

Sub-case 6: $i,k\gt r$ and $j\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int q\sqrt{1-{q}^{2}}{m}^{d,r}({\rm{d}}q)\\ &&\times {\rm{Re}}\displaystyle \int | \langle {\psi }_{2}| i\rangle {| }^{2}\langle {\psi }_{2}| k\rangle {\nu }^{d-r}({\rm{d}}{\psi }_{2})\cdot \displaystyle \int \langle j| {\psi }_{1}\rangle {\nu }^{r}({\rm{d}}{\psi }_{1})=0.\end{eqnarray*}$$

Sub-case 7: $i,j,k\gt r.$ In this case the integral is

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{dr}}=2\cdot d\cdot r\displaystyle \int \displaystyle \frac{{(1-{q}^{2})}^{2}}{{q}^{2}}{m}^{d,r}({\rm{d}}q)\\ &&\times \displaystyle \int | \langle {\psi }_{2}| i\rangle {| }^{2}{\rm{Re}}\langle j| {\psi }_{2}\rangle \langle {\psi }_{2}| k\rangle {\nu }^{d-r}({\rm{d}}{\psi }_{2})=0.\end{eqnarray*}$$

The last integral is zero for the same reason as in sub-case 1.

In conclusion all diagonal-real matrix elements are zero

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};j,k}^{{dr}}=0,\qquad \mathrm{for}\;\mathrm{all}\;i=1,...,{2}^{k},\;\mathrm{and}\;1\leqslant j\lt k\leqslant {2}^{k}.\end{eqnarray*}$$

(C) Diagonal-imaginary block.

Similarly to the case of diagonal-real elements, we obtain that all diagonal-imaginary matrix elements are zero

$$\begin{eqnarray*}&&{\bar{I}}_{{ii};{jk}}^{{di}}=0,\qquad \mathrm{for}\;\mathrm{all}\;i=1,...,{2}^{k},\;\mathrm{and}\;1\leqslant j\lt k\leqslant {2}^{k}.\end{eqnarray*}$$

(D) Real–real block.

$$\begin{eqnarray*}&&I{(\rho | {{\bf{B}}}_{U})}_{{ij};{kl}}^{{rr}}=4\displaystyle \sum _{{\bf{o}}:{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})\gt 0}\displaystyle \frac{1}{{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\rm{Re}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle )\cdot {\rm{Re}}(\langle k| {\bf{o}},U\rangle \langle {\bf{o}},U| l\rangle ).\end{eqnarray*}$$

As before the corresponding entry of the average Fisher information matrix is

$$\begin{eqnarray*}{\bar{I}}_{{ij};{kl}}^{{rr}} & = & 4\cdot d\cdot r\displaystyle \int \displaystyle \frac{{\rm{Re}}(\langle i| \psi \rangle \langle \psi | j\rangle )\cdot {\rm{Re}}(\langle k| \psi \rangle \langle \psi | l\rangle )}{\displaystyle \sum _{k=1}^{r}{| \langle \psi | k\rangle | }^{2}}{\nu }^{d}({\rm{d}}\psi )\\ & = & 4\cdot d\cdot r\displaystyle \int \displaystyle \int \displaystyle \int \displaystyle \frac{1}{{q}^{2}}{\rm{Re}}(\langle i| \psi \rangle \langle \psi | j\rangle )\cdot {\rm{Re}}(\langle k| \psi \rangle \langle \psi | l\rangle ){m}^{d,r}({\rm{d}}q){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{d-r}({\rm{d}}{\psi }_{2}).\end{eqnarray*}$$

The same reasoning as before can be applied to transform the integral into a single integral, or a product of integrals over unitaries. If a single index is larger than r while the other three are smaller, or conversely a single index is smaller than r while the other are larger, then we have two integrals over unitaries which are zero since the monomial is of odd order. Therefore the following cases remain to be analysed.

Sub-case 1: $i,j,k,l\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}{\bar{I}}_{{ij};{kl}}^{{rr}} & = & 4\cdot d\cdot r\displaystyle \int {q}^{2}{m}^{d,r}({\rm{d}}q)\displaystyle \int {\rm{Re}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{1}| j\rangle )\cdot {\rm{Re}}(\langle k| {\psi }_{1}\rangle \langle {\psi }_{1}| l\rangle ){\nu }^{r}({\rm{d}}{\psi }_{1})\\ & = & d\cdot r\cdot \displaystyle \frac{r}{d}\cdot \displaystyle \int {\mu }^{r}({\rm{d}}U)({U}_{i1}{U}_{1j}^{*}+{U}_{j1}{U}_{1i}^{*})({U}_{k1}{U}_{1\;l}^{*}+{U}_{l1}{U}_{1k}^{*}).\end{eqnarray*}$$

Using formula (A.6) we find that the integral over unitaries is zero unless we deal with a diagonal matrix element of I, i.e. i = k and j = l. For the latter we have

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{rr}} & = & d\cdot r\cdot \displaystyle \frac{r}{d}\cdot \displaystyle \int {\mu }^{r}({\rm{d}}U)({U}_{i1}{U}_{1j}^{*}+{U}_{j1}{U}_{1i}^{*})({U}_{i1}{U}_{1j}^{*}+{U}_{j1}{U}_{1i}^{*})\\ & = & 2\cdot d\cdot r\cdot \displaystyle \frac{r}{d}\left(\displaystyle \frac{1}{{r}^{2}-1}-\displaystyle \frac{1}{r({r}^{2}-1)}\right)=\displaystyle \frac{2r}{r+1}.\end{eqnarray*}$$

Sub-case 2: $i,j,k,l\gt r.$ In this case, the integral is

$$\begin{eqnarray*}{\bar{I}}_{{ij};{kl}}^{{rr}} & = & 4\cdot d\cdot r\displaystyle \int \displaystyle \frac{{(1-{q}^{2})}^{2}}{{q}^{2}}{m}^{d,r}({\rm{d}}q)\displaystyle \int {\rm{Re}}(\langle i| {\psi }_{2}\rangle \langle {\psi }_{2}| j\rangle )\cdot {\rm{Re}}(\langle k| {\psi }_{2}\rangle \langle {\psi }_{2}| l\rangle ){\nu }^{d-r}({\rm{d}}{\psi }_{2})\\ & = & d\cdot r\cdot \displaystyle \frac{(d-r)(d-r+1)}{d(r-1)}\cdot \displaystyle \int {\mu }^{d-r}({\rm{d}}U)({U}_{i1}{U}_{1j}^{*}+{U}_{j1}{U}_{1i}^{*})({U}_{k1}{U}_{1\;l}^{*}+{U}_{l1}{U}_{1k}^{*}),\end{eqnarray*}$$

where we have used formula (A.11) for the integral over q. A similar calculation as above shows that all off-diagonal elements are zero and

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{rr}} & = & 2d\cdot r\cdot \displaystyle \frac{(d-r)(d-r+1)}{d(r-1)}\cdot \left(\displaystyle \frac{1}{{(d-r)}^{2}-1}-\displaystyle \frac{1}{(d-r)({(d-r)}^{2}-1)}\right)\\ & = & \displaystyle \frac{2r}{r-1}.\end{eqnarray*}$$

Sub-case 3: ($i,j\leqslant r$ and $k,l\gt r$) or ($i,j\gt r$ and $k,l\leqslant r$). In this case we deal with a product of integrals over unitaries of the form

$$\begin{eqnarray*}&&\int {\mu }^{r}({\rm{d}}U){U}_{i1}{U}_{1j}^{*}=0.\end{eqnarray*}$$

So all matrix elements are zero.

Sub-case 4: $i,k\leqslant r$ and $j,l\gt r.$ Again, the off-diagonal elements are zero and

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{rr}}= & = & 4\cdot d\cdot r\displaystyle \int (1-{q}^{2}){m}^{d,r}({\rm{d}}q)\displaystyle \int {\rm{Re}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{2}| j\rangle )\cdot {\rm{Re}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{2}| j\rangle ){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{r}({\rm{d}}{\psi }_{2})\\ & = & 2d\cdot r\cdot (1-\displaystyle \frac{r}{d})\cdot \displaystyle \frac{1}{r}\displaystyle \frac{1}{d-r}=2.\end{eqnarray*}$$

In conclusion, the only non-zero elements of the real–real block are on the diagonal and are given by

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{rr}} & = & \displaystyle \frac{2r}{r+1},\qquad 1\leqslant i\lt j\leqslant r,\\ {\bar{I}}_{{ij};{ij}}^{{rr}} & = & 2,\qquad \qquad \;i\leqslant r\;\mathrm{and}\ \;r\lt j\leqslant {2}^{k},\\ {\bar{I}}_{{ij};{ij}}^{{rr}} & = & \displaystyle \frac{2r}{r-1},\qquad \;r\lt i\lt j\leqslant {2}^{k}.\end{eqnarray*}$$

(E) Imaginary–imaginary block. This block is similar to the real–real one. All off diagonal elements are zero, and the diagonal ones are

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{ii}} & = & \displaystyle \frac{2r}{r+1},\qquad 1\leqslant i\lt j\leqslant r,\\ {\bar{I}}_{{ij};{ij}}^{{ii}} & = & 2,\qquad \qquad \;i\leqslant r\;\mathrm{and}\ \;r\lt j\leqslant {2}^{k},\\ {\bar{I}}_{{ij};{ij}}^{{ii}} & = & \displaystyle \frac{2r}{r-1},\qquad \;r\lt i\lt j\leqslant {2}^{k}.\end{eqnarray*}$$

(F) Real-imaginary block. Next we show that all real-imaginary off-diagonal elements are equal to zero.

$$\begin{eqnarray*}&&I{(\rho | {{\bf{B}}}_{U})}_{{ij};{kl}}^{{ri}}=-4\cdot \displaystyle \sum _{{\bf{o}}:{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})\gt 0}\displaystyle \frac{1}{{p}_{\rho }({\bf{o}}| {{\bf{B}}}_{U})}{\rm{Re}}(\langle i| {\bf{o}},U\rangle \langle {\bf{o}},U| j\rangle ){\rm{Im}}(\langle k| {\bf{o}},U\rangle \langle {\bf{o}},U| l\rangle ).\end{eqnarray*}$$

By integration we obtain the corresponding matrix element of the average Fisher information matrix

$$\begin{eqnarray*}{\bar{I}}_{{ij};{kl}}^{{rr}} & = & -4\cdot d\cdot r\displaystyle \int \displaystyle \frac{{\rm{Re}}(\langle i| \psi \rangle \langle \psi | j\rangle )\cdot {\rm{Im}}(\langle k| \psi \rangle \langle \psi | l\rangle )}{{\displaystyle \sum }_{k=1}^{r}{| \langle \psi | k\rangle | }^{2}}{\nu }^{d}({\rm{d}}\psi )\\ & = & -4\cdot d\cdot r\displaystyle \int \displaystyle \int \displaystyle \int \displaystyle \frac{1}{{q}^{2}}{\rm{Re}}(\langle i| \psi \rangle \langle \psi | j\rangle )\cdot {\rm{Im}}(\langle k| \psi \rangle \langle \psi | l\rangle ){m}^{d,r}({\rm{d}}q){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{d-r}({\rm{d}}{\psi }_{2}).\end{eqnarray*}$$

Again, if one index is smaller that r while the other three are larger, or otherwise, then the matrix element in zero. We analyse the remaining cases.

Sub-case 1: $i,j,k,l\leqslant r.$ In this case the integral is

$$\begin{eqnarray*}{\bar{I}}_{{ij};{kl}}^{{ri}} & = & -4\cdot d\cdot r\displaystyle \int {q}^{2}{m}^{d,r}({\rm{d}}q)\displaystyle \int {\rm{Re}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{1}| j\rangle ){\rm{Im}}(\langle k| {\psi }_{1}\rangle \langle {\psi }_{1}| l\rangle ){\nu }^{r}({\rm{d}}{\psi }_{1})\\ & = & i\cdot d\cdot r\cdot \displaystyle \frac{r}{d}\displaystyle \int {\mu }^{r}({\rm{d}}U)({U}_{i1}{U}_{1j}^{*}+{U}_{j1}{U}_{1i}^{*})({U}_{k1}{U}_{1\;l}^{*}-{U}_{l1}{U}_{1k}^{*}).\end{eqnarray*}$$

Using formula (A.6) we find that the integral over unitaries is zero unless i = k and j = l. However even in this case, two of the four terms are zero, and the other two cancel each other.

Sub-case 2: $i,j,k,l\gt r.$ Here the integrals over the unitaries are similar as in case 1 above, with the difference that they taken with respect to ${\mu }^{d-r}({\rm{d}}U).$ Therefore the matrix elements are zero.

Sub-case 3: ($i,j\leqslant r$ and $k,l\gt r$) or ($i,j\gt r,k,l\leqslant r$). In this case we deal with a product of integrals over unitaries of the form

$$\begin{eqnarray*}&&\int {\mu }^{r}({\rm{d}}U){U}_{i1}{U}_{1j}^{*}=0.\end{eqnarray*}$$

so all matrix elements are zero.

Sub-case 4: $i,k\leqslant r$ and $j,l\gt r$ Again, the off-diagonal elements are zero and

$$\begin{eqnarray*}{\bar{I}}_{{ij};{ij}}^{{ri}}= & = & -4\cdot d\cdot r\displaystyle \int (1-{q}^{2}){m}^{d,r}({\rm{d}}q)\displaystyle \int {\rm{Re}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{2}| j\rangle )\cdot {\rm{Im}}(\langle i| {\psi }_{1}\rangle \langle {\psi }_{2}| j\rangle ){\nu }^{r}({\rm{d}}{\psi }_{1}){\nu }^{r}({\rm{d}}{\psi }_{2})\\ & = & i\cdot d\cdot r\cdot \left(1-\displaystyle \frac{r}{d}\right)\displaystyle \int \displaystyle \int \left({U}_{i1}{V}_{1j}^{*}+{V}_{j1}{U}_{1i}^{*}\right)\left({U}_{i1}{V}_{1j}^{*}-{V}_{j1}{U}_{1i}^{*}\right){\mu }^{r}({\rm{d}}U){\mu }^{d-r}({\rm{d}}V)\\ & = & 0.\end{eqnarray*}$$

In the last integral, two terms are zero and two have different signs and cancel each other. In conclusion all elements of the real-imaginary block are equal to zero.

Summary of the computation of$\bar{I}$. We found that all off-diagonal blocks of $\bar{I}$ are zero

$$\begin{eqnarray*}&&{\bar{I}}^{{dr}}={\bar{I}}^{{di}}={\bar{I}}^{{ri}}={\bar{I}}^{{rd}}={\bar{I}}^{{id}}={\bar{I}}^{{dr}}=0.\end{eqnarray*}$$

Moreover the real and imaginary diagonal blocks are diagonal, equal to each other and have three distinct values depending on the position of the indices $i\lt j$ with respect to r:

$$\begin{eqnarray*}{\bar{I}}^{{rr}/{ii}}={\rm{Diag}}\left\{\begin{array}{cccc}{\bar{I}}_{{ij};{ij}}^{{rr}/{ii}} & = & \displaystyle \frac{2r}{r+1}, & 1\leqslant i,j\leqslant r\\ {\bar{I}}_{{ij};{ij}}^{{rr}/{ii}} & = & 2, & 1\leqslant i\leqslant r,{\rm{a}}{nd}r\lt j\leqslant {2}^{k}\\ {\bar{I}}_{{ij};{ij}}^{{rr}/{ii}} & = & \displaystyle \frac{2r}{r-1} & r\lt i\lt j\leqslant {2}^{k}.\end{array}\right.\end{eqnarray*}$$

Finally, the d × d block ${\bar{I}}^{{dd}}$ is not diagonal but has a simple form

$$\begin{eqnarray*}{\bar{I}}_{{ii};{jj}}^{{dd}}=\left\{\begin{array}{cccc}{\bar{I}}_{{ii};{ii}}^{{dd}} & = & \displaystyle \frac{2r}{r+1}, & 1\leqslant i\leqslant r\\ {\bar{I}}_{{ii};{ii}}^{{dd}} & = & \displaystyle \frac{2r}{r-1}, & r\lt i\leqslant {2}^{k}\\ {\bar{I}}_{{ii};{jj}}^{{dd}} & = & \displaystyle \frac{r}{r+1}, & i,j\leqslant r,{\rm{a}}{nd}i\ne j\\ {\bar{I}}_{{ii};{jj}}^{{dd}} & = & \displaystyle \frac{r}{r+1}, & 1\leqslant i\leqslant r{\rm{a}}{nd}r\lt j\leqslant {2}^{k}\\ {\bar{I}}_{{ii};{jj}}^{{dd}} & = & \displaystyle \frac{r}{r+1}, & 1\leqslant j\leqslant r{\rm{a}}{nd}r\lt i\leqslant {2}^{k}\\ {\bar{I}}_{{ii};{jj}}^{{dd}} & = & \displaystyle \frac{r}{r-1} & r\lt i,j,{\rm{a}}{nd}i\ne j.\end{array}\right.\end{eqnarray*}$$

The model ${{\mathcal{R}}}_{d,r}$ can be seen (locally) as the restriction of the full unconstrained model ${{\mathcal{S}}}_{r,d}$ parametrized by θ, to the subset of parameters $\tilde{\theta }$ which are real and imaginary parts of matrix elements ${\rho }_{i,j}$ with $i\leqslant r\lt j.$ Therefore the corresponding average Fisher information matrix is equal to the corresponding block of $\bar{I},$ i.e. $\bar{\tilde{I}}={21}_{2r(d-r)}.$□