梯度下降--从一元到多元
什么是梯度下降法
梯度下降法(英语:Gradient descent)是一个一阶最优化算法,通常也称为最速下降法。 要使用梯度下降法找到一个函数的局部极小值,必须向函数上当前点对应梯度(或者是近似梯度)的反方向的规定步长距离点进行迭代搜索。
一元梯度下降法
为了便于理解,笔者先从一元梯度下降算法开始讲起,如图所示的损失函数为一元二次方程,为了找到函数极小值,我们从随机初始化的点向梯度反方向不断搜索,最终找到近似的极小值点。
代码实现
import numpy as np
import matplotlib.pyplot as plt
# 原函数
def J(x):
return x**2 - x + 4
# 导函数
def dJ(x):
return 2*x - 1
# 梯度搜索
def gradient_descent(theta, eta, epsilon=1e-8):
# 保存训练轨迹
theta_history = []
while True:
gradient = dJ(theta)
last_theta = theta
theta = last_theta - eta * gradient
theta_history.append(theta)
if abs(J(last_theta) - J(theta)) < epsilon:
break
return np.array(theta_history)
if __name__ == '__main__':
x = np.linspace(-10, 10, 1000)
y = J(x)
history = gradient_descent(-10, 0.1)
print(history)
# 绘制曲线和梯度下降轨迹
plt.plot(x, y)
plt.scatter(history, J(history), color='red')
plt.show()
多元梯度下降法
以上实现的一元梯度下降法只是为了方便学习和理解,接下来笔者将实现多元梯度下降法,并用一个实例直观的感受下梯度下降法的强大之处。
如图,左图是方程z=2x+3y+4的函数图像,我们将采用梯度下降法从已知数据中拟合出这个二元一次方程。右图即是拟合的过程及拟合结果。
代码实现
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 损失函数
def J(theta, x, y):
return np.sum((y - np.dot(x, theta))**2) / len(x)
# 导函数
def dJ(theta, x, y):
res = np.empty(len(theta))
res[0] = np.sum(np.dot(x, theta) - y)
for i in range(1, len(theta)):
res[i] = np.sum(np.dot(np.dot(x, theta) - y, x[:, i]))
return res * 2 / len(x)
# 梯度搜索
def gradient_descent(theta, x, y, eta, epsilon=1e-8):
theta_history = []
while True:
gradient = dJ(theta, x, y)
last_theta = theta
theta = theta - gradient * eta
esp = abs(J(last_theta, x, y) - J(theta, x, y))
if esp < epsilon:
break
if esp > 0.5:
theta_history.append(theta)
theta_history.append(theta)
return np.array(theta_history)
if __name__ == '__main__':
# 准备数据
size = 100
np.random.seed(999)
x = np.linspace(-1, 1, size)
y = np.linspace(-1, 1, size)
x_b = np.hstack([np.ones((len(x), 1)), x.reshape(-1, 1), y.reshape(-1, 1)])
z = 2*x+3*y+4+np.random.uniform(-0.2, 0.2, size)
# 通过梯度下降搜索theta
eta = 0.05
init_theta = np.zeros(x_b.shape[1])
theta_history = gradient_descent(init_theta, x_b, z, eta)
X, Y = np.meshgrid(x, y)
fig = plt.figure()
# 绘制原始图像 z=2*x+3*y+4
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
Z = 2*X + 3*Y + 4 # 原函数
ax1.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
plt.title('z=2*x+3*y+4')
# 绘制预测图像
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
for his in theta_history:
Z = his[1] * X + his[2] * Y + his[0] # 预测函数
ax2.plot_surface(X, Y, Z, rstride=1, cstride=1)
plt.show()