# Applying correction to a time series in Matlab

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]

• Asking for advice on code yet to be written or implemented is off-topic for this site. See What topics can I ask about? for reference. Once you have written working code, you're welcome to ask a new question here and we can then help you improve it! – Phrancis Apr 19 '18 at 18:02
• This question is off-topic because you have reviewed your own code, and decided that it needs to be vectorized to make it work better. Once you have vectorized the code you can bring that code back for review. Alternatively, we can review the code you currently have, but reviewers may just recommend that you vectorize it - which you already know. Note, we review the code you have, not the code you want to have. – rolfl Apr 19 '18 at 18:13
• This question has been mentioned on Meta. – 200_success Apr 19 '18 at 18:30
• It's quiet funny because I posted that question on stackoverflow first and I they said it's too broad. You are not pointing at a specific problem here. you should go on CR. And here it is more or less the counter part of this. My question is to precise. I should ask for "please what do you think of this code, generaly speaking?". No place for my question, too bad. I'll work that on my own. Thanks anyway. – Cunningham Apr 21 '18 at 10:16
• @Cuningham We've been discussing it on meta, haven't gotten a complete consensus maybe, but I've edited your question to make it match what other people thing you should ask (I personally don't). I know that there are some overlap between SO and CR, and also confusion in what's on-topic where and why. I've reopened your question after editing it, as improving performance is your ultimate goal - vectorization might be one possibility for that. – Simon Forsberg Apr 21 '18 at 15:06

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.

• Thanks a lot! I almost cried when I realized that vectorization is not a relevant practice anymore. I've spend so much time trying to vectorise my code systematically.... But all in all it's a good news :). – Cunningham Apr 26 '18 at 17:16
• @Cuningham: Vectorization is still applicable, it just doesn't yield the huge benefits it used to. Often times vectorized code is more readable, and that is the most important thing in the world. :) – Cris Luengo Apr 26 '18 at 17:27