Open access peer-reviewed chapter - ONLINE FIRST

Weighted Least Squares Perturbation Theory

Written By

Aleksandr N. Khimich, Elena A. Nikolaevskaya and Igor A. Baranov

Submitted: January 14th, 2022Reviewed: January 26th, 2022Published: April 16th, 2022

DOI: 10.5772/intechopen.102885

IntechOpen
Matrix Theory - Classics and AdvancesEdited by Mykhaylo Andriychuk

From the Edited Volume

Matrix Theory - Classics and Advances [Working Title]

Dr. Mykhaylo I. Andriychuk

Chapter metrics overview

14 Chapter Downloads

View Full Metrics

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.

Keywords

  • weighted least squares problem
  • error estimates
  • weighted matrix pseudoinverse
  • weighted condition number
  • weighted singular value decomposition

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.

Advertisement

2. Weighted least squares problem

2.1 Preliminaries

Let the set of all m×nmatrices is denoted by Rm×n. Given a matrix ARm×nlet ATis the transpose of A, rankAis the rank of A, RAis the field of values of Aand NAis the null space of A. Additionally, let denote the vector 2-norm and the consistent matrix 2-norm, and let Ibe an identity matrix.

Given an arbitrary matrix ARm×nand symmetric positive definite matrices Mand Nof orders mand n, respectively, a unique matrix XRm×n, satisfying the conditions:

AXA=A,XAX=X,MAXT=MAX,NXAT=NXA,E1

is called the weighted Moore–Penrose pseudoinverseof Aand is denoted by X=AMN+. Specifically, if M=IRm×mand N=IRn×n, then Xsatisfying conditions (1) is called the Moore–Penrose pseudoinverseand is designated as X=A+.

Let A#denote the weighted transpose of A, Pand Qbe idempotent matrices, and А¯=A+ΔAbe a perturbed matrix, i.e.,

A#=N1ATM.E2
P=AMN+A,Q=AAMN+,P¯=A¯MN+A¯,Q¯=A¯A¯MN+.E3

Let xRm, yRn. The weighted scalar products in Rmand Rnare defined as xyM=yTMx, x,yRmand xyN=yTNx, x,yRn, respectively. The weighted vector norms are defined as:

xM=xxM12=xTMx12=M12x,xRm,yN=yyN12=yTNy12=N12y,yRn.E4

Let x,yRmand xyM=0. Then the vectors xand yare called M-orthogonal, i.е. M12x- and M12y-orthogonal. It is easy to show that.

x+yM2=xM2+yM2,x,yRm.E5

The weighted matrix norms are defined as:

AMN=maxxN=1AxM=M12AN12,ARm×n,BNM=maxyM=1ByN=N12AM12,BRn×m.E6

Lemma 1(see in [37]). Let ARm×n, rankA=k, Mand Nare positive definite matrices of orders mand n, respectively. Then, there are matrices URm×mand VRn×n, satisfying UTMU=Iand VTN1V=Isuch that.

A=UD000VT,AMN+=N1VD1000UTM,E7

where D=diagμ1μ2μk, μ1μ2μk>0and μi2are the nonzero eigenvalues of the matrix A#A. The nonnegative values μiare called the weighted singular values of A, moreover, AMN=μ1, AMN+NM=1μk.

The weighted singular value decomposition of Ayields an M-orthonormal basis of the vectors of Uand an N1-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

Ax=b.E8

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

minxCxN,C=xAxbM=min,E9

where ARm×nis a rank-deficient matrix and bRm.

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

minxCx¯N,C=x¯A+ΔAx¯b+ΔbM=min,E10

where.

А¯=A+ΔA,b¯=b+Δb,x¯=x+Δx.E11

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

ΔAMNεAAMN,ΔbMεbbM.E12

The problem for the approximate solution x¯¯of a system of linear algebraic equations with approximately given initial data

A¯x¯¯=b¯+r¯,E13

where r¯=A¯x¯¯b¯is the residual vector.

The analysis of the reliability of the obtained solution includes an assessment of the hereditary error xx¯N, computational error x¯x¯¯Nand total error xx¯¯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 Lbe a nonempty subset of space R, closed with respect to the operations of addition and multiplication by a scalar (if xand yare elements of Lα,β, the αx+βyis an element of L). Vector xis N-orthogonal to the linear manifold L(xNL) if xis N-orthogonal to each vector from L.

Lemma 2(see in [38]). There exists a unique decomposition of vector x, namely x=x̂+x, where x̂L, xNL.

Let Ais an arbitrary matrix. The kernel of matrix A, denoted by N(A), is the set of vectors mapped into zero by А: NA=x:Ax=0.

The set RAof images of matrix Ais the set of vectors that are images of vectors of the space Rfrom the definition domain of A, i.e. RA=b:b=Axx.

Let Lbe a linear manifold in space R, N-orthogonal (M-orthogonal) complement to L, denoted by LN(LM), defined as the set of vectors in R, each of which is N-orthogonal (M-orthogonal) to L.

Remark 1. If xis a vector from Rand xTNy=0for any yfrom R, then x=0.

Theorem 1.Let ARm×n, then NA=RA#.

Proof. Vector xNA, if and only if Ax=0. Hence, by virtue of Remark1, we get xNA, if and only if yTMAx=0for any y. Since yTMAx=A#yTNx, we get Ax=0if and only if xis N-orthogonal to all the vectors of the formA#y. VectorsA#yform RA#. The required statement follows from here and from the definition RMA.

Theorem 2(see in [38]). If Ais an m×nmatrix and bis an m-dimensional vector, then the unique decomposition b=b̂+bholds, where b̂RAand bNA#.

Vectorb̂is a projection of bon to RA, and bis a projection of bon to NA#. Vectorsb̂and bare M-orthogonal. Hence, A#b=A#b̂.

By Theorem 1, the following relations hold for the symmetric matrix A: NA=RA,RA=NA.

Theorem 3.Let ARm×n, then RA=RAA#, RA#=RA#A, NA=NA#Aand NA#=NAA#.

Proof. It will be to establish that NA#=NAA#and NA=NA#A.

For this purpose, we will use Theorem 1. To prove the coincidence of NA#and NAA#, note that AA#x=0if A#x=0. On the other hand, if AA#x=0, then xTAA#x=0, i.е. A#xM=0, which entails equality A#x=0. So, A#x=0if and only if xTAA#x=0. We can similarly establish thatNA=NA#A.

Then let us prove the theorem about the existence and uniqueness of the solution vector that minimizes the norm of the residual AxbMby the technique proposed in [39] for the least-squares problem.

Theorem 4.Let ARm×n, bRm, bRA. Then there exists a vector x̂, that minimizes the norm of the residual AxbMand vector x̂is a unique vector from RA#, that satisfies the equation Ax=b̂, where b̂=AAMN+bis the projection of bontoRA.

Proof. By virtue of Theorem 2, we get b=b̂+b, where b=IAAMN+bis the projection of bon to NA#. Since for every x,AxRAand bRMA, then b̂AxRAand bb̂Ax. Therefore

bAxM2=b̂Ax+bM2=b̂AxM2+bM2bM2.E14

This lower bound is attained since b̂belongs to the set of images А, i.е. b̂is an image of some x0: b̂=Ax0.

Thereby, for this x0the greatest lower bound is attainable:

bAx0M2=bb̂M2=bM2.E15

It was shown earlier that

bAxM2=b̂AxM2+bM2E16

and hence, the lower bound can only be attained forx, for which Ax=b̂. According to Theorem 2, each vector x, can be presented as a sum of two orthogonal vectors: x=x̂+x, where x̂RA#, xNA.

Therefore, Ax=Ax̂and hence, bAxM2=bAx̂M2. Note that

xN2=x̂N2+xN2x̂N2,E17

where strict inequality is possible whenxx̂(i.е. if xdoes not coincide with its projection on to RA#).

It was shown above, that x0minimizes AxbM, if and only if Ax0=b̂, and among the vectors that minimize AxbM, 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 x̂andxbelong to RA#and that Ax̂=Ax=b̂. Thenxx̂RA#, howeverAxx̂=0, so xx̂NA=RNA#.

As vectorxx̂is N-orthogonal to itselfxx̂N=0, i.е. x=x̂.

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

bAA#yM=infxbAxM.E18

If

bAx0M=infxbAxM,E19

then x0NA#yNwith strict inequality forx0A#yN.

Vector ysatisfies the equation AA#y=b̂, here b̂is the projection of bonto RA.

Theorem 5. Among all the vectors xthat minimize the residual AxbM, vector x̂, which has the minimum norm x̂=minxN, is a unique vector of the form

x̂=N1ATMy=A#y,E20

satisfying the equation

A#Ax=A#b,E21

i.е. x̂can be obtained by means of any vector y0, that satisfies the equation A#AA#y=A#bby the formula x̂=A#y0.

Proof. According to the condition of Theorem 3 RA#=RA#A. Since vector A#bbelongs to the set of images A#, it should belong to the set of imagesA#Aand thus should be an image of some vector xwith respect to the transformation A#A. In other words, Eq. (21) (with respect to x) has at least one solution. If xis a solution of Eq. (21), then x̂is the projection of xon to RA#, since Ax=Ax̂according to Theorem 2. Since x̂RA#, vector x̂is an image of some vector ywith respect to the transformation A#: 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 x̂1=A#y1and x̂2=A#y2satisfy Eq. (21). Then A#AA#y1A#y2=0, therefore, A#y1y2NA#A=NA, where from the equality AA#y1y2=0follows.

Therefore y1y2NAA#=NA#; hence,x̂1=A#y1=A#y2=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 Ax=b̂, where b̂is a weighted projection of bon to RA, i.е. A#b=A#b̂.

Theorem 4 establishes that there is a unique, from RA#solution of the equation

Ax=b̂.E22

Hence, this unique solution satisfies the equation A#Ax=A#b̂.

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

AMN+=A#A#AA#+A#,E23

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

Theorem 6.Let ARm×n, then x=AMN+b—is М-weighted least squares solution with the minimum N-norm of the system Ax=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 Ax=b.

Advertisement

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,ΔARm×n, μiAandμiA¯denote the weighted singular values of Aand A¯respectively. Then,

μiAΔAMNμiA¯μiA+ΔAMN.E24

Lemma 4(see [40]). Let A,ΔARm×n, rankA¯=rankAand ΔAMNAMN+NM<1. Then

A¯MN+NMAMN+NM1ΔAMNAMN+NM.E25

Lemma 5. Let G=A¯MN+AMN+, А¯=А+ΔAand rankA¯=rankA. Then Gcan be represented as the sum of three matrices G=G1+G2+G3, where

G1=A¯MN+ΔAAMN+,E26
G2=IP¯N1ΔATAMN+TNAMN+=IP¯ΔA#AMN+#AMN+,E27
G3=A¯MN+IQ.E28

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

G=P¯+IP¯A¯MN+AMN+Q+IQ==P¯A¯MN+Q+P¯A¯MN+IQP¯AMN+QP¯AMN+IQ+IP¯A¯MN+Q++IP¯A¯MN+IQIP¯AMN+Q+IP¯AMN+IQ.E29

Since.

P¯A¯MN+=A¯MN+,IP¯A¯MN+=0,AMN+Q=AMN+,AMN+IQ=0,E30

we obtain

G=A¯MN+Q+A¯MN+IQP¯AMN++IP¯A¯MN+==A¯MN+QP¯AMN+IP¯AMN++A¯MN+IQ.E31

Consider each term in this equality separately

G1=A¯MN+QP¯AMN+=A¯MN+AAMN+A¯MN+A¯AMN+=A¯MN+AA¯AMN+=A¯MN+ΔAAMN+.E32

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

AMN+=AMN+AAMN+=N1NAMN+ATAMN+==N1ATAMN+TNAMN+=N1A¯TAMN+TNAMN+N1ΔATAMN+TNAMN+.E33

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

G2=IP¯AMN+=IP¯N1A¯TAMN+TNAMN+N1ΔATAMN+TNAMN+.E34

Since,

IP¯N1A¯TAMN+TNAMN+=N1A¯TAMN+TNAMN+A¯MN+A¯N1A¯TAMN+TNAMN+==N1A¯TAMN+TNAMN+N1A¯TAMN+TNAMN+=0E35

we obtain

G2=IP¯AMN+=IP¯N1ΔATAMN+TNAMN+=IP¯ΔA#AMN+#AMN+.E36

Finally,

G=A¯MN+AMN+=A¯MN+ΔAAMN+IP¯ΔA#AMN+#AMN++A¯MN+IQ.E37

Lemma 6(see in [41]). If rankA¯=rankA=k, then

Q¯IQMM=QIQ¯MM,E38

where Qand Q¯are defined in (3).

Lemma 7. Let A,ΔARm×n, rankA¯=rankAand ΔAMNAMN+NM<1.

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

A¯MN+AMN+NMAMN+NMCεAh1εAh,E39

where h=hA=AMNAMN+NMand the estimate of the absolute error

A¯MN+AMN+NMCεAh21εAh,E40

moreover.

if Ais not a full rank matrix, then C=3,

if m>n=kor n>m=k, then C=2,

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

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

A¯MN+AMN+=A¯MN+ΔAAMN+IP¯ΔA#AMN+#AMN++A¯MN+IQ.E41

Passing to the weighted norms, we obtain.

A¯MN+AMN+NMA¯MN+ΔAAMN+NM+ΔA#AMN+#AMN+NM+A¯MN+Q¯IQNM.E42

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

A¯MN+Q¯IQN=A¯MN+A¯A¯MN+IQNA¯MN+NMQ¯IQMM=A¯MN+NMQIQ¯MM.E43

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

A¯MN+AMN+NMA¯MN+ΔAAMN+NM+ΔA#AMN+#AMN+NM+AMN+QIQ¯NMA¯MN+NMΔAMNAMN+NM+ΔAMNAMN+NM2+A¯MN+NMAMN+NMΔAMN.E44

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

A¯MN+AMN+NMAMN+NM1ΔAMNAMN+NM(ΔAMNAMN+NM+ΔAMNAMN+NM++AMN+NMΔAMN)=hεA1hεAhεA+hεA+hεA=ChεA21hεA,C=1,2,3.E45

To estimate the relative error, we have

A¯MN+AMN+NMAMN+NM3ΔAMNAMN+NM11ΔAMNAMN+NM=ChεA1hε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:

α=ΔbMAMNxN,β=rMxNAMN,γ=r¯MAMNxNαl=ΔbMAMNxlN,βl=rMxlNAMN,γl=r¯lMAMNxlN,γk=r¯kMAMNxN.E47

Consider the following three cases.

Case 1. The rank of the original matrix A remains the same under its perturbation, i.e.,rankA=rankA¯.

Theorem 7. Assume that ΔAMNAMN+NM<1, rankA¯=rankA. Then

xx¯NxNh1hεA2εA+α+hεAβ,E48

where hA=AMNAMN+NMis the weighted condition number of A, the symbolsMNandNMdenote the weighted matrix norms defined by Eq. (4)(6), and AMN+is the weighted Moore–Penrose pseudoinverse.

Proof. The error estimate follows from the relation:

xx¯=АMN+А¯MN+b+А¯MN+bb¯.E49

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

A¯MN+AMN+=A¯MN+ΔAAMN+IP¯N1ΔATAMN+TNAMN++A¯MN+IQ.E50

Then,

xx¯=A¯MN+ΔAAMN++IP¯N1ΔATAMN+TNAMN+A¯MN+IQb+A¯MN+bb¯==A¯MN+ΔAAMN+b+IP¯N1ΔATAMN+TNAMN+bA¯MN+IQb+A¯MN+bb¯
=A¯MN+ΔAx+IP¯N1ΔATAMN+TNxA¯MN+IQb+A¯MN+bb¯E51

Thus,

x¯x=A¯MN+ΔAx+IP¯N1ΔATAMN+TNxA¯MN+IQb+A¯MN+bb¯.E52

Passing to the weighted norms yields

x¯xN=A¯MNΔAx+IP¯N1ΔATAMN+TNxA¯MN+IQb+A¯MN+bb¯NA¯MNΔAxN+IP¯N1ΔATAMN+TNxN+A¯MN+IQb+A¯MN+bb¯N.E53

By taking into account the relations

IQb=IQr=r,r=bAx,x=AMN+bE54

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

  1. A¯MN+ΔAxN=N12A¯MN+M12M12ΔAN12N12xN12A¯MN+M12M12ΔAN12N12x=A¯MN+NMΔAMNxN.(55)

  2. IP¯N1ΔATAMN+TNxN=N12IP¯N12N12ΔATM12M12AMN+TN12N12xN12IP¯N12M12ΔAN12N12AMN+M12N12x=IP¯NNΔAMNAMN+NMxN(56)

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

A¯MN+Q¯IQbN=A¯MN+A¯A¯MN+IQrNA¯MN+NMQ¯IQMMrM=A¯MN+NMQIQ¯MMrME57

where

QIQ¯MM=AAMN+IQ¯MM=M12AAMN+IQ¯M12==M12MAAMN+TIQ¯M12==M12AMN+TATA¯TM12M12IQ¯M12==M12AMN+TΔATMIQ¯M12M12AMN+TΔATM12==M12ΔAN12N12AMN+M12M12ΔAN12N12AMN+M12ΔAMNAMN+NM.E58

Substituting this result into (31) gives the inequality

A¯MN+Q¯IQbNA¯MN+NMΔAMNAMN+NMrM.E59

  1. A¯MN+bb¯N=N12A¯MN+M12M12bb¯N12A¯MN+M12M12bb¯=A¯MN+NMbb¯M(60)

Taking into account IP¯<1, and applying Lemma 4, we obtain the following weighted-norm estimate for the relative error

xx¯NxNA¯MN+NMΔAMNxNxN+ΔAMNAMN+NMxNxN++A¯MN+NMAMN+NMΔAMNrMxN+A¯MN+NMΔbMxNA¯MN+NMΔAMN+ΔAMNAMN+NM++A¯MN+NMAMN+NMΔAMNrMxN+A¯MN+NMΔbMxNAMN+NMAMN1ΔAMNAMN+NM2ΔAMNAMN+ΔbMAMNxN++AMN+NMAMNΔAMNAMNrMAMNxNhA1hAεA2εA+ΔbMAMNxN+hAεArMxNAMN.E61

as required.

Specifically, if M=IRm×mand N=IRn×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║ΔА║║А+║< 1, rankA¯=rankA=k. Then

xx¯xh1hεА2εА+εbk+hεАbbkbk,E62

where bkis the projection of the right-hand side of problem (8) onto the principal left singular subspace of the matrix A[42], i.е., bkImA, h=hA=AA+is condition number of A, the symbol , 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.rankA¯>rankA=k.

Define the idempotent matrices:

P=AMN+A,Q=AAMN+,P¯k=A¯kMN+A¯,Q¯k=A¯A¯kMN+,E63

where kis the rank of A.

Theorem 9. Assume that ΔAMNAMN+NM<12, rankA¯>rankA=k. Then

xx¯kNxNh12hεА2εА+α+hεАβ.E64

where hA=AMNAMN+NMis the weighted condition number of A, the symbols MNandNMdenote the weighted matrix norms defined Eq. (4)(6), and AMN+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, А¯is represented as a weighted singular value decomposition:

А¯=U¯D¯V¯T.E65

Along with (38), we consider the decomposition

А¯k=U¯D¯kV¯T,E66

where D¯kis a rectangular matrix whose first kdiagonal elements are nonzero and equal to the corresponding elements of 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 x¯kto the problem

minxCx¯N,C=xA¯kx¯b¯M=min.E67

The matrix A¯kis defined by (48) and has the same rank kas 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 xx¯kN/xN. The error of the weighted pseudoinverse matrix then becomes:

Gk=P¯k+IP¯kA¯kMN+AMN+Q+IQ=P¯kA¯kMN+Q+P¯kA¯kMN+IQP¯kAMN+QP¯kAMN+IQIP¯kA¯kMN+Q+IP¯kA¯kMN+IQIP¯kAMN+Q+IP¯kAMN+IQ=A¯kMN+QP¯kAMN+IP¯kAMN+++A¯kMN+IQ=A¯kMN+AAMN+A¯kMN+A¯AMN+IP¯kAMN++A¯kMN+IQ==A¯kMN+AA¯AMN+IP¯kAMN++A¯kMN+IQ,E68

Applying Lemma 5 yields

Gk=A¯kMN+AMN+=A¯kMN+ΔAAMN+IP¯kN1ΔATAMN+TNAMN++A¯kMN+IQk.E69

For the error of the WLS solution, we obtain

x¯kx=A¯kMN+ΔAx+IP¯kN1ΔATAMN+TNxA¯kMN+IQkb+A¯kMN+bb¯.E70

Passing to the weighted norms and applying Lemma 4 gives

xx¯kNxNA¯kMN+NMΔAMNxNxN+ΔAMNAMN+NMxNxN++A¯kMN+NMAMN+NMΔAMNrMxN+A¯kMN+NMΔbMxNA¯kMN+NMΔAMN+ΔAMNAMN+NM++A¯kMN+NMAMN+NMΔAMNrMxN+A¯kMN+NMΔbMxNAMN+NMAMN1ΔAkMNAMN+NM2ΔAMNAMN+ΔbMAMNxN++AMN+NMAMNΔAMNAMNrMAMNxN.E71

Let estimate ΔAk=AA¯k:

ΔAkMN=A¯kAMN=A¯kA¯+ΔAMNA¯kA¯MN+ΔAMN==U¯000Dk+1V¯TMN+ΔAMN2ΔAMN.E72

Moreover, the theorem condition ΔAMNAMN+NM<12leads to ΔAkMNAMN+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=IRm×mand N=IRn×n, then the estimates of the hereditary error of normal pseudosolutions of systems of linear algebraic equations for case rankA¯>rankA=kfollows from next theorem.

Theorem 10(see in [32]). Let ΔAA+<12, rankA¯>rankA=k. Then

xx¯kxh12hεА2εА+εbk+hεАbbkbk.E73

Case 3. The rank of the original matrix is larger than that of the perturbed matrix, i.e.,rankA>rankA¯=l.

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

Pl=AlMN+A,Ql=AAlMN+,P¯=A¯MN+A¯,Q¯=A¯A¯MN+,E74

Theorem 11. Assume that rankA>rankA¯=l, ΔAMNμl<12. Then,

xlx¯NxlNμ1/μl12ΔAMN/μl2εA+αl+μ1μlεАβl,E75

where μiare the weighted singular values of A.

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

minxCxlN,C=xAlxbM=minE76

with the matrix Al=UDlVTof rank l.

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

Gl=A¯MN+AlMN+=A¯MN+ΔAAlMN+IP¯N1ΔATAlMN+TNAlMN++A¯MN+IQl,E77
x¯xl=A¯MN+ΔAx+IP¯N1ΔATAlMN+TNxA¯MN+IQlb+A¯MN+bb¯.E78

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

xlx¯NxNA¯MN+NMΔAMN+ΔAMNAlMN+NM++A¯MN+NMAlMN+NMΔAMNrMxN+A¯MN+NMΔbMxNAlMN+NMAMN1ΔAlMNAlMN+NM2ΔAMNAMN+ΔbMAMNxN++AlMN+NMAMNΔAMNAMNrMAMNxN,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=IRm×mand N=IRn×n, then the estimates of the hereditary error of normal pseudosolutions of systems of linear algebraic equations for case rankA>rankA¯=lfollows from next theorem.

Theorem 12(see in [32]). Let rankA>rankA¯=l, ΔAμl<12. Then

xlx¯xlμ1/μl12A/μl2εA+εbl+εАμ1μlbblbl,E80

where xlis the projection of the normal pseudosolution of problem (8) onto the right principal singular subspace of the matrix Aof dimension l, blis projection of the right-hand side bonto the principal left singular subspace of dimension lof the matrix A, μiis 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 ΔAMNAMN+NM<1is met.

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

AMN+=ATMA1ATMformnandAMN+=N1ATAN1AT1fornm.E81

If mn, then problem (9) is reduced to a problem of the form

minxRnAxbM.E82

For such a problem, the following theorem is true.

Theorem 13. Let ΔAMIAMN+IM<1, m>n=k. Then

xx¯xh1hεAεA+ΔbMAMIx+hεArMxAMI,E83

where h=AMIAMN+IM.

Proof. To prove Theorem 13, as before, we will use relation (49). By (81)P¯=A¯MN+A¯=I, so that from (50) we have the equality

A¯MN+AMN+=A¯MN+ΔAAMN++A¯MN+IQ,E84

using which we obtain (83).

If nm, then problem (9) is reduced to a problem of the form

minxCxN,C=xAx=bE85

and the following theorem holds for it.

Theorem 14. Let ΔAINAMN+NI<1, n>m=k. Then

xx¯NxNh1hεA2εA+ΔbAINxN,E86

where h=AINAMN+NI.

Proof. Since in this case Q=AAMN+=I, then the expression for A¯MN+AMN+by (81) takes the form

A¯MN+AMN+=A¯MN+ΔAAMN+IP¯N1ΔATAMN+TNAMN+.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 hA¯is established by the estimates

σkΔAMNσ¯kσk+ΔAMN,σ1ΔAMNσ¯1σ1+ΔAMN,σ1ΔAMNσk+ΔAMNσ¯1σ¯kσ1+ΔAMNσkΔAMN,1εA1+εAhhАhА¯1+εA1εAh,E88

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

Lemma 7. Let A,ΔARm×n, rankA¯=rankAand ΔAMNAMN+NM<1. Then the estimate of the relative error of the condition number of the matrix Ahas the form

h¯hhεA1+h1εAhE89

where h=hA=AMNAMN+NMis weighted condition number of matrix A, h¯=h¯A=A¯MNA¯MN+NMis weighted condition number of the perturbed matrix A¯=A+ΔA.

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

Theorem 15. Let ΔAMNA¯MN+NM<1, ΔAMNεA¯A¯MN, rankA¯=rankA. Then,

x¯xNx¯NhA¯1hA¯εA¯2εA¯+ΔbMA¯MNx¯N+hA¯εA¯r¯Mx¯NA¯MN,E90

where hA¯=A¯MNA¯MN+NMis weighted matrix condition number A¯, the symbols MNand NMdenote the weighted matrix norms defined by Eq. (4)(6) and AMN+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

xx¯NxNxx¯Nx¯N1+xx¯NxNE91

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

Theorem 16.Let ΔAMNA¯MN+NM<1, ΔAMNεA¯A¯MN, rankA¯=rankA. Then

x¯xNxNβ1β,β=hA¯1hA¯εA¯2εA¯+ΔbMA¯MNx¯N+hA¯εA¯r¯Mx¯NA¯MN.E92

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

Advertisement

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 ΔA, Δ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 ΔAchange 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 rankA=minmn, and

ΔAMNAMN+NM<1,E93

Then rankA¯=rankA.

Proof. For proof, let, for example, rankA=m. Taking equalΔAMN=ε, in equality (93) can be rewritten as εμm<1, which is equivalent

μmε>0.E94

Let μ¯mm-weighted singular value of perturbed matrix A¯=A+ΔA. According to Lemma 3, we can write μ¯mμmε. Then, taking into account (94), we obtain μ¯mμmε>0.

Therefore rankA¯m, whence we come to the conclusion that rankA¯=m, i.е. rankA¯=rankA.

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

εA¯hA¯<1,E95
1.0+1hA¯1.0E96

where hA¯=A¯MNA¯MN+NMis weighted condition number of matrix 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 hA¯. 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/hA¯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 δ-rank) is

rankAδ=minABMNδrankB.E97

This means that the δ-rank of the matrix is equal to the minimum rank among all matrices in the neighborhood ABMNδ.

From [28] that if rδis the δ-rank of the matrix, then

μ1μrδ>δμrδ+1μp,p=minmn.E98

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

δμi<1,μi0,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 δcan be attributed to machine precision, for example, setting it equal machepsB.

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 ARm×nand rankA=k, M-and N-positive-defined matrices of order mand 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 Aand A¯is reduced to the following three cases.

  1. If the rank of the matrix has not changed rankA¯=rankA=k, an approximate weighted normal pseudosolution is constructed by the formula

    x¯=A¯MN+b¯,E100

where A¯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 ΔAMNAMN+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).

  • Matrix rank increased rankA¯>rankA=k. An approximate weighted normal pseudosolution is constructed by the formula

    x¯k=A¯kMN+b¯.E101

    Weighted pseudoinverse matrix A¯kMN+is defined as follows

    A¯kMN+=N1V¯D¯k+U¯TM,E102

    where D¯kis a rectangular matrix, the first kdiagonal elements of which are nonzero and coincide with the corresponding elements of the matrix 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 kof the matrix A¯and, if ΔAMNAMN+NM<12, then the error of the solution is estimated by formula (64).

  • If the rank of the matrix has decreased rankA>rankA¯=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 lof the matrix Ais approximated by the weighted normal pseudosolution of system (10) and, if ΔAMNμl<12, the projection error is estimated by formula (75).

  • Remark 5. If the rank of the original matrix is unknown, then the δ-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 ncolumns of the matrix Ucan actually contribute to the product (100). Moreover, if some of the weighted singular numbers are equal to zero, then less than ncolumns of Uare needed. If kpis the number of nonzero weighted singular numbers, then Ucan be reduced to the sizes m×kp, D+—to the sizes kp×kp, VT—up to size kp×n. Formally, such matrices Uand Vare not M-orthogonal and N1-orthogonal, respectively, since they are not square. However, their columns are weighted orthonormal systems of vectors.

    Advertisement

    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 r¯=A¯x¯¯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.,rankA=rankA¯.

    Theorem 17. Assume that ΔAMNAMN+NM<1, rankA¯=rankA=kand let x¯RA¯#. Then

    x¯x¯¯NxNh1hεА2εА+α+hεАβ+γ.E103

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

    To estimate the computational error x¯x¯¯, we use the relation

    A¯x¯x¯¯=r¯=b¯kA¯x¯¯,E104

    whereb¯kis projection of the vector b¯on the main left weighted singular subspace of the matrix A¯, i.е. b¯kRA¯.

    Considering that x¯x¯¯RA¯#and the fact that A¯MN+A¯is a projector in RA¯#, we have

    A¯MN+A¯x¯x¯¯=x¯x¯¯=A¯MN+r¯.E105

    From this, we obtain an estimate of the computational error

    x¯x¯¯Nx¯NA¯MNA¯MN+NMr¯Mb¯kM.E106

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

    xx¯¯NxNxx¯NxN+x¯x¯¯NxN,E107
    x¯x¯¯NxNAMNA¯MN+NMr¯MbkME108

    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.,rankA¯>rankA=k.

    IfΔAMNAMN+NM<1, then from [26], it follows that the rank of the perturbed matrix cannot decrease.

    Theorem 18. Assume that ΔAMNAMN+NM<12, rankA¯>rankA=kand letx¯RA¯k#. Then

    xx¯¯NxNh1hεA2εA+α+hεAβ+γk.E109

    Proof. To estimate the computational error x¯kx¯¯N, we use the fact that A¯kx¯k=b¯k. Then for arbitrary vector x¯RA¯k#

    А¯kx¯kx¯¯=r¯k=b¯kA¯kx¯¯,A¯k+А¯kx¯kx¯¯=A¯k+r¯k.E110

    Considering the fact that xkx¯¯RΑ¯k#, and operator A¯kMN+A¯kis the projection operator in RΑ¯k#, we obtain

    A¯kMN+A¯kx¯kx¯¯=x¯kx¯¯=A¯kMN+r¯k,x¯kx¯¯=A¯kMN+r¯k.E111

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

    x¯kx¯¯Nx¯kNA¯kMNA¯kMN+NMr¯kMb¯kM.E112

    The estimate of the total error follows from the inequalities

    xx¯¯NxNxx¯kNxN+x¯kx¯¯NxN,x¯kx¯¯NxNAMNA¯kMN+NMr¯kMbkME113

    and estimates (25), (64).

    Case 3. The rank of the original matrix is larger than that of the perturbed matrix, i.e.,rankA>rankA¯=l.

    Consider the case when the condition ΔAMNAMN+NM<1not satisfied and the rank of the perturbed matrix can decrease.

    Theorem 19. Assume that rankA>rankA¯=l, ΔAMNμl<12and let x¯ImA¯#. Then

    xlx¯¯NxlNμ1/μl12ΔAMN/μl2εA+αl+μ1μlεAβl+γlE114

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

    minxCxlN,C=xAlxbM=minE115

    With matrix Al=UΣlVTwith rang l.

    The estimate of the computational error in this case will be

    x¯x¯¯Nx¯NA¯MNA¯MN+NMr¯Mb¯lM.E116

    The estimate of the total error follows from the inequalities

    xlx¯¯NxlNxlx¯NxlN+x¯x¯¯NxlN,x¯x¯¯NxNAMNA¯lMN+NMr¯lMblM,E117

    obvious relationships AMN=AlMN,AlMN+NM=1μl,estimates of the hereditary error (75) and the inequality ΔAlMN2ΔAMN.

    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 AMN+is represented in accordance with the properties of the full rank matrix (81).

    Theorem 20.Let ΔAMIAMN+IM<1, m>n=kand x¯¯RA¯#.Then

    xx¯¯xh1hεАεА+ΔbMAMIx+hεАrMAMIx+r¯kMAMIx.E118

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

    x¯x¯¯x¯A¯MNA¯MN+NMr¯Mb¯kM.E119

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

    xx¯¯xxx¯x+x¯x¯¯x,x¯x¯¯xAMIA¯MN+IMr¯MbkME120

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

    Theorem 21. Let ΔAINAMN+NI<1, n>m=kand x¯¯RA¯#. Then

    xx¯¯NxNh1hεА2εA+ΔbAINxN+r¯bk,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 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 ΔAMNA¯MN+NM<1, andx¯¯RA¯#. Then, for the total error of the normal pseudosolution, the following estimate is fulfilled

    xx¯¯Nx¯NhA¯1hA¯εА¯2εА¯+ΔbMA¯MNx¯N+hA¯εA¯r¯MA¯MNx¯N+r¯Mb¯kM.E122

    Estimate (122) can be obtained from the inequality

    xx¯¯Nx¯Nxx¯Nx¯N+x¯x¯¯Nx¯NE123

    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 r¯=b¯kА¯x¯¯, the explicit form of the projection operator onto RА¯#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 μα, 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].

    Advertisement

    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.

    References

    1. 1.Chipman JS. On least squares with insufficient observation. Journal of the American Statistical Association. 1964;59(308):1078-1111
    2. 2.Milne RD. An oblique matrix pseudoinverse. SIAM Journal on Applied Mathematics. 1968;16(5):931-944
    3. 3.Ward JF, Boullion TL, Lewis TO. Weighted pseudoinverses with singular weights. SIAM Journal on Applied Mathematics. 1971;21(3):480-482
    4. 4.Galba EF, Deineka VS, Sergienko IV. Weighted pseudoinverses and weighted normal pseudosolutions with singular weights. Computational Mathematics and Mathematical Physics. 2009;49(8):1281-1297
    5. 5.Sergienko IV, Galba EF, Deineka VS. Existence and uniqueness of weighted pseudoinverse matrices and weighted normal pseudosolutions with singular weights. Ukrainian Mathematical Journal. 2011;63(1):98-124
    6. 6.Varenyuk NA, Galba EF, Sergienko IV, Khimich AN. Weighted pseudoinversion with indefinite weights. Ukrainian Mathematical Journal. 2018;70(6):866-889
    7. 7.Galba EF, Varenyuk NA. Representing weighted pseudoinverse matrices with mixed weights in terms of other pseudoinverses. Cybernetics and System Analysis. 2018;54(2):185-192
    8. 8.Galba EF, Varenyuk NA. Expansions of weighted pseudoinverses with mixed weights into matrix power series and power products. Cybernetics and System Analysis. 2019;55:760-771. DOI: 10.1007/s10559-019-00186-9
    9. 9.Goldman AJ, Zelen M. Weak generalized inverses and minimum variance linear unbiased estimation. Journal of Research of the National Bureau of Standards. 1964;68B(4):151-172
    10. 10.Rao CR, Mitra SK. Generalized Inverse of Matrices and its Applications. New York: Wiley; 1971. p. 240
    11. 11.Nashed MZ. Generalized Inverses and Applications. New York: Academic Press; 1976. p. 1068
    12. 12.Ben-Israel A, TNE G. Generalized Inverses: Theory and Applications. New York: Springer-Verlag; 2003. p. 420. DOI: 10.1007/b97366
    13. 13.Khimich AN. Perturbation bounds for the least squares problem. Cybernetics and System Analysis. 1996;32(3):434-436
    14. 14.Khimich AN, Nikolaevskaya EA. Reliability analysis of computer solutions of systems of linear algebraic equations with approximate initial data. Cybernetics and System Analysis. 2008;44(6):863-874
    15. 15.Nikolaevskaya EA, Khimich AN. Error estimation for a weighted minimum-norm least squares solution with positive definite weights. Computational Mathematics and Mathematical Physics. 2009;49(3):409-417
    16. 16.Wei Y, Wang D. Condition numbers and perturbation of the weighted Moore-Penrose inverse and weighted linear least squares problem. Applied Mathematics and Computation. 2003;145(1):45-58
    17. 17.Wei Y. A note on the sensitivity of the solution of the weighted linear least squares problem. Applied Mathematics and Computation. 2003;145(2–3):481-485
    18. 18.Molchanov IN, Galba EF. A weighed pseudoinverse for complex matrices. Ukrainian Mathematical Journal. 1983;35(1):46-50
    19. 19.Elden L. A weighted pseudoinverse, generalized singular values and constrained least squares problems. BIT. 1982;22:487-502
    20. 20.Wei Y. The weighted Moore–Penrose inverse of modified matrices. Applied Mathematics and Computation. 2001;122(1):1-13. DOI: 10.1016/S0096-3003(00)00007-2
    21. 21.Voevodin VV. On the regularization method. Zhurnal Vychislitel′noy Matematiki i Matematicheskoy Fiziki. 1969;9:673-675
    22. 22.Tikhonov AN. Regularization of ill-posed problems. Doklady Akademii Nauk SSSR. 1963;153:42-52. DOI: 10.1016/S0096-3003(00)00007-2
    23. 23.Ivanov VK, Vasin VV, Tanana VP. Theory of Linear Ill-Posed Problems and Applications. Moscow: Nauka; 1978. 206 p
    24. 24.Morozov VA. Regularization Methods for Unstable Problems. Moscow: Moscow State University; 1987 [in Russian]
    25. 25.Albert AE. Regression and the Moore–Penrose Pseudoinverse. New York: Academic Press; 1972. p. 180
    26. 26.Lawson CL, Hanson RJ. Solving Least Squares Problems. Moscow: Nauka; 1986. p. 232
    27. 27.Kirichenko NF. Analytical representation of perturbations of pseudoinverses. Kibernetika i Sistemnyi Analiz. 1997;2:98-107
    28. 28.Golub GH, Van Loan CF. Matrix Computations. Baltimore: Johns Hopkins University Press; 1996. p. 694
    29. 29.Elden L. Perturbation theory for the least squares problem with linear equality constraints. SIAM Journal on Numerical Analysis. 1980;17:338-350
    30. 30.Bjork A. Numerical Methods for Least Squares Problems. Linköping: Linköping University; 1996. p. 407
    31. 31.Voevodin VV. Computational Foundations of Linear Algebra. Moscow: Nauka; 1977
    32. 32.Khimich AN. Estimates of perturbations for least squares solutions. Kibernetika i Sistemnyi Analiz. 1996;3:142-145
    33. 33.Khimich AN, Voitsekhovskii SA, Brusnikin VN. The reliability of solutions to linear mathematical models with approximately given initial data. Mathematical Machines and Systems. 2004;3:3-17
    34. 34.Wei Y, Wu H. Expression for the perturbation of the weighted Moore–Penrose inverse. Computers & Mathematcs with Applications. 2000;39:13-18
    35. 35.Wei M. Supremum and Stability of Weighted Pseudoinverses and Weighted Least Squares Problems: Analysis and Computations. New York: Huntington; 2001. p. 180
    36. 36.Wang D. Some topics on weighted Moore–Penrose inverse, weighted least squares, and weighted regularized Tikhonov problems. Applied Mathematics and Computation. 2004;157:243-267
    37. 37.Van Loan CF. Generalizing the singular value decomposition. SIAM Journal on Numerical Analysis. 1976;13:76-83
    38. 38.Wang G, Wei Y, Qiao S. Generalized Inverses: Theory and Computations. Beijing: Science Press; 2004. p. 390
    39. 39.Albert A. Regression, Pseudoinversion, and Recurrent Estimation. Moscow: Nauka; 1977. p. 305
    40. 40.Khimich AN, Nikolaevskaya EA. Error estimation for weighted least squares solutions. Computational Mathematics. 2006;3:36-45
    41. 41.Khimich AN, Nikolaevskaya EA. Perturbation analysis of weighted least squares solutions. Theory of Optimal Solutions. 2007;6:12-18
    42. 42.Voevodin VV, Kuznecov YA. Matrices and Calculations. Мoscow: Nauka; 1984. p. 318
    43. 43.Wilkinson JH, Reinsch C. Handbook Algorithms in ALGOL. Linear Algebra. Mechanical: Moscow; 1976. p. 389
    44. 44.GNU Multiple Precision Arithmetic Library. Availabe from:www.gmplib.org
    45. 45.Nikolaevskaya EA, Khimich AN, Chistyakova TV. Programming with Multiple Precision. Berlin/Heidelberg, Germany: Springer; 2012

    Written By

    Aleksandr N. Khimich, Elena A. Nikolaevskaya and Igor A. Baranov

    Submitted: January 14th, 2022Reviewed: January 26th, 2022Published: April 16th, 2022