A comparative study of estimation methods in quantum tomography

TL;DR

This paper compares various quantum state estimation methods, analyzing their theoretical properties and practical performance, revealing insights into their asymptotic behavior, risk reduction, and computational efficiency.

Contribution

It provides a comprehensive comparison of multiple quantum tomography estimators, including new asymptotic analyses and simulation results highlighting their advantages and limitations.

Findings

Projected least squares estimator shows rank-dependent risk reduction.

Non-monotonic behavior of Bures risk observed in asymptotic analysis.

Projected estimators demonstrate high computational efficiency and strong optimality.

Abstract

As quantum tomography is becoming a key component of the quantum engineering toolbox, there is a need for a deeper understanding of the multitude of estimation methods available. Here we investigate and compare several such methods: maximum likelihood, least squares, generalised least squares, positive least squares, thresholded least squares and projected least squares. The common thread of the analysis is that each estimator projects the measurement data onto a parameter space with respect to a specific metric, thus allowing us to study the relationships between different estimators. The asymptotic behaviour of the least squares and the projected least squares estimators is studied in detail for the case of the covariant measurement and a family of states of varying ranks. This gives insight into the rank-dependent risk reduction for the projected estimator, and uncovers an interesting non-monotonic behaviour of the Bures risk. These asymptotic results complement recent non-asymptotic concentration bounds of \cite{GutaKahnKungTropp} which point to strong optimality properties, and high computational efficiency of the projected linear estimators. To illustrate the theoretical methods we present results of an extensive simulation study. An app running the different estimators has been made available online.

Figures

Tables

Table 1. Table 1: The different error functions used to measure the distance between the true state and the estimator. The Bures distance is defined only for states ρ ^ , ρ ∈ 𝒮 d ^ 𝜌 𝜌 subscript 𝒮 𝑑 \hat{\rho},\rho\in\mathcal{S}_{d} , and its classical analogue the Hellinger distance is defined between two probability distributions.

Error function	Definition
Frobenius norm squared	${\Vert \hat{\rho }-\rho \Vert }_2^2=\mathrm{Tr}\left[{\left(\hat{\rho }-\rho \right)}^2\right]$
Trace norm	${\Vert \hat{\rho }-\rho \Vert }_1=\mathrm{Tr}\left|\hat{\rho }-\rho \right|$
Operator norm	$\Vert \rho -\sigma \Vert ={\lambda }_{\max }\left(\left|\hat{\rho }-\rho \right|\right)$
Bures distance	$D_B{(\hat{\rho },\rho )}^2:=2\left(1-\mathrm{Tr}\left[\sqrt{\sqrt{\rho }\hat{\rho }\sqrt{\rho }}\right]\right)$
Hellinger distance	$D_H{(\widehat{𝝀},𝝀)}^2:=2\left(1-{\sum }_i^d\sqrt{{\hat{\lambda }}_i{\lambda }_i}\right)$

Equations

$$\mathbf{f}=X\mathbf{\beta }+\mathbf{ϵ}$$

$$P_{\mathbf{o}}^{\mathbf{s}}=|e_{o_1}^{s_1}⟩⟨e_{o_1}^{s_1}|\otimes \dots \otimes |e_{o_n}^{s_n}⟩⟨e_{o_n}^{s_n}|,{\sigma }_s|e_o^s⟩=o|e_o^s⟩.$$

$$p_{\rho }(\mathcal{D}|\mathcal{S})=\prod _{\mathbf{s}}\frac{m!}{{\prod }_{\mathbf{o}}N(\mathbf{o}|\mathbf{s})!}\prod _{\mathbf{o}}{p}_{\rho }(\mathbf{o}|\mathbf{s}{)}^{N(\mathbf{o}|\mathbf{s})}$$

$$\mathcal{M}:M({\mathbb{C}}^d)$$

$$\mathcal{M}:\rho$$

$${\hat{\mathbf{\theta }}}_{\mathrm{ML}}:=\underset{{\mathbf{\theta }}^{\prime }\in \Theta }{\mathrm{argmax}}{p}_{{\mathbf{\theta }}^{\prime }}(X_1,\dots ,X_N).$$

$${\ell }_{\mathbf{\theta }}$$

$${\hat{\mathbf{\theta }}}_{\mathrm{ML}}=\underset{{\mathbf{\theta }}^{\prime }\in \Theta }{\mathrm{argmin}}K(\mathbf{f}\Vert p_{{\mathbf{\theta }}^{\prime }})$$

$$\sqrt{N}({\hat{\mathbf{\theta }}}_{\mathrm{ML}}-\mathbf{\theta })\stackrel{\mathcal{D}}{⟶}N(0,I^{-1}(\mathbf{\theta })).$$

$$I(\mathbf{\theta }{)}_{i,j}={\mathbb{E}}_{\mathbf{\theta }}\left[\frac{\partial \log p_{\mathbf{\theta }}(X)}{\partial {\mathbf{\theta }}_i}\frac{\partial \log p_{\mathbf{\theta }}(X)}{\partial {\mathbf{\theta }}_j}\right]$$

$$d(\hat{\mathbf{\theta }},\mathbf{\theta })=(\hat{\mathbf{\theta }}-\mathbf{\theta }{)}^{T}G(\mathbf{\theta })(\hat{\mathbf{\theta }}-\mathbf{\theta })+o(\Vert \hat{\mathbf{\theta }}-\mathbf{\theta }{\Vert }^2)$$

$$N\mathbb{E}[d({\hat{\mathbf{\theta }}}_{\mathrm{ML}},\mathbf{\theta })]⟶\mathrm{Tr}(G(\mathbf{\theta })I^{-1}(\mathbf{\theta })).$$

$${\hat{\rho }}_{\mathrm{ML}}:=\underset{\tau \in {\mathcal{S}}_d}{\mathrm{argmax}}\sum _{\mathbf{o},\mathbf{s}}\mathbf{f}(\mathbf{o}|\mathbf{s})p_{\tau }(\mathbf{o}|\mathbf{s})=\underset{\tau \in {\mathcal{S}}_d}{\mathrm{argmin}}\sum _{\mathbf{s}}K(\mathbf{f}(\cdot |\mathbf{s})\Vert p_{\tau }(\cdot |\mathbf{s}))$$

$${\rho }^{\prime }={\rho }_{\mathbf{\theta }}=\rho +\delta {\rho }_{\mathbf{\theta }}$$

$$\mathbf{\theta }$$

$$=(\mathrm{Re}\delta {\rho }_{1,2}^{\prime },\dots ,\mathrm{Re}\delta {\rho }_{d-1,d}^{\prime };\mathrm{Im}\delta {\rho }_{1,2}^{\prime },\dots ,\mathrm{Im}\delta {\rho }_{d-1,d}^{\prime };\delta {\rho }_{2,2}^{\prime },\dots ,\delta {\rho }_{d,d}^{\prime })\in {\mathbb{R}}^{d^2-1}.$$

$$I(\rho |\mathcal{S})=\frac{1}{k}\sum _{\mathbf{s}\in \mathcal{S}}I(\rho |\mathbf{s}),I(\rho |\mathbf{s}{)}_{a,b}:=\sum _{\mathbf{o}:p(\mathbf{o}|\mathbf{s})>0}\frac{1}{p_{\rho }(\mathbf{o}|\mathbf{s})}\frac{\partial p_{\rho }(\mathbf{o}|\mathbf{s})}{\partial {\theta }_a}\cdot \frac{\partial p_{\rho }(\mathbf{o}|\mathbf{s})}{\partial {\theta }_b}.$$

$$N\mathbb{E}\left[d({\rho }_{\mathbf{\theta }},{\hat{\rho }}_{\mathrm{ML}})\right]⟶\mathrm{Tr}\left(I(\rho |\mathcal{S}{)}^{-1}G(\mathbf{\theta })\right).$$

$$\rho^{\prime}=\left(\begin{array}[]{cc}D&0\\ 0&0\end{array}\right)+\left(\begin{array}[]{cc}A&B\\ B^{*}&C\end{array}\right),\qquad D:={\rm Diag}(\lambda_{1},\dots,\lambda_{r})$$

$$\ell (\mathbf{\theta })-\ell ({\mathbf{\theta }}_{\mathrm{uML}})\approx -\frac{N}{2}(\mathbf{\theta }-{\hat{\mathbf{\theta }}}_{\mathrm{uML}}{)}^{T}I(\rho |\mathcal{S})(\mathbf{\theta }-{\hat{\mathbf{\theta }}}_{\mathrm{uML}}).$$

$$d_I({\rho }_{\mathbf{\theta }},{\rho }_{{\mathbf{\theta }}^{\prime }}):=(\mathbf{\theta }-{\mathbf{\theta }}^{\prime }{)}^{T}I(\rho |\mathcal{S})(\mathbf{\theta }-{\mathbf{\theta }}^{\prime }).$$

$$\rho =\sum _{i=1}^{d^2}{\beta }_i{\tau }_i$$

$$p_{\rho }(\mathbf{o}|\mathbf{s})=\mathrm{Tr}(\rho P_{\mathbf{o}}^{\mathbf{s}})=\sum _i{X}_{(\mathbf{o}|\mathbf{s})}^i{\beta }_i,X_{(\mathbf{o}|\mathbf{s})}^i:=\mathrm{Tr}({\tau }_i{P}_{\mathbf{o}}^{\mathbf{s}})$$

$$\mathbf{f}=X\mathbf{\beta }+\mathbf{ϵ},$$

$${\hat{\mathbf{\beta }}}_{\mathrm{LS}}:=\underset{v\in {\mathbb{R}}^{d^2}}{\mathrm{argmin}}\Vert Xv-\mathbf{f}{\Vert }^2$$

$$\mathcal{D}:\rho \mapsto \frac{1}{|\mathcal{S}|}\sum _{\mathbf{o},\mathbf{s}}\mathrm{Tr}(\rho P_{\mathbf{o}}^{\mathbf{s}})\cdot P_{\mathbf{o}}^{\mathbf{s}}$$

$${\hat{\rho }}_{\mathrm{LS}}={\mathcal{D}}^{-1}\circ \mathcal{P}(\mathbf{f}/|\mathcal{S}|).$$

$$\mathcal{C}:\rho \mapsto \frac{1}{3}\rho +\frac{2}{3}\frac{1}{2}.$$

$$\mathcal{D}:\rho \mapsto \frac{1}{d+1}\rho +\frac{d}{d+1}\frac{1}{d}.$$

$${\hat{\rho }}_{\mathrm{LS}}=\frac{1}{m}\sum _{i=1}^m{A}_i=\frac{1}{m}\sum _{i=1}^m\frac{1}{|\mathcal{S}|}\sum _{\mathbf{s}\in \mathcal{S}}{\mathcal{D}}^{-1}(P_{{\mathbf{o}}_{i|\mathbf{s}}}^{\mathbf{s}})$$

$$\mathrm{Pr}\left[\Vert {\hat{\rho }}_{\mathrm{LS}}-\rho \Vert \ge ϵ\right]\le d{\mathrm{e}}^{-\frac{3N{ϵ}^2}{8g(d)}}ϵ\in \left[0,1\right].$$

Full text

The asymptotic behaviour of the least squares and the projected least squares estimators is studied in detail for the case of the covariant measurement and a family of states of varying ranks. This gives insight into the rank-dependent risk reduction for the projected estimator, and uncovers an interesting non-monotonic behaviour of the Bures risk. These asymptotic results complement recent non-asymptotic concentration bounds of [36] which point to strong optimality properties, and high computational efficiency of the projected linear estimators.

To illustrate the theoretical methods we present results of an extensive simulation study. An app running the different estimators has been made available online.

1 Introduction

The problem of estimating unknown parameters of quantum systems has been at the forefront of early investigations into the statistical nature of quantum information [13, 43, 42, 79, 16]. Traditionally, key research topics have been the design of optimal measurements and estimation procedures [56, 77, 21, 31, 47, 8, 9, 41, 9], and theoretical aspects of quantum Fisher information and asymptotic estimation [11, 28, 57, 41, 6, 62, 45, 30], see also the monographs [59, 52, 40].

More recently, quantum tomography has become a crucial validation tool current quantum technology applications [38, 58, 29, 67]. The experimental challenges have stimulated research in new directions such as compressed sensing [35, 54, 27, 70, 20, 5, 3], estimation of permutationally invariant states [73], adaptive and selflearning tomography [55, 34, 60, 26, 39, 63], incomplete tomography [74], minimax bounds [25, 4], Bayesian estimation [15, 33], and confidence regions [7, 18, 71, 24, 53]. Since ‘full tomography’ becomes impossible for systems composed of even a moderate number of qubits, research has focused on the estimation of states which belong to smaller dimensional models which nevertheless capture relevant physical properties, such as low rank states [48, 17, 37, 49, 36] and matrix product states [20, 12, 50].

In this paper we analyse and compare several estimation methods for fixed (non-adaptive) measurement designs, with a focus on risks (mean errors) for different loss functions, asymptotic theory, relationships between estimators and low rank behaviour. The measurement scenarios include repeated measurements with respect to Pauli bases (as customary in multiple qubits experiments), random bases measurements, and the covariant measurement. The loss functions are given by the Frobenius, trace-norm, operator-norm, Bures and Hellinger distances. Each section deals with one class of estimators and the results of a comparative simulations study are presented at the end of the paper. While most estimators have been previously considered in the literature, our aim is to investigate them from a common perspective, as projections of data onto the parameter space. Another aim is to understand and quantify the reduction in statistical error between an unconstrained estimator such as least squares, and an estimator which take into account the physical properties of the parameter space, such as the projected least squares estimator. Among the original results, we derive the asymptotic error rates of these estimators on a class of low rank states in the covariant measurement scenario. Finally, we discuss the computational efficiency of different methods.

Below we summarise the paper using Figure 1 for illustration. A measurement $\mathcal{M}$ on $d$-dimensional system is a positive affine map from the convex set of states $\mathcal{S}_{d}$ (constrained parameter space) onto the space of outcome probabilities $\mathcal{P}_{d}\subset\mathbb{R}^{z}$, where $z$ is the number outcomes. The image of the set of trace-one selfadjoint matrices $M^{1}_{sa}(\mathbb{C}^{d})$ (unconstrained parameter space) is a hyperplane $\mathcal{L}_{d}$ of $\mathbb{R}^{z}$, which contains $\mathcal{P}_{d}$. For informationally complete measurements the map $\mathcal{M}$ is injective and we can identify all matrix estimators (whether positive or not) with their images in $\mathcal{L}_{d}$.

In section 3 we review the use of Fisher information in asymptotic normality theory and discuss to what extent it is applicable to the study of the maximum likelihood (ML) estimator in quantum tomography [10, 14, 66, 75]. We distinguish between two ML estimators: the unconstrained estimator $\hat{\rho}_{\rm uML}$, and the ‘physical’ estimator $\hat{\rho}_{\rm ML}$. The former is the maximiser of the likelihood, seen as a function over $M^{1}_{sa}(\mathbb{C}^{d})$, and may not be a density matrix. The latter performs the same optimisation but restricted to the physical space of density matrices $\mathcal{S}_{d}\subset M^{1}_{sa}(\mathbb{C}^{d})$. For large sample size, the unconstrained estimator is asymptotically normal, while the constrained estimator is equal to its projection onto $\mathcal{S}_{d}$ with respect to the metric defined by the classical Fisher information [66].

In section 4 we analyse the least squares (LS) estimator $\hat{\rho}_{\rm LS}$. This exploits the linear dependence between probabilities and states which translates into a linear regression problem

[TABLE]

where $\boldsymbol{f}$ denotes the vector of observed frequencies of measurement outcomes, $X$ is a fixed measurement design matrix, $\boldsymbol{\beta}$ is a vectorisation of the unknown state $\rho$, and $\boldsymbol{\epsilon}$ is the statistical noise. The LS estimator minimises the prediction error $\|X\hat{\boldsymbol{\beta}}-\boldsymbol{f}\|^{2}$ over all vectors $\hat{\boldsymbol{\beta}}$, and is the simplest and fastest estimator to compute. However, LS is significantly less accurate than ML, especially for low rank states and does not produce a physical state. Nevertheless, the LS has been the focus of recent investigations [17, 72, 36] as a first step towards more accurate estimators. Here we review some of the non-asymptotic concentration bounds for the operator-norm error of the LS estimator. We then study the asymptotic properties of the LS estimator in the context of covariant measurements, and rank $r$ states $\rho_{r}$ with equal non-zero eigenvalues. By exploiting the symmetry of the measurement we compute the explicit expression of the risk with respect to the Frobenius distance. Furthermore, we show that for large $d$, the eigenvalues of the ‘error matrix’ $\hat{\rho}_{\rm LS}-\rho_{r}$ are approximately distributed according to the Wigner semicircle law on the interval $[-2\sqrt{d/N},2\sqrt{d/N}]$, cf. Figure 2. This provides the asymptotic estimates of the operator-norm, and the trace-norm errors of $2\sqrt{d/N}$ and respectively $8d^{3/2}/(3\pi\sqrt{N})$, which complement the non-asymptotic bounds of [36].

In section 5 we discuss the generalised least squares (GLS) estimator. Unlike the LS estimator which minimises the prediction error, the GLS aims to optimise the estimation error, e.g. $\mathbb{E}\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|^{2}$ and needs to take into account the covariance matrix of the multinomial noise $\boldsymbol{\epsilon}$. We show that for large samplesize $N$, the GLS estimator is asymptotically normal and is equivalent to the uML estimator.

Section 6 reviews the thresholded least squares (TLS) estimator proposed in [17]. This is obtained by projecting the LS estimator onto the set of states whose non-zero eigenvalues are above a certain ‘noise level’ chosen in a data-driven fashion. The projection can be computed efficiently in terms of the eigenvalues of the LS estimator, and the truncation leads to significant error reduction for low rank states. In practice, an additional improvement is achieved by using the GLS estimator as starting point.

Section 7 discusses the positive least squares (posLS) estimator which optimises the prediction error over the physical space of density matrices, rather than over the selfadjoint matrices as is the case for the LS estimator [46]. This leads to higher accuracy, but its computational complexity is similar to that of ML. However, in the case of the covariant measurements we find that posLS is equivalent to the projected least squares estimator discussed below, which can be computed efficiently. By restricting the GLS optimisation over density matrices we obtain the positive generalised least squares estimator (posGLS), which is shown to be asymptotically equivalent to the ML estimator.

Section 8 deals with the projected least squares (PLS) estimator. This is obtained by projecting the LS estimator onto the space of density matrices with respect to the Frobenius distance [72, 36]. It is faster than TLS, and has similar statistical properties. In [36] is was shown that the PLS estimator satisfies rank-dependent concentration bounds and its trace-norm error scales as $O(r^{2}\cdot d\cdot\log d/\sqrt{N})$ for a broad class of 2-design measurements including mutually unbiased bases, stabiliser states, and symmetric informationally complete measurements, and as $O(r^{2}\cdot d^{1.6}\cdot\log d/\sqrt{N})$ for Pauli bases measurements. In this paper we focus on the asymptotic behaviour of the PLS estimator in the case of covariant measurements. Inspired by the techniques developed in [66] we show how the random matrix properties of the LS can be used to derive the asymptotic Frobenius and Bures risks of the PLS estimator for low rank states and large dimension $d$ and samplesize $N$. In particular we uncover an interesting behaviour of the Bures risk which (unlike the Frobenius or the trace norm risks) increases steeply with rank for low rank states and then decreases towards full rank, cf. Figure 3.

Section 9 presents the results of a comparative simulations study of the proposed estimators. We simulated data for a range of states of 3 and 4 atoms of different ranks, with different samplesizes, and measurement setups. For each choice we produced 100 datasets in order to estimate risks of different estimators. The measurements are chosen to be either Pauli bases measurements (as standard in ion trap experiments) or random basis measurements. The different estimators and their risks corresponding to Frobenius, Bures, Hellinger, and trace-norm distances were computed and the results are illustrated and discussed, cf. Figures 4, 5, 6, 7, 8, 9, 10, and 11. A complete set of simulation results available online via an interactive Rshiny application at https://rudhacharya.shinyapps.io/plots/. We also created an online estimation app which computes the key estimators for a range of states of certain characteristics, or user uploaded states. This is available at https://shiny.maths.nottingham.ac.uk/shiny/qt_dashboard/.

2 Quantum tomography

In quantum tomography, the aim is to estimate an unknown state from outcomes of measurements performed on an ensemble of identical prepared systems. Although a large part of our theoretical considerations hold in a general setting, we choose to formulate the tomography problem in the specific context of a system consisting of multiple qubits and projective measurements with respect to several bases. This keeps the discussion closer to realistic experimental procedures, and facilitates the understanding of our simulation results.

We consider composite systems consisting of $n$ qubits, with associated Hilbert space $\mathbb{C}^{d}$, where $d=2^{n}$. The state $\rho$ belongs to the space of $d\times d$ density matrices $\mathcal{S}_{d}\subset M(\mathbb{C}^{d})$. In our analysis we will distinguish between the constrained parameter space $\mathcal{S}_{d}$ whose elements are trace-one positive matrices, and the unconstrained parameter space $M^{1}_{sa}(\mathbb{C}^{d})$ consisting of trace-one selfadjoint matrices. This will allow us to consider procedure which produce constrained, or unconstrained estimators.

The measurement strategies consist of performing standard, von Neumann projective measurements with respect to a number of orthonormal bases (ONB), which are chosen deterministically or randomly. In particular, we focus on two scenarios: Pauli basis measurements and measurements that are drawn randomly from the uniform measure over all ONBs. While the former setup is commonly employed in experiments [38], the latter is less restrictive, more amenable to theoretical analysis, and can serve as a benchmark for current experiments. We also consider covariant measurements in the context of the asymptotic theory of the least squares and respectively projected least squares estimators in sections 4.3 and 8.1, and refer to these for further details.

In the Pauli bases setup, one measures an arbitrary Pauli observable $\sigma_{x},\sigma_{y}$, or $\sigma_{z}$ on each of the $n$ ions simultaneously. Therefore, each measurement is labelled by a sequence $\boldsymbol{s}=(s_{1},\ldots,s_{n})\in\{x,y,z\}^{n}$, and there are $3^{n}$ possible measurement bases. Such a measurement produces a $\pm 1$ outcome from each ion, and we let $\boldsymbol{o}\in\{+1,-1\}^{n}$ be the full record of outcomes from all the $n$ ions. The probability of obtaining a particular outcome $\boldsymbol{o}$ is given by $p_{\rho}(\boldsymbol{o}|\boldsymbol{s}):={\rm Tr}(\rho P_{\boldsymbol{o}}^{\boldsymbol{s}})$, where $P_{\boldsymbol{o}}^{\boldsymbol{s}}$ is the one-dimensional projection onto the tensor product of eigenvectors of Pauli matrices

[TABLE]

More generally, the measurement design is defined by a collection $\mathscr{S}=\{\boldsymbol{s}_{1},\ldots,\boldsymbol{s}_{k}\}$ of ONBs, which may be chosen deterministically or randomly. For each setting $\boldsymbol{s}$, independent measurements are performed on $m$ identical copies of the state, and the counts $N(\boldsymbol{o}|\boldsymbol{s})$ are recorded, where $N(\boldsymbol{o}|\boldsymbol{s})$ is the number of times the outcome $\boldsymbol{o}$ is observed when measuring in setting $\boldsymbol{s}$. The total number of quantum samples is therefore $N=m\times k$. The resulting dataset $\mathcal{D}:=\{N(\boldsymbol{o}|\boldsymbol{s}):\boldsymbol{o}\in\{+1,-1\}^{n},\boldsymbol{s}\in\mathscr{S}\}$ is a $2^{n}\times k$ table whose columns are independent and contain all the counts in a given setting. Its probability is given by the product of multinomials

[TABLE]

Our goal is to statistically reconstruct the density matrix $\rho$ from the counts dataset $\mathcal{D}$. This can be seen as a statistical inverse problem of reversing the measurement map

[TABLE]

in the sense that given a sample $\mathcal{D}$ from $p_{\rho}(\cdot|\mathcal{S})$, one would like to produce an estimator $\hat{\rho}(\mathcal{D})$ which is close to $\rho$ with respect to some meaningful distance measure. The next section lists several figures of merit considered here.

2.1 Error functions

Let us denote the risk or mean error of an estimator $\hat{\rho}$ as $\mathbb{E}\left[D(\hat{\rho},\rho)\right]$, where $D(\hat{\rho},\rho)$ represents a particular error function. In our theoretical analysis and simulations study we estimate the risk for several choices of the error function $D(\hat{\rho},\rho)$, which are tabulated in Table 1. Note that the Bures distance is defined only over the space of density matrices, and therefore applies only in the case where the estimators are density matrices. The classical analogue of the Bures distance is the Hellinger distance between two probability distributions. Here we consider the Hellinger distance $D_{H}(\hat{\boldsymbol{\lambda}},\boldsymbol{\lambda})^{2}$ between the eigenvalues $\{\lambda_{1},\dots,\lambda_{d}\}$ of the true state and those of the estimator $\{\hat{\lambda}_{1},\dots,\hat{\lambda}_{d}\}$, seen as probability distributions. Later on we will show that the behaviour of the Bures distance error is strongly correlated with that of the squared Hellinger distance error.

3 Fisher information, asymptotic normality and maximum likelihood

The maximum likelihood (ML) estimator is one of the most commonly used and well understood statistical tools with a universal range of applicability. Before considering its use in quantum tomography, we briefly review some of the key concepts and results related to the ML estimator and its asymptotic behaviour [51, 78]. Consider the scenario in which we are given $N$ independent samples $X_{1},\dots,X_{N}$ from a discrete probability distribution $p_{\boldsymbol{\theta}}$ over a countable set $\Omega$; the probability distribution is assumed to depend smoothly on an unknown parameter $\boldsymbol{\theta}$ which belongs to an open subset $\Theta$ of $\mathbb{R}^{p}$. The likelihood of the dataset is $p_{\boldsymbol{\theta}}(X_{1},\dots,X_{N})=\prod_{j=1}^{N}p_{\boldsymbol{\theta}}(X_{j})$. The ML estimator is the point in $\Theta$ where the likelihood of the observed data attains its maximum value

[TABLE]

The likelihood can be expressed in terms of the empirical distribution of the data $\boldsymbol{f}=\sum_{j}\delta_{X_{j}}/N$ which collects the frequencies of different outcomes in $\Omega$. Indeed the log-likelihood can be written as (cf. (1))

[TABLE]

where $K(\cdot\|\cdot)$ relative entropy (Kullback-Leibler divergence), and $C,C^{\prime}$ do not depend on $\boldsymbol{\theta}$. The ML estimator is thus the closest point to $\boldsymbol{f}$ with respect to the relative entropy

[TABLE]

Asymptotics. The appeal of the ML estimator lies in part in its asymptotic optimality properties. For large enough sample sizes $N$, a central limit behaviour emerges and the ML estimator becomes normally distributed around the true parameter, with a variance shrinking at the rate $1/N$

[TABLE]

The convergence above holds in distribution as $N\to\infty$ and the limit is a centred Gaussian distribution with covariance given by the inverse of the Fisher information matrix. The latter is a $p\times p$ positive matrix

[TABLE]

which encapsulates how informative a single sample $X$ from $p_{\boldsymbol{\theta}}$ is about $\boldsymbol{\theta}$. The limiting covariance $I^{-1}(\boldsymbol{\theta})$ is equal to the lower bound appearing in the the Cramér-Rao inequality for unbiased estimators. Estimators which satisfy equation (2) are called efficient and statistical theory shows that one cannot improve on their asymptotic accuracy, except on a measure zero set of the parameter space [51]. In particular if $d(\cdot,\cdot)$ is a locally quadratic (positive) loss function on $\Theta$,

[TABLE]

then the risk of the MLE satisfies

[TABLE]

It is important to note that the asymptotic normality property (2) relies on the smoothness of the model and the fact that $\boldsymbol{\theta}$ is an interior point of the parameter space $\Theta$. For large $N$ the ML estimator lies (with high probability) within in a small neighbourhood of $\boldsymbol{\theta}$ of size $1/\sqrt{N}$, and the parameter space looks like $\mathbb{R}^{p}$ for all practical purposes. However, when $\Theta$ has a boundary, the ML estimator will not be asymptotically normal for parameters $\boldsymbol{\theta}$ lying on the boundary, and one needs to analyse such models more carefully. This is the case in quantum tomography, when the unknown state is not full-rank and therefore lies on the boundary of $\mathcal{S}_{d}$, or is so close to the boundary that the asymptotic theory will kick in only for sample sizes that are much larger than those obtained in real experiments.

3.1 The maximum likelihood estimator in quantum tomography

We will now discuss in more detail the properties of the MLE in the specific case of quantum tomography. Given the measurement data encoded in the dataset $\mathcal{D}$, the MLE is defined as the maximum of $p_{\tau}(\mathcal{D}|\mathscr{S})$ over $\tau\in\mathcal{S}_{d}$. By passing to log-likelihood and discarding the constant factorial terms in (1), we obtain the following form of the estimator

[TABLE]

where $\boldsymbol{f}(\boldsymbol{o}|\boldsymbol{s})=N(\boldsymbol{o}|\boldsymbol{s})/m$ are the empirical frequencies of the data. The MLE is commonly used in quantum tomography [10, 14, 44, 32], and several implementation methods have been proposed including Hradil’s iterative algorithm [64]. Our specific implementation uses the CVX package for disciplined convex programming on MATLAB [2]. The general asymptotic theory guarantees that the ML estimator is asymptotically normal for full rank states, i.e. the interior of $\mathcal{S}_{d}$. To get more insight into the general asymptotic behaviour for a given state $\rho$, we will choose a local parametrisation defined in terms of the matrix elements of $\rho$ with respect to its eigenbasis. Let $\lambda_{1}\geq\lambda_{2}\dots,\geq\lambda_{d}>0$ be the eigenvalues of a full-rank state $\rho$. Then any neighbouring state can be written as

[TABLE]

where $\delta\rho_{\boldsymbol{\theta}}$ is trace zero matrix which is completely parametrised by

[TABLE]

The ML estimator $\hat{\rho}_{\rm ML}=\rho_{\hat{\boldsymbol{\theta}}_{\rm ML}}$ has parameter $\hat{\boldsymbol{\theta}}_{\rm ML}$ which is normally distributed around [math] with covariance $I^{-1}(\rho|\mathscr{S})/N$ where the Fisher information is the average of the individual informations for different sets $\mathbf{s}\in\mathscr{S}$

[TABLE]

In particular, for any locally quadratic loss function $d(\cdot,\cdot)$ (e.g. Frobenius distance, or Bures distance) with weight matrix $G(\boldsymbol{\theta})$, the associated risk has the asymptotic behaviour

[TABLE]

Now, let us assume that the unknown state $\rho$ is on the boundary of $\mathcal{S}_{d}$, and more precisely belongs to the space of rank $r$ states $\mathcal{S}_{d}(r)\subset M(\mathbb{C}^{d})$, for a fixed and known rank $r\leq d$. In its own eigenbasis, $\rho$ is the diagonal matrix of eigenvalues ${\rm Diag}(\lambda_{1},\dots,\lambda_{r},0,\dots,0)$, and any sufficiently close state is uniquely determined by its matrix elements in the first $r$ rows (or columns). Intuitively this can be understood by noting that any rank-$r$ state $\rho^{\prime}$ in the neighbourhood of $\rho$ can be obtained by perturbing the eigenvalues and performing a small rotation of the eigenbasis; in the first order of approximations these transformation leave the $(d-r)\times(d-r)$ lower-right corner unchanged so

[TABLE]

where $C=O(\|A\|^{2},\|B\|^{2})$. We therefore choose the (local) parametrisation $\rho^{\prime}=\rho_{\boldsymbol{\theta}}$ where only the matrix elements of $A$ and $B$ enter the parameter $\boldsymbol{\theta}$. For this model, any rank-$r$ state corresponds to an interior point of the parameter space $\mathcal{S}_{d}(r)$, and consequently the ML estimator obeys the asymptotic normality theory described above.

However, if the rank of $\rho$ is not known in advance, then the diagonal block $C$ of $\rho^{\prime}$ needs to be included in the model in order to describe neighbouring states of higher rank. As $\rho^{\prime}$ is a state, the block $C$ is a positive matrix, and therefore its matrix elements are constrained. This complicates the analysis of the likelihood function, and invalidates the asymptotic normality of the ML estimator.

What should be the theory replacing asymptotic normality ? At the moment there isn’t a complete answer to this question, but some important progress has been made in [66]. Following this work and [17], it is instructive to study an extended, ‘non-physical’ model in which the positivity requirement is dropped and (locally around $\rho$) the parameter space is taken to be that of selfadjoint matrices of trace-one $M_{sa}^{1}(\mathbb{C}^{d})$. Therefore, the ‘unphysical’ parameter $\boldsymbol{\theta}$ consists of the matrix elements of the blocks $A,C,B$ (except one diagonal element due to normalisation), with $\rho^{\prime}=\rho_{\boldsymbol{\theta}}$ given by equation (7). We now make the ‘mild’ assumption that $p_{\rho}(\boldsymbol{o}|\boldsymbol{s})>0$ for all pairs $(\boldsymbol{o},\boldsymbol{s})$; indeed this condition is satisfied for ‘generic’ states but fails for states whose support is orthogonal to some of the measurement basis vectors. Under this assumption, we find that locally around $\rho$, we can define a statistical model given by bona-fide probability distributions $p_{\boldsymbol{\theta}}:=p_{\rho_{\boldsymbol{\theta}}}.$ We can therefore define the unconstrained maximum likelihood (uML) estimator as in equation (3) where the optimisation is performed over the unconstrained local neighbourhood $\mathcal{O}\subset M^{1}_{sa}(\mathbb{C}^{d})$ of $\rho$, rather than over the space of states $\mathcal{S}_{d}$. Since $\rho$ is on the boundary of $\mathcal{S}_{d}$, the probability that $\hat{\rho}_{\rm uML}$ falls outside the state space is significant and does not vanish in the asymptotic limit. In this case, $\hat{\rho}_{\rm uML}$ does not coincide with the ‘constrained’ or regular ML estimator $\hat{\rho}_{\rm ML}$. In fact, each of them can be seen as the projection with respect to the relative entropy distance of the empirical distribution $\boldsymbol{f}$ onto the sets of probabilities $\mathcal{M}(\mathcal{O})$ and respectively $\mathcal{M}(\mathcal{S}_{d})$, where $\mathcal{M}:M(\mathbb{C}^{d})\to\mathbb{R}^{k\cdot d}$ is the measurement map.

Asymptotics. Since $\rho$ corresponds to an interior point of the unconstrained parameter space, the general asymptotic normality theory applies again. In particular, the uML estimator is normally distributed around $\boldsymbol{\theta}=0$ with variance $(NI(\rho|\mathscr{S}))^{-1}$, where $I(\rho|\mathscr{S})$ is the Fisher information (5) computed with respect to the unconstrained parametrisation. Moreover by Taylor expanding the log-likelihood function around the $\boldsymbol{\theta}_{\rm uML}$, and using the fact the the first derivative vanishes at this point, we obtain the second order approximation

[TABLE]

This implies tht for large $N$, the ML estimator $\hat{\rho}_{\rm ML}$ is the projection of $\hat{\rho}_{\rm uML}$ onto $\mathcal{S}_{d}$ with respect to the quadratic distance determined by the Fisher information

[TABLE]

Moreover, the probability distribution of the MLE is obtained by projecting the Gaussian distribution corresponding to $\hat{\rho}_{\rm uML}$, onto $\mathcal{S}_{d}$. Although the projection can be computed efficiently using convex optimisation, finding a general characterisation of the resulting distribution seems to be a challenging problem [71]. Nevertheless, [66] shows that the problem is tractable in those cases where the metric $d_{I}$ is (approximately) isotropic, so that random matrix results such as the emergence of Wigner semicircle law can be used to treat the asymptotic theory. In section 8.1 we will use these ideas to study the asymptotic behaviour of the projected least squares estimator.

4 The least squares estimator

The least squares (LS) estimator [63, 17] is based on the observation that the outcome probabilities $p_{\rho}(\boldsymbol{o}|\boldsymbol{s})$ depend linearly on the state $\rho$. Let

[TABLE]

be the decomposition with respect to an orthonormal basis of $M(\mathbb{C}^{d})$ consisting of selfadjoint matrices $\{\tau_{i}:1\leq i\leq d^{2}\}$, and $\boldsymbol{\beta}:=(\beta_{1},\dots,\beta_{d^{2}})^{T}$ the corresponding vectorisation. Then the probabilities can be expressed as

[TABLE]

where $X$ is a $kd\times d^{2}$ fixed matrix depending on the measurement and the chosen vectorisation. In an experimental setup we do not have access to the true probability vector. Instead from the $d\times k$ dataset of counts, we can compute the empirical probabilities $f(\boldsymbol{o}|\boldsymbol{s}):=N(\boldsymbol{o}|\boldsymbol{s})/m$, whose expectations are $\mathbb{E}f(\boldsymbol{o}|\boldsymbol{s})=p_{\rho}(\boldsymbol{o}|\boldsymbol{s})$. Replacing probabilities vector with empirical frequencies

[TABLE]

where $\boldsymbol{\epsilon}\in\mathbb{R}^{dk}$ is a mean zero vector of statistical noise. The noise distribution is irrelevant for the definition of the LS estimator, but we will return to this when discussing its error, and the generalised least squares estimator.

The least squares solution to the above system of equations is defined as the minimiser of the following optimisation problem

[TABLE]

and has the explicit form $\hat{\boldsymbol{\beta}}_{\rm LS}=(X^{T}X)^{-1}\cdot X^{T}\cdot\boldsymbol{f}$. The final estimate $\hat{\rho}_{\rm LS}$ of the density matrix is then constructed from the estimated parameter vector $\hat{\boldsymbol{\beta}}_{\rm LS}$ using equation (10). We note that the least squares estimator produces a state estimate $\hat{\rho}_{\rm LS}=\rho_{\boldsymbol{\beta}_{\rm LS}}$ that is not necessarily a density matrix.

Least squares as a projection. Let $\mathcal{M}:M(\mathbb{C}^{d})\to\mathbb{C}^{k\cdot d}$ be the map associated to the measurement, and let $\hat{p}_{\rm LS}:=\mathcal{M}(\hat{\rho}_{\rm LS})=X\hat{\boldsymbol{\beta}}_{\rm LS}$ be the ‘probability distribution’ corresponding to the LS estimator. This belongs to the hyperplane $\mathcal{L}_{d}=\mathcal{M}(M_{sa}^{1}(\mathbb{C}^{d}))$ which contains the convex set of probabilities $\mathcal{P}_{d}=\mathcal{M}(\mathcal{S}_{d})$. The LS estimator (12) can then be interpreted as the projection of the frequency vector $\boldsymbol{f}$ onto $\mathcal{L}_{d}$ with respect to the Euclidian distance in $\mathbb{R}^{k\cdot d}$.

4.1 Least squares as inversion of a measure-and-prepare channel

The LS estimator was introduced above by choosing a particular vectorisation of the density matrix. While this is useful for numerical implementations, the disadvantage is that one loses the physical interpretation of the quantum state and the measurement map. Conceptually, it is useful to define the LS estimator in a “coordinate free” way which can be interpreted as the inversion of a certain measure-and-prepare map associated to the measurement [36]. Let $\mathcal{D}:M(\mathbb{C}^{d})\to M(\mathbb{C}^{d})$ be the quantum channel

[TABLE]

which is the composition $\mathcal{P}\circ\mathcal{M}$ of the measurement with collection of bases $\mathscr{S}$, and the preparation map where the pure state $P_{\boldsymbol{o}}^{\boldsymbol{s}}$ is prepared for each outcome $\boldsymbol{o}$ of the measurement in basis $\boldsymbol{s}$. If $\rho$ is represented in its vectorised form, then the map $\mathcal{D}$ is given by the matrix $X^{T}X/|\mathscr{S}|$. On the other hand, the preparation map $\mathcal{P}$ has matrix $X^{T}$, so that the LS estimator can be expressed as

[TABLE]

From this expression we see immediately that ${\rm Tr}(\hat{\rho}_{\rm LS})=1$ as consequence of the fact that $\boldsymbol{f}/|\mathscr{S}|$ is a probability distribution and $\mathcal{D}$ is trace preserving. Additionally, we note the the accuracy of the LS estimator is linked to the properties of the channel $\mathcal{D}$. In particular, in the case of Pauli measurements the channel is given by a tensor product $\mathcal{D}=\mathcal{C}^{\otimes n}$ of qubit depolarising channels [17, 36]

[TABLE]

On the other hand, the measurements corresponding to the class of 2-designs (which includes covariant measurements, mutually unbiased bases, stabiliser states, symmetric informationally complete measurements) have associated channel given by [36]

[TABLE]

4.2 Concentration bounds and asymptotic behaviour of the LS

For more insight into the structure of the LS estimator let us unpack equation (13) and note that $\hat{\rho}_{\rm LS}$ can be written as an average of independent, identically distributed matrices $A^{(1)},\dots,A^{(m)}$ with $\mathbb{E}(A^{(i)})=\rho$

[TABLE]

where ${\bf o}_{i|{\bf s}}$ are the outcomes of the measurements with respect to setting ${\bf s}$.

In this form, the LS estimator is amenable to non-asymptotic matrix concentration techniques, as well as asymptotic normality theory. The following result [36] was obtained by applying matrix Bernstein inequalities [76] and provides a non-asymptotic confidence bound for LS with respect to the operator norm distance, see also [17].

Theorem 4.1**:**

Let $\hat{\rho}_{\rm LS}$ be the LS estimator of $\rho$ for a dataset consisting of $N=m\times k$ samples. Then

[TABLE]

where $g(d)=2d$ for $2-$design measurements and $g(d)\simeq d^{1.6}$ for Pauli measurements.

The theorem provides upper bounds for risks with respect to commonly used loss functions such as the Frobenius distance (norm-two squared) and the trace-norm distance, see also [72] for related results. Indeed using the basic inequalities $\|A\|_{2}^{2}\leq d\|A\|^{2}$ and $\|A\|_{1}\leq d\|A\|$ we obtain the upper bounds

[TABLE]

for the 2-design measurements and

[TABLE]

for the Pauli basis measurement, where $c_{1},c_{2},C_{1},C_{2}$ are a numerical constants. Moreover, in [36] it was shown that the log factor in (15) can be removed when we deal with covariant measurements. On the other hand, in section 3 we have shown that for the maximally mixed state, the optimal convergence rate for the Frobenius distance is of the order $O(d^{2}/N)$; this means that the upper bound (15) cannot be improved, when seen as a uniform bound over all states. However, this leaves open the possibility that special classes of states can be estimated with higher accuracy that that provided by the LS estimator. Indeed, we will see that simple modifications of the LS estimator which take into account the positivity of the unknown state can significantly reduce the estimation errors for low rank states.

Asymptotics. Thanks to its simplicity, the LS estimator has a tractable asymptotic theory. As in the case of the concentration Theorem 4.1 the key point is that the error $\rho_{\rm LS}$ is a sum of independent, identically distributed matrices given by equation (14). In the limit of large $m$ each matrix element of $\rho_{\rm LS}$ has a Gaussian distribution; in terms of the vectorised form we have the Central Limit for $m\to\infty$

[TABLE]

where $\Omega:=m{\rm Cov}\left[\boldsymbol{\epsilon}|\mathscr{S}\right]$ is the (renormalised) covariance of the noise $\boldsymbol{\epsilon}$. Due to the independent structure of the measurement settings, the latter has a non-trivial block diagonal form with $d\times d$ blocks corresponding to the different settings $\boldsymbol{s}$

[TABLE]

This allows us in principle to compute the asymptotic risk of an arbitrary quadratic loss function with weight matrix $G$ (such as the Frobenius distance) as

[TABLE]

Since $\boldsymbol{f}$ is a vector of frequencies, it must satisfy $k$ normalisation constraints, which is reflected in the fact each block $\Omega_{\bf s}$ is singular, with zero eigenvector $\mathbf{1}=(1,\dots,1)^{T}$. Consequently, the covariance matrix $V_{\rm LS}$ is singular with the zero eigenvector corresponding to the trace which is a fixed parameter. We will come back to this when discussing the generalised least squares estimator.

4.3 Asymptotic theory of LS for covariant measurements

While asymptotic normality tells us that the estimator $\hat{\boldsymbol{\beta}}_{\rm LS}$ lies in an ellipsoid centred at $\boldsymbol{\beta}$, it treats the estimator as a vector rather that as a matrix. For instance, it does not immediately provide an asymptotic theory for the operator norm error $\|\hat{\rho}_{\rm LS}-\rho\|$. Random matrix theory provides us with other asymptotic results which take into account the matrix structure of the estimator. To explore this avenue we will make a simplifying assumption and place ourselves in the scenario of covariant measurements. The outcome of such a measurement is a one dimensional projection $P=|\psi\rangle\langle\psi|$ and the corresponding positive operator valued measure is

[TABLE]

where ${\rm d}P$ is the uniform measure over the space of one dimensional projections; the latter is the measure induced by the Haar measure over the unitaries $U$ via the action $P\mapsto UPU^{*}$. An alternative way of obtaining a measurement sample is to choose a random basis and perform a measurement with respect to the chosen basis.

In this setup, the channel $\mathcal{D}$ is given by [36]

[TABLE]

and the LS estimator given by equation (14) can be written as

[TABLE]

where $P_{i}$ are the independent outcomes of measurements.

We will be particularly interested in the behaviour of the LS estimator for low rank states, as well as states close to the maximally mixed one. Due to covariance it suffices to choose states which are diagonal in the standard basis $\{|i\rangle:i=1,\dots,d\}$, and for simplicity we will restrict ourselves to rank-r states with equal eigenvalues $\rho_{r}=\sum_{i=1}^{r}|i\rangle\langle i|/r$. As in section 3.1 we write the LS estimator as

[TABLE]

where $A,B,C$ are block matrices of mean zero, and $I_{r}$ is the $r\times r$ identity block. We will deal with each block separately.

Block C. By covariance, the distribution of the block $C$ is invariant under unitary transformations in $\mathbb{C}^{d-r}$. As the LS estimator is unbiased, the matrix elements of $C$ are centred. The off-diagonal elements have real and imaginary parts with variances equal to

[TABLE]

with

[TABLE]

where we have written $P=U|k\rangle\langle k|U^{*}$ for $r<k\neq i,j$ and replaced the integration over projections with that over unitaries; the integral was then evaluated using Weingarten formulas [19]. Similarly, the variance and covariance of the diagonal elements are

[TABLE]

and all off-diagonal elements (of all blocks $A,B,C$) have zero covariance with other matrix elements. By the Central Limit Theorem, in the limit of large $N$, the off-diagonal elements of $\hat{\rho}_{\rm LS}$ become normally distributed and independent of all other elements. On the other hand, as the covariance of diagonal elements is non-zero, they converge to correlated Gaussian variables. However, if we also take the limit of large $d$ we find that the covariance of the diagonal elements vanishes and the matrix $C$ is distributed approximately as the Gaussian unitary ensemble (GUE). Universality results for random matrices [65] imply that the empirical distribution of the eigenvalues of $C/\sqrt{d-r}$ converges weakly to Wigner’s semicircle law on the interval $[-2,2]$, whose probability density is

[TABLE]

Indeed, panel a) of Figure 2 shows a good match between the histogram of the eigenvalues of the error block $C/\sqrt{N}$ for a rank-one state of $n=7$ atoms and $N=10^{6}$ samples, and the corresponding Wigner distribution.

Block A. Let us consider now block $A$ of equation (19). For the same symmetry reasons, its off-diagonal elements (when $r>1$) have real and imaginary parts with variance

[TABLE]

with

[TABLE]

where all indices $i,j,k,l$ are different. The diagonal elements have variance

[TABLE]

while the covariance of the diagonal elements is

[TABLE]

We note that for small $r$ the diagonal elements do not become independent when $N$ and $d$ are large. However, if we also take $r$ to be large, the covariance vanishes, as in the case of the $C$ block. By the same argument, the distribution of the eigenvalues of $A/\sqrt{r}$ converges to the Wigner semicircle law (20). This is illustrated in the right panel of Figure 2, for the case of the totally mixed state of $n=7$ atoms.

Block B. The elements of block B have variances equal to

[TABLE]

with $i\leq r<j$ and

[TABLE]

Eigenvalues distribution of the LS estimator. To better understand the LS estimator, it is instructive to study the distribution of its eigenvalues. This will also be used in the study of the projected least squares estimator furher on. Assuming $N$ is large, we note that $\hat{\rho}_{\rm LS}$ can be written as

[TABLE]

where

[TABLE]

This means that in the leading order in $N^{-1/2}$ the set of eigenvalues of $\hat{\rho}_{\rm LS}$ is the union of those of $I_{r}/r+A/\sqrt{N}$ and $C/\sqrt{N}$. In the case of large $N,d$ and small rank $r$, the first group of eigenvalues are small fluctuations of size $r/\sqrt{N}$ around $1/r$; in particular these eigenvalues are positive with high probability; the second group of eigenvalues are distributed according to the Wigner law, and have maximum absolute value approximately equal to $2\sqrt{(d-r)/N}$, and roughly half of them will be negative.

We can now evaluate the asymptotic risk of the LS estimator with respect to the chosen loss functions.

Frobenius risk. For the Frobenius distance, the asymptotic mean square error is

[TABLE]

In particular, the leading contribution to the Frobenius risk is $d^{2}/N$, and the dependence with $r$ is weak. As illustrated in panels a) and b) of Figure 3 the theoretical prediction match the simulation results for 5 and 6 atom states.

Operator norm risk. The operator-norm error of the LS estimator is

[TABLE]

Above we found that the entries of $A,B,C$ become independent in the limit of large $N$, except for correlations between the diagonal elements which vanish if we additionally take the limit of large $d$. Moreover, the variances of the elements in the $B$ block and the off-diagonal elements of $A$ differ from those of the off-diagonal elements of $C$ by factors $(r+2)/r$ and respectively $(r+1)/r$. Therefore, the error block matrix as a whole is not distributed according to GUE ensemble. However, the universality of the Wigner semicircle law, the limit holds not only for highly symmetric ensembles like GUE but also for ‘small’ perturbations of this ensemble, e.g. random matrices with independent entries, whose variances do not deviate too much from a fixed value [23]. In particular for low rank $r\ll d$, the total size of the blocks $A,B,B^{*}$ is of the order $rd$ which is much smaller than than the size $(d-r)^{2}$ of $C$. Therefore, the asymptotic behaviour of the error matrix is determined by that of $C$ and its spectrum converges to the Wigner law in the limit of large sample size and large dimension. In particular, the leading contribution to the norm error $\|\hat{\rho}_{\rm LS}-\rho_{r}\|$ for low rank states is $2\sqrt{d/N}$. A similar situation occurs for high rank states $r\approx d$ where the dominant variances are provided by the block $A$.

Indeed simulations for a pure and the fully mixed state of $n=7$ atoms with $N=10^{6}$ samples gave norm errors $0.0225$ and respectively $0.0221$, while the above theoretical prediction is $0.0226$. Further simulation results for $n=5$ and $n=6$ atoms states of different ranks show a good match with the above estimate.

Trace norm risk. As in the case of the operator norm, the trace norm error can be estimated thanks to the fact that the eigenvalues of error matrix follow an approximate Wigner semicircle distribution. Using

[TABLE]

we find that for low rank, or close to full rank states, the leading contribution to $\mathbb{E}(\|\hat{\rho}_{\rm LS}-\rho\|_{1})$ is $8d^{3/2}/(3\pi\sqrt{N})$. This rate agrees with the upper bound in equation (15), but provides an exact asymptotic constant for this rate. Indeed simulations with a pure and maximally mixed states of $n=7$ atoms and $N=10^{6}$ samples gave norm-one errors of $1.228$ and respectively $1.229$ while the theoretical prediction is $1.229$.

5 Generalised least squares estimator

Consider the generic linear regression problem

[TABLE]

with ${\bf x}\in\mathbb{R}^{a}$ and unknown vector, ${\bf y}\in\mathbb{R}^{b}$ a vector of observations and ${\bf n}$ the ‘noise’ term with fixed and known covariance matrix $C$, which is assumed to be strictly positive. While the LS estimator $\hat{\bf x}_{\bf LS}$ is optimal in the sense of minimising the prediction error $\|A\hat{\bf x}-{\bf y}\|$, this is not the case when considering an estimation error e.g. the mean square error $\mathbb{E}\|\hat{\bf x}-{\bf x}\|^{2}$, unless the covariance matrix is proportional to the identity. To recover the ‘equal noise’ situation we multiply equation (22) from the left by $C^{-1/2}$, to obtain

[TABLE]

The LS estimator of the last regression equation has the smallest covariance matrix among linear unbiased estimators

[TABLE]

This is called the generalised least squares (GLS) estimator. when the noise distribution is Gaussian, the GLS estimator coincides with the MLE. In the i.i.d. setting where $m$ samples are available, the GLS estimator is asymptotically normal and efficient

[TABLE]

GLS as a projection. Similarly to the LS estimator we consider the image $\hat{{\bf y}}_{\rm GLS}=A\hat{\bf x}_{\rm GLS}$ of the GLS estimator ${\bf x}_{\rm GLS}$. Then $\hat{{\bf y}}_{\rm GLS}$ is the projection of the data ${\bf y}$ onto the subspace ${\rm Ran}(A)\subset\mathbb{R}^{b}$ with respect to the covariance-dependent metric

[TABLE]

GLS for tomography. Let us return now to the tomography regression problem of the form given by equation (11). While the LS estimator is optimal in the sense of equation (12), in general it is not optimal in the estimation sense for any locally quadratic distance function. To remedy this we would like to construct a corresponding generalised least square estimator to take into account the nontrivial form of the noise covariance matrix $\Omega$ given by equation (17). However, there are two issues which prevent us from directly applying the GLS methodology to the tomography data. The first is that $\Omega$ is unknown since it depends on the true probabilities, and therefore on the unknown state $\rho$. We therefore propose to use an estimate of the covariance matrix instead. We describe the computation of this estimate in Appendix A.2. The second difficulty is that $\Omega$ is a singular matrix due the constraint $\sum_{\boldsymbol{o}}f(\boldsymbol{o}|\boldsymbol{s})=1$ for each setting ${\bf s}$. This means that for each setting ${\bf s}$ the vector $|{\bf 1}_{\bf s}\rangle\in\mathbb{C}^{d}$ is a zero eigenvector for $\Omega_{\bf s}$ and one should work within the orthogonal complement of $|{\bf 1}_{\bf s}\rangle$. This can be achieved by choosing a block-diagonal isometry $V:\mathbb{C}^{k(d-1)}\to\mathbb{C}^{k\cdot d}$ such that each block $V_{\bf s}$ satisfies $\langle{\bf 1}_{s}|V_{\bf s}=0$, and defining the new ‘frequencies’ vector $\tilde{\boldsymbol{f}}$ with settings components $\tilde{\boldsymbol{f}}_{\bf s}=V_{\bf s}^{*}\boldsymbol{f}_{\bf s}$. Similarly, we can also remove the trace-one constraint for states by choosing the vectorisation (10) such that the last basis element is $\tau_{d^{2}}=\mathbb{1}/\sqrt{d}$. In this case the state is uniquely determined by the free parameters $\{\beta_{i}:1\leq i\leq d^{2}-1\}$. If we denote by $J:\mathbb{C}^{d^{2}-1}\to\mathbb{C}^{d^{2}}$ the isometric embedding with respect to the standard basis, then $\tilde{\boldsymbol{\beta}}:=J^{*}\boldsymbol{\beta}$ is the truncated vector of free parameters for $\rho$.

The regression problem can now be written in terms of the ‘tilde’ vectors and matrices

[TABLE]

where $\tilde{X}=V^{*}XJ$ and $\tilde{\boldsymbol{\epsilon}}$ has covariance matrix $\tilde{\Omega}=V^{*}\Omega V$ which is non-singular. We can now apply the GLS methodology and define the estimator as

[TABLE]

where $\hat{\tilde{\Omega}}$ is the (non-singular) estimated covariance matrix. We denote $\hat{\rho}_{\rm GLS}$ as estimate of the density matrix constructed from its vectorised form $\hat{\boldsymbol{\beta}}_{\rm GLS}$. For future use, let us denote $\tilde{\mathcal{M}}:M^{1}_{sa}(\mathbb{C}^{d})\to\mathbb{R}^{k(d-1)}$ the map $\rho\mapsto\tilde{X}\tilde{\boldsymbol{\beta}}$, so that $\hat{\tilde{p}}_{\rm GLS}=\tilde{\mathcal{M}}(\hat{\rho}_{\rm GLS})$ .

5.1 Asymptotic theory of GLS

What are the asymptotic properties of $\hat{\rho}_{\rm GLS}$? For large $m$, the noise distribution becomes Gaussian, and the covariance estimator converges to the actual covariance. The estimator $\hat{\boldsymbol{\beta}}_{\rm GLS}$ becomes asymptotically normal as $m\to\infty$

[TABLE]

which means that the corresponding asymptotic average Fisher information per sample is $I_{\rm GLS}=\tilde{X}^{*}\tilde{\Omega}^{-1}\tilde{X}/k$. In appendix A.1 we show that $I_{\rm GLS}$ coincides with the Fisher information $I(\rho|\mathscr{S})$ defined as in equation (5) for the parametrisation $\tilde{\boldsymbol{\beta}}$. This equality implies that the error rates of the GLS estimator have the same asymptotic behaviour as those of the uML estimator, and both estimators satisfy asymptotic normality. In fact one can make the stronger statement that the two estimators are asymptotically close to each other.

Equivalence of GLS and uML. As stated above, the GLS can be interpreted as the projection of the data $\tilde{\boldsymbol{f}}$ onto the image of $\tilde{X}$ in $\mathbb{R}^{k(d-1)}$ with respect to the metric

[TABLE]

On the other hand, in section 3.1 we showed that for large $N$, we can define the uML estimator as the projection of $\boldsymbol{f}$ onto the hyperplane $\mathcal{L}_{d}:=\mathcal{M}(M_{sa}^{1}(\mathbb{C}^{d}))$, with respect to the relative entropy distance. Since in the first order of approximation the relative entropy is given by the quadratic form $\tilde{\Omega}^{-1}$, the two projections become identical in the asymptotic limit.

6 Thresholded least squares estimator

The LS estimator (as well as GLS) suffers from the disadvantage that it does not necessarily produce a density matrix, i.e, a positive semi-definite estimate of trace one. While the significant eigenvalues can be estimated reasonably well with enough data, the LS estimator will typically have negative eigenvalues corresponding to small or zero eigenvalues of the true state. The thresholded least squares estimator (TLS) proposed in [17], improves the LS estimator by selecting the state which is the closest in Frobenius norm to $\hat{\rho}_{\rm LS}$ and whose non-zero eigenvalues are above a certain threshold $\nu\geq 0$, i.e.

[TABLE]

with $\mathcal{S}_{d}^{(\nu)}$ the set of density matrices with spectrum in $\{0\}\cup\{[\nu,1]\}$. The choice of the statistical noise threshold is informed by the accuracy of the LS estimate, and a theoretical and a ‘data-driven’ choices for this threshold are detailed in [17], see also [69, 22]. In practice it is found that estimator’s performance improves if the threshold is allowed to be ‘data-driven’, by using cross-validation to choose the optimal value of the threshold, see Appendix A.2.

It turns out that computing $\hat{\rho}_{\rm TLS}$ is very efficient and can be easily implemented. Let

[TABLE]

be the spectral decomposition of $\hat{\rho}_{\rm LS}$ where we assume that the eigenvalues are sorted in descending order $\hat{\lambda}_{1}\geq\ldots\geq\hat{\lambda}_{d}$. The thresholded estimator

[TABLE]

has the same eigenvectors as $\hat{\rho}_{\rm TLS}$, and its eigenvalues can be computed in terms of $\hat{\lambda}_{i}$ as summarised in Algorithm 1.

In words, the algorithm checks if the smallest eigenvalue of $\hat{\rho}_{\rm LS}$ is above the noise threshold, and if it is then $\hat{\rho}_{\rm TLS}=\hat{\rho}_{\rm LS}$ ( note that by construction ${\rm Tr}(\hat{\rho}_{\rm LS})=1$, so there is no need to renormalise the LS estimator). On the other hand if the smallest eigenvalue is below the threshold, it is set to zero and the remaining eigenvalues are suitably shifted accordingly so that their sum is equal to one. The final estimate $\hat{\rho}_{\rm TLS}$ is constructed by replacing the eigenvalues of $\hat{\rho}_{\rm LS}$ with these thresholded eigenvalues $\hat{\lambda}_{i}(\nu)$. The theoretical properties of the TLS estimator are presented in [17] and are similar to those of the projected least squares estimator which is discussed in more detail below.

6.1 Thresholded Generalised Least Squares Estimator (TGLS)

This estimator is obtained by using the GLS estimate $\hat{\rho}_{\rm GLS}$ instead of the LS estimate as a starting point for the thresholding procedure. The constant for thresholding is chosen in the same way by using cross-validation. The advantage of TGLS is that it starts from a more accurate estimator than LS, which leads to smaller errors compared to TLS. The price to pay is that it is somewhat more involved computationally, although still faster than ML.

7 Positive least squares estimator

The positive least squares (posLS) estimator is obtained by restricting the minimisation in (12) to parameters that correspond to density matrices $\tau\in\mathcal{S}_{d}$. Let $\mathcal{M}:M(\mathbb{C}^{d})\mapsto\mathbb{C}^{kd}$ be the measurement map defined as $[\mathcal{M}(\tau)]_{\boldsymbol{o},\boldsymbol{s}}={\rm Tr}\left[\tau P^{\boldsymbol{s}}_{\boldsymbol{o}}\right]/|\mathcal{S}|$, Then, similarly to section 4.1 we can express the posLS estimator as

[TABLE]

By comparing with (12) we see that $\hat{\rho}_{\rm posLS}$ is the projection of $\hat{\rho}_{\rm LS}$ with respect to the distance on matrices induced by the euclidian distance on measurement probabilities $d_{\mathcal{M}}(\rho,\tau)=\|\mathcal{M}(\rho)-\mathcal{M}(\tau)\|$. To our knowledge, with the exception of the results in [46], its theoretical properties have not been studied in detail. While its statistical performance greatly improves on that of the LS, the posLS estimator has the drawback that it cannot be expressed in a closed form, and its computational complexity is comparable to that of the ML estimator.

7.1 Positive least squares estimator for covariant measurements

Let us consider the special case of the covariant measurement defined in equation (18). This measurement maps states into probability distributions in an (almost) isometric way [36]

[TABLE]

Indeed using Weingarten formulas [19] we get for any traceless $X$

[TABLE]

As a corollary, we find that the posLS estimator is the projection of the LS estimator with respect to the Frobenius distance. Therefore, for this specific measurement, the posLS estimator coincides with the projected least squares estimator which will be discussed in section 8.

7.2 Positive generalised least squares estimator

This estimator is defined in much the same way as the posLS estimator, by restricting the minimisation in (5) to parameters that correspond to density matrices. In keeping with the discussion in section 5, we consider the truncated frequency vector $\tilde{\boldsymbol{f}}$. By analogy with the GLS estimator, the positive generalised least squares (posGLS) estimator is defined as

[TABLE]

We will show that $\hat{\rho}_{\rm posGLS}$ is asymptotically equivalent to the ML estimator $\hat{\rho}_{\rm ML}$.

Asymptotic equivalence of posGLS and ML. As in the case of GLS and uML, it is easier to work with the images in $\mathbb{C}^{k(d-1)}$ of the different estimators through the (injective) map $\tilde{\mathcal{M}}$; we denote by $\hat{\tilde{p}}_{\rm posGLS}=\tilde{\mathcal{M}}(\hat{\rho}_{\rm posGLS})$ the probability vector corresponding to the posGLS estimator, and similarly for other estimators. Let us equip the the space of ‘frequencies’ $\mathbb{C}^{k(d-1)}$ with the distance

[TABLE]

Therefore, $\hat{\tilde{p}}_{\rm posGLS}$ is the projection with respect to $d_{\Omega}$ of $\tilde{\bf f}$ onto the convex set

[TABLE]

On the other hand, the GLS estimator is the projection of $\tilde{\bf f}$ onto the hyperplane $\tilde{\mathcal{M}}(M^{1}_{sa}(\mathbb{C}^{d}))$ which contains $\mathcal{P}_{d}$. Therefore, by properties of projections on convex subsets we find that $\hat{\tilde{p}}_{\rm posGLS}$ is the projection of $\hat{\tilde{p}}_{\rm GLS}$ onto $\mathcal{P}_{d}$.

In section 5.1 we showed that the GLS estimator is asymptotically equivalent to the uML, as a consequence of the fact asymptotically with $m$ both projections are determined by the Fisher information metric. Since posGLS and ML are obtained by applying the same projections onto the smaller space $\mathcal{P}_{d}$, they are also asymptotically equivalent. Unfortunately, this equivalence does not provide us with general estimation method which is more efficient than ML; the projection involved in posGLS does not seem to have closed form expression and requires a similar optimisation process as ML.

8 Projected least squares

Recently, it has been shown that the theoretical properties $\hat{\rho}_{\rm TLS}$ are preserved even if the threshold $\nu$ is chosen to be zero [36]. The projected least squares (PLS) estimator is defined as

[TABLE]

where the optimisation is performed over all states $\tau$. As in the case of the TLS estimator, this optimisation can be performed efficiently, and it only involves computing the spectral decomposition of $\hat{\rho}_{\rm LS}$ and applying Algorithm 1 with $\nu=0$. The PLS estimator is therefore faster than the data-driven TLS and turns out to have quite similar behaviour to the latter. Note that in general PLS is different from posLS , as both can be seen as projections of the LS estimator with respect to different metrics. However, as noted before, the estimators coincide in the case of covariant measurements. In the next section we will study this scenario in more detail.

Using the LS concentration bound of Theorem 4.1, the following rank dependent norm-one bound for the PLS estimator was derived in [36], where the measurement is either a 2-design or the Pauli bases measurement.

Theorem 8.1**:**

Let $\hat{\rho}_{\rm PLS}$ be the PLS estimator of $\rho$ for a dataset consisting of $N=m\times k$ samples. Then

[TABLE]

where $g(d)=2d$ for $2-$design measurements and $g(d)\simeq d^{1.6}$ for Pauli measurements.

This gives a upper bound of $O(r^{2}\cdot d\log d/\sqrt{N})$ on the convergence rate of norm-one for $2$-design measurements and $O(\log r^{2}\cdot d^{1.6}\log d/\sqrt{N})$ for Pauli measurements. In the first case the $\log d$ can be removed for covariant measurements and the resulting rate is optimal in the sense that it achives general lower bounds from [37].

8.1 The asymptotic behaviour of PLS for covariant measurements

In this section we look in more detail at the PLS estimator in the context of covariant measurements defined in equation (18). We have already seen that Theorem 8.1 provides a concentration bound on the norm-one risk of the PLS estimator. While such results are very valuable thanks to their non-asymptotic nature, it is instructive and useful to also understand the asymptotic behaviour for large sample size $N$ and dimension $d$. Indeed, in section 4.3 we showed how central limit and random matrix theory can be used to obtain tighter bounds on estimation risks of the LS estimator.

We will be particularly interested in the behaviour of PLS for low rank states. Due to covariance it suffices to choose states which are diagonal in the standard basis $\{|i\rangle:i=1,\dots,d\}$, and for simplicity we will restrict ourselves to rank-$r$ states with equal eigenvalues $\rho_{r}=\sum_{i=1}^{r}|i\rangle\langle i|/r$. As in section 4.3 we write the PLS estimator as

[TABLE]

for some error blocks $\tilde{A},\tilde{B},\tilde{C}$ whose dependence on the $A,B,C$ blocks of the LS estimator needs to be determined.

Our analysis draws on the asymptotic theory of the LS described in section 4.3, and the arguments developed in [66] for the analysis of the ML estimator. We know that $\hat{\rho}_{\rm PLS}$ has the same eigenbasis as $\hat{\rho}_{\rm LS}$, and its eigenvalues are obtained from those of $\hat{\rho}_{\rm LS}$ by the truncation procedure which involves repeatedly setting to zero negative eigenvalues and shifting the remaining ones so that they add up to one. The following 3 step procedure aims to identify the leading contributions to PLS for low rank states, and $N\gg d\gg 1$.

1. Diagonalisation. Recall that for large $N$ the LS estimator can be block diagonalised by means of a ‘small’ unitary rotation $U$, cf. equation (21)

[TABLE]

Therefore the eigenvalues of $\hat{\rho}_{\rm LS}$ can be grouped in two sets. The first is the set of eigenvalues of the block $C/\sqrt{N}$, which for large $d$ are distributed according to the Wigner semicircle law and lie between $\pm 2\sqrt{(d-r)/N}$. The second set consists of the eigenvalues of the block $I_{r}/r+A/\sqrt{N}$, which are small fluctuations of order $r/\sqrt{N}$ around $1/r$.

2. Truncation. As long as $N\gg d\gg 1$ and $r$ is small, the second set of eigenvalues is well separated from the first and it is very unlikely that any of these eigenvalues will be set to zero in the truncation process. Therefore, the cut-off point for the eigenvalues of $C/\sqrt{N}$ depends only on the sum of the larger eigenvalues

[TABLE]

Moreover, since $(d-r)\gg 1$, the eigenvalues of $C/\sqrt{N}$ are (approximately) distributed according to the Wigner distribution (20) on the interval $[-2\sqrt{(d-r)/N},2\sqrt{(d-r)/N}]$. Therefore, we can write the cut-off point $q$ as the solution of the following equation [66]

[TABLE]

where the integral corresponds to the sum of the eigenvalues of the block $C/\sqrt{N}$ above $q$, after being shifted by $q$. Writing $q=2\sqrt{d-r/N}\epsilon$ with $\epsilon=\epsilon(r,d,N,a)<1$, we get the following equation for $\epsilon$

[TABLE]

As the left side is smaller than $r$, the integral needs to be smaller than $r\pi/2(d-r)$, which means that $\epsilon$ is close to $1$ for $r\ll d$. This agrees with the intuition that a large part of the eigenvalues of the lower block will be set to zero by projecting the LS onto states. Further details on finding an (approximate) solution to (28) can be found in [66]. In particular, we will approximate $\epsilon$ by the deterministic solution of equation (28) in which $a={\rm Tr}(A)$ is set to zero; indeed, for large $d$ this will have a negligible effect on $\epsilon$ but will allow us to compute $\epsilon$ it in terms of $r$ and $d$ deterministically.

3. Rotation to original basis. Once the cut-off point $q$ has been computed, the projection of the rotated LS estimator $U^{*}\hat{\rho}_{\rm LS}U$ onto states can be written as

[TABLE]

where

[TABLE]

and $[X]_{+}$ denotes the positive part of $X$. The PLS estimator is now obtained by performing the inverse rotation

[TABLE]

We can now estimate the asymptotic risk of the PLS estimator with respect to the chosen loss functions.

Frobenius risk. The mean square error scales as $N^{-1}$ and its rescaled version is

[TABLE]

The contribution from $\tilde{B}$ is

[TABLE]

where the variance of $B$ has been computed as in section 4.3; note that the term $2r\sqrt{(d-r)/N}\epsilon$ vanishes for $N\gg d$. This error is of the order $2rd$ and can be seen as stemming from the uncertainty in estimating the eigenbasis of $\rho_{r}$. For $\tilde{A}$ we write

[TABLE]

Finally, the term $\mathbb{E}\|\tilde{C}\|_{2}^{2}$ is given by the sum of the squares of the remaining eigenvalues which can be approximated using the limit Wigner law as

[TABLE]

Although the last term appears to be of order $(d-r)^{2}$, a careful analysis of the integral [66] shows that it is of lower order than $rd$ due to the fact that $1-\epsilon$ leading term in a Taylor expansion is proportional to $(r/(d-r))^{2/5}$. Therefore, by adding the three terms, we find that the Frobenius risk scales as $6rd/N$. This agrees with the non-asymptotic results of [36], and provides the exact asymptotic constant of the Frobenius rate. By comparing with the lower bound to the asymptotic minimax rate of $2r(d-r)/N$ derived in [17] we find that PLS is optimal for such states, within a constant which is at most 3.

For the $n=7$ atoms state with $r=10$ the Frobenius error is $0.041$, compared to the theoretical one $0.043$. For rank $r=1$ state of 8 atoms with $N=10^{5}$ samples, the Frobenius error was $0.017$, while the theoretical prediction is $0.016$. The simulations results for all ranks of $n=5$ and $n=6$ atoms are presented in panels a) and respectively b) of Figure 3.

Operator norm risk. Unlike the case of the LS estimator, the error matrix of PLS does not approach a Wigner distribution. For this reason, obtaining the asymptotic operator-norm risk seems difficult. However, the following lower bound follows from equation (37)

[TABLE]

Due to the truncation, the matrix $\tilde{C}$ is positive and the largest eigenvalue is approximtely $2\sqrt{d-r}(1-\epsilon)$. On the other hand, $\tilde{A}$ is dominated by the term $2\sqrt{d-r}\epsilon I_{r}$ whose norm is $2\sqrt{d-r}\epsilon$. Since $\epsilon\approx 1$ for large $d$, the lower bound is $2\sqrt{d-r}\epsilon$. This complements the non-asymptotic upper bound of [36] which has rate $O(\sqrt{d})$. Moreover, the lower bound seems to be a good approximation to the actual risk. A simulation with a rank $r=10$ state of $n=7$ atoms and $N=10^{6}$ samples gave an operator-norm error of $0.02$ while the lower bound is $0.01$. For rank $r=1$ state of 8 atoms with $N=10^{5}$ samples, the operator-norm error was $0.12$, while the lower bound is $0.08$.

Trace-norm risk. As in the case of the norm-error, we could not derive the asymptotic expression of the trace-norm risk but we can formulate a lower bound based on the pinching inequality

[TABLE]

Note that $\tilde{A}=A-2\sqrt{d-r}\epsilon I_{r}$, and the variance of the elements of $A$ is of the order 1. Therefore, for $r\ll d$ the shift $2\sqrt{d-r}$ dominates the eigenvalues of $A$ and $\tilde{A}$ is a negative matrix. In first approximation its norm-one is then $2r\sqrt{d-r}\epsilon$. On the other hand, $\tilde{C}$ is positive and its trace can be approximated as

[TABLE]

where in last step we used equation (28) with $a=0$. Therefore the lower bound to the trace-norm error is $4r\sqrt{d-r}\epsilon$. For the $n=7$ atoms state with $r=10$ the trace-norm error was $0.33$ while the lower bound is $0.21$. For rank $r=1$ state of 8 atoms with $N=10^{5}$ samples, the operator-norm error was $0.24$, while the lower bound is $0.16$.

Bures risk. The Bures distance error can be expressed in terms of the blocks $\tilde{A},\tilde{B},\tilde{C}$ as follows

[TABLE]

where we used the Taylor expansion in the last step. The leading term of the Bures risk is then

[TABLE]

where we used the fact that the LS block $A$ has centred distribution. The second order term is

[TABLE]

The simulation results for all ranks states of $n=5$ and $n=6$ atoms are presented panels c) and respectively d) of Figure 3.

9 Comparative numerical simulations

In this section we detail the methodology and results of a general simulation study which compares the performance of the estimators presented in the previous sections for a set of states and against several estimation criteria. The states we consider are rank-$r$ states with a fixed spectrum of $r$ equal eigenvalues of magnitude $1/r$ each, and have randomly generated eigenvectors. This choice is motivated by the fact that such states are arguably harder to estimate among rank-r states (in analogy to the fact that an unbiased coin is harder to estimate than a biased one. Additionally, having such a spectrum allows for a more consistent comparison of the estimators across several ranks.

We generate the above mentioned states for $3$ or $4$ qubits, and for a particular ‘true state’ we simulate a dataset $\mathcal{D}$ of counts from which the state is to be reconstructed. The outcome statistics depend on a few variables that we may vary, namely the type of measurement design (random basis vs Pauli), the number of repetitions per settings $m$, and in the case of the random basis measurements the total number $k$ of measured settings. This allows us to study the performance of the estimators across several different combinations of variables: types of states, ranks, measurement design, number of repetitions per setting $m$, the total number of settings $k$ and the number of ions $n$.

We will present plots of the estimated mean errors $\mathbb{E}\left[D(\hat{\rho},\rho)\right]$ of the estimators for the equal eigenvalues states. For each given rank $r$ and number of qubits $n$, we generate a state with equal eigenvalues $\left(\frac{1}{r},\ldots,\frac{1}{r},0,\ldots,0\right)$ and random eigenbasis. Then for each choice of measurement design and values of $k$ and $m$, the several estimates of the true state are evaluated. The error of each resulting estimate is computed using all the error functions listed in Table 1, and the corresponding mean errors are estimated from 100 different runs of the experiment.

In order to make the results of the simulation study more accessible, we have made all plots for 3 and 4 atoms simulation available online via an interactive Rshiny application at this address: https://rudhacharya.shinyapps.io/plots/, while a selection is presented in the paper.

9.1 Squared Frobenius norm

Figures 4 and 5 show the (estimated) Frobenius risk (mean square error) for states 4 qubits measured with the Pauli and the random basis design, respectively. The states are chosen randomly from the family of rank-$r$ states with equal non-zero eigenvalues. In both cases we note that the Frobenius risk of the LS estimator has no significant dependence on the rank of the true state, and its performance is poor for small rank states. In contrast, the remaining estimators all show a scaling of the Frobenius risk with the rank of the true state. We also note that the performance of several of the estimators matches well with the Fisher-predicted risk. This is remarkable as the latter (6) was defined for a rank-$r$ parameterisation of states, while none of the estimators have any prior knowledge of the rank.

We further note that for relatively small values of $N=m\times k$ the TLS, TGLS, posLS, posGLS, ML estimators significantly better that LS and GLS at higher ranks, while the errors approach each other for larger $N$, cf. Figures 4 a) and 5 a) versus 4 b) and 5 b). This reflects the fact that in the low $N$ case the asymptotic regime has not been reached and constrained estimators have an advantage even for full rank states which lie in the interior of the parameter space. Indeed for an eigenvalue $\lambda$ of order $1/d$ and relatively small values of $N$, the unconstrained estimates $\hat{\lambda}$ will have a standard deviation of order $\sqrt{d/N}$ which may be comparable or larger than the magnitude of the eigenvalues themselves. Therefore without the constraint of positivity such estimators may produce estimates with $\hat{\lambda}<0$. In contrast, for the constrained estimators the positivity constraint provides additional information when the eigenvalues are small. For large values of $N$ however, the uncertainty in the eigenvalues is very small and the Fisher risk acts as a lower bound for all of the estimators.

Across both the Pauli and the random measurement designs we note that the performance of the posGLS and the ML estimators is very similar, and for large $m$ almost identical. This confirms our asymptotic analysis which shows that the posGLS and the ML estimators are equivalent in the limit of large $m$, cf. section 7.2.

9.2 Bures and Hellinger distance

As the Bures distance $D_{B}(\hat{\rho},\rho)^{2}$ is well defined only over density matrices, we plot the mean Bures errors only for the ML, TLS, TGLS, posLS , posGLS estimators. Figures 6 and 7 show the mean Bures errors for different sample sizes in the Pauli, and respectively random bases measurement design. For comparison, in Figures 8 and 9 we also plot the corresponding average Hellinger errors $D_{H}(\hat{\boldsymbol{\lambda}},\boldsymbol{\lambda})^{2}$, cf. Table 1.

We note that the behaviour of the Hellinger errors is very similar to that of the Bures errors, with a better match for larger values of $N$. To give some intuition about this, we look at what happens in the case of qubits. The Bures distance between neighbouring qubit states can be approximated by the sum of Hellinger distance and a quadratic form in the parameters of the unitary connecting the eigenbases of the two states [4]

[TABLE]

where $\Phi$ is the angle between the Bloch vectors of $\rho$ and $\hat{\rho}$. In a non-adaptive measurement scenario such as those considered here, the Bloch vector parameters can be estimated at rate $1/\sqrt{N}$ which means that the second term on the right side of (39) is always of the order $1/N$. However, for states which are very pure ($\lambda\approx 0$) the Hellinger component has the dominant contribution to the Bures distance, and is responsible for the “non-standard” scaling of $1/\sqrt{N}$ in the minimax risk [4]. The simulation results indicate that a similar phenomenon may occur in higher dimensional systems. For full rank states, both the Bures and the Hellinger distance have a quadratic expansion, and in this sense a relation similar to (39) can be derived by splitting the Bures distance into quadratic contributions coming from changes in the eigenvalues and small basis rotations respectively, see also [61]. Alternatively, the block-matrix techniques used for analysing the Bures risk in section 8.1 can be extended to rank deficient states of arbitrary spectrum to show that the leading $1/\sqrt{N}$ contribution comes from the Hellinger $D_{H}(\boldsymbol{\lambda},\hat{\boldsymbol{\lambda}})^{2}.$

Another noticeable feature across both the Pauli and the random bases measurement design is that for large $N$ the mean errors are seen to be larger for states of middling ranks than for the full rank states, see Figures 6b, 7b. This is however not true for smaller values of $N$, as shown in Figures 6a, 7a. More precisely, for large $N$ we see a steep increase from pure states to low rank states, followed by an almost linear decrease down to the full rank state. A similar behaviour has been uncovered in the analysis of the PLS estimator for covariant measurements in section 8.1. There, we found that for large $d$ and low rank $r$, the eigenvalues distribution of the error block $C$ of the LS estimator converges to the Wigner distribution, which allows us to compute the leading orders of the Bures risk of the PLS estimator. Since the covariant measurement arises in the large $m$ limit of the random basis measurement [3], it is therefore expected that the risks behave similarly in the two cases. Our simulations indicate that the mechanism governing the asymptotics of the Bures risk seems to be robust with respect to the details of the measurement; although the Pauli basis measurement is not expected to produce an LS estimator whose error matrix is Wigner distributed, the Bures and Hellinger risks have similar characteristics as those of the covariant measurement. These findings are in line with those of [66], and are worthy of further theoretical investigation.

While the Bures risk for rank deficient states scales as $1/\sqrt{N}$, for full rank states the Bures distance is locally quadratic, and standard asymptotic results imply that the convergence rate is in this case $1/N$, cf. section 3.1. In general,

[TABLE]

with weight matrix $G_{B}(\boldsymbol{\theta})=F(\rho_{\boldsymbol{\theta}})/4$, where $F(\rho_{\boldsymbol{\theta}})$ is the quantum Fisher information. For the maximally mixed state $\rho_{r=d}$, and the parametrisation (4), the latter has 2 uncorrelated blocks [3] corresponding to off-diagonal parameters

[TABLE]

and respectively diagonal parameters

[TABLE]

On the other hand, the classical average Fisher information for random basis measurements is given by $I=F/(d+1)$ [17]. Since the ML estimator is asymptotically normal with variance $I(\boldsymbol{\theta})^{-1}/N$, we have

[TABLE]

For 3 qubits with $N=100\cdot 1000$ samples, and 4 qubits with $N=100\cdot 1000$, the asymptotic predictions for Bures errors are $0.0014$ and $0.01$ which match closely the simulations results [1].

9.3 Trace-norm distance

The risks for the trace-norm distance exhibit a similar behaviour to those of the Frobenius distance, as illustrated in Figures 10 and 11. A noticeable feature is that all constrained estimators have close risks for large sample sizes in both the Pauli and the random bases setting.

10 Conclusions

In this paper we studied theoretical an practical aspects of quantum tomography methods. The unifying theme is that each estimator can be seen as a projection of the data onto a parameter space with respect to an appropriate metric. We considered estimators without positivity constraints (unconstrained maximum likelihood, least squares, generalised least squares) and with positivity constraints (maximum likelihood, positive least squares, thresholded least squares and projected least squares), and investigated the relationships between different estimators. While no estimator makes use of the state’s rank (which is assumed to be unknown) the constrained estimators have significantly lower errors than the unconstrained estimators, for low rank states. To better understand this behaviour we derived new asymptotic error rates for the least squares estimator and for the projected least squares estimators, for a class of given rank states and covariant measurements. These results capture the exact rate dependence on rank and dimension and complement non-asymptotic concentration bounds of [17, 72, 36], showing that PLS has strong optimality properties; for instance the leading contribution to the Frobenius risk is $6rd$ which is ‘almost optimal’, in that it is only 3 times larger than the minimax lower bound established in [17], which assumes that the rank is known. Our analysis has uncovered an interesting behaviour of the Bures error, which increases with respect to rank for small ranks and then decreases towards full rank states. The extensive simulations study shows that this behaviour (as well as the monotonic increase of other errors) is robust with respect to the measurement design.

Computationally, maximum likelihood and positive least squares involve an optimisation over states and are significantly slower than projected least squares which only requires the diagonalisation of the least squares matrix followed by a simple truncation procedure. Our results confirm and strengthen those of [36] and show that projected least squares is an attractive alternative to maximum likelihood, which is routinely used in practice. Additional improvements can be achieved by using generalised least squares as starting point, in a two steps procedure.

An interesting and practically relevant open question is whether any of the ‘fast’ estimators analysed here is statistically ‘optimal’ for more realistic measurements such as the Pauli bases. More generally, one can ask whether these methods can be adapted to non-informationally complete measurement scenarios, and other physically motivated lower dimensional statistical models.

Acknowledgements. We thank Jonas Kahn, Richard Kueng and Ian Dryden for fruitful discussions during the preparation of this manuscript.

Appendix A Appendix

A.1 Fisher information for GLS

Here we prove the equality $I(\rho|\mathscr{S})=I_{\rm GLS}$ stated in section 5.1. Note that by definition

[TABLE]

and

[TABLE]

We will show that

[TABLE]

where $c,d$ are some constants and $|p_{\rho}^{-1}\rangle\in\mathbb{C}^{k\cdot d}$ is the vector whose entries are the inverses of measurement probabilities. Indeed, since $J^{*}X^{*}|\mathbf{1}\rangle=0$, equation (41) implies (40).

Now $V(V^{*}\Omega V)^{-1}V^{*}$ is the pseudo-inverse of $\Omega$ and since $V^{*}|\mathbf{1}\rangle=0$, it satisfies

[TABLE]

To compute the inverse on the right side we use the definition (17) of $\Omega$ and apply the Sherman-Morrison formula [68]

[TABLE]

This gives

[TABLE]

where $B={\rm Diag}(p_{\rho})-|\mathbf{1}\rangle\langle\mathbf{1}|$. By applying the Sherman-Morrison formula again we get

[TABLE]

where $q=\sum_{{\bf o},{\bf s}}p^{-1}_{\rho}({\bf o}|{\bf s})$. Moreover,

[TABLE]

By plugging the last three expresssions into (42) we get

[TABLE]

which concludes the proof of (41) and of $I(\rho|\mathscr{S})=I_{\rm GLS}$.

A.2 Implementations of estimators

Here we list certain practical details about the implementation of the estimators described in the paper. In particular we discuss how we compute the covariance matrix estimator $\hat{\tilde{\Omega}}$ used in the GLS, TGLS, posGLS estimators, and we also describe the cross-validation procedure we use to select a constant $C$ for the thresholded in the TLS and TGLS estimators.

The covariance matrix estimator $\hat{\tilde{\Omega}}$ for the generalised estimators (GLS, TGLS and posGLS) is computed as follows. Given a dataset $\mathcal{D}$, we first obtain the LS estimate $\hat{\rho}_{\rm LS}$ and then construct the TLS estimate (see Algorithm 1) with threshold $\nu=0$. From this we obtain an estimate of the probabilities $\hat{p}(\boldsymbol{o}|\boldsymbol{s})={\rm Tr}\left[\hat{\rho}_{\rm TLS}P^{\boldsymbol{s}}_{\boldsymbol{o}}\right]$. The matrix $\hat{\tilde{\Omega}}$ is then constructed from these estimated probabilities via (17). The generalised estimates (GLS/TGLS/posGLS) are then evaluated using $\hat{\tilde{\Omega}}$ and the same dataset $\mathcal{D}$. 2. 2.

As mentioned briefly in section 6, the threshold for the TLS and TGLS estimators is selected using cross-validation. We describe this cross-validation method below [17].

• •

For a particular number of repetitions per setting $m$, we simulate data in 5 independent batches $m/5$ repetitions per setting in each batch. and we denote the corresponding datasets as $\mathcal{D}_{1},\ldots,\mathcal{D}_{5}$. The total dataset of counts is the sum $\mathcal{D}=\sum_{i=1}^{5}\mathcal{D}_{i}$.

• •

We choose a vector of constants $C$ forming a mesh over the interval $[0,1]$. For each value of $C$, and for each $j\in\{1,\ldots,5\}$ we compute the following estimators. We hold out the dataset $\mathcal{D}_{j}$, and compute the TLS/TGLS estimate $\hat{\rho}^{-j}_{\rm T(G)LS}(C)$ for the dataset $\mathcal{D}_{-j}=\sum_{i\neq j}\mathcal{D}_{i}$, with threshold $\nu=C\sqrt{\frac{4}{m}\log{2^{n+1}}}$, cf. [17]. For each $\mathcal{D}_{j}$ the LS estimate $\hat{\rho}_{\rm LS}^{j}$ is also evaluated.

• •

For all values of $C$, the empirical discrepancy is evaluated for a choice of error function $D(\hat{\rho},\rho)$

[TABLE]

• •

This function $CV(C)$ is then minimised over all values $C$

[TABLE]

this gives an estimate for the holding constant, which is then used to evaluate the TLS or the TGLS estimators with threshold $\nu=\hat{C}_{D}\sqrt{\frac{4}{m}\log{2^{n+1}}}$.

Notice that the cross-validation procedure picks different constants for different choices of the error function. An important caveat here is that the Bures distance is not defined for the LS estimates $\hat{\rho}^{j}_{\rm LS}$, and therefore the procedure above cannot apply. Instead in the simulations we estimate the thresholding constant $\hat{C}_{D_{B}}$ using the ML estimate as

[TABLE]