Matrix Calculus Question: a scalar-by-matrix derivative I have the following scalar-by-matrix derivative that I have completely no clue how to solve:
$f(\mathbf{R},\mathbf{S}) = \mathbf{y}^{\top}\bigg(\mathbf{1}\otimes\mathbf{R}+\mathbf{\Phi}\otimes\mathbf{S}\bigg)^{-1}\mathbf{y}$
where $f(\mathbf{R},\mathbf{S})$ is a scalar function, $\mathbf{y}$ is column vector.
Is there a closed form solution to $\dfrac{df(\mathbf{R},\mathbf{S})} {d\mathbf{R}}$ and $\dfrac{df(\mathbf{R},\mathbf{S})} {d\mathbf{S}}$?
Any guidance will be helpful.
Many thanks!
 A: For convenience, let 
$$\eqalign{
M &= 1\otimes R + \Phi\otimes S \cr
P &= -M^{-T}yy^TM^{-T} \cr
}$$
Then write your function in terms of the Frobenius (:) Inner Product and take its differential
$$\eqalign{
 f &= yy^T:M^{-1} \cr
df &= yy^T:dM^{-1} \cr
   &= -yy^T:M^{-1}\,dM\,M^{-1} \cr
   &= P:dM \cr
   &= P:1\otimes dR + P:\Phi\otimes dS \cr
}$$
At this point, assume that we can factor P as
$$\eqalign{
 P &= A\otimes B \cr
}$$
where the matrix A has the same shape as $\Phi$ and B the same shape as S.
$$ $$
Now we can use the Kronecker-Frobenius mixed product rule
$$\eqalign{
 (A\otimes B):(Z\otimes Y) &= (A:Z)\,(B:Y)  \cr
}$$
to write the differential as
$$\eqalign{
df &= (A:1)\,B:dR + (A:\Phi)\,B:dS \cr
}$$
Holding S constant yields the gradient with respect to R
$$\eqalign{
\frac{\partial f}{\partial R} &= \big(1:A\big)\,B \cr
}$$
while holding R constant yields
$$\eqalign{
\frac{\partial f}{\partial S} &= \big(\Phi:A\big)\,B \cr
}$$
You might not be able to find a Kronecker factorization of P, but you can always find a sum such that
$$\eqalign{
 P &= \sum_{k=1}^r A_k\otimes B_k \cr
}$$
Look for the classic paper "Approximation with Kronecker Products" by van Loan and Pitsianis, or Pitsianis' 1997 dissertation which contains Matlab code.
Using a sum changes the results slightly, e.g.
$$\eqalign{
\frac{\partial f}{\partial S} &= \sum_{k=1}^r \big(\Phi:A_k\big)\,B_k \cr
}$$
A: Consider a function involving the Kronecker and Frobenius products
$$\eqalign{
\psi &= P:(Y\otimes X) \\
}$$
Substitute the SVD of $\,P=\sum\sigma_ku_kv_k^T\;$
and calculate the gradient wrt $X$.
$$\eqalign{
\psi
 &= \sum_{k=1}^{\rm rank\,P}\sigma_ku_kv_k^T:(Y\otimes X) \\
 &= \sum_{k=1}^{\rm rank\,P}\sigma_ku_k^T(Y\otimes X)\,v_k \\
 &= \sum_{k=1}^{\rm rank\,P}\sigma_k\,{\rm vec}(U_k)^T(Y\otimes X)\,{\rm vec}(V_k) \\
 &= \sum_{k=1}^{\rm rank\,P}\sigma_k\,{\rm vec}(U_k):{\rm vec}(X\,V_kY^T) \\
 &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_k:X\,V_kY^T \\
 &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_kYV_k^T:X \\
d\psi
 &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_kYV_k^T:dX \\
\frac{\partial\psi}{\partial X}
 &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_kYV_k^T \\
}$$
Applying this result to John's differential relationship
$$\eqalign{
df &= P:{\tt1}\otimes dR + P:\Phi\otimes dS \\
}$$
yields
$$\eqalign{
\frac{\partial f}{\partial R}
  &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_k{\tt\large 1}V_k^T \\
\frac{\partial f}{\partial S}
  &= \sum_{k=1}^{\rm rank\,P}\sigma_kU_k{\Phi}V_k^T \\
}$$
