How to get rid of a loop that depends on its previous iterations

I have an old c++ program that I want to translate to Matlab. It has a loop in which every iteration depends on a decision in the previous iteration. I have a lot of data, so it takes a long time. Is there a way to do this with matrix operations instead of a loop? Here's the code:

data = rand(1,105);
for k = 100 : -1 : 6
    lowerAverage = mean(data(k-5:k-1));
    upperAverage = mean(data(k+1:k+5));
    surroundingAverage = mean([lowerAverage,upperAverage]);
    if data(k) > surroundingAverage
        data(k) = upperAverage
    end
end

Thanks in advance.

Answer 1

[Jan]

Start with calculating the average manually instead of the slower mean() function:

% Timings: R2018b, Win10, i7
data0 = rand(1, 100005);
data  = data0;
tic;
for k = 100000 : -1 : 6
    lowerAverage = mean(data(k-5:k-1));
    upperAverage = mean(data(k+1:k+5));
    surroundingAverage = mean([lowerAverage,upperAverage]);
    if data(k) > surroundingAverage
        data(k) = upperAverage;
    end
end
toc  % Elapsed time is 1.050304 seconds.
data_orig = data;

data = data0;
tic;
for k = 100000 : -1 : 6
    lowerAverage = sum(data(k-5:k-1)) / 5;
    upperAverage = sum(data(k+1:k+5)) / 5;
    surroundingAverage = (lowerAverage + upperAverage) / 2;
    if data(k) > surroundingAverage
        data(k) = upperAverage;
    end
end
toc  % Elapsed time is 0.050326 seconds.

isequal(data, data_orig)  % true

20 times faster already. Mayby this is some percent faster:

data = data0;
tic;
for k = 100000 : -1 : 6         % 1 multiplication, 1 division:
    low = sum(data(k-5:k-1));
    up  = sum(data(k+1:k+5));
    if 10 * data(k) > (low + up)
        data(k) = up / 5;
    end
end
toc  % Elapsed time is 0.049661 seconds.

The lower average is not influenced by overwriting the data. So this can be vectorized:

data = data0;
tic;
lowv = conv(data, ones(1, 5), 'same');
for k = 100000 : -1 : 6
    low = lowv(k - 3);
    up  = sum(data(k+1:k+5));
    if 10 * data(k) > (low + up)
        data(k) = up / 5;
    end
end
toc  % Elapsed time is 0.025510 seconds.

The last method takes significantly more time when I run it in the forums MATLAB online 2021b. So try this locally on your computer.

[EDITED] If replacing mean() increases the speed, try this with sum() also:

data = data0;
tic;
avg = conv(data, ones(1, 5), 'same');
for k = 100000 : -1 : 6
   low = avg(k - 3);
   up  = data(k+1) + data(k+2) + data(k+3) + data(k+4) + data(k+5);
   if 10 * data(k) > (low + up)
      data(k) = up / 5;
   end
end
toc  % Elapsed time is 0.007511 seconds.

A speedup of factor 140 on my R2018b and the random test data.

[JPL]

Thanks Jan,

I was afraid that would be the conclusion. I'll add your suggestions on making it faster, and just learn to live with the loop.

V/R,

John

[Jan]

My first suggestions runs at least 40 times faster than the original loop.

See [EDITED] in my answer which is 140 times faster.

[JPL]

Jan,

This will save days of processing. I'm also going to try to write a mex function and see if I can squeeze some more time out of it.

Thanks again,

JPL

[Jan]

A mex version needs about half of the processing time compared to the last M-version of my asnwer:

#include "mex.h"
mvoid mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  double *X, low, up;
  mwSize i, n;
  
  n = mxGetNumberOfElements(prhs[0]);
  plhs[0] = mxDuplicateArray(prhs[0]);
  X = mxGetDoubles(plhs[0]);
  
  for (i = n - 6; i > 4; i--) {
     low = X[i-5] + X[i-4] + X[i-3] + X[i-2] + X[i-1];
     up  = X[i+1] + X[i+2] + X[i+3] + X[i+4] + X[i+5];
     if (10 * X[i] > low + up) {
        X[i] = up / 5;
     }
  }
  
  return;
}

[JPL]

Thanks,

I think this will work out well.

John