How to speed up convolution with a million data points

I am currently doing convolution using nested for loops, for 10^6 data points in each for loop. Are there ways to speed up the following code? Thanks in advance!

% nIters = 40;
% n = 1e6; 
% mzL = rand(nIters, n); 
% gg  = rand(1, n);
mzR_temp = zeros(nIters, 1);
for c = 1:n 
    mzR_temp(:) = 0;
    for d = 1:c
        mzR_temp(:) = mzR_temp(:) + gg(c-d+1) * mzL(:,d);
    end
    
    mzR_II(:,c) = mzR_temp;
end

Answer 1

[Matt J]

Use conv,

mzR_II=conv(gg,mzL,'same');

or FFTs,

mzR_II=ifft( fft(gg,2*n) .* fft(mzL,2*n) , 'symmetric');

mzR_II=mzR_II(1:n);

[Paul]

It appears that conv has to be 'full' not 'same', after and then truncated.

n = 5;
rng(100)
mzL = rand(1,n);
gg = rand(1,n-1);
for c = 2:n 
    mzR_temp = 0;
    for d = 1:(c-1)
        mzR_temp = mzR_temp + gg(c-d) * mzL(d);
    end
    
    mzR_II(c) = mzR_temp;
end
mzR_II

mzR_II = 1x5

0 0.0661 0.3983 0.6871 0.6916

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[0 conv(gg,mzL,'same')] 

ans = 1x5

0 0.6871 0.6916 0.9559 0.7589

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[0 conv(gg,mzL,'full')]

ans = 1x9

0 0.0661 0.3983 0.6871 0.6916 0.9559 0.7589 0.1194 0.0006

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fftfilt in the Signal Processing Toolbox can be used for the fft approach.

[0 fftfilt(mzL,gg)]

ans = 1x5

0 0.0661 0.3983 0.6871 0.6916

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Or filter

[0 filter(mzL,1,gg)]

ans = 1x5

0 0.0661 0.3983 0.6871 0.6916

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Runzi Hao]

Cool! And I agree 'full' should be used.

However, due to large number of data points, my script still runs slow - aprroximately 40 min used. My computer is not super good, but is it possible to further speed up the calculation?

Besides, I would do convolution of gg(1:n) and mzL(i, 1:n) for i = 1:40, which actually will create 40 * 10^6 convolution data points.

[Matt J]

There is absolutely no way the computation should take more than 1 second on any computer made within the last 10 years.

n=1e6;
mzL = rand(1,n);
gg = rand(1,n-1);

tic
mzR_II = [0 fftfilt(mzL,gg)];
toc

Elapsed time is 0.174552 seconds.

tic;
mzR_II=ifft( fft(gg,2*n) .* fft(mzL,2*n) , 'symmetric');
mzR_II=mzR_II(1:n);
toc

Elapsed time is 0.145911 seconds.

[Runzi Hao]

Awesome! Thanks for the suggestion of using fftfilt.

It turns out that, on my computer, the convolution ran 43 s for a million data points; and with 40 iterations, the total time was 2600+ s. However, the fftfilt function ran only 8 s for the 40 million data points!

Answer 2

[Image Analyst]

Use the built-in convolutions functions: conv and conv2. They are highly optimized for speed.