Gradient Checking: MeanSquareError. Why huge epsilon improves discrepancy?
This was caused by numerical precision of floating numbers. It becomes really apparent as we change one of input values, by a tiny $\epsilon$ which (in my example above) emedially affects the cost function.

The tricky thing was - because the cost function is Mean Square Error (MSE), and there are no other layers in my network, we indeed can use any epsilon to estimate the slope. Even ridiculous epsilon will work, and will be numerically more stable, explaining why discrepancy seemed to get better. That's just how $y=x^2$ works.

In practice, when doing gradient-checking, I switched to another cost function, which is linear (instead of MSE). This improved the numerical stability, allowing me to use epsilon that is 10 times smaller:

$$C=\sum^r\sum^n(obtained-expected)$$

the gradient is then simply:

$$\frac{\partial C}{\partial (obtained)_{rn}} = 1$$

I don't use this cost in production, just for purpose of gradient-checking.

So yeah, this explains why large epsilon was improving the discrepancy.



Edit:

Notice that this linear-cost function must not be used if you have a softmax right before it. That's because the sum of softmax is always 1.0 

In that case you must use the MSE