clique problem

Annealedreplication:anewheuristicforthemaximum

cliqueproblem

ImmanuelM.Bomzea;∗,MarcoBudinichb,MarcelloPelilloc,

ClaudioRossic

aInstitut

f󰁄urStatistik&DecisionSupportSystems,Universit󰁄atWien,Universit󰁄atsstra󰀃e5,

A-1010Wien,Austria

bDipartimentodiFisica,UniversitÂadiTrieste,ViaValerio2,I-34127Trieste,ItalycDipartimentodiInformatica,UniversitÂaCa’FoscaridiVenezia,ViaTorino155,

I-30172VeneziaMestre,Italy

Received10April1998;receivedinrevisedform9April2001;accepted8May2001

Abstract

Weproposeanewheuristicforapproximatingthemaximumcliqueproblembasedonadetailedanalysisofaclassofcontinuousoptimizationmodelswhichprovideacompletecharac-terizationofsolutionstothisNP-hardcombinatorialproblem.Westartfromaknowncontinuousformulationofthemaximumclique,andtacklethesearchforlocalsolutionswithreplicatordy-namics,aclassofdynamicalsystemsdevelopedinvariousbranchesofmathematicalbiology.Hereby,weaddtotheobjectiveusedinpreviousworksaregularizationtermthatcontrolstheglobalshapeoftheenergylandscape,thatisthefunctionactuallymaximizedbythedynamics.Theparametercontrollingtheregularizationischangedduringtheevolutionofthedynamicalsystemtorenderine󰀄cientlocalsolutions(whichformerlywerestable)unstable,thusconduct-ingthesystemtoescapefromsub-optimalpoints,andsotoimprovetheÿnalresults.Theroleofthisparameteristhussuperÿciallysimilartothatoftemperatureinsimulatedannealinginthesensethatitsvariationallowstoÿndbettersolutionsfortheproblemathand.Wedemonstrateseveraltheoreticalresultsontheregularizationtermandwefurthersupportthevalidityofthisapproach,reportingonitsperformanceswhenappliedtoselectedDIMACSbenchmarkgraphs.?2002ElsevierScienceB.V.Allrightsreserved.

Keywords:Quadraticoptimization;Dynamicalsystems;Replicatordynamics;Stableset;Independentset;Combinatorialoptimization

1.Introduction

Themaximumcliqueproblem(MCP)isawell-knownproblemincombinatorialoptimizationwhichÿndsimportantapplicationsinmanydi󰀁erentdomains[11].Since

Correspondingauthor.

E-mailaddresses:immanuel.bomze@univie.ac.at(I.M.Bomze),mbh@trieste.infn.it(M.Budinich),pelillo@dsi.unive.it(M.Pelillo),rossi@dsi.unive.it(C.Rossi).

0166-218X/02/$-seefrontmatter?2002ElsevierScienceB.V.Allrightsreserved.PII:S0166-218X(01)00233-5

∗

28I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

theMCPisknowntobeNP-hard,exactalgorithmsareguaranteedtoreturnasolutiononlyinatimewhichincreasesexponentiallywiththenumberofverticesinthegraph.Thismakestheminapplicableeventomoderatelylargeprobleminstances.Moreover,aseriesofrecenttheoreticalresultsshowthattheMCPis,infact,di󰀄culttosolveevenintermsofapproximation.Strongevidenceofthisfactcamein1991,whenFeigeetal.[17](seealso[18])provedthatifthereisapolynomial-timealgorithm

1−󰀅

thatapproximatestheMCPwithinafactorof2logn,thenanyNPproblemcanbe

O(1)

solvedin“quasi-polynomial”time(i.e.,in2logntime).TheresultwasfurtherreÿnedbyAroraetal.[1,2]oneyearlater.Speciÿcally,theyprovedthatthereexistsan󰀅¿0suchthatnopolynomial-timealgorithmcanapproximatethesizeofthemaximumcliquewithinafactorofn󰀅,unlessP=NP.Morerecentdevelopmentsalongtheselinescanbefoundin[3,4,23].Inlightofthesenegativeresults,muche󰀁orthasrecentlybeendirectedtowardsdevisinge󰀄cientheuristicsfortheMCP,forwhichnoformalguaranteeofperformancemaybeprovided,butareanywayofinterestinpracticalapplications.Wereferto[25]foracollectionofpromisingheuristicsfortheMCP.Wehaverecentlyinvestigatedthee󰀁ectivenessofanapproachforapproximatingtheMCP,centeredaroundacontinuousformulationduetoMotzkinandStraus[33]anditsregularization[24,10],whichexploitsthedynamicalpropertiesoftheso-calledreplicatorequations,aclassofdynamicalsystemsdevelopedandstudiedinvariousbranchesofmathematicalbiology.Oneproblemassociatedwiththesemodels,how-ever,istheirinabilitytoescapeine󰀄cientlocalsolutions.Inthispaper,weintroduceaclassofparametrizedquadraticprograms,whichincludesboththeMotzkin–Strausprogramanditsregularizationasspecialcases,andinvestigatethepropertiesofitsso-lutionsasafunctionofitsparameter.AdetailedanalysisofthesepropertiessuggestsanewalgorithmforapproximatingtheMCPwhichisbasedontheideaofproperlyvaryingtheparameterduringthereplicatoroptimizationprocess,inmuchthesamespiritassimulatedannealingprocedures.Arelated,butdi󰀁erent,ideahasrecentlybeenproposedbyGeeandPragerintheneuralnetworkdomain[20].ExperimentalresultsconductedonvariousDIMACSbenchmarkgraphsdemonstratethevalidityoftheproposedapproach.

Theoutlineofthepaperisasfollows.InSection2,wedescribetheMotzkin–Straustheoremanditsparameterization,andpresentthereplicatordynamicalsystems.ThesedynamicsareusedtoobtainlocallyoptimalsolutionstotheMCP.Section3isde-votedtoafewresultsthatenableustoestablishboundsonaregularizationparameter󰀁whichgovernsstabilityunderthereplicatordynamics.Forillustration,weinvestigateinSection4asmall,butprototypicalexampleindetail.Inamoredetaileddynamicalanalysisexceedingtheusualperturbationtheoryapproach,wespecifyexplicitrangeswithinwhichqualitativefeaturesofthedynamicsareinvariant,andalsoobtainquanti-tativesensitivityresultsfortherelatedoptimizationproblems.Thisanalysisisdeferredtoanappendix,topromotethe󰀃owoftheargument.Thepreviouslyestablishedtheo-reticalpropertieswillleadustodevelopinSection5analgorithmforproperlyupdatingtheparameter󰀁withtheobjectiveofavoidingpoorlocalsolutions.InSection6theresultsofourexperimentsarepresented,andSection7concludesthepaper.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4929

2.Evolutiontowardsthemaximumclique

2.1.AparametrizedcontinuousformulationoftheMCP

LetG=(V;E)beanundirectedgraph,whereV={1;:::;n}isthesetofvertices,andE⊆V×Visthesetofedges.AcliqueofGisasubsetofVinwhicheverypairofverticesisconnectedbyanedge.AcliqueiscalledmaximalifnostrictsupersetofCisaclique,i.e.,novertexexternaltoCisconnectedwithmorethan|C|−1verticesofC(here,andinthesequel,|C|denotesthecardinalityofasetC).AmaximalcliqueCiscalledstrictlymaximalifnovertexiexternaltoChasthepropertythattheenlargedsetC∪{i}containsacliqueofthesamesizeasC.Inotherwords,if

dC(i)=|{j∈C:(i;j)∈E}|

denotesthedegreeofiw.r.t.C,thenamaximalcliqueCisstrictlymaximalifandonlyifdC(i)¡|C|−1foralli∈C.

Amaximumcliqueisacliquehavinglargestcardinality(notethatamaximalcliqueisnotnecessarilyamaximumone).Hence,theMCPconsistsofÿndingacliqueofmaximumsizeinagraphG.Forarecentsurveysee[11].Inthefollowing,givenasetSofverticesinG,wewilldenotebyxSitscharacteristicvector,deÿnedasxiS=1=|S|ifi∈SandxiS=0otherwise.

GivenagraphG,considerthefollowingquadraticprogramintroducedin[24,10](x󰀃alwaysdenotesthetransposeofacolumnvectorx):

maximizex󰀃(AG+12I)xsubjecttox∈Sn;

(1)

whereAG=(aij)istheadjacencymatrixofG(i.e.,aij=1if(i;j)∈E,andaij=0if(i;j)∈E),SnisthestandardsimplexofRn,thatis

󰀅󰀆n󰀈

xi=1Sn=x∈Rn:xi¿0foralli=1;:::;nand

i=1

andIisthen×nidentitymatrix.Thisturnsouttobeavariantoftheso-called

1Motzkin–Strausprogram[33],whichisobtainedfrom(1)bysimplydroppingthe2Iterm.Forcompleteness,wesummarizeheretheoriginalMotzkin–Straustheoremandsomerecentrelatedresults.1

Theorem1.LetCbeasubsetofverticesofagraphG;andletxCbeitscharacteristicvector.Then(xC)󰀃AG(xC)=1−1=|C|ifandonlyifCisaclique.Moreover:

(a)xCisastrictlocalmaximizerofx󰀃AGxoverSnifandonlyifCisastrictlymaximalclique.

(b)xCisaglobalmaximizerofx󰀃AGxoverSnifandonlyifCisamaximumclique.

TheoriginalMotzkin–Straustheorem[33]correspondstothe“if”partofTheorem1(b),whilethe“only–if”parthasbeenprovenin[39].Part(a)ofthetheoremiffrom[39,22].

1

30I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

AnimmediateconsequenceofthepreviousresultisthatanypointinSnprovidesuswithaboundonthesizeofthemaximumcliqueinG.2Infact,ifCisamaximumcliqueofG,foranyx∈Snwehavex󰀃AGx61−1=|C|,fromwhichitfollowsthat|C|¿󰀐1=(1−x󰀃AGx)󰀑.

TheMotzkin–Straustheoremhasanintriguingcomputationalsigniÿcance.Itsuggestsafundamentallynewwayofsolvingthemaximumcliqueproblem,byallowingustoshiftfromthediscretetothecontinuousdomain.Apointedoutin[35],theadvantagesofsuchareformulationaremanifold.Itnotonlyallowsustoexploitthefullarsenalofcontinuousoptimizationtechniques,therebyleadingtothedevelopmentofnewe󰀄cientalgorithms,butmayalsorevealunexpectedtheoreticalproperties.Additionally,continuousoptimizationmethodsaresometimesdescribedintermsofsetsofdi󰀁erentialequations,andarethereforepotentiallyimplementableinanalogcircuitry.TheMotzkin–Strausandrelatedtheoremshaveservedasthebasisofmanyclique-ÿndingprocedures[36,37,21,11],andhavealsobeenusedtodeterminetheoreticalboundsonthemaximumcliquesize[15].

IncontrasttotheoriginalMotzkin–Strausformulation,however,itsregularization(1)hasafurthermerit:asobservedbyPardalosandPhillips[36]andlaterformalizedbyPelilloandJagota[39],theMotzkin–Strausprogram,initsoriginalformulation,isplaguedbythepresenceof“spurious”solutions,i.e.,solutionswhicharenotintheformofcharacteristicvectors.Clearly,thisrepresentsaproblemsinceitprohibitsdirectextractionoftheverticescomprisingtheclique,andprovidesinformationonlyonitssize.Therefore,inordertodeterminethecliquevertices,onehastomakerecoursetoiterativeorrecursiveprocedure,asthosedescribedin[34,36].

Thesigniÿcanceofthefollowingresult,asharpeningofTheorem1provedin[10],isthatalocal(andhencealsoaglobal)maximumof(1)canonlybeat-tainedatacharacteristicvectorx∗=xCforsomesubsetCofverticeswhichnec-essarilythenformsamaximalclique.Thissolvesthespurioussolutionprobleminastraightforwardanddeÿnitivemannersinceitestablishesaone-to-onecorrespon-dencebetweenlocal=globalsolutionsto(1)andmaximal=maximumcliquesofG,respectively.

Theorem2.LetGbeagraphandconsiderproblem(1).Thenthefollowingassertionsareequivalent:

(a)x=xC;whereCisamaximalcliqueofsizek=|C|;(b)xisastrictlocalsolutionto(1);(c)xisalocalsolutionto(1).

Ifoneoftheaboveconditions(andthereforeall)ismet;theobjectiveisx󰀃(AG+1C

2I)x=1−1=2k.HenceCisamaximumcliqueofGifandonlyifxisaglobalsolutionto(1).

2

WethankArunJagotaforpointingthisout.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4931

Inthispaper,weconsiderthefollowingprogram,whichrepresentsalsoaregular-izationoftheMotzkin–Strausprogramandgeneralizes(1):

maximizef󰀁(x)=x󰀃(AG+󰀁I)xsubjecttox∈S:

n

(2)

ThisincludesboththeMotzkin–Straus(󰀁=0)programanditsregularization(󰀁=12)asspecialcases.Weinvestigatethepropertiesofitssolutionsasafunctionoftheparameter󰀁.Speciÿcally,weshowthatwhen󰀁∈]0;1[allthepropertiesofprogram(1)holdtrue.Fornegative󰀁,ontheotherhand,thelandscapeoff󰀁(x)changesand“󰀃atregions”canmergeinanextremumwhileotherextremacandisappear,dependingonthevaluesoftheparameter󰀁.Adetailedanalysisofthesee󰀁ectswillsuggestanewalgorithmforapproximatingtheMCPwhichisbasedontheideaofvaryingtheparameter󰀁duringanevolutionaryoptimizationprocess,insuchawayastoavoidobtainingcharacteristicvectorsofsmallcliques.

WepointoutthattheproposedparameterizationoftheMotzkin–Strausprogramiscompletelydi󰀁erent,bothincontentandmotivations,fromthatrecentlyintroducedbyGibbonsetal.[21].Theirideawastosubstitutethesignconstraintsx¿0oftheMotzkin–Strausprogramwithoneoftheformx󰀃x=1=s,sbeingaparameterintheinterval[1;n],inanattempttoavoidspurioussolutions.Withthisprogramitmayhappenthatthesolutionshavetobeprojectedontothepositiveorthant,inordertomaintainfeasibility.

2.2.ReplicatorequationsandtheirapplicationtotheMCP

LetMbeanon-negativereal-valuedn×nmatrix,andconsiderthefollowingdy-namicalsystem:

x˙i(t)=xi(t)[(Mx(t))i−x(t)󰀃Mx(t)];

(Mx(t))i

;

x(t)󰀃Mx(t)i=1;:::;n;

(3)

whereadotsigniÿesderivativew.r.t.timet,anditsdiscrete-timecounterpart

xi(t+1)=xi(t)

i=1;:::;n:

(4)

ItisreadilyseenthatthesimplexSnisinvariantunderthesedynamics,whichmeansthateverytrajectorystartinginSnwillremaininSnforallfuturetimes.Moreover,itturnsoutthattheirstationarypoints,i.e.thepointssatisfyingx˙i(t)=0for(3)orxi(t+1)=xi(t)for(4),coincideandarethesolutionsoftheequations

xi[(Mx)i−x󰀃Mx]=0;

i=1;:::;n:

(5)

Astationarypointxissaidtobeasymptoticallystableifeverysolutionto(3)or(4)whichstartscloseenoughtox,willconvergetoxast→∞.

Both(3)and(4)arecalledreplicatorequationsintheoreticalbiology,sincetheyareusedtomodelevolutionovertimeofrelativefrequenciesxi(t)ofinteracting,self-replicatingentities.Eq.(3)hasbeenintroducedinevolutionarygametheoryby

32I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

TaylorandJonker[40]tomodelevolutionofbehaviorinintraspeciÿccon󰀃ictsun-derrandompairwisematinginalarge,ideallyinÿnitepopulation.Itformalizestheideathatthegrowthratesx˙i=xiofrelativefrequencyxioftheithbehaviorpattern

󰀇󰀇

(i=1;:::;n)isequaltothe(dis)advantage(Mx)i−x󰀃Mx=jmijxj−j;kmkjxjxk,measuredbyincrementalÿtnessrelativetotheaverageperformancewithinthepopula-tioninstatex=(x1;:::;xn)󰀃.Heremijdenotesincrementalindividualÿtnessattributedtoani-individualwhenencounteringaj-individual,andM=(mij)istheresultingÿt-nessmatrix.Thebehaviorpatternsi∈{1;:::;n}areoftencalled“purestrategies”andtheinteractionmatrixMisalsotermed“payo󰀁matrix”.Similarargumentsprovidearationaleforthediscrete-timeversion(4).Surprisingly,thesedynamicalequationscanalsoberegardedasaveryspecialcaseofageneralclassofdynamicalsystemsintroducedbyBaumandEagon[5]andstudiedbyBaumandSell[6]inthecontextofMarkovchaintheory.Thiskindofprocesseshaveproventobeusefulinthespeechrecognition[31]andcomputervision[38]domains.Dynamics(3)and(4)alsoariseinpopulationgeneticsunderthenameselectionequationsinamodelassumingsep-arate(non-overlapping)generations,largepopulationsize,randomunionofgametes,andaselectionactingonlyupononechromosomallocusthroughdi󰀁erentviabilities(i.e.,survivalprobabilities),givenbythetheÿtnessmatrixMofthegenotypes,i.e.,pairsofgenesdrawnfromaset{1;:::;n}ofallelesforasinglechromosomallocus.Herexiisthegenefrequencyoftheithallele.ThematrixMisinthiscontextalwayssymmetric,sincepermutedgenepairsbelongtothesamegenotype.Models(3)and(4)asselectionequationsgowaybacktoFisher[19]andKimura[29].

Fromanoptimizationpointofview,thedi󰀁erencebetweensymmetricandnon-symmetricmatricesMiscrucial.Indeed,inthesymmetriccasethequadraticformx󰀃(t)Mx(t)isincreasingalongtrajectoriesofthereplicatordynamics;thisistheFun-damentalTheoremofNaturalSelection,see,e.g.[16,26,24].

Theorem3.IfM=M󰀃thenthefunctionx(t)󰀃Mx(t)isstrictlyincreasingwithin-creasingtalonganynon-stationarytrajectoryx(t)underbothcontinuous-time(3)anddiscrete-time(4)replicatordynamics.Furthermore;anysuchtrajectoryconvergestoastationarypoint.

ApartfromthemonotonicityresultwhichprovidesaLyapunovfunctionforbothdy-namics,theprevioustheoremalsorulesoutcomplicatedattractorslikecycles,invarianttori,orevenstrangeattractors.

Toformulatetheresultswhichrelatedynamicalpropertiestooptimality,weneedsomefurthernotionsandnotations.First,considerthegeneralquadraticoptimizationproblemoverSn,

maximizex󰀃Mxsubjecttox∈Sn

andthegeneralizedLagrangian

L(x;󰀆;󰀇)=x󰀃Mx+󰀆󰀃x+󰀇(e󰀃x−1)

(6)

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4933

of(6),wherethemultipliers󰀆iand󰀇mayhavearbitrarysign.CallacriticalpointxofthegeneralizedLagrangianageneralizedKarush–Kuhn–TuckerpointifL(x;󰀆;󰀇)=x󰀃Mxirrespectiveofthesignof󰀆i.

Finally,weneedsomenotionsfromgametheory(see,e.g.,[42]):recallthatapointx∈Snissaidtobea(symmetric)Nash(equilibrium)strategyifandonlyif

y󰀃Mx6x󰀃Mx

forally∈Sn:

(7)

Furthermore,aNashstrategyxissaidtobeaneutrallystablestrategy(NSS)ifandonlyif

y󰀃Mx=x󰀃Mx

implies

x󰀃My¿y󰀃My

(8)

andanevolutionarilystablestrategy(ESS)ifandonlyiftheinequalityin(8)isstrictfory=x.

Nowwerepeatthecharacterizationresultsfrom[12]whichlinkthreedi󰀁erentÿelds:optimizationtheory,evolutionarygametheory,andqualitativetheoryofdynamicalsystems.

Theorem4.LetM=M󰀃beanarbitrarysymmetricn×nmatrixandx∈Sn.Considerthefollowingproperties:

(a1)xisanESS;i.e.;satisÿes(8)withstrictinequality;and(7);(a2)xisastrictlocalsolutionto(6);

(a3)xisanasymptoticallystablestationarypointof(3)and(4);(b1)xisaNSS;i.e.;satisÿes(7)and(8);(b2)xisalocalsolutionof(6);

(c1)xisaNashstrategy;i.e.;satisÿes(7);(c2)xisaKarush–Kuhn–Tuckerpointfor(6);

(d1)xisastationarypointunder(3)or(4);i.e.;satisÿes(5);(d2)xisageneralizedKarush–Kuhn–Tuckerpointfor(6).

Thenthefollowingimplicationsandequivalencesholdtrue:(a1)⇔(a2)⇔(a3)⇒(b1)⇔(b2)⇒(c1)⇔(c2)⇒(d1)⇔(d2).

Thepreviousresultnaturallysuggeststheuseofreplicatorequationsforapproxi-matingtheMCP.Infact,letAGbethe(symmetric)adjacencymatrixofgraphG;byputtingM=AG+12I,thereplicatordynamicalsystemwilliterativelymaximizetheobjectivefunctionof(1)andeventuallyconverge(withprobability1)toalocalmaximizer,whichbyvirtueofTheorem2,willthencorrespondtoacharacteristicvectorofamaximalcliqueofG.OnecanalsoputM=AG,inwhichcaseweob-taintheMotzkin–Strausprogram,butduetothepresenceofspuriousmaximizers,thesesolutionscanonlyprovideanapproximationofthesizeofthemaximumclique.Theempiricalresultsobtainedin[12]overnumerousDIMACSbenchmarkgraphsareencouragingandprovethee󰀁ectivenessofthisalgorithm.Theyalsoshowthattheap-proachbasedontheoriginal(non-regularized)versionoftheMotzkin–Strausproblemperformsslightlybetterthanitsregularizedcounterpart(1),intermsofcliquesize.

34I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Thismaybeintuitivelyexplainedbyobservingthat,sincealllocalmaximaarestrict,thelandscapeofthenewobjectivefunction(1)iscertainlyless󰀃atthantheoneasso-ciatedtothenon-regularizedversionandthusadynamicsthatincreasestheobjectivefunctionateverystepwillbemorepronetoendupinacloselocalmaximum.

Finally,letusnotethatrecentempiricalinvestigations[13]indicatethatthereisnosigniÿcantgaininvaryingthestartingpointofthereplicatordynamicsbyintricatepreprocessing,orusingadiscretizationof(3)di󰀁erentfrom(4).

3.Boundsfortheannealingparameter

Inthissection,weestablishboundsfortheannealingparameter󰀁relatedtothestabilityofxSunderthereplicatordynamics.Theÿrstresultsholdforgeneralsymmetricmatrices;wethenspecializetheseÿndingstothecaseofadjacencymatrices.Proposition5.Ifx∈Snisa(local)maximizerofx󰀃(A+󰀁I)xoverSnandS={i∈V:xi¿0};thennecessarily

󰀄󰀃

(Ax)i−x󰀃Ax

:i∈S:(9)󰀁¿󰀃(x)=max

x󰀃xProof.Sincexisalocalmaximizerofx󰀃(A+󰀁I)xonSn,thenTheorem4impliestheNashequilibriumcondition(7)whichforthecaseM=A+󰀁Ientails

(Ax)i+󰀁xi=(Mx)i6x󰀃Mx=x󰀃Ax+󰀁x󰀃x

foralli∈S

(notethatequalityhastoholdifi∈S,forotherwisewewouldarriveatthecontra-󰀇

dictionx󰀃Mx=i∈Sxi(Mx)i¡x󰀃Mx).Buti∈Smeansxi=0sothat󰀁¿󰀃(x)followsreadily.

Wemovetothefollowingresult:alocalmaximizerxofbothx󰀃Axandx󰀃(A+󰀁I)xoverSnnecessarilyhastobeacharacteristicvector.

Proposition6.Ifx∈Snisa(local)maximizerofbothx󰀃Axandx󰀃(A+󰀁I)xoverSnforsome󰀁=0;thennecessarilyx=xSifS={i∈V:xi¿0}.

Proof.FromTheorem2of[12]weknowthateverylocalmaximizerhastobeastationarypointundertherespectivereplicatordynamics.Henceforalli∈Swehave,dueto(5),

[Ax]i=x󰀃Ax

andalso[Ax]i+󰀁xi=x󰀃(A+󰀁I)x=x󰀃Ax+󰀁x󰀃x;

(10)

sothatallpositivecoordinatesofxhavetobeequal.Sincetheysumuptoone,theresultfollows.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4935

Nowweneedsomeadditionalnotation.First,denotebye=(1;:::;1)󰀃∈Rnanddenotethe(n−1)-dimensionalhyperplaneofallvectorsthecoordinatesofwhichsumuptozeroby

e⊥={v∈Rn:e󰀃v=0}:

Givenanarbitraryn×nmatrixM,theactionofitsquadraticformone⊥canbefullydescribedwiththehelpoftheorthoprojectorP=I−(1=n)ee󰀃ontoe⊥:indeed,foru∈e⊥wehavePu=uwhenceu󰀃Mu=u󰀃(PMP)uresults.NowPMPissymmetricifMissymmetric,andeisaneigenvectortotheeigenvaluezeroofPMP,duetoPe=0.Henceweget

u󰀃Mu¿󰀆min(M|e⊥)u󰀃u

forallu∈e⊥;

(11)

where󰀆min(M|e⊥)denotesthesmallesteigenvalueofPMP,ifthezeroeigenvalueisignoredwithmultiplicityone,i.e.,

󰀆min(M|e⊥)=min{󰀆∈R:PMPv=󰀆vforsomev∈e⊥\\{0}}:

(12)

WerecallTheorem5of[10]accordingtowhicheverylocalmaximizerzofx󰀃Axismaximizingthisfunctionoverthewholeface{y∈Sn:yi=0ifzi=0},andthisfaceiscontainedinthebasinofattractionofzunderthereplicatordynamics.Forsimplicityofexposition,weassumeinthenextresultthatthisfaceisthewholesimplexSn,inotherwords,thatzi¿0foralli.Further,inviewofProposition6wemayanddoassumethatz=b,whereb=xVisthebarycenterofSn,i.e.,zi=1=nforalli∈V.Thisgivesusanupperboundfortheparameter󰀁asfollows:

Proposition7.Ifb=xVisaglobalmaximizerofx󰀃AxonSn;thenbisalsoa(global)maximizerofx󰀃(A+󰀁I)xoverSnprovidedthat

󰀁6ÿ=󰀆min(−A|e⊥);

(13)

where󰀆min(−A|e⊥)isthesmallesteigenvaluecorrespondingtotheactionof−Aone⊥;deÿnedin(12).

Proof.TheNashequilibriumcondition(7)forbw.r.t.Aimpliesx󰀃Ab=b󰀃Abforallx∈Sn(again,otherwisewegotatleastonestrictinequality(Ab)i¡b󰀃Ab,fromwhich

󰀇

theabsurdb󰀃Ab=bi(Ab)i¡b󰀃Abwouldresult).Butthenx󰀃Ax−b󰀃Ab=x󰀃Ax−2x󰀃Ab+b󰀃Ab=(b−x)󰀃A(b−x)6−󰀆min(−A|e⊥)(b−x)󰀃(b−x)dueto(11).Ontheotherhand,wealsoget06(b−x)󰀃(b−x)=b󰀃b−2x󰀃b+x󰀃x=1=n−21=n+x󰀃x=x󰀃x−1=n=x󰀃x−b󰀃b,whence

x󰀃(A+󰀁I)x−b󰀃(A+󰀁I)b=x󰀃Ax−b󰀃Ab+󰀁(x󰀃x−b󰀃b)

6󰀁(x󰀃x−b󰀃b)−󰀆min(−A|e⊥)(b−x)󰀃(b−x)=(󰀁−ÿ)(x󰀃x−b󰀃b)60

results,provided󰀁6ÿholds,sincethelastexpressioninparenthesesisalwaysnon-negative.Hencetheresult.

(14)

36I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Summarizing,weobtainanadmissiblerangeforourparameter:

Theorem8.LetS⊆{1;:::;n};with|S|=m.LetAS=(aij)i;j∈Sbethem×msubma-⊥

trixofAcorrespondingtoS;andeS=(1;:::;1)󰀃∈Rm.DenotebyÿS=󰀆min(−AS|eS).Ifx=xS∈Snisa(local)maximizerofx󰀃AxoverSnand󰀁∈]󰀃(xS);ÿS[;thenxSisalsoastrictlocalmaximizerofx󰀃(A+󰀁I)xoverSn.

Ontheotherhand;if󰀁¡󰀃(xS);thenxSbecomesanunstablestationarypointofthereplicatordynamicsunderA+󰀁I;andthus;withprobabilityone;cannotbeapproachedbyaninteriorpathunderthesedynamics.

Proof.FromTheorem4,theclaimedassertionfollowsifwecanestablishlocalasymp-toticstabilityofxSunderthereplicatordynamicswith,say,continuoustime,andthematrixM=A+󰀁I.NowxSliesintherelativeinterioroftheface

F={x∈Sn:xi=0foralli∈S}

ofthesimplex,whichinturnisalsotimeinvariantunderdynamics(3)and(4).Asaconsequence,wecandecomposelocalstabilityanalysisintothequestionof“internal”stability(concerningconvergenceoftrajectoriesstartingnearbywithinF)and,separately,“external”stabilitydealingwithtrajectoriesstartinginSnbuto󰀁F.Externalstabilityisgovernedbythe“external”eigenvaluesofthelinearizationofthisdynamicsaroundxS,whicharegivenbythequantities(see[8,Lemma21])

[(A+󰀁I)xS]i−(xS)󰀃(A+󰀁I)(xS)=[AxS]i−(xS)󰀃A(xS)−󰀁(xS)󰀃(xS)

fori∈S:

Nowif󰀁¿󰀃(xS),thenthelatterquantityisnegativebydeÿnition(9).Ontheotherhand,internalstabilityofxSfollowsfrom󰀁6ÿSandtheoptimalityofxSw.r.t.x󰀃(A+󰀁I)xonthefaceFasinProposition7.RecallthatoptimalityonFisguaranteedbyTheorem5of[10].Hencetheresultfor󰀁¿󰀃(xS).Theinstabilityresultfollowsbythesameargumentationasin[14,Theorem6]:allstartingpointsoftrajectoriesconvergingtothenon-asymptoticallystablepointxSlieonthecenter-stablemanifold[27,28]whichalwaysisofcodimensionatleastone.Hence,atrajectorywitharandomlychosenstartingpointwillalmostsurelynotconvergetoxS.

AfurtherresultholdswhenAistheadjacencymatrixofagraphG:

Theorem9.IfA=AGistheadjacencymatrixofgraphGandSisastrictlymaximalclique;then]−1;1[⊆]󰀃(xS);ÿS[.Moreprecisely,󰀃(xS)6−1thatbecomes󰀃(xS)60ifSisjustamaximalclique;ontheothersideÿS=1foranykindofcliqueS.Proof.IfA=AGandx=xSisthecharacteristicvectorofastrictlymaximalcliqueofsizek=|S|,weknowthat󰀃(xS)cannotexceed−1,becauseof

k([AGxS]i−(xS)󰀃AG(xS))=dS(i)−(k−1)6−1

(15)

foralli∈SduetothedeÿnitionofstrictmaximalityofS(similarly,onecanshow

󰀇

that󰀃(xC)60foranymaximalcliqueC).RecallthatdS(i)=j∈Saijdenotesthe

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4937

degreeofvertexiw.r.t.S,i.e.,thenumberofverticesinSconnectedtoi.Moreover,inthesimplifyinghypothesisofProposition7,ifxV=bisaglobalmaximizerofx󰀃Ax,thenGisacompletegraph,i.e.AG=ee󰀃−I,sothat−PAGP=PandconsequentlyÿV=1.ReturningtothegeneralsituationwherewehavetoreplaceVwithS,weseethatbyanalogy,ÿS=1musthold.

Theorem10.If0¡󰀁¡1;thentheonlystrictlocalmaximizersofx󰀃(AG+󰀁I)xoverSn(i.e.theonlyattractingstationarypointsunderthereplicatordynamicswithAG+󰀁I)arecharacteristicvectorsxSwhereSisamaximalclique.Conversely;ifSisamaximalclique;thenxSrepresentsastrictlocalmaximizer.

Proof.Weshowthateveneverylocalmaximizery∈Sn(notnecessarilystrict)isacharacteristicvector,byvirtuallythesameproofasofTheorem9in[10]:tothisend,putS={i:yi¿0}.FirstweshowthatthesubgraphofGinducedbySiscomplete.Indeed,supposethatforsomei;j∈Swithi=jwehad(i;j)∈E,i.e.aij=0wouldhold.Thenforsmall󰀄¿0,thepointx=y+󰀄(ei−ej)∈Snwhereeiistheithstandardbasisvector.Straightforwardcalculationsnowyield

x󰀃(A+󰀁I)x=y󰀃(A+󰀁I)y+󰀁󰀄2¿y󰀃(A+󰀁I)y;

acontradictiontotheoptimalityofy.HencewithASandeSasdeÿnedinTheorem

󰀃

8wegetAS=eSeS−I.NowtheKarush–Kuhn–Tuckerconditionsnecessaryforlocaloptimalityyield,inparticular,ASyS+󰀁yS+󰀇eS=oforsome󰀇∈R,whichgives,using󰀃eSyS=1,

yS=

1+󰀇

eS;1−󰀁󰀃

which,againusingeSyS=1,yieldsy=xS.ItremainstoshowthatSismaximal.Sosupposethatthereisavertexi∈SsuchthatdS(i)=|S|.Butthenasin(15),

[(AG+󰀁I)xS]i−(xS)󰀃(AG+󰀁I)(xS)=

dS(i)−|S|+1−󰀁

¿0;

|S|contradictingtheNashequilibriumproperty(7)ofxSw.r.t.M=AG+󰀁I,whichisensuredbylocaloptimalityofxSduetoTheorem4.HenceSisamaximalclique.Toshowtheconverseassertion,observethatxSisalocalmaximizerofx󰀃AGxoverSnduetoTheorem1.Now󰀁∈[0;1[⊂]󰀃(xS);ÿS[(byTheorem9),andTheorem8impliesthatxSisalsoastrictlocalmaximizerofx󰀃AG+󰀁IxoverSn.

Forthecase−1¡󰀁¡0nogeneralresulthasbeenproven,butexamplescanbeprovidedinwhichnew(spurious)localmaximaemergewhicharenotcharacteristicvectorsofanysubsetofvertices,andatthesametimelocalsolutionsintheformofcharacteristicvectorsdisappear.Inthenextsection(andinanappendix)westudysmallexampleswhichillustratethispoint.

38I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Fig.1.TheMotzkin–Strausprogramf0(x)andx󰀃xwhenxvariesoverS3.

Fig.2.f1=2(x)andf−1=2(x)onS3.

4.Aprototypicalexample

Inthissectionweinvestigateasmall,butneverthelessinteresting,example,sketchedinFig.1.Itisagraphofsize3withtwomaximumcliquesofsize2intersectinginvertex3.HenceAGhaszeroentrieswiththeexceptionofa13=a31=a23=a32=1.ThisisafrequentlyconsideredcounterexampleexhibitingspurioussolutionstotheMotzkin–Strausprogram[36].Forthissimplegraphf󰀁(x)isdeÿnedonthetwo-dimensionalsimplexS3sothatwecanactuallyplotitillustratinggraphicallytheÿndingsofthepreviousparagraphs.

Moreprecisely,S3isatrianglespannedbythevertices[1;0;0]󰀃,[0;1;0]󰀃and[0;0;1]󰀃containedintheplanedescribedbytheequationx1+x2+x3=1.Ourplotstakethisplaneastheirhorizontalplanesettingtheirorigininthevertex[0;0;1]󰀃ofS3.ToremindthereaderofthissituationtheplotsofFig.1containtheS3trianglemarkedingray;thethird,vertical,axisoftheplotsreportthevaluesoff󰀁(x).

TheplotsofFig.1containthebasicMotzkin–Strausprogramandtheregularizertermx󰀃x,respectively;Fig.2containplotsoff1=2(x)andoff−1=2(x)seenfromaslightlydi󰀁erentviewpoint.Onecanintuitivelygraspwhatishappeningbyrealizing

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4939

thatthesetwoplotsareobtainedbytheÿrstoneofFig.1whentheregularizingtermis,respectively,addedorsubtracted.

LetusnowexaminetheÿguresindetailstartingfromthebasicMotzkin–StrausprograminFig.1.InthisexampletherearetwomaximumcliquesS={1;3}and

1󰀃11󰀃T

T={2;3},andtheircharacteristicvectorsxS=[12;0;2]andx=[0;2;2]give

111󰀃

f0(xS)=f0(xT)=12.Butmoregenerally,f0(x)=2forallx=[2−s;s;2]when06s612.ThisexampleshowsexplicitlythatintheMotzkin–Strausprogramtheglobalmaximizeroff0(x)arenotnecessarilycharacteristicvectorsofmaximumcliques,asstatedbythefollowinggeneralproperty(see[39]foraproofandrelatedresults):

Theorem11.LetGbeagraphcontainingmaximumcliquesC1;:::;Cq(amongpos-sibleothers).Theneveryvectorbelongingtotheconvexhulloftheircharacteristicvectorsisaglobalmaximizerofx󰀃AGxoverSnifandonlyif;foralli;j=1;:::;q;thenumberofedgeshavingoneendpointinCi\\CjandtheotherinCj\\Ciequalsmij(mij−1);wheremij=|Ci\\Cj|=|Cj\\Ci|.

Fig.2containsf1=2(x)andshowsthatinthiscasetheonlymaximacorrespondtothecharacteristicvectorsofthemaximumcliqueswhilethesecondplot,containingf−1=2(x),showsthattheisolatedmaximizerisaninteriorpointnotcorrespondingtoanycliquevector.For󰀁¿1,e.g.f3=2(x),thesituationisessentially(apartfromtheverticalscale)thatoftheregularizerterminFig.1.Theseplotsillustratealsotheroleoftheboundsof󰀁that,forthisexample,are󰀃=0andÿ=1forbothxSandxTaspredictedbyTheorem9.Theshapesoff−1=2(x),f1=2(x)andf3=2(x),representingthethreepossiblecasesof󰀁withrespecttoitsbounds,conÿrmtheresultsofTheorems8and10.

Amoredetailedanalysistogetherwiththatofanothersimpleexampleisprovidedintheappendix.

5.Theannealedreplicationheuristic

Asdiscussedpreviously,themajordrawbackofreplicatorequationsistheirinherentinabilitytoescapefromlocalmaximizersoftheobjectivefunction.Theorem8providesuswithanimmediatestrategytoavoidunwantedlocalsolutions,i.e.,maximalcliqueswhicharenotmaximum.SupposethatSisamaximalcliqueinGthatwewanttoavoid.Byletting󰀁¡󰀃(xS),itscharacteristicvectorxSbecomesanunstablestationarypointofthereplicatordynamicsunderf󰀁,andthuswillnotbeapproachedbyanyinteriortrajectory.Ofcourse,theproblemistoobtainareasonableestimatefor󰀃(xS)withoutknowingSinadvance.Furthermore,if󰀁60,itmaywellhappenthattheprocessconvergestoavectorwhichdoesnotrepresentaclique(seebelow).Sinceweareconcernedwiththemaximumcliqueproblem,

󰀃(xS)=󰀃S=maxdS(i)−|S|+1:

i∈S

(16)

40I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Asalreadynotedin(15),󰀃S6−1ifSisstrictlymaximalwhile󰀃S=0ifSisnotstrictlymaximal.Inbothcases,󰀃S¿1−|S|withequalityattainedifSisisolatedinG.Soifonewantstoavoidcliqueswithsize|S|6m,onecouldsimplyrunthealgorithmwith󰀁¡1−m61−|S|6󰀃S60,andifthereisacliqueTsuchthatstill󰀃T¡󰀁holds,thereisa(moreorlessjustiÿed)hopetoobtaininthelimitxT,whichyieldsautomaticallyalargermaximalcliqueT.Unfortunately,twoothercasescouldoccur:

(a)noothercliqueTsatisÿes󰀃T¡󰀁,i.e.,󰀁hasatoolargevalue;

(b)evenifthereissuchaclique,otherattractorscouldemergewhicharenotchar-acteristicvectorsofaclique(notethatthisisexcludedif󰀁¿0byTheorem10).

Theproperchoiceoftheparameter󰀁isthereforeatrade-o󰀁betweenthedesiretoremoveunwantedmaximalcliquesandtheemergenceofspurioussolutions.Wepresentnowthestrategyweadoptedinthischoicestressingthat,giventhelackofpreciseindications,ourprescriptionsaresupportedmainlybynumericalresultsobtainedinextensivetestsandbytheintuitionsobtainedexaminingthesetestsandsimpleexampleslikethoseofSection4andoftheappendix.

Insteadofkeepingthevalueof󰀁ÿxed,ourapproachistostartwithasu󰀄cientlylargenegative󰀁andadaptivelyincreaseitduringtheoptimizationprocess,inmuchthesamespiritasthesimulatedannealingprocedure[30].Ofcourse,inthiscasetheannealingparameterhasnointerpretationintermsofahypotheticaltemperature,andtheresultingalgorithmiscompletelydeterministic.Therationalebehindthisideaisthatforvaluesof󰀁thataresu󰀄cientlynegativeonlythecharacteristicvectorsoflargemaximalcliqueswillbestableattractivepointsforthereplicatordynamics,togetherwithasetofspurioussolutions.Asthevalueof󰀁increases,spurioussolutionsdisappearandatthesametime(characteristicvectorsof)smallermaximalcliquesbecomestable.Weexpectthatatthebeginningoftheannealingprocessthedynamicsisattractedtoward“promising”regions,andthesearchisfurtherreÿnedastheannealingparameterincreases.Insummary,theproposedalgorithmisasfollows:

1.Startwithasufficientlylargenegative󰀁.2.LetbbethebarycenterofSnandsetx=b.

3.Runthereplicatordynamicsstartingfromx,underA+󰀁Iuntilconvergenceandletxbetheconvergedpoint.

4.Unlessastoppingconditionismet,increase󰀁andgoto3.5.Select󰀁󰀇with0¡󰀁󰀇¡1(e.g.󰀁󰀇=12),runthereplicatordynamicsstartingfromcurrentxunderA+󰀁I󰀇untilconvergence,andextractamaximalcliquefromtheconvergedsolution.

Thelaststepinthealgorithmisnecessaryifwewanttoextractalsotheverticescomprisingthecliquefound,asshowninTheorem10.

Notethatwhen󰀁¡0wearenolongerguaranteedthatthetrajectoriesoftherepli-catordynamicsinstep(3)willremaininthesimplexSn,andhencex󰀃(A+󰀁I)xwillnotnecessarilyincreaseateverystep.Admittedly,inthenumericalsimulationswe

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4941

carriedoutandwhicharereportedinthefollowingsection,thisphenomenonalmostneverhappened.Inafewcasestheÿrstiterationyieldednegativeentriesintheiteratedvector,butatthefollowingstepsthevectorwasreadilyprojectedontothesimplex.Inanycase,amatrixwithnegativeelementsisnoproblem.Itissimpletoseethat,byaddingasu󰀄cientlylargeconstanttothematrixtomakeitnon-negative,thetheoryandtheoptimizationprocessareuna󰀁ected.

Itisclearthatforthealgorithmtowork,weneedtoselectanappropriate“anneal-ing”strategy.Tothisend,onecouldemploythefollowingheuristics:supposeforthemomentthattheunderlyinggraphisarandomoneinthesensethatedgesaregener-atedindependentlyofeachotherwithacertainequalprobabilityq(inapplications,qcouldbereplacedwith|E|=(n2),theactualdensity).SupposeSisanunwantedcliqueofsizem.Take󰀄¿0small,say0.01,andconsideralowerboundwhichisexceededwithprobability1−󰀄:

Theorem12.UndertherandomgraphmodelconsideracliqueSofsize|S|=m;put󰀈=1=2(n−m)anddenoteby󰀃󰀇mthefollowinglowerboundfor󰀃(xS):

󰀉

(17)󰀃󰀇m=1−(1−q)m−mq(1−q)󰀄󰀈:Then

󰀇m)6󰀄:P(󰀃(xS)6󰀃

Moreover;󰀃󰀇mexceeds1−mforallm¿mq;󰀄where

1−q√n

mq;󰀄=󰀄:q(18)

Proof.Since󰀃(xS)=maxi∈SdS(i)−m+1,andsincefordi󰀁erenti=j,thevariatesdS(i)anddS(j)arestochasticallyindependentintherandomgraphmodel,weÿrstgettheidentity

󰀇m)=[P(dS(i)6󰀃󰀇m+m−1)]n−m:P(󰀃(xS)6󰀃

Next,observethattheexpectedvalueandvarianceofdS(i)is,accordingtotheBino-mialLaw,

EdS(i)=mq

and

VardS(i)=mq(1−q):

ÄHenceCebyÄsev’sinequalitygives

P(|dS(i)−mq|¿󰀅)6mq(1−q)=󰀅2;whichentails,putting󰀅=1−(1−q)m−󰀃󰀇m¿0,

󰀇m+m−1)=P(dS(i)−mq6󰀃󰀇m+(1−q)m−1)P(dS(i)6󰀃

=P(dS(i)−mq6−󰀅)

6P(|dS(i)−mq|6󰀅)6mq(1−q)󰀅−2:

(19)

42I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Thisgivesanlowerbound󰀃󰀇mfor󰀃Swhichisexceededwithaprobabilityofleast1−󰀄asfollows:

󰀇m)=[P(dS(i)¡󰀃󰀇m+m−1)]n−mP(󰀃(xS)6󰀃

6[mq(1−q)󰀅−2]n−m=󰀄

(20)

√󰀉

provided1−(1−q)m−󰀃󰀇m=󰀅=mq(1−q)2(n−m)󰀄,whichyields(17).Since󰀈¿1=2nbydeÿnition,obviously󰀄2󰀈6󰀄1=n.Ontheotherhand,󰀃󰀇m¿1−mifandonlyifqm−󰀉

mq(1−q)󰀄󰀈¿0,whichisequivalenttom¿[(1−q)=q]󰀄2󰀈.Hencem¿mq;󰀄¿[(1−q)=q]󰀄2󰀈yields󰀃󰀇m¿1−m.Sincemq;󰀄610ifq¿0:1forall󰀄¿0,thepreviouslyobtainedhardlowerbound1−misrelaxedby󰀃󰀇minalmostallimportantapplications.Moreoverthebound󰀃󰀇mdecreaseswithincreasingmprovidedthat

󰀊󰀁󰀂

m21−q1

√1+log󰀄¿−

(n−m)2󰀄󰀈qmholds,andthisistrueformanyimportantcasesinpractice:indeed,observethatthe

latterinequalitynecessarilyholdsiftheexpressioninbracketsispositive,whichistrue,e.g.for󰀄=0:01,whenever6m6(n−m)2.Thusitmakessensetouse󰀃󰀇masaheuristicproxyforthelowerboundof󰀃(xS),toavoidbeingattractedbyacliqueofsizem.

Furthermore,awell-knownresultduetoMatula[32]accuratelypredictsthesizeofthemaximumcliqueinrandomgraphswithsu󰀄cientlymanyvertices.Let

M(n;q)=2log1=qn−2log1=qlog1=qn+2log1=q

e

+1:2(21)

Matulaprovedthat,asn→∞,thesizeofthemaximumcliqueinann-vertexq-densityrandomgraphiseither󰀎M(n;q)󰀏or󰀐M(n;q)󰀑withprobabilitytendingto1.

Thepreviousresultssuggestusasortof“two-level”annealingstrategy:thelevelofcliquesize,whichinturninducesthatofthe“actual”annealingparameter.Moreprecisely,ifwedonothaveanyaprioriinformationabouttheexpectedsizeofthemaximumclique,wecanuseMatula’sformulaM(n;q)tohaveaninitial(moreorlessaccurate)estimateofit.Letm=󰀐M(n;q)󰀑;bysettingtheinitialvaluefor󰀁(step1ofouralgorithm)atsomeintermediatevaluebetween󰀃󰀇mand󰀃󰀇m−1,e.g.󰀁=(󰀃󰀇m+󰀃󰀇m−1)=2,weexpectthatonlythecharacteristicvectorsofmaximalcliqueshavingsizemwillsurviveinf󰀁,togetherwithmanyspurioussolutions.Aftertheinitialcycle,wedecreasem,recalculate󰀃󰀇mand󰀃󰀇m−1andupdate󰀁=(󰀃󰀇m+󰀃󰀇m−1)=2instep4asinthepreviousstep.Thewholeprocessisiterateduntileithermreaches1or󰀁becomesgreaterthanzero.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4943

6.Experimentalresults

Toassessthee󰀁ectivenessoftheproposedheuristic,extensivesimulationswerecarriedoutoveraselectionofDIMACSgraphs[25],whichrepresentastandardbenchmarkforcliqueÿndingalgorithms.3Theexperimentswereconductedusingthediscrete-timeversion(4)ofthereplicatorequations.ThecodewaswrittenintheCprogramminglanguageandrunonaDigitalAlphaStationSeries200(noattemptwasmadetooptimizethecode).Foreachgraphconsidered,theproposedalgorithmwasrunbyusingthetwo-levelannealingscheduledescribedattheendoftheprevioussection.Foreachinternalcycle(step3),thereplicatoralgorithmwasiterateduntilthe(squared)distancebetweentwosuccessivestatesbecamesmallerthan10−10.Attheÿnalcycle(i.e.,step5),theparameter󰀁wassetto12,andthereplicatordynamicswasstoppedwheneitheramaximalclique(i.e.,alocalmaximizeroff1=2onSn)wasfoundorthedistancebetweentwosuccessivepointswassmallerthanaÿxedthresh-old,whichwassetton10−15(nbeingthenumberofverticesofthegraphathand).Inthelattercasetheconvergedvectorwasrandomlyperturbed,andthealgorithmrestartedfromtheperturbedpoint.Becauseoftheone-to-onecorrespondencebetweenlocalmaximizersandmaximalcliques(seeTheorem10)thissituationcorrespondstoconvergencetoasaddlepoint.

Inordertoassesstherelativemeritsoftheproposedheuristicwecomparedouralgorithmwithplainreplicatordynamicswithÿxed󰀁,i.e.,withnoannealingstrat-egy.Speciÿcally,twocaseswereconsidered:󰀁=12,whichcorrespondstotheoriginalspurious-freequadraticprogramproposedby[24]andrecentlystudiedbyBomzeetal.[10,12],and󰀁=0whichistheoriginalMotzkin–Strausformulation[33]asstudiedbyPelillo[37].Inbothcases,thereplicatorprocesswasstartedfromthebarycenterofthesimplex,anditerateduntilthesquareddistancebetweentwosuccessivestatesbecamesmallerthan10−20.Inaddition,ourresultswerecomparedwiththosere-portedbyGibbonsetal.[21]whoproposedacontinuous-basedheuristic(CBH)alsobasedonaparameterization(completelydi󰀁erentfromours)oftheMotzkin–Strausprogram.

TheresultsofourexperimentsaresummarizedinTables1and2,whichcontainarowforeachDIMACSgraphsconsidered.Thecolumnslabeledgraph,vertices,anddens.representthenameofthecorrespondinggraph,thenumberofitsverticesanditsdensity,respectively.ThecolumnMaxClique,containsthesizeofthemaximumcliquewhenknown,oralowerboundforit(thisinformationisalreadyavailableintheÿlecontainingthegraph).ThecolumnsARH,PRD(12),PRD(0)andCBHcontainthesizeofthecliquefoundusingtheproposedannealedreplicationheuristic(ARH),theplainreplicatordynamics(PRD)appliedto(2)with󰀁=12,theplainreplicatordynamics(PRD)appliedto(2)with󰀁=0—theseresultsaretakenfrom[11]—and

Wedidnotconsidergraphswheretheplainalgorithmappliedto(1)alreadyyieldsthemaximumclique,e.g.,the“c-fat”family[11].Also,afewverylargeanddensegraphswereexcludedbecauseoftheexcessivelyhighcomputationalcostrequired.

3

44I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Table1

ResultsonDIMACSbenchmarkgraphs(partI)Graphbrock200brock200brock200brock200brock400brock400brock400brock400brock800brock800brock800brock800

1234123412341212311231Vertices2002002002004004004004008008008008001000200200200200200400400400400400200200400400

Dens.0.7450.4960.6050.6580.7480.7490.7480.7490.6490.6510.6490.6500.5010.7000.7000.9000.9000.9000.5000.7000.7000.7000.9000.7000.9000.9000.700

Maxclique2112151727293133232425261530187060441340302210018¿4213¿21

CliquesizeobtainedARHPRD(1)21910131420232323181819198151245393172015125016411321

17891221201819161516158151245363272015124014371118

PRD(0)188101321222021171718178151245353372015125516401118

CBH2012141623242324201920198151246363082015145018411220

Time(s)167.7997.33124.44150.74906.26752.69554.45937.843323.313175.442697.563181.741824.6039.6640.11106.2956.0398.76156.77232.39230.72194.12425.88131.12158.41269.64838.30

san1000san2000.7san2000.7san2000.9san2000.9san2000.9san4000.5san4000.7san4000.7san4000.7san4000.9sanr200sanr200sanr400sanr400

0.70.90.50.7theGibbonsetal.CBHalgorithm[21],respectively.Finally,thecolumnlabeledtimecontaintheCPUtimerequiredbytheprocesstoprovidetheÿnalsolution.

Ascanbeseen,theresultsareveryencouraging.Infact,inalmostallcasesweobtainedlargercliqueswithARHthanPRD(12)did(theexceptionsbeingbrock4001,san2000.93andphat700-2).Inmanycases,weobtainedthesameresultsasCBHandinafewexampleswereturnedbettersolutions,e.g.,phat1500-2,san2000.92,sanr4000.5.ARHalsoperformedbetterthanPRD(0).Onlyinsixoutof46casesPRD(0)returnedalargercliquesize,thatis:brock4001,san2000.93,san4000.91,phat500-3,phat700-2,andphat1000-3.However,asdiscussedinprevioussections,duetothepresenceofspurioussolutionsintheoriginalMotzkin–Strausprogram,PRD(0)isnotabletoalwaysreturnthenodescomprisingthecliquefound:itonlyprovidesinformationaboutitssize.ItisworthnotingthattheSan-chisgraphs(the“san”family)turnedouttobeveryhardforMotzkin–Straus-basedoptimizationalgorithmssinceneitherofthethreeheuristicsfoundgoodresults.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

Table2

ResultsonDIMACSbenchmarkgraphs(partII)GraphMANNa9MANNa27ppppppppppppppp

hat300-1hat300-2hat300-3hat500-1hat500-2hat500-3hat700-1hat700-2hat700-3hat1000-1hat1000-2hat1000-3hat1500-1hat1500-2hat1500-3Vertices45378300300300500500500700700700100010001000150015001500171776

Dens.0.9270.9900.2440.4890.7440.2530.5050.7520.2490.4980.7480.2450.4900.7440.2530.5060.7540.6490.751

Maxclique1612682536936¿491144¿6210¿46¿6512¿65¿941127

CliquesizeobtainedARHPRD(1)216117825359364794159104462106491816

121176223283347743578426196289715

PRD(0)121176243383548943598446396290715

CBH1612182536935491144601046651163941021

Time(s)

45

0.8336807.81107.63301.693221.27335.83893.211729.50739.921893.962582.181965.473010.437288.224100.548598.9116251.7734.85610.16

keller4keller5

AsfarastheCPUtimeisconcerned,itisclearthatouralgorithmturnsouttobecomputationallymoreexpensivethanplainreplicatordynamicsonÿxed󰀁(see[12]forcomparison)becausethelatterissimplyasinglestepofourheuristic.Moreover,ARHisslowerthanCBH[21]whichinturnmayhaveseriousmemoryallocationproblems.However,wenotethatthecontinuous-timeversion(3)ofreplicatorequationscannaturallybemappedontohardwarecircuitry[41],therebymakingthewholealgorithmparticularlyamenabletoparallel,distributedimplementations.

Fromtheresultsobtained,itcanbeconcludedthattheproposedannealedreplicationheuristicdoesagoodjobatÿndinglargecliques,andclearlybeatstheplainreplicatordynamicsapproach,wherenoannealingstrategyisused.Moreover,itshouldbepointedoutthattheannealingscheduleadoptedisentirelybasedontheassumptionthatthegraphsathandarerandom;clearly,DIMACSgraphscanhardlysaidtobe“random,”butneverthelesstheheuristicworkedremarkablywell.Ofcourse,betterannealingstrategiescouldbedevisedifweknewsomethingabouttheunderlyingstructureofthegraphs,butinabsenceofthiskindofinformationtherandomgraphassumptionseemstobesu󰀄cientlyrobust.7.Conclusions

Wehavepresentedanewheuristicforapproximatingthemaximumcliqueproblem.Theapproachiscenteredaroundanattractivecharacterizationoftheproblemdueto

46I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

MotzkinandStraus,whichallowsustoformulateitasalinearlyconstrainedquadraticmaximizationprogram.Speciÿcally,wehaveintroducedacontrolparameter󰀁andstudiedthepropertiesoftheobjectivefunctionas󰀁varies.Wehaveshownthatwhen󰀁ispositiveallthepropertiesenjoyedbythestandardregularizationapproach[10]holdtrue;speciÿcally,inthiscaseaone-to-onecorrespondencebetweenlocal=globalmaximizersinthecontinuousspaceandlocal=globalsolutionsinthediscretespaceex-ists.Fornegative󰀁’saninterestingpictureemerges:astheabsolutevalueof󰀁growslarger,localmaximizerscorrespondingtomaximalcliquesdisappear.Wehavederivedboundsontheparameter󰀁whicha󰀁ectthestabilityofthesesolutions.Theseresultshavesuggestedtheannealedreplicationheuristic,whichconsistsofstartingfromalargenegative󰀁andthenproperlyreducingitduringtheoptimizationprocess.Foreachvalueof󰀁standardreplicatorequationsareruninordertoobtainlocalsolu-tionsofthecorrespondingobjectivefunction.Therationalebehindthisideaisthatforvaluesof󰀁withaproperlargeabsolutevalueonlylocalsolutionscorrespondingtolargemaximalcliqueswillsurvive,togetherwithvariousspuriousmaximizers.Asthevalueof󰀁isreduced,spurioussolutionsdisappearandsmallermaximalcliquesbecomestable.Anannealingscheduleisproposedwhichisbasedontheassumptionthatthegraphsbeingconsideredarerandom.ExperimentsconductedoverseveralDI-MACSbenchmarkgraphsconÿrmthee󰀁ectivenessoftheproposedapproachandtherobustnessoftheannealingstrategy.Theoverallconclusionisthattheannealingpro-ceduredoeshelptoavoidine󰀄cientlocalsolutions,byinitiallydrivingthedynamicstowardspromisingregionsinstatespace,andthenreÿningthesearchastheannealingparameterisincreased.Acknowledgements

Theauthorswouldliketothankoneoftheanonymousrefereesforprovidingcriticalcommentsandsuggestions.

AppendixFulldynamicanalysisoftwoexamples

HereweinvestigatethedynamicbehaviorofthereplicatordynamicsofAG+󰀁Ias󰀁variesovertheentirerealline,forthegraphofsize3consideredinSection4.SinceforthereplicatordynamicsonS3acompleteclassiÿcationofthereplicator󰀃owisavailable[7,9],werefrainfromrepeatingthephaseportraitsinpictureshere,butratherrefertothenumbersinthesystemusedinthecitedarticles(seeFig.6in[7]andFig.1in[9]).Apreÿxedminussign(−)symbolizestimereversaloftherespectivephaseportrait(PP).Insomecases,thePPshavetoberotatedaccordingly.Therearethreebifurcationsat󰀁=0;1;2:for󰀁¡0thesituationisessentiallythatofthesecondplotofFig.2andthereisonlytheinteriorattractor

󰀂󰀃󰀁

1−󰀁1−󰀁2−󰀁y󰀁=:;;

4−3󰀁4−3󰀁4−3󰀁I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4947

The󰀃owisdepictedasPP7inFig.6of[7].As󰀁increasesreachingtheÿrstbifur-cationat󰀁=0thesituationchangesfromthatofthesecondplotofFig.2tothatoftheMotzkin–StrausprogramofFig.1.The󰀃owstopsnotonlyattheedgeconnectingvertices1and2,butalsoalongthetrajectoriesjoiningy󰀁withthe(former)saddlepointsxSandxT,sothatwearriveatPP1.If󰀁∈]0;1[thesituationisessentiallythatoftheÿrstplotofFig.2,andwehavethepictureofPP8,whichrendersbothxSandxTas(local)attractors,andy󰀁asasaddlepointwanderingtowardsx{3}as󰀁󰀆1.Thiscompletelybreaksdownif󰀁=1where,again,the󰀃owisstoppedattwoedges:PP20emerges.If󰀁isincreasedfurther,onlytheverticesareattractingwiththeoccurrenceofy󰀁asaninteriorrepellorif󰀁exceeds2:wegetPP−35for󰀁∈]1;2]andPP−7if󰀁¿2.TheregularizingtermofFig.1substantiallydepictsthesituationafterthelastbifurcationat󰀁=2.

Inanotherexample,Gconsistsofthreeverticeswithonlytwoofthemconnected,soAGhasonlyzeroentrieswiththeexceptionofa12=a21=1.HenceS={1;2}istheuniquemaximumclique.Inthiscase,replicatordynamicsundergoesasimpleexchange-of-stabilitybifurcationas󰀁passesthrough󰀁=−1,andamoredramatic,butsimilarphenomenonoccursas󰀁=0wherethe󰀃owontwoedgesisreversedsimultaneously:indeed,for󰀁¡−1weobtain,again,PPNo.7withinteriorglobalattractoroftheform

󰀂󰀃󰀁

󰀁󰀁+1󰀁;;x󰀁=

3󰀁+13󰀁+13󰀁+1󰀃11󰀃11S

(observethatx󰀁approaches[13;3;3]as󰀁󰀇−∞whilex󰀁→x=[2;2;0]as󰀁󰀆−1).If󰀁nowincreases,thisstablestationarypointremainsatxS.Indeed,thePPisfor󰀁∈[−1;0[qualitativelythesameasPP35whilefor󰀁=0weobtainthePP−20.StillxSistheglobalattractor.If󰀁isfurtherincreasedtopositivenumbers(e.g.󰀁=12),thenlocalstabilityofxSisretained,butthereemergesasecondlocalattractorx{3}(correspondingtothemaximalclique{3})togetherwithaninteriorsaddlepoint,againatx󰀁asabove,andthePP−8results.If󰀁=1,stabilityofx{3}isretained,butthe󰀃owattheedgecontainingxSstops:PP−1depictsthesituation,whichinsomesensecorrespondstotheoccurrenceofspurioussolutionsinthepreviousexample(emergingat󰀁=0there).Finally,if󰀁¿1,wearriveagainatPP−7wherex󰀁isnowarepellorandallverticesbecomeattractors.ThisisinaccordancewiththetheoryasinTheorem8.

References

[1]S.Arora,C.Lund,R.Motwani,M.Sudan,M.Szegedy,Proofveriÿcationandthehardnessof

approximationproblems,Proceedingsofthe33rdAnnualSymposiumonFoundationsofComputerScience,Pittsburgh,PA,1992,pp.14–23.

[2]S.Arora,S.Safra,Probabilisticcheckingofproofs:anewcharacterizationofNP,Proceedingsofthe

33rdAnnualSymposiumonFoundationsofComputerScience,Pittsburgh,PA,1992,pp.2–13.

[3]M.Bellare,S.Goldwasser,C.Lund,A.Russell,E󰀄cientprobabilisticallycheckableproofsand

applicationtoapproximation,Proceedingsofthe25thAnnualACMSymposiumonTheoryComput.,SanDiego,CA,1993,pp.294–304.

48I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–49

[4]M.Bellare,S.Goldwasser,M.Sudan,Freebits,PCPsandnon-approximability—towardstightresults,

Proceedingsofthe36thAnnualSymposiumonFoundationsofComputerScience,Milwaukee,WI,1995,pp.422–431.

[5]L.E.Baum,J.A.Eagon,Aninequalitywithapplicationstostatisticalestimationforprobabilistic

functionsofMarkovprocessesandtoamodelforecology,Bull.Amer.Math.Soc.73(1967)360–363.[6]L.E.Baum,G.R.Sell,Growthtransformationsforfunctionsonmanifolds,PaciÿcJ.Math.27(2)

(1968)211–227.

[7]I.M.Bomze,Lotka–Volterraequationandreplicatordynamics:atwo-dimensionalclassiÿcation,Biol.

Cybernet.48(1983)201–211.

[8]I.M.Bomze,Non-cooperativetwo-persongamesinbiology:aclassiÿcation,Internat.J.GameTheory

15(1986)31–57.

[9]I.M.Bomze,Lotka–Volterraequationandreplicatordynamics:newissuesinclassiÿcation,Biol.

Cybernet.72(1995)447–453.

[10]I.M.Bomze,Evolutiontowardsthemaximumclique,J.GlobalOptim.10(1997)143–164.

[11]I.M.Bomze,M.Budinich,P.M.Pardalos,M.Pelillo,Themaximumcliqueproblem,in:D.Z.Du,P.M.

Pardalos(Eds.),HandbookofCombinatorialOptimization(suppl.Vol.A),Kluwer,Dordrecht,1999,pp.1–74.

[12]I.M.Bomze,M.Pelillo,R.Giacomini,Evolutionaryapproachtothemaximumcliqueproblem:

Empiricalevidenceonalargerscale,in:I.M.Bomze,T.Csendes,R.Horst,P.M.Pardalos(Eds.),DevelopmentsinGlobalOptimization,Kluwer,Dordrecht,1997,pp.95–108.

[13]I.M.Bomze,F.Rendl,Replicatordynamicsforevolutiontowardsthemaximumclique:variations

andexperiments,in:R.DeLeone,A.Murli,P.M.Pardalos,G.Toraldo(Eds.),HighPerformanceAlgorithmsandSoftwareinNonlinearOptimization,Kluwer,Dordrecht,1998,pp.53–68.

[14]I.M.Bomze,V.Stix,Geneticengineeringvianegativeÿtness:evolutionarydynamicsforglobal

optimization,Ann.Oper.Res.89(1999)297–318.

[15]M.Budinich,Boundsonthemaximumcliqueofagraph,Disc.Appl.Math.,toappear.

[16]J.F.Crow,M.Kimura,AnIntroductiontoPopulationGeneticsTheory,Harper&Row,NewYork,

1970.

[17]U.Feige,S.Goldwasser,L.LovÃasz,S.Safra,M.Szegedy,ApproximatingcliqueisalmostNP-complete,

Proceedingsofthe32ndAnnualSymposiumonFoundationsofComputerScience,SanJuan,PuertoRico,1991,pp.2–12.

[18]U.Feige,S.Goldwasser,L.LovÃasz,S.Safra,M.Szegedy,Interactiveproofsandthehardnessof

approximatingclique,J.ACM43(2)(1996)268–292.

[19]R.A.Fisher,TheGeneticalTheoryofNaturalSelection,ClarendonPress,Oxford,1930.

[20]A.H.Gee,R.W.Prager,Polyhedralcombinatoricsandneuralnetworks,NeuralComput.6(1994)161–

180.

[21]L.E.Gibbons,D.W.Hearn,P.M.Pardalos,Acontinuousbasedheuristicforthemaximumclique

problem,in:D.S.Johnson,M.Trick(Eds.),Cliques,Coloring,andSatisÿability—SecondDIMACSImplementationChallenge,AmericanMathematicalSociety,Providence,RI,1996,pp.103–124.

[22]L.E.Gibbons,D.W.Hearn,P.M.Pardalos,M.V.Ramana,Continuouscharacterizationsofthemaximum

cliqueproblem,Math.Oper.Res.22(3)(1997)754–768.

[23]J.Hastad,Cliqueishardtoapproximatewithinn1−j,Proceedingsofthe37thAnnualSymposiumon

FoundationsofComputerScience,1996,pp.627–636.

[24]J.Hofbauer,K.Sigmund,TheTheoryofEvolutionandDynamicalSystems,CambridgeUniversity

Press,Cambridge,UK,1988.

[25]D.S.Johnson,M.A.Trick(Eds.),Cliques,Coloring,andSatisÿability:SecondDIMACS

ImplementationChallenge,DIMACSSeriesinDiscreteMathematicsandTheoreticalComputerScience,Vol.26,AmericanMathematicalSociety,Providence,RI,1996(seealsohttp:==dimacs.rutgers.edu=Volumes=Vol26.html).

[26]S.Karlin,Mathematicalmodels,problemsandcontroversiesofevolutionarytheory,Bull.Amer.Math.

Soc.10(1984)221–273.

[27]A.Kelley,Thestable,center-stable,center,center-unstable,unstablemanifolds,J.Di󰀁erentialEquations

3(1967)546–570.

[28]A.Kelley,Stabilityofthecenter-stablemanifold,J.Math.Anal.Appl.18(1967)336–344.

[29]M.Kimura,Onthechangeofpopulationÿtnessbynaturalselection,Heredity12(1958)145–167.

I.M.Bomzeetal./DiscreteAppliedMathematics121(2002)27–4949

[30]S.Kirkpatrick,C.D.GelattJr.,M.P.Vecchi,Optimizationbysimulatedannealing,Science220(4598)

(1983)671–679.

[31]S.E.Levinson,L.R.Rabiner,M.M.Sondhi,Anintroductiontotheapplicationofthetheoryof

probabilisticfunctionsofaMarkovprocesstoautomaticspeechrecognition,BellSystemTech.J.62(1983)1035–1074.

[32]D.W.Matula,Thelargestcliquesizeinarandomgraph,TechnicalReportCS7608,Departmentof

ComputerScience,SouthernMethodistUniversity,1976.

[33]T.S.Motzkin,E.G.Straus,MaximaforgraphsandanewproofofatheoremofTurÃan,Canad.J.Math.

17(1965)533–540.

[34]C.H.Papadimitriou,K.Steiglitz,CombinatorialOptimization:AlgorithmsandComplexity,

Prentice-Hall,EnglewoodCli󰀁s,NJ,1982.

[35]P.M.Pardalos,Continuousapproachestodiscreteoptimizationproblems,in:G.DiPillo,F.Giannessi

(Eds.),NonlinearOptimizationandApplications,PlenumPress,NewYork,1996,pp.313–328.

[36]P.M.Pardalos,A.T.Phillips,Aglobaloptimizationapproachforsolvingthemaximumcliqueproblem,

Internat.J.ComputerMath.33(1990)209–216.

[37]M.Pelillo,Relaxationlabelingnetworksforthemaximumcliqueproblem,J.ArtiÿcialNeuralNetworks

2(1995)313–328.

[38]M.Pelillo,Thedynamicsofnonlinearrelaxationlabelingprocesses,J.Math.ImagingVision7(4)

(1997)309–323.

[39]M.Pelillo,A.Jagota,Feasibleandinfeasiblemaximainaquadraticprogramformaximumclique,J.

ArtiÿcialNeuralNetworks2(1995)411–419.

[40]P.Taylor,L.Jonker,Evolutionarilystablestrategiesandgamedynamics,Math.Biosci.40(1978)

145–156.

[41]A.Torsello,M.Pelillo,Continuous-timerelaxationlabelingprocesses,PatternRecognition33(2000)

1897–1908.

[42]J.W.Weibull,EvolutionaryGameTheory,MITPress,Cambridge,MA,1995.