从理论到实践Python实现Grassmann流形上的子空间距离计算在机器学习与数据科学的实际应用中我们常常需要比较不同数据集或模型生成的特征子空间。传统欧氏距离无法捕捉子空间之间的几何关系而Grassmann流形提供了一套优雅的数学框架来解决这个问题。本文将带您从零开始用Python实现Grassmann流形上的五种核心距离度量。1. Grassmann流形基础与准备工作Grassmann流形G(m,D)表示D维空间中所有m维子空间的集合。想象一下当我们需要比较两个人脸识别模型提取的特征空间或者分析不同时间段用户行为特征的漂移时Grassmann流形就是理想的数学工具。环境配置要求import numpy as np from scipy.linalg import svd, orth from sklearn.datasets import load_digits import matplotlib.pyplot as plt生成测试子空间的实用函数def generate_random_subspace(D, m): 生成随机D维空间中的m维子空间 A np.random.randn(D, D) return orth(A)[:, :m] # 取前m个正交基注意所有子空间必须用标准正交基表示即Y.T Y I。实际应用中可通过SVD获取数据的正交基。2. 核心距离度量实现2.1 主角度计算所有度量的基础主角度是理解子空间关系的关键。下面函数计算两个子空间之间的主角度def principal_angles(Y1, Y2): 计算两个子空间之间的主角度(弧度) U, s, Vt svd(Y1.T Y2) return np.arccos(np.clip(s, -1, 1)) # 确保数值稳定性测试示例Y1 generate_random_subspace(100, 10) Y2 generate_random_subspace(100, 10) angles principal_angles(Y1, Y2) print(f主角度(度): {np.degrees(angles)})2.2 投影度量Projection Metric投影度量衡量的是子空间投影差异的Frobenius范数def projection_metric(Y1, Y2): 计算投影度量距离 cos_theta svd(Y1.T Y2, compute_uvFalse) return np.sqrt(min(Y1.shape[1], Y2.shape[1]) - np.sum(cos_theta**2))特性分析对小的角度变化敏感取值范围[0, √min(m1,m2)]计算复杂度O(D*m^2)2.3 最大/最小相关距离这两种度量分别关注最大和最小主角度def max_correlation(Y1, Y2): 基于最大主角度的距离 cos_theta1 svd(Y1.T Y2, compute_uvFalse)[0] return np.sqrt(1 - cos_theta1**2) def min_correlation(Y1, Y2): 基于最小主角度的距离 cos_thetam svd(Y1.T Y2, compute_uvFalse)[-1] return np.sqrt(1 - cos_thetam**2)2.4 Binet-Cauchy度量考虑所有主角度的乘积关系def binet_cauchy(Y1, Y2): Binet-Cauchy距离计算 cos_theta svd(Y1.T Y2, compute_uvFalse) return 1 - np.prod(cos_theta**2)**(1/2)2.5 Procrustes度量弦距离寻找最优旋转对齐后的距离def procrustes_metric(Y1, Y2): Procrustes距离计算 U, _, Vt svd(Y1.T Y2) return np.linalg.norm(Y1 U - Y2 Vt.T, fro)3. 实际应用案例分析3.1 MNIST数字分类中的特征空间比较比较PCA和LDA降维后的特征空间digits load_digits() X digits.data y digits.target # PCA子空间 U_pca, _, _ svd(X - X.mean(axis0)) Y_pca U_pca[:, :20] # LDA子空间简化版 class_means [] for i in range(10): class_means.append(X[yi].mean(axis0)) S_b np.cov(np.array(class_means).T) _, U_lda np.linalg.eigh(S_b) Y_lda U_lda[:, -20:] # 计算各种距离 metrics { Projection: projection_metric, Max Correlation: max_correlation, Min Correlation: min_correlation, Binet-Cauchy: binet_cauchy, Procrustes: procrustes_metric } results {name: metric(Y_pca, Y_lda) for name, metric in metrics.items()}3.2 距离度量对比表度量类型计算复杂度敏感度适用场景投影度量O(Dm^2)全局一般子空间比较最大相关O(Dm^2)最大角度异常检测最小相关O(Dm^2)最小角度保守相似性评估Binet-CauchyO(Dm^2)乘积关系整体相似性ProcrustesO(Dm^2)最优对齐需要旋转不变性的场景4. 数值稳定性与优化技巧实际应用中需要注意的常见问题1. 维度不匹配处理def safe_metric(Y1, Y2): 处理维度不匹配的安全距离计算 m min(Y1.shape[1], Y2.shape[1]) cos_theta svd(Y1[:,:m].T Y2[:,:m], compute_uvFalse) return np.sqrt(m - np.sum(cos_theta**2))2. 正则化技巧添加小的单位矩阵防止奇异regularized Y1.T Y2 1e-6 * np.eye(Y1.shape[1])3. 并行计算优化from joblib import Parallel, delayed def batch_distance(subspaces, metric): 批量计算子空间距离 return Parallel(n_jobs-1)( delayed(metric)(Y1, Y2) for i, Y1 in enumerate(subspaces) for Y2 in subspaces[i1:] )5. 进阶应用与扩展5.1 流形上的均值计算Grassmann流形上的Karcher均值实现def grassmann_mean(subspaces, max_iter100, tol1e-6): 计算多个子空间的Karcher均值 mean subspaces[0] for _ in range(max_iter): tangents [] for Y in subspaces: U, _, Vt svd(mean.T Y) tangents.append(Y Vt.T U.T - mean) update np.mean(tangents, axis0) mean orth(mean update) if np.linalg.norm(update) tol: break return mean5.2 动态子空间跟踪对于随时间变化的子空间可以实现滑动窗口跟踪class SubspaceTracker: def __init__(self, window_size10): self.window [] self.window_size window_size def update(self, new_subspace): if len(self.window) self.window_size: self.window.pop(0) self.window.append(new_subspace) return self.compute_trend() def compute_trend(self): if len(self.window) 2: return 0 distances [ projection_metric(self.window[i], self.window[i1]) for i in range(len(self.window)-1) ] return np.mean(distances)在实际项目中我发现Procrustes度量对于旋转和尺度变化具有最好的鲁棒性特别是在处理不同来源的特征提取器时。而最小相关距离则更适合检测子空间之间的根本性差异。