2. Backpropagation

I assume the reader is already familiar with backpropagation, but still, let's start by taking a step backwards and think about backpropagation in the large. The key principle underlying backpropagation is that a derivation is the arrow- part of a covariant functor from pointed differentiable manifolds to linear spaces. A.k.a, the chain rule. In the introduction I wrote that "backpropagation is a modular reformulation of the chain rule". It's time to elaborate.

Consider the two functions $g:R^{N_2}\rightarrow R^{N_1}$ and $f:R^{N_1}\rightarrow R$, and freely assume differentiability wherever it's needed: according to the chain rule, the differential of the composition $f\circ g:R^{N_2}\rightarrow R$ is given by $\nabla(f\circ g)=J_f\circ\nabla g$ (where $J_f$ is the Jacobian of $f$) - or in explicit coordinates, for all $0\le i\lt N_2$ we have $\frac{d(f\circ g)}{dx_i}=\sum_{j=0}^{N_1-1}\frac{df}{dx_j}\frac{dg}{dx_i}$.

In machine learning, typical learning algorithms involve a functional model $M:R^{N_I}\times R^{N_W}\rightarrow R^{N_O}$ and a loss-function $E:R^{N_O}\rightarrow R$, and are eventually reduced to an optimization of $E\circ M:R^{N_I}\rightarrow R$. The symbols $I,O,W$ and $E$ (for Inputs, Outputs, Weights and Error) are going to be used frequently from now on, and we treat $O:=M(I;W)$ as a function of $I\in R^{N_I}$ and $W\in R^{N_W}$, and $E(O)$ as a function of $O\in R^{N_O}$.

When the dimensions are high and the functions involved are sufficiently smooth, the best optimization strategies make use of the gradient $\frac{d(E\circ M)}{dW}$. Of course, non-smooth functions play a very important role in machine learning - but at least in principle, all of this works just as well with Lipschitz continuity and proximal-subgradients. By the chain rule, we can succinctly write: $\frac{dE}{dW}=\frac{dE}{dO}\frac{dO}{dW}$.

Now, if $M$ itself is a composition $M:=M_\text{outer}\circ M_\text{inner}$, then $\frac{dE}{dW_\text{inner}}=\frac{dE}{dO_\text{inner}}\frac{dO_\text{inner} }{dW_\text{inner}}$. But $O_\text{inner}=I_\text{outer}$, so the only contextual information $M_\text{inner}$ needs is $\frac{dE}{dI_\text{outer}}$. In its core, backpropagation is simply the recursive application of this idea. From a software-engineering perspective, it means it's possible to represent any (differentiable) real-valued composite function as a tree whose nodes are the functions' components, and for computing the evaluation of either the function or its gradient all the nodes can work locally. So each node can have many vector inputs $I_1,..., I_k$, and it should provide 2 algorithms (with a uniform interface across the nodes):

1. A function that computes (locally) $(I_1,..., I_k)\mapsto O$, referred to as the "Forward Algorithm".
2. A function that computes (locally) $(\frac{dE}{dO}, I_k)\mapsto\frac{dE}{dI_k}$ for each $k$, referred to the "Backward Algorithm".

There are several issues that deserve some more attention. First, nodes may be either parameterized or not. If they are, and those parameters are to be a subject of optimization, then the node should provide additional function that computes ($\frac{dE}{dO}, I)\mapsto\frac{dE}{dW}$. This is usually treated as part of the backward algorithm.

Secondly, the output of a node can be easily directed to several other nodes. In this case, its $\frac{dE}{dO}$ would be an accumulation of all the $\frac{dE}{dI}$s from its outgoing neighbours (since gradients are linear).

Thirdly, it is often very useful to express the backward algorithm in terms of the output $O$ instead of the inputs $I$. In practice, it's often rather simple to keep those values around after computing them in the forward-algorithm, and it can save a lot of work when executing the backward algorithm.

As an illustration, consider the very simple $R\mapsto R$ node whose forward algorithm computes $\omega x^2$, where $\omega$ is considered a parameter. Then the backward algorithm takes $\frac{dE}{dO}\in R$ as an input, and computes $\frac{dE}{dI}=2\omega\frac{dE}{dO}x=2\omega\frac{dE}{dO}I$ and $\frac{dE}{dW}=\frac{dE}{dO}x^2=\frac{dE}{dO}\omega^{-1}O$.

As for a more complicated example - well, the whole post is about one such example. Convolutional layers are functions that take an input $I$ and a kernel $W$, and in their forward algorithm computes the convolution $I\ast W$:

Next, the backward algorithms of convolutional layers will be discussed. Efficient algorithms and implementations for both the forward and backward algorithms for such functions are the subject of all of the following post.

Before going into details, a technical note: direct testing of an implementation for the backward algorithm is possible by comparing its results to a numerical differentiation of the forward algorithm. This procedure is known as gradient checking. Numerical differentiation is a pretty huge subject by its own right, but for the purpose of correctness tests, we can get by with a basic central finite-difference approximation. Python has an out-of-the-box implementation of it, given by scipy.optimize.check_grad. It works roughly as following:

An indirect sanity check for an implementation of the backward algorithm, which is certainly less reliable but more closely related to applications, is to use it for optimization. Python of course has a not-too-bad (not too-good either) optimization library as part of scipy, but I prefer using homemade toy SGD implementations for debugging. They can provide more insight on what's going on when things go south:

As for the loss function, we shall use MSE for debugging. The loss is given by $E(O,T)\propto ||O-T||$, and so its gradient is given by $\frac{dE}{dO}\propto 2(O-T)$:

Let's start with the inputs. From the chain rule, $\frac{dE}{dI}=\frac{dE}{dO}\frac{dO}{dI}$. The outputs are related to the inputs via linear convolution: the $n$-th output is $O_n=(I\ast W)_n:=\sum_{k=-M}^{+M}I_{n-k}W_{M+k}$. Thus $\frac{dO_n}{dI_k}=0$ whenever $|n-k|>M$, and otherwise, $\frac{dO_n}{dI_k}=W_{M+n-k}$. This implies that $\frac{dE}{dI_k}=\sum_{i=0}^{2M}\frac{dE}{dO_{k-M+i}}\frac{dO_{k-M+i}}{dI_k}=\sum_{i= 0}^{2M}\frac{dE}{dO_{k-M+i}}W_i$, and kinda gives an algorithm: each element of the gradient is computed via a dot-product:

As the term "kinda" gently hints, there a better perspective. Recall that $O_n$ can also be expressed as a dot-product: $O_n:=(I\ast W)_n=\langle I_{n-M}^{n+M}, W_\rho\rangle$ where $I_{n-M}^{n+M}$ is the projection of $I$ on $C^{2M+1}$ given by I[n-m:n+m+1] and $W_\rho$ is $W$ "reversed".

Looking back on what we have just done, we see that $\frac{dE}{dI}$ is actually a result of a cross-correlation, and can be computed via a convolution: $\frac{dE}{dI}=\frac{dE}{dO}\star W=\frac{dE}{dO}\ast W_\rho$. That's great! The formulation of the backwards algorithm in terms of the forward algorithm will allow us to reuse a single efficient implementation for both:

Can we do the same for the weights? Again, we start with the chain rule, $\frac{dE}{dW}=\frac{dE}{dO}\frac{dO}{dW}$, and consider $\frac{dO_n}{dW_k}$. This time we obtain $\frac{dO_n}{dW_k}=I_{n+M-k}$, thus: $\frac{dE}{dW_k}=\sum_{ i=0}^{N-1}\frac{dE}{dO_i}\frac{dO_i}{dW_k}=\sum_{i=0}^{N-1}\frac{dE}{dO_i}I_{i+M -k}$. Note that by applying a similar reasoning, then denoting by $\rho[I]$ a padded and reversed version of $I$, this means that $\frac{dE}{dW_k}=\sum_{i=0}^{N+2M}\frac{dE}{dO_i}\rho[I]_{N+k-i}$, so $\frac{dE}{dW}=\frac{dE}{dO}\star I=\frac{dE}{dO}\ast \rho[I]$. Success!

Note that this time, this is a convolution ("valid", not "full", in numpy-lingo) of two long sequences, and not a convolution of one long sequence with another short sequence as before:

To verify that we're on the right track, let's first test our functions with gradient-checking:

And finally, let's try to use those gradients in a toy-optimization problem. If our derivation is correct, we should able to make it work. First, minimizing the error with respect to the weights:

Then, minimizing the error with respect to the inputs:

A natural thought that might have occurred to you, is the following: when implementing a convolutional layer via a spectral algorithm, it seems redundant to constantly transform and then inversely transform the weights for each train sample (or a batch). It's reasonable to guess that by maintaining the weights in the frequency domain, and computing the gradient of the error with respect to the transformed weights, we can save some work at each iteration, and possibly meaningfully accelerate the whole training algorithm.

I may be wrong about it, but unfortunately, this doesn't seem to work. If you think otherwise, I'd be happy hear why and how. My reasoning is this: denoting $\hat{W}:=\mathcal{F}(W)$, we have $O_n=(I\ast W)_n=(I\ast \mathcal{F}^{-1}(\hat {W}))_n=\sum_{k=-M}^{+M}I_{n-k}\mathcal{F}^{-1}(\hat{W})_{M+k}$. Since by definition $\mathcal{F}^{-1}(\hat{W})_k=\frac{1}{2M+1}\sum_{n=0}^{2M}\hat{W}_ne^ {\frac{i2\pi kn}{2M+1}}$, we obtain $O_n=\frac{1}{2M+1}\sum_{k=-M}^{+M}I_{n-k}\sum_{r=0}^{2M}\hat{W}_re^{\frac{i2\pi (M+k)r}{2M+1}}$.

Denoting by $\rho X$ the reversal of $X$, this implies: $$\frac{dO_n}{d\hat{W}_r}=\frac{1}{2M+1}\sum_{k=-M}^{+M}I_{n-k}e^{\frac{i2\pi (M+k)r}{2M+1}}=\frac{1}{2M+1}\sum_{k=0}^{2M}I_{n-k+M}e^{\frac{i2\pi kr}{2M+1}}=\mathcal{F}^{-1}(\rho I_{n-M}^{n+M})_r$$ and we conclude: $\frac{dE}{d \hat{W}_r}=\sum_i\frac{dE}{dO_i}\frac{dO_i}{d\hat{W}_r}=\sum_i\frac{dE}{dO_i}\mathcal{F}^{-1}(\rho I_{i-M}^{i+M})_r$.

The result is kinda weird (note the application of an inverse Fourier transform for a spatial sequence), so let me assure myself that no silly mistakes - as I so often produce - were involved here, by a gradient-checking of a grotesquely unoptimized implementation of this algorithm:

So it seems correct, which is actually not happy news. Instead of 1 transformation-inversion pair of size $2M+1$ that takes $O(M\log{M})$, we now need N transformation-inversion pairs of that size, which takes $O(NM\log{M})$. Yuck. The implementation above actually does many more than that, but it is grotesquely unoptimized by design.

The root of all evil seems to be the fact the insistence of working directly with $\frac{dE}{dO}$, which is the fundamental input of a differential node. Apparently, expressing $\frac{dO_n}{d\mathcal{F}(W)_r}$ as a function of $\frac{dE}{dO}$ leads to a terribly inefficient algorithm. Unless, of course, there is some clever way I missed for quickly calculating the blob of weirdness $\sum_i\frac{dE}{dO_i}\mathcal{F}^{-1}(\sigma I_{i-M}^{i+M})_r$.

Alternatively, we can work "fully" in the frequency domain: since $\mathcal{F}(O)_n=\mathcal{F}(I)_n\cdot\mathcal{F}(W)_n$, we get that $\frac{d\mathcal{F}(O)_n}{d\mathcal{F}(W)_r}=0$ if $r\neq n$, and otherwise $\frac{d\mathcal{F}(O)_n}{d\mathcal{F}(W)_r}=\mathcal{F}(I)_n$. This implies that $\frac{dE}{\mathcal{F}(W)_r}=\sum_i\frac{dE}{d\mathcal{F}(O)_i}\frac{d\mathcal{F}(O)_i}{d\mathcal{F}(W)_r}=\frac{dE}{d\mathcal{F}(O)_r}\mathcal{F}(I)_r$.

So $\frac{dE}{d\mathcal{F}(O)_r}=\sum_n{(\frac{dE}{dO_n}\frac{dO_n}{d\mathcal{F} (O)_r})}=\frac{1}{N-2M}\sum_n{(\frac{dE}{dO_n}e^{\frac{i2\pi nr}{N-2M}})}=\mathcal{F}^{-1}(\frac{dE}{dO})_r$, and thus $\frac{dE}{\mathcal{F}(W)_r}=\mathcal{F}^{-1}(\frac{dE}{dO})_r\mathcal{F}(I)_r$. That's much better than the previous attempt, but not better than "regular" backpropagation: it is still $N\log{N}$, since we must calculate $\mathcal{F}^{-1}(\frac{dE}{dO})$.

There's nothing counter intuitive about it: $\frac{dE}{dW}$ is a function of $\frac{dE}{dO}$ and $I$. Even if the convolutional layer readily has $\mathcal{F}(I)$ (as an intermediate step of the forward algorithm), it still doesn't have $\mathcal{F}(\frac{dE}{dO})$ (if needed, it must be calculated from $\frac{dE}{dO}$). So either way, some DFT must be performed, which means the backwards algorithm for the weights will be always be superlinear with respect to $N$.

And I'm forced to conclude that the seemingly cool idea of backpropagation in the frequency domain is not helpful.