机器学习系列（十四）——梯度下降法Gradient Descen

梯度下降法

梯度下降法本身不是一个机器学习算法，而是一种基于搜索的最优化方法，它是人工智能领域非常重要的一种方法。梯度下降的作用是最小化一个损失函数。与之对应的，如果是效用函数，则是梯度上升法，它的作用是最大化一个效用函数。
为什么使用梯度下降法呢？因为并非所有机器学习算法都像线性回归那样能求得数学上的最优解（表达式），对于无法直接求得最优解的，则使用梯度下降法来最小化损失函数。
如下图所示：

Grad_Desc

若损失值当前状态对应图中圆点，依照梯度下降法，为了求得最小的损失值，它会沿着梯度的反方向移动。当前待优化变量，梯度为，有了梯度，还需要有参数变化的幅度，称为学习率，设为，于是参数更新公式为：

整个梯度下降法的过程就是不断更新参数寻找最小损失值的过程，学习率是一个超参数，它的取值影响获得最优解的过程。学习率越大参数一次更新变化幅度会越大反之则越小，更新幅度过大，目标函数在最优值附近徘徊，永远找不到极值点，幅度过小又会增加训练轮数增加计算开销影响训练效率。训练不希望参数更新到比较坏的点上，也不希望训练开销太大。
不过梯度下降算法求解的结果并不一定是全局最优解，比如函数如果有两个极小值，由于学习幅度小，参数可能在一个极小值位置徘徊。由此可见参数的最初位置会影响优化过程和结果，因此初始点也是一个超参数。不过当损失函数在一段区域二阶导数恒大于0时，可以保证求解得到全局最优解（凸规划）。

---

梯度下降法在线性回归中的应用

上一篇文章求解的线性回归模型的解是正规方程解，虽然表达式简单但是计算开销大，这里将用梯度下降法来求解线性回归模型。在此之前先用一个简单的例子模拟一下梯度下降法。

模拟梯度下降

import numpy as np
import matplotlib.pyplot as plt
plot_x = np.linspace(-1,6,141)
plot_y = (plot_x-2.5)**2-1
plt.plot(plot_x,plot_y)
plt.show()

f(x)

用梯度下降法求图示二次函数的最小值：

def dJ(theta):
    '''损失函数的导数'''
    return 2*(theta-2.5)

def J(theta):
    '''损失值计算'''
    try:
        return (theta-2.5)**2-1
    except:
        return float('inf')

theta = 0.0
epsilon = 1e-5
eta = 0.1
while True:
    gradient = dJ(theta)
    last_theta = theta
    theta = theta - eta * gradient#参数更新
    
    if (abs(J(theta)-J(last_theta))<epsilon):
        break
        
print('best theta:',theta)
print('best lost:',J(theta))

result

下面将求解过程可视化：

'''梯度下降过程可视化'''
theta = 0.0
theta_history = [theta]
while True:
    gradient = dJ(theta)
    last_theta = theta
    theta = theta - eta * gradient#参数更新
    theta_history.append(theta)
    
    if (abs(J(theta)-J(last_theta))<epsilon):
        break

plt.plot(plot_x,J(plot_x))
plt.plot(np.array(theta_history),J(np.array(theta_history)),color='r',marker='+')
plt.show()

梯度下降可视化

theta_history中保存了梯度下降过程的theta，查看它即可查看迭代次数：

迭代次数

为了调节theta和eta这两个超参数的方便，封装上述算法：

'''调节参数theta和eta，为了方便，对上述计算封装'''
def gradient_descent(initial_theta, eta, n_iters = 1e4, epsilon = 1e-5):
    theta = initial_theta
    theta_history.append(initial_theta)
    i_iter = 0
    while i_iter < n_iters:
        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

def plot_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()

接下来可以测试不同的theta初始值和不同的eta对迭代的影响：

测试

自行调整超参数，测试更多的结果，这里不再赘述。

多元线性回归中的梯度下降

可以由最优化理论证明，线性回归法的损失函数（这里取作MSE）有唯一的最优解。由于多元线性回归中有多个变量，所以求得的梯度将是一个向量：

求解过程的参数更新公式可以结合多元函数微分学由单变量情况推广过来。

'''模拟数据'''
import numpy as np
import matplotlib.pyplot as plt

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)#100行1列
plt.scatter(x,y)
plt.show()

模拟数据

为了简单并阐明原理，模拟数据的回归方程为：y=3x+4，加了微小的正态噪声。用梯度下降训练该模型：

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):
    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)

def gradient_descent(X_b, y, initial_theta, eta, n_iters = 1e4, epsilon = 1e-5):
    theta = initial_theta
    i_iter = 0
    while i_iter < n_iters:
        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
    return theta

X_b = np.hstack([np.ones((len(X),1)),X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.1
theta = gradient_descent(X_b, y, initial_theta, eta)

求得结果如下：

线性回归GD求解结果

将算法封装进play_Ml的LinearRegression.py，在线性回归中已有fit_normal（正规解）方法，这里添加fit_gd方法，即梯度下降训练过程：

    def fit_gd(self,X_train,y_train,eta=0.1, n_iters=1e4):
        assert X_train.shape[0] == y_train.shape[0],"must be the same!"

        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):
            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)
        def gradient_descent(X_b, y, initial_theta, eta, n_iters = 1e4, epsilon = 1e-5):
            theta = initial_theta
            i_iter = 0
            while i_iter < n_iters:
                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
            return theta
        X_b = np.hstack([np.ones((len(X_train),1)),X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b,y_train,initial_theta,eta,n_iters)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

使用我们的算法：

模拟sklearn梯度下降

可见得到了同样的结果。