Year: 2019

Abstract:
Learning Mahalanobis metric spaces is an important problem that has found numerous applications. Several algorithms have been designed for this problem, including Information Theoretic Metric Learning (ITML) [Davis et al. 2007] and Large Margin Nearest Neighbor (LMNN) classification [Weinberger and Saul 2009]. We consider a formulation of Mahalanobis metric learning as an optimization problem, where the objective is to minimize the number of violated similarity/dissimilarity constraints. We show that for any fixed ambient dimension, there exists a fully polynomial-time approximation scheme (FPTAS) with nearly-linear running time. This result is obtained using tools from the theory of linear programming in low dimensions. We also discuss improvements of the algorithm in practice, and present experimental results on synthetic and real-world data sets. Our algorithm is fully parallelizable and performs favorably in the presence of adversarial noise.

heading: 1 INTRODUCTION
text: Learning metric spaces is a fundamental computational primitive that has found numerous applications and has received significant attention in the literature. We refer the reader to Kulis et al. (2013); Li and Tian (2018) for detailed exposition and discussion of previous work. At the high level, the input to a metric learning problem consists of some universe of objects X , together with some similarity information on subsets of these objects. Here, we focus on pairwise similarity and dissimilarity constraints. Specifically, we are given S,D Ă `
X 2
˘
, which are sets of pairs of objects that are labeled as similar and dissimilar respectively. We are also given some u, ` ą 0, and we seek to find a mapping f : X Ñ Y , into some target metric space pY, ρq, such that for all x, y P S,
ρpfpxq, fpyqq ď u,
and for all x, y P D, ρpfpxq, fpyqq ě `.
In the case of Mahalanobis metric learning, we have X Ă Rd, with |X| “ n, for some d P N, and the mapping f : Rd Ñ Rd is linear. Specifically, we seek to find a matrix G P Rdˆd, such that for all tp, qu P S, we have
}Gp´Gq}2 ď u, (1)
and for all tp, qu P D, we have
}Gp´Gq}2 ě `. (2)

heading: 1.1 OUR CONTRIBUTION
text: In general, there might not exist any G that satisfies all constraints of type 1 and 2. We are thus interested in finding a solution that minimizes the fraction of violated constraints, which corresponds to maximizing the accuracy of the mapping. We develop a p1 ` εq-approximation algorithm for optimization problem of computing a Mahalanobis metric space of maximum accuracy, that runs in near-linear time for any fixed ambient dimension d P N. This algorithm is obtained using tools from geometric approximation algorithms and the theory of linear programming in small dimension. The following summarizes our result.
Theorem 1.1. For any d P N, ε ą 0, there exists a randomized algorithm for learning d-dimensional Mahalanobis metric spaces, which given an instance that admits a mapping with accuracy r˚, computes a mapping with accuracy at least r˚´ε, in time dOp1qnplog n{εqOpdq, with high probability.
The above algorithm can be extended to handle various forms of regularization. We also propose several modifications of our algorithm that lead to significant performance improvements in practice. The final algorithm is evaluated experimentally on both synthetic and real-world data sets, and is compared against the currently best-known algorithms for the problem.

heading: 1.2 RELATED WORK
text: Several algorithms for learning Mahalanobis metric spaces have been proposed. Notable examples include the SDP based algorithm of Xing et al. Xing et al. (2003), the algorithm of Globerson and Roweis for the fully supervised setting Globerson and Roweis (2006), Information Theoretic Metric Learning (ITML) by Davis et al. Davis et al. (2007), which casts the problem as a particular optimization minimizing LogDet divergence, as well as Large Margin Nearest Neighbor (LMNN) by Weinberger et al. Weinberger et al. (2006), which attempts to learn a metric geared towards optimizing k-NN classification. We refer the reader to the surveys Kulis et al. (2013); Li and Tian (2018) for a detailed discussion of previous work. Our algorithm differs from previous approaches in that it seeks to directly minimize the number of violated pairwise distance constraints, which is a highly non-convex objective, without resorting to a convex relaxation of the corresponding optimization problem.

heading: 1.3 ORGANIZATION
text: The rest of the paper is organized as follows. Section 2 describes the main algorithm and the proof of Theorem 1.1. Section 3 discusses practical improvements used in the implementation of the algorithm. Section 4 presents the experimental evaluation.

heading: 2 MAHALANOBIS METRIC LEARNING AS AN LP-TYPE PROBLEM
text: In this Section we present an approximation scheme for Mahalanobis metric learning in d-dimensional Euclidean space, with nearly-linear running time. We begin by recalling some prior results on the class of LP-type problems, which generalizes linear programming. We then show that linear metric learning can be cast as an LP-type problem.

heading: 2.1 LP-TYPE PROBLEMS
text: Let us recall the definition of an LP-type problem. LetH be a set of constraints, and let w : 2H Ñ RY t´8,`8u, such that for any G Ă H, wpGq is the value of the optimal solution of the instance defined by G. We say that pH, wq defines an LP-type problem if the following axioms hold:
(A1) Monotonicity. For any F Ď G Ď H, we have wpF q ď wpGq. (A2) Locality. For any F Ď G Ď H, with ´8 ă wpF q “ wpGq, and any h P H, if wpGq ă
wpGY thuq, then wpF q ă wpF Y thuq.
More generally, we say that pH, wq defines an LP-type problem on someH1 Ď H, when conditions (A1) and (A2) hold for all F Ď G Ď H1. A subset B Ď H is called a basis if wpBq ą ´8 and wpB1q ă wpBq for any proper subset B1 Ĺ B. A basic operation is defined to be one of the following:
(B0) Initial basis computation. Given some G Ď H, compute any basis for G. (B1) Violation test. For some h P H and some basis B Ď H, test whether wpB Y thuq ą wpBq
(in other words, whether B violates h).
(B2) Basis computation. For some h P H and some basis B Ď H, compute a basis of B Y thu.

heading: 2.2 AN LP-TYPE FORMULATION
text: We now show that learning Mahalanobis metric spaces can be expressed as an LP-type problem. We first note that we can rewrite (1) and (2) as
pp´ qqTApp´ qq ď u2, (3)
and
pp´ qqTApp´ qq ě `2, (4)
where A “ GTG is positive semidefinite.
We defineH “ t0, 1uˆ `Rd
2
˘
, where for each p0, tp, quq P H, we have a constraint of type (3), and for every p1, tp, quq P H, we have a constraint of type (4). Therefore, for any set of constraints F Ď H, we may associate the set of feasible solutions for F with the set AF of all positive semidefinite matrices A P Rnˆn, satisfying (3) and (4) for all constraints in F . Let w : 2H Ñ R, such that for all F P H, we have
wpF q “ " infAPAF r TAr if AF ‰ H
8 if AF “ H ,
where r P Rd is a vector chosen uniformly at random from the unit sphere from some rotationallyinvariant probability measure. Such a vector can be chosen, for example, by first choosing some r1 P Rd, where each coordinate is sampled from the normal distribution N p0, 1q, and setting r “ r1{}r1}2. Lemma 2.1. When w is chosen as above, the pair pH, wq defines an LP-type problem of combinatorial dimension Opd2q, with probability 1. Moreover, for any n ą 0, if each ri is chosen using Ωplog nq bits of precision, then for each F Ď H, with n “ |F |, the assertion holds with high probability.
Proof. Since adding constraints to a feasible instance can only make it infeasible, it follows that w satisfies the monotonicity axiom (A1).
We next argue that the locality axion (A2) also holds, with high probability. Let F Ď G Ď H, with ´8 ă wpF q “ wpGq, and let h P H, with wpGq ă wpGY thuq. Let AF P AF and AG P AG be some (not necessarily unique) infimizers of wpAq, when A ranges in AF and AG respectively. The setAF , viewed as a convex subset of Rd 2
, is the intersection of the SDP cone with n half-spaces, and thus AF has at most n facets. There are at least two distinct infimizers for wpAGq, when AG P AG, only when the randomly chosen vector r is orthogonal to a certain direction, which occurs with probability 0. When each entry of r is chosen with c log n bits of precision, the probability that r is orthogonal to any single hyperplane is at most 2´c logn “ n´c; the assertion follows by a union bound over n facets. This establishes that axiom (A2) holds with high probability.
It remains to bound the combinatorial dimension, κ. Let F Ď H be a set of constraints. For each A P AF , define the ellipsoid EA “ tv P Rd : }Av}2 “ 1u. For any A,A1 P AF , with EA “ EA1 , and A “ GTG, A1 “ G1TG1, we have that for all p, q P Rd, }Gp´Gq}2 “ pp´ qqTApp´ qq “ pp´ qqTA1pp´ qq “ }G1p´G1q}2. Therefore in order to specify a linear transformation G, up to an isometry, it suffices to specify the ellipsoid EA. Each tp, qu P S corresponds to the constraint that the point pp´ qq{u must lie in EA. Similarly each tp, qu P D corresponds to the constraint that the point pp´ qq{` must lie either on the boundary or the exterior of EA. Any ellipsoid in Rd is uniquely determined by specifying at most pd`3qd{2 “ Opd2q distinct points on its boundary (see Welzl (1991); Chazelle (2000)). Therefore, each optimal solution can be uniquely specified as the intersection of at most Opd2q constraints, and thus the combinatorial dimension is Opd2q.
Lemma 2.2. Any initial basis computation (B0), any violation test (B1), and any basis computation (B2) can be performed in time dOp1q.
Proof. The violation test (B1) can be performed by solving one SDP to compute wpBq, and another to compute wpB Y thuq. By Lemma 2.1 the combinatorial dimension is Opd2q, thus each SDP has Opd2q constraints, and be solved in time dOp1q. The basis computation step (B2) can be performed starting with the set of constraints B Y thu, and iteratively remove every constraint whose removal does not decrease the optimum cost, until we arrive at a minimal set, which is a basis. In total, we need to solve at most d SDPs, each of size Opd2q, which can be done in total time dOp1q. Finally, by the choice of w, any set containing a single constraint in S is a valid initial basis.

heading: 2.3 ALGORITHMIC IMPLICATIONS
text: Using the above formulation of Mahalanobis metric learning as an LP-type problem, we can obtain our approximation scheme. Our algorithm uses as a subroutine an exact algorithm for the problem (that is, for the special case where we seek to find a mapping that satisfies all constraints). We first present the exact algorithm and then show how it can be used to derive the approximation scheme.
An exact algorithm. Welzl (1991) obtained a simple randomized linear-time algorithm for the minimum enclosing ball and minimum enclosing ellipsoid problems. This algorithm naturally extends to general LP-type problems (we refer the reader to Har-Peled (2011); Chazelle (2000) for further details).
With the interpretation of Mahalanobis metric learning as an LP-type problem given above, we thus obtain a linear time algorithm for in Rd, for any constant d P N. The resulting algorithm on a set of constraints F Ď H is implemented by the procedure Exact-LPTMLpF ;Hq, which is presented in Algorithm 1. The procedure LPTMLpF ;Bq takes as input sets of constraints F,B Ď H. It outputs a solution A P Rdˆd to the problem induced by the set of constraints F YB, such that all constraints in B are tight (that is, they hold with equality); if no such solution solution exists, then it returns nil. The procedure Basic-LPTMLpBq computes LPTMLpH;Bq. The analysis of Welzl (1991) implies that when Basic-LPTMLpBq is called, the cardinality of B is at most the combinatorial dimension, which by Lemma 2.1 is Opd2q. Thus the procedure Basic-LPTML can be implemented using one initial basis computation (B0) and Opd2q basis computations (B2), which by Lemma 2.2 takes total time dOp1q.
Algorithm 1 An exact algorithm for Mahalanobis metric learning. procedure Exact-LPTML(F ;B)
if F “ H then AÐ Basic-LPTMLpBq else choose h P F uniformly at random AÐ Exact-LPTMLpF ´ thuq if A violates h then
A :“ Exact-LPTMLpF ´ thu;B Y thuq end if
end if return A
end procedure
An p1` εq-approximation algorithm. It is known that the above exact linear-time algorithm leads to an nearly-linear-time approximation scheme for LP-type problems. This is summarized in the following. We refer the reader to Har-Peled (2011) for a more detailed treatment.
Lemma 2.3 (Har-Peled (2011), Ch. 15). LetA be some LP-type problem of combinatorial dimension κ ą 0, defined by some pair pH, wq, and let ε ą 0. There exists a randomized algorithm which given some instance F Ď H, with |F | “ n, outputs some basis B Ď F , that violates at most p1 ` εqk constraints in F , such that wpBq ď wpB1q, for any basis B1 violating at most k constraints in F , in time O ´ t0 ` ´ n` nmin ! logκ`1 n ε2κ , logκ`2 n kε2κ`2 )¯ pt1 ` t2q ¯ , where t0 is the time needed to compute
an arbitrary initial basis of A, and t1, t2, and t3 are upper bounds on the time needed to perform the basic operations (B0), (B1) and (B2) respectively. The algorithm succeeds with high probability.
For the special case of Mahalanobis metric learning, the corresponding algorithm is given in Algorithm 2. The approximation guarantee for this algorithm is summarized in 1.1. We can now give the proof of our main result.
Proof of Theorem 1.1. Follows immediately by Lemmas 2.2 and 2.3.
Algorithm 2 An approximation algorithm for Mahalanobis metric learning. procedure LPTML(F )
for i “ 0 to log1`ε n do pÐ p1` εq´i
for j “ 1 to logOpd 2q n do
subsample Fj Ď F , where each element is chosen independently with probability p Ai,j Ð Exact-LPTMLpFjq
end for end for return a solution out of tAi,jui,j , violating the minimum number of constraints in F
end procedure
Regularization. We now argue that the LP-type algorithm described above can be extended to handle certain types of regularization on the matrix A. In methods based on convex optimization, introducing regularizers that are convex functions can often be done easily. In our case, we cannot directly introduce a regularizing term in the objective function that is implicit in Algorithm 2. More specifically, let costpAq denote the total number of constraints of type (3) and (4) that A violates. Algorithm 2 approximately minimizes the objective function costpAq. A natural regularized version of Mahalanobis metric learning is to instead minimize the objective function
cost1pAq :“ costpAq ` η ¨ regpAq,
for some η ą 0, and regularizer regpAq. One typical choice is regpAq “ trpACq, for some matrix C P Rdˆd; the case C “ I corresponds to the trace norm (see Kulis et al. (2013)). We can extend the Algorithm 2 to handle any regularizer that can be expressed as a linear function on the entries of A, such as trpAq. The following summarizes the result. Theorem 2.4. Let regpAq be a linear function on the entries of A, with polynomially bounded coefficients. For any d P N, ε ą 0, there exists a randomized algorithm for learning d-dimensional Mahalanobis metric spaces, which given an instance that admits a solution A0 with cost1pA0q “ c˚, computes a solution A with cost1pAq ď p1`εqc˚, in time dOp1qnplog n{εqOpdq, with high probability.
Proof. If η ă εt, for sufficiently large constant t ą 0, since the coefficients in regpAq are polynomially bounded, it follows that the largest possible value of η ¨ regpAq is Opεq, and can thus be omitted without affecting the result. Similarly, if η ą p1{εqnt1 , for sufficiently large constant t1 ą 0, since there are at most `
n 2
˘
constraints, it follows that the term costpAq can be omitted form the objective. Therefore, we may assume w.l.o.g. that regpA0q P rεOp1q, p1{εqnOp1qs. We can guess some i “ Oplog n` logp1{εqq, such that regpA0q P pp1` εqi´1, p1` εqis. We modify the SDP used in the proof of Lemma 2.2 by introducing the constraint regpAq ď p1` εqi. Guessing the correct value of i requires Oplog n` logp1{εqq executions of Algorithm 2, which implies the running time bound.

heading: 3 PRACTICAL IMPROVEMENTS AND PARALLELIZATION
text: We now discuss some modifications of the algorithm described in the previous section that significantly improve its performance in practical scenarios, and have been integrated in our implementation.
Move-to-front and pivoting heuristics. We use heuristics that have been previously used in algorithms for linear programming Seidel (1990); Clarkson (1995), minimum enclosing ball in R3 Megiddo (1983), minimum enclosing ball and ellipsoid is Rd, for any fixed d P N Welzl (1991), as well as in fast implementations of minimum enclosing ball algorithms Gärtner (1999). The move-to-front heuristic keeps an ordered list of constraints which gets reorganized as the algorithm runs; when the algorithm finds a violation, it moves the violating constraint to the beginning of the list of the current sub-problem. The pivoting heuristic further improves performance by choosing to add to the basis the constraint that is “violated the most”. For instance, for similarity constraints, we pick the one that is mapped to the largest distance greater than u; for dissimilarity constraints, we pick the one that is mapped to the smallest distance less than `.
Approximate counting. The main loop of Algorithm 2 involves counting the number of violated constraints in each iteration. In problems involving a large number of constraints, we use approximate counting by only counting the number of violations within a sample of Oplog 1{εq constraints. We denote by LPTMLt for the version of the algorithm that performs a total of t iterations of the inner loop.
Early termination. A bottleneck of Algorithm 2 stems from the fact that the inner loop needs to be executed for logOpd
2q n iterations. In practice, we have observed that a significantly smaller number of iterations is needed to achieve high accuracy.
Parallelization. Algorithm 2 consists of several executions of the algorithm Exact-LPTML on independently sampled sub-problems. Therefore, Algorithm 2 can trivially be parallelized by distributing a different set of sub-problems to each machine, and returning the best solution found overall.

heading: 4 EXPERIMENTAL EVALUATION
text: We have implemented Algorithm 2, incorporating the practical improvements described in Section 3, and performed experiments on synthetic and real-world data sets. Our LPTML implementation and documentation can be found at the supplementary material1. We now describe the experimental setting and discuss the main findings.

heading: 4.1 EXPERIMENTAL SETTING
text: Classification task. Each data set used in the experiments consists of a set of labeled points in Rd. The label of each point indicates its class, and there is a constant number of classes. The set of similarity constraints S (respt. dissimilarity constraints D) is formed by uniformly sampling pairs of points in the same class (resp. from different classes). We use various algorithms to learn a Mahalanobis metric for a labeled input point set in Rd, given these constraints. The values u and ` are chosen as the 90th and 10th percentiles of all pairwise distances. We used 2-fold cross-validation: At the training phase we learn a Mahalanobis metric, and in the testing phase we use k-NN classification, with k “ 4, to evaluate the performance.
Data sets. We have tested our algorithm on the following synthetic and real-world data sets:
1. Real-world: We have tested the performance of our implementation on the Iris, Wine, Ionosphere and Soybean data sets from the UCI Machine Learning Repository2.
2. Synthetic: Next, we consider a synthetic data set that is constructed by first sampling a set of 100 points from a mixture of two Gaussians in R2, with identity covariance matrices, and with means p´3, 0q and p3, 0q respectively; we then apply a linear transformation that stretches the y axis by a factor of 40. This linear transformation reduces the accuracy of k-NN on the underlying Euclidean metric with k “ 4 from 1 to 0.68.
1A copy of our implementation is also available at https://drive.google.com/drive/folders/ 1XgABqyh8E1CoRGadh1KC5or7TBLgQkdl?usp=sharing
2https://archive.ics.uci.edu/ml/datasets.php
3. Synthetic + Adversarial Noise: We modify the above synthetic data set by introducing a small fraction of points in an adversarial manner, before applying the linear transformation. Figure 3b depicts the noise added as five points labeled as one of the classes, and sampled from a Gaussian with identity covariance matrix and mean p´100, 0q (Figure 3a).
Algorithms. We compare the performance of our algorithm against ITML and LMNN. We used the implementations provided by the authors of these works, with minor modifications.

heading: 4.2 RESULTS
text: Accuracy. Algorithm 2 minimizes the number of violated pairwise distance constraints. It is interesting to examine the effect of this objective function on the accuracy of k-NN classification. Figure 1 depicts this relationship for the Wine data set. We observe that, in general, as the number of iterations of the main loop of LPTML increases, the number of violated pairwise distance constraints decreases, and the accuracy of k-NN increases. This phenomenon remains consistent when we first perform PCA to d “ 4, 8, 12 dimensions.
Comparison to ITML and LMNN. We compared the accuracy obtained by LPTMLt, for t “ 2000 iterations, against ITML and LMNN. Table 1 summarizes the findings on the real-world and data sets and the synthetic data set without adversarial noise. We observe that LPTML achieves accuracy that is comparable to ITML and LMNN.
We observe that LPTML outperforms ITML and LMNN on the Synthetic + Adversarial Noise data set. This is due to the fact that the introduction of adversarial noise causes the relaxations used in ITML and LMNN to be biased towards contracting the x-axis. In contrast, the noise does not “fool” LPTML because it only changes the optimal accuracy by a small amount. The results are summarized in Figure 2.
The effect of dimension. The running time of LPTML grows with the dimension d. This is caused mostly by the fact that the combinatorial dimension of the underlying LP-type problem is Opd2q, and thus performing each basic operation requires solving an SDP with Opd2q constraints.
Figure 4 (appearing in the Appendix) depicts the effect of dimensionality in the running time, for t “ 100, . . . , 2000 iterations of the main loop. The data set used is Wine after performing PCA to d dimensions, for d “ 2, . . . , 13.
Parallel implementation. We implemented a massively parallel version of LPTML in the MapReduce model. The program maps different sub-problems of the main loop of LPTML to different machines. In the reduce step, we keep the result with the minimum number of constraint violations. The implementation uses the mrjob Yelp and Contributors (2019) package. For these experiments, we used Amazon cloud computing instances of type m4.xlarge, AMI 5.20.0 and configured with Hadoop. As expected, the training time decreases as the number of available processors increases (Figure 5 in the Appendix). All technical details about this implementation can be found in the parallel section of the documentation of our code.

heading: A ADDITIONAL FIGURES

Review summaries:
The paper proposes a method to handle Mahalanobis metric learning thorough linear programming. The authors consider the specific setup where examples are labelled as similar or dissimilar and the task is to find a mapping such that the feature-space distance between examples is i) smaller than a certain value if the examples are labelled as 'similar' and ii) greater than a possibly different value if the examples are labelled as 'dissimilar'. Arguments from the theory of linear programming are leveraged to define exact and approximated algorithms.