Proceedings of the second"international Traveling Workshop on Interactions between Sparse models and Technology"(iTWIST'14)

The implicit objective of the biennial"international - Traveling Workshop on Interactions between Sparse models and Technology"(iTWIST) is to foster collaboration between international scientific teams by disseminating ideas through both specific oral/poster presentations and free discussions. For its second edition, the iTWIST workshop took place in the medieval and picturesque town of Namur in Belgium, from Wednesday August 27th till Friday August 29th, 2014. The workshop was conveniently located in"The Arsenal"building within walking distance of both hotels and town center. iTWIST'14 has gathered about 70 international participants and has featured 9 invited talks, 10 oral presentations, and 14 posters on the following themes, all related to the theory, application and generalization of the"sparsity paradigm": Sparsity-driven data sensing and processing; Union of low dimensional subspaces; Beyond linear and convex inverse problem; Matrix/manifold/graph sensing/processing; Blind inverse problems and dictionary learning; Sparsity and computational neuroscience; Information theory, geometry and randomness; Complexity/accuracy tradeoffs in numerical methods; Sparsity? What's next?; Sparse machine learning and inference.

iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.

Introduction
Sparse recovery has become a popular research field especially after the emergence of compressed sensing theory [1]. Consequently many natural signals have also been shown to be compressible or essentially low dimensional which enabled sparsity inducing methods to be used for the estimation of these signals in linear inverse problems. Acoustic source localization can also be formulated as a sparse recovery problem [2,3]. In the continuous space and time, the acoustic pressure field in a 2-dimensional space, p( − → s , t), induced by the pressure emitted by the sources, f ( − → s , t), is known to obey the differential wave equation at all non-boundary locations. The parameter c in (1) represents the speed of sound in the medium. For inhomogeneous mediums, this constant can be represented more generally as c( − → s , t). In addition to the wave equation, the acoustic pressure field is also restricted by the boundary conditions. Both the wave equation and the boundary equations can be discretized straightforwardly into a set of linear differential equations represented as (Ω 1 + αΩ 2 )x = z (2) in which Ω 1 and Ω 2 represent the discretized differential operators in space and time respectively, α 1/c 2 is a function of speed of sound, x and z represent the discretized pressure field and the source signals in vector form respectively. The source localization problem can be described as the estimation of z (or equivalently x) from limited number of measurements recorded by a number of microphones, i.e. from where I s is a matrix formed by keeping some of the rows of the identity matrix corresponding to microphone locations. Considering the sparse nature of the sources, source localization can be posed as a sparse recovery problem which can be solved with convex optimization such as in where the * norm is often chosen as the ℓ 1 norm or the ℓ 1,2 norm [3]. In many such approaches, the medium is assumed to be homogeneous and the medium properties (c) are assumed to be known. However, in practice this assumption may not be always valid since the behavior of the medium may vary depending on the physical conditions such as temperature, ambient pressure, or unknown inhomogeneities in the medium. These incorrect assumptions on α (and hence the sparsifying operator Ω 1 +αΩ 2 ) may lead to significantly decreased sparsity or multiplicative perturbations that would highly reduce the recovery performance [4,5].
In this abstract, we investigate the source localization problem formulated in (4) but with an unknown speed of sound c (or equivalently α). In Section 2, we formulate the source localization problem for a spatially varying c, i.e. a time invariant but inhomogeneous medium and present an algorithm to estimate the medium properties along with the sources under spatial smoothness assumptions. Preliminary experimental results with the proposed algorithm are then presented in Section 3, which is followed by the concluding remarks in Section 4.

Source Localization with Unknown Medium Parameters
Let us first assume that the medium is inhomogeneous and unknown, i.e. α = α( − → s , t) for the continuous domain. Consequently the discretized equation in (2) can be modified as where Diag(α) indicates a diagonal matrix with the entries of the vector α along the diagonal. In this setting one can attempt to recover α and x by minimizinĝ x,α = arg min x,α (Ω 1 + Diag(α)Ω 2 )x 1,2 s.t. I s x = y (6) where the ℓ 1,2 norm is defined as ℓ 1 norm along the spatial and ℓ 2 norm along the temporal dimension, assuming all the source locations are fixed along time. However the global minimum of the objective function in (5) often does not correspond to the actual x and α due to the high number of degrees of freedom introduced by the unknown α.  Figure 1: (a) The minimum sufficient number of the microphones, mx.α(n), with n = 1, 2 for successfully localizing different number of sources. The required number of microphones for the ideal case when the field α is known, mx(n), is also plotted for comparison. (b) The original and estimated signals emitted from the sources for a sample recovery.
The intrinsic degrees of freedom in the structure of a real world medium is often much smaller. In this work we will assume the medium does not change in time and exploit the fact that the spatial variation in a natural medium is often limited and assume nth order spatial difference of the field α along each spatial dimension is zero, i.e. D n α = 0 (7) where A n is an orthogonal basis for the null space of D n . It can be observed that the degrees of freedom for α are reduced to n 2 (for 2-dimensional space) from the total size of the discretized points in space, which can be a significant reduction in practice. With the assumption of (8), the optimization in (6) can be modified aŝ x,p = arg min x,p (Ω 1 + Diag(A n p)Ω 2 )x 1,2 s.t. I s x = y Unlike the case with known α, (9) is not a linear optimization when jointly considering the variables x and α, hence the well developed convex optimization methods are not applicable. It is, however, linear in terms of each of the unknowns when the other is considered a constant. Therefore we propose a form of the alternating minimization approach which is often employed in such scenarios.

Experimental Results
In order to demonstrate the joint estimation performance of the proposed approach, a rectangular grid of 10 by 10 spatial dimensions is simulated with k randomly located stationary sources, k varying from 1 to 5. The sources emit random sinusoidals for a duration of 20 time instants. The pressure field is measured by m randomly located stationary microphones for 100 time instants, m varying from 2 to 10. The boundary conditions are modeled to be Neumann boundaries as in [2,3]. The field α is generated randomly to obey (8) with parameter n set as 1 and 2 and the field values to be between 0.1 and 0.6. The field parameters p and the acoustic pressure x are estimated from the microphone measurements y = I s x with an algorithm that minimizes (9) by alternating between unknowns applying the method of Alternating Direction Method of Multipliers (ADMM [6]). The details of the algorithm are not provided here due to space constraints. Each of the randomly generated experiments is repeated 20 times with the same parameters to empirically estimate the minimum sufficient number of microphones for successful source localization for each set of n, m, k. After each simulation, k points with highest energy are chosen as estimates of the source locations. The minimum number of microphones that succeeds in localization of the sources in 80% of the 20 experiments is set as the minimum sufficient number of microphones, m x,α (n). The minimum sufficient number of the microphones for the recovery of the acoustic field when the medium parameter, α, is known, m x , is also computed the same way for comparison. Figure 1a shows the change of m x.α (n) and m x (n) with respect to the number of sources. A sample reconstruction of source signals is shown in Figure 1b. The original and the estimated field α as well as the source and the microphone locations for the same sample reconstruction can also be seen in Figure 2.
The first thing to notice in the presented results in Figure 1a is that estimating the unknown α introduces almost no performance degradation with respect to perfectly knowing the medium parameters. This can be attributed to the very small degrees of freedom in α thanks to the smoothness constraints. Hence estimating α while localizing sources is possible whenever the sources can be localized with known α. In fact it is observed in the experiments that the estimation of α is still accurate for many cases when source localization is not successful. A second observation is that the sufficient number of microphones is not significantly affected when α is spatially varying (n = 2) compared to when it is spatially constant (n = 1). However the degrees of freedom within α increases with n 2 and therefore the sufficient number of microphones is expected to increase more rapidly with increasing n.

Conclusion
In this work, the acoustic source localization problem with additional variables due to unknown medium properties is described and a sparse recovery approach is presented for estimating both the sources and the medium structure. The presented approach utilizes the assumption that the spatial variation in the medium is limited and the acoustic sources are sparse. Preliminary results that demonstrate the performance of the proposed approach are presented. The talk will include further results, the details on the recovery algorithm and discussions on the extensions of the algorithm for more realistic scenarios.

Introduction
Optimization on manifolds, or Riemannian optimization, is a fast growing research topic in the field of nonlinear optimization. Its purpose is to provide efficient numerical algorithms to solve optimization problems of the form where the search space M is a smooth space: a differentiable manifold which can be endowed with a Riemannian structure. In a nutshell, this means M can be linearized locally at each point x as a tangent space T x M and an inner product ·, · x which smoothly depends on x is available on T x M. A number of smooth search spaces arise often in applications: • The oblique manifold is a product of spheres, M = {X ∈ R n×m : diag(X X) = 1 m }. That is, X ∈ M if each column of X has unit 2-norm in R n . Absil and Gallivan (2006) show how independent component analysis can be cast on this manifold as non-orthogonal joint diagonalization. When furthermore it is only the product Y = X X which matters, matrices of the form QX are equivalent for all orthogonal Q. Quotienting out this equivalence relation yields the fixed-rank elliptope M = {Y ∈ R m×m : Y = Y 0, rank(Y ) = n, diag(Y ) = 1 m }. See the example below for application to the max-cut problem. The packing problem on the sphere, where one wishes to place m points on the unit sphere in R n such that the two closest points are as far apart as possible (Dirr et al., 2007), is another example of an optimization problem on the fixed-rank elliptope. Grubišić and Pietersz (2007) optimize over this set to produce low-rank approximations of covariance matrices. • The Grassmann manifold M = {col(X) : X ∈ R n×m * } is the manifold, where R n×m * is the set of full-rank matrices in R n×m and col(X) denotes the subspace spanned by the columns of X. That is, col(X) ∈ M is a subspace of R n of dimension m. Among other things, optimization over the Grassmann manifold proves useful in low-rank matrix completion, where it is observed that if one knows the column space spanned by the sought matrix, then completing the matrix according to a least squares criterion is easy ( • The set of symmetric, positive semidefinite, fixed-rank matrices is also a manifold, M = {X ∈ R n×n : X = X 0, rank(X) = k}. Meyer et al. (2011) exploit this to propose low-rank algorithms for metric learning. This space is tightly related to the space of Euclidean distance matrices X such that X ij is the squared distance between two fixed points x i , x j ∈ R k . Mishra et al. (2011) leverage this geometry to formulate efficient low-rank algorithms for Euclidean distance matrix completion.
• The fixed-rank spectrahedron M = {X ∈ R n×n : X = X 0, trace(X) = 1 and rank(X) = k}, without the rank constraint, is a convex set which can be used to solve relaxed (lifted) formulations of the sparse PCA problem. Journée et al. (2010) show how optimizing over the fixed-rank spectrahedron can lead to efficient algorithms for sparse PCA.
The rich geometry of Riemannian manifolds M makes it possible to define gradients and Hessians of cost functions f , as well as systematic procedures (called retractions) to move on the manifold starting at a point x, along a specified tangent direction at x. Those are sufficient ingredients to generalize standard nonlinear optimization methods such as gradient descent, conjugate-gradients, quasi-Newton, trust-regions, etc.
iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. (4) In a recent monograph, Absil et al. (2008) lay down a mature framework to analyze problems of the form (1) when f is a smooth function, with a strong emphasis on building a theory that leads to efficient numerical algorithms. In particular, they describe the necessary ingredients to design first-and second-order algorithms on Riemannian manifolds in general. These algorithms come with convergence guarantees essentially matching those of the Euclidean counterparts they generalize. For example, the Riemannian trust-region method is known to converge globally toward critical points and to converge locally quadratically when the Hessian of f is available. In many respects, this theory subsumes well-known results from an earlier paper by Edelman et al. (1998), which focused on problems of the form (1) with M either the set of orthonormal matrices (the Stiefel manifold) or the set of linear subspaces (the Grassmann manifold).
The maturity of the theory of smooth Riemannian optimization, its widespread applicability and its excellent track record performance-wise prompted us to build the Manopt toolbox: a user-friendly piece of software to help researchers and practitioners experiment with these tools. Code and documentation are available at www.manopt.org.

Architecture and features of Manopt
The toolbox architecture is based on a separation of the manifolds, the solvers and the problem descriptions. For basic use, one only needs to pick a manifold from the library, describe the cost function (and possible derivatives) on this manifold and pass it on to a solver. Accompanying tools help the user in common tasks such as numerically checking whether the cost function agrees with its derivatives up to the appropriate order, approximating the Hessian based on the gradient of the cost, etc.
Manifolds in Manopt are represented as structures and are obtained by calling a factory. The manifold descriptions include projections on tangent spaces, retractions, helpers to convert Euclidean derivatives (gradient and Hessian) to Riemannian derivatives, etc. All the manifolds mentioned in the introduction work out of the box, and more can be added (shape space (Ring and Wirth, 2012), low-rank tensors (Kressner et al., 2013), etc.). Cartesian products of known manifolds are supported too.
Solvers are functions in Manopt that implement generic Riemannian minimization algorithms. Solvers log standard information at each iteration and comply with standard stopping critera. Extra information can be logged via callbacks and, similarly, user-defined stopping criteria are allowed. Currently, Riemannian trust-regions (based on (Absil et al., 2007)) and conjugate-gradients are implemented (with preconditioning), as well as steepest-descent and a couple derivative free schemes. More solvers can be added, with an outlook toward Riemannian BFGS (Ring and Wirth, 2012), stochastic gradients (Bonnabel, 2013), nonsmooth subgradients schemes (Dirr et al., 2007), etc.
An optimization problem in Manopt is represented as a problem structure. The latter includes a field which contains a structure describing a manifold, as obtained from a factory. Additionally, the problem structure hosts function handles for the cost function f and (possibly) its derivatives. An abstraction layer at the interface between the solvers and the problem description offers great flexibility in the cost function description.
As the needs grow during the life-cycle of the toolbox and new ways of describing f become necessary (subdifferentials, partial gradients, etc.), it will be sufficient to update this interface.
Computing f (x) typically produces intermediate results which can be reused in order to compute the derivatives of f at x. To prevent redundant computations, Manopt incorporates an (optional) caching system, which becomes useful when transiting from a proof-of-concept draft of the algorithm to a convincing implementation.
3 Example: the maximum cut problem Given an undirected graph with n nodes and weights w ij ≥ 0 on the edges such that W ∈ R n×n is the weighted adjacency matrix and D ∈ R n×n is the diagonal degree matrix with D ii = j w ij , the graph Laplacian is the positive semidefinite matrix L = D − W . The max-cut problem consists in building a partition s ∈ {+1, −1} n of the nodes in two classes such that , i.e., the sum of the weights of the edges connecting the two classes, is maximum. Let X = ss . Then, max-cut is equivalent to: 0, diag(X) = 1 n and rank(X) = 1.
Goemans and Williamson (1995) proposed and analyzed the famous relaxation of this problem which consists in dropping the rank constraint, yielding a semidefinite program. Alternatively relaxing the rank constraint to be rank(X) ≤ r for some 1 < r < n yields a tighter but nonconvex relaxation. Journée et al. (2010) observe that fixing the rank with the constraint rank(X) = r turns the search space into a smooth manifold, the fixed-rank elliptope, which can be optimized over using Riemannian optimization. In Manopt, simple code for this reads (with Y ∈ R n×r such that X = Y Y ): % The problem structure hosts a manifold structure % as well as function handles to define the cost % function and its derivatives (here provided as % Euclidean derivatives, which will be converted % to their Riemannian equivalent). problem.M = elliptopefactory(n, r); problem.cost = @(Y) -trace(Y' * L * Y)/4; problem.egrad = @(Y) -(L * Y)/2; problem.ehess = @(Y, U) -(L * U)/2; % optional % These diagnostics tools help make sure the % gradient and Hessian are correct. checkgradient(problem); pause; checkhessian(problem); pause; % Minimize with trust-regions, % a random initial guess and default options. Y = trustregions(problem); Randomly projecting Y as s = sign(Y * randn(r, 1)) yields a cut. The Manopt distribution includes advanced code for this example, where the caching functionalities are used to avoid redundant computations of the product LY in the cost and the gradient, and the rank r is increased gradually to obtain a global solution of the max-cut SDP (and hence a formal upperbound), following the procedure in (Journée et al., 2010). Abstract-Emerging ultrasound phased-array technologies will soon enable the acquisition of high-resolution 3D+T images for medical applications. Processing the huge amount of spatiotemporal measurements remains a practical challenge. In this work, dynamic ultrasound images are sparsely represented by a mixture of moving speckles. We model the shape of a speckle and its locally linear motion with a weighted multivariate Gaussian kernel. Parameters of the model are estimated with online Bayesian learning from a stream of random measurements. In our preliminary experiments with a simulated phantom of a moving cylindrical structure, the optical flow of speckles is estimated for a vertical line profile and compared to the ground truth. The mean accuracy of the linear motion estimate is of 93.53%, using only a statistically sufficient random subset of the data.

Introduction
Measuring strain in structural soft tissue, such as tendons, could be an invaluable diagnostic tool for physicians [1]. For example, in orthopedics, measuring strain in ligaments will assist the placement procedure of joint implants and ultimately decrease the number of follow-up procedures [2]. Recent developments in ultrasound (US) imaging allows for the acquisition of high resolution volumetric images. This technology enables researchers to leverage existing imaging techniques such as elastography [3] or model-based biomechanical simulations of deformation [4] to the problem of assessing strain in tendons. Unfortunatelly, these techniques are not yet capable of capturing local motion occurring in small scale soft tissues.
In this work, we aim at adapting the speckle tracking echocardiography (STE) technique [5] from cardiology to the prescribed situation. Therefore, we propose a novel method for directly estimating motion of a sparse selection of speckle shape models instead of using detection and explicit tracking [6]. Speckles are forming the characteristic granular appearance in US images, which is due to deterministic interference originating from sub-wavelength scattering sites [7]. We approximate speckles and their local linear motion with a mixture of weighted multivariate Gaussian kernels. A Bayesian method is used for parametric estimation of the mixture components, which are representing the optical flow that is traced by the trajectories of speckles. Preliminary results are shown for simulated one-dimensional vertical line profiles.

Data and Image Models
We collect echography measurements in a dense discrete spatiotemporal amplitude field A(x, t), with x ∈ {1, . . . , N x } and t ∈ {1, . . . , N t }. This matrix of real positive amplitude values Spatio-temporal data Continuous reconstruction is first normalized for every time frame t such that We approximate A(x, t) with a continuous and smooth conditional Gaussian mixture model G(x, t) having the form with the set of model parameters x , σ 2 t , σ xt k where g(x, t; Θ k ) can be decomposed as the product of a marginal (temporal) and the conditional (spatial) distributions: , where the parameters of the conditional distribution arē and the Gaussian kernel g(x; µ, σ 2 ) is defined by Before estimating the mixture parameters, we first sample A(x, t) by taking a subset of M N x × N t amplitude values

Reconstruction
The set Θ of all weights and parameters of the mixture model G(x, t) is now estimated by maximizing the log-likelihood iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
No displacement 1mm displacement 2mm displacement Figure 2: Three selected frames of the moving cylinder phantom. The total course of the cylinder was 2mm long, such that the measured global displacement of the center of mass was 1.308mm.
We use the non-iterative online expectation-maximization (EM) method [8,9] for estimating the parameters. Online EM alternates between performing an expectation (E-step) which evaluates membership probabilities for the log-likelihood using the current estimate of the parameters, and a maximization (Mstep), which updates parameters for maximizing the expected log-likelihood found in the E-step. The model is first initialized with a regular grid, then the two E-and M-steps are performed successively for every sample, converging to the final model.

E-step
When updating the model from the information carried by the m-th sample (x m , t m , a m ), we first compute the membership probability p k (x m |t m ), expressing the probability that the sample is drawn from the conditional distribution of the component k at given time t m using where the parametersμ andσ 2 x of the conditional distribution are calculated as previously described.

M-step
In the second step, the parameters are updated such that the likelihood of observing the new sample is improved by applying the following update rules: Initialization We initialize each component of a given index k with the following weight w k = N t /K and parameters 12 , σ xt = 0. As presented, the weights are systematically incremented for every new sample in the M-step. In order to adapt smoothly to new data, a "forgetting" mechanism should be used to adjust parameters to the recent sufficient statistics. We simply apply a multiplicative decay factor γ = 1 − 1/N t to every component such that the total mixture weight remains normalized to N t .

Experiment
For this paper we use the Field II program [7] in order to simulate RF-signals from a linear US probe consisting of 192 (64  active) elements. To this end a 3D phantom is generated containing a 10mm long high-intensity cylinder of 4mm radius. The whole numerical phantom consists of 5000 uniformly distributed and amplitude-weighted random scatterers for a region of 30×30×10mm 3 . This cylinder is then consecutively moved downwards by 8 consecutive steps of 0.25mm along the vertical axis, while the background is left unaltered.
For the simulation, the probe is positioned 15mm from the phantom emitting a 3.5Mhz pulse focused at a depth of 30mm. Afterwards the return signals is sampled at a frequency of 100Mhz. Using this method a total of 25 RF lines are created on which we perform log envelope detection. The images shown in Figure 2 show three temporal frames of the input data after downsampling along the axial dimension by a factor 10.
For the experiment we created a 2D optical flow image A(x, t) using the central vertical line profile line at all N t = 9 temporal steps which is shown in Figure 1. The number of amplitude measurements for the profile lines is N x = 2597, the number of mixture components was set to K = 16, and the average number of samples used for statistical estimations of each component was set to was set to N = 512 such that the total number of samples was M = N K = 8192.
Reconstructed line profiles in Figure 3 show that the linear motion of speckles covering the cylinder is successfully captured with 16 mixture components. The center of mass in the data shifted by 1.308mm between the start and end of the simulation. The reconstructed displacement was of 1.224mm, reaching a relative error of 6.47% in this experiment. Similar experiments have been conducted with 8 and 32 mixture components. The motion extraction was not accurate with only 8 components, while we obtained comparable results with 32 components. We observed that the simulated speckle pattern slightly deforms over time while following the displacement path of the object. This effect has not yet been accounted for.

Conclusion
We represent dynamic ultrasound images as a mixture of moving speckle shape models that are approximated by smooth weighted multivariate Gaussian kernels. Maximum likelihood values of these kernel parameters can therefore be estimated by a conditional expectation-maximization method, using simple closed-form update rules. Experiments have been conducted on a line profile of a simple cylinder phantom. We observed that the sparse reconstruction approximates smoothly the acquired data. The subtle motion of the moving structure on a static background was successfully captured in the model. Future work will consider modeling the space-variant point spread function of speckles and extend the method to 3D+T imaging for assessing results on clinical data.  Abstract-A powerful approach to sparse representation, dictionary learning consists in finding a redundant frame in which the representation of a particular class of images is sparse. In practice, all algorithms performing dictionary learning iteratively estimate the dictionary and a sparse representation of the images using this dictionary. However, the numerical complexity of dictionary learning restricts its use to atoms with a small support.
A way to alleviate these issues is introduced in this paper, consisting in dictionary atoms obtained by translating the composition of K convolutions with S-sparse kernels of known support. The dictionary update step associated with this strategy is a nonconvex optimization problem, which we study here.
A block-coordinate descent or Gauss-Seidel algorithm is proposed to solve this problem, whose search space is of dimension KS, which is much smaller than the size of the image. Moreover, the complexity of the algorithm is linear with respect to the size of the image, allowing larger atoms to be learned (as opposed to small patches). An experiment is presented that shows the approximation of a large cosine atom with K = 7 sparse kernels, demonstrating a very good accuracy.

Introduction
The problem we introduce in this paper is motivated by the dictionary learning (DL) field. DL was pioneered by [8,9] and has received a growing attention since then. The principle behind DL is to find a representation for data that makes it simpler, sparser. We invite the reader to consult [4] for more details about sparse representations and DL. The archetype of the DL strategy is to look for a dictionary as the solution of the following optimization problem where y i are the learning database, D is the dictionary matrix, whose columns are the atoms, and f is a sparsity-inducing function. The resulting problem can be solved (or approximatively solved) by many methods including MOD [5] and K-SVD [1]. All these approaches rely on alternatively updating the codes x i and the dictionary D.
Our primary motivation for considering the observation model (1) comes from computational issues. Usually, DL is applied to small patches, because of the computational cost of repeatedly computing of the matrix-vector products Dx i , which is worse than O(N 2 ). Moreover, the cost of the dictionary update is usually worse than O(N 3 ).
We propose a model where the learned atoms are a composition of K convolutions with S-sparse kernels. The interest for such a constraint is to provide numerically effective dictionaries and allow to consider larger atoms. Indeed, the search space is only of dimension KS, which is typically smaller than the size of the target atom.
The present work focuses on the dictionary update step of one atom. In this context, the code x is known. Our goals are both to approximate a large target atom κ with our model and to obtain target atoms whose manipulation is numerically efficient. This translates into a non-convex optimization problem.

Problem formulation
Let consider an observed d-dimensional signal y of (R N ) d , assumed to result from the convolution of a known input signal x ∈ (R N ) d with an unknown target kernel κ ∈ (R N ) d , contaminated by an additive noise b following the linear model where * stands for the circular discrete convolution 1 in dimension d. For instance, the unknown target kernel κ ∈ (R N ) d may refer to the unknown impulse response of a 1D-filter or, conversely, to the point spread function of a 2D-filtering operator.
The problem addressed in this paper consists of approximating the unknown kernel κ by a composition of convolutions The kernels h 1 , . . . , h K are constrained to have less than a fixed number S of non-zero elements, i.e., they are assumed to be at most S-sparse. As stated before, this assumption aims at providing a cost-effective dictionary by reducing the computations for x * κ. Furthermore, the locations of their non-zero elements in {0, . . . , N } d are assumed to be known or pre-set. More precisely, the support of the kth kernel (i.e., the locations of the non zero elements of h k ), denoted supp h k , is constrained to a fixed set of discrete indexes S k An example of indexes for 1D convolution kernel would be Assuming that the noise vector b is an independent and identically distributed Gaussian sequence, approximating the unknown convolution kernel κ from the observed measurements y can be formulated as the following optimization problem where · 2 stands for the usual Euclidean norm in (R N ) d . The problem (P 0 ) is non convex. Thus, depending on the values of K, (S k ) 1≤k≤K , x and y, it might be difficult or impossible to find a good approximation of a global minimizer of (P 0 ). Our objectives are to study if such a problem bends itself to global optimization, and to assess the approximation power of the computed compositions of convolutions.

Block-coordinate descent
The problem (P 0 ), as formulated in the previous section, is unhandy for global optimization. As detailed in [2], it has irrelevant stationary points and is non-convex (though infinitely differentiable). To adress these issues, a scalar weight λ is introduced and kernels are constrained to have a unit norm. Moreover, we elect a block-coordinate formulation in order to solve the problem with a Gauss-Seidel type algorithm (called Alternate Least Squares, sharing many similarities with the one used in [7]).
where the kernels h k are fixed ∀k = k. The problem (P k ) is linear and can be expressed as a matrix-vector product considering only the elements of h that belong to its support: The idea is to alternatively solve (P k ) by iterating on k. The support constraint significantly reduces the search space of the problem, and thus the amount of calculations needed to solve it for a stationary point. Algorithm 1 shows an overview of the resolution of (P k ). The computational complexity associated with a passage in the while loop is O((K + S)KSN d ), i.e., it is linear with respect to the size N d of the signal. The detailed steps to solving (P k ) are given in [2].
Algorithm 1: ALS algorithm Input: y: target measurements; x: known coefficients; (S k ) 1≤k≤K : supports of the kernels (h k ) 1≤k≤K . Output: (h k ) 1≤k≤K : convolution kernels such that begin Initialize the kernels ((h k p ) p∈N ) 1≤k≤K ; while not converged do for k = 1 ,. . . , K do Update h k and λ with a minimizer of (P k ) ;

Synthetic example
In this section, we show an experiment consisting of approximating a 2D cosine atom κ in an image y of size 64 × 64 (i.e., d = 2 and N = 64). Such an atom can be seen as a large local cosine or a Fourier atom, both widely used in image processing. The interest of this atom is that it covers the whole image and is of a rather large support, making it difficult to handle with existing dictionary learning strategies.
The code x is a sparse image whose elements are chosen independent and identically distributed according to a Bernoulli-Gaussian distribution, widely used in sparse signal and image deconvolution [3,6,10]. Therefore, y contains a few weighted translations of the cosine atom κ 2 . The target y is built with additive Gaussian noise of variance σ 2 = 0.5. Kernel supports have been set to a simple 5 × 5 square, linearly dilated with k, similar to the 1-D example given in (4). Figures 1 and 2 show the cosine image y, its approximation λx * h 1 * · · · * h K , the actual atom κ and λh 1 * · · · * h K , for K = 7 and S = 25. The results obtained here are quite accurate even though the cosine image was corrupted by additive noise.

Conclusion
This work shows that simple atoms can be accurately approximated with a composition of convolutions. The kernels used in the approximation are constrained to be sparse (i.e., with sparse supports), leading to a computationally efficient algorithm, despite the non-convexity of the function to optimize. This efficiency was illustrated on a 2D-cosine function, but similar experiments conducted with archetypal kernels (e.g., wavelets or curvelets) show similar performances [2].
The proposed modeling and algorithmic schemes open new perspectives on the general problem of dictionary learning. More specifically, it seems reasonable to derive a DL technique which recovers large structured dictionary whose atoms consist of compositions of convolutions.
Finally, how to choose, set or draw the kernel supports remains a large and yet unexplored issue, that may have significant impact on the method performances.
Abstract-We consider the problem of a set of nodes which is required to collectively learn a common dictionary from noisy measurements. This distributed dictionary learning approach may be useful in several contexts including sensor networks. Diffusion cooperation schemes have been proposed to estimate a consensus solution to distributed linear regression. This work proposes a diffusion-based adaptive dictionary learning strategy. Each node receives measurements which may be shared or not with its neighbors. All nodes cooperate with their neighbors by sharing their local dictionary to estimate a common representation. In a diffusion approach, the resulting algorithm corresponds to a distributed alternate optimization. Beyond dictionary learning, this strategy could be adapted to many matrix factorization problems in various settings. We illustrate its efficiency on some numerical experiments, including the difficult problem of blind hyperspectral images unmixing.

Introduction
In a variety of contexts, huge amounts of high dimensional data are recorded from multiple sensors. When sensor networks are considered, it is desirable that computations be distributed over the network rather than centralized in some fusion unit. Indeed, centralizing all measurements lacks robustness -a failure of the central node is fatal -and scalability due to the needed energy and communication resources. In distributed computing, every node communicates with its neighbors only and processing is carried out by every node in the network. Another important remark is that relevant information from the data usually lives in a space of much reduced dimension compared to the physical space. The extraction of this relevant information calls for the identification of some adapted sparse representation of the data. Learning an adaptive sparse representation of the data using a redundant dictionary is useful for many tasks such as storing, transmitting or analyzing the data to understand its content, see [1] for an up-to-date review. Furthermore, the problem of dictionary learning belongs to the more general family of matrix factorization problems that appears in a host of applications. We study the problem of dictionary learning distributed over a sensor network in a setting where a set of nodes is required to collectively learn an adaptive sparse representation from independent observations. We consider the situation where a set of connected nodes records data from observations of the same kind of physical system: each observation is assumed to be described by a sparse representation using a common dictionary over all sensors. For instance, a set of cameras observe the same kind of scenes or a set of microphones records the same kind of sound environment.
The dictionary learning and the matrix factorization problems are connected to the linear regression problem. Indeed, the classical approach based on alternate minimization on the coefficients X and the dictionary D solves two linear regression problems knowing respectively D or X. Several recent works have proposed efficient solutions to the problem of least mean square (LMS) distributed linear regression, see [2] and references therein. The main idea is to use a so-called diffusion strategy: each node n carries out its own estimation D n of the same underlying linear regression vector D but can communicate with its neighbors as well. The information provided to some node by its neighbors is taken into account according to weights interpreted as diffusion coefficients. Under some mild conditions, the performance of such an approach in terms of mean squared error is similar to that of a centralized approach [3]. Let D c the centralized estimate which uses all the observations at once. It can be shown that the error IE D n − D 2 of the distributed estimate is of the same order as IE D c − D 2 : diffusion networks match the performance of the centralized solution.
Our work [4] gives strong indication that the classical dictionary learning technique based on block coordinate descent on the dictionary D and the coefficients X can be adapted to the distributed framework by adapting the diffusion strategy mentionned above. Our numerical experiments also strongly support this idea. Note that solving this type of matrix factorization problems is really at stake since it corresponds to many inverse problems: denoising, adaptive compression, recommendation systems... A distributed approach is highly desirable both for use in sensor networks and for parallelization of numerically expensive learning algorithms. In a second step, we may consider the more general situation where observations may also be shared between connected nodes.

Problem formulation
Many nodes, one dictionary. Consider N nodes over some region. In the following, boldfaced letters denote column vectors, and capital letters denote matrices. The node n takes q n measurements y n (i), 1 ≤ i ≤ q n from some physical system. All the observations are assumed to originate from independent realizations s n (i) of the same underlying stochastic source process s. Each measurement is a noisy measurement where z denotes the usual i.i.d. Gaussian noise with covariance matrix Σ n = σ 2 n I. Our purpose is to learn a common redundant dictionary D which carries the characteristic properties of the data. This dictionary must yield a sparse representation of s so that: ∀n, y n (i) = Dx n (i) where x n (i) features the coefficients x nk (i) associated to the contribution of atom d k , the k-th column in the dictionary ma-iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. (13) trix D, to s n (i). The sparsity of x n (i) means that only few components of x n (i) are non zero. We consider the situation where a unique dictionary D generates the observations at all nodes. On the contrary, observations will first not be shared between nodes (this is one potential generalization). Our purpose is to learn this dictionary in a distributed manner thanks to in-network computing only. As a consequence, each node will locally estimate a local dictionary D n thanks to i) its observations y n and ii) communication with its neighbors. The neighborhood of node n will be denoted by N n , including node n itself.
Dictionary learning. Various approaches to dictionary learning have been proposed [1]. Usually, in the centralized setting, the q observations are denoted by y(i) ∈ R p and grouped in a matrix Y = [y(1), ..., y(q)]. As a consequence, Y ∈ R p×q . The dictionary (associated to some linear transform) is denoted by D ∈ R p×K : each column is one atom d k of the dictionary. The coefficients associated to observations are X = [x(1), ..., x(q)]. We will consider learning methods based on block coordinate descent or alternate optimization on D and X with a sparsity constraint on X [5,1,6]. The data is represented as the sum of a linear combination of atoms and a noise term Z ∈ R p×q : In the most usual setting featuring white Gaussian noise, one wants to solve: Under some mild conditions, this problem is known to provide a solution to the underlying L0-penalized problem [7].
3 Distributed alternate optimization for dictionary learning Algorithm. The Adapt Then Combine diffusion strategy [2] for distributed estimation originates the following approach to distributed alternate optimization for dictionary learning. Diffusion is ensured by the communication between nodes sharing their dictionary estimate with neighbors in N n . Observations are taken simultaneously at each node so that a whole data matrix Y n is assumed to be available at node n. Here index i stands for iterations. The case where data arrive sequentially at each node can also be dealt with at the price of a natural adaptation of the present approach. Each node must estimate both its local dictionary D n and the coefficients X n which describe observations Y n = D n X n + Z n . At each iteration i, only the local dictionary estimates D n,i are assumed to be shared between neighbors, not observations. In summary, sparse representations are computed locally. Then each node updates its dictionary as a function of its local observations Y n (Adapt step) and its neighbors' dictionaries (Combine step). Based on known results for the ATC strategy in its usual setting, we expect that Algorithm 1 below converges to an accurate estimate of the common underlying dictionary D. Various choices can be considered for A such as some a priori fixed matrix A or with the relative degree variance (ν = degree of node ): Initialize D n,0 , ∀n (random subset of K observations y n (i)). Given a matrix A satisfying Repeat until convergence of (D n,i , X n,i ) n=1:N For each node n repeat: 1) Optimization w.r.t. X n,i (sparse coding): Given the dictionary D n,i , the coefficients X n,i are estimated using a sparse coding method (Basis Pursuit, OMP, FOCUSS,...) 2) Optimization w.r.t. D n,i (dictionary) e.g. by gradient descent: can also be updated by MOD or K-SVD,...at node n.
Numerical experiments. We present some numerical experiments to illustrate the relevance and efficiency of our approach. For instance we show the results obtained on a dataset built from a redundant random dictionary of 48 atoms of dimension 16 corresponding to image patches of size 4 × 4. Each data y n (i) is the linear combination of 3 atoms with i.i.d. coefficients uniformly distributed over [−0.5, 0.5] ; various Gaussian noise levels are considered. We show that a set of 4 nodes in a symmetrically connected network consistently learn the same dictionary of 4 × 4 patches with good accuracy (45 atoms of the initial dictionary of 48 atoms are recovered with We will also show the results of an application to the problem of blind unmixing of hyperspectral images. In this application, the network is simply made of connections between pixels which are spatially close or which carry similar spectral information. The graph underlying (hyper-)pixels makes them collaborate to learn spectral endmembers.

Conclusion
We present an original algorithm which solves the problem of distributed dictionary learning over a sensor network. This is made possible thanks to a diffusion strategy which permits local communication between neighbors. Connected nodes exchange their local dictionaries estimated from disjoint subsets of data. This algorithm adapts usual dictionary learning techniques for sparse representation to the context of in-network computing. This approach to the general problem of distributed matrix factorization paves the way towards many prospects and applications. Moreover, as far as computational complexity is concerned, distributed parallel implementations are a potentially interesting alternative to online learning techniques [8].
We may even consider a dynamical context where observations arrive over time so that the dictionary would also be learnt dynamically.
Abstract-We propose a semi-supervised crack detection method that can be used for high-dimensional and multimodal acquisitions of paintings. Our dataset consists of a recent collection of images of the Ghent Altarpiece (1432), one of Northern Europe's most important art masterpieces. We build a classifier that is able to discern crack pixels from the background consisting of non-crack pixels, making optimal use of the information that is provided by each modality. To accomplish this we employ a recently developed non-parametric Bayesian classifier, that uses tensor factorizations to characterize any conditional probability. A prior is placed on the parameters of the factorization such that every possible interaction between predictors is allowed while still identifying a sparse subset among these predictors.

Introduction
The digitization of cultural artifacts is becoming an increasingly popular practice. Consequently, their digital analysis emerged as a key non-invasive tool for art historians and preservation experts. In this work, we focus on crack detection in the Ghent Altarpiece (1432). Being able to accurately detect cracks is of major importance in the context of art conservation since cracking is one of the most common forms of deterioration found in paintings. For most 15 th century Flemish paintings on Baltic oak, humidity fluctuations cause the wooden support to shrink or expand and are the main cause for crack formation.
Previous efforts on crack detection in the Ghent Altarpiece can be found in [1] and an extension in [2], where three detection methods (filtering with elongated filters, a multi scale top hat transform and K-SVD for crack enhancement) are proposed. The images that we analyzed back then are scans of original photographic negatives. However, thanks to the project Lasting Support 1 (2010-2011), a new multimodal dataset of ultra high resolution was made available. The panels of the altarpiece were documented with digital macrophotography in both the visible and infrared parts of the electromagnetic spectrum, together with previously acquired X-rays.
This multimodal dataset requires new crack detection techniques that are able to combine the information provided by each modality. Simply applying the methods from [2], which were designed to work well on one single modality, would require choosing an additional set of parameters per modality. The goal is to build an accurate classifier and select a subset of most important predictors from a a high dimensional collection of categorical predictor vectors. We define a model, by using a carefully-structured Tucker factorization, that can characterize any conditional probability and facilitates variable selection and the modeling of higher-order interactions [3]. We follow a Bayesian approach and use a Markov chain Monte Carlo (MCMC) algorithm for posterior computation and hence allow for uncertainties in the to be included predictors. The proposed approach which we refer to as Bayesian conditional tensor factorization or BCTF, introduced in [4], shows comparable detection results to the Random Forest (RF) algorithm [5] but has the tendency to produce finer crack maps and is faster, which is an attractive property when working with ultra high-dimensional data.

Data preprocessing
The available modalities of the painting are fundamentally very different and therefore direct registration is a challenging task. Since cracks are the only components that are more or less consistent throughout the modalities, crude crack maps are initially constructed. These are obtained by first filtering the image with elongated filters [2] and a subsequent thresholding step produces binary crack maps that show an initial estimate of the crack locations. These are registered using the algorithm described in [6] and the resulting transformation is then applied to the original images.
The penetrating nature of X-rays make them an appropriate modality for the detection of cracks, since they usually appear clearly delineated. However, the presence of the wood grain in some of the panels can act as a confounding factor during the detection process. We use a texture separation technique, called morphological component analysis (MCA) to separate the texture from the natural part in images. The textural part in our images consists of the wood grain, for which we select a (Local) Discrete Cosine Transform (DCT) dictionary. The Dual Tree Complex Wavelet dictionary is used for the representation of piecewise smooth content.

Image features
Each modality is processed and filtered with a wide range of filters commonly used in image processing. The methods and filters that are used to construct the image predictor vectors include elongated filters, a Frangi Vesselness filter [7], the eigenvalues of a structure tensor, black Top Hat transforms with different structuring elements [8], Local Binary Patterns [9], pixel intensities, median filters, Laplacian of Gaussian (LoG) filters and the Leung-Malik filter bank, a multi-scale, multiorientation filter bank of 48 filters. Each pixel is now characterized by a feature vector of length 208 that is the result of all the operations described above, applied to each modality independently. Once these are extracted they are quantized into an experimentally chosen number of bins d and used in the subsequent classification step, described below.

Bayesian Conditional Tensor Factorizations
We use a recently developed framework for nonparametric Bayes classification through tensor factorizations of the conditional probability P (Y |X), with Y a categorical response and X = (X 1 , . . . , X p ) a vector of p categorical predictors [3]. In our case the values of Y are either 0 or 1, meaning a pixel belongs to the background or it is part of a crack. The proposed method, which is called Bayesian conditional tensor factorization or BCTF, is based on the principle that a conditional probability can be expressed as a multidimensional array, i.e. a d 1 × · · · × d p dimensional tensor A = [a i1,··· ,ip ] d1×···×dp , with d j , the number of quantization bins of the jth predictor X j and This conditional probability tensor has non-negative cells, i.e. a i1,··· ,ip ≥ 0 for all i 1 , . . . , i p . To determine the tensor values we perform a low rank tensor factorization, similar to higher-order singular value decomposition (HOSVD), introduced in [10]. Practically, computing an HOSVD of an nth-order tensor leads to the computation of n different matrix SVDs. The HOSVD decomposes the tensor A = [a i1,··· ,ip ] d1×···×dp as: i,j ] dm×dm are orthogonal matrices and S = [s j1,··· ,jp ] d1×···×dp is a so called all-orthogonal core tensor. In [3] it was shown that any conditional probability tensor can be decomposed in a similar (non-unique) fashion where moreover holds for every combination of (m, x m ) and where we assume that the values of k m ∈ {1, · · · , d m } are chosen as small as possible. The factorization coefficients λ j1,j2,··· ,jp (y) can be seen as the latent class allocation probabilities and π jm (x m ) as the response class probabilities and are non-negative.
A primary goal is reducing the enormous amount of covariates, as we expect that only a few features have a significant impact on the classification results, by imposing sparsity. The k m value impacts the number of parameters used to characterize the mth predictor as well as sparsity. In the special case in which k m = 1 and when taking into consideration the constraints in (2), this results in π (m) 1 (x m ) = 1, which means that P (y|x 1 , · · · , x p ) will not depend on x m and the mth predictor is excluded from the model. In other words, sparsity can be imposed by setting k m = 1 for most m's. Additionally, k m can be seen as the number of latent classes for the mth covariate. The levels of X m are clustered according to their relationship with the response variable in a soft-probabilistic manner, where k m controls the complexity of the latent structure.
To complete our Bayesian model we choose independent Dirichlet priors (commonly used in Bayesian statistics) for the parameters Λ = {λ j1,··· ,jp ; j m = 1, · · · , k m ; m = 1, · · · , p} and π = {π These priors have the advantages of imposing the nonnegativity and summing to one constraints. The hyper parameters in the Dirichlet priors are chosen to favor placing most of the probability on a few elements, including near sparsity in these vectors. To embody our prior belief that only a small number of k j 's are greater than one we set P (k j = 1) = 1 − r p and P (k j = k) = r (dj −1)p for k = 2, · · · , d j , j = 1, · · · , p, and where r is the expected number of predictors to be included. To further impose sparsity, we include a restriction that #{j : k j > 1} ≤r, wherer is a prespecified maximum number of predictors. The effective prior on the k j 's is where I A (·) is the indicator function for set A. The full conditional posterior distributions of Λ and π are of the same family as the prior distributions.  Next to the precise detection, the BCTF classifier has a number of appealing properties. For one, there is no parameter tweaking involved which is usually a time consuming process when working with more common image processing operations. The second interesting fact is that the output of the BCTF is actually the probability for each pixel of being a crack pixel. This means that we can propose a probability crack map instead of a binary crack map as a solution. A binary crack map is obtained by thresholding the corresponding probability map. A pixel is considered to be part of a crack when it has a crack probability of 0.5 or higher. Moreover, the BCTF method is significantly faster than the RF classifier, a patch of size 256 × 256 is processed within 10 seconds with BCTF while it takes 18 seconds on average with a RF classifier. Abstract-Perfusion MRI is a method in medicine used mainly for diagnosing carcinoma and cardiovascular diseases. In the method a contrast agent is injected in the patient and then its concentration is observed via MRI during time. The signal(s) captured from the affected area can be described by the curve of lognormal distribution. The standard way of obtaining the measurements is very slow and does not comply with today's challenging requirements. We propose using compressed sensing to acquire much less coefficients, having minimal effect on the signal reconstruction (using the assumption that the data is a sum of low-rank matrix and sparse matrix in row spectrum).

Introduction
Perfusion MRI Perfusion imaging of the organism is an important tool in the diagnosis of cardiovascular and oncological diseases. Currently perfusion imaging may be used in a variety of imaging modalities such as magnetic resonance imaging, computed tomography, or ultrasound.
General procedure of perfusion imaging is similar for all modalities. A suitable contrast agent is incorporated into the organism as an injection or infusion and acquisition of data is performed in area of interest. Contrast agent is distributed in the body due to the cardiovascular system and its temporal and spatial distribution can be monitored.
Areas of interest may be different organs, parts of organs or individual voxels. The time course of concentration contrast agent in the region of interest is called the perfusion curve. By the analysis of perfusion curves we obtain estimates of perfusion parameters necessary for making the diagnosis.
Perfusion curve can be mathematically modeled using the density function of lognormal distribution with appropriate parameters Here c represents the offset, t 0 is delay of bolus between injection point and area of interest, S is area under curve and µ and σ is parameters of lognormal distribution. Pefusion curve is sparse in the Fourier spectrum ( Fig. 1) which we use later. Compressed sensing Measurements in compressed sensing are linear combinations of data, and the reconstruction is nonlinear. In MRI, the signal is captured in the frequency domain, therefore a single "measurement" in this context is obtaining one Fourier coefficient. Compressed sensing helps accelerating the measurement process but this is for the price of the time for numerical reconstruction. It starts from the apriori knowledge of the signal.

Model specification Preparation of data
To be able to use compressed sensing, the MRI video is reshaped into a matrix, where each column of this matrix represents a single video frame (Fig. 2). The perfusion curves form the rows of this matrix.

Measurement process
Data is measured within each frame's spectrum (k-space) by application of mask (matrix consisting of ones/zeros). We use several types of mask (see Fig. 3) to measure: • completely random with uniform distribution • random with nonuniform distrubution • radial (randomized or not) • spiral (randomized or not) Signal recovery by L + S model We use so-called L+S model for reconstruction perfusion data. In this model, the video data M is decomposed as a superposition of a low-rank matrix L (few non-zero singular values) and a sparse matrix S (few non-zero entries). See the optimization problem in Fig. 4. For better reconstruction, the Fourier coefficients which are unmeasured but have a complex-conjugated measured counterpart, are completed in this simple manner.
iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. Numerical solution by proximal gradient method Proximal gradient method is for solving this type of convex minimization problem: where g, h are convex and g is smooth. This problem can be solved iteratively by two steps: where prox is the so-called proximity operator. Applied to our problem, function g represents the quadratic term of the objective function and h represents the nuclear and 1 norms. Then the respective proximity operators are the soft thresholding of singular values for nuclear norm and soft thresholding of coefficients for 1 norm.

Results
We create phantom of perfusion data of size 100 × 100 × 100 and 150 × 150 × 100 pixels in Matlab. Then we measured this video using a single-type mask, but with different random parameters in each frame and we reconstructed the video by the L + S algorithm. See Fig. 5. An example of the original 2D matrix, its reconstruction M and its decomposition into L and S for the random mask is depicted in Fig. 6.

Conclusion
Decomposition of the perfusion MRI matrix using the L + S model can rapidly reduce number of measurements. The best performance is obtained by using the nonuniform random mask. When the data size is bigger the results are better.

Introduction
Hyperspectral imaging is an advanced imaging technique which integrates spectroscopy into the image capturing process, providing for each pixel its light intensity as a function of wavelength. The resulting information is organized into a hyperspectral (HS) cube with two spatial and one spectral dimensions, which may be interpreted as a stack of images, one for each wavelength. Spectral cameras are currently used in research applications such as remote sensing [1] and food inspection [2], where this detailed spectral information may be used at application-level to identify objects and materials. Due to rapidly falling system complexity, size and cost, such technologies are expected to be adopted in more general purpose applications [3]. In order to acquire these vast amounts of data, spectral cameras typically use a time-consuming scanning approach, based on, e.g., line scanning with dispersive optics (such as prisms and gratings [4]) or spectral scanning using tunable filters (e.g., AOTF [5]). This problem may be overcome using snapshot acquisition, where the entire 3-D datacube is acquired during one frame period by optically multiplexing the contents of the 3-D cube onto a 2-D sensor. Whereas in scanning spectral imaging the data in a spectral cube is conveyed to multiple, consecutive frames, in snapshot spectral imaging this data is multiplexed onto a single frame of the sensor. However, due to the limited number of sensor pixels, snapshot imagers require a trade-off between spatial and spectral resolution. This trade-off has been previously tackled by using compressive sensing (CS) techniques [6] to acquire the 3-D HS cube at sub-Nyquist sampling rates [7,8]. In this work, we introduce a new compressive HS scheme that combines a snapshot imaging sensor based on monolithically integrated Fabry-Pérot [9] filter tiles with an optical front-end based on an out-of-focus spatial light modulator (SLM). By allowing the number of SLM patterns to vary at runtime, the system can reach higher SNR, thus operating as a hybrid setup between a scanning and a snapshot spectral camera. Since compressive sensing techniques avoid scanning, associated problems with motion blur and reduced SNR in high speed applications are also avoided.

Optical and Sensing Model
In this work, we extend the monochromatic compressive optical scheme developed by Björklund and Magli [10] for acquir-We thank the supercomputing facilities (CISM/UCL) and the CECI, funded by FRS-FNRS, for their computational resources. LJ and KD are funded by the FRS-FNRS. ing hyperspectral images. Our imager model is actually built around a snapshot spectral sensor, which monolithically integrates a set of Fabry-Pérot (FP) interferometers organized as tiled filters on top of a standard, off-the-shelf CMOS sensor [3] (Fig. 1). The FP filter is typically made of a transparent layer (or cavity) with two mirrors at each side of that layer. The cavity length and the mirror reflectivity determine the selected wavelength and the spectral bandwidth (or full width at half maximum) of the filter, respectively. Thanks to the monolithic integration of the filter on the sensor, the stray light in the system is heavily reduced and the sensitivity and the speed of the sensor are increased [11]. The use of CMOS process technology heavily reduces the cost and improves the compactness of the hyperspectral camera. Let us now explain how the FP filters are combined with the compressive setup. The optical model of the imager, restricted to a 2-D section for simplicity, is illustrated in Fig. 2 (a). First, in a monochromatic setting, an object is imaged by conventional objective optics O on an image plane. Classical optics tells us that every point on this plane radiates light in many directions, as initially produced by the original object point. Therefore, after the image plane, the object image is replicated in many beams of parallel light rays. By placing an out-offocus spatial light modulator (SLM) after that plane, each of such beams can be modulated differently according to its direction, before being focused (i.e., "summed") on a single pixel of the detector by a second lens F . If this lens has focal length f , this pixel is located at distance f tan θ of the optical axis. By construction, the detector actually records specific samples of the convolution (up to a kernel reversal) between the image and the modulation. In the case of polychromatic imaging with no chromatic aberrations in our optics 2 , by tiling the detector with FP spectral filters we independently convolves different slices of the HS volume with the SLM modulation. As explained hereafter, by inserting a random SLM pattern, the optical scheme can compressively image the HS volume, wavelength-wise, with a partial Toeplitz sensing [12].
Mathematically, let a HS volume x ∈ R N ×N ×L 2 with N tot = L 2 N 2 voxels to be acquired by the imager. The number of wavelengths is set to a square value L 2 since the detector is made of a 2-D grid of N s × N s pixels covered by a grid of L × L square tiles T j (1 j L 2 ) of distinct FP filters, each selecting one wavelength λ j . Assuming N s divisible by L, each tile spectrally filters the light received by a square patch of M λ × M λ pixels, with M λ = N s /L. By design (see Fig. 2 (a)), since any replication of the image must be spatially modulated at resolution N (so as to determine the spatial resolution of x) and by adjusting a shift of one pixel in the detector to a shift of one pixel for the associated image beam in the SLM plane, For limiting the SLM resolution, we arbitrary set N slm = 2N − 1 so that N s = N and M λ = N/L. Every pixel of a given tile T j is influenced by a limited part ϕ j ∈ R (N +Mλ)×(N +Mλ) of the whole SLM modulation ϕ ∈ R Nslm×Nslm . Therefore, denoting by w the vectorization of any matrix w, the recorded light intensities is a selection matrix extracting indices associated to pixels in T j , and n j accounts for possible measurement noise. Gathering all the wavelengths and vectorizing the result, the whole sensing model is y = Φx + n, with Φ := diag(Φ 1 , · · · , Φ L 2 ) and y, n ∈ R Mλ×Mλ×L 2 . The final sensing is thus realized independently for each wavelength, as observed from the block structure of Φ. This restricts the performances of our system compared to an ideal compressive imaging device, i.e., when Φ is dense and unstructured. However, conversely to the sensing operator, the reconstruction methods recovering x from y as described in Sec. 3 exploit the 3-D structure of the data. Compressed Sensing theory [6] shows that if a M × D sensing matrix A respects the restricted isometry property (RIP) with sufficiently small RIP constant, any K-sparse (or compressible) signals of R D can be recovered (or estimated) from the possibly noisy measurement vector Ax ∈ R M . Moreover, a structured sensing issued from a M -subsampled Toeplitz or circulant matrix whose generating entries are drawn from an iid Bernoulli ±1 distribution, satisfies the restricted isom- [12]. Since the matrices Φ j respect such a Toeplitz structure (even in 2-D), we can therefore expect to efficiently reconstruct each slice of the HS volume x under a sparse signal prior provided the number of observations is large. However, this work tests the regularizing capabilities of two other priors: a 3-D total variation (TV), for a "cartoon shape" HS model, and a low-rank and joint sparse prior [13]. We also consider a situation where the number M 2 λ = N 2 /L of observations per band for one modulation ϕ is fixed by the detector design. An increasing number of observations is reached by merging the measurements obtained for S different random SLM patterns. The HS volume x of N 2 L 2 voxels is thus observed with a total of M 2 λ L 2 S = N 2 S measurements, i.e., we aim at reaching 1 S L 2 with high reconstruction quality.

Simulations
The capability of our HS imaging scheme has been simulated as follows. We focus our study on two toy examples of size N = 256 and L = 8, i.e., two structured HS volumes with N tot = 2 22 voxels. The first is a "Mondrian-like" example [15] with small 3-D total variation (TV) [14] made of 8 randomly generated overlapping cubes of different sizes and values in the HS volume. The second example follows a low-rank source model (LRSM) of rank r = 4 with (joint) sparsity level K = N 2 /16 = 4096 [13] (see Fig. 1). Two kinds of reconstruction methods have been used to recover these HS volumes from their compressive observations. For the "Mondrian" example, we solve x * = argmin u u T V +ρ y− Φu 2 2 thanks to the TwIST algorithm [16] and by manually adjusting ρ. The recovery of the LRSM example is achieved by promoting a low-rank and joint sparse data model [13], i.e., where A * is the trace norm of A, Ψ T is the spatial Haar wavelet transform [17] at every wavelength, and the first term in the minimized cost is the 1 / 2 mixed norm enforcing joint sparsity at every wavelength [18]. This optimization problem is solved by proximal algorithms [19,20] with ρ = K/2r [13]. Fig. 2 (b) and Fig. 2 (c) display the reconstruction SNR for both methods as a function of the number of random SLM patterns S. Two noisy sensing scenarios were considered with measurement noise at 20 and 80 dB. The sensing associated to our out-of-focus modulation with FP filters (OFMFP) is confronted with the more efficient spread-spectrum random Fourier ensemble (SSRFE) sensing [21], also applied bandwise. For the two noise levels and the two considered HS examples, we see that our sensing quickly reaches similar reconstruction qualities as SSRFE for increasing S, with a stronger match for the LRSM model. Moreover, for these two toy examples and under a 20 dB noise, as soon as S 16 (i.e., 25% of Nyquist rate), the SNR starts to saturate to a value related to the noise power.

Discussions and Conclusion
The capability of a compressive HS imager combining a CMOS detector tiled with FP filters and an out-of-focus SLM modulation has been studied. Despite an independent sensing of each wavelength, our simulations demonstrate already good performances on toy examples for a small number of modulations compared to the number of wavelengths. In a future work, the optical aberrations of the different optics and the SLM diffraction will be integrated in the sensing model. We will also study how to modulate the cube both spatially and spectrally. Abstract-Orthogonal Matching Pursuit (OMP) algorithms enable one to recover the support of sparse vectors on the basis of indirect measurements provided that the measurement process satisfies incoherence properties. However, the problem of computing this support in a context where several measurement vectors are available is addressed to a lesser extent in the literature. This paper proposes two computationally simple approaches based on OMP to deal with these situations. Finally, the proposed algorithms are evaluated on the basis of two different problems.

Introduction
This paper examines how to recover sparse vectors β [k] from measurements y l [k]. The inverse problem to be solved is given by Equation (1).
Equation (1) is often simplified when dealing with specific applied problems.
The sparse vectors β [k] ∈ R m to be estimated are assumed to remain unchanged (the dependency upon k then disappears) or to undergo limited modifications when k varies. l [k] ∈ R n is the noise added to the measurements. The entries of refers to different sensing matrices as measurements can be obtained from different matrices X l within some particular frameworks. This paper briefly describes the two proposed algorithms and then presents the envisioned applied problems used to assess their performance.

Exponential Decay Averaging (EDA) algorithms
Both algorithms average the output of the regular Orthogonal Matching Pursuit (OMP) algorithm [1]. EDAOMPVote uses the supports returned by OMP to provide an estimate of the support while EDA-β uses the reconstructed sparse vectors provided by OMP to produce an estimate of β.
The developed strategies are easily extended to other reconstruction algorithms (e.g. 1 norm algorithms, other greedy algorithms) since the only information that is needed for the averaging procedure is the support or the reconstructed vector β regardless of the method that is used to obtain them.

.1 Noise filtering
The first problem to be addressed is that of the noise filtering described by Equation (2). (2) iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. (24) k represents a discrete time instant and β varies slowly over time. EDA algorithms can be successfully applied to that kind of problems although numerical simulations have revealed that simpler approaches yield better results. The forgetting factor g is usually chosen according to the rate of variation of β. Nodes cooperate and exchange information with the Fusion Center (FC) so as to provide an estimate of β. The FC only uses low-complexity algorithms to perform the data fusion while the nodes perform the time-consuming operations (i.e. running OMP). The scalability of the system with regards to the number of nodes is therefore promoted.

Joint estimation
The problem to be solved is given by Equation (3). The time index k has disappeared while the index l refers to the different nodes with possibly different sensing matrices X l .
EDA algorithms (with g = 1) prove to be useful for solving joint estimation problems.

Noise filtering
Numerical simulations have revealed that EDA algorithms are not well suited to Problem (2) as it is possible to obtain better results by directly filtering the measurements y [k] and then run OMP. The filter for filtering the measurements has been chosen identical to that of OMP (i.e. o [k are respectively the output and the input of the filter at time k).

Joint estimation
Numerical simulations have been carried out to determine how the proposed algorithms compare to centralized estimation architectures.
In centralized estimation, all the computationally intensive operations are performed by the FC and the nodes are merely sensors that transmit raw data to the FC. Hence, the performance that can be achieved is greater but the FC computational capabilities should be high and the scalability of the system is thereby limited.
The proposed simulation setup consists of several nodes with identical sensing matrices. The nodes send either estimations of the sparse vector β or estimations of its support depending of the algorithm that is used (i.e. EDA-β or EDAOMPVote).
Without entering too much into the details, the sensing matrix X belongs to R 250×1000 . The sensing matrix is the same for all the simulation cases and its entries have been generated on the basis of i.i.d. Gaussian random variables.
OMP is configured to choose exactly s = β 0 atoms and computes the estimates accordingly.
For comparison purposes, centralized OMP algorithms have been run on the same cases. Simultaneous OMP (S-OMP) and OMP Multiple Measurement Vectors (OMPMMV) are two similar centralized OMP algorithms. Detailed information about S-OMP and OMPMMV are available in [2,3,4]. Within the framework of centralized estimation schemes, the nodes send the raw measurements y l = X l β + l to the FC.  Figure 2 shows that S-OMP and OMPMMV deal more efficiently with high numbers of nodes although the results achieved by S-OMP/OMPMMV and EDAOMPVote/EDA-β are comparable for low numbers of nodes. Moreover, EDA-β achieves slightly better performance than EDAOMPVote. Abstract-Most linear sparse representation algorithms can be straightforwardly extended to non-linear models. Their performance however, relies on an efficient computation of the gradient of the objective function. In this paper, we focus on a particular non-linear model, defined as the nested composition of functions and propose to resort to a well-known technique from the theory of optimal control to compute the gradient. As a proof of concept, this computation is then implemented into the optimization procedure proposed by Candès et al., and applied to a geophysical dynamical model.

Introduction
Recent contributions have addressed the problem of exploiting sparse priors with non-linear observation models, that is where h : R M → R N (with M ≥ N ) is a non-linear observation operator and n stands for an observation noise. Extending the approach followed in the linear case, these contributions propose also a generalization of the penalty function, leading to an optimization problem of the form (or some variants thereof) where J(x) is some scalar function (e.g., J(x) = y −h(x) 2 2 ) accounting for discrepancies from model (1).
Noticing that many sparse representation algorithms dealing with linear observation models rely -implicitly or explicitlyon the computation of the gradient of the function J(x), nonlinear versions of them can be straightforwardly derived. Following this idea, the extensions of the well-known algorithms IHT [1], MP [2], OMP [3] and CoSaMP [4] have thus been proposed, respectively in [5], [6], [7] and [8] .
However, whereas in the linear case, the evaluation of the gradient of J(x) only involves multiplications by the dictionary and its transpose, its computational cost can be prohibitive in some non-linear cases. In this paper, we elaborate on this problem for the particular family of cost functions J(x) defined as the nested composition of some functions. Formally, we write where {f l } L l=1 are some differentiable functions and • stands for the function-composition operator. This type of model is for instance of interest in the ubiquitous situations where one collects partial information on the state of a dynamical system whose initial condition admits a sparse decomposition (see section 2.2). In particular, we emphasize that results from optimal control [9] can be exploited to provide a fast implementation of any gradient-based algorithm by taking benefit of the special structure of the non-linear model (3). We propose then a practical implementation of this computation into the optimization procedure proposed in [10].

Sparse Representations in Nested
Non-Linear Models In this section, we elaborate on the efficient evaluation of the gradient of J(x) when structured as in (3). The methodology is then applied to a particular geophysical problem.

Efficient gradient computation
We use the following definitions and notations. Considering model (3), we set, ∀l ∈ {1, . . . , L}, ∀x ∈ R M , We thus have ∀l ∈ {1, . . . , L}, with the convention s 0 = x. We also define the gradient operator as so that ∇ x applied to a vector z = [z 1 , ..., z N ] T results in the M × N matrix whose (i, j)-th element is ∂zj ∂xi . With these notations in mind, (3) evaluated at x can be rewritten as where s l is defined as in (4) with x = x . Therefore, using the chain rule of derivative, we obtain and from the dependence between s l and s l−1 , Finally, applying this expression recursively, we have The latter expression is exclusively based on the derivative of each function component. Its evaluation can then be performed through the following forward-backward procedure: iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
• The sequence {s l } L l=1 is evaluated via the forward recursion • All multiplications by a same matrix ∇ sl−1 f l (s l−1 ) T are gathered in one single operation. This is done through the backward recursion leading finally to p 0 = ∇ x J(x ). In that way, the multiplication by each matrix ∇ sl−1 f l (s l−1 ) T is only performed once during the whole recursion.
This forward-backward procedure is widely used in geophysical applications (e.g., [11]). However, to the best of our knowledge, the explicit (and motivated) use of this technique into contexts of sparsity-constrained problems has never been considered. In particular, in [12] which focuses on a similar nonlinear model, this efficient computation of the gradient is not proposed.

Application to super-resolution in SQG dynamical model
As a practical example of the proposed methodology, we focus on a super-resolution problem in a geophysical context, namely the high-resolution characterization of the state of the ocean by exploiting: i) the Surface Quasi-Geostropic (SQG) dynamical model [13]; ii) a sparse prior on the initial condition of the dynamical model; iii) low-resolution satellite images. We assume that the SQG model, of the form of (5), rules the evolution of some state variable s l . The definition of f l depends on the considered numerical integration scheme (here a 4th-order Runge-Kunta method) but is not specified hereafter for conciseness. Moreover, as a prior knowledge, the initial state is supposed to be sparse in some redundant dictionary H, for some sparse vector x. We are interested in recovering the value of {s l } L l=1 from the observation of low-dimensional images {y l } L l=1 (e.g., collected by satellites), with where G l is some known observation matrix and n is an unknown corrupting noise.
In order to solve this inverse problem, we address the following optimization problem where λ > 0 and r(x) = m log(x m + ), = 10 −1 , is some sparsity-enforcing regularizing function. A solution to (16) is searched by using the majorization-minimization optimization technique exposed in [10]; at each iteration an upper bound on the goal function is constructed by majorizing r(x) by a weighted 1 norm. We look for the minimum of each of these majorizing functions by means of descent procedures involving the gradient of J(x) (corresponding here to the first term in (16)) evaluated as presented in section 2.1.
Particularized to the SQG model, the evaluation of the forward-backward recursions (10)-(13) have a complexity of order O(M L). By comparison, using a finite-difference scheme to evaluate the gradient requires to run (at least) two forward recursions by element of x, leading to an overall complexity of O(M 2 L). This order of complexity thus precludes us from using this type of approach in moderate-to-high dimensional problems.
The simulation setup considered in this paper is as follows. The state vectors s l are assumed to live in 256-dimensional space. The initial condition is supposed to have a sparse decomposition in a dictionary H ∈ R 256×512 made up of sine and cosine functions. The observations y l ∈ R 32 are collected at four different time instants and the observation matrices G l correspond to random subsampling operators. The ratio between the number of observations and the dimension of x is therefore equal to (32 × 4)/512 = 1/4. In Fig. 1, we represent the relative mean-square error (MSE) ||x − x|| 2 2 /||x|| 2 2 achieved by the minimization of (16) via the majorization-minimization procedure described above. As a point of comparison, we run the same algorithm on a linear sparse representation problem having the same problem dimensions (namely y = GHx where G is a rate-1/2 random subsampling matrix). For each data point, we average the performance over 50 trials.
We can notice that the considered procedure can achieve an acceptable relative mean square error over a wide range of sparsity levels. We note also that the non-linear setup suffers from a reasonable degradation with respect to the linear setup.

Conclusion
In this paper, we address the problem of sparse representations in a non-linear setting. While a high computational cost of the gradient of the goal function may prevent the use of standard estimation procedures, we show that it can be overcome by applying principles from the theory of optimal control, as long as the cost function satisfies some desirable structural property. Our derivations are illustrated on a particular example dealing with the estimation of the state of a geophysical system from partial observations. Abstract-This short communication gives an overview of the main contributions from our work [6]. We put here the emphasis on the interplay between discrete (i.e. over a finite grid) and continuous (i.e. grid-free) sparse regularizations of inverse problems using 1 -type regularizations. We list results ensuring that the support of the recovered sum of Diracs matches either exactly or approximately the one of the input measure when the signal-tonoise level is large enough. We also study the convergence of the discrete problem toward the continuous one when the grid size tends to zero.

Introduction
We focus our attention to a class of super resolution problems, where one aims at recovering a sparse sum of Diracs from low-resolution noisy observations. This typical setup is an idealization of many problems in imaging sciences, from seismic sparse spikes deconvolution to source localization in medical imaging.
Discrete 1 regularization. We focus our attention in this article to 1 -type techniques. They were initially proposed in geophysics and later studied in depth by David Donoho and coworkers [4]. A popular approach to analyze the performance of discrete 1 regularization is to assess the recovery of the positions of the non-zero coefficients, as initially introduced by Fuchs [8] and exploited by Dossal and Mallat in [5] to study sparse spike deconvolution.

Continuous total variation of measures regularization.
Working over a discrete grid makes the mathematical analysis difficult. Following recent proposals (see for instance [1,3]), we consider here sparse regularization (TV) over a continuous domain, i.e. in a grid-free setting. In this setting, the discrete 1 norm is replace by the so-called total variation of a measure (not to be confounded by the total variation of a function, which is popular to recover piecewise-constant signals and images) This allows a sharper analysis of the performances because spikes are not constrained to be a on grid, and can thus slightly shift spatially.
For deconvolution from ideal low-pass measurements, the ground-breaking paper [3] shows that this TV regularization recovers exactly sum of Diracs that are well separated. This work is further refined in [2] that studies the robustness to noise. In view of the applications of superresolution, it is crucial to understand the precise location of the recovered Diracs locations when the measurements are noisy. Partial answers to this questions are given in [7] and [1], where it is shown (under different conditions on the signal-to-noise level) that the recovered spikes are clustered tightly around the initial measure's Diracs. In this article, we fully answer the question of the position of the recovered Diracs in the setting where the signal-to-noise ratio is large enough.
Contributions We first review in Section 3 a classical result which ensures perfect support recovery over a discrete grid. We explain in Section 4 that the natural generalization of this approach in the grid free setting necessitates using a more elaborate criterion, which leads to a finer performance analysis. This is because the continuous setting allows spikes to slightly move from there initial positions, which is forbidden over a fixed grid. Lastly, Section 5 gives a precise description of the discrete-tocontinuous limit, which allows a better understanding of the solution of the discrete problem when the noise is small. Note that the code to reproduce the figures of this article is available online 1 . Students and teachers might also find useful the companion numerical tour on sparse regularization over measures 2 .

Sparse Inverse Problems Regularization
We consider the problem of recovering an input Radon measure m 0 ∈ M(T), where T = R/Z, from noisy low resolution observations y = Φ(m 0 ) + w ∈ L 2 (T) where w ∈ L 2 (T) is some noise and Φ(m)(x) = T ϕ(x, x )dm(x ) is a linear operator with a smooth kernel ϕ ∈ C 2 (T × T). We note Φ * (y)(x) = T y(x )ϕ(x , x)dx the adjoint defined on L 2 (T). Note that for simplicity we focus our attention to 1-D inverse problems, but our results extend to arbitrary dimension.

Discrete Regularization
A classical way to regularize the inverse problem is to use 1 regularization over a discrete grid z = (z i ) N i=1 ⊂ T N of N points, and solve the finite dimensional inverse problem where λ > 0 should be adapted to the noise level, where Φ z a = i ϕ(·, z i )a i is the restriction of the operator to measures supported on the grid z. For a solution a of (1), one thus obtains an estimate m z,a = N i=1 a i δ zi of the initial measure m 0 .
We assume that the input measure m 0 = m z,a0 is supported on the grid, with a support I = Supp(a 0 ) of non-zero indexes. The following classical theorem shows that if the following precertificate (we assume Φ I injective) η F = Φ * Φ * ,+ I sign(a 0,I ) where Φ I = (ϕ(·, x i )) i∈I . (2) does not saturate on the grid outside the support, then the discrete regularization recovers the correct support. Here a 0,I = (a 0,i ) i∈I is the sub-vector indexed by I and Φ * ,+ I denotes the pseudo-inverse of Φ * I .

Continuous Regularization
A major issue is that Theorem 1 is overly pessimistic in practical applications, because recovering the exact positions of the spikes is not possible if the grid size is too small. One can indeed show that the condition |η F (x j )| < 1 is always violated for a grid fine enough. To allow for a finer analysis, one can replace (1) by its counterpart over the space of measures min m∈M(T) where the supremum is taken over continuous functions. Note that for a discrete measure, one has ||m z,a || TV = ||a|| 1 .
We now assume an arbitrary input measure m 0 = m x0,a0 where x 0 ∈ T s . While the discrete pre-certificate η F defined in (2) is only constrained to interpolate the sign of a 0 on the support, working over a grid-free domain enables to use refined constraints. The following Theorem 2 shows that support recovery is indeed governed by a pre-certificate which interpolates the sign of the input measure together with vanishing derivatives at the support locations where we have introduced the following operator (assumed to be injective) Theorem 2 (Duval-Peyré [6]). If ∀ i, η V (x 0,i ) = 0 and ∀ t / ∈ x 0 , |η(t)| < 1, then for ||w|| small enough and λ ∼ ||w||, the solution m λ = m xλ,aλ of (3) is unique, is composed of s Diracs, and ||a λ − a 0 || ∞ = O(||w||), ||x λ − x 0 || ∞ = O(||w||).
The important point is that (contrary to η F ) since η V has vanishing derivatives on the support, the non-saturation constraint (i.e. |η| < 1 outside the support) is much more likely to hold in the continuous setting (i.e. η = η V ) than in the discrete one (i.e. η = η F ). This is illustrated by Figure 1, on the case where ϕ(x, x ) = sin((2f c + 1)π(x − x ))/ sin(π(x − x )) is an ideal low-pass filter with cut-off frequency f c = 6.

From Discrete to Continuous
The following results shows that if the continuous problem has a stable support (i.e. Theorem 2 can be applied) then the solution of the discrete problem is supported on pairs of Diracs surrounding the initial measure. For simplicity, we suppose that the considered grids z are uniform.
Theorem 3 (Duval-Peyré [6]). Let m 0 = m x0,a0 be such that x 0 ⊂ z is on the discretization grid, i.e. x 0 = z I for some index set I. We suppose that the hypotheses of Theorem 2 hold for m 0 . Then if N is large enough, ||w|| small enough and λ ∼ ||w||, the solution of a λ of (1) is unique and Supp(a λ ) ⊂ I ∪Ĩ wherẽ This result is illustrated on Figure 2, in the case of an ideal low pass filter with no noise. It shows that for λ = 0, the discrete problem recovers exactly the input measure, but as λ increases (and similarly if some small noise w is added), two spurious spikes, located on neighboring grid points indexed byĨ, pollutes the recovered measures. This duplication of spikes is linked to the inability of the discretized representation of measure to cope with translation of the Diracs, as this is achieved by the continuous one studied in Section 4. At the same time, this gives a precise geometrical understanding of the convergence of the discrete problem toward the continuous one (duplicated spikes become closer and closer to the input spikes as N increases). Abstract-We consider in this paper an efficient classification scheme built by applying a soft-thresholding non linear map followed by a linear classifier. We propose a supervised dictionary learning algorithm tailored for this classification architecture. We cast the dictionary learning problem, that learns jointly the discriminative dictionary and the hyperplane, as a difference of convex (DC) problem and solve it efficiently with DCA. We show on texture classification experiments that our method compares favorably to standard classifiers in terms of accuracy and computational complexity.

Introduction
Object classification is one of the most challenging tasks in computer vision. Recently, techniques relying on sparse coding for feature extraction have led to state of the art results on standard datasets [1,2,3,4]. However, the prohibitive cost of computing the non-linear sparse representation of an image is a major drawback that limits the applicability of such techniques to large-scale vision problems. Indeed, the sparse representation in a dictionary D of an image x involves a non trivial optimization problem. The development of efficient techniques for solving (or approximating the solution of) Eq. (1) has been the focus of much recent research [5,6,7,8]. In this work, we consider an efficient and simple alternative to the sparse coding mapping of Eq. (1). Specifically, we consider the soft-thresholding 1 non-linear feature extraction mapping defined by where (·) + = max(0, ·) and α is a parameter that controls the sparsity of the feature vectors. The usefulness of this nonlinear mapping for classification has already been observed [3] and is popular in a number of deep learning architectures due to its efficiency [9,10]. We focus here on a classification scheme where features (obtained with the soft thresholding mapping) are fed to a linear classifier w. We design a supervised dictionary learning algorithm tailored for this classification architecture. Specifically, we learn jointly a discriminative dictionary D and a hyperplane w that separates the two classes in the feature space. We illustrate the effectiveness of our approach in terms of classification accuracy and computational efficiency on a texture classification task.

Problem formulation
We focus in this paper on a binary classification problem. Let Soft-thresholding is usually defined as sgn(x)(|x| − α). We consider only the positive part here. denote the set of training points and their associated labels. We consider the following supervised dictionary learning formulation, that optimizes jointly for the dictionary D ∈ R n×N and the linear classifier where L denotes a convex loss function that penalizes incorrect classification of a training sample and ν is a regularization parameter which prevents overfitting. The above optimization problem attempts to find a dictionary D and a linear separator w such that w T (D T x i − α) + has the same sign as y i on the training set. We adopt the hinge loss function, in this paper, which is defined by L(x) = max(0, 1 − x). The problem in Eq. (3) is difficult to solve in general. We use a suitable change of variables to put it in a nice form (difference of convex functions) that provides a practical and efficient solution. Specifically, using the change of variables u j ← |w j |d j , v j ← |w j | and s j ← sgn(w j ), the problem in Eq.(3) can be rewritten 2 in the following way: subject to v > 0.
Note that s j (i.e., the sign of w j ) essentially encodes the "class" of atom d j . In other words, an atom d j for which s j = +1 is most likely to be active in class 1 samples. Conversely, atoms that satisfy s j = −1 are most likely active for class −1 samples. We assume here that the vector s is known a priori. In other words, this means that we have a prior knowledge on the proportion of class 1 to class −1 atoms in the dictionary. For example, setting half of the entries of the vector s to be equal to +1 and the other half to −1 encodes the prior knowledge that we are searching for a dictionary with balanced class-specific atoms. Moreover, for the purpose of optimization, we approximate the term (u T j x i − αv j ) + in Eq.(4) with a smooth function q(u T j x i − αv j ) where q(x) = 1 β log (1 + exp (βx)) and β is a parameter that controls the accuracy of the approximation. Note that this approximation is done only to make the optimization easier at the dictionary learning stage. At test time, soft-thresholding is applied for feature extraction.
Thus, we end up with the optimization problem (P) 3 : Algorithm 1 DCA algorithm for solving (P) 1. Choose any initial point: 2.2 Solve the convex optimization problem: Given the optimal (U, v), D and w can be obtained by using the previously given change of variables. A test point x is then assigned to label +1 if w T (D T x − α) + > 0, and label −1 otherwise.

DC optimization
Problem (P) is typically a nonconvex optimization problem that is hard to solve using traditional methods, such as gradient descent or Newton-type methods. We show however that the problem (P) can be written as a difference of convex (DC) program [11]. A DC problem is an optimization problem which consists in minimizing a function of the form g − h, where both g and h are convex, subject to difference of convex constraints. DC problems are well studied optimization problems and efficient optimization algorithms have been proposed [11,12] and have been shown to work well in practice (see [13] and references therein). Proposition 1. For any convex loss function L and any convex function q, the problem (P ) is DC. Moreover, when L(x) = max(0, 1 − x), the objective function of problem (P) is equal to g − h with: We solve efficiently problem (P) for a local optimal solution using the DC minimization algorithm (DCA) [12] presented in Algorithm 1. In a nutshell, each iteration of DCA solves a convex optimization problem obtained by linearizing h around its current point. Although DCA is only guaranteed to reach a local minima, the authors of [12] state that DCA often converges to a global optimum.

Experimental results
We now illustrate the effectiveness of our approach on texture image classification experiments. We consider texture images from the 32 Brodatz dataset 4 . For each pair of textures under test, we build the training set by randomly selecting 200 patches per texture of size 12 × 12 from the left half of the images, and the test data is constructed by selecting 1000 patches per image from the right half of the images. For our discriminative dictionary learning approach, we fix half of the entries of s to be equal to +1 and the other half to be equal to −1. We set moreover α = 1, β = 50 and the regularization parameter ν to 0. We initialize half of the vectors in U with random patches from the training set of class 1, and the other half with random patches from the training set of class 2. We also set v 0 to 1. At test time, a datapoint x is assigned to class 1 if w T (D T x − α) + > 0 and class −1 otherwise.
We compare our approach to several classification algorithms. The first of these methods is a simple linear SVM classifier. Also, we consider a sparse representation based classifier similar to that of [2], where a dictionary is first learned (in a reconstructive manner) that sparsifies the training data 5 . A linear SVM is then trained on the sparse codes. At test time, the classification is performed by solving Eq. (1) on the test image, and then applying a linear classifier. Finally, the last comparative scheme is a nearest neighbour classifier. The dictionary is built by applying a K-Means clustering on the training samples of each class separately. At test time, the sample is assigned the label of the dictionary atom (i.e., cluster) that is closest to it. Table 1 compares these methods on 5 binary classification tasks, with different dictionary sizes. Our proposed approach clearly outperforms linear and nearest neigbour classifiers in all tasks. Moreover, our method compares favorably to sparse coding for small dictionary sizes. For larger dictionaries, sparse coding yields in some cases a better classification accuracy. This is due to the fact that the dictionary learning step in classification tasks is crucial when the dictionary is constrained to be small [3]. For larger dictionaries, there seems to be less performance gap between reconstructive and discriminative approaches for learning D. As shown in Table 2, softthresholding based features are much more efficient to compute than sparse codes, which makes our method particularly suited to large-scale problems or when computational power is limited. In a longer version of this work (in preparation), we provide more experimental results and compare the proposed approach to other algorithms.

Introduction
Nonnegative matrix factorization (NMF) is a powerful dimensionality reduction technique as it automatically extracts sparse and meaningful features from a set of nonnegative data vectors. Given n nonnegative m-dimensional vectors gathered in a nonnegative matrix M ∈ R m×n + and a factorization rank r, NMF computes two nonnegative matrices W ∈ R m×r + and H ∈ R r×n + such that M ≈ W H. Unfortunately, NMF is NP-hard in general [1]. However, if the input data matrix M is r-separable, that is, if there exists an index set K with cardinality r and an r-by-n nonnegative matrix H such that M = M (:, K)H, then the problem can be solved in polynomial time [2]. In the presence of noise, this problem is referred to as near-separable NMF. Although the separability condition is rather strong, it makes sense in several applications, e.g., in text mining [3], hyperspectral unmixing [4] and blind source separation [5]. Note that this problem is equivalent to finding a rowsparse matrix X such that ||M − M X|| 2 F is minimized [6,7]. In this talk, we present the results from [8]: we show how to make near-separable NMF algorithms more robust using preconditioning. We focus on a simple yet effective near-separable NMF algorithm, namely, the successive projection algorithm.

Successive Projection Algorithm
The successive projection algorithm is a simple but fast and robust recursive algorithm for solving near-separable NMF; see Algorithm SPA. At each step of the algorithm, the column of the input matrixM with maximum ℓ 2 norm is selected, and thenM is updated by projecting each column onto the orthogonal complement of the columns selected so far. It was first introduced in [9], and later proved to be robust in [10].

Preconditioning
Assume for simplicity that m = r. If we were given the full rank r-by-r matrix W , we could premultiply the input nearseparable matrixM = W [I r , H ′ ] + N with W −1 and obtaiñ The matrixM ′ is also near-separable, and is now perfectly conditioned, although the noise ǫ ′ = max i ||W −1 N :i || 2 might have increased but by a factor of at most σ −1 min (W ). Under the same conditions as in Theorem 2, one can check that ǫ ≤ O σmin(W ) √ r would imply that SPA identifies the columns of W up to error O (ǫ κ(W )). This is a significant improvement compared to Theorem 2, both for the noise level ǫ (improvement of a factor κ −2 (W )) and the error (improvement of a factor κ(W )).
Of course, the matrix W is unknown, otherwise the problem would be solved. However, it can be shown that a minimum volume ellipsoid problem allows to approximate W −1 (up to orthogonal transformations): Given a matrixM ∈ R r×n of iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. rank r, we can formulate the minimum volume ellipsoid centered at the origin and containing the columnsm i 1 ≤ i ≤ n of matrixM as follows (1) This problem is SDP representable [11, p.222]. In [8], it is shown that, ifM satisfies Ass. 1, then A * ≈ (W W T ) −1 , hence factoring A * (e.g., using the Cholesky decomposition) allows to recover W −1 approximately (up to orthogonal transformations). In fact, in the noiseless case (that is, for N = 0), A * = (W W T ) −1 . Finally, our new proposed algorithm, referred to as preconditioned SPA (prec-SPA), works as follows: (i) Solve Problem (1) in order to approximately compute P ≈ W −1 , (ii) PremultiplyM with the preconditioning to ob-tainM ′ = PM , and (iii) Apply SPA on the pre-conditioned matrixM ′ . We have the following robustness result:

Numerical Experiments
In this section, we illustrate the robustness of prec-SPA compared to the original SPA on some synthetic data sets. We also compare it to vertex component analysis (VCA) [12], a popular endmember extraction algorithm, and XRAY from [13]. The Matlab code is available at https://sites.google. com/site/nicolasgillis/. Please see [8] for more details and other numerical experiments. The near-separable matrixM is generated as follows: We take m = r = 20, and n = 210. The matrix W is generated using the rand(.) function of Matlab, that is, each entry is drawn uniformly at random in the interval [0, 1]. The matrix H = [I r , H ′ ] is such that H ′ has exactly two non-zero entries in each row equal to 0.5. Hence, the 190 data points are in the middle of two different columns of W . The noise is chosen such that the columns of W (that is, the first 20 columns of M ) are not perturbed, while the 190 data points are moved towards the outside of the convex hull of the columns of W : where M = W H andw = 1 r i w i is the vertex centroid of the convex hull of the columns of W . These are the same near-separable matrices as in [10]. For different noise levels δ, we generate ten such matrices and Figure 1 reports the percentage of columns of W correctly identified by the different algorithms. Clearly, prec-SPA outperforms SPA, VCA and XRAY.

Conclusion
We have presented a way to precondition near-separable NMF matrices using semidefinite programming. This is turn allowed us to robustify near-separable NMF algorithms. In particular, the preconditioning makes the popular successive projection algorithm (SPA) provably more robust to noise, which we have illustrated on some synthetic data sets.
The suboptimal nature of OMP and OLS has led many researchers to study conditions under which these procedures succeed in recovering the true sparse vector. This question has been widely addressed for OMP in the recent years, including worst-case uniform [7,8] and probabilistic analyses [9]. Although OLS has been known in the literature for a few decades (under different names [10]), exact recovery analysis for OLS have only appeared very recently, see [11,12,13].
Most of the existing works deal with uniform guarantees: these conditions ensure the success of OMP/OLS for a given sparsity level (or a given support) irrespective of the magnitude of the non-zero coefficients. In contrast with these works, we derive new guarantees of success accounting for the decay of the non-zero elements of the sparse vector. Our conditions are expressed in terms of the mutual coherence of the dictionary µ and encompass, as particular cases, some well-known results of the literature. The proofs of the results are reported in our technical report [14].

Context and Main Results
Let us assume that y ∈ R m is a (noisy) linear combination of k columns of A ∈ R m×n indexed by Q , that is where w ∈ R m denotes some noise vector and Card{·} stands for the cardinality operator. We assume that the columns a i of the dictionary are normalized: a i 2 = 1 ∀i. Let us first consider the noiseless case (w = 0). Theorem 1 provides sufficient conditions of success for OMP/OLS accounting for the possible decay of non-zero coefficients in x. In our statement, we assume without loss of generality that and and then Oxx selects atoms in Q from noiseless data during the first k iterations.
In our technical report [14], we show that this kind of result can also be extended to the characterization of the success of OMP/OLS in the noisy setting ( w 2 ≤ ) or the partial recovery of the support by the latter algorithms. In particular, we show that a well-known result by Donoho et al. [15,Th. 5.1] can be further relaxed as and then Oxx selects atoms in Q during the first k iterations.
More specifically, we see that the condition by Donoho et al. in [15,Th. 5.1] is sufficient for (7) to be satisfied and is therefore stronger than the conditions mentioned in Theorem 2. On the other hand, the two types of conditions become equivalent as soon as the nonzero elements of x have the same magnitude.

INTRODUCTION
Audio signals are prone to various degradations and clipping is among the most common ones. It may arise during recording due to the poor dynamic range of a microphone, or of analog-to-digital circuits, but it may be seen also as an extreme case of magnitude compression [1].
Clipping inevitably causes an information loss: in the time domain, the signal is "cut-off" when its amplitude is above some threshold. This manifests as an expansion (spread) in the frequency domain. An intuitive approach to de-clip clipped signals may be to "clean" their frequency spectrum and preserve only those frequencies needed to represent the original signal accurately. However, estimation of the correct support in the frequency domain is also a difficult problem. Hence, we need some additional prior information about the underlying signal.
We formalize clipping as the following element-wise operation on the original signal x ∈ R n (yielding the clipped signalx):x = sign(x) min(|x|, τ ) The threshold τ represents amplitude saturation level (here we assume symmetric clipping, i.e. the same threshold is applied to both positive and negative samples of x). Let M r (x) be the sampling operator which extracts unclipped, "reliable", samples, and M + c (x), M − c (x) be the operators which extract clipped samples with positive and negative amplitude, respectively. Knowing that (for the clipped positive/negative samples) the amplitude of the original signal has to be above (respectively, below) the threshold, we can formulate the inverse problem of signal de-clipping in the following way: This is an ill-posed inverse problem, since there are infinitely many possible solutions forx. We chose to regularize the problem by introducing a sparsity hypothesis.

SPARSE ANALYSIS DATA MODEL
The premise of sparse regularization is that many signals have low-dimensional representation, but not necessarily in their original domain. This assumption is sometimes sufficient to mitigate ill-posedness of the initial problem.
The sparse analysis, also known as the cosparse data model assumes that the signal in question can be "sparsified" by applying an adequate linear transform. Let Ω ∈ R p×n , (p ≥ n) denote the matrix form of this sparsifying transform and let the linear constraints in (2) be denoted by Γ. We introduce the following optimization problem: However, minimizing this 0 -"norm" is NP-hard, but there are convex relaxation and greedy algorithms [2] which may be used to approximate its solution. The cosparse data model has only recently attracted the attention of scientific community, and to our best knowledge, it has not yet been exploited to solve the declipping problem.

COSPARSE DECLIPPING
Our approach is based on transposing the Consistent Iterative Hard Thresholding (Consistent IHT) algorithm [3] to the cosparse data model. Consistent IHT relies on applying hard iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology. (38) thresholding operator H K (s) 1 to the transform coefficient vector s in order to approximate sparse solutions. The transform coefficients are estimated by performing a gradient descent step of an objective function that incorporates the constraints expressed in Γ.
A pragmatical issue is that computing the k-cosparse projection is proven to be NP-hard [4]. To circumvent this problem, we take a different approach, largely based on the Alternating Direction Method of Multipliers [5] framework. The resulting algorithm is presented as pseudocode in algorithm 1 and will be referred to as Cosparse Declipping by Hard Thresholding (CoDec-HT).

8:
k ← k + 1 9: go to 2 10: end if 11: returnx (k) Mimicking the heuristics used in Consistent IHT, we relax the sparsity k of z through iterations, which allows to "learn" the unknown sparsity level of a given signal. Through iterations, the vectors Ωx andẑ eventually get close to each other, thus the signal estimatex becomes approximately cosparse.
The analysis operator Ω can be chosen as any (possibly overcomplete) transform which sparsifies the audio signal. In our experiments we used concatenated Discrete Cosine Transform and Discrete Sine Transform matrices.

EXPERIMENTS
We carried declipping experiments with 10-second excerpts from two wideband music pieces, sampled at 16kHz with 16bit encoding. Audio 1 is a rock piece with vocals, drums, two electric guitars and bass, while Audio 2 is an instrumental piece performed on an analogue synthesizer. The clipped versions of these signals are divided into overlapping frames of length 1024 before processing. Afterwards, the outputs are re-synthesized by the overlap and add scheme.
Beforehand, the audio is clipped given a predefined SNR inp (signal-to-noise ratio of the input) value. The restoration performance is measured as the output signal-to-noise ratio (SNR out ): SNR inp = 20 log 10 x 2 x −x 2 , SNR out = 20 log 10 x 2 x −x 2 (4) The proposed algorithm is compared against [3] and [6] 2 , which serve as a reference. Consistent IHT (with dictionary Ψ = Ω H ) and the proposed algorithm are based on a similar algorithmic methodology and differ only by the underlying data model. Algorithm [6] exploits a structured prior termed social sparsity [7].
The results 3 presented in Figure 1 suggest that cosparsitybased hard thresholding outperforms synthesis-based one at all clipping levels. Interestingly, the SNR improvement rate of CoDec HT is roughly of the same order as for the social sparsity based algorithm, despite the fact that the latter uses a stronger prior.

CONCLUSION
Experimental results indicate that sparse analysis regularization can be successfully applied to the audio declipping inverse problem. Future work will aim at incorporating structured cosparsity in the data model to improve declipping performance, and at addressing extended magnitude corruption scenarios, such as dynamic range compression or soft clipping. Abstract-Dictionary learning is a branch of signal processing and machine learning that aims at expressing some given training data matrix as the multiplication of two factors: one dense matrix called dictionary and one sparse matrix being the representation of the data in the dictionary. The sparser the representation, the better the dictionary. However, manipulating the dictionary as a dense matrix can be computationally costly both in the learning process and later in the usage of this dictionary, thus limiting dictionary learning to relatively small-scale problems. In this paper we consider a general structure of dictionary allowing faster manipulation, and give an algorithm to learn such dictionaries over training data, as well as preliminary results showing the interest of our approach.

Introduction
The use of a dictionary to sparsely represent a certain type of data goes back to almost two centuries with the Fourier transform that was designed to sparsely represent heat flow at that time [6]. The Fourier and the Hadamard transforms, as well as the wavelets to cite just a few rely on a "simple" mathematical formula that has been shown to yield fast algorithms [5,9].
On the other hand, more recently have arisen algorithms that learn a dictionary directly over training data, without using an analytical formula (see [11] and references therein for a survey on the topic). They consider some data matrix X ∈ R d×n , which is the collection of n training vectors x i ∈ R d , and can be approximated sparsely in a certain dictionary: Γ having sparse columns. Such algorithms provide dictionaries that are well adapted to the data but they do not lead to fast algorithms because of the lack of structure of the dictionary, and storing D as a dense matrix may be impractical. Bridging the gap between analytical dictionaries that are not necessarily well adapted to the data at hand but lead to fast algorithms, and learnt dictionaries that are very well adapted to the data but do not lead to fast algorithms is an important challenge. The dictionary learning community started recently to tackle the problem, mainly in [12] and [4], where the authors introduced new dictionary structures. The model we introduce below generalizes these works, as we will see in Sec.2.
In this paper we build on the very simple observation that the fast algorithms used to apply analytical transforms like the FFT or the DWT can be seen as consecutive multiplications of the input vector by sparse matrices. This fact implies that the analytical dictionaries associated with such fast transforms can be expressed as a product of sparse matrices, that is 1 : This structure is precisely what makes possible the appearance of fast algorithms (such as the butterfly FFT for example). 1 The product being taken from left to right: Knowing that, it makes sense to learn a dictionary that has this form and is the product of several sparse matrices, making it intrinsically fast and easier to store. The problem at hand is unfortunately highly non-convex and the sparsity constraint is non-smooth, but recent advances in optimization like the algorithm proposed in [2] can be adapted to such problems with convergence guarantees as we will show in the next sections.

Model and optimization problem
Notation Throughout this paper, matrices are denoted by bold uppercase letters: A. Vectors are denoted by bold lowercase letters: a. The ith column of a matrix A is denoted by: a i . Sets are denoted by calligraphical symbols: A.
As stated in the introduction, we propose a new dictionary structure that intrinsically leads to fast algorithms although the dictionary is learnt over some training data. Our goal is to find a dictionary that sparsely represent the data, being itself a multiplication of sparse matrices.

The matrix factorization point of view
Let D ∈ R d×a be our dictionary with a atoms and Γ ∈ R a×n the corresponding sparse representation matrix such that X ≈ DΓ. In order to meet the requirements and be intrinsically fast, D must take the form of equation (2), where the S j s are sparse matrices in R aj ×aj+1 with a 1 = d and a M +1 = a. Now if we say that Γ is now called S M +1 for ease of notation, our goal is to find the sparse factors S j s such that: We see with this equation that our structured dictionary learning problem amounts to a factorization of the data matrix into M + 1 sparse factors. Actually, this model is quite general and can encompass those introduced in [12] and [4] as we will see. Such a multifactor representation of a data matrix has been introduced in the NMF framework in [8], for deep learning in [1,10] and statistics in [7,3].

Optimization objective
The proposed structured dictionary learning problem can be expressed as a constrained optimization problem. In its most general form, it can be stated as follows: where the U j s are the sets in which each factor should lie.
In [12], the authors propose to constrain each atom of the dictionary to be a sparse linear combination of the atoms of a iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
(41) so-called base dictionary. This base dictionary is assumed to be asociated with a fast algorithm (it takes the form of equation (2)), thus leading to fast manipulation of the whole learnt dictionary. It is actually equivalent to solving the problem (4) with the M − 1 leftmost factors being known (the base dictionary), S M and S M +1 = Γ to estimate, and U M and U M +1 denoting sparsity constraints. They give an algorithm to estimate jointly the dictionary and the representation. However, such a structure constrains the dictionary to be close to the base dictionary and thus does not provide full adaptability.
In [4], the authors constrain each atom in the dictionary to be the composition of several circular convolutions with sparse kernels, thus leading again to fast manipulation of the learnt dictionary. They give an algorithm to estimate the dictionary knowing the representation and the support of each kernel. It is actually equivalent to solving the problem (4) with S M +1 being known, U i , ∀i ∈ {1 · · · M } denoting the intersection of sparsity constraints with the set of circulant matrices, and the support of each factor being known.
We propose to handle problem (4) in a more general form, namely we will not assume that any factor is known, and we will consider any U i s provided that we are able to perform the projection onto it (note that in practice, only the U i s that are included in the set of sparse matrices will be of interest).

Algorithm
We propose to apply an algorithm stemming from recent advances in non-convex and non-smooth optimization to problem (4). The algorithm is introduced in [2] and called PALM (Proximal Alternating Linearized Minimization). PALM is made to handle objective functions that take the form: where each x i ∈ R ni is a block of variables, H is only assumed to be C 1 and the f i s are only assumed to be proper and lower-semicontinuous (possibly indicator functions of constraint sets). The idea behind PALM is to alternate updates between different blocks of variables (Gauss-Seidel), performing a proximal gradient step for each block. PALM is guaranteed to converge to a critical point of the objective function. Problem (4) is obviously non-convex and thus the algorithm introduced in this paper is not guaranteed to converge toward the global minimum. Adapting PALM to the structured dictionary learning problem we just introduced is quite straightforward. One option is to take each S j as a block of variables, so that H is the objective function of problem (4) and the f i s are the indicator functions of the sets U j . The proximal operator of the indicator function of a set X reducing to the projection operator onto this set P X (.), the adaptation of PALM to our problem is given in Algorithm 1, where c i j is a stepsize depending on the Lipschitz modulus of the gradient ∇ Sj H (see [2] for more details).

Algorithm 1 PALM for structured dictionary learning
for i ∈ {1 · · · N iter} do for j ∈ {1 · · · M + 1} do end for end for 3 One practical example Choosing the sets of constraints U j is crucial in order to avoid as much as possible local minima, but in order to have a fast dictionary, they have to be subsets of the following (set of k jsparse matrices in R aj ×aj+1 ): In further work, we intend to set constraints such as the combination of sparsity and orthogonality, or some sort of structured sparsity, because we observed that the factorizations of the most famous analytical dictionaries are of this type. However in this paper we will consider the simplest constraints, namely we just want the matrices S j s to have at most k j nonzero entries. This means that we will have U j = S aj ×aj+1 kj . In this simple configuration, the projection operator P Uj (.) is simply the hard thresholding operator that sets all but the k j greatest entries (in absolute value) to zero.
Complexity savings In the setting we just described, the multiplication of a vector by the dictionary or its adjoint (and its storage) would cost O( M j=1 k j ) operations instead of O(da) operations for a dictionary without structure. The interest here is twofold: first at the learning stage since this algorithm requires multiplications by the dictionary at each iteration, second at the usage stage since it involves again multiplications by the dictionary and its adjoint. It is thus necessary to tune M and the k j s in order to have a significant complexity gain, namely M j=1 k j da. For example in the FFT case with d = a = 2 P , P ∈ N, the dictionary is factorized into M = log 2 d factors that are all 2d-sparse, thus we perform multiplication in O(2d. log 2 d) operations instead of O(d 2 ).
Preliminary result We are still testing different settings for the algorithm, but it shows promising results in tasks such as denoising. To illustrate this, we present an example of result with d = 8, a = 20, n = 100, M = 3, k 4 = 100, k 1 = k 2 = k 3 = 20 and N iter = 1000. We generated randomly some data matrix X 0 following the model of equation (3) and we added white Gaussian noise B with a signal-to-noise ratio of 6dB to obtain X = X 0 +B. We ran Algorithm 1 in order to obtainX 0 . We repeated the experiment 100 times and got an output signalto-noise ratio of 7.7dB in average. We also made the interesting though expected observation that in approximation task (when we just want to factorize some matrix in sparse factors), there is a trade-off between sparsity of the factors and expressiveness of their product: the more non-zero entries we allow, the better the performances are, but the complexity is also increased.

Conclusion
In this abstract, we first presented the dictionary learning problem and the growing concern in the community to make it more computationally efficient. We then introduced a new model that is a generalization of two previously existing ones and that leads to intrinsically fast dictionaries. We presented an algorithm with convergence guarantees to a stationary point that is able to learn such type of dictionary over some training data. We finished by showing briefly some preliminary results of our approach. In the future, we intend to set new configurations for this method. We could for example modify the constraints in order to avoid as much as possible the numerous local minima inherent to the problem.

INTRODUCTION
Many convex optimization problems in image and signal processing can be solved using a relaxed and inexact iterative scheme of the form H is a nonexpansive operator on a Hilbert space H, Øλ k Ù kÈN È×0, 1×, and for x k È H, ε k È H is the error when computing T x k . One simple instance of such an algorithm is that of relaxed inexact gradient-descent for solving where f is proper convex and has β ¡1 -Lipschitz continuous gradient, in which case T Id ¡γ∇f, for γ È×0, 2βÖ. In this case, the error ε k is that when evaluating ∇f at x k .
Many structured convex optimization problems boil down to implementing the same type of iterations. Consider for instance This work has been supported by the ERC project SIGMA-Vision and l'Institut Universitaire de France.
where f È Γ 0 ÔHÕ has β ¡1 -Lipschitz continuous gradient, h i È Γ 0 ÔHÕ, and Γ 0 ÔHÕ is the class of lower semicontinuous, proper, convex functions from a Hilbert space H to × ¡ , ×. We assume that some technical domain qualification conditions are verified. Problem (1.2) was considered in [1].
In the last decades, based on the notion of the proximity operator [2], and assuming that the functions h i are simple (i.e.their proximity operators are easily computable), a wide range of proximal splitting algorithms have been proposed to solve (1.2). One can cite for instance the forward-backward splitting method (FBS) valid for n 1, the Douglas-Rachford splitting method (DRS) that applies for f 0, a generalized Forward-Backward splitting method (GFB), or primal-dual splitting methods. See [3] for a comprehensive account.
In this paper, we first establish the iteration complexity of the inexact relaxed fixed point iteration (1.1). We then build upon this result to provide iteration complexity bounds for several proximal splitting algorithms, from which an easily verifiable termination criterion for finding an approximate solution will be provided.

ITERATION COMPLEXITY OF INEXACT NONEXPANSIVE OPERATOR
Define e k ÔId ¡ T Õx k .  (i) Øe k Ù kÈN converges strongly to 0, Øx k Ù kÈN converges weakly to a point x AE È fixT ; Next we present the ergodic iteration complexity bound of (1.1), define Λ k k j 0 λ j , andē k 1 Λ k k j 0 λ j e j . Theorem 2.2 (Ergodic iteration complexity bound of (1.1)). Assume that condition Theorem 2.1 (c) holds and let C 3 jÈN λ j ε j , then, A special class of nonexpansive operators is that of αaveraged operators [3] for α È×0, 1Ö. The above two complexity bounds obviously apply, where condition Theorem 2.1 (b) can be changed to
In plain words, this means that OÔ1ßǫÕ iterations are needed to find a pair Ôx, g È h 1 ÔxÕÕ with the termination criterion g ∇f ÔxÕ 2 ǫ. This iteration-complexity improves to OÔ1ß ǫÕ in ergodic sense for the same termination criterion.
GFB Consider now n 1 in (1.2). [1] proposed a generalization in a product space H n of the FBS scheme to solve (1.2). Tn the GFB, T takes a more intricate form, omitted here for space limitation, and ε k absorbs the error when computing both the proximity operators prox γhi and the gradient. Let γ and λ k be chosen as for FBS. We prove that OÔ1ßǫÕ iterations are needed to find a pair ÔÔu i Õ 1 i n , gÕ with the termination criterion g ∇f Ô i ω i u i Õ 2 ǫ, where g È n i 1 h i Ôu i Õ, see [4] for details. Only OÔ1ß ǫÕ are however needed in ergodic sense.
Other splitting schemes Similar complexity bounds were established for DRS, ADMM and some some primal-dual splitting methods to solve more general monotone inclusion problems that encompass (1.2) as a special case; see [5] for details.

NUMERICAL EXPERIMENTS
As an illustrative example, let's consider the principal component pursuit (PCP) problem [6], and apply it to decompose a video sequence into its background and foreground components. The rationale behind this is that since the background is virtually the same in all frames, if the latter are stacked as columns of a matrix, it is likely to be low-rank (even of rank 1 for perfectly constant background). On the other hand, moving objects appear occasionally on each frame and occupy only a small fraction of it. Thus the corresponding component would be sparse. The results are depicted on Fig. 1.

Antoine Liutkus David Martina Sylvain Gigan Laurent Daudet
Institut Langevin, ESPCI ParisTech, Paris Diderot Univ., UPMC Univ. Paris 6 CNRS UMR 7587, Paris, France ABSTRACT Multiply scattering media in optics such as thin layers of ZnO -essentially white paint -behave as natural random linear multiplexers that can be exploited to yield opaque lenses or compressive imagers that beat their coherent counterparts both in terms of performance and ease of implementation. As such, they appear as a very promising research topic in many domains of optics. Still, current research is hindered by the fact that the linear model relating the input and the output of a multiply scattering material only holds for the complex optical field, which is difficult to measure without an interferometric experimental setup. Indeed, in a practical setup, only its amplitude is easily captured by a CCD camera. In this study, we show how such intensity measurements and digital micromirror devices can be exploited both for calibration, arbitrary focusing and compressive imaging using scattering materials.
Index Terms-intensity-based calibration, focusing, compressive sensing, optical imaging, scattering media

I. INTRODUCTION
Whereas optical scattering is usually seen as a time-varying nuisance, for instance when imaging through turbid media, the recent results of wave control in stable complex material have largely demonstrated that it could also be exploited, for example so as to build focusing systems that beat their coherent counterparts in terms of resolution [9]. In the context of compressive sensing for instance (CS, see [2]), such materials perform an efficient randomized multiplexing of the object into several sensors and were recently considered as analog randomizers [6]. Using such naturally occurring randomness instead of digitally-controlled pseudorandom sequences as in the single-pixel camera [3] permits both to drop the need for a designer to craft the randomizer himself, but also yields very compact sensing devices.
However, some practical difficulties come along using multiply scattering materials as multiplexing devices for focusing and imaging. Most importantly, the sampling mechanism is completely unknown a priori when using a multiply scattering material. Several studies have demonstrated that the input-output relationship of multiply scattering materials can very well be modelled through a linear Transmission Matrix (TM, see [8]), which can be used both for focusing [8], overdetermined imaging [7] and CS [6].
In this study, we go further previous work that only consider linear measurements and show how intensity-based calibration can be performed efficiently, which permits to drop the need for interferometric experimental setups and to simply use CCD cameras as sensors.
With the estimated TM, we show how focusing can be performed, where the controlling device is explicitly assumed to be a binary DMD. In comparison to the work presented in [8], [9], [1], the proposed approach has several advantages. First it does not require an interferometric mechanism nor an expensive Spatial Light Modulator (SLM) as in [8]. Second, it permits to build arbitrary output patterns to the system whereas [9], [1] only permit one single focusing point.

II-A. Theoretical background
Due to space constraints, we do not detail here the theory underlying the present study. The interested reader is refferred to [6] for more details. In short, the input to the system is a vector x 2 R N . The measurements are given as a vector y 2 R M , given by: where the M ⇥ N matrix H is called the transmission matrix (TM) and accounts for the action of the multiply scattering material. In practice, the TM H has some interesting properties. In particular, its entries can be assumed Gaussian independent and identically distributed (i.i.d.) random variables (see [6] and references therein).
With H estimated, an arbitrary output pattern y may be obtained by settling the input to x = H † y, where H † denotes the phase conjugate of H in a time-reversal interpretation, or any regularized inverse (see [6] for references). In practice, x cannot take any arbitrary values, but is often constrained to lie either on the complex unit circle (as in the case of a phase-only modulation using a SLM) or to be binary x 2 {0, 1} N as in the case of a binary DMD.

II-B. Calibration: measuring the measurement matrix
Assume that during a calibration stage, we can control the input x and measure the corresponding output y. Doing this L times and stacking the inputs and the outputs into matrices X and Y respectively, we have Y = HX. If we assume that these measurements are given with additive white Gaussian noise of variance 2 , we can estimate H throughĤ = Y X H XX H + 2 I 1 , where I is the N ⇥ N identity matrix. Correlated measurement noise can easily be taken into account, as well as a prior distribution over H in a Bayesian setting.

II-C. Experimental phase transitions under noise
All CS reconstruction algorithms exhibit some level of robustness to noise. However, when noise becomes prominent, performance of CS eventually drops. To overcome those limitations, it has been shown that it is advantageous to benefit from Multiple Measurement Vectors (MMV, [4]). In the MMV setting, P > 1 different measurements are available that correspond to the output of the system where the corresponding input all share the same support. In figure 1, we show observed performance of our compressive optical imager with linear measurements, using a varying number P of measurements per sparse object to reconstruct.

III. AMPLITUDE MEASUREMENTS III-A. Theoretical background
In the case of amplitude measurements as performed by a CCD camera for instance, the observed output is I = |y| = |Hx|. In that case, the TM H cannot be estimated simply be a least-squared error procedure as above because of the non-linearity introduced by the amplitude operator. However, it can be estimated in an iterative fashion by a dedicated algorithm that we cannot detail here due to space constraints.

III-B. Focusing with a binary DMD
Once H has been estimated, it can be used for both imaging and focusing. In the case of focusing, we want a particular output y and need the corresponding input x. When the input is to be binary as in the case of a DMD, we propose to use an antisparse coding scheme [5] that looks for x by solvinĝ Whereas the first term in 2 accounts for the focusing quality, the second term tends to minimize the maximal value of x,which in turns leads to x having only two values : ± kxk 1 , thus being binary.
In figure 2, we show an actual preliminary example of the focusing performance obtained with the proposed approach, where we want to focus on three target points at the same time. Intensitybased calibration was achieved with 4N measurements (N = 900 here). As can be seen, the proposed approach permits efficient focusing.

IV. CONCLUSION
In this paper, we have briefly summarized some important trends in optical imaging with multiply scattering materials. Those materials perform efficient multiplexing of input lights and are thus particularly efficient in the context of compressed sensing. We have shown how recent results and algorithms from the sparsity framework may be used to yield efficient optical devices that permit both focusing and imaging with a limited number of intensity calibration measurements.  Abstract-Data in high-dimensional spaces often aggregate near inherently low-dimensional structures. Manifold learning and dictionary learning motivate numerous techniques for exploiting such low intrinsic dimensionality. In manifold learning, a lowdimensional manifold constitutes the underlying structure of the data, and one "learns" this manifold from samples. In dictionary learning, the goal is to construct a dictionary (i.e. a set of vectors) such that every data point may be expressed as a sparse (compressible) linear combination of elements from this set. Manifold learning and dictionary learning are generally considered as separate problems, however, Geometric Multi-Resolution Analysis (GMRA) [2] provides a method for approximating (in a robust multiscale fashion) a low-dimensional set around which data may concentrate, while simultaneously producing a dictionary for representing the data sparsely.
In this work we introduce an estimator for low-dimensional sets supporting the data constructed from the GMRA approximations. We prove near-optimal finite sample bounds on its performance, and we also rigorously establish the robustness of this procedure with respect to both noise and model error. For data supported on a low-dimensional manifold, we demonstrate that this procedure admits error bounds primarily dependent upon the global parameters of the manifold: its intrinsic dimension and reach.

Introduction
In many high-dimensional data analysis problems, efficient representations dramatically boost the performance and efficiency of learning algorithms. In particular, it is often convenient to work with sparse representations of the data, meaning that there exists a dictionary (a large collection Φ = {φ i } N i=1 of fixed vectors) such that every observation can be represented as a linear combination of the elements of Φ with few nonzero coefficients (e.g., see [1,7,12]). In practice, the dictionary Φ is not known, and must be inferred algorithmically from the training data X = {X 1 , . . . , X n } consisting of n i.i.d. samples from some underlying distribution Π. In most cases, this approach results in a random error due to the fact that the dictionary is learned based on the empirical approximation of Π.
When the intrinsic dimension of the analyzed data is small, sparse representations generally become available. One of the simplest geometric models is the assumption that the R Dvalued observations lie close to a d-dimensional linear subspace Classical tools such as PCA [17,21] can estimate V and use its basis as a dictionary. One may generalize this model by using K ≥ 2 low-dimensional subspaces, and the resulting complexity of this generalization inspired intensive research over the past decade (e.g. [5,15,23,[25][26][27] and references therein). Another popular assumption is that the data is supported on a manifold M of dimension d isometrically embedded in R D [3, 8-11, 13, 14, 18, 19, 22, 24, 28]. Under this structural assumption, one can (locally) represent the point on the manifold M in the basis of the tangent space which has d elements. The GMRA procedure of [2,6] provides a multiscale, computationally efficient way to construct such bases, find the corresponding coordinates for a given point, and that paper demonstrates approximation guarantees in the infinite sample limit.
The main goal of this work is to obtain tight non-asymptotic bounds for the reconstruction error achieved by the GMRAbased method under the assumption that the distribution Π of the data possesses certain geometric properties. In particular, we verify the conditions for our guarantees for the case when Π is comparable with uniform distribution on a smooth submanifold of R D (or a thin tube around a submanifold), and the error rate primarily depends on the intrinsic dimension of the manifold. Detailed discussions and proofs are provided in [20].

Definitions and Notation
The following notation is used throughout the paper: · 2 denotes the standard Euclidean norm in R D , Proj V stands for the orthogonal projection on a linear subspace V ⊆ R D , dim(V ) for its dimension and V ⊥ for the orthogonal complement. Given a matrix A ∈ R D×D , tr (A) denotes its trace, · -its operator norm and A T -its transpose.
A GMRA with respect to the probability measure Π consists of a sequence of (nonlinear) operators {P j : R D → R D } j≥0 . For each "resolution level" j ≥ 0, P j is uniquely defined by a collection of pairs of subsets and affine projec- k=1 form a Borel-measurable disjoint partition of R D and P j,k (x) := c j,k + Proj Vj,k (x − c j,k ), where c j,k ∈ R D and V j,k are defined as follows. Let E j,k stand for the expectation with respect to the conditional distribution dΠ j,k (x) = dΠ(x|x ∈ C j,k ).
In other words, c j,k is the conditional mean and V j,k is the subspace spanned by eigenvectors corresponding to d largest eigenvalues of the conditional covariance matrix Given such a {(C j,k , P j,k )} where I{x ∈ C j,k } is the indicator function of the set C j,k . It was shown in [2] that if Π is supported on a smooth, compact d-dimensional submanifold M → R D , and if the partitions {C j,k } then, for any x ∈ M, x − P j (x) 2 ≤ C(M)2 −2j for all j ≥ j 0 (M). This means that the operators P j provide an efficient "compression scheme" for x ∈ M, in the sense that every x can be well-approximated by a linear combination of at most d + 1 vectors. In practice, Π is unknown and we only have access to the training data X = {X 1 , . . . , X n }, which are assumed to be i.i.d. with distribution Π. In this case, operators P j are replaced by their estimators k=1 is a suitable partition of R D obtained from the data, X j,k = C j,k ∩ X , |X j,k | is the cardinality of X j,k , (2.2) We call these P j the empirical GMRA.

L 2 (Π) error bounds for GMRA
We now state the sufficient conditions on the distribution Π and the partition {C j,k } N (j) j=1 which imply our L 2 (Π) error bound for the GMRA approximation. Suppose that for all j min ≤ j ≤ j max the following hold: x − y 2 ≤ θ 2 2 −j .

Special case: noisy manifolds
We now show that there is a large class of interesting distributions Π and a data-dependent partitioning algorithm such that (A1), (A2), and (A3) are satisfied with high probability given a sufficiently dense sample. Noisy manifold models: Let M be a smooth, compact ddimensional submanifold of R D . We recall the definition of the reach [16], an important global characteristic of M. Let Then reach(M) := sup{r ≥ 0 : M r ⊆ D(M)}. Assume that 0 < σ < τ . We shall say that the distribution Π satisfies the (τ, σ)-model assumption if there exist a compact, smooth submanifold M → R D without boundary and with reach τ such that supp(Π) = M σ , Π is absolutely continuous with respect to U Mσ -the uniform distribution on M σ , and the Radon- Data-dependent partitions: Our partitioning scheme is based on the cover trees [4] data structure. We will construct a cover tree for the collection S = {X 1 , . . . , X n } of i.i.d. samples from the distribution Π with respect to the Euclidean distance d(x, y) := x − y 2 . Assume that T j := k=1 is the set of nodes on level j of the cover tree T for S. Define the indexing map k(x) := argmin 1≤k≤N (j) x − a j,k 2 (the ties are broken by choosing the smallest value of k), and partition R D into the Voronoi regions C j,k = {x ∈ R D : k j (x) = k}. Let j(n, t) be the maximal j ∈ Z + such that n ≥ C(Π)2 jd (jd + t).

Discussion
Theorem 3.2 demonstrates that distributions supported on noisy manifolds interact in an inherently low-dimensional manner with partitions obtained via the cover trees algorithm. This is true as long as the resolution level of the partition remains above the noise level.
Theorem 3.1 supplies us with an error bound for the empirical GMRA provided that the distribution Π and the partition {C j,k } interact in an inherently low-dimensional manner. Pointwise bounds may also be obtained. This gives a rate in the noiseless case comparable to the minimax rate in [8], but with respect to the Hausdorff distance. A detailed discussion is left to [20].
For practical purposes, Theorem 3.2 may be combined with Theorem 3.1 as follows: given an i.i.d. sample {X 1 , . . . , X 2n } from Π, use the first n points {X 1 , . . . , X n } to construct the partition C j,k using the cover tree algorithm, and then the remaining {X n+1 , . . . , X 2n } are used to obtain P j (see 2.2). This makes our method entirely data-dependent.

Introduction
Four dimensional cone beam computed tomography (4D CBCT) of the free breathing thorax is important for imageguided radiation therapy (IGRT). Mainly two families of methods have been proposed to handle the problem: respiratory motion compensation, in which the motion of organs during breathing is estimated and used to deform the volume during reconstruction, and respiration-correlated reconstruction, in which several volumes are computed, each one from a subset of the projections that has been acquired during the correct respiratory phase. Advances in compressed sensing have been used to improve respiration-correlated reconstruction by enforcing spatial regularity constraints [1,2,3,4], but only a few recent methods exploit the strong correlation between successive respiratory phases [5,6,7,8]. 4D ROOSTER [8], which was recently proposed for 4D cardiac CBCT, is one such method. However, it came without a rigorous proof of convergence. In this paper, we propose a reformulation of 4D ROOSTER using the Chambolle & Pock primal-dual optimization scheme [9]. The convergence of this reformulated 4D ROOSTER is therefore guaranteed.

The original 4D ROOSTER method
This algorithm assumes that a rough segmentation of the patient's heart is available, and that movement is expected to occur only inside this segmented region. The method consists in iteratively enforcing five different constraints in an alternating manner. It starts by minimizing a quadratic data-attachment term α R α S α x − p α 2 2 , with α the projection angle, x a 4D sequence of volumes, R α the forward projection operator at angle α, S α a linear interpolator, and p α the measured projection at angle α. This data-attachment term is minimized by conjugate gradient. Then the following regularization steps are applied sequentially: positivity enforcement, averaging along time outside the heart region, spatial total-variation denoising, and temporal total-variation denoising. This constitutes one iteration of the main loop, the output of which is fed back to the conjugate gradient minimizer for the next iteration. This algorithm offers no convergence guarantees. In [10,11], using the theory of non-expansive mappings, it is only proved that if the main iteration has at least one fixed point, the algorithm converges to one of them. However, we show in this paper that each step of this method can be interpreted as a proximal operator [8], thus little effort is required to make it fit into the Chambolle & Pock framework. We remind that for f : R N → R a closed convex function, the proximal operator of f is defined as prox f (v) = arg min Rx − p 2 2 + λ 2 x T Vspace with x the 3D+t sequence of volumes, with in total N voxels, p the set of measured projections, with in total P pixels, R : R N → R P the forward projection operator, . T Vspace the spatial total-variation norm, . T Vtime the temporal totalvariation norm, R N + the set of sequences of volumes in which all voxels have non-negative values, ROI the set of sequences of volumes in which all voxels outside the heart have equal values, i R N + (x) and i ROI (x) their respective convex indicator functions, and λ 2 , λ 3 > 0 two parameters weighing the relative importance of the terms. We adopt the same formalism as in [12], and therefore use the following notations: Because the cost function has more than two terms, we have to reformulate the problem into a search in R 3N by iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
We define the total dimension W = P + 2N , the convex functions F and G and the linear operator K: With these new notations, we can write the cost function as a sum of two functions arg min which the Chambolle-Pock algorithm applies. In compact form, the Chambolle-Pock algorithm is as follows: If we expand this for our problem and simplify it since x 1 = x 2 = x 3 and y 1 = y 2 = y 3 , we obtain: → Positivity and ROI Four new operators appear in this formulation. R * : R P → R N is the back projection operator. To express the other ones, we use the relation between the proximal operator of a function and the proximal operator of its adjoint: The proximal operator for the spatial and temporal TV norms can be computed iteratively using the Chambolle algorithm [13], simply by using either the spatial or the temporal gradient and divergence operators And the proximal operator of γF * 2 (and similarly γF * 3 ) can be obtained from the same relationship as previously With this new formulation, the modified 4D ROOSTER is guaranteed to converge. The theory of proximal algorithms also helps find the suitable parameters based on the norm of R, while they can only be determined empirically in the original method.

Perspectives
4D ROOSTER can be applied to respiration-correlated 4D CBCT in order to reduce the level of streak artifacts in the reconstructions. Preliminary results are provided in Figure 1 and 2. They have been obtained using the original 4D ROOSTER method, and are compared with a standard ECG-gated FDK reconstruction. The modified 4D ROOSTER method presented in this paper will be implemented in the near future using the RTK framework [14]. Abstract-In large scale spatial surveys such as the forthcoming ESA Euclid mission, images may be undersampled due to the optical sensors sizes. Therefore, one may consider using a superresolution (SR) method to recover aliased frequencies, prior to further analysis. This is particularly relevant for point source images which provide direct measurements of the instrument point spread function (PSF). Super-resolution is a classical problem in image processing literature. It can be formulated as an inverse problem and sparsity constraint may be applied. We show that using an appropriate wavelet dictionary, this provides significant improvements over existing super-resolution methods, especially on the PSF.

Introduction
The weak lensing is one of the most promising tools to probe the dark matter distribution in the universe. The idea is to infer, from billion galaxies images, the shapes distortions due to dark matter gravitational lensing and then constraint the dark matter mass distribution. The Euclid mission will provide the data for such a purpose [1]. Nevertheless, galaxies images will be distorted due to the telescope PSF. It's critical to take accurately this PSF into account. In practice, isolated stars provide measurements of the PSF at different locations in the field of view. But these stars images will be aliased in Euclid, given the CCD sensor sizes. Nevertheless, several low resolution (LR) versions of the same PSF will be available thanks to the multiple exposures planned during the survey. This naturally leads to consider a SR approach to restore an accurate PSF.
SR techniques generally adopt the following observation model. We suppose that the desired high resolution (HR) image of size d 1 p 1 × d 2 p 2 is written in lexicographic order in a vector x = (x 1 , ..., x q ) T where q = d 1 p 1 d 2 p 2 and d 1 and d 2 are respectively the line and column downsampling factors. Let y k = (y k1 , ..., y kp ) T denote the k th LR observation written in lexicographical order, with p = p 1 p 2 and k = 1...n. So we have: where M k is a warp matrix, B k is a blur matrix, D k is a decimation matrix and n k is a noise column vector. The matrices M k traduce the motions of the observations relatively to each others and generally need to be estimated. The matrices B k account for different blurs (the atmosphere blur, which is particularly considerable for ground based telescopes, the system optics blur, the imaging system shaking, pixel response etc.). And D is simply a downsampling operator which is assumed to be the same for all the observations. In our case, we are interested in estimating the telescope PSF, in other terms the telescope contribution to the blur. We will consider that the telescope is steady so that the blur is stable in time. The Euclid telescope will be out of the atmosphere so that under our hypothesis, the blur only consists of the telescope PSF. Thus, the index k may be dropped for the blur operator. Moreover, we will consider that the PSF varies slowly in the field so that for small motions between observations, the warp and blur matrix may be commuted. So we can rewrite the model 1 as wherex = Bx and W k = DM k and we will be interested in recoveringx with x being assumed to be a point source.
In the next section, we will describe some SR techniques dedicated to astronomical images, but one may refer to [2] for more details on SR various framework in general image processing literature.

SR techniques in astronomy
The most simple super-resolution method is certainly the shiftand-add method. It performs in three steps : first the images are upsampled to the target resolution. Then they are shifted on a common grid and averaged. It has been used in astronomy for a long time, particularly for ground based telescopes. Despite its simplicity, it has been shown that this method provides an optimal solution in the sense of the maximum likelihood (ML) with additive white gaussian noise (WGN) and when only pure translation motions are considered. The interpolation operator should be D T (from Eq. 1) which comes down to a simple zero-padding interpolation, and the averaging should be performed with proper weighting [3]. Now let consider the SR method implemented in PSFEX-TRACTOR, an open source software more and more used in astronomy community for PSF modeling. This software and the associated documentation may be found on the website http://www.astromatic.net/.The desired HR image is now a matrix X of size dp × dp and the k th observation Y k is a matrix of size p × p with k = 1...n. The coordinates of the centroid of the k th observation are noted (i k , j k ). Assuming that the images are bandlimited, the samples of the LR images can be interpolated from the desired HR image, thanks to the well-known Shannon sampling theorem. In theory, this interpolation should involve a 2D sine cardinal kernel which will be approximated with a finite support function. Let h(., .) denote such a function. The estimate of the sample (i, j) of the k th observation is given by: Then we can define the cost function iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
where f k accounts for possible luminosity differences and needs to be estimated as well as σ 2 k which is related to the local background variance and the pixel values uncertainty.
PSFEXTRACTOR uses a Lanczos4 interpolant for h. Finally, PSFEXTRACTOR minimizes the cost function defined as where ∆ = X − X (0) and X (0) is an median image computed from the LR observations. This second term is meant for regularization purposes in ill-conditioned or ill-posed cases (n < d 2 ).

Sparse regularization
Now we introduce the SPRITE algorithm(SParse Recovery of InsTrumental rEsponse) which as stated before basically consists in adding a sparse penalty in the recovery process. As previously, our optimization variable will be the difference between the desired HR image and a first guess. But instead of using a quadratic regularization term, we impose the analysis coefficients of this difference to be sparse in a chosen wavelet dictionary. Let note Φ such a dictionary so that the wavelet coefficients are given by α x = Φx, then such a prior can be integrated by solving with ∆ defined as in 5 and X (0) calculated using a shift-andadd. A similar l 1 penalty has already been applied for SR. An example may be found in [4], where the cost function is minimized using variants of the alternating direction method of multipliers (ADMM). Besides advantages of such approaches over quadratic regularizers have been proved in many related problems in image and signal processing. It has been shown in [5] that under certain conditions on Φ the proximity operator of . l1 • Φ is well defined and can be estimated iteratively so that a forward-backward algorithm can be used to solve 6 ( [6]). Nevertheless, we opted for a simple iterative shrinkage scheme. LetΦ denote a reconstruction operator associated with Φ. The optimization variable is updated as follows: Proceeding this way amounts to alternately minimizing the J 1 (∆+X (0) ) and γ α ∆ l1 with α ∆ = Φ∆. This approach is inspired from the MCA decomposition algorithm in [7]. In our experiments, this gave qualitatively equivalent results as when we used the rigorously estimated analysis proximal operator, with a radical gain in complexity. Furthermore, we used the scheme introduced in [8] for FISTA to accelerate the convergence. In practice, the parameter γ, which is used as a threshold, is set at 3σ where σ is the estimated noise level at each considered wavelet scale. We used an undecimated wavelet transform with a 7/9 biorthogonal filter bank which is quite convenient since the data are nearly isotropic. Besides, we observe that reweighting the sparse penalty term as suggested in [9] tends in general to improve the final solution. Thus, the algorithm includes one "reweighted" iteration in which for each wavelet coefficient α ∆i , the threshold γ is now replaced with γ 1+|α ∆i is the wavelet coefficient previously calculated at the position i after thresholding. Thus in the reweighted iteration the thresholds will be lower on the stronger coefficients.

Numerical experiments
For Euclid mission, the actual sampling frequency is about 0.688 times the Nyquist frequency( [10]). Therefore we will target an upsampling factor of 2, which gives a sufficient bandpass to recover the high frequencies. Our data set consists of simulated optical point spread functions generated with the software Zemax, using Euclid telescope parameters. Our data set contains simulated PSF measurements on a regular 18 × 18 grid on the field of view (324 PSF). The original PSF are downsampled to almost Nyquist frequency. For each PSF, 4 randomly shifted "copies" are generated and downsampled with a factor 2 in lines and columns to achieve Euclid resolution . Different levels of WGN are added. We define the signal level as its empirical variance calculated in a 50 × 50 patch centered on the PSF main lobe. For each algorithm, we used these 4 images to reconstruct the original PSF. For quality assessment, we used the average error between the reference PSF ellipticity parameters and the super-resolved PSF ellipticity using each algorithm and these errors standard deviations. These ellipticity parameters definitions can be found in [10] and they are critical in weak lensing measurements.
As it can be seen on the figures below, smaller errors and errors dispersion are achieved with SPRITE algorithm, especially at low SNR.

Conclusion
We introduced the SPRITE algorithm which may be seen as an improvement over the super-resolution algorithm used in the software PSFEXTRACTOR. The latter method was implemented iteratively and augmented with a sparse penalty. This iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
way we achieved significantly better accuracy in terms of ellipticity error.

Introduction
In compressed sensing (CS) [1], random linear measurements are used to efficiently recover an unknown but assumed sparse signal. We propose the use of quasi-random samples as a better way of sampling. As an initial illustration of the benefits of non-random sampling methods, we reconstruct the Lena image from quasi-random point samples and compare it with the usual pseudo-random point sample reconstruction. As a regularization of the inversion, we use the total variation (TV) prior, i.e., we impose sparsity of the local image gradient.

Quasi-random sampling
The law of large number says that any sequence of randomly distributed numbers converge to a uniform distribution. This assumption is however valid only for very large numbers. When a strictly limited budget of samples should be selected for image reconstruction, it might be advantageous to rely on deterministic constructions ensuring more local uniformity in the point distribution, in contrast to clusters of points generated by true (or computer-generated pseudo-) randomness [2]. Elements of the fractional golden ratio sequence G s (i) with given seed constant s ∈ [0, 1) are given by the fractional part of the sum between s and an integer multiple of the golden ratio: where {t} is the fractional part of the real number t. Note that the conjugate golden section τ = 1 φ = φ − 1 can be subtitued to φ since only fractional parts are retained. Hence, computing G s (i + 1) given G s (i) only requires an addition and a bit mask.
Such sequences are especially interesting for distributing coordinates on the 1D unit range. A key corollary of the strong irrationality of φ (and τ ) is that the fractional parts of integer multiples will not align on any regular grid. In previous I. Loris is a research associate of the Fonds de la Recherche Scientifique-FNRS (Belgium). This work is sponsored by the iMind ICON project 3DUS.

Pseudo−random sampling
Quasi−random sampling Figure 1: Voronoi diagram for the first 512 points of a pseudo-random and quasi-random sequence. Spurious clusters emerges in the random distribution, while the quasi-random distribution ensures a more evenly uniform coverage.
work [3], the coordinates of this one-dimensional sequence are mapped to higher dimensions, using the invertible Hilbert space filling curve mapping H : R → R 2 . The Hilbert space filling curve H (t) = (x, y) , t ∈ [0, 1) defines a nested recursive grid and a locally-preserving traversal order of grid elements. Plugin-in G s (·) in H(·), we obtain the following sets of N uniformly distributed points: (2) Figure 1 shows the Voronoi tessellation for two sets of N = 512 point samples, using pseudo-random and the proposed quasi-random strategy with s = 0. Variability of Voronoi cell's areas indicates increases in local discrepancy.

Total variation reconstruction
Let S be the sampling operator S : R N1×N2 → R N : u → y = Su that corresponds to choosing a certain N pixels from the N 1 × N 2 image u. Based on samples y of some given image (y = Su orig ), the reconstructed image u rec minimizes the total variation and coincides with u orig in the chosen samples, i.e.: u rec = arg min Su=y TV(u).
Here the total variation TV(u) is defined in terms of the local differences (image gradients) D(u) ≡ (D h u, D v u) ∈ R 2 : with D h u N1,j = D v u i,N2 = 0. Thus, the total variation TV(u) in expression (3) is given explicitly by iTWIST'14: international -Traveling Workshop on Interactions between Sparse models and Technology.
Quasi−random recon. Original Figure 2: Original Lena image and reconstructions from 50.000 samples using pseudo-random and the proposed quasi-random algorithm. Visual quality and signal to noise ratio are improved when using the quasi-random method.
Total variation was introduced in imaging in 1992 [4]. This penalty promotes sparse gradients in the image. In other words, it promotes a piecewise constant reconstruction. In order to compute the solution to the convex minimization problem (3), we use the following iterative scheme: where the horizontal (h) and vertical (v) projection P λ is and it is understood that it is applied in each pixel separately. The algorithm (6) converges for positive step sizes τ 1 < S −2 2 and τ 2 < D −2 2 . In this case, the spectral norm of S is S 2 = 1 and the spectral norm of D is D 2 = 2 √ 2. The convergence speed of this algorithm depends on the positive parameter λ, but its choice is not critical. The algorithm depend on three auxiliary variables (ū, v and w). All variables are initialized (n = 0) by zero. A proof of convergence can be found in [6]. Other methods can be found e.g. in [5].

Results
In order to compare the two sampling methods, we constructed a numerical experiment in which an image is sampled at a limited number of pixels. The whole image is then reconstructed using the TV penalty. Figure 2, first row, shows the results of an initial comparison of the proposed quasi-random sampling method and standard pseudo-random sampling. Both reconstructions are based on 50.000 distinct point samples of the Lena image. In both cases we used the TV prior to regularize the inversion from 50.000 samples to a dense image of 512 × 512 pixels. Figure 2, second row, represents a close-up of the the first row images. More visual detail is discernable in the reconstruction from quasirandom samples, than from pseudo-random samples. The algorithm was run with up to 1000 iterations. The final peak signal to noise ratio (PSNR) of the reconstruction from quasi-random samples is 29.5dB whereas the final PSNR for the reconstruction from pseudo-random samples is 28.5dB. Furthermore, seven experiments were performed with 5000, 10.000, 20.000, 25.000, 30.000, 40.000 and 50.000 samples. Figure 3 plots the PSNR for the seven pairs of reconstructions of the Lena image. In both case the PSNR increases monotonically with the number of samples taken, but the quasirandom samples give better reconstruction results in all cases.

Conclusion
The law of large numbers says that randomly distributed points are expected to tend towards a uniform distribution. Experiments demonstrate that if the number of point samples is limited, spurious random clustering impairs sensing efficiency since some regions of the sampling domain are more sparsely populated. Fortunately, deterministic quasi-random methods exist for providing evenly distributed but still incoherent points.
Our first results on a natural image were encouraging and future work will quantify the potential of the method on a practical medical imaging setup. Another possibility consists of sampling the wavelet or Fourier domains in quasi-random points, as opposed to the usual pseudo-randomly placed samples, as is often done in applications of the compressed sensing framework. Abstract-Regularization plays a pivotal role when facing the challenge of solving ill-posed inverse problems, where the number of observations is smaller than the ambient dimension of the object to be estimated. In such settings, the general approach is to solve a regularized optimization problem, which combines a data fidelity term and some regularization penalty that promotes the low-dimensional/simple structure underlying the vector to recover. This work provides a general framework to capture this low-dimensional structure through what we coin piecewise regular gauges which are convex, non-negative, closed, bounded and positively homogeneous functions. This class of regularizers encompasses many popular examples such as the 1 norm, 1 − 2 norm (group sparsity), nuclear norm, as well as several others including the ∞ norm or the total variation semi-norm. Our main result presents a unified sharp analysis of exact and robust recovery of the low-dimensional subspace model associated to the object to recover from partial measurements. 1

Regularization by Gauges
We consider the forward model y = Φx 0 + w , where y ∈ R Q is the vector of observations, w ∈ R Q stands for the noise, and Φ ∈ R Q×N is a linear operator which maps the N -dimensional signal domain onto the Q-dimensional observation domain. We here consider solutions to the regularized optimization problem x ∈ Argmin x∈R N 1 2 ||y − Φx|| 2 + λJ(x), (P λ (y)) or its noiseless counterpart, when w = 0, x ∈ Argmin x∈R N J(x) subject to Φx = y.
(P 0 (y)) Here, J is a convex, non-negative, continuous (hence bounded) and positively homogeneous function. This class of regularizers J include many well-studied ones in the literature. Among them, one can think of the 1 norm used to enforce sparse solutions [10], the discrete total variation semi-norm [9], the 1 − 2 norm to induce block/group sparsity [14], the nuclear norm for low-rank matrices [5,4], or bounded polyhedral gauges [13].
Assuming that the gauges enjoy some piecewise regularity properties, our goal in this paper is to provide a unified analysis of exact and robust recovery guarantees of the subspace model underlying the object x 0 by solving (P λ (y)) from the partial measurements y. As a by-product, this entails a control on ||x − x 0 ||, where x is the unique minimizer of (P λ (y)).

Piecewise Regular Gauges
Associated to any bounded gauge, we define the model subspace at a point x as follows.
Definition 1 (Model Subspace). For any vector x ∈ R N , let S x = aff ∂J(x) be the affine hull of the subdifferential at of J at x, and e x the orthogonal projection of 0 ontoS x . We denote S x =S x − e x = span(∂J(x) − e x ) and T x = S ⊥ x .
T x is coined the model subspace of x associated to J.
When J is Gâteaux-differentiable at x, i.e. ∂J(x) = {∇J(x)}, e x = ∇J(x) and T x = R N . On the contrary, when J is non-smooth at x, the dimension of T x is strictly less that N , and the regularizing gauge J essentially promotes elements living on this model subspace.
Definition 2 (Antipromoting Gauge). Let x ∈ R N \ {0} and f x ∈ ri ∂J(x), where ri is the relative interior. The antipromoting gauge associated to f x is the gauge J • fx whose sublevel set at 1 is ∂J(x) − f x .
Toward the goal of studying robust recovery by solving (P λ (y)), we introduce a subclass of bounded gauges J for which the mappings x → e x , x → P Sx f x and x → J • fx exhibit local regularity.
Definition 3 (Piecewise Regular Gauge). A bounded gauge J is said to be piecewise regular (PRG for short) at x ∈ R N , with respect to the mapping x → f x ∈ ri ∂J(x) and Γ a coercive gauge bounded on T x = T = S ⊥ , if there exist three non-negative reals ν x , µ x , τ x and a real ξ x such that for every x ∈ T with Γ(x − x ) ν x , one has T x = T x (i.e. x ∈ T ) and Γ(e x − e x ) µ x Γ(x − x ). (1) Several examples are provided in Table 1. It turns out that the set of piecewise regular gauges is closed by summation and precomposition by a linear operator. Proposition 1. Let J and J two PRG and D a linear operator from R P to R N . Then J + J and J • D * are PRG.

Model Selection and Robustness
For some model subspace T , the restricted injectivity of Φ on T plays a central role in the sequel. This is achieved by imposing that Beside condition (C T ) stated above, the following Identifiability Criterion will play a pivotal role.   In terms of zero input symbol recovery probability P (z) r = 1 − p (z) , it is easy to show that the optimal degree distribution equals Ω(x) = x d * , where d * ≈ (log 1 1−α ) −1 . Optimal degree distributions for average input symbol recovery probability P r = 1 − p can be obtained as described in [1] and [8].

Adaptive Sudocodes
Consider the following modification to the Sudocodes scheme: Modification 1: After generating output symbol y i = 0, remove input symbols x j ∈ N (y i ) from consideration while generating following output symbols.
Modification 2: Due to the first modification, the sparsityfactor α decreases as the number of generated output symbols increase. Thus we continuously update the optimal degree d * = d * (α) before each output symbol is generated.
We analyze the Adaptive Sudocodes by following the evolution of measurement matrix with the process of generating output symbols (see Fig. 2). Starting from the initial sparsityfactor α 0 = K/N , the process runs through a sequence of increasing α-values {α 1 , α 2 , . . .} at which the optimal d * -values decrement. The total of m i measurements are generated using degree d * i during which the sparsity-factor increases from    α i to α i+1 . By tracking the parameters (dimensions and degree distributions) of the sequence of measurement matrices Φ i , i ≥ 0, and using Lemma 1, we are able to track the evolution of asymptotic recovery probabilities of LM1 reconstruction. The complete analysis is omitted due to lack of space, however, we note that due to approximations of degree distributions, the resulting analysis is not exact (see example below).
Example 3.1. Fig. 3 shows asymptotic recovery probability for Adaptive Sudocodes for initial sparsity-factor α 0 = 0.05. Simulated results for both Sudocodes and Adaptive Sudocodes scheme are provided for K = 50 and N = 1000 demonstrating that Adaptive Sudocodes performs better with lower complexity since the optimal degree d * decreases from d = 20 down to d = 2 at M/N = 0.35. We note that our approximate analysis matches well the simulated performance except that it becomes conservative in the region of large M/N -values.

Comments and Future Work
The rateless codes that apply Ω(x) with constant average degree µ = Ω ′ (1) are affected by upper bound on recovery probability that scales as e −µ M N . In LT codes, this is solved by using degree distributions whose average degree scales as O(log N ) [1], for the price of increased encoding/decoding complexity of O(N log N ). More advanced solution called Raptor codes preserves linear encoding/decoding complexity by using high-rate precoding combined with the constant average degree Ω(x) [8].
The Sudocodes scheme suffers from the same error-floor problem due to a constant-degree Ω(x) employed. This is solved by employing second phase where non-sparse measurements are combined with matrix inversion to recover small fraction of remaining input symbols [2]. However, instead of "postcoding", it is more instructive to use the Raptor-idea of "precoding" the message by adding small number of additional precoded input symbols followed by Sudocodes scheme with constant-degree Ω(x). This design is part of our ongoing work.
Finally, extending Adaptive Sudocodes to LM2 recovery has strong potential for further performance improvements and knowledge extraction from previous measurements. Note that in LM1 recovery, zero input symbols are learned only from zero-valued output symbols, while in LM2 recovery, the additional recovery rule provides new possibilities for zero input symbols identification. However, asymptotic analysis of LM2 algorithm seems to be considerably more involved [6]. The design and analysis of Adaptive Sudocodes for LM2 recovery is part of our ongoing work.