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