From ba054fa11caed05f8b12e94368b6987492c113b9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 16 May 2022 11:35:23 +0100 Subject: [PATCH] refactor: correlators._solve_GEVP simplified and optimized, test added. --- pyerrors/correlators.py | 8 +++----- tests/correlators_test.py | 15 +++++++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 233c7499..19d0c334 100644 --- a/pyerrors/correlators.py +++ b/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] diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 5242339f..c520917a 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 @@ -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):