Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp.py: 6%
390 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1# -*- coding: utf-8 -*-
2#
3# Author: Pearu Peterson, March 2002
4#
5# additions by Travis Oliphant, March 2002
6# additions by Eric Jones, June 2002
7# additions by Johannes Loehnert, June 2006
8# additions by Bart Vandereycken, June 2006
9# additions by Andrew D Straw, May 2007
10# additions by Tiziano Zito, November 2008
11#
12# April 2010: Functions for LU, QR, SVD, Schur, and Cholesky decompositions
13# were moved to their own files. Still in this file are functions for
14# eigenstuff and for the Hessenberg form.
16__all__ = ['eig', 'eigvals', 'eigh', 'eigvalsh',
17 'eig_banded', 'eigvals_banded',
18 'eigh_tridiagonal', 'eigvalsh_tridiagonal', 'hessenberg', 'cdf2rdf']
20import warnings
22import numpy
23from numpy import (array, isfinite, inexact, nonzero, iscomplexobj, cast,
24 flatnonzero, conj, asarray, argsort, empty,
25 iscomplex, zeros, einsum, eye, inf)
26# Local imports
27from scipy._lib._util import _asarray_validated
28from ._misc import LinAlgError, _datacopied, norm
29from .lapack import get_lapack_funcs, _compute_lwork
32_I = cast['F'](1j)
35def _make_complex_eigvecs(w, vin, dtype):
36 """
37 Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output
38 """
39 # - see LAPACK man page DGGEV at ALPHAI
40 v = numpy.array(vin, dtype=dtype)
41 m = (w.imag > 0)
42 m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709
43 for i in flatnonzero(m):
44 v.imag[:, i] = vin[:, i+1]
45 conj(v[:, i], v[:, i+1])
46 return v
49def _make_eigvals(alpha, beta, homogeneous_eigvals):
50 if homogeneous_eigvals:
51 if beta is None:
52 return numpy.vstack((alpha, numpy.ones_like(alpha)))
53 else:
54 return numpy.vstack((alpha, beta))
55 else:
56 if beta is None:
57 return alpha
58 else:
59 w = numpy.empty_like(alpha)
60 alpha_zero = (alpha == 0)
61 beta_zero = (beta == 0)
62 beta_nonzero = ~beta_zero
63 w[beta_nonzero] = alpha[beta_nonzero]/beta[beta_nonzero]
64 # Use numpy.inf for complex values too since
65 # 1/numpy.inf = 0, i.e., it correctly behaves as projective
66 # infinity.
67 w[~alpha_zero & beta_zero] = numpy.inf
68 if numpy.all(alpha.imag == 0):
69 w[alpha_zero & beta_zero] = numpy.nan
70 else:
71 w[alpha_zero & beta_zero] = complex(numpy.nan, numpy.nan)
72 return w
75def _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
76 homogeneous_eigvals):
77 ggev, = get_lapack_funcs(('ggev',), (a1, b1))
78 cvl, cvr = left, right
79 res = ggev(a1, b1, lwork=-1)
80 lwork = res[-2][0].real.astype(numpy.int_)
81 if ggev.typecode in 'cz':
82 alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
83 overwrite_a, overwrite_b)
84 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
85 else:
86 alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr,
87 lwork, overwrite_a,
88 overwrite_b)
89 alpha = alphar + _I * alphai
90 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
91 _check_info(info, 'generalized eig algorithm (ggev)')
93 only_real = numpy.all(w.imag == 0.0)
94 if not (ggev.typecode in 'cz' or only_real):
95 t = w.dtype.char
96 if left:
97 vl = _make_complex_eigvecs(w, vl, t)
98 if right:
99 vr = _make_complex_eigvecs(w, vr, t)
101 # the eigenvectors returned by the lapack function are NOT normalized
102 for i in range(vr.shape[0]):
103 if right:
104 vr[:, i] /= norm(vr[:, i])
105 if left:
106 vl[:, i] /= norm(vl[:, i])
108 if not (left or right):
109 return w
110 if left:
111 if right:
112 return w, vl, vr
113 return w, vl
114 return w, vr
117def eig(a, b=None, left=False, right=True, overwrite_a=False,
118 overwrite_b=False, check_finite=True, homogeneous_eigvals=False):
119 """
120 Solve an ordinary or generalized eigenvalue problem of a square matrix.
122 Find eigenvalues w and right or left eigenvectors of a general matrix::
124 a vr[:,i] = w[i] b vr[:,i]
125 a.H vl[:,i] = w[i].conj() b.H vl[:,i]
127 where ``.H`` is the Hermitian conjugation.
129 Parameters
130 ----------
131 a : (M, M) array_like
132 A complex or real matrix whose eigenvalues and eigenvectors
133 will be computed.
134 b : (M, M) array_like, optional
135 Right-hand side matrix in a generalized eigenvalue problem.
136 Default is None, identity matrix is assumed.
137 left : bool, optional
138 Whether to calculate and return left eigenvectors. Default is False.
139 right : bool, optional
140 Whether to calculate and return right eigenvectors. Default is True.
141 overwrite_a : bool, optional
142 Whether to overwrite `a`; may improve performance. Default is False.
143 overwrite_b : bool, optional
144 Whether to overwrite `b`; may improve performance. Default is False.
145 check_finite : bool, optional
146 Whether to check that the input matrices contain only finite numbers.
147 Disabling may give a performance gain, but may result in problems
148 (crashes, non-termination) if the inputs do contain infinities or NaNs.
149 homogeneous_eigvals : bool, optional
150 If True, return the eigenvalues in homogeneous coordinates.
151 In this case ``w`` is a (2, M) array so that::
153 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
155 Default is False.
157 Returns
158 -------
159 w : (M,) or (2, M) double or complex ndarray
160 The eigenvalues, each repeated according to its
161 multiplicity. The shape is (M,) unless
162 ``homogeneous_eigvals=True``.
163 vl : (M, M) double or complex ndarray
164 The normalized left eigenvector corresponding to the eigenvalue
165 ``w[i]`` is the column vl[:,i]. Only returned if ``left=True``.
166 vr : (M, M) double or complex ndarray
167 The normalized right eigenvector corresponding to the eigenvalue
168 ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
170 Raises
171 ------
172 LinAlgError
173 If eigenvalue computation does not converge.
175 See Also
176 --------
177 eigvals : eigenvalues of general arrays
178 eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
179 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
180 band matrices
181 eigh_tridiagonal : eigenvalues and right eiegenvectors for
182 symmetric/Hermitian tridiagonal matrices
184 Examples
185 --------
186 >>> import numpy as np
187 >>> from scipy import linalg
188 >>> a = np.array([[0., -1.], [1., 0.]])
189 >>> linalg.eigvals(a)
190 array([0.+1.j, 0.-1.j])
192 >>> b = np.array([[0., 1.], [1., 1.]])
193 >>> linalg.eigvals(a, b)
194 array([ 1.+0.j, -1.+0.j])
196 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
197 >>> linalg.eigvals(a, homogeneous_eigvals=True)
198 array([[3.+0.j, 8.+0.j, 7.+0.j],
199 [1.+0.j, 1.+0.j, 1.+0.j]])
201 >>> a = np.array([[0., -1.], [1., 0.]])
202 >>> linalg.eigvals(a) == linalg.eig(a)[0]
203 array([ True, True])
204 >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector
205 array([[-0.70710678+0.j , -0.70710678-0.j ],
206 [-0. +0.70710678j, -0. -0.70710678j]])
207 >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector
208 array([[0.70710678+0.j , 0.70710678-0.j ],
209 [0. -0.70710678j, 0. +0.70710678j]])
213 """
214 a1 = _asarray_validated(a, check_finite=check_finite)
215 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
216 raise ValueError('expected square matrix')
217 overwrite_a = overwrite_a or (_datacopied(a1, a))
218 if b is not None:
219 b1 = _asarray_validated(b, check_finite=check_finite)
220 overwrite_b = overwrite_b or _datacopied(b1, b)
221 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
222 raise ValueError('expected square matrix')
223 if b1.shape != a1.shape:
224 raise ValueError('a and b must have the same shape')
225 return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
226 homogeneous_eigvals)
228 geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
229 compute_vl, compute_vr = left, right
231 lwork = _compute_lwork(geev_lwork, a1.shape[0],
232 compute_vl=compute_vl,
233 compute_vr=compute_vr)
235 if geev.typecode in 'cz':
236 w, vl, vr, info = geev(a1, lwork=lwork,
237 compute_vl=compute_vl,
238 compute_vr=compute_vr,
239 overwrite_a=overwrite_a)
240 w = _make_eigvals(w, None, homogeneous_eigvals)
241 else:
242 wr, wi, vl, vr, info = geev(a1, lwork=lwork,
243 compute_vl=compute_vl,
244 compute_vr=compute_vr,
245 overwrite_a=overwrite_a)
246 t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
247 w = wr + _I * wi
248 w = _make_eigvals(w, None, homogeneous_eigvals)
250 _check_info(info, 'eig algorithm (geev)',
251 positive='did not converge (only eigenvalues '
252 'with order >= %d have converged)')
254 only_real = numpy.all(w.imag == 0.0)
255 if not (geev.typecode in 'cz' or only_real):
256 t = w.dtype.char
257 if left:
258 vl = _make_complex_eigvecs(w, vl, t)
259 if right:
260 vr = _make_complex_eigvecs(w, vr, t)
261 if not (left or right):
262 return w
263 if left:
264 if right:
265 return w, vl, vr
266 return w, vl
267 return w, vr
270def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
271 overwrite_b=False, turbo=False, eigvals=None, type=1,
272 check_finite=True, subset_by_index=None, subset_by_value=None,
273 driver=None):
274 """
275 Solve a standard or generalized eigenvalue problem for a complex
276 Hermitian or real symmetric matrix.
278 Find eigenvalues array ``w`` and optionally eigenvectors array ``v`` of
279 array ``a``, where ``b`` is positive definite such that for every
280 eigenvalue λ (i-th entry of w) and its eigenvector ``vi`` (i-th column of
281 ``v``) satisfies::
283 a @ vi = λ * b @ vi
284 vi.conj().T @ a @ vi = λ
285 vi.conj().T @ b @ vi = 1
287 In the standard problem, ``b`` is assumed to be the identity matrix.
289 Parameters
290 ----------
291 a : (M, M) array_like
292 A complex Hermitian or real symmetric matrix whose eigenvalues and
293 eigenvectors will be computed.
294 b : (M, M) array_like, optional
295 A complex Hermitian or real symmetric definite positive matrix in.
296 If omitted, identity matrix is assumed.
297 lower : bool, optional
298 Whether the pertinent array data is taken from the lower or upper
299 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
300 eigvals_only : bool, optional
301 Whether to calculate only eigenvalues and no eigenvectors.
302 (Default: both are calculated)
303 subset_by_index : iterable, optional
304 If provided, this two-element iterable defines the start and the end
305 indices of the desired eigenvalues (ascending order and 0-indexed).
306 To return only the second smallest to fifth smallest eigenvalues,
307 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
308 available with "evr", "evx", and "gvx" drivers. The entries are
309 directly converted to integers via ``int()``.
310 subset_by_value : iterable, optional
311 If provided, this two-element iterable defines the half-open interval
312 ``(a, b]`` that, if any, only the eigenvalues between these values
313 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
314 ``np.inf`` for the unconstrained ends.
315 driver : str, optional
316 Defines which LAPACK driver should be used. Valid options are "ev",
317 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
318 generalized (where b is not None) problems. See the Notes section.
319 The default for standard problems is "evr". For generalized problems,
320 "gvd" is used for full set, and "gvx" for subset requested cases.
321 type : int, optional
322 For the generalized problems, this keyword specifies the problem type
323 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
324 inputs)::
326 1 => a @ v = w @ b @ v
327 2 => a @ b @ v = w @ v
328 3 => b @ a @ v = w @ v
330 This keyword is ignored for standard problems.
331 overwrite_a : bool, optional
332 Whether to overwrite data in ``a`` (may improve performance). Default
333 is False.
334 overwrite_b : bool, optional
335 Whether to overwrite data in ``b`` (may improve performance). Default
336 is False.
337 check_finite : bool, optional
338 Whether to check that the input matrices contain only finite numbers.
339 Disabling may give a performance gain, but may result in problems
340 (crashes, non-termination) if the inputs do contain infinities or NaNs.
341 turbo : bool, optional, deprecated
342 .. deprecated:: 1.5.0
343 `eigh` keyword argument `turbo` is deprecated in favour of
344 ``driver=gvd`` keyword instead and will be removed in SciPy
345 1.12.0.
346 eigvals : tuple (lo, hi), optional, deprecated
347 .. deprecated:: 1.5.0
348 `eigh` keyword argument `eigvals` is deprecated in favour of
349 `subset_by_index` keyword instead and will be removed in SciPy
350 1.12.0.
352 Returns
353 -------
354 w : (N,) ndarray
355 The N (1<=N<=M) selected eigenvalues, in ascending order, each
356 repeated according to its multiplicity.
357 v : (M, N) ndarray
358 (if ``eigvals_only == False``)
360 Raises
361 ------
362 LinAlgError
363 If eigenvalue computation does not converge, an error occurred, or
364 b matrix is not definite positive. Note that if input matrices are
365 not symmetric or Hermitian, no error will be reported but results will
366 be wrong.
368 See Also
369 --------
370 eigvalsh : eigenvalues of symmetric or Hermitian arrays
371 eig : eigenvalues and right eigenvectors for non-symmetric arrays
372 eigh_tridiagonal : eigenvalues and right eiegenvectors for
373 symmetric/Hermitian tridiagonal matrices
375 Notes
376 -----
377 This function does not check the input array for being Hermitian/symmetric
378 in order to allow for representing arrays with only their upper/lower
379 triangular parts. Also, note that even though not taken into account,
380 finiteness check applies to the whole array and unaffected by "lower"
381 keyword.
383 This function uses LAPACK drivers for computations in all possible keyword
384 combinations, prefixed with ``sy`` if arrays are real and ``he`` if
385 complex, e.g., a float array with "evr" driver is solved via
386 "syevr", complex arrays with "gvx" driver problem is solved via "hegvx"
387 etc.
389 As a brief summary, the slowest and the most robust driver is the
390 classical ``<sy/he>ev`` which uses symmetric QR. ``<sy/he>evr`` is seen as
391 the optimal choice for the most general cases. However, there are certain
392 occasions that ``<sy/he>evd`` computes faster at the expense of more
393 memory usage. ``<sy/he>evx``, while still being faster than ``<sy/he>ev``,
394 often performs worse than the rest except when very few eigenvalues are
395 requested for large arrays though there is still no performance guarantee.
398 For the generalized problem, normalization with respect to the given
399 type argument::
401 type 1 and 3 : v.conj().T @ a @ v = w
402 type 2 : inv(v).conj().T @ a @ inv(v) = w
404 type 1 or 2 : v.conj().T @ b @ v = I
405 type 3 : v.conj().T @ inv(b) @ v = I
408 Examples
409 --------
410 >>> import numpy as np
411 >>> from scipy.linalg import eigh
412 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
413 >>> w, v = eigh(A)
414 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
415 True
417 Request only the eigenvalues
419 >>> w = eigh(A, eigvals_only=True)
421 Request eigenvalues that are less than 10.
423 >>> A = np.array([[34, -4, -10, -7, 2],
424 ... [-4, 7, 2, 12, 0],
425 ... [-10, 2, 44, 2, -19],
426 ... [-7, 12, 2, 79, -34],
427 ... [2, 0, -19, -34, 29]])
428 >>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
429 array([6.69199443e-07, 9.11938152e+00])
431 Request the second smallest eigenvalue and its eigenvector
433 >>> w, v = eigh(A, subset_by_index=[1, 1])
434 >>> w
435 array([9.11938152])
436 >>> v.shape # only a single column is returned
437 (5, 1)
439 """
440 if turbo:
441 warnings.warn("Keyword argument 'turbo' is deprecated in favour of '"
442 "driver=gvd' keyword instead and will be removed in "
443 "SciPy 1.12.0.",
444 DeprecationWarning, stacklevel=2)
445 if eigvals:
446 warnings.warn("Keyword argument 'eigvals' is deprecated in favour of "
447 "'subset_by_index' keyword instead and will be removed "
448 "in SciPy 1.12.0.",
449 DeprecationWarning, stacklevel=2)
451 # set lower
452 uplo = 'L' if lower else 'U'
453 # Set job for Fortran routines
454 _job = 'N' if eigvals_only else 'V'
456 drv_str = [None, "ev", "evd", "evr", "evx", "gv", "gvd", "gvx"]
457 if driver not in drv_str:
458 raise ValueError('"{}" is unknown. Possible values are "None", "{}".'
459 ''.format(driver, '", "'.join(drv_str[1:])))
461 a1 = _asarray_validated(a, check_finite=check_finite)
462 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
463 raise ValueError('expected square "a" matrix')
464 overwrite_a = overwrite_a or (_datacopied(a1, a))
465 cplx = True if iscomplexobj(a1) else False
466 n = a1.shape[0]
467 drv_args = {'overwrite_a': overwrite_a}
469 if b is not None:
470 b1 = _asarray_validated(b, check_finite=check_finite)
471 overwrite_b = overwrite_b or _datacopied(b1, b)
472 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
473 raise ValueError('expected square "b" matrix')
475 if b1.shape != a1.shape:
476 raise ValueError("wrong b dimensions {}, should "
477 "be {}".format(b1.shape, a1.shape))
479 if type not in [1, 2, 3]:
480 raise ValueError('"type" keyword only accepts 1, 2, and 3.')
482 cplx = True if iscomplexobj(b1) else (cplx or False)
483 drv_args.update({'overwrite_b': overwrite_b, 'itype': type})
485 # backwards-compatibility handling
486 subset_by_index = subset_by_index if (eigvals is None) else eigvals
488 subset = (subset_by_index is not None) or (subset_by_value is not None)
490 # Both subsets can't be given
491 if subset_by_index and subset_by_value:
492 raise ValueError('Either index or value subset can be requested.')
494 # Take turbo into account if all conditions are met otherwise ignore
495 if turbo and b is not None:
496 driver = 'gvx' if subset else 'gvd'
498 # Check indices if given
499 if subset_by_index:
500 lo, hi = [int(x) for x in subset_by_index]
501 if not (0 <= lo <= hi < n):
502 raise ValueError('Requested eigenvalue indices are not valid. '
503 'Valid range is [0, {}] and start <= end, but '
504 'start={}, end={} is given'.format(n-1, lo, hi))
505 # fortran is 1-indexed
506 drv_args.update({'range': 'I', 'il': lo + 1, 'iu': hi + 1})
508 if subset_by_value:
509 lo, hi = subset_by_value
510 if not (-inf <= lo < hi <= inf):
511 raise ValueError('Requested eigenvalue bounds are not valid. '
512 'Valid range is (-inf, inf) and low < high, but '
513 'low={}, high={} is given'.format(lo, hi))
515 drv_args.update({'range': 'V', 'vl': lo, 'vu': hi})
517 # fix prefix for lapack routines
518 pfx = 'he' if cplx else 'sy'
520 # decide on the driver if not given
521 # first early exit on incompatible choice
522 if driver:
523 if b is None and (driver in ["gv", "gvd", "gvx"]):
524 raise ValueError('{} requires input b array to be supplied '
525 'for generalized eigenvalue problems.'
526 ''.format(driver))
527 if (b is not None) and (driver in ['ev', 'evd', 'evr', 'evx']):
528 raise ValueError('"{}" does not accept input b array '
529 'for standard eigenvalue problems.'
530 ''.format(driver))
531 if subset and (driver in ["ev", "evd", "gv", "gvd"]):
532 raise ValueError('"{}" cannot compute subsets of eigenvalues'
533 ''.format(driver))
535 # Default driver is evr and gvd
536 else:
537 driver = "evr" if b is None else ("gvx" if subset else "gvd")
539 lwork_spec = {
540 'syevd': ['lwork', 'liwork'],
541 'syevr': ['lwork', 'liwork'],
542 'heevd': ['lwork', 'liwork', 'lrwork'],
543 'heevr': ['lwork', 'lrwork', 'liwork'],
544 }
546 if b is None: # Standard problem
547 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
548 [a1])
549 clw_args = {'n': n, 'lower': lower}
550 if driver == 'evd':
551 clw_args.update({'compute_v': 0 if _job == "N" else 1})
553 lw = _compute_lwork(drvlw, **clw_args)
554 # Multiple lwork vars
555 if isinstance(lw, tuple):
556 lwork_args = dict(zip(lwork_spec[pfx+driver], lw))
557 else:
558 lwork_args = {'lwork': lw}
560 drv_args.update({'lower': lower, 'compute_v': 0 if _job == "N" else 1})
561 w, v, *other_args, info = drv(a=a1, **drv_args, **lwork_args)
563 else: # Generalized problem
564 # 'gvd' doesn't have lwork query
565 if driver == "gvd":
566 drv = get_lapack_funcs(pfx + "gvd", [a1, b1])
567 lwork_args = {}
568 else:
569 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
570 [a1, b1])
571 # generalized drivers use uplo instead of lower
572 lw = _compute_lwork(drvlw, n, uplo=uplo)
573 lwork_args = {'lwork': lw}
575 drv_args.update({'uplo': uplo, 'jobz': _job})
577 w, v, *other_args, info = drv(a=a1, b=b1, **drv_args, **lwork_args)
579 # m is always the first extra argument
580 w = w[:other_args[0]] if subset else w
581 v = v[:, :other_args[0]] if (subset and not eigvals_only) else v
583 # Check if we had a successful exit
584 if info == 0:
585 if eigvals_only:
586 return w
587 else:
588 return w, v
589 else:
590 if info < -1:
591 raise LinAlgError('Illegal value in argument {} of internal {}'
592 ''.format(-info, drv.typecode + pfx + driver))
593 elif info > n:
594 raise LinAlgError('The leading minor of order {} of B is not '
595 'positive definite. The factorization of B '
596 'could not be completed and no eigenvalues '
597 'or eigenvectors were computed.'.format(info-n))
598 else:
599 drv_err = {'ev': 'The algorithm failed to converge; {} '
600 'off-diagonal elements of an intermediate '
601 'tridiagonal form did not converge to zero.',
602 'evx': '{} eigenvectors failed to converge.',
603 'evd': 'The algorithm failed to compute an eigenvalue '
604 'while working on the submatrix lying in rows '
605 'and columns {0}/{1} through mod({0},{1}).',
606 'evr': 'Internal Error.'
607 }
608 if driver in ['ev', 'gv']:
609 msg = drv_err['ev'].format(info)
610 elif driver in ['evx', 'gvx']:
611 msg = drv_err['evx'].format(info)
612 elif driver in ['evd', 'gvd']:
613 if eigvals_only:
614 msg = drv_err['ev'].format(info)
615 else:
616 msg = drv_err['evd'].format(info, n+1)
617 else:
618 msg = drv_err['evr']
620 raise LinAlgError(msg)
623_conv_dict = {0: 0, 1: 1, 2: 2,
624 'all': 0, 'value': 1, 'index': 2,
625 'a': 0, 'v': 1, 'i': 2}
628def _check_select(select, select_range, max_ev, max_len):
629 """Check that select is valid, convert to Fortran style."""
630 if isinstance(select, str):
631 select = select.lower()
632 try:
633 select = _conv_dict[select]
634 except KeyError as e:
635 raise ValueError('invalid argument for select') from e
636 vl, vu = 0., 1.
637 il = iu = 1
638 if select != 0: # (non-all)
639 sr = asarray(select_range)
640 if sr.ndim != 1 or sr.size != 2 or sr[1] < sr[0]:
641 raise ValueError('select_range must be a 2-element array-like '
642 'in nondecreasing order')
643 if select == 1: # (value)
644 vl, vu = sr
645 if max_ev == 0:
646 max_ev = max_len
647 else: # 2 (index)
648 if sr.dtype.char.lower() not in 'hilqp':
649 raise ValueError('when using select="i", select_range must '
650 'contain integers, got dtype %s (%s)'
651 % (sr.dtype, sr.dtype.char))
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].
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)
901def eigvalsh(a, b=None, lower=True, overwrite_a=False,
902 overwrite_b=False, turbo=False, eigvals=None, type=1,
903 check_finite=True, subset_by_index=None, subset_by_value=None,
904 driver=None):
905 """
906 Solves a standard or generalized eigenvalue problem for a complex
907 Hermitian or real symmetric matrix.
909 Find eigenvalues array ``w`` of array ``a``, where ``b`` is positive
910 definite such that for every eigenvalue λ (i-th entry of w) and its
911 eigenvector vi (i-th column of v) satisfies::
913 a @ vi = λ * b @ vi
914 vi.conj().T @ a @ vi = λ
915 vi.conj().T @ b @ vi = 1
917 In the standard problem, b is assumed to be the identity matrix.
919 Parameters
920 ----------
921 a : (M, M) array_like
922 A complex Hermitian or real symmetric matrix whose eigenvalues will
923 be computed.
924 b : (M, M) array_like, optional
925 A complex Hermitian or real symmetric definite positive matrix in.
926 If omitted, identity matrix is assumed.
927 lower : bool, optional
928 Whether the pertinent array data is taken from the lower or upper
929 triangle of ``a`` and, if applicable, ``b``. (Default: lower)
930 overwrite_a : bool, optional
931 Whether to overwrite data in ``a`` (may improve performance). Default
932 is False.
933 overwrite_b : bool, optional
934 Whether to overwrite data in ``b`` (may improve performance). Default
935 is False.
936 type : int, optional
937 For the generalized problems, this keyword specifies the problem type
938 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
939 inputs)::
941 1 => a @ v = w @ b @ v
942 2 => a @ b @ v = w @ v
943 3 => b @ a @ v = w @ v
945 This keyword is ignored for standard problems.
946 check_finite : bool, optional
947 Whether to check that the input matrices contain only finite numbers.
948 Disabling may give a performance gain, but may result in problems
949 (crashes, non-termination) if the inputs do contain infinities or NaNs.
950 subset_by_index : iterable, optional
951 If provided, this two-element iterable defines the start and the end
952 indices of the desired eigenvalues (ascending order and 0-indexed).
953 To return only the second smallest to fifth smallest eigenvalues,
954 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
955 available with "evr", "evx", and "gvx" drivers. The entries are
956 directly converted to integers via ``int()``.
957 subset_by_value : iterable, optional
958 If provided, this two-element iterable defines the half-open interval
959 ``(a, b]`` that, if any, only the eigenvalues between these values
960 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
961 ``np.inf`` for the unconstrained ends.
962 driver : str, optional
963 Defines which LAPACK driver should be used. Valid options are "ev",
964 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
965 generalized (where b is not None) problems. See the Notes section of
966 `scipy.linalg.eigh`.
967 turbo : bool, optional, deprecated
968 .. deprecated:: 1.5.0
969 'eigvalsh' keyword argument `turbo` is deprecated in favor of
970 ``driver=gvd`` option and will be removed in SciPy 1.12.0.
972 eigvals : tuple (lo, hi), optional
973 .. deprecated:: 1.5.0
974 'eigvalsh' keyword argument `eigvals` is deprecated in favor of
975 `subset_by_index` option and will be removed in SciPy 1.12.0.
977 Returns
978 -------
979 w : (N,) ndarray
980 The ``N`` (``1<=N<=M``) selected eigenvalues, in ascending order, each
981 repeated according to its multiplicity.
983 Raises
984 ------
985 LinAlgError
986 If eigenvalue computation does not converge, an error occurred, or
987 b matrix is not definite positive. Note that if input matrices are
988 not symmetric or Hermitian, no error will be reported but results will
989 be wrong.
991 See Also
992 --------
993 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
994 eigvals : eigenvalues of general arrays
995 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
996 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
997 matrices
999 Notes
1000 -----
1001 This function does not check the input array for being Hermitian/symmetric
1002 in order to allow for representing arrays with only their upper/lower
1003 triangular parts.
1005 This function serves as a one-liner shorthand for `scipy.linalg.eigh` with
1006 the option ``eigvals_only=True`` to get the eigenvalues and not the
1007 eigenvectors. Here it is kept as a legacy convenience. It might be
1008 beneficial to use the main function to have full control and to be a bit
1009 more pythonic.
1011 Examples
1012 --------
1013 For more examples see `scipy.linalg.eigh`.
1015 >>> import numpy as np
1016 >>> from scipy.linalg import eigvalsh
1017 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
1018 >>> w = eigvalsh(A)
1019 >>> w
1020 array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
1022 """
1023 return eigh(a, b=b, lower=lower, eigvals_only=True,
1024 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
1025 turbo=turbo, eigvals=eigvals, type=type,
1026 check_finite=check_finite, subset_by_index=subset_by_index,
1027 subset_by_value=subset_by_value, driver=driver)
1030def eigvals_banded(a_band, lower=False, overwrite_a_band=False,
1031 select='a', select_range=None, check_finite=True):
1032 """
1033 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
1035 Find eigenvalues w of a::
1037 a v[:,i] = w[i] v[:,i]
1038 v.H v = identity
1040 The matrix a is stored in a_band either in lower diagonal or upper
1041 diagonal ordered form:
1043 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
1044 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
1046 where u is the number of bands above the diagonal.
1048 Example of a_band (shape of a is (6,6), u=2)::
1050 upper form:
1051 * * a02 a13 a24 a35
1052 * a01 a12 a23 a34 a45
1053 a00 a11 a22 a33 a44 a55
1055 lower form:
1056 a00 a11 a22 a33 a44 a55
1057 a10 a21 a32 a43 a54 *
1058 a20 a31 a42 a53 * *
1060 Cells marked with * are not used.
1062 Parameters
1063 ----------
1064 a_band : (u+1, M) array_like
1065 The bands of the M by M matrix a.
1066 lower : bool, optional
1067 Is the matrix in the lower form. (Default is upper form)
1068 overwrite_a_band : bool, optional
1069 Discard data in a_band (may enhance performance)
1070 select : {'a', 'v', 'i'}, optional
1071 Which eigenvalues to calculate
1073 ====== ========================================
1074 select calculated
1075 ====== ========================================
1076 'a' All eigenvalues
1077 'v' Eigenvalues in the interval (min, max]
1078 'i' Eigenvalues with indices min <= i <= max
1079 ====== ========================================
1080 select_range : (min, max), optional
1081 Range of selected eigenvalues
1082 check_finite : bool, optional
1083 Whether to check that the input matrix contains only finite numbers.
1084 Disabling may give a performance gain, but may result in problems
1085 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1087 Returns
1088 -------
1089 w : (M,) ndarray
1090 The eigenvalues, in ascending order, each repeated according to its
1091 multiplicity.
1093 Raises
1094 ------
1095 LinAlgError
1096 If eigenvalue computation does not converge.
1098 See Also
1099 --------
1100 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1101 band matrices
1102 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1103 matrices
1104 eigvals : eigenvalues of general arrays
1105 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1106 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1108 Examples
1109 --------
1110 >>> import numpy as np
1111 >>> from scipy.linalg import eigvals_banded
1112 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
1113 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
1114 >>> w = eigvals_banded(Ab, lower=True)
1115 >>> w
1116 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
1117 """
1118 return eig_banded(a_band, lower=lower, eigvals_only=1,
1119 overwrite_a_band=overwrite_a_band, select=select,
1120 select_range=select_range, check_finite=check_finite)
1123def eigvalsh_tridiagonal(d, e, select='a', select_range=None,
1124 check_finite=True, tol=0., lapack_driver='auto'):
1125 """
1126 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1128 Find eigenvalues `w` of ``a``::
1130 a v[:,i] = w[i] v[:,i]
1131 v.H v = identity
1133 For a real symmetric matrix ``a`` with diagonal elements `d` and
1134 off-diagonal elements `e`.
1136 Parameters
1137 ----------
1138 d : ndarray, shape (ndim,)
1139 The diagonal elements of the array.
1140 e : ndarray, shape (ndim-1,)
1141 The off-diagonal elements of the array.
1142 select : {'a', 'v', 'i'}, optional
1143 Which eigenvalues to calculate
1145 ====== ========================================
1146 select calculated
1147 ====== ========================================
1148 'a' All eigenvalues
1149 'v' Eigenvalues in the interval (min, max]
1150 'i' Eigenvalues with indices min <= i <= max
1151 ====== ========================================
1152 select_range : (min, max), optional
1153 Range of selected eigenvalues
1154 check_finite : bool, optional
1155 Whether to check that the input matrix contains only finite numbers.
1156 Disabling may give a performance gain, but may result in problems
1157 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1158 tol : float
1159 The absolute tolerance to which each eigenvalue is required
1160 (only used when ``lapack_driver='stebz'``).
1161 An eigenvalue (or cluster) is considered to have converged if it
1162 lies in an interval of this width. If <= 0. (default),
1163 the value ``eps*|a|`` is used where eps is the machine precision,
1164 and ``|a|`` is the 1-norm of the matrix ``a``.
1165 lapack_driver : str
1166 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1167 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1168 and 'stebz' otherwise. 'sterf' and 'stev' can only be used when
1169 ``select='a'``.
1171 Returns
1172 -------
1173 w : (M,) ndarray
1174 The eigenvalues, in ascending order, each repeated according to its
1175 multiplicity.
1177 Raises
1178 ------
1179 LinAlgError
1180 If eigenvalue computation does not converge.
1182 See Also
1183 --------
1184 eigh_tridiagonal : eigenvalues and right eiegenvectors for
1185 symmetric/Hermitian tridiagonal matrices
1187 Examples
1188 --------
1189 >>> import numpy as np
1190 >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
1191 >>> d = 3*np.ones(4)
1192 >>> e = -1*np.ones(3)
1193 >>> w = eigvalsh_tridiagonal(d, e)
1194 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1195 >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
1196 >>> np.allclose(w - w2, np.zeros(4))
1197 True
1198 """
1199 return eigh_tridiagonal(
1200 d, e, eigvals_only=True, select=select, select_range=select_range,
1201 check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
1204def eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None,
1205 check_finite=True, tol=0., lapack_driver='auto'):
1206 """
1207 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1209 Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
1211 a v[:,i] = w[i] v[:,i]
1212 v.H v = identity
1214 For a real symmetric matrix ``a`` with diagonal elements `d` and
1215 off-diagonal elements `e`.
1217 Parameters
1218 ----------
1219 d : ndarray, shape (ndim,)
1220 The diagonal elements of the array.
1221 e : ndarray, shape (ndim-1,)
1222 The off-diagonal elements of the array.
1223 select : {'a', 'v', 'i'}, optional
1224 Which eigenvalues to calculate
1226 ====== ========================================
1227 select calculated
1228 ====== ========================================
1229 'a' All eigenvalues
1230 'v' Eigenvalues in the interval (min, max]
1231 'i' Eigenvalues with indices min <= i <= max
1232 ====== ========================================
1233 select_range : (min, max), optional
1234 Range of selected eigenvalues
1235 check_finite : bool, optional
1236 Whether to check that the input matrix contains only finite numbers.
1237 Disabling may give a performance gain, but may result in problems
1238 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1239 tol : float
1240 The absolute tolerance to which each eigenvalue is required
1241 (only used when 'stebz' is the `lapack_driver`).
1242 An eigenvalue (or cluster) is considered to have converged if it
1243 lies in an interval of this width. If <= 0. (default),
1244 the value ``eps*|a|`` is used where eps is the machine precision,
1245 and ``|a|`` is the 1-norm of the matrix ``a``.
1246 lapack_driver : str
1247 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1248 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1249 and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and
1250 ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is
1251 used to find the corresponding eigenvectors. 'sterf' can only be
1252 used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only
1253 be used when ``select='a'``.
1255 Returns
1256 -------
1257 w : (M,) ndarray
1258 The eigenvalues, in ascending order, each repeated according to its
1259 multiplicity.
1260 v : (M, M) ndarray
1261 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
1262 the column ``v[:,i]``.
1264 Raises
1265 ------
1266 LinAlgError
1267 If eigenvalue computation does not converge.
1269 See Also
1270 --------
1271 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1272 matrices
1273 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1274 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1275 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1276 band matrices
1278 Notes
1279 -----
1280 This function makes use of LAPACK ``S/DSTEMR`` routines.
1282 Examples
1283 --------
1284 >>> import numpy as np
1285 >>> from scipy.linalg import eigh_tridiagonal
1286 >>> d = 3*np.ones(4)
1287 >>> e = -1*np.ones(3)
1288 >>> w, v = eigh_tridiagonal(d, e)
1289 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1290 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
1291 True
1292 """
1293 d = _asarray_validated(d, check_finite=check_finite)
1294 e = _asarray_validated(e, check_finite=check_finite)
1295 for check in (d, e):
1296 if check.ndim != 1:
1297 raise ValueError('expected a 1-D array')
1298 if check.dtype.char in 'GFD': # complex
1299 raise TypeError('Only real arrays currently supported')
1300 if d.size != e.size + 1:
1301 raise ValueError('d (%s) must have one more element than e (%s)'
1302 % (d.size, e.size))
1303 select, vl, vu, il, iu, _ = _check_select(
1304 select, select_range, 0, d.size)
1305 if not isinstance(lapack_driver, str):
1306 raise TypeError('lapack_driver must be str')
1307 drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev')
1308 if lapack_driver not in drivers:
1309 raise ValueError('lapack_driver must be one of %s, got %s'
1310 % (drivers, lapack_driver))
1311 if lapack_driver == 'auto':
1312 lapack_driver = 'stemr' if select == 0 else 'stebz'
1313 func, = get_lapack_funcs((lapack_driver,), (d, e))
1314 compute_v = not eigvals_only
1315 if lapack_driver == 'sterf':
1316 if select != 0:
1317 raise ValueError('sterf can only be used when select == "a"')
1318 if not eigvals_only:
1319 raise ValueError('sterf can only be used when eigvals_only is '
1320 'True')
1321 w, info = func(d, e)
1322 m = len(w)
1323 elif lapack_driver == 'stev':
1324 if select != 0:
1325 raise ValueError('stev can only be used when select == "a"')
1326 w, v, info = func(d, e, compute_v=compute_v)
1327 m = len(w)
1328 elif lapack_driver == 'stebz':
1329 tol = float(tol)
1330 internal_name = 'stebz'
1331 stebz, = get_lapack_funcs((internal_name,), (d, e))
1332 # If getting eigenvectors, needs to be block-ordered (B) instead of
1333 # matrix-ordered (E), and we will reorder later
1334 order = 'E' if eigvals_only else 'B'
1335 m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol,
1336 order)
1337 else: # 'stemr'
1338 # ?STEMR annoyingly requires size N instead of N-1
1339 e_ = empty(e.size+1, e.dtype)
1340 e_[:-1] = e
1341 stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e))
1342 lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu,
1343 compute_v=compute_v)
1344 _check_info(info, 'stemr_lwork')
1345 m, w, v, info = func(d, e_, select, vl, vu, il, iu,
1346 compute_v=compute_v, lwork=lwork, liwork=liwork)
1347 _check_info(info, lapack_driver + ' (eigh_tridiagonal)')
1348 w = w[:m]
1349 if eigvals_only:
1350 return w
1351 else:
1352 # Do we still need to compute the eigenvalues?
1353 if lapack_driver == 'stebz':
1354 func, = get_lapack_funcs(('stein',), (d, e))
1355 v, info = func(d, e, w, iblock, isplit)
1356 _check_info(info, 'stein (eigh_tridiagonal)',
1357 positive='%d eigenvectors failed to converge')
1358 # Convert block-order to matrix-order
1359 order = argsort(w)
1360 w, v = w[order], v[:, order]
1361 else:
1362 v = v[:, :m]
1363 return w, v
1366def _check_info(info, driver, positive='did not converge (LAPACK info=%d)'):
1367 """Check info return value."""
1368 if info < 0:
1369 raise ValueError('illegal value in argument %d of internal %s'
1370 % (-info, driver))
1371 if info > 0 and positive:
1372 raise LinAlgError(("%s " + positive) % (driver, info,))
1375def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
1376 """
1377 Compute Hessenberg form of a matrix.
1379 The Hessenberg decomposition is::
1381 A = Q H Q^H
1383 where `Q` is unitary/orthogonal and `H` has only zero elements below
1384 the first sub-diagonal.
1386 Parameters
1387 ----------
1388 a : (M, M) array_like
1389 Matrix to bring into Hessenberg form.
1390 calc_q : bool, optional
1391 Whether to compute the transformation matrix. Default is False.
1392 overwrite_a : bool, optional
1393 Whether to overwrite `a`; may improve performance.
1394 Default is False.
1395 check_finite : bool, optional
1396 Whether to check that the input matrix contains only finite numbers.
1397 Disabling may give a performance gain, but may result in problems
1398 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1400 Returns
1401 -------
1402 H : (M, M) ndarray
1403 Hessenberg form of `a`.
1404 Q : (M, M) ndarray
1405 Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
1406 Only returned if ``calc_q=True``.
1408 Examples
1409 --------
1410 >>> import numpy as np
1411 >>> from scipy.linalg import hessenberg
1412 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
1413 >>> H, Q = hessenberg(A, calc_q=True)
1414 >>> H
1415 array([[ 2. , -11.65843866, 1.42005301, 0.25349066],
1416 [ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
1417 [ 0. , -1.83299243, 0.38969961, -0.51527034],
1418 [ 0. , 0. , -3.83189513, 1.07494686]])
1419 >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
1420 True
1421 """
1422 a1 = _asarray_validated(a, check_finite=check_finite)
1423 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
1424 raise ValueError('expected square matrix')
1425 overwrite_a = overwrite_a or (_datacopied(a1, a))
1427 # if 2x2 or smaller: already in Hessenberg
1428 if a1.shape[0] <= 2:
1429 if calc_q:
1430 return a1, eye(a1.shape[0])
1431 return a1
1433 gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal',
1434 'gehrd_lwork'), (a1,))
1435 ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a)
1436 _check_info(info, 'gebal (hessenberg)', positive=False)
1437 n = len(a1)
1439 lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
1441 hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1442 _check_info(info, 'gehrd (hessenberg)', positive=False)
1443 h = numpy.triu(hq, -1)
1444 if not calc_q:
1445 return h
1447 # use orghr/unghr to compute q
1448 orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,))
1449 lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
1451 q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1452 _check_info(info, 'orghr (hessenberg)', positive=False)
1453 return h, q
1456def cdf2rdf(w, v):
1457 """
1458 Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real
1459 eigenvalues in a block diagonal form ``wr`` and the associated real
1460 eigenvectors ``vr``, such that::
1462 vr @ wr = X @ vr
1464 continues to hold, where ``X`` is the original array for which ``w`` and
1465 ``v`` are the eigenvalues and eigenvectors.
1467 .. versionadded:: 1.1.0
1469 Parameters
1470 ----------
1471 w : (..., M) array_like
1472 Complex or real eigenvalues, an array or stack of arrays
1474 Conjugate pairs must not be interleaved, else the wrong result
1475 will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result,
1476 but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
1478 v : (..., M, M) array_like
1479 Complex or real eigenvectors, a square array or stack of square arrays.
1481 Returns
1482 -------
1483 wr : (..., M, M) ndarray
1484 Real diagonal block form of eigenvalues
1485 vr : (..., M, M) ndarray
1486 Real eigenvectors associated with ``wr``
1488 See Also
1489 --------
1490 eig : Eigenvalues and right eigenvectors for non-symmetric arrays
1491 rsf2csf : Convert real Schur form to complex Schur form
1493 Notes
1494 -----
1495 ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``.
1496 For example, obtained by ``w, v = scipy.linalg.eig(X)`` or
1497 ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent
1498 stacked arrays.
1500 .. versionadded:: 1.1.0
1502 Examples
1503 --------
1504 >>> import numpy as np
1505 >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
1506 >>> X
1507 array([[ 1, 2, 3],
1508 [ 0, 4, 5],
1509 [ 0, -5, 4]])
1511 >>> from scipy import linalg
1512 >>> w, v = linalg.eig(X)
1513 >>> w
1514 array([ 1.+0.j, 4.+5.j, 4.-5.j])
1515 >>> v
1516 array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
1517 [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
1518 [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
1520 >>> wr, vr = linalg.cdf2rdf(w, v)
1521 >>> wr
1522 array([[ 1., 0., 0.],
1523 [ 0., 4., 5.],
1524 [ 0., -5., 4.]])
1525 >>> vr
1526 array([[ 1. , 0.40016, -0.01906],
1527 [ 0. , 0.64788, 0. ],
1528 [ 0. , 0. , 0.64788]])
1530 >>> vr @ wr
1531 array([[ 1. , 1.69593, 1.9246 ],
1532 [ 0. , 2.59153, 3.23942],
1533 [ 0. , -3.23942, 2.59153]])
1534 >>> X @ vr
1535 array([[ 1. , 1.69593, 1.9246 ],
1536 [ 0. , 2.59153, 3.23942],
1537 [ 0. , -3.23942, 2.59153]])
1538 """
1539 w, v = _asarray_validated(w), _asarray_validated(v)
1541 # check dimensions
1542 if w.ndim < 1:
1543 raise ValueError('expected w to be at least 1D')
1544 if v.ndim < 2:
1545 raise ValueError('expected v to be at least 2D')
1546 if v.ndim != w.ndim + 1:
1547 raise ValueError('expected eigenvectors array to have exactly one '
1548 'dimension more than eigenvalues array')
1550 # check shapes
1551 n = w.shape[-1]
1552 M = w.shape[:-1]
1553 if v.shape[-2] != v.shape[-1]:
1554 raise ValueError('expected v to be a square matrix or stacked square '
1555 'matrices: v.shape[-2] = v.shape[-1]')
1556 if v.shape[-1] != n:
1557 raise ValueError('expected the same number of eigenvalues as '
1558 'eigenvectors')
1560 # get indices for each first pair of complex eigenvalues
1561 complex_mask = iscomplex(w)
1562 n_complex = complex_mask.sum(axis=-1)
1564 # check if all complex eigenvalues have conjugate pairs
1565 if not (n_complex % 2 == 0).all():
1566 raise ValueError('expected complex-conjugate pairs of eigenvalues')
1568 # find complex indices
1569 idx = nonzero(complex_mask)
1570 idx_stack = idx[:-1]
1571 idx_elem = idx[-1]
1573 # filter them to conjugate indices, assuming pairs are not interleaved
1574 j = idx_elem[0::2]
1575 k = idx_elem[1::2]
1576 stack_ind = ()
1577 for i in idx_stack:
1578 # should never happen, assuming nonzero orders by the last axis
1579 assert (i[0::2] == i[1::2]).all(),\
1580 "Conjugate pair spanned different arrays!"
1581 stack_ind += (i[0::2],)
1583 # all eigenvalues to diagonal form
1584 wr = zeros(M + (n, n), dtype=w.real.dtype)
1585 di = range(n)
1586 wr[..., di, di] = w.real
1588 # complex eigenvalues to real block diagonal form
1589 wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag
1590 wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
1592 # compute real eigenvectors associated with real block diagonal eigenvalues
1593 u = zeros(M + (n, n), dtype=numpy.cdouble)
1594 u[..., di, di] = 1.0
1595 u[stack_ind + (j, j)] = 0.5j
1596 u[stack_ind + (j, k)] = 0.5
1597 u[stack_ind + (k, j)] = -0.5j
1598 u[stack_ind + (k, k)] = 0.5
1600 # multipy matrices v and u (equivalent to v @ u)
1601 vr = einsum('...ij,...jk->...ik', v, u).real
1603 return wr, vr