Expression template to compute the Euclidean distance

Question

I was writing some geometry-related code again and had a closer look at my function supposed to compute the Euclidean distance between two points (N-dimensional points by the way, hence the N template parameter). Here is a simplified version:

template<std::size_t N, typename T>
auto distance(const Point<N, T>& lhs, const Point<N, T>& rhs)
    -> T
{
    T res{};
    for (std::size_t i = 0 ; i < N ; ++i)
    {
        auto tmp = std::abs(lhs[i] - rhs[i]);
        res += tmp * tmp;
    }
    return std::sqrt(res);
}

So far, so good. However, one very common operation is to compare the distances. Generally speaking, when comparing the distances, the sqrt is optimized away and the sum of the squares is compared instead of the distance itself. Therefore, I tried to create some kind of expression template to represent the distance between two points, so that users will benefit from both the ease of use and the "get rid of sqrt optimization" when comparing distances. Basically, the call of sqrt is not done until the exact value of the distance is needed. Here is the class:

template<typename T>
struct DistanceExpression
{
    explicit constexpr DistanceExpression(T data):
        _data(data)
    {}

    operator T() const
    {
        return std::sqrt(_data);
    }

    bool operator==(const DistanceExpression& other) const
    {
        return _data == other._data;
    }

    bool operator!=(const DistanceExpression& other) const
    {
        return !(*this == other);
    }

private:

    T _data;
};

My new distance function is implemented as such:

template<std::size_t N, typename T>
auto distance(const Point<N, T>& lhs, const Point<N, T>& rhs)
    -> DistanceExpression<T>
{
    T res{};
    for (std::size_t i = 0 ; i < N ; ++i)
    {
        auto tmp = std::abs(lhs[i] - rhs[i]);
        res += tmp * tmp;
    }
    return DistanceExpression<T>{res};
}

Here is a minimal working code at Coliru. Is such a design reasonable or is it overkill to elegantly solve this problem?

Answer 1

I find all this an overkill. Agreeing in general with Michael Urman, the programmer/user should be aware of such choices. I can clearly see at least two different functions here:

template<std::size_t N, typename T>
T squared_distance(const Point<N, T>& lhs, const Point<N, T>& rhs)
{
    T res{};
    for (std::size_t i = 0 ; i < N ; ++i)
    {
        auto tmp = lhs[i] - rhs[i];
        res += tmp * tmp;
    }
    return res;
}

template<std::size_t N, typename T>
T distance(const Point<N, T>& lhs, const Point<N, T>& rhs)
{
    return std::sqrt(squared_distance(lhs, rhs));
}

There are so many more useful cases for introducing expression templates, so why an extra burden for the compiler just for this?

For instance, given Point's x, y, vector subtraction

z = x - y

could be an expression template, so that z[i] == x[i] - y[i] for each i. Squaring each element would be another expression template, so that

z = square(x - y)

would have z[i] == w * w for each i, where w = (x - y)[i]. Then,

squared_norm(x - y)

would compute the actual reduction as

template<typename X>
T sum(const X& x) { return std::accumulate(x.begin(), x.end(), X{}); }

template<typename X>
T squared_norm(const X& x) { return sum(square(x)); }

assuming Point or a more general expression template are equipped with begin(), end() (which they should). Note I am not using std::inner_product to avoid the extra cost when x is a (lazy) expression template. Then, squared_distance would be trivial

template<typename X, typename Y>
T squared_distance(const X& x, const Y& y) { return squared_norm(x - y); }

or not needed at all, while distance generalized as

template<typename X, typename Y>
T distance(const X& x, const Y& y) { return std::sqrt(squared_distance(x, y)); }

There's more to be generalized by considering rvalue references and forwading as appropriate, but I wanted to keep this clean.

The user should explicitly use squared Euclidean norms and distances when needed. For instance, consider this scientific article on nearest-neighbor search:

In (1), the distance measures d1 and d2 in R are induced by d, so that for all a, b : d(a, b) = d1(a1, b1) + d2(a2, b2). The simplest and most important case is setting d, d1, and d2 to be squared Euclidean distances in respective spaces.

In other words, we are making use of the separability of the squared Euclidean distance in product spaces. The entire algorithm is based on this choice, and your design would not help in this case.

PS. Why would you need std::abs() at all? And why not just std::array<T, N> instead of Point<N, T>?

Answer 2

My concern now is that you make multiple calls to get the actual distance.
Every-time you do that then you incur the expense of the sqrt() operation.

I suppose it depends on the exact use case of your application which happens more often. But if this is for general code that could be used a lot in either capacity. I would add another two fields to indicate if it was dirty or not and cache the result so you don't have to recompute the actual value each time.

But that has the side affect of adding a conditional branch into the get which can also be expensive when looking at micro optimizations.

Answer 3

I agree that more context is necessary to know whether this is the right approach. When I've worked with ray tracing, you more often needed to know the closest item, or that a distance was above or below some threshold, but did not need to know what the actual distance was. In those contexts, at the time, I worked with raw numbers that were explicitly the square of the distance, returned from a function called distance2; I think I prefer that approach.

It looks like you're trying to hide the squares and square roots. Abstractions like that can be very useful, but come at a cost. Without the abstraction, the programmer has to track the complexity, and this means fewer other things will fit in the programmer's head at the same time. When hiding this optimization, however, you either have to make another time vs. space trade-off (for instance, you initially save time by delaying or avoiding the sqrt call, but then have to choose whether to cache the results; caching saves time but costs space). Essentially the class is trying to guess what the programmer needs, and this often goes badly.

If memory use is a concern, you can mitigate the extra space requirements of a cache by adding complexity. Relying on the fact that no distances are negative; using a negative distance can mark whether it is the distance or its square. But despite saving the extra storage, that complexity shows up in all operators.

How far should this abstraction go? Should you have specified comparison operators such as operator<? If so, should you specify them both for two DistanceExpression<T> values as well as for a T and a DistanceExpression<T> value? This could be useful if you want to write if (distance(p, q) < 10) without requiring the sqrt call. But again if (distance2(p, q) < 100) is also quite clear.

Reviewing beyond your explicit question, good call on marking DistanceExpression<T>'s constructor explicit, although was surprised that you initialized its _data member using parenthesis () instead of curly braces {}. Old habits die hard! I also prefer avoiding the leading underscore on member names, but it appears that only when the next letter is capital is it categorically reserved. Beware thin ice.