Aluode/PerceptionLabPortable
0
1"""Principal Component Analysis."""
2
3# Authors: The scikit-learn developers
4# SPDX-License-Identifier: BSD-3-Clause
5
6from math import lgamma, log, sqrt
7from numbers import Integral, Real
8
9import numpy as np
10from scipy import linalg
11from scipy.sparse import issparse
12from scipy.sparse.linalg import svds
13
14from ..base import _fit_context
15from ..utils import check_random_state
16from ..utils._arpack import _init_arpack_v0
17from ..utils._array_api import _convert_to_numpy, get_namespace
18from ..utils._param_validation import Interval, RealNotInt, StrOptions
19from ..utils.extmath import _randomized_svd, fast_logdet, stable_cumsum, svd_flip
20from ..utils.sparsefuncs import _implicit_column_offset, mean_variance_axis
21from ..utils.validation import check_is_fitted, validate_data
22from ._base import _BasePCA
23
24
25def _assess_dimension(spectrum, rank, n_samples):
26 """Compute the log-likelihood of a rank ``rank`` dataset.
27
28 The dataset is assumed to be embedded in gaussian noise of shape(n,
29 dimf) having spectrum ``spectrum``. This implements the method of
30 T. P. Minka.
31
32 Parameters
33 ----------
34 spectrum : ndarray of shape (n_features,)
35 Data spectrum.
36 rank : int
37 Tested rank value. It should be strictly lower than n_features,
38 otherwise the method isn't specified (division by zero in equation
39 (31) from the paper).
40 n_samples : int
41 Number of samples.
42
43 Returns
44 -------
45 ll : float
46 The log-likelihood.
47
48 References
49 ----------
50 This implements the method of `Thomas P. Minka:
51 Automatic Choice of Dimensionality for PCA. NIPS 2000: 598-604
52 <https://proceedings.neurips.cc/paper/2000/file/7503cfacd12053d309b6bed5c89de212-Paper.pdf>`_
53 """
54 xp, _ = get_namespace(spectrum)
55
56 n_features = spectrum.shape[0]
57 if not 1 <= rank < n_features:
58 raise ValueError("the tested rank should be in [1, n_features - 1]")
59
60 eps = 1e-15
61
62 if spectrum[rank - 1] < eps:
63 # When the tested rank is associated with a small eigenvalue, there's
64 # no point in computing the log-likelihood: it's going to be very
65 # small and won't be the max anyway. Also, it can lead to numerical
66 # issues below when computing pa, in particular in log((spectrum[i] -
67 # spectrum[j]) because this will take the log of something very small.
68 return -xp.inf
69
70 pu = -rank * log(2.0)
71 for i in range(1, rank + 1):
72 pu += (
73 lgamma((n_features - i + 1) / 2.0) - log(xp.pi) * (n_features - i + 1) / 2.0
74 )
75
76 pl = xp.sum(xp.log(spectrum[:rank]))
77 pl = -pl * n_samples / 2.0
78
79 v = max(eps, xp.sum(spectrum[rank:]) / (n_features - rank))
80 pv = -log(v) * n_samples * (n_features - rank) / 2.0
81
82 m = n_features * rank - rank * (rank + 1.0) / 2.0
83 pp = log(2.0 * xp.pi) * (m + rank) / 2.0
84
85 pa = 0.0
86 spectrum_ = xp.asarray(spectrum, copy=True)
87 spectrum_[rank:n_features] = v
88 for i in range(rank):
89 for j in range(i + 1, spectrum.shape[0]):
90 pa += log(
91 (spectrum[i] - spectrum[j]) * (1.0 / spectrum_[j] - 1.0 / spectrum_[i])
92 ) + log(n_samples)
93
94 ll = pu + pl + pv + pp - pa / 2.0 - rank * log(n_samples) / 2.0
95
96 return ll
97
98
99def _infer_dimension(spectrum, n_samples):
100 """Infers the dimension of a dataset with a given spectrum.
101
102 The returned value will be in [1, n_features - 1].
103 """
104 xp, _ = get_namespace(spectrum)
105
106 ll = xp.empty_like(spectrum)
107 ll[0] = -xp.inf # we don't want to return n_components = 0
108 for rank in range(1, spectrum.shape[0]):
109 ll[rank] = _assess_dimension(spectrum, rank, n_samples)
110 return xp.argmax(ll)
111
112
113class PCA(_BasePCA):
114 """Principal component analysis (PCA).
115
116 Linear dimensionality reduction using Singular Value Decomposition of the
117 data to project it to a lower dimensional space. The input data is centered
118 but not scaled for each feature before applying the SVD.
119
120 It uses the LAPACK implementation of the full SVD or a randomized truncated
121 SVD by the method of Halko et al. 2009, depending on the shape of the input
122 data and the number of components to extract.
123
124 With sparse inputs, the ARPACK implementation of the truncated SVD can be
125 used (i.e. through :func:`scipy.sparse.linalg.svds`). Alternatively, one
126 may consider :class:`TruncatedSVD` where the data are not centered.
127
128 Notice that this class only supports sparse inputs for some solvers such as
129 "arpack" and "covariance_eigh". See :class:`TruncatedSVD` for an
130 alternative with sparse data.
131
132 For a usage example, see
133 :ref:`sphx_glr_auto_examples_decomposition_plot_pca_iris.py`
134
135 Read more in the :ref:`User Guide <PCA>`.
136
137 Parameters
138 ----------
139 n_components : int, float or 'mle', default=None
140 Number of components to keep.
141 if n_components is not set all components are kept::
142
143 n_components == min(n_samples, n_features)
144
145 If ``n_components == 'mle'`` and ``svd_solver == 'full'``, Minka's
146 MLE is used to guess the dimension. Use of ``n_components == 'mle'``
147 will interpret ``svd_solver == 'auto'`` as ``svd_solver == 'full'``.
148
149 If ``0 < n_components < 1`` and ``svd_solver == 'full'``, select the
150 number of components such that the amount of variance that needs to be
151 explained is greater than the percentage specified by n_components.
152
153 If ``svd_solver == 'arpack'``, the number of components must be
154 strictly less than the minimum of n_features and n_samples.
155
156 Hence, the None case results in::
157
158 n_components == min(n_samples, n_features) - 1
159
160 copy : bool, default=True
161 If False, data passed to fit are overwritten and running
162 fit(X).transform(X) will not yield the expected results,
163 use fit_transform(X) instead.
164
165 whiten : bool, default=False
166 When True (False by default) the `components_` vectors are multiplied
167 by the square root of n_samples and then divided by the singular values
168 to ensure uncorrelated outputs with unit component-wise variances.
169
170 Whitening will remove some information from the transformed signal
171 (the relative variance scales of the components) but can sometime
172 improve the predictive accuracy of the downstream estimators by
173 making their data respect some hard-wired assumptions.
174
175 svd_solver : {'auto', 'full', 'covariance_eigh', 'arpack', 'randomized'},\
176 default='auto'
177 "auto" :
178 The solver is selected by a default 'auto' policy is based on `X.shape` and
179 `n_components`: if the input data has fewer than 1000 features and
180 more than 10 times as many samples, then the "covariance_eigh"
181 solver is used. Otherwise, if the input data is larger than 500x500
182 and the number of components to extract is lower than 80% of the
183 smallest dimension of the data, then the more efficient
184 "randomized" method is selected. Otherwise the exact "full" SVD is
185 computed and optionally truncated afterwards.
186 "full" :
187 Run exact full SVD calling the standard LAPACK solver via
188 `scipy.linalg.svd` and select the components by postprocessing
189 "covariance_eigh" :
190 Precompute the covariance matrix (on centered data), run a
191 classical eigenvalue decomposition on the covariance matrix
192 typically using LAPACK and select the components by postprocessing.
193 This solver is very efficient for n_samples >> n_features and small
194 n_features. It is, however, not tractable otherwise for large
195 n_features (large memory footprint required to materialize the
196 covariance matrix). Also note that compared to the "full" solver,
197 this solver effectively doubles the condition number and is
198 therefore less numerical stable (e.g. on input data with a large
199 range of singular values).
200 "arpack" :
201 Run SVD truncated to `n_components` calling ARPACK solver via
202 `scipy.sparse.linalg.svds`. It requires strictly
203 `0 < n_components < min(X.shape)`
204 "randomized" :
205 Run randomized SVD by the method of Halko et al.
206
207 .. versionadded:: 0.18.0
208
209 .. versionchanged:: 1.5
210 Added the 'covariance_eigh' solver.
211
212 tol : float, default=0.0
213 Tolerance for singular values computed by svd_solver == 'arpack'.
214 Must be of range [0.0, infinity).
215
216 .. versionadded:: 0.18.0
217
218 iterated_power : int or 'auto', default='auto'
219 Number of iterations for the power method computed by
220 svd_solver == 'randomized'.
221 Must be of range [0, infinity).
222
223 .. versionadded:: 0.18.0
224
225 n_oversamples : int, default=10
226 This parameter is only relevant when `svd_solver="randomized"`.
227 It corresponds to the additional number of random vectors to sample the
228 range of `X` so as to ensure proper conditioning. See
229 :func:`~sklearn.utils.extmath.randomized_svd` for more details.
230
231 .. versionadded:: 1.1
232
233 power_iteration_normalizer : {'auto', 'QR', 'LU', 'none'}, default='auto'
234 Power iteration normalizer for randomized SVD solver.
235 Not used by ARPACK. See :func:`~sklearn.utils.extmath.randomized_svd`
236 for more details.
237
238 .. versionadded:: 1.1
239
240 random_state : int, RandomState instance or None, default=None
241 Used when the 'arpack' or 'randomized' solvers are used. Pass an int
242 for reproducible results across multiple function calls.
243 See :term:`Glossary <random_state>`.
244
245 .. versionadded:: 0.18.0
246
247 Attributes
248 ----------
249 components_ : ndarray of shape (n_components, n_features)
250 Principal axes in feature space, representing the directions of
251 maximum variance in the data. Equivalently, the right singular
252 vectors of the centered input data, parallel to its eigenvectors.
253 The components are sorted by decreasing ``explained_variance_``.
254
255 explained_variance_ : ndarray of shape (n_components,)
256 The amount of variance explained by each of the selected components.
257 The variance estimation uses `n_samples - 1` degrees of freedom.
258
259 Equal to n_components largest eigenvalues
260 of the covariance matrix of X.
261
262 .. versionadded:: 0.18
263
264 explained_variance_ratio_ : ndarray of shape (n_components,)
265 Percentage of variance explained by each of the selected components.
266
267 If ``n_components`` is not set then all components are stored and the
268 sum of the ratios is equal to 1.0.
269
270 singular_values_ : ndarray of shape (n_components,)
271 The singular values corresponding to each of the selected components.
272 The singular values are equal to the 2-norms of the ``n_components``
273 variables in the lower-dimensional space.
274
275 .. versionadded:: 0.19
276
277 mean_ : ndarray of shape (n_features,)
278 Per-feature empirical mean, estimated from the training set.
279
280 Equal to `X.mean(axis=0)`.
281
282 n_components_ : int
283 The estimated number of components. When n_components is set
284 to 'mle' or a number between 0 and 1 (with svd_solver == 'full') this
285 number is estimated from input data. Otherwise it equals the parameter
286 n_components, or the lesser value of n_features and n_samples
287 if n_components is None.
288
289 n_samples_ : int
290 Number of samples in the training data.
291
292 noise_variance_ : float
293 The estimated noise covariance following the Probabilistic PCA model
294 from Tipping and Bishop 1999. See "Pattern Recognition and
295 Machine Learning" by C. Bishop, 12.2.1 p. 574 or
296 http://www.miketipping.com/papers/met-mppca.pdf. It is required to
297 compute the estimated data covariance and score samples.
298
299 Equal to the average of (min(n_features, n_samples) - n_components)
300 smallest eigenvalues of the covariance matrix of X.
301
302 n_features_in_ : int
303 Number of features seen during :term:`fit`.
304
305 .. versionadded:: 0.24
306
307 feature_names_in_ : ndarray of shape (`n_features_in_`,)
308 Names of features seen during :term:`fit`. Defined only when `X`
309 has feature names that are all strings.
310
311 .. versionadded:: 1.0
312
313 See Also
314 --------
315 KernelPCA : Kernel Principal Component Analysis.
316 SparsePCA : Sparse Principal Component Analysis.
317 TruncatedSVD : Dimensionality reduction using truncated SVD.
318 IncrementalPCA : Incremental Principal Component Analysis.
319
320 References
321 ----------
322 For n_components == 'mle', this class uses the method from:
323 `Minka, T. P.. "Automatic choice of dimensionality for PCA".
324 In NIPS, pp. 598-604 <https://tminka.github.io/papers/pca/minka-pca.pdf>`_
325
326 Implements the probabilistic PCA model from:
327 `Tipping, M. E., and Bishop, C. M. (1999). "Probabilistic principal
328 component analysis". Journal of the Royal Statistical Society:
329 Series B (Statistical Methodology), 61(3), 611-622.
330 <http://www.miketipping.com/papers/met-mppca.pdf>`_
331 via the score and score_samples methods.
332
333 For svd_solver == 'arpack', refer to `scipy.sparse.linalg.svds`.
334
335 For svd_solver == 'randomized', see:
336 :doi:`Halko, N., Martinsson, P. G., and Tropp, J. A. (2011).
337 "Finding structure with randomness: Probabilistic algorithms for
338 constructing approximate matrix decompositions".
339 SIAM review, 53(2), 217-288.
340 <10.1137/090771806>`
341 and also
342 :doi:`Martinsson, P. G., Rokhlin, V., and Tygert, M. (2011).
343 "A randomized algorithm for the decomposition of matrices".
344 Applied and Computational Harmonic Analysis, 30(1), 47-68.
345 <10.1016/j.acha.2010.02.003>`
346
347 Examples
348 --------
349 >>> import numpy as np
350 >>> from sklearn.decomposition import PCA
351 >>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
352 >>> pca = PCA(n_components=2)
353 >>> pca.fit(X)
354 PCA(n_components=2)
355 >>> print(pca.explained_variance_ratio_)
356 [0.9924 0.0075]
357 >>> print(pca.singular_values_)
358 [6.30061 0.54980]
359
360 >>> pca = PCA(n_components=2, svd_solver='full')
361 >>> pca.fit(X)
362 PCA(n_components=2, svd_solver='full')
363 >>> print(pca.explained_variance_ratio_)
364 [0.9924 0.00755]
365 >>> print(pca.singular_values_)
366 [6.30061 0.54980]
367
368 >>> pca = PCA(n_components=1, svd_solver='arpack')
369 >>> pca.fit(X)
370 PCA(n_components=1, svd_solver='arpack')
371 >>> print(pca.explained_variance_ratio_)
372 [0.99244]
373 >>> print(pca.singular_values_)
374 [6.30061]
375 """
376
377 _parameter_constraints: dict = {
378 "n_components": [
379 Interval(Integral, 0, None, closed="left"),
380 Interval(RealNotInt, 0, 1, closed="neither"),
381 StrOptions({"mle"}),
382 None,
383 ],
384 "copy": ["boolean"],
385 "whiten": ["boolean"],
386 "svd_solver": [
387 StrOptions({"auto", "full", "covariance_eigh", "arpack", "randomized"})
388 ],
389 "tol": [Interval(Real, 0, None, closed="left")],
390 "iterated_power": [
391 StrOptions({"auto"}),
392 Interval(Integral, 0, None, closed="left"),
393 ],
394 "n_oversamples": [Interval(Integral, 1, None, closed="left")],
395 "power_iteration_normalizer": [StrOptions({"auto", "QR", "LU", "none"})],
396 "random_state": ["random_state"],
397 }
398
399 def __init__(
400 self,
401 n_components=None,
402 *,
403 copy=True,
404 whiten=False,
405 svd_solver="auto",
406 tol=0.0,
407 iterated_power="auto",
408 n_oversamples=10,
409 power_iteration_normalizer="auto",
410 random_state=None,
411 ):
412 self.n_components = n_components
413 self.copy = copy
414 self.whiten = whiten
415 self.svd_solver = svd_solver
416 self.tol = tol
417 self.iterated_power = iterated_power
418 self.n_oversamples = n_oversamples
419 self.power_iteration_normalizer = power_iteration_normalizer
420 self.random_state = random_state
421
422 @_fit_context(prefer_skip_nested_validation=True)
423 def fit(self, X, y=None):
424 """Fit the model with X.
425
426 Parameters
427 ----------
428 X : {array-like, sparse matrix} of shape (n_samples, n_features)
429 Training data, where `n_samples` is the number of samples
430 and `n_features` is the number of features.
431
432 y : Ignored
433 Ignored.
434
435 Returns
436 -------
437 self : object
438 Returns the instance itself.
439 """
440 self._fit(X)
441 return self
442
443 @_fit_context(prefer_skip_nested_validation=True)
444 def fit_transform(self, X, y=None):
445 """Fit the model with X and apply the dimensionality reduction on X.
446
447 Parameters
448 ----------
449 X : {array-like, sparse matrix} of shape (n_samples, n_features)
450 Training data, where `n_samples` is the number of samples
451 and `n_features` is the number of features.
452
453 y : Ignored
454 Ignored.
455
456 Returns
457 -------
458 X_new : ndarray of shape (n_samples, n_components)
459 Transformed values.
460
461 Notes
462 -----
463 This method returns a Fortran-ordered array. To convert it to a
464 C-ordered array, use 'np.ascontiguousarray'.
465 """
466 U, S, _, X, x_is_centered, xp = self._fit(X)
467 if U is not None:
468 U = U[:, : self.n_components_]
469
470 if self.whiten:
471 # X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples)
472 U *= sqrt(X.shape[0] - 1)
473 else:
474 # X_new = X * V = U * S * Vt * V = U * S
475 U *= S[: self.n_components_]
476
477 return U
478 else: # solver="covariance_eigh" does not compute U at fit time.
479 return self._transform(X, xp, x_is_centered=x_is_centered)
480
481 def _fit(self, X):
482 """Dispatch to the right submethod depending on the chosen solver."""
483 xp, is_array_api_compliant = get_namespace(X)
484
485 # Raise an error for sparse input and unsupported svd_solver
486 if issparse(X) and self.svd_solver not in ["auto", "arpack", "covariance_eigh"]:
487 raise TypeError(
488 'PCA only support sparse inputs with the "arpack" and'
489 f' "covariance_eigh" solvers, while "{self.svd_solver}" was passed. See'
490 " TruncatedSVD for a possible alternative."
491 )
492 if self.svd_solver == "arpack" and is_array_api_compliant:
493 raise ValueError(
494 "PCA with svd_solver='arpack' is not supported for Array API inputs."
495 )
496
497 # Validate the data, without ever forcing a copy as any solver that
498 # supports sparse input data and the `covariance_eigh` solver are
499 # written in a way to avoid the need for any inplace modification of
500 # the input data contrary to the other solvers.
501 # The copy will happen
502 # later, only if needed, once the solver negotiation below is done.
503 X = validate_data(
504 self,
505 X,
506 dtype=[xp.float64, xp.float32],
507 force_writeable=True,
508 accept_sparse=("csr", "csc"),
509 ensure_2d=True,
510 copy=False,
511 )
512 self._fit_svd_solver = self.svd_solver
513 if self._fit_svd_solver == "auto" and issparse(X):
514 self._fit_svd_solver = "arpack"
515
516 if self.n_components is None:
517 if self._fit_svd_solver != "arpack":
518 n_components = min(X.shape)
519 else:
520 n_components = min(X.shape) - 1
521 else:
522 n_components = self.n_components
523
524 if self._fit_svd_solver == "auto":
525 # Tall and skinny problems are best handled by precomputing the
526 # covariance matrix.
527 if X.shape[1] <= 1_000 and X.shape[0] >= 10 * X.shape[1]:
528 self._fit_svd_solver = "covariance_eigh"
529 # Small problem or n_components == 'mle', just call full PCA
530 elif max(X.shape) <= 500 or n_components == "mle":
531 self._fit_svd_solver = "full"
532 elif 1 <= n_components < 0.8 * min(X.shape):
533 self._fit_svd_solver = "randomized"
534 # This is also the case of n_components in (0, 1)
535 else:
536 self._fit_svd_solver = "full"
537
538 # Call different fits for either full or truncated SVD
539 if self._fit_svd_solver in ("full", "covariance_eigh"):
540 return self._fit_full(X, n_components, xp, is_array_api_compliant)
541 elif self._fit_svd_solver in ["arpack", "randomized"]:
542 return self._fit_truncated(X, n_components, xp)
543
544 def _fit_full(self, X, n_components, xp, is_array_api_compliant):
545 """Fit the model by computing full SVD on X."""
546 n_samples, n_features = X.shape
547
548 if n_components == "mle":
549 if n_samples < n_features:
550 raise ValueError(
551 "n_components='mle' is only supported if n_samples >= n_features"
552 )
553 elif not 0 <= n_components <= min(n_samples, n_features):
554 raise ValueError(
555 f"n_components={n_components} must be between 0 and "
556 f"min(n_samples, n_features)={min(n_samples, n_features)} with "
557 f"svd_solver={self._fit_svd_solver!r}"
558 )
559
560 self.mean_ = xp.mean(X, axis=0)
561 # When X is a scipy sparse matrix, self.mean_ is a numpy matrix, so we need
562 # to transform it to a 1D array. Note that this is not the case when X
563 # is a scipy sparse array.
564 # TODO: remove the following two lines when scikit-learn only depends
565 # on scipy versions that no longer support scipy.sparse matrices.
566 self.mean_ = xp.reshape(xp.asarray(self.mean_), (-1,))
567
568 if self._fit_svd_solver == "full":
569 X_centered = xp.asarray(X, copy=True) if self.copy else X
570 X_centered -= self.mean_
571 x_is_centered = not self.copy
572
573 if not is_array_api_compliant:
574 # Use scipy.linalg with NumPy/SciPy inputs for the sake of not
575 # introducing unanticipated behavior changes. In the long run we
576 # could instead decide to always use xp.linalg.svd for all inputs,
577 # but that would make this code rely on numpy's SVD instead of
578 # scipy's. It's not 100% clear whether they use the same LAPACK
579 # solver by default though (assuming both are built against the
580 # same BLAS).
581 U, S, Vt = linalg.svd(X_centered, full_matrices=False)
582 else:
583 U, S, Vt = xp.linalg.svd(X_centered, full_matrices=False)
584 explained_variance_ = (S**2) / (n_samples - 1)
585
586 else:
587 assert self._fit_svd_solver == "covariance_eigh"
588 # In the following, we center the covariance matrix C afterwards
589 # (without centering the data X first) to avoid an unnecessary copy
590 # of X. Note that the mean_ attribute is still needed to center
591 # test data in the transform method.
592 #
593 # Note: at the time of writing, `xp.cov` does not exist in the
594 # Array API standard:
595 # https://github.com/data-apis/array-api/issues/43
596 #
597 # Besides, using `numpy.cov`, as of numpy 1.26.0, would not be
598 # memory efficient for our use case when `n_samples >> n_features`:
599 # `numpy.cov` centers a copy of the data before computing the
600 # matrix product instead of subtracting a small `(n_features,
601 # n_features)` square matrix from the gram matrix X.T @ X, as we do
602 # below.
603 x_is_centered = False
604 C = X.T @ X
605 C -= (
606 n_samples
607 * xp.reshape(self.mean_, (-1, 1))
608 * xp.reshape(self.mean_, (1, -1))
609 )
610 C /= n_samples - 1
611 eigenvals, eigenvecs = xp.linalg.eigh(C)
612
613 # When X is a scipy sparse matrix, the following two datastructures
614 # are returned as instances of the soft-deprecated numpy.matrix
615 # class. Note that this problem does not occur when X is a scipy
616 # sparse array (or another other kind of supported array).
617 # TODO: remove the following two lines when scikit-learn only
618 # depends on scipy versions that no longer support scipy.sparse
619 # matrices.
620 eigenvals = xp.reshape(xp.asarray(eigenvals), (-1,))
621 eigenvecs = xp.asarray(eigenvecs)
622
623 eigenvals = xp.flip(eigenvals, axis=0)
624 eigenvecs = xp.flip(eigenvecs, axis=1)
625
626 # The covariance matrix C is positive semi-definite by
627 # construction. However, the eigenvalues returned by xp.linalg.eigh
628 # can be slightly negative due to numerical errors. This would be
629 # an issue for the subsequent sqrt, hence the manual clipping.
630 eigenvals[eigenvals < 0.0] = 0.0
631 explained_variance_ = eigenvals
632
633 # Re-construct SVD of centered X indirectly and make it consistent
634 # with the other solvers.
635 S = xp.sqrt(eigenvals * (n_samples - 1))
636 Vt = eigenvecs.T
637 U = None
638
639 # flip eigenvectors' sign to enforce deterministic output
640 U, Vt = svd_flip(U, Vt, u_based_decision=False)
641
642 components_ = Vt
643
644 # Get variance explained by singular values
645 total_var = xp.sum(explained_variance_)
646 explained_variance_ratio_ = explained_variance_ / total_var
647 singular_values_ = xp.asarray(S, copy=True) # Store the singular values.
648
649 # Postprocess the number of components required
650 if n_components == "mle":
651 n_components = _infer_dimension(explained_variance_, n_samples)
652 elif 0 < n_components < 1.0:
653 # number of components for which the cumulated explained
654 # variance percentage is superior to the desired threshold
655 # side='right' ensures that number of features selected
656 # their variance is always greater than n_components float
657 # passed. More discussion in issue: #15669
658 if is_array_api_compliant:
659 # Convert to numpy as xp.cumsum and xp.searchsorted are not
660 # part of the Array API standard yet:
661 #
662 # https://github.com/data-apis/array-api/issues/597
663 # https://github.com/data-apis/array-api/issues/688
664 #
665 # Furthermore, it's not always safe to call them for namespaces
666 # that already implement them: for instance as
667 # cupy.searchsorted does not accept a float as second argument.
668 explained_variance_ratio_np = _convert_to_numpy(
669 explained_variance_ratio_, xp=xp
670 )
671 else:
672 explained_variance_ratio_np = explained_variance_ratio_
673 ratio_cumsum = stable_cumsum(explained_variance_ratio_np)
674 n_components = np.searchsorted(ratio_cumsum, n_components, side="right") + 1
675
676 # Compute noise covariance using Probabilistic PCA model
677 # The sigma2 maximum likelihood (cf. eq. 12.46)
678 if n_components < min(n_features, n_samples):
679 self.noise_variance_ = xp.mean(explained_variance_[n_components:])
680 else:
681 self.noise_variance_ = 0.0
682
683 self.n_samples_ = n_samples
684 self.n_components_ = n_components
685 # Assign a copy of the result of the truncation of the components in
686 # order to:
687 # - release the memory used by the discarded components,
688 # - ensure that the kept components are allocated contiguously in
689 # memory to make the transform method faster by leveraging cache
690 # locality.
691 self.components_ = xp.asarray(components_[:n_components, :], copy=True)
692
693 # We do the same for the other arrays for the sake of consistency.
694 self.explained_variance_ = xp.asarray(
695 explained_variance_[:n_components], copy=True
696 )
697 self.explained_variance_ratio_ = xp.asarray(
698 explained_variance_ratio_[:n_components], copy=True
699 )
700 self.singular_values_ = xp.asarray(singular_values_[:n_components], copy=True)
701
702 return U, S, Vt, X, x_is_centered, xp
703
704 def _fit_truncated(self, X, n_components, xp):
705 """Fit the model by computing truncated SVD (by ARPACK or randomized)
706 on X.
707 """
708 n_samples, n_features = X.shape
709
710 svd_solver = self._fit_svd_solver
711 if isinstance(n_components, str):
712 raise ValueError(
713 "n_components=%r cannot be a string with svd_solver='%s'"
714 % (n_components, svd_solver)
715 )
716 elif not 1 <= n_components <= min(n_samples, n_features):
717 raise ValueError(
718 "n_components=%r must be between 1 and "
719 "min(n_samples, n_features)=%r with "
720 "svd_solver='%s'"
721 % (n_components, min(n_samples, n_features), svd_solver)
722 )
723 elif svd_solver == "arpack" and n_components == min(n_samples, n_features):
724 raise ValueError(
725 "n_components=%r must be strictly less than "
726 "min(n_samples, n_features)=%r with "
727 "svd_solver='%s'"
728 % (n_components, min(n_samples, n_features), svd_solver)
729 )
730
731 random_state = check_random_state(self.random_state)
732
733 # Center data
734 total_var = None
735 if issparse(X):
736 self.mean_, var = mean_variance_axis(X, axis=0)
737 total_var = var.sum() * n_samples / (n_samples - 1) # ddof=1
738 X_centered = _implicit_column_offset(X, self.mean_)
739 x_is_centered = False
740 else:
741 self.mean_ = xp.mean(X, axis=0)
742 X_centered = xp.asarray(X, copy=True) if self.copy else X
743 X_centered -= self.mean_
744 x_is_centered = not self.copy
745
746 if svd_solver == "arpack":
747 v0 = _init_arpack_v0(min(X.shape), random_state)
748 U, S, Vt = svds(X_centered, k=n_components, tol=self.tol, v0=v0)
749 # svds doesn't abide by scipy.linalg.svd/randomized_svd
750 # conventions, so reverse its outputs.
751 S = S[::-1]
752 # flip eigenvectors' sign to enforce deterministic output
753 U, Vt = svd_flip(U[:, ::-1], Vt[::-1], u_based_decision=False)
754
755 elif svd_solver == "randomized":
756 # sign flipping is done inside
757 U, S, Vt = _randomized_svd(
758 X_centered,
759 n_components=n_components,
760 n_oversamples=self.n_oversamples,
761 n_iter=self.iterated_power,
762 power_iteration_normalizer=self.power_iteration_normalizer,
763 flip_sign=False,
764 random_state=random_state,
765 )
766 U, Vt = svd_flip(U, Vt, u_based_decision=False)
767
768 self.n_samples_ = n_samples
769 self.components_ = Vt
770 self.n_components_ = n_components
771
772 # Get variance explained by singular values
773 self.explained_variance_ = (S**2) / (n_samples - 1)
774
775 # Workaround in-place variance calculation since at the time numpy
776 # did not have a way to calculate variance in-place.
777 #
778 # TODO: update this code to either:
779 # * Use the array-api variance calculation, unless memory usage suffers
780 # * Update sklearn.utils.extmath._incremental_mean_and_var to support array-api
781 # See: https://github.com/scikit-learn/scikit-learn/pull/18689#discussion_r1335540991
782 if total_var is None:
783 N = X.shape[0] - 1
784 X_centered **= 2
785 total_var = xp.sum(X_centered) / N
786
787 self.explained_variance_ratio_ = self.explained_variance_ / total_var
788 self.singular_values_ = xp.asarray(S, copy=True) # Store the singular values.
789
790 if self.n_components_ < min(n_features, n_samples):
791 self.noise_variance_ = total_var - xp.sum(self.explained_variance_)
792 self.noise_variance_ /= min(n_features, n_samples) - n_components
793 else:
794 self.noise_variance_ = 0.0
795
796 return U, S, Vt, X, x_is_centered, xp
797
798 def score_samples(self, X):
799 """Return the log-likelihood of each sample.
800
801 See. "Pattern Recognition and Machine Learning"
802 by C. Bishop, 12.2.1 p. 574
803 or http://www.miketipping.com/papers/met-mppca.pdf
804
805 Parameters
806 ----------
807 X : array-like of shape (n_samples, n_features)
808 The data.
809
810 Returns
811 -------
812 ll : ndarray of shape (n_samples,)
813 Log-likelihood of each sample under the current model.
814 """
815 check_is_fitted(self)
816 xp, _ = get_namespace(X)
817 X = validate_data(self, X, dtype=[xp.float64, xp.float32], reset=False)
818 Xr = X - self.mean_
819 n_features = X.shape[1]
820 precision = self.get_precision()
821 log_like = -0.5 * xp.sum(Xr * (Xr @ precision), axis=1)
822 log_like -= 0.5 * (n_features * log(2.0 * np.pi) - fast_logdet(precision))
823 return log_like
824
825 def score(self, X, y=None):
826 """Return the average log-likelihood of all samples.
827
828 See. "Pattern Recognition and Machine Learning"
829 by C. Bishop, 12.2.1 p. 574
830 or http://www.miketipping.com/papers/met-mppca.pdf
831
832 Parameters
833 ----------
834 X : array-like of shape (n_samples, n_features)
835 The data.
836
837 y : Ignored
838 Ignored.
839
840 Returns
841 -------
842 ll : float
843 Average log-likelihood of the samples under the current model.
844 """
845 xp, _ = get_namespace(X)
846 return float(xp.mean(self.score_samples(X)))
847
848 def __sklearn_tags__(self):
849 tags = super().__sklearn_tags__()
850 tags.transformer_tags.preserves_dtype = ["float64", "float32"]
851 tags.array_api_support = True
852 tags.input_tags.sparse = self.svd_solver in (
853 "auto",
854 "arpack",
855 "covariance_eigh",
856 )
857 return tags
858