# Computing a sparse solution of a set of linear inequalities

Jump to: Source code Text output Plots Library index

    % Section 6.2, Boyd & Vandenberghe "Convex Optimization"
    % "Just relax: Convex programming methods for subset selection
    %  and sparse approximation" by J. A. Tropp
    % Written for CVX by Almir Mutapcic - 02/28/06
    %
    % We consider a set of linear inequalities A*x <= b which are
    % feasible. We apply two heuristics to find a sparse point x that
    % satisfies these inequalities.
    %
    % The (standard) l1-norm heuristic for finding a sparse solution is:
    %
    %   minimize   ||x||_1
    %       s.t.   Ax <= b
    %
    % The log-based heuristic is an iterative method for finding
    % a sparse solution, by finding a local optimal point for the problem:
    %
    %   minimize   sum(log( delta + |x_i| ))
    %       s.t.   Ax <= b
    %
    % where delta is a small threshold value (determines what is close to zero).
    % We cannot solve this problem since it is a minimization of a concave
    % function and thus it is not a convex problem. However, we can apply
    % a heuristic in which we linearize the objective, solve, and re-iterate.
    % This becomes a weighted l1-norm heuristic:
    %
    %   minimize sum( W_i * |x_i| )
    %       s.t. Ax <= b
    %
    % which in each iteration re-adjusts the weights W_i based on the rule:
    %
    %   W_i = 1/(delta + |x_i|), where delta is a small threshold value
    %
    % This algorithm is described in papers:
    % "An Affine Scaling Methodology for Best Basis Selection"
    %  by B. D. Rao and K. Kreutz-Delgado
    % "Portfolio optimization with linear and ï¬?xed transaction costs"
    %  by M. S. Lobo, M. Fazel, and S. Boyd

    % fix random number generator so we can repeat the experiment
    seed = 0;
    randn('state',seed);
    rand('state',seed);

    % the threshold value below which we consider an element to be zero
    delta = 1e-8;

    % problem dimensions (m inequalities in n-dimensional space)
    m = 100;
    n = 50;

    % construct a feasible set of inequalities
    % (this system is feasible for the x0 point)
    A  = randn(m,n);
    x0 = randn(n,1);
    b  = A*x0 + rand(m,1);

    % l1-norm heuristic for finding a sparse solution
    fprintf(1, 'Finding a sparse feasible point using l1-norm heuristic ...')
    cvx_begin
      variable x_l1(n)
      minimize( norm( x_l1, 1 ) )
      subject to
        A*x_l1 <= b;
    cvx_end

    % number of nonzero elements in the solution (its cardinality or diversity)
    nnz = length(find( abs(x_l1) > delta ));
    fprintf(1,['\nFound a feasible x in R^%d that has %d nonzeros ' ...
               'using the l1-norm heuristic.\n'],n,nnz);

    % iterative log heuristic
    NUM_RUNS = 15;
    nnzs = [];
    W = ones(n,1); % initial weights

    disp([char(10) 'Log-based heuristic:']);
    for k = 1:NUM_RUNS
      cvx_begin quiet
        variable x_log(n)
        minimize( sum( W.*abs(x_log) ) )
        subject to
          A*x_log <= b;
      cvx_end

      % display new number of nonzeros in the solution vector
      nnz = length(find( abs(x_log) > delta ));
      nnzs = [nnzs nnz];
      fprintf(1,'   found a solution with %d nonzeros...\n', nnz);

      % adjust the weights and re-iterate
      W = 1./(delta + abs(x_log));
    end

    % number of nonzero elements in the solution (its cardinality or diversity)
    nnz = length(find( abs(x_log) > delta ));
    fprintf(1,['\nFound a feasible x in R^%d that has %d nonzeros ' ...
               'using the log heuristic.\n'],n,nnz);

    % plot number of nonzeros versus iteration
    plot(1:NUM_RUNS, nnzs, [1 NUM_RUNS],[nnzs(1) nnzs(1)],'--');
    axis([1 NUM_RUNS 0 n])
    xlabel('iteration'), ylabel('number of nonzeros (cardinality)');
    legend('log heuristic','l1-norm heuristic','Location','SouthEast')

    Finding a sparse feasible point using l1-norm heuristic ...
    Calling sedumi: 200 variables, 100 equality constraints
       For improved efficiency, sedumi is solving the dual problem.
    ------------------------------------------------------------
    SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
    Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
    eqs m = 100, order n = 201, dim = 201, blocks = 1
    nnz(A) = 5200 + 0, nnz(ADA) = 2650, nnz(L) = 1375
     it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
      0 :            8.50E-01 0.000
      1 :   3.98E+01 3.79E-01 0.000 0.4456 0.9000 0.9000   6.27  1  1  6.6E+00
      2 :  -5.85E+00 1.47E-01 0.000 0.3882 0.9000 0.9000   1.78  1  1  1.9E+00
      3 :  -2.47E+01 6.78E-02 0.000 0.4609 0.9000 0.9000   1.49  1  1  7.5E-01
      4 :  -3.26E+01 2.55E-02 0.000 0.3769 0.9000 0.9000   1.32  1  1  2.5E-01
      5 :  -3.46E+01 1.01E-02 0.000 0.3936 0.9000 0.9000   1.13  1  1  9.5E-02
      6 :  -3.55E+01 3.53E-03 0.000 0.3513 0.9000 0.9000   1.05  1  1  3.3E-02
      7 :  -3.58E+01 1.30E-03 0.000 0.3695 0.9000 0.9000   1.02  1  1  1.2E-02
      8 :  -3.59E+01 2.58E-05 0.000 0.0198 0.9258 0.9000   1.01  1  1  1.9E-03
      9 :  -3.59E+01 6.19E-06 0.000 0.2397 0.9000 0.0000   1.00  1  1  8.5E-04
     10 :  -3.59E+01 1.39E-06 0.000 0.2251 0.9000 0.8569   1.00  1  1  2.0E-04
     11 :  -3.59E+01 1.50E-07 0.000 0.1077 0.9138 0.9000   1.00  1  1  2.9E-05
     12 :  -3.59E+01 2.89E-08 0.000 0.1924 0.9180 0.9000   1.00  1  1  6.1E-06
     13 :  -3.59E+01 3.29E-09 0.000 0.1139 0.9000 0.9112   1.00  1  1  6.8E-07
     14 :  -3.59E+01 4.33E-12 0.000 0.0013 0.9990 0.9990   1.00  1  1
    iter seconds digits       c*x               b*y
     14      0.1   Inf -3.5940776942e+01 -3.5940776942e+01
    |Ax-b| =   2.7e-15, [Ay-c]_+ =   8.9E-15, |x|=  7.1e+00, |y|=  1.0e+01

    Detailed timing (sec)
       Pre          IPM          Post
    1.000E-02    6.000E-02    0.000E+00
    Max-norms: ||b||=1, ||c|| = 2.494989e+01,
    Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.58293.
    ------------------------------------------------------------
    Status: Solved
    Optimal value (cvx_optval): +35.9408

    Found a feasible x in R^50 that has 44 nonzeros using the l1-norm heuristic.

    Log-based heuristic:
       found a solution with 44 nonzeros...
       found a solution with 44 nonzeros...
       found a solution with 43 nonzeros...
       found a solution with 43 nonzeros...
       found a solution with 41 nonzeros...
       found a solution with 40 nonzeros...
       found a solution with 39 nonzeros...
       found a solution with 39 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...
       found a solution with 37 nonzeros...

    Found a feasible x in R^50 that has 37 nonzeros using the log heuristic.