Thread Subject: vectorizing challenge

Hi all,

I'm having a hard time trying to vectorize the following code. Suppose I have a vector of interest V, and a reference vector R which is monotonically increasing (it's actually a cumulative distribution function). Both vectors contain only real numbers, except that I've modified R by tacking on -Inf at the beginning and +Inf at the end (the reason for which is hopefully obvious from the code below).

For each element of V, I want to find the first element of R that is greater than or equal to it.

To do this with a FOR loop is trivial (and yes, I realize the two lines in the loop can be combined into one; writing it this way just for clarity):

result = zeros(1,length(V));
for kk = 1:length(V);
  ind = find( V(kk) <= R, 1, 'first' );
  result(kk) = R(ind);
end;

But when V is very large, this is slow. Can anyone think of a clean way to 'vectorize' the FIND command?

Thanks very much in advance!

HJ

In article <guvhl1$qif$1@fred.mathworks.com>,
 "Heywood " <heywoodj123@yahoo.com> wrote:

> Hi all,
>
> I'm having a hard time trying to vectorize the following code. Suppose I have
> a vector of interest V, and a reference vector R which is monotonically
> increasing (it's actually a cumulative distribution function). Both vectors
> contain only real numbers, except that I've modified R by tacking on -Inf at
> the beginning and +Inf at the end (the reason for which is hopefully obvious
> from the code below).
>
> For each element of V, I want to find the first element of R that is greater
> than or equal to it.
>
> To do this with a FOR loop is trivial (and yes, I realize the two lines in
> the loop can be combined into one; writing it this way just for clarity):
>
> result = zeros(1,length(V));
> for kk = 1:length(V);
> ind = find( V(kk) <= R, 1, 'first' );
> result(kk) = R(ind);
> end;
>
> But when V is very large, this is slow. Can anyone think of a clean way to
> 'vectorize' the FIND command?
>
> Thanks very much in advance!
>
> HJ

How about this?

  % Create some sample data.
  dist = rand(1,10);
  R = [0,cumsum(dist)];
  R = R/R(end);

  % Find result.
  n = length(R);
  V = rand(1,100);
  k = interp1(R,0:n,V);
  result = R(ceil(k) + 1);

You are still finding the appropriate interval, but it's done
efficiently inside interp1.

--
Doug Schwarz
dmschwarz&ieee,org
Make obvious changes to get real email address.

"Heywood " <heywoodj123@yahoo.com> wrote in message <guvhl1$qif$1@fred.mathworks.com>...
> Hi all,
>
> I'm having a hard time trying to vectorize the following code. Suppose I have a vector of interest V, and a reference vector R which is monotonically increasing (it's actually a cumulative distribution function). Both vectors contain only real numbers, except that I've modified R by tacking on -Inf at the beginning and +Inf at the end (the reason for which is hopefully obvious from the code below).
>
> For each element of V, I want to find the first element of R that is greater than or equal to it.
------

This one is probably not as efficient as Doug's, but is nice and compact

[dummy,inds] = max( bsxfun(@ge,R(:),V(:).') );
Result = R(inds);

"Heywood " <heywoodj123@yahoo.com> wrote in message <guvhl1$qif$1@fred.mathworks.com>...

> result = zeros(1,length(V));
> for kk = 1:length(V);
> ind = find( V(kk) <= R, 1, 'first' );
> result(kk) = R(ind);
> end;

If speed is your main concern, don't use INTERP1 (slow), use HISTC or FIND_IDX on FEX

http://www.mathworks.com/matlabcentral/fileexchange/23049

Bruno

Doug Schwarz <see@sig.for.address.edu> wrote in message <see-FBBAFE.21132019052009@news.frontiernet.net>...
> In article <guvhl1$qif$1@fred.mathworks.com>,
> "Heywood " <heywoodj123@yahoo.com> wrote:
>
> > Hi all,
> >
> > I'm having a hard time trying to vectorize the following code. Suppose I have
> > a vector of interest V, and a reference vector R which is monotonically
> > increasing (it's actually a cumulative distribution function). Both vectors
> > contain only real numbers, except that I've modified R by tacking on -Inf at
> > the beginning and +Inf at the end (the reason for which is hopefully obvious
> > from the code below).
> >
> > For each element of V, I want to find the first element of R that is greater
> > than or equal to it.
> >
> > To do this with a FOR loop is trivial (and yes, I realize the two lines in
> > the loop can be combined into one; writing it this way just for clarity):
> >
> > result = zeros(1,length(V));
> > for kk = 1:length(V);
> > ind = find( V(kk) <= R, 1, 'first' );
> > result(kk) = R(ind);
> > end;
> >
> > But when V is very large, this is slow. Can anyone think of a clean way to
> > 'vectorize' the FIND command?
> >
> > Thanks very much in advance!
> >
> > HJ
>
> How about this?
>
> % Create some sample data.
> dist = rand(1,10);
> R = [0,cumsum(dist)];
> R = R/R(end);
>
> % Find result.
> n = length(R);
> V = rand(1,100);
> k = interp1(R,0:n,V);
> result = R(ceil(k) + 1);
>
> You are still finding the appropriate interval, but it's done
> efficiently inside interp1.
 
or

[dummy, k2] = histc(V,R)
result2 = R(k2+1) ;

Jos

> [dummy, k2] = histc(V,R)
> result2 = R(k2+1) ;
>
> Jos

Hi Doug, Matt, Bruno, and Jos,

Thanks for all the suggestions! I learned a lot from all your ideas (I'd never seen bsxfun before, and only vaguely remembered interp1). Ultimately Jos's suggestion was the simplest and did the trick -- resulting in a ~20X speedup in execution time!

Thanks again everyone.

Cheers,

Heywood