Convolution using non-constant kernel?

Is there a faster way of performing convolution using a non-constant kernel than with a loop? The kernel contains values completely independent from those in the first vector.

I have been able to achieve what I need by code similar to what follows, but it is being run millions of times and I would love to optimize it, if possible.

rng(42);
nr = 100000;
A = rand(nr,1);
B = rand(nr,1);
w = 100; % window length
D = NaN(nr,1); % initialize results vector

for i = w:nr
    D(i) = A((i-w+1):i)' * B((i-w+1):i);
end

The following was added on 07 Oct 2023:

This small example illustrates that the kernel of convolution, represented here by B, changes along with the sliding window.

nr = 10;
A = [5, 3, 2, 7, 8, 6, 4, 5, 8, 6];
B = [1, 2, 3, 2, 3, 2, 1, 2, 2, 3];
w = 3; % window length
D = NaN(1, nr); % initialize results vector

for i = w:nr
    D(i) = A((i-w+1):i) * B((i-w+1):i)'; % values in B change for each i
end
disp(D)

   NaN   NaN    17    26    44    50    40    26    30    44

This results in the following:

D = [NaN, NaN, 17, 26, 44, 50, 40, 26, 30, 44]

Answer 1

[Paul]

Hi goc3,

filter provides significant speedup

rng(42);
nr = 100000;
A = rand(nr,1);
B = rand(nr,1);
w = 100; % window length

tic
for count = 1:100
D = NaN(nr,1); % initialize results vector
for i = w:nr
    D(i) = A((i-w+1):i)' * B((i-w+1):i);
end
end
toc

Elapsed time is 3.781300 seconds.

tic
for count = 1:100
    % assumes A is real (original code uses conjugate transpose)
    D1 = filter(ones(1,w),1,A.*B,nan(1,w-1));
end
toc

Elapsed time is 0.058871 seconds.

isequal(isnan(D),isnan(D1))

ans = logical
   1

plot(abs(D-D1))

Answer 2

[Shivansh]

Hi goc3,

Yes, there is a faster way to perform convolution using a non-constant kernel without using a loop. You can use the conv() function in MATLAB to achieve this. Here's an example of how you can modify your code to utilize conv:

rng(42);

nr = 100000;

A = rand(nr,1);

B = rand(nr,1);

w = 100; % window length

kernel = B(w:-1:1); % Reverse the kernel to match the order used in the loop

% Perform convolution using the 'valid' option to match the loop implementation

D = conv(A, kernel, 'valid');

In this code, we create a non-constant kernel as a vector kernel. You can replace it with your actual kernel values. Then, we use the conv() function with the 'valid' option, which computes only the valid part of the convolution and discards any overlapping edges. This is equivalent to the operation you were performing in your loop.

Using conv should be significantly faster than using a loop for large arrays, as it takes advantage of optimized algorithms implemented in MATLAB. For more information on ‘conv’, you can refer to the documentation here https://in.mathworks.com/help/matlab/ref/conv.html?searchHighlight=conv&s_tid=srchtitle_support_results_1_conv.

Hope it helps!

[goc3]

Shivansh, I appreciate your response. However, it does not address the stated problem. I just added a small example to my question above to illustrate.

The answer you provided uses the 'valid' option; this serves to remove the first two (NaN) values from the D vector shown above. This is not intended.

Furthermore, the kernel definition here is extracted once from the B vector, essentially becoming a constant kernel. When I ran your code against the small example I just added, it yielded:

D = [17 28 40 41 32 29 38 39]

which is not equal to the expected answer for D above.