Assume we have a matrix $M\in \mathbb R^{t\times qt}$, a vector $p \in \mathbb R^{r}$, and a vector $z \in \mathbb R^{qt}$, .
Note that $\otimes$ is a Kronecker product and $\odot$ is a Hadamard product (Schur product), $\mathbf 1_i$ is a $i$-dimensional all-one vector, and $I_i$ is a $i\times i$ identity matrix.
How to calculate the gradient of the following expression with respect to $z$?
$$ \left( I_{t}\otimes \mathbf 1_{r}^T \right) \left( I_{tr} \odot \left( \left(Mz\otimes \mathbf1_{r} - \mathbf 1_{t}\otimes p\right)\mathbf1_{tr}^T \right) \right) \left(Mz\otimes \mathbf1_{r} - \mathbf 1_{t\times 1}\otimes p \right) $$
Note the following properties of the Diag operation $$\eqalign{ A &= {\rm Diag}(a) = I\odot a{\tt1}^T = I\odot {\tt1}a^T \\ B &= {\rm Diag}(b) \\ Ab &= a\odot b = b\odot a = Ba \\ }$$ and of the Kronecker product of two vectors $$\eqalign{ {\rm vec}(aI_1b^T) &= {\rm vec}(I_aab^T) &= {\rm vec}(ab^TI_b) \\ b\otimes a &= (b\otimes I_a)a &= (I_b\otimes a)b \\ }$$ For typing convenience, define the following variables $$\eqalign{ Q &= (I_{m}\otimes{\tt1}_n) \\ x &= (Mz\otimes{\tt1}_n - {\tt1}_m\otimes p) \;&\in {\mathbb R}^{mn\times 1} \\ X &= {\rm Diag}(x) \\ }$$ Using the above, rewrite the function and calculate its gradient. $$\eqalign{ y &= Q^T(I_{mn}\odot x{\tt1}_{mn}^T)x \;\in {\mathbb R}^{m\times 1} \\ &= Q^TXx\\ dy &= Q^T\,dX\,x + Q^TX\,dx \\ &= 2Q^TX\,dx \\ &= 2Q^TX\,(M\,dz\otimes{\tt1}_n) \\ &= 2Q^TX\,(I_m\otimes{\tt1}_n)M\,dz \\ &= 2Q^TXQM\,dz \\ \frac{\partial y}{\partial z} &= 2Q^TXQM \\ }$$