Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_eigen/_svds.py: 11%
133 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import os
2import numpy as np
4from .arpack import _arpack # type: ignore[attr-defined]
5from . import eigsh
7from scipy._lib._util import check_random_state
8from scipy.sparse.linalg._interface import LinearOperator, aslinearoperator
9from scipy.sparse.linalg._eigen.lobpcg import lobpcg # type: ignore[no-redef]
10if os.environ.get("SCIPY_USE_PROPACK"):
11 from scipy.sparse.linalg._svdp import _svdp
12 HAS_PROPACK = True
13else:
14 HAS_PROPACK = False
15from scipy.linalg import svd
17arpack_int = _arpack.timing.nbx.dtype
18__all__ = ['svds']
21def _herm(x):
22 return x.T.conj()
25def _iv(A, k, ncv, tol, which, v0, maxiter,
26 return_singular, solver, random_state):
28 # input validation/standardization for `solver`
29 # out of order because it's needed for other parameters
30 solver = str(solver).lower()
31 solvers = {"arpack", "lobpcg", "propack"}
32 if solver not in solvers:
33 raise ValueError(f"solver must be one of {solvers}.")
35 # input validation/standardization for `A`
36 A = aslinearoperator(A) # this takes care of some input validation
37 if not (np.issubdtype(A.dtype, np.complexfloating)
38 or np.issubdtype(A.dtype, np.floating)):
39 message = "`A` must be of floating or complex floating data type."
40 raise ValueError(message)
41 if np.prod(A.shape) == 0:
42 message = "`A` must not be empty."
43 raise ValueError(message)
45 # input validation/standardization for `k`
46 kmax = min(A.shape) if solver == 'propack' else min(A.shape) - 1
47 if int(k) != k or not (0 < k <= kmax):
48 message = "`k` must be an integer satisfying `0 < k < min(A.shape)`."
49 raise ValueError(message)
50 k = int(k)
52 # input validation/standardization for `ncv`
53 if solver == "arpack" and ncv is not None:
54 if int(ncv) != ncv or not (k < ncv < min(A.shape)):
55 message = ("`ncv` must be an integer satisfying "
56 "`k < ncv < min(A.shape)`.")
57 raise ValueError(message)
58 ncv = int(ncv)
60 # input validation/standardization for `tol`
61 if tol < 0 or not np.isfinite(tol):
62 message = "`tol` must be a non-negative floating point value."
63 raise ValueError(message)
64 tol = float(tol)
66 # input validation/standardization for `which`
67 which = str(which).upper()
68 whichs = {'LM', 'SM'}
69 if which not in whichs:
70 raise ValueError(f"`which` must be in {whichs}.")
72 # input validation/standardization for `v0`
73 if v0 is not None:
74 v0 = np.atleast_1d(v0)
75 if not (np.issubdtype(v0.dtype, np.complexfloating)
76 or np.issubdtype(v0.dtype, np.floating)):
77 message = ("`v0` must be of floating or complex floating "
78 "data type.")
79 raise ValueError(message)
81 shape = (A.shape[0],) if solver == 'propack' else (min(A.shape),)
82 if v0.shape != shape:
83 message = f"`v0` must have shape {shape}."
84 raise ValueError(message)
86 # input validation/standardization for `maxiter`
87 if maxiter is not None and (int(maxiter) != maxiter or maxiter <= 0):
88 message = "`maxiter` must be a positive integer."
89 raise ValueError(message)
90 maxiter = int(maxiter) if maxiter is not None else maxiter
92 # input validation/standardization for `return_singular_vectors`
93 # not going to be flexible with this; too complicated for little gain
94 rs_options = {True, False, "vh", "u"}
95 if return_singular not in rs_options:
96 raise ValueError(f"`return_singular_vectors` must be in {rs_options}.")
98 random_state = check_random_state(random_state)
100 return (A, k, ncv, tol, which, v0, maxiter,
101 return_singular, solver, random_state)
104def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None,
105 maxiter=None, return_singular_vectors=True,
106 solver='arpack', random_state=None, options=None):
107 """
108 Partial singular value decomposition of a sparse matrix.
110 Compute the largest or smallest `k` singular values and corresponding
111 singular vectors of a sparse matrix `A`. The order in which the singular
112 values are returned is not guaranteed.
114 In the descriptions below, let ``M, N = A.shape``.
116 Parameters
117 ----------
118 A : ndarray, sparse matrix, or LinearOperator
119 Matrix to decompose of a floating point numeric dtype.
120 k : int, default: 6
121 Number of singular values and singular vectors to compute.
122 Must satisfy ``1 <= k <= kmax``, where ``kmax=min(M, N)`` for
123 ``solver='propack'`` and ``kmax=min(M, N) - 1`` otherwise.
124 ncv : int, optional
125 When ``solver='arpack'``, this is the number of Lanczos vectors
126 generated. See :ref:`'arpack' <sparse.linalg.svds-arpack>` for details.
127 When ``solver='lobpcg'`` or ``solver='propack'``, this parameter is
128 ignored.
129 tol : float, optional
130 Tolerance for singular values. Zero (default) means machine precision.
131 which : {'LM', 'SM'}
132 Which `k` singular values to find: either the largest magnitude ('LM')
133 or smallest magnitude ('SM') singular values.
134 v0 : ndarray, optional
135 The starting vector for iteration; see method-specific
136 documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
137 :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
138 :ref:`'propack' <sparse.linalg.svds-propack>` for details.
139 maxiter : int, optional
140 Maximum number of iterations; see method-specific
141 documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
142 :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
143 :ref:`'propack' <sparse.linalg.svds-propack>` for details.
144 return_singular_vectors : {True, False, "u", "vh"}
145 Singular values are always computed and returned; this parameter
146 controls the computation and return of singular vectors.
148 - ``True``: return singular vectors.
149 - ``False``: do not return singular vectors.
150 - ``"u"``: if ``M <= N``, compute only the left singular vectors and
151 return ``None`` for the right singular vectors. Otherwise, compute
152 all singular vectors.
153 - ``"vh"``: if ``M > N``, compute only the right singular vectors and
154 return ``None`` for the left singular vectors. Otherwise, compute
155 all singular vectors.
157 If ``solver='propack'``, the option is respected regardless of the
158 matrix shape.
160 solver : {'arpack', 'propack', 'lobpcg'}, optional
161 The solver used.
162 :ref:`'arpack' <sparse.linalg.svds-arpack>`,
163 :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`, and
164 :ref:`'propack' <sparse.linalg.svds-propack>` are supported.
165 Default: `'arpack'`.
166 random_state : {None, int, `numpy.random.Generator`,
167 `numpy.random.RandomState`}, optional
169 Pseudorandom number generator state used to generate resamples.
171 If `random_state` is ``None`` (or `np.random`), the
172 `numpy.random.RandomState` singleton is used.
173 If `random_state` is an int, a new ``RandomState`` instance is used,
174 seeded with `random_state`.
175 If `random_state` is already a ``Generator`` or ``RandomState``
176 instance then that instance is used.
177 options : dict, optional
178 A dictionary of solver-specific options. No solver-specific options
179 are currently supported; this parameter is reserved for future use.
181 Returns
182 -------
183 u : ndarray, shape=(M, k)
184 Unitary matrix having left singular vectors as columns.
185 s : ndarray, shape=(k,)
186 The singular values.
187 vh : ndarray, shape=(k, N)
188 Unitary matrix having right singular vectors as rows.
190 Notes
191 -----
192 This is a naive implementation using ARPACK or LOBPCG as an eigensolver
193 on the matrix ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on
194 which one is smaller size, followed by the Rayleigh-Ritz method
195 as postprocessing; see
196 Using the normal matrix, in Rayleigh-Ritz method, (2022, Nov. 19),
197 Wikipedia, https://w.wiki/4zms.
199 Alternatively, the PROPACK solver can be called. ``form="array"``
201 Choices of the input matrix ``A`` numeric dtype may be limited.
202 Only ``solver="lobpcg"`` supports all floating point dtypes
203 real: 'np.single', 'np.double', 'np.longdouble' and
204 complex: 'np.csingle', 'np.cdouble', 'np.clongdouble'.
205 The ``solver="arpack"`` supports only
206 'np.single', 'np.double', and 'np.cdouble'.
208 Examples
209 --------
210 Construct a matrix ``A`` from singular values and vectors.
212 >>> import numpy as np
213 >>> from scipy.stats import ortho_group
214 >>> from scipy.sparse.linalg import svds
215 >>> from scipy.sparse import csr_matrix
216 >>> rng = np.random.default_rng()
218 Construct a dense matrix ``A`` from singular values and vectors.
220 >>> orthogonal = ortho_group.rvs(10, random_state=rng)
221 >>> s = [1e-3, 1, 2, 3, 4] # non-zero singular values
222 >>> u = orthogonal[:, :5] # left singular vectors
223 >>> vT = orthogonal[:, 5:].T # right singular vectors
224 >>> A = u @ np.diag(s) @ vT
226 With only four singular values/vectors, the SVD approximates the original
227 matrix.
229 >>> u4, s4, vT4 = svds(A, k=4)
230 >>> A4 = u4 @ np.diag(s4) @ vT4
231 >>> np.allclose(A4, A, atol=1e-3)
232 True
234 With all five non-zero singular values/vectors, we can reproduce
235 the original matrix more accurately.
237 >>> u5, s5, vT5 = svds(A, k=5)
238 >>> A5 = u5 @ np.diag(s5) @ vT5
239 >>> np.allclose(A5, A)
240 True
242 The singular values match the expected singular values.
244 >>> np.allclose(s5, s)
245 True
247 Since the singular values are not close to each other in this example,
248 every singular vector matches as expected up to a difference in sign.
250 >>> (np.allclose(np.abs(u5), np.abs(u)) and
251 ... np.allclose(np.abs(vT5), np.abs(vT)))
252 True
254 The singular vectors are also orthogonal.
256 >>> (np.allclose(u5.T @ u5, np.eye(5)) and
257 ... np.allclose(vT5 @ vT5.T, np.eye(5)))
258 True
260 If there are (nearly) multiple singular values, the corresponding
261 individual singular vectors may be unstable, but the whole invariant
262 subspace containing all such singular vectors is computed accurately
263 as can be measured by angles between subspaces via 'subspace_angles'.
265 >>> from scipy.linalg import subspace_angles as s_a
266 >>> rng = np.random.default_rng()
267 >>> s = [1, 1 + 1e-6] # non-zero singular values
268 >>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
269 >>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
270 >>> vT = v.T
271 >>> A = u @ np.diag(s) @ vT
272 >>> A = A.astype(np.float32)
273 >>> u2, s2, vT2 = svds(A, k=2)
274 >>> np.allclose(s2, s)
275 True
277 The angles between the individual exact and computed singular vectors
278 are not so small.
280 >>> s_a(u2[:, :1], u[:, :1]) + s_a(u2[:, 1:], u[:, 1:]) > 1e-3
281 True
283 >>> (s_a(vT2[:1, :].T, vT[:1, :].T) +
284 ... s_a(vT2[1:, :].T, vT[1:, :].T)) > 1e-3
285 True
287 As opposed to the angles between the 2-dimensional invariant subspaces
288 that these vectors span, which are small for rights singular vectors
290 >>> s_a(u2, u).sum() < 1e-6
291 True
293 as well as for left singular vectors.
295 >>> s_a(vT2.T, vT.T).sum() < 1e-6
296 True
298 The next example follows that of 'sklearn.decomposition.TruncatedSVD'.
300 >>> rng = np.random.RandomState(0)
301 >>> X_dense = rng.random(size=(100, 100))
302 >>> X_dense[:, 2 * np.arange(50)] = 0
303 >>> X = csr_matrix(X_dense)
304 >>> _, singular_values, _ = svds(X, k=5)
305 >>> print(singular_values)
306 [ 4.3293... 4.4491... 4.5420... 4.5987... 35.2410...]
308 The function can be called without the transpose of the input matrix
309 ever explicitly constructed.
311 >>> from scipy.linalg import svd
312 >>> from scipy.sparse import rand
313 >>> from scipy.sparse.linalg import aslinearoperator
314 >>> rng = np.random.RandomState(0)
315 >>> G = rand(8, 9, density=0.5, random_state=rng)
316 >>> Glo = aslinearoperator(G)
317 >>> _, singular_values_svds, _ = svds(Glo, k=5)
318 >>> _, singular_values_svd, _ = svd(G.toarray())
319 >>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
320 True
322 The most memory efficient scenario is where neither
323 the original matrix, nor its transpose, is explicitly constructed.
324 Our example computes the smallest singular values and vectors
325 of 'LinearOperator' constructed from the numpy function 'np.diff' used
326 column-wise to be consistent with 'LinearOperator' operating on columns.
328 >>> from scipy.sparse.linalg import LinearOperator, aslinearoperator
329 >>> diff0 = lambda a: np.diff(a, axis=0)
331 Let us create the matrix from 'diff0' to be used for validation only.
333 >>> n = 5 # The dimension of the space.
334 >>> M_from_diff0 = diff0(np.eye(n))
335 >>> print(M_from_diff0.astype(int))
336 [[-1 1 0 0 0]
337 [ 0 -1 1 0 0]
338 [ 0 0 -1 1 0]
339 [ 0 0 0 -1 1]]
341 The matrix 'M_from_diff0' is bi-diagonal and could be alternatively
342 created directly by
344 >>> M = - np.eye(n - 1, n, dtype=int)
345 >>> np.fill_diagonal(M[:,1:], 1)
346 >>> np.allclose(M, M_from_diff0)
347 True
349 Its transpose
351 >>> print(M.T)
352 [[-1 0 0 0]
353 [ 1 -1 0 0]
354 [ 0 1 -1 0]
355 [ 0 0 1 -1]
356 [ 0 0 0 1]]
358 can be viewed as the incidence matrix; see
359 Incidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU,
360 of a linear graph with 5 vertices and 4 edges. The 5x5 normal matrix
361 'M.T @ M' thus is
363 >>> print(M.T @ M)
364 [[ 1 -1 0 0 0]
365 [-1 2 -1 0 0]
366 [ 0 -1 2 -1 0]
367 [ 0 0 -1 2 -1]
368 [ 0 0 0 -1 1]]
370 the graph Laplacian, while the actually used in 'svds' smaller size
371 4x4 normal matrix 'M @ M.T'
373 >>> print(M @ M.T)
374 [[ 2 -1 0 0]
375 [-1 2 -1 0]
376 [ 0 -1 2 -1]
377 [ 0 0 -1 2]]
379 is the so-called edge-based Laplacian; see
380 Symmetric Laplacian via the incidence matrix, in Laplacian matrix,
381 (2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.
383 The 'LinearOperator' setup needs the options 'rmatvec' and 'rmatmat'
384 of multiplication by the matrix transpose 'M.T', but we want to be
385 matrix-free to save memory, so knowing how 'M.T' looks like, we
386 manually construct the following function to be used in 'rmatmat=diff0t'.
388 >>> def diff0t(a):
389 ... if a.ndim == 1:
390 ... a = a[:,np.newaxis] # Turn 1D into 2D array
391 ... d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
392 ... d[0, :] = - a[0, :]
393 ... d[1:-1, :] = a[0:-1, :] - a[1:, :]
394 ... d[-1, :] = a[-1, :]
395 ... return d
397 We check that our function 'diff0t' for the matrix transpose is valid.
399 >>> np.allclose(M.T, diff0t(np.eye(n-1)))
400 True
402 Now we setup our matrix-free 'LinearOperator' called 'diff0_func_aslo'
403 and for validation the matrix-based 'diff0_matrix_aslo'.
405 >>> def diff0_func_aslo_def(n):
406 ... return LinearOperator(matvec=diff0,
407 ... matmat=diff0,
408 ... rmatvec=diff0t,
409 ... rmatmat=diff0t,
410 ... shape=(n - 1, n))
411 >>> diff0_func_aslo = diff0_func_aslo_def(n)
412 >>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)
414 And validate both the matrix and its transpose in 'LinearOperator'.
416 >>> np.allclose(diff0_func_aslo(np.eye(n)),
417 ... diff0_matrix_aslo(np.eye(n)))
418 True
419 >>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
420 ... diff0_matrix_aslo.T(np.eye(n-1)))
421 True
423 Having the 'LinearOperator' setup validated, we run the solver.
425 >>> n = 100
426 >>> diff0_func_aslo = diff0_func_aslo_def(n)
427 >>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')
429 The singular values squared and the singular vectors are known
430 explicitly; see
431 Pure Dirichlet boundary conditions, in
432 Eigenvalues and eigenvectors of the second derivative,
433 (2022, Nov. 19), Wikipedia, https://w.wiki/5YX6,
434 since 'diff' corresponds to first
435 derivative, and its smaller size n-1 x n-1 normal matrix
436 'M @ M.T' represent the discrete second derivative with the Dirichlet
437 boundary conditions. We use these analytic expressions for validation.
439 >>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
440 >>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
441 ... np.arange(1, 4)) / n)
442 >>> np.allclose(s, se, atol=1e-3)
443 True
444 >>> print(np.allclose(np.abs(u), np.abs(ue), atol=1e-6))
445 True
446 """
447 args = _iv(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors,
448 solver, random_state)
449 (A, k, ncv, tol, which, v0, maxiter,
450 return_singular_vectors, solver, random_state) = args
452 largest = (which == 'LM')
453 n, m = A.shape
455 if n >= m:
456 X_dot = A.matvec
457 X_matmat = A.matmat
458 XH_dot = A.rmatvec
459 XH_mat = A.rmatmat
460 transpose = False
461 else:
462 X_dot = A.rmatvec
463 X_matmat = A.rmatmat
464 XH_dot = A.matvec
465 XH_mat = A.matmat
466 transpose = True
468 dtype = getattr(A, 'dtype', None)
469 if dtype is None:
470 dtype = A.dot(np.zeros([m, 1])).dtype
472 def matvec_XH_X(x):
473 return XH_dot(X_dot(x))
475 def matmat_XH_X(x):
476 return XH_mat(X_matmat(x))
478 XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype,
479 matmat=matmat_XH_X,
480 shape=(min(A.shape), min(A.shape)))
482 # Get a low rank approximation of the implicitly defined gramian matrix.
483 # This is not a stable way to approach the problem.
484 if solver == 'lobpcg':
486 if k == 1 and v0 is not None:
487 X = np.reshape(v0, (-1, 1))
488 else:
489 X = random_state.standard_normal(size=(min(A.shape), k))
491 _, eigvec = lobpcg(XH_X, X, tol=tol ** 2, maxiter=maxiter,
492 largest=largest)
493 # lobpcg does not guarantee exactly orthonormal eigenvectors
494 # until after gh-16320 is merged
495 eigvec, _ = np.linalg.qr(eigvec)
497 elif solver == 'propack':
498 if not HAS_PROPACK:
499 raise ValueError("`solver='propack'` is opt-in due "
500 "to potential issues on Windows, "
501 "it can be enabled by setting the "
502 "`SCIPY_USE_PROPACK` environment "
503 "variable before importing scipy")
504 jobu = return_singular_vectors in {True, 'u'}
505 jobv = return_singular_vectors in {True, 'vh'}
506 irl_mode = (which == 'SM')
507 res = _svdp(A, k=k, tol=tol**2, which=which, maxiter=None,
508 compute_u=jobu, compute_v=jobv, irl_mode=irl_mode,
509 kmax=maxiter, v0=v0, random_state=random_state)
511 u, s, vh, _ = res # but we'll ignore bnd, the last output
513 # PROPACK order appears to be largest first. `svds` output order is not
514 # guaranteed, according to documentation, but for ARPACK and LOBPCG
515 # they actually are ordered smallest to largest, so reverse for
516 # consistency.
517 s = s[::-1]
518 u = u[:, ::-1]
519 vh = vh[::-1]
521 u = u if jobu else None
522 vh = vh if jobv else None
524 if return_singular_vectors:
525 return u, s, vh
526 else:
527 return s
529 elif solver == 'arpack' or solver is None:
530 if v0 is None:
531 v0 = random_state.standard_normal(size=(min(A.shape),))
532 _, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter,
533 ncv=ncv, which=which, v0=v0)
534 # arpack do not guarantee exactly orthonormal eigenvectors
535 # for clustered eigenvalues, especially in complex arithmetic
536 eigvec, _ = np.linalg.qr(eigvec)
538 # the eigenvectors eigvec must be orthonomal here; see gh-16712
539 Av = X_matmat(eigvec)
540 if not return_singular_vectors:
541 s = svd(Av, compute_uv=False, overwrite_a=True)
542 return s[::-1]
544 # compute the left singular vectors of X and update the right ones
545 # accordingly
546 u, s, vh = svd(Av, full_matrices=False, overwrite_a=True)
547 u = u[:, ::-1]
548 s = s[::-1]
549 vh = vh[::-1]
551 jobu = return_singular_vectors in {True, 'u'}
552 jobv = return_singular_vectors in {True, 'vh'}
554 if transpose:
555 u_tmp = eigvec @ _herm(vh) if jobu else None
556 vh = _herm(u) if jobv else None
557 u = u_tmp
558 else:
559 if not jobu:
560 u = None
561 vh = vh @ _herm(eigvec) if jobv else None
563 return u, s, vh