Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F102836868
sparse_solution.html
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Feb 24, 16:56
Size
10 KB
Mime Type
text/html
Expires
Wed, Feb 26, 16:56 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
24438007
Attached To
R1252 EMPoWER
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:
<a href="#source">Source code</a>
<a href="#output">Text output</a>
<a href="#plots">Plots</a>
<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 & 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 <= 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 <= 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 <= 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 <= 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 �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 <= b;
cvx_end
<span class="comment">% number of nonzero elements in the solution (its cardinality or diversity)</span>
nnz = length(find( abs(x_l1) > 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 <= b;
cvx_end
<span class="comment">% display new number of nonzeros in the solution vector</span>
nnz = length(find( abs(x_log) > 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) > 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) < 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>
Event Timeline
Log In to Comment