Python中尺度矩阵（协方差矩阵）的高效计算及向量化改写

嘿，这个问题我太熟了！numpy里确实没有专门的「沿指定轴批量计算外积」的函数，但咱们完全可以用广播机制或者numpy.einsum来搞定，全程不用列表和for循环，效率还高得多。

核心思路：利用矩阵乘法或Einstein求和替代循环 #

首先明确需求：假设多变量时间序列是二维数组，比如形状为(n_time_steps, n_features)（每行对应一个时间步的特征向量），要计算时间区间(t0, t1)内的协方差（尺度）矩阵，核心是批量计算所有样本向量的外积并求和，再做归一化。

方法1：矩阵乘法（最简洁高效） #

矩阵乘法的本质就是批量外积的求和！对于截取后的子数组X_slice（形状(m, n)，m是区间内的时间步数，n是特征数），所有样本向量外积的总和直接等于X_slice.T @ X_slice（X_slice.T是(n, m)，相乘后得到(n, n)的矩阵，正好是外积和）。

如果要计算样本协方差矩阵（减去均值，除以m-1），完整代码如下：

import numpy as np

def compute_scale_matrix(X, t0, t1):
    # 截取时间区间内的子序列
    X_slice = X[t0:t1, :]
    # 中心化（减去每个特征的均值）
    X_centered = X_slice - X_slice.mean(axis=0)
    # 计算协方差矩阵：外积和除以样本数-1（样本协方差）
    scale_matrix = (X_centered.T @ X_centered) / (X_centered.shape[0] - 1)
    return scale_matrix

方法2：Einstein求和（更直观，灵活性更高） #

如果你想更清晰地控制维度运算，可以用np.einsum，它的语法直接对应外积求和的逻辑：

def compute_scale_matrix_einsum(X, t0, t1):
    X_slice = X[t0:t1, :]
    X_centered = X_slice - X_slice.mean(axis=0)
    # 'ij,ik->jk' 表示：对每个样本i，将j和k维度的元素相乘，然后对i求和
    outer_sum = np.einsum('ij,ik->jk', X_centered, X_centered)
    scale_matrix = outer_sum / (X_centered.shape[0] - 1)
    return scale_matrix

适配不同的轴方向 #

如果你的时间序列数组是(n_features, n_time_steps)（特征在轴0，时间步在轴1），只需要调整轴参数即可：

# 假设X形状是(n_features, n_time_steps)
X_slice = X[:, t0:t1]
X_centered = X_slice - X_slice.mean(axis=1, keepdims=True)
scale_matrix = (X_centered @ X_centered.T) / (X_centered.shape[1] - 1)

验证效果 #

用测试数据验证两种方法的一致性：

np.random.seed(42)
X = np.random.randn(100, 3)  # 100个时间步，3个特征
scale1 = compute_scale_matrix(X, 10, 50)
scale2 = compute_scale_matrix_einsum(X, 10, 50)
print(np.allclose(scale1, scale2))  # 输出True，说明结果一致

这两种方法都是纯numpy数组操作，完全抛弃了列表和循环，而且底层是C实现的，运行效率比Python循环高几个数量级，完美满足需求！

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