# Sparse covariance estimation for Gaussian variables

Jump to: Source code Text output Plots Library index

    % 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.
       sedumi 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 | 2.119e+00  3.337e-01  1.615e-08 | Solved
      1/  1 | 5.024e-03  1.815e-06  1.883e-08 | Solved
      1/  1 | 3.580e-04  2.791e-08  1.887e-08 | Solved
      0/  1 | 2.558e-05  1.891e-08  1.886e-08 | Solved
    -----------------------------------------------------------------
    Status: Solved
    Optimal value (cvx_optval): -31.2401