Member, IEEE Gabriel Peyré

Best Basis Compressed Sensing

Keywords: Compressed sensing, best basis, sparsity, wavelet packets, cosine packets, bandlets

This paper proposes a best basis extension of compressed sensing recovery. Instead of regularizing the compressed sensing inverse problem with a sparsity prior in a fixed basis, our framework makes use of sparsity in a tree-structured dictionary of orthogonal bases. A new iterative thresholding algorithm performs both the recovery of the signal and the estimation of the best basis. The resulting reconstruction from compressive measurements optimizes the basis to the structure of the sensed signal. Adaptivity is crucial to capture the regularity of complex natural signals. Numerical experiments on sounds and geometrical images indeed show that this best basis search improves the recovery with respect to fixed sparsity priors.

I. INTRODUCTION

C OMPRESSED sensing is a new sampling strategy that uses a fixed set of linear measurements together with a non-linear recovery process. In order for this scheme to work with a low number of measurements, compressed sensing theory requires the sensed signal to be sparse in a given orthogonal basis and the sensing vectors to be incoherent with this basis. This theory of compressive acquisition of data has been proposed jointly by Candès, Tao and Romberg [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF], [START_REF]Decoding by linear programming[END_REF] and Donoho [START_REF] Donoho | Compressed sensing[END_REF], [START_REF]For most large underdetermined systems of linear equations, the minimal 1 norm solution is also the sparsest solution[END_REF].

This paper extends the compressed sensing recovery by switching from a fixed orthogonal basis to a tree structured dictionary of orthogonal bases. The adaptivity of a best basis representation increases the sparsity of sounds and geometrical images, which in turns makes the compressed sensing recovery more efficient. The tree structure of the dictionary is used in a fast iterative thresholding algorithm that estimates both the signal or the image to recover and the best basis that optimizes the sparsity of the representation.

A. Compressed Sensing

In a series of papers, Candès, Tao and Romberg [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF], [START_REF]Decoding by linear programming[END_REF] and Donoho [START_REF] Donoho | Compressed sensing[END_REF], [START_REF]For most large underdetermined systems of linear equations, the minimal 1 norm solution is also the sparsest solution[END_REF] have proposed the idea of directly acquiring signal in a compressive form. Instead of performing the acquisition with a high sampling rate and then compressing the data in an orthogonal basis, the signal is rather projected on a reduced set of linear vectors. The compressibility of the signal is only exploited during the reconstruction phase, where one uses the sparsity of the signal in an orthogonal basis.

Compressed sensing acquisition of data might have an important impact for the design of imaging devices where data acquisition is expensive. Duarte et al. [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF] detail a single pixel camera that acquires random projections from the visual scene through a digital micromirror array. A similar acquisition Gabriel Peyré is with the Ceremade, Université Paris-Dauphine, 75775 Paris Cedex 16 France, email: gabriel.peyre@ceremade.dauphine.fr strategy can be used in MRI imaging [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid mr imaging[END_REF] to reduce the acquisition time and increase the spatial resolution.

B. Best Basis Representations

Sparse approximation in orthogonal bases is at the heart of many efficient compression and denoising algorithms. Fixed orthogonal bases are however not flexible enough to capture the complex regularity of sounds or natural images. For instance the orthogonal wavelet transform does not compress efficiently regular edges [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF], [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] and a fixed local cosine basis fails to capture transient parts of musical sounds [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

To improve the sparsity of complicated sounds or images, one can consider several orthogonal bases that compose a large dictionary of atoms. In this framework, one has to choose a best basis adapted to the signal to process. To enable fast computation, this dictionary is required to have a tree structure, so that the best basis can be optimized using a fast dynamic programming algorithm.

Local cosine bases [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] divide the time axis in segments that are adapted to the local frequency content of the sound. Other kinds of dictionaries of 1D bases have been proposed, such as the wavelet packets dictionary [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF] and non stationary wavelet packets [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF], [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF].

The set of cartoon images is a simple model that captures the sketch content of natural images [START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF]. The curvelet frame of Candès and Donoho [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] can deal with such a regularity and enjoys a better approximation rate than traditional isotropic wavelets. This result can be enhanced using a dictionary of locally elongated functions that follow the image geometry. Bandlets bases of Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF], [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF], later refined by Mallat and Peyré [START_REF] Mallat | Orthogonal bandlets bases for geometric image approximation[END_REF], [START_REF] Peyré | Surface compression with geometric bandelets[END_REF], provide such a geometric dictionary together with a fast optimization procedure to compute a basis adapted to a given image.

C. Previous Works

Compressed sensing recovery reconstructs a high resolution signal or image from low dimensional measurements. The resolution of this ill-posed linear inverse problem is regularized by introducing non-linear priors. Initial papers [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF], [START_REF]Decoding by linear programming[END_REF], [START_REF] Donoho | Compressed sensing[END_REF], [START_REF]For most large underdetermined systems of linear equations, the minimal 1 norm solution is also the sparsest solution[END_REF] give theoritical results about the performance of compressed sensing using sparsity in an orthogonal basis. Increasing the efficiency of the priors to fit the features of sounds and natural images might improve the quality of the recovered signal.

Sparsity priors over redundant dictionaries. Increasing the redundancy of the frame used to sparsify the data might be useful to remove reconstruction artifacts. One can for instance use a frame of translation invariant wavelets which has proven useful for image denoising [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Increasing the redundancy can also improve the sparsity of the representation. The curvelet frame [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] better represents edges in images than wavelets. The resulting lack of orthogonality might however deteriorate the coherence of the basis with the sensing vectors, see [START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF].

Advanced fixed priors. Other approaches to enhance compressed sensing reconstruction impose further constraints beyond sparsity. For instance, wavelets coefficients can be optimized in a scale-by-scale fashion [START_REF] Donoho | Extensions of compressed sensing[END_REF] and positivity or a total variation constraint can be enforced [START_REF] Candès | Practical signal recovery from random projections[END_REF].

One can use more advanced signals models not based on sparsity but rather on low dimensional smooth manifolds [START_REF] Baraniuk | Random projections of smooth manifolds[END_REF], [START_REF] Peyré | Manifold models for signals and images[END_REF], union of sub-spaces [START_REF] Eldar | Robust recovery of signals from a union of subspaces[END_REF], or block-sparsity [START_REF] Baraniuk | Model-based compressive sensing[END_REF], [START_REF] Duarte | Fast reconstruction of piecewise smooth signals from random projections[END_REF].

Adaptive priors. A fixed prior does not take into account the time-frequency structures of natural sound or the geometry of edges that is different in each image. Adaptive reconstruction can be obtained using an iterative non-local regularization [START_REF] Peyré | Non-local regularization of inverse problems[END_REF], [START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF]. It is possible to learn an adaptive sparsity prior together with the sensing matrix [START_REF] Carvajalino | Learning to sense sparse signals: Simultaneous sensing matrix and sparsifying dictionary optimization[END_REF], which extends the initial idea of learning the sensing vectors alone [START_REF] Elad | Optimized projections for compressed sensing[END_REF].

In these adaptive methods, an iterative algorithm computes both the recovered signal and parameters that control the regularization prior to match the patterns of the signal. In this paper, we also propose to use an adaptive reconstruction method, but within the setting of sparsity in orthogonal bases. The sparsity is increased by selecting an optimized orthogonal basis in a tree structured dictionary.

D. Contributions and Outline of the Paper

Section II recalls the basics of compressed sensing and 1 reconstruction. Section III details the main contribution of this paper, which the best basis extension of compressed sensing, together with a fast algorithm to perform an approximate minimization. This new framework minimizes an energy on both the signal to recover and on the basis that sparsifies this signal. The exploration of both the set of signals and the set of bases is however not tractable numerically. A fast iterative algorithm is thus derived from a series of surrogate functionals that progressively estimates the best basis. Section IV recalls best basis selection in a tagged tree structured dictionary, that is required to make our best basis compressed sensing fast. Sections V, VI and VII show numerical applications using 1D and 2D dictionaries which show the performance of our scheme on synthetic and natural data.

II. COMPRESSED SENSING

Compressed sensing acquisition and recovery. Compressed sensing acquisition computes a fixed set of n linear measurements of an unknown high resolution signal f ∈ R N with N n

y = Φf = { f, ϕ i } n i=1 ∈ R n with ϕ i ∈ R N .
The price to pay for this compressed sensing strategy is a non-linear reconstruction procedure to recover f from the compressed representation y = Φf .

The sparsity of f in a given orthogonal basis B = {ψ m } N -1 m=0 of R N is measured by using the 0 pseudo norm ||Ψf || 0 , where

Ψf = { f, ψ m } m and ||x|| 0 = # {m \ x m = 0} .
This sparsity can be used to recover a signal that is a solution of the following minimization

min g∈R N ||Ψg|| 0 subject to Φg = y. (1)
The minimization (1) is however combinatorial and thus intractable. It is relaxed by using the

The theory extends to signals that are well approximated with a signal that is s-sparse in B.

Incoherence: the sensing vectors {ϕ i } i should be as different as possible from the sparsity vectors {ψ m } m . This is ensured by monitoring the s-restricted isometry constant δ s , which is the smallest δ > 0 such that

(1 -δ)||x|| 2 ||ΦΨx|| 2 (1 + δ)||x|| 2 . ( 4
)
for any s-sparse signal ||x|| 0 s. One should note that although the sparity assumption (3) constrains the 0 pseudo-norm of Ψf , the actual recovery process (2) optimizes the 1 norm. This is important since the 1 norm is convex, which leads to a tractable optimization problem with fast algorithms, see [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF].

The following recovery theorem ensures the perfect recovery using the 1 minimization. [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF]) If f is s-sparse in Ψ as defined in (3) and if the sensing matrix Φ satisfies δ 2s < √ 2 -1, the solution of (2) satisfies f = f . This theorem ensures that there is an unique solution to [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF] and that it coincides with the solution of (2). It extends to the recovery of an approximately sparse signal, in which case ||f -f || is comparable with the best s-terms approximation of f in the basis B.

Theorem 1. ([1], [2],
Random sensing matrices. The allowable sparsity s for which δ 2s < √ 2 -1 needs to be large enough for Theorem 1 to be useful. For several random matrix ensembles, this is actually the case with high probability on the random matrix. For instance if the entries of Φ are drawn independently from a Gaussian distribution of variance 1/N , s Cn/ log(N/n) ensures δ 2s < √ 2-1 with high probability. This result extends to other kinds of random distributions, see for instance [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF], [START_REF]Decoding by linear programming[END_REF], [START_REF]Compressive sampling[END_REF], [START_REF] Donoho | Compressed sensing[END_REF], [START_REF]For most large underdetermined systems of linear equations, the minimal 1 norm solution is also the sparsest solution[END_REF], [START_REF] Rudelson | On sparse reconstruction from fourier and gaussian measurements[END_REF].

In this paper, following for instance [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], we consider a fast sampling operator Φf = (P 1 HP 2 f ) ↓ n , where P 1 and P 2 are realizations of a random permutation of the N entries of a vector in R N , H is a 1D or 2D orthogonal Hadamard transform, and ↓ n selects the n first entries of a vector, see [START_REF] Shanks | Computation of the fast walsh-fourier transform[END_REF] for a definition of the Hadamard transform and its fast implementation. Such a random sensing operator is computed in O(N log(N )) operations, which is important to process high dimensional data.

Robust Compressed Sensing. To deal with noisy measurements y = Φf + w, where w is a bounded noise ||w|| σ, one can turn the constrained formulation (2) into a penalized variational problem f = argmin g∈R N E(g, B, t) where ( 5)

E(g, B, t) = 1 2 ||Φg -y|| 2 + t||Ψg|| 1 . (6)
where the Lagrange multiplier t should be set so that ||Φf -y|| ≈ σ. Compressed sensing theory extends to this noisy setting, since the recovery error ||f -f || is of the order of the noise level σ, see [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF]For most large underdetermined systems of linear equations, the minimal 1 norm solution is also the sparsest solution[END_REF].

III. BEST BASIS COMPRESSED SENSING

Fixed bases are not efficient enough to sparsify sounds and natural images, and more redundancy is required. This section extends the recovery process (5) to a dictionary of orthogonal bases, the union of which contains a large collection of atoms.

A. Dictionaries and Lagrangian

Dictionaries of Orthogonal bases. A dictionary is a set

D Λ = {B λ } λ∈Λ of orthogonal bases B λ = {ψ λ m } m of R N .
Instead of using an a priori fixed basis such as the wavelet or Fourier basis, one chooses a parameter λ ∈ Λ adapted to the structures of the signal to process and use the optimized basis B λ . This dictionary D Λ defines a highly redundant set of atoms

D Λ = ψ λ m \ λ ∈ Λ, m .
One could use directly D Λ to reconstruct a signal from compressive measurements, but this is numerically intractable for large signals or images. Restricting the sparsity to orthogonal bases leads to fast algorithms for structured dictionaries such as those considered in Section IV.

Tree structure. To enable the fast optimization of a parameter λ adapted to a given signal or image f to process, we impose that each λ ∈ Λ is a tree, that is a subset of the larger tree Λ. Typical examples use the tree λ to index a segmentation of either space (local cosine bases detailed in Section V) or scale (wavelet packets detailed in Section VI). Each node in a tree λ corresponds to a refinement in the segmentation.

Lagrangian and adapted basis. A best basis B λ adapted to a signal or an image f ∈ R N is optimized to obtain the best possible approximation of f for a given sparsity, as measured by the 1 norm. This constrained formulation is turned into an Lagrangian unconstrained optimization problem by minimizing a Lagrangian E

λ = argmin λ∈Λ E(f, B λ , t) (7)
E(f, B λ , t) = min g∈R N 1 2 ||f -g|| 2 + t||Ψ λ g|| 1
and where

Ψ λ g = { g, ψ λ m } N -1 m=0
. The variable t is a Lagrange multiplier that weights the quality of approximation in the chosen basis with the sparsity of the expansion.

The following lemma characterizes the best basis together with the best 1 -penalized approximation in this basis. Its proof can be found for instance in [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Lemma 1. The minimizer (f , λ ) = argmin (g,λ)∈∈R N ×Λ 1 2 ||f -g|| 2 + t||Ψ λ g|| 1 is given by λ = argmin λ∈Λ E(f, B λ , t), f = S t (f, B λ ), (8)
where the soft thresholding operator is defined as

S t (f, B) = m s t ( f, ψ m ) ψ m (9)
where

s t (x) = max(0, 1 -t/|x|)x
The soft thresholding S t defined in (9) computes the best 1 approximation of a signal f . The Lagrangian E(f, B λ , t) can thus be written as a sum over the atoms

{ψ λ m } m of B λ E(f, B λ , t) = m γ t (| f, ψ λ m |) (10)
where

γ t (x) = x 2 /2 if |x| t, t|x| -t 2 /2 otherwise.

B. Best Basis Compressed Sensing Reconstruction

Best basis regularization. Compressed sensing is extended to a dictionary of bases D Λ by imposing that the recovered signal is sparse in at least one basis of D Λ . The original recovery procedure ( 5) is extended to this new setting by performing the minimization both on the signal to recover and on the basis

(f , λ ) = argmin (g,λ)∈R N ×Λ E(g, B λ , t), (11)
where the energy E is defined in equation [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF].

Informal discussion about best basis regularization. For several models of sounds and natural images, best basis methods have been shown either experimentally or theoretically to outperform fixed representations to sparsify the data. Although for the moment no theory is able to assess the performance of an adaptive minimization such as [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF], Sections V, VI and VII show numerically that the increase of sparsity obtained with best basis computations translates to an increase of performances for the regularization of compressed sensing reconstruction.

We believe that the effectiveness of best basis regularization for compressed sensing comes in large part from the randomness of the measurements, which translates, as already noticed by Donoho et al. [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], on reconstruction artifacts being close to Gaussian white noise during an iterative thresholding reconstruction method. This key feature is very different from the kind of artifacts obtained for non-random sensing schemes, such as Radon transform in tomography or low pass filtering in optical imaging system. This results in a high incoherence between the residual noise with the sparsity vectors of D Λ during our reconstruction scheme. This incoherence allows us to efficiently estimate an adapted basis with an iterative scheme without prior knowledge about the signal to reconstruct.

C. Surrogate Functionals

Searching in the whole dictionary D Λ for the best basis parameter λ ∈ Λ that minimizes [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF] is not feasible for large dictionaries, which typically contain of the order of 2 N bases. Furthermore, the under-determinancy of Φ creates coupling in the non-linear set of equations involved in the minimization of E.

To solve this issue, we relax the energy minimization (11) and define an approximate energy that is simpler to minimize. If one has some estimate h of the solution f , the energy E can be replaced by the following surrogate functional that depends on h

E h (g, B λ , t) = E(g, B λ , t) + µ 2 ||g -h|| 2 - 1 2 ||Φg -Φh|| 2 .
Such a surrogate functional is introduced by Daubechies et al. [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] to derive an iterative thresholding algorithm that solves the 1 regularization of inverse problems.

As long as µ > ||Φ * Φ||, the surrogate energy E h is a smooth, strictly convex, modification of the original energy E. Furthermore,

E h E and E h (h, B λ , t) = E(h, B λ , t).
It is thus reasonable to use E h as a proxy for the original minimization [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF], which leads us to consider, for a fixed h,

(f (h), λ (h)) = argmin (g,λ)∈R N ×Λ E h (g, B λ , t). ( 12
)
The following proposition shows that this surrogate minimization diminishes the energy E.

Proposition 1. If µ > ||Φ * Φ||, one has for all λ E(f (h), B λ (h) , t) E(h, B λ , t). Proof: Condition µ > ||Φ * Φ|| ensures that µ 2 ||g -h|| 2 - 1 2 ||Φg -Φh|| 2 0 so that E h (g, B λ , t) E(g, B λ , t).
One thus has

E(f (h), B λ (h) , t) E h (f (h), B λ (h) , t) E h (h, B λ , t) = E(h, B λ , t).
The following theorem ensures that the minimization ( 12) is easily solved using a best basis search and a soft thresholding.

Theorem 2. The minimization (12) has a global minimum which is given by

λ (h) = argmin λ∈Λ E( h, B λ , t), f (h) = S t/µ ( h, B λ (h) ), (13)
where h = h + 1 µ Φ T (y -Φh),
where the Lagrangian E is defined in equation ( 7) and the soft thresholding operator S t/µ is defined in equation [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Proof:

The energy E h is expanded as follow

E h (g, B λ , t) = µ 2 ||g|| 2 -Φg, y + Φg, Φh -µ g, h + t||Ψ λ g|| 1 + C h,y = µ 2 ||g|| 2 -g, µh + Φ T (y -Φh) + t||Ψ λ g|| 1 + C h,y
where C h,y is independent of λ and g. Up to multiplicative and additive constants, one has

E h (g, B λ , t) ∝ 1 2 ||h + 1 µ Φ T (y -Φh) -g|| 2 + t||Ψ λ g|| 1 .
The result of the theorem follows from lemma 1. Equation [START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF] shows that the parameter λ (h) is the best basis parameter of the modified guess h. It also shows that f (h) is obtained by thresholding h in that basis. The best parameter λ (h) is thus found by optimizing a Lagrangian E, which can be achieved with a fast algorithm, as explained in Section IV.

D. Best Basis Compressed Sensing Algorithm

The minimization of E(g, B λ , t) over (g, λ) is replaced by the minimization of a set of surrogate functionals E f k (g, B λ , t) where f k is the current estimate of the solution at iteration k.

Starting from an initial pair (f 0 , λ 0 ), this defines

(f k+1 , λ k+1 ) = argmin (g,λ)∈R N ×Λ E f k (g, B λ , t). ( 14
)
Algorithm 1 details the step by step implementation of ( 14), that is derived from ( 14) using Theorem 2. The steps of the algorithm are repeated until a user defined tolerance η is reached. In the noiseless setting y = Φf , one can still use this algorithm and decay the value of t = t k toward zero during the iterations, in a fashion similar to the Morphological Components Algorithm (MCA) algorithm [START_REF] Starck | Redundant multiscale transforms and their application for morphological component analysis[END_REF]. This algorithm extends to a best basis setting the iterative thresholding algorithm that has been proposed by several authors to solve inverse problems such as [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF], see for instance [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] and the references therein.

Convergence of the the algorithm. For a single basis B, condition µ > ||Φ * Φ|| ensures that this iterative thresholding converges to a solution of [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF], see [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. This condition can in fact be relaxed to µ > ||Φ * Φ||/2, see [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]. In contrast to the Algorithm 1: Best-basis compressed sensing algorithm.

Initialization.: Set k = 0, f 0 = 0. repeat Update the estimate: Set

fk = f k + 1 µ Φ T (y -Φ f k ). ( 15
)
Update best basis: Compute the best basis parameter

λ k+1 = argmin λ∈Λ E( fk , B λ , t). ( 16
)
This minimization is carried out with a fast procedure, as detailed in section IV.

Denoise the estimate: Compute

f k+1 = S t/µ ( fk , B λ k+1 ),
where S t is the threshold operator defined in [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

until ||f k+1 -f k || < η ;
fixed basis setting, the iterations of ( 14) that require a best basis search are difficult to analyze. Proposition 1 shows that the energy E(f k , B λ k , t) is decaying through the iterations and thus converges to some limit value. Since this energy is nonconvex in (f, λ), nothing more can however be said about the convergence of the iterates f k .

In theory, nothing prevents the best basis parameter λ k to cycle through a series of limit bases, but this was never the case in the numerical experiments performed in Sections V-B, VI-B and VII-B. We always observed that λ k converged to its final value λ after a small number of iterations, that is of the order of 10 to 20 for natural sounds and images. Once this final regime is attained, one can use classical results about the convergence of iterative thresholding algorithms, that show for instance that the energy E(f k , B λ , t) is decaying at speed O(1/k) through the iterations, see [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF].

We note that since the objective to minimize is non-convex and there is no convergence guarantee, a different initialization f 0 might lead to a different sequence of bases B λ k , and hence a different final result. Although this is indeed the case, we did not observed in our numerical tests significant improvement of using alternate initialization strategies with respect to setting f 0 = 0.

IV. BEST BASIS COMPUTATION

This section reviews a best basis framework common to the dictionaries used in the numerical experiments of this paper. These dictionaries enjoy a hierarchical tree structure, which makes the computation of the best basis in [START_REF] Peyré | Surface compression with geometric bandelets[END_REF] fast. This framework does not cover the special case of non-stationary wavelet packets, detailed in Section VI that requires a more advanced procedure, see [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF].

A. Tree structured dictionaries

Structure of the atoms. This paper focusses on dictionaries D Λ having a multiscale tree structure. The atoms of D Λ span sub-spaces V j,i of R N for scales 0 j J = log 2 (N )/d and position 0 i < 2 dj , that obey a refinement relationship

∀ j < J, V j,i = 2 d -1 ε=0 V j+1,2 d i+ε ,
the sum being orthogonal. Each sub-space V j,i has dimension N/2 dj and is equipped with one or several orthogonal bases indexed by a token ∈ Ω

B j,i = ψ j,i,s \ ∀ 0 s < N/2 dj . ( 17
)
Structure of an orthogonal basis. The parameter λ that indexes a basis B λ ∈ D Λ is a binary tree (for 1D signals, where d = 1) or a quad-tree (for 2D, where d = 2). The set of nodes of λ is denoted as N (λ) and each node (j, i) ∈ N (λ) is located at some level 0 j J and position 0 i < 2 dj . A node (j, i) is thus located in the j th row and the i th column of the tree. Each interior node (j, i) ∈ I(λ) ⊂ N (λ) has 2 d children {(j + 1, 2 d i), . . . , (j + 1, 2 d (i + 1) -1)}. The leaves nodes (j, i) ∈ L(λ) have no child.

A basis B λ is obtained by aggregating bases B j,i defined in (17) for (j, i) that are leaves of λ, and for a specific choice = j,i ∈ Ω of token at each node of the tree B λ = ψ j,i j,i,s \ ∀ (j, i) ∈ L(λ) and 0 s < N/2 dj . Figure 1 shows an example of such a tagged quad-tree λ. When one does not care about the location of the basis elements in the tree, the basis is written as B λ = {ψ λ m } m where the index is m = (j, i, s, j,i ) with (j, i) ∈ L(λ) and 0 s < N/2 dj .

Examples of dictionaries for sounds and geometrical images. Sections V, VI and VII detail three kinds of dictionaries adapted to various signals and images structures.

The local cosine dictionary [10], [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] is used to process highly oscillating signals (d = 1) such as music and speech audio data. The binary tree λ segments the time axis in order to match the variations of the local frequencies in the sound. In this case |Ω| = 1 since there is no need for additional information beside the spatial segmentation.

The non-stationary wavelet packets dictionary [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF], [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF] is used to process signals (d = 1) that require an arbitrary tiling of the scale axis. The non-stationary cascade of filterings also allows to adapt the basis functions through the scales. In this case, j,i ∈ Ω = {0, . . . , κ -1} indicates the index of a wavelet filter that is used to subdivide the frequency axis. The particular case of wavelet packets [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF] is obtained for κ = 1.

The bandlet dictionary [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF], [START_REF] Mallat | Orthogonal bandlets bases for geometric image approximation[END_REF] is used to process images (d = 2) with geometric features such as edges or directional textures. The quad-tree λ segments the square [0, 1] 2 into sub-squares S j,i ⊂ [0, 1] 2 of size 2 -j × 2 -j . A leaf node (j, i) ∈ L(λ) defines the local bandlet transform and caries a geometrical information token j,i ∈ Ω = {0, . . . , κ-1}∪∅.

A token j,i = ∅ is defined over a geometrical square S j,i in which the bandlet vectors are elongated and follow the angle π j,i /κ that approximates the orientation of the closest edge. A token j,i = ∅ is defined in isotropic square S j,i that corresponds to regions where no edge is present. Other classes of dictionaries include bases composed of atoms with rapidly varying oscillations such as the modulated bases of Coifman et al. [START_REF] Coifman | Modulated Malvar-Wilson bases[END_REF] and the chirplets dictionary of Candès [START_REF] Candès | Multiscale chirplets and near-optimal recovery of chirps[END_REF].

B. Dynamic programming for Best-basis Computation

A fast best basis search algorithm makes use of the tree structure of D Λ and the fact that the Lagrangian E(f, B λ , t) to minimize is split as a sum over the coefficients of the decomposition [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF]. It was originally presented by Coifman et al. [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF] and is a particular instance of the Classification and Regression Tree (CART) algorithm of Breidman et al. [START_REF] Breiman | Classification and Regression Trees[END_REF] as explained by Donoho [START_REF] Donoho | Cart and best-ortho-basis: A connection[END_REF].

Algorithm 2 details the implementation of this fast best basis search. This algorithm does not work for non-stationary wavelet packets, which require a more complex optimization procedure described in [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF].

This algorithm requires the decomposition of f onto each atom ψ j,i,s for all values of (i, j, s, ), and the computation of Lagrangians E j,i that depends on these inner products. For several multiscale dictionaries, such as those considered in this paper, fast algorithms perform this computation O(P ) operations, where P is the total number of atoms in D Λ . The resulting complexity is thus O(|Ω|N log 2 (N )), where |Ω| = 1 for the local cosine and stationary wavelet packets, and |Ω| is the number of orientations in the bandlet dictionary.

C. Settings for the Numerical Results

The following sections detail several dictionaries of orthogonal bases. The performance of these dictionaries is illustrated for compressed sensing recovery using the same numerical experiments. The recovery success is measured using

PSNR(f, f ) = -20 log 10 (||f -f ||),
where the signal f is assumed to take values in [0, 1]. This recovery error is measured for various values of the sensing rate n/N ∈ [0, 1]. Three kinds of recoveries are compared:

Recovery using a fixed basis (for instance fixed local DCT or fixed orthogonal basis) using the original optimization (2).

Algorithm 2: Best-basis selection algorithm.

for 0 j J, 0 i < 2 jd -1, ∈ Ω do Compute f, ψ j,i,s . for 0 j J, 0 i < 2 jd -1 do Compute j,i = argmin ∈Ω s γ t (| f, ψ j,i,s |).
for j = J, . . . , 0 do for i = 0, . . . , 2 jd -1 do Compute

E j,i = s γ t (| f, ψ j,i j,i,s |), Ẽj,i = 2 d -1 ε=0 Ẽj+1,2 d i+ε ,
with the convention that Ẽj,i = +∞ if j = J. if Ẽj,i E j,i then Declare (j, i) as interior (j, i) ∈ I(λ).

if Ẽj,i > E j,i then Declare (j, i) as leaf (j, i) ∈ L(λ). Update the cumulative Lagrangian Ẽj,i ← min( Ẽj,i , E j,i ).

Recovery using the oracle best-basis B λ (f ) estimated from the original signal f . This is an upper-bound for the performance of our method since the knowledge of this basis is not available in practice. Recovery using the algorithm of section III-D, that estimates iteratively the best basis. Experiments are performed on noisy measurements y = Φf 0 + w where w is a Gaussian white noise of variance σ 2 . The noise level σ is supposed to be known, and the regularization parameter t is adjusted so that the residual final error satisfies ||y -Φf || = √ nσ. For noiseless experiments where σ = 0, we still use our iterative thresholding algorithms, but decay linearly the threshold t = t k to 0 during the iterations, which is similar to the approach used in the Morphological Component Algorithm of Starck et al. [START_REF] Starck | Redundant multiscale transforms and their application for morphological component analysis[END_REF]. We used a relaxation parameter µ = 1 for the iteration (15).

V. BEST LOCAL COSINE BASIS COMPRESSED SENSING

A. Adapted Local Cosine Transform

A local cosine basis B λ is parameterized with a binary tree λ that segments the time axis in dyadic intervals, see [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF], [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Each leaf node (j, i) ∈ L(λ) corresponds to a selected interval [x j,i , x j,i+1 ], where x j,i = 2 -j N i -1/2. For each of these leave node, the local cosine basis vectors are defined as

∀ (j, i) ∈ L(λ), ∀ k, ψ j,i,s [k] = b 2 j (k -x j,i ) 2 2 -j N cos π s + 1 2 k -x j,i 2 -j N ,
where b is a smooth windowing function that satisfies some compatibility conditions [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

The decomposition of a signal f on the vectors of some basis B λ is computed in O(N log(N )) using fast Fourier transforms. A best basis B λ that minimizes [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] is computed in O(N log(N ) 2 ) operations using Algorithm 2. Figure 2 shows some examples of basis vectors.

j = 0 j = 1 j = 2 i = 0 i = 0 1 0 1 2 3
Fig. 2. A dyadic tree λ defining a spatial segmentation (left) ; some local cosine basis functions ψ j,i,s of the basis B λ (right).

B. Numerical Results

A synthetic sparse signal f = (Ψ λ ) -1 h is generated using a random local cosine basis B λ and a random signal of spikes h with ||h|| 0 = # {k \ h(k) = 0} = 30, see figure 3,(a). We use here noiseless measurements y = Φf . The signal recovered by the non-adaptative algorithm of section III-D in an uniform cosine basis B λ0 is significantly different from the original, figure 3,(b). This is due to the fact that f is less sparse in B λ0 , since

||Ψ λ0 f || 0 = 512 and ||Ψ λ0 f || 1 ≈ 2.8||Ψ λ f || 1 .
During the iterations of the algorithm presented in subsection III-C, the estimated best basis B λ k evolves in order to match the best basis B λ , see figure 3, (c1-c3). The recovered signal f in (c3) is nearly identical to f . We note however that the segmentation λ optimized by our algorithm differs slightly from the original segmentation λ 0 . In particular one can notice an over-segmentation of the leftmost interval, where the recovered signal f is small but not vanishing, on the contrary to the original signal f . On figure 4 one can see a sound of a tiger howling, together with the signals recovered using a fixed fully subdivided local DCT basis and the best basis recovery algorithm of section III-C. We use here noisy measurements y = Φf + w, with ||w|| = 0.03||Φf ||. Although the final adapted basis is not the same as the best basis of the original signal, it still provides an improvement of 2dB with respect to a fixed spatial subdivision. Figure 5 shows for various rates of sensing the recovery error, confirming that the iterative algorithm does not perform as good as the oracle best best basis computed from f .

VI. BEST NON-STATIONARY WAVELET PACKET COMPRESSED SENSING A. Adapted Non-stationary Wavelet Packet Transform

The non stationary (NS) wavelet transform and its extension to wavelet packets was introduced by Cohen et al. [START_REF] Cohen | Nonstationary subdivision schemes and multiresolution analysis[END_REF]. We give the definition of a non-stationary wavelet transform which corresponds to the decomposition in an orthogonal basis parameterized by a tagged binary tree λ. This tree structured dictionary is more general than the trees considered in Section IV-A because a token j,i is associated to each node (j, i) ∈ I(λ) ∪ L(λ) and not only to leaves (j, i) ∈ L(λ). The corresponding best basis algorithm is thus a generalization of Algorithm 2that is detailed in [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF].

The NS wavelet packet dictionary depends on the choice of a set H = {h 0 , . . . , h κ-1 } of low pass quadrature mirror filters, to which is associated the high pass filters g

[k] = (-1) k h [1 -k]. NS Wavelet Packet Transform. The decomposition of f = f 0,0 ∈ R N in a NS wavelet packet basis B λ iteratively computes f j,i ∈ R N/2 j
by traversing the nodes (j, i) of the tree λ from top to bottom and computing

f j+1,2i = (f j,i * h j,i ) ↓ 2, f j+1,2i+1 = (f j,i * g j,i ) ↓ 2, (18)
where ↓ 2 is the sub-sampling by two operator. The coefficients computed on the leaves correspond to the projection of f on the vectors of B λ ∀ (j, i) ∈ L(λ), f j,i [s] = f, ψ λ j,i,s . This forward NS wavelet packet transform is computed in at most O(N log 2 (N )) operations for a signal f ∈ R N . A best basis B λ defined by [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] adapted to some signal f is computed using a fast best-basis search introduced by Peyré and Ouarti [START_REF] Ouarti | Best basis denoising with non-stationary wavelet packets[END_REF]. The complexity algorithm is O(N log 2 (N )) for κ = 1 (which corresponds to the traditional wavelet packets) and O(N 1+log 2 (κ) ) for κ > 1.

B. Numerical Results

In our tests, the non-stationary wavelet packet dictionary is built using a family H D = {h } κ-1 =0 where κ = 5 and where h is the Daubechies orthogonal filter with + 1 vanishing moments.

Figure 6 (c,d,e) shows a comparison of the recovery using the Daubechies wavelets corresponding to h 3 , the wavelet packets dictionary which corresponds to H = {h 3 } and the non-stationary wavelet packets dictionary which corresponds to H = H D . The signal f (x) = f 0 (x) + ε sin(ωx) is the superposition of a piecewise-regular signal and a sinusoid with high frequency ω. We use here noiseless measurements y = Φf . Figure 6 (a) shows the index j,i of the best NS wavelet packets basis, which is able to capture both the high frequency content of sin(ωx) while minimizing the number of large coefficients created by the singularities of f 0 . (c) original signal, N = 4096 ; (d) recovery using orthogonal wavelets, n = N/3 (PSNR=23.77dB) ; (e) recovery using best wavelet packets basis, n = N/3 (PSNR=26.50dB), the corresponding dictionary tree is shown in (a) ; (f) recovery using best non-stationary wavelet packets basis, n = N/3 (PSNR=29.42dB), the corresponding dictionary tree is shown in (b). Figure 7 shows another example of recovery using wavelets, the wavelet packets and the NS wavelet packets dictionaries.

The signal f has dimension N = 4096 and is the superposition of chirps and spikes

f (x) = 3 k=1 cos(ω k (x/N ) k ) + 10 k=0 µ k (1 + |(x/N ) -x k |/ρ k ) -4
for 0

x < N , where ω 1 = π/2N, ω 2 = πN and ω 3 = 0.7πN , where the x k ∈ [0, 1] are random position, µ k ∈ [START_REF]Decoding by linear programming[END_REF][START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF] are random weights and ρ k ∈ [0, 0.03] are random widths. We use here noisy measurements y = Φf + w with ||w|| = 0.03||Φf ||. The best NS wavelet packets basis is able to better recover both the high frequency oscillations of the chirps and the localized spikes than fixed representations.

VII. BEST BANDLET BASIS COMPRESSED SENSING

A. Adapted Bandlet Transform

The bandelet bases dictionary was introduced by Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF], [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF]. Bandlets perform an efficient adaptive approximation of images with geometric singularities, such as the cartoon image in Figure 10. This transform has been refined by Mallat and Peyré [START_REF] Mallat | Orthogonal bandlets bases for geometric image approximation[END_REF], [START_REF] Peyré | Surface compression with geometric bandelets[END_REF] to obtain a dictionary of regular and orthogonal basis functions. We present a simplified bandlet transform inspired from [START_REF] Peyré | Surface compression with geometric bandelets[END_REF]. This implementation results in a decomposition similar to the wedgelets of Donoho [START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF] but within the framework of a dictionary of orthogonal bases.

The quadtree λ that parametrizes the bandelet basis B λ defines a segmentation

[0, 1] 2 = (j,i)∈L(λ) S j,i
where each S j,i is a square of size 2 -j × 2 -j . Figure 8 shows an example of such a dyadic subdivision.

The sub-image fj,i = {f (x)} x∈Sj,i ∈ R N/2 2j extracted from the square S j,i is retransformed according to the value of the token j,i . If j,i = ∅, it means that no edge is present in S j,i and an orthogonal 2D wavelet transform is applied to fj,i to obtain f j,i . If j,i = ∅, an edge passes through S j,i with a direction close to the angle θ j,i = π j,i /κ ∈ [0, π). Then a 1D directional wavelet transform is applied to fj,i to obtain f j,i . This transform is computed by ordering the Left: example of dyadic subdivision of [0, 1] 2 in squares S j,i ; right: corresponding quad-tree λ. pixels of fj,i along the direction θ j,i and then applying a 1D Haar transform, as described in [START_REF] Peyré | Surface compression with geometric bandelets[END_REF]. Figure 9 (c,d) shows this directional wavelet transform. The resulting coefficients f j,i correspond to projections of f on the bandlet vectors of the basis B λ ∀ (j, i) ∈ L(λ), f j,i [s] = f, ψ j,i j,i,s . Keeping only a few bandelet coefficients and setting the others to zero performs an approximation of the original image that follows the local direction θ j,i , see figure 9 (f).

A best basis B λ defined by ( 7) adapted to some function f is computed using a fast best-basis search introduced by Mallat and Peyré [START_REF] Mallat | Orthogonal bandlets bases for geometric image approximation[END_REF]. Figure 10 shows an example of such a segmentation adapted to a geometric image. Example of subdivision λ adapted to a geometrically regular function f . The basis B λ minimizes E(f, B λ , t).

B. Numerical Results

The geometric image f depicted in figure 12 (a), (a) is used to compare the performance of the original compressive sensing algorithm in a wavelet basis to the adaptative algorithm in a best bandlet basis. We use here noisy measurements y = Φf + w with ||w|| = 0.03||Φf ||. We use a translation invariant 7/9 wavelet tight frame, which is more efficient for inverse problems than orthogonal wavelets. Since the wavelet basis is not adapted to the geometric singularities of such an image, reconstruction (b) has ringing artifacts. The adapted reconstruction (c) exhibits fewer such artifacts since the bandlet basis functions are elongated and follow the geometry. The segmentation is depicted after the last iteration, together with the chosen directions θ j,i that closely match the direction of the edges of f . Figure 12,(a',b',c') shows the recovery results for a natural image containing complex geometric structures such as edges, junctions and sharp line features. The best bandlet recovery is able to resolve these features efficiently. Figure 11 shows the recovery error for various sensing rates, in the noiseless setting y = Φf . For low rate (n/N close to 0), the basis estimated by the iterative algorithm is not as good as the oracle basis estimated from f . For higher rates, the algorithm is able to find the underlying geometry efficiently and the algorithm performs as good as the oracle basis.

VIII. CONCLUSION

This paper tackles the problem of improving the reconstruction from compressed sensing measurements. The proposed method extends compressed sensing recovery to the setting where the sparsity is imposed in a tree structured dictionary of orthogonal bases. This best basis adaptive prior enhances the sparsity of the signal to recover with respect to fixed priors in an orthogonal basis. An iterative algorithm minimizes a variational energy over both the signal to reconstruct and the best basis parameter. Numerical results show that this approach is successful for natural sounds and geometric images that contain a broad range of sharp transitions. Best basis regularization extends to more general linear inverse problems, and could potentially improves the reconstruction beyond randomized measurements.

Fig. 1 .

 1 Fig. 1. Examples of quad-tree λ.

Fig. 3 .

 3 Fig. 3. (a) synthetic sound signal with 30 random cosine atoms N = 4096 ; (b) recovery using a fixed cosine basis ; (c1) first iteration of the best basis recovery algorithm, n = N/3 ; (c2) iteration s = 5 ; (c3) iteration s = 20.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) sound signal of a tiger howling, together with the best spatial segmentation, N = 32768 ; (b) recovery using fixed local cosine basis, n = N/3 (PSNR=19.24dB) ; (c) recovery using best cosine basis, n = N/3 (PSNR=21.32dB)

  Fig. 6.(c) original signal, N = 4096 ; (d) recovery using orthogonal wavelets, n = N/3 (PSNR=23.77dB) ; (e) recovery using best wavelet packets basis, n = N/3 (PSNR=26.50dB), the corresponding dictionary tree is shown in (a) ; (f) recovery using best non-stationary wavelet packets basis, n = N/3 (PSNR=29.42dB), the corresponding dictionary tree is shown in (b).

Fig. 7 .

 7 Fig. 7. (a) original signal, N = 4096 ; (b) recovery using orthogonal wavelets, n = N/3 (PSNR=17.59dB) ; (c) recovery using best wavelet packets basis, n = N/3 (PSNR=20.46dB) ; (d) recovery using best nonstationary wavelet packets basis, n = N/3 (PSNR=21.66dB).

Fig. 8 .

 8 Fig. 8.Left: example of dyadic subdivision of [0, 1] 2 in squares S j,i ; right: corresponding quad-tree λ.

Fig. 10 .

 10 Fig. 10.Example of subdivision λ adapted to a geometrically regular function f . The basis B λ minimizes E(f, B λ , t).

Fig. 11 .

 11 Fig. 11. Recovery results for various rate of sensing n/N , in the noiseless setting y = Φf .

Fig. 9 .

 9 Fig. 9. (a) a geometric image together with an adapted dyadic segmentation ; (b) a square S j,i together with its adapted direction θ j,i ; (c) the 1D signal obtained by re-ordering the pixels of fj,i along direction θ j,i ; (d) the 1D Haar coefficients f j,i of fj,i ; (e) the 1D approximation obtained by reconstruction from the 20 largest Haar coefficients ; (f) the corresponding square approximated in bandelet.

Fig. 12 .

 12 Fig. 12. (a,a') original image ; (b) compressed sensing reconstruction using the translation invariant wavelet frame, n = N/6 (PSNR=37.1dB) ; (c) reconstruction using iteration in a best bandlet basis (PSNR=39.3dB). (b') wavelet frame, PSNR=22.1dB, (c') bandlet basis, PSNR=23.9dB.