I am reading http://www.deeplearningbook.org/ and on chapter $4$ Numerical Computation, at page 94, we read:
Suppose we want to find the value of $\boldsymbol{x}$ that minimizes $$f(\boldsymbol{x}) = \frac{1}{2}||\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b}||_2^2$$ We can obtain the gradient $$\nabla_{\boldsymbol{x}}f(\boldsymbol{x}) = \boldsymbol{A}^T(\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b}) = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} - \boldsymbol{A}^T\boldsymbol{b}$$
So I tried to derive this myself, but didn't quite get there.
My Try
$$f(\boldsymbol{x}) = (\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b})^T(\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b}) = \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} - \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{b} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{b}^T\boldsymbol{b}$$ then since the second and third term are just scalars, their transpose is the same as the other, thus we can cancel them out. Thus we have $$\nabla_xf(\boldsymbol{x}) = \nabla_x(\boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{b}^T\boldsymbol{b}) = ?$$
I really can't continue, I have no idea how to solve that..
My attempt
From above we have $$f(\boldsymbol{x}) = \frac{1}{2} \left(\boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} - \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{b} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{b}^T\boldsymbol{b}\right)$$
From one of the answers below we calculate $$f(\boldsymbol{x} + \boldsymbol{\epsilon}) = \frac{1}{2}\left(\boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon} - \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{b} + \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon}- \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{b} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{x} -\boldsymbol{b}^T\boldsymbol{A}\boldsymbol{\epsilon}+ \boldsymbol{b}^T\boldsymbol{b}\right)$$
Now we notice that the fist is contained in the second, so we can just obtain their difference as $$f(\boldsymbol{x}+\boldsymbol{\epsilon}) - f(\boldsymbol{x}) = \frac{1}{2} \left(\boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon} + \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} + \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{b}-\boldsymbol{b}^T\boldsymbol{A}\boldsymbol{\epsilon}\right)$$
Now we look at the shapes of the matrices. Suppose $\boldsymbol{A}$ has shape (n,m), then $\boldsymbol{x}$ and $\boldsymbol{\epsilon}$ have shape (m,1) and $\boldsymbol{b}$ has shape (n,1). Then the first three terms have shape (1,1), i.e they are scalars. For scalar values, we know that they are equal to their transpose. Also, we replace $\boldsymbol{\epsilon}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon}$ by $\mathcal{O}(\epsilon^2)$. Notice that the transpose of the second term is equal to the first term. Similarly, the transpose of the penultimate term is equal to the last term. Therefore $$f(\boldsymbol{x} + \boldsymbol{\epsilon}) + f(\boldsymbol{x}) = \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\epsilon} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{\epsilon} + \mathcal{O}(\epsilon^2)$$ therefore dividing by $\boldsymbol{\epsilon}$ we have $$\nabla_{\boldsymbol{x}}f(\boldsymbol{x}) = \boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{A} - \boldsymbol{b}^T\boldsymbol{A}$$
Notice that the first term is a vector times a square matrix $\boldsymbol{M} = \boldsymbol{A}^T\boldsymbol{A}$, thus using the property suggested in the comments, we can "transpose it" and the expression is $$\nabla_{\boldsymbol{x}}f(\boldsymbol{x}) = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{x} - \boldsymbol{b}^T\boldsymbol{A}$$
Which is very similar to what I need to obtain, except that the last term is transposed. However, we cannot use the same trick we just used because $\boldsymbol{A}$ doesn't necessarily have to be square!
This is how I differentiate expressions like yours. Let $y = x+\epsilon$. Then, e.g. $$g(y) = y^TAy = x^TAx + x^TA\epsilon + \epsilon^TAx + \epsilon^TA\epsilon$$. Then $$g(x+\epsilon) - g(x) = x^TA\epsilon + x^TA^T\epsilon + O(\epsilon^2).$$ So the gradient is $$x^TA + x^TA^T.$$ The other terms in $f$ can be treated similarly. This approach works because the gradient is related to the linear approximations of a function near the base point $x$.
I looked through your work in response to my answer, and you did it exactly right, except for the transposing bit at the end. It is not actually true that for any square matrix $Mx = x^TM^T$ since the results don't even have the same shape! Also, you can't divide by epsilon, since it is a vector. The right way to finish is to go from $f(x+\epsilon) - f(x) = (x^TA^TA -b^TA)\epsilon$ to concluding that $x^TA^TA -b^TA$ is the gradient (since this is the linear function that epsilon is multiplied by). This is actually the transpose of what you are looking for, but that is just because this approach considers the gradient a row vector rather than a column vector, which is no big deal. Just go ahead and transpose it.