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