A Condition Analysis of the Weighted Linear Least Squares Problem Using Dual Norms

TL;DR

This paper introduces a new approach to analyze the sensitivity of solutions in weighted linear least squares problems using dual norms, providing explicit condition numbers and efficient estimators.

Contribution

It develops explicit normwise and componentwise condition numbers based on dual operator theory and proposes efficient estimators for these condition numbers.

Findings

Condition numbers accurately bound perturbations in solutions.

Componentwise analysis reveals individual component sensitivities.

Proposed estimators are both accurate and computationally efficient.

Abstract

In this paper, based on the theory of adjoint operators and dual norms, we define condition numbers for a linear solution function of the weighted linear least squares problem. The explicit expressions of the normwise and componentwise condition numbers derived in this paper can be computed at low cost when the dimension of the linear function is low due to dual operator theory. Moreover, we use the augmented system to perform a componentwise perturbation analysis of the solution and residual of the weighted linear least squares problems. We also propose two efficient condition number estimators. Our numerical experiments demonstrate that our condition numbers give accurate perturbation bounds and can reveal the conditioning of individual components of the solution. Our condition number estimators are accurate as well as efficient.

Tables

Table 1. Table 1: Major operations and their costs of Algorithm 1 .

Operations	flops
$\left|A^T\right|\left|𝐝\right|$	$\left(2m-1\right)n$
$\left|A\right|\left|𝐱\right|$	$\left(2n-1\right)m$
$L{𝐡}_i$	$\left(2k-1\right)n$
$Q𝐯$	$\left(2m-1\right)m$
$L_{21}^T𝐯$	$\left(2n-1\right)\left(m-n-1\right)$
${𝐬}^T𝐯$	$2n-1$
$R^{-1}𝐯$ ($R^{-T}𝐯$)	$O\left(n^2\right)$
$L_2𝐯$ ($L_2^T𝐯$)	$O\left(n\left(n+1\right)\right)$
$L_1^{-T}\left({\rho }^{-1}{\mathrm{𝐠𝐬}}^T-L_{21}^T\right)𝐯$	$O\left(mn+{\left(m-n-1\right)}^2\right)$
$\left({\rho }^{-1}{\mathrm{𝐬𝐠}}^T-L_{21}\right)L_1^{-1}𝐯$	$O\left(mn+{\left(m-n-1\right)}^2\right)$

Table 2. Table 2: Comparison of our condition numbers 𝒦 ∞ rel superscript subscript 𝒦 rel \mathcal{K}_{\infty}^{\rm rel} and 𝒦 c rel superscript subscript 𝒦 𝑐 rel \mathcal{K}_{c}^{\rm rel} and their upper bounds 𝒦 ∞ u , rel superscript subscript 𝒦 𝑢 rel \mathcal{K}_{\infty}^{u,\rm rel} and 𝒦 c u superscript subscript 𝒦 𝑐 𝑢 \mathcal{K}_{c}^{u} with their corresponding relative errors ℰ ∞ rel superscript subscript ℰ rel \mathcal{E}_{\infty}^{\rm rel} and ℰ c rel superscript subscript ℰ 𝑐 rel \mathcal{E}_{c}^{\rm rel} .

$ϵ$	$\gamma$	$L$	${ℰ}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{u,\mathrm{rel}}$	${ℰ}_c^{\mathrm{rel}}$	${𝒦}_c^{\mathrm{rel}}$	${𝒦}_c^u$
${10}^{-2}$	${10}^0$	$L_0$	4.0597e$-09$	2.0000e$+00$	2.0000e$+00$	1.3431e$-08$	4.0806e$+02$	4.0813e$+02$
		$L_1$	1.1355e$-08$	2.6350e$+02$	3.4499e$+02$	1.3431e$-08$	4.0806e$+02$	4.0813e$+02$
		$L_2$	4.0597e$-09$	2.0000e$+00$	2.0000e$+00$	4.0597e$-09$	2.0000e$+00$	2.0000e$+00$
${10}^{-2}$	${10}^{-6}$	$L_0$	8.5775e$-09$	2.0004e$+00$	2.0005e$+00$	1.2050e$-06$	5.3548e$+02$	5.3566e$+02$
		$L_1$	9.5653e$-07$	3.3222e$+02$	4.2338e$+02$	1.2050e$-06$	5.3548e$+02$	5.3566e$+02$
		$L_2$	8.5775e$-09$	2.0004e$+00$	2.0004e$+00$	8.5775e$-09$	2.0004e$+00$	2.0004e$+00$
${10}^{-6}$	${10}^0$	$L_0$	1.2311e$-07$	9.1456e$+00$	9.1456e$+00$	1.9008e$-02$	1.4121e$+06$	1.4121e$+06$
		$L_1$	1.9008e$-02$	9.9851e$+05$	1.4121e$+06$	1.9008e$-02$	1.4121e$+06$	1.4121e$+06$
		$L_2$	1.0898e$-09$	2.0000e$+00$	2.0000e$+00$	1.0898e$-09$	2.0000e$+00$	2.0000e$+00$
${10}^{-6}$	${10}^{-6}$	$L_0$	5.2538e$-09$	4.6298e$+03$	4.6298e$+03$	8.1423e$-04$	1.0280e$+09$	1.0280e$+09$
		$L_1$	8.1423e$-04$	7.2690e$+08$	1.0280e$+09$	8.1423e$-04$	1.0280e$+09$	1.0280e$+09$
		$L_2$	5.2538e$-09$	2.0000e$+00$	2.0000e$+00$	5.2538e$-09$	2.0000e$+00$	2.0000e$+00$

Table 3. Table 3: Comparison of our condition numbers 𝒦 ∞ rel superscript subscript 𝒦 rel \mathcal{K}_{\infty}^{\rm rel} and 𝒦 c rel superscript subscript 𝒦 𝑐 rel \mathcal{K}_{c}^{\rm rel} and their upper bounds 𝒦 ∞ u , rel superscript subscript 𝒦 𝑢 rel \mathcal{K}_{\infty}^{u,\rm rel} and 𝒦 c u superscript subscript 𝒦 𝑐 𝑢 \mathcal{K}_{c}^{u} with the corresponding relative errors ℰ ∞ rel superscript subscript ℰ rel \mathcal{E}_{\infty}^{\rm rel} and ℰ c rel superscript subscript ℰ 𝑐 rel \mathcal{E}_{c}^{\rm rel} , when the variances σ i 2 superscript subscript 𝜎 𝑖 2 \sigma_{i}^{2} are evenly spaced between 10 − 4 superscript 10 4 10^{-4} and 5 × 10 − 4 5 superscript 10 4 5\times 10^{-4} .

$L$	${ℰ}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{u,\mathrm{rel}}$	${ℰ}_c^{\mathrm{rel}}$	${𝒦}_c^{\mathrm{rel}}$	${𝒦}_c^u$
$L_0$	5.5085e$-09$	2.7060e$+00$	4.5323e$+00$	1.5998e$-07$	2.6098e$+02$	2.6204e$+02$
$L_1$	6.8834e$-09$	7.5328e$+00$	8.4972e$+00$	3.9877e$-08$	6.2049e$+01$	6.2311e$+01$
$L_2$	1.9161e$-08$	2.8946e$+01$	2.9075e$+01$	1.9161e$-08$	2.8946e$+01$	2.9075e$+01$

Table 4. Table 4: Comparison of our condition numbers 𝒦 ∞ rel superscript subscript 𝒦 rel \mathcal{K}_{\infty}^{\rm rel} and 𝒦 c rel superscript subscript 𝒦 𝑐 rel \mathcal{K}_{c}^{\rm rel} and their upper bounds 𝒦 ∞ u , rel superscript subscript 𝒦 𝑢 rel \mathcal{K}_{\infty}^{u,\rm rel} and 𝒦 c u superscript subscript 𝒦 𝑐 𝑢 \mathcal{K}_{c}^{u} with their corresponding relative errors ℰ ∞ rel superscript subscript ℰ rel \mathcal{E}_{\infty}^{\rm rel} and ℰ c rel superscript subscript ℰ 𝑐 rel \mathcal{E}_{c}^{\rm rel} , when the variances σ i 2 superscript subscript 𝜎 𝑖 2 \sigma_{i}^{2} are evenly spaced between 10 − 4 superscript 10 4 10^{-4} and 10 2 superscript 10 2 10^{2} .

$L$	${ℰ}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{\mathrm{rel}}$	${𝒦}_{\infty }^{u,\mathrm{rel}}$	${ℰ}_c^{\mathrm{rel}}$	${𝒦}_c^{\mathrm{rel}}$	${𝒦}_c^u$
$L_0$	9.2329e$-09$	6.4432e$+00$	8.8514e$+00$	6.0428e$-07$	1.0567e$+03$	1.2003e$+03$
$L_1$	1.2621e$-08$	1.6575e$+01$	1.7419e$+01$	1.3424e$-07$	2.4916e$+02$	2.8492e$+02$
$L_2$	4.8094e$-08$	7.6436e$+01$	8.6338e$+01$	4.8094e$-08$	7.6436e$+01$	8.6338e$+01$

Equations

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\Vert A\mathbf{x}-\mathbf{b}{\Vert }_W^2=\underset{\mathbf{x}}{\min }(A\mathbf{x}-\mathbf{b}{)}^{\mathrm{T}}W(A\mathbf{x}-\mathbf{b}),$$

$$\underset{\mathbf{x}}{\min }\Vert A\mathbf{x}-\mathbf{b}{\Vert }_2^2.$$

$$A^{\mathrm{T}}WA\mathbf{x}=A^{\mathrm{T}}W\mathbf{b}.$$

$$\left[\begin{array}[]{cc}W^{-1}&A\\ A^{\mathrm{T}}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{d}\\ \mathbf{x}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{b}\\ \mathbf{0}\end{array}\right],$$

$$\begin{array}[]{rl}\displaystyle{\min_{y\in\mathbb{R}^{m}}}&\frac{1}{2}\mathbf{y}^{\mathrm{T}}H\mathbf{y}+\mathbf{c}^{\mathrm{T}}\mathbf{y}\\ \hbox{subject to}&A^{\mathrm{T}}\mathbf{y}=\mathbf{b}\\ &\mathbf{y}\geq 0,\end{array}$$

$$\left[\begin{array}[]{cc}M&A\\ A^{\mathrm{T}}&0\end{array}\right]\,\left[\begin{array}[]{r}-\Delta\mathbf{y}\\ \Delta\mathbf{x}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{z}\\ A^{\mathrm{T}}\mathbf{y}-\mathbf{b}\end{array}\right],$$

$$\underset{h\to 0}{\lim }\frac{\Vert f(x+h)-f(x)-A_x(h){\Vert }_W}{\Vert h{\Vert }_V}=0,$$

$${\mathrm{cond}}_{\psi }(a)=\underset{ϵ\to 0}{\lim }\underset{\Vert \delta a\Vert \le ϵ\Vert a\Vert }{\sup }\frac{\Vert \psi (a+\delta a)-\psi (a)\Vert /\Vert \psi (a)\Vert }{\Vert \delta a\Vert /\Vert a\Vert }=\frac{\Vert {\psi }^{\prime }(a)\Vert \Vert a\Vert }{\Vert \psi (a)\Vert },$$

$$AXA=A,XAX=X,(MAX{)}^{\mathrm{T}}=MAX,(NXA{)}^{\mathrm{T}}=NXA,$$

$$\mathbf{x}=(A^{\mathrm{T}}WA{)}^{-1}{A}^{\mathrm{T}}W\mathbf{b}=A_{W,I_n}^{†}\mathbf{b}.$$

$$g(A,\mathbf{b})=L^{\mathrm{T}}(A^{\mathrm{T}}WA{)}^{-1}{A}^{\mathrm{T}}W\mathbf{b}=L^{\mathrm{T}}{A}_{W,I_n}^{†}\mathbf{b},$$

$$⟨\mathbf{y},J\mathbf{x}{⟩}_G=⟨J^{\ast }\mathbf{y},\mathbf{x}{⟩}_E,$$

$$\Vert \mathbf{x}{\Vert }_{E^{\ast }}=\underset{\mathbf{u}\ne 0}{\max }\frac{⟨\mathbf{x},\mathbf{u}{⟩}_E}{\Vert \mathbf{u}{\Vert }_E}.$$

$$\Vert \cdot {\Vert }_{1^{\ast }}=\Vert \cdot {\Vert }_{\infty },\Vert \cdot {\Vert }_{{\infty }^{\ast }}=\Vert \cdot {\Vert }_1,\Vert \cdot {\Vert }_{2^{\ast }}=\Vert \cdot {\Vert }_2.$$

$$\Vert J{\Vert }_{E,G}=\Vert J^{\ast }{\Vert }_{G^{\ast },E^{\ast }}.$$

$$⟨({\mathbf{u}}_1,...,{\mathbf{u}}_p),({\mathbf{v}}_1,...,{\mathbf{v}}_p)⟩=⟨{\mathbf{u}}_1,{\mathbf{v}}_1{⟩}_{E_1}+\cdots +⟨{\mathbf{u}}_p,{\mathbf{v}}_p{⟩}_{E_p}$$

$$\Vert ({\mathbf{u}}_1,...,{\mathbf{u}}_p){\Vert }_v=v(\Vert {\mathbf{u}}_1{\Vert }_{E_1},...,\Vert {\mathbf{u}}_p{\Vert }_{E_p}),$$

$$\Vert ({\mathbf{u}}_1,...,{\mathbf{u}}_p){\Vert }_{v^{\ast }}=v^{\ast }(\Vert {\mathbf{u}}_1{\Vert }_{E_1^{\ast }},...,\Vert {\mathbf{u}}_p{\Vert }_{E_p^{\ast }}).$$

$$\mathcal{K}=\Vert g^{\prime }(\mathbf{y}){\Vert }_{E,G}=\underset{\Vert \mathbf{z}{\Vert }_E=1}{\max }\Vert g^{\prime }(\mathbf{y})\cdot \mathbf{z}{\Vert }_G,$$

$${\mathcal{K}}^{\mathrm{rel}}=\mathcal{K}\frac{\Vert \mathbf{y}{\Vert }_E}{\Vert g(\mathbf{y}){\Vert }_G}$$

$$\mathcal{K}=\underset{\Vert \Delta \mathbf{y}{\Vert }_E=1}{\max }\Vert g^{\prime }(\mathbf{y})\cdot \Delta \mathbf{y}{\Vert }_G=\underset{\Vert \mathbf{z}{\Vert }_{G^{\ast }}=1}{\max }\Vert g^{\prime }(\mathbf{y}{)}^{\ast }\cdot \mathbf{z}{\Vert }_{E^{\ast }}.$$

$$\Vert \Delta \mathbf{y}{\Vert }_c=\min \{\omega ,|\Delta y_i|\le \omega |y_i|,i=1,...,n\}$$

$$\Vert \Delta \mathbf{y}{\Vert }_c=\max \{|\Delta y_i|/|y_i|\}=\Vert (|\Delta y_1|/|y_1|,...,|\Delta y_n|/|y_n|){\Vert }_{\infty },$$

$$\Vert \Delta y_i{\Vert }_{E_i^{\ast }}=\underset{z\ne 0}{\max }\frac{|\Delta y_i\cdot z|}{\Vert z{\Vert }_{E_i}}=\underset{z\ne 0}{\max }\frac{|\Delta y_i\cdot z|}{|z|/|y_i|}=|\Delta y_i||y_i|.$$

$$\Vert \Delta \mathbf{y}{\Vert }_{c^{\ast }}=\Vert (\Vert \Delta y_1{\Vert }_{E^{\ast }},...,\Vert \Delta y_n{\Vert }_{E^{\ast }}){\Vert }_{{\infty }^{\ast }}=\Vert (|\Delta y_1||y_1|,...,|\Delta y_n||y_n|){\Vert }_1.$$

$$\mathcal{K}=\underset{\Vert \mathbf{u}{\Vert }_{G^{\ast }}=1}{\max }\Vert (g^{\prime }(\mathbf{y}){)}^{\ast }\cdot \mathbf{u}{\Vert }_{c^{\ast }},$$

$$J(B,\mathbf{c})$$

$$J^{\ast }(\mathbf{u})=(\mathbf{d}{\mathbf{u}}^{\mathrm{T}}{L}^{\mathrm{T}}((A^{\mathrm{T}}WA{)}^{-1}-(A_{W,I_n}^{†}{)}^{\mathrm{T}}L\mathbf{u}{\mathbf{x}}^{\mathrm{T}},A_{W,I_n}^{†}L\mathbf{u}),$$

$$⟨\mathbf{u},J_1(B)⟩$$

$$⟨\mathbf{u},J_2(\mathbf{c})⟩=⟨\mathbf{u},L^{\mathrm{T}}{A}_{W,I_n}^{†}\mathbf{c}⟩=⟨(A_{W,I_n}^{†}{)}^{\mathrm{T}}L\mathbf{u},\mathbf{c}⟩.$$

Full text

Keywords Weighted least squares, condition number, dual norm, adjoint operator, componentwise perturbation.

Subject Classification 65F20, 65F35

1 Introduction

This paper investigates the condition of the weighted linear least squares problem using adjoint operators and dual norms.

Given $A\in\mathbb{R}^{m\times n}$, $m\geq n$, and $\mathbf{b}\in\mathbb{R}^{m}$, the weighted least squares problem (WLS)

[TABLE]

where $W\in\mathbb{R}^{m\times m}$ is symmetric and positive definite, is a generalization of the standard least squares problem (LS)

[TABLE]

When $A$ is of full column rank, the problem (1) has a unique solution that can be obtained by solving the normal equations

[TABLE]

Alternatively, the solution can be obtained by solving the augmented system:

[TABLE]

where $\mathbf{d}=W(\mathbf{b}-A\mathbf{x})$ is the weighted residual. The systems (3) and (4) are mathematically equivalent for solving $\mathbf{x}$. Methods for solving the weighted least squares problem can found in [21].

The weighted least squares problem (1) can be reduced to the standard least squares problem by the transformations $A\leftarrow W^{1/2}A$ and $\mathbf{b}\leftarrow W^{1/2}\mathbf{b}$, where $W^{1/2}$ is the square root of the symmetric and positive definite $W$, that is $W^{1/2}W^{1/2}=W$. However, the weighted least squares problem arises from applications where $W$ varies. Thus it is efficient to solve the weighted least squares problem rather than transforming it into the standard least squares problem for every $W$. For example, in the interior point method for convex quadratic programming, a convex quadratic programming problem is transformed into

[TABLE]

where $H$ is symmetric and positive semi-definite and $A$ is of full column rank. Introducing the dual variable $\mathbf{x}$ and applying the Newton’s method to the primal-dual equation, we get

[TABLE]

where $M$ is the Hessian matrix and $\Delta\mathbf{y}$ and $\Delta\mathbf{x}$ are respectively the updates for $\mathbf{y}$ and $\mathbf{x}$ in the Newton iteration. The matrix $M$ changes during the iteration while $A$ is fixed. Thus, in this problem, we need to solve a sequence of weighted least squares problems (4) with variable $W$ but fixed $A$.

Another application where the weighted least squares problem arises is the linear regression in statistics. As we know, linear least squares model is commonly used for linear regression assuming the response variables have the same error variance. In practice, however, observations may not be equally reliable. In that case, the weighted least squares model is an improvement of the standard least squares model.

Condition number plays an important role in numerical analysis. It is a measurement of the sensitivity of the solution of a problem to the perturbation of its data. In 1966, Rice [25] presented a general theory of condition based on the Fréchet derivative defined as follows.

Let $V$ and $W$ be two Banach spaces and $U$ an open subset of $V$. Considering an operator $f:\ U\rightarrow W$, if, for an $x\in U$, there exists a bounded linear operator $A_{x}:\ V\rightarrow W$ such that

[TABLE]

then $f$ is said to be Fréchet differentiable at $x$ and $A_{x}$ is called the Fréchet derivative of $f$ at $x$.

Let $\psi:\ \mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$. If $\psi$ is continuous and Fréchet differentiable in a neighbourhood of $a\in\mathbb{R}^{m}$, then, from Rice’s theory, the normwise condition number of $\psi$ at $a$ is defined by

[TABLE]

where $\psi^{\prime}(a)$ is the Fréchet derivative of $\psi$ at $a$. The condition number defined above can be interpreted as the ratio of the relative error in the solution to the relative error in the input data. Clearly, the condition number defined above is norm dependent. Moreover, it is a global measurement, which has some shortcomings. For example, the distribution of the perturbations in data is not represented. Also, as pointed out in [5], for poorly scaled (imbalanced) problems, the error can be overly estimated by the normwise condition number. To alleviate the shortcomings, componentwise perturbation analysis is introduced.

The componentwise error analysis of linear systems can be found in [6, 26]. For the linear least squares, componentwise perturbation analysis and error bounds are given in [4, 19]. In particular, for the full rank linear least squares problem, [14] and [2] present the condition number when the perturbations in both the data and solution are measured by norms. In [3], the perturbation in the data is componentwise, whereas the perturbation in the solution can be either componentwise or normwise, leading to componentwise or mixed condition numbers.

Here is a brief review of some perturbation analyses of the linear least squares problem (LS) and weighted linear least squares problem (weighted LS). For the normwise perturbation analysis, we refer the classical paper [28] and references therein. Cucker, Diao and Wei studied the mixed and componentwise condition numbers for LS in [8]. The flexible normwise condition numbers for LS was introduced in [31]. In [7], Cucker and Diao gave explicit expressions for normwise, mixed and componentwise condition numbers for LS under structured perturbations. Diao and Wei proposed and derived the weighted Frobenius normwise condition number for LS [10]. Recently, Diao et al. [9] studied the normwise, mixed and componentwise condition number for LS involving Kronecker product. For weighed LS, the perturbation analysis can be found in [17, 16, 15] and references therein. Wei and Wang studied the explicit normwise condition numbers under range conditions [32]. Wang et al. derived the results of the Frobenius normwise condition numbers for weighted LS when the coefficient matrix is of full rank [11]. In [22], Li and Sun derived explicit expressions of mixed and componentwise condition numbers of the weighted LS problem. Yang and Wang considered the flexible normwise condition numbers for weighed LS [33].

In this paper, we often use the weighted generalized inverse of $A$ defined as follows for the weighted least squares problem (1). For $A\in\mathbb{R}^{m\times n}$, let $M\in\mathbb{R}^{m\times m}$ and $N\in\mathbb{R}^{n\times n}$ be symmetric and positive definite. If there exists an $X\in\mathbb{R}^{n\times m}$ satisfying the following four equations

[TABLE]

then $X$ is called the weighted generalized inverse of $A$ corresponding to the weight matrices $M$ and $N$ and denoted by $A_{M,N}^{\dagger}$ [30].

When $A$ is of full column rank, the unique solution for (1) is $\mathbf{x}=(A^{\mathrm{T}}WA)^{-1}A^{\mathrm{T}}W\mathbf{b}$. Setting $X=(A^{\mathrm{T}}WA)^{-1}A^{\mathrm{T}}W$, $M=W$ and $N=I_{n}$, it can be verified that $X$ satisfies the four equations in (6), that is, $X$ is the weighted generalized inverse of $A$ corresponding to the weight matrices $W$ and $I_{n}$. Thus, we have $(A^{\mathrm{T}}WA)^{-1}A^{\mathrm{T}}W=A_{W,I_{n}}^{\dagger}$. The unique solution for (1) can be given by

[TABLE]

This paper studies the sensitivity of the solution $\mathbf{x}$ to the perturbation in the data $A$ and $\mathbf{b}$ in the problem (1) by applying the condition number defined in (5). So, corresponding to the function $\psi$ in the definition, we define the following function mapping the data $A$ and $\mathbf{b}$ to the solution $\mathbf{x}$:

[TABLE]

where $L$ is an $n$-by-$k$, $k\leq n$, matrix introduced for the selection of the solution components. For example, when $L=I_{n}$ ($k=n$), all the $n$ components of the solution $\mathbf{x}$ are equally selected. When $L=\mathbf{e}_{1}$ ($k=1$), the first unit vector in $\mathbb{R}^{n}$, then only the first component of the solution is selected.

In this paper, by using adjoint operators and dual norms, we derive the condition numbers of the weighted least square problem, where the perturbation in the data is componentwise whereas the perturbation in the solution is componentwise or normwise. After a brief review of adjoint operators and dual norms and their application to condition number in Section 2, the condition numbers for the weighted least squares problem are presented in Section 3. An error analysis of the augmented system (4) is performed in Section 5. Finally, Section 6 shows our numerical experiment results.

2 Dual Techniques

Consider a linear operator $J:\ E\rightarrow G$, between two Euclidean spaces $E$ and $G$ with the scalar products $\langle\cdot,\cdot\rangle_{E}$ and $\langle\cdot,\cdot\rangle_{G}$ respectively. Denote the corresponding norms by $\|\cdot\|_{E}$ and $\|\cdot\|_{G}$ respectively. We first define adjoint operator and dual norm.

Definition 1

The adjoint operator $J^{*}:\ G\rightarrow E$ of $J$ is defined by

[TABLE]

where $\mathbf{x}\in E$ and $\mathbf{y}\in G$.

Definition 2

The dual norm $\|\cdot\|_{E^{*}}$ of the norm $\|\cdot\|_{E}$ is defined by

[TABLE]

The dual norm $\|\cdot\|_{G^{*}}$ can be defined similarly.

For commonly used vector norms in $\mathbb{R}^{n}$, their dual norms are given by

[TABLE]

For matrices in $\mathbb{R}^{m\times n}$, we consider the norm corresponding to the scalar product $\langle A,B\rangle=\mathrm{trace}(A^{\mathrm{T}}B)$. Thus, we have $\|A\|_{2^{*}}=\|\sigma(A)\|_{1}$ (see [29]), where $\sigma(A)$ is the vector of the singular values of $A$. Since $\mathrm{trace}(A^{\mathrm{T}}A)=\|A\|_{F}^{2}$, we have $\|A\|_{F^{*}}=\|A\|_{F}$.

For linear operators from $E$ to $G$, $\|\cdot\|_{E,G}$ denotes the operator norm induced from $\|\cdot\|_{E}$ and $\|\cdot\|_{G}$. Similarly, for linear operators from $G$ to $E$, the norm induced from the dual norms $\|\cdot\|_{E^{*}}$ and $\|\cdot\|_{G^{*}}$, is denoted by $\|\cdot\|_{G^{*},E^{*}}$.

For the adjoint operators and dual norms, we have the following result [3]:

Lemma 1

[TABLE]

In particular, when $G$ has lower dimension than $E$, we can use the lower dimensional $\|J^{*}\|_{G^{*},E^{*}}$ instead of the higher dimensional $\|J\|_{E,G}$.

Now, we consider the product space $E=E_{1}\times...E_{p}$, where each Euclidean space $E_{i}$ is associated with a scalar product $\langle\cdot,\cdot\rangle_{E_{i}}$ and the corresponding norm $\|\cdot\|_{E_{i}}$, $1\leq i\leq p$. In $E$, we consider the scalar product defined by

[TABLE]

and the corresponding product norm defined by

[TABLE]

where $v$ is an absolute norm on $\mathbb{R}^{p}$, that is, $v(|\mathbf{x}|)=v(\mathbf{x})$, for any $\mathbf{x}\in\mathbb{R}^{p}$, where $|\mathbf{x}|=[|x_{i}|]$, see [20] for details.

Let $v^{*}$ be the dual norm of $v$ and satisfy the usual inner product on $\mathbb{R}^{p}$, we are interested in the dual norm $\|\cdot\|_{v^{*}}$ of the product norm $\|\cdot\|_{v}$ which satisfies the scalar product in $E$. The following result can be found in [3].

Lemma 2

The dual norm of a product norm can be expressed by

[TABLE]

With the necessary background in adjoint operators and dual norms, we apply them to the condition numbers for the weighted least squares problem. We can think of the Euclidean space $E$ with norm $\|\cdot\|_{E}$ as the space of the data in the weighted least squares problem and $G$ with norm $\|\cdot\|_{G}$ as the space of the solution in the weighted least squares problem. Then the function $g$ in (7) is an operator from $E$ to $G$ and the condition number is the measurement of the sensitivity of $g$ to the perturbation in its input data.

Following [25], if $g$ is Fréchet differentiable in a neighborhood of $\mathbf{y}\in E$, then the condition number $\mathcal{K}$ of $g$ at $\mathbf{y}$ is given by

[TABLE]

where $\|\cdot\|_{E,G}$ is the operator norm induced from the norms $\|\cdot\|_{E}$ and $\|\cdot\|_{G}$. If $g$ is nonzero, then we define

[TABLE]

as the relative condition number of $g$ at $\mathbf{y}\in E$. The above definition shows that $\mathcal{K}$ is dependent of the norm of the the linear operator $g^{\prime}(\mathbf{y})$. Applying Lemma 1, we have the following expression of $\mathcal{K}$ in terms of adjoint operator and dual norm:

[TABLE]

Now we consider the componentwise measurement on the data space $E=\mathbb{R}^{n}$. For any given $\mathbf{y}\in E$, $E_{Y}$ denotes the set of all the perturbations $\Delta\mathbf{y}\in\mathbb{R}^{n}$ such that $\Delta y_{i}=0$ when $y_{i}=0$, $1\leq i\leq n$. Thus in the componentwise perturbation analysis, we use the norm

[TABLE]

to measure the perturbation $\Delta\mathbf{y}\in E_{Y}$ of $\mathbf{y}$. We call $\|\cdot\|_{c}$ the componentwise relative norm. Equivalently,

[TABLE]

where $\Delta\mathbf{y}\in E_{Y}$.

Next we investigate the dual norm $\|\cdot\|_{c^{*}}$ of the componentwise norm $\|\cdot\|_{c}$. Let the product space $E$ be $\mathbb{R}^{n}$, each $E_{i}$ be $\mathbb{R}$, and the absolute norm $v$ be $\|\cdot\|_{\infty}$. Setting the norm $\|\Delta y_{i}\|_{E_{i}}$ in $E_{i}$ to $|\Delta y_{i}|/|y_{i}|$ when $y_{i}\neq 0$, from Definition 2, we have the dual norm

[TABLE]

Applying Lemma 2 and (9) and recalling $\|\cdot\|_{\infty^{*}}=\|\cdot\|_{1}$, we get the dual norm

[TABLE]

Due to the condition $\|\Delta\mathbf{y}\|_{E}=1$ in the condition number $\mathcal{K}$ in (8), whether $\Delta\mathbf{y}$ is in $E_{Y}$ or not, the expression of the condition number $\mathcal{K}$ remains valid. Indeed, if $\Delta\mathbf{y}\not\in E_{Y}$, that is, $\Delta y_{i}=0$ for some $i$, then $\|\Delta\mathbf{y}\|_{c}=\infty$. Consequently, such perturbation $\Delta\mathbf{y}$ is excluded from the calculation of $\mathcal{K}$. Following (8), we have the following lemma on the condition number in adjoint operator and dual norm.

Lemma 3

Using the above notations and the componentwise norm defined in (9), the condition number $\mathcal{K}$ can be expressed by

[TABLE]

where $\|\cdot\|_{c^{*}}$ is given by (10).

Having discussed the norms on the data space, in the next section, we study the norms on the solution space, which can be either componentwise or normwise. However, regardless of the norms chosen in the solution space, we always use the componentwise norm in the data space.

3 Condition Numbers for the Weighted Least Squares

Problem

In this section, we present an explicit expression of the condition number for the weighted least squares problem. First, we derive an explicit expression of the Fréchet derivative of the mapping $g$ in (7), when $A$ is of full column rank. Let $B\in\mathbb{R}^{m\times n}$, $\mathbf{c}\in\mathbb{R}^{m}$ and $J=g^{\prime}(A,\mathbf{b})$ be the derivative, applying the chain rule, we get

[TABLE]

recalling that $\mathbf{d}=W(\mathbf{b}-A\mathbf{x})$. Note that $J(B,\mathbf{c})$ is a mapping from the data space $\mathbb{R}^{m\times n}\times\mathbb{R}^{m}$ to $\mathbb{R}^{k}$.

From the definition of the adjoint operator and the definition of the scalar product in the data space $\mathbb{R}^{m\times n}\times\mathbb{R}^{m}$, the following lemma gives an explicit expression of the adjoint operator of the above $J(B,\mathbf{c})$.

Lemma 4

The adjoint operator of the Fréchet derivative $J(B,\mathbf{c})$ in (11) is

[TABLE]

for $\mathbf{u}\in\mathbb{R}^{k}$. Note that $J^{*}(\mathbf{u})$ is a mapping from $\mathbb{R}^{k}$ to $\mathbb{R}^{m\times n}\times\mathbb{R}^{m}$.

Proof. Let $J_{1}(B)$ and $J_{2}(\mathbf{c})$ be the first and second terms in the sum (11) respectively. By the definition of the scalar product in the matrix space, for any $\mathbf{u}\in\mathbb{R}^{k}$, we have

[TABLE]

For $J_{2}(\mathbf{c})$, we have

[TABLE]

Let

[TABLE]

and

[TABLE]

then $\langle J^{*}(\mathbf{u}),\ (B,\mathbf{c})\rangle=\langle(J_{1}^{*}(\mathbf{u}),J_{2}^{*}(\mathbf{u})),\ (B,\mathbf{c})\rangle=\langle\mathbf{u},\ J(B,\mathbf{c})\rangle$, which completes the proof. $\Box$

Having obtained an explicit expression of the adjoint operator of the Fréchet derivative, next we give an explicit expression of the condition number $\mathcal{K}$ (8) in terms the dual norm in the solution space in the following theorem, where $\mathrm{vec}(A)$ denotes the vector obtained by stacking the columns of a matrix $A$, $D_{A}$ denotes the diagonal matrix $\mathrm{diag}(\mathrm{vec}(A))$, and $\otimes$ is the Kronecker product operator [13].

Theorem 5

The condition number for the full rank weighted least squares problem can be expressed by

[TABLE]

where

[TABLE]

Proof. Let $\Delta A=[\Delta a_{i,j}]$ and $\Delta\mathbf{b}=[\Delta b_{i}]$, then, from (10), we have

[TABLE]

Applying Lemma 4, we get

[TABLE]

It can be verified that $d_{i}(A^{\mathrm{T}}WA)^{-1}\mathbf{e}_{j}$ is the $(i+(j-1)\,n)$th column of the $n\times(mn)$ matrix $(A^{\mathrm{T}}WA)^{-1}\otimes\mathbf{d}^{\mathrm{T}}$ and $x_{j}A_{W,I_{n}}^{\dagger}\mathbf{e}_{i}$ is the $(i+(j-1)\,n)$th column of the $n\times(mn)$ matrix $\mathbf{x}^{\mathrm{T}}\otimes A_{W,I_{n}}^{\dagger}$ in $V$ (12), implying that the above expression equals

[TABLE]

The theorem then follows from Lemma 3. $\Box$

The following case study discusses some commonly used norms for the norm in the solution space to obtain some specific expressions of the condition number $\mathcal{K}$.

Corollary 6

Using the above notations, when the infinity norm is chosen as the norm in the solution space $G$, we get

[TABLE]

Proof. When $\|\cdot\|_{G}=\|\cdot\|_{\infty}$, the dual norm $\|\cdot\|_{G^{*}}=\|\cdot\|_{1}$. Thus

[TABLE]

The following corollary gives an alternative expression of $\mathcal{K}_{\infty}$.

Corollary 7

Using the above notations, when the infinity norm is chosen as the norm in the solution space $G$, we get

[TABLE]

Proof. Partitioning

[TABLE]

where each $V_{j}$, $1\leq j\leq n$, is an $n\times m$ matrix, we get

[TABLE]

Recalling that $d_{i}(A^{\mathrm{T}}WA)^{-1}\mathbf{e}_{j}-x_{j}A_{W,I_{n}}^{\dagger}\mathbf{e}_{i}$ is the $(i+(j-1)n)$th column of $V$, we have

[TABLE]

The expression (14) is obtained by substituting $V_{j}$ in (15) with the above expression for $V_{j}$ and noticing that $A_{W,I_{n}}^{\dagger}=(A^{\mathrm{T}}WA)^{-1}A^{\mathrm{T}}W$. $\Box$

The advantage of the expression (14) over (13) is the absence of the Kronecker product. Consequently, its computation requires less memory. To further reduce the computational cost, we will propose efficient methods for estimating an upper bound of $\mathcal{K}_{\infty}$ in Section 4.

When $W=I_{m}$, the weighted least squares problem (1) reduces to the standard least squares problem (2) and the condition number $\mathcal{K}_{\infty}$ in (14) reduces to the condition number

[TABLE]

for the standard least squares problem given by Baboulin and Gratton [3, (3.4)], noticing that when $W=I_{m}$, we have $A_{W,I_{n}}^{\dagger}=A^{\dagger}$, $A^{\mathrm{T}}WA=A^{\mathrm{T}}A$, and $\mathbf{d}=W(\mathbf{b}-A\mathbf{x})=\mathbf{b}-AA^{\dagger}\mathbf{b}=\mathbf{r}$.

Next, we consider the 2-norm and derive an upper bound.

Corollary 8

When the 2-norm is used in the solution space, we have

[TABLE]

Proof. When $\|\cdot\|_{G}=\|\cdot\|_{2}$, then $\|\cdot\|_{G^{*}}=\|\cdot\|_{2}$. From Theorem 5,

[TABLE]

It follows from [12] that for any matrix $B$, $\|B\|_{2,1}=\max_{\|\mathbf{u}\|_{2}=1}\|B\mathbf{u}\|_{1}=\|B\hat{\mathbf{u}}\|_{1}$, where $\hat{\mathbf{u}}\in\mathbb{R}^{k}$ is a unit 2-norm vector. Applying $\|\hat{\mathbf{u}}\|_{1}\leq\sqrt{k}\,\|\hat{\mathbf{u}}\|_{2}$, we get

[TABLE]

Substituting the above $B$ with $[VD_{A}\ \ A_{W,I_{n}}^{\dagger}D_{\mathbf{b}}]^{\mathrm{T}}L$, we have

[TABLE]

which implies (16). $\Box$

The above upper bound for $\mathcal{K}_{2}$ can be obtained by computing (13) or (14).

So far, we have discussed the various mixed condition numbers, that is, componentwise norm in the data space and the infinite norm or 2-norm in the solution space. In the rest of the section, we study the case of componentwise condition number, that is, componentwise norm in the solution space as well.

Corollary 9

Considering the componentwise norm defined by

[TABLE]

in the solution space, we have the following three expressions for the componentwise condition number

[TABLE]

Proof. The expressions immediately follow from Theorem 5 and Corollaries 6 and 7. $\Box$

Similarly to $\mathcal{K}_{\infty}$, when $W=I_{m}$, the above condition number $\mathcal{K}_{c}$ reduces to the standard least squares condition number

[TABLE]

presented in [3, page 15].

4 Condition Number Estimators

In this section, we propose efficient methods for estimating $\mathcal{K}_{\infty}$ and $\mathcal{K}_{c}$, when integrated into the Paige’s method [23, 24] for solving the weighted least squares problem.

Firstly, we give upper bounds for $\mathcal{K}_{\infty}$ and $\mathcal{K}_{c}$ in the following theorem.

Theorem 10

Using the notations above, we have the upper bounds

[TABLE]

and

[TABLE]

Proof. From the monotonicity property of infinity norm and triangle inequality, we get

[TABLE]

where the last equation can be obtained by applying

[TABLE]

where ${\bf 1}=[1,...,1]^{\rm T}$.

The upper bound of $\mathcal{K}_{c}$ can be derived similarly. $\Box$

Our experiments show that the above upper bounds are tight.

The above upper bounds can be computed efficiently when the weighted least squares problem is solved by the fast numerically stable method proposed by Paige [23, 24]. To see this, we briefly describe the Paige’s method.

The weighted least squares problem (1) arises in finding the least squares estimate of the vector $\mathbf{x}$ in the linear model $\mathbf{b}=A\mathbf{x}+\mathbf{w}$, where $A$ is an $m\times n$ matrix and $\mathbf{w}$ is an unknown noise vector of zero mean and $m\times m$ covariance $Z=W^{-1}$. Usually, the factorization $Z=BB^{\rm T}$ is available. Paige considers the following form equivalent to the weighted least squares (1):

[TABLE]

By applying the plane rotations, we can get a generalized QR factorization [23, (2.1)] of the data matrix $[\mathbf{b}\ \ A\ \ B]$ of (18):

[TABLE]

where $Q,P\in{\mathbb{R}}^{m\times m}$ are orthogonal matrices and $L_{1}$, $L_{2}$, $R^{\rm T}$ are lower triangular and nonsingular, and $\rho$ is nonzero, assuming $A$ is of full column rank and $B$ is nonsingular and lower triangular. It is shown that the weighted least squares solution $\mathbf{x}$ can be obtained by solving the following nonsingular lower triangular system:

[TABLE]

The cost of Paige’s algorithm is $O(m^{2}n/2+mn^{2}-2n^{3}/3)$ [23, (4.4)].

Letting

[TABLE]

then

[TABLE]

From (19), we have

[TABLE]

Consequently, we get the factorizations:

[TABLE]

Note that $L_{1}$, $L_{2}$, $R$, and $S$ are triangular matrices. Then, using the above factorizations, the three norms in $\mathcal{K}_{\infty}^{u}$ or $\mathcal{K}_{c}^{u}$ can be efficiently estimated by the classical condition estimation method [18, Chapter 15], as shown in Algorithm 1.

Algorithm 1

Estimating $\mathcal{K}_{\infty}^{u}$.

*  Initial vectors $\mathbf{h}_{i}=k^{-1}{\bf 1}\in{\mathbb{R}}^{k}$, $i=1,2,3$; *

*  for $p=1,2,\ldots$ do

*     Using (20), calculate *

[TABLE]

*     Compute $\mathbf{s}_{i}={\rm sign}(\mathbf{y}_{i})$, $i=1,2,3$, where sign is the sign function; *

*     Using (20), calculate *

[TABLE]

*     if $\|\mathbf{z}_{i}\|_{\infty}\leq\mathbf{h}_{i}^{\rm T}\mathbf{z}_{i}$ then

*        $\gamma_{i}=\left\|\mathbf{y}_{i}\right\|_{1}$, $i=1,2,3$; *

*        break *

*     end *if

*     $\mathbf{h}_{i}=\mathbf{e}_{k_{i}}$, where $k_{i}$ is the smallest index such that $|z_{k_{i}}|=\|\mathbf{z}_{i}\|_{\infty}$; *

*  end *for

*  Return $\hat{\mathcal{K}}_{\infty}^{u}=\gamma_{1}+\gamma_{2}+\gamma_{3}$. *

Table 1 lists the major costs in Algorithm 1 when integrated into the Paige’s method, where $\mathbf{v}$ is a vector with conformal dimensions. Let $p_{\max}$ be the total number of iterations, then the total cost of Algorithm 1 is $O(p_{\max}(m^{2}+mn+n^{2}))$. Recalling that the cost of the Paige’s method for solving the weighted least squares problem is $O(m^{2}n/2+mn^{2}-2n^{3}/3)$.

5 An Error Analysis of the Augmented System

In this section we perform a componentwise perturbation analysis of the augmented system (4) for the weighted least squares problem. Our analysis is a generalization of the analysis of the standard least square problem by Arioli et al. [1] and Björck [4].

Let the perturbations $\Delta A\in\mathbb{R}^{m\times n}$ and $\mathbf{b}\in\mathbb{R}^{m}$ satisfy $|\Delta A|\leq\epsilon\,|A|$ and $|\mathbf{b}|\leq\epsilon\,|\mathbf{b}|$ for a small $\epsilon$. Suppose that the perturbed augmented system is

[TABLE]

Denoting

[TABLE]

and the perturbations

[TABLE]

When $A$ is of full column rank, $G$ is invertible. It can be verified that

[TABLE]

We know that if the spectral radius

[TABLE]

then $I_{m+n}+G^{-1}\Delta G$ is invertible. Clearly, the condition

[TABLE]

implies (21). The following results [27] are necessary for Theorem 12.

Lemma 11

For a linear system $G\mathbf{z}=\mathbf{f}$ and its perturbed system

[TABLE]

where $\mathbf{z}+\Delta\mathbf{z}$ is the solution to the perturbed system, when the perturbations $\Delta G$ and $\Delta\mathbf{f}$ are sufficiently small such that $G+\Delta G$ is invertible, the perturbation $\Delta\mathbf{z}$ in the solution $\mathbf{z}$ satisfies

[TABLE]

which implies

[TABLE]

Furthermore, when the spectral radius $\rho(|G^{-1}|\,|\Delta G|)<1$, we have

[TABLE]

Now we have the bounds for the perturbations in the weighted least squares solution and residual.

Theorem 12

Using the above notations, for any $\epsilon>0$ satisfying the condition (22), when the componentwise perturbations $|\Delta A|\leq\epsilon\,|A|$ and $|\Delta\mathbf{b}|\leq\epsilon\,|\mathbf{b}|$, the error in the solution is bounded by

[TABLE]

and error in the weighted residual is bounded by

[TABLE]

Proof. Since the condition (22) implies (21), applying (23) in Lemma 11, we get

[TABLE]

Finally, by using the conditions $|\Delta A|\leq\epsilon\,|A|$ and $|\Delta\mathbf{b}|\leq\epsilon\,|\mathbf{b}|$ and the explicit form of $G^{-1}$, the upper bounds (24) and (25) can be obtained. $\Box$

6 Numerical Experiments

In this section, we present our experimental results to demonstrate the effectiveness of our condition numbers and their estimators for the weighted LS problem. All the numerical experiments were carried out in Matlab 2015b, with the machine precision $\mu\approx 2.2\times 10^{-16}$.

Firstly, we adopted the example in Baboulin and Gratton [3] and modified $A$, $W$, and $\mathbf{b}$ as the following:

[TABLE]

where $\epsilon,\gamma>0$ and $U$ is a random orthogonal matrix obtained from the QR decomposition of a random matrix. As we can see, $A$ (or $W$) becomes ill-conditioned as $\epsilon$ (or $\gamma$) decreases to zero. The vector $\mathbf{b}$ is constructed so that the solution is imbalanced, that is, its components range widely, to show the benefit of the componentwise condition. Note that ${\bf b}_{1}\in{\rm Rang}(A)$ and ${\bf b}_{2}\in{\rm Ker}(A^{\rm T})$, where ${\rm Rang}(A)$ and ${\rm Ker}(A^{\rm T})$ denote the range space of $A$ and the null space of $A^{\rm T}$, respectively. We generated the perturbations:

[TABLE]

where entries of $E$ and ${\bf f}$ are random variables uniformly distributed in the interval $(-1,\,1)$. Thus the perturbation size $\|(\Delta A,\Delta\mathbf{b})\|_{c}\approx 10^{-8}$. For the $L$ matrix in our condition numbers, we chose

[TABLE]

That is, corresponding to the above three matrices, the whole $\mathbf{x}$, the subvector $[x_{1}\ x_{2}]^{\mathrm{T}}$, and the component $x_{3}$ are selected respectively. We called the Matlab built-in function $\tt lscov$ to compute the solutions ${\bf x}$ and $\tilde{\mathbf{x}}$ corresponding to the unperturbed WLS (1) and its perturbed WLS defined by

[TABLE]

We measured the mixed and componentwise relative errors in $L^{\rm T}\mathbf{x}$ defined by

[TABLE]

where $\|\cdot\|_{c}$ is the componentwise norm defined in (17). Since the data perturbation size is about $10^{-8}$, $\mathcal{E}_{\infty}^{\rm rel}\times 10^{8}$ and $\mathcal{E}_{c}^{\rm rel}\times 10^{8}$ are respectively indications of the mixed and componentwise condition numbers for this particular problem. Specifically, in the table, our condition numbers are

[TABLE]

where $V$ is defined in (12). Also, in the table, we define the upper bound $\mathcal{K}_{\infty}^{u,{\rm rel}}=\mathcal{K}_{\infty}^{u}/\|L^{\rm T}\mathbf{x}\|_{2}$, recalling that $\mathcal{K}_{\infty}^{u}$ and $\mathcal{K}_{c}^{u}$ are defined in Theorem 10.

Table 2 compares our condition numbers $\mathcal{K}_{\infty}^{\rm rel}$ and $\mathcal{K}_{c}^{\rm rel}$ with their corresponding relative errors $\mathcal{E}_{\infty}^{\rm rel}$ and $\mathcal{E}_{c}^{\rm rel}$. First, the table shows that our condition numbers, mixed and componentwise, are consistently close to the estimates $\mathcal{E}^{\rm rel}\times 10^{8}$. Second, our componentwise condition numbers show that the third component of the solution is better conditioned than the first two, showing the benefit of the componentwise analysis. Third, in the case when $\epsilon=\delta=10^{-6}$, our condition numbers are much larger than their corresponding estimates $\mathcal{E}_{\infty}^{\rm rel}\times 10^{8}$ and $\mathcal{E}_{c}^{\rm rel}\times 10^{8}$. Our explanation is that $\mathcal{E}_{\infty}^{\rm rel}\times 10^{8}$ and $\mathcal{E}_{c}^{\rm rel}\times 10^{8}$ give estimates of the condition numbers for this particular problem with this particular perturbation, whereas our condition numbers are upper bounds for this problem with general perturbation.

Secondly, we experimented on the linear model:

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{n}$ whose entries are random variables with standard normal distribution and $\mathbf{w}\in\mathbb{R}^{m}$ whose entries $w_{i}$ are random variables with normal distribution, mean 0, and predefined variances $\sigma_{i}^{2}$. Thus the weight matrix $W=D_{\bf z}^{-1}$, where ${\bf z}=[\sigma_{i}^{2}]$. In our experiments, we set $m=50$ and $n=10$. The $m\times n$ matrix $A$ was generated by the Matlab built-in function $\mathtt{sprandn}$ with density $0.5$. The same as before, the perturbations on $A$ and $\mathbf{b}$ were generated by (28) and both the unperturbed and perturbed weighted least squares problems were solved by the Matlab function $\mathtt{lscov}$.

For the $L$ matrix in our condition numbers, we chose

[TABLE]

Thus, corresponding to the above three matrices, the whole $\mathbf{x}$, the subvector $[x_{1}\ x_{2}]^{\mathrm{T}}$, and the last component $x_{n}$ are selected respectively.

To investigate the impact of the variances $\sigma_{i}^{2}$, we first set $\mathbf{z}=[10^{-4}:(4\times 10^{-4})/(m-1):5\times 10^{-4}]$, that is, $\sigma_{i}^{2}$ are evenly spaced between $10^{-4}$ and $5\times 10^{-4}$. We generated 1,000 samples of $A$ and $\mathbf{b}$ each. The mean values of $\mathcal{E}_{\infty}^{\rm rel}$, $\mathcal{E}_{c}^{\rm rel}$, $\mathcal{K}_{\infty}^{\rm rel}$ and $\mathcal{K}_{c}^{\rm rel}$, $\mathcal{K}_{\infty}^{u,{\rm rel}}$ and $\mathcal{K}_{c}^{u}$ are displayed in Table 3. As expected, the condition numbers are moderate, when all the variances are small, that is, the weight matrix $W$ is well-conditioned, and the data matrix $A$ is well conditioned since it is a random matrix.

We then widened the range of the variances. Specifically, $\sigma_{i}^{2}$ are evenly spaced between $10^{-4}$ and $10^{2}$. Table 4 shows the average values of $\mathcal{E}_{\infty}^{\rm rel}$, $\mathcal{E}_{c}^{\rm rel}$, $\mathcal{K}_{\infty}^{\rm rel}$ and $\mathcal{K}_{c}^{\rm rel}$, $\mathcal{K}_{\infty}^{u,{\rm rel}}$ and $\mathcal{K}_{c}^{u}$ over 1,000 samples of $A$ and $\mathbf{b}$ each. As the range of the variances is widened, that is, the condition number of the weight matrix $W$ increases, the condition numbers of the weighted least squares problem increase. However, both Tables 3 and 4, along with Table 2, show that the condition of the weighted least squares problem is more sensitive to the condition of the data matrix $A$ than that of the weight matrix $W$.

Our experiments show that our upper bounds are consistently very close to their corresponding condition numbers. In other words, our condition number estimators are accurate as well as efficient.

Conclusion

By applying adjoint operator and dual norm theory, we define the mixed and componentwise condition numbers for the linear solution function of the weighted linear least squares problem. Both the normwise and componentwise perturbation analyses of the solution are performed. Moreover, we present the componentwise perturbation analysis of both the solution and the residual of the augmented system of the weighted least squares problem. We also propose two efficient condition number estimators. Our numerical experiments show that our condition numbers are tight and can reveal the condition numbers of individual components of the solution. Moreover, our condition number estimators are accurate as well as efficient.