Weighted Least Squares Perturbation Theory

Abstract

The interest in the problem of weighted pseudoinverse matrices and the problem of weighted least squares (WLS) is largely due to their numerous applications. In particular, the problem of WLS is used in the design and optimization of building structures, in tomography, in statistics, etc. The first part of the chapter is devoted to the sensitivity of the solution to the WLS problem with approximate initial data. The second part investigates the properties of a SLAE with approximate initial data and presents an algorithm for finding a weighted normal pseudo solution of a WLS problem with approximate initial data, an algorithm for solving a WLS problem with symmetric positive semidefinite matrices and an approximate right side and also a parallel algorithm for solving a WLS problem. The third part is devoted to the analysis of the reliability of computer solutions of the WLS problem with approximate initial data. Here, estimates of the total error of the WLS problem are presented, and also software-algorithmic approaches to improving the accuracy of computer solutions.

1. Introduction

The interest in the problem of weighted pseudoinverse matrices and the WLS problem is largely due to their numerous applications. In particular, the problem of weighted least squares is used in the design and optimization of building structures, in tomography, in statistics, etc. A number of properties of weighted pseudoinverse matrices underlie the finding of weighted normal pseudosolutions. The field of application of weighted pseudoinverse matrices and weighted normal pseudosolutions is constantly expanding.

The definition of a weighted pseudoinverse matrix with positive definite weights was first introduced by Chipman in article [1]. In 1968, Milne introduced the definition of a skew pseudoinverse matrix in paper [2]. The study of the properties of weighted pseudoinverse matrices and weighted normal pseudosolutions, as well as the construction of methods for solving these and other problems, are devoted to the works of Mitra, Rao, Van Loan, Wang, Galba, Deineka, Sergienko, Ben-Israel, Elden, Wei, Wei, Ward etc. Weighted pseudoinverse matrices and weighted normal pseudosolutions with degenerate weights were studied in [3, 4, 5]. The existence and uniqueness of weighted pseudoinverse matrices with indefinite and mixed weights, as well as some of their properties, were described in [6, 7, 8]. Application of the weighted pseudoinverse matrix in statistics presented, for example, in [9, 10]. Many results on weighted generalized pseudo-inversions can be found in monographs [11, 12]. Much less work is devoted to the study of weighted pseudoinversion under conditions of approximate initial data. These issues are discussed in [13, 14, 15, 16, 17]. Analysis of the properties of weighted pseudoinverses and weighted normal pseudosolutions, as well as the creation of solution methods for these and other problems, are described in [18, 19, 20].

When solving applied problems, their mathematical models will have, as a rule, approximate initial data as a result of measurements, observations, assumptions, hypotheses, etc. Later, during discretization (‘arithmetization’) of the mathematical model, these errors are transformed into the errors of the matrix elements and the right parts of the resolving systems of equations. The input data of systems of linear algebraic equations and WLS problems can be determined directly from physical observations, and therefore they can have errors inherent in all measurements. In this case, the original data we have is an approximation of some exact data. And, finally, the initial data of mathematical models formulated in the form of linear algebra problems can be specified exactly in the form of numbers or mathematical formulas, but, given the finite length of a machine word, it is impossible to work with such an exact model on a computer. The machine model of such a problem in the general case will be approximate either due to errors in converting numbers from the decimal system to binary or due to rounding errors in the implementation of calculations on a computer.

The task is to study the properties of the machine model and to form a model of the problem and an algorithm for obtaining an approximate solution in a machine environment that will approximate the solution of a mathematical problem. The key question of numerical simulation is the reliability of the obtained machine solutions.

The most complete systematic exposition of questions related to the approximate nature of the initial data in problems of linear algebra is given in the monographs [21, 22, 23, 24]. Various approaches to the study and solution of ill-posed problems were considered, for example, in [25, 26, 27, 28]. Problems of the reliability of a machine solution for problems with approximate initial data, i.e. estimates of the proximity of the machine solution to the mathematical solution, estimates of the hereditary error in the mathematical solution and refinement of the solution were considered in the publications [12, 26, 29, 30, 31, 32, 33]. Much less work has been devoted to the study of similar questions for the WLS problem. The sensitivity analysis of a weighted normal pseudosolution under perturbation of the matrix and the right-hand side is the subject of papers [16, 34, 35, 36].

The chapter is devoted to the solution of the listed topical problems, namely the development of the perturbation theory for the WLS problem with positive definite weights and the development of numerical methods for the study and solution of mathematical models with approximate initial data.

2. Weighted least squares problem

2.1 Preliminaries

Let the set of all $m\times n$ matrices is denoted by $R^{m\times n}$. Given a matrix $A\in R^{m\times n}$ let $A^T$ is the transpose of A, $rank\left(A\right)$ is the rank of A, $\mathbb{R}\left(A\right)$ is the field of values of A and $\mathbb{N}\left(A\right)$ is the null space of A. Additionally, let $‖‖$ denote the vector 2-norm and the consistent matrix 2-norm, and let I be an identity matrix.

Given an arbitrary matrix $A\in R^{m\times n}$ and symmetric positive definite matrices M and N of orders m and n, respectively, a unique matrix $X\in R^{m\times n}$, satisfying the conditions:

$$\mathit{AXA}=A,\mathit{XAX}=X,{\left(\mathit{MAX}\right)}^T=\mathit{MAX},{\left(\mathit{NXA}\right)}^T=\mathit{NXA},$$E1

is called the weighted Moore–Penrose pseudoinverse of A and is denoted by $X=A_{\mathit{MN}}^{+}$. Specifically, if $M=I\in R^{m\times m}$ and $N=I\in R^{n\times n}$, then X satisfying conditions (1) is called the Moore–Penrose pseudoinverse and is designated as $X=A^{+}$.

Let $A^{\#}$ denote the weighted transpose of A, P and Q be idempotent matrices, and $\overline{А}=A+\mathrm{\Delta }A$ be a perturbed matrix, i.e.,

$$A^{\#}=N^{-1}{A}^{T}M.$$E2

$$P=A_{\mathit{MN}}^{+}A,Q=A{A}_{\mathit{MN}}^{+},\overline{P}={\overline{A}}_{\mathit{MN}}^{+}\overline{A},\overline{Q}=\overline{A}{\overline{A}}_{\mathit{MN}}^{+}.$$E3

Let $x\in R^m$, $y\in R^n$. The weighted scalar products in $R^m$ and $R^n$ are defined as ${\left(x,y\right)}_M=y^T\mathit{Mx}$, $x,y\in R^m$ and ${\left(x,y\right)}_N=y^T\mathit{Nx}$, $x,y\in R^n$, respectively. The weighted vector norms are defined as:

$$\begin{array}{c}{‖x‖}_M={\left(x,x\right)}_M^{\frac{1}{2}}={\left(x^T\mathit{Mx}\right)}^{\frac{1}{2}}=‖M^{\frac{1}{2}}x‖,x\in R^m,\\ {‖y‖}_N={\left(y,y\right)}_N^{\frac{1}{2}}={\left(y^T\mathit{Ny}\right)}^{\frac{1}{2}}=‖N^{\frac{1}{2}}y‖,y\in R^n.\end{array}$$E4

Let $x,y\in R^m$ and ${\left(x,y\right)}_M=0$. Then the vectors $x$ and $y$ are called M-orthogonal, i.е. $M^{\frac{1}{2}}x$- and $M^{\frac{1}{2}}y$-orthogonal. It is easy to show that.

$${‖x+y‖}_M^2={‖x‖}_M^2+{‖y‖}_M^2,x,y\in R^m.$$E5

The weighted matrix norms are defined as:

$$\begin{array}{l}{‖A‖}_{\mathit{MN}}=\underset{{‖x‖}_N=1}{max}{‖\mathit{Ax}‖}_M=‖M^{\frac{1}{2}}A{N}^{-\frac{1}{2}}‖,A\in R^{m\times n},\\ {‖B‖}_{\mathit{NM}}=\underset{{‖y‖}_M=1}{max}{‖\mathit{By}‖}_N=‖N^{\frac{1}{2}}A{M}^{-\frac{1}{2}}‖,B\in R^{n\times m}.\end{array}$$E6

Lemma 1 (see in [37]). Let $A\in R^{m\times n}$, $rank\left(A\right)=k$, M and N are positive definite matrices of orders m and n, respectively. Then, there are matrices $U\in R^{m\times m}$ and $V\in R^{n\times n}$, satisfying $U^T\mathit{MU}=I$ and $V^T{N}^{-1}V=I$ such that.

$$A=U\left(\begin{array}{cc}D& 0\\ 0& 0\end{array}\right)V^T,A_{\mathit{MN}}^{+}=N^{-1}V\left(\begin{array}{cc}{D}^{-1}& 0\\ 0& 0\end{array}\right)U^{T}M,$$E7

where $D=diag\left({\mu }_1,{\mu }_2,\dots ,{\mu }_k\right)$, ${\mu }_1\ge {\mu }_2\ge \dots \ge {\mu }_k>0$ and ${\mu }_i^2$ are the nonzero eigenvalues of the matrix $A^{\#}A$. The nonnegative values ${\mu }_i$ are called the weighted singular values of A, moreover, ${‖A‖}_{\mathit{MN}}={\mu }_1$, ${‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}=\frac{1}{{\mu }_k}$.

The weighted singular value decomposition of A yields an M-orthonormal basis of the vectors of U and an $N^{-1}$-orthonormal basis of the vectors of V.

2.2 Statement of the problem

In the study of the reliability of the obtained machine results, three linear systems are considered. A system of linear algebraic equations with exact input data

$$\mathit{Ax}=b.$$E8

We will consider the corresponding weighted least squares problem with positive definite weights M and N:

$$\underset{x\in C}{min}{‖x‖}_N,C=\left\{x|{‖\mathit{Ax}-b‖}_M=min\right\},$$E9

where $A\in R^{m\times n}$ is a rank-deficient matrix and $b\in R^m$.

Along with (9), we consider the mathematical model with approximately specified initial data.

$$\underset{x\in C}{min}{‖\overline{x}‖}_N,C=\left\{\overline{x}|{‖\left(A+\mathrm{\Delta }A\right)\overline{x}-\left(b+\mathrm{\Delta }b\right)‖}_M=min\right\},$$E10

where.

$$\overline{А}=A+\mathrm{\Delta }A,\overline{b}=b+\mathrm{\Delta }b,\overline{x}=x+\mathrm{\Delta }x.$$E11

Assume that the errors in the matrix elements and the right-hand side satisfy the relations:

$${‖\mathrm{\Delta }A‖}_{\mathit{MN}}\le {\epsilon }_A{‖A‖}_{\mathit{MN}},{‖\mathrm{\Delta }b‖}_M\le {\epsilon }_b{‖b‖}_M.$$E12

The problem for the approximate solution $\overline{\overline{x}}$ of a system of linear algebraic equations with approximately given initial data

$$\overline{A}\overline{\overline{x}}=\overline{b}+\overline{r},$$E13

where $\overline{r}=\overline{A}\overline{\overline{x}}-\overline{b}$ is the residual vector.

The analysis of the reliability of the obtained solution includes an assessment of the hereditary error ${‖x-\overline{x}‖}_N$, computational error ${‖\overline{x}-\overline{\overline{x}}‖}_N$ and total error ${‖x-\overline{\overline{x}}‖}_N$, as well as the refinement of the obtained machine solution to a given accuracy.

2.3 The existence and uniqueness of a weighted normal pseudoinverse

Let linear manifold L be a nonempty subset of space R, closed with respect to the operations of addition and multiplication by a scalar (if $x$ and y are elements of L $\forall \alpha ,\beta$, the $\mathit{\alpha x}+\mathit{\beta y}$ is an element of L). Vector $x$ is N-orthogonal to the linear manifold L ($x{\perp }_{N}L$) if $x$ is N-orthogonal to each vector from L.

Lemma 2 (see in [38]). There exists a unique decomposition of vector x, namely $x=\widehat{x}+\stackrel{\sim }{x}$, where $\widehat{x}\in L$, $\stackrel{\sim }{x}{\perp }_{N}L$.

Let A is an arbitrary matrix. The kernel of matrix A, denoted by $\mathbb{N}$ (A), is the set of vectors mapped into zero by А: $\mathbb{N}\left(\mathrm{A}\right)=\left\{x:\mathit{Ax}=0\right\}$.

The set $\mathbb{R}\left(A\right)$ of images of matrix A is the set of vectors that are images of vectors of the space R from the definition domain of A, i.e. $\mathbb{R}\left(A\right)=\left\{b:b=\mathit{Ax},\forall x\right\}$ .

Let L be a linear manifold in space R, N-orthogonal (M-orthogonal) complement to L, denoted by $L^{{\perp }_N}$ ($L^{{\perp }_M}$), defined as the set of vectors in R, each of which is N-orthogonal (M-orthogonal) to L.

Remark 1. If $x$ is a vector from R and $x^T\mathit{Ny}=0$ for any y from R, then $x=0$.

Theorem 1. Let $A\in R^{m\times n}$, then $\mathbb{N}\left(A\right)={\mathbb{R}}^{\perp }\left(A^{\#}\right)$.

Proof. Vector $x\in \mathbb{N}\left(A\right)$, if and only if $\mathit{Ax}=0$. Hence, by virtue of Remark1, we get $x\in \mathbb{N}\left(A\right)$, if and only if $y^T\mathit{MAx}=0$ for any y. Since $y^T\mathit{MAx}={\left(A^{\#}y\right)}^T\mathit{Nx}$, we get $\mathit{Ax}=0$ if and only if $x$ is N-orthogonal to all the vectors of the form$A^{\#}y$. Vectors$A^{\#}y$ form $\mathbb{R}\left(A^{\#}\right)$. The required statement follows from here and from the definition ${\mathbb{R}}^{{\perp }_M}\left(A\right)$.

Theorem 2 (see in [38]). If A is an $m\times n$ matrix and b is an m-dimensional vector, then the unique decomposition $b=\widehat{b}+\stackrel{\sim }{b}$ holds, where $\widehat{b}\in \mathbb{R}\left(A\right)$ and $\stackrel{\sim }{b}\in \mathbb{N}\left(A^{\#}\right)$.

Vector$\widehat{b}$ is a projection of b on to $\mathbb{R}\left(A\right)$, and $\stackrel{\sim }{b}$ is a projection of b on to $\mathbb{N}\left(A^{\#}\right)$. Vectors$\widehat{b}$ and $\stackrel{\sim }{b}$ are M-orthogonal. Hence, $A^{\#}b=A^{\#}\widehat{b}$.

By Theorem 1, the following relations hold for the symmetric matrix A: $\mathbb{N}\left(A\right)={\mathbb{R}}^{\perp }\left(A\right),\mathbb{R}\left(A\right)={\mathbb{N}}^{\perp }\left(A\right).$

Theorem 3. Let $A\in R^{m\times n}$, then $\mathbb{R}\left(A\right)=\mathbb{R}\left(A{A}^{\#}\right)$, $\mathbb{R}\left(A^{\#}\right)=\mathbb{R}\left(A^{\#}A\right)$, $\mathbb{N}\left(A\right)=\mathbb{N}\left(A^{\#}A\right)$ and $\mathbb{N}\left(A^{\#}\right)=\mathbb{N}\left(A{A}^{\#}\right)$ .

Proof. It will be to establish that $\mathbb{N}\left(A^{\#}\right)=\mathbb{N}\left(A{A}^{\#}\right)$ and $\mathbb{N}\left(A\right)=\mathbb{N}\left(A^{\#}A\right)$.

For this purpose, we will use Theorem 1. To prove the coincidence of $\mathbb{N}\left(A^{\#}\right)$ and $\mathbb{N}\left(A{A}^{\#}\right)$, note that $A{A}^{\#}x=0$ if $A^{\#}x=0$. On the other hand, if $A{A}^{\#}x=0$, then $x^{T}A{A}^{\#}x=0$, i.е. ${‖A^{\#}x‖}_M=0$, which entails equality $A^{\#}x=0$. So, $A^{\#}x=0$ if and only if $x^{T}A{A}^{\#}x=0$. We can similarly establish that$\mathbb{N}\left(A\right)=\mathbb{N}\left(A^{\#}A\right)$.

Then let us prove the theorem about the existence and uniqueness of the solution vector that minimizes the norm of the residual ${‖\mathit{Ax}-b‖}_M$ by the technique proposed in [39] for the least-squares problem.

Theorem 4. Let $A\in R^{m\times n}$, $b\in R^m$, $b\notin \mathbb{R}\left(A\right)$. Then there exists a vector $\widehat{x}$, that minimizes the norm of the residual ${‖\mathit{Ax}-b‖}_M$ and vector $\widehat{x}$ is a unique vector from $\mathbb{R}\left(A^{\#}\right)$, that satisfies the equation $\mathit{Ax}=\widehat{b}$, where $\widehat{b}=A{A}_{\mathit{MN}}^{+}b$ is the projection of b onto$\mathbb{R}\left(A\right)$.

Proof. By virtue of Theorem 2, we get $b=\widehat{b}+\stackrel{\sim }{b}$, where $\stackrel{\sim }{b}=\left(I-A{A}_{\mathit{MN}}^{+}\right)b$ is the projection of b on to $\mathbb{N}\left(A^{\#}\right)$. Since for every $x,\mathit{Ax}\in \mathbb{R}\left(A\right)$ and $\stackrel{\sim }{b}\in {\mathbb{R}}^{{\perp }_M}\left(A\right)$, then $\widehat{b}-\mathit{Ax}\in \mathbb{R}\left(A\right)$ and $\stackrel{\sim }{b}\perp \widehat{b}-\mathit{Ax}$. Therefore

$${‖b-\mathit{Ax}‖}_M^2={‖\widehat{b}-\mathit{Ax}+\stackrel{\sim }{b}‖}_M^2={‖\widehat{b}-\mathit{Ax}‖}_M^2+{‖\stackrel{\sim }{b}‖}_M^2\ge {‖\stackrel{\sim }{b}‖}_M^2.$$E14

This lower bound is attained since $\widehat{b}$ belongs to the set of images А, i.е. $\widehat{b}$ is an image of some $x_0$: $\widehat{b}=A{x}_0$.

Thereby, for this $x_0$ the greatest lower bound is attainable:

$${‖b-A{x}_0‖}_M^2={‖b-\widehat{b}‖}_M^2={‖\stackrel{\sim }{b}‖}_M^2.$$E15

It was shown earlier that

$${‖b-\mathit{Ax}‖}_M^2={‖\widehat{b}-\mathit{Ax}‖}_M^2+{‖\stackrel{\sim }{b}‖}_M^2$$E16

and hence, the lower bound can only be attained for$x^{\ast }$, for which $A{x}^{\ast }=\widehat{b}$. According to Theorem 2, each vector $x^{\ast }$, can be presented as a sum of two orthogonal vectors: $x^{\ast }={\widehat{x}}^{\ast }+{\stackrel{\sim }{x}}^{\ast }$, where ${\widehat{x}}^{\ast }\in \mathbb{R}\left(A^{\#}\right)$, ${\stackrel{\sim }{x}}^{\ast }\in \mathbb{N}\left(A\right)$.

Therefore, $A{x}^{\ast }=A{\widehat{x}}^{\ast }$ and hence, ${‖b-A{x}^{\ast }‖}_M^2={‖b-A{\widehat{x}}^{\ast }‖}_M^2$. Note that

$${‖x^{\ast }‖}_N^2={‖{\widehat{x}}^{\ast }‖}_N^2+{‖{\stackrel{\sim }{x}}^{\ast }‖}_N^2\ge {‖{\widehat{x}}^{\ast }‖}_N^2,$$E17

where strict inequality is possible when$x^{\ast }\ne {\widehat{x}}^{\ast }$ (i.е. if $x^{\ast }$ does not coincide with its projection on to $\mathbb{R}\left(A^{\#}\right)$).

It was shown above, that $x_0$ minimizes ${‖\mathit{Ax}-b‖}_M$, if and only if $A{x}_0=\widehat{b}$, and among the vectors that minimize ${‖\mathit{Ax}-b‖}_M$, each vector with the minimum norm should belong to the set of images $A^{\#}$. To establish the uniqueness of a minimum-norm vector, assume that $\widehat{x}$ and$x^{\ast }$ belong to $\mathbb{R}\left(A^{\#}\right)$ and that $A\widehat{x}=A{x}^{\ast }=\widehat{b}$. Then$x^{\ast }-\widehat{x}\in \mathbb{R}\left(A^{\#}\right)$, however$A\left(x^{\ast }-\widehat{x}\right)=0$, so $x^{\ast }-\widehat{x}\in \mathbb{N}\left(A\right)={\mathbb{R}}^{{\perp }_N}\left(A^{\#}\right)$.

As vector$x^{\ast }-\widehat{x}$ is N-orthogonal to itself${‖x^{\ast }-\widehat{x}‖}_N=0$, i.е. $x^{\ast }=\widehat{x}$.

Remark 2. There is another assertion that is equivalent to Theorem 4. There exists an n-dimensional vector y such that

$${‖b-A{A}^{\#}y‖}_M=i\underset{x}{\mathit{nf}}{‖b-\mathit{Ax}‖}_M.$$E18

If

$${‖b-A{x}_0‖}_M=\underset{x}{\mathit{inf}}{‖b-\mathit{Ax}‖}_M,$$E19

then ${‖x_0‖}_N\ge {‖A^{\#}y‖}_N$ with strict inequality for$x_0\ne {‖A^{\#}y‖}_N$.

Vector y satisfies the equation $A{A}^{\#}y=\widehat{b}$, here $\widehat{b}$ is the projection of b onto $\mathbb{R}\left(A\right)$.

Theorem 5. Among all the vectors $x$ that minimize the residual ${‖\mathit{Ax}-b‖}_M$, vector $\widehat{x}$, which has the minimum norm $\widehat{x}=min{‖x‖}_N$, is a unique vector of the form

$$\widehat{x}=N^{-1}{A}^T\mathit{My}=A^{\#}y,$$E20

satisfying the equation

$$A^{\#}\mathit{Ax}=A^{\#}b,$$E21

i.е. $\widehat{x}$ can be obtained by means of any vector $y_0$, that satisfies the equation $A^{\#}A{A}^{\#}y=A^{\#}b$ by the formula $\widehat{x}=A^{\#}{y}_0$.

Proof. According to the condition of Theorem 3 $\mathbb{R}\left(A^{\#}\right)=\mathbb{R}\left(A^{\#}A\right)$. Since vector $A^{\#}b$ belongs to the set of images $A^{\#}$, it should belong to the set of images$A^{\#}A$ and thus should be an image of some vector $x$ with respect to the transformation $A^{\#}A$. In other words, Eq. (21) (with respect to $x$) has at least one solution. If $x$ is a solution of Eq. (21), then $\widehat{x}$ is the projection of $x$ on to $\mathbb{R}\left(A^{\#}\right)$, since $\mathit{Ax}=A\widehat{x}$ according to Theorem 2. Since $\widehat{x}\in \mathbb{R}\left(A^{\#}\right)$, vector $\widehat{x}$ is an image of some vector y with respect to the transformation $A^{\#}$: $\widehat{x}=A^{\#}y$.

Thus, we have established that there exists at least one solution of Eq. (21) in the form of (20). To establish the uniqueness of this solution, we assume that ${\widehat{x}}_1=A^{\#}{y}_1$ and ${\widehat{x}}_2=A^{\#}{y}_2$ satisfy Eq. (21). Then $A^{\#}A\left(A^{\#}{y}_1-A^{\#}{y}_2\right)=0$, therefore, $A^{\#}\left(y_1-y_2\right)\in \mathbb{N}\left(A^{\#}A\right)=\mathbb{N}\left(A\right)$, where from the equality $A{A}^{\#}\left(y_1-y_2\right)=0$ follows.

Therefore $y_1-y_2\in \mathbb{N}\left(A{A}^{\#}\right)=\mathbb{N}\left(A^{\#}\right)$; hence,${\widehat{x}}_1=A^{\#}{y}_1=A^{\#}{y}_2={\widehat{x}}_2$.

Thus, there exists exactly one solution of Eq. (21) in the form (20). The proof of Theorem 5 will be completed if we can show that by virtue of Theorem 1 the solution found in the form (14) is also a solution of the equation $\mathit{Ax}=\widehat{b}$, where $\widehat{b}$ is a weighted projection of b on to $\mathbb{R}\left(A\right)$, i.е. $A^{\#}b=A^{\#}\widehat{b}$.

Theorem 4 establishes that there is a unique, from $\mathbb{R}\left(A^{\#}\right)$ solution of the equation

$$\mathit{Ax}=\widehat{b}.$$E22

Hence, this unique solution satisfies the equation $A^{\#}\mathit{Ax}=A^{\#}\widehat{b}$.

According to the equality $A^{\#}b=A^{\#}\widehat{b}$ the unique solution of Eq. (22) belonging to $\mathbb{R}\left(A^{\#}\right)$, should coincide with $\widehat{x}$ which is a unique solution of Eq. (21), which also belongs to $\mathbb{R}\left(A^{\#}\right)$. Finally, vector $\widehat{x}$, mentioned in the proof of Theorem 5 exactly coincides with the vector $\widehat{x}$ from Theorem 4. Using the representation of the Moore–Penrose weighted pseudoinverse from [38].

$${\mathit{A}}_{\mathit{MN}}^{+}={\mathit{A}}^{\#}{\left({\mathit{A}}^{\#}\mathit{A}{\mathit{A}}^{\#}\right)}^{+}{\mathit{A}}^{\#},$$E23

we can formulate the following theorem for problem (9) in a shorter form.

Theorem 6. Let $\mathit{A}\in {\mathit{R}}^{\mathit{m}\times \mathit{n}}$, then $\mathit{x}={\mathit{A}}_{\mathit{MN}}^{+}\mathit{b}$—is М-weighted least squares solution with the minimum N-norm of the system $\mathit{Ax}=\mathit{b}$.

Note that in [18] a slightly different mathematical apparatus was used to prove the existence and uniqueness of the M-weighted least squares solution with the minimum N-norm of the system $\mathit{Ax}=b$.

3. Error estimates for the weighted minimum-norm least squares solution

3.1 Estimates of the hereditary error of a weighted normal pseudosolution

Consider some properties of the weighted Moore–Penrose pseudoinverse.

Lemma 3 (see in [16]). Let $A,\mathrm{\Delta }A\in R^{m\times n}$, ${\mu }_i\left(A\right)$ and${\mu }_i\left(\overline{A}\right)$ denote the weighted singular values of A and $\overline{A}$ respectively. Then,

$${\mu }_i\left(A\right)-{‖\mathit{\Delta A}‖}_{\mathit{MN}}\le {\mu }_i\left(\overline{A}\right)\le {\mu }_i\left(A\right)+{‖\mathit{\Delta A}‖}_{\mathit{MN}}.$$E24

Lemma 4 (see [40]). Let $A,\mathrm{\Delta }A\in R^{m\times n}$, $rank\left(\overline{A}\right)=rank\left(A\right)$ and ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$. Then

$${‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\le \frac{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}{1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}.$$E25

Lemma 5. Let $G={\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}$, $\overline{А}=А+\mathrm{\Delta }A$ and $rank\left(\overline{A}\right)=rank\left(A\right)$. Then G can be represented as the sum of three matrices $G=G_1+G_2+G_3$, where

$$G_1=-{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+},$$E26

$$G_2=-\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=-\left(I-\overline{P}\right)\mathrm{\Delta }{A}^{\#}{\left(A_{\mathit{MN}}^{+}\right)}^{\#}{A}_{\mathit{MN}}^{+},$$E27

$$G_3={\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right).$$E28

Proof. Following [26], G can be represented as the sum of the following matrices.

$$\begin{array}{rcl}G& =& \left[\overline{P}+\left(I-\overline{P}\right)\right]\left({\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}\right)\left[Q+\left(I-Q\right)\right]=\\ & =& \overline{P}{\overline{A}}_{\mathit{MN}}^{+}Q+\overline{P}{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)-\overline{P}{A}_{\mathit{MN}}^{+}Q-\overline{P}{A}_{\mathit{MN}}^{+}\left(I-Q\right)+\left(I-\overline{P}\right){\overline{A}}_{\mathit{MN}}^{+}Q+\\ & & +\left(I-\overline{P}\right){\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)-\left(I-\overline{P}\right)A_{\mathit{MN}}^{+}Q+\left(I-\overline{P}\right)A_{\mathit{MN}}^{+}\left(I-Q\right).\end{array}$$E29

Since.

$$\overline{P}{\overline{A}}_{\mathit{MN}}^{+}={\overline{A}}_{\mathit{MN}}^{+},\left(I-\overline{P}\right){\overline{A}}_{\mathit{MN}}^{+}=0,A_{\mathit{MN}}^{+}Q=A_{\mathit{MN}}^{+},A_{\mathit{MN}}^{+}\left(I-Q\right)=0,$$E30

we obtain

$$\begin{array}{rcl}G& =& {\overline{A}}_{\mathit{MN}}^{+}Q+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)-\overline{P}{A}_{\mathit{MN}}^{+}+\left(I-\overline{P}\right){\overline{A}}_{\mathit{MN}}^{+}=\\ & =& \left({\overline{A}}_{\mathit{MN}}^{+}Q-\overline{P}{A}_{\mathit{MN}}^{+}\right)-\left(I-\overline{P}\right)A_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right).\end{array}$$E31

Consider each term in this equality separately

$$G_1={\overline{A}}_{\mathit{MN}}^{+}Q-\overline{P}{A}_{\mathit{MN}}^{+}={\overline{A}}_{\mathit{MN}}^{+}A{A}_{\mathit{MN}}^{+}-{\overline{A}}_{\mathit{MN}}^{+}\overline{A}{A}_{\mathit{MN}}^{+}={\overline{A}}_{\mathit{MN}}^{+}\left(A-\overline{A}\right)A_{\mathit{MN}}^{+}={\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}.$$E32

To estimate the second term, we use properties (1)

$$\begin{array}{rcl}{A}_{\mathit{MN}}^{+}& =& \left(A_{\mathit{MN}}^{+}A\right)A_{\mathit{MN}}^{+}=N^{-1}{\left(N{A}_{\mathit{MN}}^{+}A\right)}^T{A}_{\mathit{MN}}^{+}=\\ & =& N^{-1}{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}-N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}.\end{array}$$E33

Substituting (33) into the second term of (31) gives

$$G_2=\left(I-\overline{P}\right)A_{\mathit{MN}}^{+}=\left(I-\overline{P}\right)\left(N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}-N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}\right).$$E34

Since,

$$\begin{array}{l}\left(I-\overline{P}\right)N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}-{\overline{A}}_{\mathit{MN}}^{+}\overline{A}{N}^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=\\ =N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}-N^{-1}{\overline{A}}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=0\end{array}$$E35

we obtain

$$G_2=\left(I-\overline{P}\right)A_{\mathit{MN}}^{+}=-\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}=-\left(I-\overline{P}\right)\mathrm{\Delta }{A}^{\#}{\left(A_{\mathit{MN}}^{+}\right)}^{\#}{A}_{\mathit{MN}}^{+}.$$E36

Finally,

$$G={\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}=-{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}-\left(I-\overline{P}\right)\mathrm{\Delta }{A}^{\#}{\left(A_{\mathit{MN}}^{+}\right)}^{\#}{A}_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right).$$E37

Lemma 6 (see in [41]). If $rank\left(\overline{A}\right)=rank\left(A\right)=k$, then

$${‖\overline{Q}\left(I-Q\right)‖}_{\mathit{MM}}={‖Q\left(I-\overline{Q}\right)‖}_{\mathit{MM}},$$E38

where Q and $\overline{Q}$ are defined in (3).

Lemma 7. Let $A,\mathrm{\Delta }A\in R^{m\times n}$, $rank\left(\overline{A}\right)=rank\left(A\right)$ and ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$.

Then the relative estimate of the hereditary error of the weighted pseudoinverse matrix has the form

$$\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}\le C\frac{{\epsilon }_{A}h}{1-{\epsilon }_{A}h},$$E39

where $h=h\left(A\right)={‖A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ and the estimate of the absolute error

$${‖{\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\le C\frac{{\left({\epsilon }_{A}h\right)}^2}{1-{\epsilon }_{A}h},$$E40

moreover.

if A is not a full rank matrix, then $C=3$,

if $m>n=k$ or $n>m=k$, then $C=2$,

if $m=n=k$, then $C=1.$

Proof. To obtain estimates, we use the results of Lemma 5:

$${\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}=-{\overline{A}}_{\mathit{MN}}^{+}\mathit{\Delta A}{A}_{\mathit{MN}}^{+}-\left(I-\overline{P}\right)\Delta A^{\#}{\left(A_{\mathit{MN}}^{+}\right)}^{\#}{A}_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right).$$E41

Passing to the weighted norms, we obtain.

$${‖{\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\le {‖{\overline{A}}_{\mathit{MN}}^{+}\mathit{\Delta A}{A}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖\Delta A^{\#}{\left(A_{\mathit{MN}}^{+}\right)}^{\#}{A}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖{\overline{A}}_{\mathit{MN}}^{+}\overline{Q}\left(I-Q\right)‖}_{\mathit{NM}}.$$E42

Using the results of Lemma 6, we can estimate the last summand

$$\begin{array}{rcl}{‖{\overline{A}}_{\mathit{MN}}^{+}\overline{Q}\left(I-Q\right)‖}_N& =& {‖{\overline{A}}_{\mathit{MN}}^{+}\overline{A}{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)‖}_N\le {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\overline{Q}\left(I-Q\right)‖}_{\mathit{MM}}\\ & =& {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖Q\left(I-\overline{Q}\right)‖}_{\mathit{MM}}.\end{array}$$E43

According to (38) and (43), we can rewrite (42) in the form

$$\begin{array}{rcl}{‖{\overline{\mathit{A}}}_{\mathit{MN}}^{+}-{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}& \le & {‖{\overline{\mathit{A}}}_{\mathit{MN}}^{+}\mathit{\Delta A}{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖\mathit{\Delta }{\mathit{A}}^{\#}{\left({\mathit{A}}_{\mathit{MN}}^{+}\right)}^{\#}{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖{\mathit{A}}_{\mathit{MN}}^{+}\mathit{Q}\left(\mathit{I}-\overline{\mathit{Q}}\right)‖}_{\mathit{NM}}\le \\ & \le & {‖{\overline{\mathit{A}}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}^2& +{‖{\overline{\mathit{A}}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖{\mathit{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathit{\Delta A}‖}_{\mathit{MN}}.\end{array}$$E44

Using the results of Lemma 4, we obtain an estimate for the absolute error of the weighted pseudoinverse matrix A.

$$\begin{array}{l}{‖{\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\le \frac{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}{1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}({‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+\\ +{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}})=\frac{h{\epsilon }_A}{1-h{\epsilon }_A}\left(h{\epsilon }_A+h{\epsilon }_A+h{\epsilon }_A\right)=C\frac{{\left(h{\epsilon }_A\right)}^2}{1-h{\epsilon }_A},C=1,2,3.\end{array}$$E45

To estimate the relative error, we have

$$\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}\le 3{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{1}{1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}=C\frac{h{\epsilon }_A}{1-h{\epsilon }_A},C=1,2,3.$$E46

Let us estimate the error of the weighted minimum-norm least squares solution. Let us introduce the following notation:

$$\begin{array}{l}\alpha =\frac{{‖\mathit{\Delta b}‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N},\beta =\frac{{‖r‖}_M}{{‖x‖}_N{‖A‖}_{\mathit{MN}}},\gamma =\frac{{‖\overline{r}‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}{\alpha }_l=\frac{{‖\mathit{\Delta b}‖}_M}{{‖A‖}_{\mathit{MN}}{‖x_l‖}_N},\\ {\beta }_l=\frac{{‖r‖}_M}{{‖x_l‖}_N{‖A‖}_{\mathit{MN}}},{\gamma }_l=\frac{{‖{\overline{r}}_l‖}_M}{{‖A‖}_{\mathit{MN}}{‖x_l‖}_N},{\gamma }_k=\frac{{‖{\overline{r}}_k‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}.\end{array}$$E47

Consider the following three cases.

Case 1. The rank of the original matrix A remains the same under its perturbation, i.e., $rank\left(A\right)=rank\left(\overline{A}\right)$.

Theorem 7. Assume that ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, $rank\left(\overline{A}\right)=rank\left(A\right)$. Then

$$\frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}\le \frac{h}{1-h{\epsilon }_A}\left(2{\epsilon }_A+\alpha +h{\epsilon }_A\beta \right),$$E48

where $h\left(A\right)={‖A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is the weighted condition number of A, the symbols${\begin{array}{c}‖‖\end{array}}_{\mathit{MN}}$ and${\begin{array}{c}‖‖\end{array}}_{\mathit{NM}}$ denote the weighted matrix norms defined by Eq. (4)–(6), and $A_{\mathit{MN}}^{+}$ is the weighted Moore–Penrose pseudoinverse.

Proof. The error estimate follows from the relation:

$$\mathit{x}-\overline{\mathit{x}}=\left({\mathit{А}}_{\mathit{MN}}^{+}-{\overline{\mathit{А}}}_{\mathit{MN}}^{+}\right)\mathit{b}+{\overline{\mathit{А}}}_{\mathit{MN}}^{+}\left(\mathit{b}-\overline{\mathit{b}}\right).$$E49

For the error of the matrix pseudoinverse, we use the representation

$${\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}=-{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}-\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right).$$E50

Then,

$$\begin{array}{rcl}x-\overline{x}& =& \left[{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)\right]b+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)=\\ & =& {\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}b+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}b-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)b+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)\end{array}$$

$$\begin{array}{rc}={\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }\mathit{Ax}+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)b& +{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)\end{array}$$E51

Thus,

$$\overline{x}-x={\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }\mathit{Ax}+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)b+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right).$$E52

Passing to the weighted norms yields

$$\begin{array}{rcl}{‖\overline{x}-x‖}_N& =& {‖{\overline{A}}_{\mathit{MN}}\mathrm{\Delta }\mathit{Ax}+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)b+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)‖}_N\le \\ & \le & {‖{\overline{A}}_{\mathit{MN}}\mathrm{\Delta }\mathit{Ax}‖}_N+{‖\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}‖}_N+\Vert {\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)b{+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)\Vert }_N.\end{array}$$E53

By taking into account the relations

$$\left(I-Q\right)b=\left(I-Q\right)r=r,r=b-\mathit{A}x,x={\mathit{A}}_{\mathit{MN}}^{+}b$$E54

and applying Lemma 6, the weighted norm of each term in (27) can be rearranged as follows

1. $\begin{array}{l}{‖{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }\mathit{Ax}‖}_N=‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overline{A}}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{M}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }A{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{N}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}x‖\le \\ \le ‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overline{A}}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }A{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}x‖={‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖x‖}_N.\text{(55)}\end{array}$

2. $\begin{array}{l}{‖\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}‖}_N=‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(I-\overline{P}\right)N^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }{A}^T{M}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{A}_{\mathit{MN}}^{+T}{N}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{N}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}x‖\\ \le ‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(I-\overline{P}\right)N^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }A{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{A}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}x‖\\ ={‖\left(I-\overline{P}\right)‖}_{\mathit{NN}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖x‖}_N\text{(56)}\end{array}$

3. Using Lemma 6, and (28) we can write

$$\begin{array}{l}{‖{\overline{A}}_{\mathit{MN}}^{+}\overline{Q}\left(I-Q\right)b‖}_N={‖{\overline{A}}_{\mathit{MN}}^{+}\overline{A}{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right)r‖}_N\le {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\overline{Q}\left(I-Q\right)‖}_{\mathit{MM}}{‖r‖}_M\\ ={‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖Q\left(I-\overline{Q}\right)‖}_{\mathit{MM}}{‖r‖}_M\end{array}$$E57

where

$$\begin{array}{l}{‖Q\left(I-\overline{Q}\right)‖}_{\mathit{MM}}={‖A{A}_{\mathit{MN}}^{+}\left(I-\overline{Q}\right)‖}_{\mathit{MM}}=‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}A{A}_{\mathit{MN}}^{+}\left(I-\overline{Q}\right)M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖=\\ =‖M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\left(\mathit{MA}{A}_{\mathit{MN}}^{+}\right)}^T\left(I-\overline{Q}\right)M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖=\\ =‖M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\left(A_{\mathit{MN}}^{+}\right)}^T\left(A^T-{\overline{A}}^T\right)M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{M}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(I-\overline{Q}\right)M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖=\\ =‖M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\left(A_{\mathit{MN}}^{+}\right)}^T\mathrm{\Delta }{A}^{T}M\left(I-\overline{Q}\right)M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖\le ‖M^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\left(A_{\mathit{MN}}^{+}\right)}^T\mathrm{\Delta }{A}^T{M}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖=\\ =‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }A{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{N}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{A}_{\mathit{MN}}^{+}{M}^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖\le ‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathrm{\Delta }A{N}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{A}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖\le \\ \le {‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}.\end{array}$$E58

Substituting this result into (31) gives the inequality

$${‖{\overline{A}}_{\mathit{MN}}^{+}\overline{Q}\left(I-Q\right)b‖}_N\le {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖r‖}_M.$$E59

1. d.

   $$\begin{array}{l}{‖{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right)‖}_N=‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overline{A}}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{M}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(b-\overline{b}\right)‖\le \\ ‖N^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overline{A}}_{\mathit{MN}}^{+}{M}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}‖‖M^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(b-\overline{b}\right)‖={‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\left(b-\overline{b}\right)‖}_M\end{array}$$E60

Taking into account $‖I-\overline{P}‖<1$, and applying Lemma 4, we obtain the following weighted-norm estimate for the relative error

$$\begin{array}{l}\frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}\le \frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖x‖}_N}{{‖x‖}_N}+\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖x‖}_N}{{‖x‖}_N}+\\ +\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖r‖}_M}{{‖x‖}_N}+\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }b‖}_M}{{‖x‖}_N}\le \\ \le {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+\\ +\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖r‖}_M}{{‖x‖}_N}+\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }b‖}_M}{{‖x‖}_N}\le \\ \le \frac{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}}{1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}\left(2\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right.+\\ \left.+{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}\frac{{‖r‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right)\le \\ \le \frac{h\left(A\right)}{1-h\left(A\right){\epsilon }_A}\left(2{\epsilon }_A+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}+h\left(A\right){\epsilon }_A\frac{{‖r‖}_M}{{‖x‖}_N{‖A‖}_{\mathit{MN}}}\right).\end{array}$$E61

as required.

Specifically, if $M=I\in R^{m\times m}$ and $N=I\in R^{n\times n}$, then the estimates of the hereditary error of normal pseudosolutions of systems of linear algebraic equations follow from next theorem.

Theorem 8 (see in [32]). Let║$\mathit{\Delta А}$║║${А}^{+}$║< 1, $rank\left(\overline{A}\right)=rank\left(A\right)=k$. Then

$$\frac{‖x-\overline{x}‖}{‖x‖}\le \frac{h}{1-h{\epsilon }_{\mathrm{А}}}\left(2{\epsilon }_{\mathrm{А}}+{\epsilon }_{b_k}+h{\epsilon }_{\mathrm{А}}\frac{‖b-b_k‖}{‖b_k‖}\right),$$E62

where $b_k$ is the projection of the right-hand side of problem (8) onto the principal left singular subspace of the matrix A [42], i.е., $b_k\in Im\mathrm{A}$, $h=h\left(A\right)=‖A‖‖A^{+}‖$ is condition number of A, the symbol $\begin{array}{c}‖‖\end{array}$, unless otherwise stated, denotes the Euclidean vector norm and the corresponding spectral matrix norm, ${А}^{+}$ is the Moore–Penrose pseudoinverse.

Case 2. The rank of the perturbed matrix is larger than that of the original matrix A, i.e.$rank\left(\overline{A}\right)>\mathit{rank}\left(A\right)=k$.

Define the idempotent matrices:

$$P=A_{\mathit{MN}}^{+}A,Q=A{A}_{\mathit{MN}}^{+},{\overline{P}}_k={{\overline{A}}_k}_{\mathit{MN}}^{+}\overline{A},{\overline{Q}}_k=\overline{A}{{\overline{A}}_k}_{\mathit{MN}}^{+},$$E63

where k is the rank of A.

Theorem 9. Assume that ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<\frac{1}{2}$, $rank\left(\overline{A}\right)>rank\left(A\right)=k$. Then

$$\frac{{‖x-{\overline{x}}_k‖}_N}{{‖x‖}_N}\le \frac{h}{1-2h{\epsilon }_{\mathrm{А}}}\left(2{\epsilon }_{\mathrm{А}}+\alpha +h{\epsilon }_{\mathrm{А}}\beta \right).$$E64

where $h\left(A\right)={‖A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is the weighted condition number of A, the symbols ${\begin{array}{c}‖‖\end{array}}_{\mathit{MN}}$ and${\begin{array}{c}‖‖\end{array}}_{\mathit{NM}}$ denote the weighted matrix norms defined Eq. (4)–(6), and $A_{\mathit{MN}}^{+}$ is the weighted Moore–Penrose pseudoinverse.

Proof. The desired estimate is derived using the method of [32], which is based on the singular value decomposition of matrices. Specifically, $\overline{А}$ is represented as a weighted singular value decomposition:

$$\overline{А}=\overline{U}\overline{D}{\overline{V}}^T.$$E65

Along with (38), we consider the decomposition

$${\overline{А}}_k=\overline{U}{\overline{D}}_k{\overline{V}}^T,$$E66

where ${\overline{D}}_k$ is a rectangular matrix whose first k diagonal elements are nonzero and equal to the corresponding elements of $\overline{D}$, while all the other elements are zero.

The weighted minimum-norm least squares solution to problem (10) is approximated by the weighted minimum-norm least squares solution ${\overline{x}}_k$ to the problem

$$\underset{x\in C}{min}{‖\overline{x}‖}_N,C=\left\{x|{‖{\overline{A}}_k\overline{x}-\overline{b}‖}_M=min\right\}.$$E67

The matrix ${\overline{A}}_k$ is defined by (48) and has the same rank k as the matrix of the unperturbed problem.

Thus, the error estimation of the least-squares solution for matrices with a modified rank is reduced to the case of the same rank. This fact is used to estimate ${‖x-{\overline{x}}_k‖}_N/{‖x‖}_N$. The error of the weighted pseudoinverse matrix then becomes:

$$\begin{array}{l}{G}_k=\left[{\overline{P}}_k+\left(I-{\overline{P}}_k\right)\right]\left({{\overline{A}}_k}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}\right)\left[Q+\left(I-Q\right)\right]={\overline{P}}_k{{\overline{A}}_k}_{\mathit{MN}}^{+}Q+{\overline{P}}_k{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q\right)-\\ -{\overline{P}}_k{A}_{\mathit{MN}}^{+}Q-{\overline{P}}_k{A}_{\mathit{MN}}^{+}\left(I-Q\right)-\left(I-{\overline{P}}_k\right){{\overline{A}}_k}_{\mathit{MN}}^{+}Q+\left(I-{\overline{P}}_k\right){{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q\right)-\\ -\left(I-{\overline{P}}_k\right)A_{\mathit{MN}}^{+}Q+\left(I-{\overline{P}}_k\right)A_{\mathit{MN}}^{+}\left(I-Q\right)=\left({{\overline{A}}_k}_{\mathit{MN}}^{+}Q-{\overline{P}}_k{A}_{\mathit{MN}}^{+}\right)-\left(I-{\overline{P}}_k\right)A_{\mathit{MN}}^{+}+\\ +{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q\right)={{\overline{A}}_k}_{\mathit{MN}}^{+}A{A}_{\mathit{MN}}^{+}-{{\overline{A}}_k}_{\mathit{MN}}^{+}\overline{A}{A}_{\mathit{MN}}^{+}-\left(I-{\overline{P}}_k\right)A_{\mathit{MN}}^{+}+{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q\right)=\\ ={{\overline{A}}_k}_{\mathit{MN}}^{+}\left(A-\overline{A}\right)A_{\mathit{MN}}^{+}-\left(I-{\overline{P}}_k\right)A_{\mathit{MN}}^{+}+{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q\right),\end{array}$$E68

Applying Lemma 5 yields

$$G_k={{\overline{A}}_k}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}=-{{\overline{A}}_k}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A}_{\mathit{MN}}^{+}-\left(I-{\overline{P}}_k\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}N{A}_{\mathit{MN}}^{+}+{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q_k\right).$$E69

For the error of the WLS solution, we obtain

$${\overline{x}}_k-x={{\overline{A}}_k}_{\mathit{MN}}^{+}\mathrm{\Delta }\mathit{Ax}+\left(I-{\overline{P}}_k\right)N^{-1}\mathrm{\Delta }{A}^T{A}_{\mathit{MN}}^{+T}\mathit{Nx}-{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(I-Q_k\right)b+{{\overline{A}}_k}_{\mathit{MN}}^{+}\left(b-\overline{b}\right).$$E70

Passing to the weighted norms and applying Lemma 4 gives

$$\begin{array}{l}\frac{{‖x-{\overline{x}}_k‖}_N}{{‖x‖}_N}\le \frac{{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖x‖}_N}{{‖x‖}_N}+\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖x‖}_N}{{‖x‖}_N}+\\ +\frac{{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖r‖}_M}{{‖x‖}_N}+\frac{{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }b‖}_M}{{‖x‖}_N}\le \\ \le {‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+\\ +\frac{{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖r‖}_M}{{‖x‖}_N}+\frac{{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }b‖}_M}{{‖x‖}_N}\le \\ \le \frac{{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}}{1-{‖\mathrm{\Delta }{A}_k‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}\left(2\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right.+\\ \left.+{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}\frac{{‖r‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right)\end{array}.$$E71

Let estimate $\mathrm{\Delta }{A}_k=A-{\overline{A}}_k$:

$$\begin{array}{l}{‖\mathrm{\Delta }{A}_k‖}_{\mathit{MN}}={‖{\overline{A}}_k-A‖}_{\mathit{MN}}={‖{\overline{A}}_k-\overline{A}+\mathrm{\Delta }A‖}_{\mathit{MN}}\le {‖{\overline{A}}_k-\overline{A}‖}_{\mathit{MN}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}=\\ ={‖\overline{U}\left(\begin{array}{cc}0& 0\\ 0& D_{k+1}\end{array}\right){\overline{V}}^T‖}_{\mathit{MN}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}\le 2{‖\mathrm{\Delta }A‖}_{\mathit{MN}}.\end{array}$$E72

Moreover, the theorem condition ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<\frac{1}{2}$ leads to ${‖\mathrm{\Delta }{A}_k‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, which is necessary for expression (51) to be well defined. In view of this, (51) yields estimate (33) for the error of the minimum-norm weighted least squares solution.

Specifically, if $M=I\in R^{m\times m}$ and $N=I\in R^{n\times n}$, then the estimates of the hereditary error of normal pseudosolutions of systems of linear algebraic equations for case $rank\left(\overline{A}\right)>rank\left(A\right)=k$ follows from next theorem.

Theorem 10 (see in [32]). Let $‖\mathrm{\Delta }A‖‖A^{+}‖<\frac{1}{2}$, $rank\left(\overline{A}\right)>rank\left(A\right)=k$. Then

$$\frac{‖x-{\overline{x}}_k‖}{‖x‖}\le \frac{h}{1-2h{\epsilon }_{\mathrm{А}}}\left(2{\epsilon }_{\mathrm{А}}+{\epsilon }_{b_k}+h{\epsilon }_{\mathrm{А}}\frac{‖b-b_k‖}{‖b_k‖}\right).$$E73

Case 3. The rank of the original matrix is larger than that of the perturbed matrix, i.e., $rank\left(A\right)>rank\left(\overline{A}\right)=l$.

By analogy with (33), we define the idempotent matrices:

$$P_l={A_l}_{\mathit{MN}}^{+}A,Q_l=A{A_l}_{\mathit{MN}}^{+},\overline{P}={\overline{A}}_{\mathit{MN}}^{+}\overline{A},\overline{Q}=\overline{A}{\overline{A}}_{\mathit{MN}}^{+},$$E74

Theorem 11. Assume that $rank\left(A\right)>rank\left(\overline{A}\right)=l$, $\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{\mu }_l}<\frac{1}{2}$. Then,

$$\frac{{‖x_l-\overline{x}‖}_N}{{‖x_l‖}_N}\le \frac{{\mu }_1/{\mu }_l}{1-2{‖\mathrm{\Delta }A‖}_{\mathit{MN}}/{\mu }_l}\left(2{\epsilon }_A+{\alpha }_l+\frac{{\mu }_1}{{\mu }_l}{\epsilon }_{\mathrm{А}}{\beta }_l\right),$$E75

where ${\mu }_i$ are the weighted singular values of A.

Proof. Along with (9), we consider the problem

$$\underset{x\in C}{\mathit{min}}{‖x_l‖}_N,C=\left\{x|{‖A_{l}x-b‖}_M=min\right\}$$E76

with the matrix $A_l=U{D}_l{V}^T$ of rank l.

Similarly, writing (27) for problems (10) and (54), whose matrix ranks coincide, we obtain

$$\begin{array}{l}{G}_l={\overline{A}}_{\mathit{MN}}^{+}-{A_l}_{\mathit{MN}}^{+}=-{\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }A{A_l}_{\mathit{MN}}^{+}-\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A_l}_{\mathit{MN}}^{+T}N{A_l}_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q_l\right),\end{array}$$E77

$$\begin{array}{l}\overline{x}-x_l={\overline{A}}_{\mathit{MN}}^{+}\mathrm{\Delta }\mathit{Ax}+\left(I-\overline{P}\right)N^{-1}\mathrm{\Delta }{A}^T{A_l}_{\mathit{MN}}^{+T}\mathit{Nx}-{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q_l\right)b+{\overline{A}}_{\mathit{MN}}^{+}\left(b-\overline{b}\right).\end{array}$$E78

Applying Lemma 4 and passing to the weighted norms yields the estimate

$$\begin{array}{l}\frac{{‖x_l-\overline{x}‖}_N}{{‖x‖}_N}\le {‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}+\\ +\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖r‖}_M}{{‖x‖}_N}+\frac{{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖\mathrm{\Delta }b‖}_M}{{‖x‖}_N}\le \\ \le \frac{{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}}{1-{‖\mathrm{\Delta }{A}_l‖}_{\mathit{MN}}{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}}\left(2\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right.+\\ \left.+{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖A‖}_{\mathit{MN}}\frac{{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{‖A‖}_{\mathit{MN}}}\frac{{‖r‖}_M}{{‖A‖}_{\mathit{MN}}{‖x‖}_N}\right),\end{array}$$E79

which implies (52). This completes the proof of Theorem 11.

For approximately given initial data, the rank of the original matrix should be specified as the numerical rank of the matrix (see in [28]).

Specifically, if $M=I\in R^{m\times m}$ and $N=I\in R^{n\times n}$, then the estimates of the hereditary error of normal pseudosolutions of systems of linear algebraic equations for case $rank\left(A\right)>rank\left(\overline{A}\right)=l$ follows from next theorem.

Theorem 12 (see in [32]). Let $rank\left(A\right)>rank\left(\overline{A}\right)=l$, $\frac{‖\mathit{\Delta A}‖}{{\mu }_l}<\frac{1}{2}$. Then

$$\frac{‖x_l-\overline{x}‖}{‖x_l‖}\le \frac{{\mu }_1/{\mu }_l}{1-2‖∆A‖/{\mu }_l}\left(2{\epsilon }_A+{\epsilon }_{b_l}+{\epsilon }_{\mathrm{А}}\frac{{\mu }_1}{{\mu }_l}\frac{‖b-b_l‖}{‖b_l‖}\right),$$E80

where $x_l$ is the projection of the normal pseudosolution of problem (8) onto the right principal singular subspace of the matrix A of dimension l, $b_l$ is projection of the right-hand side b onto the principal left singular subspace of dimension l of the matrix A, ${\mu }_i$ is singular values of the matrix А.

3.2 Estimates of the hereditary error of a weighted normal pseudosolution for full rank matrices

For matrices of full rank, it is essential that their rank does not change due to the perturbation of the elements if the condition ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$ is met.

In addition, in what follows we will use the following property of matrices of full rank [28]

$$A_{\mathit{MN}}^{+}={\left(A^T\mathit{MA}\right)}^{-1}{A}^{T}M\text{for}m\ge n\text{and}{A}_{\mathit{MN}}^{+}=N^{-1}{A}^T{\left(A{N}^{-1}{A}^T\right)}^{-1}\text{for}n\ge m.$$E81

If $m\ge n$, then problem (9) is reduced to a problem of the form

$$\underset{x\in R^n}{\mathit{min}}{‖\mathit{Ax}-b‖}_M.$$E82

For such a problem, the following theorem is true.

Theorem 13. Let ${‖\mathit{\Delta A}‖}_{\mathit{MI}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{IM}}<1$, $m>n=k$. Then

$$\frac{‖x-\overline{x}‖}{‖x‖}\le \frac{h}{1-h{\epsilon }_A}\left({\epsilon }_A+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MI}}‖x‖}+h{\epsilon }_A\frac{{‖r‖}_M}{‖x‖{‖A‖}_{\mathit{MI}}}\right),$$E83

where $h={‖A‖}_{\mathit{MI}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{IM}}$.

Proof. To prove Theorem 13, as before, we will use relation (49). By (81) $\overline{P}={\overline{A}}_{\mathit{MN}}^{+}\overline{A}=I$, so that from (50) we have the equality

$${\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}=-{\overline{A}}_{\mathit{MN}}^{+}\mathit{\Delta A}{A}_{\mathit{MN}}^{+}+{\overline{A}}_{\mathit{MN}}^{+}\left(I-Q\right),$$E84

using which we obtain (83).

If $n\ge m$, then problem (9) is reduced to a problem of the form

$$\underset{\mathit{x}\in \mathit{C}}{\mathit{min}}{‖\mathit{x}‖}_{\mathit{N}},\mathit{C}=\left\{\mathit{x}|\mathit{Ax}=\mathit{b}\right\}$$E85

and the following theorem holds for it.

Theorem 14. Let ${‖\mathit{\Delta A}‖}_{\mathit{\text{IN}}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NI}}<1$, $n>m=k$. Then

$$\frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}\le \frac{h}{1-h{\epsilon }_A}\left(2{\epsilon }_A+\frac{‖\mathrm{\Delta }b‖}{{‖A‖}_{\mathit{\text{IN}}}{‖x‖}_N}\right),$$E86

where $h={‖A‖}_{\mathit{\text{IN}}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NI}}$.

Proof. Since in this case $Q=A{A}_{\mathit{MN}}^{+}=I$, then the expression for ${\overline{A}}_{\mathit{MN}}^{+}-A_{\mathit{MN}}^{+}$ by (81) takes the form

$${\overline{\mathit{A}}}_{\mathit{MN}}^{\mathit{+}}\mathit{-}{\mathit{A}}_{\mathit{MN}}^{\mathit{+}}\mathit{=}\mathit{-}{\overline{\mathit{A}}}_{\mathit{MN}}^{\mathit{+}}\mathit{\Delta A}{\mathit{A}}_{\mathit{MN}}^{\mathit{+}}\mathit{-}\left(\mathit{I}\mathit{-}\overline{\mathit{P}}\right){\mathit{N}}^{\mathit{-}1}\mathit{\Delta }{\mathit{A}}^{\mathit{T}}{\mathit{A}}_{\mathit{MN}}^{\mathit{+}\mathit{T}}\mathit{N}{\mathit{A}}_{\mathit{MN}}^{\mathit{+}}\mathit{.}$$E87

Further calculations are similar to the previous ones. As a result, we come to estimate (86).

Remark 3. The relationship between the condition number of the problem with exact initial data h(A) and the condition number of the matrix of the system with approximately given initial data $h\left(\overline{A}\right)$ is established by the estimates

$$\begin{array}{l}{\sigma }_k-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}\le {\overline{\sigma }}_k\le {\sigma }_k+{‖\mathrm{\Delta }A‖}_{\mathit{MN}},{\sigma }_1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}\le {\overline{\sigma }}_1\le {\sigma }_1+{‖\mathrm{\Delta }A‖}_{\mathit{MN}},\\ \frac{{\sigma }_1-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{\sigma }_k+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}\le \frac{{\overline{\sigma }}_1}{{\overline{\sigma }}_k}\le \frac{{\sigma }_1+{‖\mathrm{\Delta }A‖}_{\mathit{MN}}}{{\sigma }_k-{‖\mathrm{\Delta }A‖}_{\mathit{MN}}},\frac{1-{\epsilon }_A}{1+{\epsilon }_{A}h}\le \frac{h\left(\mathrm{А}\right)}{h\left(\overline{\mathrm{А}}\right)}\le \frac{1+{\epsilon }_A}{1-{\epsilon }_{A}h},\end{array}$$E88

which are easy to obtain for the weighted matrix norm based on the perturbation theory for symmetric matrices.

Lemma 7. Let $A,\mathrm{\Delta }A\in R^{m\times n}$, $rank\left(\overline{A}\right)=rank\left(A\right)$ and ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$. Then the estimate of the relative error of the condition number of the matrix A has the form

$$\left|\frac{\overline{h}-h}{h}\right|\le {\epsilon }_A\frac{1+h}{1-{\epsilon }_{A}h}$$E89

where $h=h\left(A\right)={‖A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is weighted condition number of matrix A, $\overline{h}=\overline{h}\left(A\right)={‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is weighted condition number of the perturbed matrix $\overline{A}=A+\mathrm{\Delta }A$.

Proof of Lemma 7 is easy to obtain using the inequality (25).

Theorem 15. Let ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, ${‖\mathit{\Delta A}‖}_{\mathit{MN}}\le {\epsilon }_{\overline{A}}{‖\overline{A}‖}_{\mathit{MN}}$, $rank\left(\overline{A}\right)=rank\left(A\right)$. Then,

$$\frac{{‖\overline{x}-x‖}_N}{{‖\overline{x}‖}_N}\le \frac{h\left(\overline{A}\right)}{1-h\left(\overline{A}\right){\epsilon }_{\overline{A}}}\left(2{\epsilon }_{\overline{A}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖\overline{A}‖}_{\mathit{MN}}{‖\overline{x}‖}_N}+h\left(\overline{A}\right){\epsilon }_{\overline{A}}\frac{{‖\overline{r}‖}_M}{{‖\overline{x}‖}_N{‖\overline{A}‖}_{\mathit{MN}}}\right),$$E90

where $h\left(\overline{A}\right)={‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is weighted matrix condition number $\overline{A}$, the symbols ${‖‖}_{\mathit{MN}}$ and ${‖‖}_{\mathit{NM}}$ denote the weighted matrix norms defined by Eq. (4)–(6) and $A_{\mathit{MN}}^{+}$ is the weighted Moore–Penrose pseudoinverse.

Thus, estimates of the hereditary error, the right-hand side of which is determined by approximate data, can be obtained without inequalities (88). Estimates similar to (90) can be obtained for all the previously considered cases.

Remark 4. Under the conditions of Theorem 15, using the inequality

$$\frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}\le \frac{{‖x-\overline{x}‖}_N}{{‖\overline{x}‖}_N}\left(1+\frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}\right)$$E91

and inequality (90) we arrive at the estimate in the following theorem.

Theorem 16. Let ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, ${‖\mathrm{\Delta }A‖}_{\mathit{MN}}\le {\epsilon }_{\overline{A}}{‖\overline{A}‖}_{\mathit{MN}}$, $rank\left(\overline{A}\right)=rank\left(A\right)$. Then

$$\frac{{‖\overline{x}-x‖}_N}{{‖x‖}_N}\le \frac{\beta }{1-\beta },\beta =\frac{h\left(\overline{A}\right)}{1-h\left(\overline{A}\right){\epsilon }_{\overline{A}}}\left(2{\epsilon }_{\overline{A}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖\overline{A}‖}_{\mathit{MN}}{‖\overline{x}‖}_N}+h\left(\overline{A}\right){\epsilon }_{\overline{A}}\frac{{‖\overline{r}‖}_M}{{‖\overline{x}‖}_N{‖\overline{A}‖}_{\mathit{MN}}}\right).$$E92

Estimates similar to (92) can be obtained for all the previously considered cases.

4. Research and solution of the WLS problem with approximate initial data

4.1 Investigation of the properties of WLS problem with approximate initial data

In the study of the mathematical properties of the weighted least squares problem with approximate initial data associated with computer realization as an approximate model in (10), (11) we will understand exactly the computer model of the problem. We will assume that the error of the initial data $\mathrm{\Delta }A$, $\mathrm{\Delta }b$, in this case, contains in addition to everything, the error that occurs when the matrix coefficients are written to the computer memory or its computing.

Matrix of full rank within the error of initial data, we assume a matrix that cannot change the rank of $\mathrm{\Delta }A$ change in its elements.

Matrix of full rank within the machine precision, we assume a matrix that cannot change the rank when you change the elements within the machine precision.

Lemma 8. If $rank\left(A\right)=min\left(m,n\right)$, and

$${‖\mathrm{\Delta }A‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1,$$E93

Then $rank\left(\overline{A}\right)=rank\left(A\right)$.

Proof. For proof, let, for example, $rank\left(A\right)=m$ . Taking equal${‖\mathrm{\Delta }A‖}_{\mathit{MN}}=\epsilon$, in equality (93) can be rewritten as $\frac{\epsilon }{{\mu }_m}<1$, which is equivalent

$${\mu }_m-\epsilon >0.$$E94

Let ${\overline{\mu }}_m$—m-weighted singular value of perturbed matrix $\overline{A}=A+\mathrm{\Delta }A$. According to Lemma 3, we can write ${\overline{\mu }}_m\ge {\mu }_m-\epsilon$. Then, taking into account (94), we obtain ${\overline{\mu }}_m\ge {\mu }_m-\epsilon >0$.

Therefore $rank\left(\overline{A}\right)\ge m$, whence we come to the conclusion that $rank\left(\overline{A}\right)=m$, i.е. $rank\left(\overline{A}\right)=rank\left(A\right)$.

Taking into account the results of Lemma 8, the computer algorithm for studying rank completeness is reduced to checking the two relations

$${\epsilon }_{\overline{A}}h\left(\overline{A}\right)<1,$$E95

$$1.0+\frac{1}{h\left(\overline{A}\right)}\ne 1.0$$E96

where $h\left(\overline{A}\right)={‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}$ is weighted condition number of matrix $\overline{A}$.

The fulfillment of the first condition (95) guarantees that the matrix has a full rank and is within the accuracy of the initial data, and the second (96), which is performed in floating-point arithmetic, means that the matrix has a full rank within the machine precision.

Under these conditions, the solution of the machine problem exists, it is unique and stable. Such a machine problem should be considered as correctly posed within the accuracy of initial data.

Otherwise, the matrix of the perturbed system may be a matrix, not full rank and, therefore, the machine model of the problem (10), (11) should be considered as ill-posed. A key factor in studying the properties of a machine model is the criterion of the correctness of the problem. Thereby, a useful fact is that the condition for studying the machine model of problem (96) includes the value inverse to $h\left(\overline{A}\right)$. As a result, for large condition numbers of conditionality does not occur an overflow in order. And the disappearance of the order for $1.0/h\left(\overline{A}\right)$ for large condition numbers is not fatal: the machine result is assumed to be equal to zero, which allows us to make the correct conclusion about the loss of the rank of the matrix of the machine problem.

To analyze the properties of a machine model of problems with matrices of incomplete rank under conditions of approximate initial data, a fundamental role is played definition of the rank of a matrix.

The rank of the matrix in the condition of approximate the initial data (effective rank or $\delta$ -rank) is

$$rank\left(A,\delta \right)=\underset{{‖A-B‖}_{\mathit{MN}}\le \delta }{min}rank\left(B\right).$$E97

This means that the $\delta$-rank of the matrix is equal to the minimum rank among all matrices in the neighborhood ${‖A-B‖}_{\mathit{MN}}\le \delta$.

From [28] that if $r\left(\delta \right)$ is the $\delta$-rank of the matrix, then

$${\mu }_1\ge \dots \ge {\mu }_{r\left(\delta \right)}>\delta \ge {\mu }_{r\left(\delta \right)+1}\ge \dots \ge {\mu }_p,p=min\left(m,n\right).$$E98

The practical algorithm for finding $\delta$—rank can be defined as follows: find the value of r is equal to the largest value of i, for which the inequality is fulfilled

$$\frac{\delta }{{\mu }_i}<1,{\mu }_i\ne 0,i=1,2..$$E99

Using the effective rank of a matrix, can always find the number of a stable projection that approximates the solution or projection

To analyze the rank of a matrix of values within the machine precision value $\delta$ can be attributed to machine precision, for example, setting it equal $\text{macheps}‖B‖$.

4.2 Algorithm for finding a weighted normal pseudosolution of the weighted least squares problem with approximate initial data

The algorithm is based on weighted singular value decomposition of matrices (Lemma 1).

Let $A\in R^{m\times n}$ and $rank\left(A\right)=k$, M- and N-positive-defined matrices of order m and n, respectively.

To solve the ill-posed problems in the formulation (10), (11), the algorithm for obtaining an approximate normal pseudosolution of system (9), depending on the ratio of the ranks of the matrices A and $\overline{A}$ is reduced to the following three cases.

1. If the rank of the matrix has not changed $\mathit{rank}\left(\overline{A}\right)=\mathit{rank}\left(A\right)=k$, an approximate weighted normal pseudosolution is constructed by the formula

   $$\overline{x}={\overline{A}}_{\mathit{MN}}^{+}\overline{b},$$E100

   where ${\overline{A}}_{\mathit{MN}}^{+}$ is represented as a weighted singular value decomposition (7).

   In this case, the weighted normal pseudosolution of system (9) is approximated by the weighted normal pseudosolution of system (10) and, if ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, then the error of the solution is estimated by formula (48).

   If the rank of the matrix is complete and conditions (95), (96) are satisfied, the rank of the matrix does not change, and to estimate the error, one can use formulas (100), (48).

2. Matrix rank increased $rank\left(\overline{A}\right)>rank\left(A\right)=k$ . An approximate weighted normal pseudosolution is constructed by the formula

   $${\overline{x}}_k={{\overline{A}}_k}_{\mathit{MN}}^{+}\overline{b}.$$E101

   Weighted pseudoinverse matrix ${{\overline{A}}_k}_{\mathit{MN}}^{+}$ is defined as follows

   $${{\overline{A}}_k}_{\mathit{MN}}^{+}=N^{-1}\overline{V}{\overline{D}}_k^{+}{\overline{U}}^{T}M,$$E102

   where ${\overline{D}}_k$ is a rectangular matrix, the first k diagonal elements of which are nonzero and coincide with the corresponding elements of the matrix $\overline{D}$ from (7), and all other elements are equal to zero.

   In this case, the weighted normal pseudosolution of system (9) is approximated by the projection of the weighted normal pseudosolution of system (10) onto the right principal weighted singular subspace of dimension k of the matrix $\overline{A}$ and, if ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<\frac{1}{2}$, then the error of the solution is estimated by formula (64).

3. If the rank of the matrix has decreased $rank\left(A\right)>rank\left(\overline{A}\right)=l$, an approximation to the projection of a weighted normal pseudosolution of problem (9) is constructed using formula (100). In this case, the projection of the weighted normal pseudosolution of system (9) onto the principal right weighted singular subspace of dimension l of the matrix A is approximated by the weighted normal pseudosolution of system (10) and, if $\frac{{‖\mathit{\Delta A}‖}_{\mathit{MN}}}{{\mu }_l}<\frac{1}{2}$, the projection error is estimated by formula (75).

Remark 5. If the rank of the original matrix is unknown, then the $\delta$-rank should be taken as the projection number in (101). In this case, it is guaranteed that a stable approximation is found either to a weighted normal pseudosolution or to a projection, respectively, with error estimates.

If the rank of the original matrix is known, then it is guaranteed to find an approximation to the weighted normal pseudosolution with appropriate estimates.

Remark 6. Because of the zero columns in the matrix $D^{+}$, only the largest first n columns of the matrix U can actually contribute to the product (100). Moreover, if some of the weighted singular numbers are equal to zero, then less than n columns of U are needed. If kp is the number of nonzero weighted singular numbers, then U can be reduced to the sizes $m\times \mathit{kp}$, $D^{+}$—to the sizes $\mathit{kp}\times \mathit{kp}$, VT—up to size $\mathit{kp}\times n$. Formally, such matrices U and V are not M-orthogonal and $N^{-1}$-orthogonal, respectively, since they are not square. However, their columns are weighted orthonormal systems of vectors.

5. Analysis of the reliability of computer solutions to the WLS problem with approximate initial data

5.1 Estimates of the total error of a weighted normal pseudosolution for matrices of arbitrary rank

Estimates of the total error take into account both the hereditary error due to the error in the initial data and the computational error due to an approximate method for determining the solution to the problem. In this case, the method of obtaining a solution is not taken into account. The computational error can be a consequence of both an approximate method of obtaining a solution and an error due to inaccuracy in performing arithmetic operations on a computer. The residual vector $\overline{r}=\overline{A}\overline{\overline{x}}-\overline{b}$ takes into account the overall effect of these errors.

Let us obtain estimates for the total error of the weighted normal pseudosolution using the previously introduced notation (47). Let us consider three cases.

Case 1. The rank of the original matrix A remains the same under its perturbation, i.e., $rank\left(A\right)=rank\left(\overline{A}\right)$.

Theorem 17. Assume that ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, $rank\left(\overline{A}\right)=rank\left(A\right)=k$ and let $\overline{x}\in \mathbb{R}\left({\overline{A}}^{\#}\right)$. Then

$$\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{h}{1-h{\epsilon }_{\mathrm{А}}}\left(2{\epsilon }_{\mathrm{А}}+\alpha +h{\epsilon }_{\mathrm{А}}\beta +\gamma \right).$$E103

Proof. For the hereditary error, in this case, estimate (48) holds.

To estimate the computational error $\overline{x}-\overline{\overline{x}}$, we use the relation

$$\overline{A}\left(\overline{x}-\overline{\overline{x}}\right)=\overline{r}={\overline{b}}_k-\overline{A}\overline{\overline{x}},$$E104

where${\overline{b}}_k$ is projection of the vector $\overline{b}$ on the main left weighted singular subspace of the matrix $\overline{A}$, i.е. ${\overline{b}}_k\in \mathbb{R}\left(\overline{A}\right)$.

Considering that $\overline{x}-\overline{\overline{x}}\in \mathbb{R}\left({\overline{A}}^{\#}\right)$ and the fact that ${\overline{A}}_{\mathit{MN}}^{+}\overline{A}$ is a projector in $\mathbb{R}\left({\overline{A}}^{\#}\right)$, we have

$${\overline{\mathit{A}}}_{\mathit{MN}}^{+}\overline{\mathit{A}}\left(\overline{\mathit{x}}-\overline{\overline{\mathit{x}}}\right)=\overline{\mathit{x}}-\overline{\overline{\mathit{x}}}={\overline{\mathit{A}}}_{\mathit{MN}}^{+}\overline{\mathit{r}}.$$E105

From this, we obtain an estimate of the computational error

$$\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖\overline{x}‖}_N}\le {‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{{‖\overline{r}‖}_M}{{‖{\overline{b}}_k‖}_M}.$$E106

An estimate of the total error of the normal pseudosolution follows from the relations

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{{‖x-\overline{x}‖}_N}{{‖x‖}_N}+\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖x‖}_N},$$E107

$$\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le {‖A‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{{‖\overline{r}‖}_M}{{‖b_k‖}_M}$$E108

and estimates (48), (25). The theorem is proved.

Case 2. The rank of the perturbed matrix is larger than that of the original matrix A, i.e., $rank\left(\overline{A}\right)>rank\left(A\right)=k$.

If${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, then from [26], it follows that the rank of the perturbed matrix cannot decrease.

Theorem 18. Assume that ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<\frac{1}{2}$, $rank\left(\overline{A}\right)>rank\left(A\right)=k$ and let$\overline{x}\in \mathbb{R}\left({\overline{A}}_k^{\#}\right)$. Then

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{h}{1-h{\epsilon }_A}\left(2{\epsilon }_A+\alpha +h{\epsilon }_A\beta +{\gamma }_k\right).$$E109

Proof. To estimate the computational error ${‖{\overline{x}}_k-\overline{\overline{x}}‖}_N$, we use the fact that ${\overline{A}}_k{\overline{x}}_k={\overline{b}}_k$. Then for arbitrary vector $\overline{x}\in \mathbb{R}\left({\overline{A}}_k^{\#}\right)$

$${\overline{\mathrm{А}}}_k\left({\overline{x}}_k-\overline{\overline{x}}\right)={\overline{r}}_k={\overline{b}}_k-{\overline{A}}_k\overline{\overline{x}},{{\overline{A}}_k}^{+}{\overline{\mathrm{А}}}_k\left({\overline{x}}_k-\overline{\overline{x}}\right)={{\overline{A}}_k}^{+}{\overline{r}}_k.$$E110

Considering the fact that $x_k-\overline{\overline{x}}\in \mathbb{R}\left({{\overline{{\rm A}}}_k}^{\#}\right)$, and operator ${{\overline{A}}_k}_{\mathit{MN}}^{+}{\overline{A}}_k$ is the projection operator in $\mathbb{R}\left({{\overline{{\rm A}}}_k}^{\#}\right)$, we obtain

$${{\overline{A}}_k}_{\mathit{MN}}^{+}{\overline{A}}_k\left({\overline{x}}_k-\overline{\overline{x}}\right)={\overline{x}}_k-\overline{\overline{x}}={{\overline{A}}_k}_{\mathit{MN}}^{+}{\overline{r}}_k,{\overline{x}}_k-\overline{\overline{x}}={{\overline{A}}_k}_{\mathit{MN}}^{+}{\overline{r}}_k.$$E111

Hence follows an estimate of the computational error for the projection of the normal pseudosolution

$$\frac{{‖{\overline{x}}_k-\overline{\overline{x}}‖}_N}{{‖{\overline{x}}_k‖}_N}\le {‖{\overline{A}}_k‖}_{\mathit{MN}}{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{{‖{\overline{r}}_k‖}_M}{{‖{\overline{b}}_k‖}_M}.$$E112

The estimate of the total error follows from the inequalities

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{{‖x-{\overline{x}}_k‖}_N}{{‖x‖}_N}+\frac{{‖{\overline{x}}_k-\overline{\overline{x}}‖}_N}{{‖x‖}_N},\frac{{‖{\overline{x}}_k-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{{‖A‖}_{\mathit{MN}}{‖{{\overline{A}}_k}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖{\overline{r}}_k‖}_M}{{‖b_k‖}_M}$$E113

and estimates (25), (64).

Case 3. The rank of the original matrix is larger than that of the perturbed matrix, i.e., $rank\left(A\right)>rank\left(\overline{A}\right)=l$.

Consider the case when the condition ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$ not satisfied and the rank of the perturbed matrix can decrease.

Theorem 19. Assume that $rank\left(A\right)>rank\left(\overline{A}\right)=l$, $\frac{{‖\mathit{\Delta A}‖}_{\mathit{MN}}}{{\mu }_l}<\frac{1}{2}$ and let $\overline{x}\in \mathit{Im}\left({\overline{A}}^{\#}\right)$. Then

$$\frac{{‖x_l-\overline{\overline{x}}‖}_N}{{‖x_l‖}_N}\le \frac{{\mu }_1/{\mu }_l}{1-2{‖\mathit{\Delta A}‖}_{\mathit{MN}}/{\mu }_l}\left(2{\epsilon }_A+{\alpha }_l+\frac{{\mu }_1}{{\mu }_l}{\epsilon }_A{\beta }_l+{\gamma }_l\right)$$E114

Proof. For the proof, along with problem (9), consider the problem

$$\underset{x\in C}{min}{‖x_l‖}_N,C=\left\{x|{‖A_{l}x-b‖}_M=min\right\}$$E115

With matrix $A_l=U{\Sigma }_l{V}^T$ with rang l.

The estimate of the computational error in this case will be

$$\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖\overline{x}‖}_N}\le {‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{{‖\overline{r}‖}_M}{{‖{\overline{b}}_l‖}_M}.$$E116

The estimate of the total error follows from the inequalities

$$\frac{{‖x_l-\overline{\overline{x}}‖}_N}{{‖x_l‖}_N}\le \frac{{‖x_l-\overline{x}‖}_N}{{‖x_l‖}_N}+\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖x_l‖}_N},\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{{‖A‖}_{\mathit{MN}}{‖{{\overline{A}}_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}{‖{\overline{r}}_l‖}_M}{{‖b_l‖}_M},$$E117

obvious relationships ${‖A‖}_{\mathit{MN}}={‖A_l‖}_{\mathit{MN}},{‖{A_l}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mu }_l$}\right.,$ estimates of the hereditary error (75) and the inequality ${‖\Delta A_l‖}_{\mathit{MN}}\le 2{‖\mathit{\Delta A}‖}_{\mathit{MN}}$.

5.2 Estimates of the total error of the weighted normal pseudosolution for matrices of full rank

In the following Theorems 20 and 21, the weighted pseudoinverse $A_{\mathit{MN}}^{+}$ is represented in accordance with the properties of the full rank matrix (81).

Theorem 20. Let ${‖\mathit{\Delta A}‖}_{\mathit{MI}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{IM}}<1$, $m>n=k$ and $\overline{\overline{x}}\in R\left({\overline{A}}^{\#}\right).$ Then

$$\frac{‖x-\overline{\overline{x}}‖}{‖x‖}\le \frac{h}{1-h{\epsilon }_{\mathrm{А}}}\left({\epsilon }_{\mathrm{А}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖A‖}_{\mathit{MI}}‖x‖}+h{\epsilon }_{\mathrm{А}}\frac{{‖r‖}_M}{{‖A‖}_{\mathit{MI}}‖x‖}+\frac{{‖{\overline{r}}_k‖}_M}{{‖A‖}_{\mathit{MI}}‖x‖}\right).$$E118

Proof. The estimate of the computational error is determined by formula (106), namely

$$\frac{‖\overline{x}-\overline{\overline{x}}‖}{‖\overline{x}‖}\le {‖\overline{A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}\frac{{‖\overline{r}‖}_M}{{‖{\overline{b}}_k‖}_M}.$$E119

The estimate for the total error (118) follows from the inequalities

$$\frac{‖x-\overline{\overline{x}}‖}{‖x‖}\le \frac{‖x-\overline{x}‖}{‖x‖}+\frac{‖\overline{x}-\overline{\overline{x}}‖}{‖x‖},\frac{‖\overline{x}-\overline{\overline{x}}‖}{‖x‖}\le {‖A‖}_{\mathit{MI}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{IM}}\frac{{‖\overline{r}‖}_M}{{‖b_k‖}_M}$$E120

and estimates for the pseudoinverse matrix (25) and the hereditary error (83).

Theorem 21. Let ${‖\mathit{\Delta A}‖}_{\mathit{\text{IN}}}{‖A_{\mathit{MN}}^{+}‖}_{\mathit{NI}}<1$, $n>m=k$ and $\overline{\overline{x}}\in R\left({\overline{A}}^{\#}\right)$. Then

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖x‖}_N}\le \frac{h}{1-h{\epsilon }_{\mathrm{А}}}\left(2{\epsilon }_A+\frac{‖\mathrm{\Delta }b‖}{{‖A‖}_{\mathit{\text{IN}}}{‖x‖}_N}+\frac{‖\overline{r}‖}{‖b_k‖}\right),$$E121

The proof of Theorem 21 is similar to the proof of the previous theorem, taking into account the estimate for the hereditary error (86).

Remark 7. Here, we did not indicate a method for obtaining an approximate weighted normal pseudosolution $\overline{\overline{x}}$, satisfying the conditions of the theorems. Algorithms for obtaining such approximations are considered, for example, in Section 4.2.

Remark 8. Along with estimates (103), (109), (114), (118), (121), error estimates can be obtained, the right-hand sides of which depend on the input data of systems of linear algebraic equations with approximately given initial data. For example, the following theorem holds.

Theorem 22. Let ${‖\mathit{\Delta A}‖}_{\mathit{MN}}{‖{\overline{A}}_{\mathit{MN}}^{+}‖}_{\mathit{NM}}<1$, and$\overline{\overline{x}}\in R\left({\overline{A}}^{\#}\right)$. Then, for the total error of the normal pseudosolution, the following estimate is fulfilled

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖\overline{x}‖}_N}\le \frac{h\left(\overline{A}\right)}{1-h\left(\overline{A}\right){\epsilon }_{\overline{\mathrm{А}}}}\left(2{\epsilon }_{\overline{\mathrm{А}}}+\frac{{‖\mathrm{\Delta }b‖}_M}{{‖\overline{A}‖}_{\mathit{MN}}{‖\overline{x}‖}_N}+h\left(\overline{A}\right){\epsilon }_{\overline{A}}\frac{{‖\overline{r}‖}_M}{{‖\overline{A}‖}_{\mathit{MN}}{‖\overline{x}‖}_N}+\frac{{‖\overline{r}‖}_M}{{‖{\overline{b}}_k‖}_M}\right).$$E122

Estimate (122) can be obtained from the inequality

$$\frac{{‖x-\overline{\overline{x}}‖}_N}{{‖\overline{x}‖}_N}\le \frac{{‖x-\overline{x}‖}_N}{{‖\overline{x}‖}_N}+\frac{{‖\overline{x}-\overline{\overline{x}}‖}_N}{{‖\overline{x}‖}_N}$$E123

and estimates (90), (106).

If the weighted pseudoinverse matrix is known or its weighted singular value decomposition is obtained during the process of solving the problem, then a practical estimate of the computational error can be obtained using (104). When calculating the residual $\overline{r}={\overline{b}}_k-\overline{А}\overline{\overline{x}}$, the explicit form of the projection operator onto $R\left({\overline{А}}^{\#}\right)$ is used.

In conclusion, we note that the determining factor for obtaining estimates is the use of a weighted singular value decomposition [37] and the technique of reducing the problem of estimating the error of a pseudosolution to an estimate of the error [32] for problems with matrices of the same rank. Based on the results obtained, an algorithm for finding the effective rank of matrices can be developed, as well as an algorithm for calculating stable projections of a weighted normal pseudosolution.

5.3 Software-algorithmic methods for increasing the accuracy of computer solutions

The numerical methods we have considered for solving systems of linear algebraic equations and WLS problems have one common property. Namely, the actually calculated solution (pseudosolution) is exact in accordance with the inverse analysis of errors [43] for some perturbed problem. These perturbations are very small and are often commensurate with the rounding errors of the input data. If the input data is given with an error (measurements, calculations, etc.), then usually they already contain significantly larger errors than rounding errors. In this case, any attempt to improve the machine solution (pseudosolution) without involving additional information about the exact problem or errors of the input data errors will be untenable.

The situation changes significantly if a mathematical problem with accurate input data is considered. Now the criterion of bad or good conditionality of the computer model of the problem depends on the mathematical properties of the computer model of the problem and the mathematical properties of the processor (length of the computer word), and it becomes possible in principle to achieve any given accuracy of the computer solution. In this case, as follows from estimates (48), (64), (75), (83), (86), it is obviously possible to refine the computer solution by solving a system with increased bit depth, in particular, using the GMP library [44] for implementation of computations with arbitrary bit depth.

To predict the length of the mantissa (machine word) that provides a given accuracy for a solution (joint systems), you can use the following rule of thumb: the number of correct decimal significant digits in a computer solution is $\mu -\alpha$, where μ is the decimal order of the mantissa of a floating-point number ε, α is the decimal order of the condition number. Thus, knowing the conditionality of the matrix of the system and the accuracy of calculations on a computer, it is possible to determine the required bit depth to obtain a reliable solution.

The GMP library is used to work on integers, rational numbers and floating-point numbers. The main feature of the library is the bitness of numbers (precision) is practically unlimited. Therefore, the main field of application is computer algebraic calculations, cryptography, etc. The functions of the GMP library allow not only setting the bit depth at the beginning of the program and performing calculations with this bit depth, but also changing the bit width as needed in the computation process, i.e. execute different fragments of the algorithm with different bit depths.

The library’s capabilities were tested in the study of solutions to degenerate and ill-conditioned systems in [45].

6. Conclusions

In the framework of these studies, estimates of the hereditary error of the weighted normal pseudosolution for matrices of arbitrary form and rank are obtained, including when the rank of the perturbed matrix may change. Three cases are considered: the rank of the matrix does not change when the data is disturbed, the rank increases and the rank decreases. In the first case, the weighted normal pseudosolution of the approximate problem is taken as an approximation to the weighted normal pseudosolution, in the other two, the problem is reduced to the case when the ranks of the matrices are the same. Also, the estimates of the error for the weighted pseudoinverse matrix and the weighted condition number of the matrix are obtained, the existence and uniqueness of the weighted normal pseudosolution are investigated and proved. Estimates of the total error of solving the weighted least squares problem with matrices of arbitrary form and rank are established.

The results obtained in the perturbation theory of weighted least squares problem can be a theoretical basis for further research into various aspects of the WLS problem and the development of methods for calculating weighted pseudoinverse matrices and weighted normal pseudosolutions with approximate initial data, in particular, in the design and optimization of building structures, in tomography, in the calibration of viscometers, in statistics. The results of the research can be used in the educational process when reading special courses on this section of the theory of matrices.