Eigen - 计算两组向量之间的距离矩阵

Question

我需要创建一个Eigen数组，其中包含每个源和记录器位置之间的所有距离。

我有三个Eigen::Array，分别表示源位置和记录器位置的sx、sy和sz，以及rx、ry和rz。源和记录器位置数组的长度不一定相同。

使用for循环，我可以按照以下方式计算距离矩阵

Eigen::ArrayXf  sx, sy, sz;
Eigen::ArrayXf  rx, ry, rz;
Eigen::ArrayXXf dist;

for (int s=0; s<nrSources; s++ ) {

   for (int r=0; r<nrReceivers; r++) {

      dist(h,g) = sqrt(pow(rx(h)-sx(g),2.)+
                       pow(ry(h)-sy(g),2.)+
                       pow(sz(h)-sz(g),2.));
   }
}

我需要计算500 x 1000次实验中的dist数组。使用for循环显然可以工作，但这可能不是最有效的方法，因为它没有利用向量化。

将sx、sy、sz和hx、hy和hz重写为sxyz和hxyz数组应该是可能的。

是否有可能更有效地编写方程？

- Cnoobplusplus

Answer 1

您可能希望使用表达式来丢弃内部循环。这样做可以将计算合并在一起，有助于启用矢量化。假设r和s变量被声明为Eigen :: ArrayX3f sxyz（nrSources，3），rxyz（nrReceivers，3）;，您可以将外部循环编写为:

for (int s = 0; s < nrSources; s++)
{
    dist.col(s) =
        (rxyz.rowwise() - sxyz.row(s)).matrix().rowwise().norm();
}       

这里我们使用rxyz.rowwise() - sxyz.row(s)来从所有的r中减去第s 个 s。需要使用matrix()才能访问norm()。

进行基准测试并与您的实现进行比较。

- Avi Ginsburg

---

谢谢你提供这段代码。我想到也许可以去掉内部循环，但我觉得我可能会在中途卡住。 - Cnoobplusplus

由于这个操作非常容易并行化和矢量化，我想知道用原始的C代码编写它并让英特尔编译器进行优化是否能提供更快的性能表现。 - Royi

@Royi 可能是这样（这是 OP 中的内容）。并非每个人都可以访问 ICC。在简单情况下可能有效，但即使在这种情况下，一行代码也比普通代码更具表现力和紧凑性，并且具有知道表达式模板和矢量化已经（至少在某种程度上）优化并且将在更多处理器上运行良好（例如 AMD、ARM 等）而无需修改的附加优势。 - Avi Ginsburg

@AviGinsburg，我理解你的观点。问题是我不确定Eigen是否利用了AVX / AVX2。在上述情况下，这可能非常有益。 - Royi

@Royi 从3.3版本开始，Eigen支持AVX。MSVC通常不会生成良好的汇编代码与Eigen一起使用。尝试在Linux下使用gcc（即使是虚拟机也比Windows获得更好的性能）。确保启用--march=native、_DNDEBUG和-O3以进行更公平的比较。此外，我没有看到您在项目中启用了AVX（除非VS2017默认启用它，这我不知道）。 - Avi Ginsburg

@AviGinsburg，我之前没有设置_DNDEBUG，现在已经设置了，结果几乎相同。Eigen比普通实现（使用多线程）慢了约2.5倍。如果你仔细看，你会发现我正在使用GCC 7.3（MinGW64）。MSVC打败了它。我关于AVX的评论是Eigen从中受益不大。 - Royi

Answer 2

你可以查看我的项目计算距离矩阵。
我实现了多种方法来计算两个任意向量集之间的距离矩阵。

实际上，Eigen在执行此操作时非常快（比使用MSVC的原始C实现要慢得多）。

无论如何，另一种方法是：

void CalcDistanceMatrixEigen(float* mD, float* mA, float* mB, int vecDim, int numRowsA, int numRowsB)
{

    EigenMatExt meA(mA, vecDim, numRowsA);
    EigenMatExt meB(mB, vecDim, numRowsB);

    EigenMatExt meD(mD, numRowsA, numRowsB);

    // meD = meA.colwise().squaredNorm().transpose() - (2 * meA.transpose() * meB) + meB.colwise().squaredNorm();
    // Not even close to be as fast as MATLAB (Intel MKL???)
    meD = ((-2 * meA.transpose() * meB).colwise() + meA.colwise().squaredNorm().transpose()).rowwise() + meB.colwise().squaredNorm();


}

只需要查看文件。
你会发现我在@AviGinsburg和上述方法中都比较了Eigen，似乎我的方法更快（尽管再次说明，对于该操作使用Vanilla C比Eigen更快）。

- Royi