If I remember correctly computing the gradient of relu is a bit tricky precisely because of the problem your mentioning here in that the shapes don’t match. Let us say z is what’s computed before relu and let’s say the shape is 1x1000, and h is the result after passing it through relu. What happens when we compute the gradient of h w.r.t z (dh/dz) is that you have (in theory) to form the entire jacobian matrix.
This is because in theory every value from your input (1x1000) could have influenced every value in the resulting (1x1000) vector but of course since relu is applied element wise this is not actually the case. What happens then is that the only values would be along the diagonal of the jacobian and the rest would be zero (so some values of the diagonal would 1 if z > 0 and 0 otherwise) so it wouldn’t even be a diagonal matrix. Forming the jacobian is thus a huge waste of computation, and in the example given here would be a 1000x1000 matrix that we would need to matrix multiply with z. So if we form the jacobian we have 1x1000 mm 1000x1000 and then the shapes would match.
You can simplify this and do as he does in the tutorial by simply element wise multiplying with the diagonal values directly so in our case doing 1x1000 * 1x1000. This is equivalent to also doing grad_h[h < 0] = 0.