What's the math for real world back-propagation? Considering a simple ANN:

$$x \rightarrow
f=(U_{m\times n}x^T)^T \rightarrow
g = g(f) \rightarrow
h = (V_{p \times m}g^T)^T \rightarrow
L = L(h,y)
$$

where $x\in\mathbb{R}^n$, $U$ and $V$ are matrices, $g$ is the point-wise sigmoid function, $L$ returns a real number indicating the loss by comparing the output $h$ with target $y$, and finally $\rightarrow$ represents data flow.

To minimize $L$ over $U$ and $V$ using gradient descent, we need to know $\frac{\partial L}{\partial U_{ij}}$ and $\frac{\partial L}{\partial V_{ij}}$, I know two ways to do this:


do the differentiation point wise, and having a hard time figuring out how to vectorize it
flatten $U$ and $V$ into a row vector, and use multivariate calculus (takes a vector, yields a vector) to do the differentiation


For the purpose of tutorial or illustration, the above two methods might be suffice, but say if you really want to implement back-prop by hand in the real world, what math will you use to do the derivative?
I mean, is there a branch, or method in meth, that teaches you how to take derivative of vector-valued function of matrices?
There is Matrix Calculus, (and I would recommend the very useful Matrix Cookbook as a bookmark to keep), but for the most part, when it comes to derivatives, it just boils down to pointwise differentiation and keeping your dimensionalities in check.

You might also want to look up Autodifferentiation. This is sort of a generalisation of the Chain Rule, such that it's possible to decompose any composite function, i.e. $a(x) = f(g(x))$, and calculate the gradient of the loss with respect to $g$ as a function of the gradient of the loss with respect to $f$.

This means that for every operation in your neural network, you can give it the gradient of the operation that "consumes" it, and it'll calculate its own gradient and propagate the error backwards (hence back-propagation)