-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathlinalg.py
More file actions
83 lines (69 loc) · 2.74 KB
/
linalg.py
File metadata and controls
83 lines (69 loc) · 2.74 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
"""
The linear algebra module. This module contains functions that are often
used in hyperalignment algorithms. Specifically, the robustness of
singular value decomposition (SVD) is enhanced in ``safe_svd`` to avoid
occasional crashes when the operation is performed many times (e.g., in a
searchlight algorithm), and ``svd_pca`` performs PCA based on
``safe_svd``.
"""
import numpy as np
from scipy.linalg import svd, LinAlgError
__all__ = ['safe_svd', 'svd_pca']
def safe_svd(X, remove_mean=True):
"""
Singular value decomposition without occasional LinAlgError crashes.
The default ``lapack_driver`` of ``scipy.linalg.svd`` is ``'gesdd'``,
which occassionaly crashes even if the input matrix is not singular.
This function automatically handles the ``LinAlgError`` when it's
raised and switches to the ``'gesvd'`` driver in this case.
The input matrix ``X`` is factorized as ``U @ np.diag(s) @ Vt``.
Parameters
----------
X : ndarray of shape (M, N)
The matrix to be decomposed in NumPy array format.
remove_mean : bool, default=True
Whether to subtract the mean of each column before the actual SVD
(True) or not (False). Setting `remove_mean=True` is helpful when
the SVD is used to perform PCA.
Returns
-------
U : ndarray of shape (M, K)
Unitary matrix.
s : ndarray of shape (K,)
The singular values.
Vt : ndarray of shape (K, N)
Unitary matrix.
"""
if remove_mean:
X = X - X.mean(axis=0, keepdims=True)
try:
U, s, Vt = svd(X, full_matrices=False)
except LinAlgError:
U, s, Vt = svd(X, full_matrices=False, lapack_driver='gesvd')
return U, s, Vt
def svd_pca(X, remove_mean=True):
"""
Principal component analysis (PCA) based on SVD.
This function performs a rotation and returns the transformed data in
PC space. Therefore, its behavior is similar to the ``fit_transform``
method of ``sklearn.decomposition.PCA``.
It does not throw away any PCs, and therefore there is no
dimensionality reduction in the PC space. However, the number of PCs
might be less than the number of features in ``X``, depending on the
rank of ``X``.
Parameters
----------
X : ndarray of shape (M, N)
The data matrix to be transformed into PC space.
remove_mean : bool, default=True
Whether to subtract the mean of each column before the SVD (True)
or not (False). This parameter should be set to True unless the
columns already have zero mean.
Returns
-------
X_new : ndarray of shape (M, K)
The transformed data matrix in PC space.
"""
U, s, Vt = safe_svd(X, remove_mean=remove_mean)
X_new = U * s[np.newaxis]
return X_new