Python numpy.fft changes strides
<p>Dear stackoverflow community!</p>

<p>Today I found that on a high-end cluster architecture, an elementwise multiplication of 2 cubes with dimensions 1921 x 512 x 512 takes ~ 27 s. This is much too long since I have to perform such computations at least 256 times for an azimuthal averaging of a power spectrum in the current implementation. I found that the slow performance was mainly due to different stride structures (C in one case and FORTRAN in the other). One of the two arrays was a newly generated boolean grid (C order) and the other one (FORTRAN order) came from the 3D <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftn.html#numpy.fft.fftn" rel="nofollow">numpy.fft.fftn()</a> Fourier transform of an input grid (C order). Any reasons why numpy.fft.fftn() changes the strides and ideas on how to prevent that except for reversing the axes (which would be just a workaround)? With similar strides (<a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.copy.html" rel="nofollow">ndarray.copy()</a> of the FT grid), ~ 4s are achievable, a tremendous improvement.</p>

<p>The question is therefore as following:</p>

<p>Consider the array:</p>

<pre><code>ran = np.random.rand(1921, 512, 512)
ran.strides
(2097152, 4096, 8)

a = np.fft.fftn(ran)
a.strides
(16, 30736, 15736832)
</code></pre>

<p>As we can see the stride structure is different. How can this be prevented (without using a = np.fft.fftn(ran, axes = (1,0)))? Are there any other numpy array routines that could affect stride structure? What can be done in those cases?</p>

<p>Helpful advice is as usual much appreciated!</p>

<blockquote>
  <p>Any reasons why numpy.fft.fftn() changes the strides and ideas on how to prevent that except for reversing the axes (which would be just a workaround)?</p>
</blockquote>

<p>Computing the multidimensionnal DFT of an array consists in successively computing 1D DTFs over each dimensions. There are two strategies:</p>

<ol>
<li>Restrict 1D DTF computations to contiguous 1D arrays. As the array is contiguous, problem related to latency/cache misses will be reduced. This strategy has a major drawback: the array is to be transposed between each dimension. It is likely the strategy adopted by <code>numpy.fft</code>. At the end of computations, the array has been transposed. To avoid unnecessary computations, the transposed array is returned and strides are modified.</li>
<li>Enable 1D DDFT computations for strided arrays. This might trigger some problem related to latency. It is the strategy of <code>fftw</code>, avaible through the interface <code>pyfftw</code>. As a result, the output array features the same strides as the input array.</li>
</ol>

<p>Timing <code>numpy.fftn</code> and <code>pyfftw.numpy.fftn</code> as performed <a href="http://stackoverflow.com/questions/40061307/comparatively-slow-python-numpy-3d-fourier-transformation/40064319?noredirect=1#comment67488042_40064319">here</a> and <a href="http://stackoverflow.com/questions/6365623/improving-fft-performance-in-python">there</a> or <a href="https://gist.github.com/fnielsen/99b981b9da34ae3d5035" rel="nofollow">there</a> will tell you whether FFTW is really the Fastest Fourier Transform in the West or not...</p>

<ul>
<li><p>To check that numpy uses the first strategy, take a look at <code>numpy/fft/fftpack.py</code>. At line 81-85, the call to <code>work_function(a, wsave)</code> (i.e. <code>fftpack.cfftf</code>, from <a href="http://www.netlib.org/fftpack/" rel="nofollow">FFTPACK</a>, arguments documented <a href="https://docs.oracle.com/cd/E19422-01/819-3691/cfftf.html" rel="nofollow">there</a>) is enclosed between calls to <code>numpy.swapaxes()</code> performing the transpositions.</p></li>
<li><p><code>scipy.fftpack.fftn</code> does not seem to change the strides... Nevertheless, it seems that it makes use of the first strategy. <a href="https://github.com/scipy/scipy/blob/v0.18.1/scipy/fftpack/basic.py" rel="nofollow"><code>scipy.fftpack.fftn()</code></a> calls <a href="https://github.com/scipy/scipy/blob/master/scipy/fftpack/src/zfftnd.c" rel="nofollow"><code>scipy.fftpack.zfftnd()</code></a> which calls <a href="https://github.com/scipy/scipy/blob/master/scipy/fftpack/src/zfft.c" rel="nofollow"><code>zfft()</code></a>, based on <a href="https://github.com/scipy/scipy/blob/master/scipy/fftpack/src/dfftpack/zfftf1.f" rel="nofollow"><code>zfftf1</code></a> which does not seem to handle strided DFTs. Moreover, <code>zfftnd()</code> calls many times the function <a href="https://github.com/scipy/scipy/blob/master/scipy/fftpack/src/zfftnd.c" rel="nofollow"><code>flatten()</code></a> which performs the transposition.</p></li>
<li><p>Another example: for parallel distributed memory multidimensionnal DFTs, <a href="http://www.fftw.org/fftw3_doc/Transposed-distributions.html#Transposed-distributions" rel="nofollow">FFTW-MPI uses the first strategy</a> to avoid any MPI communications between processes during 1D DTFs. Of course, <a href="http://www.fftw.org/fftw3_doc/FFTW-MPI-Transposes.html#FFTW-MPI-Transposes" rel="nofollow">functions to transpose the array</a> are not far away and a lot a MPI communications are involved in the process.</p></li>
</ul>

<blockquote>
  <p>Are there any other numpy array routines that could affect stride structure? What can be done in those cases?</p>
</blockquote>

<p>You can <a href="https://github.com/numpy/numpy/search?p=1&amp;q=swapaxes&amp;utf8=%E2%9C%93" rel="nofollow">search the github repository of numpy for <code>swapaxes</code></a>: this funtion is only used a couple of times. Hence, to my mind, this "change of strides" is particular to <code>fft.fftn()</code> and most numpy functions keep the strides unchanged.</p>

<p>Finally, the "change of strides" is a feature of the first strategy and there is no way to prevent that. The only workaround is to swap the axes at the end of the computation, which is costly. But you can rely on <code>pyfftw</code> since <code>fftw</code> implements the second strategy in a very efficient way. The DFT computations will be faster, and subsequent computations will be faster as well if the strides of the different arrays become consistent.</p>

python arrays numpy memory-management fft