Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/linalg/decomp.py : 6%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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 numpy
21from numpy import (array, isfinite, inexact, nonzero, iscomplexobj, cast,
22 flatnonzero, conj, asarray, argsort, empty,
23 iscomplex, zeros, einsum, eye, inf)
24# Local imports
25from scipy._lib._util import _asarray_validated
26from .misc import LinAlgError, _datacopied, norm
27from .lapack import get_lapack_funcs, _compute_lwork
30_I = cast['F'](1j)
33def _make_complex_eigvecs(w, vin, dtype):
34 """
35 Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output
36 """
37 # - see LAPACK man page DGGEV at ALPHAI
38 v = numpy.array(vin, dtype=dtype)
39 m = (w.imag > 0)
40 m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709
41 for i in flatnonzero(m):
42 v.imag[:, i] = vin[:, i+1]
43 conj(v[:, i], v[:, i+1])
44 return v
47def _make_eigvals(alpha, beta, homogeneous_eigvals):
48 if homogeneous_eigvals:
49 if beta is None:
50 return numpy.vstack((alpha, numpy.ones_like(alpha)))
51 else:
52 return numpy.vstack((alpha, beta))
53 else:
54 if beta is None:
55 return alpha
56 else:
57 w = numpy.empty_like(alpha)
58 alpha_zero = (alpha == 0)
59 beta_zero = (beta == 0)
60 beta_nonzero = ~beta_zero
61 w[beta_nonzero] = alpha[beta_nonzero]/beta[beta_nonzero]
62 # Use numpy.inf for complex values too since
63 # 1/numpy.inf = 0, i.e., it correctly behaves as projective
64 # infinity.
65 w[~alpha_zero & beta_zero] = numpy.inf
66 if numpy.all(alpha.imag == 0):
67 w[alpha_zero & beta_zero] = numpy.nan
68 else:
69 w[alpha_zero & beta_zero] = complex(numpy.nan, numpy.nan)
70 return w
73def _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
74 homogeneous_eigvals):
75 ggev, = get_lapack_funcs(('ggev',), (a1, b1))
76 cvl, cvr = left, right
77 res = ggev(a1, b1, lwork=-1)
78 lwork = res[-2][0].real.astype(numpy.int)
79 if ggev.typecode in 'cz':
80 alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
81 overwrite_a, overwrite_b)
82 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
83 else:
84 alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr,
85 lwork, overwrite_a,
86 overwrite_b)
87 alpha = alphar + _I * alphai
88 w = _make_eigvals(alpha, beta, homogeneous_eigvals)
89 _check_info(info, 'generalized eig algorithm (ggev)')
91 only_real = numpy.all(w.imag == 0.0)
92 if not (ggev.typecode in 'cz' or only_real):
93 t = w.dtype.char
94 if left:
95 vl = _make_complex_eigvecs(w, vl, t)
96 if right:
97 vr = _make_complex_eigvecs(w, vr, t)
99 # the eigenvectors returned by the lapack function are NOT normalized
100 for i in range(vr.shape[0]):
101 if right:
102 vr[:, i] /= norm(vr[:, i])
103 if left:
104 vl[:, i] /= norm(vl[:, i])
106 if not (left or right):
107 return w
108 if left:
109 if right:
110 return w, vl, vr
111 return w, vl
112 return w, vr
115def eig(a, b=None, left=False, right=True, overwrite_a=False,
116 overwrite_b=False, check_finite=True, homogeneous_eigvals=False):
117 """
118 Solve an ordinary or generalized eigenvalue problem of a square matrix.
120 Find eigenvalues w and right or left eigenvectors of a general matrix::
122 a vr[:,i] = w[i] b vr[:,i]
123 a.H vl[:,i] = w[i].conj() b.H vl[:,i]
125 where ``.H`` is the Hermitian conjugation.
127 Parameters
128 ----------
129 a : (M, M) array_like
130 A complex or real matrix whose eigenvalues and eigenvectors
131 will be computed.
132 b : (M, M) array_like, optional
133 Right-hand side matrix in a generalized eigenvalue problem.
134 Default is None, identity matrix is assumed.
135 left : bool, optional
136 Whether to calculate and return left eigenvectors. Default is False.
137 right : bool, optional
138 Whether to calculate and return right eigenvectors. Default is True.
139 overwrite_a : bool, optional
140 Whether to overwrite `a`; may improve performance. Default is False.
141 overwrite_b : bool, optional
142 Whether to overwrite `b`; may improve performance. Default is False.
143 check_finite : bool, optional
144 Whether to check that the input matrices contain only finite numbers.
145 Disabling may give a performance gain, but may result in problems
146 (crashes, non-termination) if the inputs do contain infinities or NaNs.
147 homogeneous_eigvals : bool, optional
148 If True, return the eigenvalues in homogeneous coordinates.
149 In this case ``w`` is a (2, M) array so that::
151 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
153 Default is False.
155 Returns
156 -------
157 w : (M,) or (2, M) double or complex ndarray
158 The eigenvalues, each repeated according to its
159 multiplicity. The shape is (M,) unless
160 ``homogeneous_eigvals=True``.
161 vl : (M, M) double or complex ndarray
162 The normalized left eigenvector corresponding to the eigenvalue
163 ``w[i]`` is the column vl[:,i]. Only returned if ``left=True``.
164 vr : (M, M) double or complex ndarray
165 The normalized right eigenvector corresponding to the eigenvalue
166 ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
168 Raises
169 ------
170 LinAlgError
171 If eigenvalue computation does not converge.
173 See Also
174 --------
175 eigvals : eigenvalues of general arrays
176 eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
177 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
178 band matrices
179 eigh_tridiagonal : eigenvalues and right eiegenvectors for
180 symmetric/Hermitian tridiagonal matrices
182 Examples
183 --------
184 >>> from scipy import linalg
185 >>> a = np.array([[0., -1.], [1., 0.]])
186 >>> linalg.eigvals(a)
187 array([0.+1.j, 0.-1.j])
189 >>> b = np.array([[0., 1.], [1., 1.]])
190 >>> linalg.eigvals(a, b)
191 array([ 1.+0.j, -1.+0.j])
193 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
194 >>> linalg.eigvals(a, homogeneous_eigvals=True)
195 array([[3.+0.j, 8.+0.j, 7.+0.j],
196 [1.+0.j, 1.+0.j, 1.+0.j]])
198 >>> a = np.array([[0., -1.], [1., 0.]])
199 >>> linalg.eigvals(a) == linalg.eig(a)[0]
200 array([ True, True])
201 >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector
202 array([[-0.70710678+0.j , -0.70710678-0.j ],
203 [-0. +0.70710678j, -0. -0.70710678j]])
204 >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector
205 array([[0.70710678+0.j , 0.70710678-0.j ],
206 [0. -0.70710678j, 0. +0.70710678j]])
210 """
211 a1 = _asarray_validated(a, check_finite=check_finite)
212 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
213 raise ValueError('expected square matrix')
214 overwrite_a = overwrite_a or (_datacopied(a1, a))
215 if b is not None:
216 b1 = _asarray_validated(b, check_finite=check_finite)
217 overwrite_b = overwrite_b or _datacopied(b1, b)
218 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
219 raise ValueError('expected square matrix')
220 if b1.shape != a1.shape:
221 raise ValueError('a and b must have the same shape')
222 return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
223 homogeneous_eigvals)
225 geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
226 compute_vl, compute_vr = left, right
228 lwork = _compute_lwork(geev_lwork, a1.shape[0],
229 compute_vl=compute_vl,
230 compute_vr=compute_vr)
232 if geev.typecode in 'cz':
233 w, vl, vr, info = geev(a1, lwork=lwork,
234 compute_vl=compute_vl,
235 compute_vr=compute_vr,
236 overwrite_a=overwrite_a)
237 w = _make_eigvals(w, None, homogeneous_eigvals)
238 else:
239 wr, wi, vl, vr, info = geev(a1, lwork=lwork,
240 compute_vl=compute_vl,
241 compute_vr=compute_vr,
242 overwrite_a=overwrite_a)
243 t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
244 w = wr + _I * wi
245 w = _make_eigvals(w, None, homogeneous_eigvals)
247 _check_info(info, 'eig algorithm (geev)',
248 positive='did not converge (only eigenvalues '
249 'with order >= %d have converged)')
251 only_real = numpy.all(w.imag == 0.0)
252 if not (geev.typecode in 'cz' or only_real):
253 t = w.dtype.char
254 if left:
255 vl = _make_complex_eigvecs(w, vl, t)
256 if right:
257 vr = _make_complex_eigvecs(w, vr, t)
258 if not (left or right):
259 return w
260 if left:
261 if right:
262 return w, vl, vr
263 return w, vl
264 return w, vr
267def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
268 overwrite_b=False, turbo=True, eigvals=None, type=1,
269 check_finite=True, subset_by_index=None, subset_by_value=None,
270 driver=None):
271 """
272 Solve a standard or generalized eigenvalue problem for a complex
273 Hermitian or real symmetric matrix.
275 Find eigenvalues array ``w`` and optionally eigenvectors array ``v`` of
276 array ``a``, where ``b`` is positive definite such that for every
277 eigenvalue λ (i-th entry of w) and its eigenvector ``vi`` (i-th column of
278 ``v``) satisfies::
280 a @ vi = λ * b @ vi
281 vi.conj().T @ a @ vi = λ
282 vi.conj().T @ b @ vi = 1
284 In the standard problem, ``b`` is assumed to be the identity matrix.
286 Parameters
287 ----------
288 a : (M, M) array_like
289 A complex Hermitian or real symmetric matrix whose eigenvalues and
290 eigenvectors will be computed.
291 b : (M, M) array_like, optional
292 A complex Hermitian or real symmetric definite positive matrix in.
293 If omitted, identity matrix is assumed.
294 lower : bool, optional
295 Whether the pertinent array data is taken from the lower or upper
296 triangle of `a`. (Default: lower)
297 eigvals_only : bool, optional
298 Whether to calculate only eigenvalues and no eigenvectors.
299 (Default: both are calculated)
300 subset_by_index : iterable, optional
301 If provided, this two-element iterable defines the start and the end
302 indices of the desired eigenvalues (ascending order and 0-indexed).
303 To return only the second smallest to fifth smallest eigenvalues,
304 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
305 available with "evr", "evx", and "gvx" drivers. The entries are
306 directly converted to integers via ``int()``.
307 subset_by_value : iterable, optional
308 If provided, this two-element iterable defines the half-open interval
309 ``(a, b]`` that, if any, only the eigenvalues between these values
310 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
311 ``np.inf`` for the unconstrained ends.
312 driver: str, optional
313 Defines which LAPACK driver should be used. Valid options are "ev",
314 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
315 generalized (where b is not None) problems. See the Notes section.
316 type : int, optional
317 For the generalized problems, this keyword specifies the problem type
318 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
319 inputs)::
321 1 => a @ v = w @ b @ v
322 2 => a @ b @ v = w @ v
323 3 => b @ a @ v = w @ v
325 This keyword is ignored for standard problems.
326 overwrite_a : bool, optional
327 Whether to overwrite data in ``a`` (may improve performance). Default
328 is False.
329 overwrite_b : bool, optional
330 Whether to overwrite data in ``b`` (may improve performance). Default
331 is False.
332 check_finite : bool, optional
333 Whether to check that the input matrices contain only finite numbers.
334 Disabling may give a performance gain, but may result in problems
335 (crashes, non-termination) if the inputs do contain infinities or NaNs.
336 turbo : bool, optional
337 *Deprecated since v1.5.0, use ``driver=gvd`` keyword instead*.
338 Use divide and conquer algorithm (faster but expensive in memory, only
339 for generalized eigenvalue problem and if full set of eigenvalues are
340 requested.). Has no significant effect if eigenvectors are not
341 requested.
342 eigvals : tuple (lo, hi), optional
343 *Deprecated since v1.5.0, use ``subset_by_index`` keyword instead*.
344 Indexes of the smallest and largest (in ascending order) eigenvalues
345 and corresponding eigenvectors to be returned: 0 <= lo <= hi <= M-1.
346 If omitted, all eigenvalues and eigenvectors are returned.
348 Returns
349 -------
350 w : (N,) ndarray
351 The N (1<=N<=M) selected eigenvalues, in ascending order, each
352 repeated according to its multiplicity.
353 v : (M, N) ndarray
354 (if ``eigvals_only == False``)
356 Raises
357 ------
358 LinAlgError
359 If eigenvalue computation does not converge, an error occurred, or
360 b matrix is not definite positive. Note that if input matrices are
361 not symmetric or Hermitian, no error will be reported but results will
362 be wrong.
364 See Also
365 --------
366 eigvalsh : eigenvalues of symmetric or Hermitian arrays
367 eig : eigenvalues and right eigenvectors for non-symmetric arrays
368 eigh_tridiagonal : eigenvalues and right eiegenvectors for
369 symmetric/Hermitian tridiagonal matrices
371 Notes
372 -----
373 This function does not check the input array for being hermitian/symmetric
374 in order to allow for representing arrays with only their upper/lower
375 triangular parts. Also, note that even though not taken into account,
376 finiteness check applies to the whole array and unaffected by "lower"
377 keyword.
379 This function uses LAPACK drivers for computations in all possible keyword
380 combinations, prefixed with ``sy`` if arrays are real and ``he`` if
381 complex, e.g., a float array with "evr" driver is solved via
382 "syevr", complex arrays with "gvx" driver problem is solved via "hegvx"
383 etc.
385 As a brief summary, the slowest and the most robust driver is the
386 classical ``<sy/he>ev`` which uses symmetric QR. ``<sy/he>evr`` is seen as
387 the optimal choice for the most general cases. However, there are certain
388 occassions that ``<sy/he>evd`` computes faster at the expense of more
389 memory usage. ``<sy/he>evx``, while still being faster than ``<sy/he>ev``,
390 often performs worse than the rest except when very few eigenvalues are
391 requested for large arrays though there is still no performance guarantee.
394 For the generalized problem, normalization with respoect to the given
395 type argument::
397 type 1 and 3 : v.conj().T @ a @ v = w
398 type 2 : inv(v).conj().T @ a @ inv(v) = w
400 type 1 or 2 : v.conj().T @ b @ v = I
401 type 3 : v.conj().T @ inv(b) @ v = I
404 Examples
405 --------
406 >>> from scipy.linalg import eigh
407 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
408 >>> w, v = eigh(A)
409 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
410 True
412 Request only the eigenvalues
414 >>> w = eigh(A, eigvals_only=True)
416 Request eigenvalues that are less than 10.
418 >>> A = np.array([[34, -4, -10, -7, 2],
419 ... [-4, 7, 2, 12, 0],
420 ... [-10, 2, 44, 2, -19],
421 ... [-7, 12, 2, 79, -34],
422 ... [2, 0, -19, -34, 29]])
423 >>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
424 array([6.69199443e-07, 9.11938152e+00])
426 Request the largest second eigenvalue and its eigenvector
428 >>> w, v = eigh(A, subset_by_index=[1, 1])
429 >>> w
430 array([9.11938152])
431 >>> v.shape # only a single column is returned
432 (5, 1)
434 """
435 # set lower
436 uplo = 'L' if lower else 'U'
437 # Set job for Fortran routines
438 _job = 'N' if eigvals_only else 'V'
440 drv_str = [None, "ev", "evd", "evr", "evx", "gv", "gvd", "gvx"]
441 if driver not in drv_str:
442 raise ValueError('"{}" is unknown. Possible values are "None", "{}".'
443 ''.format(driver, '", "'.join(drv_str[1:])))
445 a1 = _asarray_validated(a, check_finite=check_finite)
446 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
447 raise ValueError('expected square "a" matrix')
448 overwrite_a = overwrite_a or (_datacopied(a1, a))
449 cplx = True if iscomplexobj(a1) else False
450 n = a1.shape[0]
451 drv_args = {'overwrite_a': overwrite_a}
453 if b is not None:
454 b1 = _asarray_validated(b, check_finite=check_finite)
455 overwrite_b = overwrite_b or _datacopied(b1, b)
456 if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
457 raise ValueError('expected square "b" matrix')
459 if b1.shape != a1.shape:
460 raise ValueError("wrong b dimensions {}, should "
461 "be {}".format(b1.shape, a1.shape))
463 cplx = True if iscomplexobj(b1) else (cplx or False)
464 drv_args.update({'overwrite_b': overwrite_b})
466 # backwards-compatibility handling
467 subset_by_index = subset_by_index if (eigvals is None) else eigvals
469 subset = (subset_by_index is not None) or (subset_by_value is not None)
471 # Both subsets can't be given
472 if subset_by_index and subset_by_value:
473 raise ValueError('Either index or value subset can be requested.')
475 # Take turbo into account if all conditions are met otherwise ignore
476 if turbo and b is not None:
477 driver = 'gvx' if subset else 'gvd'
479 # Check indices if given
480 if subset_by_index:
481 lo, hi = [int(x) for x in subset_by_index]
482 if not (0 <= lo <= hi < n):
483 raise ValueError('Requested eigenvalue indices are not valid. '
484 'Valid range is [0, {}] and start <= end, but '
485 'start={}, end={} is given'.format(n-1, lo, hi))
486 # fortran is 1-indexed
487 drv_args.update({'range': 'I', 'il': lo + 1, 'iu': hi + 1})
489 if subset_by_value:
490 lo, hi = subset_by_value
491 if not (-inf <= lo < hi <= inf):
492 raise ValueError('Requested eigenvalue bounds are not valid. '
493 'Valid range is (-inf, inf) and low < high, but '
494 'low={}, high={} is given'.format(lo, hi))
496 drv_args.update({'range': 'V', 'vl': lo, 'vu': hi})
498 # fix prefix for lapack routines
499 pfx = 'he' if cplx else 'sy'
501 # decide on the driver if not given
502 # first early exit on incompatible choice
503 if driver:
504 if b is None and (driver in ["gv", "gvd", "gvx"]):
505 raise ValueError('{} requires input b array to be supplied '
506 'for generalized eigenvalue problems.'
507 ''.format(driver))
508 if (b is not None) and (driver in ['ev', 'evd', 'evr', 'evx']):
509 raise ValueError('"{}" does not accept input b array '
510 'for standard eigenvalue problems.'
511 ''.format(driver))
512 if subset and (driver in ["ev", "evd", "gv", "gvd"]):
513 raise ValueError('"{}" cannot compute subsets of eigenvalues'
514 ''.format(driver))
516 # Default driver is evr and gvd
517 else:
518 driver = "evr" if b is None else ("gvx" if subset else "gvd")
520 lwork_spec = {
521 'syevd': ['lwork', 'liwork'],
522 'syevr': ['lwork', 'liwork'],
523 'heevd': ['lwork', 'liwork', 'lrwork'],
524 'heevr': ['lwork', 'lrwork', 'liwork'],
525 }
527 if b is None: # Standard problem
528 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
529 [a1])
530 clw_args = {'n': n, 'lower': lower}
531 if driver == 'evd':
532 clw_args.update({'compute_v': 0 if _job == "N" else 1})
534 lw = _compute_lwork(drvlw, **clw_args)
535 # Multiple lwork vars
536 if isinstance(lw, tuple):
537 lwork_args = dict(zip(lwork_spec[pfx+driver], lw))
538 else:
539 lwork_args = {'lwork': lw}
541 drv_args.update({'lower': lower, 'compute_v': 0 if _job == "N" else 1})
542 w, v, *other_args, info = drv(a=a1, **drv_args, **lwork_args)
544 else: # Generalized problem
545 # 'gvd' doesn't have lwork query
546 if driver == "gvd":
547 drv = get_lapack_funcs(pfx + "gvd", [a1, b1])
548 lwork_args = {}
549 else:
550 drv, drvlw = get_lapack_funcs((pfx + driver, pfx+driver+'_lwork'),
551 [a1, b1])
552 # generalized drivers use uplo instead of lower
553 lw = _compute_lwork(drvlw, n, uplo=uplo)
554 lwork_args = {'lwork': lw}
556 drv_args.update({'uplo': uplo, 'jobz': _job})
558 w, v, *other_args, info = drv(a=a1, b=b1, **drv_args, **lwork_args)
560 # m is always the first extra argument
561 w = w[:other_args[0]] if subset else w
562 v = v[:, :other_args[0]] if (subset and not eigvals_only) else v
564 # Check if we had a successful exit
565 if info == 0:
566 if eigvals_only:
567 return w
568 else:
569 return w, v
570 else:
571 if info < -1:
572 raise LinAlgError('Illegal value in argument {} of internal {}'
573 ''.format(-info, drv.typecode + pfx + driver))
574 elif info > n:
575 raise LinAlgError('The leading minor of order {} of B is not '
576 'positive definite. The factorization of B '
577 'could not be completed and no eigenvalues '
578 'or eigenvectors were computed.'.format(info-n))
579 else:
580 drv_err = {'ev': 'The algorithm failed to converge; {} '
581 'off-diagonal elements of an intermediate '
582 'tridiagonal form did not converge to zero.',
583 'evx': '{} eigenvectors failed to converge.',
584 'evd': 'The algorithm failed to compute an eigenvalue '
585 'while working on the submatrix lying in rows '
586 'and columns {0}/{1} through mod({0},{1}).',
587 'evr': 'Internal Error.'
588 }
589 if driver in ['ev', 'gv']:
590 msg = drv_err['ev'].format(info)
591 elif driver in ['evx', 'gvx']:
592 msg = drv_err['evx'].format(info)
593 elif driver in ['evd', 'gvd']:
594 if eigvals_only:
595 msg = drv_err['ev'].format(info)
596 else:
597 msg = drv_err['evd'].format(info, n+1)
598 else:
599 msg = drv_err['evr']
601 raise LinAlgError(msg)
604_conv_dict = {0: 0, 1: 1, 2: 2,
605 'all': 0, 'value': 1, 'index': 2,
606 'a': 0, 'v': 1, 'i': 2}
609def _check_select(select, select_range, max_ev, max_len):
610 """Check that select is valid, convert to Fortran style."""
611 if isinstance(select, str):
612 select = select.lower()
613 try:
614 select = _conv_dict[select]
615 except KeyError:
616 raise ValueError('invalid argument for select')
617 vl, vu = 0., 1.
618 il = iu = 1
619 if select != 0: # (non-all)
620 sr = asarray(select_range)
621 if sr.ndim != 1 or sr.size != 2 or sr[1] < sr[0]:
622 raise ValueError('select_range must be a 2-element array-like '
623 'in nondecreasing order')
624 if select == 1: # (value)
625 vl, vu = sr
626 if max_ev == 0:
627 max_ev = max_len
628 else: # 2 (index)
629 if sr.dtype.char.lower() not in 'hilqp':
630 raise ValueError('when using select="i", select_range must '
631 'contain integers, got dtype %s (%s)'
632 % (sr.dtype, sr.dtype.char))
633 # translate Python (0 ... N-1) into Fortran (1 ... N) with + 1
634 il, iu = sr + 1
635 if min(il, iu) < 1 or max(il, iu) > max_len:
636 raise ValueError('select_range out of bounds')
637 max_ev = iu - il + 1
638 return select, vl, vu, il, iu, max_ev
641def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
642 select='a', select_range=None, max_ev=0, check_finite=True):
643 """
644 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
646 Find eigenvalues w and optionally right eigenvectors v of a::
648 a v[:,i] = w[i] v[:,i]
649 v.H v = identity
651 The matrix a is stored in a_band either in lower diagonal or upper
652 diagonal ordered form:
654 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
655 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
657 where u is the number of bands above the diagonal.
659 Example of a_band (shape of a is (6,6), u=2)::
661 upper form:
662 * * a02 a13 a24 a35
663 * a01 a12 a23 a34 a45
664 a00 a11 a22 a33 a44 a55
666 lower form:
667 a00 a11 a22 a33 a44 a55
668 a10 a21 a32 a43 a54 *
669 a20 a31 a42 a53 * *
671 Cells marked with * are not used.
673 Parameters
674 ----------
675 a_band : (u+1, M) array_like
676 The bands of the M by M matrix a.
677 lower : bool, optional
678 Is the matrix in the lower form. (Default is upper form)
679 eigvals_only : bool, optional
680 Compute only the eigenvalues and no eigenvectors.
681 (Default: calculate also eigenvectors)
682 overwrite_a_band : bool, optional
683 Discard data in a_band (may enhance performance)
684 select : {'a', 'v', 'i'}, optional
685 Which eigenvalues to calculate
687 ====== ========================================
688 select calculated
689 ====== ========================================
690 'a' All eigenvalues
691 'v' Eigenvalues in the interval (min, max]
692 'i' Eigenvalues with indices min <= i <= max
693 ====== ========================================
694 select_range : (min, max), optional
695 Range of selected eigenvalues
696 max_ev : int, optional
697 For select=='v', maximum number of eigenvalues expected.
698 For other values of select, has no meaning.
700 In doubt, leave this parameter untouched.
702 check_finite : bool, optional
703 Whether to check that the input matrix contains only finite numbers.
704 Disabling may give a performance gain, but may result in problems
705 (crashes, non-termination) if the inputs do contain infinities or NaNs.
707 Returns
708 -------
709 w : (M,) ndarray
710 The eigenvalues, in ascending order, each repeated according to its
711 multiplicity.
712 v : (M, M) float or complex ndarray
713 The normalized eigenvector corresponding to the eigenvalue w[i] is
714 the column v[:,i].
716 Raises
717 ------
718 LinAlgError
719 If eigenvalue computation does not converge.
721 See Also
722 --------
723 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
724 eig : eigenvalues and right eigenvectors of general arrays.
725 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
726 eigh_tridiagonal : eigenvalues and right eiegenvectors for
727 symmetric/Hermitian tridiagonal matrices
729 Examples
730 --------
731 >>> from scipy.linalg import eig_banded
732 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
733 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
734 >>> w, v = eig_banded(Ab, lower=True)
735 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
736 True
737 >>> w = eig_banded(Ab, lower=True, eigvals_only=True)
738 >>> w
739 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
741 Request only the eigenvalues between ``[-3, 4]``
743 >>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4])
744 >>> w
745 array([-2.22987175, 3.95222349])
747 """
748 if eigvals_only or overwrite_a_band:
749 a1 = _asarray_validated(a_band, check_finite=check_finite)
750 overwrite_a_band = overwrite_a_band or (_datacopied(a1, a_band))
751 else:
752 a1 = array(a_band)
753 if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
754 raise ValueError("array must not contain infs or NaNs")
755 overwrite_a_band = 1
757 if len(a1.shape) != 2:
758 raise ValueError('expected a 2-D array')
759 select, vl, vu, il, iu, max_ev = _check_select(
760 select, select_range, max_ev, a1.shape[1])
761 del select_range
762 if select == 0:
763 if a1.dtype.char in 'GFD':
764 # FIXME: implement this somewhen, for now go with builtin values
765 # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
766 # or by using calc_lwork.f ???
767 # lwork = calc_lwork.hbevd(bevd.typecode, a1.shape[0], lower)
768 internal_name = 'hbevd'
769 else: # a1.dtype.char in 'fd':
770 # FIXME: implement this somewhen, for now go with builtin values
771 # see above
772 # lwork = calc_lwork.sbevd(bevd.typecode, a1.shape[0], lower)
773 internal_name = 'sbevd'
774 bevd, = get_lapack_funcs((internal_name,), (a1,))
775 w, v, info = bevd(a1, compute_v=not eigvals_only,
776 lower=lower, overwrite_ab=overwrite_a_band)
777 else: # select in [1, 2]
778 if eigvals_only:
779 max_ev = 1
780 # calculate optimal abstol for dsbevx (see manpage)
781 if a1.dtype.char in 'fF': # single precision
782 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),))
783 else:
784 lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),))
785 abstol = 2 * lamch('s')
786 if a1.dtype.char in 'GFD':
787 internal_name = 'hbevx'
788 else: # a1.dtype.char in 'gfd'
789 internal_name = 'sbevx'
790 bevx, = get_lapack_funcs((internal_name,), (a1,))
791 w, v, m, ifail, info = bevx(
792 a1, vl, vu, il, iu, compute_v=not eigvals_only, mmax=max_ev,
793 range=select, lower=lower, overwrite_ab=overwrite_a_band,
794 abstol=abstol)
795 # crop off w and v
796 w = w[:m]
797 if not eigvals_only:
798 v = v[:, :m]
799 _check_info(info, internal_name)
801 if eigvals_only:
802 return w
803 return w, v
806def eigvals(a, b=None, overwrite_a=False, check_finite=True,
807 homogeneous_eigvals=False):
808 """
809 Compute eigenvalues from an ordinary or generalized eigenvalue problem.
811 Find eigenvalues of a general matrix::
813 a vr[:,i] = w[i] b vr[:,i]
815 Parameters
816 ----------
817 a : (M, M) array_like
818 A complex or real matrix whose eigenvalues and eigenvectors
819 will be computed.
820 b : (M, M) array_like, optional
821 Right-hand side matrix in a generalized eigenvalue problem.
822 If omitted, identity matrix is assumed.
823 overwrite_a : bool, optional
824 Whether to overwrite data in a (may improve performance)
825 check_finite : bool, optional
826 Whether to check that the input matrices contain only finite numbers.
827 Disabling may give a performance gain, but may result in problems
828 (crashes, non-termination) if the inputs do contain infinities
829 or NaNs.
830 homogeneous_eigvals : bool, optional
831 If True, return the eigenvalues in homogeneous coordinates.
832 In this case ``w`` is a (2, M) array so that::
834 w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
836 Default is False.
838 Returns
839 -------
840 w : (M,) or (2, M) double or complex ndarray
841 The eigenvalues, each repeated according to its multiplicity
842 but not in any specific order. The shape is (M,) unless
843 ``homogeneous_eigvals=True``.
845 Raises
846 ------
847 LinAlgError
848 If eigenvalue computation does not converge
850 See Also
851 --------
852 eig : eigenvalues and right eigenvectors of general arrays.
853 eigvalsh : eigenvalues of symmetric or Hermitian arrays
854 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
855 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
856 matrices
858 Examples
859 --------
860 >>> from scipy import linalg
861 >>> a = np.array([[0., -1.], [1., 0.]])
862 >>> linalg.eigvals(a)
863 array([0.+1.j, 0.-1.j])
865 >>> b = np.array([[0., 1.], [1., 1.]])
866 >>> linalg.eigvals(a, b)
867 array([ 1.+0.j, -1.+0.j])
869 >>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
870 >>> linalg.eigvals(a, homogeneous_eigvals=True)
871 array([[3.+0.j, 8.+0.j, 7.+0.j],
872 [1.+0.j, 1.+0.j, 1.+0.j]])
874 """
875 return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a,
876 check_finite=check_finite,
877 homogeneous_eigvals=homogeneous_eigvals)
880def eigvalsh(a, b=None, lower=True, overwrite_a=False,
881 overwrite_b=False, turbo=True, eigvals=None, type=1,
882 check_finite=True, subset_by_index=None, subset_by_value=None,
883 driver=None):
884 """
885 Solves a standard or generalized eigenvalue problem for a complex
886 Hermitian or real symmetric matrix.
888 Find eigenvalues array ``w`` of array ``a``, where ``b`` is positive
889 definite such that for every eigenvalue λ (i-th entry of w) and its
890 eigenvector vi (i-th column of v) satisfies::
892 a @ vi = λ * b @ vi
893 vi.conj().T @ a @ vi = λ
894 vi.conj().T @ b @ vi = 1
896 In the standard problem, b is assumed to be the identity matrix.
898 Parameters
899 ----------
900 a : (M, M) array_like
901 A complex Hermitian or real symmetric matrix whose eigenvalues will
902 be computed.
903 b : (M, M) array_like, optional
904 A complex Hermitian or real symmetric definite positive matrix in.
905 If omitted, identity matrix is assumed.
906 lower : bool, optional
907 Whether the pertinent array data is taken from the lower or upper
908 triangle of `a`. (Default: lower)
909 eigvals_only : bool, optional
910 Whether to calculate only eigenvalues and no eigenvectors.
911 (Default: both are calculated)
912 subset_by_index : iterable, optional
913 If provided, this two-element iterable defines the start and the end
914 indices of the desired eigenvalues (ascending order and 0-indexed).
915 To return only the second smallest to fifth smallest eigenvalues,
916 ``[1, 4]`` is used. ``[n-3, n-1]`` returns the largest three. Only
917 available with "evr", "evx", and "gvx" drivers. The entries are
918 directly converted to integers via ``int()``.
919 subset_by_value : iterable, optional
920 If provided, this two-element iterable defines the half-open interval
921 ``(a, b]`` that, if any, only the eigenvalues between these values
922 are returned. Only available with "evr", "evx", and "gvx" drivers. Use
923 ``np.inf`` for the unconstrained ends.
924 driver: str, optional
925 Defines which LAPACK driver should be used. Valid options are "ev",
926 "evd", "evr", "evx" for standard problems and "gv", "gvd", "gvx" for
927 generalized (where b is not None) problems. See the Notes section of
928 `scipy.linalg.eigh`.
929 type : int, optional
930 For the generalized problems, this keyword specifies the problem type
931 to be solved for ``w`` and ``v`` (only takes 1, 2, 3 as possible
932 inputs)::
934 1 => a @ v = w @ b @ v
935 2 => a @ b @ v = w @ v
936 3 => b @ a @ v = w @ v
938 This keyword is ignored for standard problems.
939 overwrite_a : bool, optional
940 Whether to overwrite data in ``a`` (may improve performance). Default
941 is False.
942 overwrite_b : bool, optional
943 Whether to overwrite data in ``b`` (may improve performance). Default
944 is False.
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 turbo : bool, optional
950 *Deprecated by ``driver=gvd`` option*. Has no significant effect for
951 eigenvalue computations since no eigenvectors are requested.
953 ..Deprecated in v1.5.0
954 eigvals : tuple (lo, hi), optional
955 *Deprecated by ``subset_by_index`` keyword*. Indexes of the smallest
956 and largest (in ascending order) eigenvalues and corresponding
957 eigenvectors to be returned: 0 <= lo <= hi <= M-1. If omitted, all
958 eigenvalues and eigenvectors are returned.
960 .. Deprecated in v1.5.0
962 Returns
963 -------
964 w : (N,) ndarray
965 The ``N`` (``1<=N<=M``) selected eigenvalues, in ascending order, each
966 repeated according to its multiplicity.
968 Raises
969 ------
970 LinAlgError
971 If eigenvalue computation does not converge, an error occurred, or
972 b matrix is not definite positive. Note that if input matrices are
973 not symmetric or Hermitian, no error will be reported but results will
974 be wrong.
976 See Also
977 --------
978 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
979 eigvals : eigenvalues of general arrays
980 eigvals_banded : eigenvalues for symmetric/Hermitian band matrices
981 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
982 matrices
984 Notes
985 -----
986 This function does not check the input array for being Hermitian/symmetric
987 in order to allow for representing arrays with only their upper/lower
988 triangular parts.
990 This function serves as a one-liner shorthand for `scipy.linalg.eigh` with
991 the option ``eigvals_only=True`` to get the eigenvalues and not the
992 eigenvectors. Here it is kept as a legacy convenience. It might be
993 beneficial to use the main function to have full control and to be a bit
994 more pythonic.
996 Examples
997 --------
998 For more examples see `scipy.linalg.eigh`.
1000 >>> from scipy.linalg import eigvalsh
1001 >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
1002 >>> w = eigvalsh(A)
1003 >>> w
1004 array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
1006 """
1007 return eigh(a, b=b, lower=lower, eigvals_only=True,
1008 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
1009 turbo=turbo, eigvals=None, type=type,
1010 check_finite=check_finite, subset_by_index=eigvals,
1011 subset_by_value=None, driver=None)
1014def eigvals_banded(a_band, lower=False, overwrite_a_band=False,
1015 select='a', select_range=None, check_finite=True):
1016 """
1017 Solve real symmetric or complex Hermitian band matrix eigenvalue problem.
1019 Find eigenvalues w of a::
1021 a v[:,i] = w[i] v[:,i]
1022 v.H v = identity
1024 The matrix a is stored in a_band either in lower diagonal or upper
1025 diagonal ordered form:
1027 a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
1028 a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
1030 where u is the number of bands above the diagonal.
1032 Example of a_band (shape of a is (6,6), u=2)::
1034 upper form:
1035 * * a02 a13 a24 a35
1036 * a01 a12 a23 a34 a45
1037 a00 a11 a22 a33 a44 a55
1039 lower form:
1040 a00 a11 a22 a33 a44 a55
1041 a10 a21 a32 a43 a54 *
1042 a20 a31 a42 a53 * *
1044 Cells marked with * are not used.
1046 Parameters
1047 ----------
1048 a_band : (u+1, M) array_like
1049 The bands of the M by M matrix a.
1050 lower : bool, optional
1051 Is the matrix in the lower form. (Default is upper form)
1052 overwrite_a_band : bool, optional
1053 Discard data in a_band (may enhance performance)
1054 select : {'a', 'v', 'i'}, optional
1055 Which eigenvalues to calculate
1057 ====== ========================================
1058 select calculated
1059 ====== ========================================
1060 'a' All eigenvalues
1061 'v' Eigenvalues in the interval (min, max]
1062 'i' Eigenvalues with indices min <= i <= max
1063 ====== ========================================
1064 select_range : (min, max), optional
1065 Range of selected eigenvalues
1066 check_finite : bool, optional
1067 Whether to check that the input matrix contains only finite numbers.
1068 Disabling may give a performance gain, but may result in problems
1069 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1071 Returns
1072 -------
1073 w : (M,) ndarray
1074 The eigenvalues, in ascending order, each repeated according to its
1075 multiplicity.
1077 Raises
1078 ------
1079 LinAlgError
1080 If eigenvalue computation does not converge.
1082 See Also
1083 --------
1084 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1085 band matrices
1086 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1087 matrices
1088 eigvals : eigenvalues of general arrays
1089 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1090 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1092 Examples
1093 --------
1094 >>> from scipy.linalg import eigvals_banded
1095 >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
1096 >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
1097 >>> w = eigvals_banded(Ab, lower=True)
1098 >>> w
1099 array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
1100 """
1101 return eig_banded(a_band, lower=lower, eigvals_only=1,
1102 overwrite_a_band=overwrite_a_band, select=select,
1103 select_range=select_range, check_finite=check_finite)
1106def eigvalsh_tridiagonal(d, e, select='a', select_range=None,
1107 check_finite=True, tol=0., lapack_driver='auto'):
1108 """
1109 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1111 Find eigenvalues `w` of ``a``::
1113 a v[:,i] = w[i] v[:,i]
1114 v.H v = identity
1116 For a real symmetric matrix ``a`` with diagonal elements `d` and
1117 off-diagonal elements `e`.
1119 Parameters
1120 ----------
1121 d : ndarray, shape (ndim,)
1122 The diagonal elements of the array.
1123 e : ndarray, shape (ndim-1,)
1124 The off-diagonal elements of the array.
1125 select : {'a', 'v', 'i'}, optional
1126 Which eigenvalues to calculate
1128 ====== ========================================
1129 select calculated
1130 ====== ========================================
1131 'a' All eigenvalues
1132 'v' Eigenvalues in the interval (min, max]
1133 'i' Eigenvalues with indices min <= i <= max
1134 ====== ========================================
1135 select_range : (min, max), optional
1136 Range of selected eigenvalues
1137 check_finite : bool, optional
1138 Whether to check that the input matrix contains only finite numbers.
1139 Disabling may give a performance gain, but may result in problems
1140 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1141 tol : float
1142 The absolute tolerance to which each eigenvalue is required
1143 (only used when ``lapack_driver='stebz'``).
1144 An eigenvalue (or cluster) is considered to have converged if it
1145 lies in an interval of this width. If <= 0. (default),
1146 the value ``eps*|a|`` is used where eps is the machine precision,
1147 and ``|a|`` is the 1-norm of the matrix ``a``.
1148 lapack_driver : str
1149 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1150 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1151 and 'stebz' otherwise. 'sterf' and 'stev' can only be used when
1152 ``select='a'``.
1154 Returns
1155 -------
1156 w : (M,) ndarray
1157 The eigenvalues, in ascending order, each repeated according to its
1158 multiplicity.
1160 Raises
1161 ------
1162 LinAlgError
1163 If eigenvalue computation does not converge.
1165 See Also
1166 --------
1167 eigh_tridiagonal : eigenvalues and right eiegenvectors for
1168 symmetric/Hermitian tridiagonal matrices
1170 Examples
1171 --------
1172 >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
1173 >>> d = 3*np.ones(4)
1174 >>> e = -1*np.ones(3)
1175 >>> w = eigvalsh_tridiagonal(d, e)
1176 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1177 >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
1178 >>> np.allclose(w - w2, np.zeros(4))
1179 True
1180 """
1181 return eigh_tridiagonal(
1182 d, e, eigvals_only=True, select=select, select_range=select_range,
1183 check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
1186def eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None,
1187 check_finite=True, tol=0., lapack_driver='auto'):
1188 """
1189 Solve eigenvalue problem for a real symmetric tridiagonal matrix.
1191 Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
1193 a v[:,i] = w[i] v[:,i]
1194 v.H v = identity
1196 For a real symmetric matrix ``a`` with diagonal elements `d` and
1197 off-diagonal elements `e`.
1199 Parameters
1200 ----------
1201 d : ndarray, shape (ndim,)
1202 The diagonal elements of the array.
1203 e : ndarray, shape (ndim-1,)
1204 The off-diagonal elements of the array.
1205 select : {'a', 'v', 'i'}, optional
1206 Which eigenvalues to calculate
1208 ====== ========================================
1209 select calculated
1210 ====== ========================================
1211 'a' All eigenvalues
1212 'v' Eigenvalues in the interval (min, max]
1213 'i' Eigenvalues with indices min <= i <= max
1214 ====== ========================================
1215 select_range : (min, max), optional
1216 Range of selected eigenvalues
1217 check_finite : bool, optional
1218 Whether to check that the input matrix contains only finite numbers.
1219 Disabling may give a performance gain, but may result in problems
1220 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1221 tol : float
1222 The absolute tolerance to which each eigenvalue is required
1223 (only used when 'stebz' is the `lapack_driver`).
1224 An eigenvalue (or cluster) is considered to have converged if it
1225 lies in an interval of this width. If <= 0. (default),
1226 the value ``eps*|a|`` is used where eps is the machine precision,
1227 and ``|a|`` is the 1-norm of the matrix ``a``.
1228 lapack_driver : str
1229 LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf',
1230 or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'``
1231 and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and
1232 ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is
1233 used to find the corresponding eigenvectors. 'sterf' can only be
1234 used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only
1235 be used when ``select='a'``.
1237 Returns
1238 -------
1239 w : (M,) ndarray
1240 The eigenvalues, in ascending order, each repeated according to its
1241 multiplicity.
1242 v : (M, M) ndarray
1243 The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
1244 the column ``v[:,i]``.
1246 Raises
1247 ------
1248 LinAlgError
1249 If eigenvalue computation does not converge.
1251 See Also
1252 --------
1253 eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal
1254 matrices
1255 eig : eigenvalues and right eigenvectors for non-symmetric arrays
1256 eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
1257 eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian
1258 band matrices
1260 Notes
1261 -----
1262 This function makes use of LAPACK ``S/DSTEMR`` routines.
1264 Examples
1265 --------
1266 >>> from scipy.linalg import eigh_tridiagonal
1267 >>> d = 3*np.ones(4)
1268 >>> e = -1*np.ones(3)
1269 >>> w, v = eigh_tridiagonal(d, e)
1270 >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
1271 >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
1272 True
1273 """
1274 d = _asarray_validated(d, check_finite=check_finite)
1275 e = _asarray_validated(e, check_finite=check_finite)
1276 for check in (d, e):
1277 if check.ndim != 1:
1278 raise ValueError('expected a 1-D array')
1279 if check.dtype.char in 'GFD': # complex
1280 raise TypeError('Only real arrays currently supported')
1281 if d.size != e.size + 1:
1282 raise ValueError('d (%s) must have one more element than e (%s)'
1283 % (d.size, e.size))
1284 select, vl, vu, il, iu, _ = _check_select(
1285 select, select_range, 0, d.size)
1286 if not isinstance(lapack_driver, str):
1287 raise TypeError('lapack_driver must be str')
1288 drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev')
1289 if lapack_driver not in drivers:
1290 raise ValueError('lapack_driver must be one of %s, got %s'
1291 % (drivers, lapack_driver))
1292 if lapack_driver == 'auto':
1293 lapack_driver = 'stemr' if select == 0 else 'stebz'
1294 func, = get_lapack_funcs((lapack_driver,), (d, e))
1295 compute_v = not eigvals_only
1296 if lapack_driver == 'sterf':
1297 if select != 0:
1298 raise ValueError('sterf can only be used when select == "a"')
1299 if not eigvals_only:
1300 raise ValueError('sterf can only be used when eigvals_only is '
1301 'True')
1302 w, info = func(d, e)
1303 m = len(w)
1304 elif lapack_driver == 'stev':
1305 if select != 0:
1306 raise ValueError('stev can only be used when select == "a"')
1307 w, v, info = func(d, e, compute_v=compute_v)
1308 m = len(w)
1309 elif lapack_driver == 'stebz':
1310 tol = float(tol)
1311 internal_name = 'stebz'
1312 stebz, = get_lapack_funcs((internal_name,), (d, e))
1313 # If getting eigenvectors, needs to be block-ordered (B) instead of
1314 # matrix-ordered (E), and we will reorder later
1315 order = 'E' if eigvals_only else 'B'
1316 m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol,
1317 order)
1318 else: # 'stemr'
1319 # ?STEMR annoyingly requires size N instead of N-1
1320 e_ = empty(e.size+1, e.dtype)
1321 e_[:-1] = e
1322 stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e))
1323 lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu,
1324 compute_v=compute_v)
1325 _check_info(info, 'stemr_lwork')
1326 m, w, v, info = func(d, e_, select, vl, vu, il, iu,
1327 compute_v=compute_v, lwork=lwork, liwork=liwork)
1328 _check_info(info, lapack_driver + ' (eigh_tridiagonal)')
1329 w = w[:m]
1330 if eigvals_only:
1331 return w
1332 else:
1333 # Do we still need to compute the eigenvalues?
1334 if lapack_driver == 'stebz':
1335 func, = get_lapack_funcs(('stein',), (d, e))
1336 v, info = func(d, e, w, iblock, isplit)
1337 _check_info(info, 'stein (eigh_tridiagonal)',
1338 positive='%d eigenvectors failed to converge')
1339 # Convert block-order to matrix-order
1340 order = argsort(w)
1341 w, v = w[order], v[:, order]
1342 else:
1343 v = v[:, :m]
1344 return w, v
1347def _check_info(info, driver, positive='did not converge (LAPACK info=%d)'):
1348 """Check info return value."""
1349 if info < 0:
1350 raise ValueError('illegal value in argument %d of internal %s'
1351 % (-info, driver))
1352 if info > 0 and positive:
1353 raise LinAlgError(("%s " + positive) % (driver, info,))
1356def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True):
1357 """
1358 Compute Hessenberg form of a matrix.
1360 The Hessenberg decomposition is::
1362 A = Q H Q^H
1364 where `Q` is unitary/orthogonal and `H` has only zero elements below
1365 the first sub-diagonal.
1367 Parameters
1368 ----------
1369 a : (M, M) array_like
1370 Matrix to bring into Hessenberg form.
1371 calc_q : bool, optional
1372 Whether to compute the transformation matrix. Default is False.
1373 overwrite_a : bool, optional
1374 Whether to overwrite `a`; may improve performance.
1375 Default is False.
1376 check_finite : bool, optional
1377 Whether to check that the input matrix contains only finite numbers.
1378 Disabling may give a performance gain, but may result in problems
1379 (crashes, non-termination) if the inputs do contain infinities or NaNs.
1381 Returns
1382 -------
1383 H : (M, M) ndarray
1384 Hessenberg form of `a`.
1385 Q : (M, M) ndarray
1386 Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
1387 Only returned if ``calc_q=True``.
1389 Examples
1390 --------
1391 >>> from scipy.linalg import hessenberg
1392 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
1393 >>> H, Q = hessenberg(A, calc_q=True)
1394 >>> H
1395 array([[ 2. , -11.65843866, 1.42005301, 0.25349066],
1396 [ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
1397 [ 0. , -1.83299243, 0.38969961, -0.51527034],
1398 [ 0. , 0. , -3.83189513, 1.07494686]])
1399 >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
1400 True
1401 """
1402 a1 = _asarray_validated(a, check_finite=check_finite)
1403 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
1404 raise ValueError('expected square matrix')
1405 overwrite_a = overwrite_a or (_datacopied(a1, a))
1407 # if 2x2 or smaller: already in Hessenberg
1408 if a1.shape[0] <= 2:
1409 if calc_q:
1410 return a1, eye(a1.shape[0])
1411 return a1
1413 gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal',
1414 'gehrd_lwork'), (a1,))
1415 ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a)
1416 _check_info(info, 'gebal (hessenberg)', positive=False)
1417 n = len(a1)
1419 lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
1421 hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1422 _check_info(info, 'gehrd (hessenberg)', positive=False)
1423 h = numpy.triu(hq, -1)
1424 if not calc_q:
1425 return h
1427 # use orghr/unghr to compute q
1428 orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,))
1429 lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
1431 q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
1432 _check_info(info, 'orghr (hessenberg)', positive=False)
1433 return h, q
1436def cdf2rdf(w, v):
1437 """
1438 Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real
1439 eigenvalues in a block diagonal form ``wr`` and the associated real
1440 eigenvectors ``vr``, such that::
1442 vr @ wr = X @ vr
1444 continues to hold, where ``X`` is the original array for which ``w`` and
1445 ``v`` are the eigenvalues and eigenvectors.
1447 .. versionadded:: 1.1.0
1449 Parameters
1450 ----------
1451 w : (..., M) array_like
1452 Complex or real eigenvalues, an array or stack of arrays
1454 Conjugate pairs must not be interleaved, else the wrong result
1455 will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result,
1456 but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
1458 v : (..., M, M) array_like
1459 Complex or real eigenvectors, a square array or stack of square arrays.
1461 Returns
1462 -------
1463 wr : (..., M, M) ndarray
1464 Real diagonal block form of eigenvalues
1465 vr : (..., M, M) ndarray
1466 Real eigenvectors associated with ``wr``
1468 See Also
1469 --------
1470 eig : Eigenvalues and right eigenvectors for non-symmetric arrays
1471 rsf2csf : Convert real Schur form to complex Schur form
1473 Notes
1474 -----
1475 ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``.
1476 For example, obtained by ``w, v = scipy.linalg.eig(X)`` or
1477 ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent
1478 stacked arrays.
1480 .. versionadded:: 1.1.0
1482 Examples
1483 --------
1484 >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
1485 >>> X
1486 array([[ 1, 2, 3],
1487 [ 0, 4, 5],
1488 [ 0, -5, 4]])
1490 >>> from scipy import linalg
1491 >>> w, v = linalg.eig(X)
1492 >>> w
1493 array([ 1.+0.j, 4.+5.j, 4.-5.j])
1494 >>> v
1495 array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
1496 [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
1497 [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
1499 >>> wr, vr = linalg.cdf2rdf(w, v)
1500 >>> wr
1501 array([[ 1., 0., 0.],
1502 [ 0., 4., 5.],
1503 [ 0., -5., 4.]])
1504 >>> vr
1505 array([[ 1. , 0.40016, -0.01906],
1506 [ 0. , 0.64788, 0. ],
1507 [ 0. , 0. , 0.64788]])
1509 >>> vr @ wr
1510 array([[ 1. , 1.69593, 1.9246 ],
1511 [ 0. , 2.59153, 3.23942],
1512 [ 0. , -3.23942, 2.59153]])
1513 >>> X @ vr
1514 array([[ 1. , 1.69593, 1.9246 ],
1515 [ 0. , 2.59153, 3.23942],
1516 [ 0. , -3.23942, 2.59153]])
1517 """
1518 w, v = _asarray_validated(w), _asarray_validated(v)
1520 # check dimensions
1521 if w.ndim < 1:
1522 raise ValueError('expected w to be at least 1D')
1523 if v.ndim < 2:
1524 raise ValueError('expected v to be at least 2D')
1525 if v.ndim != w.ndim + 1:
1526 raise ValueError('expected eigenvectors array to have exactly one '
1527 'dimension more than eigenvalues array')
1529 # check shapes
1530 n = w.shape[-1]
1531 M = w.shape[:-1]
1532 if v.shape[-2] != v.shape[-1]:
1533 raise ValueError('expected v to be a square matrix or stacked square '
1534 'matrices: v.shape[-2] = v.shape[-1]')
1535 if v.shape[-1] != n:
1536 raise ValueError('expected the same number of eigenvalues as '
1537 'eigenvectors')
1539 # get indices for each first pair of complex eigenvalues
1540 complex_mask = iscomplex(w)
1541 n_complex = complex_mask.sum(axis=-1)
1543 # check if all complex eigenvalues have conjugate pairs
1544 if not (n_complex % 2 == 0).all():
1545 raise ValueError('expected complex-conjugate pairs of eigenvalues')
1547 # find complex indices
1548 idx = nonzero(complex_mask)
1549 idx_stack = idx[:-1]
1550 idx_elem = idx[-1]
1552 # filter them to conjugate indices, assuming pairs are not interleaved
1553 j = idx_elem[0::2]
1554 k = idx_elem[1::2]
1555 stack_ind = ()
1556 for i in idx_stack:
1557 # should never happen, assuming nonzero orders by the last axis
1558 assert (i[0::2] == i[1::2]).all(),\
1559 "Conjugate pair spanned different arrays!"
1560 stack_ind += (i[0::2],)
1562 # all eigenvalues to diagonal form
1563 wr = zeros(M + (n, n), dtype=w.real.dtype)
1564 di = range(n)
1565 wr[..., di, di] = w.real
1567 # complex eigenvalues to real block diagonal form
1568 wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag
1569 wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
1571 # compute real eigenvectors associated with real block diagonal eigenvalues
1572 u = zeros(M + (n, n), dtype=numpy.cdouble)
1573 u[..., di, di] = 1.0
1574 u[stack_ind + (j, j)] = 0.5j
1575 u[stack_ind + (j, k)] = 0.5
1576 u[stack_ind + (k, j)] = -0.5j
1577 u[stack_ind + (k, k)] = 0.5
1579 # multipy matrices v and u (equivalent to v @ u)
1580 vr = einsum('...ij,...jk->...ik', v, u).real
1582 return wr, vr