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$.
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}} $$