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

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

内容概括

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

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

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

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

广而告之

可以在微信中搜索 “珍妮的算法之路” 或者 “world4458” 关注我的微信公众号；另外可以看看知乎专栏 PoorMemory-机器学习, 以后文章也会发在知乎专栏中；

LR 模型介绍

符号说明

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

• $x_i\in {\mathbb{R}}^{n\times 1}$ 表示第 $i$ 个样本
• y i ∈ { 0 , 1 } y_i\in \{0, 1\} yi​∈{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,\dots ,1{]}^T$

Sigmoid 函数

Sigmoid 函数定义为:

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

其函数图像如下:

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

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

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

具体推导过程如下:

$$\begin{array}{rl}{\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{array}$$

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,\dots ,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{array}{rl}P(y=1|x,\theta )& =h_{\theta }(x)\\ P(y=0|x,\theta )& =1-h_{\theta }(x)\end{array}$$

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

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

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

$$J(\theta )=\prod _{i=1}^m{h}_{\theta }(x_i{)}^{y_i}{\left(1-h_{\theta }(x_i)\right)}^{(1-y_i)}$$

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

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

$$\begin{array}{rl}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{array}$$

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

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

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

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

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

前面得到了损失函数为:

$$\begin{array}{rl}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{array}$$

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

$$\begin{array}{rl}\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{array}$$

那么 $\frac{\partial L(\theta )}{\partial \theta }={\left[\frac{\partial L(\theta )}{\partial {\theta }_1},\frac{\partial L(\theta )}{\partial {\theta }_2},\dots ,\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^{T}X\theta -E^T\log \left(E+\exp (X\theta )\right)\right)$$

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

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

其中需要用到的是:

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

因此推导如下:

$$\begin{array}{rl}\text{d}L(\theta )& =-\left(Y^{T}X\theta -E^T\log \left(E+\exp (X\theta )\right)\right)\\ & =-\left(\text{d}{Y}^{T}X\theta +Y^T\text{d}X\theta +Y^{T}X\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^{T}X\text{d}\theta -E^T\text{d}\log \left(E+\exp (X\theta )\right)\right)\\ & =-\left(Y^{T}X\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^{T}X\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^{T}X\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^{T}X\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^{T}X\text{d}\theta -h_{\theta }(X{)}^T\text{d}\left(X\theta \right)\right)\\ & =-\left(Y^{T}X\text{d}\theta -h_{\theta }(X{)}^T\left(\text{d}X\theta +X\text{d}\theta \right)\right)\\ & =-\left(Y^{T}X\text{d}\theta -h_{\theta }(X{)}^{T}X\text{d}\theta \right)\\ & =-\left([Y^T-h_{\theta }(X{)}^T]X\text{d}\theta \right)\\ & =[h_{\theta }(X)-Y{]}^{T}X\text{d}\theta \end{array}$$

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

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

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

参数更新

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

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

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()

x = X.numpy()
t = y.numpy()
w = model.linear.weight.data.numpy().transpose()
b = model.linear.bias.data.numpy()[0]

plot_surface(x, t, w, b)

效果:

参考文献

• 逻辑回归原理小结: 刘建平Pinard 的博客, 大佬, 博文让人受益匪浅
• 矩阵求导术（上）: 感觉这是机器学习里的内功秘籍啊, 分为上下两卷, 我今天研读了上卷, 那种心情怎么来描述呢, 大彻大悟 ? No, 还没到出家的地步 😂😂😂. 总之, 仿佛脑袋中有个灯泡💡突然亮起来了的感觉.