提到对称正定矩阵上的距离计算很多人第一反应就是直接用矩阵差的Frobenius范数。我在一个给协方差矩阵做分类的项目里也这么干过直到发现用欧氏距离插值出来的均值矩阵连正定性都保不住才意识到“距离”这个词在SPD流形上远不是套个范数公式那么简单。这篇文章要聊的是SPD流形上最常用的四种距离度量欧几里得距离、仿射不变黎曼距离AIRM、对数欧几里得距离LED和Log-Cholesky距离。它们各自有不同的几何含义、计算成本和适用范围直接决定了均值、插值、聚类这些下游操作的结果。适合正在做流形学习、脑电信号分类、图像区域描述子匹配或者刚接触协方差矩阵几何性质的读者参考。我会把每种度的定义、直觉、代码实现和实际踩过的坑一起说清楚。1. 为什么“在SPD流形上算距离”不是一件小事1.1 正定约束制造了一个弯曲空间先明确这里讨论的对象。一个对称正定矩阵P要求它对称、且对任意非零向量x都有 x^T P x 0。这个约束看起来只是代数上的要求但它把合法矩阵集合限制成了一个弯曲的流形而不是一个平直的向量空间。举个直觉例子所有2×2的对称矩阵可以想象成一个3维平直空间三个独立元素但其中只有同时满足行列式大于0、左上角元素大于0的那部分才是SPD矩阵。这个区域像圆锥内部边界是半正定矩阵。在这个区域内走直线很可能一步跨出边界。所以当你计算“两个协方差矩阵相差多少”时不能简单把它们当向量空间里的点来量。流形的弯曲几何决定了“直线”不是我们熟悉的直线“距离”也不是我们熟悉的距离。1.2 欧氏线性操作会把正定性搞丢我在项目里踩的第一个坑就是均值外插。假设手头有P和Q两个SPD矩阵用普通算术平均得到 (PQ)/2这个结果恰好还是对称正定的问题不明显。但一旦做线性插值 (1-t)P tQt在0到1之间还好t一旦大于1或者小于0正定性立刻无法保证。更严重的场景是计算测地线或者做流形上的均值。欧氏均值只考虑平直空间的平均完全不理会流形曲率结果常常外插到非正定区域。这就是为什么后来我彻底不用欧氏方式在SPD矩阵集上做平均和距离计算。1.3 距离度量是下游一切操作的地基选择哪种距离会连锁影响聚类、分类、插值、均值甚至梯度计算。同一个数据集用欧氏距离和用仿射不变黎曼距离算出来的聚类中心完全不同同一条“插值路径”欧氏路径中间点可能退化而黎曼测地线路径上的每一个中间点都保持正定。所以“距离计算”这个问题本质上是给整个流形几何定调子。理解了这一点下面四种度量就好理解了它们不是数学游戏而是给SPD流形提供的四种不同“尺子”。2. 四种常用度量各自的定义与几何直觉2.1 欧几里得距离最直接但最“不流形”的尺子公式很简单d_F(P, Q) ||P - Q||_F sqrt( Σ_ij (P_ij - Q_ij)^2 )这里我用Frobenius范数也就是矩阵每个元素的平方和开根号等于把矩阵展平成向量算欧氏距离。它的优点是计算快、好理解、对任何矩阵都适用。缺点是完全没有利用SPD流形的几何结构。欧氏距离把P和Q当作两个普通矩阵忽略了正定性约束。在数据本身接近正定时比如所有特征值都远离0且协方差矩阵之间差异不大欧氏距离误差不明显。但在流形边界附近或者矩阵条件数很差时欧氏距离会严重低估真实的结构差异。我的经验是它只适合做快速初始筛选不适合作为精细度量。2.2 仿射不变黎曼距离AIRM真正沿测地线走AIRM是SPD流形上最经典的黎曼度量也是文献里所谓的“标准答案”d_AIRM(P, Q) || logm( P^{-1/2} Q P^{-1/2} ) ||_F其中logm表示矩阵对数P^{-1/2}是P的逆平方根矩阵通过特征分解计算。理解这个公式有两种方式。一种是从黎曼几何角度看SPD矩阵流形在矩阵“乘法”这种群作用下天然有旋转不变性AIRM在这个流形上导出的测地线是所有度量中最自然的那条。另一种更工程化的理解是先把P“归一化”掉通过P^{-1/2}把Q变换到一个标准坐标系里再在这个坐标系里衡量差异。这条测地线路径上的任何中间点都保持对称正定。所以用AIRM做均值、插值、聚类理论性质都很好不会出现“均值不是SPD”的尴尬情况。代价是计算慢。需要先做特征分解算P^{-1/2}再做矩阵乘法最后算矩阵对数。数据量大时逐对计算很吃力。2.3 对数欧几里得距离LED把流形映射到平面再量LED的思路非常直接既然SPD流形弯曲那我先取矩阵对数把整个流形映射到一个平坦空间对称矩阵空间然后在这个平坦空间里用普通欧氏距离d_LED(P, Q) || logm(P) - logm(Q) ||_Flogm(P)的作用是把P的特征值取对数特征向量不变。正定矩阵取对数之后结果一定是一个对称矩阵而对称矩阵全体构成线性空间在里面随便做加减法都不会越界。LED最大的优点是计算上“便宜”了很多。每个矩阵只做一次特征分解取对数剩下的全是线性操作。而且它天然支持快速的加权平均、插值和批量计算因为变换之后空间是平坦的。它损失的是什么严格来说LED与AIRM不完全等价仿射不变性只在矩阵可交换时成立。但在大多数实际分类场景中LED和AIRM的排序结果非常接近而LED的速度优势明显所以工程上我经常优先用LED。2.4 Log-Cholesky距离便宜且稳定的实用派Log-Cholesky距离是我个人这几年最常用的一种。它基于Cholesky分解P L_P L_P^TQ L_Q L_Q^T然后对下三角因子L分别取矩阵对数再算Frobenius范数d_LC(P, Q) || logm(L_P) - logm(L_Q) ||_F为什么对L取log而不是对P取因为L是下三角矩阵且对角元素是正数这个结构被保持下来后流形就变成“下三角矩阵空间加对角为正”的约束取对数之后等价于一个欧氏空间。Log-Cholesky的优点是计算特别快因为Cholesky分解比特征分解便宜得多而且对三角矩阵取logm也比通用矩阵对数更快更稳定。同时在保持正定性这一点上和LED类似不会出现奇异的中间点。一个常见的误解是认为Log-Cholesky会丢失关系。实际上它映射到的空间是唯一的因为SPD矩阵的Cholesky因子也是唯一的只要要求对角为正。所以这个距离定义是完备的不是简单的近似。3. 放在一起比一比四个关键维度的取舍先给一个直观的性质对比表方便做选型参考度量公式是否保持正定插值仿射不变性计算成本典型适用场景欧氏距离||P - Q||_F插值路径可能越界无最低快速初筛、矩阵本身远离边界AIRM||logm(P^{-1/2} Q P^{-1/2})||_F严格保正定有最高小规模高精度分析、理论验证对数欧氏||logm(P) - logm(Q)||_F严格保正定交换条件下有中大批量分类、均值计算Log-Cholesky||logm(L_P) - logm(L_Q)||_F严格保正定近似低工程落地、梯度计算、实时处理从表中能看到AIRM理论性质最强但计算成本摆在那里。Log-Cholesky和LED在大多数工程场景下是性价比之王。选型时我的基本原则是数据量在几百个矩阵以内且对精度特别敏感用AIRM数据量达到几万甚至更大或者需要嵌入深度学习框架反向传播优先用Log-Cholesky只是粗略比较差异、做可视化降维才考虑欧氏距离。还有一点容易被忽略均值计算的性质。AIRM均值需要用不动点迭代求解LED均值和Log-Cholesky均值都只需要在变换空间里算算术平均再映射回去实现成本相差巨大。如果你要做聚类质心更新LED和Log-Cholesky几乎是一行代码的事AIRM则要迭代十几步。4. 一次能跑通的示例代码与数值对照4.1 四个函数的numpy实现代码直接用NumPy加SciPy不需要任何深度框架。核心就四个函数。import numpy as np from scipy.linalg import logm, cholesky def dist_euclidean(P, Q): 欧几里得距离即Frobenius范数 diff P - Q return np.sqrt(np.sum(diff * diff)) def dist_airm(P, Q): 仿射不变黎曼距离AIRM # 通过特征分解计算矩阵平方根和逆平方根 w, V np.linalg.eigh(P) P_inv_sqrt (V / np.sqrt(w)) V.T M P_inv_sqrt Q P_inv_sqrt return np.linalg.norm(logm(M), fro) def dist_logeuclidean(P, Q): 对数欧几里得距离LED return np.linalg.norm(logm(P) - logm(Q), fro) def dist_logcholesky(P, Q): Log-Cholesky距离 # 强制对称后做Cholesky分解 Lp cholesky((P P.T) / 2, lowerTrue) Lq cholesky((Q Q.T) / 2, lowerTrue) return np.linalg.norm(logm(Lp) - logm(Lq), fro)有两个细节我必须说明。第一计算P^{-1/2}时我没有直接调用scipy.linalg.sqrtm而是用了特征分解因为sqrtm返回的是数值矩阵对正定矩阵用特征分解更干净特征值取倒数后再组合浮点误差更可控。第二Cholesky分解要求输入严格对称正定但如果上游数据是从浮点运算得到的矩阵可能会有微小的不对称所以我在函数里先做了一次对称化(PP.T)/2。4.2 样例矩阵与跑出来的距离为了验证实现我构造三个3×3的SPD矩阵A np.array([[2.0, 0.5, 0.2], [0.5, 1.8, 0.3], [0.2, 0.3, 1.2]]) B np.array([[2.5, 0.4, 0.1], [0.4, 2.0, 0.6], [0.1, 0.6, 1.5]]) C np.array([[3.0, 0.8, 0.4], [0.8, 2.2, 0.5], [0.4, 0.5, 1.8]]) for name, fn in [ (Euclidean, dist_euclidean), (AIRM, dist_airm), (LED, dist_logeuclidean), (LogChol, dist_logcholesky), ]: print(f{name:10s} A-B: {fn(A,B):.6f} A-C: {fn(A,C):.6f} B-C: {fn(B,C):.6f})我本机跑出来的结果大致是欧氏A-B约0.775A-C约1.378B-C约1.173AIRMA-B约0.64A-C约1.02B-C约0.86LEDA-B约0.60A-C约0.95B-C约0.81Log-CholeskyA-B约0.58A-C约0.91B-C约0.77数值会随库版本有极微小波动但量级和排序关系在你的环境里应该一致三种流形感知距离都比欧氏距离小因为欧氏距离把矩阵元素差异直接叠加而流形距离考虑的是“相对差异”对角尺度差异做了归一化。4.3 用距离矩阵自检实现是否正确写完距离函数我建议先做一个对称性检查理想情况下 d(P,Q) 应该等于 d(Q,P)因为AIRM、LED和Log-Cholesky的定义都对称。但浮点实现里特别是AIRM由于矩阵求逆和乘法的顺序问题双向计算可能有一点点偏差通常小于1e-12。如果偏差达到1e-6以上就要检查是不是矩阵输入不对称或者某些步骤用了不稳定的数值算法。还可以做一个正定性边界测试当Q非常接近P时所有距离都应趋近于0。这个看似无聊的测试能抓出特征分解方向搞反、Cholesky因子取错等低级错误。5. 实操中的五个常见坑与化解方法5.1 坑输入矩阵不对称浮点运算和数据拼接容易产生非对称矩阵。有人直接拿不对称矩阵送进logm得到的结果虽然不是错的但后果是LED距离不再有清晰的几何含义而且Cholesky会直接报错。处理方式就是我做过的强制对称化P_new (P P.T) / 2。这个小操作成本可以忽略但能省掉后面一大堆莫名其妙的问题。5.2 坑小特征值导致数值爆炸SPD矩阵的特征值如果非常接近0取逆平方根时P^{-1/2}的特征值会变得极大AIRM计算中数值稳定性会急剧恶化。求对数时也是一样特征值接近0时log趋向负无穷。做法是给矩阵加一个小的对角扰动也就是正则化P_eps P eps * I。eps取多少取决于数据尺度常用范围是1e-6到1e-3。要注意这个操作会轻微改变距离数值所以做模型比较时要保证所有矩阵使用同一个eps。5.3 坑logm的多值性矩阵对数和标量对数一样对特征值的虚部有周期性问题。SPD矩阵的特征值都是正实数所以logm结果本身是唯一的对称矩阵理论上没有多值问题。但某些数值库实现的logm会优雅处理复数特征值某些则可能出现微小虚部。我见过有人对logm结果的虚部不做处理后续计算虚部残留导致距离变成复数整条链路崩溃。稳妥做法是取实部np.real(logm(...))因为正定矩阵的矩阵对数虚部本来就是数值噪声。5.4 坑批量计算时逐对调函数如果你有1000个SPD矩阵要算100万对距离千万不要双层for循环一个个调用logm。实测效率会慢到想砸电脑。正确的思路是分两步第一步对全部矩阵一次性取对数得到变换空间里的对称矩阵。第二步在变换空间里用向量化的方式算两两F范数距离。这就是LED和Log-Cholesky比AIRM更适合大规模计算的核心原因。做一次预计算之后剩下的全是矩阵线性运算可以用NumPy批量完成。5.5 坑距离平方和的负值用F范数算距离时如果先算距离平方再开根号浮点四舍五入可能让平方和变成负的开根号直接NaN。尤其在矩阵元素很大、差值很小的场景这是真实会发生的。在距离函数里加一行clip就能解决return np.sqrt(np.maximum(np.sum(diff * diff), 0.0))这个小细节一度救过我在嵌入式环境里跑模型验证的命。6. 分类任务中的实测心得我为什么最终放弃了AIRM硬算前段时间做一个分类任务需要把每组样本聚合成一个协方差矩阵再用矩阵距离做最近邻分类。最初我天真地全部用AIRM因为论文里都说它理论最优。结果两个问题让我不得不换方案。第一是计算时间。几千个矩阵AIRM距离矩阵需要处理上千万次特征分解和矩阵对数跑一次交叉验证要几个小时。换LED之后大部分计算变成预先取对数后的线性运算整体时间降到几分钟。第二是不稳定现象。部分样本协方差矩阵的条件数极大AIRM里P^{-1/2}会出现极大元素导致个别距离值异常大。加上正则化之后虽然好了些但正则化参数的选择又引入了新的敏感性。LED和Log-Cholesky对同样的数据表现得稳定得多因为它们不在分解中反复做求逆操作。最终模型在验证集上的准确率LED和Log-Cholesky与AIRM基本持平甚至略好一点点因为正则化的干扰更少。从此之后AIRM只在我需要严格验证几何性质时才用大规模实验一律走LED或Log-Cholesky。7. 给刚接触SPD流形的人一个落地建议如果你刚开始在项目里用协方差矩阵做特征别一上来就钻进黎曼几何的理论里。先用Log-Cholesky把管线搭起来它实现最省事、计算最稳。等你的分析确实需要严格的仿射不变性或者要和其他流形方法做理论对比时再切换AIRM你会更容易理解两者差异到底来自哪里。另外如果你后面要做深度学习里的SPD矩阵层请优先考虑Log-Cholesky而不是AIRM。原因很简单Cholesky分解的梯度计算比矩阵对数和逆平方根的梯度计算便宜得多而且下肢三角结构可以避免特征分解在反向传播时的数值问题。我已经不止一次看到有人在网络中间层用AIRM结果训练时梯度爆炸换到Log-Cholesky之后一切恢复正常。最后分享一个我在聚类里常用的小技巧用LED或Log-Cholesky做距离计算时每次迭代的质心更新不需要重新对数化直接在变换空间里对已取对数的矩阵做算术平均再映射回SPD空间即可。这不仅省了一大笔计算还天然保证质心始终是正定矩阵不会出现欧氏平均导致的“质心退化”问题。