逻辑回归模型 Logistic Regression 详细推导(含Numpy与PyTorch实现)

珍妮的选择

内容概括

逻辑回归模型 (Logistic Regression, LR) 是一个二分类模型, 它假设数据服从 Bernoulli 分布(也称 0 - 1 分布), 采用 Sigmoid 函数 \sigma(x) 将线性回归 Y = X\theta 的结果约束在 (0, 1) 区间内, 以表示样本属于某一类的概率. 之后通过最大似然函数的方法, 对目标函数采用梯度下降法实现对模型参数 \theta 的更新. 本质上, LR 模型仍然是一个线性模型.

下面的内容主要是对 LR 进行推导, 按照两种思路:

使用代数法进行推导
采用矩阵法进行推导

前者较为繁琐, 而后者非常简洁! 完成推导之后, 再分别使用 Numpy 或者 PyTorch 实现 LR 模型.

对自己之前分析过的文章做一个简单的总结:

机器学习基础: LR / LibFM
特征交叉: DCN / PNN / DeepMCP / xDeepFM / FiBiNet / AFM
用户行为建模: DSIN / DMR / DMIN / DIN
多任务建模: MMOE
Graph 建模: GraphSage / GAT

LR 模型介绍

符号说明

在介绍 LR 模型之前, 先对本文用到的符号进行说明;

x_i\in \mathbb{R}^{n\times 1} 表示第 i 个样本
y_i\in \{0, 1\} 表示第 i 个样本对应的标签
X\in \mathbb{R}^{m\times n} 为由 m 个样本组成的矩阵
Y\in \mathbb{R}^{m\times 1} 为 X 对应的标签组成的矩阵
\theta\in\mathbb{R}^{n\times 1} 为 LR 模型的权重参数
E\in\mathbb{R}^{m\times 1} 为全 1 1 向量, 即 [1, 1, \ldots, 1]^T

Sigmoid 函数

Sigmoid 函数定义为:

\sigma(x) = \frac{1}{1 + \exp{(-x)}} 

其函数图像如下:

当 x\rightarrow +\infty 时, \sigma(x)\rightarrow 1 ; 而当 x\rightarrow -\infty 时, \sigma(x)\rightarrow 0 . 由于 \sigma(x) 的取值范围在 (0, 1) 区间内, 因此可以用来表示概率的大小.

另外一个关于 Sigmoid 函数的有用性质是, 对其求导的结果可以用它的输出值来表示, 即:

\sigma^\prime(x) = \sigma(x)\left(1 - \sigma(x)\right)

具体推导过程如下:

\begin{aligned} \sigma^\prime(x) &= \left(\frac{1}{1 + \exp{(-x)}}\right)^\prime \\ &= -\frac{1}{\left(1 + \exp{(-x)}\right)^2}\cdot\exp{(-x)}\cdot(-x)^\prime \\ &= \frac{1}{1 + \exp{(-x)}}\cdot\frac{\exp{(-x)}}{1 + \exp{(-x)}} \\ &= \sigma(x)\left(1 - \sigma(x)\right) \end{aligned}

LR 模型

对样本 x\in\mathbb{R}^{n\times 1} , 设其类别为 y , 线性回归模型的参数设为 \theta\in\mathbb{R}^{n\times 1} , 使用 Sigmoid 函数将线性模型的结果 x^T\theta 进行转换, 便得到了二元逻辑回归的一般形式:

h_\theta(x) = \frac{1}{1 + \exp{(-x^T\theta)}}

可以用 h_\theta(x) 表示分类的概率, 如果 h_\theta(x) > 0.5 , 那么可以认为样本 x 的类别为 y=1 ; 若 h_\theta(x) < 0.5 , 则认为样本 x 的类别为 y=0 . 如果 h_\theta(x) 刚好等于 0.5 , 即此时 x 刚好等于 0 , 模型无法判断样本的具体类别, 但具体实现时, 一般将等于 0.5 的情况加入到前面两种情况之一中.

将二元逻辑回归写成矩阵的形式:

h_\theta(X) = \frac{1}{1 + \exp{(-X\theta)}} 

其中 X = [x_1, x_2, \ldots, x_m]^T\in\mathbb{R}^{m\times n} 为样本的输入特征矩阵, 参数 \theta\in\mathbb{R}^{n\times 1} , 那么输出结果 h_\theta(X)\in\mathbb{R}^{m\times 1} .

LR 模型的优化目标

似然函数与损失函数

对于输入样本 (x, y) , 利用 LR 模型可以得到它属于某一类的概率分别为:

\begin{aligned} P(y = 1 | x, \theta) &= h_{\theta}(x) \\ P(y = 0 | x, \theta) &= 1 - h_{\theta}(x) \end{aligned}

将两个式子合并为一个式子, 可以表示如下:

P(y | x, \theta) = h_{\theta}(x)^y\left(1 - h_{\theta}(x) \right)^{(1 - y)}

对于样本集 \mathcal{T}=\{(x_i, y_i)\}_{i=1}^{m} , 其似然函数可以表示为:

J(\theta) = \prod_{i=1}^{m}h_{\theta}(x_i)^{y_i}\left(1 - h_{\theta}(x_i) \right)^{(1 - y_i)}

如果采用似然函数最大化进行优化, 求的是最大值, 如果取负, 相当于优化方向是进行最小化目标, 之所以取负, 是因为一般机器学习中我们的优化目标是最小化损失函数, 通常采用梯度下降法来求解, 原因是梯度的负方向函数值下降最快. 不取负也是可以的, 但是在更新模型参数的时候, 需要做简单的修改. 这里按照惯例来.

其对数似然函数(代数形式)取负为:

\begin{aligned} L(\theta) &=-\sum_{i=1}^{m}\left[y_{i} \log h_{\theta}\left(x_{i}\right)+\left(1-y_{i}\right) \log \left(1-h_{\theta}\left(x_{i}\right)\right)\right] \\ &=-\sum_{i=1}^{m}\left[y_{i} \log \frac{h_{\theta}\left(x_{i}\right)}{1-h_{\theta}\left(x_{i}\right)}+\log \left(1-h_{\theta}\left(x_{i}\right)\right)\right] \\ &= -\sum_{i=1}^{m}\left[y_{i} \log \frac{1}{\exp{(-x_i^T\theta)}}+\log \left(1-h_{\theta}\left(x_{i}\right)\right)\right] \\ &= -\sum_{i=1}^{m}\left[y_{i} \log \frac{1}{\exp{(-x_i^T\theta)}}+\log \left(1-\frac{1}{1 + \exp{(-x_i^T\theta)}}\right)\right] \\ &= -\sum_{i=1}^{m}\left[y_{i} \log \frac{1}{\exp{(-x_i^T\theta)}}+\log \left(\frac{\exp{(-x_i^T\theta)}}{1 + \exp{(-x_i^T\theta)}}\right)\right] \\ &= -\sum_{i=1}^{m}\left[y_{i} \log \frac{1}{\exp{(-x_i^T\theta)}}+\log \left(\frac{1}{1 + \exp{(x_i^T\theta)}}\right)\right] \\ &=-\sum_{i=1}^{m}\left[y_{i}\left(x_i^T\theta\right)-\log \left(1+\exp \left(x_i^T\theta\right)\right)\right] \end{aligned}

如果用矩阵来表示 L(\theta) , 那么结果为:

L(\theta) = -\left(Y^TX\theta - E^T\log\left(E + \exp{(X\theta)}\right)\right)

其中 E\in\mathbb{R}^{m\times 1} 为全 1 向量, 即 [1, 1, \ldots, 1]^T .

模型参数更新 – 代数法求梯度

这一小节使用代数法求二元逻辑回归的梯度, 相对繁琐; 而使用矩阵法求解在形式上更为简洁, 但是理解上有一定的门槛.

前面得到了损失函数为:

\begin{aligned} L(\theta) &=-\sum_{i=1}^{m}\left[y_{i}\left(x_i^T\theta\right)-\log \left(1+\exp \left(x_i^T\theta\right)\right)\right] \end{aligned}

那么 \frac{\partial L(\theta)}{\partial \theta_j} 的结果为:

\begin{aligned} \frac{\partial L(\theta)}{\partial \theta_j} &=-\frac{\partial}{\partial \theta_j}\sum_{i=1}^{m}\left[y_{i}\left(x_i^T\theta\right)-\log \left(1+\exp \left(x_i^T\theta\right)\right)\right] \\ &= -\sum_{i=1}^{m}\left[y_{i}x_{ij} - \frac{1}{1+\exp \left(x_i^T\theta\right)}\cdot \exp \left(x_i^T\theta\right)\cdot x_{ij} \right] \\ &= -\sum_{i=1}^{m}\left[y_{i} - \frac{1}{1+\exp \left(-x_i^T\theta\right)}\right] x_{ij} \\ &= -\sum_{i=1}^{m}\left[y_{i} - h_{\theta}(x_i)\right] x_{ij} \\ &= \sum_{i=1}^{m}\left[h_{\theta}(x_i) - y_{i}\right] x_{ij} \end{aligned}

那么 \frac{\partial L(\theta)}{\partial \theta} = \left[\frac{\partial L(\theta)}{\partial \theta_1}, \frac{\partial L(\theta)}{\partial \theta_2}, \ldots, \frac{\partial L(\theta)}{\partial \theta_n}\right]^T = X^T\left(h_\theta{(X)} - Y\right)

其中 \frac{\partial L(\theta)}{\partial \theta} \in\mathbb{R}^{n\times 1} ， X^T\in{n\times m} ， \left(h_\theta{(X)} - Y\right)\in\mathbb{R}^{m\times 1}

模型参数更新 – 矩阵法求梯度

这里采用矩阵法来求梯度, 根据前面得到的损失函数的矩阵形式为:

L(\theta) = -\left(Y^TX\theta - E^T\log\left(E + \exp{(X\theta)}\right)\right)

其中 E\in\mathbb{R}^{m\times 1} 为全 1 向量.

在求解之前, 需要了解关于矩阵导数与微分以及迹的关系, 详情可以参考: 矩阵求导术（上）

其中需要用到的是:

df = \sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}\frac{\partial f}{\partial X_{ij}}dX_{ij} = \text{tr}\left({\frac{df}{dX}}^TdX\right)
若 x 为向量, 那么 df = {\frac{df}{dx}}^Tdx
逐元素函数: \text{d}\sigma(X) = \sigma^\prime(X)\odot\text{d}X , 其中 \sigma(X) = [\sigma(X_{ij})] 是逐元素标量函数运算, \sigma^\prime(X) = [\sigma^\prime(X_{ij})] 是逐元素求导数.
矩阵乘法/逐元素乘法交换: \text{tr}(A^T(B\odot C)) = \text{tr}((A\odot B)^TC) , 其中 A, B, C 尺寸相同, 两边都等于 \sum\limits_{ij}A_{ij}B_{ij}C_{ij}

因此推导如下:

\begin{aligned} \text{d}L(\theta) &= -\left(Y^TX\theta - E^T\log\left(E + \exp{(X\theta)}\right)\right) \\ &= -\left(\text{d}Y^TX\theta + Y^T\text{d}X\theta + Y^TX\text{d}\theta - \text{d}E^T\log\left(E + \exp{(X\theta)}\right) - E^T\text{d}\log\left(E + \exp{(X\theta)}\right)\right) \\ &= -\left(Y^TX\text{d}\theta - E^T\text{d}\log\left(E + \exp{(X\theta)}\right)\right) \\ &= -\left(Y^TX\text{d}\theta - E^T\left(\frac{1}{E + \exp{(X\theta)}}\odot\text{d}\left(E + \exp{(X\theta)}\right)\right)\right) \\ &= -\left(Y^TX\text{d}\theta - \left(E\odot\frac{1}{E + \exp{(X\theta)}}\right)^T\text{d}\left(E + \exp{(X\theta)}\right)\right) \\ &= -\left(Y^TX\text{d}\theta - \left(E\odot\frac{1}{E + \exp{(X\theta)}}\right)^T\left(\exp{(X\theta)}\odot\text{d}\left(X\theta\right)\right)\right) \\ &= -\left(Y^TX\text{d}\theta - \left(E\odot\frac{1}{E + \exp{(X\theta)}}\odot\exp{(X\theta)}\right)^T\text{d}\left(X\theta\right)\right) \\ &= -\left(Y^TX\text{d}\theta - h_{\theta}(X)^T\text{d}\left(X\theta\right)\right) \\ &= -\left(Y^TX\text{d}\theta - h_{\theta}(X)^T\left(\text{d}X\theta + X\text{d}\theta\right)\right) \\ &= -\left(Y^TX\text{d}\theta - h_{\theta}(X)^TX\text{d}\theta\right) \\ &= -\left([Y^T- h_{\theta}(X)^T]X\text{d}\theta\right) \\ &= [h_{\theta}(X) - Y]^TX\text{d}\theta \\ \end{aligned}

由矩阵求导公式 2, 即若 x 为向量, 那么 df = {\frac{df}{dx}}^Tdx , 那么可以得到:

\frac{\partial L(\theta)}{\partial \theta} = X^T(h_{\theta}(X) - Y)

(吐槽: 打完这些公式也太累了吧… 另外我前面说矩阵法求更简洁, 是在没有打这些公式的情况下说的, 现在弄完这些公式, 感觉也不简洁 … )

参数更新

使用梯度下降法对 \theta 的更新公式为:

\begin{aligned} \theta &= \theta - \alpha\frac{\partial L(\theta)}{\partial \theta} \\ &= \theta - \alpha X^T(h_{\theta}(X) - Y) \end{aligned}

LR Numpy 代码实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from collections import Counter
def sigmoid(x):
    x = np.array(x)
    return 1. / (1. + np.exp(-x))
## 目标函数, 极大似然
## 注意这里求取了平均值而不是直接 sum
def L(w, b, X, y):
    dot = np.dot(X, w) + b
    return np.mean(y * dot - np.log(1 + np.exp(dot)), axis=0)
## w, b 的导数
def dL(w, b, X, y):
    dot = np.dot(X, w) + b
    distance = y - sigmoid(dot)
    distance = distance.reshape(-1, 1)
    return np.mean(distance * X, axis=0), np.mean(distance, axis=0)
## 随机梯度下降? (上升)
def sgd(w, b, X, y, epoch, lr):
    for i in range(epoch):
        dw, db = dL(w, b, X, y)
        w += lr * dw
        b += lr * db
    return w, b
## 测试代码, 对于预测值, 当概率大于 0.5 时, label 属于 True
def predict(w, b, X_test):
    return sigmoid(np.dot(X_test, w) + b) >= 0.5
## 画出分类面
def plot_surface(X, y, w, b):
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01), np.arange(y_min, y_max, 0.01))
    X_test = np.c_[xx.ravel(), yy.ravel()]
    Z = predict(w, b, X_test)
    Z = Z.reshape(xx.shape)
    fig, ax = plt.subplots()
    counter = Counter(y)
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)
    ax.set_xlabel(feature_names[




    
0])
    ax.set_ylabel(feature_names[1])
    ax.contourf(xx, yy, Z, cmap=plt.cm.RdYlBu)
    ## 画出分割线
    #     i = np.linspace(x_min, x_max, 100)
    #     o = (w[0] * i + b) / -w[1]
    #     ax.plot(i, o)
    for label in counter.keys():
        ax.scatter(X[y==label, 0], X[y==label, 1])
    plt.show()
## 训练代码
iris = load_iris()
X = iris.data[:100, :2]
y = iris.target[:100] # y \in {0, 1}
feature_names = iris.feature_names[2:]
np.random.seed(123)
n = X.shape[1]
w = np.random.randn(n)
b = np.random.randn(1)
print('initial: w: {}, b: {}, L: {}'.format(w, b, L(w, b, X, y)))
w, b = sgd(w, b, X, y, 10000, 0.001)
print('final: w: {}, b: {}, L: {}'.format(w, b, L(w, b, X, y)))
plot_surface(X, y, w, b)

效果:

LR 的 PyTorch 实现

import torch
import torch.nn as nn
import torch.optim as optim
from collections import Counter
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
iris = load_iris()
X = iris.data[:100, :2]
y = iris.target[:100] # y \in {0, 1}
X = torch.tensor(X).float()
y = torch.tensor(y).float()
feature_names = iris.feature_names[2:]
class LR(nn.Module):
    def __init__(self, in_features):
        super(LR, self).__init__()
        self.linear = nn.Linear(in_features, 1, bias=True)
    def sigmoid(self, x):
        return 1. / (1 + torch.exp(-x))
    def predict(self, x):
        return (self(x) > 0.5).int()
    def forward(self, x):
        x = self.linear(x)
        x = self.sigmoid(x)
        return x
model = LR(X.size(1))
criterion = nn.BCELoss()
optimizer = optim.SGD(model.parameters(), lr=0.01)
epoch = 100
batch_size = 10
N = X.size(0)
for i in range(epoch):
    order = torch.randperm(N)
    X = X[order]
    y = y[order]
    for n in range(N // batch_size):
        input = X[n * batch_size : (n + 1) * batch_size]
        label = y[n * batch_size : (n + 1) * batch_size]
        optimizer.zero_grad()
        output = model(input)
        loss = criterion(output, label)
        loss.backward()
        optimizer.step()
        if n % 100 == 0:
            print(loss) 
def plot_surface(X, y, w, b):
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01), np.arange(y_min, y_max, 0.01))
    X_test = np.c_[xx.ravel(), yy.ravel()]
    Z = predict(w, b, X_test)
    Z = Z.reshape(xx.shape)
    fig, ax = plt.subplots()
    counter = Counter(y)
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)
    ax.set_xlabel(feature_names[0])
    ax.set_ylabel(feature_names[1])
    ax.contourf(xx, yy, Z, cmap=plt.cm.RdYlBu)
    ## 画出分割线
#     i = np.linspace(x_min, x_max, 100)
#     o = (w[0] * i + b) / -w[1]
#     ax.plot(i, o)
    for label in counter.keys():
        ax.scatter(X[y==label, 0], X[y==label, 1])
    plt.show()
def sigmoid(x):
    x = np.array(x)
    return 1. / (1. + np.exp(-x))
def predict(w, b, X_test):
    return sigmoid(np.dot(X_test, w) + b) >= 0.5
X = torch.tensor(X).float()
y = torch.tensor(y).float()