Backpropagation with Softmax / Cross Entropy I'm trying to understand how backpropagation works for a softmax/cross-entropy output layer.
The cross entropy error function is
$$E(t,o)=-\sum_j t_j \log o_j$$
with $t$ and $o$ as the target and output at neuron $j$, respectively. The sum is over each neuron in the output layer. $o_j$ itself is the result of the softmax function:
$$o_j=softmax(z_j)=\frac{e^{z_j}}{\sum_j e^{z_j}}$$
Again, the sum is over each neuron in the output layer and $z_j$ is the input to neuron $j$:
$$z_j=\sum_i w_{ij}o_i+b$$
That is the sum over all neurons in the previous layer with their corresponding output $o_i$ and weight $w_{ij}$ towards neuron $j$ plus a bias $b$.
Now, to update a weight $w_{ij}$ that connects a neuron $j$ in the output layer with a neuron $i$ in the previous layer, I need to calculate the partial derivative of the error function using the chain rule:
$$\frac{\partial E} {\partial w_{ij}}=\frac{\partial E} {\partial o_j} \frac{\partial o_j} {\partial z_{j}} \frac{\partial z_j} {\partial w_{ij}}$$
with $z_j$ as the input to neuron $j$.
The last term is quite simple. Since there's only one weight between $i$ and $j$, the derivative is:
$$\frac{\partial z_j} {\partial w_{ij}}=o_i$$
The first term is the derivation of the error function with respect to the output $o_j$:
$$\frac{\partial E} {\partial o_j} = \frac{-t_j}{o_j}$$
The middle term is the derivation of the softmax function with respect to its input $z_j$ is harder:
$$\frac{\partial o_j} {\partial z_{j}}=\frac{\partial} {\partial z_{j}} \frac{e^{z_j}}{\sum_j e^{z_j}}$$
Let's say we have three output neurons corresponding to the classes $a,b,c$ then $o_b = softmax(b)$ is:
$$o_b=\frac{e^{z_b}}{\sum e^{z}}=\frac{e^{z_b}}{e^{z_a}+e^{z_b}+e^{z_c}} $$
and its derivation using the quotient rule:
$$\frac{\partial o_b} {\partial z_{b}}=\frac{e^{z_b}*\sum e^z - (e^{z_b})^2}{(\sum_j e^{z})^2}=\frac{e^{z_b}}{\sum e^z}-\frac{(e^{z_b})^2}{(\sum e^z)^2}$$
$$=softmax(b)-softmax^2(b)=o_b-o_b^2=o_b(1-o_b)$$
Back to the middle term for backpropagation this means:
$$\frac{\partial o_j} {\partial z_{j}}=o_j(1-o_j)$$
Putting it all together I get
$$\frac{\partial E} {\partial w_{ij}}= \frac{-t_j}{o_j}*o_j(1-o_j)*o_i=-t_j(1-o_j)*o_i$$
which means, if the target for this class is $t_j=0$, then I will not update the weights for this. That does not sound right.
Investigating on this I found people having two variants for the softmax derivation, one where $i=j$ and the other for $i\ne j$, like here or here.
But I can't make any sense out of this. Also I'm not even sure if this is the cause of my error, which is why I'm posting all of my calculations. I hope someone can clarify me where I am missing something or going wrong.
 A: In place of the $\{o_i\},\,$ I want a letter whose uppercase is visually distinct from its lowercase.  So let me substitute $\{y_i\}$. Also, let's use the variable $\{p_i\}$ to designate the $\{o_i\}$ from the previous layer. 
Let $Y$ be the diagonal matrix whose diagonal equals the vector $y$, i.e.
$$Y={\rm Diag}(y)$$
Using this new matrix variable and the Frobenius Inner Product we can calculate the gradient of $E$ wrt $W$.
$$\eqalign{
 z &= Wp+b &dz= dWp \cr
 y &= {\rm softmax}(z) &dy = (Y-yy^T)\,dz \cr
 E &= -t:\log(y) &dE = -t:Y^{-1}dy \cr\cr
dE &= -t:Y^{-1}(Y-yy^T)\,dz \cr
   &= -t:(I-1y^T)\,dz \cr
   &= -t:(I-1y^T)\,dW\,p \cr
   &= (y1^T-I)tp^T:dW \cr
   &= ((1^Tt)yp^T - tp^T):dW \cr\cr
\frac{\partial E}{\partial W} &= (1^Tt)yp^T - tp^T \cr
}$$
A: 
Note: I am not an expert on backprop, but now having read a bit, I think the following caveat is appropriate. When reading papers or books on neural nets, it is not uncommon for derivatives to be written using a mix of the standard summation/index notation, matrix notation, and multi-index notation (include a hybrid of the last two for tensor-tensor derivatives). Typically the intent is that this should be "understood from context", so you have to be careful!

I noticed a couple of inconsistencies in your derivation. I do not do neural networks really, so the following may be incorrect. However, here is how I would go about the problem.
First, you need to take account of the summation in $E$, and you cannot assume each term only depends on one weight. So taking the gradient of $E$ with respect to component $k$ of $z$, we have
$$E=-\sum_jt_j\log o_j\implies\frac{\partial E}{\partial z_k}=-\sum_jt_j\frac{\partial \log o_j}{\partial z_k}$$
Then, expressing $o_j$ as
$$o_j=\tfrac{1}{\Omega}e^{z_j} \,,\, \Omega=\sum_ie^{z_i} \implies \log o_j=z_j-\log\Omega$$
we have
$$\frac{\partial \log o_j}{\partial z_k}=\delta_{jk}-\frac{1}{\Omega}\frac{\partial\Omega}{\partial z_k}$$
where $\delta_{jk}$ is the Kronecker delta. Then the gradient of the softmax-denominator is
$$\frac{\partial\Omega}{\partial z_k}=\sum_ie^{z_i}\delta_{ik}=e^{z_k}$$
which gives
$$\frac{\partial \log o_j}{\partial z_k}=\delta_{jk}-o_k$$
or, expanding the log
$$\frac{\partial o_j}{\partial z_k}=o_j(\delta_{jk}-o_k)$$
Note that the derivative is with respect to $z_k$, an arbitrary component of $z$, which gives the $\delta_{jk}$ term ($=1$ only when $k=j$).
So the gradient of $E$ with respect to $z$ is then
$$\frac{\partial E}{\partial z_k}=\sum_jt_j(o_k-\delta_{jk})=o_k\left(\sum_jt_j\right)-t_k \implies \frac{\partial E}{\partial z_k}=o_k\tau-t_k$$
where  $\tau=\sum_jt_j$ is constant (for a given $t$ vector).
This shows a first difference from your result: the $t_k$ no longer multiplies $o_k$. Note that for the typical case where $t$ is "one-hot" we have $\tau=1$ (as noted in your first link).
A second inconsistency, if I understand correctly, is that the "$o$" that is input to $z$ seems unlikely to be the "$o$" that is output from the softmax. I would think that it makes more sense that this is actually "further back" in network architecture?
Calling this vector $y$, we then have
$$z_k=\sum_iw_{ik}y_i+b_k \implies \frac{\partial z_k}{\partial w_{pq}}=\sum_iy_i\frac{\partial w_{ik}}{\partial w_{pq}}=\sum_iy_i\delta_{ip}\delta_{kq}=\delta_{kq}y_p$$
Finally, to get the gradient of $E$ with respect to the weight-matrix $w$, we use the chain rule
$$\frac{\partial E}{\partial w_{pq}}=\sum_k\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}=\sum_k(o_k\tau-t_k)\delta_{kq}y_p=y_p(o_q\tau-t_q)$$
giving the final expression (assuming a one-hot $t$, i.e. $\tau=1$)
$$\frac{\partial E}{\partial w_{ij}}=y_i(o_j-t_j)$$
where $y$ is the input on the lowest level (of your example).
So this shows a second difference from your result: the "$o_i$" should presumably be from the level below $z$, which I call $y$, rather than the level above $z$ (which is $o$).
Hopefully this helps. Does this result seem more consistent?

Update: In response to a query from the OP in the comments, here is an expansion of the first step.
  First, note that the vector chain rule requires summations (see here). Second, to be certain of getting all gradient components, you should always introduce a new subscript letter for the component in the denominator of the partial derivative. So to fully write out the gradient with the full chain rule, we have
  $$\frac{\partial E}{\partial w_{pq}}=\sum_i \frac{\partial E}{\partial o_i}\frac{\partial o_i}{\partial w_{pq}}$$ 
  and
  $$\frac{\partial o_i}{\partial w_{pq}}=\sum_k \frac{\partial o_i}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}$$
  so
  $$\frac{\partial E}{\partial w_{pq}}=\sum_i \left[ \frac{\partial E}{\partial o_i}\left(\sum_k \frac{\partial o_i}{\partial z_k}\frac{\partial z_k}{\partial w_{pq}}\right) \right]$$
  In practice the full summations reduce, because you get a lot of $\delta_{ab}$ terms. Although it involves a lot of perhaps "extra" summations and subscripts, using the full chain rule will ensure you always get the correct result.

A: While @GeoMatt22's answer is correct, I personally found it very useful to reduce the problem to a toy example and draw a picture:

I then defined the operations each node was computing, treating the $h$'s and $w$'s as inputs to a "network" ($\mathbf{t}$ is a one-hot vector representing the class label of the data point):
$$L=-t_1\log o_1 -t_2\log o_2$$
$$o_1 = \frac{\exp(y_1)}{\exp(y_1) + \exp(y_2)}$$
$$o_2 = \frac{\exp(y_2)}{\exp(y_1) + \exp(y_2)}$$
$$y_1 = w_{11}h_1 + w_{21}h_2 + w_{31}h_3$$
$$y_2 = w_{12}h_1 + w_{22}h_2 + w_{32}h_3$$
Say I want to calculate the derivative of the loss with respect to $w_{21}$. I can just use my picture to trace back the path from the loss to the weight I'm interested in (removed the second column of $w$'s for clarity):

Then, I can just calculate the desired derivatives. Note that there are two paths through $y_1$ that lead to $w_{21}$, so I need to sum the derivatives that go through each of them.
$$\frac{\partial L}{\partial o_1} = -\frac{t_1}{o_1}$$
$$\frac{\partial L}{\partial o_2} = -\frac{t_2}{o_2}$$
$$\frac{\partial o_1}{\partial y_1} = \frac{\exp(y_1)}{\exp(y_1) + \exp(y_2)} - \left(\frac{\exp(y_1)}{\exp(y_1) + \exp(y_2)}\right)^2
= o_1(1 - o_1)$$
$$\frac{\partial o_2}{\partial y_1} = \frac{-\exp(y_2)\exp(y_1)}{(\exp(y_1) + \exp(y_2))^2} = -o_2o_1$$
$$\frac{\partial y_1}{\partial w_{21}} = h_2$$
Finally, putting the chain rule together:
\begin{align}
\frac{\partial L}{\partial w_{21}} &= \frac{\partial L}{\partial o_1}\frac{\partial o_1}{\partial y_1}\frac{\partial y_1}{\partial w_{21}} + \frac{\partial L}{\partial o_2}\frac{\partial o_2}{\partial y_1}\frac{\partial y_1}{\partial w_{21}}\\
&= \frac{-t_1}{o_1}[o_1(1 - o_1)]h_2 + \frac{-t_2}{o_2}(-o_2 o_1)h_2\\
&= h_2(t_2 o_1 - t_1 + t_1 o_1)\\
&= h_2(o_1(t_1 + t_2) - t_1)\\
&= h_2(o_1 - t_1)
\end{align}
Note that in the last step, $t_1 + t_2 = 1$ because the vector $\mathbf{t}$ is a one-hot vector.
A: The original question is answered by this post Derivative of Softmax Activation -Alijah Ahmed.
However writing this out for those who have come here for the general question of Backpropagation with Softmax and Cross-Entropy.
$$
\mathbf {
\bbox[10px, border:2px solid red] { \color{red}{
\begin{aligned}
 a^0 \rightarrow
    \bbox[5px, border:2px solid black]  {
      \underbrace{\text{hidden layers}}_{a^{l-2}} }
      \,\rightarrow
    \bbox[5px, border:2px solid black]  {  
      \underbrace{w^{l-1} a^{l-2}+b^{l-1}}_{z^{l-1} }
    }
      \,\rightarrow
    \bbox[5px, border:2px solid black]  {  
      \underbrace{\sigma(z^{l-1})}_{a^{l-1}}
    }
    \,\rightarrow
    \bbox[5px, border:2px solid black]  {  
     \underbrace{w^l a^{l-1}+b^l}_{z^{l}/logits }
    }
    \,\rightarrow
    \bbox[5px, border:2px solid black]  {  
    \underbrace{P(z^l)}_{\vec P/ \text{softmax} /a^{l}}
    }
    \,\rightarrow
  \bbox[5px, border:2px solid black]  {  
    \underbrace{L ( \vec P, \vec Y)}_{\text{CrossEntropyLoss}}
  }
\end{aligned}
}}}
$$
Derivative CrossEntropy Loss wrto Weight in last layer
$$
\mathbf {
\frac {\partial L}{\partial w^l} 
=  \color{red}{\frac {\partial L}{\partial z^l}}.\color{green}{\frac {\partial z^l}{\partial w^l}} \rightarrow \quad EqA1
}
$$
Where
$$
\mathbf {
L = -\sum_k y_k \log \color{red}{p_k} \,\,and \,p_j = \frac {e^ \color{red}{z_j}} {\sum_k e^{z_k}}
}
$$
Following from Derivative of Softmax Activation -Alijah Ahmed for the first term
$$ 
\color{red}
{
\begin{aligned}
    \frac {\partial L}{\partial z_i} = \frac {\partial ({-\sum_j y_k \log {p_k})}}{\partial z_i}
   \\ \\ \text {taking the summation outside} \\ \\
   = -\sum_j y_k\frac {\partial ({ \log {p_k})}}{\partial z_i}
  \\ \\ \color{black}{ \text {since } \frac{d}{dx} (f(g(x))) = f'(g(x))g'(x) }
  \\ \\
  = -\sum_k y_k * \frac {1}{p_k} *\frac {\partial { p_k}}{\partial z_i}
\end{aligned}
}
$$
The last term $\frac {\partial { p_k}}{\partial z_i}$ is the derivative  of Softmax wrto it's inputs also called logits. This is easy to derive and there are many sites that descirbe it. Example
Dertivative of SoftMax Antoni Parellada. The more rigorous derivative via the Jacobian matrix is here The Softmax function and its derivative-Eli Bendersky
$$
\color{red}
  {
  \begin{aligned}
   \frac {\partial { p_i}}{\partial z_i} = p_i(\delta_{ij} -p_j)
   \\ \\
   \delta_{ij} = 1 \text{ when i =j}
\\   \delta_{ij} = 0 \text{ when i} \ne \text{j}
  \end{aligned}
  }
$$
Using this above and repeating as is from Derivative of Softmax Activation -Alijah Ahmed   we get the below
$$ 
 \color{red}
{
\begin{aligned}
\frac {\partial L}{\partial z_i} = -\sum_k y_k * \frac {1}{p_k} *\frac {\partial { p_k}}{\partial z_i}
\\ \\
=-\sum_k y_k * \frac {1}{p_k} * p_i(\delta_{ij} -p_j)
\\ \\ \text{these i and j are dummy indices and we can rewrite this as}
\\ \\
=-\sum_k y_k * \frac {1}{p_k} * p_k(\delta_{ik} -p_i)
\\ \\ \text{taking the two cases and adding in above equation } \\ \\
\delta_{ij} = 1  \text{ when i =k} \text{ and }
\delta_{ij} = 0  \text{ when i} \ne  \text{k}
\\ \\
= [- \sum_i y_i * \frac {1}{p_i} * p_i(1 -p_i)]+[-\sum_{k \ne i} y_k * \frac {1}{p_k} * p_k(0 -p_i) ]
\\ \\
= [- y_i * \frac {1}{p_i} * p_i(1 -p_i)]+[-\sum_{k \ne i} y_k * \frac {1}{p_k} * p_k(0 -p_i) ]
\\ \\
= [- y_i(1 -p_i)]+[-\sum_{k \ne i} y_k *(0 -p_i) ]
\\ \\
= -y_i + y_i.p_i + \sum_{k \ne i} y_k.p_i
\\ \\
= -y_i + p_i( y_i + \sum_{k \ne i} y_k)
\\ \\
= -y_i + p_i( \sum_{k} y_k)
\\ \\
\text {note that } \sum_{k} y_k = 1 \, \text{as it is a One hot encoded Vector}
\\ \\
= p_i - y_i
\\ \\
\frac {\partial L}{\partial z^l} = p_i - y_i \rightarrow  \quad  \text{EqA.1.1}
\end{aligned}
}
$$
We now need to calculate the second term, to complete the equation
$$
\begin{aligned}
\frac {\partial L}{\partial w^l}
= \color{red}{\frac {\partial L}{\partial z^l}}.\color{green}{\frac {\partial z^l}{\partial w^l}}
\\ \\
\\
\color{green}{\frac {\partial z^l}{\partial w^l} = a^{l-1}} \text{ as } z^{l} = (w^l a^{l-1}+b^l)
\\ \\ \text{Putting all together} \\ \\
\frac {\partial L}{\partial w^l} = (p_i - y_i) *a^{l-1} \quad  \rightarrow  \quad  \mathbf {EqA1}
\end{aligned}
$$
Using Gradient descent we can keep adjusting the last layer like
$$
w{^l}{_i} = w{^l}{_i} -\alpha * \frac {\partial L}{\partial w^l}
$$
Now let's do the derivation for the inner layers, which is where the Chain Rule Magic happens
Derivative of Loss wrto Weight in Inner Layers
The trick here is to derivative the Loss wrto the inner layer as a composition of the partial derivative we computed earlier.
$$
\begin{aligned}
\frac {\partial L}{\partial w^{l-1}} 
=  \color{blue}{\frac {\partial L}{\partial z^{l-1}}}.
    \color{green}{\frac {\partial z^{l-1}}{\partial w^{l-1}}} \rightarrow \text{EqA.2}
\\ \\ 
\text{the trick is to represent the first  part  in terms of what we computed earlier; in terms of } \color{blue}{\frac {\partial L}{\partial z^{l}}}
\\  \\
\color{blue}{\frac {\partial L}{\partial z^{l-1}}} =
\color{blue}{\frac {\partial L}{\partial z^{l}}}.
    \frac {\partial z^{l}}{\partial a^{l-1}}.
    \frac {\partial a^{l-1}}{\partial z^{l-1}} \rightarrow \text{ EqMagic}
\\ \\
\color{blue}{\frac {\partial L}{\partial z^{l}}} = \color{blue}{(p_i- y_i)}
\text{ from the previous layer (from EqA1.1) } 
\\ \\
 z^l = w^l a^{l-1}+b^l
 \text{ which makes }
 {\frac {\partial z^{l} }{\partial a^{l-1}} = w^l} \text{ and }
 a^{l-1} = \sigma (z^{l-1}) 
  \text{ which makes }
\frac {\partial a^{l-1}}{\partial z^{l-1}} = \sigma \color{red}{'} (z^{l-1} )
\\ \\
 \text{ Putting together we get the first part of Eq A.2 }
\\ \\
\color{blue}{\frac {\partial L}{\partial z^{l-1}}} =\color{blue}{(p_i- y_i)}.w^l.\sigma \color{red}{'} (z^{l-1} ) \rightarrow \text{EqA.2.1 }
\\ \\
\text{Value of EqA.2.1 to be used in the next layer derivation in EqMagic)}
\\ \\
 z^{l-1} = w^{l-1} a^{l-2}+b^{l-1}
    \text{ which makes }
    \color{green}{\frac {\partial z^{l-1}}{\partial w^{l-1}}=a^{l-2}}
\\ \\
\frac {\partial L}{\partial w^{l-1}} 
=  \color{blue}{\frac {\partial L}{\partial z^{l-1}}}.
     \color{green}{\frac {\partial z^{l-1}}{\partial w^{l-1}}} = \color{blue}{(p_i- y_i)}.w^l.\sigma \color{red}{'} (z^{l-1} ).
\color{green}{a^{l-2}}
\end{aligned}
$$
Disclaimer
We see that with Chain Rule we can write out an expression that looks correct; and is correct in index notation. However when we implement with an actual case, with the above equation, your weights won't match out. This is due to the fact that we need to convert from index notation to  Matrix notation, and there some Matrix products have to be written out as Hadamard product $\odot$. Wihthout having some idea of these you cannot really understand this fully. A Primer on Index Notation John Crimaldi and The Matrix Calculus You Need For Deep Learning Terence,Jermy
A: Other answers have provided the correct way of calculating the derivative, but they do not point out where you have gone wrong. In fact, $t_j$ is always 1 in your last equation, cause you have assumed that $o_j$ takes that node of target 1 in your output; $o_j$ of other nodes have different forms of probability function, thus lead to different forms of derivative, so you should now understand why other people have treated $i=j$ and $i\neq j $ differently.
