多分类问题

一、引言

本次练习中，我们将会实现一个一对多的逻辑回归以及用于识别手写数字的神经网络。当然和之前练习一样，我们需要从GitHub上下载本次练习所需的数据："ex3data1.mat","ex3data2.mat"

二、多分类问题

在本小结中，我们会使用逻辑回归和神经网络来识别手写数字(0 -9)。如今，自动识别手写数字应用广泛——从识别邮件信封上的美国邮码(邮政编码)到识别银行支票上的数字。本次练习将会为你展示我们所学到的方法是怎样运用到分类任务中去的。

首先，在练习的第一部分，我们会沿用练习二已实现的逻辑回归并将其运用到一对多分类问题中去。

1. 导入数据集

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

data = loadmat('ex3data1.mat')
data

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
 '__version__': '1.0',
 '__globals__': [],
 'X': array([[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]),
 'y': array([[10],
        [10],
        [10],
        ...,
        [ 9],
        [ 9],
        [ 9]], dtype=uint8)}

data['X'].shape, data['y'].shape

((5000, 400), (5000, 1))

可以看出，data中的X的每一行是对应手写数字图像像素的灰度(20 * 20被拉直了)，y则是对应的数字标签

好的，我们已经加载了我们的数据。图像在X矩阵中表示为400维向量（其中有5,000个图像）。 400维“特征值”是原始20 x 20图像中每个像素的灰度强度。类标签在向量y中作为表示图像中数字的数字类。

第一个任务是将我们的逻辑回归实现修改为完全向量化（即没有“for”循环）。这是因为向量化代码除了简洁外，还能够利用线性代数优化，并且通常比迭代代码快得多。

2. Sigmoid函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： $Sigmoid(z)= \frac{1}{1+e^{-z}}$

合起来，我们得到逻辑回归模型的假设函数： $z = \theta^Tx$

$h_\theta(x) = \frac{1}{1+e^{-\theta^Tx}}$

def sigmoid(z):
    return 1/(1+np.exp(-z))

3. 向量化代价函数

我们之前已经推断出正则化逻辑回归的代价函数如下： $J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]+\frac{\lambda }{2m}\sum^n_{j = 1}\theta_j^2$

定义出对应代价函数：

# 注意此处的参数顺序，约束了后续的fminix
def RegCost(theta, X , y ,  learningrate):
    # 先进性数据转换便于之后矩阵乘法
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    
    # 此处的向量和标量相乘不能直接用*，而是要用np.multiply
    
    first =np.multiply( -y , np.log(sigmoid(X * theta.T)))
    second = np.multiply((1-y) , np.log(1- sigmoid(X * theta.T)))
    # 计算正则化项
    reg = learningrate/(2 * len(X)) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    # 向量转标量
    return np.sum(first - second)/len(X)+reg

4. 向量化梯度函数

如果我们要使用梯度下降法令这个代价函数最小化，因为我们未对 $\theta_0$ 进行正则化，所以梯度下降算法将分两种情形：

$\begin{array}{l} Repeat until convergence\{\\ \theta_{0}:=\theta_{0}-a \frac{1}{m} \sum_{i=1}^{m}\left[h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right] x_{0}^{(i)} \\ \theta_{j}:=\theta_{j}-a \frac{1}{m} \sum_{i=1}^{m}\left[h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right] x_{j}^{(i)}+\frac{\lambda}{m} \theta_{j} \\\} Repeat \end{array}$

以下是使用for循环的梯度函数的原始代码(后面我们会将其替换为向量相乘)：

def gradient_with_loop(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    # 得到theta矩阵的元素个数，也就是需要求得的梯度值个数
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    # 计算误差
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            # 对添加的全一列不用加上正则化项的偏导
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad

向量化的梯度函数

A. 具体代码

def gradient(theta, X, y, learningRate):
    # 修改数据类型
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    # 得到偏导个数
    parameters = int(theta.ravel().shape[1])
    # 计算预测误差
    error = sigmoid(X * theta.T) - y
    # 利用向量计算梯度，后面一项是正则化项的偏导
    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
    
    # 额外插入的全一特征是不用正则化的
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()

B. 推理过程

我们的原始梯度下降函数利用循环来对梯度矩阵中的所有偏导依次赋值，对第i个样本的第j个特征( $x^{(i)}_j -\theta_j$ )的偏导数如下：

$\frac{\partial J}{\partial \theta_j} = \frac{1}{m} \sum^m_{i = 1}((h_\theta(x^{(i)}) - y^{(i)} )x_j^{(i)})$

为了将数据集中的求导操作向量化，我们先显式写出各偏导计算的的具体内容（本练习中的样本特征个数n = 400）：

$\begin{aligned} {\left[\begin{array}{c} \frac{\partial J}{\partial \theta_{0}} \\ \frac{\partial J}{\partial \theta_{1}} \\ \frac{\partial J}{\partial \theta_{2}} \\ \vdots \\ \frac{\partial J}{\partial \theta_{n}} \end{array}\right]=\frac{1}{m}\left[\begin{array}{c} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{0}^{(i)}\right) \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{1}^{(i)}\right) \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{2}^{(i)}\right) \\ \vdots \\ \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{n}^{(i)}\right) \end{array}\right] } \\ =\frac{1}{m} \sum_{i=1}^{m}\left(\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x^{(i)}\right)\\=\frac{1}{m} X^{T}\left(h_{\theta}(x)-y\right) . \end{aligned}$

此处的 $h_\theta(x) -y$ 是一个n维向量： $h_{\theta}(x)-y=\left[\begin{array}{c} h_{\theta}\left(x^{(1)}\right)-y^{(1)} \\ h_{\theta}\left(x^{(2)}\right)-y^{(2)} \\ \vdots \\ h_{\theta}\left(x^{(1)}\right)-y^{(m)} \end{array}\right]$

注意： $x^{(i)}$ 是一个特征向量，但是 $h_\theta(x^{(i)}) - y^{(i)}$ 是一个标量，为了便于理解求偏导的最后一步，我们令 $\beta_i = (h_\theta(x^{(i)}) - y^{(i)})$ ，可以看出：

$\sum_{i} \beta_{i} x^{(i)}=\left[\begin{array}{cccc} \mid & \mid & & \mid \\ x^{(1)} & x^{(2)} & \ldots & x^{(m)} \\ \mid & \mid & & \mid \end{array}\right]\left[\begin{array}{c} \beta_{1} \\ \beta_{2} \\ \vdots \\ \beta_{m} \end{array}\right]=X^{T} \beta$

上面的表达式可以让我们省去任何循环，加快求偏导的速度(但是需要对插入一列的特征进行特殊处理)。

现在我们已经定义了向量化的代价函数和梯度函数，现在是构建分类器的时候了。对于这个任务，我们有10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。在本练习中，我们的任务是实现一对一全分类方法，其中具有k个不同类的标签就有k个分类器，每个分类器在“类别 i”和“不是 i”之间决定。我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回为k *（n + 1）数组，其中n是参数数量(特征个数)。

5. 构建分类器(10个1对2)

from scipy.optimize import minimize

def one_vs_all(X, y, num_labels, learning_rate):
    # 得到X的样本个数和特征个数
    rows = X.shape[0]
    params = X.shape[1]
    
    # 定义一个用于存储每个分类器中的最终权重的 k*(n + 1)二维数组
    # 也就是存储10个二分类器最终参数的矩阵，添加的1列作偏移
    all_theta = np.zeros((num_labels, params + 1))
    
    # 在数据集的第一列上插入全一列，便于矩阵相乘
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # 注意，我们之前展示过标签内容，是从 1 开始而非 0开始的
    for i in range(1, num_labels + 1):
        # 每个二分类器的参数是一行401个
        theta = np.zeros(params + 1)
        # 进行“是或不是”的二分类问题
        y_i = np.array([1 if label == i else 0 for label in y])
        # 将标签向量转换为仅有标签为 i 的元素为 1 的列向量
        y_i = np.reshape(y_i, (rows, 1))
        
        # 最小化目标函数并将特征最优参数的行向量存储起来
        fmin = minimize(fun=RegCost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        # 注意预测0的二分类器的参数
        all_theta[i-1,:] = fmin.x
        
    # 将十个分类器最优参数矩阵悉数返回
    return all_theta

这里需要注意的几点：首先，我们为theta添加了一个额外的参数（与训练数据一列），以计算截距项（常数项）。其次，我们将y从类标签转换为每个分类器的二进制值（要么是类i，要么不是类i）。最后，我们使用SciPy的较新优化API来最小化每个分类器的代价函数。如果指定的话，API将采用目标函数，初始参数集，优化方法和jacobian（梯度）函数。然后将优化程序找到的参数分配给参数数组。

实现向量化代码的一个更具挑战性的部分是正确地写入所有的矩阵，保证维度正确。

# 原始数据集的尺寸
rows = data['X'].shape[0]
params = data['X'].shape[1]

# 得到的十分类矩阵(10个二分类构成)的最优参数矩阵
all_theta = np.zeros((10, params + 1))

# 原始数据集插入一列(值全一)
X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)

# 定义X特征个数（插入后），也即每个二分类器的参数个数
theta = np.zeros(params + 1)

# 分离标签为0的列向量
y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))

X.shape, y_0.shape, theta.shape, all_theta.shape

((5000, 401), (5000, 1), (401,), (10, 401))

注意，theta是一维数组，因此当它被转换为计算梯度的代码中的矩阵时，它变为（1×401）矩阵。我们还检查y中的类标签，以确保它们看起来像我们想象的一致。

# 看下有几类标签
np.unique(data['y'])

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10], dtype=uint8)

让我们确保我们的训练函数正确运行，并且得到合理的输出。

all_theta = one_vs_all(data['X'], data['y'], 10, 1)
all_theta

array([[-2.38262219e+00,  0.00000000e+00,  0.00000000e+00, ...,         1.30443025e-03, -7.48519061e-10,  0.00000000e+00],
       [-3.18328834e+00,  0.00000000e+00,  0.00000000e+00, ...,         4.45565933e-03, -5.07997883e-04,  0.00000000e+00],
       [-4.79850113e+00,  0.00000000e+00,  0.00000000e+00, ...,        -2.87680736e-05, -2.47832460e-07,  0.00000000e+00],
       ...,
       [-7.98746489e+00,  0.00000000e+00,  0.00000000e+00, ...,        -8.94565205e-05,  7.20837313e-06,  0.00000000e+00],
       [-4.57328172e+00,  0.00000000e+00,  0.00000000e+00, ...,        -1.33528996e-03,  9.98190574e-05,  0.00000000e+00],
       [-5.40500337e+00,  0.00000000e+00,  0.00000000e+00, ...,        -1.16588361e-04,  7.87229985e-06,  0.00000000e+00]])

我们现在准备好最后一步 - 使用训练完毕的分类器预测每个图像的标签。对于这一步，我们将计算每个类的类概率，对于每个训练样本（向量化代码），并将输出类标签为具有最高概率的类(极大似然)。

6. 使用分类器进行预测

# 定义预测函数，传入向量X以及十个分类器得到的最佳参数
def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # 和前面相同，插入一列便于矩阵乘法
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # 转换为矩阵
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # 为每个训练样例预测每类的概率，X-(5000,401) ,all_theta.T - (401, 10),那么 h - (5000, 10)
    h = sigmoid(X * all_theta.T)
    
    # 得到概率值最大标签的下标值(从0开始)
    h_argmax = np.argmax(h, axis=1)
    
    # 加上1偏移量作为最终结果
    h_argmax = h_argmax + 1
    
    # 返回预测结果向量
    return h_argmax

现在我们可以使用predict_all函数为每个实例生成类预测，看看我们的分类器是如何工作的。

y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))

accuracy = 94.46%

练习三 - 多分类问题