I’ve done some thinking and I think it should be differentiable wrt the mean…
My normal is simply the following:
def normal(self,mu,std,auxiliary,amplitude):
        mu = mu.view(-1)
        amplitude = amplitude.view(-1)
        a = 1/(std*torch.sqrt(2*torch.tensor(pi)))
        b = torch.exp(-0.5*(((auxiliary-mu[:,None])/std)**2))
        return amplitude[:,None]*(a*b)

Where auxiliary just refers to the 2D array with dimensions N by x^2 above.
Perhaps I should be using an inbuilt Pytorch normal distribution mechanism, or my handling of evaluating many normal distributions in one go (like in the code above) creates problems I can’t see… but as far as I can tell the differential of a normal distribution wrt the mean is possible. So I’m not sure what I’m doing wrong. 
In retrospect I now understand that the differential of a series of evaluated normal distributions with means, given as I in the post above, will simply return a zero so this is not a good way of tackling the problem.

That being said, would it not be possible to use an extension of autograd to simply backprop dL/dI = tensor-of-ones for the non-differentiable forward (in this case at the indexing step), or does that still prevent backprop functioning correctly for learning?