how to calculate errors in neural networks, when the size of the matrix dimensions of the output do not match the size of the label?

I tried to make code for image classification with neural networks without using libraries such as KERAS, etc.

but I have difficulty when I want to calculate the error, because the size of the matrix of the output feedforward does not match the size of the matrix label one hot multiclass

The following is the code for the one hot multiclass label, here I use 20 classes with each class having 15 variations (300 inputs, size = 64 x 64 pixels) :

labels_train = np.array([0] * 15 + [1] * 15 + [2] * 15 + [3] * 15 + [4] * 15 + [5] * 15 + [6] * 15 + [7] 
                        * 15 + [8] * 15 + [9] * 15 + [10] * 15 + [11] * 15 + [12] * 15 + [13] * 15 + [14] 
                        * 15 + [15] * 15 + [16] * 15 + [17] * 15 + [18] * 15 + [19] * 15)

y_train = np.zeros((300,20))

for i in range(300):
    y_train[i, labels_train[i]] = 1

The following is the code for feedforward and calculation of the error value :

for epoch in range(1):

    hasil_ih = np.dot(W_ih, x_train) + B_ih
    output_ih = sigmoid(hasil_ih)

    hasil_hh = np.dot(W_hh, output_ih) + B_hh
    output_hh = sigmoid(hasil_hh)

    hasil_ho = np.dot(W_ho, output_hh) + B_ho
    output_ho = sigmoid(hasil_ho)

    print(output_ho.shape)
    print(y_train.shape)
    error_o = labels_train - output_ho

and here are the results when the program is running :

Traceback (most recent call last):
(300, 4096)
  File "D:/Pengenalan_Induk_Aksara/NN2.py", line 71, in <module>
(300, 20)
    error_o = labels_train - output_ho
ValueError: operands could not be broadcast together with shapes (300,) (300,4096) 

Process finished with exit code 1

please help me, I'm still a beginner and I haven't found a solution. The following is the complete code :

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import os


dir = "data_train/LBP_R1/"
kategori = ["ka","ga","nga","pa","ba","ma","ta","da","na","ca","ja","nya","ya","a","la","ra","sa","wa","ha","gha"]
features = []; label = []

for y in kategori:
    path = os.path.join(dir, y)
    class_num = kategori.index(y)
    label.append(class_num)

    for img in os.listdir(path):
        im = Image.open(os.path.join(path,img))
        imgs = list(im.getdata())
        features.append(imgs)


x_train = np.vstack([features])
x_train = (x_train/255)


labels_train = np.array([0] * 15 + [1] * 15 + [2] * 15 + [3] * 15 + [4] * 15 + 
                        [5] * 15 + [6] * 15 + [7] * 15 + [8] * 15 + [9] * 15 + 
                        [10] * 15 + [11] * 15 + [12] * 15 + [13] * 15 + [14] * 
                        15 + [15] * 15 + [16] * 15 + [17] * 15 + [18] * 15 + [19] 
                        * 15)

y_train = np.zeros((300,20))

for i in range(300):
    y_train[i, labels_train[i]] = 1

#Training Phase
def sigmoid(x):
    return 1/(1+np.exp(-x)) 

def dsigmoid(y):
    return y*(1-y) 

attributes = x_train.shape[0] 
hidden_nodes = 4 
output_nodes = 300 

np.random.seed(42)
W_ih = np.random.rand(hidden_nodes,attributes) #bobot input layer 1
B_ih = np.random.rand(hidden_nodes,1) #Bias input layer 1

W_hh = np.random.rand(hidden_nodes, hidden_nodes) #weight hidden layer 2
B_hh = np.random.rand(hidden_nodes, 1) #bias hidden layer 2

W_ho = np.random.rand(output_nodes, hidden_nodes) #weight hidden output
B_ho = np.random.rand(output_nodes, 1) #bias hidden output

lr = 0.25; #Learning rate

for epoch in range(1):

    hasil_ih = np.dot(W_ih, x_train) + B_ih
    output_ih = sigmoid(hasil_ih)

    hasil_hh = np.dot(W_hh, output_ih) + B_hh
    output_hh = sigmoid(hasil_hh)

    hasil_ho = np.dot(W_ho, output_hh) + B_ho
    output_ho = sigmoid(hasil_ho)

    print(output_ho.shape)
    print(y_train.shape)

    error_o = labels_train - output_ho
    print(output_ho)
    print(y_train)
    print(error_o)