Gradient of $x^T B^T B x - x^T B^T b - b^T Bx$ I want to compute the gradient $\nabla_x f(x)$ of $f(x) = x^T B^T B x - x^T B^T b - b^T Bx$ with respect to the vector $x$. So far I have tried below. But when I try to add them together, I couldn't see they come together.
Edits
Now I think I got it. They does come together as the following answers. Since $C=B^TB$, $C$ is symmetric. Then $c_{ij} = c_{ji}$, so 
$$\frac{\partial D}{\partial x} = 2\Big[\sum_{i=1}^n c_{1i} x_i \ \ \sum_{i=1}^n c_{2i} x_i \ \ \cdots \ \ \sum_{i=1}^n c_{ni} x_i\Big] =2 C x.$$
Hence $$\frac{\partial x^TBb}{\partial x} = \Big[
\sum_{i=1}^n b_{i1} b_i \ \
\sum_{i=1}^n b_{i2} b_i \ \
\ \ \cdots \ \ \sum_{i=1}^n b_{in} b_i \ \
\Big] = B^T b.$$
So $\frac{d}{dx} (x^T B^T B x - 2x^T B^T b) = 2B^T B - 2B^Tb$.
Thank you!

Let $C = B^T B$.
\begin{align*}
x^T C x
=& 
\begin{pmatrix} x_1 & x_2 & \cdots & x_n\end{pmatrix}
\begin{pmatrix} 
c_{11} & c_{12} & \cdots & c_{1m}\\ 
c_{21} & c_{22} & \cdots & c_{2m}\\
\vdots & \vdots& \ddots & \vdots \\
 c_{n1} & c_{n2} & \cdots & c_{nm}\\
 \end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\x_n\end{pmatrix}\\
=& 
\begin{pmatrix} x_1 & x_2 & \cdots & x_n\end{pmatrix}
\begin{pmatrix} 
c_{11} x_1 + c_{12} x_2 + \cdots c_{1n} x_n\\
c_{21} x_1 + c_{22} x_2 + \cdots c_{2n} x_n\\
\vdots\\
c_{n1} x_1 + a_{n2} x_2 + \cdots a_{nn} x_n\\
\end{pmatrix}\\
=& x_1(c_{11} x_1 + c_{12} x_2 + \cdots c_{1n} x_n)+x_2(c_{21} x_1 + c_{22} x_2 + \cdots c_{2n} x_n) + \cdots + x_n(c_{n1} x_1 + a_{n2} x_2 + \cdots a_{nn} x_n)
\end{align*}
Since the derivative of the scalar $D = x^T C x$ by a vector $x$ is
$$\frac{\partial D}{\partial x} = \bigg(\frac{\partial D}{\partial x_1} \frac{\partial D}{\partial x_2} \cdots 
\frac{\partial D}{\partial x_n}\bigg),$$
I have:
\begin{align*}
\frac{\partial D}{\partial x_1} 
&= c_{11} x_1 + c_{12} x_2 + \cdots c_{1n} x_n 
+ c_{11} x_1
+ c_{21} x_2
+ \cdots + c_{n1} x_n = \sum_{i=1}^n (c_{1i} + c_{i1})x_i
\end{align*}
Hence $$\frac{\partial D}{\partial x} = \Big[\sum_{i=1}^n (c_{1i} + c_{i1})x_i \ \ \sum_{i=1}^n (c_{2i}+c_{i2})x_i \ \ \cdots \ \ \sum_{i=1}^n (c_{ni}+c_{in})x_i\Big].$$
\begin{align*}
x^T B b
=& 
\begin{pmatrix} x_1 & x_2 & \cdots & x_n\end{pmatrix}
\begin{pmatrix} 
b_{11} b_1 + b_{12} b_2 + \cdots b_{1n} b_n\\
b_{21} b_1 + b_{22} b_2 + \cdots b_{2n} b_n\\
\vdots\\
b_{n1} b_1 + b_{n2} b_2 + \cdots b_{nn} b_n\\
\end{pmatrix}\\
=& 
x_1(b_{11} b_1 + b_{12} b_2 + \cdots b_{1n} b_n)+
x_2(b_{21} b_1 + b_{22} b_2 + \cdots b_{2n} b_n)+\cdots+
x_n(b_{n1} b_1 + b_{n2} b_2 + \cdots b_{nn} b_n)
\end{align*}
Hence $$\frac{\partial x^TBb}{\partial x} = \Big[
\sum_{i=1}^n b_{1i} b_i \ \
\sum_{i=1}^n b_{2i} b_i \ \
\ \ \cdots \ \ \sum_{i=1}^n b_{ni} b_i \ \
\Big].$$

Hint: For a scalar $\epsilon$ and appropriate sized vectors $x,y$ we have: 
$f(x+\epsilon y) = (x+\epsilon y)^TB^TB(x+\epsilon y) - (x+\epsilon y)^TB^Tb - b^TB(x+\epsilon y)$ 
$= x^TB^TBx + \epsilon y^TB^TBx + \epsilon x^TB^TBy + \epsilon^2 y^TB^TBy - x^TB^Tb - \epsilon y^TB^Tb - b^TBx - \epsilon b^TBy$
$= (x^TB^TBx - x^TB^Tb - b^TBx) + \epsilon(y^TB^TBx+x^TB^TBy-y^TB^Tb-b^TBy) + \epsilon^2y^TB^TBy$
$= f(x) + \epsilon(2y^TB^TBx-2y^TBb) + \epsilon^2\text{stuff}$
where we have used the fact that $x^TB^TBy = y^TB^TBx$ and $b^TBy = y^TB^Tb$, since the transpose of a scalar is itself.
