Solved by ptrblck in post #2 
              
              
                Your manual approach is not reducing the loss in the first call: 
loss = (Y_hat - Y) ** 2

so loss will have the shape [1, 1, 6, 7] and later you will sum this loss. 
However, nn.MSELoss uses the reduction='mean' by default, which would yield a lower loss value and thus also gradients with a smaller…