Homec4science
Diffusion EMPoWER (master) Browse

R1252/simulation-code/cvx/examples/sparse_heuristics/html/sparse_solution.html
de096106547amaster

Last Change
Raw File

sparse_solution.html
View Options

<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Computing a sparse solution of a set of linear inequalities</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/sparse_heuristics/html/sparse_solution.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Computing a sparse solution of a set of linear inequalities</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.2, Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% "Just relax: Convex programming methods for subset selection</span>
<span class="comment">% and sparse approximation" by J. A. Tropp</span>
<span class="comment">% Written for CVX by Almir Mutapcic - 02/28/06</span>
<span class="comment">%</span>
<span class="comment">% We consider a set of linear inequalities A*x &lt;= b which are</span>
<span class="comment">% feasible. We apply two heuristics to find a sparse point x that</span>
<span class="comment">% satisfies these inequalities.</span>
<span class="comment">%</span>
<span class="comment">% The (standard) l1-norm heuristic for finding a sparse solution is:</span>
<span class="comment">%</span>
<span class="comment">% minimize ||x||_1</span>
<span class="comment">% s.t. Ax &lt;= b</span>
<span class="comment">%</span>
<span class="comment">% The log-based heuristic is an iterative method for finding</span>
<span class="comment">% a sparse solution, by finding a local optimal point for the problem:</span>
<span class="comment">%</span>
<span class="comment">% minimize sum(log( delta + |x_i| ))</span>
<span class="comment">% s.t. Ax &lt;= b</span>
<span class="comment">%</span>
<span class="comment">% where delta is a small threshold value (determines what is close to zero).</span>
<span class="comment">% We cannot solve this problem since it is a minimization of a concave</span>
<span class="comment">% function and thus it is not a convex problem. However, we can apply</span>
<span class="comment">% a heuristic in which we linearize the objective, solve, and re-iterate.</span>
<span class="comment">% This becomes a weighted l1-norm heuristic:</span>
<span class="comment">%</span>
<span class="comment">% minimize sum( W_i * |x_i| )</span>
<span class="comment">% s.t. Ax &lt;= b</span>
<span class="comment">%</span>
<span class="comment">% which in each iteration re-adjusts the weights W_i based on the rule:</span>
<span class="comment">%</span>
<span class="comment">% W_i = 1/(delta + |x_i|), where delta is a small threshold value</span>
<span class="comment">%</span>
<span class="comment">% This algorithm is described in papers:</span>
<span class="comment">% "An Affine Scaling Methodology for Best Basis Selection"</span>
<span class="comment">% by B. D. Rao and K. Kreutz-Delgado</span>
<span class="comment">% "Portfolio optimization with linear and &iuml;&not;?xed transaction costs"</span>
<span class="comment">% by M. S. Lobo, M. Fazel, and S. Boyd</span>
<span class="comment">% fix random number generator so we can repeat the experiment</span>
seed = 0;
randn(<span class="string">'state'</span>,seed);
rand(<span class="string">'state'</span>,seed);
<span class="comment">% the threshold value below which we consider an element to be zero</span>
delta = 1e-8;
<span class="comment">% problem dimensions (m inequalities in n-dimensional space)</span>
m = 100;
n = 50;
<span class="comment">% construct a feasible set of inequalities</span>
<span class="comment">% (this system is feasible for the x0 point)</span>
A = randn(m,n);
x0 = randn(n,1);
b = A*x0 + rand(m,1);
<span class="comment">% l1-norm heuristic for finding a sparse solution</span>
fprintf(1, <span class="string">'Finding a sparse feasible point using l1-norm heuristic ...'</span>)
cvx_begin
variable <span class="string">x_l1(n)</span>
minimize( norm( x_l1, 1 ) )
subject <span class="string">to</span>
A*x_l1 &lt;= b;
cvx_end
<span class="comment">% number of nonzero elements in the solution (its cardinality or diversity)</span>
nnz = length(find( abs(x_l1) &gt; delta ));
fprintf(1,[<span class="string">'\nFound a feasible x in R^%d that has %d nonzeros '</span> <span class="keyword">...</span>
<span class="string">'using the l1-norm heuristic.\n'</span>],n,nnz);
<span class="comment">% iterative log heuristic</span>
NUM_RUNS = 15;
nnzs = [];
W = ones(n,1); <span class="comment">% initial weights</span>
disp([char(10) <span class="string">'Log-based heuristic:'</span>]);
<span class="keyword">for</span> k = 1:NUM_RUNS
cvx_begin <span class="string">quiet</span>
variable <span class="string">x_log(n)</span>
minimize( sum( W.*abs(x_log) ) )
subject <span class="string">to</span>
A*x_log &lt;= b;
cvx_end
<span class="comment">% display new number of nonzeros in the solution vector</span>
nnz = length(find( abs(x_log) &gt; delta ));
nnzs = [nnzs nnz];
fprintf(1,<span class="string">' found a solution with %d nonzeros...\n'</span>, nnz);
<span class="comment">% adjust the weights and re-iterate</span>
W = 1./(delta + abs(x_log));
<span class="keyword">end</span>
<span class="comment">% number of nonzero elements in the solution (its cardinality or diversity)</span>
nnz = length(find( abs(x_log) &gt; delta ));
fprintf(1,[<span class="string">'\nFound a feasible x in R^%d that has %d nonzeros '</span> <span class="keyword">...</span>
<span class="string">'using the log heuristic.\n'</span>],n,nnz);
<span class="comment">% plot number of nonzeros versus iteration</span>
plot(1:NUM_RUNS, nnzs, [1 NUM_RUNS],[nnzs(1) nnzs(1)],<span class="string">'--'</span>);
axis([1 NUM_RUNS 0 n])
xlabel(<span class="string">'iteration'</span>), ylabel(<span class="string">'number of nonzeros (cardinality)'</span>);
legend(<span class="string">'log heuristic'</span>,<span class="string">'l1-norm heuristic'</span>,<span class="string">'Location'</span>,<span class="string">'SouthEast'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
Finding a sparse feasible point using l1-norm heuristic ...
Calling SDPT3: 200 variables, 100 equality constraints
------------------------------------------------------------
num. of constraints = 100
dim. of socp var = 100, num. of socp blk = 50
dim. of linear var = 100
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
NT 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|1.5e+01|1.5e+01|1.3e+05| 2.306427e+02 0.000000e+00| 0:0:00| chol 1 1
1|0.607|0.203|5.7e+00|1.2e+01|9.6e+04| 2.372055e+03 -7.798589e+01| 0:0:00| chol 1 1
2|0.527|0.755|2.7e+00|3.0e+00|4.4e+04| 2.913698e+03 -4.649185e+02| 0:0:00| chol 1 1
3|0.783|1.000|5.9e-01|2.7e-02|1.1e+04| 2.265689e+03 -7.312838e+02| 0:0:00| chol 1 1
4|0.933|1.000|3.9e-02|8.1e-03|1.2e+03| 1.752965e+02 -6.105356e+02| 0:0:00| chol 1 1
5|1.000|0.960|7.9e-07|8.9e-03|2.8e+02| 1.050573e+02 -1.760406e+02| 0:0:00| chol 1 1
6|0.842|0.872|3.2e-07|1.2e-03|5.3e+01| 5.663446e+01 3.546144e+00| 0:0:00| chol 1 1
7|0.590|0.753|1.4e-07|3.1e-04|3.3e+01| 4.833162e+01 1.570252e+01| 0:0:00| chol 1 1
8|0.864|0.581|2.2e-08|1.3e-04|1.9e+01| 4.206510e+01 2.318715e+01| 0:0:00| chol 1 1
9|1.000|0.553|1.2e-09|5.8e-05|1.1e+01| 3.937001e+01 2.867040e+01| 0:0:00| chol 1 1
10|1.000|0.851|4.1e-14|8.6e-06|3.9e+00| 3.750190e+01 3.356489e+01| 0:0:00| chol 1 1
11|0.857|0.651|2.7e-14|3.0e-06|2.2e+00| 3.675586e+01 3.455323e+01| 0:0:00| chol 1 1
12|0.885|1.000|2.9e-14|8.2e-11|1.1e+00| 3.647209e+01 3.541597e+01| 0:0:00| chol 1 1
13|1.000|0.971|5.8e-14|1.1e-11|2.2e-01| 3.604244e+01 3.581861e+01| 0:0:00| chol 1 1
14|1.000|0.667|2.8e-12|5.3e-12|9.5e-02| 3.597437e+01 3.587946e+01| 0:0:00| chol 1 1
15|0.965|0.928|6.8e-13|1.5e-12|1.7e-02| 3.594769e+01 3.593091e+01| 0:0:00| chol 1 1
16|0.873|1.000|7.3e-13|1.0e-12|5.3e-03| 3.594352e+01 3.593821e+01| 0:0:00| chol 1 2
17|0.920|0.980|1.7e-12|1.0e-12|5.9e-04| 3.594114e+01 3.594055e+01| 0:0:00| chol 2 2
18|1.000|0.977|8.6e-12|1.0e-12|5.9e-05| 3.594080e+01 3.594074e+01| 0:0:00| chol 2 1
19|0.996|0.998|3.4e-11|1.5e-12|9.1e-07| 3.594078e+01 3.594078e+01| 0:0:00|
stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
number of iterations = 19
primal objective value = 3.59407774e+01
dual objective value = 3.59407765e+01
gap := trace(XZ) = 9.08e-07
relative gap = 1.25e-08
actual relative gap = 1.24e-08
rel. primal infeas = 3.40e-11
rel. dual infeas = 1.50e-12
norm(X), norm(y), norm(Z) = 1.7e+01, 1.6e+00, 9.9e+00
norm(A), norm(b), norm(C) = 7.3e+01, 8.9e+01, 8.1e+00
Total CPU time (secs) = 0.28
CPU time per iteration = 0.01
termination code = 0
DIMACS: 1.2e-10 0.0e+00 6.1e-12 0.0e+00 1.2e-08 1.2e-08
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +35.9408
Found a feasible x in R^50 that has 46 nonzeros using the l1-norm heuristic.
Log-based heuristic:
found a solution with 46 nonzeros...
found a solution with 37 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
found a solution with 35 nonzeros...
Found a feasible x in R^50 that has 35 nonzeros using the log heuristic.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="sparse_solution__01.png" alt="">
</div>
</div>
</body>
</html>
c4science · Help