- 问题引入
- 线性回归
- 梯度下降法实现(Python)
- 多元线性回归中的梯度下降法模拟
- 多元梯度下降法实现(Python)
- 多元梯度下降法实现——修改版(Python)
- 随机梯度下降法实现(Python)
- 梯度下降中的模拟求偏导的一般方法
- 梯度下降的总结
- 梯度下降的优化
机器学习系列之梯度下降。
问题引入
在机器学习中，我们通常会将求解抽象为一个损失函数，我们的目的是最小化这个损失函数。而梯度下降法和梯度上升法则是一个很好的工具。例如，对于二次函数：
对于二次函数，我们要求极值，一般我们求导函数，然后找能使导函数为0的点即为所求。对于倒数的精确定义，我们在这里不做精确叙述(参考高数书)。我们可以明确得到的信息如下：
- 1.在曲线中，导数代表切线的斜率
- 2.导数为负，θ减小，则J增加；θ增加，则J减小
- 3.导数为正，θ减小，则J减小；θ增加，则J增大
- 4.导数可以代表方向，对应了J增大的方向(导数为负，θ减小的方向J增大；导数为正，θ增大的方向J增大)
对于梯度下降法，我们采用步进的方式，一步一步向着极值前进。每次我们都在当前位置，求出下降最快的方向，然后按照设定的步长走一步。在上述例子中，横坐标是θ，我们每次都对横坐标θ进行修改，修改量为框中所示：
其中η代表的就是步长，也叫学习率，是一个大于0的数。如果导数为负，我们在前面添加一个符号，这样得到的新方向就是J减小的方向，这样θ的变化量即为：(-1)η(dJ/dθ)。随着不断的前进，我们会达到极值点。值得注意的是，η的取值将会影响获得最优解的速度，某些时候甚至得不到最优解。η是梯度下降法的一个超参数。下面我们来看看η的取值对收敛的影响：
此外，对于一些复杂的函数，可能有很多的极值点，因此我们对初始点、学习率的选择可能会使求解结果并非全局最优，即陷入局部最优。
解决上述问题的办法就是多次运行，并随机化初始点。
线性回归
线性回归的目的
对于上述给出的一系列散点(二维空间内)，我们的目的是找出一条直线来拟合这些点，反映出这些点大致的潜在规律。
- 1.最基本的单变量线性回归：h(x)=theta0+theta1*x1
- 2.复杂点的多变量线性回归：h(x)=theta0+theta1x1+theta2x2+theta3*x3
损失函数的形式
形式如下：
其中θi(0<=i<=n)代表了一个参数，对于n=2的情况，θ0代表了截距，θ1代表斜率。n>2的情况稍微复杂，可以进行类比。损失的定义形式如下：
其中m代表点的数量。假设我们已经有了一个模型，对于每个点，我们知道其真实值和由模型预测的值，我们将其做差后平方，这样就得到了该模型在这一个点的损失，然后我们累加所有点的损失就得到了损失函数。从实际意义出发，对于所有点的损失总和，其值越小，说明我们的模型拟合度越高。从几何角度看，我们的模型刻画出的线，使得所有数据点到该线的距离和最小。此外，在线性回归法的损失函数具有唯一的最优解。
梯度下降法实现(Python)
代码如下：
import numpy as np
import matplotlib.pyplot as plt
def dJ(theta):
try:
return 2 * (theta - 2.5)
except:
return float('inf')
def J(theta):
try:
return (theta - 2.5) ** 2 - 1
except:
return float('inf')
def gradient_descent(initial_theta, eta, n_iter = 1e4, epsilon=1e-8):
theta = initial_theta
theta_history = [initial_theta]
i_iter = 0
while i_iter < n_iter:
gradient = dJ(theta)
last_theta = theta
theta = theta - eta * gradient
theta_history.append(theta)
if (abs(J(theta) - J(last_theta)) < epsilon):
break
i_iter += 1
print("theta值为：" + str(theta))
print("J(theta)值为：" + str(J(theta)))
print("theta_histoy长度为：" + str(len(theta_history)))
return theta_history
def plot_theta_history(plot_x, theta_history):
plt.plot(plot_x, J(plot_x))
plt.plot(np.array(theta_history), J(np.array(theta_history)), color='r', marker='+')
plt.show()
if __name__ == "__main__":
# 在-1到6的范围内均匀分割141个点
plot_x = np.linspace(-1, 6, 141)
print(plot_x.shape)
eta = 1.0
theta_history = gradient_descent(0.0, eta)
plot_theta_history(plot_x, theta_history)
最终我们打印出相关解：
theta值为：2.499891109642585
J(theta)值为：-0.99999998814289
执行次数为：46
可以看出，结果和实际很接近。
可以尝试自行修改学习率eta的值，查看不同eta值对收敛的影响。
多元线性回归中的梯度下降法模拟
在之前的例子中θ是一个值(y=k*x+b中的x)，对于多元情况下，θ代表一个向量(θ=(θ0，θ1，θ2，θ3…，θn))，这样我们就要修改上面的求解形式：
在上述公式中，自然地引入了偏导。这里梯度代表了方向，对应J增加最快的方向。考虑下降以及多元的情况，在每个点会有多个下降方向，梯度是所有下降方向中，下降速度最快的方向。下面这个等高线的的图可以清楚地表现出步进的过程：
当我们使用θ向量(θ=(θ0，θ1，θ2，θ3…，θn))代入后，损失函数变为：
我们对每个θi求偏导并整理后，得到：
其中X_b代表的是数据矩阵，Xn代表了一个实例(是一个行向量)，这个实例可能有许多特征，每个特征都有其特征值。X_b矩阵等同于(X1,X2,X3,…,Xn)的转置。如下图所示：
此外，我们要对这个损失函数进行调整，添加系数，最终形式为：
多元梯度下降法实现(Python)
import numpy as np
import matplotlib.pyplot as plt
def dJ(theta, X_b, y):
res = np.empty(len(theta))
# 单独设定theta0对应的偏导值
res[0] = np.sum(X_b.dot(theta) - y)
# 对于theta1-thetan，迭代设置对应的偏导值
for i in range(1, len(theta)):
# 对于矩阵来说直接相乘即可，避免了单独求和并累计和，代表的含义参考图中公式即可
res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
return res * 2 / len(X_b)
def J(theta, X_b, y):
try:
return np.sum((y - (X_b).dot(theta)) ** 2) / len(X_b)
except:
return float('inf')
def gradient_descent(X_b, y, initial_theta, eta, n_iter = 1e4, epsilon=1e-8):
theta = initial_theta
i_iter = 0
while i_iter < n_iter:
gradient = dJ(theta, X_b, y)
last_theta = theta
# 更新
theta = theta - eta * gradient
# 判断是否收敛
if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
break
i_iter += 1
print("theta值为：" + str(theta))
print("J(theta)值为：" + str(J(theta, X_b, y)))
return theta
def plot_theta_history(plot_x, theta_history):
plt.plot(plot_x, J(plot_x))
plt.plot(np.array(theta_history), J(np.array(theta_history)), color='r', marker='+')
plt.show()
if __name__ == "__main__":
np.random.seed(666)
x = 2 * np.random.random(size=100)
y = x * 3. + 4. + np.random.normal(size=100)
X = x.reshape(-1, 1)
X_b = np.hstack([np.ones((len(x), 1)), X.reshape(-1, 1)])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01
theta = gradient_descent(X_b, y, initial_theta, eta)
print(theta)
# plt.scatter(x, y)
# plt.show()
最终我们打印出相关解：
theta值为：[ 4.02145786 3.00706277]
J(theta)值为：1.09887102571
结果与我们的定义的函数：y = x*3.+4.+np.random.normal(size=100)，系数十分接近。
多元梯度下降法实现——修改版(Python)
对于上面推导出的偏导公式，我们发现J(θ)的偏导第一项最右边的系数是-1，而其他偏导最右边的系数是X(i)，这就诱导我们将第一项的系数-1转化为X(i),过程如下：
在代码中变动的部分仅为函数dJ，如下：
def dJ(theta, X_b, y):
# res = np.empty(len(theta))
# res[0] = np.sum(X_b.dot(theta) - y)
# for i in range(1, len(theta)):
# res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
# return res * 2 / len(X_b)
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)
最终我们打印出相关解：
theta值为：[ 4.02145786 3.00706277]
J(theta)值为：1.09887102571
结果与之前变化不大。
随机梯度下降法实现(Python)
在前面推导出的公式中，我们发现在求梯度时，对于各个方向的偏导，我们在求解时需要对所有的样本进行计算，对诸如这种求解方式，我们称之为批量梯度下降。在实际情况中，我们可能会遇到m十分大的情况，那么求解时的计算量将会特别大。有没有办法遏制这种情况呢，这就引出了随机梯度下降的出现。随机，就是随机选取一个方向的导数作为梯度，而不考虑这个方向是不是所有方向里下降最快的方向，因此存在向上升方向前进的情况存在。
下图展示了向非下降方向前进的情况：
实现过程如下：
import numpy as np
import matplotlib.pyplot as plt
def J(theta, X_b, y):
try:
return np.sum((y - (X_b).dot(theta)) ** 2) / len(X_b)
except:
return float('inf')
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)
def dJ_SGD(theta, X_b_i, y_i):
return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2.
def SGD(X_b, y, initial_theta, n_iters):
t0 = 5
t1 = 50
# 模拟退火
def learning_rate(t):
return t0 / (t + t1)
theta = initial_theta
for cur_iter in range(n_iters):
# 随机一个下标
rand_i = np.random.randint(len(X_b))
gradient = dJ_SGD(theta, X_b[rand_i], y[rand_i])
theta = theta - learning_rate(cur_iter) * gradient
return theta
def plot_theta_history(plot_x, theta_history):
plt.plot(plot_x, J(plot_x))
plt.plot(np.array(theta_history), J(np.array(theta_history)), color='r', marker='+')
plt.show()
if __name__ == "__main__":
m = 100000
x = 2 * np.random.random(size=m)
y = x * 3. + 4. + np.random.normal(0, 3, size=m)
X = x.reshape(-1, 1)
X_b = np.hstack([np.ones((len(x), 1)), X])
# print(X_b)
initial_theta = np.zeros(X_b.shape[1])
# print(initial_theta)
eta = 0.01
theta = SGD(X_b, y, initial_theta, len(X_b))
print(theta)
# plt.scatter(x, y)
# plt.show()
梯度下降中的模拟求偏导的一般方法
对于一个复杂的函数，如果不能直接定义其导函数，那么怎么求梯度呢？我们可以采用高数中的一些定义近似求解梯度。
这种方法在高维情况下依旧适用：
实例代码如下：
import numpy as np
import matplotlib.pyplot as plt
def J(theta, X_b, y):
try:
return np.sum((y - (X_b).dot(theta)) ** 2) / len(X_b)
except:
return float('inf')
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)
def dJ_math(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)
def dJ_debug(theta, X_b, y, epsilon=0.01):
res = np.empty(len(theta))
for i in range(len(theta)):
theta_1 = theta.copy()
theta_1[i] += epsilon
theta_2 = theta.copy()
theta_2[i] -= epsilon
res[i] = (J(theta_1, X_b, y) - J(theta_2, X_b, y)) / (2 * epsilon)
return res
def gradient_descent(dJ, X_b, y, initial_theta, eta, n_iter = 1e4, epsilon=1e-8):
theta = initial_theta
i_iter = 0
while i_iter < n_iter:
gradient = dJ(theta, X_b, y)
last_theta = theta
# 更新
theta = theta - eta * gradient
# 判断是否收敛
if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
break
i_iter += 1
print("theta值为：" + str(theta))
print("J(theta)值为：" + str(J(theta, X_b, y)))
return theta
if __name__ == "__main__":
np.random.seed(666)
X = np.random.random(size=(1000, 10))
true_theta = np.arange(1, 12, dtype=float)
X_b = np.hstack([np.ones((len(X), 1)), X])
y = X_b.dot(true_theta) + np.random.normal(size=1000)
inital_theta = np.zeros(X_b.shape[1])
eta = 0.01
theta = gradient_descent(dJ_debug, X_b, y, inital_theta, eta)
print(theta)
梯度下降的总结
一般来说，梯度下降有以下几种实现方式：
- 1.批量梯度下降法(Batch Gradient Descent, BGD)
- 2.随机梯度下降法(Stochastic Gradient Descent, SGD)
- 3.小批量梯度下降法(Mini-batch Gradient Descent, MBGD)
批量梯度下降法是梯度下降法最原始的形式，它的具体思路是在更新每一参数时都使用所有的样本来进行更新。优点：全局最优解；易于并行实现；缺点：当样本数目很多时，训练过程会很慢。
随机梯度下降法的思路是在更新每一参数时都使用一个样本来进行更新。每一次更新参数都用一个样本，更新很多次。如果样本量很大的情况（例如几十万），那么可能只用其中几万条或者几千条的样本，就已经将theta迭代到最优解了，对比上面的批量梯度下降，迭代一次需要用到十几万训练样本，一次迭代不可能最优，如果迭代10次的话就需要遍历训练样本10次，这种更新方式计算复杂度太高。但是，SGD也带来一个问题，那是噪音较BGD要多，使得SGD并不是每次迭代都向着整体最优化方向。优点：训练速度快；缺点：准确度下降，并不是全局最优；不易于并行实现。
小批量梯度下降法的具体思路是在更新每一参数时都使用一部分样本来进行更新。为了克服上面两种方法的缺点，又同时兼顾两种方法的有点。
相对于随机梯度下降，小批量梯度下降降低了收敛波动性，即降低了参数更新的方差，使得更新更加稳定。相对于批量梯度下降，其提高了每次学习的速度。并且其不用担心内存瓶颈从而可以利用矩阵运算进行高效计算。一般而言每次更新随机选择50-256个样本进行学习，但是也要根据具体问题而选择，实践中可以进行多次试验，选择一个更新速度与更次次数都较适合的样本数。小批量梯度下降可以保证收敛性，常用于神经网络中。
问题与挑战：
- 1.选择一个合理的学习速率很难。如果学习速率过小，则会导致收敛速度很慢。如果学习速率过大，那么其会阻碍收敛，即在极值点附近会振荡。
- 2.学习速率调整(又称学习速率调度，Learning rate schedules)试图在每次更新过程中，改变学习速率，如退火。一般使用某种事先设定的策略或者在每次迭代中衰减一个较小的阈值。无论哪种调整方法，都需要事先进行固定设置，这边便无法自适应每次学习的数据集特点。
- 3.模型所有的参数每次更新都是使用相同的学习速率。如果数据特征是稀疏的或者每个特征有着不同的取值统计特征与空间，那么便不能在每次更新中每个参数使用相同的学习速率，那些很少出现的特征应该使用一个相对较大的学习速率。
- 4.对于非凸目标函数，容易陷入那些次优的局部极值点中，如在神经网路中。那么如何避免呢。Dauphin指出更严重的问题不是局部极值点，而是鞍点(鞍点（saddle point）的数学含义是：目标函数在此点上的梯度（一阶导数）值为 0，但从该点出发的一个方向是函数的极大值点，而在另一个方向是函数的极小值点)。
梯度下降的优化
未完待续
- 1.Nesterov accelerated gradient（NAG）
- 2.Adagrad
- 3.Adam
- 4.Hogwild
- 5.Downpour SGD
- 6.Delay-tolerant Algorithms for SGD
- 7.Elastic Averaging SGD
- 8.Batch normalization