Friday, September 21, 2007

Compressed Sensing: How to wow your friends.


Here is a short example showing the power of Compressed Sensing. I modified a small program written by Emmanuel Candes from his short course at IMA ("The role of probability in compressive sampling", Talks(A/V) (ram)).
Before you continue on, you need to have Matlab and SeDuMi installed on your system. Let us imagine you have a function f that is unknown to you but for which you know it is sparsely composed of sine functions. Let us imagine for the sake of argument that:

f(x)=2*sin(x)+sin(2x)+21*sin(300x)

The normal way of finding what this decomposition is by computing the scalar product of this function f with sines and cosines with as large as possible of frequency content (in order to capture the high frequency content of this function). This is done iteratively because one doesn't know in advance that a frequency of 300 is part of the solution. But in effect, one is solving a least square problem.

With Compressed Sensing, we know that if we are to measure this function with an incoherent basis to sines and cosines, we are likely to find the decomposition using the L1 minimization. We know that diracs and sines are incoherent (please note we are not using random projections), so we evaluate the scalar product of f with several diracs functions centered at say 18 different locations. In effect, we query for the value of f(x) at the following points x -chosen at random by hand- : 1, 4, 6, 12, 54, 69, 75, 80, 89, 132, 133, 152, 178, 230, 300, 340, 356, 400
We then use these values to solve for the coefficients of the sine series expansion of f(x) by doing an L1 minimization using SeDuMi:

clear
% size of signal
n = 512;
x0 = zeros(n,1);
% function is f(x)=2*sin(x)+sin(2x)+21*sin(300x)
x0(1)=2;
x0(2)=1;
x0(300)=21;
% evaluating f at all sample points xx
% f(1), f(4), f(6).....f(340) f(356) f(400)
xx=[1 4 6 12 54 69 75 80 89 132 133 152 178 230 300 340 356 400];
% C is measurement matrix
for i=1:n
C(:,i)=sin(i*xx)';
end
b = C*x0;
% b is the result of evaluating f at all sample points xx
% f(1), f(4), f(6).....f(340) f(356) f(400)
% C is the measurement matrix
% let us solve for x and see if it is close to x0
% solve l1-minimization using SeDuMi
cvx_begin
variable x(n);
minimize(norm(x,1));
subject to
C*x == b;
cvx_end
figure(1)
plot(abs(x-x0),'o')
title('Error between solution solution found by L1 and actual solution')
figure(2)
plot(x,'*')
hold
plot(x0,'o')
title('Solution found for x using L1 minimization')

With only 20 function evaluations we have a near exact reconstruction. You may want to try the the L2 norm (or least-square) by replacing:
minimize(norm(x,1));

by
minimize(norm(x,2));
and then change the number of function evaluations needed to reach the same result. With 200 functions evaluations, the error between the reconstructed solution and the actual solution is several orders of magnitude larger than the L1 technique with 20 functions evaluations.
Pretty neat uh ?


Liked this entry ? subscribe to the Nuit Blanche feed, there's more where that came from

If you think this blog provides a service, please support it by ordering through the Amazon - Nuit Blanche Reference Store

38 comments:

Anonymous said...

Is the effect as impressive when noise is added?

I did

b = C*(x0+epsilon*randn(size(x0)));

and for epsilon = 0.001, it broke the algorithm down. This makes sense because the resulting vector is no longer sparse in the dictionary chosen; however, this situation happens in practice all the time. Any way to combat this theoretically or algorithmically?

Igor said...

Hello Anonymous,

Yes, it has to break down. The example is featured to show
- that you don't need to have a random projection to find a solution to compressed sensing.
- It is also there to show that how profoundly it would be different to solve an integral or a PDE using this type of technique instead of using a weak formulation (L2 based) as we have been accustomed to do for the past two hundred years.
- I also wanted to illustrate the issue of coherence. In this case, the noise is coherent with the diracs, so the sines cannot make a difference between the signal you want (diracs) and the noise.

I agree with you: noise happens all the time. But in most real cases, you don't have diracs but rather gaussians or wavelet like functions. One really is then looking for families of function that are incoherent with the gaussians or wavelets and coherent with the noise. Noiselets seem to fit that bill.

http://nuit-blanche.blogspot.com/search?q=noiselet

R. Coifman, F. Geshwind, and Y. Meyer. Noiselets. Appl. Comp. Harmonic Analysis, 10:27–44, 2001

Igor said...

Evidently, if the signal is made of gaussians or wavelets, random projections are ok. Random basis are indeed incoherent with the signal to be studied but are coherent with the noise. Results of Donoho, Candes and others show the ability to retrieve the signal of interest under random projections when noise is present.

Igor.

Jing said...

Interesting staff.

Two questions:

Where is the spot above 400 comes from? The highest frequency should be 300, isn't it? I try adding different frequency components to y, if the highest frequency in y is less than a certain value, this unexpected spot gone, otherwise, this spot will always be here. Look similar to the frequency alias when sampling a signal at a rate that is less than its Nyquist-rate.

Also, this algorithm seems to be sensitive to the frequency, I replace sin(2*x) with sin(2.01*x) in y, this algorithm fails to detect this frequency component. Any explain to this?

Thanks.

Igor said...

Jing,

I know I am answering this very late but for the second part of your question, please note that you cannot find a frequency that is not in the dictionnary in the first place. All frequencies are integer from 1 to 512 and do not include 2.01, so it cannot find it.

I think the answer to the first part of your question is similar.

Igor.

Anonymous said...

??? Undefined function or variable "cvx_begin".

Error in ==> how_to_wow_your_friends_CS at 22
cvx_begin


I am using MATLAB 7.0 and the sedumi is installed

Igor said...

Hello Anonymous,

Even though sedumi is installed, are you sure it works, or that it is in your path ? It looks to me like a sedumi specific error, you might want to contact them directly.

So that you know, I am going to try to make that example sedumi agnostic in the future since we now have at least in-house solvers ( http://igorcarron.googlepages.com/cscodes#sp )

Igor.

Anonymous said...

How can I use this 'example' if I have a 32 bit pc (which is suppose to avoid me for using Sedumi and so the optimiztion L1, L2, etc?
thanks,
Juan C.

Igor said...

Juan,

Sedumi works in 32 bit. You just need to use Sedumi 1.1.

Cheers,

Igor.

Anonymous said...

This example is somewhat confusing for me.

Your "measurement matrix" should actually be a bunch of diracs, which are multiplied by f(x) in order to derive b, right? You just didn't show this part.

You then replaced this "measurement matrix" with a new matrix C that contains sine waves (your basis vectors) and used this to derive the coefficients of f(x) in the basis of sine waves. Let's call C the "new matrix".

How did you know which sine waves to include in your new matrix C? Why does having a dirac(4) in your measurement matrix mean that you have to have a sine(4x) in your new matrix? I didn't know there was a relationship between dirac(4) in the time domain and dirac(4) in the frequency domain. Am I missing something?

Igor said...

hello Anonymous,

Generally, in CS, the system of equations is:

y = \phi * \psi * alpha

I am making the assumption that the unknown function f is really a composition of several sine functions that are listed in the \psi matrix or:

f = \psi * alpha

In other words, when alpha is sparse, f is sparse in the family of sine functions.

\phi is a dictionary of incoherent functions (incoherent with respect to sines). It is made up of diracs centered on integer points of the w-axis.

The matrix C you are refering to is just \phi * \psi

Elements of C are really the inner products of sin(w_j * x) and dirac(w - w_i).

y = C * alpha

and

f(x) = sum(alpha_j * sin(w_j * x))

alpha (made of alpha_j) is the list of coefficients (sparse) representing the frequency content of function f.

Hope this helps, let me know if this is not clear.

Cheers,

Igor.

海波 said...
This comment has been removed by the author.
Anonymous said...

Is it possible to equivalently look at the problem as that of measuring a signal that consists of impulses at different time instances? Then, C, the measurement matrix is equivalent to the Fourier measurement matrix, thus measuring the signal at random modulations. The solution of ell-1 minimization would then be the original signal, which is sparse in time.

Igor said...

Sure it could.

Igor.

Anonymous said...

Hello Igore,

measurement matrix(C) don't satisfies RIP condition :

n = 512;
x0 = zeros(n,1);
x0(1)=2;
x0(2)=1;
x0(300)=21;
xx=[1 4 6 12 54 69 75 80 89 132 133 152 178 230 300 340 356 400];
for i=1:n
C(:,i)=sin(i*xx)';
end
b = C*x0;

>> norm(x0)

ans =

21.1187

>> norm(b)

ans =

60.1531

why?,

thanks.

James said...

Regarding the cvx_begin error above...

I had this problem too. You need to install CVX:

http://www.stanford.edu/~boyd/cvx/download.html

This is not part of SeDuMi - I guess you need both installed?

It's working for me now. I use Matlab 7.8 on a macbook pro. Once again, it only started working once I installed SeDuMi and CVX.

Cheers!

Bhaskar said...

I don't know if it is too late to ask question. If I replace x0(2)=1 by x0(5)=10, it does not work. I also tried adding a fourth frequency component and it failed. Why?

With x0(5) in there, if I specify a sampling point at 30, making the xx array 19 elements long, it works again.

Igor said...

Bhaskar,

If you increase the number of samples, you will likely have convergence most of the time. 18 samples is kind of minimum but with 50 you should always be able to recover the solution.

Hope this helps,

Igor.

Anonymous said...

Hi igor. I tryed to use your algorithm and i found a matlab error, can you help me?
?? Error: File: C:\matlabR12\work\how_to_wow_your_friends_CS.m Line: 23 Column: 18
Missing operator, comma, or semicolon.

Regards
Mark

Igor said...

Mark,

Please make sure you are doing the right cut and paste. Sometimes, the procedure of cutting and pasting bring extra characters invisible to the eye.

Cheers,

Igor.

Marcos said...

Thanks Igor.
It is not working on my MLab, I typed line by line and nothing. Is that possible you send me yoir .m file?

Regards Mark

Anonymous said...

Thank you James!

Anonymous said...

Hopefully my question is not too later. This is a fantastic program. I modified the program a bit:

change sin(i*xx) to sin(pi*1*xx)


And the program cannot find the right answer. Can you explain it?

Anonymous said...

hello igor
can you tell me how can i use CS on images?

Igor said...

Think of an imgae as a matrix, decompose tha matrix as a series of vectors (i.e. Concatenate these vectors into one large vector) and tada yo have linear system of equation
Ax = b.

Anonymous said...

yeah i have already done that. but i want to do it without braking image into vectors. like if i take dct2 of an image and then multiply with sensing matrix, then how can i reconstruct it?
regarda abiha!

Igor said...

If you take the DCT2 you have already transformed your image, so i am not sure i understand. Have you looked through any of he courses listed below:
http://nuit-blanche.blogspot.com/p/teaching-compressed-sensing.html

Igor.

Anonymous said...

if the image is not sparse in time domain and it is sparse in some other domain, like i supposed that if i take dct2 and get the sparse image, then multiply with sensing matrix, recover and then take idct2.. but the problem is that how can i recover a matrix instead of a vector? and how can i define a sending matrix for a matrix?
i hope i have made it clear...
waiting for reply..

Igor said...

The least you could do is acknowledge that you have read the web page I mentioned earlier, have you ?

Anonymous said...

not the one you have mentioned,but i have read about compressive sensing from other sites.

Igor said...

Not the right ones obviously. Please read and find in the link the answer to your question. If you cannot find please come back and let me know what you have read and point to me exactly what is it you don't understand.

Cheers,

Igor.

Reuben said...

Thanks for the helpful example, Igor! I'm just wondering why you chose to take 18 samples. I read that there is a practical four-to-one rule that says you need about four incoherent samples per unknown nonzero term for exact recovery. That suggests 12 should be enough here, but I've experimented and it isn't. I usually don't get exact recovery unless I use 18 or more random incoherent samples.

Anonymous said...

Hi Igor,
I´m trying to use wavelets or splines in the psi basis, until now it doesn´t work. Do you have any suggestion? Thanks,
Sophie

Igor said...

Reuben,

This is really a rule of thumb for the constant in CKlog(N/K). the Clog(N/K) is about 4 or 5 most of the time. Your mileage **will vary** depending on the solver.

Sophie,

the matrix used here is already a combination of phi and psi. To get into random for phi and wavelet for psi, you really want to look at the examples of SPGL1 ( http://www.cs.ubc.ca/~mpf/spgl1/examples/spgexamples.html )

Igor.

Anonymous said...

Respected sir, i am the student of BSA university, and i am the beginner for this field compressive sensing, but , i am very much interested in compressive sensing, I have some doubts on the basics of cs,
i.e ., we have taken one input (either 1D or 2D) called X, then, what is sparse matrix, what is transformed matrix, what is measurment matrix, random generated matrix, psi matrix , phi matrix,,,
and what is these and why we are doing these?
sorry for the disturbance and expecting your reply soon,
thank you so much,,,,,

Igor said...

Anonymous,

Please get acquainted with the topic by reading the refeences at : http://nuit-blanche.blogspotcomr/p/teaching-compressed-sensing.html

Igor.

Anonymous said...

Thank you very much sir,

Anonymous said...

Hi. I had a question about the example code. I have run it with sedumi which now comes with cvx and I am not getting the same output as above. When I run as is the code only goes through about 11 iterations with good results. Error is small on the order of 10^-9. When I change the example to attempt the least squares suggestion the code ends after only 4 iterations and none of the results have strayed very far from zero, so the error is roughly equivalent to the actual values in the known function that we are searching for. It seems to be quitting too soon before even approaching a result. Is this an issue with modern versions of matlab and cvx? I am running Matlab2013a and cvx 2.0. I do not know which version of sedumi is included here.

Also, I was wondering if there is an extremely simple paper and pencil example of this process. Perhaps with matrices that do not stray too far into the double digits if possible. I tried working out the above, generally, by hand, but the sines make it somewhat cumbersome. Most of the examples I see are very engineering or CS centered and I have little background in either.

I've been reading through the information available on reddit, this site, and other blogs and I still do not have the process intuitively down.

Thanks,
Roo