在N维空间中计算两点之间的欧几里得距离的最快方法

Question

我需要计算N维空间中两点之间的欧几里得距离，速度至关重要。我有两个C风格的浮点数组，表示N维空间中的这两个点。

它们之间的距离公式如下（^ 表示乘方，不是异或）： sqrt(sum((p1-q1)^2 + (p2-q1)^2 + .... (pn-qn)^2))

我的当前代码如下：

sum = 0;
for(int i=0;i<N;++i){
    sum += pow(p[i]-q[i],2);
sqrt(sum)

这段代码运行速度较慢，我想知道是否有任何库可以加速它？我想象中已经有人编写了一个快速的库，在c语言中对数组执行数学运算，让我能够快速地对数组进行逐元素操作。

编辑： 回答nevsan的问题，我正在使用小N进行多个计算，大约是10或20。

- user1357607

Answer 1

一定要摆脱pow()。这在很大程度上取决于你如何使用它进行优化。你是针对非常大的N只做一次并且花费太长时间吗？还是更有可能是在一个紧密循环中多次执行？

如果你正在使用非常大的N（>1000左右），那么有高度优化的数值库可以做到这一点。例如，BLAS具有一个*nrm2函数，可以计算欧几里得范数（dnrm2、snrm2、cnrm2、znrm2，根据数据类型而定[单精度、双精度、复杂单精度、复杂双精度]）。GotoBLAS可能是某些处理器架构中最快的。MKL特色是英特尔手工调整的BLAS实现，但它不是免费的。最后，ATLAS是一个自我调整的库，实现了BLAS。

如果您有一个紧密的循环，N 很小或不是很大，那么您可能需要手动调整一些参数以使其更快。您可以使用编译器标志 -O3 或 -ftree-vectorize 来启用自动向量化。您也可以手动进行矢量化，但学习如何这样做可能会很痛苦。

您可以进行循环展开（即将 N 分成一些块，比如 4，并在 for 循环体内明确编写 4 个连续值的计算公式）。这样可以欺骗编译器使用更多寄存器进行立即计算——而寄存器是您要处理的最快形式的内存。此外，您可能还可以利用预取（使用一个内存访问调用读取一段数据）。

在这种情况下，另一件要做的事是尝试覆盖其中一个输入。也就是说，您可以将输出写入

或中，以此来帮助您计算的< p >位置仍然在缓存中，当您准备好写入时，缓存通常不会将数据写入内存，除非他们绝对必须这样做---原因之一是需要缓存行并且我们需要将上一个踢出。通过写入其中一个输入来使用更少的缓存行。

还有500,000其他要尝试的事情，但我想在这里停止。祝你好运！

- nevsan

---

除了向量化之外，如果适用于您的用例，您还可以尝试在多个核心上执行计算。 - user1118321

不幸的是，BLAS 提供的欧几里得范数函数实际上并不是我想要的。欧几里得范数或欧几里得长度相当于计算N维空间中的一个点到原点的距离，而不是两个不同点之间的距离。 - user1357607

@user1357607 是的，你必须先使用 *axpy 函数来计算两个向量的差。代码应该类似于 *axpy(N, -1, p, 1, q, 1); answer = *nrm2(N, p, 1);。 - nevsan

谢谢。冒昧问一下，是否还有一个函数可以将两个向量逐元素相乘？ - user1357607

MKL现在是免费的，而且速度相当快。 - DragonLord

Answer 2

我从不使用pow() - 我猜测，没有进行性能分析的话，这会严重拖慢你的速度。

你需要创建一个临时变量，然后对其进行平方。

double diff = p[i] - q[i];
sum += diff*diff;

sqrt函数速度有点慢，但这里唯一的选择是一些近似值。如果N大于约10，那么sqrt很可能不会成为瓶颈。

还有像boost等库，可以加快速度，但首先尝试摆脱pow()。请记住，diff*diff只需要一个浮点指令，而pow()是一个专门设计用于非整数幂的整个程序。

- Tom Andersen