ESTIMATINGCOVARIANCESTRUCTUREINHIGHDIMENSIONSByAshwiniMauryaADISSERTATIONSubmittedtoMichiganStateUniversityinpartialful˝llmentoftherequirementsforthedegreeofStatisticsDoctorofPhilosophy2016ABSTRACTESTIMATINGCOVARIANCESTRUCTUREINHIGHDIMENSIONSByAshwiniMauryaManyofscienti˝cdomainsrelyonextractingknowledgefromhigh-dimensionaldatasetstoprovideinsightsintocomplexmechanismsunderlyingthesedata.Statisticalmodelinghasbecomeubiquitousintheanalysisofhighdimensionaldataforexploringthelarge-scalegeneregulatorynetworksinhopeofdevelopingbettertreatmentsfordeadlydiseases,insearchofbetterunderstandingofcognitivesystems,andinpredictionofvolatilityinstockmarketinthehopeofavertingthepotentialrisk.Statisticalanalysisinthesehigh-dimensionaldatasetsyieldsbetterresultsonlyifanestimationprocedureexploitshiddenstructuresunderlyingthedata.Thisthesisdevelops˛exibleestimationprocedureswithprovabletheoreticalguaranteesforestimatingtheunknowncovariancestructuresunderlyingdatageneratingprocess.Ofparticularinterestareproceduresthatcanbeusedonhighdimensionaldatasetswherethenumberofsamplesnismuchsmallerthantheambientdimensionp.Duetotheimportanceofstructureestimation,themethodologyisdevelopedfortheestimationofbothcovarianceanditsinverseinparametricandaswellinnon-parametricframework.CopyrightbyASHWINIMAURYA2016Thisthesisisdedicatedtomyfamily.Fortheirendlesslove,support,andencouragement.ivACKNOWLEDGMENTSIamreallygratefultomanypeoplewhohelpedmeachievedoctorateinStatistics.ItwouldnothavebeenpossibletopursuePhDinUnitedStatesifitwasnotformyparentswhospendenormousamountoftimeande˙ortineducatingmefromtheearlieststagesofmylife.Theysupportedmeineverystepofmycareerandprovidedasafetynettofreelypursuemanydi˙erentpossibilities.AtMichiganStateUniversity,IamextremelyfortunatetobeadvisedbyProfessorHiraL.Koul,whotaughtmehowtothinkabouttheresearchproblems;Ihavebene˝tedalotfromhisclarityofthoughtandcreativeintellect.Hehasalwaysbeenconstantsourceofmotivationandencouragedmetorealizemypotential.Iamgratefultohimforprovidingthebestpossibleacademicenvironmentthatenabledmetothinkindependentlyandgrowasascienti˝cresearcher.Ialsothankhisfamilyfortheamazinghospitality,whichinmanywaysmademeathomeawayfromhomewhileatMichiganState.Ihavemuchtothankothermembersofmythesiscommitteeaswell.Dr.MarkIwen'scourseonesensingandBigandmanydiscussionsprovedveryhelpfulinmyresearchwork.IamthankfultoprofessorYuehuaCuiandDr.GraceHongforservingonmythesiscommitteeandfortakingtimeoutfromtheirbusyscheduletoteachmetheimportanceofgoodresearch.IamverygratefultoProfessorTathagataBandyopadhyayatIndianInstituteofMan-agementAhmedabad,forhissupportandencouragingmetopursueadvanceddegreefromUnitedStates.IamalsoveryfortunatetoknowProfessorArnabLahaatIndianInstituteofManagementAhmedabadandthankhimforhissupportandencouragement.AtMichiganState,IhavelearnedalotfromteachingofProfessorTapabrataMaiti,andthankhisfamilyfortheunconditionalsupport.IamthankfultoSueWatsonwhohasbeenanever-presentsourceofhelp,KimSchmuecker,andAndyHu˙ordfortheirhelpduringmanytechnicalissuesatMichiganState.vToallmyfriends,thankyouforyourunderstandingandencouragementinmymany,manymomentsofcrisis.Yourfriendshipmakesmylifeawonderfulexperience.Ican'tlistallthenameshere,butyouarealwaysonmymind.viTABLEOFCONTENTSLISTOFTABLES...................................xLISTOFFIGURES...................................xiCHAPTER1INTRODUCTION...........................11.1CovarianceStructureEstimation........................21.1.1HighDimensionalCovarianceMatrixEstimation............31.2InverseCovarianceMatrixEstimation......................31.3ThesisOverview..................................41.4Notation......................................6PartIEstimatingCovarianceStructure.................8CHAPTER2SAMPLECOVARIANCEMATRIXANDITSLIMITATIONS....92.1SampleCovarianceMatrix............................92.2WhySampleCovarianceMatrixisNOTSuitableinHighDimensions?....10CHAPTER3LOSSFUNCTIONSFORCOVARIANCEMATRIXESTIMATION.133.1LikelihoodBasedMethods............................133.2FrobeniusNormLossBasedMethods......................143.3OtherLossFunctionBasedMethods......................15CHAPTER4LEARNINGSPARSESTRUCTURE..................174.1TwoBroadClassofCovarianceMatrices....................174.2LassoTypePenalty................................194.3Discussion:.....................................21CHAPTER5ESTIMATINGAWELL-CONDITIONEDSTRUCTURE.......235.1Motivation.....................................235.2WellConditionedEstimation...........................235.3VarianceofEigenvaluesPenalty.........................24CHAPTER6LEARNINGSIMULTANEOUSSTRUCTUREWITHJOINTPENALTY266.1Motivation.....................................266.2JPENFramework.................................286.3TheoreticalAnalysisofJPENEstimators....................296.3.1ResultsonConsistency..........................316.4GeneralizedJPENEstimatorsandOptimalEstimation............336.4.1WeightedJPENEstimatorfortheCovarianceMatrixEstimation..33CHAPTER7ANALGORITHMANDITSCOMPUTATIONALCOMPLEXITY..347.1AVeryFastExactAlgorithm..........................357.1.1Derivation.................................35vii7.1.2ChoiceofRegularizationParameters..................367.1.3ChoiceofWeights.............................367.2ComputationalComplexity............................37CHAPTER8SIMULATIONS.............................398.1Preliminary....................................398.2PerformanceComparison.............................428.3RecoveryofEigen-structureandSparsity....................45PartIIEstimatingInverseCovarianceStructure..........48CHAPTER9INVERSECOVARIANCEMATRIXANDITSAPPLICATIONS...499.1Motivation.....................................499.2RelatedWork...................................499.3JointPenaltyforPrecisionMatrixEstimation.................519.4SomeApplications................................539.4.1LinearDiscriminantAnalysis......................539.4.2GaussianGraphicalModeling......................53CHAPTER10AJOINTCONVEXPENALTY(JCP)ESTIMATION.........5510.1Motivation.....................................5510.2JointConvexPenaltyEstimation........................5510.2.1ProblemFormulation...........................5610.2.2ProposedEstimator............................56CHAPTER11PROXIMALGRADIENTALGORITHMSANDITSCONVERGENCEANALYSIS...............................5811.1Introduction....................................5811.2ProximalGradientMethod............................5811.3BasicApproximationModel...........................5911.4Algorithmforoptimization............................6011.5ChoosingtheRegularizationParameter.....................6111.6ConvergenceAnalysis...............................6111.7SimulationStudy.................................6311.7.1PerformanceCriteria...........................6311.7.2StARSMethodofTuningparameterselection:.............6411.7.3SimulationResults............................6411.7.3.1ToeplitzTypePrecisionMatrix................6511.7.3.2BlockTypePrecisionMatrix.................6611.7.3.3HubGraphTypePrecisionMatrix..............6711.7.3.4NeighborhoodGraphTypePrecisionMatrix.........6811.8Summary.....................................6911.9Discussion.....................................69viiiCHAPTER12SIMULTANEOUSESTIMATIONOFSPARSEANDWELL-CONDITIONEDPRECISIONMATRIX.........................7112.1Motivation.....................................7112.2JointPenaltyEstimation:ATwoStepApproach...............7212.3WeightedJPENestimatorforprecisionmatrix.................7312.4TheoreticalAnalysisofJPENestimators....................73CHAPTER13SIMULATIONSANDANAPPLICATIONTOREALDATAANAL-YSIS...................................7513.1SimulationResults:Settings...........................7513.2PerformanceComparison.............................7613.3ColonTumorGeneExpressionDataAnalysis.................77APPENDICES......................................81APPENDIXAJPENCovarianceMatrixEstimation..............82APPENDIXBJCPforPrecisionMatrixEstimation..............88APPENDIXCJPENforPrecisionMatrixEstimation.............91BIBLIOGRAPHY....................................92ixLISTOFTABLESTable8.1Covariancematrixestimation.........................43Table8.2Covariancematrixestimation.........................44Table11.1AverageKL-Losswithstandarderrorover20replications.........65Table11.2Averagerelativeerrorwithstandarderrorover20replications.......65Table11.3AverageKL-Losswithstandarderrorover20replications.........66Table11.4Averagerelativeerrorwithstandarderrorover20replications.......66Table11.5AverageKL-Losswithstandarderrorover20replications.........67Table11.6Averagerelativeerrorwithstandarderrorover20replications.......67Table11.7AverageKL-Losswithstandarderrorover20replications.........68Table11.8Averagerelativeerrorwithstandarderrorover20replications.......68Table13.1Precisionmatrixestimation..........................76Table13.2Precisionmatrixestimation..........................77Table13.3Averagesandstandarderrorsofclassi˝cationerrorsover100replicationsin%.......................................79xLISTOFFIGURESFigure2.1Eigenvalueofsampleandpopulationcovariancematrices.........11Figure3.1Aconcavefunction..............................13Figure4.1Denseandsparsecovarianceandprecisionmatrices............18Figure6.1Comparisonofeigenvaluesofcovariancematrixestimates.........27Figure7.1TimingcomparisonofJPEN,Glasso,andPDSCE.............38Figure8.1Covariancegraphfordi˙erenttypeofmatrices...............41Figure8.2Heat-mapofzerosidenti˝edincovariancematrixoutof50realizations.Whitecoloris50/50zerosidenti˝ed,blackcoloris0/50zerosidenti˝ed.45Figure8.3Eigenvaluesplotforn=100;p=50basedon50realizationsforneigh-borhoodtypeofcovariancematrix.....................46Figure8.4Eigenvaluesplotforn=100;p=100basedon50realizationsforCov-Itypematrix..................................47Figure9.1Eigenvaluesplotofprecisionmatrix.....................52Figure9.2Illustrationofconditionalindependence..................54Figure13.1Colontumorgeneexpressiondata......................78Figure13.2Partialcorrelationnetworkofcolontumorgeneexpressiondata.....80xiCHAPTER1INTRODUCTIONTheincreasinguseoftechnologywiththedevelopmentsinstoragesystemshascreatedvastamountofhighdimensionaldataacrossmanyscienti˝cdisciplines.Examplesincludethelarge-scaleomicsdatathatenhanceourknowledgeofhumanbiology,networkdatathatexplainshowweinteractandconnectwitheachother,andthe˝nancedatathatprovidesanopportunitytobeatthemarket.Thestatisticalanalysisofthesedataischallengingduetothecurseofdimensionality.Newstatisticalmethodsareneededtomodeltheunknownstructureunderlyingthesedatasetstoleverageourscienti˝cunderstanding.Thestatisticalinferenceinhigh-dimensionaldataispossibleonlyifaninferencepro-cedureis˛exibleenoughtoexploitthehiddenstructureunderlyingthesedatasets.Thistranslatestodesigninganinferenceprocedurethatdoeswellinmodelingthestructureun-derlyingthesedatasets.Suchaninferenceprocedureoftenassumesthatmanyofhighdimensionalstructurescanberepresentedwithasmallernumberofparameterswhichisthecaseinmanyscienti˝cdisciplines.Consequently,thetheconceptofparsimonybecomescrucialinhighdimensions.Thee˙ectivenessofstatisticalestimationinhighdimensionalsettingreliesonthero-bustnessoftheprocedure,itse˚ciency,andscalability.Thelatterdependsupontheavail-abilityofscalablealgorithmictechniquesanditsabilitytoe˚cientlylearnthestructureunderlyingthesedatasets.Themaingoalofthisthesisistodevelopa˛exibleandprincipledstatisticalmethodsforuncoveringhiddenstructureunderlyingthehighdimensional,complexdatawithafocusonscienti˝cdiscovery.Inparticularthethesisaddressestwomaintasks:(i)Estimationofcovariancematricesanditsinverse,and(ii)Scalablealgorithmsforcomputingthecovariancestructurebasedontheproposedestimationprocedures,inhighdimensionalsetting.11.1CovarianceStructureEstimationInmanyscienti˝cdisciplines,astudyinvolveslargecomplexsystemswhoseoutputoftendependsonlargenumberofcomponents(variables)andtheirinteractions.Asamotivat-ingexample,insystembiology,thecellularnetworksoftenconsistsofverylargenumberofmoleculesthatinteractandexchangeinformationamongthemselves.Manyoftheexistingtechniquesdependuponthedescriptiveanalysisofmacroscopicproperties,whichincludedegreedistribution,pathlengthsandmotifpro˝leofthesemolecularnetworksordatamin-ingtoolstoidentifyclusters.Suchananalysisprovideslimitedinsightintothecomplexmechanismofthefunctionalandstructuralorganizationofbiologicalstructures.Anesti-mateoftherobustcovariancematrixcanexplainthebiologicallyimportantinteractionsandenhanceourscienti˝cunderstandingofthecomplexphenomenon.Covariancestructureestimationisoffundamentalimportanceinmultivariatedataanal-ysis.Itiswidelyusedinnumberofapplicationsincluding(i)Principalcomponentanalysis(PCA)[JohnstoneandLu[2004],Zouetal.[2006]],wherethegoalistoprojectthedataonest"k-dimensionalsubspace,andwherebestmeanstheprojecteddataexplainsasmuchofthevariationinoriginaldatawithoutincreasingk;(ii)Discriminantanalysis[Mardiaetal.[1979]]:wherethegoalistoclassifyobservationsintodi˙erentclasses,hereestimatesofcovarianceandinversecovariancematricesplayanimportantroleastheclassi˝erisoftenafunctionoftheseentities;(iii)Regressionanalysis:Ifinterestfocusesonestimationofregressioncoe˚cientswithcorrelated(orlongitudinal)data,asandwichestimatorofthecovariancematrixmaybeusedtoprovidestandarderrorsfortheestimatedcoe˚cientsthatarerobustinthesensethattheyremainconsistentundermis-speci˝cationofthecovariancestructure;and(iv)Gaussiangraphicalmodeling[YuanandLin[2007],Wainwrightetal.[2006]Wainright[2009],Yuan[2009],MeinshausenandBühlmann[2006]]:therelationshipstructureamongnodescanbeinferredfrominversecovariancematrix(alsocalledprecisionmatrix).Azeroentryintheprecisionmatriximpliesconditionalindependencebetweenthe2correspondingnodes,giventheremainingnodes.Inapplicationswheretheprobabilitydis-tributionofdataismultivariateGaussian,precisionmatrixisusedtodescribetheunderlyingdependecestructureamongthevariables.1.1.1HighDimensionalCovarianceMatrixEstimationInmanyoftheapplications,statisticianoftenencounterdatasetswheresamplesizenismuchsmallerthantheambientdimensionp,wherethelattercanbeverylargeofteninthousands,millionsorevenmore.Insuchsituations,theclassicalstatisticalestimatorstendtohavehugebiasforestimatingtheirpopulationcounterpartandmanyoftheexistingasymptotictheoriesnolongerremainvalid.Ingeneralapdimensionalcovariancematrixrequiresestimationofp(p+1)=2freeparameters.Forn<p,thisisanill-de˝nedproblem.Insuchhighdimensionalsetting,anestimationofcovariancematrixispossiblebythefactthatinmanyscienti˝cproblems,mostofthevariablesareuncorrelatedandhenceanassumptionofsparsityreducesthee˙ectivenumberofparameterstoestimate.1.2InverseCovarianceMatrixEstimationManyofscienti˝cstudiesrequireestimationofanetworkstructure,inparticulartheconditionaldependencerelationshipsamongitsnodes(variables).Inalargepopulationidentifyingthetrueinteractionamongnodesisgenerallyaveryhardproblemsincenumberofedgesscalesquadraticallytothatofnumberofnodesn.Foramotivatingexample,considertheestimationofnetworkofneurons.Withasigni˝cantimprovementintechnologyofmeasuringneuraldata,itisnowpossibletorecordthespikeactivityofhundredsofneuronsatthesametime.Acentralquestioninsuchscienti˝cstudiesishowdotheneuronscommunicateduringagiventask?Thetraditionaltimeandtrail-shiftingmethodssu˙erfromlimitationthatthesedonotaccountforbehaviorsthatencompassmultipledistinctstructuresinthebrain.Theresearch[Yatsenkoetal.[2015]]showsthattheneuralspike3patterncanbemodeledasacombinnationofsparseprecisionmatrixthataccountsforlocalinteractionandalowrankmatrixrepresentingthecommon˛uctuationsandexternalinputs.Aprecisionmatrixapproachisappealingasito˙ersfexibilityinestimatingasparseneuralnetworkandalsousefulforpredictingthefuturestatesofnetworkofneurons.AnotherimportantapplicationisinGaussiangraphicalmodeling.InGaussiangraphicalmodels,theconditionalindependencerelationshipsamongthenodes(variables)isequiva-lentlyrepresentedasamatrixofpartialcorrelationcoe˚cientsoraprecisionmatrix.ThustheestimationofnetworkisequivalenttoestimatingaprecisionmatrixwhichisobtainedbymaximizingtheGaussianlikelihoodfunctionofobservations.SincetheGaussianlikelihoodisaconcavefunctionofprecisionmatrix,itsoptimizationcanbesolvedbyanumberoffastalgorithms.Theexistingestimationmethodsofinversecovariancematricesmainlyfocusonestimat-ingtheunderlyingsparsestructureofthegivendata,anddonotaccountforminimizingtheoverdispersioninthesampleeigen-spectrum.Theproposedmethodherefortheestimationofinversecovariancematrixaddressesthisphenomenonbypenalizinganoverdispersiontermofthesampleeigenvalues.Weconsidertheestimationofprecisionmatrixinbothparamet-ricandnon-parametricframework.TheformerisbasedonGaussianlikelihoodwhereasthelatterisbasedonFrobeniusnormlossfunction.1.3ThesisOverviewThemainfocusofthisthesisistheestimationoflargedimensionalcovariancematrixanditsinversefromlimitedsampleobservations,establishthetheoreticalconsistency,andprovideaveryfastalgorithmforcomputingtheseestimates.InpartI(chapter2-chapter8),wefocusoncovariancestructureestimation.Chapter2reviewstheclassicalsamplecovariancematrixestimator,anditslimita-tionsinhighdimensionalsetting.Weaddresstheselimitationsandbuildonthesein4subsequentchapters.Chapter3reviewsthevariouslossfunctionsusedinestimationofcovariancematricesanditsinverse,therelatedoptimizationproblems,andhighlighttheiradvantagesandlimitations.Chapter4reviewsthesparsestructurelearningofcovariancematricesinhighdimen-sionalsetting,twobroadclassesofcovariancematrices,anddescribestheirestimationparadigms.Weintroduceregularizedestimationofcovariancematricesanddiscusstheirestimationframeworkwithdi˙erentlossfunctionsfromearlierchapter.Chapter5reviewstheconceptofwell-conditionedestimationanditsimportanceinhighdimensionalsettings.Someoftheeigenvaluesshrinkagemethods,theiradvantages,andlimitationsaredescribed.Weconcludethechapterbyintroducingthenewpenalty(varianceofeigenvalues)toreducetheover-dispersioninthesampleeigen-structureandthecorrespondingestimationframeworksbasedonFrobeniusnormlossfunction.InChapter6,weproposetheJointPenaltyestimationmethod.Wediscussitsasymp-toticpropertiesandratesofconvergenceinFrobeniusandspectralnorm.Wealsogivegeneralizedjointpenaltyestimatorsofcovariancematrixinhighdimensionalsettings.Chapter7containsaderivationofaveryfastalgorithmforcomputingtheproposedjointpenaltyestimatorandcomparesitscomputationaltimewithothersomeexistingmethods.InChapter8,wediscusstheextensivesimulationanalysistocomparetheperformanceoftheproposedestimatorwithsomeotherexistingmethodsforvariouschoicesofco-variancematricesinhighdimensionalsetting.Wealsoanalyzetherecoveryoftruesparseandeigen-structurebasedonjointpenaltymethods.5InpartII(chapter9-chapter13),wefocusonprecisionmatrixestimation.Chapter9reviewstheexistingliteratureandintroducestheproposedjointpenaltyframeworkforprecisionmatrixestimation.Thechapterconcludeswithdiscussionoffewapplications.InChapter10,wedescribethejointconvexpenalty(JCP)frameworkofprecisionma-trixestimationanddiscusssparseandwell-conditionedestimationunderassumptionthatthedatageneratingprocessisGaussian.Chapter11reviewstheproximalgradientalgorithm,anditsconvergenceanalysis.WegiveafastalgorithmforcomputingtheJCPestimateofprecisionmatrix.Thechapterisconcludedwithextensivesimulationanalysisforvaryingsamplesizesanddimensionsinhigh-dimensionalsetting.InChapter12,wereviewsimultaneousestimationofsparseandwell-conditionedin-versecovariancematricesinhighdimensionalsettings.WeintroduceweightedjointpenaltyestimatorsanddiscusstheirratesofconvergenceinboththeFrobeniusandspectralnorm.Chapter13reviewstheextensivesimulationanalysisusingJPENmethodforvariouschoicesofstructuredinversecovariancematrices.Weconcludethechapterwithanapplicationtocolontumorgeneexpressiondata.1.4NotationNotation:ForamatrixM,Mijdenotesits(i;j)thelement,kMk1denotesits`1normde˝nedasthesumofabsolutevaluesoftheentriesofmatrixM,kMkFdenotestheFrobeniusnormofmatrixMde˝nedassumofsquaredelementofM,kMkdenotestheoperatornorm6(alsocalledspectralnorm)de˝nedaslargestabsoluteeigenvalueofM,MdenotesmatrixMwherealldiagonalelementsaresettozero,M+denotesmatrixMwhereallo˙-diagonalelementsaresettozero,˙i(M)denotestheithlargesteigenvalueofM,tr(M)denotesitstrace,kMkdenotesitstracenormde˝nedassumofitssingularvalues,anddet(M)denotesitsdeterminant.7PartIEstimatingCovarianceStructure8CHAPTER2SAMPLECOVARIANCEMATRIXANDITSLIMITATIONS2.1SampleCovarianceMatrixGivenrandomvectors(X1;X2;;Xn)fromap-variateprobabilitydistribution,thesamplecovariancematrixisgivenby:S=[[Sij]];Sij=1n1nXk=1(XkiXi)(XkjXj);i;j=1;2;p:(2.1.1)Samplecovariancematrixisawidelyapplicableestimatorofitspopulationcounterpartandhaslowcomputationalcomplexity.Inlowdimensionalsettingwheresamplesizeissigni˝cantlylargerthanthenumberofvariables,itpossessnumberofdesiredpropertiesofagoodestimatorsuchas:Itistheoreticallyconsistent,whichmeansthatasthesamplesizedivergestoin˝nity,itconvergesalmostsurelytothepopulationcovariancematrixaslongasthedimensionpis˝xed.Itisanunbiasedestimatorofitspopulationcounterpart.Itisapproximatemaximumlikelihoodestimator.Togetherwiththevectorofsamplemeans,thesamplecovariancematrixconstitutesu˚cientstatisticsforthefamilyofGaussiandistributions.Itisinvertibleandextensivelyusedinlinearmodelsandtimeseriesanalysis.Itseigenvaluesarewellbehavedandgoodestimatorsoftheirpopulationcounterparts.Becauseofthetheseproperties,itisextensivelyusedforbothstructureestimationandpredictioninmanydataanalysisapplications.92.2WhySampleCovarianceMatrixisNOTSuitableinHighDi-mensions?Inhighdimensionalsettingwhereoftenthedimensionexceedsthesamplesize,typicallyformerisofexponentialorderoflater,manyofthesepropertiesofsamplecovariancematrixdonothold.Inhighdimensionsithasfollowinglimitations:Itisverynoisy,whichmeansthatthemanyofitsentrieshavebiases.Forn<p,itdoesnotremainspositivede˝niteandinvertible.Ithaspneigenvaluesequaltozerowhichmeansthattotalvariationindataiscontainedin˝rstneigenvaluesandthereforehighlyskewedandbiased.Infactthesampleeigenvaluesareoverdispersedinthesensethatsmallereigenvaluesarebiaseddownwardandlargereigenvaluesarebiasedupwardofthetrueeigenvalues.Figure2.1explainstheeigen-spectrumover-dispersionphenomenon.Forthisexample,wesimulatedrandomvectorsfrommultivariateGaussiandistribtuionwithmeanzeroandidentitycovariancematrix.Weconsidertwocases:case(i)Lowdimensionalsetting:n=500,p=50,andcase(ii)Highdimensionalsettingn=50,p=50.Theoverdispersioninsampleeigenvaluesarequiteapparentinthesetwosettings.Thetrueeigenvaluesareallone,whereasthesampleeigenvaluesfollowMarchenko-Pastarlaw[MarcenkoandPastur[1967]].Aresultfrom[Geman[1980]]showsthatforindependentandidentically(iid)distributedrandomvariables(thathavemeanzeroandidentityco-variancematrix)with˝nitefourthmoment,asratiopn!,thesmallestandlargestsampleeigenvaluessatisfy:10Figure2.1Eigenvalueofsampleandpopulationcovariancematricesl1!(1+p)2a:s:and(2.2.1)lp!(1p)2a:s:(2.2.2)11Inthisexample,forcase(ii),=1by(2.2.1)and(2.2.2),thesmallestandlargesteigenvalueofsamplecovariancematrixare0and4respectively.Thisshowsthatinhighdimensionthesampleeigenvaluesareoverdispersedcomparedtoitspopulationcounterpart.Becauseoftheselimitations,samplecovariancematrixisnotasuitableestimatorinhighdimensionalsettings.12CHAPTER3LOSSFUNCTIONSFORCOVARIANCEMATRIXESTIMATIONLossfunctionsarekeytoestimationproblems.Anestimatoroptimalwithrespecttoonelossfunctionmaynotbeoptimalforotherchoicesoflossfunctions.Theconsistencyandrateofconvergenceoftheestimatorsdependsuponthechoiceoflossfunction.Inthischapterwediscusssomeofthemostcommonlyusedlossfunctionsinthecontextofestimatingahighdimensionalcovariancematrix.3.1LikelihoodBasedMethodsDatalikelihood(modelbased)functionsareoneofthemostwidelyusedlossfunctionsforcovariancematrixestimation.Figure3.1AconcavefunctionThelikelihoodbasedmethodshaveadvantagethatoftentheyoutperformthenon-parametriccounterpartsinrateofconvergenceandasymptoticoptimality.Inpractice,itis13reasonabletoassumethatthedatageneratinglikelihoodfunctionissmooth.Ifthelikelihoodfunctionisstrictlyconcave(Figure3.1),theuniquemaximumlikelihoodestimatorexists.Insuchcases,maximumlikelihoodestimatorcanbeeasilycomputedusingveryfastnumericalalgorithmsuchasExpectationMaximizationalgorithm,andthealgorithmsbasedonlinearorquadraticapproximationsoflikelihoodfunctions.Anotheradvantageoflikelihoodbasedestimatoristhatthiscanbeeasilygeneralizedwhensamplescomefromamixtureofprob-abilitydistributions.Multivariatenormaldistributionisthemostwidelyusedparametricmodelforcovari-ancematrixestimation.Let(X1;X2;;Xn)followp-variatenormaldistributionwithzeromeanvectorandcovariancematrix.Thelikelihoodfunctionisgivenby:L(X1;X2;;Xn;=1(2ˇ)np=21jjp=2exp12tr(S1)g(3.1.1)Thisisconcavein1.Thereforeacommonpracticeistomaximizetheabovefunctionwithrespectto1.Let^beitssolution,thenanestimateofcovariancematrixisgivenby^1.3.2FrobeniusNormLossBasedMethodsFrobeniuslossfunctionisoneofthemostpopularalternativeestimationmethodtolikelihoodbasedmethods.Itprovidesa˛exibleestimationframeworkandhasnumberofattractivefeaturessuchas:Itisconvexandeasytosolve.Itisfullynon-parametricanddoesnotrequireanyknowledgeoffunctionalformofunderlyingdatadistributionSincetheparameterofinterestappearsinverysimpleforminthelossfunction,itiseasytointerpret.14Theconvexstructurefacilitatestheestimationprocedureanditscomputationcanbeeasilyperformedwithanumberoffastalgorithmswithlowcomputationalcomplexity.OneofmostimportantadvantageisthattheFrobeniuslossfunctionresultsinexactop-timization,unlikeinthecaseofGaussiandistribution,adirectmaximizationofthefunction(3.1)withrespecttoisverydi˚cultproblemduetoitsconvexnature.Anotherdisadvan-tageoflikelihoodbasedmodelisthatifdatadoesnotmeetthestatedmodelassumption(orifmodelisnotchoosencarefully),estimatorsbasedonlikelihoodmodelstendtoperformworsethannon-parametricestimators.Thelikelihoodfunctionmaynotalwaysbeaconcavefunction,whichmakesthecomputationverydi˚cult.Inthiscase,theestimatorsdonotremainoptimalanymore.LetSbethesamplecovariancematrix.A,un-regularizedcovariancematrixestimatorbasedonminimizationofFroebniuslossfunctionisgivenby:^=argminkSk2F(3.2.1)Theestimatoris^=S.Asdiscussedearlier,Smaynotbesuitableestimatorinhighdimensionalsetting,andanumberoftechniquesnamelyregularizationandpositivewellconditionedestimationareneededtoimprovethesamplecovarianceestimator.Formoredetailsonthistopicseechapters4and5.3.3OtherLossFunctionBasedMethodsThelikelihoodbasedmethodsrequireaprioriknowledgeabouttheunderlyingdatageneratingstochasticprocesswhichisaverystrongassumption.AlsoboththelikelihoodfunctionbasedmethodandFrobeniusnormlossfunctionbasedmethodsassumethepriorknowledgeofsamplecovariancematrixS.InpracticesamplecovariancematrixmaynotbereadilyavailablebutsomestructureofSmaybeknown.Insuchsituationsentropylossfunctionandquadraticlossfunctionaretwoothercommonlyusedlossfunctionsforcovari-15ancematrixestimation[Donohoetal.[2015]].Entropylossfunction:GivenacovariancematrixA,notnecessarilysamplecovari-ancematrix,weseekestimate^whichisobtainedmyminimizingthefollowingentropylossfunction:^=argminB˜0;B=BTtr(A1B)log(det(A1B))p(3.3.1)whereAissymmetricandpositivede˝nite.Theentropylossfunction(alsoknownbyKullback-Leiblerloss,orStein'slossfunction),isawidelyusedmethodtomeasurethediscrepancybetweentwoprobabilitydistributions.QuadraticLossFunctionAnestimatorofcovariancematrixbasedonminimizationofthequadraticlossfunctionisgivenby:^=argminB˜0;B=BT(tr(A1BI))(3.3.2)Theestimatorsbasedonentropyandquadraticlossfunctionworkwellinlowdimensionalsettingwhenn<p.Remark3.3.1.AlthoughestimatorsbasedonminimizingtheFrobenius,entropyandquadraticlossfunctionshavemanygoodproperties,toestablishtherateofconvergenceoftheseestimators,itisnecessarytoassumesomeparametricstructureontheunderlyingdatageneratingprocess.Oneofthemostcommonassumptionisthatthedatageneratingprocessissub-Gaussian.Acontinuousrandomvectorissub-GaussianifitstailsaresimilarthoseofGaussianrandomvector.Seex6:3formoredetails.16CHAPTER4LEARNINGSPARSESTRUCTUREAsdiscussedinchapter2,samplecovariancematrixisnotasuitableestimatorinhighdimensionalsettingasthisfailstoexploitthesparsestructureofthetruecovariancematrix.Inthischapter,wediscussimprovedestimationofsparsecovariancematricesbasedonregularizedlossfunctionoptimization.4.1TwoBroadClassofCovarianceMatricesInrealdataanalysisproblems,knowledgeoftruecovariancestructureisoftenunknown.Howeverinthesesituationsasuitableassumptiononthecovariancematrixstructurefacili-tatestheestimationprocedureandreducesthecomputationalcomplexity.Asamotivatingexample,giventheweathertemperatureoflocationsacrosssomegeography,weexpectthetemperatureofnearbylocationstobesimilarthanthefarawaylocations.Toeplitztypecovariancematrixwouldbesuitableformodelingthetempreturevariationsforthatgeogra-phy.Themostcommonlyusedcovariancematrixstructuresacrossmanyscienti˝cdisciplinescanbeclassi˝edintotwobroadclass:1.Naturalorderingamongvariables:Thisclassincludesthecovariancematriceswherethevariablesfarapartareweaklycorrelated.Oneoftheexampleisintimeseriesanalysis,wheretheobservationsaretypicallyautocorrelatedintime.IntheseapplicationstheresearcheroftenassumesthatthetrueunderlyingcovariancematrixhasToeplitztypeofstructure.Suchanassumptiongreatlyreducesthee˙ectivenumberofparameterstobeestimatedinthematrix.2.Nonaturalorderingamongvariables:Thisclassincludesthecovariancema-triceswherethereisnonaturalordering.Examplesincludetheanalysisofgeneexpression17datawherepriorknowledgeofanycanonicalorderingisnotavailableandsearchingoverallpermutationsofvariablesisquiteinfeasible.Inhighdimensionalsettingtypicallyn<p,andtheestimationofp(p+1)=2freeparam-etersofthecovariancematrixbasedonnobservationsisillde˝nedproblem.Theconceptofsparsestructuralassumption,whereoneoftenbelievethatonlyfewoftheentriesoftruecovariancematrixarenon-zero,greatlyreducesthee˙ectivenumberofparameterstobees-timatedandhenceimprovestheoverallestimation.Figure4.1showsthedi˙erencebetweenparsimonious(sparse)andnon-parsimonious(dense)covariancematrices.Figure4.1DenseandsparsecovarianceandprecisionmatricesThereisanextensiveliteratureontheestimationofsparsecovariancematrix[BickelandLevina[2008a],BickelandLevina[2008b],BienandTibshirani[2011],Rothman[2012],Xueetal.[2012],Dahletal.[2008].Amongearlierdevelopments,Dempster[1972]introducedtheconceptofcovarianceselectioninthecontextofprecisionmatrixestimation.Hisapproachisbasedonentrywisesparseestimationoftheprecisionmatrix.Heshowsthattheresultingestimatorcorrespondstomaximumentropymodel(maximumentropymodelismaximumsmoothmodelamongaclassofgivenmodels).Onecanfollowthesimilarprocedurefor18covariancematrixestimationinhighdimensonalsetupbysettingcertainelementsofStoequalzeroandcontinuedoingsountillthereisnosubstantialimprovementinmodel˝tting.Insuchasettingitmaynotbepossibletoderiveanexacttestofsigni˝cance,howevernumberofapproximatemethodssuchaschangein2logLikvaluecanbeusedasstoppingcriteria.Themainlimitationofsuchaprocedureisthattheresultingmatrixmaynotremainpositivede˝nite.Themethodsfortheestimationofhighdimensionalsparsecovariancematrixtendtoimposecertainstructuresassuitableonagivenclassofcovariancematrices.Fortheclassofcovariancematriceswherethevariablesareassumedtohavenaturalordering,estimatorsbasedonbandingortaperingseemtobeanaturalchoice.BickelandLevina[2008b]proposedregularizedestimationofcovariancematricesbasedonbandingwherethecorrespondingestimatorisobtainedbyselectingatmostknonzeroelementsineachrow.Achoiceofkismadebasedonre-samplingandcrossvalidation.Althoughtheirestimatorhasnaturalinterpretation,butneednotbepositivede˝nite.Toovercomethis,theyproposeataperingestimatorofcovariancematrices,whichusestheShurmatrixmultiplication.ThisisbasedonthefactthatShurmatrixmultiplicationoftwopositivede˝nitematricesisalsopositivede˝nite.FormorediscussiononthisseeBickelandLevina[2008b],Caietal.[2010],andKaroui[2008].4.2LassoTypePenaltyInsituationswherethereisnonaturalorderingamongvariables,bandingandtaperingbasedestimatorsfailtorecoverthesparsestructureofunderlyingtruecovariancematrix.Inthesesituationsthe`1regularizedcovariancematrixestimatorsaregenerallypermutationinvariant[Rothmanetal.[2008]]andbetteralternativethantheirbandingandtaperingcounterparts.The`1basedregularizedcovariancematrixestimationismotivatedbylassoinregression[Tibshiran[1996]],wherethemainideaistoshrinkthesmallerentriesof19covariancematrixtozero,whilepreservingthepositivede˝niteness.Thisprocedureisalsoknownascovariancegraphestimationofmarginallyindependentvariables.Amongthelikelihoodbasedmethods,BienandTibshirani[2011]proposedanestimatorofcovariancematrixasthesolutiontofollowingoptimizationproblem:^=argmin˜0log(det+tr(S1)+kHk1(4.2.1)whereissomepositivetuningparameter,Hisamatrixofnon-negativeweights.BothandHtogethercontrolthelevelofsparsityintheestimatedcovariancematrix.Astheaboveoptimizationproblemisconcavein,theyderiveasolutionof(4.2.1)byiterativelysolvingitsconvexapproximationusingapproach.However,suchaprocedureiscomputationalyintenseandneednotbegloballyoptimal.Amongotherrelatedworks,Chaudhurietal.[2007]considertheproblemofestimatingacovariancematrixgivenapre-speci˝edzeropattern,KhareandRajaratnam,inanunpublished2009technicalreportavailableathttp://statistics.stanford.edu/ckirby/techreports/GEN/2009/2009-01.pdf,for-mulateapriorforBayesianinferencegivenacovariancegraphstructure,andButteetal.[2000]introducetherelatednotionofrelevancenetwork,wheregeneswithpartialcorrelationexceedinggivenathresholdareconnected.AmongFrobeniusnormbasedestimationofhighdimensionalcovariancematrix,Roth-man[2012]proposedacorrelationmatrixestimatorasthesolutiontothefollowingopti-mizationproblem:^=argminTkRk2F+kk1log(det(4.2.2)whereRissamplecorrelationmatrix,issomeconstantthatensuresthepositivede˝nitenessof^.Thelog-determinantbarrierisavalidtechniquetoachievepositivede˝nitenessbutitisstillunclearwhethertheiterativeprocedureproposedinRothman[2012]actually˝ndstherightsolutiontothecorrespondingoptimizationproblem.Inanotherinterestingpaper,theauthorsinXueetal.[2012]proposedanestimatorofcovariancematrixasaminimizerof20penalizedFrobeniusnormlossfunctionoversetofpositivede˝nitematrices.Theirestimatorispositivede˝nitebutwhetheritovercomestheover-dispersionofthesampleeigenvalues,ishardtojustify.InanotherinterestinglineofworkLamandFan[2009]proposedaregularizedcovariancematrixestimatorusinganon-convexpenalty.Theyproposetheirestimatorforaclassofhard-thresholdingandSCAD(SmoothlyClippedAbsoluteDeviation)penalty.Thehard-thresholdingpenaltyisgivenby:p()=2(jj)21fg,whereastheSCADpenalyisgivenby:p()=1fg+()+1fg=(a1),forsomea>2.Themainideabehindthenon-convexpenaltyistoreducethebiaswhenthevalueofparameterhasrelativelylargermagnitude.Forexample,theSCADpenaltyremainsconstantwhenislarge,whereasthe`1penaltygrowslinearlywith.Themainlimitationsofnon-convexpenaltyisthattheproposedalgorithmsusesiterativeprocedurebasedonlocalconvexapproximationshencecomputatinalyintensive.Alsoitishardtosayiftheproposedalgorithmconvergestotheglobalminima/maxima.TheproposedJointPENalty(JPEN)methodinthisthesisusesFrobeniusnormlossfunctionandjointpenaltyof`1andvarianceofeigenvalueofunderlyingcovariancematrix.Thechoiceofsquaredlossfunctionallowssparseestimationofcovariancematrix(ratherthesparseprecisionmatrix),andresultsinveryfastalgorithm.Weintroducevarianceofeigenvaluespenaltytoensurethattheestimatedcovariancematrixispositivede˝nite.Formoredetailsonthis,seethechapter6.4.3Discussion:Assumptionofsparsityinvolvesatradeo˙betweenbene˝tandcost.Inparticularinhighdimensionaldataanalysis,whenentriesofcovariancematricesaresettozero,thenoiseduetotheerrorofestimationisgenerallyreduced.Ontheotherhand,errorsofmisspeci˝cationareintroduced.Hencethedecisionto˝tasparsemodelcomesattrade-o˙betweenover˝tting21andmodelspeci˝cation.AsnotedbyDempster[1972],oncetheparametricmodelisadopted,choiceoflevelofsparsityisoftensettleddown,speciallywhentheoptimalestimatescaneasilybecomputed.However,suchoptimalityprovidesnogauranteeagainstthecostofintroducingunecessaryparameters.22CHAPTER5ESTIMATINGAWELL-CONDITIONEDSTRUCTURE5.1MotivationInhighdimensionaldataapplications,whereaninverseofcovariancematrixisused,samplecovariancematrixcannotbeusedasgenerallythisisnotinvertible.Byawell-conditionedcovariancematrix,wemeanthatitsconditionnumber(ratioofmaximumandminimumeigenvalues)isboundedabovebyapositively˝niteconstant(heretheconstantisnottoolarge).AspointedoutbyLedoitandWolf[2004],awell-conditionedestimatorreducestheestimationerrorandisadesiredpropoertyinhighdimensionalsettings.Inthischapter,wediscusssomeoftheexistingliteratureonwell-conditionedestimation,andintroducevarianceofeigenvaluespenaltlyasane˙ectivemethodforimprovedeigen-structureestimation.5.2WellConditionedEstimationTheproblemofwell-conditionedcovariancematrixestimationisalongstudiedsubjectStein[1975,1986],LedoitandWolf[2004,2014],SheenaandGupta[2003],Wonetal.[2012].Ithasreceivedconsiderableattentioninhighdimensionalanalysisduetotheimportanceofsuchestimatorsinmanyhighdimensionaldataapplications.Amongtheearlierdevelopmentstosolvethisproblem,Stein[1975]proposedhisfamousclassofrotationinvariantshrinkageestimators.Herethemainideawastokeepthesameeigenvectorsasthatofthesamplecovariancematrixbutshrinktheeigenvaluestowardsthecenter,inordertoreducetheeigenvaluesdispersion.LetS:=UDUTbeeigen-decompositionofthesamplecovariancematrix.Stein'sestimatorisgivenby:^=UDnewUTwhereDnew=diag(dnew1;dnew2;;dnewp)(5.2.1)23withdnewii=ndiinp+1+2diipXi6=j1diidjj,where(d11;d22;;dpp)isthediagonalofD.ThisclassofestimatorsisfurtherstudiedbyHa˙[1980],LinandPerlman[1985],DeyandSrinivasan[1985],LedoitandWolf[2004,2014].AlthoughSteinestimatorisconsideredtobeStanda[Rajaratnametal.[2014]],ithasanumberoflimitationsincluding,(i)itisnotnecessarilypositivede˝nite,(ii)assumesnormality,and(iii)suitableonlyforlowdimensiondataanalysiswhensamplesizeexceedsthedimension.Amongtheearlierworkofeigenvaluesshrinkageestimationinhighdimensionalsetting,LedoitandWolf[2004]proposedanestimatorthatshrinksthesamplecovariancematrixtowardsidentity.Theirestimatorisgivenby:ˆ1S+ˆ2I;whereˆ1;ˆ2areestimatedfromdata.(5.2.2)Forˆlargeenough,theestimatorgivenby(5.2.2)iswell-conditionedbutneednotbesparseassamplecovariancematrixisgenerallynotsparse.Inanotherinterestingwork,Wonetal.[2012]considermaximumlikelhihoodestimationwithaconditionnumberconstraint.Theysolvethefollowingoptimizationproblem:MaximizeL(S;subjectto˙max=˙minmax;(5.2.3)whereL(S;:)islikelihoodfunctionofmultivariateGaussiandistributiongivenin(3.1.1),andmaxissomepositiveconstant.Theestimateof(5.2.3)isinvertibleifmaxis˝nite,andwellconditionedifmaxismoderate.Theyconsideravalueofmax<103tobemoderatebutitsvaluesalsodependsupontheeigenvaluesdispersionoftruepopulationcovariancematrix.5.3VarianceofEigenvaluesPenaltyTheestimatorsproposedbyStein[1975],DeyandSrinivasan[1985],LedoitandWolf[2004],andWonetal.[2012]arewell-conditionedandhavebeenusedinanumberofap-24plications.Howevertheseestimatorsarenotsparse,inadditiontheyhavethefollowinglimitations:Therotationinvariantestimatorsdonotchangetheeigen-vectorsandhencetheyremaininconsistent[JohnstoneandLu[2004]].Theseestimatorstendtooverestimatethenumberofnon-zerosoftheunderlyingtruecovariancematrixanditsinverse.TheestimatorgiveninLedoitandWolf[2004]resultsinlinearshrinkageofeigenvaluestowardsthoseofidentitymatrix,whichmaynotbeoptimalcriteria,aseigenvaluesfarfromcentertendtobeheavilybiasedascomparedtotheeigenvaluesinthecenter.Thechoiceofvarianceofeigenvalueshasadvantageasitallowsmoreshrinkageoftheextremeeigenvaluesthantheonesincenterandthereforenon-linearlyreducesthebias,andthequadratictermleadstoveryfastandexactalgorithm.Seechapter6formoredetailedanalysisoftheproposedmethod.25CHAPTER6LEARNINGSIMULTANEOUSSTRUCTUREWITHJOINTPENALTY6.1MotivationFromthediscussioninchapter4,itisunderstoodthatlearningasparsestructurecanbeachievedby`1regularization.Fromchapter5,itislearnedthatawell-conditionedstructurecanbeachievedbysuitablyshrinkingthesampleeigenvaluestowardsitscenter.Howevereitheroftheseregularizationdonotprovideasimultaneoustreatmentofsparseandwell-conditionedestimation.Forexample,considertheestimationofcovariancematrixbyminimizingthe`1regularizedFrobeniusnormlossfunction:^=argminT;trtr(S)jjSjj22+kk1;(6.1.1)whereissomepositiveconstant.Notethatbypenaltyfunctionkk1,weonlypenalizeo˙-diagonalelementsof.Bytheconstraint,tr=tr(S),weensurethattotalvariationintheestimatedcovariancematrixisthesameasthatinthesamplecovariancematrix.Thesolutionto(6.1.1)isthestandardsoft-thresholdingestimatoranditisgivenby(seechapter7forderivationofthisestimator):^ii=sii^ij=sign(sij)maxjsijj2;0;i6=j:(6.1.2)Itisclearfromthisexpressionthatasu˚cientlylargevalueofwillresultinsparsecovari-ancematrixestimate.Howevertheestimator^of(6.1.1)isnotnecessarilypositivede˝nite[formoredetailshereseeMaurya[2016],Xueetal.[2012]].Moreoveritishardtosaywhetheritovercomestheover-dispersioninthesampleeigenvalues.Figure6.1illustratesthisphenomenonforaneighbourhoodtypecovariancematrix.Herewesimulatedrandomvectorsfrommultivariatenormaldistributionwithsamplesizen=50anddimensionp=50.AsisevidentfromFigure6.1,eigenvaluesofsamplecovariancematrixareover-dispersedas26Figure6.1Comparisonofeigenvaluesofcovariancematrixestimatesmostofthemareeithertoolargeorclosetozero.EigenvaluesoftheproposedJointPenalty(JPEN)estimatorarewellalignedwiththoseofthetruecovariancematrix.Seechapter8fordetaileddiscussion.Thesoft-thresholdingestimator(6.1.2)issparsebutfailstorecovertheeigenstructureofthetruecovariancematrix.Toovercometheoverdispersionandachieveawell-conditionedestimator,itisnaturaltoregularizetheeigenvaluesofthesamplecovariancematrix.Considertheeigenvaluesreg-ularizedestimatorofcovariancematrixbasedonsquaredlosspenaltyasthesolutiontothe27followingoptimizationproblem:^=argminT;trtr(S)jjSjj22+pXi=1n˙i(R)˙Ro2;(6.1.3)whereissomepositiveconstant.Theminimizer^of(6.1.3)isgivenby,^=(S+tI)=(1+);(6.1.4)whereIistheidentitymatrix,andt=Ppi=1Sii=p.Toseetheadvantageofeigenvalueshrinkagepenalty,notethataftersomealgebra,forany>0,˙min(^=˙min(S+tI)=(1+)t1+>0:ThismeansthatthevarianceofeigenvaluespenaltyimprovesStoapositivede˝niteesti-mator^.Howevertheestimator(6.1.3)iswell-conditionedbutneednotbesparse.Sparsitycanbeachievedbyimposing`1penaltyontheentriesofcovariancematrix.Inthenextsectionwedescribethejointpenaltyestimationanddiscussitsadvantageasanimprovedcovariancematrixestimator.6.2JPENFrameworkWeconsiderthejointpenaltyestimatorofcovariancematrixasthesolutiontothefollowingoptimizationproblem.^=argminT;trtr(S)jjSjj22+kk1+pXi=1n˙i˙o2;(6.2.1)whereandaresomepositiveconstants.Fromhereonwardswesuppressthedependenceof^onanddenote^by^.Simulationshaveshownthat,ingeneraltheminimizerof(6.2.1)isnotpositivedef-initeforallvaluesof>0and>0.Thereforeweconsidertheoptimizationof(6.2.1)forrestrictedsetof()toensuretheresultingestimatorissparseandwell-conditioned28simultaneously.Inwhatfollows,we˝rstconsidercorrelationmatrixestimation,andlatergeneralizethemethodforcovariancematrixestimation.TheproposedJPENcovariancematrixestimatorisobtainedbyoptimizingthefollowingobjectivefunctioninRoverspeci˝cregionofvaluesof()whichdependsonthesamplecorrelationmatrixK,and.Heretheconditiontr=tr(S)reducestotr(R)=p,andthereforet=1.^RK=argminR=RT;tr(R)=pj()2^SK1jjRKjj2F+kRk1+pXi=1n˙i(R)˙Ro2;(6.2.2)where^SK1=ˆ():>0;qlogpn;8>0;˙minf(K+I)2sign(K+I)g>˙;and˙RisthemeanoftheeigenvaluesofR.InparticularifKisdiagonalmatrix,theset^SK1isgivenby,^SK1=ˆ():>0;qlogpn;8>0;<2()˙.Theminimizationin(6.2.2)overRisfor˝xed()2^SK1.FurthermoreLemmas6.2.1and6.2.2,respectivelyshowthattheobjectivefunction(6.2.2)isconvexandestimatorgivenin(6.2.3)ispositivede˝nite.Theproposedestimatorofcovariancematrix(basedonregularizedcorrelationmatrixestimator^RK)isgivenby:^K=(S+)1=2^RK(S+)1=2;(6.2.3)whereS+isthediagonalmatrixofthediagonalelementsofS.6.3TheoreticalAnalysisofJPENEstimatorsAlthoughwedonotmakeanyassumptionofdatageneratingprocessforestimation,toderiveratesofconvergencewemakeassumptionthattheunderlyingdatagenerating29processissub-Gaussian.Inthissection,wegiveratesofconvergenceoftheproposedJPENestimatorinhighdimensionalsettingwherethesamplesizeanddimensionbothdivergetoin˝nity.Def:ArandomvectorXissaidtohavesub-Gaussiandistributionifforeacht0andy2Rpwithkyk2=1,thereexist0<˝<1suchthatPfjyT(XE(X))j>tget2=2˝(6.3.1)TheJPENestimatorsexistsforany(n;p)satisfying2n<p<1.Fortheoreticalconsistencyinoperatornormwerequireslogp=o(n)andforFrobeniusnormwerequire(p+s)logp=o(n)wheresistheupperboundonthenumberofnon-zeroo˙-diagonalentriesintruecovariancematrix.Formoredetails,seetheremarkafterTheorem6.3.1.Wemakethefollowingadditionalassumptionsaboutthetruecovariancematrix0.A0.LetX:=(X1;X2;;Xp)beameanzerovectorwithcovariancematrix0suchthateachXi=p0iihassub-Gaussiandistributionwithparameter˝asde˝nedin(6.3.1).A1.WithE=f(i;j):0ij6=0;i6=jg;thejEjsforsomepositiveintegers.A2.Thereexistsa˝nitepositiverealnumberk>0suchthat1=k˙min0)˙max0)k.AssumptionA2guaranteesthatthetruecovariancematrix0iswell-conditioned(i.e.alltheeigenvaluesare˝niteandpositive).AssumptionA1ismoreofade˝nitionwhichsaysthatthenumberofnon-zeroo˙diagonalelementsareboundedbysomepositiveinteger.ThefollowingLemmas6.3.1and6.3.2,respectively,showthattheoptimizationproblemin(6.2.2)isconvexandyieldsapositivede˝nitesolution.Lemma6.3.1.Theoptimizationproblemin(6.2.2)isconvex.Lemma6.3.2.Theestimatorgivenby(6.2.2)ispositivede˝niteforany2n<1and1p<1.306.3.1ResultsonConsistencyTheorem6.3.1.Let()2^SK1and^Kbeasde˝nedin(6.2.2).UnderAssumptionsA0,A1,A2,k^RKR0kF=OPsslogpnandk^K0k=OPs(s+1)logpn;(6.3.2)whereR0istruecorrelationmatrix.Remark6.3.1.TheJPENestimator^Kismini-maxoptimalundertheoperatornorm.In(Caietal.[2015]),theauthorsobtainthemini-maxrateofconvergenceintheoperatornormoftheircovariancematrixestimatorfortheparticularconstructionofparameterspaceH0(cn;p):=ˆ:max1ipPpi=1If˙ij6=0gcn;p˙.Theyshowthatthisrateinoperatornormiscn;pqlogp=nwhichissameasthatof^Kfor1cn;p=ps.Remark6.3.2.BickelandLevina[2008b]provedthatundertheassumptionofPj=1j˙ijjqc0(p)forsome0q1,thehardthresholdingestimatorofthesamplecovariancematrixfortuningparameterq(logp)=nisconsistentinoperatornormataratenoworsethanOPc0(p)pp(logpn)(1q)=2wherec0(p)istheupperboundonthenumberofnon-zeroele-mentsineachrow.Herethetrulysparsecasecorrespondstoq=0.Therateofconvergenceof^KissameasthatofBickelandLevina[2008b]exceptinthefollowingcases:Case(i)Thecovariancematrixhasallo˙diagonalelementszeroexceptlastrowwhichhasppnon-zeroelements.Thenc0(p)=ppandps=q2pp1.TheoperatornormrateofconvergenceforJPENestimatorisOPqpp(logp)=nwhereasrateofBickelandLev-ina'sestimatorisOPqp(logp)=n.Case(ii)Whenthetruecovariancematrixistri-diagonal,wehavec0(p)=2ands=2p2,theJPENestimatorhasoperatornormrateofqplogp=nwhereasthatofBickelandLevina'sestimatorisqlogp=n.Forthecasepsc0(p)andJPENestimatorhasthesamerateofconvergenceasthatofBickelandLevina'sestimator.31Remark6.3.3.TheoperatornormrateofconvergenceismuchfasterthanFrobeniusnorm.ThisisduetothefactthatFrobeniusnormconvergenceisintermsofalleigenvaluesofthecovariancematrixwhereastheoperatornormconvergenceisintermsofthelargesteigenvalue.Remark6.3.4.Ourproposedestimatorisapplicabletoestimateanynon-negativede˝nitecovariancematrix.Notethattheestimator^Kisobtainedbyregularizationofsamplecorrelationmatrixin(6.2.2).Insomeapplicationsitisdesirabletodirectlyregularizethesamplecovariancematrix.TheJPENestimatorofthecovariancematrixbasedonregularizationofsamplecovariancematrixisobtainedbysolvingthefollowingoptimizationproblem:^S=argminT;trtr(S)j()2^SS1jjSjj2F+kk1+pXi=1f˙i˙g2;(6.3.3)where^SS1=ˆ():>0;qlogpn;8>0;˙minf(S+tI)2sign(S+tI)g>g:Theminimizationin(6.3.3)overisfor˝xed()2^SS1.Theestimator^Sispositivede˝niteandwell-conditioned.Theorem6.3.2givestherateofconvergenceoftheestimator^SinFrobeniusnorm.Theorem6.3.2.Let()2^SS1,andlet^Sbeasde˝nedin(6.3.3).UnderAssumptionsA0,A1,A2,k^S0kF=OPs(s+p)logpn(6.3.4)AsnotedinRothman[2012]theworstpartofconvergenceherecomesfromestimatingthediagonalentries.326.4GeneralizedJPENEstimatorsandOptimalEstimationTheestimatorsin(6.2.1)and(6.3.3)encourageeigenvalueshrinkagebythesameweightsforalltheeigenvalues.Howeveronemightwanttopenalizetheeigenvalueswithdi˙erentweights,especiallyifsomepriorknowledgeisavailableaboutthestructureoftrueeigenvalues.Toencouragedi˙erentlevelofshrinkagetowardsthecenter,weprovidethemoregenericestimatorsandcallitweightedJPENestimators.6.4.1WeightedJPENEstimatorfortheCovarianceMatrixEstimationAmodi˝cationofestimator^RKisobtainedbyaddingpositiveweightstotheterm(˙i(R)˙R)2.Thisleadstoweightedeigenvaluesvariancepenaltywithlargerweightsamountingtogreatershrinkagetowardsthecenterandviceversa.Notethatthechoiceoftheweightsallowsonetouseanypriorstructureoftheeigenvalues(ifknown)inestimatingthecovariancematrix.TheweightedJPENcorrelationmatrixestimator^RAisgivenby^RA=argminR=RT;tr(R)=pj()2^SK;A1jjRKjj2F+kRk1+pXi=1aif˙i(R)˙Rg2;(6.4.1)where^SK;A1=ˆ():qlogpn;(2˙min(K))(1+max(Aii)1)(1+min(Aii))1p+min(Aii)p˙;andA=diag(A11;A22;App)withAii=ai.Theproposedcovariancematrixestimatoris^K;A=(S+)1=2^RA(S+)1=2.Theoptimizationproblemin(6.4.1)isconvexandyieldsapositivede˝niteestimatorforeach()2^SK;A1.Asimpleexcerciseshowsthattheestimator^K;Ahasthesamerateofconvergenceas^S.Howtochooseweightsaiin(6.4.1),isdiscussedinnextchapter.33CHAPTER7ANALGORITHMANDITSCOMPUTATIONALCOMPLEXITYAproblemisregardedasinherentlydi˚cultifitssolutionrequiressigni˝cantresources,whateverthealgorithmused.Despiterecentambitiousdevelopmentsinsolvingconvexoptimizationproblems,e˚cientcomputationandscalabilitystillremaintwochallengingproblemsinhighdimensionsdataanalysis.Theexistingmethodsthatsolveaconvexopti-mization(herewemeanminimization)problemsoftencanbeimplementedverye˚cientlyinfarlesstimethantheconcaveoptimization.Extensiveliteratureexitsforconvexopti-mizationproblems[Bertsekas[2010]VandenbergheandBoyd[2004],Bachetal.[2011],BeckandTeboulle[2009]].ThemainchallangeincovariancematrixestimationintheGaussianliklihoodframeworkisthatthenegativeofloglikelihoodisconcavefunctionwhichmakesitaveryhardoptimizationproblem.Insuchsituations,onewaytofacilitatetheoptimiza-tionistoapproximatethenegativeofloglikelihoodfunctionbysomenon-concavefunctionandthensolvethisapproximateprobleme˚cientlyusingexistingalgorithms.Howeverthesolutionthusobtainedmaynotbeanoptimalsolutiontotheoriginalproblem.TheexistingalgorithmsofcomputingtheoptimalcovariancematrixbasedonFrobeniuslossfunctionhavecomputationalcomplexityofO(p3),wheretheconstantinO(p3)canbereallylarge,oftenmorethanthedimensionofthematrix.Themainreasonbehindsuchhighcomputationalcomplexityisthatthemethodsrequireoptimzationoverasetofpositivede˝niteconesfortheestimatortobepositivede˝nite(formoreonthistopic,seeXueetal.[2012]).TheJPENframeworkprovidesaneasysolutionforpositivede˝niteconstraintsthatdependsuponchoicesofthe().ThecompuatonalcomplexityofJPENestimatorisO(p2)andthusmuchfasterthantheotherexistingalgorithms.ThenextsectiondescribesthederivationofJPENestimatordescribedinchapter6.347.1AVeryFastExactAlgorithm7.1.1DerivationTheoptimizationproblem(6.2.2)canbewrittenas:^RK=argminR=RTj()2^SK1f(R);(7.1.1)wheref(R)=jjRKjj2F+kRk1+pXi=1f˙i(R)˙(R)g2:NotethatPpi=1f˙i(R)˙(R)g2=tr(R2)2tr(R)+p,wherewehaveusedtheconstrainttr(R)=p.Therefore,f(R)=kRKk2F+kRk1+tr(R2)2tr(R)+p=tr(R2(1+))2trfR(K+I)g+tr(KTK)+kRk1+p=(1+)ftr(R2)21+trfR(K+I)g+(1=(1+))tr(KTK)g+kRk1+p=(1+)fkR(K+I)=(1+)k2F+11+tr(KTK)g+kRk1+p:Thesolutionoftheaboveoptimizationproblemissoftthresholdingestimatorandisgivenby,^RK=11+sign(K)pmaxfabs(K+I)2;0g(7.1.2)with(^RK)ii=(Kii+)=(1+),pmax(A;b)ij:=max(Aij;b)iselementwisemaxfunctionforeachentryofthematrixA.Notethatforeach()2^SK1,^RKispositivede˝nite.357.1.2ChoiceofRegularizationParametersForagivenvalueof,wecan˝ndthevalueofsatisfying˙minf(K+I)2sign(K+I)g>0;(7.1.3)whichcanbesimpli˝edto<˙min(K+I)C12˙max(sign(K));C120:5:Then()2^SK1,andtheestimator^RKispositivede˝nite.SmallervaluesofC12yieldasolutionwhichismoresparsebutmaynotbepositivede˝nite.Theoptimalvaluesof()wereobtainedfollowingtheapproachsuggestedinBickelandLevina[2008b]byminimizingthe5foldcrossvalidationerror155Xi=1k^vivik1;where^viisJPENestimateofthecovariancematrixbasedonv=4n=5observations,viisthesamplecovariancematrixusing(n=5)observations.7.1.3ChoiceofWeightsFortheoptimizationproblemin(6.4.1),wechosetheweightsasperthefollowingcriteria:LetE=(1;2;;p)bethesetofsmallesttolargestdiagonalelementsofthesamplecovariancematrixS.LetkbethelargestintegersuchthatkthelementsofEislessthan1.Letbi=8>>><>>>:iforik1i;fork<i:A=diag(a1;a2;;ap);whereaj=bjPpi=1bi:36Suchchoiceofweightsallowsmoreshrinkageofextremesampleeigenvaluesthantheonesinthecenterofeigen-spectrum.7.2ComputationalComplexityTheJPENestimator^KhascomputationalcomplexityofO(p2).Thisisduetothefactthattherethereareatmost(p2+2p)multiplication,andatmostp2operationsforentry-wisemaximumcomputation.TheotherexistingalgorithmsGraphicalLasso(Friedmanetal.[2008]),andPDSCE(Rothman[2012])havecomputationalcomplexityofO(p3),wheretheconstantofcomplexityisoftenverylarge,mainlyduetotheiterativenatureofconvergence.AnotheradvantageofJPENestimatoristhatitisanexactsolutiontotheunderlyingoptimizationproblem.Toseethecomputingtimeperformance,weplotthecomputationaltimingofouralgorithmandsomeotherexistingalgorithmsincludingGlasso(Friedmanetal.[2008]),PDSCE(Rothman[2012]).Notethattheexacttimingofthesealgorithmalsodependsupontheimplementation,platformetc.(wedidourcomputationsinRonaAMD2.8GHzprocessor).37Figure7.1TimingcomparisonofJPEN,Glasso,andPDSCE.Figure7.1illustratesthetotalcomputationaltimetakentoestimatethecovariancematrixbyGlasso;PDSCEandJPENalgorithmsfordi˙erentvaluesofpforToeplitztypeofcovariancematrix(seechapter8forToeplitztypeofcovariancematrix).Althoughtheproposedmethodrequiresoptimizationoveragridofvaluesof()2^SK1,ouralgorithmisveryfastandeasilyscalabletolargescaledataanalysisproblems.38CHAPTER8SIMULATIONSInthischapterwecomparetheperformanceoftheproposedJPENestimatorofcovari-ancematrixforvariouschoicesofstructuredcovariancematrices.Weconsidercovariancematricesfrombothclassviz.(i)whenthereisanaturalorderingamongvariables,and(ii)whenthereisnonaturalorderingamongvariables.Inadditiontothese,wealsoincluderesultsinasettingwhentheunderlyingtruecovariancematrixisdenseandhasveryhighconditionnumber.8.1PreliminaryWeconsiderthefollowing˝vedi˙erenttypesofcovariancematricesinoursimulations.(i)HubGraph:Heretherows/columnsof0arepartitionedintoJequally-sizeddisjointgroups:fV1[V2[;:::;[VJg=f1;2;:::;pg;eachgroupisassociatedwithapivotalrowk.LetsizejV1j=s.Weset˙0i;j=˙0j;i=ˆfori2Vkand˙0i;j=˙0j;i=0otherwise.Inourexperiment,J=[p=s];k=1;s+1;2s+1;:::;andwealwaystakeˆ=1=(s+1)withJ=20.(ii)NeighborhoodGraph:We˝rstuniformlysample(y1;y2;:::;yn)fromaunitsquare.Wethenset˙0i;j=˙0j;i=ˆwithprobability(p2ˇ)1exp(4kyiyjk2).Theremainingentriesof0aresettobezero.Thenumberofnonzeroo˙-diagonalelementsofeachroworcolumnisrestrictedtobesmallerthan[1=ˆ],whereˆissettobe0.245.(iii)ToeplitzMatrix:Weset˙0i;j=2fori=j;˙0i;j=j0:75jjijj,forjijj=1;2;and˙0i;j=0,otherwise.(iv)BlockToeplitzMatrix:Inthissetting0isablockdiagonalmatrixwithvaryingblocksize.Forp=500,numberofblocksis4andforp=1000,thenumberofblocksis6.EachblockofcovariancematrixistakentobeToeplitztypematrixasinthecase(iii).(v)Cov-ItypeMatrix:Inthissetting,we˝rstsimulatearandomsample(y1;y2;:::;yp)39fromstandardnormaldistribution.Letxi=jyij3=2(1+1=p1+log(1+1=p2)).Nextwegen-eratemultivariatenormalrandomvectorsZ=(z1;z2;:::;z5p)withmeanvectorzeroandidentitycovariancematrix.LetUbeeigenvectorcorrespondingtothesamplecovariancematrixofZ.Wetake0=UDU0,whereD=diag(x1;x2;::::xp).Thisisnotasparsesettingbutthecovariancematrixhasmostofeigenvaluesclosetozeroandhenceallowsustocomparetheperformanceofvariousmethodsinasettingwheremostofeigenvaluesareclosetozeroandwidelyspreadascomparedtostructuredcovariancematricesinthecases(i)-(iv).Figure8.1showsthegraphicalcovariancestructureforthesematrices.Herewechoosep=100forbettervisualization.40Figure8.1Covariancegraphfordi˙erenttypeofmatrices41Forallthesechoicesofcovarianceandinversecovariancematrices,wegeneraterandomvectorsfrommultivariatenormaldistributionswithvaryingnandp.Wechosen=50;100andp=500;1000.Wecomparetheperformanceoftheproposedcovariancematrixestimator^Kwiththefollowingestimators.Graphicallasso[Friedmanetal.[2008]]:Graphicallassoestimatesasparseprecisionmatrix.Hereweinverttheinverse,andincludeinouranalysis.Theestimatewascom-putedusing`R'package'Glasso'.Formoredetails,refertohttp://statweb.stanford.edu/tibs/glasso/.BickelandLevina'sthresholdingestimator(BLThresh)[BickelandLevina[2008b]].Theestimatorwascomputedasperthealgorithmgivenintheirpaper.Rothman'sPositiveDe˝niteSparseCovarianceMatrixEstimator(PDSCE)[Rothman[2012]].ThePDSCEwascomputedusing`R'package'PDSCE'.Formoredetails,referto(http://cran.r-project.org/web/packages/PDSCE/index.html)LedoitandWolfestimator[LedoitandWolf[2004]]Theirestimatewascomputedusingcodefrom(http://econ.uzh.ch/faculty/wolf/publications.html#9).TheJPENestimatorwascomputedusing`R'package'JPEN'.AllthecomputationsweredoneusingRonaAMD2.8GHzprocessor.8.2PerformanceComparisonForeachofcovarianceandprecisionmatrixestimate,wecalculateAverageRelativeError(ARE)basedon50iterationsusingfollowingformula,ARE;^=jlog(f(S;^log(f(S;0))j=jlog(f(S;0))j;(8.2.1)wheref(S;)ismultivariatenormaldensitygiventhesamplecovariancematrixS,0isthetruecovariance,^istheestimateof0basedononeofthemethodsunderconsideration.42OtherchoicesofperformancecriteriaareKullback-Leibler,usedbyYuanandLin[2007]andBickelandLevina[2008b],precisionandrecall.Theoptimalvaluesoftuningparameterswereobtainedoveragridofvaluesbyminimizing5foldcross-validationasexplainedinx7:2.Table8.1CovariancematrixestimationBlocktypecovariancematrixn=50n=100p=500p=1000p=500p=1000Ledoit-Wolf1.54(0.102)2.96(0.0903)4.271(0.0394)2.18(0.11)Glasso0.322(0.0235)3.618(0.073)0.227(0.098)2.601(0.028)PDSCE3.622(0.231)4.968(0.017)1.806(0.21)2.15(0.01)BLThresh2.747(0.093)3.131(0.122)0.887(0.04)0.95(0.03)JPEN2.378(0.138)3.203(0.144)1.124(0.088)2.879(0.011)Hubtypecovariancematrixn=50n=100p=500p=1000p=500p=1000Ledoit-Wolf2.13(0.103)2.43(0.043)1.07(0.165)3.47(0.0477)Glasso0.511(0.047)0.551(0.005)0.325(0.053)0.419(0.003)PDSCE0.735(0.106)0.686(0.006)0.36(0.035)0.448(0.002)BLThresh1.782(0.047)2.389(0.036)0.875(0.102)1.82(0.027)JPEN0.732(0.111)0.688(0.006)0.356(0.058)0.38(0.007)Neighborhoodtypecovariancematrixn=50n=100p=500p=1000p=500p=1000Ledoit-Wolf1.36(0.054)2.89(0.028)1.1(0.0331)2.32(0.0262)Glasso0.608(0.054)0.63(0.005)0.428(0.047)0.419(0.038)PDSCE0.373(0.085)0.468(0.007)0.11(0.056)0.175(0.005)BLThresh1.526(0.074)2.902(0.033)0.870(0.028)1.7(0.026)JPEN0.454(0.0423)0.501(0.018)0.086(0.045)0.169(0.003)43Table8.2CovariancematrixestimationToeplitztypecovariancematrixn=50n=100p=500p=1000p=500p=1000Ledoit-Wolf1.526(0.074)2.902(0.033)1.967(0.041)2.344(0.028)Glasso2.351(0.156)3.58(0.079)1.78(0.087)2.626(0.019)PDSCE3.108(0.449)5.027(0.016)0.795(0.076)2.019(0.01)BLThresh0.858(0.040)1.206(0.059)0.703(0.039)1.293(0.018)JPEN2.517(0.214)3.205(0.16)1.182(0.084)2.919(0.011)Cov-Itypecovariancematrixn=50n=100p=500p=1000p=500p=1000Ledoit-Wolf33.2(0.04)36.7(0.03)36.2(0.03)48.0(0.03)Glasso15.4(0.25)16.1(0.4)14.0(0.03)14.9(0.02)PDSCE16.5(0.05)16.33(0.04)16.9(0.03)17.5(0.02)BLThresh15.7(0.04)17.1(0.03)13.4(0.02)17.5(0.02)JPEN7.1(0.042)11.5(0.07)8.4(0.042)7.8(0.034)Theaveragerelativeerrorandtheirstandarddeviations(inpercentage)forcovariancematrixestimatesaregiveninTable8.1andTable8.2.Thenumbersinthebracketarethestandarderrorsofrelativeerrorbasedontheestimatesusingdi˙erentmethods.AmongallthemethodsJPENandPDSCEperformsimilarformostofchoicesofnandpforall˝vetypeofcovariancematrices.ThisisduetothefactthatbothPDSCEandJPENusequadraticoptimizationfunctionwithadi˙erentpenaltyfunction.ThebehaviorofBickelandlevina'sestimatorisquitegoodinToepltizcasewhereitperformsbetterthantheothermethods.Forthistypeofcovariancematrix,theentriesawayfromthediagonaldecaytozeroandthereforesoft-thresholdingestimatorslikeBLThreshperformbetterinthissetting.Howeverforneighorhoodandhubtypecovariancematrixwhicharenotnecessarilybandedtype,BickelandLevinaestimatorisnotanaturalchoiseastheirestimatorwouldfailtorecovertheunderlyingsparsitypattern.TheperformanceofLedoit-WolfestimatorisnotveryencouragingforCov-Itypematrix,ThisisbeacuseLedoit-Wolfestimatorisgenerallynotsparseanduniformlyshrinksthesamplecovariancematrixtowardsidentitymatrix.44Figure8.2Heat-mapofzerosidenti˝edincovariancematrixoutof50realizations.Whitecoloris50/50zerosidenti˝ed,blackcoloris0/50zerosidenti˝ed.8.3RecoveryofEigen-structureandSparsityToseetheperformanceJPENestimatorinrecoveringtheeigenstructureandsparsity,weplottherecoveredeigenstructureandsparsitypatterninthebothsettingsnamely:(i)whenthetruecovariancematrixhaslowconditionnumber,and(ii)whenthetruecovariancematrixhasaveryhighconditionnumber.Theeigen-plotsinFigure8.3and8.4showthatamongallthemethods,estimatesofeigenvaluesofJPENestimatoraremostconsistentforthetrueeigenvalues.ForCov-Itypecovariancematrixwheremostofeigenvaluesareclosetozeroandwidelyspread,theperformanceofJPENestimatorisimpressive.ThisclearlyshowstheadvantageofJPENestimatorofcovariancematrixwhenthetrueeigenvaluesare45dispersedorclosetozero.TheeigenvaluesplotinFigure8.4showsthatwheneigen-spectrumofthetruecovariancematrixarenothighlydispersed,theJPENandPDSCEestimatesofeigenvaluesarealmostthesame,becausebothofthesemethodsexploitFrobeniusnormlossfunctionwithdi˙erentpenaltyfunctions.Figure8.3Eigenvaluesplotforn=100;p=50basedon50realizationsforneighborhoodtypeofcovariancematrix46Figure8.4Eigenvaluesplotforn=100;p=100basedon50realizationsforCov-Itypematrix47PartIIEstimatingInverseCovarianceStructure48CHAPTER9INVERSECOVARIANCEMATRIXANDITSAPPLICATIONS9.1MotivationInmanyofthescienti˝capplicationstheinterestistoestimatetheinverseoftheco-variancematrix(alsocalledprecisionmatrix).Precisionmatrixiswidelyusedinavarietyofapplicationsincluding:i)lineardiscriminantclassi˝cation:theclassi˝erisfunctionofpre-cisionmatrixii)Gaussiangraphicalmodeling:Azeroentryoftheprecisionmatriximpliesconditionalindependencebetweenthecorrespondingvariables,iii)regressionanalysis:theregressioncoe˚cientsarefunctionsofprecisionmatrix,andiv)Con˝denceintervalestima-tionofpopulationmeanvectorwhentheunderlyingdatadistributionisGaussian.Inhighdimensionaldataanalysisproblemswhereoneoftenhasveryfewobservationsascomparedtothenumberofvariables,theinverseofsamplecovariancematrixisnotveryuseful.Infactforn<pthesamplecovarianceissingularandtheinverseisnotde˝ned.Insuchsituationsonecanreplacetheprecisionmatrixbyitsgeneralizedinverse[Rao[1972]].AmatrixGiscalledgeneralizedinverseofA,ifAGA=A.Generalizedinversehavebeenusedinnumberofapplicationsincludingsystemofnormalequationswithsingularmatrix,leastsquarestheorytoexpressthevarianceofestimates.Althoughtheseareappealinginnumberofapplications,theirapplicationinhighdimensionalanalysisisoftenlimitedduetofactthattheystillremainsingularandnon-sparse.Inparticulartheireigenvaluesarebiased,stillremainover-dispersedcomparedtotheirpopulationcounterparts.9.2RelatedWorkRegularizationisthemostwidelyusedtechniquetoimposesomestructureontheestimationofprecisionmatrices.Inthelikelihoodbasedestimationframework,thereis49plentyofliteraturefortheestimationofsparseprecisionmatrix.Dempster[1972]introducedtheconceptofcovarianceselectionwherecertainentriesofprecisionmatricesaresettozero,locationofsuchentriesinprecisionmatrixisbasedoninformationinsamplecovarinacematrix.Inlikelihoodframework,Gaussiandistributionismostwidelyuseddatadistributionforcovariancematrixestimationduetoitsconcavityasafunctionofprecisionmatrix.Insuchframework,theoptimizationcanbecarriedoutbyanumberoffastalgorithms,viz.,interiorbasedmethods(VandenbergheandBoyd[2004],Vandenbergheetal.[1998]),subgradientbasedmethods(BeckandTeboulle[2009]),andproximalgradientmethods(Bertsekas[2010]).Amongtheearlyworkon`1regularizedhighdimensionalcovariancematrixestimation,Banerjeeetal.[2008]proposedanestimatorofsparseprecisionnatrix(basedon)asthesolutiontothefollowingoptimizationproblem:^1=argmaxlog(det+12tr(S+kk1;(9.2.1)whereisnon-negativetuningparameterthatcontrolsthelevelofsparsityintheestimate.Theyshowthattheaboveproblemisconvexandconsidertheestimationof(ratherthan1)andarguethatonecansolvetheproblembyoptimizingovereachrowandcorrespondingcolumnof=1inablockcoordinatedescentfashion.Inanotherinterestingwork,MeinshausenandBühlmann[2006]takeanotherapproachbysolvingalassoproblemtoeachvariable,usingothersaspredictors.Thecomponent^1ijisestimatedtobezeroifeithertheestimatedcoe˚cientofvariableionj,ortheestimatedcoe˚cientofvariablejoni,isnon-zero(alternativelytheyuseanANDrule).Theyshowthatasymptotically,thisconsistentlyestimatesthesetofnon-zeroelementsof1.HoweverasarguedinFriedmanetal.[2008],thisapproachdoesnotyieldmaximumlikelihoodestimator.Friedmanetal.[2008]proposedGraphicallassoestimatorassolutionto(9.1)butinsteadtheysolvefor.Usingthefactthat(9.1)isasmoothfunctionofexceptattheorigin,theyusesubgradientbasedapproachtosolvealassoleastsquareregressionofeachvariable,takingothersaspredictors.Thegraphicallassousescordinatedescentalgorithmandisveryfast.Rothman50etal.[2008]proposedanestimator(theycallitSPICE)ofprecisionmatrixassolutionto`1regularizedloglikelihoodfunction,howevertheyonlypenalizeo˙-diagonalentries.Otherauthorshaveproposedexactminimizationof`1regularizedlog-likelihood;YuanandLin[2007],Yuan[2009],Zhouetal.[2011],andPourahmadi[2007,2011].InCaietal.[2011],theauthorsproposedanestimatorofprecsionmatrixasasolutiontothefollowingoptimizationproblem.minkk1subjecttokSIk12Rpp(9.2.2)whereisatuningparameter.Theygaveasymmetricestimatorbasedonthissolutionbytaking^!ij=^!ij1fj^!ij^!jijg+^!ji1fj^!ji^!ijjgwhere^!ijisthe(i;j)thentryofthesolution^of(9.2.2).Next,wedescribethejointpenaltyestimatorofprecisionmatrix.9.3JointPenaltyforPrecisionMatrixEstimationThe`1regularizedprecisionmatrixestimatorsoftenperformverygoodinestimatingasparseprecisionmatrix,howeverwhethertheyovercometheover-dispersionineigenvalues,itishardtojustify.Figure9.3.1showstheeigenvaluesoftheestimatedinverseprecisionmatrixbasedondi˙erentmethods.Themethodsconsiderhereare:(i)JointPenalty(JPEN)(seechapter12forJPENestimator),(ii)GraphicalLasso(Glasso),(iii)SPICE,and(iv)CLIME.51Figure9.1EigenvaluesplotofprecisionmatrixAmongthesemethods,theproposedJointPenaltyestimatedeigenvaluesareclosesttothetrueeigenvalues,asitclearlyreducestheoverdispersion.Thissuggeststhatincludingapenaltyonthevarianceofeigenvaluesimprovestheoverdispersionintheeigenvaluesofprecisionmatrix.Nextwehighlightfewapplicationsofprecisionmatrixcontextofhighdimensionaldataanalysis.529.4SomeApplications9.4.1LinearDiscriminantAnalysisLineardiscriminantanalysis(LDA)isoneofthewidelyusedclassi˝cationmethodinhighdimensionalsetting.Forsimplicity,considerthetwoclassclassi˝cationproblem.Givensampleobservationsonclasslabelsandfeatures,LDAclassi˝eseachtestobservationxtoeitherclassk=0ork=1usingtherulek(x)=argmaxkˆxT^^k12^k^^k+log(^ˇk)˙;(9.4.1)where^ˇkistheproportionofclasskobservationsinthetrainingdata,^kisthesamplemeanforclasskonthetrainingdata,and^:=^1isanestimatoroftheinverseofthecommoncovariancematrix.9.4.2GaussianGraphicalModelingGivenarandomvectorY=(Y1;Y2;;Yp),Yi;Yjaresaidtobeconditionallyinde-pendent,iftheirjointdistributiongivenrestofvariablesissameastheproductoftheirconditionallymarginaldistribution.i.e.PYi;YjfYnYi;Yjg=PYifYnYi;YjgPYj;fYnYi;Yjg;wherefYnYi;YjgisvectorYexcludingrandomvariablesYiandYj.53Figure9.2IllustrationofconditionalindependenceInGaussiangraphicalmodel,Y=(Y1;Y2;;Yp)followsmultivariatenormaldistribu-tionwithmeanvectorandvariancecovariancematrix.Inthissetting,theYiandYjareconditionallyindependentgiventherestofvariablesfYk;k6=i;k6=jg,ifpartialcorrelationcoe˚cientbetweenYiandYjiszerowhichisequivalentto1ij=0.Thisparallelrelation-shipbetweenthegraphstructureandtheprecisionmatrixgreatlysimpli˝estheestimationofgraphstructuretotheproblemofsparseprecisionmatrixestimation.54CHAPTER10AJOINTCONVEXPENALTY(JCP)ESTIMATION10.1MotivationDiscussionsfromtheabovechapters5and6suggestestimationbasedonjust`1regu-larizationdoesnotyieldappropriateshrinkageoftheeigenspectrum.Thereforeitisnaturaltoimposeanadditionalpenaltytoobtainawell-conditionedprecisionmatrixestimation.Inthischapter,weextendthelikelihoodbased`1regularizedprecisionmatrixestimationbyaddingtracenormpenalty,andcallitbyJointConvexPenalty(JCP)method.Weimple-mentproximalgradientmethodforcomputingtheproposedestimator.TheproposedalgorithmisshowntoconvergeatarateO(1=k)undermildconditions.10.2JointConvexPenaltyEstimationThemethodtoovercometheover-dispersioninprecisionmatrixwaspreviouslystudiedbymanyauthorswithmajorityofworkfocusingonawell-conditionedestimation.SheenaandGupta[2003]proposeaconstrainedmaximumlikelihoodestimatorwithrestrictionsonthelowerorupperboundoftheeigenvalues.Thismethodfocusesononlyoneofthetwoendsoftheeigen-spectrumandthustheresultingestimatordoesnotcorrectfortheoveresti-mationofthelargeeigenvaluesandunderestimationofthesmalleigenvaluessimultaneously.Consequentlytheirapproachdoesnotaddressthedistortionoftheentireeigen-spectrumespeciallyinhighdimensionalsetting.Wonetal.[2012]consideramaximumlikelihoodestimationofcovariancematrixwithconditionnumberconstraint.Howeverthisapproachitselfrequiresanestimationofconditionnumber.Tocontrolthedistortionofeigen-spectrumofthecovariancematrix,weconsiderajointpenaltyofsumofsingularvalues(tracenorm)inadditionto`1norm.Byminimizing55thejointpenaltyfunctionof`1andtracenorm,theresultingestimatedprecisionmatrixissparseaswellassingularvaluesofthecorrespondingcovariancematrixaremorecenteredthanthoseofthesampleobservedcovariancematrix.10.2.1ProblemFormulationLetX˘Np(0;,˜0:Forthesimplicityofnotation,wedenote=1.Letkkbetracenormandde˝nedassumofsingularvaluesofthematrix.NextwedescribetheproposedJointConvexPenalty(JCP)estimator.10.2.2ProposedEstimatorWeconsiderthefollowingoptimizationproblemwithjointconvexpenalty:argmin˜0F:=f+g1+g2:(10.2.1)wheref=log(det+tr(S;g1=kk1;g2=˝kk(10.2.2)whereand˝arenonnegativeconstants.Theproposedalgorithmtosolveaboveproblemisgiveninchapter11andcanbeusedtosolveawidearrayofproblemsinstatisticsandmachinelearning.SomeoftheotherimportantapplicationsofthismethodincludeMatrixClassi˝cationProblems[TomiokaandAihara[2007],Bach[2008],Multi-TaskLearningAr-gyriouetal.[2008]].Notethatthefin(10.2.2)isaconvexfunction,`1normisasmoothconvexfunctionexceptatoriginandtracenormisconvexsurrogateofrankovertheunitballofspectralnormFazel[2002].Theaboveproblemisaconvexoptimizationproblemwithnon-smoothconstraints.Anaturalchoicetosolvetheaboveoptimizationproblemissubgradientmethodwhichgeneratesasequenceofestimatesfk;k=1;2;3;:::gask=k15Fk1);(10.2.3)56whereissomepositivestepsizeand5Fk1)issub-gradientindicatingdirectionofgreatestvalueincreaseofthefunctionFatk1.Thismethodhasawellknowncon-vergencerateofO(k12)fornon-smoothconvexfunctions(Nesterov[2005]).WeemployproximalgradientmethodtoobtainabetterrateofconvergenceoforderO(k1).Thismethodcanbegeneralizedtosolveanarbitrarycombinationofconvexfunctions[Bertsekas[2010]].57CHAPTER11PROXIMALGRADIENTALGORITHMSANDITSCONVERGENCEANALYSIS11.1IntroductionMuchlikeNewton'smethodisastandardtoolforsolvingunconstrainedsmoothopti-mizationproblemsofmodestsize,proximalgradientalgorithmscanbeviewedasananalo-goustoolfornon-smooth,constrained,large-scale,ordistributedversionsoftheseproblems.Theyareverygenerallyapplicable,butareespeciallywell-suitedtoproblemsofsubstantialrecentinterestinvolvinglargeorhigh-dimensionaldatasets.Themainreasonbehindthesuccessofproximalgradientalgorithmistheavailabilityofinexpensiveoperatorsofsomewellknownfunctions,Theprojectiongradientalgorithminvolvesprojectingapointontoaconvexset,andoftenadmitsaclosedformsolutionthatcanbeobtainedveryquicklywithasimplespecializedmethods.Inoursetup,werequirecomputationofproximaloperatorfor`1andtracenorm.11.2ProximalGradientMethodLethbealowersemi-continuousconvexfunctionof,whichisnotidenticallyequalto+1.Thenproximalpointalgorithm[Rockafellar[1976]]generatesasequenceofsolutionsfk;k=1;2;3;:::gtothefollowingoptimizationproblem,k=Proxhk1)=argmin˜0 h+12kk1k22!:(11.2.1)Thesequencefk;k=1;2;3;:::gweeklyconvergestotheoptimalsolutionofmin˜0h(Rockafellar[1976]).Tousethestructureoftheaboveoptimizationalgorithm,weusequadraticapproximationoff,whichisjusti˝edsincefisstrictlyconvex.5811.3BasicApproximationModelForanyL>0,considerthefollowingquadraticapproximationmodeloffat0:QL;0):=f0)+<0;5f0)>+L2k0k22(11.3.1)where<A,B>istheinnerproductofAandB,andLisapositiveconstant.Theopti-mizationproblemin(11.3.1)hastwoconvexpenalties.Proximalgradientmethodconsistsofsequentialoptimizationof(11.3.1)bytakingoneconstraintatatimeineithercyclicorrandomorder.Rewritingtheoptimizationproblem(11.3.1)withsingleconstraintProx1Lgi0)=argmin˜0 QL;0)+gi!=argmin˜0 L2kf01L5f0)gk22+gi!:(11.3.2)Ingeneral,LisunknownanditisestimatedasanupperboundofLipschitzparameterof5f)[Bachetal.[2011]].InotherwordsLsatis˝es:k5f5f0)k2Lk0k28;02Dom(f):(11.3.3)AcommonmethodofgeneratingvalueofListodoalinesearch.Fortheoptimizationproblem(11.3.2)wesequentiallygeneratenewestimatesandincreasethevalueofLbyafactor>1untilthefollowingconditionismet:fLk)fk1)+<kk1;5fk1)>+L2kkk1k22;(11.3.4)wherekisasolutionatkthiteration.InLemma11.3.1and11.3.2below,wegiveproximalgradientoperatorfor`1andtracenorm.Lemma11.3.1.LetM2Rmn.Theproximaloperatorofk:k1withconstantisgivenbyProx1(M)=argminC˜0 kCk1+12kCMk22!;(11.3.5)where59Prox1(M)=sign(M)(0;abs(M))+,>0.andabs(M)entriwisemaximumfunctionformatrixM.Proof.Proofofthelemmaisgiveninappendix.Lemma11.3.2.LetM2RmnandM=UVTbesingularvaluedecompositionofMwhereU2RmrandV2Rrnhaveorthogonalcolumns,isthediagonalmatrixofsingularvaluesofMandrisrankofmatrixM.Thenproximaloperatorofk:kwithconstant˝isgivenbyProx˝(M)=argminC˜0 ˝kCk+12kCMk22!;(11.3.6)whereProx˝(M)=U˝VT,˝isdiagonalmatrixwith˝))ii=max(0;ii˝),˝>minipii).Proof.Proofofthelemmaisgiveninappendix.For`1andtracenorm,theproximaloperatorsareinexpensivetocalculate.Thisre-sultsine˚cientoptimizationoftheobjectivefunction.Theproximaloperatorfor`1iselementwisesoft-thresholdingoperator.Theproximaloperatorfortracenormisobtainedbyshrinkingthesingularvaluesofprecisionmatrixbyregularizationparameter.Alargervalueofregularizationparameterresultsismoreshrinkageoftheeigenvalues.11.4AlgorithmforoptimizationBelowwesummarizetheoptimizationalgorithmfor(11.3.2).InitializeL0=1;>1;0=diag(1=diag(S)).Iterate:Step1:SetL=Lk1.60Step2:WhileF)>QL);k1)+g1)+g2)(where=argmin˜0QL;k1)+g1+g2)SetL=L.Step3:SetLk=L;SetZk=k1Lk5fk);SetZk+1=Prox(˝=Lk)g2(Zk);Setk+1=Prox(k)g1(Zk+1):Repeatuntilconvergence.11.5ChoosingtheRegularizationParameterThechoiceofregularizationparameterisachallengingprobleminhighdimensionaldataanalysis.Regularizationhasclearbene˝tinproducingsparsesolutionaswellreducesfalsediscoveryrate.Asmallervalueofaccountsforasparserstructureoftheprecisionmatrix.SomeofthemethodsforchoosingregularizationparameterincludeK-foldcrossvalidation(KCV),stabilityapproachtoregularizationselection(StARS)(Liuetal.[2010]),Akaike-InformationCriteria(AIC)andBayesianInformationcriteria(BIC).Experiments[Liuetal.[2010]]haveshownthatAICandBICmethodstendtogivepoorperformanceforsmallersamplesizes.AlsoK-foldcrossvalidationtendstoselectsmallervaluesofregular-izationparameterandresultsinhigherfalsediscoveryrate.WefollowStARSapproachforestimatingtheregularizationparametersand˝:Nextwepresenttheconvergenceanalysisofthealgorithmgivenin11.4.11.6ConvergenceAnalysisWeusetheproposition3.1fromBertsekas[2010],Lemma11.6.1.Letfn;Ln;n=1;2;:::gbethesequencegeneratedbyalgorithmgivenin61Letc>0beaconstantsatisfying,maxf5kfk;5kgjkgcandmaxffn)f(Zn+j1);gjn)gj(Zn+j1)gcknZn+j1k;j=1;2:Thenforacyclicorderoptimizationofcomponentsg1()andg2(),followingholds:kn1k2knk22Ln(Fn)F))+18c2=L2n(11.6.1)wheren1=Prox(˝=Ln)g2n)andisasolutionof(2.1)and(2.2).Lemma11.6.2.Letfn;Ln;n=1;2;:::gbethesequencegeneratedbyalgorithmx11:4.Letcbeaconstantasde˝nedinLemma(11.6.1).Then,Fn)F)Ln4 kn1k2knk2!+9c22Ln2<n;5knk1>(11.6.2)Proof.ProofoftheLemma11.6.2isgiveninappendix.Wegivebelowtheconvergenceofthealgorithmx11.4.Theorem11.6.1.Letfk;k=1;2;:::gbethesequencegeneratedbyalgorithmx11:4.Letcbeaconstantasde˝nedinLemma11.6.1.Inaddition,weassumethatthereexistsaconstantM<1suchthat1Pn=1j<n;5knk1>j<M,thenFk)F) Lk0k2F+18c2+M4k!;(11.6.3)whereL>0istheleastupperLipschitzconstantofthegradientoffand>1isconstantasde˝nedinalgorithmx11:4.Proof.NotethatLnLLn,foralln=1;2;:::UsingLemma11.6.2,byaddingFn)F)overn=1;2;:::,wegetthedesiredresult.62Duetothenon-smoothnessofthetracenorm,theoptimal˝rstorderblackboxmethodshaveconvergencerateofO(k12).Theproximalgradientalgorithmusesthespecialstructureofthetraceand`1normwhichimprovestheconvergencerateforjointpenaltytotheorderO(k1).11.7SimulationStudyToimplementtheproposedmethod,weperformasimulationstudyforvariouschoicesofprecisionmatrix.Weconsiderdi˙erenttypesofprecisionmatricesasoutlinedinchapter8.Hereweconsidertheunderlyingtrueprecisionmatrixsparseratherthanthecovariancematrix.Forallthesechoicesofinversecovariancematrices,wegeneraterandomnumbersfrommultivariatenormaldistributionwithvaryingnandp.Wesetn=50;100;200andp=50;100;200.TheperformanceofproposedmethodiscomparedtographicallassoandSPICEestimatesofprecisionmatrix.ThejointconvexpenaltyestimateoftheprecisionmatrixwascomputedusingRsoftwareversion3.0.1basedonthealgorithmx11:4.ThegraphicallassoestimateoftheprecisionmatrixwascomputedusingRpackage(http://statweb.stanford.edu/tibs/glasso/).Inthereisoptionofnotpenalizingthediagonalelementsbysettingtheoptionenalizing.diagonal=FALSE",thisgivesSPICEestimate.11.7.1PerformanceCriteriaForeachofprecisionmatrixestimate,wecalculateKullback-Leibler(KL)Loss,andAverageRelativeError(ARE)de˝nedbelow:KLLoss;^=log(det(^+tr1^+log(detpARE;^=jlog(f(S;^log(f(S;j=log(f(S;63wheref(;)isdensityofmultivariateGaussiandistributionandSissamplecovariancematrix.Thetuningparametersand˝wereestimatedfollowingtheLiuetal.[2010]criteria,whichwedescribebelow.11.7.2StARSMethodofTuningparameterselection:Givenasampleofsizen,methodgeneratesNsamplesofsizeb,whereb<n.Inoursettingforn<200,wechooseb=0.8nandforn200,b=10pn.ForeachoftheseNsamples,anestimateofprecisionmatrixisobtained.Foreachentryoftheprecisionmatrix,ameasureofinstabilityiscalculatedbasedonallNestimates.Finallyaregularizationparameterisselectedwhichminimizestheaverageinstabilityoverallpossibleentriesoftheestimatedprecisionmatrix.Inpracticethismethodtendstoselectleastamountofregularizationparameterthatsimultaneouslymakesestimatesoftheprecisionmatrixsparseandreplicableunderrandomsampling.StARSisusedforestimatingthepenalizationparameterforallthecompetingmethodsasgiveninsimulationviz.JCP,GraphicalLassoandSPICE.11.7.3SimulationResultsThesimulationresultsaregivenintables11.7.2.1-11.7.2.4.Thenumbersinbracketarestandarderroroftheestimatebasedon20simulations.ForToeplitztypeprecisionmatrix,theproposedmethodoutperformsothertwointermsofMSEandKL-lossforsmalln.Forlargesamplesize,theGraphicallassotendstogivebetterperformancethanothermethods.Thismaybeduetofactthatgraphicallassosolvestheconstrainedquadraticregressionproblem(dualofconstrainedobjectivefunction).Forlargesamplesize,theregressioncoe˚-cientstendstoapproximatewellthetruevalueswhichresultsinbetterestimateofprecisionmatrix.6411.7.3.1ToeplitzTypePrecisionMatrixTable11.1AverageKL-Losswithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty6.15(0.065)5.28(0.032)4.829(0.0212)GraphicalLasso6.78(0.066)5.17(0.054)4.31(0.049)SPICE6.181(0.069)5.28(0.037)4.79(0.03)p=100MixedPenalty12.37(0.0806)10.84(0.0444)9.922(0.0209)GraphicalLasso14.01(0.079)11.21(0.061)8.85(0.041)SPICE12.38(0.076)11.08(0.043)10.03(0.078)p=200MixedPenalty25.35(0.117)22.2(0.062)20.39(0.021)GraphicalLasso31.5(0.124)25.05(0.09)18.62(0.035)SPICE25.4(0.122)22.25(0.07)20.66(0.039)Table11.2Averagerelativeerrorwithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty0.0857(0.0075)0.0402(0.005)0.1106(0.0041)GraphicalLasso0.135(0.01)0.03(0.004)0.0629(0.004)SPICE0.038(0.005)0.067(0.005)0.118(0.004)p=100MixedPenalty0.0722(0.0054)0.014(0.003)0.1003(0.0031)GraphicalLasso0.24(0.006)0.125(0.003)0.024(0.003)SPICE0.031(0.01)0.0802(0.004)0.1235(0.004)p=200MixedPenalty0.1182(0.0032)0.009(0.0012)0.09(0.001)GraphicalLasso0.456(0.003)0.264(0.002)0.04(0.0014)SPICE0.013(0.002)0.1032(0.0032)0.1313(0.0014)6511.7.3.2BlockTypePrecisionMatrixTable11.3AverageKL-Losswithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty3.624(0.067)2.809(0.026)2.46(0.0165)GraphicalLasso4.219(0.0652)2.943(0.0321)2.146(0.025)SPICE3.646(0.0644)2.79(0.031)2.37(0.0262)p=100MixedPenalty7.53(0.081)6.01(0.035)5.232(0.0163)GraphicalLasso9.477(0.0717)6.791(0.0485)4.69(0.031)SPICE7.602(0.0856)6.063(0.0413)5.002(0.0372)p=200MixedPenalty15.05(0.123)12.34(0.0443)10.8(0.0264)GraphicalLasso22.46(0.2643)16.29(0.0807)10.42(0.059)SPICE15.44(0.1296)12.27(0.056)10.93(0.0444)Table11.4Averagerelativeerrorwithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty0.0865(0.005)0.013(0.003)0.0641(0.0024)GraphicalLasso0.1891(0.0076)0.092(0.0044)0.0085(0.001)SPICE0.029(0.006)0.0279(0.0036)0.0669(0.0019)p=100MixedPenalty0.1132(0.0039)0.019(0.002)0.0623(0.0013)GraphicalLasso0.3732(0.0043)0.2131(0.0028)0.0345(0.001)SPICE0.028(0.003)0.0458(0.0048)0.0729(0.0022)p=200MixedPenalty0.1844(0.0064)0.048(0.003)0.048(0.002)GraphicalLasso0.715(0.0171)0.4275(0.0023)0.1248(0.0014)SPICE0.07(0.004)0.0493(0.0051)0.0904(0.0013)6611.7.3.3HubGraphTypePrecisionMatrixTable11.5AverageKL-Losswithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty2.813(0.0618)1.74(0.027)1.066(0.0175)GraphicalLasso3.325(0.0479)1.962(0.022)1.098(0.0242)SPICE2.697(0.056)1.89(0.0373)1.399(0.0213)p=100MixedPenalty6.576(0.115)4.412(0.0545)2.91(0.029)GraphicalLasso7.993(0.1087)5.524(0.0693)3.232(0.038)SPICE5.59(0.0792)4.25(0.03)3.46(0.0317)p=200MixedPenalty10.06(0.1076)5.862(0.0628)3.163(0.0292)GraphicalLasso16.36(0.0978)11.66(0.0585)5.605(0.0358)SPICE6.88(0.098)4.53(0.072)2.83(0.0267)Table11.6Averagerelativeerrorwithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty0.0795(0.0031)0.0421(0.002)0.012(0.001)GraphicalLasso0.0786(0.0043)0.049(0.003)0.016(0.001)SPICE0.0103(0.001)0.01(0.001)0.0137(0.0008)p=100MixedPenalty0.137(0.005)0.0714(0.0031)0.021(0.0005)GraphicalLasso0.177(0.006)0.1001(0.004)0.036(0.001)SPICE0.023(0.001)0.008(0.001)0.016(0.001)p=200MixedPenalty0.229(0.0014)0.1274(0.0008)0.0415(0.0003)GraphicalLasso0.343(0.003)0.2121(0.001)0.075(0.0004)SPICE0.06(0.001)0.034(0.001)0.003(0.0003)6711.7.3.4NeighborhoodGraphTypePrecisionMatrixTable11.7AverageKL-Losswithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty3.33(0.084)2.28(0.036)1.47(0.031)GraphicalLasso3.73(0.074)2.45(0.045)1.4(0.038)SPICE3.37(0.095)2.723(0.0556)2.274(0.0576)p=100MixedPenalty6.507(0.165)4.351(0.0472)2.43(0.042)GraphicalLasso7.846(0.117)5.472(0.0601)2.598(0.0395)SPICE5.49(0.12)4.21(0.038)3.083(0.0576)p=200MixedPenalty13.04(0.1767)7.697(0.0877)4.16(0.056)GraphicalLasso18.39(0.1129)12.69(0.0815)5.79(0.0639)SPICE9.3(0.131)6.094(0.06)4.359(0.0723)Table11.8Averagerelativeerrorwithstandarderrorover20replicationsn=50n=100n=200p=50MixedPenalty0.0714(0.0029)0.0323(0.0015)0.008(0.001)GraphicalLasso0.068(0.0041)0.0363(0.0022)0.0128(0.0013)SPICE0.0056(0.001)0.018(0.002)0.0276(0.0014)p=100MixedPenalty0.139(0.006)0.069(0.0035)0.0235(0.0005)GraphicalLasso0.184(0.005)0.096(0.005)0.0386(0.001)SPICE0.019(0.002)0.005(0.001)0.016(0.001)p=200MixedPenalty0.2402(0.0029)0.131(0.0013)0.0428(0.0006)GraphicalLasso0.3591(0.0033)0.2156(0.0018)0.086(0.002)SPICE0.051(0.001)0.01(0.003)0.009(0.0004)6811.8SummaryForblocktypeprecisionmatrix,theproposedmethoddominatesthegraphicallassoandSPICEintermsofKL-lossforsmallsamplesizes.AlsoJointPenaltyandSPICEmethodshavebetterperformancethangraphicallassointermsofaveragerelativeerror.ThecorrespondingstandarderrorestimatesforJointPenaltyandSPICEaresmallerthanthoseofgraphicallasso.Forneighborhoodgraphtypeprecisionmatrix,theJointPenaltymethodhasbetterperformanceforsmallpintermsofKL-loss.HoweverSPICEseemstoperformbetterthanothermethodsforsmallnandlargepintermsofAverageRelativeError.OveralltheJointPenaltymethodhasbetterperformancethanotherthetwomethodsforsmallsamplesizeandforallchoicesoftheprecisionmatrix.Theperformanceofallthreemethodsmethodsvariesovertwochoicesoflossfunctionsastheyhavedi˙erentformulas.ThevalueofKL-lossandaveragerelativeerroraresubstantiallydi˙erentfordi˙erentchoicesoftheprecisionmatrix.Thisshowsthaterrorinestimatesdependsuponthestructureofunderlyingtrueprecisionmatrix.For˝xedptheestimatestendstoimproveforincreasingsamplesizes.Howeverfor˝xedn,asexpected,theestimatorsperformancegoesdownwithincreasingp.Foradditionalsimulations,refertoMaurya[2014].11.9DiscussionTheproposedmethodusesajointpenaltywhichismore˛exibleforpenalizingtheentriesofprecisionmatrixindi˙erentfashionthantheGraphicalLassoandtheSPICE.Theproposedproximalgradientmethodcanbeextendedtoproblemswhereonehasanarbitrarynumberofconvexpenaltyconstraints.Undermildconditionsthealgorithmachievessub-linearrateofconvergencewhichmakesitattractivechoiceformanyoptimizationproblems.Thesimulationstudyshowsthattheproposedmethodperformsbetterthantheothertwomethodsforsmallsamplesizesandlargep.Thesimulatedexamplesshowthatperformance69oftheJCPestimatesofprecisionmatrixareconsistentasAREandKL-Lossdecreasesrapidly(aswellasthecorrespondingstandarderrors)forincreasingsamplesize.70CHAPTER12SIMULTANEOUSESTIMATIONOFSPARSEANDWELL-CONDITIONEDPRECISIONMATRIXInthischapter,wediscusstheestimationoftheprecisionmatrixbasedonFrobeniusnormlossfunction.12.1MotivationEstimationofprecisionmatrixbasedontheFrobeniusnormlossfunctionisnotwellde˝nedproblemunlessSiswell-conditined.Onewaytosolvethisproblemisto˝rstob-tainagoodestimatorforthecovariancematrix,thenusethisinplaceofS.Wefollowthisapproachandconsideratwostepprocedureofprecisionmatrixestimation.Themainadvantagesofthisapproachare:(i)itallowstointroducesparsityintheprecisionmatrixitselfthaninthecovariancematrix,and(ii)thevarianceofeigenvaluespenaltyallowstoperformoptimalshrinkageintheeigenspectrumoftheprecisionmatrix.Asthecaseincovariancematrixestimation,theresultingoptimizationproblemisconvex,whichintrunyieldsanexactalgorithmwithlowcomputationalcomplexity.InthischapterweextendtheJPENapproachtoestimateawell-conditionedandsparseprecisionmatrix.Similartothecovariancematrixestimation,weproposeanestimatorfortheprecisionmatrixbasedonregularizedinversecorrelationmatrixanddiscussitsrateofconvergenceinbothFrobeniousandoperatornorm.Notation:WeshalluseZandforinversecorrelationandprecisionmatrixrespectively.7112.2JointPenaltyEstimation:ATwoStepApproachLet^RKbeaJPENestimatorforthetruecorrelationmatrix.ByLemma6.3.2,^RKispositivede˝niteandwell-conditioned.De˝netheJPENestimatorofinversecorrelationmatrixasthesolutiontothefollowingoptimizationproblem,^ZK=argminZ=ZT;tr(Z)=tr(^R1K)j()2^SK2kZ^R1Kk2+kZk1+pXi=1f˙i(Z)˙(Z)g2(12.2.1)where^SK2=ˆ():>0;slogpn;8>0;˙minf(^R1K+t1I)2sign(^R1K+t1I)g>˙;andt1istheaverageofthediagonalelementsof^R1K.Theminimizationin(12.2.1)overZisfor˝xed()2^SK2.TheproposedJPENestimatoroftheprecisionmatrix(basedonregularizedinversecorrelationmatrixestimator^ZK)isgivenby,^K=(S+)1=2^ZK(S+)1=2;whereS+isthediagonalmatrixofthediagonalelementsofS.Moreover(12.2.1)isaconvexoptimizationproblemand^ZKispositivede˝nite.Nextwegiveanotherestimateoftheprecisionmatrixbasedon^Sof6.3.3.Considerthefollowingoptimizationproblem:^S=argminT;trtr(^1S)j()2^SS2jj^1Sjj2F+kk1+pXi=1f˙i˙g2;(12.2.2)where^SS2=ˆ():>0;slogpn;8>0;˙minf(^1S+t2I)2sign(^1S+t2I)g>˙;72andt2isaverageofthediagonalelementsof^S.Theminimizationin(12.2.2)overisfor˝xed()2^SS2.Theestimatorin(12.2.2)ispositivede˝niteandwell-conditioned.12.3WeightedJPENestimatorforprecisionmatrixSimilartoweightedJPENcovariancematrixestimator^K;A,aweightedJPENestima-toroftheprecisionmatrixisobtainedbyaddingpositiveweightsaitotheterm(˙i(Z)1)2in(12.2.2).TheweightedJPENprecisionmatrixestimatoris^K;A:=(S+)1=2^ZA(S+)1=2,where^ZA=argminZ=ZT;tr(Z)=tr(^R1K)j()2^SK;A2jjZ^R1Kjj2F+kZk1+pXi=1aif˙i(Z)1g2;(12.3.1)with^SK;A2=ˆ():qlogpn;(2˙min(R1K))(1+t1max(Aii)1)(1+min(Aii)1p+min(Aii)p˙;andA=diag(A11;A22;App)withAii=ai.Theoptimizationproblemin(12.3.1)isconvexandyieldsapositivede˝niteestimatorfor()2^SK;A2.12.4TheoreticalAnalysisofJPENestimatorsInthissection,wederivetherateofconvergenceoftheproposedJPENestimatorsofprecisionmatrixinbothFrobeniusandspectralnorm.Let0=10bethetrueprecisionmatrix.LetX=(X1;X2;;Xn)besub-Gaussianrandomvectorsasde˝nedin(6.3.1).Wemakethefollowingadditionalassumptionsaboutthe0.B0.SameastheassumptionA0ofx6:3.B1.WithH=f(i;j):0ij6=0;i6=jg,thejHjs,forsomepositiveintegers.B2.Thereexist0<k<1largeenoughsuchthat(1=k)˙min0)˙max0)k.Thenexttheoremgivesconsistencyof^ZKand^K.73Theorem12.4.1.UnderAssumptionsB0,B1,B2andfor()2^SK2,k^ZKR10kF=OPsslogpnandk^K0k=OPs(s+1)logpn(12.4.1)whereR10istheinverseoftruecorrelationmatrix.Remark12.4.1.NotethattheJPENestimator^Kachievesmini-maxrateofconvergencefortheclassofcovariancematricessatisfyingassumptionB0,B1,andB2andthereforeoptimal.ThesimilarrateisobtainedinCaietal.[2015]fortheirclassofsparseinversecovariancematrices.Thenexttheoremgivesconsistencyof^S.Theorem12.4.2.Let()2^SS2.UnderAssumptionsB0,B1,andB2,the^Sof(12.2.1)satis˝es,k^S0kF=OPs(s+p)logpn:(12.4.2)Consistencyof^ZA:Asimpleexerciseshowsthattheestimator^ZAhassimilarrateofconvergenceasthatof^ZK.74CHAPTER13SIMULATIONSANDANAPPLICATIONTOREALDATAANALYSISInthischapter,wecomparetheperformanceoftheproposedmethodtovariousothermethodsforanumberofstructuredprecisionmatrices.13.1SimulationResults:SettingsWechosesimilarstructuresasinchapter8forprecisionmatrix0i.e.wereplace0by0andforallthesechoicesofinversecovariancematrices,generaterandomvectorsfrommultivariatenormaldistributionswithvaryingnandp.Wechosen=50;100andp=500;1000.Theperformanceofproposedcovariancematrixestimator^Kiscomparedtothefollowingmethods.Graphicallasso[Friedmanetal.[2008]]:Graphicallassoestimatesasparsepre-cisionmatrix.Hereweinverttheinverse,andincludeinouranalysis.Theestimatewascomputedusing`R'package'Glasso'.Formoredetails,refertohttp://statweb.stanford.edu/tibs/glasso/.SparsePermutationInvariantEstimation(SPICE)[Rothman[2012]]:theSPICEwascomputedusing`R'package'Glasso'wherewedonotpenalizediago-nalelements.Formoredetails,referto(http://cran.r-project.org/web/pack-ages/PDSCE/index.html)Constrained`1minimizationforinversecovariancematrices(CLIME)[Caietal.[2011]]:Theirestimatewascomputedusing`R'package'clime'.Formoredetailssee(https://cran.r-project.org/web/packages/clime/index.html).AllthecomputationsweredoneusingstatisticalsoftwareRonaAMD2.8GHzprocessor.75Table13.1PrecisionmatrixestimationBlocktypeprecisionmatrixn=50n=100p=500p=1000p=500p=1000Glasoo4.144(0.523)1.202(0.042)0.168(0.136)1.524(0.028)PDSCE1.355(0.497)1.201(0.044)0.516(0.196)0.558(0.032)CLIME4.24(0.23)6.56(0.25)6.88(0.802)10.64(0.822)JPEN1.248(0.33)1.106(0.029)0.562(0.183)0.607(0.03)Hubtypeprecisionmatrixn=50n=100p=500p=1000p=500p=1000Glasoo1.122(0.082)0.805(0.007)0.07(0.038)0.285(0.004)PDSCE0.717(0.108)0.702(0.007)0.358(0.046)0.356(0.005)CLIME10.5(0.329)10.6(0.219)6.98(0.237)10.8(0.243)JPEN0.684(0.051)0.669(0.003)0.34(0.024)0.337(0.002)Neighborhoodtypeprecisionmatrixn=50n=100p=500p=1000p=500p=1000Glasoo1.597(0.109)0.879(0.013)1.29(0.847)0.428(0.007)PDSCE0.587(0.13)0.736(0.014)0.094(0.058)0.288(0.01)CLIME10.5(0.535)11.5(0.233)10.5(0.563)11.5(0.245)JPEN0.551(0.075)0.691(0.008)0.066(0.042)0.201(0.007)13.2PerformanceComparisonForeachoftheprecisionmatrixestimate,wecalculateAverageRelativeError(ARE)basedon50iterationsbasedontheformulagivenin(8.2.1).Theoptimalvaluesoftuningparameterswereobtainedoveragridofvaluesbyminimizing5foldcross-validationasexplainedinx7:1:2.TheJPENestimator^Koutperformsothermethodsforthemostofthechoicesofnandpforall˝vetypesoftheprecisionmatrices.Additionalsimulations(notincludedhere)showthatfornˇp,alltheunderlyingmethodsperformsimilarlyandtheestimatesoftheireigenvaluesarealsowellalignedwithtruevalues.Howeverinhighdimen-sionalsetting,forlargepandsmalln,theirperformanceisdi˙erentasseeninsimulationsofTable13.1andTable13.2.76Table13.2PrecisionmatrixestimationToeplitztypeprecisionmatrixn=50n=100p=500p=1000p=500p=1000Glasoo2.862(0.475)2.89(0.048)2.028(0.267)2.073(0.078)PDSCE1.223(0.5)1.238(0.065)0.49(0.269)0.473(0.061)CLIME4.91(0.22)7.597(0.34)5.27(1.14)8.154(1.168)JPEN1.151(0.333)2.718(0.032)0.607(0.196)2.569(0.057)Cov-Itypeprecisionmatrixn=50n=100p=500p=1000p=500p=1000Glasoo54.0(0.19)190.(5.91)14.7(0.37)49.9(0.08)PDSCE28.8(0.19)45.8(0.32)16.9(0.04)26.9(0.08)CLIME59.8(0.82)207.5(3.44)15.4(0.03)53.7(0.69)JPEN26.3(0.36)7.0(0.07)15.7(0.08)23.5(0.3)13.3ColonTumorGeneExpressionDataAnalysisInthissection,wecomparetheperformanceofJPENestimatorofprecisionmatrixfortumorclassi˝cationusingLinearDiscriminantAnalysis(LDA).Thegeneexpressiondata(Alonetal.[1999]consistsof40tumorousand22non-tumorousadenocarcinomatis-sues.Afterpreprocessing,datawasreducedtoasubsetof2,000geneexpressionvalueswiththelargestminimalintensityoverthe62tissuesamples(source:http://genomics-pubs.princeton.edu/oncology/a˙ydata/index.html).Figure13.1showsthegeneexpressiondataforbothtumorousandnon-tumoroustissues.Tumoroustissuesaremarkedwitharrowsontheleftwhereasnormaltissuesareunmarked.Thereisclearseparationofnormalandtumortissues.77Figure13.1ColontumorgeneexpressiondataInouranalysis,wereducedthenumberofgenesbyselectingpmostsigni˝cantgenesbasedonell1regularizedlogisticregression.Weobtainestimatesofprecisionmatrixforp=50;100;200andthenusetheLDAtoclassifythesetissuesaseithertumorousornon-tumorous(normal).Classifyeachtestobservationxtoeitherclassk=0ork=1usingtheLDArulek(x)=argmaxkˆxT^^k12^k^^k+log(^ˇk)˙;(13.3.1)where^ˇkistheproportionofclasskobservationsinthetrainingdata,^kisthesamplemeanforclasskonthetrainingdata,and^:=^1isanestimatoroftheinverseofthecommoncovariancematrixcomputedfromthetrainingdata.Tuningparametersandwerechosenusing5-foldcrossvalidation.Tocreatetrainingandtestsets,werandomlysplitthedataintoatrainingandtestsetofsizes42and20,respectively;followingtheapproachusedbyWangetal.[2007],thetrainingsethas27tumoroustissuesand15non-tumoroustissues.Sincewedonothaveseparatevalidationset,wedothe5-foldcrossvalidationon78Table13.3Averagesandstandarderrorsofclassi˝cationerrorsover100replicationsin%.Methodp=50p=100p=200LogisticRegression21.0(0.84)19.31(0.89)21.5(0.85)SVM16.70(0.85)16.76(0.97)18.18(0.96)NaiveBayes13.3(0.75)14.33(0.85)14.63(0.75)GraphicalLasso10.9(1.3)9.4(0.89)9.8(0.90)JointPenalty9.9(0.98)8.9(0.93)8.2(0.81)trainingdataasfollowing.Ateachsplit,wedividethetrainingdatainto5subsets(fold)where4subsetsareusedtoestimatetheprecisionmatrixandonesubsetisusedtomeasuretheclassi˝er'sperformance,andforeachsplit.Thisprocedureisrepeated5timesbytakingoneofthe5subsetsasvalidationdata.Anoptimalcombinationofandisobtainedbyminimizingthe5-foldcrossvalidationerror.Theaverageclassi˝cationerrorswithstandarderrorsoverthe100splitsarepresentedinTable13.3.Sincethesamplesizeislessthanthenumberofgenes,weomittheinversesamplecovariancematrixasitisnotwellde˝nedandinsteadincludethenaiveBayes',andsupportvectormachineclassi˝ers.NaiveBayeshasbeenshowntoperformbetterthanthesamplecovariancematrixinhigh-dimensionalsettings(BickelandLevina[2004].SupportVectorMachine(SVM)isanotherpopularchoiceforhighdimensionalclassi˝cation.Amongallthemethodsintable13.3,theprecisionmatrixbasedLDAclassi˝ersperformfarbetterthanNaiveBayes,SVMandLogisticRegression.Forallotherclassi˝erstheclassi˝cationperformancedeterioratesforincreasingp.Forlargerp,i.e.,whenmoregenesareaddedtothedataset,theperformanceofJPENestimatebasedLDAclassi˝erinitiallyimprovesbutitdeterioratesforlargep.Forp=2000,whenallthegenesareusedinanalysis,theclassi˝erbasedonprecisionmatrixhasaccuracyof30%.Thisisduetothefactthatasdimensionofcovariancematrixincreases,theestimatordoesnotremainveryinformative.79Figure13.2PartialcorrelationnetworkofcolontumorgeneexpressiondataToseetheunderlyinggene-geneinteractionininversecorrelationmatrix,wedothenetworkplotthegenesthatwithnon-zeropartialcorrelationintheJPENestimatedinversecorrelationmatrices.Thenetworkgraphshowsthesparsestructureunderlyingthecolontumordata,astherearefewconnectededgesoutoftotal1225possibleedges.80APPENDICES81APPENDIXAJPENCOVARIANCEMATRIXESTIMATIONProof.[Lemma6.3.1]Letf(R)=kRKk2+kRk1+pXi=1f˙i(R)˙Rg2:(A.0.1)where˙RisthemeanofeigenvaluesofR.Duetotheconstrainttr(R)=p,wehave˙R=1.Thethirdtermof(A.0.1)canbewrittenaspXi=1f˙i(R)˙Rg2=tr(R2)2tr(R)+pThen,f(R)=tr(R2)2tr(RK)+tr(K2)+kRk1+ftr(R2)2tr(R)+pg=tr(R2(1+))2tr(K+I)+tr(K2)+kRk1+p=(1+)kR(K+I)=(1+)k2+tr(K2)+kRk1+p(A.0.2)ThisisquadraticinRwithan`1penaltytotheo˙-diagonalentriesofR,thereforeaconvexfunctioninR.Proof.[ProofofLemma6.3.2]Thesolutionto(A.0.2)satis˝es:2(R(K+I))(1+)1+@kRk1@R=0(A.0.3)where@kRk1@Risgivenby:@kRk1@R=8>>>>>><>>>>>>:1:ifRij>01:ifRij<0˝2(1;1):ifRij=0NotethatkRk1hasthesamevalueirrespectiveofsignofR,thereforetherighthandsideof(A.0.2)isminimumif:82sign(R)=sign(K+I)=sign(K)forevery>0,using(A.0.3),˙minf(K+I)2sign(K)g>givesa()2^SK1andsuchachoiceof()guaranteestheestimatortobepositivede˝nite.RemarkA.0.1.Intuitively,alargershrinkstheeigenvaluestowardscenterwhichis1,alargerwouldresultinpositivede˝niteestimator,whereasalargerresultsinsparseestimate.Acombinationof()resultsinasparseandwell-conditionedestimator.Inparticularcase,whenKisdiagonalmatrix,the<2.Proof.[Theorem6.3.1]LetQ(R)=f(R)f(R0),whereR0isthetruecorrelationmatrixandRisanyothercorrelationmatrix.LetR=UDUTbeeigenvaluedecompositionofR,whereDisdiagonalmatrixofeigenvaluesandUismatrixofeigen-vectors.R0=U0D0UT0iseigenvaluedecompositionofR0.Wehave,Q(R)=kRKk2F+kRk1+tr(D22D+p)kR0Kk2FkR0k1tr(D202D0+p)(A.0.4)Letn(M):=f:=T;kk2=Mrn;0<M<1g.Theestimate^RminimizestheQ(R)orequivalently^=^RR0minimizestheG=Q(R0+.NotethatGisconvexandif^isitssolution,thenG(^G(0)=0.Therefore,ifweshowthatGisnon-negativefor2n(M),then^willbewithinsphereofradiusMrn.Werequirern=oq(p+s)logp=n.Consider,kRKk2FkR0Kk2F=tr(RTR2RTK+KTK)tr(RT0R02R0S+KTK)=tr(RTRRT0R0)2tr((RR0)TK)=tr((R0+T(R0+RT0R0)2trTK)=trT2trT(KR0))83Next,weboundterminvolvingKintheaboveexpression.WehavejtrT(R0K))jXi6=jjij(R0ijKij)jmaxi6=j(jR0ijKijj)kk1C0(1+˝)slogpnkk1C1slogpnkk1holdswithhighprobability,byaresult(Lemma1fromRavikumaretal.[2011])onthetailinequalityforthesamplecovariancematrixofsub-GaussianrandomvectorsandwhereC1=C0(1+˝);C0>0.Nextweobtainanupperboundonthetermsinvolvingin(A.0.4).ByCauchy-Schwarzinequality,tr(D22D)tr(D202D0)=trfR2R20g2trfRR0)g=tr(R0+2tr(R20)=2tr(R0+trT2pskkF+kk2F:Toboundtheterm(kRk1kR0k1)=(k+R0k1kR0k1),letEbetheindexsetasde˝nedinAssumptionA.2ofTheorem6.3.1.Thenusingthetriangleinequality,weobtain,(k+R0k1kR0k1)=(kE+R0k1+kEk1kR0k1)(kR0k1kEk1+kEk1kR0k1)(kEk1kEk1)Let=(C1)qlogp=n,=(C11)qlogp=n;where()2^SK1,weobtain,GtrT+)2C1ˆslogpn(kk1)+11sslogpnkkF˙+C1slogpnkEk1Ek1kk2F(1+)2C1slogpnkEk1+kEk1+C1slogpnkEk1Ek12C11sslogpnkkF:84AlsobecausekEk1=P(i;j)2E;i6=jijpskkF,2C1slogpnkEk1+C1slogpnkEk1slogpnkEk12C1+C10forsu˚cientlysmall.Therefore,Gkk2F1+C11slogpnC1sslogpnk+kFf1+1+21gkk2F1+C11slogpnC1Mf1+1+21g0;forallsu˚cientlylargenandM.Whichprovesthe˝rstpartoftheorem.Toprovetheoperatornormconsistency,bysub-multiplicativenormpropertykABkkAkkBk,k^K0k=k^W^R^WWKWkk^WWkk^RKkk^WWk+k^WWk(k^RkkWk+k^WkkKk)+k^RKkk^WkkWk:SincekKk=O(1)andk^RKkF=O(qslogpn)thesetogetherimpliesthatk^Rk=Op(1).Also,k^W2W2k=maxkxk2=1pXi=1j(^w2iw2i)jx2imax1ipj(^w2iw2i)jpXi=1x2i=max1ipj(^w2iw2i)j=Opslogpn;byusingaresult(Lemma1fromRavikumaretal.[2011]).Nextweshallshowsthatk^WWkk^W2W2k,(whereABmeansA=OP(B)andB=OP(A)).Wehave,k^WWk=maxkxk2=1pXi=1j(^wiwi)jx2i=maxkxk2=1pXi=1j^w2iw2i^wi+wijx2ipXi=1j(^w2iw2i)jx2i=C3k^W2W2k:85wherewehaveusedthefactthatthetruestandarddeviationsarewellabovezero,i.e.,90<C3<1suchthat1=C3w1iC38i=1;2;;p,andthesamplestandarddevia-tionsareallpositive,i.e,^wi>08i=1;2;;p:Nowsincek^W2W2kk^WWk,itfollowsthatk^Wk=Op(1)andwehavek^K0k2=Opslogpn+logpn.Thiscompletestheproof.Proof.[Theorem3.2]Letf=jjSjj2F+kk1+pXi=1f˙i˙g2;Similartotheproofoftheorem6.3.1,de˝nethefunction,Q1=ff0)where0isthetruecovariancematrixandisanyothercovariancematrix.Let=UDUTbetheeigenvaluedecompositionof,whereDisdiagonalmatrixofeigenvaluesandUismatrixofeigen-vectors.Let0=U0D0UT0iseigenvaluedecompositionof0.Then,Q1=kSk2F+kk1+tr(D2)(tr(D))2=pk0Sk2Fk0k1tr(D20)(tr(D0))2=p(A.0.5)whereA=diag(a1;a2;;ap).Write=0,andletn(M):=f:=T;kk2=Mrn;0<M<1g.Theestimate^minimizesQorequivalently^=^0minimizesG=Q0+.NotethatGisconvexandif^beitssolution,thenwehaveG(^G(0)=0.ThereforeifwecanshowthatGisnon-negativefor2n(M),then^willliewithinthesphereofradiusMrn.Werequireq(p+s)logp=opn.kSk2Fk0Sk2F=trTTS+STS)trT000S+STS)=trTT00)2tr0)S)=tr0+T0+T00)2trTS)=trT2trT(S0))86Next,weboundtheterminvolvingSintheaboveexpression,wehavejtr0S))jXi6=jjij0ijSij)j+Xi=1jii0iiSii)jmaxi6=j(j0ijSijj)kk1+ppmaxi=1(j0iiSiij)sXi=12iiC0(1+˝)maxi0ii)ˆslogpnkk1+splogpnk+k2˙C1ˆslogpnkk1+splogpnk+k2˙holdswithhighprobability,byaresult(Lemma1fromRavikumaretal.[2011])whereC1=C0(1+˝)maxi0ii);C0>0and+ismatrixwithallo˙-diagonalelementssetequaltozero.Nextweobtainupperboundonthetermsinvolving.wehave,tr(D2)(tr(D))2=ptr(D20)(tr(D))2=p=tr2)tr20)(tr2=p+(tr0))2=p(i)tr2)20))tr0+2tr0)2=tr2+2tr20)tr2+C1pskkF(ii)tr2(tr0))2=(tr0+2(tr0))2(tr2+2tr0)trpkk2F+2kpppk+kF:Thereforethetermcanbeboundedby2kk2F+(C1ps+2ppk)kkF.WeboundtheterminvolvingasinsimilartotheproofofTheorem6.3.1.Forqlogpn,therestoftheprooffollowsverysimilartoTheorem6.3.1.87APPENDIXBJCPFORPRECISIONMATRIXESTIMATIONProof.[Lemma11.3.1]LetMbeasolutionof(11.3.6)whichexitsbecause(11.3.6)isaconvexoptimizationproblem.Thenthefollowingsub-gradientoptimalityconditionholds[VandenbergheandBoyd[2004]]02MM+kMk1:(B.0.1)where((@kMk1))ij=@jmijjandgivenby@jmijj=8>>>>>><>>>>>>:+1ifmij>01ifmij<0i=1;2;:::m;j=1;2::::n:2[1;1]ifmi;j=0Notethat(B.0.1)issatis˝edifandonlyifjmijjandthereforeoptimialsolutionisgivenbymij=sign(mij)(mij)+.Thiscompletestheproof.Proof.[Lemma11.3.2]LetLbeasolutionto(11.3.7).Thenthefollowingsub-gradientoptimalityconditionholds:02LM+˝@kLk(B.0.2)LetW=U˝VT.WeshallshowthatthischoiceofWsatis˝estheaboveoptimalitycondi-tion.Thesub-di˙erential@kWkofkWkisgivenbyBach[2008]as@kWk=(UVT+HsuchthatH2Rmn;kHk21;UTH=0andHV=0):ThereforeWM+˝@kWk=U˝VTUVT+˝(UVT+H):88MultiplyingbothsidesbyUUT,andnotingthatUUT=IweobtainWM+@kWk=UUT(U˝VTUVT+˝(UVT+H))=U˝VTU˝I)VT+UUTH=0:Therefore,W=U˝VTisasolutionto(11.3.7),thiscompletestheproof.Proof.[Lemma11.3.2]ForsuchchoiceofLnwehaveF(Wn)=f(Wn)+kWnk1+kWnk=QLn(Wn;Wn1)+kWnk1+˝kWnk=f(Wn1)+Ln2kWnWn1k2+<WnWn1;5f(Wn1)>+kWnk1+˝kWnkAlsowehave,F(W)=f(W)+kWk1+kWkf(W)f(Wn1)+<WWn1;5f(Wn1)>kWk1kWnk1+<WWn;5kWnk1>kWkkWnk+<WWn;5kWnk>Weget,F(W)F(Wn)Ln2kWnWn1k2+<WWn;5f(Wn1)+5kWnk1+˝5kWnk>(B.0.3)NotethatWnissolutionof5f(Wn1)+Ln(WnWn1)+˝5kWnk=0(using11.3.1)Therefore(11.6.2)becomesF(W)F(Wn)Ln2 kWn1Wnk2+2<Wn1Wn;WnW>2Ln<WWn;5kWnk1>!WeknowthatforanythreematricesA;B;CkBAk2+2<BA;AC>=kBCk2kACk289Usingthis,weobtainF(Wn)F(W)Ln2 kWn1Wk2kWnWk22Ln<WWn;5kWnk1>!:Usinglemma(11.3.2),wegetF(Wn)F(W)Ln4 kWn1Wk2kWnWk2!+9c22Ln2<WWn;5kWnk1>:Thiscompletestheproof.90APPENDIXCJPENFORPRECISIONMATRIXESTIMATIONProof.[Theorem12.4.1]Toboundthecrossproductterminvolvingand^R1K,wehave,jtr((R10^R1Kj=jtr(R10(^RKR0)^R1Kj˙1(R10)jtr((^RKR0)^R1Kjk˙1(^R1K)jtr((^RKR0jkk1jtr((^RKR0j:where˙min(^RK)(1=k1)>0,isapositivelowerboundontheeigenvaluesofJPENestimate^RKofthecorrelationmatrixR0.SuchaconstantexistsbyLemma6.3.2.RestoftheproofcloselyfollowsthatofTheorem6.3.1.Proof.[Theorem12.4.2]Weboundthetermtr((^S0similartothatinproofofTheorem12.4.1.RestoftheproofcloselyfollowstothatTheorem12.4.1.91BIBLIOGRAPHY92BIBLIOGRAPHYI.JohnstoneandY.Lu.Sparseprincipalcomponentsanalysis.UnpublishedManuscript,2004.H.Zou,HastieT.,andTibshiraniR.Sparseprincipalcomponentsanalysis.JournalofComputationalandGraphicalStatistics,2006.K.Mardia,KentJ.,andBibbyJ.MultivariateAnalysis.,volume1.AcademicPress,NewYork,NY,1979.M.YuanandY.Lin.Modelselectionandestimationinthegaussiangraphicalmodel.Biometrika,2007.M.Wainwright,RavikumarP.,andLa˙ertyJ.High-dimensionalgraphicalmodelselectionusingl1-regularizedlogisticregression.ProceedingsofAdvancesinNeuralInformationProcessingSystems.,2006.M.Wainright.Sharpthresholdsforhigh-dimensionalandnoisysparsityrecoveryusingl1-constrainedquadraticprogrammming(lasso).IEEETransactionsonInformationTheoryarchive,55,2009.M.Yuan.Sparseinversecovariancematrixestimationvialinearprogramming.JournalofMachineLearningResearch,2009.MeinshausenandP.Bühlmann.Highdimensionalgraphsandvariableselectionwiththelasso.AnnalsofStatistics,2006.D.Yatsenko,K.Josic,rA.S.Ecke,E.Froudarakis,R.J.Cotton,andA.STolias.Improvedestimationandinterpretationofcorrelationsinneuralcircuits.PLoSComput.Biol,11,2015.V.MarcenkoandL.Pastur.Distributionsofeigenvaluesofsomesetsofrandommatrices.Math.USSR-Sb,1967.S.Geman.Alimittheoremforthenormofrandommatrices.TheAnnalsofStatistics,8(2):1980.D.L.Donoho,M.Gavish,andI.M.Johnstone.Optimalshrinkageofeigenvaluesinthespikedcovariancemodel.http://arxiv.org/pdf/1311.0851.pdf,2015.P.BickelandE.Levina.Regulatizedestimationoflargecovariancematrices.AnnalsofStatistics,2008a.P.BickelandE.Levina.Covarianceregularizationbythresholding.TheAnnalsofStatistics,2008b.J.BienandR.Tibshirani.Sparseestimationofacovariancematrix.Biometrica,2011.93A.Rothman.Positivede˝niteestimatorsoflargecovariancematrices.Biometrica,740,2012.L.Xue,MaS.,andZouHui.Positive-de˝nitel1-penalizedestimationoflargecovariancematrices.JournalofAmericanStatisticalAssociation,2012.J.Dahl,L.Vandenberghe,andV.Roychowdhury.Covarianceselectionfornon-chordalgraphsviachordalembedding.OptimizationMethodsandSoftware,2008.A.Dempster.Covarianceselection.Biometrika,1972.T.Cai,C.Zhang,andH.Zhou.Optimalratesofconvergenceforcovariancematrixestima-tion.TheAnnalsofStatistics,2010.N.Karoui.Operatornormconsistentestimationoflargedimensionalsparsecovariancematrices.TheAnnalsofStatistics,2008.A.Rothman,BickelP.J.,LevinaE.,andZhuJ.Sparsepermutationinvariantcovarianceestimation.ElectronicJournalofStatistics,2008.R.Tibshiran.Regressionshrinkageandselectionviathelasso.JournaloftheRoyalStatisticalSociety,SeriesB(StatisticalMethodology),pages1996.S.Chaudhuri,M.Drton,andT.S.Richardson.Estimationofacovariancematrixwithzeros.Biometrika,2007.A.J.Butte,P.Tamayo,D.Slonim,T.R.Golub,andI.S.Kohane.Discoveringfunctionalrelationshipsbetweenrnaexpressionandchemotherapeuticsusceptibilityusingrelevancenetworks.ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica,2000.C.LamandJ.Fan.Sparsistencyandratesofconvergenceinlargecovariancematricesestimation.AnnalsofStatistics,2009.O.LedoitandM.Wolf.Awell-conditionedestimatorforlarge-dimensionalcovariancema-trices.JournalofMultivariateAnalysis,2004.C.Stein.Estimationofacovariancematrix.Rietzlecture,39thAnnualMeetingIMS.Atlanta,Georgia,1975.C.Stein.Lecturesonthetheoryofestimationofmanyparameters.JournalofMathematicalSciences.,1986.O.LedoitandM.Wolf.Optimalestimationofalarge-dimensionalcovariancematrixunderstein'sloss.http://papers.ssrn.com/,2014.Y.SheenaandA.Gupta.Estimationofthemultivariatenormalcovariancematrixundersomerestrictions.StatisticsandDecisions,2003.94J.H.Won,J.Lim,S.J.Kim,andB.Rajaratnam.Condition-numberregularizedcovarianceestimation.JournaloftheRoyalStatisticalSocietyB,75,2012.L.R.Ha˙.Empiricalbayesestimationofthemultivariatenormalcovariancematrix.AnnalsofStatistics,1980.S.LinandM.D.Perlman.Amontecarlocomparisonoffourestimatorsofacovariancematrix.MultivariateAnalysis,1985.D.DeyandC.Srinivasan.Estimationofacovariancematrixunderstein'sloss.AnnalsofStatistics,1985.B.Rajaratnam,D.Vincenzi,andB.Naul.Atheoreticalstudyofstein'scovarianceestimator.Technicalreport,DepartmentofStatistics,StanfordUniversity,2014.A.Maurya.Asparseandwell-conditionedestimationofcovarianceandinversecovariancematricesusingajointpenalty.JournalofMachineLearningResearch,15-345,2016.T.Cai,Z.Ren,andH.Zhou.Estimatingstructuredhigh-dimensionalcovarianceandpre-cisionmatrices:Optimalratesandadaptiveestimation.ElectronicJournalofStatistics,2015.D.P.Bertsekas.Incrementalgradient,subgradient,andproximalmethodsforconvexopti-mization,asurvey.LabaratoryforInformationandDecisionSystemsReportLIDS-P-2848.MIT,2010.L.VandenbergheandS.Boyd.Convexoptimization.CambridgeUniversityPress,2004.F.Bach,R.Jenatton,J.Mairal,andG.Obozinski.Convexoptimizationwithsparsity-inducingnorms.InOptimizationforMachineLearning,MITpress,2011.A.BeckandM.Teboulle.Afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems.SIAMJournalofImagingScience,2009.J.Friedman,HastieT.,andTibshiraniR.Sparseinversecovarianceestimationwiththegraphicallasso.Biostatistics.,2008.C.R.Rao.Generalizedinverseofamatrixanditsapplications.Proc.SixthBerkeleySymp.onMath.Statist.andProb.,Univ.ofCalif.Press,1972.L.Vandenberghe,S.Boyd,andS.-P.Wu.Determinantmaximizationwithlinearmatrixinequalityconstraints.SIAMJournalonMatrixAnalysisandApplications,1998.O.Banerjee,E.L.Ghaoui,andd'AspremontA.Modelselectionthroughsparsemaximumlikelihoodestimationformultivariategaussianorbinarydata.JournalofMachineLearn-ingResearch,2008.S.Zhou,P.Rutimann,andBuhlmannP.XuM.High-dimensionalcovarianceestimationbasedongaussiangraphicalmodels.JournalofMachineLearningResearch,2011.95M.Pourahmadi.Choleskydecompositionsandestimationofacovariancematrix:orthogo-nalityofvariance-correlationparameters.Biometrika94,2007.M.Pourahmadi.Modelingcovariancematrices:Theglmandregularizationperspectives.StatisticalScience,2011.T.Cai,W.Liu,andX.Luo.,aconstrained`1minimizationapproachtosparseprecisionmatrixestimation.JournalofAmericanStatisticalAssociation,2011.R.TomiokaandK.Aihara.Classifyingmatriceswithaspectralregularization.Proc.24thInt.Conf.MachineLearning,pages2007.F.Bach.Consistencyoftracenormminimization.JournalofMachineLearningResearch,2008.A.Argyriou,T.EvgeniouT.,andM.Pontil.Convexmulti-taskfeaturelearning.MachineLearning,SpecialIssueonInductiveTransferLearning,2008.M.Fazel.Matrixrankminimizationwithapplications.,phdthesis.Elec.Eng.Dept,StanfordUniversity,2002.Y.Nesterov.Smoothminimizationofnon-smoothfunctions.Math.Program.,pages12005.R.T.Rockafellar.Monotoneoperatorsandproximalpointalgorithm.SIAMJournalofControlandOptimization,14,1976.H.Liu,K.Roede,andL.Wasserman.Stabilityapproachtoregularizationselection(stars)forhighdimensionalgraphicalmodels.InProceedingsoftheTwenty-ThirdAnnualConferenceonNeuralInformationProcessingSystems(NIPS),2010.A.Maurya.Ajointconvexpenaltyforinversecovariancematrixestimation.ComputationalStatisticsandDataAnalysis,2014.U.Alon,BarkaiN.,NottermanD.,GishK.,YbarraS.,MackD.,andLevineA.Broadpatternsofgeneexpressionrevealedbyclusteringanalysisoftumorandnormalcolontissuesprobedbyoligonucleotidearrays.ProceedingofNationalAcademyofScienceUSA,1999.L.Wang,J.Zhu,andH.Zou.Hybridhuberizedsupportvectormachinesformicroarrayclassi˝cation.Proceedingsofthe24thInternationalConferenceonMachineLearning.,pages2007.P.BickelandE.Levina.Sometheoryfor˝sher'slineardiscriminantfunction,ebayandsomealternativeswhentherearemanymorevariablesthanobservations.Bernoulli,2004.P.Ravikumar,WainwrightM.andRaskuttiG.,andYuB.High-dimensionalcovarianceestimationbyminimizingl1-penalizedlog-determinantdivergence.ElectronicJournalofStatistics,2011.96