多变量滞后时间序列分量重构技术实现求助

Hey there! Let's work through this multivariate SSA reconstruction problem you're stuck on. The main issue with your current code is that the Reconstruction function is built exclusively for univariate time series—it doesn't account for the multiple channels stacked in your trajectory matrix. Let's break down the key math first, then update the code to handle multivariate data properly.

多变量SSA重构的核心逻辑 #

你的lagged_covariance_matrix函数已经正确生成了多变量轨迹矩阵：对于L个通道，它为每个通道创建滞后矩阵，然后按列堆叠成N × (M*L)的矩阵（其中N = T-M+1，T是时间步数）。

重构的关键步骤：

1. 将轨迹矩阵投影到ST-PCs（排序后的特征向量）后，需要把每个分量的轨迹矩阵拆回L个独立的N × M滞后矩阵（每个通道对应一个）。
2. 对每个通道的滞后矩阵执行对角线平均（和单变量SSA逻辑一致），得到该分量对应的通道子序列。
3. 将所有分量的子序列累加，得到最终的多变量重构序列。

修改后的代码实现 #

以下是适配多变量场景的代码，同时调整了函数间的元数据传递逻辑：

import numpy as np

def lagged_covariance_matrix(data, M):
    """
    Computes the lagged covariance matrix using the Broomhead & King method
    Background: Plaut, G., & Vautard, R. (1994). Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere. Journal of the atmospheric sciences, 51(2), 210-236.
    Arguments:
        data : pxn time series, where p denotes the length of the time series and n the number of channels
        M : window length
    Returns:
        covmat : lagged covariance matrix
        X : multivariate trajectory matrix (N × M*L)
        L : number of channels
    """
    # 为单序列输入添加通道维度
    if np.ndim(data) == 1:
        data = np.reshape(data,(len(data),1))
    T = data.shape[0]
    L = data.shape[1]
    N = T - M + 1
    X = np.zeros((T, L, M))
    for i in range(M):
        X[:,:,i] = np.roll(data, -i, axis = 0)
    X = X[:N]
    # 堆叠为N × M*L的轨迹矩阵
    X = np.reshape(X, (N, M*L), order = 'C')
    # 选择最小的投影基计算协方差矩阵
    if M*L >= N:
        return 1/(M*L) * X.dot(X.T), X, L
    else:
        return 1/N * X.T.dot(X), X, L

def sort_by_eigenvalues(eigenvalues, PCs):
    """
    Sorts the PCs and eigenvalues by descending size of the eigenvalues
    """
    desc = np.argsort(-eigenvalues)
    return eigenvalues[desc], PCs[:,desc]

def multivariate_ssa_reconstruction(M, L, E, X):
    """
    Reconstructs multivariate time series from the lagged trajectory matrix and ST-PCs.
    Arguments:
        M : window length
        L : number of channels (time series)
        E : eigenvector basis (ST-PCs)
        X : multivariate trajectory matrix (N × M*L)
    Returns:
        recons_data : T × L matrix of reconstructed time series (T = N + M -1)
    """
    T = X.shape[0] + M - 1
    recons_data = np.zeros((T, L))
    
    # 遍历每个特征向量分量
    for i in range(E.shape[1]):
        # 提取第i个特征向量并投影轨迹矩阵
        pc = E[:, i].reshape(-1, 1)
        component_trajectory = X @ pc @ pc.T
        
        # 将分量轨迹矩阵拆分为L个通道专属的滞后矩阵
        channel_lagged_matrices = np.split(component_trajectory, L, axis=1)
        
        # 对每个通道的滞后矩阵执行对角线平均
        for channel_idx in range(L):
            lag_mat = channel_lagged_matrices[channel_idx]
            Q = np.flipud(lag_mat)
            # 计算每个时间步的对角线平均值
            for k in range(T):
                offset = -(T - M - k)
                diag = np.diagonal(Q, offset=offset)
                recons_data[k, channel_idx] += diag.mean()
    
    return recons_data

# -----------------------------------------------------------------------------
# 示例使用
# -----------------------------------------------------------------------------
# 模拟输入数据：288个时间步 × 3个通道
np.random.seed(42)
data = np.random.randn(288, 3)
M = 45  # 窗口长度

# 计算协方差矩阵、轨迹矩阵和通道数
covmat, X, L = lagged_covariance_matrix(data, M)
# 获取排序后的特征值和特征向量
vals, vecs = np.linalg.eig(covmat)
eig_vals, eig_vecs = sort_by_eigenvalues(vals, vecs)
# 执行多变量重构
recons_data = multivariate_ssa_reconstruction(M, L, eig_vecs, X)

# 验证输出形状与输入一致（288 × 3）
print(f"Reconstructed data shape: {recons_data.shape}")

关键修改说明 #

1. lagged_covariance_matrix更新：新增返回通道数L，避免在重构函数中重复计算，减少维度对齐错误。
2. 多变量重构逻辑：
   • 用np.split将投影后的分量轨迹矩阵拆回通道专属的滞后矩阵（和原始堆叠逻辑匹配）。
   • 对每个通道独立执行对角线平均，再累加所有分量的结果。
3. 维度一致性：最终输出recons_data的形状与输入data完全一致（T × L），符合你要求的数据格式。

额外实用提示 #

• 实值特征向量：由于协方差矩阵是对称矩阵，建议用np.linalg.eigh替代np.linalg.eig，确保得到实值特征值和正交特征向量，避免数值噪声导致的伪复数结果。
• 蒙特卡洛检验集成：多变量SSA的蒙特卡洛检验可通过生成替代数据（如相位随机化的输入序列），重复分解流程，将真实数据的特征值与替代数据的特征值分布对比，评估分量的显著性。

内容的提问来源于stack exchange，提问作者Daan