Nicolas Courty
email: nicolas.courty@univ-ubs.fr

Xing Gong

Jimmy Vandel

Thomas Burger

This paper presents a new non-negative matrix factorization technique which (1) allows the decomposition of the original data on multiple latent factors accounting for the geometrical structure of the manifold embedding the data;

(2) provides an optimal representation with a controllable level of sparsity; (3) has an overall linear complexity allowing handling in tractable time large and high dimensional datasets. It operates by coding the data with respect to local neighbors with non-linear weights. This locality is obtained as a consequence of the simultaneous sparsity and convexity constraints. Our method is demonstrated over several experiments, including a feature extraction and classification task, where it achieves better performances than the state-of-the-art factorization methods, with a shorter computational time.

Introduction

Context: Non-negative matrix factorization (or NMF for short) has long been studied and used as a powerful data analysis tool providing a basis for numerous processing, such as dimensionality reduction, clustering, denoising, unmixing, etc. This article is concerned with a variant of NMF, where one tries to decompose a given matrix X = [x 1 , . . . , x n ] (x i ∈R p ) formed by a set of n samples described with p variables (∈ R p ) as a product of two matrices F ∈ R p×ℓ and G ∈ R n×ℓ such that ||X -FG ⊤ || 2 is minimized, under the constraints that the elements of G are positive or nil, and that its rows sum to one. In other words, G serves as a convex

SAGA: Sparse And Geometry-Aware non-negative matrix factorization through non-linear local embedding surrogate description of X in the reduced embedding formed by ℓ prototypes (the columns of F).

This convexity property on G is desirable in a lot of applications for its relations with the physics of the underlying observations: in this acceptation, the datapoints are usually observations of a mixing process where the components of the mixture are not known. Hyperspectral images are a good example, as each pixel describes a spectrum that can be defined as a combination of pure materials spectra (trees, concrete, water, etc.); the combination occurring because of the captor spatial resolution and different scattering effects [START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF][START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. Other examples of applications are found in archetypal analysis [START_REF] Mørup | Archetypal analysis for machine learning and data mining[END_REF], or biology [START_REF] Kersting | Pre-symptomatic prediction of plant drought stress using dirichlet-aggregation regression on hyperspectral images[END_REF].

Notations: Bold capital letters, such as X or M refer to matrices. The transposition operator is denoted . ⊤ (so that M ⊤ refers to the transpose of M). The set of m × n real matrices (respectively real matrices with nonnegative entries) is denoted R m×n (respectively R m×n + ). I n is the n × n identity matrix, and ||.|| refers to the Frobenius norm. For any matrix M, M •i (respectively M i• ) corresponds to the ith column (respectively ith row) of M. However, in the particular case of X, each column X •i corresponds to the datum x i , so that this latter more intuitive notation is preferred. ∆ n (M) is a (n -1)-simplicial polytope formed by n columns of M. The unit (n -1)-simplex ∆ n (I n ) is simply referred to as ∆ n .

Related work: Despite appealing properties, NMF presents a number of difficulties: First, while real data are often embedded on complex manifolds, the seminal convex NMF formulation prohibits a large number of latent factors that would explain the data in a geometry preserving way; Second, the level of sparsity in the final embedding (the matrix G) is rarely controllable, and far from data analyst expectations; finally, the computational time is prohibitive when it comes to handling very large size matrices. Even if to date, various state-of-the-art methods already tackle one of or two of these issues, no method addresses these three issues all together.

The convexity constraint we consider on G differs from that of what is classically referred to as Convex NMF [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF], i.e. NMF where the columns of F are convex combinations of columns of X. However, our problem is also not unheard of in the literature, as, beyond its easiness of interpretation, it provides a very appealing computational framework: if G encodes convex combinations, the ℓ columns of F are expected to characterize a hull of X in R p , as data points lying within are perfectly reconstructed. Based on this idea, it is possible to separate the NMF problem into two steps, namely the computation of F (in order to determine, or to best approximate the hull of the dataset) and of G (which corresponds to the projection of the dataset onto the region bound by the hull). From a computational point of view, this is really interesting, as it is now possible to decompose the NMF (which amounts to finding two matrices with close relationship) into two simple problems, each focused on a particular matrix: first one computes F regardless G, and second, one computes G by projection of the columns of X onto F. With such an approach, it is possible to expect a lower complexity than that of elder methods based on singular value decomposition (with a o(n 3 ) complexity), or than that of demanding procedures that iteratively alternates between the minimization of the two matrices, such as in [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF].

Initially, F was first related to the convex hull of X [START_REF] Cutler | Archetypal analysis[END_REF]: Notably, numerous works in the remote sensing and hyperspectral imaging community (see [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] as well as its references) have pointed out that F should be chosen so that ∆ ℓ (F) best encompasses X, while other investigations focused on how to find this simplicial convex hull [START_REF] Wang | Nonnegative least-correlated component analysis for separation of dependent sources by volume maximization[END_REF][START_REF] Thurau | Yes we can: simplex volume maximization for descriptive web-scale matrix factorization[END_REF][START_REF] Çivril | On selecting a maximum volume sub-matrix of a matrix and related problems[END_REF]. More recently, a series of works [START_REF] Arora | Computing a nonnegative matrix factorizationprovably[END_REF][START_REF] Kumar | Fast conical hull algorithms for near-separable non-negative matrix factorization[END_REF][START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF][START_REF] Recht | Factoring nonnegative matrices with linear programs[END_REF] focused on the conical hull of X: Here, the idea is to find F so that it spans a cone encompassing the dataset. The reason of the recent focus on the conical hull is that it is the geometric translation of an important assumption in NMF, the separability assumptions proposed by Donoho in [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts?[END_REF]. This separability assumption reads that (1) the residue X -FG ⊤ is nil, (2) columns of F are collinear to columns of X. This second condition regarding the collinearity of the columns of F and X is interesting beyond the separability assumption: whatever the type of hull defined by F, it is possible to assume that the components of the mixture belong to the data, and that any datum is a convex combination of a restricted number of particular selected datapoints (those forming F). This approach drastically reduces the complexity of finding F as this latter is simply defined by a column subset selection (CSS) of X.

Fig. 1 illustrates a dataset lying on a non-linear manifold resulting from some hidden factors that could be of interest, as well as its conical and convex hulls. It clearly appears that the geometries of these hulls are not adapted to that of the dataset. So far, there has been little interest in trying to characterize the boundary of the manifold dataset in spite of its non-convexity; in the sequel, we shall address this boundary with the shorter and imaged name manifold hull, that is illustrated on Fig. 1(d). Naturally, precisely characterizing such a manifold hull would require to increase ℓ, the number of datapoints involved in the CSS; which stresses to a larger extend the need for an adapted control of the sparsity level that is already sought for in numerous applications: Each point should be described as a convex combination of a restricted number λ of prototype points among the ℓ of the CSS.

Whatever the type of hull (convex or conical), separating F and G computations has drastically reduced the overall NMF complexity of the state-of-the-art methods: [START_REF] Thurau | Yes we can: simplex volume maximization for descriptive web-scale matrix factorization[END_REF][START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF] compute F in linear time and few works [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF][START_REF] Recht | Factoring nonnegative matrices with linear programs[END_REF] even reach overall linear-complex NMF. However, among them, none allows characterizing the manifold hull, and none allows controlling the solution sparsity. Even if to date, more computationally efficient methods than adding the classical L 1 penalty [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF][START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF] have been developed, such as [START_REF] Kim | Sparse non-negative matrix factorizations via alternating nonnegativity-constrained least squares for microarray data analysis[END_REF][START_REF] Gillis | Sparse and unique nonnegative matrix factorization through data preprocessing[END_REF], none reaches a linear complexity.

Even if the characterization of the manifold hull has never been addressed so far, it is well-known that the kernel trick is an efficient way to provide a manifold preserving description of a dataset. In fact, several pre-existing works applied the kernel trick to NMF [START_REF] Da | Non-negative matrix factorization on kernels[END_REF][START_REF] Buciu | Nonnegative matrix factorization in polynomial feature space[END_REF][START_REF] Cai | Graph regularized nonnegative matrix factorization for data representation[END_REF]. However, the convexity constraint on G (and consequently the notion of manifold hull) does not appear. Moreover, neither sparsity nor any linear complexity is achieved.

Proposal:

In this article, we present SAGA (Sparse And Geometry-Aware), a new NMF method avoiding the aforementioned limits. It operates in a Reproducing Kernel Hilbert Space (RKHS) [START_REF] Schölkopf | Learning with kernels: Support vector machines, regularization, optimization, and beyond[END_REF], where F is defined by a CSS. Then, according to an expected sparsity level, G, the best dataset projection onto the simplex formed by F is computed. The advantages of our method are:

First, kernelization is interesting to several extents: i) it makes the algorithm compliant with dataset where only relationships among objects are available; ii) Fig. 1 (a) A dataset X embedded in a U-shaped manifold; (b) the separability assumption assumes the dataset is encompassed in a cone spanned by the columns of F; a datum is described in a column of G through its non-negative coordinates in the cone; (c) Convex hull: the dataset is embedded in a simplex, and each datum is described with barycentric coordinates; (d) Manifold hull: one uses a great enough number of reference points to precisely characterize the manifold, while forcing the sparsity to achieve some local coding. Note that, in order to make these imaged representations clearer, we did not represent F as a CSS of X.

its regularization property improves robustness to noise; iii) it allows using more latent factors than the dimensions of the input data, providing an insightful tool to consider the geometric structure of the data manifold. As both the number of latent factors and the sparsity level are increasing, it appears that the locality of the support of G is increasing, turning the factorization problem in a non-linear local embedding of the original data. As shown in some recent works [START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF][START_REF] Guillemot | Neighbor Embedding with Non-negative Matrix Factorization for image prediction[END_REF], this kind of embedding is powerful to describe the non-linear structure of the data manifold, and as such serve as a very good feature extraction framework.

Second, the CSS is defined thanks to a manifold subsampling method [START_REF] Shroff | Manifold precis: An annealing technique for diverse sampling of manifolds[END_REF], which reaches a linear complexity with respect to n, the size of the dataset.

Third, the computation of G corresponds to a sparse RKHS simplex projection, which is a non-linear optimization problem. Based on recent advances on sparse projected gradient methods, the projection is solved with an algorithm of linear complexity (with respect to n), while naturally embedding sparsity mechanism control.

Forth and finally, since both computations of F and G are linear, the overall complexity of SAGA is linear. This makes this algorithm perfectly suitable for very big data.

To the best of our knowledge, no state-of-the-art method simultaneously compiles all these advantages.

Contributions: SAGA has been developed on the top of several state-of-theart algorithms, thanks to several technical extensions that are listed here. The computation of matrix F is largely inspired by the simplex volume maximization (SiVM) approach of [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]. However, the method itself has been extended in several manners:

1. It is slightly generalized in order to operate in the Hilbert space reproducing the Gaussian kernel, rather than in R p . 2. Most importantly, while SiVM proposes an approximate optimization (for several simplifications are made to reach a linear complexity), we propose another solution to the simplex volume maximization which performs an exact optimization.

In spite of these two extensions, our method remains of linear complexity, so that finally, the computation of F with SiVM is both a particular case and an approximation of the one produced by SAGA. Then, the computation of matrix G is not inspired by any other NMF technique, but is based on a recent sparse projected gradient strategy [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF]:

3. This latter is adapted to the projection over the CSS in the RKHS, and its linear complexity is kept. 4. We provide with theoretical bounds on the convergence of the projector (linked to the kernel bandwidth, the minimum pairwise distance between the CSS and the sparsity level).

Outline: We solve our NMF problem by separating the computations of the CSS matrix and of the projection matrix. This is why, Sections 2 and 3 focus on the computations of F and G respectively. More specifically, the structure of Sections 2 is the following: First a small introductive paragraph recalls the basics of the kernel trick, as well as why it is interesting to operate in a RKHS to fulfill our objectives. Then, Section 2.1 investigates the consequences of working in such RKHS; They lead us to a particular strategy, which is to characterize the manifold hull of the dataset in R p via a simplicial convex hull in the RKHS; and to define this latter with a manifold sampling strategy. At this point, computational constraints direct us toward SiVM-like procedures rather than toward more resource-demanding ones. Section 2.2 jointly presents SiVM, such as defined in the literature [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF], as well as the kernel generalization we propose, while Section 2.3 describes our modifications to reach exact maximal simplex volume. In a similar way, Section 3 is divided into two parts: the first one (Section 3.1) explains how projecting the image in the RKHS of any datum x i onto the image of the CSS in the RKHS, along with sparsity constraints; Section 3.2 provides with a proof that the projector defined in Section 3.1 converges. As it is established in the literature [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] that the Restricted Isometry Property (RIP, [START_REF] Candes | The restricted isometry property and its implications for compressed sensing[END_REF]) implies the convergence, we prove that our projector respect the RIP for some particular tuning of its parameters. At this stage, most of the required mathematics is established. Then, in Section 4, one summarizes the entire procedure made of the concatenation of the two matrices computation through an easy to implement algorithm, completed by some theoretical assessments of the linear complexity of the algorithm. Finally, Section 5 is devoted to experimental validations. In the first series of experiments, one focuses on toy examples, in order to illustrate the behavior of the manifold hull. Then, follow several experiments on simulated datasets, in order to compare the computational efficiency and the computational precision of SAGA with respect to the state-of-the-art. Finally, we consider real datasets through challenging image classification tasks.

2 Geometry aware CSS procedure Despite living in R p , X may span a nonlinear manifold of intrinsic dimensionality lower than p. A major problem is thus to extend to this nonlinear manifold the classical statistics that are used to work in a vector space. To do so, an essential element is to replace the Euclidean distances by geodesic distances. Depending on the manifold, the associated metric may be difficult to formally define. However, it is possible to characterize it through the time-scale of the well-studied heat diffusion process, the kernel formulation of which is well approximated by the Gaussian kernel [START_REF] Lafferty | Diffusion maps and coarse-graining: A unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF]27]: The geometry of the manifold in which X lies is captured by

K ij = k(x i , x j ) = exp -||x i -x j || 2 /2σ 2 with variance σ 2 .
Let us call φ(.) the implicit feature map from R p onto H, the RKHS associated to k(., .). We use the shorthand notation Φ = φ(X) = [φ(x 1 ), . . . , φ(x n )] (φ(x i )∈H) . Then, following the notations of [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF], the SAGA solution amounts to finding the indicator matrix W (ℓ) ∈ {0, 1} n×ℓ defining the CSS1 and the projection matrix G where φ(F) = φ(XW (ℓ) ) = ΦW (ℓ) , such that ||Φ -ΦW (ℓ) G ⊤ || 2 is minimized under convexity constraints.

CSS as a manifold sampling procedure

Let us start by a basic remark, Remark 1 In the input space, we must have ℓ ≤ n and ℓ ≤ p. In the Gaussian RKHS, one still has ℓ ≤ n, however, ℓ > p becomes possible, for each sample spans its own dimension.

which leads to the following property:

Property 1 (Non-Separability) Separability assumption does not hold for NMF in the Gaussian RKHS.

Proof. According to Remark 1, any datum not included in XW cannot be expressed as a linear combination of elements of φ(XW).

Thus, one should not consider conical hull in the RKHS. However, Remark 1 leads to:

Corollary 1 In the Gaussian RKHS, it is possible to use more than p points in the CSS. The latter forms a non-simplicial polytope in R p while their image in the Gaussian RKHS is a simplex.

In other words, the manifold hull of X can be characterized through the simplicial convex hull ∆ ℓ (φ(F)). Then, it follows that:

Corollary 2 No sample lies in ∆ ℓ (φ(F)) and all the samples will be projected on hyperfaces of ∆ ℓ (φ(F)), leading to approximate reconstructions. Yet, such an approximation comes with the appealing property of sparsity, discussed later in the article.

At this point, finding ∆ ℓ (φ(F)) amounts to finding W, which turns out to subsample the boundary of the data manifold. Most of the methods from the literature address it with the objective of maximizing the representation of the dataset, while here, we are interested in its boundary, which makes the sampling completely different: For instance, a kernel k-means sampling, although very efficient to subsample a given manifold [27], leads to Convex NMF of [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF], the aim of which is completely different of ours.

However, our problem is not completely unheard of: In [START_REF] Shroff | Manifold precis: An annealing technique for diverse sampling of manifolds[END_REF], the authors consider maximizing the diversity of the selected samples using a Karcher variance2 maximization criterion. Alternative formulation exists, where one seeks for the maximum volume parallelepiped in the data matrix [START_REF] Çivril | On selecting a maximum volume sub-matrix of a matrix and related problems[END_REF]. Interestingly enough, whatever the interpretation (diversity or volume criterion), the corresponding computation can reduce to recursive QR decompositions of the data matrix [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF][START_REF] Shroff | Manifold precis: An annealing technique for diverse sampling of manifolds[END_REF]. However, in the RKHS, since no explicit coordinates are available, those methods cannot be transposed. The recent proposal of [START_REF] Courty | A kernel view on manifold sub-sampling based on karcher variance optimization[END_REF] regarding a kernel rank revealing Cholesky technique is also of interest, unfortunately, it fails to scale up to big data, because it implies a computationally demanding decomposition at each selection step.

Finally, among all the methods available in the literature, if one discards those (1) which do not sample the boundary, (2) which cannot be conducted in a RKHS, (3) which do not have a linear complexity, we are aware of a single remaining method: the Simplex Volume Maximization (SiVM) [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF].

Original Simplex Volume Maximization

We begin with a review of the original formulation of Thurau and co-workers [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF] and its direct transposition to the kernel framework. SiVM tries to maximize Vol(∆ ℓ (F)), the volume of the simplex spanned by ∆ ℓ (F), which reads:

Vol(∆ ℓ (F)) = -1 ℓ • cmd(F) 2 ℓ-1 (ℓ -1)! , with cmd(F) = det       0 1 1 1 . . . 1 1 0 d 2 1,2 d 2 1,3 . . . d 2 1,ℓ 1 d 2 2,1 0 d 2 2,3 . . . d 2 2,ℓ 1 d 2 3,1 d 2 3,2 0 . . . d 2 3,ℓ . . . . . . . . . . . . . . . . . . 1 d 2 ℓ,1 d 2 ℓ,2 d 2 ℓ,3 . . . 0       . (1)
cmd(F) is the Cayley-Menger determinant of F (i.e. the determinant of the matrix accounting for the pairwise distances in F) and d 2 i,j is the square Euclidean distance between elements x i and x j . This obviously transposes to the RKHS: The volume of ∆ ℓ (ΦW) can be expressed by replacing the Cayley-Menger determinant by that of a matrix accounting for pairwise distances in H:

Vol(∆ ℓ (ΦW)) = -1 ℓ 2 ℓ-1 (ℓ -1)! det(A) (2)
with A similar to the matrix of Eq. ( 1) except that ∀i, j ≤ ℓ, d 2 i,j is replaced by:

A (i+1)(j+1) = ||φ(x i ) -φ(x j )|| 2 = k(x i , x i ) + k(x j , x j ) -2k(x i , x j )
The comparison of the volume of the simplices spanned by all the possible CSS is computationally prohibitive. Thus, an approximate search in linear time of the best simplex is proposed in [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]. This is possible thanks to two tricks: The first one is a strong result (Theorem 1 of [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]), which states that if one adds an element to the simplex, it can only make the reconstruction error smaller or equal. Thus, a simple greedy algorithm can be used to compute W: Starting with the best 1-simplex (i.e. W has 2 columns), one adds the best third sample (according to the maximization of Eq. ( 1)), then the fourth, and so on until ℓ samples are selected. Practically, this procedure can be transposed to the RKHS: At iteration p -1 one selects the element φ(x i ) of Φ such that i = arg max q Vol (∆ p (ΦW) ∪ φ(x q )) .

(3)

However, this procedure still requires the computation of several Cayley-Menger determinants (each of them being computationally intensive). At this point shows up the second trick: If one makes the simplifying assumption that the distances between the elements of the CSS defined in the previous iteration are constant and noted a, and if α j,q = d 2 j,q /2, then, Eq. 3 amounts to finding φ(x i ) such that

i = arg max q   p k=1 α kq •   a 2 + 2 p j=k+1 α jq   -(p -1) • p k=1 α 2 kq   , (4)
where d j,q refers to the distance between a point of the CSS x j and a point out of the CSS x q which is considered for adjunction to the CSS (thus, in the RKHS, one has α j,q = 1 -k(x j , x q )). Finally, the computation of Eq. 4 is sped up by considering that d 2 i,q ≈ d i,q and a 2 ≈ a, leading to

i = arg max q   p k=1 d kq •   a + p j=k+1 d jq   - p -1 2 • p k=1 d 2 kq   . (5)
Remark 2 Naturally, the magnitude to which the constant distance assumption is violated strongly depends on the dataset. As a consequence, the instantiation of Eq. 3 into Eq. 5 may lead to some approximations as already noted in [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF], where Gillis and Vavasis showed that SiVM can underperform on ill-conditionned datasets.

Despite this remark, the general principle advocated in Eq. 3 remains unalterably valid. This is why, we rely on it to propose an alternative to SiVM, which provides exact volume computation, with a similar linear complexity.

Exact simplex volume maximization

Let us consider a simplex ∆ p of dimensionality p-1 spanned by a subset XW (p) = {x 1 , . . . , x p } of p points of X. If we add a (p + 1)th point x i ∈ X \ XW (p) to the simplex, the new volume is given by:

Vol(∆ p+1 ) = Vol(∆ p ) × dist(x i , XW (p) ) p (6)
where dist(x i , XW (p) ) is the distance between x i and its projection onto the subspace spanned by XW (p) . According to [START_REF] Courty | PerTurbo: a new classification algorithm based on the spectrum perturbations of the laplace-beltrami operator[END_REF], in H, this distance reads:

dist(φ(x i ), ΦW (p) ) = 1 -k ⊤ x i • K -1 p • k x i (7)
where K -1 p is the inverse of the kernel matrix of the elements of the CSS, i.e.

K -1 p = (W (p) ⊤ Φ ⊤ ΦW (p)
) -1 , and where k x i is a vector of length p such that

k x i = [k(x j , x i )] x j ∈XW (p) .
Then, it is possible to use a greedy procedure similar to that of SiVM, where Eq. 3 translates into:

i = arg max q Vol(∆ p (ΦW (p) )) × dist(φ(x q ), ΦW (p) ) p = arg min q k ⊤ x q • K -1 p • k x q . ( 8
)
Remark 3 This procedure tends to add datum that most changes the geometry of the manifold spanned by the CSS. As such, it acts as the spectral sampling procedure proposed in [START_REF] Öztireli | Spectral sampling of manifolds[END_REF] for sampling a 3D mesh in a computer graphics context.

Even if K -1 p
is computed a single time at each iteration (it does not depend on x q ), a matrix inversion remains a resource demanding operation. Moreover, if ℓ elements are to be selected, then ℓ inversions of matrices of increasing size are expected. Fortunately, it is possible to bypass this inversion, by iteratively constructing K -1 p on the basis of the Schur complement [START_REF] Boyd | Convex optimization[END_REF]. Once i, the index of the best (p + 1)th point to add to the CSS is defined (Eq. 8), one computes

K -1 p+1 = K p k x i k ⊤ x i 1 -1 = K• 1-k ⊤ x i • K -1 p • k x i -1 0 ⊤ p 0 p K p -1 •K ⊤ (9)
with

K = -K -1 p • k x i I p 1 0 p ⊤
where I p is the identity matrix of size p and 0 p is a vector of p zeros. The computation works as long as k ⊤ x i • K -1 p • k x i differs from 1, which is always true in the Gaussian RKHS as long as the data points are separated.

Sparse projections onto the RKHS simplex

In this section, we focus on the computation of G. We give the formulation of our sparse RKHS simplex projector, and then we discuss its convergence. We notably show some analytical bounds required for the convergence of the method.

Projection on the RKHS simplex

We search for the projection of any point φ(x i ) onto the simplex ∆ ℓ (ΦW), i.e. the point of the simplex which minimizes the Euclidean distance to φ(x i ). It amounts to solving the n independent problems of computing the rows of G:

G i• = arg min G i• ||φ(x i ) -ΦWG ⊤ i• || 2 s. t. j G ij = 1, G ij ≥ 0, ∀j (10)
The

constraint j G ij = 1, G ij ≥ 0, ∀j is equivalent to have G i• in the unit standard simplex ∆ ℓ , (G encodes the barycentric coordinates of X in ∆ ℓ (ΦW)).
At this point, it is possible to force the sparsity of the projection to λ < ℓ, without any extra computational cost: We only need to replace the previous constraint by G i• ∈ ∆ λ . Thus, Eq. ( 10) reads:

G i• = arg min G i• ||φ(x i ) -ΦWG ⊤ i• || 2 s. t. G i• ∈ ∆ λ (11)
Instead of considering quadratic programming, such as in [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF] or [START_REF] Kumar | Fast conical hull algorithms for near-separable non-negative matrix factorization[END_REF], we follow some recent work on the projection on the unit standard simplex [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF], and we propose to use a simple projected gradient descent algorithm to solve Eq. ( 11), which amounts to iterating through different possible solutions of

G t+1 i• = P λ G t i• -ε t ∇(||φ(x i ) -ΦWG t i• ⊤ || 2 ) (12)
t being the iteration index, ε t a (possibly varying) step size, ∇(.) the gradient operator and P λ (.) the projector onto ∆ λ . This kind of projected gradient descent method has recently shown its computational efficiency and is also endowed with theoretical convergence guarantees [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF]. The gradient reads (we omit the iteration index t for clarity):

∇(||φ(x i ) -ΦWG ⊤ i• || 2 ) = ∇((φ(x i ) -ΦWG ⊤ i• ) ⊤ (φ(x i ) -ΦWG ⊤ i• )) = ∇(G i• K ℓ G ⊤ i• -2k x i G ⊤ i• + k(x i , x i )), = 2(G ⊤ i• K ℓ -k x i ). (13)
As for P λ (.), we rely on the Greedy Selector and Simplex Projector (GSSP) algorithm of [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF] which can be summarized as a two-step procedure: firstly the coordinates of the vector are sorted by magnitude, and then the λ greatest values are projected on the unit simplex ∆ λ (while the other vector entries are set to zero).

Remark 4 Since the GSSP procedure projects G i• on the subspace spanned by the ℓ columns of F, the coordinates of G i• embeds the projection over each of the ℓ selected elements in the feature space, so that the sparsity support is chosen in the closest elements in the feature space.

The Gaussian kernel is monotonically decreasing according to the neighboring distance; This implies that for each datum, the sparsity support is made of the closest CSS elements in the input space. It follows that:

Property 2 (Non-linear local coding)The sparse RKHS simplex projector describes any element x i of X with G i• , which interprets as its non-linear (because of the kernel non-linearity) barycentric coordinates according to λ prototype points. These prototype points are found in the closest elements of the CSS, thus providing a non-linear local coding for x i .

Convergence of the projector

Finally, we establish the convergence of that projector, which ensures [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF] that the final vector G i• is the best λ-sparse solution. To do so, we rely on [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF], which states that for Eq. 11 to be minimized via the projected gradient approach, ΦW has to satisfy λ-restricted isometry property (or λ-RIP for short), with δ λ ∈ [0, 1[. This latter reads:

Definition 1 The linear operator ΦW respects the λ-restricted isometry property [START_REF] Candes | The restricted isometry property and its implications for compressed sensing[END_REF] with constant δ λ if

(1 -δ λ )||x|| 2 2 ≤ ||ΦWx|| 2 2 ≤ (1 + δ λ )||x|| 2 2 ( 14
)
where ||.|| 2 refers to the L 2 norm, and for every λ-sparse vector x.

Equivalently, the constant δ λ can also be defined as:

δ λ := max L⊆{1,...,ℓ}, |L|=λ ||(ΦW) ⊤ L (ΦW) L -I λ || 2 (15)
where (ΦW) L denotes a subset matrix of ΦW, with λ columns corresponding to a subset L of cardinality λ picked up among the ℓ indices of the CSS. Thus, (ΦW) ⊤ L (ΦW) L is simply the related Gram matrix, noted K λ . δ λ is defined according to the subset providing a maximum among all the possible combinations of those columns.

As a matter of fact, such convergence holds for particular values of λ and σ, such as stated by the following proposition:

Proposition 1 If λ > 2, and if σ < d min √ 2 ln(λ-1)
, then, ΦW satisfies the λ-RIP

with constant δ λ ∈ [0, 1[. If λ ≤ 2,
there is no particular bound to σ.

Proof. Let us first note that since we are working in the Gaussian RKHS, the columns of ΦW have unit norms. The diagonal entries of K λ are 1, and ∀i, j ∈ L the off-diagonal element (i, j) is k(x i , x j ).

Let ν i be an eigenvalue of a matrix A. By the Gershgorin circle theorem, we know that:

|ν i -A(i, i)| < j≤λ,i =j |A(i, j)| (16)
Thus, for A = K λ -I λ , we obtain |ν i | < j≤λ,i =j k(x i , x j ). Let µ be the greatest dot product between the elements of (ΦW) L , i.e. µ := max k(x i , x j ), ∀i, j ∈ L, i = j. We can write:

j≤λ,i =j k(x i , x j ) ≤ (λ -1)µ (17)
and we have a bound for every eigenvalue ν i and every subset L of cardinality λ. Thus, ΦW follows the λ-RIP with δ λ = (λ -1)µ [START_REF] Bandeira | The road to deterministic matrices with the restricted isometry property[END_REF]. Let d 2 min be the minimum squared distance between two distinct elements of XW, i.e. d 2 min = min i,j,i =j ||x ix j || 2 . As µ = exp -d 2 min /2σ 2 , one has:

δ λ < 1 ⇔ (λ -1) exp -d 2 min 2σ 2 < 1 ⇔ σ < d min 2 ln(λ -1) , λ > 2 (18)
This allows deriving conditions on σ for different convergence or approximation guarantees, such as δ 2λ < 1/3 for convergence [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] (Theorem 2.1). This point is discussed in [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF] (Section 2). We note that similar developments using the Gershgorin circle theorem have been used for the deterministic creation of projection matrices in the domain of compressed sensing [3] (Section 2.1).

Complete SAGA procedure

We give in Alg. 1 the complete SAGA matrix factorization procedure 3 . We then discuss its computational complexity.

Algorithm

In addition to the data matrix X, the user needs to tune the following parameters: σ the Gaussian kernel bandwidth, λ the expected sparsity level, and ℓ the number of prototypes in the CSS (with λ ≤ ℓ). Additional parameters for the gradient descent can be tuned to optimize the computation load, yet it is not mandatory. At first, one computes the CSS, then the matrix of projections. Regarding the CSS, the procedure is initialized (Lines 1 to 4) as in [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]: First, one randomly selects a datum. Then, one finds the most distant datum to that first one. Finally, the most distant datum to this second datum is selected, and is considered as the first element of the CSS (see Lines 2 and 3 in Alg. 1). After initialization, the iterative increment of the simplex dimensionality is implemented in the loop from Line 5 to 8. Regarding the projection, each datum is processed separately thanks to the loop from Line 9 to 16. Within this loop, another loop deals with the gradient descent up to a stopping criterion (loop from Line 11 to 16).

Computational complexity

Property 3 (Computational complexity for W -noted as First step in Alg. 1 )

The selection of the CSS defining W based on the exact and incremental simplex volume maximization has complexity of o(nℓ 3 ), i.e. it has a linear complexity with respect to n.

Proof. The entire CSS construction is based on Eq. 8 and 9. At each step of the selection, the procedure amounts to a linear scanning procedure where the volume increment is computed for each element in the dataset (size n). For one element x i , with i ∈ [1, n], this requires to compute the associated bilinear form k ⊤ x q •K -1 i-1 •k x q k x q vector. with an assymptotical computational complexity of o(p 2 ). As p varies 1 to ℓ, one ends up with a cubical complexity term regarding ℓ.

We note that thanks to Eq. 9, almost all the values (except the one corresponding to the kernel evaluation with the last chosen element) have already been computed in the previous iteration and do not need to be computed again. This makes the overall CSS computation that efficient. However, for exceptionally large data matrices, it is possible to improve it with the randomized approximate search Algorithm 1: The SAGA Matrix Factorization algorithm input : X = [x 1 , . . . , xn] (x i ∈R p ) the data matrix σ the Gaussian kernel bandwidth λ the expected sparsity level ℓ the number of prototypes E = {ǫ, (εt) (t∈N) } additional parameters for the gradient descent output: W the indicator matrix G the reduced sparse convex embedding of X // First step: Column Subset Selection as a manifold subsampling

1 W (i=1) ← 0 // W (i) ∈ {0, 1} n×i 2 t ← arg min q k(xq, x rand[1,n] ) 3 e ← arg min q [k(xq, xt)] // xe is the first element 4 W (i=1) et ← 1 5 for i ← 2 to ℓ do 6 e ← arg min q k ⊤ xq • K -1 i-1 • kx q 7 W (i) •i ← 0; W (i) ei ← 1 8 compute K -1 i from K -1
i-1 using Eq. 9

// Second step: Sparse Projection over the defined simplex 9 for

x i ∈ X do 10 G (k=0) i• ← [1/ℓ, . . . , 1/ℓ] 11 repeat 12 G (k+1) i• ← G (k) i• -εt(G (k) i• ⊤ K ℓ -kx i ) 13 find the indices of the λ greatest values of G (k+1) i• 14 project the corresponding elements of G (k+1) i• onto ∆ λ [30]
15 set the other elements to 0

16 until ||G (k+1) i• -G (k) i• || 2 < ǫ;
of [START_REF] Thurau | Yes we can: simplex volume maximization for descriptive web-scale matrix factorization[END_REF], that can be directly adapted. Finally, let us remark that, in spite of being linear in terms of n, the number of elements in the dataset, the procedure is not linear with respect to ℓ. However, as ℓ is classically several order of magnitude smaller than n, so that it is seldom important. This is why, in a similar way, SiVM is also not linear with respect to ℓ.

The complexity of the computation of the projection matrix G is by construction linear with respect to n as in the algorithm, it is decomposed into a succession of n independent projections. However, the complexity of the projection with respect to ℓ is linear: Property 4 (Computational complexity for G -noted as Second step in Alg. 1 ) Each projection has a linear complexity regarding λ.

Proof. Each projection is computed through a gradient descent, each iteration of which has a complexity dominated by that of the GSSP (Eq. 12). The latter requires getting the λ greatest values of G t i• . To do so, the GSSP classically relies on sorting the elements of G t i• (a vector of size ℓ), with a o(ℓ log(ℓ)) complexity. However, it is possible to be more efficient by achieving this task in o(λ) thanks to the median-finding algorithm [START_REF] Maculan | A linear-time median-finding algorithm for projecting a vector on the simplex of rn[END_REF]. Finally, if q iterations are needed in the descent, then the total complexity of one projection is o(qλ).

Let us note here that q typically depends on the choice of the magnitude of the gradient step ε t , which can be efficiently set following the results of [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF]. In practice, only a few tens of iterations are necessary and if one has n >> qλ, the influence of the number of iterations is immaterial.

Finally, as both the definition of the CSS and the projection have a linear complexity with respect to the dataset size, the overall SAGA procedure also has. This allows factorizing very large matrices in tractable time.

Experiments and results

In this Section, we first observe the behavior of SAGA on some toy datasets. Our goal is to verify the behavior of our algorithm with respect to the theoretical properties given in the previous section, with a special focus on the nature of the subsampling occurring the RKHS. Then, the volume maximization strategy is discussed, as well as its impact on the reconstruction errors over a toy and a real dataset. We finally discuss the potential use of SAGA in a feature extraction context for classification purpose. The performances of our method are then compared to a selection of state-of-the-art methods performing NMF with characteristics shared by our method (sparsity, kernels, convexity, etc.). The SAGA paradigm is first illustrated on toy datasets (Fig. 2, 3 and4). In the first examples we consider simulated datasets (ring and S shaped respectively) made of 600 points in R 2 : In the input space, the points of the CSS form a nonsimplicial polytope which intuitively fits with the manifold hull idea (both inner and outer contours are displayed for the ring). If one projects points onto the corresponding CSS (ℓ = 30 for the ring, and ℓ = 50 for the S), the magnitude of the reconstruction error fits with the non-linear geometry of the simplex in the input space (Fig. 2). The sparsity and the locality of the reconstruction are displayed on Fig. 3: It appears that for each point, the number of non-null components is smaller than or equal to λ = 3 (ring) or λ = 5 (S). Moreover, these non-nil components are all located in the very close neighborhood of the projected point, as expected for a non-linear local embedding. Figure 4(a) shows the evolution of the reconstruction error and the locality of the samples used for the reconstruction. This last term is measured as the radius of the minimum volume enclosing ball (computed as a smallest enclosing ball problem [START_REF] Gärtner | Fast and robust smallest enclosing balls[END_REF]). As expected, this mean radius is decreasing as more samples are taken from the dataset. Next, the reconstruction error of SAGA is evaluated with respect to the chosen sparsity level and compared to the original version of SiVM [START_REF] Thurau | Yes we can: simplex volume maximization for descriptive web-scale matrix factorization[END_REF], yet in the Gaussian RKHS, i.e. a simplex volume maximization of Eq. 5, followed by a projection based on quadratic programming. We draw 2, 000 points in R 30 according to a Gaussian distribution with covariance matrix σ = 0.5I. We measure for SiVM and SAGA the normalized reconstruction error ||Φ -ΦWG ⊤ || 2 /||Φ|| 2 . We remark here that this formula is correct for SiVM if it is implemented through a kernel form with the linear kernel. Results are displayed in Figure 4(b). When ℓ ≤ p = 30, SAGA (with or without sparsity constraint) as well as SiVM directly operates as a dimensionality reduction method: As the SiVM curve is below that of SAGA λ = ℓ, SiVM appears as more reliable, which makes sense, as the Gaussian distribution provides a rather convex dataset. However, even on such a dataset, if the reconstruction error is considered along with sparsity, it appears that, whatever the value λ * chosen for parameter λ, it is always possible to find a value ℓ * for ℓ, such that the reconstruction error with SAGA tuned with (ℓ = ℓ * , λ = λ * ) is smaller than with SiVM tuned with (ℓ = λ * ). This illustrates well that it is possible to reduce the reconstruction error while constraining the solution sparsity. If one considers the case where ℓ ≥ p = 30, SiVM fails in producing reconstruction error which decreases when ℓ increases: From Eq. ( 1), the addition of a p+1 vertex leads to a simplex of null volume, which is impossible to maximize, and turning SiVM into a random projection method. Thus, the reconstruction error becomes greater than with SAGA, the latter enhancing the reconstruction quality, despite strong sparsity constraints.

Experiments on toy datasets

Comparison on subsampling strategies

The quality of the overall matrix factorization procedure relies on two elements: the ability of the CSS to correctly define the manifold hull, and the precision of the projection. A particular focus is given here on the first one, as various sampling strategies are compared and discussed. The goal of our first comparison is to confirm that, among the methods based on volume maximization, ours is the most accurate. To this end, we compare three volume maximization methods (operating in the RKHS, to fit the objectives of this work): The first one is very method from [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF], reported through Eq. 5, including the constant distance assumption as well as the replacement of squared distances by simple distances. As this point is only supported by computational considerations, we also consider the method summarized by Eq. 4, where the constant distance assumption still holds, but where the squared distances are kept. These two methods are referred to as SiVM-approx-(4) and SiVM-approx- [START_REF] Boyd | Convex optimization[END_REF]. Naturally, the third one is that of SAGA, based on Eq. 8 and 9. As a reference, we consider the exact volume computation based on the Cayley-Menger determinant (CMD), and we compare the ratios of the volumes derived by the three methods over the reference one. Due to the computation cost of CMD, we have restrict the size of the CSS to ℓ = 8.

Three types of datasets are used for this comparison. First datasets are composed of 2, 000 points in R 30 according to a uniform distribution. The second type of datasets is used to test the robustness of the methods to ill-conditioned data; so, 2, 000 elements in R 50 are generated and their singular values are transformed such as described in [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF]. Finally, our last dataset is obtained through a random selection of 160 images from the real dataset COIL-20 (Columbia University Image Library). COIL-20 contains 128 × 128 gray images of 20 objects at different view angles, for a total number of sample of 1, 440 [START_REF] Nene | Columbia object image library (coil-20)[END_REF]. We use 30 datasets of each type and we compute the mean ratio (and variance) of the volumes, as described above (see Table 1). Some conclusions can be drawn from Table 1. First, the lower ratios of SiVMapprox-(4) and SiVM-approx-(5) for ill-conditioned datasets confirm the conclusion of [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF] regarding the constant distance approximation, while SAGA is not Table 1 Mean ratio and variance (in %, over 30 runs) of maximum simplex volumes found by both approximations from Eq. 4 and Eq. 5, as well as SAGA over the exact CMD approach. For each dataset, the most accurate result is in bold font.

Uniform

Ill-conditioned COIL-20 SiVM-approx-( 4 bothered. Second, it is possible to notice the slight decrement of the performances due to the supplemental approximation in SiVM-approx- [START_REF] Boyd | Convex optimization[END_REF] where the squared distances are not considered anymore for computational reasons. Also, for all datasets, SAGA finds the most similar simplex volumes to the reference ones, as the ratios are the closest to 1. Oddly enough, the ratios are sometimes even slightly greater than 1, due to the numerical imprecisions of the determinant computations which may induce different choices for the simplex vertices. In the meantime, while volume differences between SAGA and CMD are low, their respective computational performances are quite different. As an example, with COIL dataset, each CMD computation requires ≈ 10 sec. whereas SAGA takes ≈ 0.01 seconds in an unoptimized implementation. As a conclusion, SAGA provides a volume as large as what of CMD approach, yet in much less time. However, the computational accuracy of SAGA is fully investigated in the next subsection. Now that it is established that SAGA provides a better hull for the manifold than classical methods based on a kernelization of the simplex volume maximization, let us compare it to other manifold sampling strategies. As it is not meaningful to use the simplex volume as a criterion, we consider the normalized reconstruction errors, such as with toy datasets. We can also consider larger CSS values, ranging from 10 to more than 100, as the computation of the CMD is not an issue anymore. In the comparison, we keep SiVM-approx- [START_REF] Boyd | Convex optimization[END_REF], as it corresponds to a kernelized version of the original paper [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]. We also consider random sampling (the NMF reducing to a random projection algorithm), kernel k-means and kernel rank-revealing Cholesky of [START_REF] Courty | A kernel view on manifold sub-sampling based on karcher variance optimization[END_REF], as a surrogate for the rank-revealing QR method of [START_REF] Shroff | Manifold precis: An annealing technique for diverse sampling of manifolds[END_REF] which cannot be kernelized.

Two datasets are used, namely the UCI Wine dataset [START_REF] Bache | UCI machine learning repository[END_REF] composed of 178 instances with 13 attributes and a simulated dataset which consist of 2, 000 points in R 50 according to a Gaussian distribution. Reconstruction error curves are presented in Fig 5 and demonstrate the best overall performances of SAGA with increasing ℓ values. However this general trend differs according to the considered dataset and the ℓ values. For example if SAGA clearly outperforms SiVM with the Wine dataset, their performances are very close with the Gaussian one. With lowest ℓ values, kernel k-means and kernel rank-revealing Cholesky get better results on the Wine dataset, however, the reconstruction error remains high whatever the strategy. both linear and Gaussian kernel. Each time the experiment is repeated 20 times to stabilize the measure, and the computation times are reported in Fig. 6(a). As expected, we observe a near linear trend with both different versions of kernel. The computational differences occurring between the two kernels are mostly due to the evaluation of the kernel: the linear kernel results in a simple dot product operation, whereas evaluating the Gaussian kernel involves computations of transcendental operators.

Computational complexity

If we now turn to the complexity of the projection, it is useless to consider it with respect to n as the projection is dealt datum by datum. However, it is interesting to consider the complexity regarding ℓ. To do so, we use the same experimental setting, yet, n is fixed and ℓ varies from 1 to 120. Fig. 6(b) clearly highlights the outperformance of our projector compared to the quadratic programming approach used in SiVM.

Application to classification on real datasets

In order to propose a more application-oriented test of SAGA, we propose to use the result of various NMF methods as feature extraction processes in a classification problem. Experiments are conducted on four publicly available datasets. The previously described COIL-20 dataset, the CMU PIE dataset, the MNIST dataset and the CIFAR-10 dataset: The CMU PIE dataset contains the 32 gray scale facial images of 68 people, any of each having 42 images at different illumination and facial expression conditions [START_REF] Sim | The cmu pose, illumination, and expression (pie) database[END_REF]. The MNIST dataset [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] contains the 28 × 28 gray scale images of handwritten numbers (from 0 to 9). Its total sample number is 70,000. Finally, CIFAR-10 is made of 60,000 32 × 32 color images in 10 classes, with 6,000 images per class.

Each dataset is partitioned into two sets: one for training (one tenth of the dataset), and the other for testing. On the training dataset X t , a factorization is conducted and leads to derive F t and G t . F t plays the role of a visual coding dictionary, and is used to reconstruct the test set. We note here that F t could have been constructed with respect to the entire dataset, in a unsupervised learning manner, but it was not the case. G t is used to train a SVM classifier (with Gaussian kernel), with parameters optimized by standard cross-validation. The testing set X s is then projected over F t which allows deriving G s , used for testing with the SVM classifier. This process is repeated 10 times for each value of ℓ ∈ {10, 20, 30, 40, 50}, and for each ℓ/λ ratio ∈ {1, 1.5, 2}, in order to stabilize the performances. In this setting, we have compared the result of SAGA to other stateof-the-art algorithms: Sparse NMF [START_REF] Kim | Sparse non-negative matrix factorizations via alternating nonnegativity-constrained least squares for microarray data analysis[END_REF] (with ℓ/λ = 2), Kernel NMF [START_REF] Li | A new kernel non-negative matrix factorization and its application in microarray data analysis[END_REF], Convex NMF [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF] and Kernel Convex NMF [START_REF] Ding | Convex and semi-nonnegative matrix factorizations[END_REF] (or KC NMF for short). We have used the same kernel variance for all kernel-based techniques, on the basis of a standard rule of thumb [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. The mean results over the 10 repetitions are displayed in Fig. 7. It shows that SAGA produces the most accurate results on three datasets out of four. In the case of CIFAR-10 dataset, one notices first the very low performances of all the methods with respect to the state-of-the-art works focusing on classification performances. The reason is that the size of the CSS remains rather low in our experimental setting: the point of this comparison is not to exhibit the highest possible accuracies, but rather to provide a sound experimental setting across various datasets of heterogeneous difficulty. However, this does not explain why SAGA does not compete with Sparse NMF and Kernel NMF on this dataset. A possible explanation stems from the complexity of the manifold hull which cannot be described efficiently by so few elements. In this particular case, SAGA would perform slightly worse than other state-of-the-art methods. Interestingly enough, whatever the dataset, when λ stays the same, the SAGA classification accuracy improves with the increasing number ℓ of simplex vertices. However, ℓ exerts little influence when λ is already large. As discussed in [START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF], we can relate the optimal value of ℓ to the intrinsic dimensionality of the manifold where the data live.

Table 2 summarizes the results of Fig. 7 by averaging the performances ∀ℓ ∈ {10, 20, 30, 40, 50}, (for Sparse NMF and SAGA, with ℓ/λ = 2): The variances which are not displayed on Fig. 7 for clarity sakes are given here. SAGA variance is sometimes important due to the strong increment of the performances when ℓ increases. To allow for more complete comparisons, we have also added a kernel version of the original SiVM with our projection method, noted Sparse Kernel SiVM (SK SiVM for short), which basically amounts to using SAGA, yet with Eq. 5 instead of Eq. 8 and 9; it appears that it is always less accurate than SAGA. Finally, on the majority of the datasets, the superiority of SAGA is established. Also, it is the fastest of all, as can be seen in Fig. 8, where the factorization performances are reported for SAGA and the considered state-of-the-art NMF methods for the whole PIE dataset when varying the size of the CSS. This illustrates the computational benefits of our approach.

Conclusion

SAGA (Sparse and Geometry Aware) is a new matrix factorization algorithm which has the following properties: (1) it operates in the Gaussian RKHS, which accounts for potential nonlinearity in the dataset geometry; (2) it provides sparse and convex projections onto a reduced embedding spanned by selected typical samples, which facilitates the human interpretation, and leads to a non-linear local representation of the data; (3) it has a complexity linear with the number of data entries, which allows dealing with big data. SAGA relies on both a manifold sampling strategy and a data projection. This latter has been proved to converge under some conditions regarding the Gaussian kernel variance. Finally, SAGA has been tested on toy, simulated and real datasets. The following conclusions can be drawn from the experiments: we observed in accordance with the theory that SAGA encodes the data as convex combinations of neighbor samples; the proposed volume maximization heuristic leads to better subsampling of the original data with respect to the volume of the simplex formed by the CSS; and the performances of its feature extraction have proved to outperform the selected state-of-the-art other NMF methods on classification tasks.

Fig. 2 A

 2 Fig. 2 A ring and an S-shape datasets (a-c) It is possible to have ℓ > p elements in the CSS forming a non-simplicial polytope (here illustrated by a 4-NN graph) which approximates well the contour of the shape. (b-d) Each pixel of the image is projected in the RKHS onto the 29-simplex (b) or 49-simplex (d). Each pixel is colored according to its reconstruction error (black for low values, white for high ones).

Fig. 3

 3 Fig. 3 (a-b) The 29-simplex on the ring shape dataset: The red point is projected in the RKHS as a 3-sparse vector. Each vertex color accounts for the corresponding convex coordinates (from blue= 0 to red= 1) (c-d) The S dataset with 5-sparse vectors (in (d), only 4 projections are non-nil).

Fig. 4

 4 Fig. 4 (a) Reconstruction error (in red) and mean radius of smallest enclosing ball (blue) for the ring dataset. (b) Normalized reconstruction error.

First, we considerFig. 5

 5 Fig. 5 Comparison of the reconstruction error according to the size of the CSS: (a) Wine dataset; (b) simulated dataset (50 dimensional Gaussian).

Fig. 6

 6 Fig. 6 Comparison of the computational times for CSS (a) and the projection (b) for SiVM and SAGA.

Fig. 7

 7 Fig. 7 Feature extraction evaluation on (a) COIL-20, (b) PIE, (c) MNIST and (d) CIFAR-10 datasets at different reduced dimension.

Fig. 8

 8 Fig.8Computational performances between SAGA and the other considered NMF methods on the whole PIE dataset.

Table 2

 2 Mean accuracy and variance NMF-based classification. For each dataset, the most accurate method is in bold font. This property is assessed thanks to a paired Student's T-test between the SAGA and the best other method (apart from SK SiVM which can be seen as a particular case of SAGA)

		COIL-20	PIE	MNIST	CIFAR-10
	Sparse NMF	88.04 (0.66)	80.24 (1.25)	86.08 (1.99)	38.34 (0.56)
	Kernel NMF	94.38 (1.75)	79.89 (7.25)	83.61 (2.86)	38.90 (1.87)
	Convex NMF	83.46 (2.84)	54.95 (9.07)	80.66 (2.80)	31.31 (0.45)
	KC NMF	85.18 (0.39)	59.80 (1.75)	81.94 (0.35)	25.76 (0.16)
	SK SiVM	96.79 (1.34)	83.73 (9.32)	89.71 (3.78)	35.50 (2.20)
	SAGA	97.23 (1.58) 84.44 (13.62) 89.92 (3.18)	35.59 (2.44)
	t-statistics	3.2327	1.8872	2.3257	X
	Confidence Level	≥ 99%	≥ 95%	≥ 97.5%	X

We write W instead of W (ℓ) when ℓ does not matter, or is implicit regarding the context.

The Karcher variance is a variance accounting for Riemannian distances over the manifold, rather than Euclidean ones.

Acknowledgements This work was partially supported by the ANR fundings ANR-10-INBS-08 (ProFI project, "Infrastructures Nationales en Biologie et Santé", "Investissements d'Avenir"), ANR-13-JS02-0005-01 (Asterix project). and the Prospectom project (Mastodons 2012 CNRS challenge).


Julien Rebetez

Devis Tuia
email: devis.tuia@epfl.ch

Nicolas Courty
email: ncourty@irisa.fr

Network-based correlated correspondence for unsupervised domain adaptation of hyperspectral satellite images

Adapting a model to changes in the data distribution is a relevant problem in machine learning and pattern recognition since such changes degrade the performances of classifiers trained on undistorted samples. This paper tackles the problem of domain adaptation in the context of hyperspectral satellite image analysis. We propose a new correlated correspondence algorithm based on network analysis. The algorithm finds a matching between two distributions, which preserves the geometrical and topological information of the corresponding graphs. We evaluate the performance of the algorithm on a shadow compensation problem in hyperspectral image analysis: the land use classification obtained with the compensated data is improved.

I. INTRODUCTION

Domain adaptation problems occur naturally in many applications of machine learning to real-world datasets [START_REF] Quiñonero-Candela | Dataset shift in machine learning, ser. Neural information processing series[END_REF]. In remote sensing image analysis this problem arises frequently, since the acquisition conditions of the images (cloud cover, acquisition angle, seasonal variations) are most often different. As a consequence, even if the images contain the same type of objects, the observed data distribution undergoes a ddimensional and often nonlinear spectral distortion, i.e. a distortion that is local, class-specific and that impacts differently each region of the electromagnetic spectrum [START_REF] Camps-Valls | Advances in hyperspectral image classification: Earth monitoring with statistical learning methods[END_REF], [START_REF] Lunga | Manifold-learningbased feature extraction for classification of hyperspectral data[END_REF]. One way to solve this problem is to perform an adaptation between the two d-dimensional image domains, in order to achieve a relative compensation of the shift by matching the data clouds to each other. Provided that the data are expressed as graphs and embed a topological structure, this problem can be seen as a graph matching problem [START_REF] Bunke | Error correcting graph matching: on the influence of the underlying cost function[END_REF]. In hyperspectral remote sensing, the adaptation problem has been tackled equivalently in [START_REF] Tuia | Graph matching for adaptation in remote sensing[END_REF], where the authors use the distance between nodes to constrain the possible assignments and enforce structure preservation of the graph. Recent alternative methodologies consider the domains as manifold to be aligned [START_REF] Wang | Manifold alignment[END_REF] and use labeled examples [START_REF] Wang | Heterogeneous domain adaptation using manifold alignment[END_REF] or local geometric matching [START_REF]Manifold alignment without correspondence[END_REF] to perform the alignment. If the first option has been successfully applied to remote sensing data in a multimodal setting [START_REF] Tuia | Semisupervised manifold alignment of multimodal remote sensing images[END_REF], local geometric matching must be handled with care, as the datasets to be matched are of large size and comparing all possible arrangements of local patterns is computationally really expensive.

To find the best matching, one can resort to energy minimization on Markov Random Fields (MRF). In the last years, this topic has received a lot of attention in computer vision (e.g. for dense stereo matching, segmentation and image stitching) [START_REF] Szeliski | A comparative study of energy minimization methods for Markov Random Fields[END_REF] and efficient algorithms can be used to solve the optimization, especially if it consists solely of unary and pairwise constraints. In hyperspectral remote sensing, MRF are widely used to enforce spatial constraints in the classifiers [START_REF] Moser | Land-cover mapping by Markov modeling of spatial-contextual information[END_REF] but, to our knowledge, the only attempt to use them for domain adaptation is found in [START_REF] Jacobs | Domain adaptation with hidden Markov Random Fields[END_REF]. In that study, the two domains are considered as observations of a single hidden MRF. Nodes in the source and target graphs are matched by finding a hidden node that is likely to have generated both. To compute unary and pairwise potentials, a Gaussian distribution of the distances between nodes and of the edges lengths is assumed. In this paper, we propose a MRF formulation of the graph matching problem. As a first contribution, we extend the algorithm in [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces[END_REF], originally developed to register a pair of 3D meshes, to the d-dimensional case by using descriptors issued from network theory [START_REF] Crucitti | Centrality measures in spatial networks of urban streets[END_REF] and secondly propose an alternative optimization to reduce the number of candidate matches and allow working on larger graphs. Contrarily to [START_REF] Jacobs | Domain adaptation with hidden Markov Random Fields[END_REF], we rely on node descriptors that are more robust to deformations between source and target domains. We also avoid the requirement of having the same number of nodes in both domains. Experiments on relative shadowing compensation show important benefits for land use pixel classification of a dataset comprising a 144-channels hyperspectral image and a LiDAR digital surface model.

II. THE CORRELATED CORRESPONDENCE ALGORITHM (CC)

Consider two graphs, one for the source data, S, and one for the target data, T . Each graph is composed of a set of nodes, V , connected by a set of edges, E. The total number of nodes, K, can differ from one graph to another. We will refer to the source graph as G S = (V S , E S ) and to the target graph as G T = (V T , E T ), with K S and K T nodes, respectively.

The correlated correspondence algorithm (CC, [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces[END_REF]) solves a three-dimensional graph matching problem that consists in finding correspondences between the nodes of both graphs. Each node in the source graph {V S k } K S k=1 is associated with a correspondence variable, c k , whose value represents the index of the node in V T to which V S k has been matched. The CC solution to the graph matching problem is given by the (K S × 1) vector c containing all the nodes assignments. CC solves this graph matching problem as an MRF energy minimization problem. Each node on the source graph is associated with a unary potential and a number of pairwise potentials. Unary potentials ψ(c k ) encode the dissimilarity between a source node and each of the target nodes. Pairwise potentials ψ(c k , c l ) encode the cost of assigning two neighboring source nodes to two specific target nodes, penalizing the assignments producing strong changes in the topology of G S . There are two such potentials: ψ n (c k , c l ), which penalizes an increase of the distance between neighboring nodes in G S and ψ f (c k , c l ), which penalizes a decrease of the distance between non-neighboring nodes in G S . CC minimizes the energy function

E(c) = k∈K S ψ(c k ) + k,l∈N S ψn(c k , c l ) + k,l∈F S ψ f (c k , c l ), (1)
where N S is the set of neighboring nodes on the source graph and F S is the set of nodes which are far from each other.

Although minimizing E(c) is in general NP-hard, efficient approximate algorithms exist such as the Tree-Reweighted Message Passing algorithm (TRWS [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF]) used in this paper and Loopy Belief Propagation [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF].

III. THE NETWORK-BASED CORRELATED CORRESPONDENCE ALGORITHM (NETCC)

In this section, we detail the proposed network-based correlated correspondence (netCC) algorithm.

A. Node descriptors

The first step necessary to compute all the potentials in Eq. ( 1) is to obtain a set of descriptors used to match the nodes. These descriptors must reflect proximity for nodes to be matched and dissimilarity for nodes not to be matched. If the graphs are relatively similar in shape and position in the input space (as in [START_REF] Tuia | Graph matching for adaptation in remote sensing[END_REF]), the Euclidean distance between the d-dimensional vectors can be used. However, when strong distortions (scaling, local transformations) are observed, like in Fig. 2, this distance is likely to produce too many mismatches, as the nodes in the source graph would be matched to their direct neighbors. Authors in [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces[END_REF] use local surface signatures based on local 2D histograms. These local histograms are obtained after projecting the neighboring points onto the tangent plane associated with the point. Although this works well for 3D meshes, it is not straightforward to extend this method to the d-dimensional case (where we do not have surface information) and we therefore propose an alternative method.

In the proposed netCC, we compute a vector of descriptors, which includes several node characteristics invariant to rotation and translation of the graph (Fig. 1) :

• Closeness centrality [START_REF] Crucitti | Centrality measures in spatial networks of urban streets[END_REF], V c i : how far (in term of shortest path) a node i is from all the other nodes.

• Eccentricity, V e i : the maximum shortest path distance between node i and all other nodes.

• Median distance to the manifold, V m : the median Euclidean distance between node i and all the points in the domain. • Density, V d : the average Euclidean distance between node i and its neighbors in the domain.

• Distance to center, V dc : the Euclidean distance between node i and the graph's gravity center. To obtain a descriptor vector for each node, these measures are stacked as

x i = [V c i V e i V m i V d i V dc i ] ⊤ .

B. Unary potentials

Each node in the source graph {V S k } K S k=1 yields an unary potential that encodes the cost of assigning it to each node in the target graph

{V T i } K T i=1 : ψ(c k = i) = d(x S k , x T i ), (2)
where d is the Euclidean distance between the vectors of descriptors of the two nodes.

C. Pairwise potentials

For the pairwise potentials, we follow [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces[END_REF] and define potentials that favor assignments preserving distances between nodes of the source graph G S . Such distances are computed as the shortest path distances on the graph between the two nodes, d sp (V 1 , V 2 ). We also normalize all the distances on both graphs by the maximum distance on each graph in order to be able to compare the distances between pairs of nodes computed on G S with those on G T . Such a normalization only deals with global scaling problems (as those considered in the experiments presented), but should be modified for nonuniform local scalings.

Nearness preservation potentials: The nearness preservation potentials ψ n favor assignments in which two neighboring nodes V S k and V S l on G S will stay close after matching with G T . Two nodes V T i and V T j on G T are considered close if their shortest path distance is lower than a threshold α.

ψn(c k = i, c l = j) = 0 if dsp(V T i , V T j ) < α 1 otherwise (3)
The threshold α is to be chosen by the analyst. Since its value corresponds to the radius of a local neighborhood sphere, it should remain quite small (e.g. α ∈ [0.15...0.35], in order to avoid considering all the nodes as neighbors.

Farness preservation potentials: The farness preservation potential ψ f encodes the constraint that two nodes, which are far on G S , should be assigned to nodes that are far on G T . Two nodes are considered far if their shortest path distance is greater than a threshold β. In principle, we admit that both farness and nearness potentials are considered for disjoint set of pairs, i.e. β > α.

ψ f (c k = i, c l = j) = 0 if dsp(V T i , V T j ) > β 1 otherwise (4)

D. Computational complexity

We apply a number of optimizations to reduce the computational cost of the minimization of Eq.( 1). One shortcoming of this formulation is that each node V S k of the source graph can be matched to any target node V T i . This means that unary potentials will be stored as vectors of length K T and pairwise potentials as matrices of size (K T × K T ). Since we have K S such potentials, this leads to a O(K S K T 2 ) memory requirement for each pairwise potential.

To reduce the number of farness potentials, we follow the strategy in [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces[END_REF]. We first solve the energy minimization problem of Eq. ( 1) without the farness potentials and we then check which potentials are violated by the solution. These potentials are added to the objective function and the optimization is iterated until no additional potentials are required.

We also reduce the number of candidate matches among the target nodes. We do so by considering only the m closest target nodes in descriptor space. In this work, we empirically choose m = 30, which seems to be a good tradeoff between the number of missed matches and the computational gain. By using only m candidate target nodes, the memory requirement for pairwise potentials reduces to O(K S m 2 ). In our tests, increasing this value did not yield any performance improvement, but this might depend on the dataset.

E. Full algorithm

Algorithm 1 contains the detailed steps to perform domain adaptation using netCC. Creating a graph where each point in the domain is a node would result in a very large graph that would be impractical to work with. Instead, we use vector quantization and build a graph on the centroids found by the clustering algorithm for each domain. We use two levels of clustering :

• A fine level with C clusters. We use these graphs to approximate geodesic distance between two data points on the manifold [START_REF] Bernstein | Graph approximations to geodesics on embedded manifolds[END_REF]. • A coarse level where we subsample the fine level to K S < C, respectively K T < C clusters. The resulting graphs are used in the netCC algorithm, but all distances are computed on the fine graph. The algorithm defines a local transformation that can be used to transfer points between the two domains. Therefore, given a d-dimensional source dataset, netCC creates an adapted d-dimensional dataset composed of each source point transferred to the target domain. A classifier can then be trained on the target dataset and applied directly to the adapted source points.

IV. DATA AND SETUP OF EXPERIMENTS

A. Dataset

We use the dataset of the 2013 IEEE GRSS Data Fusion Contest [START_REF] Pacifici | Report on the 2013 IEEE GRSS data fusion contest: Fusion of hyperspectral and LiDAR data[END_REF] which consists of an hyperspectral image and a Digital Surface Model (DSM) derived from airborne LiDAR at a 2.5m spatial resolution with an image size of 1903 × 329. The hyperspectral image has 144 spectral bands. The goal is to classify each pixel into a land use class using its spectral Algorithm 1 Domain adaptation with network-based correlated correspondence Inputs -Source dataset X S ∈ R d -Target dataset X T ∈ R d 1: Eigendecompose X T 2: Project both datasets on the npc top eigenvectors 3: Cluster each dataset using k-means with C clusters 4: Build a kNN graph per dataset using the cluster centers as nodes 5: Compute geodesic distances d sp within each graph 6: Subsample graph G S to keep K S nodes 7: Subsample graph G T to keep K T nodes 8: Apply the netCC algorithm to find the optimal correspondence vector c. 9: for each node k ∈ K S : do

10: Compute the displacement D k = V T c k -V S k 11:
Apply the displacement to all the points in X S assigned to node V S k 12: end for 13: Reconstruction : invert the projection to retrieve the ddimensional data. signature and its altitude on the DSM. Labels are provided for around 15000 pixels as shown in table I. The shadow of a cloud affects the right part of the image, leading to important shifts in the spectral signatures of the hyperspectral image.

Figure 2 shows the lit and shadowed domains on the same scale. Thus, adaptation is required to allow a single classifier to deal correctly with both shadowed and lit areas. We obtained a cloud mask (fig. 3c) by thresholding the magnitude of the bands corresponding to RGB and Near-Infrared and using manual intervention to remove some of the spurious points.

We can therefore divide our image in two domains : the source domain contains all the pixels affected by the cloud shadow while the target domain is the remainder of the image.

B. Experimental setup

The labels provided for this dataset are divided into two groups : one for training and one for testing. The training set only contains pixels that are not affected by the shadow while the testing set contains both lit and shadowed pixels.

To evaluate the performance of our adaptation procedure, we use a one-vs-one SVM with a Gaussian kernel. After model selection by cross-validation in the ranges γ = [10 -3 , ..., 10 2 ] and C SVM = [10 0 , ..., 10 3 ], we obtained the optimal SVM parameters. For each pixel, the following features are used by the SVM : The value of C depends on the dataset and should be set large enough to yield a good representation of data. Since the fine graphs are used only for the computation of the geodesic distances, which are computed only once at the beginning of the process, C can be increased without impacting the performances.

We compare the proposed netCC algorithm with histogram matching applied on all the spectral bands. Since we are interested in the classification improvement after adaptation, we consider two separate testing scenarios : testing pixels under the shadow only and all testing pixels. In the experiments, we study the variability with regard to 2 parameters : the dimensionality of the space where the matching is performed (npc) and the number of nodes in the coarse graphs (K S , K T ). For each choice of parameters, we run the whole graph matching procedure (clustering, matching) 10 times and report the average performances.

V. RESULTS AND DISCUSSIONS

With the setting described above, the proposed netCC algorithm converges in a limited number of iteration, typically around 100. This corresponds to an average running time of a few minutes on a laptop (Intel i5).

Table II shows the results when using only pixels under the shadow as test samples (around 1400 points). Without adaptation, we get a κ score of 0.131, which is very low. We can see that by using histogram matching or netCC, we respectively reach a performance of 0.36-0.47, which is acceptable considering that this is a problem involving 15 classes, and some of them are very similar (e.g. there are three types of grass, two types of parking lots, two types of buildings). We obtain an improvement of 0.01 to 0.11 κ points by using netCC over histogram matching. This corresponds to an increase in accuracy of about 11%. The choice of the number of centroids, K S and K T , has an important influence on the results. This is of course dependent on the problem, but on our testing dataset, it seems K S = K T = 50 is enough to get a good representation of the data manifold and increasing to 100 mostly increases the number of clusters near the center of the graph. This makes it harder to find a set of correspondences that satisfy all nearness potentials. On the other hand, we see that increasing the number of principal components (npc) has a positive effect on the performances. This is because with less principal components the reconstruction will be of lower quality.

Table III reports the results obtained when using of all the testing labels (about 12000 points). We get a small increase in performance when using netCC over histogram matching. The effect of DA is less visible due to the large number of test points in the lit area that do not benefit from adaptation.

Figure 4 shows classification maps for the right part of the image with the tested methods. Visually, the benefits of domain adaptation are obvious. We can see that without domain adaptation, the classifier fails to produce a meaningful result and classifies everything as water. The problem disappears when adapting the spectra prior to classification.

Figure 5 provides a magnified view of two areas where histogram matching and netCC result in different classifications. In the second row, we can see that histogram matching causes a grass area (visible on the labels and the flat lidar) to be classified as trees while netCC correctly classify it as grass. In the third row, a zoom on an urban area is shown. From the RGB and lidar images, it is clear this area contains buildings, some of which are misclassified as grass after histogram matching. With netCC, they are correctly classified as buildings, albeit not of the correct class. The same observation can be made for some roads which are misclassified as soil after histogram matching.

VI. CONCLUSION

In this paper, we propose an extension of the correlated correspondence algorithm for graph matching to handle d dimensional datasets. We propose the use of descriptors used in network analysis to guide the matching process. To reduce the computational costs, we use a two-levels graph and only allow matches between neighbors in the descriptor space. Experiments on a challenging real-world dataset show good performances and a significant improvement over histogram matching. Future research directions could focus on the choice of descriptors and on the use of different potentials. Deeper exploration and sensitivity analysis of the various parameters would also be interesting, as well as the application to other datasets.

Another area for improvement is to study better ways to transport the points for a given transformation of the nodes.

In our current implementation, the local transformation is a simple translation, but more generic classes of transformation (e.g. affine) could be considered and possibly improve the performance of the classifier trained in the adapted domain.

Fig. 1 :

 1 Fig. 1: Node descriptors for the G S (top) and G T (bottom) graphs.

  (a) RGB (three channels among the 144 available).

Fig. 3 :

 3 Fig. 3: The 2013 GRSS DFC dataset.

•

  144 spectral bands • 1 LiDAR band • 6 bands obtained by applying opening by reconstruction and closing by reconstruction on the LiDAR band with three different square structuring elements of sizes : {7, 19, 31}. Morphological filters are often used in remote sensing pixel classification to enforce the smoothness of the decision function in the spatial domain [19]. The spectral bands are normalized by their overall maximal value and the lidar band is normalized by its maximal value. We set the nearness threshold α = 0.2 and the farness threshold β = 0.3. We use C = 200 for the fine graphs and experiment with two different values for K S = K T = [50, 100].

Fig. 2 :

 2 Fig. 2: Projection of the two domains on the two first principal components. Colors indicate class membership (black is unlabeled).

Fig. 4 :Fig. 5 :

 45 Fig. 4: Comparison of classification maps for the three methods on the part of the image affected by the shadow of the cloud.

TABLE I :

 I Labeled pixels in the dataset

TABLE II :

 II Results with the shadowed test points

			Overall Accuracy (%)	Kappa (κ)
	npc	K S , K T	µ	σ	µ	σ
	10	100	53.3	4.1	0.471 0.044
		50	53.9	4.4	0.476	0.047
	5	100	44.4	4.3	0.375 0.046
		50	49.1	3.2	0.424 0.035
	2	100	43.6	3.3	0.368 0.036
		50	48.0	7.4	0.412 0.078
	Histogram matching 42.2	0.0	0.360 0.000
	No adaptation	18.6	0.0	0.131 0.000

TABLE III :

 III Results with all the test points

			Overall Accuracy (%)	Kappa (κ)
	npc	K S , K T	µ	σ	µ	σ
	10	100	89.3	0.5	0.884 0.005
		50	89.4	0.5	0.885	0.006
	5	100	88.3	0.5	0.873 0.005
		50	88.8	0.4	0.879 0.004
	2	100	88.2	0.4	0.872 0.004
		50	88.7	0.9	0.877 0.009
	Histogram matching 88.0	0.0	0.870 0.000
	No adaptation	85.2	0.0	0.841 0.000

ACKNOWLEDGMENT

This work has been partially funded by the Swiss NSF (grant PZ00P2-136827) and by the Belgian Science Office (BELSPO, grant SR/00/154). We acknowledge Prof. Prasad and the IEEE IADF TC for providing the 2013 Data Fusion Contest Data.