梯度下降的原理、代码、调试

大家好，我是东哥。

作者简介：前审计 | 现风控模型工程师 | 数据挖掘

目  录

• 一、基本原理

• （一）什么是梯度？

• （二）为什么梯度方向（正、反）是函数值局部变化（上升、下降）最快的方向？

• （三）梯度下降法

• （四）补充：梯度上升法

• 二、模拟实现与原理可视化展示

• （一）单一参数

• （二）双参数

• 三、线性回归中的梯度下降法

• （一）基本原理

• （二）分步实现

• （三）封装

• 四、梯度下降法中数据归一化的必要性

• 五、梯度下降法相对于正规方程法的优势

• 六、随机“梯度”下降法

• （一）基本原理

• （二）分步实现

• （三）封装

• （四）scikit - learn 中的随机梯度下降法

• 七、梯度的调试

• 八、总结

• （一）梯度下降法分类

• （二）随机的力量

# 导入基本包
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import *
from mpl_toolkits import mplot3d #用于绘制3D图形

# 显示所有过程结果
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'

plt.rcParams['font.sans-serif']=['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号


一、基本原理

（一）什么是梯度？

1、从偏导数说起

（1）一元函数情况

对于一元函数  来说，偏导数实际上就是指的导数，其定义如下：

几何意义：在点  处的切线（  关于  的变化率），如下图：

（2）多元函数情况（以二元函数为例，二元以上可类推）

对于二元函数 ，我们同样也可以定义：

①  关于  的偏导数可以定义如下：

几何意义：  关于  的变化率，即：  不变，  变化，  随之发生的变化量。

②  关于  的偏导数可以定义如下：

几何意义： 关于  的变化率，即：  不变，  变化，  随之发生的变化量。

2、方向导数

上述定义了多元函数的偏导数，但是偏导数只是描述了函数在某一个变量方向上的的变化率，但是对于多元函数来说，某一点上有无数个变化方向，那么我们如何衡量多元函数向任意一个方向变化的变化率呢？答案就是方向导数。

还是以二元函数  为例。如下图：

点  沿着  方向变化，变化量为向量  的长度。那么此时  的变化率就为：

令 ，则 ，。因此定义式  可以变换为：

那么我们如何计算方向导数呢？首先函数在某一点  处方向导数存在的充分条件是函数  在该点可微。根据可微的定义可知：

3、梯度

进一步引入向量内积的概念。向量  与  的内积，其基本定义为：

几何定义为：

因此，我们方向导数  式可以写成向量内积的形式，如下：

我们将上述  式中  称为 梯度向量 ，记作：，或者 。这里切记，梯度是个向量。向量  其实就是变化方向上的单位向量。

（二）为什么梯度方向（正、反）是函数值局部变化（上升、下降）最快的方向？

根据上述定义，方向导数其实就是梯度向量和变化方向上的单位向量的内积。即：

如下图：

根据内积几何定义可得：

当 ，即：  时，也就是说当变化方向与梯度同向时， 方向导数取正数最大值，因此，函数上升最快。

当 ，即：  时，也就是说当变化方向与梯度反向时， 方向导数取负数最小值，因此，函数下降最快。

（三）梯度下降法

1、基本原理

那么什么是梯度下降法呢？首先我们要明确，梯度下降法（Gradient Descent）并不是一个机器学习算法，而是一种基于搜索的最优化方法。通过梯度下降法，能够最小化一个损失函数。

以线性回归为例，我们的损失函数为：

我们的目标是最优化损失函数 ，即：找到一组参数  使得  最小。

由于损失函数  是关于  的函数，所以我们只要控制  朝着  下降最快（梯度负方向）的方向进行变化，就最终能找到  的局部最小值。

如何控制参数向量  向着  下降最快的方向进行变化呢？方法就是在梯度的负方向上进行逐步增量，用公式表示就是：

其中：

（1）负号表示“负方向”；

（2） 为学习率，是控制  向着  的局部最小值进行变化的步长，学习率  不宜过大（太大容易跑过局部最小值），也不宜太小（太小前进的速度太缓慢）；

（3） 为梯度向量，梯度向量控制变化方向。

注意： 一定要形成， 和梯度都是向量的概念，向量是。

我们以单一参数的损失函数  为例，以图形进行直观地展示，如下：

如上图，我们首先初始化一个参数向量 ，由于这里只是单参数，所以向量  就是一个标量（为了使向量的概念深深引入脑海，这里也将其当作一个特殊向量），梯度 。我们将初始化的参数向量  向着梯度的负方向进行更新，得到新的参数向量。即：  ，继续更新，只到损失函数  达到（局部）最小。

那我们什么时候停止迭代呢？因为在沿着梯度的方向进行迭代时，  逐渐接近局部最小值，相应地，梯度也会逐渐趋向0，因此前进的步子也会越来越小。根据这一特性，我们可以设计参数迭代的收敛条件如下：

(1) 

(2) 

2、学习率的选择

（1）学习率过小

如图，学习率太小，会导致收敛速度缓慢。

（2）学习率过大

如图，学习率太大，可能会越过最小值点，导致最终无法收敛。

在实践中，一般选择 。

3、全局最小值与局部最小值的问题

如上图，损失函数  是一个具有两个极值点的函数，当我们初始点设置在“初始点1”时，我们就会滑入一个局部最优解，而不能找到真正的全局最优解。如何解决这个问题呢？解决方法就是多次随机化初始点，以期望能找到类似“初始点2”这样的能够滑入全局最优解的点。

（四）补充：梯度上升法

二、模拟实现与原理可视化展示

（一）单一参数

1、分步验证

# 生成模拟数据
plot_theta = np.linspace(-2., 6., 50)
plot_theta

plot_J = (plot_theta - 2) ** 2 + 5
plot_J

# 绘图
plt.figure()
plt.plot(plot_theta, plot_J, color='steelblue')
plt.show()

# 定义损失函数
def Func_J(theta):
return (theta - 2) ** 2 + 5

# 定义求导（求梯度）函数
def grad(theta):
return 2 * (theta - 2)

# 以theta变化极小作为收敛条件
eta = 0.1 # 学习率
epsilon = 1e-8 # 收敛阈值
theta = 0.0 # 初始theta
i = 0
while True:
gradient = grad(theta)
theta_old = theta
theta = theta_old - eta * gradient
i += 1

if abs(theta - theta_old) < epsilon:
break

print('一共迭代了{}次'.format(i))
print('以theta变化极小作为收敛方式时得到的参数theta为：{}'.format(theta))
print('以theta变化极小作为收敛方式时得到的最小值为：{}'.format(Func_J(theta)))


# 以J(theta)变化极小作为收敛条件
eta = 0.1 # 学习率
epsilon = 1e-8 # 收敛阈值
theta = 0.0 # 初始theta
i = 0
while True:
gradient = grad(theta)
theta_old = theta
theta = theta_old - eta * gradient
i += 1
if abs(Func_J(theta) - Func_J(theta_old)) < epsilon:
break

print('一共迭代了{}次'.format(i))
print('以J变化极小作为收敛方式时得到的参数theta为：{}'.format(theta))
print('以J变化极小作为收敛方式时得到的最小值为：{}'.format(Func_J(theta)))


输出：

一共迭代了80次 

以theta变化极小作为收敛方式时得到的参数theta为：1.9999999646630588

以theta变化极小作为收敛方式时得到的最小值为：5.000000000000001

一共迭代了44次 

以J变化极小作为收敛方式时得到的参数theta为：1.999891109642585

以J变化极小作为收敛方式时得到的最小值为：5.00000001185711

# 绘图展示梯度下降过程

def Func_J(theta):
    return (theta - 2) ** 2 + 5

def grad(theta):
    return 2 * (theta - 2)

eta = 0.1
epsilon = 1e-8
theta = 0.0
theta_list = [theta]
while True:
    gradint = grad(theta)
    theta_old = theta
    theta = theta_old - eta * gradint  
    theta_list.append(theta)
    
    if abs(Func_J(theta) - Func_J(theta_old)) < epsilon:
        break

theta
Func_J(theta)
        
plt.figure()
plt.plot(plot_theta, plot_J, color='steelblue')
plt.plot(np.array(theta_list), Func_J(np.array(theta_list)), color='indianred', marker='>')
plt.show()
print('一共迭代了{}次'.format(len(theta_list)))  # 输出：一共迭代了45次

2、简单封装

class GradientDescent(object):

def __init__(self, eta=0.1, epsilon=1e-8, initial_theta=0.0, n_iters=1e4):
self.eta = eta
self.epsilon = epsilon
self.theta = initial_theta
self.theta_list = [initial_theta]

self.n_iters = n_iters
self.i_iter = 0
self.i_iter_list = []

def Func_J(self, theta):
return (theta - 2) ** 2 + 5

def grad(self, theta):
return 2 * (theta - 2)

def gradient_descent(self):

while self.i_iter < self.n_iters:
theta_old = self.theta
self.theta = self.theta - self.eta * self.grad(theta_old)
self.theta_list.append(self.theta)
self.i_iter_list.append(self.i_iter)

if abs(self.Func_J(self.theta) - self.Func_J(theta_old)) < self.epsilon:
break

self.i_iter += 1

return self

def plot_theta_history(self):
plot_theta = np.linspace(-2., 6., 50)
plot_J = (plot_theta - 2) ** 2 + 5

plt.figure()
plt.title('梯度下降过程图')
plt.plot(plot_theta, plot_J, color='steelblue')
plt.plot(np.array(self.theta_list), Func_J(np.array(self.theta_list)), color='indianred', marker='>')
plt.show()

print('一共迭代了{}次'.format(len(self.theta_list) - 1))

def plot_J_history(self):
plt.figure()
plt.title('损失函数变化图')
plt.plot(np.array(self.i_iter_list), Func_J(np.array(self.theta_list[1:])), color='steelblue')
plt.show()

def __repr__(self):
return 'GradientDescent()'


# 这是将参数调整可视化的一段代码，你可以得到图中一样的调整控件
%matplotlib inline
@interact(eta=(0.0, 2.0, 0.0001), epsilon=(1e-8,1,1e-9), n_iters=(1,1e4,1), initial_theta=(-2,6,0.01), continuous_update=False)
# eta为学习率（步长）；n_iters为最大迭代次数；init_theta为初始参数
def visualize(eta=0.1, epsilon=1e-8, initial_theta=0.0, n_iters=1e4):
    try:
        GD = GradientDescent(eta, epsilon, initial_theta, n_iters)
        GD.gradient_descent()
        GD.plot_theta_history()
        
        print('最终参数为：{}'.format(GD.theta))
        print('学习率为：{}'.format(eta))
        print('收敛阈值为：{}'.format(epsilon))
        print('初始参数为：{}'.format(initial_theta))
        print('可迭代最大次数为：{}'.format(n_iters))
        
        GD.plot_J_history()
        
    except OverflowError:
        print('(34, \'Result too large\')')

3、学习率的选择

# 学习率为0.1 （合适）
GD_1 = GradientDescent(eta=0.1, epsilon=1e-8, initial_theta=0.0)
GD_1.gradient_descent()
GD_1.plot_theta_history()
GD_1.theta
GD_1.plot_J_history()


输出：

一共迭代了44次

1.999891109642585

# 学习率为0.001 （太小，收敛缓慢）
GD_2 = GradientDescent(eta=0.001, epsilon=1e-8, initial_theta=4.0)
GD_2.gradient_descent()
GD_2.plot_theta_history()
GD_2.theta
GD_2.plot_J_history()

输出：

一共迭代了3569次

2.0015773637295693

# 学习率为1.1 （过大，不收敛）
GD_3 = GradientDescent(eta=1.1, epsilon=1e-8, initial_theta=0.0, n_iters=8)
GD_3.gradient_descent()
GD_3.plot_theta_history()
GD_3.theta
GD_3.plot_J_history()


输出：

一共迭代了8次

-6.599633920000013

（二）双参数

# 这是双参数下的一个梯度下降过程可视化图形
#梯度函数
def Func_J(x, y):
    return 2 * np.power(x, 2) + np.power(y, 2)

#梯度函数的导数
def grad_theta_1(theta):
    return 4 * theta
def grad_theta_2(theta):
    return 2 * theta

def gradient_descent(eta, n_iters, theta1, theta2, up, dirc):
    t1 = [theta1]
    t2 = [theta2]
    for i in range(n_iters):
        gradient = grad_theta_1(theta1)
        theta1 = theta1 - eta * gradient
        t1.append(theta1)
        gradient = grad_theta_2(theta2)
        theta2 = theta2 - eta * gradient
        t2.append(theta2)
        
    plt.figure(figsize=(10,10))  #设置画布大小
    # 生成虚拟数据
    x = np.linspace(-3,3,30)
    y = np.linspace(-3,3,30)
    
    # 转换成网格数据
    X, Y = np.meshgrid(x, y)
    Z = Func_J(X, Y)    
    ax = plt.axes(projection='3d')
    fig =plt.figure()
    ax.contour3D(X, Y, Func_J(X,Y), 50, cmap='binary') #等高线图
    
    ax.scatter3D(t1, t2, Func_J(t1,t2), c='r',marker = 'o')

    ax.view_init(up, dirc)
    
    return t1, t2

%matplotlib inline
@interact(eta=(0, 2, 0.0002),n_iters=(1,100,1),initial_theta1=(-3,3,0.1),initial_theta2=(-3,3,0.1),up=(-180,180,1),dirc=(-180,180,1),continuous_update=False)
#lr为学习率（步长） epoch为迭代次数   init_theta为初始参数的设置 up调整图片上下视角 dirc调整左右视角
def visualize_gradient_descent(eta=0.05,n_iters=10,initial_theta1=-2,initial_theta2=-3,up=45,dirc=100):
    gradient_descent(eta,n_iters,initial_theta1,initial_theta2,up,dirc)

三、线性回归中的梯度下降法

（一）基本原理

对于一元线性回归，其损失函数如下：

令 ，则上式可以写成下面的形式

进一步考虑到样本量的影响，将上式除以一个样本量 ，得到：

我们将  式作为我们最终的损失函数。因此，我们的目标就是最优化  式。

对损失函数求关于  的梯度：

令：

则：

求出梯度之后，我们开始进行迭代，过程如下：

step1： 初始化一组参数 

Step2：

开始迭代

只到：

则停止，并返回 

（二）分步实现

1、使用 for 循环方式实现梯度下降

# 生成数据
np.random.seed(666)
x_1 = 2 * np.random.random(size=100)
y = x_1 * 3. + 4. + np.random.normal(size=100)

# 绘图
plt.figure()
plt.scatter(x_1, y, color='steelblue')
plt.show()


# 定义损失函数
def Func_J1(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta))**2) / len(X_b)
    except:
        return float('inf')

# 定义求梯度函数（用for的方法）
def grad1(theta, X_b, y):
grad_j = np.empty(len(theta))
grad_j[0] = np.sum(X_b.dot(theta) - y)
for i in range(1, len(theta)):
grad_j[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
return grad_j * 2 / len(X_b)


# 定义梯度下降迭代函数
def gradient_descent1(X_b, y, initial_theta, eta=0.1, n_iters=1e4, epsilon=1e-8):
    
    theta = initial_theta
    theta_list = [initial_theta]
    i_iter = 0
    
    while i_iter < n_iters:
        theta_old = theta
        theta = theta - eta * grad1(theta_old, X_b, y)
        theta_list.append(theta)

        if abs(Func_J1(theta, X_b, y) - Func_J1(theta_old, X_b, y)) < epsilon:
            break

        i_iter += 1
        
    return theta, np.array(theta_list)

x_0 = np.ones((len(x_1), 1))
X_b = np.hstack([x_0, x_1.reshape(-1,1)])
initial_theta = np.zeros(X_b.shape[1])

theta, theta_list = gradient_descent1(X_b, y, initial_theta, eta=0.01, n_iters=1e4, epsilon=1e-8)
theta # 输出：array([4.02145786, 3.00706277])


2、使用向量方式实现梯度下降

# 生成数据
np.random.seed(666)
x_1 = 2 * np.random.random(size=100)
y = x_1 * 3. + 4. + np.random.normal(size=100)

# 绘图
plt.figure()
plt.scatter(x_1, y, color='steelblue')
plt.show()

def Func_J2(theta, X_b, y):
try:
return (y - X_b.dot(theta)).T.dot(y - X_b.dot(theta))/ len(X_b)
except:
return float('inf')


# 用向量方法
def grad2(theta, X_b, y):
    return (X_b).T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

# 定义梯度下降迭代函数
def gradient_descent2(X_b, y, initial_theta, eta=0.1, n_iters=1e4, epsilon=1e-8):

theta = initial_theta
theta_list = [initial_theta]
i_iter = 0

while i_iter < n_iters:
theta_old = theta
theta = theta - eta * grad2(theta_old, X_b, y)
theta_list.append(theta)

if abs(Func_J2(theta, X_b, y) - Func_J2(theta_old, X_b, y)) < epsilon:
break

i_iter += 1

return theta, np.array(theta_list)


x_0 = np.ones((len(x_1), 1))
X_b = np.hstack([x_0, x_1.reshape(-1,1)])
initial_theta = np.zeros(X_b.shape[1])

theta, theta_list = gradient_descent2(X_b, y, initial_theta, eta=0.01, n_iters=1e4, epsilon=1e-8)
theta  # 输出：array([4.02145786, 3.00706277])

（三）封装

import numpy as np

class MY_LinearRegression:

def __init__(self):
'''初始化'''
self.intercept_ = None
self.coef_ = None
self._theta =None

def fit_normal(self, X_train, y_train):
'''使用正规方程方法训练'''
assert X_train.shape[0] == y_train.shape[0], 'X_train 和 y_train 的行数要一致'

X_0 = np.ones(shape=(len(X_train), 1))
X_b = np.hstack([X_0, X_train])

self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

self.intercept_ = self._theta[0]
self.coef_ = self._theta[1:]

return self

def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4, epsilon=1e-8, method='f'):
'使用梯度下降法进行训练'
assert X_train.shape[0] == y_train.shape[0], 'X_train 和 y_train 的行数要一致'

def Func_J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta))**2) / len(X_b)
except:
return float('inf')

def grad(theta, X_b, y):
# for循环
if method == 'f':
grad_j = np.empty(len(theta))
grad_j[0] = np.sum(X_b.dot(theta) - y)
for i in range(1, len(theta)):
grad_j[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
return grad_j * 2 / len(X_b)
# 向量法
if method == 'm':
return (X_b).T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

# 定义梯度下降迭代函数
def gradient_descent(X_b, y, initial_theta, eta, n_iters, epsilon):

theta = initial_theta
i_iter = 0

while i_iter < n_iters:
theta_old = theta
theta = theta - eta * grad(theta_old, X_b, y)

if abs(Func_J(theta, X_b, y) - Func_J(theta_old, X_b, y)) < epsilon:
break

i_iter += 1

return theta

x_0 = np.ones((len(X_train), 1))

if len(X_train.shape) == 1:
X_train = X_train.reshape(-1,1)

X_b = np.hstack([x_0, X_train])
initial_theta = np.zeros(X_b.shape[1])

self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters, epsilon)
self.intercept_ = self._theta[0]
self.coef_ = self._theta[1:]

return self

def predict(self, X_predict):
'''预测'''
assert self.intercept_ is not None and self.coef_ is not None, '请先进行训练'
assert X_predict.shape[1] == len(self.coef_), '特证个数必须与参数（非截距）个数相等'

X_0 = np.ones(shape=(len(X_predict), 1))
X_b = np.hstack([X_0, X_predict])

return X_b.dot(self._theta) # y = X_b.dot(𝜃)

def r2_score(self, X_test, y_test):
'''模型评价R2'''
y_test_predict = self.predict(X_test)
mean_squared_error = np.sum((y_test - y_test_predict) ** 2) / len(y_test)
return 1 - mean_squared_error / np.var(y_test)

def __repr__(self):
return 'MY_LinearRegression()'


# 导入波士顿房价数据
from sklearn import datasets

boston = datasets.load_boston()

# print(boston.DESCR)
# boston.feature_names

X = boston.data
y = boston.target

# 去除天花板值
X = X[y < 50.0]
y = y[y < 50.0]

X.shape  # 输出：(490, 13)
y.shape  # 输出：(490,)

# 切分数据集
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=666)

X_train.shape # 输出：(392, 13)
X_test.shape # 输出：(98, 13)
y_train.shape # 输出：(392,)
y_test.shape # 输出：(98,)


1、正规方程训练

reg_n = MY_LinearRegression()
reg_n.fit_normal(X_train, y_train)
reg_n._theta  
reg_n.r2_score(X_test, y_test)  # 输出：0.8129794056212729

2、梯度下降法（for循环方法训练）

reg_gd1 = MY_LinearRegression()
reg_gd1.fit_gd(X_train, y_train, eta=0.000001, method='f', n_iters=1e6)
reg_gd1._theta
reg_gd1.r2_score(X_test, y_test) # 输出：0.7542932581943915


3、梯度下降法（向量法训练）

reg_gd2 = MY_LinearRegression()
reg_gd2.fit_gd(X_train, y_train, eta=0.000001, method='m', n_iters=1e6)
reg_gd2._theta
reg_gd2.r2_score(X_test, y_test)  # 输出：0.7542932581943915

四、梯度下降法中数据归一化的必要性

通过上面的实际训练可以看出，梯度下降法相对于正规方程法，模型训练时间较长，原因之一就是因为我们在训练之前没有进行数据归一化。

为什么不进行数据归一化会导致梯度下降法的训练时间时间较长呢？我们以两参数模型的损失函数  为例。

首先我们必须知道，梯度下降的方向是与损失函数的等高线 （ 为常数）的切线垂直的。为什么垂直？我们做如下证明：

已知曲面 ，其等高线方程可以表示为：（ 为常数）。我们对  两边同时进行微分，得：

如上推导，我们可以看出梯度向量  与切向量  的内积为0，因此两个向量垂直。又由于梯度下降的方向是梯度的反方向，因此梯度下降的方向也是与等高线切线垂直的。

知道梯度下降的方向与等高线切线垂直后，我们可以绘制两种情况下的损失函数等高线图，如下：

图一：

图二：

图一表示的是特征  与特征  的取值范围相近的情况，在这种情况下，等高线比较接近正圆。因此梯度下降的时候直奔最小值点去，因此迭代次数较少。

图二表示的是特征  与特征  的取值范围量级相差较大的情况，在这种情况下，等高线比较“消瘦”。梯度下降的时候会在下降过程中反复横跳，因此导致迭代次数较多。

由此上述可见，我们使用梯度下降法进行寻找最优解的时候，进行数据归一化是很有必要的！

# 对特征数据进行标准化
from sklearn.preprocessing import StandardScaler

StandardScaler = StandardScaler()
StandardScaler.fit(X_train)
X_train_standard = StandardScaler.transform(X_train)


# 训练
reg_gd3 = MY_LinearRegression()
reg_gd3.fit_gd(X_train_standard, y_train, method='m')

X_test_standard = StandardScaler.transform(X_test)
reg_gd3._theta
reg_gd3.r2_score(X_test_standard, y_test) # 输出：0.8129873310487505


通过运行时间和  来看，在进行数据归一化后，整个训练时间缩短了很多（没进行数据归于化之前需要48s，归一化之后直用了238ms），并且从  来看，模型的训练效果也有一定的提升。

五、梯度下降法相对于正规方程法的优势

当特征维度  比较大时，正规方程的时间复杂度较大，但是梯度下降法训练速度就很快。以下以一个虚拟数据进行二者的比较：

# 生成虚拟数据
n = 10000
m = 1000

X = np.random.normal(size=(m, n))
true_theta = np.random.uniform(0.0, 50.0, size=n+1)
y = X.dot(true_theta[1:]) + true_theta[0] + np.random.normal(0.0, 10.0,size=m)

# 正规方程方法训练 --花了17s，不同计算机时间不一样，仅作参考
reg_n = MY_LinearRegression()
reg_n.fit_normal(X, y)


# 梯度下降法训练  --花了1.7s，不同计算机时间不一样，仅作参考
reg_gd = MY_LinearRegression()
reg_gd.fit_gd(X, y, method='m')

通过以上测试，梯度下降法的训练速度明显快于正规方程

六、随机“梯度”下降法

（一）基本原理

1、“梯度”的计算

前述梯度下降法中，我们的梯度表示为：

这种梯度下降法在每次计算梯度时间，将所有的样本都进行了使用，我们称之为批量梯度下降法。批量梯度下降法在样本量很大时，使用矩阵计算梯度时十分耗时，因此我们考虑在每次计算梯度时，只考虑其中一个样本进行计算。即：

迭代式为：

这种方法我们称之为随机梯度下降法。标题之所以打上引号，其实是因为这里的计算方式，使得我们参数迭代的方向并不是沿着梯度的反方向，而是随机方向，因此并不能称为严格意义的“梯度下降”。

2、学习率的计算

正如以上讨论，随机梯度下降法下，我们参数迭代的方向并不是严格的梯度反方向，因此并不能保证每次迭代都能使  减小，也不能保证每次迭代都能使损失函数  减小。因此可能存在当我们迭代到损失函数局部最小值附近时，由于迭代步长太大了，结果慢慢地又远离了最优解。因此我们需要人为地定义一个递减的的学习率，使得迭代量逐渐递减（批量梯度下降法下，每次的迭代量是天然递减的）。

因此我们定义学习率为：

（二）分步实现

# 生成模拟数据
m = 100000
np.random.seed(666)
x_1 = np.random.normal(size=m)
y = 4.*x_1 + 3. + np.random.normal(0, 3, size=m)


1、使用批量梯度下降法训练

def Func_J(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta))**2) / len(X_b)
    except:
        return float('inf')

def Bgrad(theta, X_b, y):
    return (X_b).T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

# 定义梯度下降迭代函数
def B_gradient_descent(X_b, y, initial_theta, eta=0.01, n_iters=1e4, epsilon=1e-8):

    theta = initial_theta
    i_iter = 0
    i_iter_list = []
    J_list = []

    while i_iter < n_iters:
        theta_old = theta
        theta = theta - eta * Bgrad(theta_old, X_b, y)
        i_iter_list.append(i_iter)
        J = Func_J(theta, X_b, y)
        J_list.append(J)
        if abs(Func_J(theta, X_b, y) - Func_J(theta_old, X_b, y)) < epsilon:
            break

        i_iter += 1

    return theta, i_iter_list, J_list

x_0 = np.ones((len(x_1), 1))
x_1 = x_1.reshape(-1,1)
X_b = np.hstack([x_0, x_1])
initial_theta = np.zeros(X_b.shape[1])
theta, i_iter_list, J_list = B_gradient_descent(X_b, y, initial_theta)
theta # 输出：array([3.00545329, 4.00684163])


plt.figure()
plt.title('批量梯度下降法损失函数变化图示')
plt.plot(i_iter_list, J_list, color='steelblue')
plt.show()

2、使用随机梯度下降法训练

def Func_J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta))**2) / len(X_b)
except:
return float('inf')

def Sgrad(theta, X_b_i, y_i):
return (X_b_i).T.dot(X_b_i.dot(theta) - y_i) * 2

def S_gradient_descent(X_b, y, initial_theta, n_iters=100, t0=5, t1=50):

theta = initial_theta
i_iter_list = []
J_list = []

# 定义学习率函数
def learning_rate(i_iter):
return t0 / (i_iter + t1)

for i_iter in range(n_iters):
theta_old = theta
eta = learning_rate(i_iter)
rand_i = np.random.randint(len(X_b)) # 建立一个随机索引
grad = Sgrad(theta_old, X_b[rand_i], y[rand_i])
theta = theta - eta * grad

i_iter_list.append(i_iter)
J = Func_J(theta, X_b, y)
J_list.append(J)

return theta, i_iter_list, J_list


x_0 = np.ones((len(x_1), 1))
x_1 = x_1.reshape(-1,1)
X_b = np.hstack([x_0, x_1])
initial_theta = np.zeros(X_b.shape[1])
n_iters = len(X_b)//3  # // 取整

theta, i_iter_list, J_list = S_gradient_descent(X_b, y, initial_theta, n_iters)
theta  # 输出：array([2.91924481, 4.0059984 ])

plt.figure()
plt.title('随机梯度下降法损失函数变化图示')
plt.plot(i_iter_list, J_list, color='steelblue')
plt.show()


比较“批量梯度下降”和“随机梯度下降”的损失函数变化曲线可以看出，“随机梯度下降”的损失函数会出现反复横跳的现象，这也是我们要人为定义一个递减的学习率的重要原因。

（三）封装

import numpy as np

# 使用for循环
class MY_LinearRegression:
    
    def __init__(self):
        '''初始化'''
        self.intercept_ = None
        self.coef_ = None
        self._theta =None
     
    
    # 正规方程法    
    def fit_normal(self, X_train, y_train):
        '''使用正规方程方法训练'''
        assert X_train.shape[0] == y_train.shape[0], 'X_train 和 y_train 的行数要一致'
        
        X_0 = np.ones(shape=(len(X_train), 1))
        X_b = np.hstack([X_0, X_train])
        
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
        
        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]
        
        return self
    
    
    # 批量梯度下降法
    def fit_bgd(self, X_train, y_train, eta=0.01, n_iters=1e4, epsilon=1e-8, method='f'):
        '使用批量梯度下降法进行训练'
        assert X_train.shape[0] == y_train.shape[0], 'X_train 和 y_train 的行数要一致'
        
        def Func_J(theta, X_b, y):
            try:
                return np.sum((y - X_b.dot(theta))**2) / len(X_b)
            except:
                return float('inf')
                
        def b_grad(theta, X_b, y):
            # for循环
            if method == 'f':
                grad_j = np.empty(len(theta))
                grad_j[0] = np.sum(X_b.dot(theta) - y)
                for i in range(1, len(theta)):
                    grad_j[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
                return grad_j * 2 / len(X_b)
            # 向量法
            if method == 'm':
                    return (X_b).T.dot(X_b.dot(theta) - y) * 2 / len(X_b)
                
        # 定义梯度下降迭代函数
        def b_gradient_descent(X_b, y, initial_theta, eta, n_iters, epsilon):

            theta = initial_theta
            i_iter = 0

            while i_iter < n_iters:
                theta_old = theta
                theta = theta - eta * b_grad(theta_old, X_b, y)

                if abs(Func_J(theta, X_b, y) - Func_J(theta_old, X_b, y)) < epsilon:
                    break

                i_iter += 1

            return theta
        
        x_0 = np.ones((len(X_train), 1))
        
        if len(X_train.shape) == 1:
            X_train = X_train.reshape(-1,1)
            
        X_b = np.hstack([x_0, X_train])
        initial_theta = np.zeros(X_b.shape[1])

        self._theta = b_gradient_descent(X_b, y_train, initial_theta, eta, n_iters, epsilon)
        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]
        
        return self
    
    
    # 随机梯度下降法
    def fit_sgd(self, X_train, y_train, n_iters=5, epsilon=1e-8, t0=5, t1=50):
        '使用随机梯度下降法进行训练'
        assert X_train.shape[0] == y_train.shape[0], 'X_train 和 y_train 的行数要一致'
        
        def Func_J(theta, X_b, y):
            try:
                return np.sum((y - X_b.dot(theta))**2) / len(X_b)
            except:
                return float('inf')
                
        def s_grad(theta, X_b_i, y_i):
            return (X_b_i).T.dot(X_b_i.dot(theta) - y_i) * 2
                
        # 定义梯度下降迭代函数
        def s_gradient_descent(X_b, y, initial_theta, n_iters, t0, t1):

            theta = initial_theta

            # 定义学习率函数
            def learning_rate(i_iter):
                return t0 / (i_iter + t1)
            
            m = len(X_b)
            
            for i in range(n_iters):  # n_iters表示看几遍样本
                # 将样本打乱顺序
                indexs = np.random.permutation(m) 
                X_b_new = X_b[indexs,:]
                y_new = y[indexs]
                
                # 遍历当前样本集
                for i_iter in range(m):
                    theta_old = theta
                    eta = learning_rate(i * m + i_iter)
                    grad = s_grad(theta_old, X_b[i_iter], y[i_iter])
                    theta = theta - eta * grad

            return theta
        
        x_0 = np.ones((len(X_train), 1))
        
        if len(X_train.shape) == 1:
            X_train = X_train.reshape(-1,1)
            
        X_b = np.hstack([x_0, X_train])
        initial_theta = np.zeros(X_b.shape[1])

        self._theta = s_gradient_descent(X_b, y_train, initial_theta, n_iters, t0, t1)
        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]
        
        return self
    

    def predict(self, X_predict):
        '''预测'''
        assert self.intercept_ is not None and self.coef_ is not None, '请先进行训练'
        assert X_predict.shape[1] == len(self.coef_), '特证个数必须与参数（非截距）个数相等'
        
        X_0 = np.ones(shape=(len(X_predict), 1))
        X_b = np.hstack([X_0, X_predict])
        
        return X_b.dot(self._theta) # y = X_b.dot(𝜃)
    
    def r2_score(self, X_test, y_test):
        '''模型评价R2'''
        y_test_predict = self.predict(X_test)
        mean_squared_error = np.sum((y_test - y_test_predict) ** 2) / len(y_test)
        return 1 - mean_squared_error / np.var(y_test)
    
    def __repr__(self):
        return 'MY_LinearRegression()'

1、模拟数据验证

# 生成模拟数据
m = 100000
np.random.seed(666)
X_train = np.random.normal(size=m)
y_train = 4.*X_train + 3. + np.random.normal(0, 3, size=m)


# 批量梯度下降
reg_bgd = MY_LinearRegression()
reg_bgd.fit_bgd(X_train, y_train, method='m')
reg_bgd._theta  # 输出：array([3.00545329, 4.00684163])

# 随机梯度下降
reg_sgd = MY_LinearRegression()
reg_sgd.fit_sgd(X_train, y_train, n_iters=1)
reg_sgd._theta # 输出：array([3.01080423, 4.02101668])


2、波士顿房价数据

# 导入波士顿房价数据
from sklearn import datasets

boston = datasets.load_boston()

# print(boston.DESCR)
# boston.feature_names

X = boston.data
y = boston.target

# 去除天花板值
X = X[y < 50.0]
y = y[y < 50.0]

X.shape  # 输出：(490, 13)
y.shape  # 输出：(490,)

# 切分数据集
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=666)

X_train.shape  # 输出：(392, 13)
X_test.shape  # 输出：(98, 13)
y_train.shape  # 输出：(392,)
y_test.shape  # 输出：(98,)

# 数据归一化
from sklearn.preprocessing import StandardScaler

StandardScaler = StandardScaler()
StandardScaler.fit(X_train)
X_train_standard = StandardScaler.transform(X_train)
X_test_standard = StandardScaler.transform(X_test)

# 批量梯度下降
reg_bgd = MY_LinearRegression()
reg_bgd.fit_bgd(X_train_standard, y_train, eta=0.000001, method='m', n_iters=1e6)
reg_bgd._theta
reg_bgd.r2_score(X_test_standard, y_test) # 输出：0.6656195606376432


# 随机梯度下降
reg_sgd = MY_LinearRegression()
reg_sgd.fit_sgd(X_train_standard, y_train, n_iters=100)
reg_sgd._theta
reg_sgd.r2_score(X_test_standard, y_test)  # 输出：0.813295063363921

（四）scikit - learn 中的随机梯度下降法

from sklearn.linear_model import SGDRegressor

# sklearn中的随机梯度下降 浏览样本遍数 n_iter_no_change 选择默认
s_reg_sgd = SGDRegressor()
s_reg_sgd.fit(X_train_standard, y_train)
s_reg_sgd.score(X_test_standard, y_test) # 输出：0.8129254242161921


# sklearn中的随机梯度下降 浏览样本遍数 n_iter_no_change=100
s_reg_sgd_100 = SGDRegressor(n_iter_no_change=100)
s_reg_sgd_100.fit(X_train_standard, y_train)
s_reg_sgd_100.score(X_test_standard, y_test)  # 输出：0.8129364925974142

七、梯度的调试

前述我们都是通过数学推导的方式得出计算梯度的公式，但是对于一些复杂的损失函数来说，推导梯度计算公式并不是那么容易的事情，很可能会出现推导错误的情况。那么我们怎样验证推导的计算方法的准确性呢？

首先，我们引入一个近似表示函数在某点的导数的方法，如图：

如图，我们可以使用  近似表示损失函数  在点  处的导数。

上述方法对于高维同样适用。

如下

令：

则：

同理可以推出  关于其他参数的偏导数。因此梯度可以表示为：

def Func_J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta))**2) / len(X_b)
except:
return float('inf')

def b_grad_math(theta, X_b, y):
return (X_b).T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

def b_grad_debug(theta, X_b, y, epsilon=0.001):
grad = 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
grad[i] = (Func_J(theta_1, X_b, y) - Func_J(theta_2, X_b, y)) / (2 * epsilon)

return grad


# 定义梯度下降迭代函数
def b_gradient_descent(grad_method, X_b, y, initial_theta, eta=0.01, n_iters=1e4, epsilon=1e-8):

theta = initial_theta
i_iter = 0

while i_iter < n_iters:
theta_old = theta
theta = theta - eta * grad_method(theta_old, X_b, y)

if abs(Func_J(theta, X_b, y) - Func_J(theta_old, X_b, y)) < epsilon:
break

i_iter += 1

return theta


# 生成虚拟数据
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)

# 使用推导方式计算的梯度
initial_theta = np.zeros(X_b.shape[1])
theta = b_gradient_descent(b_grad_math, X_b, y, initial_theta)
theta


输出：

array([ 1.1251597 ,  2.05312521,  2.91522497,  4.11895968,  5.05002117, 5.90494046,  6.97383745,  8.00088367,  8.86213468,  9.98608331, 10.90529198])

# 使用近似的方式计算梯度
initial_theta = np.zeros(X_b.shape[1])
theta = b_gradient_descent(b_grad_debug, X_b, y, initial_theta)
theta

输出：

array([ 1.1251597 ,  2.05312521,  2.91522497,  4.11895968,  5.05002117, 5.90494046,  6.97383745,  8.00088367,  8.86213468,  9.98608331, 10.90529198])

如上所示，近似方式计算的梯度训练效果也比较好，但是其时间复杂度较高，一次一般不采用这种方法。但是我们可以将其应用到梯度验证上。我们在处理一些比较复杂的损失函数的时候，可以先用近似的方法计算梯度，然后进行训练。再对我们推导出的梯度计算方法进行验证。

八、总结

（一）梯度下降法分类

梯度下降法主要包括：

1、批量梯度下降法 （Batch Gradient Descent）

• 该种梯度下降法比较稳定，一定能够走向局部最小值，但是时间复杂度较大。

2、随机梯度下降法 （Stochastic Gradient Descent）

• 该种梯度下降法由于每次并不是朝着梯度的反方向进行前进，所以不一定能够走到局部最小值，因此不是很稳定，但是时间复杂度较小。

3、小批量梯度下降法 （Mini-Batch Gradient Descent）

• 该种梯度下降法是结合了以上两种方法，每次计算“梯度”时考虑k个样本。降低了两者单独呈现的缺点。

（二）随机的力量

1、跳出局部最优解

2、更快的运行速度

3、随机搜索；随机森林