Aluode/PerceptionLabPortable

1# Authors: The scikit-learn developers
2# SPDX-License-Identifier: BSD-3-Clause
3
4
5import numpy as np
6import pytest
7from numpy.testing import assert_array_equal
8
9from sklearn.datasets import make_low_rank_matrix
10from sklearn.decomposition import PCA, MiniBatchSparsePCA, SparsePCA
11from sklearn.utils import check_random_state
12from sklearn.utils._testing import (
13    assert_allclose,
14    assert_array_almost_equal,
15    if_safe_multiprocessing_with_blas,
16)
17from sklearn.utils.extmath import svd_flip
18
19
20def generate_toy_data(n_components, n_samples, image_size, random_state=None):
21    n_features = image_size[0] * image_size[1]
22
23    rng = check_random_state(random_state)
24    U = rng.randn(n_samples, n_components)
25    V = rng.randn(n_components, n_features)
26
27    centers = [(3, 3), (6, 7), (8, 1)]
28    sz = [1, 2, 1]
29    for k in range(n_components):
30        img = np.zeros(image_size)
31        xmin, xmax = centers[k][0] - sz[k], centers[k][0] + sz[k]
32        ymin, ymax = centers[k][1] - sz[k], centers[k][1] + sz[k]
33        img[xmin:xmax][:, ymin:ymax] = 1.0
34        V[k, :] = img.ravel()
35
36    # Y is defined by : Y = UV + noise
37    Y = np.dot(U, V)
38    Y += 0.1 * rng.randn(Y.shape[0], Y.shape[1])  # Add noise
39    return Y, U, V
40
41
42# SparsePCA can be a bit slow. To avoid having test times go up, we
43# test different aspects of the code in the same test
44
45
46def test_correct_shapes():
47    rng = np.random.RandomState(0)
48    X = rng.randn(12, 10)
49    spca = SparsePCA(n_components=8, random_state=rng)
50    U = spca.fit_transform(X)
51    assert spca.components_.shape == (8, 10)
52    assert U.shape == (12, 8)
53    # test overcomplete decomposition
54    spca = SparsePCA(n_components=13, random_state=rng)
55    U = spca.fit_transform(X)
56    assert spca.components_.shape == (13, 10)
57    assert U.shape == (12, 13)
58
59
60def test_fit_transform(global_random_seed):
61    alpha = 1
62    rng = np.random.RandomState(global_random_seed)
63    Y, _, _ = generate_toy_data(3, 10, (8, 8), random_state=rng)  # wide array
64    spca_lars = SparsePCA(
65        n_components=3, method="lars", alpha=alpha, random_state=global_random_seed
66    )
67    spca_lars.fit(Y)
68
69    # Test that CD gives similar results
70    spca_lasso = SparsePCA(
71        n_components=3, method="cd", random_state=global_random_seed, alpha=alpha
72    )
73    spca_lasso.fit(Y)
74    assert_array_almost_equal(spca_lasso.components_, spca_lars.components_)
75
76
77@if_safe_multiprocessing_with_blas
78def test_fit_transform_parallel(global_random_seed):
79    alpha = 1
80    rng = np.random.RandomState(global_random_seed)
81    Y, _, _ = generate_toy_data(3, 10, (8, 8), random_state=rng)  # wide array
82    spca_lars = SparsePCA(
83        n_components=3, method="lars", alpha=alpha, random_state=global_random_seed
84    )
85    spca_lars.fit(Y)
86    U1 = spca_lars.transform(Y)
87    # Test multiple CPUs
88    spca = SparsePCA(
89        n_components=3,
90        n_jobs=2,
91        method="lars",
92        alpha=alpha,
93        random_state=global_random_seed,
94    ).fit(Y)
95    U2 = spca.transform(Y)
96    assert not np.all(spca_lars.components_ == 0)
97    assert_array_almost_equal(U1, U2)
98
99
100def test_transform_nan(global_random_seed):
101    # Test that SparsePCA won't return NaN when there is 0 feature in all
102    # samples.
103    rng = np.random.RandomState(global_random_seed)
104    Y, _, _ = generate_toy_data(3, 10, (8, 8), random_state=rng)  # wide array
105    Y[:, 0] = 0
106    estimator = SparsePCA(n_components=8, random_state=global_random_seed)
107    assert not np.any(np.isnan(estimator.fit_transform(Y)))
108
109
110def test_fit_transform_tall(global_random_seed):
111    rng = np.random.RandomState(global_random_seed)
112    Y, _, _ = generate_toy_data(3, 65, (8, 8), random_state=rng)  # tall array
113    spca_lars = SparsePCA(n_components=3, method="lars", random_state=rng)
114    U1 = spca_lars.fit_transform(Y)
115    spca_lasso = SparsePCA(n_components=3, method="cd", random_state=rng)
116    U2 = spca_lasso.fit(Y).transform(Y)
117    assert_array_almost_equal(U1, U2)
118
119
120def test_initialization(global_random_seed):
121    rng = np.random.RandomState(global_random_seed)
122    U_init = rng.randn(5, 3)
123    V_init = rng.randn(3, 4)
124    model = SparsePCA(
125        n_components=3, U_init=U_init, V_init=V_init, max_iter=0, random_state=rng
126    )
127    model.fit(rng.randn(5, 4))
128
129    expected_components = V_init / np.linalg.norm(V_init, axis=1, keepdims=True)
130    expected_components = svd_flip(u=expected_components.T, v=None)[0].T
131    assert_allclose(model.components_, expected_components)
132
133
134def test_mini_batch_correct_shapes():
135    rng = np.random.RandomState(0)
136    X = rng.randn(12, 10)
137    pca = MiniBatchSparsePCA(n_components=8, max_iter=1, random_state=rng)
138    U = pca.fit_transform(X)
139    assert pca.components_.shape == (8, 10)
140    assert U.shape == (12, 8)
141    # test overcomplete decomposition
142    pca = MiniBatchSparsePCA(n_components=13, max_iter=1, random_state=rng)
143    U = pca.fit_transform(X)
144    assert pca.components_.shape == (13, 10)
145    assert U.shape == (12, 13)
146
147
148def test_scaling_fit_transform(global_random_seed):
149    alpha = 1
150    rng = np.random.RandomState(global_random_seed)
151    Y, _, _ = generate_toy_data(3, 1000, (8, 8), random_state=rng)
152    spca_lars = SparsePCA(n_components=3, method="lars", alpha=alpha, random_state=rng)
153    results_train = spca_lars.fit_transform(Y)
154    results_test = spca_lars.transform(Y[:10])
155    assert_allclose(results_train[0], results_test[0])
156
157
158def test_pca_vs_spca(global_random_seed):
159    rng = np.random.RandomState(global_random_seed)
160    Y, _, _ = generate_toy_data(3, 1000, (8, 8), random_state=rng)
161    Z, _, _ = generate_toy_data(3, 10, (8, 8), random_state=rng)
162    spca = SparsePCA(alpha=0, ridge_alpha=0, n_components=2, random_state=rng)
163    pca = PCA(n_components=2, random_state=rng)
164    pca.fit(Y)
165    spca.fit(Y)
166    results_test_pca = pca.transform(Z)
167    results_test_spca = spca.transform(Z)
168    assert_allclose(
169        np.abs(spca.components_.dot(pca.components_.T)), np.eye(2), atol=1e-4
170    )
171    results_test_pca *= np.sign(results_test_pca[0, :])
172    results_test_spca *= np.sign(results_test_spca[0, :])
173    assert_allclose(results_test_pca, results_test_spca, atol=1e-4)
174
175
176@pytest.mark.parametrize("SPCA", [SparsePCA, MiniBatchSparsePCA])
177@pytest.mark.parametrize("n_components", [None, 3])
178def test_spca_n_components_(SPCA, n_components):
179    rng = np.random.RandomState(0)
180    n_samples, n_features = 12, 10
181    X = rng.randn(n_samples, n_features)
182
183    model = SPCA(n_components=n_components).fit(X)
184
185    if n_components is not None:
186        assert model.n_components_ == n_components
187    else:
188        assert model.n_components_ == n_features
189
190
191@pytest.mark.parametrize("SPCA", (SparsePCA, MiniBatchSparsePCA))
192@pytest.mark.parametrize("method", ("lars", "cd"))
193@pytest.mark.parametrize(
194    "data_type, expected_type",
195    (
196        (np.float32, np.float32),
197        (np.float64, np.float64),
198        (np.int32, np.float64),
199        (np.int64, np.float64),
200    ),
201)
202def test_sparse_pca_dtype_match(SPCA, method, data_type, expected_type):
203    # Verify output matrix dtype
204    n_samples, n_features, n_components = 12, 10, 3
205    rng = np.random.RandomState(0)
206    input_array = rng.randn(n_samples, n_features).astype(data_type)
207    model = SPCA(n_components=n_components, method=method)
208    transformed = model.fit_transform(input_array)
209
210    assert transformed.dtype == expected_type
211    assert model.components_.dtype == expected_type
212
213
214@pytest.mark.parametrize("SPCA", (SparsePCA, MiniBatchSparsePCA))
215@pytest.mark.parametrize("method", ("lars", "cd"))
216def test_sparse_pca_numerical_consistency(SPCA, method, global_random_seed):
217    # Verify numericall consistentency among np.float32 and np.float64
218    n_samples, n_features, n_components = 20, 20, 5
219    input_array = make_low_rank_matrix(
220        n_samples=n_samples,
221        n_features=n_features,
222        effective_rank=n_components,
223        random_state=global_random_seed,
224    )
225
226    model_32 = SPCA(
227        n_components=n_components,
228        method=method,
229        random_state=global_random_seed,
230    )
231    transformed_32 = model_32.fit_transform(input_array.astype(np.float32))
232
233    model_64 = SPCA(
234        n_components=n_components,
235        method=method,
236        random_state=global_random_seed,
237    )
238    transformed_64 = model_64.fit_transform(input_array.astype(np.float64))
239    assert_allclose(transformed_64, transformed_32)
240    assert_allclose(model_64.components_, model_32.components_)
241
242
243@pytest.mark.parametrize("SPCA", [SparsePCA, MiniBatchSparsePCA])
244def test_spca_feature_names_out(SPCA):
245    """Check feature names out for *SparsePCA."""
246    rng = np.random.RandomState(0)
247    n_samples, n_features = 12, 10
248    X = rng.randn(n_samples, n_features)
249
250    model = SPCA(n_components=4).fit(X)
251    names = model.get_feature_names_out()
252
253    estimator_name = SPCA.__name__.lower()
254    assert_array_equal([f"{estimator_name}{i}" for i in range(4)], names)
255
256
257def test_spca_early_stopping(global_random_seed):
258    """Check that `tol` and `max_no_improvement` act as early stopping."""
259    rng = np.random.RandomState(global_random_seed)
260    n_samples, n_features = 50, 10
261    X = rng.randn(n_samples, n_features)
262
263    # vary the tolerance to force the early stopping of one of the model
264    model_early_stopped = MiniBatchSparsePCA(
265        max_iter=100, tol=0.5, random_state=global_random_seed
266    ).fit(X)
267    model_not_early_stopped = MiniBatchSparsePCA(
268        max_iter=100, tol=1e-3, random_state=global_random_seed
269    ).fit(X)
270    assert model_early_stopped.n_iter_ < model_not_early_stopped.n_iter_
271
272    # force the max number of no improvement to a large value to check that
273    # it does help to early stop
274    model_early_stopped = MiniBatchSparsePCA(
275        max_iter=100, tol=1e-6, max_no_improvement=2, random_state=global_random_seed
276    ).fit(X)
277    model_not_early_stopped = MiniBatchSparsePCA(
278        max_iter=100, tol=1e-6, max_no_improvement=100, random_state=global_random_seed
279    ).fit(X)
280    assert model_early_stopped.n_iter_ < model_not_early_stopped.n_iter_
281
282
283def test_equivalence_components_pca_spca(global_random_seed):
284    """Check the equivalence of the components found by PCA and SparsePCA.
285
286    Non-regression test for:
287    https://github.com/scikit-learn/scikit-learn/issues/23932
288    """
289    rng = np.random.RandomState(global_random_seed)
290    X = rng.randn(50, 4)
291
292    n_components = 2
293    pca = PCA(
294        n_components=n_components,
295        svd_solver="randomized",
296        random_state=0,
297    ).fit(X)
298    spca = SparsePCA(
299        n_components=n_components,
300        method="lars",
301        ridge_alpha=0,
302        alpha=0,
303        random_state=0,
304    ).fit(X)
305
306    assert_allclose(pca.components_, spca.components_)
307
308
309def test_sparse_pca_inverse_transform(global_random_seed):
310    """Check that `inverse_transform` in `SparsePCA` and `PCA` are similar."""
311    rng = np.random.RandomState(global_random_seed)
312    n_samples, n_features = 10, 5
313    X = rng.randn(n_samples, n_features)
314
315    n_components = 2
316    spca = SparsePCA(
317        n_components=n_components,
318        alpha=1e-12,
319        ridge_alpha=1e-12,
320        random_state=global_random_seed,
321    )
322    pca = PCA(n_components=n_components, random_state=global_random_seed)
323    X_trans_spca = spca.fit_transform(X)
324    X_trans_pca = pca.fit_transform(X)
325    assert_allclose(
326        spca.inverse_transform(X_trans_spca), pca.inverse_transform(X_trans_pca)
327    )
328
329
330@pytest.mark.parametrize("SPCA", [SparsePCA, MiniBatchSparsePCA])
331def test_transform_inverse_transform_round_trip(SPCA, global_random_seed):
332    """Check the `transform` and `inverse_transform` round trip with no loss of
333    information.
334    """
335    rng = np.random.RandomState(global_random_seed)
336    n_samples, n_features = 10, 5
337    X = rng.randn(n_samples, n_features)
338
339    n_components = n_features
340    spca = SPCA(
341        n_components=n_components,
342        alpha=1e-12,
343        ridge_alpha=1e-12,
344        random_state=global_random_seed,
345    )
346    X_trans_spca = spca.fit_transform(X)
347    assert_allclose(spca.inverse_transform(X_trans_spca), X)
348