I want to start by taking an example for a normal neural network with 2 input nodes, 3 hidden nodes and 2 output nodes. Let the weights between input and hidden nodes are $W_i{_j}$ (2x3) and weights between hidden nodes and output (3x2) are $W_j{_k}$. Last layer has linear activation. The gradient of error wrt:
$$W_i{_j}= W_j{_k} * g'(W_i{_j} * x) * x * err$$
But this does not satisfy the weight update equation as the product gives the size 2x1. But it should be 2x3.
It's important to note that you are storing weights towards previous layer in columns. So my example will be specifically for that case.
If you are using a different notation (weights from neuron towards previous layer are in rows of your matrix) - swap the words "rows" and "columns" below:
For your case, any incoming gradient (from the layers above or from the Error aka Cost function) must have dimension 2. It must be a vector. Otherwise it will not physically "match" your output nodes.
Section A
The gradient wrt to your weights (hidden-to-output layer) will be a 3x2 matrix. You can subtract this matrix from your weights, to correct them. Don't forget about learning rate.
To build this matrix, copy the incoming gradient (vec of dimension 2, already passed through the activation function) into every row of a new, empty matrix (3x2), then component-wise multiply its every column by the fwdprop output of your 3-node hidden layer. This is your GradForWeights matrix
Section B
The gradient wrt hidden nodes will be a 3d vector.
To make it, use dot product between your Weights (3x2) and the incoming gradient vector (2d). The resulting vector will be of dimension 3, and can also be used as the "incoming gradient vector" for the layer below (of course, after you pass it through the activation function of that sub-layer).
