# Application of chain rule, and some recursion

Consider the differentiable functions $$L^1(x,\theta^1),L^2(x^2,\theta^2),...,L^l(x^l,\theta^l)$$, where every $$x_k,\theta^k$$ are real vectors, for $$k=1,...,l$$. Also define $$\theta=(\theta^1,...,\theta^l)$$.

Define the composite function $$f(x,\theta)=x^{l+1}$$ recursively by doing $$x^k= L^{k-1}(x^{k-1},\theta^{k-1})$$, $$x^1=x$$.

Compute $$J_\theta f$$, the jacobian of $$f$$ with respect to $$\theta$$

For some context, I'm trying to implement gradient descent for optimizing the loss function of a neural network, and if my computations are correct I don't understand why we do back-propagation, instead of, say, forward-propagation... Here is my attempt, is there any mistake?

1. Compute $$J f$$: using the chain rule: $$Jf=JL^l(x^l,\theta^l)= \left ( J_{x^l}L^l\cdot J_{x,\theta^1,...,\theta^{l-1}}x^l \middle| J_{\theta^l}L^l\right )= \left ( J_{x^l}L^l\cdot J_{x,\theta^1,...,\theta^{l-1}}L^{l-1} \middle| J_{\theta^l}L^l\right )$$ Hence we can write $$Jf=J^l$$, where $$J^l$$ is given by the following recursive rule: $$J^k=\left ( J_{x^k}L^k\cdot J^{k-1}\middle| J_{\theta^k}L^k\right ), \quad J^1=J_{x,\theta^1}L^1$$

2. Obtain $$J_\theta f$$: we want to obtain the last columns of $$Jf$$, corresponding to the derivatives with respect to $$\theta^1,...,\theta^l$$. Clearly $$J_\theta f=\left ( J_{x^l}L^l\cdot J_{\theta^1,...,\theta^{l-1}}L^{l-1} \middle| J_{\theta^l}L^l\right )$$ Hence $$J_\theta f=G^l$$, where: $$G^k=\left ( J_{x^k}L^k\cdot G^{k-1}\middle| J_{\theta^k}L^k\right ), \quad G^1=J_{\theta^1}L^1$$

1. You wondered why backpropagation and not "forward-propagation". Khue gave a great answer, to which there is not much to add. As he said, automatic differentiation can be done in the forward mode or in the reverse mode. One way may require fewer arithmetic operations than the other, depending on the dimensions of the free parameters and the output. It is further explained in this answer.

As for the terminology, backpropagation stands for "backward propagation of errors", which is a name for the backward-mode differentiation in the context of neural networks. Calling a forward-mode differentiation a "forward-propagation" would be a bit inappropriate, since the error is the function's output and can be only propagated from that end.

2. Your derivations look correct to me. I am not sure whether you were simply asking for a verification or you were trying to derive the backpropagation in your own way, but got stuck. In the latter case, what you are missing is perhaps the right interpretation of your last line:

$$G^k=\left ( J_{x^k}L^k\cdot G^{k-1}\middle| J_{\theta^k}L^k\right ), \quad G^1=J_{\theta^1}L^1.\tag{1}\label{eq1}$$

This recursive relation indeed prompts us to start the computation with $$k=1,2,\dots$$, because $$G^1$$ is known and $$G^k$$ on the left-hand side depends on $$G^{k-1}$$ on the right-hand side; the computation is then straightforward.

This does not mean, however, that we can't start from the other end, $$k=l,l-1,\dots$$. Recall that we are interested not in $$G^k$$, but in the $$k$$-th columns of $$G^l$$. The last ($$l$$th) column of $$G^l$$ is readily available, as it does not depend on $$G^{l-1}$$:

$$G^l=\left ( J_{x^l}L^l\cdot G^{l-1}\middle| J_{\theta^l}L^l\right ).$$

For $$k=l-1$$ we need to take the second-to-last column. It does depend on $$G^{l-1}$$, but to be precise, it depends on the last column of $$G^{l-1}$$, which, in turn, does not depend on $$G^{l-2}$$. So we can pull it out, as follows:

$$G^{l}=\left(J_{x^{l}}L^{l}\cdot J_{x^{l-1}}L^{l-1}\cdot G^{l-2}|J_{x^{l}}L^{l}\cdot J_{\theta^{l-1}}L^{l-1}|J_{\theta^{l}}L^{l}\right),$$ which becomes $$G^{l}=\left(J_{x^{l-1}}L^{l}\cdot G^{l-2}|J_{\theta^{l-1}}L^{l}|J_{\theta^{l}}L^{l}\right).$$

At this point, it should be clear how to continue.

Update. In the above transition, the second-to-last column was computed as $$J_{\theta^{l-1}}L^{l}=J_{x^{l}}L^{l}\cdot J_{\theta^{l-1}}L^{l-1}$$. By analogy, we will observe that the consequent columns (moving from last to first) are computed as $$J_{\theta^{k-1}}L^{l}=J_{x^{k}}L^{l}\cdot J_{\theta^{k-1}}L^{k-1},\tag{2a}\label{eq3}$$

where $$J_{x^{k}}L^{l}$$ can be obtained through $$J_{x^{k}}L^{l}=J_{x^{k+1}}L^{l}\cdot J_{x^{k}}L^{k}.\tag{2b}\label{eq4}$$

The left-hand sides of \eqref{eq3}, \eqref{eq4} have $$k-1$$ and $$k$$, while the right-hand sides have $$k$$, $$k+1$$, and the terms we can know directly. So now you can use relations \eqref{eq3}, \eqref{eq4} recursively starting from $$k=l,l-1,\dots$$. This corresponds to the reverse-mode AD.

Of course, you could obtain \eqref{eq3}, \eqref{eq4} directly, without relying on your previous computations with $$G^k$$. I just wanted to show that where you stopped was not the dead end. If you were to start over, you would go like

Compute $$J_{\theta^{1}\dots\theta^{l}}f=\left(J_{\theta^{1}}f\mid\dots\mid J_{\theta^{l}}f\right)$$

where you'd carefully apply the chain rule for full derivatives in each column and would notice that the columns have common sub-expressions. I suppose, instead of going column by column you could formulate the same in a matrix form, like you did in \eqref{eq1}, but I don't see a point in such an exercise.

• +1 for the detailed explanation!
– f10w
Commented Oct 12, 2020 at 21:08
• I got a bit confused by your notation, since I use $J_zg$ if $g$ directly depends on the variable $z$, i.e. if I can write $g=g(z,...)$. I hope I have sorted it out. I have posted a summary-answer to check if I got the idea. Other than that, thank you very much for your time. Commented Oct 12, 2020 at 22:08
• To be less confusing, I tried as hard as I could to use your notation, which I personally would abandon if provided a chance. When $g$ directly depends on $z$ and we consider infinitesimal effects of this direct dependence, we call it a partial derivative $\partial g/\partial z$. When we are interested in all (direct and indirect) dependencies, we speak of total derivatives $\mathrm{d} g / \mathrm{d} z$. There is no way you can avoid total derivatives in this question: you yourself write $J_{x,\theta^1,...,\theta^{l-1}}L^{l-1}$, where $L^{l-1}$ is a function of only $x^{l-1}, \theta^{l-1}$. Commented Oct 12, 2020 at 23:55
• Your summary answer looks correct. [I cannot comment on it due to low reputation score, omg] Is it this the algorithm that is implemented in the backward pass of every layer? Conceptually, yes. Real implementations of automatic differentiation are based on computational graphs, which are built during the forward evaluation (at least, this is roughly how torch works). Computational graphs are flexible and can be used for many more things, compared to a hard-coded backpropagation in a multi-layer perceptron. Folks on Cross Validated or SO would be qualified to explain the implementation. Commented Oct 13, 2020 at 0:34
• @warm_fish You asked for a reference on backpropagation. While I never taught this subject and don't know which source is particularly good, I can suggest Chapter 6.5 of Deep Learning by Goodfellow, Bengio, and Courville deeplearningbook.org . This book is very well known in the machine learning community. Commented Oct 13, 2020 at 0:57

It's straightforward to see that the gradient of the output with respect to all the parameters can be computed in a recursive, forward manner (as you have shown above). This procedure is called the forward-mode differentiation. The well-known backpropagation algorithm, on the other hand, is a special case of the reverse-mode differentiation, which is much harder to see (that's why its invention is appreciated).

The question is, if the forward-mode differentiation is straightforward, why do people keep using the reverse mode?

The answer lies in the computational efficiency of the reverse mode. Indeed, for a general computational graph, if the dimension of the input is much larger than the output's, then the reverse mode is much more efficient (and vice-versa). This is a well-known result in automatic differentiation (see e.g. Who Invented the Reverse Mode of Differentiation? by Griewank).

It turns out that, in machine learning, the so-called training task often involves the gradient of a scalar-valued objective function with respect to a large number of parameters, i.e. the dimension of the output (1d) is much smaller than the dimension of the parameter vector (as well as the dimension of the input features), and thus the reverse-mode differentiation is much more efficient in this case.

(Try deriving the backpropagation algorithm yourself, then you'll see that the computation of the gradient of the loss will involve a lot of matrix-vector multiplications, which are much less expensive than the many matrix-matrix multiplications in the forward mode. I believe that you are able to see this yourself, but let me know if you need additional help.)

• Hey, thank you very much for your answer, and for your comment to the other question, which I bountied before receiving this answer. If you could point me towards a good source explaining backpropagation, it would be very helpful. So far, I've computed the gradient of the network (and derived the forward-mode differentiation), so I only need to derive the backward mode. Actually, I don't really have clear what backpropagation does, my current understanding is: use the chain rule on the $L^k$ above and then just multiply stuff together (in a backward fashion) Commented Oct 12, 2020 at 15:15
• @warm_fish Sorry I won’t be available for the next 10-12 hours. Please see the other great answer by paperskilltrees and let us know if you have further questions.
– f10w
Commented Oct 12, 2020 at 21:13

So, as far as I can understand, backwards differentiation is the following. After initializing $$D=I$$:

for $$k$$ from $$l$$ to $$1$$:

1. Save $$D\cdot J_{\theta^{k}}L^{k}$$ as $$J_{\theta^{k}}f$$
2. $$D=D\cdot J_{x^{k}}L^{k}$$

Is it this the algorithm that is implemented in the backward pass of every layer?