SoftMax回归详解

SoftMax 回归

• 在观看此文章之前，推荐先去仔细了解 Logistic 回归：Logistic 回归详解

推导

数据说明

• 数据集： { ( x i , y i ) } i = 1 m \{(x_i, y_i)\}_{i=1}^m {(xi​,yi​)}i=1m​（$m$ 个样本）
• $x_i$ 可能是多维（假设有 $d$ 维）
• $y_i$ 是类别 { 0 , 1 , 2 , . . . , k } \{0,1,2,...,k\} {0,1,2,...,k}（分为 $k+1$ 类）的 $OneHot$ 编码
• $y_i=j$（属于第 $j$ 类）简称 $y_i^j$，且 $y_i^j=[0,0,...,1,...,0]$（下标 $j$ 的位置为 $1$）

本质

• 使用归一化 $SoftMax$ 函数：$p=\frac{e^{-({\omega }^{T}x+b)}}{{\sum }_{j=0}^k{e}^{-({\omega }_j^{T}x+b_j)}}=\frac{e^{-y}}{{\sum }_{j=0}^k{e}^{-y^j}}$，将线性回归 $y={\omega }^{T}x+b$ 与多分类任务 { 0 , 1 , 2 , . . . , k } \{0,1,2,...,k\} {0,1,2,...,k} 联系起来，得到了 $p=\frac{e^{-y}}{{\sum }_{j=0}^k{e}^{-y^j}}$
• 通过梯度下降法得到参数 $\omega$ 和 $b$（每一类对应一组参数，总共 $k+1$ 组 $\omega$ 和 $b$），代入要分类的样本 $x$，计算出对应的 $y$ 值（相当于该样本的类后验概率 $p(y|x)$，简称属于第 $j$ 类的后验概率 $p(y=j|x)$ 为 $p^j$），并与多个阈值比较，完成分类

与 Logistic 回归的关系（重点）

关系

• Logistic 回归是 SoftMax 回归的特殊形式
• 若 y ∈ { 0 , 1 } y\in\{0,1\} y∈{0,1}，则 $SoftMax$ 函数计算类后验概率 $p(y=0|x)=\frac{e^{-y^0}}{e^{-y^0}+e^{-y^1}}=\frac{1}{1+e^{-y^1}}$，即 Logistic 回归的形式

求导的关系（重点）

• 便于推导，令 $\theta =(b,\omega )$（ $d+1$ 维），对应的 $x=(1,x)$，（也就是增广）
• 在 Logistic 回归中，我们实际是对 $p^0=\frac{1}{1+e^{-y^1}}=\frac{e^{-y^0}}{e^{-y^0}+e^{-y^1}}$ 中的 $\theta =({\theta }_0,{\theta }_1)$ 求导，即对每个类对应的一组 $\theta$ 求导，结果再求和，如下
• $\frac{\partial p^0}{\partial \theta }=\frac{\partial p^0}{\partial {\theta }_0}+\frac{\partial p^0}{\partial {\theta }_1}=p^0(1-p^0)+0=p^0(1-p^0)$
• 于是求导的情况就分为两类：
  • 要计算类后验概率 $p^a$，对类 $a$ 对应的参数 ${\theta }_a$ 求导（如上面的 $\frac{\partial p^0}{\partial {\theta }_0}$）
  • 要计算类后验概率 $p^a$，却对类 $b$ 对应的参数 ${\theta }_b$ 求导（如上面的 $\frac{\partial p^0}{\partial {\theta }_1}$）

损失函数

• 由于预测结果是多类概率，这里 $p=[p^0,p^1,...,p^k]$
• 交叉熵损失函数：$J(\theta )=-\frac{1}{m}{\sum }_{i=1}^m{\sum }_{j=0}^k[y_i^{j}ln{p}_i]$

梯度下降法求参数 $\omega$, $b$

• $SoftMax$ 函数求导性质：对于 $p^a=\frac{e^{-y^a}}{{\sum }_{j=0}^k{e}^{-y^j}}$，有$\frac{\partial p^a}{\partial y}=\{\begin{array}{rlr}{p}^a(1-p^a),& & y=y^a\\ -p^a{p}^b,& & y=y^b\end{array}$（对应上面求导的两类情况）
• 求偏导：$\frac{\partial J(\theta )}{\partial \theta }=-\frac{1}{m}{\sum }_{i=1}^m[\frac{y_i{p}_i^a(1-p_i^a)x_i}{p_i^a}+{\sum }_{j=0,j\ne a}^k\frac{y_i^j(-p_i^a{p}_i^j)x_i}{p_i^j}]$
• 继续：$=-\frac{1}{m}{\sum }_{i=1}^m[y_i^a(1-p_i^a)+{\sum }_{j=0,j\ne a}^k{y}_i^j(-p_i^a)]x_i$
• 继续：$=-\frac{1}{m}{\sum }_{i=1}^m[y_i^a{\sum }_{j=0,j\ne a}^k{p}_i^j-p_i^a{\sum }_{j=0,j\ne a}^k{y}_i^j]x_i$ ，（注意 $1-p_i^a={\sum }_{j=0,j\ne a}^k{p}_i^j$）
• 得到：$=\frac{1}{m}{\sum }_{i=1}^m[p_i-y_i^j]x_i$（与 Logistic 完全一样）

实现

• 在 jupyter notebook 中代码均可运行

读取鸢尾花数据集

from sklearn import datasets
import numpy as np

iris = datasets.load_iris()    #加载鸢尾花数据集

X = iris["data"][:, (2, 3)]     #花瓣的长，宽
y = iris["target"]

X_with_bias = np.c_[np.ones([len(X), 1]), X]    #增广矩阵

切分数据集

test_ratio = 0.2
validation_ratio = 0.2
total_size = len(X_with_bias)

test_size = int(total_size * test_ratio)
validation_size = int(total_size * validation_ratio)
train_size = total_size - test_size - validation_size

rnd_indices = np.random.permutation(total_size)

X_train = X_with_bias[rnd_indices[:train_size]]
y_train = y[rnd_indices[:train_size]]
X_valid = X_with_bias[rnd_indices[train_size:-test_size]]
y_valid = y[rnd_indices[train_size:-test_size]]
X_test = X_with_bias[rnd_indices[-test_size:]]
y_test = y[rnd_indices[-test_size:]]

$OneHot$ 编码

def to_one_hot(y):    #类索引转换为矩阵，例如 0,1,2 类中属于 1 类，则 [0,1,0]
    n_classes = y.max() + 1
    m = len(y)
    Y_one_hot = np.zeros((m, n_classes))    #全 0 向量
    Y_one_hot[np.arange(m), y] = 1
    return Y_one_hot

Y_train_one_hot = to_one_hot(y_train)    #转换
Y_valid_one_hot = to_one_hot(y_valid)
Y_test_one_hot = to_one_hot(y_test)

$SoftMax$ 函数

def softmax(logits):    #与 Logistic 函数完全一样
    exps = np.exp(logits)
    exp_sums = np.sum(exps, axis=1, keepdims=True)
    return exps / exp_sums

梯度下降法求参数

eta = 0.01    #学习率
n_iterations = 5001    #迭代次数
m = len(X_train)    #训练样本数
epsilon = 1e-7    #平滑，为了防止 log(0) 出错

Theta = np.random.randn(n_inputs, n_outputs)    #随便生成 3 个 3 维的 theta

for iteration in range(n_iterations):
    logits = X_train.dot(Theta)
    Y_proba = softmax(logits)
    
    loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))    #计算损失
    
    error = Y_proba - Y_train_one_hot
    
    # if iteration % 500 == 0:    #输出损失测试一下
    #    print(iteration, loss)
        
    gradients = 1/m * X_train.T.dot(error)    #偏导数
    Theta = Theta - eta * gradients

验证精度

logits = X_valid.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)

accuracy_score = np.mean(y_predict == y_valid)

• 可以看到精度达到了 $0.966666...$，至于如何再提高精度，可以逐步减小步长，或者提前停止，这里就暂不讨论了