Positive bias causes the calculation of incorrect gradients
Ok, so I was right in assuming it was the discontinuity at x=0 for relu(x). The positive bias actually had nothing to do with the incorrect gradients. As for why in my edit I was getting such massive differences between gradient values, it was because I initialised the input to be integers, and upon adding a decimal, no matter how small (> 0 though), numpy would jump entire integers, so -2 + 0.0001 = -1.

This is reproduced by:

&gt;&gt;&gt; a = np.random.randint(-3,3,(1,3,3))
&gt;&gt;&gt; a
array([[[-2,  1, -2],
        [-2, -1, -3],
        [-3,  0,  1]]])
&gt;&gt;&gt; a[0][0][0] += 0.0001
&gt;&gt;&gt; a
array([[[-1,  1, -2],
        [-2, -1, -3],
        [-3,  0,  1]]])


The problem with my input is large amounts of it were 0, so it produced a bunch of discontinuous gradients. As such, I simply check both sides of the graph (have epsilon positive then negative), and if they disagree, I set the value to 0, as in that case I'd just ignore any gradient passing through that node.