Thursday, December 13, 2007

Compressed Sensing: Single-Pixel Imaging via Compressive Sampling, Stable Manifold Embedding

In an upcoming issue of the IEEE Signal Processing Magazine, we will be treated with a Special Issue on Compressive Sampling. It is scheduled to come out in March 2008. In it, is a new well written review of the trade studies performed on the current single pixel camera development using Compressed Sensing at Rice. The review is Single-Pixel Imaging via Compressive Sampling and was written by Marco Duarte, Mark Davenport, Dharmpal Takhar, Jason Laska, Ting Sun, Kevin Kelly, and Richard Baraniuk. The abstract reads:
Humans are visual animals, and imaging sensors that extend our reach – cameras – have improved dramatically in recent times thanks to the introduction of CCD and CMOS digital technology. Consumer digital cameras in the mega-pixel range are now ubiquitous thanks to the happy coincidence that the semiconductor material of choice for large-scale electronics integration (silicon) also happens to readily convert photons at visual wavelengths into electrons. On the contrary, imaging at wavelengths where silicon is blind is considerably more complicated, bulky, and expensive. Thus, for comparable resolution, a $500 digital camera for the visible becomes a $50,000 camera for the infrared. In this paper, we present a new approach to building simpler, smaller, and cheaper digital cameras that can operate efficiently across a much broader spectral range than conventional silicon-based cameras. Our approach fuses a new camera architecture based on a digital micromirror device with the new mathematical theory and algorithms of compressive sampling . CS combines sampling and compression into a single nonadaptive linear measurement process [1–4]. Rather than measuring pixel samples of the scene under view, we measure inner products between the scene and a set of test functions. Interestingly, random test functions play a key role, making each measurement a random sum of pixel values taken across the entire image. When the scene under view is compressible by an algorithm like JPEG or JPEG2000, the CS theory enables us to stably reconstruct an image of the scene from fewer measurements than the number of reconstructed pixels. In this manner we achieve sub-Nyquist image acquisition. Our “single-pixel” CS camera architecture is basically an optical computer (comprising a DMD, two lenses, a single photon detector, and an analog-to-digital (A/D) converter) that computes random linear measurements of the scene under view. The image is then recovered or processed from the measurements by a digital computer. The camera design reduces the required size, complexity, and cost of the photon detector array down to a single unit, which enables the use of exotic detectors that would be impossible in a conventional digital camera. The random CS measurements also enable a tradeoff between space and time during image acquisition. Finally, since the camera compresses as it images, it has the capability to efficiently and scalably handle high-dimensional data sets from applications like video and hyperspectral imaging. This article is organized as follows. After describing the hardware, theory, and algorithms of the single-pixel camera in detail, we analyze its theoretical and practical performance and compare it to more conventional cameras based on pixel arrays and raster scanning. We also explain how the camera is information scalable in that its random measurements can be used to directly perform simple image processing tasks, such as target classification, without first reconstructing the underlying imagery. We conclude with a review of related camera architectures and a discussion of ongoing and future work.
When I last explained the single pixel camera, I didn't realize the first option is also called the raster scan. The noteworthy result of the trade study compares different current scanning capabilities and that provided by compressed sensing.They look at dynamic range, quantization error, photon counting noise and remind us of the fact that with CS, one can obtain much more light on the sensitive detector (pixel) than conventional cameras and that it should help in low light conditions as presented in the IMA presentation (Multiscale Geometric Analysis).



X-rays could be implemented the same way if the DMD would be able to reflect at these wavelengths, for instance by using a well thought multilayer system for instance. They also give a good summary of the other hardware technologies implementing Compressed Sensing. It looks like I have not missed any (part 1, part 2, part 3, part 4)

In another good summary of the current investigations of Richard Baraniuk's group, I found the following slide most illuminating (as presented in Single-Pixel Imaging via Compressive Sampling,) explaining the geometric reason why L1 works

The same figure could be used to explained the relative better ability of non-convex Lp regularization techniques to do even better than L1. Because the Lp (for p less than 1) ball shape tend to be squashed compared to the L1 ball, it is easier to see that there are more randomly oriented N-M planes anchored on a K-face of the Lp ball which would not interesect with the ball. [Update: Muthu MuthuKrishnan just posted a nice picture of an Lp ball, the picture is due to Petros Boufounos ].




With regards to extension of compressed sensing or compressed sensing ideas, Richard Baraniuk and Mike Wakin [1] have these three illuminating slides where they shows how Random projections can reduce Manifolds of low dimension in high dimensional space or how the condition for which random projection can be used as dimensionality reduction techniques. It shows how CS provides better bounds than the older Whitney's embedding theorem.




The current drawback of this approach (random projections) is the inability to provide some type of parametrization on the manifold. Never say never.

Source: Rice Compressed Sensing Repository.
Reference : [1] R. G. Baraniuk and M. B. Wakin, Random Projections of Smooth Manifolds, to appear in Foundations of Computational Mathematics.

Tuesday, December 11, 2007

Complex Systems that Matter: an unfortunate Nuclear reactor shutdown and the Space Shuttle ECO malfunction


Soon after September 11th, I recall a short discussion with Alan Waltar who was the head of our department on the amount on some of the hate mail he was getting because of his nuclear engineering leadership. I recall specifically one e-mail where he was touted as a terrorist. I cannot say how extremely annoyed I was when I read it. Alan is a truly nice guy and is very passionate about a subject that is becoming increasingly relevant everyday. In America the Powerless: Facing our Nuclear Energy Dilemma he was mentioning the important issue on our dramatic reliance on Canadian and other foreign reactors to get medical isotopes ( isotopes of molybdenum, technetium ). Then I read this yesterday about a Nuclear plant shutdown brings hospital delays where it explains that because of a unduly conservative regulatory environment, the shutdown of one single Canadian ( Chalk River nuclear facility ) will delay significantly the wait for medical tests in the U.S. and other places. Let us not kid ourselves, it will have some dramatic consequences on some patients.
In line with dealing with very complex systems, here is the internal NASA e-mails from the director of Shuttle Safety at the Johnson Space Center and the Shuttle Program Manager show how they struggle to understand and address the failure of the engine cutoff (ECO) sensor malfunction (Rainer Gerhards has a more detailed view of the system) of the Space Shuttle (In NASA's lingo, the Shuttle is not just the plane, it is the plane and the boosters) that postponed last week's launch. This is not the first time we see these internal passionate discussions about such a detailed system, but as usual, I am always left with the bitter taste that if NASA were to provide some type of funding for what they call Red team reviews with an emphasis of funding people outside of the space business, analyses would be a lot more fundamental. This tiny investment would go a long way toward sustaining and improving the program.




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

Monday, December 10, 2007

Compressed Sensing: Blog entries and website updates.

Terry Tao has an interesting entry on random matrices in Additive combinatorics and random matrices. At the same time, Muthu Muthukrishnan has an entry on Compressed Sensing, entitled: Compressed Sensing Enriched. Also, Rick Chartrand just tells me the article mentioned in the Restricted Isometry Properties and Nonconvex Compressive Sensing entry is now available here. Michael Grant, Stephen Boyd, and Yinyu Ye released a new build (Version 1.1, Build 577) for CVX ( Matlab Software for Disciplined Convex Programming), they also seem to have a new introduction page providing a good explanation of what the software does to newcomers.

And remember kids, Cassini will flyby Titan on Dec. 20 or ten days from now. It will fly above Titan at an altitude of 970 km above it. Not the 1 million km of this photo:


Photo credit: NASA/JPL/Space Science Institute

Monday Morning Algorithm Part 6: The Fast Johnson-Lindenstrauss Transform

Last year, Nir Ailon and Bernard Chazelle described a new algorithm that improves the ability to perform the Johnson-Lindenstrauss transform. It is presented in Approximate Nearest Neighbors and the Fast Johnson-Lindenstrauss Transform (slides (ppt) are here). The summary can be described in the Fast Johnson-Lindenstrauss Transform Lemma:

Fix any set X of n vectors in R^d, epsilon < 1, and p is either 1 or 2. Let phi be the FJLT. With probability at least 2/3, we have for all x element of X: (1 + epsilon) * alpha_p * ||x||_2 < || Phi x||_p < (1 + epsilon) * alpha_p * ||x||_2 This dimension reduction technique improves a classic algorithm by Johnson and Lindenstrauss from the mid 80’s. Their technique is the first to offer an asymptotic improvement, and has already been used in design of efficient algorithms for nearest neighbor searching and high dimensional linear algebraic numerical computations. A better technique has recently surfaced as well. A good description of this new technique can be found in this presentation by Edo Liberty.

I have implemented a version of it by following the instructions on the paper but I am not convinced I have done the right implementation as the results are not so sharp to me. I haven't had much time to debug it so here it is, if you find a mistake please let me know:

clear
% Fast Johnson-Lindenstrauss Transform
% as described by Nir Ailon and Bernard Chazelle
% in their paper:
% "Approximate Nearest Neighbors and the Fast Johnson-Lindenstrauss
% Transform"
% The algorithm uses the PHD transform of paragraph 2
% The transform is a k x d projection
% for every set S of n points in R^d
% Let us define an example case
% All columns of G are sets of points
% i.e. we have 1000 points of dimension 128
n = 1000;
d = 256;
% let us examine a reduction from d = 256 to k = 150
%
k = 150;
G = 180 * rand(d,n);
% k = O(epsilon^(-2) * log(n))
epsilon = sqrt(log(n)/k)
%
% epsilon is really O( sqrt(log(n)/k) )
%
%
% Projection in dim k << d % let us take d = 2^w w = log2(d); % % % Defining P (k x d) % first define q % q = min((log(n))^2/d,1); % P = rand(k,d); P = (P < q); P1 = 1/q * randn(k,d); P = P1 .* P; % % Defining H (d x d) H = [1]; for i=1:w H = [ H H;H -H]; end H = 1/d^(0.5)* H; % % Defining D ( d x d) vec1=rand(d,1); v1 = ((vec1 > 0.5) - 1*(vec1<=0.5)); D = diag(v1); % % % FJLT = P * H * D; u = FJLT * G(:,5); v = FJLT * G(:,36); norm(G(:,5)-G(:,36))^2 * k * (1-epsilon) norm( u - v)^2 norm(G(:,5)-G(:,36))^2 * k * (1+epsilon)

Thursday, December 06, 2007

Compressed Sensing: Neighborly polytopes, Acoustic Tomography, Sampling in union of linear subspaces, ROMP stability

Suresh Venkatasubramanian posts about polytopes, compressed sensing, and lower bounds for convex hulls in Neighborly polytopes and compressed sensing. He has an interesting wiki from his seminar course.

Recently, I mentioned a paper on Efficient and Stable Acoustic Tomography Using Sparse Reconstruction Methods by Ivan Jovanović, Ali Hormati, Luciano Sbaiz and Martin Vetterli. It is now here. They solve an inverse problem finding a temperature profile by an inference from the speed of sound. In effect they solve an integral equation like Richard Baraniuk and Philippe Steeghs in Compressive radar imaging and Lawrence Carin, Dehong Liu and Ya Xue in In Situ Compressive Sensing

From the Rice repository, two new items:

Thomas Blumensath and Mike Davies, Sampling theorems for signals from the union of linear subspaces.
The abstract reads
Compressed sensing is an emerging signal acquisition technique that enables signals to be sampled well below the Nyquist rate, given that the signal has a sparse representation in an orthonormal basis. In fact, sparsity in an orthonormal basis is only one possible signal model that allows for sampling strategies below the Nyquist rate. In this paper we consider a more general signal model and assume signals that live on or close to the union of linear subspaces of low dimension. We present sampling theorems for this model that are in the same spirit as the Nyquist-Shannon sampling theorem in that they connect the number of required samples to certain model parameters. Contrary to the Nyquist-Shannon sampling theorem, which gives a necessary and sufficient condition for the number of required samples as well as a simple linear algorithm for signal reconstruction, the model studied here is more complex. We therefore concentrate on two aspects of the signal model, the existence of one to one maps to lower dimensional observation spaces and the smoothness of the inverse map. We show that almost all linear maps are one to one when the observation space is at least twice the dimension of the largest subspace in the union of subspaces. However, we also show that in order for the inverse map to have certain smoothness properties such as a given finite Lipschitz constant, the required observation dimension necessarily depends on the number of subspaces in the signal model. These results are then applied to two examples, the standard compressed sensing signal model in which the signal has a sparse representation in an orthonormal basis and to a sparse signal model with additional tree structure.


and Nam H. Nguyen and Trac D. Tran, The stability of regularized orthogonal matching pursuit. The abstract reads:
This paper studies a fundamental problem that arises in sparse representation and compressed sensing community: can greedy algorithms give us a stable recovery from incomplete and contaminated observations ? Using the Regularized Orthogonal Matching Pursuit (ROMP) algorithm, a modified version of Orthogonal Matching Pursuit (OMP) [1], which was recently introduced by Deanna Needell and Roman Vershynin [ref: Signal recovery from incomplete and inaccurate measurements via Regularized Orthogonal Matching Pursuit ], we assert that ROMP is stable and guarantees approximate recovery of non-sparse signals, as good as the Basis Pursuit algorithm [16]. We also will set up criterions at which the algorithm halts, and the upper bounds for the reconstructed error. It will be proved in the paper that these upper bounds are proportional to the noise energy.

Credit Photo: NASA/JPL/Space Science Institute, Saturn moon Tethys taken on October 29, 2007.

Wednesday, December 05, 2007

Compressed Sensing: Solving the Basis Pursuit problem with a Bregman Iterative Algorithm

While TwIST is solving (IST) compressed sensing reconstruction problems with a two-step Iterative shrinkage/thresholding (IST) algorithm for different regularization terms. Wotao Yin, Stanley Osher, Donald Goldfarb and Jerome Darbon are using Bregman iterative regularization to attack this Basis Pursuit problem in Bregman Iterative Algorithms for Compressed Sensing and Related Problems. The abstract reads:
We propose simple and extremely efficient methods for solving the Basis Pursuit problem



which is used in compressed sensing. Our methods are based on Bregman iterative regularization and they give a very accurate solution after solving only a very small number of instances of the unconstrained problem



for given matrix A and vector fk. We show analytically that this iterative approach yields exact solutions in a finite number of steps, and present numerical results that demonstrate that as few as two to six iterations are sufficient in most cases. Our approach is especially useful for many compressed sensing applications where matrix-vector operations involving A and A' can be computed by fast transforms. Utilizing a fast fixed-point continuation solver that is solely based on such operations for solving the above unconstrained sub-problem, we were able to solve huge instances of compressed sensing problems quickly on a standard PC.

They also seem to aim for very large reconstruction problems as well. The webpage for this effort is here and the code is here. The main difference would certainly be in how much time is spent on both algorithms to solve the same problem but also the fact that the TwIST would have the advantage of also dealing with non-convex regularization lp terms with p =0.7 or 0.8. If anybody has done the comparison between the two methods with p=1, I'd love to hear about it.

Source: Rice Compressed Sensing page.

Compressed Sensing: A new TwIST

After introducing us to the very fast GPSR (new version 5.0), Mário Figueiredo is at it again. José Bioucas-Dias and Mário Figueiredo introduces us to TwIST (Two-step Iterative Shrinkage/Thresholding Algorithm for Linear Inverse Problems) that seems to be both speedy and capable of handling very large unknown vectors found when dealing with images. The main page has examples that include compressed sensing reconstruction examples.


References that explain the algorithms are A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration and Two-step algorithms for linear inverse problems with non-quadratic regularization" both by José Bioucas-Dias, Mário Figueiredo. MATLAB code is available here.

Source: Rice Compressed Sensing page.

Tuesday, December 04, 2007

Compressed Sensing: Random Undersampling in Geophysics and ROMP, a stable greedy algorithm

Gilles Hennenfent and Felix Herrmann just came out with a new paper on the use of Compressed Sensing and Geophysics entitled: Simply denoise: wavefield reconstruction via jittered undersampling. The abstract reads:


In this paper, we present a new discrete undersampling scheme designed to favor wavefield reconstruction by sparsity-promoting inversion with transform elements that are localized in the Fourier domain. Our work is motivated by empirical observations in the seismic community, corroborated by recent results from compressive sampling, which indicate favorable (wavefield) reconstructions from random as opposed to regular undersampling. As predicted by theory, random undersampling renders coherent aliases into harmless incoherent random noise, effectively turning the interpolation problem into a much simpler denoising problem.

A practical requirement of wavefield reconstruction with localized sparsifying transforms is the control on the maximum gap size. Unfortunately, random undersampling does not provide such a control and the main purpose of this paper is to introduce a sampling scheme, coined jittered undersampling, that shares the benefits of random sampling, while offering control on the maximum gap size. Our contribution of jittered sub-Nyquist sampling proofs to be key in the formulation of a versatile wavefield sparsity-promoting recovery scheme that follows the principles of compressive sampling.

After studying the behavior of the jittered-undersampling scheme in the Fourier domain, its performance is studied for curvelet recovery by sparsity-promoting inversion (CRSI). Our findings on synthetic and real seismic data indicate an improvement of several decibels over recovery from regularly-undersampled data for the same amount of data collected.

Deanna Needell and Roman Vershynin just made available a paper entitled:Signal recovery from incomplete and inaccurate measurements via Regularized Orthogonal Matching Pursuit. Roman Vershynin writes the following commentary on it:


This is our second paper with my student Deanna Needell, where we continue to study alternatives to linear programming methods in Compressed Sensing. In our first paper, we developed a greedy algorihtm called ROMP, which can perfectly recover any sparse signal from a small number of linear non-adaptive measurements. In this paper, we prove the stability of ROMP. It seems to be the first greedy algorithm in Compressed Sensing that is provably stable in the sense that the recovery error is proportional to the level of the noise. In particular, the recovery is exact in the absence of noise. Also, one can use ROMP to accurately recover signals, not necessarily sparse. This kind of stability was known for the linear programming methods (a result of Candes, Romberg and Tao), but has been previously unknown for any greedy method (like OMP).

ROMP is available here. Let us note that like the reweighted Lp algorithm of Rick Chartrand and Wotao Yin, reconstruction is performed by a series of iterations solving a Least Squares problem.

Sunday, December 02, 2007

Monday Morning Algorithm part 5: The 1-D Experimental Probabilistic Hypersurface


The Experimental Probabilistic Hypersurface is a tool initially designed by Bernard Beauzamy. Let me first try to say why it is useful first. Imagine you have a function of several variables, say 100, and it takes 2 hours to evaluate that function. At the end of his/her 8 hours day, the engineer will have processed at best 4 function evaluations. It will take a long time of doing any type of sensitivity studies.

In order to have a faster turn around, one could foresee a scheme where one would do a linear interpolation between points (of dimension 100) that we already know. This would provide an idea of what the function looks like in a larger region. Unfortunately, this interpolation scheme in dimension 100 would by itself take too much time. Irrespective to this time constraint, there is nothing that says this function should be linear or even continuous.

In 2004, Bernard Beauzamy made this original construction whereby the function to be evaluated is a multidimensional probabilistic distribution. For every set of (100) parameters that have been already evaluated, the result of the EPH for this function is a dirac. For all other sets of parameters that have not been evaluated, the result of the EPH for this function is a distribution whose shape will be as peaky (resp. flat) as one is close (resp. far) from the sets of parameters for which the function has already been evaluated. Olga Zeydina, a Ph.D candidate, is writing her thesis on the subject. Most information on this construction can be found in the document section of the Robust Mathematical Modeling program at SCM. Of particular interest are: the initial practical implementation [ Report 1 ] and a newer implementation [Report 2].

Today's algorithm will be featured in two files. One that implements the EPH, the other using the EPH function and produce graphs from this document. Listed in this entry, I will list the EPH function file. The second file can be found in the MMA repository.

function [t,pxx]=EPH(xi,theta,x,t1,tau,xe)
% Algorithm developed by Bernard Beauzamy
% and Olga Zeydina
%
% This implementation was written by
% Igor Carron
%
% Inputs
% xi: a parameter set where the function
% has been evaluated
%
% theta: the result of the
% function evaluation for parameter set xi
%
% x: bounding value for parameter set xi
%
% t1: bounding value for result of
% function evaluation
%
% tau: is the increment in the interval
% defined by the bounding value of t1
%
%
% xe: Query parameter set. We want to
% know the value of the function at
% this point.
%
% Outputs
%
% t: interval where the distribution is
% defined
%
% pxx: probability distribution sought
% sampled over t.
%
N =size(xi,2)
%
x1_min= x(1,1)
x1_max= x(1,2)
%
t_min= t1(1,1);
t_max= t1(1,2);
%
%
nu=size([t_min:tau:t_max],2)-1
t=[t_min:tau:t_max];
%
dmax(1)=max(abs(x(1,1)-xi(1)),abs(x(1,2)-xi(1)));
lambd(1)=log(nu)/dmax(1);
for i=2:N
dmax(i)=max(abs(x(1,1)-xi(i)),abs(x(1,2)-xi(i)));
if dmax(i)>dmax(i-1)
lambd(i)=log(nu)/dmax(i);
else
lambd(i)=lambd(i-1);
end
end
nxe=size(xe,2);
options = optimset('TolFun',1e-3);
for jxx=1:nxe
for i=1:N
d(i)=abs(xe(jxx)-xi(i));
end
for i=1:N
a =pi/tau.^2*exp(1-2*lambd(i)*d(i));
lamn=lambd(i)*d(i);
at2=sum(exp(-a.*t.^2));
at3=sum(t.^2.*exp(-a.*t.^2));
diffe=abs(at2-1/tau*(pi/a).^(1/2))/(1/tau*(pi/a).^(1/2));
if abs(diffe)<1e-10
areal(i) = a;
else
areal(i) = fzero(@(x) fina(x,t,lamn),a,options);
end
b = 1/2-lambd(i)*d(i);
breal(i) = log(1/(sum(exp(-areal(i).*t.^2))));
aa(i)=areal(i);
bb(i)=breal(i);
end
for i=1:N
pro(i,:)=d;
end
prd=cumprod(d,2);
for i=1:N
pro(i,i)=1.0;
end
xpro=cumprod(pro,2);
pioverdi=xpro(:,N);
d11s=sum(pioverdi);
for j=1:nu+1
for i=1:N
pij1(i,j)=exp(-aa(i)*(t(j)-theta(i)).^2 + bb(i));
end
pxx(jxx,j)=1/d11s*(pioverdi'*pij1(:,j));
end
end

To run this function, one needs to run the MMA5.m file. This is a simple 1-D implementation. One can obviously foresee its application to 2-D or n-D functions. But more important, one could use the resident expert knowledge of the engineer and use the algorithm of last week [1] to produce a new metric in the parameter space and then use the EPH. If you are using it for some serious work, you may want to join us in the Robust Mathematical Modeling program. This probabilistic approach to function evaluation can also be used as a storage mechanism. Another issue that most engineers face is that after so many function evaluations, it generally is difficult to summarize all the findings. This approach solves this problem. One can even implement a way for the engineer to add parameters as he/she goes.


Compressed Sensing: Restricted Isometry Properties and Nonconvex Compressive Sensing

Rick Chartrand and Valentina Staneva just submitted a paper on the subject of "Restricted isometry properties and nonconvex compressive sensing"

In previous work, numerical experiments showed that lp minimization with p above 0 and p less than 1 recovers sparse signals from fewer linear measurements than does l1 minimization. It was also shown that a weaker restricted isometry property is sufficient to guarantee perfect recovery in the lp case. In this work, we generalize this result to an lp variant of the restricted isometry property, and then determine how many random, Gaussian measurements are sufficient for the condition to hold with high probability. The resulting sufficient condition is met by fewer measurements for smaller p.
In order to prove that a weaker restricted isometry property is sufficient to guarantee perfect recovery in the lp case, they used an improvement of a proof of Dvoretzky's theorem presented in Gilles Pisier 's book. The algorithm exploiting this weaker RIP was featured in Iteratively Reweighted Algorithms for Compressive Sensing and was implemented in the Monday Morning Algorithm Part deux: Reweighted Lp for Non-convex Compressed Sensing Reconstruction. The paper is available from the authors and will be on the LANL website in about a week. [ Update: it is now available here ]

Besides the impressive result on the weaker RIP, one of the experimentally noticeable item featured in the article is that below p=0.8, using the reweighted Lp algorithm, perfect recovery below a certain number of measurements is not possible. That number 0.8 is also the coefficient used for efficient recovery of natural images, in Deconvolution using natural image priors by Anat Levin, Rob Fergus, Fredo Durand, and Bill Freeman on the reconstruction from coded aperture measurements. Naaaaah... it's just a coincidence....

Reference:
[1] Pisier Gilles, The volume of convex bodies and Banach space geometry.

Printfriendly