Using Linear Regression to Learn Polynomial Regression
It is quite simple to understand (and to implement using matrices).
Consider a specific example (to generalise later). You have a polynomial function of a single feature $x$):
$$ f(x) = \omega_0 x^0 + \omega_1 x^1 + \ldots \omega_n x^n $$
You can organise coefficients and features in vectors and get $f$ by a scalar product:
$$ \mathbf{\omega} = \begin{pmatrix} \omega_0, \\ \vdots \\ \omega_n \end{pmatrix}, \qquad \mathbf{x} =  \begin{pmatrix} 1, \\ x \\ x^2 \\ \vdots \\ x^n \end{pmatrix}$$
Hence $$ f(x) = \omega^T\mathbf{x}$$.
This is nothing else than a multi-feature linear regression where the $i$-th feature is now the $i$-th power of $x$.
In numpy, imagine you have an array of data x.
To create the vector $\mathbf{x}$ above, you can do (for $n=3$, for instance)
X = np.ones((len(x), 4))
X[:,1] = x
X[:,2] = np.power(x,2)
X[:,3] = np.power(x,3)

And then using sklearn LinearRegression,
model = LinearRegression()
model.fit(X, y)

UPDATE: In sklearn has been recently introduced PolynomialFeatures that precisely performs the transformation I described in numpy (you asked in numpy, but this might be useful as well).