Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_eigen/_svds.py: 11%
132 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +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.
201 Choices of the input matrix `A` numeric dtype may be limited.
202 Only ``solver="lobpcg"`` supports all floating point dtypes
203 real: 'np.float32', 'np.float64', 'np.longdouble' and
204 complex: 'np.complex64', 'np.complex128', 'np.clongdouble'.
205 The ``solver="arpack"`` supports only
206 'np.float32', 'np.float64', and 'np.complex128'.
208 Examples
209 --------
210 Construct a matrix `A` from singular values and vectors.
212 >>> import numpy as np
213 >>> from scipy import sparse, linalg, stats
214 >>> from scipy.sparse.linalg import svds, aslinearoperator, LinearOperator
216 Construct a dense matrix `A` from singular values and vectors.
218 >>> rng = np.random.default_rng(258265244568965474821194062361901728911)
219 >>> orthogonal = stats.ortho_group.rvs(10, random_state=rng)
220 >>> s = [1e-3, 1, 2, 3, 4] # non-zero singular values
221 >>> u = orthogonal[:, :5] # left singular vectors
222 >>> vT = orthogonal[:, 5:].T # right singular vectors
223 >>> A = u @ np.diag(s) @ vT
225 With only four singular values/vectors, the SVD approximates the original
226 matrix.
228 >>> u4, s4, vT4 = svds(A, k=4)
229 >>> A4 = u4 @ np.diag(s4) @ vT4
230 >>> np.allclose(A4, A, atol=1e-3)
231 True
233 With all five non-zero singular values/vectors, we can reproduce
234 the original matrix more accurately.
236 >>> u5, s5, vT5 = svds(A, k=5)
237 >>> A5 = u5 @ np.diag(s5) @ vT5
238 >>> np.allclose(A5, A)
239 True
241 The singular values match the expected singular values.
243 >>> np.allclose(s5, s)
244 True
246 Since the singular values are not close to each other in this example,
247 every singular vector matches as expected up to a difference in sign.
249 >>> (np.allclose(np.abs(u5), np.abs(u)) and
250 ... np.allclose(np.abs(vT5), np.abs(vT)))
251 True
253 The singular vectors are also orthogonal.
255 >>> (np.allclose(u5.T @ u5, np.eye(5)) and
256 ... np.allclose(vT5 @ vT5.T, np.eye(5)))
257 True
259 If there are (nearly) multiple singular values, the corresponding
260 individual singular vectors may be unstable, but the whole invariant
261 subspace containing all such singular vectors is computed accurately
262 as can be measured by angles between subspaces via 'subspace_angles'.
264 >>> rng = np.random.default_rng(178686584221410808734965903901790843963)
265 >>> s = [1, 1 + 1e-6] # non-zero singular values
266 >>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
267 >>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
268 >>> vT = v.T
269 >>> A = u @ np.diag(s) @ vT
270 >>> A = A.astype(np.float32)
271 >>> u2, s2, vT2 = svds(A, k=2, random_state=rng)
272 >>> np.allclose(s2, s)
273 True
275 The angles between the individual exact and computed singular vectors
276 may not be so small. To check use:
278 >>> (linalg.subspace_angles(u2[:, :1], u[:, :1]) +
279 ... linalg.subspace_angles(u2[:, 1:], u[:, 1:]))
280 array([0.06562513]) # may vary
281 >>> (linalg.subspace_angles(vT2[:1, :].T, vT[:1, :].T) +
282 ... linalg.subspace_angles(vT2[1:, :].T, vT[1:, :].T))
283 array([0.06562507]) # may vary
285 As opposed to the angles between the 2-dimensional invariant subspaces
286 that these vectors span, which are small for rights singular vectors
288 >>> linalg.subspace_angles(u2, u).sum() < 1e-6
289 True
291 as well as for left singular vectors.
293 >>> linalg.subspace_angles(vT2.T, vT.T).sum() < 1e-6
294 True
296 The next example follows that of 'sklearn.decomposition.TruncatedSVD'.
298 >>> rng = np.random.RandomState(0)
299 >>> X_dense = rng.random(size=(100, 100))
300 >>> X_dense[:, 2 * np.arange(50)] = 0
301 >>> X = sparse.csr_matrix(X_dense)
302 >>> _, singular_values, _ = svds(X, k=5, random_state=rng)
303 >>> print(singular_values)
304 [ 4.3293... 4.4491... 4.5420... 4.5987... 35.2410...]
306 The function can be called without the transpose of the input matrix
307 ever explicitly constructed.
309 >>> rng = np.random.default_rng(102524723947864966825913730119128190974)
310 >>> G = sparse.rand(8, 9, density=0.5, random_state=rng)
311 >>> Glo = aslinearoperator(G)
312 >>> _, singular_values_svds, _ = svds(Glo, k=5, random_state=rng)
313 >>> _, singular_values_svd, _ = linalg.svd(G.toarray())
314 >>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
315 True
317 The most memory efficient scenario is where neither
318 the original matrix, nor its transpose, is explicitly constructed.
319 Our example computes the smallest singular values and vectors
320 of 'LinearOperator' constructed from the numpy function 'np.diff' used
321 column-wise to be consistent with 'LinearOperator' operating on columns.
323 >>> diff0 = lambda a: np.diff(a, axis=0)
325 Let us create the matrix from 'diff0' to be used for validation only.
327 >>> n = 5 # The dimension of the space.
328 >>> M_from_diff0 = diff0(np.eye(n))
329 >>> print(M_from_diff0.astype(int))
330 [[-1 1 0 0 0]
331 [ 0 -1 1 0 0]
332 [ 0 0 -1 1 0]
333 [ 0 0 0 -1 1]]
335 The matrix 'M_from_diff0' is bi-diagonal and could be alternatively
336 created directly by
338 >>> M = - np.eye(n - 1, n, dtype=int)
339 >>> np.fill_diagonal(M[:,1:], 1)
340 >>> np.allclose(M, M_from_diff0)
341 True
343 Its transpose
345 >>> print(M.T)
346 [[-1 0 0 0]
347 [ 1 -1 0 0]
348 [ 0 1 -1 0]
349 [ 0 0 1 -1]
350 [ 0 0 0 1]]
352 can be viewed as the incidence matrix; see
353 Incidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU,
354 of a linear graph with 5 vertices and 4 edges. The 5x5 normal matrix
355 ``M.T @ M`` thus is
357 >>> print(M.T @ M)
358 [[ 1 -1 0 0 0]
359 [-1 2 -1 0 0]
360 [ 0 -1 2 -1 0]
361 [ 0 0 -1 2 -1]
362 [ 0 0 0 -1 1]]
364 the graph Laplacian, while the actually used in 'svds' smaller size
365 4x4 normal matrix ``M @ M.T``
367 >>> print(M @ M.T)
368 [[ 2 -1 0 0]
369 [-1 2 -1 0]
370 [ 0 -1 2 -1]
371 [ 0 0 -1 2]]
373 is the so-called edge-based Laplacian; see
374 Symmetric Laplacian via the incidence matrix, in Laplacian matrix,
375 (2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.
377 The 'LinearOperator' setup needs the options 'rmatvec' and 'rmatmat'
378 of multiplication by the matrix transpose ``M.T``, but we want to be
379 matrix-free to save memory, so knowing how ``M.T`` looks like, we
380 manually construct the following function to be
381 used in ``rmatmat=diff0t``.
383 >>> def diff0t(a):
384 ... if a.ndim == 1:
385 ... a = a[:,np.newaxis] # Turn 1D into 2D array
386 ... d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
387 ... d[0, :] = - a[0, :]
388 ... d[1:-1, :] = a[0:-1, :] - a[1:, :]
389 ... d[-1, :] = a[-1, :]
390 ... return d
392 We check that our function 'diff0t' for the matrix transpose is valid.
394 >>> np.allclose(M.T, diff0t(np.eye(n-1)))
395 True
397 Now we setup our matrix-free 'LinearOperator' called 'diff0_func_aslo'
398 and for validation the matrix-based 'diff0_matrix_aslo'.
400 >>> def diff0_func_aslo_def(n):
401 ... return LinearOperator(matvec=diff0,
402 ... matmat=diff0,
403 ... rmatvec=diff0t,
404 ... rmatmat=diff0t,
405 ... shape=(n - 1, n))
406 >>> diff0_func_aslo = diff0_func_aslo_def(n)
407 >>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)
409 And validate both the matrix and its transpose in 'LinearOperator'.
411 >>> np.allclose(diff0_func_aslo(np.eye(n)),
412 ... diff0_matrix_aslo(np.eye(n)))
413 True
414 >>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
415 ... diff0_matrix_aslo.T(np.eye(n-1)))
416 True
418 Having the 'LinearOperator' setup validated, we run the solver.
420 >>> n = 100
421 >>> diff0_func_aslo = diff0_func_aslo_def(n)
422 >>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')
424 The singular values squared and the singular vectors are known
425 explicitly; see
426 Pure Dirichlet boundary conditions, in
427 Eigenvalues and eigenvectors of the second derivative,
428 (2022, Nov. 19), Wikipedia, https://w.wiki/5YX6,
429 since 'diff' corresponds to first
430 derivative, and its smaller size n-1 x n-1 normal matrix
431 ``M @ M.T`` represent the discrete second derivative with the Dirichlet
432 boundary conditions. We use these analytic expressions for validation.
434 >>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
435 >>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
436 ... np.arange(1, 4)) / n)
437 >>> np.allclose(s, se, atol=1e-3)
438 True
439 >>> print(np.allclose(np.abs(u), np.abs(ue), atol=1e-6))
440 True
442 """
443 args = _iv(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors,
444 solver, random_state)
445 (A, k, ncv, tol, which, v0, maxiter,
446 return_singular_vectors, solver, random_state) = args
448 largest = (which == 'LM')
449 n, m = A.shape
451 if n >= m:
452 X_dot = A.matvec
453 X_matmat = A.matmat
454 XH_dot = A.rmatvec
455 XH_mat = A.rmatmat
456 transpose = False
457 else:
458 X_dot = A.rmatvec
459 X_matmat = A.rmatmat
460 XH_dot = A.matvec
461 XH_mat = A.matmat
462 transpose = True
464 dtype = getattr(A, 'dtype', None)
465 if dtype is None:
466 dtype = A.dot(np.zeros([m, 1])).dtype
468 def matvec_XH_X(x):
469 return XH_dot(X_dot(x))
471 def matmat_XH_X(x):
472 return XH_mat(X_matmat(x))
474 XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype,
475 matmat=matmat_XH_X,
476 shape=(min(A.shape), min(A.shape)))
478 # Get a low rank approximation of the implicitly defined gramian matrix.
479 # This is not a stable way to approach the problem.
480 if solver == 'lobpcg':
482 if k == 1 and v0 is not None:
483 X = np.reshape(v0, (-1, 1))
484 else:
485 X = random_state.standard_normal(size=(min(A.shape), k))
487 _, eigvec = lobpcg(XH_X, X, tol=tol ** 2, maxiter=maxiter,
488 largest=largest)
490 elif solver == 'propack':
491 if not HAS_PROPACK:
492 raise ValueError("`solver='propack'` is opt-in due "
493 "to potential issues on Windows, "
494 "it can be enabled by setting the "
495 "`SCIPY_USE_PROPACK` environment "
496 "variable before importing scipy")
497 jobu = return_singular_vectors in {True, 'u'}
498 jobv = return_singular_vectors in {True, 'vh'}
499 irl_mode = (which == 'SM')
500 res = _svdp(A, k=k, tol=tol**2, which=which, maxiter=None,
501 compute_u=jobu, compute_v=jobv, irl_mode=irl_mode,
502 kmax=maxiter, v0=v0, random_state=random_state)
504 u, s, vh, _ = res # but we'll ignore bnd, the last output
506 # PROPACK order appears to be largest first. `svds` output order is not
507 # guaranteed, according to documentation, but for ARPACK and LOBPCG
508 # they actually are ordered smallest to largest, so reverse for
509 # consistency.
510 s = s[::-1]
511 u = u[:, ::-1]
512 vh = vh[::-1]
514 u = u if jobu else None
515 vh = vh if jobv else None
517 if return_singular_vectors:
518 return u, s, vh
519 else:
520 return s
522 elif solver == 'arpack' or solver is None:
523 if v0 is None:
524 v0 = random_state.standard_normal(size=(min(A.shape),))
525 _, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter,
526 ncv=ncv, which=which, v0=v0)
527 # arpack do not guarantee exactly orthonormal eigenvectors
528 # for clustered eigenvalues, especially in complex arithmetic
529 eigvec, _ = np.linalg.qr(eigvec)
531 # the eigenvectors eigvec must be orthonomal here; see gh-16712
532 Av = X_matmat(eigvec)
533 if not return_singular_vectors:
534 s = svd(Av, compute_uv=False, overwrite_a=True)
535 return s[::-1]
537 # compute the left singular vectors of X and update the right ones
538 # accordingly
539 u, s, vh = svd(Av, full_matrices=False, overwrite_a=True)
540 u = u[:, ::-1]
541 s = s[::-1]
542 vh = vh[::-1]
544 jobu = return_singular_vectors in {True, 'u'}
545 jobv = return_singular_vectors in {True, 'vh'}
547 if transpose:
548 u_tmp = eigvec @ _herm(vh) if jobu else None
549 vh = _herm(u) if jobv else None
550 u = u_tmp
551 else:
552 if not jobu:
553 u = None
554 vh = vh @ _herm(eigvec) if jobv else None
556 return u, s, vh