Hamiltonian Monte Carlo Acceleration Using Surrogate Functions with Random Bases

 We just had a meetup on Stan where Eric mentioned Michael Betancourt's video presentation on Hamiltonian Monte Carlo (see below) and I note the Random Features expansion can speed some of these computations:

Hamiltonian Monte Carlo Acceleration Using Surrogate Functions with Random Bases by Cheng Zhang, Babak Shahbaba, Hongkai Zhao

For big data analysis, high computational cost for Bayesian methods often limits their applications in practice. In recent years, there have been many attempts to improve computational efficiency of Bayesian inference. Here we propose an efficient and scalable computational technique for a state-of-the-art Markov Chain Monte Carlo (MCMC) methods, namely, Hamiltonian Monte Carlo (HMC). The key idea is to explore and exploit the structure and regularity in parameter space for the underlying probabilistic model to construct an effective approximation of its geometric properties. To this end, we build a surrogate function to approximate the target distribution using properly chosen random bases and an efficient optimization process. The resulting method provides a flexible, scalable, and efficient sampling algorithm, which converges to the correct target distribution. We show that by choosing the basis functions and optimization process differently, our method can be related to other approaches for the construction of surrogate functions such as generalized additive models or Gaussian process models. Experiments based on simulated and real data show that our approach leads to substantially more efficient sampling algorithms compared to existing state-of-the art methods.

The Geometric Foundations of Hamiltonian Monte Carlo
M. J. Betancourt, Simon Byrne, Samuel Livingstone, Mark Girolami

Although Hamiltonian Monte Carlo has proven an empirical success, the lack of a rigorous theoretical understanding of the algorithm has in many ways impeded both principled developments of the method and use of the algorithm in practice. In this paper we develop the formal foundations of the algorithm through the construction of measures on smooth manifolds, and demonstrate how the theory naturally identifies efficient implementations and motivates promising generalizations.

Join the CompressiveSensing subreddit or the Google+ Community or the Facebook page and post there !
Liked this entry ? subscribe to Nuit Blanche's feed, there's more where that came from. You can also subscribe to Nuit Blanche by Email, explore the Big Picture in Compressive Sensing or the Matrix Factorization Jungle and join the conversations on compressive sensing, advanced matrix factorization and calibration issues on Linkedin.

REMBO : Bayesian Optimization in a Billion Dimensions via Random Embeddings - implementation -

I am a little late on this one, but here is a Random Features approach used to enable large scale bayesian computations:

Bayesian Optimization in a Billion Dimensions via Random Embeddings by Ziyu Wang, Masrour Zoghi, Frank Hutter, David Matheson, Nando de Freitas

Bayesian optimization techniques have been successfully applied to robotics, planning, sensor placement, recommendation, advertising, intelligent user interfaces and automatic algorithm configuration. Despite these successes, the approach is restricted to problems of moderate dimension, and several workshops on Bayesian optimization have identified its scaling to high-dimensions as one of the holy grails of the field. In this paper, we introduce a novel random embedding idea to attack this problem. The resulting Random EMbedding Bayesian Optimization (REMBO) algorithm is very simple, has important invariance properties, and applies to domains with both categorical and continuous variables. We present a thorough theoretical analysis of REMBO, including regret bounds that only depend on the problem's intrinsic dimensionality. Empirical results confirm that REMBO can effectively solve problems with billions of dimensions, provided the intrinsic dimensionality is low. They also show that REMBO achieves state-of-the-art performance in optimizing the 47 discrete parameters of a popular mixed integer linear programming solver.

Bayesian Optimization in High Dimensions via Random Embeddings by Ziyu Wang, Masrour Zoghi, Frank Hutter, David Matheson, Nando de Freitas

Bayesian optimization techniques have been successfully applied to robotics, planning, sensor placement, recommendation, advertising, intelligent user interfaces and automatic algorithm configuration. Despite these successes, the approach is restricted to problems of moderate dimension, and several workshops on Bayesian optimization have identified its scaling to high dimensions as one of the holy grails of the field. In this paper, we introduce a novel random embedding idea to attack this problem. The resulting Random EMbedding Bayesian Optimization (REMBO) algorithm is very simple and applies to domains with both categorical and continuous variables. The experiments demonstrate that REMBO can effectively solve high-dimensional problems, including automatic parameter configuration of a popular mixedinteger linear programming solver.

An implementation of REMBO can be found at: https://github.com/ziyuw/rembo

Let us note the potential use of this technique for synthetic gene design

  Join the CompressiveSensing subreddit or the Google+ Community and post there !
Liked this entry ? subscribe to Nuit Blanche's feed, there's more where that came from. You can also subscribe to Nuit Blanche by Email, explore the Big Picture in Compressive Sensing or the Matrix Factorization Jungle and join the conversations on compressive sensing, advanced matrix factorization and calibration issues on Linkedin.

A Unified Algorithmic Framework for Multi-Dimensional Scaling, Philosophy and the practice of Bayesian statistics

Here are some of the reading I took with me in a place away from the interwebs.

As we all know, many issues in machine learning depends crucially on an initial good distance measure, From Suresh's tweet, here is: A Unified Algorithmic Framework for Multi-Dimensional Scaling by Arvind Agarwal, Jeff M. Phillips, Suresh Venkatasubramanian. The abstract reads:

In this paper, we propose a unified algorithmic framework for solving many known variants of \mds. Our algorithm is a simple iterative scheme with guaranteed convergence, and is \emph{modular}; by changing the internals of a single subroutine in the algorithm, we can switch cost functions and target spaces easily. In addition to the formal guarantees of convergence, our algorithms are accurate; in most cases, they converge to better quality solutions than existing methods, in comparable time. We expect that this framework will be useful for a number of \mds variants that have not yet been studied. 

Our framework extends to embedding high-dimensional points lying on a sphere to points on a lower dimensional sphere, preserving geodesic distances. As a compliment to this result, we also extend the Johnson-Lindenstrauss Lemma to this spherical setting, where projecting to a random $O((1/\eps^2) \log n)$-dimensional sphere causes $\eps$-distortion.

A press release can be found here.The Matlab code implementing the Euclidian MDS is here.

Philosophy and the practice of Bayesian statistics by Andrew Gelman, Cosma Rohilla Shalizi. The abstract reads:

A substantial school in the philosophy of science identifies Bayesian inference with inductive inference and even rationality as such, and seems to be strengthened by the rise and practical success of Bayesian statistics. We argue that the most successful forms of Bayesian statistics do not actually support that particular philosophy but rather accord much better with sophisticated forms of hypothetico-deductivism. We examine the actual role played by prior distributions in Bayesian models, and the crucial aspects of model checking and model revision, which fall outside the scope of Bayesian confirmation theory. We draw on the literature on the consistency of Bayesian updating and also on our experience of applied work in social science.

Clarity about these matters should benefit not just philosophy of science, but also statistical practice. At best, the inductivist view has encouraged researchers to fit and compare models without checking them; at worst, theorists have actively discouraged practitioners from performing model checking because it does not fit into their framework.

Cosma wrote a small blog entry on his paper here.

Other items to read include:

* Paul Graham's The top idea in your mind

* From Suresh 's blog here are some interesting links: to read and ponder:

• Measure Concentration by Alexander Barvinok
• The Homepage of Nearest Neighbors and Similarity Search by  Yury Lifshits and its successor page The Similarity Search Webpage maintained by Arnoldo Muller

CS: Seeking CS data where the signal prior is not sparse or noise is non-gaussian

Danny Bickson, a reader of this blog, asks the following interesting question:

Does anyone have/know of measurement data which captures signals with either non-gaussian noise, skewed noise, skewed signal prior, or non-unimodal signal (signal has several peeks not just at zero). Most of the works until now, assume sparse signal, so a Laplacian prior captures sparsity very well.

I initially replied that an answer could include most datasets that depends on power-laws (not strictly sparse) such as natural images but this is obviously not what Danny is after. There are also datasets coming out of hyperspectral work photon counting follows Poisson distribution (see for instance: Performance bounds on compressed sensing with Poisson noise by Rebecca Willett and Maxim Raginsky). This happens not just in hyperspectral work as other types of radiation (i.e. radioactive decay for instance,...) would also fall in that category, but I know of no CS radiation sensing system other than those for photons so far. In a different direction , there are different kinds of noises as well as shown in General Deviants: An Analysis of Perturbations in Compressed Sensing by Matthew A. Herman and Thomas Strohmer, the authors talk about a multiplicative noise and how it relates to the physics studied in different areas:

....It is important to consider this kind of noise since it can account for precision errors when applications call for physically implementing the measurement matrix A in a sensor. In other CS scenarios, such as when A represents a system model, E can absorb errors in assumptions made about the transmission channel. This can be realized in radar [7], remote sensing [8], telecommunications, source separation [5], [6], and countless other problems....

as Matt later explained here:

It is important to consider perturbations to the measurement matrix A since the physical implementation of sensors is never perfect in the real world (thus the matrix E can represent precision errors or other non-ideal phenomena). Viewed in a different way, the matrix A can also model a channel that a signal passes through such as in radar, telecommunications, or in source separation problems. The models of these types of channels always involve assumptions and/or approximations of their physical characteristics. In that way, the matrix E can absorb errors in these assumptions. Factors such as these must be accounted for in real-world devices.

Finally, the only bimodal priors I have ever seen were in sensorimotor studies [2] and they were used to see how humans could learn them and how sleep was essential in that learning process:

It could be that subjects need a consolidation period to adequately learn the distribution. Such improvements in learning contingent upon sleep have also been observed in visual learning [7].

But I am not sure this really fit what Danny was asking for and beyond these examples, I am a little bit at loss. Can any expert help ?




While we are on the subject, Aaron Clauset mentioned that his paper on power laws [1] has finally been accepted in SIAM review. I note from his blog entry the following:

Here's a brief summary of the 24 data sets we looked at, and our conclusions as to how much statistical support there is in the data for them to follow a power-law distribution:

Good:
frequency of words (Zipf's law)

Moderate:
frequency of bird sightings
size of blackouts
book sales
population of US cities
size of religions
severity of inter-state wars
number of citations
papers authored
protein-interaction degree distribution
severity of terrorist attacks

With an exponential cut-off:
size of forest fires
intensity of solar flares
intensity of earthquakes (Gutenberg-Richter law)
popularity of surnames
number of web hits
number of web links, with cut-off
Internet (AS) degree distribution
number of phone calls
size of email address book
number of species per genus

None:
HTTP session sizes
wealth
metabolite degree distribution

One should note that while power-laws are interesting for the purpose of spotting phenomena with potentially sparse sets (larger probability elements dominating the counting process and lower probability elements being compared to noise with the same counting process), those phenomena who do not follow power laws can evidently produce population of sparse sets as well. In the end, we are left with whether a CS counting process would be more appropriate than a more conventional counting process.


[1] Power-law distributions in empirical data by Aaron Clauset, Cosma Rohilla Shalizi, M. E. J. Newman
[2] Körding, KP. and Wolpert, D. (2003) Bayesian Integration with Multimodal priors.

CS: Bayesian Compressive Sensing using Laplace Priors, Distilled Sensing

I do not know if it is a clue of things to come but over the course of two days, two graduate students ( Derin Babacan and Jarvis Haupt ) have set up pages for their projects and wanted to make an announcement on Nuit Blanche. I am gladly obliging with the added word that I am very impressed to see students taking these initiatives. It used to be only Assistant Professors who would go for some self promotion but I am glad students are doing it as it shows their interest in providing the full strength of their ideas to the community.

First, Derin Babacan let me know that his research team has made available two papers and attendant source code on their Bayesian approach to compressive sensing. The two papers are:
Bayesian Compressive Sensing using Laplace Priors by S. Derin Babacan, Rafael Molina, and Aggelos Katsaggelos. The abstract reads:

In this paper we model the components of the compressive sensing (CS) problem, i.e., the signal acquisition process, the unknown signal coefficients and the model parameters for the signal and noise using the Bayesian framework. We utilize a hierarchical form of the Laplace prior to model sparsity of the unknown signal. We describe the relationship among a number of sparsity priors proposed in the literature, and show the advantages of the proposed model including its high degree of sparsity. Moreover, we show that some of the existing models are special cases of the proposed model. We present two algorithms resulting from our model; one global optimization algorithm and one constructive (greedy) algorithm designed for fast reconstruction useful in practical settings. Unlike most existing CS reconstruction methods, both algorithms are fully-automated, i.e., the unknown signal coefficients and all necessary parameters are estimated solely from the observation and therefore no user-intervention is needed. Additionally, the proposed algorithms provide estimates of the uncertainty of the reconstructions. We provide experimental results with synthetic 1D signals and images, and compare with the state-of-the-art CS reconstruction algorithms demonstrating the superior performance of the proposed approach.

and Fast Bayesian Compressive Sensing using Laplace Priors. by S. Derin Babacan, Rafael Molina, and Aggelos Katsaggelos. The abstract reads:

In this paper we model the components of the compressive sensing (CS) problem using the Bayesian framework by utilizing a hierarchical form of the Laplace prior to model sparsity of the unknown signal. This signal prior includes some of the existing models as special cases and achieves a high degree of sparsity. We develop a constructive (greedy) algorithm resulting from this formulation where necessary parameters are estimated solely from the observation and therefore no user-intervention is needed. We provide experimental results with synthetic 1D signals and images, and compare with the state-of-the-art CS reconstruction algorithms demonstrating the superior performance of the proposed approach.

The source code can be found on the project webpage.



The algorithm is compared with BCS, BP, OMP, StOMP with CFAR thresholding (denoted by FAR) and GPSR.


Jarvis Haupt put together a page describing (at a relatively high-level) what Distilled Sensing is and how it works:

http://homepages.cae.wisc.edu/~jhaupt/ds.html

Jarvis also let me know that the latest DS paper is Distilled Sensing: Selective Sampling for Sparse Signal Recovery by Jarvis Haupt, Rui Castro and Robert Nowak. The abstract readss;

A selective sampling procedure called distilled sensing (DS) is proposed, and shown to be an effective method for recovering sparse signals in noise. Based on the notion that it is often easier to rule out locations that do not contain signal than it is to directly identify non-zero signal components, DS is a sequential method that systematically focuses sensing resources towards the signal subspace. This adaptivity in sensing results in rather surprising gains in sparse signal recovery dramatically weaker sparse signals can be recovered using DS compared with conventional non-adaptive sensing procedures.

Jarvis also told me that some toy model illustrating the DS algorithm will be available soon.

Thank you Jarvis and Derin.

CS: Bayesian Compressive Sensing via Belief Propagation

Dror Baron just let me know of the release of a new paper by him, Shriram Sarvotham and Rich Baraniuk entitled Bayesian Compressive Sensing via Belief Propagation whose abstract reads:

This project presents an O(N log^2(N)) decoder using two contributions. First, we use CS-LDPC encoding matrices, which consist mostly of zeros, where nonzero entries are {-1,+1}. Second, we incorporate a prior on the signal model (and possible measurement noise model) via belief propagation (BP), a well-known technique for performing approximate Bayesian inference. Our specific example mostly use a two-state mixture Gaussian signal model, but other models can be incorporated in the code.

Dror also kindly provided me also some context to the paper and provided the attendant Matlab code implementing the CS-BP algorithm as well as a much detailed listing of how the functions should be used. Here is the context:

Linear programming CS decoding received much initial attention, yet requires at least quadratic (or even cubic) computation, which is excessive for decoding large signals such as million-pixel images. Decoding is computationally challenging, because we need to (i) take products of the measurement matrix with the signal vector, and (ii) many algorithms iterate over the data numerous times.

To reduce computation, we propose to use sparse CS encoding matrices that offer fast matrix-vector multiplication. In a previous paper, we used these matrices for decoding noiseless measurements of strictly sparse data. The current paper considers problems where the data is not strictly sparse and measurements may be noisy. To approach these problems, we use a Bayesian setting. We focus on a two-state mixture Gaussian signal model, but more sophisticated signal and noise models can be incorporated. Our decoder relies on belief propagation (BP), a technique that is commonly used for signal estimation in various signal processing applications and for channel decoding of turbo and low density parity check (LDPC) codes. The CS-BP decoder represents the sparse CS encoding matrix by a sparse bipartite graph. In each iteration, nodes that correspond to signal coefficients and measurements estimate their corresponding conditional probabilities and then convey this information to neighboring nodes. Because the encoding matrix is sparse, the bipartite graph is also sparse, and so each CS-BP iteration runs in O(N log(N)) time, where the number of iterations is logarithmic in N. Numerical results are encouraging, and the Matlab code is available online.

Thanks Dror!.

Credit Photo: NASA/JPL-Caltech/Cornell, the road not taken by Opportunity.

CS: Solving Helmholtz, Linear Codes, Nonlinear Signal Modeling, Jointly Sparse Signals, Sparse PCA and RIP, Geometrical Study of MP, Cameras tradeoffs

Tim Lin, an undergraduate student at UBC contacted my last week to let me know about a paper I had not covered and how he and his advisor went about discovering what they did, this is fascinating:

... That said, I'm writing you because I've recently published a paper on Geophysics which you might be interested in mentioning on Nuit Blanche. It involves the concept of using the eigenbasis of a matrix operator as the measurement basis for compressive sensing.

I was looking for ways to improve the efficiency of explicit methods of wavefield propagation (via a matrix operator), working as a research assistant in seismic processing for my undergraduate degree in Physics at UBC. My professor Felix Herrmann and I started exploring the one-way propagator which is diagonalized in the eigenbasis of the Helmholtz operator, and it turns out that under certain boundary conditions the transforms into this eigenbasis greatly resemble the discrete (co)sine transform. Out of curiosity we "compressively sampled" the wavefield under a randomly chosen subset of this basis, applied the diagonalized operation, and then attempted to recover the signal with l1-regularization.

To our surprise, we found that we could successfully recover a propagated wavefield. This is upto the compressive-sensing defined limit. We're able to get away using only 30% of the eigenvectors of the original operator in the worst case. We believe with further research this can be a fresh alternative to "largest eigenvalue" or "largest SVD" based approach to compressing matrix operators. This is of course provided that the wavefield is suitably represented in a sparsifying basis such as curvelets.

This is what I understand, the eigenfunctions of the Helmholtz kernel are widespread (exponential) and the trick is to use the incoherence of this type of basis with the general signal seen in geophysics (concentrated wavelets in 1-d/2-d/3-d) to perform some compressed sensing. This is outstanding as it potentially opens the door to other ways of solving PDEs and other integral equations. Tim Lin and Felix Herrmann describe this in detail in Compressed wavefield extrapolation, the abstract reads:

An explicit algorithm for the extrapolation of one-way wavefields is proposed which combines recent developments in information theory and theoretical signal processing with the physics of wave propagation. Because of excessive memory requirements, explicit formulations for wave propagation have proven to be a challenge in 3-D. By using ideas from “compressed sensing”, we are able to formulate the (inverse) wavefield extrapolation problem on small subsets of the data volume ,thereby reducing the size of the operators. Compressed sensing entails a new paradigm for signal recovery that provides conditions under which signals can be recovered from incomplete samplings by nonlinear recovery methods that promote sparsity of the to-be-recovered signal. According to this theory, signals can successfully be recovered when the measurement basis is incoherent with the representation in which the wavefield is sparse. In this new approach, the eigenfunctions of the Helmholtz operator are recognized as a basis that is incoherent with curvelets that are known to compress seismic wavefields. By casting the wavefield extrapolation problem in this framework, wavefields can successfully be extrapolated in the modal domain, despite evanescent wave modes. The degree to which the wavefield can be recovered depends on the number of missing (evanescent) wave modes and on the complexity of the wavefield. A proof of principle for the “compressed sensing” method is given for wavefield extrapolation in 2-D, together with a pathway to 3-D during which the multiscale and multiangular properties of curvelets in relation to the Helmholtz operator are exploited. The results show that our method is stable, has reduced dip limitations and handles evanescent waves in inverse extrapolation.

Also solving the Helmholtz equation, Lawrence Carin, Dehong Liu and B. Guo just released In Situ Compressive Sensing for Multi-Static Scattering : Imaging and the Restricted Isometry Property,

Compressive sensing (CS) is a framework in which one attempts to measure a signal in a compressive mode, implying that fewer total measurements are required vis-`a-vis direct sampling methods. Compressive sensing exploits the fact that the signal of interest is compressible in some basis, and the CS measurements correspond to projections performed on the basis-function coefficients. In this paper we demonstrate that when a target is situated in the presence of a complicated background medium, the frequency-dependent multi-static fields scattered from the target may be measured in a CS manner; these scattered fields are projected onto the two-way Green’s function of the environment. With knowledge of the Green’s function, and using a mutual coherence principle, one may invert for the multi-static scattered fields of the target based on a relatively small number of CS measurements. This is termed in situ CS, because the medium in which the target is employed is used to constitute the CS projections, via the medium Green’s function. We develop a hierarchical statistical algorithm to jointly invert multiple in situ CS measurements performed across a range of frequencies. This framework is demonstrated based on simulated scattering data. In addition, we examine the “restricted isometry property” associated with proper CS projections, and discuss how the observed CS data may be used as a sufficient statistic for a detector or classifier, without the need to know the background Green’s function.

In a different area, Henry Pfister, Fan Zhang released Compressed Sensing and Linear Codes over Real Numbers which makes a connection between compressed sensing and syndrome source coding. The abstract reads:

Compressed sensing (CS) is a relatively new area of signal processing and statistics that focuses on signal reconstruction from a small number of linear (e.g., dot product) measurements. In this paper, we analyze CS using tools from coding theory because CS can also be viewed as syndrome-based source coding of sparse vectors using linear codes over real numbers. While coding theory does not typically deal with codes over real numbers, there is actually a very close relationship between CS and error-correcting codes over large discrete alphabets. This connection leads naturally to new reconstruction methods and analysis. In some cases, the resulting methods provably require many fewer measurements than previous approaches.

The presentation slides are here.

There is also Gradient Pursuit for Non-Linear Sparse Signal Modelling by Thomas Blumensath and Mike Davies. The abstract reads:

In this paper the linear sparse signal model is extended to allow more general, non-linear relationships and more general measures of approximation error. A greedy gradient based strategy is presented to estimate the sparse coefficients. This algorithm can be understood as a generalisation of the recently introduced Gradient Pursuit framework. Using the presented approach with the traditional linear model but with a different cost function is shown to outperform OMP in terms of recovery of the original sparse coefficients. A second set of experiments then shows that for the nonlinear model studied and for highly sparse signals, recovery is still possible in at least a percentage of cases.

This solver is part of version 0.3 of the ‘Sparsify’ Matlab toolbox. Let us note the following use of this technique for coherent image reconstruction:

This experiment was motivated by coherent imaging applications. Assume complex valued data is measured and the images of interest are the magnitude and phase of the complex data. Further assume that the magnitude and phase is sparse in some transform domain. In compressed sensing for imaging applications, instead of observing the complex valued image data, a subset of coefficients is observed in another domain, such as, for example, the Fourier domain.

Marco F. Duarte, Shriram Sarvotham, Dror Baron, Michael Wakin and Richard Baraniuk just released Performace Limits for Jointly Sparse Signals via Graphical Models (and the attendant Poster). The abstract reads:

The compressed sensing (CS) framework has been proposed for efficient acquisition of sparse and compressible signals through incoherent measurements. In our recent work, we introduced a new concept of joint sparsity of a signal ensemble. For several specific joint sparsity models, we demonstrated distributed CS schemes. This paper considers joint sparsity via graphical models that link the sparse underlying coefficient vector, signal entries, and measurements. Our converse and achievable bounds establish that the number of measurements required in the noiseless measurement setting is closely related to the dimensionality of the sparse coefficient vector. Single signal and joint (single-encoder) CS are special cases of joint sparsity, and their performance limits fit into our graphical model framework for distributed (multi-encoder) CS.

In A Direct Formulation for Sparse PCA Using Semidefinite Programming by Alexandre d'Aspremont, Laurent El Ghaoui, Michael I. Jordan, Gert Lanckriet ( DSPCA source code). I noted the following slides: In the first slide, one can find an unexpected way to compute the restricted isometry constant.

and then one discovers how the Sparse PCA problem is transformed into a semidefinite problem:


When looking at this formulation I could not help myself trying to draw conclusion between this transformation and a previous entry entitled:Trace: What is it good for ? - How Compressed Sensing relates to Dimensionality Reduction

In the big picture, the manifold signal processing is not as popular a subject (red arrows) as the one that touches on the reconstruction. Laurent Jacques and Christophe De Vleeschouwer
contribute to our enlightment by looking at the relationship between dictionaries and (smooth) manifolds in A Geometrical Study of Matching Pursuit Parametrization. The abstract reads:

This paper studies the effect of discretizing the parametrization of a dictionary used for Matching Pursuit decompositions of signals. Our approach relies on viewing the continuously parametrized dictionary as an embedded manifold in the signal space on which the tools of differential (Riemannian) geometry can be applied. The main contribution of this paper is twofold. First, we prove that if a discrete dictionary reaches a minimal density criterion, then the corresponding discrete MP (dMP) is equivalent in terms of convergence to a weakened hypothetical continuous MP. Interestingly, the corresponding weakness factor depends on a density measure of the discrete dictionary. Second, we show that the insertion of a simple geometric gradient ascent optimization on the atom dMP selection maintains the previous comparison but with a weakness factor at least two times closer to unity than without optimization. Finally, we present numerical experiments confirming our theoretical predictions for decomposition of signals and images on regular discretizations of dictionary parametrizations.

Alexandre d'Aspremont also released Subsampling Algorithms for Semidefinite Programming. The abstract reads:

We derive two stochastic gradient algorithms for semidefinite optimization using randomization techniques. One is based on the robust stochastic approximation method and uses random sparsifications of the current iterate to both accelerate eigenvalue computations and reduce memory requirements. The other relies on gradient sampling to reduce the per iteration cost of smooth semidefinite optimization algorithms.

at the conclusion, here is an item of interest:

Overall, progress on these issues is likely to come from a better understanding of the measure concentration phenomenon on eigenvectors. At this point, a lot is known about concentration of eigenvalues of random matrices with independent coefficients but random matrix theory is somewhat silent on eigenvectors.

My little understanding of this is: the sensitivity of eigenvalues to small changes in the value of matrix elements is directly related to how nearly colinear eigenvectors are from each other. Therefore, aren't bounds on eigenvalues of random matrices providing some type of statement of the pseudospectra of these matrices ? and if so, what types of results does it imply on the pseudospectra and by the same token on eigenvectors ? Inquiring mind wants to know.

I also noted the following paper as it tries to give some unifying framework to different types of cameras. It is stands it is a useful paper when it comes to understanding how compressed sensing can provide a new way of reconstructing 4D lightfields. Understanding camera trade-offs through a Bayesian analysis of light field projections by Anat Levin, William T. Freeman, and Fredo Durand. The abstract reads:

Computer vision has traditionally focused on extracting structure, such as depth, from images acquired using thin-lens or pinhole optics. The development of computational imaging is broadening this scope; a variety of unconventional cameras do not directly capture a traditional image anymore, but instead require the joint reconstruction of structure and image information. For example, recent coded aperture designs have been optimized to facilitate the joint reconstruction of depth and intensity. The breadth of imaging designs requires new tools to understand the tradeoffs implied by different strategies. This paper introduces a unified framework for analyzing computational imaging approaches. Each sensor element is modeled as an inner product over the 4D light field. The imaging task is then posed as Bayesian inference: given the observed noisy light field projections and a new prior on light field signals, estimate the original light field. Under common imaging conditions, we compare the performance of various camera designs using 2D light field simulations. This framework allows us to better understand the tradeoffs of each camera type and analyze their limitations.

Compressed Sensing: There is no theory behind it right now, but does it work ? and Lena in 3D

This is the question Raghunandan Kainkaryam and Radu Berinde are debating in the comments section of Raghunandan's blog entry entitled: Shifted Transversal Designs and Expander Graphs. Is there anyone who wants to try it out on their signals ? (i.e. reconstructing signals using this measurement matrix). I put the measurement matrix function in this m.file.

Unrelated to the previous statement, here is Lena in 3D using Make3d that was tried before on out of this world photographs before. You can see the full view here as well as some tornado views (here and here).


What does this say about this picture, but more importantly what does it say about the technique developed by Ashutosh Saxena, Min Sun and Andrew Ng.

Predicting NROL-21 's fall: When and How.

Here is a different kind of debris. Aviation Week reports that an NRO spacecraft is uncontrollably descending into the atmosphere "at a rate of about 2,310 feet per day, according to Canadian analyst Ted Molczan." The reporter seems to be saying that the Columbia fateful reentry radar and recovery data are aiding in the analysis of the potential for accidents with the debris of that satellite. I am a little bit skeptical on several counts. When people tried to evaluate when a small disco ball would reenter the atmosphere, few guessed the right time and day it eventually happened. Also, when I last looked into it, most of the modeling that used specific small scale experiments did not seem to have some real predictive capability. It seems that this conclusion was also reached by some of the SCM folks who eventually did a new kind of statistical study on this very problem for the French Space Agency. Finally, the two photos below show before and after photographs of our camera that was on-board the Columbia Space Shuttle. It is in aluminum and even though the fusion temperature of Aluminum is 600 C, well below the temperature felt outside of the craft on reentry, the base plate was nearly intact on impact (our camera was outside). The rest melted away except for the optical component in Invar.

Make3D: Depth Information From A Single Shot. Stressing the Algorithm With Other World Views

In a previous installment, we saw how coded aperture could use specific bayesian priors to infer depth. Ashutosh Saxena went a little further to make a similar technology more accessible using normal images. He just made available to everybody the ability to produce 3-D models from single snapshots. This is an extension of his work with Min Sun and Andrew Ng. The site is at:

http://Make3d.stanford.edu

If you register with the site, you can upload your own images and see the results either in some Quicktime format or in VRML which can be viewed using the Cortona plugin for Windows, and this extension for Linux.

In order to push a little bit the boundaries of the algorithm, I specifically used imagery that is a little "weird" where the background has very sharp discontinuities: Planetary exploration.

The landscape as seen from Huygens when it landed on Titan is here:


Not bad for a single view (the Huygens probe survived only 30 minutes on Titan as expected and did not move) that no human has ever seen.

I am less impressed by the 3-D view rendered from the Mars rover Opportunity taken a week ago.


(the VRML can be seen here)
or Spirit



Some of the very sharp contrast images from Cassini looking a Phoebe did not converge.


But the beautiful view from the (amateur) HALO flight 2 is very interesting.


When using the 3D viewer, it looks like you are flying over clouds.

You can try it here.
Indoor scenes seem to provide also some good estimate, starting with the original shot:

while the 3D scene produces this:

one can view it with the 3D viewer.

I am half surprised by some of these results, if one recalls how the technique works, for outdoor photography, haze is an important part of the process that allows the method to differentiate between different depths. And so I am only half surprised that it does not converge for Phoebe (no atmosphere) but somewhat well on Mars (low atmosphere) and well on Titan (with an atmosphere) and on Earth at 28 km altitude looking at clouds and indoors.

Credit: ESA/NASA/JPL/University of Arizona/Alexei Karpenko

Building a Roaming Dinosaur: Why is Policy Learning Not Used ?

Dubai is going to host a Jurassic Park. It is about time: the current display of Animatronics are very much underwhelming and certainly do not yield the type of magic moment as displayed by Laura Dern's face in Jurassic Park. Yet, all over the world, kids and families line up to pay from $10 to $100 for the privilege of being wowed. The most interesting feature of the Dubai 'Resteless Planet' undertaking will be the claimed ability for the dinosaurs to be roaming. This is no small feat as none of the current animatronics are able to do that. I have a keen interest in this as you probably have noticed from the different entries on muscles, scaling and various autonomous robotics undertaking.

So over the years, I have kept an eye on the current understanding of dinosaur gait. Examples on the web can be found here and here. A sentence seems to be pretty much a good summary:

When the American Museum of Natural History wanted to create a digital walking Tyrannosaurus rex for a new dinosaur exhibit, it turned to dinosaur locomotion experts John Hutchinson and Stephen Gatesy for guidance.

The pair found the process humbling.

With powerful computers and sophisticated modeling software, animators can take a pile of digital bones and move them in any way they want. This part is easy; but choosing the most likely motion from the myriad of possibilities proved difficult......

The researchers think that one way to narrow down the possibilities of dinosaur movement is to use more rigorous physical constraints in computer models. These constraints fall into two broad categories: kinematic (motion-based) and kinetic (force-based). One simple kinematic constraint, for example, is that the ankle and knee cannot bend backwards.

So in short, it is one thing to know the skeleton, it is another one to devise how the movement goes. Yet, the current methods devised to figure gait are relying on not so sophisticated methods. As it turns out, we have a similar problem in robotics and machine learning. Because robots are becoming increasingly complex, there needs to be new methods of collecting data and summarizing them in what are called 'policies'. New methods are able to learn behavior for robots even though they have many degrees of freedom though some type of supervised learning. Some of the techniques include Non-Negative Matrix Factorization (NMF), diffusion processes and some of the techniques we tried in our unsuccesful attempt in DARPA's race in the future.

[ Update: a dinosaur finding showing preserved organic parts shows us that basing our intuition on just bones is not enough. It looks as though dinosaurs may have been much larger. Ona different note, it is one thing to model human behavior (and by extension dinosaur behavior) using Differential equations,but the problem you are trying to solve is 'given a certain behavior, how can it fit the model set forth by the differential equations?'. This is what is called an inverse problem and while a set of differential equations may give you a sentiment that you are modeling everything right, they generally are simplification of the real joint behavior and their interaction with the environment (soil,...). In short, to give a sense of realness, you have to go beyond a description with differential equations alone, for these reasons alone. For this reason, building a real roaming dinosaur need the type of undertaking mentioned above in this entry ]

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.

L1 -- What is it good for ?

In a previous entry, I mentioned the connection between how L1 and maximum entropy considerations. As we now know, in Compressed Sensing, the reconstruction of the signal from incoherent measurements involves a Linear Programming (L1) step featured in the Basis Pursuit approach. We also know that there is a bayesian approach to the reconstruction which obviously involves Laplace priors (because a maxent approach to the problem involving the L1 norm point to Laplace distribution as ideal priors).

In The Price of Privacy and the Limits of LP Decoding by Cynthia Dwork, Frank McSherry and Kunal Talwar show the following:

Our principal result is the discovery of a sharp threshhold ρ∗ ≈ 0.239, so that if ρ < ρ∗ and A is a random m × n encoding matrix of independently chosen standard Gaussians, where m = O(n), then with overwhelming probability over choice of A, for all x ∈ Rn, LP decoding corrects ρm arbitrary errors in the encoding Ax, while decoding can be made to fail if the error rate exceeds ρ∗.

and then in the conclusion

Comparing the results of section 7 with those of section 4, we note that while near-perfect decoding is information theoretically possible for error rates up to a half, LP decoding fails at much lower error rates. It is thus natural to look for other efficient decoding procedures for the case of adversarial error.

We already somehow knew this from previous investigations (David Donoho and Victoria Stodden in Breakdown Point of Model When the Number of Variables Exceeds the Observations ) but we now have a number: 24 percent. Interestingly, Rick Chartrand in Exact reconstruction of sparse signals via nonconvex minimization makes the case that since L1 might not be optimal and that looking at Lp with p less than 1 might lead to better results

. The idea being that asymptotically, we want p to go to zero in order to fulfill the real L0 requirement. From the abstract:

We show that by replacing the L1 norm with the Lp norm with p < 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.

Another paper by Chartrand, Nonconvex compressed sensing and error correction shows similar results. His computations seem to go very fast compared to BP/LP so this is noteworthy (we could always go for L0 directly but since it is combinatorial, the time spent on the computation is just too enormous). In light of the reconstruction error shown in that article, one cannot but recall the statement made by Aleks Jakulin and Andrew Gelman in this comments section on using log(1+d^2) norm.

As an log(1+d^2) ball does not have (thank you Yaroslav) the same shape as the Lp norm ball for p less than 1, how can we reconcile the findings that Aleks seems to find Cauchy/t-distributions are doing a good job in sparse decomposition (better than L-1)?

For other reasons, I'd like to think that Cauchy distributions are indeed grounded in serious theoretical work (besides the observation that they are not sensitive to outliers). We already know that random projections exist when modeling the primary visual cortex. We may eventually figure that some type of Autism is related to an equivalent phase change between L0 and L1 or L_log(1+d^2). There is a nice parallel between the metabolic constraints on the cortex when doing sparse coding and CPU requirements to do L0 or L1 or Lp.

Resources:
[1] Loss function semantics.
[2] Rice Compressed sensing page.
[3] Object recognition in the cortex

Adding Search and Rescue Capabilities (part II): Modeling what we see and do not see

One of the concern that one has during a search and rescue operation (part I is here) is whether or not, the item of interest was seen or detected. I am not entirely sure for instance that SAROPS includes this, so here is the result of some of the discussions I have had with some friends on this. While the discussions were about the Tenacious, one should keep an eye on how it applies to other types of mishap that may lead to a similar undertaking.

In the search for the Tenacious, there were several sensors used at different times:

• Mark One Eyeball from Coast Guards or from some private parties or from onlookers from the coast
  
• sensors used by the Coast Guard in Planes and Boats
• sensors (Radar, visual, IR, multispectral) from satellites or high altitude planes
• webcams looking at the SF port and bay.

Each and every one of these sensors give some information about their field of view but they are limited by their capabilities. The information from the sensor is dependent on its resolution and other elements. While the issue of resolution is well understood, at least spatially, sensor visibility is dependent on:

• cloud cover (high altitude, satellites), haze (low altitude)
• the calmness of the sea
• the orientation of the sensor (was the object of interest in the sensor cone ?)
• the ability of the sensor to discriminate the target of interest from the background (signature of the target)
  
• the size of the target (are we looking for a full boat or debris ?)
  

And so whenever there is a negative sighting over an area, the statement is really about the inability of the detector to detect the target of interest due the elements listed above. And so the probability of the target of interest not being there is not zero (except in very specific circumstances). In effect, when the data fusion occurs when merging information from all these sensors, it is important to be able to quantify what we don't know as much as what we know. It is also important to realize that different maps are really needed for each scenario. A scenario about searching for debris is different from that of searching for a full size boat. What the detectors/sensors see is different in these two scenarios. While one can expect to have a good signal when searching for a full size boat, most sensors are useless when it comes to detecting minutes debris.

In the aerospace business, some of us use software like STK that provides different modules in order to schedule and understand information about specific satellite trajectories and so forth. It may be a good add-on to the current SAROPS capabilities in terms of quantifying the field of view.



But the main issue is really about building the right probability distribution as the search goes on and how one can add any heterogenous information into a coherent view of the search.

Time is also a variable that becomes more and more important as the search goes. In particular it is important to figure out the ability to do data fusion with time stamped data. One can see in this presentation, that while the search grid is regular, one can see some elements drifting out of the field of view as the search is underway. So the issue is really about quantifying data fusion with sensors input as well as maritime currents and provide a probability of escaping the search grid. SAROPS already does some of this, but I am not sure the timing element of the actual search (made by CG planes, boat) is entered in the software as the search go on. It was difficult for us to get back that timing from the search effort (it was rightfully not their priority) and one simply wonders if this is an input to SAROPS when iterating on the first empty searches. If one thinks along the lines of the 8000 containers scenario, this is important as it has been shown that some of these containers have different lifespan at sea level and right under the surface. In this case, the correlation between time stamped sensor outputs become central as a submerged but within a few feet underwater containers may not be viewable from specific sensors (but would remain dangerous to navigation). Also this is not because we did not see anything on the second path at the same location (provided no current) that the object is not here anymore, rather the sensor did not detect it. In the illustration below one can see the different targets found by the Radarsat/John Hopkins team for the Tenacious. Without time stamp it is nearly impossible to make a correlation between hits on the first and the second satellite path.

The bayesian framework seems to have already been adopted by SAROPS and previous versions. It may need some additional capabilities to take into account most the issues mentioned above (sensor network or EPH). In either case, a challenge of some kind, with real data might be a way to advance the current state of the art.