diff --git a/metric_learn/lsml.py b/metric_learn/lsml.py
index 1d66cbc0..7d81fdf5 100644
--- a/metric_learn/lsml.py
+++ b/metric_learn/lsml.py
@@ -39,7 +39,7 @@ def __init__(self, tol=1e-3, max_iter=1000, prior=None, verbose=False,
     tol : float, optional
     max_iter : int, optional
     prior : (d x d) matrix, optional
-        guess at a metric [default: inv(covariance(X))]
+        guess at a metric [default: pinvh(covariance(X))]
     verbose : bool, optional
         if True, prints information while learning
     preprocessor : array-like, shape=(n_samples, n_features) or callable
@@ -70,10 +70,10 @@ def _fit(self, quadruplets, weights=None):
       X = np.vstack({tuple(row) for row in
                      quadruplets.reshape(-1, quadruplets.shape[2])})
       prior_inv = np.atleast_2d(np.cov(X, rowvar=False))
-      M = np.linalg.inv(prior_inv)
+      M = scipy.linalg.pinvh(prior_inv)
     else:
       M = self.prior
-      prior_inv = np.linalg.inv(self.prior)
+      prior_inv = scipy.linalg.pinvh(self.prior)
 
     step_sizes = np.logspace(-10, 0, 10)
     # Keep track of the best step size and the loss at that step.
@@ -126,7 +126,7 @@ def _total_loss(self, metric, vab, vcd, prior_inv):
     return self._comparison_loss(metric, vab, vcd) + reg_loss
 
   def _gradient(self, metric, vab, vcd, prior_inv):
-    dMetric = prior_inv - np.linalg.inv(metric)
+    dMetric = prior_inv - scipy.linalg.pinvh(metric)
     dabs = np.sum(vab.dot(metric) * vab, axis=1)
     dcds = np.sum(vcd.dot(metric) * vcd, axis=1)
     violations = dabs > dcds
diff --git a/test/metric_learn_test.py b/test/metric_learn_test.py
index a785d60d..518d9c1e 100644
--- a/test/metric_learn_test.py
+++ b/test/metric_learn_test.py
@@ -74,6 +74,18 @@ def test_deprecation_num_labeled(self):
            'removed in 0.6.0')
     assert_warns_message(DeprecationWarning, msg, lsml_supervised.fit, X, y)
 
+  def test_singular_covariance_does_not_diverge(self):
+    # Test that LSML does not diverge when using the covariance prior and
+    # when this covariance has null eigenvalues (See
+    # https://github.com/metric-learn/metric-learn/issues/202)
+    # TODO: remove when # 195 is merged
+    rng = np.random.RandomState(42)
+    X, y = load_iris(return_X_y=True)
+    # we add a feature that is a linear combination of the two first ones
+    X = np.concatenate([X, X[:, :2].dot(rng.randn(2, 1))], axis=-1)
+    lsml = LSML_Supervised()
+    lsml.fit(X, y, random_state=rng)
+
 
 class TestITML(MetricTestCase):
   def test_iris(self):