[Numpy-discussion] (some) generalized eigenvalues in numpy

[Sven Schreiber]

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