提到對(duì)稱正定矩陣上的距離計(jì)算很多人第一反應(yīng)就是直接用矩陣差的Frobenius范數(shù)。我在一個(gè)給協(xié)方差矩陣做分類的項(xiàng)目里也這么干過直到發(fā)現(xiàn)用歐氏距離插值出來的均值矩陣連正定性都保不住才意識(shí)到“距離”這個(gè)詞在SPD流形上遠(yuǎn)不是套個(gè)范數(shù)公式那么簡(jiǎn)單。這篇文章要聊的是SPD流形上最常用的四種距離度量歐幾里得距離、仿射不變黎曼距離AIRM、對(duì)數(shù)歐幾里得距離LED和Log-Cholesky距離。它們各自有不同的幾何含義、計(jì)算成本和適用范圍直接決定了均值、插值、聚類這些下游操作的結(jié)果。適合正在做流形學(xué)習(xí)、腦電信號(hào)分類、圖像區(qū)域描述子匹配或者剛接觸協(xié)方差矩陣幾何性質(zhì)的讀者參考。我會(huì)把每種度的定義、直覺、代碼實(shí)現(xiàn)和實(shí)際踩過的坑一起說清楚。1. 為什么“在SPD流形上算距離”不是一件小事1.1 正定約束制造了一個(gè)彎曲空間先明確這里討論的對(duì)象。一個(gè)對(duì)稱正定矩陣P要求它對(duì)稱、且對(duì)任意非零向量x都有 x^T P x 0。這個(gè)約束看起來只是代數(shù)上的要求但它把合法矩陣集合限制成了一個(gè)彎曲的流形而不是一個(gè)平直的向量空間。舉個(gè)直覺例子所有2×2的對(duì)稱矩陣可以想象成一個(gè)3維平直空間三個(gè)獨(dú)立元素但其中只有同時(shí)滿足行列式大于0、左上角元素大于0的那部分才是SPD矩陣。這個(gè)區(qū)域像圓錐內(nèi)部邊界是半正定矩陣。在這個(gè)區(qū)域內(nèi)走直線很可能一步跨出邊界。所以當(dāng)你計(jì)算“兩個(gè)協(xié)方差矩陣相差多少”時(shí)不能簡(jiǎn)單把它們當(dāng)向量空間里的點(diǎn)來量。流形的彎曲幾何決定了“直線”不是我們熟悉的直線“距離”也不是我們熟悉的距離。1.2 歐氏線性操作會(huì)把正定性搞丟我在項(xiàng)目里踩的第一個(gè)坑就是均值外插。假設(shè)手頭有P和Q兩個(gè)SPD矩陣用普通算術(shù)平均得到 (PQ)/2這個(gè)結(jié)果恰好還是對(duì)稱正定的問題不明顯。但一旦做線性插值 (1-t)P tQt在0到1之間還好t一旦大于1或者小于0正定性立刻無法保證。更嚴(yán)重的場(chǎng)景是計(jì)算測(cè)地線或者做流形上的均值。歐氏均值只考慮平直空間的平均完全不理會(huì)流形曲率結(jié)果常常外插到非正定區(qū)域。這就是為什么后來我徹底不用歐氏方式在SPD矩陣集上做平均和距離計(jì)算。1.3 距離度量是下游一切操作的地基選擇哪種距離會(huì)連鎖影響聚類、分類、插值、均值甚至梯度計(jì)算。同一個(gè)數(shù)據(jù)集用歐氏距離和用仿射不變黎曼距離算出來的聚類中心完全不同同一條“插值路徑”歐氏路徑中間點(diǎn)可能退化而黎曼測(cè)地線路徑上的每一個(gè)中間點(diǎn)都保持正定。所以“距離計(jì)算”這個(gè)問題本質(zhì)上是給整個(gè)流形幾何定調(diào)子。理解了這一點(diǎn)下面四種度量就好理解了它們不是數(shù)學(xué)游戲而是給SPD流形提供的四種不同“尺子”。2. 四種常用度量各自的定義與幾何直覺2.1 歐幾里得距離最直接但最“不流形”的尺子公式很簡(jiǎn)單d_F(P, Q) ||P - Q||_F sqrt( Σ_ij (P_ij - Q_ij)^2 )這里我用Frobenius范數(shù)也就是矩陣每個(gè)元素的平方和開根號(hào)等于把矩陣展平成向量算歐氏距離。它的優(yōu)點(diǎn)是計(jì)算快、好理解、對(duì)任何矩陣都適用。缺點(diǎn)是完全沒有利用SPD流形的幾何結(jié)構(gòu)。歐氏距離把P和Q當(dāng)作兩個(gè)普通矩陣忽略了正定性約束。在數(shù)據(jù)本身接近正定時(shí)比如所有特征值都遠(yuǎn)離0且協(xié)方差矩陣之間差異不大歐氏距離誤差不明顯。但在流形邊界附近或者矩陣條件數(shù)很差時(shí)歐氏距離會(huì)嚴(yán)重低估真實(shí)的結(jié)構(gòu)差異。我的經(jīng)驗(yàn)是它只適合做快速初始篩選不適合作為精細(xì)度量。2.2 仿射不變黎曼距離AIRM真正沿測(cè)地線走AIRM是SPD流形上最經(jīng)典的黎曼度量也是文獻(xiàn)里所謂的“標(biāo)準(zhǔn)答案”d_AIRM(P, Q) || logm( P^{-1/2} Q P^{-1/2} ) ||_F其中l(wèi)ogm表示矩陣對(duì)數(shù)P^{-1/2}是P的逆平方根矩陣通過特征分解計(jì)算。理解這個(gè)公式有兩種方式。一種是從黎曼幾何角度看SPD矩陣流形在矩陣“乘法”這種群作用下天然有旋轉(zhuǎn)不變性AIRM在這個(gè)流形上導(dǎo)出的測(cè)地線是所有度量中最自然的那條。另一種更工程化的理解是先把P“歸一化”掉通過P^{-1/2}把Q變換到一個(gè)標(biāo)準(zhǔn)坐標(biāo)系里再在這個(gè)坐標(biāo)系里衡量差異。這條測(cè)地線路徑上的任何中間點(diǎn)都保持對(duì)稱正定。所以用AIRM做均值、插值、聚類理論性質(zhì)都很好不會(huì)出現(xiàn)“均值不是SPD”的尷尬情況。代價(jià)是計(jì)算慢。需要先做特征分解算P^{-1/2}再做矩陣乘法最后算矩陣對(duì)數(shù)。數(shù)據(jù)量大時(shí)逐對(duì)計(jì)算很吃力。2.3 對(duì)數(shù)歐幾里得距離LED把流形映射到平面再量LED的思路非常直接既然SPD流形彎曲那我先取矩陣對(duì)數(shù)把整個(gè)流形映射到一個(gè)平坦空間對(duì)稱矩陣空間然后在這個(gè)平坦空間里用普通歐氏距離d_LED(P, Q) || logm(P) - logm(Q) ||_Flogm(P)的作用是把P的特征值取對(duì)數(shù)特征向量不變。正定矩陣取對(duì)數(shù)之后結(jié)果一定是一個(gè)對(duì)稱矩陣而對(duì)稱矩陣全體構(gòu)成線性空間在里面隨便做加減法都不會(huì)越界。LED最大的優(yōu)點(diǎn)是計(jì)算上“便宜”了很多。每個(gè)矩陣只做一次特征分解取對(duì)數(shù)剩下的全是線性操作。而且它天然支持快速的加權(quán)平均、插值和批量計(jì)算因?yàn)樽儞Q之后空間是平坦的。它損失的是什么嚴(yán)格來說LED與AIRM不完全等價(jià)仿射不變性只在矩陣可交換時(shí)成立。但在大多數(shù)實(shí)際分類場(chǎng)景中LED和AIRM的排序結(jié)果非常接近而LED的速度優(yōu)勢(shì)明顯所以工程上我經(jīng)常優(yōu)先用LED。2.4 Log-Cholesky距離便宜且穩(wěn)定的實(shí)用派Log-Cholesky距離是我個(gè)人這幾年最常用的一種。它基于Cholesky分解P L_P L_P^TQ L_Q L_Q^T然后對(duì)下三角因子L分別取矩陣對(duì)數(shù)再算Frobenius范數(shù)d_LC(P, Q) || logm(L_P) - logm(L_Q) ||_F為什么對(duì)L取log而不是對(duì)P取因?yàn)長(zhǎng)是下三角矩陣且對(duì)角元素是正數(shù)這個(gè)結(jié)構(gòu)被保持下來后流形就變成“下三角矩陣空間加對(duì)角為正”的約束取對(duì)數(shù)之后等價(jià)于一個(gè)歐氏空間。Log-Cholesky的優(yōu)點(diǎn)是計(jì)算特別快因?yàn)镃holesky分解比特征分解便宜得多而且對(duì)三角矩陣取logm也比通用矩陣對(duì)數(shù)更快更穩(wěn)定。同時(shí)在保持正定性這一點(diǎn)上和LED類似不會(huì)出現(xiàn)奇異的中間點(diǎn)。一個(gè)常見的誤解是認(rèn)為L(zhǎng)og-Cholesky會(huì)丟失關(guān)系。實(shí)際上它映射到的空間是唯一的因?yàn)镾PD矩陣的Cholesky因子也是唯一的只要要求對(duì)角為正。所以這個(gè)距離定義是完備的不是簡(jiǎn)單的近似。3. 放在一起比一比四個(gè)關(guān)鍵維度的取舍先給一個(gè)直觀的性質(zhì)對(duì)比表方便做選型參考度量公式是否保持正定插值仿射不變性計(jì)算成本典型適用場(chǎng)景歐氏距離||P - Q||_F插值路徑可能越界無最低快速初篩、矩陣本身遠(yuǎn)離邊界AIRM||logm(P^{-1/2} Q P^{-1/2})||_F嚴(yán)格保正定有最高小規(guī)模高精度分析、理論驗(yàn)證對(duì)數(shù)歐氏||logm(P) - logm(Q)||_F嚴(yán)格保正定交換條件下有中大批量分類、均值計(jì)算Log-Cholesky||logm(L_P) - logm(L_Q)||_F嚴(yán)格保正定近似低工程落地、梯度計(jì)算、實(shí)時(shí)處理從表中能看到AIRM理論性質(zhì)最強(qiáng)但計(jì)算成本擺在那里。Log-Cholesky和LED在大多數(shù)工程場(chǎng)景下是性價(jià)比之王。選型時(shí)我的基本原則是數(shù)據(jù)量在幾百個(gè)矩陣以內(nèi)且對(duì)精度特別敏感用AIRM數(shù)據(jù)量達(dá)到幾萬甚至更大或者需要嵌入深度學(xué)習(xí)框架反向傳播優(yōu)先用Log-Cholesky只是粗略比較差異、做可視化降維才考慮歐氏距離。還有一點(diǎn)容易被忽略均值計(jì)算的性質(zhì)。AIRM均值需要用不動(dòng)點(diǎn)迭代求解LED均值和Log-Cholesky均值都只需要在變換空間里算算術(shù)平均再映射回去實(shí)現(xiàn)成本相差巨大。如果你要做聚類質(zhì)心更新LED和Log-Cholesky幾乎是一行代碼的事AIRM則要迭代十幾步。4. 一次能跑通的示例代碼與數(shù)值對(duì)照4.1 四個(gè)函數(shù)的numpy實(shí)現(xiàn)代碼直接用NumPy加SciPy不需要任何深度框架。核心就四個(gè)函數(shù)。import numpy as np from scipy.linalg import logm, cholesky def dist_euclidean(P, Q): 歐幾里得距離即Frobenius范數(shù) diff P - Q return np.sqrt(np.sum(diff * diff)) def dist_airm(P, Q): 仿射不變黎曼距離AIRM # 通過特征分解計(jì)算矩陣平方根和逆平方根 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): 對(duì)數(shù)歐幾里得距離LED return np.linalg.norm(logm(P) - logm(Q), fro) def dist_logcholesky(P, Q): Log-Cholesky距離 # 強(qiáng)制對(duì)稱后做Cholesky分解 Lp cholesky((P P.T) / 2, lowerTrue) Lq cholesky((Q Q.T) / 2, lowerTrue) return np.linalg.norm(logm(Lp) - logm(Lq), fro)有兩個(gè)細(xì)節(jié)我必須說明。第一計(jì)算P^{-1/2}時(shí)我沒有直接調(diào)用scipy.linalg.sqrtm而是用了特征分解因?yàn)閟qrtm返回的是數(shù)值矩陣對(duì)正定矩陣用特征分解更干凈特征值取倒數(shù)后再組合浮點(diǎn)誤差更可控。第二Cholesky分解要求輸入嚴(yán)格對(duì)稱正定但如果上游數(shù)據(jù)是從浮點(diǎn)運(yùn)算得到的矩陣可能會(huì)有微小的不對(duì)稱所以我在函數(shù)里先做了一次對(duì)稱化(PP.T)/2。4.2 樣例矩陣與跑出來的距離為了驗(yàn)證實(shí)現(xiàn)我構(gòu)造三個(gè)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})我本機(jī)跑出來的結(jié)果大致是歐氏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數(shù)值會(huì)隨庫版本有極微小波動(dòng)但量級(jí)和排序關(guān)系在你的環(huán)境里應(yīng)該一致三種流形感知距離都比歐氏距離小因?yàn)闅W氏距離把矩陣元素差異直接疊加而流形距離考慮的是“相對(duì)差異”對(duì)角尺度差異做了歸一化。4.3 用距離矩陣自檢實(shí)現(xiàn)是否正確寫完距離函數(shù)我建議先做一個(gè)對(duì)稱性檢查理想情況下 d(P,Q) 應(yīng)該等于 d(Q,P)因?yàn)锳IRM、LED和Log-Cholesky的定義都對(duì)稱。但浮點(diǎn)實(shí)現(xiàn)里特別是AIRM由于矩陣求逆和乘法的順序問題雙向計(jì)算可能有一點(diǎn)點(diǎn)偏差通常小于1e-12。如果偏差達(dá)到1e-6以上就要檢查是不是矩陣輸入不對(duì)稱或者某些步驟用了不穩(wěn)定的數(shù)值算法。還可以做一個(gè)正定性邊界測(cè)試當(dāng)Q非常接近P時(shí)所有距離都應(yīng)趨近于0。這個(gè)看似無聊的測(cè)試能抓出特征分解方向搞反、Cholesky因子取錯(cuò)等低級(jí)錯(cuò)誤。5. 實(shí)操中的五個(gè)常見坑與化解方法5.1 坑輸入矩陣不對(duì)稱浮點(diǎn)運(yùn)算和數(shù)據(jù)拼接容易產(chǎn)生非對(duì)稱矩陣。有人直接拿不對(duì)稱矩陣送進(jìn)logm得到的結(jié)果雖然不是錯(cuò)的但后果是LED距離不再有清晰的幾何含義而且Cholesky會(huì)直接報(bào)錯(cuò)。處理方式就是我做過的強(qiáng)制對(duì)稱化P_new (P P.T) / 2。這個(gè)小操作成本可以忽略但能省掉后面一大堆莫名其妙的問題。5.2 坑小特征值導(dǎo)致數(shù)值爆炸SPD矩陣的特征值如果非常接近0取逆平方根時(shí)P^{-1/2}的特征值會(huì)變得極大AIRM計(jì)算中數(shù)值穩(wěn)定性會(huì)急劇惡化。求對(duì)數(shù)時(shí)也是一樣特征值接近0時(shí)log趨向負(fù)無窮。做法是給矩陣加一個(gè)小的對(duì)角擾動(dòng)也就是正則化P_eps P eps * I。eps取多少取決于數(shù)據(jù)尺度常用范圍是1e-6到1e-3。要注意這個(gè)操作會(huì)輕微改變距離數(shù)值所以做模型比較時(shí)要保證所有矩陣使用同一個(gè)eps。5.3 坑logm的多值性矩陣對(duì)數(shù)和標(biāo)量對(duì)數(shù)一樣對(duì)特征值的虛部有周期性問題。SPD矩陣的特征值都是正實(shí)數(shù)所以logm結(jié)果本身是唯一的對(duì)稱矩陣?yán)碚撋蠜]有多值問題。但某些數(shù)值庫實(shí)現(xiàn)的logm會(huì)優(yōu)雅處理復(fù)數(shù)特征值某些則可能出現(xiàn)微小虛部。我見過有人對(duì)logm結(jié)果的虛部不做處理后續(xù)計(jì)算虛部殘留導(dǎo)致距離變成復(fù)數(shù)整條鏈路崩潰。穩(wěn)妥做法是取實(shí)部np.real(logm(...))因?yàn)檎ň仃嚨木仃噷?duì)數(shù)虛部本來就是數(shù)值噪聲。5.4 坑批量計(jì)算時(shí)逐對(duì)調(diào)函數(shù)如果你有1000個(gè)SPD矩陣要算100萬對(duì)距離千萬不要雙層for循環(huán)一個(gè)個(gè)調(diào)用logm。實(shí)測(cè)效率會(huì)慢到想砸電腦。正確的思路是分兩步第一步對(duì)全部矩陣一次性取對(duì)數(shù)得到變換空間里的對(duì)稱矩陣。第二步在變換空間里用向量化的方式算兩兩F范數(shù)距離。這就是LED和Log-Cholesky比AIRM更適合大規(guī)模計(jì)算的核心原因。做一次預(yù)計(jì)算之后剩下的全是矩陣線性運(yùn)算可以用NumPy批量完成。5.5 坑距離平方和的負(fù)值用F范數(shù)算距離時(shí)如果先算距離平方再開根號(hào)浮點(diǎn)四舍五入可能讓平方和變成負(fù)的開根號(hào)直接NaN。尤其在矩陣元素很大、差值很小的場(chǎng)景這是真實(shí)會(huì)發(fā)生的。在距離函數(shù)里加一行clip就能解決return np.sqrt(np.maximum(np.sum(diff * diff), 0.0))這個(gè)小細(xì)節(jié)一度救過我在嵌入式環(huán)境里跑模型驗(yàn)證的命。6. 分類任務(wù)中的實(shí)測(cè)心得我為什么最終放棄了AIRM硬算前段時(shí)間做一個(gè)分類任務(wù)需要把每組樣本聚合成一個(gè)協(xié)方差矩陣再用矩陣距離做最近鄰分類。最初我天真地全部用AIRM因?yàn)檎撐睦锒颊f它理論最優(yōu)。結(jié)果兩個(gè)問題讓我不得不換方案。第一是計(jì)算時(shí)間。幾千個(gè)矩陣AIRM距離矩陣需要處理上千萬次特征分解和矩陣對(duì)數(shù)跑一次交叉驗(yàn)證要幾個(gè)小時(shí)。換LED之后大部分計(jì)算變成預(yù)先取對(duì)數(shù)后的線性運(yùn)算整體時(shí)間降到幾分鐘。第二是不穩(wěn)定現(xiàn)象。部分樣本協(xié)方差矩陣的條件數(shù)極大AIRM里P^{-1/2}會(huì)出現(xiàn)極大元素導(dǎo)致個(gè)別距離值異常大。加上正則化之后雖然好了些但正則化參數(shù)的選擇又引入了新的敏感性。LED和Log-Cholesky對(duì)同樣的數(shù)據(jù)表現(xiàn)得穩(wěn)定得多因?yàn)樗鼈儾辉诜纸庵蟹磸?fù)做求逆操作。最終模型在驗(yàn)證集上的準(zhǔn)確率LED和Log-Cholesky與AIRM基本持平甚至略好一點(diǎn)點(diǎn)因?yàn)檎齽t化的干擾更少。從此之后AIRM只在我需要嚴(yán)格驗(yàn)證幾何性質(zhì)時(shí)才用大規(guī)模實(shí)驗(yàn)一律走LED或Log-Cholesky。7. 給剛接觸SPD流形的人一個(gè)落地建議如果你剛開始在項(xiàng)目里用協(xié)方差矩陣做特征別一上來就鉆進(jìn)黎曼幾何的理論里。先用Log-Cholesky把管線搭起來它實(shí)現(xiàn)最省事、計(jì)算最穩(wěn)。等你的分析確實(shí)需要嚴(yán)格的仿射不變性或者要和其他流形方法做理論對(duì)比時(shí)再切換AIRM你會(huì)更容易理解兩者差異到底來自哪里。另外如果你后面要做深度學(xué)習(xí)里的SPD矩陣層請(qǐng)優(yōu)先考慮Log-Cholesky而不是AIRM。原因很簡(jiǎn)單Cholesky分解的梯度計(jì)算比矩陣對(duì)數(shù)和逆平方根的梯度計(jì)算便宜得多而且下肢三角結(jié)構(gòu)可以避免特征分解在反向傳播時(shí)的數(shù)值問題。我已經(jīng)不止一次看到有人在網(wǎng)絡(luò)中間層用AIRM結(jié)果訓(xùn)練時(shí)梯度爆炸換到Log-Cholesky之后一切恢復(fù)正常。最后分享一個(gè)我在聚類里常用的小技巧用LED或Log-Cholesky做距離計(jì)算時(shí)每次迭代的質(zhì)心更新不需要重新對(duì)數(shù)化直接在變換空間里對(duì)已取對(duì)數(shù)的矩陣做算術(shù)平均再映射回SPD空間即可。這不僅省了一大筆計(jì)算還天然保證質(zhì)心始終是正定矩陣不會(huì)出現(xiàn)歐氏平均導(dǎo)致的“質(zhì)心退化”問題。