Thread Subject: Vectorization : Min/Max Function

I am just getting onboard with the notion of vectorizing my for loops, and while I have done quite a bit of reading on it, for some reason I may be too daft to figure out why I can't determine how to vectorize the following for loop:

var(1:length(data)) = 0;
for i = 1:length(data)
     var(i) = max(data(i-30:i-01));
end

Thanks!

"Eric " <elb@windwardinv.com> wrote in message <h33055$eah$1@fred.mathworks.com>...
> I am just getting onboard with the notion of vectorizing my for loops, and while I have done quite a bit of reading on it, for some reason I may be too daft to figure out why I can't determine how to vectorize the following for loop:
>
> var(1:length(data)) = 0;
> for i = 1:length(data)
> var(i) = max(data(i-30:i-01));
> end
>
> Thanks!

As an additional thought I found that I could perform averaging in this context no problem using the filter command...not sure that it applies with the min and max function though.

Before:
var(1:length(data)) = 0;
for i = 1:length(data)
     var(i) = mean(data(i-30:i-01));
end

After:
var = filter(ones(1,30)/3,1,data);

"Eric " <elb@windwardinv.com> wrote in message <h330vd$8rd$1@fred.mathworks.com>...
> "Eric " <elb@windwardinv.com> wrote in message <h33055$eah$1@fred.mathworks.com>...
> > I am just getting onboard with the notion of vectorizing my for loops, and while I have done quite a bit of reading on it, for some reason I may be too daft to figure out why I can't determine how to vectorize the following for loop:
> >
> > var(1:length(data)) = 0;
> > for i = 1:length(data)
> > var(i) = max(data(i-30:i-01));
> > end
> >
> > Thanks!
>
> As an additional thought I found that I could perform averaging in this context no problem using the filter command...not sure that it applies with the min and max function though.
>
> Before:
> var(1:length(data)) = 0;
> for i = 1:length(data)
> var(i) = mean(data(i-30:i-01));
> end
>
> After:
> var = filter(ones(1,30)/3,1,data);

% A few comments:
1. You may want to avoid using var as a variable
since it is the MATLAB "builtin" function for variance
(which you might be interested in with these types of problems).
Also same goes for i (MATLAB's sqrt(-1)).
2. Am I missing something, since when i < 30
you will be using a negative index for data?

"Eric " <elb@windwardinv.com> wrote in message <h33055$eah$1@fred.mathworks.com>...
> I am just getting onboard with the notion of vectorizing my for loops, and while I have done quite a bit of reading on it, for some reason I may be too daft to figure out why I can't determine how to vectorize the following for loop:
>
> var(1:length(data)) = 0;
> for i = 1:length(data)
> var(i) = max(data(i-30:i-01));
> end
>
> Thanks!

% data is row vector
win = 30;
var = max(hankel(data(1:win),[data(win:end) -inf(1,win-1)]),1)

Bruno

Sorry there is a typo, it should read:

% data is row vector
win = 30;
maxdata = max(hankel(data(1:3),[data(3:end) -inf(1,2)]), [], 1)

An explanation would help: maxdata will have the same size as data. The
- 1st element of maxdata is max(data(1:30))
- 2nd element of maxdata is max(data(2:31))
...
- before last of maxdata is max(data(end-1:end))
- last element of maxdata is data(end)

% Bruno

Jamais deux sans trois!
 
> % data is row vector
win = 30;
maxdata = max(hankel(data(1:win),[data(win:end) -inf(1,win-1)]), [], 1)

Bruno

On Jul 8, 2:43 pm, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Jamais deux sans trois!
>
> > % data is row vector
>
> win = 30;
> maxdata = max(hankel(data(1:win),[data(win:end) -inf(1,win-1)]), [], 1)
>
> Bruno

This is really good, but I doubt in such case vectorized code can be
any faster than the for loop.

Siyi Deng <mr.siyi.deng@gmail.com> wrote in message <690dd0d7-2130-4d08-
> This is really good, but I doubt in such case vectorized code can be
> any faster than the for loop.

Especially since the given For loop will error out on it's first pass. (As someone noticed earlier.)

Ideally what we need is a Mex implementation of a MAX/MIN filter using a better algorithm (one of the fancy one that does use the sliding nature of the problem). Is there any expert on this problem can make a comment?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h33tie$9c9$1@fred.mathworks.com>...
> Ideally what we need is a Mex implementation of a MAX/MIN filter using a better algorithm (one of the fancy one that does use the sliding nature of the problem). Is there any expert on this problem can make a comment?
>
> Bruno

I was thinking a similar thing with filter. Then my mind drifted to accumarray. Of course I think it would have to use the @max, and so be slow. But it let me to think this could be vectorized another way. I have only worked out the case for when win = 2 here, just to see if it might be worth pursuing the general case. It seems fairly fast for large data vectors, so it might be worth developing. Of course using only C-code stock functions helps.


subs = ones(length(data)*2-1,1);
subs(3:2:end) = 0; % Can you tell this started as an accumarray arg?
A = max(reshape([data(cumsum(subs)),-inf],2,[]));


I'm sure this still wouldn't be faster than a MEX, but it was interesting. It is late here, signing off till tomorrow....

"Eric " <elb@windwardinv.com> wrote in message <h33055$eah$1@fred.mathworks.com>...
> I am just getting onboard with the notion of vectorizing my for loops, and while I have done quite a bit of reading on it, for some reason I may be too daft to figure out why I can't determine how to vectorize the following for loop:
>
> var(1:length(data)) = 0;
> for i = 1:length(data)
> var(i) = max(data(i-30:i-01));
> end
>
> Thanks!

You might be interested in my SLIDEFUN function (which does use for-loops as well):
http://www.mathworks.com/matlabcentral/fileexchange/12550

For instance, one of the examples is quite close to what you want

% Example 1) Sliding max filter - return the maximum of every three
% consecutive elements:
% V = [1 2 3 9 4 2 1 1 5 6] ;
% R = slidefun(@max, 3, V)
% % -> [2 3 9 9 9 4 2 5 6 6]
% % So R(i) = max(R(i-1:i+1)) ;
% % and R(1) = max(V(1:2))

hth
Jos

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h333t5$k5b$1@fred.mathworks.com>...
> Jamais deux sans trois!
>
> > % data is row vector
> win = 30;
> maxdata = max(hankel(data(1:win),[data(win:end) -inf(1,win-1)]), [], 1)
>
> Bruno


First off...thank you to everyone who wrote in!!!
Bruno, this seems to work perfectly....my assumption is that hankel is totally vector based with no loops?? I'll obviously run the profiler for pre-post comparison....My data vectors are quite long and operated on multiple times, so this should really provide a very significant speed boost!
I changed it such that I ignore the last few max functions where there is not the full window size and through some zero padding in the front....Is this the best way for me to rewrite your code to achieve that?

win=4; maxdata = ([zeros(1,win-1) max(hankel(b(1:win),b(win:end)), [], 1)])

"Eric " <elb@windwardinv.com> wrote in message <h34qvk$r1b$1@fred.mathworks.com>...

> Bruno, this seems to work perfectly....my assumption is that hankel is totally vector based with no loops?? I'll obviously run the profiler for pre-post comparison....My data vectors are quite long and operated on multiple times, so this should really provide a very significant speed boost!
> I changed it such that I ignore the last few max functions where there is not the full window size and through some zero padding in the front....Is this the best way for me to rewrite your code to achieve that?
>
> win=4; maxdata = ([zeros(1,win-1) max(hankel(b(1:win),b(win:end)), [], 1)])

Yes Hankel is vectorized. I have some doubt about it runs faster than for-loop. If it's faster, I'm happy for you.

Your zero padding look all right to me.

Bruno

It seems to me the padding should be at the back, since this is where the maxing isn't occurring. But that is your choice of course. This is a bit faster than the hankel solution, but I still have my doubts about a good For loop not being just as fast or faster.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function benchmovingmax
% Data
B = randperm(90000);
win = 6;

% Engine
tic
subs = ones(win*length(B)-(win-1)*win,1);
subs(win+1:win:end) = 2-win;
maxdata3 = [zeros(1,win-1) max(reshape(B(cumsum(subs)),win,[]))];
toc % Elapsed time is 0.015388 seconds.

% Compare
tic
maxdata2 = ([zeros(1,win-1) max(hankel(B(1:win),B(win:end)), [], 1) ]);
toc % Elapsed time is 0.023616 seconds.

isequal(maxdata2,maxdata3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With win=200 I get:
Elapsed time is 1.246890 seconds.
Elapsed time is 1.930484 seconds.

Thanks Matt. Great solution...the reshape method is certainly a speed improvement on the hankel method...I wonder if there is a way to utilize the filter function instead as this method for doing a similar mean function (as oppose to min/max) is far and away faster (almost 2 orders of magnitude) than both the reshape and hankel approaches.
Either way, my guess is that reshape should suffice compared to the nested for loops I am currently using.
Thanks again.

Eric

"Matt Fig" <spamanon@yahoo.com> wrote in message <h35ee6$def$1@fred.mathworks.com>...
> It seems to me the padding should be at the back, since this is where the maxing isn't occurring. But that is your choice of course. This is a bit faster than the hankel solution, but I still have my doubts about a good For loop not being just as fast or faster.
>
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> function benchmovingmax
> % Data
> B = randperm(90000);
> win = 6;
>
> % Engine
> tic
> subs = ones(win*length(B)-(win-1)*win,1);
> subs(win+1:win:end) = 2-win;
> maxdata3 = [zeros(1,win-1) max(reshape(B(cumsum(subs)),win,[]))];
> toc % Elapsed time is 0.015388 seconds.
>
> % Compare
> tic
> maxdata2 = ([zeros(1,win-1) max(hankel(B(1:win),B(win:end)), [], 1) ]);
> toc % Elapsed time is 0.023616 seconds.
>
> isequal(maxdata2,maxdata3)
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> With win=200 I get:
> Elapsed time is 1.246890 seconds.
> Elapsed time is 1.930484 seconds.

I have coded a simple for loop that can be faster than the vectorized solution (both Matt's and Hankel) for medium window size. According to my benchmark test when applied on "random" vector, it's worth to use from win>10.

function benchrunningmax()

% Data
B = randperm(90000);
win = 30;

% Matt's Engine
tic
subs = ones(win*length(B)-(win-1)*win,1);
subs(win+1:win:end) = 2-win;
maxdata1 = [zeros(1,win-1) max(reshape(B(cumsum(subs)),win,[]))];
toc % Elapsed time is 0.132831 seconds.

tic
maxdata2 = zeros(size(B));
m = max(B(1:win));
for k=win+1:length(B)
    maxdata2(k-1) = m;
    if B(k-win) < m
        m = max(m, B(k));
    else
        m = max(B(1+k-win:k));
    end
end
maxdata2(k) = m;
toc % Elapsed time is 0.029430 seconds.

isequal(maxdata1,maxdata2)

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h382o9$f5d$1@fred.mathworks.com>...
> I have coded a simple for loop that can be faster than the vectorized solution (both Matt's and Hankel) for medium window size. According to my benchmark test when applied on "random" vector, it's worth to use from win>10.


We were right to doubt (see both our comments above)! This kinda ties in nicely with the question Loren asks at the end of this weeks blog, eh?

Hey Bruno,
You don't mind if I piggy back off your hard work do you?


function benchrunningmax()
% Bruno!!

% Data
B = randperm(90000);
win = 2;
t = zeros(1,3);

% Matt's Engine
tic
subs = ones(win*length(B)-(win-1)*win,1);
subs(win+1:win:end) = 2-win;
maxdata1 = [zeros(1,win-1) max(reshape(B(cumsum(subs)),win,[]))];
t(1) = toc;

tic % Bruno's For loop
maxdata2 = zeros(size(B));
m = max(B(1:win));
for k=win+1:length(B)
    maxdata2(k-1) = m;
    if B(k-win) < m
        m = max(m, B(k));
    else
        m = max(B(1+k-win:k));
    end
end
maxdata2(k) = m;
t(2) = toc;


tic % Fig's Modification
maxdata3 = zeros(size(B));
m = max(B(1:win));
for k=win+1:length(B)
    maxdata3(k-1) = m;
    if B(k-win) < m
        m = max(m, B(k));
    else
        m = B(1+k-win);
        for ii = 1+k-win+1:k
            if B(ii)>m
                m = B(ii);
            end
        end
    end
end
maxdata3(k) = m;
t(3) = toc;

isequal(maxdata1,maxdata2,maxdata3)
t./min(t)
% For win = 100: 151.38 4.7816 1
% For win = 2: 1.2156 8.1736 1

"Matt Fig" <spamanon@yahoo.com> wrote in message <h385po$663$1@fred.mathworks.com>...
> Hey Bruno,
> You don't mind if I piggy back off your hard work do you?

Never. NWG is mainly about sharing and mutual stimulation. You just did that.

Bruno

On 10 Juli, 22:03, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> "Matt Fig" <spama...@yahoo.com> wrote in message <h385po$66...@fred.mathworks.com>...
> > Hey Bruno,
> > You don't mind if I piggy back off your hard work do you?
>
> Never. NWG is mainly about sharing and mutual stimulation. You just did that.
>
> Bruno

This subject was discussed in sep/oct 2006. See e.g.
http://groups.google.se/group/comp.soft-sys.matlab/browse_thread/thread/
3c33975d3646e6ea/e2ed5ef4d75a8445?hl=sv#e2ed5ef4d75a8445
or search for "running max" in this NG. See the mex-function rmax.c
for the fastest version.

Best regards

>
> This subject was discussed in sep/oct 2006. See e.g.
> http://groups.google.se/group/comp.soft-sys.matlab/browse_thread/thread/
> 3c33975d3646e6ea/e2ed5ef4d75a8445?hl=sv#e2ed5ef4d75a8445
> or search for "running max" in this NG. See the mex-function rmax.c
> for the fastest version.
>

I see the same ideas of limiting the comparison is pup up twice. What I do not like is: in the worse case, this trick will not save the number of comparison to the brute force sliding/max fore example when applying on the sequence that is monotonically decrease. I'm pretty sure there are better algorithm out there.

When google; I see this article has been often cited

M. Brookes “Algorithms for max and min filters with improved worst-case performance” IEEE Transactions on Circuits and Systems, Volume 47, Issue 9, Sep 2000, pp. 930–935

Does anyone has a copy of it?

Bruno

I found this paper: http://www.archipel.uqam.ca/309/1/webmaximinalgo.pdf

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h39fme$ing$1@fred.mathworks.com>...
> I found this paper: http://www.archipel.uqam.ca/309/1/webmaximinalgo.pdf

I have tried the algorithm in this paper (Matlab coding) it seems working all right. The next step would be MEX it.

Bruno

Bellow is the first attempt of MEX code of the algorithm of the Lemire's paper.

This algorithm can of course be extended on ND problems (min/max while sliding on a hypercube) by applying successively in different dimensions.

It will be part of soon submitted to FEX. I would like to hear any comment. Thanks.

Bruno

/*************************************************************************
 * MATLAB MEX ROUTINE LEMIRE_ENGINE.C
 * [minval maxval] = LEMIRE_ENGINE(a, window)
 *
 * INPUTS
 * A: double real vector
 * window: size of the sliding window, >2
 *
 * OUTPUTS
 * minval, maxval: running min/max, dimension length(A)-window+1
 * minval(1) is min(A(1:win))
 * ...
 * minval(end) is min(A(end-window+1:end))
 * Same indexing scheme for maxval
 *
 * Algorithm: Lemire's "STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN
 * THREE COMPARISONS PER ELEMENT" Nordic Journal of Computing, Volume 13,
 * Number 4, pages 328-339, 2006.
 *
 * Compilation:
 * >> mex -O -v lemire_engine.c % add -largeArrayDims on 64-bit computer
 *
 * Author: Bruno Luong <brunoluong@yahoo.com>
 * History
 * Original: 12/July/2009
 ************************************************************************/

#include "mex.h"
#include "matrix.h"

/* Define the name for Input/Output ARGUMENTS */
#define A prhs[0]
#define WINDOW prhs[1]
#define MINVAL plhs[0]
#define MAXVAL plhs[1]

/* Gateway of LEMIRE_ENGINE */
void mexFunction(int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[]) {
    
    mxClassID ClassID;
    double *a, *minval, *maxval;
    mwSize i, n, window, size, left;
    mwSize *U, *L; /* wedge */
    mwSize nU, nL; /* wedge number of elements */
    mwSize Ufirst, Lfirst, Ulast, Llast;
    
    /* Check number of arguments */
    if (nrhs!=2)
        mexErrMsgTxt("LEMIRE_ENGINE: two arguments are required.");
    
    /* Get class of input matrix */
    ClassID = mxGetClassID(A);
   
   /* TO DO other classes */
    switch (ClassID) {
        case mxDOUBLE_CLASS:
            break;
        default:
            mexErrMsgTxt("LEMIRE_ENGINE: Class not supported.");
    }
            
    a = mxGetPr(A);
    
    /* Get the number of elements of A */
    n = mxGetM(A)*mxGetN(A);
    
    /* Get the window size */
    window = (mwSize)(*mxGetPr(WINDOW));
    if (window<3) /* Check if it's valid */
        mexErrMsgTxt("LEMIRE_ENGINE: windows must be 2 or greater.");
    
    /* Allocate wedges buffers for L and U, each is size (window+1)*/
    size = (window+1);
    L = mxMalloc((2*size)*sizeof(mwSize));
    if (L==NULL)
        mexErrMsgTxt("LEMIRE_ENGINE: out of memory.");
    U = L + size;
    /* Initialize empty wedge */
    nU = nL = 0;
    Lfirst = Ufirst = 0;
    Llast = Ulast = -1;
    
    /* Create output arrays */
    MINVAL = mxCreateNumericMatrix(1, n-window+1, ClassID, mxREAL);
    MAXVAL = mxCreateNumericMatrix(1, n-window+1, ClassID, mxREAL);
    
    minval = mxGetPr(MINVAL);
    maxval = mxGetPr(MAXVAL);
    
    /* Let's go */
    for (i=1; i<n; i++) {
        
        left = i-window;
        /* Assign the output */
        if (left >= 0) {
            maxval[left] = a[(nU? U[Ufirst] : i-1)];
            minval[left] = a[(nL? L[Lfirst] : i-1)];
        }
        
        if (a[i] > a[i-1])
        {
            /* L = pushback(L, i-1); */
            nL++;
            if ((++Llast) == size) Llast = 0;
            L[Llast] = i-1;
            /* if i==window+getfirst(L), L = popfront(L); end */
            if (left == L[Lfirst]) {
                nL--;
                if ((++Lfirst) == size) Lfirst = 0;
            }
            while (nU) {
                /* if a(i)<=a(getlast(U)) */
                if (a[i] <= a[U[Ulast]]) {
                    /* if i==window+getfirst(U), U = popfront(U); end */
                    if (left == U[Ufirst]) {
                        nU--;
                        if ((++Ufirst) == size) Ufirst = 0;
                    }
                    break;
                } /* if */
                /* U = popback(U); */
                nU--;
                if ((--Ulast) < 0) Ulast += size;
            } /* while */
        }
        else {
            /* U = pushback(U, i-1); */
            nU++;
            if ((++Ulast) == size) Ulast = 0;
            U[Ulast] = i-1;
            /* if i==window+getfirst(U), U = popfront(U); end */
            if (left == U[Ufirst]) {
                nU--;
                if ((++Ufirst) == size) Ufirst = 0;
            }
            while (nL) {
                /* if a(i)>=a(getlast(L)) */
                if (a[i] >= a[L[Llast]]) {
                    /* if i==window+getfirst(L), L = popfront(L); end */
                    if (left == L[Lfirst]) {
                        nL--;
                        if ((++Lfirst) == size) Lfirst = 0;
                    }
                    break;
                } /* if */
                /* L = popback(L); */
                nL--;
                if ((--Llast) < 0) Llast += size;
            } /* while */
        } /* if */
                        
    } /* for-loop */
    
    /* Assign the last output */
    left = n-window;
    maxval[left] = a[(nU? U[Ufirst] : n-1)];
    minval[left] = a[(nL? L[Lfirst] : n-1)];
    
    /* Free the buffer */
    mxFree(L);
    
    return;
} /* LEMIRE_ENGINE */