Stuck implementing simple neural network
<p>I've been bashing my head against this brick wall for what seems like an eternity, and I just can't seem to wrap my head around it. I'm trying to implement an autoencoder using only numpy and matrix multiplication. No theano or keras tricks allowed.</p>

<p>I'll describe the problem and all its details. It is a bit complex at first since there are a lot of variables, but it really is quite straightforward.</p>

<p><strong>What we know</strong></p>

<p>1) <code>X</code> is an <code>m</code> by <code>n</code> matrix which is our inputs. The inputs are rows of this matrix. Each input is an <code>n</code> dimensional row vector, and we have <code>m</code> of them.</p>

<p>2)The number of neurons in our (single) hidden layer, which is <code>k</code>.</p>

<p>3) The activation function of our neurons (sigmoid, will be denoted as <code>g(x)</code>) and its derivative <code>g'(x)</code></p>

<p><strong>What we don't know and want to find</strong></p>

<p>Overall our goal is to find 6 matrices: <code>w1</code> which is <code>n</code> by <code>k</code>, <code>b1</code> which is <code>m</code> by <code>k</code>, <code>w2</code> which is <code>k</code> by <code>n</code>, b2 which is <code>m</code> by <code>n</code>, <code>w3</code> which is <code>n</code> by <code>n</code> and <code>b3</code> which is <code>m</code> by <code>n</code>.</p>

<p>They are initallized randomly and we find the best solution using gradient descent.</p>

<p><strong>The process</strong></p>

<p>The entire process looks something like this
<a href="http://i.stack.imgur.com/OQl3g.jpg"><img src="http://i.stack.imgur.com/OQl3g.jpg" alt="enter image description here"></a></p>

<p>First we compute <code>z1 = Xw1+b1</code>. It is <code>m</code> by <code>k</code> and is the input to our hidden layer. We then compute <code>h1 = g(z1)</code>, which is simply applying the sigmoid function to all elements of <code>z1</code>. naturally it is also <code>m</code> by <code>k</code> and is the output of our hidden layer.</p>

<p>We then compute <code>z2 = h1w2+b2</code> which is <code>m</code> by <code>n</code> and is the input to the output layer of our neural network. Then we compute <code>h2 = g(z2)</code> which again is naturally also <code>m</code> by <code>n</code> and is the output of our neural network.</p>

<p>Finally, we take this output and perform some linear operator on it: <code>Xhat = h2w3+b3</code> which is also <code>m</code> by <code>n</code> and is our final result.</p>

<p><strong>Where I am stuck</strong></p>

<p>The cost function I want to minimize is the mean squared error. I already implemented it in numpy code</p>

<pre><code>def cost(x, xhat):
    return (1.0/(2 * m)) * np.trace(np.dot(x-xhat,(x-xhat).T))
</code></pre>

<p>The problem is finding the derivatives of cost with respect to <code>w1,b1,w2,b2,w3,b3</code>. Let's call the cost <code>S</code>.</p>

<p>After deriving myself <strong>and checking myself numerically</strong>, I have established the following facts:</p>

<p>1) <code>dSdxhat = (1/m) * np.dot(xhat-x)</code></p>

<p>2) <code>dSdw3 = np.dot(h2.T,dSdxhat)</code></p>

<p>3) <code>dSdb3 = dSdxhat</code></p>

<p>4) <code>dSdh2 = np.dot(dSdxhat, w3.T)</code></p>

<p>But I can't for the life of me figure out dSdz2. It's a brick wall.</p>

<p>From chain-rule, it should be that dSdz2 = dSdh2 * dh2dz2 but the dimensions don't match. </p>

<p>What is the formula to compute the derivative of S with respect to z2?</p>

<p><strong>Edit</strong> - This is my code for the entire feed forward operation of the autoencoder.</p>

<pre><code>import numpy as np

def g(x): #sigmoid activation functions
    return 1/(1+np.exp(-x)) #same shape as x!

def gGradient(x): #gradient of sigmoid
    return g(x)*(1-g(x)) #same shape as x!

def cost(x, xhat): #mean squared error between x the data and xhat the output of the machine
    return (1.0/(2 * m)) * np.trace(np.dot(x-xhat,(x-xhat).T))

#Just small random numbers so we can test that it's working small scale
m = 5 #num of examples
n = 2 #num of features in each example
k = 2 #num of neurons in the hidden layer of the autoencoder
x = np.random.rand(m, n) #the data, shape (m, n)

w1 = np.random.rand(n, k) #weights from input layer to hidden layer, shape (n, k)
b1 = np.random.rand(m, k) #bias term from input layer to hidden layer (m, k)
z1 = np.dot(x,w1)+b1 #output of the input layer, shape (m, k)
h1 = g(z1) #input of hidden layer, shape (m, k)

w2 = np.random.rand(k, n) #weights from hidden layer to output layer of the autoencoder, shape (k, n)
b2 = np.random.rand(m, n) #bias term from hidden layer to output layer of autoencoder, shape (m, n)
z2 = np.dot(h1, w2)+b2 #output of the hidden layer, shape (m, n)
h2 = g(z2) #Output of the entire autoencoder. The output layer of the autoencoder. shape (m, n)

w3 = np.random.rand(n, n) #weights from output layer of autoencoder to entire output of the machine, shape (n, n)
b3 = np.random.rand(m, n) #bias term from output layer of autoencoder to entire output of the machine, shape (m, n)
xhat = np.dot(h2, w3)+b3 #the output of the machine, which hopefully resembles the original data x, shape (m, n)
</code></pre>

<p>OK, here's a suggestion.  In the vector case, if you have <em>x</em> as a vector of length <code>n</code>, then <code>g(x)</code> is also a vector of length <code>n</code>.  However, <code>g'(x)</code> is not a vector, it's the <a href="https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant" rel="nofollow">Jacobian matrix</a>, and will be of size <code>n X n</code>.  Similarly, in the minibatch case, where <em>X</em> is a matrix of size <code>m X n</code>, <code>g(X)</code> is <code>m X n</code> but <code>g'(X)</code> is <code>n X n</code>.  Try:</p>

<pre><code>def gGradient(x): #gradient of sigmoid
    return np.dot(g(x).T, 1 - g(x))
</code></pre>

<p>@Paul is right that the bias terms should be vectors, not matrices.  You should have:</p>

<pre><code>b1 = np.random.rand(k) #bias term from input layer to hidden layer (k,)
b2 = np.random.rand(n) #bias term from hidden layer to output layer of autoencoder, shape (n,)
b3 = np.random.rand(n) #bias term from output layer of autoencoder to entire output of the machine, shape (n,)
</code></pre>

<p>Numpy's broadcasting means that you don't have to change your calculation of <code>xhat</code>.</p>

<p>Then (I think!) you can compute the derivatives like this:</p>

<pre><code>dSdxhat = (1/float(m)) * (xhat-x)
dSdw3 = np.dot(h2.T,dSdxhat)
dSdb3 = dSdxhat.mean(axis=0)
dSdh2 = np.dot(dSdxhat, w3.T)
dSdz2 = np.dot(dSdh2, gGradient(z2))
dSdb2 = dSdz2.mean(axis=0)
dSdw2 = np.dot(h1.T,dSdz2)
dSdh1 = np.dot(dSdz2, w2.T)
dSdz1 = np.dot(dSdh1, gGradient(z1))
dSdb1 = dSdz1.mean(axis=0)
dSdw1 = np.dot(x.T,dSdz1)
</code></pre>

<p>Does this work for you?</p>

<p><strong>Edit</strong></p>

<p>I've decided that I'm not at all sure that <code>gGradient</code> is supposed to be a matrix.  How about:</p>

<pre><code>dSdxhat = (xhat-x) / m
dSdw3 = np.dot(h2.T,dSdxhat)
dSdb3 = dSdxhat.sum(axis=0)
dSdh2 = np.dot(dSdxhat, w3.T)
dSdz2 = h2 * (1-h2) * dSdh2
dSdb2 = dSdz2.sum(axis=0)
dSdw2 = np.dot(h1.T,dSdz2)
dSdh1 = np.dot(dSdz2, w2.T)
dSdz1 = h1 * (1-h1) * dSdh1
dSdb1 = dSdz1.sum(axis=0)
dSdw1 = np.dot(x.T,dSdz1)
</code></pre>

python algorithm numpy matrix derivative