The assumption underlying the linear model approach is not that the residual is mean-zero. That statement is tautological. The assumption is that the data, itself, obeys a linear distribution,
$$
y = a x + b + \epsilon,
$$
where $\epsilon$ is a mean-zero residual term and $b$ is a constant. The assumption is in the linear functional form of $y$, not in the fact that $\epsilon$ is mean-zero. The reason I call that latter statement "tautological", is because you can always redefine the component of the equation that's independent of $x$ as follows:
$$
b + \epsilon = c + \delta,
$$
where $c$ is a constant, and $\delta$ is a residual with non-zero mean, but rather, a mean value equal to $b - c$.
Any random variable $\delta$ can be re-written as $(\delta - E[\delta]) + E[\delta]$. We know that $\delta -E[\delta]$ is a RV with zero mean by definition, and $E[\delta]$ is just a constant. In choosing the $x$-independent part to read as "$b + \epsilon$", where $\epsilon$ is mean-zero RV and $b$ is a constant, that's just another way of saying "the portion of $y$ that's independent of $x$ is a random variable, $\delta$, and we'll refer to $\delta - E[\delta]$ as "$\epsilon$", and $E[\delta]$ as "$b$".
It's true that a least-squares optimization is guaranteed to have zero-mean residual on the training set on which you did the optimization. However, that is independent of the assumption. When you say that your assumption is $y = a x + b + \epsilon$, that's not an assumption about your model, but rather, an assumption about the data itself. You're only training a linear model to learn what those parameters are.