HoroPCA学习笔记

1. 背景介绍

本篇博客学习的文章为：HoroPCA: Hyperbolic Dimensionality Reduction via Horospherical Projections。

传统欧氏空间（Euclidean spaces）的主成分分析（Principal Component Analysis，PCA）是一种基本的降维技术，它寻找最能解释原始数据的一系列方向。

下面将PCA推广到一般的双曲空间。

给定一个方向，PCA依赖于下述几个核心点：

1. 这些方向所张成的仿射子空间（affine subspaces）的嵌套序列；
2. 将这些点映射到子空间上的一种投影方法，可以沿着每个方向保留信息；
3. 以方差最大化作为优化目标，选择最优方向。

总体的优化过程为（后面的代码也是基于此思想来进行编写）：给定一个数据集，选择一系列方向，使这些方向张成子空间上的投影方差最大化，从而产生的方向序列最优地解释了数据。比较关键的一点是，算法只依赖于仿射子空间的方向，而不依赖于仿射子空间在空间中的位置。这点可以通过下图(a)看出。

图(a)中黑色的点表示真实空间中的点，红色箭头表示仿射子空间的方向，绿色与蓝色点分别表示在不同仿射空间上的投影点（两个仿射子空间是同方向的）。

目前双曲空间上PCA的几种著名推广分别为：

1. 主测地线分析（Principal Geodesic Analysis, PGA）;
2. 重心子空间分析（Barycentric Subspace Analysis, BSA）.

其中，测地线（geodesic）是指：双曲空间中的最短路径。在Poincaré模型中，它们由穿过原点的直线段和垂直于单位球边界的圆弧表示。

PGA和BSA都使用最近点或测地线投影将点映射到子流形上，这两种方法均不是沿主方向保存信息，例如，它们不能等距地保持任意点之间的双曲距离，并以指数因子缩小所有路径长度。这两种方法只是寻找子空间，而不是在寻找解释数据的方向。这两种方法假设所有最优子空间都经过一个选定的基点，但与欧氏空间的设置不同，这个假设不能保证原始数据在转换后的子流形上保持投影之间的等距（图(b)，绿色点所在空间与蓝色点所在空间中，内部点间的距离在两个子空间中不一致）。

基于此，文章提出了一种只依赖于子空间方向，不依赖其位置的方法（location-independence property）。沿着正交方向(极限球面，horospheres)，也就是图©中的红色虚线，平移目标子流形进行投影，这样可以保持固定的投影距离，从而保证位置独立性。

---

文章中的双曲空间考虑的是 d d d维Poincaré模型，对应的空间为 H d = { x ∈ R d : ∥ x ∥ < 1 } \mathbb{H}^{d}=\left\{x \in \mathbb{R}^{d}:\|x\|<1\right\} Hd={x∈Rd:∥x∥<1}。为了说明双曲空间中的PCA，下面引入理想点（ideal points）的概念。

1）理想点（Ideal points）

与欧氏空间中的平行射线一样， H d \mathbb{H}^{d} Hd空间中相互靠近的测地线射线可以被视为在无穷远处共享一个共同的端点，也称为理想点。直观地说，理想点代表了 H d \mathbb{H}^{d} Hd中的点可以向无穷远处移动的方向。下图则为一个理想点，黑色的线表示双曲测地线，每一个红色的圈表示极限球面。每两个极限球面之间的测地线是等长的。

2）Busemann 坐标

在欧氏空间中，每个方向都可以用一个单位向量 w w w来表示。点 x x x在 w w w方向上投影的长度就是点乘运算： w ⋅ x w \cdot x w⋅x。

在双曲几何中，方向可以用理想点来表示，但点乘并没有明确的定义。

我们可以采用基于射线（ray）的视角来看这个问题。欧氏空间中，如果我们从原点向方向 w w w发射一条光线，坐标 w ⋅ x w \cdot x w⋅x可以被视为该光线方向上到无穷远的标准化距离。换句话说，作为点 y = t w , ( t > 0 ) y = tw, (t>0) y=tw,(t>0)沿着方向 w w w向无穷远移动：

w ⋅ x = lim ⁡ t → ∞ ( d ( 0 , t w ) − d ( x , t w ) ) w \cdot x=\lim _{t \rightarrow \infty}(d(0, t w)-d(x, t w)) w⋅x=t→∞lim​(d(0,tw)−d(x,tw))

这种方法可以推广到其他几何空间，给定一个单位速度的测地线射线 γ ( t ) \gamma(t) γ(t)， γ \gamma γ的Busemann函数 B γ ( x ) B_{\gamma}(x) Bγ​(x)定义为：

B γ ( x ) = lim ⁡ t → ∞ ( d ( x , γ ( t ) ) − t ) B_{\gamma}(x)=\lim _{t \rightarrow \infty}(d(x, \gamma(t))-t) Bγ​(x)=t→∞lim​(d(x,γ(t))−t)

这个函数只依赖于测地线射线无穷远处的端点，而不是起始点 γ ( 0 ) \gamma(0) γ(0)。因此，给定一个理想点 p p p（方向），我们将 p p p的Busemann函数 B p ( x ) B_{p}(x) Bp​(x)定义为从单位球模型原点开始，端点为 p p p的测地线射线的Busemann函数。直观上，它表示 x x x在 p p p方向上的坐标。在Poincaré模型中，有显示表达式:

B p ( x ) = ln ⁡ ∥ p − x ∥ 2 1 − ∥ x ∥ 2 B_{p}(x)=\ln \frac{\|p-x\|^{2}}{1-\|x\|^{2}} Bp​(x)=ln1−∥x∥2∥p−x∥2​

对应的代码定义：

def busemann(x, p, keepdim=True):
    xnorm = x.norm(dim=-1, p=2, keepdim=True)
    pnorm = p.norm(dim=-1, p=2, keepdim=True)
    p = p / pnorm.clamp_min(MIN_NORM) # clamp_min(MIN_NORM) 表示最小值阶段
    num = torch.norm(p - x, dim=-1, keepdim=True) ** 2
    den = (1 - xnorm ** 2).clamp_min(MIN_NORM) 
    ans = torch.log((num / den).clamp_min(MIN_NORM))
    if not keepdim:
        ans = ans.squeeze(-1)
    return ans

总结：

	欧式空间	双曲空间
成分	单位向量 w w w	理想点 p p p
坐标	内积 x ⋅ w x \cdot w x⋅w	Busemann函数 B p ( x ) B_p(x) Bp​(x)

3）测地线投影

PCA使用的是正交投影将数据投影到子空间。在双曲空间中，给定目标子流形 M M M，空间中的每个点 x x x都被映射到 M M M中与之最近的点:

π M G ( x ) = argmin ⁡ y ∈ M d M ( x , y ) \pi_{M}^{\mathrm{G}}(x)=\underset{y \in M}{\operatorname{argmin}} d_{M}(x, y) πMG​(x)=y∈Margmin​dM​(x,y)

π M G ( x ) \pi_{M}^{\mathrm{G}}(x) πMG​(x)为 x x x投影到的子流形上的点。

4）双曲标记 (Hyperbolic Flags)

flag的含义就是一个一个递归嵌套的线性子空间序列，给定一系列点 S S S（可以在空间 H d \mathbb{H}^{d} Hd内，也可以在边界球 S ∞ d − 1 \mathbb{S}_{\infty}^{d-1} S∞d−1​内），包含所有 S S S的 H d \mathbb{H}^{d} Hd中的最小测地线子流形称为 S S S的geodesic hull，用 G H ( S ) \mathrm{GH}(S) GH(S)表示。

因此，给定 K K K个理想点 p 1 , p 2 , … , p K p_{1}, p_{2}, \ldots, p_{K} p1​,p2​,…,pK​以及一个基点 b ∈ H d b \in \mathbb{H}^{d} b∈Hd，我们可以得到一系列嵌套的测地线子流形 GH ⁡ ( b , p 1 ) ⊂ G H ( b , p 1 , p 2 ) ⊂ ⋯ ⊂ G H ( b , p 1 , … , p K ) \operatorname{GH}\left(b, p_{1}\right) \subset \mathrm{GH}\left(b, p_{1}, p_{2}\right) \subset \cdots \subset \mathrm{GH}\left(b, p_{1}, \ldots, p_{K}\right) GH(b,p1​)⊂GH(b,p1​,p2​)⊂⋯⊂GH(b,p1​,…,pK​)。

这里的基点 b ∈ H d b \in \mathbb{H}^{d} b∈Hd其实就是欧氏空间中的原点 o o o，可以随意选择一个点，对我们最后的运算结果不影响。同时我们也我们假设 b , p 1 , … , p K b, p_{1}, \ldots, p_{K} b,p1​,…,pK​都不在其他 K − 1 K - 1 K−1个点的测地壳中（例如： p k p_{k} pk​不在 G H ( b , p 1 , … , p k − 1 ) \mathrm{GH}\left(b, p_{1}, \ldots, p_{k-1}\right) GH(b,p1​,…,pk−1​)中）。这个类似于欧氏空间中的线性无关。

5）通过极限球面投影（Projections via Horospheres）

a） K = 1 K=1 K=1

对于 K = 1 K=1 K=1, 有一个理想点 p p p 与一个基点 b b b，geodesic hull G H ( b , p ) \mathrm{GH}(b, p) GH(b,p) 为测地线 γ \gamma γ。 我们的目标是将 π b , p H ( x ) \pi_{b, p}^{\mathrm{H}}(x) πb,pH​(x) 中的每一个 x ∈ H d x \in \mathbb{H}^{d} x∈Hd 映射到 γ \gamma γ 在方向 p p p 上具有相同Busemann坐标 :
B p ( x ) = B p ( π b , p H ( x ) ) B_{p}(x)=B_{p}\left(\pi_{b, p}^{\mathrm{H}}(x)\right) Bp​(x)=Bp​(πb,pH​(x))
由方向 p p p与点 x x x构成的极限球面记为 S ( p , x ) S(p, x) S(p,x)， π b , p H ( x ) \pi_{b, p}^{\mathrm{H}}(x) πb,pH​(x)为其上的点， p p p为中心，球面穿过 x x x，因此可以定义：
π b , p H ( x ) : = γ ∩ S ( p , x ) \pi_{b, p}^{\mathrm{H}}(x):=\gamma \cap S(p, x) πb,pH​(x):=γ∩S(p,x)

对任意 x ∈ H d x \in \mathbb{H}^{d} x∈Hd，若 y ∈ G H ( x , p ) y \in \mathrm{GH}(x, p) y∈GH(x,p)，则有：
d H ( π b , p H ( x ) , π b , p H ( y ) ) = d H ( x , y ) d_{\mathbb{H}}\left(\pi_{b, p}^{\mathrm{H}}(x), \pi_{b, p}^{\mathrm{H}}(y)\right)=d_{\mathbb{H}}(x, y) dH​(πb,pH​(x),πb,pH​(y))=dH​(x,y)

下图中的 O O O即为前面定义的基点 b b b， x ′ x^\prime x′为 π b , p H ( x ) \pi_{b, p}^{\mathrm{H}}(x) πb,pH​(x)，也就是 x x x沿着极限球面（绿色圈）在测地线 γ \gamma γ上的投影。

b） K > 1 K>1 K>1

下面考虑在更高维子流形上的投影。

首先固定一个基点 b ∈ H d b \in \mathbb{H}^{d} b∈Hd与 K > 1 K>1 K>1个理想点 { p 1 , … , p K } \left\{p_{1}, \ldots, p_{K}\right\} {p1​,…,pK​}。下面需要定义一个从 H d \mathbb{H}^{d} Hd 到 M : = G H ( b , p 1 , … , p K ) M:=\mathrm{GH}\left(b, p_{1}, \ldots, p_{K}\right) M:=GH(b,p1​,…,pK​)的映射，从而在 p 1 , … , p K p_{1}, \ldots, p_{K} p1​,…,pK​每个方向上均有一个对应的：Busemann 坐标值：
B p j ( x ) = B p j ( π b , p 1 , … , p K H ( x ) )  for every  j = 1 , … , K B_{p_{j}}(x)=B_{p_{j}}\left(\pi_{b, p_{1}, \ldots, p_{K}}^{\mathrm{H}}(x)\right) \text { for every } j=1, \ldots, K Bpj​​(x)=Bpj​​(πb,p1​,…,pK​H​(x)) for every j=1,…,K
类似 K = 1 K=1 K=1的情形，也能类似定义投影点：
π b , p 1 , … , p K H : H d → M x ↦ M ∩ S ( p 1 , x ) ∩ ⋯ ∩ S ( p K , x ) \begin{aligned} \pi_{b, p_{1}, \ldots, p_{K}}^{\mathrm{H}}: \mathbb{H}^{d} & \rightarrow M \\ & x \mapsto M \cap S\left(p_{1}, x\right) \cap \cdots \cap S\left(p_{K}, x\right) \end{aligned} πb,p1​,…,pK​H​:Hd​→Mx↦M∩S(p1​,x)∩⋯∩S(pK​,x)​
需要注意的是，这里的投影点通常可能会有两个，我们选择离基点更近的一个定义为最终的投影点 π b , p 1 , … , p K H \pi_{b, p_{1}, \ldots, p_{K}}^{\mathrm{H}} πb,p1​,…,pK​H​。

下面以一个 d = 3 , K = 2 d=3, K=2 d=3,K=2的例子，来说明高维情形的投影。

Horospherical投影 π b , p 1 , p 2 H \pi_{b, p_{1}, p_{2}}^{\mathrm{H}} πb,p1​,p2​H​ 是从 H 3 \mathbb{H}^{3} H3 空间投影到2维子流形中得到的。 p 1 , p 2 p_{1}, p_{2} p1​,p2​ 为理想点（即方向），基点 b ∈ H 3 b \in \mathbb{H}^{3} b∈H3 为Poincaré 球的中心。两个理想点与基点构成的geodesic hull G H ( b , p 1 , p 2 ) \mathrm{GH}\left(b, p_{1}, p_{2}\right) GH(b,p1​,p2​) 为黄色圆表示的双曲平面（圆内区域构成的平面均为此区域）。 p 1 p_{1} p1​ 与 p 2 p_{2} p2​ 之间的测地线 P = GH ⁡ ( p 1 , p 2 ) P=\operatorname{GH}\left(p_{1}, p_{2}\right) P=GH(p1​,p2​) 标记为蓝色，称之为投影的脊（spine）。对任意输入点 x ∈ H 3 x \in \mathbb{H}^{3} x∈H3，我们可以画出两个极限球面 S ( p 1 , x ) S\left(p_{1}, x\right) S(p1​,x) 与 S ( p 2 , x ) S\left(p_{2}, x\right) S(p2​,x) ，由红色线条表示（两个红色圆表示两个球面，需要立体看待）。两个极限球面相交与一个圆 S ( x ) S(x) S(x) ，由绿色线条表示。注意，此圆上的集合均满足： { y ∈ H 3 : B p j ( y ) = B p j ( x ) , j = 1 , 2 } \left\{y \in \mathbb{H}^{3}: B_{p_{j}}(y)=B_{p_{j}}(x), j=1,2\right\} {y∈H3:Bpj​​(y)=Bpj​​(x),j=1,2}，即投影到 geodesic hull G H ( b , p 1 , p 2 ) \mathrm{GH}\left(b, p_{1}, p_{2}\right) GH(b,p1​,p2​) 上的坐标值均相等，且此圆关于脊 P P P对称。绿色圆环与黄色圆面 G H ( b , p 1 , p 2 ) \mathrm{GH}\left(b, p_{1}, p_{2}\right) GH(b,p1​,p2​) 相交于两点 x ′ x^{\prime} x′ 与 x ′ ′ x^{\prime \prime} x′′，他们同样关于脊对称。由于 x ′ x^{\prime} x′ 距离基点 b b b更近，因此定义 π b , p 1 , p 2 H ( x ) \pi_{b, p_{1}, p_{2}}^{\mathrm{H}}(x) πb,p1​,p2​H​(x) 为 x ′ x^{\prime} x′。

对应的投影代码如下：

def _project(self, x, Q):
    if self.n_components == 1:
        proj = project_kd(Q, x)[0]
    else:
        if self.hyperboloid:
            hyperboloid_ideals = hyperboloid.from_poincare(Q, ideal=True)
            hyperboloid_x = hyperboloid.from_poincare(x)
            hyperboloid_proj = hyperboloid.horo_projection(hyperboloid_ideals, hyperboloid_x)[0]
            proj = hyperboloid.to_poincare(hyperboloid_proj)
        else:
            proj = project_kd(Q, x)[0]
    return proj

若hyperboloid=True，则会将目前的庞加莱圆盘模型转化到双曲空间中（理想点与数据点分别进行转换）。而后在双曲空间中进行投影，最后再将投影的结果映射回原本的庞加莱圆盘中。

而当hyperboloid=False时，则直接利用庞加莱圆盘中，对前面描述的投影方式进行计算。下面代码为d维投影到k维。

def project_kd(p, x, keep_ambient=True):
    """ Project n points in dimension d onto 'direction' spanned by k ideal points
    p: (..., k, d) ideal points
    x: (..., n, d) points to project

    Returns:
    projection_1: (..., n, s) where s = d if keep_ambient==True otherwise s = k
    projection_2: same as projection_1. this is guaranteed to be the ideal point in the case k = 1
    p: the ideal points
    """

    if len(p.shape) < 2:
        p = p.unsqueeze(0)
    if len(x.shape) < 2:
        x = x.unsqueeze(0)
    k = p.size(-2)
    d = x.size(-1)
    assert d == p.size(-1)
    busemann_distances = busemann(x.unsqueeze(-2), p.unsqueeze(-3), keepdim=False)  # (..., n, k)
    c, r = busemann_to_horocycle(p.unsqueeze(-3), busemann_distances)  # (..., n, k, d) (..., n, k)
    c, r, ortho = sphere_intersections(c, r)  # (..., n, d) (..., n) (..., n, d, k-1)
    # we are looking for a vector spanned by the k ideal points, orthogonal to k-1 given vectors
    # i.e. x @ p @ ortho = 0
    if ortho is None:
        direction = torch.ones_like(busemann_distances)  # (..., n, k)
    else:
        a = torch.matmul(p.unsqueeze(-3), ortho)  # (..., n, k, k-1) = (..., n, k, d) @ (..., n, d, k-1)
        u, s, v = torch.svd(a, some=False)  # a = u s v^T
        direction = u[..., -1]  # (..., n, k)
    direction = direction @ p  # (..., n, d)
    direction = direction / torch.norm(direction, dim=-1, keepdim=True).clamp_min(MIN_NORM)

    projection_1 = c - r.unsqueeze(-1) * direction
    projection_2 = c + r.unsqueeze(-1) * direction
    if not keep_ambient:
        _, _, v = torch.svd(p, some=False)  # P = USV^T => PV = US so last d-k columns of PV are 0
        projection_1 = (projection_1 @ v)[..., :k]
        projection_2 = (projection_2 @ v)[..., :k]
        p = (p @ v)[..., :k]

    return projection_1, projection_2, p

hyperboloid=False时，总体的代码思想是和前面高维投影的想法是一致的，首先根据样本点x与方向坐标p，将坐标转换为busemann坐标，而后再对每个方向作极限球。得到极限球后，针对每个极限球的交集区域，计算中心点与半径。最后得到所有样本投影的坐标。

6）方差

在进行PCA分析时，需要选择合适的方向，从而能够最大限度地保留了原始数据中的信息。欧氏空间中的衡量指标是使投影数据的方差尽可能大，而在双曲几何中，可以类比定义一个衡量方差的指标，其仅依赖于数据点之间的距离。

根据欧氏空间中方差的定义与距离的关系，

1 n 2 ∑ x , y ∈ S ∥ x − y ∥ 2 = 1 n 2 ∑ x , y ∈ S ∥ ( x − μ ( S ) ) − ( y − μ ( S ) ) ∥ 2 = 2 n 2 ∑ x , y ∈ S ∥ x − μ ( S ) ∥ 2 − 2 n 2 ∑ x , y ∈ S < x − μ ( S ) , y − μ ( S ) > = 2 n ∑ x ∈ S ∥ x − μ ( S ) ∥ 2 − 2 n 2 < ∑ x ∈ S ( x − μ ( S ) ) , ∑ y ∈ S ( y − μ ( S ) ) > = 2 n ∑ x ∈ S ∥ x − μ ( S ) ∥ 2 = 2 σ 2 ( S ) \begin{aligned} \frac{1}{n^{2}} \sum_{x, y \in S}\|x-y\|^{2}& = \frac{1}{n^{2}} \sum_{x, y \in S}\|(x-\mu(S)) - (y-\mu(S))\|^{2}\\ & =\frac{2}{n^{2}} \sum_{x, y \in S}\|x-\mu(S) \|^{2} - \frac{2}{n^{2}} \sum_{x, y \in S} \left<x-\mu(S), y-\mu(S)\right> \\ & =\frac{2}{n} \sum_{x \in S}\|x-\mu(S) \|^{2} - \frac{2}{n^{2}}\left< \sum_{x \in S} (x-\mu(S)) , \sum_{y \in S} (y-\mu(S))\right> \\ & =\frac{2}{n} \sum_{x \in S}\|x-\mu(S) \|^{2} \\ & =2 \sigma^{2}(S)\\ \end{aligned} n21​x,y∈S∑​∥x−y∥2​=n21​x,y∈S∑​∥(x−μ(S))−(y−μ(S))∥2=n22​x,y∈S∑​∥x−μ(S)∥2−n22​x,y∈S∑​⟨x−μ(S),y−μ(S)⟩=n2​x∈S∑​∥x−μ(S)∥2−n22​⟨x∈S∑​(x−μ(S)),y∈S∑​(y−μ(S))⟩=n2​x∈S∑​∥x−μ(S)∥2=2σ2(S)​

类比到其他空间，这里定义一个广义的方差：

σ H 2 ( S ) = 1 n 2 ∑ x , y ∈ S d H ( x , y ) 2 \sigma_{\mathbb{H}}^{2}(S)=\frac{1}{n^{2}} \sum_{x, y \in S} d_{\mathbb{H}}(x, y)^{2} σH2​(S)=n21​x,y∈S∑​dH​(x,y)2

上述定义的直观想法就是使得在投影平面上，样本与样本之间距离尽可能的大。代码中方差的计算如下（frechet_variance默认为False）：

def compute_variance(self, x):
    """ x are projected points. """
    if self.frechet_variance:
        Q = [self.mean_weights[i] * self.components[i] for i in range(self.n_components)]  # (k, d)
        mean = sum(Q).squeeze(0)
        distances = poincare.distance(mean, x)
        var = torch.mean(distances ** 2)
    else:
        distances = poincare.pairwise_distance(x)
        var = torch.mean(distances ** 2)
    return var

其中，距离的定义如下：

def distance(x, y, keepdim=True):
    """Hyperbolic distance on the Poincare ball with curvature c.
    Args:
        x: torch.Tensor of size B x d with hyperbolic points
        y: torch.Tensor of size B x d with hyperbolic points
    Returns: torch,Tensor with hyperbolic distances, size B x 1
    """
    pairwise_norm = mobius_add(-x, y).norm(dim=-1, p=2, keepdim=True)
    dist = 2.0 * torch.atanh(pairwise_norm.clamp(-1 + MIN_NORM, 1 - MIN_NORM))
    if not keepdim:
        dist = dist.squeeze(-1)
    return dist

def pairwise_distance(x, keepdim=False):
    """All pairs of hyperbolic distances (NxN matrix)."""
    return distance(x.unsqueeze(-2), x.unsqueeze(-3), keepdim=keepdim)

代码中，unsqueeze(-2)表示在倒数第二个维度上增加一个维度，如：（5，6）变为（5，1，6）；squeeze(-1)则表示将倒数第一个维度去掉（但只有维度为1的时候才能去掉）。

注意到，这里的距离与我们前面定义的距离形式不一致，常见的庞加莱圆盘中距离的定义为：

d ( u , v ) = arccosh ⁡ ( 1 + 2 ∣ ∣ u − v ∣ ∣ 2 ( 1 − ∣ ∣ u ∣ ∣ 2 ) ( 1 − ∣ ∣ v ∣ ∣ 2 ) ) . d(u,v)=\operatorname {arccosh} (1+2{\frac {||u-v||^{2}}{(1-||u||^{2})(1-||v||^{2})}}). d(u,v)=arccosh(1+2(1−∣∣u∣∣2)(1−∣∣v∣∣2)∣∣u−v∣∣2​).

而这里使用的是一种更广义的庞加莱度量，不仅仅局限于庞加莱圆盘上，在半径为1的庞加莱圆盘下，两者的计算其实是等价的。庞加莱度量定义如下。

首先定义单位圆盘：
U = { z = x + i y : ∣ z ∣ = ( x 2 + y 2 ) ≤ 1 } U=\{z=x+iy:|z|={\sqrt {(x^{2}+y^{2})}}\leq 1\} U={z=x+iy:∣z∣=(x2+y2) ​≤1}

对 z 1 , z 2 ∈ U z_{1},z_{2}\in U z1​,z2​∈U 的庞加莱度量为：

d ( z 1 , z 2 ) = 2 tanh ⁡ − 1 ∣ z 1 − z 2 1 − z 1 z 2 ‾ ∣ . d(z_{1},z_{2})=2\tanh ^{-1}\left|{\frac {z_{1}-z_{2}}{1-z_{1}{\overline {z_{2}}}}}\right|. d(z1​,z2​)=2tanh−1∣∣∣∣​1−z1​z2​​z1​−z2​​∣∣∣∣​.

由于实际中两个样本正交，因此这里的 z 1 z 2 ‾ = 0 z_{1}{\overline {z_{2}}}=0 z1​z2​​=0，从而得到代码中的计算公式。

此外，需要注意：欧式空间的加和在庞加莱球空间中就是莫比乌斯加和（Möbius addition），同理也有对应的莫比乌斯环乘（Möbius multiplication）等，公式定义如下：

对于半径为 s s s的球，Möbius addition：
u ⊕ M v = ( 1 + 2 s 2 u ⋅ v + 1 s 2 ∥ v ∥ 2 ) u + ( 1 − 1 s 2 ∥ u ∥ 2 ) v 1 + 2 s 2 u ⋅ v + 1 s 4 ∥ u ∥ 2 ∥ v ∥ 2 \mathbf{u} \oplus_{M} \mathbf{v}=\frac{\left(1+\frac{2}{s^{2}} \mathbf{u} \cdot \mathbf{v}+\frac{1}{s^{2}}\|\mathbf{v}\|^{2}\right) \mathbf{u}+\left(1-\frac{1}{s^{2}}\|\mathbf{u}\|^{2}\right) \mathbf{v}}{1+\frac{2}{s^{2}} \mathbf{u} \cdot \mathbf{v}+\frac{1}{s^{4}}\|\mathbf{u}\|^{2}\|\mathbf{v}\|^{2}} u⊕M​v=1+s22​u⋅v+s41​∥u∥2∥v∥2(1+s22​u⋅v+s21​∥v∥2)u+(1−s21​∥u∥2)v​

Möbius multiplication:
r ⊗ M v = s ( 1 + ∥ v ∥ s ) r − ( 1 − ∥ v ∥ s ) r ( 1 + ∥ v ∥ s ) r + ( 1 − ∥ v ∥ s ) r v ∥ v ∥ = s tanh ⁡ ( r tanh ⁡ − 1 ∥ v ∥ s ) v ∥ v ∥ \begin{aligned} r \otimes_{M} \mathbf{v} &=s \frac{\left(1+\frac{\|\mathbf{v}\|}{s}\right)^{r}-\left(1-\frac{\|\mathbf{v}\|}{s}\right)^{r}}{\left(1+\frac{\|\mathbf{v}\|}{s}\right)^{r}+\left(1-\frac{\|\mathbf{v}\|}{s}\right)^{r}} \frac{\mathbf{v}}{\|\mathbf{v}\|} \\ &=s \tanh \left(r \tanh ^{-1} \frac{\|\mathbf{v}\|}{s}\right) \frac{\mathbf{v}}{\|\mathbf{v}\|} \end{aligned} r⊗M​v​=s(1+s∥v∥​)r+(1−s∥v∥​)r(1+s∥v∥​)r−(1−s∥v∥​)r​∥v∥v​=stanh(rtanh−1s∥v∥​)∥v∥v​​

对应的代码如下：

def mobius_add(x, y):
    """Mobius addition."""
    x2 = torch.sum(x * x, dim=-1, keepdim=True)
    y2 = torch.sum(y * y, dim=-1, keepdim=True)
    xy = torch.sum(x * y, dim=-1, keepdim=True)
    num = (1 + 2 * xy + y2) * x + (1 - x2) * y
    denom = 1 + 2 * xy + x2 * y2
    return num / denom.clamp_min(MIN_NORM)

def mobius_mul(x, t):
    """Mobius multiplication."""
    normx = x.norm(dim=-1, p=2, keepdim=True).clamp(min=MIN_NORM, max=1. - 1e-5)
    return torch.tanh(t * torch.atanh(normx)) * x / normx

def midpoint(x, y):
    """Computes hyperbolic midpoint beween x and y."""
    t1 = mobius_add(-x, y)
    t2 = mobius_mul(t1, 0.5)
    return mobius_add(x, t2)

后续方法与代码请参见：HoroPCA: Hyperbolic Dimensionality Reduction via Horospherical Projections 学习笔记——2.方法与代码