refactor: correlators._solve_GEVP simplified and optimized, test added.

pyerrors/correlators.py

@@ -1219,8 +1219,6 @@ 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."""
+    return scipy.linalg.eigh(Gt, G0)[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
 
@@ -290,6 +291,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):