填充曲线之美-理解希尔伯特曲线

填充曲线之美：理解希尔伯特曲线

原文：towardsdatascience.com/the-beauty-of-space-filling-curves-understanding-the-hilbert-curve/

0. 简介

填充曲线（SFC）是迷人的数学结构，在数据科学和数据工程中有许多实际应用。虽然它们可能听起来很抽象，但它们通常隐藏在显而易见的地方——例如在 Z 排序或液态聚类（例如在 Databricks 等平台上使用）等术语背后。如果你已经与大规模数据平台合作过，那么你很可能已经在不经意间使用了 SFCs。

尽管在现代系统中具有重要意义，但关于这个主题的信息往往是零散的，这使得理论与实践之间的桥梁难以搭建。本文旨在弥合这一差距，同时聚焦于希尔伯特曲线。

我的目的是提供一个浓缩且易于理解的 SFCs 概述：从它们的数学起源开始，经过实际实施技术，最后到数据处理和优化的实际应用。这不是要取代现有的资源，而是为了提供更详细信息的参考。术语和细节的进一步来源将在全文中引用。

你可能会问：曲线有什么如此迷人的地方？ 毕竟，一条常规曲线很容易理解，可能不是我会首先选择的书籍主题。但是 SFCs 是不同的。它们穿越连续空间中的每一个点，具有分形特性，当在 2D 或 3D 中绘制时会产生视觉上引人注目的图案——尤其是在低迭代次数时。因此，让我们更深入地了解一下。

（如果你想直接从可视化动画开始，请查看我的GitHub 仓库）

1. 填充曲线的历史与理论

SFCs 的研究可以追溯到 19 世纪，当时乔治·康托尔做出了开创性的发现。他表明，“两个有限维光滑流形具有相同的基数，无论它们的维度如何。” [1]

为了说明这一点，考虑单位区间[0, 1] ⊂ R 和单位正方形[0, 1]² ⊂ R²。直观上，人们可能会预期正方形的基数大于线段。然而，康托尔通过交错小数的方法证明了这两个集合实际上具有相同的基数。

这个结果意味着区间和正方形之间存在双射，这意味着它们的元素之间存在一一对应关系。在康托尔的发现之后，一个自然的问题出现了：这些集合之间也存在连续的双射吗？尤根·内托对此问题给出了否定的答案。

在这个背景下，连续性可以几何地解释：一个连续映射将允许一个人在 2D 或 3D 中“绘制”图像，而无需抬起笔——形成一个曲线。这一洞察为后来 SFCs 的发展奠定了基础——这些曲线虽然连续，但可以任意接近填充更高维的空间。

空间。

2. Peano 曲线：空间填充曲线的发现

在 Netto 令人清醒的发现之后，人们开始质疑，如果这样的映射不是双射，它是否可以是满射。第一个能够定义这种映射的人是 G. Peano，他构建了所谓的 Peano 曲线。

Peano 曲线是递归定义的。其定义域是单位区间 [0, 1] ⊂ R，其像位于单位正方形 [0, 1]² ⊂ R²。通过将区间 [0, 1] 重复三等分，并相应地将 R² 中的正方形划分为 3 × 3 的网格，随着迭代次数趋向于无穷大，构造会收敛到实际的空间填充曲线。[1]

图 1：一阶、二阶和三阶 Peano 曲线（从左到右）。

一阶 Peano 曲线的像在更高阶中被复制并镜像。可以观察到，一阶 Peano 曲线的基本模式在更高阶中重新出现，但在每第二次迭代时被镜像。这种交替的镜像和旋转基本元素的过程是其他 SFCs 共有的特征。

(来自维基百科的图片 在公共领域许可下，作者修改)

因此，有限迭代次数的 Peano 曲线图（图 1）并不代表“最终”的 SFC。只有当这个递归映射的迭代次数趋向于无穷大时，实际的 SFC 才会出现，它遍历 [0, 1]² 中的每一个点。从视觉上看，在这个极限下，曲线基本上会呈现出一个从 (0, 0) 到 (1, 1) 的填充正方形。

这个观察引发了一个最初令人反直觉的问题：根据定义，曲线是一维的。虽然它可以嵌入到更高维的空间（n > 1）中，但其内在参数域仍然是一维的。然而，如果 Peano 曲线穿过 [0, 1]² 中的每一个点并因此完全填满平面，它的像是否还可以被视为一维的？答案是：不：Peano 曲线的像有 豪斯多夫维度 2。SFC 的另一个特征是其像有正的 Jordan 内容（Peano-Jordan 度量）。这些事实可能看起来令人惊讶，但它们与分形属性相一致：许多这样的集合具有大于 1 的豪斯多夫维度，甚至有些具有非整数的豪斯多夫维度。

3. 希尔伯特曲线——至今仍受欢迎！

虽然佩亚诺是第一个构造 SFC 的人，但一个更为著名的例子是希尔伯特曲线，由大卫·希尔伯特在 1891 年定义。其定义稍微简单一些，从 2x2 网格开始。像佩亚诺曲线一样，希尔伯特曲线的映射递归地将 [0, 1] 区间内的每个区间以及 [0, 1]² 中的每个正方形在每一步划分为四个更小的区间/正方形。与佩亚诺曲线一样，希尔伯特曲线在迭代次数趋于无穷大时收敛到真正的 SFC。

图 2：左侧的基本单元（阶数为 1）被重复用于构建更高阶的希尔伯特曲线。然而，必要的变换（如镜像和旋转）比佩亚诺曲线的情况更复杂。

(图片由作者提供)

为了本文的目的，我们将重点关注希尔伯特曲线，因为其属性使其成为现代数据平台中的宝贵工具。

3.1 希尔伯特曲线的正式定义

从希尔伯特曲线的定义域区间 [0,1] 开始，每个递归步骤将当前区间划分为四个相等的子区间：a 是左端点，h 是区间宽度，子区间如下：

在 [0,1] 区间内分割区间。（公式来自 [2]，图片由作者提供）

对于 [0, 1] 中任意选择的点，恰好有一个子区间包含该点。然后可以使用相同的规则再次细分该区间，产生一个更细的区间，该区间仍然包含该点。这个过程可以无限进行下去，从而得到曲线上的任意精确位置。相同的递归细分也应用于 [0, 1]²，将每个正方形分割成四个更小的正方形：

在 [0,1]² 中分割象限。（公式来自 [2]，图片由作者提供）

一般属性：

• 满射: 从其递归定义中可以得出，希尔伯特曲线是满射的：在极限情况下，[0, 1]² 中的每个点都被覆盖。嵌套区间是紧致的，相邻区间共享边界点（例如，a + h/4 既是第一个子区间的右端点，也是第二个子区间的左端点）。

  因此整个正方形被填满。然而，映射不是单射——尝试强制双射（例如，通过打开区间）会破坏连续性。

• 连续: 从直观表示中可以清楚地看出这个属性：曲线可以不提笔地绘制。形式上，可以通过证明希尔伯特曲线是连续函数的均匀极限来建立，而均匀收敛保持连续性。

• 不可微: 通过观察希尔伯特曲线的图形，很明显这条曲线是不可微的。

  可微。H.Sagan 使用差商给出了这个属性的证明。

• 局部保持性：与 Z-order 曲线等更简单的映射相比，希尔伯特曲线倾向于保持局部性：在一维参数中靠近的点通常会被映射到附近。这一特性对于大数据平台的应用至关重要。

• 正 Jordan 内容：在无限多次迭代的过程中，希尔伯特曲线的图像具有正 Jordan 测度，这意味着它占据了平面上的非零面积。（佩亚诺-乔丹测度）

• 二维的豪斯多夫维度：相应地，希尔伯特曲线不像通常的一维线那样表现，而是具有Hausdorff 维度 2，这反映了它完全填充了单位正方形。

尽管希尔伯特曲线的早期定义是在二维空间中进行的，但更高维度也是可行的。我们将在下一节讨论的算法可以在任何有限维度中工作。

4 使用斯金林算法计算希尔伯特曲线

希尔伯特曲线的定义是以几何方式给出的，没有给出在给定网格上计算坐标的代数定义，对于 I 中的给定点。在希尔伯特发布他的想法后近 100 年，数学家们才开始考虑如何计算给定希尔伯特索引的点。谁又能责怪他们呢？毕竟，长期以来，没有计算机能够绘制包含数百或数千个点的曲线。在研究过程中，我发现有多种计算希尔伯特曲线的方法——从复数到 L-系统。有些方法非常复杂，而有些则保留了计算曲线单一点的迭代方法。我所寻找的是简单的东西：

• 一个函数，它接受一个希尔伯特索引（即 1D 空间中的任何数字，如 1、2、3）并返回其坐标。你可以将希尔伯特索引视为从左到右的区间编号，对于无限阶的希尔伯特曲线。

• 一个执行逆映射的函数，将坐标映射回其希尔伯特索引。

在互联网上搜索可能的实现方法时，我偶然发现了一个普林斯顿大学的GitHub 仓库，该仓库实现了约翰·斯金林在 2004 年发表的一篇名为编程希尔伯特曲线的论文中的算法。不幸的是，这篇论文对公众并不免费，所以我决定分析普林斯顿仓库中的代码。

4.1 斯金林算法 – 概述

Skilling 观察到，将希尔伯特索引映射到坐标可以用二进制运算优雅地表达。例如，考虑一维中的索引 0，1，2，3。这些对应于 2×2 网格中的坐标（0，0），（1，0），（1，1），（0，1）。在这里，值 0，1，2，3 不再代表单位区间中的分数点（如 1/3），而是离散的区间数。在 2×2 网格中，[0, 1]中有四个区间，以及四个对应的[0, 1]²中的正方形。Skilling 算法将这个想法推广。它使用二进制变换计算从希尔伯特索引到其对应坐标的映射（反之亦然），适用于任何有限维度。基本步骤如下：

1. 将希尔伯特索引从十进制转换为二进制。

2. 将二进制数转换为其格雷码表示。

3. 将格雷码解开成坐标结构。

4. 使用 XOR 操作应用旋转和反射。

5. 将二进制坐标转换回十进制

4.2 二进制表示

为了理解为什么二进制更适合从希尔伯特索引和逆过来计算希尔伯特曲线的点，以下示例可能会有所帮助（我们讨论的是二维空间中的所有内容，但算法适用于任何维度的空间）：

希尔伯特曲线定义在 2×2，4×4，8×8，16×16 等网格上。（记住上面的定义及其递归方法）。

通过观察数字，人们可能会发现间隔的数量随着 2^n 增长，其中 n 是曲线的阶数。这与二进制编码完美匹配：对于一个 n 阶曲线，我们

每个轴需要恰好 n 位来描述网格。

以 4×4 网格（第二阶）为例。每个轴上两个比特就足够了：

1. 第一个比特标识主要象限（左下，左上，右下或右上）。

2. 第二个比特指定该象限内的位置。

例如，希尔伯特索引 2 的二进制形式为 0010。解释如下：

• 00 选择左下象限。

• 10 选择其内部的右上子象限。

图 3：将二进制映射到网格单元。前两个比特编码主要象限，最后两个比特编码

子象限。考虑每个象限中重复出现的模式 00，01，10，11，形成一个希尔伯特曲线。

阶数 1. （图片由作者提供）

然而，如果我们继续对大于 3 的索引进行此过程，我们会遇到一个挑战：曲线的方向从一个象限变为下一个象限。正确处理这些旋转和反射正是格雷码和 XOR 操作（如 Skilling 算法中所示）变得至关重要的地方。

4.3 格雷码表示

Skilling 算法的下一步是从二进制到格雷码的转换。关键区别在于格雷码中，连续的数字只在一个位上不同。这一属性至关重要：它确保曲线可以平滑地从一个个象限移动到下一个（即使曲线在每个象限中的取向仍然不正确）。

通过观察二进制数和曲线不同部分的取向，我们可以看到曲线仍然不正确，但现在每个象限的末端现在连接到下一个象限的起始处。

图 4：将二进制值转换为格雷码后，当前象限的最后一个单元格具有相同的值

作为下一个单元格的第一个单元格（图片由作者提供）

4.4 位的解缠

Skilling 方法的真正“魔法”始于对格雷编码位的重新排序——这一步被称为解缠。在我们的 4×4 示例中，我们最初将四个位解释为(bitx[1], bity[1], bitx[2], bity[2])，其中第一对编码主象限，第二对编码子象限。然而，对于坐标计算，我们需要一种形式为(bitx[1], bitx[2], bity[1], bity[2]*)的结构，以便所有 x 位和 y 位都可以稍后组合成相应的十进制坐标(x, y)。这一步被称为位的解缠。

图 5：格雷码解缠后 4×4 网格中子象限的取向（图片由作者提供）

4.5 纠正变换

在位解缠之后，Skilling 算法的最终一步是在每个象限内旋转和镜像子曲线，以便它们无缝地连接到 n 阶 Hilbert 曲线。

图 6 说明了 4×4 情况下的这个过程。左侧的表格显示了如何通过应用简单的变换（交换和位翻转）将格雷编码坐标转换为标准二进制数。

右侧的图表可视化了这个效果：上象限旋转了 180°，下象限沿对角线镜像，在某些情况下（例如，黄色象限）甚至不需要任何变换。

关键的洞察力在于，经过这些纠正变换后，坐标再次处于标准二进制形式。这意味着 Skilling 算法的输出可以直接转换为(x, y)格式的十进制坐标，无需进一步调整。

图 6：将灰度码转换为二进制坐标的最终转换（图片由作者提供）

斯金林算法的关键转换：输入：以灰度码格式化（bitx[1]，bitx[2]，bity[1]，bity[2]）。在 Python 中，格式将是：[-1，ndims，nbits]。例如，数字 4 将被表示为以下列表/np-array：[[01]，[10]]。对于 x 维度，1 是最不重要的位（LSB），而 0 是最重要的位（MSB）。

（最高有效位）。

1. 从最高有效位（MSB）到最低有效位（LSB）循环。

2. 从最高维度（2D 中的 y）到最低维度的内循环

3. : 查看当前位。如果是 1：翻转维度 0（通常是 x）中所有较低的位。如果是 0：交换当前维度和维度 0（如果它们不同）之间的值。

   在当前维度和维度 0（如果它们不同）的较低位之间交换值。

第 3 步可以通过使用 numpy 的 XOR 操作轻松计算。每次迭代中翻转和交换位的整个过程在以下动画中进行了可视化。

图 7：使用约翰·斯金林算法创建 2D 希尔伯特曲线的过程（图片由作者提供）

图 8：使用约翰·斯金林算法创建 3D 希尔伯特曲线的过程（图片由作者提供）

如果你想更详细地分析算法或简单地生成自己的 2D 或 3D 动画，请查看我的 GitHub 仓库。

5 应用空间填充曲线

在讨论了希尔伯特曲线的理论方面和实现细节之后，问题随之而来，即它可以在哪里应用。在实现过程中，我们看到了如何将希尔伯特索引转换为坐标。对于以下应用，这一过程的逆操作更有趣。

希尔伯特曲线的一个有价值的特点是它将 1D 有序集（即 1，2，3……）映射到 n 维空间中的坐标。它为其遍历的点提供顺序，并且它可以存在于任意大小的向量空间中。因此，希尔伯特曲线用于数据分区和聚类、图像压缩，以及在处理空间数据时用于机器学习中的特征构建。

5.1 使用 SFCs 进行数据分区/聚类

SFCs（空间填充曲线）最突出的应用之一是数据分区。例如，在 Databricks 中，Z 排序基于 Z 曲线，而液体聚类依赖于希尔伯特曲线。原因很简单：

希尔伯特曲线比 Z 曲线更好地保留了局部性，这在索引和分区多维数据时至关重要。在图 9 中，你可以看到一些示例数据点是如何通过将每个点分配给曲线给出的一个分区来映射到希尔伯特曲线点的。

图 9：数据映射到希尔伯特曲线点的过程。红色虚线箭头指示一些映射

例示（图片由作者提供）

当查询应用于数据时（例如，SELECT * FROM table WHERE x in (1,2) and y in (2,3)，所有在这个范围内的点（（1,2），（1,3），（2,2），（2,3））都被转换为希尔伯特索引，系统可以直接检索所有匹配的条目。关键优势是这种映射使得数据检索既快又灵活。与传统索引不同，基于希尔伯特的分区可以自然地适应数据集的更新或增长，而无需重新计算整个索引。

5.2 数据索引：希尔伯特曲线与 Z 曲线

为了突出希尔伯特曲线的实际优势，我将其性能与一组合成范围查询中的 Z 曲线进行了比较。

对于实验，我生成了 100 个固定大小的随机范围查询。对于每个查询，我计算了希尔伯特曲线和 Z 曲线的索引，并计算了簇的数量，其中簇是一组连续的索引。例如，如果查询返回的索引为[1,2,3,5,6,8,9]，这将形成三个簇：[1,2,3]，[5,6]和[8,9]。

如果数据按索引顺序存储，簇对应于顺序读取，而簇之间的间隙意味着昂贵的跳转到新的存储地址。

图 10：使用希尔伯特曲线和 Z 曲线的 2D 设置的 100 个随机查询。如图所示，您什么也看不到！😉

（图片由作者提供）

为了量化性能，我使用了两个指标：

1. 簇计数：簇越少，意味着碎片化越少，存储跳跃越少。

2. 簇内扩散：每个簇的平均索引数

最坏的情况将是极端碎片化：每个点都形成一个自己的簇。图 11 比较了 Z 曲线和希尔伯特曲线在二维、三维和四维中的性能，查询大小为 7（2D 中的 7×7，3D 中的 7x7x7 等），每轴 6 位（即每轴 64 个值）

图 11：基于簇数和簇内扩散的希尔伯特曲线与 Z 曲线的比较（2，3 和 4 维）。结果清楚地表明，希尔伯特曲线比 Z 曲线更好地保留了局部性（图片由作者提供）

结果清楚地表明，希尔伯特曲线比 Z 曲线更好地保留了局部性。在所有测试的维度中，查询结果导致簇的数量更少，因此使用希尔伯特索引的簇内密度更高。在实践中，这转化为更有效的数据检索和降低 I/O 成本，尤其是对于多维范围查询。

6 超越空间填充曲线

本文的目标是展示 SFCs 的优雅性，并对其在数据索引中的应用给出一个概览。然而，该领域的最新研究已经超越了经典 SFCs。

所有空间填充曲线的主要局限性在于它们的固定机制。一旦定义，它们的结构为适应不同的数据集或工作负载模式提供了很少的适应空间。在实践中，这种刚性可能会限制性能。

为了克服这一点，像陈等研究人员（中国电子科技大学与华为）已经提出了 AdaCurve，这是一种基于机器学习的方法。AdaCurve 不是依赖于预定的映射，而是训练一个模型，直接从高维数据点生成一维索引，根据数据集和查询负载进行优化。[3]

这个想法非常有前景：虽然希尔伯特和其他 SFCs 提供了优雅但刚性的映射，但 AdaCurve 能够动态适应，生成一个针对当前数据和查询定制的索引系统。这种适应性可能为未来大规模数据平台上的显著更高效索引铺平道路。

参考文献

[1] H. Sagan, 空间填充曲线. Springer-Verlag, 1994.

[2] M. Bader, 空间填充曲线 – 科学计算中的应用导论. Springer-Verlag, 2013.

[3] X. CHEN, “使用学习自适应曲线优化高维数据的块跳过,” SIGMOD, 第 3 卷, 2025. [在线]. 可用:zheng-kai.com/paper/2025_sigmod_chen.pdf