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, 19 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