岭回归算法详解（Python实现）

岭回归与多项式回归唯一的不同在于代价函数上的差别。岭回归的代价函数如下：


为了方便计算导数，通常也写成下面的形式：


其中，w 是长度为 n 的向量，不包括截距项的系数 θ0；θ 是长度为 n+1 的向量，包括截距项的系数 θ0；m 为样本数；n 为特征数。

岭回归的代价函数仍然是一个凸函数，因此可以利用梯度等于 0 的方式求得全局最优解（正规方程）：

θ=(XTX+λI)-1(XTy)

上述正规方程与一般线性回归的正规方程相比，多了一项 λI，其中 I 表示单位矩阵。假如 X^TX 是一个奇异矩阵（不满秩），添加这一项后可以保证该项可逆。由于单位矩阵的形状是对角线上为 1，其他地方都为 0，看起来像一条山岭，因此而得名。

除了上述正规方程之外，还可以使用梯度下降的方式求解：


为了计算方便，在式中添加 θ0=0 到 w。因此在梯度下降的过程中，参数的更新可以表示为：


其中，α 为学习率，λ 为正则化项的参数。

岭回归的实现

'''数据以及相关函数'''
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

data = np.array([[-2.95507616, 10.94533252],
                 [-0.44226119, 2.96705822],
                 [-2.13294087, 6.57336839],
                 [1.84990823, 5.44244467],
                 [0.35139795, 2.83533936],
                 [-1.77443098, 5.6800407],
                 [-1.8657203, 6.34470814],
                 [1.61526823, 4.77833358],
                 [-2.38043687, 8.51887713],
                 [-1.40513866, 4.18262786]])
m = data.shape[0]                          # 样本大小
X = data[:, 0].reshape(-1, 1)             # 将array转换为矩阵
y = data[:, 1].reshape(-1, 1)

'''岭回归的实现'''
# 代价函数
def L_theta(theta, X_x0, y, lamb):
    h = np.dot(X_x0, theta)                # np.dot表示矩阵乘法
    theta_without_t0 = theta[1:]
    L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(theta_without_t0))
    return L_theta

# 梯度下降
def GD(lamb, X_x0, theta, y, alpha):
    for i in range(T):
        h = np.dot(X_x0, theta)
        theta_with_t0_0 = np.r_[np.zeros([1, 1]), theta[1:]]   # 设theta[0]=0
        theta -= (alpha * 1/m * np.dot(X_x0.T, h - y) + lamb * (theta_with_t0_0))
        # 添加正则化项的梯度
        if i % 50000 == 0:
            print(L_theta(theta, X_x0, y, lamb))
    return theta

T = 1200000                                # 迭代次数
degree = 11
theta = np.ones((degree + 1, 1))          # 参数的初始化，degree=11
alpha = 0.0000000006                      # 学习率
lamb = 0.0001
poly_features_d = PolynomialFeatures(degree=degree, include_bias=False)
X_poly_d = poly_features_d.fit_transform(X)
X_x0 = np.c_[np.ones((m, 1)), X_poly_d]   # 向每个实例添加X0=1
theta = GD(lamb=lamb, X_x0=X_x0, theta=theta, y=y, alpha=alpha)

由于自由度比较大，此时利用梯度下降的方法训练模型较困难，学习率稍微大些就会出现损失函数的值越过最低点时不断增长的情况。下面是训练结束后的参数以及代价函数值：

185842996.9976393
3.3865692058065915
3.6003632848879725
3.602155809370964
...
3.599879516788523
3.5997451046573166
3.599610698056696
3.599476296986373

从结果看，截距项的参数最大，高阶项的参数都比较小。以下代码比较原始数据和训练出来的模型之间的关系：

X_plot=np.linspace(-2.99, 1.9, 1000).reshape(-1, 1)
poly_features_d_with_bias=PolynomialFeatures(degree=degree, include_bias=True)
X_plot_poly=poly_features_d_with_bias.fit_transform(X_plot)
y_plot=np.dot(X_plot_poly, theta)
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

效果如下图所示：


图 5 显示模型与原始数据的匹配度不是太好，但是过拟合的情况极大地改善了，模型变得更简单了。

正则方程

下面代码使用正则方程求解，其中 λ=0。

theta2=np.linalg.inv(np.dot(X_x0.T, X_x0)+  10 * np.identity(X_x0.shape[1])).dot(X_x0.T)
.dot(y)
print(theta2)
print(L_theta(theta2, X_x0, y, lamb))

运行程序，参数即代价函数的值为：

[[0.56502654]
[-0.12459547]
[0.26772443]
[-0.15642405]
[0.29249514]
[-0.10084392]
[0.22791769]
[0.1648667]
[-0.05686718]
[-0.03906615]
[-0.00111673]
[0.00101724]]
0.604428719059618

从参数来看，截距项的系数减小了，1~7 阶都有比较大的参数，后面更高阶项的参数越来越小：

X_plot=np.linspace(-3, 2, 1000).reshape(-1, 1)
poly_features_d_with_bias=PolynomialFeatures(degree=degree, include_bias=True)
X_plot_poly=poly_features_d_with_bias.fit_transform(X_plot)
y_plot=np.dot(X_plot_poly, theta2)
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.show()　　　#效果如图 6 所示

下图为函数图像：


从图 6 中可以看到，虽然模型的自由度没变，还是 11，但是过拟合的程度得到了改善。