Why is using python scipy's convolve function on a series of vector pairs in a for loop faster than using fftconvolve on two equivalent matrices?

To give some background, I am performing matched filtering between two batches of signals. The transmitted signals (i.e. the kernel) and the received signals (the 'data') are then represented by two numpy arrays, each of size: rows = 50,000 columns = 64. (That is to say, both the transmit and receive arrays each contain 64 signals, where each signal is sampled 50,000 times, and placed in a column in one of the arrays, depending on whether it was a transmitted signal or a received signal).

For a given transmitted and received signal, (i.e. two vectors, each of length 50,000) the matched filter operation I'm using would look like:

convolve(reversed(transmit signal_vector), received_signal_vector)

In order to codify this (using python), for my two arrays (which contain the data from 64 signals each), I have done the following:

from scipy import signal

for column in range(64):
    output_array[:,column] = signal.convolve(
        np.conj(np.flipud(transmitted_signal_array[:,column])), 
        received_signal_array[:,column], mode='full', method='auto')

So essentially, every column in the transmitted_signal_array is (flipped and conjugated) and convolved with the corresponding column from received_signal_array. Note also, the method is always the fft method, not direct. The for loop sequentially goes through each column in the arrays such that each of the 64 pairs of transmitted and received signals are convolved without regard for the other signals contained in the arrays.

This seems to work well, but I wanted to get some speed up since it's a clear bottleneck when running a profiler. I imagined that performing the convolution using some in-built functionality rather than my for loop would probably be a good place to start but as far as I can tell the scipy.signal.convolve function does not allow convolution between two arrays and specification about which axis to concentrate on. However, the scipy.signal.fftconvolve function does. So I tried:

output = signal.fftconvolve(np.conj(np.flipud(transmitted_signal_array)),received_signal_array,mode='full',axis=0)

However, this performs slower than the original for loop method I was using. I understand that fft based convolution is not always the best way to implement it; however, as I previously mentioned, the scipy.signal.convolve function was ALWAYS choosing the fft method over the direct method anyway. Not only this, but if I directly substitute fftconvolve into my for loop operation, I get approximately the same time as when using convolve, which leads me to believe the reason fftconvolve is performing slower when being used with the entire arrays rather than sequentially on each pair of columns via the for loop is because something other than convolving each column from array1 with each column from array2 is occurring.

So why is my for loop method working faster than the fftconvolve method (with an axis specified as an argument to the function) and is there a faster python convolution method to do what I want (perhaps from a different library)?

And, perhaps I should ask this in a separate question, but, is there a smarter way to do what I want to do? Ultimately, this convolution operation is part of a simulation application which will use large quantities of repetitions of this operation. On the scales I'm working with, would it make sense to move arrays to a GPU, perform the convolution, and move the result back again?

EDIT:

The methods used are entirely from numpy: numpy.conj and numpy.flipud or scipy: scipy.signal.convolve and scipy.signal.fftconvolve

For inputs, the transmitted_signal_array could be:

X = 1
transmitted_signal_array = numpy.random.normal(0,X,size=(50000,64))+1j*numpy.random.normal(0,X,size=(50000,64))

and the received_signal_array could then be:

noise_array = numpy.random.normal(0,X,size=(50000,64))+1j*numpy.random.normal(0,X,size=(50000,64))
received_signal_array  = transmitted_signal_array + noise_array

For timing, I was using a profiler:

import cProfile
import pstats
import io 

def Profile(function):
   def inner(*args,**kwargs):
       pr = cProfile.Profile()
       pr.enable()
       retval = function(*args,**kwargs)
       pr.disable()
       s = io.StringIO()
       sortby = 'cumulative'
       ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
       ps.print_stats('My_Convolution_Method')
       with open('timing_results.txt', 'w+') as f:
          f.write(s.getvalue())
       return retval
   return inner

This is then used as a decorator for the function 'My_Convolution_Method'