Efficient algorithm to produce closest triplet from 3 arrays?

Question

I need to implement an algorithm in C++ that, when given three arrays of unequal sizes, produces triplets a,b,c (one element contributed by each array) such that max(a,b,c) - min(a,b,c) is minimized. The algorithm should produce a list of these triplets, in order of size of max(a,b,c)-min(a,b,c). The arrays are sorted.

I've implemented the following algorithm (note that I now use arrays of type double), however it runs excruciatingly slow (even when compiled using GCC with -03 optimization, and other combinations of optimizations). The dataset (and, therefore, each array) has potentially tens of millions of elements. Is there a faster/more efficient method? A significant speed increase is necessary to accomplish the required task in a reasonable time frame.

void findClosest(vector<double> vec1, vector<double> vec2, vector<double> vec3){

 //calculate size of each array
    int len1 = vec1.size();
    int len2 = vec2.size();
    int len3 = vec3.size();

    int i = 0; int j = 0; int k = 0; int res_i, res_j, res_k;
    int diff = INT_MAX;

    int iter = 0; int iter_bound = min(min(len1,len2),len3);
    while(iter < iter_bound)
        while(i < len1 && j < len2 && k < len3){

            int minimum = min(min(vec1[i], vec2[j]), vec3[k]);
            int maximum = max(max(vec1[i], vec2[j]), vec3[k]);

            //if new difference less than previous difference, update difference, store
            //resultants
            if(fabs(maximum - minimum) < diff){ diff = maximum-minimum; res_i = i; res_j = j; res_k = k;}

            //increment minimum value
            if(vec1[i] == minimum) ++i;
            else if(vec2[j] == minimum) ++j;
            else ++k;

    }

    //"remove" triplet
    vec1.erase(vec1.begin() + res_i);
    vec2.erase(vec2.begin() + res_j);
    vec3.erase(vec3.begin() + res_k);

    --len1; --len2; --len3;

    ++iter_bound;

}

Answer 1

OK, you're going to need to be clever in a few ways to make this run well.

The first thing that you need is a priority queue, which is usually implemented with a heap. With that, the algorithm in pseudocode is:

Make a priority queue for possible triples in order of max - min, then how close median is to their average.
Make a pass through all 3 arrays, putting reasonable triples for every element into the priority queue
While the priority queue is not empty:
    Pull a triple out
    If all three of the triple are not used:
        Add triple to output
        Mark the triple used
    else:
        If you can construct reasonable triplets for unused elements:
            Add them to the queue

Now for this operation to succeed, you need to efficiently find elements that are currently unused. Doing that at first is easy, just keep an array of bools where you mark off the indexes of the used values. But once a lot have been taken off, your search gets long.

The trick for that is to have a vector of bools for individual elements, a second for whether both in a pair have been used, a third for where all 4 in a quadruple have been used and so on. When you use an element just mark the individual bool, then go up the hierarchy, marking off the next level if the one you're paired with is marked off, else stopping. This additional data structure of size 2n will require an average of marking 2 bools per element used, but allows you to find the next unused index in either direction in at most O(log(n)) steps.

The resulting algorithm will be O(n log(n)).