gdbackprop.py

import numpy as np

def sigmoid(X):
    return 1. / (1. + np.exp(-X))

def sigmoid_deriv(X):
    return sigmoid(X) * (1. - sigmoid(X))

def vec(X):
    return X.T.reshape(-1, 1)

def vecinv(m, n, X):
    return X.reshape(m, n, order='F') # Fortran-like indexing to reshape by column order

X = np.array([[1, 1], [1, 0], [0, 1], [0, 0]])
Y = np.array([[0], [1], [1], [0]])

W = 2 * np.random.random((2, 4)) - 1
V = 2 * np.random.random((4, 1)) - 1

print ((sigmoid(sigmoid(X.dot(W)).dot(V)) - Y) ** 2.) / (2 * 4.)

for i in range(1000):
    forward_hidden = sigmoid(X.dot(W))
    forward_output = sigmoid(forward_hidden.dot(V))

    dcost = np.diag(vec(forward_output - Y)[:, 0])
    dsderiv_ol = np.diag(vec(sigmoid_deriv(forward_hidden.dot(V)))[:, 0])
    err = dcost * dsderiv_ol
    dv = np.kron(np.eye(V.shape[1]).T, forward_hidden)
    delta_V = err.dot(dv)

    dmatmulderiv_ol = np.kron(V.T, np.eye(forward_hidden.shape[0]))
    dsderiv_hl = np.diag(vec(sigmoid_deriv(X.dot(W)))[:, 0])
    err = err.dot(dmatmulderiv_ol).dot(dsderiv_hl)
    dw = np.kron(np.eye(W.shape[1]).T, X) #transpose useless
    delta_W = err.dot(dw)

    ####### MATRIX VERSION IN ORDER TO SEE IF THE GRADIENTS ARE THE SAME
    merr = (forward_output - Y) * (sigmoid_deriv(forward_hidden.dot(V)))
    mdelta_V = forward_hidden.T.dot(merr)
    merr = merr.dot(V.T) * sigmoid_deriv(X.dot(W))
    mdelta_W = X.T.dot(merr)
    mV = V - mdelta_V
    mW = W - mdelta_W
    #####

    V = V -  vecinv(V.shape[0], V.shape[1], (np.ones((1, delta_V.shape[0])).dot(delta_V)).T)
    W = W -  vecinv(W.shape[0], W.shape[1], (np.ones((1, delta_W.shape[0])).dot(delta_W)).T)

    print 'Error at step %d : %f' % (i, np.sum(((sigmoid(sigmoid(X.dot(W)).dot(V)) - Y) ** 2.) / (2 * 4.)))
    print sigmoid(sigmoid(X.dot(W)).dot(V))