The “simplex of probabilty measures” and semi-continuous compressed sensing

Let be a compact subset of and consider the space of continuous functions with the usual supremum norm. The Riesz Representation Theorem states that the dual space of is in this case the set of all Radon measures, denoted by and the canonical duality pairing is given by

We can equip with the usual notion of weak* convergence which read as

We call a measure positive if implies that . If a positive measure satisfies (i.e. it integrates the constant function with unit value to one), we call it a probability measure and we denote with the set of all probability measures.

Example 1 Every non-negative integrable function with induces a probability measure via

Quite different probability measures are the -measures: For every there is the -measure at this point, defined by

In some sense, the set of probability measure is the generalization of the standard simplex in to infinite dimensions (in fact uncountably many dimensions): The -measures are the extreme points of and since the set is compact in the weak* topology, the Krein-Milman Theorem states that is the weak*-closure of the set of convex combinations of the -measures – similarly as the standard simplex in is the convex combination of the canonical basis vectors of .

Remark 1 If we drop the positivity assumption and form the set

we have the is the set of convex combinations of the measures (). Hence, resembles the hyper-octahedron (aka cross polytope or -ball).

I’ve taken the above (with almost similar notation) from the book “ A Course in Convexity” by Alexander Barvinok. I was curious to find (in Chapter III, Section 9) something which reads as a nice glimpse on semi-continuous compressed sensing: Proposition 9.4 reads as follows

Proposition 1 Let , and suppose that the subset of consisting of the probability measures such that for

is not empty. Then there exists such that

1. and are convex combinations of at most -measures, and
2. it holds that for all we have

   

In terms of compressed sensing this says: Among all probability measures which comply with the data measured by linear measurements, there are two extremal ones which consists of -measures.

Note that something similar to “support-pursuit” does not work here: The minimization problem does not make much sense, since for all .

 

ISMP over – non-convex and non-smooth minimization, l^1 and l^p

ISMP is over now and I’m already home. I do not have many things to report on from the last day. This is not due the lower quality of the talks but due to the fact that I was a little bit exhausted, as usual at the end of a five-day conference. However, I collect a few things for the record:

• In the morning I visited the semi-planary by Xiaojun Chenon non-convex and non-smooth minimization with smoothing methods. Not surprisingly, she treated the problem

  

  with convex and smooth and . She proposed and analyzed smoothing methods, that is, to smooth the problem a bit to obtain a Lipschitz-continuous objective function , minimizing this and then gradually decreasing . This works, as she showed. If I remember correctly, she also treated “iteratively reweighted least squares” as I described in my previous post. Unfortunately, she did not include the generalized forward-backward methods based on -functions for non-convex functions. Kristian and I pursued this approach in our paper Minimization of non-smooth, non-convex functionals by iterative thresholding and some special features of our analysis include:

  • A condition which excludes some (but not all) local minimizers from being global.
  • An algorithm which avoids this non-global minimizers by carefully adjusting the steplength of the method.
  • A result that the number of local minimizers is still finite, even if the problem is posed in and not in .

  Most of our results hold true, if the -quasi-norm is replaced by functions of the form

  

  with special non-convex , namely fulfilling a list of assumptions like

  • for (infinite slope at ) and for (mild coercivity),
  • strictly convex on and for ,
  • for each there is such that for it holds that , and
  • local integrability of some section of .

  As one easily sees, -quasi-norms fulfill the assumptions and some other interesting functions as well (e.g. some with very steep slope at like ).

• Jorge Nocedalgave a talk on second-order methods for non-smooth problems and his main example was a functional like

  

  with a convex and smooth , but different from Xiaojun Chen, he only considered the -norm. His talked is among the best planary talks I have ever attended and it was a great pleasure to listen to him. He carefully explained things and put them in perspective. In the case he skipped slides, he made me feel that I either did not miss an important thing, or understood them even though he didn’t show them He argued that it is not necessarily more expensive to use second order information in contrast to first order methods. Indeed, the -norm can be used to reduce the number of degrees of freedom for a second order step. What was pretty interesting is, that he advocated semismooth Newton methods for this problem. Roland and I pursued this approach some time ago in our paper A Semismooth Newton Method for Tikhonov Functionals with Sparsity Constraints and, if I remember correctly (my notes are not complete at this point), his family of methods included our ssn-method. The method Roland and I proposed worked amazingly well in the cases in which it converged but the method suffered from non-global convergence. We had some preliminary ideas for globalization, which we could not tune enough to retain the speed of the method, and abandoned the topic. Now, that the topic will most probably be revived by the community, I am looking forward to fresh ideas here.

 

ISMP – alternatingly projecting on non-convex sets and demixing by convex optimization

Today I report on two things I came across here at ISMP:

• The first is a talk by Russell Luke on Constraint qualifications for nonconvex feasibility problems. Luke treated the NP-hard problem of sparsest solutions of linear systems. In fact he did not tackle this problem but the problem to find an -sparse solution of an system of equations. He formulated this as a feasibility-problem (well, Heinz Bauschke was a collaborator) as follows: With the usual malpractice let us denote by the number of non-zero entries of . Then the problem of finding an -sparse solution to is:

  

  In other words: find a feasible point, i.e. a point which lies in the intersection of the two sets. Well, most often feasibility problems involve convex sets but here, the first one given by this “-norm” is definitely not convex. One of the simplest algorithms for the convex feasibility problem is to alternatingly project onto both sets. This algorithm dates back to von Neumann and has been analyzed in great detail. To make this method work for non-convex sets one only needs to know how to project onto both sets. For the case of the equality constraint one can use numerical linear algebra to obtain the projection. The non-convex constraint on the number of non-zero entries is in fact even easier: For the projection onto consists of just keeping the largest entries of while setting the others to zero (known as the “best -term approximation”). However, the theory breaks down in the case of non-convex sets. Russell treated problem in several papers (have a look at his publication page) and in the talk he focused on the problem of constraint qualification, i.e. what kind of regularity has to be imposed on the intersection of the two sets. He could shows that (local) linear convergence of the algorithm (which is observed numerically) can indeed be justified theoretically. One point which is still open is the phenomenon that the method seems to be convergent regardless of the initialization and that (even more surprisingly) that the limit point seems to be independent of the starting point (and also seems to be robust with respect to overestimating the sparsity ). I wondered if his results are robust with respect to inexact projections. For larger problems the projection onto the equality constraint are computationally expensive. For example it would be interesting to see what happens if one approximates the projection with a truncated CG-iteration as Andreas, Marc and I did in our paper on subgradient methods for Basis Pursuit.

• Joel Tropp reported on his paper Sharp recovery bounds for convex deconvolution, with applications together with Michael McCoy. However, in his title he used demixing instead of deconvolution (which, I think, is more appropriate and leads to less confusion). With “demixing” they mean the following: Suppose you have two signals and of which you observe only the superposition of and a unitarily transformed , i.e. for a unitary matrix you observe

  

  Of course, without further assumptions there is no way to recover and from the knowledge of and . As one motivation he used the assumption that both and are sparse. After the big bang of compressed sensing nobody wonders that one turns to convex optimization with -norms in the following manner:

  

  This looks a lot like sparse approximation: Eliminating one obtains the unconstraint problem \begin{equation*} \min_y \|z_0-Uy\|_1 + \lambda \|y\|_1. \end{equation*}

  Phrased differently, this problem aims at finding an approximate sparse solution of such that the residual (could also say “noise”) is also sparse. This differs from the common Basis Pursuit Denoising (BPDN) by the structure function for the residual (which is the squared -norm). This is due to the fact that in BPDN one usually assumes Gaussian noise which naturally lead to the squared -norm. Well, one man’s noise is the other man’s signal, as we see here. Tropp and McCoy obtained very sharp thresholds on the sparsity of and which allow for exact recovery of both of them by solving (1). One thing which makes their analysis simpler is the following reformulation: The treated the related problem \begin{equation*} \min_{x,y} \|x\|_1 \text{such that} \|y\|_1\leq\alpha, x+Uy=z_0 \end{equation*} (which I would call this the Ivanov version of the Tikhonov-problem (1)). This allows for precise exploitation of prior knowledge by assuming that the number is known.

  First I wondered if this reformulation was responsible for their unusual sharp results (sharper the results for exact recovery by BDPN), but I think it’s not. I think this is due to the fact that they have this strong assumption on the “residual”, namely that it is sparse. This can be formulated with the help of -norm (which is “non-smooth”) in contrast to the smooth -norm which is what one gets as prior for Gaussian noise. Moreover, McCoy and Tropp generalized their result to the case in which the structure of and is formulated by two functionals and , respectively. Assuming a kind of non-smoothness of and the obtain the same kind of results and especially matrix decomposition problems are covered.

 

ISMP, second day

The second day of ISMP started (for me) with the session I organized and chaired.

The first talk was by Michael Goldman on Continuous Primal-Dual Methods in Image Processing. He considered the continuous Arrow-Hurwitz method for saddle point problems

with convex in the first and concave in the second variable. The continuous Arrow-Hurwitz method consists of solving

His talk evolved around the problem if comes from a functional which contains the total variation, namely he considered

with the additional constraints and . For the case of he presented a nice analysis of the problem including convergence of the method to a solution of the primal problem and some a-posteriori estimates. This reminded me of Showalters method for the regularization of ill-posed problems. The Arrow-Hurwitz method looks like a regularized version of Showalters method and hence, early stopping does not seem to be necessary for regularization. The related paper is Continuous Primal-Dual Methods for Image Processing.

The second talk was given by Elias Helou and was on Incremental Algorithms for Convex Non-Smooth Optimization with Applications in Image Reconstructions. He presented his work on a very general framework for problems of the class

with a convex function and a convex set . Basically, he abstracted the properties of the projected subgradient method. This consists of taking subgradient descent steps for followed by projection onto iteratively: With a subgradient this reads as

he extracted the conditions one needs from the subgradient descent step and from the projection step and formulated an algorithm which consists of successive application of an “optimality operator” (replacing the subgradient step) and a feasibility operator (replacing the projection step). The algorithm then reads as

and he showed convergence under the extracted conditions. The related paper is , Incremental Subgradients for Constrained Convex Optimization: a Unified Framework and New Methods.

The third talk was by Jerome Fehrenbach on Stripes removel in images, apllications in microscopy. He considered the problem of very specific noise which is appear in the form of stripes (and appears, for example, “single plane illumination microscopy”). In fact he considered a little more general case and the model he proposed was as follows: The observed image is

i.e. the usual sum of the true image and noise . However, for the noise he assumed that it is given by

i.e. it is a sum of different convolutions. The are kind of shape-functions which describe the “pattern of the noise” and the are samples of noise processes, following specific distributions (could be white noise realizations, impulsive noise or something else)-. He then formulated a variational method to identify the variables which reads as

Basically, this is the usual variational approach to image denoising, but nor the optimization variable is the noise rather than the image. This is due to the fact that the noise has a specific complicated structure and the usual formulation with is not feasible. He used the primal-dual algorithm by Chambolle and Pock for this problem and showed that the method works well on real world problems.

Another theme which caught my attention here is “optimization with variational inequalities as constraints”. At first glance that sounds pretty awkward. Variational inequalities can be quite complicated things and why on earth would somebody considers these things as side conditions in optimization problems? In fact there are good reasons to do so. One reason is, if you have to deal with bi-level optimization problems. Consider an optimization problem

 

with convex and (omitting regularity conditions which could be necessary to impose) depending on a parameter . Now consider the case that you want to choose the parameter in an optimal way, i.e. it solves another optimization problem. This could look like

 

Now you have an optimization problem as a constraint. Now we use the optimality condition for the problem~(1): For differentiable , solves~(1) if and only if

In other words: We con reformulate (2) as

 

And there it is, our optimization problem with a variational inequality as constraint. Here at ISMP there are entire sessions devoted to this, see here and here.

 

ISMP first day

The scientific program at ISMP started today and I planned to write a small personal summary of each day. However, it is a very intense meeting. Lot’s of excellent talks, lot’s of people to meet and little spare time. So I’m afraid that I have to deviate from my plan a little bit. Instead of a summary of every day I just pick out a few events. I remark that these picks do not reflect quality, significance or something like this in any way. I just pick things for which I have something to record for personal reasons.

My day started after the first plenary which the session Testing environments for machine learning and compressed sensing in which my own talk was located. The session started with the talk by Michael Friedlander of the SPOT toolbox. Haven’t heard of SPOT yet? Take a look! In a nutshell its a toolbox which turns MATLAB into “OPLAB”, i.e. it allows to treat abstract linear operators like matrices. By the way, the code is on github.

The second talk was by Katya Scheinberg (who is giving a semi-planary talk on derivative free optimization at the moment…). She talked about speeding up FISTA by cleverly adjusting step-sizes and over-relaxation parameters and generalizing these ideas to other methods like alternating direction methods. Notably, she used the “SPEAR test instances” from our project homepage! (And credited them as “surprisingly hard sparsity problems”.)

My own talk was the third and last one in that session. I talked about the issue of constructing test instance for Basis Pursuit Denoising. I argued that the naive approach (which takes a matrix , a right hand side and a parameter and let some great solver run for a while to obtain a solution ) may suffer from “trusted method bias”. I proposed to use “reverse instance construction” which is: First choose , and the solution and the construct the right hand side (I blogged on this before here).

Last but not least, I’d like to mention the talk by Thomas Pock: He talked about parameter selection on variational models (think of the regularization parameter in Tikhonov, for example). In a paper with Karl Kunisch titled A bilevel optimization approach for parameter learning in variational models they formulated this as a bi-level optimization problem. An approach which seemed to have been overdue! Although they treat somehow simple inverse problems (well, denoising) (but with not so easy regularizers) it is a promising first step in this direction.

 

Estimating unknown sparsity, sparse signal processing and proximal Newton-type methods

In this post I just collect a few papers that caught my attention in the last moth.

I begin with Estimating Unknown Sparsity in Compressed Sensing by Miles E. Lopes. The abstract reads

Within the framework of compressed sensing, many theoretical guarantees for signal reconstruction require that the number of linear measurements exceed the sparsity of the unknown signal . However, if the sparsity is unknown, the choice of remains problematic. This paper considers the problem of estimating the unknown degree of sparsity of with only a small number of linear measurements. Although we show that estimation of is generally intractable in this framework, we consider an alternative measure of sparsity , which is a sharp lower bound on , and is more amenable to estimation. When is a non-negative vector, we propose a computationally efficient estimator , and use non-asymptotic methods to bound the relative error of in terms of a finite number of measurements. Remarkably, the quality of estimation is dimension-free, which ensures that is well-suited to the high-dimensional regime where . These results also extend naturally to the problem of using linear measurements to estimate the rank of a positive semi-definite matrix, or the sparsity of a non-negative matrix. Finally, we show that if no structural assumption (such as non-negativity) is made on the signal , then the quantity cannot generally be estimated when .

It’s a nice combination of the observation that the quotient is a sharp lower bound for and that it is possible to estimate the one-norm and the two norm of a vector (with additional properties) from carefully chosen measurements. For a non-negative vector you just measure with the constant-one vector which (in a noisy environment) gives you an estimate of . Similarly, measuring with Gaussian random vector you can obtain an estimate of .

Then there is the dissertation of Dustin Mixon on the arxiv: Sparse Signal Processing with Frame Theory which is well worth reading but too long to provide a short overview. Here is the abstract:

Many emerging applications involve sparse signals, and their processing is a subject of active research. We desire a large class of sensing matrices which allow the user to discern important properties of the measured sparse signal. Of particular interest are matrices with the restricted isometry property (RIP). RIP matrices are known to enable efficient and stable reconstruction of sfficiently sparse signals, but the deterministic construction of such matrices has proven very dfficult. In this thesis, we discuss this matrix design problem in the context of a growing field of study known as frame theory. In the first two chapters, we build large families of equiangular tight frames and full spark frames, and we discuss their relationship to RIP matrices as well as their utility in other aspects of sparse signal processing. In Chapter 3, we pave the road to deterministic RIP matrices, evaluating various techniques to demonstrate RIP, and making interesting connections with graph theory and number theory. We conclude in Chapter 4 with a coherence-based alternative to RIP, which provides near-optimal probabilistic guarantees for various aspects of sparse signal processing while at the same time admitting a whole host of deterministic constructions.

By the way, the thesis is dedicated “To all those who never dedicated a dissertation to themselves.”

Further we have Proximal Newton-type Methods for Minimizing Convex Objective Functions in Composite Form by Jason D Lee, Yuekai Sun, Michael A. Saunders. This paper extends the well explored first order methods for problem of the type with Lipschitz-differentiable or simple to second order Newton-type methods. The abstract reads

We consider minimizing convex objective functions in composite form

where is convex and twice-continuously differentiable and is a convex but not necessarily differentiable function whose proximal mapping can be evaluated efficiently. We derive a generalization of Newton-type methods to handle such convex but nonsmooth objective functions. Many problems of relevance in high-dimensional statistics, machine learning, and signal processing can be formulated in composite form. We prove such methods are globally convergent to a minimizer and achieve quadratic rates of convergence in the vicinity of a unique minimizer. We also demonstrate the performance of such methods using problems of relevance in machine learning and high-dimensional statistics.

With this post I say goodbye for a few weeks of holiday.

 

Problems solved: RIP and NSP are NP-hard, Homotopy for l1 has exponential complexity

In this post I gladly announce that three problems that bothered me have been solved: The computational complexity of certifying RIP and NSP and the number of steps the homotopy method needs to obtain a solution of the Basis Pursuit problem.

1. Complexity of RIP and NSP

On this issue we have two papers:

• The Computational Complexity of RIP, NSP, and Related Concepts in Compressed Sensing by Andreas M. Tillmann and Marc E. Pfetsch, arXiv/1205.2081
• Certifying the restricted isometry property is hard by Afonso S. Bandeira, Edgar Dobriban, Dustin G. Mixon and William F. Sawin, arXiv/1204.1580

The first paper has the more general results and hence, we start with the second one: The main result of the second paper is this:

Theorem 1 Let a matrix , a positive integer and some be given. It is hard for NP under randomized polynomial-time reductions to check if satisfies the restricted isometry property.

That does not yet say that it’s NP-hard to check if is an RIP constant for -sparse vectors but it’s close. I think that Dustin Mixon has explained this issue better on his blog than I could do here.

In the first paper (which is, by the way, on outcome of the SPEAR-project in which I am involved…) the main result is indeed the conjectured NP-hardness of calculating RIP constants:

Theorem 2 For a given matrix and a positive integer , it is NP-hard to compute the restricted isometry constant.

Moreover, this is just a corollary to the main theorem of that paper which reads as

Theorem 3 For a given matrix and a positive integer , the problem to decide whether satisfies the restricted isometry property of order for some constant is coNP-complete.

They also provide a slightly strengthened version of Theorem~1:

Theorem 4 Let a matrix , a positive integer and some be given. It is coNP-complete to check if satisfies the restricted isometry property.

Moreover, the paper by Pfetsch and Tillmann also proves something about the null space property (NSP):

Definition 5 A matrix satisfies the null space property of order if there is a constant such that for all elements in the null space of it holds that the sum of the largest absolute values of is smaller that times the 1-norm of . The smallest such constant is called the null space constant of order .

Their main result is as follows:

Theorem 6 F or a given matrix and a positive integer , the problem to decide whether satisfies the null space property order for some constant is coNP-complete. Consequently, it is NP-hard to compute the null space constant of .

2. Complexity of the homotopy method for Basis Pursuit

The second issue is about the basis pursuit problem

which can be approximated by the “denoising variant”

What is pretty interesting about the denoising variant is, that the solution (if it is unique throughout) depends on in a piecewise linear way and converges to the solution of basis pursuit for . This leads to an algorithm for the solution of basis pursuit: Start with (for which the unique solution is ), calculate the direction of the “solution path”, follow it until you reach a “break point”, calculate the next direction and so on until hits zero. This is for example implemented for MATLAB in L1Homotopy (the SPAMS package also seems to have this implemented, however, I haven’t used it yet). In practice, this approach (usually called homotopy method) is pretty fast and moreover, only detects a few break points. However, an obvious upper bound on the number of break point is exponential in the number of entries in . Hence, it seemed that one was faced with a situation similar to the simplex method for linear programming: The algorithms performs great an average but the worst case complexity is bad. That this is really true for linear programming is known since some time by the Klee-Minty example, an example for which the simplex method takes an exponential number of steps. What I asked myself for some time: Is there a Klee-Minty example for the homotopy method?

Now the answer is there: Yes, there is!

The denoising variant of basis pursuit is also known as LASSO regularization in the statistics literature and this explains the title of the paper which comes up with the example:

• Complexity Analysis of the Lasso Regularization Path by Julien Mairal and Bin Yu, arxiv.org/1205.0079

Julien and Bin investigate the number of linear segments in the regularization path and first observe that this is upper bounded by is is the number of entries in (i.e. the number of variables of the problem). Then they try to construct an instance that matches this upper bound. They succeed in a clever way: For a given instance with a path with linear segments they try to construct an instance which has one more variable such that the number of linear segments in increased by a factor. Their result goes like this:

Theorem 7 Let have full rank and let be in the range of . Assume that the homotopy path has linear segments and denote by the regularization parameter which corresponds to the smallest kink in the path. Now choose and such that

and define and by

Then the homotopy path for the basis pursuit problem with matrix and right hand side has linear segments.

With this theorem at hand, it is straightforward to recursively build a “Mairal-Yu”-example which matches the upper bound for the number of linear segments. The idea is to start with a example and let it grow by one row and one column according to Theorem~7. We start with the simplest example, namely and . To move to the next bigger example you can choose the next entry and we always choose for convenience. Moreover, you need the next and you need to know the smallest kink in the path. I calculated the paths and kinks with L1Packv2 by Ignace Loris because it is written in Mathematica and can use exact arithmetics with rational numbers (and you will see, that accuracy will be an issue even for small instances) and seemed bullet proof for me. Let’s see where this idea brings us:

Example 1 (Mairal-Yu example)

• Stage 1: We start with , and . The homotopy path has one kink at (with corresponding solution ) and hence, two linear segments. Now let’s go to the next larger instance:
• Stage 2: We can choose the entry as we like and choose it equals to 1, i.e. our new is

  

  Now we have to choose according to (1), i.e

  

  and we can choose, e.g., which gives our new matrix

  

  The calculation of the new regularization path shows that it has exactly the announced number of 5 segments and the parameter of the smallest kink is .

• Stage 2: Again we choose giving

  

  For the choice of we need that

  

  and we may choose

  

  which gives the next matrix

  

  We calculate the regularization path, observe that it has the predicted 14 segments and that the parameter of the smallest kink is .

• Stage 3: Again we choose giving

  

  For the choice of we need that

  

  and we see that things are getting awkward here…

Proceeding in this way we always increase the number of linear segments for the -case from to in each step and one checks easily that this leads to which is the worst case! If you are interested in the regularization path: I produced picture for the first three dimensions (well, I could not draw a 4d -ball) and here they are:

1d Mairal-Yu example

2d Mairal-Yu example

3d Mairal-Yu example

It is not really easy to perceive the whole paths from the pictures because the magnitude of the entries vary strongly. I’ve drawn the path in red, each kink marked with a small circle. Moreover, I have drawn the according -balls of the respective radii to provide more geometric information.

The paper by Mairal and Yu has more results of the paths if one looks for approximate solutions of the linear system but I will not go into detail about them here.

At least two questions come to mind:

• The Mairal-Yu example is . What is the worst case complexity for the true rectangular case? In other words: What is the complexity for in terms of and ?
• The example and the construction leads to matrices that does not have normed columns and moreover, they are far from being equal in norm. But matrices with normed columns seem to be more “well behaved”. Does the worst case complexity lowers if the consider matrices with unit-norm columns? Probably one can construct a unit-norm example by proper choice of …

 

The elastic-net as augmentation and super-resolution by semi-continuous compressed sensing

Today I would like to comment on two arxiv-preprints I stumbled upon:

1. “Augmented L1 and Nuclear-Norm Models with a Globally Linearly Convergent Algorithm” – The Elastic Net rediscovered

The paper “Augmented L1 and Nuclear-Norm Models with a Globally Linearly Convergent Algorithm” by Ming-Jun Lai and Wotao Yin is another contribution to a field which is (or was?) probably the fastest growing field in applied mathematics: Algorithms for convex problems with non-smooth -like terms. The “mother problem” here is as follows: Consider a matrix , try to find a solution of

or, for

which appeared here on this blog previously. Although this is a convex problem and even has a reformulation as linear program, some instances of this problem are notoriously hard to solve and gained a lot of attention (because their applicability in sparse recovery and compressed sensing). Very roughly speaking, a part of its hardness comes from the fact that the problem is neither smooth nor strictly convex.

The contribution of Lai and Yin is that they analyze a slight perturbation of the problem which makes its solution much easier: They add another term in the objective; for they consider

or

This perturbation does not make the problem smooth but renders it strongly convex (which usually makes the dual more smooth). It turns out that this perturbation makes life with this problem (and related ones) much easier – recovery guarantees still exists and algorithms behave better.

I think it is important to note that the “augmentation” of the objective with an additional squared -term goes back to Zou and Hastie from the statistics community. There, the motivation was as follows: They observed that the pure objective tends to “overpromote” sparsity in the sense that if there are two columns in which are almost equally good in explaining some component of then only one of them is used. The “augmented problem”, however, tends to use both of them. They coined the method as “elastic net” (for reasons which I never really got).

I also worked on elastic-net problems for problems in the form

in this paper (doi-link). Here it also turns out that the problem gets much easier algorithmically. I found it very convenient to rewrite the elastic-net problem as

which turns the elastic-net problem into just another -penalized problem with a special matrix and right hand side. Quite convenient for analysis and also somehow algorithmically.

2. Towards a Mathematical Theory of Super-Resolution

The second preprint is “Towards a Mathematical Theory of Super-Resolution” by Emmanuel Candes and Carlos Fernandez-Granda.

The idea of super-resolution seems to pretty old and, very roughly speaking, is to extract a higher resolution of a measured quantity (e.g. an image) than the measured data allows. Of course, in this formulation this is impossible. But often one can gain something by additional knowledge of the image. Basically, this also is the idea behind compressed sensing and hence, it does not come as a surprise that the results in compressed sensing are used to try to explain when super-resolution is possible.

The paper by Candes and Fernandez-Granada seems to be pretty close in spirit to Exact Reconstruction using Support Pursuit on which I blogged earlier. They model the sparse signal as a Radon measure, especially as a sum of Diracs. However, different from the support-pursuit-paper they use complex exponentials (in contrast to real polynomials). Their reconstruction method is basically the same as support pursuit: The try to solve

i.e. they minimize over the set of Radon measures under the constraint that certain measurements result in certain given values . Moreover, they make a thorough analysis of what is “reconstructable” by their ansatz and obtain a lower bound on the distance of two Diracs (in other words, a lower bound in the Prokhorov distance). I have to admit that I do not share one of their claims from the abstract: “We show that one can super-resolve these point sources with infinite precision—i.e. recover the exact locations and amplitudes—by solving a simple convex program.” My point is that I can not see to what extend the problem (1) is a simple one. Well, it is convex, but it does not seem to be simple.

I want to add that the idea of “continuous sparse modelling” in the space of signed measures is very appealing to me and appeared first in Inverse problems in spaces of measures by Kristian Bredies and Hanna Pikkarainen.

 

Some Sloan-Fellowships 2012 related to sparsity and signal processing

Recently the recipients of Sloan Fellowships for 2012 has been announced. This is a kind grant/price awarded to young scientists, usually people who are assistant professors on the tenure track, from the U.S. and Canada. While the actual award is not exorbitant (but still large) the Sloan Fellowships seem to be a good indicator for further success in the career and, more importantly, for awaited breakthroughs by the fellows.

This year there are 20 mathematicians among the recipients and it appeared that I knew three of them:

• Rachel Ward (UTexas at Austin), a student of Ingrid Daubechies, has written a very nice PhD Thesis “Freedom through imperfection” on signal processing based on redundancy. Among several interesting works she has a very recent preprint Stable image reconstruction using total variation minimization which I plan to cover in a future post (according to the abstract the paper provides rigorous theory for exact recovery with discrete total variation which is cool).
• Ben Recht (University of Wisconsin, Madison) also works in the field of signal processing (among other fields) and is especially known for his work on exact matrix completion with compressed sensing techniques (together with Emmanuel Candes); I also liked his paper “The Convex Geometry of Linear Inverse Problems” (with Venkat Chandrasekaran, Pablo A. Parrilo, and Alan Willsky) with elaborated on the generalization of -minimization and nuclear norm minimization.
• Greg Blekherman (Georgia Tech) works (among other things) on algebraic geometry and especially on the problem which non-negative polynomials (in more the one variable) can be written as sums of squares of polynomials, a question which is related to one of Hilbert’s 23 problems, namely Hilbert’s seventeenth problem. An interesting thing about non-negative polynomials and sums of squares is that: 1. Checking if a polynomial is non-negative is NP-hard. 2. Checking is a polynomial is a sum of squares can be done fast. 3. It seem that “most” non-negative polynomials are actually sums of squares but it unclear “how many of them”, see Greg’s chapter of the forthcoming book “Semidefinite Optimization and Convex Algebraic Geometry” here.

Congratulations!

 

Sparse recovery of multidimensional signal with Kronecker products

Today I write about sparse recovery of “multidimensional” signal. With “multidimensional” I mean something like this: A one-dimensional signal is a vector while a two-dimensional signal is a matrix . Similarly, -dimensional signal is . Of course, images as two-dimensional signals, come time mind. Moreover, movies are three-dimensional, a hyperspectral 2D image (which has a whole spectrum attached to any pixel) is also three-dimensional), and time-dependent volume-data is four-dimensional.

Multidimensional data is often a challenge due the large amount of data. While the size of the signals is usually not the problem it is more the size of the measurement matrices. In the context of compressed sensing or sparse recovery the signal is measured with a linear operator, i.e. one applies a number of linear functionals to the signal. In the -dimensional case this can be encoded as a matrix and this is where the trouble with the data comes in: If you have megapixel image (which is still quite small) the matrix has a million of columns and if you have a dense matrix, storage becomes an issue.

One approach (which is indeed quite old) to tackle this problem is, to consider special measurement matrices: If the signal has a sparse structure is every slice, i.e. every vectors of the form where we fix all but the -th component, then the Kronecker product of measurement matrices for each dimension is the right thing.

The Kronecker product of two matrices and is the matrix

This has a lot to do with the tensor product and you should read the Wikipedia entry. Moreover, it is numerically advantageous not to build the Kronecker product of dense matrices if you only want to apply it to a given signal. To see this, we introduce the vectorization operator which takes a matrix and stacks its columns into a tall column vector. For matrices and (of fitting sizes) it holds that

So, multiplying from the left and from the right gives the application of the Kronecker product.

The use of Kronecker products in numerical linear algebra is fairly old (for example they are helpful for multidimensional finite difference schemes where you can build Kronecker products of sparse difference operators in respective dimensions). Recently, they have been discovered for compressed sensing and sparse recovery in these two papers: Sparse solutions to underdetermined Kronecker product systems by Sadegh Jokar and Volker Mehrmann and the more recent Kronecker Compressed Sensing by Marco Duarte and Rich Baraniuk.

From these papers you can extract some interestingly simple and nice theorems:

Theorem 1 For matrices with restricted isometry constant of order it holds that the restricted isometry constant of the Kronecker product fulfills

Basically, the RIP constant of a Kronecker product is not better than the worst one but still not too large.

Theorem 2 For matrices with columns normalized to one, it hold that the spark of their Kronecker product fulfills

Theorem 3 For matrices with columns normalized to one, it hold that the mutual coherence of their Kronecker product fulfills