BP neural network algorithm realizes handwritten numeral recognition

import numpy as np
from sklearn.datasets import load_digits
from sklearn.preprocessing import LabelBinarizer
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt

# load dataset
digits = load_digits()
print(digits.images.shape)
# show image
plt.imshow(digits.images[20],cmap='gray')
plt.show()

# data
X = digits.data
# label
y = digits.target
print(X.shape)
print(y.shape)

#build a network model 64-100-10
# definite input_layer to hidden_layer
V = np.random.random((64, 100))*2-1
# definite hidden_layer to output_layer
W = np.random.random((100, 10))*2-1

# split dataset
# default 1/4 test dataset, 3/4 train dataset
X_train, X_test, y_train, y_test = train_test_split(X,y)

# one-hot labels
labels_train = LabelBinarizer().fit_transform(y_train)
print(y_train[:5])
print(labels_train[:5])

#activation function
def sigmoid(x):
    return 1/(1+np.exp(-x))

# delta
def dsigmoid(x):
    return x*(1-x)

# train model
def train(X,y,steps=10000,learning_rate=0.11):
    global V,W
    for n in range(steps):
        # stochastic choice one data
        i = np.random.randint(X.shape[0])
        # get one data
        x = X[i]
        x = np.atleast_2d(x)
        # BP algorithm
        # compute the output of hidden_layer
        L1 = sigmoid(np.dot(x,V))
        # compute the output of output_layer
        L2 = sigmoid(np.dot(L1,W))
        dL2 = (y[i] - L2)*dsigmoid(L2)
        dL1 = dL2.dot(W.T)*dsigmoid(L1)

        # update weights
        W += learning_rate * L1.T.dot(dL2)
        V += learning_rate * x.T.dot(dL1)

        # predict accuracy once of 1000 steps
        if (n+1)%1000 == 0:
            output = predict(X_test)
            predictions = np.argmax(output, axis=1)
            acc = np.mean(np.equal(predictions, y_test))
            print("steps: {}, accuracy: {}".format(n+1, acc))

def predict(x):
    L1 = sigmoid(np.dot(x,V))
    L2 = sigmoid(np.dot(L1,W))
    return L2

train(X_train, labels_train, 30000)
predictions = np.argmax(predict(X_test),axis=1)
print(classification_report(predictions, y_test))
print(confusion_matrix(predictions,y_test))