diff --git a/examples/06_gevp.ipynb b/examples/06_gevp.ipynb index 9bb9ed49..e969008c 100644 --- a/examples/06_gevp.ipynb +++ b/examples/06_gevp.ipynb @@ -168,9 +168,9 @@ } ], "source": [ - "vec = matrix_V1V1.GEVP(t0=3, ts=6 ,state=0, sorted_list=None)\n", - "assert len(vec) == matrix_V1V1.N\n", - "print(vec)" + "vecs = matrix_V1V1.GEVP(t0=3, ts=6, sort=None)\n", + "assert len(vecs[0]) == matrix_V1V1.N\n", + "print(vecs[0])" ] }, { @@ -202,7 +202,7 @@ } ], "source": [ - "matrix_V1V1.projected(vec).m_eff().show(comp=single_smearing.m_eff(), auto_gamma=True)" + "matrix_V1V1.projected(vecs[0]).m_eff().show(comp=single_smearing.m_eff(), auto_gamma=True)" ] }, { @@ -266,7 +266,7 @@ "\n", "# We then solve the GEVP to get eigenvectors corresponding to the ground state.\n", "\n", - "vec_ground = matrix_VnVn.GEVP(t0=3, ts=6, state=0, sorted_list=None)\n", + "vec_ground = matrix_VnVn.GEVP(t0=3, ts=6, sort=None)[0]\n", "\n", "# Now we project the matrix-correlators to get new correlators belonging to the ground state.\n", "\n", diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 233c7499..ab4f6f50 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -241,34 +241,49 @@ class Corr: if self.N == 1: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") - def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"): - """Solve the general eigenvalue problem on the current correlator + def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs): + r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors. + + The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the + largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing + ```python + C.GEVP(t0=2)[0] # Ground state vector(s) + C.GEVP(t0=2)[:3] # Vectors for the lowest three states + ``` Parameters ---------- t0 : int - The time t0 for G(t)v= lambda G(t_0)v + The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$ ts : int - fixed time G(t_s)v= lambda G(t_0)v if return_list=False - If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. + fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None. + If sort="Eigenvector" it gives a reference point for the sorting method. + sort : string + If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned. + - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. + - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state. + The reference state is identified by its eigenvalue at $t=t_s$. + + Other Parameters + ---------------- state : int - The state one is interested in, ordered by energy. The lowest state is zero. - sorted_list : string - if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. - "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. - "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. - The reference state is identified by its eigenvalue at t=ts - """ + Returns only the vector(s) for a specified state. The lowest state is zero. + ''' if self.N == 1: raise Exception("GEVP methods only works on correlator matrices and not single correlators.") - - symmetric_corr = self.matrix_symmetric() - if sorted_list is None: - if (ts is None): - raise Exception("ts is required if sorted_list=None.") + if ts is not None: if (ts <= t0): raise Exception("ts has to be larger than t0.") + + if "sorted_list" in kwargs: + warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning) + sort = kwargs.get("sorted_list") + + symmetric_corr = self.matrix_symmetric() + if sort is None: + if (ts is None): + raise Exception("ts is required if sort=None.") 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") @@ -277,11 +292,10 @@ class Corr: G0[i, j] = symmetric_corr[t0][i, j].value Gt[i, j] = symmetric_corr[ts][i, j].value - sp_vecs = _GEVP_solver(Gt, G0) - sp_vec = sp_vecs[state] - return sp_vec - elif sorted_list in ["Eigenvalue", "Eigenvector"]: - if sorted_list == "Eigenvalue" and ts is not None: + reordered_vecs = _GEVP_solver(Gt, G0) + + elif sort in ["Eigenvalue", "Eigenvector"]: + if sort == "Eigenvalue" and ts is not None: warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) all_vecs = [None] * (t0 + 1) for t in range(t0 + 1, self.T): @@ -292,43 +306,34 @@ class Corr: G0[i, j] = symmetric_corr[t0][i, j].value Gt[i, j] = symmetric_corr[t][i, j].value - sp_vecs = _GEVP_solver(Gt, G0) - if sorted_list == "Eigenvalue": - sp_vec = sp_vecs[state] - all_vecs.append(sp_vec) - else: - all_vecs.append(sp_vecs) + all_vecs.append(_GEVP_solver(Gt, G0)) except Exception: all_vecs.append(None) - if sorted_list == "Eigenvector": + if sort == "Eigenvector": if (ts is None): raise Exception("ts is required for the Eigenvector sorting method.") all_vecs = _sort_vectors(all_vecs, ts) - all_vecs = [a[state] if a is not None else None for a in all_vecs] + + reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] else: - raise Exception("Unkown value for 'sorted_list'.") + raise Exception("Unkown value for 'sort'.") - return all_vecs + if "state" in kwargs: + return reordered_vecs[kwargs.get("state")] + else: + return reordered_vecs - def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): + def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): """Determines the eigenvalue of the GEVP by solving and projecting the correlator Parameters ---------- - t0 : int - The time t0 for G(t)v= lambda G(t_0)v - ts : int - fixed time G(t_s)v= lambda G(t_0)v if return_list=False - If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. state : int The state one is interested in ordered by energy. The lowest state is zero. - sorted_list : string - if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. - "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. - "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. - The reference state is identified by its eigenvalue at t=ts + + All other parameters are identical to the ones of Corr.GEVP. """ - vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) + vec = self.GEVP(t0, ts=ts, sort=sort)[state] return self.projected(vec) def Hankel(self, N, periodic=False): @@ -1176,9 +1181,7 @@ class Corr: if basematrix.N != self.N: raise Exception('basematrix and targetmatrix have to be of the same size.') - evecs = [] - for i in range(Ntrunc): - evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None)) + evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) rmat = [] @@ -1219,8 +1222,10 @@ def _sort_vectors(vec_set, ts): return sorted_vec_set -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.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 +def _GEVP_solver(Gt, G0): + """Helper function for solving the GEVP and sorting the eigenvectors. + + The helper function assumes that both provided matrices are symmetric and + only processes the lower triangular part of both matrices. In case the matrices + are not symmetric the upper triangular parts are effectively discarded.""" + return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 5242339f..97bbebb1 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -1,5 +1,6 @@ import os import numpy as np +import scipy import pyerrors as pe import pytest @@ -228,28 +229,70 @@ def test_matrix_corr(): corr_ab = 0.5 * corr_aa corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.gamma_method() corr_mat.item(0, 0) - vec_0 = corr_mat.GEVP(0, 1, sorted_list=None) - vec_1 = corr_mat.GEVP(0, 1, state=1, sorted_list=None) + for (ts, sort) in zip([None, 1, 1], ["Eigenvalue", "Eigenvector", None]): + vecs = corr_mat.GEVP(0, ts=ts, sort=sort) - corr_0 = corr_mat.projected(vec_0) - corr_1 = corr_mat.projected(vec_1) + corr_0 = corr_mat.projected(vecs[0]) + corr_1 = corr_mat.projected(vecs[1]) assert np.all([o == 0 for o in corr_0 - corr_aa]) assert np.all([o == 0 for o in corr_1 - corr_aa]) - corr_mat.GEVP(0, sorted_list="Eigenvalue") - corr_mat.GEVP(0, 1, sorted_list="Eigenvector") - corr_mat.matrix_symmetric() + corr_mat.GEVP(0, state=0) + corr_mat.Eigenvalue(2, state=0) + + +def test_GEVP_warnings(): + corr_aa = _gen_corr(1) + corr_ab = 0.5 * corr_aa + + corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.item(0, 0) with pytest.warns(RuntimeWarning): - corr_mat.GEVP(0, 1, sorted_list="Eigenvalue") + corr_mat.GEVP(0, 1, sort="Eigenvalue") + + with pytest.warns(DeprecationWarning): + corr_mat.GEVP(0, sorted_list="Eigenvalue") + +def test_GEVP_exceptions(): + corr_aa = _gen_corr(1) + corr_ab = 0.5 * corr_aa + + corr_mat = pe.Corr(np.array([[corr_aa, corr_ab], [corr_ab, corr_aa]])) + corr_mat.item(0, 0) + + with pytest.raises(Exception): + corr_mat.item(0, 0).projected() + + with pytest.raises(Exception): + corr_mat.item(0, 0).GEVP(2) + + with pytest.raises(Exception): + corr_mat.item(0, 0).matrix_symmetric() + + with pytest.raises(Exception): + corr_mat.GEVP(0, 0, sort=None) + + with pytest.raises(Exception): + corr_mat.GEVP(0, sort=None) + + with pytest.raises(Exception): + corr_mat.GEVP(1, 0, sort="Eigenvector") + + with pytest.raises(Exception): + corr_mat.GEVP(0, 1, sort="This sorting method does not exist.") with pytest.raises(Exception): corr_mat.plottable() + with pytest.raises(Exception): + corr_mat.spaghetti_plot() + with pytest.raises(Exception): corr_mat.show() @@ -257,7 +300,7 @@ def test_matrix_corr(): corr_mat.m_eff() with pytest.raises(Exception): - corr_mat.Hankel() + corr_mat.Hankel(2) with pytest.raises(Exception): corr_mat.plateau() @@ -290,6 +333,20 @@ def test_matrix_symmetric(): assert np.all([np.all(o == o.T) for o in sym_corr_mat]) +def test_GEVP_solver(): + + mat1 = np.random.rand(15, 15) + mat2 = np.random.rand(15, 15) + mat1 = mat1 @ mat1.T + mat2 = mat2 @ mat2.T + + sp_val, sp_vecs = scipy.linalg.eigh(mat1, mat2) + 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 @ mat2 @ v)) for v in sp_vecs] + + assert np.allclose(sp_vecs, pe.correlators._GEVP_solver(mat1, mat2), atol=1e-14) + + def test_hankel(): corr_content = [] for t in range(8):