Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add Eigenvalue Range Functionality for Symmetric Eigenvalue Problems #6510

Closed
wants to merge 7 commits into from
Closed

Add Eigenvalue Range Functionality for Symmetric Eigenvalue Problems #6510

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
48cf326
Update -syevr routines to take in vl and vu
apapanico Aug 22, 2016
21d97f6
Monkey patch for -evr methods
apapanico Aug 22, 2016
73bba9f
Tests for monkey patched -evr methods
apapanico Aug 22, 2016
bd80c76
Fixing data type in -heevr
apapanico Aug 24, 2016
b66593a
Compute eigenvalues within a specified range
apapanico Aug 25, 2016
ca11239
Tests for parameter `eigrng` in `eigh`
apapanico Aug 25, 2016
da1f57e
PEP8 formatting
apapanico Aug 25, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
34 changes: 27 additions & 7 deletions scipy/linalg/decomp.py
View file
Open in desktop
Expand Up @@ -7,6 +7,7 @@
# additions by Bart Vandereycken, June 2006
# additions by Andrew D Straw, May 2007
# additions by Tiziano Zito, November 2008
# additions by Alex Papanicolaou, August 2016
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You don't necessarily need to add a comment here --- who modified what is kept track of by the version contorl.
In particular, the file has been modified by many people since 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
Expand Down Expand Up @@ -240,7 +241,7 @@ def eig(a, b=None, left=False, right=True, overwrite_a=False,

def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
overwrite_b=False, turbo=True, eigvals=None, type=1,
check_finite=True):
check_finite=True, eigrng=None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer to write eig_range or just range instead of abbreviations.

"""
Solve an ordinary or generalized eigenvalue problem for a complex
Hermitian or real symmetric matrix.
Expand Down Expand Up @@ -272,7 +273,8 @@ def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
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.
If eigvals and eigrng are omitted, all eigenvalues and eigenvectors
are returned.
type : int, optional
Specifies the problem type to be solved:

Expand All @@ -289,6 +291,10 @@ def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
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.
eigrng : tuple (vl, vu), optional
The half-open interval (vl,vu] on which eigenvalues and corresponding
eigenvectors are to be returned: vl < vu. If eigrng and eigvals are
omitted, all eigenvalues and eigenvectors are returned.

Returns
-------
Expand Down Expand Up @@ -352,14 +358,21 @@ def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
_job = (eigvals_only and 'N') or 'V'

# port eigenvalue range from python to fortran convention
if eigvals is not None and eigrng is not None:
raise ValueError("Only one of eigvals or eigrng may be set")
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'
raise ValueError('The eigenvalue index range specified is not valid.\n'
'Valid range is [%s,%s]' % (0, a1.shape[0]-1))
lo += 1
hi += 1
eigvals = (lo, hi)
if eigrng is not None:
vl, vu = eigrng
if vl >= vu:
raise ValueError('The eigenvalue value range specified is not valid:\n'
'(%s,%s]' % (vl, vu))

# set lower
if lower:
Expand All @@ -379,14 +392,21 @@ def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
# for all lapack routines
if b1 is None:
(evr,) = get_lapack_funcs((pfx+'evr',), (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:
if eigvals is not None:
(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]
elif eigrng is not None:
(vl, vu) = eigrng
w_full, v_full, info = evr(a1, uplo=uplo, jobz=_job, range="V", vl=vl,
vu=vu, overwrite_a=overwrite_a)
ix = nonzero(w_full)[0].tolist()
Copy link
Member

@pv pv Sep 17, 2016 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if there are zero eigenvalues?

Is the number of eigenvalues found returned by the LAPACK routine in the M parameter?
If not, how is it possible to know how many were found?
What is the purpose of ISUPPZ parameter in the dsyevr?

The nonzero and tolist are probably not needed.

w = w_full[ix]
v = v_full[:, ix]
else:
w, v, info = evr(a1, uplo=uplo, jobz=_job, range="A", il=1,
iu=a1.shape[0], overwrite_a=overwrite_a)

# Generalized Eigenvalue Problem
else:
Expand Down
105 changes: 25 additions & 80 deletions scipy/linalg/flapack.pyf.src
View file
Open in desktop
Expand Up @@ -6,7 +6,7 @@
! $Revision$ $Date$
!
! Additions by Travis Oliphant, Tiziano Zito, Collin RM Stocks, Fabian Pedregosa
! Skipper Seabold
! Skipper Seabold, Alex Papanicolaou
!
! <prefix2=s,d> <ctype2=float,double> <ftype2=real,double precision> <wrap2=ws,d>
! <prefix2c=c,z> <ftype2c=complex,double complex> <ctype2c=complex_float,complex_double> <wrap2c=wc,z>
Expand Down Expand Up @@ -2531,122 +2531,67 @@ end subroutine <prefix>gbtrs
!
! RRR routines for standard eigenvalue problem
!
subroutine ssyevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)

subroutine <prefix2>syevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
! Standard Eigenvalue Problem
! simple/expert driver: all eigenvectors or optionally selected eigenvalues
! algorithm: Relatively Robust Representation
! matrix storage
! Real - Single precision
! Real - Single/Double precision
character intent(in) :: jobz='V'
character intent(in) :: range='A'
character intent(in) :: uplo='L'
integer intent(hide) :: n=shape(a,0)
real intent(in,copy,aligned8),dimension(n,n) :: a
<ftype2> intent(in,copy,aligned8),dimension(n,n) :: a
integer intent(hide),depend(n,a) :: lda=n
real intent(hide) :: vl=0
real intent(hide) :: vu=1
<ftype2> optional,intent(in) :: vl=0
<ftype2> optional,intent(in) :: vu=1
integer optional,intent(in) :: il=1
integer optional,intent(in),depend(n) :: iu=n
real intent(hide) :: abstol=0.
<ftype2> intent(hide) :: abstol=0.
integer intent(hide),depend(iu) :: m=iu-il+1
real intent(out),dimension(n),depend(n) :: w
real intent(out),dimension(n,m),depend(n,m) :: z
<ftype2> intent(out),dimension(n),depend(n) :: w
<ftype2> intent(out),dimension(n,m),depend(n,m) :: z
integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
integer intent(hide),dimension(2*m) :: isuppz
integer intent(in),depend(n) :: lwork=26*n
real intent(hide),dimension(lwork) :: work
<ftype2> intent(hide),dimension(lwork) :: work
integer intent(hide),depend(n):: liwork=10*n
integer intent(hide),dimension(liwork) :: iwork
integer intent(out) :: info
end subroutine ssyevr
subroutine dsyevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
end subroutine <prefix2>syevr

subroutine <prefix2c>heevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,rwork,lrwork,iwork,liwork,info)
! Standard Eigenvalue Problem
! simple/expert driver: all eigenvectors or optionally selected eigenvalues
! algorithm: Relatively Robust Representation
! matrix storage
! Real - Double precision
! Complex - Single/Double precision
character intent(in) :: jobz='V'
character intent(in) :: range='A'
character intent(in) :: uplo='L'
integer intent(hide) :: n=shape(a,0)
double precision intent(in,copy,aligned8),dimension(n,n) :: a
<ftype2c> intent(in,copy,aligned8),dimension(n,n) :: a
integer intent(hide),depend(n,a) :: lda=n
double precision intent(hide) :: vl=0
double precision intent(hide) :: vu=1
<ftype2> optional,intent(in) :: vl=0
<ftype2> optional,intent(in) :: vu=1
integer optional,intent(in) :: il=1
integer optional,intent(in),depend(n) :: iu=n
double precision intent(hide) :: abstol=0.
<ftype2> intent(hide) :: abstol=0.
integer intent(hide),depend(iu) :: m=iu-il+1
double precision intent(out),dimension(n),depend(n) :: w
double precision intent(out),dimension(n,m),depend(n,m) :: z
<ftype2> intent(out),dimension(n),depend(n) :: w
<ftype2c> intent(out),dimension(n,m),depend(n,m) :: z
integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
integer intent(hide),dimension(2*m) :: isuppz
integer intent(in),depend(n) :: lwork=26*n
double precision intent(hide),dimension(lwork) :: work
integer intent(hide),depend(n):: liwork=10*n
integer intent(hide),dimension(liwork) :: iwork
integer intent(out) :: info
end subroutine dsyevr
subroutine cheevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,rwork,lrwork,iwork,liwork,info)
! Standard Eigenvalue Problem
! simple/expert driver: all eigenvectors or optionally selected eigenvalues
! algorithm: Relatively Robust Representation
! matrix storage
! Complex - Single precision
character intent(in) :: jobz='V'
character intent(in) :: range='A'
character intent(in) :: uplo='L'
integer intent(hide) :: n=shape(a,0)
complex intent(in,copy,aligned8),dimension(n,n) :: a
integer intent(hide),depend(n,a) :: lda=n
real intent(hide) :: vl=0
real intent(hide) :: vu=1
integer optional,intent(in) :: il=1
integer optional,intent(in),depend(n) :: iu=n
real intent(hide) :: abstol=0.
integer intent(hide),depend(iu) :: m=iu-il+1
real intent(out),dimension(n),depend(n) :: w
complex intent(out),dimension(n,m),depend(n,m) :: z
integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
integer intent(hide),dimension(2*m) :: isuppz
integer intent(in),depend(n) :: lwork=18*n
complex intent(hide),dimension(lwork) :: work
<ftype2c> intent(hide),dimension(lwork) :: work
integer intent(hide),depend(n) :: lrwork=24*n
real intent(hide),dimension(lrwork) :: rwork
<ftype2> intent(hide),dimension(lrwork) :: rwork
integer intent(hide),depend(n):: liwork=10*n
integer intent(hide),dimension(liwork) :: iwork
integer intent(out) :: info
end subroutine cheevr
subroutine zheevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,rwork,lrwork,iwork,liwork,info)
! Standard Eigenvalue Problem
! simple/expert driver: all eigenvectors or optionally selected eigenvalues
! algorithm: Relatively Robust Representation
! matrix storage
! Complex - Double precision
character intent(in) :: jobz='V'
character intent(in) :: range='A'
character intent(in) :: uplo='L'
integer intent(hide) :: n=shape(a,0)
complex*16 intent(in,copy,aligned8),dimension(n,n) :: a
integer intent(hide),depend(n,a) :: lda=n
double precision intent(hide) :: vl=0
double precision intent(hide) :: vu=1
integer optional,intent(in) :: il=1
integer optional,intent(in),depend(n) :: iu=n
double precision intent(hide) :: abstol=0.
integer intent(hide),depend(iu) :: m=iu-il+1
double precision intent(out),dimension(n),depend(n) :: w
complex*16 intent(out),dimension(n,m),depend(n,m) :: z
integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
integer intent(hide),dimension(2*m) :: isuppz
integer intent(in),depend(n) :: lwork=18*n
complex*16 intent(hide),dimension(lwork) :: work
integer intent(hide),depend(n) :: lrwork=24*n
double precision intent(hide),dimension(lrwork) :: rwork
integer intent(hide),depend(n):: liwork=10*n
integer intent(hide),dimension(liwork) :: iwork
integer intent(out) :: info
end subroutine zheevr
end subroutine <prefix2c>heevr

subroutine ssygv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,info)
! Generalized Eigenvalue Problem
! simple driver (all eigenvectors)
Expand Down
67 changes: 67 additions & 0 deletions scipy/linalg/lapack.py
View file
Open in desktop
Expand Up @@ -404,6 +404,73 @@
_flapack.sgegv = sgegv
_flapack.zgegv = zgegv

# Patching EVR Methods to handle value ranges
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should say that this is used to retain backward compatibility.

_evr_doc = """
w,z,info = {name}(a,[jobz,range,uplo,il,iu,lwork,overwrite_a,vl,vu])

Wrapper for ``{name}``.

Parameters
----------
a : input rank-2 array({type!r}) with bounds (n,n)

Other Parameters
----------------
jobz : input string(len=1), optional
Default: 'V'
range : input string(len=1), optional
Default: 'A'
uplo : input string(len=1), optional
Default: 'L'
overwrite_a : input int, optional
Default: 0
il : input int, optional
Default: 1
iu : input int, optional
Default: n
lwork : input int, optional
Default: 26*n
vl : input float, optional
Default: 0
vu : input float, optional
Default: 1

Returns
-------
w : rank-1 array({w_type!r}) with bounds (n)
z : rank-2 array({type!r}) with bounds (n,m)
info : int
"""


def evr_decorator(evr_fun, name, type):

_arg_list = ['a', 'jobz', 'range', 'uplo', 'il',
'iu', 'lwork', 'overwrite_a', 'vl', 'vu']

def wrapper(*args, **kwargs):
_new_kwargs = kwargs.copy()
for k, arg in zip(_arg_list, args):
_new_kwargs[k] = arg
# Fix bug in -EVR where il/iu are not ignored
if _new_kwargs['range'] != 'I':
_new_kwargs['il'] = 1
_new_kwargs['iu'] = min(_new_kwargs['a'].shape)
return evr_fun(**_new_kwargs)

wrapper.__module__ = 'scipy.linalg._flapack'
wrapper.__doc__ = _evr_doc.format(
name=name, w_type=type.lower(), type=type)
wrapper.__name__ = name
Copy link
Member

@pv pv Sep 17, 2016 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should also copy the _cpointer attribute for full backward compat.


return wrapper

ssyevr = evr_decorator(ssyevr, 'ssyevr', 'f')
dsyevr = evr_decorator(dsyevr, 'dsyevr', 'd')
cheevr = evr_decorator(cheevr, 'cheevr', 'F')
zheevr = evr_decorator(zheevr, 'zheevr', 'D')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need also replace the function in _flapack:

_flapack.ssyevr = ssyevr



# some convenience alias for complex functions
_lapack_alias = {
'corghr': 'cunghr', 'zorghr': 'zunghr',
Expand Down
50 changes: 50 additions & 0 deletions scipy/linalg/tests/test_decomp.py
View file
Open in desktop
Expand Up @@ -653,6 +653,25 @@ def test_eigh():
turbo, eigenvalues)


def test_eigh_rng():
DIM = 6
v = {'dim': (DIM,),
'dtype': ('f','d','F','D'),
'overwrite': (True, False),
'lower': (True, False),
'eigrng': (None, (-.5, .5), (-.75, .75))}

for dim in v['dim']:
for typ in v['dtype']:
for overwrite in v['overwrite']:
for eigenrange in v['eigrng']:
for lower in v['lower']:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that itertools.product can be used to avoid deeply nested loops:

dims = (DIM,)
dtypes = ('f', 'd', 'F', 'D')
...
parts = [dims, dtypes, ...]
for dim, typ, overwrite, eigenrange, lower in itertools.product(*parts):
    yield (...)

yield (eigenhproblem_range,
'ordinary',
dim, typ, overwrite, lower,
eigenrange)


def test_eigh_of_sparse():
# This tests the rejection of inputs that eigh cannot currently handle.
import scipy.sparse
Expand Down Expand Up @@ -715,6 +734,37 @@ def eigenhproblem_general(desc, dim, dtype,
assert_array_almost_equal(diag2_, ones(diag2_.shape[0]), DIGITS[dtype])


def eigenhproblem_range(desc, dim, dtype,
overwrite, lower,
eigenrange):
"""Solve a standard eigenvalue problem."""
if iscomplex(empty(1, dtype=dtype)):
a = _complex_symrand(dim, dtype)
else:
a = symrand(dim).astype(dtype)

if overwrite:
a_c = a.copy()
else:
a_c = a
w, z = eigh(a, overwrite_a=overwrite, lower=lower, eigrng=eigenrange)
assert_dtype_equal(z.dtype, dtype)
assert_equal(len(np.nonzero(w)[0]), len(w))
assert_equal(z.shape[1], len(w))
w = w.astype(dtype)
diag_ = diag(dot(z.T.conj(), dot(a_c, z))).real
assert_array_almost_equal(diag_, w, DIGITS[dtype])


def test_eigh_range_invalid():
# This tests the rejection of inputs that eigh cannot currently handle.
dim = 6
dtype = 'd'
a = symrand(dim).astype(dtype)
assert_raises(ValueError, eigh, a, eigvals=True, eigrng=True)
assert_raises(ValueError, eigh, a, eigrng=(1, 0))


def test_eigh_integer():
a = array([[1,2],[2,7]])
b = array([[3,1],[1,5]])
Expand Down