Ivan Markovsky Low-Rank Approximation Algorithms, Implementation, Applications September2,2014 Springer vi Preface Preface mationproblems.Thetheoryandalgorithmsofthesenewclassesofproblemsare interestingintheirownrightandbeingapplicationdrivenarepracticallyrelevant. Stochastic estimation and deterministic approximation are two complementary aspectsofdatamodeling.Theformeraimstofindfromnoisydata,generatedbya low-complexitysystem,anestimateofthatdatageneratingsystem.Thelatteraims tofindfromexactdata,generatedbyahighcomplexitysystem,alow-complexity approximation of the data generating system. In applications both the stochastic estimation and deterministic approximation aspects are likely to be present. The dataislikelytobeimpreciseduetomeasurementerrorsandislikelytobegener- atedbyacomplicatedphenomenonthatisnotexactlyrepresentablebyamodelin theconsideredmodelclass.Thedevelopmentofdatamodelingmethodsinsystem identification andsignalprocessing,however,hasbeendominatedbythestochas- tic estimation point of view. If considered, the approximation error is represented inthemainstreamdatamodelingliteratureasarandomprocess.Thisisnotnatural Mathematicalmodelsareobtainedfromfirstprinciples(naturallaws,interconnec- becausethe approximation error is by definition deterministic and even if consid- tion,etc.)andexperimentaldata.Modelingfromfirstprinciplesiscommoninnat- ered as a random process, it is not likely to satisfy standard stochastic regularity ural sciences, while modeling from data is common in engineering. In engineer- conditionssuchaszeromean,stationarity,ergodicity,andGaussianity. ing, often experimental data is available and a simple approximate model is pre- Anexceptiontothestochasticparadigmindatamodelingisthebehavioralap- ferredtoacomplicateddetailedone.Indeed,althoughoptimalpredictionandcon- proach,initiatedbyJ.C.Willemsinthemid80’s.Althoughthebehavioralapproach trolofacomplex(high-order,nonlinear,time-varying)systemiscurrentlydifficult is motivated by the deterministic approximation aspect of data modeling, it does toachieve,robustanalysisanddesignmethods,basedonasimple(low-order,lin- not exclude the stochastic estimation approach. In this book, we use the behav- ear,time-invariant)approximatemodel,mayachievesufficientlyhighperformance. ioralapproachasalanguagefordefiningdifferentmodelingproblemsandpresent- Thisbookaddressestheproblemofdataapproximationbylow-complexitymodels. ingtheirsolutions.Weemphasizetheimportanceofdeterministicapproximationin A unifying theme of the book is low-rank approximation: a prototypical data datamodeling,however,weformulateandsolvestochasticestimationproblemsas modelingproblem.Therankofamatrixconstructedfromthedatacorrespondsto low-rankapproximationproblems. thecomplexityofalinearmodelthatfitsthedataexactly.Thedatamatrixbeingfull Manywellknownconceptsandproblemsfromsystemsandcontrol,signalpro- rankimpliesthatthereisnoexactlowcomplexitylinearmodelforthatdata.Inthis cessing, and machine learning reduce to low-rank approximation. Generic exam- case,theaimistofindanapproximatemodel.Oneapproachforapproximatemod- plesinsystemtheoryaremodelreductionandsystemidentification.Theprincipal eling,consideredinthebook,istofindsmall(insomespecifiedsense)modification componentanalysismethodinmachinelearningisequivalenttolow-rankapproxi- ofthedatathatrendersthemodifieddataexact.Theexactmodelforthemodified mation,whichsuggeststhatrelateddimensionalityreduction,classification,andin- dataisanoptimal(inthespecifiedsense)approximatemodelfortheoriginaldata. formationretrievalproblemscanbephrasedaslow-rankapproximationproblems. Thecorrespondingcomputationalproblemislow-rankapproximation.Itallowsthe Sylvesterstructuredlow-rankapproximationhasapplicationsincomputationswith usertotradeoffaccuracyvscomplexitybyvaryingtherankoftheapproximation. polynomialsandisrelatedtomethodsfromcomputeralgebra. Thedistancemeasureforthedatamodificationisauserchoicethatspecifiesthe The developed ideas lead to algorithms, which are implemented in software. desiredapproximationcriterionorreflectspriorknowledgeabouttheaccuracyofthe The algorithms clarify the ideas and the software implementation clarifies the al- data.Inaddition,theusermayhavepriorknowledgeaboutthesystemthatgenerates gorithms.Indeed,thesoftwareistheultimateunambiguousdescriptionofhowthe thedata.Suchknowledgecanbeincorporatedinthemodelingproblembyimposing ideas are put to work. In addition, the provided software allows the reader to re- constraintsonthemodel.Forexample,ifthemodelisknown(orpostulated)tobe producetheexamplesinthebookandtomodifythem.Theexpositionreflectsthe alineartime-invariantdynamicalsystem,thedatamatrixhasHankelstructureand sequence the approximating matrix should have the same structure. This leads to a Hankel theory algorithms implementation. structuredlow-rankapproximationproblem. 7→ 7→ A tenet of the book is: the estimation accuracy of the basic low-rank approx- Correspondingly,thetextisinterwovenwithcodethatgeneratesthenumericalex- imation method can be improved by exploiting prior knowledge, i.e., by adding amplesbeingdiscussed. constraintsthatareknowntoholdforthedatageneratingsystem.Thispathofde- velopment leads to weighted, structured, and other constrained low-rank approxi- v Preface vii viii Preface Prerequisites andpracticeproblems in MATLABin order to generate the figures and numerical results. Finally, LATEX, bibTEX,dvips, andps2pdf arerunonthefiles createdintheprevious stepsto generate the final pdf file of the book. The described process guarantees that the Acommonfeatureofthecurrentresearchactivityinallareasofscienceandengi- codeinthetextistheactualcodethathasgeneratedtheresultsshowninthebook. neeringisthenarrowspecialization.Inthisbook,wepickapplicationsinthebroad Alaptopandadesktopcomputerswereusedforworkingonthebook.Keeping area of data modeling, posing and solving them as low-rank approximation prob- the files synchronized betweenthetwo computers was done usingthe Unisonfile lems.Thisunifiesseeminglyunrelatedapplicationsandsolutiontechniquesbyem- synchronizer. phasisingtheircommonaspects(e.g.,complexity–accuracytrade-off)andabstract- ing from theapplicationspecificdetails,terminology, andimplementation details. Despiteofthefactthatapplicationsinsystemsandcontrol,signalprocessing,ma- chinelearning,andcomputervisionareusedasexamples,theonlyrealprerequisites Acknowledgements forfollowingthepresentationisknowledgeoflinearalgebra. Thebookisintendedtobeusedforselfstudybyresearchersintheareaofdata AnumberofindividualsandtheEuropeanResearchCouncilcontributedandsup- modelingandbyadvancedundergraduate/graduatelevelstudentsasacomplemen- ported me during the preparation of the book. Oliver Jackson—Springer’s edi- tary text for acourse on systemidentification or machine learning. In either case, tor (engineering)—encouraged me to embark on the project. My colleagues in theexpectedknowledgeisundergraduatelevellinearalgebra.Inaddition,MATLAB ESAT/SISTA,K.U.LeuvenandECS/ISIS,Southampton,U.K.createdtherighten- codeisused,sothatfamiliaritywithMATLABprogramminglanguageisrequired. vironmentfordevelopingtheideasinthebook.Inparticular,IamindebttoJanC. Passivereadingofthebookgivesabroadperspectiveonthesubject.Deeperun- Willems (SISTA) for his personal guidance and example of critical thinking. The derstanding,however,requiresactiveinvolvement,suchassupplyingmissingjusti- behavioralapproachthatJaninitiatedintheearly1980’sispresentinthisbook. ficationofstatementsandspecificexamplesofthegeneralconcepts,applicationand Maarten De Vos, Diana Sima, Konstantin Usevich, and Jan Willems proofread modificationofpresentedideas,andsolutionoftheprovidedexercisesandpractice chaptersofthebookandsuggestedimprovements.Igratefullyacknowledgefunding problems.Therearetwotypesofpracticeproblems:analytical,askingforaproof fromtheEuropeanResearchCouncilundertheEuropeanUnion’sSeventhFrame- ofastatementclarifyingorexpandingthematerialinthebook,andcomputational, workProgramme(FP7/2007-2013)/ERCGrantagreementnumber258581“Struc- askingforexperimentswithrealorsimulateddataofspecificapplications.Mostof turedlow-rankapproximation:Theory,algorithms,andapplications”. theproblemsareeasytomediumdifficulty.Afewproblems(markedwithstars)can beusedassmallresearchprojects. Thecodeinthebook,availablefrom Southampton, IvanMarkovsky http://extra.springer.com/ September2,2014 hasbeentestedwithMATLAB7.9,runningunderLinux,andusestheOptimization Toolbox4.3,ControlSystemToolbox8.4,andSymbolicMathToolbox5.3.Aver- sion of thecodethat is compatible with Octave(afree alternative to MATLAB)is alsoavailablefromthebook’swebpage. Softwaretoolsusedfortypesettingthebook ThebookistypesetusingLATEXandanumberofextensionpackages.Diagramsare createdusingthePSTricksandxypackagesandthestandalonexfigprogram.All editingisdoneintheGNUEmacseditor,usingorg-modeandlatex-mode.MATLAB code,presentedinaliterateprogrammingstyleisanintegralpartofthetext.Thisis achievedbytheNormanRamsey’snowebsystem. The process of compiling the book involves a number of steps automated by a GNUMakefile.FirsttheLATEXsourcecodeandtheMATLABfunctionsandscripts are extracted from noweb’s source files. The obtained scripts are then executed x Contents Contents PartII Miscellaneousgeneralizations 5 Missingdata,centering,andconstraints ..........................137 5.1 Weightedlow-rankapproximationwithmissingdata .............138 5.2 Affinedatamodeling........................................150 5.3 Complexleastsquaresproblemwithconstrainedphase...........159 5.4 Approximatelowrankfactorizationwithstructuredfactors........165 5.5 Notesandreferences........................................177 6 Nonlinearstaticdatamodeling ..................................181 6.1 Aframeworkfornonlinearstaticdatamodeling .................181 6.2 Nonlinearlow-rankapproximation ............................186 6.3 Algorithms................................................189 6.4 Examples .................................................193 6.5 Notesandreferences........................................198 1 Introduction................................................... 1 1.1 Classicalandbehavioralparadigmsfordatamodeling............ 1 7 Fastmeasurementsofslowprocesses .............................201 1.2 Motivatingexampleforlow-rankapproximation................. 3 7.1 Introduction ...............................................201 1.3 Overviewofapplications .................................... 7 7.2 Estimationwithknownmeasurementprocessdynamics ..........204 1.4 Overviewofalgorithms ..................................... 21 7.3 Estimationwithunknownmeasurementprocessdynamics ........206 1.5 Literateprogramming....................................... 24 7.4 Examplesandreal-lifetesting ................................214 1.6 Notesandreferences........................................ 27 7.5 Notesandreferences........................................225 References.....................................................227 PartI Linearmodelingproblems A Approximatesolutionofanoverdeterminedsystemofequations.....229 2 Fromdatatomodels............................................ 37 2.1 Linearstaticmodelrepresentations............................ 37 B Proofs ........................................................237 2.2 Lineartime-invariantmodelrepresentations .................... 47 2.3 Exactandapproximatedatamodeling ......................... 54 P Problems......................................................241 2.4 Unstructuredlow-rankapproximation ......................... 62 2.5 Structuredlow-rankapproximation............................ 68 S Solutions......................................................247 2.6 Notesandreferences........................................ 72 Notation...........................................................263 3 Algorithms .................................................... 75 3.1 Subspacemethods.......................................... 75 Listofcodechunks .................................................265 3.2 Algorithmsbasedonlocaloptimization ........................ 83 Functionsandscriptsindex..........................................269 3.3 Datamodelingusingthenuclearnormheuristic ................. 98 3.4 Notesandreferences........................................106 Index .............................................................271 4 Applicationsinsystem,control,andsignalprocessing ..............109 4.1 Introduction ...............................................109 4.2 Modelreduction ...........................................110 4.3 Systemidentification........................................117 4.4 Analysisandsynthesis ......................................124 4.5 Simulationexamples........................................127 4.6 Notesandreferences........................................130 ix 2 1 Introduction Chapter 1 Toseethatexistenceofalowcomplexityexactlinearmodelisequivalenttorank Introduction deficiencyofthedatamatrix,letthecolumnsd1,...,dNofDbetheobservationsand theelementsd ,...,d ofd betheobservedvariables.Weassumethatthereareat 1j qj j leastasmanyobservationsasobservedvariables,i.e.,q N.AlinearmodelforD ≤ declaresthattherearelinearrelationsamongthevariables,i.e.,therearevectorsr , k suchthat rk⊤dj=0, for j=1,...,N. If there are p independent linear relations, then D has rank less than or equal to m:=q pandtheobservationsbelongtoatmostm-dimensionalsubspaceBofRq. Theveryartofmathematicsistosaythesamething Weide−ntifythemodelforD,definedbythelinearrelationsr1,...,rp Rq,withthe anotherway. set B Rq. Oncea model B is obtained from the data,all possibl∈einput/output ⊂ partitionscanbeenumerated,whichisananalysisproblemfortheidentifiedmodel. Unknown Therefore,thechoiceofaninput/outputpartitioninthebehavioralparadigmtodata modeling can be incorporated, if desired, in the modeling problem and thus need 1.1 Classicalandbehavioralparadigmsfordatamodeling notbehypothesizedasnecessarilydoneintheclassicalparadigm. The classical and behavioral paradigms for data modeling are related but not Fittinglinearmodelstodatacanbeachieved,bothconceptuallyandalgorithmically, equivalent.AlthoughexistenceofsolutionofthesystemAX=Bimpliesthatthema- bysolvingapproximatelyasystemoflinearequations trix AB islowrank,itisnottruethat AB havingasufficientlylowrankimplies that the system AX =B is solvable. This lack of equivalence causes ill-posed (or (cid:2) (cid:3) (cid:2) (cid:3) AX B, (LSE) numericallyill-conditioned)datafittingproblemsintheclassicalparadigm,which ≈ have no solution (or are numerically difficult to solve). In terms of the data fit- where thematrices A andB are constructedfrom the given dataandthematrix X ting problem, ill-conditioning of the problem (LSE) means that the a priori fixed parametrizesthemodel.Inthisclassicalparadigm,themaintoolsaretheordinary input/output partition of the variables is not corroborated by the data. In the be- linearleastsquaresmethodanditsvariations—regularizedleastsquares,totalleast havioral setting without the a priori fixed input/output partition of the variables, squares, robust least squares, etc. The least squares method and its variations are ill-conditioning of the data matrix D implies that the data approximately satisfies mainlymotivatedbytheirapplicationsfordatafitting,buttheyinvariablyconsider linearrelations,sothatnearlyrankdeficiencyisagoodfeatureofthedata. solvingapproximatelyanoverdeterminedsystemofequations. Theclassicalparadigmisincludedinthebehavioralparadigmasaspecialcase Theunderlyingpremiseintheclassicalparadigmisthatexistenceofanexactlin- because approximate solution of an overdetermined system of equations (LSE) is earmodelforthedataisequivalenttoexistenceofsolutionX toasystemAX =B. a possible approach to achieve low-rank approximation. Alternatively, low-rank Suchamodelisalinearmap:thevariablescorrespondingtotheAmatrixareinputs approximation can be achieved by approximating the data matrix with a matrix (or causes)andthevariables correspondingtotheBmatrix areoutputs (or conse- thathasatleastp-dimensionalnullspace,oratmostm-dimensionalcolumnspace. quences) in the sensethat they are determined by the inputs and the model. Note Parametrizingthenullspaceandthecolumnspacebysetsofbasisvectors,theal- thatintheclassicalparadigmtheinput/outputpartitionofthevariablesispostulated ternativeapproachesare: a priori. Unless the model is required to have the a priori specified input/output partition, imposing such structure in advance is ad hoc and leads to undesirable 1. kernelrepresentation thereisafullrowrankmatrixR Rp×q,suchthat ∈ theoreticalandnumericalfeaturesofthemodelingmethodsderived. RD=0, An alternative to the classical paradigm that does not impose an a priori fixed input/outputpartitionisthebehavioralparadigm.Inthebehavioralparadigm,fitting 2. imagerepresentation therearematricesP Rq mandL Rm N,suchthat linear models to data is equivalent to the problem of approximating a matrix D, ∈ × ∈ × constructed from the data, by a matrix D of lower rank. Indeed, existence of an D=PL. exactlinearmodelforDisequivalenttoDbeingrankdeficient.Moreover,therank ofDisrelatedtothecomplexityofthemobdel.Thisfactisthetenetofthebookand Theapproachesusingkernelandimagerepresentationsareequivalenttotheorig- is revisited in the following chapters in the context of applications from systems inallow-rankapproximationproblem.Next,theuseofAX =B,kernel,andimage and control, signal processing, computer algebra, and machine learning. Also its representationsisillustratedonthemostsimpledatafittingproblem—linefitting. implicationtothedevelopmentofnumericalalgorithmsfordatafittingisexplored. 1 1.2 Motivatingexampleforlow-rankapproximation 3 4 1 Introduction 1.2 Motivatingexampleforlow-rankapproximation Thisfeatureoftheclassicalapproachisundesirable:itismorenaturaltospecifya desiredfittingcriterionindependentlyofhowthemodelhappenstobeparametrized. Given a (multi)set of points d ,...,d R2 in the plane, the aim of the line Inmanydatamodelingmethods,however,amodelrepresentationisapriorifixed 1 N { }⊂ andimplicitlycorrespondstoaparticularfittingcriterion. fitting problem is to findaline passingthrough theorigin that “best”matches the givenpoints.Theclassicalapproachforlinefittingistodefine axls=bfit a=bx′lsfit a col(aj,bj):= bj :=dj 6 6 j (cid:20) (cid:21) 4 4 (“:=”standsfor“bydefinition”,seepage263foralistofnotation) andsolveap- proximatelytheoverdeterminedsystem 2 2 col(a ,...,a )x=col(b ,...,b ) (lse) b 0 b 0 1 N 1 N −2 −2 bytheleastsquaresmethod.Letx betheleastsquaressolutionto(lse).Thenthe ls leastsquaresfittinglineis −4 −4 Bls:= d=col(a,b) R2 axls=b . −6 −6 { ∈ | } −2 0 2 −2 0 2 a a Geometrically, Bls minimizes the sum of the squared vertical distances from the datapointstothefittingline. Fig.1.1 Leastsquaresfits(solidlines)minimizingvertical(leftplot)andhorizontal(rightplot) distances. TheleftplotinFigure1.1showsaparticularexamplewithN=10datapoints. Thedatapointsd1,...,d10 arethecirclesinthefigure,thefitBls isthesolidline, andthefittingerrorse:=ax barethedashedlines.Visuallyoneexpectsthebest Thetotalleastsquaresmethodisanalternativetoleastsquaresmethodforsolv- ls − fittobetheverticalaxis,sominimizingverticaldistancesdoesnotseemappropriate ing approximately anoverdetermined system of linear equations.In terms of data inthisexample. fitting,thetotalleastsquaresmethodminimizesthesumofthesquaredorthogonal Notethatbysolving(lse),a(thefirstcomponentsofthed)istreateddifferently distancesfromthedatapointstothefittingline.Usingthesystemofequations(lse), from b (the second components): b is assumed to be a function of a. This is an linefittingbythetotalleastsquaresmethodleadstotheproblem arbitrarychoice;thedatacanbefittedalsobysolvingapproximatelythesystem a b col(a1,...,aN)=col(b1,...,bN)x, (lse′) minimize overx R, ..1 RN,and ..1 RN ∑N d aj 2 inwhichcaseaisassumedtobeafunctionofb.Letx′lsbetheleastsquaressolution ∈ ab.N∈ bb.N∈ j=1(cid:13)(cid:13)(cid:13) j−(cid:20)bbj(cid:21)(cid:13)(cid:13)(cid:13)2 (tls) to(lse′).Itgivesthefittingline subjectto a x=b , for j=1,...,N.   (cid:13) b (cid:13) j j b b Bl′s:={d=col(a,b)∈R2|a=bx′ls}, However,forthbedatabinFigure1.1thetotalleastsquaresproblemhasnosolution. Informally, the approximate solution is x =∞, which corresponds to a fit by a whichminimizesthesumofthesquaredhorizontaldistances(seetherightplotin tls verticalline.Formally, Figure1.1).ThelineBl′shappenstoachievethedesiredfitintheexample. thetotalleastsquaresproblem(tls)mayhavenosolutionandthereforefailto Intheclassicalapproachfordatafitting,i.e.,solvingapproximatelyasystem giveamodel. oflinearequationsintheleastsquaressense,thechoiceofthemodelrepre- sentationaffectsthefittingcriterion. The use of (lse) in the definition of the total least squares line fitting problem restricts the fitting line to be a graph of a function ax=b for some x R. Thus, ∈ theverticallineisaprioriexcludedasapossiblesolution.Intheexample,theline 1.2 Motivatingexampleforlow-rankapproximation 5 6 1 Introduction minimizingthesumofthesquaredorthogonaldistanceshappenstobethevertical rank [d d ] 1. 1 N ··· ≤ line.Forthisreason,x doesnotexist. tls Any line B passing through the origin can be represented as an image and a Thus,(lraP)and(lraR)areinstances(cid:0)obfoneabnd(cid:1)thesame kernel,i.e.,thereexistmatricesP R2 1andR R1 2,suchthat × × ∈ ∈ B=image(P):= d=Pℓ R2 ℓ R abstractproblem:approximatethedatamatrixDbyalowrankmatrixD. { ∈ | ∈ } and b B=ker(R):= d R2 Rd=0 . { ∈ | } InChapter2,theobservationsmadeinthelinefittingexamplearegeneralizedto Usingtheimagerepresentationofthemodel,thelinefittingproblemofminimizing modelingofq-dimensionaldata.Theunderlyinggoalis: thesumofthesquaredorthogonaldistancesis N minimize overP R2 1and ℓ ℓ R1 N ∑ d d 2 GivenasetofpointsinRq (thedata),findasubspaceofRq ofboundeddi- ∈ × 1 ··· N ∈ × j=1k j− jk2 (lra′P) mension(amodel)thathastheleast(2-norm)distancetothedatapoints. subjectto d =Pℓ , for j=(cid:2)1,...,N. (cid:3) j j b With b Such asubspaceis a(2-norm) optimal fitting model. Generalrepresentationsof a D:= d1 ···dN , D:= d1 ···dN , subspaceinRq arethekernelortheimageofamatrix.Theclassicalleastsquares and total least squares formulations of the data modeling problem exclude some and theFrobeniusno(cid:2)rm, (cid:3) h i k·kF b b b subspaces.TheequationsAX =BandA=BX′,usedintheleastsquaresandtotal leastsquaresproblemformulationstorepresentthesubspace,mightfailtorepresent kEkF:= vec(E) 2= e11 ···eq1 ···e1N ···eqN ⊤ 2, forallE∈Rq×N theoptimalsolution,whilethekernelandimagerepresentationsdonothavesuch (lra′P)ismo(cid:13)(cid:13)recompa(cid:13)(cid:13)ctlyw(cid:13)(cid:13)(cid:13)r(cid:2)ittenas (cid:3) (cid:13)(cid:13)(cid:13) dfoerfidcaietancmy.oTdehliisnsgu.ggeststhatthekernelandimagerepresentationsarebettersuited TheequationsAX=BandA=BX wereintroducedfromanalgorithmicpointof ′ minimize overP∈R2×1andL∈R1×N kD−Dk2F (lra ) view—byusingthem,thedatafittingproblemisturnedintothestandardproblemof subjectto D=PL. P solvingapproximatelyanoverdeterminedsystemoflinearequations.Aninterpreta- b tionoftheseequationsinthedatamodelingcontextisthatinthemodelrepresented Similarly,usingakernelrepresentation,thelinefittingproblem,minimizingthesum bytheequationAX =B,thevariableAisaninputandthevariableBisanoutput. b ofsquaresoftheorthogonaldistancesis Similarly,inthemodelrepresentedbytheequationA=BX ,AisanoutputandBis ′ aninput.Theinput/outputinterpretationhasanintuitiveappealbecauseitimpliesa minimize overR∈R1×2,R6=0,andD∈R2×N kD−Dk2F (lra ) causaldependenceofthevariables:theinputiscausingtheoutput. R Representing the model by an equation AX =B and A=BX , as done in the subjectto RD=0. ′ b b classicalapproach,oneaprioriassumesthattheoptimalfittingmodelhasacertain Contrarytothetotalleastsquaresproblem(tls),problems(lra )and(lra )always input/outputstructure.Theconsequencesare: b P R have(nonunique) solutions.In theexample,solutions are,e.g.,P =col(0,1) and ∗ existenceofexceptional(nongeneric)cases,whichcomplicatethetheory, R = 10 ,whichdescribetheverticalline • ∗ ill-conditioningcausedby“nearly”exceptionalcases,whichleadstolackofnu- • (cid:2) (cid:3) B∗:=image(P∗)=ker(R∗). mericalrobustnessofthealgorithms,and needofregularization,whichleadstoachangeofthespecifiedfittingcriterion. • Theconstraints Theseaspectsoftheclassicalapproacharegenerallyconsideredasinherenttothe data modeling problem. By choosing the alternative image and kernel model rep- D=PL, withP∈R2×1, L∈R1×N and RD=0, withR∈R1×2, R6=0 resentations, the problem of solving approximately an overdetermined system of equationsbecomesalow-rankapproximationproblem,wherethenongenericcases areebquivalenttotheconstraintrank(D)≤1,whichbshowsthatthepoints{d1,...,dN} (andtherelatedissuesofill-conditioningandneedofregularization)areavoided. beingfittedexactlybyalinepassingthroughtheoriginisequivalentto b b b 1.3 Overviewofapplications 7 8 1 Introduction 1.3 Overviewofapplications exactdeterministic approximatedeterministic → ↓ ↓ exactstochastic approximatestochastic → Inthissection,examplesoflow-rankapproximationdrawnfromdifferentapplica- tionareasarelisted.Thefactthatamatrixconstructedfromexactdataislowrank Fig.1.2 Transitionsamongexactdeterministic,approximatedeterministic,exactstochastic,and approximatestochasticmodelingproblems.Thearrowsshowprogressionfromsimpletocomplex. and the approximate modeling problem is low-rank approximation is sometimes wellknown(e.g.,inrealizationtheory,modelreduction,andapproximategreatest commondivisor).Inothercases(e.g.,naturallanguageprocessingandconicsection fitting),thelinktolow-rankapproximationislesswellknownandisnotexploited. Theapplicationscanbereadinanyorderorskippedwithoutlossofcontinuity. Commonpatternindatamodeling Applications insystemsandcontrol Themottoofthebookis: Deterministicsystemrealizationandmodelreduction Behindeverydatamodelingproblemthereisa(hidden)low-rankapproxima- tionproblem:themodelimposesrelationsonthedatawhichrenderamatrix Realizationtheoryaddressestheproblemoffindingastaterepresentationofalinear constructedfromexactdatarankdeficient. time-invariantdynamicalsystemdefinedbyatransferfunctionorimpulseresponse representation.Thekeyresultinrealizationtheoryisthatasequence Althoughanexactdatamatrixislowrank,amatrixconstructedfromobserved dataisgenericallyfullrankduetomeasurementnoise,unaccountedeffects,andas- H= H(0),H(1),...,H(t),... sumptionsaboutthedatageneratingsystemthatarenotsatisfiedinpractice.There- isanimpulseresponseofadis(cid:0)crete-timelineartime-inv(cid:1)ariantsystemofordernif fore,generically,theobserveddatadoesnothaveanexactlowcomplexitymodel. andonlyifthetwosidedinfiniteHankelmatrix Thisleadstotheproblemofapproximatemodeling,whichcanbeformulatedasa low-rank approximation problem as follows. Modify the data as little as possible, H(1)H(2)H(3) sothatthematrixconstructedfromthemodifieddatahasaspecifiedlowrank.The ··· modified data matrix being low rank implies that there is an exact model for the H(2)H(3) ...  modifieddata.Thismodelisbydefinitionanapproximatemodelforthegivendata. H(H):=H(3) ... , Thetransitionfromexacttoapproximatemodelingisanimportantstepinbuilding   acoherenttheoryfordatamodelingandisemphasizedinthisbook.  ...    Inallapplications,theexactmodelingproblemisdiscussedbeforethepractically   moreimportantapproximatemodelingproblem.Thisisdonebecause1)exactmod- constructedfromH hasrankn,i.e., elingissimplerthanapproximatemodeling,sothatitistherightstartingplace,and 2)exactmodelingisapartofoptimalapproximatemodelingandsuggestswaysof rank H(H) =orderofaminimalrealizationofH. solvingsuchproblemssuboptimally.Indeed,smallmodificationsofexactmodeling Therefore,existence(cid:0)ofafin(cid:1)itedimensionalrealizationofH (exactlowcomplexity algorithms leadtoeffectiveapproximatemodelingalgorithms. Wellknownexam- lineartime-invariantmodelforH)isequivalenttorankdeficiencyofaHankelma- ples of the transition from exact to approximate modeling in systems theory are trixconstructedfromthedata.Aminimalstaterepresentationcanbeobtainedfrom theprogressionsfromrealizationtheorytomodelreductionandfromdeterministic arankrevealingfactorizationofH(H). subspaceidentificationtoapproximateandstochasticsubspaceidentification. When there is no exact finite dimensional realization of the data or the exact The estimator consistency question in stochastic estimation problems corre- realization is of high order, one may want to find an approximate realization of a spondstoexactdatamodelingbecauseasymptoticallythetruedatageneratingsys- specified low order n. These, respectively, approximate realization and model re- temisrecoveredfrom observeddata.Estimationwithfinitesamplesize,however, ductionproblemsnaturallyleadtoHankelstructuredlow-rankapproximation. necessarily involves approximation. Thus in stochastic estimation theory there is The deterministic system realization and model reduction problems are further alsoastepoftransitionfromexacttoapproximate,seeFigure1.2. consideredinSections2.2,3.1,and4.2. 1.3 Overviewofapplications 9 10 1 Introduction Stochasticsystemrealization w(1) w(2) w(T nmax) ··· − w(2) w(3) w(T n +1) max Let y be the output of an nth order linear time-invariant system, driven by white Hnmax+1(w):= .. .. ··· − ..  (Hi) . . . noise(astochasticsystem)andletEbetheexpectationoperator.Thesequence   w(nmax+1)w(nmax+1) w(T)   ···    R= R(0),R(1),...,R(t),... withn +1blockrows,constructedfromthetrajectoryw,isrankdeficient: max definedby (cid:0) (cid:1) R(τ):=E y(t)y⊤(t τ) rank Hnmax+1(w) ≤rank Hnmax+1(u) +orderofthesystem. (SYSID) − is called the autocorrelation sequence(cid:0)of y. Stochas(cid:1)tic realization theory is con- Conversely(cid:0),iftheHank(cid:1)elmatrix(cid:0)Hnmax+1(w)(cid:1)hasrank(nmax+1)m+nandthema- cerned with the problem of finding a state representation of a stochastic system trix H2nmax+1(u) is full row rank (persistency of excitation of u), then w is a tra- jectory of a controllable linear time-invariant system of order n. Under the above thatcouldhavegeneratedtheobservedoutputy,i.e.,alineartime-invariantsystem assumptions, the data generating system can be identified from a rank revealing drivenbywhitenoise,whoseoutputcorrelationsequenceisequaltoR. An important result in stochastic realization theory is that R is the output cor- factorizationofthematrixHnmax+1(w). When there are measurement errors or the data generating system is not a low relation sequence of an nth order stochastic system if and only if the Hankel ma- trixH(R)constructedfromRhasrankn,i.e., complexity linear time-invariant system, the data matrix Hnmax+1(w) is generi- callyfullrank.Insuchcases,anapproximatelow-complexitylineartime-invariant rank H(R) =orderofaminimalstochasticrealizationofR. modelforwcanbederivedbyfindingaHankelstructuredlow-rankapproximation Therefore,stoch(cid:0)asticre(cid:1)alizationofarandomprocessyisequivalenttodeterministic coafnHbnemaapx+p1li(ewd)a.lTsoheforerfaoprep,rothxeimHaatenkseylstsetmrucidtuernetdifilcoawti-ornan.Lkianpeparrotximime-aitniovnarpiarnotbsleyms- realizationofitsautocorrelationsequenceR.Whenitexists,thefinitedimensional temidentificationisamaintopicofthebookandappearsfrequentlyinthefollowing stochastic realizations can be obtained from a rank revealing factorization of the chapters. matrixH(R). Similarly,totheanalogybetweendeterministicandstochasticsystemrealization, Inpractice,onlyafinitenumberoffinitelengthrealizationsoftheoutputyare thereisananalogybetweendeterministicandstochasticsystemidentification.The available, so that the autocorrelation sequence is estimated from y. With an esti- latteranalogysuggestsanapplicationofHankelstructuredlow-rankapproximation mate R of the autocorrelation R, the Hankel matrix H(R) is almost certainly full tostochasticsystemidentification. rank,whichimpliesthatafinitedimensionalstochasticrealizationcannotbefound. Therefbore,theproblemoffindinganapproximatestochastbicrealizationoccurs.This problemisagainHankelstructuredlow-rankapproximation. Applications incomputeralgebra Systemidentification Greatestcommondivisoroftwopolynomials Realizationtheoryconsidersasystemrepresentationproblem:passfromonerepre- Thegreatestcommondivisorofthepolynomials sentationofasystemtoanother.Alternatively,itcanbeviewedasaspecialexact identificationproblem:findfromimpulseresponsedata(aspecialtrajectoryofthe p(z)=p0+p1z+ +pnzn and q(z)=q0+q1z+ +qmzm ··· ··· system)astatespacerepresentationofthedatageneratingsystem.Theexactidenti- isapolynomialcofmaximaldegreethatdividesbothpandq,i.e.,amaximaldegree ficationproblem(alsocalleddeterministicidentificationproblem)istofindfroma polynomialc,forwhichtherearepolynomialsrands,suchthat generalresponseofasystem,arepresentationofthatsystem.Let p=rc and q=sc. w=col(u,y), where u= u(1),...,u(T) and y= y(1),...,y(T) be an input/output trajectory of a d(cid:0)iscrete-time lin(cid:1)ear time-inva(cid:0)riant system of(cid:1)or- DefinetheSylvestermatrixofthepolynomials pandq dernwithminputsandpoutputsandletn beagivenupperboundonn.Then max theHankelmatrix 1.3 Overviewofapplications 11 12 1 Introduction ℓ1 p0 q0 • p p q q  ...1 p01 ... ...1 q10 ...  ··· ℓm R(p,q):=pn ... ... p0 qm ... ... q0∈R(n+m)×(n+m). (R) •  pn p1 qm q1  ... ... ... ...  w1 ··· wq    pn qm Fig.1.3 Antennaarrayprocessingsetup.   (Byconvention,inthisbook,allmissingentriesinamatrixareassumedtobezeros.) Awellknownfactinalgebraisthatthedegreeofthegreatestcommondivisorof p where τ is the time needed for the wave to travel from the source to the array 1 andqisequaltotherankdeficiency(co-rank)ofR(p,q),i.e., and p Rq is the array’s response to the source emitting at a unit intensity. The 1 ∈ vector p dependsonlyonthearraygeometryandthesourcelocationandisthere- 1 degree(c)=n+m rank R(p,q) . (GCD) − foreconstantintime.Measurementsoftheantennaattimeinstantst=1,...,T give Suppose that p and q havea greatestcommon(cid:0) divisor(cid:1)of degree d>0, but the adatamatrix coefficientsofthepolynomialspandqareimprecise,resultinginperturbedpolyno- D:= w(1) w(T) =p1 ℓ1(1 τ) ℓ1(T τ) =p1ℓ1, mials pd andqd.Generically,thematrixR(pd,qd),constructedfromtheperturbed ··· − ··· − polynomials,isfullrank,implyingthatthegreatestcommondivisorofpdandqdhas (cid:2) (cid:3) (cid:2) ℓ1 (cid:3) degreezero.Theproblemoffindinganapproximatecommondivisor of p andq d d whichhasrankequaltoone. | {z } withdegreed,canbeformulatedasfollows.Modifythecoefficientsof p andq , aslittleaspossible,sothattheresultingpolynomials,say, pandqhaveadgreatedst Consider now m<q distant sources emitting with intensities ℓ1,...,ℓm. Let pk be the response of the array to the kth source emitting alone with unit intensity. common divisor of degree d. This problem is a Sylvester structured low-rank ap- Assumingthatthearrayrespondslinearlytoamixtureofsources,wehave proximationproblem.Therefore,Sylvesterstructuredlow-rankapproximationcan b b beappliedforcomputinganapproximatecommondivisorwithaspecifieddegree. m Theapproximategreatestcommondivisorcfortheperturbedpolynomialspdandqd D= w(1)···w(T) = ∑pk ℓk(1−τk)···ℓk(T−τk) =PL, istheexactgreatestcommondivisorof pandq. k=1 (cid:2) (cid:3) (cid:2) ℓk (cid:3) TheapproximategreatestcommondivisorproblemisconsideredinSection3.2. b | {z } b b whereP:= p1 pm ,L:=col(ℓ1,...,ℓm),andτkisthedelayofthewavecoming ··· fromthekthsource.ThisshowsthattherankofDislessthanorequaltothenumber (cid:2) (cid:3) of sources m. If the number of sources m is less than the number of antennas q Applications insignalprocessing and m is less than the number of samples T, the sources intensities ℓ ,...,ℓ are 1 m linearly independent, and the unit intensity array patterns p ,...,p are linearly 1 m Arraysignalprocessing independent,thenwehavethat Anarrayofantennasorsensorsisusedfordirectionofarrivalestimationandadap- rank(D)=thenumberofsourcestransmittingtothearray. tivebeamforming.Considerqantennasinafixedconfigurationandawavepropa- gatingfromdistantsources,seeFigure1.3. Moreover,thefactorsPandLinarankrevealingfactorizationPLofDcarryinfor- Consider,first,thecaseofasinglesource.Thesourceintensityℓ (thesignal)is mationaboutthesourcelocations. 1 afunctionoftime.Letw(t) Rqbetheresponseofthearrayattimet(w beingthe With noisy observations, the matrix D is generically a full rank matrix. Then, i responseoftheithantenna).∈Assumingthatthesourceisfarfromthearray(relative assumingthatthearray’sgeometryisknown,low-rankapproximationcanbeused tothearray’slength),thearray’sresponseisproportionaltothesourceintensity toestimatethenumberofsourcesandtheirlocations. w(t)=p ℓ (t τ), 1 1 1 −