Kuan Liu
email: kuanl@usc.edu.

Aurélien Bellet
email: aurelien.bellet@inria.fr

Escaping the Curse of Dimensionality in Similarity Learning: Efficient Frank-Wolfe Algorithm and Generalization Bounds

Similarity and metric learning provides a principled approach to construct a task-specific similarity from weakly supervised data. However, these methods are subject to the curse of dimensionality: as the number of features grows large, poor generalization is to be expected and training becomes intractable due to high computational and memory costs. In this paper, we propose a similarity learning method that can efficiently deal with high-dimensional sparse data. This is achieved through a parameterization of similarity functions by convex combinations of sparse rank-one matrices, together with the use of a greedy approximate Frank-Wolfe algorithm which provides an efficient way to control the number of active features. We show that the convergence rate of the algorithm, as well as its time and memory complexity, are independent of the data dimension. We further provide a theoretical justification of our modeling choices through an analysis of the generalization error, which depends logarithmically on the sparsity of the solution rather than on the number of features. Our experiments on datasets with up to one million features demonstrate the ability of our approach to generalize well despite the high dimensionality as well as its superiority compared to several competing methods.

Introduction

High-dimensional and sparse data are commonly encountered in many applications of machine learning, such as computer vision, bioinformatics, text mining and behavioral targeting. To classify, cluster or rank data points, it is important to be able to compute semantically meaningful similarities between them. However, defining an appropriate similarity measure for a given task is often difficult as only a small and unknown subset of all features are actually relevant. For instance, in drug discovery studies, chemical compounds are typically represented by a large number of sparse features describing their 2D and 3D properties, and only a few of them play in role in determining whether the compound will bind to a particular target receptor [START_REF] Leach | An Introduction to Chemoinformatics[END_REF]. In text classification and clustering, a document is often represented as a sparse bag of words, and only a small subset of the dictionary is generally useful to discriminate between documents about different topics. Another example is targeted advertising, where ads are selected based on fine-grained user history [START_REF] Chen | Large-scale behavioral targeting[END_REF].

Similarity and metric learning (Bellet et al., 2015) offers principled approaches to construct a taskspecific similarity measure by learning it from weakly supervised data, and has been used in many application domains. The main theme in these methods is to learn the parameters of a similarity (or distance) function such that it agrees with task-specific similarity judgments (e.g., of the form "data point x should algorithm and its analysis. Generalization bounds are established in Section 4. Finally, Section 5 describes our experimental results, and we conclude in Section 6.

Background and Related Work

In this section, we review some background and related work in metric and similarity learning (Section 2.1) and the Frank-Wolfe algorithm (Section 2.2).

Metric and Similarity Learning

Metric and similarity learning has attracted a lot of interest over the past ten years. The great majority of work has focused on learning either a Mahalanobis distance d M (x, x ) = (x-x ) T M (x-x ) where M is a symmetric positive semi-definite (PSD) matrix, or a bilinear similarity S M (x, x ) = x T M x where M is often taken to be an arbitrary d × d matrix. A comprehensive survey of existing approaches can be found in [START_REF] Bellet | A Survey on Metric Learning for Feature Vectors and Structured Data[END_REF]. We focus below on the two topics most relevant to our work: (i) efficient algorithms for the high-dimensional setting, and (ii) the derivation of generalization guarantees for metric and similarity learning.

Metric learning in high dimensions. Both Mahalanobis distance metric learning and bilinear similarity learning require estimating O(d 2 ) parameters, which is undesirable in the high-dimensional setting for the reasons mentioned earlier. In practice, it is thus customary to resort to dimensionality reduction (such as PCA, SVD or random projections) to preprocess the data when it has more than a few hundred dimensions (see e.g., [START_REF] Davis | Information-theoretic metric learning[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Guillaumin | Is that you? Metric learning approaches for face identification[END_REF][START_REF] Ying | Distance Metric Learning with Eigenvalue Optimization[END_REF][START_REF] Wang | Parametric Local Metric Learning for Nearest Neighbor Classification[END_REF][START_REF] Lim | Robust Structural Metric Learning[END_REF][START_REF] Qi Qian | An Integrated Framework for High Dimensional Distance Metric Learning and Its Application to Fine-Grained Visual Categorization[END_REF]Liu et al., 2015b;[START_REF] Yao | High-dimensional Similarity Learning via Dual-sparse Random Projection[END_REF]. Although this strategy can be justified formally in some cases (Liu et al., 2015b;[START_REF] Qi Qian | Towards making high dimensional distance metric learning practical[END_REF], the projection may intertwine useful features and irrelevant/noisy ones and thus hurt the performance of the resulting similarity function. It also makes it hard to interpret and use for data exploration, preventing the discovery of knowledge that can be valuable to domain experts.

There have been very few satisfactory solutions to this essential limitation. The most drastic strategy is to learn a diagonal matrix M [START_REF] Schultz | Learning a Distance Metric from Relative Comparisons[END_REF][START_REF] Gao | SOML: Sparse Online Metric Learning with Application to Image Retrieval[END_REF], which is very restrictive as it amounts to a simple weighting of the features. Instead, some approaches assume an explicit low-rank decomposition M = L T L and learn L ∈ R r×d in order to reduce the number of parameters [START_REF] Goldberger | Neighbourhood Components Analysis[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Kedem | Non-linear Metric Learning[END_REF]. This results in nonconvex formulations with many local optima [START_REF] Kulis | Metric Learning: A Survey[END_REF], and requires to tune r carefully. Moreover, the training complexity still depends on d and can thus remain quite large. Another direction is to learn M as a combination of rank-one matrices. In particular, [START_REF] Shi | Sparse Compositional Metric Learning[END_REF] generate a set of rank-one matrices from the training data and then learn a metric as a sparse combination. However, as the dimension increases, a larger dictionary is needed and can be expensive to generate. Some other work has studied sparse and/or low-rank regularization to reduce overfitting in high dimensions [START_REF] Rosales | Learning Sparse Metrics via Linear Programming[END_REF][START_REF] Qi | An Efficient Sparse Metric Learning in High-Dimensional Space via l1-Penalized Log-Determinant Regularization[END_REF][START_REF] Ying | Sparse Metric Learning via Smooth Optimization[END_REF] but this does not in itself reduce the training complexity of the algorithm. [START_REF] Zhang | Efficient Stochastic Optimization for Low-Rank Distance Metric Learning[END_REF] proposed a stochastic gradient descent solver together with low-rank regularization in an attempt to keep the intermediate solutions low-rank. The complexity per iteration of their approach is linear in d but cubic in the rank of the current solution, which quickly becomes intractable unless the regularization is very strong.

Finally, some greedy algorithms for metric learning have been proposed in the literature to guarantee a tighter bound on the rank of intermediate solutions. [START_REF] Atzmon | Learning sparse metrics, one feature at a time[END_REF] use a block coordinate descent algorithm to update the metric one feature at a time. [START_REF] Shen | Positive Semidefinite Metric Learning Using Boosting-like Algorithms[END_REF] selects rank-one updates in a boosting manner, while DML-eig [START_REF] Ying | Distance Metric Learning with Eigenvalue Optimization[END_REF] and its extension DML-ρ (Cao et al., 2012b) rely on a greedy Frank-Wolfe algorithm to optimize over the set of PSD matrices with unit trace. However, these greedy methods still suffer from a computational cost of O(d 2 ) per iteration and are thus unsuitable for the highdimensional setting we consider in this work. In contrast, we will propose an algorithm which is linear in the number of nonzero features and can thus be efficiently applied to high-dimensional sparse data.

Generalization bounds for metric learning. The derivation of generalization guarantees for metric and similarity learning has been investigated in the supervised setting, where the metric or similarity is learned from a labeled dataset of n points by (regularized) empirical risk minimization. For a given family of loss functions, the results generally bound the maximal deviation between the expected risk (where the expectation is taken over the unknown data distribution) and the empirical risk of the learned metric.1 These bounds are generally of order O(1/ √ n). Several technical tools have been used to address the challenge of learning from dependent pairs/triplets, leading to different trade-offs in terms of tightness, generality, and dependence on the feature dimension d. The results of [START_REF] Jin | Regularized Distance Metric Learning: Theory and Algorithm[END_REF] apply only under Frobenius norm regularization of M and have a

√

d factor in the rate. Using an adaptation of algorithmic robustness, Bellet and Habrard (2015) obtain bounds which hold also for sparsity-inducing regularizers but with a covering number term that can be exponential in the dimension. [START_REF] Bian | Learning a Distance Metric by Empirical Loss Minimization[END_REF] rely on assumptions on the data distribution and do not show an explicit dependence on the dimension. Cao et al. (2012a) derive bounds based on Rademacher complexity and maximal deviation results for U -statistics [START_REF] Clémenc ¸on | Ranking and Empirical Minimization of Ustatistics[END_REF]. Depending on the regularization used, the dependence on the dimension d ranges from logarithmic to linear. [START_REF] Verma | Sample complexity of learning mahalanobis distance metrics[END_REF] show that the √ d factor of [START_REF] Jin | Regularized Distance Metric Learning: Theory and Algorithm[END_REF] is in fact unavoidable in the worst case without some form of regularization (or restriction of the hypothesis class). They derive bounds which do not depend on the dimension d but on the Frobenius norm of the optimal parameter M . Note however that their analysis assumes that the metrics are learned from a set of i.i.d. pairs or triplets, which is rarely seen in practice.

In all the above work, generalization in metric learning is studied independently of the algorithm used to solve the empirical risk minimization problem, and none of the bounds are adaptive to the actual sparsity of the solution. In contrast, we will show that one can use early stopping in our algorithm to control the complexity of the hypothesis class so as to make the bounds independent of the dimension d, effectively balancing the empirical (optimization) error and the generalization error.

Frank-Wolfe Algorithms

The Frank-Wolfe (FW) algorithm was originally introduced by [START_REF] Frank | An algorithm for quadratic programming[END_REF] and further generalized by [START_REF] Kenneth | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF] and [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF]. FW aims at solving constrained optimization problems of the following general form: min

M ∈D f (M ), (1)
where f is a convex and continuously differentiable function, and the feasible domain D is a convex and compact subset of some Hilbert space equipped with inner product •, • .

Starting from a feasible initial point M (0) ∈ D, the standard FW algorithm iterates over the following steps. First, it finds the feasible point S (k) ∈ D which minimizes the linearization of f at the current point Algorithm 1 Standard Frank-Wolfe algorithm

Input: Initial point M (0) ∈ D for k = 0, 1, 2, . . . do S (k) ← arg min S∈D S, ∇f (M (k) ) γ (k) ← 2 k+2 (or determined by line search) M (k+1) ← (1 -γ (k) )M (k) + γ (k) S (k) end for M (k) : S (k) ∈ arg min S∈D S, ∇f (M (k) ) . (2)
The next iterate M (k+1) is then constructed as a convex combination of M (k) and S (k) , where the relative weight of each component is given by a step size γ (k) . The step size can be decreasing with the iteration number k or set by line search. The overall algorithm is summarized in Algorithm 1. FW is guaranteed to converge to an optimal solution of (1) at rate O(1/k), see for instance [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF] for a generic and concise proof. Unlike projected gradient, FW is a projection-free algorithm: each iterate M (k) is feasible by construction since it is a convex combination of elements of D. Instead of computing projections onto the feasible domain D, FW solves the linear optimization subproblem (2). The linearity of the objective (2) implies that a solution S (k) always lies at an extremal point of D. This leads to the interpretation of FW as a greedy algorithm which adds an extremal point to the current solution at each iteration [START_REF] Kenneth | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF]. In other words, M (k) can be written as a sparse convex combination of extremal points: k) , where

M (k) = S (k) ∈S (k) α (k) S (k) S (
S (k) ∈S (k) α (k) S (k) = 1 and α (k) S (k) ≥ 0, (3)
where S (k) denotes the set of "active" extremal points that have been added up to iteration k. When the extremal points of D have specific structure (such as sparsity, or low-rankness), this structure can be leveraged to compute a solution of (2) much more efficiently than the projection operator, see [START_REF] Jaggi | Sparse Convex Optimization Methods for Machine Learning[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF] for compelling examples.

A drawback of the standard FW algorithm is that "removing" an extremal point S (k) from the current iterate (or significantly reducing its weight α (k) S (k) ) can only be done indirectly by adding (increasing the weight of) other extremal points. The variant of FW with away steps [START_REF] Guélat | Some comments on Wolfe's away step[END_REF] addresses this issue by allowing the algorithm to choose between adding a new extremal point (forward step) or reducing the weight of an existing one (away step), as shown in Algorithm 2. This can lead to sparser solutions [START_REF] Guélat | Some comments on Wolfe's away step[END_REF][START_REF] Kenneth | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF][START_REF] Jaggi | Sparse Convex Optimization Methods for Machine Learning[END_REF] and faster convergence in some cases [START_REF] Guélat | Some comments on Wolfe's away step[END_REF][START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF].

In the present work, we will introduce a FW algorithm with away steps to efficiently perform similarity learning for high-dimensional sparse data. One of our key ingredients will be the design of a feasible domain with appropriate sparsity structure.

Proposed Approach

This section introduces HDSL (High-Dimensional Similarity Learning), the approach proposed in this paper. We first describe our problem formulation (Section 3.1), then derive and analyze an efficient FW algorithm Algorithm 2 Frank-Wolfe algorithm with away steps

Input: Initial point M (0) ∈ D for k = 0, 1, 2, . . . do S (k) F ← arg min S∈D S, ∇f (M (k) ) , D (k) F = S (k) F -M (k) // forward direction S (k) A ← arg max S∈S (k) S, ∇f (M (k) ) , D (k)
A = M (k) -S (k) A // away direction if D (k) F , ∇f (M (k) ) ≤ D (k) A , ∇f (M (k) ) then D (k) ← D (k) F and γ max ← 1 // choose forward step else D (k) ← D (k)
A and γ max ← α

(k) S (k) A /(1 -α (k) S (k) A ) // choose away step end if γ (k) ← 2 k+2 (or determined by line search) M (k+1) ← M (k) + γ (k) D (k)
end for to solve it in Section 3.2.

Problem Formulation

In this work, our goal is to learn a similarity function for high-dimensional sparse data. We assume the data points lie in some space X ⊆ R d , where d is large (d > 10 4 ) but points are s-sparse on average (s d). In other words, their number of nonzero entries is typically much smaller than d. We focus on learning a similarity function S M : X × X → R of the form

S M (x, x ) = x T M x = xx , M ,
where M ∈ R d×d and •, • denotes the Frobenius inner product. Notice that for any M , S M can be computed in O(s 2 ) time on average if data points are stored in a sparse format.

Feasible domain. We will derive an algorithm to learn a very sparse M with time and memory requirements that depend on s but not on d. To this end, given a scale λ > 0 which will play the role of a regularization parameter, we parameterize M as a convex combination of rank-one, 4-sparse d × d bases:

M ∈ D λ = conv(B λ ), with B λ = ij P (ij) λ , N (ij) λ ,
where for any pair of features i, j ∈ {1, . . . , d}, i = j,

P (ij) λ = λ(e i + e j )(e i + e j ) T =   • • • • • • λ • λ • • • • • • • λ • λ • • • • • •   , N (ij) λ = λ(e i -e j )(e i -e j ) T =   • • • • • • λ • -λ • • • • • • • -λ • λ • • • • • •   .
The use of such sparse matrices was first suggested by [START_REF] Jaggi | Sparse Convex Optimization Methods for Machine Learning[END_REF]. Besides the fact that they are instrumental to the efficiency of our algorithm (see Section 3.2), we give some additional motivation for their use in the context of similarity learning.

First, any M ∈ D λ is a convex combination of symmetric PSD matrices and is thus also symmetric PSD. Unlike many metric learning algorithms, we thus avoid the O(d 3 ) cost of projecting onto the PSD cone. Constraining M to be symmetric PSD provides useful regularization to prevent overfitting [START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF] and ensures that S M can be interpreted as a dot product after a linear transformation of the inputs:

S M (x, x ) = x T M x = (Lx) T (Lx ),
where M = LL T with L ∈ R d×k . Because the bases in B λ are rank-one, the dimensionality k of the transformed space is at most the number of bases composing M . Second, each basis operates on two features only. In particular, S P (ij)

λ (x, x ) = λ(x i x i + x j x j + x i x j +
x j x i ) assigns a higher similarity score when feature i appears jointly in x and x (likewise for j), as well as when feature i in x and feature j in y co-occur (and vice versa). Conversely, S N (ij) λ penalizes the crossoccurrences of features i and j. In the context of text data represented as bags-of-words (or other count data), the semantic behind the bases in B λ is quite natural: they can be intuitively thought of as encoding the fact that a term i or j present in both documents makes them more similar, and that two terms i and j are associated with the same/different class or topic.

Optimizing over the convex hull D λ of B λ will allow us to easily control the number of active features, thereby learning a very compact representation with efficient similarity computations.

Optimization problem. We now describe the optimization problem to learn the similarity parameters. Following previous work (see for instance [START_REF] Schultz | Learning a Distance Metric from Relative Comparisons[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF], our training data consists of weak supervision in the form of triplet constraints: T = {x t should be more similar to y t than to z t } T t=1 .

Such constraints can be built from a labeled training sample (see Section 4), provided directly by domain experts or crowdsourcing campaign, or obtained through implicit feedback such as clicks on search engine results. For notational convenience, we denote A t = x t (y t -z t ) T ∈ R d×d for each constraint t = 1, . . . , T so that we can concisely write S M (x t , y t ) -S M (x t , z t ) = A t , M . We measure the degree of violation of each constraint t with the smoothed hinge loss : R → R + defined as

A t , M =    0 if A t , M ≥ 1 1 2 -A t , M if A t , M ≤ 0 1 2 1 -A t , M 2 otherwise .
This convex loss is a continuously differentiable version of the standard hinge loss which tries to enforce a margin constraint of the form S M (x t , y t ) ≥ S M (x t , z t ) + 1. When this constraint is satisfied, the value of the loss is zero. On the other hand, when the margin is negative, i.e. S M (x t , y t ) ≤ S M (x t , z t ), the penalty is linear in the margin violation. A quadratic interpolation is used to bridge between these two cases to ensure that the loss is differentiable everywhere.

Remark 1 (Choice of loss). One may use any other convex and continuously differentiable loss function in our framework, such as the squared hinge loss, the logistic loss or the exponential loss.

Algorithm 3 Frank Wolfe algorithm for problem (4)

1: initialize M (0) to an arbitrary B ∈ B λ 2: for k = 0, 1, 2, . . . do 3:

B (k) F ← arg min B∈B λ B, ∇f (M (k) ) , D (k) F ← B (k) F -M (k)
// forward dir.

4:

B (k) A ← arg max B∈S (k) B, ∇f (M (k) ) , D (k) A ← M (k) -B (k) A
// away dir.

5:

if

D (k) F , ∇f (M (k) ) ≤ D (k) A , ∇f (M (k) ) then 6: D (k) ← D (k)
F and γ max ← 1 // choose forward step 7:

else 8:

D (k) ← D (k) A and γ max ← α (k) B (k) A /(1 -α (k) B (k) A ) // choose away step 9:
end if 10:

γ (k) ← arg min γ∈[0,γmax] f (M (k) + γD (k) )
// perform line search 11:

M (k+1) ← M (k) + γ (k) D (k)
// update iterate towards direction 12: end for Given λ > 0, our similarity learning formulation aims at finding the matrix M ∈ D λ that minimizes the average margin penalty (as measured by ) over the triplet constraints in T :

min M ∈R d×d f (M ) = 1 T T t=1 A t , M s.t. M ∈ D λ . (4)
Due to the convexity of the smoothed hinge loss, (4) involves minimizing a convex function over the convex domain D λ . Note that the gradient of the objective is given by

∇f (M ) = 1 T T t=1 G t (M ), with G t (M ) =    0 if A t , M ≥ 1 -A t if A t , M ≤ 0 A t , M -1 A t otherwise . (5)
In the next section, we propose a greedy algorithm to efficiently find sparse approximate solutions to this problem.

Algorithm

Exact Frank-Wolfe Algorithm

We propose to use a Frank-Wolfe algorithm with away steps (see Section 2.2) to learn the similarity. We will exploit the fact that in our formulation (4), the extremal points (vertices) of the feasible domain D λ are the elements of B λ and have special structure. Our algorithm is shown in Algorithm 3. During the course of the algorithm, we explicitly maintain a representation of each iterate M (k) as a convex combination of basis elements as previously discussed in Section 2.2:

M (k) = B∈B λ α (k) B B,
where

B∈B λ α (k) B = 1 and α (k) B ≥ 0.
We denote the set of active basis elements in

M (k) as S (k) = {B ∈ B λ : α (k) B > 0}.
The algorithm goes as follows. We initialize M (0) to a random basis element. Then, at each iteration, we greedily choose between moving towards a (possibly) new basis (forward step) or reducing the weight of an active one (away step). The extent of the step is determined by line search. As a result, Algorithm 3 adds only one basis (at most 2 new features) at each iteration, which provides a convenient way to control the number of active features and maintains a compact representation of M (k) for a memory cost of O(k). Furthermore, away steps provide a way to reduce the importance of a potentially "bad" basis element added at an earlier iteration (or even remove it completely when γ (k) = γ max ). Recall that throughout the execution of the FW algorithm, all iterates M (k) remain convex combinations of basis elements and are thus feasible. The following proposition shows that the iterates of Algorithm 3 converge to an optimal solution of (4) with a rate of O(1/k).

Proposition 1. Let λ > 0, M * be an optimal solution to (4) and

L = 1 T T t=1 A t 2 F . At any iteration k ≥ 1 of Algorithm 3, the iterate M (k) ∈ D λ satisfies f (M (k) ) -f (M * ) ≤ 16Lλ 2 /(k + 2).
Furthermore, it has at most rank k + 1 with 4(k + 1) nonzero entries, and uses at most 2(k + 1) distinct features.

Proof. We first show that ∇f is L-Lipschitz continuous on D λ with respect to the Frobenius norm, i.e. for any

M 1 , M 2 ∈ D λ , ∇f (M 1 ) -∇f (M 2 ) F ≤ L M 1 -M 2 F (6)
for some L ≥ 0. Note that

∇f (M 1 ) -∇f (M 2 ) F = 1 T T t=1 G t (M 1 ) - 1 T T t=1 G t (M 2 ) F ≤ 1 T T t=1 G t (M 1 ) -G t (M 2 ) F .
Let t ∈ {1, . . . , T }. We will now bound

∆ t = G t (M 1 ) -G t (M 2 ) F for any M 1 , M 2 ∈ D λ .
The form of the gradient (5) requires to consider several cases:

(i) If A t , M 1 ≥ 1 and A t , M 2 ≥ 1, we have ∆ t = 0. (ii) If A t , M 1 ≤ 0 and A t , M 2 ≤ 0, we have ∆ t = 0. (iii) If 0 < A t , M 1 < 1 and 0 < A t , M 2 < 1, we have: ∆ t = A t , M 1 -M 2 A t F = A t F | A t , M 1 -M 2 | ≤ A t 2 F M 1 -M 2 F . (iv) If A t , M 1 ≥ 1 and A t , M 2 ≤ 0, we have ∆ t = A t F ≤ A t F | A t , M 1 -M 2 | ≤ A t 2 F M 1 -M 2 F . (v) If A t , M 1 ≥ 1 and 0 < A t , M 2 < 1, we have: ∆ t = ( A t , M 2 -1)A t F = A t F (1 -A t , M 2 ) ≤ A t F (1 -A t , M 2 ) + A t F ( A t , M 1 -1) = A t F A t , M 1 -M 2 ≤ A t 2 F M 1 -M 2 F . (vi) If A t , M 1 ≤ 0 and 0 < A t , M 2 < 1, we have: ∆ t = -A t -( A t , M 2 -1)A t F = A t A t , M 2 F = A t F A t , M 2 ≤ A t F A t , M 2 -A t F A t , M 1 = A t F A t , M 2 -M 1 ≤ A t 2 F M 1 -M 2 F .
The remaining cases are also bounded by

A t 2 F M 1 -M 2 F by symmetry to cases (iv)-(v)-(vi). Hence ∇f is L-Lipschitz continuous with L = A t 2 F . It is easy to see that diam • F (D λ ) = √ 8λ.
The convergence rate then follows from the general analysis of the FW algorithm [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF].

The second part of the proposition follows directly from the structure of the bases and the greedy nature of the algorithm.

Note that the optimality gap in Proposition 1 is independent of d. Indeed, A t has O(s 2 ) nonzero entries on average, hence the term A t 2 F in the Lipschitz constant L can be bounded by

s 2 A t ∞ , where A ∞ = max d i,j=1 |A i,j |.
This means that Algorithm 3 is able to find a good approximate solution based on a small number of features in only a few iterations, which is very appealing in the high-dimensional setting we consider.

Complexity Analysis

We now analyze the time and memory complexity of Algorithm 3. The form of the gradient (5) along with the structure of the algorithm's updates are crucial to its efficiency. Since M (k+1) is a convex combination of M (k) and a 4-sparse matrix B (k) , we can efficiently compute most of the quantities of interest through careful book-keeping.

In particular, storing M (k) at iteration k requires O(k) memory. We can also recursively compute A t , M (k+1) for all constraints in only O(T ) time and O(T ) memory based on A t , M (k) and A t , B (k) . This allows us, for instance, to efficiently compute the objective value as well as to identify the set of satisfied constraints (those with A t , M (k) ≥ 1) which are ignored in the computation of the gradient. Finding the away direction at iteration k can be done in O(T k) time. For the line search, we use a bisection algorithm to find a root of the gradient of the 1-dimensional function of γ, which only depends on A t , M (k) and A t , B (k) , both of which are readily available. Its time complexity is O(T log 1 ) where is the precision of the line-search, with a memory cost of O(1).

The bottleneck is to find the forward direction. Indeed, sequentially considering each basis element is intractable as it takes O(T d 2 ) time. A more efficient strategy is to sequentially consider each constraint, which requires O(T s 2 ) time and O(T s 2 ) memory. The overall iteration complexity of Algorithm 3 is given in Table 1.

Approximate Forward Step

Finding the forward direction can be expensive when T and s are both large. We propose two strategies to alleviate this cost by finding an approximately optimal basis (see Table 1 for iteration complexity).

Mini-batch approximation. Instead of finding the forward and away directions based on the full gradient at each iteration, we can estimate it on a mini-batch of M T constraints drawn uniformly at random (without replacement). The complexity of finding the forward direction is thus reduced to O(M s 2 ) time Variant Time complexity Memory complexity Exact (Algorithm 3)

Õ(T s 2 + T k) Õ(T s 2 + k) Mini-batch Õ(M s 2 + T k) Õ(T + M s 2 + k) Mini-batch + heuristic Õ(M s + T k) Õ(T + M s + k)
Table 1: Complexity of iteration k (ignoring logarithmic factors) for different variants of the algorithm.

and O(M s 2 ) memory. Consider the deviation between the "value" of any basis element B ∈ B λ on the full set of constraints and its estimation on the mini-batch, namely

1 M t∈M B, G t - 1 T T t=1 B, G t , ( 7
)
where M is the set of M constraint indices drawn uniformly and without replacement from the set {1, . . . , T }. Under mild assumptions, concentration bounds such as Hoeffding's inequality for sampling without replacement [START_REF] Serfling | Probability inequalities for the sum in sampling without replacement[END_REF][START_REF] Bardenet | Concentration inequalities for sampling without replacement[END_REF] can be used to show that the probability of ( 7) being larger than some constant decreases exponentially fast with M . The FW algorithm is known to be robust to inexact gradients, and convergence guarantees similar to Proposition 1 can be obtained directly from [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF][START_REF] Freund | New Analysis and Results for the Conditional Gradient Method[END_REF].

Fast heuristic. To avoid the quadratic dependence on s, we propose to use the following heuristic to find a good forward basis. We first pick a feature i ∈ [d] uniformly at random, and solve the linear problem over the restricted set j {P

(ij) λ , N (ij)
λ }. We then solve it again over the set k {P

(kj) λ , N (kj)
λ } and use the resulting basis for the forward direction. This can be done in only O(M s) time and O(M s) memory and gives good performance in practice, as we shall see in Section 5.

Generalization Analysis

In this section, we derive generalization bounds for the proposed method. Our main goal is to give a theoretical justification of our approach, in particular by (i) showing that our choice of feasible domain D λ helps to reduce overfitting in high dimensions, and (ii) showing that the proposed greedy Frank-Wolfe algorithm provides a simple way to balance between optimization and generalization errors through early stopping.

Setup and Notations

As in previous work on generalization bounds for metric learning, we consider the supervised learning setting where the training sample is a set of labeled points S = {z i = (x i , y i )} n i=1 drawn i.i.d. from a probability distribution µ over the space Z = X × Y, where X ⊆ R d and Y = {1, . . . , C} is the label set. We assume that B X = sup x,x ,x ∈X x(xx ) T is bounded for some convenient matrix norm • .

For simplicity, we assume that the univariate loss function : R → R + is 1-Lipschitz, which is the case for the smoothed hinge loss used in our algorithm. Given a triplet (z, z , z ) ∈ Z 3 , we say that it is admissible if y = y = y . Since we only want to consider admissible triplets, we will use the triplet-wise loss function

L M (z, z , z ) = I[y = y = y ] • ( x(x -x ) T , M ) indexed by M ∈ D λ ,
which is equal to zero for non-admissible triplets.

Given a matrix M ∈ D λ , we define its empirical risk associated on the training set S as follows:

L S (M ) = 1 n(n -1)(n -2) i =j =k L M (z i , z j , z k ). (8)
Similarly, its expected risk is defined as

L(M ) = E z,z ,z ∼µ L M (z, z , z ) . (9)
In contrast to the standard supervised classification setting, note that the empirical risk (8) takes the form of an average of dependent terms known as a U -statistic [START_REF] Lee | U-Statistics: Theory and Practice[END_REF].

From our feasible domain D λ = conv(B λ ), we can define a sequence of nested sets as follows:

D (k) λ = k i=1 α i B i : B i ∈ B λ , α i ≥ 0, k i=1 α i = 1 , k = 1, . . . , 2d(d -1). (10)
In other words, D

λ consists of all d × d matrices which can be decomposed as a convex combination of at most k elements of the basis set B λ . Clearly, we have (d-1)) λ = D λ . Note also that since is 1-Lipschitz, by Holder's inequality we have ∀k:

D (1) λ ⊂ D (2) λ ⊂ • • • ⊂ D (2d
sup z,z ,z ∈Z,M ∈D (k) λ |L M (z, z , z )| ≤ sup x,x ,x ∈X ,M ∈D (k) λ | ( x(x -x ) T , M )| ≤ B X sup M ∈D (k) λ M * , (11)
where • * is the dual norm of • .

In the following, we derive theoretical results that take advantage of the structural properties of our algorithm, namely that the matrix M (k) returned after k ≥ 1 iterations of Algorithm 3 belongs to D λ . We then use these results to derive bounds on the excess risk L(M (k) ) -L(M * ), where M * ∈ arg min M ∈D λ L(M ) is the expected risk minimizer. All proofs can be found in the appendix.

Main Results

We first characterize the Rademacher complexity of the loss functions indexed by elements of

D (k) λ . Given k ∈ {1, . . . , 2d(d -1)}, consider the family F (k) = {L M : M ∈ D (k)
λ } of functions mapping from Z 3 to R + . We will consider the following definition of the Rademacher complexity of F (k) with respect to distribution µ and sample size n ≥ 3, adapted from [START_REF] Clémenc ¸on | Ranking and Empirical Minimization of Ustatistics[END_REF]Cao et al., 2012a):

R n F (k) = E σ,S∼µ n sup M ∈D (k) λ 1 n/3 n/3 i=1 σ i L M (z i , z i+ n/3 , z i+2× n/3 ) , (12)
where σ = (σ 1 , . . . , σ n/3 ) are independent uniform random variables taking values in {-1, 1}. The following lemma gives a bound on the above Rademacher complexity.

Lemma 1 (Bounded Rademacher complexity). Let n ≥ 3, λ > 0 and 1 ≤ k ≤ 2d(d -1). We have

R n (F (k) ) ≤ 8λB X 2 log k n/3 . Proof. See Appendix B.
There are two important consequences to Lemma 1. First, restricting the set of feasible matrices M to D λ = D (2d(d-1)) λ instead of R d×d leads to a Rademacher complexity with a very mild O( √ log d) dependence in the dimension. This validates our design choice for the feasible domain in the high-dimensional setting we consider. Second, the Rademacher complexity can actually be made independent of d by further restricting the number of bases k.

Using this result, we derive a bound for the deviation between the expected risk L(M ) and the empirical risk

L S (M ) of any M ∈ D (k) λ .
Theorem 1 (Maximal deviations). Let S be a set of of n points drawn i.i.d. from µ, λ > 0 and 1 ≤ k ≤ 2d(d -1). For any δ > 0, with probability 1 -δ we have

sup M ∈D (k) λ [L(M ) -L S (M )] ≤ 16λB X 2 log k n/3 + 3B X B D (k) λ 2 ln (2/δ) n , (13)
where

B D (k) λ = sup M ∈D (k) λ M * . Proof. See Appendix C.
The generalization bounds given by Theorem 1 exhibit a standard O(1 √ n) rate. They also confirm that restricting the number k of bases is a good strategy to guard against overfitting when the feature dimension d is high. Interestingly, note that due to the convex hull structure of our basis set, B D (k)

λ = sup M ∈D (k) λ M *
can be easily bounded by a quantity independent of d for any k ≥ 1 and any dual norm • * . We thus have complete freedom to choose the primal norm • so as to make B X = sup x,x ,x ∈X x(xx ) T as small as possible. A good choice of primal norm is the infinity norm A ∞ = max d i,j=1 |A i,j |, which is independent of d. For instance, if X = [0, 1] d we have B X = 1. The dual norm of the infinity norm being the L 1 norm, we then have for any k ≥ 1:

B D (k) λ = sup M ∈D (k) λ M 1 = sup M ∈D (k) λ d i,j=1 |M i,j | ≤ 4λ. ( 14
)
Theorem 1 is directly comparable to the results of Cao et al. (2012a), who derived generalization bounds for similarity learning under various norm regularizers. Their bounds have a similar form, but exhibit a dependence on the feature dimension d which is at least logarithmic (sometimes even linear, depending on the norm used to regularize the empirical risk). In contrast, our bounds depend logarithmically on k d. This offers more flexibility in the high-dimensional setting because k can be directly controlled by stopping our algorithm after k d iterations to guarantee that the output is in

D (k)
λ . This is highlighted by the following corollary, which combines the generalization bounds of Theorem 1 with the O(1/k) convergence rate of our Frank-Wolfe optimization algorithm (Proposition 1).

Corollary 1 (Excess risk bound). Let S be a set of n points drawn i.i.d. from µ, λ > 0. Given k ∈ {1, . . . , 2d(d -1)}, let M (k) be the solution returned after k iterations of Algorithm 3 applied to the problem min M ∈D λ L S (M ), and let M * ∈ arg min M ∈D λ L(M ) be the expected risk minimizer over D λ . For any δ > 0, with probability 1 -δ we have

L(M (k) ) -L(M * ) ≤ 16Lλ 2 k + 2 + 16λB X 2 log k n/3 + 5B X B D (k) λ ln (4/δ) n .
Proof. See Appendix D.

Corollary 1 shows that the excess risk with respect to the expected risk minimizer M * depends on a trade-off between optimization error and complexity of the hypothesis class. Remarkably, this trade-off is ruled by the number k of iterations of the algorithm: as k increases, the optimization error term decreases but the Rademacher complexity terms gets larger. We thus obtain an excess risk bound which adapts to the actual sparsity of the solution output by our algorithm. This is in accordance with our overall goal of reducing overfitting by allowing a strict control on the complexity of the learned similarity, and justifies an early-stopping strategy to achieve a good reduction in empirical risk by selecting the most useful bases while keeping the solution complexity small enough. Again, the excess risk is independent of the feature dimension d, suggesting that in the high-dimensional setting it is possible to find sparse solutions with small excess risk. To the best of our knowledge, this is the first result of this nature for metric or similarity learning.

Remark 2 (Approximation of empirical risk by subsampling). The empirical risk (8) is a sum of O(n 3 ) term, which can be costly to minimize in the large-scale setting. To reduce the computational cost, an alternative to the mini-batch strategy described in Section 3.2.3 is to randomly subsample M terms of the sum (e.g., uniformly without replacement) and to solve the resulting approximate empirical risk minimization problem. For general problems involving U -statistics, Clémenc ¸on et al. (2016) showed that sampling only M = O(n) terms is sufficient to maintain the O(1/ √ n) rate. These arguments can be adapted to our setting to obtain results similar to Theorem 1 and Corollary 1 for this subsampled empirical risk.

Experiments

In this section, we present experiments to evaluate the performance and robustness of HDSL. In Section 5.1, we use synthetic data to study the performance of our approach in terms of similarity recovery and generalization in high dimensions in a controlled environment. Section 5.2 evaluates our algorithm against competing approaches on classification and dimensionality reduction using real-world datasets.

Experiments on Synthetic Data

We first conduct experiments on synthetic datasets in order to address two questions:

1. Is the algorithm able to recover the ground truth sparse similarity function from (potentially weak) similarity judgments?

2. How well does the algorithm generalize as the dimensionality increases?

Similarity Recovery

To investigate the algorithm's ability to recover the underlying similarity, we generate a ground truth similarity metric M ∈ R d×d where d = 2000. M is constructed as a convex combination of 100 randomly selected rank-one 4-sparse bases as specified in Section 3.1. The combination coefficients are drawn from a Dirichlet distribution with shape parameter 9 and scale parameter 0.5. Without loss of generality, we choose the metric to be block structured by restricting the basis selection from two blocks. This makes the resulting matrix easier to visualize, as show in Figure 2(a).

We then generate 5000 training samples from the uniform distribution on [0, 1] with 2% sparsity. From this sample, we create 30, 000 training triplets {(x 1 , x 2 , x 3 )} where x 1 is randomly picked and x 2 (or x 3 ) is sampled among x 1 's top α% similar (or dissimilar) samples as measured by the ground truth metric M . The parameter α controls the "quality" of the triplet constraints: a larger α leads to less similar (or dissimilar) samples in the triplets, thereby providing a weaker signal about the underlying similarity. We experiment with various α (10%, 20%, 25%, 30%) to investigate the robustness of HDSL to the quality of the supervision. In all our experiments, we use λ = 100.

Results. We aim to measure how accurately we recover the entries (i.e., pairs of features) that are active in the ground truth similarity as training proceeds. To do so, at each iteration k of HDSL, we rank each pair of features by descending order of the absolute value of the corresponding entry in the current matrix M (k) . We then compute the Area under the ROC Curve (AUC) of the ranking induced by the similarity with respect to the list of active entries in the ground truth similarity. The AUC is well-suited to the imbalanced setting (as active entries in the ground truth are a small subset of all entries). It can be interpreted as the probability that a random entry that is active in the ground truth is ranked higher than a random inactive one. Following a similar process, we also compute an AUC score for individual features: this is done by ranking each feature by the L 1 norm of its associated row in the matrix.

The AUC scores for feature and entry recovery along the iterations are reported in Figure 1 for different values of α. When the quality of the triplet constraints is high (α =10%,20%), the AUC increases quickly to converge very close to 1.0, indicating an almost perfect recovery of relevant features/entries. This confirms that HDSL is able to accurately identify the small number of correct features and pairs of features. As α increases (i.e., the similarity constraints become noisy and less informative), the AUC increases at a slower pace and the final value decreases. This is expected as the quality of the information carried by the similarity judgments is key to recover the ground truth similarity. Yet, even for α =30%, the final AUC score is still very high (above 0.85 for both feature and entry recovery). This good recovery behavior is confirmed by the visual representations of the ground truth and learned similarity matrices shown in Figure 2. We observe that the learned similarity (when α = 20%) clearly recovers the block structure of the true similarity, and is able to correctly identify most individual entries with very few false positives.

Link Prediction

We now investigate the ability of our algorithm to generalize well as the feature dimensionality increases by conducting a signed link prediction experiment, which is the task of distinguishing positive and negative interactions in a network (see e.g. [START_REF] Agrawal | Link label prediction in signed social networks[END_REF].

We generate 500 samples with different number of features d ranging from 5, 000 to 1, 000, 000. As the dimension d increases, we decrease the average sparsity of data (from 0.02 to 0.002) to limit running time. In real high-dimensional datasets, features typically do not appear in a uniform frequency: instead, a small portion of features tends to dominate the others. Following this observation, we generate features whose frequency follow a power law style distribution, as shown in Figure 3(a). The ground truth similarity is then a convex combination of randomly selected bases as in the previous experiment, except that we restrict the selected bases to those involving features that are frequent enough (a frequency of at least 0.1 was chosen for this experiment). This is needed to ensure that the features involved in the ground truth similarity will occur at least a few times in our small dataset, but we emphasize that the algorithm is exposed to the entire feature set and does not know which features are relevant.

Based on the samples and the ground truth similarity, we generate signed link observations of the form

{x i 1 , x i 2 , y i } N i (y i ∈ {-1, 1}
). We associate the label y i = 1 (positive link) to pairs for which the similarity between x 1 and x 2 ranks in the top 5% of x 1 's (or x 2 's) neighbors according to the ground truth similarity measure. On the other hand, y i = -1 (negative link) indicates that the similarity ranks in the bottom 5% of x 1 's (or x 2 's) neighbors. We split these link observations into training, validation and test sets of 1, 000 observations each. Triplets constraints are generated from training links -given a pair x 1 , x 2 , y, we randomly sample x 3 as a similar (if y = -1) or dissimilar (if y = 1) node. The validation set is used to tune the hyperparameter λ and for early stopping.

Results. We measure the generalization ability of HDSL by the AUC score of link prediction on the test set. Figure 3(b) reports these AUC scores across different dimensions. We also show results for different numbers of constraints per training link. The results are averaged over 5 random runs. As one would expect, the task becomes increasingly difficult as the dimension becomes larger, since the size of the training set is fixed (1, 000 training links generated from 500 nodes). However, the performance decreases slowly (roughly logarithmically) with the dimension, and we achieve very high AUC scores (larger than 0.9) even for one million features. We also see that training from more constraints tends to improve the prediction performance.

Experiments on Real Datasets

We now present comparative experiments on several high-dimensional real datasets, evaluating our approach against several baselines and competing methods. Competing methods. We compare the proposed approach (HDSL) to several methods:

• DOT: The standard dot product, which is equivalent to setting M = I.

• DIAG: Diagonal similarity learning (i.e., a weighting of the features), as done in [START_REF] Gao | SOML: Sparse Online Metric Learning with Application to Image Retrieval[END_REF]. We obtain it by minimizing the same loss as in HDSL with 2 and 1 regularization, i.e., min

w∈R d f (w) = 1 T T t=1 A t , diag(w) + λΩ(w),
where Ω(w) ∈ { w2 2 , w 1 } and λ is the regularization parameter. Optimization was done using (proximal) gradient descent.

• RP+OASIS: Similarity learning in random projected space. Given r d, let R ∈ R d×r be a matrix where each entry r ij is randomly drawn from N (0, 1). For each data instance x ∈ R d , we generate x = 1 √ r R T x ∈ R r and use this reduced data in OASIS [START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF], a fast online method to learn a bilinear similarity from triplet constraints.

• PCA+OASIS: Similarity learning in PCA space. Same as RP+OASIS, except that PCA is used instead of random projections to project the data into R r .

• SVM: Support Vector Machines. We use linear SVM, which is known to perform well for sparse highdimensional data [START_REF] Caruana | An empirical evaluation of supervised learning in high dimensions[END_REF], with 2 and 1 regularization. We also use nonlinear SVM with the polynomial kernel (2nd and 3rd degree) popular in text classification [START_REF] Chang | Training and Testing Low-degree Polynomial Data Mappings via Linear SVM[END_REF].

The SVM models are trained using liblinear [START_REF] Fan | LIBLINEAR: A Library for Large Linear Classification[END_REF] and libsvm [START_REF] Chang | LIBSVM : a library for support vector machines[END_REF] with 1vs1 paradigm for multiclass. We have also tried to compare our method with COMET [START_REF] Atzmon | Learning sparse metrics, one feature at a time[END_REF], which also learns a bilinear similarity in a greedy fashion with rank-1 updates. However, as mentioned in Section 2.1 their coordinate descent algorithm has a time complexity of O(d 2 ) per iteration, as well as overall memory complexity of O(d 2 ). We run the sparse version of code provided by the authors3 on a machine with a 2.3GHz Intel Core i7 and 16GB memory. On the dexter dataset (which has the smallest dimensionality in our benchmark), a single pass over the features took more than 4 hours, while the authors reported that about 10 passes are generally needed for COMET to converge [START_REF] Atzmon | Learning sparse metrics, one feature at a time[END_REF]. On the dorothea dataset, COMET returned a memory error. As a result, we did not include COMET to our empirical comparison. In contrast, on the same hardware, our approach HDSL takes less than 1 minute on dexter and less than 1 hour on dorothea.

Training Procedure. For all similarity learning algorithms, we generate 15 training constraints for each instance by identifying its 3 target neighbors (nearest neighbors with same label) and 5 impostors (nearest neighbors with different label), following [START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF]. Due to the very small number of training instances in dexter, we found that better performance is achieved across all methods when using 20 training constraints per instance, drawn at random based on its label. All parameters are tuned using the accuracy on the validation set. For HDSL, we use the fast heuristic described in Section 3.2.3 and tune the scale parameter λ ∈ {1, 10, . . . , 10 9 }. The regularization parameters of DIAG and SVM are tuned in {10 -9 , . . . , 10 8 } and the "aggressiveness" parameter of OASIS is tuned in {10 -9 , . . . , 10 2 }.

Results

Classification Performance. We first investigate the performance of each similarity learning approach in k-NN classification (k was set to 3 for all experiments). For RP+OASIS and PCA+OASIS, we choose the dimension r of the reduced space based on the accuracy of the learned similarity on the validation set, limiting our search to r ≤ 2000 because OASIS is extremely slow beyond this point. 4 Similarly, we use the performance on validation data to do early stopping in HDSL, which also has the effect of restricting the number of features used by the learned similarity.

Table 3 shows the k-NN classification performance. We can first observe that RP+OASIS often performs worse than DOT, which is consistent with previous observations showing that a large number of random projections may be needed to obtain good performance [START_REF] Fradkin | Experiments with random projections for machine learning[END_REF]. PCA+OASIS gives much better results, but is generally outperformed by a simple diagonal similarity learned directly in the original high-dimensional space. HDSL, however, outperforms all other algorithms on these datasets, including DIAG. This shows the good generalization performance of the proposed approach, even though the number of training samples is sometimes very small compared to the number of features, as in dexter and dorothea. It also shows the relevance of encoding "second order" information (pairwise interactions between the original features) in the similarity instead of considering a simple weighting of features as in DIAG.

Table 4 shows the comparison with SVMs. Interestingly, HDSL outperforms all SVM variants on dexter and dorothea, both of which have a large proportion of irrelevant features. This shows that its greedy strategy and early stopping mechanism achieves better feature selection and generalization than the 1 version of linear SVM. On the other two datasets, HDSL is competitive with SVM, although it is outperformed slightly by one variant (SVM-poly-3 on rcv1 2 and SVM-linear-2 on rcv1 4), both of which rely on all features.

Feature selection and sparsity. We now focus on the ability of HDSL to perform feature selection and more generally to learn sparse similarity functions. To better understand the behavior of HDSL, we show in Figure 4 the number of selected features as a function of the iteration number for two of the datasets. Remember that at most two new features can be added at each iteration. Figure 4 shows that HDSL incorporates many features early on but tends to eventually converge to a modest fraction of features (the same observation holds for the other two datasets). This may explain why HDSL does not suffer much from overfitting even when training data is scarce as in dexter.

Another attractive characteristic of HDSL is its ability to learn a matrix that is sparse not only on the diagonal but also off-diagonal (the proportion of nonzero entries is in the order of 0.0001% for all datasets). In other words, the learned similarity only relies on a few relevant pairwise interactions between features. Figure 5 shows two examples, where we can see that HDSL is able to exploit the product of two features as either a positive or negative contribution to the similarity score. This opens the door to an analysis of the importance of pairs of features (for instance, word co-occurrence) for the application at hand. Finally, the extreme sparsity of the matrices allows very fast similarity computation. Together with the superior accuracy brought by HDSL, it makes our approach potentially useful in a variety of contexts (k-NN, clustering, ranking, etc).

Finally, it is also worth noticing that HDSL uses significantly less features than DIAG-1 (see numbers in brackets in Table 3). We attribute this to the extra modeling capability brought by the non-diagonal similarity observed in Figure 5. 5Dimension reduction. We now investigate the potential of HDSL for dimensionality reduction. Recall that HDSL learns a sequence of PSD matrices M (k) . We can use the square root of M (k) to project the data into a new space where the dot product is equivalent to S M (k) in the original space. The dimension of the projection space is equal to the rank of M (k) , which is upper bounded by k + 1 (see Section 3.1). A single run of HDSL can thus be seen as incrementally building projection spaces of increasing dimensionality.

To assess the dimensionality reduction quality of HDSL (measured by k-NN classification error on the test set), we plot its performance at various iterations during the runs that generated the results of Table 3. We compare it to two standard dimensionality reduction techniques: random projection and PCA. We also evaluate RP+OASIS and PCA+OASIS, i.e., learn a similarity with OASIS on top of the RP and PCA features.6 Note that OASIS was tuned separately for each projection size, making the comparison a bit unfair to HDSL. The results are shown in Figure 6. As observed earlier, random projection-based approaches achieve poor performance. When the features are not too noisy (as in rcv1 2 and rcv1 4), PCA-based methods are better than HDSL at compressing the space into very few dimensions, but HDSL eventually catches up. On the other hand, PCA suffers heavily from the presence of noise (dexter and dorothea), while HDSL is able to quickly improve upon the standard similarity in the original space. Finally, on all datasets, we observe that HDSL converges to a stationary dimension without overfitting, unlike PCA+OASIS which exhibits signs of overfitting on dexter and rcv1 4 especially.

Concluding Remarks

In this work, we proposed an efficient approach to learn similarity functions from high-dimensional sparse data. This is achieved by forming the similarity as a combination of simple sparse basis elements that operate on only two features and the use of an (approximate) Frank-Wolfe algorithm. Our algorithm is completed by a novel generalization analysis which validates the design choices and highlights the robustness of our approach to high dimensions. Experiments on synthetic and real datasets confirmed the good practical behavior of our method for classification and dimensionality reduction. The learned similarity may be applied to other algorithms that rely on a similarity function (clustering, ranking), or as a way to preprocess the data before applying another learning algorithm. We also note that [START_REF] St | Sparse Compositional Local Metric Learning[END_REF] have recently extended our HDSL algorithm to learn local metrics for different regions of the space in addition to the global metric.

We leave several fundamental questions for future work. In particular, our framework could be extended to optimize a loss function related to a linear classification objective. We could then attempt to adapt our analysis to obtain generalization bounds directly for the classification error. Such bounds exist in the literature (see [START_REF] Bellet | Similarity Learning for Provably Accurate Sparse Linear Classification[END_REF][START_REF] Guo | Guaranteed Classification via Regularized Similarity Learning[END_REF] but exhibit a classic dependence on the data dimension that could be avoided with our approach. Another interesting, though challenging direction is to formally study the conditions under which a sparse ground truth similarity can be accurately recovered from similarity judgments. Inspiration could be drawn from the related problem of sparse recovery in the compressed sensing literature [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF].

The next technical lemma is based on the above representation. Lemma 4 (McDiarmid, 1989). Let Z be some set and let f : Z n → R be a function of n variables such that for some c > 0, for all i ∈ {1, . . . , n} and for all z 1 , . . . , z n , z i ∈ Z, we have |f (z 1 , . . . , z i-1 , z i , z i+1 , . . . , , z n ) -f (z 1 , . . . , z i-1 , z i , z i+1 , . . . , z n )| ≤ c.

Let Z 1 , . . . , Z n be n independent random variables taking values in Z. Then, with probability at least 1 -δ, we have 15) follows from the contraction property (see Shalev-Shwartz and Ben-David, 2014, Lemma 26.9). The inequality (16) follows from the fact M is a convex combination of set of k bases combined with the properties in Shalev-Shwartz and Ben-David (2014, Lemma 26.7, 26.8). Finally, inequality (17) follows from the sparsity structure of the bases and the fact that x i (x j -x k ) T has no entries with absolute value greater than B X .

Let S = {z 1 , . . . , z q-1 , z q , z q+1 , . . . , z n } and S = {z 1 , . . . , z q-1 , z q , z q+1 , . . . , z n } be two samples differing by exactly one point. We have:

Φ(S ) -Φ(S) ≤ sup M ∈D (k) λ [L S (M ) -L S (M )] ≤ 1 n(n -1)(n -2) sup M ∈D (k) λ i =j =k |L M (z i , z j , z k ) -L M (z i , z j , z k )| ≤ 1 n(n -1)(n -2) 6(n -1)(n -2)B X B D (k) λ = 6 n B X B D (k) λ .
The first inequality comes from the fact that the difference of suprema does not exceed the supremum of the difference. The last inequality makes use of ( 11). Similarly, we can obtain Φ(S) -Φ(S ) ≤ 6B X B D (k) ).

  the Rademacher complexity of D (k) λ and derive bounds on the maximal deviation between L(M ) and L S (M ) for any M ∈ D (k)

Figure 1 :

 1 Figure 1: Similarity recovery experiment on synthetic data. Figure 1(a) and Figure 1(b) show the AUC scores (for feature recovery and entry recovery respectively) along the iterations of the algorithm for different values of α.

  Figure 2: Similarity recovery experiment on synthetic data. Figure 2(a) shows the underlying ground truth similarity, where blue dots represent positive entries and red dots represent negative entries (combination coefficients are not displayed). Figure 2(b) shows the similarity learned by HDSL (α =20%), which is visually very close to the ground truth.

  Figure3: Link prediction experiment on synthetic data. Figure3(a)shows the feature frequency distribution, which follows a power law as in many real high-dimensional datasets. Figure3(b)shows AUC scores on the test set for different number of features (in log scale) and number of training constraints per link.

Figure 4 :

 4 Figure 4: Number of active features learned by HDSL as a function of the iteration number.

  Figure 5: Sparsity structure of the matrix M learned by HDSL. Positive and negative entries are shown in blue and red, respectively (best seen in color).

  Figure 6: k-NN test error as a function of the dimensionality of the space (in log scale). Best seen in color.

Lemma 3 .

 3 Let Q be a set of functions from Z 3 to R. If z 1 , z 2 , ..., z n ∈ Z are i.i.d., then we haveE[sup q∈Q 1 n(n -1)(n -2) i =j =k q(z i , z j , z k )] ≤ E[sup i , z i+ n/3 , z i+2× n/3 )].Proof. From Lemma 2, we observe that π(i) , z π(i+ n/3 ) , z π(i+2× n/3 ) )] π(i) , z π(i+ n/3 ) , z π(i+2× n/3 ) π(i) , z π(i+ n/3 ) , z π(i+2× n/3 ) i , z i+ n/3 , z i+2× n/3 )],which proves the result.Finally, we recall McDiarmid's inequality.

  |f (Z 1 , . . . , Z n ) -E[f (Z 1 , . . . , Z n )]| ≤ c n log(2/δ) 2 .Appendix B Proof of Lemma 1Proof. Given a training sample S = {z i = (x i , y i ) : i ∈ 1, . . . , n} ∼ µ n , we denote the set of admissible triplets involved in the Rademacher complexity byA S = i : y i = y i+ n/3 = y i+2× n/3 , i = 1, . . . , n/3 ,and let m = |A S | ≤ n/3 . We have:R n (F (k) ) = E σ,S∼µ n sup i ( x i (x i+ n/3 -x i+2× n/3 ) T , M ) ≤ E σ,S∼µ n sup i (x i+ n/3 -x i+2× n/3 ) T , M σ i x i (x i+ n/3 -x i+2× n/3 ) T , U = {u τ ∈ R m : τ = 1, . . . , k, (u τ ) i = x γ(i) (x γ(i)+ n/3 -x γ(i)+2× n/3 ) T , B τ , γ : {1, . . . , m} → A S isbijective, B τ ∈ B λ }, and ū = 1 k k τ =1 u τ . The inequality (

σ

  |Φ(S) -Φ(S )| ≤ 6B X B D (k) λ /n. We can therefore apply McDiarmid's inequality (see Lemma 4 in Appendix A) to Φ(S): for any δ > 0, with probability at least 1 -δ we have:sup M ∈D (k) λ [L(M ) -L S (M )] ≤ E S sup M ∈D (k) λ [L(M ) -L S (M )] + 3B X B D (k) to bound E S sup M ∈D (k) λ [L(M ) -L S (M )]. Applying Lemma 3 (see Appendix A) with q M (z, z , z ) = L(M ) -L M (z, z , z ) gives (z i , z i+ n/3 , z i+2× n/3). Let S = {z 1 , ..., zn } be an i.i.d. sample independent of S. ThenE S sup M ∈D (k) λ [L(M ) -LS (M )] = E S sup M ∈D (k) λ [E S L S (M ) -LS (M )] ≤ E S, S sup M ∈D (k) λ [ L S (M ) -LS (M )].Let σ 1 , . . . , σ n 3 ∈ {-1, 1} be a collection of i.i.d. Rademacher variables. By standard symmetrization techniques, we have that E S, S sup M ( zi , zi+ n/3 , zi+2× n/3 ) -L M (z i , z i+ n/3 , z i+2× n/3 ) M ( zi , zi+ n/3 , zi+2× n/3 ) + E σ,S sup M (z i , z i+ n/3 , z i+2× n/3 i L M (z i , z i+ n/3 , z i+2× n/3 ) = 2R n (F

Table 2 :

 2 Datasets used in the experiments5.2.1 SetupDatasets. We report experimental results on several real-world classification datasets with up to 100,000 features. Dorothea and dexter come from the NIPS 2003 feature selection challenge[START_REF] Guyon | Result Analysis of the NIPS 2003 Feature Selection Challenge[END_REF] and are respectively pharmaceutical and text data with predefined splitting into training, validation and test sets. They both contain a large proportion of noisy/irrelevant features. Reuters CV1 is a popular text classification dataset with bag-of-words representation. We use the binary classification version from the LIBSVM dataset collection 2 (with 60%/20%/20% random splits) and the 4-classes version (with 40%/30%/30% random splits) introduced by[START_REF] Cai | Manifold Adaptive Experimental Design for Text Categorization[END_REF]. Detailed information on the datasets and splits is given in Table2. All datasets are normalized such that each feature takes values in [0, 1].

Table 3 :

 3 k-NN test error (%) of the similarities learned with each method. The number of features used by each similarity (when smaller than d) is given in brackets. Best accuracy on each dataset is shown in bold.

	Dataset	SVM-poly-2	SVM-poly-3	SVM-linear-2	SVM-linear-1	HDSL
	dexter	9.4	9.2	8.9	8.9 [281]	6.5 [183]
	dorothea	7	6.6	8.1	6.6 [366]	6.5 [731]
	rcv1 2	3.4	3.3	3.5	4.0 [1915]	3.4 [2126]
	rcv1 4	5.7	5.7	5.1	5.7 [2770]	5.7 [1888]

Table 4 :

 4 Test error (%) of several SVM variants compared to HDSL. As in Table3, the number of features is given in brackets and best accuracies are shown in bold.

This is in contrast to a different line of work, inspired by the problem of ordinal embedding, which aims to learn a metric which correctly orders a fixed set of known points (see for instanceJain et al.,

2017)

http://www.csie.ntu.edu.tw/ ˜cjlin/libsvmtools/datasets/

https://github.com/yuvalatzmon/COMET

Note that the number of PCA dimensions is at most the number of training examples. Therefore, for dexter and dorothea, r is at most 300 and 800 respectively.

Note that HDSL uses roughly the same number of features as SVM-linear-1 (Table4), but it is difficult to draw any solid conclusion because the objective and training data for each method are different, and SVM is a combination of binary models.

Again, we were not able to run OASIS beyond a certain dimension due to computational complexity.

Acknowledgments. This work was partially supported by a grant from CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020. It was also supported in part by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Defense U.S. Army Research Laboratory (DoD / ARL) contract number W911NF-12-C-0012, a NSF IIS-1065243, an Alfred. P. Sloan Research Fellowship, DARPA award D11AP00278, and an ARO YIP Award (W911NF-12-1-0241). The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoD/ARL, or the U.S. Government.

Appendix A Technical Lemmas

The following classic result, known as the first Hoeffding's decomposition, allows to represent a U -statistic as a sum of i.i.d. blocks.

Lemma 2 [START_REF] Hoeffding | A Class of Statistics with Asymptotically Normal Distribution[END_REF]. Let q : Z × Z × Z → R be a real-valued function. Given the i.i.d. random variables z 1 , z 2 , ..., z n ∈ Z, we have

Proof. Observe that ∀i = j = k, q(z i , z j , z k ) appears once on the left hand side and U n (q) has

of its value, while on the right hand side it appears (n -3)! × n/3 times, because for each of the n/3 positions there are (n -3)! possible permutations. Thus the right hand side also has 1 n(n-1)(n-2) of its function value. We thus have the equality.

Appendix C Proof of Theorem 1

Proof. Let us consider the function

We have thus shown:

Plugging ( 19) into (18) and using Lemma 1, we get the desired result.

Appendix D Proof of Corollary 1

Proof. The excess risk of M (k) with respect to M * can be decomposed as follows:

generalization error

where M S ∈ arg min M ∈D λ L S (M ) is an empirical risk minimizer.

The generalization error term in (20) can be bounded using Theorem 1 (recalling that M (k) ∈ D

(k) λ by construction), while the optimization error term is bounded by the convergence rate of our Frank-Wolfe algorithm (Proposition 1). In the last term, M * does not depend on S, hence we can use Hoeffding's inequality together with ( 11) and ( 14) to obtain that for any δ > 0, with probability at least 1 -δ/2:

We get the corollary by combining the above results using the union bound.


PFT3D

(Plateforme Technologique 3D)

Archéotransfert

•

  Prev. : Integrated High-Quality Acquisition of Geometry and Appearance for Cultural Heritage -Schwartz, Weinmann,Ruiters, and Klein -VAST 2011

  and Transmission of Knowledge 3. Life cycle, interoperability & archiving 4. HSS meetings: exchanges of uses, communication & formation consortium3d-coordination@services.cnrs.fr