_incremental_pca.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406
  1. """Incremental Principal Components Analysis."""
  2. # Author: Kyle Kastner <kastnerkyle@gmail.com>
  3. # Giorgio Patrini
  4. # License: BSD 3 clause
  5. from numbers import Integral
  6. import numpy as np
  7. from scipy import linalg, sparse
  8. from ..base import _fit_context
  9. from ..utils import gen_batches
  10. from ..utils._param_validation import Interval
  11. from ..utils.extmath import _incremental_mean_and_var, svd_flip
  12. from ._base import _BasePCA
  13. class IncrementalPCA(_BasePCA):
  14. """Incremental principal components analysis (IPCA).
  15. Linear dimensionality reduction using Singular Value Decomposition of
  16. the data, keeping only the most significant singular vectors to
  17. project the data to a lower dimensional space. The input data is centered
  18. but not scaled for each feature before applying the SVD.
  19. Depending on the size of the input data, this algorithm can be much more
  20. memory efficient than a PCA, and allows sparse input.
  21. This algorithm has constant memory complexity, on the order
  22. of ``batch_size * n_features``, enabling use of np.memmap files without
  23. loading the entire file into memory. For sparse matrices, the input
  24. is converted to dense in batches (in order to be able to subtract the
  25. mean) which avoids storing the entire dense matrix at any one time.
  26. The computational overhead of each SVD is
  27. ``O(batch_size * n_features ** 2)``, but only 2 * batch_size samples
  28. remain in memory at a time. There will be ``n_samples / batch_size`` SVD
  29. computations to get the principal components, versus 1 large SVD of
  30. complexity ``O(n_samples * n_features ** 2)`` for PCA.
  31. Read more in the :ref:`User Guide <IncrementalPCA>`.
  32. .. versionadded:: 0.16
  33. Parameters
  34. ----------
  35. n_components : int, default=None
  36. Number of components to keep. If ``n_components`` is ``None``,
  37. then ``n_components`` is set to ``min(n_samples, n_features)``.
  38. whiten : bool, default=False
  39. When True (False by default) the ``components_`` vectors are divided
  40. by ``n_samples`` times ``components_`` to ensure uncorrelated outputs
  41. with unit component-wise variances.
  42. Whitening will remove some information from the transformed signal
  43. (the relative variance scales of the components) but can sometimes
  44. improve the predictive accuracy of the downstream estimators by
  45. making data respect some hard-wired assumptions.
  46. copy : bool, default=True
  47. If False, X will be overwritten. ``copy=False`` can be used to
  48. save memory but is unsafe for general use.
  49. batch_size : int, default=None
  50. The number of samples to use for each batch. Only used when calling
  51. ``fit``. If ``batch_size`` is ``None``, then ``batch_size``
  52. is inferred from the data and set to ``5 * n_features``, to provide a
  53. balance between approximation accuracy and memory consumption.
  54. Attributes
  55. ----------
  56. components_ : ndarray of shape (n_components, n_features)
  57. Principal axes in feature space, representing the directions of
  58. maximum variance in the data. Equivalently, the right singular
  59. vectors of the centered input data, parallel to its eigenvectors.
  60. The components are sorted by decreasing ``explained_variance_``.
  61. explained_variance_ : ndarray of shape (n_components,)
  62. Variance explained by each of the selected components.
  63. explained_variance_ratio_ : ndarray of shape (n_components,)
  64. Percentage of variance explained by each of the selected components.
  65. If all components are stored, the sum of explained variances is equal
  66. to 1.0.
  67. singular_values_ : ndarray of shape (n_components,)
  68. The singular values corresponding to each of the selected components.
  69. The singular values are equal to the 2-norms of the ``n_components``
  70. variables in the lower-dimensional space.
  71. mean_ : ndarray of shape (n_features,)
  72. Per-feature empirical mean, aggregate over calls to ``partial_fit``.
  73. var_ : ndarray of shape (n_features,)
  74. Per-feature empirical variance, aggregate over calls to
  75. ``partial_fit``.
  76. noise_variance_ : float
  77. The estimated noise covariance following the Probabilistic PCA model
  78. from Tipping and Bishop 1999. See "Pattern Recognition and
  79. Machine Learning" by C. Bishop, 12.2.1 p. 574 or
  80. http://www.miketipping.com/papers/met-mppca.pdf.
  81. n_components_ : int
  82. The estimated number of components. Relevant when
  83. ``n_components=None``.
  84. n_samples_seen_ : int
  85. The number of samples processed by the estimator. Will be reset on
  86. new calls to fit, but increments across ``partial_fit`` calls.
  87. batch_size_ : int
  88. Inferred batch size from ``batch_size``.
  89. n_features_in_ : int
  90. Number of features seen during :term:`fit`.
  91. .. versionadded:: 0.24
  92. feature_names_in_ : ndarray of shape (`n_features_in_`,)
  93. Names of features seen during :term:`fit`. Defined only when `X`
  94. has feature names that are all strings.
  95. .. versionadded:: 1.0
  96. See Also
  97. --------
  98. PCA : Principal component analysis (PCA).
  99. KernelPCA : Kernel Principal component analysis (KPCA).
  100. SparsePCA : Sparse Principal Components Analysis (SparsePCA).
  101. TruncatedSVD : Dimensionality reduction using truncated SVD.
  102. Notes
  103. -----
  104. Implements the incremental PCA model from:
  105. *D. Ross, J. Lim, R. Lin, M. Yang, Incremental Learning for Robust Visual
  106. Tracking, International Journal of Computer Vision, Volume 77, Issue 1-3,
  107. pp. 125-141, May 2008.*
  108. See https://www.cs.toronto.edu/~dross/ivt/RossLimLinYang_ijcv.pdf
  109. This model is an extension of the Sequential Karhunen-Loeve Transform from:
  110. :doi:`A. Levy and M. Lindenbaum, Sequential Karhunen-Loeve Basis Extraction and
  111. its Application to Images, IEEE Transactions on Image Processing, Volume 9,
  112. Number 8, pp. 1371-1374, August 2000. <10.1109/83.855432>`
  113. We have specifically abstained from an optimization used by authors of both
  114. papers, a QR decomposition used in specific situations to reduce the
  115. algorithmic complexity of the SVD. The source for this technique is
  116. *Matrix Computations, Third Edition, G. Holub and C. Van Loan, Chapter 5,
  117. section 5.4.4, pp 252-253.*. This technique has been omitted because it is
  118. advantageous only when decomposing a matrix with ``n_samples`` (rows)
  119. >= 5/3 * ``n_features`` (columns), and hurts the readability of the
  120. implemented algorithm. This would be a good opportunity for future
  121. optimization, if it is deemed necessary.
  122. References
  123. ----------
  124. D. Ross, J. Lim, R. Lin, M. Yang. Incremental Learning for Robust Visual
  125. Tracking, International Journal of Computer Vision, Volume 77,
  126. Issue 1-3, pp. 125-141, May 2008.
  127. G. Golub and C. Van Loan. Matrix Computations, Third Edition, Chapter 5,
  128. Section 5.4.4, pp. 252-253.
  129. Examples
  130. --------
  131. >>> from sklearn.datasets import load_digits
  132. >>> from sklearn.decomposition import IncrementalPCA
  133. >>> from scipy import sparse
  134. >>> X, _ = load_digits(return_X_y=True)
  135. >>> transformer = IncrementalPCA(n_components=7, batch_size=200)
  136. >>> # either partially fit on smaller batches of data
  137. >>> transformer.partial_fit(X[:100, :])
  138. IncrementalPCA(batch_size=200, n_components=7)
  139. >>> # or let the fit function itself divide the data into batches
  140. >>> X_sparse = sparse.csr_matrix(X)
  141. >>> X_transformed = transformer.fit_transform(X_sparse)
  142. >>> X_transformed.shape
  143. (1797, 7)
  144. """
  145. _parameter_constraints: dict = {
  146. "n_components": [Interval(Integral, 1, None, closed="left"), None],
  147. "whiten": ["boolean"],
  148. "copy": ["boolean"],
  149. "batch_size": [Interval(Integral, 1, None, closed="left"), None],
  150. }
  151. def __init__(self, n_components=None, *, whiten=False, copy=True, batch_size=None):
  152. self.n_components = n_components
  153. self.whiten = whiten
  154. self.copy = copy
  155. self.batch_size = batch_size
  156. @_fit_context(prefer_skip_nested_validation=True)
  157. def fit(self, X, y=None):
  158. """Fit the model with X, using minibatches of size batch_size.
  159. Parameters
  160. ----------
  161. X : {array-like, sparse matrix} of shape (n_samples, n_features)
  162. Training data, where `n_samples` is the number of samples and
  163. `n_features` is the number of features.
  164. y : Ignored
  165. Not used, present for API consistency by convention.
  166. Returns
  167. -------
  168. self : object
  169. Returns the instance itself.
  170. """
  171. self.components_ = None
  172. self.n_samples_seen_ = 0
  173. self.mean_ = 0.0
  174. self.var_ = 0.0
  175. self.singular_values_ = None
  176. self.explained_variance_ = None
  177. self.explained_variance_ratio_ = None
  178. self.noise_variance_ = None
  179. X = self._validate_data(
  180. X,
  181. accept_sparse=["csr", "csc", "lil"],
  182. copy=self.copy,
  183. dtype=[np.float64, np.float32],
  184. )
  185. n_samples, n_features = X.shape
  186. if self.batch_size is None:
  187. self.batch_size_ = 5 * n_features
  188. else:
  189. self.batch_size_ = self.batch_size
  190. for batch in gen_batches(
  191. n_samples, self.batch_size_, min_batch_size=self.n_components or 0
  192. ):
  193. X_batch = X[batch]
  194. if sparse.issparse(X_batch):
  195. X_batch = X_batch.toarray()
  196. self.partial_fit(X_batch, check_input=False)
  197. return self
  198. @_fit_context(prefer_skip_nested_validation=True)
  199. def partial_fit(self, X, y=None, check_input=True):
  200. """Incremental fit with X. All of X is processed as a single batch.
  201. Parameters
  202. ----------
  203. X : array-like of shape (n_samples, n_features)
  204. Training data, where `n_samples` is the number of samples and
  205. `n_features` is the number of features.
  206. y : Ignored
  207. Not used, present for API consistency by convention.
  208. check_input : bool, default=True
  209. Run check_array on X.
  210. Returns
  211. -------
  212. self : object
  213. Returns the instance itself.
  214. """
  215. first_pass = not hasattr(self, "components_")
  216. if check_input:
  217. if sparse.issparse(X):
  218. raise TypeError(
  219. "IncrementalPCA.partial_fit does not support "
  220. "sparse input. Either convert data to dense "
  221. "or use IncrementalPCA.fit to do so in batches."
  222. )
  223. X = self._validate_data(
  224. X, copy=self.copy, dtype=[np.float64, np.float32], reset=first_pass
  225. )
  226. n_samples, n_features = X.shape
  227. if first_pass:
  228. self.components_ = None
  229. if self.n_components is None:
  230. if self.components_ is None:
  231. self.n_components_ = min(n_samples, n_features)
  232. else:
  233. self.n_components_ = self.components_.shape[0]
  234. elif not self.n_components <= n_features:
  235. raise ValueError(
  236. "n_components=%r invalid for n_features=%d, need "
  237. "more rows than columns for IncrementalPCA "
  238. "processing" % (self.n_components, n_features)
  239. )
  240. elif not self.n_components <= n_samples:
  241. raise ValueError(
  242. "n_components=%r must be less or equal to "
  243. "the batch number of samples "
  244. "%d." % (self.n_components, n_samples)
  245. )
  246. else:
  247. self.n_components_ = self.n_components
  248. if (self.components_ is not None) and (
  249. self.components_.shape[0] != self.n_components_
  250. ):
  251. raise ValueError(
  252. "Number of input features has changed from %i "
  253. "to %i between calls to partial_fit! Try "
  254. "setting n_components to a fixed value."
  255. % (self.components_.shape[0], self.n_components_)
  256. )
  257. # This is the first partial_fit
  258. if not hasattr(self, "n_samples_seen_"):
  259. self.n_samples_seen_ = 0
  260. self.mean_ = 0.0
  261. self.var_ = 0.0
  262. # Update stats - they are 0 if this is the first step
  263. col_mean, col_var, n_total_samples = _incremental_mean_and_var(
  264. X,
  265. last_mean=self.mean_,
  266. last_variance=self.var_,
  267. last_sample_count=np.repeat(self.n_samples_seen_, X.shape[1]),
  268. )
  269. n_total_samples = n_total_samples[0]
  270. # Whitening
  271. if self.n_samples_seen_ == 0:
  272. # If it is the first step, simply whiten X
  273. X -= col_mean
  274. else:
  275. col_batch_mean = np.mean(X, axis=0)
  276. X -= col_batch_mean
  277. # Build matrix of combined previous basis and new data
  278. mean_correction = np.sqrt(
  279. (self.n_samples_seen_ / n_total_samples) * n_samples
  280. ) * (self.mean_ - col_batch_mean)
  281. X = np.vstack(
  282. (
  283. self.singular_values_.reshape((-1, 1)) * self.components_,
  284. X,
  285. mean_correction,
  286. )
  287. )
  288. U, S, Vt = linalg.svd(X, full_matrices=False, check_finite=False)
  289. U, Vt = svd_flip(U, Vt, u_based_decision=False)
  290. explained_variance = S**2 / (n_total_samples - 1)
  291. explained_variance_ratio = S**2 / np.sum(col_var * n_total_samples)
  292. self.n_samples_seen_ = n_total_samples
  293. self.components_ = Vt[: self.n_components_]
  294. self.singular_values_ = S[: self.n_components_]
  295. self.mean_ = col_mean
  296. self.var_ = col_var
  297. self.explained_variance_ = explained_variance[: self.n_components_]
  298. self.explained_variance_ratio_ = explained_variance_ratio[: self.n_components_]
  299. # we already checked `self.n_components <= n_samples` above
  300. if self.n_components_ not in (n_samples, n_features):
  301. self.noise_variance_ = explained_variance[self.n_components_ :].mean()
  302. else:
  303. self.noise_variance_ = 0.0
  304. return self
  305. def transform(self, X):
  306. """Apply dimensionality reduction to X.
  307. X is projected on the first principal components previously extracted
  308. from a training set, using minibatches of size batch_size if X is
  309. sparse.
  310. Parameters
  311. ----------
  312. X : {array-like, sparse matrix} of shape (n_samples, n_features)
  313. New data, where `n_samples` is the number of samples
  314. and `n_features` is the number of features.
  315. Returns
  316. -------
  317. X_new : ndarray of shape (n_samples, n_components)
  318. Projection of X in the first principal components.
  319. Examples
  320. --------
  321. >>> import numpy as np
  322. >>> from sklearn.decomposition import IncrementalPCA
  323. >>> X = np.array([[-1, -1], [-2, -1], [-3, -2],
  324. ... [1, 1], [2, 1], [3, 2]])
  325. >>> ipca = IncrementalPCA(n_components=2, batch_size=3)
  326. >>> ipca.fit(X)
  327. IncrementalPCA(batch_size=3, n_components=2)
  328. >>> ipca.transform(X) # doctest: +SKIP
  329. """
  330. if sparse.issparse(X):
  331. n_samples = X.shape[0]
  332. output = []
  333. for batch in gen_batches(
  334. n_samples, self.batch_size_, min_batch_size=self.n_components or 0
  335. ):
  336. output.append(super().transform(X[batch].toarray()))
  337. return np.vstack(output)
  338. else:
  339. return super().transform(X)