Caching for complex comparison functions with STL algorithms?

Viewed 97

I have a long to compute function f(T) -> int, and a vector v of T. The task is to find the minimum element after applying f.

Generally, I would use

std::min_element(begin(v), end(v), [&f](auto a, auto b){ return f(a) < f(b); }

Do compilers try to store the computed values (if that makes sense) or do I have to do that by hand? In the second case: Is there a good solution using STL algorithms to do that by hand?

EDIT: Note that the implementation cannot ensure this optimization as it only gets the comparison function and not the structure of that (in this case a unary function).

This is also a more general question because this arises all the time when using STL algorithms. If I measure one example, I only get the insurance for one specific example. I am interested whether the compiler tries to fix this in general (with reasonable optimization enabled). And if not, can you fix this without rewriting the algorithm?

EDIT 2: I think that all questions are well-answered except for the replacement method. This should satisfy that it has the same running time as a (not implemented) function

std::min_element(begin(v), end(v), f);

that stores the value of the last accessed element. Furthermore I would like a solution that is applicable to all algorithms where this optimization can be made.

With c++20, we get the possibility to use projections, but as far as I see, the suggested implementation https://en.cppreference.com/w/cpp/algorithm/ranges/min_element is not optimized for caching (I wonder why, it would not make anything slower, right?).

4 Answers

Do compilers try to store the computed values (if that makes sense) or do I have to do that by hand?

In general this is called memoization and in general the compiler can't do it for you. In specific cases, inlining might allow the optimizer to do something clever.

You can write an automatic memoizing function wrapper if you're expecting to do this a lot - it just needs to keep some (maybe bounded) amount of storage tracking the outputs for previously-encountered inputs. Then you could write something like

auto mf = memoize(f);
std::min_element(begin(v), end(v), [&mf](auto a, auto b){ return mf(a) < mf(b); }

(note that you're still making an explicit decision about how long this particular memo lasts, assuming its cached values are lost when mf goes out of scope).

The language won't do it for you because there are lots of implementation tradeoffs (how much storage is reasonable? what happens if the return value's copy constructor throws an exception when copying it into your memo storage?) that don't have any obvious good default.

In order to check that the compiler was not able to cache the calculated values, I implemented and benchmarked three possibilities:

  • The direct min implemention
  • Min calculation after first calculating the array of f(.) values
  • A specific min_element function, with an unitary function as an argument

Output:

min = 9193
distance = 120
822913 micro-s

New version:
min = 9193
distance = 120
425393 micro-s

With rewriting of min_element function:
min = 9193
distance = 120
416941 micro-s

The second version effectively brings a real speed improvement.

The third version brings a small speed improvement only, but has the advantage to avoid an increase of the memory used.

#include <iostream>
#include <vector>
#include <algorithm>
#include <chrono>

const int param = 10000;
const int N = 10000;

int slowF (int x) {
    int y = x;
    for (int i = 0; i < param; ++i) {
        y = y*y % N;
    }
    return y+2;
}


template<typename T, class ForwardIt, class Funct>
ForwardIt min_element_fct(ForwardIt first, ForwardIt last, Funct f)
{
    if (first == last) return last;
 
    ForwardIt smallest = first;
    T val_smallest = f(*first);
    ++first;
    for (; first != last; ++first) {
        T val;
        if ((val=f(*first)) < val_smallest) {
            smallest = first;
            val_smallest = val;
        }
    }
    return smallest;
}

int main() {
    std::vector<int> v(N);
    v[0] = N/2 - 7;
    for (int i = 1; i < N; ++i) v[i] = (13*v[i-1] + 27) % N;
    
    auto comp = [] (int a, int b) {return slowF(a) < slowF(b);};
    auto t1 = std::chrono::high_resolution_clock::now();
    auto adr_min = std::min_element (v.begin(), v.end(), comp);
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "min = " << *adr_min << std::endl;
    std::cout << "distance = " << std::distance (v.begin(), adr_min) << std::endl;

    auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    std::cout << duration << " micro-s" << std::endl;
    std::cout << "\nNew version: \n";
        
    v[0] = N/2 - 7;
    for (int i = 1; i < N; ++i) v[i] = (13*v[i-1] + 27) % N;
    
    t1 = std::chrono::high_resolution_clock::now();
    std::vector<int> val(N);
    for (int i = 0; i < N; ++i) val[i] = slowF(v[i]);
    
    auto adr_min2 = std::min_element (val.begin(), val.end());
    adr_min = v.begin() + std::distance (val.begin(), adr_min2);
    t2 = std::chrono::high_resolution_clock::now();
    std::cout << "min = " << *adr_min << std::endl;
    std::cout << "distance = " << std::distance (val.begin(), adr_min2) << std::endl;

    duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    std::cout << duration << " micro-s" << std::endl;   
    
    std::cout << "\nWith rewriting of min_element function: \n";
    
    v[0] = N/2 - 7;
    for (int i = 1; i < N; ++i) v[i] = (13*v[i-1] + 27) % N;
    
    t1 = std::chrono::high_resolution_clock::now();
    
    adr_min = min_element_fct<int> (v.begin(), v.end(), slowF);
    t2 = std::chrono::high_resolution_clock::now();
    std::cout << "min = " << *adr_min << std::endl;
    std::cout << "distance = " << std::distance (v.begin(), adr_min) << std::endl;

    duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    std::cout << duration << " micro-s" << std::endl;   

    return 0;
}

The implementation of min_element can't cache f(a), as it works on a predicate. But looking at the larger implementation (compiler plus its standard library), because you pass a lambda the compiler will likely be able to inline the lambda body into min_element. As a result, the optimizer does get a chance to look at that f(a) inside a loop.

A good optimizer might just cache f(*returnvalue) for simple iterator types such as pointers and vector iterators. That would not be specially programmed for min_element, it would be a result of a generic optimization that tries to reuse subexpressions inside loops.

A replacement method is not that easy to do in general without adding some runtime overhead that you would not have for single hand-crafted solutions for min_element, max_element etc. The problem is that for memoizing you always have to know what value you should memoize. "The last accessed" does not really make sense in this context since you will always access two elements for a compare.

So if your memoization only works on the projection (agnostic with regards to how it is used in a comparison) it would need two slots for cached return values, two slots for saving the arguments that create these return values and an indication for which should be overridden next.

And even that may not be enough. I am not sure in which way there sequence points in calls to <, but I think that there aren't any in a call to a comp(_,_) function. So comparing a, b, c can use comp(f(a), f(b)) and comp(f(c), f(a)) and may execute that in order f(a), f(b), f(c), f(a). Here, memoizing with two slots would not be enough.

Having established that we cannot really be agnostic with regards to the comparison function, how about we try it with access to it? Well, then we have another problem. Take this simple implementation of a memoizer. std::min_element always packs the current minimum in the right hand side of the comparison and it stays there if the result is false (standard quote).

#include <type_traits>
#include <optional>
#include <utility>
#include <algorithm>
#include <iostream>

template<class ArgType, class Func, class Comp>
struct Memoizer {
    Memoizer(Func func, Comp comp) : m_func(std::move(func)), m_comp(std::move(comp)) {}

    [[nodiscard]] constexpr bool operator()(const ArgType &lhs, const ArgType &rhs) {
        if (!m_cache) {
            m_cache = std::invoke(m_func, rhs);
        }

        auto temp = std::invoke(m_func, lhs);

        if (m_comp(temp, *m_cache)) {
            m_cache = std::move(temp);
            return true;
        }
        return false;
    }

    using result_t = std::invoke_result_t<Func, ArgType>;
    Func m_func;
    Comp m_comp;

  private:
    std::optional<result_t> m_cache;
};

template<class ArgType, class Func, class Comp>
auto make_memoizer(Func&& f, Comp&& c) -> Memoizer<ArgType, Func, Comp> {
    return {std::forward<Func>(f), std::forward<Comp>(c)};
}

int main() {
    int arr[] = {2,3,1,5,4};

    auto memo = make_memoizer<int>(std::negate{}, std::less{});

    auto min = std::min_element(std::begin(arr), std::end(arr), memo);
    auto max = std::max_element(std::begin(arr), std::end(arr), memo);

    std::cout << "Min: " << -*min << " (Should be: " << -*std::max_element(std::begin(arr), std::end(arr)) << ")\n";
    std::cout << "Max: " << -*max << " (Should be: " << -*std::min_element(std::begin(arr), std::end(arr)) << ")\n";
}

As you can see, max_element does it the other way round and hence we obtain a false result. Now you can template the whole thing again on which side to memoize in which case, but that is just a nicely written invitation for bugs because you will use the wrong variant for the wrong algorithm. Further that just won't work for minmax_element. And at that point writing seperate implementations for min_element and max_element is just easier and safer.

There is still the option to write an iterator adaptor that calls the function on the first *iter and memoizes it for further *iter calls. But I does not seem that the standard gives any guarantees to how the iterators are called so I am not exactly sure if we can guarantee that the function calls are as small as they can be.

This is a very bare bones implementation of such an memoizing iter. You should probably use boost or a similar library to write iterator adaptors cause otherwise they are a pain in the neck with everything they have to foward.

#include <type_traits>
#include <optional>
#include <utility>
#include <algorithm>
#include <iostream>
#include <iterator>

template<class Func, class Iter>
struct MemoIter {
    using result_t = std::invoke_result_t<Func, decltype(*std::declval<Iter>())>;

    using iterator_category = typename std::iterator_traits<Iter>::iterator_category;

    auto operator*() {
        if(!m_cache) {
            m_cache = m_func(*m_iter);
        }
        return *m_cache;
    }

    auto operator++() {
        invalidate_cache();
        ++m_iter;
        return *this;
    }

    auto operator++(int) {
        auto ret = MemoIter(m_iter, invalidate_cache(), m_func);
        ++m_iter;
        return ret;
    }

    MemoIter(Iter iter, Func func) : m_iter(iter), m_func(func) {}
    explicit MemoIter(Iter iter) : MemoIter(iter, {}) {}

    friend bool operator==(const MemoIter& lhs, const MemoIter& rhs) {
        return lhs.m_iter == rhs.m_iter;
    }


    friend bool operator!=(const MemoIter& lhs, const MemoIter& rhs) {
        return lhs.m_iter != rhs.m_iter;
    }


private:
    auto invalidate_cache() {
        auto ret = std::move(m_cache);
        m_cache.reset();
        return ret;
    }
    MemoIter(Iter iter, std::optional<result_t>&& cache, Func func) : m_iter(iter), m_cache(std::move(cache)), m_func(func) {}
    Iter m_iter;
    std::optional<result_t> m_cache;
    Func m_func;

};

int main() {
    int arr[] = {2,3,1,5,4};
    auto begin = MemoIter(std::begin(arr), std::negate<>{});
    auto end = MemoIter(std::end(arr), std::negate<>{});

    auto min = std::min_element(begin, end);
    auto max = std::max_element(begin, end);

    std::cout << "Min: " << *min << " (Should be: " << -*std::max_element(std::begin(arr), std::end(arr)) << ")\n";
    std::cout << "Max: " << *max << " (Should be: " << -*std::min_element(std::begin(arr), std::end(arr)) << ")\n";
}

ADDENDUM: Concerning C++20. The standard explicitly specifies that the number of projection operations is exactly twice the number of comparisons, so the not only the example implementation but also the standard does not give the possibility of caching for this. (Of course, if the compiler could prove that there is no side effect to the projection, it could cache due to the as-if-rule, but that is a big if). As to the why: As you see, memoization is not that easy. You have still further problems if the iterator returns some kind of proxy that may be invalidated in between comparisons. So, just calling the projection more often is more general, easier and also more predictable (you know the number of calls to your projection to be 2*range.size() - 1; if you memoize, it can depend on the order of the elements in your range and whatnot).

Related