I need to solve symmetric generalized eigenvalue problems with large, sparse stiffness and mass matrices, say A and B. The problem is of the form A*v = lambda *BV.
They are both Hermitian, and B is positive definite. Both are approximately of size 70 000 x 70 000. A has about 5 million non- zero values and B has about 1.6 million. I'm using them both in Matlab's sparse matrix format. A has condition number on the order of 10^9 and B has a condition number on the order of 10^6.
I had no problems getting the correct values/vectors with eigs. However, I would like to solve this problem using the lobpcg function, available from Matlab Central's File Exchange website. This is partly because I'm interested in solving this problem in Python, where I have to use a function identical to this lobpcg.
I tried this first on a scaled down eigenvalue problem, using a = A(1:10000, 1:10000) and b = B(1:10000,1:10000), finding the 10 smallest eigenvalues. eigs and lobpcg returned similar results for these values;
Here is the syntax I used for the two:
X = rand(68394, 10)
[v, lambda] = lobpcg(X, a, b, 1e-5, 20,0)
[V, D] = eigs(a,b, 10, 'sm')
I could not get lobpcg to work for full sized A and B, no matter how many iterations I used. lobpcg requires an initial "guess" for the eigenvectors ( X ), a matrix whose columns determine the number of eigenvalues desired. I made a plot (link below) for an eigenvalue computed correctly with eigs ("Matlab" in the legend), and how the number of iterations affects the result using lobpcg. For lobpcg, I tried 4 different guesses for X, where X1 to X3 are random matrices and Xorthornormal is an orthonormalized random matrix.
The plot shows that an unreasonable number of iterations would be necessary to get the correct values from lobpcg.