# Back-propagation through cross entropy or logistic loss function

I have neural network which ends with softmax function and I want to minimize cross-entropy cost function which takes output of this network and one-hot labels as arguments.

To calculate partial derivatives of cost function with respect to pre-softmax signal (output of fully connected layer before activation function), I may use chain rule of algebra and multiply derivatives of cost function with respect to network output with derivative of softmax function:

$$\begin{equation} \frac{\partial J(a, y)}{\partial z} = \frac{\partial J(a, y)}{\partial a} \circ \frac{\partial a(z)}{\partial z} = \frac{a - y}{a \circ (1 - a)} \circ a \circ (1 - a) = a - y \end{equation}$$

But if I want to “see” how gradient flows between loss function and softmax function, I often have a = 1 or a = 0 which leads to division by zero. Why it just works with a - y? Is there a way to see gradient between last layer activation and cost functions in computational graph? Is my calculations correct? What about other output activations / cost functions pairs, is there a way to separate them?