Making lobpcg support complex matrix
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
Hi Greg, that is already a ticket http://projects.scipy.org/scipy/ticket/452 Cheers, Nils On 10/11/12, Gregory Crosswhite <g.crosswhite@uq.edu.au> wrote:
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 _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
Gregory Crosswhite <g.crosswhite <at> uq.edu.au> writes: [clip]
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.
The original LOBPCG paper IIRC talked only about real symmetric problems, but it seems other existing implementations accept also Hermitian input. So yes, a patch fixing the implementation in Scipy for complex Hermitian matrices (plus adding tests) would be welcome. The easiest way (for us, not necessarily for you at least the first time :) to submit it would be via a pull request on Github; not a requirement though. As Nils said, there's a ticket with at least some work done in this direction: http://projects.scipy.org/scipy/ticket/452 I'm not sure what is the exact status with that, but you may want to take a look if it is of some help. -- Pauli Virtanen
participants (3)
-
Gregory Crosswhite -
Nils Wagner -
Pauli Virtanen