Attempting to implement a vectorized version of one of the backpropagation equations I'm trying to implement a neural net that can perform backpropagation on a whole minibatch rather than iterating through each element.
Here is what I want to do:
Given two matrices

*

*$\delta_l$ of shape (batch size, # of neurons in in layer $l$). Where each element $(k, i)$ corresponds to the error for the neuron $i$ in layer $l$ at minibatch $k$. Each element looks like this $\dfrac{\partial C}{\partial z^l_{i}}$.


*$a_{l-1}$ of shape (batch size, # of neurons in layer $l-1$). Where the each element $(k, j)$ corresponds to the activation of neuron $j$ in layer $l-1$ at minibatch $k$. Each element $a^{l-1}_j$ is equal to (calculations not shown) $\dfrac{\partial z^l_i}{\partial w^l_{ij}}$ for some arbitrary neuron $i$ in layer $l$. Note that the notation I am using for weights is $w^{\text{layer}}_{\text{where it is going, where it is coming from}}$.
The reason these two matrices are used is because using the chain rule, we can see that for any weight $w^l_{ij}$, we can compute it's derivative with respect to the cost by $ \dfrac{\partial C}{\partial z^l_{i}} \dfrac{\partial z^l_i}{\partial w^l_{ij}}$
Let $m$ be the number of neurons in layer $l$ and $n$ be the number of neurons in layer $l-1$. Using these two matrices, I want to compute a matrix that represents the gradient of the weights, which looks like the following:
$$\nabla W = \begin{bmatrix}
\dfrac{\partial C}{\partial z^l_1} a^{l-1}_1 & \cdots & \dfrac{\partial C}{\partial z^l_1} a^{l-1}_n \\
\vdots &\ddots & \vdots \\ 
\dfrac{\partial C}{\partial z^l_m} a^{l-1}_1 & \cdots & \dfrac{\partial C}{\partial z^l_m} a^{l-1}_n \\
\end{bmatrix}$$
$$  = \begin{bmatrix}
\dfrac{\partial C}{\partial z^l_1} \dfrac{\partial z^l_1}{\partial w_{1,1}^l} & \cdots &\dfrac{\partial C}{\partial z^l_1} \dfrac{\partial z^l_1}{\partial w_{1,n}^l} \\
\vdots & \ddots & \vdots \\
\dfrac{\partial C}{\partial z^l_m} \dfrac{\partial z^l_m}{\partial w_{m,1}^l}& \cdots & \dfrac{\partial C}{\partial z^l_m} \dfrac{\partial z^l_m}{\partial w_{m,n}^l}\\
\end{bmatrix}$$
Is there way I can use these two matrices and compute the gradient of the weights that connect layers $l$ and $l-1$ together? I can think of a way to do it with NumPy operations (Using tiling and transposing), but I can't think of simpler, more elegant way to do it. Perhaps I'm just looking at this whole thing wrong and should go around representing my error and activations in another way.
 A: Yes, by using those two matrices, $\delta_l$ and $a_{l-1}$ , it's possible to do a vectorized version of the backpropagation equation that gives the updates of the weights, $\Delta W$, using a mini-batch of size $K$.
Let's see why is this the case, and also just for clarity the notation that is going to be used along this post.
Notation
For the sake of clarity, the main notation that I'm going to use (using also the notation of the question) is:

*

*$\Delta\to$ to express updates

*$K\to$ size of the mini-batch

*$m\to$ number of neurons at layer $l$

*$n\to$ number of neurons at layer $l-1$

*$\text{mb}\to$ Subscript to denote that a variable contains information of the whole minibatch.

Using this notation, I've taken the liberty of changing a bit the original notation given by the question $\to$ the updates of the weights that connect a layer $l-1$ to a layer $l$ that we want to achieve in a vectorized way, will be contained in $\Delta W^l_{\text{mb}}$.
Besides, by using this notation, the matrix of errors, $\delta$, at layer $l$ and the matrix of activations, $a$, at layer $l-1$ will be expressed as $\delta_{\text{mb}}^l$ and $a_{\text{mb}}^{l-1}$ (instead of $\delta_l$ and $a_{l-1}$).
As said earlier, by using the matrices $\delta_{\text{mb}}^l$ and $a_{\text{mb}}^{l-1}$, we can express $\Delta W^l_{\text{mb}}$ by a quantity proportional to:
$$ \underbrace{\Delta W^l_{\text{mb}}}_{m\times n} \propto \underbrace{(\delta_{\text{mb}}^l)^T}_{m\times K} \,\,\,\underbrace{a_{\text{mb}}^{l-1}}_{K \times n}$$
Which represent a vectorized implementation.
Let's see why this works.
What is contained in $\Delta W^l_{\text{mb}}$
As said in the question, the weight updates for a single sample (not a mini-batch), $\Delta W^l$, are given by:
$$\Delta W^l \propto \frac{\partial C}{\partial W^l} = 
\pmatrix{\frac{\partial C}{\partial w^l_{11}} & \frac{\partial C}{\partial w^l_{12}}  & \cdots & \frac{\partial C}{\partial w^l_{1n}} \\
\frac{\partial C}{\partial w^l_{21}} & \frac{\partial C}{\partial w^l_{22}}  & \cdots & \frac{\partial C}{\partial w^l_{2n}} \\
\vdots & \vdots & \ddots & \vdots\\
\frac{\partial C}{\partial w^l_{m1}} & \frac{\partial C}{\partial w^l_{m2}}  & \cdots & \frac{\partial C}{\partial w^l_{mn}}}$$
Because of this, what is contained in $\Delta W^l_{\text{mb}}$ will be given by:
$$\Delta W^l_{\text{mb}} \propto \sum_{k=1}^K \left(\frac{\partial C}{\partial W^l}\right)_k = 
\sum_{k=1}^K \pmatrix{\frac{\partial C}{\partial w^l_{11}} & \frac{\partial C}{\partial w^l_{12}}  & \cdots & \frac{\partial C}{\partial w^l_{1n}} \\
\frac{\partial C}{\partial w^l_{21}} & \frac{\partial C}{\partial w^l_{22}}  & \cdots & \frac{\partial C}{\partial w^l_{2J}} \\
\vdots & \vdots & \ddots & \vdots\\
\frac{\partial C}{\partial w^l_{m1}} & \frac{\partial C}{\partial w^l_{I2}}  & \cdots & \frac{\partial C}{\partial w^l_{mn}}}_k$$
Which is the same as:
$$\Delta W^l_{\text{mb}} \propto 
\pmatrix{\sum_k^K (\frac{\partial C}{\partial w^l_{11}})_k & \sum_k^K (\frac{\partial C}{\partial w^l_{12}})_k  & \cdots & \sum_k^K (\frac{\partial C}{\partial w^l_{1n}})_k \\
\sum_k^K (\frac{\partial C}{\partial w^l_{21}})_k & \sum_k^K (\frac{\partial C}{\partial w^l_{22}})_k  & \cdots & \sum_k^K (\frac{\partial C}{\partial w^l_{2n}})_k \\
\vdots & \vdots & \ddots & \vdots\\
\sum_k^K (\frac{\partial C}{\partial w^l_{m1}})_k & \sum_k^K (\frac{\partial C}{\partial w^l_{m2}})_k  & \cdots & \sum_k^K (\frac{\partial C}{\partial w^l_{mn}})_k}$$

Achieving $\sum_k^K (\frac{\partial C}{\partial w^l_{ij}})_k$
To compute each element of the previous matrix, $\sum_k^K (\frac{\partial C}{\partial w^l_{ij}})_k$, we can make use of the columns of $\delta_{\text{mb}}^l$ and $a_{\text{mb}}^{l-1}$:

*

*$\delta^l_{\text{mb}} (:,i) = (\delta^l_{\text{mb}} (1,i), \delta^l_{\text{mb}} (2,i), \cdots, \delta^l_{\text{mb}} (K,i))^T \Rightarrow$ It has $K\times1$ dimensions.

*$a_{\text{mb}}^{l-1} (:,j) = (a_{\text{mb}}^{l-1} (1,j), a_{\text{mb}}^{l-1} (2,j), \cdots, a_{\text{mb}}^{l-1} (K,j))^T \Rightarrow$ It also has $K\times1$ dimensions.

This way, each element will be given by $ \Rightarrow \sum_k^K (\frac{\partial C}{\partial w^l_{ij}})_k = (\delta^l_{\text{mb}} (:,i))^T \,\,a_{\text{mb}}^{l-1} (:,j)$

Achieving $\Delta W^l_{\text{mb}}$
To extend the above reasoning of a certain element to the whole matrix itself, we have to make use of mini-batch matrices mentioned in the question:
$$\begin{aligned}
(\delta_{\text{mb}}^l)^T &= \pmatrix{
\delta_{1,1}^l & \delta_{1,2}^l & \cdots & \delta_{1,K}^l\\
\delta_{2,1}^l & \delta_{2,2}^l & \cdots & \delta_{2,K}^l\\
\vdots & \vdots & \ddots & \vdots\\
\delta_{m,1}^l & \delta_{m,2}^l & \cdots & \delta_{m,K}^l
} &&  \text{Matrix of mini-batch errors}\\
\\
a_{\text{mb}}^{l-1} &= \pmatrix{
a_{1,1}^{l-1} & a_{2,1}^{l-1} & \cdots & a_{n,1}^{l-1}\\
a_{1,2}^{l-1} & a_{2,2}^{l-1} & \cdots & a_{n,2}^{l-1}\\
\vdots & \vdots & \ddots & \vdots\\
a_{1,K}^{l-1} & a_{2,K}^{l-1} & \cdots & a_{n,K}^{l-1}\\
} && \text{Matrix of mini-batch activations}
\end{aligned}$$
This way, we have proven that we can achieve a vectorized way of obtanining the updates of the weights, because:
$$\underbrace{\Delta W^l_{\text{mb}}}_{m\times n} \propto \underbrace{(\delta_{\text{mb}}^l)^T}_{m\times K} \,\,\,\underbrace{a_{\text{mb}}^{l-1}}_{K \times n}$$
