Gradient Checking: MeanSquareError. Why huge epsilon improves discrepancy?

I am using custom C++ code, and coded a simple “Mean Squared Error” layer.
Temporarily using it for the ‘classification task’, not a simple regression. …maybe this causes the issues?

I don’t have anything else before this layer – not even a simple Dense layer. It’s just MSE on its own.
Its input is a collection of rows of input features. For example, here are 8 rows of input features that will be passed to MSE all at once:

{ a0, a1, a2, a3, a4, a5, a6, a7 }
{ b0, b1, b2, b3, b4, b5, b6, b7 }
{ c0, c1, c2, c3, c4, c5, c6, c7 }
{ d0, d1, d2, d3, d4, d5, d6, d7 }
{ e0, e1, e2, e3, e4, e5, e6, e7 }
{ f0, f1, f2, f3, f4, f5, f6, f7 }
{ g0, g1, g2, g3, g4, g5, g6, g7 }
{ h0, h1, h2, g3, h4, h5, h6, h7 }    //8x8 matrix (contains 64 different values)

Every row of this matrix gets passed into my “Mean Square Error” layer, returning a single scalar for such a row: “Cost”.

I then compute a “final error” scalar, which is the average of such Costs.

When doing Gradient Checking, I am looking at how this “final error” quantity changes as I perturb each of the 64 input values, seen above. The idea is that the changes in finalError must correspond to the gradient computed by formula, with respect to my 64 input values. If they match, then I’ve coded backprop correctly.

Here is the forward prop:

$$finalError = frac{1}{r}sum^r{ left( frac{1}{2n}sum^n{(input_i-target)^2} right) } $$

where $n$ is the number of features per row, and $r$ is the number of rows.

Here is the gradient wrt one of input values, that my backprop is using:

$$frac{partial finalError}{partial input_i} = frac{1}{rn}(input_i – target)$$

Question:

I peturb each input value ‘up’, then ‘down’, running forward prop 64*2 = 128 times. This gives me numerical estimate of gradient for my 64 input values.

However this numerical estimate and the actual analytical gradient become less similar when smaller epsilon is used. This is counterintuitive to me. On the contrary, my vectors match almost exactly, when I use a ridiculously large value for epsilon, such as $1$

Is this expected, or do I have an error in C++ code?

Here is the pseudocode

for every input value i:
   i -= EPSILON
   finalCost_down =  fwdprop( inputMatrix )//very simple - just computes final cost via MSE layer.  finalCost_down is a scalar.
   i += EPSILON
   finalCost_up   =  fwdprop( inputMatrix ) 
   gradientEstimate[i] = (finalCost_up - finalCost_down) / (2*EPSILON)

//after the loop, some time later, just one invocation of backprop:
trueGradientVec = backprop( vec )

//some time later:

discrepancyScalar =  (gradientEstimate - trueGradientVec).magnitude / gradientEstimate.magnitude + trueGradientVec.magnitude)

//somehow discrepancyScalar decreases the larger the EPSILON was used:
// discrepancy is 0.00275, if EPSILON is 0.0001
// discrepancy is 0.00025, if EPSILON is 0.001
// discrepancy is 2.198e-05, if EPSILON is 0.01
// discrepancy is 3.149e-06, if EPSILON is 0.1
// discrepancy is 2.751e-07, if EPSILON is 1

I would expect discrepancy to decrease when epsilon is decreased, because finer perturbations should give more precise slope estimate…

Andrew NG explanation of GradientChecking