title: Convolution of NumPy arrays of arbitrary dimension for Cauchy product of multivariate power series
body: <p>I'm trying to implement the idea I have suggested <a href="https://math.stackexchange.com/questions/2877478/cauchy-product-of-multivariate-formal-power-series">here</a>, for <a href="https://en.wikipedia.org/wiki/Cauchy_product" rel="nofollow noreferrer">Cauchy product</a> of multivariate finite power series (i.e. polynomials) represented as NumPy ndarrays. <code>numpy.convolve</code> does the job for 1D arrays, respectively. But to my best knowledge there is no implementations of convolution for arbitrary dimensional arrays. In the above link, I have suggested the equation:</p>
  
  <p><a href="https://i.stack.imgur.com/agkOk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/agkOk.png" alt="enter image description here"></a></p>
  
  <p>for convolution of two <code>n</code> dimensional arrays <code>Phi</code> of shape <code>P=[p1,...,pn]</code> and <code>Psi</code> of the shape <code>Q=[q1,...,qn]</code>, where:</p>
  
  <ol>
  <li><code>omega</code>s are the elements of <code>n</code> dimensional array <code>Omega</code> of the shape <code>O=P+Q-1</code></li>
  <li><code>&lt;A,B&gt;_F</code> is the generalization of <a href="https://en.wikipedia.org/wiki/Frobenius_inner_product" rel="nofollow noreferrer">Frobenius inner product</a> for arbitrary dimensional arrays <code>A</code> and <code>B</code> of the same shape</li>
  <li><code>A^F</code> is <code>A</code> flipped in all <code>n</code> directions</li>
  <li><code>{A}_[k1,...,kn]</code> is a slice of <code>A</code> starting from <code>[0,...,0]</code> to <code>[k1,...,kn]</code></li>
  <li><code>Psi'</code> is <code>Psi</code> extended with zeros to have the shape <code>O</code> as defined above</li>
  </ol>
  
  <p>I tried implementing the above functions one by one: </p>
  
  <pre><code>import numpy as np
  </code></pre>
  
  
  
  <pre><code>def crop(A,D1,D2):
      return A[tuple(slice(D1[i], D2[i]) for i in range(D1.shape[0]))]
  </code></pre>
  
  <p>as was suggested <a href="https://stackoverflow.com/a/51788969/4999991">here</a>, slices/crops <code>A</code> from <code>D1</code> to <code>D2</code>,</p>
  
  <pre><code>def sumall(A):
      sum1=A
      for k in range(A.ndim):
          sum1 = np.sum(sum1,axis=0)
      return sum1 
  </code></pre>
  
  <p>is a generalization of <code>numpy.sum</code> for multidimensional ndarrays,</p>
  
  <pre><code>def flipall(A):
      A1=A
      for k in range(A.ndim):
          A1=np.flip(A1,k)
      return A1
  </code></pre>
  
  <p>flips <code>A</code> is all existing axises, and finally </p>
  
  <pre><code>def conv(A,B,K):
      D0=np.zeros(K.shape,dtype=K.dtype)
      return sumall(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,K-B.shape)) \
                             ,np.minimum(A.shape,K)), \
                        flipall(crop(B,np.maximum(D0,np.minimum(B.shape,K-A.shape)) \
                             ,np.minimum(B.shape,K)))))
  </code></pre>
  
  <p>where <code>K=[k1,...,kn]</code> and for all <code>0&lt;=kj&lt;=oj</code>, is a modified version of formula above which only calculate the non-zero multiplications to be more efficient. Now I'm trying to populate the <code>Omega</code> array using <code>fromfunction</code> or <code>meshgrid</code> in combination to  <code>vectorize</code> as suggested <a href="https://stackoverflow.com/questions/51782224/populating-a-numpy-ndarray-of-an-arbitrary-shape-preferably-without-using-for-l">here</a>, but I have failed so far. Now my questions in prioritized order are: </p>
  
  <ul>
  <li>how can I implement the final step and populate the final array in an efficient and pythonic way?</li>
  <li>are there more efficient implementations of the functions above? Or how would you implement the formula?</li>
  <li>is my equation correct? does this represent multiplication of multivariate finite power series?</li>
  <li>haven't really others implemented this before in NumPy or am I reinventing the wheel here? I would appreciate if you could point me towards other solutions. </li>
  </ul>
  
  <p>I would appreciate if you could help me with these questions. Thanks for your help in advance. </p>
  
  <p><strong>P.S.1</strong> You may find some examples and other information in this <a href="https://gist.github.com/Foadsf/80abba07ea892e50da3df7e318844f33" rel="nofollow noreferrer">GitHub Gist</a></p>
  
  <p><strong>P.S.2</strong> <a href="https://mail.python.org/pipermail/astropy/2018-August/004440.html" rel="nofollow noreferrer">Here in the AstroPy mailing list</a> I was told that <code>scipy.signal.convolve</code> and/or <code>scipy.ndimage.convolve</code> do the job for higher dimensions as well. There is also a <code>scipy.ndimage.filters.convolve</code>. <a href="https://mail.python.org/pipermail/astropy/2018-August/004442.html" rel="nofollow noreferrer">Here</a> I have explained why they are not what I'm looking for. </p>