# Ridge regression - varying alpha and observing the residual

I am trying to reproduce this figure from Bishop:

Residual vs. Alpha (lambda in figure)

The code is pasted below:

import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model

from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error

x_train = np.linspace(0, 1, 10)
noise_10 = np.random.normal(0, 0.3, 10)
y_train = np.sin(2*np.pi*x_train) + noise_10
x_train = x_train[:, np.newaxis]

x_test = np.linspace(0, 1, 100)
noise_100 = np.random.normal(0, 0.3, 100)
y_test = np.sin(2*np.pi*x_test) + noise_100
x_test = x_test[:, np.newaxis]

n_alphas = 200
alphas = np.logspace(-40, -18, n_alphas)

errors = []
for a in alphas:
ridge = make_pipeline(PolynomialFeatures(degree = 9),
Ridge(alpha=a))
ridge.fit(x_train, y_train)
mse = mean_squared_error(y_test, ridge.predict(x_test))
errors.append(np.sqrt(mse))

print(errors)


However, the errors array has the same value for all values of alpha. It's taking the first value of alpha = np.exp(-40) and all the other values seem to be the same for all future iterations of the for loop. How can I correct this error?

I still haven't figured out what the previously posted code is doing wrong. However, manually populating the alphas array gives me results that are close to the original figure.

import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error
np.random.seed(12344)
x_1000 = np.linspace(0, 1, 1000)
x_train = np.linspace(0, 1, 10)
noise_10 = np.random.normal(0, 0.3, 10)
y_train = np.sin(2*np.pi*x_train) + noise_10
x_train = x_train[:, np.newaxis]
x_test = np.linspace(0, 1, 100)
noise_100 = np.random.normal(0, 0.3, 100)
y_test = np.sin(2*np.pi*x_test) + noise_100
x_test = x_test[:, np.newaxis]
low_alpha = make_pipeline(PolynomialFeatures(degree = 9), Ridge(alpha=np.exp(-18)))
low_alpha.fit(x_test, y_test)
plt.figure(1)
plt.plot(x_train, low_alpha.predict(x_train), label = 'alpha = ln(-18)')
plt.plot(x_1000, np.sin(2*np.pi*x_1000))
plt.legend()
plt.show
high_alpha = make_pipeline(PolynomialFeatures(degree = 9),
Ridge(alpha=np.exp(0)))
high_alpha.fit(x_test, y_test)
plt.figure(2)
plt.plot(x_train, high_alpha.predict(x_train), label = 'alpha = 1')
plt.plot(x_1000, np.sin(2*np.pi*x_1000))
plt.legend()
plt.show
alphas = np.array([np.exp(-30), np.exp(-29), np.exp(-28), np.exp(-27), np.exp(-26), np.exp(-25), np.exp(-24), np.exp(-23), np.exp(-22), np.exp(-21), np.exp(-20), np.exp(-19), np.exp(-18), np.exp(-17), np.exp(-16), np.exp(-15), np.exp(-14), np.exp(-13), np.exp(-12), np.exp(-11), np.exp(-10), np.exp(-9), np.exp(-8), np.exp(-7), np.exp(-6), np.exp(-5), np.exp(-4), np.exp(-3), np.exp(-2), np.exp(-1), np.exp(-1), np.exp(0), np.exp(1)])

test_errors = []
train_errors = []

for a in np.nditer(alphas):
ridge = make_pipeline(PolynomialFeatures(degree = 9),
Ridge(alpha=a))
ridge.fit(x_train, y_train)
mse_train = mean_squared_error(y_train, ridge.predict(x_train))
mse_test = mean_squared_error(y_test, ridge.predict(x_test))
train_errors.append(np.sqrt(mse_train))
test_errors.append(np.sqrt(mse_test))
plt.figure(3)
plt.plot(alphas, test_errors, 'g^', label = 'Test Error')
plt.plot(alphas, train_errors, 'bs', label = 'Train Error')
plt.xscale('log')
plt.xlabel('Regression coefficient Lambda')
plt.ylabel('Residuals')
plt.legend()
plt.show()


The output (only for the third figure) is shown below:

I have taken a look on your code. You obtain same errors results for each alpha value because your regularization strength is too small. Replacing :

alphas = np.logspace(-40, -18, n_alphas)


with :

alphas = np.logspace(-40, -1, n_alphas)


will yields different errors values for alpha values large enough. Are you sure about figure alpha values? Do you have a link to this hands-on?

Also, I would like to highlight the fact that you have to standardize your features before using regularization. Reason is, by creating polynomial features, polynomial features will have different magnitudes. Therefore when fitting model, coefficients to be estimated won't have the same magnitudes neither and so regularization will highly penalize coefficients with large values. Standardization / Normalization is a strong prerequisite to regularization.