Graphical-Lasso

GraphicalLassoMatlabimplementationofthegraphicalLassomodelforestimatingsparseinversecovariancematrix(a.k.a.precisionorconcentrationmatrix)用于估计逆协方差矩阵的graphicalLasso模型的Matlab实现。minimizetr(Theta*S)-logdet(Theta)+ρ*||Theta||1overallpositive-definiteandsymmetricmatricesTheta.Sisanestimateofthecovariancematrix(usuallysamplecovariancematrix)andρisaregularizationparameter.Theta为正定对称矩阵，S是协方差矩阵的估计（通常为样本协方差矩阵），ρ是正则化参数I/O:Input:samplecovariancematrixS,penaltyparameterρ.Output:theestimatedprecisionmatrixandtheregularizedcovariancematrix.输入：样本协方差矩阵S和惩罚参数ρ输出：估计的精度矩阵和正则化协方差矩阵Example:Wesimulateanexamplewith50variablesand200observations.Left:trueprecisionmatrixpattern.Right:estimatedprecisionmatrixpatternbygraphicalLasso.我们用50个变量，200个观测值，仿真了一个案例，Upper：真正的协方差矩阵Down:通过graphicallasso估计的协方差矩阵GraphicalLasso代码%GraphicalLassofunction%Author:XiaohuiChen([e-mail])%Version:2012-Febfunction[ThetaW]=graphicalLasso(S,rho,maxIt,tol)%SolvethegraphicalLasso%minimize_{Theta0}tr(S*Theta)-logdet(Theta)+rho*||Theta||_1%Ref:Friedmanetal.(2007)Sparseinversecovarianceestimationwiththe%graphicallasso.Biostatistics.%Note:ThisfunctionneedstocallanalgorithmthatsolvestheLasso%problem.Here,wechoosetousetothefunction*lassoShooting*(shooting%algorithm)forthispurpose.However,anyLassoalgorithminthe%penelizedformwillwork.%%Input:%S--samplecovariancematrix%rho--regularizationparameter%maxIt--maximumnumberofiterations%tol--convergencetolerancelevel%%Output:%Theta--inversecovariancematrixestimate%W--regularizedcovariancematrixestimate,W=Theta^-1p=size(S,1);ifnargin4,tol=1e-6;endifnargin3,maxIt=1e2;end%InitializationW=S+rho*eye(p);%diagonalofWremainsunchangedW_old=W;i=0;%GraphicalLassoloopwhileimaxIt,i=i+1;forj=p:-1:1,jminus=setdiff(1:p,j);[VD]=eig(W(jminus,jminus));d=diag(D);X=V*diag(sqrt(d))*V';%W_11^(1/2)Y=V*diag(1./sqrt(d))*V'*S(jminus,j);%W_11^(-1/2)*s_12b=lassoShooting(X,Y,rho,maxIt,tol);W(jminus,j)=W(jminus,jminus)*b;W(j,jminus)=W(jminus,j)';end%Stopcriterionifnorm(W-W_old,1)tol,break;endW_old=W;endifi==maxIt,fprintf('%s\n','Maximumnumberofiterationreached,glassomaynotconverge.');endTheta=W^-1;%ShootingalgorithmforLasso(unstandardizedversion)functionb=lassoShooting(X,Y,lambda,maxIt,tol),ifnargin4,tol=1e-6;endifnargin3,maxIt=1e2;end%Initialization[n,p]=size(X);ifpn,b=zeros(p,1);%Fromthenullmodel,ifpnelseb=X\Y;%FromtheOLSestimate,ifp=nendb_old=b;i=0;%PrecomputeX'XandX'YXTX=X'*X;XTY=X'*Y;%ShootingloopwhileimaxIt,i=i+1;forj=1:p,jminus=setdiff(1:p,j);S0=XTX(j,jminus)*b(jminus)-XTY(j);%S0=X(:,j)'*(X(:,jminus)*b(jminus)-Y)ifS0lambda,b(j)=(lambda-S0)/norm(X(:,j),2)^2;elseifS0-lambda,b(j)=-(lambda+S0)/norm(X(:,j),2)^2;elseb(j)=0;endenddelta=norm(b-b_old,1);%Normchangeduringsuccessiveiterationsifdeltatol,break;endb_old=b;endifi==maxIt,fprintf('%s\n','Maximumnumberofiterationreached,shootingmaynotconverge.');endReference~xiaohuic/code/glasso/glasso.htm