title: 2d numpy.power for polynomial expansion
body: <p>I am trying to write a function that maps 2d-ndarray to 2d-ndarray. The rows of the input array can be processed independently and there shall be a 1-to-1 correspondence between rows of the input and rows of the output. For each row of the input, the polynomial expansion of a given order for the row shall be computed (see docstring for an example). The current implementation works; however it requires an explicit loop over the rows and duplication of rows in the "powerMatrix"). Is it possible to get the same result with a single call to numpy.power? Btw.: the order of the entries in the result's rows doesn't matter to me.</p>
  
  <pre><code>import numpy
  def polynomialFeatures(x, order):
      """ Generate polynomial features of given order for data x.
  
      For each row of ndarray x, the polynomial expansions are computed, i.e
      for row [x1, x2] and order 2, the following row of the result matrix is
      computed: [1, x1, x1**2, x2, x1*x2, x1**2*x2, x2**2, x1*x2**2, x1**2*x2**2]
  
      Parameters
      ----------
      x : array-like
          2-D array; for each of its rows, the polynomial features are created
  
      order : int
          The order of the polynomial features
  
      Returns
      -------
      out : ndarray
          2-D array of shape (x.shape[0], (order+1)**x.shape[1]) containing the 
          polynomial features computed for the rows of the array x
  
      Examples
      --------
      &gt;&gt;&gt; polynomialFeatures([[1, 2, 3], [-1, -2, -3]], 2)
      array([[  1   3   9   2   6  18   4  12  36   1   3   9   2   6  18   4  12  
               36   1   3   9   2   6  18   4  12  36]
             [  1  -3   9  -2   6 -18   4 -12  36  -1   3  -9   2  -6  18  -4  12 
              -36   1  -3   9  -2   6 -18   4 -12  36]])
      """
      x = numpy.asarray(x)
      # TODO: Avoid duplication of rows
      powerMatrix = numpy.array([range(order+1)] * x.shape[1]).T
      # TODO: Avoid explicit loop, and use numpy's broadcasting
      F = []
      for i in range(x.shape[0]):
          X = numpy.power(x[i], powerMatrix).T
          F.append(numpy.multiply.reduce(cartesian(X), axis=1))
  
      return numpy.array(F)
  
  print numpy.all(polynomialFeatures([[1, 2, 3], [-1, -2, -3]], 2) ==
                  numpy.array([[1,   3,   9,   2,   6,  18,   4,  12,  36,   1,
                                3,   9,   2,   6,  18,   4,  12,  36,   1,   3, 
                                9,   2,   6,  18,   4,  12,  36],
                               [1,  -3,   9,  -2,   6, -18,   4, -12,  36,  -1,
                                3,  -9,   2,  -6,  18,  -4,  12, -36,   1,  -3,
                                9,  -2,   6, -18,   4, -12,  36]]))
  </code></pre>
  
  <p>Thanks,
  Jan</p>
  
  <p>EDIT: The missing function cartesian is defined here: <a href="https://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays">Using numpy to build an array of all combinations of two arrays</a></p>