From 8cd96c584e0404480cc69064bfb71f2a57ee746a Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 7 Mar 2022 14:10:02 +0000 Subject: [PATCH] feat: implemented eigenvalue smoothing method of hep-lat/9412087 --- pyerrors/fits.py | 2 +- pyerrors/obs.py | 27 ++++++++++++++++++++++++++- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 6db0b1dc..123356cc 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -461,7 +461,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs): x0 = [0.1] * n_parms if kwargs.get('correlated_fit') is True: - corr = covariance(y, correlation=True) + corr = covariance(y, correlation=True, **kwargs) covdiag = np.diag(1 / np.asarray(dy_f)) condn = np.linalg.cond(corr) if condn > 1e8: diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 2439bc03..5e0d4fa1 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1335,7 +1335,7 @@ def correlate(obs_a, obs_b): return o -def covariance(obs, visualize=False, correlation=False, **kwargs): +def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): r'''Calculates the covariance matrix of a set of observables. The gamma method has to be applied first to all observables. @@ -1348,6 +1348,11 @@ def covariance(obs, visualize=False, correlation=False, **kwargs): If True plots the corresponding normalized correlation matrix (default False). correlation : bool If True the correlation instead of the covariance is returned (default False). + smooth : None or int + If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue + smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the + largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely + small ones. Notes ----- @@ -1372,6 +1377,9 @@ def covariance(obs, visualize=False, correlation=False, **kwargs): corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov))) + if isinstance(smooth, int): + corr = _smooth_eigenvalues(corr, smooth) + errors = [o.dvalue for o in obs] cov = np.diag(errors) @ corr @ np.diag(errors) @@ -1391,6 +1399,23 @@ def covariance(obs, visualize=False, correlation=False, **kwargs): return cov +def _smooth_eigenvalues(corr, E): + """Eigenvalue smoothing as described in hep-lat/9412087 + + corr : np.ndarray + correlation matrix + E : integer + Number of eigenvalues to be left substantially unchanged + """ + if not (2 < E < corr.shape[0] - 1): + raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).") + vals, vec = np.linalg.eigh(corr) + lambda_min = np.mean(vals[:-E]) + vals[vals < lambda_min] = lambda_min + vals /= np.mean(vals) + return vec @ np.diag(vals) @ vec.T + + def _covariance_element(obs1, obs2): """Estimates the covariance of two Obs objects, neglecting autocorrelations."""