title: What is a vectorized way to create multiple powers of a NumPy array?
body: <p>I have a NumPy array:</p>
  
  <pre><code>arr = [[1, 2],
         [3, 4]]
  </code></pre>
  
  <p>I want to create a new array that contains powers of <code>arr</code> up to a power <code>order</code>:</p>
  
  <pre><code># arr_new = [arr^0, arr^1, arr^2, arr^3,...arr^order]
  arr_new = [[1, 1, 1, 2, 1, 4, 1, 8],
             [1, 1, 3, 4, 9, 16, 27, 64]]
  </code></pre>
  
  <p>My current approach uses <code>for</code> loops:</p>
  
  <pre><code># Pre-allocate an array for powers
  arr = np.array([[1, 2],[3,4]])
  order = 3
  rows, cols = arr.shape
  arr_new = np.zeros((rows, (order+1) * cols))
  
  # Iterate over each exponent
  for i in range(order + 1):
      arr_new[:, (i * cols) : (i + 1) * cols] = arr**i
  print(arr_new)
  </code></pre>
  
  <p>Is there a faster (i.e. vectorized) approach to creating powers of an array?</p>
  
  <hr>
  
  <h3>Benchmarking</h3>
  
  <p>Thanks to @hpaulj and @Divakar and @Paul Panzer for the answers. I benchmarked the loop-based and broadcasting-based operations on the following test arrays.</p>
  
  <pre><code>arr = np.array([[1, 2],
                  [3,4]])
  order = 3
  
  arrLarge = np.random.randint(0, 10, (100, 100))  # 100 x 100 array
  orderLarge = 10
  </code></pre>
  
  <p>The <code>loop_based</code> function is:</p>
  
  <pre><code>def loop_based(arr, order):
      # pre-allocate an array for powers
      rows, cols = arr.shape
      arr_new = np.zeros((rows, (order+1) * cols))
      # iterate over each exponent
      for i in range(order + 1):
          arr_new[:, (i * cols) : (i + 1) * cols] = arr**i
      return arr_new
  </code></pre>
  
  <p>The <code>broadcast_based</code> function using <code>hstack</code> is:</p>
  
  <pre><code>def broadcast_based_hstack(arr, order):
      # Create a 3D exponent array for a 2D input array to force broadcasting
      powers = np.arange(order + 1)[:, None, None]
      # Generate values (third axis contains array at various powers)
      exponentiated = arr ** powers
      # Reshape and return array
      return np.hstack(exponentiated)   # &lt;== using hstack function
  </code></pre>
  
  <p>The <code>broadcast_based</code> function using <code>reshape</code> is:</p>
  
  <pre><code>def broadcast_based_reshape(arr, order):
      # Create a 3D exponent array for a 2D input array to force broadcasting
      powers = np.arange(order + 1)[:, None]
      # Generate values (3-rd axis contains array at various powers)
      exponentiated = arr[:, None] ** powers
      # reshape and return array
      return exponentiated.reshape(arr.shape[0], -1)  # &lt;== using reshape function
  </code></pre>
  
  <p>The <code>broadcast_based</code> function using cumulative product <code>cumprod</code> and <code>reshape</code>:</p>
  
  <pre><code>def broadcast_cumprod_reshape(arr, order):
      rows, cols = arr.shape
      # Create 3D empty array where the middle dimension is
      # the array at powers 0 through order
      out = np.empty((rows, order + 1, cols), dtype=arr.dtype)
      out[:, 0, :] = 1   # 0th power is always 1
      a = np.broadcast_to(arr[:, None], (rows, order, cols))
      # Cumulatively multiply arrays so each multiplication produces the next order
      np.cumprod(a, axis=1, out=out[:,1:,:])
      return out.reshape(rows, -1)
  </code></pre>
  
  <p>On Jupyter notebook, I used the <code>timeit</code> command and got these results:</p>
  
  <p><strong>Small arrays (2x2)</strong>:</p>
  
  <pre><code>%timeit -n 100000 loop_based(arr, order)
  7.41 µs ± 174 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
  
  %timeit -n 100000 broadcast_based_hstack(arr, order)
  10.1 µs ± 137 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
  
  %timeit -n 100000 broadcast_based_reshape(arr, order)
  3.31 µs ± 61.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
  
  %timeit -n 100000 broadcast_cumprod_reshape(arr, order)
  11 µs ± 102 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
  </code></pre>
  
  <p><strong>Large arrays (100x100)</strong>:</p>
  
  <pre><code>%timeit -n 1000 loop_based(arrLarge, orderLarge)
  261 µs ± 5.82 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
  
  %timeit -n 1000 broadcast_based_hstack(arrLarge, orderLarge)
  225 µs ± 4.15 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
  
  %timeit -n 1000 broadcast_based_reshape(arrLarge, orderLarge)
  223 µs ± 2.16 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
  
  %timeit -n 1000 broadcast_cumprod_reshape(arrLarge, orderLarge)
  157 µs ± 1.02 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
  </code></pre>
  
  <h3>Conclusions:</h3>
  
  <p>It seems that the broadcast based approach using <code>reshape</code> is faster for smaller arrays. However, for large arrays, the <code>cumprod</code> approach scales better and is faster.</p>