N次元のポイントの膨大なセットがあります（数千万、Nは100に近い）。
空間的な局所性を維持しながら、これらの点を 1 つの次元にマッピングする必要があります。ヒルベルトの空間充填曲線を使用してそれを行いたいです。
各点について、曲線上の最も近い点を選びたいと思います。点のヒルベルト値 (曲線の始点から選択した点までの曲線の長さ) は、私が求める単一次元の値です。
計算は即時である必要はありませんが、まともな最新の家庭用 PC ハードウェアでは数時間以内で済むと思います。
実装に関する提案はありますか？私を助けるライブラリはありますか？（言語はあまり関係ありません。）
私はついに故障し、いくらかのお金を払いました。AIP (American Institute of Physics) には、C のソース コードを含むすばらしい短い記事があります。John Skilling による「ヒルベルト曲線のプログラミング」(AIP Conf. Proc. 707, 381 (2004) から) には、次のコードを含む付録があります。双方向のマッピング。1 を超える任意の数の次元に対して機能し、再帰的ではなく、大量のメモリを消費する状態遷移ルックアップ テーブルを使用せず、主にビット操作を使用します。したがって、適度に高速で、メモリ フットプリントが良好です。
記事を購入することを選択した場合、ソース コードにエラーが見つかりました。
次のコード行 (関数 TransposetoAxes にあります) にはエラーがあります。
for( i = n-1; i >= 0; i-- ) X[i] ^= X[i-1];
修正は、より大きいか等しい (>=) をより大きい (>) に変更することです。この修正を行わないと、変数「i」がゼロになると X 配列が負のインデックスを使用してアクセスされ、プログラムが失敗します。
この記事 (コードを含めて 7 ページの長さ) を読むことをお勧めします。この記事では、アルゴリズムがどのように機能するかが説明されていますが、これは明らかではありません。
私は彼のコードを自分で使用するために C# に翻訳しました。コードは次のとおりです。Skilling はその場で変換を実行し、渡されたベクトルを上書きします。入力ベクトルのクローンを作成し、新しいコピーを返すことにしました。また、メソッドを拡張メソッドとして実装しました。
Skilling のコードは、ヒルベルト インデックスを転置として表し、配列として格納します。ビットをインターリーブして単一の 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 = A B C D E F G H I J K L M N O is stored
/// as its Transpose ^
/// X[0] = A D G J M X[2]| 7
/// X[1] = B E H K N <-------> | /X[1]
/// X[2] = C F I L O 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を参照してください。
更新: 私はちょうど (2019 年秋) crates.io で "hilbert" と呼ばれる Rust クレートを公開しました。また、Skilling のアルゴリズムも使用します。https://crates.io/crates/hilbertを参照
ここで与えられた n->1 および 1->n からのマッピングのアルゴリズム 「ヒルベルト空間充填曲線を使用した 1 次元値と n 次元値の間のマッピングの計算」JK Lawder
「SFC モジュールと Kademlia オーバーレイ」を Google で検索すると、それをシステムの一部として使用していると主張するグループが見つかります。ソースを表示すると、おそらく関連する関数を抽出できます。
これがあなたが望むことをどのように行うかは私にはわかりません。この些細な3Dの場合を考えてみましょう。
001 ------ 101
|\ |\
| \ | \
| 011 ------ 111
| | | |
| | | |
000 -|---- 100 |
\ | \ |
\ | \ |
010 ------ 110
これは、次のパスで「Hilbertized」できます。
001 -----> 101
\ \
\ \
011 111
^ |
| |
000 | 100 |
\ | \ |
\ | \ V
010 110
1Dオーダーに：
000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100
これが厄介な部分です。以下のペアと1D距離のリストを検討してください。
000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1
すべての場合において、左側と右側の値は互いに同じ3D距離（最初の位置で+/- 1）であり、これは同様の「空間的局所性」を意味しているように見えます。ただし、次元の順序（上記の例では、y、z、zの順に）を選択して線形化すると、その局所性が失われます。
別の言い方をすれば、開始点を取得し、残りの点をその開始点からの距離で並べ替えると、結果が大幅に異なります。出発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
この効果は、次元の数とともに指数関数的に増大します（各次元の「サイズ」が同じであると想定）。
少し時間をかけて、Paul Chernoch のコードを Java に翻訳し、クリーンアップしました。コードにバグがある可能性があります。特に、元の論文にアクセスできないためです。しかし、それは私が書くことができた単体テストに合格します。以下です。
大規模なデータ セットの空間インデックス付けについて、 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;
}
もう 1 つの可能性は、データにkd ツリーを構築し、ツリーを順番に走査して順序を取得することです。kd ツリーを構築するには、優れた中央値検出アルゴリズムが必要なだけであり、その中には多数あります。
ヒルベルト曲線を 1 次元で使用する方法がわかりません。
距離を維持しながら (最小のエラーで) ポイントをより低い次元にマッピングすることに関心がある場合は、「多次元スケーリング」アルゴリズムを調べることができます。
シミュレーテッド アニーリングは 1 つのアプローチです。
編集：コメントありがとうございます。ヒルベルト曲線アプローチの意味がわかりました。しかし、これは難しい問題であり、N=100 と 1000 万のデータ ポイントを考えると、局所性を適切に維持し、妥当な時間内に実行できるアプローチはないと思います。ここでは kd-trees は機能しないと思います。
全体の順序を見つけることが重要でない場合は、局所性ベースのハッシングやその他のおおよその最近傍スキームを調べることができます。入力サイズを減らすためのポイントのバケットを使用した階層的多次元スケーリングは、適切な順序付けを提供する可能性がありますが、そのような高次元では疑わしいです。