1. 项目概述为什么我们需要双曲空间里的统计模型如果你在机器学习领域特别是几何深度学习方向摸索过一阵子大概率听说过“双曲空间”这个词。从词嵌入到社交网络分析再到最近的生物序列建模双曲表示法因其能优雅地捕捉数据的层次或树状结构而备受青睐。简单来说想象一下你要画一棵巨大的家族树或者一个复杂的公司组织架构图。在平坦的欧几里得纸上随着层级和分支的增多图会变得异常拥挤边与边之间交叉混乱。但如果你有一张具有“负曲率”的、像马鞍面一样的纸你就能把越往深处树根或高层的节点画得“越小”而把外围的叶子节点舒展开从而用更少的空间、更自然的几何关系来容纳整个结构。这就是双曲空间的核心魅力——它天生适合表示具有指数增长关系的层次数据。然而理论上的美好愿景在工程落地时总会遇到现实的骨感。过去几年大多数双曲机器学习实验都依赖于一个叫做“双曲空间上的广义正态分布”的模型来编码不确定性。这个模型听起来很通用但实际用起来尤其是在数据维度升高时问题就暴露出来了它的表达能力不足采样效率低下参数估计也不稳定。这就像你有一把号称“万能”的扳手结果拧螺丝滑丝拧螺母又卡不住严重制约了更复杂模型比如双曲变分自编码器、双曲正态化流的性能天花板。这正是我们这次要深入探讨的核心为双曲空间特别是庞加莱圆盘和球模型建立一个更坚实、更实用的统计建模与优化框架。这个框架的基石不是别的正是数学中优美而强大的群论。我们将看到通过双曲空间的等距变换群主要是莫比乌斯变换群我们可以构造出一族性质极佳的莫比乌斯分布。这族分布不仅具有我们梦寐以求的共形不变性即在双曲等距变换下分布的基本形态保持不变还带来了两大实操福音极其简便的随机采样方法和高效的最大似然参数估计。本文将带你从理论到实践完整拆解这个基于群论的双曲机器学习框架手把手展示如何实现双曲空间中的统计建模与优化。2. 理论基石双曲空间、群作用与不变性在跳进公式和代码之前我们必须先打好地基理解双曲空间到底是什么以及“群”在其中扮演了何等关键的角色。这并非单纯的数学炫技而是理解后续所有算法设计“为什么这么做”的根本。2.1 双曲空间的两种常见模型庞加莱与贝尔格曼双曲空间有多种等价的模型最常用的两个是庞加莱圆盘/球模型和贝尔格曼球模型。它们就像地图的不同投影方式本质描述的是同一个几何空间但坐标表示和计算公式不同各有优劣。庞加莱模型是最直观的。对于二维情况就是单位圆盘B^2 {z ∈ C: |z| 1}。两点z1, z2之间的双曲距离公式为ρ(z1, z2) arcosh(1 2 * |z1 - z2|^2 / ((1 - |z1|^2)(1 - |z2|^2)))。 这个模型下双曲直线表现为与单位圆边界垂直的圆弧或直径。它的优点是可视化极好所有点都落在有限区域内。贝尔格曼模型则是复向量空间C^m中的单位球D^m {z ∈ C^m: |z| 1}。其度量更复杂一些与复结构相容。当m1时贝尔格曼球就退化成了庞加莱圆盘。对于偶数维实空间中的双曲球我们可以选择将其视为实向量空间采用庞加莱度量或复向量空间采用贝尔格曼度量。本文的精彩之处在于它证明了对于机器学习中的统计建模这两种选择在计算效率和表示能力上没有本质区别这给了我们根据具体问题选择计算更方便的模型的自由。注意选择哪种模型往往取决于你的数据形式和已有代码库。庞加莱模型在实数运算和可视化上更简单而贝尔格曼模型在涉及复数的场景如某些信号处理中可能更自然。实践中庞加莱模型的使用更为广泛。2.2 核心主角莫比乌斯变换群双曲空间的“刚性”结构由其等距变换群决定。在庞加莱圆盘中这个群就是莫比乌斯变换群G2其元素是形如g(z) (αz β) / (β̄z ᾱ)的映射其中|α|^2 - |β|^2 1。这些变换有三个关键几何性质保角性保持曲线间的夹角不变。保圆性将圆或直线映射为圆或直线。等距性保持双曲距离不变。正是这个群的对称性为我们构建具有不变性的统计模型提供了可能。共形不变性是我们模型的核心要求如果我们有一个定义在双曲空间上的概率分布对它施加任何一个莫比乌斯变换g那么变换后的数据应该服从一个同族但参数不同的分布。这意味着我们的模型不会因为坐标系观测视角的旋转、平移或缩放而失效其统计特性在群作用下是协变的。这极大地增强了模型的鲁棒性和理论美感。3. 莫比乌斯分布族定义、性质与直观理解现在我们迎来本文的第一个重磅贡献莫比乌斯分布族Moeb_d(a, s)。它是为解决双曲空间中缺乏良好统计模型这一问题而量身定制的。3.1 分布密度函数解析在d维庞加莱球B^d中莫比乌斯分布的概率密度函数为p(x; a, s) C_d * [ (1 - |x|^2)(1 - |a|^2) / ρ(a, x)^2 ]^s。 其中a ∈ B^d位置参数即分布的“均值”点。在双曲几何中它被称为共形重心也是该分布的黎曼中心。s d-1集中度参数。s值越大样本点就越紧密地聚集在均值点a周围。C_d归一化常数具体形式为π^(d/2) * Γ(1s-d/2) / Γ(1s-d)确保密度函数在整个球上的积分为1。ρ(a, x)a与x之间的双曲距离函数。这个公式是怎么来的它并非凭空猜想而是源于一个深刻的几何事实我们希望密度函数在莫比乌斯变换g下具有简单的变换规律。具体来说如果我们从原点0处的某个对称分布出发通过一个将原点映射到a点的莫比乌斯变换g_a就能自然地得到以a为中心的分布。公式中的(1 - |g_a(x)|^2)^s正是这一思想的体现经过化简就得到了上面的形式。3.2 关键性质与优势为什么说这个分布族是双曲统计建模的“利器”它拥有以下几个令人心动的性质共形不变性核心对于固定的集中度参数s分布族Moeb_d(a, s)在莫比乌斯变换群G_d作用下是封闭的。即若x ~ Moeb_d(a, s)则对任意g ∈ G_d有g(x) ~ Moeb_d(g(a), s)。这意味着整个分布族在群作用下是“齐性”的极大简化了理论分析。单峰与对称性在双曲度量下该分布是单峰且对称的。其等高线等概率密度线是以a为中心的双曲圆。这为它作为许多模型如潜在变量模型中的先验或噪声分布奠定了基础。参数意义清晰参数a和s具有明确的几何和统计意义位置和集中度不像某些复杂分布的参数那样难以解释。计算友好性正如我们将在下一节看到的它支持高效的随机采样和参数估计。这是其能投入实际应用的关键。实操心得不要被复杂的密度函数吓到。在实现时我们几乎不要直接计算这个密度值除了在最大似然估计中。更重要的是利用其性质来设计采样和优化算法。记住a是“中心点”s控制“胖瘦”这就够了。4. 高效算法一双曲空间中的随机采样一个实用的概率分布必须能方便地从中生成随机样本。莫比乌斯分布的一大优势就在于其采样过程可以分解为简单的步骤。4.1 从中心在原点的分布采样我们首先解决最简单的情况如何从Moeb_d(0, s)采样即均值为原点0的分布。 由于分布的旋转对称性采样可以分解为两步采样方向从d维单位球面S^(d-1)上均匀采样一个单位向量u。这在高维统计和图形学中有标准算法如通过标准化一个d维标准正态随机向量。采样半径采样径向距离r |x|。其累积分布函数P(|x| b)已由论文给出见公式43。我们需要从中反解出半径。对于二维庞加莱圆盘 (d2)情况特别简单。半径r的分布函数为P(|z| b) 1 - (1 - b^2)^(s-1)。因此采样过程为生成一个均匀随机数κ ~ U[0,1]。计算半径r sqrt(1 - (1 - κ)^(1/(s-1)))。生成一个均匀随机角度φ ~ U[0, 2π]。则样本点z r * e^(iφ)。对于更高维度半径分布涉及超几何函数但论文指出由于维度d是整数这些级数会退化为关于b的多项式因此求解方程P(|x| b) κ来得到b*在计算上是可行的。4.2 变换到任意均值点一旦我们能从Moeb_d(0, s)采样要得到来自Moeb_d(a, s)的样本就轻而易举了。这正是共形不变性带来的红利采样y ~ Moeb_d(0, s)。计算x g_a(y)其中g_a是那个将原点0映射到目标点a的莫比乌斯变换。则x ~ Moeb_d(a, s)。在庞加莱圆盘中这个变换有显式公式g_a(z) (a z) / (1 \bar{a} z)。在高维庞加莱球中也有类似的变换公式。这意味着我们只需要实现一次从中心分布的采样然后通过一个确定性的、计算代价很低的变换就能得到任意中心点的样本。import numpy as np import matplotlib.pyplot as plt def sample_moebius_2d(a, s, n_samples): 从二维庞加莱圆盘上的 Moeb2(a, s) 分布采样。 参数: a: 复数均值点 (|a| 1) s: 浮点数集中度参数 (s 1) n_samples: 整数样本数量 返回: samples: 形状为 (n_samples,) 的复数数组 # 1. 从中心在原点的分布采样半径和角度 kappa np.random.rand(n_samples) r np.sqrt(1 - np.power(1 - kappa, 1/(s-1))) phi np.random.rand(n_samples) * 2 * np.pi z0 r * np.exp(1j * phi) # 中心化样本 # 2. 通过莫比乌斯变换移动到目标均值 a # g_a(z) (a z) / (1 conj(a) * z) # 注意处理数值稳定性当分母接近0时 denominator 1 np.conj(a) * z0 # 避免除以零但理论上在单位圆盘内不会发生 samples (a z0) / denominator return samples # 示例可视化不同参数的采样 fig, axes plt.subplots(2, 3, figsize(12, 8)) params [(0, 2), (0, 4), (0, 8), (0.9*np.exp(1j*3*np.pi/4), 2), (0.9*np.exp(1j*3*np.pi/4), 4), (0.9*np.exp(1j*3*np.pi/4), 8)] for ax, (a, s) in zip(axes.flat, params): samples sample_moebius_2d(a, s, 1000) ax.scatter(samples.real, samples.imag, s1, alpha0.6) # 画出单位圆边界 circle plt.Circle((0, 0), 1, fillFalse, colorr, linestyle--, linewidth0.5) ax.add_patch(circle) ax.set_aspect(equal) ax.set_xlim(-1.1, 1.1) ax.set_ylim(-1.1, 1.1) ax.set_title(fa{a:.2f}, s{s}) ax.grid(True, alpha0.3) plt.tight_layout() plt.show()代码说明这段代码实现了二维庞加莱圆盘中莫比乌斯分布的采样。它清晰地展示了两个步骤先采样中心化分布再进行莫比乌斯变换。可视化结果能直观看到参数a中心位置和s集中度对样本分布的影响。5. 高效算法二参数的最大似然估计当我们有一组观测数据{x_1, ..., x_N}假设来自莫比乌斯分布时如何估计参数(a, s)最大似然估计MLE是最自然的方法。幸运的是对于莫比乌斯分布MLE 问题具有非常优雅且可分解的结构。5.1 似然函数与对数似然给定观测数据似然函数为所有样本点概率密度的乘积。代入密度公式并取对数后我们得到对数似然函数忽略常数项log L(a, s) N log(s-1) s * Σ_{i1}^N log[ (1-|a|^2)(1-|x_i|^2) / ρ(a, x_i)^2 ]。 为了简化我们定义函数H(a) - (1/N) * Σ_{i1}^N log[ (1-|a|^2)(1-|x_i|^2) / ρ(a, x_i)^2 ]。 那么最大化log L(a, s)就等价于最大化N log(s-1) - N * s * H(a)。5.2 参数估计的可分解性这里出现了美妙的一幕优化问题可以完全分解为两个独立的子问题估计位置参数a观察公式s只作为H(a)的乘数出现。为了最大化整体似然对于任何正的s我们都必须最小化H(a)。而H(a)的最小值点â正是我们观测数据点{x_i}的共形重心。估计集中度参数s在找到â后将其代入关于s的优化问题变为一个单变量凹函数的最大化问题max_{sd-1} N log(s-1) - N * s * H(â)。通过求导置零我们可以得到s的解析解或一个非常容易数值求解的方程。这种可分解性极大地降低了计算复杂度。我们不需要用复杂的联合优化算法去同时折腾a和s而是可以分而治之。5.3 计算核心双曲梯度下降求共形重心那么如何计算共形重心â即最小化H(a)呢这需要用到双曲空间中的优化。论文中提出了一种基于双曲梯度下降的优雅算法。在欧氏空间中梯度下降是a_{t1} a_t - η * ∇_{Eucl} H(a_t)。在双曲空间中我们需要遵循双曲几何的“弯曲”规则进行更新。庞加莱球中的双曲梯度∇_{hyp}与欧氏梯度∇_{Eucl}的关系为∇_{hyp} H(a) (1 - |a|^2)^2 / 4 * ∇_{Eucl} H(a)。但直接使用这个关系进行更新a_{t1} a_t - η * ∇_{hyp} H(a_t)并不方便因为要确保更新后的点仍在单位球内。论文采用了更几何的方法双曲梯度流。它本质上描述了一个动力系统使得数据点x_i(t)在双曲等距变换h_t下演化x_i(t) h_t(x_i(0))。而所有数据点的共形重心a(t)的演化遵循一个简单的常微分方程ODEda/dt (K/2N) * (1 - |a(t)|^2) * Σ_{i1}^N h_a(x_i(0))。 其中K是负的学习率h_a是特定的莫比乌斯变换。这个 ODE 正是H(a)在双曲度量下的梯度流。我们可以用数值方法如欧拉法来求解它从而迭代地找到â。一个更工程化的简化实现方式是使用指数映射和对数映射import torch def poincare_distance(u, v): 计算庞加莱球中两点间的双曲距离 sqrt torch.sqrt norm_u torch.norm(u, dim-1, keepdimTrue) norm_v torch.norm(v, dim-1, keepdimTrue) diff_norm torch.norm(u - v, dim-1, keepdimTrue) denominator (1 - norm_u**2) * (1 - norm_v**2) return torch.acosh(1 2 * diff_norm**2 / denominator.clamp_min(1e-10)) def exp_map(x, v): 在点x处将切空间向量v指数映射回庞加莱球 norm_v torch.norm(v, dim-1, keepdimTrue).clamp_min(1e-10) norm_x torch.norm(x, dim-1, keepdimTrue) # 避免数值问题 lambda_x 2 / (1 - norm_x**2).clamp_min(1e-10) direction v / norm_v return torch.tanh(lambda_x * norm_v / 2) * direction def log_map(x, y): 将庞加莱球中的点y映射到点x处的切空间 diff y - x norm_diff torch.norm(diff, dim-1, keepdimTrue).clamp_min(1e-10) norm_x torch.norm(x, dim-1, keepdimTrue) lambda_x 2 / (1 - norm_x**2).clamp_min(1e-10) scale torch.atanh(norm_diff) / (lambda_x * norm_diff) return scale * diff def hyperbolic_gradient_descent(data, lr0.01, max_iters1000, tol1e-6): 使用双曲梯度下降法估计莫比乌斯分布的位置参数a共形重心。 参数: data: 张量形状 (N, d)N个d维样本点位于庞加莱单位球内。 lr: 学习率。 max_iters: 最大迭代次数。 tol: 收敛容忍度。 返回: a_est: 估计的共形重心。 N, d data.shape # 初始化a为数据的欧氏均值并投影到球内确保|a| 1 a torch.mean(data, dim0) a_norm torch.norm(a) if a_norm 1: a 0.99 * a / a_norm # 轻微缩放到球内 for i in range(max_iters): # 计算H(a)的欧氏梯度 (简化版基于公式推导) # ∇H(a) ∝ Σ_i [ (a - x_i) / (1 - 2a, x_i |a|^2|xi|^2) - ... ] # 这里使用一个更稳定的数值实现计算对数似然关于a的梯度 # 为了稳定我们使用自动微分 a.requires_grad_(True) # 计算 H(a) -1/N * Σ_i log( (1-|a|^2)(1-|x_i|^2) / ρ(a, x_i)^2 ) # 忽略与a无关的项我们最小化 F(a) -1/N * Σ_i [ log(1-|a|^2) - log(ρ(a, x_i)^2) ] term1 torch.log(1 - torch.norm(a)**2) dist_sq poincare_distance(a.unsqueeze(0), data)**2 # (N, 1) term2 torch.log(dist_sq) F - (term1 - torch.mean(term2)) # 计算梯度 F.backward() grad_eucl a.grad a.requires_grad_(False) # 将欧氏梯度转换为双曲梯度: ∇_hyp (1-|a|^2)^2 / 4 * ∇_eucl lambda_a (1 - torch.norm(a)**2) grad_hyp (lambda_a**2 / 4) * grad_eucl # 在双曲空间中沿负梯度方向更新a_new exp_a(-lr * grad_hyp) # 首先将梯度向量移动到a点的切空间它已经在切空间了因为是通过对数/指数映射定义的 # 然后沿该方向进行指数映射 update_vec -lr * grad_hyp a_new exp_map(a, update_vec) # 检查收敛 if torch.norm(a_new - a) tol: a a_new print(f在第 {i1} 次迭代收敛。) break a a_new return a # 示例生成一些样本数据并估计其参数 true_a torch.tensor([0.5, 0.2]) true_s 4.0 # 这里需要调用之前定义的采样函数假设已适配PyTorch # samples sample_moebius_nd(true_a, true_s, 1000) # 为演示我们使用一个简单的替代在true_a附近生成一些扰动点 torch.manual_seed(42) data true_a 0.1 * torch.randn(1000, 2) # 确保所有点都在单位圆盘内 norms torch.norm(data, dim1) if (norms 1).any(): data data / (norms.unsqueeze(1) 1e-5) * 0.99 estimated_a hyperbolic_gradient_descent(data, lr0.1, max_iters500) print(f真实 a: {true_a}) print(f估计 a: {estimated_a})代码说明这是一个简化的双曲梯度下降实现。在实际应用中计算H(a)的解析梯度可能更高效稳定。关键在于使用指数映射exp_map来确保更新后的点始终位于双曲空间庞加莱球内。注意对于大规模数据需要采用随机梯度下降SGD或小批量处理。一旦得到â计算ŝ就非常直接了。对于二维圆盘有解析解ŝ 1 N / Σ_i log[ (1-|â|^2)(1-|x_i|^2) / |1 - \bar{â} x_i|^2 ]。对于高维情况需要解一个关于s的简单方程通常可以用牛顿法快速求解。6. 从理论到实践应用场景与扩展理解了莫比乌斯分布和其参数估计方法我们来看看它能用在哪些地方以及如何扩展。6.1 主要应用场景双曲变分自编码器HVAE的潜在空间先验在标准的 VAE 中潜在空间通常使用标准正态分布。在双曲 VAE 中我们可以用莫比乌斯分布Moeb_d(0, s)作为先验分布它能更好地匹配层次化数据的几何特性。编码器输出分布参数(a, s)解码器从该分布采样并进行重构。双曲嵌入的不确定性量化当我们将节点、单词等实体嵌入到双曲空间时一个点估计如通过梯度下降学到的嵌入向量可能不足以表达置信度。我们可以为每个实体学习一个莫比乌斯分布其中均值a表示嵌入的中心位置集中度s表示不确定性s越大越确定。双曲混合模型单个莫比乌斯分布是单峰的。为了对多模态的双曲数据进行密度估计我们可以使用莫比乌斯分布的混合模型。这类似于高斯混合模型GMM但基分布换成了莫比乌斯分布。期望最大化EM算法可以适配过来其中 E 步计算后验概率M 步就用我们上面介绍的最大似然估计来更新每个分量的参数(a_k, s_k)和混合权重。6.2 模型局限性与扩展方向莫比乌斯分布并非万能。它的主要局限在于其对称性和单峰性。现实世界的数据分布可能是偏斜的、多峰的或有更复杂的结构。对此有几个自然的扩展方向混合模型如上所述这是增加模型表达能力最直接的方法。更复杂的参数化可以考虑让集中度参数s也依赖于位置x或者引入偏斜参数。但这可能会破坏共形不变性这一优美性质并增加计算复杂度。基于流的模型在双曲空间上构建归一化流将一个简单的基分布如原点处的莫比乌斯分布通过一系列可逆的双曲变换映射到复杂的目标分布。这能提供极强的表达能力是当前研究的前沿。6.3 实现中的注意事项与调参经验数值稳定性双曲距离和莫比乌斯变换的公式中涉及1 - |x|^2项当点靠近边界|x| → 1时此项趋于零容易导致数值溢出或下溢。在实现时需要对范数进行裁剪如clamp_max(0.999)并使用对数空间计算以避免中间项溢出。梯度下降的学习率双曲梯度下降对学习率非常敏感。由于双曲空间的几何特性越靠近边界局部看起来“越开阔”建议使用自适应优化器如 Riemannian Adam或者采用学习率衰减策略。一个经验法则是初始学习率可以设得比欧氏空间稍大一些。参数初始化对于a的初始化使用数据的欧氏均值并投影到球内是一个不错的起点。对于s的初始化可以考虑基于数据点到初始a的平均双曲距离来启发式地设置一个值距离小则s大。维度诅咒虽然双曲空间本身是为了缓解高维欧氏空间的维度诅咒而引入的但我们的模型在非常高维时采样步骤中涉及的超几何函数计算可能变得繁琐。需要检查是否有针对高维的近似或简化计算。7. 常见问题与排查实录在实际编码和调试基于此框架的模型时你可能会遇到以下典型问题。这里记录了我的排查思路和解决方法。问题现象可能原因排查步骤与解决方案采样点“溢出”单位圆盘1. 采样算法中计算半径r时数值错误导致r1。2. 莫比乌斯变换g_a(z)中分母1 \bar{a}z计算接近零放大误差。1. 检查半径计算公式确保κ在 (0,1) 开区间并验证r sqrt(1 - (1-κ)^(1/(s-1)))在s1时是否小于1。2. 在变换后对点的范数进行强制裁剪x x / (|x| eps) * 0.999。添加微小常数eps防止除零。双曲梯度下降不收敛或震荡1. 学习率过大。2. 梯度计算有误特别是双曲梯度与欧氏梯度的转换。3. 数据点过于靠近边界导致曲率极大优化不稳定。1. 大幅降低学习率如从0.1试到0.001或换用Riemannian Adam。2. 用自动微分如PyTorch的autograd计算H(a)对a的梯度与手推公式对比验证。3. 考虑对输入数据做一个温和的收缩变换如x - 0.9*x将所有点拉离边界。估计的集中度参数ŝ异常大或为负1. 在计算H(â)时出现数值错误如对数项为负。2. 数据点â估计不准导致H(â)计算错误。3. 公式ŝ 1 N / Σ_i [...]中分母接近或小于零。1. 确保在计算log( (1-|a|^2)(1-|x_i|^2) / ρ(a, x_i)^2 )时内部项大于0。双曲距离ρ总是 ≥0但需防止舍入误差导致为零。2. 先仔细检查并可视化â的估计结果看其是否合理地位于数据点中心。3. 为分母添加一个小的正数epsilon防止非正或检查数据是否严重违背了莫比乌斯分布的假设。在高维d10下采样或似然计算速度极慢1. 高维球面均匀采样开销大。2. 高维下双曲距离和梯度计算涉及大量向量运算。3. 超几何函数计算复杂。1. 使用高效的球面采样算法如从标准正态分布采样后归一化。2. 利用线性代数库如NumPy, PyTorch的向量化操作避免循环。3. 对于特定维度d超几何函数2F1可能退化为初等函数查阅数学手册或使用符号计算预先化简。混合模型EM算法中某个分量的s_k迭代中趋于无穷大该分量可能只匹配到了非常少的几个点导致其试图用无限集中的分布去拟合。1. 引入s_k的先验分布如伽马分布进行正则化转为最大后验估计。2. 设置s_k的上界。3. 检查分量初始化避免空簇。可以使用K-Means在双曲空间中的变体进行初始化。最后我想分享一点个人体会。双曲机器学习是一个充满几何直觉与数学美感的领域。将群论的思想引入不仅让模型有了坚实的理论基础和良好的不变性更常常能导出像本文中这样计算高效的算法。莫比乌斯分布族是一个漂亮的起点但它绝不是终点。当你真正开始用它来建模数据时可能会发现其对称性的限制。这时不要急于抛弃它而是可以思考如何以它为基元去构建更复杂的模型例如混合模型或流模型。理解了这个框架的精髓——利用对称性来简化问题——你就能更自如地驾驭双曲空间甚至将其思想迁移到其他具有对称性的流形上。
