Using Linear Regression to Learn Polynomial Regression

Question

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?

Answer 1

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