动手学机器学习(第二版)-第四章训练模型

第四章 训练模型

线性模型(Linear Regression)

线性模型的一般形式如下

$$\hat{y}={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_n{x}_n$$

• $\hat{y}$ 预测结果
• $n$ 特征数量
• $x_i$ 第 $i$ 个特征值
• ${\theta }_j$ 是第 $j$ 个模型参数(包括偏置项和特征圈子 ${\theta }_1,{\theta }_2,\dots ,{\theta }_n)$

用矢量形式描述如下

$$\hat{y}=h_{\theta }(X)=\theta \cdot X$$

• $\theta$ 模型的参数矢量，包含偏置项和特征值权重 ${\theta }_1-{\theta }_n$
• $X$ 实例的特征矢量，包含 $x_0$ 到 $x_n$，其中 $x_0$ 始终置为1
• $\theta \times X$ 是两个矢量的乘积，按照线性代数计算展开是 ${\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_n{x}_n$
• $h_{\theta }$是假设函数

正则方程(the Normal Equation)

求解最小化代价函数的 $\theta$ 值，有一个闭合解(closed-form solution)，换句话说，可以通过数学公式直接计算得到结果，这个公式叫做正规方程(Normal Equation)

$$\hat{\theta }=(X^{T}X{)}^{-1}{X}^{T}y$$

• $\hat{\theta }$ 是最小化代价函数的 $\theta$ 值
• $y$ 是包含 $y^{(1)}$ 到 $y^{(n)}$ 的矢量

让我们生成一些线性的数据来测试这个公式

import numpy as np

X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

X_b = np.c_[np.ones((100, 1)), X] # 添加x0 = 1列
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) # 正规方程计算

theta_best

array([[3.97884872],
       [2.92718847]])

X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)
y_predict

array([[3.97884872],
       [9.83322566]])

import matplotlib.pyplot as plt

plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0, 2, 0, 15])
plt.show()

使用Scikit-Learn来实现线性回归就很简单了

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(X, y)
lin_reg.intercept_, lin_reg.coef_

(array([3.97884872]), array([[2.92718847]]))

lin_reg.predict(X_new)

array([[3.97884872],
       [9.83322566]])

LinearRegression类是基于scipy.linalg.lstsq()函数(名字是指“标准二乘法(least squares)")，可以直接这样调用

theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd

array([[3.97884872],
       [2.92718847]])

这个函数计算 $\hat{\theta }=X^{+}y$ ，其中 $X^{+}$ 是 $X$ 的伪逆矩阵(穆尔-彭罗斯逆矩阵Moore-Penrose inverse)。你可以使用 np.linalg.pinv()直接计算伪逆矩阵

np.linalg.pinv(X_b).dot(y)

array([[3.97884872],
       [2.92718847]])

求伪逆矩阵本身采用了一种标准的矩阵因式分解方法叫做奇异值分解(Singular Value Decomposition)，这种方法可以把训练数据矩阵 $X$ 分解为三个矩阵 $U,\sum ,V^T$ 的乘积(更多的信息可以查看numpy.linalg.svg())。求伪逆矩阵的计算是 $X^{+}=V\sum U^T$。当 $X^{T}X$ 不可逆时，正规方程无法求解，而求伪逆矩阵仍然可以在这种情况运行。

计算复杂度(Computational Complexity)

正规方程需要计算 $X^{T}X$ 的转置，计算复杂度是 $(n+1)\times (n+1)$(其中 $n$ 是特征数量)。转置一个这种矩阵的计算复杂度一般是 $O(n^{2.4})$ 到 $O(n^3)$。

Scikit-Learn的LinearRegression类采用的SVD方法计算逆的复杂度是 $O(n^2)$，所以如果特征数变为两倍，则计算耗时变为原来的四倍

正规方程和SVD方法在特征数量变的很大时计算都非常缓慢。但是优点是计算复杂度与训练集实例数量的关系是线性增长的，计算复杂度是 $O(m)$ ，所以它们可以高效地处理训练集

梯度下降(Gradient Descent)

梯度下降是一种常用的优化算法，它能够在大范围的问题中找到最优解。梯度下降的思想是迭代微调参数来最小化代价函数。

梯度下降背后的思想是：开始时我们随机选择一个参数的组合$\left({\theta }_0,{\theta }_1,......,{\theta }_n\right)$，计算代价函数，然后我们寻找下一个能让代价函数值下降最多的参数组合。我们持续这么做直到找到一个局部最小值（local minimum），因为我们并没有尝试完所有的参数组合，所以不能确定我们得到的局部最小值是否便是全局最小值（global minimum），选择不同的初始参数组合，可能会找到不同的局部最小值。

使用梯度下降方法时，需要确保所以的特征尺度相似，否则梯度下降需要花很长的时间收敛。

批量梯度下降(Batch Gradient Descent)

为了实现梯度下降方法，需要计算每个模型参数 ${\theta }_j$ 对应代价函数的梯度。这叫做偏导数(partial derivative)。

$$\frac{\partial }{\partial {\theta }_j}MSE(\theta )=\frac{2}{m}\sum _{i=1}^m({\theta }^T{x}^{(i)}-y^{(i)})x_j^{(i)}$$

为可以使用以下公式直接计算所以的偏导数，而不是一个个地计算

$${\nabla }_{\theta }MSE(\theta )=\left\{\begin{array}{ccccc}a& b& c& d& e\\ f& g& h& i& j\\ k& l& m& n& o\\ p& q& r& s& t\end{array}\right\}=\frac{2}{m}{X}^T(X\theta -y)$$

注意到，这个公式在每个梯度下降迭代需要计算训练集中所有的 $X$！结果就是它在大数据集上十分缓慢。但是，梯度下降在特征维数很多时效果很好，用梯度下降训练一个有数百或数千维数特征的线性回归模型比使用正规方程或SVD(奇异值分解)快很多。

有了梯度下降向量后，当前点在上方，就直接往下方移动。意思就是直接从 $\theta$ 中减去 ${\nabla }_{\theta }MSE(\theta )$。这里就是学习率登场的时候了:将学习率 $\eta$ 与梯度向量相乘来决定下降步长。

$${\theta }^{(nextstep)}=\theta -\eta {\partial }_{\theta }MSE(\theta )$$

用以上的思想来快速地实现算法

eta = 0.1   # 学习率
n_iterations = 1000 
m = 100

theta = np.random.randn(2, 1)

for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta -eta * gradients

theta

array([[3.97884872],
       [2.92718847]])

• 找到最佳的学习率可以使用网格搜索的方法
• 迭代次数可以选择很大一个数量，当梯度向量变的很小时终止迭代

收敛速度(Convergence Rate):当代价函数是凸函数，斜率没有突变(例如MSE代价函数), 具有固定学习率的批量梯度下降算法最终会收敛在最优解，但是你或许需要等待一会儿:在一系列根据代价函数的形状定义的 $ϵ$ ，需要迭代 $O(1/ϵ)$ 的次数。如果为了更精确的解将界限(指 $ϵ)$ 变小10倍，那么算法或许需要运行10倍于之前的时间。

随机梯度下降(Stochastic Gradient Descent)

随机梯度下降(Stochastic Gradient Descent, SGD)算法每步迭代随机选择一个实例然后仅根据这个实例计算梯度。

• 因为每次迭代需要操作的数据很少,算法运行得更快
• 因为每次迭代只需要一个实例，所以可以在大数据集上训练(SGD可以作为外存算法)
• SGD运行时比较随机，在整体上是逐渐接近最小值，但是它会反复弹跳，随着时间的变化，它会到达最小值附加，但是它不会稳定下来，仍会反复横跳。所以在算法停止运行时，最终的结果很好，但是不是最优解。
• 当代价函数不平整时，这会让算法跳出局部最优解，所以比起BGD,SGD有更多的机会找到全局最优解

因此随机性帮助跳出局部最优，但是也意味着算法无法稳定在最小值。一种解决的方法是逐渐地减小学习率。步长开始的时候设置的很大(这有助于快速迭代和跳出局部最优)，然后步长慢慢变小，然算法在全局最优解稳定。这个步骤类似于模拟退火。每次迭代控制学习率的函数叫做学习计划(learning schedule)

• 学习率降低的太快，可能会一直陷在局部最优或者甚至未达到局部最优就停滞了
• 学习率降低的太慢，可能很长一段时间在最优解附近弹跳，如果你很早就终止算法可能会在次优解终止

以下代码是使用了一个简单学习计划的随机梯度下降算法的实现

n_epochs = 50
t0, t1 = 5, 50

def learning_schedule(t):
    return t0 / (t+t1)

theta = np.random.randn(2, 1)

for epoch in range(n_epochs):
    for i in range(m):
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta)-yi)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients

theta

array([[4.01896884],
       [2.968826  ]])

需要注意的一点是：实例的选择是随机的，一些实例在迭代过程中可能被多次选择而一些实例完全没有被选中。

• 为了保证参数趋向全局最优解，简单的方法是打乱训练集，然后随机的选择实例。否则可能出现一个标签一个标签的分类。

采用Scikit-Learn的SGDRegressor类实现SGD，例子如下所示

from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=1000, tol=1e-3, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())

SGDRegressor(eta0=0.1, penalty=None)

sgd_reg.intercept_, sgd_reg.coef_

(array([3.89348504]), array([2.86385718]))

sgd_reg

SGDRegressor(eta0=0.1, penalty=None)

小批量梯度下降(Mini-batch Gradient Descent)

Mini-batch GD在小部分的随机数据集实例上计算梯度。Mini-batch GD比起随机梯度下降的主要优点是可以通过硬件矩阵操作的优化获得性能增益，尤其在使用GPU的时候。

• Mini-batch GD在参数空间的变化比SGD变化少一些不稳定性，尤其是在大数据集下
• Mini-batch GD比起SGD会更近的靠近最小值
• Mini-batch GD会更难跳出局部最小值

多元回归(Polynomial Regression)

如果数据实际上比起简单的直线要复杂，也可以使用线性模型来拟合非线性数据。一个简单的方法就是增加特征的幂数，然后在这些扩增的数据特征上训练线性模型，这种技术叫做多元回归。

一个简单的平方等式生成的例子如下

m = 100
X = 6 * np.random.rand(m, 1) -3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)

数据生成了，画图查看数据的分布

plt.scatter(X, y)
plt.xlabel('X_1')
plt.ylabel('y')
plt.show()

很显然，直线不能很好地拟合数据。所以利用Scikit-Learn的PolynomialFeatures类来转化和训练数据，为每个特征添加平方项作为新的特征

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
X[0]

array([-2.80664181])

X_poly[0]

array([-2.80664181,  7.87723826])

X_poly现在包含了原来X特征和其对于的平方项特征。现在可以使用LinearRegression模型来拟合这个扩展后的数据集

lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_

(array([2.11424994]), array([[0.98144922, 0.4598956 ]]))

x_pre = np.arange(-4,3,0.01).reshape(-1,1)
x_ppre = poly_features.fit_transform(x_pre)
y_ppre = lin_reg.predict(x_ppre)
plt.scatter(X, y)
plt.plot(x_pre, y_ppre,'r-', label="Predictions")
plt.legend()
plt.xlim(-3, 3)
plt.xlabel("X1")
plt.ylabel("y")
plt.show()

得到的模型是 $\hat{y}=0.56x_1^2+0.93x_1+1.78$，实际的模型函数是 $y=0.5x^2+1.0x_1+2.0+Gaussiannoise$

学习曲线(Learning Curves)

实际情况下，模型过于复杂会导致过拟合，模型太简单会出现欠拟合。

第二章中利用交叉验证可以判断模型是过拟合或者欠拟合。

另一种方式是通过学习曲线(learning curves)：绘制模型在训练集和测试验证集上的表现。为了生成曲线，简单地在不同大小的训练子集上训练模型。

from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))
    plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")
    plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")

lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
plt.legend()
plt.show()

# 10th degree 多元线性模型
from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline([
    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
    ("lin_reg", LinearRegression()),
])

plot_learning_curves(polynomial_regression, X, y)
plt.ylim(0,3)
plt.legend()
plt.show()

该学习曲线和前者很相似，但是两者有两个很重要的不同：

• 这次的训练数据误差远低于线性回归模型
• 在两条曲线中间有间隔。这表示模型在训练数据上表现比校验误差好，这是过拟合的标志。但是，如果你用一个更大的数据集，两条曲线之间的间隔会逐渐变小

一种提高过拟合模型的方法是给模型喂更多的训练数据直到验证误差解决训练误差

偏差与方差的权衡: 统计学和机器学习领域很重要的一个理论结果是-模型的泛化误差可以表示为三个不同误差的和

• 偏差(Bias)：这部分的泛化误差来源于错误的假设，例如假设数据是线性的结果数据实际上是二次型的。高偏差的模型绝大多数欠拟合训练数据
• 方差(Variance)：这部分是由于模型对训练数据的小变化的过分敏感。有很高次方自由项的模型(例如高次方的多元线性模型)更有可能有高的方差，然后可能过拟合数据
• 不可避免误差(Irreducible error)：这部分是由于数据本身的噪声。减少这部分误差的唯一方法是清理数据(比如，修复数据源，确实的张量(tensor)，或检测和清楚异常值)

增加模型的复杂度通常会提升其方差并降低偏差。相反地，减少模型的复杂度会增加偏差和减少方差，这就是需要考虑权衡的原因。

正则化线性模型(Regularized Linear Models)

降低过拟合的一种方法的是正则化模型：自由度越低，过拟合程度越低。以下是三种不同的限制权重来降低过拟合的方法。

岭回归(Ridge Regression)

代价函数中添加了一个正则项 $a{\sum }_{i=1}^n{\theta }_i^2$. 正则化项只有在训练中添加，模型训练之后，评估模型的性能则使用未正则化的性能评估方法

岭回归代价函数的公式表示为：

$$J(\theta )=MSE(\theta )+a\frac{1}{2}\sum _{i=1}^n{\theta }_i^2$$

需要注意的是　${\theta }_0$ 没有正则化
在实施岭回归之前，标准化数据很重要，因为岭回归对与输入特征的尺度很敏感

岭回归的闭合解

$$\hat{\theta }=(X^{T}X+aA{)}^{-1}{X}^{Ｔ}y$$

以下是如何使用Scikit-Learn通过闭合解实现岭回归

from sklearn.linear_model import Ridge

ridge_reg = Ridge(alpha=1, solver="cholesky")
ridge_reg.fit(X, y)

Ridge(alpha=1, solver='cholesky')

ridge_reg.predict([[1.5]])

array([[5.0048275]])

# 使用随机梯度下降法
sgd_reg = SGDRegressor(penalty="l2")
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])

array([4.99518164])

套索回归(Lasso Regression)

最小绝对值收敛和选择运算回归(简称Lasso Regression)是线性回归的另一种正则化方法，它使用的是 $l_1$ 范数

*Lasso Regression代价函数

$$J(\theta )=MSE(\theta )+a\sum _{i=1}^n|{\theta }_i|$$

from sklearn.linear_model import Lasso
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
lasso_reg.predict([[1.5]])

array([4.95391855])

弹性网络(Elastic Net)

弹性网络是岭回归和Lasso回归的折中

弹性网络代价函数公式

$$J(\theta )=MSE(\theta )+ra\sum _{i=1}^n|{\theta }_i|+\frac{1-r}{2}a\sum _{i=1}^n{\theta }_i^2$$

• 岭回归一般作为默认的正则化，但是如果只有很少的特征真的有用，可以使用Lasso或Elastic，因为它们倾向减少无用的特征权重至0
• 通常情况下，Elastic比Lasso好，因为当数据的特征数比数据实例数多或数个特征有较强关联时，Lasso可能表现的不太稳定。

以下是使用Scikit-Learn的ElasticNet

from sklearn.linear_model import ElasticNet

elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5)
elastic_net.fit(X, y)
elastic_net.predict([[1.5]])

array([4.95584998])

及时停止(Early Stopping)

在一些正则化迭代学习算法中，在验证集误差达到最小时，及时停止是一个很好的方式

逻辑回归(Logistic Regression)

一些回归算法也可以被用于分类，如果评估到可能性高于50%，那么模型将这个实例划分为这个类，否则划分为另一类

评估可能性

逻辑线性回归不是直接输出结果，而是结果的逻辑

逻辑线性回归评估可能性的公式

$$\hat{p}=h_{\theta }(x)=\sigma (x^T\theta )$$

$\sigma$ 是sigmod function，输出的是0到1的数字，定义如下：

$$\sigma (t)=\frac{1}{1+exp(-t)}$$

训练和损失函数(Training and Cost Function)

逻辑回归的损失函数如下(log loss)

$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m[y^{(i)}log({\hat{p}}^{(i)})+(1-y^{(i)})log(1-\widehat{p^{(i)}})]$$

其关于 ${\theta }_j$ 的偏导式如下所示

$$\frac{\partial }{\partial {\theta }_j}J(\theta )=\frac{1}{m}\sum _{i=1}^m(\sigma ({\theta }^T{x}^{(i)})-y^{(i)})x_j^{(i)}$$

通过这个公式，你也可以实现批量梯度下降算法，SGD，Mini-batch GD。

决策边界(Decision Boundaries)

使用鸢尾花数据集来阐述逻辑回归

from sklearn import datasets

iris = datasets.load_iris()

list(iris.keys())

['data',
 'target',
 'frame',
 'target_names',
 'DESCR',
 'feature_names',
 'filename']

X = iris["data"][:, 3:] # petal width
y = (iris["target"] == 2).astype(np.int)

获取了数据集，然后训练逻辑回归模型

from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression()
log_reg.fit(X, y)

LogisticRegression()

查看模型对花瓣宽度从0到3cm的评估

X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
plt.plot(X_new, y_proba[:, 1], "g-", label="Iris-Varigica")
plt.plot(X_new, y_proba[:, 0], label="Not Iris-Virginica")

[<matplotlib.lines.Line2D at 0x7f422f7510f0>]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uIROTzgb-1596532786757)(chap4_files/chap4_72_1.svg)]

Softmax Regression

$$S_k(x)=x^T{\theta }^{(k)}$$

Equation 4-20.Softmax function

$${\hat{p}}_k=\sigma (s(x){)}_k=\frac{exp(s_k(x))}{{\sum }_{j=1}^{k}exp(s_k(x))}$$

X = iris["data"][:, (2, 3)]
y = iris["target"]

softmax_reg = LogisticRegression(multi_class="multinomial", solver="lbfgs", C=10)
softmax_reg.fit(X, y)

LogisticRegression(C=10, multi_class=‘multinomial’)

softmax_reg.predict([[5, 2]])

array([2])

softmax_reg.predict_proba([[5, 2]])

array([[6.38014896e-07, 5.74929995e-02, 9.42506362e-01]])

参考文献

1. 《Hands-on Machine Learning with Scikit-Learn, Keras & TensorFlow》(2nd Edition)
2. github 代码地址