对我来说最好的最近邻算法是什么?
What is the best nearest neighbor algorithm for my case?
我有一个预定义的 gps 位置列表,它基本上构成了一个预定义的汽车轨迹。列表中大约有 15000 个点。整个列表是事先已知的,之后不需要插入任何点。然后我得到大约 1 百万 个额外的采样 gps 位置,我需要在预定义列表中找到最近的邻居。我需要在一次迭代中处理所有 100 万个项目,并且我需要尽快完成。这种情况下最好的最近邻算法是什么?
我可以根据需要尽可能多地预处理预定义列表,但是处理 100 万个项目应该尽可能快。
我已经测试了 KDTree c# 实现,但性能似乎很差,也许存在更适合我的 2D 数据的算法。 (在我的例子中 gps 高度被忽略)
感谢您的任何建议!
K-D 树确实很适合这个问题。你应该先用 known-good 实现再试一次,如果性能不够好,你可以轻松地并行化查询——因为每个查询都完全独立于其他查询,你可以通过处理 N 个查询来实现 N 的加速并行,如果你有足够的硬件。
我推荐OpenCV的implementation, as mentioned in this answer
Performance-wise,您插入的点的顺序 可以 对查询时间有影响,因为实现可能会选择是否重新平衡不平衡的树(和,例如,OpenCV 不会这样做)。一个简单的安全措施是以随机顺序插入点:首先打乱列表,然后以打乱的顺序插入所有点。虽然不是最优的,但这确保了以压倒性的概率,结果顺序不会是病态的。
CGAL 有一个 2d point library 用于基于 Delaunay 三角剖分数据结构的最近邻和范围搜索。
这是他们的库针对您的用例的基准测试:
// file: cgal_benchmark_2dnn.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_2.h>
#include <chrono>
#include <list>
#include <random>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Point_set_2<K>::Vertex_handle Vertex_handle;
typedef K::Point_2 Point_2;
/**
* @brief Time a lambda function.
*
* @param lambda - the function to execute and time
*
* @return the number of microseconds elapsed while executing lambda
*/
template <typename Lambda>
std::chrono::microseconds time_lambda(Lambda lambda) {
auto start_time = std::chrono::high_resolution_clock::now();
lambda();
auto end_time = std::chrono::high_resolution_clock::now();
return std::chrono::duration_cast<std::chrono::microseconds>(end_time -
start_time);
}
int main() {
const int num_index_points = 15000;
const int num_trials = 1000000;
std::random_device
rd; // Will be used to obtain a seed for the random number engine
std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
std::uniform_real_distribution<> dis(-1, 1.);
std::list<Point_2> index_point_list;
{
auto elapsed_microseconds = time_lambda([&] {
for (int i = 0; i < num_index_points; ++i) {
index_point_list.emplace_back(dis(gen), dis(gen));
}
});
std::cout << " Generating " << num_index_points << " random points took "
<< elapsed_microseconds.count() << " microseconds.\n";
}
CGAL::Point_set_2<K> point_set;
{
auto elapsed_microseconds = time_lambda([&] {
point_set.insert(index_point_list.begin(), index_point_list.end());
});
std::cout << " Building point set took " << elapsed_microseconds.count()
<< " microseconds.\n";
}
{
auto elapsed_microseconds = time_lambda([&] {
for (int j = 0; j < num_trials; ++j) {
Point_2 query_point(dis(gen), dis(gen));
Vertex_handle v = point_set.nearest_neighbor(query_point);
}
});
auto rate = elapsed_microseconds.count() / static_cast<double>(num_trials);
std::cout << " Querying " << num_trials << " random points took "
<< elapsed_microseconds.count()
<< " microseconds.\n >> Microseconds / query :" << rate << "\n";
}
}
在我的系统 (Ubuntu 18.04) 上可以用
编译
g++ cgal_benchmark_2dnn.cpp -lCGAL -lgmp -O3
当 运行 产生性能时:
Generating 15000 random points took 1131 microseconds.
Building point set took 11469 microseconds.
Querying 1000000 random points took 2971201 microseconds.
>> Microseconds / query :2.9712
速度相当快。请注意,使用 N 个处理器,您可以将其速度提高大约 N 倍。
最快的实施
如果以下两项或多项为真:
- 您有一个用于 150000 个索引点的小边界框
- 您只关心小数点后几位的精度(请注意,对于小数点后 6 位以上的纬度和经度坐标,会产生 centimeter/millimeter 比例精度)
- 您的系统内存充足
然后缓存所有内容! 您可以 pre-compute 在索引点的边界框上创建所需精度的网格。将每个网格单元格映射到一个唯一的地址,该地址可以在知道查询点的二维坐标的情况下进行索引。
然后简单地使用任何最近邻算法(例如我提供的算法)将每个网格单元映射到最近的索引点。请注意,此步骤只需执行一次即可初始化网格中的网格单元。
对于 运行 一个查询,这将需要一个 2D 坐标到网格单元格坐标计算,然后是一个内存访问,这意味着您不能真正希望有一种更快的方法(可能是 2-3 CPU 个查询周期。)
我怀疑(有一些洞察力)这就是像 Google 或 Facebook 这样的大公司会如何解决这个问题(因为 #3 对他们来说甚至对整个世界来说都不是问题。)更小的 non-profit 组织使用这样的方案(例如 NASA。)尽管如此,NASA 使用的方案要复杂得多,具有 resolution/precision.
的多个尺度
澄清
从下面的评论来看,很明显最后一节没有被很好地理解,所以我会添加更多细节。
假设您的点集由两个向量 x
和 y
给出,它们包含数据的 x 和 y 坐标(或纬度和经度或您正在使用的任何坐标。)
然后数据的边界框定义为维度 width = max(x)-min(x)
& height=max(y)-min(y)
。
现在使用一组测试点 (x_t,y_t)
的映射,使用 NxM 个点创建一个精细网格来表示整个边界框
u(x_t) = round((x_t - min(x)) / double(width) * N)
v(y_t) = round((y_t - min(y)) / double(height) * M)
然后只需使用 indices = grid[u(x_t),v(y_t)]
,其中 indices
是最接近 [x_t,y_t]
的索引点的索引,grid
是预先计算的查找 table将网格中的每个项目映射到最近的索引点 [x,y]
.
例如,假设您的索引点是 [0,0]
和 [2,2]
(按此顺序)。您可以将网格创建为
grid[0,0] = 0
grid[0,1] = 0
grid[0,2] = 0 // this is a tie
grid[1,0] = 0
grid[1,1] = 0 // this is a tie
grid[1,2] = 1
grid[2,0] = 1 // this is a tie
grid[2,1] = 1
grid[2,2] = 1
上面的右侧是索引 0
(映射到点 [0,0]
)或 1
(映射到点 [2,2]
)。注意:由于这种方法的离散性,您会遇到与一个点的距离恰好等于到另一个索引点的距离的关系,您将不得不想出一些方法来确定如何打破这些关系。请注意,grid
中的条目数决定了您要达到的精确度。显然在我上面给出的例子中精度很差。
我有一个预定义的 gps 位置列表,它基本上构成了一个预定义的汽车轨迹。列表中大约有 15000 个点。整个列表是事先已知的,之后不需要插入任何点。然后我得到大约 1 百万 个额外的采样 gps 位置,我需要在预定义列表中找到最近的邻居。我需要在一次迭代中处理所有 100 万个项目,并且我需要尽快完成。这种情况下最好的最近邻算法是什么?
我可以根据需要尽可能多地预处理预定义列表,但是处理 100 万个项目应该尽可能快。
我已经测试了 KDTree c# 实现,但性能似乎很差,也许存在更适合我的 2D 数据的算法。 (在我的例子中 gps 高度被忽略)
感谢您的任何建议!
K-D 树确实很适合这个问题。你应该先用 known-good 实现再试一次,如果性能不够好,你可以轻松地并行化查询——因为每个查询都完全独立于其他查询,你可以通过处理 N 个查询来实现 N 的加速并行,如果你有足够的硬件。
我推荐OpenCV的implementation, as mentioned in this answer
Performance-wise,您插入的点的顺序 可以 对查询时间有影响,因为实现可能会选择是否重新平衡不平衡的树(和,例如,OpenCV 不会这样做)。一个简单的安全措施是以随机顺序插入点:首先打乱列表,然后以打乱的顺序插入所有点。虽然不是最优的,但这确保了以压倒性的概率,结果顺序不会是病态的。
CGAL 有一个 2d point library 用于基于 Delaunay 三角剖分数据结构的最近邻和范围搜索。
这是他们的库针对您的用例的基准测试:
// file: cgal_benchmark_2dnn.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_2.h>
#include <chrono>
#include <list>
#include <random>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Point_set_2<K>::Vertex_handle Vertex_handle;
typedef K::Point_2 Point_2;
/**
* @brief Time a lambda function.
*
* @param lambda - the function to execute and time
*
* @return the number of microseconds elapsed while executing lambda
*/
template <typename Lambda>
std::chrono::microseconds time_lambda(Lambda lambda) {
auto start_time = std::chrono::high_resolution_clock::now();
lambda();
auto end_time = std::chrono::high_resolution_clock::now();
return std::chrono::duration_cast<std::chrono::microseconds>(end_time -
start_time);
}
int main() {
const int num_index_points = 15000;
const int num_trials = 1000000;
std::random_device
rd; // Will be used to obtain a seed for the random number engine
std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
std::uniform_real_distribution<> dis(-1, 1.);
std::list<Point_2> index_point_list;
{
auto elapsed_microseconds = time_lambda([&] {
for (int i = 0; i < num_index_points; ++i) {
index_point_list.emplace_back(dis(gen), dis(gen));
}
});
std::cout << " Generating " << num_index_points << " random points took "
<< elapsed_microseconds.count() << " microseconds.\n";
}
CGAL::Point_set_2<K> point_set;
{
auto elapsed_microseconds = time_lambda([&] {
point_set.insert(index_point_list.begin(), index_point_list.end());
});
std::cout << " Building point set took " << elapsed_microseconds.count()
<< " microseconds.\n";
}
{
auto elapsed_microseconds = time_lambda([&] {
for (int j = 0; j < num_trials; ++j) {
Point_2 query_point(dis(gen), dis(gen));
Vertex_handle v = point_set.nearest_neighbor(query_point);
}
});
auto rate = elapsed_microseconds.count() / static_cast<double>(num_trials);
std::cout << " Querying " << num_trials << " random points took "
<< elapsed_microseconds.count()
<< " microseconds.\n >> Microseconds / query :" << rate << "\n";
}
}
在我的系统 (Ubuntu 18.04) 上可以用
编译g++ cgal_benchmark_2dnn.cpp -lCGAL -lgmp -O3
当 运行 产生性能时:
Generating 15000 random points took 1131 microseconds.
Building point set took 11469 microseconds.
Querying 1000000 random points took 2971201 microseconds.
>> Microseconds / query :2.9712
速度相当快。请注意,使用 N 个处理器,您可以将其速度提高大约 N 倍。
最快的实施
如果以下两项或多项为真:
- 您有一个用于 150000 个索引点的小边界框
- 您只关心小数点后几位的精度(请注意,对于小数点后 6 位以上的纬度和经度坐标,会产生 centimeter/millimeter 比例精度)
- 您的系统内存充足
然后缓存所有内容! 您可以 pre-compute 在索引点的边界框上创建所需精度的网格。将每个网格单元格映射到一个唯一的地址,该地址可以在知道查询点的二维坐标的情况下进行索引。
然后简单地使用任何最近邻算法(例如我提供的算法)将每个网格单元映射到最近的索引点。请注意,此步骤只需执行一次即可初始化网格中的网格单元。
对于 运行 一个查询,这将需要一个 2D 坐标到网格单元格坐标计算,然后是一个内存访问,这意味着您不能真正希望有一种更快的方法(可能是 2-3 CPU 个查询周期。)
我怀疑(有一些洞察力)这就是像 Google 或 Facebook 这样的大公司会如何解决这个问题(因为 #3 对他们来说甚至对整个世界来说都不是问题。)更小的 non-profit 组织使用这样的方案(例如 NASA。)尽管如此,NASA 使用的方案要复杂得多,具有 resolution/precision.
的多个尺度澄清
从下面的评论来看,很明显最后一节没有被很好地理解,所以我会添加更多细节。
假设您的点集由两个向量 x
和 y
给出,它们包含数据的 x 和 y 坐标(或纬度和经度或您正在使用的任何坐标。)
然后数据的边界框定义为维度 width = max(x)-min(x)
& height=max(y)-min(y)
。
现在使用一组测试点 (x_t,y_t)
u(x_t) = round((x_t - min(x)) / double(width) * N)
v(y_t) = round((y_t - min(y)) / double(height) * M)
然后只需使用 indices = grid[u(x_t),v(y_t)]
,其中 indices
是最接近 [x_t,y_t]
的索引点的索引,grid
是预先计算的查找 table将网格中的每个项目映射到最近的索引点 [x,y]
.
例如,假设您的索引点是 [0,0]
和 [2,2]
(按此顺序)。您可以将网格创建为
grid[0,0] = 0
grid[0,1] = 0
grid[0,2] = 0 // this is a tie
grid[1,0] = 0
grid[1,1] = 0 // this is a tie
grid[1,2] = 1
grid[2,0] = 1 // this is a tie
grid[2,1] = 1
grid[2,2] = 1
上面的右侧是索引 0
(映射到点 [0,0]
)或 1
(映射到点 [2,2]
)。注意:由于这种方法的离散性,您会遇到与一个点的距离恰好等于到另一个索引点的距离的关系,您将不得不想出一些方法来确定如何打破这些关系。请注意,grid
中的条目数决定了您要达到的精确度。显然在我上面给出的例子中精度很差。