ML-梯度下降法求解线性回归

作者: 倪桦 | 来源:发表于2023-01-26 01:34 被阅读0次

梯度法 是基于搜索来最优化一个目标函数的方法。
分为梯度下降法 和 梯度上升法 :

梯度下降法 用来最小化一个损失函数；

梯度上升法，用作最大化效用函数。

对于很多无法求取数学解(类似线性回归的正规数学方程解)的机器学习算法模型，就需要用到 梯度法 这种搜索最优化方法来找到其最优解。

1、梯度含义

对于一个存在极值点的损失函数 $J$ ，比如其与参数 $\theta$ 的关系在平面中描述为二次曲线:

对于损失函数曲线上的一个点的导数 $\frac {\mathrm{d} J}{\mathrm{d} \theta}$ （切线斜率）其本质刻画了使该点增大的移动方向。如上图点A处的导数为负数，意味这使A点函数值增大对应的是x轴负方向上的 $\theta$ ，或者说通过减小 $\theta$ 来增大函数值；导数值的绝对值描述了该点切线的倾斜程度--梯度。函数极值点处的切线斜率为零(梯度为零)，当曲线上任意点滚落到曲线极值点的时候，这段轨迹上的梯度是逐渐变小的过程(落在极值点梯度变为0)-- 梯度下降。从而通过使曲线上一点滚落到函数的极值点来取得损失函数最优解的方法称为梯度下降法。
$\therefore$ 要使得曲线上任意点往函数值变小的方向移动一段距离，结合导数代表的移动方向和移动步长 $\eta$ 表示为:
$-\eta \frac {\mathrm{d} J}{\mathrm{d} \theta}$
在多维函数中，通过对各个方向的分量分别求导，最终得到的方向就是梯度。

关于步长 $\eta$ 的取值问题：
步长 $\eta$ 在机器学习领域又称 学习率(learning rate)，它的取值会影响获得最优解的速度；当 $\eta$ 取值不合适的时候，甚至会得不到最优解； $\eta$ 也是梯度下降法的一个超参数。

当 $\eta$ 取值太小，意味着损失函数曲线上的一点向函数最低点的滚动速度很慢，影响算法得到最优解的速度，也就是减慢收敛学习速度；
当 $\eta$ 取值太大，意味着损失函数曲线上的一点每次移动的距离很长，这种长距离移动甚至会导致该点无法滚落到函数的最低点，出现不收敛问题；
对于大多数ML算法， $\eta = 0.01$ 是一个相对合适的取值

2、梯度下降搜索注意事项

对于类似二次函数这种明显有唯一极值点的损失函数，它的最值点可以用数学解直接求得。梯度下降法 更多地是应用在 复杂损失函数 的极值点求解。

2.1 多极值点目标函数最优解

多极值点函数
面对存在多极值点的函数，使用 梯度下降法，如果滚落起始点选的不合适，一次搜索很可能只搜索到一个局部最优解 ，而不得到全局最优解。
解决方案：

多次运行，随机化初始点，并比较模型性能

梯度下降法的初始点也是一个超参数

该方案无法保证一定能取得全局最优解 ，但能取得相对较好的 局部最优解 。

3、梯度下降python实现

3.1 算法基本过程

import numpy as np
### Prepare dataset, 定义二次关系函数J，生成x和y
def f_J(x):
    try:
        return (x - 2.5)**2 -1
    except:
        return float('inf') ### 异常检测，如果出现异常返回浮点数最大值float('inf')

x = np.linspace(-1,6,141)
y = [f_J(x_i) for x_i in x]

### Derivative function
def Derivate(x):
    return 2*(x - 2.5)

### Gradient Descent
eta = .01 ### 学习率
epsilon = 1e-8 ### 定义一个误差范围，当初始点开始滚动，如果当前点位与一点位函数值误差小于定义的epsilon误差范围，认为当前点为就是最优解。
theta_lst = [float(np.random.randint(10,size = 1) + np.random.normal(size = 1))] ### 随机化一个滚落初始点
while True:
    gradient = Derivate(theta_lst[-1])
    theta_lst.append(theta_lst[-1] - eta * gradient)
    if(abs(f_J(theta_lst[-1]) - f_J(theta_lst[-2])) < epsilon):break
    if(len(theta_lst) > 1e4): break       ### 循环1e4次以后如果没有取得最优解将退出循环
        
print(theta_lst[-1]) ### 输出最优参数

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

从简单线性回归的梯度推广到多元线性回归 $\hat y^{(i)} = \theta_{0}x_{0}^{(i)} + \theta_{1}x_{1}^{(i)} ... + \theta_{n}x_{n}^{(i)},x_{0}^{(i)} = 1$ ，求解其目标函数 $J = \sum^{m}_{i} {(\hat y^{(i)} - y^{(i)})^{2}}$ 的最优解，通过对目标函数求每个 $\theta$ 的偏导获得该 $\theta$ 上的 梯度函数，综合成目标函数的梯度，描述为向量形式 :
$-\eta \nabla, \nabla J(\theta) = (\frac{\partial J}{\partial \theta_{0}},\frac{\partial J}{\partial \theta_{1}},...,\frac{\partial J}{\partial \theta_{n}})$

偏导结果：

目标函数 $\sum_{i=1}^{m}(y^{i} - \theta_{0}\cdot 1 - \theta_{1}x_{1}^{(i)} ... - \theta_{n}x_{n}^{(i)})^2$ , $m$ 为样本量。

令 $X_{b}^{(i)} = (1,x^{(i)}_1,...,x^{(i)}_n),\theta = (\theta_{0},\theta_{1},...,\theta_{n})^{T}$ ;

$\nabla J(\theta)= \begin{bmatrix} \partial J/ \partial \theta _{0} \\ \partial J/ \partial \theta _{1} \\ ... \\ \partial J/ \partial \theta _{n} \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{m}{2(y^{i} - X_{b}^{(i)}\theta)\cdot ( -1 \ \ )} \\ \sum_{i=1}^{m}{2(y^{i} - X_{b}^{(i)}\theta)\cdot (-x^{i}_1)} \\... \\ \sum_{i=1}^{m}{2(y^{i} - X_{b}^{(i)}\theta)\cdot (-x^{i}_n)} \end{bmatrix} = 2 \begin{bmatrix} \sum_{i=1}^{m}{(X_{b}^{(i)}\theta - y^{i})\cdot 1 \ \ } \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_1} \\... \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_n} \end{bmatrix}$

上式中，由于最后每个参数上 梯度值 的结果都受到了样本数量 $m$ 的影响 $\sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i})}$ ，为了脱离这种影响，通过除以样本量来实现:

$\nabla J(\theta)= \frac{2}{m} \cdot \begin{bmatrix} \sum_{i=1}^{m}{(X_{b}^{(i)}\theta - y^{i})\cdot 1 \ \ } \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_1} \\... \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_n} \end{bmatrix}$

这种改变，其对应的目标函数背后意味着使用MSE作为目标函数： $J(\theta) = MES(y,\hat y)= \frac{1}{m}\sum^{m}_{i} {(\hat y^{(i)} - y^{(i)})^{2}}$

令样本的特征矩阵为 $X$ ，
定义 $X_b = \begin{bmatrix} x_{0}^{(1)}&x_{1}^{(1)}&...&x_{n}^{(1)} \\ x_{0}^{(2)}&x_{2}^{(2)}&...&x_{n}^{(2)}\\ ...&...&...&... \\ x_{0}^{(m)}&x_{2}^{(m)}&...&x_{n}^{(m)}\end{bmatrix}$ ,其中 $x^{(i)}_{0} = 1$
则有：
$\nabla J(\theta) = \frac{2}{m} (X_b\theta - y)^{T} \cdot X_b \rightarrow$
使计算结果按列向量返回
$\nabla J(\theta) = \frac{2}{m} (X_b\theta - y)^{T} \cdot X_b = (\frac{2}{m} (X_b\theta - y)^{T} \cdot X_b ) ^{T} = \frac{2}{m}\cdot X_b^{T} \cdot (X_b\theta - y)$

算法过程的python 实现

### Prepare data
import numpy as np
from sklearn import datasets
boston = datasets.load_boston()
x = boston.data[ boston.target < 50 ]
y = boston.target[boston.target < 50]
from sklearn.model_selection import train_test_split
Train_X,Test_X,Train_Y,Test_Y = train_test_split(x,y,test_size= 0.2,random_state=666)
X_b = np.hstack([np.ones((len(Train_X),1)),Train_X])
y = Train_Y

### Loss Function
def J(theta,X_b,y):
    try:
        return np.sum((y - X_b.dot(theta))**2)/len(X_b)
    except:
        return float('inf') ### 异常检测，如果出现异常返回浮点数最大值float('inf')

### Derivative function
def dJ(X_b,y,theta):
    return (2./len(X_b)) * (X_b.T.dot(X_b.dot(theta) - y))

### Gradient Descent
eta = .000001 ### 学习率
epsilon = 1e-8 ### 定义一个误差范围，当初始点开始滚动，如果当前点位与一点位函数值误差小于定义的epsilon误差范围，认为当前点为就是最优解。
theta_lst = [np.random.randint(10,size = X_b.shape[1])+ np.random.normal(size = X_b.shape[1])] ### 随机化滚落初始点
while True:
    gradient = dJ(X_b,y,theta_lst[-1])
    theta_lst.append(theta_lst[-1] - eta * gradient)
    
    if(abs(J(theta_lst[-1],X_b,y) - J(theta_lst[-2],X_b,y)) < epsilon):break
    if(len(theta_lst) > 1e4): break       ### 循环1e4次以后如果没有取得最优解将退出循环
        
print(theta_lst[-1])

4、梯度下降法与数据归一化

在使用梯度下降法进行搜索的时候，由 $\eta$ 学习率和偏导结果决定了损失函数上初始点往极值滚落的运动方向和步长。这种方案下就会存在由于特征的数量级的不同，从而损失函数上各维度的滚动步长不一致，有些维度滚动幅度特别大，而有些特别小，出现收敛速度奇慢或者无法收敛的问题。

解决方案就是把样本特征规约到同一个尺度 数据归一化 。

from sklearn.preprocessing import StandardScaler
standScaler = StandardScaler()
standScaler.fit(Train_X)
Train_X_scale = standScaler.transform(Train_X)
X_b = np.hstack([np.ones((len(Train_X),1)),Train_X_scale])
y = Train_Y

5、随机梯度下降法

批量梯度下降法(Batch Gradient Descent)，上述通过纳入所有训练样本来计算目标函数上滚落初始点准确 梯度值 的方法。该方法存在缺点，在当样本量急剧膨胀的时候，就会出现效能下降的问题。
随机梯度下降法(Stochastic Gradient Descent)，该方法仅从训练样本中取一个样本来计算损失函数上滚落初始点的搜索方向，不综合所有样本来计算最优前进方向，参考等式的意思就是基于一个样本在当前 $\theta$ 参数下真值与预测值的变化来推测 $\theta$ 的变化方向。对于这样一个粗略的方向，每次的搜索运动虽然不能保证一定朝着损失函数上极小值点的方向前进，但大概率是向着这个方向前进的，随着搜索次数增加也能够来到损失函数上的极小值范围。这种方法牺牲一定的精度来节约时间成本，同时能够获得一个不错的解，是面对样本量急剧膨胀使用 批量梯度下降法 出现效能下降或运算困难情况时的可选解决方案。
在使用 随机梯度下降法 的时候，如果使用一个固定值的学习率，就可能会出现当滚落点来到了目标函数的极小值区域又随即跳出该区域的情况，所以应该随着搜索次数的增加相应地减小学习率使得滚落点稳定在目标函数的极小值区域 (模拟退火思想) 。这样学习率 $\eta$ 与搜索次数的关系函数可以表达为 $\eta = \frac {t_0}{i\_iters + t_1}$ ，式中常量 $t_0,t_1$ 用来防止学习率下降太快（如迭代次数从1到2，学习率就出现了50%的变化），如让 $\eta$ 从0.01开始变化，可取( $t_0 = 2, t_1 = 50$ )。

$\frac{2}{m} \cdot \begin{bmatrix} \sum_{i=1}^{m}{(X_{b}^{(i)}\theta - y^{i})\cdot 1 \ \ } \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_1} \\... \\ \sum_{i=1}^{m}{(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_n} \end{bmatrix} \rightarrow 2 \cdot \begin{bmatrix} {(X_{b}^{(i)}\theta - y^{i})\cdot 1 \ \ } \\ {(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_1} \\... \\ {(X_{b}^{(i)}\theta -y^{i} )\cdot x^{i}_n} \end{bmatrix}$
此时的 $X_b^{{i}}$ 仅是 $X_b$ 矩阵中的一行，也即一个样本。

5.2 随机梯度下降python实现

### Prepare data
import numpy as np
from sklearn import datasets
boston = datasets.load_boston()
x = boston.data[ boston.target < 50 ]
y = boston.target[boston.target < 50]
from sklearn.model_selection import train_test_split
Train_X,Test_X,Train_Y,Test_Y = train_test_split(x,y,test_size= 0.2,random_state=666)

### Scale data
from sklearn.preprocessing import StandardScaler
standScaler = StandardScaler()
standScaler.fit(Train_X)
Train_X_scale = standScaler.transform(Train_X)
X_b = np.hstack([np.ones((len(Train_X),1)),Train_X_scale])
y = Train_Y


### Derivative function
def sgd(X_b_i,y_i,theta):
    return 2 * (X_b_i.T.dot(X_b_i.dot(theta) - y_i))

### Stochastic Gradient Descent >>>
theta = np.random.randint(10,size = X_b.shape[1])+ np.random.normal(size = X_b.shape[1]) ### 随机化滚落初始点
for cur_iters in range(150): ### 使用随机梯度下降法遍历样本的时候，至少保证每个样本都要被纳入计算5次
    ### 打乱样本
    shuffle_index = np.random.permutation(len(X_b))
    X_b_shuffle = X_b[shuffle_index]
    y_shuffle = y[shuffle_index]
    
    for i in range(len(X_b_shuffle)): ### 遍历乱序样本
            gradient = sgd(X_b_shuffle[i],y_shuffle[i],theta)
            theta = theta - ( 5./(cur_iters*len(X_b) + i + 50) ) * gradient ### 更新每次迭代的 theta
    
### 随机梯度下降法由于搜索步长是随机的，与梯度值本身无关
### 即使运动两点的损失函数值误差很小，并不代表搜索运动来到了损失函数的极值区域，所以无法在损失函数上根据运动前后两点的差值来定义收敛域
### 退出随机梯度下降搜索是在搜索运动迭代一定次数后，接受搜索最后一次搜索的结果

print(theta)

5.3 scikit-learn 中的SGD

### Prepare data
import numpy as np
from sklearn import datasets
boston = datasets.load_boston()
from sklearn.model_selection import train_test_split
Train_X,Test_X,Train_Y,Test_Y = train_test_split(boston.data[ boston.target < 50 ],boston.target[boston.target < 50],test_size= 0.2,random_state=666)

### Scale data
from sklearn.preprocessing import StandardScaler
standScaler = StandardScaler()
standScaler.fit(Train_X)
Train_X_scale = standScaler.transform(Train_X)
Test_X_scale = standScaler.transform(Test_X)
### SGD regressor
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(n_iter_no_change = 100)
sgd_reg.fit(Train_X_scale,Train_Y)
sgd_reg.score(Test_X_scale,Test_Y) ### 模型准确度

除了样本极端的批量梯度下降法和随机梯度下降法，比较折中的方法是结合两者的小批量梯度下降法(Mini-Batch Gradient Descent)，使用全样本中的 $k$ 个样本来判断下降方向。

随机梯度下降法的优势还有：

一定概率可以跳出 局部最优解，找到整体的 全局最优解；

机器学习领域本身求解的是系列不确定性问题，随机性在ML中扮演重要角色。

6、任意形式目标函数的梯度

对于复杂目标函数，在数学上求解其梯度可能会变得非常苦难。这种情况下可以从 导数的定义 出发来近似出目标函数上一点的梯度。

切线斜率

切线的近似
导数的实质就是 曲线上切线的斜率 。根据微分思想，在曲线上任意点 $A$ 的两侧各区一点，当这两点相对 $A$ 点的偏移 $\epsilon$ 趋于无穷小，那么这两点连线的斜率就是点 $A$ 处切线的斜率，也就是该点的导数。所以，当取一个特别小的 $\epsilon$ ，这两点连线的斜率可以近似为点 $A$ 处的导数。
$\frac {\mathrm {d} J}{\mathrm {d} \theta} \approx \frac {J(\theta + \epsilon) - J(\theta - \epsilon)}{2\epsilon}$
推广到高维特征 $\theta = (\theta_0,\theta_1,..,\theta_n)$
$\nabla J(\theta) = (\frac{\partial J}{\partial \theta_{0}},\frac{\partial J}{\partial \theta_{1}},...,\frac{\partial J}{\partial \theta_{n}})$

设置 $\theta_{i}^{+} = (\theta_0,\theta_1,...,\theta_i+\epsilon,...,\theta_n)$
设置 $\theta_{i}^{-} = (\theta_0,\theta_1,...,\theta_i-\epsilon,...,\theta_n)$
目标函数 $J$ 对 $\theta_i$ 求导的结果即 $\frac {\partial J}{\partial \theta_i} = \frac {J(\theta_{i}^{+}) - J(\theta_{i}^{+})}{2\epsilon}$

6.2 近似梯度的 python实现

### Loss Function
def J(theta,X_b,y):
    try:
        return np.sum((y - X_b.dot(theta))**2)/len(X_b)
    except:
        return float('inf') ### 异常检测，如果出现异常返回浮点数最大值float('inf')

def approx_dJ(main_theta,X_b,y,epsilon = 1e-5):
    res = np.empty(len(main_theta))
    for i in range(len(main_theta)):
        upper_theta = main_theta.copy() ; upper_theta[i] += epsilon
        lower_theta = main_theta.copy() ; upper_theta[i] -= epsilon
        
        res[i] = (J(upper_theta,X_b,y) - J(lower_theta,X_b,y)) / (2 * epsilon)