Monday, November 05, 2007

Monday Morning Algorithm Part 1: Fast Low Rank Approximation using Random Projections

I am not sure I will be able to do this often, but it is always interesting to start the week with an interesting algorithm. So I decided to implement one every beginning of the week. Generally, these algorithms are listed in papers and sometimes when one is lucky they are implemented in Matlab, Octave or equivalent. The beauty of these is certainly underscored when they are readily implemented and you can wow yourself by using them right on the spot. I will try to make it so that one can cut and paste the rest of the entry and run it directly on their Matlab or Octave implementation.
I will generally try to focus on small algorithms (not more than 10 important lines) relevant to some of the issues tackled and mentioned in this blog namely Compressed Sensing, randomized algorithms, bayesian computations, dimensionality reduction....
This week, the script I wrote in matlab is relevant to finding a low-rank approximation to a matrix with more rows than columns (overdetermined systems are like that). One could readily use the SVD command in Matlab/Octave but it can take a long time for large matrices. This algorithm uses properties from Random matrices to do the bidding in less time than the traditional SVD based method. I initially found it in [3], the original reference is [1] and it is used in the context of LSI is [2].

clear
% preparing the problem
% trying to find a low approximation to A, an m x n matrix
% where m >= n
m = 1000;
n = 900;
% first let's produce example A
A = rand(m,n);
%
% beginning of the algorithm designed to find alow rank matrix of A
% let us define that rank to be equal to k
k = 50;
% R is an m x l matrix drawn from a N(0,1)
% where l is such that l > c log(n)/ epsilon^2
%
l = 100;
% timing the random algorithm
trand =cputime;
R = randn(m,l);
B = 1/sqrt(l)* R' * A;
[a,s,b]=svd(B);
Ak = A*b(:,1:k)*b(:,1:k)';
trandend = cputime-trand;
% now timing the normal SVD algorithm
tsvd = cputime;
% doing it the normal SVD way
[U,S,V] = svd(A,0);
Aksvd= U(1:m,1:k)*S(1:k,1:k)*V(1:n,1:k)';
tsvdend = cputime -tsvd;
%
%
% relative error between the two computations in percent
rel = norm(Ak-Aksvd)/norm(Aksvd)*100
% gain in time
gain_in_time = tsvdend - trandend
% the random algorithm is faster by
tsvdend/trandend


If you find any mistake, please let me know.

References:
[1] The Random Projection Method, Santosh S. Vempala
[2] Latent Semantic Indexing: A Probabilistic Analysis (1998) Christos H. Papadimitriou, Prabhakar Raghavan, Hisao Tamaki, Santosh Vempala
[3] The Random Projection Method; chosen chapters from DIMACS vol.65 by Santosh S. Vempala, presentation by Edo Liberty October 13, 2006
Photo Credit: ESA/ASI/NASA/Univ. of Rome/JPL/Smithsonian, This image shows the topographic divide between the Martian highlands and lowlands. The mysterious deposits of the Medusae Fossae Formation are found in the lowlands along the divide. The radar sounder on ESA's Mars Express orbiter, MARSIS, has revealed echoes from the lowland plains buried by these mysterious deposits.

Sunday, November 04, 2007

Compressed Sensing: High Resolution Radar, Identification of Sparse Matrices


Sometimes I get to find preprints that are not yet on the Rice Resource page. Here are some:

Matthew Herman and Thomas Strohmer just came out with a preprint on High Resolution Radar via Compressed Sensing. The abstract reads:
A stylized compressed sensing radar is proposed in which the time-frequency plane is discretized into an N×N grid. Assuming the number of targets K is small (i.e., K  N2), then we can transmit a sufficiently “incoherent” pulse and employ the techniques of compressed sensing to reconstruct the target scene. A theoretical upper bound on the sparsity K is presented. Numerical simulations verify that even better performance can be achieved in practice. This novel compressed sensing approach offers great potential for better resolution over classical radar.
As for any new technique, when being compared with techniques that have been used for decades, it is hard to make a breakthrough. Here this difficulty seems to lie more in the sophistication of the target matching currently used than the detection per se. For the detection by this "stylised" approach, one seems to have an enormous gain in the resolution. The compressed sensing radar has a resolution of 1/(2 sqrt( N)) whereas a similar classical radar has a resolution of 1/2.


Identification of Matrices having a Sparse Representation by Jared Tanner, Gotz E. Pfander and Holger Rauhut. The abstract reads:
We consider the problem of recovering a matrix from its action on a known vector in the setting where the matrix can be represented efficiently in a known matrix dictionary. Connections with sparse signal recovery allows for the use of efficient reconstruction techniques such as Basis Pursuit (BP). Of particular interest is the dictionary of time–frequency shift matrices and its role for channel estimation and identification in communications engineering. We present recovery results for BP with the time–frequency shift dictionary and various dictionaries of random matrices.

Saturday, November 03, 2007

DARPA Urban Challenge: It's today

We did not qualify for the semi-finals but others did. DARPA has finally selected the teams that will compete in the Urban Challenge Final Event on Saturday, November 3 at the former George Air Force Base in Victorville, Calif. The teams will attempt to complete a complex 60-mile urban course with live traffic in less than six hours. The finalists will operate on the course roads with approximately 50 human-driven traffic vehicles. Speed is not the only factor in determining the winners, as vehicles must also meet the same standards required to pass the California DMV road test. The Urban Challenge Final Event is to take place today (November 3) and will be webcast live at http://www.grandchallenge.org starting at 7:30 am PT. Grounds will open at 6:00 AM PT for spectators, but the opening ceremony will begin at 7:30 AM, and vehicles will begin to launch at 8:00 AM.

Watch the Race here, Status Board is here, Live Video is here, Tracking Map is here.

A brief highlight video of National Qualifying Event (semi-final) operations is available for download here: High (27MB) Med (7.5MB) Low (3MB)

Vehicles are being tested in three test areas to evaluate their ability to operate with live traffic, make safe left turns across moving traffic, and pull out at T-intersections with cars arriving from both directions. Vehicles also have to follow narrow winding roads, avoiding parked cars and other obstructions, maneuver into a designated parking spot and negotiate 4-way intersections and road-blocks.

List of teams that did qualify for the Final Event can be found here.
List of teams that did not qualify for the Final Event can be found here.

I personally will not watch, I am still sore from this.

Friday, November 02, 2007

Compressed Sensing: Toeplitz, Radar, SAR, Bayesian Approach, Non Convex Minimization, Angle Preserving RIP




There is a new batch of papers showing up on the Compressed Sensing Repository at Rice.

* Toeplitz-structured compressed sensing matrices by Waheed Bajwa, Jarvis Haupt, Gil Raz, Stephen Wright, and Robert Nowak. In the current Compressed Sensing Literature, random matrices used to produce CS measurements are made by sampling each element out of a specific probability distribution. This paper shows that one can build Toeplitz-structured matrices with entries drawn independently from probability distributions that yield good CS measurement matrices. The abstract reads:
The problem of recovering a sparse signal x ∈ Rn from a relatively small number of its observations of the form y = Ax ∈ Rk, where A is a known matrix and k ≪ n, has recently received a lot of attention under the rubric of compressed sensing (CS) and has applications in many areas of signal processing such as data compression, image processing, dimensionality reduction, etc. Recent work has established that if A is a random matrix with entries drawn independently from certain probability distributions then exact recovery of x from these observations can be guaranteed with high probability. In this paper, we show that Toeplitz-structured matrices with entries drawn independently from the same distributions are also sufficient to recover x from y with high probability, and we compare the performance of such matrices with that of fully independent and identically distributed ones. The use of Toeplitz matrices in CS applications has several potential advantages: (i) they require the generation of only O(n) independent random variables; (ii) multiplication with Toeplitz matrices can be efficiently implemented using fast Fourier transform, resulting in faster acquisition and reconstruction algorithms; and (iii) Toeplitz-structured matrices arise naturally in certain application areas such as system identification.
I wonder how building these matrices would fit into the work of Benjamin Recht, Maryam Fazel, Pablo Parrilo on minimum rank approximation ( Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization )

* In Situ Compressive Sensing (short version) by Lawrence Carin, Dehong Liu, and Ya Xue. The abstract reads:
It has been demonstrated that the fields observed at a point in a complex propagation medium may be viewed as a linear combination the fields from the same source, as radiated into a homogeneous environment. In the context of compressive sensing (CS), the former fields represent the data v and the latter u. Exploiting the compressible nature of most fields, as a function of position and frequency, one may use CS to infer u from v, with the number of required measurements for v small relative to those typically required for u. This has been referred to as in situ CS, because the complex propagation medium is itself exploited to perform the compressive measurements. The basic concept has been demonstrated for the problem of a source radiating in the presence of parallel plates, with encouraging results demonstrated. The linkages between CS and the superresolution properties of time-reversal [8] have been made, and it has been informally asserted that media that are appropriate for timereversal will likely also be appropriate for in situ CS. There is significant fundamental research required to formalize the requirements on the propagation media to optimize the performance of in situ CS.
* Exact reconstructions of sparse signals via nonconvex minimization. by Rick Chartrand is a useful complement to Iteratively reweighted algorithms for compressive sensing by Rick Chartrand and Wotao Yin, that was the subject of a previous entry. The abstract reads:

Several authors have shown recently that is possible to reconstruct exactly a sparse signal from fewer linear measurements than would be expected from traditional sampling theory. The methods used involve computing the signal of minimum l1 norm among those having the given measurements. We show that by replacing the l1 norm with the lp norm with p less than 1, exact reconstruction is possible with substantially fewer measurements. We give a theorem in this direction, and many numerical examples, both in one complex dimension, and larger scale examples in two real dimensions.

* There is also: Stable sparse approximations via nonconvex optimization by Rayan Saab, Rick Chartrand and Özgür Yilmaz.

We present theoretical results pertaining to the ability of lp minimization to recover sparse and compressible signals from incomplete and noisy measurements. In particular, we extend the results of Candes, Romberg and Tao [1] to the p less than 1 case. Our results indicate that depending on the restricted isometry constants (see, e.g., [2] and [3]) and the noise level, lp minimization with certain values of p less than 1 provides better theoretical guarantees in terms of stability and robustness than l1 minimization does. This is especially true when the restricted isometry constants are relatively large.
The paper takes on the stable recovery issue of the signal in the presence of noise with the Lp approach.

*Compressive sampling for signal detection by Jarvis Haupt and Robert Nowak.

Compressive sampling (CS) refers to a generalized sampling paradigm in which observations are inner products between an unknown signal vector and user-specified test vectors. Among the attractive features of CS is the ability to reconstruct any sparse (or nearly sparse) signal from a relatively small number of samples, even when the observations are corrupted by additive noise. However, the potential of CS in other signal processing applications is still not fully known. This paper examines the performance of CS for the problem of signal detection. A generalized restricted isometry property (GRIP) is introduced, which guarantees that angles are preserved, in addition to the usual norm preservation, by CS. The GRIP is leveraged to derive error bounds for a CS matched filtering scheme, and to show that the scheme is robust to signal mismatch.
This approach introduces GRIP which is a restricted isometric property that also preserves angles. I wonder how it will influence results in manifold learning. I can think of many different problems where this is of interest.

* Detecting signal structure from randomly-sampled data by Frank Boyle, Gerald Fudge, Jarvis Haupt and Robert Nowak. The abstract reads:
Recent theoretical results in Compressive Sensing (CS) show that sparse (or compressible) signals can be accurately reconstructed from a reduced set of linear measurements in the form of projections onto random vectors. The associated reconstruction consists of a nonlinear optimization that requires knowledge of the actual projection vectors. This work demonstrates that random time samples of a data stream could be used to identify certain signal features, even when no time reference is available. Since random sampling suppresses aliasing, a small (sub-Nyquist) set of samples can represent high-bandwidth signals. Simulations were carried out to explore the utility of such a procedure for detecting and classifying signals of interest.
* Bayesian compressive sensing by Shihao Ji, Ya Xue, and Lawrence Carin (I covered this one here) and Bayesian Compressive Sensing and Projection Optimization by Shihao Ji, Lawrence Carin,

The abstract reads:
This paper introduces a new problem for which machine-learning tools may make an impact. The problem considered is termed “compressive sensing”, in which a real signal of dimension N is measured accurately based on KN real measurements. This is achieved under the assumption that the underlying signal has a sparse representation in some basis (e.g., wavelets). In this paper we demonstrate how techniques developed in machine learning, specifically sparse Bayesian regression and active learning, may be leveraged to this new problem. We also point out future research directions in compressive sensing of interest to the machine learning community.
This approach is indeed very interesting as they seem to have faster convergence in the reconstruction stage than the best greedy algorithm. They also claim to have sparser results as well.


* In a previous entry, we were shown that one could do a better job with radars. The paper is now out:
Compressive radar imaging by Richard Baraniuk and Philippe Steeghs.

The abstract reads:
We introduce a new approach to radar imaging based on the concept of compressive sensing (CS). In CS, a low-dimensional, nonadaptive, linear projection is used to acquire an efficient representation of a compressible signal directly using just a few measurements. The signal is then reconstructed by solving an inverse problem either through a linear program or a greedy pursuit. We demonstrate that CS has the potential to make two significant improvements to radar systems: (i) eliminating the need for the pulse compression matched filter at the receiver, and (ii) reducing the required receiver analog-to digital conversion bandwidth so that it need operate only at the radar reflectivity’s potentially low “information rate” rather than at its potentially high Nyquist rate. These ideas could enable the design of new, simplified radar systems, shifting the emphasis from expensive receiver hardware to smart signal recovery algorithms.

This is the other paper, along with the work of Carin (see above on In-situ Compressed Sensing by Lawrence Carin, Dehong Liu, and Ya Xue) that solves an integral equation using strictly CS. It is aimed at reducing the amount of electronics that goes into the processing stage of radars in order to be able to deal with objects of different complexity.

* Fast encoding of synthetic aperture radar raw data using compressed sensing by Sujit Bhattacharya, Thomas Blumensath, Bernard Mulgrew, and Mike Davies.

The abstract reads:
Synthetic Aperture Radar (SAR) is active and coherent microwave high resolution imaging system, which has the capability to image in all weather and day-night conditions. SAR transmits chirp signals and the received echoes are sampled into In-phase (I) and Quadrature (Q) components, generally referred to as raw SAR data. The various modes of SAR coupled with the high resolution and wide swath requirements result in a huge amount of data, which will easily exceed the on-board storage and downlink bandwidth of a satellite. This paper addresses the compression of the raw SAR data by sampling the signal below Nyquist rate using ideas from Compressed Sensing (CS). Due to the low computational resources available onboard satellite, the idea is to use a simple encoder, with a 2D FFT and a random sampler. Decoding is then based on convex optimization or uses greedy algorithms such as Orthogonal Matching Pursuit (OMP).

It looks like they are using OMP because they deal with images that are inherently large. Since current reconstruction l1 reconstruction capabilities are limited one has to resort to greedy algorithms. One wonders how the Lp-L2 of Chartrand (see above) would fare....


Reference: Rice Compressed Sensing Repository
Image Credit: NASA/JPL/Space Science Institute, Rings of Saturn taken on June 12, 2007 by the Cassini probe.

Monday, October 29, 2007

Compressed Sensing: Sparco helps you test all these spiffy reconstruction algorithms


With the large number of reconstruction techniques, it was just a question of time before somebody would produce a toolbox for unifying the testing of these sparse reconstruction algorithms. This is what Ewout van den Berg, Michael P. Friedlander, Gilles Hennenfent, Felix J. Herrmann, Rayan Saab, and Ozgur Yilmaz have done by releasing Sparco.

Sparco: A Testing Framework for Sparse Reconstruction

Sparco is a framework for testing and benchmarking algorithms for sparse reconstruction. It includes a large collection of sparse reconstruction problems drawn from the imaging, compressed sensing, and geophysics literature. Sparco is also a framework for implementing new test problems and can be used as a tool for reproducible research. Sparco is implemented entirely in Matlab, and is released as open-source software under the GNU Public License.
The technical report is here.

Making sense of two findings

Two wows this past week:

First, Scott Parazynski and Daniel Tani found that there are metal shavings in one of the joint of the rotating solar arrays of the International Space Station (ISS). Since you can't change them, this is major news as it may critically impact the ISS in the future.

Second, the news that Uranium isotope ratios have been shown to not be invariant. This is the first time since the Oklo natural reactor that we show that natural reactors have occurred in many other places (though not at the rate of the Oklo). This is significant because it shows that deep burying of current nuclear waste is not unlike what is currently done in nature. Furthermore, it maybe another element showing us that Earth is a giant nuclear reactor. In this latest theory, I have always been amazed at the Origen-S and SAS2 jobs performed to compare to results found in nature. These codes are generally used in nuclear engineering to determine parameters of interest in nuclear reactor cores.

Credit: NASA, Earth-Moon system as viewed from a Mars probe on its way to Mars. 2003.

Saturday, October 27, 2007

Compressed Sensing: RandomFaces prove that faces lie in a low dimensional manifold

It was only a matter of time, the Machine Learning community is finding ways to use compressed sensing following the work of Richard Baraniuk and Michael Wakin on Random projections of smooth manifolds.
There were a laplacianfaces, eigenfaces, we now have randomfaces, the projection of faces onto random bases. The work of Yi Ma ([1] [2] [3])seems to show that one can get similar results in terms of classification as those obtain from the more complex techniques producing the Laplacianfaces and other Eigenfaces. Some people would say that it may parallel what is happening in the primary visual cortex.

[ Des presentations sur le compressed sensing (acquisition comprimee) se feront a l'Institut Henri Poincaré (Paris, France) en Mars 2008, Méthodes variationnelles et parcimonieuses en traitement des signaux et des images , Merci Laurent]

References:
[1] Robust Face Recognition via Sparse Representation, John Wright, Arvind Ganesh, Allen Yang and Yi Ma. Submitted to IEEE Transactions on Pattern Analysis and Machine Intelligence. [Updated January 2008]
[2] Feature Selection in Face Recognition: A Sparse Representation Perspective, Allen Yang, John Wright, Yi Ma and Shankar Sastry. Submitted to IEEE Transactions on Pattern Analysis and Machine Intelligence.
[3] Computation and Relaxation of Conditions for Equivalence between L1 and L0 Minimization, Yoav Sharon, John Wright and Yi Ma. Submitted to IEEE Transactions on Information Theory.

Compressed Sensing: Reweighted Lp through Least square (L2), Forget Europa, let's shoot for Titan


In a previous entry, I mentioned the impressive improvement of the L1 reconstruction capability for signals that were not so sparse. We also knew that an Lp (p less than 1) exact reconstruction of sparse signals via nonconvex minimization was also possible.

Based on the work of Bhaskar Rao and Kenneth Kreutz-Delgado [2], Rick Chartrand and Wotao Yin produced an astounding reconstruction technique in this article [2]:

ITERATIVELY REWEIGHTED ALGORITHMS FOR COMPRESSIVE SENSING

The theory of compressive sensing has shown that sparse signals can be reconstructed exactly from remarkably few measurements. In [1], it was shown empirically that using Lp minimization with p less than 1 can do so with fewer measurements than with p = 1. In this paper we consider the use of iteratively reweighted algorithms for computing local minima of the nonconvex problem. In particular, a particular regularization strategy is found to greatly improve the ability of a reweighted least-squares algorithm to recover sparse signals, with exact recovery being observed for signals that are much less sparse than required by an unregularized version (such as FOCUSS, [2]). Improvements are also observed for the reweighted-l1 approach of [3].


Why is this is an important paper ? Because there is
  • no need to use SDP programming, the scheme is very simple computationally. It can be written in ten lines of Matlab (Update: This algorithm was featured in this Monday Morning Algorithm)
  • everybody in engineering knows least squares
  • some of us now have a stab at the very large problems.

The major feature of this paper is showing that $\epsilon$ should probably be not that small in the first place. Let us hope that this heuristic will be proven in the future. In order to understand the graph above, one needs to know the acronyms:

  • IRLS stands for Iteratively Reweighted or weighted Least Square (the algorithm of this paper)
  • IRL1 stands for Iteratively Reweighted L1 (the case of p=0 is the same as the reweighted L1 algorithm in [3]).
K stands for the sparsity number of the signal, M is the number of measurements and N is the signal size.

Please note comparable numbers between IRL1/p=0 (linear programming algorithm in [3]) and IRLS/p=0 (least square algorithm in [1]).

Reference:
[1] ITERATIVELY REWEIGHTED ALGORITHMS FOR COMPRESSIVE SENSING, Rick Chartrand, Wotao Yin, ICASSP 2008

[2] An affine scaling methodology for best basis selection.
Rao, B.D.; Kreutz-Delgado, K. IEEE Transactions on Signal Processing, vol.47, (no.1), IEEE, Jan. 1999. p.187-200. [Rao1999.pdf]

[3] Enhancing sparsity by reweighted L1 minimization, Emmanuel Candès, Mike Wakin and Stephen Boyd

Compressive Sensing Radiation Detector Part II

In the case where radiation does not suffer much scattering, the principle of Coded Aperture used for light can be directly applied to these radiation fields. Berthold Horn, Richard Lanza, Roberto Accorsi,
Klaus Ziock, and Lorenzo Fabris at MIT have been following this principle to do imaging by performing Computational Imaging. As pointed out before, Coded Aperture is also what the DISP group does at Duke University (but with radiation that are in the visible range). The big difference resides in the Mask, in radiation with smaller wavelengths (X-rays,...), radiation goes through the mask, but more generally you cannot guide the radiation to your detector and so it becomes a little bit more complex. Spatial resolution is extremely difficult to do because the detector are expensive and need to be where the radiation goes. So one has to resort to batteries of detectors such as the ones below.
In order to provide spatial resolution, one has to either to have many detectors or move the detector or the source as can be seen in this animation (“Synthetic Aperture” radiography). An introduction to instrumentation using coded aperture for astronomy can be found here. This is a great introduction to the subject, especially in regards to the current inversion techniques.

When one deals with radiation scattering (neutrons in nuclear reactor cores or infra-red/visible light in tissues) , radiation sensors have a much more difficult kernel to invert (because the diffusion/scattering of the rays in matter) and requires methods like Bayesian Compressed Sensing (i.e. In-situ detector of Lawrence Carin, Dehong Liu, and Ya Xue at Duke University). All in all, it looks like some of these techniques for Coded Aperture could be improved at the software and hardware level using Compressed Sensing.

Wednesday, October 24, 2007

New satellite imagery of the California fires

ESA just released new satellite imagery from the California fires using the MERIS camera on ENVISAT. These shots were taken 3 hours ago. IMagery from two days ago can be found in a previous entry.

I made a close up version centered on L.A.

The full imagery is shown below. The cloud at the bottom is clearly different in nature than the ones at the top. i.e. they are the result of the fire.

Printfriendly