A centroid approach to the joint approximate diagonalisation problem - PDF document

J.H.Manton/DigitalSignalProcessing16(2006)468–478,wheretherepresentrandomerrormatrices,cannotbefoundbyastraightforwardeigendecomposition.Iftheerrordistributionsareknown,themaximumlikelihoodestimateofisobtainedbysolvingtheconstrainedoptimisationproblem:Given,ndtheX,Cwhichmaximisethejointlikelihoodof,where,subjecttotheequalityconstraintsthatthearediagonalmatricesand.Thisiscalledjointapproximatediagonalisation,emphasisingthatonlyapproximatelydiagonalisestheobservedmatrices.Differentdistributionsforthe,andestimatorsotherthanthemaximum-likelihoodestimator,leadtodifferentdenitionsofanoptimaljointapprox-imatediagonaliser.Allalgorithmsforcomputingthesevariousdiagonalisersiterativelyminimiseacostfunctionandmayconvergetoalocalminimum,hencetheirperformancecannotbeguaranteed[3,6,7].Thispaperproposesanewdenitionofjointapproximatediagonalisationwhichisgeo-metricallymeaningful,irrespectiveofanyparticularsignalprocessingapplication.Thisgeometricalinsightshedslightonthestrengthsandlimitationsofjointapproximatediag-onalisationtechniquesingeneral.Forinstance,justlikeantipodalpointsonaspheredonothaveauniquecentroid,incertaincasesthejointapproximatediagonalisercannotbedeneduniquelyeither.Moreover,theproposeddenitionischosentoensuretheoptimaldiagonaliserwithrespecttothisdenitioncanbecomputedreliablyinpractice.Althoughtheproposedalgorithmcomputesadiagonaliserwhichisoptimalaccordingtothegeometricdenition,thediagonalisermaynotbeoptimalforagivensignalprocess-ingapplication.Nevertheless,insuchsituations,theintentionisfortheproposedalgorithmeithertobeuseddirectlyasasub-optimalbutreliablesolution—reliablealbeitsub-optimalalgorithmsareoftenpreferable—orasameansofcomputingagoodinitialguessforasub-sequentlocaloptimisationroutine,suchasifthemaximum-likelihoodestimateissought.Itisemphasisedtheideasinthispaperareintheirinfancyand,asnotedthroughoutthepaper,therearenumerouswaystheresultscanbeeithergeneralisedortailoredtospe-cicapplications.Therstofsuchnotesisthatthecomplex-valuedcasewillbestudiedelsewhere;whilethegeometricdenitionofjointapproximatediagonalisationeasilygen-eralises,differenttechniquesarerequiredforamathematicalanalysis.2.GeometricjointapproximatediagonalisationThejointapproximatediagonalisationproblemistondanorthonormalmatrixsuchthattheareapproximatelydiagonal,wherethearesymmetricma-trices.Therearenumerouswaysofdeningwhatapproximatelydiagonalmeans.Thissectionintroducesageometricdenition,themotivationbeingthatunderthisdenition,theoptimaldiagonalisercanbefoundreliablyinpracticebyminimisingaconvexfunctiononamanifold.Intuitively,previousdenitionsofoptimaljointdiagonalisationledtodifcultoptimi-sationproblemsbecausethejointdiagonalisationproblemisageneralisationofthematrixeigenvalueproblem,andtodate,attemptingtosolvethematrixeigenvalueproblembyminimisingacostfunction(withoutanyformofimplicitorexplicitdeationtakingplace)isnotcompetitivewithcurrentbestalgorithms.Thissuggeststheinherentmatrixeigen- J.H.Manton/DigitalSignalProcessing16(2006)468–478onit,inducedfromtheuniquebi-invariantmetric[9,10],whichdoesindeedmeasuretheamountofrotation:(X,Y)  ,X,YwherelogistheprincipalmatrixlogarithmandtheFrobeniusnorm.bethesub-groupofsuchthatifandonlyifthereexistsapermuta-tionmatrixandadiagonalmatrix.Totakeintoaccountthecorrectpermutationandsignchange,dened(X,Y)P,Q(XP,YQ)(X,YQ).DeÞnition1.Givensymmetric,letbesuchthatisdiagonal.Ageometricjointapproximatediagonaliseroftheisanyimising(1),withdenedin(3).Section3usesthegeometryoftoestablishconditionsforthegeometricjointapproximatediagonalisertobeunique.2.1.OrderedjointapproximatediagonalisationIfanhasamultipleeigenvaluethenDenition1isapparentlyunsatisfactorybecausethenon-uniquenessofextendsbeyondpermutationsandsignchangesofitscolumns.Themathematicalsolutionisstraightforward;includethisextraambiguitywhensearchingfortheminimumin(3).Inpractice,thisextralevelofdifcultyisnotwarranted.First,ifarenoise-corruptedversionsofthematrices,thentypicallytheprobabilitythatanhasamultipleeigenvalueiszero.Moreimportantly,thejointdiagonalisationproblemisill-conditionedifanhasamultipleeigenvalue.Forexample,ifarbitrary,arbitrarilysmallperturbationsofcanleadtoitseigenvectorsbeingfarawayfromtheeigenvectorsofThus,ajointdiagonalisationalgorithmdetectinganwithtwocloseeigenvaluesshouldsignalthattheproblemisill-conditioned.Furthermore,alljointdiagonalisationroutinesshouldbeusedcautiouslyinapplicationswherethenoisefreematricescontaincloselyspacedeigenvalues.Thislineofthinkingmotivatesseveralnewgeometricproblemswhichshouldbeofinteresttothesignalprocessingcommunity.TherstistondajointdiagonaliserwheretheorderofthecolumnsofDeÞnition2.Givensymmetric,letbesuchthatthethcolumnofistheeigenvectorassociatedwiththethlargesteigenvalueof.Anorderedjointapproximatediagonaliseroftheisanyminimising(1),within(3)butwhereisinsteadthesub-groupofalldiagonalmatricesin Promptedbyareviewer,wepointoutthatill-conditionedisdifferentfromill-dened;ourclaimdoesnotcontradictTheorem3of[2]. J.H.Manton/DigitalSignalProcessing16(2006)468–478 Fig.3.Thedottedlinesareageodesictriangleconnectingthreepoints(projectedontoatwo-dimensionalplane).Thesolidlinesaregeodesicsjoiningthemidpointofeachsideofthetriangletotheoppositecorner.Thecurvatureofpreventsthethreesolidlinesmeetingatauniquepoint.Althoughtheprojectionismany-to-one,locallyitisone-to-one.Thatis,thegeometry(curvature,geodesics,etc.)ofcoincideslocallywiththegeometryof,andthelatteriswellunderstood:equiptwithitsbi-invariantmetricisacompactRiemannianmanifoldwithnon-negativecurvature.Geodesicsonoftheform(t),where,whereisthesetofskew-symmetricmatrices.,andhencehaspositivecurvature,ifthearetoofarapartthecentroidneednotbeunique.Tovisualisethis,considerthenorthandsouthpolesonasphereandnoteanypointontheequatorisequidistantfromthem.Inpractice,thismeansthatiftoomuchnoiseispresent,sothearefarfrombeingjointlydiagonalisable,thenthejointapproximatediagonaliserisill-dened.Sincealljointdiagonalisationalgorithmsminimiseacostfunctionon(oraquotientspacethereof),thisisaninherentfea-tureofthejointapproximatediagonalisationproblemandnotofanyspecicdenitionorItisremarkedthatalthoughthereareseveraldifferentdenitionsofcentroidwhichturnouttobeequivalentinEuclideanspace,onthisisnolongertrue.See,forexample,Fig.3.Thedenitionin(1)ofacentroidistakenfrom[11]andhasbeenstudiedinthespecialcaseofin[12].Thegloballyconvergentalgorithmin[13]forcomputingthiscentroidformspartofthealgorithminSection4.Furtherdetailsofthegeometryofcanbefoundinthesepapers.Thelargestdomainonwhichisinjectiveisthefundamentaldomain(X,I)d(X,Q)DenotetheopenballinB(X(X,Y) J.H.Manton/DigitalSignalProcessing16(2006)468–478Lemma3..TheopenballB(IiscontainedinthefundamentaldomainifandonlyifProof. (Q,I).ClearlyB(Iifandonlyifbethematrixwithallelementszeroexceptfor3.Then(Q,I)2,proving4.Theproofthat4istediousandhenceomitted.Itisintuitivethoughthatistheelementinclosesttobecauseitdiffersfrominonlytwocolumns.Theorem4..Referringto,ifthereexistsasuchthatB(Ythenthereispreciselyoneminimumoff(X)intheregionB(Y.More-over,restrictedtoB(Yisstrictlyconvex,andisthelargestradiusforwhichthisistrue.Proof.f(X)(X,X.IfB(Ythenitisprovedin[11,12]restrictedtoB(YisstrictlyconvexandhencehasauniqueminimumB(Y,whichispreciselytheKarchermeanofthepoints.(Infact,f(X)isconvexonamuchlargerball.)Now,Lemma3impliesthatif(X,X)/4thend(X,X(X,X,soifX,XB(Yd(X,X)/4,henceagreeonB(Y,thusprovingrestrictedtoB(Yisstrictlyconvexwithmini-Toprove8isthelargestradius,deneQ()tobetheidentitymatrixexceptforthetopleft22blockwhichisscos,Šsin;sin,cos],arotationmatrix.NoteQ(/.LetQ(/.Sinceisequivalentto)Q(/Q()[Š,///8,3/8](5)andclearlyf(X)isnotconvexatQ(/Admittedly,theradius8issmall;ifamoderateamountofnoiseispresent,itmightnotbepossibletochoosesuchthatB(Y,eventhoughthegeometricjointdiagonaliseriswelldenedandcorrectlycomputedbythealgorithminSection4.TheproblemisprimarilycausedbythevolumeofB(Ibeingconsiderablysmallerthanthevolumeofeventhoughcontainsnolargerball.Thewaytoovercomethisistouseadifferentmetricfordeningopenballs.Thisiscurrentlyunderinvestigation.4.Thealgorithm,anumericalexample,andadiscussionAnalgorithmforcomputingthegeometricjointdiagonaliser(Denition1)isnowdescribed.Givensymmetric,nd,byeigen-decomposition,suchthatthearediagonal.Ifanyofthehavemulti-pleeigenvalues,warnaboutill-conditioning.Multiplytherstcolumnofeach1ifnecessary,sodet1.Recallingthedenitionofin(3),compute J.H.Manton/DigitalSignalProcessing16(2006)468–478argminbyexhaustivesearch.(Thecardinalityof,sothisisonlypracticalforsmall.Fastersolutionsnotrequiringanexhaustivesearcharecurrentlyunderinvestigation.Seetoothenotebelow.)Set2thenreturnasthejointdiagonaliser.Otherwise,foreach,computeargmin(X,Xandset.(Thisensurestheareclustered.)Ifthereexistsansuchthat,X)8,warnthatreliabilitycan-notbeguaranteed.Repeatthesteps:(1)Set .Toreducenumericalroundofferrors,set(2)Ifthenreturnasthejointdiagonaliser.(3)Setandgotostep1.(Thesethreestepsimplementthesteepestdescentalgorithmproposedin[13].)Notethattheexhaustivesearchisnotnecessaryifthearesufcientlyclosetogether,suchaswhentheareclosetobeingexactlydiagonalisableandhavewellseparatedeigenvalues.Inthiscase,ifthethcolumnofistheeigenvectorofatedwiththethlargesteigenvalue,thenminachievesitsminimumisdiagonalwith1iftheinnerproductofthethcolumnsofpositive,otherwise1.Similarly,computingtheorderedjointdiagonaliser(Den-ition2)eliminatestheneedtosearchthroughallpermutationmatrices.ThealgorithmwasimplementedandcomparedwithCardoso’swell-knownJacobimethodforcomputingjointapproximatediagonalisers[7].SinceCardoso’smethodmin-imisesadifferentcostfunction(namely,thesumofthesquaresoftheoffdiagonalele-ments),whetherCardoso’smethodortheproposedmethodperformsbettershoulddependonthenoisemodel.Twodifferentnoisemodelsareused.TherstisUDUisanarbitraryorthonormalmatrix,wheretheelementsofhaveaGaussiandistribution.ItisanticipatedCardoso’smethodwilloutperformthegeometricjointdiagonaliserbecausethelatteres-sentiallyassumeseacheigenvectorofisperturbedbyapproximatelythesameamountfromthecorrespondingeigenvectorofUDU,whereastheadditivenoisewillperturbsomeeigenvectorsmorethanothersbecausetheeigenvaluesrangefrom1to5.Figure4illustratesthiswhen5.Here,ifisthejointapproximatediagonaliser,theerrorisgraphedasd(U,V).Notethat5islargeenoughforthe8distanceruletobeviolatedfrequently,yetthealgorithmstillworkswell.Notetoothatresultsnotpresentedindicatethatif2thentheperformancedifferencebetweenthetwomethodsisneg-ligible.(Whensimulatingthisandthefollowingmodel,50trialswereperformed,withthegeneratedatrandomforeachtrial.Eachpointintheguresrepresentsonetrial.)ThesecondmodelisD(e4,whereisanarbitraryorthonormalmatrix,wheretheelementsofhaveaGaussiandistribution.Sincerepresentsauniformrandomperturba-tionofthecolumnsof,itisanticipatedthegeometricjointdiagonaliserwilloutperformCardoso’smethod.Figure5showsthatwhen5,thisisindeedthecase;thegeo- J.H.Manton/DigitalSignalProcessing16(2006)468–478withalargereigenvalueshouldincuralargerpenaltysinceadditivenoiseislesslikelytoaffecttheseeigenvectors.However,anad-hocchangeto(3)caneasilydestroytheconvex-ityofthecostfunctionf(X)in(1).Initsfullgenerality,thekeyideainthispaperistotobeadistancefunctiononaRiemannianmanifold,sincethiswillensuref(X)isconvexprovidedthearesufcientlyclosetogether.Theonlyotherrestrictionbeleft-invariant,thatisd(X,Y)d(UX,UY)forall,sincethisensuresthegeometricjointdiagonaliserisequivarianttoorthonormalchangesofbasisofthe(requirementR3inSection3).denotethetangentspaceattheidentityof.(ItistheLiealgebraassociatedandconsistsofskew-symmetricmatrices.)Aleft-invariantdistancefunctiononisfoundbyassigninganinnerproducttoandthenextendingittothewholetangentbundleofbylefttranslation.Fornoisemodelone,asuitablechoiceofinnerproductisfoundasfollows.If(sothathatD,AAD,Aisdiagonal.Providedhasdistincteigenvalues,givenasymmetricthereexistsansuchthatthatD,A.Thus,(6)saystorstorderthatadiagonalmatrixperturbedbyhasitseigenvectorsperturbedfrom,wherehereD,A.Aleast-squaresestimatorseekstominimisetrtrD,A.ThissuggestsdeningthenormoftobebeD,A,andinfact,thisnormcomesfromtheinnerproductA,BBD,AAD,B,A,BAsalreadymentioned,(7)inducesaleft-invariantdistancefunction,andgeometricjointdiagonalisationwithrespecttothisdistancefunctionisexpectedtoperformwellfornoisemodelone.Beforethiscanbeveriedthough,thegeometryinducedby(7)mustbede-rived,andthisisbeyondthescopeofthepresentpaper.5.ConclusionAjointapproximatediagonalisationalgorithmshouldfulllrequirementsR1–R4inSection3andgivereliableresults.Toachievethis,thispaperproposedtofactoroutthedifcultmatrixeigenvalueproblemfromthejointdiagonalisationproblemandthenusegeometrytoconstructaconvexcostfunctiononaRiemannianmanifoldwhoseminimumdenesthejointapproximatediagonaliser.Assuch,thisistherstjointdiagonalisationalgorithmtouseaconvexcostfunctionandhencenotsufferfromconvergencetoalocalAcknowledgmentsTheauthoracknowledgesthesupportoftheAustralianResearchCouncilandtheARCSpecialResearchCentreforUltra-BroadbandInformationNetworks(CUBIN).