% 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
rand('state',0);
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.
   SDPT3 will be called several times to refine the solution.
   Original size: 503 variables, 280 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 1.843e+00  3.337e-01  0.000e+00 | Solved
  1/  1 | 2.766e-01  6.363e-03  0.000e+00 | Solved
  1/  1 | 4.484e-03  1.630e-06  0.000e+00 | Solved
  0/  0 | 1.487e-05  0.000e+00  0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.2401