[MRG] FIX LMNN gradient and cost function

test/metric_learn_test.py

@@ -2,7 +2,7 @@
 import re
 import pytest
 import numpy as np
-from scipy.optimize import check_grad
+from scipy.optimize import check_grad, approx_fprime
 from six.moves import xrange
 from sklearn.metrics import pairwise_distances
 from sklearn.datasets import load_iris, make_classification, make_regression
@@ -21,7 +21,7 @@
                           RCA_Supervised, MMC_Supervised, SDML)
 # Import this specially for testing.
 from metric_learn.constraints import wrap_pairs
-from metric_learn.lmnn import python_LMNN
+from metric_learn.lmnn import python_LMNN, _sum_outer_products
 
 
 def class_separation(X, labels):
@@ -120,6 +120,61 @@ def test_iris(self):
                               self.iris_labels)
       self.assertLess(csep, 0.25)
 
+  def test_loss_grad_lbfgs(self):
+    """Test gradient of loss function
+    Assert that the gradient is almost equal to its finite differences
+    approximation.
+    """
+    rng = np.random.RandomState(42)
+    X, y = make_classification(random_state=rng)
+    L = rng.randn(rng.randint(1, X.shape[1] + 1), X.shape[1])
+    lmnn = LMNN()
+
+    k = lmnn.k
+    reg = lmnn.regularization
+
+    X, y = lmnn._prepare_inputs(X, y, dtype=float,
+                                ensure_min_samples=2)
+    num_pts, num_dims = X.shape
+    unique_labels, label_inds = np.unique(y, return_inverse=True)
+    lmnn.labels_ = np.arange(len(unique_labels))
+    lmnn.transformer_ = np.eye(num_dims)
+
+    target_neighbors = lmnn._select_targets(X, label_inds)
+    impostors = lmnn._find_impostors(target_neighbors[:, -1], X, label_inds)
+
+    # sum outer products
+    dfG = _sum_outer_products(X, target_neighbors.flatten(),
+                              np.repeat(np.arange(X.shape[0]), k))
+    df = np.zeros_like(dfG)
+
+    # storage
+    a1 = [None]*k
+    a2 = [None]*k
+    for nn_idx in xrange(k):
+      a1[nn_idx] = np.array([])
+      a2[nn_idx] = np.array([])
+
+    # initialize L
+
+    def fun(L):
+        return lmnn._loss_grad(X, L.reshape(-1, X.shape[1]), dfG, impostors, 1,
+                               k, reg,
+                               target_neighbors, df, a1, a2)[1]
+
+    def grad(L):
+        return lmnn._loss_grad(X, L.reshape(-1, X.shape[1]), dfG, impostors,
+                               1, k, reg,
+          target_neighbors, df, a1, a2)[0].ravel()
+
+    # compute relative error
+    epsilon = np.sqrt(np.finfo(float).eps)
+    rel_diff = (check_grad(fun, grad, L.ravel()) /
+                           np.linalg.norm(approx_fprime(L.ravel(), fun,
+                                                       epsilon)))
+               # np.linalg.norm(grad(L))
+    np.testing.assert_almost_equal(rel_diff, 0., decimal=5)
+
 
 def test_convergence_simple_example(capsys):
   # LMNN should converge on this simple example, which it did not with
@@ -421,8 +476,10 @@ def grad(M):
       return nca._loss_grad_lbfgs(M, X, mask)[1].ravel()
 
     # compute relative error
-    rel_diff = check_grad(fun, grad, M.ravel()) / np.linalg.norm(grad(M))
-    np.testing.assert_almost_equal(rel_diff, 0., decimal=6)
+    epsilon = np.sqrt(np.finfo(float).eps)
+    rel_diff = (check_grad(fun, grad, M.ravel()) /
+                np.linalg.norm(approx_fprime(M.ravel(), fun, epsilon)))
+    np.testing.assert_almost_equal(rel_diff, 0., decimal=10)
 
   def test_simple_example(self):
     """Test on a simple example.