Merge pull request #100 from fjosw/feature/GEVP_return_all_vectors

GEVP rehaul 2

examples/06_gevp.ipynb

@@ -168,9 +168,9 @@
     }
    ],
    "source": [
-    "vec = matrix_V1V1.GEVP(t0=3, ts=6 ,state=0, sorted_list=None)\n",
+    "vecs = matrix_V1V1.GEVP(t0=3, ts=6, sort=None)\n",
-    "assert len(vec) == matrix_V1V1.N\n",
+    "assert len(vecs[0]) == matrix_V1V1.N\n",
-    "print(vec)"
+    "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",

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"):
+    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
-        """Solve the general eigenvalue problem on the current correlator
+        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
+            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
-            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+            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.
+           Returns only the vector(s) for a specified state. 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
-        """
 
         if self.N == 1:
             raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
-
+        if ts is not None:
-        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 <= 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)
+            reordered_vecs = _GEVP_solver(Gt, G0)
-            sp_vec = sp_vecs[state]
+
-            return sp_vec
+        elif sort in ["Eigenvalue", "Eigenvector"]:
-        elif sorted_list in ["Eigenvalue", "Eigenvector"]:
+            if sort == "Eigenvalue" and ts is not None:
-            if sorted_list == "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)
+                    all_vecs.append(_GEVP_solver(Gt, G0))
-                    if sorted_list == "Eigenvalue":
-                        sp_vec = sp_vecs[state]
-                        all_vecs.append(sp_vec)
-                    else:
-                        all_vecs.append(sp_vecs)
                 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.
+        All other parameters are identical to the ones of Corr.GEVP.
-             "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
         """
-        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 = []
+        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-        for i in range(Ntrunc):
-            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
 
         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
+def _GEVP_solver(Gt, G0):
-    sp_val, sp_vecs = scipy.linalg.eigh(Gt, G0)
+    """Helper function for solving the GEVP and sorting the eigenvectors.
-    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]
+    The helper function assumes that both provided matrices are symmetric and
-    return sp_vecs
+    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]

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)
+    for (ts, sort) in zip([None, 1, 1], ["Eigenvalue", "Eigenvector", None]):
-    vec_1 = corr_mat.GEVP(0, 1, state=1, sorted_list=None)
+        vecs = corr_mat.GEVP(0, ts=ts, sort=sort)
 
-    corr_0 = corr_mat.projected(vec_0)
+        corr_0 = corr_mat.projected(vecs[0])
-    corr_1 = corr_mat.projected(vec_1)
+        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):