Gradient of sum of products of matrix traces For a matrix $X \in \Re_{n\times d}$ find the gradient of
$\sum_{i,j}[\langle X_{i.},X_{j.} \rangle\operatorname{tr}(X^TA_{ij}X)]$ w.r.t $X$, 
where $A_{ij}=(e_i-e_j)(e_i-e_j)^T$ using the basis vectors while $X_{i.}$ denotes the $i$'th row. Do note that, $\langle X_{i.},X_{j.} \rangle$ can be written as $\langle X_{i.},X_{j.} \rangle = \operatorname{Tr}(X^Te_ie_j^TX)$ making the original question a sum of products of trace functions.
Hint: The gradient of $\operatorname{tr}(X^TMX)$ w.r.t $X$ for any real matrix $M$ is given by $MX+M^TX$. 
 A: First, let's normalize the notation by rewriting the inner product.
$$\langle X_{i.},X_{j.} \rangle = \mathop{\textrm{Tr}}(X^Te_ie_j^TX) = 
\tfrac{1}{2}\mathop{\textrm{Tr}}(X^T(e_ie_j^T+e_je_i^T)X)$$
Let $B_{ij}\triangleq e_ie_j^T+e_je_i^T$. The expression, which we will call $f(X)$, simplifies to
$$
f(X)=\tfrac{1}{2}\sum_{ij} \mathop{\textrm{Tr}}(X^TB_{ij}X)\mathop{\textrm{Tr}}(X^TA_{ij}X)
$$
The gradient follows from a combination of the standard product rule and the fact that $\nabla X\mathop{\textrm{Tr}}(X^TPX)=2PX$ when $P$ is a symmetric constant. (This is why we took the extra symmetrizing step above.) The gradient is
$$
\nabla_Xf(X)=\sum_{ij} \mathop{\textrm{Tr}}(X^TB_{ij}X)A_{ij}X + \mathop{\textrm{Tr}}(X^TA_{ij}X)B_{ij}X = (Q + R) X,
$$
where $Q$ and $R$ are defined as follows:
$$
Q \triangleq \sum_{ij} \mathop{\textrm{Tr}}(X^TB_{ij}X)A_{ij}, \quad
R \triangleq \sum_{ij} \mathop{\textrm{Tr}}(X^TA_{ij}X)B_{ij}.
$$
Let's find some clean expressions for $Q$ and $R$. For $Q$, we have
$$
\mathop{\textrm{Tr}}(X^TB_{ij}X) 
= \mathop{\textrm{Tr}}(X^T(e_ie_j^T+e_je_i^T)X)
= 2e_i^TXX^Te_j = 2Z_{ij}
$$
where $Z\triangleq XX^T$. Continuing:
$$
Q = \sum_{ij} 2Z_{ij}A_{ij} = \sum_{ij} 2Z_{ij}(e_i-e_j)(e_i-e_j)^T
= \sum_{ij} 2Z_{ij} ( e_ie_i^T + e_je_j^T - e_ie_j^T - e_je_i^T )
$$
For each $(i,j)$, this expression adds $Z_{ij}$ to elements $Q_{ii}$ and $Q_{jj}$, and subtracts $Z_{ij}$ from $Q_{ij}$ and $Q_{ji}$. (When $i=j$, these steps cancel.) The total is then multiplied by two. This will do that:
$$
Q = 2(\mathop{\textrm{diag}}(Z\textbf{1})+\mathop{\textrm{diag}}((\textbf{1}^TZ)^T)-Z-Z^T)=4\mathop{\textrm{diag}}(Z\textbf{1})-4Z
$$
The $\mathop{\textrm{diag}}$ operator constructs a diagonal matrix from a column vector. You can verify this result by substituting $Z\rightarrow Z_{ij}e_ie_j^T$ into the first form for $Q$ above and simplifying; the result should equal the $(i,j)$ summand.
Now consider the second term in the summation:
$$\begin{aligned}
\mathop{\textrm{Tr}}(X^TA_{ij}X) &= \mathop{\textrm{Tr}}(X^T(e_i-e_j)(e_i-e_j)^TX) \\& = (e_i-e_j)^TXX^T(e_i-e_j) = Z_{ii}+Z_{jj}-Z_{ij}-Z_{ji}\end{aligned}$$
$$
R=\sum_{ij} \mathop{\textrm{Tr}}(X^TA_{ij}X)B_{ij} = \sum_{ij} (Z_{ii}+Z_{jj}-Z_{ij}-Z_{ji})(e_ie_j^T+e_ie_j^T)
$$
This summation copies each quantity $Z_{ii}+Z_{jj}-Z_{ij}-Z_{ji}$ to the $(i,j)$ and $(j,i)$ positions. Thus $R$ is
$$R=2(\mathop{\textrm{diag}^*}(Z)\textbf{1}^T+\textbf{1}\mathop{\textrm{diag}^*}(Z)^T-Z-Z^T)=2\mathop{\textrm{diag}^*}(Z)\textbf{1}^T+2\textbf{1}\mathop{\textrm{diag}^*}(Z)^T-4Z.$$
The $\mathop{\textrm{diag}^*}$ operator extracts the diagonal elements of a matrix into a column vector.
The final result, therefore, is
$$
\boxed{
\begin{aligned}
\nabla_X f(X) &= 
(4\cdot\mathop{\textrm{diag}}(XX^T\textbf{1})
+2\cdot\mathop{\textrm{diag}^*}(XX^T)\textbf{1}^T \\ &\quad
+2\cdot\textbf{1}\mathop{\textrm{diag}^*}(XX^T)^T
-8\cdot XX^T)X.
\end{aligned}}
$$
A: Given the standard basis vectors $(e_k)$ and matrices $(E_{ij}=e_ie_j^T)$ construct an all-ones vector and matrix.
$$f = \sum_k e_k, \quad F = \sum_{ij} E_{ij}\quad\implies F=ff^T$$
Define some auxiliary variables
$$\eqalign{
M &= XX^T = M^T &\implies dM = 2\,{\rm Sym}(dX\,X^T) \\
\mu_{ij} &= (E_{ij}:M) = \mu_{ji}
  &\implies M=\sum_{ij} \mu_{ij}E_{ij} \\
B &= I\odot M &\implies dB = I\odot dM \\
b &= {\rm diag}(M) = Bf &\implies bf^T = Bff^T = BF \\
\beta_{k} &= b:e_k \;\doteq \mu_{kk} 
  &\implies b=\sum_{k} \beta_ke_k \\
A_{ij} &= E_{ii} + E_{jj} - E_{ji} - E_{ij} \\
\lambda_{ij} &= (A_{ij}:M) \\
 &= \mu_{ii} + \mu_{jj} - \mu_{ji} - \mu_{ij}  \\
 &= \beta_{i}{\tt1}_j + {\tt1}_i\beta_{j} - 2\mu_{ij}  \\
 &= (E_{ij}:L) 
  &\implies L=\sum_{ij} \lambda_{ij}E_{ij} \\
L &= bf^T + fb^T - 2M \\
}$$
where
$\quad{\rm Sym}$ denotes the symmetrization function  $\;{\rm Sym}(X)=\tfrac{1}{2}(X+X^T)$
 $\quad(\odot)$ denotes the elementwise/Hadamard product
 $\quad(\;:\;)$ denotes the trace/Frobenius product, i.e. 
$\;M\!:\!L={\rm Tr}(M^TL)=\sum\mu_{ij}\lambda_{ij}$
The objective function is now easy to handle.
$$\eqalign{
\phi
 &= \sum_{ij}\mu_{ij}\lambda_{ij} \;= M:L \\
 &= M:(bf^T + fb^T - 2M) \\
 &= 2M:(bf^T - M) \\
 &= 2M:(BF - M) \\
d\phi
 &= 2\,dM:(BF - M) + 2M:(dB\,F - dM) \\
 &= 2(BF - M):dM + 2MF:(I\odot dM)  - 2M:dM \\
 &= \big(2BF - 2M + 2I\odot MF - 2M\big):dM \\
 &= 2\,{\rm Sym}\big(2BF + 2I\odot MF - 4M\big)X:dX \\
\frac{\partial \phi}{\partial X} &= \big(2BF + 2FB + 4I\odot MF - 8M\big)X \\
}$$
