向量自相关加权和高效计算方法

本文介绍一种基于 NumPy 矩阵运算的向量化方法，用于高效计算给定二维数组各行与全数组点积之和，并按位移步长 $ t $ 进行归一化，避免 Python 显式循环，显著提升性能。

本文介绍一种基于 NumPy 矩阵运算的向量化方法，用于高效计算给定二维数组各行与全数组点积之和，并按位移步长 $ t $ 进行归一化，避免 Python 显式循环，显著提升性能。

在时间序列分析、信号处理或统计建模中，常需对数据矩阵按“滞后”（lag）维度执行类似自相关结构的聚合计算：对每个起始行索引 $ t $，取第 $ t $ 行作为权重向量，与整个数据矩阵做逐行点积，再对所有点积结果求和，并除以有效样本数 $ n - t $。若用传统 for 循环实现，不仅代码冗长，且在大数据量下效率低下。

幸运的是，该操作可完全向量化。核心洞察在于：

• 第 $ t $ 行与全矩阵的点积等价于矩阵乘法 my_data @ my_data[t, :]；
• 将所有 $ t $ 对应的点积结果统一表达，等价于计算完整 Gram 矩阵 my_data @ my_data.T（形状为 $ n \times n $），其中第 $ t $ 列即为 my_data @ my_data[t, :] 的结果；
• 因此，对 Gram 矩阵按列求和（sum(axis=0)），再逐元素除以对应分母 $ n - t $，即可一次性获得全部结果。

以下是简洁、高效的实现：

import numpy as np

my_data = np.array([
    [1, 1, 1],
    [2, 2, 2],
    [3, 3, 3],
    [4, 4, 4],
    [5, 5, 5]
])
n = len(my_data)

# 向量化计算：t ∈ [0, 1, ..., n-1]
t_values = np.arange(n)
gram_matrix = my_data @ my_data.T      # shape: (n, n)
result = gram_matrix.sum(axis=0) / (n - t_values)

print(result)
# 输出: [  9.   22.5  45.   90.  225. ]

✅ 关键优势：

• 零显式循环，充分利用 NumPy 底层优化（BLAS 加速）；
• 时间复杂度仍为 $ O(n^2 d) $（$ d $ 为列数），但常数因子远低于解释型循环；
• 代码简短、可读性强，易于嵌入 pipeline 或批量处理。

⚠️ 注意事项：

• 此方法要求 my_data 为二维数组（ndim == 2），若输入为一维需先 reshape(-1, 1)；
• 分母 n - t_values 在 $ t = n $ 时为零，故 t 范围必须严格限制在 [0, n-1]，代码中使用 np.arange(n) 已天然保证安全；
• 内存占用为 $ O(n^2) $，当 $ n $ 超过 $ 10^4 $ 量级时，需评估 Gram 矩阵是否适合内存——此时可考虑分块计算或迭代近似策略。

综上，该向量化方案是兼顾简洁性、性能与可维护性的最佳实践，适用于大多数中等规模科学计算场景。