回归算法学习笔记——线性回归、随机梯度（SGD、BGD）、逻辑回归（牛顿法）、Softmax回归算法、L1/L2正则化、Ridge、Lasso、ElasticNet

线性回归

• 回归分析
• 目标函数：线性回归方程 $y=wx+b$
• 一个或多个自变量和因变量之间的关系进行建模(其中${\theta }_i$为权重，${\theta }_0$为bias偏置值)：
  • 一维特征：$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1$
  • 二维特征：$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2$
  • N维特征：${\sum }_{i=0}^N{\theta }_i{x}_i={\theta }^{T}X$

梯度下降算法

构建损失函数

• 未得到目标线性方程，得到上述公式中的${\theta }^T$参数
• 为确定所选的${\theta }_T$效果好坏使用损失函数（loss function）来评估$h(x)$函数的好坏
• 损失函数如下(MSE误差平方和）：
  $$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h({\theta }^T{x}_i)-y_i{)}^2\text{\hspace{1em}(1)}$$

梯度下降法

• 使用了导数的概念，对点$x_0$的导数反映了函数在该点处的瞬时变化速率
  $$f^{\prime }(x_0)=\underset{△x\to \infty }{\lim }=\frac{△y}{△x}=\underset{△x\to \infty }{\lim }=\frac{f(x_0+△x)-f(x_0)}{△x}\text{\hspace{1em}(2)}$$
• 推广到多维函数中，梯度反映了多维图形中变化速率最快的方向
• 起始时导数为正，$\theta$减小后并以新的$\theta$为基点重新求导，一直迭代就会找到最小的$J(\theta )$
  
• StochasticGradientDescent(SGD)随机梯度下降
  Repeat { ⋯ θ j = θ j − α ( ( h ( θ T x i ) − y i ) x i ) ⋯ } \text{Repeat} \{ \\ \dotsb\\ \theta_j = \theta_j -\alpha((h(\theta^Tx_i)-y_i)x_i)\\ \dotsb\\ \} Repeat{ ⋯θj​=θj​−α((h(θTxi​)−yi​)xi​)⋯}
  其中：$x_i$为随机的样本点
• 代码实现

import random
import matplotlib.pyplot as plt

x = [(0.,3), (1.,3) ,(2.,3), (3.,2), (4.,4), (0.,3) , (1.,3.1) ,(2.,3.5), (3.,2.1) , (4.,4.2)]
#y[i] is the output of y = theta0 * x[0] + theta1 * x[1] +theta2 * x[2]
y = [95.364,97.217205,75.195834,60.105519,49.342380, 100.364,100.217205,100.195834,100.105519,12.342380]

epsilon = 0.001
#learning rate
alpha = 0.001
diff = [0,0]
error1 = 0
error0 =0
m = len(x)

#init the parameters to zero
theta0 = 0
theta1 = 0
theta2 = 0
epoch = 0

error_array = []
epoch_array = []
while True:
    #calculate the parameters
    # 线性回归：h(x) = theta0  + theta1 * x[i][0] + theta2 * x[i][1] 
    # 损失函数(评测指标MSE)：累和 (1/2) *  (y - h(x)) ^ 2
    # d((1/2) *  (y - h(x)) ^ 2) / dtheta0 = (y - h(x)) * (-1) * (1) 
    # theta0 = theta0 - (-alpha * (y - h(x))* 1 )
    # theta1 = theta1 - (-alpha * (y - h(x))* x[i][0])
    # theta2 = theta2 - (-alpha * (y - h(x))* x[i][1])
    # 1. 随机梯度下降算法在迭代的时候，每迭代一个新的样本，就会更新一次所有的theta参数。
    i = random.randint(0, m - 1)

    # (y - h(x))
    diff[0] = y[i]-( theta0 * 1 + theta1 * x[i][0] + theta2 * x[i][1] )
    # - (y - h(x))x
    gradient0 = - diff[0]* 1
    gradient1 = - diff[0]* x[i][0]
    gradient2 = - diff[0]* x[i][1]
    # theta = theta - (  - alpha * (y - h(x))x )
    theta0 = theta0 - alpha * gradient0
    theta1 = theta1 - alpha * gradient1
    theta2 = theta2 - alpha * gradient2
    #theta3
    #calculate the cost function
    error1 = 0
    # 此处error为一个相对的Error值。
    for i in range(len(x)):
        error1 += (y[i]-( theta0 * 1 + theta1 * x[i][0] + theta2 * x[i][1]))**2 
    
    error1 = error1 / m    
    print("delta  error {}".format(abs(error1-error0)))
    error_array.append(error1)
    epoch_array.append(epoch)
    if epoch == 500:
        break
    else:
        error0 = error1
    epoch += 1
    print(' theta0 : %f, theta1 : %f, theta2 : %f, sgd error1 : %f, epoch : % f'%(theta0,theta1,theta2,error1, epoch))

print('Done: theta0 : %f, theta1 : %f, theta2 : %f'%(theta0,theta1,theta2))

plt.plot(epoch_array, error_array, color='blue', linewidth=3)
plt.xlabel('epochs')
plt.ylabel('errors')
plt.show()

• BatchGradientDescent(BGD)批量随机梯度下降
  Repeat { ⋯ θ j = θ j − α ( 1 m ∑ i = 0 m ( h ( θ T x i ) − y i ) x i ) ⋯ } \text{Repeat} \{ \\ \dotsb\\ \theta_j = \theta_j -\alpha(\frac{1}{m}\sum_{i=0}^m(h(\theta^Tx_i)-y_i)x_i)\\ \dotsb\\ \} Repeat{ ⋯θj​=θj​−α(m1​i=0∑m​(h(θTxi​)−yi​)xi​)⋯}
• 代码实现

import matplotlib.pyplot as plt

#Training data set
#each element in x represents (x0,x1,x2)
x = [(0.,3) , (1.,3) ,(2.,3), (3.,2) , (4.,4), (0.,3) , (1.,3.1) ,(2.,3.5), (3.,2.1) , (4.,4.2)]
#y[i] is the output of y = theta0 * x[0] + theta1 * x[1] +theta2 * x[2]
y = [95.364,97.217205,75.195834,60.105519,49.342380, 100.364,100.217205,100.195834,100.105519,12.342380]


#learning rate
alpha = 0.001 
diff = [0,0]
error1 = 0
error0 =0
m = len(x)

#init the parameters to zero
theta0 = 0
theta1 = 0
theta2 = 0
sum0 = 0
sum1 = 0
sum2 = 0

epoch = 0
error_array = []
epoch_array = []
while True:
    #calculate the parameters
    # 线性回归：hi(x) = theta0 + theta1 * x[i][1] + theta2 * x[i][2]  
    # 损失函数：(1/2) 累加 * (y - h(x)) ^ 2
    # theta = theta - 累和(  - alpha * (y - h(x))x )
    # 1. 随机梯度下降算法在迭代的时候，每迭代一个新的样本，就会更新一次所有的theta参数。
    #calculate the parameters
    # 2. 批梯度下降算法在迭代的时候，是完成所有样本的迭代后才会去更新一次theta参数
    print(m)
    for i in range(m):
        #begin batch gradient descent
        diff[0] = y[i]-( theta0 + theta1 * x[i][0] + theta2 * x[i][1] )
        sum0 = sum0 + ( -alpha * diff[0]* 1)
        sum1 = sum1 + ( -alpha * diff[0]* x[i][0])
        sum2 = sum2 + ( -alpha * diff[0]* x[i][1])
        #end  batch gradient descent
    
    theta0 = theta0 - sum0 / m
    theta1 = theta1 - sum1 / m
    theta2 = theta2 - sum2 / m

    sum0 = 0
    sum1 = 0
    sum2 = 0
    #calculate the cost function
    error1 = 0
    for i in range(m):
        error1 += ( y[i]-( theta0 + theta1 * x[i][0] + theta2 * x[i][1] ) )**2
        
    error1 = error1 / m
    error_array.append(error1)
    epoch_array.append(epoch)
    if epoch == 200:
        break
    else:
        error0 = error1
    epoch += 1
    print(' theta0 : %f, theta1 : %f, theta2 : %f, bgd error1 : %f, epoch: %f'%(theta0,theta1,theta2,error1,epoch))

print('Done: theta0 : %f, theta1 : %f, theta2 : %f'%(theta0,theta1,theta2))
plt.plot(epoch_array, error_array, color='blue', linewidth=3)
plt.xlabel('epochs')
plt.ylabel('errors')
plt.show()

对比两类梯度下降法各有利弊，SGD的计算量小，但随机性强，有时可以快速收敛，有时则收敛的很慢，loss曲线抖动下降；MGD的计算量大，但包含批量的样本，收敛效果比较好，loss曲线平滑下降。

Logistic Regression算法

• 针对分类问题，利用回归问题的思路，解决分类问题。

sigmoid函数

• sigmoid函数是一个s形的曲线
• 取值在(0, 1)之间
• 在远离0的地方函数的值会很快接近0或1
• $g(z)=\frac{1}{1+e^{-z}}$

构造目标函数

• ${\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n={\sum }_{i=0}^n{\theta }_i{x}_i={\theta }^{T}X$
• $h_{\theta }(x)=g({\theta }^{T}X)=\frac{1}{1+e^{-{\theta }^{T}X}}$
• $$\begin{gathered} p(y=1|x,\theta )=h_{\theta }(x) \\ p(y=0|x,\theta )=1-h_{\theta }(x)\text{\hspace{1em}(3)} \end{gathered}$$
• 将预测结果转化成预测某类别的概率公式(3)中y为真实的类别

构造损失函数-极大似然估计

• 根据上述公式构建损失函数，目的：通过${\theta }^T$的取值使预测正确类别的概率最大化
• $p(y|x,\theta )=(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$
• Lossfunction: $L(\theta )={\prod }_{i=0}^{m}p(y_i|x_i,\theta )={\prod }_{i=0}^m(h_{\theta }(x_i){)}^{y_i}(1-h_{\theta }(x_i){)}^{1-y_i}$
• 由于累乘的式子难以求导，因此对上式进行对数处理，转化为累加
• $l(\theta )=\log L(\theta )={\sum }_{i=0}^m(y_i(\log h_{\theta }(x_i)+(1-y_i)\log (1-h_{\theta })))$
• 转化为类似线性回归问题中求最小损失函数的问题：
  • $J(\theta )=-\frac{1}{m}l(\theta )$

梯度下降

• 首先求解损失函数梯度
  $$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}J(\theta )& =-\frac{1}{m}\sum _{i=1}^m\left[y_i\frac{1}{h_{\theta }(x_i)}\frac{\partial }{\partial {\theta }_j}{h}_{\theta }(x_i)-(1-y_i)\frac{1}{1-h_{\theta }(x_i)}\frac{\partial }{\partial {\theta }_j}{h}_{\theta }(x_i)\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left(h_{\theta }(x_i)-y_i\right)x_i^j\end{array}$$
• 使用梯度下降求解最优参数
  $${\theta }_j={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m\left(h_{\theta }(x_i)-y_i\right)x_i^j$$

多分类问题

• 解决方案：one vs rest
  • 适合场景：类别可以同时出现，互相不互斥的情况
  • 参考资料：one-vs-rest与one-vs-one以及sklearn的实现

优化算法：牛顿法

• 牛顿法比梯度下降法收敛的要快（迭代次数更少）
  • 梯度下降法每次只从当前位置选一个上升速度最大的方向走一步
  • 牛顿法在选择方向时，不仅会考虑上升速度是否够大，还会考虑你走了一步之后，上升速度是否会变得更大（即二阶导数，类似物理上的加速度）
• 红色的牛顿法的迭代路径，绿色的是梯度下降法的迭代路径。
  
• 由于最速梯度下降收敛速度并不“最速”，局部搜索中最速下降的方向和全局的最小值方向并不一致，所以就有后来改进的方法，包括牛顿法以及拟牛顿法。

切线法

• 牛顿法的几何意义
  
• 从上图可以总结出如下公式：
  $${\theta }_j={\theta }_j-\frac{J^{\prime }({\theta }_j)}{J^{\prime \prime }({\theta }_j)}$$
  其中分母上的二阶导函数可以理解为正则项，当一阶导数下降方向不是最优方向，进行修改，而且由于二阶导数包含数值，因此也不用设置学习率，下降的步长是最优的。

另一种理解方式

• 随时函数最小值 = 损失函数导数为0
• 将损失函数进行泰勒展开：
  $$f(x)\approx f\left(x^{(k)}\right)+g_k^T\left(x-x^{(k)}\right)+\frac{1}{2}{\left(x-x^{(k)}\right)}^{T}H\left(x^{(k)}\right)\left(x-x^{(k)}\right)\text{\hspace{1em}(4)}$$
  其中$g(x)$为一阶导数，$H$为二阶导数
• 对(4)式等式两端进行求导，转化为梯度下降：
  $$\begin{gathered} ▽f(x)=g_k+H_k\left(x-x^{(k)}\right) \\ {g}_k+H_k\left(x-x^{(k)}\right)=0 \\ {x}^{(x+1)}=x^{(x)}-H_k^{-1}{g}_k \end{gathered}$$

  最终的结果与切线法得到的结果表达式是一样的

改进：拟牛顿法

• 牛顿法中需要计算$H_k$二阶导数矩阵，其尺寸为$N\times N$的，其中N为特征的数量，计算代价很大
• 而拟牛顿法不在使用$H_k$二阶导数矩阵，而是选取相似矩阵

Softmax Regression算法

• Softmax解决多分类的分类器，即类标签y的取值大于等于2
• 类标记为: y ( i ) ∈ { 1 , 2 , ⋯   , k } y(i) \in \left\{ 1,2,\cdots,k \right\} y(i)∈{ 1,2,⋯,k}
• 假设函数为对于每一个样本估计其所属的类别的概率$p(y=j\mid x)$
• 每一个样本估计所属类别概率为：
  $$p\left(y^{(i)}|x^{(x)};\theta \right)=\frac{e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}}$$
  即对预测数的各类别的值取$log$后除以所有类别概率值取$log$的总和

Softmax回归代价函数

• 类似于Logistic回归，在Softmax的代价函数中引入指示函数 I { ⋅ } I\{⋅\} I{ ⋅}，其具体形式为：
  I { expression } = { 0 if expression = false 1 if expression = true I\left \{ \text{expression}\right\} = $$\begin{cases} 0 \quad \text{if}\ \text{expression}=\text{false}\\ 1 \quad \text{if}\ \text{expression}=\text{true} \end{cases}$$ I{ expression}={ 0if expression=false1if expression=true​
• 那么，对于Softmax回归的代价函数为交叉熵：
  J ( θ ) = − 1 m [ ∑ i = 1 m ∑ j = 1 k I { y i = j } log ⁡ e θ j T x ( i ) ∑ l = 1 k e θ l T x ( i ) ] J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^kI\{y^{i}=j \}\log \frac{e^{\theta^T_jx^{(i)}}}{\sum_{l=1}^ke^{\theta^T_lx^{(i)}}} \right] J(θ)=−m1​[i=1∑m​j=1∑k​I{ yi=j}log∑l=1k​eθlT​x(i)eθjT​x(i)​]
• Sofrmax回归求解
  • 对上述的代价函数，可以使用梯度下降法进行求解，首先对其进行求梯度：
    $$\begin{array}{rl}{▽}_{{\theta }_j}J(\theta )& =-\frac{1}{m}\sum _{i=1}^m\left[{▽}_{{\theta }_j}\sum _{j=1}^{k}I\left\{y^{(i)}=j\right\}\log \frac{e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[I\left\{y^{(i)}=j\right\}\frac{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}}{e^{{\theta }_j^T{x}^{(i)}}}\cdot \frac{e^{{\theta }_j^T{x}^{(i)}}\cdot x^{(i)}\cdot {\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}-e^{{\theta }_j^T{x}^{(i)}}\cdot x^{(i)}\cdot e^{{\theta }_j^T{x}^{(i)}}}{{\left({\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}\right)}^2}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[I\left\{y^{(i)}=j\right\}\cdot \frac{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}-e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}}\cdot x^{(i)}\right]\end{array}$$
    $${\theta }_j={\theta }_j-{▽}_{{\theta }_j}J(\theta )$$

L1/L2正则化

• 为了防止过拟合，拟合过程中通常都倾向于让权值尽可能小，最后构造一个所有参数都比较小的模型
• 参数值小的模型比较简单，能适应不同的数据集，也在一定程度上避免了过拟合现象
• 对于一个线性回归方程，若参数很大，那么只要数据偏移一点点，就会对结果
  造成很大的影响
• 但如果参数足够小，数据偏移得多一点也不会对结果造成什么影响

L1

• L1正则化是指权值向量$w$中各个元素的绝对值之和，通常表示为${\Vert w\Vert }_1$
• 如下图所示，图中等值线是J0的等值线，黑色方形是L函数的图形。在图中，当J0等值线与L首次相交的地方就是最优解。上图中J0与L在L的一个顶点处相交，这个顶点就是最优解。
  
• 正则化和“带约束的目标函数”是等价的对L1正则化建立数学模型如下:
  $$\begin{gathered} J(\theta )=\sum _{i=1}^n{\left(y_i-{\theta }^T{x}_i\right)}^2 \\ s.t.{\Vert \theta \Vert }_1⩽m \end{gathered}$$
• 通过拉格朗日乘子法，可以变为如下形式：
  $$J(\theta )=\sum _{i=1}^n{\left(y_i-{\theta }^T{x}_i\right)}^2+\lambda \left({\Vert \theta \Vert }_1-m\right)$$

L2

• L2正则化是指权值向量$w$中各个元素的平方和然后再求平方根（在Ridge回归的L2正则化项有平方符号），通常表示为${\Vert w\Vert }_2$
• 如下图所示，二维平面下L2正则化的函数图形是个圆，与方形相比，被磨去了棱角。
  
• 对比L1和L2两图，采用L1范数时平方误差项等值线与正则化项等值线的交点出现在坐标轴上，即$w_1$或$w_2$为0，而采用L2范数时，两者的交点常出现在某个象限中，即$w_1$或$w_2$均非0；换而言之，采用L1范数比L2范数更容易于得到稀疏解。
• 正则化和“带约束的目标函数”是等价的对L1正则化建立数学模型如下:
  $$\begin{gathered} J(\theta )=\sum _{i=1}^n{\left(y_i-{\theta }^T{x}_i\right)}^2 \\ s.t.{\Vert \theta \Vert }_2^2⩽m \end{gathered}$$
• 通过拉格朗日乘子法，可以变为如下形式：
  $$J(\theta )=\sum _{i=1}^n{\left(y_i-{\theta }^T{x}_i\right)}^2+\lambda \left({\Vert \theta \Vert }_2^2-m\right)$$

L1和L2对比

其中c为超参数，控制L1和L2的限制区域大小，c越小对损失函数的约束越大；图中像素点表示[0, 1]的取值，0为白色1为黑色。可以看出L1下是很多权值为0，做到了稀疏性。

参考资料：L1与L2范数

正则化目的

• 预防或修复过拟合情况
• 其他方式：
  • 数据层面：
    • 增加数据（行上扩展）
    • 减少feature，特征选择
  • 模型层面：
    • 模型融合

Ridge与Lasso

• L1正则化和L2正则化可以看做是损失函数的惩罚项
• Lasso回归
  • 对于线性回归模型，使用L1正则化的模型建叫做Lasso回归
  • Lasso回归的损失函数，式中加号后面一项即为L1正则化项
• Ridge
  • Ridge回归的损失函数，式中加号后面一项$\alpha {\Vert \theta \Vert }_2^2$即为L2正则化项

ElasticNet

对于ElasticNet回归，可以看做是Lasso和ridge回归的组合
$$\begin{gathered} J(\theta )=\sum _{i=1}^n{\left(y_i-{\theta }^T{x}_i\right)}^2+{\lambda }_1({\Vert \theta \Vert }_1)+{\lambda }_2({\Vert \theta \Vert }_2^2) \\ \underset{\theta }{\min }J(\theta ) \end{gathered}$$

转载：https://blog.csdn.net/weixin_43870329/article/details/106604834