1# Authors: Pierre Ablin <pierre.ablin@inria.fr>
2#
3# License: BSD (3-clause)
4import warnings
5
6import numpy as np
7
8from scipy import linalg
9
10from sklearn.decomposition import FastICA
11from sklearn.utils import check_random_state, as_float_array
12from sklearn.utils.validation import FLOAT_DTYPES
13
14from .solver import picard
15
16
17class Picard(FastICA):
18    """Picard: a **very** fast algorithm for Independent Component Analysis.
19
20    Parameters
21    ----------
22    n_components : int, default=None
23        Number of components to use. If None is passed, all are used.
24    ortho : bool, optional
25        If True, uses Picard-O and enforce an orthogonal constraint.
26        Otherwise, uses the standard Picard.
27    extended : None or bool, optional
28        If True, uses the extended algorithm to separate sub and super-gaussian
29        sources. If ``None`` (default), it is set to True if `ortho == True`,
30        and `False` otherwise. With `extended=True` we recommend you keep the
31        different density to `'tanh'`. See notes below.
32    whiten : bool, default=True
33        If whiten is false, the data is already considered to be
34        whitened, and no whitening is performed.
35    fun : str or class, optional
36        Either a built-in density model ('tanh', 'exp' and 'cube'), or a custom
37        density.
38        A custom density is a class that should contain two methods called
39        'log_lik' and 'score_and_der'. See examples in the densities.py file.
40    max_iter : int, default=200
41        Maximum number of iterations during fit.
42    tol : float, default=1e-7
43        Tolerance on update at each iteration.
44    w_init : ndarray of shape (n_components, n_components), default=None
45        The mixing matrix to be used to initialize the algorithm.
46    m : int, optional
47        Size of L-BFGS's memory.
48    ls_tries : int, optional
49        Number of attempts during the backtracking line-search.
50    lambda_min : float, optional
51        Threshold on the eigenvalues of the Hessian approximation. Any
52        eigenvalue below lambda_min is shifted to lambda_min.
53    random_state : int, RandomState instance or None, default=None
54        Used to initialize ``w_init`` when not specified, with a
55        normal distribution. Pass an int, for reproducible results
56        across multiple function calls.
57
58    Attributes
59    ----------
60    components_ : ndarray of shape (n_components, n_features)
61        The linear operator to apply to the data to get the independent
62        sources. This is equal to the unmixing matrix when ``whiten`` is
63        False, and equal to ``np.dot(unmixing_matrix, self.whitening_)`` when
64        ``whiten`` is True.
65    mixing_ : ndarray of shape (n_features, n_components)
66        The pseudo-inverse of ``components_``. It is the linear operator
67        that maps independent sources to the data.
68    mean_ : ndarray of shape(n_features,)
69        The mean over features. Only set if `self.whiten` is True.
70    whitening_ : ndarray of shape (n_components, n_features)
71        Only set if whiten is 'True'. This is the pre-whitening matrix
72        that projects data onto the first `n_components` principal components.
73
74    Examples
75    --------
76    >>> from sklearn.datasets import load_digits
77    >>> from picard import Picard
78    >>> X, _ = load_digits(return_X_y=True)
79    >>> transformer = Picard(n_components=7,
80    ...         random_state=0)
81    >>> X_transformed = transformer.fit_transform(X)
82    >>> X_transformed.shape
83    (1797, 7)
84
85    Notes
86    -----
87    Using a different density than `'tanh'` may lead to erratic behavior of
88    the algorithm: when `extended=True`, the non-linearity used by the
89    algorithm is `x +/- fun(x)`. The non-linearity should correspond to a
90    density, hence `fun` should be dominated by `x ** 2`. Further,
91    `x + fun(x)` should separate super-Gaussian sources and `x-fun(x)`
92    should separate sub-Gaussian sources. This set of requirement is met by
93    `'tanh'`.
94    """
95    def __init__(self, n_components=None, *, ortho=True, extended=None,
96                 whiten=True, fun='tanh', max_iter=200, tol=1e-7,
97                 w_init=None, m=7,  ls_tries=10, lambda_min=0.01,
98                 random_state=None):
99        super().__init__()
100        if max_iter < 1:
101            raise ValueError("max_iter should be greater than 1, got "
102                             "(max_iter={})".format(max_iter))
103        self.n_components = n_components
104        self.ortho = ortho
105        self.extended = extended
106        self.whiten = whiten
107        self.fun = fun
108        self.max_iter = max_iter
109        self.tol = tol
110        self.w_init = w_init
111        self.m = m
112        self.ls_tries = ls_tries
113        self.lambda_min = lambda_min
114        self.random_state = random_state
115
116    def _fit(self, X, compute_sources=False):
117        """Fit the model
118
119        Parameters
120        ----------
121        X : array-like of shape (n_samples, n_features)
122            Training data, where n_samples is the number of samples
123            and n_features is the number of features.
124        compute_sources : bool, default=False
125            If False, sources are not computes but only the rotation matrix.
126            This can save memory when working with big data. Defaults to False.
127
128        Returns
129        -------
130        X_new : ndarray of shape (n_samples, n_components)
131            The estimated sources.
132        """
133
134        X = self._validate_data(X, copy=self.whiten, dtype=FLOAT_DTYPES,
135                                ensure_min_samples=2).T
136        random_state = check_random_state(self.random_state)
137
138        n_samples, n_features = X.shape
139
140        n_components = self.n_components
141        if not self.whiten and n_components is not None:
142            n_components = None
143            warnings.warn('Ignoring n_components with whiten=False.')
144
145        if n_components is None:
146            n_components = min(n_samples, n_features)
147        if (n_components > min(n_samples, n_features)):
148            n_components = min(n_samples, n_features)
149            warnings.warn(
150                'n_components is too large: it will be set to %s'
151                % n_components
152            )
153
154        if self.whiten:
155            # Centering the columns (ie the variables)
156            X_mean = X.mean(axis=-1)
157            X -= X_mean[:, np.newaxis]
158
159            # Whitening and preprocessing by PCA
160            u, d, _ = linalg.svd(X, full_matrices=False, check_finite=False)
161
162            del _
163            K = (u / d).T[:n_components]
164            del u, d
165            X1 = np.dot(K, X)
166            X1 *= np.sqrt(n_features)
167        else:
168            # X must be casted to floats to avoid typing issues with numpy
169            # 2.0 and the line below
170            X1 = as_float_array(X, copy=False)  # copy has been taken care of
171
172        w_init = self.w_init
173
174        kwargs = {'tol': self.tol,
175                  'fun': self.fun,
176                  'ortho': self.ortho,
177                  'extended': self.extended,
178                  'max_iter': self.max_iter,
179                  'w_init': w_init,
180                  'whiten': None,  # Already done above
181                  'm': self.m,
182                  'ls_tries': self.ls_tries,
183                  'lambda_min': self.lambda_min,
184                  'random_state': random_state}
185        _, W, _ = picard(X1, **kwargs)
186        del X1
187
188        if compute_sources:
189            if self.whiten:
190                S = np.linalg.multi_dot([W, K, X]).T
191            else:
192                S = np.dot(W, X).T
193        else:
194            S = None
195
196        if self.whiten:
197            self.components_ = np.dot(W, K)
198            self.mean_ = X_mean
199            self.whitening_ = K
200        else:
201            self.components_ = W
202
203        self.mixing_ = linalg.pinv(self.components_, check_finite=False)
204        self._unmixing = W
205
206        return S
207