|
# # Author: Pearu Peterson, March 2002 # # additions by Travis Oliphant, March 2002 # additions by Eric Jones, June 2002 # additions by Johannes Loehnert, June 2006 # additions by Bart Vandereycken, June 2006 # additions by Andrew D Straw, May 2007 # additions by Tiziano Zito, November 2008 # # April 2010: Functions for LU, QR, SVD, Schur and Cholesky decompositions were # moved to their own files. Still in this file are functions for eigenstuff # and for the Hessenberg form.
'eig_banded', 'eigvals_banded', 'eigh_tridiagonal', 'eigvalsh_tridiagonal', 'hessenberg', 'cdf2rdf']
flatnonzero, conj, asarray, argsort, empty, newaxis, argwhere, iscomplex, eye, zeros, einsum) # Local imports
""" Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output """ # - see LAPACK man page DGGEV at ALPHAI v = numpy.array(vin, dtype=dtype) m = (w.imag > 0) m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709 for i in flatnonzero(m): v.imag[:, i] = vin[:, i+1] conj(v[:, i], v[:, i+1]) return v
if homogeneous_eigvals: if beta is None: return numpy.vstack((alpha, numpy.ones_like(alpha))) else: return numpy.vstack((alpha, beta)) else: if beta is None: return alpha else: w = numpy.empty_like(alpha) alpha_zero = (alpha == 0) beta_zero = (beta == 0) beta_nonzero = ~beta_zero w[beta_nonzero] = alpha[beta_nonzero]/beta[beta_nonzero] # Use numpy.inf for complex values too since # 1/numpy.inf = 0, i.e. it correctly behaves as projective # infinity. w[~alpha_zero & beta_zero] = numpy.inf if numpy.all(alpha.imag == 0): w[alpha_zero & beta_zero] = numpy.nan else: w[alpha_zero & beta_zero] = complex(numpy.nan, numpy.nan) return w
homogeneous_eigvals): ggev, = get_lapack_funcs(('ggev',), (a1, b1)) cvl, cvr = left, right res = ggev(a1, b1, lwork=-1) lwork = res[-2][0].real.astype(numpy.int) if ggev.typecode in 'cz': alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork, overwrite_a, overwrite_b) w = _make_eigvals(alpha, beta, homogeneous_eigvals) else: alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork, overwrite_a, overwrite_b) alpha = alphar + _I * alphai w = _make_eigvals(alpha, beta, homogeneous_eigvals) _check_info(info, 'generalized eig algorithm (ggev)')
only_real = numpy.all(w.imag == 0.0) if not (ggev.typecode in 'cz' or only_real): t = w.dtype.char if left: vl = _make_complex_eigvecs(w, vl, t) if right: vr = _make_complex_eigvecs(w, vr, t)
# the eigenvectors returned by the lapack function are NOT normalized for i in xrange(vr.shape[0]): if right: vr[:, i] /= norm(vr[:, i]) if left: vl[:, i] /= norm(vl[:, i])
if not (left or right): return w if left: if right: return w, vl, vr return w, vl return w, vr
overwrite_b=False, check_finite=True, homogeneous_eigvals=False): """ Solve an ordinary or generalized eigenvalue problem of a square matrix.
Find eigenvalues w and right or left eigenvectors of a general matrix::
a vr[:,i] = w[i] b vr[:,i] a.H vl[:,i] = w[i].conj() b.H vl[:,i]
where ``.H`` is the Hermitian conjugation.
Parameters ---------- a : (M, M) array_like A complex or real matrix whose eigenvalues and eigenvectors will be computed. b : (M, M) array_like, optional Right-hand side matrix in a generalized eigenvalue problem. Default is None, identity matrix is assumed. left : bool, optional Whether to calculate and return left eigenvectors. Default is False. right : bool, optional Whether to calculate and return right eigenvectors. Default is True. overwrite_a : bool, optional Whether to overwrite `a`; may improve performance. Default is False. overwrite_b : bool, optional Whether to overwrite `b`; may improve performance. Default is False. check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. homogeneous_eigvals : bool, optional If True, return the eigenvalues in homogeneous coordinates. In this case ``w`` is a (2, M) array so that::
w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
Default is False.
Returns ------- w : (M,) or (2, M) double or complex ndarray The eigenvalues, each repeated according to its multiplicity. The shape is (M,) unless ``homogeneous_eigvals=True``. vl : (M, M) double or complex ndarray The normalized left eigenvector corresponding to the eigenvalue ``w[i]`` is the column vl[:,i]. Only returned if ``left=True``. vr : (M, M) double or complex ndarray The normalized right eigenvector corresponding to the eigenvalue ``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
Raises ------ LinAlgError If eigenvalue computation does not converge.
See Also -------- eigvals : eigenvalues of general arrays eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays. eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices eigh_tridiagonal : eigenvalues and right eiegenvectors for symmetric/Hermitian tridiagonal matrices
Examples -------- >>> from scipy import linalg >>> a = np.array([[0., -1.], [1., 0.]]) >>> linalg.eigvals(a) array([0.+1.j, 0.-1.j])
>>> b = np.array([[0., 1.], [1., 1.]]) >>> linalg.eigvals(a, b) array([ 1.+0.j, -1.+0.j])
>>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]]) >>> linalg.eigvals(a, homogeneous_eigvals=True) array([[3.+0.j, 8.+0.j, 7.+0.j], [1.+0.j, 1.+0.j, 1.+0.j]])
>>> a = np.array([[0., -1.], [1., 0.]]) >>> linalg.eigvals(a) == linalg.eig(a)[0] array([ True, True]) >>> linalg.eig(a, left=True, right=False)[1] # normalized left eigenvector array([[-0.70710678+0.j , -0.70710678-0.j ], [-0. +0.70710678j, -0. -0.70710678j]]) >>> linalg.eig(a, left=False, right=True)[1] # normalized right eigenvector array([[0.70710678+0.j , 0.70710678-0.j ], [0. -0.70710678j, 0. +0.70710678j]])
""" a1 = _asarray_validated(a, check_finite=check_finite) if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: raise ValueError('expected square matrix') overwrite_a = overwrite_a or (_datacopied(a1, a)) if b is not None: b1 = _asarray_validated(b, check_finite=check_finite) overwrite_b = overwrite_b or _datacopied(b1, b) if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]: raise ValueError('expected square matrix') if b1.shape != a1.shape: raise ValueError('a and b must have the same shape') return _geneig(a1, b1, left, right, overwrite_a, overwrite_b, homogeneous_eigvals)
geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,)) compute_vl, compute_vr = left, right
lwork = _compute_lwork(geev_lwork, a1.shape[0], compute_vl=compute_vl, compute_vr=compute_vr)
if geev.typecode in 'cz': w, vl, vr, info = geev(a1, lwork=lwork, compute_vl=compute_vl, compute_vr=compute_vr, overwrite_a=overwrite_a) w = _make_eigvals(w, None, homogeneous_eigvals) else: wr, wi, vl, vr, info = geev(a1, lwork=lwork, compute_vl=compute_vl, compute_vr=compute_vr, overwrite_a=overwrite_a) t = {'f': 'F', 'd': 'D'}[wr.dtype.char] w = wr + _I * wi w = _make_eigvals(w, None, homogeneous_eigvals)
_check_info(info, 'eig algorithm (geev)', positive='did not converge (only eigenvalues ' 'with order >= %d have converged)')
only_real = numpy.all(w.imag == 0.0) if not (geev.typecode in 'cz' or only_real): t = w.dtype.char if left: vl = _make_complex_eigvecs(w, vl, t) if right: vr = _make_complex_eigvecs(w, vr, t) if not (left or right): return w if left: if right: return w, vl, vr return w, vl return w, vr
overwrite_b=False, turbo=True, eigvals=None, type=1, check_finite=True): """ Solve an ordinary or generalized eigenvalue problem for a complex Hermitian or real symmetric matrix.
Find eigenvalues w and optionally eigenvectors v of matrix `a`, where `b` is positive definite::
a v[:,i] = w[i] b v[:,i] v[i,:].conj() a v[:,i] = w[i] v[i,:].conj() b v[:,i] = 1
Parameters ---------- a : (M, M) array_like A complex Hermitian or real symmetric matrix whose eigenvalues and eigenvectors will be computed. b : (M, M) array_like, optional A complex Hermitian or real symmetric definite positive matrix in. If omitted, identity matrix is assumed. lower : bool, optional Whether the pertinent array data is taken from the lower or upper triangle of `a`. (Default: lower) eigvals_only : bool, optional Whether to calculate only eigenvalues and no eigenvectors. (Default: both are calculated) turbo : bool, optional Use divide and conquer algorithm (faster but expensive in memory, only for generalized eigenvalue problem and if eigvals=None) eigvals : tuple (lo, hi), optional Indexes of the smallest and largest (in ascending order) eigenvalues and corresponding eigenvectors to be returned: 0 <= lo <= hi <= M-1. If omitted, all eigenvalues and eigenvectors are returned. type : int, optional Specifies the problem type to be solved:
type = 1: a v[:,i] = w[i] b v[:,i]
type = 2: a b v[:,i] = w[i] v[:,i]
type = 3: b a v[:,i] = w[i] v[:,i] overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance) overwrite_b : bool, optional Whether to overwrite data in `b` (may improve performance) check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- w : (N,) float ndarray The N (1<=N<=M) selected eigenvalues, in ascending order, each repeated according to its multiplicity. v : (M, N) complex ndarray (if eigvals_only == False)
The normalized selected eigenvector corresponding to the eigenvalue w[i] is the column v[:,i].
Normalization:
type 1 and 3: v.conj() a v = w
type 2: inv(v).conj() a inv(v) = w
type = 1 or 2: v.conj() b v = I
type = 3: v.conj() inv(b) v = I
Raises ------ LinAlgError If eigenvalue computation does not converge, an error occurred, or b matrix is not definite positive. Note that if input matrices are not symmetric or hermitian, no error is reported but results will be wrong.
See Also -------- eigvalsh : eigenvalues of symmetric or Hermitian arrays eig : eigenvalues and right eigenvectors for non-symmetric arrays eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays eigh_tridiagonal : eigenvalues and right eiegenvectors for symmetric/Hermitian tridiagonal matrices
Notes ----- This function does not check the input array for being hermitian/symmetric in order to allow for representing arrays with only their upper/lower triangular parts.
Examples -------- >>> from scipy.linalg import eigh >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]]) >>> w, v = eigh(A) >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4))) True
""" a1 = _asarray_validated(a, check_finite=check_finite) if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: raise ValueError('expected square matrix') overwrite_a = overwrite_a or (_datacopied(a1, a)) if iscomplexobj(a1): cplx = True else: cplx = False if b is not None: b1 = _asarray_validated(b, check_finite=check_finite) overwrite_b = overwrite_b or _datacopied(b1, b) if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]: raise ValueError('expected square matrix')
if b1.shape != a1.shape: raise ValueError("wrong b dimensions %s, should " "be %s" % (str(b1.shape), str(a1.shape))) if iscomplexobj(b1): cplx = True else: cplx = cplx or False else: b1 = None
# Set job for fortran routines _job = (eigvals_only and 'N') or 'V'
# port eigenvalue range from python to fortran convention if eigvals is not None: lo, hi = eigvals if lo < 0 or hi >= a1.shape[0]: raise ValueError('The eigenvalue range specified is not valid.\n' 'Valid range is [%s,%s]' % (0, a1.shape[0]-1)) lo += 1 hi += 1 eigvals = (lo, hi)
# set lower if lower: uplo = 'L' else: uplo = 'U'
# fix prefix for lapack routines if cplx: pfx = 'he' else: pfx = 'sy'
# Standard Eigenvalue Problem # Use '*evr' routines # FIXME: implement calculation of optimal lwork # for all lapack routines if b1 is None: driver = pfx+'evr' (evr,) = get_lapack_funcs((driver,), (a1,)) if eigvals is None: w, v, info = evr(a1, uplo=uplo, jobz=_job, range="A", il=1, iu=a1.shape[0], overwrite_a=overwrite_a) else: (lo, hi) = eigvals w_tot, v, info = evr(a1, uplo=uplo, jobz=_job, range="I", il=lo, iu=hi, overwrite_a=overwrite_a) w = w_tot[0:hi-lo+1]
# Generalized Eigenvalue Problem else: # Use '*gvx' routines if range is specified if eigvals is not None: driver = pfx+'gvx' (gvx,) = get_lapack_funcs((driver,), (a1, b1)) (lo, hi) = eigvals w_tot, v, ifail, info = gvx(a1, b1, uplo=uplo, iu=hi, itype=type, jobz=_job, il=lo, overwrite_a=overwrite_a, overwrite_b=overwrite_b) w = w_tot[0:hi-lo+1] # Use '*gvd' routine if turbo is on and no eigvals are specified elif turbo: driver = pfx+'gvd' (gvd,) = get_lapack_funcs((driver,), (a1, b1)) v, w, info = gvd(a1, b1, uplo=uplo, itype=type, jobz=_job, overwrite_a=overwrite_a, overwrite_b=overwrite_b) # Use '*gv' routine if turbo is off and no eigvals are specified else: driver = pfx+'gv' (gv,) = get_lapack_funcs((driver,), (a1, b1)) v, w, info = gv(a1, b1, uplo=uplo, itype=type, jobz=_job, overwrite_a=overwrite_a, overwrite_b=overwrite_b)
# Check if we had a successful exit if info == 0: if eigvals_only: return w else: return w, v _check_info(info, driver, positive=False) # triage more specifically if info > 0 and b1 is None: raise LinAlgError("unrecoverable internal error.")
# The algorithm failed to converge. elif 0 < info <= b1.shape[0]: if eigvals is not None: raise LinAlgError("the eigenvectors %s failed to" " converge." % nonzero(ifail)-1) else: raise LinAlgError("internal fortran routine failed to converge: " "%i off-diagonal elements of an " "intermediate tridiagonal form did not converge" " to zero." % info)
# This occurs when b is not positive definite else: raise LinAlgError("the leading minor of order %i" " of 'b' is not positive definite. The" " factorization of 'b' could not be completed" " and no eigenvalues or eigenvectors were" " computed." % (info-b1.shape[0]))
'all': 0, 'value': 1, 'index': 2, 'a': 0, 'v': 1, 'i': 2}
"""Check that select is valid, convert to Fortran style.""" if isinstance(select, string_types): select = select.lower() try: select = _conv_dict[select] except KeyError: raise ValueError('invalid argument for select') vl, vu = 0., 1. il = iu = 1 if select != 0: # (non-all) sr = asarray(select_range) if sr.ndim != 1 or sr.size != 2 or sr[1] < sr[0]: raise ValueError('select_range must be a 2-element array-like ' 'in nondecreasing order') if select == 1: # (value) vl, vu = sr if max_ev == 0: max_ev = max_len else: # 2 (index) if sr.dtype.char.lower() not in 'hilqp': raise ValueError('when using select="i", select_range must ' 'contain integers, got dtype %s (%s)' % (sr.dtype, sr.dtype.char)) # translate Python (0 ... N-1) into Fortran (1 ... N) with + 1 il, iu = sr + 1 if min(il, iu) < 1 or max(il, iu) > max_len: raise ValueError('select_range out of bounds') max_ev = iu - il + 1 return select, vl, vu, il, iu, max_ev
select='a', select_range=None, max_ev=0, check_finite=True): """ Solve real symmetric or complex hermitian band matrix eigenvalue problem.
Find eigenvalues w and optionally right eigenvectors v of a::
a v[:,i] = w[i] v[:,i] v.H v = identity
The matrix a is stored in a_band either in lower diagonal or upper diagonal ordered form:
a_band[u + i - j, j] == a[i,j] (if upper form; i <= j) a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
where u is the number of bands above the diagonal.
Example of a_band (shape of a is (6,6), u=2)::
upper form: * * a02 a13 a24 a35 * a01 a12 a23 a34 a45 a00 a11 a22 a33 a44 a55
lower form: a00 a11 a22 a33 a44 a55 a10 a21 a32 a43 a54 * a20 a31 a42 a53 * *
Cells marked with * are not used.
Parameters ---------- a_band : (u+1, M) array_like The bands of the M by M matrix a. lower : bool, optional Is the matrix in the lower form. (Default is upper form) eigvals_only : bool, optional Compute only the eigenvalues and no eigenvectors. (Default: calculate also eigenvectors) overwrite_a_band : bool, optional Discard data in a_band (may enhance performance) select : {'a', 'v', 'i'}, optional Which eigenvalues to calculate
====== ======================================== select calculated ====== ======================================== 'a' All eigenvalues 'v' Eigenvalues in the interval (min, max] 'i' Eigenvalues with indices min <= i <= max ====== ======================================== select_range : (min, max), optional Range of selected eigenvalues max_ev : int, optional For select=='v', maximum number of eigenvalues expected. For other values of select, has no meaning.
In doubt, leave this parameter untouched.
check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- w : (M,) ndarray The eigenvalues, in ascending order, each repeated according to its multiplicity. v : (M, M) float or complex ndarray The normalized eigenvector corresponding to the eigenvalue w[i] is the column v[:,i].
Raises ------ LinAlgError If eigenvalue computation does not converge.
See Also -------- eigvals_banded : eigenvalues for symmetric/Hermitian band matrices eig : eigenvalues and right eigenvectors of general arrays. eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays eigh_tridiagonal : eigenvalues and right eiegenvectors for symmetric/Hermitian tridiagonal matrices
Examples -------- >>> from scipy.linalg import eig_banded >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]]) >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]]) >>> w, v = eig_banded(Ab, lower=True) >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4))) True >>> w = eig_banded(Ab, lower=True, eigvals_only=True) >>> w array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
Request only the eigenvalues between ``[-3, 4]``
>>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4]) >>> w array([-2.22987175, 3.95222349])
""" if eigvals_only or overwrite_a_band: a1 = _asarray_validated(a_band, check_finite=check_finite) overwrite_a_band = overwrite_a_band or (_datacopied(a1, a_band)) else: a1 = array(a_band) if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all(): raise ValueError("array must not contain infs or NaNs") overwrite_a_band = 1
if len(a1.shape) != 2: raise ValueError('expected two-dimensional array') select, vl, vu, il, iu, max_ev = _check_select( select, select_range, max_ev, a1.shape[1]) del select_range if select == 0: if a1.dtype.char in 'GFD': # FIXME: implement this somewhen, for now go with builtin values # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1) # or by using calc_lwork.f ??? # lwork = calc_lwork.hbevd(bevd.typecode, a1.shape[0], lower) internal_name = 'hbevd' else: # a1.dtype.char in 'fd': # FIXME: implement this somewhen, for now go with builtin values # see above # lwork = calc_lwork.sbevd(bevd.typecode, a1.shape[0], lower) internal_name = 'sbevd' bevd, = get_lapack_funcs((internal_name,), (a1,)) w, v, info = bevd(a1, compute_v=not eigvals_only, lower=lower, overwrite_ab=overwrite_a_band) else: # select in [1, 2] if eigvals_only: max_ev = 1 # calculate optimal abstol for dsbevx (see manpage) if a1.dtype.char in 'fF': # single precision lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),)) else: lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),)) abstol = 2 * lamch('s') if a1.dtype.char in 'GFD': internal_name = 'hbevx' else: # a1.dtype.char in 'gfd' internal_name = 'sbevx' bevx, = get_lapack_funcs((internal_name,), (a1,)) w, v, m, ifail, info = bevx( a1, vl, vu, il, iu, compute_v=not eigvals_only, mmax=max_ev, range=select, lower=lower, overwrite_ab=overwrite_a_band, abstol=abstol) # crop off w and v w = w[:m] if not eigvals_only: v = v[:, :m] _check_info(info, internal_name)
if eigvals_only: return w return w, v
homogeneous_eigvals=False): """ Compute eigenvalues from an ordinary or generalized eigenvalue problem.
Find eigenvalues of a general matrix::
a vr[:,i] = w[i] b vr[:,i]
Parameters ---------- a : (M, M) array_like A complex or real matrix whose eigenvalues and eigenvectors will be computed. b : (M, M) array_like, optional Right-hand side matrix in a generalized eigenvalue problem. If omitted, identity matrix is assumed. overwrite_a : bool, optional Whether to overwrite data in a (may improve performance) check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. homogeneous_eigvals : bool, optional If True, return the eigenvalues in homogeneous coordinates. In this case ``w`` is a (2, M) array so that::
w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
Default is False.
Returns ------- w : (M,) or (2, M) double or complex ndarray The eigenvalues, each repeated according to its multiplicity but not in any specific order. The shape is (M,) unless ``homogeneous_eigvals=True``.
Raises ------ LinAlgError If eigenvalue computation does not converge
See Also -------- eig : eigenvalues and right eigenvectors of general arrays. eigvalsh : eigenvalues of symmetric or Hermitian arrays eigvals_banded : eigenvalues for symmetric/Hermitian band matrices eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal matrices
Examples -------- >>> from scipy import linalg >>> a = np.array([[0., -1.], [1., 0.]]) >>> linalg.eigvals(a) array([0.+1.j, 0.-1.j])
>>> b = np.array([[0., 1.], [1., 1.]]) >>> linalg.eigvals(a, b) array([ 1.+0.j, -1.+0.j])
>>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]]) >>> linalg.eigvals(a, homogeneous_eigvals=True) array([[3.+0.j, 8.+0.j, 7.+0.j], [1.+0.j, 1.+0.j, 1.+0.j]])
""" return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a, check_finite=check_finite, homogeneous_eigvals=homogeneous_eigvals)
overwrite_b=False, turbo=True, eigvals=None, type=1, check_finite=True): """ Solve an ordinary or generalized eigenvalue problem for a complex Hermitian or real symmetric matrix.
Find eigenvalues w of matrix a, where b is positive definite::
a v[:,i] = w[i] b v[:,i] v[i,:].conj() a v[:,i] = w[i] v[i,:].conj() b v[:,i] = 1
Parameters ---------- a : (M, M) array_like A complex Hermitian or real symmetric matrix whose eigenvalues and eigenvectors will be computed. b : (M, M) array_like, optional A complex Hermitian or real symmetric definite positive matrix in. If omitted, identity matrix is assumed. lower : bool, optional Whether the pertinent array data is taken from the lower or upper triangle of `a`. (Default: lower) turbo : bool, optional Use divide and conquer algorithm (faster but expensive in memory, only for generalized eigenvalue problem and if eigvals=None) eigvals : tuple (lo, hi), optional Indexes of the smallest and largest (in ascending order) eigenvalues and corresponding eigenvectors to be returned: 0 <= lo < hi <= M-1. If omitted, all eigenvalues and eigenvectors are returned. type : int, optional Specifies the problem type to be solved:
type = 1: a v[:,i] = w[i] b v[:,i]
type = 2: a b v[:,i] = w[i] v[:,i]
type = 3: b a v[:,i] = w[i] v[:,i] overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance) overwrite_b : bool, optional Whether to overwrite data in `b` (may improve performance) check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- w : (N,) float ndarray The N (1<=N<=M) selected eigenvalues, in ascending order, each repeated according to its multiplicity.
Raises ------ LinAlgError If eigenvalue computation does not converge, an error occurred, or b matrix is not definite positive. Note that if input matrices are not symmetric or hermitian, no error is reported but results will be wrong.
See Also -------- eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays eigvals : eigenvalues of general arrays eigvals_banded : eigenvalues for symmetric/Hermitian band matrices eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal matrices
Notes ----- This function does not check the input array for being hermitian/symmetric in order to allow for representing arrays with only their upper/lower triangular parts.
Examples -------- >>> from scipy.linalg import eigvalsh >>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]]) >>> w = eigvalsh(A) >>> w array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
""" return eigh(a, b=b, lower=lower, eigvals_only=True, overwrite_a=overwrite_a, overwrite_b=overwrite_b, turbo=turbo, eigvals=eigvals, type=type, check_finite=check_finite)
select='a', select_range=None, check_finite=True): """ Solve real symmetric or complex hermitian band matrix eigenvalue problem.
Find eigenvalues w of a::
a v[:,i] = w[i] v[:,i] v.H v = identity
The matrix a is stored in a_band either in lower diagonal or upper diagonal ordered form:
a_band[u + i - j, j] == a[i,j] (if upper form; i <= j) a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
where u is the number of bands above the diagonal.
Example of a_band (shape of a is (6,6), u=2)::
upper form: * * a02 a13 a24 a35 * a01 a12 a23 a34 a45 a00 a11 a22 a33 a44 a55
lower form: a00 a11 a22 a33 a44 a55 a10 a21 a32 a43 a54 * a20 a31 a42 a53 * *
Cells marked with * are not used.
Parameters ---------- a_band : (u+1, M) array_like The bands of the M by M matrix a. lower : bool, optional Is the matrix in the lower form. (Default is upper form) overwrite_a_band : bool, optional Discard data in a_band (may enhance performance) select : {'a', 'v', 'i'}, optional Which eigenvalues to calculate
====== ======================================== select calculated ====== ======================================== 'a' All eigenvalues 'v' Eigenvalues in the interval (min, max] 'i' Eigenvalues with indices min <= i <= max ====== ======================================== select_range : (min, max), optional Range of selected eigenvalues check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- w : (M,) ndarray The eigenvalues, in ascending order, each repeated according to its multiplicity.
Raises ------ LinAlgError If eigenvalue computation does not converge.
See Also -------- eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal matrices eigvals : eigenvalues of general arrays eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays eig : eigenvalues and right eigenvectors for non-symmetric arrays
Examples -------- >>> from scipy.linalg import eigvals_banded >>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]]) >>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]]) >>> w = eigvals_banded(Ab, lower=True) >>> w array([-4.26200532, -2.22987175, 3.95222349, 12.53965359]) """ return eig_banded(a_band, lower=lower, eigvals_only=1, overwrite_a_band=overwrite_a_band, select=select, select_range=select_range, check_finite=check_finite)
check_finite=True, tol=0., lapack_driver='auto'): """ Solve eigenvalue problem for a real symmetric tridiagonal matrix.
Find eigenvalues `w` of ``a``::
a v[:,i] = w[i] v[:,i] v.H v = identity
For a real symmetric matrix ``a`` with diagonal elements `d` and off-diagonal elements `e`.
Parameters ---------- d : ndarray, shape (ndim,) The diagonal elements of the array. e : ndarray, shape (ndim-1,) The off-diagonal elements of the array. select : {'a', 'v', 'i'}, optional Which eigenvalues to calculate
====== ======================================== select calculated ====== ======================================== 'a' All eigenvalues 'v' Eigenvalues in the interval (min, max] 'i' Eigenvalues with indices min <= i <= max ====== ======================================== select_range : (min, max), optional Range of selected eigenvalues check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. tol : float The absolute tolerance to which each eigenvalue is required (only used when ``lapack_driver='stebz'``). An eigenvalue (or cluster) is considered to have converged if it lies in an interval of this width. If <= 0. (default), the value ``eps*|a|`` is used where eps is the machine precision, and ``|a|`` is the 1-norm of the matrix ``a``. lapack_driver : str LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf', or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'`` and 'stebz' otherwise. 'sterf' and 'stev' can only be used when ``select='a'``.
Returns ------- w : (M,) ndarray The eigenvalues, in ascending order, each repeated according to its multiplicity.
Raises ------ LinAlgError If eigenvalue computation does not converge.
See Also -------- eigh_tridiagonal : eigenvalues and right eiegenvectors for symmetric/Hermitian tridiagonal matrices
Examples -------- >>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh >>> d = 3*np.ones(4) >>> e = -1*np.ones(3) >>> w = eigvalsh_tridiagonal(d, e) >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1) >>> w2 = eigvalsh(A) # Verify with other eigenvalue routines >>> np.allclose(w - w2, np.zeros(4)) True """ return eigh_tridiagonal( d, e, eigvals_only=True, select=select, select_range=select_range, check_finite=check_finite, tol=tol, lapack_driver=lapack_driver)
check_finite=True, tol=0., lapack_driver='auto'): """ Solve eigenvalue problem for a real symmetric tridiagonal matrix.
Find eigenvalues `w` and optionally right eigenvectors `v` of ``a``::
a v[:,i] = w[i] v[:,i] v.H v = identity
For a real symmetric matrix ``a`` with diagonal elements `d` and off-diagonal elements `e`.
Parameters ---------- d : ndarray, shape (ndim,) The diagonal elements of the array. e : ndarray, shape (ndim-1,) The off-diagonal elements of the array. select : {'a', 'v', 'i'}, optional Which eigenvalues to calculate
====== ======================================== select calculated ====== ======================================== 'a' All eigenvalues 'v' Eigenvalues in the interval (min, max] 'i' Eigenvalues with indices min <= i <= max ====== ======================================== select_range : (min, max), optional Range of selected eigenvalues check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. tol : float The absolute tolerance to which each eigenvalue is required (only used when 'stebz' is the `lapack_driver`). An eigenvalue (or cluster) is considered to have converged if it lies in an interval of this width. If <= 0. (default), the value ``eps*|a|`` is used where eps is the machine precision, and ``|a|`` is the 1-norm of the matrix ``a``. lapack_driver : str LAPACK function to use, can be 'auto', 'stemr', 'stebz', 'sterf', or 'stev'. When 'auto' (default), it will use 'stemr' if ``select='a'`` and 'stebz' otherwise. When 'stebz' is used to find the eigenvalues and ``eigvals_only=False``, then a second LAPACK call (to ``?STEIN``) is used to find the corresponding eigenvectors. 'sterf' can only be used when ``eigvals_only=True`` and ``select='a'``. 'stev' can only be used when ``select='a'``.
Returns ------- w : (M,) ndarray The eigenvalues, in ascending order, each repeated according to its multiplicity. v : (M, M) ndarray The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is the column ``v[:,i]``.
Raises ------ LinAlgError If eigenvalue computation does not converge.
See Also -------- eigvalsh_tridiagonal : eigenvalues of symmetric/Hermitian tridiagonal matrices eig : eigenvalues and right eigenvectors for non-symmetric arrays eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices
Notes ----- This function makes use of LAPACK ``S/DSTEMR`` routines.
Examples -------- >>> from scipy.linalg import eigh_tridiagonal >>> d = 3*np.ones(4) >>> e = -1*np.ones(3) >>> w, v = eigh_tridiagonal(d, e) >>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1) >>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4))) True """ d = _asarray_validated(d, check_finite=check_finite) e = _asarray_validated(e, check_finite=check_finite) for check in (d, e): if check.ndim != 1: raise ValueError('expected one-dimensional array') if check.dtype.char in 'GFD': # complex raise TypeError('Only real arrays currently supported') if d.size != e.size + 1: raise ValueError('d (%s) must have one more element than e (%s)' % (d.size, e.size)) select, vl, vu, il, iu, _ = _check_select( select, select_range, 0, d.size) if not isinstance(lapack_driver, string_types): raise TypeError('lapack_driver must be str') drivers = ('auto', 'stemr', 'sterf', 'stebz', 'stev') if lapack_driver not in drivers: raise ValueError('lapack_driver must be one of %s, got %s' % (drivers, lapack_driver)) if lapack_driver == 'auto': lapack_driver = 'stemr' if select == 0 else 'stebz' func, = get_lapack_funcs((lapack_driver,), (d, e)) compute_v = not eigvals_only if lapack_driver == 'sterf': if select != 0: raise ValueError('sterf can only be used when select == "a"') if not eigvals_only: raise ValueError('sterf can only be used when eigvals_only is ' 'True') w, info = func(d, e) m = len(w) elif lapack_driver == 'stev': if select != 0: raise ValueError('stev can only be used when select == "a"') w, v, info = func(d, e, compute_v=compute_v) m = len(w) elif lapack_driver == 'stebz': tol = float(tol) internal_name = 'stebz' stebz, = get_lapack_funcs((internal_name,), (d, e)) # If getting eigenvectors, needs to be block-ordered (B) instead of # matirx-ordered (E), and we will reorder later order = 'E' if eigvals_only else 'B' m, w, iblock, isplit, info = stebz(d, e, select, vl, vu, il, iu, tol, order) else: # 'stemr' # ?STEMR annoyingly requires size N instead of N-1 e_ = empty(e.size+1, e.dtype) e_[:-1] = e stemr_lwork, = get_lapack_funcs(('stemr_lwork',), (d, e)) lwork, liwork, info = stemr_lwork(d, e_, select, vl, vu, il, iu, compute_v=compute_v) _check_info(info, 'stemr_lwork') m, w, v, info = func(d, e_, select, vl, vu, il, iu, compute_v=compute_v, lwork=lwork, liwork=liwork) _check_info(info, lapack_driver + ' (eigh_tridiagonal)') w = w[:m] if eigvals_only: return w else: # Do we still need to compute the eigenvalues? if lapack_driver == 'stebz': func, = get_lapack_funcs(('stein',), (d, e)) v, info = func(d, e, w, iblock, isplit) _check_info(info, 'stein (eigh_tridiagonal)', positive='%d eigenvectors failed to converge') # Convert block-order to matrix-order order = argsort(w) w, v = w[order], v[:, order] else: v = v[:, :m] return w, v
"""Check info return value.""" if info < 0: raise ValueError('illegal value in argument %d of internal %s' % (-info, driver)) if info > 0 and positive: raise LinAlgError(("%s " + positive) % (driver, info,))
""" Compute Hessenberg form of a matrix.
The Hessenberg decomposition is::
A = Q H Q^H
where `Q` is unitary/orthogonal and `H` has only zero elements below the first sub-diagonal.
Parameters ---------- a : (M, M) array_like Matrix to bring into Hessenberg form. calc_q : bool, optional Whether to compute the transformation matrix. Default is False. overwrite_a : bool, optional Whether to overwrite `a`; may improve performance. Default is False. check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- H : (M, M) ndarray Hessenberg form of `a`. Q : (M, M) ndarray Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``. Only returned if ``calc_q=True``.
Examples -------- >>> from scipy.linalg import hessenberg >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]]) >>> H, Q = hessenberg(A, calc_q=True) >>> H array([[ 2. , -11.65843866, 1.42005301, 0.25349066], [ -9.94987437, 14.53535354, -5.31022304, 2.43081618], [ 0. , -1.83299243, 0.38969961, -0.51527034], [ 0. , 0. , -3.83189513, 1.07494686]]) >>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4))) True """ a1 = _asarray_validated(a, check_finite=check_finite) if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]): raise ValueError('expected square matrix') overwrite_a = overwrite_a or (_datacopied(a1, a))
# if 2x2 or smaller: already in Hessenberg if a1.shape[0] <= 2: if calc_q: return a1, numpy.eye(a1.shape[0]) return a1
gehrd, gebal, gehrd_lwork = get_lapack_funcs(('gehrd', 'gebal', 'gehrd_lwork'), (a1,)) ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a) _check_info(info, 'gebal (hessenberg)', positive=False) n = len(a1)
lwork = _compute_lwork(gehrd_lwork, ba.shape[0], lo=lo, hi=hi)
hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1) _check_info(info, 'gehrd (hessenberg)', positive=False) h = numpy.triu(hq, -1) if not calc_q: return h
# use orghr/unghr to compute q orghr, orghr_lwork = get_lapack_funcs(('orghr', 'orghr_lwork'), (a1,)) lwork = _compute_lwork(orghr_lwork, n, lo=lo, hi=hi)
q, info = orghr(a=hq, tau=tau, lo=lo, hi=hi, lwork=lwork, overwrite_a=1) _check_info(info, 'orghr (hessenberg)', positive=False) return h, q
""" Converts complex eigenvalues ``w`` and eigenvectors ``v`` to real eigenvalues in a block diagonal form ``wr`` and the associated real eigenvectors ``vr``, such that::
vr @ wr = X @ vr
continues to hold, where ``X`` is the original array for which ``w`` and ``v`` are the eigenvalues and eigenvectors.
.. versionadded:: 1.1.0
Parameters ---------- w : (..., M) array_like Complex or real eigenvalues, an array or stack of arrays
Conjugate pairs must not be interleaved, else the wrong result will be produced. So ``[1+1j, 1, 1-1j]`` will give a correct result, but ``[1+1j, 2+1j, 1-1j, 2-1j]`` will not.
v : (..., M, M) array_like Complex or real eigenvectors, a square array or stack of square arrays.
Returns ------- wr : (..., M, M) ndarray Real diagonal block form of eigenvalues vr : (..., M, M) ndarray Real eigenvectors associated with ``wr``
See Also -------- eig : Eigenvalues and right eigenvectors for non-symmetric arrays rsf2csf : Convert real Schur form to complex Schur form
Notes ----- ``w``, ``v`` must be the eigenstructure for some *real* matrix ``X``. For example, obtained by ``w, v = scipy.linalg.eig(X)`` or ``w, v = numpy.linalg.eig(X)`` in which case ``X`` can also represent stacked arrays.
.. versionadded:: 1.1.0
Examples -------- >>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]]) >>> X array([[ 1, 2, 3], [ 0, 4, 5], [ 0, -5, 4]])
>>> from scipy import linalg >>> w, v = linalg.eig(X) >>> w array([ 1.+0.j, 4.+5.j, 4.-5.j]) >>> v array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j], [ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j], [ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
>>> wr, vr = linalg.cdf2rdf(w, v) >>> wr array([[ 1., 0., 0.], [ 0., 4., 5.], [ 0., -5., 4.]]) >>> vr array([[ 1. , 0.40016, -0.01906], [ 0. , 0.64788, 0. ], [ 0. , 0. , 0.64788]])
>>> vr @ wr array([[ 1. , 1.69593, 1.9246 ], [ 0. , 2.59153, 3.23942], [ 0. , -3.23942, 2.59153]]) >>> X @ vr array([[ 1. , 1.69593, 1.9246 ], [ 0. , 2.59153, 3.23942], [ 0. , -3.23942, 2.59153]]) """ w, v = _asarray_validated(w), _asarray_validated(v)
# check dimensions if w.ndim < 1: raise ValueError('expected w to be at least one-dimensional') if v.ndim < 2: raise ValueError('expected v to be at least two-dimensional') if v.ndim != w.ndim + 1: raise ValueError('expected eigenvectors array to have exactly one ' 'dimension more than eigenvalues array')
# check shapes n = w.shape[-1] M = w.shape[:-1] if v.shape[-2] != v.shape[-1]: raise ValueError('expected v to be a square matrix or stacked square ' 'matrices: v.shape[-2] = v.shape[-1]') if v.shape[-1] != n: raise ValueError('expected the same number of eigenvalues as ' 'eigenvectors')
# get indices for each first pair of complex eigenvalues complex_mask = iscomplex(w) n_complex = complex_mask.sum(axis=-1)
# check if all complex eigenvalues have conjugate pairs if not (n_complex % 2 == 0).all(): raise ValueError('expected complex-conjugate pairs of eigenvalues')
# find complex indices idx = nonzero(complex_mask) idx_stack = idx[:-1] idx_elem = idx[-1]
# filter them to conjugate indices, assuming pairs are not interleaved j = idx_elem[0::2] k = idx_elem[1::2] stack_ind = () for i in idx_stack: # should never happen, assuming nonzero orders by the last axis assert (i[0::2] == i[1::2]).all(), "Conjugate pair spanned different arrays!" stack_ind += (i[0::2],)
# all eigenvalues to diagonal form wr = zeros(M + (n, n), dtype=w.real.dtype) di = range(n) wr[..., di, di] = w.real
# complex eigenvalues to real block diagonal form wr[stack_ind + (j, k)] = w[stack_ind + (j,)].imag wr[stack_ind + (k, j)] = w[stack_ind + (k,)].imag
# compute real eigenvectors associated with real block diagonal eigenvalues u = zeros(M + (n, n), dtype=numpy.cdouble) u[..., di, di] = 1.0 u[stack_ind + (j, j)] = 0.5j u[stack_ind + (j, k)] = 0.5 u[stack_ind + (k, j)] = -0.5j u[stack_ind + (k, k)] = 0.5
# multipy matrices v and u (equivalent to v @ u) vr = einsum('...ij,...jk->...ik', v, u).real
return wr, vr |