Moving window Autocorrelation. One pass algorithm. Speed and Stability concerns.

Hello,

I have a timeseries for which I want to generate a 'moving window autocorrelation at lag 1' series using a window of size window_size and a step of size window_step.

I have written the following one pass algorithm code for evaluating AC at lag 1:

function AC_timeseries = rolling_autocorr_lag1(x, window_size, window_step)
    
    N = length(x);
    
    sum_x = zeros(1, N);
    sum_x(1) = x(1);
    sum_xy = zeros(1, N);

    for n = 2: N
        sum_x(n) = sum_x(n - 1) + x(n);
        sum_xy(n) = sum_xy(n - 1) + x(n) * x(n - 1);
    end
    
    % Get the indices at which to calculate autocorrelation lag 1 values
    window_ends_idx = window_size: window_step: length(x);
    
    % Preallocate array for storing autocorrelation lag 1 time series
    AC_timeseries = zeros(1, length(window_ends_idx));

    % Generate autocorrelation lag 1 time series
    for i = 1: length(window_ends_idx)
        idx = window_ends_idx(i);
        if i == 1
            sum_use = sum_x(idx);
        else
            sum_use = sum_x(idx) - sum_x(idx - window_size);
        end
        sum_x_lag1 = sum_use - x(idx);
        sum_x_fwd1 = sum_use - x(idx - window_size + 1);
        sum_xy_window = sum_xy(idx) - sum_xy(idx - window_size + 1);
        std_dev_window = std(x(idx - window_size + 1: idx));
        N = window_size - 1;
        % Calculate autocorrelation lag 1 for the window
        AC_timeseries(i) = (N * sum_xy_window - sum_x_lag1 * sum_x_fwd1) / (N * std_dev_window * N * std_dev_window);
    end

end

My Concerns:

I would first like to know if this code can be improved further to speed up computations. Also, if a better method is available please let me know.

Also, I am sceptical that the variable sum_xy can overflow. How can I tackle this issue.

Is normalization of the data something that could help here? Can I make changes to the calculations/assignments to make this faster?

Answer 1

[Bruno Luong]

Do not put if inside loop, start the loop from i=2 and do i=1 outside

Do not compute inside the 2nd loop std(x(idx - window_size + 1: idx))

Rather compute sum(x.^2) in the first loop the use the formula

std(xw) = sqrt((sum(xw.^2) - sum(xw).^2 /window_size)/(window_size-1))

to get the std, with sum(xw) sum(xw.^2) values get from cumulative results and obtained with 2 substractions.

Normalization won't help.

You might want to quantify x as integer multiple of some fixed and small quantity. The result is approximative but round off error won't accumulate for long sequence.

You did not tell which release you are using. In recent version MATLAB for-loop with basic arithmetics is quite fast and usually C++-mex can be notably slower (useless at the moment), and C-mex won't have much chance to beat MATLAB.

[Bruno Luong]

This function rolling_autocorr_lag1_BLU seems much faster (27 times) and returns same result

x = randi([1e3 1e4],[1 10000]);
t = zeros(1,2);

window_size = 100;
window_step = 10;
t(1) = timeit(@() rolling_autocorr_lag1(x, window_size, window_step), 1);
t(2) = timeit(@() rolling_autocorr_lag1_BLU(x, window_size, window_step), 1);

% Check if the result match
AC = rolling_autocorr_lag1(x, window_size, window_step);
AC_BLU = rolling_autocorr_lag1_BLU(x, window_size, window_step);
relerr = norm(AC-AC_BLU) / norm(AC)

relerr = 2.7716e-16

tms = t*1e3

tms = 1×2

2.7576 0.0948

function AC_timeseries = rolling_autocorr_lag1(x, window_size, window_step)
    
    N = length(x);
    
    sum_x = zeros(1, N);
    sum_x(1) = x(1);
    sum_xy = zeros(1, N);
    for n = 2: N
        sum_x(n) = sum_x(n - 1) + x(n);
        sum_xy(n) = sum_xy(n - 1) + x(n) * x(n - 1);
    end
    
    % Get the indices at which to calculate autocorrelation lag 1 values
    window_ends_idx = window_size: window_step: length(x);
    
    % Preallocate array for storing autocorrelation lag 1 time series
    AC_timeseries = zeros(1, length(window_ends_idx));
    % Generate autocorrelation lag 1 time series
    for i = 1: length(window_ends_idx)
        idx = window_ends_idx(i);
        if i == 1
            sum_use = sum_x(idx);
        else
            sum_use = sum_x(idx) - sum_x(idx - window_size);
        end
        sum_x_lag1 = sum_use - x(idx);
        sum_x_fwd1 = sum_use - x(idx - window_size + 1);
        sum_xy_window = sum_xy(idx) - sum_xy(idx - window_size + 1);
        std_dev_window = std(x(idx - window_size + 1: idx));
        N = window_size - 1;
        % Calculate autocorrelation lag 1 for the window
        AC_timeseries(i) = (N * sum_xy_window - sum_x_lag1 * sum_x_fwd1) / (N * std_dev_window * N * std_dev_window);
    end
end


function AC_timeseries = rolling_autocorr_lag1_BLU(x, window_size, window_step)

    sum_x = [0 cumsum(x)];
    sum_x2 = [0 cumsum(x.^2)];
    sum_xy = [0, cumsum(x(2:end).*x(1:end-1))];
    
    % Get the indices at which to calculate autocorrelation lag 1 values
    window_ends_idx = window_size: window_step: length(x);
    
    % Generate autocorrelation lag 1 time series
    N = window_size - 1;
    idx = window_ends_idx;
    iwin1 = idx - window_size + 1;
    sum_use = sum_x(idx+1) - sum_x(iwin1);
    sum_x_lag1 = sum_use - x(idx);
    sum_x_fwd1 = sum_use - x(iwin1);
    sum_xy_window = sum_xy(idx) - sum_xy(iwin1);
    sum_x2_window = sum_x2(idx+1) - sum_x2(iwin1);
    var_dev_window = (sum_x2_window - sum_use.^2/window_size); % N * square of std_dev_window
    % Calculate autocorrelation lag 1 for the window
    AC_timeseries = (N*sum_xy_window - sum_x_lag1.*sum_x_fwd1) ./ (N * var_dev_window);

end

[atharva aalok]

@Bruno Luong Unbelievable stuff! This is blazingly fast. Using my above function I was getting around a 150x boost in speed from my original naive implementation. I am getting another 20x through your function over mine.

So I am now able to get a whopping 3000x+ boost thanks to you!!!

[Bruno Luong]

Code on wikipedia is harder to beat. ;-)

PS: don't forget to look at my last edit where the entire loop is removed.

[atharva aalok]

@Bruno Luong I asked the same question on stackoverflow as well. If you would like to, you can provide your answer there as well and I will accept that. For now, I have linked your answer there.

Answer 2

[Matt J]

y=circshift(x,-window_step);

w=normalize(ones(1,window_size),'n',1);

Ex=cyconv(x,w);
Ey=cyconv(y,w);
Exy=cyconv(x.*y,w);

AC_timeseries=Exy-Ex.*Ey;


function z=cyconv(x,y)
%Non-Fourier domain cyclic convolution
%
%  z=cyconv(x,y)

 siz=num2cell(size(x));

 subs=cellfun(@(n)[2:n,1:n],siz,'uni',0);

 x=x(subs{:});

 z=convn(x,y,'valid');
end