OpenMP SIMD Euclidean Distance

Question

I'm using cv::norm() to compute the L2 distance between row vectors in thousands of dimensions (at least 4096 dimensions). I'm using an AVX2 machine and icpc version 2017.

Since I'm computing this function millions of times, I need the possible most efficient version for this function and I highly doubt that cv::norm() exploits vectorization properly.

So I implemnted this function:

#include <opencv2/opencv.hpp>

float euclDist (const cv::Mat1f &first, const cv::Mat1f &second){
    float dist = 0;
    const float *firstp = first.ptr<float>(0);
    const float *secondp = second.ptr<float>(0);
    #pragma omp simd reduction(+:dist)
    for(size_t i=0; i<first.cols; i++)
        dist += std::pow(firstp[i]-secondp[i],2);
    return std::sqrt(dist);
}

int main (int argc, const char *argv[]){
    int codeSize = 4096;
    cv::Mat1f mat(1, codeSize); 
    float low = 0;
    float high = 2.0;
     srand (0);
     cv::theRNG().state = cv::getTickCount() ;
    randu(mat, cv::Scalar(low), cv::Scalar(high));
    euclDist(mat, mat);
    return 0;
}

The example above is similar what the actually program does: it generates a row vector and it computes the distance with itself (i.e. 0). In the actual program generates a matrix and it randomly compute the distance between some rows. I don't think that I can get only that part of the code easily.

I already noticed an improvement w.r.t. cv::norm(), but scales badly when computing multiple distances in parallel. According to Intel VTune I fully use all the processors, but I get poor memory performance, as it shown from this general exploration analysis:

This is the result from: sudo lshw -short -C memory

H/W path       Device      Class          Description
=====================================================
/0/0                       memory         128KiB BIOS
/0/4/b                     memory         32KiB L1 cache
/0/4/c                     memory         256KiB L2 cache
/0/4/d                     memory         6MiB L3 cache
/0/a                       memory         32KiB L1 cache
/0/2a                      memory         12GiB System Memory
/0/2a/0                    memory         4GiB SODIMM DDR3 Synchronous 1600 MHz (0.6 ns)
/0/2a/1                    memory         DIMM [empty]
/0/2a/2                    memory         8GiB SODIMM DDR3 Synchronous 1600 MHz (0.6 ns)
/0/2a/3                    memory         DIMM [empty]

What I could do?

Answer 1

I had to adapt the test program to get meaningful results using /usr/bin/time:

#include <iostream>
int main()
{
    static const int codeSize = 4096;
    cv::Mat1f mat1(1, codeSize);
    cv::Mat1f mat2(1, codeSize);
    float low = 0;
    float high = 2.0;
    srand(0);
    cv::theRNG().state = cv::getTickCount() ;

    randu(mat1, cv::Scalar(low), cv::Scalar(high));
    randu(mat2, cv::Scalar(low), cv::Scalar(high));

    double d = 0;
    for (auto i = 0;  i < 1000000;  ++i) {
        d += euclDist(mat1, mat2);
    }

    std::cout << d << std::endl;
}

(I used two different inputs to avoid a predictable result of 0, made sure I used the result to prevent optimising out, and executed a million iterations to get the run time over a second).

Here's my compilation flags (used with my standard Makefile):

167501: PKGS += opencv
167501: CXXFLAGS += -O3 -march=native
167501: CXXFLAGS += -fopenmp

---

Having done that, I tried a couple of variations, all compiled with GCC version 7.1.0, and executed them on an Intel i7-6700K with 128KiB/1MiB/8MiB L1d/2/3 caches.

The original clocks in at 14.6 seconds user time without the OpenMP directive and 4.5 seconds with it.

Changing the std::pow() call to a plain multiply improved those times to 4.2 and 0.6 seconds respectively.

I got another factor of three improvement in wall-clock time, at the expense of much more processor time, by distributing across 8 cores with

#pragma omp parallel for simd reduction(+:dist)

For this code, the loop is a bit small for parallelization, it seems.

---

float euclDist(const cv::Mat1f& first, const cv::Mat1f& second)
{
    float dist = 0;
    const float *const firstp = first.ptr<float>(0);
    const float *const secondp = second.ptr<float>(0);
#pragma omp simd reduction(+:dist)
    for(auto i = 0;  i < first.cols;  ++i) {
        auto const d = firstp[i] - secondp[i];
        dist += d * d;
    }
    return std::sqrt(dist);
}

Answer 2

Avoid sqrt if possible

One thing that can give a large speed improvement is if you can avoid computing the square root, and just work with squared distances instead. This works fine for things like finding which point is closest to a specified point, or all that within a given distance of a specified point.

This won't work if (for example) you're trying to find the shortest route between two points, so you're adding up distances between points.

Summary: you can't always do this, and makes the biggest difference with lower dimensions (because the square roots are a higher percentage of the computations) but will probably still give some improvement if you can use it.

Consider a different algorithm

If you really need to compute the distance instead of the squared distance, there are other algorithms for computing it that are often faster than sqrt(X * X + Y * Y).

One example was published years ago in More Programming Pearls: Confessions of a Coder. C++ code for that algorithm looks like this:

#include <math.h>

/* Iterations   Accuracy
 *  2          6.5 digits
 *  3           20 digits
 *  4           62 digits
 * assuming a numeric type able to maintain that degree of accuracy in
 * the individual operations.
 */
#define ITER 2

double dist(double P, double Q) {
/* A reasonably robust method of calculating `sqrt(P*P + Q*Q)'
 *
 * Transliterated from _More Programming Pearls, Confessions of a Coder_
 * by Jon Bentley, pg. 156.
 */

    double R;
    int i;

    P = fabs(P);
    Q = fabs(Q);

    if (P<Q) {
        R = P;
        P = Q;
        Q = R;
    }

/* The book has this as:
 *  if P = 0.0 return Q; # in AWK
 * However, this makes no sense to me - we've just insured that P>=Q, so
 * P==0 only if Q==0;  OTOH, if Q==0, then distance == P...
 */
    if ( Q == 0.0 )
        return P;

    for (i=0;i<ITER;i++) {
        R = Q / P;
        R = R * R;
        R = R / (4.0 + R);
        P = P + 2.0 * R * P;
        Q = Q * R;
    }
    return P;
}

Do note that initial comment about accuracy though--basically, you need two iterations for single precision, and 3 for full double precision (obviously, two iterations is generally faster than 3).

People have sometimes misinterpreted the "reasonably robust" comment. It's trying to point out that this is substantially more robust than the usual formula, not that it's any less so.

As it stands, this is written for 2D distances though--I haven't extended it to higher dimensions, but believe it's supposed to be fairly easy to do so.

Use of AVX

The easiest way to use AVX well (for this task) is probably to just unroll the loop, and compute a number of distances in parallel. Basically:

for all the X/Y values:
    load 8 X values into register 0
    load 8 Y values into register 1
    square values in register 0
    square values in register 1
    add values in register 1 to register 0
    add results to accumulated distances
next 8 values
horizontal add values in accumulated distance
compute square root of result

With AVX 2, this will let you compute 8 values at a time (per core). With AVX 512 (if you have a new enough processor to support it) you can double that to 16 values at a time (per core).