Repository: scikit-learn/scikit-learn
Pull request: [MRG+1] Fixes issue #3367 -> MCD fails on data with singular covariance matrix
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/master/CONTRIBUTING.md#Contributing-Pull-Requests
-->
#### Reference Issue
<!-- Example: Fixes #1234 -->

Fixes #3367 

See Issue #4635 (my previous fix two years ago)

I am attempting to fix #3367 again since it seems the original code diverged some since back when I originally attempted to fix it (pull request #4635). Before I think I was also patching to fix an issue with array shapes / sizes (see pull request #4517) that was reported in issue #4512 and eventually fixed in pull request #5234. 

#### What does this implement/fix? Explain your changes.

This fixes an issue in the Minimum Covariance Determinant. The implementation would throw a ValueError if data that produced a singular covariance matrix was passed in. In the original FastMCD paper (Rousseeuw and Van Driessen, 1999), there is no such requirement that the data not produce a singular matrix. In fact, a singular matrix means that you've found some set or subset of your data that produces a covariance determinant of zero, i.e. the absolute minimum covariance determinant. Since this implementation uses the `fast_logdet` procedure to get the determinant (by calculating the log of the determinant), that means that if the determinant is infinity (or negative infinity), we have found the MCD for that sample. The fix was basically to just return immediately instead of throwing a ValueError.

#### Any other comments?

I also added a test that should fail on master (see `sklearn/covariance/tests/test_robust_covariance/` and the procedure `test_issue3367`). I added a bunch of comments explaining what I was doing, but basically we can test that the fix works by:

1. Generating a set of data with a singular covariance matrix. 
2. Running this through the MCD fit.

I did this by generating a bunch of points that fit within a plane (X, Y, Z=0). It should fail on master by throwing a ValueError exception, and should just pass on my branch. Note that a UserWarning is annoyingly generated that says something along the lines of "The covariance matrix associated to your dataset is not full rank". I wasn't sure if I should remove this, since it could be a helpful error to tell people not to bother with the MCD if they already have clean data. The MCD is especially useful if you have outliers in your data (in my example, that would be points that are not on the plane I generated). In such a case, the above UserWarning won't be generated, and the outliers will actually be filtered out by the MCD process since outlier points will not produce a covariance matrix with the minimum determinant. 

Note that the fix I applied is still necessary, because it also solves the case where any subset of your original data has a singular covariance matrix. I just don't know if the UserWarning should be removed or changed to better reflect what I described above.

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Robust covaraince (fast_mcd) does not handle singular covariance matrices
When the c_step of fast_mcd encounters a covariance matrix with a 0 determinant, it raises an error. According to the paper (Rousseeuw & Van Driessen 1999), it shouldn't be the case. Technically, 0-determinant covariance is the optimal solution. A singular case requires a special treatment. (Section 5 of the paper). I'm new here, but I can contribute to this.

Hints from issue discussion:
Sorry for the long silence. If you are still interested a fix would be very welcome.

Hey, sorry for the long wait, I haven't taken the time to commit to this issue since I have been busy with final projects / exams. I submitted a PR which I believe fixes this issue (PR #4635), but does manage to produce some warnings with one of the sample datasets I used to test the fix. 

I believe the issue with that sample dataset (described in the PR) is different than what caused this issue. Overall, I think this issue is solved, as despite the warnings on some of the trials, this still produces the correct answer for exact fit scenarios.

I'm starting fresh to try and tackle this again, two years later (sorry, it seems I really delayed my old fix, which is now incompatible to an extent with the old code). I'm gonna try and write some test cases and submit a new PR. Sorry for being such a wreck about actually fixing this, I think I was quite close to getting everything correct last time.

For the record, I will also test whatever data I generate against the reference implementation in [LIBRA](http://wis.kuleuven.be/stat/robust/LIBRA) as well, which should let me know if I change things too much.

Commit message: Style corrections to pass flake runner

Implements the style corrections as mentioned in pull request #8328

Patch to review:
diff --git a/sklearn/covariance/robust_covariance.py b/sklearn/covariance/robust_covariance.py
index 29cbd52e183d3..a19a0ac036eb7 100644
--- a/sklearn/covariance/robust_covariance.py
+++ b/sklearn/covariance/robust_covariance.py
@@ -94,6 +94,7 @@ def _c_step(X, n_support, random_state, remaining_iterations=30,
             initial_estimates=None, verbose=False,
             cov_computation_method=empirical_covariance):
     n_samples, n_features = X.shape
+    dist = np.inf
 
     # Initialisation
     support = np.zeros(n_samples, dtype=bool)
@@ -114,11 +115,13 @@ def _c_step(X, n_support, random_state, remaining_iterations=30,
     X_support = X[support]
     location = X_support.mean(0)
     covariance = cov_computation_method(X_support)
+    precision = pinvh(covariance)
 
     # Iterative procedure for Minimum Covariance Determinant computation
     det = fast_logdet(covariance)
     previous_det = np.inf
-    while (det < previous_det) and (remaining_iterations > 0):
+    while (det < previous_det and remaining_iterations > 0
+            and not np.isinf(det)):
         # save old estimates values
         previous_location = location
         previous_covariance = covariance
@@ -140,14 +143,9 @@ def _c_step(X, n_support, random_state, remaining_iterations=30,
 
     previous_dist = dist
     dist = (np.dot(X - location, precision) * (X - location)).sum(axis=1)
-    # Catch computation errors
+    # Check if best fit already found (det => 0, logdet => -inf)
     if np.isinf(det):
-        raise ValueError(
-            "Singular covariance matrix. "
-            "Please check that the covariance matrix corresponding "
-            "to the dataset is full rank and that MinCovDet is used with "
-            "Gaussian-distributed data (or at least data drawn from a "
-            "unimodal, symmetric distribution.")
+        results = location, covariance, det, support, dist
     # Check convergence
     if np.allclose(det, previous_det):
         # c_step procedure converged
diff --git a/sklearn/covariance/tests/test_robust_covariance.py b/sklearn/covariance/tests/test_robust_covariance.py
index 27e423b410210..92b2229173809 100644
--- a/sklearn/covariance/tests/test_robust_covariance.py
+++ b/sklearn/covariance/tests/test_robust_covariance.py
@@ -4,6 +4,8 @@
 #
 # License: BSD 3 clause
 
+import itertools
+
 import numpy as np
 
 from sklearn.utils.testing import assert_almost_equal
@@ -92,6 +94,38 @@ def test_mcd_issue1127():
     mcd.fit(X)
 
 
+def test_mcd_issue3367():
+    # Check that MCD completes when the covariance matrix is singular
+    # i.e. one of the rows and columns are all zeros
+    rand_gen = np.random.RandomState(0)
+
+    # Think of these as the values for X and Y -> 10 values between -5 and 5
+    data_values = np.linspace(-5, 5, 10).tolist()
+    # Get the cartesian product of all possible coordinate pairs from above set
+    data = np.array(list(itertools.product(data_values, data_values)))
+
+    # Add a third column that's all zeros to make our data a set of point
+    # within a plane, which means that the covariance matrix will be singular
+    data = np.hstack((data, np.zeros((data.shape[0], 1))))
+
+    # The below line of code should raise an exception if the covariance matrix
+    # is singular As a further test, since we have points in XYZ, the principle
+    # components (Eigenvectors) of these directly relate to the geometry of the
+    # points. Since it's a plane, we should be able to test that the
+    # Eigenvector that corresponds to the smallest Eigenvalue is the plane
+    # normal, specifically [0, 0, 1], since everything is in the XY plane (as
+    # I've set it up above). To do this one would start by:
+    #
+    #     evals, evecs = np.linalg.eigh(mcd_fit.covariance_)
+    #     normal = evecs[:, np.argmin(evals)]
+    #
+    # After which we need to assert that our `normal` is equal to [0, 0, 1].
+    # Do note that there is floating point error associated with this, so it's
+    # best to subtract the two and then compare some small tolerance (e.g.
+    # 1e-12).
+    MinCovDet(random_state=rand_gen).fit(data)
+
+
 def test_outlier_detection():
     rnd = np.random.RandomState(0)
     X = rnd.randn(100, 10)


Code review comments:
diff_hunk: @@ -114,11 +115,13 @@ def _c_step(X, n_support, random_state, remaining_iterations=30,
     X_support = X[support]
     location = X_support.mean(0)
     covariance = cov_computation_method(X_support)
+    precision = pinvh(covariance)
path: sklearn/covariance/robust_covariance.py
text: @user1:
Can we please avoid duplicating this calculation in the usual case?

@author:
This causes problems if you put in data that has a singular covariance matrix. The following while loop is skipped, precision is never defined, and then it is referenced on line 145, which throws an error. To avoid the extra `pinvh` call, I can wrap the `precision = pinvh(covariance)` in a conditional after the initial determinant is calculated. I'll commit this, please let me know if that's sufficient.

diff_hunk: @@ -92,6 +94,38 @@ def test_mcd_issue1127():
     mcd.fit(X)
 
 
+def test_mcd_issue3367():
+    # Check that MCD completes when the covariance matrix is singular
+    # i.e. one of the rows and columns are all zeros
+    rand_gen = np.random.RandomState(0)
+
+    # Think of these as the values for X and Y -> 10 values between -5 and 5
+    data_values = np.linspace(-5, 5, 10).tolist()
+    # Get the cartesian product of all possible coordinate pairs from above set
+    data = np.array(list(itertools.product(data_values, data_values)))
+
+    # Add a third column that's all zeros to make our data a set of point
+    # within a plane, which means that the covariance matrix will be singular
+    data = np.hstack((data, np.zeros((data.shape[0], 1))))
+
+    # The below line of code should raise an exception if the covariance matrix
+    # is singular As a further test, since we have points in XYZ, the principle
path: sklearn/covariance/tests/test_robust_covariance.py
text: @user1:
. after singular.