# Using Matrix Calculus in Backpropagation derivation. Rules of order of matmul and transposition when taking derivatives in different layouts.

I define a neural network with $$L$$ layers ($$L-1$$ hidden layers). The forward pass is as follows:

$$\mathbf{a}^{(l)} = f(\mathbf{W}^{(l)}\mathbf{a}^{(l-1)}+\mathbf{b}^{(l)})$$ Where $$l \in [0,L]$$ and $$\mathbf{a}^{(0)}=\mathbf{x}$$.

We have a dataset $$D = \{(\mathbf{X}_i, \mathbf{Y}_i) \mid \mathbf{X}_i = \mathbf{x} \in \mathbf{X}, \mathbf{Y}_i = \mathbf{y} \in \mathbf{Y}, i = 1, \ldots, N\}$$, where $$N$$ is the size of the dataset and $$\mathbf{x}$$ and $$\mathbf{y}$$ are the $$i$$th input and label of the dataset.

We have a cost/loss function that for example is defined as half the mean square error. $$N_L$$ is the size of the output layer. For a single data point (SGD) we have: $$C(\mathbf{a}^{(L)},\mathbf{y}) = \frac{1}{N_{L}}\sum_{i=1}^{N_{L}}{\frac{1}{2}(a^{(L)}_{i}-y_i)^2}$$

We want to minimise the cost function to train the network.

$$\min_{\mathbf{W^{(l)}}, \mathbf{b^{(l)}} \space \forall{l \in [0,L]}} C(\mathbf{a}^{(L)},\mathbf{y})$$

We derive backpropagation algorithm without matrix calculus as follows (where$$\space\mathbf{a}^{(l)}=f(\mathbf{z}^{(l)})\space$$).

For the final layer ($$l=L$$): $$\frac{\partial C}{\partial W^{(L)}_{ij}} = \frac{\partial C}{\partial a^{(L)}_i}\frac{\partial a^{(L)}_i}{\partial z^{(L)}_i}\frac{\partial z^{(L)}_i}{\partial W^{(L)}_{ij}} = (a^{(L)}_i - y_i)f'(z_i)a^{(L-1)}_j$$ For the penultimate layer ($$l=L-1$$) using $$j$$ and $$k$$ as indexes to the $$L-1$$ and $$L-2$$ neuron layers: $$\frac{\partial C}{\partial W^{(L-1)}_{ij}} = (\sum_{i}{\frac{\partial C}{\partial a^{(L)}_i}\frac{\partial a^{(L)}_i}{\partial z^{(L)}_i}\frac{\partial z^{(L)}_i}{\partial a^{(L-1)}_{j}}}) \frac{\partial a^{(L-1)}_{j}}{\partial z^{(L-1)}_j}\frac{\partial z^{(L-1)}_{j}}{\partial W^{(L-1)}_{jk}}$$

Now for a general layer $$l$$, where $$k$$, $$i$$ and $$j$$ are the indexes for the next, current and previous layers respectively. We now generalise the summation term in the parentheses above ($$l=L-1$$) for any $$l \in [0,L-1]$$: $$\frac{\partial C}{\partial a^{(l)}_i} = \sum_{k}{\frac{\partial C}{\partial a^{(l+1)}_k}\frac{\partial a^{(l+1)}_k}{\partial z^{(l+1)}_k}\frac{\partial z^{(l+1)}_k}{\partial a^{(l)}_{i}}} = \sum_{k}{\frac{\partial C}{\partial a^{(l+1)}_k} f(z^{(l+1)}_k) W^{(l)}_{ki}}$$ And for any $$l$$: $$\frac{\partial C}{\partial a^{(l)}_i} = \begin{cases} \sum_{k}{\frac{\partial C}{\partial a^{(l+1)}_k}\frac{\partial a^{(l+1)}_k}{\partial z^{(l+1)}_k}\frac{\partial z^{(l+1)}_k}{\partial a^{(l)}_{i}}} & \text{if} \space l < L \\ \frac{\partial C}{\partial a^{(l)}_i} & \text{if} \space l = L \end{cases}$$

We now define $$\delta$$:

$$\delta^{(l)}_i = \frac{\partial C}{\partial a^{(l)}_i} \frac{\partial a^{(l)}_i}{\partial z^{(l)}_i}$$

For any $$l$$ we have the gradient for the weights (and a similar derivtion for biases): $$\frac{\partial C}{\partial W^{(l)}_{ij}} = \delta^{(l)}_i a^{(l)}_j \space , \space \frac{\partial C}{\partial b^{(l)}_{i}} = \delta^{(l)}_i$$

I now derive backpropagation using matrix calculus.

In matrix calculus there are denominator and numerator layouts.

Denominator: $$\frac{\partial f}{\partial \mathbf{x}} = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_n} \end{bmatrix}$$ Numerator: $$\frac{\partial f}{\partial \mathbf{x}} = \begin{bmatrix} \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \dots, \frac{\partial f}{\partial x_n}\end{bmatrix}$$

Using denominator layout:

$$\frac{\partial C}{\partial \mathbf{W}^{(L)}} = \frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial C}{\partial \mathbf{a}^{(L)}}\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{W}^{(L)}} = \mathrm{diag}(f'(\mathbf{z}))(\mathbf{a}^{(L)} - \mathbf{y})(\mathbf{a}^{(L-1)})^T$$

$$\frac{\partial C}{\partial \mathbf{W}^{(L-1)}} = \frac{\partial \mathbf{a}^{(L-1)}}{\partial \mathbf{z}^{(L-1)}}(\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{a}^{(L-1)}}\frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial C}{\partial \mathbf{a}^{(L)}}) \frac{\partial \mathbf{z}^{(L-1)}}{\partial \mathbf{W}^{(L-1)}} = (\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{a}^{(L-1)}}\frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial C}{\partial \mathbf{a}^{(L)}}) \frac{\partial \mathbf{z}^{(L-1)}}{\partial \mathbf{W}^{(L-1)}}\frac{\partial \mathbf{a}^{(L-1)}}{\partial \mathbf{z}^{(L-1)}}$$ We now generalise the summation term in the parentheses above ($$l=L-1$$) for any $$l \in [0,L-1]$$ $$\frac{\partial C}{\partial \mathbf{a}^{(l)}} = \frac{\partial \mathbf{z}^{(l+1)}}{\partial \mathbf{a}^{(l)}}\frac{\partial \mathbf{a}^{(l+1)}}{\partial \mathbf{z}^{(l+1)}}\frac{\partial C}{\partial \mathbf{a}^{(l+1)}} = (\mathbf{W}^{(l)})^T \mathrm{diag}(f'(\mathbf{z}^{(l+1)})) \frac{\partial C}{\partial \mathbf{a}^{(l+1)}}$$ And for any $$l$$: $$\frac{\partial C}{\partial \mathbf{a}^{(l)}} = \begin{cases} \frac{\partial \mathbf{z}^{(l+1)}}{\partial \mathbf{a}^{(l)}}\frac{\partial \mathbf{a}^{(l+1)}}{\partial \mathbf{z}^{(l+1)}}\frac{\partial C}{\partial \mathbf{a}^{(l+1)}} & \text{if} \space l < L \\ \frac{\partial C}{\partial \mathbf{a^{(l)}}} & \text{if} \space l = L \end{cases}$$

We now define $$\delta$$:

$$\mathbf{\delta}^{(l)} = \frac{\partial \mathbf{a}^{(l)}}{\partial \mathbf{z}^{(l)}}\frac{\partial C}{\partial \mathbf{a}^{(l)}}$$ For any $$l$$ we have the gradient for the weights (and a similar derivtion for biases): $$\frac{\partial C}{\partial \mathbf{W}^{(l)}} = \mathbf{\delta}^{(l)} (\mathbf{a}^{(l)})^T \space , \space \frac{\partial C}{\partial \mathbf{b}^{(l)}} = \mathbf{\delta}^{(l)}$$

Using numerator layout:

$$\frac{\partial C}{\partial \mathbf{W}^{(L)}} = \frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{W}^{(L)}}\frac{\partial C}{\partial \mathbf{a}^{(L)}} = \mathrm{diag}(f'(\mathbf{z}))\space\mathbf{a}^{(L-1)}(\mathbf{a}^{(L)} - \mathbf{y})^T$$

$$\frac{\partial C}{\partial \mathbf{W}^{(L-1)}} = \frac{\partial \mathbf{z}^{(L-1)}}{\partial \mathbf{W}^{(L-1)}}(\frac{\partial C}{\partial \mathbf{a}^{(L)}}\frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{a}^{(L-1)}})\frac{\partial \mathbf{a}^{(L-1)}}{\partial \mathbf{z}^{(L-1)}} = (\frac{\partial C}{\partial \mathbf{a}^{(L)}}\frac{\partial \mathbf{a}^{(L)}}{\partial \mathbf{z}^{(L)}}\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{a}^{(L-1)}})\frac{\partial \mathbf{a}^{(L-1)}}{\partial \mathbf{z}^{(L-1)}}\frac{\partial \mathbf{z}^{(L-1)}}{\partial \mathbf{W}^{(L-1)}}$$ We now generalise the summation term in the parentheses above ($$l=L-1$$) for any $$l \in [0,L-1]$$ $$\frac{\partial C}{\partial \mathbf{a}^{(l)}} = \frac{\partial C}{\partial \mathbf{a}^{(l+1)}}\frac{\partial \mathbf{a}^{(l+1)}}{\partial \mathbf{z}^{(l+1)}}\frac{\partial \mathbf{z}^{(l+1)}}{\partial \mathbf{a}^{(l)}} = \frac{\partial C}{\partial \mathbf{a}^{(l+1)}} \mathrm{diag}(f'(\mathbf{z}^{(l+1)})) \mathbf{W}^{(l)}$$ And for any $$l$$: $$\frac{\partial C}{\partial \mathbf{a}^{(l)}} = \begin{cases} \frac{\partial C}{\partial \mathbf{a}^{(l+1)}}\frac{\partial \mathbf{a}^{(l+1)}}{\partial \mathbf{z}^{(l+1)}}\frac{\partial \mathbf{z}^{(l+1)}}{\partial \mathbf{a}^{(l)}} & \text{if} \space l < L \\ \frac{\partial C}{\partial \mathbf{a^{(l)}}} & \text{if} \space l = L \end{cases}$$

We now define $$\delta$$ (a row vector):

$$\mathbf{\delta}^{(l)} = \frac{\partial C}{\partial \mathbf{a}^{(l)}}\frac{\partial \mathbf{a}^{(l)}}{\partial \mathbf{z}^{(l)}}$$ For any $$l$$ we have the gradient for the weights (and a similar derivtion for biases): $$\frac{\partial C}{\partial \mathbf{W}^{(l)}} = \mathbf{a}^{(l)} \mathbf{\delta}^{(l)} \space , \space \frac{\partial C}{\partial \mathbf{b}^{(l)}} = \mathbf{\delta}^{(l)}$$

As you can see when I do chain rule for matrix calculus in denominator and numerator layouts I find multiple ways I can order multiplications. Are there rules for ordering the chain rule multiplications in matrix calculus? Is there a situation agnostic way/rule of doing chain rule in matrix calculus which you can follow to order the multiplications and transpositions without having to do it in index form?

• This question is a bit verbose, do you mind tagging the equations which you are comparing when you say "As you can see when I do chain rule for matrix calculus in denominator and numerator layouts I find multiple ways I can order multiplications." Commented Feb 21 at 7:30

There is in fact a chain-rule for matrix calculus: https://ccrma.stanford.edu/~dattorro/matrixcalc.pdf

However, from an implementation perspective, there is no "right" answer, it depends mostly on how the deep learning framework is optimized for the underlying hardware. The idea is to maximize the amount of parallelism achieved, because ML training is done on GPUs (which execute in parallel).

In frameworks like PyTorch/Tensorflow, the dominating paradigm is autograd, which constructs an execution graph in a very granular way (at the level of individual calculations like summation, subtraction, etc) and stores gradients for each granular calculation. Here is an overview: https://pytorch.org/blog/overview-of-pytorch-autograd-engine/. The execution thus proceeds layer-by-layer, with each layer running the matmuls in paralle.

You can probably find the exact ordering but it requires sifting through a lot of CUDA code. This detail is not important in terms of understanding neural nets, unless you are proposing a change to the underlying PyTorch execution graph or something similar.

I don't think it matters because on back-propagation you would need to update the weights. meaning the dimensions of the derivatives need to match the dimensions. So, practically, I think the col/row representation of the derivatives is fixed.

Because instead of taking the derivative form and you write out exactly what the derivative values are you will be constrained by the order of operations of the forward propagation.

• Yeah but im asking about specific ordering of operations when it comes to matrix calculus applied to backpropagation. Are there ordering of operations that look like the normal single variable chain rule (or a flipped version).
– xTom
Commented Feb 20 at 22:13