pyerrors.mpm

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