From 6b7aa29cdbc3ca350371e0e5f4672a92a83255f0 Mon Sep 17 00:00:00 2001 From: JanNeuendorf Date: Mon, 21 Feb 2022 18:24:40 +0100 Subject: [PATCH 1/4] GEVP method now uses `scipy.eigh()` and applies `matrix_symmetric()`. `Eigenvalue` takes the same arguments as `GEVP` and projects automatically. --- pyerrors/correlators.py | 29 +++++++---------------------- 1 file changed, 7 insertions(+), 22 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 2274f7e1..c1ff9417 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -7,7 +7,7 @@ import scipy.linalg from .obs import Obs, reweight, correlate, CObs from .misc import dump_object, _assert_equal_properties from .fits import least_squares -from .linalg import eigh, inv, cholesky +# from .linalg import eigh, inv, cholesky from .roots import find_root @@ -268,8 +268,8 @@ class Corr: G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") for i in range(self.N): for j in range(self.N): - G0[i, j] = self.content[t0][i, j].value - Gt[i, j] = self.content[ts][i, j].value + G0[i, j] = self.matrix_symmetric().content[t0][i, j].value + Gt[i, j] = self.matrix_symmetric().content[ts][i, j].value sp_vecs = _GEVP_solver(Gt, G0) sp_vec = sp_vecs[state] @@ -301,24 +301,9 @@ class Corr: return all_vecs - def Eigenvalue(self, t0, state=1): - G = self.matrix_symmetric() - G0 = G.content[t0] - L = cholesky(G0) - Li = inv(L) - LT = L.T - LTi = inv(LT) - newcontent = [] - for t in range(self.T): - if self.content[t] is None: - newcontent.append(None) - else: - Gt = G.content[t] - M = Li @ Gt @ LTi - eigenvalues = eigh(M)[0] - eigenvalue = eigenvalues[-state] - newcontent.append(eigenvalue) - return Corr(newcontent) + def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): + vec = self.GEVP(self, t0=t0, ts=ts, state=state, sorted_list=sorted_list) + return self.projected(vec) def Hankel(self, N, periodic=False): """Constructs an NxN Hankel matrix @@ -1081,7 +1066,7 @@ def _sort_vectors(vec_set, ts): def _GEVP_solver(Gt, G0): # Just so normalization an sorting does not need to be repeated. Here we could later put in some checks - sp_val, sp_vecs = scipy.linalg.eig(Gt, G0) + sp_val, sp_vecs = scipy.linalg.eigh(Gt, G0) sp_vecs = [sp_vecs[:, np.argsort(sp_val)[-i]] for i in range(1, sp_vecs.shape[0] + 1)] sp_vecs = [v / np.sqrt((v.T @ G0 @ v)) for v in sp_vecs] return sp_vecs From ab960ee7ad3c110acb50f1b36a073bd41b37f4ff Mon Sep 17 00:00:00 2001 From: JanNeuendorf Date: Mon, 21 Feb 2022 18:34:19 +0100 Subject: [PATCH 2/4] typo --- pyerrors/correlators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index c1ff9417..187ba9b1 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -302,7 +302,7 @@ class Corr: return all_vecs def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): - vec = self.GEVP(self, t0=t0, ts=ts, state=state, sorted_list=sorted_list) + vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) return self.projected(vec) def Hankel(self, N, periodic=False): From c2bc1c0bed1ad70b77e79e98a793341f8adbe93f Mon Sep 17 00:00:00 2001 From: JanNeuendorf Date: Mon, 21 Feb 2022 18:48:20 +0100 Subject: [PATCH 3/4] called matrix_symmetric only once --- pyerrors/correlators.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 187ba9b1..c7078fa3 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -266,10 +266,11 @@ class Corr: if (self.content[t0] is None) or (self.content[ts] is None): raise Exception("Corr not defined at t0/ts") G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") + symmetric_corr = self.matrix_symmetric() for i in range(self.N): for j in range(self.N): - G0[i, j] = self.matrix_symmetric().content[t0][i, j].value - Gt[i, j] = self.matrix_symmetric().content[ts][i, j].value + G0[i, j] = symmetric_corr.content[t0][i, j].value + Gt[i, j] = symmetric_corr[ts][i, j].value sp_vecs = _GEVP_solver(Gt, G0) sp_vec = sp_vecs[state] From 280dbb5d984a2fb3a9f6c91f9a77bc84f351b977 Mon Sep 17 00:00:00 2001 From: JanNeuendorf Date: Tue, 22 Feb 2022 08:25:45 +0100 Subject: [PATCH 4/4] comment deleted --- pyerrors/correlators.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index c7078fa3..1da8f3bc 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -7,7 +7,6 @@ import scipy.linalg from .obs import Obs, reweight, correlate, CObs from .misc import dump_object, _assert_equal_properties from .fits import least_squares -# from .linalg import eigh, inv, cholesky from .roots import find_root