Matrix calculus : Derivative with respect to two vectors

I have a function of two variables, $$\vec{x_i} \in \mathbb{R}^n$$ and $$\vec{y_j} \in \mathbb{R}^n$$:

$$f(\vec{x_i},\vec{y_j}) =\sum_\limits{i,j} \left( \vec{x_i} \cdot \vec{y_j} + z_{ij} \right)^2$$

I am trying to find the second derivative of this expression. I understand how to accomplish the task when keeping the $$\vec{x}_i$$ or $$\vec{y}_j$$ constant but I cannot figure out a way to do it when both of the variables are not constant. I would really appreciate being pointed to a certain direction.

Let $$(e_k)$$ be the standard basis vectors and $$(X,Y)$$ be the matrices whose columns are the $$(x_i,\,y_j)$$ vectors, i.e. \eqalign{ x_i &= Xe_i, &y_j = Ye_j \cr x_i\cdot y_j &= e_i^TX^TY&e_j = X^TY:e_ie_j^T \cr } where colon represents the trace/Frobenius product, i.e. $$\,\,A:B={\rm Tr}(A^TB)$$
Define $$K=\begin{bmatrix}0&I\\0&0\end{bmatrix}$$ and concatenate the XY-matrices into a single matrix $$\,\,W^T=[X^T\,\,Y^T]$$
Note that $$X^TY=W^TKW$$.
Using these new variables, rewrite the cost function and find its differential and gradient. \eqalign{ f &= (X^TY+Z):(X^TY+Z) \cr &= (W^TKW+Z):(W^TKW+Z) \cr df &= 2(W^TKW+Z):d(W^TKW) \cr &= 2(W^TKW+Z):(dW^TKW+W^TK\,dW) \cr &= 2\Big(KW(W^TKW+Z)^T+K^TW(W^TKW+Z)\Big):dW \cr G=\frac{\partial f}{\partial W} &= 2KW\big(W^TKW+Z\big)^T+2K^TW\big(W^TKW+Z\big) \cr } Now find the differential of the gradient \eqalign{ dG &= 2KdW\big(W^TKW+Z\big)^T+2K^TdW\big(W^TKW+Z\big) \cr &+\,2KW\big(dW^TKW\big)^T+2K^TW\big(dW^TKW\big) \cr &+\,2KW\big(W^TKdW\big)^T+2K^TW\big(W^TKdW\big) \cr \cr } Shuffling these terms around to produce the hessian is a tedious exercise, which I'll leave to you.
Keep in mind that, since the gradient is a matrix, the resulting hessian will be a 4th order tensor. $${\mathcal H} = \frac{\partial G}{\partial W} = \frac{\partial^2f}{\partial W\partial W}$$ One way to handle the term-shuffling is to flatten the expressions using vectorization, i.e. $$A\,dX\,B = (B^T\otimes A)\,{\rm vec}(dX)$$ where $$\otimes$$ is the Kronecker product.
You'll also need a special permutation matrix $$P\,$$ called the Commutation matrix to reduce some of the intermediate terms. $${\rm vec}(dX^T)=P\,{\rm vec}(dX)$$