Processing math: 0%

统计学习:逻辑回归与交叉熵损失(Pytorch实现)

1 Logistic 分布和对率回归

监督学习的模型可以是概率模型或非概率模型,由条件概率分布P(Y|\bm{X})或决 策函数(decision function)Y=f(\bm{X})表示,随具体学习方法而定。对具体的输入\bm{x}进行相应的输出预测并得到某个结果时,写作P(y|\bm{x})y=f(\bm{x})

我们这里的 Logistic 分类模型是概率模型,模型P(Y|\bm{X})表示给定随机向量\bm{X}下,输出类别Y的条件概率分布。这里我们只讨论二分类问题,后面我们介绍多层感知机的时候会介绍多分类问题,道理是一样的。

首先,我们先来介绍 Logistic 分布。

Logistic 的分布函数为:

F(x) = P(X<=x) = \frac{1}{1+e^{-(x-u)/\gamma}}

该分布函数的图像为:

超平面几何示意图

该函数以点(u,1/2)中心对称,在中心附近增长速度较快,在两端增长速度较慢。

我们后面对于未归一化的数z,采用\text{sigmoid}函数

σ(z) = \frac{1}{1 + \text{exp}(−z)}

将其“压”在(0, 1)之间。这样就可以使z归一化,一般我们使归一化后的值表示置信度(belief),可以理解成概率,但是与概率有着细微的差别,更体现 “属于××类的把握”的意思。

对于二分类,输出Y=0Y=1,我们将P(Y=0| \bm{x})P(Y=1|\bm{x})表示如下,也就定义了二项逻辑回归模型:

\begin{aligned} P(Y=1|\bm{x}) &= \sigma(\bm{w}^T\bm{x}+b) \\ &= \frac{1}{1+\text{exp}(-(\bm{w}^T\bm{x}+b))} \end{aligned}\space (1) \\ \begin{aligned} P(Y=0|\bm{x}) &= 1-\sigma(\bm{w}^T\bm{x}+b) \\ &= \frac{1}{1+\text{exp}(\bm{w}^T\bm{x}+b)} \end{aligned} \qquad (2)

现在考察 Logistic 回归模型的特点,一个事件的几率(odds)是指该事件发生的概率与不发生的概率的比值。如果事件发生的概率是p,那么该事件的几率是p/(1-p),该事件的对数几率(log odds,简称对率)或 logit 函数是

\text{logit}(p) = \text{log}\space p/(1-p)

对 Logistic 回归而言,由式子(1)和式子(2)得到:

\text{log}\frac{P(Y = 1|\bm{x})}{1 − P(Y = 1|\bm{x})} = \bm{w}^T\bm{x} + b

这玩意在统计学里面称之为“对率回归”,其实就是“Logistic regression 名称”的由来。这里的 Logistic 和“逻辑”没有任何关系,和对率才是有关系的。 可以看出,输出Y=1的对数几率是由输入\bm{x}的线性函数表示的模型,即 Logistic 回归模型。

2 Logistic 回归模型的参数估计

我们在《统计推断:极大似然估计、贝叶斯估计与方差偏差分解》)中说过,从概率模型的视角来看,机器模型学习的过程就是概率分布参数估计的过程,这里就是估计条件概率分布P(Y|\bm{X})的参数\bm{θ}。对于给定的训练数据集T=\{\{\bm{x}_1, y_1\},\{\bm{x}_2, y_2\}, \dots, \{\bm{x}_N, y_N\}\},其中,\bm{x}_i∈\mathbb{R}^ny_i∈\{0, 1\},可以应用极大似然估计法估计模型参数,从而得到Logistic回归模型。

我们设P(Y=1|x)=π(\bm{x})P(Y=0|\bm{x})=1-π(\bm{x}),很显然P(Y|\bm{x})服从一个伯努利分 布,不过这个伯努利分布很复杂,它的参数p是一个函数π(\bm{x})。我们这里采用一个抽象的观念,将函数π(\bm{x})看做一个整体,这样P(Y|\bm{x})服从二项分布,可以方便我们理解。我们的参数\bm{θ}=(w, b),不过我们还是采用之前的技巧,将权值向量和输入向量加以扩充,直接记做\bm{θ}=\bm{w}

这样,对于N个样本,定义似然函数:

L(\bm{w}|\bm{x}_1,\bm{x}_2,...,\bm{x}_N) = \Pi_{i=1}^N[\pi(\bm{x}_i)]^{y_i}[1-\pi(\bm{x}_i)]^{1-y_i}

因为函数比较复杂,我们采用数值优化进行求解。然而这个函数是非凸的, 如果直接用数值迭代算法作用于次函数,可能陷入局部最优解不能保证到收敛到全局最优解。我们为了将其变为负对数似然函数(想一想为什么要取负号),也就是一个凸函数,方便用梯度下降法、牛顿法进行求解:

-L(\bm{w}|\bm{x}_1,\bm{x}_2,...,\bm{x}_N) = -\sum_{i=1}^N[y_i\text{log}\space \pi(\bm{x}_i)+(1-y_i)\text{log}(1-\pi(\bm{x}_i))]

配凑 1/N转换为关于数据经验分布期望的无偏估计 ,即\mathbb{E}_{\bm{x}\sim \hat{p}_{data}}\text{log}\space p_{model}(\bm{x}_i| \bm{\theta})而不改变目标函数的优化方向。

-\frac{1}{N}\sum_{i=1}^N[y_i\text{log}\space \pi(\bm{x}_i)+(1-y_i)\text{log}(1-\pi(\bm{x}_i))]

上面这个式子,就是我们在深度学习里面常用的交叉熵损失函数(cross entropy loss function)的特殊情况。(交叉熵损失函数一般是多分类,这里是二分类)。后面我们会学习,在深度学习领域中,我们用f(\cdot)表示神经网络对输入\bm{x}加以的 映射(这里没有激活函数,是线性的映射),f(\bm{x})就是输出的条件概率分布P(Y|\bm{X}=\bm{x})。 设类别总数为K\bm{x}_i为第i个样本的特征向量(特征维度为D),f(\bm{x}_i)输出维度也为Kf(\bm{x}_i)_k表示P(y_k | \bm{x}_i)。因为是多分类,P(y_k | \bm{x}_i)的计算方法就不能通过\text{sigmoid}函数来得到了,取而代之是更通用的\text{softmax}函数。我们给定输入样本\bm{x},设其被预测为第k类的logit z_k=\sum_{j=1}^Dw_{kj}x_{j},其中\textbf{W} ∈\mathbb{R}^{K\times D}神经网络的权重矩阵D为特征向量维度,K为类别总数。则我们有:

f(\bm{x})_k = \text{softmax}(\bm{z})_k=\frac{\text{exp}(z_k)}{\sum_{k^{\prime}}\text{exp}(z_{k^{\prime}})}

我们规定第i个样本的训练标签\bm{y}_iK维one-hot向量。则交叉熵函数表述如下:

-\frac{1}{N}\sum_{i=1}^N\sum_{k=1}^K y_{ik}\text{log}f(\bm{x}_i)_k

因为\bm{y}_i为 one-hot 向量,只有一个维度为 1,设这个维度为c(即表示类别c),则交叉熵函数可以化简为:

-\frac{1}{N}\sum_{i=1}^N \text{log}f(\bm{x}_i)_c

这里f(\bm{x}_i)_c=\frac{\text{exp}(z_c)}{\sum_{k}\text{exp}(z_k)}

很明显,如果f(\bm{x}_i)_c越大,表示神经网络预测\bm{x}_i样本类别为第c类的概率越大。这也是目标函数优化的方向——即对于类别为c的样本\bm{x}_i,优化器尽量使样本\bm{x}_i被预测为第c类的logits(也即z_c)更大,而使其被预测为其它类别的logits(也即z_{k\neq c})更小。如果你对这点有疑惑,可以进一步将f(\bm{x}_i)_c写成\frac{1}{1 + \sum_{k \neq c}\text{exp}(z_{k} - z_c)},这样看的更明显一些,这就是糖水不等式\frac{b + c}{a + c} > \frac{b}{a} (a>b>0, c>0)的原理233。

关于sigmoid函数和softmax函数的联系和区别这里额外强调一下。首先softmax函数是sigmoid函数在多分类情况下的扩展。不过,即使在二分类的情况下,二者也有一些区别。比如给定一个正类样本\bm{x}(其训练标签y=1),设z = \bm{w}^T\bm{x}+b\in\mathbb{R},sigmoid函数会使得下列概率最大化:

P(Y=1|\bm{x}) = σ(z) = \frac{1}{1 + \text{exp}(−{z})}

这里本质上是在使得z最大化。

而对softmax函数而言,设\bm{z} = \bm{W}\bm{x}+b\in\mathbb{R}^2。则对于正类样本(其训练标签为一个one-hot向量\bm{y},满足y_1=1, y_2 = 0),softmax函数会使得下列概率最大化:

P(Y=1|\bm{x}) = \text{softmax}(\bm{z})_1=\frac{\text{exp}(z_1)}{\text{exp}(z_1) + \text{exp}(z_2)} =\frac{1}{1 + \text{exp}(z_2 - z_1)}

这里不但使得z_1最大化,而且还使得z_2最小化。
因此我们可以认为对softmax而言,参与其中的单元之间形成了竞争,这和神经科学中的现象很吻合。事实上,softmax的输出总满足和为1,所以一个单元的值增加必然对应着其他单元值的减少。这与皮质中相邻神经元间的侧抑制类似。在极端情况下,它变成了赢者通吃(winner-take-all) 的形式(其中一个输出接近1,其它的接近0)。

如下为使用梯度下降法求解二分类逻辑回归问题(并采用breast_cancer数据集进行测试,该数据集是二分类数据集):

from locale import normalize
import numpy as np
import random
import torch
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# batch_size表示单批次用于参数估计的样本个数
# y_pred大小为(batch_size, 1)
# y大小为(batch_size, ),为类别型变量
def cross_entropy(y_pred, y):
    y = y.reshape(-1, 1)
    return -(torch.mul(y, torch.log(y_pred)) + torch.mul(1-y, torch.log(1-y_pred))).sum()/y.shape[0]

# 前向函数
def logistic_f(X, w): 
    z = torch.matmul(X, w).reshape(-1, 1) 
    return 1/(torch.exp(-z) + 1)

# 之前实现的梯度下降法,做了一些小修改
def gradient_descent(X, w, y, n_iter, eta, loss_func, f):
    for i in range(1, n_iter+1):
        y_pred = f(X, w)
        #print(y_pred)
        loss_v = loss_func(y_pred, y)
        #print(loss_v)
        loss_v.backward() 
        with torch.no_grad(): 
            w -= eta*w.grad
        w.grad.zero_()  
    w_star = w.detach()
    return w_star 

# 迭代次数
n_iter = 600

# 学习率
eta = 10 # 因为我们前面loss除了样本总数,这里调大些

if __name__ == '__main__':

    X, y = load_breast_cancer(return_X_y=True)

    D = X.shape[1]
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=0)
    n_train, n_test = X_train.shape[0], X_test.shape[0]

    # 初始化权重(含偏置)
    w = torch.tensor(2 * np.random.ranf(size=D + 1) - 1, requires_grad=True)
    print(w)
    
    # 扩充1维
    X_train = torch.tensor(np.concatenate([X_train, np.ones([X_train.shape[0], 1])], axis=1))
    X_train = torch.nn.functional.normalize(X_train) # 先归一化,否则后面sigmoid部分输出无穷趋近于0,导致log()得-inf
    X_test = torch.tensor(np.concatenate([X_test, np.ones([X_test.shape[0], 1])], axis=1))
    X_test = torch.nn.functional.normalize(X_test) # 先归一化,否则后面sigmoid部分输出无穷趋近于0,导致log()得-inf

 
    y_train, y_test = torch.tensor(y_train), torch.tensor(y_test)

    # 调用梯度下降法对函数进行优化
    # 这里采用单次迭代对所有样本进行估计,后面我们会介绍小批量法减少时间复杂度
    w_star = gradient_descent(X_train, w, y_train, n_iter, eta, cross_entropy, logistic_f)
    print(w_star)
    
    y_pred = logistic_f(X_test, w_star)
    pred_label = torch.where(y_pred < 0.5, 0, 1)
    acc = accuracy_score(y_test, pred_label)
    print("accuracy: %f" % (acc))


注意输入数据要先归一化到(0, 1),否则后面|Xw + b|过大,以致Xw + b<0时使sigmoid函数的输出无穷趋近于0,最终导致log函数得-inf。我们来看一下算法的输出结果。初始权重向量、最终得到的权重向量和模型精度如下:

Initial w: tensor([-0.3109, -0.4637, -0.0326,  0.0081, -0.0108, -0.9996, -0.2989,  0.2311,
        -0.9197,  0.7227,  0.2692,  0.4146,  0.6281, -0.8207,  0.3153,  0.9207,
        -0.4418,  0.6558, -0.1684,  0.9044, -0.7847, -0.1003, -0.2369, -0.0608,
        -0.3722,  0.3982, -0.7417, -0.8493,  0.0326,  0.1718, -0.3145],
       dtype=torch.float64, requires_grad=True)
Final w: tensor([  3.0186,   5.7167,  19.8512,  18.5805,   0.0232,  -1.0004,  -0.3366,
          0.2130,  -0.8558,   0.7494,   0.2977,   0.9121,   0.7568,  -6.9429,
          0.3183,   0.9233,  -0.4391,   0.6570,  -0.1606,   0.9056,   2.5354,
          7.5881,  19.1409, -19.0107,  -0.3286,   0.3691,  -0.8191,  -0.8705,
          0.1200,   0.1985,   0.1274], dtype=torch.float64)
Final acc: 0.918129

可见模型收敛正常。

如下是多分类逻辑回归问题,我们需要将\text{sigmoid}函数修改为\text{softmax}函数,损失函数修改为更为通用的交叉熵函数,根据输出维度变化将权重向量修改为权重矩阵等。我们采用iris数据集进行测试,该数据集是一个3分类数据集。

from locale import normalize
import numpy as np
import random
import torch
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# batch_size表示单批次用于参数估计的样本个数
# n_feature为特征向量维度
# n_class为类别个数
# y_pred大小为(batch_size, n_class)
# y大小为(batch_size, ),为类别型变量
def cross_entropy(y_pred, y):
    # 这里y是创建新的对象,这里将y变为(batch_size. 1)形式
    y = y.reshape(-1, 1)
    return -torch.log(torch.gather(y_pred, 1, y)).sum()/ y_pred.shape[0]

# 前向函数,注意,这里的sigmoid改为多分类的softmax函数
def logistic_f(X, W): 
    z = torch.matmul(X, W)
    return torch.exp(z)/torch.exp(z).sum()

# 之前实现的梯度下降法,做了一些小修改
def gradient_descent(X, W, y, n_iter, eta, loss_func, f):
    for i in range(1, n_iter+1):
        y_pred = f(X, W)
        loss_v = loss_func(y_pred, y)
        loss_v.backward() 
        with torch.no_grad(): 
            W -= eta*W.grad
        W.grad.zero_()  
    W_star = W.detach()
    return W_star 

# 迭代次数
n_iter = 1200

# 学习率
eta = 1 # 因为我们前面loss除了样本总数,这里也要调大些

if __name__ == '__main__':

    X, y = load_iris(return_X_y=True)

    n_classes = y.max() - y.min() # 分类类别总数
    D = X.shape[1]
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=0)
    n_train, n_test = X_train.shape[0], X_test.shape[0]

    # 初始化权重(含偏置)
    W = torch.tensor(2 * np.random.ranf(size=(D + 1, n_classes)), requires_grad=True)
    print("Initial W: %s" % W)
    
    # 扩充1维
    X_train = torch.tensor(np.concatenate([X_train, np.ones([X_train.shape[0], 1])], axis=1))
    X_train = torch.nn.functional.normalize(X_train) # 先归一化,否则后面sigmoid输出无穷趋近于0,导致log()得-inf
    X_test = torch.tensor(np.concatenate([X_test, np.ones([X_test.shape[0], 1])], axis=1))
    X_test = torch.nn.functional.normalize(X_test) # 先归一化,否则后面sigmoid输出无穷趋近于0,导致log()得-inf

    print(y_test)
    y_train, y_test = torch.tensor(y_train), torch.tensor(y_test)

    # 调用梯度下降法对函数进行优化
    # 这里采用单次迭代对所有样本进行估计,后面我们会介绍小批量法减少时间复杂度
    W_star = gradient_descent(X_train, W, y_train, n_iter, eta, cross_entropy, logistic_f)
    print("Final w: %s" % W_star)
    
    y_pred = logistic_f(X_test, W_star)
    pred_label = torch.argmax(y_pred, axis=1)

    acc = accuracy_score(y_test, pred_label)
    print("Final acc: %f" % (acc))

同样需要注意输入数据要先归一化到(0, 1)。我们来看一下算法的输出结果。初始权重向量、最终得到的权重向量和模型精度如下(同样,权重最后一行为偏置向量\bm{b}。因为是多分类,每一个类别维度k都会对应一个偏置b_k):

Initial W: tensor([[0.2088, 0.7352, 1.3713],
        [0.6877, 0.1925, 0.7291],
        [0.4245, 1.7149, 1.8838],
        [0.5134, 1.7090, 0.8319],
        [1.8033, 1.8256, 0.1155]], dtype=torch.float64, requires_grad=True)
Final w: tensor([[ 3.8652,  4.2337, -2.7884],
        [ 5.9122, -3.3024, -3.2518],
        [-7.2260,  3.5287, 10.9562],
        [-3.4692, -1.2400,  7.3598],
        [ 3.5502,  3.5274, -2.3503]], dtype=torch.float64)
Final acc: 0.933333

可见模型收敛正常。

参考

posted @   orion-orion  阅读(1238)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 实操Deepseek接入个人知识库
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
点击右上角即可分享
微信分享提示