Iterative Hard Thresholding IOP Instance Optimality Property KL Kullback-Leibler KLSH Kernelized Locality Sensitive Hashing KPCA Kernel Principal Component Analysis LSH Locality Sensitive Hashing NSP Null Space Property PCA Principal Component Analysis RIP Restricted Isometry Property RMMH Random Maximum Margin Hashing SVM Support Vector Machine

Mes remerciements se dirigent en premier lieu vers mon directeur de thèse Rémi Gribonval et mon encadrant en entreprise Patrick Pérez. Je ne saurais exprimer pleinement ma gratitude pour leur encadrement de rare qualité. Leur pédagogie, leur acuité scientifique et leur passion communicative n'ont cessé de m'impressionner, de m'inspirer et de me donner goût à la recherche pendant ces trois années. J'ai particulièrement apprécié la capacité de Rémi à dégager avec une étonnante rapidité le sens profond des concepts alors que j'étais aveuglé par certains aspects techniques, ainsi que sa grande patience face à mes questions stupides, vagues ou mal formulées. Patrick m'a maintes fois passionné par sa culture gargantuesque des sciences numériques, et son sens du détail m'a été extrêmement profitable. Je leur adresse donc mes plus sincères remerciements pour ces trois années au cours desquelles ils m'ont tant appris, aussi bien au niveau technique que méthodologique.

Merci également à Hervé Jégou, avec qui j'ai eu le plaisir de travailler et qui m'a fait découvrir et comprendre a plusieurs techniques. Sa vision très pragmatique des concepts m'a inspiré. J'ai également une pensée amicale pour mes autres coauteurs Mike Davies, Tomer Peleg et Florent Perronnin.

Je remercie bien sûr les rapporteurs Francis Bach et Gilles Blanchard, et plus généralement l'ensemble du jury, pour l'intérêt qu'il porte à ces travaux et le temps consacré à l'écriture des rapports et à la soutenance.

J'ai enfin une pensée chaleureuse pour mes collègues de Technicolor et de l'IRISA b pour leur bonne humeur, leurs répliques hilarantes et les discussions de très haute volée c autour d'un café d .

a Enfin, je crois... b Le monde rennais étant petit, certains ont d'ailleurs endossé les deux casquettes successivement. c Ou pas ? d Ou autre boisson chaude pour les quelques braves qui ne succombent pas à l'appel de la caféine.

iii iv

REMERCIEMENTS

Résumé des contributions

Cette section rédigée en français reprend la motivation initiale de la thèse et en résume les contributions principales. Le manuscrit principal commence à la suite de ce résumé et est rédigé en anglais.

En parallèle, des réductions similaires de dimension tout en gardant l'"information utile" ont été exploitées dans plusieurs outils de traitement de données à grande échelle qui opèrent sur les données en réduisant leur taille de façon analogue à ce qui est fait en échantillonnage compressé. Il semble ainsi que l'échantillonnage compressé ait un réel potentiel en apprentissage. Le point de départ de la thèse est donc l'exploration de ce potentiel, en examinant si des méthodes d'échantillonnage compressé peuvent s'appliquer avec succès à des problèmes d'apprentissage statistique.

La thèse comprend trois contributions principales, faisant chacune l'objet d'un chapitre du présent manuscrit.

Échantillonnage compressé pour l'estimation de densité

La première contribution est un apport immédiatement relié à la question évoquée précédemment : des méthodes d'échantillonnage compressé peuvent-elles s'appliquer dans un cadre d'apprentissage statistique ? En échantillonnage compressé, on considère généralement un vecteur x s'écrivant comme combinaison linéaire d'un nombre réduit de vecteurs pris dans une famille {u j } (la famille étant généralement orthogonale ou composée de vecteurs quasiment orthogonaux). On cherche alors à reconstruire x à partir de mesures linéaires Mx, où M est un opérateur linéaire réduisant la dimension.

Or, des modèles analogues de "signaux" apparaissent en apprentissage statistique : un problème classique d'apprentissage non supervisé est par exemple l'estimation d'un mélange de densités à partir d'un ensemble de vecteurs. Dans ce cadre, on peut considérer que les vecteurs sont indépendemment et identiquement distribués selon une loi de probabilité de densité p qui s'écrit comme combinaison linéaire d'un nombre réduit de densités prises dans une famille {p θ }. Le but de la tâche d'estimation est alors de trouver lesquelles de ces densités apparaissent dans la décomposition de p et d'estimer les coefficients de cette décomposition.

La ressemblance entre ces deux cadres de travail suggère que l'on pourrait considérer le problème d'estimation de paramètres de mélange comme un problème d'échantillonnage compressé : il suffit en effet de construire un opérateur de mesure M qui, appliqué à une densité de probabilité p, en calcule une représentation compressée Mp. Le problème d'estimation des paramètres de mélange pourrait alors être interprété comme la reconstruction de la densité p comme combinaison linéaire des {p θ }.

Dans notre première contribution (Bourrier et al., 2013b,c) :

• Nous proposons un cadre de travail analogue à celui de l'échantillonnage compressé pour l'estimation de paramètres de mélange.

• Nous identifions des opérateurs M qui permettent de calculer une estimation de Mp à partir d'un ensemble de vecteurs tirés selon p.

• Nous instancions ce cadre de travail sur un exemple simple, où les densités considérées sont des Gaussiennes isotropes et l'opérateur M est un opérateur d'échantillonnage de la transformée de Fourier. Nous prouvons qu'un choix judicieux de fraquences assure l'injectivité de l'opérateur de mesure associé sur l'ensemble des mélanges parcimonieux de Gaussiennes.

• Nous proposons un algorithme inspiré d'une méthode utilisée en échantillonnage compressé permettant d'estimer les paramètres de mélange à partir de la représentation compressée estimée sur les données.

• Nous proposons des résultats numériques d'expériences effectuées dans le cadre cité ci-dessus où les densités sont des Gaussiennes isotropes, et nous mesurons la qualité de l'estimation par rapport à un algorithme classique.

L'estimation de paramètres de mélange à l'aide de Gaussiennes isotropes peut être interprété comme un problème classique d'apprentissage non supervisé : le partitionnement de données. Nos expériences montrent, sur un exemple simple, qu'il est possible d'estimer des paramètres de mélange de façon précise par rapport à un algorithme d'estimation classique, tout en travaillant à partir d'une représentation compressée, et donc moins volumineuse que l'ensemble des données de départ.

Cette analogie entre échantillonnage compressé et estimation de paramètres de mélange suggère la mise au point d'outils théoriques permettant d'étudier les problèmes inverses de reconstruction et d'estimation dans un cadre plus général que celui de l'échantillonnage compressé classique. Il s'agit du point de départ de notre deuxième contribution.

Performances théoriques des décodeurs en grande dimension

En échantillonnage compressé, et plus généralement dans le cadre de problèmes inverses, on cherche à reconstruire un signal x à partir de mesures Mx, où M est un opérateur linéaire donné. En toute généralité, et particulièrement si M réduit la dimension, cette "inversion" ne va être possible qu'à l'aide d'un a priori sur le signal x. L'hypothèse généralement faite est que x est approximativement parcimonieux, c'est-à-dire proche (au sens d'une certaine distance) d'un ensemble Σ k de signaux parcimonieux (typiquement les signaux qui n'ont au maximum que k composantes non nulles dans la base canonique).

L'étape de reconstruction du signal est effectuée par un "décodeur" ∆ qui va servir de pseudo-inverse à l'opérateur M relativement au modèle Σ k : on veut décoder précisément les vecteurs x situés au voisinage de Σ k , sans accorder d'importance à la précision du décodage si x est éloigné de Σ k . Dans ce cadre, une propriété usuelle requise pour un décodeur est l'instance optimality : on veut majorer, pour tout signal x, l'erreur de décodage x -∆(Mx) par un terme de la forme d(x, Σ k ) représentant la distance du vecteur au modèle. Ainsi, le décodage sera d'autant plus précis que le signal est proche du modèle. Des résultats étudiant l'existence de tels décodeurs ont été proposés pour les modèles de vecteurs parcimonieux et, plus récemment, pour les unions finies de sous-espaces.

Dans notre deuxième contribution (Bourrier et al., 2013a) :

• Nous proposons un cadre général relaxant plusieurs hypothèses utilisées dans les travaux précédents sur l'instance optimality. Nous généralisons dans ce cadre les notions de ces travaux, ainsi que leurs relations entre elles. Les généralisations effectuées couvrent notamment le cadre où Σ est un modèle quelconque de signaux et où les mesures sont éventuellement bruitées.

• Nous généralisons un résultat classique d'instance optimality avec normes

• Nous proposons une généralisation de la Propriété d'Isométrie Restreinte (PIR), usuellement utilisée pour majorer l'erreur de reconstruction d'un vecteur en échantillonnage compressé, à des modèles généraux. Nous relions cette PIR à l'existence d'un décodeur instance optimal au sens d'une certaine norme que nous définissons. Nous proposons une majoration de cette norme sous une hypothèse de PIR. Cette majoration est étudiée dans deux cas classiques, pour lesquels nous prouvons qu'elle est équivalente à des normes usuelles.

Ce deuxième groupe de contributions fournit ainsi des outils permettant d'étudier la précision théorique que l'on peut attendre en terme d'instance optimality sur un problème donné. Le cadre développé englobe notamment le cas où ce que l'on cherche à reconstruire n'est pas forcément le signal de départ mais une caractéristique du signal Ax, où A est un opérateur linéaire. L'étude du cas où ce A est non trivial (c'est-à-dire différent de l'identité) soulève de nombreuses questions ouvertes, englobant l'apprentissage sur données compressées.

L'apprentissage non supervisé sur données compressées a déjà été évoqué dans la première contribution, dans laquelle était étudiée une méthode d'estimation de paramètres de mélange sur une base de données de vecteurs compressée globalement. Un autre moyen, plus usuel, de compresser des vecteurs est de réduire individuellement leur dimension, tout en conservant leur nombre. La dernière partie du travail traite de telles méthodes dans le cadre de la recherche de plus proche voisin, problème courant en analyse statistique.

Plongements explicites pour la recherche de plus proche voisin

En recherche de plus proche voisin, on considère un espace muni d'une distance d et un ensemble X = {x 1 , . . . , x L } d'éléments de cet espace. Étant donné un élément-requête q, la recherche de plus proche voisin de q consiste à trouver argmin x∈X d(x, q). Si L est grand et/ou si d est une distance dont le calcul est coûteux, le recours à une recherche approximative est obligatoire. De tels schémas approximatifs de recherche sont nombreux lorsque d est la distance Euclidienne usuelle, mais sont plus rares si d sort de ce cadre. Or, des distances fréquemment utilisées pour les données textuelles, sonores ou bien visuelles, sont des distances définies à partir d'un noyau k, qui ne sont en général pas des distances Euclidiennes dans l'espace de départ.

Pour de telles distances, des méthodes approximatives de recherche de plus proche voisin ont été dérivées, par analogie avec le cas Euclidien. La majorité s'appuie sur des méthodes de hachage, ou plus généralement d'indexation : chaque vecteur x i est compressé en une représentation s(x i ), bien moins volumineuse que celle de départ. Cette compression est typiquement faite de façon à conserver les voisinages, et donc faire en sorte que deux voisins "proches" aient des images par s similaires tout en s'assurant que des voisins "éloignés" aient des images différentes. La recherche de plus proche voisin se fait alors de manière moins coûteuse dans l'espace compressé : la requête q est également réduite en une représentation s(q), que l'on compare ensuite aux quantités s(x i ). La recherche de plus proche voisin se fait donc dans le domaine compressé des signatures s(x i ), potentiellement de façon bien plus économe en terme de mémoire et de temps de calcul : d'une part les représentations utilisées sont plus compactes que celles de départ, d'autre part la distance entre signatures est généralement plus efficace à calculer que la distance de départ. Cette recherche approximative permet d'identifier un certain nombre de voisins proches potentiels. La recherche se termine usuellement par un réordonnancement de ces voisins approchés. ix En particulier, plusieurs méthodes de hachage binaire visant à conserver les voisinages au sens d'une distance dérivée d'un noyau ont été proposées.

Dans notre troisième contribution [START_REF] Bourrier | Nearest neighbor search for arbitrary kernels with explicit embeddings[END_REF] :

• Nous proposons, pour la recherche de plus proche voisin au sens d'une distance dérivant d'un noyau, de s'appuyer sur une approximation Euclidienne du noyau, typiquement appelée plongement explicite. De telles méthodes sont éprouvées, notamment l'Analyse en Composantes Principales (ACP) à noyaux, qui a fait l'objet de travaux à la fois théoriques et pratiques sur la précision et la réduction de la complexité de calcul de l'approximation.

• L'approximation via l'ACP à noyaux nous permet de dériver un schéma de recherche exacte de plus proche voisin pour les noyaux normalisés, s'appuyant sur l'inégalité de Cauchy-Schwarz bornant la précision de l'approximation. Nous proposons des résultats montrant le gain de cette approche dans un cas particulier, et ses limites en général. Nous évoquons une relaxation possible de l'inégalité permettant d'apprendre des bornes théoriques de précision de l'approximation en général.

• Suite à cette étape de plongement explicite, il est également possible d'appliquer, si l'on veut effectuer une recherche approximative, un schéma de compression adapté à la distance Euclidienne. En particulier, nous montrons à l'aide d'expériences que la combinaison de l'ACP à noyaux et d'un hachage Euclidien offre de meilleurs résultats que des méthodes usuelles de hachage à noyaux. Il est en outre possible, à l'aide de la méthode proposée, d'appliquer suite à l'étape de plongement explicite un schéma de compression offrant une meilleure précision que le hachage binaire : la Quantification Produit.

Utiliser un plongement explicite afin de se ramener au cas Euclidien semble donc être une étape utile, permettant notamment de dériver des schémas de recherche approximative plus précis que les méthodes de hachage existantes. Les algorithmes de recherche de plus proche voisin devraient ainsi se comparer à cette succession d'étapes simples afin de mesurer le réel apport de l'algorithme en question.

Introduction

For more than a decade, numerical sciences have known a fast evolution, reflecting the growing need for data processing tools. Practical interests for such tools are numerous: some examples of general tasks where statistical learning on data is critical are contentbased search, automatic classification of data, or collaborative filtering. This type of general tasks is addressed by many software or online services.

The gain in interest for such tools in the past years is tied to several phenomena, which led to a large increase in the size of available data, especially concerning the Internet. This increase in the number of sources and in the volume of content contributed to the elaboration of databases comprised of billions of elements e , and on which traditional statistical analysis tools cannot apply, principally because of hardware constraints, implying the computers or servers cannot cope with memory or computation time requirements of these algorithms on such voluminous data.

Particular statistical processing methods of such data have therefore arisen, principally consisting in a volume reduction of the data and/or in an approximative but faster way of dealing with computations. At the theoretical level, the compromise between computational time, precision of the statistical analysis and data size has also been studied.

In parallel to this evolution in data processing, has been developed in the signal processing community a theory of compressed sensing, a technique ensuring that under certain hypotheses on a signal, it is possible to reconstruct it precisely from a number of measures which is a priori far smaller than the dimension of the signal. Compressed sensing theory has developed with both theoretical results controlling the necessary and sufficient hypotheses on the signal for such a reconstruction to be possible, and practical algorithms to perform such a reconstruction, as well as additional results on the precision of the reconstruction provided by these algorithms.

Interestingly, analogous ideas to those used in compressed sensing have been successfully applied for statistical analysis of voluminous data. The goal was principally to reduce the size of the considered data (similarly to the measurement step in compressed sensing), thus reducing the memory and/or computational time required for an algorithm to perform, while still being able to approximately extract the desired information (similarly to the reconstruction step). Therefore, compressed sensing seems to have a genuine potential for statistical analysis of voluminous data.

The starting point of the Ph.D. thesis presented in this manuscript is therefore the study of the potential of compressed sensing for learning on data. More precisely, we will consider unsupervised learning, in which the considered data is not divided into categories before processing.

The document layout is the following: Chapter 1 provides a more thorough description of the notions and ideas at the source of this work. In particular, it draws links between compressed sensing and some techniques which have been proposed in the learning come Databases built from Facebook or Flickr image data can comprise 10 11 elements.

xvii INTRODUCTION munity. Then come three chapters, each one presenting a main contribution of the thesis. These chapters all begin by a quick review of the state of the art relative to the corresponding contribution. Then come a summary of the contributions relative to the corresponding chapter. Finally, the contributions are described.

Chapter 2 contains the first contribution: a parameter estimation method for probability density mixtures on compressed data, in a similar way as what is done in compressed sensing, and which can be interpreted as a generalized compressed sensing problem. We propose to compress data to a fixed-size representation called a sketch, then derive an algorithm to induce mixture parameters from this sketch. We provide numerical experiments evaluating the experimental performance of this compressed framework in the case of isotropic Gaussian mixture estimation.

Chapter 3 describes the second contribution : a generalization of results relative to the theoretical performance of reconstruction methods in compressed sensing, and more generally in inverse problems, to models beyond the models to which the previous results apply. More particularly, we provide necessary and sufficient conditions to the existence of an instance optimal reconstruction function, which is optimal in a certain sense with respect to the model. We study this instance optimality property in a general case and link it with a generalized Restricted Isometry Property, which is a widely-used property in compressed sensing.

Finally, Chapter 4 presents the third main contribution: an approximate nearest neighbor search of a certain element with respect to certain distances, relying on the approximation of these distances by a Euclidean distance. We provide a simple framework combining two existing dimension-reducing techniques: explicit embedding and Product Quantization.

We provide experiments showing that these frameworks allow a better precision than hashing techniques designed to work with kernels.

Chapter 1

Motivation

We introduce in this section the underlying motivation to the thesis work depicted in the present manuscript. The initial idea is to make a parallel between two a priori different domains: on the one hand, inverse problems, and more particularly the recent theory of compressed sensing, which are typically used in signal processing; on the other hand, somewhat general data analysis methods which exploit ideas similar to the paradigms of compressed sensing.

We will first, in Section 1.1, give a review of inverse problems with a particular focus on compressed sensing, describing the ideas at the core of this concept. Then we will talk, in Section 1.2, about some data processing methods related to statistical learning and relying on analogous schemes to those used in compressed sensing. This will suggest links between compressed sensing and learning, which will be developed in our first contribution in Chapter 2.

Inverse problems and Compressed sensing 1.Sparse signal representations

The linear abstraction of signal processing. Physical signals typically obey linear partial differential equations. Thanks to this linear behavior, the abstract linear algebra framework is particularly well-suited for signal processing. In this sense, a signal is typically represented:

• either as a function defined on an infinite domain (typically spatial or temporal), and belonging to a space of functions sharing a certain regularity property;

• or more frequently as a quantized signal, i.e., a finite-dimensional vector of R n or C n , where each entry represents the value of a quantization cell (in time for audio signals or in space for images for instance).

Many transformations undergone by signals can be modeled by linear operators, such as convolution operators. This allows the usage of powerful mathematical tools to model signal data and derive algorithms to process it.

Sparsity of a signal. The "canonical" representation of a finite-dimensional signal as a vector, such as the representation of an image as a concatenation of pixel values, is not necessarily well-suited for processing. An interesting property on the representation of a signal in a certain basis is sparsity [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]. When linearly decomposing a signal x on a certain basis B = {b 1 , . . . , b n }, we generically say that its representation is sparse if the decomposition involves only a few nonzero coefficients, that is x = n i=1 x i b i with x i = 0 for most subscripts i.

In this case, the basis B is well-suited to represent x: it provides a compact representation (the set of couples {(i, x i ) : x i = 0} fully describes the signal), and can give a simple interpretation by describing the signal with only a few components. This sparse representation of the signal x also allows quick linear computations, and the a priori knowledge that x is sparse when decomposed on B allows one to solve linear inverse problems, as will be discussed in the following.

Classical sparsity-inducing bases. In signal processing, given a specific task, one usually considers a certain model of signals, which can contain very specific signals or on the contrary a wide range of them. To efficiently exploit sparsity, a desirable property is that there exists a basis B which provides an approximately sparse decomposition for every vector of the considered model. Such a basis B allows one to precisely approximate all (or most of) the signals in the considered model by sparse vectors.

Such classical "sparsity-inducing" bases used in signal processing include many trigonometric-based families, initially proposed by Fourier at the beginning of the XIX th century for periodic functions [START_REF] Fourier | Mémoire sur la propagation de la chaleur dans les corps solides[END_REF]. Indeed, natural signals often present a form of periodicity associated to particular frequencies, which will be the frequencies around which the energy of the Fourier decomposition will be concentrated. Over the years, this type of frequency decomposition has been applied to other linear objects like nonperiodic functions or discrete vectors. Such trigonometric families are typically denoted as Fourier, or frequential, bases.

The main drawback of Fourier representations is the so-called Gibbs phenomenon, which designates their difficulty to adequately represent a discontinuity in the signal. To circumvent this difficulty, another family of functions is widely used as a representation basis for signals: the wavelet functions [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF], which typically are piecewise regular functions with localized discontinuities or strong variations. With such wavelet bases, the decomposition of a large class of natural signals, such as natural images, is concentrated on a small subset of vectors in the family, while the rest of the decomposition involves small coefficients.

Both Fourier and wavelet families, which include many different variants, are usually well-structured sets of vectors, which allows one to derive efficient algorithms to compute the corresponding decomposition of a vector.

More recently, other sparsity-inducing families of vectors have been considered and are denoted as dictionaries. They are typically families of vectors learned from a data set to adaptively provide a sparse decomposition on data which is similar to the training set. Unlike Fourier and wavelet families, orthogonality or even linear independence is not usually required, hence their usual denomination of redundant dictionaries. This type of representation, while being potentially much better suited for specific data than a nonadpative family of vectors, suffers from harder learning and decomposition steps compared to the previously mentioned representations, which is mainly due to the lessstructured form of the family.

Inverse problems

The concision and interpretability given by sparse representations make them useful for many signal processing tasks. A common considered framework in signal processing is the case where a vector s * ∈ R n undergoes the action of a linear operator M, with a possible additional noise e. This results in a vector Ms * + e, from which one would want to recover s * . This framework applies to multiple problems, among which appear denoising, deconvolution or super-resolution a .

In general, the operator M will not be invertible, so that even in the absence of noise, the problem of recovering s * from the quantity y = Ms * + e is ill-posed: several vectors can indeed have the same image by M, impeding from recovering the correct preimage. To correctly define the problem, one must impose additional constraints on s * . One of the typical considered constraints is that s * has a sparse linear representation in a certain basis, or more generally a certain dictionary D. This assumption implies that the vector s * can be written as Dx * , where D is a matrix containing the vectors of D as columns, and x * is a sparse vector containing an unknown but "small" number of nonzero entries.

In this case, the problem of "inverting" the observation y = MDx * + e to find x * can be expressed as the search for the sparsest vector x such that MDx is close to the observation y in a certain sense. More precisely, the problem can be formulated as:

argmin x x 0 , subject to y -MDx 2 ≤ σ, (1.1)
where x 0 is the pseudo-norm counting the number of nonzero entries of x, and σ is the estimated order of magnitude of the noise e, which is typically taken in the 2 -norm sense.

Note that there may be multiple solutions to the problem (1.1), so that it may be considered as ill-posed. However, if σ = 0, there exist sufficient conditions so that the solution is unique. In this case, the problem is expressed as the search for the sparsest vector of the affine space {x ∈ R n : MDx = y}. If x * is sparse enough b , then the solution to (1.1) is unique (Chapter 2 of [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF]).

If σ > 0, uniqueness is never achieved (unless in the trivial case where x * = 0): indeed, if x is a solution of (1.1), the continuity of the operator associated to MD ensures that any vector sufficiently close to x is also solution. Since there exists an infinite number of vectors with the same sparsity as x and arbitrarily close to x, there is an infinite number of solutions to (1.1). Therefore, in this case where σ > 0, one simply aims at finding one of the many solutions of (1.1) (Chapter 5 of [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF]).

However, even when the problem (1.1) has a unique solution, it is in general NP-hard to solve it [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF][START_REF] Davis | Adaptive greedy approximations[END_REF], therefore the problem is asymptotically intractable. Two main approaches can be considered to try and solve this problem approximately:

• Greedy approaches: Since one searches for a sparse vector x such that y ≈ MDx, y will be similar to a linear combination of a few columns of MD. Therefore, one can find columns of this matrix which are well-correlated with y, so that there is a linear combination of these columns which approximately yields y, hopefully leading to a good solution of the problem. The research for such columns can be done in one step by finding the most correlated columns (Gribonval et al., 2007), or in an iterative fashion with a pursuit algorithm [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] which updates at each step a sparse approximate solution by adding a nonzero entry corresponding to a column of MD which is well-correlated with the residual.

a We consider real vectors for simplicity, but all considered results can easily be applied to complex vectors.

b Namely, if x * 0 is less than half the spark of the matrix MD, the spark being the minimal number of linearly dependent columns of a matrix.

Even if these approaches can yield an exact solution of the problem (1.1) in particular cases, the conditions are usually very restrictive and not satisfied in practice (Chapter 4 of [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF]). Nonetheless, these greedy algorithms can still perform well in practice and yield a good approximate solution (Chapter 3 of [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF]).

• 1 relaxation: One can replace in (1.1) the term x 0 by the convex term x 1 . The problem then becomes

x * = argmin x x 1 , subject to y -MDx 2 ≤ σ. (1.2)
This formulation is typically known as Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] and can be solved using standard optimization techniques by casting it as a linear program (Chapters 13 and 14 of [START_REF] Nocedal | Numerical Optimization[END_REF]) or a second order cone program (Chapter 4 of [START_REF] Boyd | Convex Optimization[END_REF]). They can typically be solved using interior point methods. When σ = 0, under stronger assumptions on M, D and x * 0 than in the • 0 formulation, it has been proven that the solution to the problem (1.2) is unique and the same as problem (1.1) [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF][START_REF] Donoho | For most large underdetermined systems of linear equations the minimal 1 -norm solution is also the sparsest solution[END_REF]. When σ > 0, there also exist robustness results upper bounding the discrepancy between the solutions of (1.1) and (1.2) [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF].

Both types of methods therefore provide feasible alternatives to the intractable problem (1.1), as well as theoretical results controlling the precision of the approximate solutions with respect to the solutions of the initial problem under certain conditions.

We have introduced inverse problems in the case where the operator M is fixed and one aims at reconstructing a signal from its (possibly corrupted) measure by M. For such inverse problems, one looks for an adequate model and/or an adequate optimization formulation in order to perform the inversion. The next section focuses on another way of considering inverse problems, leading to the theory of compressed sensing.

Compressed sensing paradigms

Motivation for compressed sensing. Typically, the acquisition of a signal is performed by sampling it at regularly spaced points. In this framework, the Nyquist-Shannon theorem is a well-known result precising the number of regular measurements one needs to be able to theoretically recover the signal. It states more precisely that if a function f has a low-frequency spectrum, that is the Fourier transform of f is contained in a segment [-ω, ω], then f can be (theoretically) perfectly recovered from a discrete regular sampling of period π ω . In this case, the reconstruction process is a simple convolution step with a sinc function.

In some practical acquisition cases of a signal, there are benefits to reducing the number of measurements, even if this causes the reconstruction step to be more costly. For instance, in medical imaging, the number of measurements should be minimal to avoid a negative impact on the patient. Another example where a reduced number of measurements would be beneficial is when sensing a signal using a device with limited energy, such as a mobile device or a satellite.

This objective of reducing the number of measurements allows one to look at linear inverse problems under another angle. One now considers a given model Σ ⊂ R n , which is a set of signals of interest. Given these particular signals, one aims at finding linear operators M : R n → R m , with m "small" (hopefully way smaller than n) and such that for any x ∈ Σ, one can recover x from measurements Mx with an algorithmic procedure (typically numerical optimization). As well as in the usual inverse problem discussed previously, one can also study a noisy case where the measurements are corrupted with a noise e, and where one aims at approximately recovering x from measurements Mx + e.

The somewhat recent but already well developed theory of compressed sensing [START_REF] Kutyniok | Compressed sensing: Theory and applications[END_REF][START_REF] Qaisar | Compressive sensing: From theory to applications, a survey[END_REF] gives interesting and surprising answers to this measurement reduction objective.

The initial compressed sensing setup. Early works in the compressed sensing theory [START_REF] Donoho | Compressed sensing[END_REF]Candès and Tao, 2006) mainly consider the noiseless recovery of sparse vectors in the canonical basis of R n . Let's denote by Σ k the set containing the k-sparse vectors: .3) In this compressed sensing framework, the goal is to find operators M satisfying several conditions:

Σ k = {x ∈ R n : x 0 ≤ k} . ( 1
• Dimension reduction: The operators must map R n in a low-dimensional space R m , that is m < n, and hopefully m n.

• Theoretical recovery of Σ k : The operators must yield theoretical noiseless recovery of all vectors in Σ k . As we have seen in the previous section, the typical linear inverse problem way of casting the recovery of a signal

x * ∈ Σ k from y = Mx * is as argmin x∈R n x 0 , subject to Mx = y, (1.4)
which is the problem (1.1) reformulated in the case where D is the identity matrix of dimension n and σ = 0. The problem (1.4) should have x * as a unique solution for the theoretical recovery condition to be satisfied.

• Practical recovery of Σ k : The operators must yield practical noiseless recovery of all vectors in Σ k . We have seen that the problem (1.1) is practically intractable in general. Therefore, an additional practical recovery condition is that this problem can be recast as an equivalent problem, which is easier to solve and still yields the same solution. This equivalent problem is typically the 1 relaxation argmin

x∈R n x 1 , subject to Mx = y. (1.5)
An optional extra constraint is to enforce that problems (1.4) and (1.5) yield precise reconstruction (in a certain sense) of approximately sparse vectors x * , that is vectors which are close to Σ k in a certain sense. This additional constraint comes from the fact that practical signals are not usually exactly sparse but rather have fast decaying coefficients in specific basis: they are therefore approximately sparse and not exactly sparse, and the "precise recovery" property should encompass these signals.

Early results of compressed sensing. Initial results of compressed sensing have been stated in a quick succession of papers from 2004 to 2006, and many theoretical results have kept on being found ever since. In the next few sections, we propose to summarily unknot the braid of this dense succession of fundamental papers, presenting the main flavor of the results in the approximate order they were proposed.

All the initial papers about compressed sensing mainly give examples of linear operators M satisfying the conditions mentioned in the previous section. The proposed operators have two aspects in common, even if there are slight variations between the results:

• Drastic dimension reduction: The results prove that there are a lot of operators M which drastically reduce the initial dimension n while satisfying the other compressed sensing prerequisites. More precisely, the exact reconstruction of Σ k is usually possible with m = O(k ln(n)), which is negligible with respect to n when n goes to infinity c . Moreover, still with m = O(k ln(n)), the reconstruction error of certain approximately sparse vectors x * is upper bounded by the 2 distance from x * to Σ k , that is the residual 2 -norm if one approximates x * by its best k-term approximation.

• Randomized choice of the operator: Interestingly, the proposed operators are not constructed in a deterministic manner. Instead, all the works consider a simple class M of operators, supplied with a probability distribution. The operator used to perform compressed sensing is randomly drawn in M, and theoretical results prove that with "high probability" on the drawing of M in M, M will satisfy the compressed sensing constraints. The term "high probability" usually refers to the fact that when n → ∞, if m is chosen as O(k ln(n)), then the probability that M satisfies the reconstruction conditions is 1-O(n -α ), where α is a parameter depending on the considered signals, the model M and the desired reconstruction guarantees.

In particular, the considered classes M of operators, along with the corresponding probability distributions leading to compressed sensing guarantees, were:

• In (Candès et al., 2006), the authors considered the case where M was a submatrix of the Fourier matrix F, defined by

Fr, s = 1 √ n exp -2πi (r -1)(s -1) n . ( 1.6)
M was obtained by uniformly extracting m different lines of F. In this case, the authors proved the theoretical and practical recovery of Σ k with high probability provided m = O(k ln(n)).

• The even more groundbreaking papers [START_REF] Donoho | Compressed sensing[END_REF] and (Candès and Tao, 2006) prove similar results, and are often viewed as the founding works in compressed sensing. They consider the reconstruction of vectors of R n with coefficients obeying a power-law decay, that is signals x * for which the q th largest entry is upper bounded by Cq -α with C, α > 0. They prove that if m = O(k ln(n)), then with high probability on the drawing of the linear operator M, which will be detailed in the following, a signal x * belonging to this class of signals can be reconstructed up to a precision of the order of the best k-term approximation of x * . The reconstruction can once again be cast as a 1 minimization problem (1.5).

Even if their choices of operators is similar, there are slight differences between the two papers. In [START_REF] Donoho | Compressed sensing[END_REF], the author considers a set M of m × n matrices with unit-norm columns, provided with uniform distribution. In (Candès and Tao, 2006), the authors again consider Fourier subsampling matrices, as well as matrices with entries drawn i.i.d. either from a symmetric Bernoulli distribution (taking values ±1 with probability 1 2 ) or from a normal distribution.

These surprising results show that compressed sensing is feasible with particular dimension reducing operators. As has been mentioned, all these operators share in common the interesting property of being chosen randomly among a certain family. In (Candès and Tao, 2006), the authors also precise some common properties shared between the different considered families M of operators, which pins down the reason why they are well-suited for compressed sensing. These conditions will then be simplified into a powerful property described in the next section.

The Restricted Isometry Property The direct legacy of the aforementioned founding results is a series of two papers [START_REF] Candès | Decoding by linear programming[END_REF]Candès et al., 2006) which mainly introduce a simple sufficient condition on an operator M for it to satisfy the compressed sensing constraints (and even more). This property is called the Restricted Isometry Property (RIP) which is a property concerning a certain sparsity level t. This RIP states that for all t-sparse vectors x,

(1 -δ) x 2 2 ≤ Mx 2 2 ≤ (1 + δ) x 2 2 (1.7)
for a certain 0 < δ < 1. The RIP therefore says that even if M is a highly dimensionalreducing operator, so that the column vectors which appear in M are highly linearly dependent, any subfamily of t such column vectors is nonetheless approximately orthogonal, so that M acts approximately as an isometry on t 2 -sparse vectors (supposing t is even), hence the name "Restricted Isometry".

This property was first precised and used in [START_REF] Candès | Decoding by linear programming[END_REF] in a noncompressive setting, then applied to compressed sensing in (Candès et al., 2006). The success of the recovery procedure in a usual compressed sensing setting can be linked to the RIP:

• The previously mentioned random matrices asymptotically satisfy the RIP with high probability: this includes Gaussian, Bernoulli and Fourier subsample matrices.

• A matrix M satisfying some kind of RIP is well-suited for compressed sensing. Intuitively, let's suppose that M satisfies the RIP on the set of 2k-sparse vectors. If

x and y are k-sparse vectors, then the distances xy 2 and Mx -My 2 will be close: the ratio between the two distances will be between 1 -δ and 1 + δ. This makes intuitively understandable the recovery of k-sparse vectors: if δ is small, two k-sparse vectors cannot be mixed up in the compressed domain.

This intuitive power of the RIP is theoretically justified: the results of (Candès et al., 2006) prove that if M satisfies certain RIPs, it satisfies the previously mentioned compressed sensing conditions, and moreover permits reconstruction from noisy measurements. Indeed, under some assumptions on RIP constants, if x * ∈ Σ k and y = Mx + e is a noisy measurement with e 2 ≤ σ, the minimization problem

argmin x∈R n x 1 , subject to y -Mx 2 ≤ σ (1.8) yields a solution x # such that x * -x # 2 ≤
Cσ, that is the reconstruction error is of magnitude similar to the noise of the measurement. Moreover, the reconstruction error of an approximately sparse vector with noisy measurements can also be upper bounded under RIP assumptions.

The main advantages of the RIP are its simplicity and the fact that it separates the probabilistic choice of an operator from the theoretical implications of the choice of a well-suited operator. Because of this, it has been widely used in the compressed sensing theory to prove theoretical reconstruction results. The simplicity of this property is particularly well illustrated in [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF], where the author concisely proves noiseless and noisy compressed reconstruction bounds for approximately sparse vectors under RIP assumptions.

Note that more recently, compressed sensing has been considered directly in a probabilistic framework, allowing one to get reconstruction results under "RIPless" conditions on M (Candès and Plan, 2011). The corresponding results directly give conditions on a family of operators M and on the random choice of these operators for compressed sensing to be feasible with high probability.

The importance of randomness. The main characteristic of compressed sensing is the random choice of the measurement operator. Indeed, even if the previously mentioned results prove the existence of multiple well-suited operators for compressed sensing, they do not provide a way to construct them in a deterministic way. There is actually no known deterministic way of building compressed sensing operators in polynomial time which have guarantees comparable to those obtained in a probabilistic manner.

The random choice of the operator, although strange at first glance, is a good way of finding a suitable operator for the initial purpose of compressively acquiring the signal: it is fairly useless to really aim for an operator which will be guaranteed to work, and one can instead pick randomly an operator among a simple model in which most of the elements are suitable. This good behavior of randomly drawn operators is theoretically linked to the concentration of measure phenomenon, which designates a collection of results underlying the fact that when a large number of random variables is drawn, the variability of a wide range of statistical quantities computed from these variables will be very low, so that the observed statistics will very likely be close to their expected value.

This idea of picking an element at random instead of struggling to find the "optimal" element while still having theoretical guarantees will be discussed in the next section, through examples of methods used in statistical analysis or optimization and aiming at reducing complexity and/or memory usage for different tasks.

"CS-like" techniques in statistical analysis

In parallel to the development of compressed sensing, and actually before, the idea of using randomness to save computational time or memory has been exploited in techniques more related to statistical analysis than signal processing.

We give some examples of such techniques, shedding some light on the similarities with compressed sensing. Apart from the following mentioned techniques, sketch-based methods share common points with compressed sensing and some of them will be discussed in Chapter 2.

Locality-Sensitive Hashing

Databases are often comprised of vectors living in the same vector space. When statistically manipulating a database, one often needs to perform a similarity search among these vectors, i.e., looking for the nearest neighbors of a certain vector in the database. When the dimension and/or the number of vectors is large, it may become compulsory to rely on approximate but faster schemes. To this end, hashing can be exploited.

Hashing for approximate neighbor search. The core idea of hashing for computational savings is to replace the initial vectors x i , belonging to a metric space with metric d, by shorter signatures h(x i ) belonging to a simpler space with a metric d which is easier

x 1 . . . x L Initial vectors ∈ (R n , d) Hashed vectors ∈ (S k , d ) h 1 (x 1 )
. . .

h k (x 1 )
. . .

h 1 (x L )
. . .

h k (x L ) = h(x 1 ) = h(x L ) H Hash functions : R n → S h 1 h k . . .
(drawn i.i.d. )

Figure 1.1: Core idea of LSH. Each initial vector is mapped to a signature for which each entry corresponds to the image of the vector by a randomly drawn hash function.

and faster to compute. Even though there is no immediate link between d and d , the hash functions are chosen so that d approximately preserves neighborhoods: if x and y are "close neighbors" in the initial space, we should have h(x) ≈ h(y) and conversely, if they are far apart, so should h(x) and h(y) be.

Locality Sensitive Hashing (LSH) [START_REF] Indyk | Approximate nearest neighbors: Towards removing the curse of dimensionality[END_REF][START_REF] Gionis | Similarity search in high dimensions via hashing[END_REF][START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF]) is a particular hashing scheme which was introduced to provide a fast approximate similarity search method for distances corresponding to the p-norm in dimension n:

x p = n i=1 |x i | p 1 p
.

(1.9)

LSH relies on a family H of elementary hash functions mapping the initial space R n to a finite set of elements S. The signature of a vector x is a k-tuple h(x) = (h 1 (x), . . . , h k (x)) ∈ S k , where the h j are elements of H drawn i.i.d. with respect to a certain probability distribution on H. This general framework is illustrated in Figure 1.1.

For the signature to be well-designed, the family H should be locality-sensitive, i.e., for any couple of vectors of R n , the probability that they are mapped to the same element in S (with respect to the drawing of the hash function) should be high if the vectors are close to one another, and low if they are far apart. More precisely, this property is usually stated as follows: for all x, y ∈ R n ,

• If d(x, y) ≤ r 1 , then P(h(x) = h(y)) ≥ p 1 . • If d(x, y) ≥ r 2 , then P(h(x) = h(y)) ≤ p 2 .
(1.10)

Here, r 1 , r 2 , p 1 and p 2 are parameters of the locality sensitive property and P denotes the probability with respect to the drawing of the hash function h among the family H.

Obviously, for this property to have a true meaning, one must have r 1 < r 2 and p 1 > p 2 .

To make a parallel with compressed sensing, one may compare this property to the RIP: indeed, several hashing families satisfy this property, which is useful in showing theoretical guarantees for algorithms which use LSH, as we will mention in the following. Among the families which satisfy this locality-sensitive property appears a family of binary functions h θ indexed on θ ∈ R n , with values in {0, 1} and defined as h θ (x) = 1 if θ, x ≥ 0 and h θ (x) = 0 otherwise [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]Charikar, 2002a). This family will be further discussed in Chapter 4, in a more general hashing setting. Other usually considered hash functions families are comprised of projections on a line followed by partitioning into segments [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF].

Note that LSH, and more generally hashing techniques aimed at reducing complexity, can be used in two different ways. The first basic way is to compute a single signature of respectable size for every database vector, then compute a full nearest neighbor search in the signature space, getting a list of nearest neighbors. This yields an easily implementable research scheme. This first way of exploiting hashing will be used in Chapter 4. On the other hand, one may compute several smaller signatures for a vector, each belonging to a different bucket. In this second case, a lot of vectors may have similar signatures in a given bucket. The nearest neighbor list of a "query" vector is computed by identifying the vectors which have similar signatures as the query in most buckets. This method relies on a harder implementation to be efficient, but can yield sublinear research times [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF][START_REF] Andoni | Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions[END_REF]. This second framework is usually the considered framework when one aims at proving theoretical guarantees for approximate nearest neighbor search.

In both cases, when a list of probable nearest neighbors is formed, an additional step of reranking is performed, where the true distances between the query and the considered vectors are computed in order to better identify the true nearest neighbors.

Theoretical results about LSH can be proved: if a family H satisfies the localitysensitive property stated above, the main result one can obtain is that with suitable choices for the number of buckets and the dimension of the hashed vectors, LSH can be used to efficiently solve (with high probability on the drawing of the hash functions) the so-called -nearest neighbor search problem, that is given a query q and its nearest neighbor x in a database X , find a vector y ∈ X such that d(q, y) ≤ (1 + )d(q, x).

(1.11)

Similarity with compressed sensing. Even though the frameworks of compressed sensing and LSH are very different, there are interesting similarities in the procedures and the theoretical study: in both cases, one randomly draws a dimension reduction operator among a certain model, and it can be proven that in suitable conditions, with high probability on the drawing of the operator, the compressed representation can be used to precisely perform a certain task.

Compressed methods for matrix decomposition

Complexity of matrix decompositions. Linear algebra is a core tool in statistical analysis. In particular, data can often be represented as a set of vectors of the same dimension, concatenated into a matrix. Typically, each column represents an individual and each line is a feature, which is a numerical quantity characterizing some aspect of the individual. When considering such a matrix, it is often useful to simplify its representation through a certain decomposition. One of the most commonly considered decomposition in the Singular Value Decomposition (SVD), which provides optimal low-rank decompositions in the Frobenius sense. SVD and other matrix decompositions are numerically performed by iterative algorithms. They typically require extensive access to the matrix, so that the matrix should be small enough to fit in main memory. Moreover, if the matrix is of size m × n, the usual complexity of the decomposition of the matrix is O (mn min(m, n)). If one only needs a rank-k approximation of the matrix without performing a full decomposition, other iterative algorithms can be used which have complexity O(mnk). These algorithms with reduced complexity still require extensive access to the matrix.

Randomized scheme for matrix decomposition. Finding more efficient ways (either in complexity or in memory usage) than the usual iterative algorithms to perform matrix decomposition has been the topic of a lot of works ( [START_REF] Mahoney | Randomized algorithms for matrices and data[END_REF][START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] compile numerous references). In particular, [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] describes several randomized methods for computing matrix decompositions such as Singular Value Decomposition (SVD) or QR decomposition. Typically, randomness is used to quickly compute a low-rank estimate of the matrix, and then a standard decomposition algorithm is applied to this estimate.

More precisely, if the matrix N to be decomposed approximately is of size m × n, the randomized layout for computing the decomposition is composed of two stages:

1. Approximation of the action of the matrix: The first goal is to efficiently "capture" the action of the matrix, that is compute a low-rank estimate of the action of the matrix. Typically, a matrix will not have singular values of the same magnitude, so that the energy of a random vector will be mostly mapped by the matrix on the first singular directions. This phenomenon is used for randomized approximation: if one searches for a rank-k approximation, the matrix N is applied to a set of k +p random vectors (typically Gaussian vectors), where p > 0 is an oversampling parameter, usually taken constant and k. The complexity of this step is O(mnk), but can be reduced to O(mn ln(k)) by applying N to a structured set of random vectors, such as Subsample Random Fourier Transform (SRFT) [START_REF] Woolfe | A fast randomized algorithm for the approximation of matrices[END_REF]. The output is an m × (k + p) orthogonal matrix Q such that .12) where • typically denotes the Frobenius or the 2 operator norm and N k is the best rank-k approximation of N with respect to the norm . .

N -QQ T N ≈ N -N k , ( 1

Decomposition of the matrix:

Since the matrix Q allows to compute a good low-rank approximation of the main singular components of N, the desired decomposition of N can directly be performed on the much smaller matrix Q T N. For SVD, it is for instance sufficient to compute the SVD of Q T N = UΣV * and the approximate SVD of N (up to rank k) is then given by N k ≈ (QU)ΣV * .

This framework can be adapted to other decompositions such as QR. The complexity is not different from a standard partial SVD algorithm for general dense matrices, but since this layout relies mainly on matrix-vectors multiplications to reduce the size of the data, the complexity can be significantly dropped if these multiplications can be quickly computed, typically if the matrix is sparse. Finally, the scheme can also be adapted to aim at reducing the number of passes on the matrix N, in the case where N is too large to fit in main memory and the bottleneck of a decomposition algorithm is the disk access time.

Links with compressed sensing. These randomized decompositions have provable precision bounds, relying on concentration of measure tools. The results are various, but they mainly consist either in an average case analysis, upper bounding the expectation of the decomposition error over the drawing of the random vectors, or in a probabilistic bound on the error, upper bounding it with high probability on the drawing of the random vectors, similarly to compressed sensing. In these upper bounds usually appear the singular values of the matrix N for ranks > k, so that the theoretical results confirm the intuition that the randomized scheme will be useful only if the singular values of N decrease at a sufficient rate.

These randomized decomposition methods share certain common points with the compressed sensing framework: knowing that a matrix has sufficiently decreasing singular values (which can be compared to the sparse model of compressed sensing), it is possible to precisely estimate a low-rank decomposition from a compressed representation of the matrix obtained via its image from a randomly picked operator. However, a singular difference is that one still needs to access the initial matrix to perform the decomposition in the second step, whereas compressed sensing ensures the reconstruction of the initial vector simply from the compressed representation and the projection operator.

Other randomized statistical procedures

In this section, we will mention other randomized schemes that have been used in various learning frameworks, but which have fewer immediate links with compressed sensing. However, the core idea is to save computational time by "picking randomly an item" instead of "searching for the best item" to perform a task, which draws links to compressed sensing.

Linear kernel approximations Kernels (Scholkopf and Smola, 2001) are widely used functions in machine learning: they provide a considerable range of metrics between elements of a database and have nice properties leading to standard algorithms to perform learning, such as Support Vector Machines (SVM). However, since these metrics can be elaborate, the cost of computing kernel values is usually substantial, so that kernel methods cannot be applied at a very large scale.

In the past few years, there has been in the machine learning community a gain of interest for explicit embedding methods, which consist in finding an approximation of a certain kernel K as a scalar product, that is finding a certain (nonlinear) transformation Φ of the data satisfying K(x, y) ≈ Φ(x), Φ(y) . The benefit is that once the data has been transformed by Φ, the approximate kernel values are fast to compute, since the computation of scalar products relies on very optimized standard linear algebra routines. This kind of approximation will be further discussed in Chapter 4.

In [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF], the authors consider two random schemes for defining such a function Φ. They consider the case where the kernel K to approximate is a shiftinvariant kernel on R n , that is K(x, y) only depends on the difference xy.

Their first method relies on Bochner's theorem, which allows to write K(x, y) as the Fourier transform of a probability measure, that is K(x, y) can be seen as an expected value of cosines. Therefore, they propose to approximate K(x, y) as an empirical mean of cosines drawn from the proper probability distribution, which can be considered as a scalar product. Concentration of measure results provide a uniform bound on the approximation error of the kernel on a compact subset of R n with high probability on the drawing of the cosines.

The second method is less general and applies only to certain type of shift-invariant kernels. It is similar to binary hashing and consists in randomly partitioning the space into bins, ensuring the probability of collision between two vectors x, y on a compact subset of R n is proportional to K(x, y). By averaging several such hashed representations, one again gets an approximation of the kernel that can be seen as a scalar product. Moreover, a uniform bound on the approximation error can still be derived, upper bounding it with high probability on the drawing of the hash functions.

In both cases, the probability distributions used to draw the function Φ are adapted to the kernel one wants to approximate. These linear random approximations were shown to yield substantial savings on training time for several learning databases and problems [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF].

Random Kitchen Sinks. In [START_REF] Rahimi | Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning[END_REF], the authors consider a binary classification problem where the classifier is sought as the sign of a linear combination of some simple classifiers taken from a family Ω = {Φ θ : θ ∈ Θ}. This framework encompasses classification models such as SVM or boosting.

Given a training set of labeled points {(x i , y i )} L i=1 ⊂ R n × {-1, 1}, the search for such a classifier is usually expressed as the following minimization problem:

argmin θ 1 ,...,θm,α 1 ,...,αm L i=1   m q=1 α q Φ θq (x i ), y i   , (1.13)
where m is a predetermined parameter enforcing the maximal number of terms in the classifier and (•, •) is a loss function measuring the discrepancy between the predicted class and the observed class. This optimization problem is complex, especially when finding correct values for θ q , since the relationship between θ and Φ θ can be complicated.

The authors therefore propose to simplify this optimization problem by relaxing the optimization on the parameters θ q : instead of finding the best m such parameters, one could just pick m parameters θ 1 , . . . , θ m at random in Θ with respect to a probability density p. The optimization problem simply becomes a simple matter of minimizing over the weights α ∈ R m the quantity

L i=1 α T z i , y i , (1.14)
where

z i is the concatenated feature vector [Φ θ 1 (x i ), . . . , Φ θm (x i )] T .
Under regularity conditions on the elements of Ω and the loss function , with high probability on the drawing of the parameters θ q , one can upper bound the difference between the loss induced by the solution of (1.14) and the loss induced by the best classifier in a certain regular family of linear combinations of elements in Ω.

Atop this theoretical bound, experiments show that this randomized formulation yields substantial savings in training time for classification using AdaBoost.

Stochastic gradient

Another interesting example where randomness can proficiently be used to reduce learning costs is the so-called stochastic gradient method, which relies on a method initially proposed in [START_REF] Robbins | A stochastic approximation method[END_REF]. A long time after, the method has been readjusted to online learning and studied statistically [START_REF] Bottou | Online algorithms and stochastic approximations[END_REF][START_REF] Murata | A statistical study on on-line learning[END_REF].

The stochastic gradient is typically applied to problems where one aims at minimizing with respect to x an objective function of the form f (x) = L i=1 i (x), where the terms i (x) are typically loss terms of individual data of a training set. Provided all the functions i are differentiable with respect to x a usual gradient descent algorithm would need to compute all terms ∇ i to find the gradient descent direction ∇f . However, it has been proven that in a large-scale setting, it is sometimes preferable to reduce the cost of computing this direction by iteratively considering a descent direction of the form ∇ i for a particular i, which is chosen randomly among the training data in the case of stochastic gradient. While such a method obviously yields a cost reduction of a single iteration of the optimization algorithm, it has been proven theoretically that it can also offer a better time/precision tradeoff at a large scale [START_REF] Bottou | The tradeoffs of large scale learning[END_REF].

Compressed least-squares regression. In [START_REF] Maillard | Compressed least squares regression[END_REF], the authors propose a compressive framework to perform linear regression. Given data X = {(x r , y r )} L r=1 ⊂ Ω × R and a family {f s : Ω → R} n s=1 of functions called features, the initial goal is to find a vector α ∈ R n such that the function f α = n s=1 α s f s is a viable regression function, that is for all r, y r ≈ f α (x r ).

(1.15)

The authors consider the case where α is chosen in a penalized least-squares fashion, that is minimize a quantity of the type

L r=1 |y r -f α (x r )| 2 + J(α), (1.16)
where J is a regularization function.

Instead of searching for a regression function of the form f α , the authors propose to compress the set {f s } by replacing it with a set of m features {g t } m t=1 , obtained as random linear combinations of f s . The proposed random combinations follow analogous probability distributions to the distributions of [START_REF] Achlioptas | Database-friendly random projections[END_REF], aimed at keeping the scalar products between a certain number of vectors while ensuring the linear combinations are somewhat sparse in order to fasten the projections computations. The regression function is then sought as a function g β = m t=1 β t g t , with β ∈ R m . The authors derive theoretical bounds on the average risk of the optimal regression function of the compressed model. For a suitable choice of the parameter m, one can obtain similar error bounds on the approximation error as other noncompressive methods, while learning the regressing function with a reduced complexity.

Conclusion

The success of compressed sensing shows us that randomization can be very robust if applied to the right model and problem. The provided examples of computational methods using randomization to solve a learning task act as an incentive to express such a task in a framework analog to compressed sensing. This will hopefully allow to develop learning methods which combine numerical efficiency and theoretical soundness.

From this observation, we will build our first contribution in the following chapter by proposing and instantiating a compressive learning framework on a usual estimation problem.

Chapter 2

Compressive density mixture estimation

This chapter will appear in a slightly modified version in the Springer/Birkhäuser book entitled Compressed sensing and its applications.

As we have seen in the previous chapter, methods that are close to compressed sensing have been exploited to reduce complexity and/or memory costs of several learning tasks. To further study the potential of compressed sensing to learning problems, let's first present a very conceptual view of what would be a "compressive learning" framework.

Conceptual compressive learning outline. When considering a learning problem, one is interested in performing a certain task on a certain type of data. A learning procedure will usually consist in trying to fit an underlying model to the type of data one is interested in by picking a model in a parametrized set M = {M θ : θ ∈ Θ}, where Θ is a parameter set which is used to index the models in M. In order to achieve this model fitting, a learning procedure will consist in finding a parameter θ * so that the model M θ * is in some sense adequate to a training set X = {x 1 , . . . , x L } containing some data one is interested in. The computational cost of estimating such a parameter θ * will depend on the size of X , on the size of the models in M and on the algorithm used for the estimation.

A compressive framework to perform such an estimation is outlined in Figure 2.1, which represents two main ways of compressing training data in order to apply a learning algorithm to data of reduced size. The top scheme represents the case where each vector of X is compressed individually: this is performed for instance in [START_REF] Calderbank | Compressed learning : Universal sparse dimensionality reduction and learning in the measurement domain[END_REF], with a method which will be described in the next section. The bottom scheme represents the case where X will be compressed into a single representation usually called sketch, of size m which should not depend on the number L of elements in the database but rather on the complexity of the model one wants to fit (represented by the parameter set Θ) and on the task one aims at achieving after the learning. This second scheme has been instantiated in a simple estimation problem in [START_REF] Thaper | Dynamic multidimensional histograms[END_REF], which will also be discussed in the next section.

Density mixture estimation. In this chapter, we will focus more particularly on a classical unsupervised learning problem: density mixture estimation. Let's model the problem as follows: suppose the data set X comprises vectors of R n , and consider a set P of parametrized probability densities defined on R n , that is

P = p θ : R n → R + R n p θ dµ = 1, θ ∈ Θ , (2.1)
where µ is the Lebesgue measure on R n . The goal of density mixture estimation is to find a density function p on R n which satisfies two prerequisites:

• p must be written as a linear combination, or mixture, of functions in P, that is p = k s=1 α s p θs , where the quantities α s are typically positive weights satisfying k s=1 α s = 1 and θ s are the parameters of the mixture. Let's notice the similarity of this condition with the sparsity condition discussed in the previous chapter.

• X can "reasonably well" be considered as a set of vectors drawn i.i.d. with respect to a probability law of density p. This condition is typically enforced by optimizing a certain quantity representing the consistency between the data X and the probability p. Said quantity is usually closely related to the parametrization of the probability densities, such as the likelihood of the parameters P(X |(α s , θ s ) k s=1 ), representing the probability of the drawing of the data in X knowing the probability density p. Alternatively, this condition can be formulated by supposing the data X is drawn i.i.d. from a probability distribution of density f , and that one searches for the best approximation p of f (in a certain sense), p being an element of the set of sparse linear combinations of functions in P.

Since p must be written as a linear combination of a few functions of P, one can consider it as a "sparse vector" over the family P. In that sense, several works have considered the density mixture estimation problem in a linear inverse problem fashion, and will be discussed in the next section. Drawing our inspiration from these works, our goal in this chapter will be to propose an instantiation of the compressive scheme represented in bottom part of Figure 2.1, applying it to a density mixture estimation problem.

In the next section, on top of discussing the works relative to density mixture estimation in a linear inverse problem fashion, we will also mention some other "compressive learning" works which can also be considered as instantiations of the outline in Figure 2.1.

State of the art: density mixture estimation and compressed learning

In this section, we present different sets of works that motivate our contribution: on the one hand, works on density mixture estimation expressed as a linear inverse problem; on the other hand, works on "compressed learning" aiming at performing learning tasks on compressed data, some of which coming from the database literature with a concern in analyzing data streams thanks to global compressed representations called sketches. Our approach will combine both aspects.

Density mixture estimation as a linear inverse problem

To express the density mixture estimation problem in a linear inverse fashion, the works [START_REF] Bunea | Spades and mixture models[END_REF][START_REF] Bertin | Adaptive Dantzig density estimation[END_REF] consider P as a finite set {p 1 , . . . , p k }. Moreover, they also consider the alternative formulation of density mixture estimation mentioned earlier, that is X is drawn i.i.d. from a probability distribution of density f . Finally, they consider the case where all the densities p s and the density f belong to the Hilbert space L 2 (R n ).

In this case, the density mixture estimation consists in finding a sparse vector α ∈ R k , such that f ≈ f α , where f α = k s=1 α s p s . This objective can be specified in different ways in a linear inverse problem fashion. They all share the following observation: if E is the expectation with respect to the density f , then the scalar product between f and p s , defined by

f, p s = R n p s f dµ = E[p s (X)], (2.2)
can be approximated by an empirical mean obtained from the data X . This empirical estimator is defined as

Ê[p s (X)] = 1 L L r=1 p s (x r ), (2.3)
where Ê is the expectation with respect to the empirical probability distribution of the data X .

In [START_REF] Bunea | Spades and mixture models[END_REF], the authors aim at minimizing over α the quantity

f -f α 2 2 + J(α), (2.4)
where J is a penalty function that promotes sparsity, defined as a weighted version of 1 -norm depending on the family P.

The first term can be developed using the scalar product, and can be replaced without changing the solution to problem (2.4) by the term

-2E[f α (X)] + f α 2 2 .
(2.5)

From there, the authors simply replace the theoretical expectation E in (2.5) by its empirical counterpart Ê. This defines the "SPADES" estimator as

α = argmin α∈R k - 2 L L r=1 f α (x r ) + f α 2 2 + J(α). (2.6)
They further propose oracle inequalities for the estimation error under different assumptions on the family P and the choice of the penalty term J.

In [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF], the authors apply the Dantzig selector [START_REF] Candès | The Dantzig selector: Statistical estimation when p is much larger than n[END_REF] to the problem by expressing it as

argmin α∈R k α 1 , subject to |(Gα) s -Ê[p s (X)]| ≤ η s for all 1 ≤ s ≤ k.
(2.7)

In this formulation, the matrix G is the Gram matrix of the family P, for which the (s 1 , s 2 ) entry is equal to p s 1 , p s 2 . The quantity η s is an adaptive threshold relative to each density p s . The objective aims at finding a sparse solution via 1 relaxation, while the constraints ensure that the theoretical correlation p s , f α = (Gα) s is close to the empirical estimator of p s , f defined as Ê[p s (X)]. These two methods have been successfully applied to density mixture estimation for several models in dimension n = 1, 2 in [START_REF] Bunea | Spades and mixture models[END_REF] and n = 1 in [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF]. However, they suffer from two major drawbacks if one aims at applying them to higher-dimensional models. Let's give somewhat intuitive insights into the reasons for these drawbacks:

• Finiteness of the density model: The family P is supposed finite. However, several density models typically used in estimation are infinite, such as Gaussian Mixture Models (GMMs), which consider Gaussian densities indexed on a certain number of parameters depending on n. If one aims for instance at applying these methods to GMM estimation, one needs to discretize the continuous model of Gaussian densities. Since a number of R h n is required to sample a cube of side R with a mesh step h in each direction, the number of centers in a discrete representation will grow exponentially fast with n, which will not be viable computationally.

• Incoherence of the density model: Both methods rely on the fact that P is a family of incoherent densities, that is the quantities

p s 1 , p s 2 p s 1 2 p s 2 2
(2.8)

are not too close to 1 if s 1 = s 2 . This incoherence necessity prevents from decomposing f on a refined model containing many similar densities.

Therefore, if these approaches are viable and theoretically sound for density mixture estimation in small dimension (n = 1, 2, 3), they may not be applied in all generality for moderate dimensions (say, even n = 5, 10). For such dimensions, it would be more interesting to consider a continuous model of densities p θ indexed by a continuous parameter θ, such as a vector of a certain space R d . However, such a continuous model cannot be exploited by the aforementioned methods.

Keeping these limitations in mind, let's present the other set of inspiring contributions for our work, which can be viewed as compressive learning instances.

Learning with randomized dimensionality reduction and sketching

Learning with individual dimension reduction. In [START_REF] Calderbank | Compressed learning : Universal sparse dimensionality reduction and learning in the measurement domain[END_REF], the authors address the classical linear Support Vector Machine (SVM) classification problem : labeled data (X , Y) = {(x r , y r )} L r=1 ⊂ R n × {-1, 1} must be separated by a hyperplane with normal vector u satisfying y r = sign( u, x r ) for all r ∈ 1, L . Since this condition usually cannot be satisfied, one often considers a relaxation of this problem where u is sought as the minimizer of a loss

L r=1 ( v, x r , y r ) (2.9)
over all vectors v of R n . A reformulation of this problem in the case where is the socalled Hinge loss proves that v can be looked for as a linear combination of the vectors x r . Therefore, the objective function is expressed in terms of scalar products x r 1 , x r 2 .

In order to reduce the learning complexity of optimizing the objective function (2.9), the authors propose to reduce the dimension of the vectors x r by replacing them with vectors Mx r , with M : R n → R m being a dimensionality-reducing linear operator. Since the vectors x r are supposedly sparse, we have seen in Chapter 1 that it is possible to randomly choose M in order to get with high probability a RIP so that M acts approximately as an isometry on X . Therefore, we will have

x r 1 , x r 2 ≈ Mx r 1 , Mx r 2 .
(2.10)

The learning problem can thus approximately be tackled with the lower-dimensional vectors Mx r instead of the vectors x r . The authors derive upper bounds on the error incurred by the low-dimensional projection.

Data stream sketches. Data streams consist in a flow of data (x r ) r≥1 . In some models, the data stream is important in itself; in others, the data x r act simply as modification step of an underlying item x, which is the object of interest (for instance, x could be a vector and the elements x r could encode modifications of this vector such as adding or subtracting a unit in a particular entry).

Standard statistical problems involving data streams include the search for frequent items among the x r , usually called heavy hitters [START_REF] Gilbert | One sketch for all: fast algorithms for compressed sensing[END_REF][START_REF] Cormode | Methods for finding frequent items in data streams[END_REF], or more generally estimation of quantiles of the content of the stream. When aiming at obtaining such information at a given time without having to store all the data flow up to this point, it is necessary to maintain a compressed representation of the data stream which will allow one to perform the desired estimations. Such compressed representations can be deterministically built, but are sometimes built similarly to hash functions. In this case, they are called sketches.

A sketch z is usually updated each time a new element (x r ) is streamed. Examples of usual sketches include the Count Sketch [START_REF] Charikar | Finding frequent items in data streams[END_REF] and the Count-Min Sketch [START_REF] Cormode | An improved data stream summary: The count-min sketch and its applications[END_REF]. They both rely on updating z thanks to randomly chosen hash functions. They are mainly used to estimate the heavy hitters.

Compressed histogram estimation. Interestingly, [START_REF] Thaper | Dynamic multidimensional histograms[END_REF] proposes a sketching procedure which can be linked to compressed sensing and to density mixture estimation. In this work, the authors consider a data stream (x r ) r≥1 , where each x r is a n-dimensional vector taken in a finite set A ⊂ R n . The goal is to obtain, at a given point r 0 , a histogram H r 0 approximating the distribution of vectors x r for r ≤ r 0 .

In order to avoid storing and updating a complete histogram of the data stream as it flows, the authors propose instead to build and update a sketch of such a histogram. This sketch is obtained by considering a low-dimensional projection of a histogram H by a randomly built linear operator M designed to approximately keep distances between histograms, still in a way similar to [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space[END_REF][START_REF] Achlioptas | Database-friendly random projections[END_REF]. The sketch can be updated at each time r by considering x r as a histogram H r which is null everywhere except in the bin corresponding to x r , where it is 1.

The benefit of this framework is to reduce memory costs, since the whole histogram needs not be updated, while still being able to compute at any time a good approximation of the histogram. However, the main drawback is the complexity of the recovery procedure, which is exponential in the dimension n of the data. This prevents from applying this method to even moderate dimensions (say, n = 10).

Layout of the chapter

In Section 2.3, we build an instantiation of the second scheme of Figure 2.1 to a density mixture estimation problem. We identify several conditions that such a compressive framework must satisfy to succeed in proposing a correct answer to the problem of compressive mixture estimation. We provide a precise framework to tackle this problem in the case where the considered probability densities are isotropic Gaussians. We will also prove that it is possible to (deterministically) linearly compress any sparse mixture of k isotropic Gaussians in dimension n with less than 8k 3 n measurements.

In Section 2.4, we propose an algorithm derived from Iterative Hard Thresholding IHT (Blumensath and Davies, 2009a) to solve the considered compressive estimation problem. We propose an analysis of the memory usage of the algorithm and discuss its computational cost.

In Section 2.5, we apply our compressive algorithm to synthetic data. We compare the precision of the compressive estimation to the estimation obtained with a standard EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] with respect to two usual pseudo-distance between densities. We experimentally show that the compressive scheme achieves similar precision as an EM algorithm while not requiring to store the learning data, thus leading to memory savings when the data is numerous.

Compressive estimation framework

Let's recall that we want to consider the following problem in a compressive way: let X = {x r } L r=1 be vectors of R n , supposedly i.i.d. from a certain probability distribution of density p ∈ L 1 (R n ). In the rest of this chapter, the notation Σ + k (P) will denote the positive linear combinations of k densities in P, that is

Σ + k (P) = k s=1 λ s p θs : λ s ∈ R + , θ s ∈ Θ . (2.11)
Let's note that P ⊂ Σ + k (P). Our goal is to find a good estimate of p in the set Σ + k (P). To simplify the problem, we will suppose that p is an exact k-sparse mixture of densities taken in P, defined in (2.1), that is p = k s=1 α s p θs , with α s ≥ 0, k s=1 α s = 1 and θ s ∈ Θ. In this case, the most natural way to approximate p with a density of Σ + k (P) is to try and estimate α = (α 1 , . . . , α k ) and θ 1 , . . . , θ k from X . We further want to perform this estimation in a compressive fashion. By analogy with compressed sensing, let's derive a conceptual method.

The compressive operator.

The unknown "signal" p is a sparse linear combination of the p θ 's. One would like to reconstruct p from a "compressive" measure Mp, where M is a linear operator which transforms a density function into a finite-dimensional representation, so that one is able to manipulate it. To reconstruct p from Mp, one will then look for an element q of Σ + k (P) satisfying Mq ≈ Mp. The main raised issue for now is to design an adequate linear measurement operator M. This operator should satisfy the following requirements:

• Estimation of Mp: The empirical representation of the density p is the discrete density associated to the collection of vectors X . Since p is unknown, one cannot compute directly Mp, and so the operator M must be such that Mp can be estimated through this empirical distribution.

• Computation of MP: Intuitively, the reconstruction algorithm will aim at finding densities of Σ + k (P) which will have an image by M similar to the empirical value of Mp computed from X . To reconstruct the density, one should therefore be able to compute "easily" the value of Mf for any f ∈ Σ + k (P) (or for any f ∈ P, which is equivalent). "Easily" essentially means one must have a closed-form expression of Mf .

Suppose the operator M transforms a function into a compressed representation of dimension m. Operator M can be seen as the concatenation of m linear forms M 1 , . . . , M m . Since we made the simplifying assumption that p belonged to Σ + k (P) ⊂ P , where P denotes the complex span of P, one only needs to define the linear forms M j on P .

Considering the complex span instead of the real positive span will allow us to simplify the expressions of the linear forms we consider in the following, which are based on Fourier transform.

These linear forms must satisfy the two above conditions. Simple linear forms on P can be defined as

M g : f → R n f g dµ, (2.12)
where g is a bounded measurable function on R n . In particular, the required conditions are easily interpreted for M g :

• Estimation of M g p: Since p is a probability density on R n ,

M g p = R n gp dµ = E[g(X)], (2.13)
where E[•] is the expectation with respect to the probability law of density p. Therefore, the value of M g p can be approximated by the empirical estimate

Mg (X ) = 1 L L r=1 g(x r ) = Ê[g(X)], (2.14)
where Ê[•] is the expectation with respect to the empirical distribution with L equal masses 1 L at each vector of X = {x 1 , . . . , x L }. For such a choice of linear forms, it is therefore possible to estimate M g p for virtually any g. Moreover, concentration of measure results such as Hoeffding's inequality or McDiarmid's inequality provide confidence intervals on the estimation error.

• Computation of MP: The functions g should be chosen so that the value of M g p θ is computable in closed form for any θ ∈ Θ.

A compression scheme analog to compressed sensing will consist in considering a family G of functions so that g.p θ is integrable for any θ ∈ Θ and g ∈ G. Having defined a probability distribution on the family G, one will be able to randomly choose a compressive operator M by drawing i.i.d. m functions g 1 , . . . , g m ∈ G and defining M = (M g 1 , . . . , M gm ).

Proposed instantiation: isotropic Gaussians.

Let's now propose a particular instantiation of this framework which we will use in the rest of this chapter. Let's consider σ ∈ R * + , which will be fixed for the rest of the chapter. Given n ∈ * , let's define the considered family of densities as

P n = p µ : x → 1 (2π) n 2 σ n exp - x -µ 2 2 2σ 2 , µ ∈ R n . (2.15)
This family contains all isotropic Gaussians of R n with variance σ 2 , indexed by their ndimensional mean µ, which uniquely characterizes them. As before, we define Σ + k (P n ) as the linear combinations of k (or less) functions of P n and we adopt for the rest of the chapter the simplifying notation Σ + k (P n ) = Σ + k,n . Natural linear forms associated to this type of functions are Fourier measurements. Each Fourier measurement can be indexed by a frequency vector ω ∈ R n and the corresponding function g can be defined as g ω (x) = exp (-i x, ω ). The corresponding linear form, denoted as M ω , is therefore:

M ω : q → R n q(x)e -i x,ω dx.
(2.16)

Given data X = {x 1 , . . . , x L }, the empirical counterpart of this linear form is

Mω (X ) = 1 L L r=1 exp (-i x r , ω ) .
(2.17)

A compressive operator can therefore be defined by choosing m frequencies ω 1 , . . . , ω m and posing M = (M ω 1 , . . . , M ωm ). These frequencies will typically be randomly drawn, so that for a probability density f ∈ P n , Mf can be interpreted as the sampling of the characteristic function of f at m random frequencies. We will make precise the random choice of the frequencies we retained in section 2.5.

The value of M ω p µ is explicitly computable: one has

M ω p µ = exp - σ 2 2 ω 2 2 exp (-i µ, ω ) .
(2.18)

Injectivity of the compressive operator.

Let's now prove that for some deterministic choices of less than 8k 3 n frequencies, the corresponding sketching operator M is injective on Σ + k,n , and even on Σ k,n , which is the set of all linear combinations of k densities of P n (non-necessarily with positive coefficients).

Even though the frequencies we will consider in our experiments will not be deterministically chosen (and will be far less numerous), the following results prove that reconstructing a function of Σ + k,n from a finite number of Fourier measurements is a theoretically well-posed problem provided the frequencies are well-chosen. The proofs of the following theorems can be found in Appendix A.

In dimension 1. We first consider the case where n = 1, from which we can deduce the generalization for larger dimensions. In this case, the injectivity can be obtained with 4k 2 well-chosen frequencies (which is indeed inferior to the number of 8k 3 announced earlier): Theorem 1. Let ω 1 , . . . , ω 2k ∈ R\{0} be such that ∀p = q, ωp ωq = Q. Let's define M be the linear operator on L 1 (R) which performs a Fourier sampling at the following 4k 2 frequencies:

• 2k first multiples of ω 1 : ω 1 , 2ω 1 , . . . , 2kω 1 , • 2k first multiples of ω 2 : ω 2 , 2ω 2 , . . . , 2kω 2 , • . . . • 2k first multiples of ω 2k : ω 2k , 2ω 2k , . . . , 2kω 2k .
Then M is injective on Σ k,1 , and in particular on Σ + k,1 .

Let's note that for any k, there exists 2k frequencies ω s satisfying the condition ωp ωq = Q for all p = q. This comes from the fact that R is a vector space of infinite dimension over Q (if it were finite-dimensional, R would be isomorphic to Q L as a Q vector space for a certain integer L, and thus would be countable).

In dimension n > 1. For n > 1, let's keep the choice of ω 1 , . . . , ω 2k of Theorem 1. The following theorem states that there exists m ≤ 2kn vectors (u j ) m j=1 of R n such that taking the measurements at frequencies ω s u j with (s, j) ∈ 1, 2k × 1, m yields injectivity on Σ k,n . The total number of measurements is therefore 4k 2 m ≤ 8k 3 n. Theorem 2. Let n > 1 and ω 1 , . . . , ω 2k ∈ R\{0} be such that ∀p = q, ωp ωq = Q. There exists m ≤ 2kn vectors (u j ) m j=1 of R n such that if M is the linear operator on L 1 (R n ) which performs a Fourier sampling at frequencies ω s u j with (s, j) ∈ 1, 2k × 1, m , then M is injective on Σ k,n , and in particular on Σ + k,n .

Recovery problem formulation.

Given the data X , we will denote ẑ the empirical sketch of X , which is the m-dimensional vector ẑ = ( Mω 1 (X ), . . . , Mωm (X )). One can express the recovery of the initial density p from ẑ as the following minimization problem:

p = argmin q∈Σ k 1 2 ẑ -Mq 2 2 . (2.19)
Despite being nonconvex, this kind of formulation of the problem is addressed in a regular compressed sensing setting by greedy algorithms. In the next section, we will derive from such a standard method an algorithm aimed at solving (2.19).

Compressive reconstruction algorithm

To address the estimation problem (2.19), we propose an algorithm analogous to Iterative Hard Thresholding (IHT) (Blumensath and Davies, 2009a).

Reminder of Iterative Hard Thresholding

IHT is a standard greedy method aimed at solving sparse inverse problems, as mentioned in Section 1.1.2. Consider a k-sparse signal x of dimension n and a measurement matrix M of size m × n (with m < n). Denoting y = Mx the measurement of x, IHT considers the minimization argmin

z∈Σ k y -Mz 2 2 , (2.20)
where Σ k is this time the set of k-sparse vectors of R n . At each iteration, IHT updates an estimate x of x, decreasing the objective function φ : x → 1 2 y -Mx 2 2 while ensuring the k-sparsity of x. The quantity r = y -Mx is named the residual. The update step is performed in two steps:

1. The n-dimensional gradient of φ at current iterate x, noted ∇φ, is computed.

2. The update is given by x ← H k (x -λ∇φ), where λ is a descent step and H k is a hard thresholding operator which keeps only the k entries of the vector with largest module and sets the others to 0.

This algorithm has been proven to converge to the global minimum of (2.20) under a RIP assumption on M (Blumensath and Davies, 2009a). We now adapt this scheme of work to our reconstruction problem (2.19).

Proposed continuous case algorithm

We also adopt an iterative greedy method to perform the reconstruction of p. We therefore iteratively update an estimate p, which is parametrized by a vector â ∈ R k of positive weights and by the support Γ = { μ1 , . . . , μk } ⊂ R n corresponding to the means of the current estimated Gaussians. The current residual is defined by r = ẑ -Mp. In our case, the function φ takes p as an argument and is defined as

φ(p) = 1 2 ẑ -Mp 2 2 .
(2.21)

There are some differences between our density estimation problem and the problem addressed by IHT which require modifications of the procedure. They are explained in the following sections. The algorithm is then more precisely described.

The "gradient": continuous version. In IHT, the signal one wants to reconstruct is supposed to be sparse in a finite basis of vectors. The gradient computed in the first step is a finite-dimensional vector, and each entry corresponds to the infinitesimal shift of the objective function φ when a certain entry of the vector is shifted.

In our case, the density is supposed to be sparse in the infinite basis P σ , which is parametrized by R n . The "canonical" directions in which p can be shifted are therefore also parametrized by R n , and the "gradient" is the collection of these possible shifts, which are noted ∇ µ φ for all µ ∈ R n and defined as follows:

∇ µ φ(p) = ∂ ∂t 1 2 ẑ -M(p + tp µ ) 2 2 t=0 = -Mp µ , r . (2.22)
Again, this quantity represents the local variation of the objective function (2.19) when an infinitesimal fraction of the density p µ is added to the current estimate. Since we cannot compute these values for every µ ∈ R n , we must only choose a finite number of µ for which we will compute ∇ µ φ.

Since we aim at decreasing φ, these directions should be chosen so that ∇ µ φ(p) is negatively minimal, so that p µ is a seemingly good candidate to be added to the current estimate p. Therefore, we seek instead a certain number M of local minima of µ → ∇ µ φ(p), which will typically be chosen as O(k). These local minima parametrize elements of P which are the best correlated elements to the residual r. They are the best directions in which to "move" locally the estimate p in order to decrease the objective function φ. These local minima are searched for by a randomly initialized minimization algorithm. When they are found, they are added to the current support Γ, increasing its size up to M + k elements.

Hard Thresholding. The second step in IHT consists in choosing a descent step used to shift the current estimate in the direction of the gradient, and enforcing sparsity through hard thresholding. In our case, we have at this point an updated collection of candidate means Γ, and we want to keep only k of these means.

In order to do this, we aim at decomposing the empirical sketch ẑ as a positive linear combination of vectors in Mp ν , with ν ∈ Γ, that is to project the sketch ẑ on the cone generated by the sketches of the functions {p ν : ν ∈ Γ}. This cone can be noted C(M Γ) and defined by

C(M Γ) =    K j=1 λ j Mp ν j : N > 0, λ j ≥ 0, ν j ∈ Γ   . (2.23)
The aforementioned projection of ẑ on C(M Γ) is expressed as the following minimization problem, supposing Γ contains K vectors ν 1 , . . . , ν K :

argmin β∈R K + ẑ -Nβ 2 , (2.24)
where N is the concatenation of the sketches of the functions parametrized by Γ, that is

N = [Mp ν 1 . . . Mp ν K ] . (2.25)
The hard thresholding step is then performed by keeping the k largest coefficients and the k corresponding parameters of Γ found in (2.24). Note that in the framework we consider (isotropic Gaussians with Fourier measurements), the quantity Mp ν 2 do not depend on ν, that is the sketches all have the same energy. In the case where they do not have the same energy, one should keep the k coefficients such that β j Mp ν j 2 is maximal.

Gradient descent step. In IHT, an iteration stops when hard thresholding is performed. In our case, we can still perform an additional step, which consists in decreasing further the objective function ϕ.

At this point in the iteration, p is defined as

k s=1 αs p μs , ( 2.26)
where the parameters âlpha s and μs hopefully estimate the real parameters α s and µ s of p.

Since the family P σ is extremely coherent, the local minima we found in the previous steps may be shifted from the true mean vectors because of the imprecision induced by the other components of the mixture. This imprecision on the µ s obviously also imply imprecision on the α s . However, there may exist a better estimate for p in the vicinity of p.

To find it, we simply consider ϕ as a slightly different version of φ: ϕ represents the same cost function, but takes the parameters α 1 , . . . , α k and µ 1 , µ k as arguments. Therefore, it is defined as:

ϕ : R k × (R n ) k → R (α, µ 1 , . . . , µ k ) → 1 2 ẑ -[Mp µ 1 . . . Mp µ k ]α 2 2 .
(2.27)

Initializing the parameters to the current estimators α and μ1 , . . . , μk , we can apply a gradient descent algorithm on ϕ to find, in the vicinity of p, a better estimate in the sense that it has a smaller image by φ.

Algorithmic scheme. Algorithm 1 summarizes the overall procedure. It mainly consists of three steps by iteration:

1. M local minima of µ → ∇ µ (p) are sought with a gradient descent algorithm with random initialization and are added to the current support Γ.

2. The sketch ẑ is projected on M Γ with a positivity constraint on the coefficients. Only the k highest coefficients and the corresponding vectors of the support are kept.

3.

A gradient descent algorithm is applied to further decrease the objective function with respect to the weights and support vectors.

Algorithms 2, 3, 4 and 5 describe the subfunctions.

Algorithm 1 Compressive isotropic Gaussian mixture parameter estimation Input: Sketch ẑ, operator M, target sparsity k, integer M .

Initialize Γ = ∅, r = ẑ. repeat Set ν 1 , . . . , ν M ← Find min(M, r, M). Set Γ ← Γ ∪ {ν 1 , . . . , ν M }. Set â ← Proj cone(ẑ, Γ ). Set â, Γ ← Hard threshold(â , Γ , k). Set â, Γ ← Shift support(ẑ, â, Γ). Set r ← ẑ -k j=1 αj Mp μj . until Stopping criterion is satisfied Return â, Γ. Algorithm 2 Find min(M, r, M) For i = 1 to M Find a local minimum ν i of the function: ν ∈ R n → -Mp ν , r
with a gradient descent algorithm, initialized randomly. End For Return ν 1 , . . . , ν M .

Memory usage

Let's now estimate the order of magnitude of the memory required by the compressive algorithm to estimate p from ẑ. Let's consider that n, k and m are much larger than 1. If we suppose that optimization algorithms only use first-order quantities, their memory costs are dominated by O(kn). The computation of the cost function of Algorithm 5 requires O(km). The storage of the operator M (via the frequencies ω j ) requires O(mn).

Therefore, the total memory usage is O((k + n)m + kn) and does not depend on the number L of vectors. In comparison, the memory requirement of EM is O(L(k + n)) to store both the vectors and their probabilities to belong to each current component of the mixture. The compressed algorithm allows memory savings as soon as m+ kn k+n L. Since kn m, this condition is nearly equivalent to m L.

This suggests that one will be able to make memory savings if the number of vectors in the training set X is larger than the size of the sketch required to performed the reconstruction.

Computational complexity

Computational complexity is the main drawback of the compressed procedure, since this procedure relies on several optimization steps, which can involve many variables for large k and n.

More precisely, the computational bottleneck is the last step of the iteration where a gradient descent is performed. This optimization procedure involves k(n+1) variables and the cost for computing the function at a certain point is O(mk). Therefore, since a first order optimization algorithm requires the computation of the gradient for each variable, the complexity of a simple gradient descent implementation is O(mk 2 (n + 1)). Moreover, the cost of computing the sketch from the training data is O(mL), and must be taken into account unless the data is streamed so that the sketch can be computed "on the fly" before

Algorithm 3 Proj cone(v, Γ = {ν 1 . . . ν K })
Solve the following convex optimization problem: a = argmin

β∈R K + ||v -Nβ|| 2 , with N = [Mp ν 1 , . . . , Mp ν K ].
Return a.

Algorithm 4 Hard threshold(a, Γ = {u 1 , . . . , u K }, k) Let a i 1 , . . . , a i k be the k highest entries of a.

Return (a i 1 , . . . , a i k ), {u i 1 , . . . , u i k }. Algorithm 5 Shift support(ẑ, a, Γ = {µ 1 , . . . , µ k }) Find a local minimum (a , µ 1 , . . . , µ k ) of the function: R k × (R n ) k → R + (b, ν 1 , . . . , ν k ) → ẑ -[Mp ν 1 , . . . , Mp ν k ]b 2 , using a gradient descent algorithm initialized at a, µ 1 , . . . , µ k . Return a , {µ 1 , . . . , µ k }.
the estimation procedure per se. The overall cost is virtually O(mL + mk 2 (n + 1)n comp ), where n comp is the number of iterations performed in the compressed estimation procedure.

This must be compared to a standard EM algorithm, which has complexity O(kLn EM ), where n EM is again the number of iterations performed in the EM algorithm. In the case the sketch is not computed on the fly, there will be a gain in complexity in the compressed case only if m kn EM . This was not the case in our experiments (described in the next section), since the EM algorithm converged quickly enough so that m ∼ kn EM . We will further discuss the computational outlooks in Section 2.6.

Experiments

Experimental setup

To evaluate the behavior of the compressive reconstruction algorithm, we conducted experiments on vectors drawn from a mixture of k isotropic Gaussians with identity covariance matrices (σ = 1). In each case, we drew weights uniformly on the simplex a , and we chose the Gaussian means µ j by drawing random vectors, each entry being drawn from a probability law of density N (0, 1).

The experiments were performed in the following way: after the choice of the probability distribution p, we drew L random vectors from this probability distribution and computed the empirical sketch of the distribution in one pass of the data. The training samples were then discarded from hard memory. We chose the sketching operator M randomly, following the scheme described in Section 2.5.2. We then applied the reconstruction algorithm to the empirical sketch ẑ to get an approximated mixture p. The random initialization for the reconstruction algorithm is detailed in section 2.5.3.

To evaluate the quality of the estimation, we relied on two usual discrepancy measures between probability density functions. They were used to quantify the difference between the true mixture p and the estimated mixture p. The two considered quantities are defined a We also performed experiments where all the weights were equal to 1 k and this didn't alter the conclusions drawn from the experiments.

by integrals, which in our case could not be computed explicitly. Therefore, we approximated the integrals by empirical means: we drew N = 10 5 points (y i ) N i=1 i.i.d. from p and computed the empirical estimates described below. The two chosen measures were:

• Kullback-Leibler (KL) divergence: A symmetrized version of KL divergence can be defined as

D KL (p, p) = R n ln p(x) p(x) p(x) + ln p(x) p(x) p(x) dx.
(2.28)

The empirical estimate we considered is defined as

DKL (p, p) = 1 N N r=1 ln p(y r ) p(y r ) + p(y r ) p(y r ) ln p(y r ) p(y r ) . ( 2

.29)

The KL divergence ranges from 0 to +∞, lower values meaning closer distributions.

• Hellinger distance: The Hellinger distance can be defined as

D H (p, p) = 1 - R n p(x)p(x)dx.
(2.30)

The empirical estimate we considered is defined as

DH (p, p) = 1 - 1 N N r=1
p(y r ) p(y r ) .

(2.31)

The Hellinger distance ranges from 0 to 1. Here again, lower values mean closer distributions.

Choice of the frequencies

Let's now describe the heuristic we considered to randomly choose the compressive operator M. Let's recall that p ∈ Σ k , so that p = k s=1 α s p µ s , with the α s positive which sum to 1. Denoting by F(f ).ω the Fourier transform of a function f taken at frequency ω, we have .32) This upper bound on the value of F(p) gives a hint on the way to design the random choice of frequencies. Indeed, we want to sample frequencies which are likely to be "energetic", so that |F(p).ω| is not "too low". Designing a random choice with a density proportional to the upper bound found in (2.32) seems like a reasonable choice. Since σ = 1, the random sampling we chose for the frequencies ω j is from a probability law of density N (0, Id). Real and reconstructed centroids (respectively represented as circles and squares) of 4 Gaussians in dimension 2 from 10 3 points drawn from the mixture. To estimate the 12 real parameters of the mixture, the data was compressed to a complexvalued sketch of dimension 30, represented to the right as a 60-dimensional real signal.

|F(p).ω| = k s=1 α s p µ s ≤ k s=1 α s |F(p µ s ).ω| = k s=1 α s exp - σ 2 2 ω 2 2 = exp - σ 2 2 ω 2 2 . ( 2

Heuristic for random initialization

The search for local minima in Algorithm 2 was initialized randomly by exploiting a measure performed during the construction of the sketch: during the single pass on the data, the norms of the vectors x r are computed and the maximum of the norms, R = max x∈X ||x|| 2 , is computed. These calculations have a negligible impact on the computation time of the sketch, and on its size (it only adds one component which "completes" the sketch). The knowledge of R allows us to delimit a ball in which the centers of the Gaussians are very probably contained.

We performed the random initialization by drawing a direction uniformly on the unit sphere and multiplying this unit vector by a scalar uniformly drawn in [0; R]. For each choice of parameters, 10 experiments were performed and the depicted value is the Hellinger distance such that 80% of the experiments lead to a smaller Hellinger distance. We can essentially observe that the Hellinger distance gradually decreases as the number of mixture components rises. For the considered parameters range, choosing m = 10kn, i.e., choosing m so that it contains 10 times more values than the number of parameters to estimate, leads to a Hellinger distance smaller than 0.03 for 80% of the cases. Note that the number of measurements we took for the experiments is way below the deterministic number of 8k 3 n provided in Section 2.3.3, suggesting that a random choice of frequencies can be proven as a robust choice allowing the drastic reduction of the number of measurements compared to the proposed deterministic choice. Table 2.1 compares our algorithm with a standard EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] in the case where n = 20, k = 10, m = 1000 for values of dataset size L ranging from 10 3 to 10 5 . For each case, we can see that the precision of the estimation increases with the number of samples. In the compressed case, this can be explained by the fact that the components of the sketch are better estimated with more points. We notice that the memory used for EM is proportional to the number L of samples in the dataset, while the memory required by the compressed algorithm does not depend on this parameter, which leads to a substantial improvement in memory usage for L ≥ 10 4 . Even with this reduced memory cost, the compressed algorithm is able to provide a precision comparable to the precision of the EM algorithm.

Results

Conclusion and outlooks

In this chapter, we first proposed a review of techniques reminiscent of inverse problems and compressed sensing applied to certain learning tasks, mainly to density estimation. Our contribution consisted in a framework for density mixture estimation, which was instantiated to isotropic Gaussians and random Fourier sampling. We could derive an L Compressed KL div.

Hell. Mem. 10 3 0.68 ± 0.28 0.06 ± 0.01 0.6 10 4 0.24 ± 0.31 0.02 ± 0.02 0.6 10 5 0.13 ± 0.15 0.01 ± 0.02 0.6 L EM KL div.

Hell. Mem. 10 3 0.68 ± 0.44 0.07 ± 0.03 0.24 10 4 0.19 ± 0.21 0.01 ± 0.02 2.4 10 5 0.13 ± 0.21 0.01 ± 0.02 24 Table 2.1: Comparison between our compressed estimation algorithm and an EM algorithm in terms of precision of the estimation and memory usage (in megabytes). Experiments were performed with n = 20, k = 10, m = 1000. In each cell, the value is a median on 10 experiments with the standard deviation for the precision measures.

algorithm which experimentally shows good reconstruction properties with respect to a standard estimation algorithm. Let's now mention some outlooks related to these results.

Extension to richer families of densities

Density mixture estimation using isotropic Gaussians can be seen as a clustering problem.

It would be particularly interesting to extend the experimental results to more general families of Gaussians, for instance with diagonal covariance matrices, to allow for variations in the form of the "clusters".

Considering larger families of densities would probably require finer choices for the sketching operator M, in order to be able to separate the compressed representations of all vectors of the enriched family. For this purpose, it seems interesting to investigate multiscale sketches, where coarser frequencies may be used at the beginning of the reconstruction algorithm to avoid getting too many local minima, and higher frequencies can be added throughout the algorithm to approach the solution.

Computational savings

If such a compressive framework has the potential to reduce memory requirements by computing the sketch on the fly with streamed data, the computational complexity of the density reconstruction is still large, especially due to the last stage of the algorithm (Algorithm 5), which is an optimization step involving (n + 1)k variables and which is performed at each iteration in the current state of the algorithm. Algorithmic savings could be made by finding a faster alternative to this step. The cost of computing the sketch may also be reduced by finding better sketching operators or procedures. In particular, adopting a "multiscale" paradigm, where several sketches are computed at different ranges of frequencies may provide faster convergence to an acceptable solution.

Let's also note that the size of the sketch does not depend on the number of training vectors, that is the complexity of the estimation do not increase in the case where the sketch is updated using more training vectors (this does not include the time which is necessary to build the sketch). Moreover, the algorithm should be more precise when the number L of training vectors increases, since the empirical sketch ẑ is in this case a more reliable estimate of the theoretical sketch Mp. The general approach described in this chapter may therefore be compared to conceptual results suggesting that for a certain learning task, if an estimation algorithm is well chosen in an abstract collection of algorithms, then the computational complexity required to achieve a given estimation error ratio should decrease when the size of the training set increases [START_REF] Shalev-Shwartz | Using more data to speed-up training time[END_REF].

Theoretical analysis

We have proved that some deterministic choices of the frequencies yielded the injectivity of the sketching operator M on Σ k , so that the reconstruction of a density from a sketch is theoretically well-posed. The next step in studying the well-posedness of such a problem would be a finer study of Σ k and the action of M when the frequencies are chosen at random, with a Gaussian distribution such as the distribution chosen for the frequencies.

Investigating the geometry of Σ k and using "union bound" techniques such as in classical compressed sensing, it seems possible to study the average amount of energy of the function one is able to keep by randomly sampling the Fourier transform. This type of result typically allows deduction of robustness results about the operator M.

In particular, the regularity of Gaussians seems interesting to study this "conservation of energy" and it may be conclusive to study the behavior of M on the Schwartz space.

Chapter 3

Performance of decoders in linear inverse problems

Note: The main part of this chapter is taken from the preprint Fundamental performance limits for ideal decoders in high-dimensional linear inverse problems, which was submitted to IEEE Transactions on Information Theory on November 2013. A version of this preprint can be found at: http://arxiv.org/abs/1311.6239.

Inverse problems have been introduced in Chapter 1 in the case where one wants to invert a linear operator applied to a sparse vector in a certain basis. On top of this standard case, inverse problems concern a wide range of other signals: one can be interested in adding hypotheses to the usual sparse model if they have more information than just sparsity about the considered signals or conversely relax some hypotheses on the sparse model to build a model containing more objects. Inverse problems can even be considered for models which are structured differently than the set of sparse vectors.

Models beyond sparsity for inverse problems. Even if the sparse vectors model still attracts a lot of attention, other works on inverse problems or compressed sensing considering more exotic signal models have multiplied over the years [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF]. In these "generalized" models, the signals of interest live in or close to a subset Σ of the space. Such models are more adequate to describe particular type of signals or objects, and some of them are represented in Figure 3.1. A list of various models which have been considered in inverse problem or compressed sensing setting is:

• Block-sparse signals [START_REF] Eldar | Block-sparse signals: uncertainty relations and efficient recovery[END_REF] are vectors divided into blocks and for which only a few blocks are nonzero. In practice, this model can be used to model signals with a multi-band spectrum [START_REF] Mishali | Blind multiband signal reconstruction: Compressed sensing for analog signals[END_REF], or the difference of gene expression levels between two samples [START_REF] Parvaresh | Recovering sparse signals using sparse measurement matrices in compressed dna microarrays[END_REF].

• Dictionary-sparse signals [START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF]. Sparse signals models can be considered relatively to a certain basis, but also to a dictionary, with possibly linearly redundant vectors but which represent more adequately the considered signals. In particular, these models can be used for various signal processing tasks such as denoising or deblurring.

• Cosparse signals [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF] are signals which are sparse after undergoing a linear transform Ω, typically dimension-increasing. This framework can model overcomplete transforms of a signal which supposedly yields sparsity, such as the shift invariant wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF] or the finite difference operator for an image.

• Low-rank matrices [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]Candès and Plan, 2011) often appear when considering data which live near a subspace of limited dimension relatively to the ambient dimension. Moreover, instead of only considering matrices which are approximately low-rank, one can model the signals of interest as the sum of a low-rank matrix and of a sparse matrix [START_REF] Zhou | Stable principal component pursuit[END_REF]. In order to circumvent the ambiguity of such a decomposition (a matrix can be low-rank and sparse), a nonsparsity constraint on the low-rank matrix can be considered (Candès et al., 2011a).

• Unions of subspaces [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF][START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF] can be considered in a general setting, which encompasses a lot of other models. In particular, the usual k-sparse model in dimension n is the union of n k subspaces corresponding to all the possible supports of a k-sparse vector. Most other sparse models can be considered as finite union of subspaces. Furthermore, some models of low-rank matrices can be viewed as an infinite union of subspaces [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF].

• Low-dimensional manifolds have been considered in a compressed sensing setting [START_REF] Baraniuk | Random projections of smooth manifolds[END_REF][START_REF] Eftekhari | New analysis of manifold embeddings and signal recovery from compressive measurements[END_REF]. Such models can apply to the representation of small variations of an image [START_REF] Wakin | The multiscale structure of non-differentiable image manifolds[END_REF][START_REF] Wakin | A manifold lifting algorithm for multi view compressive imaging[END_REF].

• Symmetric definite positive square matrices with k-sparse inverse. Highdimensional Gaussian graphical models have a covariance matrix of this type: the numerous pairwise conditional independences that characterize the structure of such models, and make them tractable, translate into zeros entries of the inverse covariance matrix (the concentration matrix). Combining sparsity prior on the concentration matrix with maximum likelihood estimation of covariance from data, permits to learn jointly the structure and the parameters of Gaussian graphical models (so called "covariance selection" problem) [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF][START_REF] Yuan | High dimensional inverse covariance matrix estimation via linear programming[END_REF].

• Low-dimensional embedding of a point cloud consists in finding a projection of the points onto a low-dimensional subspace while approximately keeping distances between points. This is asymptotically feasible while substantially reducing the dimension, as proved by Johnson-Lindenstrauss lemma [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space[END_REF]. As in compressed sensing, such a projection can be achieved with high probability with usual random drawings of matrices [START_REF] Achlioptas | Database-friendly random projections[END_REF].

All these models share a common point: they typically contain far fewer vectors than the whole space. Indeed, these models serve as a priori information on the considered signals in order to enable their reconstruction from a non-injective operator M. The well-posedness of an inverse problem thus comes from the information that the signal to reconstruct x is sufficiently close to a subset Σ of the signal space. For this information to be useful, the model Σ needs to be a substantially small part of the space for the inverse problem to make sense, otherwise the number of solutions of the problem will be too large to have a meaning. Since these models generalize the sparse model, the following question arises: can they be considered under a general framework, sharing common reconstruction properties? The following discussion precises what type of results would be expected in this generalized framework.

What should we expect for an inverse problem to be well-posed? The wellposedness of an inverse problem is not a universal property. There are numerous ways to define when a problem is well-posed: either theoretically or practically, either uniformly or probabilistically. A common requirement for well-posedness is that M is injective on Σ, so that the problem is theoretically well-posed: if x ∈ Σ, the knowledge of Mx yields the knowledge of x since there is one, and one only, element of Σ having the image Mx by M.

However, this theoretical well-posedness given by the single injectivity on Σ is usually not sufficient practically for several reasons:

• Stability to the model : The model Σ do not contain exactly the vectors one is interested in recovering. It is only comprised of approximations of the signals of interest. Therefore, M should be such that one can recover sufficiently precisely vectors which do not exactly belong to Σ but "live near" Σ in a certain sense.

• Robustness to noise: The measure Mx cannot be known with infinite precision. Aside from the obvious fact that this quantity is usually considered numerically, so that a quantization step occurs and involves a precision loss, linear operators M considered in inverse problems represent in general physical filters applied to a signal, such as transformations that a signal undergoes when captured or processed by a device. In this case, chances are that the measure is corrupted by a certain amount of noise, so that the measure is of the form Mx + e, where e is the additive noise term.

• Practical recovery algorithm: Probably the most important issue, one needs to be capable of practically solve the inverse problem with all the aforementioned constraints. The existence of a unique solution is therefore not sufficient to ensure the reconstruction: one needs an algorithm which will find the signal x within a given precision.

Expressing conditions under which a signal model can be recovered with these additional constraints is much more difficult and involves a part of subjectivity. In this chapter, we first review a well-posedness formulation previously proposed [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] and called Instance Optimality. This condition was originally defined for the model of sparse vectors and characterizes a uniform stability to the model. The contributions of this chapter concern the generalization of these results to other models as those previously mentioned, as well as the study of Instance Optimality in the general case.

State of the art: Instance Optimality in the sparse case

Let's consider the vector space E = R n and the set of k-sparse vectors Σ k . Define on E a linear measurement operator M mapping the signal space in R m . The problem of inverting M for signals in Σ k can be seen as the definition of a decoder ∆ : R m → R n such that ∀x ∈ Σ k , ∆(Mx) = x, thus making M a linear encoder associated to the (typically nonlinear) decoder ∆.

As has already been mentioned in the previous section, a good decoder ∆ is certainly expected to have nicer properties than simply reconstructing Σ k , the first of it being the stability to the model: the signal x to be reconstructed may not belong exactly in Σ k but "live near" Σ k under a distance d, meaning that d(x, Σ k ) = inf z∈Σ k d(x, z) is "small" in a certain sense. In this case, one wants to be able to build a sufficiently precise estimate of x from Mx, that is a quantity ∆(Mx) such that x -∆(Mx) is "small" for a certain norm • . This stability to the model has been formalized into the so-called Instance Optimality assumption on ∆. ∆ is said to be instance optimal if:

∀x ∈ R n , x -∆(Mx) ≤ Cd(x, Σ k ), (3.1)
for a certain choice of norm • and distance d. For this property to be meaningful, the constant C must not scale with n and typically "good" instance optimal decoders are decoders which involve a constant which is the same for all n (note that this implicitly relies on the fact that a sparse set Σ k ⊂ R n can be defined for any n). When the norm is 2 or 1 and the distance is 1 , such good instance optimal decoders exist and can be implemented as the minimization of a convex objective [START_REF] Donoho | Compressed sensing[END_REF]Candès and Tao, 2006;[START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF] under assumptions on M such as the Restricted Isometry Property (RIP). Note that Instance Optimality is a uniform upper bound on the reconstruction error, and that other types of bounds on decoders can be studied, particularly from a probabilistic point of view [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF]. Other early work include upper bounds on the reconstruction error from noisy measurements with a regularizing function when the signal belongs exactly to the model [START_REF] Engl | Regularization of Inverse Problems[END_REF].

In [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], the authors considered the following question: Given the encoder M, is there a simple characterization of the existence of an instance optimal decoder? Their goal was not to find implementable decoders that would have this property, but rather to identify conditions on M and Σ k under which the reconstruction problem is ill-posed if one aims at finding an instance optimal decoder with small constant. The existence of a decoder ∆ which satisfies (3.1) will be called the Instance Optimality Property (IOP). The authors proved that this IOP is closely related to a property of the kernel of M with respect to Σ 2k , called the Null Space Property (NSP). This relation allowed them to study the existence of stable decoders under several choices of norm . and distance d(., .).

More precisely, the authors considered two norms • G and • E defined on the signal space R n , the distance derived from • E being denoted d E . Instance Optimality with respect to these two norms rewrites as follows: a decoder ∆ : R m → R n is said to be instance optimal for k-sparse signals if

∀x ∈ R n , x -∆(Mx) G ≤ Cd E (x, Σ k ), (3.2)
for some constant C > 0. This property on ∆ upper bounds the reconstruction error of a vector, measured by • G , by the distance from the vector to the model, measured by d E . The authors prove that the existence of an instance optimal decoder, called IOP, is closely related to the NSP of M with respect to the set Σ 2k of 2k-sparse vectors. Noting N = ker(M), this NSP states

∀h ∈ N , h G ≤ Dd E (h, Σ 2k ) (3.3)
for some constant D.

The relationship between the IOP and the NSP is the following: if there exists an instance optimal decoder ∆ satisfying (3.2), then (3.3) holds with D = C. Conversely, if (3.3) holds, then there exists a decoder ∆ such that (3.2) holds with C = 2D. Such a decoder can be defined as follows, supposing M is onto:

∆(Mx) = argmin z∈(x+N ) d E (z, Σ k ), (3.4)
x + N denoting the set {x + h, h ∈ N }. The well-posedness of this definition is discussed in Appendix B.1, in the more general setting where the model is a finite union of subspaces in finite dimension. Note that for generalized models, such a decoder may not necessarily exist since the infimum of d E (z, Σ) may not be achieved, as we will discuss in the next section. This result can be seen as an "equivalence" between the IOP and the NSP, with similar constants.

On top of this fundamental relationship between IOP and NSP, a question addressed in [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] is that of the fundamental limits of dimension reduction: Given the target dimension m and desired constant C, is there an encoder M with an associated instance optimal decoder? They particularly showed that there is a fundamental tradeoff between the size of the constant C in (3.1) (with 2 norm and 2 distance) and the dimension reduction ratio m/n.

In [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], the theoretical results of [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] are generalized in the case where one aims at stably decoding a vector living near a finite union of subspaces (UoS). They also show in this case the impossibility of getting a good 2 / 2 instance optimal decoder with substantial dimensionality reduction. Their extension also covers the case where the quantity one wants to decode is not the signal itself but a linear measure of the signal.

In our thesis work, we further extended the study of the IOP to general models of signals. The global summary of our contributions is the subject of the next section.

Summary of the contributions

In our work, we consider signals of interest living in or near a subset Σ of a vector space E, without further restriction, and show that instance optimality can be generalized for such models. In fact, we consider the following generalizations:

• Robustness to noise: noise-robust instance optimality is characterized, showing somewhat surprisingly the equivalence between the existence of two flavors of noiserobust decoders (noise-aware and noise-blind);

• Infinite dimension: signal spaces E that may be infinite dimensional are considered. For example E may be a Banach space such as an L p space, the space of signed measures, etc. This is motivated by recent work on infinite dimensional compressed sensing [START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF][START_REF] Hansen | Generalized sampling and infinite dimensional compressed sensing[END_REF] or compressive density estimation [START_REF] Bourrier | Compressive Gaussian Mixture Estimation[END_REF], discussed in Chapter 2;

• Task-oriented decoders: the decoder is not constrained to approximate the signal x itself but rather a linear feature derived from the signal, Ax, as in [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF]; in the usual inverse problem framework, A is the identity. Examples of problems where A = I include:

-Medical imaging of a particular region of the body: as in Magnetic Resonance Imaging, one may acquire Fourier coefficients of a function defined on the body, but only want to reconstruct properly a particular region. In this case, A would be the orthogonal projection on this region.

-Partial source separation: given an audio signal mixed from several sources whose positions are known, as well as the microphone filters, the task of isolating one of the sources from the mixed signal is a reconstruction task where E is the space of concatenated sources, and A orthogonally projects such a signal in a single source signal space.

Instance optimality for inverse problems with general models

In the noiseless case, we express a concept of Instance Optimality which does not necessarily involve homogeneous norms and distances but some pseudo-norms instead. Such a generalized Instance Optimality can be expressed as follows:

∀x ∈ E, Ax -∆(Mx) G ≤ Cd E (x, Σ), (3.5)
where • G is a pseudo-norm and d E is a distance the properties of which will be specified in due time, and A is a linear operator representing the feature one wants to estimate from Mx. Our first contribution is to prove that the existence of a decoder ∆ satisfying (3.5), which is a generalized IOP, can be linked with a generalized NSP, similarly to the sparse case. This generalized NSP can be stated as:

∀h ∈ ker(M), Ah G ≤ Dd E (h, Σ -Σ), (3.6)
where the set Σ -Σ is comprised of all differences of elements in Σ, that is Σ -Σ = {z 1 -z 2 |z 1 , z 2 ∈ Σ}. The constants C and D are related by a factor no more than 2, as will be stated in Theorems 3 and 4 characterizing the relationships between these two properties. In particular, all previously mentioned low-dimensional models can fit in this generalized framework.

Noise-robust instance optimality

Our second contribution (Theorems 5 and 6) is to provide a noise-robust extension of instance optimality and link it to a property called the Robust NSP. Section 3.4 regroups these noiseless and noise-robust results after a review of the initial IOP/NSP results of [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF]. We show somewhat surprisingly that the existence of noise-aware instance optimal decoders for all noise levels implies the existence of a noise-blind decoder.

Infinite-dimensional inverse problems

The generalization to arbitrary vector spaces allows us to consider infinite dimensional inverse problems. Here things are not always as straightforward as in the finite dimensional setting. For example, in the theory of generalized sampling [START_REF] Adcock | A generalized sampling theorem for stable reconstructions in arbitrary bases[END_REF], even when the signal model Σ is simply a finite dimensional subspace, it can be necessary to oversample by some factor in order to guarantee stable recovery. In fact Theorem 4.1 of (Adcock and Hansen, 2011) can be read as a statement of 2 / 2 instance optimality for a specific (linear) decoder given in terms of the NSP constant of the measurement operator.

The results presented here therefore provide an extension of generalized sampling for linear models beyond 2 .

In [START_REF] Adcock | Generalized sampling and infinite-dimensional compressed sensing[END_REF] the authors also combine generalized sampling theory with compressed sensing. However, rather than developing a uniform instance optimality theory which seems the most natural extension of generalized sampling, they adopt a non-uniform approach based on (Candès and Plan, 2011). Our results should enable the development of a uniform infinite dimensional CS theory.

Limits of dimensionality reduction with generalized models

The reformulation of IOP as an NSP allows us to consider the 2 / 2 instance optimality for general models in Section 3.5. In this case, the NSP can be interpreted in terms of scalar product and we precise the necessity of the NSP for the existence of an instance optimal decoder. This leads to the proof of Theorem 8 stating that, just as in the sparse case, one cannot expect to build an 2 / 2 instance optimal decoder if M reduces substantially the dimension and the model is "too large" in a precise sense. In particular, we will see that the model is "too large" when the set Σ -Σ contains an orthonormal basis. This encompasses a wide range of standard models where a consequence of our results is that 2 / 2 IOP with dimensionality reduction is impossible:

• k-sparse vectors. In the case where Σ = Σ k is the set of k-sparse vectors, Σ contains the null vector and the canonical basis, so that Σ -Σ contains the canonical basis. Note that the impossibility of good 2 / 2 IOP has been proved in [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF].

• Block-sparse vectors [START_REF] Eldar | Block-sparse signals: uncertainty relations and efficient recovery[END_REF]. The same argument as above applies in this case as well, implying that imposing a block structure on sparsity does not improve 2 / 2 feasibility.

• Low-rank matrices [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]Candès and Plan, 2011). In the case where E = M n (R) and Σ is the set of matrices of rank ≤ k, Σ also contains the null matrix and the canonical basis.

• Low-rank + sparse matrices [START_REF] Zhou | Stable principal component pursuit[END_REF]Candès et al., 2011a). The same argument applies to the case where the model contains all matrices that are jointly low-rank and sparse, which appear in phase retrieval [START_REF] Oymak | Simultaneously structured models with application to sparse and low-rank matrices[END_REF][START_REF] Ohlsson | Compressive phase retrieval from squared output measurements via semidefinite programming[END_REF][START_REF] Candès | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF].

• Low-rank matrices with non-sparsity constraints. In order to reduce the ambivalence of the low-rank + sparse decomposition of a matrix, (Candès et al., 2011a) introduced non-sparsity constraints on the low-rank matrix in order to enforce its entries to have approximately the same magnitude. However, as shown in Lemma 3, an orthonormal Fourier basis of the matrix space can be written as differences of matrices which belong to this model.

• Reduced union of subspace models [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF] obtained by pruning out the combinatorial collection of k-dimensional subspaces associated to k-sparse vectors. This covers block-sparse vectors [START_REF] Eldar | Block-sparse signals: uncertainty relations and efficient recovery[END_REF], tree-structured sparse vectors, and more. Despite the fact that these unions of subspaces may contain much fewer k-dimensional subspaces than the combinatorial number of subspaces of the standard k-sparse model, the same argument as in the k-sparse model applies to these signal models, provided they contain the basis collection of 1-sparse signals.

This contradicts the naive intuition that 2 / 2 IOP could be achievable at the price of substantially reducing the richness of the model through a drastic pruning of its subspaces.

• k-sparse expansions in a dictionary model [START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF]. More generally, if the model is the set of vectors which a linear combination of at most k elements of a dictionary D which contains an orthogonal family or a tight frame, then Theorem 8 applies.

• Cosparse vectors with respect to the finite difference operator [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF][START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF]. As shown in [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], the canonical basis is highly cosparse with respect to the finite difference operator, hence it is contained in the corresponding union of subspaces.

• As shown in Lemma 2, this is also the case for symmetric definite positive square matrices with k-sparse inverse. The covariance matrix of high-dimensional Gaussian graphical models is of this type: the numerous pairwise conditional independences that characterize the structure of such models, and make them tractable, translate into zeros entries of the inverse covariance matrix (the concentration matrix). Combining sparsity prior on the concentration matrix with maximum likelihood estimation of covariance from data, permits to learn jointly the structure and the parameters of Gaussian graphical models (so called "covariance selection" problem) [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF][START_REF] Yuan | High dimensional inverse covariance matrix estimation via linear programming[END_REF]. In very high-dimensional cases, compressive solutions to this problem would be appealing.

• Johnson-Lindenstrauss embedding of point clouds [START_REF] Achlioptas | Database-friendly random projections[END_REF]. Given a set X of L vectors in R n and > 0, there exists a linear mapping f : R n → R m , with m = O(ln(L)/2 ) and

(1 -) x -y 2 ≤ f (x) -f (y) 2 ≤ (1 + ) x -y 2 (3.7)
holds for all x, y ∈ X . The fact that the point cloud contains a tight frame is satisfied if it "spreads" in a number of directions which span the space. In this case, one cannot guarantee precise out-of-sample reconstruction of the points in R n in the

Generalized Restricted Isometry Property

Our last contribution, in Section 3.6, is to study the relations between the NSP and a generalized version of the Restricted Isometry Property (RIP), which encompasses classical or recent RIP formulations, such as

• the D-RIP (Candès et al., 2011) for the dictionary model;

• the RIP for low-rank matrices (Candès and Plan, 2011);

• the Ω-RIP [START_REF] Giryes | Greedy-Like Algorithms for the Cosparse Analysis Model[END_REF] for the cosparse model.

This generalized RIP bounds Mx F from below and/or above on a certain set V , and can be decomposed in:

Lower -RIP : ∀x ∈ V, α x G ≤ Mx F (3.8) Upper -RIP : ∀x ∈ V, Mx F ≤ β x G , (3.9)
where • G and • F are norms defined respectively on the signal space and on the measure space, and 0 < α ≤ β < +∞. We prove particularly in Theorem 9 that a generalized lower-RIP on Σ -Σ implies the existence of instance optimal decoders in the noiseless and the noisy cases for a certain norm • E we call the "M -norm" a . Furthermore, we prove that under an upper-RIP assumption on Σ, this M -norm can be upper bounded by an atomic norm [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF] defined using Σ and denoted • Σ . This norm is easier to interpret than the M -norm: it can in particular be upper bounded by usual norms for the k-sparse vectors and low-rank matrices models. We have the following general result relating generalized RIP and IOP (Theorem 10): if M satisfies a lower-RIP (3.8) for V = Σ -Σ and an upper-RIP (3.9) for V = Σ, then for all δ > 0, there exists a decoder ∆ δ satisfying ∀x ∈ E, ∀e ∈ F ,

x -∆ δ (Mx + e) G ≤ 2 1 + β α d Σ (x, Σ) + 2 α e E + δ, (3.10)
which is a particular case of Robust instance optimality, as described in Section 3.4.

Structure of the chapter

The structure of the chapter is as follows: Section 3.4 first contains a quick review of the relationship between IOP and NSP in the usual sparse case, then exposes the more general setting considered in this paper, for which these properties and their relationship are extended, both in noiseless and noisy settings. Section 3.5 then focuses on the particular case of 2 / 2 IOP, proving the impossibility for a certain class of models to achieve such IOP with decent precision in dimension reducing scenarii. In particular, we show that this encompasses a wide range of usual models. Finally, in Section 3.6, we get back to the problem of IO with general norms and prove that a generalized version of the lower-RIP implies the existence of an instance optimal decoder for a certain norm we call the "M -norm". We propose an upper-bound on this norm under a generalized upper-RIP assumption to get an IOP with simpler norms, illustrating the result in standard cases.

The proofs of most of the results presented in this chapter are found in Appendix B.

a The prefix "M-" should be thought as "Measurement-related norm" since in other works the measurement matrix may be denoted by other letters.

Signal space E

Pseudo-norm • E (pseudo-distance d E ) Measure space F Pseudo-norm • F Σ x d E (x, Σ) Mx Mx + e MΣ Ax ∆(Mx + e) Ax -∆(Mx + e) G AΣ M ∆ A Feature space G Pseudo-norm • G Figure 3.2:
Illustration of the proposed generalized setting. The signals belong to the space E, supplied with a pseudo-norm • E used to measure the distance from a vector to the model Σ containing the signals of interest. E is mapped in the measure space F by the operator M and the measure is perturbed by an additive noise e. The space F is supplied with a pseudo-norm • F . The feature space G, supplied with a norm • G , is composed of vectors obtained by applying a linear operator A to the signals in E. These feature vectors are the vectors one wants to reconstruct from the measures in M by applying a decoder ∆. The reconstruction error for the vector x and noise e is therefore Ax -∆(Mx + e) G . Note that in the case where E = G and A = I, the decoder is aimed at reconstructing exactly the signals.

Generalized IOP and NSP equivalences

In this section, we extend the initial IOP/NSP relationship in several ways.

Proposed extensions

The framework we consider is more general. The signal space is a vector space E, possibly infinite-dimensional. In particular, E may be a Banach space such as an L p space, the space of signed measures, etc. On this space is defined a linear operator M : E → F , where F is the measurement space, which will most likely be finite-dimensional in practice. We assume that M is onto. We further define a signal model Σ ⊂ E comprising the signals which we want to be able to "reconstruct" from their images by M. In the framework we consider, this "reconstruction" is not necessarily an inverse problem where we want to recover x from Mx. More precisely, as in [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], we consider a case where we want to recover from Mx a quantity Ax, where A is a linear operator mapping E into a space G. When G = E and A = I, we are brought back to the usual case where we want to reconstruct x. This generalized framework is illustrated in Figure 3.2.

In this generalized framework, we are now interested in the concepts of IOP and NSP, as well as their relationship. A decoder ∆ : F → G will aim at approximating Ax from Mx.

The approximation error will be measured by a function

• G : G → R + . This Triangle Inequality Symmetry 0 = 0 Definiteness Homogeneity • E X X X - - • F X X X - - • G X X - - -
Symmetry : x G = -x G (3.11) Triangle inequality : x + y G ≤ x G + y G . (3.
12)

The differences with a regular norm is that neither definiteness nor homogeneity is required:

x G = 0 needs not imply x = 0 and λx G needs not equal |λ| x G . We provide two examples of such pseudo-norms in the case where G = R n : [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF].

• • G can be defined as a "non-normalized" p -quasinorm for 0 ≤ p ≤ 1, that is x G = n i=1 |x i | p . In this case, λx G = |λ| p x G . • More generally, if f : R + → R + is a concave function such that f (x) = 0 ⇔ x = 0, then • G can be defined as the f -(pseudo-)norm x f = n i=1 f (|x i |), see
In order to measure the distance from a vector to the model, we also endow E with a pseudo-norm • E : E → R + which satisfies the same properties as • G with the additional requirement that 0 E = 0. The pseudo-distance d E is defined on E 2 by d(x, y) = xy E . Yet again, • E can be defined as a non-normalized p -norm or an f -norm.

We will also consider a noisy framework where the measure Mx is perturbed by an additive noise term e. To consider IOP and NSP in this context, we measure the amount of noise with a pseudo-norm in the measurement space F , which we will denote by • F . The assumptions we make on • F are the same as the assumptions on • E .

To sum up, here are the extensions we propose compared to the framework of [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF][START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF] :

• The measure Mx can be perturbed by an additive noise e.

• The model set Σ can be any subset of E.

• E is not necessarily R n but can be any vector space, possibly infinite-dimensional.

• The reconstruction of Ax is targeted rather than that of x.

• The functions • E , • F and • G need not be norms but can be pseudo-norms with relaxed hypotheses. In particular, Table 3.1 summarizes the requirements on these functions.

48CHAPTER 3. PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

N Σ x + N • x • x 1 d(x 1 , Σ) • x 2 d(x 2 , Σ) • x 3 d(x 3 , Σ)
Figure 3.3: Necessity of the additive δ term in a simple case. For each x in the blue half-plane, the distance d E (x + N , Σ) is never reached at a particular point of x + N : the distance strictly decreases as one goes right along the affine plane x + N (d(x 1 , Σ) < d(x 2 , Σ) < d(x 3 , Σ)), so that the minimal distance is reached "at infinity".

The noiseless case

We first consider the same framework as [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF][START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], where one measures Mx with infinite precision. In our generalized framework, instance optimality for a decoder ∆ reads:

∀x ∈ E, Ax -∆(Mx) G ≤ Cd E (x, Σ).
We will prove that if IOP holds, i.e., if the above holds for a certain decoder ∆, then a generalized NSP is satisfied, that is:

∀h ∈ N , Ah G ≤ Dd E (h, Σ -Σ), with D = C. Note that the set Σ 2k has been replaced by Σ -Σ = {x -y, x ∈ Σ, y ∈ Σ}. When Σ = Σ k , we have indeed Σ -Σ = Σ 2k .
The construction of an instance optimal decoder from the NSP is more complicated and the form of the instance optimality we get depends on additional assumptions on Σ and M. Let's first suppose that for all x ∈ E, there exists z ∈ (x + N ) such that d E (z, Σ) = d E (x + N , Σ). Then the NSP (3.6) implies the existence of an instance optimal decoder satisfying (3.5) with C = 2D. If this assumption is not true anymore, then the NSP implies a slightly modified IOP, which states, for any δ > 0, the existence of a decoder ∆ δ such that: (3.13) reflecting the fact that one cannot necessarily consider the exact quantity

∀x ∈ E, Ax -∆ δ (Mx) G ≤ Cd E (x, Σ) + δ,
argmin z∈(x+N ) d E (z, Σ) but rather a certain vector z ∈ (x + N ) satisfying d E (z, Σ) ≤ d E (x + N , Σ) + δ.
A similar positive "projection error" appears in [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF].

Remark 1. To understand the necessity of such an additive error term when Σ is a general set, we can consider the following toy example depicted in Figure 3.3 where

E = R 2 , N = R × {0}, Σ = {(x 1 , x 2 ) ∈ (R + ) 2 : x 2 = 1 x 1 } and • G / • E are the 2 norm.
In this case, the minimal distance between x + N and Σ is not reached at any point, making it necessary to add the δ term for the decoder to be well-defined.

In this setting, the NSP (3.6) implies the existence of instance optimal decoders in the sense of (3.13) for all δ > 0. Moreover, this weak IOP formulation still implies the regular NSP with D = C. This is summarized in Theorems 3 and 4.

Theorem 3. Suppose ∀δ > 0, there exists a decoder ∆ δ satisfying (3.13):

∀x ∈ E, Ax -∆ δ (Mx) G ≤ Cd E (x, Σ) + δ.
Then M satisfies the NSP (3.6):

∀h ∈ N , Ah G ≤ Dd E (h, Σ -Σ), with constant D = C. Theorem 4. Suppose that M satisfies the NSP (3.6): ∀h ∈ N , Ah G ≤ Dd E (h, Σ -Σ).
Then ∀δ > 0, there exists a decoder ∆ δ satisfying (3.13):

∀x ∈ E, Ax -∆ δ (Mx) G ≤ Cd E (x, Σ) + δ, with C = 2D.
If we further assume that

∀x ∈ E, ∃z ∈ (x + N ), d E (z, Σ) = d E (x + N , Σ), (3.14)
then there exists a decoder ∆ satisfying (3.5):

∀x ∈ E, Ax -∆(Mx) G ≤ Cd E (x, Σ) (3.15) with C = 2D.
Note that this result is similar to the result proven in [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], which was stated in the case where Σ is a finite union of subspaces in finite dimension. In this framework, condition (3.14) is always satisfied as soon as • E is a norm, by the same argument as in usual CS (see Appendix B.1).

Let's also note the following property: if

• E is definite, that is x E = 0 ⇒ x = 0, then d E is a distance.
In the following proposition, we prove that if we further suppose that the set Σ+N is a closed set with respect to d E , then the NSP (3.6) implies for any δ > 0 the existence of a decoder ∆ δ satisfying (3.5) with C = (2 + δ)D. This assumption therefore allows us to suppress the additive constant in (3.13) and replace it by an arbitrarily small increase in the multiplicative constant of (3.5).

Proposition 1. Suppose that M satisfies the NSP (3.6), that d E is a distance and that Σ+N is a closed set with respect to d E . Then ∀δ > 0, there exists a decoder ∆ δ satisfying:

∀x ∈ E, Ax -∆ δ (Mx) G ≤ (2 + δ)Dd E (x, Σ).
(3.16) 50CHAPTER

PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

The noisy case

In practice, it is not likely that one can measure with infinite precision the quantity Mx. This measure is likely to be contaminated with some noise, which will be considered in the following as an additive term e ∈ F , so that the measure one gets is y = Mx + e. In this case, a good decoder should be robust to noise, so that moderate values of e should not have a severe impact on the approximation error. We are interested in the existence of similar results as before in this noisy setting. We first need to define a noise-robust version of instance optimality. The robustness to noise of practical decoders is in fact a problem that has been considered by many authors. A first type of result considers noise-aware decoders, where given the noise level ≥ 0 a decoder ∆ fulfills the following property: ∀x ∈ E, ∀e ∈ F,

e F ≤ ⇒ Ax -∆(Mx + e) G ≤ C 1 d E (x, Σ) + C 2 .
(3.17)

Here, the upper bound on the approximation error gets a new term measuring the amplitude of the noise. For example, this noise-robust instance optimality holds for a noise-aware 1 decoder in the sparse case with bounded noise [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF] for

• G = • 2 and • E = . 1 / √ k, provided M satisfies the RIP on Σ 2k .
In practical settings, it is hard to assume that one knows precisely the noise level. To exploit the above guarantee with a noise-aware decoder, one typically needs to overestimate the noise level. This loosens the effective performance guarantee and potentially degrades the actual performance of the decoder. An apparently stronger property for a decoder is to be robust even without knowledge of the noise level: ∀x ∈ E, ∀e ∈ F,

Ax -∆(Mx + e) G ≤ C 1 d E (x, Σ) + C 2 e F .
(3.18)

Further on, such decoders will be referred to as noise-blind. Guarantees of this type have been obtained under a RIP assumption for practical decoders such as iterative hard thresholding, CoSAMP, or hard thresholding pursuit, see e.g. (Foucart, 2011, Corollary 3.9). Of course, the existence of a noise-blind noise-robust decoder in the sense of (3.18) implies the existence of a noise-aware noise-robust decoder in the sense of (3.17) for any noise level . We will see that, somewhat surprisingly, the converse is true in a sense, for both are equivalent to a noise-robust NSP.

Just as in the noiseless case, dealing with an arbitrary model Σ and possibly infinite dimensional E requires some caution. For δ > 0, the noise-robust (and noise-blind) instance optimality of a decoder ∆ δ is defined as

: ∀x ∈ E, ∀e ∈ F , Ax -∆ δ (Mx + e) G ≤ C 1 d E (x, Σ) + C 2 e F + δ. (3.19)
One can see that ∆ δ necessarily also satisfies the noiseless instance optimality (3.13) by setting e = 0. As we show below, if for every δ > 0 there exists a noise-robust instance optimal decoder ∆ δ satisfying (3.19), then a generalized NSP for M relatively to Σ -Σ, referred to as Robust NSP, must hold: (3.20) with D 1 = C 1 and D 2 = C 2 . This property appears e.g. in (Foucart and Rauhut, 2013) (Chap. 4) with

∀h ∈ E, Ah G ≤ D 1 d E (h, Σ -Σ) + D 2 Mh F ,
• G = • E = • 1 and • F any norm.
Note that this Robust NSP concerns every vector of E and not just the vectors of the null space N = ker(M) b . In the case where h ∈ N , one retrieves the regular NSP. For other vectors h, another additive term, measuring the "size" of Mh, appears in the upper bound.

Conversely, the Robust NSP implies the existence of noise-robust instance optimal decoders ∆ δ satisfying (3.19) with C 1 = 2D 1 and C 2 = 2D 2 for all δ > 0. These results are summarized in Theorems 5 and 6.

Theorem 5. Suppose ∀δ > 0, there exists a decoder ∆ δ satisfying (3.19)

: ∀x ∈ E, ∀e ∈ F, Ax -∆ δ (Mx + e) G ≤ C 1 d E (x, Σ) + C 2 e F + δ.
Then M satisfies the Robust NSP (3.20):

∀h ∈ E, Ah G ≤ D 1 d E (h, Σ -Σ) + D 2 Mh F , with constants D 1 = C 1 and D 2 = C 2 .
Theorem 6. Suppose that M satisfies the Robust NSP (3.20):

∀h ∈ E, Ah G ≤ D 1 d E (h, Σ -Σ) + D 2 Mh F .
Then ∀δ > 0, there exists a decoder ∆ δ satisfying (3.19)

: ∀x ∈ E, ∀e ∈ F, Ax -∆ δ (Mx + e) G ≤ C 1 d E (x, Σ) + C 2 e F + δ, with constants C 1 = 2D 1 and C 2 = 2D 2 .
We conclude this section by discussing the relation between noise-aware and noise-blind decoders. A noise-aware version of noise-robust instance optimality can be defined where for ≥ 0, δ > 0 we require ∀x ∈ E, ∀e ∈ F,

e F ≤ ⇒ Ax -∆ δ, (Mx + e) G ≤ C 1 d E (x, Σ) + C 2 + δ. (3.21)
Of course, the existence of a noise-blind instance optimal decoder implies that of noiseaware decoders for every ≥ 0. The converse is indeed essentially true, up to the value of the constants C i :

Theorem 7. Suppose ∀ , δ > 0, there exists a noise-aware decoder ∆ δ, satisfying (3.21):

∀x ∈ E, ∀e ∈ F, e F ≤ ⇒ Ax -∆ δ, (Mx + e) G ≤ C 1 d E (x, Σ) + C 2 + δ.
Then M satisfies the Robust NSP (3.20) with constants D 1 = C 1 and D 2 = 2C 2 . Therefore, by Theorem 6, there exists an instance optimal noise-blind decoder satisfying:

∀x ∈ E, ∀e ∈ F, Ax -∆ δ (Mx + e) G ≤ 2C 1 d E (x, Σ) + 4C 2 e F + δ.

2 / 2 Instance Optimality

In this section, we suppose that E is a Hilbert space equipped with the norm • 2 and scalar product •, • , that F = R m and we consider a finite-dimensional subspace V of dimension n, on which we define the measure operator M. We are interested in the following question in the noiseless framework: Is it possible to have a "good" noiseless instance optimal decoder with • G = • E = • 2 in a dimensionality reducing context where m n? A result of [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] states that in the usual sparse setting, one cannot expect to get a good instance optimal decoder if M performs a substantial dimensionality reduction, the best corresponding constant being n m . In [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF], the authors prove that this lower bound on the constant holds in the case where Σ is a finite union of subspaces in finite dimension. Here, we are interested in a version of this result for the general case where Σ can be a more general subset of E. More precisely, we will give a sufficient condition on Σ under which the optimal 2 / 2 instance optimalityconstant is of the order of n m , thus preventing the existence of a 2 / 2 instance optimal decoder with small constant if m n.

Homogeneity of the NSP

In the case where • G , • E and • F are actual norms, the general NSP can be rewritten as an NSP holding on the cone R(Σ -Σ) generated by Σ -Σ, i.e., the set {λz|λ ∈ R, z ∈ Σ -Σ}.

Lemma 1. If • G and • E are norms, we have an equivalence between the NSP on Σ -Σ: ∀h ∈ N , Ah G ≤ Dd E (h, Σ -Σ), (3.22) and the NSP on R(Σ -Σ):

∀h ∈ N , Ah G ≤ Dd E (h, R(Σ -Σ)). (3.23)
Similarly, if • G , • E and • F are norms, we have an equivalence between the robust NSP on Σ -Σ:

∀h ∈ E, Ah G ≤ D 1 d E (h, Σ -Σ) + D 2 Mh F , (3.24)
and the robust NSP on R(Σ -Σ):

∀h ∈ E, Ah G ≤ D 1 d E (h, R(Σ -Σ)) + D 2 Mh F . (3.25)
This lemma, which is valid even in the case where A is not the identity, shows that the NSP imposes a constraint on the whole linear cone spanned by the elements of Σ -Σ and not only on the elements themselves. Note that this equivalence is trivial in the case where Σ is a union of subspaces since Σ -Σ is already a cone in this case.

3.5.2

The optimal 2 / 2 NSP constant Remark 2. In the subsequent sections of the paper, we will assume that A = I (this implies G = E), so that one aims at reconstructing the actual signal.

In the 2 / 2 case, one can give a simple definition of the optimal NSP constant D * , that is the minimal real positive number D such that the 2 / 2 NSP is satisfied with constant D: (3.26) This definition assumes that there exists some constant so that the NSP is satisfied. Using the NSP definition and Lemma 1, we get that

D * = inf {D ∈ R + |∀h ∈ N , h 2 ≤ Dd 2 (h, Σ -Σ)}.
D * = sup h∈N z∈R(Σ-Σ) h 2 h -z 2 = sup h∈N z∈R(Σ-Σ) 1 h h 2 -z h 2 2
.

(3.27)

Denoting N 1 and Σ 1 respectively the set of unit-norm vectors (in the 2 sense) of N and Σ -Σ, we can rewrite the expression above as:

D * = sup h∈N 1 z∈R(Σ-Σ) 1 h -z 2 = sup h∈N 1 z∈Σ 1 λ∈R 1 h -λz 2 .
(3.28)

A simple study gives that if

h 2 = z 2 = 1, then sup λ∈R 1 h-λz 2 = 1 √ 1-h,z 2 , so that: D * = sup h∈N 1 z∈Σ 1 1 1 -h, z 2 . (3.29)
The contraposition of Theorem 3 gives the following result : if the NSP (3.6) is not satisfied for a certain constant D, then no decoder ∆ δ can satisfy instance optimality (3.13) with constant D. In the 2 / 2 case, considering D < D * , h ∈ N ∩B 2 and z ∈ R(Σ-Σ)∩B 2 such that h, z 2 ≥ 1 -1 D 2 , we can construct two vectors such that for any decoder, instance optimality with constant < √ D 2 -1 can only be satisfied for at most one of them. This will shed light on the link between NSP and IOP. We have z

= z 1 -z 2 z 1 -z 2 2
for some z 1 , z 2 ∈ Σ. Let ∆ be a decoder. If ∆(Mz 1 ) = z 1 , then this vector prevents ∆ from being instance optimal. The same goes for z 2 if ∆(Mz 2 ) = z 2 . Now, let's suppose that z 1 and z 2 are correctly decoded. In this case, (z 1 + z 2 )/2 is decoded with a constant worse than √ D 2 -1, as depicted in Figure 3.4. Indeed, noting p = (z 1 + z 2 )/2 and defining the vectors p 1 and p 2 respectively as the orthogonal projections of z 1 and z 2 on the affine plane p + N , we must have ∆(Mp 1 ) = ∆(Mp 2 ). Denoting as p N ⊥ the orthogonal projection on N ⊥ , we have

d 2 (p 1 , Σ) ≤ d 2 (p 1 , z 1 ) = p N ⊥ (z 2 -z 1 ) 2 /2. Similarly, d 2 (p 2 , Σ) ≤ p N ⊥ (z 2 -z 1 ) 2 /2. The fact that ∆(Mp 1 ) = ∆(Mp 2 ) implies that there exists i ∈ {1, 2} such that p i -∆(Mp i ) 2 ≥ p 1 -p 2 2 /2 = p N (z 1 -z 2 ) 2 /2. Therefore, p i -∆(Mp i ) 2 d 2 (p i ,Σ) ≥ p N (z 2 -z 1 ) 2 p N ⊥ (z 2 -z 1 ) 2 ≥ D 1 -1 D 2 = √ D 2 -1.
This illustrates the closeness between NSP and IOP: a vector of R(Σ -Σ) which is correlated with N can be used to define a couple of vectors such that for any decoder, one of the vectors will not be well decoded.

3.5.3

2 / 2 IO with dimensionality reduction

Main theorem

Let's now exploit the expression of D * to state the main result of this section: if R(Σ -Σ) contains an orthonormal basis of the finite-dimensional subspace V ⊂ E (or even a family of vectors that is sufficiently correlated with every vector of V ), then one cannot expect to

N p 1 z 1 z 2 p 2 • p p + N Figure 3
.4: Illustration of the impact of the correlation between N and Σ -Σ on instance optimality. Here, z 1 and z 2 are two vectors in Σ such that z 1 -z 2 is well correlated with N , implying that at least one of the two vectors p 1 and p 2 , which are close to Σ but far from one another, will not be well decoded.

get a 2 / 2 instance optimal decoder with a small constant while M substantially reduces the dimension of V . The fact that R(Σ -Σ) contains such a tight frame implies that the dimension of N cannot be too big without N being strongly correlated with Σ -Σ, thus yielding the impossibility of a good instance optimal decoder.

Before showing examples where this theorem applies, let's first state it and prove it.

Theorem 8. Suppose V is of dimension n and Σ -Σ contains a family z 1 , . . . , z n of unit-norm vectors of E satisfying ∀x ∈ V , n i=1 z i , x 2 ≥ K x 2 2 .
Then to satisfy the NSP on V , M must map V into a space of dimension at least

1 -1 K 1 -1 D 2 * n.
If the number of measurements m is fixed, then an 2 / 2 IO decoder must have a constant at least

1 1-K(1-m n ) .
In particular, if Σ-Σ contains an orthonormal basis of V , then K = 1 and the minimal number of measures to achieve NSP with constant D * is n/D 2 * . Similarly, if m is fixed so that m n, then a 2 / 2 instance optimal decoder has constant at least n m .

Examples

As discussed in the introduction, there is a wide range of standard models where Σ -Σ contains an orthonormal basis, and so where 2 / 2 IOP with dimensionality reduction is impossible. We provide here less trivial examples, where E = V is finite-dimensional.

Symmetric definite positive matrices with sparse inverse.

Lemma 2. Consider E is the space of symmetric n-dimensional matrices, and Σ ⊂ E the subset of symmetric positive-definite matrices with sparse inverse and with sparsity constant k ≥ n + 2 (note that k ≥ n is necessary for the matrix to be invertible). The set Σ -Σ contains an orthonormal basis of E.

Proof. This orthonormal basis we consider is made of the n(n + 1)/2 matrices: E i,i and

1 √ 2 (E i,j + E j,i ) i =j
, where E i,j is the matrix where the only nonzero entry is the (i, j) entry which has value 1.

First, consider B

i = I + E i,i , where I is the identity matrix. Since B -1 i = I -1 2 E i,i is n-sparse, we have B i ∈ Σ. Since, I ∈ Σ, we have E i,i = B i -I ∈ Σ -Σ.
Now, consider the matrix C i,j = 2I + E i,j + E j,i . This matrix is symmetric and for

x = (x 1 , . . . , x n ) ∈ R n , we have x T C i,j x = 2( x 2
2 -x i x j ) ≥ 0, so that C i,j is semi-definite positive. We can remark that C i,j is invertible and that its inverse is

1 2 I + 1 6 (E i,i + E j,j ) - 1 3 (E i,j + E j,i ), which is n + 2-sparse.
The fact that C i,j is invertible implies that it is definite, so that C i,j ∈ Σ. Therefore, we can write

E i,i + E j,j = C i,j -2I ∈ Σ -Σ.
Since Σ is a positive cone, multiplying this equality by 1 √ 2 yields the desired result.

Low-rank and nonsparse matrices. In (Candès et al., 2011a), the authors consider a matrix decomposition of the form L + S, where L is low-rank and S is sparse. In order to give meaning to this decomposition, one must avoid L to be sparse. To this end, a "nonsparsity model" for low-rank matrices was introduced.

Let E be the space of complex matrices of size n 1 × n 2 . Given µ ≥ 1 and r ≤ min(n 1 , n 2 ), let Σ µ,r be the set of matrices of E of rank ≤ r satisfying the two following conditions (denoting the SVD of such a matrix by r k=1 σ k u k v * k , where σ k > 0 and the u k and v k are unit-norm vectors) :

1. ∀k, u k ∞ ≤ µr n 1 and v k ∞ ≤ µr n 2 .
2. Denoting U and V the matrices obtained by concatenating the vectors

u k and v k , UV * ∞ ≤ µr n 1 n 2 .
These two conditions aim at "homogenizing" the entries of U and V. Note that we necessarily have µ ≥ 1.

Lemma 3. Let E = M n 1 ,n 2 (C) and Σ µ,r be the subset of E containing the matrices satisfying the two above conditions (with µ ≥ 1 and r ≥ 1). Then Σ µ,r -Σ µ,r contains an orthonormal basis.

Proof. Since Σ µ,r contains the null matrix, it is sufficient to prove that Σ µ,r contains an orthonormal basis. Let {e k } n 1 k=1 and {f } n 2 =1 be the discrete Fourier bases of C n 1 and C n 2 , that is

e k = 1 √ n 1 1, exp 2iπk n 1 , . . . , exp 2iπ(n 1 -1)k n 1 T f = 1 √ n 2 1, exp 2iπ n 2 , . . . , exp 2iπ(n 2 -1) n 2 T .
Then the n 1 n 2 rank-1 matrices of the form e k f * are elements of Σ µ,r since they obviously satisfy the two above conditions. But they also form an orthonormal basis of E, since each entry of e k f * is of module 1

√ n 1 n 2 and that, denoting •, • the Hermitian scalar product on E, we have

e k f * , e k f * = n 1 -1 u=0 exp 2iπu k -k n 1 n 2 -1 v=0 exp 2iπv - n 2 = δ k k δ ,
proving that these matrices form an orthonormal basis of E.

56CHAPTER 3. PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

3.6 The NSP and its relationship with the RIP

In the following section, we consider E = G (and A = I). As we have seen, one cannot expect to get 2 / 2 instance optimality in a dimensionality reduction context. This raises the following question: given pseudo-norms • G and • F defined respectively on E and F , is there a pseudo-norm • E such that IOP holds? We will see that this property is closely related to the RIP on M.

Generalized RIP and its necessity for robustness

The Restricted Isometry Property is a widely-used property on the operator M which yields nice stability and robustness results on the recovery of vectors from their compressive measurements. In the usual CS framework, the RIP provides a relation of the form (1

- δ) x G ≤ Mx F ≤ (1 + δ)
x G for any vector x in Σ 2k . The norms • G and • F are usually both taken as the 2 -norm. A form of RIP can easily be stated in a generalized framework: we will say that M satisfies the RIP on Σ -Σ if there exists positive constants

α, β such that ∀z ∈ Σ -Σ, α z G ≤ Mz F ≤ β z G . (3.30)
Similarly to the sparse case, it is possible to make a distinction between lower-RIP (left inequality) and upper-RIP (right inequality). Let's remark that this definition has been stated for vectors of Σ -Σ: this choice is justified by the links between this formulation and the NSP, which will be discussed later in this section. Let's also note that this form of RIP encompasses several generalized RIP previously proposed: the Ω-RIP [START_REF] Giryes | Greedy-Like Algorithms for the Cosparse Analysis Model[END_REF], the D-RIP (Candès et al., 2011) and the Union of Subspaces RIP [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF]). Let's now suppose the existence of decoders robust to noise, that is for all δ > 0, (3.19) is satisfied for a certain ∆ δ . This property implies the Robust NSP (3.20) with the same constants according to Theorem 5. By considering h ∈ Σ -Σ, the Robust NSP reads:

∀h ∈ Σ -Σ, h G ≤ D 2 Mh F . (3.31)
This is the lower-RIP on Σ -Σ, with constant 1/D 2 . The stability to noise therefore implies the lower-RIP on the set of differences of vectors of Σ, which is therefore necessary if one seeks the existence of a decoder robust to noise.

M -norm instance optimality with the RIP

The lower-RIP is necessary for the existence of a Robust instance optimal decoder, but what can we say this time if we suppose that M satisfies the lower-RIP on Σ -Σ with constant α, that is ∀z ∈ Σ -Σ, α z G ≤ Mz F ? We will prove that in both the noiseless and the noisy cases, this implies the IOP with norms • G and • M , the latter being called "M -norm" c and involving • G and • F . Let's define the M -norm on E as the following quantity, extending its definition for 2 norms in [START_REF] Peleg | Compressed sensing and best approximation from union of subspaces: Beyond dictionaries[END_REF] and its implicit appearance in the proof of early results of the field [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF]:

∀x ∈ E, x M = x G + 1 α Mx F . (3.32)
c to highlight its dependency on the Measurement operator Note that the term M -norm should be understood as M -pseudo-norm in the general case: if • F and • G satisfy the properties listed in Table 3.1, then • M satisfies the same properties as • G . However, when • G and • F are norms, • M is also a norm. We will note d M (., .) its associated (pseudo-)distance. The following theorem states that this • M allows one to derive an NSP from the lower-RIP on Σ -Σ.

Theorem 9. Let's suppose that M satisfies the lower-RIP on Σ -Σ with constant α (left inequality of (3.30)). Then the following Robust NSP is satisfied:

∀h ∈ E, h G ≤ d M (h, Σ -Σ) + 1 α Mh F . (3.33)
In particular, the following regular NSP is satisfied:

∀h ∈ N , h G ≤ d M (h, Σ -Σ). (3.34)
Therefore, if M satisfies the lower-RIP on Σ -Σ with constant α, then for all δ > 0, there exists a noise-robust instance optimal decoder ∆ δ satisfying the following property (Theorem 6): ∀x ∈ E, ∀e ∈ F,

x -∆ δ (Mx + e) G ≤ 2d M (x, Σ) + 2 α e F + δ. (3.35)
Note that in [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF], the author explored the implication of a lower-RIP on Σ -Σ for the case where Σ is an arbitrary UoS and . G / . F are the 2 norm. He proved that this generalized lower-RIP implies the following IOP: for all δ > 0, there exists a decoder ∆ δ such that ∀x ∈ E, ∀e ∈ F, ∀z ∈ Σ, .36) In this set-up, the instance optimality in equation (3.35) can be reformulated as:

x -∆ δ (Mx + e) 2 ≤ x -z 2 + 2 α M(x -z) + e 2 + δ. ( 3
∀x ∈ E, ∀e ∈ F, ∀z ∈ Σ, x -∆ δ (Mx + e) 2 ≤ 2 x -z 2 + 2 α M(x -z) 2 + 2 α e 2 + δ. (3.37)
Comparing these two instance optimality results, we can remark that the one in (Blumensath, 2011) is slightly tighter. This is merely a consequence of the difference in our method of proof, as we add the NSP as an intermediate result to prove instance optimality. The upper bound in [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF] can also be derived in our case with the same proof layout if we suppose the lower-RIP. Compared to [START_REF] Blumensath | Sampling and reconstructing signals from a union of linear subspaces[END_REF], our theory deals with general (pseudo-)norms and sets Σ beyond Union of Subspaces.

Upper-bound on the M -norm by an atomic norm

As we have seen, provided a lower-RIP on Σ -Σ, an NSP can be derived with the Mnorm as • E . However, this may look like a tautology since the M -norm explicitly depends on M. Hence, one may wonder if this NSP is of any use. We will prove in the following that provided an upper-RIP on a certain cone Σ (which can be taken as RΣ), a more natural upper bound can be derived by bounding the M -norm with an atomic norm [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF]. In particular, this type of inequality applied to the usual k-sparse vectors and low-rank matrices models give, under standard RIP conditions, instance optimality upper bounds with typical norms.

We will suppose in this section that • G is a norm.

58CHAPTER 3. PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

The atomic norm • Σ Let Σ be a subset of E and let E be the closure of span(Σ ) with respect to the norm • G . For x ∈ E , one can define the "norm" x Σ by:

x Σ = inf +∞ k=0 x k G | ∀k, x k ∈ RΣ and x - K k=0 x k G → K→+∞ 0 . (3.38)
Remark that there may be some vectors x for which x Σ = +∞, if +∞ k=0 x k G = +∞ for any decomposition of x as an infinite sum of elements of RΣ . However, the set [START_REF] Devore | Some remarks on greedy algorithms[END_REF]. In the following, we assume that V = E. Note that this norm can be linked to atomic norms defined in [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF] by considering A as the set of normalized elements of Σ with respect to • G . Now suppose M satisfies an upper-RIP on Σ , so that

V = {x ∈ E| x Σ < +∞} is a normed subspace of E which contains Σ (DeVore
∀x ∈ Σ , Mx F ≤ β x G . (3.39) For x ∈ E admitting a decomposition +∞ k=0 x k on RΣ , we can therefore upper bound Mx F by +∞ k=0 Mx k F ≤ β +∞ k=0 x k G .
This inequality is valid for any decomposition of x as a sum of elements of RΣ , so that Mx F ≤ β x Σ . Therefore, under these hypotheses,

∀x ∈ E, x M ≤ x G + β α x Σ ≤ 1 + β α x Σ . ( 3.40)
In particular, we have the following result:

Theorem 10. Suppose M satisfies the lower-RIP on Σ -Σ with constant α and the upper-RIP on Σ with constant β, that is

∀x ∈ Σ -Σ, α x G ≤ Mx F (3.41) and ∀x ∈ Σ, Mx F ≤ β x G . (3.42)
Then for all δ > 0, there exists a decoder ∆ δ satisfying: ∀x ∈ E, ∀e ∈ F,

x -∆ δ (Mx + e) G ≤ 2 1 + β α d Σ (x, Σ) + 2 α e E + δ, (3.43)
where d Σ is the distance associated to the norm • Σ .

Remark 3. Note that these results can be extended with relative ease to the case where

• G is not necessarily homogeneous but p-homogeneous, that is λx G = |λ| p x G .

Study of • Σ in two usual cases

We now provide a more thorough analysis of the norm • Σ for usual models which are sparse vectors and low-rank matrices. In particular, we give a simple equivalent of this norm involving usual norms in the case where • G = • 2 (for matrices, this is the Frobenius norm).

The norm • Σ relies on the decomposition of a vector as a sum of elements of RΣ. When Σ is the set of k-sparse vectors or the set or matrices of rank k, there are particular decompositions of this type:

• In the case where Σ is the set of k-sparse vectors, a vector x can be decomposed as ∞ j=1 x j , where all x j are k-sparse vectors with disjoint supports, which are eventually zero, and such that any entry of x j does not exceed any entry of x j-1 (in magnitude). This is a decomposition of x into disjoint supports of size k with a nonincreasing constraint on the coefficients.

• Similarly, in the case where Σ is the set of matrices of rank k and N is a matrix, the SVD of N gives a decomposition of the form N = ∞ j=1 N j , where the N j are rank k, eventually zero matrices such that any singular value of N j does not exceed any singular value of N j-1 .

For j ≥ 2, we can upper bound the quantity x j 2 using the assumption on the particular decomposition:

x j 2 ≤ √ k x j ∞ ≤ √ k x j-1 1 k = x j-1 1 √ k . Similarly, N j 2 ≤ N j-1 * √ k
, where • * is the trace norm, defined as the sum of singular values. We can therefore, in both cases, upper bound the norm • Σ . In the case of k-sparse vectors, this gives:

x Σ ≤ x 1 2 + j≥1 x j 1 √ k ≤ x 2 + x 1 √ k . (3.44)
In the case of matrices of rank k, this gives:

N Σ ≤ N 1 2 + j≥1 N j 1 √ k ≤ N F + N * √ k . ( 3.45)
We can also upper bound the right hand side of these equations by O( • Σ ) with a small constant, which will prove that the norms defined in these equations are of the same order. Indeed, a simple application of the triangle inequality gives us first that x 2 ≤ x Σ and N F ≤ N Σ . Then, considering a decomposition of x as a sum of k-sparse vectors j≥1 x j , we get

x 1 √ k ≤ j≥1 x j 1 √ k ≤ j≥1 x j 2 (3.46)
(indeed, each x j can be viewed as a k-dimensional vector and we have for such a vector

x j 1 ≤ √ k x j 2 ). Similarly, N * √ k ≤ j≥1 N j F . (3.47)
Since these upper bounds are satisfied for any decomposition, they can be replaced respectively by x Σ and N Σ . Finally, we have

x 2 + x 1 √ k ≤ 2 x Σ and N F + N * √ k ≤ 2 N Σ .
(3.48)

In these two cases, the norm • Σ is therefore equivalent (with constants bounded by 2) to the norm

• 2 + • 1 √ k (which is • F + • * √ k for matrices).
We have thus shown: Lemma 4. When Σ is the set of k-sparse vectors, the norm • Σ satisfies

• Σ ≤ • 2 + • 1 √ k ≤ 2 • Σ . (3.49) 60CHAPTER 3

. PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

When Σ is the set of rank-k matrices, it satisfies

• Σ ≤ • F + • * √ k ≤ 2 • Σ . (3.50)
We can therefore remark that for these two standard models, the norm • Σ can easily be upper bounded by usual norms under RIP conditions, yielding an IOP with a usual upper bound. We can also note that stronger RIP conditions can yield a stronger result: in [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF], the author proves that under upper and lower-RIP on Σ -Σ with Σ being the set of k-sparse vectors, an instance optimal decoder can be defined as the minimization of a convex objective: the 1 norm, which appears as strongly connected to the norm • Σ . One may then wonder if a generalization of such a result is possible: when can an instance optimal decoder be obtained by solving a convex minimization problem with a norm related to • Σ ?

Discussion and outlooks on Instance Optimality

Let's now review the results presented in this chapter and give some insights on interesting future work. As has been detailed throughout the chapter, Instance Optimality is a property presenting several benefits:

• It can be defined in a very general framework, for any signal space, signal model and pseudo-norms, as well as for both noiseless and noisy settings.

• It is a nice uniform formulation of the "good behavior" of a decoder and thus of the well-posedness of an inverse problem.

• It can be linked to Null Space Property and Restricted Isometry Property, which provide necessary and/or sufficient conditions for the existence of an instance optimal decoder.

We now present some immediate outlooks and interesting open questions related to Instance Optimality and to the results presented in this paper.

3.7.1 Condition for the well-posedness of the "optimal" decoder.

We have seen that for general models Σ, an additionnal term δ appears in the right hand side term of the instance optimality inequality ((3.13), (3.19)), reflecting the fact that the minimal distance of the "optimal" decoder (B.7) may not be reached at a specific point. However, as mentioned in Property 1, this additive constant can be dropped in the noiseless case provided Σ+N is a closed set. One can then wonder if there exists a similar condition (e.g., a sort of local compactness property) in the noisy case for which one can drop the constant δ and get a more usual instance optimality result.

Compressed graphical models.

As has been mentioned in Section 3.2.4, the case where Σ is the set of symmetric definite positive square matrices with sparse inverse is related to high-dimensional Gaussian graphical models. In Lemma 2, we showed this type of models fits in our theory since we could apply Theorem 8 in this case, proving the impossibility of 2 / 2 IOP in a dimensionreduction case. Yet, as for other signal models, can Gaussian graphical models satisfy some IOP/NSP with different norms in a compressive framework? 3.7.3 Guarantees for signal-space reconstructions and more.

When D is a redundant dictionary of size d×n and the signals of interest are vectors of the form z = Dx, where x is a sparse vector, traditional reconstruction guarantees from y = Mz assume the RIP on the matrix MD. This is often too restrictive: for example when D has strongly correlated columns, failure to identify x from y does not necessarily prevent one from correctly estimating z. Recent work on signal-space algorithms [START_REF] Davenport | Signal space cosamp for sparse recovery with redundant dictionaries[END_REF] has shown that the D-RIP assumption on M is in fact sufficient.

The framework presented in this paper offers two ways to approach this setting:

• Considering Σ = Σ k as the set of k-sparse vectors of dimension n and A = D, the upper bound on the reconstruction error is of the form d E (x, Σ k ). Signal-space guarantees can be envisioned by choosing a metric

• E = D • .
• Considering Σ = DΣ k as the set of d-dimensional vectors that have a k-sparse representation in the dictionary D and A = I, the upper bound is of the form d (z, DΣ k ).

In [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF], the authors propose a result similar to instance optimality by upper bounding, for a Total Variation decoder, the reconstruction error of an image X from compressive measurement by a quantity involving d 1 (∇X, Σ k ), where ∇ is the gradient operator, Σ k the k-sparse union of subspaces (in the gradient space) and d 1 is the 1 distance. This quantity is therefore the distance between the gradient of the image and the k-sparse vectors model. Can such a bound be interpreted in our framework, and possibly be generalized to other types of signals?

3.7.4 Task-oriented decoders versus general purpose decoders.

We already mentioned two very different application set-ups, in medical imaging and audio source separation, where only parts of the original signals need to be recovered. One can think of other, more dramatic, cases where only task-oriented linear features should be reconstructed. One such situation is met in current image classification work-flows. Indeed, most recent state-of-art image classification methods rely on very high-dimensional image representation (e.g., so called Fisher vectors, of dimension ranging from 10,000 to 200,000) and conduct supervised learning on such labeled signals by means of linear SVMs [START_REF] Sanchez | Image classification with the fisher vector: Theory and practice[END_REF]. Not only this approach yields top-ranking performance in terms of classification accuracy on challenging image classification benchmarks, but it also permits very large scale learning thanks to the low complexity of linear SVM training and its efficient implementations, e.g., with stochastic gradient descent. For each visual category to recognize, a linear classifier w is learned, which associates to an input image with representation x the score w T x. The single or multiple labels that are finally assigned to x by the system depend on the scores provided by all trained classifiers (typically from 10 to 100), hence on a vector of the form Ax, where each row of A is one one-vs-all linear SVM. In this set-up, the operator A implies a dramatic dimension reduction. For very large scale problems of this type, storing and manipulating original image signatures in the database can become intractable. The theoretical framework proposed in this paper might help designing new solutions to this problem in the future. In particular, it will provide tools to answer the following questions:

• A being given (learned on a labeled subset of the database): can one design a compressive measurement operator M such that the "classifiers" scores can be recovered 62CHAPTER

PERFORMANCE OF DECODERS IN LINEAR INVERSE PROBLEMS

directly from the compressed image signature Mx, hence avoiding the prior reconstruction of the high-dimensional signal x?

• M being given ("legacy" compressed storing of image signatures): what are the linear classifier collections that can be precisely emulated in the compressed domain thanks to a good decoder ∆?

Note that this classification-oriented set-up might call for a specific norm • G on the output of a linear score bank.

Another important domain of application that might benefit from both aspects (general purpose and task-oriented) of our work is data analysis under privacy constraints. Two scenarii can be envisioned, where our framework could help decide whether or not such constraints are compatible with the analysis of interest:

• General purpose scenario: given a linear measurement operator M of interest for further analysis, can one guarantee that there is no decoder permitting good enough recovery of original signals?

• Task-oriented scenario: the operator M serving as a means to obfuscate original signals such that critical information can't be recovered, let's consider a specific analysis task on original signals requiring the application of the feature extractor A.

Can this task be implemented on obfuscated signals instead, via a good decoder ∆, hence in a privacy-preserving fashion?

3.7.5 Worst case versus average case instance optimality.

The raw concept of Instance Optimality has a major drawback: the uniformity of the bound may impose, in some settings, a large global instance optimality constant whereas the inverse problem is well posed for the vast majority of signals. Let's consider the example depicted in Figure 3.5, where the signal space E is of dimension 2, the signal model Σ is a point cloud mostly concentrated along the line D and the measurement operator M is the orthogonal projection on D. The figure depicts the ratio (approximation error)/(distance to model) for each x ∈ R 2 . The optimal constant, which is the supremum of these ratios, is infinite: the ratio actually goes to infinity in the vicinity of the point p. However, for the vast majority of vectors, the ratio is rather low (the blue section covers most of the space).

An interesting outlook to circumvent this pessimistic "worst-case" phenomenon is to consider a probabilistic formulation of instance optimality, as in [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF]: given Ω a probability space with probability measure P , and considering M as a random variable on Ω, is there a decoder ∆(.|M) (which computes an estimate given the observation and the particular draw of the measurement operator M) such that for any x ∈ E, the instance optimality inequality

x -∆(Mx|M) G ≤ Cd E (x, Σ) (3.51) holds with high probability on the drawing of M? A particular challenge would be to understand in which dimension reduction scenarii there exists both a probability measure and a decoder with the above property. Another possible formulation of probabilistic instance optimality is to define a probability distribution on the signal space and to upper bound the average reconstruction error of the vectors, as in [START_REF] Yu | Statistical compressed sensing of gaussian mixture models[END_REF]. Chapter 4

Explicit embeddings for nearest neighbor search

As has been mentioned in Chapters 1 and 2 and pictured in Figure 2.1, individual dimension reduction of vectors can be used to reduce the cost of a learning task. In this chapter, we will consider a problem which is not itself a learning problem as the beginning of Chapter 2 depicts them a , but which can constitute a particular step in some global learning task. The envisioned task is nearest neighbor (NN) search. Let Ω be a set supplied with a distance d and X = {x 1 , . . . , x L } be a database of elements of Ω. Given a new element q ∈ Ω, typically called a query, a nearest neighbor of q in X is a solution of the problem argmin x∈X d(q, x).

(4.1)

Usually, this solution will be unique. Such a problem arises in applications such as image search [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF][START_REF] Torralba | Small codes and large databases for recognition[END_REF] and classification [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF][START_REF] Boiman | In defense of nearest neighbor based image classification[END_REF][START_REF] Deng | ImageNet: A Large-Scale Hierarchical Image Database[END_REF][START_REF] Perronnin | Large-scale image categorization with explicit data embedding[END_REF]. In such practical applications, one is usually interested in computing not only the nearest neighbor of a query but its k nearest neighbors. In this chapter, we will develop methods and perform experiments in a framework where we only seek the nearest neighbor of a query for simplicity, but the contributions described in this chapter can be easily extended to the case where one searches for several nearest neighbors of a query.

Approximate search. The most basic way to find a nearest neighbor is to compute all the distances d(q, x) for all x ∈ X and to find the minimal values. However, these computations come at a cost of O(LC d ), where C d is the cost of computing one distance value between two elements of Ω. This cost can be prohibitive when L is large, which motivates the search for less costly approximate nearest neighbor (ANN) search schemes. These schemes typically rely on the computation of small-dimensional signatures, which are compressed representations of the elements of Ω. The signatures will be outputs of an encoding function s : Ω → S, where the signature space S is provided with a distance which is much faster to compute than the distance of Ω.

The ANN search of a query q is therefore performed in this signature space. Note that as discussed in Section 1.2.1, an obvious way of using these signatures is to perform a full a That is a problem where one wants to infer some parameters from a training data set.

NN search in the signature space by computing all the distances between the signature of the query and the signatures of the database elements. In this chapter, we will consider methods applied in this "full-search" framework.

An ANN search pipeline will typically divide into offline steps, which will regroup the "fixed costs" and prepare the database X to the actual search step, and online steps, which are the steps where one gets a query q and needs to perform the search. The offline meta-procedure will usually be as follows:

• Choice of the encoding function: The encoding function s which will compute signatures is determined. In methods such as LSH, this step is generally quick since the function is a concatenation of randomly chosen hash functions. However, as we will see in Section 4.1, the choice of the signatures can be the result of a learning step involving the data, in order to choose adequate signatures.

• Computation of the signatures: The signatures s(x 1 ), . . . , s(x L ) are computed and stored. For a fixed database, this operation is only required to be performed once, since the same signatures can be used to perform as many NN searches as needed.

Given a query q for which one wants to perform an ANN search among the elements in X , the online meta-procedure will be as follows:

• Full NN search in the signature space: The signature s(q) is computed, and the distances between s(q) and every s(x r ) is computed. The N nearest neighbors in the signature space are identified, where N is a predetermined integer, supposedly far smaller than the size of X . This step may be costly since it linearly scales with the size of X , but the signatures are chosen so that it is far less costly than a full search in the space Ω.

As an alternative, in the case where S ⊂ Ω, one can perform the full NN search in an asymmetric way, that is by computing the distances d(q, s(x r )) instead of the distances d(s(q), s(x r )). This case, although non-intuitive at first glance, will be instantiated in this chapter via the Product Quantization (Jégou et al., 2011), in which the authors argue that an asymmetric scheme yields better precision than a symmetric scheme.

• Reranking: The N approximate nearest neighbors identified in the previous step are then reranked, that is the effective distances between q and the corresponding x i are computed. If the signature function s and the number N are well designed, the nearest neighbor of q is likely to be among these N elements.

In this case, if C d is the cost of computing a distance in Ω and C s is the cost of computing a distance in the signature space (one should have C s C d ), then the fullsearch cost in Ω is O(LC d ) whereas the cost of the described approximate scheme is O(LC s + NC d ). The approximate scheme will yield computational savings if

LC s + NC d LC d = C s C d + N L 1. (4.2)
In practice, these two fractions are 1, so that the approximate scheme yields substantial complexity reduction.

ANN search methods. Significant progress has been achieved on ANN search in the last decade, in particular with the substantial contribution of the LSH approach [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF], described in Chapter 1. Variants of this method have been developed, consisting of different types of low-dimensional signatures [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF][START_REF] Weiss | Spectral hashing[END_REF]Jégou et al., 2011). These methods typically allows ANN search systems to scale to datasets comprising 10 9 elements.

However, the aforementioned methods only consider simple metrics such as p norms, and mainly the Euclidean norm. In this chapter, we will consider the case where the distance d derives from a Mercer kernel, which is a broad generalization of the case where the distance is Euclidean. This case has less been studied in the ANN literature, although kernels provide a wide range of distances which are more elaborate than the Euclidean distance. Kernels are reviewed in the next paragraph.

Mercer kernels and nearest neighbors. Let's first recall standard definitions and properties of Mercer kernels [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization and Beyond[END_REF]. A positive semi-definite (PSD) kernel on a space Ω is an application K : Ω 2 → R such that for any collection X = {x 1 , . . . , x L } of elements of Ω, the Gram matrix

K(X ) =    K(x 1 , x 1 ) . . . K(x 1 , x L ) . . . . . . K(x L , x 1 ) . . . K(x L , x L )    (4.3)
is symmetric semi-definite positive. For any such kernel, there exists a unique (up to isomorphism) Hilbert space H and application Ψ : Ω → H such that ∀x, y ∈ Ω, K(x, y) = Ψ(x), Ψ(y) H , (4.4) where •, • H is the scalar product of H. H is usually named implicit space, since it is a space which allows the implicit interpretation of the action of K as a scalar product, provided one knows the implicit function Ψ which embeds Ω into H. Usually, one does not know this function but some interesting properties can still be derived, relying on the so-called kernel trick : one does not need to explicitly know the values of Φ(x) and Φ(y) to compute their scalar product in H, since this scalar product is equal to K(x, y).

In particular, K induces a distance d K on Ω by posing

d K (x, y) = K(x, x) + K(y, y) -2K(x, y). (4.5)
Let's notice that if we define

S K (1) = {x ∈ Ω : K(x, x) = 1}, (4.6)
then for all x, y ∈ S K (1),

d K (x, y) = 2 (1 -K(x, y)) . (4.7)
This property ensures that when considering data X ⊂ S K (1) and a query q ∈ S K (1), the problem (4.1) with d = d K is equivalent to the following problem:

argmax x∈X K(q, x), (4.8)
that is searching for the nearest neighbor in the sense of d K is equivalent to searching for the most correlated vector in the sense of K. This constitutes a good incentive to normalize data, that is project it onto S K (1), when it is possible: it allows one to consider kernels instead of distances. This will be used in our contribution, where we will perform our experiments on normalized data.

In the particular case where Ω is a compact Hausdorff space provided with a finite measure µ with support Ω, the kernel K is a Mercer kernel [START_REF] König | Eigenvalues of compact operators with applications to integral operators[END_REF], that is H = 2 (R) and for all i ∈ N, there exists λ i ∈ R and Ψ i : Ω → R satisfying

Ψ(x) = λ i Ψ i (x) ∞ i=0 , ( 4.9)
with λ i ↓ 0 and Ω Ψ i Ψ j dµ = δ j i . The λ i 's and Ψ i 's are respectively the eigenvalues and the eigenfunctions of K. In this case, the embedding Ψ therefore transforms elements of Ω into square-summable sequences with decreasing average energy. This is particularly important in order to approximate the kernel values, as we will review when describing our contribution.

Quick overview of the contribution. In this chapter, we propose two methods for NN search when the considered distance is derived from a Mercer kernel K. The described methods rely on explicit embeddings, that is approximating the unknown function Ψ by an explicit function Ψ which maps Ω to the Euclidean space R E (E > 0 is the embedding dimension), and satisfying Ψ(q), Ψ(x i ) H ≈ Ψ(q), Ψ(x i ) , (4.10) where •, • is the scalar product in R E . We will underline two benefits from exploiting explicit embeddings when designing ANN search methods with kernel distances:

• Exact search scheme: As shown in our first proposed method, one can derive error bounds from an explicit embedding such as KPCA, which can be used to derive an exact search scheme while performing the NN search in R E instead of performing it in Ω. This can yield slightly lower search times while still be ensured to obtain the nearest neighbor of the query.

• Approximate search scheme: More importantly, our second proposed method aims at computing low-dimensional signatures for efficient large-scale ANN search in two steps: first the data X is mapped from Ω to R E with an explicit embedding Ψ, then small-dimensional signatures are computed using methods adapted to an Euclidean space, and which are therefore not applicable directly in Ω. This two-step procedure is experimentally shown to perform better than previously proposed methods which aim at finding signatures directly in the implicit space, and are detailed in the following section.

State of the art: ANN search for kernels and explicit embeddings

In this section, we first describe two previously existing methods developed in order to compute signatures specifically aimed at solving problem (4.1) when d is a Mercer kernel distance. The two presented methods both produce binary signatures, which are then compared with respect to the Hamming distance in the signature space. Note that other methods aimed at producing small-size signatures for kernels exist, such as (He et al., 4.1. STATE OF THE ART: ANN SEARCH FOR KERNELS AND EXPLICIT EMBEDDINGS69 2010;[START_REF] Gorisse | Locality-sensitive hashing for chi2 distance. Pattern Analysis and Machine Intelligence[END_REF], but were not considered in our work because the framework presented in these papers is different from ours and closer to a classification context. Then is reviewed the main explicit embedding method: the Kernel Principal Component Analysis (KPCA), which will be used in our approach. We also briefly mention several other explicit embeddings procedures, which may also be used in the framework we consider, provided they are applied to the proper kernels.

Kernel ANN methods

Kernelized Locality Sensitive Hashing

In LSH, a typical choice for the hash functions for n-dimensional vectors are functions of the type h : x → sign( x, r ), where r is usually chosen as a random draw with respect to a normal distribution N (0, I). Kernelized Locality Sensitive Hashing (KLSH) [START_REF] Kulis | Kernelized locality-sensitive hashing for scalable image search[END_REF] aims at imitating this Gaussian random choice in the implicit space H.

To this end, the authors select M data points y 1 , . . . , y M among X and approximate the space H by the finite-dimensional space F spanned by Ψ(y 1 ), . . . , Ψ(y M ). The hash functions are determined by randomly choosing subsets of Q < M points among the Ψ(y i ) and taking the empirical mean of these Q points. The image of a vector of Ω by such a hash function is computable thanks to the kernel trick. After a proper rescaling aimed at correcting variance deviations, this random choice yields a vector of F which is approximately drawn with respect to an isotropic Gaussian distribution N (0, I), according to the Central Limit Therorem.

Moreover, since the considered directions belong to the space F spanned by elements of X , the image by the corresponding hash functions of any vector Ψ(x) with x ∈ Ω can explicitly be computable since the calculation only requires to compute kernel values between x and the y i .

Random Maximum Margin Hashing

Random Maximum Margin Hashing (RMMH) [START_REF] Joly | Random maximum margin hashing[END_REF] aims at finding hash functions of the form h : x → sign( Ψ(x), r H + b), where r is once again a vector of H and b is a real number. The authors express the problem of finding such hash functions as an SVM problem.

Each hash function is chosen as follows: an even number M of vectors y 1 , . . . , y M are randomly drawn from X and half of these vectors are labeled with label +1, the other half with label -1. The corresponding hash function is given by (4.11) where i is the label of y i and (α i ) M i=1 and b maximize the usal C-SVM cost:

h(x) = sign M i=1 α i i K(x, y i ) + b ,
1 2 M i,j=1 i j α i α j K(y i , y j ) + C M i=1 ξ i , (4.12) subject to ∀i, ξ i ≥ 0 and 1 -ξ i ≤ i   M j=1 j α j K(y i , y j ) + b   . (4.13)
Intuitively, this formulation finds hash functions which prevent close neighbors from having a different image by the hash function, since the hyperplane is chosen by margin maximization. This method has been presented as yielding better precision than KLSH for approximate search.

Explicit embeddings

In this section, we provide a quick review of the standard KPCA method, which can be considered as the standard explicit embedding procedure. We also give a brief list of other embedding methods which have been proposed.

Kernel Principal Component Analysis

KPCA was introduced in [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] as a way to perform explicitly PCA in the implicit space H, thus being a generic way to compute a finite-dimensional embedding of the data satisfying the condition (4.10).

Performing PCA on the set of vectors Ψ(X ) = (Ψ(x i )) M i=1 in implicit space H consists in projecting them onto the E-dimensional subspace V that minimizes the following meansquared error:

M i,j=1 (K(x i , x j ) -P W Ψ(x i ), P W Ψ(x j ) H ) 2 (4.14)
over all E-dimensional subspaces W (P W being the orthogonal projection onto W ). It is possible to compute explicitly this projection for any vector of the form Ψ(x) by exploiting the Gram matrix K(X ) = (K(x i , x j )) i,j=1,...M . Let us denote its E largest eigenvalues by (σ 2 i ) E i=1 and by (u i ) E i=1 the corresponding eigenvectors. Let K(x, :) be the vector [K(x, x 1 ), . . . , K(x, x M )] T . Then if (v i ) E i=1 denote the E first eigenvectors of the covariance operator of Ψ(X ), one has [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]: .15) This result allows us to define by PCA approximation an explicit embedding Ψ : Ω

ϕ i = P V (Ψ(x)), v i H = 1 σ i K(x, :), u i . ( 4
→ R E such that Ψ(x) = [ϕ i ] E
i=1 . This result allows us to define by PCA approximation an explicit embedding Ψ(x) = P V (Ψ(x)) and an approximate kernel K(x, y) = Ψ(x), Ψ(y) . We therefore have

K(x, y) = P V (Ψ(x)), P V (Ψ(y)) H . (4.16)
Note that KPCA is adaptive, that is the explicit embedding Φ it produces is learned on training data X , and thus depend on it.

Other explicit embeddings

Other embedding methods have been proposed for specific kernels in a classification context. In [START_REF] Maji | Max-margin additive models for detection[END_REF], the authors describe an explicit embedding for the histogram intersection kernel, defined for x, y ∈ R n by K(x, y) = n i=1 min(x i , y i ).

(4.17) In (Vedaldi andZisserman, 2010, 2012), the authors propose a nonadaptive (thus requiring no input data nor learning step) embedding algorithm for additive or multiplicative kernels, that is kernels defined on

R n as K(x, y) = n i=1 k(x i , y i ) or K(x, y) = n i=1 k(x i , y i ),
where k is a kernel on R.

In [START_REF] Perronnin | Large-scale image categorization with explicit data embedding[END_REF], the authors suggest applying KPCA independently in each dimension if the kernel K is additive. This yields a computational gain at the cost of a loss of precision in the embedding approximation.

As we will see, the methods we propose heavily relies on explicit embeddings. In all our experiments, we used KPCA as the embedding method, mainly because we are aiming at a generic kernel-independent method. In that sense, the KPCA has the benefit that it does not make any assumption on K other than the fact it is a Mercer kernel. Moreover, KPCA allows one to derive an embedding in R E without constraints on E, whereas for instance, [START_REF] Maji | Max-margin additive models for detection[END_REF][START_REF] Vedaldi | Efficient additive kernels via explicit feature maps[END_REF][START_REF] Perronnin | Large-scale image categorization with explicit data embedding[END_REF] require the output dimensionality to be greater than or equal to that of the initial space. Finally, KPCA yields in practice good approximations on real datasets.

Layout of the chapter

Our contribution consists in two methods relying on explicit embedding for ANN search. Since the description of each method is directly followed by practical experiments, the datasets used in the experiments are first described in Section 4.3.

The first proposed method is an exact search method and is described in Section 4.4, along with an analysis of its benefits and an experimental illustration.

The second proposed method is an approximate scheme and is described in Section 4.5. It is experimentally compared to the previously mentioned KLSH and RMMH schemes in terms of precision.

Datasets and experimental protocol

For the sake of reproducibility, we have used publicly available datasets of vectors that are large enough to evaluate the interest of NN search techniques. Table 4.1 gives a brief description of these benchmarks, which all consist of image descriptors extracted from real images.

In our case, each of these datasets is separated into three disjoint sets: a query set, a training set and a validation set, this separation being random for SIFT1M and BIGANN which do not have such explicit prior separation of the full dataset.

• SIFT1M : This benchmark introduced in (Jégou et al., 2011) consists of one million local SIFT descriptors [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] computed on patches extracted with the Hessian-Affine detector [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF].

• BIGANN (Jégou et al., 2011): This set provides a NN ground-truth for subsets of increasing size, from one million to one billion vectors. These descriptors were extracted using the Difference of Gaussians (DoG) detector and are therefore not equivalent to those of the SIFT1M benchmark.

• Imagenet [START_REF] Deng | ImageNet: A Large-Scale Hierarchical Image Database[END_REF]: This set contains the pre-computed 1,000-dimensional bag-of-words (BOW) vectors provided with the images of the 2010 Large-Scale Visual Recognition (ILSVRC'2010). For our search problem, the larger "train" set is used as the database to be searched, the "val" set to learn the parameters, and we randomly selected a subset of 1,000 queries from the "test" set to act as queries to be processed by our methods.

For these three datasets, we are looking for the most similar vectors with respect to the χ 2 kernel, as done in related works [START_REF] Joly | Random maximum margin hashing[END_REF][START_REF] Vedaldi | Efficient additive kernels via explicit feature maps[END_REF]. The χ 2 kernel is defined on Ω = R n by

K(x, y) = K χ 2 (x, y) = n i=1 x i y i x i + y i . ( 4

.18)

In our experiments, all the data vectors were 1 -normalized. This is the normalization for which K χ 2 (x, x) = 1, and thus it is the required normalization to get the equivalence between nearest neighbor and highest kernel value, as discussed at the beginning of the chapter.

In the experiments we describe, the search performance is measured by the Recall@R, which measures the rate of queries for which the nearest neighbor is present among the R closest neighbors in the signature/embedded space.

Exact search method

The first method we propose relies on an explicit embedding in order to cast the search problem in a Euclidean space. The benefit is that in such a space, the scalar products are matrix multiplications and such computations can be performed by heavily optimized routines.

In order to derive an exact NN search scheme, one needs to be able to retain every potential nearest neighbor, that is upper bounding the precision of the explicit embedding. Since we are using KPCA, the next section is dedicated on the derivation of a standard error bound for KPCA.

Error bound for the KPCA embedding

Precision bounds on the KPCA provided in the literature are usually probabilistic bounds on the error between the empirical eigenvalues/eigenfunctions computed from the matrix K(X ) (denoted σ 2 i and u i in Section 4.1.2) and their continuous counterparts, defined in the preliminary section of this chapter. Such bounds are difficult to obtain and they are not tight or require oracle knowledge about the implicit embedding [START_REF] Shawe-Taylor | On the eigenspectrum of the gram matrix and the generalization error of kernel pca[END_REF][START_REF] Braun | Accurate error bounds for the eigenvalues of the kernel matrix[END_REF].

However, it is possible to get tighter precision bounds for our problem by exploiting the data we consider. As shown in Section 4.1.2, KPCA is an orthogonal projection onto a subspace V in the implicit space H. The approximation error (x, y) = K(x, y) -K(x, y) can thus be written as a scalar product in the implicit space: (4.19) where P V ⊥ is the orthogonal projection on the subspace orthogonal to V . Let us denote .20) This quantity can be easily computed in practice. The Cauchy-Schwarz inequality provides a bound on (x, y) in terms of R(x) and R(y):

(x, y) = K(x, y) -K(x, y) = Ψ(x), Ψ(y) H -P V Ψ(x), P V Ψ(y) H = P V ⊥ Ψ(x), P V ⊥ Ψ(y) H ,
R(x) = ||P V ⊥ Ψ(x)|| H = K(x, x) -|| Ψ(x)|| 2 1 2 . ( 4
| (x, y)| ≤ R(x)R(y). ( 4

.21)

This bound can be used for nearest neighbor search. Despite its simplicity, we believe it has never been used in such a context. Let us assume that we look for the nearest neighbor of a query q in the vector dataset X and let us define

x 1 = argmax x∈X K(q, x) and x1 = argmax x∈X K(q, x), (4.22)
that is respectively the nearest neighbor and approximate nearest neighbor of the query q. We can get a lower bound for K(q, x 1 ) by noticing:

K(q, x 1 ) ≥ K(q, x1 ) ≥ K(q, x1 ) -R(q)R(x 1 ). (4.23)
For any vector x ∈ X , we can get an upper bound on K(q, x):

K(q, x) ≤ K(q, x) + R(q)R(x). ( 4

.24)

Combining (4.23) and ( 4.24), we obtain that, if for some x we have:

K(q, x) + R(q)R(x) < K(q, x1 ) -R(q)R(x 1 ), (4.25)
then K(q, x) < K(q, x 1 ) and x cannot be the nearest neighbor of q. The following section describes a procedure that uses this result to avoid kernel computations while still being sure to find the nearest neighbor of a query.

Exact search procedure: description and illustration

Assuming that the quantities R(x) are computed and stored while performing the embedding of all vectors in dataset X , the nearest neighbor of the query q is retrieved exactly by using the following procedure, which principle is illustrated in Figure 4.1 (left):

1. Compute Ψ(q) and R(q) using Equation (4.15).

2. Compute K(q, x) = Ψ(q), Ψ(x) for all x ∈ X .

Find x1 = argmax

x∈X K(q, x)

4. Find all elements y ∈ X such that: K(q, x1 ) -R(q)R(x 1 ) ≤ K(q, y) + R(q)R(y).

Let us denote by Y this set of vectors.

5. The true nearest neighbor of the query is x 1 = argmax y∈Y K(q, y).

Complexity Analysis and Experiments

Let L be the size of the base X , E be the embedding dimension, M be the number of vectors used for KPCA, N be the cardinality of Y, C K be the cost of computing the kernel between two vectors in the initial space and C K be the cost of computing the approximate kernel between two embedded vectors (which is the cost of a scalar product in the embedded space). The cost of the exact search in the initial space has complexity O(C K L).

The exact search in the embedded space as described above has complexity O(MC K ) for step 1 (scalar products are negligible if E M ), O(NC K ) for steps 2 to 4 (steps 3 and 4 have constant complexity per data vector, which is negligible compared to C K for high-dimensional data), and O(NC K ) for step 5.

Hence this method is guaranteed to find the nearest neighbor of q with complexity

O (M + N )C K + LC K . (4.26) If C K C K (
which is likely to be the case, especially for complicated kernels), and if M + N L, we can be sure to retrieve the nearest neighbor while performing much fewer calculations.

Figure 4.1 (right) illustrates this method on 10 6 SIFT descriptors taken from the dataset BIGANN and shows that it reduces the number of kernel values computed. Indeed, for M = 2, 048, 90% of the tested queries require less than 15,000 kernel computations, which is a comparison with less than 0.15% of the base with respect to the distance we are interested in. These comparisons are of negligible complexity, the exact search process in the embedded space thus roughly requires 10 6 scalar products in dimension 128 instead of 10 6 χ 2 kernel computations in dimension 128 for the exact search in the initial space, which leads to a speedup considering that scalar products have more optimized implementations than the χ 2 kernel.

A tradeoff has to be found on M : a high value will tighten the Cauchy-Schwarz bound but will result in more kernel computations to perform the embedding on the query. In this example, M = 2, 048 performs on average better than M = 8, 192 (8,000 average kernel computations vs. more than 11,000). See also Figures 4.2 and 4.3 for an analysis of the impact of parameters M and E on the search quality (independently of the number of kernel computations).

On the other hand, this approach does not perform well on the Imagenet database. For this dataset, the Cauchy-Schwarz inequality is too pessimistic to provide a good selection of the vectors to be reranked. Indeed, experiments showed that the number of reranked vectors reprensented 30% of the database. In this case, it is still possible to select fewer vectors by relaxing the Cauchy-Schwarz inequality and considering that αR(x)R(y) ≤ (x, y) ≤ βR(x)R(y), (4.27) where α and β are simply lower and upper bounds for cos(x, y). By empirically estimating values for α and β that work for many vectors, one can use concentration of measure inequalities to derive theoretical guarantees learned from data. However, one is not able then to precisely ensure the retrieval of the true nearest neighbor of a vector.

Approximate search method

Even though the previous method allows one to perform an exact search at a slightly reduced cost, it may still be too costly to perform the approximate search in the embedded space. In this case, it is preferable to compute very small-dimensional signatures from the database elements, as described in the preliminary section of this chapter.

The binary schemes described in Section 4.1.1 are implicit schemes: they aim at finding the encoding function directly as a hashing function in the implicit space H, expressing it as a function of the kernel K. Instead, the approximate scheme we propose relies on two separate steps:

1. The data is embedded in a Euclidean space R E using a function Ψ corresponding to KPCA b .

2. The signatures are computed by applying an Euclidean signature method on the Ψ(x i ) in R E .

This scheme is illustrated in Figure 4.4. Note that at step 2, if the data is normalized with respect to the kernel, any compression technique which aims at approximating the dot-product, the Euclidean distance or the cosine similarity could be employed. This includes among others LSH [START_REF] Charikar | Similarity estimation techniques from rounding algorithms[END_REF]Datar et al., 2004a) which approximates the cosine or Spectral Hashing [START_REF] Weiss | Spectral hashing[END_REF] and Locality Sensitive Binary Code [START_REF] Raginsky | Locality-sensitive binary codes from shift-invariant kernels[END_REF] which approximate the Euclidean distance. In this section, we particularly focus on a recently proposed technique known as Product Quantization (PQ) which has been shown to be state-of-the-art for the Euclidean distance (Jégou et al., 2011). It is discussed briefly in the next section.

Product quantization and variations

Most of the techniques based on short signatures for ANN search consider binary codes [START_REF] Kulis | Kernelized locality-sensitive hashing for scalable image search[END_REF][START_REF] Joly | Random maximum margin hashing[END_REF][START_REF] Gong | Iterative quantization: A procrustean approach to learning binary codes[END_REF]. In (Jégou et al., 2011), the authors proposed instead the PQ algorithm which produces non-binary codes allowing the estimation of Euclidean distances in the signature domain. For this purpose, a product quantizer is used to decompose the feature space into a Cartesian product of D subspaces. A given database vector x is decomposed into D subvectors of equal lengths, each of which is quantized separately using a k-means quantizer with a limited number of centroids. This number is typically set to 256 such that each subvector is quantized b Any other explicit embedding method can be used if applicable, but we used KPCA in our experiments. .5: Illustration of the PQ procedure on one vector. The vector is considered as a product of D subvectors of equal length, and each subvector is quantized separately using a k-means algorithm on the database X .

as a 1-byte code (8 bits) and therefore the vector is represented by a D-bytes code. This quantization is illustrated in Figure 4.5. The estimation of d(q, x) is done in an asymmetrical manner, i.e., the query is not quantized to prevent it from being approximated. Computing a distance estimate requires D table look-ups and additions, which remains competitive compared with operations on binary vectors at a given signature size.

Experiments on approximate search

In this section, we provide and discuss experimental results of our approximate search scheme compared to the other previously described binary methods. Let's first discuss a particular pre-processing which has been chosen for the PQ technique in our experiments.

Pre-processing for energy balance.

Since we are applying PQ after a KPCA, the energy is not equally distributed in the embedded vectors Ψ(x i ): the first components have high energy while the last components have smaller values. This has a negative impact on PQ since the last subvectors will have less variance than the first while the k-means clustering still produces the same number of centroids for each subvector. Intuitively, one should find a way to balance energy between the subvectors so that all the quantizers have a subvector of comparable energy to encode.

In [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF], the authors faced the same phenomenon after a regular PCA. In this context, they argue that the PQ technique is more effective if a rotation matrix is applied on the vectors produced by the PCA and prior to PQ encoding. The drawback is that the individual product quantizers partially encode the same information, since the decorrelation performed by PCA is lost.

In our contribution, we consider the use of a random permutation of the components of the vectors Ψ(x i ) as an alternative to applying a random rotation. The random permutation is cheaper to compute and it ensures that the components remain uncorrelated.

In order to choose the appropriate pre-processing for PQ, let's study the impact of three pre-processing techniques:

(i) No pre-processing is done at the output of the KPCA.

(ii) A random rotation is applied to the Ψ(x i ).

(iii) The components of the Ψ(x i ) are randomly permuted. Table 4.2 shows that the choice of the pre-processing stage has an impact on the search quality in terms of Recall@R. On SIFT1M, the random permutation is clearly the best choice. On Imagenet, the difference is less clear but the random permutation still yields better results. We will therefore adopt this choice in all subsequent experiments. Overall, the improvement of KPCA+PQ with respect to KLSH and RMMH can be considerable, especially for small codes. For instance, for SIFT1M and for B = 64, the recall@100 is about 80% for KPCA+PQ while it is 20% for KLSH. For Imagenet and for the same number B of bits, the recall@1000 of RMMH is somewhat above 30% while KPCA+PQ achieves on the order of 70%.

Large scale experiments.

Conclusion

This study shows that when dealing with a distance derived from a Mercer kernel, the use of explicit embeddings can benefit nearest-neighbor search both for exact and more generally for approximate search.

The proposed exact search algorithm based on the Cauchy-Schwarz inequality works very well on the SIFT descriptors and allows to reduce the cost of the exact search roughly to the cost of computing the scalar products in the embedded space. However, it is not applicable to the Imagenet high-dimensional descriptors due to the rough precision of the Cauchy-Schwarz inequality in a high-dimensional space. Nonetheless, it may be possible in this context to sacrifice the exact search for a quasi-exact search with probability bounds by learning the discrepancy between Cauchy-Schwarz inequality and the actual value of the approximation error for a range of training vectors.

The proposed approximate scheme, despite its simplicity, performs better than other encoding methods in this framework. Performing KPCA followed by a standard binary LSH encoding yields better precision than standard methods which directly mimic binary LSH in the implicit space, thus giving a first incentive to use explicit embeddings. Furthermore, mapping the initial elements into an Euclidean space allows one to use a wide range of schemes designed at computing signatures, and especially Product Quantization, which performs way better than binary schemes at a given signature size. Therefore, using an explicit embedding prior to an efficient Euclidean coding scheme should be considered as a baseline when measuring the performance of ANN search methods with small signatures in a kernelized framework.

Chapter 5

Conclusion

In this short chapter, we first propose to wrap up the principal contributions presented in this manuscript and to draw additional links between them, deriving some lessons from these different works. We then suggest some global research directions which seem interesting to pursue in order to benefit from the most conceptual part of this thesis work.

Summary of the contributions

The proof of concept of Chapter 2 focuses on estimating the parameters of a mixture of densities, instantiating it on isotropic Gaussian mixture models. Even though these are a simple problem and model, the estimation of such Gaussian mixture parameters acts as a step in many practical learning or data processing mechanisms. Therefore, this proof of concept is meaningful and acts as an incentive to come up with extensions of the method to more general Gaussian mixture models.

Aside from the experimental results, the originality of the approach itself motivates the search for alternate algorithms designed to perform learning with more data and less work [START_REF] Shalev-Shwartz | Using more data to speed-up training time[END_REF]. The standard algorithmic procedure, updating in an iterative way some parameters using all available data, is replaced in our approach by a two-step procedure in which the data is compressed to a sketch of fixed-size representation prior to the estimation step. The way this sketching step is built may heavily depend on the task one wishes to perform afterwards: one could either consider building a sketch particularly fitted to a specific parameter estimation or more generally building a (probably higher-dimensional) sketch aimed at different estimation or processing tasks. In all cases, this sketching procedure can be viewed as a "storage step", where the only amount of information that is kept on the data is the amount one needs to be able to exploit it.

Since the estimation procedure we proposed can be interpreted as a generalized compressed sensing problem, it seems natural to be inspired by the rich signal processing literature on the subject. In particular, we showed in Chapter 3 that standard theoretical results can be generalized to a broad class of models. The assumptions we made to derive our results are very mild, so that the studied guarantees are potentially applicable to the study of many linear compression schemes for learning. More particularly, the stated theorems give tools aimed at studying both the performance one will not be able to achieve and the performance one can achieve if a generalized RIP is satisfied. This gives an angle for tackling the theoretical study of linear sketching tools: one may aim at finding frameworks where these tools satisfy the RIP in order to get the developed instance optimality guarantees.

The work in Chapter 4 shows that a simple scheme for nearest neighbor search exploiting a well-used technique (namely, the KPCA) can still achieve better performance than more elaborated techniques. One may see there an incentive to make the most of existing techniques.

Perspectives

Some short-term perspectives have already been presented at the end of Chapters 2 and 3. In this section, we will elaborate on some of them and evoke more general lines of research (possibly very conceptual in their current form) inspired by the work presented in this manuscript. These perspectives mainly focus on compressed density estimation (and more generally compressed learning), linear inverse problems for general models and algorithmic methods to solve such problems.

Short-term perspectives

Robustness analysis of compressed Gaussian mixture estimation. In Chapter 2, we proved a theorem stating the injectivity of the sketching operator on the set of sparse mixtures of isotropic Gaussians, provided the frequencies were (deterministically) wellchosen. However, the operators we used in our practical experiments were randomly drawn, and the method was experimentally shown to work with far fewer measurements than stated in the theorem. In order to better estimate the number of required measurements needed to reconstruct accurately the mixture, we certainly need, as in compressed sensing theory, to consider a probabilistic framework where M is randomly drawn. This should allow us to find a number of measurements so that with high probability on the draw of M, a robustness result is satisfies (such as a RIP, so that the results of Chapter 3 are satisfied). In particular, it seems interesting to consider a method analog to [START_REF] Baraniuk | A simple proof of the restricted isometry property for random matrices[END_REF] by studying the average conservation of the 2 norm between f and Mf for f ∈ Σ + k,n , then using concentration of measure and union bound arguments to obtain a RIP-like result.

Reduction of the complexity of the compressed estimation algorithm. The compressed approach we developed in Chapter 2 allows memory savings in the case of numerous data, but the complexity of the algorithm is still high due to the numerous optimization steps required by the estimation. It is therefore important to overcome these complexity limitations in the case of numerous parameters to estimate. The bottleneck of the algorithm is the last step of the iteration, which performs a gradient descent on the objective function taken as a function of the parameters. To reduce the computational cost of this step, one may consider trying not to perform it at each iteration, or update in this fashion only a fraction of the parameters. The complexity will also be reduced if the algorithm is modified so that it converges faster: in that sense, it may be interesting to consider a hierarchical procedure where only low-frequency components of the sketch are used at the beginning of the optimization, leading to a rough estimation of the parameters, then high-frequency components are considered to refine the estimation.

Practical extension of the compressed framework to a broader class of Gaussians. As has been mentioned in the summary of the contributions, Gaussian models are widely used in learning frameworks because of their simplicity. It is therefore an interesting outlook to extend the proof of concept of Chapter 2 to richer families of Gaussians than isotropic Gaussians. In particular, Gaussians with diagonal covariance matrices are often used since they allow more variability than isotropic Gaussians without changing the order of magnitude of the number of parameters needed to describe them, which is O(kn) for k Gaussians in dimension n. An experimental study of the compressed framework instantiated on this enriched class of Gaussians is therefore motivated by potentially important applications. Moreover, the practical challenges can help finding a better way to build sketching operators and design algorithms to solve the compressed estimation problem.

Theoretical analysis of the difference between explicit and implicit hashing. We have experimentally shown in Chapter 4 that performing an explicit embedding then applying a Euclidean hashing scheme yielded better precision than using an implicit embedding which directly tries to mimic a hashing technique in the implicit space and use the kernel trick to define the hashing functions. It may be insightful to perform a theoretical analysis of this behavior. For instance, under assumptions on the probability distribution of the data and the decrease of the eigenvalues of K, can theoretical results be exploited to compare the expected precision of the different methods?

Mid-term perspectives

Generalized theoretical recovery results for compressive density mixture estimation. The compressive framework of Chapter 2 could be refined in the case where the mixture model is generic. In particular, a crucial question is the existence of key properties on the probability family P and the linear sketching operator M such that there exist theoretical recovery guarantees of a mixture of Σ + k (P) given its image by M. More precisely, finding conditions ensuring the generalized RIP for M stated in Chapter 3 would yield an IOP for M. If a probabilistic analysis of the robustness of the sketching operator has already been performed in the case of isotropic Gaussians, it may serve as a starting point to identify the prerequisites on the density families allowing to generalize the robustness results.

Characterization of the existence of "well-behaved" decoders. In Chapter 3, we studied the existence of instance optimal decoders without constraints on these decoders. The theorems in Section 3.4 provide decoders which are minimizers of possibly nonconvex and irregular functions. In order to link the theory to practical application, we need to identify the cases when it is possible to regularize these decoders (in particular, make them minimizers of convex objectives) and still getting instance optimality. A starting point can be the study of the difference between the decoders mentioned in the proofs of Section 3.4 theorems and the 1 decoder used in classical compressed sensing. Can the method used to prove the 1 ⇔ 0 minimization equivalence for sparse vectors be adapted to more general norms and models?

Long-term perspectives

Randomized methods for non-convex optimization and sparse decomposition. In order to solve the non-convex problem (2.19), we designed an algorithm based on random initialization in the parameter space to search for good candidates to add to the current estimate. Such methods, aimed at reconstructing a sparse signal on a very coherent (and even continuous) dictionary, may be studied more thoroughly and improved. An important theoretical question is the existence of probabilistic recovery guarantees for such algorithms, which would certainly need to identify some particular required properties on the random step for the algorithm to be robust. Given a dictionary D, is there an "optimal" parametrization of D which essentially makes decomposition algorithms more robust?

The enriched tradeoffs of large-scale learning. In [START_REF] Bottou | The tradeoffs of large scale learning[END_REF], the authors studied the impact of modifying the expression of learning objectives to sacrifice optimization precision for computational savings. In particular, they proved that at a given target training time, the required time to reach an asymptotic generalization error can be reduced by sacrificing the precision of the iteration steps for computational gains. However, they do not take account the memory used by the training algorithm. Moreover, they consider algorithms aimed at minimizing a single objective.

A deep question is therefore the tradeoff between generalization error, required memory and computational time in a learning procedure. It would be particularly appealing to come with the relationship between memory and computational time at a target precision for a class of algorithms which compute low-dimensional representations of training data as a step before the actual estimation of the training parameters. The consideration of computational tradeoffs for a class of algorithms aimed at performing a certain task is linked to [START_REF] Shalev-Shwartz | Using more data to speed-up training time[END_REF] and is a thrilling problematic. It has the potential to answer some profound questions about the optimal way to perform learning depending on the limitations of a particular practical framework in terms of available memory, computational time, and number of available training samples.

Denoting c j = u, ν j , the hypothesis on u implies that the c j are distinct. We then have for any x ∈ R, by replacing ω by xu in (A.3), N j=1 λ j e -ic j x = 0.

(A.4)

The functions x → e -ic j x defined on R are eigenfunctions of the linear derivation operator on C ∞ (R) associated to distinct eigenvalues -ic 1 , . . . , -ic N . Therefore, they form a linearly independent family and the λ j are all zero. Finally, the family P n is linearly independent.

Lemma 5 will principally be used in the following way: if a linear combination f = k s=1 λ s p µ s of distinct functions of P n is the zero function, then all the λ s are equal to 0. The next lemma relates the injectivity of an operator M on Σ k,n with the preimage of {0} by M in Σ 2k,n . Lemma 6. Let n > 0 and M be a linear operator on L 1 (R n ). Then the following properties are equivalent for any integer k > 0:

1. M is injective on Σ k,n .

Az -Ah -∆ δ (Mz) G ≤ C 1 d Y (z -h, Σ) + C 2 Mh F + δ. Az -Ah -∆ δ, (Mz) G ≤ C 1 d Y (z -h, Σ) + C 2 Mh F + δ.
Let's remark that (3.11) and (3.12) imply for all x, y ∈ G y G ≤ xy G + x G . Therefore, since Az -∆ δ (Mz) G ≤ δ (resp. Az -∆ δ, (Mz) G ≤ C 2 Mh F + δ), we have:

Ah G ≤ C 1 d Y (z -h, Σ) + C 2 Mh F + 2δ. (resp. Ah G ≤ C 1 d Y (z -h, Σ) + 2C 2 Mh F + 2δ.)

B.8 Proof of Theorem 8

Let's note M = M |V and N = N ∩ V . Let m be the dimension of the range of M, so that N is of dimension n -m. Let h 1 , . . . , h n-m be an orthonormal basis of N . We have:

n -m = n-m j=1 h j 2 2 ≤ 1 K n-
m ≥ n 1 - 1 K 1 - 1 D 2 * . (B.21)
We get the lower bound on D 2 * by isolating it in the inequality.

B.9 Proof of Theorem 9

Let h ∈ E and z ∈ Σ -Σ. We have the following inequalities: (B.22) where we have used the lower-RIP for the second inequality.

h G ≤ h -z G + z G ≤ h -z G + 1 α Mz F ,
A similar consideration on Mz yields:

Mz F ≤ M(z -h) F + Mh
h G ≤ h -z G + 1 α M(h -z) F + 1 α Mh F = h -z M + 1 α Mh F . (B.24)
Taking the infimum of the right hand-side quantity over all z ∈ Σ -Σ, one gets the desired Robust NSP: (B.25)

h G ≤ d M (h, Σ -Σ) + 1 α Mh F .

CHAPTER 2 .Figure 2 . 1 :

 221 Figure 2.1: Compressive learning outline. The learning data X = {x 1 , . . . , x L } is compressed into a smaller representation, which can either consist in reducing the dimensions of each individual entry x r or in computing a more global compressed representation of the data, called sketch. Parameters θ are then inferred from such a compressed representation by an algorithm adapted to it.

  Figure 2.2: Real and reconstructed centroids (respectively represented as circles and squares) of 4 Gaussians in dimension 2 from 10 3 points drawn from the mixture. To estimate the 12 real parameters of the mixture, the data was compressed to a complexvalued sketch of dimension 30, represented to the right as a 60-dimensional real signal.

Figure 2 .

 2 Figure 2.2 visually illustrates the behavior of the algorithm on a simple mixture of 4 Gaussians in dimension 2. L = 10 3 points were drawn from this mixture and used to compute a m = 30-dimensional sketch. As shown in the figure, the mixture parameters are precisely estimated without referring to the initial data. The symmetric KL divergence and Hellinger distance are respectively 0.026 and 0.003.Figure 2.3 illustrates the reconstruction quality of our algorithm in dimension 10 for different values of mixture components k and sketch sizes m in terms of Hellinger distance. For each sketch size m ranging from 200 to 2000 with step size 200, k was chosen to range from m/200 to m/10 with step m/200.For each choice of parameters, 10 experiments were performed and the depicted value is the Hellinger distance such that 80% of the experiments lead to a smaller Hellinger distance. We can essentially observe that the Hellinger distance gradually decreases as the number of mixture components rises. For the considered parameters range, choosing m = 10kn, i.e., choosing m so that it contains 10 times more values than the number of parameters to estimate, leads to a Hellinger distance smaller than 0.03 for 80% of the cases. Note that the number of measurements we took for the experiments is way below the deterministic number of 8k 3 n provided in

Figure 2 .

 2 Figure 2.2 visually illustrates the behavior of the algorithm on a simple mixture of 4 Gaussians in dimension 2. L = 10 3 points were drawn from this mixture and used to compute a m = 30-dimensional sketch. As shown in the figure, the mixture parameters are precisely estimated without referring to the initial data. The symmetric KL divergence and Hellinger distance are respectively 0.026 and 0.003.Figure 2.3 illustrates the reconstruction quality of our algorithm in dimension 10 for different values of mixture components k and sketch sizes m in terms of Hellinger distance. For each sketch size m ranging from 200 to 2000 with step size 200, k was chosen to range from m/200 to m/10 with step m/200.For each choice of parameters, 10 experiments were performed and the depicted value is the Hellinger distance such that 80% of the experiments lead to a smaller Hellinger distance. We can essentially observe that the Hellinger distance gradually decreases as the number of mixture components rises. For the considered parameters range, choosing m = 10kn, i.e., choosing m so that it contains 10 times more values than the number of parameters to estimate, leads to a Hellinger distance smaller than 0.03 for 80% of the cases. Note that the number of measurements we took for the experiments is way below the deterministic number of 8k 3 n provided in

Figure 2 . 3 :

 23 Figure 2.3: Quality of reconstruction in dimension n = 10, with N = 10 4 points, measured as a Hellinger distance. Each square corresponds to 10 experiments, and the depicted values are the values of the Hellinger distance under which 80% of performed experiments are placed.

Figure 3 . 1 :

 31 Figure 3.1: Illustration of several CS models. From left to right: k-sparse vectors, union of subspaces, smooth manifold and point cloud.

Figure 3 . 5 :

 35 Figure 3.5: Drawback of uniform instance optimality in a simple case: the model Σ (Left) is the set of black points including those on D and the point p and the operator M is the 1-dimensional orthogonal projection on the horizontal axis. If we choose ∆ as the pseudo-inverse of M, the depicted IO ratio (Right) is low on most of the space, but the uniform constant is infinite.

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: Impact of the size M of the learning set on the ranking performance for the exact search procedure, measured using recall@R curves. For the two datasets considered (Imagenet and SIFT1M), the vector output by the explicit embedding is about one half of the original descriptor dimensionality (E=512 and E=64, respectively).

  Figure 4.5: Illustration of the PQ procedure on one vector. The vector is considered as a product of D subvectors of equal length, and each subvector is quantized separately using a k-means algorithm on the database X .

Figure 4 .

 4 Figure 4.6: Comparison with the state-of-the-art: search quality for the χ 2 kernel with coded representations for 1M SIFT. Left: recall@R curves for a fixed size of B=256 bits. Right: recall@100 as a function of the number of bits (32-256). Parameters: KPCA+LSH → E = B and M =1,024; KPCA+PQ → E=128 and M =1,024.

Figure 4 . 7 :

 47 Figure 4.7: Comparison with the state-of-the-art: search quality for the χ 2 kernel with coded representations for 1.2M BOW. Left: recall@R curves for a fixed size of B=256 bits. Right: recall@1000 as a function of the number of bits (32-256). Parameters: KPCA+LSH → E = B and M = 1, 024; KPCA+PQ → E = 128 and M = 1, 024.

Figure 4 . 8 :

 48 Figure 4.8: Performance as a function of the database size (1M to 200M vectors). Parameters for KPCA+PQ: D = 16, E = 64 and M = 1, 024.

Figure 4 .

 4 Figure 4.8 reports the Recall@R measured on subsets of 1M-200M SIFT descriptors taken from the BIGANN dataset. These results show that the proposed KPCA+PQ scheme can scale to very large datasets without a drastic loss in performance.

2 .

 2 ker(M) ∩ Σ 2k,n = {0}. Proof. 1. ⇒ 2. Suppose M is injective on Σ k,n and let f ∈ ker(M)∩Σ 2k,n , that is Mf = 0 and f = 2k s=1 λ s p µ s , (A.5)with λ s ∈ R and p µ s ∈ R n for all s. Since Mf = 0, we can apply M to (A.5) to obtainM k s=1 λ s p µ s = M 2k s=k+1 -λ s p µ s . (A.6)The left and right-hand side of (A.6) are images by M of functions of Σ k,n . Since M is injective on Σ k,n , these functions are equal, which yields (M)∩ Σ 2k,n = {0}. 2. ⇒ 1. Let's suppose ker(M) ∩ Σ 2k,n = {0}. Let µ 1 , . . . , µ k , ν 1 , . . . , ν k be vectors of R n and α 1 , . . . , α k , β 1 , . . . , β k be real numbers such that M k s=1 α s p µ s = M k s=1 β s p νs . (A.8)By subtracting the right hand side in (A.8), we getM k s=1 α s p µ s -k s=1 β s p νs = 0. (A.9)Note that this set may not contain a unique element and thus this definition relies on the axiom of choice. For x ∈ E, we have again x -∆ ν (Mx) ∈ N , so that by(3.6):Ax -A∆ ν (Mx) G ≤ Dd E (x -∆ ν (Mx), Σ -Σ) ≤ Dd E (x, Σ) + Dd E (∆ ν (Mx), Σ) ≤ 2Dd E (x, Σ) + Dν, (B.10)where we have used (3.12) again for the second inequality. The last inequality comes from (B.9), which yields d E (∆ ν (Mx), Σ) ≤ d E (x, Σ) + ν. Therefore, by posing ∆ δ = A∆ δ/D , we get (3.13).B.4 Proof of Proposition 1Let x ∈ E and ν > 0. If 0 = d E (x + N , Σ) = d E (x, Σ + N ), then since Σ + N is a closed set,x ∈ Σ + N , and therefore (x + N ) ∩ Σ = ∅. In this case, we define ∆ ν (Mx) as any element of (x + N ) ∩ Σ.Ifd E (x + N , Σ) > 0, then we define ∆ ν (Mx) ∈ {u ∈ (x + N )|d E (u, Σ) ≤ (1 + ν)d E (x + N , Σ)}.This provides a consistent definition of ∆ ν . Let's remark that for all x ∈ E,d E (∆ ν (Mx), Σ) ≤ (1 + ν)d E (x + N , Σ) ≤ (1 + ν)d E (x, Σ)For x ∈ E, x -∆ ν (Mx) ∈ N , so that(3.6) gives:Ax -A∆ ν (Mx) G ≤ Dd E (x -∆ ν (Mx), Σ -Σ) ≤ Dd E (x, Σ) + Dd E (∆ ν (Mx), Σ) ≤ (2 + ν)Dd E (x,Σ). (B.11) Defining ∆ δ = A∆ ν , we get the desired result. B.5 Proof of Theorem 5 and Theorem 7 Let's first remark that applying (3.19) (resp. (3.21)) with x = z ∈ Σ and e = 0 yields Az -∆ δ (Mz) G ≤ δ and Az -∆ δ, (Mz) G ≤ C 2 + δ for any z ∈ Σ, ≥ 0, where we have used the fact that 0 F = 0. Let h ∈ E and z ∈ Σ. We apply (3.19) (resp. (3.21)) with x = zh, e = Mh, and = Mh F , which yields:

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . State of the art: density mixture estimation and compressed learning . . . . 2.1.1 Density mixture estimation as a linear inverse problem . . . . . . . . 2.1.2 Learning with randomized dimensionality reduction and sketching . 2.2 Layout of the chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Compressive estimation framework . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 The compressive operator. . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Proposed instantiation: isotropic Gaussians. . . . . . . . . . . . . . . 2.3.3 Injectivity of the compressive operator. . . . . . . . . . . . . . . . . 2.3.4 Recovery problem formulation. . . . . . . . . . . . . . . . . . . . . . 2.4 Compressive reconstruction algorithm . . . . . . . . . . . . . . . . . . . . . 2.4.1 Reminder of Iterative Hard Thresholding . . . . . . . . . . . . . . . 2.4.2 Proposed continuous case algorithm . . . . . . . . . . . . . . . . . . 2.4.3 Memory usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.4 Computational complexity . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of the contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Short-term perspectives . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Mid-term perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 Long-term perspectives . . . . . . . . . . . . . . . . . . . . . . . . . Preliminary lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	xii A.1 xiv 5 Conclusion 5.1 Bibliography List of Figures A Proofs of Chapter 2 theorems Notations Mathematical notations	CONTENTS CONTENTS
	x	Vector (bold lowercase letter)
	A	Linear operator (bold uppercase letter)
	M L I	Vector or operator family (curved uppercase letter) Identity matrix
	∇f •, • E[•]	Gradient of a functional f Integer interval Expected value of a random variable
	P(•)	Probability of an event
	2 Compressive density mixture estimation
	Manuscript conventions n Signal or data dimension m Compressed dimension / Number of measurements k Sparsity level L Number of elements in a database M Measurement operator Σ k Set of k-sparse vectors Σ General model of signals E Signal space F Measurement space G Feature space 2.1 xi xv

Introduction xvii 1 Motivation 1.1 Inverse problems and Compressed sensing . . . . . . . . . . . . . . . . . . . 1.1.1 Sparse signal representations . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Inverse problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Compressed sensing paradigms . . . . . . . . . . . . . . . . . . . . . 1.2 "CS-like" techniques in statistical analysis . . . . . . . . . . . . . . . . . . . 1.2.1 Locality-Sensitive Hashing . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Compressed methods for matrix decomposition . . . . . . . . . . . . 1.2.3 Other randomized statistical procedures . . . . . . . . . . . . . . . . 1.3 2.5.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Choice of the frequencies . . . . . . . . . . . . . . . . . . . . . . . . 2.5.3 Heuristic for random initialization . . . . . . . . . . . . . . . . . . . 2.5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Proofs of Chapter 3 theorems B.1 Well-posedness of the finite UoS decoder . . . . . . . . . . . . . . . . . . . . B.2 Proof of Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Proof of Theorem 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Proof of Proposition 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Proof of Theorem 5 and Theorem 7 . . . . . . . . . . . . . . . . . . . . . . . B.6 Proof of Theorem 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.7 Proof of Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.8 Proof of Theorem 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.9 Proof of Theorem 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p

(R n ) Lebesgue p-integrable complex functions on R n • p p -norm in R n or L p (R n ) •, •

Scalar product in a pre-Hilbert space F(f ).ω Fourier transform of f taken at frequency ω M m,n (K) Matrices of size m × n over the field K

Table 3 .

 3 1: Summary of the hypotheses on the pseudo-norms • E , • F and • G . A cross means the property is required, a horizontal bar means it is not. function needs not be a norm in order to state the following results. It still must satisfy the following properties:

Table 4 .

 4 1: Datasets used for the evaluation of nearest neighbor search

  Figure 4.1: Exact nearest neighbor search in a database of 10M 128-dimensional SIFT vectors (extracted from the BIGANN database), embedded in dimension 128 with KPCA on χ 2 kernel. Left: Procedure for a given query. The reranked vectors are the neighbors which have an upper bound value higher than the lower bound of the approximate nearest neighbor. Right: one minus the cumulative density function of the number of kernel values computed to find the exact nearest neighbor of a query (for different values of M ). The M kernel calculations required to compute the embedding of the query are taken into account.

	Kernel values and bounds	0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9	True K Approximate K Upper bound for database vectors (Equ. (4.23)) Lower bound for ANN (Equ. (4.22))	Complementary cumulative density function	0.1 0.5 1			M=512 M=1,024 M=2,048 M=4,096 M=8,192
		0	2000	4000	6000	8000	10000	0	0	10000	20000	30000
				Approximate neighbor rank					Number of kernel values computed

  Proposed approximate search scheme. As opposed as the existing implicit hashing methods, we first compute explicit vectors belonging to an Euclidean space, then apply an encoding method designed for Euclidean distances.

	Initial	Implicit hashing	Implicit signatures	Proposed scheme
	vectors		Explicit		Explicit
		Explicit embedding	vectors	Euclidean encoding	signatures
	Figure 4.4:

  Using(3.29), we get that, for all h ∈ N and unit-norm vectorz ∈ Σ -Σ, h, z 2 ≤ 1 -1 D 2 * h 2 2 .If we denote by p N the orthogonal projection on N and apply this inequality withh = p N (z i ) = n-m j=1 h j , z i h j and z = z i , we get that p N (z i ) 4 2 ≤ 1 -1 D 2 * p N (z i ) 2 2 , which can be simplified to p N (z i ) 2 2 = n-m j=1 h j , z i 2 ≤ 1 -1 D 2 * even if p N (z i ) 2 = 0.Using this relation in (B.19), we get:

				m	n
					h j , z i	2 .	(B.19)
			j=1	i=1
	n -m ≤	n K	1 -	1 * D 2	,		(B.20)
	so that:

, initialement énoncé pour un modèle de vecteurs k-parcimonieux puis généralisé à des unions de sous-espaces, à un cadre tout à fait quelconque de modèle. Nous montrons que ce résultat s'applique à de nombreux modèles classiques utilisés en problèmes inverses.

We now summarize more precisely our main contributions which fall under two categories: on the one hand the characterization of the existence of an instance optimal decoder given a linear measurement operator M and a model Σ, on the other hand the limits of dimensionality reduction.

-sense, except for a very limited neighborhood of the point cloud. This is further discussed in Section

3.7.

b In fact, unlike the NSP (3.6),(3.20) is not purely a property of the null space N even though it implies the NSP. The name Robust NSP is thus somewhat improper, but has become a standard for this type of property.

Remerciements

Abbreviations

Dataset SIFT1M Imagenet Pre-processing None RR RP None RR RP Recall@1 0.17 0.15 0.19 0.04 0.04 0.05 Recall@10 0.42 0.39 0.51 0.12 0.12 0.14 Recall@100 0.75 0.74 0.85 0.30 0.33 0.37 Recall@1000 0.96 0.97 0.99 0.60 0.67 0.69 Table 4.2: Choice of the pre-processing: performance of KPCA+PQ (Recall@R) with no pre-processing (Jégou et al., 2011), Random rotation (RR) [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF] and Random permutation (RP). Parameters: M = 1, 024, E = 64 and D = 8 for the two datasets considered. Averages on 10 experiments.

Comparison with the state of the art.

We now compare the proposed coding scheme based on KPCA+PQ with the state-of-theart techniques described earlier: KLSH [START_REF] Kulis | Kernelized locality-sensitive hashing for scalable image search[END_REF] and RMMH [START_REF] Joly | Random maximum margin hashing[END_REF]. We note that the proposed approach differs in two significant ways from KLSH and RMMH:

(i) As previously mentioned, KLSH and RMMH rely on an implicit embedding while we rely on an explicit KPCA embedding.

(ii) The signatures produced in KLSH and RMMH are binary while PQ encodes an embedded vector Ψ(x i ) as a succession of bytes.

To better understand the respective impact of these two differences, we also experimented with a coding scheme based on KPCA+LSH which performs an explicit embedding of the data (as in the first step of the proposed scheme) but performs a standard binary coding (as in the second step). Comparing this scheme with KLSH and RMMH allows one to measure exactly the gain of adding the explicit embedding step since the produced binary codes are very similar to those produce by KLSH and RMMH.

We performed experiments for different signature sizes (noted B in bits). The number of learning vectors for the KPCA was chosen as M = 1, 024, so that the number of kernel computations when processing a query is limited. A similar number of learning vectors was chosen as an input to KLSH and RMMH algorithms. For the embedding dimension E, we set it differently depending on the scheme. For KPCA+PQ, we set E = 128. For KPCA+LSH, we set E = B which led to optimal or near optimal results in preliminary experiments.

We report results for SIFT1M on Figure 4.6 and for Imagenet on Figure 4.7. We repeated all experiments 10 times, with different random draws for any step that relies on randomness. We show the average and standard deviation.

We can first observe that, as reported in [START_REF] Joly | Random maximum margin hashing[END_REF], RMMH outperforms KLSH on the BOW Imagenet features. However, for SIFT features, KLSH performs significantly better than RMMH. For instance, for B = 128 bits, the improvement in recall@100 is on the order of 10% absolute. Second, we can see that the simple KPCA+LSH scheme outperforms KLSH and RMMH at all ranks and code sizes. This difference can be significant (on the order of 10% absolute on SIFT1M for 128 bits). This shows the practical superiority of explicit embedding over implicit embedding for approximate search on these datasets. Finally, we can observe that KPCA+PQ performs significantly better than KPCA+LSH, which confirms the superiority of a non-binary encoding scheme over a binary coding scheme (Jégou et al., 2011).

Appendix A

Proofs of Chapter 2 theorems

This appendix is aimed at proving Theorems 1 and 2. Let's recall that the standard deviation of the Gaussians, denoted σ, is a fixed positive real number in all the subsequent sections of this appendix. For ω ∈ R n , M ω will denote the Fourier sampling operator on L 1 (R n ) at frequency ω. Recall that L 1 (R n ) denotes the complex vector space of Lebesgueintegrable functions mapping R n into C.

A.1 Preliminary lemmas

In this section, we consider n as any positive integer and prove two elementary lemmas which will be used in the proofs of Theorems 1 and 2.

The first lemma states a simple result:

Lemma 5. For any integer n > 0, the family

Proof. Let N > 0, (λ j ) N j=1 be complex numbers and (ν j ) N j=1 be real n-dimensional distinct vectors such that f = N j=1 λ j p ν j = 0. Let's prove that we necessarily have λ j = 0 for all j, which will prove the lemma.

Denoting by F the Fourier transform, we get that for any ω ∈ R n :

Since F(p 0 ) takes no zero value, this is equivalent to

Denoting f j the linear form ω → ω, ν j on R n , there exists u ∈ R n such that u / ∈ ker(f r -f s ) for all distinct r, s ∈ 1, N . A straightforward argument to justify this claim is to invoke Baire's theorem, since R n is complete and a hyperplane is a closed set of empty interior.

The argument is a function of ker(M) ∩ Σ 2k,n , so it is zero by hypothesis. Therefore, we necessarily have

and M is injective on Σ k,n .

In the proofs of Theorems 1 and 2, we will actually prove that ker(M) ∩ Σ 2k,n = {0}. Lemma 6 states that this implies the injectivity result we want to obtain. The converse will also be used to prove Theorem 2 from Theorem 1.

A.2 Proof of Theorem 1

In this section, we fix n = 1. The numbers ω 1 , . . . , ω 2k defined in the statement of the theorem are also fixed, and we note ω = (ω 1 , . . . , ω 2k ) ∈ R 2k . We define M as the linear Fourier sampling operator on L 1 (R) associated to the 4k 2 frequencies defined as the 2k first nonzero multiples of the ω j .

Let's now define a family of finite-dimensional linear operators which will be injective iff M is injective, as stated by the following lemma. Given real numbers ω and µ, define the function

, this new notation will be used to simplify the expressions. Given ω ∈ R and an -dimensional vector µ = (µ 1 , . . . , µ ), we denote by N ω,µ the following matrix:

We note N µ the matrix:

The following lemma links the injectivity of M and the injectivity of the N µ .

Lemma 7. The operator M is injective on Σ k,n if and only if the operator N µ is injective for any vector µ ∈ R 2k composed of distinct entries.

Proof. Let's first remark that for all coefficients α

The matrix N µ thus allows us to interpret M as a finite-dimensional linear map for linear combinations of Gaussians in the support µ. Let's now prove the equivalence.

Necessary condition. Suppose M is injective on Σ k,n . Let µ ∈ R 2k be a vector composed of distinct entries and γ ∈ ker(N µ ). This implies

Since the µ s are distinct, Lemma 5 implies that all γ s are 0, so that γ = 0. Finally, ker(N µ ) = {0} and N µ is injective.

With a proper reordering of the terms of the two sums, we can put this last inequality in the form:

where the first sum is the gathering of the common functions between the two sums in (A.16), the second and third are the distinct functions, and the fourth is composed of any other functions p ξs chosen so that the linear combination comprises 2k terms composed of distinct functions. The injectivity of N µ for µ = (θ 1 , . . . , θ k , ν +1 , ν k , ξ 1 , . . . , ξ ) gives that α s = β s for all s ∈ 1, and α s = β s = 0 for all s ∈ + 1, k , so that We are now ready to prove Theorem 1. Lemma 7 casts the initial claim as the injectivity of N µ for all µ ∈ R 2k with distinct entries. For such a vector µ, we have

Let's therefore study ker(N ω,µ ) for ω ∈ R. For s ∈ {1, . . . , 2k}, in order to simplify notations, let's note β s = exp -σ 2 2 s 2 ω 2 and x s = exp(-iωµ s ). N ω,µ can be written as:

The determinant of this matrix is proportional (within a nonzero factor) to a Vandermonde determinant. Therefore, this matrix has a nonzero kernel if and only if ∃s = t such that x s = x t , which is equivalent to

The structure of the matrix gives the form of ker(N ω,µ ) which is determined by the equivalence classes of the µ s in R\ where

Let's remark that U (I j ) = {0} if and only if I j is a singleton. Let's now show that the choice taken for the frequencies is sufficient to yield the injectivity of N µ .

For t = 1, . . . , 2k, we note I t 1 , . . . , I t Pt the partition induced on 1, 2k by the equivalence relation of equality modulo 2π ωt for the µ s , s ∈ 1, 2k . To prove the theorem, it is sufficient to prove that:

Indeed, this will prove that ∀s ∈ 1, 2k , ∃t ∈ 1, 2k , ker(N ωt,µ ) ⊥ e s , where e s denotes the s th vector of the canonical basis of R 2k . We will thus have: Let's therefore prove (A.24). The proof relies on the following fact: if x = y are in I s for some and s, then for all t = s, x and y cannot belong to a same I t j . Indeed, by contraposition, if x and y belong to I s and I t j for s = t, we have µ x -µ y = 0 2π ωs = 0 2π ωt . By hypothesis, ωs ωt / ∈ Q, so that 2π ωs Z ∩ 2π ωt Z = {0}. Therefore, µ x = µ y and x = y since µ is composed of distinct entries.

Let s ∈ 1, 2k . The above result tells us that for any t = s, t and s are in the same subset in at most one of these partitions. Since there are 2k partitions and only 2k -1 elements t distinct from s, there is necessarily at least one partition in which {s} is a subset. This yields (A.24), and therefore the result.

A.3 Proof of Theorem 2

We fix ω 1 , . . . , ω 2k as real frequencies satisfying the hypotheses of Theorem 1.

Given a family U = {u 1 , . . . , u m } ⊂ R n , we note M(U) the functional operator which samples the Fourier transform of a function at every frequency xω y u j for (x, y, j) ∈ 1, 2k

Proof. The assumption on U means that for any collection X of 2k vectors and any vector x ∈ X , there exists a direction u in U such that the orthogonal projection of x on u is different from the orthogonal projection of any other vector of X on u. Since U is fixed in the current proof, let's note M = M(U ). Without loss of generality, let's suppose the vectors of U are normalized in the 2 sense.

Using Lemma 6, it is sufficient to prove that for all collection of vectors µ 1 , . . . , µ 2k ∈ R n and all real numbers λ 1 , . . . , λ 2k ,

Without loss of generality, we can suppose the µ s are distinct (if some of them are equal, we can gather the common values and add any other p µ with a coefficient 0 until we get 2k terms). We note p = 2k s=1 λ s p µ s . We will now prove that every λ t is 0, which will prove the result. Let t ∈ 1, 2k . From the hypothesis, there exists u ∈ U such that (A.27) holds. In the expression of Mp as a vector of dimension 4k 2 m, let's consider the 4k 2 entries M xω j p, with (x, j) ∈ 1, 2k 2 . Using (2.18), we can get their expressions in terms of p:

x 2 ω 2 j exp(-ixω j u, µ s ).

(A.29) Let's now remark that the above quantity is the Fourier transform at frequency xω j of the mixture of 1-dimensional Gaussians q = 2k s=1 λ s p u,µ s .

(A.30)

Let N be the operator defined on L 1 (R) which samples the Fourier transform of a function at frequencies xω j for (x, j) ∈ 1, 2k 2 . From the above reasoning, we have Nq = 0, so q ∈ ker(N).

Theorem 1 ensures that N is injective on Σ k,1 and therefore Lemma 6 ensures that ker(N) ∩ Σ 2k,1 = {0}. Thus, q = 0. But the hypothesis on u ensures that u, µ t = u, µ s for all s = t. By regrouping the same functions p u,µ s in the expression (A.30), we get that q = λ t p u,µ t + s ξ s p u,µ s , (A.31) where the second sum contains at most 2k -1 terms and the µ s and µ t are distinct. Lemma 5 ensures that, since q = 0, we have λ t = 0. This is valid for all t ∈ 1, 2k , which proves the result.

Let's now prove a final lemma, which states the existence of certain remarkable families of vectors in R n which will be used in the rest of the proof.

Lemma 9. For any integer q ≥ n, there exists a family B q of vectors of R n such that any subfamily of n vectors of B q is a basis of R n .

Proof. We prove this by induction on q. For q = n, the property is trivial: it is sufficient to define B q as a basis of R n .

Let's suppose the property is satisfied for q ≥ n. Consider the hyperplanes of R n spanned by all the collections of n -1 vectors of B q : there are finitely many such hyperplanes. The union of these hyperplanes cannot recover R n : a straightforward argument to justify this claim is to invoke Baire's theorem, since R n is complete and a hyperplane is a closed set of empty interior. Therefore, there exists a vector x ∈ R n which is not in this union.

Let's define B q+1 = B q ∪ {x}, and let U be a subfamily of n vectors of B q+1 . If U ⊂ B q , then U is a basis of R n by hypothesis. If not, then U is composed of x and n -1 vectors of B q . These n -1 vectors are linearly independent: if not, adding one other vector of B q to them would produce a subfamily of n vectors of B q which is not a basis. Moreover, by hypothesis on x, it is not included in the hyperplane spanned by these n -1 vectors, so that U is a linearly independent family of R n composed of n vectors, that is a basis of R n . This proves the result.

Let's now finish the proof of the theorem. Let m be an integer > (2k -1)(n -1) and pose U = B m , where B m is the corresponding set in Lemma 9. We will prove that U satisfies the hypothesis of Lemma 8, which will prove the theorem.

Let x 1 , . . . , x 2k be a family of distinct vectors of R n , and t ∈ 1, 2k . For any s = t, there are at most n -1 vectors u of U such that u, x t = u, x s : indeed, if there are n such vectors, they form a basis and we would therefore have x t = x s , which contradicts the fact that the vectors are distinct. Therefore, the set [START_REF]Core idea of LSH[END_REF]. Since m > (2k -1)(n -1), there exists a vector u ∈ U such that for all s = t, u, x t = u, x s . This is valid for all t, so that U satisfies the hypothesis of Lemma 8, which proves the result, since we can pose m = (2k -1)(n -1) + 1 ≤ 2kn. (v). For all u such that u > M, we have Φ(u) ≥ x ≥ Φ(v), so that v is a global minimizer of Φ. We therefore have ∀(u, w) ∈ V × W, xv ≤ x -(u + w) . The vector v can be written f + g with f ∈ V and g ∈ W , so that the vector y = xf , which belongs to

Therefore, ∆(Mx) can be defined as

is properly defined. In particular, this applies to the decoder (3.4).

B.2 Proof of Theorem 3

Let δ > 0 and ∆ δ and C be such that (3.13)

Applying (3.13) to x = h 2 ∈ Σ and using the fact that 0 E = 0, we get:

Let's now find an upper bound for Ah G :

where we have used (3.11) and (3.12) for the last inequality. Combining (B.1) and (B.2), we get that:

Applying (3.13) to x = h 1 + h 3 , we get:

B.3 Proof of Theorem 4

Let's first assume that (3.14) holds and define the following decoder on F :

Note that the decoder is well defined, since

where we have used (3.12) for the second inequality. The last inequality comes from (B.7), which yields d E (∆ (Mx), Σ) ≤ d E (x, Σ). Therefore, by posing ∆ = A∆ , we get (3.5).

Let's return to the general case, and consider ν > 0. We define the following decoder on This last inequality is valid for all z ∈ Σ, therefore (3.19) implies: (B.12) where we have used (3.11) for the last inequality. Similarly, (3.21) implies

We conclude by using the fact that (B.12) and (B.13) hold for all δ > 0.

B.6 Proof of Theorem 6

Let's suppose (3.20) and define for δ > 0 the decoder ∆ δ : F → E such that ∀y ∈ F : where we have used (3.11) and (3.12) for the second inequality and the last inequality is a consequence of (B.14).

Posing ∆ δ = A∆ δ proves (3.19) with C 1 = 2D 1 and C 2 = 2D 2 .

B.7 Proof of Lemma 1

The two equivalences are very similar to prove, so that we will only prove the first. (3.23) ⇒ (3.22) is obvious. Let's now suppose (3.22), so that:

We also have: (B.17) so that: (B.18) This last inequality yields (3.23).

B.8. PROOF OF THEOREM 8