sparse_covariance_est_tradeoff.m
de096106547amaster

sparse_covariance_est_tradeoff.m
View Options

			% Sparse covariance estimation for Gaussian variables
			% Joëlle Skaf - 04/24/08
			% (a figure is generated)
			%
			% Suppose y \in\reals^n is a Gaussian random variable with zero mean and
			% covariance matrix R = \Expect(yy^T), with sparse inverse S = R^{-1}
			% (S_ij = 0 means that y_i and y_j are conditionally independent).
			% We want to estimate the covariance matrix R based on N independent
			% samples y1,...,yN drawn from the distribution, and using prior knowledge
			% that S is sparse
			% A good heuristic for estimating R is to solve the problem
			% maximize logdet(S) - tr(SY) - lambda*sum(sum(abs(S)))
			% subject to S >= 0
			% where Y is the sample covariance of y1,...,yN, and lambda is a sparsity
			% parameter to be chosen or tuned.
			% A figure showing the sparsity (number of nonzeros) of S versus lambda
			% is generated.

			% Input data
			randn('state',0);
			n = 10;
			N = 100;
			Strue = sprandsym(n,0.5,0.01,1);
			nnz_true = sum(Strue(:)>1e-4);
			R = inv(full(Strue));
			y_sample = sqrtm(R)*randn(n,N);
			Y = cov(y_sample');
			Nlambda = 20;
			lambda = logspace(-2, 3, Nlambda);
			nnz = zeros(1,Nlambda);

			for i=1:Nlambda
			disp(['i = ' num2str(i) ', lambda(i) = ' num2str(lambda(i))]);
			% Maximum likelihood estimate of R^{-1}
			cvx_begin sdp quiet
			variable S(n,n) symmetric
			maximize log_det(S) - trace(SY) - lambda(i)sum(sum(abs(S)))
			S >= 0
			cvx_end
			nnz(i) = sum(S(:)>1e-4);
			end

			figure;
			semilogx(lambda, nnz);
			hold on;
			semilogx(lambda, nnz_true*ones(1,Nlambda),'r');
			xlabel('\lambda');
			legend('nonzeros in S', 'nonzeros in R^{-1}');