Suppose a multi-layer feed-forward neural network, e.g.:

[![Network Structure][1]][1] 

Using matrix form to account for all training samples $(i)$, the forward propagation can be written as follows:

$Z^{[l]}=W^{[l]}A^{[l-1]}+\overline{b}^{[l]}$

$A^{[l]}=g^{[l]}(Z^{[l]})$

where $g^{[l]}$ is the activation function used at layer ${[l]}$. 


Let $L$ denote the loss function. For the backpropagation, we want to compute partial derivatives of $L$ with respect $z^{[l](i)}_j$ for all nodes $j$ of the layer $[l]$ and all training examples $(i)$. Many tutorials (e.g. [this][2]) call the resulting matrix a Jacobian. I do not understand how this is the case.

In particular, we can view $L$ as a function of the inputs at the layer $[l]$, i.e. $L=L(z^{[l]}_1, z^{[l]}_2,\ldots,z^{[l]}_{n^{[l]}})$. For each training sample, the output of this function is a **scalar**, whereas the definition of Jacobian requires that the function's output be a **vector**. So, it seems to me that what we have here (i.e. when we join the derivatives for all the training samples into one matrix form) is not a Jacobian, but a vector of gradients, each computed at a different point. What am I missing?

  [1]: https://i.sstatic.net/FeNRC.png
  [2]: https://www.jasonosajima.com/backprop