Is Conv2 based on direct convolution?

Is the conv2 function in Matlab 2023a, in particular the way it operates on 1d vectors, based on direct convolution or fft? I understand that you do not want to share all details on how this function is implemented, but for making use of this function in a scientific work where round-off errors and computational cost enter, it would be good to state with confidence that the method is either direct-based or fft-based.

Answer 1

[Matt J]

It is not FFT-based.

To do FFT-based convolution (using the overlap-and-add method) you could use fftfilt.

[Håkon Hoel]

Thank you for the reply. Does that mean that it is direct, in the sense that for

v = conv(f,g) (which again calls conv2(f,g))

the method computes the i-th element by v(i) = sum_{j=1}^i f(j) g(i+1-j) (*)

through additions and multiplications of the involved vectors? Surely avoiding unneccessary repetitions etc, but I mean only through using these two operations and the values in the vectors, no other transformation methods or similar that may make a discussion on round-off and complexity more difficult.

[Matt J]

As I've always understood it yes, the evaluation is through implementation of equation (*). However, the implementation is multithreaded by Matlab's internal parallelization, so that could complicate a dicussion of error accumulation and complexity. In particular, I don't know you would find that the runtime of conv is O(N^2), as an example, though we could certainly try it.

N=[1000:1000:100000];
for i=1:numel(N)

    x=rand(1,N(i));
    tic
     conv(x,x);
    t(i)=toc;
    
end

p=polyfit(N,t,2)

p = 1×3

1.0e-03 * 0.0000 0.0001 -0.2466

NN=linspace(min(N),max(N),1000);
plot(N,t,'o',NN, polyval(p,NN)); xlabel N; ylabel time(sec)

[Håkon Hoel]

I agree the cpu-timing shows that conv is not implemented using FFT. I have observed similarly, and I have also observed that the result of calling conv(f,g) is almost identical to the result of a direct computation of (*) with a "naive" double for-loop, while it is not always almost identical to using FFT-based convolution. That is, when you convolve vectors with all elements smaller than the double precision eps_double in f and g, then FFT seems to introduce round-off errors of magintude eps_double, making the relative error of the convolution very large. This does not happen with conv(f,g), which in my experience manages to keep the relative error small also in such settings, as one would expect from a direct convolution method. My issue is that these are all indications and I would like a definite yes or no, but your insight was nonetheless helpful!

[Christine Tobler]

I can confirm that conv2 uses a direct convolution method and not fft, and have passed along the suggestion to make this explicit in the documentation.

Answer 2

[Walter Roberson]

Historically, MATLAB implemented several different cases according to size and data type (and matters such as complex-valued or not.) The general case was done by direct convolution -- but possibly by calls into high speed libraries. The special cases involved decomposition into kernels and some kind of transform.

The source code for the direct convolution used to be visible. The decomposition etc., was something I do not recall ever seeing the source code for (that is, it might have been present in MATLAB source form in sufficiently old versions in the time before I thought to look for it.)

[Walter Roberson]

Note: historically, conv2() was implemented as a series of calls to conv() with transposes between.

[Håkon Hoel]

Thank you for complementing the above answer! I would accept this answer too, if possible. I believe that parallelization would not produce the round-off errors I have in mind, as they relate to the relative error blowing up when eps_double \approx 10^{-16} is large in magnitude relative to all elements in the vectors being convolved. In such cases my observation is that round-off errors from FFT-based convolution can lead to O(1) relative error while conv2() manages to maintain a low relative error.