Code optimization (3 line function)

I need a function that efficiently applies a reflective bound to a signal: when the signal reaches the bound, it stays there until the signal increases again (example below).

I wrote a short function that applies such a bound to all signals in a matrix (code below).

My question is: does someone know a way to speed up this code? I'm applying it to matrices of size 5000x5000, which takes about a second per matrix. Since i have to do this on thousands of matrices, it would be nice to speed up the function. Thanks!

function Y = reflect(Y,b)
% apply reflecting bound to signals in matrix Y; columns are signals, rows time steps
mhist = zeros(1,size(Y,2));        % history of shifts applied in prev time steps
for ii=1:size(Y,1)                 % loop over time steps
    m = max(b-Y(ii,:)-mhist,0);       % compute required shift for current time step 
    Y(ii,:) = Y(ii,:) + mhist + m;    % apply shift
    mhist = mhist + m;                % add shift to history
end

Here is an example of its output (red=original, black=after applying reflective bound at Y=-50):

The plot was produced using the following code:

rng(1); 
X=0:.1:100; 
Y=cumsum(normrnd(0,5,size(X)))-X; 
plot(X,Y,'r'); 
hold on; 
plot(X,reflect(Y',-50),'k-'); 
plot([X(1) X(end)],[-50 -50],'b');

[luc]

Pre-allocating matrices with nan vallues.

m, Y and mhist change within the loop, you know their final size so before entering the loop try to create them and fill them with nan vallues.

function:nan(3,3)=[nan nan nan;nan nan nan;nan nan nan]

[Ronald van den Berg]

Luc, thanks for your response, but i'm afraid this won't help. Preallocation helps when variables change size within a loop, but m, nor Y, nor mhist does so (m and mhist are vectors, Y is the input matrix which is changed row by row in the loop).

Answer 1

[Sean de Wolski]

If you have MATLAB Coder, this could be a potentially good candidate for C code generation and MEXing. It involves a loop with a persistent state that's updated on each iteration, something compiled languages are often faster with.

If you don't, post a zip file with the matrix and I'll benchmark it for you.

[Ronald van den Berg]

Thanks, this is exactly what we're trying now. And it seems to be pretty fast (about to do a comparison now) -- i don't think any clever .m function would be able to beat this good old C solution :)

Answer 2

[Titus Edelhofer]

Hi,

one thing that comes to my mind: can you swap rows and columns, i.e., transpose your input matrix?

function Y = reflect(Y,b)
% apply reflecting bound to signals in matrix Y; columns are signals, rows time steps
mhist = zeros(size(Y,1),1);        % history of shifts applied in prev time steps
for ii=1:size(Y,2)                 % loop over time steps
  m = max(b-Y(:,ii)-mhist,0);       % compute required shift for current time step 
  Y(:,ii) = Y(:,ii) + mhist + m;    % apply shift
  mhist = mhist + m;                % add shift to history
end

That should be faster, since MATLAB likes columns better than rows ...

Titus

[Ronald van den Berg]

Oh, i didn't know that column operations are faster than row (this may have to do with the order in which matrix elements are stored in memory?). I'll try this. Even if it doesn't help here (since we have a mex solution now), this could be very useful to keep in mind for future functions. Thanks!

[Titus Edelhofer]

I'm curious to see the speed up you gain with the MEX/MATLAB Coder approach. And yes, matrices are stored columnwise in MATLAB (you see this e.g. with

A = magic(4);
A(:)

Therefore both from programming standpoint but also from memory access (cache!) working on columns is faster.

Titus