Thread Subject: vectorizing challenge

Subject: vectorizing challenge

From: Heywood

Date: 20 May, 2009 00:08:01

Message: 1 of 6

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

Subject: vectorizing challenge

From: Doug Schwarz

Date: 20 May, 2009 01:13:20

Message: 2 of 6

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.

Subject: vectorizing challenge

From: Matt

Date: 20 May, 2009 01:30:20

Message: 3 of 6

"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);

Subject: vectorizing challenge

From: Bruno Luong

Date: 20 May, 2009 07:04:02

Message: 4 of 6

"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

Subject: vectorizing challenge

From: Jos

Date: 20 May, 2009 09:17:16

Message: 5 of 6

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

Subject: vectorizing challenge

From: Heywood

Date: 20 May, 2009 13:08:01

Message: 6 of 6

> [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

Tags for this Thread

Everyone's Tags:

Add a New Tag:

Separated by commas
Ex.: root locus, bode

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Tag Activity for This Thread
Tag Applied By Date/Time
vectorize Heywood 19 May, 2009 20:09:04
rssFeed for this Thread

Contact us at files@mathworks.com