Additional GEVP method with errors #195

JanNeuendorf · 2023-06-15T17:30:23Z

This adds an additional method to the correlator class called error_gevp().
It uses cholesky decomposition to give the eigenvectors with errors in the same form as the output of Corr.GEVP().
This could have been implemented as an extra flag for Corr.GEVP(), but I think this would not be optimal for the following reasons.

It would need to work with the different sorting methods and all changes to them.
The output might differ due to numerical issues and the different solver. The user might turn on errors and suddenly see a different result.
The presence of an error flag suggests to the user that they should always leave it on for correct error propagation. That is likely not what they want.

If you do not want this as a method on Corr , this could just be moved to misc?

…errors.

fjosw · 2023-06-15T17:34:48Z

Thanks Jan, this looks good. I will have a closer look tomorrow. Maybe @s-kuberski also has an opinion?

As mentioned in #191 this feature will break with numpy 1.25 due to an incompatability in autograd. Another reason to fix this behavior upstream.

fjosw

I left a few small comments in the code.

Can we also have some consistency checks versus the standard GEVP method? You mention in the comments that the results do not necessarily agree but can we have some check within a tolerance?

pyerrors/correlators.py

…ic gamma method

JanNeuendorf · 2023-06-16T13:29:58Z

I added a test that shows that the mean values of the vectors with error are close to the ones from the other method. It also checks that there are errors, that they increase the errors of the projected correlator and, finally, that the result is None for a ts where the correlator is None.

fjosw · 2023-06-16T13:34:05Z

Very good, thanks a lot. I would wait for @s-kuberski's feedback but from my side, this is ready to merge.

Another question is maybe how to interpret the fact that the error of the projected correlator is larger when including the error of the eigenvectors. Do you have a feeling if one neglects anything when using the original GEVP method?

JanNeuendorf · 2023-06-16T13:45:38Z

That is an interesting question. We had a discussion about this before. My opinion is that it is correct to project without errors. If I measure a correlator with some source, I get some linear combination of states. If I measure with multiple sources and project with some (error-less) vector, I get another linear combination of states. Ideally, i got this vector from the GEVP and it produces a better overlap with the state I care about, but this is something that will be reflected in my later analysis. For example, I get longer plateaus. The GEVP-error does not help to quantify state assignment. I either found a projection to the state I want, or I did not. I feel like this error should only be used when we are taking about the vector or source itself.

s-kuberski · 2023-06-19T15:54:49Z

Hi! Thanks for implementing this. I try to have a look as soon as possible, sorry for being slow.

s-kuberski · 2023-06-23T10:25:14Z

Hi,

first, thank you very much for taking the time to implement this. Also, I think that the perspective that you take in the last comment is an interesting one and I would agree concerning the fact that one does not need to include the uncertainties on the eigenvectors if one is just interested in having a nice projection.

Concerning the changes:

Conceptually, I would really prefer to have more overlap between the GEVP and the error_GEVP methods or even a combination. You write that the user could be tempted to switch on the errors, but I think that is not necessarily the case if you set the default to no error and carefully explain the flag in the docstring.
I agree that in this case one should have matching central values. But why not using the same algorithm to solve the GEVP in both cases? As I see it, you could just implement your new setup in _GEVP_solver, taking the correct routines for either scalars or Obs, and end up with the same central values in both cases. I think this would be a cleaner solution. Do you have a feeling if one of the two variants is more stable?
I also think that the sorting could be applied to both variants by just using the central values (floats) to sort the vectors, which might be Obs. If one is interested in the eigenvectors themselves, the correct sorting would be mandatory.

To conclude: I very much like the addition. For smooth integration in the existing workflow, I would really like to have an implementation where the user can work with or without errors on the eigenvectors, just by changing a flag, with the same possibilities in both cases. I would be really interested in your opinion. I think I could also contribute to this, if you like.

fjosw · 2023-06-23T10:48:41Z

Just a GitHub technicality: As a maintainer of the repo @s-kuberski should be able to push to @JanNeuendorf's feature branch in case you want to work on this together.

JanNeuendorf · 2023-06-23T17:25:46Z

Okay, I think there is nothing wrong with having errors as an option in the GEVP method, if explained well.
I agree, that both methods would need to use the same algorithm under the hood. I do not know, what the current scipy solver is doing under the hood; I believe, it is calling some LAPACK routines. Whatever it is doing, I assume it is better than the straightforward implementation used here. We should test this extensively, or give the user the choice which solver to use. There could be a solver option with three possible values: scipy, cholesky and cholesky_errors, where the last two are guaranteed to agree in the mean. I think that just switching the default behavior of the GEVP method is a breaking change. Someone might be using this in a scenario with lot of numerical instability and if the solver now behaves differently, the results could change completely. So the default flag would be the one used now.
Yes, feel very free to contribute!

s-kuberski · 2023-10-20T15:28:03Z

Hi,
I finally implemented the changes that I suggested some time ago:

The error propagation for the eigenvectors can now be chosen via a keyword in GEVP. Everything works as in the old case, including the sorting algorithms.
The solver for the case without error propagation may be chosen (more on this later).
I got a significant speedup (~ factor 7) compared to the earlier version in this pr. If I use the new implementation via Cholesky decomposition without uncertainties, I am about 40% percent slower than with scipy.linalg.eigh, which I consider a success since some python work is done here (compared to the fancy fortran implementation). When switching on the error propagation, it takes O(40) times longer to compute the eigenvectors.

Things to discuss:

I did some checks and the new implementation method='cholesky' seems to agree to machine precision with scipy.linalg.eigh. Therefore, we could abolish the whole idea of choosing the solver and just using the new one, if one wants to propagate uncertainties. Did you, @JanNeuendorf have an example, where the two methods differ?
The sorting can be significantly improved in terms of computer time. I'd open a new pull request, as soon as we have settle this one.

I am happy to discuss and to do changes, if you come up with ideas (also considering the keywords etc.)

fjosw · 2023-10-20T16:19:44Z

pyerrors/misc.py

@@ -183,3 +183,18 @@ def _assert_equal_properties(ol, otype=Obs):
            if hasattr(ol[0], attr):
                if not getattr(ol[0], attr) == getattr(o, attr):
                    raise Exception(f"All Obs in list have to have the same state '{attr}'.")
+
+
+def obsval(o):


I think the built-in float function should be able to do this as well.

float(o) is o.value

Very good point. In general terms, float(o) could not cope with o being None, so obsval would be a bit more flexible. But that should not matter in this case. I can remove the function.

fjosw · 2023-10-27T13:19:46Z

The interface looks good to me, but I didn't run any thorough tests. Would be interesting to hear @JanNeuendorf's opinion.

JanNeuendorf · 2023-11-13T08:58:14Z

I apologize for the late reply.
As for your question. No, I have not found a scenario where the results between this solver and scipy.linalg.eigh differ wildly. But, as you said, the mean values are extremely close but not identical. This means that there could be a scenario, where the ordering of states changes. This is certainly unlikely and I would agree that this might be worth it to have a unified solver.
I like the parameter name vector_obs a lot because it does not sound like it promises anything about the nature of the error propagation. As for your notes on performance, none of them sound too surprising. Is there anything I should do now?

s-kuberski · 2023-11-13T11:35:14Z

Thanks for you feedback.

If you did not see any other difference than machine precision between the two ways of solving the GEVP, I would not worry too much about the sorting of states because it would then also depend on every other step of the analysis chain... Now that I think about this, I would still leave everything as implemented as the scipy solver is still faster in the standard computation. If the user wants to have bit-identical results when switching vector_obs on/off, they can choose to use the same solver. Since the possibility is there, we do not have to remove it.... Is this fine with you?

I'll have another look at the code and push a small update (with obsval removed) that could then hopefully be merged.

s-kuberski · 2023-11-13T14:44:45Z

Sorry for the mess, I'll try to find out why my test fails (randomly or on specific versions).

s-kuberski · 2023-11-13T15:38:46Z

Turns out there is an ambiguity in the overall sign of single eigenvectors that might differ when changing the solver. This could also appear in production but would not have any impact for the physics because it cancels in the projection.
Therefore, it might be good to have the possibility to fix the solver to cholesky, when not using vector_obs.

JanNeuendorf · 2023-11-13T15:49:18Z

Okay, the sign difference might be important. I have a scenario where I project different sides of a matrix with different vectors. The end-results are independent from the sign, but maybe somewhere in between one might plot the correlator without taking an abs. This plot would then be flipped, which would be pretty major. Is there a simple convention that scipy.linalg.eigh uses for the sign, that we can just copy?

s-kuberski · 2023-11-13T17:15:29Z

Just to remind you: The current version of the pull request does not change any behavior of existing code, since you would have to turn on vector_obs or method='cholesky' to use the manual solver. I that sense, if one likes to be insensitive to this in future projects, one could set the method manually.

I don't know if there can be any simple convention as the whole problem is invariant under flipping the sign of an eigenvector. I can think a bit about this, but I would not be sure if there is a way to streamline this (after all, everything is just passed to Fortran code). The routines act on different data in the two cases. Funnily, when I use scipy.linalg.eigh instead of np.linalg.eigh for the ordinary EVP in the Cholesky implementation, I get sign flips for basically every case, whereas, in the current implementation, this only happens in about 1% of the cases.

Added the function error_gevp() to compute the gevp with statistical …

df3836c

…errors.

JanNeuendorf requested a review from fjosw as a code owner June 15, 2023 17:30

fjosw requested a review from s-kuberski June 15, 2023 17:32

fjosw reviewed Jun 16, 2023

View reviewed changes

pyerrors/correlators.py Outdated Show resolved Hide resolved

pyerrors/correlators.py Outdated Show resolved Hide resolved

pyerrors/correlators.py Outdated Show resolved Hide resolved

JanNeuendorf added 6 commits June 16, 2023 12:41

Changed method name from error_gevp to error_GEVP and removed automat…

b97db0f

…ic gamma method

added auto_gamma to error_GEVP

f0bffcb

Specified exceptions in Corr.error_GEVP

6f9cd06

Fixed a wrong path. It should be np.linalg.LinAlgError

ce0794b

Added a test for error_GEVP

cdaaf1b

The tests of error_gevp loads a test matrix

a56f872

fjosw linked an issue Jun 16, 2023 that may be closed by this pull request

GEVP eigenvectors with errors #183

Closed

fjosw approved these changes Jun 16, 2023

View reviewed changes

Incorporated eigenvectors with uncertainties in GEVP routine

9d83e57

fjosw reviewed Oct 20, 2023

View reviewed changes

Cleaned up GEVP routines

08d323a

s-kuberski force-pushed the develop branch from 883f32d to 08d323a Compare November 13, 2023 13:09

s-kuberski added 2 commits November 13, 2023 14:11

Merge branch 'develop' into develop

f2762a6

Cleaned up breaking change from merge

5ac5187

Released tolerance in test of GEVP

7d545df

Repaired broken GEVP test

ed07b87

fjosw self-requested a review November 17, 2023 17:56

fjosw approved these changes Nov 17, 2023

View reviewed changes

fjosw merged commit e1a4d0c into fjosw:develop Nov 17, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Additional GEVP method with errors #195

Additional GEVP method with errors #195

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Additional GEVP method with errors #195

Additional GEVP method with errors #195

Uh oh!

Conversation

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!