Scikit-learn与Scipy余弦距离计算结果不一致问题排查

处理N×M的观测-特征矩阵时，计算观测间两两距离得到N×N距离矩阵是常见需求。Scipy和Scikit-learn都提供了相关实现，理论上同度量下结果应该一致，我编写了如下测试函数验证：

import numpy as np
from sklearn.metrics import pairwise_distances
from scipy.spatial.distance import pdist, squareform

def test_equivalence(arr: np.array, metric="cosine") -> bool:
    scipy_result = squareform(pdist(arr, metric=metric))
    sklearn_result = pairwise_distances(arr, metric=metric)
    return np.isclose(scipy_result, sklearn_result).all()

我有一个1219×37652的归一化矩阵（每行和为1）：

• 完整矩阵测试时，两个库的结果一致；
• 截取矩阵最后i列测试时，仅当i≤25676时结果一致，超过该维度阈值后结果不符；
• 测试其他度量时，除correlation外，其余度量的结果均一致。

简化测试矩阵（1219×96）的加载方式为np.load("tf_matrix.npz")["arr_0"]。

现需明确：

1. 该现象的原因是什么？
2. 应该选用哪种实现？

原因分析与实现选择建议 #

原因分析 #

• 数值计算路径与误差累积差异：Scipy和Scikit-learn的余弦距离底层实现细节不同。Scipy的pdist直接基于公式1 - (u·v)/(||u|| ||v||)计算，高维度下浮点数累加的舍入误差会被逐步放大；而Scikit-learn的pairwise_distances采用了分块计算、归一化顺序调整等优化策略，误差累积的方向和程度与Scipy不同，当维度超过某个阈值时，误差累积的差异达到了np.isclose的检测阈值，导致结果不一致。
• 归一化矩阵的L2范数计算误差：你的矩阵是L1归一化（每行和为1），但余弦距离依赖L2范数。高维度下，L2范数的计算涉及大量平方和累加，两种实现的浮点运算顺序不同，进一步放大了误差差异。
• correlation度量的固有差异：correlation距离的计算涉及均值、协方差的计算，两个库的实现逻辑差异更大，因此即使低维度也容易出现不一致，本质是同一类数值计算路径差异导致的问题。

实现选择建议 #

• 机器学习场景优先选Scikit-learn：如果你的任务属于机器学习生态（如聚类、分类前的距离计算），pairwise_distances与Scikit-learn的其他组件兼容性更好，且针对大规模数据做了内存优化（支持稀疏矩阵、分块处理），更适合高维度场景。
• Scipy生态联动选Scipy：如果你的工作流以Scipy为核心，或需要与Scipy的其他距离分析工具配合，可选择pdist+squareform。但建议使用np.float64精度计算，缓解高维度下的误差累积问题。
• 根据精度需求决策：若仅需距离的相对排序（如聚类），两种实现的差异不会影响结果；若要求高精度绝对距离值，建议对比两者的误差范围，选择符合需求的实现，或自行实现严格的余弦距离计算逻辑。

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