Quantum Tomography by Regularized Linear Regression

TL;DR

This paper introduces regularized linear regression methods for quantum state tomography, improving estimation accuracy through bias-variance trade-offs and providing asymptotic optimality with practical tuning strategies.

Contribution

It develops a unified framework for regularized linear regression in quantum tomography, including explicit formulas for optimal regularization and asymptotic risk minimization.

Findings

Weighted least squares estimate effectively reconstructs quantum states.

Regularization reduces mean-square error in state estimation.

Proposed tuning methods achieve asymptotic optimality.

Abstract

In this paper, we study extended linear regression approaches for quantum state tomography based on regularization techniques. For unknown quantum states represented by density matrices, performing measurements under certain basis yields random outcomes, from which a classical linear regression model can be established. First of all, for complete or over-complete measurement bases, we show that the empirical data can be utilized for the construction of a weighted least squares estimate (LSE) for quantum tomography. Taking into consideration the trace-one condition, a constrained weighted LSE can be explicitly computed, being the optimal unbiased estimation among all linear estimators. Next, for general measurement bases, we show that $\ell_2$-regularization with proper regularization gain provides even lower mean-square error under a cost in bias. The regularization parameter is tuned by two estimators in terms of a risk characterization. Finally, a concise and unified formula is established for the regularization parameter with complete measurement basis under an equivalent regression model, which proves that the proposed tuning estimators are asymptotically optimal as the number of samples grows to infinity under the risk metric. Additionally, numerical examples are provided to validate the established results.

Figures

Equations

$$\rho =\sum _{i=1}^{d^2}{\theta }_i{\mathsf{B}}_i$$

$${\mathsf{E}}_m=\sum _{i=1}^{d^2}{\beta }_{mi}{\mathsf{B}}_i$$

$$\mathbf{p}=\mathbf{A}\mathbf{\theta }.$$

$${\hat{p}}_m=\frac{\#m}{n}$$

$${\hat{p}}_m={\mathbf{\beta }}_m^{\top }\theta +e_m,$$

$$\mathbf{y}=\mathbf{A}\mathbf{\theta }+\mathbf{e}$$

$$e_m={\hat{p}}_m-p_m=\frac{{\sum }_{l=1}^n{b}_l^{(m)}}{n}-p_m=\sum _{l=1}^n\frac{b_l^{(m)}-p_m}{n}.$$

$$\mathbb{E}(e_m)=0$$

$$\mathbb{V}(e_m)=\mathbb{E}(e_m{)}^2=(p_m-p_m^2)/n.$$

$$\mathbb{P}\left(e_m={\left(\frac{1-p_m}{n}\right)}^k{\left(-\frac{p_m}{n}\right)}^{n-k}\right)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p_m^k(1-p_m{)}^{n-k}.$$

$${\bar{p}}_m=\frac{{\#}^{\prime }m}{n}$$

$${\bar{p}}_m={\mathbf{\beta }}_m^{\top }\theta +{\bar{e}}_m,m=1,\dots ,M.$$

$${\hat{\mathbf{\theta }}}^{\mathrm{LS}}$$

$$=({\mathbf{A}}^{\top }\mathbf{A}{)}^{-1}{\mathbf{A}}^{\top }\mathbf{y}$$

$$\displaystyle\mathbb{MSE}\big{(}\widehat{\bm{\theta}}^{\rm LS}\big{)}$$

$$=({\mathbf{A}}^{\top }\mathbf{A}{)}^{-1}{\mathbf{A}}^{\top }\mathrm{P}\mathbf{A}({\mathbf{A}}^{\top }\mathbf{A}{)}^{-1}$$

$${\hat{\mathbf{\theta }}}^{\mathrm{WLS}}$$

$$=({\mathbf{A}}^{\top }\mathrm{W}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\top }\mathrm{W}\mathbf{y}$$

$$\displaystyle\mathbb{MSE}\big{(}\widehat{\bm{\theta}}^{\rm WLS}\big{)}=(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}.$$

$$\mathbb{MSE}(\hat{\mathbf{\theta }})\ge \mathbb{MSE}({\hat{\mathbf{\theta }}}^{\mathrm{WLS}}).$$

$${\hat{\mathbf{\theta }}}^{\mathrm{AWLS}}=({\mathbf{A}}^{\top }\hat{\mathrm{W}}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\top }\hat{\mathrm{W}}\mathbf{y},$$

$$\displaystyle\widehat{\mathrm{W}}=n\cdot\mathrm{diag}\big{(}[1/(\widehat{p}_{1}-\widehat{p}_{1}^{2}),\cdots,1/(\widehat{p}_{M}-\widehat{p}_{M}^{2})]\big{)}$$

$${\hat{\mathbf{\theta }}}^{\mathrm{AWLS}}-{\hat{\mathbf{\theta }}}^{\mathrm{WLS}}$$

$$\displaystyle=\big{(}(\mathbf{A}^{\top}\widehat{\mathrm{W}}\mathbf{A})^{-1}-(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}\big{)}\mathbf{A}^{\top}\widehat{\mathrm{W}}\mathbf{e}+(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}\mathbf{A}^{\top}\big{(}\widehat{\mathrm{W}}-\mathrm{W}\big{)}\mathbf{e}$$

$$\displaystyle=O_{p}(1/\sqrt{n})(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathrm{W}\big{(}1+O_{p}(1/\sqrt{n})\big{)}\mathbf{e}+(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}\mathbf{A}^{\top}O_{p}\big{(}1/\sqrt{n}\big{)}\mathrm{W}\mathbf{e}$$

$$=O_p(1/\sqrt{n})({\mathbf{A}}^{\top }\mathrm{W}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\top }\mathrm{W}\mathbf{e}$$

$$e_m=O_p(1/\sqrt{n}),1\le m\le M$$

$$\displaystyle\widehat{\mathrm{W}}=\mathrm{W}\big{(}1+O_{p}(1/\sqrt{n})\big{)}$$

$$\displaystyle\widehat{\mathrm{W}}-\mathrm{W}=O_{p}\big{(}1/\sqrt{n}\big{)}\mathrm{W}.$$

$$\mathrm{Tr}(\rho )=1.$$

$${\mathbf{\theta }}^{\top }\mathrm{Tr}(\mathsf{B})=1$$

$$\mathrm{Tr}(\mathsf{B})\stackrel{△}{=}[\mathrm{Tr}({\mathsf{B}}_1),\cdots ,\mathrm{Tr}({\mathsf{B}}_{d^2}){]}^{\top }.$$

$${\hat{\mathbf{\theta }}}^{\mathrm{CLS}}=\underset{{\mathbf{\theta }}^{\top }\mathrm{Tr}(\mathsf{B})=1}{\mathrm{argmin}}(\mathbf{y}-\mathbf{A}\mathbf{\theta }{)}^{\top }(\mathbf{y}-\mathbf{A}\mathbf{\theta }).$$

$$\displaystyle\widehat{\bm{\theta}}^{\rm CLS}=\widehat{\bm{\theta}}^{\rm LS}-\frac{C\operatorname*{Tr}(\mathsf{B})}{\operatorname*{Tr}(\mathsf{B})^{\top}C\operatorname*{Tr}(\mathsf{B})}\big{(}\operatorname*{Tr}(\mathsf{B})^{\top}\widehat{\bm{\theta}}^{\rm LS}-1\big{)}$$

Full text

1 Introduction

Since Feynman pointed out the possibility of using quantum resources to carry out computation in the early 1980s, significant progresses have been made in both the theoretical understanding and the real-world implementations for computing and communication mechanisms based on quantum states (Nielsen & Chuang, 2001). Underpinning such efforts lies in the development of quantum tomography (Wootters & Fields, 1989; James et al., 2001; Artiles et al., 2005; Senko et al., 2014), where reliable quantum state estimation (Bisio et al., 2009; Blume-Kohout, 2010; Teo et al., 2011) and system identification methods (Bonnabel et al., 2009; Leghtas et al., 2012; Wang et al., 2018; Xue et al., 2018; Wang et al., 2019) provide basic assurance for the validity of quantum systems that we intend to work on. The fundamental quantum measurement postulate indicates that any form of quantum information probe would have an inherent probabilistic nature. The exponentially growing complexity of quantum systems along with the increasing scale further adds to the challenging reality: only partial information can be made available via measurements for uncertain quantum systems; processing the measurement data faces enormously high computation barrier for large-scale quantum systems.

One primary task of quantum tomography is to determine an unknown quantum state from a number of identical copies (Bisio et al., 2009; Blume-Kohout, 2010; Teo et al., 2011). Performing measurement on those copies along certain observables, i.e., measurement bases, yields independent realizations of some hidden random variable whose statistics encode the quantum state and the observable. Therefore, utilizing the outcomes of the measurements we can build estimations of the unknown quantum state since the observables are known (which can be selected and designed). Apparently the choice of the estimation method is not unique since the estimation error metrics can be characterized by different metrics, and the resulting computational feasibility also raises constraint on the potentially viable estimation approaches. Therefore, there is often a tradeoff between estimation quality and computational efficiency (Smolin et al., 2012).

Linear regression, as a universal estimator (Ljung, 1999), becomes a natural and important quantum state tomography approach due to its simplicity and practicability. A thorough comparison was made in Qi et al. (2013) between linear regression method and maximum likelihood estimation for quantum state tomography under complete or over-complete measurement bases. Linear regression was also applied to quantum tomography with incomplete measurement bases for low-rank, e.g., Gross et al. (2010); Gross (2011); Alquier et al. (2013) or sparse, e.g., Cai et al. (2016) quantum states in view of the insights from compressed sensing, where a small number of measurement bases was proven to be enough for the recovery of a high-dimensional quantum state with high probability as the dimension increases. Recently, linear regression method was also generalized to the adaptive measurement case where selection of the measurement basis depends on the previous measurement outcomes (Qi et al., 2017).

In this paper, we study the role of regularization for linear regression-based quantum state tomography. In recent years, the power of regularization has been well noted in the literature of optimization, machine learning, and system identification, e.g., Shalev-Shwartz (2012); Goodfellow et al. (2016); Chiuso (2016) for the purpose of avoiding overfitting in empirical learning. Noting any quantum state can be represented as a trace-one positive Hermitian density matrix, which is of low rank if it is a combination of a small number of pure states, we establish the following results.

• •

For complete or over-complete measurement basis, the empirical data can be utilized for constructing of a weighted least squares estimate (LSE) for quantum tomography. The weighted LSE provides reduced mean-square error compared to standard LSE. Taking into consideration the trace-one condition, the constrained weighted LSE can be explicitly computed, which is the optimal unbiased estimation that is linear in the measurement data.

• •

For any (complete, over-complete, or under-complete) measurement basis, a closed form solution is established for tomography with $\ell_{2}$-regularized weighted linear regression. It is shown that with proper regularization parameter, this regularized regression always provides even lower mean-square error subject to, of course, a price of additional bias.

• •

The regularization parameter can be further optimized subject to a risk characterization. An explicit formula is established for the regularization parameter under an equivalent regression model, which proves that the proposed tuning estimators are asymptotically optimal for complete bases as the number of samples grows to infinity for the risk metric.

Numerical examples are provided for the validation of the established theoretical results, which confirms the potential usefulness of the proposed linear regression methods in quantum state tomography.

The remainder of the paper is organized as follows. In Section 2, we present the standard linear regression model for quantum state tomography, and review some preliminary knowledge on the underlying rationale. In Section 3, we present the extended quantum tomography methods based on weighted LSE, constrained weighted LSE, and constrained regularized LSE, respectively, whose performances in terms of mean-square error are thoroughly investigated. Section 4 further presents the asymptotically optimal regularization gain under an equivalent model. Numerical examples are presented in Section 5, and finally, some concluding remarks are given in Section 6.

2 Problem Definition and Preliminaries

2.1 Linear Regression for Quantum State Tomography

Let $\mathcal{H}$ be a $d$-dimensional Hilbert space that characterizes the state space of a quantum system. Denote the space of linear Hermitian operators over $\mathcal{H}$ by $\mathcal{L}(\mathcal{H})$. Suppose that \big{\{}\mathsf{B}_{i}\big{\}}_{i=1}^{d^{2}} is an orthonormal basis of $\mathcal{L}(\mathcal{H})$ with $\operatorname*{Tr}(\mathsf{B}_{i}^{\dagger}\mathsf{B}_{j})=\delta_{ij}$ and $\mathsf{B}_{i}^{\dagger}=\mathsf{B}_{i}$, where $\operatorname*{Tr}(\cdot)$ means the trace of a square matrix, $(\cdot)^{\dagger}$ represents the Hermitian conjugate of a complex matrix, and $\delta_{ij}$ is the Kronecker function. A quantum state $\rho$ as a density operator over $\mathcal{H}$ can then be expressed by

[TABLE]

where $\theta_{i}=\operatorname*{Tr}(\rho\mathsf{B}_{i})\in\mathbb{R}$ is the coordinate of $\rho$ under the given basis \big{\{}\mathsf{B}_{i}\big{\}}_{i=1}^{d^{2}}. Let there be a positive operator-valued measurement (POVM) over the space $\mathcal{H}$, denoted by $\{\mathsf{M}_{m}\}_{m=1}^{M}$ with $\sum_{m=1}^{M}\mathsf{M}_{m}^{\dagger}\mathsf{M}_{m}=\mathsf{I}$, where $\mathsf{I}$ is the identity operator. Then $\mathsf{E}_{m}\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\mathsf{M}_{m}^{\dagger}\mathsf{M}_{m}$ can be expressed as a linear combination of the orthogonal basis of \big{\{}\mathsf{B}_{i}\big{\}}_{i=1}^{d^{2}}:

[TABLE]

for each $1\leq m\leq M$, where $\beta_{mi}=\operatorname*{Tr}(\mathsf{E}_{m}\mathsf{B}_{i})$. When the quantum state $\rho$ is being measured under the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$, the probability of observing outcome $m$ is $p_{m}=\operatorname*{Tr}(\mathsf{E}_{m}\rho)=\bm{\beta}_{m}^{\top}\bm{\theta}$, where $\bm{\beta}_{m}=[\beta_{m1},\cdots,\beta_{md^{2}}]^{\top}$ and $\bm{\theta}=[\theta_{1},\cdots,\theta_{d^{2}}]^{\top}$. Denoting $\mathbf{p}=[p_{1},\dots,p_{M}]^{\top}\in\mathbb{R}^{M}$, $\mathbf{A}=[\bm{\beta}_{1},\dots,\bm{\beta}_{M}]^{\top}\in\mathbb{R}^{M\times d^{2}}$, we have the following fundamental quantum measurement description in the form of a linear algebraic equation:

[TABLE]

The tomography of an unknown quantum state $\rho$ is therefore equivalent to identifying the vector $\bm{\theta}$, where $\mathbf{A}$ is known and $\mathbf{p}$ is estimated by experimental realizations of measuring $\rho$ from the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$. The POVM can be in general represented under Pauli matrices, see e.g, Cai et al. (2016); Wang (2013).

A standard quantum state tomography process is as follows: (i) Prepare $N=nM$ identical copies of an uncertain quantum state $\rho$; (ii) Perform measurement along each $\mathsf{M}_{m}$ within the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$ independently for $n$ copies; (iii) For each $1\leq m\leq M$, record the number of times that the outcome $\mathsf{M}_{m}$ is observed among those $n$ experiments, denoted by $\#m$, from the $n$ experiments. Then,

[TABLE]

is a natural estimator of the probability $p_{m}$, leading to

[TABLE]

where $e_{m}=\widehat{p}_{m}-p_{m}$ is the estimation error. The distribution of $e_{m}$ depends on the sample size $n$, as it is the sum of $n$ identical and independently distributed (i.i.d.) Bernoulli random variables with mean $p_{m}$. This naturally yields the following linear regression problem:

[TABLE]

with $\mathbf{y}=[\widehat{p}_{1},\cdots,\widehat{p}_{M}]^{\top}$ and $\mathbf{e}=(e_{1},\cdots,e_{M})^{\top}$.

2.2 The Noise Distribution

Define i.i.d. random variables $b_{l}^{(m)}$ for $1\leq l\leq n$, which takes value $1$ with probability $p_{m}$ and [math] with probability $1-p_{m}$. Then there holds

[TABLE]

Note that $(b_{l}^{(m)}-p_{m})/n$ takes value $(1-p_{m})/n$ with probability $p_{m}$ and $-p_{m}/n$ with probability $1-p_{m}$. It follows that

[TABLE]

As a result, the distribution of $e_{m}$ is as follows:

[TABLE]

As $n$ tends to infinity, each $e_{m}$ will converge to a Gaussian random variable with mean [math] and variance $(p_{m}-p_{m}^{2})/n$.

2.3 Simultaneous Measurements

In the tomography process described above, each $\mathsf{M}_{m}$ is separately measured, i.e., a binary outcome is recorded for any copy of $\rho$, where $1$ represents $\mathsf{M}_{m}$, and [math] represents $\mathsf{I}-\mathsf{M}_{m}$. An alternative quantum tomography process can be described based on $n$ copies of $\rho$, where we perform measurement by the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$ collectively. To be precise, the outcome associated with each copy of the quantum state now takes value in $\{1,\dots,m\}$, and then the number of times that the outcome $m$ is observed among those $n$ experiments, denoted by $\#^{\prime}m$, is recorded from the $n$ experiments for each $1\leq m\leq M$. Consequently,

[TABLE]

is still an estimator of the probability $p_{m}$, leading to

[TABLE]

The estimation error $\bar{e}_{m}$ as a random variable has the same distribution of $e_{m}$. However, the $\bar{e}_{m}$ are no longer independent since now $\sum_{m=1}^{n}\bar{p}_{m}=1$ is a sure event. Except for this minor difference, this new formulation of quantum tomography procedure remains the same.

2.4 Standard Least Squares

For the estimation problem (4), the least squares (LS) solution

[TABLE]

is a common choice provided that $\mathbf{A}$ has full column rank. The estimate $\widehat{\bm{\theta}}^{\rm LS}$ admits the following properties:

• •

$\widehat{\bm{\theta}}^{\rm LS}$ is unbiased, namely, \mathbb{E}\big{(}\widehat{\bm{\theta}}^{\rm LS}\big{)}=\bm{\theta};

• •

The mean squared error (MSE) matrix of $\widehat{\bm{\theta}}^{\rm LS}$ is

[TABLE]

where \mathrm{P}=\mathrm{diag}\big{(}[p_{1}-p_{1}^{2},\cdots,p_{M}-p_{M}^{2}]\big{)}/n.

However, this standard least squares neglected the fact that the $e_{m}$ have different variances, although they are all zero mean. As a result, the above covariance is not optimal. Furthermore, the condition that $\mathbf{A}$ be full column rank means the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$ is informationally complete, i.e., any two density operators are distinguishable under the POVM given sufficiently large number of samples. This is not practical for large-scale quantum systems.

3 Regularized Linear Regressions

In this section, we present a few generalizations to standard LSE for the considered quantum state tomography problem, and investigate their performances in terms of mean-square errors.

3.1 Weighted Regression

Noticing $\mathbb{V}(e_{m})=(p_{m}-p_{m}^{2})/n$, we can instead use the following weighted least squares (WLS) estimate

[TABLE]

with \mathrm{W}=\mathrm{P}^{-1}=n\mathrm{diag}\big{(}[1/(p_{1}-p_{1}^{2}),\cdots,1/(p_{M}-p_{M}^{2})]\big{)} penalizing the difference in variances for the noises $e_{m}$. This weighted least square $\widehat{\bm{\theta}}^{\rm WLS}$ continues to be unbiased since \mathbb{E}\big{(}\widehat{\bm{\theta}}^{\rm WLS}\big{)}=\bm{\theta} is easily verifiable and its MSE is

[TABLE]

Suppose $\operatorname*{rank}(\mathbf{A})=d^{2}$ and let $\widehat{\bm{\theta}}$ be any linear unbiased estimate for $\bm{\theta}$. Then we have

[TABLE]

This means it is the best estimator of $\bm{\theta}$ among all unbiased linear estimators in the sense that it achieves the minimal covariance.

In practice, the matrix $\mathrm{W}$ in (12) is unknown and a feasible solution is to use the estimate

[TABLE]

where $\mathrm{W}$ in (12) is replaced by its consistent estimate

[TABLE]

with $\widehat{p}_{m},1\leq m\leq M$ given by (2). In the following, it is shown that the estimate (14) is accurate enough and asymptotically coincides with (12).

For a random sequence $\xi_{n}$, we define $\xi_{n}=O_{p}(a_{n})$ by that $\{\xi_{n}/a_{n}\}$ is bounded in probability, i.e., $\forall\epsilon>0,\exists L>0$ such that $\mathbb{P}(|\xi_{n}/a_{n}|>L)<\epsilon,~{}\forall n$. Then there holds for large $n$ that

[TABLE]

in terms of

[TABLE]

and further

[TABLE]

This means $\widehat{\bm{\theta}}^{\rm AWLS}$ is a consistently practical approximation for the weighted LSE $\widehat{\bm{\theta}}^{\rm WLS}$. Actually, the approximation (15) and resulting conclusions (17) also hold for the following introduced estimators.

3.2 Constrained Weighted Regression

The standard or weighted least squares solutions might lead to estimates that are not legitimate quantum states. In fact, the quantum state has an essential requirement

[TABLE]

This becomes for the model (1) that

[TABLE]

where ${\rm Tr}(\mathsf{B})$ is defined by

[TABLE]

This inspires us to define the constrained least squares (CLS) estimate

[TABLE]

For the estimate (21), we have the following proposition to characterize its property.

Proposition 1

Suppose $\operatorname*{rank}(\mathbf{A})=d^{2}$. The CLS estimate $\widehat{\bm{\theta}}^{\rm CLS}$ has the following closed-form solution

[TABLE]

where $\widehat{\bm{\theta}}^{\rm LS}$ is the least squares estimate given by (10) and $C=(\mathbf{A}^{\top}\mathbf{A})^{-1}$, and its MSE matrix is

[TABLE]

*where $F\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}C-\frac{C\operatorname*{Tr}(\mathsf{B})\operatorname*{Tr}(\mathsf{B})^{\top}C}{\operatorname*{Tr}(\mathsf{B})^{\top}C\operatorname*{Tr}(\mathsf{B})}.$ *

To make the notation simple, we will a little abuse the symbols $F$ and $C$ for different cases in the following.

To reduce the MSE of the estimate (21), we can similarly introduce the constrained weighted least squares (CWLS) estimate

[TABLE]

Theorem 1

Suppose $\operatorname*{rank}(\mathbf{A})=d^{2}$ and $p_{m}\in(0,1)$ for $m=1,\cdots,M$. The estimate $\widehat{\bm{\theta}}^{\rm CWLS}$ can be explicitly written as

[TABLE]

where $\widehat{\bm{\theta}}^{\rm WLS}$ is the WLS estimate (12) and $C=(\mathbf{A}^{\top}\mathrm{W}\mathbf{A})^{-1}$. The resulting MSE

[TABLE]

*where $F\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}C-\frac{C\operatorname*{Tr}(\mathsf{B})\operatorname*{Tr}(\mathsf{B})^{\top}C}{\operatorname*{Tr}(\mathsf{B})^{\top}C\operatorname*{Tr}(\mathsf{B})}$, is optimal in the sense that *

[TABLE]

where $\widehat{\bm{\theta}}$ is any unbiased estimate for $\bm{\theta}$ that is affine in $\mathbf{y}$ and $\bm{\theta}$ satisfies the constraint $\bm{\theta}^{\top}\operatorname*{Tr}(\mathsf{B})=1$.

3.3 Regularized Weighted Regression

Further, we introduce the following weighted regression with $\ell_{2}$-regularization:

[TABLE]

where $\gamma\geq 0$ is a regularization parameter and $\|\cdot\|$ represents the 2-norm of a vector. The motivation for introducing (25) may arise from the following two aspects:

• (i)

When the POVM $\{\mathsf{M}_{m}\}_{m=1}^{M}$ is under-determinate, the matrix $\mathbf{A}$ in (4) might not have full column rank. As a result, the $\widehat{\bm{\theta}}^{\rm LS}$, $\widehat{\bm{\theta}}^{\rm CLS}$, and $\widehat{\bm{\theta}}^{\rm CWLS}$ will all fail to produce a unique estimate to the quantum state. The additional $\ell_{2}$-regularization term in the prediction error function will resolve this non-uniqueness challenge.

• (ii)

In practice, the quantum state $\rho$ is often a combination of some finite number of pure states. As a result, a significant prior knowledge on $\rho$ would be that it is of low rank. Since the rank minimization optimization problem with convex constraints is NP-hard (Recht et al., 2010), the nuclear norm is a common alternative as an approximation of the rank constraint for matrices in various matrix optimization problems. Note that $\rho^{\dagger}\rho$ has the same rank as that of $\rho$. As a result, $\rho^{\dagger}\rho$ is still of low rank and the nuclear norm of $\rho^{\dagger}\rho$ is

[TABLE]

Therefore $\ell_{2}$ regularization can be a good rank penalty as well.

The two aspects are certainly connected in practice, where reconstruction of unknown low-rank quantum state is desried with a small number of measurement basis.

Remark 1

For the positive semidefinite quantum state $\rho$, penalizing the nuclear norm $\rho$ (see, e.g., Gross et al. (2010)) is not quite well-defined because

[TABLE]

where $\sigma_{i}$ and $\lambda_{i}$ are the singular values and eigenvalues of $\rho$, respectively. Note that (25) is essentially the regularized optimization approach adopted in Gross et al. (2010) for the numerical study of quantum state reconstruction problems.

Remark 2

The optimization problem (25) can be equivalently represented as

[TABLE]

where $c>0$ corresponds to $\gamma$. In (27b), it is clear that the $\ell_{2}$ norm of the $\bm{\theta}$ serves as a constraint from the two aspects of motivations for such regularization.

For convenience and consistence of the results displayed in the paper, here we first introduce the regularized weighted least squares (RWLS) estimate

[TABLE]

where the constraint $\bm{\theta}^{\top}\operatorname*{Tr}(\mathsf{B})=1$ is neglected.

The problem (25) also has a closed-from solution, which is stated in the following theorem.

Theorem 2

The optimal weighted regularized quantum state estimate, denoted $\widehat{\bm{\theta}}^{\rm CRWLS}$, as the solution to (25) is given by

[TABLE]

where $C=(\mathbf{A}^{\top}\mathrm{W}\mathbf{A}+\gamma\mathsf{I})^{-1}$. The resulting MSE matrix of $\widehat{\bm{\theta}}^{\rm CRWLS}$ is

[TABLE]

*where $F=C-\frac{C\operatorname*{Tr}(\mathsf{B})\operatorname*{Tr}(\mathsf{B})^{\top}C}{\operatorname*{Tr}(\mathsf{B})^{\top}C\operatorname*{Tr}(\mathsf{B})}$. *

It is worth noting that Theorem 2 does not depend on the rank of $\mathbf{A}$. The next theorem shows that the CRWLS estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$ yields immediate improvement in terms of mean squred error if the regularization parameter $\gamma$ is well chosen.

Theorem 3

There holds

[TABLE]

if 0<\gamma<2/\big{(}\|\bm{\theta}\|^{2}-\frac{1}{\|{\rm Tr}(\mathsf{B})\|^{2}}\big{)}.

Remark 3

There holds from the Cauchy–Schwarz inequality that

[TABLE]

for all quantum states $\rho$. Moreover, when strict equality takes place, there is $\lambda\in\mathbb{R}$ such that $\theta_{i}=\lambda{\rm Tr}(\mathsf{B}_{i})$ for all $i=1,\dots,d^{2}$. As a result,

[TABLE]

which implies $\rho=\lambda\mathsf{I}$, and hence $\lambda$ must be $1/d$. Therefore, we have just established that

[TABLE]

*for all $\rho$ as quantum states except for $\rho=\mathsf{I}/d$. *

Theorem 3 shows that regularization that considers the low rank property of the quantum state $\rho$ can further improve the estimate for the parameter vector $\bm{\theta}$ if we can choose a proper $\gamma$.

Remark 4

Theorem 1 indicates $\widehat{\bm{\theta}}^{\rm CWLS}$ has the smallest MSE among all the unbiased estimate of $\bm{\theta}$ linear with $\mathbf{y}$ while Theorem 3 shows that $\widehat{\bm{\theta}}^{\rm CRWLS}$ has a smaller MSE than $\widehat{\bm{\theta}}^{\rm CWLS}$ even if $\widehat{\bm{\theta}}^{\rm CRWLS}$ is also linear with $\mathbf{y}$. The reason is that regularization introduces a small bias but decreases the variance more such that the total MSE is smaller.

The estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$ is a function of the regularization parameter $\gamma$, the selection of which needs to be determined carefully to achieve a better performance. The essence of tuning $\gamma$ is to choose a proper model complexity for the estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$ given the data. Here we provide a method of tuning $\gamma$ by the measurements based on the risk definition of the estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$. Also, we will prove that the tuning method is asymptotically optimal in the risk sense. For convenience of derivation, rewrite the estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$ as the affine form with respective the output $\mathbf{y}$

[TABLE]

where with $C=(\mathbf{A}^{\top}\mathrm{W}\mathbf{A}+\gamma\mathsf{I})^{-1}$,

[TABLE]

Let us introduce the risk for the estimate $\widehat{\bm{\theta}}^{\rm CRWLS}$ defined by (Rao et al., 2008):

[TABLE]

which is a reference measure to characterize how well the estimate (29) can achieve, namely, gives an upper bound of the estimate (29) in the risk sense (32). Thus the regularization parameter $\gamma$ tuned by the risk (32)

[TABLE]

is the optimal regularization parameter of $\gamma$ for any given data in the risk sense. Unfortunately, the cost function (32) of (33) requires the access to the true parameter $\bm{\theta}$, which is usually unavailable for a system to be identified.

In the following, we use an unbiased estimate for (32) as the cost function of an implementable tuning estimator in terms of data to estimate $\gamma$, which is given by

[TABLE]

It can be verified that the expectation of the cost function (34) over the estimation error $\mathbf{e}$ is exactly the risk (32).

The properties of $\widehat{\gamma}_{R}(\widehat{\bm{\theta}}^{\rm CRWLS})$ and $\widehat{\gamma}_{u}(\widehat{\bm{\theta}}^{\rm CRWLS})$ will be given in the next section under an alternative regression model.

4 An Equivalent Regression Model

Up to now our discussions on the quantum state tomography problem have been around the linear model with an equality constraint:

[TABLE]

In this section, we first present a way of transforming this standing linear regression model into an unconstrained version. Then, under the new but equivalent model we establish some important asymptotic properties of the regularized regression solutions as sample size grows.

4.1 Eliminating Equality Constraint

Let us first construct an orthogonal matrix $\mathbf{Q}$ of size $d^{2}\times d^{2}$ as follows: The first row of $\mathbf{Q}$ is ${\rm Tr}(\mathsf{B})^{\top}/\|{\rm Tr}(\mathsf{B})\|$ and the remaining rows are chosen such that $\mathbf{Q}$ is orthogonal. Thus, we have from (35) that

[TABLE]

where

[TABLE]

The constraint (35b) on $\bm{\theta}$ is forced by the fact that the first element $\beta_{1}$ of $\bm{\beta}$ is $1/\|{\rm Tr}(\mathsf{B})\|$. As a result, the problem (35) is equivalent to the unconstrained linear model

[TABLE]

Clearly, by (37b)

[TABLE]

Thus, regularization (low rank property) on $\rho^{\dagger}\rho$ can also be embedded into $\bm{\alpha}$.

For the model (38), we can produce the corresponding LS, WLS, RWLS estimates. Here, we consider the RWLS estimate for (38) since other estimates (LS, WLS) are special cases by setting $\gamma=0$ and/or $\mathrm{W}=\mathsf{I}$. The RWLS estimate for (38) is defined as

[TABLE]

where

[TABLE]

Intuitively, for an estimate $\widehat{\bm{\alpha}}$ of (38), the vector defined by

[TABLE]

should be the corresponding estimate for (35) and independent of the choice of $\mathbf{Q}$. However, this is not obvious. Now, we intend to show that the hypothesis above is true.

Proposition 2

For any regularization parameter $\gamma\geq 0$, there holds

[TABLE]

*Moreover, *

[TABLE]

Remark 5

When $\gamma=0$, the estimate $\widehat{\bm{\alpha}}^{\rm RWLS}$ is reduced to the WLS estimate of (38). Meanwhile, we have

[TABLE]

for $0<\gamma<2/\bm{\alpha}^{\top}\bm{\alpha}$, an equivalent statement as Theorem 3.

4.2 Asymptotically Optimal Regularization Gain

For the estimate (40), it also needs to well tune the regularization parameter $\gamma$. The risk for the estimate $\widehat{\bm{\alpha}}^{\rm RWLS}$ can be similarly defined as

[TABLE]

and the resulting optimal regularization parameter is

[TABLE]

Let us construct an unbiased estimate

[TABLE]

for (47) and it can straightforwardly check its expectation with respect to $\mathbf{e}$ is $R(\widehat{\bm{\alpha}}^{\rm RWLS})$ up to a constant. Consequently, we propose the tuning estimator for $\gamma$ as

[TABLE]

which gives a way to estimate $\gamma$ directly by the data.

The following proposition illustrates that the tuning estimators (33) and (34) as well as (48) and (50) developed for the constrained model (34) and its unconstrained counterpart (50), respectively, are equivalent.

Proposition 3

There hold

[TABLE]

Denote

[TABLE]

We can establish the asymptotically optimal selection of the regularization parameter $\gamma$ explicitly for the regularized regression estimate of the quantum state in the risk senses (32) and (47).

Theorem 4

Suppose $\operatorname*{rank}(\mathbf{A})=d^{2}$. The limits take place as the sample size $n\xrightarrow{}\infty$ by

[TABLE]

where

[TABLE]

is the asymptotically optimal selection of $\gamma$ for the estimate (29) of the quantum state in the risk senses (32) and (47). Moreover, there hold as $n\xrightarrow{}\infty$

[TABLE]

deterministically and

[TABLE]

in distribution.

Remark 6

Theorem 4 shows that the implementable estimators $\widehat{\gamma}_{u}(\widehat{\bm{\theta}}^{\rm CRWLS})$ and $\widehat{\gamma}_{u}(\widehat{\bm{\alpha}}^{\rm RWLS})$ converges to the asymptotically optimal $\gamma^{\star}$ as the optimal estimators $\widehat{\gamma}_{R}(\widehat{\bm{\theta}}^{\rm CRWLS})$ and $\widehat{\gamma}_{R}(\widehat{\bm{\alpha}}^{\rm RWLS})$ for any finite sample data do. On the other hand, $\widehat{\gamma}_{u}(\widehat{\bm{\theta}}^{\rm CRWLS})$ and $\widehat{\gamma}_{u}(\widehat{\bm{\alpha}}^{\rm RWLS})$ have a slower rate of convergence than that of $\widehat{\gamma}_{R}(\widehat{\bm{\theta}}^{\rm CRWLS})$ and $\widehat{\gamma}_{R}(\widehat{\bm{\alpha}}^{\rm RWLS})$.

5 Numerical Examples

5.1 Overdeterminate Measurement Basis

Example 1

We consider the following quantum Werner state tomography for a two-qubit system as studied in Qi et al. (2013):

[TABLE]

where $|{\boldmath\Psi}^{-}\rangle=({|01\rangle-|10\rangle})/{\sqrt{2}}$ and $q\in[0,1]$ is a parameter associated with the state. Take an orthonormal basis \big{\{}\mathsf{B}_{i}\big{\}}_{i=1}^{16} as

[TABLE]

for $j,k=0,1,2,3$, where

[TABLE]

are the Pauli matrices. Let

[TABLE]

Then

[TABLE]

for $j,k=1,2,\dots,6$ form our measurement basis \big{\{}\mathsf{M}_{m}\big{\}}_{m=1}^{36} with $\mathsf{M}_{m}=|\varphi_{j}\rangle\otimes|\varphi_{k}\rangle$. We can verify that the measurement set \big{\{}\mathsf{M}_{m}\big{\}}_{m=1}^{36} is overcomplete and the matrix $\mathbf{A}=[{\boldmath\beta}_{1},\dots,{\boldmath\beta}_{36}]^{\top}\in\mathbb{R}^{36\times 16}$ has full column rank.

We first sample the parameter $q$ to identify the quantum states $\rho_{q}$, and then for any $\rho_{q}$ carry out the tomography procedure for $1000$ rounds based on $n=110,1100,11000$ copies, respectively. The measurement process is simulated by i.i.d. random numbers $\lambda\in\{1,\dots,36\}$ according to $\mathbb{P}(\lambda=m)=p_{m}=\operatorname*{Tr}(\mathsf{E}_{m}\rho_{q})$. Then for the experimental round $k$, $k=1,\dots,1000$, $\hat{p}_{m}{(k)}$ is recorded as the estimation of $p_{m}$. Based on the $\hat{p}_{m}{(k)}$, we derive the standard, weighted, and constrained weighted estimates $\widehat{\bm{\theta}}^{\mathrm{LS}}{(k)}$, $\widehat{\bm{\theta}}^{\mathrm{WLS}}{(k)}$, $\widehat{\bm{\theta}}^{\mathrm{CWLS}}{(k)}$, and $\widehat{\bm{\theta}}^{\mathrm{CRWLS}}{(k)}$ (with \gamma=1/\big{(}\|\bm{\theta}\|^{2}-\frac{1}{\|{\rm Tr}(\mathsf{B})\|^{2}}\big{)}) according to (10), (12), (23), and (29) respectively, where $\mathrm{W}$ is replaced by its estimate (15). Note that some of the $\hat{p}_{m}$ might inevitably take value zero, and whenever that happens we set $\hat{p}_{m}=10^{-8}$ in the weight matrix $\mathrm{W}$ for the sake of computation. The experimental mean-square error $\mathbb{MSE}_{\rm exp}$ is then computed by averaging the square error from each round of experiments, which are plotted with their theoretical predictions from (11), (13), (24), and (30), for each state and their estimates in Figs. 1–3, where the true $\mathrm{W}$ instead of its estimate is used.

From these figures one easily sees that the experimental estimates are approaching the theoretical ones as the number of samples $n$ grows large for all four estimates, LS, WLS, CWLS, and CRWLS, which validates Theorem 1, Theorem 2, and Theorem 3. For small sample size ($n=110$), the WLS, CWLS, and CRWLS are apparently producing worse experimental mean-square error compared to LS. The reason for that might that the $\widehat{\mathrm{W}}$ constructed from the $\hat{p}_{m}$ may not be accurate enough for approximating the true $\mathrm{W}$. For relatively larger sample size ($n=11000$), the WLS, CWLS, and CRWLS all provide significant improvments compared to LS. It is worth noting that even with small sample size, the CRWLS may lead to drastically reduced error for small $q$, where $\rho_{q}$ tends to be closer to a separable quantum state.

5.2 Small Sample-size and Optimal Regularizer

As we have seen from Example 1, when the sample size $n$ is small, the weighted estimates $\widehat{\bm{\theta}}^{\mathrm{WLS}}$, $\widehat{\bm{\theta}}^{\mathrm{CWLS}}$, and $\widehat{\bm{\theta}}^{\mathrm{CRWLS}}$ may lead to lower accuracy compared to $\widehat{\bm{\theta}}^{\mathrm{LS}}$. In the following example, we show that in this case forcing $\mathrm{W}=\mathrm{I}$ in $\widehat{\bm{\theta}}^{\mathrm{CRWLS}}{(k)}$ to obtain a constrained regularized LS estimate (CRLS) would resolve the issue under the tuning methods in (33) and (34).

Example 2

Consider exactly the same quantum state and tomography setup as in Example 1. Let $\mathrm{W}=\mathrm{I}$ in $\widehat{\bm{\theta}}^{\mathrm{CRWLS}}$ so that we define

[TABLE]

as the unweighted CRLS estimate. The regularization gain $\gamma$ is selected under the optimal value $\widehat{\gamma}_{R}$ from (33) in the risk sense and its unbiased estimate $\widehat{\gamma}_{u}$ from (34), under which for any $\rho_{q}$ we carry out the tomography procedure for $1000$ rounds based on $n=110,1100$ copies, respectively. The resulting experimental mean-square errors $\mathbb{MSE}_{\rm exp}(\widehat{\bm{\theta}}^{\mathrm{CRLS}},\widehat{\gamma}_{R})$ and $\mathbb{MSE}_{\rm exp}(\widehat{\bm{\theta}}^{\mathrm{CRLS}},\widehat{\gamma}_{u})$ are then computed and plotted in Figs. 4–5, respectively, in comparison to the theoretical and experimental MSE of standard LS estimate $\widehat{\bm{\theta}}^{\mathrm{LS}}$.

As we can see from the numerical results, with $n=110$, the regularizer for $\widehat{\bm{\theta}}^{\mathrm{CRLS}}$ significantly improves the estimation accuracy compared to $\widehat{\bm{\theta}}^{\mathrm{LS}}$ under both $\widehat{\gamma}_{R}$ and $\widehat{\gamma}_{u}$. While with $n=1100$, for relatively large $q$, the advantage of $\widehat{\bm{\theta}}^{\mathrm{CRLS}}$ is no longer obvious compared to $\widehat{\bm{\theta}}^{\mathrm{LS}}$ since in this case, the use of the weight ${\rm W}$ becomes essential for the performance.

5.3 Under-determinate Measurement Basis

Example 3

We consider a 6-qubit quantum state. We use the $6$-qubit Pauli matrices to form our basis $\{\mathsf{B}_{j},j=1,\dots,d^{2}\}$ with $d=2^{6}=64$. Let $u=\begin{bmatrix}\frac{\sqrt{3}}{2}&\frac{1}{2}\\ -\frac{1}{2}i&\frac{\sqrt{3}}{2}i\end{bmatrix}$ be a $2\times 2$ unitary matrix. Let

[TABLE]

be three pure states, where $i_{1}=42,i_{2}=8,i_{3}=59,i_{4}=30$. Then $\rho_{p}=\left(u^{\otimes 6}\right)^{\dagger}(\frac{1}{3}|\psi_{1}\rangle\langle\psi_{1}|+\frac{1}{3}|\psi_{2}\rangle\langle\psi_{2}|+\frac{1}{3}|\psi_{3}\rangle\langle\psi_{3}|)u^{\otimes 6}$ is a rank-$3$ density matrix for a $6$-qubit system for all $p\in[0,1]$. Note that $\rho_{p}$ is low-rank but not sparse due to the existence of $u$. We index the set of $6$-fold tensor product of Pauli matrices \big{\{}\sigma_{l_{1}}\otimes\dots\otimes\sigma_{l_{6}}:(l_{1},\dots,l_{6})\in\{1,2,3\}^{6}\big{\}} by $\{P_{j},j=1,\dots,3^{6}\}$. The $P_{j}$’s are of full rank and have eigenvalues $\pm 1$. Denote by $Q_{j\pm}$ the projection onto the eigenspaces of $P_{j}$ with respect to $\pm 1$ respectively. We randomly choose and then fix $\{Q_{j_{1}+},\dots,Q_{j_{200}+}\}$ from $\{Q_{j+},j=1,\dots,3^{6}\}$. Then $\mathcal{M}^{+}\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\{M_{k}=\sqrt{Q_{j_{k}+}/200},k=1,\dots,200\}$ forms an under-complete measurement basis, and the resulting $\mathbf{A}=[\beta_{1},\dots,\beta_{200}]^{\top}\in\mathbb{R}^{200\times 4096}$ becomes under-determinate.

We use $n=1100,11000,110000$ copies for each $\rho_{p}$ and perform independent measurements over each copy along any element in the basis $\mathcal{M}^{+}$. The parameter $p$ is sampled at $p=0,0.1,\dots,1$, and for each $p$, we carry out the tomography procedure for $1000$ rounds. For each round, we set $\gamma=1,10,100,1000$ for $\widehat{\bm{\theta}}^{\mathrm{CRWLS}}$, whose experimental mean-square errors $\mathbb{MSE}_{\mathrm{exp}}$ with respect to the parametrized $p$ are plotted with $n=1100,11000,110000$ in Figs. 6–8, respectively. From these plots we see that the MSE is fundamentally lower bounded by the $\mathcal{M}^{+}$ instead of heavily relying on the sample size $n$. Moreover, the MSE and the risk are not quite sensitive with respect to the regularization gain $\gamma$. It is expected that these estimation results can be improved by utilizing the optimal regularizer selection from (34), but would require significantly higher computation cost.

6 Conclusions

We have studied a series of linear regression methods for quantum state tomography based on regularization. With complete or over-complete measurement bases, the empirical data was shown to be useful for the construction of a weighted LSE from the measurement outcomes of an unknown quantum state. It was proven that the trace-constrained weighted LSE is the optimal unbiased estimation among all linear estimators. For general measurement bases, either complete or incomplete, we showed that $\ell_{2}$-regularization with proper regularization parameter could yield even lower mean-square error under a penalty in bias. An explicit formula was established for the regularization parameter under an equivalent regression model, which shows that the proposed tuning estimators are asymptotically optimal as the number of samples grows to infinity for the risk metric. An interesting future direction lies in regularization-based approach for Hamiltonian identification of quantum dynamical systems.

Appendix A

The proofs of Proposition 1 and Theorem 1 are standard and can be found in the books (Theil, 1971; Amemiya, 1985). So, they are omitted for saving space.

A.1 Proof of Theorem 2

The constrained optimization problem (25) is transformed to an unconstrained one by introducing Lagrange multiplier $\lambda$ and the resulting Lagrange function of (25) is

[TABLE]

Thus the optimal solution $(\widehat{\bm{\theta}}^{\rm CRWLS},\lambda^{*})$ of the problem (25) satisfies the first optimality condition

[TABLE]

It follows from (A.1a) that

[TABLE]

Furthermore, by using (A.1b), we have

[TABLE]

which yields

[TABLE]

Substituting (A.3) into (A.2) obtains

[TABLE]

Next, we compute the MSE matrix of $\widehat{\bm{\theta}}^{\rm CRWLS}$. By the constraint $\operatorname*{Tr}(\mathsf{B})^{\top}\bm{\theta}=1$, we have

[TABLE]

The matrix $F$ has the following properties:

[TABLE]

As a result, the MSE matrix of $\widehat{\bm{\theta}}^{\rm CRWLS}$ is

[TABLE]

This completes the proof.

A.2 Proof of Theorem 3

The proof requires the equivalent model (36) as well as resulting notations and conclusion in Section 4. Denote the orthogonal matrix appearing in (37) by

[TABLE]

Thus by (37) and (45) we have

[TABLE]

When 0<\gamma<2/\big{(}\|\bm{\theta}\|^{2}-\frac{1}{\|{\rm Tr}(\mathsf{B})\|^{2}}\big{)}, there holds

[TABLE]

where the inequality is obtained by (46) in Section 4 and $\widehat{\bm{\theta}}^{\rm CRWLS}(0)$ is exactly $\widehat{\bm{\theta}}^{\rm CWLS}$.

A.3 Proof of Proposition 2 and Remark 5

It follows from Lemma B1 that

[TABLE]

Pre-multiplying $\mathbf{Q}^{\top}$ on both sides of (A.22) proves (45).

By (40), we have

[TABLE]

which derives

[TABLE]

When $0<\gamma<2/\bm{\alpha}^{\top}\bm{\alpha}$, we have $\gamma\bm{\alpha}\bm{\alpha}^{\top}-2\mathsf{I}<0$ and further

[TABLE]

since $\gamma(\mathbf{K}^{\top}\mathrm{W}\mathbf{K})^{-1}$ is always positive definite. Thus we obtain

[TABLE]

if $0<\gamma<2/\bm{\alpha}^{\top}\bm{\alpha}$.

A.4 Proof of Proposition 3

By (A.22), we have

[TABLE]

which means that $R(\widehat{\bm{\theta}}^{\rm CRWLS})=R(\widehat{\bm{\alpha}}^{\rm RWLS})$. So, the assertion (51a) holds.

Similarly, it gives

[TABLE]

The first equation in Lemma (B1) derives

[TABLE]

Combining (A.26) with (A.29) proves (51b).

A.5 Proof of Theorem 4

We first prove the convergence (53). It follows from (A.23) that

[TABLE]

the two terms of which are

[TABLE]

where $\mathbf{S}$ is defined in (B.7). Define

[TABLE]

Thus

[TABLE]

Substituting (A.30) into (A.31) yields

[TABLE]

as $n\xrightarrow{}\infty$. It is clear that

[TABLE]

which can also be expressed by $\bm{\theta}$ and $\Upsilon$ in terms of Lemma B3.

By Lemma B4, the limit $\widehat{\gamma}_{R}(\widehat{\bm{\alpha}}^{\rm RWLS})\xrightarrow{}\gamma^{\star}$ holds as $n\xrightarrow{}\infty$ since the convergence (A.33) is uniform over a compact subset of $[0,+\infty)$ that includes $\gamma^{\star}$.

For convenience of proving (55), denote

[TABLE]

It follows that

[TABLE]

We first consider the decomposition of the first term of the cost function of (50)

[TABLE]

The second term of (A.37) is

[TABLE]

The third term of (A.37) is

[TABLE]

Further, we have

[TABLE]

Define

[TABLE]

Thus

[TABLE]

since \big{(}\mathbf{z}-\mathbf{K}\widehat{\bm{\alpha}}^{\rm WLS}\big{)}^{\top}\mathrm{W}\big{(}\mathbf{z}-\mathbf{K}\widehat{\bm{\alpha}}^{\rm WLS}\big{)} is independent of $\gamma$. Substituting (A.38)–(A.40) into (A.5) turns out

[TABLE]

as $n\xrightarrow{}\infty$ since $\widehat{\bm{\alpha}}^{\rm WLS}\xrightarrow{}\bm{\alpha}$ almost surely as $n\xrightarrow{}\infty$. It follows from Lemma B4 that the limit (55) is true.

It remains to show the rates of convergence (54) and (55).

For this, we first calculate the first and second order derivatives of $\mathscr{C}_{R}(\gamma)$ with respective to $\gamma$ with the help of Lemma B2:

[TABLE]

By using a Taylor expansion, we have

[TABLE]

where $\overline{\gamma}$ is a real number between $\widehat{\gamma}_{R}(\widehat{\bm{\alpha}}^{\rm RWLS})$ and $\gamma^{\star}$, which implies that

[TABLE]

As $n\xrightarrow{}\infty$, we have

[TABLE]

which yields

[TABLE]

For proving (55), the procedure is similar. By Lemma B2, the first and second derivatives of $\mathscr{C}_{U}(\gamma)$ are

[TABLE]

Applying the Taylor expansion obtains

[TABLE]

where $\overline{\gamma}$ is a real number between $\widehat{\gamma}_{u}(\widehat{\bm{\alpha}}^{\rm RWLS})$ and $\gamma^{\star}$. By a straightforward calculation, we have

[TABLE]

by further using the Delta method,

[TABLE]

It follows that

[TABLE]

As a result,

[TABLE]

Appendix B

This appendix contains the technical lemmas used in the proof in Appendix A.

Lemma B1

We have

[TABLE]

Proof. For convenience of proof, we denote

[TABLE]

as well as the column vector $\bm{e_{1}}$ of length $d^{2}-1$ and matrix $\bm{e_{-1}}$ of size $d^{2}$ by $d^{2}-1$ that make up the identity matrix

[TABLE]

Thus the following identities hold.

[TABLE]

The identities (B.6a)–(B.6b) can be verified by a straightforward calculation.

The equation (B.6c) is obtained by choosing the $(2,1)$-block submatrix of the identity

[TABLE]

By using (B.6a)–(B.6b), we have

[TABLE]

It is clear that

[TABLE]

and further it follows from (B.6c) that

[TABLE]

Using (B.6a)–(B.6b) again, one yields

[TABLE]

Thus by (B.6c) we have

[TABLE]

This completes the proof.

Lemma B2

Denote

[TABLE]

We have

[TABLE]

where $\bm{g}$ is any column vector.

The proof is carried out by making use of the matrix differentiation formulas in Petersen & Pedersen (2012, Chapter 2) and is omitted due to limited space.

Lemma B3

We have

[TABLE]

Proof. By letting $\gamma=0$ in Lemma B1, we have

[TABLE]

where

[TABLE]

It follows that

[TABLE]

which derives that

[TABLE]

This completes the proof.

Lemma B4

(Ljung, 1999, Theorem 8.2)*

Assume that*

$\mathscr{C}(\gamma)$* is a deterministic function that is continuous in $\gamma\in\Omega$ and minimized at the point $\gamma^{\star}$, where $\Omega$ is a compact subset of $\mathbb{R}$.* 2. 2)

A sequence of functions $\{\mathscr{C}_{n}(\gamma)\}$ converges to $\mathscr{C}(\gamma)$ almost surely and uniformly in $\Omega$ as $n$ goes to $\infty$.

Then $\widehat{\gamma}_{n}=\operatorname*{arg\,min}_{\gamma\in\Omega}\mathscr{C}_{n}(\gamma)$ converges to $\gamma^{\star}$ almost surely, namely,

[TABLE]