In this contribution we extend the result of (Markovsky et. al, SIAM J. of Matrix Anal. and Appl., 2005) to the case of weighted cost function. It is shown that the computational complexity of the proposed algorithm is preserved linear in the sample size when the weight matrix is banded with bandwidth that is independent of the sample size

G. Golub

I. Markovsky

S. Huffel Van

English

Crossref

On Weighted Structured Total Least Squares

Markovsky, Ivan

Van Huffel, Sabine

Southampton (e-Prints Soton)

On Weighted Structured Total Least SquaresIvan Markovsky and Sabine Van HuﬀelK.U. Leuven, ESAT-SCD, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium{ivan.markovsky, sabine.vanhuffel}@esat.kuleuven.ac.bewww.esat.kuleuven.ac.be/~imarkovswww.esat.kuleuven.ac.be/sista/members/vanhuffel.htmlAbstract. In this contribution we extend our previous results on thestructured total least squares problem to the case of weighted cost func-tions. It is shown that the computational complexity of the proposedalgorithm is preserved linear in the sample size when the weight matrixis banded with bandwidth that is independent of the sample size.1 IntroductionThe total least squares (TLS) method (Golub and Van Loan, [1], Van Huﬀel andVandewalle, [2])min∆A,∆B,X   ∆A ∆B   2F subject to (A − ∆A)X = B − ∆B , (1)is a solution technique for an overdetermined system of equations AX ≈ B,A ∈ IRm×n, B ∈ IRm×d. It is a natural generalization of the least squaresapproximation method when the data in both A and B is perturbed. The methodhas been generalized in two directions:– weighted total least squaresmin∆A,∆B,X   ∆A ∆B   2W subject to (A − ∆A)X = B − ∆B , (2)where  ∆C 2W := vec (∆C )Wvec(∆C ), W>0, and– structured total least squares (STLS)min∆A,∆B,X   ∆A ∆B   2F subject to (A − ∆A)X = B − ∆B and ∆A ∆B has the same structure as AB .(3)While the basic TLS problem allows for an analytic solution in terms of thesingular value decomposition of the data matrix C := AB ,t h ew e i g h t e da n dstructured TLS problems are solved numerically via local optimization methods.In [4] we show that under a general assumption (see Assumption 1) aboutthe structure, the cost function and ﬁrst derivative of the STLS problem canbe evaluated in O(m) ﬂoating point operations (ﬂops). This allows for eﬃcientI. Lirkov, S. Margenov, and J. Wa´ sniewski (Eds.): LSSC 2005, LNCS 3743, pp. 695–702, 2006.c   Springer-Verlag Berlin Heidelberg 2006696 I. Markovsky and S. Van Huﬀelcomputational algorithms based on standard methods for local optimization. Viaa similar approach, see (Markovsky et al., [5]), the weighted TLS problem canbe solved eﬃciently when the weight matrix W is block-diagonal with blocks ofsize n + d.In this paper, we extend our earlier results on the STLS problem by account-ing for weighted cost function. Thus the weighted TLS problem becomes a specialcase of the considered weighted STLS problem when the data matrix is unstruc-tured. In Sect. 2 we review the results of Markovsky, Van Huﬀel, Pintelon [4].Section 3 presents the necessary modiﬁcations for the weighted STLS problemand Sect. 4 discusses the implementation of the algorithm.2 Review of Results for the STLS ProblemLet S :I Rnp → IRm×(n+d) be an injective function. A matrix C ∈ IRm×(n+d)is said to be S-structured if C ∈ image(S). The vector p for which C = S(p)is called the parameter vector of the structured matrix C. Respectively, IRnp iscalled the parameter space of the structure S. The aim of the STLS problemis to perturb as little as possible a given parameter vector p by a vector ∆p,so that the perturbed structured matrix S(p+∆p) becomes rank deﬁcient withrank at most n.Problem 1 (STLS). Given a data vector p ∈ IRnp, a structure speciﬁcation S :IRnp → IRm×(n+d), and a rank speciﬁcation n, solve the optimization problemˆ X =a r gm i nX,∆p ∆p 22 subject to S(p − ∆p) X−Id =0 . (4)Let AB := S(p). Problem 1 makes precise the STLS problem formulation (3)from the introduction. In what follows, we often use the notationXext := X−I .The STLS problem is said to be aﬃne structured if thefunction S is aﬃne,i.e.,S(p)=S0 +np  i=1Sipi, for all p ∈ IRnp and for some Si,i =1 ,...,n p . (5)In an aﬃne STLS problem, the constraint S(p − ∆p)Xext = 0 is bilinear in thedecision variables X and ∆p.Lemma 1. Let S :I Rnp → IRm×(n+d) be an aﬃne function. ThenS(p − ∆p)Xext =0 ⇐⇒ G(X)∆p = r(X) ,whereG(X): = vec (S1Xext)  ··· vec (SnpXext)   ∈ IRmd×np , (6)andr(X): =v e c  S(p)Xext   ∈ IRmd .On Weighted Structured Total Least Squares 697Using Lemma 1, we rewrite the aﬃne STLS problem as followsminX min∆p ∆p 22 subject to G(X)∆p = r(X) . (7)The inner minimization problem has an analytic solution, which allows to derivean equivalent optimization problem.Theorem 1 (Equivalent optimization problem for aﬃne STLS). Assum-ing that np ≥ md, the aﬃne STLS problem (7) is equivalent tominXf(X) where f(X): =r (X)Γ †(X)r(X) and Γ(X): =G(X)G (X) .The signiﬁcance of Theorem 1 is that the constraint and the decision variable ∆pin problem (7) are eliminated. Typically the number of elements nd in X is muchsmaller than the number of elements np in the correction ∆p. Thus the reductionin the complexity is signiﬁcant.The equivalent optimization problem (1) is a nonlinear least squares prob-lem, so that classical optimization methods can be used for its solution. Theoptimization methods require a cost function and ﬁrst derivative evaluation.In order to evaluate the cost function f for a given value of the argument X,we need to form the weight matrix Γ(X) and to solve the system of equationsΓ(X)y(X)=r(X). This straightforward implementation requires O(m3)ﬂ o p s .For large m (the applications that we aim at) this computational complexitybecomes prohibitive.It turns out, however, that for a special case of aﬃne structures S,t h ew e i g h tmatrix Γ(X) has a block-Toeplitz and block-banded structure, which can beexploited for eﬃcient cost function and ﬁrst derivative evaluations.Assumption 1 (Flexible structure speciﬁcation). The structure speciﬁ-cation S :I Rnp → IRm×(n+d) is such that for all p ∈ IRnp, the data matrixS(p)= :C =: AB is of the typeS(p)= C1 ··· Cq , where Cl, for l =1 ,...,q, is block-Toeplitz,block-Hankel, unstructured, or exact and all block-Toeplitz/Hankelstructured blocks Cl have equal row dimension K of the blocks.Assumption 1 says that S(p) is composed of blocks, each one of which is block-Toeplitz, block-Hankel, unstructured, or exact. A block Cl that is exact is notmodiﬁed in the solution ˆ C := S(p − ∆p), i.e., ˆ Cl = Cl. Assumption 1 is theessential structural assumption that we impose on the STLS problem. It is fairlygeneral and covers many applications.We use the notation nl for the number of block columns of the block Cl.F o runstructured and exact blocks nl := 1.Theorem 2 (Structure of the weight matrix Γ). Consider the equivalentoptimization problem (1) from Theorem 1. If in addition to the assumptions of698 I. Markovsky and S. Van HuﬀelTheorem 1, the structure S is such that Assumption 1 holds, then the weightmatrix Γ(X) has the block-Toeplitz and block-banded structure,Γ(X)=⎡⎢ ⎢⎢⎢⎢⎢⎢ ⎢⎢⎢⎣Γ0 Γ  1 ··· Γ  s 0Γ1............. . ............. Γ  sΓs............. . ............. Γ  10 Γs ··· Γ1 Γ0⎤⎥ ⎥⎥⎥⎥⎥⎥ ⎥⎥⎥⎦∈ IRmd×md , (8)where Γk ∈ IRdK×dK,f o rk =0 ,1,...,s,a n ds =m a x l=1,...,q(nl − 1).3 Modiﬁcations for the Weighted STLS ProblemNext we consider the generalization of the STLS problem where the cost functionis weighted.Problem 2 (Weighted STLS). Given a data vector p ∈ IRnp, a positive deﬁniteweight matrix W ∈ IRnp×np, a structure speciﬁcation S :I Rnp → IRm×(n+d),and a rank speciﬁcation n, solve the optimization problemˆ Xw =a r gm i nX,∆p∆p W∆p subject to S(p − ∆p) X−Id =0 . (9)The counterpart of Theorem 1 for the case at hand is the following one.Theorem 3 (Equivalent optimization problem for weighted STLS).Assuming that np ≥ md, the aﬃne weighted STLS problem (9) is equivalent tominXfw(X) where fw(X): =r (X)Γ†w(X)r(X)and Γw(X): =G(X)W −1G (X) . (10)Proof. The equivalent optimization problem in Theorem 1 is obtained by solvinga least squares problem. In the weighted STLS case, we solve the weighted leastsquares problemmin∆p∆p W∆p subject to G(X)∆p = r(X) .The optimal parameter correction as a function of X is∆pw(X)=W −1G (X) G(X)W −1G (X) †r(X) ,so thatfw(X)=∆p w(X)W∆pw(X)=r (X) G(X)W −1G (X)      Γw †r(X) .    On Weighted Structured Total Least Squares 699In general neither the block-Toeplitz nor the block-banded properties of Γ =GG  are present in Γw = GW −1G . In the rest of this section, we show that incertain special cases these properties are preserved.Assumption 2 (Block-diagonal weight matrix). Consider the ﬂexiblestructure speciﬁcation of Assumption 1, let the blocks Cl, l =1 ,...,q be param-eterized by parameter vectors pl ∈ IRnp,l, and assume without loss of generalitythat p =c o l ( p1,...,p q). The weight matrix W is assumed to be block-diagonalW = blkdiag(W1,...,Wq), where W l ∈ IRnp,l×np,l .Assumption 2 forbids cross-weighting among the parameters of the blocksC1,...,Cq. Under Assumption 2 the eﬀect of Cl on Γw is independent fromthose of the other blocks. Thus the problem of determining the structure of Γw,resulting from the ﬂexible structure speciﬁcation of C decouples into three in-dependent problems: what is the structure of Γw, resulting from respectively anunstructured matrix C, a block-Hankel matrix C, and a block-Toeplitz matrix C.In what follows “i-block-Toeplitz matrix” stands for block-Toeplitz matrixwith i× i block size and “s-block-banded matrix” stands for a block-symmetricand block-banded matrix with upper/lower block-bandwidth i.L e tV := W−1and Vi := W−1i .Proposition 1. Let G be deﬁned as in (6) and let Assumptions 1 and 2 hold.If all blocks V l corresponding to unstructured blocks Cl are (n + d)K-block-Toeplitz and all blocks V l corresponding to block-Toeplitz/Hankel blocks Cl aredK-block-Toeplitz, Γw = GV G  is dK-block-Toeplitz.Proof. See the Appendix.    For a particular type of weight matrices, the block-Toeplitz structure of Γ ispreserved. More important, however, is the implication of the following propo-sition.Proposition 2. Let G be deﬁned as in (6) and let Assumptions 1 and 2 hold.If W is p-block-banded, then Γw = GV G  is (s + p)-block-banded, where s isgiven in Theorem 2.Proof. See the Appendix.    For block-banded weight matrix W, the block-banded structure of Γ is preserved,however, the block-bandwidth is increased by the block-bandwidth of W.I nt h efollowing section, the block-banded structure of Γ (and Γw) is utilized for O(m)cost function and ﬁrst derivative evaluation.Summary: We have established the following special cases:V block-Toeplitz =⇒Γw block-Toeplitz (generally not block-banded),V l p-block-banded=⇒Γw (s + p)-block-banded (generally not block-Toeplitz),W block-diagonal =⇒Γw s-block-banded (generally not block-Toeplitz).The case W block-diagonal, i.e., Wl = blkdiag(Wl1,...,Wlm), for l =1 ,...,q,covers most applications of interest and will be considered in the next section.700 I. Markovsky and S. Van Huﬀel4 Algorithm for Solving Weighted STLS ProblemIn [3] we have proposed an algorithm for solving the STLS problem (4) with theﬂexible structure speciﬁcation of Assumption 1. The structure of S(·) is speciﬁedby the integer K, the number of rows in a block of a block-Toeplitz/Hankelstructured block Cl, and the array S ∈ ({T,H,U,E}×IN × IN)q that describesthe structure of the blocks {Cl}ql=1.T h elth element Sl of the array S speciﬁesthe block Cl by giving its type Sl(1), the number of columns nl = Sl(2), and (ifCl is block-Hankel or block-Toeplitz) the column dimension tl = Sl(3) of a blockin Cl. Therefore, the input data for the STLS problem is the data matrix S(p)(alternatively the parameter vector p) and the structure speciﬁcation K and S.It is shown that the blocks Γk of Γ are quadratic functions of XΓk(X)=( IK ⊗ X ext)Sk(IK ⊗ X ext)  ,k =0 ,1,...,s , (11)where the matrices Sk ∈ IRK(n+d)×K(n+d) depend on the structure S. The ﬁrststep of the algorithm is to translate the structure speciﬁcation S to the set ofmatrices Sk, k =0 ,1,...,s.T h e nf o rag i v e nX,t h eΓk matrices can be formed,which speciﬁes the Γ matrix.For cost function evaluation, the structured system of equations Γ(X)y(X)=r(X) is solved and the product f(X)=r (X)y(X) is computed. Eﬃciency isachieved by exploiting the structure of Γ in solving the system of equations.Moreover, as shown in [3], the ﬁrst derivative f (X) can also be evaluatedfrom y(X)w i t hO(m) extra computations. The resulting solution method isoutlined in Algorithm 1.Algorithm 1. Algorithm for solving the STLS problemInput: structure speciﬁcation K,S and matrices A and B, such that￿AB￿= S(p).1: Form the matrices {Sk}.2: Compute the TLS solution Xini of AX ≈ B.3: Execute a standard optimization algorithm, e.g., the BFGS quasi-Newtonmethod, for the minimization of f0 over X with initial approximation Xini andwith eﬃcient cost function and ﬁrst derivative evaluation.Output: ˆ X the approximation found by the optimization algorithm uponconvergence.The changes for the case of weighted STLS problem are only in (11). Nowthe matrix Γ is replaced by Γw, which is no longer block-Toeplitz but is s-block-banded with block-elementsΓij(X)=( IK ⊗ X ext)Sij(IK ⊗ X ext)  , (12)whereSij :=⎧⎪ ⎨⎪ ⎩blkdiag(V 1i ,...,Vqi )Si−j if 0 ≤ i − j ≤ sS ji if − s ≤ i − j<00o t h e r w i s e.On Weighted Structured Total Least Squares 701For an exact block Cl, with some abuse of notation, we deﬁne Wli =0f o ra l li.(Our previous deﬁnition of Wl is an empty matrix since pl = 0 in this case.)5 ConclusionsWe have extended the theory of Markovsky, Van Huﬀel, Pintelon ([4]) for thecase of weighted STLS problems. The main question of interest is what proper-ties of the Γ matrix in the equivalent optimization problem are preserved whenthe cost function is weighted. Block-Toeplitz inverse weight matrix V ,r e s u l t sin corresponding Γw matrix that is also block-Toeplitz. More important for fastcomputational methods, however, is the fact that block-banded weight matrix Wwith block-bandwidth p leads to increase of the block-bandwidth of Γ with p.In particular W block-diagonal, results in Γw block-banded with the same band-width as Γ. This observation was used for eﬃcient solution of weighted STLSproblems.AcknowledgmentsDr. Ivan Markovsky is a postdoctoral researcher and Dr. Sabine Van Huﬀel is afull professor at the Katholieke Universiteit Leuven, Belgium.Our research is supported by Research Council KUL: GOA-AMBioRICS,GOA-Meﬁsto 666, several PhD/postdoc & fellow grants; Flemish Govern-ment: FWO: PhD/postdoc grants, projects, G.0078.01 (structured matri-ces), G.0407.02 (support vector machines), G.0269.02 (magnetic resonancespectroscopic imaging), G.0270.02 (nonlinear Lp approximation), G.0360.05(EEG signal processing), research communities (ICCoS, ANMMM); IWT: PhDGrants; Belgian Federal Science Policy Oﬃce IUAP P5/22 (‘Dynamical Systemsand Control: Computation, Identiﬁcation and Modelling’); EU: PDT-COIL,BIOPATTERN, ETUMOUR; HPC-EUROPA (RII3-CT-2003-506079), with thesupport of the European Community — Research Infrastructure Action underthe FP6 “Structuring the European Research Area” Program.References1. Golub, G., Van Loan, C.: An analysis of the total least squares problem. SIAM J.Numer. Anal. 17 (1980) 883–8932. Van Huﬀel, S., Vandewalle, J.: The total least squares problem: Computationalaspects and analysis. SIAM, Philadelphia (1991)3. Markovsky, I., Van Huﬀel, S.: High-performance numerical algorithms and softwarefor structured total least squares. J. of Comput. and Appl. Math. 180 2 (2005)311–3314. Markovsky, I., Van Huﬀel, S., Pintelon, R.: Block-Toeplitz/Hankel structured totalleast squares. SIAM J. Matrix Anal. Appl. 26 4 (2005) 1083–10995. Markovsky, I., Rastello, M.L., Premoli, A., Kukush, A., Van Huﬀel, S.: The element-wise weighted total least squares problem. Comp. Stat. & Data Anal. 51 1 (2005)181–209702 I. Markovsky and S. Van HuﬀelA Proof of Propositions 1 and 2The G matrix, see (6), has the following structure G = G1 ··· Gq ,w h e r eGl ∈IRmd×pl depends only on the structure of Cl (see Lemma 3.2 of [4]). For anunstructured block Cl,Gl = Im ⊗ X ext,l , (13)where Xext =: col(Xext,1,...,X ext,q), Xext,l ∈ IRnl×d and ⊗ is the Kroneckerproduct. For a block-Toeplitz block Cl,Gl =⎡⎢⎢⎢⎢⎣X1 X2 ··· Xnl 0 ··· 00 X1 X2 ··· Xnl.... . .. . ............. 00 ··· 0 X1 X2 ··· Xnl⎤⎥⎥⎥⎥⎦, (14)where Xk := IK ⊗ Xext,k, and for a block-Hankel block Cl,Gl =⎡⎢⎢⎢⎢⎣Xnl Xnl−1 ··· X1 0 ··· 00 Xnl Xnl−1 ··· X1.... . .. . ............. 00 ··· 0 Xnl Xnl−1 ··· X1⎤⎥⎥⎥⎥⎦. (15)Due to Assumption 2, we haveΓw = GV G  =q  l=1GlV l(Gl)       Γ lw, (16)so that we need to consider the three independent problems: structure of Γ lw forrespectively unstructured, block-Toeplitz, and block-Hankel block Cl.T h es t a t e -ments of Propositions 1 and 2 are now easy to see by substituting respectively(13), (14), and (15) in (16) and doing the matrix products.

An analysis of the total least squares problem.

On weighted structured total least squares

On weighted structured total least squares

Abstract

Similar works

Full text

Available Versions

Crossref

Southampton (e-Prints Soton)