Computing the gradient of $x^TW^TWW^TWx$ I'm working on an independent research project, and need the gradient of an equation for stochastic gradient. I was able to compute the gradient $\nabla_W x^TW^TWx$ which turned out to be 
$
2\begin{pmatrix}
w_1^\top \cdot x\\
w_2^\top \cdot x\\
\vdots \\
w_n^\top \cdot x
\end{pmatrix}
x^\top$
where $x \in \mathbb{R}^n$ and $w_i^\top$ is the $i^{\rm th}$ row of $W$, an $n\times n$ matrix. This is an outer product and produces an $n\times n$ matrix, as required. (I also realized I can write it as $2Wxx^\top$.)
I am now attempting to compute the following more challenging gradient:
$\nabla_W x^TW^TWW^TWx = \nabla_W x^T(W^TW)^2x$
I can convert the gradient in that matrix form to the following summation form:
$\nabla_W \left[\sum_{i=1}^{n}\sum_{j=1}^{n} x_ix_j \sum_{k=1}^{n} (w_i \cdot w_k) (w_k \cdot w_j) \right]$
where here $w_i \in \mathbb{R}^n$ is a column of $W$. To compute gradients, I am trying to do $\partial/\partial w_{ij}$ to get the pattern for one component. However, even for that, the algebra becomes a complete slog here. For instance, with $\partial/\partial w_{ij}$ I am getting
$\sum_{i=1}^{n}\sum_{j=1}^n x_ix_j \frac{\partial}{\partial w_{ij}} \left( (w_i\cdot w_j)\|w_j\|_2^2 + \sum_{k=1, k\ne j}^{n} (w_i\cdot w_k)(w_k \cdot w_j) \right)$
Do you know of a better way to compute gradients that can take advantage of the pattern from the first gradient here, or perhaps another concept from gradient algebra? Or is doing brute-force the way to go? (I've gone through the algebra, but it's difficult for me to see how to condense it nicely.) Maybe it's $2(W^\top W)Wxx^\top$? Thanks for any help you may wish to provide.
 A: There are various ways to do this, I think the simplest is just to expand the expression and pick out the linear parts.
Let $\phi(W) = x^T (W^T W)^2 x$. Then
$\phi(W+H) = x^T ((W+H)^T (W+H))^2 x$.
Expanding and collecting the terms in $H$ we have
$D \phi(W)(H) = x^T (W^T W H^T W + W^T W W^T H + H^T W W^T W + W^T H W^T W) x$.
This can be simplified slightly to $D \phi(W)(H) = 2 x^T (W^T W H^T W + W^T W W^T H) x$.
By looking at $D \phi(W)(e_i e_j^T) $, we obtain the gradient
$\nabla \phi(W) = 2(W^T W x x^T W^T + W W^T W x x^T)$.
Clarifications:
(i) The $H$ is analogous to $h$ in $f(x+h) = f(x) + Df(x)(h) + o(h)$.
One computes $\phi(W+H)-\phi(W)$ and looks for the linear terms in $H$ (and checks the the remaining terms are $o(H)$).
(ii) $D \phi(W)$ (or ${ \partial \phi(W) \over \partial W}$) is the derivative of $\phi$, the only complications are that we are dealing with
matrices.
The derivative evaluated in the direction $H$ is $D\phi(W)(H)$. In the usual
$f: \mathbb{R}^n \to \mathbb{R}$ sort of stuff, we can represent $Df(x)$ as a dual vector (transpose) and $Df(x)(h)$ is just the multiplication by $h$.
The matrices complicate this slightly.
(iii) To talk about a gradient, we need an inner product. In this case
we use $\langle A, B \rangle = \sum_{ij} [A]_{ij} [B]_{ij}$. Then  I need to
find some element $G$ such that $\langle G, H \rangle = D \phi(W)(H)$ and
then $\nabla \phi(W) = G$. We see that
$[G]_{ij} = \langle G, e_i e_j^T \rangle = D \phi(W)(e_i e_j^T)$.
