Aluode/PerceptionLabPortable
0
1"""Tests for Incremental PCA."""
2
3import itertools
4import warnings
5
6import numpy as np
7import pytest
8from numpy.testing import assert_allclose, assert_array_equal
9
10from sklearn import datasets
11from sklearn.decomposition import PCA, IncrementalPCA
12from sklearn.utils._testing import (
13 assert_allclose_dense_sparse,
14 assert_almost_equal,
15 assert_array_almost_equal,
16)
17from sklearn.utils.fixes import CSC_CONTAINERS, CSR_CONTAINERS, LIL_CONTAINERS
18
19iris = datasets.load_iris()
20
21
22def test_incremental_pca():
23 # Incremental PCA on dense arrays.
24 X = iris.data
25 batch_size = X.shape[0] // 3
26 ipca = IncrementalPCA(n_components=2, batch_size=batch_size)
27 pca = PCA(n_components=2)
28 pca.fit_transform(X)
29
30 X_transformed = ipca.fit_transform(X)
31
32 assert X_transformed.shape == (X.shape[0], 2)
33 np.testing.assert_allclose(
34 ipca.explained_variance_ratio_.sum(),
35 pca.explained_variance_ratio_.sum(),
36 rtol=1e-3,
37 )
38
39 for n_components in [1, 2, X.shape[1]]:
40 ipca = IncrementalPCA(n_components, batch_size=batch_size)
41 ipca.fit(X)
42 cov = ipca.get_covariance()
43 precision = ipca.get_precision()
44 np.testing.assert_allclose(
45 np.dot(cov, precision), np.eye(X.shape[1]), atol=1e-13
46 )
47
48
49@pytest.mark.parametrize(
50 "sparse_container", CSC_CONTAINERS + CSR_CONTAINERS + LIL_CONTAINERS
51)
52def test_incremental_pca_sparse(sparse_container):
53 # Incremental PCA on sparse arrays.
54 X = iris.data
55 pca = PCA(n_components=2)
56 pca.fit_transform(X)
57 X_sparse = sparse_container(X)
58 batch_size = X_sparse.shape[0] // 3
59 ipca = IncrementalPCA(n_components=2, batch_size=batch_size)
60
61 X_transformed = ipca.fit_transform(X_sparse)
62
63 assert X_transformed.shape == (X_sparse.shape[0], 2)
64 np.testing.assert_allclose(
65 ipca.explained_variance_ratio_.sum(),
66 pca.explained_variance_ratio_.sum(),
67 rtol=1e-3,
68 )
69
70 for n_components in [1, 2, X.shape[1]]:
71 ipca = IncrementalPCA(n_components, batch_size=batch_size)
72 ipca.fit(X_sparse)
73 cov = ipca.get_covariance()
74 precision = ipca.get_precision()
75 np.testing.assert_allclose(
76 np.dot(cov, precision), np.eye(X_sparse.shape[1]), atol=1e-13
77 )
78
79 with pytest.raises(
80 TypeError,
81 match=(
82 "IncrementalPCA.partial_fit does not support "
83 "sparse input. Either convert data to dense "
84 "or use IncrementalPCA.fit to do so in batches."
85 ),
86 ):
87 ipca.partial_fit(X_sparse)
88
89
90def test_incremental_pca_check_projection(global_random_seed):
91 # Test that the projection of data is correct.
92 rng = np.random.RandomState(global_random_seed)
93 n, p = 100, 3
94 X = rng.randn(n, p) * 0.1
95 X[:10] += np.array([3, 4, 5])
96 Xt = 0.1 * rng.randn(1, p) + np.array([3, 4, 5])
97
98 # Get the reconstruction of the generated data X
99 # Note that Xt has the same "components" as X, just separated
100 # This is what we want to ensure is recreated correctly
101 Yt = IncrementalPCA(n_components=2).fit(X).transform(Xt)
102
103 # Normalize
104 Yt /= np.sqrt((Yt**2).sum())
105
106 # Make sure that the first element of Yt is ~1, this means
107 # the reconstruction worked as expected
108 assert_almost_equal(np.abs(Yt[0][0]), 1.0, 1)
109
110
111def test_incremental_pca_inverse(global_random_seed):
112 # Test that the projection of data can be inverted.
113 rng = np.random.RandomState(global_random_seed)
114 n, p = 50, 3
115 X = rng.randn(n, p) # spherical data
116 X[:, 1] *= 0.00001 # make middle component relatively small
117 X += [5, 4, 3] # make a large mean
118
119 # same check that we can find the original data from the transformed
120 # signal (since the data is almost of rank n_components)
121 ipca = IncrementalPCA(n_components=2, batch_size=10).fit(X)
122 Y = ipca.transform(X)
123 Y_inverse = ipca.inverse_transform(Y)
124 assert_almost_equal(X, Y_inverse, decimal=3)
125
126
127def test_incremental_pca_validation():
128 # Test that n_components is <= n_features.
129 X = np.array([[0, 1, 0], [1, 0, 0]])
130 n_samples, n_features = X.shape
131 n_components = 4
132 with pytest.raises(
133 ValueError,
134 match=(
135 "n_components={} invalid"
136 " for n_features={}, need more rows than"
137 " columns for IncrementalPCA"
138 " processing".format(n_components, n_features)
139 ),
140 ):
141 IncrementalPCA(n_components, batch_size=10).fit(X)
142
143 # Test that n_components is also <= n_samples in first call to partial fit.
144 n_components = 3
145 with pytest.raises(
146 ValueError,
147 match=(
148 f"n_components={n_components} must be less or equal to the batch "
149 f"number of samples {n_samples} for the first partial_fit call."
150 ),
151 ):
152 IncrementalPCA(n_components=n_components).partial_fit(X)
153
154
155def test_n_samples_equal_n_components():
156 # Ensures no warning is raised when n_samples==n_components
157 # Non-regression test for gh-19050
158 ipca = IncrementalPCA(n_components=5)
159 with warnings.catch_warnings():
160 warnings.simplefilter("error", RuntimeWarning)
161 ipca.partial_fit(np.random.randn(5, 7))
162 with warnings.catch_warnings():
163 warnings.simplefilter("error", RuntimeWarning)
164 ipca.fit(np.random.randn(5, 7))
165
166
167def test_n_components_none():
168 # Ensures that n_components == None is handled correctly
169 rng = np.random.RandomState(1999)
170 for n_samples, n_features in [(50, 10), (10, 50)]:
171 X = rng.rand(n_samples, n_features)
172 ipca = IncrementalPCA(n_components=None)
173
174 # First partial_fit call, ipca.n_components_ is inferred from
175 # min(X.shape)
176 ipca.partial_fit(X)
177 assert ipca.n_components_ == min(X.shape)
178
179 # Second partial_fit call, ipca.n_components_ is inferred from
180 # ipca.components_ computed from the first partial_fit call
181 ipca.partial_fit(X)
182 assert ipca.n_components_ == ipca.components_.shape[0]
183
184
185def test_incremental_pca_set_params():
186 # Test that components_ sign is stable over batch sizes.
187 rng = np.random.RandomState(1999)
188 n_samples = 100
189 n_features = 20
190 X = rng.randn(n_samples, n_features)
191 X2 = rng.randn(n_samples, n_features)
192 X3 = rng.randn(n_samples, n_features)
193 ipca = IncrementalPCA(n_components=20)
194 ipca.fit(X)
195 # Decreasing number of components
196 ipca.set_params(n_components=10)
197 with pytest.raises(ValueError):
198 ipca.partial_fit(X2)
199 # Increasing number of components
200 ipca.set_params(n_components=15)
201 with pytest.raises(ValueError):
202 ipca.partial_fit(X3)
203 # Returning to original setting
204 ipca.set_params(n_components=20)
205 ipca.partial_fit(X)
206
207
208def test_incremental_pca_num_features_change():
209 # Test that changing n_components will raise an error.
210 rng = np.random.RandomState(1999)
211 n_samples = 100
212 X = rng.randn(n_samples, 20)
213 X2 = rng.randn(n_samples, 50)
214 ipca = IncrementalPCA(n_components=None)
215 ipca.fit(X)
216 with pytest.raises(ValueError):
217 ipca.partial_fit(X2)
218
219
220def test_incremental_pca_batch_signs(global_random_seed):
221 # Test that components_ sign is stable over batch sizes.
222 rng = np.random.RandomState(global_random_seed)
223 n_samples = 100
224 n_features = 3
225 X = rng.randn(n_samples, n_features)
226 all_components = []
227 batch_sizes = np.arange(10, 20)
228 for batch_size in batch_sizes:
229 ipca = IncrementalPCA(n_components=None, batch_size=batch_size).fit(X)
230 all_components.append(ipca.components_)
231
232 for i, j in itertools.pairwise(all_components):
233 assert_almost_equal(np.sign(i), np.sign(j), decimal=6)
234
235
236def test_incremental_pca_partial_fit_small_batch():
237 # Test that there is no minimum batch size after the first partial_fit
238 # Non-regression test
239 rng = np.random.RandomState(1999)
240 n, p = 50, 3
241 X = rng.randn(n, p) # spherical data
242 X[:, 1] *= 0.00001 # make middle component relatively small
243 X += [5, 4, 3] # make a large mean
244
245 n_components = p
246 pipca = IncrementalPCA(n_components=n_components)
247 pipca.partial_fit(X[:n_components])
248 for idx in range(n_components, n):
249 pipca.partial_fit(X[idx : idx + 1])
250
251 pca = PCA(n_components=n_components)
252 pca.fit(X)
253
254 assert_allclose(pca.components_, pipca.components_, atol=1e-3)
255
256
257def test_incremental_pca_batch_values(global_random_seed):
258 # Test that components_ values are stable over batch sizes.
259 rng = np.random.RandomState(global_random_seed)
260 n_samples = 100
261 n_features = 3
262 X = rng.randn(n_samples, n_features)
263 all_components = []
264 batch_sizes = np.arange(20, 40, 3)
265 for batch_size in batch_sizes:
266 ipca = IncrementalPCA(n_components=None, batch_size=batch_size).fit(X)
267 all_components.append(ipca.components_)
268
269 for i, j in itertools.pairwise(all_components):
270 assert_almost_equal(i, j, decimal=1)
271
272
273def test_incremental_pca_batch_rank():
274 # Test sample size in each batch is always larger or equal to n_components
275 rng = np.random.RandomState(1999)
276 n_samples = 100
277 n_features = 20
278 X = rng.randn(n_samples, n_features)
279 all_components = []
280 batch_sizes = np.arange(20, 90, 3)
281 for batch_size in batch_sizes:
282 ipca = IncrementalPCA(n_components=20, batch_size=batch_size).fit(X)
283 all_components.append(ipca.components_)
284
285 for components_i, components_j in itertools.pairwise(all_components):
286 assert_allclose_dense_sparse(components_i, components_j)
287
288
289def test_incremental_pca_partial_fit(global_random_seed):
290 # Test that fit and partial_fit get equivalent results.
291 rng = np.random.RandomState(global_random_seed)
292 n, p = 50, 3
293 X = rng.randn(n, p) # spherical data
294 X[:, 1] *= 0.00001 # make middle component relatively small
295 X += [5, 4, 3] # make a large mean
296
297 # same check that we can find the original data from the transformed
298 # signal (since the data is almost of rank n_components)
299 batch_size = 10
300 ipca = IncrementalPCA(n_components=2, batch_size=batch_size).fit(X)
301 pipca = IncrementalPCA(n_components=2, batch_size=batch_size)
302 # Add one to make sure endpoint is included
303 batch_itr = np.arange(0, n + 1, batch_size)
304 for i, j in itertools.pairwise(batch_itr):
305 pipca.partial_fit(X[i:j, :])
306 assert_almost_equal(ipca.components_, pipca.components_, decimal=3)
307
308
309def test_incremental_pca_against_pca_iris():
310 # Test that IncrementalPCA and PCA are approximate (to a sign flip).
311 X = iris.data
312
313 Y_pca = PCA(n_components=2).fit_transform(X)
314 Y_ipca = IncrementalPCA(n_components=2, batch_size=25).fit_transform(X)
315
316 assert_almost_equal(np.abs(Y_pca), np.abs(Y_ipca), 1)
317
318
319def test_incremental_pca_against_pca_random_data(global_random_seed):
320 # Test that IncrementalPCA and PCA are approximate (to a sign flip).
321 rng = np.random.RandomState(global_random_seed)
322 n_samples = 100
323 n_features = 3
324 X = rng.randn(n_samples, n_features) + 5 * rng.rand(1, n_features)
325
326 Y_pca = PCA(n_components=3).fit_transform(X)
327 Y_ipca = IncrementalPCA(n_components=3, batch_size=25).fit_transform(X)
328
329 assert_almost_equal(np.abs(Y_pca), np.abs(Y_ipca), 1)
330
331
332def test_explained_variances():
333 # Test that PCA and IncrementalPCA calculations match
334 X = datasets.make_low_rank_matrix(
335 1000, 100, tail_strength=0.0, effective_rank=10, random_state=1999
336 )
337 prec = 3
338 n_samples, n_features = X.shape
339 for nc in [None, 99]:
340 pca = PCA(n_components=nc).fit(X)
341 ipca = IncrementalPCA(n_components=nc, batch_size=100).fit(X)
342 assert_almost_equal(
343 pca.explained_variance_, ipca.explained_variance_, decimal=prec
344 )
345 assert_almost_equal(
346 pca.explained_variance_ratio_, ipca.explained_variance_ratio_, decimal=prec
347 )
348 assert_almost_equal(pca.noise_variance_, ipca.noise_variance_, decimal=prec)
349
350
351def test_singular_values(global_random_seed):
352 # Check that the IncrementalPCA output has the correct singular values
353
354 rng = np.random.RandomState(global_random_seed)
355 n_samples = 1000
356 n_features = 100
357
358 X = datasets.make_low_rank_matrix(
359 n_samples, n_features, tail_strength=0.0, effective_rank=10, random_state=rng
360 )
361
362 pca = PCA(n_components=10, svd_solver="full", random_state=rng).fit(X)
363 ipca = IncrementalPCA(n_components=10, batch_size=150).fit(X)
364 assert_array_almost_equal(pca.singular_values_, ipca.singular_values_, 2)
365
366 # Compare to the Frobenius norm
367 X_pca = pca.transform(X)
368 X_ipca = ipca.transform(X)
369 assert_array_almost_equal(
370 np.sum(pca.singular_values_**2.0), np.linalg.norm(X_pca, "fro") ** 2.0, 12
371 )
372 assert_array_almost_equal(
373 np.sum(ipca.singular_values_**2.0), np.linalg.norm(X_ipca, "fro") ** 2.0, 2
374 )
375
376 # Compare to the 2-norms of the score vectors
377 assert_array_almost_equal(
378 pca.singular_values_, np.sqrt(np.sum(X_pca**2.0, axis=0)), 12
379 )
380 assert_array_almost_equal(
381 ipca.singular_values_, np.sqrt(np.sum(X_ipca**2.0, axis=0)), 2
382 )
383
384 # Set the singular values and see what we get back
385 rng = np.random.RandomState(global_random_seed)
386 n_samples = 100
387 n_features = 110
388
389 X = datasets.make_low_rank_matrix(
390 n_samples, n_features, tail_strength=0.0, effective_rank=3, random_state=rng
391 )
392
393 pca = PCA(n_components=3, svd_solver="full", random_state=rng)
394 ipca = IncrementalPCA(n_components=3, batch_size=100)
395
396 X_pca = pca.fit_transform(X)
397 X_pca /= np.sqrt(np.sum(X_pca**2.0, axis=0))
398 X_pca[:, 0] *= 3.142
399 X_pca[:, 1] *= 2.718
400
401 X_hat = np.dot(X_pca, pca.components_)
402 pca.fit(X_hat)
403 ipca.fit(X_hat)
404 assert_array_almost_equal(pca.singular_values_, [3.142, 2.718, 1.0], 14)
405 assert_array_almost_equal(ipca.singular_values_, [3.142, 2.718, 1.0], 14)
406
407
408def test_whitening(global_random_seed):
409 # Test that PCA and IncrementalPCA transforms match to sign flip.
410 X = datasets.make_low_rank_matrix(
411 1000, 10, tail_strength=0.0, effective_rank=2, random_state=global_random_seed
412 )
413 atol = 1e-3
414 for nc in [None, 9]:
415 pca = PCA(whiten=True, n_components=nc).fit(X)
416 ipca = IncrementalPCA(whiten=True, n_components=nc, batch_size=250).fit(X)
417
418 # Since the data is rank deficient, some components are pure noise. We
419 # should not expect those dimensions to carry any signal and their
420 # values might be arbitrarily changed by implementation details of the
421 # internal SVD solver. We therefore filter them out before comparison.
422 stable_mask = pca.explained_variance_ratio_ > 1e-12
423
424 Xt_pca = pca.transform(X)
425 Xt_ipca = ipca.transform(X)
426 assert_allclose(
427 np.abs(Xt_pca)[:, stable_mask],
428 np.abs(Xt_ipca)[:, stable_mask],
429 atol=atol,
430 )
431
432 # The noisy dimensions are in the null space of the inverse transform,
433 # so they are not influencing the reconstruction. We therefore don't
434 # need to apply the mask here.
435 Xinv_ipca = ipca.inverse_transform(Xt_ipca)
436 Xinv_pca = pca.inverse_transform(Xt_pca)
437 assert_allclose(X, Xinv_ipca, atol=atol)
438 assert_allclose(X, Xinv_pca, atol=atol)
439 assert_allclose(Xinv_pca, Xinv_ipca, atol=atol)
440
441
442def test_incremental_pca_partial_fit_float_division():
443 # Test to ensure float division is used in all versions of Python
444 # (non-regression test for issue #9489)
445
446 rng = np.random.RandomState(0)
447 A = rng.randn(5, 3) + 2
448 B = rng.randn(7, 3) + 5
449
450 pca = IncrementalPCA(n_components=2)
451 pca.partial_fit(A)
452 # Set n_samples_seen_ to be a floating point number instead of an int
453 pca.n_samples_seen_ = float(pca.n_samples_seen_)
454 pca.partial_fit(B)
455 singular_vals_float_samples_seen = pca.singular_values_
456
457 pca2 = IncrementalPCA(n_components=2)
458 pca2.partial_fit(A)
459 pca2.partial_fit(B)
460 singular_vals_int_samples_seen = pca2.singular_values_
461
462 np.testing.assert_allclose(
463 singular_vals_float_samples_seen, singular_vals_int_samples_seen
464 )
465
466
467def test_incremental_pca_fit_overflow_error():
468 # Test for overflow error on Windows OS
469 # (non-regression test for issue #17693)
470 rng = np.random.RandomState(0)
471 A = rng.rand(500000, 2)
472
473 ipca = IncrementalPCA(n_components=2, batch_size=10000)
474 ipca.fit(A)
475
476 pca = PCA(n_components=2)
477 pca.fit(A)
478
479 np.testing.assert_allclose(ipca.singular_values_, pca.singular_values_)
480
481
482def test_incremental_pca_feature_names_out():
483 """Check feature names out for IncrementalPCA."""
484 ipca = IncrementalPCA(n_components=2).fit(iris.data)
485
486 names = ipca.get_feature_names_out()
487 assert_array_equal([f"incrementalpca{i}" for i in range(2)], names)
488