# scipy.linalg.sparse.eigsh does not work for generalised eigenvalues

I asked this question over at StackOverflow and someone told me that I'd get a better answer here. So here's my problem:

I'm working on a machine learning project which involves doing a Principal Component Analysis on some labeled data and using those labels to extract more valuable information from the data.

To do that, I'm calculating a scatter matrix for each class, and for each pair of classes I need to solve a generalised eigenvalue problem for their scatter matrices, as follows:

$$S_i v = w (S_j + \beta I) v,$$

where $\beta$ is a multiplier and $I$ is the identity matrix. Now, this is the code in python:

        jeigenvalues = eigsh(scatter_j, k=10, return_eigenvectors=False, maxiter=100)
beta = betaMult*mean(jeigenvalues)
print(beta)
print(scatter_j+beta*eye(shape(x_data)[1]))
w, v = eigsh(scatter_i,M=scatter_j+beta*eye(shape(x_data)[1]),k=int(numberOfEVs/45), maxiter=100)
print(i,j,'done')


numberOfEVs is 90 in my current code (so that it's divisible by 45).

But the problem is, at the line where I use the eigsh for the aforementioned formula, it never gives me an answer. It keeps eating more and more memory without even completing a single iteration (I set its maxiter input to 1, and it still didn't give an answer). When I don't give the eigsh function the M argument (which is the matrix on the right side of the generalised EV problem and it is assumed to be "I" when not specified), it works correctly. But when M is provided, it becomes unresponsive.

Any ideas?

EDIT: The scatter matrices have rather small entries, mostly around 10^-5. I've also tried multiplying the left hand side by the inverse of the RHS matrix, and again it's having the same issue (goes on for a long time without an answer). Is the smallness of these entries the issue? How can I solve it, then?

• Welcome to SciComp.SE! How large are your matrices $S_i$? Are they indeed sparse and symmetric? For a generalized eigenvalue problem, the matrix M must be symmetric and positive definite. Maybe your shift $\beta$ is not large enough if $S_j$ is not s.p.d.? – Christian Clason Jan 21 '16 at 12:31
• The algorithm behind eigsh uses implicitly restarted Arnoldi for the standard eigenvalue problem $M^{-1}Ax = \lambda x$, i.e., it requires the solution of $Mx = y$ in every iteration. For this, eigsh computes an LU factorization of $M$, which is probably what's taking so long. You might get better performance (or at least an idea what's going wrong) by precomputing such a factorization and passing a function that uses this factorization to solve the system as the Minv argument. – Christian Clason Jan 21 '16 at 12:39
• Thank you! The matrices are scatter matrices, so they're both symmetric and positive definite. I'll try out the solution you told me, but I also tried the inv function (for sparse) and that too displayed the same kind of behaviour. I've tried changing beta (it's been both ~0.001 and ~0.00001, with the same behaviour from the eigsh function). – Ramtin Yazdanian Jan 21 '16 at 14:31
• inv is never a good idea for sparse matrices. Your range of $\beta$ might still be too small. You could try to compute the smallest eigenvalue of $S_j$ (e.g., using this trick) to see how large of a shift you need. – Christian Clason Jan 21 '16 at 14:36
• Thank you for your advice. Increasing beta finally worked: the algorithm takes ~25 minutes for each eigsh and takes a couple of GBs of memory, but it finally gives an answer. – Ramtin Yazdanian Jan 22 '16 at 17:37