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

Thread Subject:
fast way to get second output of min

Subject: fast way to get second output of min

From: Christian

Date: 6 Mar, 2014 22:35:08

Message: 1 of 6

Hi all,

I'm looking for a faster way for the following code:

    for k1=1:K1
        for k21=1:K2
            for k22=1:K2
                [~, xind(k1,k21,k22)] = min(abs(yp(k1,k21,k22)-x)); %nearest neighbor of yp in x
            end
        end
    end

where x is a vector of dimension 1:K1. I have to do this operation zillions of times, so it might be nice to have a speedup. Converting it to a mex file helps a bit, but still it's too slow for what I need. Maybe one could leverage on the fact that I only need the second output element of min.

Thanks a lot,
Christian

Subject: fast way to get second output of min

From: dpb

Date: 7 Mar, 2014 02:12:43

Message: 2 of 6

On 3/6/2014 4:35 PM, Christian wrote:
> Hi all,
>
> I'm looking for a faster way for the following code:
>
> for k1=1:K1
> for k21=1:K2
> for k22=1:K2
> [~, xind(k1,k21,k22)] ...

1) Have you preallocated xind?

2) Reverse the order of the loops and traverse in column-major order.
Your loop order could be causing cache misses continually depending on
size of the array. This could be a reason the mex didn't help any
significant amount.

I don't see anything to be done just because you're only looking at the
location.

--

Subject: fast way to get second output of min

From: Christian

Date: 7 Mar, 2014 04:52:11

Message: 3 of 6

Hi,

thanks for your reply. Yes, I have preallocated xind as xind=NaN(K1,K2,K3). For column-major order I would put the loop as:

for k21=1:K2
    for k1=1:K1
         for k22=1:K2
? That is, I start with the column, then row, and then all the other dimensions of the matrix?

Anyway, I have tried it in different ways and the results were all very similar. When my dimensions of K1 and K2 are 500 and 7, respectively, the normal matlab file needs 43 seconds for 440 iterations, while the mex file actually needs 52 seconds. For smaller dimensions, like 100 and 3, the mex file is twice as fast as the normal matlab file. Does anybody know why mex is slow?


dpb <none@non.net> wrote in message <lfba03$49a$1@speranza.aioe.org>...
> On 3/6/2014 4:35 PM, Christian wrote:
> > Hi all,
> >
> > I'm looking for a faster way for the following code:
> >
> > for k1=1:K1
> > for k21=1:K2
> > for k22=1:K2
> > [~, xind(k1,k21,k22)] ...
>
> 1) Have you preallocated xind?
>
> 2) Reverse the order of the loops and traverse in column-major order.
> Your loop order could be causing cache misses continually depending on
> size of the array. This could be a reason the mex didn't help any
> significant amount.
>
> I don't see anything to be done just because you're only looking at the
> location.
>
> --

Subject: fast way to get second output of min

From: Roger Stafford

Date: 7 Mar, 2014 07:49:11

Message: 4 of 6

"Christian " <proechri@umich.edu> wrote in message <lfat6s$te$1@newscl01ah.mathworks.com>...
> I'm looking for a faster way for the following code:
>
> for k1=1:K1
> for k21=1:K2
> for k22=1:K2
> [~, xind(k1,k21,k22)] = min(abs(yp(k1,k21,k22)-x)); %nearest neighbor of yp in x
> end
> end
> end
>
> where x is a vector of dimension 1:K1.
- - - - - - - - -
  It seems to me inefficient to have the 'min' function repeatedly scan over all of the x elements inside the innermost for-loop for every element in yp. If you first sort x and then use 'histc' appropriately to find the nearest value of x, it ought to go much faster. Seeking a nearest element in a sorted list of n elements can be done in O(log(n)) as opposed to O(n), which I believe is the way 'histc' works.

  Try the following. The p vector below serves to transform the indices relative to the sorted x2 back to those relative to the original x. I assume here that x is a column vector. If it is a row vector, make the obvious change in the second argument in 'histc'. I also assume that all the elements in x and yp are finite.

 [x2,p] = sort(x);
 [~,xind] = histc(yp,[-inf;(x2(1:K1-1)+x2(2:K1))/2;inf]); % Use x2 midpoints
 xind = reshape(p(xind),size(yp));

Roger Stafford

Subject: fast way to get second output of min

From: Yair Altman

Date: 7 Mar, 2014 10:13:07

Message: 5 of 6

"Roger Stafford" wrote in message <lfbtln$df6$1@newscl01ah.mathworks.com>...
> "Christian " <proechri@umich.edu> wrote in message <lfat6s$te$1@newscl01ah.mathworks.com>...
> > I'm looking for a faster way for the following code:
> >
> > for k1=1:K1
> > for k21=1:K2
> > for k22=1:K2
> > [~, xind(k1,k21,k22)] = min(abs(yp(k1,k21,k22)-x)); %nearest neighbor of yp in x
> > end
> > end
> > end
> >
> > where x is a vector of dimension 1:K1.
> - - - - - - - - -
> It seems to me inefficient to have the 'min' function repeatedly scan over all of the x elements inside the innermost for-loop for every element in yp. If you first sort x and then use 'histc' appropriately to find the nearest value of x, it ought to go much faster. Seeking a nearest element in a sorted list of n elements can be done in O(log(n)) as opposed to O(n), which I believe is the way 'histc' works.
>
> Try the following. The p vector below serves to transform the indices relative to the sorted x2 back to those relative to the original x. I assume here that x is a column vector. If it is a row vector, make the obvious change in the second argument in 'histc'. I also assume that all the elements in x and yp are finite.
>
> [x2,p] = sort(x);
> [~,xind] = histc(yp,[-inf;(x2(1:K1-1)+x2(2:K1))/2;inf]); % Use x2 midpoints
> xind = reshape(p(xind),size(yp));
>
> Roger Stafford


Minor correction to Roger's excellent answer:

   [x2,p] = sort(x);
   [~,xind] = histc(yp,[-inf;(x2(1:end-1)+x2(2:end))/2;inf]); % Use x2 midpoints
   xind = reshape(p(xind),size(yp));

i.e., use end (not K1), for the x2 midpoints in the call to histc

Yair Altman
http://UndocumentedMatlab.com
 

Subject: fast way to get second output of min

From: Christian

Date: 7 Mar, 2014 15:07:08

Message: 6 of 6

That's wonderful!
My x vector already is sorted and strictly monotonically increasing so I could skip line 1 and 3. The normal matlab file now requires 5.5 seconds compared to 43 before, and the mex file is actually down to 3.5 seconds. There are some minor other operations in the code that take up the chunk of the calculation time. So I might post another question on that topic.
For histc I know that it's not available on gpu. But I think I've already asked for a manual workaround for histc, so I might want to revive that thread.

"Yair Altman" wrote in message <lfc63j$2mf$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <lfbtln$df6$1@newscl01ah.mathworks.com>...
> > "Christian " <proechri@umich.edu> wrote in message <lfat6s$te$1@newscl01ah.mathworks.com>...
> > > I'm looking for a faster way for the following code:
> > >
> > > for k1=1:K1
> > > for k21=1:K2
> > > for k22=1:K2
> > > [~, xind(k1,k21,k22)] = min(abs(yp(k1,k21,k22)-x)); %nearest neighbor of yp in x
> > > end
> > > end
> > > end
> > >
> > > where x is a vector of dimension 1:K1.
> > - - - - - - - - -
> > It seems to me inefficient to have the 'min' function repeatedly scan over all of the x elements inside the innermost for-loop for every element in yp. If you first sort x and then use 'histc' appropriately to find the nearest value of x, it ought to go much faster. Seeking a nearest element in a sorted list of n elements can be done in O(log(n)) as opposed to O(n), which I believe is the way 'histc' works.
> >
> > Try the following. The p vector below serves to transform the indices relative to the sorted x2 back to those relative to the original x. I assume here that x is a column vector. If it is a row vector, make the obvious change in the second argument in 'histc'. I also assume that all the elements in x and yp are finite.
> >
> > [x2,p] = sort(x);
> > [~,xind] = histc(yp,[-inf;(x2(1:K1-1)+x2(2:K1))/2;inf]); % Use x2 midpoints
> > xind = reshape(p(xind),size(yp));
> >
> > Roger Stafford
>
>
> Minor correction to Roger's excellent answer:
>
> [x2,p] = sort(x);
> [~,xind] = histc(yp,[-inf;(x2(1:end-1)+x2(2:end))/2;inf]); % Use x2 midpoints
> xind = reshape(p(xind),size(yp));
>
> i.e., use end (not K1), for the x2 midpoints in the call to histc
>
> Yair Altman
> http://UndocumentedMatlab.com
>

Tags for this Thread

No tags are associated with this thread.

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.

Contact us