search

gabrielasd / pyeom

GNU General Public License v3.0
2 stars 0 forks source link

Eigenvalue problem solvers #53

Closed gabrielasd closed 2 years ago

gabrielasd commented 3 years ago

Currently, we have two options to transform the EOM generalized eigenvalue problem into a standard one (both defined in the solve_dense method of the base module). A c = w B c ---> A' d = w d

One is to use symmetric orthogonalization and the other is to simply multiply by the inverse of the right-hand-side matrix (RHS) on both sides of the equation. In both cases, numerical instabilities are handled by setting to zero the inverse of close to zero eigenvalues (by some threshold value). Singular value decomposition is used to find the right-hand-side matrix eigenvalues.

This is one way of neglecting the eigenvectors in the RHS with close to zero eigenvalues. However, the final solutions (w) have multiple zero values that either correspond to real zero energy transitions or originate from the RHS. To distinguish between these zero values we can modify the solve_dense method to explicitly discard these eigenvectors of the RHS and solve the "new" eigenvalue problem expanded on the reduced basis.

Additionally, the nature of the RHS (positive definite/semidefinite, other) depends on the type of excited state that is being modeled and the form of the EOM equation used. It may be appropriate to include an option in the implemented EOM methods that sets a default solver.

PaulWAyers commented 3 years ago

I think having an option to choose a solver. It would allow (in the long run) a move to sparse solvers or other similar methods also.

Choosing between the "reduced space" and "whole space" methods is probably something we should do based on the results when we can do a more systematic study but then we can set up a default option for the way the eigenproblem is approached.

msricher commented 3 years ago

Sometimes (when the inappropriate solver is being used, most likely) numpy complains about dividing by zero.

We should make use of Python's warnings module and np.seterr to handle this gracefully. @gabrielasd I can help, if you can decide what behaviour you want in the case of div/0.

msricher commented 3 years ago

I think having an option to choose a solver. It would allow (in the long run) a move to sparse solvers or other similar methods also.

Choosing between the "reduced space" and "whole space" methods is probably something we should do based on the results when we can do a more systematic study but then we can set up a default option for the way the eigenproblem is approached.

@PaulWAyers The dense method is about as good as it gets, because it's a numpy call to C. We should keep it though, and give the user the option to choose, although the default should probably be the sparse solver.

If SciPy's sparse solver ever becomes too slow (it's a lot of Python wrapping Arpack), I can always rip out the eigensolver from PyCI (Spectralib) and use it here too.

PaulWAyers commented 3 years ago

Sounds good. I think the dense solver will get us as far as we need for now, but having the flexibility is great.

msricher commented 3 years ago

branch floaterror should be the final commit we need to merge to close this issue @gabrielasd

msricher commented 3 years ago

I'm actually not sure if the default behaviour should be "warn". It's a known behaviour that's accounted for, so maybe it should just be "ignore" so it proceeds quietly. Any feedback @gabrielasd ?

gabrielasd commented 3 years ago

@msricher I think having "ignore" as the default is OK

gabrielasd commented 3 years ago

branch floaterror should be the final commit we need to merge to close this issue @gabrielasd

commits 1cc11d25fb4f1b815b44cdb725fc03805582b8ef and 47625acfaebc15fd932c68ac669d535e6fff8845