Correct way of handling the output layer deltas in neural networks

According to Andrew Ng's Coursera Machine Learning course (notes here and here), a training algorithm for neural nets is this:

1. $\delta^{(L)} = a^{(L)} - y$
2. $\text{for }l = L-1, ..., 2\text{ do}$
3. $\quad \delta^{(l)} = (\Theta^{(l)})^T \delta^{(l + 1)} .* g'(z^{(l)})$
4. $\text{end}$

Then, for computing the actual gradients:

1. $\text{Set }\Delta^{(l)} = 0,\text{ for }l = 1, ..., L-1$
2. $\text{for }i = 1, ..., m\text{ do}$
3. $\quad \text{Set }a^{(1)} = x^{(i)}$
4. $\quad \text{Compute }a^{(l)}\text{ for }l=2,...,L \text{ using forward propagation (Algorithm 1)}$
5. $\quad \text{Using }y^{(i)}, \text{compute } \delta^{(l)} \text{ for } l = L, ..., 2 \\ \qquad\text{ using the backpropagation algorithm (Algorithm 2)}$
6. $\quad \Delta^{(l)} = \Delta^{(l)} + \delta^{(l+1)}(a^{(l)})^T \text{ for }l=L-1, ..., 1$
7. $\text{end}$

However, this seems to contradict this (page 7, taken from Mitchell's book I believe, but also in various other tutorials: here and here for example), where the output layer delta is defined as:

$${\delta^{(L)} = (a^{(L)} - y}).*(a^{(L)}.*(1-a^{(L)}))$$

So the derivative of the output layer outputs is also factored in. In the above, $a$ always represents the output of a layer and $z$ is the "raw computation" of a layer, before the sigmoid function is applied to it. $.*$ is element-wise multiplication.

According to what I learned and understand, you should factor in the derivative as well. According to my experiments, things work a little better if you do, but I haven't been able to test too much.

Which is the correct way? Or are they equivalent but I misunderstood something about the first notation (and therefore also implemented it wrong)? Or it doesn't matter and I'll just need proper hyperparameters for each to obtain similar results?

Here is some code that shows what I mean:

X = np.array([ [0,0,1],[0,1,1],[1,0,1],[1,1,1] ])
y = np.array([[0,1,1,0]]).T
syn0 = 2*np.random.random((3,4)) - 1
syn1 = 2*np.random.random((4,1)) - 1
for j in xrange(60000):
l1 = 1/(1+np.exp(-(np.dot(X,syn0))))
l2 = 1/(1+np.exp(-(np.dot(l1,syn1))))
l2_delta = (y - l2)*(l2*(1-l2)) #######################
l1_delta = l2_delta.dot(syn1.T) * (l1 * (1-l1))
syn1 += l1.T.dot(l2_delta)
syn0 += X.T.dot(l1_delta)


According to my understanding, Andrew Ng's version suggests replacing the line marked with hashtags with just this:

l2_delta = y - l2


Also according to my understanding, doing that would be wrong.

I think I figured out my mistake.

Both versions are correct under the right loss function.

Andrew Ng's version uses log loss:

$$J(\theta)=-\frac{1}{N}\sum_{i=1}^{N}\Big(y_ilog\big(h_\theta(x_i)\big)+(1-y_i)log\big(1-h_\theta(x_i)\big)\Big)$$

Where:

$$h_\theta(x) = g(x\theta) = \frac{1}{1+exp(-x\theta)}$$

Let's take the derivative of this (I'm shortening this a little by dropping the sum) with respect to $\theta$:

\begin{align} &\frac{\partial}{\partial \theta} \Big(ylog\big(h_\theta(x)\big)+(1-y)log\big(1-h_\theta(x)\big)\Big) \\ &=y\frac{1}{h_\theta(x)}g(x\theta)\big(1-g(x\theta)\big)x+(1-y)\frac{1}{1-h_\theta(x)}g(x\theta)\big(g(x\theta)-1\big)x \\ &=y\big(1-g(x\theta)\big)x-(1-y)g(x\theta)x \\ &=yx-yg(x\theta)x-g(x\theta)x+yg(x\theta)x \\ &=\big(y-g(x\theta)\big)x \end{align}

Mitchell's version uses squared loss, which is defined along these lines:

$$J(\theta)=-\frac{1}{2N}\sum_{i=1}^{N}\Big(y_i-h_\theta(x_i)\Big)^2$$

Taking the derivative with respect to $\theta$ as we did above, we get:

\begin{align} &\frac{\partial}{\partial \theta}\Big(y-h_\theta(x)\Big)^2\\ &=2\big(y-h_\theta(x)\big)g(x\theta)\big(g(x\theta)-1\big)x \\ &=2\big(h_\theta(x)-y\big)g(x\theta)\big(1-g(x\theta)\big)x \end{align}

I missed the fact that they were using different loss functions. Both versions are correct.