梯度下降算法在逻辑回归中的应用

逻辑回归简介

• sigmoid函数：

$$g(z)=\frac{1}{1+e^{-z}}$$

• 逻辑回归假设函数：

$$\hat{y}=h_{\theta }(x)=g({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$

其中${\theta }^T$和$x$延续在梯度下降算法在线性回归中的应用一文中的表达，即：
$$\begin{gathered} {\theta }^T=\left[\begin{array}{c}{\theta }_0,{\theta }_1,{\theta }_2,\cdots ,{\theta }_n\end{array}\right] \\ x=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ ⋮\\ x_n\end{array}\right] \end{gathered}$$

• 损失函数：

$$L(\hat{y},y)=-ylog(\hat{y})-(1-y)log(1-\hat{y})$$

• 损失函数可以这样理解，当某一样本的真实标签$y=1$（正类）时，损失函数变为：

$$L(\hat{y},y)=-log(\hat{y})$$

• 我们希望损失函数越小越好，那么也就是说希望$\hat{y}$越大越好，而$\hat{y}$表示是的概率，定义域是$[0,1]$，$\hat{y}$最大只能是1，与$y=1$相同。
• 同样的，当某一样本的真实标签$y=0$（负类）时，损失函数变为：

$$L(\hat{y},y)=log(1-\hat{y})$$

我们希望损失函数越小越好，那么也就是说希望$\hat{y}$越小越好，而$\hat{y}$表示是的概率，定义域是$[0,1]$，$\hat{y}$最小只能是0，与$y=0$相同。

• 代价函数：

$$\begin{array}{rl}J(\theta )& =\frac{1}{m}\sum _{i=1}^{m}L({\hat{y}}^{(i)},y^{(i)})\\ & =\frac{1}{m}\sum _{i=1}^m[-y^{(i)}log({\hat{y}}^{(i)})-(1-y^{(i)})log(1-{\hat{y}}^{(i)})]\\ & =\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}log\left(\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\right)-(1-y^{(i)})log\left(\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\right)\right]\end{array}$$

梯度下降算法求解逻辑回归

• $J(\theta )$对${\theta }_j$求偏导，${\theta }_j$表示第$j$个$\theta$的值：

$$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}J(\theta )& =\frac{\partial }{\partial {\theta }_j}\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}log\left(\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\right)-(1-y^{(i)})log\left(\frac{1}{1+e^{{\theta }^T{x}^{(i)}}}\right)\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}log\left(1+e^{-{\theta }^T{x}^{(i)}}\right)+(1-y^{(i)})log\left(1+e^{{\theta }^T{x}^{(i)}}\right)\right]\end{array}$$

• 对$log\left(1+e^{{\theta }^T{x}^{(i)}}\right)$​求偏导，注意这里采用复合函数求导法则

$$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}log\left(1+e^{-{\theta }^T{x}^{(i)}}\right)& =\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\times e^{-{\theta }^T{x}^{(i)}}\times -x_j^{(i)}\\ & =\frac{-x_j^{(i)}{e}^{-{\theta }^T{x}^{(i)}}}{1+e^{-{\theta }^T{x}^{(i)}}}\end{array}$$

• 则，原式为：

$$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}J(\theta )& =\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\times \frac{-x_j^{(i)}{e}^{-{\theta }^T{x}^{(i)}}}{1+e^{-{\theta }^T{x}^{(i)}}}+(1-y^{(i)})\times \frac{x_j^{(i)}{e}^{{\theta }^T{x}^{(i)}}}{1+e^{{\theta }^T{x}^{(i)}}}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\times \frac{-x_j^{(i)}}{e^{{\theta }^T{x}^{(i)}}+1}+(1-y^{(i)})\times \frac{x_j^{(i)}{e}^{{\theta }^T{x}^{(i)}}}{1+e^{{\theta }^T{x}^{(i)}}}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[\frac{-x_j^{(i)}{y}^{(i)}+x_j^{(i)}{e}^{{\theta }^T{x}^{(i)}}-y^{(i)}{x}_j^{(i)}{e}^{{\theta }^T{x}^{(i)}}}{1+e^{{\theta }^T{x}^{(i)}}}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[\frac{-y^{(i)}(1+e^{{\theta }^T{x}^{(i)}})+e^{{\theta }^T{x}^{(i)}}}{1+e^{{\theta }^T{x}^{(i)}}}{x}_j^{(i)}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[\left(-y^{(i)}+\frac{e^{{\theta }^T{x}^{(i)}}}{1+e^{{\theta }^T{x}^{(i)}}}\right)x_j^{(i)}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[\left(-y^{(i)}+\frac{1}{e^{-{\theta }^T{x}^{(i)}}+1}\right)x_j^{(i)}\right]\\ & =\frac{1}{m}\sum _{i=1}^m\left[\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)x_j^{(i)}\right]\end{array}$$

根据梯度下降算法公式，迭代公式为：
$${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m\left[\left(h_{\theta }(x^{(i)})-y^{(i)}\right)x_j^{(i)}\right]$$
上面的$:=$符号表示先算右边的式子，算完后赋值给左边的变量。注意：左边和右边的${\theta }_i$值并不相等，右边是迭代前的${\theta }_i$值，左边的迭代后${\theta }_i$的值。

python代码实现

• 导入必要包

import numpy as np
import statsmodels.api as sm
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing, make_regression

• 读取数据函数

def load_data(data_path):
    data = np.loadtxt(data_path,delimiter=',')
    n = data.shape[1]-1
    data_x = data[:,0:n]
    data_y = data[:,-1].reshape(-1,1)
    return data_x,data_y

• sigmoid函数

def sigmoid(z):
    r = 1/(1 + np.exp(-z))
    return r

• 逻辑回归函数

def logic_lr(data_x, theta):
    z = np.dot(data_x,theta)
    return sigmoid(z)

• 代价函数

def compute_loss(data_x, data_y, theta):
    row_num, col_num = data_x.shape
    l = -1 * data_y * np.log(logic_lr(data_x, theta)) - (1-data_y)*np.log(1-logic_lr(data_x, theta))
    return np.sum(l)/row_num

• 梯度下降算法求解逻辑回归函数

def solve_logic_lr(data_x, data_y, theta, alpha, steps):
    temp_ones = np.ones(data_x.shape[0]).transpose()
    data_x = np.insert(data_x, 0, values=temp_ones, axis=1)
    row_num, col_num = data_x.shape
    loss_list = []
    for step in range(steps):
        loss_list.append(compute_loss(data_x, data_y, theta))
        for i in range(col_num):
            theta[i] = theta[i] - (alpha/row_num) * np.sum((logic_lr(data_x, theta) - data_y) * data_x[:,i].reshape(-1,1))
    return theta, loss_list

• 预测函数

def predict(data_x, theta):
    temp_ones = np.ones(data_x.shape[0]).transpose()
    data_x = np.insert(data_x, 0, values=temp_ones, axis=1)
    p = logic_lr(data_x, theta)
    p[p >= 0.5] = 1
    p[p < 0.5] = 0
    return p

• 函数调用

data_x, data_y = load_data('/kaggle/input/studyml/cls.txt')
data_x, mu, sigma = data_std(data_x)
# 变成列向量
theta = np.zeros(data_x.shape[1]+1).reshape(-1,1)
steps = 100
alpha = 0.0001
theta, loss = solve_logic_lr(data_x, data_y, theta, alpha, steps)
print(theta)
# 打印输出：[0.0009, 0.0027, 0.00249]

• 绘制Loss曲线（左）和拟合效果图（右）

plt.figure(figsize=(12,5),dpi=600)
plt.subplot(1,2,1)
plt.plot(loss)
plt.title("Loss Curve")
plt.xlabel("steps")
plt.ylabel("loss")
plt.subplot(1,2,2)
plt.scatter(data_x[:,0],data_x[:,1],c=data_y)
temp_x1 = np.arange(min(data_x[:,0]), max(data_x[:,0]), 0.1)
temp_x2 = -(theta[1] * temp_x1 + theta[0])/theta[2]
plt.plot(temp_x1, temp_x2)
plt.title("Fitting Effect Diagram")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

带L2正则化的逻辑回归

• 引入L2正则化，用于惩罚大的回归系数，减轻过拟合。则代价函数变为：

$$\begin{array}{rl}J(\theta )& =\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}log\left(\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\right)-(1-y^{(i)})log\left(\frac{1}{1+e^{-{\theta }^T{x}^{(i)}}}\right)\right]+\frac{\lambda }{2m}\Vert \theta {\Vert }_2^2\end{array}$$

• 根据上述的推导，可以得到$J(\theta )$对${\theta }_j$的偏导数为：

$$\frac{\partial }{\partial {\theta }_j}J(\theta )=\frac{1}{m}\sum _{i=1}^m\left[\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)x_j^{(i)}\right]+\frac{1}{m}\times \lambda {\theta }_j$$

根据梯度下降算法公式，迭代公式为：
$${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m\left[\left(h_{\theta }(x^{(i)})-y^{(i)}\right)x_j^{(i)}+\lambda {\theta }_j\right]$$

python代码实现

• 在上面的基础上，只需要修改代价函数和梯度下降求解函数就可以了
• 代价函数变为：

def compute_loss(data_x, data_y, theta, lambd):
    row_num, col_num = data_x.shape
    l = -1 * data_y * np.log(logic_lr(data_x, theta)) - (1-data_y)*np.log(1-logic_lr(data_x, theta))
    return np.sum(l)/row_num + (lambd/2*row_num) + np.sum(np.power(theta,2))

• 梯度下降求解函数变为：

def solve_logic_lr(data_x, data_y, theta, alpha, steps, lambd=0.01):
    temp_ones = np.ones(data_x.shape[0]).transpose()
    data_x = np.insert(data_x, 0, values=temp_ones, axis=1)
    row_num, col_num = data_x.shape
    loss_list = []
    for step in range(steps):
        loss_list.append(compute_loss(data_x, data_y, theta, lambd))
        for i in range(col_num):
            theta[i] = theta[i] - (alpha/row_num) * np.sum((logic_lr(data_x, theta) - data_y) * data_x[:,i].reshape(-1,1) + lambd * theta[i])
    return theta, loss_list

• 函数调用

data_x, data_y = load_data('/kaggle/input/studyml/cls.txt')
data_x, mu, sigma = data_std(data_x)
# 变成列向量
theta = np.zeros(data_x.shape[1]+1).reshape(-1,1)
steps = 100
alpha = 0.0001
theta, loss = solve_logic_lr(data_x, data_y, theta, alpha, steps)
print(theta)

• 绘制Loss曲线（左）和拟合效果图（右）

plt.figure(figsize=(12,5),dpi=600)
plt.subplot(1,2,1)
plt.plot(loss)
plt.title("Loss Curve")
plt.xlabel("steps")
plt.ylabel("loss")
plt.subplot(1,2,2)
plt.scatter(data_x[:,0],data_x[:,1],c=data_y)
temp_x1 = np.arange(min(data_x[:,0]), max(data_x[:,0]), 0.1)
temp_x2 = -(theta[1] * temp_x1 + theta[0])/theta[2]
plt.plot(temp_x1, temp_x2)
plt.title("Fitting Effect Diagram")
plt.xlabel("x")
plt.ylabel("y")
plt.show()