Performance: repeated sweep through large binary dataset

Question

As part of a bigger mcmc sampling algorithm, I'm repeatedly looping through my entire dataset. When the data grows large this becomes a pretty large bottleneck in my algorithm and I'm looking for ways to make it faster.

Some things I'm considering are:

• What are the most efficient datatypes? For example, for the variable x I use std::vector< std::vector<bool> > because data are binary; 0 or 1. However, this leads to many type conversions. Is there a better way?
• Can I use iterators to loop over the data more efficiently?

All suggestions are welcome!

Edit: the performance critical function does the following:

The goal of the code is to simulate a dataset of size n_p by n_i, while simultaneously evaluating a likelihood ratio. The simulated new data is binary, and created by rlogis(0.0, 1.0) <= val ? 1 : 0. The subsequent evaluation of the ratio happens via acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];. These two steps seem optimal from a performance perspective to me; my question is primarily about the surrounding structure.

Code:

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

using namespace std; // <-- only here for the example; I don't normally do this

typedef unsigned int uint;
typedef std::vector<double> vec;
typedef std::vector<vec> mat;

// will be replaced by actual random number generator
double rlogis() {
    return rand() % 5 - 2;
}


// improving performance of this function is my goal
double orig(const uint np, const uint ni, const uint nd,
               const mat& theta, const mat& C,
               const vec& delta, const vec& lambdaNew,
               const int d, std::vector<std::vector<bool> >& x) {

    double acc = 0.0;
    bool xnew;
    vec tmp(nd);
    for (uint j = 0; j < ni; ++j) {

        for (uint d2 = 0; d2 < nd; ++d2) {
            tmp[d2] = lambdaNew[d2] * C[j][d2];
        }

        for (uint i = 0; i < np; ++i) {

            double val = delta[j];
            for (uint d2 = 0; d2 < nd; ++d2) {
                val += tmp[d2] * theta[i][d2];
            }

            xnew = rlogis() <= val;
            acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];

        }
    }
    return acc;
}

// wrapper function for timing
int main() {

    // dimensions of the data
    const uint np = 2E4; // this could be larger, easily 1E5
    const uint ni = 200; // this stays constant
    const uint nd = 3;   // could be 4-5, but not 10

    // objects
    std::vector<double> delta(ni);
    std::vector<double> lambdaNew(nd);
    vector< std::vector<double> > theta(np, std::vector<double>(nd));
    vector< std::vector<double> > C(ni, std::vector<double>(nd));
    std::vector< std::vector<bool> > x(np, std::vector<bool>(ni));

    // fill objects with values (exact values are irrelevant)
    srand(1);
    int d = rand() % nd;
    for (uint i = 0; i < np; ++i) {
        for (uint d = 0; d < nd; ++d) {
            theta[i][d] = 2 * rand() - 1;
        }
        for (uint j = 0; j < ni; ++j) {
            x[i][j] = bool(int(rand() + .5));
        }
    }
    for (uint j = 0; j < ni; ++j) {
        delta[j] = 2 * rand() - 1;
        for (uint d = 0; d < nd; ++d) {
            C[j][d] = 2 * rand() - 1;
        }
    }
    for (uint d = 0; d < nd; ++d) {
        lambdaNew[d] = 5 * rand() + 1;
    }

    // start timing
    {
        srand(2);
        auto start = chrono::steady_clock::now();
        auto ans = orig(np, ni, nd, theta, C, delta, lambdaNew, d, x);
        auto time1 = chrono::steady_clock::now() - start;

        cout << "Orig: " << ans << "\n"
             << "Time: " << chrono::duration <double, milli> (time1).count() << " ms\n" << endl;
    }

    return 0;
}

Answer 1

You did not provide any performance measurements, so every answer here is basically guesswork, but I suspect that the lack of speed you are experiencing could have something to do with std::vector<bool>.

The problem is that std::vector<bool> is not really what it appears to be. While a vector of any other type actually works in units of bytes, the specialization for bool works with individual bits, packing them as tightly as possible. This is, of course, great if you want to reduce memory usage, but it also is really slow since you have to mask and shift bits around every time you want to access a value.

For this reason, in combination with the fact that it doesn't behave like a normal vector in certain cases (which are, sadly, not exactly rare), std::vector<bool> is mostly seen as a design flaw in the standard library and does not see much use.

You should try to replace std::vector<bool> with a similar fitting type, e.g. std::vector<char>, and measure that to determine whether and by how much this alleviates your issue.

If you are not convinced yet, here is a quick-and-dirty benchmark to show just how much slower std::vector<bool> is.

Answer 2

I'll compact the code for easier reading, and add the const values:

    // loop1, ni = 200; this stays constant
    for (uint j = 0; j < ni; ++j) {

        // loop2, nd = 3-5
        for (uint d2 = 0; d2 < nd; ++d2) tmp[d2] = lambdaNew[d2] * C[j][d2];

        //loop3 np = 2E4, easily 1E5
        for (uint i = 0; i < np; ++i) {

            //loop4, nd = 3-5
            double val = delta[j];
            for (uint d2 = 0; d2 < nd; ++d2) val += tmp[d2] * theta[i][d2];

            xnew = rlogis() <= val;
            acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];    
        }
    }

This will be memory-intensive, so here are object sizes:

double delta[200];
double lambdaNew[3...5]
double C[200][3...5]

double theta[2e4...1e5][3...5]    4.0 MB
bool   x[2e4...1e5][200]          2.5 MB

Since theta is read 200 times we really want it to fit in cache (otherwise it will have to be fetched from RAM every time) but its large size may be a problem.

Additionally, theta is a vector of vectors and I believe the implementation will allocate each line of the matrix separately with a new[], which will add some headers and pad it to the nearest legal block size. This is problematic here as the lines are very small (3-5 doubles) so the heap memory management overhead risks being significant. See here for some numbers.

DDR RAM is optimized for high throughput burst or sequential reads, and utterly sucks at random accesses. What you want is read your matrix sequentially, but that won't happen here, instead the cpu will read each line, then skip to the next pointer. Most likely due to prefetching etc, accesses will still be sequential, but perhaps half of the RAM bandwidth will be wasted on useless stuff like headers and padding. Worse, because cache memory works by cache lines, precious cache will be wasted on the headers and padding. Cache is vital, you don't want to fill it with useless stuff, instead you want to make maximum use of it.

Therefore:

• Use a proper matrix allocated as a single contiguous block of memory which will be read sequentially without any holes or jumps or messing with pointer indirection
• Use floats instead of doubles if you can, if that makes it fit into L2 cache then it will be much faster.

Now, you want loop4 to be unrolled and use proper SIMD vector instructions. If nd is a known compile time constant, template the function and put it as a parameter, then check your compiler flags and make sure the compiler unrolls the loop properly by looking at the disassembled output. I believe the newest flavor of SSE has a "dot product" instruction which would be perfect here. If nd can vary in a tiny range like 3..5, it would probably be worth it to compile templated versions with nd=3,4,5 with the loops unrolled. This should be really fast. Again check if you can use floats instead of doubles.

This should turn loop4 into a ridiculously low number of SSE instructions.

Then you should profile, and your next bottleneck will probably be the random number generator in rlogis. Try to find a fast one. You should check its period. The period of rand() may not be enough for your needs. If you're a sneaky bastard, you can have another core generate the random numbers and use a FIFO, but this needs a bit of skill to make sure the synchronization overhead doesn't negate the gains.

Now, the boolean array.

Since it is accessed as x[i][j] but i is the innermost index, you should transpose it and put the innermost loop index last as in x[j][i]. This is important as the array is large, and the [j] goes to 200 so your current code will not access memory sequentially, rather it will make random skips.

Random DDR RAM access cost is HUGE! It can be up to 100ns. L1-L2 caches take 1-3 ns. 100ns on a 3GHz cpu is 300 cycles! Your CPU can run several instructions per cycle... One random RAM access can cost as much as 1000 instructions or more if you are really unlucky... sequential accesses are your friend!

Needless to say, a bool array with one value per bit needs a few bit twiddling operations per access, so using one byte per value is faster... unless the now 8x larger array no longer fits in cache. If DDR RAM needs to be accessed, then it will be a lot slower than twiddling bits.

loop2 can be unrolled also.

loop3 could be parallelized with OpenMP but it is so short I'm not sure there will be any gains. That depends how costly rlogis() is. You can always try and benchmark.

The whole thing could also run on a GPU. In this case loop3 can be parallelized to the max.

Remember stuff like using iterators or indexes can matter a little bit, but memory access paterns, caching, random/sequential, etc, are the real issues in your scenario. Just transposing a matrix can have a huge impact on performance if the memory accesses become sequential instead of random/skips.

DDR4 may give you 60 GB/s sequential access... but random byte access is 10M/s. Ouch.

Answer 3

Your interface is a real mess between C and C++. Why do you provide ni, nd and np when they are actually sizes of the respective arrays? Also you should really not use std::vector<bool> especially not if you use it as a fake integer

That leads to the following:

double orig(const mat& theta, const mat& C,
            const vec& delta, const vec& lambdaNew,
            const int d, std::vector<std::vector<int> >& x) {

    double acc = 0.0;
    bool xnew;
    vec tmp(lambdaNew.size());
    for (std::size_t j = 0; j < C.size(); ++j) {

        for (std::size_t d2 = 0; d2 < C[j].size(); ++d2) {
            tmp[d2] = lambdaNew[d2] * C[j][d2];
        }

        for (std::size_t i = 0; i < theta.size(); ++i) {    
            double val = delta[j];
            for (uint d2 = 0; d2 < nd; ++d2) {
                val += tmp[d2] * theta[i][d2];
            }

            xnew = rlogis() <= val;
            acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];    
        }
    }
    return acc;
}

Why is xnew defined outside the loop it is used in although you never use it anymore?

The next thing to notice, is that tmp is actually more of a copy of lambdaNew. One could be tempted to just write vec tmp = lambdaNew; rather than vec tmp(lambdaNew.size());, however that prevents a possible optimization by the compiler if lambdaNew is passed in a move constructible way to the function. Therfore it is better to let the compiler do the copy

double orig(const mat& theta, const mat& C,
            const vec& delta, vec lambdaNew,
            const int d, std::vector<std::vector<int> >& x) {

    double acc = 0.0;
    for (std::size_t j = 0; j < C.size(); ++j) {

        for (std::size_t d2 = 0; d2 < C[j].size(); ++d2) {
            lambdaNew[d2] *= C[j][d2];
        }

        for (std::size_t i = 0; i < theta.size(); ++i) {    
            double val = delta[j];
            for (uint d2 = 0; d2 < nd; ++d2) {
                val += lambdaNew[d2] * theta[i][d2];
            }    

            bool xnew = rlogis() <= val;
            acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];    
        }
    }
    return acc;
}

If you want to do the multiplication i a C++ way you can use std::transform for it.

for (std::size_t d2 = 0; d2 < C[j].size(); ++d2) {
    lambdaNew[d2] *= C[j][d2];
}

becomes

std::transform(lambdaNew.begin(), 
               lambdaNew.end(), 
               C[j].begin(), 
               lambdaNew.begin(),
               std::multiplies<double>());

Not that its pretty but if you really like STL... More interesting are other STL functions for example std::inner_product. It does what it says, by adding the elementwise products of two ranges to an initial value. So this piece of code:

double val = delta[j];
for (uint d2 = 0; d2 < nd; ++d2) {
    val += lambdaNew[d2] * theta[i][d2];
}

becomes:

double val = std::inner_product(lambdaNew.begin(),
                                lambdaNew.end(),
                                theta[i].begin(),
                                delta[j]);

Now together we have this

double orig(const mat& theta, const mat& C,
            const vec& delta, vec lambdaNew,
            const int d, std::vector<std::vector<int> >& x) {

    double acc = 0.0;
    for (std::size_t j = 0; j < C.size(); ++j) {
        std::transform(lambdaNew.begin(), 
                       lambdaNew.end(), 
                       C[j].begin(),        
                       lambdaNew.begin(),
                       std::multiplies<double>());

        for (std::size_t i = 0; i < theta.size(); ++i) {    
           const  double val = std::inner_product(lambdaNew.begin(),
                                                  lambdaNew.end(),
                                                  theta[i].begin(),
                                                  delta[j]);

            bool xnew = rlogis() <= val;
            acc += theta[i][d] * (x[i][j] - xnew) * C[j][d];    
        }
    }
    return acc;
}

I would rather put x-xnew init its own temporary than defining xnew

bool xnew = rlogis() <= val;
acc += theta[i][d] * (x[i][j] - xnew) * C[j][d]; 

becomes:

const int factor = x[i][j] - (rlogis() <= val ? 0 : 1);
acc += factor * theta[i][d] * C[j][d]; 

And we can put it together:

double orig(const mat& theta, const mat& C,
            const vec& delta, vec lambdaNew,
            const int d, std::vector<std::vector<int> >& x) {

    double acc = 0.0;
    for (std::size_t j = 0; j < C.size(); ++j) {
        std::transform(lambdaNew.begin(), 
                       lambdaNew.end(), 
                       C[j].begin(),        
                       lambdaNew.begin(),
                       std::multiplies<double>());

        for (std::size_t i = 0; i < theta.size(); ++i) {    
            const double val = std::inner_product(lambdaNew.begin(),
                                                  lambdaNew.end(),
                                                  theta[i].begin(),
                                                  delta[j]);

            const int factor = x[i][j] - (rlogis() <= val ? 0 : 1);
            acc += factor * theta[i][d] * C[j][d];   
        }
    }
    return acc;
}

I would guess, that there will only be two really improvements. The first is the switch away from vecbool. The second is the direct copy of lambdaNew that ommits the zeroing of tmp and the possibility for move construction.