# Multi-layer perceptron program is super slow

I'm writing a multi-layer perceptron from scratch and I think it's way slower than it should be. the culprit seems to be my compute_gradients-function, which according to my investigation answers for most of the execution time. It looks like this:

def compute_gradients(X, Y, S1, H, P, W1, W2, lamb):
# Y must be one-hot
# Y and P must be (10, n)
# X and H must be (3072, n)
# P is the softmax layer
if not (Y.shape[0] == 10 and P.shape[0] == 10 and Y.shape == P.shape):
raise ValueError("Y and P must have shape (10, k). Y: {}, P: {}".format(Y.shape, P.shape))
if not X.shape[0] == n_input:
raise ValueError("X must have shape ({}, k), has {}".format(n_input, X.shape))
if not H.shape[0] == n_hidden:
raise ValueError("H must have shape ({}, k)".format(n_hidden))

X, Y, H, P = X.T, Y.T, H.T, P.T
for x, y, s1, h, p in zip(X, Y, S1, H, P):
h = np.reshape(h, [1, n_hidden])
y = np.reshape(y, [10, 1])
p = np.reshape(p, [10, 1])

# Second layer
g = -(y-p).T
# First layer
g = np.matmul(g, W2)
ind = np.zeros(h.shape[1])
for i, val in enumerate(s1):
if val > 0:
ind[i] = 1
diag = np.diag(ind)

g = np.matmul(g, diag)

grad_W1 += np.matmul(g.T, np.reshape(x, [1, n_input]))
# Divide by batch size


If X, Y, H, P are all 10 rows long (n=10), the computations take about 1 second. This is way too much compared to my friends who are doing the same task. But I can't see any obvious inefficiencies in my code. What can I do to speed the computations up?

EDIT: Here's a runnable example. You also need the CIFAR dataset, found here: https://www.cs.toronto.edu/~kriz/cifar-100-python.tar.gz .

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

n_input = 3072
n_hidden = 3072

def one_hot(Y):
# assume Y = [1, 4, 9, 0, ...]
result = [None]*len(Y)
for i, cls in enumerate(Y):
onehot = {
0: lambda: [1,0,0,0,0,0,0,0,0,0],
1: lambda: [0,1,0,0,0,0,0,0,0,0],
2: lambda: [0,0,1,0,0,0,0,0,0,0],
3: lambda: [0,0,0,1,0,0,0,0,0,0],
4: lambda: [0,0,0,0,1,0,0,0,0,0],
5: lambda: [0,0,0,0,0,1,0,0,0,0],
6: lambda: [0,0,0,0,0,0,1,0,0,0],
7: lambda: [0,0,0,0,0,0,0,1,0,0],
8: lambda: [0,0,0,0,0,0,0,0,1,0],
9: lambda: [0,0,0,0,0,0,0,0,0,1],
}[cls]()
result[i] = onehot

result = np.array(result).T
return result

def unpickle(file):
import pickle
with open(file, "rb") as fo:
return d

path = "../../datasets/cifar/"

names = ["data_batch_1",
"data_batch_2",
"data_batch_3",
"data_batch_4",
"data_batch_5",
]

dataset_raw = unpickle(os.path.join(path, names[0]))
dataset_small = {}
dataset_small["data"] = dataset_raw[b"data"]/255
dataset_small["labels"] = one_hot(dataset_raw[b"labels"]).T

validation_set_raw = unpickle(os.path.join(path, names[1]))
validation_small = {}
validation_small["data"] = validation_set_raw[b"data"]/255
validation_small["labels"] = one_hot(validation_set_raw[b"labels"]).T

print("small dataset data shape: {}".format(dataset_small["data"].shape))
print("small dataset labels shape: {}".format(dataset_small["labels"].shape))
print("small validation data shape: {}".format(validation_small["data"].shape))
print("small validation labels shape: {}".format(validation_small["labels"].shape))

test_set_raw = unpickle(os.path.join(path, "test_batch"))
X_test = test_set_raw[b"data"]
Y_test = one_hot(test_set_raw[b"labels"])

# All data sets

dataset_large = {"data": np.zeros([0, 3072]), "labels": np.array([])}
validation_large = {}

## All data batches

for name in names[1:4]:
raw = unpickle(os.path.join(path, name))
dataset_large["data"] = np.append(dataset_large["data"], raw[b"data"], axis = 0)
dataset_large["labels"] = np.append(dataset_large["labels"], raw[b"labels"], axis = 0)
raw = unpickle(os.path.join(path, names[4]))
dataset_large["data"] = np.append(dataset_large["data"], raw[b"data"][0: -1000], axis = 0)
dataset_large["labels"] = np.append(dataset_large["labels"], raw[b"labels"][0: -1000], axis = 0)
validation_large["data"] = raw[b"data"][-1000: ]
validation_large["labels"] = raw[b"labels"][-1000: ]

# Make one-hot
dataset_large["labels"] = one_hot(dataset_large["labels"]).T
validation_large["labels"] = one_hot(validation_large["labels"]).T

# Normalize
dataset_large["data"] = dataset_large["data"]/255
validation_large["data"] = validation_large["data"]/255

print("large dataset data: {}".format(dataset_large["data"].shape))
print("large dataset labels: {}".format(dataset_large["labels"].shape))
print("large validation set data: {}".format(validation_large["data"].shape))
print("large validation set labels: {}".format(validation_large["labels"].shape))

def softmax(s):
numerator = np.exp(s)
denominator = np.matmul(np.ones(10).T, np.exp(s))
return numerator/denominator

def relu(s1):
h = np.maximum(np.zeros_like(s1), s1)
return h

def layer_1(X, W1, b1):
# WARNING: X must have shape (3072, n)
# W1 must have shape (3072, 3072)
# b1 must have shape (3072, 1)
# s will get shape (10, n)
# Return p with shape (10, n)

if not X.shape[0] == n_input:
raise ValueError("Got wrong shape of s1: {}".format(X.shape))

s1 = np.matmul(W1, X) + b1

if not s1.shape == (n_hidden, X.shape[1]):
raise ValueError("Got wrong shape of s1: {}".format(s1.shape))

h = relu(s1)

if not h.shape == (n_hidden, X.shape[1]):
raise ValueError("Got wrong shape of h: {}".format(h.shape))

return s1, h

def layer_2(h, W2, b2):
# WARNING: h must have shape (3072, n)
# W2 must have shape (10, 3072)
# b2 must have shape (10, 1)
# s will get shape (10, n)
# Returns p with shape (10, n)

s = np.matmul(W2, h) + b2

if not s.shape == (10, h.shape[1]):
raise ValueError("Got wrong shape of s: {}".format(s.shape))

p = softmax(s)

if not p.shape == (10, h.shape[1]):
raise ValueError("Got wrong shape of p: {}".format(p.shape))

return p

def evaluate_classifier(X, W1, W2, b1, b2):
if not X.shape[0] == n_input:
ValueError("Wrong shape of X: {}".format(X.shape))
if not len(X.shape) == 2:
ValueError("Wrong shape of X: {}".format(X.shape))
if not W1.shape == (n_hidden, n_input):
raise ValueError("Wrong shape of W1: {}".format(W1.shape))
if not b1.shape == (n_hidden, 1):
raise ValueError("Wrong shape of b1: {}".format(b1.shape))
if not W2.shape == (10, n_hidden):
raise ValueError("Wrong shape of W2: {}".format(W2.shape))
if not b2.shape == (10, 1):
raise ValueError("Wrong shape of b2: {}".format(b2.shape))

s1, h = layer_1(X, W1, b1)

p = layer_2(h, W2, b2)

return s1, h, p

def compute_cost(X, Y, W1, W2, b1, b2, lamb = 0):
# WARNING! Y must be one hot!
# X must have shape (3072, -1)
# Y and P must have shape (-1, 10)

if not Y.shape[1] == 10 and len(Y.shape) == 2:
raise ValueError("Wrong shape of Y: {}. Must be (-1, 10)".format(Y.shape))
if not X.shape[0] == 3072:
raise ValueError("Wrong shape of X: {}. Must be (3072, -1)".format(X.shape))
_, _, P = evaluate_classifier(X, W1, W2, b1, b2)
P = P.T

if not P.shape[1] == 10 and len(P.shape) == 2:
raise ValueError("Wrong shape of P: {}. Must be (-1, 10)".format(P.shape))
if not P.shape == Y.shape:
raise ValueError("Y and P must have the same shape. {} != {}".format(Y.shape, P.shape))

l_cross = 0
for i in range(X.shape[1]):
y = Y[i]
p = P[i]
y = np.reshape(y, [1, 10])
p = np.reshape(p, [10, 1])
l_cross += -np.log(np.matmul(y, p))[0][0]
avg_loss = l_cross/X.shape[1]
reg = lamb * (np.sum(W1) + np.sum(W2))

return avg_loss + reg

def compute_accuracy(X, Y, W1, W2, b1, b2):
# X must have shape (3072, n)
# Y must have shape (10, n)
if(not X.shape[0] == 3072):
raise ValueError("X must have shape (3072, n)")
try:
X.shape[1]
except IndexError:
raise(ValueError("X must have shape (3072, n)"))
if not Y.shape[0] == 10:
raise ValueError("Y must have shape (10, n)")
if not len(Y.shape) == 2:
raise ValueError("Y mus have shape (10, n)")
correct = 0
for x, y in zip(X.T, Y.T):
_, _, p = evaluate_classifier(np.reshape(x, (3072,1)), W1, W2, b1, b2)
pred = np.argmax(p)
corr = np.argmax(y)
correct += 1 if pred == corr else 0
return correct/X.shape[1]

def compute_gradients(X, Y, S1, H, P, W1, W2, lamb):
# Y must be one-hot
# Y and P must be (10, n)
# X and H must be (3072, n)
# P is the softmax layer
t = time.clock()
if not (Y.shape[0] == 10 and P.shape[0] == 10 and Y.shape == P.shape):
raise ValueError("Y and P must have shape (10, k). Y: {}, P: {}".format(Y.shape, P.shape))
if not X.shape[0] == n_input:
raise ValueError("X must have shape ({}, k), has {}".format(n_input, X.shape))
if not H.shape[0] == n_hidden:
raise ValueError("H must have shape ({}, k)".format(n_hidden))

X, Y, H, P = X.T, Y.T, H.T, P.T

for x, y, s1, h, p in zip(X, Y, S1, H, P):
h = np.reshape(h, [1, n_hidden])
y = np.reshape(y, [10, 1])
p = np.reshape(p, [10, 1])

# Second layer
g = -(y-p).T
# First layer
g = np.matmul(g, W2)
print(g.shape)
ind = np.zeros(h.shape[1])
for i, val in enumerate(s1):
if val > 0:
ind[i] = 1

diag = np.diag(ind)

g = np.matmul(g, diag)

grad_W1 += np.matmul(g.T, np.reshape(x, [1, n_input]))
# Divide by batch size
dur = time.clock() - t

def compute_gradient_num(X, Y, W1, W2, b1, b2, lamb = 0, h = 0.000001):
Y = Y.T

for i in range(b1.shape[0]):
b1_try = np.copy(b1)
b1_try[i,0] = b1[i,0] - h
c1 = compute_cost(X, Y, W1, W2, b1_try, b2, lamb)
b1_try = np.copy(b1)
b1_try[i, 0] = b1[i, 0] + h
c2 = compute_cost(X, Y, W1, W2, b1_try, b2, lamb)
grad_b1[i,0] = (c2-c1) / (2 * h)

for i in range(W1.shape[0]):
for j in range(W1.shape[1]):
W1_try = np.copy(W1)
W1_try[i,j] = W1[i,j] - h
c1 = compute_cost(X, Y, W1_try, W2, b1, b2, lamb)
W1_try = np.copy(W1)
W1_try[i, j] = W1[i, j] + h
c2 = compute_cost(X, Y, W1_try, W2, b1, b2, lamb)
grad_W1[i, j] = (c2 - c1) / (2 * h)

for i in range(b2.shape[0]):
b2_try = np.copy(b2)
b2_try[i,0] = b2[i,0] - h
c1 = compute_cost(X, Y, W1, W2, b1, b2_try, lamb)
b2_try = np.copy(b2)
b2_try[i, 0] = b2[i, 0] + h
c2 = compute_cost(X, Y, W1, W2, b1, b2_try, lamb)
grad_b2[i,0] = (c2-c1) / (2 * h)

for i in range(W2.shape[0]):
for j in range(W2.shape[1]):
W2_try = np.copy(W2)
W2_try[i,j] = W2[i,j] - h
c1 = compute_cost(X, Y, W1, W2_try, b1, b2, lamb)
W2_try = np.copy(W2)
W2_try[i, j] = W2[i, j] + h
c2 = compute_cost(X, Y, W1, W2_try, b1, b2, lamb)
grad_W2[i, j] = (c2 - c1) / (2 * h)

W1 = np.random.normal(0, 0.01, [n_hidden, n_input])
W2 = np.random.normal(0, 0.01, [10, n_hidden])
b1 = np.random.normal(0, 0.1, [n_hidden, 1])
b2 = np.random.normal(0, 0.1, [10, 1])

def train_batch(X, Y, W1, W2, b1, b2, gd_params, old_grads):
S1, H, P = evaluate_classifier(X.T, W1, W2, b1, b2)
try:
except KeyError:

def minibatch_GD(X, Y, W1, W2, b1, b2, gd_params, validation=None, verbose = True):
# X must have shape (n, 3072)
# Y must have shape (k, 10)
# validation["data"] must have shape (k, 3072)
# validation["labels"] must have shape (k, 10)
training_loss = [None]*gd_params["epochs"]
validation_loss = training_loss[:] # copying the list
training_accuracy = training_loss[:]
validation_accuracy = training_loss[:]

for i in range(gd_params["epochs"]):

if i % 10 == 0 and verbose:
print("Training epoch {} started".format(i+1))

# For each batch
for j in range(X.shape[0]//gd_params["batch_size"]):
start = j*gd_params["batch_size"]
end = start + gd_params["batch_size"]
batch = X[start: end]
targets = Y[start:end]

# Train the datapoints that didn't fit into last minibatch
n_rest = X.shape[0]%gd_params["batch_size"]
if n_rest > 0:
start = X.shape[0]-n_rest
end = X.shape[0]
batch = X[start: end]
targets = Y[start: end]

training_loss[i] = compute_cost(X.T, Y, W1, W2, b1, b2, gd_params["lambda"])
training_accuracy[i] = compute_accuracy(X.T, Y.T, W1, W2, b1, b2)
if validation:
validation_loss[i] = compute_cost(validation["data"].T, validation["labels"], W1, W2, b1, b2, gd_params["lambda"])
validation_accuracy[i] = compute_accuracy(validation["data"].T, validation["labels"].T, W1, W2, b1, b2)
return training_loss, None, training_accuracy, None

dataset_mini = {}
dataset_mini["data"] = dataset_small["data"][0:100]
dataset_mini["labels"] = dataset_small["labels"][0:100]

W1 = np.random.normal(0, 0.01, [n_hidden, n_input])
W2 = np.random.normal(0, 0.01, [10, n_hidden])
b1 = np.random.normal(0, 0.1, [n_hidden, 1])
b2 = np.random.normal(0, 0.1, [10, 1])

gd_params = {"epochs": 200, "batch_size":10, "lambda": 0, "eta": 0.01, "rho": 0.5}
X = dataset_mini["data"]
Y = dataset_mini["labels"]
training_loss, _, training_accuracy, _ = minibatch_GD(X, Y, W1, W2, b1, b2, gd_params)

plt.plot(training_loss)
plt.figure()
plt.plot(training_accuracy)

• Could you modify the question to add a runnable example? It would make it easier to investigate, because we can python -m cProfile your_example.py Apr 18, 2018 at 4:57
• On first glance the the culprit is that you don't use vectorized computation, but rather compute the gradient row by row in a for loop. Also you are using a linear output layer and a ReLu hidden layer correct? You can vectorize gradient computation of the ReLu instead of nesting for loops Apr 18, 2018 at 5:03
• @FirefoxMetzger, added runnable example Apr 18, 2018 at 6:59
• You are linking to the wrong data set. Your code seems to use CIFAR-10, but you link to CIFAR-100 which has a slightly different structure and won't work. Apr 18, 2018 at 9:57

Here is a faster version (around 300 ms per batch). Its essentially refactored for readability and often in vectorized code, more readability / code beauty makes code run faster:

def compute_gradients(X, Y, S1, H, P, W1, W2, lamb):
# Y must be one-hot
# Y and P must be (10, n)
# X and H must be (3072, n)
# P is the softmax layer

batch_size = X.shape[1]

t = time.time()

dp = -(Y-P).T
dl2 = dp
drelu = np.matmul(dl2, W2)
dl1 = drelu * (S1 > 0).T

dW2 = np.matmul(H,dl2).T/batch_size + 2*lamb*W2
dW1 = np.matmul(X,dl1).T/batch_size + 2*lamb*W1

db2 = np.sum(dl2)/batch_size
db1 = np.sum(dl1)/batch_size

dur = time.time() - t
return dW1, dW2, db1, db2


Things I changes:

• I threw out the assertions that you did in the beginning. Numpy will do the same and it really just slows you down. It also adds visual clutter and makes code hard to read
• I vectorized computing ReLu (comparing the value with 0 directly)
• I vectorized the gradient of the ReLu (matmul(a,diag(b)) is the same as a*b element wise)
• I vectorized the accumulation of the gradient update (this one took me a long time to learn as a Bachelor's; if you don't see it, try writing down the sum over an update batch and the dyadic product for each update on paper and refactor. The idea is that you do an outer product along one rank and an inner product along another)
• I changed time.clock() for time.time() to compute the runtime in seconds instead of 1/10th a second.

There is a lot more you could do to the rest of the code, but that is out of scope for this question.

• I don't think it's a bug actually, we're using the cross entropy loss and softmax for the last layer, relu for the hidden layer. The math ends up looking like that after some simplifications. Apr 18, 2018 at 12:03
• Ah, yes that's true. I didn't see the cross entropy term. Apr 18, 2018 at 12:58