CODEWITHSUNDEEP
matlabimage-processingsparse-matrix
P.Bisoyi
P.Bisoyi

Reputation: 133

Implementing matching pursuit algorithm

I have implemented matching pursuit algorithm but i m unable to get the required result.

Here is my code:

D=[1 6 11 16 21 26 31 36 41 46
    2 7 12 17 22 27 32 37 42 47
    3 8 13 18 23 28 33 38 43 48
    4 9 14 19 24 29 34 39 44 49
    5 10 15 20 25 30 35 40 45 50];
b=[6;7;8;9;10];
n=size(D);
A1=zeros(n);
R=b;
H=10;
if(H <= 0)
    error('The number of iterations needs to be greater then 0')
end;
for k=1:1:H
    [c,d] = max(abs(D'*R));    %//'
    A1(:,d)=D(:,d);
    D(:,d)=0;
    y = A1\b;
    R = b-A1*y;
end

Output 

y=

0.8889  
     0  
     0  
     0  
     0  
     0  
     0  
     0  
     0  
0.1111

I should get only non-zero value at (2,1) and other values should be zero but I'm getting 2 non-zero value. Can you please help me find out where the error is?

Thanks.

Upvotes: 2

Views: 2832

Answers (3)

Shailesh Kumar
Shailesh Kumar

Reputation: 6957

Yes the atoms in the dictionary must be normalized so that the inner products of the current residual with different atoms can be compared fairly.

You may want to check my OMP implementation which also includes incremental Cholesky updates for the least square step in OMP at https://github.com/indigits/sparse-plex/blob/master/library/%2Bspx/%2Bpursuit/%2Bsingle/omp_chol.m

I have written detailed tutorial notes on OMP in my library documentation at https://sparse-plex.readthedocs.io/en/latest/book/pursuit/omp/index.html

My library sparse-plex contains a C implementation of OMP which is close to 4 times faster than fastest MATLAB implementations. See the discussion here https://sparse-plex.readthedocs.io/en/latest/book/pursuit/omp/fast_omp.html

Upvotes: 0

Xcross
Xcross

Reputation: 2132

Here is the updated code. This is based on the algorithm provided at https://en.wikipedia.org/wiki/Matching_pursuit

clc;
clear all;

D=[1 6 11 16 21 26 31 36 41 46
    2 7 12 17 22 27 32 37 42 47
    3 8 13 18 23 28 33 38 43 48
    4 9 14 19 24 29 34 39 44 49
    5 10 15 20 25 30 35 40 45 50];

b=[6;7;8;9;10];

H=10;

for index=1:10
   G(:,index)=D(:,index)./norm(D(:,index)); 
end
G1=G;

n=size(G);
R=b;
if(H <= 0)
    error('The number of iterations needs to be greater then 0')
end;
if(H >size(D,2))
    error('The number of iterations needs to be less than dictionary size')
end;

bIndex=1:size(G,2);
for k=H:-1:1
    innerProduct=[];
    for index=1:size(G,2)
        innerProduct(index)=dot(R,G(:,index));
    end
    [c,d] = max(abs(innerProduct));    
    An(H-k+1)=innerProduct(d);
    R = R-(An(H-k+1)*G(:,d));
    G(:,d)=[];
    strong(H-k+1)=bIndex(d);
    bIndex(d)=[];
end
G_new=G1(:,strong);


%% reconstruction
bReconstructed=zeros(size(G_new,1),1);
for index=1:size(G_new,2)
    bReconstructed(:,index) = (An(index)*G_new(:,index));
end
b_new=sum(bReconstructed,2)

Upvotes: 1

Steffen
Steffen

Reputation: 2431

I checked with: http://www.scholarpedia.org/article/Matching_pursuit

Your functions need to be normalized!

D = D./repmat(sum(D,1),5,1);

I get the following algorithm:

D=[1 6 11 16 21 26 31 36 41 46
    2 7 12 17 22 27 32 37 42 47
    3 8 13 18 23 28 33 38 43 48
    4 9 14 19 24 29 34 39 44 49
    5 10 15 20 25 30 35 40 45 50];
D = D./repmat(sum(D,1),5,1);
b=[6;7;8;9;10];
n=size(D);
A1=zeros(n);
R=b;
H=100;
if(H <= 0)
    error('The number of iterations needs to be greater then 0')
end;
a = zeros(1,H);
G = zeros(size(D,1),H);
for k=1:1:H
    ip = D'*R;
    [~,d] = max(abs(ip));    %//'
    G(:,k) = D(:,d);
    a(k) = ip(d);
    R = R-a(k)*G(:,k);
end

% recover signal:
Rrec = zeros(size(R));
for i=1:H
    Rrec = Rrec + a(i)*G(:,i);
end
figure();
plot(b);
hold on;
plot(Rrec)

It approximates the signal quite well. But not with D(:,2) at first as expected. Maybe it is a starting point...

Upvotes: 1

Related Questions