梯度检测没有实现。有借鉴网上的部分

导入相关库，读取数据

因为这次的数据是mat文件，需要使用scipy库中的loadmat进行读取数据。

通过对数据类型的分析，发现是字典类型，查看该字典的键，可以发现又X，y等关键字。

import numpy as np
import scipy.io as sio
from scipy.optimize import minimize
from sklearn.preprocessing import OneHotEncoder#读取数据
path = "./ex4data1.mat"
data = sio.loadmat(path)
# print(type(data))
# print(data.keys())
X = data.get("X")
Y = data.get("y").flatten()
# X = np.insert(X,0,values=1,axis=1)
# print(Y.shape)
# print(Y)

one-hot编码

在之前没有涉及神经网络的方向传播等时，可以直接将Y降维。

而在神经网络中需要对Y进行编码，使每一个y值都是一个10维的向量

#one-hot编码
encoder = OneHotEncoder(sparse=False)
Y_onehot = encoder.fit_transform(Y.reshape(-1,1))

随机初始化

当使用高级优化算法或者梯度下降算法时，需要对向量进行赋值。

第一种想法是全部设为0，在逻辑回归中，是完全被允许的，但在训练神经网络时，这种初始值起不到任何作用。因为这会导致第二层激活项的值完全相同（因为前一层的权重相同），这个问题称为对称权重问题。

为了解决这个问题，我们采用随机初始化的方法，也就是将每一个都初始化为一个范围在中的一个随机数。

#初始化参数
input_size = 400
hidden_size = 25
num_labels = 10
lamda = 1size = (input_size+1)*hidden_size+(hidden_size+1)*num_labels
params = np.random.uniform(-1.2,1.2,size)

序列化

因为到后面优化的时候fun和jac参数要求代价函数和梯度的第一个参数必须为一维向量，所以这里需要序列化。

#序列化
def serialize(theta1,theta2):return np.r_[theta1.flatten(),theta2.flatten()]
def deserialize(params):return params[:(input_size+1)*hidden_size].reshape(hidden_size,input_size+1),params[(input_size+1)*hidden_size:].reshape(num_labels,hidden_size+1)theta1,theta2 = deserialize(params)

前向传播

同上一题一致的前向传播

#前向传播
def feed_forward(params,X):theta1, theta2 = deserialize(params)a1 = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)z2 = a1 @ theta1.Ta2 = 1/(1+np.exp(-z2))a2 = np.insert(a2,0,values=1,axis=1)z3 = a2 @ theta2.Th = 1/(1+np.exp(-z3))return a1,z2,a2,z3,h

代价函数

不需要正则化

公式如下：

#代价函数
def cost_func(params,X,Y,lamda):theta1, theta2 = deserialize(params)a1,z2, a2, z3, h = feed_forward(params,X)m = len(X)cost = -np.sum(Y * np.log(h) + (1-Y) * np.log(1-h))/mreg = (np.sum(np.power(theta1[:,1:],2)) + np.sum(np.power(theta2[:,1:],2)))* lamda/(2*m)return reg+cost

反向传播

相关公式如下：

#反向传播
def sigmoid_gradient(z):return 1/(1+np.exp(-z))*(1-1/(1+np.exp(-z)))def gradientReg(params,X,Y,lamda=1):theta1, theta2 = deserialize(params)m = len(X)a1, z2, a2, z3, h = feed_forward(params,X)delta3 = h - Ydelta2 = delta3 @ theta2[:,1:] * sigmoid_gradient(z2)Delta2 = delta3.T @ a2 / mDelta1 = delta2.T @ a1 / mtheta1[:, 0] = 0theta2[:, 0] = 0regDelta1 = Delta1 + (lamda / m) * theta1regDelta2 = Delta2 + (lamda / m) * theta2return serialize(regDelta1,regDelta2)

优化

fmin = minimize(fun=cost_func,x0=params,args=(X,Y_onehot,lamda),method="TNC",jac=gradientReg,options={'maxiter': 400})

预测

a1,z2,a2,z3,h = feed_forward(fmin.x,X)
y_pred = np.argmax(h,axis=1)
y_pred = y_pred + 1
acc = np.mean(y_pred==Y)
print(acc)

准确率达到了98.64%

---------------------------------------------------------------------------------------------------------------------------------

可以先将数据可视化，展示手写数字。

def one_image(X):pick_one = np.random.randint(5000)image = X[pick_one,:]fig,ax = plt.subplots()ax.imshow(image.reshape(20,20).T,cmap="gray_r")plt.xticks([])plt.yticks([])plt.show()one_image(X)
def more_image(X):pick_more = np.random.choice(len(X),100)images = X[pick_more,:]fig,ax = plt.subplots(nrows=10,ncols=10, figsize=(8, 8), sharex=True, sharey=True)for i in range(10):for j in range(10):ax[i,j].imshow(images[10*i+j].reshape(20,20).T,cmap="gray_r")plt.xticks([])plt.yticks([])plt.show()
more_image(X)

最后可以使用上面写出的函数将隐藏层可视化。

thetafinal1, thetafinal2 = deserialize(fmin.x)
hidden_layer = thetafinal1[:, 1:] 
more_image(hidden_layer)