Quicker way to subtract large vectors using a sliding window?
3 views (last 30 days)
Show older comments
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
%% 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
%% 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.'])
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
%% 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
%% 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.'])
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.
0 Comments
Accepted Answer
Paul
on 31 Dec 2023
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
tic
for t_loop = 1:10
D3 = filter(ones(1,w),1,A) - w*B;
D3(1:w-1) = 0;
end
t3 = toc
max(abs(D1-D3))
4 Comments
Paul
on 31 Dec 2023
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.
More Answers (1)
Nipun
on 29 Dec 2023
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:
Hope this helps.
Regards,
Nipun
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!