如何用更线性代数的方式通过矩阵计算回归预测标准误？

当然有更贴合线性代数逻辑的实现方式，核心是利用二次型的矩阵表示替代逐元素乘法与求和：

方法1：利用矩阵乘积的对角线元素 #

对于形状为 (p, n) 的协变量矩阵 X（每一列对应一个待预测观测），所有观测的预测方差本质是每个列向量 x_i 的二次型 x_i^T @ sigma @ x_i。而这些二次型恰好是矩阵 X.T @ sigma @ X 的对角线元素。

对应的代码实现：

import numpy as np

# 假设 sigma 是 (p,p) 的协方差矩阵，X 是 (p,n) 的协变量矩阵
pred_se = np.sqrt(np.diag(X.T @ sigma @ X))

这个写法完全基于矩阵乘法的线性代数定义，没有额外的逐元素操作，逻辑直观——X.T @ sigma @ X 是一个 n×n 的矩阵，其对角线就是每个观测的预测方差，取对角线后开根号即得标准误。

方法2：利用爱因斯坦求和约定（更高效） #

如果观测数 n 很大，生成 n×n 的中间矩阵会占用较多内存，此时可以用 np.einsum 直接计算每个观测的二次型，避免生成冗余的非对角线元素：

pred_se = np.sqrt(np.einsum('pi,pq,qj->i', X, sigma, X))

这里的 einsum 表达式 'pi,pq,qj->i' 表示：对每个 i，计算 sum_p sum_q X[p,i] * sigma[p,q] * X[q,i]，这正好是 x_i^T @ sigma @ x_i 的展开形式，直接得到每个观测的预测方差数组，再开根号即可。

和现有方法的对比 #

你当前的方法 np.sqrt((X.T * (sigma@X.T)).sum(0)) 本质是把二次型展开为逐元素相乘再求和，而上面两种方法直接利用线性代数的矩阵/张量运算特性，代码更简洁，逻辑也更贴合回归分析中预测方差的数学定义。

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