LINEAR ALGEBRA AND ITS APPLICATIONS ELSEVIER Linear Algebra and its Applications 291(1999)185-199 Inverse updating and downdating for weighted linear least squares using M-invariant reflections Weiguo Wang a, Jinxi Zhao b. * Department of Mathematics, Qufu Normal University, Qufu, Shandong 273165, People's Republic of China b Department of Mathematics, Nanjing University, Nanjing 210008, People's Republic of China Received 23 May 1996; accepted 30 November 1998 Submitted by L. Elsner Abstract A new method for the weighted linear least squares problem min, M-12(b - Ax), is presented by introducing a row M-invariant matrix (i.e., QMQ' = M). Our purpose in this paper is to introduce new row M-invariant and row hyperbolic M-invariant re- flections. We then show how these row M-invariant reflections can be used to design efficient sliding-date-window recursive weighted linear least squares covariance algo- rithms, which are based upon rank-k modifications to the inverse like-Cholesky factor R of the covariance matrix. The algorithms are rich in matrix-matrix BLAS-3 com- putations. We also provide computational experiments indicating the numerical stability of the methods. 1999 Elsevier Science Inc. All rights reserved. AMS classification: 65F25; 65F35; 65G05 Keywords:Linear weighted least squares;Row M-invariant reflection;Updating:Downdating 1.Introduction Consider the weighted linear least squares problem [1,7-9] min(b-Ax)'M1(b-Ax), (1) x∈R Corresponding author. Project supported by the State 863-high science technology plan of China. 0024-3795/99/S- see front matter 1999 Elsevier Science Inc. All rights reserved. PII:S0024-3795(99)00007-5
186 W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 where A E R"x,b Rm and rank (A)=n,and M diag(:...,m)>0. It is assumed that m>n.An equivalent formulation of Eq.(1)is (6)()-(6) (2) where E R,xE R",see Ref.[9]for more details. Since rank(4)=n and M is positive definite (2)has a unique solution. Assume that E Rxm and MEx,then is said to be M-invariant if it is nonsingular and OMOT =M. In a pioneering paper by Gulliksson and Wedin [2],it was shown how Eq.(2)could be solved by using the column M-invariant matrices.In this paper we introduce new row M-invariant and row hyperbolic M-invariant matrices,that play the same role as row householder reflections [3]and M- invariant rotations [2]. We will develop efficient algorithms for the recursive weighted least squares problem of the sliding-window type (see Refs.[5,61). The outline of this paper is as follows.In Section 2 we introduce the new row M-invariant methods.In Sections 3 and 4 we show how these matrices can be used to efficiently modify weighted least-squares solutions when observations are added to and/or deleted from the linear system.In Sec- tion 5 some computational experiments and some concluding remarks are provided. 2.Row M-invariant methods In this section we introduce a row M-invariant method which is a rank-1 modification to the identity matrix eliminating all elements in a row of a ma- trix. 2.1.Row M-invariant reflections Lemma 2.1 [2].Assume that =I-2cd,dMd >0.where Q is M-invariant with M nonsingular.Then O=1-2Mdd/dMd,with =I, i.e.,O is a reflector.The matrix O is called an M-invariant reflection
W.Wang,J.Zhao Linear Algebra and its Applications 291 (1999)185-199 187 Theorem 2.2.Let B be an (n +1)x n matrix B }1 in. (3) where D is nonsingular,and M=diag(,...)with >0.Then there is an M-invariant reflection O such that e=(8) (4) Proof.Let c=(e)a=()h We will construct =I-2cdT such that Eq.(4)is satisfied.Then we obtain the relation () 2 cidib+cid D ed b+ed D =(8) (5) giving DTd ub where =(1-2cid)/2c1. By Lemma 2.1,we have c=Md/d Md and c=ud/d Md that inserted in the expression for u gives 2p1d+21d-dMd=0. (6) By expanding dMd and using d=uD-Tb we have d"Md=diu+u(D-Tb)'M(D-Tb) enabling us to rewrite Eg.(6)as 41d+241d-u(DTb)M(D-Tb)=0. Solving this equation for di we get . By choosing u.we obtain d and =I-2MddT/dMd. In order to avoid rounding errors we choose the negative sign of the square root in di. We may choose u =1/lbll2 and if =0,we set I. We have the follwing algorithm for determining the M-invariant reflection
188 W.Wang.J.Zhao I Linear Algebra and its Applications 291 (1999)185-199 Algorithm 1 (RowMR). INPUT:An n-vector b,a nonsingular n x n-matrix D,and the weight matrix M=diag(41,,4n,na+i】 if=0 then d=0,2=I. else Put u=1/llbll2.and solve DTd ub. Compute d=--+aMd.Let d=(d.) en OUTPUT:=I-2MddT/dTMd having the property that the first row of OB is all zeros. Algorithm RowMR will have good numerical properties as long as d is solved by a numerically stable method. In Section 3 we consider the problem of annihilating r rows.This is easily done by applying a sequence of M-invariant reflections described above as 0=2Q-1…Q1. 2.2.Row hyperbolic M-invariant reflections LetΦ=diag(士l),and assume that2∈Rmxm and M∈Rmxm,then is said to be hyperbolic M-invariant if it is nonsingular and OMor=M(see Ref. [10]) Lemma2.3.Assume that Q=中-2cdr,dTMΦd≠0，whereΦ=diag(士l),and O is hyperbolic M-invariant with M nonsingular.Then 2=Φ-2Mddr/d'Mpd,with O2=Φ， i.e.,O is a hyperbolic reflector.We call o a hyperbolic M-invariant reflection. Theorem 2.4.Let B be an (n+1)x n matrix of the form B D》 where D is nonsingular,and M=diag(,....)with >0.Assume that DTM,D-1/u bb>0 then there is a hyperbolic M-invariant reflection 2=Φ-2cd,whereΦ=diag(-l,l,,I),such that QB= (8) (7) where DE Raxn
W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 189 Proof.Let c=()a-() where c,d4∈R,c,deR". From Eq.(7)we get QB=ΦB-2cdrB=（ -bT-2(cid bT +cidD) D-2(cd bT +cd D) and D'd=ub,where 4=-J-2cd (8) 2c1 By Lemma 2.4,we have c=Md/dM d,and c=ud/dMod which, inserted in Eq.(8),gives 241d+24d+dMΦd=0. (9) By expanding dMod=-du+(D-b)'Mi(D-Tb), where M=diag(2,...,),we can rewrite Eq.(9)as ud2+2uud+(D-Tb)Mi(D-Tb)=0. (10) Solving for di gives d=-2A±V公 241 where we have used the assumption DTMD-(1/u)bbT20 giving bDMD-Tb≥0and△=4412(41-(D-Tb)'M1(D-Tb)≥0. By choosing u we attain di,d and =-2MddT/dMd. We arrive at the following algorithm for determining a hyperbolic M-invariant reflection of the kind described above. Algorithm 2(RowHMR). INPUT:An n-vector b,a nonsingular n x n-matrix D,and M=diag(,...,). if lb2=0,then d=0;Q=Φ. else Let u=1/lbll2 and solve DTd ub
190 W.Wang,J.Zhao Linear Algebra and its Applications 291 (1999)185-199 Compute di =u-udMid. Let d=(di,dT). end OUTPUT:=-2Mdd/dMd that has the property that the first row of OB is all zeros. 3.Inverse updating In this section we will solve the updating problem by M-invariant matrix methods. Consider the weighted linear least squares problem below m9aM2(s-Kwl2, (11) where Mi diag(,...u),u>0 and X is an m x n matrix with rank(X)=n. Let X=OR,where o is an M-invariant matrix with columns, R1=(0 ∈代mxn and R is an n x n upper triangular matrix.Then the solution to (11)is given by w=(R10)2's. Suppose k new observations (YT u),where YT e Rex",and ueR,be added to the dating defining the weighted linear least squares problem(11).We then show how the solution w to (11)can be updated to the solution w to m-a(a)-()儿 (12) where M=diag(M,M2). Lemma 3.1.Let ()-(8) (13) where is an M-invariant matrix,R.R is upper triangular,then ()=() (1 where ET E Rkxn
W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 191 Proof.Let (Q()'=w then Eq.(13)gives 1=w(⑧) UR, andU=R-l.▣ Lemma 3.2.Assume that V=-R-TY,R given in Eq.(13),Q is M-invariant,and O=On..0,where Oi are M-invariant reflection with O2=I,such that e()=(8) (15) where Ik is the k x k identity matrix and D is a k xk matrix.Then ()-(8) If U is upper triangular,then U=R,and ()=() Proof.Since is M-invariant,then is M--invariant.We choose row M--invariant matrix O-T,such that e()(8) where D is nonsingular. For =I and the definition of V,we obtain e()() Let e-(888)adg-(88)
192 W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 We get (-RTY)+2=0,O1-YTR-2,=0, 02191+Q2YTR01=0. Since -T is row M-1-invariant matrix,u is nonsingular.We obtain 2+z2YTR-=0 and y=02R+022YT=0. Hence e(R)-(6) If U is upper triangular for e()-(8) then it is easy from Lemma 3.1 to see e()=() Theorem 3.3.Let O satisfy the same assumptions as in Lemma 3.2,if w is the solution to (11),then the solution to (12)is given by w=w-ED-T(u-YTw): (16) with E and D given in Lemmas 3.1 and 3.2. Proof.Let ex=(6)c=() where s E R",S2E R-",is Mi-invariant matrix. Then Eq.(12)can be rewritten as mmw-r()-()儿 =mw(食()儿 (17)
W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 193 which is equivalent to 〔台任)g (18) Hence we have 创 (19) where e=diag(O,1). Since w=R-s1,moreover let V=-R-TY.Consequently,from Lemma 3.2. there exist M-invariant matrix O,such that ()-(8) ndg'(。)-(T further let (a)-() hence w=R-151. Note that w=R-s and hence from the definition of r,we have the relation ()-(）()） -(）eoa)-()） (20) Then we get
194 W.Wang.J.Zhao Linear Algebra and its Applications 291 (1999)185-199 W=w-E元=w-EDT(u-YTw).▣ (21) This algorithm computes w,where w solves Eq.(12). Algorithm 3. INPUT:An upper triangular matrix R-T and w,where QX=R and w solves (11).A new set of k observations (YT u)and M=diag(,...,um). 1.Compute V=-R-TY. 2.Find -T=OTOr..Of.where 2 is a row M-invariant reflection, such that o()-(8) 3.Update R-T to R-T,i.e, (0)-() 4.Update w to w,i.e, w=w-ED-T(u-YTw). The cost in flops for each step is:1.kn2/2,2.15k2n,3.kn2,4.2+2kn with a total cost for Algorithm 3 as kn2+15k2n+2kn+flops.A straightforward implementation of the rank-1 method of Pan and Plemmons [4]would require 3kn2+O(kn)multiplications.Thus,roughly speaking,Algorithm 3 requires less flops when n≥l5k. 4.Inverse downdating Let the matrix Y and the vector s be given by the partition =(@)x=() where 2T∈Rxn,d∈R,Then the problem min llM 2(s-w)ll2 (22) is our downdating problem.Thus,we assume that we have the solution to(11) where Mi=diag(M3,M4)and want the the solution to (22)by our row hy- perbolic M-invariant method. Assume that =O.,where are hyperbolic M-invariant reflections. Then we define o=1...O to be used in the sequel. We have the following lemma that is used to construct the downdating al- gorithm