# 题目

利用神经网络识别手写数字。

# 任务

# 多级回归

在这个练习中，你将使用逻辑回归和神经网络去识别 0-9 的手写数字。手写数字自动识别目前应用已经十分广泛，从邮政编码到银行支票都有它的身影。这个练习将向你表示你所学的方法是如何被应用在多级回归任务中的。

练习的第一部分，你需要扩展在逻辑回归中的代码并将它应用到一对多的分类任务中。

# 数据集

数据集存放在 ex3data1.mat 中，它包含了 5000 个手写数字的训练样本。 “.mat” 文件代表着数据以矩阵形式被存放，而非 csv 格式的文件那样以文本形式存放。通过 load 命令，就可以将他们在你的程序中调用。矩阵已被命名，所以你无需再为他们指定名字。

在 ex3data1.mat 中有 5000 个训练样本，其中每个训练样本都是一个 20*20 像素的灰度图像。每个像素由一个浮点数表示，表示该位置的灰度强度。将每一个 20*20 的像素网格 “展开” 成一个 400 维的向量，组成数据矩阵 X 中的一行。我们将得到一个 5000*400 的矩阵 X，其中每一行都是手写数字图像的训练示例。

$$X =\begin{bmatrix}(-(x^{(1)})^T-) \\ (-(x^{(2)})^T-) \\ \cdots \\ (-(x^{(k)})^T-) \end{bmatrix}$$

训练集的第二部分是一个 5000 维的包含训练集标签的 向量y 。为了使数据和程序的索引更匹配，其中不包含索引 “0” 。我们将数字 0 标记为 10，而 1-9 则按顺序标记。

# 可视化

你首先要可视化训练集的子集。在第 1 部分中，代码从 X 中随机选择 100 行，并将这些行传递给一个可以显示图像的函数。 该函数将每一行映射为 20*20 像素空间的灰度图像显示在屏幕上。在你操作之后，你将得到类似下面的图像。

# 向量化逻辑回归

你将使用多个一对多逻辑回归模型来构建多级分类器。由于现在有 10 个类，所以我们需要训练 10 个独立的逻辑回归分类器。为了使这种训练有效，重要的是确保您的代码很好地向量化。在这一部分，你将实现一个不使用任何 for 循环的向量化逻辑循环。

# 向量化代价函数

首先我们要写一个向量化版本的代价函数。回想一下，非正则化的逻辑回归函数是这样的：

$$J(\theta)=-\frac{1} {m}\left[\sum_{i=1}^{m}{y}^{(i)}\log{h_\theta({x}^{(i)})}+\left(1-{y}^{(i)}\right)log\left(1-h_\theta\left({x}^{(i)}\right)\right)\right]$$

我们需要计算每个样例的假设函数，S 形函数。如果我们将其转换为向量形式，就会方便很多：

$$X =\begin{bmatrix}(-(x^{(1)})^T-) \\ (-(x^{(2)})^T-) \\ \cdots \\ (-(x^{(k)})^T-) \end{bmatrix} \, and \: \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \cdots \\ \theta_n \end{bmatrix}$$

然后我们计算矩阵$\theta \cdot X$, 我们将得到：

$$\theta \cdot X =\begin{bmatrix}(-\theta^T(x^{(1)})-) \\ (-\theta^T(x^{(2)})-) \\ \cdots \\ (-\theta^T(x^{(k)})-) \end{bmatrix}$$

你需要像上面这样编写完成你的代价函数。

# 向量化梯度函数

未正则化的梯度函数是这样定义的：

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

对其进行向量化操作，我们要写出关于$\theta$ 的所有偏导数项，

$$\begin{aligned} \begin{bmatrix}\frac{\partial J}{\partial { {\theta }_0 } } \\ \frac{\partial J}{\partial { {\theta }_0 } } \\ \cdots \\ \frac{\partial J}{\partial { {\theta }_0 } } \end{bmatrix} &= \begin{bmatrix}\sum\limits_{i=1}^{m}{ {\left( {h_\theta}\left( {x}^{\left( i \right)} \right)-{y}^{\left( i \right)} \right)} } {x}_{0}^{(i)} \\ \sum\limits_{i=1}^{m}{ {\left( {h_\theta}\left( {x}^{\left( i \right)} \right)-{y}^{\left( i \right)} \right)} } {x}_{2}^{(i)} \\ \cdots \\ \sum\limits_{i=1}^{m}{ {\left( {h_\theta}\left( {x}^{\left( i \right)} \right)-{y}^{\left( i \right)} \right)} } {x}_{n}^{(i)} \end{bmatrix} \\&=\frac{1}{m} \sum\limits_{i=1}^{m}{ {\left( {h_\theta}\left( {x}^{\left( i \right)} \right)-{y}^{\left( i \right)} \right)} } {x}_{0}^{(i)} \\&=\frac{1}{m}X^T(h_{\theta}(x)-y) \end{aligned}$$

$$h_{\theta}(x)-y=\begin{bmatrix}{h_\theta}\left( {x}^{\left( 1 \right)} \right)-{y}^{\left( 1 \right)} \\ {h_\theta}\left( {x}^{\left( 2 \right)} \right)-{y}^{\left( 2 \right)} \\ \cdots \\ {h_\theta}\left( {x}^{\left( m \right)} \right)-{y}^{\left( m \right)} \end{bmatrix}$$

注意$x^{(i)}$ 是向量，而$h_{\theta}(x)-y$ 是标量，为了理解这一步的由来，我们令$\beta_i =h_{\theta}(x)-y$ ，则

$$\sum_{i}\beta_ix^{(i)}=\begin{bmatrix} x^{(1)} & x^{(2)} & \cdot \cdot \cdot & x^{(m)}\end{bmatrix} \begin{bmatrix}\beta_1 \\ \beta_2 \\ \cdot \cdot \cdot \\ \beta_m \end{bmatrix}=X^T\beta$$

上面的表达式能让我们在不适用循环的方式下计算所有的偏导数。

# 向量化正则逻辑回归

当你将逻辑回归实现向量化后，你需要将正则项加入到代价函数。我们正则化后的代价函数长这样：

$$J(\theta)=-\frac{1}{m}[\sum_{i=1}^{m}{y}^{(i)}\log{h_\theta({x}^{(i)})}+(1-{y}^{(i)})\log(1-h_\theta({x}^{(i)}))]+\frac{\lambda}{2m}\sum_{j=1}^{n}\theta_j^2$$

由于我们的$\theta_0$ 项没有偏差项目，所以代价函数的偏导数应该这样定义:

$$\frac{\partial }{\partial { {\theta }_{j} } }J(\theta)=\frac{1}{m}\sum\limits_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{i})x_j^{(i)})\hspace{1em}for\;j=0 \\ \frac{\partial }{\partial { {\theta }_{j} } }J(\theta)=\frac{1}{m}\sum\limits_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{i})x_j^{(i)})+\frac{\lambda}{m}\theta_j\hspace{1em}for\;j\geq1$$

# 一对多分类

在这一部分中，你将通过训练多个正则化逻辑回归分类器来实现一对多分类。

你的代码应该返回矩阵$\Theta\in R^{K\times(N1)}$ 中的所有分类器参数，其中每行的$\Theta$ 对应于一个类的学习逻辑回归参数。 你可以使用从 1 到 K 的 “for” 循环来完成此操作，独立地训练每个分类器。 请注意，此函数的 y 参数是标签为 1 到 10 的向量，其中我们将数字 0 映射到标签 10（以避免与索引混淆）。

当训练类$k\in{\{1,\dots,K\} }$ 的分类器时，你需要已经标号的 m 维的向量标签 y，它将高速你第 j 个训练实例是否属于该类，属于则$y_j=1$, 不属于则$y_j=0$。

# 一对多预测

在一对多分类器训练好之后，就可以使用它来预测给定图像中的数字。对于每个输入，你应该使用逻辑回归分类器计算它属于每个类的概率，一对多预测函数将选择逻辑回归分类器输出最高概率并返回类标签 (1,2，…，K) 作为输入示例的预测。输出的准群率应该在 94.9% 左右。

# 神经网络

在本部分练习中，您将实现一个神经网络，使用与以前相同的训练集来识别手写数字。 神经网络将能够形成非线性假设的复杂模型。

# 模型呈现

我们的神经网络如图 2 所示。 它具有三层 —— 输入层、隐藏层和输出层。 回想一下，我们的输入是数字图像的像素值。 由于图像大小为 20×20，这给了我们 400 个输入层单元（不包括总是输出 1 的额外偏置单元）。 和以前一样，训练数据将被加载到变量 X 和 y 中。 你已经获得了一组参数，这些参数被存储在 ex3weights.mat 中，这些参数将用于第二层中的 25 个单元和 10 个输出单元。

# 前馈传播与预测

现在，你将实现神经网络的前馈传播。 您需要完成代码来返回神经网络的预测。 你应该实现前馈计算，计算每个示例 i 的$(h_\theta(x))_k$，并返回相关的预测。 类似于一对多分类策略，预测将输出概率最大的标签。最后输出的准确率应该在 97.5% 左右。

# 代码

# 多级回归

#exercise-3 neural network

from scipy.io import loadmat

import scipy.optimize as opt

import matplotlib.pyplot as plt

import matplotlib

import numpy as np

from sklearn.metrics import classification_report #用于模型评价

# 加载数据

def load_data(path,transpose=True):

    data = loadmat(path)

    y = data.get('y')   #(5000,1)

    y = y.reshape(y.shape[0])

    x = data.get('X')   #(5000,400)

    if transpose: #转置矩阵，不然图像会横着

        x = np.array([im.reshape((20,20)).T for im in x])

        x = np.array([im.reshape(400) for im in x])

    return x, y

X,y = load_data("ex3data1.mat") # shape of X is (5000,400), y is (5000,1)

raw_X, raw_y = load_data('ex3data1.mat')

# 数据可视化

def plot_100_images(X):

    'sample 100 images and show them in a square'

    size = int(np.sqrt(X.shape[1]))

    sample_idx = np.random.choice(np.arange(X.shape[0]),100)

    sample_images = X[sample_idx, :]

    fig, ax = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(8,8))

    for r in range(10):

        for c in range(10):

            ax[r,c].matshow(sample_images[10*r+c].reshape(size,size), cmap=matplotlib.cm.binary)

            plt.xticks(np.array([]))

            plt.yticks(np.array([]))

plot_100_images(X)

plt.show()

# 处理数据

X = np.insert(raw_X, 0, values=np.ones(raw_X.shape[0]), axis=1)

y_matrix = []

for k in range(1,11):

    y_matrix.append((y == k).astype(int))

y_matrix = [y_matrix[-1]] + y_matrix[:-1]

y = np.array(y_matrix)

由于最后一列是用 10 代替的 0，这里将最后一列移到第一列。for 循环被用于向量化标签，处理之后 y 的 shape 为（10,5000）

# 训练一维模型

def sigmoid(z):

    return 1 / (1 + np.exp(-z))

def cost(theta, X, y):

    return np.mean((-y) *np.log(sigmoid(X @ theta))- (1-y)*np.log(1-sigmoid(X@theta)))

def regularized_cost(theta, X, y, l=1):

    theta_j1_to_n = theta[1:]

    regularized_term = l/(2*len(X))*np.power(theta_j1_to_n, 2).sum()

    return cost(theta, X, y) + regularized_term

def gradient(theta, X, y):

    return 1/len(X) * X.T @ (sigmoid(X @ theta) - y)

def regularized_gradient(theta, X, y, l=1):

    theta_j1_to_n = theta[1:]

    regularized_theta = 1/ len(X)* theta_j1_to_n

    regularized_term = np.concatenate([np.array([0]),regularized_theta])

    return gradient(theta, X, y) + regularized_term

def logistic_regression(X, y, l=1):

    theta = np.zeros(X.shape[1])

    res = opt.minimize(fun=regularized_cost, 

                       x0=theta, 

                       args=(X, y, l), 

                       method='TNC', 

                       jac=regularized_gradient, 

                       #options={'disp':True} #打印信息

                       )

    final_theta = res.x

    return final_theta

def predict(x, theta):

    prob = sigmoid(x @ theta)

    return ( prob >= 0.5).astype(int)

t0 = logistic_regression(X, y[0])

y_pred = predict(X, t0)

print('Accuracy={}'.format(np.mean(y[0]==y_pred)))

这里是我们需要的各种函数，在之前的练习中都做过。$t_0$ 的 shape 为（401,1），拿 y [0] 进行训练，得到的准确率为 0.9974。

# 训练 K 维模型

k_theta = np.array([logistic_regression(X, y[k]) for k in range(10)])	

prob_matrix = sigmoid(X @ k_theta.T) 

#k_theta.shape->(10,401)

#X.shape->(5000,401)

#prob_matrix.shape->(5000,10)

np.set_printoptions(suppress = True) 	#输出方式控制，小数无需以科学计数法输出

y_pred = np.argmax(prob_matrix,axis=1)  #返回最大值的索引

y_answer = raw_y.copy()

y_answer[y_answer==10] = 0

这部分代码是将数据喂入，得到每组数据属于各组的概率并取概率最高者。

# 评估结果

print(classification_report(y_answer, y_pred))

# 神经网络

def load_weight(path):

    data = loadmat(path)

    return data['Theta1'], data['Theta2']

theta1, theta2 = load_weight('ex3weights.mat')

#theta1.shape->(25,401)

#theta2.shape->(10,26)

X,y = load_data('ex3data1.mat',transpose=False)

X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)  # intercept

要注意的是，这里所给的权重并没有作转置，所以我们重新调用未经转置的原始数据。

# 前馈预测

a1 = X					#a1.shape->(5000,401)

z2 = a1 @ theta1.T		#z2.shape->(5000,25)

z2 = np.insert(z2, 0,	 values=np.ones(z2.shape[0]), axis=1)

a2 = sigmoid(z2)		#a2.shape->(5000,26)

z3 = a2 @ theta2.T		#z3.shape->(5000,10)

a3 = sigmoid(z3)		#a3.shape->(5000,10)

y_pred = np.argmax(a3, axis=1) +1

# 评估结果

print(classification_report(y, y_pred))