Repository: scikit-learn/scikit-learn
Pull request: MAINT Make fastica call FastICA
~~NOMRG, I think there's a bug in fastica that needs to be fixed first~~

Closes #14988 

This PR makes `fastica` call `FastICA`. Following the discussion in #14988 I have created a private `_unmixing` attribute and clarified the docstring of the `components_` attribute.

Problem statement:
Inconsistency between fastica and FastICA
`fastica` returns `W` which is, according to the docs "the unmixing matrix"

`FastICA` has a `components_` attribute which is also "the unmixing matrix"

But they are not equal when `whiten` is True:

```py
if whiten == True:
	self.components_ = np.dot(unmixing, whitening)
else:
    self.components_ = unmixing
```

I don't know whether this is a bug (I hope it is), but if it's not the doc should be updated.

ping @agramfort

Hints from issue discussion:
sources = np.dot(self.components_, X.T) with the current code (always and
no matter what is the value of whiten)

It does apply the unmixing transformation to the data (taking into account
the whitening done prior to calling the fastica code).

are you still confused?

Sorry I'm still confused. This is unrelated to `sources` as far as I understand.

The problem is that both `W` from `fastica()` and `components_` from `FastICA` are documented as "the unmixing matrix", but they are clearly not the same thing. They don't even have the same shape.

`est.components_` is equal to `W.dot(est.whitening_)` when whiten is True.

<details>

```py
from sklearn.datasets import load_digits
from sklearn.decomposition import FastICA, fastica
import numpy as np

X, _ = load_digits(return_X_y=True)
whiten = True

est = FastICA(n_components=7, random_state=0, whiten=whiten).fit(X)
print(est.components_[0])  # print first comp

K, W, _ = fastica(X, n_components=7, random_state=0, whiten=whiten)
print(W.dot(est.whitening_)[0])  # equal to first comp

# according to the doc W is supposed to be the same as components_ but that's
# clearly not the case.
```
</details>
fastica estimates an unmixing matrix (square) from whitened data
(K.dot(X.T))

now if I give you X and you want to know what is the linear operator to go
from X
to independant sources you are asking me what is the unmixing matrix (it
unmixes
the X to get unmixed / independant sources). This matrix that we call
components_
in FastICA contains the whitener times the fastica output obtained from
whtened data.

if you want you can document components_ as the "linear operator to apply
to the data
to get the independent sources"

OK so:

1. the docs are wrong and must be updated

2. `fastica` and `FastICA` do not have the same API: it is impossible to retrieve `W` (as returned by `fastica`) if you use `FastICA`.

Feel free to disagree but I see that as a problem.
>
> OK so:
>
>    1.
>
>    the docs are wrong and must be updated
>
>
let's say they must be clarified....


>
>    1.
>
>    fastica and FastICA do not have the same API: it is impossible to
>    retrieve W (as returned by fastica) if you use FastICA.
>
>
ok. If you want to use FastICA in fastica you need to add 2 attributes to
FastICA:

unmixing_matrix_
whitener_

and clarify that components_ is the "linear operator to apply to the data
to get the independent sources. Made using unmixing_matrix_.dot(whitener_)
if whiten is set to True"

or something like that

>

Commit message: removed comments

Patch to review:
diff --git a/sklearn/decomposition/fastica_.py b/sklearn/decomposition/fastica_.py
index dffce0dc0d8bc..85682b84883f5 100644
--- a/sklearn/decomposition/fastica_.py
+++ b/sklearn/decomposition/fastica_.py
@@ -264,120 +264,37 @@ def my_g(x):
     pp. 411-430*
 
     """
-    random_state = check_random_state(random_state)
-    fun_args = {} if fun_args is None else fun_args
-    # make interface compatible with other decompositions
-    # a copy is required only for non whitened data
-    X = check_array(X, copy=whiten, dtype=FLOAT_DTYPES,
-                    ensure_min_samples=2).T
-
-    alpha = fun_args.get('alpha', 1.0)
-    if not 1 <= alpha <= 2:
-        raise ValueError('alpha must be in [1,2]')
-
-    if fun == 'logcosh':
-        g = _logcosh
-    elif fun == 'exp':
-        g = _exp
-    elif fun == 'cube':
-        g = _cube
-    elif callable(fun):
-        def g(x, fun_args):
-            return fun(x, **fun_args)
-    else:
-        exc = ValueError if isinstance(fun, str) else TypeError
-        raise exc("Unknown function %r;"
-                  " should be one of 'logcosh', 'exp', 'cube' or callable"
-                  % fun)
-
-    n, p = X.shape
 
-    if not whiten and n_components is not None:
-        n_components = None
-        warnings.warn('Ignoring n_components with whiten=False.')
-
-    if n_components is None:
-        n_components = min(n, p)
-    if (n_components > min(n, p)):
-        n_components = min(n, p)
-        warnings.warn('n_components is too large: it will be set to %s' % n_components)
+    est = FastICA(n_components=n_components, algorithm=algorithm,
+                  whiten=whiten, fun=fun, fun_args=fun_args,
+                  max_iter=max_iter, tol=tol, w_init=w_init,
+                  random_state=random_state)
+    sources = est._fit(X, compute_sources=compute_sources)
 
     if whiten:
-        # Centering the columns (ie the variables)
-        X_mean = X.mean(axis=-1)
-        X -= X_mean[:, np.newaxis]
-
-        # Whitening and preprocessing by PCA
-        u, d, _ = linalg.svd(X, full_matrices=False)
-
-        del _
-        K = (u / d).T[:n_components]  # see (6.33) p.140
-        del u, d
-        X1 = np.dot(K, X)
-        # see (13.6) p.267 Here X1 is white and data
-        # in X has been projected onto a subspace by PCA
-        X1 *= np.sqrt(p)
-    else:
-        # X must be casted to floats to avoid typing issues with numpy
-        # 2.0 and the line below
-        X1 = as_float_array(X, copy=False)  # copy has been taken care of
-
-    if w_init is None:
-        w_init = np.asarray(random_state.normal(size=(n_components,
-                            n_components)), dtype=X1.dtype)
-
-    else:
-        w_init = np.asarray(w_init)
-        if w_init.shape != (n_components, n_components):
-            raise ValueError('w_init has invalid shape -- should be %(shape)s'
-                             % {'shape': (n_components, n_components)})
-
-    kwargs = {'tol': tol,
-              'g': g,
-              'fun_args': fun_args,
-              'max_iter': max_iter,
-              'w_init': w_init}
-
-    if algorithm == 'parallel':
-        W, n_iter = _ica_par(X1, **kwargs)
-    elif algorithm == 'deflation':
-        W, n_iter = _ica_def(X1, **kwargs)
-    else:
-        raise ValueError('Invalid algorithm: must be either `parallel` or'
-                         ' `deflation`.')
-    del X1
-
-    if whiten:
-        if compute_sources:
-            S = np.dot(np.dot(W, K), X).T
-        else:
-            S = None
         if return_X_mean:
             if return_n_iter:
-                return K, W, S, X_mean, n_iter
+                return (est.whitening_, est._unmixing, sources, est.mean_,
+                        est.n_iter_)
             else:
-                return K, W, S, X_mean
+                return est.whitening_, est._unmixing, sources, est.mean_
         else:
             if return_n_iter:
-                return K, W, S, n_iter
+                return est.whitening_, est._unmixing, sources, est.n_iter_
             else:
-                return K, W, S
+                return est.whitening_, est._unmixing, sources
 
     else:
-        if compute_sources:
-            S = np.dot(W, X).T
-        else:
-            S = None
         if return_X_mean:
             if return_n_iter:
-                return None, W, S, None, n_iter
+                return None, sources, None, est.n_iter_
             else:
-                return None, W, S, None
+                return None, sources, None
         else:
             if return_n_iter:
-                return None, W, S, n_iter
+                return None, est._unmixing, sources, est.n_iter_
             else:
-                return None, W, S
+                return None, est._unmixing, sources
 
 
 class FastICA(TransformerMixin, BaseEstimator):
@@ -431,7 +348,10 @@ def my_g(x):
     Attributes
     ----------
     components_ : 2D array, shape (n_components, n_features)
-        The unmixing matrix.
+        The linear operator to apply to the data to get the independent
+        sources. This is equal to the unmixing matrix when ``whiten`` is
+        False, and equal to ``np.dot(unmixing_matrix, whitening_)`` when
+        ``whiten`` is True.
 
     mixing_ : array, shape (n_features, n_components)
         The mixing matrix.
@@ -502,26 +422,121 @@ def _fit(self, X, compute_sources=False):
             X_new : array-like, shape (n_samples, n_components)
         """
         fun_args = {} if self.fun_args is None else self.fun_args
-        whitening, unmixing, sources, X_mean, self.n_iter_ = fastica(
-            X=X, n_components=self.n_components, algorithm=self.algorithm,
-            whiten=self.whiten, fun=self.fun, fun_args=fun_args,
-            max_iter=self.max_iter, tol=self.tol, w_init=self.w_init,
-            random_state=self.random_state, return_X_mean=True,
-            compute_sources=compute_sources, return_n_iter=True)
+        random_state = check_random_state(self.random_state)
+
+        # make interface compatible with other decompositions
+        # a copy is required only for non whitened data
+        X = check_array(X, copy=self.whiten, dtype=FLOAT_DTYPES,
+                        ensure_min_samples=2).T
+
+        alpha = fun_args.get('alpha', 1.0)
+        if not 1 <= alpha <= 2:
+            raise ValueError('alpha must be in [1,2]')
+
+        if self.fun == 'logcosh':
+            g = _logcosh
+        elif self.fun == 'exp':
+            g = _exp
+        elif self.fun == 'cube':
+            g = _cube
+        elif callable(self.fun):
+            def g(x, fun_args):
+                return self.fun(x, **fun_args)
+        else:
+            exc = ValueError if isinstance(self.fun, str) else TypeError
+            raise exc(
+                "Unknown function %r;"
+                " should be one of 'logcosh', 'exp', 'cube' or callable"
+                % self.fun
+            )
+
+        n, p = X.shape
+
+        n_components = self.n_components
+        if not self.whiten and n_components is not None:
+            n_components = None
+            warnings.warn('Ignoring n_components with whiten=False.')
+
+        if n_components is None:
+            n_components = min(n, p)
+        if (n_components > min(n, p)):
+            n_components = min(n, p)
+            warnings.warn(
+                'n_components is too large: it will be set to %s'
+                % n_components
+            )
+
+        if self.whiten:
+            # Centering the columns (ie the variables)
+            X_mean = X.mean(axis=-1)
+            X -= X_mean[:, np.newaxis]
+
+            # Whitening and preprocessing by PCA
+            u, d, _ = linalg.svd(X, full_matrices=False)
+
+            del _
+            K = (u / d).T[:n_components]  # see (6.33) p.140
+            del u, d
+            X1 = np.dot(K, X)
+            # see (13.6) p.267 Here X1 is white and data
+            # in X has been projected onto a subspace by PCA
+            X1 *= np.sqrt(p)
+        else:
+            # X must be casted to floats to avoid typing issues with numpy
+            # 2.0 and the line below
+            X1 = as_float_array(X, copy=False)  # copy has been taken care of
+
+        w_init = self.w_init
+        if w_init is None:
+            w_init = np.asarray(random_state.normal(
+                size=(n_components, n_components)), dtype=X1.dtype)
+
+        else:
+            w_init = np.asarray(w_init)
+            if w_init.shape != (n_components, n_components):
+                raise ValueError(
+                    'w_init has invalid shape -- should be %(shape)s'
+                    % {'shape': (n_components, n_components)})
+
+        kwargs = {'tol': self.tol,
+                  'g': g,
+                  'fun_args': fun_args,
+                  'max_iter': self.max_iter,
+                  'w_init': w_init}
+
+        if self.algorithm == 'parallel':
+            W, n_iter = _ica_par(X1, **kwargs)
+        elif self.algorithm == 'deflation':
+            W, n_iter = _ica_def(X1, **kwargs)
+        else:
+            raise ValueError('Invalid algorithm: must be either `parallel` or'
+                             ' `deflation`.')
+        del X1
+
+        if compute_sources:
+            if self.whiten:
+                S = np.dot(np.dot(W, K), X).T
+            else:
+                S = np.dot(W, X).T
+        else:
+            S = None
+
+        self.n_iter_ = n_iter
 
         if self.whiten:
-            self.components_ = np.dot(unmixing, whitening)
+            self.components_ = np.dot(W, K)
             self.mean_ = X_mean
-            self.whitening_ = whitening
+            self.whitening_ = K
         else:
-            self.components_ = unmixing
+            self.components_ = W
 
         self.mixing_ = linalg.pinv(self.components_)
+        self._unmixing = W
 
         if compute_sources:
-            self.__sources = sources
+            self.__sources = S
 
-        return sources
+        return S
 
     def fit_transform(self, X, y=None):
         """Fit the model and recover the sources from X.


Code review comments:
diff_hunk: @@ -431,7 +348,10 @@ def my_g(x):
     Attributes
     ----------
     components_ : 2D array, shape (n_components, n_features)
-        The unmixing matrix.
+        The linear operator to apply to the data to get the independent
+        sources. This is equal to the unmixing matrix when ``whiten`` is
+        False, and equal to ``np.dot(unmixing_matrix, whitening_)`` when
path: sklearn/decomposition/fastica_.py
text: @user1:
```suggestion
        False, and equal to ``np.dot(unmixing_matrix_, whitening_)`` when
```

@author:
The unmixing matrix is a private attribute

diff_hunk: @@ -502,26 +422,121 @@ def _fit(self, X, compute_sources=False):
             X_new : array-like, shape (n_samples, n_components)
         """
         fun_args = {} if self.fun_args is None else self.fun_args
-        whitening, unmixing, sources, X_mean, self.n_iter_ = fastica(
-            X=X, n_components=self.n_components, algorithm=self.algorithm,
-            whiten=self.whiten, fun=self.fun, fun_args=fun_args,
-            max_iter=self.max_iter, tol=self.tol, w_init=self.w_init,
-            random_state=self.random_state, return_X_mean=True,
-            compute_sources=compute_sources, return_n_iter=True)
+        random_state = check_random_state(self.random_state)
+
+        # make interface compatible with other decompositions
+        # a copy is required only for non whitened data
+        X = check_array(X, copy=self.whiten, dtype=FLOAT_DTYPES,
+                        ensure_min_samples=2).T
+
+        alpha = fun_args.get('alpha', 1.0)
+        if not 1 <= alpha <= 2:
+            raise ValueError('alpha must be in [1,2]')
+
+        if self.fun == 'logcosh':
+            g = _logcosh
+        elif self.fun == 'exp':
+            g = _exp
+        elif self.fun == 'cube':
+            g = _cube
+        elif callable(self.fun):
+            def g(x, fun_args):
+                return self.fun(x, **fun_args)
+        else:
+            exc = ValueError if isinstance(self.fun, str) else TypeError
+            raise exc(
+                "Unknown function %r;"
+                " should be one of 'logcosh', 'exp', 'cube' or callable"
+                % self.fun
+            )
+
+        n, p = X.shape
path: sklearn/decomposition/fastica_.py
text: @user1:
a good occasion to change this to 

`n_samples, n_features = X.shape`

@author:
OK but I'd rather not change to code, this is out of scope for this PR