Homec4science
Diffusion EMPoWER (master) Browse

R1252/simulation-code/cvx/examples/cvxbook/Ch06_approx_fitting/html/basispursuit.html
de096106547amaster

Last Change
Raw File

basispursuit.html
View Options

<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Figures 6.21-6.23: Basis pursuit using Gabor functions</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/basispursuit.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Figures 6.21-6.23: Basis pursuit using Gabor functions</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.5.4</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Argyris Zymnis - 11/27/2005</span>
<span class="comment">%</span>
<span class="comment">% Here we find a sparse basis for a signal y out of</span>
<span class="comment">% a set of Gabor functions. We do this by solving</span>
<span class="comment">% minimize ||A*x-y||_2 + ||x||_1</span>
<span class="comment">%</span>
<span class="comment">% where the columns of A are sampled Gabor functions.</span>
<span class="comment">% We then fix the sparsity pattern obtained and solve</span>
<span class="comment">% minimize ||A*x-y||_2</span>
<span class="comment">%</span>
<span class="comment">% NOTE: The file takes a while to run</span>
clear
<span class="comment">% Problem parameters</span>
sigma = 0.05; <span class="comment">% Size of Gaussian function</span>
Tinv = 500; <span class="comment">% Inverse of sample time</span>
Thr = 0.001; <span class="comment">% Basis signal threshold</span>
kmax = 30; <span class="comment">% Number of signals are 2*kmax+1</span>
w0 = 5; <span class="comment">% Base frequency (w0 * kmax should be 150 for good results)</span>
<span class="comment">% Build sine/cosine basis</span>
fprintf(1,<span class="string">'Building dictionary matrix...'</span>);
<span class="comment">% Gaussian kernels</span>
TK = (Tinv+1)*(2*kmax+1);
t = (0:Tinv)'/Tinv;
A = exp(-t.^2/(sigma^2));
ns = nnz(A&gt;=Thr)-1;
A = A([ns+1:-1:1,2:ns+1],:);
ii = (0:2*ns)';
jj = ones(2*ns+1,1)*(1:Tinv+1);
oT = ones(1,Tinv+1);
A = sparse(ii(:,oT)+jj,jj,A(:,oT));
A = A(ns+1:ns+Tinv+1,:);
<span class="comment">% Sine/Cosine basis</span>
k = [ 0, reshape( [ 1 : kmax ; 1 : kmax ], 1, 2 * kmax ) ];
p = zeros(1,2*kmax+1); p(3:2:end) = -pi/2;
SC = cos(w0*t*k+ones(Tinv+1,1)*p);
<span class="comment">% Multiply</span>
ii = 1:numel(SC);
jj = rem(ii-1,Tinv+1)+1;
A = sparse(ii,jj,SC(:)) * A;
A = reshape(A,Tinv+1,(Tinv+1)*(2*kmax+1));
fprintf(1,<span class="string">'done.\n'</span>);
<span class="comment">% Construct example signal</span>
a = 0.5*sin(t*11)+1;
theta = sin(5*t)*30;
b = a.*sin(theta);
<span class="comment">% Solve the Basis Pursuit problem</span>
disp(<span class="string">'Solving Basis Pursuit problem...'</span>);
tic
cvx_begin
variable <span class="string">x(30561)</span>
minimize(sum_square(A*x-b)+norm(x,1))
cvx_end
disp(<span class="string">'done'</span>);
toc
<span class="comment">% Reoptimize problem over nonzero coefficients</span>
p = find(abs(x) &gt; 1e-5);
A2 = A(:,p);
x2 = A2 \ b;
<span class="comment">% Constants</span>
M = 61; <span class="comment">% Number of different Basis signals</span>
sk = 250; <span class="comment">% Index of s = 0.5</span>
<span class="comment">% Plot example basis functions;</span>
<span class="comment">%if (0) % to do this, re-run basispursuit.m to create A</span>
figure(1); clf;
subplot(3,1,1); plot(t,A(:,M*sk+1)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 1'</span>);
subplot(3,1,2); plot(t,A(:,M*sk+31)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 2'</span>);
subplot(3,1,3); plot(t,A(:,M*sk+61)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 3'</span>);
<span class="comment">%print -deps bp-dict_helv.eps</span>
<span class="comment">% Plot reconstructed signal</span>
figure(2); clf;
subplot(2,1,1);
plot(t,A2*x2,<span class="string">'--'</span>,t,b,<span class="string">'-'</span>); axis([0 1 -1.5 1.5]);
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y_{hat} and y'</span>);
title(<span class="string">'Original and Reconstructed signals'</span>)
subplot(2,1,2);
plot(t,A2*x2-b); axis([0 1 -0.06 0.06]);
title(<span class="string">'Reconstruction error'</span>)
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y - y_{hat}'</span>);
<span class="comment">%print -deps bp-approx_helv.eps</span>
<span class="comment">% Plot frequency plot</span>
figure(3); clf;
subplot(2,1,1);
plot(t,b); xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y'</span>); axis([0 1 -1.5 1.5]);
title(<span class="string">'Original Signal'</span>)
subplot(2,1,2);
plot(t,150*abs(cos(w0*t)),<span class="string">'--'</span>);
hold <span class="string">on</span>;
<span class="keyword">for</span> k = 1:length(t);
<span class="keyword">if</span>(abs(x((k-1)*M+1)) &gt; 1e-5), plot(t(k),0,<span class="string">'o'</span>); <span class="keyword">end</span>;
<span class="keyword">for</span> j = 2:2:kmax*2
<span class="keyword">if</span>((abs(x((k-1)*M+j)) &gt; 1e-5) | (abs(x((k-1)*M+j+1)) &gt; 1e-5)),
plot(t(k),w0*j/2,<span class="string">'o'</span>);
<span class="keyword">end</span>;
<span class="keyword">end</span>;
<span class="keyword">end</span>;
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'w'</span>);
title(<span class="string">'Instantaneous frequency'</span>)
hold <span class="string">off</span>;
</pre>
<a id="output"></a>
<pre class="codeoutput">
Building dictionary matrix...done.
Solving Basis Pursuit problem...
Calling SDPT3: 61625 variables, 502 equality constraints
------------------------------------------------------------
num. of constraints = 502
dim. of socp var = 61625, num. of socp blk = 30562
*******************************************************************
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+00|1.7e+02|7.6e+06| 4.324221e+04 0.000000e+00| 0:0:02| chol 1 1
1|1.000|0.995|3.6e-07|1.0e+00|8.6e+04| 4.270270e+04 -7.964108e+01| 0:0:04| chol 1 1
2|1.000|0.989|4.6e-07|2.1e-02|1.1e+04| 1.102337e+04 -7.582646e-01| 0:0:07| chol 1 1
3|0.970|0.989|2.1e-08|1.2e-03|3.4e+02| 3.397756e+02 3.356358e-01| 0:0:10| chol 1 1
4|0.484|0.886|1.2e-08|2.3e-04|2.1e+02| 2.206064e+02 6.421842e+00| 0:0:14| chol 1 1
5|0.772|0.915|2.8e-09|2.8e-05|1.0e+02| 1.106558e+02 1.003458e+01| 0:0:17| chol 1 1
6|0.507|0.828|1.4e-09|5.7e-06|6.6e+01| 7.735771e+01 1.171877e+01| 0:0:20| chol 1 1
7|0.514|0.678|6.7e-10|1.9e-06|4.0e+01| 5.172854e+01 1.210713e+01| 0:0:23| chol 1 1
8|0.362|0.850|4.3e-10|3.0e-07|3.0e+01| 4.235935e+01 1.241274e+01| 0:0:26| chol 1 1
9|0.447|0.604|2.4e-10|1.2e-07|2.0e+01| 3.222978e+01 1.253063e+01| 0:0:30| chol 1 1
10|0.527|0.932|1.1e-10|8.1e-09|1.2e+01| 2.482939e+01 1.268651e+01| 0:0:33| chol 1 1
11|0.624|0.687|4.2e-11|2.6e-09|6.6e+00| 1.932119e+01 1.274802e+01| 0:0:36| chol 1 1
12|0.599|0.980|1.7e-11|6.1e-11|3.7e+00| 1.649674e+01 1.281034e+01| 0:0:39| chol 1 1
13|0.543|0.647|7.7e-12|2.5e-11|2.2e+00| 1.507030e+01 1.282581e+01| 0:0:42| chol 1 1
14|0.541|0.964|3.5e-12|2.5e-12|1.3e+00| 1.413866e+01 1.283901e+01| 0:0:45| chol 1 1
15|0.728|0.656|1.1e-12|1.8e-12|5.3e-01| 1.337669e+01 1.284169e+01| 0:0:49| chol 1 1
16|0.458|0.608|6.3e-13|1.7e-12|3.4e-01| 1.318123e+01 1.284339e+01| 0:0:52| chol 1 1
17|0.590|0.955|5.3e-13|1.1e-12|1.9e-01| 1.303175e+01 1.284468e+01| 0:0:55| chol 1 1
18|0.961|0.987|2.5e-12|1.0e-12|4.6e-02| 1.289055e+01 1.284500e+01| 0:0:58| chol 1 1
19|0.761|0.883|4.2e-12|1.1e-12|1.6e-02| 1.286137e+01 1.284512e+01| 0:1:01| chol 1 1
20|0.882|0.884|9.5e-12|1.1e-12|2.8e-03| 1.284795e+01 1.284515e+01| 0:1:05| chol 1 1
21|0.943|0.968|4.8e-11|1.5e-12|2.8e-04| 1.284543e+01 1.284516e+01| 0:1:08| chol 2 2
22|0.628|0.970|3.8e-11|2.3e-12|1.5e-04| 1.284530e+01 1.284516e+01| 0:1:11| chol 1 2
23|0.631|1.000|2.2e-11|3.4e-12|7.8e-05| 1.284523e+01 1.284516e+01| 0:1:14| chol 1 2
24|0.623|1.000|1.1e-11|4.3e-12|4.2e-05| 1.284520e+01 1.284516e+01| 0:1:17| chol 1 2
25|0.619|1.000|6.8e-12|2.3e-12|2.3e-05| 1.284518e+01 1.284516e+01| 0:1:21| chol 2 2
26|0.602|1.000|5.0e-12|1.4e-12|1.3e-05| 1.284517e+01 1.284516e+01| 0:1:24| chol 2 2
27|0.591|1.000|3.4e-12|1.0e-12|7.0e-06| 1.284516e+01 1.284516e+01| 0:1:27| chol 2 1
28|0.599|1.000|1.1e-11|1.0e-12|3.9e-06| 1.284516e+01 1.284516e+01| 0:1:30| chol 2 2
29|0.606|1.000|5.3e-12|1.5e-12|2.2e-06| 1.284516e+01 1.284516e+01| 0:1:33| chol 2 2
30|0.612|1.000|3.0e-12|1.1e-12|1.2e-06| 1.284516e+01 1.284516e+01| 0:1:37| chol 2 2
31|0.618|1.000|2.4e-12|1.0e-12|6.4e-07| 1.284516e+01 1.284516e+01| 0:1:40| chol 2 2
32|0.622|1.000|2.1e-12|1.0e-12|3.4e-07| 1.284516e+01 1.284516e+01| 0:1:43|
stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
number of iterations = 32
primal objective value = 1.28451560e+01
dual objective value = 1.28451557e+01
gap := trace(XZ) = 3.44e-07
relative gap = 1.29e-08
actual relative gap = 1.29e-08
rel. primal infeas = 2.06e-12
rel. dual infeas = 1.00e-12
norm(X), norm(y), norm(Z) = 3.5e+00, 9.9e-01, 1.9e+02
norm(A), norm(b), norm(C) = 6.9e+02, 1.9e+01, 1.8e+02
Total CPU time (secs) = 102.93
CPU time per iteration = 3.22
termination code = 0
DIMACS: 1.6e-11 0.0e+00 8.8e-11 0.0e+00 1.3e-08 1.3e-08
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +12.8452
done
Elapsed time is 106.049858 seconds.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="basispursuit__01.png" alt=""> <img src="basispursuit__02.png" alt=""> <img src="basispursuit__03.png" alt="">
</div>
</div>
</body>
</html>
c4science · Help