# Numerical gradient checking (best practices)

I've implemented a neural network and am using numerical gradient checking to validate the back-propagation algorithm is working correctly.

I'm using the standard method to calculate the numerical gradient:

$f'(x) \approx \frac{J(\theta + \epsilon) - J(\theta - \epsilon)}{2 \epsilon}$ where $\epsilon = 10^{-4}$, and

norm(gradients - numericalGradients)/norm(gradients + numericalGradients) < 10e-8


to check back-propagation is operating correctly.

I'm currently checking throughout training for testing purposes and I think that it is functioning correctly, however, as the network approaches its solution the gradients become very small and I think this is why the check fails.

So my question is what are the best practices for when to perform gradient checking?

Once at the start only, only when the minimum gradients are N times bigger than epsilon?

#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

class NeuralNet:

EPSILSON = 10e-4

def __init__(self, _maxIter=250, _nHidden=3):
# Set hyperparameters
self.maxIter = _maxIter
self.nHidden = _nHidden
self.learningRate = np.linspace(0.5, 0.05, self.maxIter)
self.momentum = 0.5
self.enNesterov = False  # Nesterov accelerated gradient (https://arxiv.org/pdf/1212.0901v2.pdf)
self.cost = []
self.dW1 = 0
self.dW2 = 0

def initNet(self, _nInput, _nOutput):
self.nInput = _nInput
self.nOutput = _nOutput
self.W1 = np.random.rand(self.nInput, self.nHidden)/100
self.W2 = np.random.rand(self.nHidden, self.nOutput)/100

def train(self, X, y):
# Initialise network structure
self.initNet(X.shape[1], y.shape[1])

# Train - full batch
for m in range(self.maxIter):

# Check gradient descent is operating as expected

# Train net
yHat = self.feedforward(X)
J = self.costFunction(y, yHat)
self.cost.append(J)
dJdW1, dJdW2 = self.backprop(X, y, yHat)
self.updateWeights(dJdW1, dJdW2, self.learningRate[m])

def feedforward(self, X):
self.z2 = np.dot(X, self.W1)
self.a2 = self.sigmoid(self.z2)
self.z3 = np.dot(self.a2, self.W2)
yHat = self.sigmoid(self.z3)
return yHat

def costFunction(self, y, yHat):
J = 0.5 * np.sum((y - yHat)**2)
return J

def backprop(self, X, y, yHat):
# NB: delta 3 depends on the cost function applied
delta3 = np.multiply(-(y - yHat), self.sigmoidPrime(self.z3))
dJdW2 = np.dot(self.a2.T, delta3)
delta2 = np.dot(delta3, self.W2.T) * self.sigmoidPrime(self.z2)
dJdW1 = np.dot(X.T, delta2)
return dJdW1, dJdW2

def updateWeights(self, dJdW1, dJdW2, learnRate):
if self.enNesterov:
dW1_prev = self.dW1
dW2_prev = self.dW2
self.dW1 = learnRate*dJdW1 + self.momentum*self.dW1
self.dW2 = learnRate*dJdW2 + self.momentum*self.dW2
self.W1 = self.W1 - (1+self.momentum)*self.dW1 - self.momentum*dW1_prev
self.W2 = self.W2 - (1+self.momentum)*self.dW2 - self.momentum*dW2_prev
else:
self.dW1 = learnRate*dJdW1 + self.momentum*self.dW1
self.dW2 = learnRate*dJdW2 + self.momentum*self.dW2
self.W1 = self.W1 - self.dW1
self.W2 = self.W2 - self.dW2

def getWeights(self):
return np.concatenate((self.W1.ravel(), self.W2.ravel()))

def setWeights(self, weights):
W1_start = 0
W1_end = self.nInput * self.nHidden
self.W1 = np.reshape(weights[W1_start:W1_end], (self.nInput, self.nHidden))
W2_end = W1_end + self.nHidden * self.nOutput
self.W2 = np.reshape(weights[W1_end:W2_end], (self.nHidden, self.nOutput))

yHat = self.feedforward(X)
dJdW1, dJdW2 = self.backprop(X, y, yHat)

# Compare
str = 'PASS'
else:
str = 'FAIL'
print('[{0}] Gradient checking. diff = {1}'.format(str, diff))

weights = self.getWeights()
perturb = np.zeros(weights.shape)

for p in range(len(weights)):
# Set pertubation for this weight only
perturb[p] = self.EPSILSON
# Positive perturbation
self.setWeights(weights + perturb)
yHat = self.feedforward(X)
Jpos = self.costFunction(y, yHat)
# Negative perturbation
self.setWeights(weights - perturb)
yHat = self.feedforward(X)
Jneg = self.costFunction(y, yHat)
numGrad[p] = (Jpos - Jneg) / (2 * self.EPSILSON)
# Reset perturbation for next iteration
perturb[p] = 0

# Reset weights
self.setWeights(weights)

def sigmoid(self, z):
return 1 / (1 + np.exp(-z))

def sigmoidPrime(self, z):
x = np.exp(-z)
return (x / ((1 + x)**2))

#  Data
X = np.array(([3, 5], [5, 1], [10, 2]), dtype=float)
y = np.array(([75], [82], [93]), dtype=float)

# Normalize
X = X/np.amax(X, axis=0)
y = y/100

# Train a network
NN = NeuralNet()
NN.train(X, y)

h = NN.feedforward(X)

# PLotting
f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=False)
ax1.plot(NN.cost)
ax1.set(xlabel='Iteration', ylabel='Cost', title='')
ax1.grid()

• Have you tried smaller values for $\epsilon$ (e.g. $10^{-8}$)? Near a minimum the gradient usually changes more quickly so a large step size might yield too coarse of a numerical approximation. – Ruben van Bergen Mar 6 '18 at 20:54
• @RubenvanBergen Thanks, yes I have played with the epsilon value, and indeed $10^{-8}$ produces a pass to a tolerance of $10^{-8}$ throughout training as the gradients tend towards their minima. Is the test just as valid with a much smaller epsilon (within numerical tolerances)? I read $10^{-4}$ throughout the literature and am wondering why this value? – CatsLoveJazz Mar 6 '18 at 21:04
• Personally, as a general rule I always try to use as small a step-size as possible (i.e. without causing numerical underflow issues) to get the best approximation. I don't know why $10^{-4}$ would be special - perhaps if you use some kind of standardization for your variables/cost function this is a step-size that is usually good enough? But yes, as far as I know the test is just as valid, if not more so, with smaller $\epsilon$. – Ruben van Bergen Mar 7 '18 at 7:57