FIX numerical stability in GMM with eigh sampling

sklearn/mixture/gmm.py

@@ -98,8 +98,10 @@ def sample_gaussian(mean, covar, covariance_type='diag', n_samples=1,
         rand = np.dot(np.diag(np.sqrt(covar)), rand)
     else:
         s, U = linalg.eigh(covar)
-        sqrt_covar = U * np.sqrt(s)
-        rand = np.dot(sqrt_covar, rand)
+        np.abs(s, out=s)        # get rid of tiny negatives
+        np.sqrt(s, out=s)
+        U *= s
+        rand = np.dot(U, rand)
 
     return (rand.T + mean).T
 

sklearn/mixture/tests/test_gmm.py

@@ -45,6 +45,14 @@ def test_sample_gaussian():
     assert_true(np.allclose(samples.mean(axis), mu, atol=1.3))
     assert_true(np.allclose(np.cov(samples), cv, atol=2.5))
 
+    # Numerical stability check: in SciPy 0.12.0 at least, eigh may return
+    # tiny negative values in its second return value.
+    from sklearn.mixture import sample_gaussian
+    x = sample_gaussian([0, 0], [[4, 3], [1, .1]],
+                        covariance_type='full', random_state=42)
+    print(x)
+    assert_true(np.isfinite(x).all())
+
 
 def _naive_lmvnpdf_diag(X, mu, cv):
     # slow and naive implementation of lmvnpdf