高维数据相似度计算实战用Python实现Grassmann流形度量当我们需要比较两组高维数据的相似性时传统的欧氏距离往往力不从心。想象一下你手上有两组图像特征向量每组包含数百个样本每个样本都是上千维的向量。简单的点对点距离计算不仅计算量大更重要的是它无法捕捉数据背后的整体结构信息。这就是Grassmann流形大显身手的地方——它让我们能够比较两个高维子空间之间的形状差异。1. 为什么需要子空间相似度计算在机器学习和数据分析领域我们经常遇到需要比较两组高维数据的情况。比如跨域图像识别比较源域和目标域的特征分布时间序列分析评估不同时间段数据特征的稳定性模型评估对比不同模型提取的特征空间结构推荐系统衡量用户群体兴趣子空间的相似度传统方法如欧氏距离或余弦相似度只能计算单个向量之间的距离而Grassmann流形则允许我们比较整个子空间的几何结构。这种更高层次的比较往往能揭示数据间更本质的关系。import numpy as np from scipy.linalg import svd # 生成两组随机高维数据作为示例 np.random.seed(42) data1 np.random.randn(100, 50) # 第一组数据100个样本每个50维 data2 np.random.randn(80, 50) # 第二组数据80个样本每个50维2. Grassmann流形基础从理论到代码实现Grassmann流形G(m,D)表示D维空间中所有m维子空间的集合。简单来说如果我们有一组高维数据可以通过其主成分分析(PCA)得到一个低维子空间这个子空间就是Grassmann流形上的一个点。2.1 准备工作数据预处理与子空间提取在比较子空间之前我们需要从原始数据中提取代表性的子空间。通常使用PCA降维def extract_subspace(data, subspace_dim): 从数据中提取子空间基 # 中心化数据 data_centered data - np.mean(data, axis0) # 计算SVD _, s, vh svd(data_centered, full_matricesFalse) # 选择前subspace_dim个右奇异向量作为子空间基 subspace vh[:subspace_dim].T return subspace # 为两组数据提取10维子空间 subspace_dim 10 subspace1 extract_subspace(data1, subspace_dim) subspace2 extract_subspace(data2, subspace_dim)2.2 主角度计算子空间相似度的核心指标主角度(Principal Angles)是两个子空间之间的一系列角度类似于向量间的夹角但推广到了高维空间。计算主角度的Python实现def principal_angles(subspace1, subspace2): 计算两个子空间之间的主角度 # 计算两个子空间基的乘积 M subspace1.T subspace2 # 对M进行SVD分解 U, s, Vh svd(M, full_matricesFalse) # 奇异值就是主角度的余弦值 cos_angles s # 确保余弦值在[-1,1]范围内 cos_angles np.clip(cos_angles, -1.0, 1.0) # 计算角度(弧度) angles np.arccos(cos_angles) return angles # 计算两组数据子空间的主角度 angles principal_angles(subspace1, subspace2) print(f主角度(弧度): {angles}) print(f主角度(度): {np.degrees(angles)})3. Grassmann流形上的五种距离度量基于主角度我们可以定义多种Grassmann流形上的距离度量每种都有其特定的应用场景。3.1 投影度量(Projection Metric)投影度量是最常用的Grassmann距离之一定义为所有主角度正弦值的2-范数def projection_metric(angles): 计算投影度量距离 return np.linalg.norm(np.sin(angles)) proj_dist projection_metric(angles) print(f投影度量距离: {proj_dist:.4f})适用场景当我们需要考虑子空间在所有方向上的整体差异时投影度量是一个很好的选择。它在域适应(Domain Adaptation)和子空间聚类中表现良好。3.2 Binet-Cauchy度量Binet-Cauchy度量基于主角度余弦值的乘积def binet_cauchy_metric(angles): 计算Binet-Cauchy度量距离 cos_sq np.cos(angles)**2 return 1 - np.sqrt(np.prod(cos_sq)) bc_dist binet_cauchy_metric(angles) print(fBinet-Cauchy度量距离: {bc_dist:.4f})适用场景当子空间的主成分重要性不同时Binet-Cauchy度量能更好地反映主导方向的相似性。3.3 最大相关距离(Max Correlation Distance)只考虑最小的主角度(即最大的余弦值)def max_correlation_metric(angles): 计算最大相关距离 return np.sin(angles[0]) max_corr_dist max_correlation_metric(angles) print(f最大相关距离: {max_corr_dist:.4f})适用场景当我们只关心两个子空间最相似的方向时比如在特征选择或模型稳定性评估中。3.4 最小相关距离(Min Correlation Distance)考虑最大的主角度(即最小的余弦值)def min_correlation_metric(angles): 计算最小相关距离 return np.sin(angles[-1]) min_corr_dist min_correlation_metric(angles) print(f最小相关距离: {min_corr_dist:.4f})适用场景检测子空间之间的最不相似方向适用于异常检测或鲁棒性分析。3.5 Procrustes度量(Chordal Distance)Procrustes度量考虑了两个子空间基之间的最小Frobenius范数距离def procrustes_metric(subspace1, subspace2): 计算Procrustes度量距离 M subspace1.T subspace2 U, _, Vh svd(M, full_matricesFalse) R U Vh dist np.linalg.norm(subspace1 - subspace2 R.T, fro) return dist proc_dist procrustes_metric(subspace1, subspace2) print(fProcrustes度量距离: {proc_dist:.4f})适用场景当我们需要考虑子空间基的具体表示而不仅仅是它们的张成空间时Procrustes度量特别有用。4. 实际应用案例图像特征空间比较让我们通过一个实际例子来展示Grassmann流形距离的应用。假设我们有两组来自不同领域的图像特征源域自然场景图像(如COCO数据集)目标域医学X光图像我们想评估这两个领域的特征分布差异以决定是否需要域适应技术。# 假设我们已经提取了特征(实际应用中需要真实数据) source_features np.random.randn(500, 2048) # 500张自然图像每张2048维特征 target_features np.random.randn(300, 2048) # 300张医学图像每张2048维特征 # 提取子空间(降维到64维) subspace_source extract_subspace(source_features, 64) subspace_target extract_subspace(target_features, 64) # 计算各种距离 angles principal_angles(subspace_source, subspace_target) print(\n域适应场景下的子空间距离:) print(f- 投影度量: {projection_metric(angles):.4f}) print(f- Binet-Cauchy度量: {binet_cauchy_metric(angles):.4f}) print(f- Procrustes度量: {procrustes_metric(subspace_source, subspace_target):.4f})4.1 距离度量的选择指南不同的距离度量适用于不同的场景度量类型数学表达式适用场景计算复杂度投影度量‖sinθ‖₂整体子空间差异评估中等Binet-Cauchy1-∏cos²θ主导方向相似性中等最大相关sinθ₁最相似方向评估低最小相关sinθₘ最不相似方向评估低Procrustes‖Y₁-Y₂R‖F基对齐评估高5. 高级技巧与优化建议5.1 计算效率优化对于大规模高维数据直接计算SVD可能很耗时。可以采用以下优化策略def efficient_principal_angles(subspace1, subspace2, k10): 使用随机SVD加速主角度计算 from sklearn.utils.extmath import randomized_svd M subspace1.T subspace2 # 只计算前k个奇异值 U, s, Vh randomized_svd(M, n_componentsk) cos_angles s cos_angles np.clip(cos_angles, -1.0, 1.0) angles np.arccos(cos_angles) return angles # 使用随机SVD加速计算 angles_fast efficient_principal_angles(subspace1, subspace2, k5)5.2 子空间维度选择子空间维度的选择对结果影响很大。常用的选择方法包括方差解释率保留足够的主成分以解释大部分方差(如95%)稳定性分析通过bootstrap采样评估子空间估计的稳定性任务驱动基于下游任务(如分类)性能选择最优维度def select_subspace_dim(data, var_explained0.95): 基于方差解释率选择子空间维度 data_centered data - np.mean(data, axis0) _, s, _ svd(data_centered, full_matricesFalse) explained_variance np.cumsum(s**2) / np.sum(s**2) dim np.argmax(explained_variance var_explained) 1 return dim optimal_dim select_subspace_dim(data1) print(f基于95%方差解释率的最优子空间维度: {optimal_dim})5.3 流形上的统计学习Grassmann流形上的统计操作(如均值计算)需要特殊处理。一个常见方法是使用Karcher均值def grassmann_mean(subspaces, max_iter100, tol1e-6): 计算Grassmann流形上的Karcher均值 # 初始猜测第一个子空间 mean subspaces[0] for _ in range(max_iter): # 计算对数映射(切空间向量) tangents [] for sub in subspaces: M mean.T sub U, s, Vh svd(M, full_matricesFalse) R U Vh tangent sub R.T - mean tangents.append(tangent) # 计算平均切向量 avg_tangent np.mean(tangents, axis0) # 检查收敛 if np.linalg.norm(avg_tangent) tol: break # 指数映射更新均值 U, s, Vh svd(avg_tangent, full_matricesFalse) mean mean Vh.T np.diag(np.cos(s)) Vh U np.diag(np.sin(s)) Vh return mean # 示例计算多个子空间的均值 subspace_list [extract_subspace(np.random.randn(100,50), 10) for _ in range(5)] mean_subspace grassmann_mean(subspace_list)