Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py: 12%
637 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2Find a few eigenvectors and eigenvalues of a matrix.
5Uses ARPACK: http://www.caam.rice.edu/software/ARPACK/
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.
38__docformat__ = "restructuredtext en"
40__all__ = ['eigs', 'eigsh', 'ArpackError', 'ArpackNoConvergence']
42from . import _arpack
43arpack_int = _arpack.timing.nbx.dtype
45import numpy as np
46import warnings
47from scipy.sparse.linalg._interface import aslinearoperator, LinearOperator
48from scipy.sparse import eye, issparse, isspmatrix, isspmatrix_csr
49from scipy.linalg import eig, eigh, lu_factor, lu_solve
50from scipy.sparse._sputils import isdense, is_pydata_spmatrix
51from scipy.sparse.linalg import gmres, splu
52from scipy._lib._util import _aligned_zeros
53from scipy._lib._threadsafety import ReentrancyLock
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 = "%s [%s]" % (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 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info =\
724 self._arpack_solver(self.ido, self.bmat, self.which, self.k,
725 self.tol, self.resid, self.v, self.iparam,
726 self.ipntr, self.workd, self.workl,
727 self.info)
728 else:
729 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info =\
730 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)
735 xslice = slice(self.ipntr[0] - 1, self.ipntr[0] - 1 + self.n)
736 yslice = slice(self.ipntr[1] - 1, self.ipntr[1] - 1 + self.n)
737 if self.ido == -1:
738 # initialization
739 self.workd[yslice] = self.OP(self.workd[xslice])
740 elif self.ido == 1:
741 # compute y = Op*x
742 if self.mode in (1, 2):
743 self.workd[yslice] = self.OP(self.workd[xslice])
744 else:
745 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n)
746 self.workd[yslice] = self.OPa(self.workd[Bxslice])
747 elif self.ido == 2:
748 self.workd[yslice] = self.B(self.workd[xslice])
749 elif self.ido == 3:
750 raise ValueError("ARPACK requested user shifts. Assure ISHIFT==0")
751 else:
752 self.converged = True
754 if self.info == 0:
755 pass
756 elif self.info == 1:
757 self._raise_no_convergence()
758 else:
759 raise ArpackError(self.info, infodict=self.iterate_infodict)
761 def extract(self, return_eigenvectors):
762 k, n = self.k, self.n
764 ierr = 0
765 howmny = 'A' # return all eigenvectors
766 sselect = np.zeros(self.ncv, 'int') # unused
767 sigmar = np.real(self.sigma)
768 sigmai = np.imag(self.sigma)
769 workev = np.zeros(3 * self.ncv, self.tp)
771 if self.tp in 'fd':
772 dr = np.zeros(k + 1, self.tp)
773 di = np.zeros(k + 1, self.tp)
774 zr = np.zeros((n, k + 1), self.tp)
775 dr, di, zr, ierr = \
776 self._arpack_extract(return_eigenvectors,
777 howmny, sselect, sigmar, sigmai, workev,
778 self.bmat, self.which, k, self.tol, self.resid,
779 self.v, self.iparam, self.ipntr,
780 self.workd, self.workl, self.info)
781 if ierr != 0:
782 raise ArpackError(ierr, infodict=self.extract_infodict)
783 nreturned = self.iparam[4] # number of good eigenvalues returned
785 # Build complex eigenvalues from real and imaginary parts
786 d = dr + 1.0j * di
788 # Arrange the eigenvectors: complex eigenvectors are stored as
789 # real,imaginary in consecutive columns
790 z = zr.astype(self.tp.upper())
792 # The ARPACK nonsymmetric real and double interface (s,d)naupd
793 # return eigenvalues and eigenvectors in real (float,double)
794 # arrays.
796 # Efficiency: this should check that return_eigenvectors == True
797 # before going through this construction.
798 if sigmai == 0:
799 i = 0
800 while i <= k:
801 # check if complex
802 if abs(d[i].imag) != 0:
803 # this is a complex conjugate pair with eigenvalues
804 # in consecutive columns
805 if i < k:
806 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1]
807 z[:, i + 1] = z[:, i].conjugate()
808 i += 1
809 else:
810 #last eigenvalue is complex: the imaginary part of
811 # the eigenvector has not been returned
812 #this can only happen if nreturned > k, so we'll
813 # throw out this case.
814 nreturned -= 1
815 i += 1
817 else:
818 # real matrix, mode 3 or 4, imag(sigma) is nonzero:
819 # see remark 3 in <s,d>neupd.f
820 # Build complex eigenvalues from real and imaginary parts
821 i = 0
822 while i <= k:
823 if abs(d[i].imag) == 0:
824 d[i] = np.dot(zr[:, i], self.matvec(zr[:, i]))
825 else:
826 if i < k:
827 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1]
828 z[:, i + 1] = z[:, i].conjugate()
829 d[i] = ((np.dot(zr[:, i],
830 self.matvec(zr[:, i]))
831 + np.dot(zr[:, i + 1],
832 self.matvec(zr[:, i + 1])))
833 + 1j * (np.dot(zr[:, i],
834 self.matvec(zr[:, i + 1]))
835 - np.dot(zr[:, i + 1],
836 self.matvec(zr[:, i]))))
837 d[i + 1] = d[i].conj()
838 i += 1
839 else:
840 #last eigenvalue is complex: the imaginary part of
841 # the eigenvector has not been returned
842 #this can only happen if nreturned > k, so we'll
843 # throw out this case.
844 nreturned -= 1
845 i += 1
847 # Now we have k+1 possible eigenvalues and eigenvectors
848 # Return the ones specified by the keyword "which"
850 if nreturned <= k:
851 # we got less or equal as many eigenvalues we wanted
852 d = d[:nreturned]
853 z = z[:, :nreturned]
854 else:
855 # we got one extra eigenvalue (likely a cc pair, but which?)
856 if self.mode in (1, 2):
857 rd = d
858 elif self.mode in (3, 4):
859 rd = 1 / (d - self.sigma)
861 if self.which in ['LR', 'SR']:
862 ind = np.argsort(rd.real)
863 elif self.which in ['LI', 'SI']:
864 # for LI,SI ARPACK returns largest,smallest
865 # abs(imaginary) (complex pairs come together)
866 ind = np.argsort(abs(rd.imag))
867 else:
868 ind = np.argsort(abs(rd))
870 if self.which in ['LR', 'LM', 'LI']:
871 ind = ind[-k:][::-1]
872 elif self.which in ['SR', 'SM', 'SI']:
873 ind = ind[:k]
875 d = d[ind]
876 z = z[:, ind]
877 else:
878 # complex is so much simpler...
879 d, z, ierr =\
880 self._arpack_extract(return_eigenvectors,
881 howmny, sselect, self.sigma, workev,
882 self.bmat, self.which, k, self.tol, self.resid,
883 self.v, self.iparam, self.ipntr,
884 self.workd, self.workl, self.rwork, ierr)
886 if ierr != 0:
887 raise ArpackError(ierr, infodict=self.extract_infodict)
889 k_ok = self.iparam[4]
890 d = d[:k_ok]
891 z = z[:, :k_ok]
893 if return_eigenvectors:
894 return d, z
895 else:
896 return d
899def _aslinearoperator_with_dtype(m):
900 m = aslinearoperator(m)
901 if not hasattr(m, 'dtype'):
902 x = np.zeros(m.shape[1])
903 m.dtype = (m * x).dtype
904 return m
907class SpLuInv(LinearOperator):
908 """
909 SpLuInv:
910 helper class to repeatedly solve M*x=b
911 using a sparse LU-decomposition of M
912 """
914 def __init__(self, M):
915 self.M_lu = splu(M)
916 self.shape = M.shape
917 self.dtype = M.dtype
918 self.isreal = not np.issubdtype(self.dtype, np.complexfloating)
920 def _matvec(self, x):
921 # careful here: splu.solve will throw away imaginary
922 # part of x if M is real
923 x = np.asarray(x)
924 if self.isreal and np.issubdtype(x.dtype, np.complexfloating):
925 return (self.M_lu.solve(np.real(x).astype(self.dtype))
926 + 1j * self.M_lu.solve(np.imag(x).astype(self.dtype)))
927 else:
928 return self.M_lu.solve(x.astype(self.dtype))
931class LuInv(LinearOperator):
932 """
933 LuInv:
934 helper class to repeatedly solve M*x=b
935 using an LU-decomposition of M
936 """
938 def __init__(self, M):
939 self.M_lu = lu_factor(M)
940 self.shape = M.shape
941 self.dtype = M.dtype
943 def _matvec(self, x):
944 return lu_solve(self.M_lu, x)
947def gmres_loose(A, b, tol):
948 """
949 gmres with looser termination condition.
950 """
951 b = np.asarray(b)
952 min_tol = 1000 * np.sqrt(b.size) * np.finfo(b.dtype).eps
953 return gmres(A, b, tol=max(tol, min_tol), atol=0)
956class IterInv(LinearOperator):
957 """
958 IterInv:
959 helper class to repeatedly solve M*x=b
960 using an iterative method.
961 """
963 def __init__(self, M, ifunc=gmres_loose, tol=0):
964 self.M = M
965 if hasattr(M, 'dtype'):
966 self.dtype = M.dtype
967 else:
968 x = np.zeros(M.shape[1])
969 self.dtype = (M * x).dtype
970 self.shape = M.shape
972 if tol <= 0:
973 # when tol=0, ARPACK uses machine tolerance as calculated
974 # by LAPACK's _LAMCH function. We should match this
975 tol = 2 * np.finfo(self.dtype).eps
976 self.ifunc = ifunc
977 self.tol = tol
979 def _matvec(self, x):
980 b, info = self.ifunc(self.M, x, tol=self.tol)
981 if info != 0:
982 raise ValueError("Error in inverting M: function "
983 "%s did not converge (info = %i)."
984 % (self.ifunc.__name__, info))
985 return b
988class IterOpInv(LinearOperator):
989 """
990 IterOpInv:
991 helper class to repeatedly solve [A-sigma*M]*x = b
992 using an iterative method
993 """
995 def __init__(self, A, M, sigma, ifunc=gmres_loose, tol=0):
996 self.A = A
997 self.M = M
998 self.sigma = sigma
1000 def mult_func(x):
1001 return A.matvec(x) - sigma * M.matvec(x)
1003 def mult_func_M_None(x):
1004 return A.matvec(x) - sigma * x
1006 x = np.zeros(A.shape[1])
1007 if M is None:
1008 dtype = mult_func_M_None(x).dtype
1009 self.OP = LinearOperator(self.A.shape,
1010 mult_func_M_None,
1011 dtype=dtype)
1012 else:
1013 dtype = mult_func(x).dtype
1014 self.OP = LinearOperator(self.A.shape,
1015 mult_func,
1016 dtype=dtype)
1017 self.shape = A.shape
1019 if tol <= 0:
1020 # when tol=0, ARPACK uses machine tolerance as calculated
1021 # by LAPACK's _LAMCH function. We should match this
1022 tol = 2 * np.finfo(self.OP.dtype).eps
1023 self.ifunc = ifunc
1024 self.tol = tol
1026 def _matvec(self, x):
1027 b, info = self.ifunc(self.OP, x, tol=self.tol)
1028 if info != 0:
1029 raise ValueError("Error in inverting [A-sigma*M]: function "
1030 "%s did not converge (info = %i)."
1031 % (self.ifunc.__name__, info))
1032 return b
1034 @property
1035 def dtype(self):
1036 return self.OP.dtype
1039def _fast_spmatrix_to_csc(A, hermitian=False):
1040 """Convert sparse matrix to CSC (by transposing, if possible)"""
1041 if (isspmatrix_csr(A) and hermitian
1042 and not np.issubdtype(A.dtype, np.complexfloating)):
1043 return A.T
1044 elif is_pydata_spmatrix(A):
1045 # No need to convert
1046 return A
1047 else:
1048 return A.tocsc()
1051def get_inv_matvec(M, hermitian=False, tol=0):
1052 if isdense(M):
1053 return LuInv(M).matvec
1054 elif isspmatrix(M) or is_pydata_spmatrix(M):
1055 M = _fast_spmatrix_to_csc(M, hermitian=hermitian)
1056 return SpLuInv(M).matvec
1057 else:
1058 return IterInv(M, tol=tol).matvec
1061def get_OPinv_matvec(A, M, sigma, hermitian=False, tol=0):
1062 if sigma == 0:
1063 return get_inv_matvec(A, hermitian=hermitian, tol=tol)
1065 if M is None:
1066 #M is the identity matrix
1067 if isdense(A):
1068 if (np.issubdtype(A.dtype, np.complexfloating)
1069 or np.imag(sigma) == 0):
1070 A = np.copy(A)
1071 else:
1072 A = A + 0j
1073 A.flat[::A.shape[1] + 1] -= sigma
1074 return LuInv(A).matvec
1075 elif isspmatrix(A) or is_pydata_spmatrix(A):
1076 A = A - sigma * eye(A.shape[0])
1077 A = _fast_spmatrix_to_csc(A, hermitian=hermitian)
1078 return SpLuInv(A).matvec
1079 else:
1080 return IterOpInv(_aslinearoperator_with_dtype(A),
1081 M, sigma, tol=tol).matvec
1082 else:
1083 if ((not isdense(A) and not isspmatrix(A) and not is_pydata_spmatrix(A)) or
1084 (not isdense(M) and not isspmatrix(M) and not is_pydata_spmatrix(A))):
1085 return IterOpInv(_aslinearoperator_with_dtype(A),
1086 _aslinearoperator_with_dtype(M),
1087 sigma, tol=tol).matvec
1088 elif isdense(A) or isdense(M):
1089 return LuInv(A - sigma * M).matvec
1090 else:
1091 OP = A - sigma * M
1092 OP = _fast_spmatrix_to_csc(OP, hermitian=hermitian)
1093 return SpLuInv(OP).matvec
1096# ARPACK is not threadsafe or reentrant (SAVE variables), so we need a
1097# lock and a re-entering check.
1098_ARPACK_LOCK = ReentrancyLock("Nested calls to eigs/eighs not allowed: "
1099 "ARPACK is not re-entrant")
1102def eigs(A, k=6, M=None, sigma=None, which='LM', v0=None,
1103 ncv=None, maxiter=None, tol=0, return_eigenvectors=True,
1104 Minv=None, OPinv=None, OPpart=None):
1105 """
1106 Find k eigenvalues and eigenvectors of the square matrix A.
1108 Solves ``A @ x[i] = w[i] * x[i]``, the standard eigenvalue problem
1109 for w[i] eigenvalues with corresponding eigenvectors x[i].
1111 If M is specified, solves ``A @ x[i] = w[i] * M @ x[i]``, the
1112 generalized eigenvalue problem for w[i] eigenvalues
1113 with corresponding eigenvectors x[i]
1115 Parameters
1116 ----------
1117 A : ndarray, sparse matrix or LinearOperator
1118 An array, sparse matrix, or LinearOperator representing
1119 the operation ``A @ x``, where A is a real or complex square matrix.
1120 k : int, optional
1121 The number of eigenvalues and eigenvectors desired.
1122 `k` must be smaller than N-1. It is not possible to compute all
1123 eigenvectors of a matrix.
1124 M : ndarray, sparse matrix or LinearOperator, optional
1125 An array, sparse matrix, or LinearOperator representing
1126 the operation M@x for the generalized eigenvalue problem
1128 A @ x = w * M @ x.
1130 M must represent a real symmetric matrix if A is real, and must
1131 represent a complex Hermitian matrix if A is complex. For best
1132 results, the data type of M should be the same as that of A.
1133 Additionally:
1135 If `sigma` is None, M is positive definite
1137 If sigma is specified, M is positive semi-definite
1139 If sigma is None, eigs requires an operator to compute the solution
1140 of the linear equation ``M @ x = b``. This is done internally via a
1141 (sparse) LU decomposition for an explicit matrix M, or via an
1142 iterative solver for a general linear operator. Alternatively,
1143 the user can supply the matrix or operator Minv, which gives
1144 ``x = Minv @ b = M^-1 @ b``.
1145 sigma : real or complex, optional
1146 Find eigenvalues near sigma using shift-invert mode. This requires
1147 an operator to compute the solution of the linear system
1148 ``[A - sigma * M] @ x = b``, where M is the identity matrix if
1149 unspecified. This is computed internally via a (sparse) LU
1150 decomposition for explicit matrices A & M, or via an iterative
1151 solver if either A or M is a general linear operator.
1152 Alternatively, the user can supply the matrix or operator OPinv,
1153 which gives ``x = OPinv @ b = [A - sigma * M]^-1 @ b``.
1154 For a real matrix A, shift-invert can either be done in imaginary
1155 mode or real mode, specified by the parameter OPpart ('r' or 'i').
1156 Note that when sigma is specified, the keyword 'which' (below)
1157 refers to the shifted eigenvalues ``w'[i]`` where:
1159 If A is real and OPpart == 'r' (default),
1160 ``w'[i] = 1/2 * [1/(w[i]-sigma) + 1/(w[i]-conj(sigma))]``.
1162 If A is real and OPpart == 'i',
1163 ``w'[i] = 1/2i * [1/(w[i]-sigma) - 1/(w[i]-conj(sigma))]``.
1165 If A is complex, ``w'[i] = 1/(w[i]-sigma)``.
1167 v0 : ndarray, optional
1168 Starting vector for iteration.
1169 Default: random
1170 ncv : int, optional
1171 The number of Lanczos vectors generated
1172 `ncv` must be greater than `k`; it is recommended that ``ncv > 2*k``.
1173 Default: ``min(n, max(2*k + 1, 20))``
1174 which : str, ['LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'], optional
1175 Which `k` eigenvectors and eigenvalues to find:
1177 'LM' : largest magnitude
1179 'SM' : smallest magnitude
1181 'LR' : largest real part
1183 'SR' : smallest real part
1185 'LI' : largest imaginary part
1187 'SI' : smallest imaginary part
1189 When sigma != None, 'which' refers to the shifted eigenvalues w'[i]
1190 (see discussion in 'sigma', above). ARPACK is generally better
1191 at finding large values than small values. If small eigenvalues are
1192 desired, consider using shift-invert mode for better performance.
1193 maxiter : int, optional
1194 Maximum number of Arnoldi update iterations allowed
1195 Default: ``n*10``
1196 tol : float, optional
1197 Relative accuracy for eigenvalues (stopping criterion)
1198 The default value of 0 implies machine precision.
1199 return_eigenvectors : bool, optional
1200 Return eigenvectors (True) in addition to eigenvalues
1201 Minv : ndarray, sparse matrix or LinearOperator, optional
1202 See notes in M, above.
1203 OPinv : ndarray, sparse matrix or LinearOperator, optional
1204 See notes in sigma, above.
1205 OPpart : {'r' or 'i'}, optional
1206 See notes in sigma, above
1208 Returns
1209 -------
1210 w : ndarray
1211 Array of k eigenvalues.
1212 v : ndarray
1213 An array of `k` eigenvectors.
1214 ``v[:, i]`` is the eigenvector corresponding to the eigenvalue w[i].
1216 Raises
1217 ------
1218 ArpackNoConvergence
1219 When the requested convergence is not obtained.
1220 The currently converged eigenvalues and eigenvectors can be found
1221 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
1222 object.
1224 See Also
1225 --------
1226 eigsh : eigenvalues and eigenvectors for symmetric matrix A
1227 svds : singular value decomposition for a matrix A
1229 Notes
1230 -----
1231 This function is a wrapper to the ARPACK [1]_ SNEUPD, DNEUPD, CNEUPD,
1232 ZNEUPD, functions which use the Implicitly Restarted Arnoldi Method to
1233 find the eigenvalues and eigenvectors [2]_.
1235 References
1236 ----------
1237 .. [1] ARPACK Software, http://www.caam.rice.edu/software/ARPACK/
1238 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE:
1239 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted
1240 Arnoldi Methods. SIAM, Philadelphia, PA, 1998.
1242 Examples
1243 --------
1244 Find 6 eigenvectors of the identity matrix:
1246 >>> import numpy as np
1247 >>> from scipy.sparse.linalg import eigs
1248 >>> id = np.eye(13)
1249 >>> vals, vecs = eigs(id, k=6)
1250 >>> vals
1251 array([ 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])
1252 >>> vecs.shape
1253 (13, 6)
1255 """
1256 if A.shape[0] != A.shape[1]:
1257 raise ValueError('expected square matrix (shape=%s)' % (A.shape,))
1258 if M is not None:
1259 if M.shape != A.shape:
1260 raise ValueError('wrong M dimensions %s, should be %s'
1261 % (M.shape, A.shape))
1262 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower():
1263 warnings.warn('M does not have the same type precision as A. '
1264 'This may adversely affect ARPACK convergence')
1266 n = A.shape[0]
1268 if k <= 0:
1269 raise ValueError("k=%d must be greater than 0." % k)
1271 if k >= n - 1:
1272 warnings.warn("k >= N - 1 for N * N square matrix. "
1273 "Attempting to use scipy.linalg.eig instead.",
1274 RuntimeWarning)
1276 if issparse(A):
1277 raise TypeError("Cannot use scipy.linalg.eig for sparse A with "
1278 "k >= N - 1. Use scipy.linalg.eig(A.toarray()) or"
1279 " reduce k.")
1280 if isinstance(A, LinearOperator):
1281 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator "
1282 "A with k >= N - 1.")
1283 if isinstance(M, LinearOperator):
1284 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator "
1285 "M with k >= N - 1.")
1287 return eig(A, b=M, right=return_eigenvectors)
1289 if sigma is None:
1290 matvec = _aslinearoperator_with_dtype(A).matvec
1292 if OPinv is not None:
1293 raise ValueError("OPinv should not be specified "
1294 "with sigma = None.")
1295 if OPpart is not None:
1296 raise ValueError("OPpart should not be specified with "
1297 "sigma = None or complex A")
1299 if M is None:
1300 #standard eigenvalue problem
1301 mode = 1
1302 M_matvec = None
1303 Minv_matvec = None
1304 if Minv is not None:
1305 raise ValueError("Minv should not be "
1306 "specified with M = None.")
1307 else:
1308 #general eigenvalue problem
1309 mode = 2
1310 if Minv is None:
1311 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol)
1312 else:
1313 Minv = _aslinearoperator_with_dtype(Minv)
1314 Minv_matvec = Minv.matvec
1315 M_matvec = _aslinearoperator_with_dtype(M).matvec
1316 else:
1317 #sigma is not None: shift-invert mode
1318 if np.issubdtype(A.dtype, np.complexfloating):
1319 if OPpart is not None:
1320 raise ValueError("OPpart should not be specified "
1321 "with sigma=None or complex A")
1322 mode = 3
1323 elif OPpart is None or OPpart.lower() == 'r':
1324 mode = 3
1325 elif OPpart.lower() == 'i':
1326 if np.imag(sigma) == 0:
1327 raise ValueError("OPpart cannot be 'i' if sigma is real")
1328 mode = 4
1329 else:
1330 raise ValueError("OPpart must be one of ('r','i')")
1332 matvec = _aslinearoperator_with_dtype(A).matvec
1333 if Minv is not None:
1334 raise ValueError("Minv should not be specified when sigma is")
1335 if OPinv is None:
1336 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1337 hermitian=False, tol=tol)
1338 else:
1339 OPinv = _aslinearoperator_with_dtype(OPinv)
1340 Minv_matvec = OPinv.matvec
1341 if M is None:
1342 M_matvec = None
1343 else:
1344 M_matvec = _aslinearoperator_with_dtype(M).matvec
1346 params = _UnsymmetricArpackParams(n, k, A.dtype.char, matvec, mode,
1347 M_matvec, Minv_matvec, sigma,
1348 ncv, v0, maxiter, which, tol)
1350 with _ARPACK_LOCK:
1351 while not params.converged:
1352 params.iterate()
1354 return params.extract(return_eigenvectors)
1357def eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None,
1358 ncv=None, maxiter=None, tol=0, return_eigenvectors=True,
1359 Minv=None, OPinv=None, mode='normal'):
1360 """
1361 Find k eigenvalues and eigenvectors of the real symmetric square matrix
1362 or complex Hermitian matrix A.
1364 Solves ``A @ x[i] = w[i] * x[i]``, the standard eigenvalue problem for
1365 w[i] eigenvalues with corresponding eigenvectors x[i].
1367 If M is specified, solves ``A @ x[i] = w[i] * M @ x[i]``, the
1368 generalized eigenvalue problem for w[i] eigenvalues
1369 with corresponding eigenvectors x[i].
1371 Note that there is no specialized routine for the case when A is a complex
1372 Hermitian matrix. In this case, ``eigsh()`` will call ``eigs()`` and return the
1373 real parts of the eigenvalues thus obtained.
1375 Parameters
1376 ----------
1377 A : ndarray, sparse matrix or LinearOperator
1378 A square operator representing the operation ``A @ x``, where ``A`` is
1379 real symmetric or complex Hermitian. For buckling mode (see below)
1380 ``A`` must additionally be positive-definite.
1381 k : int, optional
1382 The number of eigenvalues and eigenvectors desired.
1383 `k` must be smaller than N. It is not possible to compute all
1384 eigenvectors of a matrix.
1386 Returns
1387 -------
1388 w : array
1389 Array of k eigenvalues.
1390 v : array
1391 An array representing the `k` eigenvectors. The column ``v[:, i]`` is
1392 the eigenvector corresponding to the eigenvalue ``w[i]``.
1394 Other Parameters
1395 ----------------
1396 M : An N x N matrix, array, sparse matrix, or linear operator representing
1397 the operation ``M @ x`` for the generalized eigenvalue problem
1399 A @ x = w * M @ x.
1401 M must represent a real symmetric matrix if A is real, and must
1402 represent a complex Hermitian matrix if A is complex. For best
1403 results, the data type of M should be the same as that of A.
1404 Additionally:
1406 If sigma is None, M is symmetric positive definite.
1408 If sigma is specified, M is symmetric positive semi-definite.
1410 In buckling mode, M is symmetric indefinite.
1412 If sigma is None, eigsh requires an operator to compute the solution
1413 of the linear equation ``M @ x = b``. This is done internally via a
1414 (sparse) LU decomposition for an explicit matrix M, or via an
1415 iterative solver for a general linear operator. Alternatively,
1416 the user can supply the matrix or operator Minv, which gives
1417 ``x = Minv @ b = M^-1 @ b``.
1418 sigma : real
1419 Find eigenvalues near sigma using shift-invert mode. This requires
1420 an operator to compute the solution of the linear system
1421 ``[A - sigma * M] x = b``, where M is the identity matrix if
1422 unspecified. This is computed internally via a (sparse) LU
1423 decomposition for explicit matrices A & M, or via an iterative
1424 solver if either A or M is a general linear operator.
1425 Alternatively, the user can supply the matrix or operator OPinv,
1426 which gives ``x = OPinv @ b = [A - sigma * M]^-1 @ b``.
1427 Note that when sigma is specified, the keyword 'which' refers to
1428 the shifted eigenvalues ``w'[i]`` where:
1430 if mode == 'normal', ``w'[i] = 1 / (w[i] - sigma)``.
1432 if mode == 'cayley', ``w'[i] = (w[i] + sigma) / (w[i] - sigma)``.
1434 if mode == 'buckling', ``w'[i] = w[i] / (w[i] - sigma)``.
1436 (see further discussion in 'mode' below)
1437 v0 : ndarray, optional
1438 Starting vector for iteration.
1439 Default: random
1440 ncv : int, optional
1441 The number of Lanczos vectors generated ncv must be greater than k and
1442 smaller than n; it is recommended that ``ncv > 2*k``.
1443 Default: ``min(n, max(2*k + 1, 20))``
1444 which : str ['LM' | 'SM' | 'LA' | 'SA' | 'BE']
1445 If A is a complex Hermitian matrix, 'BE' is invalid.
1446 Which `k` eigenvectors and eigenvalues to find:
1448 'LM' : Largest (in magnitude) eigenvalues.
1450 'SM' : Smallest (in magnitude) eigenvalues.
1452 'LA' : Largest (algebraic) eigenvalues.
1454 'SA' : Smallest (algebraic) eigenvalues.
1456 'BE' : Half (k/2) from each end of the spectrum.
1458 When k is odd, return one more (k/2+1) from the high end.
1459 When sigma != None, 'which' refers to the shifted eigenvalues ``w'[i]``
1460 (see discussion in 'sigma', above). ARPACK is generally better
1461 at finding large values than small values. If small eigenvalues are
1462 desired, consider using shift-invert mode for better performance.
1463 maxiter : int, optional
1464 Maximum number of Arnoldi update iterations allowed.
1465 Default: ``n*10``
1466 tol : float
1467 Relative accuracy for eigenvalues (stopping criterion).
1468 The default value of 0 implies machine precision.
1469 Minv : N x N matrix, array, sparse matrix, or LinearOperator
1470 See notes in M, above.
1471 OPinv : N x N matrix, array, sparse matrix, or LinearOperator
1472 See notes in sigma, above.
1473 return_eigenvectors : bool
1474 Return eigenvectors (True) in addition to eigenvalues.
1475 This value determines the order in which eigenvalues are sorted.
1476 The sort order is also dependent on the `which` variable.
1478 For which = 'LM' or 'SA':
1479 If `return_eigenvectors` is True, eigenvalues are sorted by
1480 algebraic value.
1482 If `return_eigenvectors` is False, eigenvalues are sorted by
1483 absolute value.
1485 For which = 'BE' or 'LA':
1486 eigenvalues are always sorted by algebraic value.
1488 For which = 'SM':
1489 If `return_eigenvectors` is True, eigenvalues are sorted by
1490 algebraic value.
1492 If `return_eigenvectors` is False, eigenvalues are sorted by
1493 decreasing absolute value.
1495 mode : string ['normal' | 'buckling' | 'cayley']
1496 Specify strategy to use for shift-invert mode. This argument applies
1497 only for real-valued A and sigma != None. For shift-invert mode,
1498 ARPACK internally solves the eigenvalue problem
1499 ``OP @ x'[i] = w'[i] * B @ x'[i]``
1500 and transforms the resulting Ritz vectors x'[i] and Ritz values w'[i]
1501 into the desired eigenvectors and eigenvalues of the problem
1502 ``A @ x[i] = w[i] * M @ x[i]``.
1503 The modes are as follows:
1505 'normal' :
1506 OP = [A - sigma * M]^-1 @ M,
1507 B = M,
1508 w'[i] = 1 / (w[i] - sigma)
1510 'buckling' :
1511 OP = [A - sigma * M]^-1 @ A,
1512 B = A,
1513 w'[i] = w[i] / (w[i] - sigma)
1515 'cayley' :
1516 OP = [A - sigma * M]^-1 @ [A + sigma * M],
1517 B = M,
1518 w'[i] = (w[i] + sigma) / (w[i] - sigma)
1520 The choice of mode will affect which eigenvalues are selected by
1521 the keyword 'which', and can also impact the stability of
1522 convergence (see [2] for a discussion).
1524 Raises
1525 ------
1526 ArpackNoConvergence
1527 When the requested convergence is not obtained.
1529 The currently converged eigenvalues and eigenvectors can be found
1530 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
1531 object.
1533 See Also
1534 --------
1535 eigs : eigenvalues and eigenvectors for a general (nonsymmetric) matrix A
1536 svds : singular value decomposition for a matrix A
1538 Notes
1539 -----
1540 This function is a wrapper to the ARPACK [1]_ SSEUPD and DSEUPD
1541 functions which use the Implicitly Restarted Lanczos Method to
1542 find the eigenvalues and eigenvectors [2]_.
1544 References
1545 ----------
1546 .. [1] ARPACK Software, http://www.caam.rice.edu/software/ARPACK/
1547 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE:
1548 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted
1549 Arnoldi Methods. SIAM, Philadelphia, PA, 1998.
1551 Examples
1552 --------
1553 >>> import numpy as np
1554 >>> from scipy.sparse.linalg import eigsh
1555 >>> identity = np.eye(13)
1556 >>> eigenvalues, eigenvectors = eigsh(identity, k=6)
1557 >>> eigenvalues
1558 array([1., 1., 1., 1., 1., 1.])
1559 >>> eigenvectors.shape
1560 (13, 6)
1562 """
1563 # complex Hermitian matrices should be solved with eigs
1564 if np.issubdtype(A.dtype, np.complexfloating):
1565 if mode != 'normal':
1566 raise ValueError("mode=%s cannot be used with "
1567 "complex matrix A" % mode)
1568 if which == 'BE':
1569 raise ValueError("which='BE' cannot be used with complex matrix A")
1570 elif which == 'LA':
1571 which = 'LR'
1572 elif which == 'SA':
1573 which = 'SR'
1574 ret = eigs(A, k, M=M, sigma=sigma, which=which, v0=v0,
1575 ncv=ncv, maxiter=maxiter, tol=tol,
1576 return_eigenvectors=return_eigenvectors, Minv=Minv,
1577 OPinv=OPinv)
1579 if return_eigenvectors:
1580 return ret[0].real, ret[1]
1581 else:
1582 return ret.real
1584 if A.shape[0] != A.shape[1]:
1585 raise ValueError('expected square matrix (shape=%s)' % (A.shape,))
1586 if M is not None:
1587 if M.shape != A.shape:
1588 raise ValueError('wrong M dimensions %s, should be %s'
1589 % (M.shape, A.shape))
1590 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower():
1591 warnings.warn('M does not have the same type precision as A. '
1592 'This may adversely affect ARPACK convergence')
1594 n = A.shape[0]
1596 if k <= 0:
1597 raise ValueError("k must be greater than 0.")
1599 if k >= n:
1600 warnings.warn("k >= N for N * N square matrix. "
1601 "Attempting to use scipy.linalg.eigh instead.",
1602 RuntimeWarning)
1604 if issparse(A):
1605 raise TypeError("Cannot use scipy.linalg.eigh for sparse A with "
1606 "k >= N. Use scipy.linalg.eigh(A.toarray()) or"
1607 " reduce k.")
1608 if isinstance(A, LinearOperator):
1609 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator "
1610 "A with k >= N.")
1611 if isinstance(M, LinearOperator):
1612 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator "
1613 "M with k >= N.")
1615 return eigh(A, b=M, eigvals_only=not return_eigenvectors)
1617 if sigma is None:
1618 A = _aslinearoperator_with_dtype(A)
1619 matvec = A.matvec
1621 if OPinv is not None:
1622 raise ValueError("OPinv should not be specified "
1623 "with sigma = None.")
1624 if M is None:
1625 #standard eigenvalue problem
1626 mode = 1
1627 M_matvec = None
1628 Minv_matvec = None
1629 if Minv is not None:
1630 raise ValueError("Minv should not be "
1631 "specified with M = None.")
1632 else:
1633 #general eigenvalue problem
1634 mode = 2
1635 if Minv is None:
1636 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol)
1637 else:
1638 Minv = _aslinearoperator_with_dtype(Minv)
1639 Minv_matvec = Minv.matvec
1640 M_matvec = _aslinearoperator_with_dtype(M).matvec
1641 else:
1642 # sigma is not None: shift-invert mode
1643 if Minv is not None:
1644 raise ValueError("Minv should not be specified when sigma is")
1646 # normal mode
1647 if mode == 'normal':
1648 mode = 3
1649 matvec = None
1650 if OPinv is None:
1651 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1652 hermitian=True, tol=tol)
1653 else:
1654 OPinv = _aslinearoperator_with_dtype(OPinv)
1655 Minv_matvec = OPinv.matvec
1656 if M is None:
1657 M_matvec = None
1658 else:
1659 M = _aslinearoperator_with_dtype(M)
1660 M_matvec = M.matvec
1662 # buckling mode
1663 elif mode == 'buckling':
1664 mode = 4
1665 if OPinv is None:
1666 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1667 hermitian=True, tol=tol)
1668 else:
1669 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec
1670 matvec = _aslinearoperator_with_dtype(A).matvec
1671 M_matvec = None
1673 # cayley-transform mode
1674 elif mode == 'cayley':
1675 mode = 5
1676 matvec = _aslinearoperator_with_dtype(A).matvec
1677 if OPinv is None:
1678 Minv_matvec = get_OPinv_matvec(A, M, sigma,
1679 hermitian=True, tol=tol)
1680 else:
1681 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec
1682 if M is None:
1683 M_matvec = None
1684 else:
1685 M_matvec = _aslinearoperator_with_dtype(M).matvec
1687 # unrecognized mode
1688 else:
1689 raise ValueError("unrecognized mode '%s'" % mode)
1691 params = _SymmetricArpackParams(n, k, A.dtype.char, matvec, mode,
1692 M_matvec, Minv_matvec, sigma,
1693 ncv, v0, maxiter, which, tol)
1695 with _ARPACK_LOCK:
1696 while not params.converged:
1697 params.iterate()
1699 return params.extract(return_eigenvectors)