Quicker way to subtract large vectors using a sliding window?

Below is code with comments for two versions of vector subtraction using a sliding window. I got the math to work (first method for each version) and then was able to optimize it (second method for each version). However, I wonder if there is an even faster way to accomplish this sort of calculation, perhaps a vectorized method. I realize that this post is lengthy, but I figured that it would be helpful to clearly state the entire problem. Any suggestions or feedback are appreciated.

First, variables are initialized for both versions:

% initialize variables
rng(42);
nr = 100000;
A = rand(nr,1);
B = rand(nr,1);
w = 10; % window length
D1 = zeros(nr,1); % initialize results vector
D2 = zeros(nr,1); % initialize results vector

This is the first example, wherein a window slides over the first vector to perform subtraction at each value in the second vector, which results are then summed to populate the ith value. Two loop methods are presented and timing results are displayed on running the code.

%% Example 1: sliding window over one vector
%% First loop method
tic
for t_loop = 10; % is typically 1,000's to 1,000,000s for my application

for i = w:nr
    D1(i) = sum(A((i-w+1):i) - B(i));
end

end
t1_1st = toc

t1_1st = 0.0233

%% Second loop method
tic
for t_loop = 10; % is typically 1,000's to 1,000,000s for my application

A_arr = zeros(nr,w);
for i = 1:w
    A_arr(w:end,i) = A(i:(end-w+i));
end
D2 = sum( bsxfun(@minus, A_arr, B), 2);
D2(1:(w-1)) = 0; % mimic initialization in D1

end
t1_2nd = toc

t1_2nd = 0.0085

%% Check that the results are the same
assert(sum(abs(D1-D2)>1e-6)==0)
disp(['Method 2 is ',num2str(t1_1st / t1_2nd),' times faster than method 1.'])

Method 2 is 2.7574 times faster than method 1.

The second method is two to three times faster than the first method (the amount varies each time I run it), which is very helpful when running thousands to millions of iterations for combinatoric studies. But, as stated above, I would like to know if this can be optimized further.

This is the second example, wherein windows slide over both vectors to perform similar subtractions.

%% Example 2: sliding window over both vectors
D1 = zeros(nr,1); % initialize results vector
D2 = zeros(nr,1); % initialize results vector

%% First loop method
tic
for i = w:nr
    s = 0; % initialize sum for the window at each iteration
    for j = 1:(w-1)
        s = s + sum(A(i-j) - B((i-j+1):i));
    end
    D1(i) = s;
end
t2_1st = toc

t2_1st = 0.1753

%% Second loop method
tic
B_arr = zeros(nr,w-1);
for i = 1:(w-1)
    B_arr_slice = zeros(nr-w+1,i);
    for j = 1:i
        B_arr_slice(1:end,j) = B((w-j+1):(end-j+1));
    end
    B_arr(w:end,i) = sum( bsxfun(@minus, A((w-i):(end-i)), B_arr_slice), 2);
end
D2 = sum(B_arr,2);
t2_2nd = toc

t2_2nd = 0.0255

%% Check that the results are the same
assert(sum(abs(D1-D2)>1e-6)==0)
disp(['Method 2 is ',num2str(t2_1st / t2_2nd),' times faster than method 1.'])

Method 2 is 6.8863 times faster than method 1.

This second method is much faster, as expected, given the improvement for the single-window case and that the sliding window operates over both vectors here.

However, as mentioned previously, these calculations are performed thousand or millions of times. Every additional optimization would be very helpful. And, I can't help but wonder if there is a faster and/or vectorized method of performing these calculations.

Answer 1

[Paul]

Hi goc3,

Using filter seems to be a lot faster than the loop method 1 for Example 1. I didn't compare to the other loop method for Example 1, or the different operation in Example 2.

rng(42);
nr = 100000;
A = rand(nr,1);
B = rand(nr,1);
w = 10; % window length
D1 = zeros(nr,1); % initialize results vector
D2 = zeros(nr,1); % initialize results vector

tic
% corrected this line, loop should be 1:10, not just 10
for t_loop = 1:10; % is typically 1,000's to 1,000,000s for my application 
for i = w:nr
    D1(i) = sum(A((i-w+1):i) - B(i));
end
end
t1_1st = toc

t1_1st = 0.1981

tic
for t_loop = 1:10
    D3 = filter(ones(1,w),1,A) - w*B;
    D3(1:w-1) = 0;
end
t3 = toc

t3 = 0.0077

max(abs(D1-D3))

ans = 3.5527e-15

[goc3]

@Paul: your solution is the fastest option so far for Example 1 by an order of magnitude.

Could you also provide a solution for Example 2? That one calculates differences using two sliding windows instead of just one. It is not immediately apparent to me how to extend your solution for Example 1 to Example 2.

[Paul]

I think this does it.

rng(42);
nr = 100000;
A = rand(nr,1);
B = rand(nr,1);
w = 10; % window length
D1 = zeros(nr,1); % initialize results vector
D2 = zeros(nr,1); % initialize results vector

%% First loop method
tic
for t_loop = 1:10
for i = w:nr
    s = 0; % initialize sum for the window at each iteration
    for j = 1:(w-1)
        s = s + sum(A(i-j) - B((i-j+1):i));
    end
    D1(i) = s;
end
end
t2_1st = toc

t2_1st = 1.6746

tic
for t_loop = 1:10
D3 = filter(0:(w-1),1,A) - filter((w-1):-1:1,1,B); % corrected for general case based on comment below
D3(1:w-1) = 0;
end
t3 = toc

t3 = 0.0159

max(abs(D3-D1))

ans = 1.5987e-14

[goc3]

Your answer to Example 2 provides even more of an optimization than the one for Example 1. Thanks, for your prompt responses and answers, @Paul.

For anyone else who this might help, the two 9's in the code above (0:9 and 9:-1:1) should be replaced by w-1 to support the general case [i.e., 0:(w-1) and (w-1):-1:1].

[Paul]

I was so focused on identifying the pattern for the specific case that I forgot to go back to the general case. I'll correct it.

Answer 2

[Nipun]

Hi,

I understand that you are looking for a quicker way to calculate the difference of large vectors using sliding window in MATLAB.

Based on the provided code and the attached code snippet, there is indeed a potential for further optimization through vectorization, which is one of MATLAB's strengths. I will split this answer into two parts: vectorized method for example 1 and vectorized method for example 2.

Vectorized method for example 1

The operation being performed in Example 1 is essentially a moving sum of the differences between a window in vector A and a corresponding value in vector B. This can be vectorized using the cumsum function (cumulative sum) and then subtracting the cumulative sum at two points to get the sum over a window.

Here's a vectorized approach for Example 1:

%% Example 1: Vectorized approach
tic
cumsumA = cumsum(A); % Cumulative sum of A
cumsumA = [zeros(w-1, 1); cumsumA]; % Pad with zeros to align with the window size
D1_vect = cumsumA(w:end) - cumsumA(1:end-w+1) - w * B(w:end);
D1_vect = [zeros(w-1, 1); D1_vect]; % Pad with zeros for the initial values
t_vect = toc;

% Check that the results are the same
assert(sum(abs(D1-D1_vect)>1e-6)==0)
disp(['Vectorized method for Example 1 is ',num2str(t1 / t_vect),' times faster.'])

Vectorized method for example 2

In Example 2, a double sliding window operation is being performed. The vectorization of such an operation can be quite tricky, but it's still possible to optimize it further by eliminating the innermost loop and performing operations on larger chunks of data.

Here's a vectorized approach for Example 2:

%% Example 2: Vectorized approach
tic
% Precompute the cumulative sums for A and B
cumsumA = cumsum([0; A]);
cumsumB = cumsum([0; B]);

% Initialize the result vector
D2_vect = zeros(nr, 1);

% Compute the result using vectorized indexing and operations
for i = w:nr
    A_sum = cumsumA(i:-1:i-w+1) - cumsumA(i-1:-1:i-w);
    B_sum = cumsumB(i:-1:i-w+2) - cumsumB(i-1:-1:i-w+1);
    D2_vect(i) = sum(A_sum(2:end) - B_sum(1:end-1));
end
t_vect = toc;

% Check that the results are the same
assert(sum(abs(D2-D2_vect)>1e-6)==0)
disp(['Vectorized method for Example 2 is ',num2str(t2 / t_vect),' times faster.'])

Link to documentation:

1. cumsum : Cumulative sum - MATLAB cumsum - MathWorks India

Hope this helps.

Regards,

Nipun

[goc3]

@Nipun, this does look like a good approach. However, I copied the code from your answer and could not get it to run or produce correct results. Could you edit the code in your answer so that it runs without error?

I made a few minor edits in my question so that all four timing variables are accessible. (Before, I used two variables that were overwritten for Example 2.)