update
This commit is contained in:
@@ -1,22 +0,0 @@
|
||||
"""
|
||||
Sparse Eigenvalue Solvers
|
||||
-------------------------
|
||||
|
||||
The submodules of sparse.linalg._eigen:
|
||||
1. lobpcg: Locally Optimal Block Preconditioned Conjugate Gradient Method
|
||||
|
||||
"""
|
||||
from .arpack import *
|
||||
from .lobpcg import *
|
||||
from ._svds import svds
|
||||
|
||||
from . import arpack
|
||||
|
||||
__all__ = [
|
||||
'ArpackError', 'ArpackNoConvergence',
|
||||
'eigs', 'eigsh', 'lobpcg', 'svds'
|
||||
]
|
||||
|
||||
from scipy._lib._testutils import PytestTester
|
||||
test = PytestTester(__name__)
|
||||
del PytestTester
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -1,563 +0,0 @@
|
||||
import os
|
||||
import numpy as np
|
||||
|
||||
from .arpack import _arpack # type: ignore[attr-defined]
|
||||
from . import eigsh
|
||||
|
||||
from scipy._lib._util import check_random_state
|
||||
from scipy.sparse.linalg._interface import LinearOperator, aslinearoperator
|
||||
from scipy.sparse.linalg._eigen.lobpcg import lobpcg # type: ignore[no-redef]
|
||||
if os.environ.get("SCIPY_USE_PROPACK"):
|
||||
from scipy.sparse.linalg._svdp import _svdp
|
||||
HAS_PROPACK = True
|
||||
else:
|
||||
HAS_PROPACK = False
|
||||
from scipy.linalg import svd
|
||||
|
||||
arpack_int = _arpack.timing.nbx.dtype
|
||||
__all__ = ['svds']
|
||||
|
||||
|
||||
def _herm(x):
|
||||
return x.T.conj()
|
||||
|
||||
|
||||
def _iv(A, k, ncv, tol, which, v0, maxiter,
|
||||
return_singular, solver, random_state):
|
||||
|
||||
# input validation/standardization for `solver`
|
||||
# out of order because it's needed for other parameters
|
||||
solver = str(solver).lower()
|
||||
solvers = {"arpack", "lobpcg", "propack"}
|
||||
if solver not in solvers:
|
||||
raise ValueError(f"solver must be one of {solvers}.")
|
||||
|
||||
# input validation/standardization for `A`
|
||||
A = aslinearoperator(A) # this takes care of some input validation
|
||||
if not (np.issubdtype(A.dtype, np.complexfloating)
|
||||
or np.issubdtype(A.dtype, np.floating)):
|
||||
message = "`A` must be of floating or complex floating data type."
|
||||
raise ValueError(message)
|
||||
if np.prod(A.shape) == 0:
|
||||
message = "`A` must not be empty."
|
||||
raise ValueError(message)
|
||||
|
||||
# input validation/standardization for `k`
|
||||
kmax = min(A.shape) if solver == 'propack' else min(A.shape) - 1
|
||||
if int(k) != k or not (0 < k <= kmax):
|
||||
message = "`k` must be an integer satisfying `0 < k < min(A.shape)`."
|
||||
raise ValueError(message)
|
||||
k = int(k)
|
||||
|
||||
# input validation/standardization for `ncv`
|
||||
if solver == "arpack" and ncv is not None:
|
||||
if int(ncv) != ncv or not (k < ncv < min(A.shape)):
|
||||
message = ("`ncv` must be an integer satisfying "
|
||||
"`k < ncv < min(A.shape)`.")
|
||||
raise ValueError(message)
|
||||
ncv = int(ncv)
|
||||
|
||||
# input validation/standardization for `tol`
|
||||
if tol < 0 or not np.isfinite(tol):
|
||||
message = "`tol` must be a non-negative floating point value."
|
||||
raise ValueError(message)
|
||||
tol = float(tol)
|
||||
|
||||
# input validation/standardization for `which`
|
||||
which = str(which).upper()
|
||||
whichs = {'LM', 'SM'}
|
||||
if which not in whichs:
|
||||
raise ValueError(f"`which` must be in {whichs}.")
|
||||
|
||||
# input validation/standardization for `v0`
|
||||
if v0 is not None:
|
||||
v0 = np.atleast_1d(v0)
|
||||
if not (np.issubdtype(v0.dtype, np.complexfloating)
|
||||
or np.issubdtype(v0.dtype, np.floating)):
|
||||
message = ("`v0` must be of floating or complex floating "
|
||||
"data type.")
|
||||
raise ValueError(message)
|
||||
|
||||
shape = (A.shape[0],) if solver == 'propack' else (min(A.shape),)
|
||||
if v0.shape != shape:
|
||||
message = f"`v0` must have shape {shape}."
|
||||
raise ValueError(message)
|
||||
|
||||
# input validation/standardization for `maxiter`
|
||||
if maxiter is not None and (int(maxiter) != maxiter or maxiter <= 0):
|
||||
message = "`maxiter` must be a positive integer."
|
||||
raise ValueError(message)
|
||||
maxiter = int(maxiter) if maxiter is not None else maxiter
|
||||
|
||||
# input validation/standardization for `return_singular_vectors`
|
||||
# not going to be flexible with this; too complicated for little gain
|
||||
rs_options = {True, False, "vh", "u"}
|
||||
if return_singular not in rs_options:
|
||||
raise ValueError(f"`return_singular_vectors` must be in {rs_options}.")
|
||||
|
||||
random_state = check_random_state(random_state)
|
||||
|
||||
return (A, k, ncv, tol, which, v0, maxiter,
|
||||
return_singular, solver, random_state)
|
||||
|
||||
|
||||
def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None,
|
||||
maxiter=None, return_singular_vectors=True,
|
||||
solver='arpack', random_state=None, options=None):
|
||||
"""
|
||||
Partial singular value decomposition of a sparse matrix.
|
||||
|
||||
Compute the largest or smallest `k` singular values and corresponding
|
||||
singular vectors of a sparse matrix `A`. The order in which the singular
|
||||
values are returned is not guaranteed.
|
||||
|
||||
In the descriptions below, let ``M, N = A.shape``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : ndarray, sparse matrix, or LinearOperator
|
||||
Matrix to decompose of a floating point numeric dtype.
|
||||
k : int, default: 6
|
||||
Number of singular values and singular vectors to compute.
|
||||
Must satisfy ``1 <= k <= kmax``, where ``kmax=min(M, N)`` for
|
||||
``solver='propack'`` and ``kmax=min(M, N) - 1`` otherwise.
|
||||
ncv : int, optional
|
||||
When ``solver='arpack'``, this is the number of Lanczos vectors
|
||||
generated. See :ref:`'arpack' <sparse.linalg.svds-arpack>` for details.
|
||||
When ``solver='lobpcg'`` or ``solver='propack'``, this parameter is
|
||||
ignored.
|
||||
tol : float, optional
|
||||
Tolerance for singular values. Zero (default) means machine precision.
|
||||
which : {'LM', 'SM'}
|
||||
Which `k` singular values to find: either the largest magnitude ('LM')
|
||||
or smallest magnitude ('SM') singular values.
|
||||
v0 : ndarray, optional
|
||||
The starting vector for iteration; see method-specific
|
||||
documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
|
||||
:ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
|
||||
:ref:`'propack' <sparse.linalg.svds-propack>` for details.
|
||||
maxiter : int, optional
|
||||
Maximum number of iterations; see method-specific
|
||||
documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
|
||||
:ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
|
||||
:ref:`'propack' <sparse.linalg.svds-propack>` for details.
|
||||
return_singular_vectors : {True, False, "u", "vh"}
|
||||
Singular values are always computed and returned; this parameter
|
||||
controls the computation and return of singular vectors.
|
||||
|
||||
- ``True``: return singular vectors.
|
||||
- ``False``: do not return singular vectors.
|
||||
- ``"u"``: if ``M <= N``, compute only the left singular vectors and
|
||||
return ``None`` for the right singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
- ``"vh"``: if ``M > N``, compute only the right singular vectors and
|
||||
return ``None`` for the left singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
|
||||
If ``solver='propack'``, the option is respected regardless of the
|
||||
matrix shape.
|
||||
|
||||
solver : {'arpack', 'propack', 'lobpcg'}, optional
|
||||
The solver used.
|
||||
:ref:`'arpack' <sparse.linalg.svds-arpack>`,
|
||||
:ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`, and
|
||||
:ref:`'propack' <sparse.linalg.svds-propack>` are supported.
|
||||
Default: `'arpack'`.
|
||||
random_state : {None, int, `numpy.random.Generator`,
|
||||
`numpy.random.RandomState`}, optional
|
||||
|
||||
Pseudorandom number generator state used to generate resamples.
|
||||
|
||||
If `random_state` is ``None`` (or `np.random`), the
|
||||
`numpy.random.RandomState` singleton is used.
|
||||
If `random_state` is an int, a new ``RandomState`` instance is used,
|
||||
seeded with `random_state`.
|
||||
If `random_state` is already a ``Generator`` or ``RandomState``
|
||||
instance then that instance is used.
|
||||
options : dict, optional
|
||||
A dictionary of solver-specific options. No solver-specific options
|
||||
are currently supported; this parameter is reserved for future use.
|
||||
|
||||
Returns
|
||||
-------
|
||||
u : ndarray, shape=(M, k)
|
||||
Unitary matrix having left singular vectors as columns.
|
||||
s : ndarray, shape=(k,)
|
||||
The singular values.
|
||||
vh : ndarray, shape=(k, N)
|
||||
Unitary matrix having right singular vectors as rows.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This is a naive implementation using ARPACK or LOBPCG as an eigensolver
|
||||
on the matrix ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on
|
||||
which one is smaller size, followed by the Rayleigh-Ritz method
|
||||
as postprocessing; see
|
||||
Using the normal matrix, in Rayleigh-Ritz method, (2022, Nov. 19),
|
||||
Wikipedia, https://w.wiki/4zms.
|
||||
|
||||
Alternatively, the PROPACK solver can be called. ``form="array"``
|
||||
|
||||
Choices of the input matrix ``A`` numeric dtype may be limited.
|
||||
Only ``solver="lobpcg"`` supports all floating point dtypes
|
||||
real: 'np.single', 'np.double', 'np.longdouble' and
|
||||
complex: 'np.csingle', 'np.cdouble', 'np.clongdouble'.
|
||||
The ``solver="arpack"`` supports only
|
||||
'np.single', 'np.double', and 'np.cdouble'.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Construct a matrix ``A`` from singular values and vectors.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats import ortho_group
|
||||
>>> from scipy.sparse.linalg import svds
|
||||
>>> from scipy.sparse import csr_matrix
|
||||
>>> rng = np.random.default_rng()
|
||||
|
||||
Construct a dense matrix ``A`` from singular values and vectors.
|
||||
|
||||
>>> orthogonal = ortho_group.rvs(10, random_state=rng)
|
||||
>>> s = [1e-3, 1, 2, 3, 4] # non-zero singular values
|
||||
>>> u = orthogonal[:, :5] # left singular vectors
|
||||
>>> vT = orthogonal[:, 5:].T # right singular vectors
|
||||
>>> A = u @ np.diag(s) @ vT
|
||||
|
||||
With only four singular values/vectors, the SVD approximates the original
|
||||
matrix.
|
||||
|
||||
>>> u4, s4, vT4 = svds(A, k=4)
|
||||
>>> A4 = u4 @ np.diag(s4) @ vT4
|
||||
>>> np.allclose(A4, A, atol=1e-3)
|
||||
True
|
||||
|
||||
With all five non-zero singular values/vectors, we can reproduce
|
||||
the original matrix more accurately.
|
||||
|
||||
>>> u5, s5, vT5 = svds(A, k=5)
|
||||
>>> A5 = u5 @ np.diag(s5) @ vT5
|
||||
>>> np.allclose(A5, A)
|
||||
True
|
||||
|
||||
The singular values match the expected singular values.
|
||||
|
||||
>>> np.allclose(s5, s)
|
||||
True
|
||||
|
||||
Since the singular values are not close to each other in this example,
|
||||
every singular vector matches as expected up to a difference in sign.
|
||||
|
||||
>>> (np.allclose(np.abs(u5), np.abs(u)) and
|
||||
... np.allclose(np.abs(vT5), np.abs(vT)))
|
||||
True
|
||||
|
||||
The singular vectors are also orthogonal.
|
||||
|
||||
>>> (np.allclose(u5.T @ u5, np.eye(5)) and
|
||||
... np.allclose(vT5 @ vT5.T, np.eye(5)))
|
||||
True
|
||||
|
||||
If there are (nearly) multiple singular values, the corresponding
|
||||
individual singular vectors may be unstable, but the whole invariant
|
||||
subspace containing all such singular vectors is computed accurately
|
||||
as can be measured by angles between subspaces via 'subspace_angles'.
|
||||
|
||||
>>> from scipy.linalg import subspace_angles as s_a
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> s = [1, 1 + 1e-6] # non-zero singular values
|
||||
>>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
|
||||
>>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
|
||||
>>> vT = v.T
|
||||
>>> A = u @ np.diag(s) @ vT
|
||||
>>> A = A.astype(np.float32)
|
||||
>>> u2, s2, vT2 = svds(A, k=2)
|
||||
>>> np.allclose(s2, s)
|
||||
True
|
||||
|
||||
The angles between the individual exact and computed singular vectors
|
||||
are not so small.
|
||||
|
||||
>>> s_a(u2[:, :1], u[:, :1]) + s_a(u2[:, 1:], u[:, 1:]) > 1e-3
|
||||
True
|
||||
|
||||
>>> (s_a(vT2[:1, :].T, vT[:1, :].T) +
|
||||
... s_a(vT2[1:, :].T, vT[1:, :].T)) > 1e-3
|
||||
True
|
||||
|
||||
As opposed to the angles between the 2-dimensional invariant subspaces
|
||||
that these vectors span, which are small for rights singular vectors
|
||||
|
||||
>>> s_a(u2, u).sum() < 1e-6
|
||||
True
|
||||
|
||||
as well as for left singular vectors.
|
||||
|
||||
>>> s_a(vT2.T, vT.T).sum() < 1e-6
|
||||
True
|
||||
|
||||
The next example follows that of 'sklearn.decomposition.TruncatedSVD'.
|
||||
|
||||
>>> rng = np.random.RandomState(0)
|
||||
>>> X_dense = rng.random(size=(100, 100))
|
||||
>>> X_dense[:, 2 * np.arange(50)] = 0
|
||||
>>> X = csr_matrix(X_dense)
|
||||
>>> _, singular_values, _ = svds(X, k=5)
|
||||
>>> print(singular_values)
|
||||
[ 4.3293... 4.4491... 4.5420... 4.5987... 35.2410...]
|
||||
|
||||
The function can be called without the transpose of the input matrix
|
||||
ever explicitly constructed.
|
||||
|
||||
>>> from scipy.linalg import svd
|
||||
>>> from scipy.sparse import rand
|
||||
>>> from scipy.sparse.linalg import aslinearoperator
|
||||
>>> rng = np.random.RandomState(0)
|
||||
>>> G = rand(8, 9, density=0.5, random_state=rng)
|
||||
>>> Glo = aslinearoperator(G)
|
||||
>>> _, singular_values_svds, _ = svds(Glo, k=5)
|
||||
>>> _, singular_values_svd, _ = svd(G.toarray())
|
||||
>>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
|
||||
True
|
||||
|
||||
The most memory efficient scenario is where neither
|
||||
the original matrix, nor its transpose, is explicitly constructed.
|
||||
Our example computes the smallest singular values and vectors
|
||||
of 'LinearOperator' constructed from the numpy function 'np.diff' used
|
||||
column-wise to be consistent with 'LinearOperator' operating on columns.
|
||||
|
||||
>>> from scipy.sparse.linalg import LinearOperator, aslinearoperator
|
||||
>>> diff0 = lambda a: np.diff(a, axis=0)
|
||||
|
||||
Let us create the matrix from 'diff0' to be used for validation only.
|
||||
|
||||
>>> n = 5 # The dimension of the space.
|
||||
>>> M_from_diff0 = diff0(np.eye(n))
|
||||
>>> print(M_from_diff0.astype(int))
|
||||
[[-1 1 0 0 0]
|
||||
[ 0 -1 1 0 0]
|
||||
[ 0 0 -1 1 0]
|
||||
[ 0 0 0 -1 1]]
|
||||
|
||||
The matrix 'M_from_diff0' is bi-diagonal and could be alternatively
|
||||
created directly by
|
||||
|
||||
>>> M = - np.eye(n - 1, n, dtype=int)
|
||||
>>> np.fill_diagonal(M[:,1:], 1)
|
||||
>>> np.allclose(M, M_from_diff0)
|
||||
True
|
||||
|
||||
Its transpose
|
||||
|
||||
>>> print(M.T)
|
||||
[[-1 0 0 0]
|
||||
[ 1 -1 0 0]
|
||||
[ 0 1 -1 0]
|
||||
[ 0 0 1 -1]
|
||||
[ 0 0 0 1]]
|
||||
|
||||
can be viewed as the incidence matrix; see
|
||||
Incidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU,
|
||||
of a linear graph with 5 vertices and 4 edges. The 5x5 normal matrix
|
||||
'M.T @ M' thus is
|
||||
|
||||
>>> print(M.T @ M)
|
||||
[[ 1 -1 0 0 0]
|
||||
[-1 2 -1 0 0]
|
||||
[ 0 -1 2 -1 0]
|
||||
[ 0 0 -1 2 -1]
|
||||
[ 0 0 0 -1 1]]
|
||||
|
||||
the graph Laplacian, while the actually used in 'svds' smaller size
|
||||
4x4 normal matrix 'M @ M.T'
|
||||
|
||||
>>> print(M @ M.T)
|
||||
[[ 2 -1 0 0]
|
||||
[-1 2 -1 0]
|
||||
[ 0 -1 2 -1]
|
||||
[ 0 0 -1 2]]
|
||||
|
||||
is the so-called edge-based Laplacian; see
|
||||
Symmetric Laplacian via the incidence matrix, in Laplacian matrix,
|
||||
(2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.
|
||||
|
||||
The 'LinearOperator' setup needs the options 'rmatvec' and 'rmatmat'
|
||||
of multiplication by the matrix transpose 'M.T', but we want to be
|
||||
matrix-free to save memory, so knowing how 'M.T' looks like, we
|
||||
manually construct the following function to be used in 'rmatmat=diff0t'.
|
||||
|
||||
>>> def diff0t(a):
|
||||
... if a.ndim == 1:
|
||||
... a = a[:,np.newaxis] # Turn 1D into 2D array
|
||||
... d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
|
||||
... d[0, :] = - a[0, :]
|
||||
... d[1:-1, :] = a[0:-1, :] - a[1:, :]
|
||||
... d[-1, :] = a[-1, :]
|
||||
... return d
|
||||
|
||||
We check that our function 'diff0t' for the matrix transpose is valid.
|
||||
|
||||
>>> np.allclose(M.T, diff0t(np.eye(n-1)))
|
||||
True
|
||||
|
||||
Now we setup our matrix-free 'LinearOperator' called 'diff0_func_aslo'
|
||||
and for validation the matrix-based 'diff0_matrix_aslo'.
|
||||
|
||||
>>> def diff0_func_aslo_def(n):
|
||||
... return LinearOperator(matvec=diff0,
|
||||
... matmat=diff0,
|
||||
... rmatvec=diff0t,
|
||||
... rmatmat=diff0t,
|
||||
... shape=(n - 1, n))
|
||||
>>> diff0_func_aslo = diff0_func_aslo_def(n)
|
||||
>>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)
|
||||
|
||||
And validate both the matrix and its transpose in 'LinearOperator'.
|
||||
|
||||
>>> np.allclose(diff0_func_aslo(np.eye(n)),
|
||||
... diff0_matrix_aslo(np.eye(n)))
|
||||
True
|
||||
>>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
|
||||
... diff0_matrix_aslo.T(np.eye(n-1)))
|
||||
True
|
||||
|
||||
Having the 'LinearOperator' setup validated, we run the solver.
|
||||
|
||||
>>> n = 100
|
||||
>>> diff0_func_aslo = diff0_func_aslo_def(n)
|
||||
>>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')
|
||||
|
||||
The singular values squared and the singular vectors are known
|
||||
explicitly; see
|
||||
Pure Dirichlet boundary conditions, in
|
||||
Eigenvalues and eigenvectors of the second derivative,
|
||||
(2022, Nov. 19), Wikipedia, https://w.wiki/5YX6,
|
||||
since 'diff' corresponds to first
|
||||
derivative, and its smaller size n-1 x n-1 normal matrix
|
||||
'M @ M.T' represent the discrete second derivative with the Dirichlet
|
||||
boundary conditions. We use these analytic expressions for validation.
|
||||
|
||||
>>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
|
||||
>>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
|
||||
... np.arange(1, 4)) / n)
|
||||
>>> np.allclose(s, se, atol=1e-3)
|
||||
True
|
||||
>>> print(np.allclose(np.abs(u), np.abs(ue), atol=1e-6))
|
||||
True
|
||||
"""
|
||||
args = _iv(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors,
|
||||
solver, random_state)
|
||||
(A, k, ncv, tol, which, v0, maxiter,
|
||||
return_singular_vectors, solver, random_state) = args
|
||||
|
||||
largest = (which == 'LM')
|
||||
n, m = A.shape
|
||||
|
||||
if n >= m:
|
||||
X_dot = A.matvec
|
||||
X_matmat = A.matmat
|
||||
XH_dot = A.rmatvec
|
||||
XH_mat = A.rmatmat
|
||||
transpose = False
|
||||
else:
|
||||
X_dot = A.rmatvec
|
||||
X_matmat = A.rmatmat
|
||||
XH_dot = A.matvec
|
||||
XH_mat = A.matmat
|
||||
transpose = True
|
||||
|
||||
dtype = getattr(A, 'dtype', None)
|
||||
if dtype is None:
|
||||
dtype = A.dot(np.zeros([m, 1])).dtype
|
||||
|
||||
def matvec_XH_X(x):
|
||||
return XH_dot(X_dot(x))
|
||||
|
||||
def matmat_XH_X(x):
|
||||
return XH_mat(X_matmat(x))
|
||||
|
||||
XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype,
|
||||
matmat=matmat_XH_X,
|
||||
shape=(min(A.shape), min(A.shape)))
|
||||
|
||||
# Get a low rank approximation of the implicitly defined gramian matrix.
|
||||
# This is not a stable way to approach the problem.
|
||||
if solver == 'lobpcg':
|
||||
|
||||
if k == 1 and v0 is not None:
|
||||
X = np.reshape(v0, (-1, 1))
|
||||
else:
|
||||
X = random_state.standard_normal(size=(min(A.shape), k))
|
||||
|
||||
_, eigvec = lobpcg(XH_X, X, tol=tol ** 2, maxiter=maxiter,
|
||||
largest=largest)
|
||||
# lobpcg does not guarantee exactly orthonormal eigenvectors
|
||||
# until after gh-16320 is merged
|
||||
eigvec, _ = np.linalg.qr(eigvec)
|
||||
|
||||
elif solver == 'propack':
|
||||
if not HAS_PROPACK:
|
||||
raise ValueError("`solver='propack'` is opt-in due "
|
||||
"to potential issues on Windows, "
|
||||
"it can be enabled by setting the "
|
||||
"`SCIPY_USE_PROPACK` environment "
|
||||
"variable before importing scipy")
|
||||
jobu = return_singular_vectors in {True, 'u'}
|
||||
jobv = return_singular_vectors in {True, 'vh'}
|
||||
irl_mode = (which == 'SM')
|
||||
res = _svdp(A, k=k, tol=tol**2, which=which, maxiter=None,
|
||||
compute_u=jobu, compute_v=jobv, irl_mode=irl_mode,
|
||||
kmax=maxiter, v0=v0, random_state=random_state)
|
||||
|
||||
u, s, vh, _ = res # but we'll ignore bnd, the last output
|
||||
|
||||
# PROPACK order appears to be largest first. `svds` output order is not
|
||||
# guaranteed, according to documentation, but for ARPACK and LOBPCG
|
||||
# they actually are ordered smallest to largest, so reverse for
|
||||
# consistency.
|
||||
s = s[::-1]
|
||||
u = u[:, ::-1]
|
||||
vh = vh[::-1]
|
||||
|
||||
u = u if jobu else None
|
||||
vh = vh if jobv else None
|
||||
|
||||
if return_singular_vectors:
|
||||
return u, s, vh
|
||||
else:
|
||||
return s
|
||||
|
||||
elif solver == 'arpack' or solver is None:
|
||||
if v0 is None:
|
||||
v0 = random_state.standard_normal(size=(min(A.shape),))
|
||||
_, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter,
|
||||
ncv=ncv, which=which, v0=v0)
|
||||
# arpack do not guarantee exactly orthonormal eigenvectors
|
||||
# for clustered eigenvalues, especially in complex arithmetic
|
||||
eigvec, _ = np.linalg.qr(eigvec)
|
||||
|
||||
# the eigenvectors eigvec must be orthonomal here; see gh-16712
|
||||
Av = X_matmat(eigvec)
|
||||
if not return_singular_vectors:
|
||||
s = svd(Av, compute_uv=False, overwrite_a=True)
|
||||
return s[::-1]
|
||||
|
||||
# compute the left singular vectors of X and update the right ones
|
||||
# accordingly
|
||||
u, s, vh = svd(Av, full_matrices=False, overwrite_a=True)
|
||||
u = u[:, ::-1]
|
||||
s = s[::-1]
|
||||
vh = vh[::-1]
|
||||
|
||||
jobu = return_singular_vectors in {True, 'u'}
|
||||
jobv = return_singular_vectors in {True, 'vh'}
|
||||
|
||||
if transpose:
|
||||
u_tmp = eigvec @ _herm(vh) if jobu else None
|
||||
vh = _herm(u) if jobv else None
|
||||
u = u_tmp
|
||||
else:
|
||||
if not jobu:
|
||||
u = None
|
||||
vh = vh @ _herm(eigvec) if jobv else None
|
||||
|
||||
return u, s, vh
|
||||
@@ -1,398 +0,0 @@
|
||||
|
||||
def _svds_arpack_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
|
||||
maxiter=None, return_singular_vectors=True,
|
||||
solver='arpack', random_state=None):
|
||||
"""
|
||||
Partial singular value decomposition of a sparse matrix using ARPACK.
|
||||
|
||||
Compute the largest or smallest `k` singular values and corresponding
|
||||
singular vectors of a sparse matrix `A`. The order in which the singular
|
||||
values are returned is not guaranteed.
|
||||
|
||||
In the descriptions below, let ``M, N = A.shape``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : sparse matrix or LinearOperator
|
||||
Matrix to decompose.
|
||||
k : int, optional
|
||||
Number of singular values and singular vectors to compute.
|
||||
Must satisfy ``1 <= k <= min(M, N) - 1``.
|
||||
Default is 6.
|
||||
ncv : int, optional
|
||||
The number of Lanczos vectors generated.
|
||||
The default is ``min(n, max(2*k + 1, 20))``.
|
||||
If specified, must satistify ``k + 1 < ncv < min(M, N)``; ``ncv > 2*k``
|
||||
is recommended.
|
||||
tol : float, optional
|
||||
Tolerance for singular values. Zero (default) means machine precision.
|
||||
which : {'LM', 'SM'}
|
||||
Which `k` singular values to find: either the largest magnitude ('LM')
|
||||
or smallest magnitude ('SM') singular values.
|
||||
v0 : ndarray, optional
|
||||
The starting vector for iteration:
|
||||
an (approximate) left singular vector if ``N > M`` and a right singular
|
||||
vector otherwise. Must be of length ``min(M, N)``.
|
||||
Default: random
|
||||
maxiter : int, optional
|
||||
Maximum number of Arnoldi update iterations allowed;
|
||||
default is ``min(M, N) * 10``.
|
||||
return_singular_vectors : {True, False, "u", "vh"}
|
||||
Singular values are always computed and returned; this parameter
|
||||
controls the computation and return of singular vectors.
|
||||
|
||||
- ``True``: return singular vectors.
|
||||
- ``False``: do not return singular vectors.
|
||||
- ``"u"``: if ``M <= N``, compute only the left singular vectors and
|
||||
return ``None`` for the right singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
- ``"vh"``: if ``M > N``, compute only the right singular vectors and
|
||||
return ``None`` for the left singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
|
||||
solver : {'arpack', 'propack', 'lobpcg'}, optional
|
||||
This is the solver-specific documentation for ``solver='arpack'``.
|
||||
:ref:`'lobpcg' <sparse.linalg.svds-lobpcg>` and
|
||||
:ref:`'propack' <sparse.linalg.svds-propack>`
|
||||
are also supported.
|
||||
random_state : {None, int, `numpy.random.Generator`,
|
||||
`numpy.random.RandomState`}, optional
|
||||
|
||||
Pseudorandom number generator state used to generate resamples.
|
||||
|
||||
If `random_state` is ``None`` (or `np.random`), the
|
||||
`numpy.random.RandomState` singleton is used.
|
||||
If `random_state` is an int, a new ``RandomState`` instance is used,
|
||||
seeded with `random_state`.
|
||||
If `random_state` is already a ``Generator`` or ``RandomState``
|
||||
instance then that instance is used.
|
||||
options : dict, optional
|
||||
A dictionary of solver-specific options. No solver-specific options
|
||||
are currently supported; this parameter is reserved for future use.
|
||||
|
||||
Returns
|
||||
-------
|
||||
u : ndarray, shape=(M, k)
|
||||
Unitary matrix having left singular vectors as columns.
|
||||
s : ndarray, shape=(k,)
|
||||
The singular values.
|
||||
vh : ndarray, shape=(k, N)
|
||||
Unitary matrix having right singular vectors as rows.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This is a naive implementation using ARPACK as an eigensolver
|
||||
on ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on which one is more
|
||||
efficient.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Construct a matrix ``A`` from singular values and vectors.
|
||||
|
||||
>>> from scipy.stats import ortho_group
|
||||
>>> from scipy.sparse import csc_matrix, diags
|
||||
>>> from scipy.sparse.linalg import svds
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> orthogonal = csc_matrix(ortho_group.rvs(10, random_state=rng))
|
||||
>>> s = [0.0001, 0.001, 3, 4, 5] # singular values
|
||||
>>> u = orthogonal[:, :5] # left singular vectors
|
||||
>>> vT = orthogonal[:, 5:].T # right singular vectors
|
||||
>>> A = u @ diags(s) @ vT
|
||||
|
||||
With only three singular values/vectors, the SVD approximates the original
|
||||
matrix.
|
||||
|
||||
>>> u2, s2, vT2 = svds(A, k=3, solver='arpack')
|
||||
>>> A2 = u2 @ np.diag(s2) @ vT2
|
||||
>>> np.allclose(A2, A.toarray(), atol=1e-3)
|
||||
True
|
||||
|
||||
With all five singular values/vectors, we can reproduce the original
|
||||
matrix.
|
||||
|
||||
>>> u3, s3, vT3 = svds(A, k=5, solver='arpack')
|
||||
>>> A3 = u3 @ np.diag(s3) @ vT3
|
||||
>>> np.allclose(A3, A.toarray())
|
||||
True
|
||||
|
||||
The singular values match the expected singular values, and the singular
|
||||
vectors are as expected up to a difference in sign.
|
||||
|
||||
>>> (np.allclose(s3, s) and
|
||||
... np.allclose(np.abs(u3), np.abs(u.toarray())) and
|
||||
... np.allclose(np.abs(vT3), np.abs(vT.toarray())))
|
||||
True
|
||||
|
||||
The singular vectors are also orthogonal.
|
||||
|
||||
>>> (np.allclose(u3.T @ u3, np.eye(5)) and
|
||||
... np.allclose(vT3 @ vT3.T, np.eye(5)))
|
||||
True
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
def _svds_lobpcg_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
|
||||
maxiter=None, return_singular_vectors=True,
|
||||
solver='lobpcg', random_state=None):
|
||||
"""
|
||||
Partial singular value decomposition of a sparse matrix using LOBPCG.
|
||||
|
||||
Compute the largest or smallest `k` singular values and corresponding
|
||||
singular vectors of a sparse matrix `A`. The order in which the singular
|
||||
values are returned is not guaranteed.
|
||||
|
||||
In the descriptions below, let ``M, N = A.shape``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : sparse matrix or LinearOperator
|
||||
Matrix to decompose.
|
||||
k : int, default: 6
|
||||
Number of singular values and singular vectors to compute.
|
||||
Must satisfy ``1 <= k <= min(M, N) - 1``.
|
||||
ncv : int, optional
|
||||
Ignored.
|
||||
tol : float, optional
|
||||
Tolerance for singular values. Zero (default) means machine precision.
|
||||
which : {'LM', 'SM'}
|
||||
Which `k` singular values to find: either the largest magnitude ('LM')
|
||||
or smallest magnitude ('SM') singular values.
|
||||
v0 : ndarray, optional
|
||||
If `k` is 1, the starting vector for iteration:
|
||||
an (approximate) left singular vector if ``N > M`` and a right singular
|
||||
vector otherwise. Must be of length ``min(M, N)``.
|
||||
Ignored otherwise.
|
||||
Default: random
|
||||
maxiter : int, default: 20
|
||||
Maximum number of iterations.
|
||||
return_singular_vectors : {True, False, "u", "vh"}
|
||||
Singular values are always computed and returned; this parameter
|
||||
controls the computation and return of singular vectors.
|
||||
|
||||
- ``True``: return singular vectors.
|
||||
- ``False``: do not return singular vectors.
|
||||
- ``"u"``: if ``M <= N``, compute only the left singular vectors and
|
||||
return ``None`` for the right singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
- ``"vh"``: if ``M > N``, compute only the right singular vectors and
|
||||
return ``None`` for the left singular vectors. Otherwise, compute
|
||||
all singular vectors.
|
||||
|
||||
solver : {'arpack', 'propack', 'lobpcg'}, optional
|
||||
This is the solver-specific documentation for ``solver='lobpcg'``.
|
||||
:ref:`'arpack' <sparse.linalg.svds-arpack>` and
|
||||
:ref:`'propack' <sparse.linalg.svds-propack>`
|
||||
are also supported.
|
||||
random_state : {None, int, `numpy.random.Generator`,
|
||||
`numpy.random.RandomState`}, optional
|
||||
|
||||
Pseudorandom number generator state used to generate resamples.
|
||||
|
||||
If `random_state` is ``None`` (or `np.random`), the
|
||||
`numpy.random.RandomState` singleton is used.
|
||||
If `random_state` is an int, a new ``RandomState`` instance is used,
|
||||
seeded with `random_state`.
|
||||
If `random_state` is already a ``Generator`` or ``RandomState``
|
||||
instance then that instance is used.
|
||||
options : dict, optional
|
||||
A dictionary of solver-specific options. No solver-specific options
|
||||
are currently supported; this parameter is reserved for future use.
|
||||
|
||||
Returns
|
||||
-------
|
||||
u : ndarray, shape=(M, k)
|
||||
Unitary matrix having left singular vectors as columns.
|
||||
s : ndarray, shape=(k,)
|
||||
The singular values.
|
||||
vh : ndarray, shape=(k, N)
|
||||
Unitary matrix having right singular vectors as rows.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This is a naive implementation using LOBPCG as an eigensolver
|
||||
on ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on which one is more
|
||||
efficient.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Construct a matrix ``A`` from singular values and vectors.
|
||||
|
||||
>>> from scipy.stats import ortho_group
|
||||
>>> from scipy.sparse import csc_matrix, diags
|
||||
>>> from scipy.sparse.linalg import svds
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> orthogonal = csc_matrix(ortho_group.rvs(10, random_state=rng))
|
||||
>>> s = [0.0001, 0.001, 3, 4, 5] # singular values
|
||||
>>> u = orthogonal[:, :5] # left singular vectors
|
||||
>>> vT = orthogonal[:, 5:].T # right singular vectors
|
||||
>>> A = u @ diags(s) @ vT
|
||||
|
||||
With only three singular values/vectors, the SVD approximates the original
|
||||
matrix.
|
||||
|
||||
>>> u2, s2, vT2 = svds(A, k=3, solver='lobpcg')
|
||||
>>> A2 = u2 @ np.diag(s2) @ vT2
|
||||
>>> np.allclose(A2, A.toarray(), atol=1e-3)
|
||||
True
|
||||
|
||||
With all five singular values/vectors, we can reproduce the original
|
||||
matrix.
|
||||
|
||||
>>> u3, s3, vT3 = svds(A, k=5, solver='lobpcg')
|
||||
>>> A3 = u3 @ np.diag(s3) @ vT3
|
||||
>>> np.allclose(A3, A.toarray())
|
||||
True
|
||||
|
||||
The singular values match the expected singular values, and the singular
|
||||
vectors are as expected up to a difference in sign.
|
||||
|
||||
>>> (np.allclose(s3, s) and
|
||||
... np.allclose(np.abs(u3), np.abs(u.todense())) and
|
||||
... np.allclose(np.abs(vT3), np.abs(vT.todense())))
|
||||
True
|
||||
|
||||
The singular vectors are also orthogonal.
|
||||
|
||||
>>> (np.allclose(u3.T @ u3, np.eye(5)) and
|
||||
... np.allclose(vT3 @ vT3.T, np.eye(5)))
|
||||
True
|
||||
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
def _svds_propack_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
|
||||
maxiter=None, return_singular_vectors=True,
|
||||
solver='propack', random_state=None):
|
||||
"""
|
||||
Partial singular value decomposition of a sparse matrix using PROPACK.
|
||||
|
||||
Compute the largest or smallest `k` singular values and corresponding
|
||||
singular vectors of a sparse matrix `A`. The order in which the singular
|
||||
values are returned is not guaranteed.
|
||||
|
||||
In the descriptions below, let ``M, N = A.shape``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : sparse matrix or LinearOperator
|
||||
Matrix to decompose. If `A` is a ``LinearOperator``
|
||||
object, it must define both ``matvec`` and ``rmatvec`` methods.
|
||||
k : int, default: 6
|
||||
Number of singular values and singular vectors to compute.
|
||||
Must satisfy ``1 <= k <= min(M, N)``.
|
||||
ncv : int, optional
|
||||
Ignored.
|
||||
tol : float, optional
|
||||
The desired relative accuracy for computed singular values.
|
||||
Zero (default) means machine precision.
|
||||
which : {'LM', 'SM'}
|
||||
Which `k` singular values to find: either the largest magnitude ('LM')
|
||||
or smallest magnitude ('SM') singular values. Note that choosing
|
||||
``which='SM'`` will force the ``irl`` option to be set ``True``.
|
||||
v0 : ndarray, optional
|
||||
Starting vector for iterations: must be of length ``A.shape[0]``.
|
||||
If not specified, PROPACK will generate a starting vector.
|
||||
maxiter : int, optional
|
||||
Maximum number of iterations / maximal dimension of the Krylov
|
||||
subspace. Default is ``10 * k``.
|
||||
return_singular_vectors : {True, False, "u", "vh"}
|
||||
Singular values are always computed and returned; this parameter
|
||||
controls the computation and return of singular vectors.
|
||||
|
||||
- ``True``: return singular vectors.
|
||||
- ``False``: do not return singular vectors.
|
||||
- ``"u"``: compute only the left singular vectors; return ``None`` for
|
||||
the right singular vectors.
|
||||
- ``"vh"``: compute only the right singular vectors; return ``None``
|
||||
for the left singular vectors.
|
||||
|
||||
solver : {'arpack', 'propack', 'lobpcg'}, optional
|
||||
This is the solver-specific documentation for ``solver='propack'``.
|
||||
:ref:`'arpack' <sparse.linalg.svds-arpack>` and
|
||||
:ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`
|
||||
are also supported.
|
||||
random_state : {None, int, `numpy.random.Generator`,
|
||||
`numpy.random.RandomState`}, optional
|
||||
|
||||
Pseudorandom number generator state used to generate resamples.
|
||||
|
||||
If `random_state` is ``None`` (or `np.random`), the
|
||||
`numpy.random.RandomState` singleton is used.
|
||||
If `random_state` is an int, a new ``RandomState`` instance is used,
|
||||
seeded with `random_state`.
|
||||
If `random_state` is already a ``Generator`` or ``RandomState``
|
||||
instance then that instance is used.
|
||||
options : dict, optional
|
||||
A dictionary of solver-specific options. No solver-specific options
|
||||
are currently supported; this parameter is reserved for future use.
|
||||
|
||||
Returns
|
||||
-------
|
||||
u : ndarray, shape=(M, k)
|
||||
Unitary matrix having left singular vectors as columns.
|
||||
s : ndarray, shape=(k,)
|
||||
The singular values.
|
||||
vh : ndarray, shape=(k, N)
|
||||
Unitary matrix having right singular vectors as rows.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This is an interface to the Fortran library PROPACK [1]_.
|
||||
The current default is to run with IRL mode disabled unless seeking the
|
||||
smallest singular values/vectors (``which='SM'``).
|
||||
|
||||
References
|
||||
----------
|
||||
|
||||
.. [1] Larsen, Rasmus Munk. "PROPACK-Software for large and sparse SVD
|
||||
calculations." Available online. URL
|
||||
http://sun.stanford.edu/~rmunk/PROPACK (2004): 2008-2009.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Construct a matrix ``A`` from singular values and vectors.
|
||||
|
||||
>>> from scipy.stats import ortho_group
|
||||
>>> from scipy.sparse import csc_matrix, diags
|
||||
>>> from scipy.sparse.linalg import svds
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> orthogonal = csc_matrix(ortho_group.rvs(10, random_state=rng))
|
||||
>>> s = [0.0001, 0.001, 3, 4, 5] # singular values
|
||||
>>> u = orthogonal[:, :5] # left singular vectors
|
||||
>>> vT = orthogonal[:, 5:].T # right singular vectors
|
||||
>>> A = u @ diags(s) @ vT
|
||||
|
||||
With only three singular values/vectors, the SVD approximates the original
|
||||
matrix.
|
||||
|
||||
>>> u2, s2, vT2 = svds(A, k=3, solver='propack')
|
||||
>>> A2 = u2 @ np.diag(s2) @ vT2
|
||||
>>> np.allclose(A2, A.todense(), atol=1e-3)
|
||||
True
|
||||
|
||||
With all five singular values/vectors, we can reproduce the original
|
||||
matrix.
|
||||
|
||||
>>> u3, s3, vT3 = svds(A, k=5, solver='propack')
|
||||
>>> A3 = u3 @ np.diag(s3) @ vT3
|
||||
>>> np.allclose(A3, A.todense())
|
||||
True
|
||||
|
||||
The singular values match the expected singular values, and the singular
|
||||
vectors are as expected up to a difference in sign.
|
||||
|
||||
>>> (np.allclose(s3, s) and
|
||||
... np.allclose(np.abs(u3), np.abs(u.toarray())) and
|
||||
... np.allclose(np.abs(vT3), np.abs(vT.toarray())))
|
||||
True
|
||||
|
||||
The singular vectors are also orthogonal.
|
||||
|
||||
>>> (np.allclose(u3.T @ u3, np.eye(5)) and
|
||||
... np.allclose(vT3 @ vT3.T, np.eye(5)))
|
||||
True
|
||||
|
||||
"""
|
||||
pass
|
||||
@@ -1,45 +0,0 @@
|
||||
|
||||
BSD Software License
|
||||
|
||||
Pertains to ARPACK and P_ARPACK
|
||||
|
||||
Copyright (c) 1996-2008 Rice University.
|
||||
Developed by D.C. Sorensen, R.B. Lehoucq, C. Yang, and K. Maschhoff.
|
||||
All rights reserved.
|
||||
|
||||
Arpack has been renamed to arpack-ng.
|
||||
|
||||
Copyright (c) 2001-2011 - Scilab Enterprises
|
||||
Updated by Allan Cornet, Sylvestre Ledru.
|
||||
|
||||
Copyright (c) 2010 - Jordi Gutiérrez Hermoso (Octave patch)
|
||||
|
||||
Copyright (c) 2007 - Sébastien Fabbro (gentoo patch)
|
||||
|
||||
Redistribution and use in source and binary forms, with or without
|
||||
modification, are permitted provided that the following conditions are
|
||||
met:
|
||||
|
||||
- Redistributions of source code must retain the above copyright
|
||||
notice, this list of conditions and the following disclaimer.
|
||||
|
||||
- Redistributions in binary form must reproduce the above copyright
|
||||
notice, this list of conditions and the following disclaimer listed
|
||||
in this license in the documentation and/or other materials
|
||||
provided with the distribution.
|
||||
|
||||
- Neither the name of the copyright holders nor the names of its
|
||||
contributors may be used to endorse or promote products derived from
|
||||
this software without specific prior written permission.
|
||||
|
||||
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
||||
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
||||
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
||||
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
||||
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
||||
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
||||
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
||||
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
||||
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
||||
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
||||
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||||
@@ -1,20 +0,0 @@
|
||||
"""
|
||||
Eigenvalue solver using iterative methods.
|
||||
|
||||
Find k eigenvectors and eigenvalues of a matrix A using the
|
||||
Arnoldi/Lanczos iterative methods from ARPACK [1]_,[2]_.
|
||||
|
||||
These methods are most useful for large sparse matrices.
|
||||
|
||||
- eigs(A,k)
|
||||
- eigsh(A,k)
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] ARPACK Software, http://www.caam.rice.edu/software/ARPACK/
|
||||
.. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE:
|
||||
Solution of Large Scale Eigenvalue Problems by Implicitly Restarted
|
||||
Arnoldi Methods. SIAM, Philadelphia, PA, 1998.
|
||||
|
||||
"""
|
||||
from .arpack import *
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
@@ -1,725 +0,0 @@
|
||||
__usage__ = """
|
||||
To run tests locally:
|
||||
python tests/test_arpack.py [-l<int>] [-v<int>]
|
||||
|
||||
"""
|
||||
|
||||
import threading
|
||||
import itertools
|
||||
|
||||
import numpy as np
|
||||
|
||||
from numpy.testing import assert_allclose, assert_equal, suppress_warnings
|
||||
from pytest import raises as assert_raises
|
||||
import pytest
|
||||
|
||||
from numpy import dot, conj, random
|
||||
from scipy.linalg import eig, eigh
|
||||
from scipy.sparse import csc_matrix, csr_matrix, diags, rand
|
||||
from scipy.sparse.linalg import LinearOperator, aslinearoperator
|
||||
from scipy.sparse.linalg._eigen.arpack import (eigs, eigsh, arpack,
|
||||
ArpackNoConvergence)
|
||||
|
||||
|
||||
from scipy._lib._gcutils import assert_deallocated, IS_PYPY
|
||||
|
||||
|
||||
# precision for tests
|
||||
_ndigits = {'f': 3, 'd': 11, 'F': 3, 'D': 11}
|
||||
|
||||
|
||||
def _get_test_tolerance(type_char, mattype=None, D_type=None, which=None):
|
||||
"""
|
||||
Return tolerance values suitable for a given test:
|
||||
|
||||
Parameters
|
||||
----------
|
||||
type_char : {'f', 'd', 'F', 'D'}
|
||||
Data type in ARPACK eigenvalue problem
|
||||
mattype : {csr_matrix, aslinearoperator, asarray}, optional
|
||||
Linear operator type
|
||||
|
||||
Returns
|
||||
-------
|
||||
tol
|
||||
Tolerance to pass to the ARPACK routine
|
||||
rtol
|
||||
Relative tolerance for outputs
|
||||
atol
|
||||
Absolute tolerance for outputs
|
||||
|
||||
"""
|
||||
|
||||
rtol = {'f': 3000 * np.finfo(np.float32).eps,
|
||||
'F': 3000 * np.finfo(np.float32).eps,
|
||||
'd': 2000 * np.finfo(np.float64).eps,
|
||||
'D': 2000 * np.finfo(np.float64).eps}[type_char]
|
||||
atol = rtol
|
||||
tol = 0
|
||||
|
||||
if mattype is aslinearoperator and type_char in ('f', 'F'):
|
||||
# iterative methods in single precision: worse errors
|
||||
# also: bump ARPACK tolerance so that the iterative method converges
|
||||
tol = 30 * np.finfo(np.float32).eps
|
||||
rtol *= 5
|
||||
|
||||
if mattype is csr_matrix and type_char in ('f', 'F'):
|
||||
# sparse in single precision: worse errors
|
||||
rtol *= 5
|
||||
|
||||
if (
|
||||
which in ('LM', 'SM', 'LA')
|
||||
and D_type.name == "gen-hermitian-Mc"
|
||||
):
|
||||
if type_char == 'F':
|
||||
# missing case 1, 2, and more, from PR 14798
|
||||
rtol *= 5
|
||||
|
||||
if type_char == 'D':
|
||||
# missing more cases, from PR 14798
|
||||
rtol *= 7
|
||||
|
||||
return tol, rtol, atol
|
||||
|
||||
|
||||
def generate_matrix(N, complex_=False, hermitian=False,
|
||||
pos_definite=False, sparse=False):
|
||||
M = np.random.random((N, N))
|
||||
if complex_:
|
||||
M = M + 1j * np.random.random((N, N))
|
||||
|
||||
if hermitian:
|
||||
if pos_definite:
|
||||
if sparse:
|
||||
i = np.arange(N)
|
||||
j = np.random.randint(N, size=N-2)
|
||||
i, j = np.meshgrid(i, j)
|
||||
M[i, j] = 0
|
||||
M = np.dot(M.conj(), M.T)
|
||||
else:
|
||||
M = np.dot(M.conj(), M.T)
|
||||
if sparse:
|
||||
i = np.random.randint(N, size=N * N // 4)
|
||||
j = np.random.randint(N, size=N * N // 4)
|
||||
ind = np.nonzero(i == j)
|
||||
j[ind] = (j[ind] + 1) % N
|
||||
M[i, j] = 0
|
||||
M[j, i] = 0
|
||||
else:
|
||||
if sparse:
|
||||
i = np.random.randint(N, size=N * N // 2)
|
||||
j = np.random.randint(N, size=N * N // 2)
|
||||
M[i, j] = 0
|
||||
return M
|
||||
|
||||
|
||||
def generate_matrix_symmetric(N, pos_definite=False, sparse=False):
|
||||
M = np.random.random((N, N))
|
||||
|
||||
M = 0.5 * (M + M.T) # Make M symmetric
|
||||
|
||||
if pos_definite:
|
||||
Id = N * np.eye(N)
|
||||
if sparse:
|
||||
M = csr_matrix(M)
|
||||
M += Id
|
||||
else:
|
||||
if sparse:
|
||||
M = csr_matrix(M)
|
||||
|
||||
return M
|
||||
|
||||
|
||||
def _aslinearoperator_with_dtype(m):
|
||||
m = aslinearoperator(m)
|
||||
if not hasattr(m, 'dtype'):
|
||||
x = np.zeros(m.shape[1])
|
||||
m.dtype = (m * x).dtype
|
||||
return m
|
||||
|
||||
|
||||
def assert_allclose_cc(actual, desired, **kw):
|
||||
"""Almost equal or complex conjugates almost equal"""
|
||||
try:
|
||||
assert_allclose(actual, desired, **kw)
|
||||
except AssertionError:
|
||||
assert_allclose(actual, conj(desired), **kw)
|
||||
|
||||
|
||||
def argsort_which(eigenvalues, typ, k, which,
|
||||
sigma=None, OPpart=None, mode=None):
|
||||
"""Return sorted indices of eigenvalues using the "which" keyword
|
||||
from eigs and eigsh"""
|
||||
if sigma is None:
|
||||
reval = np.round(eigenvalues, decimals=_ndigits[typ])
|
||||
else:
|
||||
if mode is None or mode == 'normal':
|
||||
if OPpart is None:
|
||||
reval = 1. / (eigenvalues - sigma)
|
||||
elif OPpart == 'r':
|
||||
reval = 0.5 * (1. / (eigenvalues - sigma)
|
||||
+ 1. / (eigenvalues - np.conj(sigma)))
|
||||
elif OPpart == 'i':
|
||||
reval = -0.5j * (1. / (eigenvalues - sigma)
|
||||
- 1. / (eigenvalues - np.conj(sigma)))
|
||||
elif mode == 'cayley':
|
||||
reval = (eigenvalues + sigma) / (eigenvalues - sigma)
|
||||
elif mode == 'buckling':
|
||||
reval = eigenvalues / (eigenvalues - sigma)
|
||||
else:
|
||||
raise ValueError("mode='%s' not recognized" % mode)
|
||||
|
||||
reval = np.round(reval, decimals=_ndigits[typ])
|
||||
|
||||
if which in ['LM', 'SM']:
|
||||
ind = np.argsort(abs(reval))
|
||||
elif which in ['LR', 'SR', 'LA', 'SA', 'BE']:
|
||||
ind = np.argsort(np.real(reval))
|
||||
elif which in ['LI', 'SI']:
|
||||
# for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
|
||||
if typ.islower():
|
||||
ind = np.argsort(abs(np.imag(reval)))
|
||||
else:
|
||||
ind = np.argsort(np.imag(reval))
|
||||
else:
|
||||
raise ValueError("which='%s' is unrecognized" % which)
|
||||
|
||||
if which in ['LM', 'LA', 'LR', 'LI']:
|
||||
return ind[-k:]
|
||||
elif which in ['SM', 'SA', 'SR', 'SI']:
|
||||
return ind[:k]
|
||||
elif which == 'BE':
|
||||
return np.concatenate((ind[:k//2], ind[k//2-k:]))
|
||||
|
||||
|
||||
def eval_evec(symmetric, d, typ, k, which, v0=None, sigma=None,
|
||||
mattype=np.asarray, OPpart=None, mode='normal'):
|
||||
general = ('bmat' in d)
|
||||
|
||||
if symmetric:
|
||||
eigs_func = eigsh
|
||||
else:
|
||||
eigs_func = eigs
|
||||
|
||||
if general:
|
||||
err = ("error for %s:general, typ=%s, which=%s, sigma=%s, "
|
||||
"mattype=%s, OPpart=%s, mode=%s" % (eigs_func.__name__,
|
||||
typ, which, sigma,
|
||||
mattype.__name__,
|
||||
OPpart, mode))
|
||||
else:
|
||||
err = ("error for %s:standard, typ=%s, which=%s, sigma=%s, "
|
||||
"mattype=%s, OPpart=%s, mode=%s" % (eigs_func.__name__,
|
||||
typ, which, sigma,
|
||||
mattype.__name__,
|
||||
OPpart, mode))
|
||||
|
||||
a = d['mat'].astype(typ)
|
||||
ac = mattype(a)
|
||||
|
||||
if general:
|
||||
b = d['bmat'].astype(typ)
|
||||
bc = mattype(b)
|
||||
|
||||
# get exact eigenvalues
|
||||
exact_eval = d['eval'].astype(typ.upper())
|
||||
ind = argsort_which(exact_eval, typ, k, which,
|
||||
sigma, OPpart, mode)
|
||||
exact_eval = exact_eval[ind]
|
||||
|
||||
# compute arpack eigenvalues
|
||||
kwargs = dict(which=which, v0=v0, sigma=sigma)
|
||||
if eigs_func is eigsh:
|
||||
kwargs['mode'] = mode
|
||||
else:
|
||||
kwargs['OPpart'] = OPpart
|
||||
|
||||
# compute suitable tolerances
|
||||
kwargs['tol'], rtol, atol = _get_test_tolerance(typ, mattype, d, which)
|
||||
# on rare occasions, ARPACK routines return results that are proper
|
||||
# eigenvalues and -vectors, but not necessarily the ones requested in
|
||||
# the parameter which. This is inherent to the Krylov methods, and
|
||||
# should not be treated as a failure. If such a rare situation
|
||||
# occurs, the calculation is tried again (but at most a few times).
|
||||
ntries = 0
|
||||
while ntries < 5:
|
||||
# solve
|
||||
if general:
|
||||
try:
|
||||
eigenvalues, evec = eigs_func(ac, k, bc, **kwargs)
|
||||
except ArpackNoConvergence:
|
||||
kwargs['maxiter'] = 20*a.shape[0]
|
||||
eigenvalues, evec = eigs_func(ac, k, bc, **kwargs)
|
||||
else:
|
||||
try:
|
||||
eigenvalues, evec = eigs_func(ac, k, **kwargs)
|
||||
except ArpackNoConvergence:
|
||||
kwargs['maxiter'] = 20*a.shape[0]
|
||||
eigenvalues, evec = eigs_func(ac, k, **kwargs)
|
||||
|
||||
ind = argsort_which(eigenvalues, typ, k, which,
|
||||
sigma, OPpart, mode)
|
||||
eigenvalues = eigenvalues[ind]
|
||||
evec = evec[:, ind]
|
||||
|
||||
try:
|
||||
# check eigenvalues
|
||||
assert_allclose_cc(eigenvalues, exact_eval, rtol=rtol, atol=atol,
|
||||
err_msg=err)
|
||||
check_evecs = True
|
||||
except AssertionError:
|
||||
check_evecs = False
|
||||
ntries += 1
|
||||
|
||||
if check_evecs:
|
||||
# check eigenvectors
|
||||
LHS = np.dot(a, evec)
|
||||
if general:
|
||||
RHS = eigenvalues * np.dot(b, evec)
|
||||
else:
|
||||
RHS = eigenvalues * evec
|
||||
|
||||
assert_allclose(LHS, RHS, rtol=rtol, atol=atol, err_msg=err)
|
||||
break
|
||||
|
||||
# check eigenvalues
|
||||
assert_allclose_cc(eigenvalues, exact_eval, rtol=rtol, atol=atol, err_msg=err)
|
||||
|
||||
|
||||
class DictWithRepr(dict):
|
||||
def __init__(self, name):
|
||||
self.name = name
|
||||
|
||||
def __repr__(self):
|
||||
return "<%s>" % self.name
|
||||
|
||||
|
||||
class SymmetricParams:
|
||||
def __init__(self):
|
||||
self.eigs = eigsh
|
||||
self.which = ['LM', 'SM', 'LA', 'SA', 'BE']
|
||||
self.mattypes = [csr_matrix, aslinearoperator, np.asarray]
|
||||
self.sigmas_modes = {None: ['normal'],
|
||||
0.5: ['normal', 'buckling', 'cayley']}
|
||||
|
||||
# generate matrices
|
||||
# these should all be float32 so that the eigenvalues
|
||||
# are the same in float32 and float64
|
||||
N = 6
|
||||
np.random.seed(2300)
|
||||
Ar = generate_matrix(N, hermitian=True,
|
||||
pos_definite=True).astype('f').astype('d')
|
||||
M = generate_matrix(N, hermitian=True,
|
||||
pos_definite=True).astype('f').astype('d')
|
||||
Ac = generate_matrix(N, hermitian=True, pos_definite=True,
|
||||
complex_=True).astype('F').astype('D')
|
||||
Mc = generate_matrix(N, hermitian=True, pos_definite=True,
|
||||
complex_=True).astype('F').astype('D')
|
||||
v0 = np.random.random(N)
|
||||
|
||||
# standard symmetric problem
|
||||
SS = DictWithRepr("std-symmetric")
|
||||
SS['mat'] = Ar
|
||||
SS['v0'] = v0
|
||||
SS['eval'] = eigh(SS['mat'], eigvals_only=True)
|
||||
|
||||
# general symmetric problem
|
||||
GS = DictWithRepr("gen-symmetric")
|
||||
GS['mat'] = Ar
|
||||
GS['bmat'] = M
|
||||
GS['v0'] = v0
|
||||
GS['eval'] = eigh(GS['mat'], GS['bmat'], eigvals_only=True)
|
||||
|
||||
# standard hermitian problem
|
||||
SH = DictWithRepr("std-hermitian")
|
||||
SH['mat'] = Ac
|
||||
SH['v0'] = v0
|
||||
SH['eval'] = eigh(SH['mat'], eigvals_only=True)
|
||||
|
||||
# general hermitian problem
|
||||
GH = DictWithRepr("gen-hermitian")
|
||||
GH['mat'] = Ac
|
||||
GH['bmat'] = M
|
||||
GH['v0'] = v0
|
||||
GH['eval'] = eigh(GH['mat'], GH['bmat'], eigvals_only=True)
|
||||
|
||||
# general hermitian problem with hermitian M
|
||||
GHc = DictWithRepr("gen-hermitian-Mc")
|
||||
GHc['mat'] = Ac
|
||||
GHc['bmat'] = Mc
|
||||
GHc['v0'] = v0
|
||||
GHc['eval'] = eigh(GHc['mat'], GHc['bmat'], eigvals_only=True)
|
||||
|
||||
self.real_test_cases = [SS, GS]
|
||||
self.complex_test_cases = [SH, GH, GHc]
|
||||
|
||||
|
||||
class NonSymmetricParams:
|
||||
def __init__(self):
|
||||
self.eigs = eigs
|
||||
self.which = ['LM', 'LR', 'LI'] # , 'SM', 'LR', 'SR', 'LI', 'SI']
|
||||
self.mattypes = [csr_matrix, aslinearoperator, np.asarray]
|
||||
self.sigmas_OPparts = {None: [None],
|
||||
0.1: ['r'],
|
||||
0.1 + 0.1j: ['r', 'i']}
|
||||
|
||||
# generate matrices
|
||||
# these should all be float32 so that the eigenvalues
|
||||
# are the same in float32 and float64
|
||||
N = 6
|
||||
np.random.seed(2300)
|
||||
Ar = generate_matrix(N).astype('f').astype('d')
|
||||
M = generate_matrix(N, hermitian=True,
|
||||
pos_definite=True).astype('f').astype('d')
|
||||
Ac = generate_matrix(N, complex_=True).astype('F').astype('D')
|
||||
v0 = np.random.random(N)
|
||||
|
||||
# standard real nonsymmetric problem
|
||||
SNR = DictWithRepr("std-real-nonsym")
|
||||
SNR['mat'] = Ar
|
||||
SNR['v0'] = v0
|
||||
SNR['eval'] = eig(SNR['mat'], left=False, right=False)
|
||||
|
||||
# general real nonsymmetric problem
|
||||
GNR = DictWithRepr("gen-real-nonsym")
|
||||
GNR['mat'] = Ar
|
||||
GNR['bmat'] = M
|
||||
GNR['v0'] = v0
|
||||
GNR['eval'] = eig(GNR['mat'], GNR['bmat'], left=False, right=False)
|
||||
|
||||
# standard complex nonsymmetric problem
|
||||
SNC = DictWithRepr("std-cmplx-nonsym")
|
||||
SNC['mat'] = Ac
|
||||
SNC['v0'] = v0
|
||||
SNC['eval'] = eig(SNC['mat'], left=False, right=False)
|
||||
|
||||
# general complex nonsymmetric problem
|
||||
GNC = DictWithRepr("gen-cmplx-nonsym")
|
||||
GNC['mat'] = Ac
|
||||
GNC['bmat'] = M
|
||||
GNC['v0'] = v0
|
||||
GNC['eval'] = eig(GNC['mat'], GNC['bmat'], left=False, right=False)
|
||||
|
||||
self.real_test_cases = [SNR, GNR]
|
||||
self.complex_test_cases = [SNC, GNC]
|
||||
|
||||
|
||||
def test_symmetric_modes():
|
||||
params = SymmetricParams()
|
||||
k = 2
|
||||
symmetric = True
|
||||
for D in params.real_test_cases:
|
||||
for typ in 'fd':
|
||||
for which in params.which:
|
||||
for mattype in params.mattypes:
|
||||
for (sigma, modes) in params.sigmas_modes.items():
|
||||
for mode in modes:
|
||||
eval_evec(symmetric, D, typ, k, which,
|
||||
None, sigma, mattype, None, mode)
|
||||
|
||||
|
||||
def test_hermitian_modes():
|
||||
params = SymmetricParams()
|
||||
k = 2
|
||||
symmetric = True
|
||||
for D in params.complex_test_cases:
|
||||
for typ in 'FD':
|
||||
for which in params.which:
|
||||
if which == 'BE':
|
||||
continue # BE invalid for complex
|
||||
for mattype in params.mattypes:
|
||||
for sigma in params.sigmas_modes:
|
||||
eval_evec(symmetric, D, typ, k, which,
|
||||
None, sigma, mattype)
|
||||
|
||||
|
||||
def test_symmetric_starting_vector():
|
||||
params = SymmetricParams()
|
||||
symmetric = True
|
||||
for k in [1, 2, 3, 4, 5]:
|
||||
for D in params.real_test_cases:
|
||||
for typ in 'fd':
|
||||
v0 = random.rand(len(D['v0'])).astype(typ)
|
||||
eval_evec(symmetric, D, typ, k, 'LM', v0)
|
||||
|
||||
|
||||
def test_symmetric_no_convergence():
|
||||
np.random.seed(1234)
|
||||
m = generate_matrix(30, hermitian=True, pos_definite=True)
|
||||
tol, rtol, atol = _get_test_tolerance('d')
|
||||
try:
|
||||
w, v = eigsh(m, 4, which='LM', v0=m[:, 0], maxiter=5, tol=tol, ncv=9)
|
||||
raise AssertionError("Spurious no-error exit")
|
||||
except ArpackNoConvergence as err:
|
||||
k = len(err.eigenvalues)
|
||||
if k <= 0:
|
||||
raise AssertionError("Spurious no-eigenvalues-found case") from err
|
||||
w, v = err.eigenvalues, err.eigenvectors
|
||||
assert_allclose(dot(m, v), w * v, rtol=rtol, atol=atol)
|
||||
|
||||
|
||||
def test_real_nonsymmetric_modes():
|
||||
params = NonSymmetricParams()
|
||||
k = 2
|
||||
symmetric = False
|
||||
for D in params.real_test_cases:
|
||||
for typ in 'fd':
|
||||
for which in params.which:
|
||||
for mattype in params.mattypes:
|
||||
for sigma, OPparts in params.sigmas_OPparts.items():
|
||||
for OPpart in OPparts:
|
||||
eval_evec(symmetric, D, typ, k, which,
|
||||
None, sigma, mattype, OPpart)
|
||||
|
||||
|
||||
def test_complex_nonsymmetric_modes():
|
||||
params = NonSymmetricParams()
|
||||
k = 2
|
||||
symmetric = False
|
||||
for D in params.complex_test_cases:
|
||||
for typ in 'DF':
|
||||
for which in params.which:
|
||||
for mattype in params.mattypes:
|
||||
for sigma in params.sigmas_OPparts:
|
||||
eval_evec(symmetric, D, typ, k, which,
|
||||
None, sigma, mattype)
|
||||
|
||||
|
||||
def test_standard_nonsymmetric_starting_vector():
|
||||
params = NonSymmetricParams()
|
||||
sigma = None
|
||||
symmetric = False
|
||||
for k in [1, 2, 3, 4]:
|
||||
for d in params.complex_test_cases:
|
||||
for typ in 'FD':
|
||||
A = d['mat']
|
||||
n = A.shape[0]
|
||||
v0 = random.rand(n).astype(typ)
|
||||
eval_evec(symmetric, d, typ, k, "LM", v0, sigma)
|
||||
|
||||
|
||||
def test_general_nonsymmetric_starting_vector():
|
||||
params = NonSymmetricParams()
|
||||
sigma = None
|
||||
symmetric = False
|
||||
for k in [1, 2, 3, 4]:
|
||||
for d in params.complex_test_cases:
|
||||
for typ in 'FD':
|
||||
A = d['mat']
|
||||
n = A.shape[0]
|
||||
v0 = random.rand(n).astype(typ)
|
||||
eval_evec(symmetric, d, typ, k, "LM", v0, sigma)
|
||||
|
||||
|
||||
def test_standard_nonsymmetric_no_convergence():
|
||||
np.random.seed(1234)
|
||||
m = generate_matrix(30, complex_=True)
|
||||
tol, rtol, atol = _get_test_tolerance('d')
|
||||
try:
|
||||
w, v = eigs(m, 4, which='LM', v0=m[:, 0], maxiter=5, tol=tol)
|
||||
raise AssertionError("Spurious no-error exit")
|
||||
except ArpackNoConvergence as err:
|
||||
k = len(err.eigenvalues)
|
||||
if k <= 0:
|
||||
raise AssertionError("Spurious no-eigenvalues-found case") from err
|
||||
w, v = err.eigenvalues, err.eigenvectors
|
||||
for ww, vv in zip(w, v.T):
|
||||
assert_allclose(dot(m, vv), ww * vv, rtol=rtol, atol=atol)
|
||||
|
||||
|
||||
def test_eigen_bad_shapes():
|
||||
# A is not square.
|
||||
A = csc_matrix(np.zeros((2, 3)))
|
||||
assert_raises(ValueError, eigs, A)
|
||||
|
||||
|
||||
def test_eigen_bad_kwargs():
|
||||
# Test eigen on wrong keyword argument
|
||||
A = csc_matrix(np.zeros((8, 8)))
|
||||
assert_raises(ValueError, eigs, A, which='XX')
|
||||
|
||||
|
||||
def test_ticket_1459_arpack_crash():
|
||||
for dtype in [np.float32, np.float64]:
|
||||
# This test does not seem to catch the issue for float32,
|
||||
# but we made the same fix there, just to be sure
|
||||
|
||||
N = 6
|
||||
k = 2
|
||||
|
||||
np.random.seed(2301)
|
||||
A = np.random.random((N, N)).astype(dtype)
|
||||
v0 = np.array([-0.71063568258907849895, -0.83185111795729227424,
|
||||
-0.34365925382227402451, 0.46122533684552280420,
|
||||
-0.58001341115969040629, -0.78844877570084292984e-01],
|
||||
dtype=dtype)
|
||||
|
||||
# Should not crash:
|
||||
evals, evecs = eigs(A, k, v0=v0)
|
||||
|
||||
|
||||
@pytest.mark.skipif(IS_PYPY, reason="Test not meaningful on PyPy")
|
||||
def test_linearoperator_deallocation():
|
||||
# Check that the linear operators used by the Arpack wrappers are
|
||||
# deallocatable by reference counting -- they are big objects, so
|
||||
# Python's cyclic GC may not collect them fast enough before
|
||||
# running out of memory if eigs/eigsh are called in a tight loop.
|
||||
|
||||
M_d = np.eye(10)
|
||||
M_s = csc_matrix(M_d)
|
||||
M_o = aslinearoperator(M_d)
|
||||
|
||||
with assert_deallocated(lambda: arpack.SpLuInv(M_s)):
|
||||
pass
|
||||
with assert_deallocated(lambda: arpack.LuInv(M_d)):
|
||||
pass
|
||||
with assert_deallocated(lambda: arpack.IterInv(M_s)):
|
||||
pass
|
||||
with assert_deallocated(lambda: arpack.IterOpInv(M_o, None, 0.3)):
|
||||
pass
|
||||
with assert_deallocated(lambda: arpack.IterOpInv(M_o, M_o, 0.3)):
|
||||
pass
|
||||
|
||||
def test_parallel_threads():
|
||||
results = []
|
||||
v0 = np.random.rand(50)
|
||||
|
||||
def worker():
|
||||
x = diags([1, -2, 1], [-1, 0, 1], shape=(50, 50))
|
||||
w, v = eigs(x, k=3, v0=v0)
|
||||
results.append(w)
|
||||
|
||||
w, v = eigsh(x, k=3, v0=v0)
|
||||
results.append(w)
|
||||
|
||||
threads = [threading.Thread(target=worker) for k in range(10)]
|
||||
for t in threads:
|
||||
t.start()
|
||||
for t in threads:
|
||||
t.join()
|
||||
|
||||
worker()
|
||||
|
||||
for r in results:
|
||||
assert_allclose(r, results[-1])
|
||||
|
||||
|
||||
def test_reentering():
|
||||
# Just some linear operator that calls eigs recursively
|
||||
def A_matvec(x):
|
||||
x = diags([1, -2, 1], [-1, 0, 1], shape=(50, 50))
|
||||
w, v = eigs(x, k=1)
|
||||
return v / w[0]
|
||||
A = LinearOperator(matvec=A_matvec, dtype=float, shape=(50, 50))
|
||||
|
||||
# The Fortran code is not reentrant, so this fails (gracefully, not crashing)
|
||||
assert_raises(RuntimeError, eigs, A, k=1)
|
||||
assert_raises(RuntimeError, eigsh, A, k=1)
|
||||
|
||||
|
||||
def test_regression_arpackng_1315():
|
||||
# Check that issue arpack-ng/#1315 is not present.
|
||||
# Adapted from arpack-ng/TESTS/bug_1315_single.c
|
||||
# If this fails, then the installed ARPACK library is faulty.
|
||||
|
||||
for dtype in [np.float32, np.float64]:
|
||||
np.random.seed(1234)
|
||||
|
||||
w0 = np.arange(1, 1000+1).astype(dtype)
|
||||
A = diags([w0], [0], shape=(1000, 1000))
|
||||
|
||||
v0 = np.random.rand(1000).astype(dtype)
|
||||
w, v = eigs(A, k=9, ncv=2*9+1, which="LM", v0=v0)
|
||||
|
||||
assert_allclose(np.sort(w), np.sort(w0[-9:]),
|
||||
rtol=1e-4)
|
||||
|
||||
|
||||
def test_eigs_for_k_greater():
|
||||
# Test eigs() for k beyond limits.
|
||||
A_sparse = diags([1, -2, 1], [-1, 0, 1], shape=(4, 4)) # sparse
|
||||
A = generate_matrix(4, sparse=False)
|
||||
M_dense = np.random.random((4, 4))
|
||||
M_sparse = generate_matrix(4, sparse=True)
|
||||
M_linop = aslinearoperator(M_dense)
|
||||
eig_tuple1 = eig(A, b=M_dense)
|
||||
eig_tuple2 = eig(A, b=M_sparse)
|
||||
|
||||
with suppress_warnings() as sup:
|
||||
sup.filter(RuntimeWarning)
|
||||
|
||||
assert_equal(eigs(A, M=M_dense, k=3), eig_tuple1)
|
||||
assert_equal(eigs(A, M=M_dense, k=4), eig_tuple1)
|
||||
assert_equal(eigs(A, M=M_dense, k=5), eig_tuple1)
|
||||
assert_equal(eigs(A, M=M_sparse, k=5), eig_tuple2)
|
||||
|
||||
# M as LinearOperator
|
||||
assert_raises(TypeError, eigs, A, M=M_linop, k=3)
|
||||
|
||||
# Test 'A' for different types
|
||||
assert_raises(TypeError, eigs, aslinearoperator(A), k=3)
|
||||
assert_raises(TypeError, eigs, A_sparse, k=3)
|
||||
|
||||
|
||||
def test_eigsh_for_k_greater():
|
||||
# Test eigsh() for k beyond limits.
|
||||
A_sparse = diags([1, -2, 1], [-1, 0, 1], shape=(4, 4)) # sparse
|
||||
A = generate_matrix(4, sparse=False)
|
||||
M_dense = generate_matrix_symmetric(4, pos_definite=True)
|
||||
M_sparse = generate_matrix_symmetric(4, pos_definite=True, sparse=True)
|
||||
M_linop = aslinearoperator(M_dense)
|
||||
eig_tuple1 = eigh(A, b=M_dense)
|
||||
eig_tuple2 = eigh(A, b=M_sparse)
|
||||
|
||||
with suppress_warnings() as sup:
|
||||
sup.filter(RuntimeWarning)
|
||||
|
||||
assert_equal(eigsh(A, M=M_dense, k=4), eig_tuple1)
|
||||
assert_equal(eigsh(A, M=M_dense, k=5), eig_tuple1)
|
||||
assert_equal(eigsh(A, M=M_sparse, k=5), eig_tuple2)
|
||||
|
||||
# M as LinearOperator
|
||||
assert_raises(TypeError, eigsh, A, M=M_linop, k=4)
|
||||
|
||||
# Test 'A' for different types
|
||||
assert_raises(TypeError, eigsh, aslinearoperator(A), k=4)
|
||||
assert_raises(TypeError, eigsh, A_sparse, M=M_dense, k=4)
|
||||
|
||||
|
||||
def test_real_eigs_real_k_subset():
|
||||
np.random.seed(1)
|
||||
|
||||
n = 10
|
||||
A = rand(n, n, density=0.5)
|
||||
A.data *= 2
|
||||
A.data -= 1
|
||||
|
||||
v0 = np.ones(n)
|
||||
|
||||
whichs = ['LM', 'SM', 'LR', 'SR', 'LI', 'SI']
|
||||
dtypes = [np.float32, np.float64]
|
||||
|
||||
for which, sigma, dtype in itertools.product(whichs, [None, 0, 5], dtypes):
|
||||
prev_w = np.array([], dtype=dtype)
|
||||
eps = np.finfo(dtype).eps
|
||||
for k in range(1, 9):
|
||||
w, z = eigs(A.astype(dtype), k=k, which=which, sigma=sigma,
|
||||
v0=v0.astype(dtype), tol=0)
|
||||
assert_allclose(np.linalg.norm(A.dot(z) - z * w), 0, atol=np.sqrt(eps))
|
||||
|
||||
# Check that the set of eigenvalues for `k` is a subset of that for `k+1`
|
||||
dist = abs(prev_w[:,None] - w).min(axis=1)
|
||||
assert_allclose(dist, 0, atol=np.sqrt(eps))
|
||||
|
||||
prev_w = w
|
||||
|
||||
# Check sort order
|
||||
if sigma is None:
|
||||
d = w
|
||||
else:
|
||||
d = 1 / (w - sigma)
|
||||
|
||||
if which == 'LM':
|
||||
# ARPACK is systematic for 'LM', but sort order
|
||||
# appears not well defined for other modes
|
||||
assert np.all(np.diff(abs(d)) <= 1e-6)
|
||||
@@ -1,16 +0,0 @@
|
||||
"""
|
||||
Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG)
|
||||
|
||||
LOBPCG is a preconditioned eigensolver for large symmetric positive definite
|
||||
(SPD) generalized eigenproblems.
|
||||
|
||||
Call the function lobpcg - see help for lobpcg.lobpcg.
|
||||
|
||||
"""
|
||||
from .lobpcg import *
|
||||
|
||||
__all__ = [s for s in dir() if not s.startswith('_')]
|
||||
|
||||
from scipy._lib._testutils import PytestTester
|
||||
test = PytestTester(__name__)
|
||||
del PytestTester
|
||||
Binary file not shown.
Binary file not shown.
@@ -1,982 +0,0 @@
|
||||
"""
|
||||
Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] A. V. Knyazev (2001),
|
||||
Toward the Optimal Preconditioned Eigensolver: Locally Optimal
|
||||
Block Preconditioned Conjugate Gradient Method.
|
||||
SIAM Journal on Scientific Computing 23, no. 2,
|
||||
pp. 517-541. :doi:`10.1137/S1064827500366124`
|
||||
|
||||
.. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov (2007),
|
||||
Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX)
|
||||
in hypre and PETSc. :arxiv:`0705.2626`
|
||||
|
||||
.. [3] A. V. Knyazev's C and MATLAB implementations:
|
||||
https://github.com/lobpcg/blopex
|
||||
"""
|
||||
|
||||
import warnings
|
||||
import numpy as np
|
||||
from scipy.linalg import (inv, eigh, cho_factor, cho_solve,
|
||||
cholesky, LinAlgError)
|
||||
from scipy.sparse.linalg import LinearOperator
|
||||
from scipy.sparse import isspmatrix
|
||||
from numpy import block as bmat
|
||||
|
||||
__all__ = ["lobpcg"]
|
||||
|
||||
|
||||
def _report_nonhermitian(M, name):
|
||||
"""
|
||||
Report if `M` is not a Hermitian matrix given its type.
|
||||
"""
|
||||
from scipy.linalg import norm
|
||||
|
||||
md = M - M.T.conj()
|
||||
nmd = norm(md, 1)
|
||||
tol = 10 * np.finfo(M.dtype).eps
|
||||
tol = max(tol, tol * norm(M, 1))
|
||||
if nmd > tol:
|
||||
warnings.warn(
|
||||
f"Matrix {name} of the type {M.dtype} is not Hermitian: "
|
||||
f"condition: {nmd} < {tol} fails.",
|
||||
UserWarning, stacklevel=4
|
||||
)
|
||||
|
||||
def _as2d(ar):
|
||||
"""
|
||||
If the input array is 2D return it, if it is 1D, append a dimension,
|
||||
making it a column vector.
|
||||
"""
|
||||
if ar.ndim == 2:
|
||||
return ar
|
||||
else: # Assume 1!
|
||||
aux = np.array(ar, copy=False)
|
||||
aux.shape = (ar.shape[0], 1)
|
||||
return aux
|
||||
|
||||
|
||||
def _makeMatMat(m):
|
||||
if m is None:
|
||||
return None
|
||||
elif callable(m):
|
||||
return lambda v: m(v)
|
||||
else:
|
||||
return lambda v: m @ v
|
||||
|
||||
|
||||
def _applyConstraints(blockVectorV, factYBY, blockVectorBY, blockVectorY):
|
||||
"""Changes blockVectorV in place."""
|
||||
YBV = np.dot(blockVectorBY.T.conj(), blockVectorV)
|
||||
tmp = cho_solve(factYBY, YBV)
|
||||
blockVectorV -= np.dot(blockVectorY, tmp)
|
||||
|
||||
|
||||
def _b_orthonormalize(B, blockVectorV, blockVectorBV=None,
|
||||
verbosityLevel=0):
|
||||
"""in-place B-orthonormalize the given block vector using Cholesky."""
|
||||
normalization = blockVectorV.max(axis=0) + np.finfo(blockVectorV.dtype).eps
|
||||
blockVectorV = blockVectorV / normalization
|
||||
if blockVectorBV is None:
|
||||
if B is not None:
|
||||
try:
|
||||
blockVectorBV = B(blockVectorV)
|
||||
except Exception as e:
|
||||
if verbosityLevel:
|
||||
warnings.warn(
|
||||
f"Secondary MatMul call failed with error\n"
|
||||
f"{e}\n",
|
||||
UserWarning, stacklevel=3
|
||||
)
|
||||
return None, None, None, normalization
|
||||
if blockVectorBV.shape != blockVectorV.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorV.shape} "
|
||||
f"of the orthogonalized matrix not preserved\n"
|
||||
f"and changed to {blockVectorBV.shape} "
|
||||
f"after multiplying by the secondary matrix.\n"
|
||||
)
|
||||
else:
|
||||
blockVectorBV = blockVectorV # Shared data!!!
|
||||
else:
|
||||
blockVectorBV = blockVectorBV / normalization
|
||||
VBV = blockVectorV.T.conj() @ blockVectorBV
|
||||
try:
|
||||
# VBV is a Cholesky factor from now on...
|
||||
VBV = cholesky(VBV, overwrite_a=True)
|
||||
VBV = inv(VBV, overwrite_a=True)
|
||||
blockVectorV = blockVectorV @ VBV
|
||||
# blockVectorV = (cho_solve((VBV.T, True), blockVectorV.T)).T
|
||||
if B is not None:
|
||||
blockVectorBV = blockVectorBV @ VBV
|
||||
# blockVectorBV = (cho_solve((VBV.T, True), blockVectorBV.T)).T
|
||||
return blockVectorV, blockVectorBV, VBV, normalization
|
||||
except LinAlgError:
|
||||
if verbosityLevel:
|
||||
warnings.warn(
|
||||
"Cholesky has failed.",
|
||||
UserWarning, stacklevel=3
|
||||
)
|
||||
return None, None, None, normalization
|
||||
|
||||
|
||||
def _get_indx(_lambda, num, largest):
|
||||
"""Get `num` indices into `_lambda` depending on `largest` option."""
|
||||
ii = np.argsort(_lambda)
|
||||
if largest:
|
||||
ii = ii[:-num - 1:-1]
|
||||
else:
|
||||
ii = ii[:num]
|
||||
|
||||
return ii
|
||||
|
||||
|
||||
def _handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel):
|
||||
if verbosityLevel:
|
||||
_report_nonhermitian(gramA, "gramA")
|
||||
_report_nonhermitian(gramB, "gramB")
|
||||
|
||||
|
||||
def lobpcg(
|
||||
A,
|
||||
X,
|
||||
B=None,
|
||||
M=None,
|
||||
Y=None,
|
||||
tol=None,
|
||||
maxiter=None,
|
||||
largest=True,
|
||||
verbosityLevel=0,
|
||||
retLambdaHistory=False,
|
||||
retResidualNormsHistory=False,
|
||||
restartControl=20,
|
||||
):
|
||||
"""Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).
|
||||
|
||||
LOBPCG is a preconditioned eigensolver for large symmetric positive
|
||||
definite (SPD) generalized eigenproblems.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : {sparse matrix, dense matrix, LinearOperator, callable object}
|
||||
The symmetric linear operator of the problem, usually a
|
||||
sparse matrix. Often called the "stiffness matrix".
|
||||
X : ndarray, float32 or float64
|
||||
Initial approximation to the ``k`` eigenvectors (non-sparse). If `A`
|
||||
has ``shape=(n,n)`` then `X` should have shape ``shape=(n,k)``.
|
||||
B : {dense matrix, sparse matrix, LinearOperator, callable object}
|
||||
Optional.
|
||||
The right hand side operator in a generalized eigenproblem.
|
||||
By default, ``B = Identity``. Often called the "mass matrix".
|
||||
M : {dense matrix, sparse matrix, LinearOperator, callable object}
|
||||
Optional.
|
||||
Preconditioner to `A`; by default ``M = Identity``.
|
||||
`M` should approximate the inverse of `A`.
|
||||
Y : ndarray, float32 or float64, optional.
|
||||
An n-by-sizeY matrix of constraints (non-sparse), sizeY < n.
|
||||
The iterations will be performed in the B-orthogonal complement
|
||||
of the column-space of Y. Y must be full rank.
|
||||
tol : scalar, optional.
|
||||
Solver tolerance (stopping criterion).
|
||||
The default is ``tol=n*sqrt(eps)``.
|
||||
maxiter : int, optional.
|
||||
Maximum number of iterations. The default is ``maxiter=20``.
|
||||
largest : bool, optional.
|
||||
When True, solve for the largest eigenvalues, otherwise the smallest.
|
||||
verbosityLevel : int, optional
|
||||
Controls solver output. The default is ``verbosityLevel=0``.
|
||||
retLambdaHistory : bool, optional.
|
||||
Whether to return eigenvalue history. Default is False.
|
||||
retResidualNormsHistory : bool, optional.
|
||||
Whether to return history of residual norms. Default is False.
|
||||
restartControl : int, optional.
|
||||
Iterations restart if the residuals jump up 2**restartControl times
|
||||
compared to the smallest ones recorded in retResidualNormsHistory.
|
||||
The default is ``restartControl=20``, making the restarts rare for
|
||||
backward compatibility.
|
||||
|
||||
Returns
|
||||
-------
|
||||
w : ndarray
|
||||
Array of ``k`` eigenvalues.
|
||||
v : ndarray
|
||||
An array of ``k`` eigenvectors. `v` has the same shape as `X`.
|
||||
lambdas : ndarray, optional
|
||||
The eigenvalue history, if `retLambdaHistory` is True.
|
||||
rnorms : ndarray, optional
|
||||
The history of residual norms, if `retResidualNormsHistory` is True.
|
||||
|
||||
Notes
|
||||
-----
|
||||
The iterative loop in lobpcg runs maxit=maxiter (or 20 if maxit=None)
|
||||
iterations at most and finishes earler if the tolerance is met.
|
||||
Breaking backward compatibility with the previous version, lobpcg
|
||||
now returns the block of iterative vectors with the best accuracy rather
|
||||
than the last one iterated, as a cure for possible divergence.
|
||||
|
||||
The size of the iteration history output equals to the number of the best
|
||||
(limited by maxit) iterations plus 3 (initial, final, and postprocessing).
|
||||
|
||||
If both ``retLambdaHistory`` and ``retResidualNormsHistory`` are True,
|
||||
the return tuple has the following format
|
||||
``(lambda, V, lambda history, residual norms history)``.
|
||||
|
||||
In the following ``n`` denotes the matrix size and ``k`` the number
|
||||
of required eigenvalues (smallest or largest).
|
||||
|
||||
The LOBPCG code internally solves eigenproblems of the size ``3k`` on every
|
||||
iteration by calling the dense eigensolver `eigh`, so if ``k`` is not
|
||||
small enough compared to ``n``, it makes no sense to call the LOBPCG code.
|
||||
Moreover, if one calls the LOBPCG algorithm for ``5k > n``, it would likely
|
||||
break internally, so the code calls the standard function `eigh` instead.
|
||||
It is not that ``n`` should be large for the LOBPCG to work, but rather the
|
||||
ratio ``n / k`` should be large. It you call LOBPCG with ``k=1``
|
||||
and ``n=10``, it works though ``n`` is small. The method is intended
|
||||
for extremely large ``n / k``.
|
||||
|
||||
The convergence speed depends basically on two factors:
|
||||
|
||||
1. Relative separation of the seeking eigenvalues from the rest
|
||||
of the eigenvalues. One can vary ``k`` to improve the absolute
|
||||
separation and use proper preconditioning to shrink the spectral spread.
|
||||
For example, a rod vibration test problem (under tests
|
||||
directory) is ill-conditioned for large ``n``, so convergence will be
|
||||
slow, unless efficient preconditioning is used. For this specific
|
||||
problem, a good simple preconditioner function would be a linear solve
|
||||
for `A`, which is easy to code since `A` is tridiagonal.
|
||||
|
||||
2. Quality of the initial approximations `X` to the seeking eigenvectors.
|
||||
Randomly distributed around the origin vectors work well if no better
|
||||
choice is known.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] A. V. Knyazev (2001),
|
||||
Toward the Optimal Preconditioned Eigensolver: Locally Optimal
|
||||
Block Preconditioned Conjugate Gradient Method.
|
||||
SIAM Journal on Scientific Computing 23, no. 2,
|
||||
pp. 517-541. :doi:`10.1137/S1064827500366124`
|
||||
|
||||
.. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov
|
||||
(2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers
|
||||
(BLOPEX) in hypre and PETSc. :arxiv:`0705.2626`
|
||||
|
||||
.. [3] A. V. Knyazev's C and MATLAB implementations:
|
||||
https://github.com/lobpcg/blopex
|
||||
|
||||
Examples
|
||||
--------
|
||||
Solve ``A x = lambda x`` with constraints and preconditioning.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy.sparse import spdiags, issparse
|
||||
>>> from scipy.sparse.linalg import lobpcg, LinearOperator
|
||||
|
||||
The square matrix size:
|
||||
|
||||
>>> n = 100
|
||||
>>> vals = np.arange(1, n + 1)
|
||||
|
||||
The first mandatory input parameter, in this test
|
||||
a sparse 2D array representing the square matrix
|
||||
of the eigenvalue problem to solve:
|
||||
|
||||
>>> A = spdiags(vals, 0, n, n)
|
||||
>>> A.toarray()
|
||||
array([[ 1, 0, 0, ..., 0, 0, 0],
|
||||
[ 0, 2, 0, ..., 0, 0, 0],
|
||||
[ 0, 0, 3, ..., 0, 0, 0],
|
||||
...,
|
||||
[ 0, 0, 0, ..., 98, 0, 0],
|
||||
[ 0, 0, 0, ..., 0, 99, 0],
|
||||
[ 0, 0, 0, ..., 0, 0, 100]])
|
||||
|
||||
Initial guess for eigenvectors, should have linearly independent
|
||||
columns. The second mandatory input parameter, a 2D array with the
|
||||
row dimension determining the number of requested eigenvalues.
|
||||
If no initial approximations available, randomly oriented vectors
|
||||
commonly work best, e.g., with components normally disrtibuted
|
||||
around zero or uniformly distributed on the interval [-1 1].
|
||||
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> X = rng.normal(size=(n, 3))
|
||||
|
||||
Constraints - an optional input parameter is a 2D array comprising
|
||||
of column vectors that the eigenvectors must be orthogonal to:
|
||||
|
||||
>>> Y = np.eye(n, 3)
|
||||
|
||||
Preconditioner in the inverse of A in this example:
|
||||
|
||||
>>> invA = spdiags([1./vals], 0, n, n)
|
||||
|
||||
The preconditiner must be defined by a function:
|
||||
|
||||
>>> def precond( x ):
|
||||
... return invA @ x
|
||||
|
||||
The argument x of the preconditioner function is a matrix inside `lobpcg`,
|
||||
thus the use of matrix-matrix product ``@``.
|
||||
|
||||
The preconditioner function is passed to lobpcg as a `LinearOperator`:
|
||||
|
||||
>>> M = LinearOperator(matvec=precond, matmat=precond,
|
||||
... shape=(n, n), dtype=np.float64)
|
||||
|
||||
Let us now solve the eigenvalue problem for the matrix A:
|
||||
|
||||
>>> eigenvalues, _ = lobpcg(A, X, Y=Y, M=M, largest=False)
|
||||
>>> eigenvalues
|
||||
array([4., 5., 6.])
|
||||
|
||||
Note that the vectors passed in Y are the eigenvectors of the 3 smallest
|
||||
eigenvalues. The results returned are orthogonal to those.
|
||||
"""
|
||||
blockVectorX = X
|
||||
bestblockVectorX = blockVectorX
|
||||
blockVectorY = Y
|
||||
residualTolerance = tol
|
||||
if maxiter is None:
|
||||
maxiter = 20
|
||||
|
||||
bestIterationNumber = maxiter
|
||||
|
||||
sizeY = 0
|
||||
if blockVectorY is not None:
|
||||
if len(blockVectorY.shape) != 2:
|
||||
warnings.warn(
|
||||
f"Expected rank-2 array for argument Y, instead got "
|
||||
f"{len(blockVectorY.shape)}, "
|
||||
f"so ignore it and use no constraints.",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
blockVectorY = None
|
||||
else:
|
||||
sizeY = blockVectorY.shape[1]
|
||||
|
||||
# Block size.
|
||||
if blockVectorX is None:
|
||||
raise ValueError("The mandatory initial matrix X cannot be None")
|
||||
if len(blockVectorX.shape) != 2:
|
||||
raise ValueError("expected rank-2 array for argument X")
|
||||
|
||||
n, sizeX = blockVectorX.shape
|
||||
|
||||
# Data type of iterates, determined by X, must be inexact
|
||||
if not np.issubdtype(blockVectorX.dtype, np.inexact):
|
||||
warnings.warn(
|
||||
f"Data type for argument X is {blockVectorX.dtype}, "
|
||||
f"which is not inexact, so casted to np.float32.",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
blockVectorX = np.asarray(blockVectorX, dtype=np.float32)
|
||||
|
||||
if retLambdaHistory:
|
||||
lambdaHistory = np.zeros((maxiter + 3, sizeX),
|
||||
dtype=blockVectorX.dtype)
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory = np.zeros((maxiter + 3, sizeX),
|
||||
dtype=blockVectorX.dtype)
|
||||
|
||||
if verbosityLevel:
|
||||
aux = "Solving "
|
||||
if B is None:
|
||||
aux += "standard"
|
||||
else:
|
||||
aux += "generalized"
|
||||
aux += " eigenvalue problem with"
|
||||
if M is None:
|
||||
aux += "out"
|
||||
aux += " preconditioning\n\n"
|
||||
aux += "matrix size %d\n" % n
|
||||
aux += "block size %d\n\n" % sizeX
|
||||
if blockVectorY is None:
|
||||
aux += "No constraints\n\n"
|
||||
else:
|
||||
if sizeY > 1:
|
||||
aux += "%d constraints\n\n" % sizeY
|
||||
else:
|
||||
aux += "%d constraint\n\n" % sizeY
|
||||
print(aux)
|
||||
|
||||
if (n - sizeY) < (5 * sizeX):
|
||||
warnings.warn(
|
||||
f"The problem size {n} minus the constraints size {sizeY} "
|
||||
f"is too small relative to the block size {sizeX}. "
|
||||
f"Using a dense eigensolver instead of LOBPCG iterations."
|
||||
f"No output of the history of the iterations.",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
|
||||
sizeX = min(sizeX, n)
|
||||
|
||||
if blockVectorY is not None:
|
||||
raise NotImplementedError(
|
||||
"The dense eigensolver does not support constraints."
|
||||
)
|
||||
|
||||
# Define the closed range of indices of eigenvalues to return.
|
||||
if largest:
|
||||
eigvals = (n - sizeX, n - 1)
|
||||
else:
|
||||
eigvals = (0, sizeX - 1)
|
||||
|
||||
try:
|
||||
if isinstance(A, LinearOperator):
|
||||
A = A(np.eye(n, dtype=int))
|
||||
elif callable(A):
|
||||
A = A(np.eye(n, dtype=int))
|
||||
if A.shape != (n, n):
|
||||
raise ValueError(
|
||||
f"The shape {A.shape} of the primary matrix\n"
|
||||
f"defined by a callable object is wrong.\n"
|
||||
)
|
||||
elif isspmatrix(A):
|
||||
A = A.toarray()
|
||||
else:
|
||||
A = np.asarray(A)
|
||||
except Exception as e:
|
||||
raise Exception(
|
||||
f"Primary MatMul call failed with error\n"
|
||||
f"{e}\n")
|
||||
|
||||
if B is not None:
|
||||
try:
|
||||
if isinstance(B, LinearOperator):
|
||||
B = B(np.eye(n, dtype=int))
|
||||
elif callable(B):
|
||||
B = B(np.eye(n, dtype=int))
|
||||
if B.shape != (n, n):
|
||||
raise ValueError(
|
||||
f"The shape {B.shape} of the secondary matrix\n"
|
||||
f"defined by a callable object is wrong.\n"
|
||||
)
|
||||
elif isspmatrix(B):
|
||||
B = B.toarray()
|
||||
else:
|
||||
B = np.asarray(B)
|
||||
except Exception as e:
|
||||
raise Exception(
|
||||
f"Secondary MatMul call failed with error\n"
|
||||
f"{e}\n")
|
||||
|
||||
try:
|
||||
vals, vecs = eigh(A,
|
||||
B,
|
||||
subset_by_index=eigvals,
|
||||
check_finite=False)
|
||||
if largest:
|
||||
# Reverse order to be compatible with eigs() in 'LM' mode.
|
||||
vals = vals[::-1]
|
||||
vecs = vecs[:, ::-1]
|
||||
|
||||
return vals, vecs
|
||||
except Exception as e:
|
||||
raise Exception(
|
||||
f"Dense eigensolver failed with error\n"
|
||||
f"{e}\n"
|
||||
)
|
||||
|
||||
if (residualTolerance is None) or (residualTolerance <= 0.0):
|
||||
residualTolerance = np.sqrt(np.finfo(blockVectorX.dtype).eps) * n
|
||||
|
||||
A = _makeMatMat(A)
|
||||
B = _makeMatMat(B)
|
||||
M = _makeMatMat(M)
|
||||
|
||||
# Apply constraints to X.
|
||||
if blockVectorY is not None:
|
||||
|
||||
if B is not None:
|
||||
blockVectorBY = B(blockVectorY)
|
||||
if blockVectorBY.shape != blockVectorY.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorY.shape} "
|
||||
f"of the constraint not preserved\n"
|
||||
f"and changed to {blockVectorBY.shape} "
|
||||
f"after multiplying by the secondary matrix.\n"
|
||||
)
|
||||
else:
|
||||
blockVectorBY = blockVectorY
|
||||
|
||||
# gramYBY is a dense array.
|
||||
gramYBY = np.dot(blockVectorY.T.conj(), blockVectorBY)
|
||||
try:
|
||||
# gramYBY is a Cholesky factor from now on...
|
||||
gramYBY = cho_factor(gramYBY)
|
||||
except LinAlgError as e:
|
||||
raise ValueError("Linearly dependent constraints") from e
|
||||
|
||||
_applyConstraints(blockVectorX, gramYBY, blockVectorBY, blockVectorY)
|
||||
|
||||
##
|
||||
# B-orthonormalize X.
|
||||
blockVectorX, blockVectorBX, _, _ = _b_orthonormalize(
|
||||
B, blockVectorX, verbosityLevel=verbosityLevel)
|
||||
if blockVectorX is None:
|
||||
raise ValueError("Linearly dependent initial approximations")
|
||||
|
||||
##
|
||||
# Compute the initial Ritz vectors: solve the eigenproblem.
|
||||
blockVectorAX = A(blockVectorX)
|
||||
if blockVectorAX.shape != blockVectorX.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorX.shape} "
|
||||
f"of the initial approximations not preserved\n"
|
||||
f"and changed to {blockVectorAX.shape} "
|
||||
f"after multiplying by the primary matrix.\n"
|
||||
)
|
||||
|
||||
gramXAX = np.dot(blockVectorX.T.conj(), blockVectorAX)
|
||||
|
||||
_lambda, eigBlockVector = eigh(gramXAX, check_finite=False)
|
||||
ii = _get_indx(_lambda, sizeX, largest)
|
||||
_lambda = _lambda[ii]
|
||||
if retLambdaHistory:
|
||||
lambdaHistory[0, :] = _lambda
|
||||
|
||||
eigBlockVector = np.asarray(eigBlockVector[:, ii])
|
||||
blockVectorX = np.dot(blockVectorX, eigBlockVector)
|
||||
blockVectorAX = np.dot(blockVectorAX, eigBlockVector)
|
||||
if B is not None:
|
||||
blockVectorBX = np.dot(blockVectorBX, eigBlockVector)
|
||||
|
||||
##
|
||||
# Active index set.
|
||||
activeMask = np.ones((sizeX,), dtype=bool)
|
||||
|
||||
##
|
||||
# Main iteration loop.
|
||||
|
||||
blockVectorP = None # set during iteration
|
||||
blockVectorAP = None
|
||||
blockVectorBP = None
|
||||
|
||||
smallestResidualNorm = np.abs(np.finfo(blockVectorX.dtype).max)
|
||||
|
||||
iterationNumber = -1
|
||||
restart = True
|
||||
forcedRestart = False
|
||||
explicitGramFlag = False
|
||||
while iterationNumber < maxiter:
|
||||
iterationNumber += 1
|
||||
|
||||
if B is not None:
|
||||
aux = blockVectorBX * _lambda[np.newaxis, :]
|
||||
else:
|
||||
aux = blockVectorX * _lambda[np.newaxis, :]
|
||||
|
||||
blockVectorR = blockVectorAX - aux
|
||||
|
||||
aux = np.sum(blockVectorR.conj() * blockVectorR, 0)
|
||||
residualNorms = np.sqrt(np.abs(aux))
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory[iterationNumber, :] = residualNorms
|
||||
residualNorm = np.sum(np.abs(residualNorms)) / sizeX
|
||||
|
||||
if residualNorm < smallestResidualNorm:
|
||||
smallestResidualNorm = residualNorm
|
||||
bestIterationNumber = iterationNumber
|
||||
bestblockVectorX = blockVectorX
|
||||
elif residualNorm > 2**restartControl * smallestResidualNorm:
|
||||
forcedRestart = True
|
||||
blockVectorAX = A(blockVectorX)
|
||||
if blockVectorAX.shape != blockVectorX.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorX.shape} "
|
||||
f"of the restarted iterate not preserved\n"
|
||||
f"and changed to {blockVectorAX.shape} "
|
||||
f"after multiplying by the primary matrix.\n"
|
||||
)
|
||||
if B is not None:
|
||||
blockVectorBX = B(blockVectorX)
|
||||
if blockVectorBX.shape != blockVectorX.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorX.shape} "
|
||||
f"of the restarted iterate not preserved\n"
|
||||
f"and changed to {blockVectorBX.shape} "
|
||||
f"after multiplying by the secondary matrix.\n"
|
||||
)
|
||||
|
||||
ii = np.where(residualNorms > residualTolerance, True, False)
|
||||
activeMask = activeMask & ii
|
||||
currentBlockSize = activeMask.sum()
|
||||
|
||||
if verbosityLevel:
|
||||
print(f"iteration {iterationNumber}")
|
||||
print(f"current block size: {currentBlockSize}")
|
||||
print(f"eigenvalue(s):\n{_lambda}")
|
||||
print(f"residual norm(s):\n{residualNorms}")
|
||||
|
||||
if currentBlockSize == 0:
|
||||
break
|
||||
|
||||
activeBlockVectorR = _as2d(blockVectorR[:, activeMask])
|
||||
|
||||
if iterationNumber > 0:
|
||||
activeBlockVectorP = _as2d(blockVectorP[:, activeMask])
|
||||
activeBlockVectorAP = _as2d(blockVectorAP[:, activeMask])
|
||||
if B is not None:
|
||||
activeBlockVectorBP = _as2d(blockVectorBP[:, activeMask])
|
||||
|
||||
if M is not None:
|
||||
# Apply preconditioner T to the active residuals.
|
||||
activeBlockVectorR = M(activeBlockVectorR)
|
||||
|
||||
##
|
||||
# Apply constraints to the preconditioned residuals.
|
||||
if blockVectorY is not None:
|
||||
_applyConstraints(activeBlockVectorR,
|
||||
gramYBY,
|
||||
blockVectorBY,
|
||||
blockVectorY)
|
||||
|
||||
##
|
||||
# B-orthogonalize the preconditioned residuals to X.
|
||||
if B is not None:
|
||||
activeBlockVectorR = activeBlockVectorR - (
|
||||
blockVectorX @
|
||||
(blockVectorBX.T.conj() @ activeBlockVectorR)
|
||||
)
|
||||
else:
|
||||
activeBlockVectorR = activeBlockVectorR - (
|
||||
blockVectorX @
|
||||
(blockVectorX.T.conj() @ activeBlockVectorR)
|
||||
)
|
||||
|
||||
##
|
||||
# B-orthonormalize the preconditioned residuals.
|
||||
aux = _b_orthonormalize(
|
||||
B, activeBlockVectorR, verbosityLevel=verbosityLevel)
|
||||
activeBlockVectorR, activeBlockVectorBR, _, _ = aux
|
||||
|
||||
if activeBlockVectorR is None:
|
||||
warnings.warn(
|
||||
f"Failed at iteration {iterationNumber} with accuracies "
|
||||
f"{residualNorms}\n not reaching the requested "
|
||||
f"tolerance {residualTolerance}.",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
break
|
||||
activeBlockVectorAR = A(activeBlockVectorR)
|
||||
|
||||
if iterationNumber > 0:
|
||||
if B is not None:
|
||||
aux = _b_orthonormalize(
|
||||
B, activeBlockVectorP, activeBlockVectorBP,
|
||||
verbosityLevel=verbosityLevel
|
||||
)
|
||||
activeBlockVectorP, activeBlockVectorBP, invR, normal = aux
|
||||
else:
|
||||
aux = _b_orthonormalize(B, activeBlockVectorP,
|
||||
verbosityLevel=verbosityLevel)
|
||||
activeBlockVectorP, _, invR, normal = aux
|
||||
# Function _b_orthonormalize returns None if Cholesky fails
|
||||
if activeBlockVectorP is not None:
|
||||
activeBlockVectorAP = activeBlockVectorAP / normal
|
||||
activeBlockVectorAP = np.dot(activeBlockVectorAP, invR)
|
||||
restart = forcedRestart
|
||||
else:
|
||||
restart = True
|
||||
|
||||
##
|
||||
# Perform the Rayleigh Ritz Procedure:
|
||||
# Compute symmetric Gram matrices:
|
||||
|
||||
if activeBlockVectorAR.dtype == "float32":
|
||||
myeps = 1
|
||||
else:
|
||||
myeps = np.sqrt(np.finfo(activeBlockVectorR.dtype).eps)
|
||||
|
||||
if residualNorms.max() > myeps and not explicitGramFlag:
|
||||
explicitGramFlag = False
|
||||
else:
|
||||
# Once explicitGramFlag, forever explicitGramFlag.
|
||||
explicitGramFlag = True
|
||||
|
||||
# Shared memory assingments to simplify the code
|
||||
if B is None:
|
||||
blockVectorBX = blockVectorX
|
||||
activeBlockVectorBR = activeBlockVectorR
|
||||
if not restart:
|
||||
activeBlockVectorBP = activeBlockVectorP
|
||||
|
||||
# Common submatrices:
|
||||
gramXAR = np.dot(blockVectorX.T.conj(), activeBlockVectorAR)
|
||||
gramRAR = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorAR)
|
||||
|
||||
gramDtype = activeBlockVectorAR.dtype
|
||||
if explicitGramFlag:
|
||||
gramRAR = (gramRAR + gramRAR.T.conj()) / 2
|
||||
gramXAX = np.dot(blockVectorX.T.conj(), blockVectorAX)
|
||||
gramXAX = (gramXAX + gramXAX.T.conj()) / 2
|
||||
gramXBX = np.dot(blockVectorX.T.conj(), blockVectorBX)
|
||||
gramRBR = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorBR)
|
||||
gramXBR = np.dot(blockVectorX.T.conj(), activeBlockVectorBR)
|
||||
else:
|
||||
gramXAX = np.diag(_lambda).astype(gramDtype)
|
||||
gramXBX = np.eye(sizeX, dtype=gramDtype)
|
||||
gramRBR = np.eye(currentBlockSize, dtype=gramDtype)
|
||||
gramXBR = np.zeros((sizeX, currentBlockSize), dtype=gramDtype)
|
||||
|
||||
if not restart:
|
||||
gramXAP = np.dot(blockVectorX.T.conj(), activeBlockVectorAP)
|
||||
gramRAP = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorAP)
|
||||
gramPAP = np.dot(activeBlockVectorP.T.conj(), activeBlockVectorAP)
|
||||
gramXBP = np.dot(blockVectorX.T.conj(), activeBlockVectorBP)
|
||||
gramRBP = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorBP)
|
||||
if explicitGramFlag:
|
||||
gramPAP = (gramPAP + gramPAP.T.conj()) / 2
|
||||
gramPBP = np.dot(activeBlockVectorP.T.conj(),
|
||||
activeBlockVectorBP)
|
||||
else:
|
||||
gramPBP = np.eye(currentBlockSize, dtype=gramDtype)
|
||||
|
||||
gramA = bmat(
|
||||
[
|
||||
[gramXAX, gramXAR, gramXAP],
|
||||
[gramXAR.T.conj(), gramRAR, gramRAP],
|
||||
[gramXAP.T.conj(), gramRAP.T.conj(), gramPAP],
|
||||
]
|
||||
)
|
||||
gramB = bmat(
|
||||
[
|
||||
[gramXBX, gramXBR, gramXBP],
|
||||
[gramXBR.T.conj(), gramRBR, gramRBP],
|
||||
[gramXBP.T.conj(), gramRBP.T.conj(), gramPBP],
|
||||
]
|
||||
)
|
||||
|
||||
_handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel)
|
||||
|
||||
try:
|
||||
_lambda, eigBlockVector = eigh(gramA,
|
||||
gramB,
|
||||
check_finite=False)
|
||||
except LinAlgError as e:
|
||||
# raise ValueError("eigh failed in lobpcg iterations") from e
|
||||
if verbosityLevel:
|
||||
warnings.warn(
|
||||
f"eigh failed at iteration {iterationNumber} \n"
|
||||
f"with error {e} causing a restart.\n",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
# try again after dropping the direction vectors P from RR
|
||||
restart = True
|
||||
|
||||
if restart:
|
||||
gramA = bmat([[gramXAX, gramXAR], [gramXAR.T.conj(), gramRAR]])
|
||||
gramB = bmat([[gramXBX, gramXBR], [gramXBR.T.conj(), gramRBR]])
|
||||
|
||||
_handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel)
|
||||
|
||||
try:
|
||||
_lambda, eigBlockVector = eigh(gramA,
|
||||
gramB,
|
||||
check_finite=False)
|
||||
except LinAlgError as e:
|
||||
# raise ValueError("eigh failed in lobpcg iterations") from e
|
||||
warnings.warn(
|
||||
f"eigh failed at iteration {iterationNumber} with error\n"
|
||||
f"{e}\n",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
break
|
||||
|
||||
ii = _get_indx(_lambda, sizeX, largest)
|
||||
_lambda = _lambda[ii]
|
||||
eigBlockVector = eigBlockVector[:, ii]
|
||||
if retLambdaHistory:
|
||||
lambdaHistory[iterationNumber + 1, :] = _lambda
|
||||
|
||||
# Compute Ritz vectors.
|
||||
if B is not None:
|
||||
if not restart:
|
||||
eigBlockVectorX = eigBlockVector[:sizeX]
|
||||
eigBlockVectorR = eigBlockVector[sizeX:
|
||||
sizeX + currentBlockSize]
|
||||
eigBlockVectorP = eigBlockVector[sizeX + currentBlockSize:]
|
||||
|
||||
pp = np.dot(activeBlockVectorR, eigBlockVectorR)
|
||||
pp += np.dot(activeBlockVectorP, eigBlockVectorP)
|
||||
|
||||
app = np.dot(activeBlockVectorAR, eigBlockVectorR)
|
||||
app += np.dot(activeBlockVectorAP, eigBlockVectorP)
|
||||
|
||||
bpp = np.dot(activeBlockVectorBR, eigBlockVectorR)
|
||||
bpp += np.dot(activeBlockVectorBP, eigBlockVectorP)
|
||||
else:
|
||||
eigBlockVectorX = eigBlockVector[:sizeX]
|
||||
eigBlockVectorR = eigBlockVector[sizeX:]
|
||||
|
||||
pp = np.dot(activeBlockVectorR, eigBlockVectorR)
|
||||
app = np.dot(activeBlockVectorAR, eigBlockVectorR)
|
||||
bpp = np.dot(activeBlockVectorBR, eigBlockVectorR)
|
||||
|
||||
blockVectorX = np.dot(blockVectorX, eigBlockVectorX) + pp
|
||||
blockVectorAX = np.dot(blockVectorAX, eigBlockVectorX) + app
|
||||
blockVectorBX = np.dot(blockVectorBX, eigBlockVectorX) + bpp
|
||||
|
||||
blockVectorP, blockVectorAP, blockVectorBP = pp, app, bpp
|
||||
|
||||
else:
|
||||
if not restart:
|
||||
eigBlockVectorX = eigBlockVector[:sizeX]
|
||||
eigBlockVectorR = eigBlockVector[sizeX:
|
||||
sizeX + currentBlockSize]
|
||||
eigBlockVectorP = eigBlockVector[sizeX + currentBlockSize:]
|
||||
|
||||
pp = np.dot(activeBlockVectorR, eigBlockVectorR)
|
||||
pp += np.dot(activeBlockVectorP, eigBlockVectorP)
|
||||
|
||||
app = np.dot(activeBlockVectorAR, eigBlockVectorR)
|
||||
app += np.dot(activeBlockVectorAP, eigBlockVectorP)
|
||||
else:
|
||||
eigBlockVectorX = eigBlockVector[:sizeX]
|
||||
eigBlockVectorR = eigBlockVector[sizeX:]
|
||||
|
||||
pp = np.dot(activeBlockVectorR, eigBlockVectorR)
|
||||
app = np.dot(activeBlockVectorAR, eigBlockVectorR)
|
||||
|
||||
blockVectorX = np.dot(blockVectorX, eigBlockVectorX) + pp
|
||||
blockVectorAX = np.dot(blockVectorAX, eigBlockVectorX) + app
|
||||
|
||||
blockVectorP, blockVectorAP = pp, app
|
||||
|
||||
if B is not None:
|
||||
aux = blockVectorBX * _lambda[np.newaxis, :]
|
||||
else:
|
||||
aux = blockVectorX * _lambda[np.newaxis, :]
|
||||
|
||||
blockVectorR = blockVectorAX - aux
|
||||
|
||||
aux = np.sum(blockVectorR.conj() * blockVectorR, 0)
|
||||
residualNorms = np.sqrt(np.abs(aux))
|
||||
# Use old lambda in case of early loop exit.
|
||||
if retLambdaHistory:
|
||||
lambdaHistory[iterationNumber + 1, :] = _lambda
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory[iterationNumber + 1, :] = residualNorms
|
||||
residualNorm = np.sum(np.abs(residualNorms)) / sizeX
|
||||
if residualNorm < smallestResidualNorm:
|
||||
smallestResidualNorm = residualNorm
|
||||
bestIterationNumber = iterationNumber + 1
|
||||
bestblockVectorX = blockVectorX
|
||||
|
||||
if np.max(np.abs(residualNorms)) > residualTolerance:
|
||||
warnings.warn(
|
||||
f"Exited at iteration {iterationNumber} with accuracies \n"
|
||||
f"{residualNorms}\n"
|
||||
f"not reaching the requested tolerance {residualTolerance}.\n"
|
||||
f"Use iteration {bestIterationNumber} instead with accuracy \n"
|
||||
f"{smallestResidualNorm}.\n",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
|
||||
if verbosityLevel:
|
||||
print(f"Final iterative eigenvalue(s):\n{_lambda}")
|
||||
print(f"Final iterative residual norm(s):\n{residualNorms}")
|
||||
|
||||
blockVectorX = bestblockVectorX
|
||||
# Making eigenvectors "exactly" satisfy the blockVectorY constrains
|
||||
if blockVectorY is not None:
|
||||
_applyConstraints(blockVectorX,
|
||||
gramYBY,
|
||||
blockVectorBY,
|
||||
blockVectorY)
|
||||
|
||||
# Making eigenvectors "exactly" othonormalized by final "exact" RR
|
||||
blockVectorAX = A(blockVectorX)
|
||||
if blockVectorAX.shape != blockVectorX.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorX.shape} "
|
||||
f"of the postprocessing iterate not preserved\n"
|
||||
f"and changed to {blockVectorAX.shape} "
|
||||
f"after multiplying by the primary matrix.\n"
|
||||
)
|
||||
gramXAX = np.dot(blockVectorX.T.conj(), blockVectorAX)
|
||||
|
||||
blockVectorBX = blockVectorX
|
||||
if B is not None:
|
||||
blockVectorBX = B(blockVectorX)
|
||||
if blockVectorBX.shape != blockVectorX.shape:
|
||||
raise ValueError(
|
||||
f"The shape {blockVectorX.shape} "
|
||||
f"of the postprocessing iterate not preserved\n"
|
||||
f"and changed to {blockVectorBX.shape} "
|
||||
f"after multiplying by the secondary matrix.\n"
|
||||
)
|
||||
|
||||
gramXBX = np.dot(blockVectorX.T.conj(), blockVectorBX)
|
||||
_handle_gramA_gramB_verbosity(gramXAX, gramXBX, verbosityLevel)
|
||||
gramXAX = (gramXAX + gramXAX.T.conj()) / 2
|
||||
gramXBX = (gramXBX + gramXBX.T.conj()) / 2
|
||||
try:
|
||||
_lambda, eigBlockVector = eigh(gramXAX,
|
||||
gramXBX,
|
||||
check_finite=False)
|
||||
except LinAlgError as e:
|
||||
raise ValueError("eigh has failed in lobpcg postprocessing") from e
|
||||
|
||||
ii = _get_indx(_lambda, sizeX, largest)
|
||||
_lambda = _lambda[ii]
|
||||
eigBlockVector = np.asarray(eigBlockVector[:, ii])
|
||||
|
||||
blockVectorX = np.dot(blockVectorX, eigBlockVector)
|
||||
blockVectorAX = np.dot(blockVectorAX, eigBlockVector)
|
||||
|
||||
if B is not None:
|
||||
blockVectorBX = np.dot(blockVectorBX, eigBlockVector)
|
||||
aux = blockVectorBX * _lambda[np.newaxis, :]
|
||||
else:
|
||||
aux = blockVectorX * _lambda[np.newaxis, :]
|
||||
|
||||
blockVectorR = blockVectorAX - aux
|
||||
|
||||
aux = np.sum(blockVectorR.conj() * blockVectorR, 0)
|
||||
residualNorms = np.sqrt(np.abs(aux))
|
||||
|
||||
if retLambdaHistory:
|
||||
lambdaHistory[bestIterationNumber + 1, :] = _lambda
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory[bestIterationNumber + 1, :] = residualNorms
|
||||
|
||||
if retLambdaHistory:
|
||||
lambdaHistory = lambdaHistory[
|
||||
: bestIterationNumber + 2, :]
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory = residualNormsHistory[
|
||||
: bestIterationNumber + 2, :]
|
||||
|
||||
if np.max(np.abs(residualNorms)) > residualTolerance:
|
||||
warnings.warn(
|
||||
f"Exited postprocessing with accuracies \n"
|
||||
f"{residualNorms}\n"
|
||||
f"not reaching the requested tolerance {residualTolerance}.",
|
||||
UserWarning, stacklevel=2
|
||||
)
|
||||
|
||||
if verbosityLevel:
|
||||
print(f"Final postprocessing eigenvalue(s):\n{_lambda}")
|
||||
print(f"Final residual norm(s):\n{residualNorms}")
|
||||
|
||||
if retLambdaHistory:
|
||||
lambdaHistory = np.vsplit(lambdaHistory, np.shape(lambdaHistory)[0])
|
||||
lambdaHistory = [np.squeeze(i) for i in lambdaHistory]
|
||||
if retResidualNormsHistory:
|
||||
residualNormsHistory = np.vsplit(residualNormsHistory,
|
||||
np.shape(residualNormsHistory)[0])
|
||||
residualNormsHistory = [np.squeeze(i) for i in residualNormsHistory]
|
||||
|
||||
if retLambdaHistory:
|
||||
if retResidualNormsHistory:
|
||||
return _lambda, blockVectorX, lambdaHistory, residualNormsHistory
|
||||
else:
|
||||
return _lambda, blockVectorX, lambdaHistory
|
||||
else:
|
||||
if retResidualNormsHistory:
|
||||
return _lambda, blockVectorX, residualNormsHistory
|
||||
else:
|
||||
return _lambda, blockVectorX
|
||||
Binary file not shown.
Binary file not shown.
@@ -1,534 +0,0 @@
|
||||
""" Test functions for the sparse.linalg._eigen.lobpcg module
|
||||
"""
|
||||
|
||||
import itertools
|
||||
import platform
|
||||
import sys
|
||||
import pytest
|
||||
import numpy as np
|
||||
from numpy import ones, r_, diag
|
||||
from numpy.testing import (assert_almost_equal, assert_equal,
|
||||
assert_allclose, assert_array_less)
|
||||
|
||||
from scipy.linalg import eig, eigh, toeplitz, orth
|
||||
from scipy.sparse import spdiags, diags, eye, csr_matrix
|
||||
from scipy.sparse.linalg import eigs, LinearOperator
|
||||
from scipy.sparse.linalg._eigen.lobpcg import lobpcg
|
||||
|
||||
_IS_32BIT = (sys.maxsize < 2**32)
|
||||
|
||||
|
||||
def ElasticRod(n):
|
||||
"""Build the matrices for the generalized eigenvalue problem of the
|
||||
fixed-free elastic rod vibration model.
|
||||
"""
|
||||
L = 1.0
|
||||
le = L/n
|
||||
rho = 7.85e3
|
||||
S = 1.e-4
|
||||
E = 2.1e11
|
||||
mass = rho*S*le/6.
|
||||
k = E*S/le
|
||||
A = k*(diag(r_[2.*ones(n-1), 1])-diag(ones(n-1), 1)-diag(ones(n-1), -1))
|
||||
B = mass*(diag(r_[4.*ones(n-1), 2])+diag(ones(n-1), 1)+diag(ones(n-1), -1))
|
||||
return A, B
|
||||
|
||||
|
||||
def MikotaPair(n):
|
||||
"""Build a pair of full diagonal matrices for the generalized eigenvalue
|
||||
problem. The Mikota pair acts as a nice test since the eigenvalues are the
|
||||
squares of the integers n, n=1,2,...
|
||||
"""
|
||||
x = np.arange(1, n+1)
|
||||
B = diag(1./x)
|
||||
y = np.arange(n-1, 0, -1)
|
||||
z = np.arange(2*n-1, 0, -2)
|
||||
A = diag(z)-diag(y, -1)-diag(y, 1)
|
||||
return A, B
|
||||
|
||||
|
||||
def compare_solutions(A, B, m):
|
||||
"""Check eig vs. lobpcg consistency.
|
||||
"""
|
||||
n = A.shape[0]
|
||||
rnd = np.random.RandomState(0)
|
||||
V = rnd.random((n, m))
|
||||
X = orth(V)
|
||||
eigvals, _ = lobpcg(A, X, B=B, tol=1e-2, maxiter=50, largest=False)
|
||||
eigvals.sort()
|
||||
w, _ = eig(A, b=B)
|
||||
w.sort()
|
||||
assert_almost_equal(w[:int(m/2)], eigvals[:int(m/2)], decimal=2)
|
||||
|
||||
|
||||
def test_Small():
|
||||
A, B = ElasticRod(10)
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
compare_solutions(A, B, 10)
|
||||
A, B = MikotaPair(10)
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
compare_solutions(A, B, 10)
|
||||
|
||||
|
||||
def test_ElasticRod():
|
||||
A, B = ElasticRod(20)
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
compare_solutions(A, B, 2)
|
||||
|
||||
|
||||
def test_MikotaPair():
|
||||
A, B = MikotaPair(20)
|
||||
compare_solutions(A, B, 2)
|
||||
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:Exited at iteration 0")
|
||||
@pytest.mark.filterwarnings("ignore:Exited postprocessing")
|
||||
def test_nonhermitian_warning(capsys):
|
||||
"""Check the warning of a Ritz matrix being not Hermitian
|
||||
by feeding a non-Hermitian input matrix.
|
||||
Also check stdout since verbosityLevel=1 and lack of stderr.
|
||||
"""
|
||||
n = 10
|
||||
X = np.arange(n * 2).reshape(n, 2).astype(np.float32)
|
||||
A = np.arange(n * n).reshape(n, n).astype(np.float32)
|
||||
with pytest.warns(UserWarning, match="Matrix gramA"):
|
||||
_, _ = lobpcg(A, X, verbosityLevel=1, maxiter=0)
|
||||
out, err = capsys.readouterr() # Capture output
|
||||
assert out.startswith("Solving standard eigenvalue") # Test stdout
|
||||
assert err == '' # Test empty stderr
|
||||
# Make the matrix symmetric and the UserWarning dissappears.
|
||||
A += A.T
|
||||
_, _ = lobpcg(A, X, verbosityLevel=1, maxiter=0)
|
||||
out, err = capsys.readouterr() # Capture output
|
||||
assert out.startswith("Solving standard eigenvalue") # Test stdout
|
||||
assert err == '' # Test empty stderr
|
||||
|
||||
|
||||
def test_regression():
|
||||
"""Check the eigenvalue of the identity matrix is one.
|
||||
"""
|
||||
# https://mail.python.org/pipermail/scipy-user/2010-October/026944.html
|
||||
n = 10
|
||||
X = np.ones((n, 1))
|
||||
A = np.identity(n)
|
||||
w, _ = lobpcg(A, X)
|
||||
assert_allclose(w, [1])
|
||||
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:The problem size")
|
||||
@pytest.mark.parametrize('n, m, m_excluded', [(100, 4, 3), (4, 2, 0)])
|
||||
def test_diagonal(n, m, m_excluded):
|
||||
"""Test ``m - m_excluded`` eigenvalues and eigenvectors of
|
||||
diagonal matrices of the size ``n`` varying matrix formats:
|
||||
dense array, spare matrix, and ``LinearOperator`` for both
|
||||
matrixes in the generalized eigenvalue problem ``Av = cBv``
|
||||
and for the preconditioner.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
|
||||
# Define the generalized eigenvalue problem Av = cBv
|
||||
# where (c, v) is a generalized eigenpair,
|
||||
# A is the diagonal matrix whose entries are 1,...n,
|
||||
# B is the identity matrix.
|
||||
vals = np.arange(1, n+1, dtype=float)
|
||||
A_s = diags([vals], [0], (n, n))
|
||||
A_a = A_s.toarray()
|
||||
|
||||
def A_f(x):
|
||||
return A_s @ x
|
||||
|
||||
A_lo = LinearOperator(matvec=A_f,
|
||||
matmat=A_f,
|
||||
shape=(n, n), dtype=float)
|
||||
|
||||
B_a = eye(n)
|
||||
B_s = csr_matrix(B_a)
|
||||
|
||||
def B_f(x):
|
||||
return B_a @ x
|
||||
|
||||
B_lo = LinearOperator(matvec=B_f,
|
||||
matmat=B_f,
|
||||
shape=(n, n), dtype=float)
|
||||
|
||||
# Let the preconditioner M be the inverse of A.
|
||||
M_s = diags([1./vals], [0], (n, n))
|
||||
M_a = M_s.toarray()
|
||||
|
||||
def M_f(x):
|
||||
return M_s @ x
|
||||
|
||||
M_lo = LinearOperator(matvec=M_f,
|
||||
matmat=M_f,
|
||||
shape=(n, n), dtype=float)
|
||||
|
||||
# Pick random initial vectors.
|
||||
X = rnd.normal(size=(n, m))
|
||||
|
||||
# Require that the returned eigenvectors be in the orthogonal complement
|
||||
# of the first few standard basis vectors.
|
||||
if m_excluded > 0:
|
||||
Y = np.eye(n, m_excluded)
|
||||
else:
|
||||
Y = None
|
||||
|
||||
for A in [A_a, A_s, A_lo]:
|
||||
for B in [B_a, B_s, B_lo]:
|
||||
for M in [M_a, M_s, M_lo]:
|
||||
eigvals, vecs = lobpcg(A, X, B, M=M, Y=Y,
|
||||
maxiter=40, largest=False)
|
||||
|
||||
assert_allclose(eigvals, np.arange(1+m_excluded,
|
||||
1+m_excluded+m))
|
||||
_check_eigen(A, eigvals, vecs, rtol=1e-3, atol=1e-3)
|
||||
|
||||
|
||||
def _check_eigen(M, w, V, rtol=1e-8, atol=1e-14):
|
||||
"""Check if the eigenvalue residual is small.
|
||||
"""
|
||||
mult_wV = np.multiply(w, V)
|
||||
dot_MV = M.dot(V)
|
||||
assert_allclose(mult_wV, dot_MV, rtol=rtol, atol=atol)
|
||||
|
||||
|
||||
def _check_fiedler(n, p):
|
||||
"""Check the Fiedler vector computation.
|
||||
"""
|
||||
# This is not necessarily the recommended way to find the Fiedler vector.
|
||||
col = np.zeros(n)
|
||||
col[1] = 1
|
||||
A = toeplitz(col)
|
||||
D = np.diag(A.sum(axis=1))
|
||||
L = D - A
|
||||
# Compute the full eigendecomposition using tricks, e.g.
|
||||
# http://www.cs.yale.edu/homes/spielman/561/2009/lect02-09.pdf
|
||||
tmp = np.pi * np.arange(n) / n
|
||||
analytic_w = 2 * (1 - np.cos(tmp))
|
||||
analytic_V = np.cos(np.outer(np.arange(n) + 1/2, tmp))
|
||||
_check_eigen(L, analytic_w, analytic_V)
|
||||
# Compute the full eigendecomposition using eigh.
|
||||
eigh_w, eigh_V = eigh(L)
|
||||
_check_eigen(L, eigh_w, eigh_V)
|
||||
# Check that the first eigenvalue is near zero and that the rest agree.
|
||||
assert_array_less(np.abs([eigh_w[0], analytic_w[0]]), 1e-14)
|
||||
assert_allclose(eigh_w[1:], analytic_w[1:])
|
||||
|
||||
# Check small lobpcg eigenvalues.
|
||||
X = analytic_V[:, :p]
|
||||
lobpcg_w, lobpcg_V = lobpcg(L, X, largest=False)
|
||||
assert_equal(lobpcg_w.shape, (p,))
|
||||
assert_equal(lobpcg_V.shape, (n, p))
|
||||
_check_eigen(L, lobpcg_w, lobpcg_V)
|
||||
assert_array_less(np.abs(np.min(lobpcg_w)), 1e-14)
|
||||
assert_allclose(np.sort(lobpcg_w)[1:], analytic_w[1:p])
|
||||
|
||||
# Check large lobpcg eigenvalues.
|
||||
X = analytic_V[:, -p:]
|
||||
lobpcg_w, lobpcg_V = lobpcg(L, X, largest=True)
|
||||
assert_equal(lobpcg_w.shape, (p,))
|
||||
assert_equal(lobpcg_V.shape, (n, p))
|
||||
_check_eigen(L, lobpcg_w, lobpcg_V)
|
||||
assert_allclose(np.sort(lobpcg_w), analytic_w[-p:])
|
||||
|
||||
# Look for the Fiedler vector using good but not exactly correct guesses.
|
||||
fiedler_guess = np.concatenate((np.ones(n//2), -np.ones(n-n//2)))
|
||||
X = np.vstack((np.ones(n), fiedler_guess)).T
|
||||
lobpcg_w, _ = lobpcg(L, X, largest=False)
|
||||
# Mathematically, the smaller eigenvalue should be zero
|
||||
# and the larger should be the algebraic connectivity.
|
||||
lobpcg_w = np.sort(lobpcg_w)
|
||||
assert_allclose(lobpcg_w, analytic_w[:2], atol=1e-14)
|
||||
|
||||
|
||||
def test_fiedler_small_8():
|
||||
"""Check the dense workaround path for small matrices.
|
||||
"""
|
||||
# This triggers the dense path because 8 < 2*5.
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
_check_fiedler(8, 2)
|
||||
|
||||
|
||||
def test_fiedler_large_12():
|
||||
"""Check the dense workaround path avoided for non-small matrices.
|
||||
"""
|
||||
# This does not trigger the dense path, because 2*5 <= 12.
|
||||
_check_fiedler(12, 2)
|
||||
|
||||
|
||||
@pytest.mark.skipif(platform.machine() == 'aarch64',
|
||||
reason="issue #15935")
|
||||
def test_failure_to_run_iterations():
|
||||
"""Check that the code exists gracefully without breaking. Issue #10974.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
X = rnd.standard_normal((100, 10))
|
||||
A = X @ X.T
|
||||
Q = rnd.standard_normal((X.shape[0], 4))
|
||||
with pytest.warns(UserWarning, match="Failed at iteration"):
|
||||
eigenvalues, _ = lobpcg(A, Q, maxiter=40, tol=1e-12)
|
||||
assert(np.max(eigenvalues) > 0)
|
||||
|
||||
|
||||
def test_failure_to_run_iterations_nonsymmetric():
|
||||
"""Check that the code exists gracefully without breaking
|
||||
if the matrix in not symmetric.
|
||||
"""
|
||||
A = np.zeros((10, 10))
|
||||
A[0, 1] = 1
|
||||
Q = np.ones((10, 1))
|
||||
with pytest.warns(UserWarning, match="Exited at iteration 2"):
|
||||
eigenvalues, _ = lobpcg(A, Q, maxiter=20)
|
||||
assert np.max(eigenvalues) > 0
|
||||
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:The problem size")
|
||||
def test_hermitian():
|
||||
"""Check complex-value Hermitian cases.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
|
||||
sizes = [3, 10, 50]
|
||||
ks = [1, 3, 10, 50]
|
||||
gens = [True, False]
|
||||
|
||||
for s, k, gen in itertools.product(sizes, ks, gens):
|
||||
if k > s:
|
||||
continue
|
||||
|
||||
H = rnd.random((s, s)) + 1.j * rnd.random((s, s))
|
||||
H = 10 * np.eye(s) + H + H.T.conj()
|
||||
|
||||
X = rnd.standard_normal((s, k))
|
||||
X = X + 1.j * rnd.standard_normal((s, k))
|
||||
|
||||
if not gen:
|
||||
B = np.eye(s)
|
||||
w, v = lobpcg(H, X, maxiter=5000)
|
||||
w0, _ = eigh(H)
|
||||
else:
|
||||
B = rnd.random((s, s)) + 1.j * rnd.random((s, s))
|
||||
B = 10 * np.eye(s) + B.dot(B.T.conj())
|
||||
w, v = lobpcg(H, X, B, maxiter=5000, largest=False)
|
||||
w0, _ = eigh(H, B)
|
||||
|
||||
for wx, vx in zip(w, v.T):
|
||||
# Check eigenvector
|
||||
assert_allclose(np.linalg.norm(H.dot(vx) - B.dot(vx) * wx)
|
||||
/ np.linalg.norm(H.dot(vx)),
|
||||
0, atol=5e-4, rtol=0)
|
||||
|
||||
# Compare eigenvalues
|
||||
j = np.argmin(abs(w0 - wx))
|
||||
assert_allclose(wx, w0[j], rtol=1e-4)
|
||||
|
||||
|
||||
# The n=5 case tests the alternative small matrix code path that uses eigh().
|
||||
@pytest.mark.filterwarnings("ignore:The problem size")
|
||||
@pytest.mark.parametrize('n, atol', [(20, 1e-3), (5, 1e-8)])
|
||||
def test_eigs_consistency(n, atol):
|
||||
"""Check eigs vs. lobpcg consistency.
|
||||
"""
|
||||
vals = np.arange(1, n+1, dtype=np.float64)
|
||||
A = spdiags(vals, 0, n, n)
|
||||
rnd = np.random.RandomState(0)
|
||||
X = rnd.random((n, 2))
|
||||
lvals, lvecs = lobpcg(A, X, largest=True, maxiter=100)
|
||||
vals, _ = eigs(A, k=2)
|
||||
|
||||
_check_eigen(A, lvals, lvecs, atol=atol, rtol=0)
|
||||
assert_allclose(np.sort(vals), np.sort(lvals), atol=1e-14)
|
||||
|
||||
|
||||
def test_verbosity():
|
||||
"""Check that nonzero verbosity level code runs.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
X = rnd.standard_normal((10, 10))
|
||||
A = X @ X.T
|
||||
Q = rnd.standard_normal((X.shape[0], 1))
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
_, _ = lobpcg(A, Q, maxiter=3, verbosityLevel=9)
|
||||
|
||||
|
||||
@pytest.mark.xfail(_IS_32BIT and sys.platform == 'win32',
|
||||
reason="tolerance violation on windows")
|
||||
@pytest.mark.xfail(platform.machine() == 'ppc64le',
|
||||
reason="fails on ppc64le")
|
||||
@pytest.mark.filterwarnings("ignore:Exited postprocessing")
|
||||
def test_tolerance_float32():
|
||||
"""Check lobpcg for attainable tolerance in float32.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
n = 50
|
||||
m = 3
|
||||
vals = -np.arange(1, n + 1)
|
||||
A = diags([vals], [0], (n, n))
|
||||
A = A.astype(np.float32)
|
||||
X = rnd.standard_normal((n, m))
|
||||
X = X.astype(np.float32)
|
||||
eigvals, _ = lobpcg(A, X, tol=1.25e-5, maxiter=50, verbosityLevel=0)
|
||||
assert_allclose(eigvals, -np.arange(1, 1 + m), atol=2e-5, rtol=1e-5)
|
||||
|
||||
|
||||
def test_random_initial_float32():
|
||||
"""Check lobpcg in float32 for specific initial.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
n = 50
|
||||
m = 4
|
||||
vals = -np.arange(1, n + 1)
|
||||
A = diags([vals], [0], (n, n))
|
||||
A = A.astype(np.float32)
|
||||
X = rnd.random((n, m))
|
||||
X = X.astype(np.float32)
|
||||
eigvals, _ = lobpcg(A, X, tol=1e-3, maxiter=50, verbosityLevel=1)
|
||||
assert_allclose(eigvals, -np.arange(1, 1 + m), atol=1e-2)
|
||||
|
||||
|
||||
def test_maxit():
|
||||
"""Check lobpcg if maxit=maxiter runs maxiter iterations and
|
||||
if maxit=None runs 20 iterations (the default)
|
||||
by checking the size of the iteration history output, which should
|
||||
be the number of iterations plus 3 (initial, final, and postprocessing)
|
||||
typically when maxiter is small and the choice of the best is passive.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
n = 50
|
||||
m = 4
|
||||
vals = -np.arange(1, n + 1)
|
||||
A = diags([vals], [0], (n, n))
|
||||
A = A.astype(np.float32)
|
||||
X = rnd.standard_normal((n, m))
|
||||
X = X.astype(np.float64)
|
||||
for maxiter in range(1, 4):
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
_, _, l_h, r_h = lobpcg(A, X, tol=1e-8, maxiter=maxiter,
|
||||
retLambdaHistory=True,
|
||||
retResidualNormsHistory=True)
|
||||
assert_allclose(np.shape(l_h)[0], maxiter+3)
|
||||
assert_allclose(np.shape(r_h)[0], maxiter+3)
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
l, _, l_h, r_h = lobpcg(A, X, tol=1e-8,
|
||||
retLambdaHistory=True,
|
||||
retResidualNormsHistory=True)
|
||||
assert_allclose(np.shape(l_h)[0], 20+3)
|
||||
assert_allclose(np.shape(r_h)[0], 20+3)
|
||||
# Check that eigenvalue output is the last one in history
|
||||
assert_allclose(l, l_h[-1])
|
||||
# Make sure that both history outputs are lists
|
||||
assert isinstance(l_h, list)
|
||||
assert isinstance(r_h, list)
|
||||
# Make sure that both history lists are arrays-like
|
||||
assert_allclose(np.shape(l_h), np.shape(np.asarray(l_h)))
|
||||
assert_allclose(np.shape(r_h), np.shape(np.asarray(r_h)))
|
||||
|
||||
|
||||
@pytest.mark.slow
|
||||
@pytest.mark.parametrize("n", [20])
|
||||
@pytest.mark.parametrize("m", [1, 3])
|
||||
@pytest.mark.filterwarnings("ignore:Exited at iteration")
|
||||
@pytest.mark.filterwarnings("ignore:Exited postprocessing")
|
||||
def test_diagonal_data_types(n, m):
|
||||
"""Check lobpcg for diagonal matrices for all matrix types.
|
||||
"""
|
||||
rnd = np.random.RandomState(0)
|
||||
# Define the generalized eigenvalue problem Av = cBv
|
||||
# where (c, v) is a generalized eigenpair,
|
||||
# and where we choose A and B to be diagonal.
|
||||
vals = np.arange(1, n + 1)
|
||||
|
||||
# list_sparse_format = ['bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil']
|
||||
list_sparse_format = ['coo']
|
||||
sparse_formats = len(list_sparse_format)
|
||||
for s_f_i, s_f in enumerate(list_sparse_format):
|
||||
|
||||
As64 = diags([vals * vals], [0], (n, n), format=s_f)
|
||||
As32 = As64.astype(np.float32)
|
||||
Af64 = As64.toarray()
|
||||
Af32 = Af64.astype(np.float32)
|
||||
|
||||
def As32f(x):
|
||||
return As32 @ x
|
||||
As32LO = LinearOperator(matvec=As32f,
|
||||
matmat=As32f,
|
||||
shape=(n, n),
|
||||
dtype=As32.dtype)
|
||||
|
||||
listA = [Af64, As64, Af32, As32, As32f, As32LO, lambda v: As32 @ v]
|
||||
|
||||
Bs64 = diags([vals], [0], (n, n), format=s_f)
|
||||
Bf64 = Bs64.toarray()
|
||||
Bs32 = Bs64.astype(np.float32)
|
||||
|
||||
def Bs32f(x):
|
||||
return Bs32 @ x
|
||||
Bs32LO = LinearOperator(matvec=Bs32f,
|
||||
matmat=Bs32f,
|
||||
shape=(n, n),
|
||||
dtype=Bs32.dtype)
|
||||
listB = [Bf64, Bs64, Bs32, Bs32f, Bs32LO, lambda v: Bs32 @ v]
|
||||
|
||||
# Define the preconditioner function as LinearOperator.
|
||||
Ms64 = diags([1./vals], [0], (n, n), format=s_f)
|
||||
|
||||
def Ms64precond(x):
|
||||
return Ms64 @ x
|
||||
Ms64precondLO = LinearOperator(matvec=Ms64precond,
|
||||
matmat=Ms64precond,
|
||||
shape=(n, n),
|
||||
dtype=Ms64.dtype)
|
||||
Mf64 = Ms64.toarray()
|
||||
|
||||
def Mf64precond(x):
|
||||
return Mf64 @ x
|
||||
Mf64precondLO = LinearOperator(matvec=Mf64precond,
|
||||
matmat=Mf64precond,
|
||||
shape=(n, n),
|
||||
dtype=Mf64.dtype)
|
||||
Ms32 = Ms64.astype(np.float32)
|
||||
|
||||
def Ms32precond(x):
|
||||
return Ms32 @ x
|
||||
Ms32precondLO = LinearOperator(matvec=Ms32precond,
|
||||
matmat=Ms32precond,
|
||||
shape=(n, n),
|
||||
dtype=Ms32.dtype)
|
||||
Mf32 = Ms32.toarray()
|
||||
|
||||
def Mf32precond(x):
|
||||
return Mf32 @ x
|
||||
Mf32precondLO = LinearOperator(matvec=Mf32precond,
|
||||
matmat=Mf32precond,
|
||||
shape=(n, n),
|
||||
dtype=Mf32.dtype)
|
||||
listM = [None, Ms64, Ms64precondLO, Mf64precondLO, Ms64precond,
|
||||
Ms32, Ms32precondLO, Mf32precondLO, Ms32precond]
|
||||
|
||||
# Setup matrix of the initial approximation to the eigenvectors
|
||||
# (cannot be sparse array).
|
||||
Xf64 = rnd.random((n, m))
|
||||
Xf32 = Xf64.astype(np.float32)
|
||||
listX = [Xf64, Xf32]
|
||||
|
||||
# Require that the returned eigenvectors be in the orthogonal complement
|
||||
# of the first few standard basis vectors (cannot be sparse array).
|
||||
m_excluded = 3
|
||||
Yf64 = np.eye(n, m_excluded, dtype=float)
|
||||
Yf32 = np.eye(n, m_excluded, dtype=np.float32)
|
||||
listY = [Yf64, Yf32]
|
||||
|
||||
tests = list(itertools.product(listA, listB, listM, listX, listY))
|
||||
# This is one of the slower tests because there are >1,000 configs
|
||||
# to test here, instead of checking product of all input, output types
|
||||
# test each configuration for the first sparse format, and then
|
||||
# for one additional sparse format. this takes 2/7=30% as long as
|
||||
# testing all configurations for all sparse formats.
|
||||
if s_f_i > 0:
|
||||
tests = tests[s_f_i - 1::sparse_formats-1]
|
||||
|
||||
for A, B, M, X, Y in tests:
|
||||
eigvals, _ = lobpcg(A, X, B=B, M=M, Y=Y, tol=1e-4,
|
||||
maxiter=100, largest=False)
|
||||
assert_allclose(eigvals,
|
||||
np.arange(1 + m_excluded, 1 + m_excluded + m),
|
||||
atol=1e-5)
|
||||
Binary file not shown.
Binary file not shown.
@@ -1,907 +0,0 @@
|
||||
import os
|
||||
import re
|
||||
import copy
|
||||
import numpy as np
|
||||
|
||||
from numpy.testing import assert_allclose, assert_equal, assert_array_equal
|
||||
import pytest
|
||||
|
||||
from scipy.linalg import svd, null_space
|
||||
from scipy.sparse import csc_matrix, isspmatrix, spdiags, random
|
||||
from scipy.sparse.linalg import LinearOperator, aslinearoperator
|
||||
if os.environ.get("SCIPY_USE_PROPACK"):
|
||||
has_propack = True
|
||||
else:
|
||||
has_propack = False
|
||||
from scipy.sparse.linalg import svds
|
||||
from scipy.sparse.linalg._eigen.arpack import ArpackNoConvergence
|
||||
|
||||
|
||||
# --- Helper Functions / Classes ---
|
||||
|
||||
|
||||
def sorted_svd(m, k, which='LM'):
|
||||
# Compute svd of a dense matrix m, and return singular vectors/values
|
||||
# sorted.
|
||||
if isspmatrix(m):
|
||||
m = m.toarray()
|
||||
u, s, vh = svd(m)
|
||||
if which == 'LM':
|
||||
ii = np.argsort(s)[-k:]
|
||||
elif which == 'SM':
|
||||
ii = np.argsort(s)[:k]
|
||||
else:
|
||||
raise ValueError("unknown which=%r" % (which,))
|
||||
|
||||
return u[:, ii], s[ii], vh[ii]
|
||||
|
||||
|
||||
def svd_estimate(u, s, vh):
|
||||
return np.dot(u, np.dot(np.diag(s), vh))
|
||||
|
||||
|
||||
def _check_svds(A, k, u, s, vh, which="LM", check_usvh_A=False,
|
||||
check_svd=True, atol=1e-10, rtol=1e-7):
|
||||
n, m = A.shape
|
||||
|
||||
# Check shapes.
|
||||
assert_equal(u.shape, (n, k))
|
||||
assert_equal(s.shape, (k,))
|
||||
assert_equal(vh.shape, (k, m))
|
||||
|
||||
# Check that the original matrix can be reconstituted.
|
||||
A_rebuilt = (u*s).dot(vh)
|
||||
assert_equal(A_rebuilt.shape, A.shape)
|
||||
if check_usvh_A:
|
||||
assert_allclose(A_rebuilt, A, atol=atol, rtol=rtol)
|
||||
|
||||
# Check that u is a semi-orthogonal matrix.
|
||||
uh_u = np.dot(u.T.conj(), u)
|
||||
assert_equal(uh_u.shape, (k, k))
|
||||
assert_allclose(uh_u, np.identity(k), atol=atol, rtol=rtol)
|
||||
|
||||
# Check that vh is a semi-orthogonal matrix.
|
||||
vh_v = np.dot(vh, vh.T.conj())
|
||||
assert_equal(vh_v.shape, (k, k))
|
||||
assert_allclose(vh_v, np.identity(k), atol=atol, rtol=rtol)
|
||||
|
||||
# Check that scipy.sparse.linalg.svds ~ scipy.linalg.svd
|
||||
if check_svd:
|
||||
u2, s2, vh2 = sorted_svd(A, k, which)
|
||||
assert_allclose(np.abs(u), np.abs(u2), atol=atol, rtol=rtol)
|
||||
assert_allclose(s, s2, atol=atol, rtol=rtol)
|
||||
assert_allclose(np.abs(vh), np.abs(vh2), atol=atol, rtol=rtol)
|
||||
|
||||
|
||||
def _check_svds_n(A, k, u, s, vh, which="LM", check_res=True,
|
||||
check_svd=True, atol=1e-10, rtol=1e-7):
|
||||
n, m = A.shape
|
||||
|
||||
# Check shapes.
|
||||
assert_equal(u.shape, (n, k))
|
||||
assert_equal(s.shape, (k,))
|
||||
assert_equal(vh.shape, (k, m))
|
||||
|
||||
# Check that u is a semi-orthogonal matrix.
|
||||
uh_u = np.dot(u.T.conj(), u)
|
||||
assert_equal(uh_u.shape, (k, k))
|
||||
error = np.sum(np.abs(uh_u - np.identity(k))) / (k * k)
|
||||
assert_allclose(error, 0.0, atol=atol, rtol=rtol)
|
||||
|
||||
# Check that vh is a semi-orthogonal matrix.
|
||||
vh_v = np.dot(vh, vh.T.conj())
|
||||
assert_equal(vh_v.shape, (k, k))
|
||||
error = np.sum(np.abs(vh_v - np.identity(k))) / (k * k)
|
||||
assert_allclose(error, 0.0, atol=atol, rtol=rtol)
|
||||
|
||||
# Check residuals
|
||||
if check_res:
|
||||
ru = A.T.conj() @ u - vh.T.conj() * s
|
||||
rus = np.sum(np.abs(ru)) / (n * k)
|
||||
rvh = A @ vh.T.conj() - u * s
|
||||
rvhs = np.sum(np.abs(rvh)) / (m * k)
|
||||
assert_allclose(rus, 0.0, atol=atol, rtol=rtol)
|
||||
assert_allclose(rvhs, 0.0, atol=atol, rtol=rtol)
|
||||
|
||||
# Check that scipy.sparse.linalg.svds ~ scipy.linalg.svd
|
||||
if check_svd:
|
||||
u2, s2, vh2 = sorted_svd(A, k, which)
|
||||
assert_allclose(s, s2, atol=atol, rtol=rtol)
|
||||
A_rebuilt_svd = (u2*s2).dot(vh2)
|
||||
A_rebuilt = (u*s).dot(vh)
|
||||
assert_equal(A_rebuilt.shape, A.shape)
|
||||
error = np.sum(np.abs(A_rebuilt_svd - A_rebuilt)) / (k * k)
|
||||
assert_allclose(error, 0.0, atol=atol, rtol=rtol)
|
||||
|
||||
|
||||
class CheckingLinearOperator(LinearOperator):
|
||||
def __init__(self, A):
|
||||
self.A = A
|
||||
self.dtype = A.dtype
|
||||
self.shape = A.shape
|
||||
|
||||
def _matvec(self, x):
|
||||
assert_equal(max(x.shape), np.size(x))
|
||||
return self.A.dot(x)
|
||||
|
||||
def _rmatvec(self, x):
|
||||
assert_equal(max(x.shape), np.size(x))
|
||||
return self.A.T.conjugate().dot(x)
|
||||
|
||||
|
||||
# --- Test Input Validation ---
|
||||
# Tests input validation on parameters `k` and `which`.
|
||||
# Needs better input validation checks for all other parameters.
|
||||
|
||||
class SVDSCommonTests:
|
||||
|
||||
solver = None
|
||||
|
||||
# some of these IV tests could run only once, say with solver=None
|
||||
|
||||
_A_empty_msg = "`A` must not be empty."
|
||||
_A_dtype_msg = "`A` must be of floating or complex floating data type"
|
||||
_A_type_msg = "type not understood"
|
||||
_A_ndim_msg = "array must have ndim <= 2"
|
||||
_A_validation_inputs = [
|
||||
(np.asarray([[]]), ValueError, _A_empty_msg),
|
||||
(np.asarray([[1, 2], [3, 4]]), ValueError, _A_dtype_msg),
|
||||
("hi", TypeError, _A_type_msg),
|
||||
(np.asarray([[[1., 2.], [3., 4.]]]), ValueError, _A_ndim_msg)]
|
||||
|
||||
@pytest.mark.parametrize("args", _A_validation_inputs)
|
||||
def test_svds_input_validation_A(self, args):
|
||||
A, error_type, message = args
|
||||
with pytest.raises(error_type, match=message):
|
||||
svds(A, k=1, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("k", [-1, 0, 3, 4, 5, 1.5, "1"])
|
||||
def test_svds_input_validation_k_1(self, k):
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((4, 3))
|
||||
|
||||
# propack can do complete SVD
|
||||
if self.solver == 'propack' and k == 3:
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not enabled")
|
||||
res = svds(A, k=k, solver=self.solver)
|
||||
_check_svds(A, k, *res, check_usvh_A=True, check_svd=True)
|
||||
return
|
||||
|
||||
message = ("`k` must be an integer satisfying")
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(A, k=k, solver=self.solver)
|
||||
|
||||
def test_svds_input_validation_k_2(self):
|
||||
# I think the stack trace is reasonable when `k` can't be converted
|
||||
# to an int.
|
||||
message = "int() argument must be a"
|
||||
with pytest.raises(TypeError, match=re.escape(message)):
|
||||
svds(np.eye(10), k=[], solver=self.solver)
|
||||
|
||||
message = "invalid literal for int()"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), k="hi", solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("tol", (-1, np.inf, np.nan))
|
||||
def test_svds_input_validation_tol_1(self, tol):
|
||||
message = "`tol` must be a non-negative floating point value."
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), tol=tol, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("tol", ([], 'hi'))
|
||||
def test_svds_input_validation_tol_2(self, tol):
|
||||
# I think the stack trace is reasonable here
|
||||
message = "'<' not supported between instances"
|
||||
with pytest.raises(TypeError, match=message):
|
||||
svds(np.eye(10), tol=tol, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("which", ('LA', 'SA', 'ekki', 0))
|
||||
def test_svds_input_validation_which(self, which):
|
||||
# Regression test for a github issue.
|
||||
# https://github.com/scipy/scipy/issues/4590
|
||||
# Function was not checking for eigenvalue type and unintended
|
||||
# values could be returned.
|
||||
with pytest.raises(ValueError, match="`which` must be in"):
|
||||
svds(np.eye(10), which=which, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("transpose", (True, False))
|
||||
@pytest.mark.parametrize("n", range(4, 9))
|
||||
def test_svds_input_validation_v0_1(self, transpose, n):
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((5, 7))
|
||||
v0 = rng.random(n)
|
||||
if transpose:
|
||||
A = A.T
|
||||
k = 2
|
||||
message = "`v0` must have shape"
|
||||
|
||||
required_length = (A.shape[0] if self.solver == 'propack'
|
||||
else min(A.shape))
|
||||
if n != required_length:
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(A, k=k, v0=v0, solver=self.solver)
|
||||
|
||||
def test_svds_input_validation_v0_2(self):
|
||||
A = np.ones((10, 10))
|
||||
v0 = np.ones((1, 10))
|
||||
message = "`v0` must have shape"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(A, k=1, v0=v0, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("v0", ("hi", 1, np.ones(10, dtype=int)))
|
||||
def test_svds_input_validation_v0_3(self, v0):
|
||||
A = np.ones((10, 10))
|
||||
message = "`v0` must be of floating or complex floating data type."
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(A, k=1, v0=v0, solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("maxiter", (-1, 0, 5.5))
|
||||
def test_svds_input_validation_maxiter_1(self, maxiter):
|
||||
message = ("`maxiter` must be a positive integer.")
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), maxiter=maxiter, solver=self.solver)
|
||||
|
||||
def test_svds_input_validation_maxiter_2(self):
|
||||
# I think the stack trace is reasonable when `k` can't be converted
|
||||
# to an int.
|
||||
message = "int() argument must be a"
|
||||
with pytest.raises(TypeError, match=re.escape(message)):
|
||||
svds(np.eye(10), maxiter=[], solver=self.solver)
|
||||
|
||||
message = "invalid literal for int()"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), maxiter="hi", solver=self.solver)
|
||||
|
||||
@pytest.mark.parametrize("rsv", ('ekki', 10))
|
||||
def test_svds_input_validation_return_singular_vectors(self, rsv):
|
||||
message = "`return_singular_vectors` must be in"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), return_singular_vectors=rsv, solver=self.solver)
|
||||
|
||||
# --- Test Parameters ---
|
||||
|
||||
@pytest.mark.parametrize("k", [3, 5])
|
||||
@pytest.mark.parametrize("which", ["LM", "SM"])
|
||||
def test_svds_parameter_k_which(self, k, which):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
# check that the `k` parameter sets the number of eigenvalues/
|
||||
# eigenvectors returned.
|
||||
# Also check that the `which` parameter sets whether the largest or
|
||||
# smallest eigenvalues are returned
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((10, 10))
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
res = svds(A, k=k, which=which, solver=self.solver,
|
||||
random_state=0)
|
||||
else:
|
||||
res = svds(A, k=k, which=which, solver=self.solver,
|
||||
random_state=0)
|
||||
_check_svds(A, k, *res, which=which, atol=8e-10)
|
||||
|
||||
# loop instead of parametrize for simplicity
|
||||
def test_svds_parameter_tol(self):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
return # TODO: needs work, disabling for now
|
||||
# check the effect of the `tol` parameter on solver accuracy by solving
|
||||
# the same problem with varying `tol` and comparing the eigenvalues
|
||||
# against ground truth computed
|
||||
n = 100 # matrix size
|
||||
k = 3 # number of eigenvalues to check
|
||||
|
||||
# generate a random, sparse-ish matrix
|
||||
# effect isn't apparent for matrices that are too small
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((n, n))
|
||||
A[A > .1] = 0
|
||||
A = A @ A.T
|
||||
|
||||
_, s, _ = svd(A) # calculate ground truth
|
||||
|
||||
# calculate the error as a function of `tol`
|
||||
A = csc_matrix(A)
|
||||
|
||||
def err(tol):
|
||||
if self.solver == 'lobpcg' and tol == 1e-4:
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
_, s2, _ = svds(A, k=k, v0=np.ones(n),
|
||||
solver=self.solver, tol=tol)
|
||||
else:
|
||||
_, s2, _ = svds(A, k=k, v0=np.ones(n),
|
||||
solver=self.solver, tol=tol)
|
||||
return np.linalg.norm((s2 - s[k-1::-1])/s[k-1::-1])
|
||||
|
||||
tols = [1e-4, 1e-2, 1e0] # tolerance levels to check
|
||||
# for 'arpack' and 'propack', accuracies make discrete steps
|
||||
accuracies = {'propack': [1e-12, 1e-6, 1e-4],
|
||||
'arpack': [2e-15, 1e-10, 1e-10],
|
||||
'lobpcg': [1e-11, 1e-3, 10]}
|
||||
|
||||
for tol, accuracy in zip(tols, accuracies[self.solver]):
|
||||
error = err(tol)
|
||||
assert error < accuracy
|
||||
assert error > accuracy/10
|
||||
|
||||
def test_svd_v0(self):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
# check that the `v0` parameter affects the solution
|
||||
n = 100
|
||||
k = 1
|
||||
# If k != 1, LOBPCG needs more initial vectors, which are generated
|
||||
# with random_state, so it does not pass w/ k >= 2.
|
||||
# For some other values of `n`, the AssertionErrors are not raised
|
||||
# with different v0s, which is reasonable.
|
||||
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((n, n))
|
||||
|
||||
# with the same v0, solutions are the same, and they are accurate
|
||||
# v0 takes precedence over random_state
|
||||
v0a = rng.random(n)
|
||||
res1a = svds(A, k, v0=v0a, solver=self.solver, random_state=0)
|
||||
res2a = svds(A, k, v0=v0a, solver=self.solver, random_state=1)
|
||||
assert_equal(res1a, res2a)
|
||||
_check_svds(A, k, *res1a)
|
||||
|
||||
# with the same v0, solutions are the same, and they are accurate
|
||||
v0b = rng.random(n)
|
||||
res1b = svds(A, k, v0=v0b, solver=self.solver, random_state=2)
|
||||
res2b = svds(A, k, v0=v0b, solver=self.solver, random_state=3)
|
||||
assert_equal(res1b, res2b)
|
||||
_check_svds(A, k, *res1b)
|
||||
|
||||
# with different v0, solutions can be numerically different
|
||||
message = "Arrays are not equal"
|
||||
with pytest.raises(AssertionError, match=message):
|
||||
assert_equal(res1a, res1b)
|
||||
|
||||
def test_svd_random_state(self):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
# check that the `random_state` parameter affects the solution
|
||||
# Admittedly, `n` and `k` are chosen so that all solver pass all
|
||||
# these checks. That's a tall order, since LOBPCG doesn't want to
|
||||
# achieve the desired accuracy and ARPACK often returns the same
|
||||
# singular values/vectors for different v0.
|
||||
n = 100
|
||||
k = 1
|
||||
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((n, n))
|
||||
|
||||
# with the same random_state, solutions are the same and accurate
|
||||
res1a = svds(A, k, solver=self.solver, random_state=0)
|
||||
res2a = svds(A, k, solver=self.solver, random_state=0)
|
||||
assert_equal(res1a, res2a)
|
||||
_check_svds(A, k, *res1a)
|
||||
|
||||
# with the same random_state, solutions are the same and accurate
|
||||
res1b = svds(A, k, solver=self.solver, random_state=1)
|
||||
res2b = svds(A, k, solver=self.solver, random_state=1)
|
||||
assert_equal(res1b, res2b)
|
||||
_check_svds(A, k, *res1b)
|
||||
|
||||
# with different random_state, solutions can be numerically different
|
||||
message = "Arrays are not equal"
|
||||
with pytest.raises(AssertionError, match=message):
|
||||
assert_equal(res1a, res1b)
|
||||
|
||||
@pytest.mark.parametrize("random_state", (0, 1,
|
||||
np.random.RandomState(0),
|
||||
np.random.default_rng(0)))
|
||||
def test_svd_random_state_2(self, random_state):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
|
||||
n = 100
|
||||
k = 1
|
||||
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((n, n))
|
||||
|
||||
random_state_2 = copy.deepcopy(random_state)
|
||||
|
||||
# with the same random_state, solutions are the same and accurate
|
||||
res1a = svds(A, k, solver=self.solver, random_state=random_state)
|
||||
res2a = svds(A, k, solver=self.solver, random_state=random_state_2)
|
||||
assert_equal(res1a, res2a)
|
||||
_check_svds(A, k, *res1a)
|
||||
|
||||
@pytest.mark.parametrize("random_state", (None,
|
||||
np.random.RandomState(0),
|
||||
np.random.default_rng(0)))
|
||||
def test_svd_random_state_3(self, random_state):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
|
||||
n = 100
|
||||
k = 5
|
||||
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((n, n))
|
||||
|
||||
# random_state in different state produces accurate - but not
|
||||
# not necessarily identical - results
|
||||
res1a = svds(A, k, solver=self.solver, random_state=random_state)
|
||||
res2a = svds(A, k, solver=self.solver, random_state=random_state)
|
||||
_check_svds(A, k, *res1a, atol=2e-10, rtol=1e-6)
|
||||
_check_svds(A, k, *res2a, atol=2e-10, rtol=1e-6)
|
||||
|
||||
message = "Arrays are not equal"
|
||||
with pytest.raises(AssertionError, match=message):
|
||||
assert_equal(res1a, res2a)
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:Exited postprocessing")
|
||||
def test_svd_maxiter(self):
|
||||
# check that maxiter works as expected: should not return accurate
|
||||
# solution after 1 iteration, but should with default `maxiter`
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
A = np.diag(np.arange(9)).astype(np.float64)
|
||||
k = 1
|
||||
u, s, vh = sorted_svd(A, k)
|
||||
|
||||
if self.solver == 'arpack':
|
||||
message = "ARPACK error -1: No convergence"
|
||||
with pytest.raises(ArpackNoConvergence, match=message):
|
||||
svds(A, k, ncv=3, maxiter=1, solver=self.solver)
|
||||
elif self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="Exited at iteration"):
|
||||
svds(A, k, maxiter=1, solver=self.solver)
|
||||
elif self.solver == 'propack':
|
||||
message = "k=1 singular triplets did not converge within"
|
||||
with pytest.raises(np.linalg.LinAlgError, match=message):
|
||||
svds(A, k, maxiter=1, solver=self.solver)
|
||||
|
||||
ud, sd, vhd = svds(A, k, solver=self.solver) # default maxiter
|
||||
_check_svds(A, k, ud, sd, vhd, atol=1e-8)
|
||||
assert_allclose(np.abs(ud), np.abs(u), atol=1e-8)
|
||||
assert_allclose(np.abs(vhd), np.abs(vh), atol=1e-8)
|
||||
assert_allclose(np.abs(sd), np.abs(s), atol=1e-9)
|
||||
|
||||
@pytest.mark.parametrize("rsv", (True, False, 'u', 'vh'))
|
||||
@pytest.mark.parametrize("shape", ((5, 7), (6, 6), (7, 5)))
|
||||
def test_svd_return_singular_vectors(self, rsv, shape):
|
||||
# check that the return_singular_vectors parameter works as expected
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random(shape)
|
||||
k = 2
|
||||
M, N = shape
|
||||
u, s, vh = sorted_svd(A, k)
|
||||
|
||||
respect_u = True if self.solver == 'propack' else M <= N
|
||||
respect_vh = True if self.solver == 'propack' else M > N
|
||||
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
if rsv is False:
|
||||
s2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert_allclose(s2, s)
|
||||
elif rsv == 'u' and respect_u:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert_allclose(np.abs(u2), np.abs(u))
|
||||
assert_allclose(s2, s)
|
||||
assert vh2 is None
|
||||
elif rsv == 'vh' and respect_vh:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert u2 is None
|
||||
assert_allclose(s2, s)
|
||||
assert_allclose(np.abs(vh2), np.abs(vh))
|
||||
else:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
if u2 is not None:
|
||||
assert_allclose(np.abs(u2), np.abs(u))
|
||||
assert_allclose(s2, s)
|
||||
if vh2 is not None:
|
||||
assert_allclose(np.abs(vh2), np.abs(vh))
|
||||
else:
|
||||
if rsv is False:
|
||||
s2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert_allclose(s2, s)
|
||||
elif rsv == 'u' and respect_u:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert_allclose(np.abs(u2), np.abs(u))
|
||||
assert_allclose(s2, s)
|
||||
assert vh2 is None
|
||||
elif rsv == 'vh' and respect_vh:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
assert u2 is None
|
||||
assert_allclose(s2, s)
|
||||
assert_allclose(np.abs(vh2), np.abs(vh))
|
||||
else:
|
||||
u2, s2, vh2 = svds(A, k, return_singular_vectors=rsv,
|
||||
solver=self.solver, random_state=rng)
|
||||
if u2 is not None:
|
||||
assert_allclose(np.abs(u2), np.abs(u))
|
||||
assert_allclose(s2, s)
|
||||
if vh2 is not None:
|
||||
assert_allclose(np.abs(vh2), np.abs(vh))
|
||||
|
||||
# --- Test Basic Functionality ---
|
||||
# Tests the accuracy of each solver for real and complex matrices provided
|
||||
# as list, dense array, sparse matrix, and LinearOperator.
|
||||
|
||||
A1 = [[1, 2, 3], [3, 4, 3], [1 + 1j, 0, 2], [0, 0, 1]]
|
||||
A2 = [[1, 2, 3, 8 + 5j], [3 - 2j, 4, 3, 5], [1, 0, 2, 3], [0, 0, 1, 0]]
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:k >= N - 1",
|
||||
reason="needed to demonstrate #16725")
|
||||
@pytest.mark.parametrize('A', (A1, A2))
|
||||
@pytest.mark.parametrize('k', range(1, 5))
|
||||
# PROPACK fails a lot if @pytest.mark.parametrize('which', ("SM", "LM"))
|
||||
@pytest.mark.parametrize('real', (True, False))
|
||||
@pytest.mark.parametrize('transpose', (False, True))
|
||||
# In gh-14299, it was suggested the `svds` should _not_ work with lists
|
||||
@pytest.mark.parametrize('lo_type', (np.asarray, csc_matrix,
|
||||
aslinearoperator))
|
||||
def test_svd_simple(self, A, k, real, transpose, lo_type):
|
||||
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
|
||||
A = np.asarray(A)
|
||||
A = np.real(A) if real else A
|
||||
A = A.T if transpose else A
|
||||
A2 = lo_type(A)
|
||||
|
||||
# could check for the appropriate errors, but that is tested above
|
||||
if k > min(A.shape):
|
||||
pytest.skip("`k` cannot be greater than `min(A.shape)`")
|
||||
if self.solver != 'propack' and k >= min(A.shape):
|
||||
pytest.skip("Only PROPACK supports complete SVD")
|
||||
if self.solver == 'arpack' and not real and k == min(A.shape) - 1:
|
||||
pytest.skip("#16725")
|
||||
|
||||
if self.solver == 'propack' and (np.intp(0).itemsize < 8 and not real):
|
||||
pytest.skip('PROPACK complex-valued SVD methods not available '
|
||||
'for 32-bit builds')
|
||||
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
u, s, vh = svds(A2, k, solver=self.solver)
|
||||
else:
|
||||
u, s, vh = svds(A2, k, solver=self.solver)
|
||||
_check_svds(A, k, u, s, vh, atol=3e-10)
|
||||
|
||||
def test_svd_linop(self):
|
||||
solver = self.solver
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not available")
|
||||
|
||||
nmks = [(6, 7, 3),
|
||||
(9, 5, 4),
|
||||
(10, 8, 5)]
|
||||
|
||||
def reorder(args):
|
||||
U, s, VH = args
|
||||
j = np.argsort(s)
|
||||
return U[:, j], s[j], VH[j, :]
|
||||
|
||||
for n, m, k in nmks:
|
||||
# Test svds on a LinearOperator.
|
||||
A = np.random.RandomState(52).randn(n, m)
|
||||
L = CheckingLinearOperator(A)
|
||||
|
||||
if solver == 'propack':
|
||||
v0 = np.ones(n)
|
||||
else:
|
||||
v0 = np.ones(min(A.shape))
|
||||
if solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
U1, s1, VH1 = reorder(svds(A, k, v0=v0, solver=solver))
|
||||
U2, s2, VH2 = reorder(svds(L, k, v0=v0, solver=solver))
|
||||
else:
|
||||
U1, s1, VH1 = reorder(svds(A, k, v0=v0, solver=solver))
|
||||
U2, s2, VH2 = reorder(svds(L, k, v0=v0, solver=solver))
|
||||
|
||||
assert_allclose(np.abs(U1), np.abs(U2))
|
||||
assert_allclose(s1, s2)
|
||||
assert_allclose(np.abs(VH1), np.abs(VH2))
|
||||
assert_allclose(np.dot(U1, np.dot(np.diag(s1), VH1)),
|
||||
np.dot(U2, np.dot(np.diag(s2), VH2)))
|
||||
|
||||
# Try again with which="SM".
|
||||
A = np.random.RandomState(1909).randn(n, m)
|
||||
L = CheckingLinearOperator(A)
|
||||
|
||||
# TODO: arpack crashes when v0=v0, which="SM"
|
||||
kwargs = {'v0': v0} if solver not in {None, 'arpack'} else {}
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
U1, s1, VH1 = reorder(svds(A, k, which="SM", solver=solver,
|
||||
**kwargs))
|
||||
U2, s2, VH2 = reorder(svds(L, k, which="SM", solver=solver,
|
||||
**kwargs))
|
||||
else:
|
||||
U1, s1, VH1 = reorder(svds(A, k, which="SM", solver=solver,
|
||||
**kwargs))
|
||||
U2, s2, VH2 = reorder(svds(L, k, which="SM", solver=solver,
|
||||
**kwargs))
|
||||
|
||||
assert_allclose(np.abs(U1), np.abs(U2))
|
||||
assert_allclose(s1 + 1, s2 + 1)
|
||||
assert_allclose(np.abs(VH1), np.abs(VH2))
|
||||
assert_allclose(np.dot(U1, np.dot(np.diag(s1), VH1)),
|
||||
np.dot(U2, np.dot(np.diag(s2), VH2)))
|
||||
|
||||
if k < min(n, m) - 1:
|
||||
# Complex input and explicit which="LM".
|
||||
for (dt, eps) in [(complex, 1e-7), (np.complex64, 1e-3)]:
|
||||
if self.solver == 'propack' and np.intp(0).itemsize < 8:
|
||||
pytest.skip('PROPACK complex-valued SVD methods '
|
||||
'not available for 32-bit builds')
|
||||
rng = np.random.RandomState(1648)
|
||||
A = (rng.randn(n, m) + 1j * rng.randn(n, m)).astype(dt)
|
||||
L = CheckingLinearOperator(A)
|
||||
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning,
|
||||
match="The problem size"):
|
||||
U1, s1, VH1 = reorder(svds(A, k, which="LM",
|
||||
solver=solver))
|
||||
U2, s2, VH2 = reorder(svds(L, k, which="LM",
|
||||
solver=solver))
|
||||
else:
|
||||
U1, s1, VH1 = reorder(svds(A, k, which="LM",
|
||||
solver=solver))
|
||||
U2, s2, VH2 = reorder(svds(L, k, which="LM",
|
||||
solver=solver))
|
||||
|
||||
assert_allclose(np.abs(U1), np.abs(U2), rtol=eps)
|
||||
assert_allclose(s1, s2, rtol=eps)
|
||||
assert_allclose(np.abs(VH1), np.abs(VH2), rtol=eps)
|
||||
assert_allclose(np.dot(U1, np.dot(np.diag(s1), VH1)),
|
||||
np.dot(U2, np.dot(np.diag(s2), VH2)),
|
||||
rtol=eps)
|
||||
|
||||
SHAPES = ((100, 100), (100, 101), (101, 100))
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:Exited at iteration")
|
||||
@pytest.mark.filterwarnings("ignore:Exited postprocessing")
|
||||
@pytest.mark.parametrize("shape", SHAPES)
|
||||
# ARPACK supports only dtype float, complex, or np.float32
|
||||
@pytest.mark.parametrize("dtype", (float, complex, np.float32))
|
||||
def test_small_sigma_sparse(self, shape, dtype):
|
||||
# https://github.com/scipy/scipy/pull/11829
|
||||
solver = self.solver
|
||||
# 2do: PROPACK fails orthogonality of singular vectors
|
||||
# if dtype == complex and self.solver == 'propack':
|
||||
# pytest.skip("PROPACK unsupported for complex dtype")
|
||||
if solver == 'propack':
|
||||
pytest.skip("PROPACK failures unrelated to PR")
|
||||
rng = np.random.default_rng(0)
|
||||
k = 5
|
||||
(m, n) = shape
|
||||
S = random(m, n, density=0.1, random_state=rng)
|
||||
if dtype == complex:
|
||||
S = + 1j * random(m, n, density=0.1, random_state=rng)
|
||||
e = np.ones(m)
|
||||
e[0:5] *= 1e1 ** np.arange(-5, 0, 1)
|
||||
S = spdiags(e, 0, m, m) @ S
|
||||
S = S.astype(dtype)
|
||||
u, s, vh = svds(S, k, which='SM', solver=solver, maxiter=1000)
|
||||
c_svd = False # partial SVD can be different from full SVD
|
||||
_check_svds_n(S, k, u, s, vh, which="SM", check_svd=c_svd, atol=1e-1)
|
||||
|
||||
# --- Test Edge Cases ---
|
||||
# Checks a few edge cases.
|
||||
|
||||
@pytest.mark.parametrize("shape", ((6, 5), (5, 5), (5, 6)))
|
||||
@pytest.mark.parametrize("dtype", (float, complex))
|
||||
def test_svd_LM_ones_matrix(self, shape, dtype):
|
||||
# Check that svds can deal with matrix_rank less than k in LM mode.
|
||||
k = 3
|
||||
n, m = shape
|
||||
A = np.ones((n, m), dtype=dtype)
|
||||
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
U, s, VH = svds(A, k, solver=self.solver)
|
||||
else:
|
||||
U, s, VH = svds(A, k, solver=self.solver)
|
||||
|
||||
_check_svds(A, k, U, s, VH, check_usvh_A=True, check_svd=False)
|
||||
|
||||
# Check that the largest singular value is near sqrt(n*m)
|
||||
# and the other singular values have been forced to zero.
|
||||
assert_allclose(np.max(s), np.sqrt(n*m))
|
||||
s = np.array(sorted(s)[:-1]) + 1
|
||||
z = np.ones_like(s)
|
||||
assert_allclose(s, z)
|
||||
|
||||
@pytest.mark.filterwarnings("ignore:k >= N - 1",
|
||||
reason="needed to demonstrate #16725")
|
||||
@pytest.mark.parametrize("shape", ((3, 4), (4, 4), (4, 3), (4, 2)))
|
||||
@pytest.mark.parametrize("dtype", (float, complex))
|
||||
def test_zero_matrix(self, shape, dtype):
|
||||
# Check that svds can deal with matrices containing only zeros;
|
||||
# see https://github.com/scipy/scipy/issues/3452/
|
||||
# shape = (4, 2) is included because it is the particular case
|
||||
# reported in the issue
|
||||
k = 1
|
||||
n, m = shape
|
||||
A = np.zeros((n, m), dtype=dtype)
|
||||
|
||||
if (self.solver == 'arpack' and dtype is complex
|
||||
and k == min(A.shape) - 1):
|
||||
pytest.skip("#16725")
|
||||
|
||||
if self.solver == 'propack':
|
||||
pytest.skip("PROPACK failures unrelated to PR #16712")
|
||||
|
||||
if self.solver == 'lobpcg':
|
||||
with pytest.warns(UserWarning, match="The problem size"):
|
||||
U, s, VH = svds(A, k, solver=self.solver)
|
||||
else:
|
||||
U, s, VH = svds(A, k, solver=self.solver)
|
||||
|
||||
# Check some generic properties of svd.
|
||||
_check_svds(A, k, U, s, VH, check_usvh_A=True, check_svd=False)
|
||||
|
||||
# Check that the singular values are zero.
|
||||
assert_array_equal(s, 0)
|
||||
|
||||
@pytest.mark.parametrize("shape", ((20, 20), (20, 21), (21, 20)))
|
||||
# ARPACK supports only dtype float, complex, or np.float32
|
||||
@pytest.mark.parametrize("dtype", (float, complex, np.float32))
|
||||
def test_small_sigma(self, shape, dtype):
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not enabled")
|
||||
# https://github.com/scipy/scipy/pull/11829
|
||||
if dtype == complex and self.solver == 'propack':
|
||||
pytest.skip("PROPACK unsupported for complex dtype")
|
||||
rng = np.random.default_rng(179847540)
|
||||
A = rng.random(shape).astype(dtype)
|
||||
u, _, vh = svd(A, full_matrices=False)
|
||||
if dtype == np.float32:
|
||||
e = 10.0
|
||||
else:
|
||||
e = 100.0
|
||||
t = e**(-np.arange(len(vh))).astype(dtype)
|
||||
A = (u*t).dot(vh)
|
||||
k = 4
|
||||
u, s, vh = svds(A, k, solver=self.solver, maxiter=100)
|
||||
t = np.sum(s > 0)
|
||||
assert_equal(t, k)
|
||||
# LOBPCG needs larger atol and rtol to pass
|
||||
_check_svds_n(A, k, u, s, vh, atol=1e-3, rtol=1e0, check_svd=False)
|
||||
|
||||
# ARPACK supports only dtype float, complex, or np.float32
|
||||
@pytest.mark.filterwarnings("ignore:The problem size")
|
||||
@pytest.mark.parametrize("dtype", (float, complex, np.float32))
|
||||
def test_small_sigma2(self, dtype):
|
||||
if self.solver == 'propack':
|
||||
if not has_propack:
|
||||
pytest.skip("PROPACK not enabled")
|
||||
elif dtype == np.float32:
|
||||
pytest.skip("Test failures in CI, see gh-17004")
|
||||
elif dtype == complex:
|
||||
# https://github.com/scipy/scipy/issues/11406
|
||||
pytest.skip("PROPACK unsupported for complex dtype")
|
||||
|
||||
rng = np.random.default_rng(179847540)
|
||||
# create a 10x10 singular matrix with a 4-dim null space
|
||||
dim = 4
|
||||
size = 10
|
||||
x = rng.random((size, size-dim))
|
||||
y = x[:, :dim] * rng.random(dim)
|
||||
mat = np.hstack((x, y))
|
||||
mat = mat.astype(dtype)
|
||||
|
||||
nz = null_space(mat)
|
||||
assert_equal(nz.shape[1], dim)
|
||||
|
||||
# Tolerances atol and rtol adjusted to pass np.float32
|
||||
# Use non-sparse svd
|
||||
u, s, vh = svd(mat)
|
||||
# Singular values are 0:
|
||||
assert_allclose(s[-dim:], 0, atol=1e-6, rtol=1e0)
|
||||
# Smallest right singular vectors in null space:
|
||||
assert_allclose(mat @ vh[-dim:, :].T, 0, atol=1e-6, rtol=1e0)
|
||||
|
||||
# Smallest singular values should be 0
|
||||
sp_mat = csc_matrix(mat)
|
||||
su, ss, svh = svds(sp_mat, k=dim, which='SM', solver=self.solver)
|
||||
# Smallest dim singular values are 0:
|
||||
assert_allclose(ss, 0, atol=1e-5, rtol=1e0)
|
||||
# Smallest singular vectors via svds in null space:
|
||||
n, m = mat.shape
|
||||
if n < m: # else the assert fails with some libraries unclear why
|
||||
assert_allclose(sp_mat.transpose() @ su, 0, atol=1e-5, rtol=1e0)
|
||||
assert_allclose(sp_mat @ svh.T, 0, atol=1e-5, rtol=1e0)
|
||||
|
||||
# --- Perform tests with each solver ---
|
||||
|
||||
|
||||
class Test_SVDS_once():
|
||||
@pytest.mark.parametrize("solver", ['ekki', object])
|
||||
def test_svds_input_validation_solver(self, solver):
|
||||
message = "solver must be one of"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.ones((3, 4)), k=2, solver=solver)
|
||||
|
||||
|
||||
class Test_SVDS_ARPACK(SVDSCommonTests):
|
||||
|
||||
def setup_method(self):
|
||||
self.solver = 'arpack'
|
||||
|
||||
@pytest.mark.parametrize("ncv", list(range(-1, 8)) + [4.5, "5"])
|
||||
def test_svds_input_validation_ncv_1(self, ncv):
|
||||
rng = np.random.default_rng(0)
|
||||
A = rng.random((6, 7))
|
||||
k = 3
|
||||
if ncv in {4, 5}:
|
||||
u, s, vh = svds(A, k=k, ncv=ncv, solver=self.solver)
|
||||
# partial decomposition, so don't check that u@diag(s)@vh=A;
|
||||
# do check that scipy.sparse.linalg.svds ~ scipy.linalg.svd
|
||||
_check_svds(A, k, u, s, vh)
|
||||
else:
|
||||
message = ("`ncv` must be an integer satisfying")
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(A, k=k, ncv=ncv, solver=self.solver)
|
||||
|
||||
def test_svds_input_validation_ncv_2(self):
|
||||
# I think the stack trace is reasonable when `ncv` can't be converted
|
||||
# to an int.
|
||||
message = "int() argument must be a"
|
||||
with pytest.raises(TypeError, match=re.escape(message)):
|
||||
svds(np.eye(10), ncv=[], solver=self.solver)
|
||||
|
||||
message = "invalid literal for int()"
|
||||
with pytest.raises(ValueError, match=message):
|
||||
svds(np.eye(10), ncv="hi", solver=self.solver)
|
||||
|
||||
# I can't see a robust relationship between `ncv` and relevant outputs
|
||||
# (e.g. accuracy, time), so no test of the parameter.
|
||||
|
||||
|
||||
class Test_SVDS_LOBPCG(SVDSCommonTests):
|
||||
|
||||
def setup_method(self):
|
||||
self.solver = 'lobpcg'
|
||||
|
||||
def test_svd_random_state_3(self):
|
||||
pytest.xfail("LOBPCG is having trouble with accuracy.")
|
||||
|
||||
|
||||
class Test_SVDS_PROPACK(SVDSCommonTests):
|
||||
|
||||
def setup_method(self):
|
||||
self.solver = 'propack'
|
||||
|
||||
def test_svd_LM_ones_matrix(self):
|
||||
message = ("PROPACK does not return orthonormal singular vectors "
|
||||
"associated with zero singular values.")
|
||||
# There are some other issues with this matrix of all ones, e.g.
|
||||
# `which='sm'` and `k=1` returns the largest singular value
|
||||
pytest.xfail(message)
|
||||
|
||||
def test_svd_LM_zeros_matrix(self):
|
||||
message = ("PROPACK does not return orthonormal singular vectors "
|
||||
"associated with zero singular values.")
|
||||
pytest.xfail(message)
|
||||
Reference in New Issue
Block a user