[Warning: Before you read this entry, read this one first ]I know this is not Monday, but if you read this on Monday, it'll look like it was posted on Monday. Today's algorithm has to do with Compressed Sensing. In Compressed Sensing, the idea is that you can sample a signal much below the Nyquist rate if you have the prior knowledge that it is sparse. Since most signal are sparse, this technique should eventually become mainstream. If you want to know more about Compressed Sensing check the Big Picture out.
If you have been obtaining Compressed Sensing measurements, the most difficult task is then to decode that information into the actual signal. You can use any of the amazing reconstruction codes produced so far or you could implement a simple one proposed by Wei Dai and Olgica Milenkovic in Subspace Pursuit for Compressive Sensing: Closing the Gap Between Performance and Complexity. [Update: Wei Dai provided me with the original code used in the paper and asked me to put it on this site ]
The algorithm is pretty straightforward from the explanation given in the paper. Here is the SP.m file but the full set of m files implementing this technique can be found either in the local CS code section or in the Monday Morning Algorithm archive. Once you have downloaded all the functions in one directory, you can run directly trialSP from the prompt. It runs one simple example of a vector of 290 elements and 3 non-zero elements. As expected by the log(K) asymptotic, I am pretty much amazed at the speed of this thing. If you find a bug please do let me know.
function [xfinal,That]=SP(K, phi, y,nn)
%
% This is an example using the Subspace Pursuit
% Algorithm of Wei Dai and Olgica Milenkovic
% "Subspace Pursuit for Compressive Sensing: Closing the
% Gap Between Performance and Complexity"
%
% Written by Igor Carron
% http://nuit-blanche.blogspot.com
% Creative Commons Licence
% Initialization
%
u = phi' * y;
[That,phiThat]=deft(phi,u,K);
y_r = resid(y,phiThat);
%
%
% Iteration
%
while y_r ~= 0
u1 = phi' * y_r;
[Tprime1,phiTprime1]=deft(phi,u1,K);
[Tprime,phiTprime]=uni(Tprime1,That,phi);
phiTprimecross = cros(phiTprime);
x_pprime = phiTprimecross * y;
[Ttilde,phiTtilde] = deft(phi,x_pprime,K);
ytilde_r = resid(y,phiTtilde);
if gt(norm(ytilde_r) ,norm(y_r))
break
else
That = Ttilde;
y_r = ytilde_r;
end
end
%
% Output
%
phiThatcross = cros(phiThat);
x_That = phiThatcross * y;
That;
xfinal = zeros(nn,1);
for ii = 1:K
xfinal(That(ii))=x_That(ii);
end
Hi,
ReplyDeleteI am new to your weblog.
I have a silly question and would be thankful if you answer it:
Why should somebody do that, we can simply sort the signal values which has the cost of ~O(Nlog N) and threshols it with zero (in the absence of a noise) which has the cost of O(N). In the presence of noise, we can simply sort the signal and find an estimate of variance of the noise from the tail of the sorted signal and pick some threshold base on that. Do you mean this algorithm improve the complexity from this to O(K)? but both of them are modestly fast and not that expensive (of course O(K) is faster) but is it the only reason?
I know there must be a reason people doing that but I don't know why.
Sorry for such a silly question,
would you please give me some elementary references to start reading about it?
Xiaoying
Hello Xiaoying,
ReplyDeleteif you acquire all the samples, sort them and then threshold them, you are in effect doing traditional sampling. In compressed sensing, sampling is reduced so you don't have to sort, threshold and throw away. You can do this because you have some prior information that the signal is sparse. The sample size is limited to about O(Klog(K)) instead of N. But in order to get the full signal back from O(klog(k)) samples to N data, you have to perform the reconstruction following this algorithm.
You can find more information on the subject of compressed sensing (tutorial/videos, etc...) here:
http://igorcarron.googlepages.com/cs#understanding
Hope this helps,
Igor.
Hello Igor,
ReplyDeleteI ran into a problem using your code.
It seemed that the algorithm got stucked somewhere.
The problem was the line:
if gt(norm(ytilde_r) ,norm(y_r))
For my special dataset, norm(ytilde_r) and norm(y_r) where equal sometimes, therefore the while loop was repeated infinitely. I replaced the > by >= and it works.
Maybe you want to change that too.
Have a nice day,
Sebastian