Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

How to find vector elements between two values, efficiently

Asked by Gonzalo on 31 May 2011

I need to find all elements that fall between 2 values (L,U) in a matrix (A) with 2.8 million rows. I'm currently doing this:

[ind,~] = find(A(:,1) >= L & A(:,1) < U);

Variables U and L are updated 41 times inside a loop. This single line of code is slowing down the program by 20 minutes! Any ideas how can I speed this up?

2 Comments

Angus on 20 Jun 2011

The fastest way I can think of is the following:

data = randn(3000000,3);
inds = not(abs(sign(sign(L - data) + sign(U - data))));

this will give you a matrix of 1s and 0s indicating the indices where the values are between your two bounds (L, U); note not [L, U].

This should be lightening fast, for 3 columns and 3 million rows it computes in a fraction of a second.

chirag on 18 May 2013

just type this

op=A(A>l & A<U);
Gonzalo

Products

5 Answers

Answer by Matt Fig on 31 May 2011
Accepted answer

Do you need the actual indices or the values? If you only need the values, then it would probably be faster to do:

A(A(:,1) >= L & A(:,1) < U)

3 Comments

Gonzalo on 31 May 2011

I actually need the indices...

Matt Fig on 31 May 2011

Then you may be stuck unless there are other efficiencies you can make. One alternative that I can think of to get the indices would be to use a dummy variable. I am not sure if this would be faster or not. Make IDX before hand, if you are looping....

IDX = uint32(1:size(A,1));
ind = IDX(A(:,1) >= L & A(:,1) < U);

Gonzalo on 31 May 2011

This method reduces the time by 13,7%, which is good. Thanks.

Matt Fig
Answer by James Tursa on 1 Jun 2011

You can try a mex approach. The following file does the exact calculation shown above. If you need to modify it for different columns etc let me know. To mex it, simply put the file someplace on your path, make that directory your current directory, then type

mex findrange.c

The mex program does the calculations fast at the expense of memory. It always allocates enough to hold the indexes of the entire column and then just sets the return size to the amount that it found without reallocating & copying. But from the looks of things this may be a relatively minor temporary memory waste.

/* File: findrange.c                                                */
/* IND = findrange(A,L,U) returns the same result as the following: */
/* IND = find(A(:,1)>=L & A(:,1)<U)                                 */
/* Programmer: James Tursa                                          */
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize i, k, m;
    double L, U;
    double *pr, *ind;
      if( nlhs > 1 ) {
          mexErrMsgTxt("Too many outputs.");
      }
      if( nrhs != 3 ) {
          mexErrMsgTxt("Need exactly 3 inputs.");
      }
      if( !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxGetNumberOfDimensions(prhs[0]) > 2 ) {
          mexErrMsgTxt("First argument must be full double 2D matrix.");
      }
      if( !mxIsNumeric(prhs[1]) || mxGetNumberOfElements(prhs[1]) != 1 ) {
          mexErrMsgTxt("2nd argument must be a scalar.");
      }
      if( !mxIsNumeric(prhs[2]) || mxGetNumberOfElements(prhs[2]) != 1 ) {
          mexErrMsgTxt("3rd argument must be a scalar.");
      }
      L = mxGetScalar(prhs[1]);
      U = mxGetScalar(prhs[2]);
      pr = mxGetPr(prhs[0]);
      m = mxGetM(prhs[0]);
      plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL);
      ind = mxGetPr(plhs[0]);
      for( i=0; i<m; i++ ) {
          if( *pr >= L && *pr < U ) {
              *ind++ = i+1;
          }
          pr++;
      }
      mxSetM(plhs[0],ind-mxGetPr(plhs[0]));
  }

2 Comments

Jan Simon on 1 Jun 2011

Instead of the DOUBLE output, an UINT32 vector would use the half memory:
uint32_T *ind, i, m;
plhs[0] = mxCreateNumericMatrix(m,1,mxUINT32_CLASS, mxREAL);
ind = (uint32_T *) mxGetData(plhs[0]);

With "for (i = 1; i <= m; i++), ... *ind++=i;" you can save some milliseconds.

James Tursa on 2 Jun 2011

DOUBLE vs UINT32 memory savings:
True. Whether this is an actual savings overall depends on what is done with the result downstream (which OP doesn't show). If operations are done with it that turn it into a double then the memory savings will turn into a memory waste instead.

Changing the loop index to start from 1 instead of 0:
Yeah, I saw that right after I posted it (I habitually start my loops from 0 in C), but then I thought ... "Why don't I just let Jan find it ..."

James Tursa
Answer by Walter Roberson on 31 May 2011

How "dense" are the values? It might be faster to use

[ind1, ~] = find(A(:,1)>=L);
ind = ind1(A(ind1,:)<U);

Or reversing the order of the test, putting the test less likely to succeed first.

This could be cost-effective if relatively few matches are found, reducing the number of ind1 subscripts that need to be looked up in the second step.

2 Comments

Gonzalo on 31 May 2011

It's very dense... I'll give this a try..

Gonzalo on 31 May 2011

This approach doesn't reduce the time.. but thanks anyway

Walter Roberson
Answer by Angus on 22 Jun 2011

The fastest way I can think of is the following:

data = randn(3000000,3);
inds = not(abs(sign(sign(L - data) + sign(U - data))));

this will give you a matrix of 1s and 0s indicating the indices where the values are between your two bounds (L, U); note not [L, U].

This should be lightening fast, for 3 columns and 3 million rows it computes in a fraction of a second.

2 Comments

Matt Fig on 22 Jun 2011

But surely this is not faster than:

idx = (A >= L & A < U);

You are calling SIGN 3 times, MINUS twice, PLUS once, ABS once - all on double data, then a logical negation. All that is involved above expression is two logical comparisons and a logical AND.

Angus on 22 Jun 2011

Matt. You're right it shouldn't be and isn't. It appeared and I believe I misread that this solution only lowered the execution time by 13.7% on a 20 minute execution time for merely the find statement. The example I gave above with 3 million rows and on 3 columns runs in sub 1 sec. Your answer definitely makes more sense.

Angus
Answer by chirag on 18 May 2013

just type this

op=A(A>l & A<U);

0 Comments

chirag

Contact us