scipy
ENH: improved accuracy and orthonormality of output eigenvectors in lobpcg
Reference issue
Closes gh-16319, gh-16408. gh-15461, gh-15935
What does this implement/fix?
- Introduced a free and simple tunable restart to potentially increase the attainable accuracy for edge cases
- Internal postprocessing in lobpcg that runs one final exact RR giving more accurate and orthonormal eigenvectors
- Output lists of history are now in np array format, to improve performance a bit
- In cases of iteration failure or not reaching the required tolerance in the given number of iterations, the code now outputs the computed iterate with the smallest max norm of the residual and drops the history of subsequent iterations.
- Some minor code beatification and restructuring, not affecting the iterates.
- Replace the internal matrix converter into the LinearOperator with a simple function, removing the unnecessary dependence.
- Remove the check for LinearOperator format input and thus allow a simple function handle of a callable object as an input.
- Better consistency in using dtypes for different objects
- Multiple warnings added for the user to easier see how to fix edge cases or violating the assumptions crucial for the algorithm health
- Fool-proof enhancements, attempting to guess and fix user errors in the input, rather than letting the execution to break
@tylerjereddy - thanks for reviewing!
- I don't know that I saw tests added for some of the new error handling scenarios (or how important those tests would be); some new warnings do seem covered I think
Yes. It's lots of nontrivial work to come up with examples initiating every new error handling code. I have added only those unit tests that are relatively simple to grenerate. It's still an improvement compared to the original code without this error handling.
- it looked like there was some positive feedback on the changes here from the JAX team at Google here: #16408 (comment); but then they noticed a tradeoff here: #16408 (comment); am I reading that correctly? can you summarize that situation for non-expert?
There is no tradeoff. The reviewer there and I in my tests get different accuracy for one of the examples of this PR, but the original code simply fails on this example to my recollection, so the PR brings an improvement either way.
One of the main goals of this PR is very simple - without changing the algorithm make the code as fool-proof as possible, e.g., even in the worst-case scenario of immediate failure just returning the unchanged initial approximations with a proper warning rather than an error. Allowing the user to call function handles instead of LinearOperators is extremely flexible and convenient but opens a whole new can of warms allowing gross abuses of the underlying assumptions on the input matrices making the algorithm to fail in different parts of the code.
I'm not super familiar with the algorithm but this looks like an improvement to me. I'm wondering if there are any changes in performance, but also not sure how easy those would be to check. Overall I'm +1 - thanks @lobpcg !
@melissawm thanks for reviewing!
The accuracy/stability improvement due to outputting the best rather than the last approximation is documented in https://github.com/scipy/scipy/issues/16408#issuecomment-1171165201
The accuracy improvement due to adding the postprocessing step should hopefully be seen for example in another edit - removing orthogonalization of eigenvectors in _svds.py calling lobpcg and still passing all the unit tests for svds, which I am about to make.
The postprocessing step adds a bit of new operations to run but the extra time cost is smaller than a cost of running an extra iteration so negligible and no new memory is needed.
@lobpcg it would also be nice if you could propose a note on the improvements to lobpcg to add to the 1.10.0 release notes, in a comment here.
@lobpcg it would also be nice if you could propose a note on the improvements to
lobpcgto add to the 1.10.0 release notes, in a comment here.
@rgommers What does this implement/fix summary of main new features: A simple tunable restart potentially increases the attainable accuracy for edge cases. Internal postprocessing runs one final exact Rayleigh-Ritz method giving more accurate and orthonormal eigenvectors. Output the computed iterate with the smallest max norm of the residual and drop the history of subsequent iterations. Remove the check for LinearOperator format input and thus allow a simple function handle of a callable object as an input. Fool-proof enhancements, attempting to guess and fix user errors in the input, rather than letting the execution to break.