Assume an n-dimensional array of observations that are reshaped to be a 2d-array with each row being one observation set. Using this reshape approach, `np.polyfit` can compute 2nd order fit coefficients for the entire ndarray (vectorized):

```
fit = np.polynomial.polynomialpolyfit(X, Y, 2)

```

where Y is shape (304000, 21) and X is a vector. This results in a (304000,3) array of coefficients, fit. 

Using an iterator it is possible to call `np.polyval(fit, X)` for each row. This is inefficient when a vectorized approach may exist. Could the `fit` result be applied to the entire observation array without iterating? If so, how?

This is along the lines of [this SO](https://stackoverflow.com/questions/19282429/regression-along-a-dimension-in-a-numpy-array) question.