File Exchange

image thumbnail

CoSaMP and OMP for sparse recovery

version 1.7 (10.7 KB) by

Orthogonal Matching Pursuit (OMP) and Compressive Sampling Matched Pursuit (CoSaMP).



View License

Orthogonal matching Pursuit (OMP) and Compressive Sampling Matched Pursuit (CoSaMP) algorithm (see Needell and Tropp's 2008 paper ). This implementation allows several variants, and it also allows you to specify a matrix via function handles (useful if your matrix represents an FFT or similar).
A demo code shows how to use both the OMP.m and CoSaMP.m functions.
OMP and CoSaMP are useful for sparse recovery problems; in particular, they can be used for compressed sensing (aka compressive sampling), image denoising and deblurring, seismic tomography problems, MRI, etc.

Another good OMP implementation (C++, Matlab) is here:
(Updated, March 2012: SPAMS now has python and R bindings as well)

And a CoSaMP implementation (I haven't tested):
Edit: that CoSaMP implementation mentioned above is buggy. Read this:

Update, Feb 2012: for a blog discussion of several way to implement CoSaMP, see this website:

Comments and Ratings (22)

Ahmad Mousavi


VMat (view profile)

Thanks Stephen. Testing some cs greedy techniques, ideally I wish to plot the input signal and the reconstructed one: could you please tell what is the input signal variable in OMP.m and CoSaMP.m codes ?


mpvae (view profile)

good. thank you for sharing

Stephen Becker

Stephen Becker (view profile)

Thanks for noticing the bug Pham. I've updated the code.

Pham Dang

Dear M. Becker,
I suspect a small mistake at line 65 of test_OMP_and_CoSaMP.m. The line is :
z = sigma*randn(M,1); if COMPLEX, z = (z + sigma*randn(M,1))/sqrt(2); end

Should not it be :
z = sigma*randn(M,1); if COMPLEX, z = (z + i1*sigma*randn(M,1))/sqrt(2); end

in order to get a complex value for z ?
Best regards


Joo (view profile)

how can i use this code with images and get the output does any one know

fatemeh nj

I just have a bit problem, I receive this warning and my audio signal has not been recovered well, my warning is:
Warning: did not reach target size of residual

I also increase the number of 'K' but at the end the error was repeated!

my first question is that I know the k determined sparsity but how can i know the sparsity of unknown signal?

second question, could you also mention which paper you use for writing omp function?


Yang (view profile)

Thanks for this source code


TDJIO (view profile)

Thanks for sharing!


sheng (view profile)

Thank you very much for your share.


Jason (view profile)

How to estimate sparsity?

Works great if I know the sparse size exactly (by creating test data, for example), otherwise not such a great improvement over least squares (for my particular problem).


Jason (view profile)

Stephen Becker

Stephen Becker (view profile)

Israa, I'm not exactly sure what you want to do, but if you wish to work with 2D images, then use vectorize (e.g. use vec = @(x) x(:) ) and reshape operators, so that the code only sees vectors. If done right, it will work just fine. Hope that helps a bit.

israa tawfic

can i use this source code with an image,

Stephen Becker

Stephen Becker (view profile)

Robin, that's a typo. It's an unimportant default value that doesn't affect the code, so I will change it the next time I have a major update.

Why do you use round with 2 nargin in CoSaMP line 97?

It doesn't work.

Thanks for this source code.


Bo (view profile)

I just get it work.
Sorry for the 4 stars. It should be a five stars.
Thank you,
Bo Yuan


Bo (view profile)

Hi Stephan,

It seems your functions can only accept A which is square matrix? Do I miss anything or this function will not work for overcomplete frames?

Thank you!
Bo Yuan

Stephen Becker

Stephen Becker (view profile)

good catch Noam! I'm updating the fixed code. The new version also displays some more information for CoSaMP when using function handles and using LSQR.

sasikala mr


Noam Wagner

Hi Stephan,
Thank you very much for the very nice and useful OMP code.
I believe there might be a mistake in the orthogonalization loop of "atom_new", when working in fast mode. I think the equation:

atom_new = atom_new - (atom_new'*A_T(:,j))*A_T(:,j);

should be modified to:

atom_new = atom_new - (A_T(:,j)'*atom_new)*A_T(:,j);

This turns out significant when working with complex numbers. Otherwise, one may notice different behavior in slow and fast modes.
Thanks again!

Arunima c.v.




Fixed bug with complex numbers (Aug 2016)


Fixing the bug for complex mode. Minor changes to the solvers. Added complex data test mode to the test script.


Fixing a bug that affected versions of Matlab prior to 2009b. See


editing description text a bit


Adding new links in the description, and updated the demo file slightly.

MATLAB Release
MATLAB 7.10 (R2010a)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video