将N维值映射到希尔伯特曲线上的点

我有一个庞大的N维点（数千万，N接近100）。

我需要将这些点映射到一个维度，同时保留空间局部性。我想用希尔伯特空间填充曲线来做到这一点。

对于每个点，我想select曲线上最近的点。该点的Hilbert值（从曲线起点到拾取点的曲线长度）是I seek的单维值。

计算不一定是即时的，但我期望在不错的现代家庭PC硬件上不超过几个小时。

任何关于实施的build议？有没有任何图书馆可以帮助我？（语言不重要。）

我终于打破了一些钱。 AIP（美国物理研究所）有一篇很好的短文，其源代码为C. John Skilling（来自AIP Conf。Proc。707,381（2004））的“编程希尔伯特曲线”有一个附录，其代码为在两个方向映射。它适用于任何维度> 1，不recursion，不使用状态转换查找表吞噬大量的内存，并且大多使用位操作。因此它速度相当快，并具有良好的内存占用。

如果您select购买文章，我发现源代码中有错误。

以下代码行（在函数TransposetoAxes中find）有错误：

对于（i = n-1; i> = 0; i–）X [i] ^ = X [i-1]

更正是将大于或等于（> =）更改为大于（>）。如果没有这个更正，当variables“i”变为零时，使用负索引访问X数组，导致程序失败。

我build议阅读这篇文章（这是长达七页，包括代码），因为它解释了algorithm是如何工作的，这是不明显的。

我将他的代码翻译成C＃供我自己使用。代码如下。 Skilling会执行转换，覆盖您传入的向量。我select复制input向量并返回新副本。另外，我实现了作为扩展方法的方法。

Skilling的代码将Hilbert索引表示为一个转置，作为一个数组存储。我发现交叉位和形成单个BigInteger（在字典中更有用，在循环中更容易迭代等）更方便，但是我使用幻数，位操作等来优化该操作及其反转。代码很长，所以我省略了它。

namespace HilbertExtensions { /// <summary> /// Convert between Hilbert index and N-dimensional points. /// /// The Hilbert index is expressed as an array of transposed bits. /// /// Example: 5 bits for each of n=3 coordinates. /// 15-bit Hilbert integer = ABCDEFGHIJKLMNO is stored /// as its Transpose ^ /// X[0] = ADGJMX[2]| 7 /// X[1] = BEHKN <-------> | /X[1] /// X[2] = CFILO axes |/ /// high low 0------> X[0] /// /// NOTE: This algorithm is derived from work done by John Skilling and published in "Programming the Hilbert curve". /// (c) 2004 American Institute of Physics. /// /// </summary> public static class HilbertCurveTransform { /// <summary> /// Convert the Hilbert index into an N-dimensional point expressed as a vector of uints. /// /// Note: In Skilling's paper, this function is named TransposetoAxes. /// </summary> /// <param name="transposedIndex">The Hilbert index stored in transposed form.</param> /// <param name="bits">Number of bits per coordinate.</param> /// <returns>Coordinate vector.</returns> public static uint[] HilbertAxes(this uint[] transposedIndex, int bits) { var X = (uint[])transposedIndex.Clone(); int n = X.Length; // n: Number of dimensions uint N = 2U << (bits - 1), P, Q, t; int i; // Gray decode by H ^ (H/2) t = X[n - 1] >> 1; // Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index. for (i = n - 1; i > 0; i--) X[i] ^= X[i - 1]; X[0] ^= t; // Undo excess work for (Q = 2; Q != N; Q <<= 1) { P = Q - 1; for (i = n - 1; i >= 0; i--) if ((X[i] & Q) != 0U) X[0] ^= P; // invert else { t = (X[0] ^ X[i]) & P; X[0] ^= t; X[i] ^= t; } } // exchange return X; } /// <summary> /// Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve. /// That distance will be transposed; broken into pieces and distributed into an array. /// /// The number of dimensions is the length of the hilbertAxes array. /// /// Note: In Skilling's paper, this function is called AxestoTranspose. /// </summary> /// <param name="hilbertAxes">Point in N-space.</param> /// <param name="bits">Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.</param> /// <returns>The Hilbert distance (or index) as a transposed Hilbert index.</returns> public static uint[] HilbertIndexTransposed(this uint[] hilbertAxes, int bits) { var X = (uint[])hilbertAxes.Clone(); var n = hilbertAxes.Length; // n: Number of dimensions uint M = 1U << (bits - 1), P, Q, t; int i; // Inverse undo for (Q = M; Q > 1; Q >>= 1) { P = Q - 1; for (i = 0; i < n; i++) if ((X[i] & Q) != 0) X[0] ^= P; // invert else { t = (X[0] ^ X[i]) & P; X[0] ^= t; X[i] ^= t; } } // exchange // Gray encode for (i = 1; i < n; i++) X[i] ^= X[i - 1]; t = 0; for (Q = M; Q > 1; Q >>= 1) if ((X[n - 1] & Q)!=0) t ^= Q - 1; for (i = 0; i < n; i++) X[i] ^= t; return X; } } }

我已经把C＃中的工作代码发布到github上。

请参阅https://github.com/paulchernoch/HilbertTransformation

从n> 1和1 – > n映射的algorithm在这里给出“使用希尔伯特空间填充曲线计算一维和n维值之间的映射”JK Lawder

如果您是Google的“SFC模块和Kademlia覆盖”，您可以find一个声称将其用作系统一部分的组。如果你查看源代码，你可能会提取相关的function。

我不清楚这将如何做你想要的。考虑这个trival 3D案例：

 001 ------ 101 |\ |\ | \ | \ | 011 ------ 111 | | | | | | | | 000 -|---- 100 | \ | \ | \ | \ | 010 ------ 110

可以通过以下path“Hilbertized”：

 001 -----> 101 \ \ \ \ 011 111 ^ | | | 000 | 100 | \ | \ | \ | \ V 010 110

进入1D命令：

 000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100

这是讨厌的一点。考虑以下对的列表和一维距离：

 000 : 100 -> 7 010 : 110 -> 5 011 : 111 -> 3 001 : 101 -> 1

在所有情况下，左手和右手的距离都是相同的3D距离（在第一个位置+/- 1），这似乎意味着相似的“空间局部性”。但是，通过维度sorting（y，然后z，然后在上面的例子中，z）的任何select线性化打破该地点。

另一种说法是，以出发点为起点，按距出发点的距离来sorting剩余点，将会产生明显不同的结果。以000为起点，例如：

 1D ordering : distance 3D ordering : distance ---------------------- ---------------------- 010 : 1 001,010,100 : 1 011,101,110 : sqrt(2) 111 : sqrt(3) 011 : 2 001 : 3 101 : 4 111 : 5 110 : 6 100 : 7

这种效应随着维数的增加呈指数增长（假设每个维度具有相同的“大小”）。

另一种可能性是在你的数据上构build一个kd-tree ，然后到树的顺序遍历来获得sorting。构buildkd树只需要你有一个很好的中值searchalgorithm，其中有很多。

我花了一点时间把Paul Chernoch的代码翻译成Java并清理了它。我的代码中可能存在一个错误，尤其是因为我无法访问它原来的文件。但是，它通过了我能写的unit testing。在下面。

请注意，我已经评估了大数据集上的空间索引的Z-阶和希尔伯特曲线。我不得不说，Z-Order提供了更好的质量。但随时为自己尝试。

  /** * Convert the Hilbert index into an N-dimensional point expressed as a vector of uints. * * Note: In Skilling's paper, this function is named TransposetoAxes. * @param transposedIndex The Hilbert index stored in transposed form. * @param bits Number of bits per coordinate. * @return Point in N-space. */ static long[] HilbertAxes(final long[] transposedIndex, final int bits) { final long[] result = transposedIndex.clone(); final int dims = result.length; grayDecode(result, dims); undoExcessWork(result, dims, bits); return result; } static void grayDecode(final long[] result, final int dims) { final long swap = result[dims - 1] >>> 1; // Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index. for (int i = dims - 1; i > 0; i--) result[i] ^= result[i - 1]; result[0] ^= swap; } static void undoExcessWork(final long[] result, final int dims, final int bits) { for (long bit = 2, n = 1; n != bits; bit <<= 1, ++n) { final long mask = bit - 1; for (int i = dims - 1; i >= 0; i--) if ((result[i] & bit) != 0) result[0] ^= mask; // invert else swapBits(result, mask, i); } } /** * Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve. * That distance will be transposed; broken into pieces and distributed into an array. * * The number of dimensions is the length of the hilbertAxes array. * * Note: In Skilling's paper, this function is called AxestoTranspose. * @param hilbertAxes Point in N-space. * @param bits Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve. * @return The Hilbert distance (or index) as a transposed Hilbert index. */ static long[] HilbertIndexTransposed(final long[] hilbertAxes, final int bits) { final long[] result = hilbertAxes.clone(); final int dims = hilbertAxes.length; final long maxBit = 1L << (bits - 1); inverseUndo(result, dims, maxBit); grayEncode(result, dims, maxBit); return result; } static void inverseUndo(final long[] result, final int dims, final long maxBit) { for (long bit = maxBit; bit != 0; bit >>>= 1) { final long mask = bit - 1; for (int i = 0; i < dims; i++) if ((result[i] & bit) != 0) result[0] ^= mask; // invert else swapBits(result, mask, i); } // exchange } static void grayEncode(final long[] result, final int dims, final long maxBit) { for (int i = 1; i < dims; i++) result[i] ^= result[i - 1]; long mask = 0; for (long bit = maxBit; bit != 0; bit >>>= 1) if ((result[dims - 1] & bit) != 0) mask ^= bit - 1; for (int i = 0; i < dims; i++) result[i] ^= mask; } static void swapBits(final long[] array, final long mask, final int index) { final long swap = (array[0] ^ array[index]) & mask; array[0] ^= swap; array[index] ^= swap; }

~~我不知道如何在一个维度上使用希尔伯特曲线。~~

如果您有兴趣将点映射到较低的维度，同时保留距离（误差最小），那么您可以查看“多维比例”algorithm。

模拟退火是一种方法。

编辑：感谢您的评论。我明白了你现在希尔伯特曲线的含义。然而，这是一个很难的问题，在N = 100和1000万个数据点的情况下，我认为任何方法都不能很好地保持局部性，并在合理的时间内运行。我不认为kd-trees会在这里工作。

如果find一个总的顺序对你来说并不重要，那么你可以看看基于局部的散列和其他近似最近邻居scheme。分级多维缩放与点的桶减lessinput大小可能会给你一个很好的sorting，但在这么高的维度再次是可疑的。