Looking for a faster way of finding the first element larger than a given number in a sorted array

Hi,

I have a vector

a = [1 1 1.01 1 1.05 1 .... 1]

which is basically all numbers >=1. I also have a vector b that is the cummulative sum of this vector (and thus sorted), so each increment in this vector is at least 1 and sometimes (often) more than 1.

I also have a random number r between 0 and the total sum of a. Now, I need to find the first element k in b that is larger than r.

I am currently using a simple find comment:

k = find(b>r, 1)

However, I feel like it should be possible for this to be faster. Since I know that k should be at least floor(r) as increments in the cumsum vector are at least 1, I have tried

n = floor(r)
k = find(b(n:end)>r, 1) + n-1

But this seems to be slower despite skipping searching the first n enties of b. Is there an build-in option in matlab for starting the search at a certain index? If not, is there any other faster way? This function takes about 65% of my total computing time. a Is a vector of size 10.000 but it changes very often and I have to make this call millions of times.

[dpb]

The alternative requires creating the temporary vector b(n:end) that is undoubtedly the overhead problem besides the additional calculation of floor() and the extra arithmetic.

Nothing earth-shaking comes to mind -- perhaps one could test for whether the value of r is past the halfway point and use the reverse test from back to front for those cases, but again might be more overhead introduced in the test than gained in the result.

[Fangjun Jiang]

not to create b, do a while-loop subtracting r with values in a, starting from the first element.

with the resulting k could be first or last, on average, this might be faster.

[Tom van den Bosch]

Thank you! But creating b is not the problem here, as I am able to update that in another way that does not cost much time. The slowest step is purely finding the first element of b larger than k. I'd imagine doing a while loop one-by-one is slower than the single find command but I'll try it.

I'll also try the reverse find for if the value is past halfway, hopefully that helps a bit.

[Fangjun Jiang]

sum(b<r) might give you the answer. index=b<r is true or false. You know the size of b, so don't use find(), you can simply count the zeros or ones in "index" to derive "k".

[James Tursa]

A mex routine with a binary search using truncated starting points will probably be fastest. Do you have C/C++ compiler installed for this?

[Tom van den Bosch]

Thanks for the suggestions all!

I tried sum(b<r) but this was about 30% slower.

I am unsure how the matlab C/C++ compiler works, although I do not doubt that the fastest possible method is through a C++. However, the accepted answer by dpb is sufficient for now as this step is no longer a bottleneck and only takes about 7% of total runtime.

Answer 1

[dpb]

With the JIT compiler, loops are pretty quick any more...in fact, with the starting point as the beginning loop index, that may, indeed, be the fastest option. I'd use a counted loop and break though, not the while I think, although you can compare...

for k=floor(r):numel(b)
  if b(k)>r, break, end
end

[Bruno Luong]

dpb method is very fast

n = 10000;
a=ones(1,n)+0.2*rand(1,n);
b=cumsum(a);
ntests=100000;
rtab = rand(1,ntests)*b(end);
rtab = max(rtab, b(1));

tic
for i=1:ntests
    r = rtab(i);
    k = find(b>r, 1);
end
toc 

Elapsed time is 1.188909 seconds.

tic
for i=1:ntests
    r = rtab(i);
    k = discretize(r,b);
    if b(k) == r
        k = k+1;
    end
end
toc % 1.782346 seconds.

Elapsed time is 1.852914 seconds.

tic
for i=1:ntests
    r = rtab(i);
    for k=floor(r):numel(b)
        if b(k)>r
            break
        end
    end
end
toc

Elapsed time is 0.008726 seconds.

[Tom van den Bosch]

This is amazing, thank you! By far the fastest method of them all and it completely removes this step in my code as a bottleneck.

[dpb]

I had a machine failure so have been offline since -- glad to hear that seems to help; I had one additional idea that might be better than the straight loop which would be a binary search starting at the known minimum location. I haven't tried coding it to test; I just got MATLAB reloaded on the new machine and have a lot more to go to get it rebuilt but I just came to check out briefly before going on...

Answer 2

[Bruno Luong]

You might try

k = discretize(r,b);
if b(k) == r
    k = k+1;
end

[Tom van den Bosch]

Thanks! I tried but it seems to be about 3-4 times as slow.

Answer 3

[Zhunping Xue]

If it's an ordered cumsum vector and about 10,000 entries in size, you probably can simply implement the binary search algorith (you can find it anywhere, say, https://en.wikipedia.org/wiki/Binary_search_algorithm). Matlab's find option is a linear search and does not take advantage of the fact that your vector is ordered. It should come out to be orders faster than what you're doing!

James Tursa had the best solution I can think of, a binary search with mex -- but you don't have to do it in mex.

My really, really primitive code (I'm a biologist, don't laugh):

function m=find_ordered_array(array, value)
% finds value in ordered array, else returns the closest value larger than
% the value sought for


% if(value<array(1) || value > array(end))
%     disp('value outside array');
%     m=-1;
%     return;
% else
% l=1;
tol=1e-15;
h=length(array);
if(array(1)>=value)
    m=1;
    return;
elseif(array(h)<=value)
    m=h;
    return;
else
    l=1;
    h=length(array);
    m=floor((l+h)/2);
    while(l<h-1 && abs(array(m)-value)>tol)

        if(array(m)>value)
            h=m;
        else
            l=m;
        end
        m=floor((l+h)/2);
    end
    
    if(abs(array(m)-value)<tol)
    elseif(abs(array(h)-value)<tol)
        m=h;
    elseif(abs(array(l)-value)<tol)
        m=l;
    else
        m=l;
            
    end
    m=m+1;
end
% end

[Bruno Luong]

@Zhunping Xue discretize in my solution is binary search