机器学习第五周 最优化方法：梯度下降

机器学习第五周 最优化方法：梯度下降

一、学习目标

• 梯度下降：实现梯度下降、线性回归中的梯度下降
• 随机梯度下降：相关代码及调用

二、学习内容

《还不了解梯度下降法？看完这篇就懂了！》

《手动实现梯度下降（可视化）》

《线性回归中的梯度下降》

《速度更快的随机梯度下降法》

《梯度下降番外：非常有用的调试方式及总结》

三、开动

概括总结：

机器学习于千万里寻找那一个函数f(x)，对这个函数进行优化，而这个函数可以进行预测、分类、生成等工作。

对函数的优化，也就是取得最优的值，怎么让损失函数最小。下面将要介绍的这种方法无疑是最核心、最广泛的一种方法。梯度下降

简单线性回归，计算用最小二乘法推导出参数计算公式；多元线性回归也是推导出损失函数。除了能进行推导，有没有一种通用的寻找最优化的过程方法，而这种方法对大多数函数都适用？

梯度下降： gradient descent GD.这是一种基于搜索的最优化方法，而不是一个机器学习算法。其作用是用来对原始模型的损失函数进行优化，以便寻找到最优的参数，使得损失函数的值最小。

最简单的我们可以想到选择1000个随机值，得到1000个结果，选取损失值最小的作为最优的参数值。这是一种方法。有没有一种方式通过每次取得的结果去更新参数，而且能够大幅降低计算次数。

通过参数与损失值的导数，找到朝什么方向运动，速度快慢，安全高效搜索最优值。

从高等数据微积分知识，梯度，可以用多元函数的导数，即分别对每个变量进行微分，然后用逗号分隔开。
$$▽L=(\frac{\partial L}{\partial a},\frac{\partial L}{\partial b})$$
梯度是向量。

举例：二元函数
$$f(x_1,x_2)=x_1^2+x_1{x}_2-3x_2$$

$$▽f=(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2})=(2x_1+x_2,x_1-3)$$

某个点带进去就可以求得对应点的梯度。梯度是函数的微分，代表函数在给定的某个点的切线或者变化上升或者下降最快的方向。

利用自然界中模拟某个过程抽象而来的算法，如遗传算法、下山法，模拟退火法等。
$$针对损失函数L，假设选取的初始值(a_0,b_0),现在该点稍微一动一点点，得到(a_1,b_1),根据泰勒展开式$$

$$f(x_1,x_2,...,x_n)=f(a_1,...,a_n)+\sum _{i=1}^n\frac{\partial f(a_1,a_2,...,a_n)}{\partial x_i}(x_i-a_i)+o(\sum (|x_i-a_i|))$$

$$△L=L(a_1,b_1)-L(a_0,b_0)\approx \frac{\partial L}{\partial a}△a+\frac{\partial L}{\partial b}△b$$

$$△a=-\eta \frac{\partial L}{\partial a},△b=-\eta \frac{\partial L}{\partial b},\eta >0$$

$$(a_{k+1},b_{k+1})=(a_{k+1}-\eta \frac{\partial L}{\partial a},b_{k+1}-\eta \frac{\partial L}{\partial b})$$

在搜索过程中，容易陷入局部最优解。一种解决方案是初始化多个参数及，即多组a0,b0,分别进行最优解的梯度下降算法，指导函数收敛于某个值；最后从其中找到最小值。该种方法不一定能找到全局最优解，但起码能找到较好的点。

注：这里的泰勒公式需要在注意一下。

from scipy.misc import  derivative
def f(x):
    return x**3 + x**2
derivative(f, 1.0, dx =1e-6)
#scipy.misc.derivative (func, x0 , dx=1.0 ,n=1 ,args(),order==3)
#func求导的函数，x0要求导的点，dx间距，n阶导数，order 奇数

import sympy  as sy
x = sy.Symbol('x')
func = 1/(1 + x**2)
type(x)
func
diff(func,x)
diff(func,x).subs(x,3)
diff(func,x).subs(x,3).evalf()

import numpy  as  np
import matplotlib.pyplot as  plt
from scipy.misc import derivative

#定义损失函数
def lossFunction(x):
    return (x-2.5)**2-1
plot_x = np.linspace(-1,6,141)
plot_y = lossFunction(plot_x)
plt.plot(plot_x, plot_y)
plt.show()

#定义导数
def dLF(theta):
    return  derivative(lossFunction, theta , dx =1e-6)

theta = 0.0 
eta = 0.1
epsilon =1e-6
while True:
    #沿着梯度下降
    gradient = dLF(theta)
    last_theta = theta
    #下一步走到
    theta = theta - eta * gradient
    #判断是否达到最小值
    if(abs(lossFunction(theta) - lossFunction(last_theta)) < epsilon):
        break
theta , lossFunction(theta)
#相当于逼近与一个最小值，因为最小值处导数为0，两次迭代的函数值差值应该为0   

def gradient_descent(initial_theta , eta , epsilon):
    theta = initial_theta
    #记录中间过程最佳,[]记录中间寻找的过程....
    history_theta.append(theta)
    while True:
        last_theta = theta 
        gradient = dLF(theta)
        theta = theta - eta * gradient
        history_theta.append(theta)
        if (abs(lossFunction(theta)-lossFunction(last_theta))<epsilon):
            break
    print('执行次数：%f',len(history_theta))
def plot_theta_history():
        plt.plot(plot_x,plot_y)
        plt.plot(np.array(history_theta),lossFunction(np.array(history_theta)),color='red',marker='+')
        plt.show()

#测试输出将结果
eta = 0.1
history_theta = []
gradient_descent(1.0 , eta , 1e-6)
plot_theta_history()

这是一种往最优值方向寻找的过程，当eta大小调整后，会影响收敛速度，甚至是导致不收敛。还有一种是不收敛搜索，最终趋近INF。

eta=0.05时，次数及结果：

所以要对其进行改进，设置搜索次数

优化修改搜索次数：

def gradient_descent2(initial_theta,eta,epsilon,iter_num):
    theta = initial_theta
    #记录中间过程最佳
    history_theta.append(theta)
    iter = 0
    while iter < iter_num:
        iter = iter + 1
        last_theta = theta 
        gradient = dLF(theta)
        theta = theta - eta * gradient
        history_theta.append(theta)
        if (abs(lossFunction(theta)-lossFunction(last_theta))<epsilon):
            break
    print('执行次数：%f',len(history_theta))
def plot_theta_history():
        plt.plot(plot_x,plot_y)
        plt.plot(np.array(history_theta),lossFunction(np.array(history_theta)),color='red',marker='+')
        plt.show()

eta = 0.01
history_theta = []
gradient_descent(1.0 , eta ,1e-6)
history_theta = []
eta = 0.01
gradient_descent2(1.0 , eta ,1e-6,100)

分别执行了285次 和101次

同时对于损失函数加个异常处理try: except:处理

线性回归中的梯度实现：
$$复习一下多元线性回归：{\hat{y}}_{(i)}={\theta }_0{X}_0^{(i)}+{\theta }_1{X}_1^{(i)}+{\theta }_2{X}_2^{(i)}+...+{\theta }_n{X}_n^{(i)}$$

$$已知训练数据样本x,y，找到{\theta }_0,{\theta }_1,...,{\theta }_n，使得损失函数尽可能小J=\sum _{i=1}^m(y^{(i)}-{\hat{y}}^{(i)}{)}^2$$

之前是使用推导公式直接求解，现在利用梯度下降法来求损失函数最小值。

将预测值带入损失函数J中，得到
$$J=\sum _{i=1}^m(y^{(i)}-{\theta }_0{X}_0^{(i)}-{\theta }_1{X}_1^{(i)}-{\theta }_2{X}_2^{(i)}-...-{\theta }_n{X}_n^{(i)}{)}^2最小$$

$${\hat{y}}_{(i)}=X_b\theta ,变换后有J=\sum _{i=1}^m(y^{(i)}-X_b\theta {)}^2$$

$$X_b=\left(\begin{array}{ccccc}1& X_1^{(1)}& X_2^{(1)}& ...& X_n^{(1)}\\ 1& X_1^{(2)}& X_2^{(2)}& ...& X_n^{(2)}\\ ...& & & & ...\\ 1& X_1^{(m)}& X_2^{(m)}& ...& X_n^{(m)}\end{array}\right)$$

$$\theta =\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ ...\\ {\theta }_n\end{array}\right)$$

$$\theta 这里是列向量n+1行，X_b为m\ast (n+1)矩阵$$

$$▽J(\theta )=\left(\begin{array}{c}\frac{\partial J}{\partial {\theta }_0}\\ \frac{\partial J}{\partial {\theta }_1}\\ \frac{\partial J}{\partial {\theta }_2}\\ ...\\ \frac{\partial J}{\partial {\theta }_n}\end{array}\right)$$

$$▽J(\theta )=\left(\begin{array}{c}{\sum }_{i=1}^{m}2(y^{(i)}-X_b^{(i)}\theta )(-X_0^{(i)})\\ {\sum }_{i=1}^{m}2(y^{(i)}-X_b^{(i)}\theta )(-X_1^{(i)})\\ {\sum }_{i=1}^{m}2(y^{(i)}-X_b^{(i)}\theta )(-X_2^{(i)})\\ ...\\ {\sum }_{i=1}^{m}2(y^{(i)}-X_b^{(i)}\theta )(-X_n^{(i)})\end{array}\right)$$

$$▽J(\theta )=\left(\begin{array}{c}{\sum }_{i=1}^{m}2(X_b^{(i)}\theta -y^{(i)})(X_0^{(i)})\\ {\sum }_{i=1}^{m}2(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})\\ {\sum }_{i=1}^{m}2(X_b^{(i)}\theta -y^{(i)})(X_2^{(i)})\\ ...\\ {\sum }_{i=1}^{m}2(X_b^{(i)}\theta -y^{(i)})(X_n^{(i)})\end{array}\right)$$

消除m的影响，每一列都除以m,
$$▽J(\theta )=\frac{2}{m}\left(\begin{array}{c}{\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_0^{(i)})\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_2^{(i)})\\ ...\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_n^{(i)})\end{array}\right)$$
这里看每一行则看做
$$C=\frac{2}{m}(X_b^{(1)}-y^{(1)},X_b^{(2)}-y^{(2)},X_b^{(3)}-y^{(3)},...,X_b^{(m)}-y^{(m)})$$
矩阵C点乘矩阵Xb ===》C矩阵列向量 (m*1) ； Xb矩阵 （m*(n+1)）

这里可以直接根据Xb的转置
$$▽J(\theta )=\frac{2}{m}{X}_b^{T}C=\frac{2}{m}{X}_b^T(X_b\theta -y)$$

$$X_b^T=\left(\begin{array}{ccccc}{X}_0^{(1)}& X_0^{(2)}& X_0^{(3)}& ...& X_0^{(m)}\\ X_1^{(1)}& X_1^{(2)}& X_1^{(3)}& ...& X_1^{(m)}\\ X_2^{(1)}& X_2^{(2)}& X_2^{(3)}& ...& X_2^{(m)}\\ ...& & & & ...\\ X_n^{(1)}& X_n^{(2)}& X_n^{(3)}& ...& X_n^{(m)}\end{array}\right)$$

import numpy as  np
import pandas  as  pd

'''计算均方误差MSE'''
def mse(train_y ,hat_y):
    assert len(train_y)==len(hat_y),"估计值和实际值的行数应该相同"
    assert len(train_y) > 0,"训练标签维度必须要大于0"
    return  sum((train_y-hat_y)**2)/len(train_y)
def fit_gd(train_X,train_Y, initial_theta ,eta ,epsilon, iter_max):
    '''根据训练数据集train_X ,train_Y,使用梯度下降算法训练线性回归'''
    assert train_X.shape[0] == len(train_Y) ,"确保训练数据行数相同"
    assert train_X.shape[1] + 1 == len(initial_theta),"确保训练参数行数与训练数据X的列数+1相等"
    
    theta = initial_theta
    '''给train_X添加一列全1 ''' 
    m = train_X.shape[0]
    train_x = np.hstack([np.ones(m).reshape(-1,1),train_X])
    iter = 1
    while iter < iter_max:
        iter = iter + 1
        last_theta = theta
        theta = theta -train_x.T.dot(train_x.dot(theta)-train_Y) * 2. / m * eta
        '''计算此时损失函数与y值之间的均方误差'''
        hat_y = np.dot(train_x,theta)
        mse_epsilon = mse(train_Y,hat_y)
        if  mse_epsilon < epsilon :
            break
    return  theta
def score(y_true,y_pred):
    u = sum((y_true - y_pred)**2)
    v = sum((y_true - y_true.mean())**2)
    return 1- u / v

from  sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn import model_selection

import random
boston = datasets.load_boston()
standardScaler = StandardScaler()
X = boston.data
y = boston.target
#归一化
standardScaler.fit(X)
x_std = standardScaler.transform(X)
#划分训练集与测试集
x_train ,x_test ,y_train,y_test = model_selection.train_test_split(x_std,y,test_size=0.2,random_state=666)
x_train = x_train[y_train < 50.0]
y_train = y_train[y_train < 50.0]
#x_train,y_train
initial_theta = np.zeros(x_train.shape[1] + 1)
eta = 0.005
epsilon = 1e-6
iter_max = 2000
theta = fit_gd(x_train,y_train,initial_theta , eta, epsilon , iter_max)
#出现了异常，溢出
#测试数据集
x_test = np.hstack((np.ones(x_test.shape[0]).reshape(-1,1),x_test))
y_pred = np.dot(x_test, theta)
score(y_test, y_pred) ,theta

这里与使用公式直接推导出的参数的结果类似：

(0.6443606885710433,
 array([ 2.20380509e+01, -6.56825615e-01,  8.36945559e-01, -3.51688663e-03,
         8.76600256e-02, -1.35159650e+00,  2.85093629e+00, -5.07274323e-01,
        -2.38713802e+00,  2.08562386e+00, -2.58930990e+00, -1.82071566e+00,
         9.02164783e-01, -2.45314781e+00]))

#公式推导方法 二：
import numpy  as  np
class  LinearRegression:
    def __init__(self):
        '''初始化线性回归模型'''
        self.coef_ = None  #
        self.interception = None  #截距
        self._theta = None
    
    def fit(self,x_train,y_train):
        assert x_train.shape[0]==y_train.shape[0],"特征值与标签要个数对应"
        x0 = np.ones(x_train.shape[0]).reshape(-1,1)
        x_b = np.hstack((x0,x_train))
        #self._theta = (((x_b.T).dot(x_b)).I).dot(x_b.T).dot(y_train)
        try:
            self._theta = np.linalg.inv(x_b.T.dot(x_b)).dot(x_b.T).dot(y_train)
            self.coef_ = self._theta[1:]
            self.interception = self._theta[0]
        except:
            print("矩阵不可逆")
       
        return  self
    
    def predict(self,x_predict):
        assert self._theta is not None ,"清先训练再预测"
        assert len(self.coef_) == x_predict.shape[1],"特征数要一致" 
        x_b = np.hstack((np.ones(x_predict.shape[0]).reshape(-1,1),x_predict))
        return x_b.dot(self._theta)
    def normalization(train_x):
        return None
    def score(self,y_true,y_pred):
        u = sum((y_true - y_pred)**2)
        v = sum((y_true - y_true.mean())**2)
        return 1- u / v
        
    def __repr__(self):
        return "LinearRegression()"

from  sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn import model_selection

import random
boston = datasets.load_boston()
linearRegression = LinearRegression()

standardScaler = StandardScaler()
X = boston.data
y = boston.target
#归一化
standardScaler.fit(X)
x_std = standardScaler.transform(X)
#划分训练集与测试集
x_train ,x_test ,y_train,y_test = model_selection.train_test_split(x_std,y,test_size=0.2,random_state=666)
x_train = x_train[y_train < 50.0]
y_train = y_train[y_train < 50.0]
#x_train,y_train
linearRegression.fit(x_train,y_train)
y_pred = linearRegression.predict(x_test)
#y_pred.shape ,y_test.shape
linearRegression.score(y_test, y_pred) 

0.64436044965542

上述算法是，批量梯度下降法 BGD batch gradient descent 每次计算时，是将整个数据都计算在内，在遇到数据量很大时需要计算很久。

针对该种缺点，随机梯度下降法 SGD stochastic gradient descent .随机梯度下降是每次迭代使用一个样本对参数进行更新。虽然不是每次迭代得到的损失函数都是朝向全局最优方向，但是大方向是向全局最优解，最终结果也应该在全局附近。其速度更快。

随机取出i个向量计算参与计算，同时eta 参数的变化，
$$eta=\frac{n_0}{(1+iter)}，iter=0,1,2,...,ite{r}_{max}$$

import matplotlib.pyplot as plt
iter_max  = 1000
x0 = np.linspace(0,1000,iter_max)
y = 1. / (1+ x0)
plt.plot(x0,y)
plt.show()

可以看出 eta的变化，刚开始变化迅速，后面就会缓慢，防止寻找最优解过程的跳动。

调试参数的过程这里就不涉及了，总结一下：

本文先从下山法开始进行原理解析，随后从简单代码实现；之后讲解最优化搜索的梯度下降法，及其代码实现；接着从梯度下降整体数据参与的优缺点，引出了随机梯度下降优化，在居士和饼干的文章中还有对于小批量梯度下降MBGD的描述，这里不再赘述。

四、打卡链接

【腾讯文档】机器学习小组第五周（20200316-20200322)

https://blog.csdn.net/wait_for_eva/article/details/84307306 各种typora公式符号