Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py: 12%
639 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1"""
2Find a few eigenvectors and eigenvalues of a matrix.
5Uses ARPACK: https://github.com/opencollab/arpack-ng
7"""
8# Wrapper implementation notes
9#
10# ARPACK Entry Points
11# -------------------
12# The entry points to ARPACK are
13# - (s,d)seupd : single and double precision symmetric matrix
14# - (s,d,c,z)neupd: single,double,complex,double complex general matrix
15# This wrapper puts the *neupd (general matrix) interfaces in eigs()
16# and the *seupd (symmetric matrix) in eigsh().
17# There is no specialized interface for complex Hermitian matrices.
18# To find eigenvalues of a complex Hermitian matrix you
19# may use eigsh(), but eigsh() will simply call eigs()
20# and return the real part of the eigenvalues thus obtained.
22# Number of eigenvalues returned and complex eigenvalues
23# ------------------------------------------------------
24# The ARPACK nonsymmetric real and double interface (s,d)naupd return
25# eigenvalues and eigenvectors in real (float,double) arrays.
26# Since the eigenvalues and eigenvectors are, in general, complex
27# ARPACK puts the real and imaginary parts in consecutive entries
28# in real-valued arrays. This wrapper puts the real entries
29# into complex data types and attempts to return the requested eigenvalues
30# and eigenvectors.
33# Solver modes
34# ------------
35# ARPACK and handle shifted and shift-inverse computations
36# for eigenvalues by providing a shift (sigma) and a solver.
38import numpy as np
39import warnings
40from scipy.sparse.linalg._interface import aslinearoperator, LinearOperator
41from scipy.sparse import eye, issparse
42from scipy.linalg import eig, eigh, lu_factor, lu_solve
43from scipy.sparse._sputils import isdense, is_pydata_spmatrix
44from scipy.sparse.linalg import gmres, splu
45from scipy._lib._util import _aligned_zeros
46from scipy._lib._threadsafety import ReentrancyLock
48from . import _arpack
49arpack_int = _arpack.timing.nbx.dtype
51__docformat__ = "restructuredtext en"
53__all__ = ['eigs', 'eigsh', 'ArpackError', 'ArpackNoConvergence']
56_type_conv = {'f': 's', 'd': 'd', 'F': 'c', 'D': 'z'}
57_ndigits = {'f': 5, 'd': 12, 'F': 5, 'D': 12}
59DNAUPD_ERRORS = {
60 0: "Normal exit.",
61 1: "Maximum number of iterations taken. "
62 "All possible eigenvalues of OP has been found. IPARAM(5) "
63 "returns the number of wanted converged Ritz values.",
64 2: "No longer an informational error. Deprecated starting "
65 "with release 2 of ARPACK.",
66 3: "No shifts could be applied during a cycle of the "
67 "Implicitly restarted Arnoldi iteration. One possibility "
68 "is to increase the size of NCV relative to NEV. ",
69 -1: "N must be positive.",
70 -2: "NEV must be positive.",
71 -3: "NCV-NEV >= 2 and less than or equal to N.",
72 -4: "The maximum number of Arnoldi update iterations allowed "
73 "must be greater than zero.",
74 -5: " WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'",
75 -6: "BMAT must be one of 'I' or 'G'.",
76 -7: "Length of private work array WORKL is not sufficient.",
77 -8: "Error return from LAPACK eigenvalue calculation;",
78 -9: "Starting vector is zero.",
79 -10: "IPARAM(7) must be 1,2,3,4.",
80 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
81 -12: "IPARAM(1) must be equal to 0 or 1.",
82 -13: "NEV and WHICH = 'BE' are incompatible.",
83 -9999: "Could not build an Arnoldi factorization. "
84 "IPARAM(5) returns the size of the current Arnoldi "
85 "factorization. The user is advised to check that "
86 "enough workspace and array storage has been allocated."
87}
89SNAUPD_ERRORS = DNAUPD_ERRORS
91ZNAUPD_ERRORS = DNAUPD_ERRORS.copy()
92ZNAUPD_ERRORS[-10] = "IPARAM(7) must be 1,2,3."
94CNAUPD_ERRORS = ZNAUPD_ERRORS
96DSAUPD_ERRORS = {
97 0: "Normal exit.",
98 1: "Maximum number of iterations taken. "
99 "All possible eigenvalues of OP has been found.",
100 2: "No longer an informational error. Deprecated starting with "
101 "release 2 of ARPACK.",
102 3: "No shifts could be applied during a cycle of the Implicitly "
103 "restarted Arnoldi iteration. One possibility is to increase "
104 "the size of NCV relative to NEV. ",
105 -1: "N must be positive.",
106 -2: "NEV must be positive.",
107 -3: "NCV must be greater than NEV and less than or equal to N.",
108 -4: "The maximum number of Arnoldi update iterations allowed "
109 "must be greater than zero.",
110 -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.",
111 -6: "BMAT must be one of 'I' or 'G'.",
112 -7: "Length of private work array WORKL is not sufficient.",
113 -8: "Error return from trid. eigenvalue calculation; "
114 "Informational error from LAPACK routine dsteqr .",
115 -9: "Starting vector is zero.",
116 -10: "IPARAM(7) must be 1,2,3,4,5.",
117 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
118 -12: "IPARAM(1) must be equal to 0 or 1.",
119 -13: "NEV and WHICH = 'BE' are incompatible. ",
120 -9999: "Could not build an Arnoldi factorization. "
121 "IPARAM(5) returns the size of the current Arnoldi "
122 "factorization. The user is advised to check that "
123 "enough workspace and array storage has been allocated.",
124}
126SSAUPD_ERRORS = DSAUPD_ERRORS
128DNEUPD_ERRORS = {
129 0: "Normal exit.",
130 1: "The Schur form computed by LAPACK routine dlahqr "
131 "could not be reordered by LAPACK routine dtrsen. "
132 "Re-enter subroutine dneupd with IPARAM(5)NCV and "
133 "increase the size of the arrays DR and DI to have "
134 "dimension at least dimension NCV and allocate at least NCV "
135 "columns for Z. NOTE: Not necessary if Z and V share "
136 "the same space. Please notify the authors if this error"
137 "occurs.",
138 -1: "N must be positive.",
139 -2: "NEV must be positive.",
140 -3: "NCV-NEV >= 2 and less than or equal to N.",
141 -5: "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'",
142 -6: "BMAT must be one of 'I' or 'G'.",
143 -7: "Length of private work WORKL array is not sufficient.",
144 -8: "Error return from calculation of a real Schur form. "
145 "Informational error from LAPACK routine dlahqr .",
146 -9: "Error return from calculation of eigenvectors. "
147 "Informational error from LAPACK routine dtrevc.",
148 -10: "IPARAM(7) must be 1,2,3,4.",
149 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
150 -12: "HOWMNY = 'S' not yet implemented",
151 -13: "HOWMNY must be one of 'A' or 'P' if RVEC = .true.",
152 -14: "DNAUPD did not find any eigenvalues to sufficient "
153 "accuracy.",
154 -15: "DNEUPD got a different count of the number of converged "
155 "Ritz values than DNAUPD got. This indicates the user "
156 "probably made an error in passing data from DNAUPD to "
157 "DNEUPD or that the data was modified before entering "
158 "DNEUPD",
159}
161SNEUPD_ERRORS = DNEUPD_ERRORS.copy()
162SNEUPD_ERRORS[1] = ("The Schur form computed by LAPACK routine slahqr "
163 "could not be reordered by LAPACK routine strsen . "
164 "Re-enter subroutine dneupd with IPARAM(5)=NCV and "
165 "increase the size of the arrays DR and DI to have "
166 "dimension at least dimension NCV and allocate at least "
167 "NCV columns for Z. NOTE: Not necessary if Z and V share "
168 "the same space. Please notify the authors if this error "
169 "occurs.")
170SNEUPD_ERRORS[-14] = ("SNAUPD did not find any eigenvalues to sufficient "
171 "accuracy.")
172SNEUPD_ERRORS[-15] = ("SNEUPD got a different count of the number of "
173 "converged Ritz values than SNAUPD got. This indicates "
174 "the user probably made an error in passing data from "
175 "SNAUPD to SNEUPD or that the data was modified before "
176 "entering SNEUPD")
178ZNEUPD_ERRORS = {0: "Normal exit.",
179 1: "The Schur form computed by LAPACK routine csheqr "
180 "could not be reordered by LAPACK routine ztrsen. "
181 "Re-enter subroutine zneupd with IPARAM(5)=NCV and "
182 "increase the size of the array D to have "
183 "dimension at least dimension NCV and allocate at least "
184 "NCV columns for Z. NOTE: Not necessary if Z and V share "
185 "the same space. Please notify the authors if this error "
186 "occurs.",
187 -1: "N must be positive.",
188 -2: "NEV must be positive.",
189 -3: "NCV-NEV >= 1 and less than or equal to N.",
190 -5: "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'",
191 -6: "BMAT must be one of 'I' or 'G'.",
192 -7: "Length of private work WORKL array is not sufficient.",
193 -8: "Error return from LAPACK eigenvalue calculation. "
194 "This should never happened.",
195 -9: "Error return from calculation of eigenvectors. "
196 "Informational error from LAPACK routine ztrevc.",
197 -10: "IPARAM(7) must be 1,2,3",
198 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
199 -12: "HOWMNY = 'S' not yet implemented",
200 -13: "HOWMNY must be one of 'A' or 'P' if RVEC = .true.",
201 -14: "ZNAUPD did not find any eigenvalues to sufficient "
202 "accuracy.",
203 -15: "ZNEUPD got a different count of the number of "
204 "converged Ritz values than ZNAUPD got. This "
205 "indicates the user probably made an error in passing "
206 "data from ZNAUPD to ZNEUPD or that the data was "
207 "modified before entering ZNEUPD"
208 }
210CNEUPD_ERRORS = ZNEUPD_ERRORS.copy()
211CNEUPD_ERRORS[-14] = ("CNAUPD did not find any eigenvalues to sufficient "
212 "accuracy.")
213CNEUPD_ERRORS[-15] = ("CNEUPD got a different count of the number of "
214 "converged Ritz values than CNAUPD got. This indicates "
215 "the user probably made an error in passing data from "
216 "CNAUPD to CNEUPD or that the data was modified before "
217 "entering CNEUPD")
219DSEUPD_ERRORS = {
220 0: "Normal exit.",
221 -1: "N must be positive.",
222 -2: "NEV must be positive.",
223 -3: "NCV must be greater than NEV and less than or equal to N.",
224 -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.",
225 -6: "BMAT must be one of 'I' or 'G'.",
226 -7: "Length of private work WORKL array is not sufficient.",
227 -8: ("Error return from trid. eigenvalue calculation; "
228 "Information error from LAPACK routine dsteqr."),
229 -9: "Starting vector is zero.",
230 -10: "IPARAM(7) must be 1,2,3,4,5.",
231 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
232 -12: "NEV and WHICH = 'BE' are incompatible.",
233 -14: "DSAUPD did not find any eigenvalues to sufficient accuracy.",
234 -15: "HOWMNY must be one of 'A' or 'S' if RVEC = .true.",
235 -16: "HOWMNY = 'S' not yet implemented",
236 -17: ("DSEUPD got a different count of the number of converged "
237 "Ritz values than DSAUPD got. This indicates the user "
238 "probably made an error in passing data from DSAUPD to "
239 "DSEUPD or that the data was modified before entering "
240 "DSEUPD.")
241}
243SSEUPD_ERRORS = DSEUPD_ERRORS.copy()
244SSEUPD_ERRORS[-14] = ("SSAUPD did not find any eigenvalues "
245 "to sufficient accuracy.")
246SSEUPD_ERRORS[-17] = ("SSEUPD got a different count of the number of "
247 "converged "
248 "Ritz values than SSAUPD got. This indicates the user "
249 "probably made an error in passing data from SSAUPD to "
250 "SSEUPD or that the data was modified before entering "
251 "SSEUPD.")
253_SAUPD_ERRORS = {'d': DSAUPD_ERRORS,
254 's': SSAUPD_ERRORS}
255_NAUPD_ERRORS = {'d': DNAUPD_ERRORS,
256 's': SNAUPD_ERRORS,
257 'z': ZNAUPD_ERRORS,
258 'c': CNAUPD_ERRORS}
259_SEUPD_ERRORS = {'d': DSEUPD_ERRORS,
260 's': SSEUPD_ERRORS}
261_NEUPD_ERRORS = {'d': DNEUPD_ERRORS,
262 's': SNEUPD_ERRORS,
263 'z': ZNEUPD_ERRORS,
264 'c': CNEUPD_ERRORS}
266# accepted values of parameter WHICH in _SEUPD
267_SEUPD_WHICH = ['LM', 'SM', 'LA', 'SA', 'BE']
269# accepted values of parameter WHICH in _NAUPD
270_NEUPD_WHICH = ['LM', 'SM', 'LR', 'SR', 'LI', 'SI']
273class ArpackError(RuntimeError):
274 """
275 ARPACK error
276 """
278 def __init__(self, info, infodict=_NAUPD_ERRORS):
279 msg = infodict.get(info, "Unknown error")
280 RuntimeError.__init__(self, "ARPACK error %d: %s" % (info, msg))
283class ArpackNoConvergence(ArpackError):
284 """
285 ARPACK iteration did not converge
287 Attributes
288 ----------
289 eigenvalues : ndarray
290 Partial result. Converged eigenvalues.
291 eigenvectors : ndarray
292 Partial result. Converged eigenvectors.
294 """
296 def __init__(self, msg, eigenvalues, eigenvectors):
297 ArpackError.__init__(self, -1, {-1: msg})
298 self.eigenvalues = eigenvalues
299 self.eigenvectors = eigenvectors
302def choose_ncv(k):
303 """
304 Choose number of lanczos vectors based on target number
305 of singular/eigen values and vectors to compute, k.
306 """
307 return max(2 * k + 1, 20)
310class _ArpackParams:
311 def __init__(self, n, k, tp, mode=1, sigma=None,
312 ncv=None, v0=None, maxiter=None, which="LM", tol=0):
313 if k <= 0:
314 raise ValueError("k must be positive, k=%d" % k)
316 if maxiter is None:
317 maxiter = n * 10
318 if maxiter <= 0:
319 raise ValueError("maxiter must be positive, maxiter=%d" % maxiter)
321 if tp not in 'fdFD':
322 raise ValueError("matrix type must be 'f', 'd', 'F', or 'D'")
324 if v0 is not None:
325 # ARPACK overwrites its initial resid, make a copy
326 self.resid = np.array(v0, copy=True)
327 info = 1
328 else:
329 # ARPACK will use a random initial vector.
330 self.resid = np.zeros(n, tp)
331 info = 0
333 if sigma is None:
334 #sigma not used
335 self.sigma = 0
336 else:
337 self.sigma = sigma
339 if ncv is None:
340 ncv = choose_ncv(k)
341 ncv = min(ncv, n)
343 self.v = np.zeros((n, ncv), tp) # holds Ritz vectors
344 self.iparam = np.zeros(11, arpack_int)
346 # set solver mode and parameters
347 ishfts = 1
348 self.mode = mode
349 self.iparam[0] = ishfts
350 self.iparam[2] = maxiter
351 self.iparam[3] = 1
352 self.iparam[6] = mode
354 self.n = n
355 self.tol = tol
356 self.k = k
357 self.maxiter = maxiter
358 self.ncv = ncv
359 self.which = which
360 self.tp = tp
361 self.info = info
363 self.converged = False
364 self.ido = 0
366 def _raise_no_convergence(self):
367 msg = "No convergence (%d iterations, %d/%d eigenvectors converged)"
368 k_ok = self.iparam[4]
369 num_iter = self.iparam[2]
370 try:
371 ev, vec = self.extract(True)
372 except ArpackError as err:
373 msg = f"{msg} [{err}]"
374 ev = np.zeros((0,))
375 vec = np.zeros((self.n, 0))
376 k_ok = 0
377 raise ArpackNoConvergence(msg % (num_iter, k_ok, self.k), ev, vec)
380class _SymmetricArpackParams(_ArpackParams):
381 def __init__(self, n, k, tp, matvec, mode=1, M_matvec=None,
382 Minv_matvec=None, sigma=None,
383 ncv=None, v0=None, maxiter=None, which="LM", tol=0):
384 # The following modes are supported:
385 # mode = 1:
386 # Solve the standard eigenvalue problem:
387 # A*x = lambda*x :
388 # A - symmetric
389 # Arguments should be
390 # matvec = left multiplication by A
391 # M_matvec = None [not used]
392 # Minv_matvec = None [not used]
393 #
394 # mode = 2:
395 # Solve the general eigenvalue problem:
396 # A*x = lambda*M*x
397 # A - symmetric
398 # M - symmetric positive definite
399 # Arguments should be
400 # matvec = left multiplication by A
401 # M_matvec = left multiplication by M
402 # Minv_matvec = left multiplication by M^-1
403 #
404 # mode = 3:
405 # Solve the general eigenvalue problem in shift-invert mode:
406 # A*x = lambda*M*x
407 # A - symmetric
408 # M - symmetric positive semi-definite
409 # Arguments should be
410 # matvec = None [not used]
411 # M_matvec = left multiplication by M
412 # or None, if M is the identity
413 # Minv_matvec = left multiplication by [A-sigma*M]^-1
414 #
415 # mode = 4:
416 # Solve the general eigenvalue problem in Buckling mode:
417 # A*x = lambda*AG*x
418 # A - symmetric positive semi-definite
419 # AG - symmetric indefinite
420 # Arguments should be
421 # matvec = left multiplication by A
422 # M_matvec = None [not used]
423 # Minv_matvec = left multiplication by [A-sigma*AG]^-1
424 #
425 # mode = 5:
426 # Solve the general eigenvalue problem in Cayley-transformed mode:
427 # A*x = lambda*M*x
428 # A - symmetric
429 # M - symmetric positive semi-definite
430 # Arguments should be
431 # matvec = left multiplication by A
432 # M_matvec = left multiplication by M
433 # or None, if M is the identity
434 # Minv_matvec = left multiplication by [A-sigma*M]^-1
435 if mode == 1:
436 if matvec is None:
437 raise ValueError("matvec must be specified for mode=1")
438 if M_matvec is not None:
439 raise ValueError("M_matvec cannot be specified for mode=1")
440 if Minv_matvec is not None:
441 raise ValueError("Minv_matvec cannot be specified for mode=1")
443 self.OP = matvec
444 self.B = lambda x: x
445 self.bmat = 'I'
446 elif mode == 2:
447 if matvec is None:
448 raise ValueError("matvec must be specified for mode=2")
449 if M_matvec is None:
450 raise ValueError("M_matvec must be specified for mode=2")
451 if Minv_matvec is None:
452 raise ValueError("Minv_matvec must be specified for mode=2")
454 self.OP = lambda x: Minv_matvec(matvec(x))
455 self.OPa = Minv_matvec
456 self.OPb = matvec
457 self.B = M_matvec
458 self.bmat = 'G'
459 elif mode == 3:
460 if matvec is not None:
461 raise ValueError("matvec must not be specified for mode=3")
462 if Minv_matvec is None:
463 raise ValueError("Minv_matvec must be specified for mode=3")
465 if M_matvec is None:
466 self.OP = Minv_matvec
467 self.OPa = Minv_matvec
468 self.B = lambda x: x
469 self.bmat = 'I'
470 else:
471 self.OP = lambda x: Minv_matvec(M_matvec(x))
472 self.OPa = Minv_matvec
473 self.B = M_matvec
474 self.bmat = 'G'
475 elif mode == 4:
476 if matvec is None:
477 raise ValueError("matvec must be specified for mode=4")
478 if M_matvec is not None:
479 raise ValueError("M_matvec must not be specified for mode=4")
480 if Minv_matvec is None:
481 raise ValueError("Minv_matvec must be specified for mode=4")
482 self.OPa = Minv_matvec
483 self.OP = lambda x: self.OPa(matvec(x))
484 self.B = matvec
485 self.bmat = 'G'
486 elif mode == 5:
487 if matvec is None:
488 raise ValueError("matvec must be specified for mode=5")
489 if Minv_matvec is None:
490 raise ValueError("Minv_matvec must be specified for mode=5")
492 self.OPa = Minv_matvec
493 self.A_matvec = matvec
495 if M_matvec is None:
496 self.OP = lambda x: Minv_matvec(matvec(x) + sigma * x)
497 self.B = lambda x: x
498 self.bmat = 'I'
499 else:
500 self.OP = lambda x: Minv_matvec(matvec(x)
501 + sigma * M_matvec(x))
502 self.B = M_matvec
503 self.bmat = 'G'
504 else:
505 raise ValueError("mode=%i not implemented" % mode)
507 if which not in _SEUPD_WHICH:
508 raise ValueError("which must be one of %s"
509 % ' '.join(_SEUPD_WHICH))
510 if k >= n:
511 raise ValueError("k must be less than ndim(A), k=%d" % k)
513 _ArpackParams.__init__(self, n, k, tp, mode, sigma,
514 ncv, v0, maxiter, which, tol)
516 if self.ncv > n or self.ncv <= k:
517 raise ValueError("ncv must be k<ncv<=n, ncv=%s" % self.ncv)
519 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
520 self.workd = _aligned_zeros(3 * n, self.tp)
521 self.workl = _aligned_zeros(self.ncv * (self.ncv + 8), self.tp)
523 ltr = _type_conv[self.tp]
524 if ltr not in ["s", "d"]:
525 raise ValueError("Input matrix is not real-valued.")
527 self._arpack_solver = _arpack.__dict__[ltr + 'saupd']
528 self._arpack_extract = _arpack.__dict__[ltr + 'seupd']
530 self.iterate_infodict = _SAUPD_ERRORS[ltr]
531 self.extract_infodict = _SEUPD_ERRORS[ltr]
533 self.ipntr = np.zeros(11, arpack_int)
535 def iterate(self):
536 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info = \
537 self._arpack_solver(self.ido, self.bmat, self.which, self.k,
538 self.tol, self.resid, self.v, self.iparam,
539 self.ipntr, self.workd, self.workl, self.info)
541 xslice = slice(self.ipntr[0] - 1, self.ipntr[0] - 1 + self.n)
542 yslice = slice(self.ipntr[1] - 1, self.ipntr[1] - 1 + self.n)
543 if self.ido == -1:
544 # initialization
545 self.workd[yslice] = self.OP(self.workd[xslice])
546 elif self.ido == 1:
547 # compute y = Op*x
548 if self.mode == 1:
549 self.workd[yslice] = self.OP(self.workd[xslice])
550 elif self.mode == 2:
551 self.workd[xslice] = self.OPb(self.workd[xslice])
552 self.workd[yslice] = self.OPa(self.workd[xslice])
553 elif self.mode == 5:
554 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n)
555 Ax = self.A_matvec(self.workd[xslice])
556 self.workd[yslice] = self.OPa(Ax + (self.sigma *
557 self.workd[Bxslice]))
558 else:
559 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n)
560 self.workd[yslice] = self.OPa(self.workd[Bxslice])
561 elif self.ido == 2:
562 self.workd[yslice] = self.B(self.workd[xslice])
563 elif self.ido == 3:
564 raise ValueError("ARPACK requested user shifts. Assure ISHIFT==0")
565 else:
566 self.converged = True
568 if self.info == 0:
569 pass
570 elif self.info == 1:
571 self._raise_no_convergence()
572 else:
573 raise ArpackError(self.info, infodict=self.iterate_infodict)
575 def extract(self, return_eigenvectors):
576 rvec = return_eigenvectors
577 ierr = 0
578 howmny = 'A' # return all eigenvectors
579 sselect = np.zeros(self.ncv, 'int') # unused
580 d, z, ierr = self._arpack_extract(rvec, howmny, sselect, self.sigma,
581 self.bmat, self.which, self.k,
582 self.tol, self.resid, self.v,
583 self.iparam[0:7], self.ipntr,
584 self.workd[0:2 * self.n],
585 self.workl, ierr)
586 if ierr != 0:
587 raise ArpackError(ierr, infodict=self.extract_infodict)
588 k_ok = self.iparam[4]
589 d = d[:k_ok]
590 z = z[:, :k_ok]
592 if return_eigenvectors:
593 return d, z
594 else:
595 return d
598class _UnsymmetricArpackParams(_ArpackParams):
599 def __init__(self, n, k, tp, matvec, mode=1, M_matvec=None,
600 Minv_matvec=None, sigma=None,
601 ncv=None, v0=None, maxiter=None, which="LM", tol=0):
602 # The following modes are supported:
603 # mode = 1:
604 # Solve the standard eigenvalue problem:
605 # A*x = lambda*x
606 # A - square matrix
607 # Arguments should be
608 # matvec = left multiplication by A
609 # M_matvec = None [not used]
610 # Minv_matvec = None [not used]
611 #
612 # mode = 2:
613 # Solve the generalized eigenvalue problem:
614 # A*x = lambda*M*x
615 # A - square matrix
616 # M - symmetric, positive semi-definite
617 # Arguments should be
618 # matvec = left multiplication by A
619 # M_matvec = left multiplication by M
620 # Minv_matvec = left multiplication by M^-1
621 #
622 # mode = 3,4:
623 # Solve the general eigenvalue problem in shift-invert mode:
624 # A*x = lambda*M*x
625 # A - square matrix
626 # M - symmetric, positive semi-definite
627 # Arguments should be
628 # matvec = None [not used]
629 # M_matvec = left multiplication by M
630 # or None, if M is the identity
631 # Minv_matvec = left multiplication by [A-sigma*M]^-1
632 # if A is real and mode==3, use the real part of Minv_matvec
633 # if A is real and mode==4, use the imag part of Minv_matvec
634 # if A is complex and mode==3,
635 # use real and imag parts of Minv_matvec
636 if mode == 1:
637 if matvec is None:
638 raise ValueError("matvec must be specified for mode=1")
639 if M_matvec is not None:
640 raise ValueError("M_matvec cannot be specified for mode=1")
641 if Minv_matvec is not None:
642 raise ValueError("Minv_matvec cannot be specified for mode=1")
644 self.OP = matvec
645 self.B = lambda x: x
646 self.bmat = 'I'
647 elif mode == 2:
648 if matvec is None:
649 raise ValueError("matvec must be specified for mode=2")
650 if M_matvec is None:
651 raise ValueError("M_matvec must be specified for mode=2")
652 if Minv_matvec is None:
653 raise ValueError("Minv_matvec must be specified for mode=2")
655 self.OP = lambda x: Minv_matvec(matvec(x))
656 self.OPa = Minv_matvec
657 self.OPb = matvec
658 self.B = M_matvec
659 self.bmat = 'G'
660 elif mode in (3, 4):
661 if matvec is None:
662 raise ValueError("matvec must be specified "
663 "for mode in (3,4)")
664 if Minv_matvec is None:
665 raise ValueError("Minv_matvec must be specified "
666 "for mode in (3,4)")
668 self.matvec = matvec
669 if tp in 'DF': # complex type
670 if mode == 3:
671 self.OPa = Minv_matvec
672 else:
673 raise ValueError("mode=4 invalid for complex A")
674 else: # real type
675 if mode == 3:
676 self.OPa = lambda x: np.real(Minv_matvec(x))
677 else:
678 self.OPa = lambda x: np.imag(Minv_matvec(x))
679 if M_matvec is None:
680 self.B = lambda x: x
681 self.bmat = 'I'
682 self.OP = self.OPa
683 else:
684 self.B = M_matvec
685 self.bmat = 'G'
686 self.OP = lambda x: self.OPa(M_matvec(x))
687 else:
688 raise ValueError("mode=%i not implemented" % mode)
690 if which not in _NEUPD_WHICH:
691 raise ValueError("Parameter which must be one of %s"
692 % ' '.join(_NEUPD_WHICH))
693 if k >= n - 1:
694 raise ValueError("k must be less than ndim(A)-1, k=%d" % k)
696 _ArpackParams.__init__(self, n, k, tp, mode, sigma,
697 ncv, v0, maxiter, which, tol)
699 if self.ncv > n or self.ncv <= k + 1:
700 raise ValueError("ncv must be k+1<ncv<=n, ncv=%s" % self.ncv)
702 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
703 self.workd = _aligned_zeros(3 * n, self.tp)
704 self.workl = _aligned_zeros(3 * self.ncv * (self.ncv + 2), self.tp)
706 ltr = _type_conv[self.tp]
707 self._arpack_solver = _arpack.__dict__[ltr + 'naupd']
708 self._arpack_extract = _arpack.__dict__[ltr + 'neupd']
710 self.iterate_infodict = _NAUPD_ERRORS[ltr]
711 self.extract_infodict = _NEUPD_ERRORS[ltr]
713 self.ipntr = np.zeros(14, arpack_int)
715 if self.tp in 'FD':
716 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
717 self.rwork = _aligned_zeros(self.ncv, self.tp.lower())
718 else:
719 self.rwork = None
721 def iterate(self):
722 if self.tp in 'fd':
723 results = self._arpack_solver(self.ido, self.bmat, self.which, self.k,
724 self.tol, self.resid, self.v, self.iparam,
725 self.ipntr, self.workd, self.workl, self.info)
726 self.ido, self.tol, self.resid, self.v, \
727 self.iparam, self.ipntr, self.info = results
729 else:
730 results = self._arpack_solver(self.ido, self.bmat, self.which, self.k,
731 self.tol, self.resid, self.v, self.iparam,
732 self.ipntr, self.workd, self.workl,
733 self.rwork, self.info)
734 self.ido, self.tol, self.resid, self.v, \
735 self.iparam, self.ipntr, self.info = results
738 xslice = slice(self.ipntr[0] - 1, self.ipntr[0] - 1 + self.n)
739 yslice = slice(self.ipntr[1] - 1, self.ipntr[1] - 1 + self.n)
740 if self.ido == -1:
741 # initialization
742 self.workd[yslice] = self.OP(self.workd[xslice])
743 elif self.ido == 1:
744 # compute y = Op*x
745 if self.mode in (1, 2):
746 self.workd[yslice] = self.OP(self.workd[xslice])
747 else:
748 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n)
749 self.workd[yslice] = self.OPa(self.workd[Bxslice])
750 elif self.ido == 2:
751 self.workd[yslice] = self.B(self.workd[xslice])
752 elif self.ido == 3:
753 raise ValueError("ARPACK requested user shifts. Assure ISHIFT==0")
754 else:
755 self.converged = True
757 if self.info == 0:
758 pass
759 elif self.info == 1:
760 self._raise_no_convergence()
761 else:
762 raise ArpackError(self.info, infodict=self.iterate_infodict)
764 def extract(self, return_eigenvectors):
765 k, n = self.k, self.n
767 ierr = 0
768 howmny = 'A' # return all eigenvectors
769 sselect = np.zeros(self.ncv, 'int') # unused
770 sigmar = np.real(self.sigma)
771 sigmai = np.imag(self.sigma)
772 workev = np.zeros(3 * self.ncv, self.tp)
774 if self.tp in 'fd':
775 dr = np.zeros(k + 1, self.tp)
776 di = np.zeros(k + 1, self.tp)
777 zr = np.zeros((n, k + 1), self.tp)
778 dr, di, zr, ierr = \
779 self._arpack_extract(return_eigenvectors,
780 howmny, sselect, sigmar, sigmai, workev,
781 self.bmat, self.which, k, self.tol, self.resid,
782 self.v, self.iparam, self.ipntr,
783 self.workd, self.workl, self.info)
784 if ierr != 0:
785 raise ArpackError(ierr, infodict=self.extract_infodict)
786 nreturned = self.iparam[4] # number of good eigenvalues returned
788 # Build complex eigenvalues from real and imaginary parts
789 d = dr + 1.0j * di
791 # Arrange the eigenvectors: complex eigenvectors are stored as
792 # real,imaginary in consecutive columns
793 z = zr.astype(self.tp.upper())
795 # The ARPACK nonsymmetric real and double interface (s,d)naupd
796 # return eigenvalues and eigenvectors in real (float,double)
797 # arrays.
799 # Efficiency: this should check that return_eigenvectors == True
800 # before going through this construction.
801 if sigmai == 0:
802 i = 0
803 while i <= k:
804 # check if complex
805 if abs(d[i].imag) != 0:
806 # this is a complex conjugate pair with eigenvalues
807 # in consecutive columns
808 if i < k:
809 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1]
810 z[:, i + 1] = z[:, i].conjugate()
811 i += 1
812 else:
813 #last eigenvalue is complex: the imaginary part of
814 # the eigenvector has not been returned
815 #this can only happen if nreturned > k, so we'll
816 # throw out this case.
817 nreturned -= 1
818 i += 1
820 else:
821 # real matrix, mode 3 or 4, imag(sigma) is nonzero:
822 # see remark 3 in <s,d>neupd.f
823 # Build complex eigenvalues from real and imaginary parts
824 i = 0
825 while i <= k:
826 if abs(d[i].imag) == 0:
827 d[i] = np.dot(zr[:, i], self.matvec(zr[:, i]))
828 else:
829 if i < k:
830 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1]
831 z[:, i + 1] = z[:, i].conjugate()
832 d[i] = ((np.dot(zr[:, i],
833 self.matvec(zr[:, i]))
834 + np.dot(zr[:, i + 1],
835 self.matvec(zr[:, i + 1])))
836 + 1j * (np.dot(zr[:, i],
837 self.matvec(zr[:, i + 1]))
838 - np.dot(zr[:, i + 1],
839 self.matvec(zr[:, i]))))
840 d[i + 1] = d[i].conj()
841 i += 1
842 else:
843 #last eigenvalue is complex: the imaginary part of
844 # the eigenvector has not been returned
845 #this can only happen if nreturned > k, so we'll
846 # throw out this case.
847 nreturned -= 1
848 i += 1
850 # Now we have k+1 possible eigenvalues and eigenvectors
851 # Return the ones specified by the keyword "which"
853 if nreturned <= k:
854 # we got less or equal as many eigenvalues we wanted
855 d = d[:nreturned]
856 z = z[:, :nreturned]
857 else:
858 # we got one extra eigenvalue (likely a cc pair, but which?)
859 if self.mode in (1, 2):
860 rd = d
861 elif self.mode in (3, 4):
862 rd = 1 / (d - self.sigma)
864 if self.which in ['LR', 'SR']:
865 ind = np.argsort(rd.real)
866 elif self.which in ['LI', 'SI']:
867 # for LI,SI ARPACK returns largest,smallest
868 # abs(imaginary) (complex pairs come together)
869 ind = np.argsort(abs(rd.imag))
870 else:
871 ind = np.argsort(abs(rd))
873 if self.which in ['LR', 'LM', 'LI']:
874 ind = ind[-k:][::-1]
875 elif self.which in ['SR', 'SM', 'SI']:
876 ind = ind[:k]
878 d = d[ind]
879 z = z[:, ind]
880 else:
881 # complex is so much simpler...
882 d, z, ierr =\
883 self._arpack_extract(return_eigenvectors,
884 howmny, sselect, self.sigma, workev,
885 self.bmat, self.which, k, self.tol, self.resid,
886 self.v, self.iparam, self.ipntr,
887 self.workd, self.workl, self.rwork, ierr)
889 if ierr != 0:
890 raise ArpackError(ierr, infodict=self.extract_infodict)
892 k_ok = self.iparam[4]
893 d = d[:k_ok]
894 z = z[:, :k_ok]
896 if return_eigenvectors:
897 return d, z
898 else:
899 return d
902def _aslinearoperator_with_dtype(m):
903 m = aslinearoperator(m)
904 if not hasattr(m, 'dtype'):
905 x = np.zeros(m.shape[1])
906 m.dtype = (m * x).dtype
907 return m
910class SpLuInv(LinearOperator):
911 """
912 SpLuInv:
913 helper class to repeatedly solve M*x=b
914 using a sparse LU-decomposition of M
915 """
917 def __init__(self, M):
918 self.M_lu = splu(M)
919 self.shape = M.shape
920 self.dtype = M.dtype
921 self.isreal = not np.issubdtype(self.dtype, np.complexfloating)
923 def _matvec(self, x):
924 # careful here: splu.solve will throw away imaginary
925 # part of x if M is real
926 x = np.asarray(x)
927 if self.isreal and np.issubdtype(x.dtype, np.complexfloating):
928 return (self.M_lu.solve(np.real(x).astype(self.dtype))
929 + 1j * self.M_lu.solve(np.imag(x).astype(self.dtype)))
930 else:
931 return self.M_lu.solve(x.astype(self.dtype))
934class LuInv(LinearOperator):
935 """
936 LuInv:
937 helper class to repeatedly solve M*x=b
938 using an LU-decomposition of M
939 """
941 def __init__(self, M):
942 self.M_lu = lu_factor(M)
943 self.shape = M.shape
944 self.dtype = M.dtype
946 def _matvec(self, x):
947 return lu_solve(self.M_lu, x)
950def gmres_loose(A, b, tol):
951 """
952 gmres with looser termination condition.
953 """
954 b = np.asarray(b)
955 min_tol = 1000 * np.sqrt(b.size) * np.finfo(b.dtype).eps
956 return gmres(A, b, rtol=max(tol, min_tol), atol=0)
959class IterInv(LinearOperator):
960 """
961 IterInv:
962 helper class to repeatedly solve M*x=b
963 using an iterative method.
964 """
966 def __init__(self, M, ifunc=gmres_loose, tol=0):
967 self.M = M
968 if hasattr(M, 'dtype'):
969 self.dtype = M.dtype
970 else:
971 x = np.zeros(M.shape[1])
972 self.dtype = (M * x).dtype
973 self.shape = M.shape
975 if tol <= 0:
976 # when tol=0, ARPACK uses machine tolerance as calculated
977 # by LAPACK's _LAMCH function. We should match this
978 tol = 2 * np.finfo(self.dtype).eps
979 self.ifunc = ifunc
980 self.tol = tol
982 def _matvec(self, x):
983 b, info = self.ifunc(self.M, x, tol=self.tol)
984 if info != 0:
985 raise ValueError("Error in inverting M: function "
986 "%s did not converge (info = %i)."
987 % (self.ifunc.__name__, info))
988 return b
991class IterOpInv(LinearOperator):
992 """
993 IterOpInv:
994 helper class to repeatedly solve [A-sigma*M]*x = b
995 using an iterative method
996 """
998 def __init__(self, A, M, sigma, ifunc=gmres_loose, tol=0):
999 self.A = A
1000 self.M = M
1001 self.sigma = sigma
1003 def mult_func(x):
1004 return A.matvec(x) - sigma * M.matvec(x)
1006 def mult_func_M_None(x):
1007 return A.matvec(x) - sigma * x
1009 x = np.zeros(A.shape[1])
1010 if M is None:
1011 dtype = mult_func_M_None(x).dtype
1012 self.OP = LinearOperator(self.A.shape,
1013 mult_func_M_None,
1014 dtype=dtype)
1015 else:
1016 dtype = mult_func(x).dtype
1017 self.OP = LinearOperator(self.A.shape,
1018 mult_func,
1019 dtype=dtype)
1020 self.shape = A.shape
1022 if tol <= 0:
1023 # when tol=0, ARPACK uses machine tolerance as calculated
1024 # by LAPACK's _LAMCH function. We should match this
1025 tol = 2 * np.finfo(self.OP.dtype).eps
1026 self.ifunc = ifunc
1027 self.tol = tol
1029 def _matvec(self, x):
1030 b, info = self.ifunc(self.OP, x, tol=self.tol)
1031 if info != 0:
1032 raise ValueError("Error in inverting [A-sigma*M]: function "
1033 "%s did not converge (info = %i)."
1034 % (self.ifunc.__name__, info))
1035 return b
1037 @property
1038 def dtype(self):
1039 return self.OP.dtype
1042def _fast_spmatrix_to_csc(A, hermitian=False):
1043 """Convert sparse matrix to CSC (by transposing, if possible)"""
1044 if (A.format == "csr" and hermitian
1045 and not np.issubdtype(A.dtype, np.complexfloating)):
1046 return A.T
1047 elif is_pydata_spmatrix(A):
1048 # No need to convert
1049 return A
1050 else:
1051 return A.tocsc()
1054def get_inv_matvec(M, hermitian=False, tol=0):
1055 if isdense(M):
1056 return LuInv(M).matvec
1057 elif issparse(M) or is_pydata_spmatrix(M):
1058 M = _fast_spmatrix_to_csc(M, hermitian=hermitian)
1059 return SpLuInv(M).matvec
1060 else:
1061 return IterInv(M, tol=tol).matvec
1064def get_OPinv_matvec(A, M, sigma, hermitian=False, tol=0):
1065 if sigma == 0:
1066 return get_inv_matvec(A, hermitian=hermitian, tol=tol)
1068 if M is None:
1069 #M is the identity matrix
1070 if isdense(A):
1071 if (np.issubdtype(A.dtype, np.complexfloating)
1072 or np.imag(sigma) == 0):
1073 A = np.copy(A)
1074 else:
1075 A = A + 0j
1076 A.flat[::A.shape[1] + 1] -= sigma
1077 return LuInv(A).matvec
1078 elif issparse(A) or is_pydata_spmatrix(A):
1079 A = A - sigma * eye(A.shape[0])
1080 A = _fast_spmatrix_to_csc(A, hermitian=hermitian)
1081 return SpLuInv(A).matvec
1082 else:
1083 return IterOpInv(_aslinearoperator_with_dtype(A),
1084 M, sigma, tol=tol).matvec
1085 else:
1086 if ((not isdense(A) and not issparse(A) and not is_pydata_spmatrix(A)) or
1087 (not isdense(M) and not issparse(M) and not is_pydata_spmatrix(A))):
1088 return IterOpInv(_aslinearoperator_with_dtype(A),
1089 _aslinearoperator_with_dtype(M),
1090 sigma, tol=tol).matvec
1091 elif isdense(A) or isdense(M):
1092 return LuInv(A - sigma * M).matvec
1093 else:
1094 OP = A - sigma * M
1095 OP = _fast_spmatrix_to_csc(OP, hermitian=hermitian)
1096 return SpLuInv(OP).matvec
1099# ARPACK is not threadsafe or reentrant (SAVE variables), so we need a
1100# lock and a re-entering check.
1101_ARPACK_LOCK = ReentrancyLock("Nested calls to eigs/eighs not allowed: "
1102 "ARPACK is not re-entrant")
1105def eigs(A, k=6, M=None, sigma=None, which='LM', v0=None,
1106 ncv=None, maxiter=None, tol=0, return_eigenvectors=True,
1107 Minv=None, OPinv=None, OPpart=None):
1108 """
1109 Find k eigenvalues and eigenvectors of the square matrix A.
1111 Solves ``A @ x[i] = w[i] * x[i]``, the standard eigenvalue problem
1112 for w[i] eigenvalues with corresponding eigenvectors x[i].
1114 If M is specified, solves ``A @ x[i] = w[i] * M @ x[i]``, the
1115 generalized eigenvalue problem for w[i] eigenvalues
1116 with corresponding eigenvectors x[i]
1118 Parameters
1119 ----------
1120 A : ndarray, sparse matrix or LinearOperator
1121 An array, sparse matrix, or LinearOperator representing
1122 the operation ``A @ x``, where A is a real or complex square matrix.
1123 k : int, optional
1124 The number of eigenvalues and eigenvectors desired.
1125 `k` must be smaller than N-1. It is not possible to compute all
1126 eigenvectors of a matrix.
1127 M : ndarray, sparse matrix or LinearOperator, optional
1128 An array, sparse matrix, or LinearOperator representing
1129 the operation M@x for the generalized eigenvalue problem
1131 A @ x = w * M @ x.
1133 M must represent a real symmetric matrix if A is real, and must
1134 represent a complex Hermitian matrix if A is complex. For best
1135 results, the data type of M should be the same as that of A.
1136 Additionally:
1138 If `sigma` is None, M is positive definite
1140 If sigma is specified, M is positive semi-definite
1142 If sigma is None, eigs requires an operator to compute the solution
1143 of the linear equation ``M @ x = b``. This is done internally via a
1144 (sparse) LU decomposition for an explicit matrix M, or via an
1145 iterative solver for a general linear operator. Alternatively,
1146 the user can supply the matrix or operator Minv, which gives
1147 ``x = Minv @ b = M^-1 @ b``.
1148 sigma : real or complex, optional
1149 Find eigenvalues near sigma using shift-invert mode. This requires
1150 an operator to compute the solution of the linear system
1151 ``[A - sigma * M] @ x = b``, where M is the identity matrix if
1152 unspecified. This is computed internally via a (sparse) LU
1153 decomposition for explicit matrices A & M, or via an iterative
1154 solver if either A or M is a general linear operator.
1155 Alternatively, the user can supply the matrix or operator OPinv,
1156 which gives ``x = OPinv @ b = [A - sigma * M]^-1 @ b``.
1157 For a real matrix A, shift-invert can either be done in imaginary
1158 mode or real mode, specified by the parameter OPpart ('r' or 'i').
1159 Note that when sigma is specified, the keyword 'which' (below)
1160 refers to the shifted eigenvalues ``w'[i]`` where:
1162 If A is real and OPpart == 'r' (default),
1163 ``w'[i] = 1/2 * [1/(w[i]-sigma) + 1/(w[i]-conj(sigma))]``.
1165 If A is real and OPpart == 'i',
1166 ``w'[i] = 1/2i * [1/(w[i]-sigma) - 1/(w[i]-conj(sigma))]``.
1168 If A is complex, ``w'[i] = 1/(w[i]-sigma)``.
1170 v0 : ndarray, optional
1171 Starting vector for iteration.
1172 Default: random
1173 ncv : int, optional
1174 The number of Lanczos vectors generated
1175 `ncv` must be greater than `k`; it is recommended that ``ncv > 2*k``.
1176 Default: ``min(n, max(2*k + 1, 20))``
1177 which : str, ['LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'], optional
1178 Which `k` eigenvectors and eigenvalues to find:
1180 'LM' : largest magnitude
1182 'SM' : smallest magnitude
1184 'LR' : largest real part
1186 'SR' : smallest real part
1188 'LI' : largest imaginary part
1190 'SI' : smallest imaginary part
1192 When sigma != None, 'which' refers to the shifted eigenvalues w'[i]
1193 (see discussion in 'sigma', above). ARPACK is generally better
1194 at finding large values than small values. If small eigenvalues are
1195 desired, consider using shift-invert mode for better performance.
1196 maxiter : int, optional
1197 Maximum number of Arnoldi update iterations allowed
1198 Default: ``n*10``
1199 tol : float, optional
1200 Relative accuracy for eigenvalues (stopping criterion)
1201 The default value of 0 implies machine precision.
1202 return_eigenvectors : bool, optional
1203 Return eigenvectors (True) in addition to eigenvalues
1204 Minv : ndarray, sparse matrix or LinearOperator, optional
1205 See notes in M, above.
1206 OPinv : ndarray, sparse matrix or LinearOperator, optional
1207 See notes in sigma, above.
1208 OPpart : {'r' or 'i'}, optional
1209 See notes in sigma, above
1211 Returns
1212 -------
1213 w : ndarray
1214 Array of k eigenvalues.
1215 v : ndarray
1216 An array of `k` eigenvectors.
1217 ``v[:, i]`` is the eigenvector corresponding to the eigenvalue w[i].
1219 Raises
1220 ------
1221 ArpackNoConvergence
1222 When the requested convergence is not obtained.
1223 The currently converged eigenvalues and eigenvectors can be found
1224 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
1225 object.
1227 See Also
1228 --------
1229 eigsh : eigenvalues and eigenvectors for symmetric matrix A
1230 svds : singular value decomposition for a matrix A
1232 Notes
1233 -----
1234 This function is a wrapper to the ARPACK [1]_ SNEUPD, DNEUPD, CNEUPD,
1235 ZNEUPD, functions which use the Implicitly Restarted Arnoldi Method to
1236 find the eigenvalues and eigenvectors [2]_.
1238 References
1239 ----------
1240 .. [1] ARPACK Software, https://github.com/opencollab/arpack-ng
1241 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE:
1242 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted
1243 Arnoldi Methods. SIAM, Philadelphia, PA, 1998.
1245 Examples
1246 --------
1247 Find 6 eigenvectors of the identity matrix:
1249 >>> import numpy as np
1250 >>> from scipy.sparse.linalg import eigs
1251 >>> id = np.eye(13)
1252 >>> vals, vecs = eigs(id, k=6)
1253 >>> vals
1254 array([ 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])
1255 >>> vecs.shape
1256 (13, 6)
1258 """
1259 if A.shape[0] != A.shape[1]:
1260 raise ValueError(f'expected square matrix (shape={A.shape})')
1261 if M is not None:
1262 if M.shape != A.shape:
1263 raise ValueError(f'wrong M dimensions {M.shape}, should be {A.shape}')
1264 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower():
1265 warnings.warn('M does not have the same type precision as A. '
1266 'This may adversely affect ARPACK convergence',
1267 stacklevel=2)
1269 n = A.shape[0]
1271 if k <= 0:
1272 raise ValueError("k=%d must be greater than 0." % k)
1274 if k >= n - 1:
1275 warnings.warn("k >= N - 1 for N * N square matrix. "
1276 "Attempting to use scipy.linalg.eig instead.",
1277 RuntimeWarning, stacklevel=2)
1279 if issparse(A):
1280 raise TypeError("Cannot use scipy.linalg.eig for sparse A with "
1281 "k >= N - 1. Use scipy.linalg.eig(A.toarray()) or"
1282 " reduce k.")
1283 if isinstance(A, LinearOperator):
1284 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator "
1285 "A with k >= N - 1.")
1286 if isinstance(M, LinearOperator):
1287 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator "
1288 "M with k >= N - 1.")
1290 return eig(A, b=M, right=return_eigenvectors)
1292 if sigma is None:
1293 matvec = _aslinearoperator_with_dtype(A).matvec
1295 if OPinv is not None:
1296 raise ValueError("OPinv should not be specified "
1297 "with sigma = None.")
1298 if OPpart is not None:
1299 raise ValueError("OPpart should not be specified with "
1300 "sigma = None or complex A")
1302 if M is None:
1303 #standard eigenvalue problem
1304 mode = 1
1305 M_matvec = None
1306 Minv_matvec = None
1307 if Minv is not None:
1308 raise ValueError("Minv should not be "
1309 "specified with M = None.")
1310 else:
1311 #general eigenvalue problem
1312 mode = 2
1313 if Minv is None:
1314 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol)
1315 else:
1316 Minv = _aslinearoperator_with_dtype(Minv)
1317 Minv_matvec = Minv.matvec
1318 M_matvec = _aslinearoperator_with_dtype(M).matvec
1319 else:
1320 #sigma is not None: shift-invert mode
1321 if np.issubdtype(A.dtype, np.complexfloating):
1322 if OPpart is not None:
1323 raise ValueError("OPpart should not be specified "
1324 "with sigma=None or complex A")
1325 mode = 3
1326 elif OPpart is None or OPpart.lower() == 'r':
1327 mode = 3
1328 elif OPpart.lower() == 'i':
1329 if np.imag(sigma) == 0:
1330 raise ValueError("OPpart cannot be 'i' if sigma is real")
1331 mode = 4
1332 else:
1333 raise ValueError("OPpart must be one of ('r','i')")
1335 matvec = _aslinearoperator_with_dtype(A).matvec
1336 if Minv is not None:
1337 raise ValueError("Minv should not be specified when sigma is")
1338 if OPinv is None:
1339 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1340 hermitian=False, tol=tol)
1341 else:
1342 OPinv = _aslinearoperator_with_dtype(OPinv)
1343 Minv_matvec = OPinv.matvec
1344 if M is None:
1345 M_matvec = None
1346 else:
1347 M_matvec = _aslinearoperator_with_dtype(M).matvec
1349 params = _UnsymmetricArpackParams(n, k, A.dtype.char, matvec, mode,
1350 M_matvec, Minv_matvec, sigma,
1351 ncv, v0, maxiter, which, tol)
1353 with _ARPACK_LOCK:
1354 while not params.converged:
1355 params.iterate()
1357 return params.extract(return_eigenvectors)
1360def eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None,
1361 ncv=None, maxiter=None, tol=0, return_eigenvectors=True,
1362 Minv=None, OPinv=None, mode='normal'):
1363 """
1364 Find k eigenvalues and eigenvectors of the real symmetric square matrix
1365 or complex Hermitian matrix A.
1367 Solves ``A @ x[i] = w[i] * x[i]``, the standard eigenvalue problem for
1368 w[i] eigenvalues with corresponding eigenvectors x[i].
1370 If M is specified, solves ``A @ x[i] = w[i] * M @ x[i]``, the
1371 generalized eigenvalue problem for w[i] eigenvalues
1372 with corresponding eigenvectors x[i].
1374 Note that there is no specialized routine for the case when A is a complex
1375 Hermitian matrix. In this case, ``eigsh()`` will call ``eigs()`` and return the
1376 real parts of the eigenvalues thus obtained.
1378 Parameters
1379 ----------
1380 A : ndarray, sparse matrix or LinearOperator
1381 A square operator representing the operation ``A @ x``, where ``A`` is
1382 real symmetric or complex Hermitian. For buckling mode (see below)
1383 ``A`` must additionally be positive-definite.
1384 k : int, optional
1385 The number of eigenvalues and eigenvectors desired.
1386 `k` must be smaller than N. It is not possible to compute all
1387 eigenvectors of a matrix.
1389 Returns
1390 -------
1391 w : array
1392 Array of k eigenvalues.
1393 v : array
1394 An array representing the `k` eigenvectors. The column ``v[:, i]`` is
1395 the eigenvector corresponding to the eigenvalue ``w[i]``.
1397 Other Parameters
1398 ----------------
1399 M : An N x N matrix, array, sparse matrix, or linear operator representing
1400 the operation ``M @ x`` for the generalized eigenvalue problem
1402 A @ x = w * M @ x.
1404 M must represent a real symmetric matrix if A is real, and must
1405 represent a complex Hermitian matrix if A is complex. For best
1406 results, the data type of M should be the same as that of A.
1407 Additionally:
1409 If sigma is None, M is symmetric positive definite.
1411 If sigma is specified, M is symmetric positive semi-definite.
1413 In buckling mode, M is symmetric indefinite.
1415 If sigma is None, eigsh requires an operator to compute the solution
1416 of the linear equation ``M @ x = b``. This is done internally via a
1417 (sparse) LU decomposition for an explicit matrix M, or via an
1418 iterative solver for a general linear operator. Alternatively,
1419 the user can supply the matrix or operator Minv, which gives
1420 ``x = Minv @ b = M^-1 @ b``.
1421 sigma : real
1422 Find eigenvalues near sigma using shift-invert mode. This requires
1423 an operator to compute the solution of the linear system
1424 ``[A - sigma * M] x = b``, where M is the identity matrix if
1425 unspecified. This is computed internally via a (sparse) LU
1426 decomposition for explicit matrices A & M, or via an iterative
1427 solver if either A or M is a general linear operator.
1428 Alternatively, the user can supply the matrix or operator OPinv,
1429 which gives ``x = OPinv @ b = [A - sigma * M]^-1 @ b``.
1430 Note that when sigma is specified, the keyword 'which' refers to
1431 the shifted eigenvalues ``w'[i]`` where:
1433 if mode == 'normal', ``w'[i] = 1 / (w[i] - sigma)``.
1435 if mode == 'cayley', ``w'[i] = (w[i] + sigma) / (w[i] - sigma)``.
1437 if mode == 'buckling', ``w'[i] = w[i] / (w[i] - sigma)``.
1439 (see further discussion in 'mode' below)
1440 v0 : ndarray, optional
1441 Starting vector for iteration.
1442 Default: random
1443 ncv : int, optional
1444 The number of Lanczos vectors generated ncv must be greater than k and
1445 smaller than n; it is recommended that ``ncv > 2*k``.
1446 Default: ``min(n, max(2*k + 1, 20))``
1447 which : str ['LM' | 'SM' | 'LA' | 'SA' | 'BE']
1448 If A is a complex Hermitian matrix, 'BE' is invalid.
1449 Which `k` eigenvectors and eigenvalues to find:
1451 'LM' : Largest (in magnitude) eigenvalues.
1453 'SM' : Smallest (in magnitude) eigenvalues.
1455 'LA' : Largest (algebraic) eigenvalues.
1457 'SA' : Smallest (algebraic) eigenvalues.
1459 'BE' : Half (k/2) from each end of the spectrum.
1461 When k is odd, return one more (k/2+1) from the high end.
1462 When sigma != None, 'which' refers to the shifted eigenvalues ``w'[i]``
1463 (see discussion in 'sigma', above). ARPACK is generally better
1464 at finding large values than small values. If small eigenvalues are
1465 desired, consider using shift-invert mode for better performance.
1466 maxiter : int, optional
1467 Maximum number of Arnoldi update iterations allowed.
1468 Default: ``n*10``
1469 tol : float
1470 Relative accuracy for eigenvalues (stopping criterion).
1471 The default value of 0 implies machine precision.
1472 Minv : N x N matrix, array, sparse matrix, or LinearOperator
1473 See notes in M, above.
1474 OPinv : N x N matrix, array, sparse matrix, or LinearOperator
1475 See notes in sigma, above.
1476 return_eigenvectors : bool
1477 Return eigenvectors (True) in addition to eigenvalues.
1478 This value determines the order in which eigenvalues are sorted.
1479 The sort order is also dependent on the `which` variable.
1481 For which = 'LM' or 'SA':
1482 If `return_eigenvectors` is True, eigenvalues are sorted by
1483 algebraic value.
1485 If `return_eigenvectors` is False, eigenvalues are sorted by
1486 absolute value.
1488 For which = 'BE' or 'LA':
1489 eigenvalues are always sorted by algebraic value.
1491 For which = 'SM':
1492 If `return_eigenvectors` is True, eigenvalues are sorted by
1493 algebraic value.
1495 If `return_eigenvectors` is False, eigenvalues are sorted by
1496 decreasing absolute value.
1498 mode : string ['normal' | 'buckling' | 'cayley']
1499 Specify strategy to use for shift-invert mode. This argument applies
1500 only for real-valued A and sigma != None. For shift-invert mode,
1501 ARPACK internally solves the eigenvalue problem
1502 ``OP @ x'[i] = w'[i] * B @ x'[i]``
1503 and transforms the resulting Ritz vectors x'[i] and Ritz values w'[i]
1504 into the desired eigenvectors and eigenvalues of the problem
1505 ``A @ x[i] = w[i] * M @ x[i]``.
1506 The modes are as follows:
1508 'normal' :
1509 OP = [A - sigma * M]^-1 @ M,
1510 B = M,
1511 w'[i] = 1 / (w[i] - sigma)
1513 'buckling' :
1514 OP = [A - sigma * M]^-1 @ A,
1515 B = A,
1516 w'[i] = w[i] / (w[i] - sigma)
1518 'cayley' :
1519 OP = [A - sigma * M]^-1 @ [A + sigma * M],
1520 B = M,
1521 w'[i] = (w[i] + sigma) / (w[i] - sigma)
1523 The choice of mode will affect which eigenvalues are selected by
1524 the keyword 'which', and can also impact the stability of
1525 convergence (see [2] for a discussion).
1527 Raises
1528 ------
1529 ArpackNoConvergence
1530 When the requested convergence is not obtained.
1532 The currently converged eigenvalues and eigenvectors can be found
1533 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
1534 object.
1536 See Also
1537 --------
1538 eigs : eigenvalues and eigenvectors for a general (nonsymmetric) matrix A
1539 svds : singular value decomposition for a matrix A
1541 Notes
1542 -----
1543 This function is a wrapper to the ARPACK [1]_ SSEUPD and DSEUPD
1544 functions which use the Implicitly Restarted Lanczos Method to
1545 find the eigenvalues and eigenvectors [2]_.
1547 References
1548 ----------
1549 .. [1] ARPACK Software, https://github.com/opencollab/arpack-ng
1550 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE:
1551 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted
1552 Arnoldi Methods. SIAM, Philadelphia, PA, 1998.
1554 Examples
1555 --------
1556 >>> import numpy as np
1557 >>> from scipy.sparse.linalg import eigsh
1558 >>> identity = np.eye(13)
1559 >>> eigenvalues, eigenvectors = eigsh(identity, k=6)
1560 >>> eigenvalues
1561 array([1., 1., 1., 1., 1., 1.])
1562 >>> eigenvectors.shape
1563 (13, 6)
1565 """
1566 # complex Hermitian matrices should be solved with eigs
1567 if np.issubdtype(A.dtype, np.complexfloating):
1568 if mode != 'normal':
1569 raise ValueError("mode=%s cannot be used with "
1570 "complex matrix A" % mode)
1571 if which == 'BE':
1572 raise ValueError("which='BE' cannot be used with complex matrix A")
1573 elif which == 'LA':
1574 which = 'LR'
1575 elif which == 'SA':
1576 which = 'SR'
1577 ret = eigs(A, k, M=M, sigma=sigma, which=which, v0=v0,
1578 ncv=ncv, maxiter=maxiter, tol=tol,
1579 return_eigenvectors=return_eigenvectors, Minv=Minv,
1580 OPinv=OPinv)
1582 if return_eigenvectors:
1583 return ret[0].real, ret[1]
1584 else:
1585 return ret.real
1587 if A.shape[0] != A.shape[1]:
1588 raise ValueError(f'expected square matrix (shape={A.shape})')
1589 if M is not None:
1590 if M.shape != A.shape:
1591 raise ValueError(f'wrong M dimensions {M.shape}, should be {A.shape}')
1592 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower():
1593 warnings.warn('M does not have the same type precision as A. '
1594 'This may adversely affect ARPACK convergence',
1595 stacklevel=2)
1597 n = A.shape[0]
1599 if k <= 0:
1600 raise ValueError("k must be greater than 0.")
1602 if k >= n:
1603 warnings.warn("k >= N for N * N square matrix. "
1604 "Attempting to use scipy.linalg.eigh instead.",
1605 RuntimeWarning, stacklevel=2)
1607 if issparse(A):
1608 raise TypeError("Cannot use scipy.linalg.eigh for sparse A with "
1609 "k >= N. Use scipy.linalg.eigh(A.toarray()) or"
1610 " reduce k.")
1611 if isinstance(A, LinearOperator):
1612 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator "
1613 "A with k >= N.")
1614 if isinstance(M, LinearOperator):
1615 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator "
1616 "M with k >= N.")
1618 return eigh(A, b=M, eigvals_only=not return_eigenvectors)
1620 if sigma is None:
1621 A = _aslinearoperator_with_dtype(A)
1622 matvec = A.matvec
1624 if OPinv is not None:
1625 raise ValueError("OPinv should not be specified "
1626 "with sigma = None.")
1627 if M is None:
1628 #standard eigenvalue problem
1629 mode = 1
1630 M_matvec = None
1631 Minv_matvec = None
1632 if Minv is not None:
1633 raise ValueError("Minv should not be "
1634 "specified with M = None.")
1635 else:
1636 #general eigenvalue problem
1637 mode = 2
1638 if Minv is None:
1639 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol)
1640 else:
1641 Minv = _aslinearoperator_with_dtype(Minv)
1642 Minv_matvec = Minv.matvec
1643 M_matvec = _aslinearoperator_with_dtype(M).matvec
1644 else:
1645 # sigma is not None: shift-invert mode
1646 if Minv is not None:
1647 raise ValueError("Minv should not be specified when sigma is")
1649 # normal mode
1650 if mode == 'normal':
1651 mode = 3
1652 matvec = None
1653 if OPinv is None:
1654 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1655 hermitian=True, tol=tol)
1656 else:
1657 OPinv = _aslinearoperator_with_dtype(OPinv)
1658 Minv_matvec = OPinv.matvec
1659 if M is None:
1660 M_matvec = None
1661 else:
1662 M = _aslinearoperator_with_dtype(M)
1663 M_matvec = M.matvec
1665 # buckling mode
1666 elif mode == 'buckling':
1667 mode = 4
1668 if OPinv is None:
1669 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1670 hermitian=True, tol=tol)
1671 else:
1672 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec
1673 matvec = _aslinearoperator_with_dtype(A).matvec
1674 M_matvec = None
1676 # cayley-transform mode
1677 elif mode == 'cayley':
1678 mode = 5
1679 matvec = _aslinearoperator_with_dtype(A).matvec
1680 if OPinv is None:
1681 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1682 hermitian=True, tol=tol)
1683 else:
1684 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec
1685 if M is None:
1686 M_matvec = None
1687 else:
1688 M_matvec = _aslinearoperator_with_dtype(M).matvec
1690 # unrecognized mode
1691 else:
1692 raise ValueError("unrecognized mode '%s'" % mode)
1694 params = _SymmetricArpackParams(n, k, A.dtype.char, matvec, mode,
1695 M_matvec, Minv_matvec, sigma,
1696 ncv, v0, maxiter, which, tol)
1698 with _ARPACK_LOCK:
1699 while not params.converged:
1700 params.iterate()
1702 return params.extract(return_eigenvectors)