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)
%           subject to  sum(sum(abs(S))) <= alpha
%                       S >= 0
% where Y is the sample covariance of y1,...,yN, and alpha is a sparsity
% parameter to be chosen or tuned.

% Input data
randn('state',0);
n = 10;
N = 100;
Strue = sprandsym(n,0.5,0.01,1);
R = inv(full(Strue));
y_sample = sqrtm(R)*randn(n,N);
Y = cov(y_sample');
alpha = 50;

% Computing sparse estimate of R^{-1}
cvx_begin sdp
    variable S(n,n) symmetric
    maximize log_det(S) - trace(S*Y)
    sum(sum(abs(S))) <= alpha
    S >= 0
cvx_end
R_hat = inv(S);

S(find(S<1e-4)) = 0;
figure;
subplot(121);
spy(Strue);
title('Inverse of true covariance matrix')
subplot(122);
spy(S)
title('Inverse of estimated covariance matrix')
 
Successive approximation method to be employed.
   For improved efficiency, sedumi is solving the dual problem.
   sedumi will be called several times to refine the solution.
   Original size: 501 variables, 221 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
          Errors   
Act Centering    Conic    Status
-----------------------------------
1   2.734e+00  6.978e-01  Solved
1   7.454e-02  4.548e-04  Solved
1   9.197e-03  6.694e-06  Solved
1   1.106e-03  9.935e-09  Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.24