[SciPy-Dev] Making lobpcg support complex matrix
Wed Oct 10 22:46:33 CDT 2012
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.
More information about the SciPy-Dev