高维数据相似度计算实战：用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流形距离的应用。假设我们有两组来自不同领域的图像特征：

1. 源域 ：自然场景图像(如COCO数据集)
2. 目标域 ：医学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 子空间维度选择

子空间维度的选择对结果影响很大。常用的选择方法包括：

1. 方差解释率 ：保留足够的主成分以解释大部分方差(如95%)
2. 稳定性分析 ：通过bootstrap采样评估子空间估计的稳定性
3. 任务驱动 ：基于下游任务(如分类)性能选择最优维度

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)