[SciPy-Dev] Making lobpcg support complex matrix

[Gregory Crosswhite]

Dear all,

At the moment lobpcg seems to have been designed with only real numbers 
in mind;  this is unfortunate because I have been having some trouble 
with eigs and was hoping to try it out as an alternative.

Fortunately, it looks to me like the solution is rather simple.  I was 
able to get it to work by replacing all instances of ".T" with 
".T.conj()" --- which meant that all of the dot products were now 
correct for complex numbers --- and by replacing the cast to "float64" 
with a cast to "complex128".  Changing the second cast is of course less 
than ideal in the cases where "float64" is indeed what is being used, 
but something like it is needed to make the answers not be garbage for 
complex input matrices.

With these changes, I generated 10x10 complex Hermitian matrices A and B 
to serve as respectively the problem matrix and the normalization/metric 
matrix, and ran lobpcg with a tolerance of 1e-10.  For the standard 
eigenvalue problem with k=2 (and random X) lobpcg got essentially the 
same answers as eigh for the lowest two eigenvalue/eigenvector pairs, 
and for the generalized eigenvalue problem it got the same eigenvalues 
as eigh and it got eigenvectors whose entries were within ~10^-4 of the 
eigenvectors returned by eigh (after dividing by the first entry to make 
the vectors proportionately the same) and such that the normalized 
overlap between the eigenvectors of lobpcg and eigh (using B as the 
metric) was within ~ 10^-8 of 1 (not surprising as this is the square of 
the first number).  I ran a quick check where I dropped the imaginary 
parts of A and B and re-ran this analysis and say the errors fall to 
respectively ~ 10^-6 and 10^-12, so the algorithm gets less accurate 
results for complex numbers than for real numbers, though I don't know 
enough about how it works to speculate on which this would be.

Anyway, I don't know much about lobpcg so there might be some issue that 
I've missed in simply adding ".conj()" everywhere a ".T" appeared to fix 
the dot products.  I do think it would be very nice, though, to be able 
to use lobpcg for complex matrices, and so I would be willing to submit 
a patch towards this end.

Thoughts?

Cheers,
Greg