Machine Learning（吴恩达＜二＞）神经网络

-小透明-

已于 2022-06-28 09:45:18 修改

阅读量4.1k

点赞数 2

分类专栏：吴恩达机器学习笔记文章标签：机器学习神经网络人工智能深度学习

于 2021-11-12 18:12:50 首次发布

本文链接：https://blog.csdn.net/qq_54809548/article/details/121196101

版权

吴恩达机器学习笔记专栏收录该内容

6 篇文章 5 订阅

订阅专栏

本文详细探讨了引入非线性假设在计算机视觉中的必要性，通过神经元模型理解神经网络工作原理，展示了从模型构建到反向传播优化的完整流程，重点介绍了手写体识别中的神经网络应用，以及如何利用神经网络进行多元分类。

摘要由CSDN通过智能技术生成

如何高效进行计算并展现一个向量化的实现方法？

非线性假设

这个时候用线性回归是不能很好的去拟合，并且就算能拟合，表达式也会很复杂。

计算机视觉中一个问题：假设你想用机器学习算法训练一个分类器，来检测图像，来判断图像是否是一辆汽车。对人来说，一看就知道；但对机器来说，是很艰难的一件事。

因为，当人看一辆车时，计算机看到的是一个数据矩阵（或者说一个表示像素强度值得网络）告诉我们图片中每个像素点的亮度值，计算机需要通过这些像素亮度值矩阵告诉我们这些数值代表一个什么（汽车把手）

当我们要训练一个汽车分类器时，我们需要提供一个带标签的样本集，其中一些样本是各类汽车，另外一些不是汽车。

引入非线性假设的必要性

：逻辑回归当n（特征features）过大时，是不在适合的。

而神经网路被证明是（无论n多大，都能很轻松的解决）的一种方法。

神经元与大脑

背景

神经网络的起源是人们想设计出模仿大脑的算法，其理念是：如果想要建立学习系统，何不去模仿人的大脑呢。在19世纪左右就流行了一段时间，后来慢慢不提了，最近几年又火了。

可以假设大脑做的一些事情是不需要写一堆复杂的程序一个一个执行的，相反，只是通过一种学习算法就可以了(让大脑自己学习）。

模型展示1

如何用神经网络表示我们的假设和模型

首先看看单个的神经元（大脑中的细胞）它有一个细胞体，和许多"输入线"（即树突）来接收其他神经元的消息，他还有一个"输出通道"（即轴突）给其他神经元传递信号。

单个神经元是一个计算单元，他从输入通道接收一定数目的信息，并做出一些计算，通过输出通道传递信息。下图其逻辑单元图中x0=1成为“偏执单元或偏执神经元”；另外，神经网络中的一个术语“激活函数”指的就是逻辑回归中的g(z)sigmoid函数；另外，将θ称为“权重”或“参数”。

下面这幅图是多个神经元组成的一个网络，术语又来了：通常神经网路的第一层称为输入层（因为我们在这一层输入特征）；最后一层称为输出层（因为我们通过这一层神经元输出假设的最终结果）；中间第二层称为“隐藏层”，图中只有三层，实际上，隐藏层会有多个，所有非输入输出层都称作隐藏层。

神经网络具体的计算步骤

下图一些符号的意思：

$a_{i}^{(j)}$ :表示第j层的第i个神经(单)元的激活项(指由一个具体的神经元计算并输出的值)

$\theta ^{(j)}$ :权重矩阵（神经网络被这些矩阵参数化）他控制从某一层（如第一层或第二层）到（如第三层）的映射

如果一个神经网络在第s_ j层有j个单元，第s_(j+1)层有j+1个单元，则 $\theta ^{(j)}$ 则控制第j层到第j+1层的映射的矩阵，其维度为s_(j+1)*(s_ j + 1)

最后一个输出层的神经元计算 $h_{(\theta )}(x)$ 注意 $a^{(3)}_{1}$ 的表达式中 $\theta ^{(2)}$ 这个权重矩阵控制的是从第二层的三个单元到第三层的一个单元的映射

总结：上图通解释了如何定义一个人工神经网络，其中神经网路定义了函数h，从输入x到输出y的映射，这些假设被参数化，将参数记为大写的θ，这样，改变θ就可以得到不同的假设。我们有不同的函数比如从x到y的映射。 以上是如何从数学上定义和计算神经网络的假设函数。

模型展示2

如何高效进行计算并展现一个向量化的实现方法？

下图计算h(x)的过程也称为向前传播（从输入单元的激活项开始进行向前传播给隐藏层，计算隐藏层的激活项，然后继续进行向前传播计算输出层的激活项，这个过程叫做向前传播）

下图推导是这一过程的向量化实现方法，如果用图中右侧的公式来实现这一过程，会得到一个相对有效的计算h(x)的方法

下图，如果挡住第一层，只看蓝色的公式，它的hx就是逻辑回归，不同的是，它的特征不是原始的特征值x1，x2，x3，而是用新的隐藏层计算的a1，a2，a3作为新的特征。有趣的是，ai是从学习得到的输入值（就是从第一层到第二层的映射函数，该函数有 $\theta^{(1)}$ 决定）。

在神经网络中，它没有用输入的特征x1，x2，x3来训练逻辑回归，而是自己训练逻辑回归的输入a1，a2，a3，可以根据为 $\theta^{(1)}$ 选择的不同参数有时可以学习到一些复杂的有趣的特征，就可以得到一个比直接使用输入特征更好的假设函数。也可以选择多项式项作为输入，这个算法可以快速灵活的学习任意特征项。

神经网路的架构

神经网络中不同神经元的连接方式称为神经网络的架构。下图有三层隐藏层，通过三层能够训练出更加复杂的特征作为输入，以此得到一个复杂的非线性假设函数。

例子与直觉理解1

AND：逻辑与；当x1=1ANDx2=1时，z=-30+20+20=10；g(z)≈1

例子与直觉理解2

h(x)计算了非线性分类的(not x1) and (not x2)

手写体的识别过程

利用神经网络 输入的原始像素点（图片）通过隐藏层1得到一些特征，进而在通过一层隐藏层，计算出更复杂的特征...,最终，这些特将被用于逻辑回归分类器的最后一层，来准确预测出神经网络所看到的数字。如下图

多元分类

上面的手写体识别就是一个多类别识别问题，因为它需要识别 0~9 10个类别。

多类别识别具体实现细节

要在神经网络中实现多类别分类，采用的方法本质上是一对多法的扩展。

下图例子，我们需要识别出行人，汽车，摩托车，卡车四种类别。

以上是对如何构建假设模型来表示神经网络的细节补充，下面将说明如何构建训练集，以及如何让神经网络自动学习这些参数。

代价函数

一个栗子

下图中最上面（x，y）表示训练样本；

L 表示神经网络的层数

sl 表示每一层神经元的个数（当然不包括偏执单元1）

定义代价函数

使我们在逻辑回归中用的代价函数J(θ)的一种形式，不同的是，不再是一个逻辑回归输出的单元，而是K个。神经网络输出了属于R^K的向量；如果处理的是二元分类问题，则K可能=1。其中 $(h_{\theta}(x))_{i}$ 表示第i个输出。 $h_{\theta}(x)$ 是一个K维向量。

J(θ)中的第二项，就是在逻辑回归中所用到的正则化项。这里要除（θji中i=0，可以看出式子中都是从1开始求和的）对应于偏差值的项。这是因为，当我们计算神经元的激励值时，会有像θ_i0这样的项*x0+θ_i1*x1+...,含0的项类似于偏执项，在逻辑回归正则化时，我们是不想正则化这些项的，因此在这里，也不用加到正则项中，通常约定把他们置为0，当然，如果加入进去了，结果也是有效的，并且不会有很大的区别，只是不常见。

反向传播法优化代价函数

最小化代价函数

反向传播法

直观的讲就是对每一个结点计算出这样的一项 $\delta _{j}^{(l)}$ 表示了第l层的第j个结点的误差；前面 $a_{j}^{(l)}$ 代表着第l层第j个结点的激活值， $\delta$ 项在一定程度上捕捉到了我们在这个神经结点的激活值的误差。因此我们可能希望这个节点的激活值稍微不一样。下面的式子中，每一项也可是一个向量，并且每个向量的维度为输出单元的数目。下图中.*表示点乘。

反向传播的意思就是：我们从输出层开始计算 $\delta$ 项，然后返回上一层计算最后一个隐藏层的误差项，依次向后计算。因此类似于将输出值得误差反向传播给了一层层的隐藏层。

如何实现向量反传播算法来计算这些参数的偏导项，当我们有一个很大的（m个）训练样本时，如下图所示，下图中△是大写的 $\delta$ ，a^(1)输入层的激活函数=x^(i)第i个训练样本的输入值，然后用正向传播来计算第二层激活值，一直到最后一层；用反向传播计算δ(L), (L1)...(2)。强调：不用计算输入层的误差，因此没有δ(1)。用△ij来累计前面的偏导数项，如果把它写成向量形式，可以把δij看做矩阵ij为矩阵中的一个位置。写成红色的式子，这样实现了对所有i和j的自动更新值。

理解反向传播

是一个很复杂的“黑箱子”，有很复杂的算法，但很神奇，确实很难理解。

进一步理解向前传播的过程，代价函数cost(i)看做是一个"方差"。

再看看反向传播的过程，其中δ实际上是代价函数cost(i)关于z^(i)_j的偏导，也就是计算出z的项的加权和；cost(i)是关于标签y和神经网络中h(x)的输出值的函数

使用：将参数矩阵展开为向量

在函数中都假设θ为一个n+1维的向量，同样得到的梯度也是一个n+1维的向量

神经网络中θ不再是向量而是参数矩阵，同样有梯度下降矩阵

将矩阵展开为向量，以便最终成为恰当的格式能够传入到方程的参数

其中θ1矩阵中有110个数，thetaVec为三个矩阵组成的只有一列的向量（很多行）

reshape：即取出（例如110个元素）<重组函数>成一个(10*11)的矩阵

通过重组得到初始参数θ1 2 3，这样可以执行向前传播，反向传播，来计算导数和代价函数；最后我们可以取出这些导数值，然后展开它们，让他们保持和展开的θ值同样的顺序，D通过代价函数以向量的形式返回gradientVec。

梯度检测（gradient check）

神经网络算法使用反向传播算目标函数关于每个参数的梯度，可以看做解析梯度。由于计算过程中涉及到的参数很多，反向传播计算的梯度很容易出现误差，导致最后迭代得到效果很差的参数值。梯度检测从数值方法上，会计算梯度（或导数）的近似值，然后和程序计算出的梯度（导数）的值对比，以判断梯度（导数）值是否正确。

下图左面红色的公式是双侧差分，右面蓝色的为单侧差分，而双侧差分更精确。

总结：梯度检测过程

first，通过反向传播计算Devc

second，实现数值上的梯度检测

third ，保证上面两个值非常接近

last，最重要的一点，当你用反向传播训练之前，一定要关掉梯度检测！！<这是因为，检测算法计算量非常大，速度非常慢算导数的方法；而反向传播是一个高性能的算导数的方法。因此，一旦你通过梯度检测确定你的反向传播算法是正确的，就可以直接使用反向传播了。>

在优化前，先使用梯度检测先来确保你的算法是正确的。

随机初始化

如何对参数θ设初值

假设将所有参数初始化为0，参数同时改变并且最后参数同时相等，这时，每次输入隐藏单元都以相同的函数作为输入，这是一种冗余的现象，这阻止了神经网络学习更有趣的东西。未解决这一问题（对称权重(即权重都相等)），引入随机初始化。

将参数随机初始化为一个范围内接近0但不能为0的一个值

这里范围的 $\varepsilon$ 与梯度检测的 $\varepsilon$ 无关

总结：首先随机（打破对称性）初始化参数，然后使用反向传播计算梯度，做梯度检测，然后使用梯度下降或者更高级的优化算法最小化代价函数。

总结整体实现过程

选择

训练一个神经网络第一件事就是选择一种网络架构（架构：神经元之间的连接方式）

输入层的神经元个数取决于你即将输入的特征数量

输出层的神经元个数取决于你的分类问题最终分的类别数

下图中只有一个隐藏层是比较合理的默认结构；通常隐藏层的额单元越多越好（计算量也越大）

一般隐藏单元的个数还应与输入的维度相匹配（可以与特征数相同也可以是整数2、3、4倍）

训练

1. 随机初始化权重（一般初始化为很小的接近0的值）

2. 执行向前传播算法（即计算尺每一个输入算出对应的y）

3. 通过代码计算代价函数J(θ)

4.执行反向传播算法计算出偏导数项。（在第一次使用反向传播时，建议使用for循环实现，当然有更高级的算法，但目前还没有了解）

5. 使用梯度检测来比较计算的偏导数项，保证基本接近后不再执行检测代码

6. 使用一个最优化算法（如梯度下降等（目前只知道这一个））与反向传播算法结合来最小化参数θ，这里的代价函数是一个非凸函数，也就是可以达到局部最小位置

梯度下降在神经网络中的应用

神经网络中有很多θ，由于绘图限制，为更直观的感受，假设只有两个。

神经网络学习的例子

无人驾驶

涉及路径规划，图像识别（分类）等，神经网络算法的一个应用。

无人驾驶基本原理_人工智能AI技术-CSDN博客_无人驾驶技术原理

应用（copy的）

手写体识别

逻辑回归算法

# 数据集 ex3data1.mat
# 这是一个MATLAB格式的.mat文件，其中包含5000个20*20像素的手写字体图像，以及他对应的数字。另外，数字0的y值，对应的是10
# 用Python读取我们需要使用SciPy
# 数字0的y值，对应的是10
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from scipy.io import loadmat
from sklearn.metrics import classification_report#  评价报告

data = loadmat('E:/Python/machine learning/data/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)}

print(data['X'].shape, data['y'].shape)# 维度

# 随机展示100个数据
sample_idx = np.random.choice(np.arange(data['X'].shape[0]), 100)
sample_images = data['X'][sample_idx, :]
print(sample_images)

fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(12, 12))
for r in range(10):
    for c in range(10):
        ax_array[r, c].matshow(np.array(sample_images[10 * r + c].reshape((20, 20))).T,cmap=matplotlib.cm.binary)
        plt.xticks(np.array([]))
        plt.yticks(np.array([]))

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

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.]])

# 用多分类逻辑回归做一个分类器。因为现在有10个数字类别，所以需要训练10个不同的逻辑回归分类器
# 为了让训练效率更高，将逻辑回归向量化

# 向量化sigmod函数：Z为向量
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# 代价函数
def cost(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    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

# 向量化梯度
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)
    
    # intercept gradient is not regularized
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()

# 一对多分类器
from scipy.optimize import minimize

def one_vs_all(X, y, num_labels, learning_rate):
    rows = X.shape[0]
    params = X.shape[1]
    
    # k * (n + 1) array for the parameters of each of the k classifiers
    all_theta = np.zeros((num_labels, params + 1))
    
    # insert a column of ones at the beginning for the intercept term
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # labels are 1-indexed instead of 0-indexed
    for i in range(1, num_labels + 1):
        theta = np.zeros(params + 1)
        y_i = np.array([1 if label == i else 0 for label in y])
        y_i = np.reshape(y_i, (rows, 1))
        
        # minimize the objective function
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        all_theta[i-1,:] = fmin.x
    
    return all_theta


rows = data['X'].shape[0]
params = data['X'].shape[1]

all_theta = np.zeros((10, params + 1))

X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)

theta = np.zeros(params + 1)

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))

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.38208691e+00,  0.00000000e+00,  0.00000000e+00, ...,
         1.30384042e-03, -6.14427666e-10,  0.00000000e+00],
       [-3.18352842e+00,  0.00000000e+00,  0.00000000e+00, ...,
         4.46123909e-03, -5.08642939e-04,  0.00000000e+00],
       [-4.79735416e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -2.87309679e-05, -2.47481807e-07,  0.00000000e+00],
       ...,
       [-7.98467966e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -8.97640711e-05,  7.23641521e-06,  0.00000000e+00],
       [-4.57003568e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.33433208e-03,  1.00011405e-04,  0.00000000e+00],
       [-5.40502829e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.16485909e-04,  7.85363055e-06,  0.00000000e+00]])

# 一对多预测
def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # same as before, insert ones to match the shape
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # convert to matrices
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # compute the class probability for each class on each training instance
    h = sigmoid(X * all_theta.T)
    
    # create array of the index with the maximum probability
    h_argmax = np.argmax(h, axis=1)
    
    # because our array was zero-indexed we need to add one for the true label prediction
    h_argmax = h_argmax + 1
    
    return h_argmax

# 为每个实例生成类预测，看看分类器是如何工作的
y_pred = predict_all(data['X'], all_theta)
print(classification_report(data['y'], y_pred))

              precision    recall  f1-score   support

           1       0.95      0.99      0.97       500
           2       0.95      0.92      0.93       500
           3       0.95      0.91      0.93       500
           4       0.95      0.95      0.95       500
           5       0.92      0.92      0.92       500
           6       0.97      0.98      0.97       500
           7       0.95      0.95      0.95       500
           8       0.93      0.92      0.92       500
           9       0.92      0.92      0.92       500
          10       0.97      0.99      0.98       500

    accuracy       0.94      0.94      0.94      5000
   macro avg       0.94      0.94      0.94      5000
weighted avg       0.94      0.94      0.94      5000

逻辑回归并不能承载更复杂的假设，因为它是一个线性分类器。

神经网络

而神经网络可以表示一些非线性复杂的模型

神经网络预测

# 实现前馈神经网络预测手写数字

weight = loadmat("E:/Python/machine learning/data/ex3weights.mat")
theta1, theta2 = weight['Theta1'], weight['Theta2']
                 
print(theta1.shape, theta2.shape)
                 
# 插入偏执项
X2 = np.matrix(np.insert(data['X'], 0, values=np.ones(X.shape[0]), axis=1))
y2 = np.matrix(data['y'])
print(X2.shape, y2.shape)

a1 = X2
z2 = a1 * theta1.T
print(z2.shape)

a2 = sigmoid(z2)
print(a2.shape)

a2 = np.insert(a2, 0, values=np.ones(a2.shape[0]), axis=1)
z3 = a2 * theta2.T
print(z3.shape)

a3 = sigmoid(z3)
print(a3)

y_pred2 = np.argmax(a3, axis=1) + 1
print(y_pred2.shape)

print(classification_report(y2, y_pred))

(25, 401) (10, 26)
(5000, 401) (5000, 1)
(5000, 25)
(5000, 25)
(5000, 10)
[[1.12661530e-04 1.74127856e-03 2.52696959e-03 ... 4.01468105e-04
  6.48072305e-03 9.95734012e-01]
 [4.79026796e-04 2.41495958e-03 3.44755685e-03 ... 2.39107046e-03
  1.97025086e-03 9.95696931e-01]
 [8.85702310e-05 3.24266731e-03 2.55419797e-02 ... 6.22892325e-02
  5.49803551e-03 9.28008397e-01]
 ...
 [5.17641791e-02 3.81715020e-03 2.96297510e-02 ... 2.15667361e-03
  6.49826950e-01 2.42384687e-05]
 [8.30631310e-04 6.22003774e-04 3.14518512e-04 ... 1.19366192e-02
  9.71410499e-01 2.06173648e-04]
 [4.81465717e-05 4.58821829e-04 2.15146201e-05 ... 5.73434571e-03
  6.96288990e-01 8.18576980e-02]]
(5000, 1)
              precision    recall  f1-score   support

           1       0.95      0.99      0.97       500
           2       0.95      0.92      0.93       500
           3       0.95      0.91      0.93       500
           4       0.95      0.95      0.95       500
           5       0.92      0.92      0.92       500
           6       0.97      0.98      0.97       500
           7       0.95      0.95      0.95       500
           8       0.93      0.92      0.92       500
           9       0.92      0.92      0.92       500
          10       0.97      0.99      0.98       500
   micro avg       0.94      0.94      0.94      5000
   macro avg       0.94      0.94      0.94      5000
weighted avg       0.94      0.94      0.94      5000

使用反向传播的前馈神经网络，自动学习神经网络的参数。

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

X = data['X']
y = data['y']

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

# 前向传播函数
def forward_propagate(X, theta1, theta2):
    m = X.shape[0]
    
    a1 = np.insert(X, 0, values=np.ones(m), axis=1)
    z2 = a1 * theta1.T
    a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)
    z3 = a2 * theta2.T
    h = sigmoid(z3)
    
    return a1, z2, a2, z3, h

def cost(theta1, theta2 , input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # compute the cost
    J = 0
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    return J

# 原本的y是5000*1维的向量，需把它编码成5000*10的矩阵
encoder = OneHotEncoder(sparse=False)
y_onehot = encoder.fit_transform(y)
print(y_onehot.shape)

print(y[0], y_onehot[0,:]) # y0是数字0

(5000, 10)
[10] [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]

# 初始化设置
input_size = 400
hidden_size = 25
num_labels = 10
learning_rate = 1

# 用ex4weights.mat中给定的theta1 和 theta2 计算初始代价
cost(theta1, theta2, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

0.2876291651613188

# 正则化代价函数
def costReg(theta1, theta2, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)

    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # compute the cost
    J = 0
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    return J
# 计算初始代价
costReg(theta1, theta2, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

0.38376985909092354

# 作用在矩阵上时，应该是计算每个元素的梯度
def sigmoid_gradient(z):
    return np.multiply(sigmoid(z), (1 - sigmoid(z)))
# z=0的梯度
sigmoid_gradient(0)

0.25

# 随机初始，设置范围值为-0.12~0.12
# np.random.random(size) 返回size大小的0-1随机浮点数
params = (np.random.random(size=hidden_size * (input_size + 1) + num_labels * (hidden_size + 1)) - 0.5) * 0.24

# 反向传播
def backprop(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # initializations
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m

    # perform backpropagation
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        d3t = ht - yt  # (1, 10)
        
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        
        delta1 = delta1 + (d2t[:,1:]).T * a1t
        delta2 = delta2 + d3t.T * a2t
        
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    return J, delta1, delta2

# 梯度检测运行一次太慢了，因此跳过了

# 正则化神经网络
def backpropReg(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # initializations
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    # perform backpropagation
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        d3t = ht - yt  # (1, 10)
        
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        
        delta1 = delta1 + (d2t[:,1:]).T * a1t
        delta2 = delta2 + d3t.T * a2t
        
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    # add the gradient regularization term
    delta1[:,1:] = delta1[:,1:] + (theta1[:,1:] * learning_rate) / m
    delta2[:,1:] = delta2[:,1:] + (theta2[:,1:] * learning_rate) / m
    
    # unravel the gradient matrices into a single array
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
    
    return J, grad

# 使用工具库计算参数最优解
from scipy.optimize import minimize

# minimize the objective function
fmin = minimize(fun=backpropReg, x0=(params), args=(input_size, hidden_size, num_labels, X, y_onehot, learning_rate), 
                method='TNC', jac=True, options={'maxiter': 250})
fmin

     fun: 0.41235724537034724
     jac: array([-9.72188479e-04,  1.49551756e-06,  6.61280427e-06, ...,
        7.91347749e-05, -5.32445291e-04, -1.99904227e-04])
 message: 'Linear search failed'
    nfev: 219
     nit: 16
  status: 4
 success: False
       x: array([-0.29219349,  0.00747759,  0.03306402, ..., -0.50501111,
        0.56041664, -1.53064174])

X = np.matrix(X)
thetafinal1 = np.matrix(np.reshape(fmin.x[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
thetafinal2 = np.matrix(np.reshape(fmin.x[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))

# 计算使用优化后的θ得出的预测
a1, z2, a2, z3, h = forward_propagate(X, thetafinal1, thetafinal2 )
y_pred = np.array(np.argmax(h, axis=1) + 1)
y_pred

# 预测值与实际值比较
from sklearn.metrics import classification_report#这个包是评价报告
print(classification_report(y, y_pred))

             precision    recall  f1-score   support

           1       0.98      0.97      0.98       500
           2       0.98      0.97      0.97       500
           3       0.98      0.96      0.97       500
           4       0.97      0.98      0.97       500
           5       0.96      0.97      0.96       500
           6       0.99      0.98      0.98       500
           7       0.96      0.98      0.97       500
           8       0.97      0.98      0.98       500
           9       0.97      0.96      0.97       500
          10       0.99      1.00      0.99       500

    accuracy       0.97      0.97      0.97      5000
   macro avg       0.97      0.97      0.97      5000
weighted avg       0.97      0.97      0.97      5000

#可视化隐藏层

hidden_layer = thetafinal1[:, 1:] 
print(hidden_layer.shape)

fig, ax_array = plt.subplots(nrows=5, ncols=5, sharey=True, sharex=True, figsize=(12, 12))
for r in range(5):
    for c in range(5):
        ax_array[r, c].matshow(np.array(hidden_layer[5 * r + c].reshape((20, 20))),cmap=matplotlib.cm.binary)
        plt.xticks(np.array([]))
        plt.yticks(np.array([]))