LogisticRegression -逻辑回归/对率回归与三种优化算法（梯度下降/牛顿/拟牛顿与sklearn实现

文章目录

• 1. 简介
• 2. 实现思路
• 2.1 迭代计算系数
• a. 梯度下降法
• b. 牛顿迭代法
• c. 拟牛顿法（BFGS
• 2.2 Sigmoid 函数转化
• 3. 数据尝试
• 4. sklearn 实现
• 4.1 参数介绍
• 4.2 常用调用方法

1. 简介

逻辑回归即考察在样本各属性值前加上一个系数后的和（类似于加权平均），通过与阈值的比较实现分类

• 对于数据集D，样本由d个属性描述，此时我们试图学得：$f(x_i) = w^Tx_i + b \tag{1.1}$f(xi​)=wTxi​+b(1.1)

  • 也就是：$a_0+a_1x_1 + a_2x_2 + \cdots + a_nx_n = f(x_i)$a0​+a1​x1​+a2​x2​+⋯+an​xn​=f(xi​)

• 目标：求得最加拟合的一串系数+常数项

• 适用范围
  • 只适用于二分类问题
  • 所有的变量必须是数值型变量
  • 残差和因变量都要服从二项分布（是多次重复的伯努利试验
    • 注意二项分布对应的是分类变量，所以不是正态分布
    • 进而不是用最小二乘法，而是最大似然法来解决方程估计和检验问题
  • 自变量与概率之间要是线性关系，此时概率才能表示为自变量的线性组合
  • 各个观测对象之间互相独立

2. 实现思路

主要分为两个步骤

• 迭代计算系数
• 通过 Sigmoid 函数将$f(x_i)$f(xi​)进行转换并与阈值比较得出结论

2.1 迭代计算系数

• 令$\hat{\beta}=(w;b)$β^​=(w;b) , $\hat{x}= (x;1)$x^=(x;1)（注意此时的x和beta都是列向量

  • 即此时的目标计算函数转化为：$f(x_i) = \hat{\beta} \hat{x}$f(xi​)=β^​x^

  • 以下对参数 beta 求解

  
  

• 此时迭代的目标：即损失函数最小化（极大似然法

• 注意此时可以对似然函数两边取ln，将连乘转为连加，从而得到代价函数：
  $l(\beta) = \sum^m_{i=1}\left(-y_i\beta^T\hat{x_i}+ln\left(1+e^{\beta^T\hat{x_i}}\right)\right) \tag{2.1}$l(β)=i=1∑m​(−yi​βTxi​^​+ln(1+eβTxi​^​))(2.1)

• 令代价函数最小的过程实际上属于无约束优化问题

  • 以下考虑三种方法：
    • 梯度下降法
    • 牛顿法
    • 拟牛顿法

a. 梯度下降法

• 可以理解为每次按照梯度下降的方向下降一定的步长，以寻找最优解的位置
• 同理此时求的是代价函数的最小值 → 梯度下降；求最大值时用梯度上升法，即系数 1/m 为正

• 对于第 i 个样本的第 j 个属性（注意此时的j包括x^ 最后的1:$\theta_j^{new} = \theta_j - \alpha \frac{1}{m} \sum^m_{i=1}(h_\theta(x_i)-y_i)x_i^j\tag{2.2.1}$θjnew​=θj​−αm1​i=1∑m​(hθ​(xi​)−yi​)xij​(2.2.1)

  • 其中 $\theta_j$θj​ 表示θ内第j个值
  • $h_\theta(x)$hθ​(x) 表示y=1 时的后验概率：$h_\theta(x) = P(y = 1 | x;\theta) = \frac{e^{w^Tx+b}}{1+e^{w^Tx+b}} \tag{2.2.2}$hθ​(x)=P(y=1∣x;θ)=1+ewTx+bewTx+b​(2.2.2)
  • $x_i^j$xij​ 表示第i个样本第j个属性的取值
  • α 表示学习步长，取小一点比较精确同时也会比较慢

特点

• 梯度下降法比较万能，基本都能用
• 此处使用梯度下降，也就是说求似然函数的最小值（似然函数是代价函数），将1/m前改为 + 则梯度上升法，可以求最大值
• 梯度下降法一般很慢

python实现

def sag(old_beta,alpha):
    
    '''
    使用梯度下降法得到更新的系数值
    
    param old_beta: 初始系数值
    param alpha: 梯度下降学习步长
    return cur_beta: 迭代结束后计算得到的系数值
    '''
    cur_beta = old_beta
    cur_lb = 0
    
    while 1:
        for j in range(num_a):
            for i in range(number):
                wtxi = np.dot(datax[i,:],old_beta.reshape(-1,1))
                hx = np.exp(wtxi)/(1+np.exp(wtxi))    # 2.2.2
                cur_beta[0,j]-= alpha/number*(hx-datay[i])*datax[i,j]    # 2.2.1
    
    # 计算代价函数，判断迭代终止条件
        lb = 0
        for i in range(number):
            wtxi_cur = np.dot(datax[i,:],cur_beta.reshape(-1,1)) 
            # 即当前迭代后beta对number号样本的预测值
            lb += -datay[i]*wtxi_cur + np.log(1+np.exp(wtxi_cur))
        if np.abs(cur_lb-lb)<1e-5:
            return cur_beta
        old_beta = cur_beta
        cur_lb = lb
        # 不满足结束条件，继续进行迭代

b. 牛顿迭代法

思路：

• 对 $\beta$β 给定一个初始值，进行迭代

• t+1 轮迭代公式为：$\beta_{t+1}=\beta_t-\left(\frac{\partial^2l(\beta)}{\partial\beta\partial\beta^T}\right)^{-1}\frac{\partial l(\beta)}{\partial\beta} \tag{2.2.1}$βt+1​=βt​−(∂β∂βT∂2l(β)​)−1∂β∂l(β)​(2.2.1)

• 其中：关于$\beta$β的一阶导数为：$\frac{\partial l(\beta)}{\partial\beta} = -\sum^m_{i=1}\hat{x_i}(y_i-p_1(\hat{x_i};\beta)) \tag{2.2.2}$∂β∂l(β)​=−i=1∑m​xi​^​(yi​−p1​(xi​^​;β))(2.2.2)

  关于$\beta$β的二阶导数为：
  $\frac{\partial^2 l(\beta)}{\partial\beta\partial\beta^T}=\sum^m_{i=1}\hat{x_i}\hat{x_i}^Tp_1(\hat{x_i};\beta)(1-p_1(\hat{x_i};\beta)) \tag{2.2.3}$∂β∂βT∂2l(β)​=i=1∑m​xi​^​xi​^​Tp1​(xi​^​;β)(1−p1​(xi​^​;β))(2.2.3)

  其中的 $p_1$p1​ 表示：$p(y=1|x) = \frac{e^{w^Tx+b}}{1+e^{w^Tx+b}} \tag{2.2.4}$p(y=1∣x)=1+ewTx+bewTx+b​(2.2.4)

特点

• 很快，比梯度下降法快得多
• 注意以上 $\beta$β 的二阶导(黑塞矩阵 Hessian )不一定是可逆的, 迭代公式可能根本进行不下去（注意matlab的伪逆也没用
• 同理由于迭代步长固定为1，容易出现无法收敛的情况
• Hessian 矩阵不可导时，考虑拟牛顿法；无法收敛时考虑阻尼牛顿法

def newton_cg(old_beta):
    
    '''
    使用牛顿迭代法得到系数 beta
    param old_beta: 系数的初始值
    return cur_beta: 迭代结束得到的系数值
    '''
    
    beta_d1 = 0    # 一阶导
    beta_d2 = 0    # 二阶导（Hessian矩阵
    cur_lb = 0    # 代价损失函数（用于判断迭代是否终止
    
    while 1:
        
        for i in range(number):
            wtxi = np.dot(datax[i,:],old_beta.reshape(-1,1))
            # 即当前beta对number号样本的预测值
            p1 = np.exp(wtxi)/(1+np.exp(wtxi))   # 2.2.4
            beta_d1 -= datax[i,:].reshape(-1,1)*(datay[i]-p1)    # 2.2.2
            beta_d2 += np.dot(datax[i,:].reshape(-1,1),datax[i,:].reshape(1,-1))*p1*(1-p1)    # 2.2.3
        
        # 进行beta更新
        cur_beta = old_beta - np.dot(np.linalg.inv(beta_d2),beta_d1).reshape(1,-1)    # 2.2.1

        # 计算代价函数，判断迭代终止条件
        lb = 0
        for i in range(number):
            wtxi_cur = np.dot(datax[i,:],cur_beta.reshape(-1,1)) 
            # 即当前迭代后beta对number号样本的预测值
            lb += -datay[i]*wtxi_cur + np.log(1+np.exp(wtxi_cur))
        if np.abs(cur_lb-lb)<1e-5:
            return cur_beta
        old_beta = cur_beta
        cur_lb = lb
        # 不满足结束条件，继续进行迭代
            

c. 拟牛顿法（BFGS

• 考虑到梯度下降法很慢，而牛顿法存在黑塞矩阵不可逆的问题，考虑拟牛顿法
• 思路（BFGS）：通过迭代的方式，求得一个黑塞矩阵逆的近似矩阵，在每次进行目标值的迭代时，沿着迭代方向进行armijo搜索确定的步长的迭代，余下思路与牛顿法相同

a) 考虑使用另一不断迭代的矩阵表示黑塞矩阵逆的近似
使用 G_k 表示第k次迭代更新beta时的黑塞矩阵的逆的近似

• 则多次迭代中 G_k 的更新公式为：
  $G_{k+1} = (I- \frac{\delta _k y_k^T}{y_k^T\delta _k})G_k(I- \frac{y_k\delta _k^T}{y_k^T\delta _k}) + \frac{\delta _k\delta _k^T}{y_k^T\delta _k}\tag{2.3.1}$Gk+1​=(I−ykT​δk​δk​ykT​​)Gk​(I−ykT​δk​yk​δkT​​)+ykT​δk​δk​δkT​​(2.3.1)

  • 其中：G_0初始设为单位矩阵 I
  • $\delta _k$δk​ 即表示多次迭代之间的差值：$\delta _k = \beta_{k+1} - \beta_k\tag{2.3.2}$δk​=βk+1​−βk​(2.3.2)
  • y_k 表示多次迭代求一阶导的差值$y_k = g_{k+1}-g_k\tag{2.3.3}$yk​=gk+1​−gk​(2.3.3)

b) 注意具体迭代时与牛顿法存在差异：拟牛顿法在迭代时使用阻尼牛顿法的思路

• 即增加一个步长因子，每次beta的迭代是向着牛顿方向，以步长因子为步长迭代

• 此时可以缓解牛顿法迭代点列发散的情况，即牛顿法是步长因子为1时的特殊的阻尼牛顿法

• 此时迭代：$\beta_{k+1} = \beta_k - \lambda_kH_k^{-1}g_k\tag{2.3.3}$βk+1​=βk​−λk​Hk−1​gk​(2.3.3)

  • 其中 H_k 即二阶导（黑塞矩阵Hessian：$H_k^{-1} = G_k$Hk−1​=Gk​
  • g_k 表示一阶导：$g_k = \frac{\partial l(\beta)}{\partial\beta} = -\sum^m_{i=1}\hat{x_i}(y_i-p_1(\hat{x_i};\beta)) \tag{2.3.4}$gk​=∂β∂l(β)​=−i=1∑m​xi​^​(yi​−p1​(xi​^​;β))(2.3.4)
    • 其中：$p(y=1|x) = \frac{e^{w^Tx+b}}{1+e^{w^Tx+b}} \tag{2.3.5}$p(y=1∣x)=1+ewTx+bewTx+b​(2.3.5)

• 对于步长因子 λ：为一个迭代的过程（Armijo搜索

  • 考虑调参的思路 $\delta \in (0,1), \sigma \in (0,0.5)$δ∈(0,1),σ∈(0,0.5)
  • $\lambda = \delta^{m_k}$λ=δmk​, m从0开始迭代（m是整数，从0开始往上加），找到第一个符合下列式子的m_k即结束，根据参数 δ 和 σ 得到步长因子的值
    • 根据此时已有的m计算出步长因子，用当下步长因子对 beta 进行更新，更新后beta满足下式：$l(\beta_{new})<= l(\beta_{old})+\sigma g_k^T \beta_{new}\tag{2.3.6}$l(βnew​)<=l(βold​)+σgkT​βnew​(2.3.6)
      • l(β) 即似然函数的对数 → 最小化的目标代价函数

• BFGS法整体步骤
  • 用所有下标带k的先给出 beta_k+1 （得到的过程中步长因子需要迭代
  • 判断 此时的 l(b) 目标函数下降足够少, 即退出迭代
  • 计算 y_k
  • 迭代黑塞矩阵的近似值 G 矩阵
  • 重新迭代计算新的beta

在每次迭代beta的时候迭代一次黑塞矩阵
在每次迭代beta的时候都要完整迭代完步长因子

def BFGS(old_beta,sigma,delta):
    
    '''
    使用拟牛顿法得到更新迭代的系数值
    param old_beta : 初始系数值
    param sigma, delta : 步长因子迭代部分调参
    return cur_beta
    '''
    
    G = np.identity(num_a)    # 使用 G 逼近黑塞矩阵的逆，初始从单位矩阵开始
    I = np.identity(num_a)
    g_k = 0    # 一阶导
    cur_lb = 0
    for i in range(number):    # 计算迭代开始前用到的 g_k
            wtxi = np.dot(old_beta,datax[i,:].reshape(-1,1))
            p1 = np.exp(wtxi)/(1+np.exp(wtxi))
            g_k -= datax[i,:].reshape(-1,1)*(datay[i]-p1)
            
            
    while 1:
        cur_beta = old_beta-np.dot(G,g_k).reshape(1,-1)    # 默认步长为1时的迭代
        
        # 寻找最佳步长因子
        for m in range(30):
            lam = np.power(delta,m)
            beta_new = old_beta-lam*np.dot(G,g_k).reshape(1,-1)
            left = 0
            right = 0
            for i in range(number):    # armijo搜索
                wtxi_new = np.dot(beta_new,datax[i,:].reshape(-1,1))
                wtxi_old = np.dot(old_beta,datax[i,:].reshape(-1,1))
                left += -datay[i]*wtxi_new + np.log(1+np.exp(wtxi_new))
                right += -datay[i]*wtxi_old + np.log(1+np.exp(wtxi_old))
            if left<= right+sigma*np.dot(beta_new,g_k.reshape(-1,1)):
                cur_beta = beta_new
                break
       
        # 计算 g_k+1
        g_k1 = 0
        for i in range(number):
            wtxi = np.dot(cur_beta,datax[i,:].reshape(-1,1))
            p1 = np.exp(wtxi)/(1+np.exp(wtxi))
            g_k1 -= datax[i,:].reshape(-1,1)*(datay[i]-p1)
        # 判断迭代结束条件
        lb = 0
        for i in range(number):
            wtxi_cur = np.dot(datax[i,:],cur_beta.reshape(-1,1)) 
            # 即当前迭代后beta对number号样本的预测值
            lb += -datay[i]*wtxi_cur + np.log(1+np.exp(wtxi_cur))
        if np.abs(cur_lb-lb)<1e-5:
            return cur_beta
            
        # 更新黑塞矩阵的近似矩阵
        yk = g_k1-g_k
        sk = old_beta-cur_beta
        cur_lb = lb
        G = np.dot(np.dot((I-np.dot(sk.reshape(-1,1),yk.reshape(1,-1))/np.dot(yk.reshape(1,-1),\
                                                                                 sk.reshape(-1,1))),G),\
                   (I-np.dot(yk,sk.reshape(1,-1)/np.dot(yk.reshape(1,-1),sk.reshape(-1,1)))))+\
            np.dot(sk.reshape(-1,1),sk.reshape(1,-1))/np.dot(yk.reshape(1,-1),sk.reshape(-1,1))

2.2 Sigmoid 函数转化

考虑到此时计算得到的结果本身是不能用作判断的
考虑一个 Sigmoid 函数将计算出的预测值转化为 0-1 之间的数字

• Sigmoid 函数：$y = \frac{1}{1+e^{-z}}$y=1+e−z1​

• 此时不仅将计算出的值转化至0-1之间，同时利用Sigmoid函数的连续可微性

def predict(beta,test_data):
    '''
    结合sigmoid函数对测试数据进行预测
    param beta: 即计算出的系数集
    param test_data: 即需要进行预测的数据，每行为一个样本，ndarray形式，注意此时最后一列并没有置1
    return result: 一个长为待预测样本数的列表，记录分类结果 0/1
    '''
    
    num,num_a = test_data.shape()
    result = [0]*num
    for i in range(num):
        for j in range(num_a):
            result[i] += test_data[i,j]*beta[j]
        result[i] = 1/1+np.exp(-result[i])
        if result[i]>0.5:
            result[i] = 1
        else:
            result[i] = 0
            
    return result

3. 数据尝试

• 此时采用《机器学习》_周志华 中 西瓜数据集 3.0a 进行试验
• 此时的数据只有两个数值型属性

密度 含糖率 好瓜
0.697 0.460 1
0.774 0.376 1
0.634 0.264 1
0.608 0.318 1
0.556 0.215 1
0.403 0.237 1
0.481 0.149 1
0.437 0.211 1
0.666 0.091 0
0.243 0.267 0
0.245 0.057 0
0.343 0.099 0
0.639 0.161 0
0.657 0.198 0
0.360 0.370 0
0.593 0.042 0
0.719 0.103 0

# 预先载入与处理数据
# 此时使用西瓜数据集 3.0a

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import time 

data = pd.read_csv('data+wm3a.csv')
datax = np.array(data.iloc[:,1:].copy())    #注意第一列为编号
datax[:,-1]=1    # 处理为 x^
datay = np.array(data.iloc[:,-1].copy())

number,num_a = datax.shape
# number：样本数
# num_a: 属性个数+1

inibeta = np.zeros((1,num_a))
inibeta[0,-1]=1
# 初始化系数

alpha = 0.5
sigma = 0.3
delta = 0.5

# 梯度下降法
start = time.time()
old_beta = inibeta
beta_sag = sag(inibeta,alpha)
end = time.time()
print('beta_sag:',beta_sag)
print('sag cost time:',end-start)

# 牛顿迭代法
start = time.time()
old_beta = inibeta
beta_newton = newton_cg(inibeta)
end = time.time()
print('\nbeta_newton:',beta_newton)
print('newton cost time:',end-start)

# BFGS法
start = time.time()
old_beta = inibeta
beta_BFGS = BFGS(inibeta,sigma,delta)
end = time.time()
print('\nbeta_BFGS:',beta_BFGS)
print('BFGS cost time:',end-start)

输出结果

beta_sag: [[ 3.60304639 11.98140328 -4.5915188 ]]
sag cost time: 5.659311532974243

beta_newton: [[ 2.97700622 12.75249236 -4.36482012]]
newton cost time: 0.06464886665344238

beta_BFGS: [[ 4.14063382 11.7053009  -4.85396699]]
BFGS cost time: 0.8432471752166748

# 作图表示
for i in range(17):
    if i<8:
        plt.plot(datax[i,0],datax[i,1],'*r')
        continue
    plt.plot(datax[i,0],datax[i,1],'+b')

listx = np.arange(0,1,0.002)
listy = np.arange(0,0.6,0.0012)
X,Y = np.meshgrid(listx,listy)

def predict_draw(beta,x,y):
    result = 0
    result+= beta[0,0]*x+beta[0,1]*y+beta[0,2]
    return result

# 用不同计算结果展示一下
Z = predict_draw(beta_BFGS,X,Y)
cline = plt.contour(X,Y,Z,[0.5],colors = 'orange')
plt.clabel(cline)

plt.xlabel('density')
plt.ylabel('sugar ratio')

4. sklearn 实现

参考链接

4.1 参数介绍

• 调用:

from sklearn import linear_model
# linear_model.LogisticRegression

• 参数：

  • penalty=‘l2’ : 字符串‘l1’或‘l2’,默认‘l2’。

    • 用来指定惩罚的基准（正则化参数）
    • 只有‘l2’支持‘newton-cg’、‘sag’和‘lbfgs’这三种算法。
    • 选择‘l2’，solver参数可选‘liblinear’、‘newton-cg’、‘sag’和‘lbfgs’; 选择‘l1’只能用‘liblinear’算法
    
    

  • C=1.0 : C为正则化系数λ的倒数，必须为正数，默认为1

    • 和SVM中的C一样，值越小，代表正则化越强
  
  

  • solver=‘liblinear’ : 对逻辑回归损失函数的优化方法

    • a) liblinear：使用了开源的liblinear库实现，内部使用了坐标轴下降法来迭代优化损失函数。

    • b) lbfgs：拟牛顿法的一种，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。

    • c) newton-cg：也是牛顿法家族的一种，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。

    • d) sag：即随机平均梯度下降，是梯度下降法的变种

      • 和普通梯度下降法的区别是每次迭代仅仅用一部分的样本来计算梯度，适合于样本数据较多的情况
      
      

  • multi_class=‘ovr’ : 分类方式

    • 即对于多分类问题的分类方式
    • ovr → one-vs-rest(OvR)
    • multinomial → many-vs-many(MvM)
      • 注意选择了 multinomial 时不能使用 liblinear 进行优化
  
  
  • class_weight=None : 类型权重参数
    • 用于标示分类模型中各种类型的权重。默认不输入，即所有的分类的权重一样。
    • 选择‘balanced’自动根据y值计算类型权重。
    • 自己设置权重，格式：{class_label: weight}
      • 例如0,1分类的er’yuan二元模型，设置class_weight={0:0.9, 1:0.1}，这样类型0的权重为90%，而类型1的权重为10%
  
  

  • max_iter=100 : 算法收敛的最大迭代次数

  • tol=0.0001 : 迭代终止判据的误差范围

4.2 常用调用方法

• fit(X, y, sample_weight=None)

  • 拟合模型，用来训练LR分类器，其中X是训练样本，y是对应的标记向量
  • 返回对象，self
  
  

• predict(X)

  • 用来预测样本，也就是分类（0/1值），X是测试集。返回array
  
  

• score(X, y, sample_weight=None)

  • 返回给定测试集合的平均准确率（mean accuracy），浮点型数值。
  • 对于多个分类返回，则返回每个类别的准确率组成的哈希矩阵

from sklearn import linear_model
clf = linear_model.LogisticRegression(solver = 'lbfgs',tol = 1e-5, C = 1e5)
clf.fit(datax[:,:-1],datay)

Z = clf.predict(datax[:,:-1])

此时 Z 直接是分类结果 0/1
但这里得到的分类结果并不是很好：

[1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 0 0]

——————————————————————————
*存在不足欢迎指正

转载请保留作者信息

• 点赞
• 收藏
• 分享
• • 文章举报

Blacksheep103221 发布了1 篇原创文章 · 获赞 0 · 访问量 32 私信 关注

标签：LogisticRegression,old,迭代,回归,牛顿,beta,xi,np,2.2
来源： https://blog.csdn.net/weixin_44914063/article/details/104512292