Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_decomp.py: 7%
393 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +0000
1#
2# Author: Pearu Peterson, March 2002
3#
4# additions by Travis Oliphant, March 2002
5# additions by Eric Jones, June 2002
6# additions by Johannes Loehnert, June 2006
7# additions by Bart Vandereycken, June 2006
8# additions by Andrew D Straw, May 2007
9# additions by Tiziano Zito, November 2008
10#
11# April 2010: Functions for LU, QR, SVD, Schur, and Cholesky decompositions
12# were moved to their own files. Still in this file are functions for
13# eigenstuff and for the Hessenberg form.
15__all__ = ['eig', 'eigvals', 'eigh', 'eigvalsh',
16 'eig_banded', 'eigvals_banded',
17 'eigh_tridiagonal', 'eigvalsh_tridiagonal', 'hessenberg', 'cdf2rdf']
19import warnings
21import numpy
22from numpy import (array, isfinite, inexact, nonzero, iscomplexobj,
23 flatnonzero, conj, asarray, argsort, empty,
24 iscomplex, zeros, einsum, eye, inf)
25# Local imports
26from scipy._lib._util import _asarray_validated
27from ._misc import LinAlgError, _datacopied, norm
28from .lapack import get_lapack_funcs, _compute_lwork
29from scipy._lib.deprecation import _NoValue, _deprecate_positional_args
32_I = numpy.array(1j, dtype='F')
35def _make_complex_eigvecs(w, vin, dtype):
36 """
37 Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output
38 """
39 # - see LAPACK man page DGGEV at ALPHAI
40 v = numpy.array(vin, dtype=dtype)
41 m = (w.imag > 0)
42 m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709
43 for i in flatnonzero(m):
44 v.imag[:, i] = vin[:, i+1]
45 conj(v[:, i], v[:, i+1])
46 return v
49def _make_eigvals(alpha, beta, homogeneous_eigvals):
50 if homogeneous_eigvals:
51 if beta is None:
52 return numpy.vstack((alpha, numpy.ones_like(alpha)))
53 else:
54 return numpy.vstack((alpha, beta))
55 else:
56 if beta is None:
57 return alpha
58 else:
59 w = numpy.empty_like(alpha)
60 alpha_zero = (alpha == 0)
61 beta_zero = (beta == 0)
62 beta_nonzero = ~beta_zero
63 w[beta_nonzero] = alpha[beta_nonzero]/beta[beta_nonzero]
64 # Use numpy.inf for complex values too since
65 # 1/numpy.inf = 0, i.e., it correctly behaves as projective
66 # infinity.
67 w[~alpha_zero & beta_zero] = numpy.inf
68 if numpy.all(alpha.imag == 0):
69 w[alpha_zero & beta_zero] = numpy.nan
70 else:
71 w[alpha_zero & beta_zero] = complex(numpy.nan, numpy.nan)
72 return w
75def _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
76 homogeneous_eigvals):
77 ggev, = get_lapack_funcs(('ggev',), (a1, b1))
78 cvl, cvr = left, right
79 res = ggev(a1, b1, lwork=-1)
80 lwork = res[-2][0].real.astype(numpy.int_)
81 if ggev.typecode in 'cz':
82 alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
83 overwrite_a, overwrite_b)
84 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
85 else:
86 alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr,
87 lwork, overwrite_a,
88 overwrite_b)
89 alpha = alphar + _I * alphai
90 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
91 _check_info(info, 'generalized eig algorithm (ggev)')
93 only_real = numpy.all(w.imag == 0.0)
94 if not (ggev.typecode in 'cz' or only_real):
95 t = w.dtype.char
96 if left:
97 vl = _make_complex_eigvecs(w, vl, t)
98 if right:
99 vr = _make_complex_eigvecs(w, vr, t)
101 # the eigenvectors returned by the lapack function are NOT normalized
102 for i in range(vr.shape[0]):
103 if right:
104 vr[:, i] /= norm(vr[:, i])
105 if left:
106 vl[:, i] /= norm(vl[:, i])
108 if not (left or right):
109 return w
110 if left:
111 if right:
112 return w, vl, vr
113 return w, vl
114 return w, vr
117def eig(a, b=None, left=False, right=True, overwrite_a=False,
118 overwrite_b=False, check_finite=True, homogeneous_eigvals=False):
119 """
120 Solve an ordinary or generalized eigenvalue problem of a square matrix.
122 Find eigenvalues w and right or left eigenvectors of a general matrix::
124 a vr[:,i] = w[i] b vr[:,i]
125 a.H vl[:,i] = w[i].conj() b.H vl[:,i]
127 where ``.H`` is the Hermitian conjugation.
129 Parameters
130 ----------
131 a : (M, M) array_like
132 A complex or real matrix whose eigenvalues and eigenvectors
133 will be computed.
134 b : (M, M) array_like, optional
135 Right-hand side matrix in a generalized eigenvalue problem.
136 Default is None, identity matrix is assumed.
137 left : bool, optional
138 Whether to calculate and return left eigenvectors. Default is False.
139 right : bool, optional
140 Whether to calculate and return right eigenvectors. Default is True.
141 overwrite_a : bool, optional
142 Whether to overwrite `a`; may improve performance. Default is False.
143 overwrite_b : bool, optional
144 Whether to overwrite `b`; may improve performance. Default is False.
145 check_finite : bool, optional
146 Whether to check that the input matrices contain only finite numbers.
147 Disabling may give a performance gain, but may result in problems
148 (crashes, non-termination) if the inputs do contain infinities or NaNs.
149 homogeneous_eigvals : bool, optional
150 If True, return the eigenvalues in homogeneous coordinates.
151 In this case ``w`` is a (2, M) array so that::
153 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
155 Default is False.
157 Returns
158 -------
159 w : (M,) or (2, M) double or complex ndarray
160 The eigenvalues, each repeated according to its
161 multiplicity. The shape is (M,) unless
162 ``homogeneous_eigvals=True``.
163 vl : (M, M) double or complex ndarray
164 The normalized left eigenvector corresponding to the eigenvalue
165 ``w[i]`` is the column ``vl[:,i]``. Only returned if ``left=True``.
166 vr : (M, M) double or complex ndarray
167 The normalized right eigenvector corresponding to the eigenvalue
168 ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
170 Raises
171 ------
172 LinAlgError
173 If eigenvalue computation does not converge.
175 See Also
176 --------
177 eigvals : eigenvalues of general arrays
178 eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
179 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
180 band matrices
181 eigh_tridiagonal : eigenvalues and right eiegenvectors for
182 symmetric/Hermitian tridiagonal matrices
184 Examples
185 --------
186 >>> import numpy as np
187 >>> from scipy import linalg
188 >>> a = np.array([[0., -1.], [1., 0.]])
189 >>> linalg.eigvals(a)
190 array([0.+1.j, 0.-1.j])
192 >>> b = np.array([[0., 1.], [1., 1.]])
193 >>> linalg.eigvals(a, b)
194 array([ 1.+0.j, -1.+0.j])
196 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
197 >>> linalg.eigvals(a, homogeneous_eigvals=True)
198 array([[3.+0.j, 8.+0.j, 7.+0.j],
199 [1.+0.j, 1.+0.j, 1.+0.j]])
201 >>> a = np.array([[0., -1.], [1., 0.]])
202 >>> linalg.eigvals(a) == linalg.eig(a)[0]
203 array([ True, True])
204 >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector
205 array([[-0.70710678+0.j , -0.70710678-0.j ],
206 [-0. +0.70710678j, -0. -0.70710678j]])
207 >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector
208 array([[0.70710678+0.j , 0.70710678-0.j ],
209 [0. -0.70710678j, 0. +0.70710678j]])
213 """
214 a1 = _asarray_validated(a, check_finite=check_finite)
215 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
216 raise ValueError('expected square matrix')
217 overwrite_a = overwrite_a or (_datacopied(a1, a))
218 if b is not None:
219 b1 = _asarray_validated(b, check_finite=check_finite)
220 overwrite_b = overwrite_b or _datacopied(b1, b)
221 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
222 raise ValueError('expected square matrix')
223 if b1.shape != a1.shape:
224 raise ValueError('a and b must have the same shape')
225 return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
226 homogeneous_eigvals)
228 geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
229 compute_vl, compute_vr = left, right
231 lwork = _compute_lwork(geev_lwork, a1.shape[0],
232 compute_vl=compute_vl,
233 compute_vr=compute_vr)
235 if geev.typecode in 'cz':
236 w, vl, vr, info = geev(a1, lwork=lwork,
237 compute_vl=compute_vl,
238 compute_vr=compute_vr,
239 overwrite_a=overwrite_a)
240 w = _make_eigvals(w, None, homogeneous_eigvals)
241 else:
242 wr, wi, vl, vr, info = geev(a1, lwork=lwork,
243 compute_vl=compute_vl,
244 compute_vr=compute_vr,
245 overwrite_a=overwrite_a)
246 t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
247 w = wr + _I * wi
248 w = _make_eigvals(w, None, homogeneous_eigvals)
250 _check_info(info, 'eig algorithm (geev)',
251 positive='did not converge (only eigenvalues '
252 'with order >= %d have converged)')
254 only_real = numpy.all(w.imag == 0.0)
255 if not (geev.typecode in 'cz' or only_real):
256 t = w.dtype.char
257 if left:
258 vl = _make_complex_eigvecs(w, vl, t)
259 if right:
260 vr = _make_complex_eigvecs(w, vr, t)
261 if not (left or right):
262 return w
263 if left:
264 if right:
265 return w, vl, vr
266 return w, vl
267 return w, vr
270@_deprecate_positional_args(version="1.14.0")
271def eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False,
272 overwrite_b=False, turbo=_NoValue, eigvals=_NoValue, type=1,
273 check_finite=True, subset_by_index=None, subset_by_value=None,
274 driver=None):
275 """
276 Solve a standard or generalized eigenvalue problem for a complex
277 Hermitian or real symmetric matrix.
279 Find eigenvalues array ``w`` and optionally eigenvectors array ``v`` of
280 array ``a``, where ``b`` is positive definite such that for every
281 eigenvalue λ (i-th entry of w) and its eigenvector ``vi`` (i-th column of
282 ``v``) satisfies::
284 a @ vi = λ * b @ vi
285 vi.conj().T @ a @ vi = λ
286 vi.conj().T @ b @ vi = 1
288 In the standard problem, ``b`` is assumed to be the identity matrix.
290 Parameters
291 ----------
292 a : (M, M) array_like
293 A complex Hermitian or real symmetric matrix whose eigenvalues and
294 eigenvectors will be computed.
295 b : (M, M) array_like, optional
296 A complex Hermitian or real symmetric definite positive matrix in.
297 If omitted, identity matrix is assumed.
298 lower : bool, optional
299 Whether the pertinent array data is taken from the lower or upper
300 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
301 eigvals_only : bool, optional
302 Whether to calculate only eigenvalues and no eigenvectors.
303 (Default: both are calculated)
304 subset_by_index : iterable, optional
305 If provided, this two-element iterable defines the start and the end
306 indices of the desired eigenvalues (ascending order and 0-indexed).
307 To return only the second smallest to fifth smallest eigenvalues,
308 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
309 available with "evr", "evx", and "gvx" drivers. The entries are
310 directly converted to integers via ``int()``.
311 subset_by_value : iterable, optional
312 If provided, this two-element iterable defines the half-open interval
313 ``(a, b]`` that, if any, only the eigenvalues between these values
314 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
315 ``np.inf`` for the unconstrained ends.
316 driver : str, optional
317 Defines which LAPACK driver should be used. Valid options are "ev",
318 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
319 generalized (where b is not None) problems. See the Notes section.
320 The default for standard problems is "evr". For generalized problems,
321 "gvd" is used for full set, and "gvx" for subset requested cases.
322 type : int, optional
323 For the generalized problems, this keyword specifies the problem type
324 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
325 inputs)::
327 1 => a @ v = w @ b @ v
328 2 => a @ b @ v = w @ v
329 3 => b @ a @ v = w @ v
331 This keyword is ignored for standard problems.
332 overwrite_a : bool, optional
333 Whether to overwrite data in ``a`` (may improve performance). Default
334 is False.
335 overwrite_b : bool, optional
336 Whether to overwrite data in ``b`` (may improve performance). Default
337 is False.
338 check_finite : bool, optional
339 Whether to check that the input matrices contain only finite numbers.
340 Disabling may give a performance gain, but may result in problems
341 (crashes, non-termination) if the inputs do contain infinities or NaNs.
342 turbo : bool, optional, deprecated
343 .. deprecated:: 1.5.0
344 `eigh` keyword argument `turbo` is deprecated in favour of
345 ``driver=gvd`` keyword instead and will be removed in SciPy
346 1.14.0.
347 eigvals : tuple (lo, hi), optional, deprecated
348 .. deprecated:: 1.5.0
349 `eigh` keyword argument `eigvals` is deprecated in favour of
350 `subset_by_index` keyword instead and will be removed in SciPy
351 1.14.0.
353 Returns
354 -------
355 w : (N,) ndarray
356 The N (N<=M) selected eigenvalues, in ascending order, each
357 repeated according to its multiplicity.
358 v : (M, N) ndarray
359 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
360 the column ``v[:,i]``. Only returned if ``eigvals_only=False``.
362 Raises
363 ------
364 LinAlgError
365 If eigenvalue computation does not converge, an error occurred, or
366 b matrix is not definite positive. Note that if input matrices are
367 not symmetric or Hermitian, no error will be reported but results will
368 be wrong.
370 See Also
371 --------
372 eigvalsh : eigenvalues of symmetric or Hermitian arrays
373 eig : eigenvalues and right eigenvectors for non-symmetric arrays
374 eigh_tridiagonal : eigenvalues and right eiegenvectors for
375 symmetric/Hermitian tridiagonal matrices
377 Notes
378 -----
379 This function does not check the input array for being Hermitian/symmetric
380 in order to allow for representing arrays with only their upper/lower
381 triangular parts. Also, note that even though not taken into account,
382 finiteness check applies to the whole array and unaffected by "lower"
383 keyword.
385 This function uses LAPACK drivers for computations in all possible keyword
386 combinations, prefixed with ``sy`` if arrays are real and ``he`` if
387 complex, e.g., a float array with "evr" driver is solved via
388 "syevr", complex arrays with "gvx" driver problem is solved via "hegvx"
389 etc.
391 As a brief summary, the slowest and the most robust driver is the
392 classical ``<sy/he>ev`` which uses symmetric QR. ``<sy/he>evr`` is seen as
393 the optimal choice for the most general cases. However, there are certain
394 occasions that ``<sy/he>evd`` computes faster at the expense of more
395 memory usage. ``<sy/he>evx``, while still being faster than ``<sy/he>ev``,
396 often performs worse than the rest except when very few eigenvalues are
397 requested for large arrays though there is still no performance guarantee.
400 For the generalized problem, normalization with respect to the given
401 type argument::
403 type 1 and 3 : v.conj().T @ a @ v = w
404 type 2 : inv(v).conj().T @ a @ inv(v) = w
406 type 1 or 2 : v.conj().T @ b @ v = I
407 type 3 : v.conj().T @ inv(b) @ v = I
410 Examples
411 --------
412 >>> import numpy as np
413 >>> from scipy.linalg import eigh
414 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
415 >>> w, v = eigh(A)
416 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
417 True
419 Request only the eigenvalues
421 >>> w = eigh(A, eigvals_only=True)
423 Request eigenvalues that are less than 10.
425 >>> A = np.array([[34, -4, -10, -7, 2],
426 ... [-4, 7, 2, 12, 0],
427 ... [-10, 2, 44, 2, -19],
428 ... [-7, 12, 2, 79, -34],
429 ... [2, 0, -19, -34, 29]])
430 >>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
431 array([6.69199443e-07, 9.11938152e+00])
433 Request the second smallest eigenvalue and its eigenvector
435 >>> w, v = eigh(A, subset_by_index=[1, 1])
436 >>> w
437 array([9.11938152])
438 >>> v.shape # only a single column is returned
439 (5, 1)
441 """
442 if turbo is not _NoValue:
443 warnings.warn("Keyword argument 'turbo' is deprecated in favour of '"
444 "driver=gvd' keyword instead and will be removed in "
445 "SciPy 1.14.0.",
446 DeprecationWarning, stacklevel=2)
447 if eigvals is not _NoValue:
448 warnings.warn("Keyword argument 'eigvals' is deprecated in favour of "
449 "'subset_by_index' keyword instead and will be removed "
450 "in SciPy 1.14.0.",
451 DeprecationWarning, stacklevel=2)
453 # set lower
454 uplo = 'L' if lower else 'U'
455 # Set job for Fortran routines
456 _job = 'N' if eigvals_only else 'V'
458 drv_str = [None, "ev", "evd", "evr", "evx", "gv", "gvd", "gvx"]
459 if driver not in drv_str:
460 raise ValueError('"{}" is unknown. Possible values are "None", "{}".'
461 ''.format(driver, '", "'.join(drv_str[1:])))
463 a1 = _asarray_validated(a, check_finite=check_finite)
464 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
465 raise ValueError('expected square "a" matrix')
466 overwrite_a = overwrite_a or (_datacopied(a1, a))
467 cplx = True if iscomplexobj(a1) else False
468 n = a1.shape[0]
469 drv_args = {'overwrite_a': overwrite_a}
471 if b is not None:
472 b1 = _asarray_validated(b, check_finite=check_finite)
473 overwrite_b = overwrite_b or _datacopied(b1, b)
474 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
475 raise ValueError('expected square "b" matrix')
477 if b1.shape != a1.shape:
478 raise ValueError("wrong b dimensions {}, should "
479 "be {}".format(b1.shape, a1.shape))
481 if type not in [1, 2, 3]:
482 raise ValueError('"type" keyword only accepts 1, 2, and 3.')
484 cplx = True if iscomplexobj(b1) else (cplx or False)
485 drv_args.update({'overwrite_b': overwrite_b, 'itype': type})
487 # backwards-compatibility handling
488 subset_by_index = subset_by_index if (eigvals in (None, _NoValue)) else eigvals
490 subset = (subset_by_index is not None) or (subset_by_value is not None)
492 # Both subsets can't be given
493 if subset_by_index and subset_by_value:
494 raise ValueError('Either index or value subset can be requested.')
496 # Take turbo into account if all conditions are met otherwise ignore
497 if turbo not in (None, _NoValue) and b is not None:
498 driver = 'gvx' if subset else 'gvd'
500 # Check indices if given
501 if subset_by_index:
502 lo, hi = (int(x) for x in subset_by_index)
503 if not (0 <= lo <= hi < n):
504 raise ValueError('Requested eigenvalue indices are not valid. '
505 'Valid range is [0, {}] and start <= end, but '
506 'start={}, end={} is given'.format(n-1, lo, hi))
507 # fortran is 1-indexed
508 drv_args.update({'range': 'I', 'il': lo + 1, 'iu': hi + 1})
510 if subset_by_value:
511 lo, hi = subset_by_value
512 if not (-inf <= lo < hi <= inf):
513 raise ValueError('Requested eigenvalue bounds are not valid. '
514 'Valid range is (-inf, inf) and low < high, but '
515 'low={}, high={} is given'.format(lo, hi))
517 drv_args.update({'range': 'V', 'vl': lo, 'vu': hi})
519 # fix prefix for lapack routines
520 pfx = 'he' if cplx else 'sy'
522 # decide on the driver if not given
523 # first early exit on incompatible choice
524 if driver:
525 if b is None and (driver in ["gv", "gvd", "gvx"]):
526 raise ValueError('{} requires input b array to be supplied '
527 'for generalized eigenvalue problems.'
528 ''.format(driver))
529 if (b is not None) and (driver in ['ev', 'evd', 'evr', 'evx']):
530 raise ValueError('"{}" does not accept input b array '
531 'for standard eigenvalue problems.'
532 ''.format(driver))
533 if subset and (driver in ["ev", "evd", "gv", "gvd"]):
534 raise ValueError('"{}" cannot compute subsets of eigenvalues'
535 ''.format(driver))
537 # Default driver is evr and gvd
538 else:
539 driver = "evr" if b is None else ("gvx" if subset else "gvd")
541 lwork_spec = {
542 'syevd': ['lwork', 'liwork'],
543 'syevr': ['lwork', 'liwork'],
544 'heevd': ['lwork', 'liwork', 'lrwork'],
545 'heevr': ['lwork', 'lrwork', 'liwork'],
546 }
548 if b is None: # Standard problem
549 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
550 [a1])
551 clw_args = {'n': n, 'lower': lower}
552 if driver == 'evd':
553 clw_args.update({'compute_v': 0 if _job == "N" else 1})
555 lw = _compute_lwork(drvlw, **clw_args)
556 # Multiple lwork vars
557 if isinstance(lw, tuple):
558 lwork_args = dict(zip(lwork_spec[pfx+driver], lw))
559 else:
560 lwork_args = {'lwork': lw}
562 drv_args.update({'lower': lower, 'compute_v': 0 if _job == "N" else 1})
563 w, v, *other_args, info = drv(a=a1, **drv_args, **lwork_args)
565 else: # Generalized problem
566 # 'gvd' doesn't have lwork query
567 if driver == "gvd":
568 drv = get_lapack_funcs(pfx + "gvd", [a1, b1])
569 lwork_args = {}
570 else:
571 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
572 [a1, b1])
573 # generalized drivers use uplo instead of lower
574 lw = _compute_lwork(drvlw, n, uplo=uplo)
575 lwork_args = {'lwork': lw}
577 drv_args.update({'uplo': uplo, 'jobz': _job})
579 w, v, *other_args, info = drv(a=a1, b=b1, **drv_args, **lwork_args)
581 # m is always the first extra argument
582 w = w[:other_args[0]] if subset else w
583 v = v[:, :other_args[0]] if (subset and not eigvals_only) else v
585 # Check if we had a successful exit
586 if info == 0:
587 if eigvals_only:
588 return w
589 else:
590 return w, v
591 else:
592 if info < -1:
593 raise LinAlgError('Illegal value in argument {} of internal {}'
594 ''.format(-info, drv.typecode + pfx + driver))
595 elif info > n:
596 raise LinAlgError('The leading minor of order {} of B is not '
597 'positive definite. The factorization of B '
598 'could not be completed and no eigenvalues '
599 'or eigenvectors were computed.'.format(info-n))
600 else:
601 drv_err = {'ev': 'The algorithm failed to converge; {} '
602 'off-diagonal elements of an intermediate '
603 'tridiagonal form did not converge to zero.',
604 'evx': '{} eigenvectors failed to converge.',
605 'evd': 'The algorithm failed to compute an eigenvalue '
606 'while working on the submatrix lying in rows '
607 'and columns {0}/{1} through mod({0},{1}).',
608 'evr': 'Internal Error.'
609 }
610 if driver in ['ev', 'gv']:
611 msg = drv_err['ev'].format(info)
612 elif driver in ['evx', 'gvx']:
613 msg = drv_err['evx'].format(info)
614 elif driver in ['evd', 'gvd']:
615 if eigvals_only:
616 msg = drv_err['ev'].format(info)
617 else:
618 msg = drv_err['evd'].format(info, n+1)
619 else:
620 msg = drv_err['evr']
622 raise LinAlgError(msg)
625_conv_dict = {0: 0, 1: 1, 2: 2,
626 'all': 0, 'value': 1, 'index': 2,
627 'a': 0, 'v': 1, 'i': 2}
630def _check_select(select, select_range, max_ev, max_len):
631 """Check that select is valid, convert to Fortran style."""
632 if isinstance(select, str):
633 select = select.lower()
634 try:
635 select = _conv_dict[select]
636 except KeyError as e:
637 raise ValueError('invalid argument for select') from e
638 vl, vu = 0., 1.
639 il = iu = 1
640 if select != 0: # (non-all)
641 sr = asarray(select_range)
642 if sr.ndim != 1 or sr.size != 2 or sr[1] < sr[0]:
643 raise ValueError('select_range must be a 2-element array-like '
644 'in nondecreasing order')
645 if select == 1: # (value)
646 vl, vu = sr
647 if max_ev == 0:
648 max_ev = max_len
649 else: # 2 (index)
650 if sr.dtype.char.lower() not in 'hilqp':
651 raise ValueError('when using select="i", select_range must '
652 'contain integers, got dtype %s (%s)'
653 % (sr.dtype, sr.dtype.char))
654 # translate Python (0 ... N-1) into Fortran (1 ... N) with + 1
655 il, iu = sr + 1
656 if min(il, iu) < 1 or max(il, iu) > max_len:
657 raise ValueError('select_range out of bounds')
658 max_ev = iu - il + 1
659 return select, vl, vu, il, iu, max_ev
662def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
663 select='a', select_range=None, max_ev=0, check_finite=True):
664 """
665 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
667 Find eigenvalues w and optionally right eigenvectors v of a::
669 a v[:,i] = w[i] v[:,i]
670 v.H v = identity
672 The matrix a is stored in a_band either in lower diagonal or upper
673 diagonal ordered form:
675 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
676 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
678 where u is the number of bands above the diagonal.
680 Example of a_band (shape of a is (6,6), u=2)::
682 upper form:
683 * * a02 a13 a24 a35
684 * a01 a12 a23 a34 a45
685 a00 a11 a22 a33 a44 a55
687 lower form:
688 a00 a11 a22 a33 a44 a55
689 a10 a21 a32 a43 a54 *
690 a20 a31 a42 a53 * *
692 Cells marked with * are not used.
694 Parameters
695 ----------
696 a_band : (u+1, M) array_like
697 The bands of the M by M matrix a.
698 lower : bool, optional
699 Is the matrix in the lower form. (Default is upper form)
700 eigvals_only : bool, optional
701 Compute only the eigenvalues and no eigenvectors.
702 (Default: calculate also eigenvectors)
703 overwrite_a_band : bool, optional
704 Discard data in a_band (may enhance performance)
705 select : {'a', 'v', 'i'}, optional
706 Which eigenvalues to calculate
708 ====== ========================================
709 select calculated
710 ====== ========================================
711 'a' All eigenvalues
712 'v' Eigenvalues in the interval (min, max]
713 'i' Eigenvalues with indices min <= i <= max
714 ====== ========================================
715 select_range : (min, max), optional
716 Range of selected eigenvalues
717 max_ev : int, optional
718 For select=='v', maximum number of eigenvalues expected.
719 For other values of select, has no meaning.
721 In doubt, leave this parameter untouched.
723 check_finite : bool, optional
724 Whether to check that the input matrix contains only finite numbers.
725 Disabling may give a performance gain, but may result in problems
726 (crashes, non-termination) if the inputs do contain infinities or NaNs.
728 Returns
729 -------
730 w : (M,) ndarray
731 The eigenvalues, in ascending order, each repeated according to its
732 multiplicity.
733 v : (M, M) float or complex ndarray
734 The normalized eigenvector corresponding to the eigenvalue w[i] is
735 the column v[:,i]. Only returned if ``eigvals_only=False``.
737 Raises
738 ------
739 LinAlgError
740 If eigenvalue computation does not converge.
742 See Also
743 --------
744 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
745 eig : eigenvalues and right eigenvectors of general arrays.
746 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
747 eigh_tridiagonal : eigenvalues and right eigenvectors for
748 symmetric/Hermitian tridiagonal matrices
750 Examples
751 --------
752 >>> import numpy as np
753 >>> from scipy.linalg import eig_banded
754 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
755 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
756 >>> w, v = eig_banded(Ab, lower=True)
757 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
758 True
759 >>> w = eig_banded(Ab, lower=True, eigvals_only=True)
760 >>> w
761 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
763 Request only the eigenvalues between ``[-3, 4]``
765 >>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4])
766 >>> w
767 array([-2.22987175, 3.95222349])
769 """
770 if eigvals_only or overwrite_a_band:
771 a1 = _asarray_validated(a_band, check_finite=check_finite)
772 overwrite_a_band = overwrite_a_band or (_datacopied(a1, a_band))
773 else:
774 a1 = array(a_band)
775 if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
776 raise ValueError("array must not contain infs or NaNs")
777 overwrite_a_band = 1
779 if len(a1.shape) != 2:
780 raise ValueError('expected a 2-D array')
781 select, vl, vu, il, iu, max_ev = _check_select(
782 select, select_range, max_ev, a1.shape[1])
783 del select_range
784 if select == 0:
785 if a1.dtype.char in 'GFD':
786 # FIXME: implement this somewhen, for now go with builtin values
787 # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
788 # or by using calc_lwork.f ???
789 # lwork = calc_lwork.hbevd(bevd.typecode, a1.shape[0], lower)
790 internal_name = 'hbevd'
791 else: # a1.dtype.char in 'fd':
792 # FIXME: implement this somewhen, for now go with builtin values
793 # see above
794 # lwork = calc_lwork.sbevd(bevd.typecode, a1.shape[0], lower)
795 internal_name = 'sbevd'
796 bevd, = get_lapack_funcs((internal_name,), (a1,))
797 w, v, info = bevd(a1, compute_v=not eigvals_only,
798 lower=lower, overwrite_ab=overwrite_a_band)
799 else: # select in [1, 2]
800 if eigvals_only:
801 max_ev = 1
802 # calculate optimal abstol for dsbevx (see manpage)
803 if a1.dtype.char in 'fF': # single precision
804 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),))
805 else:
806 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),))
807 abstol = 2 * lamch('s')
808 if a1.dtype.char in 'GFD':
809 internal_name = 'hbevx'
810 else: # a1.dtype.char in 'gfd'
811 internal_name = 'sbevx'
812 bevx, = get_lapack_funcs((internal_name,), (a1,))
813 w, v, m, ifail, info = bevx(
814 a1, vl, vu, il, iu, compute_v=not eigvals_only, mmax=max_ev,
815 range=select, lower=lower, overwrite_ab=overwrite_a_band,
816 abstol=abstol)
817 # crop off w and v
818 w = w[:m]
819 if not eigvals_only:
820 v = v[:, :m]
821 _check_info(info, internal_name)
823 if eigvals_only:
824 return w
825 return w, v
828def eigvals(a, b=None, overwrite_a=False, check_finite=True,
829 homogeneous_eigvals=False):
830 """
831 Compute eigenvalues from an ordinary or generalized eigenvalue problem.
833 Find eigenvalues of a general matrix::
835 a vr[:,i] = w[i] b vr[:,i]
837 Parameters
838 ----------
839 a : (M, M) array_like
840 A complex or real matrix whose eigenvalues and eigenvectors
841 will be computed.
842 b : (M, M) array_like, optional
843 Right-hand side matrix in a generalized eigenvalue problem.
844 If omitted, identity matrix is assumed.
845 overwrite_a : bool, optional
846 Whether to overwrite data in a (may improve performance)
847 check_finite : bool, optional
848 Whether to check that the input matrices contain only finite numbers.
849 Disabling may give a performance gain, but may result in problems
850 (crashes, non-termination) if the inputs do contain infinities
851 or NaNs.
852 homogeneous_eigvals : bool, optional
853 If True, return the eigenvalues in homogeneous coordinates.
854 In this case ``w`` is a (2, M) array so that::
856 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
858 Default is False.
860 Returns
861 -------
862 w : (M,) or (2, M) double or complex ndarray
863 The eigenvalues, each repeated according to its multiplicity
864 but not in any specific order. The shape is (M,) unless
865 ``homogeneous_eigvals=True``.
867 Raises
868 ------
869 LinAlgError
870 If eigenvalue computation does not converge
872 See Also
873 --------
874 eig : eigenvalues and right eigenvectors of general arrays.
875 eigvalsh : eigenvalues of symmetric or Hermitian arrays
876 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
877 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
878 matrices
880 Examples
881 --------
882 >>> import numpy as np
883 >>> from scipy import linalg
884 >>> a = np.array([[0., -1.], [1., 0.]])
885 >>> linalg.eigvals(a)
886 array([0.+1.j, 0.-1.j])
888 >>> b = np.array([[0., 1.], [1., 1.]])
889 >>> linalg.eigvals(a, b)
890 array([ 1.+0.j, -1.+0.j])
892 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
893 >>> linalg.eigvals(a, homogeneous_eigvals=True)
894 array([[3.+0.j, 8.+0.j, 7.+0.j],
895 [1.+0.j, 1.+0.j, 1.+0.j]])
897 """
898 return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a,
899 check_finite=check_finite,
900 homogeneous_eigvals=homogeneous_eigvals)
903@_deprecate_positional_args(version="1.14.0")
904def eigvalsh(a, b=None, *, lower=True, overwrite_a=False,
905 overwrite_b=False, turbo=_NoValue, eigvals=_NoValue, type=1,
906 check_finite=True, subset_by_index=None, subset_by_value=None,
907 driver=None):
908 """
909 Solves a standard or generalized eigenvalue problem for a complex
910 Hermitian or real symmetric matrix.
912 Find eigenvalues array ``w`` of array ``a``, where ``b`` is positive
913 definite such that for every eigenvalue λ (i-th entry of w) and its
914 eigenvector vi (i-th column of v) satisfies::
916 a @ vi = λ * b @ vi
917 vi.conj().T @ a @ vi = λ
918 vi.conj().T @ b @ vi = 1
920 In the standard problem, b is assumed to be the identity matrix.
922 Parameters
923 ----------
924 a : (M, M) array_like
925 A complex Hermitian or real symmetric matrix whose eigenvalues will
926 be computed.
927 b : (M, M) array_like, optional
928 A complex Hermitian or real symmetric definite positive matrix in.
929 If omitted, identity matrix is assumed.
930 lower : bool, optional
931 Whether the pertinent array data is taken from the lower or upper
932 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
933 overwrite_a : bool, optional
934 Whether to overwrite data in ``a`` (may improve performance). Default
935 is False.
936 overwrite_b : bool, optional
937 Whether to overwrite data in ``b`` (may improve performance). Default
938 is False.
939 type : int, optional
940 For the generalized problems, this keyword specifies the problem type
941 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
942 inputs)::
944 1 => a @ v = w @ b @ v
945 2 => a @ b @ v = w @ v
946 3 => b @ a @ v = w @ v
948 This keyword is ignored for standard problems.
949 check_finite : bool, optional
950 Whether to check that the input matrices contain only finite numbers.
951 Disabling may give a performance gain, but may result in problems
952 (crashes, non-termination) if the inputs do contain infinities or NaNs.
953 subset_by_index : iterable, optional
954 If provided, this two-element iterable defines the start and the end
955 indices of the desired eigenvalues (ascending order and 0-indexed).
956 To return only the second smallest to fifth smallest eigenvalues,
957 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
958 available with "evr", "evx", and "gvx" drivers. The entries are
959 directly converted to integers via ``int()``.
960 subset_by_value : iterable, optional
961 If provided, this two-element iterable defines the half-open interval
962 ``(a, b]`` that, if any, only the eigenvalues between these values
963 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
964 ``np.inf`` for the unconstrained ends.
965 driver : str, optional
966 Defines which LAPACK driver should be used. Valid options are "ev",
967 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
968 generalized (where b is not None) problems. See the Notes section of
969 `scipy.linalg.eigh`.
970 turbo : bool, optional, deprecated
971 .. deprecated:: 1.5.0
972 'eigvalsh' keyword argument `turbo` is deprecated in favor of
973 ``driver=gvd`` option and will be removed in SciPy 1.14.0.
975 eigvals : tuple (lo, hi), optional
976 .. deprecated:: 1.5.0
977 'eigvalsh' keyword argument `eigvals` is deprecated in favor of
978 `subset_by_index` option and will be removed in SciPy 1.14.0.
980 Returns
981 -------
982 w : (N,) ndarray
983 The N (N<=M) selected eigenvalues, in ascending order, each
984 repeated according to its multiplicity.
986 Raises
987 ------
988 LinAlgError
989 If eigenvalue computation does not converge, an error occurred, or
990 b matrix is not definite positive. Note that if input matrices are
991 not symmetric or Hermitian, no error will be reported but results will
992 be wrong.
994 See Also
995 --------
996 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
997 eigvals : eigenvalues of general arrays
998 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
999 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1000 matrices
1002 Notes
1003 -----
1004 This function does not check the input array for being Hermitian/symmetric
1005 in order to allow for representing arrays with only their upper/lower
1006 triangular parts.
1008 This function serves as a one-liner shorthand for `scipy.linalg.eigh` with
1009 the option ``eigvals_only=True`` to get the eigenvalues and not the
1010 eigenvectors. Here it is kept as a legacy convenience. It might be
1011 beneficial to use the main function to have full control and to be a bit
1012 more pythonic.
1014 Examples
1015 --------
1016 For more examples see `scipy.linalg.eigh`.
1018 >>> import numpy as np
1019 >>> from scipy.linalg import eigvalsh
1020 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
1021 >>> w = eigvalsh(A)
1022 >>> w
1023 array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
1025 """
1026 return eigh(a, b=b, lower=lower, eigvals_only=True,
1027 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
1028 turbo=turbo, eigvals=eigvals, type=type,
1029 check_finite=check_finite, subset_by_index=subset_by_index,
1030 subset_by_value=subset_by_value, driver=driver)
1033def eigvals_banded(a_band, lower=False, overwrite_a_band=False,
1034 select='a', select_range=None, check_finite=True):
1035 """
1036 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
1038 Find eigenvalues w of a::
1040 a v[:,i] = w[i] v[:,i]
1041 v.H v = identity
1043 The matrix a is stored in a_band either in lower diagonal or upper
1044 diagonal ordered form:
1046 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
1047 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
1049 where u is the number of bands above the diagonal.
1051 Example of a_band (shape of a is (6,6), u=2)::
1053 upper form:
1054 * * a02 a13 a24 a35
1055 * a01 a12 a23 a34 a45
1056 a00 a11 a22 a33 a44 a55
1058 lower form:
1059 a00 a11 a22 a33 a44 a55
1060 a10 a21 a32 a43 a54 *
1061 a20 a31 a42 a53 * *
1063 Cells marked with * are not used.
1065 Parameters
1066 ----------
1067 a_band : (u+1, M) array_like
1068 The bands of the M by M matrix a.
1069 lower : bool, optional
1070 Is the matrix in the lower form. (Default is upper form)
1071 overwrite_a_band : bool, optional
1072 Discard data in a_band (may enhance performance)
1073 select : {'a', 'v', 'i'}, optional
1074 Which eigenvalues to calculate
1076 ====== ========================================
1077 select calculated
1078 ====== ========================================
1079 'a' All eigenvalues
1080 'v' Eigenvalues in the interval (min, max]
1081 'i' Eigenvalues with indices min <= i <= max
1082 ====== ========================================
1083 select_range : (min, max), optional
1084 Range of selected eigenvalues
1085 check_finite : bool, optional
1086 Whether to check that the input matrix contains only finite numbers.
1087 Disabling may give a performance gain, but may result in problems
1088 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1090 Returns
1091 -------
1092 w : (M,) ndarray
1093 The eigenvalues, in ascending order, each repeated according to its
1094 multiplicity.
1096 Raises
1097 ------
1098 LinAlgError
1099 If eigenvalue computation does not converge.
1101 See Also
1102 --------
1103 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1104 band matrices
1105 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1106 matrices
1107 eigvals : eigenvalues of general arrays
1108 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1109 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1111 Examples
1112 --------
1113 >>> import numpy as np
1114 >>> from scipy.linalg import eigvals_banded
1115 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
1116 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
1117 >>> w = eigvals_banded(Ab, lower=True)
1118 >>> w
1119 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
1120 """
1121 return eig_banded(a_band, lower=lower, eigvals_only=1,
1122 overwrite_a_band=overwrite_a_band, select=select,
1123 select_range=select_range, check_finite=check_finite)
1126def eigvalsh_tridiagonal(d, e, select='a', select_range=None,
1127 check_finite=True, tol=0., lapack_driver='auto'):
1128 """
1129 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1131 Find eigenvalues `w` of ``a``::
1133 a v[:,i] = w[i] v[:,i]
1134 v.H v = identity
1136 For a real symmetric matrix ``a`` with diagonal elements `d` and
1137 off-diagonal elements `e`.
1139 Parameters
1140 ----------
1141 d : ndarray, shape (ndim,)
1142 The diagonal elements of the array.
1143 e : ndarray, shape (ndim-1,)
1144 The off-diagonal elements of the array.
1145 select : {'a', 'v', 'i'}, optional
1146 Which eigenvalues to calculate
1148 ====== ========================================
1149 select calculated
1150 ====== ========================================
1151 'a' All eigenvalues
1152 'v' Eigenvalues in the interval (min, max]
1153 'i' Eigenvalues with indices min <= i <= max
1154 ====== ========================================
1155 select_range : (min, max), optional
1156 Range of selected eigenvalues
1157 check_finite : bool, optional
1158 Whether to check that the input matrix contains only finite numbers.
1159 Disabling may give a performance gain, but may result in problems
1160 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1161 tol : float
1162 The absolute tolerance to which each eigenvalue is required
1163 (only used when ``lapack_driver='stebz'``).
1164 An eigenvalue (or cluster) is considered to have converged if it
1165 lies in an interval of this width. If <= 0. (default),
1166 the value ``eps*|a|`` is used where eps is the machine precision,
1167 and ``|a|`` is the 1-norm of the matrix ``a``.
1168 lapack_driver : str
1169 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1170 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1171 and 'stebz' otherwise. 'sterf' and 'stev' can only be used when
1172 ``select='a'``.
1174 Returns
1175 -------
1176 w : (M,) ndarray
1177 The eigenvalues, in ascending order, each repeated according to its
1178 multiplicity.
1180 Raises
1181 ------
1182 LinAlgError
1183 If eigenvalue computation does not converge.
1185 See Also
1186 --------
1187 eigh_tridiagonal : eigenvalues and right eiegenvectors for
1188 symmetric/Hermitian tridiagonal matrices
1190 Examples
1191 --------
1192 >>> import numpy as np
1193 >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
1194 >>> d = 3*np.ones(4)
1195 >>> e = -1*np.ones(3)
1196 >>> w = eigvalsh_tridiagonal(d, e)
1197 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1198 >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
1199 >>> np.allclose(w - w2, np.zeros(4))
1200 True
1201 """
1202 return eigh_tridiagonal(
1203 d, e, eigvals_only=True, select=select, select_range=select_range,
1204 check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
1207def eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None,
1208 check_finite=True, tol=0., lapack_driver='auto'):
1209 """
1210 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1212 Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
1214 a v[:,i] = w[i] v[:,i]
1215 v.H v = identity
1217 For a real symmetric matrix ``a`` with diagonal elements `d` and
1218 off-diagonal elements `e`.
1220 Parameters
1221 ----------
1222 d : ndarray, shape (ndim,)
1223 The diagonal elements of the array.
1224 e : ndarray, shape (ndim-1,)
1225 The off-diagonal elements of the array.
1226 eigvals_only : bool, optional
1227 Compute only the eigenvalues and no eigenvectors.
1228 (Default: calculate also eigenvectors)
1229 select : {'a', 'v', 'i'}, optional
1230 Which eigenvalues to calculate
1232 ====== ========================================
1233 select calculated
1234 ====== ========================================
1235 'a' All eigenvalues
1236 'v' Eigenvalues in the interval (min, max]
1237 'i' Eigenvalues with indices min <= i <= max
1238 ====== ========================================
1239 select_range : (min, max), optional
1240 Range of selected eigenvalues
1241 check_finite : bool, optional
1242 Whether to check that the input matrix contains only finite numbers.
1243 Disabling may give a performance gain, but may result in problems
1244 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1245 tol : float
1246 The absolute tolerance to which each eigenvalue is required
1247 (only used when 'stebz' is the `lapack_driver`).
1248 An eigenvalue (or cluster) is considered to have converged if it
1249 lies in an interval of this width. If <= 0. (default),
1250 the value ``eps*|a|`` is used where eps is the machine precision,
1251 and ``|a|`` is the 1-norm of the matrix ``a``.
1252 lapack_driver : str
1253 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1254 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1255 and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and
1256 ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is
1257 used to find the corresponding eigenvectors. 'sterf' can only be
1258 used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only
1259 be used when ``select='a'``.
1261 Returns
1262 -------
1263 w : (M,) ndarray
1264 The eigenvalues, in ascending order, each repeated according to its
1265 multiplicity.
1266 v : (M, M) ndarray
1267 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
1268 the column ``v[:,i]``. Only returned if ``eigvals_only=False``.
1270 Raises
1271 ------
1272 LinAlgError
1273 If eigenvalue computation does not converge.
1275 See Also
1276 --------
1277 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1278 matrices
1279 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1280 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1281 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1282 band matrices
1284 Notes
1285 -----
1286 This function makes use of LAPACK ``S/DSTEMR`` routines.
1288 Examples
1289 --------
1290 >>> import numpy as np
1291 >>> from scipy.linalg import eigh_tridiagonal
1292 >>> d = 3*np.ones(4)
1293 >>> e = -1*np.ones(3)
1294 >>> w, v = eigh_tridiagonal(d, e)
1295 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1296 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
1297 True
1298 """
1299 d = _asarray_validated(d, check_finite=check_finite)
1300 e = _asarray_validated(e, check_finite=check_finite)
1301 for check in (d, e):
1302 if check.ndim != 1:
1303 raise ValueError('expected a 1-D array')
1304 if check.dtype.char in 'GFD': # complex
1305 raise TypeError('Only real arrays currently supported')
1306 if d.size != e.size + 1:
1307 raise ValueError('d (%s) must have one more element than e (%s)'
1308 % (d.size, e.size))
1309 select, vl, vu, il, iu, _ = _check_select(
1310 select, select_range, 0, d.size)
1311 if not isinstance(lapack_driver, str):
1312 raise TypeError('lapack_driver must be str')
1313 drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev')
1314 if lapack_driver not in drivers:
1315 raise ValueError('lapack_driver must be one of %s, got %s'
1316 % (drivers, lapack_driver))
1317 if lapack_driver == 'auto':
1318 lapack_driver = 'stemr' if select == 0 else 'stebz'
1319 func, = get_lapack_funcs((lapack_driver,), (d, e))
1320 compute_v = not eigvals_only
1321 if lapack_driver == 'sterf':
1322 if select != 0:
1323 raise ValueError('sterf can only be used when select == "a"')
1324 if not eigvals_only:
1325 raise ValueError('sterf can only be used when eigvals_only is '
1326 'True')
1327 w, info = func(d, e)
1328 m = len(w)
1329 elif lapack_driver == 'stev':
1330 if select != 0:
1331 raise ValueError('stev can only be used when select == "a"')
1332 w, v, info = func(d, e, compute_v=compute_v)
1333 m = len(w)
1334 elif lapack_driver == 'stebz':
1335 tol = float(tol)
1336 internal_name = 'stebz'
1337 stebz, = get_lapack_funcs((internal_name,), (d, e))
1338 # If getting eigenvectors, needs to be block-ordered (B) instead of
1339 # matrix-ordered (E), and we will reorder later
1340 order = 'E' if eigvals_only else 'B'
1341 m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol,
1342 order)
1343 else: # 'stemr'
1344 # ?STEMR annoyingly requires size N instead of N-1
1345 e_ = empty(e.size+1, e.dtype)
1346 e_[:-1] = e
1347 stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e))
1348 lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu,
1349 compute_v=compute_v)
1350 _check_info(info, 'stemr_lwork')
1351 m, w, v, info = func(d, e_, select, vl, vu, il, iu,
1352 compute_v=compute_v, lwork=lwork, liwork=liwork)
1353 _check_info(info, lapack_driver + ' (eigh_tridiagonal)')
1354 w = w[:m]
1355 if eigvals_only:
1356 return w
1357 else:
1358 # Do we still need to compute the eigenvalues?
1359 if lapack_driver == 'stebz':
1360 func, = get_lapack_funcs(('stein',), (d, e))
1361 v, info = func(d, e, w, iblock, isplit)
1362 _check_info(info, 'stein (eigh_tridiagonal)',
1363 positive='%d eigenvectors failed to converge')
1364 # Convert block-order to matrix-order
1365 order = argsort(w)
1366 w, v = w[order], v[:, order]
1367 else:
1368 v = v[:, :m]
1369 return w, v
1372def _check_info(info, driver, positive='did not converge (LAPACK info=%d)'):
1373 """Check info return value."""
1374 if info < 0:
1375 raise ValueError('illegal value in argument %d of internal %s'
1376 % (-info, driver))
1377 if info > 0 and positive:
1378 raise LinAlgError(("%s " + positive) % (driver, info,))
1381def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
1382 """
1383 Compute Hessenberg form of a matrix.
1385 The Hessenberg decomposition is::
1387 A = Q H Q^H
1389 where `Q` is unitary/orthogonal and `H` has only zero elements below
1390 the first sub-diagonal.
1392 Parameters
1393 ----------
1394 a : (M, M) array_like
1395 Matrix to bring into Hessenberg form.
1396 calc_q : bool, optional
1397 Whether to compute the transformation matrix. Default is False.
1398 overwrite_a : bool, optional
1399 Whether to overwrite `a`; may improve performance.
1400 Default is False.
1401 check_finite : bool, optional
1402 Whether to check that the input matrix contains only finite numbers.
1403 Disabling may give a performance gain, but may result in problems
1404 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1406 Returns
1407 -------
1408 H : (M, M) ndarray
1409 Hessenberg form of `a`.
1410 Q : (M, M) ndarray
1411 Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
1412 Only returned if ``calc_q=True``.
1414 Examples
1415 --------
1416 >>> import numpy as np
1417 >>> from scipy.linalg import hessenberg
1418 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
1419 >>> H, Q = hessenberg(A, calc_q=True)
1420 >>> H
1421 array([[ 2. , -11.65843866, 1.42005301, 0.25349066],
1422 [ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
1423 [ 0. , -1.83299243, 0.38969961, -0.51527034],
1424 [ 0. , 0. , -3.83189513, 1.07494686]])
1425 >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
1426 True
1427 """
1428 a1 = _asarray_validated(a, check_finite=check_finite)
1429 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
1430 raise ValueError('expected square matrix')
1431 overwrite_a = overwrite_a or (_datacopied(a1, a))
1433 # if 2x2 or smaller: already in Hessenberg
1434 if a1.shape[0] <= 2:
1435 if calc_q:
1436 return a1, eye(a1.shape[0])
1437 return a1
1439 gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal',
1440 'gehrd_lwork'), (a1,))
1441 ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a)
1442 _check_info(info, 'gebal (hessenberg)', positive=False)
1443 n = len(a1)
1445 lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
1447 hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1448 _check_info(info, 'gehrd (hessenberg)', positive=False)
1449 h = numpy.triu(hq, -1)
1450 if not calc_q:
1451 return h
1453 # use orghr/unghr to compute q
1454 orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,))
1455 lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
1457 q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1458 _check_info(info, 'orghr (hessenberg)', positive=False)
1459 return h, q
1462def cdf2rdf(w, v):
1463 """
1464 Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real
1465 eigenvalues in a block diagonal form ``wr`` and the associated real
1466 eigenvectors ``vr``, such that::
1468 vr @ wr = X @ vr
1470 continues to hold, where ``X`` is the original array for which ``w`` and
1471 ``v`` are the eigenvalues and eigenvectors.
1473 .. versionadded:: 1.1.0
1475 Parameters
1476 ----------
1477 w : (..., M) array_like
1478 Complex or real eigenvalues, an array or stack of arrays
1480 Conjugate pairs must not be interleaved, else the wrong result
1481 will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result,
1482 but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
1484 v : (..., M, M) array_like
1485 Complex or real eigenvectors, a square array or stack of square arrays.
1487 Returns
1488 -------
1489 wr : (..., M, M) ndarray
1490 Real diagonal block form of eigenvalues
1491 vr : (..., M, M) ndarray
1492 Real eigenvectors associated with ``wr``
1494 See Also
1495 --------
1496 eig : Eigenvalues and right eigenvectors for non-symmetric arrays
1497 rsf2csf : Convert real Schur form to complex Schur form
1499 Notes
1500 -----
1501 ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``.
1502 For example, obtained by ``w, v = scipy.linalg.eig(X)`` or
1503 ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent
1504 stacked arrays.
1506 .. versionadded:: 1.1.0
1508 Examples
1509 --------
1510 >>> import numpy as np
1511 >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
1512 >>> X
1513 array([[ 1, 2, 3],
1514 [ 0, 4, 5],
1515 [ 0, -5, 4]])
1517 >>> from scipy import linalg
1518 >>> w, v = linalg.eig(X)
1519 >>> w
1520 array([ 1.+0.j, 4.+5.j, 4.-5.j])
1521 >>> v
1522 array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
1523 [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
1524 [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
1526 >>> wr, vr = linalg.cdf2rdf(w, v)
1527 >>> wr
1528 array([[ 1., 0., 0.],
1529 [ 0., 4., 5.],
1530 [ 0., -5., 4.]])
1531 >>> vr
1532 array([[ 1. , 0.40016, -0.01906],
1533 [ 0. , 0.64788, 0. ],
1534 [ 0. , 0. , 0.64788]])
1536 >>> vr @ wr
1537 array([[ 1. , 1.69593, 1.9246 ],
1538 [ 0. , 2.59153, 3.23942],
1539 [ 0. , -3.23942, 2.59153]])
1540 >>> X @ vr
1541 array([[ 1. , 1.69593, 1.9246 ],
1542 [ 0. , 2.59153, 3.23942],
1543 [ 0. , -3.23942, 2.59153]])
1544 """
1545 w, v = _asarray_validated(w), _asarray_validated(v)
1547 # check dimensions
1548 if w.ndim < 1:
1549 raise ValueError('expected w to be at least 1D')
1550 if v.ndim < 2:
1551 raise ValueError('expected v to be at least 2D')
1552 if v.ndim != w.ndim + 1:
1553 raise ValueError('expected eigenvectors array to have exactly one '
1554 'dimension more than eigenvalues array')
1556 # check shapes
1557 n = w.shape[-1]
1558 M = w.shape[:-1]
1559 if v.shape[-2] != v.shape[-1]:
1560 raise ValueError('expected v to be a square matrix or stacked square '
1561 'matrices: v.shape[-2] = v.shape[-1]')
1562 if v.shape[-1] != n:
1563 raise ValueError('expected the same number of eigenvalues as '
1564 'eigenvectors')
1566 # get indices for each first pair of complex eigenvalues
1567 complex_mask = iscomplex(w)
1568 n_complex = complex_mask.sum(axis=-1)
1570 # check if all complex eigenvalues have conjugate pairs
1571 if not (n_complex % 2 == 0).all():
1572 raise ValueError('expected complex-conjugate pairs of eigenvalues')
1574 # find complex indices
1575 idx = nonzero(complex_mask)
1576 idx_stack = idx[:-1]
1577 idx_elem = idx[-1]
1579 # filter them to conjugate indices, assuming pairs are not interleaved
1580 j = idx_elem[0::2]
1581 k = idx_elem[1::2]
1582 stack_ind = ()
1583 for i in idx_stack:
1584 # should never happen, assuming nonzero orders by the last axis
1585 assert (i[0::2] == i[1::2]).all(),\
1586 "Conjugate pair spanned different arrays!"
1587 stack_ind += (i[0::2],)
1589 # all eigenvalues to diagonal form
1590 wr = zeros(M + (n, n), dtype=w.real.dtype)
1591 di = range(n)
1592 wr[..., di, di] = w.real
1594 # complex eigenvalues to real block diagonal form
1595 wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag
1596 wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
1598 # compute real eigenvectors associated with real block diagonal eigenvalues
1599 u = zeros(M + (n, n), dtype=numpy.cdouble)
1600 u[..., di, di] = 1.0
1601 u[stack_ind + (j, j)] = 0.5j
1602 u[stack_ind + (j, k)] = 0.5
1603 u[stack_ind + (k, j)] = -0.5j
1604 u[stack_ind + (k, k)] = 0.5
1606 # multipy matrices v and u (equivalent to v @ u)
1607 vr = einsum('...ij,...jk->...ik', v, u).real
1609 return wr, vr