Why averaging the gradient works in Gradient Descent?
Each training sample ends up in a distant, completely separate location on the error-surface


That is not a correct visualisation of what is going on. The error surface plot is tied to the value of the network parameters, not to the values of the data inputs. During back-propagation of an individual item in a mini-batch or full batch, each example gives an estimate of the gradient in the same location in parameter space. The more examples you use, the better the estimate will be (more on that below).

A more accurate representation of what is going on would be this:



Your question here is still valid though:


  But why does averaging the gathered gradient work? 


In other words, why do you expect that taking all these individual gradients from separate examples should combine into a better approximation of the average gradient over the error surface? 

This is entirely to do with how the error surface is itself constructed as the average of individual loss functions. If we note cost function for the error surface as $C$, then 

$$C(X, \theta) = \frac{1}{|X|}\sum_{x \in X} L(x, \theta)$$

Where $X$ represents the whole dataset, $\theta$ are your model's trainable parameters, $L$ is an individual loss function for $x$. Note I have rolled labels into $X$ here, it doesn't matter for this argument whether loss is due to comparison of model output with some part of the training data - all we care about is finding a gradient to the error surface.

The error gradient that you want to calculate for gradient descent is $\nabla_{\theta} C(X, \theta)$, which you can therefore write as:

$$\nabla_{\theta} C(X, \theta) = \nabla_{\theta}(\frac{1}{|X|}\sum_{x \in X} L(x, \theta))$$

The derivative of the sum of any two functions is the sum of the derivatives, i.e.

$$\frac{d}{dx}(y+z) = \frac{dy}{dx} + \frac{dz}{dx}$$

In addition, any fixed multiplier that doesn't depend on the parameters you are taking the gradient with (in this case, the size of the dataset) can just be treated as an external factor:

$$\nabla_{\theta} C(X, \theta) = \frac{1}{|X|}\sum_{x \in X} \nabla_{\theta} L(x, \theta)$$

So . . . the gradient of an average of many functions, is equal to the average of the gradients of those functions taken separately. Taking any completely random subset of $X$ will result in an unbiased estimate of the mean gradient, same as taking a random subset of any variable and taking its mean will give you an unbiased estimate of the population's mean.

This will not work if your samples are somehow correlated, hence why you will often see recommendations to shuffle the data prior to training, in order to make it i.i.d.

This will also not work if your cost function combines examples in any way other than addition. However, it would be an unusual cost function that combined separate training examples by multiplying their loss functions, or some other non-linear combination.