最简洁高效的方式是用 np.cov 计算完整协方差矩阵后接 np.triu 提取上三角；若需一维向量，用 np.triu_indices 索引提取；注意数据格式：N×D 输入须设 rowvar=False。

直接用 np.cov 计算完整协方差矩阵，再用布尔索引或 np.triu 提取上三角部分，是最简洁高效的方式。不需要手动循环或分块计算，NumPy 的底层实现已高度优化。

用 np.cov + n

p.triu 一步提取上三角

默认 np.cov(X) 返回对称的 D×D 协方差矩阵（X 是 D 维特征 × N 样本的二维数组，注意 rowvar=True 是默认）。只需再套一层 np.triu 即可保留上三角（含对角线），其余置零：

import numpy as np
假设 X.shape = (D, N)，例如 5 个特征、1000 个样本
X = np.random.randn(5, 1000)
cov_full = np.cov(X)           # D×D 对称矩阵
cov_upper = np.triu(cov_full)  # 上三角：对角线及以上为原值，下方为 0

若只需上三角非零值（一维数组），用 np.triu_indices

当后续只处理上三角元素（如聚类、图构建、存储压缩），可直接索引提取一维向量，节省内存：

i, j = np.triu_indices(cov_full.shape[0])
upper_vec = cov_full[i, j]  # 长度为 D*(D+1)//2 的一维数组
# 例如 D=3 → indices: (0,0),(0,1),(0,2),(1,1),(1,2),(2,2)

注意数据布局和参数设置

np.cov 默认把每行当作一个变量（即 rowvar=True）。如果你的数据是 N×D（样本数×特征数），必须显式设 rowvar=False，否则结果会出错：

X_nd = np.random.randn(1000, 5)  # N×D 格式
cov_nd = np.cov(X_nd, rowvar=False)  # 正确：输出 5×5
# 错误写法（不加 rowvar=False）会得到 1000×1000 矩阵！

无需手动优化：避免常见低效操作

以下做法既慢又易错，不推荐：

• 用双重 for 循环逐个计算 cov(X[i], X[j]) —— 完全放弃 NumPy 向量化优势
• 先计算下三角再转置 —— 多余且破坏对称性假设
• 用 scipy.linalg.eigh 等分解后重构 —— 过度设计，协方差本身已足够稳定