Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_decomp.py: 7%
402 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-14 06:37 +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 left eigenvector corresponding to the eigenvalue
165 ``w[i]`` is the column ``vl[:,i]``. Only returned if ``left=True``.
166 The left eigenvector is not normalized.
167 vr : (M, M) double or complex ndarray
168 The normalized right eigenvector corresponding to the eigenvalue
169 ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
171 Raises
172 ------
173 LinAlgError
174 If eigenvalue computation does not converge.
176 See Also
177 --------
178 eigvals : eigenvalues of general arrays
179 eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
180 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
181 band matrices
182 eigh_tridiagonal : eigenvalues and right eiegenvectors for
183 symmetric/Hermitian tridiagonal matrices
185 Examples
186 --------
187 >>> import numpy as np
188 >>> from scipy import linalg
189 >>> a = np.array([[0., -1.], [1., 0.]])
190 >>> linalg.eigvals(a)
191 array([0.+1.j, 0.-1.j])
193 >>> b = np.array([[0., 1.], [1., 1.]])
194 >>> linalg.eigvals(a, b)
195 array([ 1.+0.j, -1.+0.j])
197 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
198 >>> linalg.eigvals(a, homogeneous_eigvals=True)
199 array([[3.+0.j, 8.+0.j, 7.+0.j],
200 [1.+0.j, 1.+0.j, 1.+0.j]])
202 >>> a = np.array([[0., -1.], [1., 0.]])
203 >>> linalg.eigvals(a) == linalg.eig(a)[0]
204 array([ True, True])
205 >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector
206 array([[-0.70710678+0.j , -0.70710678-0.j ],
207 [-0. +0.70710678j, -0. -0.70710678j]])
208 >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector
209 array([[0.70710678+0.j , 0.70710678-0.j ],
210 [0. -0.70710678j, 0. +0.70710678j]])
214 """
215 a1 = _asarray_validated(a, check_finite=check_finite)
216 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
217 raise ValueError('expected square matrix')
218 overwrite_a = overwrite_a or (_datacopied(a1, a))
219 if b is not None:
220 b1 = _asarray_validated(b, check_finite=check_finite)
221 overwrite_b = overwrite_b or _datacopied(b1, b)
222 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
223 raise ValueError('expected square matrix')
224 if b1.shape != a1.shape:
225 raise ValueError('a and b must have the same shape')
226 return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
227 homogeneous_eigvals)
229 geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
230 compute_vl, compute_vr = left, right
232 lwork = _compute_lwork(geev_lwork, a1.shape[0],
233 compute_vl=compute_vl,
234 compute_vr=compute_vr)
236 if geev.typecode in 'cz':
237 w, vl, vr, info = geev(a1, lwork=lwork,
238 compute_vl=compute_vl,
239 compute_vr=compute_vr,
240 overwrite_a=overwrite_a)
241 w = _make_eigvals(w, None, homogeneous_eigvals)
242 else:
243 wr, wi, vl, vr, info = geev(a1, lwork=lwork,
244 compute_vl=compute_vl,
245 compute_vr=compute_vr,
246 overwrite_a=overwrite_a)
247 t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
248 w = wr + _I * wi
249 w = _make_eigvals(w, None, homogeneous_eigvals)
251 _check_info(info, 'eig algorithm (geev)',
252 positive='did not converge (only eigenvalues '
253 'with order >= %d have converged)')
255 only_real = numpy.all(w.imag == 0.0)
256 if not (geev.typecode in 'cz' or only_real):
257 t = w.dtype.char
258 if left:
259 vl = _make_complex_eigvecs(w, vl, t)
260 if right:
261 vr = _make_complex_eigvecs(w, vr, t)
262 if not (left or right):
263 return w
264 if left:
265 if right:
266 return w, vl, vr
267 return w, vl
268 return w, vr
271@_deprecate_positional_args(version="1.14.0")
272def eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False,
273 overwrite_b=False, turbo=_NoValue, eigvals=_NoValue, type=1,
274 check_finite=True, subset_by_index=None, subset_by_value=None,
275 driver=None):
276 """
277 Solve a standard or generalized eigenvalue problem for a complex
278 Hermitian or real symmetric matrix.
280 Find eigenvalues array ``w`` and optionally eigenvectors array ``v`` of
281 array ``a``, where ``b`` is positive definite such that for every
282 eigenvalue λ (i-th entry of w) and its eigenvector ``vi`` (i-th column of
283 ``v``) satisfies::
285 a @ vi = λ * b @ vi
286 vi.conj().T @ a @ vi = λ
287 vi.conj().T @ b @ vi = 1
289 In the standard problem, ``b`` is assumed to be the identity matrix.
291 Parameters
292 ----------
293 a : (M, M) array_like
294 A complex Hermitian or real symmetric matrix whose eigenvalues and
295 eigenvectors will be computed.
296 b : (M, M) array_like, optional
297 A complex Hermitian or real symmetric definite positive matrix in.
298 If omitted, identity matrix is assumed.
299 lower : bool, optional
300 Whether the pertinent array data is taken from the lower or upper
301 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
302 eigvals_only : bool, optional
303 Whether to calculate only eigenvalues and no eigenvectors.
304 (Default: both are calculated)
305 subset_by_index : iterable, optional
306 If provided, this two-element iterable defines the start and the end
307 indices of the desired eigenvalues (ascending order and 0-indexed).
308 To return only the second smallest to fifth smallest eigenvalues,
309 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
310 available with "evr", "evx", and "gvx" drivers. The entries are
311 directly converted to integers via ``int()``.
312 subset_by_value : iterable, optional
313 If provided, this two-element iterable defines the half-open interval
314 ``(a, b]`` that, if any, only the eigenvalues between these values
315 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
316 ``np.inf`` for the unconstrained ends.
317 driver : str, optional
318 Defines which LAPACK driver should be used. Valid options are "ev",
319 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
320 generalized (where b is not None) problems. See the Notes section.
321 The default for standard problems is "evr". For generalized problems,
322 "gvd" is used for full set, and "gvx" for subset requested cases.
323 type : int, optional
324 For the generalized problems, this keyword specifies the problem type
325 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
326 inputs)::
328 1 => a @ v = w @ b @ v
329 2 => a @ b @ v = w @ v
330 3 => b @ a @ v = w @ v
332 This keyword is ignored for standard problems.
333 overwrite_a : bool, optional
334 Whether to overwrite data in ``a`` (may improve performance). Default
335 is False.
336 overwrite_b : bool, optional
337 Whether to overwrite data in ``b`` (may improve performance). Default
338 is False.
339 check_finite : bool, optional
340 Whether to check that the input matrices contain only finite numbers.
341 Disabling may give a performance gain, but may result in problems
342 (crashes, non-termination) if the inputs do contain infinities or NaNs.
343 turbo : bool, optional, deprecated
344 .. deprecated:: 1.5.0
345 `eigh` keyword argument `turbo` is deprecated in favour of
346 ``driver=gvd`` keyword instead and will be removed in SciPy
347 1.14.0.
348 eigvals : tuple (lo, hi), optional, deprecated
349 .. deprecated:: 1.5.0
350 `eigh` keyword argument `eigvals` is deprecated in favour of
351 `subset_by_index` keyword instead and will be removed in SciPy
352 1.14.0.
354 Returns
355 -------
356 w : (N,) ndarray
357 The N (N<=M) selected eigenvalues, in ascending order, each
358 repeated according to its multiplicity.
359 v : (M, N) ndarray
360 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
361 the column ``v[:,i]``. Only returned if ``eigvals_only=False``.
363 Raises
364 ------
365 LinAlgError
366 If eigenvalue computation does not converge, an error occurred, or
367 b matrix is not definite positive. Note that if input matrices are
368 not symmetric or Hermitian, no error will be reported but results will
369 be wrong.
371 See Also
372 --------
373 eigvalsh : eigenvalues of symmetric or Hermitian arrays
374 eig : eigenvalues and right eigenvectors for non-symmetric arrays
375 eigh_tridiagonal : eigenvalues and right eiegenvectors for
376 symmetric/Hermitian tridiagonal matrices
378 Notes
379 -----
380 This function does not check the input array for being Hermitian/symmetric
381 in order to allow for representing arrays with only their upper/lower
382 triangular parts. Also, note that even though not taken into account,
383 finiteness check applies to the whole array and unaffected by "lower"
384 keyword.
386 This function uses LAPACK drivers for computations in all possible keyword
387 combinations, prefixed with ``sy`` if arrays are real and ``he`` if
388 complex, e.g., a float array with "evr" driver is solved via
389 "syevr", complex arrays with "gvx" driver problem is solved via "hegvx"
390 etc.
392 As a brief summary, the slowest and the most robust driver is the
393 classical ``<sy/he>ev`` which uses symmetric QR. ``<sy/he>evr`` is seen as
394 the optimal choice for the most general cases. However, there are certain
395 occasions that ``<sy/he>evd`` computes faster at the expense of more
396 memory usage. ``<sy/he>evx``, while still being faster than ``<sy/he>ev``,
397 often performs worse than the rest except when very few eigenvalues are
398 requested for large arrays though there is still no performance guarantee.
401 For the generalized problem, normalization with respect to the given
402 type argument::
404 type 1 and 3 : v.conj().T @ a @ v = w
405 type 2 : inv(v).conj().T @ a @ inv(v) = w
407 type 1 or 2 : v.conj().T @ b @ v = I
408 type 3 : v.conj().T @ inv(b) @ v = I
411 Examples
412 --------
413 >>> import numpy as np
414 >>> from scipy.linalg import eigh
415 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
416 >>> w, v = eigh(A)
417 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
418 True
420 Request only the eigenvalues
422 >>> w = eigh(A, eigvals_only=True)
424 Request eigenvalues that are less than 10.
426 >>> A = np.array([[34, -4, -10, -7, 2],
427 ... [-4, 7, 2, 12, 0],
428 ... [-10, 2, 44, 2, -19],
429 ... [-7, 12, 2, 79, -34],
430 ... [2, 0, -19, -34, 29]])
431 >>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
432 array([6.69199443e-07, 9.11938152e+00])
434 Request the second smallest eigenvalue and its eigenvector
436 >>> w, v = eigh(A, subset_by_index=[1, 1])
437 >>> w
438 array([9.11938152])
439 >>> v.shape # only a single column is returned
440 (5, 1)
442 """
443 if turbo is not _NoValue:
444 warnings.warn("Keyword argument 'turbo' is deprecated in favour of '"
445 "driver=gvd' keyword instead and will be removed in "
446 "SciPy 1.14.0.",
447 DeprecationWarning, stacklevel=2)
448 if eigvals is not _NoValue:
449 warnings.warn("Keyword argument 'eigvals' is deprecated in favour of "
450 "'subset_by_index' keyword instead and will be removed "
451 "in SciPy 1.14.0.",
452 DeprecationWarning, stacklevel=2)
454 # set lower
455 uplo = 'L' if lower else 'U'
456 # Set job for Fortran routines
457 _job = 'N' if eigvals_only else 'V'
459 drv_str = [None, "ev", "evd", "evr", "evx", "gv", "gvd", "gvx"]
460 if driver not in drv_str:
461 raise ValueError('"{}" is unknown. Possible values are "None", "{}".'
462 ''.format(driver, '", "'.join(drv_str[1:])))
464 a1 = _asarray_validated(a, check_finite=check_finite)
465 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
466 raise ValueError('expected square "a" matrix')
467 overwrite_a = overwrite_a or (_datacopied(a1, a))
468 cplx = True if iscomplexobj(a1) else False
469 n = a1.shape[0]
470 drv_args = {'overwrite_a': overwrite_a}
472 if b is not None:
473 b1 = _asarray_validated(b, check_finite=check_finite)
474 overwrite_b = overwrite_b or _datacopied(b1, b)
475 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
476 raise ValueError('expected square "b" matrix')
478 if b1.shape != a1.shape:
479 raise ValueError(f"wrong b dimensions {b1.shape}, should be {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 f'Valid range is [0, {n-1}] and start <= end, but '
506 f'start={lo}, end={hi} is given')
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 f'low={lo}, high={hi} is given')
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(f'{driver} requires input b array to be supplied '
527 'for generalized eigenvalue problems.')
528 if (b is not None) and (driver in ['ev', 'evd', 'evr', 'evx']):
529 raise ValueError(f'"{driver}" does not accept input b array '
530 'for standard eigenvalue problems.')
531 if subset and (driver in ["ev", "evd", "gv", "gvd"]):
532 raise ValueError(f'"{driver}" cannot compute subsets of eigenvalues')
534 # Default driver is evr and gvd
535 else:
536 driver = "evr" if b is None else ("gvx" if subset else "gvd")
538 lwork_spec = {
539 'syevd': ['lwork', 'liwork'],
540 'syevr': ['lwork', 'liwork'],
541 'heevd': ['lwork', 'liwork', 'lrwork'],
542 'heevr': ['lwork', 'lrwork', 'liwork'],
543 }
545 if b is None: # Standard problem
546 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
547 [a1])
548 clw_args = {'n': n, 'lower': lower}
549 if driver == 'evd':
550 clw_args.update({'compute_v': 0 if _job == "N" else 1})
552 lw = _compute_lwork(drvlw, **clw_args)
553 # Multiple lwork vars
554 if isinstance(lw, tuple):
555 lwork_args = dict(zip(lwork_spec[pfx+driver], lw))
556 else:
557 lwork_args = {'lwork': lw}
559 drv_args.update({'lower': lower, 'compute_v': 0 if _job == "N" else 1})
560 w, v, *other_args, info = drv(a=a1, **drv_args, **lwork_args)
562 else: # Generalized problem
563 # 'gvd' doesn't have lwork query
564 if driver == "gvd":
565 drv = get_lapack_funcs(pfx + "gvd", [a1, b1])
566 lwork_args = {}
567 else:
568 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
569 [a1, b1])
570 # generalized drivers use uplo instead of lower
571 lw = _compute_lwork(drvlw, n, uplo=uplo)
572 lwork_args = {'lwork': lw}
574 drv_args.update({'uplo': uplo, 'jobz': _job})
576 w, v, *other_args, info = drv(a=a1, b=b1, **drv_args, **lwork_args)
578 # m is always the first extra argument
579 w = w[:other_args[0]] if subset else w
580 v = v[:, :other_args[0]] if (subset and not eigvals_only) else v
582 # Check if we had a successful exit
583 if info == 0:
584 if eigvals_only:
585 return w
586 else:
587 return w, v
588 else:
589 if info < -1:
590 raise LinAlgError('Illegal value in argument {} of internal {}'
591 ''.format(-info, drv.typecode + pfx + driver))
592 elif info > n:
593 raise LinAlgError(f'The leading minor of order {info-n} of B is not '
594 'positive definite. The factorization of B '
595 'could not be completed and no eigenvalues '
596 'or eigenvectors were computed.')
597 else:
598 drv_err = {'ev': 'The algorithm failed to converge; {} '
599 'off-diagonal elements of an intermediate '
600 'tridiagonal form did not converge to zero.',
601 'evx': '{} eigenvectors failed to converge.',
602 'evd': 'The algorithm failed to compute an eigenvalue '
603 'while working on the submatrix lying in rows '
604 'and columns {0}/{1} through mod({0},{1}).',
605 'evr': 'Internal Error.'
606 }
607 if driver in ['ev', 'gv']:
608 msg = drv_err['ev'].format(info)
609 elif driver in ['evx', 'gvx']:
610 msg = drv_err['evx'].format(info)
611 elif driver in ['evd', 'gvd']:
612 if eigvals_only:
613 msg = drv_err['ev'].format(info)
614 else:
615 msg = drv_err['evd'].format(info, n+1)
616 else:
617 msg = drv_err['evr']
619 raise LinAlgError(msg)
622_conv_dict = {0: 0, 1: 1, 2: 2,
623 'all': 0, 'value': 1, 'index': 2,
624 'a': 0, 'v': 1, 'i': 2}
627def _check_select(select, select_range, max_ev, max_len):
628 """Check that select is valid, convert to Fortran style."""
629 if isinstance(select, str):
630 select = select.lower()
631 try:
632 select = _conv_dict[select]
633 except KeyError as e:
634 raise ValueError('invalid argument for select') from e
635 vl, vu = 0., 1.
636 il = iu = 1
637 if select != 0: # (non-all)
638 sr = asarray(select_range)
639 if sr.ndim != 1 or sr.size != 2 or sr[1] < sr[0]:
640 raise ValueError('select_range must be a 2-element array-like '
641 'in nondecreasing order')
642 if select == 1: # (value)
643 vl, vu = sr
644 if max_ev == 0:
645 max_ev = max_len
646 else: # 2 (index)
647 if sr.dtype.char.lower() not in 'hilqp':
648 raise ValueError(
649 f'when using select="i", select_range must '
650 f'contain integers, got dtype {sr.dtype} ({sr.dtype.char})'
651 )
652 # translate Python (0 ... N-1) into Fortran (1 ... N) with + 1
653 il, iu = sr + 1
654 if min(il, iu) < 1 or max(il, iu) > max_len:
655 raise ValueError('select_range out of bounds')
656 max_ev = iu - il + 1
657 return select, vl, vu, il, iu, max_ev
660def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
661 select='a', select_range=None, max_ev=0, check_finite=True):
662 """
663 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
665 Find eigenvalues w and optionally right eigenvectors v of a::
667 a v[:,i] = w[i] v[:,i]
668 v.H v = identity
670 The matrix a is stored in a_band either in lower diagonal or upper
671 diagonal ordered form:
673 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
674 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
676 where u is the number of bands above the diagonal.
678 Example of a_band (shape of a is (6,6), u=2)::
680 upper form:
681 * * a02 a13 a24 a35
682 * a01 a12 a23 a34 a45
683 a00 a11 a22 a33 a44 a55
685 lower form:
686 a00 a11 a22 a33 a44 a55
687 a10 a21 a32 a43 a54 *
688 a20 a31 a42 a53 * *
690 Cells marked with * are not used.
692 Parameters
693 ----------
694 a_band : (u+1, M) array_like
695 The bands of the M by M matrix a.
696 lower : bool, optional
697 Is the matrix in the lower form. (Default is upper form)
698 eigvals_only : bool, optional
699 Compute only the eigenvalues and no eigenvectors.
700 (Default: calculate also eigenvectors)
701 overwrite_a_band : bool, optional
702 Discard data in a_band (may enhance performance)
703 select : {'a', 'v', 'i'}, optional
704 Which eigenvalues to calculate
706 ====== ========================================
707 select calculated
708 ====== ========================================
709 'a' All eigenvalues
710 'v' Eigenvalues in the interval (min, max]
711 'i' Eigenvalues with indices min <= i <= max
712 ====== ========================================
713 select_range : (min, max), optional
714 Range of selected eigenvalues
715 max_ev : int, optional
716 For select=='v', maximum number of eigenvalues expected.
717 For other values of select, has no meaning.
719 In doubt, leave this parameter untouched.
721 check_finite : bool, optional
722 Whether to check that the input matrix contains only finite numbers.
723 Disabling may give a performance gain, but may result in problems
724 (crashes, non-termination) if the inputs do contain infinities or NaNs.
726 Returns
727 -------
728 w : (M,) ndarray
729 The eigenvalues, in ascending order, each repeated according to its
730 multiplicity.
731 v : (M, M) float or complex ndarray
732 The normalized eigenvector corresponding to the eigenvalue w[i] is
733 the column v[:,i]. Only returned if ``eigvals_only=False``.
735 Raises
736 ------
737 LinAlgError
738 If eigenvalue computation does not converge.
740 See Also
741 --------
742 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
743 eig : eigenvalues and right eigenvectors of general arrays.
744 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
745 eigh_tridiagonal : eigenvalues and right eigenvectors for
746 symmetric/Hermitian tridiagonal matrices
748 Examples
749 --------
750 >>> import numpy as np
751 >>> from scipy.linalg import eig_banded
752 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
753 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
754 >>> w, v = eig_banded(Ab, lower=True)
755 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
756 True
757 >>> w = eig_banded(Ab, lower=True, eigvals_only=True)
758 >>> w
759 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
761 Request only the eigenvalues between ``[-3, 4]``
763 >>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4])
764 >>> w
765 array([-2.22987175, 3.95222349])
767 """
768 if eigvals_only or overwrite_a_band:
769 a1 = _asarray_validated(a_band, check_finite=check_finite)
770 overwrite_a_band = overwrite_a_band or (_datacopied(a1, a_band))
771 else:
772 a1 = array(a_band)
773 if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
774 raise ValueError("array must not contain infs or NaNs")
775 overwrite_a_band = 1
777 if len(a1.shape) != 2:
778 raise ValueError('expected a 2-D array')
779 select, vl, vu, il, iu, max_ev = _check_select(
780 select, select_range, max_ev, a1.shape[1])
781 del select_range
782 if select == 0:
783 if a1.dtype.char in 'GFD':
784 # FIXME: implement this somewhen, for now go with builtin values
785 # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
786 # or by using calc_lwork.f ???
787 # lwork = calc_lwork.hbevd(bevd.typecode, a1.shape[0], lower)
788 internal_name = 'hbevd'
789 else: # a1.dtype.char in 'fd':
790 # FIXME: implement this somewhen, for now go with builtin values
791 # see above
792 # lwork = calc_lwork.sbevd(bevd.typecode, a1.shape[0], lower)
793 internal_name = 'sbevd'
794 bevd, = get_lapack_funcs((internal_name,), (a1,))
795 w, v, info = bevd(a1, compute_v=not eigvals_only,
796 lower=lower, overwrite_ab=overwrite_a_band)
797 else: # select in [1, 2]
798 if eigvals_only:
799 max_ev = 1
800 # calculate optimal abstol for dsbevx (see manpage)
801 if a1.dtype.char in 'fF': # single precision
802 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),))
803 else:
804 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),))
805 abstol = 2 * lamch('s')
806 if a1.dtype.char in 'GFD':
807 internal_name = 'hbevx'
808 else: # a1.dtype.char in 'gfd'
809 internal_name = 'sbevx'
810 bevx, = get_lapack_funcs((internal_name,), (a1,))
811 w, v, m, ifail, info = bevx(
812 a1, vl, vu, il, iu, compute_v=not eigvals_only, mmax=max_ev,
813 range=select, lower=lower, overwrite_ab=overwrite_a_band,
814 abstol=abstol)
815 # crop off w and v
816 w = w[:m]
817 if not eigvals_only:
818 v = v[:, :m]
819 _check_info(info, internal_name)
821 if eigvals_only:
822 return w
823 return w, v
826def eigvals(a, b=None, overwrite_a=False, check_finite=True,
827 homogeneous_eigvals=False):
828 """
829 Compute eigenvalues from an ordinary or generalized eigenvalue problem.
831 Find eigenvalues of a general matrix::
833 a vr[:,i] = w[i] b vr[:,i]
835 Parameters
836 ----------
837 a : (M, M) array_like
838 A complex or real matrix whose eigenvalues and eigenvectors
839 will be computed.
840 b : (M, M) array_like, optional
841 Right-hand side matrix in a generalized eigenvalue problem.
842 If omitted, identity matrix is assumed.
843 overwrite_a : bool, optional
844 Whether to overwrite data in a (may improve performance)
845 check_finite : bool, optional
846 Whether to check that the input matrices contain only finite numbers.
847 Disabling may give a performance gain, but may result in problems
848 (crashes, non-termination) if the inputs do contain infinities
849 or NaNs.
850 homogeneous_eigvals : bool, optional
851 If True, return the eigenvalues in homogeneous coordinates.
852 In this case ``w`` is a (2, M) array so that::
854 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
856 Default is False.
858 Returns
859 -------
860 w : (M,) or (2, M) double or complex ndarray
861 The eigenvalues, each repeated according to its multiplicity
862 but not in any specific order. The shape is (M,) unless
863 ``homogeneous_eigvals=True``.
865 Raises
866 ------
867 LinAlgError
868 If eigenvalue computation does not converge
870 See Also
871 --------
872 eig : eigenvalues and right eigenvectors of general arrays.
873 eigvalsh : eigenvalues of symmetric or Hermitian arrays
874 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
875 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
876 matrices
878 Examples
879 --------
880 >>> import numpy as np
881 >>> from scipy import linalg
882 >>> a = np.array([[0., -1.], [1., 0.]])
883 >>> linalg.eigvals(a)
884 array([0.+1.j, 0.-1.j])
886 >>> b = np.array([[0., 1.], [1., 1.]])
887 >>> linalg.eigvals(a, b)
888 array([ 1.+0.j, -1.+0.j])
890 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
891 >>> linalg.eigvals(a, homogeneous_eigvals=True)
892 array([[3.+0.j, 8.+0.j, 7.+0.j],
893 [1.+0.j, 1.+0.j, 1.+0.j]])
895 """
896 return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a,
897 check_finite=check_finite,
898 homogeneous_eigvals=homogeneous_eigvals)
901@_deprecate_positional_args(version="1.14.0")
902def eigvalsh(a, b=None, *, lower=True, overwrite_a=False,
903 overwrite_b=False, turbo=_NoValue, eigvals=_NoValue, type=1,
904 check_finite=True, subset_by_index=None, subset_by_value=None,
905 driver=None):
906 """
907 Solves a standard or generalized eigenvalue problem for a complex
908 Hermitian or real symmetric matrix.
910 Find eigenvalues array ``w`` of array ``a``, where ``b`` is positive
911 definite such that for every eigenvalue λ (i-th entry of w) and its
912 eigenvector vi (i-th column of v) satisfies::
914 a @ vi = λ * b @ vi
915 vi.conj().T @ a @ vi = λ
916 vi.conj().T @ b @ vi = 1
918 In the standard problem, b is assumed to be the identity matrix.
920 Parameters
921 ----------
922 a : (M, M) array_like
923 A complex Hermitian or real symmetric matrix whose eigenvalues will
924 be computed.
925 b : (M, M) array_like, optional
926 A complex Hermitian or real symmetric definite positive matrix in.
927 If omitted, identity matrix is assumed.
928 lower : bool, optional
929 Whether the pertinent array data is taken from the lower or upper
930 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
931 overwrite_a : bool, optional
932 Whether to overwrite data in ``a`` (may improve performance). Default
933 is False.
934 overwrite_b : bool, optional
935 Whether to overwrite data in ``b`` (may improve performance). Default
936 is False.
937 type : int, optional
938 For the generalized problems, this keyword specifies the problem type
939 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
940 inputs)::
942 1 => a @ v = w @ b @ v
943 2 => a @ b @ v = w @ v
944 3 => b @ a @ v = w @ v
946 This keyword is ignored for standard problems.
947 check_finite : bool, optional
948 Whether to check that the input matrices contain only finite numbers.
949 Disabling may give a performance gain, but may result in problems
950 (crashes, non-termination) if the inputs do contain infinities or NaNs.
951 subset_by_index : iterable, optional
952 If provided, this two-element iterable defines the start and the end
953 indices of the desired eigenvalues (ascending order and 0-indexed).
954 To return only the second smallest to fifth smallest eigenvalues,
955 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
956 available with "evr", "evx", and "gvx" drivers. The entries are
957 directly converted to integers via ``int()``.
958 subset_by_value : iterable, optional
959 If provided, this two-element iterable defines the half-open interval
960 ``(a, b]`` that, if any, only the eigenvalues between these values
961 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
962 ``np.inf`` for the unconstrained ends.
963 driver : str, optional
964 Defines which LAPACK driver should be used. Valid options are "ev",
965 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
966 generalized (where b is not None) problems. See the Notes section of
967 `scipy.linalg.eigh`.
968 turbo : bool, optional, deprecated
969 .. deprecated:: 1.5.0
970 'eigvalsh' keyword argument `turbo` is deprecated in favor of
971 ``driver=gvd`` option and will be removed in SciPy 1.14.0.
973 eigvals : tuple (lo, hi), optional
974 .. deprecated:: 1.5.0
975 'eigvalsh' keyword argument `eigvals` is deprecated in favor of
976 `subset_by_index` option and will be removed in SciPy 1.14.0.
978 Returns
979 -------
980 w : (N,) ndarray
981 The N (N<=M) selected eigenvalues, in ascending order, each
982 repeated according to its multiplicity.
984 Raises
985 ------
986 LinAlgError
987 If eigenvalue computation does not converge, an error occurred, or
988 b matrix is not definite positive. Note that if input matrices are
989 not symmetric or Hermitian, no error will be reported but results will
990 be wrong.
992 See Also
993 --------
994 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
995 eigvals : eigenvalues of general arrays
996 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
997 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
998 matrices
1000 Notes
1001 -----
1002 This function does not check the input array for being Hermitian/symmetric
1003 in order to allow for representing arrays with only their upper/lower
1004 triangular parts.
1006 This function serves as a one-liner shorthand for `scipy.linalg.eigh` with
1007 the option ``eigvals_only=True`` to get the eigenvalues and not the
1008 eigenvectors. Here it is kept as a legacy convenience. It might be
1009 beneficial to use the main function to have full control and to be a bit
1010 more pythonic.
1012 Examples
1013 --------
1014 For more examples see `scipy.linalg.eigh`.
1016 >>> import numpy as np
1017 >>> from scipy.linalg import eigvalsh
1018 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
1019 >>> w = eigvalsh(A)
1020 >>> w
1021 array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
1023 """
1024 return eigh(a, b=b, lower=lower, eigvals_only=True,
1025 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
1026 turbo=turbo, eigvals=eigvals, type=type,
1027 check_finite=check_finite, subset_by_index=subset_by_index,
1028 subset_by_value=subset_by_value, driver=driver)
1031def eigvals_banded(a_band, lower=False, overwrite_a_band=False,
1032 select='a', select_range=None, check_finite=True):
1033 """
1034 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
1036 Find eigenvalues w of a::
1038 a v[:,i] = w[i] v[:,i]
1039 v.H v = identity
1041 The matrix a is stored in a_band either in lower diagonal or upper
1042 diagonal ordered form:
1044 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
1045 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
1047 where u is the number of bands above the diagonal.
1049 Example of a_band (shape of a is (6,6), u=2)::
1051 upper form:
1052 * * a02 a13 a24 a35
1053 * a01 a12 a23 a34 a45
1054 a00 a11 a22 a33 a44 a55
1056 lower form:
1057 a00 a11 a22 a33 a44 a55
1058 a10 a21 a32 a43 a54 *
1059 a20 a31 a42 a53 * *
1061 Cells marked with * are not used.
1063 Parameters
1064 ----------
1065 a_band : (u+1, M) array_like
1066 The bands of the M by M matrix a.
1067 lower : bool, optional
1068 Is the matrix in the lower form. (Default is upper form)
1069 overwrite_a_band : bool, optional
1070 Discard data in a_band (may enhance performance)
1071 select : {'a', 'v', 'i'}, optional
1072 Which eigenvalues to calculate
1074 ====== ========================================
1075 select calculated
1076 ====== ========================================
1077 'a' All eigenvalues
1078 'v' Eigenvalues in the interval (min, max]
1079 'i' Eigenvalues with indices min <= i <= max
1080 ====== ========================================
1081 select_range : (min, max), optional
1082 Range of selected eigenvalues
1083 check_finite : bool, optional
1084 Whether to check that the input matrix contains only finite numbers.
1085 Disabling may give a performance gain, but may result in problems
1086 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1088 Returns
1089 -------
1090 w : (M,) ndarray
1091 The eigenvalues, in ascending order, each repeated according to its
1092 multiplicity.
1094 Raises
1095 ------
1096 LinAlgError
1097 If eigenvalue computation does not converge.
1099 See Also
1100 --------
1101 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1102 band matrices
1103 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1104 matrices
1105 eigvals : eigenvalues of general arrays
1106 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1107 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1109 Examples
1110 --------
1111 >>> import numpy as np
1112 >>> from scipy.linalg import eigvals_banded
1113 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
1114 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
1115 >>> w = eigvals_banded(Ab, lower=True)
1116 >>> w
1117 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
1118 """
1119 return eig_banded(a_band, lower=lower, eigvals_only=1,
1120 overwrite_a_band=overwrite_a_band, select=select,
1121 select_range=select_range, check_finite=check_finite)
1124def eigvalsh_tridiagonal(d, e, select='a', select_range=None,
1125 check_finite=True, tol=0., lapack_driver='auto'):
1126 """
1127 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1129 Find eigenvalues `w` of ``a``::
1131 a v[:,i] = w[i] v[:,i]
1132 v.H v = identity
1134 For a real symmetric matrix ``a`` with diagonal elements `d` and
1135 off-diagonal elements `e`.
1137 Parameters
1138 ----------
1139 d : ndarray, shape (ndim,)
1140 The diagonal elements of the array.
1141 e : ndarray, shape (ndim-1,)
1142 The off-diagonal elements of the array.
1143 select : {'a', 'v', 'i'}, optional
1144 Which eigenvalues to calculate
1146 ====== ========================================
1147 select calculated
1148 ====== ========================================
1149 'a' All eigenvalues
1150 'v' Eigenvalues in the interval (min, max]
1151 'i' Eigenvalues with indices min <= i <= max
1152 ====== ========================================
1153 select_range : (min, max), optional
1154 Range of selected eigenvalues
1155 check_finite : bool, optional
1156 Whether to check that the input matrix contains only finite numbers.
1157 Disabling may give a performance gain, but may result in problems
1158 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1159 tol : float
1160 The absolute tolerance to which each eigenvalue is required
1161 (only used when ``lapack_driver='stebz'``).
1162 An eigenvalue (or cluster) is considered to have converged if it
1163 lies in an interval of this width. If <= 0. (default),
1164 the value ``eps*|a|`` is used where eps is the machine precision,
1165 and ``|a|`` is the 1-norm of the matrix ``a``.
1166 lapack_driver : str
1167 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1168 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1169 and 'stebz' otherwise. 'sterf' and 'stev' can only be used when
1170 ``select='a'``.
1172 Returns
1173 -------
1174 w : (M,) ndarray
1175 The eigenvalues, in ascending order, each repeated according to its
1176 multiplicity.
1178 Raises
1179 ------
1180 LinAlgError
1181 If eigenvalue computation does not converge.
1183 See Also
1184 --------
1185 eigh_tridiagonal : eigenvalues and right eiegenvectors for
1186 symmetric/Hermitian tridiagonal matrices
1188 Examples
1189 --------
1190 >>> import numpy as np
1191 >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
1192 >>> d = 3*np.ones(4)
1193 >>> e = -1*np.ones(3)
1194 >>> w = eigvalsh_tridiagonal(d, e)
1195 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1196 >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
1197 >>> np.allclose(w - w2, np.zeros(4))
1198 True
1199 """
1200 return eigh_tridiagonal(
1201 d, e, eigvals_only=True, select=select, select_range=select_range,
1202 check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
1205def eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None,
1206 check_finite=True, tol=0., lapack_driver='auto'):
1207 """
1208 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1210 Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
1212 a v[:,i] = w[i] v[:,i]
1213 v.H v = identity
1215 For a real symmetric matrix ``a`` with diagonal elements `d` and
1216 off-diagonal elements `e`.
1218 Parameters
1219 ----------
1220 d : ndarray, shape (ndim,)
1221 The diagonal elements of the array.
1222 e : ndarray, shape (ndim-1,)
1223 The off-diagonal elements of the array.
1224 eigvals_only : bool, optional
1225 Compute only the eigenvalues and no eigenvectors.
1226 (Default: calculate also eigenvectors)
1227 select : {'a', 'v', 'i'}, optional
1228 Which eigenvalues to calculate
1230 ====== ========================================
1231 select calculated
1232 ====== ========================================
1233 'a' All eigenvalues
1234 'v' Eigenvalues in the interval (min, max]
1235 'i' Eigenvalues with indices min <= i <= max
1236 ====== ========================================
1237 select_range : (min, max), optional
1238 Range of selected eigenvalues
1239 check_finite : bool, optional
1240 Whether to check that the input matrix contains only finite numbers.
1241 Disabling may give a performance gain, but may result in problems
1242 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1243 tol : float
1244 The absolute tolerance to which each eigenvalue is required
1245 (only used when 'stebz' is the `lapack_driver`).
1246 An eigenvalue (or cluster) is considered to have converged if it
1247 lies in an interval of this width. If <= 0. (default),
1248 the value ``eps*|a|`` is used where eps is the machine precision,
1249 and ``|a|`` is the 1-norm of the matrix ``a``.
1250 lapack_driver : str
1251 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1252 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1253 and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and
1254 ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is
1255 used to find the corresponding eigenvectors. 'sterf' can only be
1256 used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only
1257 be used when ``select='a'``.
1259 Returns
1260 -------
1261 w : (M,) ndarray
1262 The eigenvalues, in ascending order, each repeated according to its
1263 multiplicity.
1264 v : (M, M) ndarray
1265 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
1266 the column ``v[:,i]``. Only returned if ``eigvals_only=False``.
1268 Raises
1269 ------
1270 LinAlgError
1271 If eigenvalue computation does not converge.
1273 See Also
1274 --------
1275 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1276 matrices
1277 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1278 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1279 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1280 band matrices
1282 Notes
1283 -----
1284 This function makes use of LAPACK ``S/DSTEMR`` routines.
1286 Examples
1287 --------
1288 >>> import numpy as np
1289 >>> from scipy.linalg import eigh_tridiagonal
1290 >>> d = 3*np.ones(4)
1291 >>> e = -1*np.ones(3)
1292 >>> w, v = eigh_tridiagonal(d, e)
1293 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1294 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
1295 True
1296 """
1297 d = _asarray_validated(d, check_finite=check_finite)
1298 e = _asarray_validated(e, check_finite=check_finite)
1299 for check in (d, e):
1300 if check.ndim != 1:
1301 raise ValueError('expected a 1-D array')
1302 if check.dtype.char in 'GFD': # complex
1303 raise TypeError('Only real arrays currently supported')
1304 if d.size != e.size + 1:
1305 raise ValueError(f'd ({d.size}) must have one more element than e ({e.size})')
1306 select, vl, vu, il, iu, _ = _check_select(
1307 select, select_range, 0, d.size)
1308 if not isinstance(lapack_driver, str):
1309 raise TypeError('lapack_driver must be str')
1310 drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev')
1311 if lapack_driver not in drivers:
1312 raise ValueError(f'lapack_driver must be one of {drivers}, '
1313 f'got {lapack_driver}')
1314 if lapack_driver == 'auto':
1315 lapack_driver = 'stemr' if select == 0 else 'stebz'
1317 # Quick exit for 1x1 case
1318 if len(d) == 1:
1319 if select == 1 and (not (vl < d[0] <= vu)): # request by value
1320 w = array([])
1321 v = empty([1, 0], dtype=d.dtype)
1322 else: # all and request by index
1323 w = array([d[0]], dtype=d.dtype)
1324 v = array([[1.]], dtype=d.dtype)
1326 if eigvals_only:
1327 return w
1328 else:
1329 return w, v
1331 func, = get_lapack_funcs((lapack_driver,), (d, e))
1332 compute_v = not eigvals_only
1333 if lapack_driver == 'sterf':
1334 if select != 0:
1335 raise ValueError('sterf can only be used when select == "a"')
1336 if not eigvals_only:
1337 raise ValueError('sterf can only be used when eigvals_only is '
1338 'True')
1339 w, info = func(d, e)
1340 m = len(w)
1341 elif lapack_driver == 'stev':
1342 if select != 0:
1343 raise ValueError('stev can only be used when select == "a"')
1344 w, v, info = func(d, e, compute_v=compute_v)
1345 m = len(w)
1346 elif lapack_driver == 'stebz':
1347 tol = float(tol)
1348 internal_name = 'stebz'
1349 stebz, = get_lapack_funcs((internal_name,), (d, e))
1350 # If getting eigenvectors, needs to be block-ordered (B) instead of
1351 # matrix-ordered (E), and we will reorder later
1352 order = 'E' if eigvals_only else 'B'
1353 m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol,
1354 order)
1355 else: # 'stemr'
1356 # ?STEMR annoyingly requires size N instead of N-1
1357 e_ = empty(e.size+1, e.dtype)
1358 e_[:-1] = e
1359 stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e))
1360 lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu,
1361 compute_v=compute_v)
1362 _check_info(info, 'stemr_lwork')
1363 m, w, v, info = func(d, e_, select, vl, vu, il, iu,
1364 compute_v=compute_v, lwork=lwork, liwork=liwork)
1365 _check_info(info, lapack_driver + ' (eigh_tridiagonal)')
1366 w = w[:m]
1367 if eigvals_only:
1368 return w
1369 else:
1370 # Do we still need to compute the eigenvalues?
1371 if lapack_driver == 'stebz':
1372 func, = get_lapack_funcs(('stein',), (d, e))
1373 v, info = func(d, e, w, iblock, isplit)
1374 _check_info(info, 'stein (eigh_tridiagonal)',
1375 positive='%d eigenvectors failed to converge')
1376 # Convert block-order to matrix-order
1377 order = argsort(w)
1378 w, v = w[order], v[:, order]
1379 else:
1380 v = v[:, :m]
1381 return w, v
1384def _check_info(info, driver, positive='did not converge (LAPACK info=%d)'):
1385 """Check info return value."""
1386 if info < 0:
1387 raise ValueError('illegal value in argument %d of internal %s'
1388 % (-info, driver))
1389 if info > 0 and positive:
1390 raise LinAlgError(("%s " + positive) % (driver, info,))
1393def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
1394 """
1395 Compute Hessenberg form of a matrix.
1397 The Hessenberg decomposition is::
1399 A = Q H Q^H
1401 where `Q` is unitary/orthogonal and `H` has only zero elements below
1402 the first sub-diagonal.
1404 Parameters
1405 ----------
1406 a : (M, M) array_like
1407 Matrix to bring into Hessenberg form.
1408 calc_q : bool, optional
1409 Whether to compute the transformation matrix. Default is False.
1410 overwrite_a : bool, optional
1411 Whether to overwrite `a`; may improve performance.
1412 Default is False.
1413 check_finite : bool, optional
1414 Whether to check that the input matrix contains only finite numbers.
1415 Disabling may give a performance gain, but may result in problems
1416 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1418 Returns
1419 -------
1420 H : (M, M) ndarray
1421 Hessenberg form of `a`.
1422 Q : (M, M) ndarray
1423 Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
1424 Only returned if ``calc_q=True``.
1426 Examples
1427 --------
1428 >>> import numpy as np
1429 >>> from scipy.linalg import hessenberg
1430 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
1431 >>> H, Q = hessenberg(A, calc_q=True)
1432 >>> H
1433 array([[ 2. , -11.65843866, 1.42005301, 0.25349066],
1434 [ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
1435 [ 0. , -1.83299243, 0.38969961, -0.51527034],
1436 [ 0. , 0. , -3.83189513, 1.07494686]])
1437 >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
1438 True
1439 """
1440 a1 = _asarray_validated(a, check_finite=check_finite)
1441 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
1442 raise ValueError('expected square matrix')
1443 overwrite_a = overwrite_a or (_datacopied(a1, a))
1445 # if 2x2 or smaller: already in Hessenberg
1446 if a1.shape[0] <= 2:
1447 if calc_q:
1448 return a1, eye(a1.shape[0])
1449 return a1
1451 gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal',
1452 'gehrd_lwork'), (a1,))
1453 ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a)
1454 _check_info(info, 'gebal (hessenberg)', positive=False)
1455 n = len(a1)
1457 lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
1459 hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1460 _check_info(info, 'gehrd (hessenberg)', positive=False)
1461 h = numpy.triu(hq, -1)
1462 if not calc_q:
1463 return h
1465 # use orghr/unghr to compute q
1466 orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,))
1467 lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
1469 q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1470 _check_info(info, 'orghr (hessenberg)', positive=False)
1471 return h, q
1474def cdf2rdf(w, v):
1475 """
1476 Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real
1477 eigenvalues in a block diagonal form ``wr`` and the associated real
1478 eigenvectors ``vr``, such that::
1480 vr @ wr = X @ vr
1482 continues to hold, where ``X`` is the original array for which ``w`` and
1483 ``v`` are the eigenvalues and eigenvectors.
1485 .. versionadded:: 1.1.0
1487 Parameters
1488 ----------
1489 w : (..., M) array_like
1490 Complex or real eigenvalues, an array or stack of arrays
1492 Conjugate pairs must not be interleaved, else the wrong result
1493 will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result,
1494 but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
1496 v : (..., M, M) array_like
1497 Complex or real eigenvectors, a square array or stack of square arrays.
1499 Returns
1500 -------
1501 wr : (..., M, M) ndarray
1502 Real diagonal block form of eigenvalues
1503 vr : (..., M, M) ndarray
1504 Real eigenvectors associated with ``wr``
1506 See Also
1507 --------
1508 eig : Eigenvalues and right eigenvectors for non-symmetric arrays
1509 rsf2csf : Convert real Schur form to complex Schur form
1511 Notes
1512 -----
1513 ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``.
1514 For example, obtained by ``w, v = scipy.linalg.eig(X)`` or
1515 ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent
1516 stacked arrays.
1518 .. versionadded:: 1.1.0
1520 Examples
1521 --------
1522 >>> import numpy as np
1523 >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
1524 >>> X
1525 array([[ 1, 2, 3],
1526 [ 0, 4, 5],
1527 [ 0, -5, 4]])
1529 >>> from scipy import linalg
1530 >>> w, v = linalg.eig(X)
1531 >>> w
1532 array([ 1.+0.j, 4.+5.j, 4.-5.j])
1533 >>> v
1534 array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
1535 [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
1536 [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
1538 >>> wr, vr = linalg.cdf2rdf(w, v)
1539 >>> wr
1540 array([[ 1., 0., 0.],
1541 [ 0., 4., 5.],
1542 [ 0., -5., 4.]])
1543 >>> vr
1544 array([[ 1. , 0.40016, -0.01906],
1545 [ 0. , 0.64788, 0. ],
1546 [ 0. , 0. , 0.64788]])
1548 >>> vr @ wr
1549 array([[ 1. , 1.69593, 1.9246 ],
1550 [ 0. , 2.59153, 3.23942],
1551 [ 0. , -3.23942, 2.59153]])
1552 >>> X @ vr
1553 array([[ 1. , 1.69593, 1.9246 ],
1554 [ 0. , 2.59153, 3.23942],
1555 [ 0. , -3.23942, 2.59153]])
1556 """
1557 w, v = _asarray_validated(w), _asarray_validated(v)
1559 # check dimensions
1560 if w.ndim < 1:
1561 raise ValueError('expected w to be at least 1D')
1562 if v.ndim < 2:
1563 raise ValueError('expected v to be at least 2D')
1564 if v.ndim != w.ndim + 1:
1565 raise ValueError('expected eigenvectors array to have exactly one '
1566 'dimension more than eigenvalues array')
1568 # check shapes
1569 n = w.shape[-1]
1570 M = w.shape[:-1]
1571 if v.shape[-2] != v.shape[-1]:
1572 raise ValueError('expected v to be a square matrix or stacked square '
1573 'matrices: v.shape[-2] = v.shape[-1]')
1574 if v.shape[-1] != n:
1575 raise ValueError('expected the same number of eigenvalues as '
1576 'eigenvectors')
1578 # get indices for each first pair of complex eigenvalues
1579 complex_mask = iscomplex(w)
1580 n_complex = complex_mask.sum(axis=-1)
1582 # check if all complex eigenvalues have conjugate pairs
1583 if not (n_complex % 2 == 0).all():
1584 raise ValueError('expected complex-conjugate pairs of eigenvalues')
1586 # find complex indices
1587 idx = nonzero(complex_mask)
1588 idx_stack = idx[:-1]
1589 idx_elem = idx[-1]
1591 # filter them to conjugate indices, assuming pairs are not interleaved
1592 j = idx_elem[0::2]
1593 k = idx_elem[1::2]
1594 stack_ind = ()
1595 for i in idx_stack:
1596 # should never happen, assuming nonzero orders by the last axis
1597 assert (i[0::2] == i[1::2]).all(), \
1598 "Conjugate pair spanned different arrays!"
1599 stack_ind += (i[0::2],)
1601 # all eigenvalues to diagonal form
1602 wr = zeros(M + (n, n), dtype=w.real.dtype)
1603 di = range(n)
1604 wr[..., di, di] = w.real
1606 # complex eigenvalues to real block diagonal form
1607 wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag
1608 wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
1610 # compute real eigenvectors associated with real block diagonal eigenvalues
1611 u = zeros(M + (n, n), dtype=numpy.cdouble)
1612 u[..., di, di] = 1.0
1613 u[stack_ind + (j, j)] = 0.5j
1614 u[stack_ind + (j, k)] = 0.5
1615 u[stack_ind + (k, j)] = -0.5j
1616 u[stack_ind + (k, k)] = 0.5
1618 # multiply matrices v and u (equivalent to v @ u)
1619 vr = einsum('...ij,...jk->...ik', v, u).real
1621 return wr, vr