derivative of the trace of a quantity involving Kronecker product. Let $A\in\mathcal{M}(m\times m)$ and $B\in \mathcal{M}(n\times n)$ and define the Kronecker product $\Sigma=A\otimes B.$ Denote $a\in\Bbb{R.}$
I would like to compute the derivative of $\log(a+(y-\mu)^T\Sigma^{-1}(y-\mu))$ respect to $A.$
There is the ${\rm tr}$ trick, so we have $(y-\mu)^T\Sigma^{-1}(y-\mu)={\rm tr}((y-\mu)^T(y-\mu)\Sigma^{-1}).$
Now we have $\Sigma^{-1}=A^{-1}\otimes B^{-1}$ so that 
$\log(a+(y-\mu)^T\Sigma^{-1}(y-\mu))={\rm tr}((y-\mu)^T(y-\mu)A^{-1}\otimes B^{-1}).$

How can get $$\frac{\partial {\rm tr}((y-\mu)^T(y-\mu)A^{-1}\otimes B^{-1})}{\partial A}$$ ?

 A: For consistency, I'll use the convention were uppercase Latin letters are matrices, lowercase Latin are vectors, and Greek letters are scalars.
I'll also make frequent use of the trace/Frobenius product, i.e.
$$A:B = {\rm tr}(A^TB)$$
Finally, define new variables which are easier to type 
$$\eqalign{
 w &= y-u \cr
 S &= A\otimes B \cr
 \beta &= \alpha + ww^T:S^{-1} \cr
 \lambda &= \log\beta\cr
}$$
Now find the differential of $\lambda$, which is the function that you asked about
$$\eqalign{
d\lambda &= \beta^{-1}d\beta \cr
\beta\,d\lambda
 &= d\beta \cr &= ww^T:dS^{-1} \cr
 &= -ww^T:S^{-1}\,dS\,S^{-1} \cr
 &= -S^{-T}ww^TS^{-T}:dS \cr
 &= -S^{-T}ww^TS^{-T}:dA\otimes B \cr
}$$
At this point, assume we can find a Kronecker factorization of the LHS, i.e.
$$\eqalign{
F\otimes G &= S^{-T}ww^TS^{-T} \cr
}$$
where the sizes of the matrices $(F,G)$ are the same as $(A,B)$, respectively.
Then we can proceed to find the gradient 
$$\eqalign{
\beta\,d\lambda
 &= -(F\otimes G):(dA\otimes B) \cr
 &= -(G:B)(F:dA) \cr
\frac{\partial\lambda}{\partial A} &= -\beta^{-1}(B:G)F \cr
}$$
Our assumed Kronecker factorization may not exist, but we can always find a decomposition by summing over multiple factors
$$\eqalign{
S^{-T}ww^TS^{-T} &= \sum_{k=1}^r F_k\otimes G_k
}$$
Search for papers by vanLoan & Pitsianis on "Kronecker approximation". Pitsianis's 1997 dissertation contains Matlab code to calculate this decomposition.
It's basically (yet another) a clever use of the SVD decomposition.
In this case, since the LHS is a rank-1 matrix and the Kronecker factors are square matrices of sizes $(m\times m)$ and $(n\times n)$, we know exactly how many terms are needed for the decomposition:  $r=\min(m^2,n^2)$ 
This more complicated Kronecker decomposition changes our gradient expression to 
$$\eqalign{
\frac{\partial\lambda}{\partial A} &= -\beta^{-1}\sum_{k=1}^r \big(B:G_k\big)F_k \cr\cr
}$$
NB: There are lots of ways to rearrange the terms in a Frobenius product, all of which follow from the cyclic properties of the trace function.
For example, all of the following are equivalent
$$\eqalign{
A:BC
 &= BC:A  \cr
 &= A^T:(BC)^T \cr
 &= B^TA:C \cr
 &= AC^T:B \cr
}$$
