5.梯度下降法
5.1什么是梯度下降法
5.1.1梯度下降法的实现思想
最优化方法是寻找函数(目标函数)极值点的数值方法,包括梯度下降法、牛顿法、坐标下降法、拉格朗日乘数法、凸优化等。
梯度下降法是一种基于搜索的最优化方法,它不是机器学习算法;其作用是最小化一个损失函数。
梯度是导数对多元函数的推广,它是多元函数对各个自变量偏导数形成的向量,常表示为▽f。
一元函数y = f(x)的梯度是df/dx,其在x0处的梯度为df/dx|x=x0。
多元函数y = f(x),x = (x1,x2,x3,...xn)的梯度是:
需要注意的是,在机器学习中梯度都是列向量。
梯度下降法就是沿着梯度向量反方向进行迭代以达到函数极值点的方法。在这里直接给出公式,不再推导。从初始点x0开始,使用如下的迭代公式:
其中,η是一个接近于0的正数(一般设置为0.01),称为学习率。其取值会影响获得最优解的速度;若η取值不合适,甚至得不到最优解。η是梯度下降法中的一个超参数。
当f(xk+1)-f(xk)无限接近于0时,可以认为xk+1为极小值点;否则函数值xk+1会继续沿着序列xk+2的方向进行搜索,最终找到梯度为0的点。
当一个函数有多个极值点时,获得的极值点可能是局部极值点;为了获得全局极值点(最小值点),可以多次随机化初始点来重新运行程序。因此,初始点也是梯度下降法的一个超参数。
5.1.2简单实现梯度下降法
#程序5-1:求其y = (x-1)**2 + 2的极小值点
import numpy as np
def dy(x):
return 2*(x-1)
def y(x):
return (x-1)**2 + 2
#设置初始点x_0
x = 0.0
#设置eta
eta = 0.01
#设置近0值
near_zero = 1e-10
while True:
grad_x = dy(x)
last_x = x
x = x - eta*grad_x
if abs(y(last_x)-y(x)) < near_zero:
break
print(x)
print(y(x))
运行结果:
0.9999507956878463
2.0000000024210642
在梯度下降法中,η的取值非常关键,若η取值过大,可能执行x = x - eta*grad_x时导致x值不断增大/减小,使y(last_x)和y(x)趋向于无穷大/小。因此,修改代码如下。
#程序5-2:修改eta
import numpy as np
def dy(x):
return 2*(x-1)
def y(x):
try:
return (x-1)**2 + 2
except:
return float('inf') #返回最大的浮点数
#设置初始点x_0
x = 0.0
#设置eta
eta = 1.11
#设置近0值
near_zero = 1e-10
cur_iters = 0
max_iters = 1e4
while cur_iters < max_iters:
grad_x = dy(x)
last_x = x
x = x - eta*grad_x
if abs(y(last_x)-y(x)) < near_zero:
break
cur_iters += 1
print(x)
print(y(x))
运行结果:
nan
nan
当x和y(x)取值过大时,会返回nan。
5.2梯度下降法和线性回归
梯度下降法分为批量梯度下降法(BGD)、随机梯度下降法(SGD)和小批量梯度下降法(MBSD)。
5.2.1批量梯度下降法
在前面多元线性回归算法中,需要求损失函数的极小值,即:
为了使损失函数和样本数无关,因此修改上式为:
式中,X(i)表示样本特征的数据集,y(i)表示样本标签的数据集,都是已知量;我们需要求当θ为何值时,损失函数J最小。
上式还可以用向量的形式来表示,在这里不再写推导过程,如下所示:
#程序5-3
import numpy as np
from sklearn import datasets
#波士顿房价
boston = datasets.load_boston()
print(boston.keys())
#以每个住宅的平均房间数(RM)为特征
X = boston.data.copy()
y = boston.target.copy()
#取出y的边界值
X = X[y < np.max(y)]
y = y[y < np.max(y)]
#划分训练集和测试集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=123)
print(X_train.shape)
print(y_train.shape)
#数据归一化
from sklearn.preprocessing import StandardScaler
std_sca = StandardScaler()
std_sca.fit(X_train)
X_train_std = std_sca.transform(X_train)
X_test_std = std_sca.transform(X_test)
#求X_b
X_train_std_b = np.hstack([np.ones(shape=(len(X_train_std),1)),X_train_std])
X_test_std_b = np.hstack([np.ones(shape=(len(X_test_std),1)),X_test_std])
print(X_train_std_b.shape)
print(X_test_std_b.shape)
#求theta
#方法1
#设置theta初始值
# theta = np.zeros(X_train_std_b.shape[1])
# #设置eta
# eta = 0.01
# #设置近0值
# near_zero = 1e-10
# def J(theta,X_b,y):
# return np.sum((y-X_b.dot(theta))**2)/len(X_b)
# def dJ(theta,X_b,y):
# ret = np.zeros(len(theta))
# for i in range(0,len(theta)):
# ret[i] = np.sum((X_b.dot(theta)-y).dot(X_b[:,i]))
# return ret*2/len(X_b)
# while True:
# grad_J = dJ(theta,X_train_std_b,y_train)
# last_theta = theta
# theta = theta - eta*grad_J
# if abs(J(last_theta,X_train_std_b,y_train)-J(theta,X_train_std_b,y_train)) < near_zero:
# break
# print(theta)
#方法2
#设置theta初始值
theta = np.zeros(X_train_std_b.shape[1])
#设置eta
eta = 0.01
#设置近0值
near_zero = 1e-10
def J(theta,X_b,y):
return np.sum((y-X_b.dot(theta))**2)/len(X_b)
def dJ(theta,X_b,y):
return X_b.T.dot(X_b.dot(theta)-y)*2/len(X_b)
while True:
grad_J = dJ(theta,X_train_std_b,y_train)
last_theta = theta
theta = theta - eta*grad_J
if abs(J(last_theta,X_train_std_b,y_train)-J(theta,X_train_std_b,y_train)) < near_zero:
break
print(theta)
#预测
y_predict = X_test_std_b.dot(theta)
#R**2
from sklearn.metrics import r2_score
R_2 = r2_score(y_test,y_predict)
print(R_2)
运行结果:
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])
(392, 13)
(392,)
(392, 14)
(98, 14)
[21.62933673 -0.95284124 0.55299072 -0.30574477 -0.03944411 -1.3706543
2.61353192 -0.82432509 -2.36439855 2.02745115 -2.18335363 -1.76904847
0.74384624 -2.25717465]
0.8007705984940655
由于步长η和量纲之间有关联,因此梯度下降法需要对数据进行归一化处理。
多元线性回归算法处理样本特征多、样本数据集大的情况,其效率低下。批量梯度下降法适合处理样本特征较多的情况;而随机梯度下降法则可以处理样本特征多且样本数多的情况。
5.2.2随机梯度下降法
在批量梯度下降法中,需要对所有的样本求取损失函数J的梯度,因此当样本数据集较大时,运算效率低。
因此在随机梯度下降法中,对样本数据集进行随机抽取,使用抽取的随机样本来求取“梯度”。由于随机样本和数据集的样本总数并无关系,因此修改批量梯度下降法中的梯度公式,如下所示:
可以简单的认为,对于任何一个样本都可以用上式来求随机梯度,虽然并不是100%准确。
随机梯度下降法的本质就是对于一个初始值θ,不断的往梯度的反方向变化,每次变化的长度称为步长η;而步长η模拟退火的思想,其等于a/(iters+b),其中iters表示迭代的次数,正常情况下a取值为5,b取值为50(经验法)。
#程序5-4
import numpy as np
from sklearn import datasets
#波士顿房价
boston = datasets.load_boston()
print(boston.keys())
#以每个住宅的平均房间数(RM)为特征
X = boston.data.copy()
y = boston.target.copy()
#取出y的边界值
X = X[y < np.max(y)]
y = y[y < np.max(y)]
#划分训练集和测试集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=123)
print(X_train.shape)
print(y_train.shape)
#数据归一化
from sklearn.preprocessing import StandardScaler
std_sca = StandardScaler()
std_sca.fit(X_train)
X_train_std = std_sca.transform(X_train)
X_test_std = std_sca.transform(X_test)
#求X_b
X_train_std_b = np.hstack([np.ones(shape=(len(X_train_std),1)),X_train_std])
X_test_std_b = np.hstack([np.ones(shape=(len(X_test_std),1)),X_test_std])
print(X_train_std_b.shape)
print(X_test_std_b.shape)
#求theta的方法:
def dJ_sgd(theta, X_b_i, y_i):
#随机梯度
return 2. * X_b_i.T.dot(X_b_i.dot(theta) - y_i)
def sgd(X_b, y, initial_theta, n_iters=5, t0=5, t1=50):
#学习率eta = t0/(t+t1),由于t不断增大,eta不断变小
def learn_rate(t):
return t0 / (t + t1)
#设置初始值theta,其不断往梯度的反方向变化,每次变化步长为eta
theta = initial_theta
#求样本数
m = len(X_b)
#n_iters表示对样本数据集的遍历总次数,i_iter表示对样本数据集遍历到第几次
for i_iter in range(n_iters):
#将训练数据集进行乱序处理
indexes = np.random.permutation(m)
X_b_new = X_b[indexes,:]
y_new = y[indexes]
for i in range(m):
#对数据集中的每个样本都求取梯度
#使用批量梯度下降法,遍历次数为m*(n+1)
#使用随机梯度下降法,遍历次数为m*n_iters
grad = dJ_sgd(theta, X_b_new[i], y_new[i])
theta = theta - learn_rate(i_iter * m + i) * grad
return theta
#设置theta初始值
initial_theta = np.zeros(X_train_std_b.shape[1])
#获取theta的数值
theta = sgd(X_train_std_b, y_train,initial_theta, n_iters=100)
print(theta)
#预测
y_predict = X_test_std_b.dot(theta)
#R**2
from sklearn.metrics import r2_score
R_2 = r2_score(y_test,y_predict)
print(R_2)
运行结果:
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])
(392, 13)
(392,)
(392, 14)
(98, 14)
[21.62973169 -0.9612318 0.56728355 -0.26750109 -0.04311416 -1.38761527
2.59674183 -0.81831922 -2.36458287 2.12518371 -2.2918619 -1.77474864
0.7440291 -2.26267114]
0.8012609768467776
随着n_iters数值越来越大,其θ值会越来越准确。在sklearn.linear_model中封装了SGD,其算法的实现更加复杂,θ值也更加准确。官方文档:https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html#sklearn.linear_model.SGDRegressor
#程序5-4
import numpy as np
from sklearn import datasets
#波士顿房价
boston = datasets.load_boston()
print(boston.keys())
#以每个住宅的平均房间数(RM)为特征
X = boston.data.copy()
y = boston.target.copy()
#取出y的边界值
X = X[y < np.max(y)]
y = y[y < np.max(y)]
#划分训练集和测试集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=123)
print(X_train.shape)
print(y_train.shape)
#数据归一化
from sklearn.preprocessing import StandardScaler
std_sca = StandardScaler()
std_sca.fit(X_train)
X_train_std = std_sca.transform(X_train)
X_test_std = std_sca.transform(X_test)
#使用sklearn.
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=100, tol=1e-4)
sgd_reg.fit(X_train_std,y_train)
score = sgd_reg.score(X_test_std,y_test)
print(score)
运行结果:
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])
(392, 13)
(392,)
0.7939722013627373
在最近版本中,将使用max_iter和tol代替n_iters。