Applying correction to a time series in Matlab

Question

I face a computation efficiency problem. I have a time series of a non-monotonically drifting variable, that is measurements of objects going through a machine (where the measurement is made) in a production line. The job to do consists of simulating what this time series would yield if there was a correction made to the object each time the measurement drifts above or below a threshold.

To do that I could simply make a for loop, and each time the thresholds are crossed, apply the correction to the rest of the time series. However the time series is very long and the for loop would take too much time to compute. I would like to improve the performance of my code. Does anyone see a way to do this?

Here is a working example using a for loop:

ts = [1 2 3 4 3 2 1 0 -1 -2 -3];
threshold = [-3.5 3];
correction = [1 -1];
for i = 1:numel(ts)
    if ts(i) > threshold(2)
       ts(i:end) = ts(i:end) + correction(2);
   elseif ts(i) < threshold(1)
       ts(i:end) = ts(i:end) + correction(1);
   end
end
disp(ts)

The result is and should be the following array:

[1 2 3 3 2 1 0 -1 -2 -3 -3]

Answer 1

There is really nothing wrong with using a loop in MATLAB. For decades, they have been teaching us to vectorize our code. But in recent years the differences between vectorized code and the equivalent loop code has been shrinking dramatically. And I've run into examples where vectorized code is actually slower!

Nonetheless, there is a simple improvement to be made to your code. Note that what you implemented is an algorithm of \$O(n^2)\$. This simplified bit:

for i = 1:numel(ts)
   ts(i:end) = ts(i:end) + correction(2);
end

does \$n(n+1)/2\$ additions (hence \$O(n^2)\$). Also the indexing is expensive, since you're copying half the array (on average) every loop iteration.

I suggest that, instead of adding correction(2) to all subsequent elements, you store the "current correction", update it every loop iteration, and add it to the current value only.

Below is the test function I wrote. method1 is your code, method2 is my suggestion. At the top is a function that exercises these two methods, compares their output, and times them:

function test_methods
ts = [1 2 3 4 3 2 1 0 -1 -2 -3];
ts1 = method1(ts);
ts2 = method2(ts);
if any(abs(ts1-ts2)>1e-6)
   error('the methods differ');
end
timeit(@()method1(ts))
timeit(@()method2(ts))

function ts = method1(ts)
threshold = [-3.5 3];
correction = [1 -1];
for i = 1:numel(ts)
   if ts(i) > threshold(2)
      ts(i:end) = ts(i:end) + correction(2);
   elseif ts(i) < threshold(1)
      ts(i:end) = ts(i:end) + correction(1);
   end
end

function ts = method2(ts)
threshold = [-3.5 3];
correction = [1 -1];
current = 0;
for i = 1:numel(ts)
   value = ts(i) + current;
   if value > threshold(2)
      current = current + correction(2);
      value = value + correction(2);
   elseif value < threshold(1)
      current = current + correction(1);
      value = value + correction(1);
   end
   ts(i) = value;
end

For the short example input, these functions both run too fast for accurate timing. I see 3.1 μs and 1.4 μs, my version is only twice as fast as yours. But for larger inputs the differences become more important (I figured that the cumulative sum of a random process would imitate appropriately your drifting variable):

ts = cumsum(randn(1,1000));

0.29 ms and 8.33 μs, an order of magnitude difference.

ts = cumsum(randn(1,100000));

Now I see 1.66 s and 0.812 ms, 3 orders of magnitude difference.

Because method1 is quadratic in the input length, and method2 is linear, the time difference grows quadratically.