Installing toolboxes and setting up the path.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

In order to recover sparse measures (i.e. sums of Diracs), it makes sense to consider the following regularization (sometimes
called BLASSO for Beurling LASSO in [deCastroGamboa12]) \[ \umin{m \in \Mm(\mathbb{T})} \frac{1}{2}\norm{y-\Phi(m)}^2 + \la \abs{m}(\mathbb{T}) \] where \(\abs{m}(\mathbb{T})\)
is the total variation of the measure \(m\), defined as \[ \abs{m}(\mathbb{T}) = \sup \enscond{ \int_{\mathbb{T}} g(x) d
m(x) }{ g \in \Cc(\mathbb{T}), \: \normi{g} \leq 1 }.\] This formulation of the recovery of sparse Radon measures has recently
received lots of attention in the literature, see for instance the works of [CandesFernandezGranda13,deCastroGamboa12,DuvalPeyre13].

The BLASSO optimization problem is convex but infinite dimensional, and while there exists solvers when \(\Phi\) is measuring
a finite number of Fourier frequency (see [CandesFernandezGranda13], and the related Numerical Tour "Sparse Spikes Deconvolution over the Space of Measures"), they are do not scale well with
the number of frequencies. The vast majority of practitioners thus approximate the recovery over arbitrary measures by a finite
dimensional problem computed over a finite grid \(z = (z_i)_i \in \mathbb{T}^N\), by restricting their attention to measures
of the form \[ m_{a,z} = \sum_{i=1}^N a_i \de_{z_i} \in \Mm(\mathbb{T}). \] For such a discrete measure, one has \( \abs{m_{a,z}}(\mathbb{T})
= \norm{a}_1\), which can be interpreted as the fact that \(\abs{\cdot}(\mathbb{T})\) is the natural extension of the \(\ell^1\)
norm from finite dimensional vectors to infinite dimensional space of measures. Inserting this parameterization inside the
BLASSO leads to the celebrated Basis-Pursuit problem [Chen99], \[ \umin{a \in \RR^N} \frac{1}{2}\norm{y-\Phi_z a}^2 + \la \norm{a}_1. \] where in the following we make use of the notations
\[ \Phi_z a = \sum_{i=1}^N a_i \phi(\cdot-z_i), \] \[ \Phi_z' b = \sum_{i=1}^N b_i \phi'(\cdot-z_i). \] One can understand
the BP problem as performing a nearest neighbor interpolation of the Dirac's locations.

Continuous Basis-Pursuit

To obtain a better approximation of the infinite dimensional problem, [Ekanadham11] proposes to perform a first order approximation of the kernel. This method assumes that the measures are positive.

To ease the exposition, we consider an uniform grid \(z_i=i/N\) of \(N\) points, so that the grid size is \(\De=1/N\). The
C-BP method solves \[ \umin{(a,b) \in \RR^N \times \RR^N} \norm{ y - \Phi_z a - \frac{\De}{2} \Phi'_z b }^2 + \la \norm{a}_1
\] subject to \[ \abs{b} \leq a \qandq a \geq 0 \] where the constraint inequalities should be understood to hold coordinate-wise.
This is a convex optimization problem, which can be solved using traditional conic optimization methods. We detail in this
tour a simple proximal splitting scheme which can be used.

If \((a^\star,b^\star)\) are solutions of C-BP, one recovers an output discrete measure defined by \[ m^\star = \sum_{a^\star_i
\neq 0} a^\star_i \de_{x_i^\star} \qwhereq x_i^\star = z_i + \frac{\De b_i^\star}{2 a_i^\star} . \] The rationale behind C-BP
is to perform a first order Taylor approximation of the kernel \(\Phi\), where the variable \(\tau_i = \De b_i /(2 a_i) \in
[-\De/2,\De/2]\) encodes the horizontal shift of the Dirac location with respect to the grid sample \(z_i\). The landmark
idea introduced in [Ekanadham11] is that, while the optimization is non-convex with respect to \(\tau\), it is convex with respect to \(b\).

We discribed the method using continuous observation \(y \in L^2(\mathbb{T})\). In fact, it can be applied to a discretized
obervation space as well. We set \(P\) to be the number of observation points.

Generate the input signal \(m_{a_0,x_0}\) composed of three spikes. The variable \(\kappa \in [0,1]\) controls the amount
of shift between the locations \(x_0\) and the sampling grid. The larger \(\kappa\), the worse is the approximation performed
by C-BP. Setting \(\kappa\) close to 1 will make C-BP fail to estimate correctly the number of spikes of the input measure.

Generate the true observations \(y=\Phi(m_{a_0,x_0})\). Note that here we assume there is no noise, i.e. \(w_0=0\). It is
important to realize that, even if there is no noise, when the spikes are off the grid, the linear approximation made by C-BP
(and of course also by BP) create an error term which should be interpreted as noise. This thus necessitate the use of a well-chosen
value of \(\lambda\), which should scale like \(\De^2\) (magnitude of the approximation error).

It is thus the perturbation of \(\iota_{\Cc}\) with a linear form. This allows one to compute the proximal operator of the
\(J\) functional as \[ \forall u \in \Hh, \quad \text{Prox}_{\ga J}(u) = \text{Proj}_\Cc\pa{ u - \la \xi } \] where \(\text{Proj}_\Cc\)
is the orthogonal projector on the cone \(\Cc\).

Exercice 1: (check the solution) Compute the Lipschitz constant \(L=\norm{\Ga_z}^2\). Hint: you can use power iteration to estimate the largest eigenvalues of \(\Ga_z^*\Ga_z\). We display here the evolution of the
estimate of this eigenvalue with the number of power iterations.

Exercice 2: (check the solution) Implement the forward-backward algorithm. Monitor the decay of the energy \( \log_{10}( E(u^{(\ell)})/ E^\star -1 ) \)
where \(E^\star\) is an estimate of the minimum energy (obtained by running the algorithm with a lot of iterations.

exo2;

Recovered parameter \(a,b,x\) that are intended to approximate \(a_0,b_0,x_0\).

Exercice 3: (check the solution) Compute the full homotopy path \(\la \mapsto (a_\la,b_\la) \) where \( (a_\la,b_\la) \) is the solution of C-BP with regularization
parameter \(\la\). The display bellow shows in color the evolution of the correctly estimated spikes and in black wrong perturbating
spikes. Test with different values of \(\kappa\). What can you conclude about the performances of C-BP ?