Softmax Regression --【python代码实现】

上次的Logistic Regression 只能处理二元分类问题，而在其推广的Softmax回归能处理多元分类问题->称为多元Logistic Regression。

---

有 K 个 类别就有 K 个 w_j列矩阵。

• ​ 因为对于每个实例x，都要计算ta成为K个特征的概率

  • $$Z_j=g_j(x)=w_j^{T}x$$

  • X = \begin{pmatrix}
        1 & x_1 & x_1^2 & \cdots & x_1^n \\
        1 & x_2 & x_2^2 & \cdots & x_2^n \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        1 & x_m & x_m^2 & \cdots & x_m^n \\
    \end{pmatrix}
    \\
    

    $$\begin{gathered} W=\left(\begin{array}{c}{w}_1^T\\ w_2^T\\ ⋮\\ w_K^T\end{array}\right) \\ z=g(x)=Wx=\left(\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_K\end{array}\right) \end{gathered}$$

softmax 函数预测概率

• softmax函数输出为一个向量（每一个分量 j 是x预测为 j 类别的概率）
  $$\sigma (z)=\left(\begin{array}{c}\sigma (z_1)\\ \sigma (z_2)\\ ⋮\\ \sigma (z_K)\end{array}\right)$$

  $$\sigma (z_j)=\frac{e^{z_j}}{{\sum }_{k=1}^K{e}^{z_k}}$$

  当k=2，分子分母同时除以e^zj就是logistic函数

Softmax回归模型（分子是列向量，展开）（概率）

• 回归模型假设函数

$$h_w(x)=\sigma (g(x))=\frac{1}{{\sum }_{k=1}^K{e}^{w_k^{T}x}}\left(\begin{array}{c}{e}^{w_1^{T}x}\\ e^{w_2^{T}x}\\ ⋮\\ e^{w_k^{T}x}\end{array}\right)$$

• 如果通过训练确定了参数模型W，可构建多选分类函数如下
  $$H(h_w(x))=argma{x}_k{h}_w(x{)}_k=argmax(w_k^{T}x)->K$$

损失函数，交叉熵:

$$J(W)=-\frac{1}{m}\sum _{i=1}^m\sum _{j=1}^{K}I(y_i=j)\ln h_w(x_i{)}_j$$

注: I 是 指示函数 y_j== j 的时候为1

梯度下降

• If w.shape in Logistic Regression is (n+1, 1), here it must be a matrix (n+1, K).

  要更新W就要更新每一个w_j,需要计算J(W)对每个w_j的梯度

  公式：
  $${\nabla }_{w_j}J(W)=\frac{1}{m}\sum _{i=1}^m(h_w(x_i{)}_j-I(y_i=j))x_i$$

  更新参数W的公式

  $$W：=W-\eta \left(\begin{array}{c}{\nabla }_{w_1}J(W{)}^T\\ {\nabla }_{w_2}J(W{)}^T\\ ⋮\\ {\nabla }_{w_K}J(W{)}^T\end{array}\right)$$

  ---

  感觉脑子不够用了…

代码实现（代码来自书籍）

import numpy as np

class SoftmaxRegression:
    def __init__(self, n_iter=200, eta=1e-3, tol=None):
        # 训练迭代次数
        self.n_iter = n_iter
        # 学习率
        self.eta = eta
        # 误差变化阈值
        self.tol = tol
        # 模型参数W(训练时初始化)
        self.W = None

    def _z(self, X, W):
        '''g(x)函数: 计算x与w的内积.'''
        if X.ndim == 1:
            return np.dot(W, X)   # 如果维度是一维的，就同书中公式那样进行乘积
        return np.matmul(X, W.T)  # 如果不是,则 (m, n+1) (n+1, K)  -> (m, K)

    def _softmax(self, Z):  # (m, K)
        '''softmax函数'''
        E = np.exp(Z)  # 上边公式中的分子
        if Z.ndim == 1:  # 维度为1，则直接类似于归一化的操作
            return E / np.sum(E)  # 归一化
        return E / np.sum(E, axis=1, keepdims=True)  # 保持维度的，同行相加 ->对应所有种类的值求sum

    def _predict_proba(self, X, W):
        '''h(x)函数: 预测y为各个类别的概率.'''
        Z = self._z(X, W)
        return self._softmax(Z)

    def _loss(self, y, y_proba):
        '''计算损失'''
        m = y.size  # (m, k)
        p = y_proba[range(m), y]  # 选去所有行的和？？？无法理解。可能是进行指示函数的步骤，但是难看。
        print(">>y:", y)
        print(">>p:", p)
        return -np.sum(np.log(p)) / m  #

    def _gradient(self, xi, yi, yi_proba):
        '''计算梯度'''
        K = yi_proba.size  #
        y_bin = np.zeros(K)
        y_bin[yi] = 1

        return (yi_proba - y_bin)[:, None] * xi  # 不行了, 这个操作好像了升维的。

    def _stochastic_gradient_descent(self, W, X, y):
        '''随机梯度下降算法'''

        # 若用户指定tol, 则启用早期停止法.
        if self.tol is not None:
            loss_old = np.inf
            end_count = 0

        # 使用随机梯度下降至多迭代n_iter次, 更新w.
        m = y.size
        idx = np.arange(m)
        for step_i in range(self.n_iter):
            # 计算损失
            y_proba = self._predict_proba(X, W)
            loss = self._loss(y, y_proba)
            print('%4i Loss: %s' % (step_i, loss))

            # 早期停止法
            if self.tol is not None:
                # 随机梯度下降的loss曲线不像批量梯度下降那么平滑(上下起伏),
                # 因此连续多次(而非一次)下降不足阈值, 才终止迭代.
                if loss_old - loss < self.tol:
                    print('haha')
                    end_count += 1
                    if end_count == 5:
                        break
                else:
                    end_count = 0

                loss_old = loss

            # 每一轮迭代之前, 随机打乱训练集.
            np.random.shuffle(idx)
            for i in idx:
                # 预测xi为各类别概率
                yi_proba = self._predict_proba(X[i], W)
                # 计算梯度
                grad = self._gradient(X[i], y[i], yi_proba)
                # 更新参数w
                W -= self.eta * grad


    def _preprocess_data_X(self, X):
        '''数据预处理'''

        # 扩展X, 添加x0列并置1.
        m, n = X.shape  # (m, n+1)
        X_ = np.empty((m, n + 1))
        X_[:, 0] = 1
        X_[:, 1:] = X

        return X_

    def train(self, X_train, y_train):
        '''训练'''

        # 预处理X_train(添加x0=1)
        X_train = self._preprocess_data_X(X_train)

        # 初始化参数向量W
        k = np.unique(y_train).size  # 获得种类数量
        _, n = X_train.shape  # (m, n+1)
        self.W = np.random.random((k, n)) * 0.05  # 初始化参数W (k种类, n+1特征)

        # 执行随机梯度下降训练W
        self._stochastic_gradient_descent(self.W, X_train, y_train)

    def predict(self, X):
        '''预测'''

        # 预处理X_test(添加x0=1)
        X = self._preprocess_data_X(X)

        # 对每个实例计算向量z.
        Z = self._z(X, self.W)

        # 向量z中最大分量的索引即为预测的类别.
        return np.argmax(Z, axis=1)