cumulative maximum - vectorised
gg
I would like to calculate the cumulative maximum of a vector. This should work exactly like cumsum but taking max instead of sum. (E.g. cummax([2 1 3 5 4])is [2 2 3 5 5]) Of course the challenge is to do this in a vectorised way.
My first attempt was
%-------------------------------------------------------------
function cmax = cummax(A)
%find cumulative maximum of the rows of matrix A
%cummax is to max what cumsum is to sum
A=A'; % transpose the input
[n m] = size(A);
onetri = shiftdim(triu(ones(m,m)),-1);
onevec = repmat(onetri,[n,1,1]);
tmp = repmat(A,[1,1,m]);
tmp(~onevec) = -inf;
% transpose output, because input has been transposed.
cmax= squeeze(max(tmp,[],2))';
%-------------------------------------------------------------
The problem with above function is the use of repmat, which uses up too much memory for larger matrices.
Has anyone a clever idea?
Thank you
gg
Bruno Luong
> Dear All,
>
> I would like to calculate the cumulative maximum of a vector. This should work exactly like cumsum but taking max instead of sum. (E.g. cummax([2 1 3 5 4])is [2 2 3 5 5]) Of course the challenge is to do this in a vectorised way.
>
Why it has to be vectorized?
I'm curious to see any vectorized algorithm one can beat this one:
function b=cummaxfor(a)
b=zeros(size(a));
b(1)=a(1);
for k=2:length(a)
b(k)=max(a(k),b(k-1));
end
To gain speed, make a Mex.
Bruno
John D'Errico
This is something I've wanted before. I looked
at the idea of a descending sort, then using
those tags to create the cumulative max. Its
too unwieldy to be more efficient as Matlab
code. And since a sort does more work than
this requires, why bother? Finally, while a mex
would surely gain some speed, an accelerated,
lean matlab loop is probably the best.
John
John D'Errico
I doubt this is worth posting on the
FEX, but it does work on any specified
dimension.
function res = cummax(V,dim)
% cummax: cumulative maximum
% usage: res = cummax(V)
% usage: res = cummax(V,dim)
%
% Arguments: (input)
% V - any numeric vector or array
%
% dim - (OPTIONAL) - scalar positive integer
% Indicates the dimension along which the
% cumulative mamimum is to be taken.
%
% If dim is not supplied, or is left empty,
% then the first non-singleton dimension is
% chosen as the default.
%
% Arguments: (output)
% res - cumulative maximum
%
%
% Example:
% V = rand(2,5)
% V =
% 0.90588 0.5073 0.23138 0.56156 0.090554
% 0.087153 0.35399 0.77111 0.29053 0.52427
% cummax(V,2)
% ans =
% 0.90588 0.90588 0.90588 0.90588 0.90588
% 0.087153 0.35399 0.77111 0.77111 0.77111
%
%
% See also: min, max, cumsum, cumprod
%
%
% Author: John D'Errico
% e-mail: wood...@rochester.rr.com
% Release: 1.0
% Release date: 10/04/08
% empty begets empty, scalars also remain unchanged
if (numel(V) <= 1)
res = V;
return
end
% default for the dimension: first non-singleton
% There must be at least one such dimension, since
% we have already kicked out the scalars. Shiftdim
% does much of the work.
sizeV = size(V);
if (nargin<2) || isempty(dim)
dim = min(find(sizeV~=1));
elseif (numel(dim)>1) || (dim<0)
error('dim must be scalar, positive, numeric, integer')
end
[Vs,nshift] = shiftdim(V,dim-1);
% we may now presume to be always working on the
% first dimension of Vs.
res = zeros(size(Vs));
ns = size(Vs);
n1 = ns(1);
n2 = prod(ns(2:end));
% simple loop is the best
res(1,:) = Vs(1,:);
for i = 2:n1
res(i,:) = max(res(i-1,:),Vs(i,:));
end
% All done, so unshift the dimensions
res = shiftdim(res,length(sizeV) - nshift);
Jos
Take a look at SLIDEFUN
a = [2 1 3 5 4]
slidefun(@max,numel(a),a,'backward')
% ans = 2 2 3 5 5
SLIDEFUN can be found here:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=12550&objectType=FILE
Although SLIDEFUN is not vectorized, it is quite efficient.
Possible oneliners for this particular problem include:
A = [2 1 3 5 4]
fliplr(max(hankel(fliplr(A))))
max(triu(toeplitz(A)))
which only works for positive values in A. If you have the stats toolbox, the following can be used to circumvent that problem:
A = [-2 3 1 5 4]
fliplr(nanmax(hankel(fliplr(A),[A(1) nan(1,numel(A)-1)])))
hth
Jos
Bruno Luong
>
> Possible oneliners for this particular problem include:
>
> A = [2 1 3 5 4]
> fliplr(max(hankel(fliplr(A))))
> max(triu(toeplitz(A)))
>
> which only works for positive values in A. If you have the stats toolbox, the following can be used to circumvent that problem:
>
> A = [-2 3 1 5 4]
> fliplr(nanmax(hankel(fliplr(A),[A(1) nan(1,numel(A)-1)])))
>
Hah hah very funny Jos, but why try to make us by statistic toolbox:
A = [2 1 3 5 4]
fliplr(max(toeplitz(A,-inf(size(A)))))
Bruno
Bruno Luong
Typo:
> Hah hah very funny Jos, but why try to make us *buying* statistic toolbox
Bruno
Matt Fig
>
>
> I'm curious to see any vectorized algorithm one can beat this one:
>
> function b=cummaxfor(a)
>
> b=zeros(size(a));
> b(1)=a(1);
> for k=2:length(a)
> b(k)=max(a(k),b(k-1));
> end
>
> To gain speed, make a Mex.
>
> Bruno
I am with Bruno that for very large vectors, a loop is the way to go here. However, I would suggest a slight change to his algorithm. This results in an average increase in speed for large vectors.
function b = cummaxfor2(a)
b = zeros(size(a));
[val,idx] = max(a);
b(idx:end) = val;
b(1) = a(1);
for k = 2:idx-1
b(k) = max(a(k),b(k-1));
end
I tested the relative speeds like this:
niter = 300;
t = zeros(niter,2);
for ii = 1:niter
A = round(rand(1,20000)*2000);
f = @() cummaxfor(A);
t(ii,1) = timeit(f);
f = @() cummaxfor2(A);
t(ii,2) = timeit(f);
end
>> sum(t)/min(sum(t))
ans =
1.2876 1.0000
Jos
> Hah hah very funny Jos, but why try to make us by statistic toolbox:
The humor was not intended, as it was not my goal to promote the stats toolbox. On the contrary, I think that functions like nanmax, nanmean etc should be part of basic matlab.
Jos
Jan Simon
So I tried it with a temporary variable and got a speedup of 50% (Matlab 6.5, vector length: 1000).
> "Bruno Luong" <b.l...@fogale.findmycountry> wrote in
> > function b=cummaxfor(a)
> >
> > b=zeros(size(a));
> > b(1)=a(1);
> > for k=2:length(a)
> > b(k)=max(a(k),b(k-1));
> > end
function b = cummaxfor3(a)
b = zeros(size(a));
t = a(1);
for k = 1:length(a)
if a(k) > t
t = a(k);
end
b(k) = t;
end
Nevertheless, accessing two vectors [a] and [b] occupies the processor cache twice. So further 10% can be saved:
function a = cummaxfor4(a)
t = a(1);
for k = 1:length(a)
if a(k) > t
t = a(k);
end
a(k) = t;
end
Good luck, Jan
Bruno Luong
function a = cummaxfor5(a)
t = a(1);
for k = 2:length(a)
if a(k) > t
t = a(k);
else
a(k) = t;
end
end
% Bruno
Jan Simon
> function a = cummaxfor5(a)
> t = a(1);
> for k = 2:length(a)
> if a(k) > t
> t = a(k);
> else
> a(k) = t;
> end
> end
function a = cummaxfor6(a)
t = a(1);
for k = 2:length(a)
if a(k) < t
a(k) = t;
else
t = a(k);
end
end
But if racing is really needed (is it ever?), surely a MEX script is preferable (Matlab 6.5, BCC 5.5, vector length: 1000 ==> Mex is 65% faster than cummaxfor6, speedup grows with vector length):
------------ 8< ----------------
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *P, *Pend, M;
plhs[0] = mxDuplicateArray(prhs[0]);
if (mxIsEmpty(plhs[0])) {
return;
}
P = mxGetPr(plhs[0]);
Pend = P + mxGetNumberOfElements(plhs[0]);
M = *P++;
while (P < Pend) {
if (*P < M) {
*P++ = M;
} else {
M = *P++;
}
}
return;
}
-------------------------- >8 ------------
Jan
Bruno Luong
> For not strictly ascending values, the number of calls to the ELSE branch of the IF statement is higher than to the IF branch. Then the number of jumps to ELSE can be reduced by swapping the branchs (further 3% acceleration, but depending on the input data!):
But there must be a JUMP machine instruction is at the end of the first branch of IF, so I'm not sure this swapping can reduce the CPU time. In fact I ran six times of 1000 random data (on 10^3-10^5 vector length) in length, and the cummaxfor6 is slighly slower (about 0-2%) than cummaxfor5. Put a "continue" statement in the branch makes time even worse!
Bruno
Matt
> Take a look at SLIDEFUN
>
> a = [2 1 3 5 4]
> slidefun(@max,numel(a),a,'backward')
> % ans = 2 2 3 5 5
>
> SLIDEFUN can be found here:
> http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=12550&objectType=FILE
>
> Although SLIDEFUN is not vectorized, it is quite efficient.
It would be good to have slidefun() operate columnwise on matrix inputs as well as just vectors. That way, you could also make 2D sliding window functions as in
A=slidefun(@max,3, A ,'center');
A=slidefun(@max,3, A.' ,'center');
Of course, it would be great to have it as a mex also, but I guess that might be hard...
Jos
Nice idea, Matt. I will think about how to implement it into SLIDEFUN.
Jos
Bryan
It works if the vector is strictly positive.
a = [2 1 3 5 4]';
n = length(a);
A = kron(a,ones(1,n)); % creates a matrix whose columns are repeats of a.
b = max(triu(A))'; % b is the cumulative max
"gg" wrote in message <gc7isd$4a3$1...@fred.mathworks.com>...