How to correctly calculate gradients in neural network with numpy

Viewed 2560

I am trying to build a simple neural network class from scratch using numpy, and test it using the XOR problem. But the backpropagation function (backprop) does not seem to be working correctly.

In the class, I construct instances by passing in the size of each layer, and the activation functions to use at each layer. I assume that the final activation function is softmax, so that I can calculate the derivative of cross-entropy loss wrt to Z of the last layer. I also do not have a separate set of bias matrices in my class. I just include them in the weight matrices as an extra column at the end.

I know that my backprop function is not working correctly, because the neural network does not ever converge on a somewhat correct output. I also created a numerical gradient function, and when comparing the results of both. I get drastically different numbers.

My understanding from what I have read is that the delta values of each layer (with L being the last layer, and i representing any other layer) should be:

delta_L = (dE/dZ)

delta_i = WT(i+1)dotdelta(i+1)*sigmoid_prime(zi)

And the respective gradients/weight-update of those layers should be:

gradient_l = (delta_l)(activations(l-1).T)

W = W - alpha(W)(dE/dW)

Where * is the hardamard product, a represents the activation of some layer, and z represents the nonactivated output of some layer.

The sample data that I am using to test this is at the bottom of the file.

This is my first time trying to implement the backpropagation algorithm from scratch. So I am a bit lost on where to go from here.

import numpy as np

def sigmoid(n, deriv=False):
    if deriv:
        return np.multiply(n, np.subtract(1, n))
    return 1 / (1 + np.exp(-n))

def softmax(X, deriv=False):
    if not deriv:
        exps = np.exp(X - np.max(X))
        return exps / np.sum(exps)
    else:
        raise Error('Unimplemented')

def cross_entropy(y, p, deriv=False):
    """
    when deriv = True, returns deriv of cost wrt z
    """
    if deriv:
        ret = p - y
        return ret
    else:
        p = np.clip(p, 1e-12, 1. - 1e-12)
        N = p.shape[0]
        return -np.sum(y*np.log(p))/(N)

class NN:
    def __init__(self, layers, activations):
        """random initialization of weights/biases
        NOTE - biases are built into the standard weight matrices by adding an extra column
        and multiplying it by one in every layer"""
        self.activate_fns = activations
        self.weights = [np.random.rand(layers[1], layers[0]+1)]
        for i in range(1, len(layers)):
            if i != len(layers)-1:
                self.weights.append(np.random.rand(layers[i+1], layers[i]+1))

                for j in range(layers[i+1]):
                    for k in range(layers[i]+1):
                        if np.random.rand(1,1)[0,0] > .5:
                            self.weights[-1][j,k] = -self.weights[-1][j,k]

    def ff(self, X, get_activations=False):
         """Feedforward"""
         activations, zs = [], []
         for activate, w in zip(self.activate_fns, self.weights):
             X = np.vstack([X, np.ones((1, 1))]) # adding bias
             z = w.dot(X)
             X = activate(z)
             if get_activations:
                 zs.append(z)
                 activations.append(X)
         return (activations, zs) if get_activations else X

    def grad_descent(self, data, epochs, learning_rate):
        """gradient descent
        data - list of 2 item tuples, the first item being an input, and the second being its label"""
        grad_w = [np.zeros_like(w) for w in self.weights]
        for _ in range(epochs):
            for x, y in data:
                grad_w = [n+o for n, o in zip(self.backprop(x, y), grad_w)]
            self.weights = [w-(learning_rate/len(data))*gw for w, gw in zip(self.weights, grad_w)]

    def backprop(self, X, y):
        """perfoms backprop for one layer of a NN with softmax/cross_entropy output layer"""
        (activations, zs) = self.ff(X, True)
        activations.insert(0, X)

        deltas = [0 for _ in range(len(self.weights))]
        grad_w = [0 for _ in range(len(self.weights))]
        deltas[-1] = cross_entropy(y, activations[-1], True) # assumes output activation is softmax
        grad_w[-1] = np.dot(deltas[-1], np.vstack([activations[-2], np.ones((1, 1))]).transpose())
        for i in range(len(self.weights)-2, -1, -1):
            deltas[i] = np.dot(self.weights[i+1][:, :-1].transpose(), deltas[i+1]) * self.activate_fns[i](zs[i], True)
            grad_w[i] = np.hstack((np.dot(deltas[i], activations[max(0, i-1)].transpose()), deltas[i]))

        # check gradient
        num_gw = self.gradient_check(X, y, i)
        print('numerical:', num_gw, '\nanalytic:', grad_w)

        return grad_w

    def gradient_check(self, x, y, i, epsilon=1e-4):
        """Numerically calculate the gradient in order to check analytical correctness"""
        grad_w = [np.zeros_like(w) for w in self.weights]
        for w, gw in zip(self.weights, grad_w):
            for j in range(w.shape[0]):
                for k in range(w.shape[1]):
                    w[j,k] += epsilon
                    out1 = cross_entropy(self.ff(x), y)
                    w[j,k] -= 2*epsilon
                    out2 = cross_entropy(self.ff(x), y)
                    gw[j,k] = np.float64(out1 - out2) / (2*epsilon)
                    w[j,k] += epsilon # return weight to original value
        return grad_w

##### TESTING #####
X = [np.array([[0],[0]]), np.array([[0],[1]]), np.array([[1],[0]]), np.array([[1],[1]])]
y = [np.array([[1], [0]]), np.array([[0], [1]]), np.array([[0], [1]]), np.array([[1], [0]])]
data = []
for x, t in zip(X, y):
    data.append((x, t))

def nn_test():
    c = NN([2, 2, 2], [sigmoid, sigmoid, softmax])
    c.grad_descent(data, 100, .01)
    for x in X:
        print(c.ff(x))
nn_test()

UPDATE: I found one small bug in the code, but it still does not converge correctly. I calculated/derived the gradients for both matrices by hand and found no errors in my implementation, so I still do not know what is wrong with it.

UPDATE #2: I created a procedural version of what I was using above with the following code. Upon testing I discovered that the NN was able to learn the correct weights for classifying each of the 4 cases in XOR separately, but when I try to train using all the training examples at once (as shown), the resultant weights almost always output something around .5 for both output nodes. Could someone please tell me why this is occurring?

X = [np.array([[0],[0]]), np.array([[0],[1]]), np.array([[1],[0]]), np.array([[1],[1]])]
y = [np.array([[1], [0]]), np.array([[0], [1]]), np.array([[0], [1]]), np.array([[1], [0]])]
weights = [np.random.rand(2, 3) for _ in range(2)]
for _ in range(1000):
    for i in range(4):
        #Feedforward
        a0 = X[i]
        z0 = weights[0].dot(np.vstack([a0, np.ones((1, 1))]))
        a1 = sigmoid(z0)
        z1 = weights[1].dot(np.vstack([a1, np.ones((1, 1))]))
        a2 = softmax(z1)
        # print('output:', a2, '\ncost:', cross_entropy(y[i], a2))

        #backprop
        del1 = cross_entropy(y[i], a2, True)
        dcdw1 = del1.dot(np.vstack([a1, np.ones((1, 1))]).T)
        del0 = weights[1][:, :-1].T.dot(del1)*sigmoid(z0, True)
        dcdw0 = del0.dot(np.vstack([a0, np.ones((1, 1))]).T)

        weights[0] -= .03*weights[0]*dcdw0
        weights[1] -= .03*weights[1]*dcdw1
i = 0
a0 = X[i]
z0 = weights[0].dot(np.vstack([a0, np.ones((1, 1))]))
a1 = sigmoid(z0)
z1 = weights[1].dot(np.vstack([a1, np.ones((1, 1))]))
a2 = softmax(z1)
print(a2)
1 Answers

Softmax doesn't look right

Using cross entropy loss, the derivative for softmax is really nice (assuming you are using a 1 hot vector, where "1 hot" essentially means an array of all 0's except for a single 1, ie: [0,0,0,0,0,0,1,0,0])

For node y_n it ends up being y_n-t_n. So for a softmax with output:

[0.2,0.2,0.3,0.3]

And desired output:

[0,1,0,0]

The gradient at each of the softmax nodes is:

[0.2,-0.8,0.3,0.3]

It looks as if you are subtracting 1 from the entire array. The variable names aren't very clear, so if you could possibly rename them from L to what L represents, such as output_layer I'd be able to help more.

Also, for the other layers just to clear things up. When you say a^(L-1) as an example, do you mean "a to the power of (l-1)" or do you mean "a xor (l-1)"? Because in python ^ means xor.

EDIT:

I used this code and found the strange matrix dimensions (modified at line 69 in the function backprop)

deltas = [0 for _ in range(len(self.weights))]
grad_w = [0 for _ in range(len(self.weights))]
deltas[-1] = cross_entropy(y, activations[-1], True) # assumes output activation is softmax
print(deltas[-1].shape)
grad_w[-1] = np.dot(deltas[-1], np.vstack([activations[-2], np.ones((1, 1))]).transpose())
print(self.weights[-1].shape)
print(activations[-2].shape)
exit()
Related