Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_svdp.py: 19%
107 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1"""
2Python wrapper for PROPACK
3--------------------------
5PROPACK is a collection of Fortran routines for iterative computation
6of partial SVDs of large matrices or linear operators.
8Based on BSD licensed pypropack project:
9 http://github.com/jakevdp/pypropack
10 Author: Jake Vanderplas <vanderplas@astro.washington.edu>
12PROPACK source is BSD licensed, and available at
13 http://soi.stanford.edu/~rmunk/PROPACK/
14"""
16__all__ = ['_svdp']
18import numpy as np
20from scipy._lib._util import check_random_state
21from scipy.sparse.linalg import aslinearoperator
22from scipy.linalg import LinAlgError
24from ._propack import _spropack # type: ignore[attr-defined]
25from ._propack import _dpropack # type: ignore[attr-defined]
26from ._propack import _cpropack # type: ignore[attr-defined]
27from ._propack import _zpropack # type: ignore[attr-defined]
30_lansvd_dict = {
31 'f': _spropack.slansvd,
32 'd': _dpropack.dlansvd,
33 'F': _cpropack.clansvd,
34 'D': _zpropack.zlansvd,
35}
38_lansvd_irl_dict = {
39 'f': _spropack.slansvd_irl,
40 'd': _dpropack.dlansvd_irl,
41 'F': _cpropack.clansvd_irl,
42 'D': _zpropack.zlansvd_irl,
43}
45_which_converter = {
46 'LM': 'L',
47 'SM': 'S',
48}
51class _AProd:
52 """
53 Wrapper class for linear operator
55 The call signature of the __call__ method matches the callback of
56 the PROPACK routines.
57 """
58 def __init__(self, A):
59 try:
60 self.A = aslinearoperator(A)
61 except TypeError:
62 self.A = aslinearoperator(np.asarray(A))
64 def __call__(self, transa, m, n, x, y, sparm, iparm):
65 if transa == 'n':
66 y[:] = self.A.matvec(x)
67 else:
68 y[:] = self.A.rmatvec(x)
70 @property
71 def shape(self):
72 return self.A.shape
74 @property
75 def dtype(self):
76 try:
77 return self.A.dtype
78 except AttributeError:
79 return self.A.matvec(np.zeros(self.A.shape[1])).dtype
82def _svdp(A, k, which='LM', irl_mode=True, kmax=None,
83 compute_u=True, compute_v=True, v0=None, full_output=False, tol=0,
84 delta=None, eta=None, anorm=0, cgs=False, elr=True,
85 min_relgap=0.002, shifts=None, maxiter=None, random_state=None):
86 """
87 Compute the singular value decomposition of a linear operator using PROPACK
89 Parameters
90 ----------
91 A : array_like, sparse matrix, or LinearOperator
92 Operator for which SVD will be computed. If `A` is a LinearOperator
93 object, it must define both ``matvec`` and ``rmatvec`` methods.
94 k : int
95 Number of singular values/vectors to compute
96 which : {"LM", "SM"}
97 Which singular triplets to compute:
98 - 'LM': compute triplets corresponding to the `k` largest singular
99 values
100 - 'SM': compute triplets corresponding to the `k` smallest singular
101 values
102 `which='SM'` requires `irl_mode=True`. Computes largest singular
103 values by default.
104 irl_mode : bool, optional
105 If `True`, then compute SVD using IRL (implicitly restarted Lanczos)
106 mode. Default is `True`.
107 kmax : int, optional
108 Maximal number of iterations / maximal dimension of the Krylov
109 subspace. Default is ``10 * k``.
110 compute_u : bool, optional
111 If `True` (default) then compute left singular vectors, `u`.
112 compute_v : bool, optional
113 If `True` (default) then compute right singular vectors, `v`.
114 tol : float, optional
115 The desired relative accuracy for computed singular values.
116 If not specified, it will be set based on machine precision.
117 v0 : array_like, optional
118 Starting vector for iterations: must be of length ``A.shape[0]``.
119 If not specified, PROPACK will generate a starting vector.
120 full_output : bool, optional
121 If `True`, then return sigma_bound. Default is `False`.
122 delta : float, optional
123 Level of orthogonality to maintain between Lanczos vectors.
124 Default is set based on machine precision.
125 eta : float, optional
126 Orthogonality cutoff. During reorthogonalization, vectors with
127 component larger than `eta` along the Lanczos vector will be purged.
128 Default is set based on machine precision.
129 anorm : float, optional
130 Estimate of ``||A||``. Default is `0`.
131 cgs : bool, optional
132 If `True`, reorthogonalization is done using classical Gram-Schmidt.
133 If `False` (default), it is done using modified Gram-Schmidt.
134 elr : bool, optional
135 If `True` (default), then extended local orthogonality is enforced
136 when obtaining singular vectors.
137 min_relgap : float, optional
138 The smallest relative gap allowed between any shift in IRL mode.
139 Default is `0.001`. Accessed only if ``irl_mode=True``.
140 shifts : int, optional
141 Number of shifts per restart in IRL mode. Default is determined
142 to satisfy ``k <= min(kmax-shifts, m, n)``. Must be
143 >= 0, but choosing 0 might lead to performance degradation.
144 Accessed only if ``irl_mode=True``.
145 maxiter : int, optional
146 Maximum number of restarts in IRL mode. Default is `1000`.
147 Accessed only if ``irl_mode=True``.
148 random_state : {None, int, `numpy.random.Generator`,
149 `numpy.random.RandomState`}, optional
151 Pseudorandom number generator state used to generate resamples.
153 If `random_state` is ``None`` (or `np.random`), the
154 `numpy.random.RandomState` singleton is used.
155 If `random_state` is an int, a new ``RandomState`` instance is used,
156 seeded with `random_state`.
157 If `random_state` is already a ``Generator`` or ``RandomState``
158 instance then that instance is used.
160 Returns
161 -------
162 u : ndarray
163 The `k` largest (``which="LM"``) or smallest (``which="SM"``) left
164 singular vectors, ``shape == (A.shape[0], 3)``, returned only if
165 ``compute_u=True``.
166 sigma : ndarray
167 The top `k` singular values, ``shape == (k,)``
168 vt : ndarray
169 The `k` largest (``which="LM"``) or smallest (``which="SM"``) right
170 singular vectors, ``shape == (3, A.shape[1])``, returned only if
171 ``compute_v=True``.
172 sigma_bound : ndarray
173 the error bounds on the singular values sigma, returned only if
174 ``full_output=True``.
176 """
177 random_state = check_random_state(random_state)
179 which = which.upper()
180 if which not in {'LM', 'SM'}:
181 raise ValueError("`which` must be either 'LM' or 'SM'")
182 if not irl_mode and which == 'SM':
183 raise ValueError("`which`='SM' requires irl_mode=True")
185 aprod = _AProd(A)
186 typ = aprod.dtype.char
188 try:
189 lansvd_irl = _lansvd_irl_dict[typ]
190 lansvd = _lansvd_dict[typ]
191 except KeyError:
192 # work with non-supported types using native system precision
193 if np.iscomplexobj(np.empty(0, dtype=typ)):
194 typ = np.dtype(complex).char
195 else:
196 typ = np.dtype(float).char
197 lansvd_irl = _lansvd_irl_dict[typ]
198 lansvd = _lansvd_dict[typ]
200 m, n = aprod.shape
201 if (k < 1) or (k > min(m, n)):
202 raise ValueError("k must be positive and not greater than m or n")
204 if kmax is None:
205 kmax = 10*k
206 if maxiter is None:
207 maxiter = 1000
209 # guard against unnecessarily large kmax
210 kmax = min(m + 1, n + 1, kmax)
211 if kmax < k:
212 raise ValueError(
213 "kmax must be greater than or equal to k, "
214 f"but kmax ({kmax}) < k ({k})")
216 # convert python args to fortran args
217 jobu = 'y' if compute_u else 'n'
218 jobv = 'y' if compute_v else 'n'
220 # these will be the output arrays
221 u = np.zeros((m, kmax + 1), order='F', dtype=typ)
222 v = np.zeros((n, kmax), order='F', dtype=typ)
224 # Specify the starting vector. if v0 is all zero, PROPACK will generate
225 # a random starting vector: the random seed cannot be controlled in that
226 # case, so we'll instead use numpy to generate a random vector
227 if v0 is None:
228 u[:, 0] = random_state.uniform(size=m)
229 if np.iscomplexobj(np.empty(0, dtype=typ)): # complex type
230 u[:, 0] += 1j * random_state.uniform(size=m)
231 else:
232 try:
233 u[:, 0] = v0
234 except ValueError:
235 raise ValueError(f"v0 must be of length {m}")
237 # process options for the fit
238 if delta is None:
239 delta = np.sqrt(np.finfo(typ).eps)
240 if eta is None:
241 eta = np.finfo(typ).eps ** 0.75
243 if irl_mode:
244 doption = np.array((delta, eta, anorm, min_relgap), dtype=typ.lower())
246 # validate or find default shifts
247 if shifts is None:
248 shifts = kmax - k
249 if k > min(kmax - shifts, m, n):
250 raise ValueError('shifts must satisfy '
251 'k <= min(kmax-shifts, m, n)!')
252 elif shifts < 0:
253 raise ValueError('shifts must be >= 0!')
255 else:
256 doption = np.array((delta, eta, anorm), dtype=typ.lower())
258 ioption = np.array((int(bool(cgs)), int(bool(elr))), dtype='i')
260 # If computing `u` or `v` (left and right singular vectors,
261 # respectively), `blocksize` controls how large a fraction of the
262 # work is done via fast BLAS level 3 operations. A larger blocksize
263 # may lead to faster computation at the expense of greater memory
264 # consumption. `blocksize` must be ``>= 1``. Choosing blocksize
265 # of 16, but docs don't specify; it's almost surely a
266 # power of 2.
267 blocksize = 16
269 # Determine lwork & liwork:
270 # the required lengths are specified in the PROPACK documentation
271 if compute_u or compute_v:
272 lwork = m + n + 9*kmax + 5*kmax*kmax + 4 + max(
273 3*kmax*kmax + 4*kmax + 4,
274 blocksize*max(m, n))
275 liwork = 8*kmax
276 else:
277 lwork = m + n + 9*kmax + 2*kmax*kmax + 4 + max(m + n, 4*kmax + 4)
278 liwork = 2*kmax + 1
279 work = np.empty(lwork, dtype=typ.lower())
280 iwork = np.empty(liwork, dtype=np.int32)
282 # dummy arguments: these are passed to aprod, and not used in this wrapper
283 dparm = np.empty(1, dtype=typ.lower())
284 iparm = np.empty(1, dtype=np.int32)
286 if typ.isupper():
287 # PROPACK documentation is unclear on the required length of zwork.
288 # Use the same length Julia's wrapper uses
289 # see https://github.com/JuliaSmoothOptimizers/PROPACK.jl/
290 zwork = np.empty(m + n + 32*m, dtype=typ)
291 works = work, zwork, iwork
292 else:
293 works = work, iwork
295 if irl_mode:
296 u, sigma, bnd, v, info = lansvd_irl(_which_converter[which], jobu,
297 jobv, m, n, shifts, k, maxiter,
298 aprod, u, v, tol, *works, doption,
299 ioption, dparm, iparm)
300 else:
301 u, sigma, bnd, v, info = lansvd(jobu, jobv, m, n, k, aprod, u, v, tol,
302 *works, doption, ioption, dparm, iparm)
304 if info > 0:
305 raise LinAlgError(
306 f"An invariant subspace of dimension {info} was found.")
307 elif info < 0:
308 raise LinAlgError(
309 f"k={k} singular triplets did not converge within "
310 f"kmax={kmax} iterations")
312 # info == 0: The K largest (or smallest) singular triplets were computed
313 # successfully!
315 return u[:, :k], sigma, v[:, :k].conj().T, bnd