统计学习：逻辑回归与交叉熵损失（Pytorch实现）
在Logistic 回归模型中，一个事件的几率(odds)是指该事件发生的概率与不发生的概率的比值。如果事件发生的概率是p，那么该事件的几率是p/(1-p)，该事件的对数几率(log odds，简称对率)或 logit 函数是log(p/1-p)。这玩意在统计学里面称之为“对率回归”，其实就是“Logistic regression 名称”的由来。这里的 Logistic 和“逻辑”没有任何关系，和对率才是有关系的。 可以看出，输出Y=1的对数几率是由输入x的线性函数表示的模型，即 Logistic回归模型。
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 的分布函数为：
该分布函数的图像为：
该函数以点\((u,1/2)\)中心对称，在中心附近增长速度较快，在两端增长速度较慢。
我们后面对于未归一化的数\(z\)，采用\(\text{sigmoid}\)函数
将其“压”在(0, 1)之间。这样就可以使\(z\)归一化，一般我们使归一化后的值表示置信度（belief），可以理解成概率，但是与概率有着细微的差别，更体现 “属于××类的把握”的意思。
对于二分类，输出\(Y=0\)或\(Y=1\)，我们将\(P(Y=0| \bm{x})\)和\(P(Y=1|\bm{x})\)表示如下，也就定义了二项逻辑回归模型：
现在考察 Logistic 回归模型的特点，一个事件的几率(odds)是指该事件发生的概率与不发生的概率的比值。如果事件发生的概率是\(p\)，那么该事件的几率是\(p/(1-p)\)，该事件的对数几率(log odds，简称对率)或 logit 函数是
对 Logistic 回归而言，由式子\((1)\)和式子\((2)\)得到:
这玩意在统计学里面称之为“对率回归”，其实就是“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}^n\)，\(y_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\)个样本，定义似然函数：
因为函数比较复杂，我们采用数值优化进行求解。然而这个函数是非凸的， 如果直接用数值迭代算法作用于次函数，可能陷入局部最优解不能保证到收敛到全局最优解。我们为了将其变为负对数似然函数（想一想为什么要取负号），也就是一个凸函数，方便用梯度下降法、牛顿法进行求解：
配凑 \(1/N\)转换为关于数据经验分布期望的无偏估计 ，即\(\mathbb{E}_{\bm{x}\sim \hat{p}_{data}}\text{log}\space p_{model}(\bm{x}_i| \bm{\theta})\)而不改变目标函数的优化方向。
上面这个式子，就是我们在深度学习里面常用的交叉熵损失函数(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)\)输出维度也为\(K\)。\(f(\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\)为类别总数。则我们有：
我们规定第\(i\)个样本的训练标签\(\bm{y}_i\)为\(K\)维one-hot向量。则交叉熵函数表述如下：
因为\(\bm{y}_i\)为 one-hot 向量，只有一个维度为 1，设这个维度为\(c\)（即表示类别\(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
可见模型收敛正常。
参考
- [1] 李航. 统计学习方法(第2版)[M]. 清华大学出版社, 2019.
- [2] 周志华. 机器学习[M]. 清华大学出版社, 2016.
- [3] Goodfellow I, Bengio Y, Courville A. Deep learning[M]. MIT press, 2016.
- [4] Pytorch官网： torch.nn.functional.normalize 函数
- [5] Scikit Learn 文档： Toy datasets
