c++performancelinear-algebradistanceeigen

Eigen - Calculate the Distance Matrix between 2 Sets of Vectors


I need to create an Eigen array with all distances between each source and recorder position.

I have three Eigen::Array sx, sy and sz denoting the source positions and three Eigen::Array rx, ry and rz denoting the recorder positions. The length of source and recording position arrays are not necessary the same.

Using for-loops I can calculate the distance matrix as follow

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.));
   }
}

I need to compute dist array like 500 x 1000 times for the experiment. Using the for-loop obvious works but it is probably not the most efficient method since it does not make use of the vectorization.

Rewriting the sx, sx, sz and hx, hy and hz into sxyz and hxyz arrays should be possible.

Is it possible to write the equation more efficiently!?


Solution

  • You probably want to discard the inner loop in favor of an expression. By doing so, we allow calculations to be lumped together and hopefully enable vectorization. Assuming the r and s variables were declared as Eigen::ArrayX3f sxyz(nrSources,3), rxyz(nrReceivers,3);, you could then write the outer loop as:

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

    Here we use rxyz.rowwise() - sxyz.row(s) to subtract the s'th s from all rs. The matrix() is needed to give access to norm().

    Benchmark and compare to your implementation.