Chapter_07_ShowTheGap1.m

% Figure - 7.1
% =========================================
% This program shows the actual performance of BP and OMP versus 
% the weak bounds we got for them - the noiseless case


function [ErT,A]=Chapter_07_ShowTheGap1(n,m,Smin,Sd,Smax,Experiments)

%==========================================
%  The figure in Chapter 7 was generated by
%           [ErT,A]=Chapter_07_ShowTheGap1(100,200,1,1,70,100);
%==========================================

A=randn(n,m);
W=sqrt(diag(A'*A));
for k=1:1:m, 
    A(:,k)=A(:,k)/W(k); 
end; 
ErT=zeros(Smax-Smin+1,2);
h=waitbar(0,'experimenting ...');

countT=length(Smin:Sd:Smax)*Experiments;
count=1;
for S=Smin:Sd:Smax,
    disp(['=======> ',num2str(S)]);
    Er=zeros(Experiments,2);
    for k=1:1:Experiments
        waitbar(count/countT);
        x=GenerateX(m,S);
        b=A*x;
        xBP=BasisPursuit(A,b,x);
        Er(k,1)=mean((xBP-x).^2)/mean(x.^2);
        xMP=MatchingPursuit(A,b);
        Er(k,2)=mean((xMP-x).^2)/mean(x.^2);
        count=count+1;
    end;
    ErT(S-Smin+1,1)=mean(Er(:,1)<1e-5);
    ErT(S-Smin+1,2)=mean(Er(:,2)<1e-5);
    CreateGraph(Smin,Sd,Smax,ErT); 
    drawnow;
end;
close(h); 
% print -depsc2 Chapter_07_Comp1.eps

%=================================================

function x=GenerateX(m,S)

x=zeros(m,1);
pos=randperm(m);
x(pos(1:S))=randn(S,1);
return;

%=================================================

function xBP=BasisPursuit(A,b,x)

[n,m]=size(A);
V=ones(2*m,1);
xBP=linprog(V,[],[],[A, -A],b,0*V,V*10);
xBP=xBP(1:m)-xBP(m+1:end);
return;

%=================================================

function xMP=MatchingPursuit(A,b)

[n,m]=size(A);
r=b;
SS=[];
while r'*r>1e-8,
    Z=(A'*r).^2;
    pos=find(Z==max(Z));
    SS=sort([SS,pos(1)]);
    r=b-A(:,SS)*pinv(A(:,SS))*b;    
end;
xMP=zeros(m,1);
xMP(SS)=pinv(A(:,SS))*b;
return;

%=================================================

function []=CreateGraph(Smin,Sd,Smax,ErT)

figure(1); clf; 
h=plot(Smin:Sd:Smax,ErT(1:Sd:end,1),'r'); hold on; 
set(h,'LineWidth',2);
h=plot(Smin:Sd:Smax,ErT(1:Sd:end,2),'b'); 
set(h,'LineWidth',2);
h=xlabel('Cardinality of the solution'); set(h,'FontSize',14);
h=ylabel('Probability of success'); set(h,'FontSize',14);
h=legend({'Basis Pursuit','Matching Pursuit'},3);
set(h,'FontSize',14);
h=plot(Smin:Sd:Smax,ErT(1:Sd:end,2),'bo');
set(h,'MarkerFaceColor','b');
h=plot(Smin:Sd:Smax,ErT(1:Sd:end,1),'ro');
set(h,'MarkerFaceColor','r');
set(gca,'FontSize',14);
grid on;
axis([0 70 -0.1 1.1]);