Additional GEVP method with errors (#195)

* Added the function error_gevp() to compute the gevp with statistical errors.

* Changed method name from error_gevp to error_GEVP and removed automatic gamma method

* added auto_gamma to error_GEVP

* Specified exceptions in Corr.error_GEVP

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

* Added a test for error_GEVP

* The tests of error_gevp loads a test matrix

* Incorporated eigenvectors with uncertainties in GEVP routine

* Cleaned up GEVP routines

* Cleaned up breaking change from merge

* Released tolerance in test of GEVP

* Repaired broken GEVP test

---------

Co-authored-by: Simon Kuberski <simon.kuberski@uni-muenster.de>

Author: JanNeuendorf

pyerrors/correlators.py

@@ -8,6 +8,7 @@
 from .misc import dump_object, _assert_equal_properties
 from .fits import least_squares
 from .roots import find_root
+from . import linalg
 
 
 class Corr:
@@ -298,7 +299,7 @@ def matrix_symmetric(self):
             transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
             return 0.5 * (Corr(transposed) + self)
 
-    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
+    def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **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
@@ -317,14 +318,21 @@ def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
             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.
+            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. (default)
             - "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$.
+            - None: The GEVP is solved only at ts, no sorting is necessary
+        vector_obs : bool
+            If True, uncertainties are propagated in the eigenvector computation (default False).
 
         Other Parameters
         ----------------
         state : int
            Returns only the vector(s) for a specified state. The lowest state is zero.
+        method : str
+           Method used to solve the GEVP.
+           - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
+           - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.
         '''
 
         if self.N == 1:
@@ -342,25 +350,43 @@ def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
         else:
             symmetric_corr = self.matrix_symmetric()
 
-        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
-        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
+        def _get_mat_at_t(t, vector_obs=vector_obs):
+            if vector_obs:
+                return symmetric_corr[t]
+            else:
+                return np.vectorize(lambda x: x.value)(symmetric_corr[t])
+        G0 = _get_mat_at_t(t0)
+
+        method = kwargs.get('method', 'eigh')
+        if vector_obs:
+            chol = linalg.cholesky(G0)
+            chol_inv = linalg.inv(chol)
+            method = 'cholesky'
+        else:
+            chol = np.linalg.cholesky(_get_mat_at_t(t0, vector_obs=False))  # Check if matrix G0 is positive-semidefinite.
+            if method == 'cholesky':
+                chol_inv = np.linalg.inv(chol)
+            else:
+                chol_inv = None
 
         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.")
-            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
-            reordered_vecs = _GEVP_solver(Gt, G0)
+            Gt = _get_mat_at_t(ts)
+            reordered_vecs = _GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv)
+            if kwargs.get('auto_gamma', False) and vector_obs:
+                [[o.gm() for o in ev if isinstance(o, Obs)] for ev in reordered_vecs]
 
         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):
                 try:
-                    Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t])
-                    all_vecs.append(_GEVP_solver(Gt, G0))
+                    Gt = _get_mat_at_t(t)
+                    all_vecs.append(_GEVP_solver(Gt, G0, method=method, chol_inv=chol_inv))
                 except Exception:
                     all_vecs.append(None)
             if sort == "Eigenvector":
@@ -369,15 +395,17 @@ def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
                 all_vecs = _sort_vectors(all_vecs, ts)
 
             reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
+            if kwargs.get('auto_gamma', False) and vector_obs:
+                [[[o.gm() for o in evn] for evn in ev if evn is not None] for ev in reordered_vecs]
         else:
-            raise Exception("Unkown value for 'sort'.")
+            raise Exception("Unknown value for 'sort'. Choose 'Eigenvalue', 'Eigenvector' or None.")
 
         if "state" in kwargs:
             return reordered_vecs[kwargs.get("state")]
         else:
             return reordered_vecs
 
-    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
+    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue", **kwargs):
         """Determines the eigenvalue of the GEVP by solving and projecting the correlator
 
         Parameters
@@ -387,7 +415,7 @@ def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
 
         All other parameters are identical to the ones of Corr.GEVP.
         """
-        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
+        vec = self.GEVP(t0, ts=ts, sort=sort, **kwargs)[state]
         return self.projected(vec)
 
     def Hankel(self, N, periodic=False):
@@ -1386,8 +1414,13 @@ def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
         return Corr(newcontent)
 
 
-def _sort_vectors(vec_set, ts):
+def _sort_vectors(vec_set_in, ts):
     """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
+
+    if isinstance(vec_set_in[ts][0][0], Obs):
+        vec_set = [anp.vectorize(lambda x: float(x))(vi) if vi is not None else vi for vi in vec_set_in]
+    else:
+        vec_set = vec_set_in
     reference_sorting = np.array(vec_set[ts])
     N = reference_sorting.shape[0]
     sorted_vec_set = []
@@ -1406,9 +1439,9 @@ def _sort_vectors(vec_set, ts):
                 if current_score > best_score:
                     best_score = current_score
                     best_perm = perm
-            sorted_vec_set.append([vec_set[t][k] for k in best_perm])
+            sorted_vec_set.append([vec_set_in[t][k] for k in best_perm])
         else:
-            sorted_vec_set.append(vec_set[t])
+            sorted_vec_set.append(vec_set_in[t])
 
     return sorted_vec_set
 
@@ -1418,10 +1451,63 @@ def _check_for_none(corr, entry):
     return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
 
 
-def _GEVP_solver(Gt, G0):
-    """Helper function for solving the GEVP and sorting the eigenvectors.
+def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
+    r"""Helper function for solving the GEVP and sorting the eigenvectors.
+
+    Solves $G(t)v_i=\lambda_i G(t_0)v_i$ and returns the eigenvectors v_i
 
     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]
+    are not symmetric the upper triangular parts are effectively discarded.
+
+    Parameters
+    ----------
+    Gt : array
+        The correlator at time t for the left hand side of the GEVP
+    G0 : array
+        The correlator at time t0 for the right hand side of the GEVP
+    Method used to solve the GEVP.
+       - "eigh": Use scipy.linalg.eigh to solve the GEVP.
+       - "cholesky": Use manually implemented solution via the Cholesky decomposition.
+    chol_inv : array, optional
+        Inverse of the Cholesky decomposition of G0. May be provided to
+        speed up the computation in the case of method=='cholesky'
+
+    """
+    if isinstance(G0[0][0], Obs):
+        vector_obs = True
+    else:
+        vector_obs = False
+
+    if method == 'cholesky':
+        if vector_obs:
+            cholesky = linalg.cholesky
+            inv = linalg.inv
+            eigv = linalg.eigv
+            matmul = linalg.matmul
+        else:
+            cholesky = np.linalg.cholesky
+            inv = np.linalg.inv
+
+            def eigv(x, **kwargs):
+                return np.linalg.eigh(x)[1]
+
+            def matmul(*operands):
+                return np.linalg.multi_dot(operands)
+        N = Gt.shape[0]
+        output = [[] for j in range(N)]
+        if chol_inv is None:
+            chol = cholesky(G0)  # This will automatically report if the matrix is not pos-def
+            chol_inv = inv(chol)
+
+        try:
+            new_matrix = matmul(chol_inv, Gt, chol_inv.T)
+            ev = eigv(new_matrix)
+            ev = matmul(chol_inv.T, ev)
+            output = np.flip(ev, axis=1).T
+        except (np.linalg.LinAlgError, TypeError, ValueError):  # The above code can fail because of linalg-errors or because the entry of the corr is None
+            for s in range(N):
+                output[s] = None
+        return output
+    elif method == 'eigh':
+        return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1]

pyerrors/linalg.py

@@ -271,6 +271,12 @@ def eig(obs, **kwargs):
     return w
 
 
+def eigv(obs, **kwargs):
+    """Computes the eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh."""
+    v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs)
+    return v
+
+
 def pinv(obs, **kwargs):
     """Computes the Moore-Penrose pseudoinverse of a matrix of Obs."""
     return derived_observable(lambda x, **kwargs: anp.linalg.pinv(x), obs)

tests/correlators_test.py

@@ -200,17 +200,17 @@ def test_padded_correlator():
 
 def test_corr_exceptions():
     obs_a = pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test'])
-    obs_b= pe.Obs([np.random.normal(0.1, 0.1, 99)], ['test'])
+    obs_b = pe.Obs([np.random.normal(0.1, 0.1, 99)], ['test'])
     with pytest.raises(Exception):
         pe.Corr([obs_a, obs_b])
 
     obs_a = pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test'])
-    obs_b= pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test'], idl=[range(1, 200, 2)])
+    obs_b = pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test'], idl=[range(1, 200, 2)])
     with pytest.raises(Exception):
         pe.Corr([obs_a, obs_b])
 
     obs_a = pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test'])
-    obs_b= pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test2'])
+    obs_b = pe.Obs([np.random.normal(0.1, 0.1, 100)], ['test2'])
     with pytest.raises(Exception):
         pe.Corr([obs_a, obs_b])
 
@@ -436,6 +436,7 @@ def test_GEVP_solver():
     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)
+    assert np.allclose(np.abs(sp_vecs), np.abs(pe.correlators._GEVP_solver(mat1, mat2, method='cholesky')))
 
 
 def test_GEVP_none_entries():
@@ -552,7 +553,7 @@ def test_corr_no_filtering():
     li = [-pe.pseudo_Obs(.2, .1, 'a', samples=10) for i in range(96)]
     for i in range(len(li)):
         li[i].idl['a'] = range(1, 21, 2)
-    c= pe.Corr(li)
+    c = pe.Corr(li)
     b = pe.pseudo_Obs(1, 1e-11, 'a', samples=30)
     c *= b
     assert np.all([c[0].idl == o.idl for o in c])
@@ -572,6 +573,28 @@ def test_corr_symmetric():
         assert scorr[0] == corr[0]
 
 
+def test_error_GEVP():
+    corr = pe.input.json.load_json("tests/data/test_matrix_corr.json.gz")
+    t0, ts, state = 3, 6, 0
+    vec_regular = corr.GEVP(t0=t0)[state][ts]
+    vec_errors = corr.GEVP(t0=t0, vector_obs=True, auto_gamma=True)[state][ts]
+    vec_regular_chol = corr.GEVP(t0=t0, method='cholesky')[state][ts]
+    print(vec_errors)
+    print(type(vec_errors[0]))
+    assert(np.isclose(np.asarray([e.value for e in vec_errors]), vec_regular).all())
+    assert(all([e.dvalue > 0. for e in vec_errors]))
+
+    projected_regular = corr.projected(vec_regular).content[ts + 1][0]
+    projected_errors = corr.projected(vec_errors).content[ts + 1][0]
+    projected_regular.gamma_method()
+    projected_errors.gamma_method()
+    assert(projected_errors.dvalue > projected_regular.dvalue)
+    assert(corr.GEVP(t0, vector_obs=True)[state][42] is None)
+
+    assert(np.isclose(vec_regular_chol, vec_regular).all())
+    assert(np.isclose(corr.GEVP(t0=t0, state=state)[ts], vec_regular).all())
+
+
 def test_corr_array_ndim1_init():
     y = [pe.pseudo_Obs(2 + np.random.normal(0.0, 0.1), .1, 't') for i in np.arange(5)]
     cc1 = pe.Corr(y)
@@ -688,7 +711,6 @@ def test_matrix_trace():
     for el in corr.trace():
         assert el == 0
 
-
     with pytest.raises(ValueError):
         corr.item(0, 0).trace()
 

tests/data/test_matrix_corr.json.gz

@@ -297,6 +297,10 @@ def test_matrix_functions():
         for j in range(dim):
             assert tmp[j].is_zero()
 
+    # Check eigv
+    v2 = pe.linalg.eigv(sym)
+    assert(np.all(v - v2).is_zero())
+
     # Check eig function
     e2 = pe.linalg.eig(sym)
     assert np.all(np.sort(e) == np.sort(e2))