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

Question

atharva aalok 2023-8-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010292-moving-window-autocorrelation-one-pass-algorithm-speed-and-stability-concerns

评论： atharva aalok 2023-8-19

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2023-8-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010292-moving-window-autocorrelation-one-pass-algorithm-speed-and-stability-concerns#answer_1289707

编辑：Bruno Luong 2023-8-19

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Bruno Luong 2023-8-19

编辑：Bruno Luong 2023-8-19

在 MATLAB Online 中打开

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

Bruno Luong 2023-8-19

编辑：Bruno Luong 2023-8-19

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 2023-8-19

@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 2023-8-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010292-moving-window-autocorrelation-one-pass-algorithm-speed-and-stability-concerns#answer_1289737

编辑：Matt J 2023-8-19

在 MATLAB Online 中打开

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