A FAST ALGORITHM FOR MATRIX BALANCING

AFASTALGORITHMFORMATRIXBALANCINGPHILIPA.KNIGHT∗ANDDANIELRUIZ†Abstract.AslongasasquarenonnegativematrixAcontainssufficientnonzeroelements,thenthematrixcanbebalanced,thatiswecanfindadiagonalscalingofAthatisdoublystochastic.Anum-berofalgorithmshavebeenproposedtoachievethebalancing,themostwellknownofthesebeingtheSinkhorn-Knoppalgorithm.Inthispaperwederivenewalgorithmsbasedoninner-outeriterationschemes.WeshowthattheSinkhorn-Knoppalgorithmbelongstothisfamily,butothermemberscanconvergemuchmorequickly.Inparticular,weshowthatwhilestationaryiterativemethodsofferlittleornoimprovementinmanycases,aschemeusingapreconditionedconjugategradientmethodastheinneriterationcangivequadraticconvergenceatlowcost.Keywords.Matrixbalancing,Sinkhorn-Knoppalgorithm,doublystochasticmatrix,conjugategradi-entiteration.AMSsubjectclassifications.15A48,15A51,65F10,65H10.1.Introduction.Foratleast70years,scientistsinawidevarietyofdisciplineshaveattemptedtotransformsquarenonnegativematricesintodoublystochasticformbyapplyingdiagonalscalings.Thatis,givenA∈Rn×n,A≥0,finddi-agonalmatricesD1andD2sothatP=D1AD2isdoublystochastic.Motivationforachievingthisbalanceincludeinterpretingeconomicdata[2],preconditioningsparsematrices[13],understandingtrafficcirculation[11]andorderingnodesinagraph[9].Inalloftheseapplications,oneofthemainmethodsconsideredistheSinkhorn-Knoppalgorithm1.ThisisaniterativeprocessthatattemptstofindD1andD2byalternatelynormalisingcolumnsandrowsinasequenceofmatricesstart-ingwithA.Convergenceconditionsforthisalgorithmarewellknown:ifAhasfullsupport2thenthealgorithmwillconvergelinearlywithasymptoticrateequaltothesquareofthesubdominantsingularvalueofP[18,19,9].Clearly,insomecasestheconvergencewillbepainfullyslow.Theprincipalaimofthispaperistoderivesomenewalgorithmsforthematrixbalancingproblemwithaneyeonspeed,especiallyforlargesystems.FirstwelookatasimpleNewtonmethodforsymmetricmatrices,closelyrelatedtoamethodproposedbyKhachiyanandKalantari[8]forpositivedefinite(butnotnecessarilynonnegative)matrices.We∗DepartmentofMathematics,UniversityofStrathclyde,26RichmondStreet,GlasgowG11XHScot-land([e-mail]).ThisworkwassupportedinpartbyagrantfromtheCarnegieTrustfortheUniversitiesofScotland.†ENSEEIHT2,rueCharlesCamichel,Toulouse,France([e-mail])1Thealgorithmhasbeengivenmanydifferentnames,Sinkhorn-Knoppisthenameusuallyadoptedbythelinearalgebracommunity2Fullsupportisdefinedin§3.1DagstuhlSeminarProceedings07071WebInformationRetrievalandLinearAlgebraAlgorithms’smethodproducesasequenceofpositiveiterates,theJacobianswegeneratewillbepositivesemi-definite.ToapplyNewton’smethodexactlywerequirealinearsystemsolveateachstep,andthisisusuallyprohibitivelyexpensive.Wethereforelookatiterativetechniquesforapproximatingthesolutionateachstep.FirstwelookatsplittingmethodsandweseethattheSinkhorn-Knoppalgorithmisamemberofthisfamilyofmethods.Wegiveanasymptoticboundonthe(linear)rateofconvergenceofthesemethods.Nextwelookatapreconditionedconjugategradientmethodforsolvingtheinneriteration.Wediscussimplementationdetailsandshowthatasymptotically,super-linearconvergenceisachievable.ByimplementinganArmijo-typeruleweensurethatouriteratesretainpositivityandwedemonstratethereliabilityandspeedofourmethodintests.Anumberofauthorshavepresentedalternativetechniquesforbalancingma-tricesthatcanconvergefasterthantheSinkhorn-Knoppalgorithm.Forexample,ParlettandLandis[16]lookatsomesimplewaysoftryingtoacceleratethecon-vergencebyfocusingonreducingstatisticssuchasthestandarddeviationbetweenrowsumsoftheiterates.Incertaincases,theyshowgreatimprovementispossible,suggestingthattherateonconvergenceisnotdependentonthesingularvaluesofP;buttheyalsogiveexampleswheretheiralternativesperformsignificantlyworse.Weincludeacomparisonofoneofthesealgorithmsagainstourproposedapproachin§6.Linialetal.[12]useasimilarapproachto[16]inthecontextofestimatingmatrixpermanents,althoughtheirupperboundoniterationcounts(O(n7))makesthemdistinctlyunappealingforlargeproblems!Acompletelydifferentapproachistoviewmatrixbalancingasanoptimisationproblem.Therearemanyalternativeformulations,perhapsthefirstbeingthatofMarshallandOlkhin[14]whoshowthatbalancingisequivalenttominimizingthebilinearformxTAysubjecttotheconstraintsΠxi=Πyi=1.However,theexper-imentalresultswehaveseenforoptimisationtechniquesforbalancing[14,17,3]suggestthatthesemethodsarenotparticularlycheaptoimplement.2.Newton’smethod.LetD:Rn→Rn×nrepresenttheoperatorthatturnsavectorintoadiagonalmatrix,D(x)=diag(x),andleterepresentavectorofones.Thentobalanceanonnegativematrix,A,weneedtofindpositivevectorsrandcsuchthatD(r)AD(c)e=D(r)Ac=e,D(c)ATr=e.(2.1)Rearrangingtheseidentitiesgivesc=D(ATr)−1e,r=D(Ac)−1e,FASTMATRIXBALANCING3andonewayofwritingtheSinkhorn-Knoppalgorithm[6,9]isasck+1=D(ATrk)−1e,rk+1=D(Ack+1)−1e.(2.2)IfAissymmetricthen(2.1)canbesimplified.Toachievebalancingweneedavectorxthatsatisfiesf(x)=D(x)Ax−e=0.(2.3)Thisleadstotheiterativestepxk+1=D(Axk)−1e.(2.4)Exactlythes