Coverage for /usr/lib/python3/dist-packages/scipy/linalg/_decomp.py: 6%
390 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-06-14 15:25 +0200
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, cast,
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
31_I = cast['F'](1j)
34def _make_complex_eigvecs(w, vin, dtype):
35 """
36 Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output
37 """
38 # - see LAPACK man page DGGEV at ALPHAI
39 v = numpy.array(vin, dtype=dtype)
40 m = (w.imag > 0)
41 m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709
42 for i in flatnonzero(m):
43 v.imag[:, i] = vin[:, i+1]
44 conj(v[:, i], v[:, i+1])
45 return v
48def _make_eigvals(alpha, beta, homogeneous_eigvals):
49 if homogeneous_eigvals:
50 if beta is None:
51 return numpy.vstack((alpha, numpy.ones_like(alpha)))
52 else:
53 return numpy.vstack((alpha, beta))
54 else:
55 if beta is None:
56 return alpha
57 else:
58 w = numpy.empty_like(alpha)
59 alpha_zero = (alpha == 0)
60 beta_zero = (beta == 0)
61 beta_nonzero = ~beta_zero
62 w[beta_nonzero] = alpha[beta_nonzero]/beta[beta_nonzero]
63 # Use numpy.inf for complex values too since
64 # 1/numpy.inf = 0, i.e., it correctly behaves as projective
65 # infinity.
66 w[~alpha_zero & beta_zero] = numpy.inf
67 if numpy.all(alpha.imag == 0):
68 w[alpha_zero & beta_zero] = numpy.nan
69 else:
70 w[alpha_zero & beta_zero] = complex(numpy.nan, numpy.nan)
71 return w
74def _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
75 homogeneous_eigvals):
76 ggev, = get_lapack_funcs(('ggev',), (a1, b1))
77 cvl, cvr = left, right
78 res = ggev(a1, b1, lwork=-1)
79 lwork = res[-2][0].real.astype(numpy.int_)
80 if ggev.typecode in 'cz':
81 alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
82 overwrite_a, overwrite_b)
83 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
84 else:
85 alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr,
86 lwork, overwrite_a,
87 overwrite_b)
88 alpha = alphar + _I * alphai
89 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
90 _check_info(info, 'generalized eig algorithm (ggev)')
92 only_real = numpy.all(w.imag == 0.0)
93 if not (ggev.typecode in 'cz' or only_real):
94 t = w.dtype.char
95 if left:
96 vl = _make_complex_eigvecs(w, vl, t)
97 if right:
98 vr = _make_complex_eigvecs(w, vr, t)
100 # the eigenvectors returned by the lapack function are NOT normalized
101 for i in range(vr.shape[0]):
102 if right:
103 vr[:, i] /= norm(vr[:, i])
104 if left:
105 vl[:, i] /= norm(vl[:, i])
107 if not (left or right):
108 return w
109 if left:
110 if right:
111 return w, vl, vr
112 return w, vl
113 return w, vr
116def eig(a, b=None, left=False, right=True, overwrite_a=False,
117 overwrite_b=False, check_finite=True, homogeneous_eigvals=False):
118 """
119 Solve an ordinary or generalized eigenvalue problem of a square matrix.
121 Find eigenvalues w and right or left eigenvectors of a general matrix::
123 a vr[:,i] = w[i] b vr[:,i]
124 a.H vl[:,i] = w[i].conj() b.H vl[:,i]
126 where ``.H`` is the Hermitian conjugation.
128 Parameters
129 ----------
130 a : (M, M) array_like
131 A complex or real matrix whose eigenvalues and eigenvectors
132 will be computed.
133 b : (M, M) array_like, optional
134 Right-hand side matrix in a generalized eigenvalue problem.
135 Default is None, identity matrix is assumed.
136 left : bool, optional
137 Whether to calculate and return left eigenvectors. Default is False.
138 right : bool, optional
139 Whether to calculate and return right eigenvectors. Default is True.
140 overwrite_a : bool, optional
141 Whether to overwrite `a`; may improve performance. Default is False.
142 overwrite_b : bool, optional
143 Whether to overwrite `b`; may improve performance. Default is False.
144 check_finite : bool, optional
145 Whether to check that the input matrices contain only finite numbers.
146 Disabling may give a performance gain, but may result in problems
147 (crashes, non-termination) if the inputs do contain infinities or NaNs.
148 homogeneous_eigvals : bool, optional
149 If True, return the eigenvalues in homogeneous coordinates.
150 In this case ``w`` is a (2, M) array so that::
152 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
154 Default is False.
156 Returns
157 -------
158 w : (M,) or (2, M) double or complex ndarray
159 The eigenvalues, each repeated according to its
160 multiplicity. The shape is (M,) unless
161 ``homogeneous_eigvals=True``.
162 vl : (M, M) double or complex ndarray
163 The normalized left eigenvector corresponding to the eigenvalue
164 ``w[i]`` is the column vl[:,i]. Only returned if ``left=True``.
165 vr : (M, M) double or complex ndarray
166 The normalized right eigenvector corresponding to the eigenvalue
167 ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
169 Raises
170 ------
171 LinAlgError
172 If eigenvalue computation does not converge.
174 See Also
175 --------
176 eigvals : eigenvalues of general arrays
177 eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
178 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
179 band matrices
180 eigh_tridiagonal : eigenvalues and right eiegenvectors for
181 symmetric/Hermitian tridiagonal matrices
183 Examples
184 --------
185 >>> import numpy as np
186 >>> from scipy import linalg
187 >>> a = np.array([[0., -1.], [1., 0.]])
188 >>> linalg.eigvals(a)
189 array([0.+1.j, 0.-1.j])
191 >>> b = np.array([[0., 1.], [1., 1.]])
192 >>> linalg.eigvals(a, b)
193 array([ 1.+0.j, -1.+0.j])
195 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
196 >>> linalg.eigvals(a, homogeneous_eigvals=True)
197 array([[3.+0.j, 8.+0.j, 7.+0.j],
198 [1.+0.j, 1.+0.j, 1.+0.j]])
200 >>> a = np.array([[0., -1.], [1., 0.]])
201 >>> linalg.eigvals(a) == linalg.eig(a)[0]
202 array([ True, True])
203 >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector
204 array([[-0.70710678+0.j , -0.70710678-0.j ],
205 [-0. +0.70710678j, -0. -0.70710678j]])
206 >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector
207 array([[0.70710678+0.j , 0.70710678-0.j ],
208 [0. -0.70710678j, 0. +0.70710678j]])
212 """
213 a1 = _asarray_validated(a, check_finite=check_finite)
214 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
215 raise ValueError('expected square matrix')
216 overwrite_a = overwrite_a or (_datacopied(a1, a))
217 if b is not None:
218 b1 = _asarray_validated(b, check_finite=check_finite)
219 overwrite_b = overwrite_b or _datacopied(b1, b)
220 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
221 raise ValueError('expected square matrix')
222 if b1.shape != a1.shape:
223 raise ValueError('a and b must have the same shape')
224 return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
225 homogeneous_eigvals)
227 geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
228 compute_vl, compute_vr = left, right
230 lwork = _compute_lwork(geev_lwork, a1.shape[0],
231 compute_vl=compute_vl,
232 compute_vr=compute_vr)
234 if geev.typecode in 'cz':
235 w, vl, vr, info = geev(a1, lwork=lwork,
236 compute_vl=compute_vl,
237 compute_vr=compute_vr,
238 overwrite_a=overwrite_a)
239 w = _make_eigvals(w, None, homogeneous_eigvals)
240 else:
241 wr, wi, vl, vr, info = geev(a1, lwork=lwork,
242 compute_vl=compute_vl,
243 compute_vr=compute_vr,
244 overwrite_a=overwrite_a)
245 t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
246 w = wr + _I * wi
247 w = _make_eigvals(w, None, homogeneous_eigvals)
249 _check_info(info, 'eig algorithm (geev)',
250 positive='did not converge (only eigenvalues '
251 'with order >= %d have converged)')
253 only_real = numpy.all(w.imag == 0.0)
254 if not (geev.typecode in 'cz' or only_real):
255 t = w.dtype.char
256 if left:
257 vl = _make_complex_eigvecs(w, vl, t)
258 if right:
259 vr = _make_complex_eigvecs(w, vr, t)
260 if not (left or right):
261 return w
262 if left:
263 if right:
264 return w, vl, vr
265 return w, vl
266 return w, vr
269def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
270 overwrite_b=False, turbo=False, eigvals=None, type=1,
271 check_finite=True, subset_by_index=None, subset_by_value=None,
272 driver=None):
273 """
274 Solve a standard or generalized eigenvalue problem for a complex
275 Hermitian or real symmetric matrix.
277 Find eigenvalues array ``w`` and optionally eigenvectors array ``v`` of
278 array ``a``, where ``b`` is positive definite such that for every
279 eigenvalue λ (i-th entry of w) and its eigenvector ``vi`` (i-th column of
280 ``v``) satisfies::
282 a @ vi = λ * b @ vi
283 vi.conj().T @ a @ vi = λ
284 vi.conj().T @ b @ vi = 1
286 In the standard problem, ``b`` is assumed to be the identity matrix.
288 Parameters
289 ----------
290 a : (M, M) array_like
291 A complex Hermitian or real symmetric matrix whose eigenvalues and
292 eigenvectors will be computed.
293 b : (M, M) array_like, optional
294 A complex Hermitian or real symmetric definite positive matrix in.
295 If omitted, identity matrix is assumed.
296 lower : bool, optional
297 Whether the pertinent array data is taken from the lower or upper
298 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
299 eigvals_only : bool, optional
300 Whether to calculate only eigenvalues and no eigenvectors.
301 (Default: both are calculated)
302 subset_by_index : iterable, optional
303 If provided, this two-element iterable defines the start and the end
304 indices of the desired eigenvalues (ascending order and 0-indexed).
305 To return only the second smallest to fifth smallest eigenvalues,
306 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
307 available with "evr", "evx", and "gvx" drivers. The entries are
308 directly converted to integers via ``int()``.
309 subset_by_value : iterable, optional
310 If provided, this two-element iterable defines the half-open interval
311 ``(a, b]`` that, if any, only the eigenvalues between these values
312 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
313 ``np.inf`` for the unconstrained ends.
314 driver : str, optional
315 Defines which LAPACK driver should be used. Valid options are "ev",
316 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
317 generalized (where b is not None) problems. See the Notes section.
318 The default for standard problems is "evr". For generalized problems,
319 "gvd" is used for full set, and "gvx" for subset requested cases.
320 type : int, optional
321 For the generalized problems, this keyword specifies the problem type
322 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
323 inputs)::
325 1 => a @ v = w @ b @ v
326 2 => a @ b @ v = w @ v
327 3 => b @ a @ v = w @ v
329 This keyword is ignored for standard problems.
330 overwrite_a : bool, optional
331 Whether to overwrite data in ``a`` (may improve performance). Default
332 is False.
333 overwrite_b : bool, optional
334 Whether to overwrite data in ``b`` (may improve performance). Default
335 is False.
336 check_finite : bool, optional
337 Whether to check that the input matrices contain only finite numbers.
338 Disabling may give a performance gain, but may result in problems
339 (crashes, non-termination) if the inputs do contain infinities or NaNs.
340 turbo : bool, optional, deprecated
341 .. deprecated:: 1.5.0
342 `eigh` keyword argument `turbo` is deprecated in favour of
343 ``driver=gvd`` keyword instead and will be removed in SciPy
344 1.12.0.
345 eigvals : tuple (lo, hi), optional, deprecated
346 .. deprecated:: 1.5.0
347 `eigh` keyword argument `eigvals` is deprecated in favour of
348 `subset_by_index` keyword instead and will be removed in SciPy
349 1.12.0.
351 Returns
352 -------
353 w : (N,) ndarray
354 The N (1<=N<=M) selected eigenvalues, in ascending order, each
355 repeated according to its multiplicity.
356 v : (M, N) ndarray
357 (if ``eigvals_only == False``)
359 Raises
360 ------
361 LinAlgError
362 If eigenvalue computation does not converge, an error occurred, or
363 b matrix is not definite positive. Note that if input matrices are
364 not symmetric or Hermitian, no error will be reported but results will
365 be wrong.
367 See Also
368 --------
369 eigvalsh : eigenvalues of symmetric or Hermitian arrays
370 eig : eigenvalues and right eigenvectors for non-symmetric arrays
371 eigh_tridiagonal : eigenvalues and right eiegenvectors for
372 symmetric/Hermitian tridiagonal matrices
374 Notes
375 -----
376 This function does not check the input array for being Hermitian/symmetric
377 in order to allow for representing arrays with only their upper/lower
378 triangular parts. Also, note that even though not taken into account,
379 finiteness check applies to the whole array and unaffected by "lower"
380 keyword.
382 This function uses LAPACK drivers for computations in all possible keyword
383 combinations, prefixed with ``sy`` if arrays are real and ``he`` if
384 complex, e.g., a float array with "evr" driver is solved via
385 "syevr", complex arrays with "gvx" driver problem is solved via "hegvx"
386 etc.
388 As a brief summary, the slowest and the most robust driver is the
389 classical ``<sy/he>ev`` which uses symmetric QR. ``<sy/he>evr`` is seen as
390 the optimal choice for the most general cases. However, there are certain
391 occasions that ``<sy/he>evd`` computes faster at the expense of more
392 memory usage. ``<sy/he>evx``, while still being faster than ``<sy/he>ev``,
393 often performs worse than the rest except when very few eigenvalues are
394 requested for large arrays though there is still no performance guarantee.
397 For the generalized problem, normalization with respect to the given
398 type argument::
400 type 1 and 3 : v.conj().T @ a @ v = w
401 type 2 : inv(v).conj().T @ a @ inv(v) = w
403 type 1 or 2 : v.conj().T @ b @ v = I
404 type 3 : v.conj().T @ inv(b) @ v = I
407 Examples
408 --------
409 >>> import numpy as np
410 >>> from scipy.linalg import eigh
411 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
412 >>> w, v = eigh(A)
413 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
414 True
416 Request only the eigenvalues
418 >>> w = eigh(A, eigvals_only=True)
420 Request eigenvalues that are less than 10.
422 >>> A = np.array([[34, -4, -10, -7, 2],
423 ... [-4, 7, 2, 12, 0],
424 ... [-10, 2, 44, 2, -19],
425 ... [-7, 12, 2, 79, -34],
426 ... [2, 0, -19, -34, 29]])
427 >>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
428 array([6.69199443e-07, 9.11938152e+00])
430 Request the second smallest eigenvalue and its eigenvector
432 >>> w, v = eigh(A, subset_by_index=[1, 1])
433 >>> w
434 array([9.11938152])
435 >>> v.shape # only a single column is returned
436 (5, 1)
438 """
439 if turbo:
440 warnings.warn("Keyword argument 'turbo' is deprecated in favour of '"
441 "driver=gvd' keyword instead and will be removed in "
442 "SciPy 1.12.0.",
443 DeprecationWarning, stacklevel=2)
444 if eigvals:
445 warnings.warn("Keyword argument 'eigvals' is deprecated in favour of "
446 "'subset_by_index' keyword instead and will be removed "
447 "in SciPy 1.12.0.",
448 DeprecationWarning, stacklevel=2)
450 # set lower
451 uplo = 'L' if lower else 'U'
452 # Set job for Fortran routines
453 _job = 'N' if eigvals_only else 'V'
455 drv_str = [None, "ev", "evd", "evr", "evx", "gv", "gvd", "gvx"]
456 if driver not in drv_str:
457 raise ValueError('"{}" is unknown. Possible values are "None", "{}".'
458 ''.format(driver, '", "'.join(drv_str[1:])))
460 a1 = _asarray_validated(a, check_finite=check_finite)
461 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
462 raise ValueError('expected square "a" matrix')
463 overwrite_a = overwrite_a or (_datacopied(a1, a))
464 cplx = True if iscomplexobj(a1) else False
465 n = a1.shape[0]
466 drv_args = {'overwrite_a': overwrite_a}
468 if b is not None:
469 b1 = _asarray_validated(b, check_finite=check_finite)
470 overwrite_b = overwrite_b or _datacopied(b1, b)
471 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
472 raise ValueError('expected square "b" matrix')
474 if b1.shape != a1.shape:
475 raise ValueError("wrong b dimensions {}, should "
476 "be {}".format(b1.shape, a1.shape))
478 if type not in [1, 2, 3]:
479 raise ValueError('"type" keyword only accepts 1, 2, and 3.')
481 cplx = True if iscomplexobj(b1) else (cplx or False)
482 drv_args.update({'overwrite_b': overwrite_b, 'itype': type})
484 # backwards-compatibility handling
485 subset_by_index = subset_by_index if (eigvals is None) else eigvals
487 subset = (subset_by_index is not None) or (subset_by_value is not None)
489 # Both subsets can't be given
490 if subset_by_index and subset_by_value:
491 raise ValueError('Either index or value subset can be requested.')
493 # Take turbo into account if all conditions are met otherwise ignore
494 if turbo and b is not None:
495 driver = 'gvx' if subset else 'gvd'
497 # Check indices if given
498 if subset_by_index:
499 lo, hi = (int(x) for x in subset_by_index)
500 if not (0 <= lo <= hi < n):
501 raise ValueError('Requested eigenvalue indices are not valid. '
502 'Valid range is [0, {}] and start <= end, but '
503 'start={}, end={} is given'.format(n-1, lo, hi))
504 # fortran is 1-indexed
505 drv_args.update({'range': 'I', 'il': lo + 1, 'iu': hi + 1})
507 if subset_by_value:
508 lo, hi = subset_by_value
509 if not (-inf <= lo < hi <= inf):
510 raise ValueError('Requested eigenvalue bounds are not valid. '
511 'Valid range is (-inf, inf) and low < high, but '
512 'low={}, high={} is given'.format(lo, hi))
514 drv_args.update({'range': 'V', 'vl': lo, 'vu': hi})
516 # fix prefix for lapack routines
517 pfx = 'he' if cplx else 'sy'
519 # decide on the driver if not given
520 # first early exit on incompatible choice
521 if driver:
522 if b is None and (driver in ["gv", "gvd", "gvx"]):
523 raise ValueError('{} requires input b array to be supplied '
524 'for generalized eigenvalue problems.'
525 ''.format(driver))
526 if (b is not None) and (driver in ['ev', 'evd', 'evr', 'evx']):
527 raise ValueError('"{}" does not accept input b array '
528 'for standard eigenvalue problems.'
529 ''.format(driver))
530 if subset and (driver in ["ev", "evd", "gv", "gvd"]):
531 raise ValueError('"{}" cannot compute subsets of eigenvalues'
532 ''.format(driver))
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('The leading minor of order {} of B is not '
594 'positive definite. The factorization of B '
595 'could not be completed and no eigenvalues '
596 'or eigenvectors were computed.'.format(info-n))
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('when using select="i", select_range must '
649 'contain integers, got dtype %s (%s)'
650 % (sr.dtype, sr.dtype.char))
651 # translate Python (0 ... N-1) into Fortran (1 ... N) with + 1
652 il, iu = sr + 1
653 if min(il, iu) < 1 or max(il, iu) > max_len:
654 raise ValueError('select_range out of bounds')
655 max_ev = iu - il + 1
656 return select, vl, vu, il, iu, max_ev
659def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
660 select='a', select_range=None, max_ev=0, check_finite=True):
661 """
662 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
664 Find eigenvalues w and optionally right eigenvectors v of a::
666 a v[:,i] = w[i] v[:,i]
667 v.H v = identity
669 The matrix a is stored in a_band either in lower diagonal or upper
670 diagonal ordered form:
672 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
673 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
675 where u is the number of bands above the diagonal.
677 Example of a_band (shape of a is (6,6), u=2)::
679 upper form:
680 * * a02 a13 a24 a35
681 * a01 a12 a23 a34 a45
682 a00 a11 a22 a33 a44 a55
684 lower form:
685 a00 a11 a22 a33 a44 a55
686 a10 a21 a32 a43 a54 *
687 a20 a31 a42 a53 * *
689 Cells marked with * are not used.
691 Parameters
692 ----------
693 a_band : (u+1, M) array_like
694 The bands of the M by M matrix a.
695 lower : bool, optional
696 Is the matrix in the lower form. (Default is upper form)
697 eigvals_only : bool, optional
698 Compute only the eigenvalues and no eigenvectors.
699 (Default: calculate also eigenvectors)
700 overwrite_a_band : bool, optional
701 Discard data in a_band (may enhance performance)
702 select : {'a', 'v', 'i'}, optional
703 Which eigenvalues to calculate
705 ====== ========================================
706 select calculated
707 ====== ========================================
708 'a' All eigenvalues
709 'v' Eigenvalues in the interval (min, max]
710 'i' Eigenvalues with indices min <= i <= max
711 ====== ========================================
712 select_range : (min, max), optional
713 Range of selected eigenvalues
714 max_ev : int, optional
715 For select=='v', maximum number of eigenvalues expected.
716 For other values of select, has no meaning.
718 In doubt, leave this parameter untouched.
720 check_finite : bool, optional
721 Whether to check that the input matrix contains only finite numbers.
722 Disabling may give a performance gain, but may result in problems
723 (crashes, non-termination) if the inputs do contain infinities or NaNs.
725 Returns
726 -------
727 w : (M,) ndarray
728 The eigenvalues, in ascending order, each repeated according to its
729 multiplicity.
730 v : (M, M) float or complex ndarray
731 The normalized eigenvector corresponding to the eigenvalue w[i] is
732 the column v[:,i].
734 Raises
735 ------
736 LinAlgError
737 If eigenvalue computation does not converge.
739 See Also
740 --------
741 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
742 eig : eigenvalues and right eigenvectors of general arrays.
743 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
744 eigh_tridiagonal : eigenvalues and right eigenvectors for
745 symmetric/Hermitian tridiagonal matrices
747 Examples
748 --------
749 >>> import numpy as np
750 >>> from scipy.linalg import eig_banded
751 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
752 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
753 >>> w, v = eig_banded(Ab, lower=True)
754 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
755 True
756 >>> w = eig_banded(Ab, lower=True, eigvals_only=True)
757 >>> w
758 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
760 Request only the eigenvalues between ``[-3, 4]``
762 >>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4])
763 >>> w
764 array([-2.22987175, 3.95222349])
766 """
767 if eigvals_only or overwrite_a_band:
768 a1 = _asarray_validated(a_band, check_finite=check_finite)
769 overwrite_a_band = overwrite_a_band or (_datacopied(a1, a_band))
770 else:
771 a1 = array(a_band)
772 if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
773 raise ValueError("array must not contain infs or NaNs")
774 overwrite_a_band = 1
776 if len(a1.shape) != 2:
777 raise ValueError('expected a 2-D array')
778 select, vl, vu, il, iu, max_ev = _check_select(
779 select, select_range, max_ev, a1.shape[1])
780 del select_range
781 if select == 0:
782 if a1.dtype.char in 'GFD':
783 # FIXME: implement this somewhen, for now go with builtin values
784 # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
785 # or by using calc_lwork.f ???
786 # lwork = calc_lwork.hbevd(bevd.typecode, a1.shape[0], lower)
787 internal_name = 'hbevd'
788 else: # a1.dtype.char in 'fd':
789 # FIXME: implement this somewhen, for now go with builtin values
790 # see above
791 # lwork = calc_lwork.sbevd(bevd.typecode, a1.shape[0], lower)
792 internal_name = 'sbevd'
793 bevd, = get_lapack_funcs((internal_name,), (a1,))
794 w, v, info = bevd(a1, compute_v=not eigvals_only,
795 lower=lower, overwrite_ab=overwrite_a_band)
796 else: # select in [1, 2]
797 if eigvals_only:
798 max_ev = 1
799 # calculate optimal abstol for dsbevx (see manpage)
800 if a1.dtype.char in 'fF': # single precision
801 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),))
802 else:
803 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),))
804 abstol = 2 * lamch('s')
805 if a1.dtype.char in 'GFD':
806 internal_name = 'hbevx'
807 else: # a1.dtype.char in 'gfd'
808 internal_name = 'sbevx'
809 bevx, = get_lapack_funcs((internal_name,), (a1,))
810 w, v, m, ifail, info = bevx(
811 a1, vl, vu, il, iu, compute_v=not eigvals_only, mmax=max_ev,
812 range=select, lower=lower, overwrite_ab=overwrite_a_band,
813 abstol=abstol)
814 # crop off w and v
815 w = w[:m]
816 if not eigvals_only:
817 v = v[:, :m]
818 _check_info(info, internal_name)
820 if eigvals_only:
821 return w
822 return w, v
825def eigvals(a, b=None, overwrite_a=False, check_finite=True,
826 homogeneous_eigvals=False):
827 """
828 Compute eigenvalues from an ordinary or generalized eigenvalue problem.
830 Find eigenvalues of a general matrix::
832 a vr[:,i] = w[i] b vr[:,i]
834 Parameters
835 ----------
836 a : (M, M) array_like
837 A complex or real matrix whose eigenvalues and eigenvectors
838 will be computed.
839 b : (M, M) array_like, optional
840 Right-hand side matrix in a generalized eigenvalue problem.
841 If omitted, identity matrix is assumed.
842 overwrite_a : bool, optional
843 Whether to overwrite data in a (may improve performance)
844 check_finite : bool, optional
845 Whether to check that the input matrices contain only finite numbers.
846 Disabling may give a performance gain, but may result in problems
847 (crashes, non-termination) if the inputs do contain infinities
848 or NaNs.
849 homogeneous_eigvals : bool, optional
850 If True, return the eigenvalues in homogeneous coordinates.
851 In this case ``w`` is a (2, M) array so that::
853 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
855 Default is False.
857 Returns
858 -------
859 w : (M,) or (2, M) double or complex ndarray
860 The eigenvalues, each repeated according to its multiplicity
861 but not in any specific order. The shape is (M,) unless
862 ``homogeneous_eigvals=True``.
864 Raises
865 ------
866 LinAlgError
867 If eigenvalue computation does not converge
869 See Also
870 --------
871 eig : eigenvalues and right eigenvectors of general arrays.
872 eigvalsh : eigenvalues of symmetric or Hermitian arrays
873 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
874 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
875 matrices
877 Examples
878 --------
879 >>> import numpy as np
880 >>> from scipy import linalg
881 >>> a = np.array([[0., -1.], [1., 0.]])
882 >>> linalg.eigvals(a)
883 array([0.+1.j, 0.-1.j])
885 >>> b = np.array([[0., 1.], [1., 1.]])
886 >>> linalg.eigvals(a, b)
887 array([ 1.+0.j, -1.+0.j])
889 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
890 >>> linalg.eigvals(a, homogeneous_eigvals=True)
891 array([[3.+0.j, 8.+0.j, 7.+0.j],
892 [1.+0.j, 1.+0.j, 1.+0.j]])
894 """
895 return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a,
896 check_finite=check_finite,
897 homogeneous_eigvals=homogeneous_eigvals)
900def eigvalsh(a, b=None, lower=True, overwrite_a=False,
901 overwrite_b=False, turbo=False, eigvals=None, type=1,
902 check_finite=True, subset_by_index=None, subset_by_value=None,
903 driver=None):
904 """
905 Solves a standard or generalized eigenvalue problem for a complex
906 Hermitian or real symmetric matrix.
908 Find eigenvalues array ``w`` of array ``a``, where ``b`` is positive
909 definite such that for every eigenvalue λ (i-th entry of w) and its
910 eigenvector vi (i-th column of v) satisfies::
912 a @ vi = λ * b @ vi
913 vi.conj().T @ a @ vi = λ
914 vi.conj().T @ b @ vi = 1
916 In the standard problem, b is assumed to be the identity matrix.
918 Parameters
919 ----------
920 a : (M, M) array_like
921 A complex Hermitian or real symmetric matrix whose eigenvalues will
922 be computed.
923 b : (M, M) array_like, optional
924 A complex Hermitian or real symmetric definite positive matrix in.
925 If omitted, identity matrix is assumed.
926 lower : bool, optional
927 Whether the pertinent array data is taken from the lower or upper
928 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
929 overwrite_a : bool, optional
930 Whether to overwrite data in ``a`` (may improve performance). Default
931 is False.
932 overwrite_b : bool, optional
933 Whether to overwrite data in ``b`` (may improve performance). Default
934 is False.
935 type : int, optional
936 For the generalized problems, this keyword specifies the problem type
937 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
938 inputs)::
940 1 => a @ v = w @ b @ v
941 2 => a @ b @ v = w @ v
942 3 => b @ a @ v = w @ v
944 This keyword is ignored for standard problems.
945 check_finite : bool, optional
946 Whether to check that the input matrices contain only finite numbers.
947 Disabling may give a performance gain, but may result in problems
948 (crashes, non-termination) if the inputs do contain infinities or NaNs.
949 subset_by_index : iterable, optional
950 If provided, this two-element iterable defines the start and the end
951 indices of the desired eigenvalues (ascending order and 0-indexed).
952 To return only the second smallest to fifth smallest eigenvalues,
953 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
954 available with "evr", "evx", and "gvx" drivers. The entries are
955 directly converted to integers via ``int()``.
956 subset_by_value : iterable, optional
957 If provided, this two-element iterable defines the half-open interval
958 ``(a, b]`` that, if any, only the eigenvalues between these values
959 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
960 ``np.inf`` for the unconstrained ends.
961 driver : str, optional
962 Defines which LAPACK driver should be used. Valid options are "ev",
963 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
964 generalized (where b is not None) problems. See the Notes section of
965 `scipy.linalg.eigh`.
966 turbo : bool, optional, deprecated
967 .. deprecated:: 1.5.0
968 'eigvalsh' keyword argument `turbo` is deprecated in favor of
969 ``driver=gvd`` option and will be removed in SciPy 1.12.0.
971 eigvals : tuple (lo, hi), optional
972 .. deprecated:: 1.5.0
973 'eigvalsh' keyword argument `eigvals` is deprecated in favor of
974 `subset_by_index` option and will be removed in SciPy 1.12.0.
976 Returns
977 -------
978 w : (N,) ndarray
979 The ``N`` (``1<=N<=M``) selected eigenvalues, in ascending order, each
980 repeated according to its multiplicity.
982 Raises
983 ------
984 LinAlgError
985 If eigenvalue computation does not converge, an error occurred, or
986 b matrix is not definite positive. Note that if input matrices are
987 not symmetric or Hermitian, no error will be reported but results will
988 be wrong.
990 See Also
991 --------
992 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
993 eigvals : eigenvalues of general arrays
994 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
995 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
996 matrices
998 Notes
999 -----
1000 This function does not check the input array for being Hermitian/symmetric
1001 in order to allow for representing arrays with only their upper/lower
1002 triangular parts.
1004 This function serves as a one-liner shorthand for `scipy.linalg.eigh` with
1005 the option ``eigvals_only=True`` to get the eigenvalues and not the
1006 eigenvectors. Here it is kept as a legacy convenience. It might be
1007 beneficial to use the main function to have full control and to be a bit
1008 more pythonic.
1010 Examples
1011 --------
1012 For more examples see `scipy.linalg.eigh`.
1014 >>> import numpy as np
1015 >>> from scipy.linalg import eigvalsh
1016 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
1017 >>> w = eigvalsh(A)
1018 >>> w
1019 array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
1021 """
1022 return eigh(a, b=b, lower=lower, eigvals_only=True,
1023 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
1024 turbo=turbo, eigvals=eigvals, type=type,
1025 check_finite=check_finite, subset_by_index=subset_by_index,
1026 subset_by_value=subset_by_value, driver=driver)
1029def eigvals_banded(a_band, lower=False, overwrite_a_band=False,
1030 select='a', select_range=None, check_finite=True):
1031 """
1032 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
1034 Find eigenvalues w of a::
1036 a v[:,i] = w[i] v[:,i]
1037 v.H v = identity
1039 The matrix a is stored in a_band either in lower diagonal or upper
1040 diagonal ordered form:
1042 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
1043 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
1045 where u is the number of bands above the diagonal.
1047 Example of a_band (shape of a is (6,6), u=2)::
1049 upper form:
1050 * * a02 a13 a24 a35
1051 * a01 a12 a23 a34 a45
1052 a00 a11 a22 a33 a44 a55
1054 lower form:
1055 a00 a11 a22 a33 a44 a55
1056 a10 a21 a32 a43 a54 *
1057 a20 a31 a42 a53 * *
1059 Cells marked with * are not used.
1061 Parameters
1062 ----------
1063 a_band : (u+1, M) array_like
1064 The bands of the M by M matrix a.
1065 lower : bool, optional
1066 Is the matrix in the lower form. (Default is upper form)
1067 overwrite_a_band : bool, optional
1068 Discard data in a_band (may enhance performance)
1069 select : {'a', 'v', 'i'}, optional
1070 Which eigenvalues to calculate
1072 ====== ========================================
1073 select calculated
1074 ====== ========================================
1075 'a' All eigenvalues
1076 'v' Eigenvalues in the interval (min, max]
1077 'i' Eigenvalues with indices min <= i <= max
1078 ====== ========================================
1079 select_range : (min, max), optional
1080 Range of selected eigenvalues
1081 check_finite : bool, optional
1082 Whether to check that the input matrix contains only finite numbers.
1083 Disabling may give a performance gain, but may result in problems
1084 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1086 Returns
1087 -------
1088 w : (M,) ndarray
1089 The eigenvalues, in ascending order, each repeated according to its
1090 multiplicity.
1092 Raises
1093 ------
1094 LinAlgError
1095 If eigenvalue computation does not converge.
1097 See Also
1098 --------
1099 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1100 band matrices
1101 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1102 matrices
1103 eigvals : eigenvalues of general arrays
1104 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1105 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1107 Examples
1108 --------
1109 >>> import numpy as np
1110 >>> from scipy.linalg import eigvals_banded
1111 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
1112 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
1113 >>> w = eigvals_banded(Ab, lower=True)
1114 >>> w
1115 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
1116 """
1117 return eig_banded(a_band, lower=lower, eigvals_only=1,
1118 overwrite_a_band=overwrite_a_band, select=select,
1119 select_range=select_range, check_finite=check_finite)
1122def eigvalsh_tridiagonal(d, e, select='a', select_range=None,
1123 check_finite=True, tol=0., lapack_driver='auto'):
1124 """
1125 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1127 Find eigenvalues `w` of ``a``::
1129 a v[:,i] = w[i] v[:,i]
1130 v.H v = identity
1132 For a real symmetric matrix ``a`` with diagonal elements `d` and
1133 off-diagonal elements `e`.
1135 Parameters
1136 ----------
1137 d : ndarray, shape (ndim,)
1138 The diagonal elements of the array.
1139 e : ndarray, shape (ndim-1,)
1140 The off-diagonal elements of the array.
1141 select : {'a', 'v', 'i'}, optional
1142 Which eigenvalues to calculate
1144 ====== ========================================
1145 select calculated
1146 ====== ========================================
1147 'a' All eigenvalues
1148 'v' Eigenvalues in the interval (min, max]
1149 'i' Eigenvalues with indices min <= i <= max
1150 ====== ========================================
1151 select_range : (min, max), optional
1152 Range of selected eigenvalues
1153 check_finite : bool, optional
1154 Whether to check that the input matrix contains only finite numbers.
1155 Disabling may give a performance gain, but may result in problems
1156 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1157 tol : float
1158 The absolute tolerance to which each eigenvalue is required
1159 (only used when ``lapack_driver='stebz'``).
1160 An eigenvalue (or cluster) is considered to have converged if it
1161 lies in an interval of this width. If <= 0. (default),
1162 the value ``eps*|a|`` is used where eps is the machine precision,
1163 and ``|a|`` is the 1-norm of the matrix ``a``.
1164 lapack_driver : str
1165 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1166 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1167 and 'stebz' otherwise. 'sterf' and 'stev' can only be used when
1168 ``select='a'``.
1170 Returns
1171 -------
1172 w : (M,) ndarray
1173 The eigenvalues, in ascending order, each repeated according to its
1174 multiplicity.
1176 Raises
1177 ------
1178 LinAlgError
1179 If eigenvalue computation does not converge.
1181 See Also
1182 --------
1183 eigh_tridiagonal : eigenvalues and right eiegenvectors for
1184 symmetric/Hermitian tridiagonal matrices
1186 Examples
1187 --------
1188 >>> import numpy as np
1189 >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
1190 >>> d = 3*np.ones(4)
1191 >>> e = -1*np.ones(3)
1192 >>> w = eigvalsh_tridiagonal(d, e)
1193 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1194 >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
1195 >>> np.allclose(w - w2, np.zeros(4))
1196 True
1197 """
1198 return eigh_tridiagonal(
1199 d, e, eigvals_only=True, select=select, select_range=select_range,
1200 check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
1203def eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None,
1204 check_finite=True, tol=0., lapack_driver='auto'):
1205 """
1206 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1208 Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
1210 a v[:,i] = w[i] v[:,i]
1211 v.H v = identity
1213 For a real symmetric matrix ``a`` with diagonal elements `d` and
1214 off-diagonal elements `e`.
1216 Parameters
1217 ----------
1218 d : ndarray, shape (ndim,)
1219 The diagonal elements of the array.
1220 e : ndarray, shape (ndim-1,)
1221 The off-diagonal elements of the array.
1222 select : {'a', 'v', 'i'}, optional
1223 Which eigenvalues to calculate
1225 ====== ========================================
1226 select calculated
1227 ====== ========================================
1228 'a' All eigenvalues
1229 'v' Eigenvalues in the interval (min, max]
1230 'i' Eigenvalues with indices min <= i <= max
1231 ====== ========================================
1232 select_range : (min, max), optional
1233 Range of selected eigenvalues
1234 check_finite : bool, optional
1235 Whether to check that the input matrix contains only finite numbers.
1236 Disabling may give a performance gain, but may result in problems
1237 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1238 tol : float
1239 The absolute tolerance to which each eigenvalue is required
1240 (only used when 'stebz' is the `lapack_driver`).
1241 An eigenvalue (or cluster) is considered to have converged if it
1242 lies in an interval of this width. If <= 0. (default),
1243 the value ``eps*|a|`` is used where eps is the machine precision,
1244 and ``|a|`` is the 1-norm of the matrix ``a``.
1245 lapack_driver : str
1246 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1247 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1248 and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and
1249 ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is
1250 used to find the corresponding eigenvectors. 'sterf' can only be
1251 used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only
1252 be used when ``select='a'``.
1254 Returns
1255 -------
1256 w : (M,) ndarray
1257 The eigenvalues, in ascending order, each repeated according to its
1258 multiplicity.
1259 v : (M, M) ndarray
1260 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
1261 the column ``v[:,i]``.
1263 Raises
1264 ------
1265 LinAlgError
1266 If eigenvalue computation does not converge.
1268 See Also
1269 --------
1270 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1271 matrices
1272 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1273 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1274 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1275 band matrices
1277 Notes
1278 -----
1279 This function makes use of LAPACK ``S/DSTEMR`` routines.
1281 Examples
1282 --------
1283 >>> import numpy as np
1284 >>> from scipy.linalg import eigh_tridiagonal
1285 >>> d = 3*np.ones(4)
1286 >>> e = -1*np.ones(3)
1287 >>> w, v = eigh_tridiagonal(d, e)
1288 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1289 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
1290 True
1291 """
1292 d = _asarray_validated(d, check_finite=check_finite)
1293 e = _asarray_validated(e, check_finite=check_finite)
1294 for check in (d, e):
1295 if check.ndim != 1:
1296 raise ValueError('expected a 1-D array')
1297 if check.dtype.char in 'GFD': # complex
1298 raise TypeError('Only real arrays currently supported')
1299 if d.size != e.size + 1:
1300 raise ValueError('d (%s) must have one more element than e (%s)'
1301 % (d.size, e.size))
1302 select, vl, vu, il, iu, _ = _check_select(
1303 select, select_range, 0, d.size)
1304 if not isinstance(lapack_driver, str):
1305 raise TypeError('lapack_driver must be str')
1306 drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev')
1307 if lapack_driver not in drivers:
1308 raise ValueError('lapack_driver must be one of %s, got %s'
1309 % (drivers, lapack_driver))
1310 if lapack_driver == 'auto':
1311 lapack_driver = 'stemr' if select == 0 else 'stebz'
1312 func, = get_lapack_funcs((lapack_driver,), (d, e))
1313 compute_v = not eigvals_only
1314 if lapack_driver == 'sterf':
1315 if select != 0:
1316 raise ValueError('sterf can only be used when select == "a"')
1317 if not eigvals_only:
1318 raise ValueError('sterf can only be used when eigvals_only is '
1319 'True')
1320 w, info = func(d, e)
1321 m = len(w)
1322 elif lapack_driver == 'stev':
1323 if select != 0:
1324 raise ValueError('stev can only be used when select == "a"')
1325 w, v, info = func(d, e, compute_v=compute_v)
1326 m = len(w)
1327 elif lapack_driver == 'stebz':
1328 tol = float(tol)
1329 internal_name = 'stebz'
1330 stebz, = get_lapack_funcs((internal_name,), (d, e))
1331 # If getting eigenvectors, needs to be block-ordered (B) instead of
1332 # matrix-ordered (E), and we will reorder later
1333 order = 'E' if eigvals_only else 'B'
1334 m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol,
1335 order)
1336 else: # 'stemr'
1337 # ?STEMR annoyingly requires size N instead of N-1
1338 e_ = empty(e.size+1, e.dtype)
1339 e_[:-1] = e
1340 stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e))
1341 lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu,
1342 compute_v=compute_v)
1343 _check_info(info, 'stemr_lwork')
1344 m, w, v, info = func(d, e_, select, vl, vu, il, iu,
1345 compute_v=compute_v, lwork=lwork, liwork=liwork)
1346 _check_info(info, lapack_driver + ' (eigh_tridiagonal)')
1347 w = w[:m]
1348 if eigvals_only:
1349 return w
1350 else:
1351 # Do we still need to compute the eigenvalues?
1352 if lapack_driver == 'stebz':
1353 func, = get_lapack_funcs(('stein',), (d, e))
1354 v, info = func(d, e, w, iblock, isplit)
1355 _check_info(info, 'stein (eigh_tridiagonal)',
1356 positive='%d eigenvectors failed to converge')
1357 # Convert block-order to matrix-order
1358 order = argsort(w)
1359 w, v = w[order], v[:, order]
1360 else:
1361 v = v[:, :m]
1362 return w, v
1365def _check_info(info, driver, positive='did not converge (LAPACK info=%d)'):
1366 """Check info return value."""
1367 if info < 0:
1368 raise ValueError('illegal value in argument %d of internal %s'
1369 % (-info, driver))
1370 if info > 0 and positive:
1371 raise LinAlgError(("%s " + positive) % (driver, info,))
1374def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
1375 """
1376 Compute Hessenberg form of a matrix.
1378 The Hessenberg decomposition is::
1380 A = Q H Q^H
1382 where `Q` is unitary/orthogonal and `H` has only zero elements below
1383 the first sub-diagonal.
1385 Parameters
1386 ----------
1387 a : (M, M) array_like
1388 Matrix to bring into Hessenberg form.
1389 calc_q : bool, optional
1390 Whether to compute the transformation matrix. Default is False.
1391 overwrite_a : bool, optional
1392 Whether to overwrite `a`; may improve performance.
1393 Default is False.
1394 check_finite : bool, optional
1395 Whether to check that the input matrix contains only finite numbers.
1396 Disabling may give a performance gain, but may result in problems
1397 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1399 Returns
1400 -------
1401 H : (M, M) ndarray
1402 Hessenberg form of `a`.
1403 Q : (M, M) ndarray
1404 Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
1405 Only returned if ``calc_q=True``.
1407 Examples
1408 --------
1409 >>> import numpy as np
1410 >>> from scipy.linalg import hessenberg
1411 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
1412 >>> H, Q = hessenberg(A, calc_q=True)
1413 >>> H
1414 array([[ 2. , -11.65843866, 1.42005301, 0.25349066],
1415 [ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
1416 [ 0. , -1.83299243, 0.38969961, -0.51527034],
1417 [ 0. , 0. , -3.83189513, 1.07494686]])
1418 >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
1419 True
1420 """
1421 a1 = _asarray_validated(a, check_finite=check_finite)
1422 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
1423 raise ValueError('expected square matrix')
1424 overwrite_a = overwrite_a or (_datacopied(a1, a))
1426 # if 2x2 or smaller: already in Hessenberg
1427 if a1.shape[0] <= 2:
1428 if calc_q:
1429 return a1, eye(a1.shape[0])
1430 return a1
1432 gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal',
1433 'gehrd_lwork'), (a1,))
1434 ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a)
1435 _check_info(info, 'gebal (hessenberg)', positive=False)
1436 n = len(a1)
1438 lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
1440 hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1441 _check_info(info, 'gehrd (hessenberg)', positive=False)
1442 h = numpy.triu(hq, -1)
1443 if not calc_q:
1444 return h
1446 # use orghr/unghr to compute q
1447 orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,))
1448 lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
1450 q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1451 _check_info(info, 'orghr (hessenberg)', positive=False)
1452 return h, q
1455def cdf2rdf(w, v):
1456 """
1457 Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real
1458 eigenvalues in a block diagonal form ``wr`` and the associated real
1459 eigenvectors ``vr``, such that::
1461 vr @ wr = X @ vr
1463 continues to hold, where ``X`` is the original array for which ``w`` and
1464 ``v`` are the eigenvalues and eigenvectors.
1466 .. versionadded:: 1.1.0
1468 Parameters
1469 ----------
1470 w : (..., M) array_like
1471 Complex or real eigenvalues, an array or stack of arrays
1473 Conjugate pairs must not be interleaved, else the wrong result
1474 will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result,
1475 but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
1477 v : (..., M, M) array_like
1478 Complex or real eigenvectors, a square array or stack of square arrays.
1480 Returns
1481 -------
1482 wr : (..., M, M) ndarray
1483 Real diagonal block form of eigenvalues
1484 vr : (..., M, M) ndarray
1485 Real eigenvectors associated with ``wr``
1487 See Also
1488 --------
1489 eig : Eigenvalues and right eigenvectors for non-symmetric arrays
1490 rsf2csf : Convert real Schur form to complex Schur form
1492 Notes
1493 -----
1494 ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``.
1495 For example, obtained by ``w, v = scipy.linalg.eig(X)`` or
1496 ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent
1497 stacked arrays.
1499 .. versionadded:: 1.1.0
1501 Examples
1502 --------
1503 >>> import numpy as np
1504 >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
1505 >>> X
1506 array([[ 1, 2, 3],
1507 [ 0, 4, 5],
1508 [ 0, -5, 4]])
1510 >>> from scipy import linalg
1511 >>> w, v = linalg.eig(X)
1512 >>> w
1513 array([ 1.+0.j, 4.+5.j, 4.-5.j])
1514 >>> v
1515 array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
1516 [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
1517 [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
1519 >>> wr, vr = linalg.cdf2rdf(w, v)
1520 >>> wr
1521 array([[ 1., 0., 0.],
1522 [ 0., 4., 5.],
1523 [ 0., -5., 4.]])
1524 >>> vr
1525 array([[ 1. , 0.40016, -0.01906],
1526 [ 0. , 0.64788, 0. ],
1527 [ 0. , 0. , 0.64788]])
1529 >>> vr @ wr
1530 array([[ 1. , 1.69593, 1.9246 ],
1531 [ 0. , 2.59153, 3.23942],
1532 [ 0. , -3.23942, 2.59153]])
1533 >>> X @ vr
1534 array([[ 1. , 1.69593, 1.9246 ],
1535 [ 0. , 2.59153, 3.23942],
1536 [ 0. , -3.23942, 2.59153]])
1537 """
1538 w, v = _asarray_validated(w), _asarray_validated(v)
1540 # check dimensions
1541 if w.ndim < 1:
1542 raise ValueError('expected w to be at least 1D')
1543 if v.ndim < 2:
1544 raise ValueError('expected v to be at least 2D')
1545 if v.ndim != w.ndim + 1:
1546 raise ValueError('expected eigenvectors array to have exactly one '
1547 'dimension more than eigenvalues array')
1549 # check shapes
1550 n = w.shape[-1]
1551 M = w.shape[:-1]
1552 if v.shape[-2] != v.shape[-1]:
1553 raise ValueError('expected v to be a square matrix or stacked square '
1554 'matrices: v.shape[-2] = v.shape[-1]')
1555 if v.shape[-1] != n:
1556 raise ValueError('expected the same number of eigenvalues as '
1557 'eigenvectors')
1559 # get indices for each first pair of complex eigenvalues
1560 complex_mask = iscomplex(w)
1561 n_complex = complex_mask.sum(axis=-1)
1563 # check if all complex eigenvalues have conjugate pairs
1564 if not (n_complex % 2 == 0).all():
1565 raise ValueError('expected complex-conjugate pairs of eigenvalues')
1567 # find complex indices
1568 idx = nonzero(complex_mask)
1569 idx_stack = idx[:-1]
1570 idx_elem = idx[-1]
1572 # filter them to conjugate indices, assuming pairs are not interleaved
1573 j = idx_elem[0::2]
1574 k = idx_elem[1::2]
1575 stack_ind = ()
1576 for i in idx_stack:
1577 # should never happen, assuming nonzero orders by the last axis
1578 assert (i[0::2] == i[1::2]).all(),\
1579 "Conjugate pair spanned different arrays!"
1580 stack_ind += (i[0::2],)
1582 # all eigenvalues to diagonal form
1583 wr = zeros(M + (n, n), dtype=w.real.dtype)
1584 di = range(n)
1585 wr[..., di, di] = w.real
1587 # complex eigenvalues to real block diagonal form
1588 wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag
1589 wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
1591 # compute real eigenvectors associated with real block diagonal eigenvalues
1592 u = zeros(M + (n, n), dtype=numpy.cdouble)
1593 u[..., di, di] = 1.0
1594 u[stack_ind + (j, j)] = 0.5j
1595 u[stack_ind + (j, k)] = 0.5
1596 u[stack_ind + (k, j)] = -0.5j
1597 u[stack_ind + (k, k)] = 0.5
1599 # multipy matrices v and u (equivalent to v @ u)
1600 vr = einsum('...ij,...jk->...ik', v, u).real
1602 return wr, vr