pyerrors.mpm

View Source

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

def matrix_pencil_method(corrs, k=1, p=None, **kwargs) View Source

 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    if isinstance(corrs[0], Obs):
27        data = [corrs]
28    else:
29        data = corrs
30
31    lengths = [len(d) for d in data]
32    if lengths.count(lengths[0]) != len(lengths):
33        raise Exception('All datasets have to have the same length.')
34
35    data_sets = len(data)
36    n_data = len(data[0])
37
38    if p is None:
39        p = max(n_data // 2, k)
40    if n_data <= p:
41        raise Exception('The pencil p has to be smaller than the number of data samples.')
42    if p < k or n_data - p < k:
43        raise Exception('Cannot extract', k, 'energy levels with p=', p, 'and N-p=', n_data - p)
44
45    # Construct the hankel matrices
46    matrix = []
47    for n in range(data_sets):
48        matrix.append(scipy.linalg.hankel(data[n][:n_data - p], data[n][n_data - p - 1:]))
49    matrix = np.array(matrix)
50    # Construct y1 and y2
51    y1 = np.concatenate(matrix[:, :, :p])
52    y2 = np.concatenate(matrix[:, :, 1:])
53    # Apply SVD to y2
54    u, s, vh = svd(y2, **kwargs)
55    # Construct z from y1 and SVD of y2, setting all singular values beyond the kth to zero
56    z = np.diag(1. / s[:k]) @ u[:, :k].T @ y1 @ vh.T[:, :k]
57    # Return the sorted logarithms of the real eigenvalues as Obs
58    energy_levels = np.log(np.abs(eig(z, **kwargs)))
59    return sorted(energy_levels, key=lambda x: abs(x.value))

Matrix pencil method to extract k energy levels from data

Implementation of the matrix pencil method based on eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990)

Parameters

data (list): can be a list of Obs for the analysis of a single correlator, or a list of lists of Obs if several correlators are to analyzed at once.
k (int): Number of states to extract (default 1).
p (int): matrix pencil parameter which filters noise. The optimal value is expected between len(data)/3 and 2*len(data)/3. The computation is more expensive the closer p is to len(data)/2 but could possibly suppress more noise (default len(data)//2).