# Polynomial Regression

## Setup


```python
%matplotlib inline
```


```python
import numpy as np
import matplotlib.pyplot as plt
```

Define the convenience function *polynomial_features* that implements the polynomial feature mapping/transformation $\boldsymbol{\phi}(\cdot)$ for a $D$-degree polynomial

\begin{align}
    \boldsymbol{\phi}(x) & \triangleq
    \begin{bmatrix}
        x^D     \\
        x^{D-1} \\
        \vdots  \\
        x^2     \\
        x       \\
        1
    \end{bmatrix}    
\end{align}

For a vector $\mathbf{x} \in \mathbb{R}^N$ as input argument, it constructs and returns the matrix

\begin{align}
    \boldsymbol{\Phi} & \triangleq
    \begin{bmatrix}
        \boldsymbol{\phi}^T(x_1) \\
        \boldsymbol{\phi}^T(x_2) \\
        \vdots  \\
        \boldsymbol{\phi}^T(x_N) \\ 
    \end{bmatrix} =    
    \begin{bmatrix}
        x_1^D & x_1^{D-1} & \cdots & x_1^2 & x_1 & 1 \\
        x_2^D & x_2^{D-1} & \cdots & x_2^2 & x_2 & 1 \\
        \vdots  \\
        x_N^D & x_N^{D-1} & \cdots & x_N^2 & x_N & 1 
    \end{bmatrix} \in \mathbb{R}^{N \times (D+1)}   
\end{align}




```python
def polynomial_features(x, degree):
    """
    Construct design matrix for polynomial regression of degree `polydegree`
    
    Inputs
    ------
    x: numpy.ndarray, ndims=1
        A 1-dimensional numpy array of length `N` contianing `x` values
    
    degree: int
        An integer specifying the maximum degree of polynomial to be constructed
        in the design matrix

    Returns
    -------
    Phi: np.ndarray, ndims=2
        A 2-dimensional numpy array of shape `(N, polydegree+1)` representing the 
        design matrix. Columns are: 
        `[x**polydegree, x**(polydegree-1) ... x**2, x**1, 1]`
    """
    N = x.shape[0]
    Phi = np.ones((N, degree + 1))
    for deg in range(1, degree+1):
        Phi[:, -deg-1] = Phi[:, -deg] * x
    
    return Phi
```

## Specification of true model

We want the true model to be the 3rd degree polynomial in $x$

\begin{align}
    & f(x|\mathbf{w}_{\mathrm{true}}) = \frac{A}{3} x^3 - \frac{A(a+b)}{2} x^2 + A a b x = \boldsymbol{\phi}^T(x)  \mathbf{w}_{\mathrm{true}}
\end{align}

where the polynomial coefficients $\mathbf{w}_{\mathrm{true}}$ are given as

\begin{align}
    \mathbf{w}_{\mathrm{true}} & \triangleq
    \begin{bmatrix}
        \frac{A}{3} \\
        -\frac{A (a+b)}{2} \\
        a b A \\
        0
    \end{bmatrix}
\end{align}

and the polynomial feature map $\boldsymbol{\phi}(\cdot)$ is given as

\begin{align}
    \boldsymbol{\phi}(x) & \triangleq
    \begin{bmatrix}
        x^3 \\
        x^2 \\
        x   \\
        1
    \end{bmatrix}    
\end{align}


```python
## True model (polynomial) specification

# The parameters that follow influence the shape of the polynomial (feel free to change)
a = 1.0    # First stationary point of polynomial
b = 1.5    # Second stationary point of polynomial
A = 2.0    # Scaling factor

# Polynomial coefficients (do not change)
w_true = np.zeros(4)
w_true[0] = A / 3.0;
w_true[1] = - A * (a + b) / 2.0
w_true[2] = A * a * b
w_true[3] = 0.2

w_true
```

## True Data Generation Process

We will feed the polynomial regression models noisy data obtained from the model above

For any $x$ value $x_i$ we will compute the corresponding target $y_i$ as:

\begin{align*}
y_i &= \phi(x_i) \cdot w_i  + \epsilon_i, \quad \epsilon_i \stackrel{\text{iid}}{\sim} N(0, \sigma)
\end{align*}

The parameter $\sigma$ controls the strength of the noise

Below we sample oa grid of x on $[0, 3]$, compute $\phi(x)$, then add random noise to compute $y$


```python
np.random.seed(42)  # for reproducible resultss

N = 100
xmin, xmax = 0, 3
x = xmin + (xmax - xmin) * np.random.rand(N)


sigma = 0.05
epsilon = np.random.randn(100) * np.sqrt(sigma)

## Construct the target values
Phi = polynomial_features(x, 3)
y = Phi @ w_true + epsilon

# also construct data for plotting later
x_plot = np.linspace(xmin, xmax, 40)
y_true_plot = polynomial_features(x_plot, 3) @ w_true
```

## Data Partitioning

The available data will be partitioned into a *training* set and a *hold-out* set (i.e. samples not used for training/fitting). In particular, the hold-out set will be used to obtain an honest estimate of how good the model performs on previously unseen samples (i.e. samples other than ones used for training).


```python
## Partition data into training and holdout set

# Number of training samples (feel free to change)
N_train = 10

N_holdout = N - N_train

x_train = x[:N_train]
y_train = y[:N_train]

x_holdout = x[N_train:]
y_holdout = y[N_train:]
```

## Model Training/Fitting


```python
# Specify the degree of the polynomial that we plan to fit (feel free to change)
polydegree = 3   # must be integer >=0 

# Construct design matrix for polynomial regression (with intercept term)
# for the training data
Phi_train = polynomial_features(x_train, polydegree)

# Fit/train the model by computing the optimal coefficient (weight) vector
Phi_train_dagger = np.linalg.pinv(Phi_train)
w_star = np.dot(Phi_train_dagger, y_train)

# Compute the (optimal) predicted outputs for the training inputs
y_hat_train = np.dot(Phi_train, w_star)

# Compute the (minimum) MSE on the training set
res_train = y_train - y_hat_train   # these are the optimal regression residuals
MSE_train = np.linalg.norm(res_train,2)**2 / N_train
```

## Hold-Out Set Predictions


```python
# Construct design matrix for polynomial regression (with intercept term)
# for the hold-out data
Phi_holdout = polynomial_features(x_holdout, polydegree)

# Compute the (optimal) predicted outputs for the hould-out inputs
# using the optimal coefficients obtained by fitting the model
y_hat_holdout = np.dot(Phi_holdout, w_star)

# Compute the (minimum) MSE on the hold-out set
res_holdout = y_holdout - y_hat_holdout  
MSE_holdout = np.linalg.norm(res_holdout,2)**2 / N_holdout
```

## Reporting


```python
print('The optimal coefficients are:\n')
print(w_star)

print('\n')

print('The (minimum) MSE of the fitted model as computed on the training set is\nMSE_train={:.3f}\n'.format(MSE_train))
print('The MSE of the fitted model as computed on the hold-out set is\nMSE_holdout={:.3f}\n'.format(MSE_holdout))
print('The noise variance of the true model is\nsigmaSquared_e={:.3f}\n'.format(sigma))
```

## Graphing of Results


```python
## Create data for plotting the fitted polynomial curve
Phi_plot = polynomial_features(x_plot, polydegree)
y_hat_plot_fitted = np.dot(Phi_plot, w_star)

```


```python
# Resize figure to make it bigger
fig, ax = plt.subplots(figsize=(13, 8))

ax.plot(x_plot, y_hat_plot_fitted, 'r', lw=3, label="Fitted model")
ax.plot(x_plot, y_true_plot, 'g', lw=3, label="True model")
ax.scatter(x_holdout, y_holdout, s=100, facecolors='none', edgecolors='b', label="Hold-Out data")
ax.scatter(x_train, y_train, s=100, facecolors='r', edgecolors='k', label="train data")

# Annotate graph
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='lower right')
ax.set_title((
    "Fitted vs. True Model Responses: "
    f"Ntrain={N_train:d} deg={polydegree:d} MSEtrain={MSE_train:.3f}, "
    f"MSEholdout={MSE_holdout:.3f}, NoiseVar={sigma:.3f}"
))
```

## Using sklearn

Below we will show an example of how to use scikit-learn to:

1. construct the polynomial feature matrix
2. Fit the polynomial regression
3. Compute the MSE on training and holdout sets

We will wrap this process up in a function that takes as an input the degree of the polynomial as well as the number of points to use for training

This function will then be used to construct an interactive GUI where we will get sliders for the degree and number of training points, and as we adjust the sliders a plot of the model fit updates automatically

Feel free to experiment with the two sliders and try to gain some intutions for when overfitting will and won't occur


```python
from sklearn import preprocessing, linear_model, pipeline, metrics, model_selection
from ipywidgets import widgets
```


```python
def polyreg_demo(degree=8, n_train=10):
    # define model
    model = pipeline.make_pipeline(
        preprocessing.PolynomialFeatures(degree=degree),
        linear_model.LinearRegression()
    )
    
    X = x[:, None]  # convert to 2d
    test_size = X.shape[0] - n_train
    split = model_selection.train_test_split(X, y, train_size=n_train, test_size=test_size, random_state=12)
    X_train, X_test, y_train, y_test = split
   
    # fit model
    model.fit(X_train, y_train)
    yhat = model.predict(x_plot[:, None])
    
    # compute metrics
    MSE_train = metrics.mean_squared_error(y_train, model.predict(X_train))
    MSE_holdout = metrics.mean_squared_error(y_test, model.predict(X_test))

    # make the plot
    fig, ax = plt.subplots(figsize=(11, 8))
    ax.plot(x_plot, yhat, "k-", lw=3, label="Fitted Model")
    ax.scatter(x_holdout, y_holdout, color="b", s=60, alpha=0.5, label="Hold-Out Data")
    ax.scatter(X_train.flatten(), y_train, color="r", s=80, alpha=0.7, label="Training Data") 
    ax.set_ylim((-20, 20))
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.legend(loc='upper left')
    ax.set_title((
        "Fitted vs. True Model Responses: "
        f"Ntrain={N_train:d} deg={degree:d} MSEtrain={MSE_train:.3f}, "
        f"MSEholdout={MSE_holdout:.3f}, NoiseVar={sigma:.3f}"
    ))
    
    return fig


widgets.interactive(polyreg_demo, degree=(1, 10, 1), n_train=(5, 95, 1))
```


```python
fig = polyreg_demo(0, 10)
fig.set_size_inches((11, 5))
ax = fig.axes[0]
ax.set_ylim((-2, 5))
fig.savefig("underfit.pdf")
fig.savefig("../overfitting_model_selection/underfitting.pdf")
```

## Using Keras (Optional)

We have had a few questions about how Linear Regression relates to a neural network in keras

In this section we will show how you can do Linear Regression using keras

Note that this is like using a blowtorch when a match would do, but it will illustrate the following:

- How you can view a LinearRegression as a `Dense` neural network layer

- This will be the first example of an _iterative_ approach to finding coefficients for a linear regression. As mentioned in the lectures this the most numerically stable approach and often the only computationaly reasonable approach as either the number of samples or the number of features gets very large

### Theory

Recall that a multi-layer-perceptron (MLP) is specific kind of neural network that is constructed by stacking one or more `Dense` layers, separated by non-linear activation functions

An example of a 2 hidden layer MLP in keras is:

```python
from tensorflow import keras

model = keras.Sequential()
model.add(keras.layers.Dense(10), activation="relu")  # hidden layer 1
model.add(keras.layers.Dense(5), activation="relu")   # hidden layer 2
model.add(keras.layers.Dense(1), activation=None)     # output layer -- assuming we have 1 target to predict
```

Each dense layer in the network has the following structure: 

$$\text{output} = \mathbf{input} \mathbf{W} + \mathbf{b}$$ 

where $\mathbf{input}$ is the input to the layer, $\mathbf{W}$ is a matrix of weights or coefficients and $\mathbf{b}$ is a vector of intercepts

The number of columns in the matrix $\mathbf{W}$ and the number of elements in the vector $\mathbf{b}$ are the number of neurons in the hidden layer

In our example, the $\mathbf{W}$ for the first hidden layer would have shape $n_{\text{features}} \times 10$

The output of this first hidden layer would he $n_{\text{samples}} \times 10$

The $\mathbf{W}$ for the second hidden layer would have shape $10 \times 5$ -- there are $5$ rows because that is how many neurons there are in the hidden layer

The $\mathbf{b}$ fot the second layer would have $5$ elements and the output would be of shape $n_{\text{samples}} \times 5$, and so on until the last layer...

Given an input sample $\mathbf{x}$, we can compute the prediction of the network as follows:

$$\begin{align*}
\text{output}_1 &= \text{relu}(\mathbf{x} \mathbf{W}_1 + \mathbf{b}_1) \\
\text{output}_2 &= \text{relu}(\text{output}_1 \mathbf{W}_2 + \mathbf{b}_2) \\
\hat{y} &= \text{output}_2 \mathbf{W}_{\text{final}} + \mathbf{b}_{\text{final}} \\
        &= (\text{relu}(\text{relu}(\mathbf{x} \mathbf{W}_1 + \mathbf{b}_1)\mathbf{W}_2 + \mathbf{b}_2) \mathbf{W}_{\text{out}} + \mathbf{b}_{\text{out}}
\end{align*}$$

### Special Case: Linear Regression

We can consider a special case where there are no hidden layers and no activation function

The keras model will be written as:

```python
model = keras.Sequential()
model.add(keras.layers.Dense(1), activation=None)   # output layer -- assuming we have 1 target to predict
```

and the prediction for an input sample $\mathbf{x}$ is

$$\hat{y} = \mathbf{x}\mathbf{W} + \mathbf{b},$$

which is exactly our linear regression equation!

Let's test this out with our polynomial regression


```python
from tensorflow import keras

model = keras.Sequential()
model.add(keras.layers.Dense(1, use_bias=False, activation=None))   # output layer -- assuming we have 1 target to predict

model.compile(optimizer="rmsprop", loss="mse")  # note we use mse as loss!
```

Now let's fit a degree 3 polynomial regression for 10 epochs


```python
X3 = polynomial_features(x_train, 3)
X3_holdout = polynomial_features(x_holdout, 3)
```


```python
model.fit(X3, y_train, epochs=10, validation_data=(X3_holdout, y_holdout));
```

Notice that keras is reporting the `loss` after each pass through the training data

This is the MSE with the current parameters

It is changing because the parameters are being updated (stay tuned for later in the course and you'll learn _how_ theose parameters are being updated)

Also, notice that the MSE after 10 epocs is much larger than the optimal MSE of 0.015 we found using the closed form optimal solution from linear regression 

Let's train the model for many more epochs and see if we can close that gap


```python
history = model.fit(X3, y_train, epochs=4000, verbose=False, initial_epoch=11, validation_data=(X3_holdout, y_holdout))
```


```python
fig, ax = plt.subplots(figsize=(8, 6))
x = np.arange(len(history.history["loss"]))
ax.plot(x, history.history["loss"], label="train loss")
ax.plot(x, history.history["val_loss"], label="Validation loss")
ax.hlines(MSE_train, xmin=x[0], xmax=x[-1], alpha=0.5, lw=4, color="red")
ax.set_ylabel("MSE")
ax.set_xlabel("Epoch")
ax.legend()

msg = "The final MSE_train and MSE_validation are {} and {}"
print(msg.format(history.history["loss"][-1], history.history["val_loss"][-1]))
```

Now the training MSE is much closer to the optimum than before

b75d8864f77e8e86b6a589b15e874e1daaec0b18
ipynb
Jupyter Notebook
2020-10-02_validation/polynomial-regression.ipynb
sglyon/UCF-MSDA-ML-Seminar
0f62fa2f3ee190e70c83a9dc52a4749bb4ce7d60
MIT
2020-10-02_validation/polynomial-regression.ipynb
sglyon/UCF-MSDA-ML-Seminar
0f62fa2f3ee190e70c83a9dc52a4749bb4ce7d60
MIT
2020-10-02_validation/polynomial-regression.ipynb
sglyon/UCF-MSDA-ML-Seminar
0f62fa2f3ee190e70c83a9dc52a4749bb4ce7d60
MIT
Qwen/Qwen-72B
1. YES
2. YES