[MRG] FIX: make proposal for sdml formulation

.travis.yml

@@ -4,8 +4,14 @@ cache: pip
 python:
   - "2.7"
   - "3.4"
+  - "3.6"
 before_install:
+  - sudo apt-get install liblapack-dev
   - pip install --upgrade pip
   - pip install wheel
-  - pip install numpy scipy scikit-learn
+  - pip install cython numpy scipy scikit-learn
+  - if [[ ($TRAVIS_PYTHON_VERSION == "3.6") ||
+          ($TRAVIS_PYTHON_VERSION == "2.7")]]; then
+              pip install git+https://github.com/skggm/skggm.git@a0ed406586c4364ea3297a658f415e13b5cbdaf8;
+    fi
 script: pytest test

README.rst

@@ -20,7 +20,7 @@ Metric Learning algorithms in Python.
 **Dependencies**
 
 -  Python 2.7+, 3.4+
--  numpy, scipy, scikit-learn
+-  numpy, scipy, scikit-learn, and skggm (commit `a0ed406 <https://github.com/skggm/skggm/commit/a0ed406586c4364ea3297a658f415e13b5cbdaf8>`_) for `SDML`
 -  (for running the examples only: matplotlib)
 
 **Installation/Setup**

doc/getting_started.rst

@@ -15,7 +15,7 @@ Alternately, download the source repository and run:
 **Dependencies**
 
 -  Python 2.7+, 3.4+
--  numpy, scipy, scikit-learn
+-  numpy, scipy, scikit-learn, and skggm (commit `a0ed406 <https://github.com/skggm/skggm/commit/a0ed406586c4364ea3297a658f415e13b5cbdaf8>`_) for `SDML`
 -  (for running the examples only: matplotlib)
 
 **Notes**

metric_learn/_util.py

@@ -15,6 +15,14 @@ def vector_norm(X):
     return np.linalg.norm(X, axis=1)
 
 
+def has_installed_skggm():
+  try:
+    import inverse_covariance
+    return True
+  except ImportError:
+    return False
+
+
 def check_input(input_data, y=None, preprocessor=None,
                 type_of_inputs='classic', tuple_size=None, accept_sparse=False,
                 dtype='numeric', order=None,

metric_learn/constraints.py

@@ -96,6 +96,6 @@ def wrap_pairs(X, constraints):
   c = np.array(constraints[2])
   d = np.array(constraints[3])
   constraints = np.vstack((np.column_stack((a, b)), np.column_stack((c, d))))
-  y = np.vstack([np.ones((len(a), 1)), - np.ones((len(c), 1))])
+  y = np.hstack([np.ones((len(a),)), - np.ones((len(c),))])
   pairs = X[constraints]
   return pairs, y

metric_learn/sdml.py

@@ -12,12 +12,15 @@
 import warnings
 import numpy as np
 from sklearn.base import TransformerMixin
-from sklearn.covariance import graph_lasso
-from sklearn.utils.extmath import pinvh
+from scipy.linalg import pinvh
+from sklearn.covariance import graphical_lasso
+from sklearn.exceptions import ConvergenceWarning
 
 from .base_metric import MahalanobisMixin, _PairsClassifierMixin
 from .constraints import Constraints, wrap_pairs
-from ._util import transformer_from_metric
+from ._util import transformer_from_metric, has_installed_skggm
+if has_installed_skggm():
+  from inverse_covariance import quic
 
 
 class _BaseSDML(MahalanobisMixin):
@@ -52,24 +55,50 @@ def __init__(self, balance_param=0.5, sparsity_param=0.01, use_cov=True,
     super(_BaseSDML, self).__init__(preprocessor)
 
   def _fit(self, pairs, y):
+    if not has_installed_skggm():
+      msg = ("Warning, skggm is not installed, so SDML will use "
+             "scikit-learn's graphical_lasso method. It can fail to converge"
+             "on some non SPD matrices where skggm would converge. If so, "
+             "try to install skggm. (see the README.md for the right "
+             "version.)")
+      warnings.warn(msg)
+    else:
+      print("SDML will use skggm's solver.")
     pairs, y = self._prepare_inputs(pairs, y,
                                     type_of_inputs='tuples')
 
     # set up prior M
     if self.use_cov:
       X = np.vstack({tuple(row) for row in pairs.reshape(-1, pairs.shape[2])})
-      M = pinvh(np.atleast_2d(np.cov(X, rowvar = False)))
+      prior = pinvh(np.atleast_2d(np.cov(X, rowvar=False)))
     else:
-      M = np.identity(pairs.shape[2])
+      prior = np.identity(pairs.shape[2])
     diff = pairs[:, 0] - pairs[:, 1]
     loss_matrix = (diff.T * y).dot(diff)
-    P = M + self.balance_param * loss_matrix
-    emp_cov = pinvh(P)
-    # hack: ensure positive semidefinite
-    emp_cov = emp_cov.T.dot(emp_cov)
-    _, M = graph_lasso(emp_cov, self.sparsity_param, verbose=self.verbose)
-
-    self.transformer_ = transformer_from_metric(M)
+    emp_cov = pinvh(prior) + self.balance_param * loss_matrix
+
+    # our initialization will be the matrix with emp_cov's eigenvalues,
+    # with a constant added so that they are all positive (plus an epsilon
+    # to ensure definiteness). This is empirical.
+    w, V = np.linalg.eigh(emp_cov)
+    if any(w < 0.):
+      warnings.warn("Warning, the input matrix of graphical lasso is not "
+                    "positive semi-definite (PSD). The algorithm may diverge, "
+                    "and lead to degenerate solutions. "
+                    "To prevent that, try to decrease the balance parameter "
+                    "`balance_param` and/or to set use_covariance=False.",
+                    ConvergenceWarning)
+    sigma0 = (V * (w - min(0, np.min(w)) + 1e-10)).dot(V.T)
+    if has_installed_skggm():
+      theta0 = pinvh(sigma0)
+      M, _, _, _, _, _ = quic(emp_cov, lam=self.sparsity_param,
+                              msg=self.verbose,
+                              Theta0=theta0, Sigma0=sigma0)
+    else:
+      _, M = graphical_lasso(emp_cov, alpha=self.sparsity_param,
+                             verbose=self.verbose,
+                             cov_init=sigma0)
+    self.transformer_ = transformer_from_metric(np.atleast_2d(M))
     return self
 
 

setup.py

@@ -38,6 +38,7 @@
       extras_require=dict(
           docs=['sphinx', 'shinx_rtd_theme', 'numpydoc'],
           demo=['matplotlib'],
+          sdml=['skggm']
       ),
       test_suite='test',
       keywords=[

test/metric_learn_test.py

@@ -11,10 +11,11 @@
 from sklearn.exceptions import ConvergenceWarning
 from sklearn.utils.validation import check_X_y
 
-from metric_learn import (
-    LMNN, NCA, LFDA, Covariance, MLKR, MMC,
-    LSML_Supervised, ITML_Supervised, SDML_Supervised, RCA_Supervised, MMC_Supervised)
+from metric_learn import (LMNN, NCA, LFDA, Covariance, MLKR, MMC,
+                          LSML_Supervised, ITML_Supervised, SDML_Supervised,
+                          RCA_Supervised, MMC_Supervised, SDML)
 # Import this specially for testing.
+from metric_learn._util import has_installed_skggm
 from metric_learn.constraints import wrap_pairs
 from metric_learn.lmnn import python_LMNN
 
@@ -148,27 +149,142 @@ def test_no_twice_same_objective(capsys):
 
 
 class TestSDML(MetricTestCase):
-  def test_iris(self):
-    # Note: this is a flaky test, which fails for certain seeds.
-    # TODO: un-flake it!
-    rs = np.random.RandomState(5555)
-
-    sdml = SDML_Supervised(num_constraints=1500)
-    sdml.fit(self.iris_points, self.iris_labels, random_state=rs)
-    csep = class_separation(sdml.transform(self.iris_points), self.iris_labels)
-    self.assertLess(csep, 0.25)
 
-  def test_deprecation_num_labeled(self):
-    # test that a deprecation message is thrown if num_labeled is set at
-    # initialization
-    # TODO: remove in v.0.6
-    X = np.array([[0, 0], [0, 1], [2, 0], [2, 1]])
-    y = np.array([1, 0, 1, 0])
-    sdml_supervised = SDML_Supervised(num_labeled=np.inf)
-    msg = ('"num_labeled" parameter is not used.'
-           ' It has been deprecated in version 0.5.0 and will be'
-           'removed in 0.6.0')
-    assert_warns_message(DeprecationWarning, msg, sdml_supervised.fit, X, y)
+    @pytest.mark.skipif(has_installed_skggm(),
+                        reason="The warning will be thrown only if skggm is "
+                               "not installed.")
+    def test_raises_warning_msg_not_installed_skggm(self):
+      """Tests that the right warning message is raised if someone tries to
+      use SDML but has not installed skggm"""
+      # TODO: remove if we don't need skggm anymore
+      pairs = np.array([[[-10., 0.], [10., 0.]], [[0., -55.], [0., -60]]])
+      y_pairs = [1, -1]
+      X, y = make_classification(random_state=42)
+      sdml = SDML()
+      sdml_supervised = SDML_Supervised(use_cov=False, balance_param=1e-5)
+      msg = ("Warning, skggm is not installed, so SDML will use "
+             "scikit-learn's graphical_lasso method. It can fail to converge"
+             "on some non SPD matrices where skggm would converge. If so, "
+             "try to install skggm. (see the README.md for the right "
+             "version.)")
+      with pytest.warns(None) as record:
+        sdml.fit(pairs, y_pairs)
+      assert str(record[0].message) == msg
+      with pytest.warns(None) as record:
+        sdml_supervised.fit(X, y)
+      assert str(record[0].message) == msg
+
+    @pytest.mark.skipif(not has_installed_skggm(),
+                        reason="It's only in the case where skggm is installed"
+                               "that no warning should be thrown.")
+    def test_raises_no_warning_installed_skggm(self):
+      # otherwise we should be able to instantiate and fit SDML and it
+      # should raise no warning
+      pairs = np.array([[[-10., 0.], [10., 0.]], [[0., -55.], [0., -60]]])
+      y_pairs = [1, -1]
+      X, y = make_classification(random_state=42)
+      with pytest.warns(None) as record:
+        sdml = SDML()
+        sdml.fit(pairs, y_pairs)
+      assert len(record) == 0
+      with pytest.warns(None) as record:
+        sdml = SDML_Supervised(use_cov=False, balance_param=1e-5)
+        sdml.fit(X, y)
+      assert len(record) == 0
+
+    def test_iris(self):
+      # Note: this is a flaky test, which fails for certain seeds.
+      # TODO: un-flake it!
+      rs = np.random.RandomState(5555)
+
+      sdml = SDML_Supervised(num_constraints=1500, use_cov=False,
+                             balance_param=5e-5)
+      sdml.fit(self.iris_points, self.iris_labels, random_state=rs)
+      csep = class_separation(sdml.transform(self.iris_points),
+                              self.iris_labels)
+      self.assertLess(csep, 0.22)
+
+    def test_deprecation_num_labeled(self):
+      # test that a deprecation message is thrown if num_labeled is set at
+      # initialization
+      # TODO: remove in v.0.6
+      X, y = make_classification(random_state=42)
+      sdml_supervised = SDML_Supervised(num_labeled=np.inf, use_cov=False,
+                                        balance_param=5e-5)
+      msg = ('"num_labeled" parameter is not used.'
+             ' It has been deprecated in version 0.5.0 and will be'
+             'removed in 0.6.0')
+      assert_warns_message(DeprecationWarning, msg, sdml_supervised.fit, X, y)
+
+    def test_sdml_raises_warning_non_psd(self):
+      """Tests that SDML raises a warning on a toy example where we know the
+      pseudo-covariance matrix is not PSD"""
+      pairs = np.array([[[-10., 0.], [10., 0.]], [[0., 50.], [0., -60]]])
+      y = [1, -1]
+      sdml = SDML(use_cov=True, sparsity_param=0.01, balance_param=0.5)
+      msg = ("Warning, the input matrix of graphical lasso is not "
+             "positive semi-definite (PSD). The algorithm may diverge, "
+             "and lead to degenerate solutions. "
+             "To prevent that, try to decrease the balance parameter "
+             "`balance_param` and/or to set use_covariance=False.")
+      with pytest.warns(ConvergenceWarning) as raised_warning:
+        try:
+          sdml.fit(pairs, y)
+        except Exception:
+          pass
+      # we assert that this warning is in one of the warning raised by the
+      # estimator
+      assert msg in list(map(lambda w: str(w.message), raised_warning))
+
+    def test_sdml_converges_if_psd(self):
+      """Tests that sdml converges on a simple problem where we know the
+      pseudo-covariance matrix is PSD"""
+      pairs = np.array([[[-10., 0.], [10., 0.]], [[0., -55.], [0., -60]]])
+      y = [1, -1]
+      sdml = SDML(use_cov=True, sparsity_param=0.01, balance_param=0.5)
+      sdml.fit(pairs, y)
+      assert np.isfinite(sdml.get_mahalanobis_matrix()).all()
+
+    @pytest.mark.skipif(not has_installed_skggm(),
+                        reason="sklearn's graphical_lasso can sometimes not "
+                               "work on some non SPD problems. We test that "
+                               "is works only if skggm is installed.")
+    def test_sdml_works_on_non_spd_pb_with_skggm(self):
+      """Test that SDML works on a certain non SPD problem on which we know
+      it should work, but scikit-learn's graphical_lasso does not work"""
+      X, y = load_iris(return_X_y=True)
+      sdml = SDML_Supervised(balance_param=0.5, sparsity_param=0.01,
+                             use_cov=True)
+      sdml.fit(X, y)
+
+
+@pytest.mark.skipif(not has_installed_skggm(),
+                    reason='The message should be printed only if skggm is '
+                           'installed.')
+def test_verbose_has_installed_skggm_sdml(capsys):
+  # Test that if users have installed skggm, a message is printed telling them
+  # skggm's solver is used (when they use SDML)
+  # TODO: remove if we don't need skggm anymore
+  pairs = np.array([[[-10., 0.], [10., 0.]], [[0., -55.], [0., -60]]])
+  y_pairs = [1, -1]
+  sdml = SDML()
+  sdml.fit(pairs, y_pairs)
+  out, _ = capsys.readouterr()
+  assert "SDML will use skggm's solver." in out
+
+
+@pytest.mark.skipif(not has_installed_skggm(),
+                    reason='The message should be printed only if skggm is '
+                           'installed.')
+def test_verbose_has_installed_skggm_sdml_supervised(capsys):
+  # Test that if users have installed skggm, a message is printed telling them
+  # skggm's solver is used (when they use SDML_Supervised)
+  # TODO: remove if we don't need skggm anymore
+  X, y = make_classification(random_state=42)
+  sdml = SDML_Supervised()
+  sdml.fit(X, y)
+  out, _ = capsys.readouterr()
+  assert "SDML will use skggm's solver." in out
 
 
 class TestNCA(MetricTestCase):

test/test_base_metric.py

@@ -4,6 +4,7 @@
 import numpy as np
 from sklearn import clone
 from sklearn.utils.testing import set_random_state
+
 from test.test_utils import ids_metric_learners, metric_learners
 
 
@@ -55,7 +56,8 @@ def test_lsml(self):
   def test_sdml(self):
     self.assertEqual(str(metric_learn.SDML()),
                      "SDML(balance_param=0.5, preprocessor=None, "
-                     "sparsity_param=0.01, use_cov=True,\n   verbose=False)")
+                     "sparsity_param=0.01, use_cov=True,\n   "
+                     "verbose=False)")
     self.assertEqual(str(metric_learn.SDML_Supervised()), """
 SDML_Supervised(balance_param=0.5, num_constraints=None,
         num_labeled='deprecated', preprocessor=None, sparsity_param=0.01,

test/test_fit_transform.py

@@ -1,11 +1,13 @@
+import pytest
 import unittest
 import numpy as np
 from sklearn.datasets import load_iris
 from numpy.testing import assert_array_almost_equal
 
 from metric_learn import (
     LMNN, NCA, LFDA, Covariance, MLKR,
-    LSML_Supervised, ITML_Supervised, SDML_Supervised, RCA_Supervised, MMC_Supervised)
+    LSML_Supervised, ITML_Supervised, SDML_Supervised, RCA_Supervised,
+    MMC_Supervised)
 
 
 class TestFitTransform(unittest.TestCase):
@@ -62,12 +64,14 @@ def test_lmnn(self):
 
   def test_sdml_supervised(self):
     seed = np.random.RandomState(1234)
-    sdml = SDML_Supervised(num_constraints=1500)
+    sdml = SDML_Supervised(num_constraints=1500, balance_param=1e-5,
+                           use_cov=False)
     sdml.fit(self.X, self.y, random_state=seed)
     res_1 = sdml.transform(self.X)
 
     seed = np.random.RandomState(1234)
-    sdml = SDML_Supervised(num_constraints=1500)
+    sdml = SDML_Supervised(num_constraints=1500, balance_param=1e-5,
+                           use_cov=False)
     res_2 = sdml.fit_transform(self.X, self.y, random_state=seed)
 
     assert_array_almost_equal(res_1, res_2)

test/test_mahalanobis_mixin.py

@@ -96,7 +96,7 @@ def check_is_distance_matrix(pairwise):
   assert np.array_equal(pairwise, pairwise.T)  # symmetry
   assert (pairwise.diagonal() == 0).all()  # identity
   # triangular inequality
-  tol = 1e-15
+  tol = 1e-12
   assert (pairwise <= pairwise[:, :, np.newaxis] +
           pairwise[:, np.newaxis, :] + tol).all()
 
@@ -281,5 +281,15 @@ def test_transformer_is_2D(estimator, build_dataset):
 
   # test that it works for 1 feature
   trunc_data = input_data[..., :1]
+  # we drop duplicates that might have been formed, i.e. of the form
+  # aabc or abcc or aabb for quadruplets, and aa for pairs.
+  slices = {4: [slice(0, 2), slice(2, 4)], 2: [slice(0, 2)]}
+  if trunc_data.ndim == 3:
+    for slice_idx in slices[trunc_data.shape[1]]:
+      pairs = trunc_data[:, slice_idx, :]
+      diffs = pairs[:, 1, :] - pairs[:, 0, :]
+      to_keep = np.nonzero(diffs.ravel())
+      trunc_data = trunc_data[to_keep]
+      labels = labels[to_keep]
   model.fit(trunc_data, labels)
   assert model.transformer_.shape == (1, 1)  # the transformer must be 2D