% 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
        variable S(n,n) symmetric
        maximize log_det(S) - trace(S*Y) - 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}');
i = 1, lambda(i) = 0.01
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.125e+00  3.355e-01  1.255e-07 | Solved
  1/  1 | 5.487e-03  2.254e-06  9.933e-08 | Solved
  0/  1 | 3.910e-04  1.089e-07  9.895e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.3012
i = 2, lambda(i) = 0.01833
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.130e+00  3.370e-01  7.234e-08 | Solved
  1/  1 | 5.818e-03  2.430e-06  2.271e-08 | Solved
  1/  1 | 4.145e-04  3.492e-08  2.271e-08 | Solved
  0/  1 | 2.961e-05  2.277e-08  2.270e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.3504
i = 3, lambda(i) = 0.033598
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.138e+00  3.396e-01  4.734e-08 | Solved
  1/  1 | 6.455e-03  3.080e-06  1.118e-07 | Solved
  0/  1 | 4.597e-04  1.269e-07  1.122e-07 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.4368
i = 4, lambda(i) = 0.061585
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.151e+00  3.441e-01  4.352e-08 | Solved
  1/  1 | 7.595e-03  4.110e-06  6.682e-09 | Solved
  1/  1 | 5.405e-04  2.740e-08  6.681e-09 | Solved
  0/  1 | 3.862e-05  6.836e-09  6.725e-09 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.5845
i = 5, lambda(i) = 0.11288
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.174e+00  3.515e-01  6.697e-08 | Solved
  1/  1 | 9.467e-03  6.486e-06  1.061e-07 | Solved
  1/  1 | 6.733e-04  1.377e-07  1.058e-07 | Solved
  0/  1 | 4.811e-05  1.060e-07  1.059e-07 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.8272
i = 6, lambda(i) = 0.20691
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.209e+00  3.633e-01  6.798e-09 | Solved
  1/  1 | 1.248e-02  1.109e-05  1.180e-08 | Solved
  1/  1 | 8.860e-04  6.712e-08  1.195e-08 | Solved
  0/  1 | 6.331e-05  1.209e-08  1.176e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -32.2113
i = 7, lambda(i) = 0.37927
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.264e+00  3.818e-01  5.682e-08 | Solved
  1/  1 | 1.727e-02  2.127e-05  5.134e-08 | Solved
  1/  1 | 1.223e-03  1.616e-07  5.588e-08 | Solved
  0/  1 | 8.742e-05  5.606e-08  5.546e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -32.8011
i = 8, lambda(i) = 0.69519
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.344e+00  4.098e-01  1.613e-07 | Solved
  1/  1 | 2.464e-02  4.333e-05  1.018e-07 | Solved
  1/  1 | 1.740e-03  3.169e-07  1.025e-07 | Solved
  0/  1 | 1.244e-04  1.033e-07  1.021e-07 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -33.6657
i = 9, lambda(i) = 1.2743
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.454e+00  4.503e-01  3.628e-08 | Solved
  1/  1 | 3.573e-02  9.096e-05  2.945e-08 | Solved
  1/  1 | 2.510e-03  4.726e-07  3.051e-08 | Solved
  0/  1 | 1.795e-04  3.355e-08  3.081e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -34.876
i = 10, lambda(i) = 2.3357
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.601e+00  5.072e-01  2.229e-08 | Solved
  1/  1 | 5.185e-02  1.915e-04  1.271e-07 | Solved
  1/  1 | 3.612e-03  1.106e-06  1.789e-07 | Solved
  0/  1 | 2.550e-04  2.526e-08  2.071e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -36.4981
i = 11, lambda(i) = 4.2813
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 2.786e+00  5.839e-01  1.381e-08 | Solved
  1/  1 | 7.449e-02  3.953e-04  1.529e-08 | Solved
  1/  1 | 5.131e-03  1.861e-06  1.158e-08 | Solved
  0/  1 | 3.671e-04  1.401e-08  1.178e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -38.5618
i = 12, lambda(i) = 7.8476
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 3.011e+00  6.842e-01  7.251e-09 | Solved
  1/  1 | 1.053e-01  7.906e-04  3.331e-09 | Solved
  1/  1 | 7.155e-03  3.630e-06  5.011e-09 | Solved
  1/  1 | 5.127e-04  2.168e-08  5.055e-09 | Solved
  0/  1 | 3.661e-05  2.136e-09  5.052e-09 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -41.092
i = 13, lambda(i) = 14.3845
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 3.273e+00  8.124e-01  1.292e-10 | Solved
  1/  1 | 1.463e-01  1.528e-03  9.899e-08 | Solved
  1/  1 | 9.745e-03  6.735e-06  7.207e-08 | Solved
  0/  1 | 7.002e-04  3.429e-08  6.008e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -44.1061
i = 14, lambda(i) = 26.3665
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 3.584e+00  9.786e-01  1.254e-09 | Solved
  1/  1 | 2.014e-01  2.896e-03  1.615e-08 | Solved
  1/  1 | 1.300e-02  1.199e-05  1.282e-10 | Solved
  1/  1 | 9.410e-04  6.470e-08  7.387e-09 | Solved
  0/  0 | 0.000e+00  0.000e+00  0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -47.7294
i = 15, lambda(i) = 48.3293
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 3.940e+00  1.190e+00  2.178e-08 | Solved
  1/  1 | 2.737e-01  5.358e-03  4.829e-10 | Solved
  1/  1 | 1.705e-02  2.062e-05  1.092e-09 | Solved
  1/  1 | 1.227e-03  1.110e-07  1.033e-09 | Solved
  1/  1 | 8.758e-05  1.404e-09  1.043e-09 | Solved
  0/  1 | 6.256e-06  1.075e-09  1.038e-09 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -51.9812
i = 16, lambda(i) = 88.5867
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 4.335e+00  1.448e+00  9.696e-09 | Solved
  1/  1 | 3.644e-01  9.520e-03  3.968e-08 | Solved
  1/  1 | 2.159e-02  3.315e-05  3.598e-08 | Solved
  1/  1 | 1.560e-03  1.707e-07  4.037e-08 | Solved
  0/  1 | 1.112e-04  6.868e-10  3.981e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -56.7913
i = 17, lambda(i) = 162.3777
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 4.755e+00  1.752e+00  4.613e-09 | Solved
  1/  1 | 4.738e-01  1.613e-02  2.491e-08 | Solved
  1/  1 | 2.631e-02  4.924e-05  5.624e-08 | Solved
  0/  1 | 1.856e-03  2.248e-07  1.052e-06 | Inaccurate/Solved
  0/  1 | 9.017e-05  1.853e-11  1.060e-07 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -62.0441
i = 18, lambda(i) = 297.6351
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 5.190e+00  2.097e+00  1.453e-07 | Solved
  1/  1 | 6.002e-01  2.595e-02  1.851e-10 | Solved
  1/  1 | 3.072e-02  6.695e-05  2.359e-09 | Solved
  1/  1 | 2.226e-03  3.612e-07  2.405e-09 | Solved
  0/  1 | 1.589e-04  6.270e-10  2.401e-09 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -67.6118
i = 19, lambda(i) = 545.5595
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 5.629e+00  2.479e+00  1.728e-06 | Solved
  1/  1 | 7.412e-01  3.970e-02  1.001e-07 | Solved
  1/  1 | 3.437e-02  8.404e-05  1.021e-07 | Solved
  1/  1 | 2.497e-03  4.346e-07  4.581e-08 | Solved
  0/  1 | 1.646e-04  1.906e-09  1.348e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -73.3839
i = 20, lambda(i) = 1000
 
Successive approximation method to be employed.
   sedumi will be called several times to refine the solution.
   Original size: 502 variables, 279 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 6.068e+00  2.891e+00  8.466e-07 | Solved
  1/  1 | 8.947e-01  5.802e-02  6.623e-07 | Solved
  1/  1 | 3.671e-02  9.585e-05  1.067e-07 | Solved
  1/  1 | 2.682e-03  5.072e-07  8.768e-08 | Solved
  0/  1 | 2.031e-04  1.601e-09  9.269e-08 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -79.2801