Does the popular packages like PyTorch, Tensorflow, Keras, etc., use this or a variant of this algorithm to find the gradients at a particular point?
Yes. This is effectively what back-propagation is. However, there are a couple of important details:
Using a loss function flattens the matrix form you have to a vector, because with a loss function the dimension of the output $n=1$. This is important, since it is not possible to directly search for a maximum or minimum value for multi-dimensional output.
Back-propagation is resolved as follows:
- Set up by constructing a computation graphcomputation graph that composes multiple functions. Some libraries set this graph up directly e.g. TensorFlow, whilst others like PyTorch can resolve it dynamically at run time from requested computations (this is convenient and easier to write, but there is a performance cost).
- Gradients are calculated numerically by resolving function composition in reverse using the chain rule.
- There is an implied analytic form for the forward calculation, but it never used directly.
If yes, will those packages be able to write the analytical forms of all the output variables in terms of input variables?
In theory, yes, as the data to do so is in the computational graphs that are used to run neural networks forward, and is the same data that is used to resolve back propagation. In practice I have not seen this done. There may be an add-on or library that could print out the analytic forms of neural network outputs at any layer, or for the loss function. If there is not one, it would not be too hard to write one.
However, other than as a teaching aid for very small networks, then the analytical view of the outputs or loss function are not a practical or useful description of what is going on. The graph view showing how the linear algebra and non-linear activations are composed into a sequence of computations (and reversed step-by-step for gradient calculations) is much easier to comprehend. Many neural network libraries can generate views of this compuational graph.
One important take-away though, is that despite the apparent complexity of back-propagation, it is doing exactly what you expect in the question: Calculating the gradient components in a Jacobian.