Using Linear Regression to Learn Polynomial Regression

Let's start by considering one-dimensional data, i.e., $$d=1$$. In OLS regression, we would learn the function $$f(x)=w_{0}+w_{1} x,$$ where $$x$$ is the data point and $$\mathbf{w}=\left(w_{0}, w_{1}\right)$$ is the weight vector. To achieve a polynomial fit of degree $$p$$, we will modify the previous expression into $$f(x)=\sum_{j=0}^{p} w_{j} x^{j}$$ where $$p$$ is the degree of the polynomial. We will rewrite this expression using a set of basis functions as $$f(x)=\sum_{j=0}^{p} w_{j} \phi_{j}(x)=\mathbf{w}^{\top} \boldsymbol{\phi}$$ where $$\phi_{j}(x)=x^{j}$$ and $$\phi=\left(\phi_{0}(x), \phi_{1}(x), \ldots, \phi_{p}(x)\right)$$. We simply apply this transformation to every data point $$x_{i}$$ to get a new dataset $$\left\{\left(\phi\left(x_{i}\right), y_{i}\right)\right\}$$. Then we use linear regression on this dataset, to get the weights $$\mathbf{w}$$ and the nonlinear predictor $$f(x)=\sum_{j=0}^{p} w_{j} \phi_{j}(x),$$ which is a polynomial (nonlinear) function in the original observation space. Notes

How does this work? Could anyone give me an example and explain it to me in simple terms? How would I go about implementing this in Numpy?

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).