高维数据相似度计算实战：用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, axis=0) # 计算SVD _, s, vh = svd(data_centered, full_matrices=False) # 选择前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_matrices=False) # 奇异值就是主角度的余弦值 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(f"Binet-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_matrices=False) R = U @ Vh dist = np.linalg.norm(subspace1 - subspace2 @ R.T, 'fro') return dist proc_dist = procrustes_metric(subspace1, subspace2) print(f"Procrustes度量距离: {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-Cauchy | 1-∏cos²θ | 主导方向相似性 | 中等 |
| 最大相关 | sinθ₁ | 最相似方向评估 | 低 |
| 最小相关 | sinθₘ | 最不相似方向评估 | 低 |
| Procrustes | ‖Y₁-Y₂R‖F | 基对齐评估 | 高 |
5. 高级技巧与优化建议
5.1 计算效率优化
对于大规模高维数据，直接计算SVD可能很耗时。可以采用以下优化策略：
def efficient_principal_angles(subspace1, subspace2, k=10): """使用随机SVD加速主角度计算""" from sklearn.utils.extmath import randomized_svd M = subspace1.T @ subspace2 # 只计算前k个奇异值 U, s, Vh = randomized_svd(M, n_components=k) 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, k=5)5.2 子空间维度选择
子空间维度的选择对结果影响很大。常用的选择方法包括：
- 方差解释率：保留足够的主成分以解释大部分方差(如95%)
- 稳定性分析：通过bootstrap采样评估子空间估计的稳定性
- 任务驱动：基于下游任务(如分类)性能选择最优维度
def select_subspace_dim(data, var_explained=0.95): """基于方差解释率选择子空间维度""" data_centered = data - np.mean(data, axis=0) _, s, _ = svd(data_centered, full_matrices=False) 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_iter=100, tol=1e-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_matrices=False) R = U @ Vh tangent = sub @ R.T - mean tangents.append(tangent) # 计算平均切向量 avg_tangent = np.mean(tangents, axis=0) # 检查收敛 if np.linalg.norm(avg_tangent) < tol: break # 指数映射更新均值 U, s, Vh = svd(avg_tangent, full_matrices=False) 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)
