scikit-learn · ogrisel · Nov 13, 2020 · Nov 3, 2020 · Nov 3, 2020 · Nov 4, 2020
diff --git a/doc/whats_new/v0.24.rst b/doc/whats_new/v0.24.rst
@@ -132,6 +132,9 @@ Changelog
   predictions for `est.transform(Y)` when the training data is single-target.
   :pr:`17095` by `Nicolas Hug`_.
 
+- |Fix| Increases the stability of :class:`cross_decomposition.CCA` :pr:`18746`
 Xk -= np.outer(x_scores, x_loadings) 
 Xk -= np.outer(x_scores, x_loadings) 
+  by `Thomas Fan`_.
+
 - |API| For :class:`cross_decomposition.NMF`,
   the `init` value, when 'init=None' and
   n_components <= min(n_samples, n_features) will be changed from

diff --git a/sklearn/cross_decomposition/_pls.py b/sklearn/cross_decomposition/_pls.py
@@ -12,7 +12,6 @@
 from scipy.linalg import pinv2, svd
 
 from ..base import BaseEstimator, RegressorMixin, TransformerMixin
-from ..base import _UnstableArchMixin
 from ..base import MultiOutputMixin
 from ..utils import check_array, check_consistent_length
 from ..utils.extmath import svd_flip
@@ -45,8 +44,8 @@ def _get_first_singular_vectors_power_method(X, Y, mode="A", max_iter=500,
         # As a result, and as detailed in the Wegelin's review, CCA (i.e. mode
         # B) will be unstable if n_features > n_samples or n_targets >
         # n_samples
-        X_pinv = pinv2(X, check_finite=False)
-        Y_pinv = pinv2(Y, check_finite=False)
+        X_pinv = pinv2(X, check_finite=False, cond=10*eps)
+        Y_pinv = pinv2(Y, check_finite=False, cond=10*eps)
 
     for i in range(max_iter):
         if mode == "B":
@@ -683,7 +682,7 @@ def __init__(self, n_components=2, *, scale=True, algorithm="nipals",
             max_iter=max_iter, tol=tol, copy=copy)
 
 
-class CCA(_UnstableArchMixin, _PLS):
+class CCA(_PLS):
     """Canonical Correlation Analysis, also known as "Mode B" PLS.
 
     Read more in the :ref:`User Guide <cross_decomposition>`.

diff --git a/sklearn/cross_decomposition/tests/test_pls.py b/sklearn/cross_decomposition/tests/test_pls.py
@@ -382,54 +382,57 @@ def test_copy(Est):
                               pls.predict(X.copy(), copy=False))
 
 
-@pytest.mark.xfail
-@pytest.mark.parametrize('Est', (CCA, PLSCanonical, PLSRegression, PLSSVD))
-def test_scale_and_stability(Est):
-    # scale=True is equivalent to scale=False on centered/scaled data
-    # This allows to check numerical stability over platforms as well
-
+def _generate_test_scale_and_stability_datasets():
+    """Generate dataset for test_scale_and_stability"""
+    # dataset for non-regression 7818
     rng = np.random.RandomState(0)
-
-    d = load_linnerud()
-    X1 = d.data
-    Y1 = d.target
-    # causes X[:, -1].std() to be zero
-    X1[:, -1] = 1.0
-
-    # From bug #2821
-    # Test with X2, Y2 s.t. clf.x_score[:, 1] == 0, clf.y_score[:, 1] == 0
-    # This test robustness of algorithm when dealing with value close to 0
-    X2 = np.array([[0., 0., 1.],
-                   [1., 0., 0.],
-                   [2., 2., 2.],
-                   [3., 5., 4.]])
-    Y2 = np.array([[0.1, -0.2],
-                   [0.9, 1.1],
-                   [6.2, 5.9],
-                   [11.9, 12.3]])
-
-    # Non-regression for https://github.com/scikit-learn/scikit-learn/pull/7819
     n_samples = 1000
     n_targets = 5
     n_features = 10
     Q = rng.randn(n_targets, n_features)
-    Y3 = rng.randn(n_samples, n_targets)
-    X3 = np.dot(Y3, Q) + 2 * rng.randn(n_samples, n_features) + 1
-    X3 *= 1000
-
-    for (X, Y) in [(X1, Y1), (X2, Y2), (X3, Y3)]:
-        X_std = X.std(axis=0, ddof=1)
-        X_std[X_std == 0] = 1
-        Y_std = Y.std(axis=0, ddof=1)
-        Y_std[Y_std == 0] = 1
-        X_s = (X - X.mean(axis=0)) / X_std
-        Y_s = (Y - Y.mean(axis=0)) / Y_std
-
-        X_score, Y_score = Est(scale=True).fit_transform(X, Y)
-        X_s_score, Y_s_score = Est(scale=False).fit_transform(X_s, Y_s)
-
-        assert_array_almost_equal(X_s_score, X_score)
-        assert_array_almost_equal(Y_s_score, Y_score)
+    Y = rng.randn(n_samples, n_targets)
+    X = np.dot(Y, Q) + 2 * rng.randn(n_samples, n_features) + 1
+    X *= 1000
+    yield X, Y
+
+    # Data set where one of the features is constaint
+    X, Y = load_linnerud(return_X_y=True)
+    # causes X[:, -1].std() to be zero
+    X[:, -1] = 1.0
+    yield X, Y
+
+    X = np.array([[0., 0., 1.],
+                  [1., 0., 0.],
+                  [2., 2., 2.],
+                  [3., 5., 4.]])
+    Y = np.array([[0.1, -0.2],
+                  [0.9, 1.1],
+                  [6.2, 5.9],
+                  [11.9, 12.3]])
+    yield X, Y
+
+    # Seeds that provide a non-regression test for #18746, where CCA fails
+    seeds = [530, 741]
+    for seed in seeds:
+        rng = np.random.RandomState(seed)
+        X = rng.randn(4, 3)
+        Y = rng.randn(4, 2)
+        yield X, Y
+
+
+@pytest.mark.parametrize('Est', (CCA, PLSCanonical, PLSRegression, PLSSVD))
+@pytest.mark.parametrize('X, Y', _generate_test_scale_and_stability_datasets())
+def test_scale_and_stability(Est, X, Y):
+    """scale=True is equivalent to scale=False on centered/scaled data
+    This allows to check numerical stability over platforms as well"""
+
+    X_s, Y_s, *_ = _center_scale_xy(X, Y)
+
+    X_score, Y_score = Est(scale=True).fit_transform(X, Y)
+    X_s_score, Y_s_score = Est(scale=False).fit_transform(X_s, Y_s)
+
+    assert_allclose(X_s_score, X_score, atol=1e-4)
+    assert_allclose(Y_s_score, Y_score, atol=1e-4)
 
 
 @pytest.mark.parametrize('Est', (PLSSVD, PLSCanonical, CCA))