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