Aluode/PerceptionLabPortable
0
1"""Matrix factorization with Sparse PCA."""
2
3# Authors: The scikit-learn developers
4# SPDX-License-Identifier: BSD-3-Clause
5
6from numbers import Integral, Real
7
8import numpy as np
9
10from ..base import (
11 BaseEstimator,
12 ClassNamePrefixFeaturesOutMixin,
13 TransformerMixin,
14 _fit_context,
15)
16from ..linear_model import ridge_regression
17from ..utils import check_random_state
18from ..utils._param_validation import Interval, StrOptions
19from ..utils.extmath import svd_flip
20from ..utils.validation import check_array, check_is_fitted, validate_data
21from ._dict_learning import MiniBatchDictionaryLearning, dict_learning
22
23
24class _BaseSparsePCA(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator):
25 """Base class for SparsePCA and MiniBatchSparsePCA"""
26
27 _parameter_constraints: dict = {
28 "n_components": [None, Interval(Integral, 1, None, closed="left")],
29 "alpha": [Interval(Real, 0.0, None, closed="left")],
30 "ridge_alpha": [Interval(Real, 0.0, None, closed="left")],
31 "max_iter": [Interval(Integral, 0, None, closed="left")],
32 "tol": [Interval(Real, 0.0, None, closed="left")],
33 "method": [StrOptions({"lars", "cd"})],
34 "n_jobs": [Integral, None],
35 "verbose": ["verbose"],
36 "random_state": ["random_state"],
37 }
38
39 def __init__(
40 self,
41 n_components=None,
42 *,
43 alpha=1,
44 ridge_alpha=0.01,
45 max_iter=1000,
46 tol=1e-8,
47 method="lars",
48 n_jobs=None,
49 verbose=False,
50 random_state=None,
51 ):
52 self.n_components = n_components
53 self.alpha = alpha
54 self.ridge_alpha = ridge_alpha
55 self.max_iter = max_iter
56 self.tol = tol
57 self.method = method
58 self.n_jobs = n_jobs
59 self.verbose = verbose
60 self.random_state = random_state
61
62 @_fit_context(prefer_skip_nested_validation=True)
63 def fit(self, X, y=None):
64 """Fit the model from data in X.
65
66 Parameters
67 ----------
68 X : array-like of shape (n_samples, n_features)
69 Training vector, where `n_samples` is the number of samples
70 and `n_features` is the number of features.
71
72 y : Ignored
73 Not used, present here for API consistency by convention.
74
75 Returns
76 -------
77 self : object
78 Returns the instance itself.
79 """
80 random_state = check_random_state(self.random_state)
81 X = validate_data(self, X)
82
83 self.mean_ = X.mean(axis=0)
84 X = X - self.mean_
85
86 if self.n_components is None:
87 n_components = X.shape[1]
88 else:
89 n_components = self.n_components
90
91 return self._fit(X, n_components, random_state)
92
93 def transform(self, X):
94 """Least Squares projection of the data onto the sparse components.
95
96 To avoid instability issues in case the system is under-determined,
97 regularization can be applied (Ridge regression) via the
98 `ridge_alpha` parameter.
99
100 Note that Sparse PCA components orthogonality is not enforced as in PCA
101 hence one cannot use a simple linear projection.
102
103 Parameters
104 ----------
105 X : ndarray of shape (n_samples, n_features)
106 Test data to be transformed, must have the same number of
107 features as the data used to train the model.
108
109 Returns
110 -------
111 X_new : ndarray of shape (n_samples, n_components)
112 Transformed data.
113 """
114 check_is_fitted(self)
115
116 X = validate_data(self, X, reset=False)
117 X = X - self.mean_
118
119 U = ridge_regression(
120 self.components_.T, X.T, self.ridge_alpha, solver="cholesky"
121 )
122
123 return U
124
125 def inverse_transform(self, X):
126 """Transform data from the latent space to the original space.
127
128 This inversion is an approximation due to the loss of information
129 induced by the forward decomposition.
130
131 .. versionadded:: 1.2
132
133 Parameters
134 ----------
135 X : ndarray of shape (n_samples, n_components)
136 Data in the latent space.
137
138 Returns
139 -------
140 X_original : ndarray of shape (n_samples, n_features)
141 Reconstructed data in the original space.
142 """
143 check_is_fitted(self)
144 X = check_array(X)
145
146 return (X @ self.components_) + self.mean_
147
148 @property
149 def _n_features_out(self):
150 """Number of transformed output features."""
151 return self.components_.shape[0]
152
153 def __sklearn_tags__(self):
154 tags = super().__sklearn_tags__()
155 tags.transformer_tags.preserves_dtype = ["float64", "float32"]
156 return tags
157
158
159class SparsePCA(_BaseSparsePCA):
160 """Sparse Principal Components Analysis (SparsePCA).
161
162 Finds the set of sparse components that can optimally reconstruct
163 the data. The amount of sparseness is controllable by the coefficient
164 of the L1 penalty, given by the parameter alpha.
165
166 Read more in the :ref:`User Guide <SparsePCA>`.
167
168 Parameters
169 ----------
170 n_components : int, default=None
171 Number of sparse atoms to extract. If None, then ``n_components``
172 is set to ``n_features``.
173
174 alpha : float, default=1
175 Sparsity controlling parameter. Higher values lead to sparser
176 components.
177
178 ridge_alpha : float, default=0.01
179 Amount of ridge shrinkage to apply in order to improve
180 conditioning when calling the transform method.
181
182 max_iter : int, default=1000
183 Maximum number of iterations to perform.
184
185 tol : float, default=1e-8
186 Tolerance for the stopping condition.
187
188 method : {'lars', 'cd'}, default='lars'
189 Method to be used for optimization.
190 lars: uses the least angle regression method to solve the lasso problem
191 (linear_model.lars_path)
192 cd: uses the coordinate descent method to compute the
193 Lasso solution (linear_model.Lasso). Lars will be faster if
194 the estimated components are sparse.
195
196 n_jobs : int, default=None
197 Number of parallel jobs to run.
198 ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
199 ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
200 for more details.
201
202 U_init : ndarray of shape (n_samples, n_components), default=None
203 Initial values for the loadings for warm restart scenarios. Only used
204 if `U_init` and `V_init` are not None.
205
206 V_init : ndarray of shape (n_components, n_features), default=None
207 Initial values for the components for warm restart scenarios. Only used
208 if `U_init` and `V_init` are not None.
209
210 verbose : int or bool, default=False
211 Controls the verbosity; the higher, the more messages. Defaults to 0.
212
213 random_state : int, RandomState instance or None, default=None
214 Used during dictionary learning. Pass an int for reproducible results
215 across multiple function calls.
216 See :term:`Glossary <random_state>`.
217
218 Attributes
219 ----------
220 components_ : ndarray of shape (n_components, n_features)
221 Sparse components extracted from the data.
222
223 error_ : ndarray
224 Vector of errors at each iteration.
225
226 n_components_ : int
227 Estimated number of components.
228
229 .. versionadded:: 0.23
230
231 n_iter_ : int
232 Number of iterations run.
233
234 mean_ : ndarray of shape (n_features,)
235 Per-feature empirical mean, estimated from the training set.
236 Equal to ``X.mean(axis=0)``.
237
238 n_features_in_ : int
239 Number of features seen during :term:`fit`.
240
241 .. versionadded:: 0.24
242
243 feature_names_in_ : ndarray of shape (`n_features_in_`,)
244 Names of features seen during :term:`fit`. Defined only when `X`
245 has feature names that are all strings.
246
247 .. versionadded:: 1.0
248
249 See Also
250 --------
251 PCA : Principal Component Analysis implementation.
252 MiniBatchSparsePCA : Mini batch variant of `SparsePCA` that is faster but less
253 accurate.
254 DictionaryLearning : Generic dictionary learning problem using a sparse code.
255
256 Examples
257 --------
258 >>> import numpy as np
259 >>> from sklearn.datasets import make_friedman1
260 >>> from sklearn.decomposition import SparsePCA
261 >>> X, _ = make_friedman1(n_samples=200, n_features=30, random_state=0)
262 >>> transformer = SparsePCA(n_components=5, random_state=0)
263 >>> transformer.fit(X)
264 SparsePCA(...)
265 >>> X_transformed = transformer.transform(X)
266 >>> X_transformed.shape
267 (200, 5)
268 >>> # most values in the components_ are zero (sparsity)
269 >>> np.mean(transformer.components_ == 0)
270 np.float64(0.9666)
271 """
272
273 _parameter_constraints: dict = {
274 **_BaseSparsePCA._parameter_constraints,
275 "U_init": [None, np.ndarray],
276 "V_init": [None, np.ndarray],
277 }
278
279 def __init__(
280 self,
281 n_components=None,
282 *,
283 alpha=1,
284 ridge_alpha=0.01,
285 max_iter=1000,
286 tol=1e-8,
287 method="lars",
288 n_jobs=None,
289 U_init=None,
290 V_init=None,
291 verbose=False,
292 random_state=None,
293 ):
294 super().__init__(
295 n_components=n_components,
296 alpha=alpha,
297 ridge_alpha=ridge_alpha,
298 max_iter=max_iter,
299 tol=tol,
300 method=method,
301 n_jobs=n_jobs,
302 verbose=verbose,
303 random_state=random_state,
304 )
305 self.U_init = U_init
306 self.V_init = V_init
307
308 def _fit(self, X, n_components, random_state):
309 """Specialized `fit` for SparsePCA."""
310
311 code_init = self.V_init.T if self.V_init is not None else None
312 dict_init = self.U_init.T if self.U_init is not None else None
313 code, dictionary, E, self.n_iter_ = dict_learning(
314 X.T,
315 n_components,
316 alpha=self.alpha,
317 tol=self.tol,
318 max_iter=self.max_iter,
319 method=self.method,
320 n_jobs=self.n_jobs,
321 verbose=self.verbose,
322 random_state=random_state,
323 code_init=code_init,
324 dict_init=dict_init,
325 return_n_iter=True,
326 )
327 # flip eigenvectors' sign to enforce deterministic output
328 code, dictionary = svd_flip(code, dictionary, u_based_decision=True)
329 self.components_ = code.T
330 components_norm = np.linalg.norm(self.components_, axis=1)[:, np.newaxis]
331 components_norm[components_norm == 0] = 1
332 self.components_ /= components_norm
333 self.n_components_ = len(self.components_)
334
335 self.error_ = E
336 return self
337
338
339class MiniBatchSparsePCA(_BaseSparsePCA):
340 """Mini-batch Sparse Principal Components Analysis.
341
342 Finds the set of sparse components that can optimally reconstruct
343 the data. The amount of sparseness is controllable by the coefficient
344 of the L1 penalty, given by the parameter alpha.
345
346 For an example comparing sparse PCA to PCA, see
347 :ref:`sphx_glr_auto_examples_decomposition_plot_faces_decomposition.py`
348
349 Read more in the :ref:`User Guide <SparsePCA>`.
350
351 Parameters
352 ----------
353 n_components : int, default=None
354 Number of sparse atoms to extract. If None, then ``n_components``
355 is set to ``n_features``.
356
357 alpha : int, default=1
358 Sparsity controlling parameter. Higher values lead to sparser
359 components.
360
361 ridge_alpha : float, default=0.01
362 Amount of ridge shrinkage to apply in order to improve
363 conditioning when calling the transform method.
364
365 max_iter : int, default=1_000
366 Maximum number of iterations over the complete dataset before
367 stopping independently of any early stopping criterion heuristics.
368
369 .. versionadded:: 1.2
370
371 callback : callable, default=None
372 Callable that gets invoked every five iterations.
373
374 batch_size : int, default=3
375 The number of features to take in each mini batch.
376
377 verbose : int or bool, default=False
378 Controls the verbosity; the higher, the more messages. Defaults to 0.
379
380 shuffle : bool, default=True
381 Whether to shuffle the data before splitting it in batches.
382
383 n_jobs : int, default=None
384 Number of parallel jobs to run.
385 ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
386 ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
387 for more details.
388
389 method : {'lars', 'cd'}, default='lars'
390 Method to be used for optimization.
391 lars: uses the least angle regression method to solve the lasso problem
392 (linear_model.lars_path)
393 cd: uses the coordinate descent method to compute the
394 Lasso solution (linear_model.Lasso). Lars will be faster if
395 the estimated components are sparse.
396
397 random_state : int, RandomState instance or None, default=None
398 Used for random shuffling when ``shuffle`` is set to ``True``,
399 during online dictionary learning. Pass an int for reproducible results
400 across multiple function calls.
401 See :term:`Glossary <random_state>`.
402
403 tol : float, default=1e-3
404 Control early stopping based on the norm of the differences in the
405 dictionary between 2 steps.
406
407 To disable early stopping based on changes in the dictionary, set
408 `tol` to 0.0.
409
410 .. versionadded:: 1.1
411
412 max_no_improvement : int or None, default=10
413 Control early stopping based on the consecutive number of mini batches
414 that does not yield an improvement on the smoothed cost function.
415
416 To disable convergence detection based on cost function, set
417 `max_no_improvement` to `None`.
418
419 .. versionadded:: 1.1
420
421 Attributes
422 ----------
423 components_ : ndarray of shape (n_components, n_features)
424 Sparse components extracted from the data.
425
426 n_components_ : int
427 Estimated number of components.
428
429 .. versionadded:: 0.23
430
431 n_iter_ : int
432 Number of iterations run.
433
434 mean_ : ndarray of shape (n_features,)
435 Per-feature empirical mean, estimated from the training set.
436 Equal to ``X.mean(axis=0)``.
437
438 n_features_in_ : int
439 Number of features seen during :term:`fit`.
440
441 .. versionadded:: 0.24
442
443 feature_names_in_ : ndarray of shape (`n_features_in_`,)
444 Names of features seen during :term:`fit`. Defined only when `X`
445 has feature names that are all strings.
446
447 .. versionadded:: 1.0
448
449 See Also
450 --------
451 DictionaryLearning : Find a dictionary that sparsely encodes data.
452 IncrementalPCA : Incremental principal components analysis.
453 PCA : Principal component analysis.
454 SparsePCA : Sparse Principal Components Analysis.
455 TruncatedSVD : Dimensionality reduction using truncated SVD.
456
457 Examples
458 --------
459 >>> import numpy as np
460 >>> from sklearn.datasets import make_friedman1
461 >>> from sklearn.decomposition import MiniBatchSparsePCA
462 >>> X, _ = make_friedman1(n_samples=200, n_features=30, random_state=0)
463 >>> transformer = MiniBatchSparsePCA(n_components=5, batch_size=50,
464 ... max_iter=10, random_state=0)
465 >>> transformer.fit(X)
466 MiniBatchSparsePCA(...)
467 >>> X_transformed = transformer.transform(X)
468 >>> X_transformed.shape
469 (200, 5)
470 >>> # most values in the components_ are zero (sparsity)
471 >>> np.mean(transformer.components_ == 0)
472 np.float64(0.9)
473 """
474
475 _parameter_constraints: dict = {
476 **_BaseSparsePCA._parameter_constraints,
477 "max_iter": [Interval(Integral, 0, None, closed="left")],
478 "callback": [None, callable],
479 "batch_size": [Interval(Integral, 1, None, closed="left")],
480 "shuffle": ["boolean"],
481 "max_no_improvement": [Interval(Integral, 0, None, closed="left"), None],
482 }
483
484 def __init__(
485 self,
486 n_components=None,
487 *,
488 alpha=1,
489 ridge_alpha=0.01,
490 max_iter=1_000,
491 callback=None,
492 batch_size=3,
493 verbose=False,
494 shuffle=True,
495 n_jobs=None,
496 method="lars",
497 random_state=None,
498 tol=1e-3,
499 max_no_improvement=10,
500 ):
501 super().__init__(
502 n_components=n_components,
503 alpha=alpha,
504 ridge_alpha=ridge_alpha,
505 max_iter=max_iter,
506 tol=tol,
507 method=method,
508 n_jobs=n_jobs,
509 verbose=verbose,
510 random_state=random_state,
511 )
512 self.callback = callback
513 self.batch_size = batch_size
514 self.shuffle = shuffle
515 self.max_no_improvement = max_no_improvement
516
517 def _fit(self, X, n_components, random_state):
518 """Specialized `fit` for MiniBatchSparsePCA."""
519
520 transform_algorithm = "lasso_" + self.method
521 est = MiniBatchDictionaryLearning(
522 n_components=n_components,
523 alpha=self.alpha,
524 max_iter=self.max_iter,
525 dict_init=None,
526 batch_size=self.batch_size,
527 shuffle=self.shuffle,
528 n_jobs=self.n_jobs,
529 fit_algorithm=self.method,
530 random_state=random_state,
531 transform_algorithm=transform_algorithm,
532 transform_alpha=self.alpha,
533 verbose=self.verbose,
534 callback=self.callback,
535 tol=self.tol,
536 max_no_improvement=self.max_no_improvement,
537 )
538 est.set_output(transform="default")
539 est.fit(X.T)
540
541 self.components_, self.n_iter_ = est.transform(X.T).T, est.n_iter_
542
543 components_norm = np.linalg.norm(self.components_, axis=1)[:, np.newaxis]
544 components_norm[components_norm == 0] = 1
545 self.components_ /= components_norm
546 self.n_components_ = len(self.components_)
547
548 return self
549