Sébastien Bourguignon

Jordan Ninin

Hervé Carfantan

Member, IEEE Marcel Mongeau

Exact Sparse Approximation Problems via Mixed-Integer Programming: Formulations and Computational Performance

Keywords: sparse approximation, 0-norm-based problems, optimization, mixed-integer programming, deconvolution

Sparse approximation addresses the problem of approximately fitting a linear model with a solution having as few non-zero components as possible. While most sparse estimation algorithms rely on suboptimal formulations, this work studies the performance of exact optimization of 0-norm-based problems through Mixed-Integer Programs (MIPs). Nine different sparse optimization problems are formulated based on 1, 2 or ∞ data misfit measures, and involving whether constrained or penalized formulations. For each problem, MIP reformulations allow exact optimization, with optimality proof, for moderate-size yet difficult sparse estimation problems. Algorithmic efficiency of all formulations is evaluated on sparse deconvolution problems. This study promotes error-constrained minimization of the 0 norm as the most efficient choice when associated with 1 and ∞ misfits, while the 2 misfit is more efficiently optimized with sparsity-constrained and sparsity-penalized problems. Then, exact 0-norm optimization is shown to outperform classical methods in terms of solution quality, both for over-and underdetermined problems. Finally, numerical simulations emphasize the relevance of the different p fitting possibilities as a function of the noise statistical distribution. Such exact approaches are shown to be an efficient alternative, in moderate dimension, to classical (suboptimal) sparse approximation algorithms with 2 data misfit. They also provide an algorithmic solution to less common sparse optimization problems based on 1 and ∞ misfits. For each formulation, simulated test problems are proposed where optima have been successfully computed. Data and optimal solutions are made available as potential benchmarks for evaluating other sparse approximation methods.

I. INTRODUCTION

A. Sparse estimation for inverse problems

The problem of sparse representation of data y ∈ R N in a dictionary H ∈ R N × R Q consists in finding a solution x ∈ R Q to the system y = Hx with the fewest non-zero components, i.e., with the lowest sparsity level. In sparse approximation, in order to account for noise and model errors, the equality constraint is relaxed through the minimization of the data misfit measure y -Hx , where • generally stands for the standard Euclidean norm in R N . Such sparsest representation and approximation problems are essentially combinatorial. Finding the best K-sparse solution (the solution with K non-zero components) is usually considered This work was partially supported by the French Groupement de Recherche ISIS (Information, Signal, Image, viSion) through its "young researchers" program. One of the authors of this work has been supported by French National Research Agency (ANR) through JCJC program (project ATOMIC n • ANR 12-JS02-009-01).

too difficult in practical large-scale instances. Indeed, the brute-force approach that amounts to exploring all the Q K possible combinations, is computationally prohibitive. In the abundant literature on sparse approximation, much work has been dedicated to the relaxation approach that replaces the 0 -"norm" sparsity measure, x 0 := Card{i | x i = 0}, with the 1 norm x 1 := i |x i |. Many specific convex optimization algorithms have been proposed in the past decade, see for example [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF], [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. In addition, conditions were established for which both the 0 and the relaxed 1 problems yield the same solution support (the set of non-zero components). These mostly rely on a low sparsity level assumption and on structural hypotheses on the matrix H, such as low correlation of its columns (see [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF] and references therein). Alternatively, greedy algorithms build a sparse solution by iteratively adding non-zero components to an initially emptysupport solution [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF], [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF]. More complex forward-backward methods [START_REF] Miller | Subset Selection in Regression[END_REF], [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] may show better performance in practice but with higher computation time. Tree-search-based methods also try to improve the classical greedy algorithms using heuristics to reduce the complexity of exhaustive combinatorial exploration (see e.g., [START_REF] Karahanoglu | A* orthogonal matching pursuit: Best-first search for compressed sensing signal recovery[END_REF] and references therein). Other support exploration strategies maintain the desired sparsity level at each iteration, and perform local combinatorial exploration steps [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]. Optimality proofs for all such strategies also rely on very restrictive hypotheses [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF], [START_REF] Soussen | Joint k-step analysis of orthogonal matching pursuit and orthogonal least squares[END_REF]. More " 0 -oriented" approaches were proposed, e.g., by successive continuous approximations of the 0 norm [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed 0 norm[END_REF], by descent-based Iterative Hard Thresholding (IHT) [START_REF] Herrity | Sparse approximation via iterative thresholding[END_REF], [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF] and by penalty decomposition methods [START_REF] Lu | Sparse approximation via penalty decomposition methods[END_REF]. However, without additional assumptions on H, one can only prove that the solution found is a local minimum of the optimization problem. Moreover, for IHT, optimality conditions suffer from the same restrictions as the aforementioned greedy methods [START_REF] Herrity | Sparse approximation via iterative thresholding[END_REF], [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF].

In many inverse problems, the model y Hx results from the discretization of an intrinsically continuous physical model. A typical example is sparse deconvolution, where Hx models the convolution of a spike train (in one-dimensional signals) or of point sources (in imaging problems) by the impulse response of the acquisition device [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], [START_REF] O'brien | Recovery of a sparse spike time series by 1 norm deconvolution[END_REF], [START_REF] Mendel | Some modeling problems in reflection seismology[END_REF]. A similar problem concerns nonlinear parameter identification, where parameters are discretized on arbitrarily thin grids [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], [START_REF] Bourguignon | A sparsity-based method for the estimation of spectral lines from irregularly sampled data[END_REF], [START_REF] Tang | Sparse recovery over continuous dictionaries -just discretize[END_REF] and estimation amounts to finding a sparse solution to a linear problem of high dimension. In such cases, the columns of H can be highly correlated, so no optimality guaranty can be obtained for greedy and 1 -norm-based methods.

Similar problems also arise for variable selection in machine learning and statistics [START_REF] Miller | Subset Selection in Regression[END_REF], where the set of features (the columns of H) is not designed to satisfy any recovery property.

The aforementioned problems essentially focus on the correct estimation of the support of x. In deconvolution, for example, support identification corresponds to detection and localization of the sources. Since the true sparsity measure is indeed the 0 norm, a global optimum of 0 -based formulations is more likely to yield exact support identification than approximate solutions. Consequently, our interest focuses on optimization methods for 0 -norm-based criteria providing optimality guarantees. Such exact approaches are usually discarded, based on the argument that sparse optimization problems are NP hard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. It is also commonly considered that exact optimization amounts to combinatorial exploration of all possible supports, which is nothing but a worst-casescenario argument. Note however that, in order to reduce the number of explored supports, Tropp and Wright [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF] mentioned the possible use of cutting-plane methods, which are one of the basic elements of resolution of the mixed-integer programs explored in the present paper.

Here, we focus on sparse optimization occurring in certain inverse problems with moderate size yet with a complexity sufficient to make the usual methods fail in estimating the sparsest solutions. Examples include spike-train deconvolution in ultrasonic nondestructive testing (NDT) [START_REF] O'brien | Recovery of a sparse spike time series by 1 norm deconvolution[END_REF] or Geophysics [START_REF] Mendel | Some modeling problems in reflection seismology[END_REF], sparse hyperspectral unmixing [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], and spectral analysis with short data sets [START_REF] Klein | A sparsity-based method for the analysis of magnetic fluctuations in unevenly-spaced Mirnov coils[END_REF]. A first objective of this contribution is to show the viability of exact resolution approaches for such problems.

B. Global optimization via Mixed-Integer Programming

We focus on Mixed-Integer Programs (MIP), that is, optimization problems involving both continuous and integer variables. Such problems are well suited to 0 -norm-based optimization, since the 0 norm naturally introduces a binary decision variable for each component (zero or non-zero?). In this paper, MIP refers to the minimization of linear or quadratic criteria subject to linear or quadratic inequality constraints. It is commonly claimed that, in the past fifteen years, a factor 10 9 was gained in the required computing time for solving such problems. This gain is due in (roughly) equal parts to (i) hardware improvement, (ii) progress in the resolution of linear programs, and (iii) implementation efficiency of advanced mathematical techniques [START_REF] Bixby | A brief history of linear and mixed-integer programming computation[END_REF]. Therefore, as it will be shown in this paper, some exact approaches can now be advantageously used to address the moderate-size, yet difficult, applications enumerated at the end of Section I-A.

To our knowledge, the first MIP reformulation of a sparse optimization problem is proposed in [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF]. However, the authors argue that the assumption that |x| is upper bounded, which is required for the MIP reformulation, leads to computational inefficiency. Therefore, they choose to consider only a related problem: the maximum feasible subsystem problem, for which exact solutions can be found only for very small instances (N = 16, Q = 32) and no result is given concerning the MIP approach. A similar formulation with binary variables appears in [START_REF] Tosic | Learning joint intensity-depth sparse representations[END_REF], but binary variables are replaced by continuous variables in [0, 1] in order to yield a convex problem, which is obviously not equivalent to the original one. In [START_REF] Sun | Recent advances in mathematical programming with semi-continuous variables and cardinality constraint[END_REF], some exact and approximate reformulations of 0 -based problems are surveyed. The authors deplore the inherent combinatorial difficulty of such MIP problems but no practical result is provided. Finally, in [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF], noisefree sparse representation problems are formulated as MIP. Here, we address the noisy case, which opens perspectives to different possible formulations of the optimization problem. Establishing MIP reformulations for such problems, studying their computational efficiency, investigating properties of optimal solutions and comparing them with the results of standard methods are the core of this paper.

C. Objectives of the paper

This paper shows that different possible reformulations of sparse approximation problems can be tackled by MIP solvers. Sparse approximation is intrinsically a bi-objective optimization problem, where both the sparsity measure and the data misfit measure are optimized. In inverse problems, it is usually formulated through the optimization of a weighted sum of the two terms. However, constrained formulations (involving one criterion to be minimized and the other subject to a constraint) may also be well suited to MIP reformulations, which are constrained programs by nature. Therefore, we study the efficiency of MIP solving techniques applied to the following three formulations:

• minimize the 0 norm under a bounded-data-misfit constraint, • minimize the data misfit under an 0 -boundedness constraint, • minimize a weighted sum of the two criteria. Additionally, we also consider non-quadratic data misfit measures, which may be appropriate if the error term y -Hx is non-Gaussian. Moreover, piecewise-linear alternatives to the 2 norm x 2 := i |x i | 2 , may also prove to be more attractive computationally, because MIPs essentially rely on the resolution of linear subproblems. In particular, the 1 and ∞ ( x ∞ := max i |x i |) norms are easily linearized; formulations based on those norms naturally boil down to optimization problems involving linear inequality constraints or linear objective functions. Therefore, we also consider formulations involving 1 and ∞ misfits, for which much fewer algorithms have been proposed.

Our work establishes MIP reformulations of nine different sparse approximation problems, which are all evaluated in terms of computational efficiency, depending on the sparsity level and on the noise level. Then, these formulations are compared in their ability to identify the exact support in the presence of noise, depending on the noise statistical distribution. Our experimental results additionally show that the classical methods are far from reaching acceptable results in such cases, whereas the 0 -norm formulations do yield more satisfactory solutions-but with much higher computing time. Note that experiments in both [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF] and [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF] involve random matrices, which are more likely to satisfy the conditions ensuring the optimality of 1 -norm-based and greedy approaches. In most ill-posed inverse problems, the columns of H are highly correlated, so that such conditions certainly do not hold.

The remainder of the paper is organized as follows. Section II introduces nine optimization formulations of the sparse approximation problem, and discusses their statistical interpretation and the structure of the solution sets. MIP reformulations are then established in Section III. Then, basic elements concerning the resolution of MIP problems are given in Section IV. Experimental results in Section V are dedicated to the evaluation of computational costs. Section VI compares the solutions obtained through MIP optimization with those of classical sparse approximation algorithms, on both overdetermined and underdetermined sparse deconvolution problems. Simulations in Section VII evaluate the support identification performance of p -misfit-constrained formulations as a function of p and of the noise statistical distribution. Finally, a discussion is given in Section VIII.

II. SPARSE OPTIMIZATION PROBLEMS

In sparse approximation, both the sparsity of the solution and the fidelity of its corresponding data approximation are optimized. Therefore, the generic sparse approximation problem, which we are interested in, is the following unconstrained bi-objective optimization problem:

min x ( x 0 , y -Hx p ), (1)
where p is either 1, 2 or ∞. This section presents nine formulations of this problem and discusses their statistical interpretation and the structure of the sets of solutions.

A. Taxonomy

Various mono-objective optimization problems can be formulated to address the bi-criterion problem [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. For p ∈ {1, 2, ∞}, the bounded-error problems read P 0/p : min x x 0 s.t. y -Hx p ≤ α p , and the sparsity-constrained problems read

P p/0 : min x y -Hx p s.t. x 0 ≤ K p ,
where α p and K p are user-defined threshold parameters. Finally, the penalized problems read P 0+p : min x µ p x 0 + y -Hx p p , for p ∈ {1, 2}, and P 0+∞ : min

x µ ∞ x 0 + y -Hx ∞ ,
where µ p and µ ∞ are user-defined penalty parameters.

In this paper, we propose a reformulation of each of these problems as MIPs. To the best of our knowledge, P 0/∞ is the only sparse approximation problem for which a MIP reformulation is mentioned [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF]. Remark that the sparse representation case (noise-free data), which was recently tackled via MIPs in [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF] with equality constraint y = Hx, is the special case of P 0/p with α p set to 0, for which P 0/1 , P 0/2 , and P 0/∞ are obviously equivalent. Recall that the sparsity-based inverse problems considered here are sparse approximation problems: data are always contaminated by measurement noise and the model may be inexact, so that α p = 0.

Choosing one of the nine formulations and the values of the parameters (α p , K p or µ p ) amounts to selecting some particular solutions among the wide variety of Pareto-optimal solutions of problem [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. Note that, for a given p misfit, no equivalence between the three problems P 0/p , P p/0 and P 0+p can be obtained because the 0 norm is not convex. In particular, solutions in the non-convex part of the Pareto frontier cannot be reached by solving the penalized formulation [START_REF] Das | A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems[END_REF].

B. Statistical interpretations and parameter tuning

In practice, one has to choose one among the nine optimization problems and must set a value for the corresponding parameter. Such choices can be based on statistical arguments.

The p data-misfit measures, with p ∈ {1, 2, ∞}, can be interpreted in terms of likelihood functions. Let p be the statistical distribution of the additive noise term = y -Hx. The likelihood function is defined as: L(y; x) := p (y -Hx). If noise samples n , n = 1, . . . , N are independent and identically distributed (i.i.d.) according to a centered Gaussian distribution, thenlog L(y; x) is proportional to y -Hx up to an additive constant. Similarly,log L(y; x) is proportional to y -Hx 1 (up to an additive constant) if n are i.i.d. according to a centered Laplace distribution. Such a heavy-tailed distribution assumption may be appropriate in the presence of impulsive noise [START_REF] Marks | Detection in Laplace noise[END_REF], [START_REF] O'brien | Recovery of a sparse spike time series by 1 norm deconvolution[END_REF]. The ∞ misfit is connected to an i.i.d. uniform noise distribution assumption. Suppose that is uniformly distributed on [-a, a] N , for some given a > 0. Then, the likelihood function is constant for any x such that y -Hx ∞ ≤ a, otherwise it is zero. Consequently, x = arg min x y -Hx ∞ is a maximumlikelihood estimator if noise samples are uniformly distributed on [-a, a], for any a > y -H x ∞ [START_REF] Boyd | Convex Optimization[END_REF]Ch. 7.1]. In this case, which arises for example when accounting for quantization noise, ∞ data fitting may be a relevant choice-see [START_REF] Clason | L ∞ fitting for inverse problems with uniform noise[END_REF] for both theoretical and numerical arguments.

Consequently, P p/0 is a sparsity-constrained maximum likelihood estimation problem with the aforementioned corresponding noise distribution assumption. In a Bayesian setting, P 0+p defines a Maximum A Posteriori (MAP) estimate, where the 0 term results from a Bernoulli-Gaussian prior model with infinite variance Gaussian distribution [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF]. Note that, within such a MAP interpretation, solving P 0+∞ reduces to solving P 0/∞ . Indeed, sincelog L(y; x) is constant if y -Hx ∞ ≤ a and equals +∞ otherwise, minimizing µ ∞ x 0 -log L(y; x) amounts to minimizing x 0 subject to y -Hx ∞ ≤ a. Finally, P 0/p considers a maximal tolerance on the approximation error and cannot be interpreted as a maximum likelihood or MAP estimation problem.

The choice of the formulation also depends on available prior information about the considered practical problem. If reasonable bounds on the acceptable approximation error can be inferred, e.g., from the knowledge of the signal-to-noise ratio or from a desired approximation quality, then P 0/p may be preferred. In particular, the parameter α p can be fixed according to the statistics of p , which can be obtained for any noise distribution (analytically or numerically). If the sparsity level is fixed or can be upper bounded, e.g., in a compression context, then P p/0 may be appropriate. In P 0+p , the parameter µ p trades off between the noise level and the sparsity level.

With the previous MAP interpretation, for p ∈ {1, 2}, it is an explicit function of the noise variance and of the rate of non-zero values in the Bernoulli process. Therefore, tuning µ p requires more information than tuning the parameters of the two other formulations. When too little prior information is available, a practical solution consists in computing optimal solutions corresponding to different parameter tuningswhatever the considered formulation-and then selecting a posteriori the most appropriate one, according to some expert supervision or to model order selection criteria [START_REF] Stoica | Model-order selection: A review of information criterion rules[END_REF].

C. Structure of the solution sets

We investigate hereafter the structure of the solution sets of the different problems, for fixed values of the corresponding parameters α p , K p and µ p . In the following, an optimal support refers to a set of indices which supports at least one optimal solution. The 0 norm is a piecewise constant function, where each value is attained on a finite number of supports. Hence, for any problem P defined in Section II-A, the set of minimizers can be defined as the finite union of sets of minimizers on each optimal support: if S denotes the set of optimal supports, then:

Arg P = s∈S Arg P s ,
where P s denotes the restriction of problem P to the support s.

Let us characterize the solution set of P s . We assume that the sparsity level of all solutions is lower than N and that the matrix H satisfies the Unique Representation Property (URP) [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm[END_REF], that is, any N columns of H are linearly independent. For any x supported by s, x 0 is constant, hence P s p/0 and P s 0+p are solved by minimizing y -H s x s p , where x s (respectively, H s ) collects the non-zero components in x (respectively, the corresponding columns of H). Thanks to the URP, H s has full column rank and, for p = 2, P s admits a unique solution (the least-squares solution). Consequently, the solution sets of P 2/0 and P 0+2 are both (finite) unions of singletons. The p norms for p = 1 and p = ∞ are not strictly convex, therefore one can only claim that min xs y -H s x s p is attained on a convex set, such that H s x s lies in an p -sphere in dimension K s -1, centered at y. Consequently, for p = 1 and p = ∞, the solution sets of P p/0 and P 0+p are (finite) unions of convex sets of the form y -H s x s p = constant . Now, consider the solution set of P s 0/p , which is formed by all vectors x s such that y -H s x s p ≤ α p . In the particular case where min xs y -H s x s p = α p , it comes from the previous arguments that the solution set is a singleton for p = 2, and a convex set for p = 1 and p = ∞. But, in the most frequent case where min xs y -H s x s p < α p , the solution set x s y -H s x s p ≤ α p is such that H s x s lies in an p -ball of dimension K s , centered at y, and the solution set of P 0/p is a finite union of such sets. Consequently, for a given p, the solution set of P 0/p is generally "larger" than the solution sets of P p/0 and P 0+p . In particular, the minimizers of these last two problems may be unique. For example, with some additional assumptions on the data y and on the matrix H, the solution of P 0+2 is unique [START_REF] Nikolova | Description of the minimizers of least squares regularized with 0 -norm. Uniqueness of the global minimizer[END_REF]. On the contrary, the minimizer of P 0/p is certainly not unique, except in very specific cases.

III. MIXED-INTEGER REFORMULATIONS

In this section, we establish the reformulations of optimization problems P 0/p , P p/0 and P 0+p for p ∈ {1, 2, ∞}, as Mixed-Integer Linear Programs (MILPs), Mixed-Integer Quadratic Programs (MIQPs) or Mixed-Integer Quadratically Constrained (linear) Programs (MIQCPs).

A. Definitions of MILP, MIQP and MIQCP

The general form of an MILP is

min v c v, subject to (s.t.)          A in v ≤ b in , A eq v = b eq , l b ≤ v ≤ u b , v j ∈ Z, ∀j ∈ I, where v ∈ R J is the vector of optimization variables; c ∈ R J defines the linear objective function; b in ∈ R Pin , b eq ∈ R Peq , A in ∈ R Pin × R J and A eq ∈ R Peq × R J define
the inequality and equality constraints; l b and u b ∈ R J are respectively the vectors of lower and upper bounds of the optimization variables; I is the index set corresponding to the components of v that are constrained to be integer-valued.

An MIQP has the general form:

min v 1 2 v Fv + c v, s.t.          A in v ≤ b in , A eq v = b eq , l b ≤ v ≤ u b , v j ∈ Z, ∀j ∈ I,
where F is a J × J matrix.

Finally, the form of an MIQCP that is of interest in this paper is:

min v c v, s.t.                  1 2 v Bv + d v ≤ e, A in v ≤ b in , A eq v = b eq , l b ≤ v ≤ u b , v j ∈ Z, ∀j ∈ I,
where B is a J × J matrix, d ∈ R J and e ∈ R.

B. Equivalent reformulation techniques

We now present standard reformulation techniques that enables to express each of the nine optimization problems introduced in Section II-A as an MILP, an MIQP or an MIQCP, without any approximation.

1)

Boundedness assumption and "big-M " reformulation of the 0 norm: For each q = 1, 2, . . . , Q, let us introduce an additional binary optimization variable b q , such that b q = 0 ⇔ x q = 0.

(2)

Then, the non-linear sparsity measure x 0 is equal to the linear term 1 Q b (= q b q ), where 1 Q is the Q-dimensional allones column vector. The logical constraint (2) must however be translated into (in)equality constraints compatible with MIP, ideally through linear constraints. One standard way to achieve this (see e.g. [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF], [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF]) is to assume that a solution x of the problem under consideration satisfies the following constraints for some sufficiently large pre-defined value M > 0:

-M 1 Q < x < M 1 Q . ( 3
)
This assumption supposes that the problem admits bounded optimal solutions, which is not restrictive in our practical applications. The parameter M has to be large enough so that x ∞ < M at any desirable optimal solution x. On the other hand, the bound M must be as tight as possible in order to improve computational efficiency; therefore tuning the value of parameter M may be a critical issue [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF]. In the problems addressed in this paper, satisfactory results are obtained with a rather simple empirical rule discussed in Section V-B. Note that specific lower and upper bounds for each component x q could also be advantageously considered [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF] if corresponding prior information is available.

The reformulations of 0 -norm-based constraints and objective functions are obtained through the two following respective lemmas.

Lemma 1: Considering the boundedness assumption (3),

x 0 ≤ K ⇔            ∃b ∈ {0; 1} Q such that Q q=1 b q ≤ K, (i) -M b ≤ x ≤ M b. (ii)
Proof: The ⇒ implication is straightforward by considering b defined by Equation (2). Now, let b ∈ {0; 1} Q satisfy (i) and (ii), and suppose x 0 > K. From (ii), one has

(b q = 0) ⇒ (x q = 0), that is, (x q = 0) ⇒ (b q = 1). Hence b q = 1 for at least K + 1 indices q, which contradicts (i). Consequently, x 0 ≤ K.
Lemma 2: Considering the boundedness assumption (3),

min x∈F x 0 ⇔ min x∈F b∈{0,1} Q Q q=1 b q s.t. -M b ≤ x ≤ M b,
where F represents the feasible domain of the problem under consideration.

Proof: Similar to that of Lemma 1. Such a reformulation technique is commonly referred to as "big-M " reformulation. Remark finally that another reformulation of the cumbersome logical constraint (2) consists in introducing the equality constraint x q (b q -1) = 0. However, the latter is a bi-linear constraint, typically less interesting in terms of computation time for off-the-self MIP solvers than linear constraints [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF].

2) Reformulation of the 1 data misfit measure: The 1 misfit term can be written linearly as y -Hx 1 = n w n , with additional constraints w n = |y nh r n x|, n = 1, . . . , N , where h r n denotes the nth row of H. Then, these constraints can be relaxed (exactly) by the linear inequalities: -w ≤ y-Hx ≤ w, with column vector w = [w 1 , . . . , w N ] , thanks to the two following lemmas.

Lemma 3:

min x∈F y -Hx 1 ⇔ min x∈F ,w∈R N n w n s.t. -w ≤ y -Hx ≤ w. Proof: Let f n (x) = y n -h r n x.
The following optimization problems are equivalent:

P a : min x∈F n |f n (x)| P b : min x∈F ,w∈R N n w n s.t. |f n (x)| = w n , ∀n P c : min x∈F ,w∈R N n w n s.t. |f n (x)| ≤ w n , ∀n.
Indeed, P a ⇔ P b is trivial. In order to show that P b ⇔ P c , one can simply remark that if an optimal solution (x , w ) of P c is such that |f n0 (x )| < w n0 for some index n 0 , then one can straightforwardly construct a better feasible solution for P c , which yields a contradiction. Lemma 3 will be used to obtain a MIP reformulation of P 1/0 and P 0+1 , which involve the 1 -misfit term in the objective function. For P 0/1 , which involves the 1 -misfit term as a constraint, we use the following lemma:

Lemma 4: Let (x , w ) solve the optimization problem:

P d : min x∈R Q ,w∈R N x 0 s.t.    n w n ≤ α 1 -w ≤ y -Hx ≤ w.
Then, x is a solution of P 0/1 . Proof: Suppose that (x , w ) solves (P d ) and let w n := |yh r n x |, ∀n. Then, n w n ≤ n w n ≤ α 1 and (x , w ) is a solution of:

P e : min x∈R Q ,w∈R N x 0 s.t.    n w n ≤ α 1
|y nh r n x| = w n , ∀n. Indeed, since P d is a relaxation of P e -the feasible set of P e is a subset of that of P d -, its optimal value, x 0 , is a lower bound for the optimal value of P e . The solution (x , w ) is clearly feasible for P e and it attains the lower bound x 0 . Hence, it is optimal for P e . Finally, P e is clearly equivalent to P 0/1 . Let us remark finally that problems P d and P e are not strictly equivalent because they are not minimized by the same couple of vectors (x, w), but they share the same solution set for x.

3) Reformulation of the ∞ data misfit measure: P 0/∞ naturally brings linear inequality constraints as

y -Hx ∞ ≤ α ∞ ⇔ -α ∞ 1 N ≤ y -Hx ≤ α ∞ 1 N .
For both P ∞/0 and P 0+∞ , involving the ∞ norm in the objective function, one can simply introduce an additional scalar variable z such that minimizing y -Hx ∞ amounts to minimizing z under the constraints

-z1 N ≤ y -Hx ≤ z1 N .
Here again, one can easily show that an optimal solution of the original problem will necessarily satisfy y -Hx ∞ = z.

C. Mixed-integer programming reformulations

Given the reformulation techniques of Section III-B, the nine problem reformulations proposed in Table I are straightforward.

IV. MIP RESOLUTION: BASIC ELEMENTS

Mixed-integer programming problems are not easy problems. MILP problems are already NP-hard [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. As a consequence, reducing sparse approximation problems to MILP, MIQP or MIQCP problems does not per se reduce the complexity. Nevertheless, such MIP reformulations not only open up possibilities to prove the optimality (or quantify the sub-optimality) of solutions, but also allows one to benefit from decades of exponential progress in terms of required computing time to solve a given MIP problem. This progress does not simply reflect the doubling of computing power every 18 months, it is also a consequence of fast progress in both the theory and practice of linear programming and discrete optimization (duality, numerical linear algebra, interiorpoint methods, semi-definite positive relaxations, branch and bound/cut/price methods, decomposition approaches, global optimization, etc.) Once the sparse approximation problem is recast as a MIP problem, then state-of-the-art off-theshelf software can be used, such as BARON, COUENNE, CPLEX, GloMIQO, GUROBI, MOSEK or Xpress-MP-see for example [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF] and references therein. We chose to use CPLEX [START_REF] Bixby | A brief history of linear and mixed-integer programming computation[END_REF] because it is unanimously considered among the best MIP solvers. CPLEX has been developed over the last thirty years and includes the best strategies developed by the MIP community. Moreover, it is freely available for research and teaching purposes.

The main method behind the CPLEX MIP solver is a branch-and-cut algorithm. Globally, it implements a branchand-bound strategy (i.e., a tree-structured implicit enumeration algorithm) based on successive continuous relaxations of the integer variables [START_REF] Lodi | Mixed integer programming computation[END_REF]. Each branch generates a subproblem by fixing some integer (in our case, binary) variables, and a set of branches is constructed corresponding to the different possible configurations of such variables. Then, the aim is to discard (huge) parts of the remaining combinatorial tree by lower bounding the objective function on the corresponding subproblems. To obtain such lower bounds, a continuous relaxation of each subproblem is formed by relaxing the integer variables. Linear constraints, such as Gomory cutting planes [START_REF] Gomory | Outline of an algorithm for integer solutions to linear programs[END_REF], are added to each generated (continuous relaxation) subproblem, so that the continuous solution converges to an integer solution. Such cutting planes remove part of the feasible domain of the subproblem that does not contain any integer solution. This approach amounts to attempting to construct the convex hull of the set of integer feasible solutions of each subproblem. CPLEX incorporates several techniques in order to improve performance, such as constraint propagation techniques [START_REF] Benhamou | Revising hull and box consistency[END_REF], linear algebra techniques [START_REF] Bixby | A brief history of linear and mixed-integer programming computation[END_REF] and heuristic techniques to find rapidly a good integer solution. Doing so, parts of the research space are eliminated only if it is proved that they do not contain the global minimum.

The best current integer solution provides an upper bound of the global minimum of the entire problem. The solution of the current relaxed (continuous) subproblem gives a lower bound of the global minimum of the current subproblem with integer constraints under consideration. The worst solution of all relaxed subproblems-the one that achieves the lowest lower bound-gives a certified lower bound of the global minimum of the entire problem. If such a lower bound is attained by the best current integer solution, then a global minimizer is found and optimality is proved. Otherwise, the entire process is iterated by creating new branches. The algorithm converges towards such a certified optimum in a finite number of steps. If the solver reaches the time limit, the duality gap (the difference between the best solution found and the certified lower bound of the global minimum) provides a measure of sub-optimality of the solution found. Note that tree-search based greedy algorithms such as in [START_REF] Karahanoglu | A* orthogonal matching pursuit: Best-first search for compressed sensing signal recovery[END_REF] also rely on treebased (local) exploration and on lower-bounding the objective function on sets of solutions, which is used inside a greedy procedure. Therefore, they do not come with any optimality guarantee.

Remark that in most sparsity-inspired signal processing problems, the matrix H satisfies specific properties that are exploited for efficient computations. In particular, if H represents a redundant set of transforms based on multi-scale representations such as wavelets [START_REF] Mallat | A wavelet tour of signal processing: The sparse way[END_REF], matrix-vector products Hv and H v (where v is some given vector of an appropriate dimension) can be computed by fast algorithms using Fast Fourier Transforms (FFT). In 1D (respectively, 2D) deconvolution problems, H is a Toeplitz (respectively, Toeplitz-block-Toeplitz) matrix, so that matrix-vector products can also be computed by FFT. No such matrix structure is exploited within the standard mixed-integer optimization algorithms, as they are implemented as general-purpose MIP solvers. The only matrix property that is exploited here is the possible sparsity of the matrix H to compute fast vector products. This is particularly the case for the finite impulse response (FIR) deconvolution problems that are considered hereafter.

V. EXPERIMENTAL RESULTS: OPTIMIZATION PERFORMANCE This section presents the test problems and the computational framework. Then, the computational efficiency of the nine MIP reformulations is studied.

A. Definition of test problems

The different MIP reformulations are evaluated on onedimensional sparse deconvolution problems. Such problems are typically encountered for example in ultrasonic NDT [START_REF] O'brien | Recovery of a sparse spike time series by 1 norm deconvolution[END_REF] Problem

Equivalent MIP reformulation MIP class Bounded-misfit formulations (parameters: α 2 , α 1 and α ∞ )

P 0/2 min x∈R Q ,b∈{0,1} Q 1 Q b s.t. x H Hx -2y Hx ≤ α 2 2 -y y -M b ≤ x ≤ M b MIQCP P 0/1 min x∈R Q ,b∈{0,1} Q ,w∈R N 1 Q b s.t.      1 N w ≤ α 1 -w ≤ y -Hx ≤ w -M b ≤ x ≤ M b MILP P 0/∞ min x∈R Q ,b∈{0,1} Q 1 Q b s.t. -α ∞ 1 N ≤ y -Hx ≤ α ∞ 1 N -M b ≤ x ≤ M b MILP Sparsity-constrained formulations (parameters: K 2 , K 1 and K ∞ ) P 2/0 min x∈R Q ,b∈{0,1} Q x H Hx -2y Hx s.t. 1 Q b ≤ K 2 -M b ≤ x ≤ M b MIQP P 1/0 min x∈R Q ,b∈{0,1} Q ,w∈R N 1 N w s.t.      1 Q b ≤ K 1 -w ≤ y -Hx ≤ w -M b ≤ x ≤ M b MILP P ∞/0 min x∈R Q ,b∈{0,1} Q ,z∈R z s.t.      1 Q b ≤ K ∞ -z1 N ≤ y -Hx ≤ z1 N -M b ≤ x ≤ M b MILP Penalized formulations (parameters: µ 2 , µ 1 and µ ∞ ) P 0+2 min x∈R Q ,b∈{0,1} Q µ 2 1 Q b + x H Hx -2y Hx s.t. -M b ≤ x ≤ M b MIQP P 0+1 min x∈R Q ,b∈{0,1} Q ,w∈R N µ 1 1 Q b + 1 N w s.t. -w ≤ y -Hx ≤ w -M b ≤ x ≤ M b MILP P 0+∞ min x∈R Q ,b∈{0,1} Q ,z∈R µ ∞ 1 Q b + z s.t. -z1 N ≤ y -Hx ≤ z1 N -M b ≤ x ≤ M b MILP TABLE I MIXED-INTEGER PROGRAMMING REFORMULATIONS OF NINE SPARSE APPROXIMATION PROBLEMS.
and in seismic reflection in Geophysics [START_REF] Mendel | Some modeling problems in reflection seismology[END_REF]. In the following, x is a K-sparse sequence in R 100 (i.e., Q = 100), with uniformly distributed spike locations, where the sparsity level, K, is varying. In order to avoid arbitrary small spike values, each non-zero amplitude is drawn as sign(u) + u, where u is a centered Gaussian sample with unit variance. The matrix H is the discrete convolution matrix corresponding to the 21-sample impulse response shown in Fig. 1 (left). With the boundary assumption that x is zero outside its domain, y is a 120-sample signal (i.e., N = 120). White noise is added with variable signal-to-noise ratio (SNR). In this section, noise samples are generated according to a centered normal distribution N (0, σ 2 ), with σ such that SNR dB = 10 log 10 Hx 2 2 /(N σ 2 ) . We name such problems SA SNR K . Note that they are slightly overdetermined (Q < N ), whereas typical sparse approximation problems deal with largely underdetermined systems. However, trivial inversion is not satisfac-tory here, because of the presence of noise and of the illconditioned nature of H (cond(H) 10 3 ).

We consider problems SA SNR K with K varying between 5 and 11, and SNR varying from +∞ (noise-free data) downto 10 dB. One example of data is given in Fig. 1 (right) for K = 7 and SNR = 20 dB. It illustrates the difficulty of sparse deconvolution problems arising for example in ultrasonic NDT [START_REF] O'brien | Recovery of a sparse spike time series by 1 norm deconvolution[END_REF]. The oscillating impulse response and the proximity of the spikes produce overlapping echoes in the available data. As the echoes cannot be distinguished visually, numerical techniques are required. All data and optimization results of this section are available online as supplementary multimedia material.

B. Machine configuration and implementation details

Optimization is run with IBM ILOG CPLEX V12.6.0 from a Matlab interface on a computer with eight Intel Xeon X5472 processors with Central Processing Units (CPU) clocked at 3 GHz. The maximum time allowed for each resolution is set to T max = 1000s. The other CPLEX parameters are set to their default value. For each problem, the "big-M " constant is set to M = 1.1x

C. Evaluation of computational costs

Each of the nine MIP reformulations of Table I is run for fifty random instances (for both spike distributions and noise realizations) of each SA SNR K problem. In this section, in order to ensure a fair comparison (in terms of computational efficiency), the parameters K p are set to the true sparsity level K, and parameters α p and µ p are set by trial and error, until the sparsity level of each solution equals K. Note that in our case, the matrix H has full column rank, hence 0 inequality constraints in P p/0 will yield the same results as if they were equality constraints. That is, when imposing x 0 ≤ K, all the K columns of matrix H contribute to the reduction of the data misfit. However, formulations with inequality constraints yielded lower computing times. Sparse representation problems (noise-free data, SNR = ∞) are addressed through P 0/p for p ∈ {1, 2, ∞}, with threshold α p = 10 -8 . Remark that in the noise-free case, no sparsityenhancing algorithm is indeed necessary, since the solution can simply be computed by least squares.

Average CPU times obtained for MIP reformulations are given in Table II. The figures on the left-hand side of each column is the time required to prove the global optimality of the solution found. The figures on the right-hand side indicate the time after which the support of the solution was found, which is generally much lower. The figures in parentheses indicate the number of instances for which optimality was not proved within 1000s.

All CPU times increase with the sparsity level K, but also with the noise level. In particular, for SNR = 10 dB and K = 11, for each formulation, optimality of the solutions was obtained in less than 1000s only on a fraction of the fifty instances. In order to explain such behavior, let us remark that the sparsity level (respectively, the noise level) increases the size of the feasible domain of P p/0 (respectively, of P 0/p ).

More generally, for all problems, if either the sparsity level or the noise level increases, then the branch-and-bound strategy becomes less efficient in discriminating between concurrent candidates and in eliminating some of them.

With SNR = 30 dB, the lowest CPU times are achieved by solving P 0+2 -although optimization did not terminate within 1000s for one instance of SA 30 11 . When the noise level increases, solving P 0/1 and P 0/∞ problems becomes more efficient computationally, and their superiority over other problems increases with both the sparsity level and the noise level. In particular, P 0/1 problems were always solved exactly in less than 1000s, except for (SNR = 20 dB, K = 11). Results are slightly not as good for P 0/∞ , where optimization did not terminate within 1000s for three more instances with SNR = 20 dB and K = 11. Note that even a small amount of noise severely degrades the resolution performance of problems P 0/p . Optimization of P 0/2 , the only problem with a quadratic constraint, is the most time-consuming among all proposed formulations. Non-linear constraints are known to make the MIP really difficult to solve [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF]. An element of explanation can be found by comparing the Lagrangians of formulations involving 2 misfits. Indeed, the Lagrangian of P 0/2 contains trilinear terms. On the contrary, the Lagrangians of P 2/0 and of P 0+2 are quadratic functions. Therefore, optimizing a linear function under quadratic constraints is more complex than optimizing a quadratic function under linear constraints.

For any p ∈ {1, 2, ∞}, solving problems P 0+p generally performs better than solving problems P p/0 at high SNR. On the contrary, as the noise level increases, sparsity-constrained formulations outperform penalized versions. For both formulations, which involve the data misfit in the objective function, using an 2 misfit measure is the most efficient choice, and both 1 -and ∞ -misfit optimizations behave similarly.

We also note a high dispersion of the required CPU times among the fifty realizations of each problem. For example, the average time for the resolution of P 0+1 on SA 20 9 problems was approximately 11s on forty-eight instances, whereas optimization did not terminate after 1000s on the two other instances. We also remark that, for P 0/2 , optimality was not proved within 1000s for two instances of the simplest test problem SA +∞ 5 . Two other instances of SA +∞ 5 also required a much higher CPU time than the others, which leads to an atypically high average time of 53s, reflecting once again the difficulty of the quadratically-constrained problem P 0/2 .

Finally, let us evaluate the CPU time of exhaustive combinatorial exploration for P 2/0 . Using notations introduced in Section II-C, for a given support s with K components, the minimizer of P s 2/0 has the closed-form expression x s = (H s H s ) -1 H s y. Then, the least-squares misfit value is computed by: y -H s x s 2 = y 2 -H s x s 2 . In practice,

x s can be computed by performing Cholesky decomposition of H s H s , so that one computation of the objective function mainly amounts to two K × K triangular inversions. The CPU time, denoted c K , of one such computation is estimated by averaging over 10 5 inversions. Then, neglecting the cost of Cholesky factorizations, the cost for testing all K-sparse solutions is extrapolated as 399| 96 (18) 18 | 4.1 19 | 2.1 29 |4.0 (1) 72 | 11 (2) 78 | 17 (9) 11 |2.7 (1) 60 | 13 (1) (3) 166| 21 (6) 244| 22 (11) 355|260 (20) 104| 16 (5) 179| 58 (13) 299| 38 (17) 485| 63 (34) 24 | 3.9 26 | 1.2 54 | 12 (1) 189| 28 (4) 220|218 (11) 96 | 25 (11) 230| 11 (18) 215| 32 (16) 50) 287 | 44 (12) 211 | 9.0 (9) 358| 26 (24) 670| 75 (42) 600|598 (39) 340| 34 (32) 300| 13 (43) 508|101 (44) 1500s for K = 5, 4 days for K = 7, and more than one year for K = 9. Problems based on both 1 and ∞ misfits require the resolution of a linear program for each support, therefore the corresponding exhaustive search yields still higher computational costs than in the 2 case. Consequently, exhaustive search cannot be considered as a practical solution, even for such moderate-size problems. In order to emphasize the ability of the MIP solver to remove important parts of the full combinatorial tree search, we give a last indicator. For P 2/0 with SNR = 30 dB and K = 9, for which all instances were successfully solved, the average number of combinations that were explicitly considered in the branch-and-bound procedure is about 4.5 10 5 , to be compared with the total number of 100 9

K = 11 -| -(
1.9 10 12 combinations.

VI. EXPERIMENTAL RESULTS: EVALUATION OF SOLUTIONS

We now compare the solutions obtained via MIP optimization with those of classical sparse approximation methods: Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], Single Best Replacement (SBR) [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], Iterative Hard Thresholding (IHT)1 [START_REF] Herrity | Sparse approximation via iterative thresholding[END_REF], [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF] and the minimization of the 1 -norm-penalized leastsquares criterion ( 1 -relax) using homotopy continuation [START_REF] Malioutov | Homotopy continuation for sparse signal representation[END_REF]. All algorithms are tuned so that all solutions have the correct number of spikes. Therefore, all methods are put in a favorable setting for evaluating their ability to retrieve the correct support of the solution.

Recall that for sparse deconvolution, none of the classical methods are theoretically guaranteed to solve the 0norm problem. As seen in the previous section, the MIP approach can compute an exact solution, but requires a larger computation time. In order to evaluate intermediate solutions found by the MIP solver, we consider several values of the maximum time allowed for each MIP resolution: T max =1s, 10s, 100s and 1000s. If the maximum time is reached, then the current solution is considered-which is the best solution found, without any optimality guarantee. The parameter M is tuned as explained in Section V-B when running optimization with T max = 1000s. Then, this value is also used for lower values of T max .

In this section, we consider two types of sparse deconvolution problems. The first problems are similar to those of Section V-A, with N = 120 and Q = 100, and are therefore slightly overdetermined. The second ones are underdetermined problems, the most frequent case in sparse approximation. Such problems may arise in high resolution sparse deconvolution. Indeed, the true spike locations are generally continuousvalued (e.g., representing times of flight of reflected waves), and the model Hx is a discrete approximation of a continuous integral equation. The discretization step is usually chosen equal to the sampling period of the data, thus H is a discrete convolution matrix. Such a choice may be too rough and can lead to erroneous spike locations [START_REF] Carcreff | Resolution enhancement of ultrasonic signals by up-sampled sparse deconvolution[END_REF]. In order to improve the model, one may consider an upsampled convolution model, where both the impulse response and the sparse sequence are upsampled by an integer factor, UF. Then, H can be viewed as the concatenation of UF discrete convolution matrices. Detail about the corresponding matrix structure can be found in [START_REF] Carcreff | Resolution enhancement of ultrasonic signals by up-sampled sparse deconvolution[END_REF]. In the following, we consider UF = 2, so that N = 120 and Q = 200. Recall however that the intrinsic difficulty of sparse deconvolution is mostly due to the ill-conditioned nature of matrix H and to the presence of noise, even in the overdetermined case.

Fifty random instances are run with UF = 1 and UF = 2, SNR = ∞ and SNR = 20 dB, and sparsity levels K varying between 5 and 11. In the noise-free case, for the MIP approach, only P 0/p problems are relevant. Therefore, solutions are computed via the optimization of P 0/∞ (with threshold α ∞ = 10 -8 ), which was the most efficient computationally (see Section V-C). All data and MIP optimization results of this section are available online as supplementary multimedia material.

In the noisy case, we consider P 2/0 , since the classical methods rely on 2 data misfit measures. Two quality indices are computed for each simulated data set:

• the exact recovery rate, which is the average number of simulations for which a given algorithm correctly locates the support of the true sequence; • the average spike distance between the estimated and the true sequences. To this aim, a distance similar to that used in neuroscience [START_REF] Van Rossum | A novel spike distance[END_REF] is defined: both estimated and true spike trains are convolved with a continuous-domain Gaussian kernel (with standard deviation chosen equal to the discretization step of the spikes). Then, the 2norm error between both signals is computed. Such a criterion is less severe than the exact recovery rate if the estimated spikes are slightly shifted with respect to the true ones, and also gives more importance to spikes with high amplitudes. Note that with SNR = ∞, the true sequence is a minimizer of P 0/p , with p = 1, 2, ∞. Hence the expected exact recovery rate for the MIP approach is 100%. On the other hand, in the noisy case, the minimizer of P 2/0 may not be obtained on the true support. Therefore, successful global optimization does not always produce exact recovery.

Results are summarized in Fig. 2. The top row shows the average exact recovery rate as a function of the number of spikes, and the bottom row plots the average spike distance to the true sequence as a function of the CPU time.

Let us first focus on the noise-free case (left columns). Recall that, in the overdetermined noise-free case, the solution can simply be computed by least squares. Simulations are still of interest, however, in order to compare the algorithms in an ideal and simple context. For all classical algorithms, the exact recovery rate is lower than 40% in the overdetermined case, and decreases as the sparsity level increases (top row). Their performance is still worse in the underdetermined case, where the exact recovery rate is close to zero, except for the simplest problems (K=5). Their average spike distance to the true sequence (bottom row) is logically lower for algorithms requiring more computation time. We note in particular the bad results obtained by IHT in the underdetermined case. Such bad performance of IHT was already attested when theoretical optimality conditions do not hold [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF]. This is particularly true in the underdetermined case, where the columns of H are strongly correlated. On the contrary, the MIP strategy correctly retrieves the support in nearly 100% of the noisefree instances, even with the computation time limited to 1s. Actually, only one instance led to erroneous support identification (for UF = 2 and K = 5), meaning that the solution was not found, even within 1000s. The MIP approach also gives an average spike distance close to zero, which means that both the supports and the amplitudes of the solutions have been correctly recovered, even in the underdetermined case, but with a larger computation time (from 0.03s to 0.2s for UF = 1, and from 0.1s to 20s for UF = 2). Note however that all classical algorithms are still much faster on such relatively small problems (between 10 -3 s and 10 -2 s).

In the more realistic noisy case (SNR = 20 dB), the results of the classical algorithms are very similar to those obtained in the noise-free case, both in terms of exact recovery rate, average spike distance and CPU time. In contrast, the MIP performance deteriorates, and the exact recovery rate quickly decreases as the number of spikes increases. Recall however that, in the presence of noise, the minimizer of P 2/0 may not retrieve the correct support: in the overdetermined case, for example, the MIP solver returns the optimal solution in less than 1000s in 94% of the instances, whereas the average exact recovery rate is much lower. However, it is still better than that of the classical methods with T max = 1s, and even much better if T max is increased to 100s or 1000s. The MIP approach also outperforms classical methods in terms of average spike distance, in particular if T max is high enough. In the overdetermined case, the average computing time ranges from 0.25s (for K = 5) to 350s (for K = 11) and, as mentioned earlier, global optimality was obtained in less than 1000s for most simulations. In the underdetermined case, however, an optimum was not proved to be found within 1000s in 51 % of the instances, that mostly correspond to the cases where K = 9 and K = 11. This analysis corroborates the results in Section V-C: the presence of noise strongly impacts the computing time of the MIP solver, and therefore the quality of the solutions obtained by early stopping.

A typical data set and estimated sparse sequences are shown in Fig. 3. It corresponds to the overdetermined case, the true sequence is 9-sparse and SNR = 20 dB. In this example, the MIP approach is the only algorithm that correctly identifies the support. Note that the resulting 2 misfit at the MIP solution is lower than the 2 norm of the noise.

VII. EXPERIMENTAL RESULTS

: RELEVANCE OF 1 -, 2 - AND ∞ -NORM DATA MISFITS
In this section, the impact of the data misfit measure (through 1 , 2 and ∞ norms) on the quality of the solution is studied, as a function of the noise distribution-as motivated by the discussion in Section II-B. To this aim, data are simulated in a manner similar to Section V-A and Fig. 1. The 7-sparse spike sequence x true is fixed, and 200 noise realizations are drawn, where noise samples are i.i.d. according to Gaussian, Laplacian and uniform distributions, successively. The SNR here is set to 15 dB. The three errorconstrained problems P 0/p are considered here. We focus on these formulations because, in practical cases, tuning the parameter α p requires less prior information than tuning the parameter K p for P p/0 or the parameter µ p for P 0+p . Indeed, for any given noise distribution p , the parameters For each algorithm, the residual norm and the average spike distance to the true sequence (ASD) are given.

α p , p ∈ {1, 2, ∞}, can be naturally estimated from a common statistical rule. More precisely, setting α p to a value satisfying Pr p ≤ α p | ∼p = 95 % amounts to considering that the approximation error cannot be distinguished from noise with probability 95%. Doing so, fair comparisons can be performed between the solutions of the three problems P 0/p . These three problems share the same objective function, hence they can be compared through their minimum value, that is, the estimated number of spikes. In order to evaluate the quality of the solutions, we also consider the support error e supp ( x, x true ) := || bb true || 0 , where b q (respectively, b true q ) equals 1 if x q = 0 (respectively, x true q = 0), and 0 otherwise. Optimization terminated in less than 1000s for all simulations. Fig. 4 shows the distribution of the 0 norms (left column) and of the support errors (right column) for the three noise distributions, and for the three p data misfits.

We first note that the correct value of the 0 norm is the most frequently reached in all cases, except if an ∞ misfit is used with Laplacian noise. When the estimated 0 norm is wrong, it is generally lower than the true value. Indeed, in most cases (statistically, in 95% of the cases from the definition of α p ), the noise realization satisfies p < α p . That is, the true sequence x satisfies the misfit-bound constraint. Thus, the threshold α p allows a higher approximation error than the noise in the data, which enables the possibility of solutions with fewer spikes. Consequently, in such cases, x 0 ≤ x 0 . This is particularly true if the noise distribution is heavy-tailed, where the p -norm of the noise samples may be much below the threshold α p , as show the 0 norm estimation statistics in the case of Laplacian noise. On the other hand, if p ≥ α p , then the true sequence does not correspond to a solution of the optimization problem, and y -H x p < y -Hx p . In Gaussian noise Occurrences (%) Support error e supp ( x, x true ) Fig. 4. Estimation results obtained for the deconvolution of a 7-sparse sequence with SNR = 15 dB, averaged over 200 noise realizations, in the case of Gaussian (top), Laplacian (center), and uniform (bottom) noise distributions. Error-constrained problems P 0/p are considered, with p = 2 (white), p = 1 (gray) and p = ∞ (black). Left: distributions of the 0 norm of the solutions for the three p misfits. For each problem, K indicates the average estimated 0 norm value. Right: distributions of the support errors. For each problem, e indicates the average support error between the estimated and the true sequences. such cases, one may have x 0 > x 0 . In our simulations, only very few instances led to such an overestimation of the 0 norm. As it could be expected, the lowest support errors are achieved by using the 2 (respectively, 1 and ∞ ) misfit in the case of Gaussian (respectively, Laplacian and uniform) noise. For each noise distribution, the corresponding misfit yields the smallest average support error, and more frequently achieves correct support identification-even if, for Gaussian or Laplacian noise, it is only obtained in a few cases (respectively, 15% and 11%). We also remark that switching from 2 to 1 misfits with Gaussian noise only slightly degrades the support identification performance, whereas optimization is computationally more efficient in the 1 case-see Section V-C. Much better support identification is achieved with uniform noise combined with the ∞ misfit, which yields exact identification in 90 % of the cases, whereas 1 and 2 data fitting achieve much worse results. Note that with both ∞ and 2 misfits, the 0 norm is correctly estimated in 98% of the cases. However, support recovery performance is much worse in the 2 case, as some spikes are misplaced. Such superiority of ∞ data fitting for uniform noise was already attested in [START_REF] Clason | L ∞ fitting for inverse problems with uniform noise[END_REF] in a non-sparse context. Fig. 5 displays typical results obtained for one particular realization of the noise process. For Gaussian (respectively, Laplacian and uniform) noise distributions, one example is shown such that 2 (respectively, 1 and ∞ ) data fitting yields a solution with the most frequent support error obtained among the 200 realizations. Note that, in each case, the solution shown corresponds to one solution of the considered optimization problem, that is, with the lowest 0 norm that satisfies the boundedp -misfit constraint. Recall indeed that, for most values of threshold parameters α p , problems P 0/p feature an infinite number of solutions-see Section II-C. Consequently, the presented solution is almost certainly not the solution with minimal p misfit. With Gaussian noise, the minimizer of P 0/2 has the correct 0 norm, but with two misplaced spikes, that leads to a support error equal to 4. With the minimizer of P 0/1 , two spikes are slightly misplaced, and a third one is not detected. The minimizer of P 0/∞ also has the correct 0 norm, but its spikes are very badly located. For Laplacian noise, the most frequent support error for the minimizer of P 0/1 is 3, which corresponds to one misplaced spike and the non-detection of one spike. Note that on the presented example, both minimizers of P 0/1 and of P 0/2 identify the same support, whereas the solution of P 0/∞ features only four spikes (among which one is erroneous). In the case of uniform noise, the solution of P 0/∞ correctly locates all spikes. The solution of P 0/1 misplaces one spike and misses another one, and the solution of P 0/2 is still worse, with three misplaced spikes-although with the correct sparsity level.

VIII. DISCUSSION

In this paper, nine sparse approximation problems involving the 0 norm were considered and reformulated as mixedinteger programs (MIP). Bounded-error, sparsity-constrained and penalized formulations were studied, involving p -norm data misfit measures, for p ∈ {1, 2, ∞}. Thanks to efficient dedicated MIP solvers, we demonstrated that moderate-size sparse approximation problems can be solved exactly, whereas exhaustive search remains computationally prohibitive for such instances. In particular, the use of a branch-and-bound strategy, coupled with efficient cutting planes methods, allows most combinations to be discarded without being evaluated.

Computational costs were evaluated on simulated difficult sparse deconvolution problems. Simulated data and corresponding optimal solutions are made available as potential benchmarks for evaluating other (potentially suboptimal) sparse approximation algorithms 2 . Our experiments show that misfit-constrained minimization of the 0 norm is the most efficient optimization formulation when associated with 1 and ∞ misfit measure. Conversely, the 2 misfit measure is advantageously used as an objective function, not as a Gaussian noise constraint. All CPU times increase with the number of nonzero components in the true solution, and also with the amount of noise in the data. Our encouraging numerical results tend to indicate however that such optimization formulations may be appropriate for tackling sparse approximation problems with several hundreds of unknowns, as long as the solution is highly sparse and/or the noise level is low. In particular, they do represent an alternative to 1 -norm-based and greedy methods for difficult estimation problems with highly correlated dictionaries, both of which are likely to fail. Simulations revealed, in particular, that exact solutions of 0 problems almost always achieve perfect support recovery for underdetermined noise-free problems, whereas classical methods perform relatively badly. In the presence of noise, the MIP solutions still outperform that of classical methods (in both over-and underdetermined cases), although the required computing time for obtaining exact solutions dramatically increases.

The 0 sparse approximation problem with 2 data misfit measure has been used in a huge quantity of works in signal processing, statistics, machine learning, etc. To the best of our knowledge, the methods presented in this paper are the only guaranteed-optimality alternatives to exhaustive search that do not rely on any strong assumption on the dictionary structure.

With the introduced MIP reformulations, we also proposed to solve exactly less common sparse optimization problems based on 1 and ∞ misfits. Such problems may be of interest from an informational point of view. Simulations illustrated this point, and confirmed a rather intuitive fact: choosing an p misfit with p = 2 (respectively, p = 1 and p = ∞) is relevant if the noise distribution is Gaussian (respectively, Laplacian and uniform) as far as support identification is concerned. In particular, with uniformly distributed noise, introducing an ∞ misfit constraint frequently achieves correct support identification, which is not the case for any other combination of data misfit and noise distribution.

Several points in the MIP reformulations could be considered in order to improve computational efficiency. First, as acknowledged in previous works on MIP reformulations of sparsity-based problems [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF], [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF], tuning the value, M , in the "big-M " reformulation impacts algorithmic performance. For a given problem, statistical rules may be used in order to infer reasonable M values. Then, new constraints in the optimization formulations may be added in order to reduce the feasible domain. For example, in [START_REF] Karahanoglu | A mixed integer linear programming formulation for the sparse recovery problem in compressed sensing[END_REF], an upper bound on the 0 norm of the solution sought is considered. Furthermore, many signal processing problems naturally involve linear constraints such as positivity or sum-to-one requirements. The proposed MIPbased approaches can easily be adapted to such cases, for which exact solutions can still be obtained. Adding such extra constraints may also contribute to reducing the computational time, whereas it generally penalizes the efficiency of classical (convex or greedy) sparse approximation algorithms. One may also consider directly the bi-objective optimization problem with multi-criterion optimization methods [START_REF] Ehrgott | Multicriteria optimization[END_REF] in order to propose a whole range of trade-off (sparsity vs. data fitting) solutions.

Global optimization of criteria involving structured sparsity would also be worth being studied, where (possibly overlapping) subsets of coefficients are jointly zero or non-zero. Such problems are generally tackled by convex optimization approaches involving mixed norms [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF] or by extensions of greedy algorithms [START_REF] Huang | Learning with structured sparsity[END_REF]. Both suffer from similar limitations than their scalar 1 -norm relaxation and greedy counterparts, as far as optimality with respect to the 0 -based problem is concerned. We believe that exact optimization of such problems through MIP should also be possible for moderatesize problems. For example, MIP-like formulations of some structured sparsity problems are shown in [START_REF] Tosic | Learning joint intensity-depth sparse representations[END_REF]-although the authors finally resort to (inexact) continuous relaxation of the binary variables-and in [START_REF] Baldassarre | Group-sparse model selection: Hardness and relaxation[END_REF], where specific structured sparsity problems defined through totally unimodular systems allow exact optimization in polynomial time.

Fig. 1 .

 1 Fig. 1. Example of sparse deconvolution data. Left: impulse response. Right: 7-sparse sequence x (circles) and noisy convolved data y, with SNR = 20 dB.

Fig. 2 .

 2 Fig.2. Performance of classical sparse approximation algorithms and MIP optimization on sparse deconvolution problems, as a function of the signal-to-noise ratio and of the size of the unknown vector. Overdetermined cases correspond to the standard deconvolution problem (N = 120, Q = 100). Underdetermined cases correspond to an upsampling factor UF = 2 such that N = 120 and Q = 200. Top row: exact recovery rate as a function of the sparsity level K. Bottom row: average spike distance to the true sequence as a function of the computation time. For each algorithm, the four marks correspond to K = 5, 7, 9 and 11. Note the log-scale and the different scalings of the time axes. The notation MIP Tmax corresponds to the MIP optimization strategy running in a maximum of Tmax seconds. All results are averaged over 50 random realizations.

1 3 SBR, r 2 26 Fig. 3 .

 32263 Fig.3. Data (top-left) and deconvolution results obtained by different sparse approximation algorithms in the overdetermined case (N = 120, Q = 100), with SNR = 20 dB and K = 9 spikes. Circles locate the true spikes. On the top-left panel, the dashed (respectively, solid) line represents the data y = Hx+ (respectively, the noise ). On the other panels, crosses show the obtained minimizer x and the solid line represents the residual r = y -H x. For each algorithm, the residual norm and the average spike distance to the true sequence (ASD) are given.

45 P0/ 2 : 0 P0/ 2 :

 45202 e = 3.17 P0/1 : e = 3.39 P0/∞ : e = 4.24x Support error e supp ( x, x true ) K = 6.7 P0/1 : K = 6.58 P0/∞ : K = 5.17

45 P0/ 2 :P0/ 2 :

 4522 e = 3.98 P0/1 : e = 3.04 P0/∞ : e = 5.75 x 0 Support error e supp ( x, x true ) e = 2.97 P0/1 : e = 3.69 P0/∞ : e = 0.43x 0

3Fig. 5 .

 5 Fig.5. Solutions of deconvolution problems P 0/p with Gaussian (top), Laplacian (center) and uniform (bottom) noises, for one particular noise realization, with SNR = 15 dB. Circles locate the true spikes. On the left column, the dashed (respectively, solid) line represents the data y (respectively, the noise ). On the three other columns, crosses show the obtained minimizer x and the solid line represents the residual y -H x.

TABLE II CPU

 II TIMES (IN SECONDS) OBTAINED BY THE NINE MIP REFORMULATIONS, AS A FUNCTION OF THE SIGNAL-TO-NOISE RATIO (SNR) AND OF THE SPARSITY LEVEL (K), AVERAGED OVER 50 INSTANCES OF EACH PROBLEM. THE FIRST NUMBER INDICATES THE CPU TIME FOR ESTABLISHING OPTIMALITY OF THE SOLUTION. THE SECOND NUMBER GIVES THE CPU TIME AT WHICH THE CORRESPONDING OPTIMUM WAS FOUND. NUMBERS IN PARENTHESES COUNT THE INSTANCES FOR WHICH OPTIMALITY WAS NOT PROVEN IN LESS THAN 1000S. FOR REFERENCE, THE COMPUTATION TIME OF CLASSICAL SPARSE APPROXIMATION ALGORITHMS ON OUR HARDWARE CAN BE FOUND IN FIG. 2.

We used T. Blumensath's implementation of IHT (program AIHT.m) available at: http://www.personal.soton.ac.uk/tb1m08/sparsify/sparsify.html

Matlab implementations of the nine formulations are available at http://www.irccyn.ec-nantes.fr/ ∼ bourguignon