So, after a couple dozen tries I finally implemented a standalone nice and flashy softmax layer for my neural network in numpy. All works well, but I have a question regarding the maths part because there's just one tiny point I can't understand, like at all.
Having any kind of activation function in the output layer, back-propagation looks like:
$$ z = x\cdot w+b \\ a = f(z) \\ E = \frac{1}{N} \sum^{N}_{i=1} (a_{i} - y_{i})^{2} \\ \frac{\partial E}{\partial w} := \frac{\partial E}{\partial a_{i}} \cdot\frac{\partial a_{i}}{\partial z} \cdot \frac{\partial z}{\partial w} \\ \frac{\partial E}{\partial w} := (a_{i} - y_{i}) \cdot f^{'}(z) \cdot x $$
So by differentiating $ a_{l} $ with respect to $ z_{l} $, the result is the derivative of the activation function with $ z_{l} $ itself. Now, with Softmax in the final layer, this does not apply. If I use $ Softmax'(z_{l}) $ I get incorrect results, but I rather need $ Softmax'(a_{l}) $. Looking at a couple online materials like LINK in the Backpropagation phase section in the python code, the author also uses the activation itself in the argument for softmax's derivative and I've been differentiating there and back for the last 2 weeks, ... I don't understand why it works like that.
In short, my question is: With specifically Softmax in the output layer, why is my $ \frac{\partial a_{l}}{\partial z_{l}} $ becomes $ Softmax'(a_{l}) $ instead of $ Softmax'(z_{l}) $ like with pretty much any other normal activation function that's out there?