# Candes.m

(Difference between revisions)
 Revision as of 23:29, 15 November 2008 (edit) (wIdHeYggpGXTenP)← Previous diff Revision as of 14:19, 16 November 2008 (edit) (undo)m (Reverted edits by 84.244.189.99 (Talk); changed back to last version by Ranjelin)Next diff → Line 1: Line 1: - jGxq35 zxfvgjjlvqej, [url=http://xerodzoxagxo.com/]xerodzoxagxo[/url], [link=http://zkasnuzddvmk.com/]zkasnuzddvmk[/link], http://uuspekpdbwgk.com/ + This Matlab demonstration of ''compressive sampling'' (a.k.a. ''compressed sensing'') by [http://www.acm.caltech.edu/~emmanuel Emmanuel Candes]
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. + +
+                                                                                                                                                     %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 [http://www.stanford.edu/~boyd/cvx 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')
+

## 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')
```