Hi,
in econometrics/statistics it's relatively common to solve a special
generalized eigenvalue problem, where the matrices are real, symmetric,
and positive definite. It bothered me to depend on scipy just because of
that. Also, scipy.linalg.eig is much more general and returns complex
types which is inconvenient for sorting etc. So I've written the
workaround function below, but it's probably not very efficient (at
least I hope the algebra is right).
I'd be happy to hear your comments, suggestions; first with respect to
the function itself, but also whether you agree that the implemented
functionality would be useful enough to include in official numpy (in a
more professional way than I have done).
cheers,
Sven
def geneigsympos(A, B):
""" Solves symmetric-positive-definite generalized
eigenvalue problem Az=lBz.
Takes two real-valued symmetric matrices A and B (B must also be
positive-definite) and returns the corresponding generalized (but
also real-valued) eigenvalues and eigenvectors. Return format: as
in scipy.linalg.eig, tuple (l, Z); l is directly taken from eigh
output (a 1-dim array of length A.shape[0] ?), and Z is an array
or matrix (depending on type of input A) with the corresponding
eigenvectors in columns (hopefully).
Eigenvalues are ordered ascending (thanks to eigh).
"""
from numpy import asmatrix, asarray, linalg
#fixme: input checks on the matrices
LI = asmatrix(linalg.cholesky(B)).I
C = LI * asmatrix(A) * LI.T
evals, evecs = linalg.eigh(C)
if type(A) == type(asarray(A)): output = "array"
# A was passed as numpy-array
else: output = "matrix"
#but the evecs need to be transformed back:
evecs = LI.T * asmatrix(evecs)
if output == "array": return evals, asarray(evecs)
else: return asmatrix(evals), evecs