Candes.m

From Wikimization

(Difference between revisions)
Jump to: navigation, search
(wIdHeYggpGXTenP)
m (Reverted edits by 84.244.189.99 (Talk); changed back to last version by Ranjelin)
Line 1: Line 1:
-
jGxq35 <a href="http://zxfvgjjlvqej.com/">zxfvgjjlvqej</a>, [url=http://xerodzoxagxo.com/]xerodzoxagxo[/url], [link=http://zkasnuzddvmk.com/]zkasnuzddvmk[/link], http://uuspekpdbwgk.com/
+
This Matlab demonstration of ''compressive sampling'' (<tt>a.k.a.</tt> ''compressed sensing'') by [http://www.acm.caltech.edu/~emmanuel Emmanuel Candes]<br> comes from his [http://www.ima.umn.edu/recordings/New_Directions_Short_Course/ND6.4-15.07/candes6-6-07.ram June 6 2007] video on the [[Optimization Videos]] page.
 +
 
 +
<pre>
 +
%Emmanuel Candes, California Institute of Technology, June 6 2007, IMA Summerschool.
 +
%Transcribed by Jon Dattorro.
 +
%Fails using SDP solver SDPT3 on 7th consecutive run after Matlab R2007b startup. CVX version 1.2 (build 656).
 +
%Fails using SDP solver Sedumi on 4th consecutive run after Matlab R2007b startup. CVX version 1.2 (build 656).
 +
clear all, close all
 +
n = 512; % Size of signal
 +
m = 64; % Number of samples (undersample by a factor 8)
 +
 
 +
k = 0:n-1; t = 0:n-1;
 +
F = exp(-i*2*pi*k'*t/n)/sqrt(n); % Fourier matrix
 +
freq = randsample(n,m);
 +
A = [real(F(freq,:));
 +
imag(F(freq,:))]; % Incomplete Fourier matrix
 +
 
 +
S = 28;
 +
support = randsample(n,S);
 +
x0 = zeros(n,1); x0(support) = randn(S,1);
 +
b = A*x0;
 +
 
 +
% Solve l1 using CVX
 +
cvx_quiet(true);
 +
%cvx_solver('sedumi');
 +
cvx_begin
 +
variable x(n);
 +
minimize(norm(x,1));
 +
A*x == b;
 +
cvx_end
 +
 
 +
norm(x - x0)/norm(x0)
 +
figure, plot(1:n,x0,'b*',1:n,x,'ro'), legend('original','decoded')
 +
</pre>
 +
 
 +
Code between <code>cvx_begin</code> and <code>cvx_end</code> requires [http://www.stanford.edu/~boyd/cvx CVX].
 +
 
 +
<code>randsample()</code> is from Matlab Statistics Toolbox.
 +
 
 +
Failure modes are reparable by [[Convex Iteration]]:
 +
<pre>
 +
%Emmanuel Candes, California Institute of Technology, June 6 2007, IMA Summerschool.
 +
%Convex Iteration implementation by Jon Dattorro.
 +
%Failure modes repaired.
 +
clear all, close all
 +
n = 512; % Size of signal
 +
m = 64; % Number of samples (undersample by a factor 8)
 +
 
 +
k = 0:n-1; t = 0:n-1;
 +
F = exp(-i*2*pi*k'*t/n)/sqrt(n); % Fourier matrix
 +
freq = randsample(n,m);
 +
A = [real(F(freq,:));
 +
imag(F(freq,:))]; % Incomplete Fourier matrix
 +
 
 +
S = 28;
 +
support = randsample(n,S);
 +
x0 = zeros(n,1); x0(support) = randn(S,1);
 +
b = A*x0;
 +
 
 +
cvx_quiet(true);
 +
%cvx_solver('sedumi');
 +
 
 +
%convex iteration
 +
y = ones(n,1);
 +
while 1
 +
% Solve l0 using CVX and Convex Iteration
 +
cvx_begin
 +
variable x(n);
 +
minimize(norm(y.*x,1));
 +
A*x == b;
 +
cvx_end
 +
 
 +
% update search direction y
 +
[x_sorted, indices] = sort(abs(x), 'descend');
 +
y = ones(n,1);
 +
y(indices(1:S)) = 0;
 +
 
 +
cardx = sum(abs(x) > 1e-6)
 +
if cardx <= S, break, end
 +
end
 +
norm(x - x0)/norm(x0)
 +
figure, plot(1:n,x0,'b*',1:n,x,'ro'), legend('original','decoded')
 +
</pre>

Revision as of 14:19, 16 November 2008

This Matlab demonstration of compressive sampling (a.k.a. compressed sensing) by Emmanuel Candes
comes from his June 6 2007 video on the Optimization Videos page.

%Emmanuel Candes, California Institute of Technology, June 6 2007, IMA Summerschool.
%Transcribed by Jon Dattorro.
%Fails using SDP solver SDPT3  on 7th consecutive run after Matlab R2007b startup.  CVX version 1.2 (build 656).
%Fails using SDP solver Sedumi on 4th consecutive run after Matlab R2007b startup.  CVX version 1.2 (build 656).
clear all, close all                
n = 512;                            % Size of signal
m = 64;                             % Number of samples (undersample by a factor 8)

k = 0:n-1;  t = 0:n-1;
F = exp(-i*2*pi*k'*t/n)/sqrt(n);    % Fourier matrix
freq = randsample(n,m);
A = [real(F(freq,:)); 
     imag(F(freq,:))];              % Incomplete Fourier matrix

S = 28;
support = randsample(n,S);
x0 = zeros(n,1);  x0(support) = randn(S,1);
b = A*x0;

% Solve l1 using CVX
cvx_quiet(true);
%cvx_solver('sedumi'); 
cvx_begin
    variable x(n);
    minimize(norm(x,1));
    A*x == b;
cvx_end

norm(x - x0)/norm(x0)
figure, plot(1:n,x0,'b*',1:n,x,'ro'), legend('original','decoded')

Code between cvx_begin and cvx_end requires CVX.

randsample() is from Matlab Statistics Toolbox.

Failure modes are reparable by Convex Iteration:

%Emmanuel Candes, California Institute of Technology, June 6 2007, IMA Summerschool.
%Convex Iteration implementation by Jon Dattorro.
%Failure modes repaired.
clear all, close all
n = 512;                            % Size of signal
m = 64;                             % Number of samples (undersample by a factor 8)

k = 0:n-1;  t = 0:n-1;
F = exp(-i*2*pi*k'*t/n)/sqrt(n);    % Fourier matrix
freq = randsample(n,m);
A = [real(F(freq,:));
     imag(F(freq,:))];              % Incomplete Fourier matrix

S = 28;
support = randsample(n,S);
x0 = zeros(n,1);  x0(support) = randn(S,1);
b = A*x0;

cvx_quiet(true);
%cvx_solver('sedumi');

%convex iteration
y = ones(n,1);
while 1
    % Solve l0 using CVX and Convex Iteration
    cvx_begin
        variable x(n);
        minimize(norm(y.*x,1));
        A*x == b;
    cvx_end

    % update search direction y
    [x_sorted, indices] = sort(abs(x), 'descend');  
    y = ones(n,1);
    y(indices(1:S)) = 0;

    cardx = sum(abs(x) > 1e-6)
    if cardx <= S, break, end
end
norm(x - x0)/norm(x0)
figure, plot(1:n,x0,'b*',1:n,x,'ro'), legend('original','decoded')
Personal tools