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

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) 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 is not a Jacobian, but a vector of gradients, each computed at a different point. What am I missing?