醉后骑龙吹铁笛

酒醒不知何处

概述

一般情况下，我们不需要了解一个模型内部的工作细节。但是理解一个系统是如何工作的，有助于找到合适的模型、正确的训练算法，以及一套合适的超参数，而且能让我们更高效地执行错误调试和错误分析。

训练模型方法

闭式方程

直接计算出最适合训练集的模型参数（也就是使训练集上的成本函数最小化的模型参数）

梯度下降（迭代优化）

通过迭代的方法逐渐调整模型参数直至训练集上的成本函数调至最低，最终趋于同第一种方法计算出来的模型参数。梯度下降（GD）包括：批量梯度下降、小批量梯度下降以及随机梯度下降。

线性回归

基本理论

线性模型就是对输入特征的加权求和，再加上一个偏置相（截距项）的常数，以此进行预测。

线性模型公式：

• \(\hat{y}\)是预测值
• n 是特征数量
• \(x_i\)是第i个特征值
• \(\theta_j\)是第j个模型参数

向量化线性回归模型：

• \(\theta\)是模型的参数向量，包括偏置项、特征权重
• \(\theta^T\) 是\(\theta\)的转置向量（为行向量，不再是列向量）
• x是实例的特征向量，其中\(x_0\)始终是1
• \(\theta^T·x\)是\(\theta^T\)和x的点积
• \(h_\theta\)是使用模型参数\(\theta\)的假设函数

训练模型的过程就是设置模型参数直到模型最适应训练集的过程，即找到最小化的RMSE（均方根误差）或MSE（均方误差）

线性回归模型的MSE成本函数：

标准方程

为得到使成本函数最小的\(\theta\)值，有一个闭式解方法：

• \(\hat{\theta}\)是使成本函数最小的参数值
• y则是目标值向量

生成一些随机数，再由这些随机数生成线性数据, y = 4+3x
计算出其闭式解

可见非常接近真实值，模型如下
Scikit-Learn的结果也一样

X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.rand(100,1)
X_b = np.c_[np.ones((100,1)), X]
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

plt.scatter(y=y,x=X)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('random.png')

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)
plt.plot(X,y,'b.')
plt.plot(X_new,y_predict,'r-')
plt.savefig('model.png')

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

计算复杂度
标准方程求逆的矩阵\(X^T·X\)，是一个\(n \times n\)的矩阵。对这种矩阵求逆的时间复杂度通常在特征数量的5.3～8倍之间，特征数量比较大时，计算时间将及其缓慢。

梯度下降

基本原理

梯度下降是一种通用的优化算法，能够为大范围的问题找到最优解。其核心思想就是迭代地调整参数从而使成本函数最小化。具体来说，就是首先使用一个随机的\(\theta\)值（随机初始化），然后逐步改进，知道算法收敛出一个最小值。

梯度中的一个重要参数是每一步的步长，这取决于超参数的学习率，如果学习率太低，算法需要大量时间才能收敛
如果学习率太高，则会导致算法发散，无法找到最优解

PS:线性回归模型的MSE成本函数是一个凸函数，这意味着连接曲线上任意两个点的线段永远不会跟曲线相交，即不存在局部最小，只有全局最小。它同时还是一个连续函数，所以斜率也不会产生陡峭的变化。这保证了在学习率不是太高的情况下，即便是乱走，等的时间够长，都可以找到全局最小值，这规避了梯度下降陷阱
梯度下降陷阱

批量梯度下降

要实现梯度下降，我们需要计算\(\theta_j\)改变，成本函数会改变多少，这是关于 \(\theta_j\)的偏导数

如不想单独计算这些梯度，可以利用公式进行一次性计算

PS:此公式每一次计算梯度时，都是基于完整训练集的
那么此时步长就等于学习率（\(\eta\)）乘以梯度向量

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

学习率为0.1时的学习情况

学习率为0.02时的学习情况
学习率为0.5时的学习情况

可知学习率为0.02时太低，但还是可以找到最优解；学习率为0.5时就太高了，每次迭代都在偏离最优解。
要找到合适的学习率，可以通过网格搜索，但需要限制迭代次数。通常方法：设置一个非常大的迭代次数，但当梯度下降的值变得很微小时中断算法——即其泛数低于\(\epsilon\)（容差）时中断。

随机梯度下降

批量梯度下降每次都需要整个训练集进行计算，训练集很大时，算法会特别慢。而随机梯度下降每一步则是从训练集中随机选取一个实例，基于该单个实例来计算梯度。但是由于算法性质，其成本函数将不再是缓缓降低直到达到最小值，而是不断上上下下，但从整体来看，最终还是会接近最小值，不过即使达到了最小值，算法依旧还会反弹，不会停止。此算法得到的参数是足够好的，但不是最优的。
当成本函数不规则时，随机梯度下降其实可以帮助算法跳出局部最小值，因此相较于批量梯度下降，它对找到全局最小值更有优势。
随机性虽然可以逃离局部最优，但始终找不到最小值，有一个办法可以解决这个困境——逐步降低学习率。开始时步长比较大，后来越来越小，让算法尽量靠近全局最小值，此过程叫做模拟退火。确定每个迭代学习率的函数叫做学习计划。如果学习率降得太快，可能会陷入局部最小值，甚至停留在途中；如果学习率降的太慢，耗时比较长。

下面是一个简单的学习计划：

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

才短短50次，结果就已经很不错了，学习情况如下

Scikit-Learn中使用SGDRegressor

sgd_reg = SGDRegressor(n_iter=50, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())
sgd_reg.intercept_, sgd_reg.coef_

小批量梯度下降

每一步的梯度计算既不是基于整个训练集，也不是基于单个实例，而是基于一小部分随机实例。小批量梯度下降的优势在于可以从矩阵运算的硬件优化中获得显著的性能提升，特别是要用到GPU时。
但是，这个算法可能更难匆局部最小值陷阱中逃脱。线性回归算法比较：

多项式回归

基本原理

可以用线性模型来拟合非线性数据，一个简单的方法是将每个特征的幂指次方添加为一个新特征，然后在这个拓展过的数据集上训练线性模型，此即为多项式回归

基于简单二次方程(\(y=0.5x^2+x+2\))制造一些非线性数据

m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m,1)
plt.scatter(y=y,x=X)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('aitiLinear.png')

ploy_features = PolynomialFeatures(degree=2, include_bias=False)
X_ploy = ploy_features.fit_transform(X)

现在X_ploy包含原本的特征X和X的平方，对此进行线性回归
结果还算理想

学习曲线

高阶多项式回归对训练数据的拟合，很可能会比简单线性回归要好。但是也容易过拟合，我们可以使用交叉验证来判断，也可以观察学习曲线进行判断。学习曲线是绘制模型在训练集和验证集上，关于“训练集大小”的性能函数。要生成这个函数，只要在不同大小的训练子集上多次训练模型就行了。

def plot_learning_curves(model, X, y, file_name):
	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_predict, y_train[:m]))
		val_errors.append(mean_squared_error(y_val_predict, y_val))
	plt.plot(np.sqrt(train_errors), 'r-+', linewidth=2, label='train')
	plt.plot(np.sqrt(val_errors), 'b-', linewidth=3, label='val')
	plt.legend()
	plt.savefig(file_name)

一开始，训练集中只有1个实例，模型可以完美拟合，因此曲线从0开始。但是随着新实例的加入，模型不再能够完美拟合，误差一路上升，直到达到一个高地，从这一点开始，添加新示例，不再能够让平均误差上升或下降。验证集亦是相对应的，训练集很少时，模型不能很好地泛化，所以验证集一开始误差非常大，逐渐降低，直到靠近一条直线，不再变化。
再看同一数据集上5阶多项式的学习曲线

PS:如果两条曲线间有一定距离，则是过拟合表现，说明此模型在训练数据上表现好于验证数据

正则线性回归

基本原理

减少过拟合的一个好办法就是对模型进行正则化（约束模型）处理：拥有的自由度越低，越不容易过拟合。比如将多项式模型正则化的简单方法就是降低多项式阶数。

岭回归

吉洪诺夫正则化，在成本函数中添加一个\(\alpha\sum^{n}_{i=1}\theta^{2}_{i}\)的正则项，这使得学习中的算法不仅需要拟合数据，同时还需要让模型权重保持最小。

PS：正则项只能在训练的时候添加到成本函数中，一旦训练完成，需要使用未经正则化的性能指标来评估模型性能

\(\alpha\)控制模型正则化程度，如果其为0，那么模型就是线性模型；如果其值过大，那么所有权重都将非常接近0，结果是一条穿过数据平均值的水平线。

岭回归成本函数：

这里偏置项\(\theta_0\)没有正则化，求和从1而非0开始，正则项即是\(1/2(||w||_2)^2\),\(||w||_2\)是权重向量的\(l_2\)范数

PS：执行岭回归之前，必须对数据进行缩放，因为它对特征大小非常敏感

岭回归的闭式解：

A为nxn的单位矩阵

rigge_reg = Ridge(alpha=1, solver='cholesky')
rigge_reg.fit(X, y.ravel())
rigge_reg.predict([[1.5]])

sgd_reg = SGDRegressor(penalty='l2')
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])

套索回归

最小绝对收缩和选择算子回归（Lasso回归），也是向成本函数增加一个正则项，但其增加的是权重向量的\(l_1\)范数，而不是\(l_2\)范数平方的一半，Lasso回归成本函数：

Lasso回归的重要特点是倾向于完全消除掉最不重要的特征的权重（设为0）

右图中虚线看起来像二次函数，但接近线性，因为其高次项的特征权重为0，Lasso回归会自动进行特征筛选并输出一个稀疏模型（只有很少的特征有非0权重）
当\(\theta_i=0(i=1,2,…,n)\)时，Lasso成本函数是不可微的，但当任意\(\theta_i=0\)时，可以使用次梯度向量作为替代，如下公式：

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

sgd_reg = SGDRegressor(penalty='l1')
sgd_reg.fit(X, y)
sgd_reg.predict([[1.5]])

弹性网络

介于岭回归和Lasso回归之间，是两者的组合，组合比例按r来控制，r=0时，为岭回归；r=1时，为Lasso回归。弹性网络成本函数：

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

PS:通常来说，有正则化总比没有好，哪怕很小。岭回归是不错的默认选择，如果实际用到的特征比较少，倾向于Lasso回归和弹性网络，其中弹性网络优于Lasso回归，因为特征数超过训练实例数量或几个特征强相关时，Lasso回归很不稳定

早期停止法

验证误差达到最小值时停止训练

sgd_reg = SGDRegressor(n_iter=1, warm_start=True, penalty=None, learning_rate='constant', eta0=0.0005)
minimum_val_error = float("inf")
best_epoch = None
best_model = None
for epoch in range(1000):
	sgd_reg.fit(X_tarin, y_train)
	y_val_predict = sgd_reg.predict(X_val)
	val_error = mean_squared_error(y,y_val_predict)
	if val_error < minimum_val_error:
		minimum_val_error = val_error
		best_epoch = epoch
		best_model = clone(sgd_reg)

逻辑回归

估算一个实例属于某一类别的概率

概率估算

逻辑回归也是计算输入特征的加权和（加上偏置项），但其输出的是结果的数理逻辑：

逻辑模型是一个sigmoid函数（S形），写作\(\sigma(·)\)，它的输出是0～1之间的一个数字，逻辑函数如下：


预测逻辑如下：

训练和成本函数

单个实例成本函数：

当t接近0时，-log(t)会非常大，如果模型预测一个实例结果为正类结果接近0或负类结果接近1时，成本会非常大。
逻辑回归成本函数（log损失函数）：

坏消息是此函数没有闭式解，但好在是凸函数，可以通过梯度下降或其他优化算法求得最优解，逻辑回归成本函数偏导数：

决策边界

利用鸢尾花数据集进行说明

iris = datasets.load_iris()
list(iris.keys())
X = iris['data'][:,3:]
y = (iris['target'] == 2).astype(np.int)

log_reg = LogisticRegression()
log_reg.fit(X, y)

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-Virginica')
plt.plot(X_new, y_proba[:,0], 'b--', label='Not Iris-Virginica')
plt.legend()
plt.savefig('logistic.png')

横坐标为花瓣宽度，纵坐标为概率，大约1.6cm处是决策边界

决策边界处即为概率50%处，Logistic回归也可正则化，默认为使用\(l_2\)惩罚函数

Softmax回归

逻辑回归模型推广后可以直接支持多个类别，而不需要训练多个组合二元分类器，即Softmax回归（多元逻辑回归）
原理是首先计算出每个类别k的分数\(s_k(x)\)，然后对这些分数应用Softmax函数（归一化函数），估算出每个类别的概率。\(s_k(x)\)公式：

每个类别都有自己特定的参数向量\(\theta_k\)，所有这些向量通常存储在参数矩阵\(\Theta\)中，之后应用Softmax函数得到概率：

• K是类别数量
• \(s(x)\)是实例x每个类别的分数向量
• \(\sigma(s(x))_k\)是给定类别分数下，实例x属于类别k的概率

Softmax回归分类器预测：

argmax返回使函数最大化时对应变量的值，即返回类别k

PS：Softmax回归分类器一次只会预测一个类，它是多类别分类器，但不是多输出分类器

利用使交叉熵成本函数最小化的方法对模型进行训练，交叉熵成本函数：

• 如果第i个实例的目标类别为k，则\(y^{(i)}_k\)等于1，否则为0
• 当只有两个类别时（K=2），其等价于逻辑回归成本函数

对于类别k的交叉熵梯度向量：

X = iris['data'][:,(2,3)]
y = iris['target']
softmax_reg = LogisticRegression(multi_class='multinomial', solver='lbfgs', C=10)
softmax_reg.fit(X, y)
softmax_reg.predict([[5,2]])
softmax_reg.predict_proba([[5,2]])

后记

范数

• \(l_0\)范数
  向量中非零元素个数
• \(l_1\)范数
  向量中非零元素绝对值之和
• \(l_2\)范数
  向量元素的平方和再开平方

参考

Markdown数学公式语法

• Previous
  SSH远程连接Linux图形用户界面（GUI）
• Next
  机器学习（五）