std::function has performances issues, how to avoid it?

Question

I have classes which allow to compound covariance functions (also called Kernels see https://stats.stackexchange.com/questions/228552/covariance-functions-or-kernels-what-exactly-are-they) and then compute the covariance given the new kernel for example :

auto C = GaussianKernel(50,60) + GaussianKernel(100,200);
auto result = C.covarianceFunction(30.0,40.0);

But The problem is that I call a std::function when I want to compute the covariance, Is there a simple way to avoid it ?
Note that I want to compute a big covariance matrix (approximately 50K*50K) which means that performances matters.

Here is the code

class Kernel {
public: 
    /*
    Covariance function : return the covariance between two R.V. for the entire kernel's domain definition. 
    */
    virtual double covarianceFunction(
        double   X,
        double   Y
    )const = 0 ;
    ~Kernel() = default;
};

class FooKernel : public Kernel {
public:
    FooKernel(std::function<double(double, double)> fun) : fun_(fun) {}
    double covarianceFunction(
        double   X,
        double   Y
    ) const {
        return fun_(X, Y);
    }
    template<class T>
    auto operator+(const T b) const {
        return FooKernel([b, this](double X, double Y) -> double {
            return this->covarianceFunction(X, Y) + b.covarianceFunction(X, Y);
        });
    }
    FooKernel operator=(const FooKernel other) const {
        return other;
    }
private:
    std::function<double(double, double)> fun_;
};

class GaussianKernel : public Kernel {
public:
    GaussianKernel(double sigma, double scale) : m_sigma(sigma), m_scale(scale) {}
    GaussianKernel(double sigma) : m_sigma(sigma), m_scale(1) {}
    /*
    A well known covariance function that enforces smooth deformations
    Ref : Shape modeling using Gaussian process Morphable Models, Luethi et al.
    */
    double covarianceFunction(
        double   X,
        double   Y
    ) const 
    {
        //use diagonal matrix
    doulbe result;
    result = m_scale  *  exp(-std::norm(X - Y) / (m_sigma*m_sigma));
    return result;      
    }
    template<class T>
    auto operator+(const T b) const {
        return FooKernel([b, this](double X, double Y) -> double {
            auto debugBval = b.covarianceFunction(X, Y);
            auto debugAval = this->covarianceFunction(X, Y);
            auto test = debugBval + debugAval;
            return test;
        });
    }
private:
    double m_sigma;
    double m_scale;
};

Answer 1

by templating FooKernel you can avoid the need for std::function.

#include <iostream>
#include <complex>
#include <functional>


class Kernel {
public: 
    /*
    Covariance function : return the covariance between two R.V. for the entire kernel's domain definition. 
    */
    virtual double covarianceFunction(
        double   X,
        double   Y
    )const = 0 ;
    ~Kernel() = default;
};


template <typename Func>
class FooKernel : public Kernel {
public:

    FooKernel(Func&& fun) : fun_(std::forward<Func>(fun)) {}
    double covarianceFunction(
        double   X,
        double   Y
    ) const {
        return fun_(X, Y);
    }
    template<class T>
    auto operator+(const T b) const {
        return make_foo_kernel([b, this](double X, double Y) -> double {
            return this->covarianceFunction(X, Y) + b.covarianceFunction(X, Y);
        });
    }
    FooKernel operator=(const FooKernel other) const {
        return other;
    }
private:
   Func fun_;
};

template <typename Func>
auto make_foo_kernel(Func&& fun)
{
    return FooKernel<Func>(std::forward<Func>(fun));
}


class GaussianKernel : public Kernel {
public:
    GaussianKernel(double sigma, double scale) : m_sigma(sigma), m_scale(scale) {}
    GaussianKernel(double sigma) : m_sigma(sigma), m_scale(1) {}
    /*
    A well known covariance function that enforces smooth deformations
    Ref : Shape modeling using Gaussian process Morphable Models, Luethi et al.
    */
    double covarianceFunction(
        double   X,
        double   Y
    ) const 
    {
        //use diagonal matrix
    double result;
    result = m_scale  *  exp(-std::norm(X - Y) / (m_sigma*m_sigma));
    return result;      
    }
    template<class T>
    auto operator+(const T b) const {
        return make_foo_kernel([b, this](double X, double Y) -> double {
            auto debugBval = b.covarianceFunction(X, Y);
            auto debugAval = this->covarianceFunction(X, Y);
            auto test = debugBval + debugAval;
            return test;
        });
    }
private:
    double m_sigma;
    double m_scale;
};

int main()
{
    auto C = GaussianKernel(50,60) + GaussianKernel(100,200);
    auto result = C.covarianceFunction(30.0,40.0);

    return 0;
}

Demo

Answer 2

With this design, the only improvement over using std::function is to template-parameterize your classes, which can create other undesired problems.

template<class Fun>
class FooKernel : public Kernel {
public:
    FooKernel(Fun&& fun) : fun_(std::forward<Fun>(fun)) {}
...
private:
    Fun fun_;
};

If you don't want to template your class and If you need your classes to own a stateful function object, pretty much std::function is the only way to go.

However, if you don't need ownership or the the function or the function object is stateless (e.g. a free function) and you state that in your question I can give you an alternative option.

Answer 3

As you said you liked the clarity of std::function you could try this non-owning function reference class:

#include <utility>

template<typename TSignature> class function_ref;

template<typename TRet, typename ...TParams>
class function_ref<TRet(TParams...)> final
{
    using refptr_t = void*;
    using callback_t = TRet (*)(refptr_t, TParams&&...);

    callback_t m_callback = nullptr;
    refptr_t m_callable = nullptr;

public:
    constexpr function_ref() noexcept = default;
    constexpr function_ref(const function_ref&) noexcept = default;
    constexpr function_ref& operator=(const function_ref&) noexcept = default;
    constexpr function_ref(function_ref&&) noexcept = default;
    constexpr function_ref& operator=(function_ref&&) noexcept = default;
    ~function_ref() noexcept = default;

    template <
        typename T,
        typename = typename std::enable_if_t<
            std::is_invocable_r_v<TRet, T(TParams...), TParams...> &&
            !std::is_convertible_v<std::decay_t<T>, function_ref>
        >
    >
    constexpr function_ref(T &&_callable) noexcept :
        m_callback(
            [](refptr_t callable, TParams&& ...params)
            {return (*reinterpret_cast<std::remove_reference_t<T>*>(callable))(std::forward<TParams>(params)...);}
            ),
        m_callable(reinterpret_cast<refptr_t>(std::addressof(_callable))) 
    {}

    constexpr decltype(auto) operator()(TParams&& ...params) noexcept 
    {
        return m_callback(m_callable, std::forward<TParams>(params)...);
    }

    constexpr operator bool() noexcept { return m_callback; }
};

This doesn't have the overhead of std::function as it doesn't need to own the callable, and with my tests it is usually completely inlined with -O3 optimisations. This is my modified implementation of the class discussed by Vittorio Romeo in this talk. You still need to watch the lifetimes of the functions you pass to the constructor, however it is perfect to take function parameters.

Example usage:

void func(int x)
{
    std::cout<<x<< " I'm a free func!\n";
}

class Obj 
{
    public: 
    void member(int x) { std::cout<<x<< " I'm a member func!\n";}
};

int main()
{
    // Define the signature
    using func_ref_t = function_ref<void(int)>;

    // Can be used with stateful lambdas
    int bar = 1;
    auto lambda = [&bar](int x){std::cout<<x<< " I'm a lambda!\n"; ++bar;};

    // Copy and move
    func_ref_t lref(lambda);
    auto cpy = lref;
    auto mv = std::move(lref);
    cpy(1);
    mv(2);

    // See the modified var from the lambda
    std::cout<<bar<<'\n';

    // Use with free functions
    auto fref = func_ref_t{func};
    fref(4);

    // We can wrap member functions with stateful lamdas
    Obj obj;
    auto mem = [&obj](int x) { obj.member(x); };

    auto mref = func_ref_t{mem};
    mref(5);

}