Tuesday, January 20, 2009

CS: GPSR 6.0, ILRS Minimization for Sparse Recovery Code and On Sparse Solutions of Underdetermined Linear Systems, Data Streams

The GPSR page was updated showing that now GPSR is at version 6.0: Download it at GPSR_6.0. Thank you Mário Figueiredo, Robert D. Nowak, Stephen J. Wright!

I featured the paper before, but Massimo Fornasier just released the attendant poster and the code featured in Iteratively Re-weighted Least Squares Minimization for Sparse Recovery by Ingrid Daubechies, Ronald DeVore, Massimo Fornasier, and C. Sinan Güntürk. The codes are here:

They are also in the reconstruction section of the big picture in compressive sensing.

Here is also a survey article on l_p areconstruction techniques in On Sparse Solutions of Underdetermined Linear Systems by Ming-Jun Lai. The abstract reads:
We first explain the research problem of finding the sparse solution of underdetermined linear systems with some applications. Then we explain three different approaches how to solve the sparse solution: the ℓ1 approach, the orthogonal greedy approach, and the ℓq approach with 0 less q less 1. We mainly survey recent results and present some new or simplified proofs. In particular, we give a good reason why the orthogonal greedy algorithm converges and why it can be used to find the sparse solution. About the restricted isometry property (RIP) of matrices, we provide an elementary proof to a known result that the probability that the random matrix with iid Gaussian variables possesses the PIP is strictly positive.


The 2009 Annual Workshop on Computational Complexity in Barbados will feature S. Muthu Muthukrishnan who will be the main speaker. The course is titled Data Streams:
Proposed outline

In order to deal with massive data that arrives rapidly and has to be processed instantly, the theoretical computer science community has developed the sublinear space, polylog time-per item, one-pass "data stream" model. In this course, we will study the theory of algorithms in the data stream model. Lectures will include (from here):

* Data stream models and applications
* Forward distributions and Sketches
* Inverse distributions and Samples.
* Advanced Topic: Lower bounds in Data Stream model.
* Advanced Topic: Algorithms for graph, geometric, matrix streams.
* Advanced Topic: Second generation models including multipass, MapReduce and Probablistic streams.
* Applications beyond data streams: compressed sensing, machine learning.

We start with the basic material from the book "Data Streams and Applications", but the advanced material is from the ongoing Stream 2.0 survey with Andrew McGregor.

Background material is a document entitled: Data Streams: Algorithms and Applications.

Sunday, January 18, 2009

CS: Near-Optimal Bayesian Localization via Incoherence and Sparsity, Lp/RIP/Pessimism, Two jobs and a call for paper



Everybody has seen the Youtube video, but for those interested in using it as an example in a paper, the avi video is on the Coast Guards' site. Maybe it could be used as a benchmark of some kind. In the meantime, we have one paper, two job announcements for two postdocs, an internship and a CfP.

Volkan Cevher just sent me this freshly accepted and very interesting paper entitled Near-Optimal Bayesian Localization via Incoherence and Sparsity by himself, Petros Boufounos, Richard Baraniuk, Anna Gilbert, Martin Strauss. The abstract reads:

Source localization using a network of sensors is a classical problem with applications in tracking, habitat monitoring, etc. A solution to this estimation problem must satisfy a number of competing resource constraints, such as estimation accuracy, communication and energy costs, signal sampling requirements and computational complexity. This paper exploits recent developments in sparse approximation and compressive sensing to efficiently perform localization in a sensor network. We introduce a Bayesian framework for the localization problem and provide sparse approximations to its optimal solution. By exploiting the spatial sparsity of the posterior density, we demonstrate that the optimal solution can be computed using fast sparse approximation algorithms. We show that exploiting the signal sparsity can reduce the sensing and computational cost on the sensors, as well as the communication bandwidth. We further illustrate that the sparsity of the source locations can be exploited to decentralize the computation of the source locations and reduce the sensor communications even further. We also discuss how recent results in 1-bit compressive sensing can impact the sensor communications by transmitting only the timing information relevant to the problem. Finally, we develop a computationally efficient algorithm for bearing estimation using a network of sensors with provable guarantees.

I'll have to dwell into it further, but the mixing of sensor network and 1-bit compressive sensing makes it a unique and promising approach.

Remi Gribonval released the slides of his lecture at the Cambridge Workshop on Sparsity and Large-Scale inverse problems, on Lp minimisation and the role of restricted isometry constants. The presentation is entitled: Some stories about Lp-minimization, the Restricted Isometry Property, and excessive pessimism.. I should come back to this later and will probably integrate it in the Big Picture Section.

Remi also made me aware of this first announcement as featured in the Compressive Sensing Jobs listing:

January 16th, 2009: Two Post-Docs at INRIA Rennes, France (IRISA). The postdocs will work on theoretical, algorithmic and practical aspects of sparse representations of large-dimensional data, with a particular emphasis on acoustic fields, for various applications such as compressed sensing, source separation and localization, and signal classification. Previous experience in sparse representations (time-frequency and time-scale transforms, pursuit algorithms, support vector machines and related approaches) is desirable, as well as a strong taste for the mathematical aspects of signal processing. For additional technical information, please contact :Rémi GRIBONVAL - SMALL/ECHANGE project leader - METISS Project-Team - INRIA-Bretagne Atlantique - Email: remi.gribonval@inria.fr - phone: +33 2 99 84 25 06. The positions, funded for at least 2 years (up to three years), will be renewed on a yearly basis depending on scientific progress and achievement. The gross minimum salary will be 28287 € annually (~ 1923 € net per month) and will be adjusted according to experience. The usual funding support of any French institution (medical insurance, etc.) will be provided. More information can be found here

and then I found this one on the internets:

January 16th, 2009: Intern at Microsoft Research Asia. Beijing. MRA is recruiting interns for Stars of Tomorrow Internship Program.. For more information see here.
Position: Full-time Intern
Group: Visual Computing Group
Quantity: 1
Work Location: Beijing. The interns will be working with the VC manager, professor Yi Ma, on applications of sparse representation and compressive sensing to problems in computer vision and pattern recognition (e.g., face recognition, image/video segmentation and analysis). The interns will be required to conduct programming, simulation, and experiment for various research projects, or if capable, conduct mathematical analysis of algorithms and solve new problems. Qualified applicants please fill in the application form and send it together with a full resume in both English and Chinese (PDF/Word/Txt/Html format) to: MSRAih@microsoft.com, and please note you are applying for MSRA VC Group. Know more about Stars of Tomorrow Internship Program, please visit http://www.msra.cn/ur/intern.aspx. If any question, please email to msrastar@microsoft.com.


Finally, here is a Call for Paper in the IEEE Signal Processing Magazine (SPM) Convex Optimization for Signal Processing (SPM 2009)
  • Submission Deadline: May 5, 2009
  • Notification Due: Dec 1, 2009
  • Final Version Due: Jan 15, 2010
link: http://apollo.ee.columbia.edu/spm/?i=cfp/May10

Call For Papers

In recent years, we have witnessed technical breakthroughs in a wide variety of topics where the key to success is the use of convex optimization. In fact, convex optimization has now emerged as a major signal processing technique and has made significant impact on numerous problems previously considered intractable. Today, innovative applications of convex optimization in signal processing range from those in adaptive filtering, detection and estimation, sensor array processing, MIMO communications, sensor networks, sampling theory, and more recently, image processing, speech processing and cognitive radios - and the scope of its applications is still expanding. Considering the foundational nature and potential impact of convex optimization in signal processing, there appears to be a clear need for a special issue that introduces convex optimization to the broad signal processing community, gives insights into how convex optimization can make a difference, and showcases some notable successes. This special issue aims to solicit papers that provide tutorials of convex optimization techniques (including available software) and various successful signal processing applications. Also welcome are tutorial papers that deal with emerging, meaningful applications; or that give friendly overviews of certain theoretically advanced convex optimization techniques relevant to signal processing.

To enhance readability and appeal for a broad signal processing audience, prospective authors are encouraged to use an intuitive approach in their presentation; e.g., by using simple instructive examples, considering special cases that show insights into the ideas, and using illustrations to the extent possible.

Examples of topics that will be addressed in this special issue include, but are not limited to:

* Adaptive filtering
* Beamforming
* Convex optimization fundamentals including relaxation techniques and software toolboxes
* Image processing applications, including denoising, MRI image processing, phase unwrapping
* Compressed sensing
* Detection and estimation
* MIMO communications
* Sensor networks
* Cognitive radios
* Speech applications

Submission Procedure:

Prospective authors should submit their white papers through the web submission system at www.ee.columbia.edu/spm. The white paper should be no more than 6 pages in the IEEE double-space one-column 11-point format.

Schedule (all deadlines are firm no exceptions)
  • White paper due: May 5, 2009
  • Invitation notification: June 1, 2009
  • Manuscript submission: September 1, 2009
  • Notification of acceptance: December 1, 2009
  • Final manuscript decision: January 15, 2010
  • Publication date: May, 2010

Guest Editors:
  • Yonina Eldar, Technion, Israel Institute of Technology, yonina@ee.technion.ac.il
  • Zhi-Quan Luo, University of Minnesota, luozq@umn.edu
  • Wing-Kin (Ken) Ma, The Chinese University of Hong Kong, wkma@ee.cuhk.edu.hk
  • Daniel Palomar, Hong Kong University of Science and Technology, palomar@ust.hk
  • Nikos Sidiropoulos, Technical University of Crete, nikos@telecom.tuc.gr

Friday, January 16, 2009

CS: SparSA, Q/A, CS and Random Arrays, KF for reconstruction, Inferring Rankings Under Constrained Sensing, Post-silicon timing, a course and SSP09

And you thought, you'd go home for the week-end with nothing to stuff your brain with. Tough luck!

First, after a request, Mario Figueiredo kindly released the SparSA code featured in “Sparse reconstruction by separable approximation” (new version) by Stephen Wright, Robert Nowak, Mario Figueiredo. The MATLAB code available here:

http://www.lx.it.pt/~mtf/SpaRSA

It has been added to the reconstruction section of the Big Picture. By the way, I have tried to clean up both the reconstruction and the measurement/encoding matrix section, if you find I missed an item or I am wrong, please let me know. Thanks. In particular, the measurement/encoding matrix section also includes a blurb on the different properties (RIP, JN,....) that a matrix should have to allow for recovery of sparse solutions.

After the CS in radio interferometry entry of last week . I went ahead and asked dumb questions to the authors:
...What I am gathering out of it is that CS-spgl1 and Clean are kind of close to each other except for an advantage for CS for the lower number of measurements....
Laurent Jacques responded with :
This is not exactly what we shown. Let me summarize our approach here and Yves will possibly complement my explanation.

First, CLEAN is a Matching Pursuit with an additional degree of freedom : a parameter gamma that weight the removal of the selected atom from the current residual. (for usual MP, gamma=1 since the whole atom is removed from the residual at each iteration). Since astronmers took this parameter very small (typically gamma=1/10), they got a greedy method very efficient for deconvolving the true sky from the frequency coverage seen as a convolutive kernel. Or, using CS terminology, CLEAN is good to minimize the number of measurements while still having a good reconstruction. Therefore, in its results, CLEAN is closer to OMP than MP. My interpretation is that CLEAN is a "kind of" Taylor development of OMP making CLEAN very close to OMP when gamma is small.

Second, what we observed is that CLEAN and Basis Pursuit BP (solved by any methods like spgl1 or proximal methods) are working (in average) equally well, i.e. they are similar in the number of measurements they required to reach a certain quality of reconstruction. Remember that there is no exact recovery since the signal is compressible and since the measurements can be noisy. In my mind, we tried also with OMP and the quality reached was also very similar for the same number of measurements. This could be reconfirmed.

So, for us, BP does not beat CLEAN for the quality (SNR) of the reconstruction for this particular application where the sensing is performed randomly in the Fourier plane and where the sparsity basis is the canonical basis. BP has however the advantage to see clearly the global minimization to reach. Therefore, the number of iterations in a BP solver is well controlled compared to CLEAN where the stopping criterion is not always clearly set (and a small gamma provides perhaps a better quality but also a huge number of iterations).

I also asked :
... Except that I am under the weird impression that you are not comparing similar techniques since you put an additional constraint of positivity on the CS scheme....

Laurent Jacques answered :
BP+, i.e. BP plus the positivity constraint, is outside of the CS theory in our presentation of CS. (I know that they are some recent progress on that point in the CS theory, e.g working with the null space property.) So here, we have observed that BP+ provides better quality than CLEAN/BP for the same frequency coverage (number of measurements). We do not prove why. We just assume that the addition of this new a priori information(that correspond to a convex constraint in the problem) makes the problem more stable in its resolution.

One could mention that Clean or MEM are not using specifically the sparsity of the signal and on could say that indeed you are not comparing similar techniques. CLEAN assumes the sparsity as MP does. So, we can compare it to BP and BP+.

However, MEM assumes the smoothness/constancy of the signal (the converse of the sparsity) since by definition it maximizes the entropy. So, we didn't compare MEM to CLEAN/BP/BP+ because of this.

Yves Wiaux further added:
...Astronomers use both CLEAN and MEM. In that regard we might want to compare our methods to both. But indeed we only compare BP and our prior-enhanced-BP to CLEAN (which somehow already assumes sparsity, while MEM does not). And for both kinds of signals analyzed (a positive signal without noise, and a general signal with noise and whose statistics were preliminarily studied) simulations show that BP and CLEAN provide more or less the same performance. On the contrary prior-enhanced-BP does significantly better (here the accuracy is proven from simulations only, as CS gives no theoretical result for these enhanced problems, for the moment)...
while Pierre Vandergheynst finally added:
...another advantage of BP or global minimization algos with respect to greedy (or local minimization) is this possibility of adding external constrains (like positivity). You can't do that too easily with MP for example. Same thing for plugging in deeper statistical models : you can more easily weight BP with a priori than MP (though we've done that in a paper 2 years ago and compared both, see: On the Use of A Priori Information for Sparse Signal Approximations...
This is all great, thank you Laurent, Yves, Pierre for being kind enough to answer. As one can see older algorithms will eventually be included in the CS framework and we can expect that information to lead into improved algorithms. Lawrence Carin is making that sentiment crystal clear with On the Relationship Between Compressive Sensing and Random Sensor Arrays. The abstract reads:
Random sensor arrays are examined from a compressive sensing (CS) perspective. It is demonstrated that the natural random-array projections manifested by the media Green’s function are consistent with the projection-type measurements associated with CS. This linkage allows the use of existing CS theory to quantify the performance of random arrays, of interest for array design. The analysis demonstrates that the CS theory is applicable to arrays in vacuum as well as in the presence of a surrounding media; further, the presence of a surrounding media with known properties may be used to improve array performance.

A very nice paper. Echoing the previous Q/A, we can read:
We demonstrate that these algorithms, as well as RELAX from array processing [15], while developed independently, are highly related to one another (in fact, OMP and CLEAN are essentially the same algorithm).
or
This paper makes the explicit connection between decades-old random sensor arrays and the much newer CS, demonstrating that the former is a special case of the latter. It is therefore not surprising that the aforementioned independently developed algorithms are highly inter-related.
and finally,
A challenge for CS inversion involves a requirement for knowledge of psi, which is sensitive to the exact placement of the antennas and of the surrounding scattering media (if a design like that in Figure 3 is employed). However, a given structure may be “calibrated” to infer psi, by simply performing far-field measurements of e, with a source antenna placed at large R within the aforementioned plane; the array response e is then measured as a function of the source angle . By performing this one-time calibration for a sufficient set of angles  element of [0, 2pi (or desired subset of angles), one may directly measure psi. Once this one-time measurement of psi is performed, “standard” CS inversion algorithms may be used to recover x for an arbitrary source J().
This aspect of CS will clearly become tremendously important in the future. Be it the Random Lens Imager at MIT or the Hyperspectral Imager at Duke, whenever we are designing or spreading the Point Spread Functions or PSF of these systems, we are going to need a simple and not so hard way to implement efficient calibration procedures. In some way this is a simpler problem than the question being asked by Yaron Rachlin and Dror Baron on The Secrecy of Compressive Sensing Measurements ( featured earlier) since here, we have some x's and corresponding y's and are trying to figure out the measurement matrix. This is clearly a simple problem if one has enough samples, but the question is, how can the physics help in determining a low bound for the number of measurements necessary to know all the measurement matrix with a certain accuracy.



In other news three papers showed up on my radar screen: A Simple Method for Sparse Signal Recovery from Noisy Observation Using Kalman Filtering by Avishy Carmi, Pini Gurfil, Dimitri Kanevsky. The abstract reads:

We present a simple method for recovering sparse signals from a series of noisy observations. Our algorithm is a Kalman filter (KF) that utilize a so-called pseudo-measurement technqiue for optimizing the convex minimization problem following from the theory of compressed sensing (CS). Compared to the recently introduced CS-KF in [1] which involves the implementation of an additional CS optimization algorithm (e.g., the Dantzig selector), our method is remarkebly easy to implement as it is exclusively based on the KF formulation. The results of an extensive numerical study are provided demonstrating the performance and viability of the new method.

In a primary research report we have introduced a simple method for recovering sparse signals for a series of noisy observations using a Kalman filter (KF). The new method utilizes a so-called pseudo-measuremnet technique for optimizing the convex minimization problem following from the theory of compressed sensing (CS). The CS-embedded KF approach has shown promising results when applied for reconstruction of linear Gaussian sparse processes.
This report presents an imporved version fo the KF algorithm for the recovery of sparse signals. In this work we substitute the l_1 norm by an appropriate quasi-norm l_p, 0 p 1. This modification which better suits the original combinatorial problem, greatly improves the accuracy of the resulting KF algorithm. This however involves the implementation of an extended KF (EKF) stage for properly computing the state statistics.
If we are talking EKF, I wonder how better performing an algorithms based on the Unscented Kalman Filter (an a-ah moment for me) would do.

Finally, Post-silicon Timing Characterization by Compressed Sensing by Davood Shamsi, Petros Boufounos, Farinaz Koushanfar. The abstract reads

We address post-silicon timing characterization of the unique gate delays and their distributions on each manufactured IC. Our proposed approach is based upon the new theory of compressed sensing that enables efficient sampling and reconstruction of sparse signals using significantly fewer measurements than previously thought possible. The first step in performing timing measurements is to find the sensitizable paths via traditional testing methods. Next, we show that variations are sparse in wavelet domain. Then, using the compressed sensing theory, our method estimates the delay distributions by a small number of timing measurements. We discuss how the post-silicon characterization method can enable a range of new and emerging applications including improved simulation models, post-silicon optimization, and IC fingerprinting. Experimental results on benchmark circuits show that using compressed sensing theory can characterize the post-silicon variations with a mean accurately of 95% in the pertinent sparse basis.
Mike Wakin has written a course on Connexions on Compressive Sensing. Check it out.

A student conference at MIT will feature Srikanth Jagabathula who will talk about Inferring Rankings Under Constrained Sensing. The abstract of his talk reads:

Motivated by applications like elections, web-page ranking, revenue maximization etc., we consider the question of inferring popular rankings using constrained data. More specifically, we consider the problem of inferring a probability distribution over the group of permutations using its first order marginals. We first prove that it is not possible to recover more than O(n) permutations over n elements with the given information. We then provide a simple and novel algorithm that can recover up to O(n) permutations under a natural stochastic model; in this sense, the algorithm is optimal. In certain applications, the interest is in recovering only the most popular (or mode) ranking. As a second result, we provide an algorithm based on the Fourier Transform over the symmetric group to recover the mode under a natural majority condition; the algorithm turns out to be a maximum weight matching on an appropriately defined weighted bipartite graph. The questions considered are also thematically related to Fourier Transforms over the symmetric group and the currently popular topic of compressed sensing.

There is a IEEE Workshop on Statistical Signal Processing in Cardiff, Wales, UK on Aug.31-Sept. 3, 2009. In the submission section, there is a track for compressive sensing and sparse approximation.

  • Submission of Proposals for Special Sessions February 9, 2009
  • Paper submissions April 14, 2009
  • Notification of acceptance June 15, 2009
  • Camera ready paper submission June 29, 2009
  • Advanced Registration June 29, 2009
Credit: Grego!'s flicker photostream, Janis Krum. Images first appeared on Twitter. Plane Crash on Hudson river. I'll definitely consider flying US Airways from now on. Who would have known that Airbus planes floated ?!

Thursday, January 15, 2009

CS: Q/A, Accelerating SENSE Using CS, Regularized SENSE thru Bregman Iterations, SparseSENSE, Bregmanized Nonlocal Regularization for Deconvolution

In light of monday's post featuring Michael Lustig presentation on the combination of Parallel Imaging and Compressive Sensing, here are the comments made on that post and other preprints/papers in relation to the same topic followed by a paper from UCLA. In the comment section of Monday's entry one can read an anonymous commenter asking:


I wonder how random sampling is better than poisson-disc sampling in case of single coil with regard to CS point of view as we have local high and low density areas in case of random sampling?

Miki responded with:
I am not yet sure if random sampling is better. Random sampling has uniform coherence (or lack of) between pixels. Poisson disc has higher coherence between pixels at a distance larger than the disc radius, and lower coherence with closer pixels.
Then Davide, another semi-anonymous reader asked:
Which method you are using to generate poisson disc sampling pattern as it is quite computationally expensive.
Miki kindly responded with:
I used dart throwing using the code available at: 
The computational complexity is negligible. We generated a very large set of a poisson disc sampling pattern (supports a Nyquist rate equivalent of 1024x1024 matrix). Whenever we need to scan, we crop the sampling pattern and stretch it according to the desired field of view and resolution. This is implemented on an MRI scanner and has no delay what so ever.  Originally, I think it took about 20 min to generate that sampling pattern, but we did it only once.
Thanks Miki!

The other papers/preprints I found are: Accelerating SENSE Using Compressed Sensing (also viewable here) by Dong Liang, Bo Liu, JiunJie Wang, Leslie Ying. The abstract reads:
Both parallel magnetic resonance imaging (pMRI) and compressed sensing (CS) are emerging techniques to accelerate conventional MRI by reducing the number of acquired data. The combination of pMRI and CS for further acceleration is of great interests. In this paper, we propose two methods to combine SENSE, one of the standard methods for pMRI, and SparseMRI, a recently proposed method for CS-MRI with Cartesian trajectories. The first method, named SparseSENSE, directly formulates the reconstruction from multi-channel reduced k-space data as the same nonlinear convex optimization problem as SparseMRI, except that the encoding matrix is the Fourier transform of the channel-specific sensitivity modulation. The second method, named CS-SENSE, first employs SparseMRI to reconstruct a set of aliased reduced-field-of-view images in each channel, and then applies Cartesian SENSE to reconstruct the final image. The results from simulations, phantom and in vivo experiments demonstrate that both SparseSENSE and CS-SENSE can achieve a reduction factor higher than those achieved by SparseMRI and SENSE individually, and CS-SENSE outperforms SparseSENSE in most cases.

also of interest in light of yesterday's entry, Regularized SENSE Reconstruction Using Bregman Iterations by Bo Liu, Kevin King, Michael Steckner, Jun Xie, Jinhua Sheng and Leslie Ying. The abstract reads: 
In parallel imaging, the signal-to-noise ratio of SENSE reconstruction is usually degraded by the illconditioning problem, which becomes especially serious at large acceleration factors. Existing regularization methods have been shown to alleviate the problem. However, they usually suffer from image artifacts at high acceleration factors due to the large data inconsistency resulting from heavy regularization. In this paper, we propose Bregman iteration for SENSE regularization. Unlike the existing regularization methods where the regularization function is fixed, the method adaptively updates the regularization function using the Bregman distance at different iterations, such that the iteration gradually removes the aliasing artifacts and recovers fine structures before the noise finally comes back. With a discrepancy principle as the stopping criterion, our results demonstrate that the reconstructed image using Bregman iteration preserves both sharp edges lost in Tikhonov regularization and fines structures missed in total variation regularization, while reducing more noise and aliasing artifacts.
Finally, I don't why I did not cover it before, but here it is: SparseSENSE: Randomly-Sampled Parallel Imaging using Compressed Sensing by Bo Liu, Florian Sebert, Yi Ming Zou, and Leslie Ying. The introduction reads:
Since the advent of compressed Sensing (CS) (1), much effort has been made to apply this new concept to various applications (2,3).The most desirable property of CS in MRI application is that it allows sampling of k-space well below Nyquist sampling rate, while still being able to reconstruct the image if certain conditions are satisfied. Recent work (4,5) applied CS to reduce scanning time in conventional Fourier imaging and demonstrated impressive results. In this abstract, we investigate the structure of parallel imaging encoding matrix, and apply CS to parallel imaging to achieve an even higher reduction in scanning time than what can be achieved by each individual method alone. Our experiments show the the combined method, named SparseSENSE, can achieve a reduction factor higher than the number of channels.
Let us note that:
The proposed reconstruction algorithm is computationally intensive, with a running time of 45 minutes for the phantom data on a 2.8GHz CPU/512MB RAM PC.

And therefore it makes sense to consider Multicores or GPUs to speed these computations as the UCLA folks did for a Bregman based algorithm.

While we are talking about UCLA, here is a tech report from these guys entitled: Bregmanized Nonlocal Regularization for Deconvolution and Sparse Reconstruction by Xiaoqun Zhang, Martin Burger, Xavier Bresson, and Stanley Osher. The abstract reads:
We propose two algorithms based on Bregman iteration and operator splitting technique for nonlocal TV regularization problems. The convergence of the algorithms is analyzed and applications to deconvolution and sparse reconstruction are presented.
I note in the end of the paper that
As expected, the standard TV regularization is not capable of recovering texture patterns presented in these images. The results based on wavelet are obtained by using a daubqf(8) wavelet with maximum decomposition level and optimal thresholding parameter setting in a wide range. Since there is no noise considered in these two examples, we solve the equality constrained problem by activating the continuation option in the GPSR code. The nonlocal regularization schemes with Bregman iteration (BOS/PBOS) achieve the best reconstruction result. Surprisingly, with only few measurements, the image textures are almost perfectly reconstructed by the nonlocal TV regularization. This is because image structures are expressed implicitly in the nonlocal weight function, and the nonlocal regularization process with Bregman iteration provides an efficient way to recover textures without explicitly construct a basis.
Stanley Osher tells me that they intend on making their implementation available at some point in time.

And by the way, good luck Sina

Credit: NASA, Appollo 11 LEM

Wednesday, January 14, 2009

CS: Bayesian Compressive Sensing, Availability of TSW-CS and BCS-VB codes, Coordinate Gradient Descent Method for l_1-regularized Convex Minimization


Back in October, I mentioned, the following preprint entitled Exploiting structure in wavelet-based bayesian compressed sensing by Lihan He and Lawrence Carin. There is now an implementation of the TSW-CS algorithm available as well as another implementation of the Bayesian Compressive Sensing package at Duke.

  • TSW-CS: The TSW-CS implemented by an MCMC approach. The package includes the core TSW-CS code, one example of a 1-dimensional signal and two examples of 2-dimensional images.
    Download: tswcs.zip (Last Updated: Dec. 09, 2008)
  • VB-BCS: The basic BCS implemented via a variational Bayesian approach. The package includes the core VB-BCS code, one example of a 1-dimensional signal and two examples of 2-dimensional images.

    Download: bcs_vb.zip (Last Updated: Dec. 09, 2008)



In applications such as signal processing and statistics, many problems involve finding sparse solutions to under-determined linear systems of equations. These problems can be formulated as a structured nonsmooth optimization problems, i.e., the problem of minimizing l_1-regularized linear least squares problems. In this paper, we propose a block coordinate gradient descent method (abbreviated as CGD) to solve the more general l_1-regularized convex minimization problems, i.e., the problem of minimizing an l_1-regularized convex smooth function. We establish a Q-linear convergence rate for our method when the coordinate block is chosen by a Gauss-Southwell-type rule to ensure sufficient descent. We propose efficient implementations of the CGD method and report numerical results for solving large-scale l_1-regularized linear least squares problems arising in compressed sensing and image deconvolution as well as large-scale l_1-regularized logistic regression problems for feature selection in data classification. Comparison with several state-of-the-art algorithms specifically designed for solving large-scale l_1-regularized linear least squares or logistic regression problems suggests that an efficiently implemented CGD method may outperform these algorithms despite the fact that the CGD method is not specifically designed just to solve these special classes of problems.
Let us note how faster these CGD codes are compared to FPC, GPSR_BB. The authors also made these codes available:
All the links to these codes can be found in the reconstruction section of the big picture.

Tuesday, January 13, 2009

CS: Two reconstruction codes: RecPF and FTVd (v 3.0)


Yin Zhang is at it again, he just released two reconstruction codes with attendant papers. First a new code:

RecPF (Version 1.0)
A MATLAB code for image reconstruction from partial Fourier data that solves models with total-variation and l_1 regularization and an l_2-norm fidelity to fit the available incomplete Fourier data. Co-developed with Junfeng Yang and Wotao Yin. RecPF solved the following modelwhere
-- u is the signal/image to be reconstructed 
-- TV(u) is the total variation regularization term 
-- Ψ is a sparsifying basis 
-- Fp is a partial Fourier matrix 
-- fp is a vector of partial Fourier coefficients 
and a new version of FTVd (version 3.0)
A MATLAB code for image deblurring and denoising that solves the model with total-variation regularization and l_2-norm fidelity. Co-developed with Junfeng Yang, Yilun Wang and Wotao Yin. To recall:
This is a Matlab package for recovering images, gray scale or color, from blurry and noisy observations based on solving one of the following 2 problems: 

min TV(u) + (p/2) ||h*u -f||2 or
 min TV(u) + p ||h*u -f||1,
where f is an input blurry and noise image, u is the output image, h is a blurring kernel, and p>0 is a regularization parameter. The noise can be either Gaussian or impulsive like salt-and-pepper.

Both codes are now listed in the reconstruction section of the big picture.

On top of making the code available, the authors also have created a course on Connexions entitled A Class of Fast Algorithms for Total Variation Image Restoration.

FYI, Connexions is a place to view and share educational material made of small knowledge chunks called modules that can be organized as courses, books, reports, etc. Anyone may view or contribute. Richard Baraniuk makes a case for it in this TED talk:








Credit: NASA/JPL/Space Science Institute, Shadows and Gores.

Monday, January 12, 2009

CS: Combining Parallel Imaging and Compressed Sensing, A short discussion with Michael Lustig


Michael Lustig just pointed my attention to the work he and his co-authors were presenting in th upcoming workshop on data sampling and image reconstruction. There is 1 page abstract on-line and a related presentation on the subject entitled:





I went ahead and first checked the pdf version and asked him the following question:
In the incoherent sampling section (slide 107 in the pdf presentation), one thing I'd like to understand, why are you saying that it is has "too many holes", do you mean to say you cannot control the accuracy of the signal in certain region of the frequency space because you don't gather information there ?
Miki responded with:
In the figure, only the white pixels are collected.

There are several aspects here. We are trying to use two separate sources of prior information on the signal in order to recover it.
The first source is that we use a multiple receiver array (slide 7,8), where each array element receives the image information weighted by receive sensitivity of the element, so there is redundancy there. The second, is the compressibility of the signal (slide 9,10).

If the undersampling is much less than the number of elements, the first source is enough to reconstruct, as the measurement matrix is over-determined. But, as the undersampling gets close to the number of elements and beyond, the system of equations becomes ill-conditioned (see slide 88, that shows the noise amplification due to the ill-conditionedness) and then under-determined.  -

The rule of thumb is that when points are relatively close in frequency, the reconstruction is well conditioned, but when they are far away, the reconstruction becomes ill-conditioned.

Now, in order to exploit  sparsity, we need to have an incoherent sampling pattern. But it turns out, that just choosing samples at random is not a good idea, because the distances between samples are not preserved. Globally,  such a sampling scheme has uniform density, but locally, you will get high density areas and "holes". These "holes" really mess up the receiver array source of redundancy. On the other hand, using a Poisson disc sampling distribution provides local uniform sampling everywhere, and also incoherency. So, both sources of redundancy are maximally exploited.

Slides 86-102 explain how parallel imaging reconstruction is done in practice. But I guess if you are not familiar with the concept it is harder to understand without the animation. Try the quicktime version, it is cool.


The conditioning of the parallel imaging reconstruction is actually more than a rule of thumb, It can  be calculated analytically. But I wanted to simplify the explanation.

I watched the Quicktime version and then asked:

One more thing, since MRI is ahead with regards to other technologiesimplementing CS, can you poinpoint for me a formula/bound in your work that lays out thw analytical calculation you just mentioned ?

Miki kindly responded with 

Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P., SENSE: Sensitivity Encoding for Fast MRI , Magn Reson Med. 1999 Nov;42(5):952-62.

Pruessmann KP, Weiger M, Börnert P, Boesiger P., Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med. 2001 Oct;46(4):638-51.

Thanks Miki!

The first reference pinpoints the issue at hand in gathering a signal from multiple sensors at the same time 

There are actually two kinds of noise that affect SENSE images, i.e., noise in sample values and noise in sensitivity data. The latter, however, can usually be reduced to a negligible level by smoothing. Then Eq. [8] for the calculation of image noise holds. This equation illustrates two important aspects of noise propagation in SENSE reconstruction.

First, with multiple channels the diagonal entries in \Psi vary from channel to channel and there is noise correlation between samples taken simultaneously, i.e., there are non-zero cross-terms. Second, unlike a matrix representation of FFT, a SENSE reconstruction matrix generally is not unitary. As a consequence, unlike standard Fourier images the noise level in a SENSE image varies from pixel to pixel and there is noise correlation between pixels.

Thursday, January 08, 2009

CS: A Theoretical Analysis of Joint Manifolds and the release of the Ann Arbor Fast Fourier Transform

If you recall,  Richard Baraniuk made a presentation on the subject entitled "Manifold models for signal acquisition, compression, and processing" where the slides are here and the attendant video are here. There is now an attendant preprint entitled: A Theoretical Analysis of Joint Manifolds by Mark Davenport, Chinmay Hegde, Marco Duarte, and Richard Baraniuk. The abstract reads:
The emergence of low-cost sensor architectures for diverse modalities has made it possible to deploy sensor arrays that capture a single event from a large number of vantage points and using multiple modalities. In many scenarios, these sensors acquire very high-dimensional data such as audio signals, images, and video. To cope with such high-dimensional data, we typically rely on low-dimensional models. Manifold models provide a particularly powerful model that captures the structure of high-dimensional data when it is governed by a low-dimensional set of parameters. However, these models do not typically take into account dependencies among multiple sensors. We thus propose a new joint manifold framework for data ensembles that exploits such dependencies. We show that simple algorithms can exploit the joint manifold structure to improve their performance on standard signal processing applications. Additionally, recent results concerning dimensionality reduction for manifolds enable us to formulate a network-scalable data compression scheme that uses random projections of the sensed data. This scheme efficiently fuses the data from all sensors through the addition of such projections, regardless of the data modalities and dimensions.
as the authors note:

This method enables a novel scheme for compressive, multi-modal data fusion; in addition, the number of random projections required by this scheme is only logarithmic in the number of sensors J.

Before the Christmas break, I talked to Mark Iwen about whether he would make his FFT algorithm available. As you all know, FFT is a cornerstone of most of the scientific revolution of the past 40 years. An algorithm that potentially improves on it by using a sparsity prior is ground breaking. Mark said that some work needed to be done to release it and that he would try to make it available at the end of January. It looks like he made his Ann Arbor Fast Fourier Transform available for download at sourceforge at the beginning of this week! This is the code implemented in Empirical Evaluation of a Sub-Linear Time Sparse DFT Algorithm. Since it is on Sourceforge, I am sure that Mark would not mind having contributors to this project. One of the nice thing to have would include a compiled version for several platforms and even a .dll for those of us using matlab....

Thank you Mark!

CS: Construction of a Large Class of Deterministic Sensing Matrices that Satisfy a Statistical Isometry Property

Here is a new addition at the Rice repository by Robert Calderbank, Stephen Howard and Sina Jafarpour entitled Construction of a Large Class of Deterministic Sensing Matrices that Satisfy a Statistical Isometry Property. The abstract reads:
Compressed Sensing aims to capture attributes of a signal using very few measurements. The Restricted Isometry Property is the condition that the sensing matrix acts as as near isometry on all k-sparse signals. Cand`es and Tao showed that this condition is sufficient for sparse reconstruction and that random matrices, where the entries are generated by an iid Gaussian or Bernoulli process, satisfy the RIP with high probability. This approach treats all k-sparse signals equally likely, in contrast to mainstream signal processing where the filtering is deterministic, and the signal is described probabilistically. In the mainstream framework the sensing matrix is deterministic and it is required to act as a near-isometry on k-sparse vectors with high probability. This paper provides weak conditions that are sufficient to show that a deterministic sensing matrix satisfies this Statistical Restricted Isometry Property (STRIP). The proof is elementary and avoids intricate combinatorial arguments involving coherence of orthonormal bases. The new framework encompases many families of deterministic sensing matrices, including those formed from discrete chirps, Delsarte-Goethals codes, and Extended BCH codes. It is resilient to noise, and generalizes to k-compressible signals, where only k entries are significant, and the magnitude of all remaining entries is close to zero.

As some of you have guessed I am interested in the engineering perspective where one is given a measurement matrix by the physics and the challenge is to find out if that measurement system can fit into an underdetermined system satisfying some property where Compressive Sensing apply. After reading this paper, I went ahead and ask dumb questions to Sina Jafarpour, one of the authors:
Question:
....My question entails the verification that a particular matrix satisfies the STRIP sufficient condition as featured in your Th 2.4:

"1) The columns of Φ form a group UC under pointwise multiplication.
2) the rows of Φ are orthogonal, and all row sums are equal to zero."

Before doing any type of reconstruction or acquiring any signal with it,

* the verification of item 2 means O(N^3) or more operations, as it could evaluated through an SVD.

* for item 1, since a group is a set that is closed under a binary associative operation, contains an identity element, and has an inverse for every element. I am guessing it's a also a O(C^2) maybe even O(C^3) operations but I may be mistaken.

Sina kindly responded with:
...That's generally true, if I give you a (possibly new) matrix you need this amount of calculation to figure out if the conditions for Th 2.4 holds. But if the matrix (or better to say family of matrices) comes with a well-definied structure (like chirps, 2nd order reed-mullers etc) there is no need to check, the structure guarantees that the matrix satisfies conditions of Th 2.4, and I believe the main purpose of StRIP is to finish the proof of recovery of those matrices.

Igor here, I realize there is a need for known matrices like chirps, 2nd order reed-mullers to fit a sufficient condition however, I'd like to think it could be used for other, not yet found families of measurement matrices. This new condition is very attractive compared to the NP-hard RIP. I also wonder how the set of families that is StRIP fare with the families satisfying the JN condition. Thanks Sina!

Printfriendly