嵌入到非欧几里得空间
默认情况下，UMAP将数据嵌入到欧几里得空间中。对于2D可视化，这意味着数据被嵌入到适合散点图的2D平面中。然而，实际上并没有任何主要的限制阻止算法与其他更有趣的嵌入空间一起工作。在本教程中，我们将探讨如何让UMAP嵌入到其他空间中，如何嵌入到您自己的自定义空间中，以及这种方法可能为何有用。
首先，我们将加载常用的库。在这种情况下，我们不会使用umap.plot功能，而是直接使用matplotlib，因为我们将为一些更独特的嵌入空间生成一些自定义的可视化。
import numpy as np
import numba
import sklearn.datasets
import matplotlib.pyplot as plt
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
import umap
%matplotlib inline
sns.set(style='white', rc={'figure.figsize':(10,10)})
作为测试数据集，我们将使用来自sklearn的PenDigits数据集——嵌入到异质空间可能会显著增加计算负担，因此一个相对简单的小数据集将会很有用。
digits = sklearn.datasets.load_digits()
平面嵌入
普通的平面嵌入非常简单——这是UMAP的默认设置。这里我们将再次运行示例，只是为了确保你熟悉这是如何工作的，以及在平面嵌入的简单情况下，PenDigits数据集的UMAP嵌入结果是什么样子的。
plane_mapper = umap.UMAP(random_state=42).fit(digits.data)
plt.scatter(plane_mapper.embedding_.T[0], plane_mapper.embedding_.T[1], c=digits.target, cmap='Spectral')
球形嵌入
sphere_mapper = umap.UMAP(output_metric='haversine', random_state=42).fit(digits.data)
结果是基于球面上的半正矢距离嵌入的pendigits数据。问题在于，如果我们天真地可视化它，那么我们将得到无意义的结果。
plt.scatter(sphere_mapper.embedding_.T[0], sphere_mapper.embedding_.T[1], c=digits.target, cmap='Spectral')
出错的地方在于，在嵌入距离度量下，位于\((0, \pi)\)的点与位于\((0, 3\pi)\)的点的距离为零，因为这会绕赤道一圈。你会注意到，上述图的x轴和y轴的比例远远超出了\((-\pi, \pi)\)和\((0, 2\pi)\)的范围，所以这不是数据的正确表示。然而，我们可以使用简单的公式将这些数据映射到嵌入在三维空间中的球体上。
x = np.sin(sphere_mapper.embedding_[:, 0]) * np.cos(sphere_mapper.embedding_[:, 1])
y = np.sin(sphere_mapper.embedding_[:, 0]) * np.sin(sphere_mapper.embedding_[:, 1])
z = np.cos(sphere_mapper.embedding_[:, 0])
现在 x, y, 和 z 提供了每个嵌入点的三维坐标，这些点位于球体的表面上。我们可以使用 matplotlib 的三维绘图功能来可视化这一点，并看到我们实际上已经将数据嵌入到了球体的表面上。
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=digits.target, cmap='Spectral')
如果你更喜欢二维图，我们可以将这些数据转换为适当范围内的经纬度坐标，并获得相当于球体数据的地图投影。
x = np.arctan2(x, y)
y = -np.arccos(z)
plt.scatter(x, y, c=digits.target.astype(np.int32), cmap='Spectral')
在自定义度量空间上嵌入
如果你有一些其他的自定义度量空间概念，想要将数据嵌入其中，该怎么办？就像UMAP可以支持输入数据的自定义距离度量（只要它们可以用numba编译），output_metric参数也可以接受自定义距离函数。一个需要注意的是，为了支持梯度下降优化，距离函数需要返回距离和距离梯度的向量。后者可能需要用户进行一些微积分计算。第二个需要注意的是，以没有坐标约束的方式参数化嵌入空间是非常有益的——否则梯度下降可能会使一个点超出嵌入空间，导致不良后果。这就是为什么，例如，球体示例中的点只是环绕而不是将坐标限制在适当的范围内。
让我们通过一个例子来构建一个不同空间的距离度量和梯度：一个环面。环面本质上就是一个甜甜圈的外表面。我们可以用x、y坐标来参数化环面，但需要注意的是我们可以“环绕”（类似于球体）。在这样的模型中，距离大多是欧几里得距离，我们只需要检查哪个方向更短——穿过还是环绕——并确保我们考虑到多次环绕的等价性。我们可以编写一个简单的函数来计算这个距离。
@numba.njit(fastmath=True)
def torus_euclidean_grad(x, y, torus_dimensions=(2*np.pi,2*np.pi)):
"""Standard euclidean distance.
..math::
D(x, y) = \sqrt{\sum_i (x_i - y_i)^2}
"""
distance_sqr = 0.0
g = np.zeros_like(x)
for i in range(x.shape[0]):
a = abs(x[i] - y[i])
if 2*a < torus_dimensions[i]:
distance_sqr += a ** 2
g[i] = (x[i] - y[i])
else:
distance_sqr += (torus_dimensions[i]-a) ** 2
g[i] = (x[i] - y[i]) * (a - torus_dimensions[i]) / a
distance = np.sqrt(distance_sqr)
return distance, g/(1e-6 + distance)
请注意，梯度仅源自标准欧几里得梯度，我们只需根据我们计算距离的方式检查方向。我们现在可以直接将该函数插入到output_metric参数中，并最终将数据嵌入到环面上。
torus_mapper = umap.UMAP(output_metric=torus_euclidean_grad, random_state=42).fit(digits.data)
R = 3 # Size of the doughnut circle
r = 1 # Size of the doughnut cross-section
x = (R + r * np.cos(torus_mapper.embedding_[:, 0])) * np.cos(torus_mapper.embedding_[:, 1])
y = (R + r * np.cos(torus_mapper.embedding_[:, 0])) * np.sin(torus_mapper.embedding_[:, 1])
z = r * np.sin(torus_mapper.embedding_[:, 0])
现在我们可以使用matplotlib可视化结果，并看到数据确实已经被适当地嵌入到了一个环面上。
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=digits.target, cmap='Spectral')
ax.set_zlim3d(-3, 3)
ax.view_init(35, 70)
与环面一样，我们可以进行一些几何操作，并将环面展开成一个具有适当边界的平面。
u = np.arctan2(x,y)
v = np.arctan2(np.sqrt(x**2 + y**2) - R, z)
plt.scatter(u, v, c=digits.target, cmap='Spectral')
虽然到目前为止给出的例子可能有一些用途（因为一些数据确实具有适合的周期性或循环结构，我们期望这些结构在球体或环面上能更好地表示），但大多数数据并不真正属于用户可以事先期望存在于奇异流形上的领域。那么，嵌入到其他空间的能力是否有更实际的用途呢？事实证明是有的。一个有趣的例子是考虑由二维高斯分布形成的空间。我们可以通过两个高斯分布之间的PDF的内积的负对数来测量它们之间的距离（因为这是一个很好的闭式解，并且是合理可计算的）。这为我们提供了一个可以嵌入的度量空间，其中样本不是表示为二维中的点，而是表示为二维中的高斯分布，编码了高维空间中每个样本如何嵌入的一些不确定性。
当然，我们仍然有适合SGD的参数化问题——要求协方差矩阵是对称且正定的，这是具有挑战性的。相反，我们可以用宽度、高度和角度来参数化协方差，并在需要时从这些参数中恢复协方差矩阵。这为我们提供了总共5个组件来嵌入（两个用于均值，3个用于描述协方差的参数）。我们可以简单地这样做，因为已经定义了适当的度量。请注意，我们必须特别传递n_components=5，因为我们需要显式地嵌入到一个5维空间，以支持与2d高斯相关的所有协方差参数。
gaussian_mapper = umap.UMAP(output_metric='gaussian_energy',
n_components=5,
random_state=42).fit(digits.data)
由于我们将数据嵌入到了一个五维空间中，可视化不像之前那样简单。我们可以通过仅查看均值来开始可视化结果，这些均值是高斯分布模式的二维位置。传统的散点图足以满足这一需求。
plt.scatter(gaussian_mapper.embedding_.T[0], gaussian_mapper.embedding_.T[1], c=digits.target, cmap='Spectral')
我们看到我们得到了一个类似于标准嵌入到欧几里得空间的结果，但聚类不太清晰，且聚类之间有更多的点。为了更清楚地了解发生了什么，有必要设计一种方法来显示包含在额外3个维度中的一些额外信息，这些信息提供了协方差数据。为此，能够绘制对应于2D高斯PDF的超水平集的椭圆将是有帮助的。我们可以通过编写一个简单的函数来开始，该函数根据位置、宽度、高度和角度在图上绘制椭圆（因为这是嵌入计算数据的格式）。
from matplotlib.patches import Ellipse
def draw_simple_ellipse(position, width, height, angle,
ax=None, from_size=0.1, to_size=0.5, n_ellipses=3,
alpha=0.1, color=None,
**kwargs):
ax = ax or plt.gca()
angle = (angle / np.pi) * 180
width, height = np.sqrt(width), np.sqrt(height)
# Draw the Ellipse
for nsig in np.linspace(from_size, to_size, n_ellipses):
ax.add_patch(Ellipse(position, nsig * width, nsig * height,
angle, alpha=alpha, lw=0, color=color, **kwargs))
现在我们可以通过提供中心的散点图（如前所述）来绘制数据，但将其叠加在相关高斯分布的超水平集椭圆上。明显的缺点是这会导致大量的重叠绘制，但它至少提供了一种开始理解我们生成的嵌入的方法。
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
colors = plt.get_cmap('Spectral')(np.linspace(0, 1, 10))
for i in range(gaussian_mapper.embedding_.shape[0]):
pos = gaussian_mapper.embedding_[i, :2]
draw_simple_ellipse(pos, gaussian_mapper.embedding_[i, 2],
gaussian_mapper.embedding_[i, 3],
gaussian_mapper.embedding_[i, 4],
ax, color=colors[digits.target[i]],
from_size=0.2, to_size=1.0, alpha=0.05)
ax.scatter(gaussian_mapper.embedding_.T[0],
gaussian_mapper.embedding_.T[1],
c=digits.target, cmap='Spectral', s=3)
现在我们可以看到，点的协方差结构在绝对大小和形状上都有很大的变化。我们注意到，许多落在簇之间的点具有更大的方差，在某种意义上代表了嵌入位置的不确定性更大。还值得注意的是，椭圆的形状可以显著变化——有几个非常拉伸的椭圆，与许多非常圆的椭圆截然不同；在某种意义上，这代表了不确定性更多地沿着一条线分布的情况。
虽然这个图突出了一些异常点中的协方差结构，但在实践中，这里的过度绘制掩盖了集群本身中许多更有趣的结构。我们可以尝试通过每个点只绘制一个椭圆并使用较低的alpha通道值来更好地看到这种结构，使它们更加半透明。
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
for i in range(gaussian_mapper.embedding_.shape[0]):
pos = gaussian_mapper.embedding_[i, :2]
draw_simple_ellipse(pos, gaussian_mapper.embedding_[i, 2],
gaussian_mapper.embedding_[i, 3],
gaussian_mapper.embedding_[i, 4],
ax, n_ellipses=1,
color=colors[digits.target[i]],
from_size=1.0, to_size=1.0, alpha=0.01)
ax.scatter(gaussian_mapper.embedding_.T[0],
gaussian_mapper.embedding_.T[1],
c=digits.target, cmap='Spectral', s=3)
这让我们可以看到集群密度相对于协方差结构的变化——一些集群的协方差一直非常紧密，而其他集群则更加分散（因此在某种意义上，具有更大的相关不确定性）。当然，即使在这里，我们仍然存在一定程度的重叠绘图，调整alpha通道以使内容可见将变得越来越困难。相反，我们想要的是一个实际的密度图，显示所有这些高斯分布总和的密度。
为此，我们需要定义一些函数，这些函数的执行将使用numba进行加速：在给定点评估二维高斯分布的密度；评估给定点的密度，该密度是对一组多个高斯分布求和的结果；以及一个函数，用于生成某个网格中每个点的密度（仅对附近的高斯分布求和，以使这种简单方法更易于计算）。
from sklearn.neighbors import KDTree
@numba.njit(fastmath=True)
def eval_gaussian(x, pos=np.array([0, 0]), cov=np.eye(2, dtype=np.float32)):
det = cov[0,0] * cov[1,1] - cov[0,1] * cov[1,0]
if det > 1e-16:
cov_inv = np.array([[cov[1,1], -cov[0,1]], [-cov[1,0], cov[0,0]]]) * 1.0 / det
diff = x - pos
m_dist = cov_inv[0,0] * diff[0]**2 - \
(cov_inv[0,1] + cov_inv[1,0]) * diff[0] * diff[1] + \
cov_inv[1,1] * diff[1]**2
return (np.exp(-0.5 * m_dist)) / (2 * np.pi * np.sqrt(np.abs(det)))
else:
return 0.0
@numba.njit(fastmath=True)
def eval_density_at_point(x, embedding):
result = 0.0
for i in range(embedding.shape[0]):
pos = embedding[i, :2]
t = embedding[i, 4]
U = np.array([[np.cos(t), np.sin(t)], [np.sin(t), -np.cos(t)]])
cov = U @ np.diag(embedding[i, 2:4]) @ U
result += eval_gaussian(x, pos=pos, cov=cov)
return result
def create_density_plot(X, Y, embedding):
Z = np.zeros_like(X)
tree = KDTree(embedding[:, :2])
for i in range(X.shape[0]):
for j in range(X.shape[1]):
nearby_points = embedding[tree.query_radius([[X[i,j],Y[i,j]]], r=2)[0]]
Z[i, j] = eval_density_at_point(np.array([X[i,j],Y[i,j]]), nearby_points)
return Z / Z.sum()
现在我们只需要一个合适的点网格。我们可以使用上面看到的绘图边界，并选择一个适合计算的网格大小。numpy的meshgrid函数可以提供实际的网格。
X, Y = np.meshgrid(np.linspace(-7, 9, 300), np.linspace(-8, 8, 300))
现在我们可以使用上面定义的函数来计算网格中每个点的密度，给定嵌入产生的高斯分布。
Z = create_density_plot(X, Y, gaussian_mapper.embedding_)
现在我们可以使用imshow将结果查看为密度图。
plt.imshow(Z, origin='lower', cmap='Reds', extent=(-7, 9, -8, 8), vmax=0.0005)
plt.colorbar()
额外内容：在双曲空间中的嵌入
作为一个额外的例子，让我们看看如何将数据嵌入到双曲空间中。 最流行的可视化模型是Poincare的圆盘模型。 下面展示了一个在Poincare圆盘模型中对双曲空间进行规则镶嵌的例子；你可能会注意到它与M.C. Escher的著名图像相似。
hyperbolic_mapper = umap.UMAP(output_metric='hyperboloid',
random_state=42).fit(digits.data)
一个简单的可视化选项是简单地查看我们得到的x和y坐标：
plt.scatter(hyperbolic_mapper.embedding_.T[0],
hyperbolic_mapper.embedding_.T[1],
c=digits.target, cmap='Spectral')
我们也可以求解z坐标，并查看位于三维空间中的双曲面上的数据。
x = hyperbolic_mapper.embedding_[:, 0]
y = hyperbolic_mapper.embedding_[:, 1]
z = np.sqrt(1 + np.sum(hyperbolic_mapper.embedding_**2, axis=1))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=digits.target, cmap='Spectral')
ax.view_init(35, 80)
但我们还可以做得更多——既然我们已经成功地将数据嵌入到双曲空间中，我们可以将数据映射到庞加莱圆盘模型中。实际上，这是一个直接的计算。
disk_x = x / (1 + z)
disk_y = y / (1 + z)
现在我们可以按照最初的想法，在庞加莱圆盘模型嵌入中可视化数据。为此，我们只需生成数据的散点图，然后绘制无限远处的边界圆。
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(disk_x, disk_y, c=digits.target, cmap='Spectral')
boundary = plt.Circle((0,0), 1, fc='none', ec='k')
ax.add_artist(boundary)
ax.axis('off');
希望这提供了一个有用的示例，说明如何嵌入到非欧几里得空间中。最后一个示例理想地突出了这种方法的局限性（我们确实需要一个合适的参数化），以及一些潜在的解决方法：我们可以使用替代的参数化进行嵌入，然后将数据转换为所需的表示形式。