File Exchange

image thumbnail

CoSaMP and OMP for sparse recovery

version 1.7 (10.7 KB) by Stephen Becker
Orthogonal Matching Pursuit (OMP) and Compressive Sampling Matched Pursuit (CoSaMP).


Updated 05 Aug 2016

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:

Cite As

Stephen Becker (2020). CoSaMP and OMP for sparse recovery (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (25)

feng yanfei

There is a wrong in CoSaMP function.
In the CoSaMP, the 123rd code should be ' if k > M '.
Anyway thank you~!


Lei Chen


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 ?


good. thank you for sharing

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


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?


Thanks for this source code


Thanks for sharing!


Thank you very much for your share.


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).


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.

can i use this source code with an image,

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.


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


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

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!




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 Compatibility
Created with R2010a
Compatible with any release
Platform Compatibility
Windows macOS Linux