# Candes.m

### From Wikimization

(Difference between revisions)

(wIdHeYggpGXTenP) |
m (Reverted edits by 84.244.189.99 (Talk); changed back to last version by Ranjelin) |
||

Line 1: | Line 1: | ||

- | + | 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')