Matrix : joined ops - blas

[pfeatherstone]

In dlib's expression template suite, there isn't an overload of operator* for joined matrices.
This could be an optimization when performing dot products with circular buffers for example, which are essentially two joined arrays. The optimization would be to split the dot product into a sum of two dot products. You could then use a BLAS routine. At the moment this isn't possible.

Here is some example code:

#include <vector>
#include <algorithm>
#include <numeric>
#include <cstdio>
#include <cassert>
#include <cstdint>
#include <chrono>

#include <dlib/sliding_buffer.h>
#include <dlib/matrix.h>
#include <dlib/matrix/matrix_mat.h>

using namespace std::chrono;
using namespace dlib;

namespace p
{
    template<typename T>
    class circular_buffer
    {
    public:
        circular_buffer(
            size_t size = 1
        ) : buf(size, T(0))
        {
        }

        void reset(size_t size)
        {
            buf.resize(size, T(0));
            offset = 0;
            count  = 0;
        }

        void push(T x)
        {
            buf[offset] = x;
            offset = (offset + 1) % buf.size();
            count++;
        }

        size_t size() const
        {
            return buf.size();
        }

        template<typename It>
        auto dot_product (
            It hbeg,
            It hend
        ) const
        {
            const auto hlen = std::distance(hbeg, hend);
            assert(hlen == buf.size() || 0 != "filter_size != buf.size()");

            const auto sum_oldest = std::inner_product(beg(),  mid1(), hbeg,          T{0});
            const auto sum_newest = std::inner_product(mid2(), end(),  hend - offset, T{0}); 
            return sum_newest + sum_oldest;
        }

        template<typename Container>
        auto dot_product (
            const Container& y
        ) const
        {
            return dot_product(y.begin(), y.end());
        }

        template<typename Container>
        auto dot_product3 (
            const Container& y
        ) const
        {
            return dot(join_cols(mat(buf.data() + offset, buf.size() - offset),
                                 mat(buf.data(),          offset)),
                       mat(y));
        }

        template<typename Container>
        auto dot_product_reversed (
            const Container& y
        ) const
        {
            return dot_product(y.rbegin(), y.rend());
        }

        template<typename Container>
        auto dot_product_reversed3 (
            const Container& y
        ) const
        {
            return dot(flip(join_cols(mat(buf.data() + offset, buf.size() - offset),
                                      mat(buf.data(),          offset))),
                       mat(y));
        }

    private:

        const auto beg()    const {return buf.begin() + offset;}
        const auto mid1()   const {return buf.end();}
        const auto mid2()   const {return buf.begin();}
        const auto end()    const {return beg();}
        
        std::vector<T>  buf;
        size_t          offset = 0;
        size_t          count  = 0;
    };
}

const auto time_it = [](const size_t ntests, auto func)
{
    const auto start = high_resolution_clock::now();
    auto ret{0};

    for (size_t i = 0 ; i < ntests ; ++i)
        ret = func();

    const auto stop = high_resolution_clock::now(); 
    return std::make_tuple(ret, (duration_cast<nanoseconds>(stop - start).count() * 1e-3) / ntests);
};                                                       

int main()
{
    p::circular_buffer<float> a1(100);
    dlib::circular_buffer<float> a2(100);

    for (int i = 0 ; i < 160 ; ++i)
    {
        a1.push(i);
        a2.push_back(i);
    }

    std::vector<float> b(a1.size());
    std::iota(b.begin(), b.end(), 1.0);
    int ret = 1;
    double delta_us;

    size_t LOOP_LIMIT = 10000000;

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product(b); });
    printf("STL      dot       : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product_reversed(b); });
    printf("STL      rdot      : ret = %i , computed in : %f us\n", ret, delta_us);

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product3(b); });
    printf("dlib mat dot       : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product_reversed3(b); });
    printf("dlib mat rdot      : ret = %i , computed in : %f us\n", ret, delta_us);

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(mat(a2), mat(b)); });
    printf("dlib circ buf dot  : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(flip(mat(a2)), mat(b)); });
    printf("dlib circ buf rdot : ret = %i , computed in : %f us\n", ret, delta_us);
}

Also, for whatever reason, when using dlib::circular_buffer directly, it's nearly 10X slower...

Here is the output on my machine:

STL      dot       : ret = 636300 , computed in : 0.050482 us
STL      rdot      : ret = 469650 , computed in : 0.055836 us
dlib mat dot       : ret = 636300 , computed in : 0.082418 us
dlib mat rdot      : ret = 469650 , computed in : 0.085437 us
dlib circ buf dot  : ret = 636300 , computed in : 0.557763 us
dlib circ buf rdot : ret = 469650 , computed in : 0.573191 us

Note, I have assert statements disabled and MKL enabled. Even so, BLAS functions aren't used anywhere above.

[davisking]

BLAS usually doesn't get you much for dot products. Although depends how big they are. But yeah, some overload for operator* for mat(some_circular_buffer) would be cool. The more overloads that make things faster the better :)

[pfeatherstone]

By the way, replacing:

return data[(i+offset)%data.size()]; 

with

const auto index = i + offset;
return data[index >= data.size() ? index % data.size() : index];

in dlib::circular_buffer<T>::operator[] gave me

dlib  dot  : ret = 636300 , computed in : 0.343111 us
dlib  rdot : ret = 469650 , computed in : 0.342680 us

That's a 40% improvement. That's using a trick from Chandler Carruth talk at CppCon.

[pfeatherstone]

That's a quick win but we are nowhere near manually calling std::inner_product twice on two sub-arrays.

[davisking]

Yeah nice.

[pfeatherstone]

I've added dlib::mat() overload for p::circular_buffer

Update code:

#include <vector>
#include <algorithm>
#include <numeric>
#include <cstdio>
#include <cassert>
#include <cstdint>
#include <chrono>

#include <dlib/sliding_buffer.h>
#include <dlib/matrix.h>
#include <dlib/matrix/matrix_mat.h>

using namespace std::chrono;
using namespace dlib;

namespace p
{
    template<typename T>
    class circular_buffer
    {
    public:
        using value_type = T;

        circular_buffer(
            size_t size = 1
        ) : buf(size, T(0))
        {
        }

        void reset(size_t size)
        {
            buf.resize(size, T(0));
            offset = 0;
            count  = 0;
        }

        void push(T x)
        {
            buf[offset] = x;
            offset = (offset + 1) % buf.size();
            count++;
        }

        size_t size() const
        {
            return buf.size();
        }

        const T& operator[](size_t i) const
        {
            const auto index = i + offset;
            return buf[index >= buf.size() ? index % buf.size() : index];
 	    //return buf[(i+offset)%buf.size()];
        }

        template<typename It>
        auto dot_product (
            It hbeg,
            It hend
        ) const
        {
            const auto hlen = std::distance(hbeg, hend);
            assert(hlen == buf.size() || 0 != "filter_size != buf.size()");

            const auto sum_oldest = std::inner_product(beg(),  mid1(), hbeg,          T{0});
            const auto sum_newest = std::inner_product(mid2(), end(),  hend - offset, T{0}); 
            return sum_newest + sum_oldest;
        }

        template<typename Container>
        auto dot_product (
            const Container& y
        ) const
        {
            return dot_product(y.begin(), y.end());
        }

        template<typename Container>
        auto dot_product2 (
            const Container& y
        ) const
        {
            return dot(join_cols(mat(buf.data() + offset, buf.size() - offset),
                                 mat(buf.data(),          offset)),
                       mat(y));
        }

        template<typename Container>
        auto dot_product_reversed (
            const Container& y
        ) const
        {
            return dot_product(y.rbegin(), y.rend());
        }

        template<typename Container>
        auto dot_product_reversed2 (
            const Container& y
        ) const
        {
            return dot(flip(join_cols(mat(buf.data() + offset, buf.size() - offset),
                                      mat(buf.data(),          offset))),
                       mat(y));
        }

    private:

        const auto beg()    const {return buf.begin() + offset;}
        const auto mid1()   const {return buf.end();}
        const auto mid2()   const {return buf.begin();}
        const auto end()    const {return beg();}
        
        std::vector<T>  buf;
        size_t          offset = 0;
        size_t          count  = 0;
    };

    template <class T>
    auto mat (
        const circular_buffer<T>& m 
    )
    {
        using namespace dlib;
        using op = op_std_vect_to_mat<circular_buffer<T>>;
        return matrix_op<op>(op(m));
    }
}

const auto time_it = [](const size_t ntests, auto func)
{
    const auto start = high_resolution_clock::now();
    auto ret{0};

    for (size_t i = 0 ; i < ntests ; ++i)
        ret = func();

    const auto stop = high_resolution_clock::now(); 
    return std::make_tuple(ret, (duration_cast<nanoseconds>(stop - start).count() * 1e-3) / ntests);
};                                                       

int main()
{
    using dlib::mat;
    using dlib::flip;
    using dlib::dot;

    p::circular_buffer<float> a1(100);
    dlib::circular_buffer<float> a2(100);

    for (int i = 0 ; i < 160 ; ++i)
    {
        a1.push(i);
        a2.push_back(i);
    }

    std::vector<float> b(a1.size());
    std::iota(b.begin(), b.end(), 1.0);
    int ret = 1;
    double delta_us;

    size_t LOOP_LIMIT = 10000000;

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product(b); });
    printf("STL      dot       : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product_reversed(b); });
    printf("STL      rdot      : ret = %i , computed in : %f us\n", ret, delta_us);

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product2(b); });
    printf("dlib concat dot    : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return a1.dot_product_reversed2(b); });
    printf("dlib concat rdot   : ret = %i , computed in : %f us\n", ret, delta_us);

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(mat(a1), mat(b)); });
    printf("dlib mat dot       : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(flip(mat(a1)), mat(b)); });
    printf("dlib mat rdot      : ret = %i , computed in : %f us\n", ret, delta_us);

    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(mat(a2), mat(b)); });
    printf("dlib circ buf dot  : ret = %i , computed in : %f us\n", ret, delta_us);
    std::tie(ret, delta_us) = time_it(LOOP_LIMIT, [&]{ return dot(flip(mat(a2)), mat(b)); });
    printf("dlib circ buf rdot : ret = %i , computed in : %f us\n", ret, delta_us);
}

I've compiled the above with -O3 -march=native -ffast-math using gcc 12 and I made the suggested change to dlib::circular_buffer<T>::operator[].
I get the following results:

STL      dot       : ret = 636300 , computed in : 0.007123 us
STL      rdot      : ret = 469650 , computed in : 0.008296 us
dlib concat dot    : ret = 636300 , computed in : 0.036760 us
dlib concat rdot   : ret = 469650 , computed in : 0.039871 us
dlib mat dot       : ret = 636300 , computed in : 0.085077 us
dlib mat rdot      : ret = 469650 , computed in : 0.078495 us
dlib circ buf dot  : ret = 636300 , computed in : 0.344959 us
dlib circ buf rdot : ret = 469650 , computed in : 0.347123 us

Fundamentally, I don't understand why dot products on mat(p::circular_buffer) are a lot faster than on mat(dlib::circular_buffer). They are doing the exact same thing now and DLIB_ASSERT is disabled.

[davisking]

Yeah I was confused about that. Superficially it looks like what you pasted above is identical but the timings are way different. You 100% sure the assert is disabled?

[pfeatherstone]

Yeah certain. So I changed operator[] for both p::circular_buffer and dlib::circular_buffer to:

auto index = i + offset;
if (index >= data.size())
    index -= data.size();
return data[index];

while keeping the DLIB_ASSERT statements in dlib::circular_buffer<T>::operator[]
and I get the results:

STL      dot       : ret = 636300 , computed in : 0.007418 us
STL      rdot      : ret = 469650 , computed in : 0.007805 us
dlib mat dot       : ret = 636300 , computed in : 0.034434 us
dlib mat rdot      : ret = 469650 , computed in : 0.037156 us
dlib circ buf dot  : ret = 636300 , computed in : 0.038685 us
dlib circ buf rdot : ret = 469650 , computed in : 0.041538 us

So that's equivalent behaviour now. But in both cases, we are way off.