How does neural network training work, if there are A HUGE number of points that not differentiable? When I first saw ReLu function, I would not guess it will work in neural network because there is a point that is not differentiable. But it seems works very well on modern neural network. 
My question is that, if we have a neural network with ReLu everywhere, then there are HUGE number of points that is not differentiable, then how to update the weights when we do not have the gradient at those points?
 A: The simple and perhaps unsatisfying answer is that we arbitrarily choose a gradient at 0.
Typically deep learning libraries will choose to have a gradient of 0. We can see this using the python libraries of Keras (with tensorflow backend) and numpy:
import keras
import keras.backend as K
import numpy as np

# Define simple model with only ReLU activation.
inp = keras.layers.Input([1])
x = keras.layers.Activation('relu')(inp)
m = keras.models.Model(inp, x)
m.compile(loss='mse', optimizer='sgd')

# get the gradients of the input with respect to the loss.
# This gradient is the gradient through the ReLU activation.
gradients = m.optimizer.get_gradients(m.total_loss, inp)

# create a function we can use to retrieve the gradient
input_tensors = m.inputs + m.sample_weights + m.targets + [K.learning_phase()]
get_gradients = K.function(inputs=input_tensors, outputs=gradients)

# Let's see the gradient of ReLU when passing in [-1, 0, 1, 2]
# as the input, with a target value of always 1
data_input = np.array([[-1.0, 0.0, 1.0, 2.0]]).T
sample_weights = [1.0] * 4 # just scales the gradient
target = np.ones([4, 1]) # The target is 1.0
training_phase = 1 # doesn't effect anything in this example
inputs = [data_input, sample_weights, target, training_phase]
grads = get_gradients(inputs)

X = data_input.flatten()
G = grads[0].flatten()
Y = np.array([1] * len(X))
print('\n'.join(['x = {:+.2f}, loss = {:+.2f}, (∂ loss)/(∂ x) = {:+.8f}'.format(*x)
                 for x in zip(X, (Y - X) ** 2.0, G)]))

This prints the following gradients:
x = -1.00, loss = +4.00, (∂ loss)/(∂ x) = -0.00000000
x = +0.00, loss = +1.00, (∂ loss)/(∂ x) = -0.00000000
x = +1.00, loss = +0.00, (∂ loss)/(∂ x) = +0.00000000
x = +2.00, loss = +1.00, (∂ loss)/(∂ x) = +0.50000000

