What exactly is averaged when doing batch gradient descent?
Introduction

First of all, it's completely normal that you are confused because nobody really explains this well and accurately enough. Here's my partial attempt to do that. So, this answer doesn't completely answer the original question. In fact, I leave some unanswered questions at the end (that I will eventually answer).

The gradient is a linear operator

The gradient operator $\nabla$ is a linear operator, because, for some $f : \mathbb{R} \rightarrow \mathbb{R} $ and $g: \mathbb{R} \rightarrow \mathbb{R}$, the following two conditions hold.


$\nabla(f + g)(x) = (\nabla f)(x) + (\nabla g)(x),\; \forall x \in \mathbb{R}$
$\nabla(kf)(x) = k(\nabla f)(x),\; \forall k, x \in \mathbb{R}$


In other words, the restriction, in this case, is that the functions are evaluated at the same point $x$ in the domain. This is a very important restriction to understand the answer to your question below!

The linearity of the gradient directly follows from the linearity of the derivative. See a simple proof here. 

Example

For example, let $f(x) = x^2$, $g(x) = x^3$ and $h(x) = f(x) + g(x) = x^2 + x^3$, then $\frac{dh}{dx} = \frac{d (x^2 + x^3)}{d x} = \frac{d x^2}{d x} + \frac{d x^3}{d x} = \frac{d f}{d x} + \frac{d g}{d x} = 2x + 3x$. 

Note that both $f$ and $g$ are not linear functions (i.e. straight-lines), so the linearity of the gradients is not just applicable in the case of straight-lines.

Straight-lines are not necessarily linear maps

Before proceeding, I want to note that there are at least two notions of linearity. 


There's the notion of a linear map (or linear operator), i.e. which is the definition above (i.e. the gradient operator is a linear operator because it satisfies the two conditions, i.e. it preserves addition and scalar multiplication). 
There's the notion of a straight-line function: $f(x) = c*x + k$. A function can be a straight-line and not be a linear map. For example, $f(x) = x+1$ is a straight-line but it doesn't satisfy the conditions above. More precisely, in general, $f(x+y) \neq f(x) + f(y)$, and you can easily verify that this is the case if $x = 2$ and $y=3$ (i.e. $f(2+3) = 6$, $f(2) = 3$, $f(3) = 4$, but $f(2) + f(3) = 7 \neq f(2+3)$.


Neural networks

A neural network is a composition of (typically) non-linear functions (let's ignore the case of linear functions), which can thus be represented as $$y'_{\theta}= f^{L}_{\theta_L} \circ f^{L-1}_{\theta_{L-1}} \circ \dots \circ f_{\theta_1},$$ where 


$f^{l}_{\theta_l}$ is the $i$th layer of your neural network and it 
computes a non-linear function
${\theta_l}$ is a vector of parameters associated with the $l$th layer
$L$ is the number of layers, 
$y'_{\theta}$ is your neural network, 
$\theta$ is a vector containing all parameters of the neural network
$y'_{\theta}(x)$ is the output of your neural network
$\circ $ means the composition of functions


Given that $f^l_{\theta}$ are non-linear, $y'_{\theta}$ is also a non-linear function of the input $x$. This notion of linearity is the second one above (i.e. $y'_{\theta}$ is not a straight-line). In fact, neural networks are typically composed of sigmoids, ReLUs, and hyperbolic tangents, which are not straight-lines.

Sum of squared errors

Now, for simplicity, let's consider the sum of squared error (SSE) as the loss function of your neural network, which is defined as

$$
\mathcal{L}_{\theta}(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^N  \mathcal{S}_{\theta}(\mathbf{x}_i, \mathbf{y}_i) = \sum_{i=1}^N (\mathbf{y}_i - y'_{\theta}(\mathbf{x}_i))^2
$$
where 


$\mathbf{x} \in \mathbb{R}$ and $\mathbf{y} \in \mathbb{R}$ are vectors of inputs and labels, respectively
$\mathbf{y}_i$ is the label for the $i$th input $\mathbf{x}_i$
$\mathcal{S}_{\theta}(\mathbf{x}_i, \mathbf{y}_i) = (\mathbf{y}_i - y'_{\theta}(\mathbf{x}_i))^2$


Sum of gradients vs gradient of a sum

Given the gradient is a linear operator, one could think that computing the sum of the gradients is equal to the gradient of the sums.

However, in our case, we are summing $\mathcal{S}_{\theta}(\mathbf{x}_i, \mathbf{y}_i)$ and, in general, $\mathbf{x}_i \neq \mathbf{x}_j$, for $i \neq j$. So, essentially, the SSE is the sum of the same function, i.e. $S_{\theta}$, evaluated at different points of the domain. However, the definition of a linear map applies when the functions are evaluated at the same point in the domain, as I said above.

So, in general, in the case of neural networks with SSE, the gradient of the sum may not be equal to the sum of gradients, i.e. the definition of the linear operator for the gradient doesn't apply here because we are evaluating every squared error at different points of their domains.

Stochastic gradient descent

The idea of stochastic gradient descent is to approximate the true gradient (i.e. the gradient that would be computed with all training examples) with a noisy gradient (which is an approximation of the true gradient).

How does the noisy gradient approximate the true gradient?

In the case of mini-batch ($M \leq N$, where $M$ is the size of the mini-batch and $N$ is the total number of training examples), this is actually a sum of the gradients, one for each example in the mini-batch. 

The papers Bayesian Learning via Stochastic Gradient Langevin Dynamics (equation 1) or Auto-Encoding Variational Bayes (in section 2.2) use this type of approximation. See also these slides.

Why?

To give you some intuition of why we sum the gradients of the error of each input point $\mathbf{x}_i$, let's consider the case $M=1$, which is often referred to as the (actual) stochastic gradient descent algorithm. 

Let's assume we uniformly sample an arbitrary tuple $(\mathbf{x}_j, \mathbf{y}_j)$ from the dataset $\mathcal{D} = \{ (\mathbf{x}_i, \mathbf{y}_i) \}_{i=1}^N$. 

Formally, we want to show that 

\begin{align}
\nabla_{\theta} \mathcal{L}_{\theta}(\mathbf{x}, \mathbf{y}) 
&amp;=
\mathbb{E}_{(\mathbf{x}_j, \mathbf{y}_j) \sim \mathbb{U}}\left[ \nabla_{\theta} \mathcal{S}_{\theta} \right] \label{1} \tag{1}
\end{align}

where 


$\nabla_{\theta} \mathcal{S}_{\theta}$ is the gradient of $\mathcal{S}_{\theta}$ with respect to the parameters $\theta$
$\mathbb{E}_{(\mathbf{x}_j, \mathbf{y}_j) \sim \mathbb{U}}$ is the expectation with respect to the random variable associated with a sample $(\mathbf{x}_j, \mathbf{y}_j)$ from the uniform distribution $\mathbb{U}$


Under some conditions (see this), we can exchange the expectation and gradient operators, so \ref{1} becomes
\begin{align}
\nabla_{\theta} \mathcal{L}_{\theta}(\mathbf{x}, \mathbf{y}) 
&amp;=
\nabla_{\theta} \mathbb{E}_{(\mathbf{x}_j, \mathbf{y}_j) \sim \mathbb{U}}\left[  \mathcal{S}_{\theta} \right] 
\label{2} \tag{2}
\end{align}
Given that we uniformly sample, the probability of sampling an arbitrary $(\mathbf{x}_j, \mathbf{y}_j)$ is $\frac{1}{N}$. So, equation \ref{2} becomes
\begin{align}
\nabla_{\theta} \mathcal{L}_{\theta} (\mathbf{x}, \mathbf{y}) 
&amp;=
\nabla_{\theta} \sum_{i=1}^N \frac{1}{N} \mathcal{S}_{\theta}(\mathbf{x}_i, \mathbf{y}_i) \\
&amp;=
\nabla_{\theta} \frac{1}{N}  \sum_{i=1}^N \mathcal{S}_{\theta}(\mathbf{x}_i, \mathbf{y}_i) 
\end{align}

Note that $\frac{1}{N}$ is a constant with respect to the summation variable $i$ and so it can be taken out of the summation.

This shows that the gradient with respect to $\theta$ of the loss function $\mathcal{L}_{\theta}$ that includes all training examples is equivalent, in expectation, to the gradient of $\mathcal{S}_{\theta}$ (the loss function of one training example). 

Questions


How can we extend the previous proof to the case $1 &lt; M \leq N$?
Which conditions need exactly to be satisfied so that we can exchange the gradient and the expectation operators? And are they satisfied in the case of typical loss functions, or sometimes they aren't (but in which cases)?
What is the relationship between the proof above and the linearity of the gradient?


In the proof above, we are dealing with expectations and probabilities!

What would the gradient of a sum of errors represent? Can we still use it in place of the sum of gradients?