pyerrors.mpm

 0import numpy as np
 1import scipy.linalg
 2from .obs import Obs
 3from .linalg import svd, eig
 4
 5
 6def matrix_pencil_method(corrs, k=1, p=None, **kwargs):
 7    """Matrix pencil method to extract k energy levels from data
 8
 9    Implementation of the matrix pencil method based on
10    eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990)
11
12    Parameters
13    ----------
14    data : list
15        can be a list of Obs for the analysis of a single correlator, or a list of lists
16        of Obs if several correlators are to analyzed at once.
17    k : int
18        Number of states to extract (default 1).
19    p : int
20        matrix pencil parameter which filters noise. The optimal value is expected between
21        len(data)/3 and 2*len(data)/3. The computation is more expensive the closer p is
22        to len(data)/2 but could possibly suppress more noise (default len(data)//2).
23    """
24    if isinstance(corrs[0], Obs):
25        data = [corrs]
26    else:
27        data = corrs
28
29    lengths = [len(d) for d in data]
30    if lengths.count(lengths[0]) != len(lengths):
31        raise Exception('All datasets have to have the same length.')
32
33    data_sets = len(data)
34    n_data = len(data[0])
35
36    if p is None:
37        p = max(n_data // 2, k)
38    if n_data <= p:
39        raise Exception('The pencil p has to be smaller than the number of data samples.')
40    if p < k or n_data - p < k:
41        raise Exception('Cannot extract', k, 'energy levels with p=', p, 'and N-p=', n_data - p)
42
43    # Construct the hankel matrices
44    matrix = []
45    for n in range(data_sets):
46        matrix.append(scipy.linalg.hankel(data[n][:n_data - p], data[n][n_data - p - 1:]))
47    matrix = np.array(matrix)
48    # Construct y1 and y2
49    y1 = np.concatenate(matrix[:, :, :p])
50    y2 = np.concatenate(matrix[:, :, 1:])
51    # Apply SVD to y2
52    u, s, vh = svd(y2, **kwargs)
53    # Construct z from y1 and SVD of y2, setting all singular values beyond the kth to zero
54    z = np.diag(1. / s[:k]) @ u[:, :k].T @ y1 @ vh.T[:, :k]
55    # Return the sorted logarithms of the real eigenvalues as Obs
56    energy_levels = np.log(np.abs(eig(z, **kwargs)))
57    return sorted(energy_levels, key=lambda x: abs(x.value))