Let $M \in \mathbb{C}^{vr \times vc}$ and $P \in \mathbb{R}^{v \times v}$, where $P \succ 0 $
then,
$F = Tr(M (P^{-1} \otimes I_{c}) M^H (P \otimes I_{r}))$,
where,
$I_{r} \in \mathbb{R}^{r \times r}$ and $I_{c} \in \mathbb{R}^{c \times c}$ are identity matrices, respectively. $Tr()$ is the trace of a matrix.
We want to compute $\frac{dF}{dP}$. I know that we can decompose $M = \sum_i B_i \otimes E_i$ and then use the properties of a Kronecker product to distribute the trace, which would yield $F = \sum_{ij} Tr(B_iP^{-1}B_j^HP) Tr(E_i E_j^H)$ and then take the derivative $\frac{dF}{dP}$, which is easily done.
However, is it possible to compute $\frac{dF}{dP}$ without the decomposition of $M$?
The end goal is to solve $\frac{dF}{dP} = 0$
$ \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\e{\varepsilon} \def\l{j} \def\I{I_{\tt1}} \def\C{C^{-\tt1}} \def\Ct{C^{-T}} \def\bbR#1{{\mathbb R}^{#1}} \def\L{\left}\def\R{\right}\def\LR#1{\L(#1\R)} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\p{\partial}\def\o{{\,\big[\tt1}\big]\,} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\qiq{\quad\implies\quad} $Use a colon to denote the Frobenius product, which is a concise notation for the trace $$\eqalign{ A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij} \;=\; \trace{A^TB} \\ A:A &= \|A\|^2_F \\ }$$ The properties of the underlying trace function allow the terms in a Frobenius product to be rearranged in many different but equivalent ways, e.g. $$\eqalign{ A:B &= B:A \\ A:B &= A^T:B^T \\ C:\LR{AB} &= \LR{CB^T}:A = \LR{A^TC}:B \\ }$$ Use $\{e_\l,\;\e_k\}$ to denote the euclidean basis vectors for $\{\bbR{r},\bbR{c}\}$ respectively, which can be used to write identity matrices as rank-one sums $$\eqalign{ \I &= \o \\ I_r &= \sum_{\l=1}^r e_\l e_\l^T = \sum_{\l=1}^r e_\l\,I_1\,e_\l^T \\ I_c &= \sum_{k=1}^c \e_k \e_k^T = \sum_{k=1}^c \e_k \,I_1\,\e_k^T \\ }$$ The above can be used to rearrange the trace of a Kronecker product $$\eqalign{ \trace{A^T\LR{P\otimes I_c}} &= A:\LR{P\otimes I_c} \\ &= \sum_{k=1}^c A:\LR{P\otimes \e_k\I\e_k^T} \\ &= \sum_{k=1}^c A :\LR{\LR{I_v\otimes\e_k}\LR{P\otimes\I}\LR{I_v\otimes\e_k}^T} \\ &= \sum_{k=1}^c A :\LR{\LR{I_v\otimes\e_k}\,P\,\LR{I_v\otimes\e_k}^T} \\ &= \LR{\sum_{k=1}^c \LR{I_v\otimes\e_k}^TA\LR{I_v\otimes\e_k}}:P \\ }$$
For typing convenience, define the matrix variables $$\eqalign{ C &= P\otimes I_c &\qiq dC = dP\otimes I_c \\ R &= P\otimes I_r &\qiq dR = dP\otimes I_r \\ A &= \LR{\C M^HRM\C}^T \\ B &= \LR{M\C M^H}^T \\ }$$ Calculate the differential and gradient of the function. $$\eqalign{ F &= M^T:\LR{\C M^HR} \\ dF &= M^T:\LR{\C M^H\,dR} + M^T:\LR{d\C\,M^HR} \\ &= \LR{M^*\Ct M^T}:{dR} - {M^T}:\LR{\C\,dC\,\C M^HR} \\ &= \LR{M^*\Ct M^T}:{dR} - \LR{\Ct M^TR^TM^*\Ct}:{dC} \\ &= B:\LR{dP\otimes I_r} - A:\LR{dP\otimes I_c} \\ &= \LR{\sum_{\l=1}^r \LR{I_v\otimes e_\l}^TB\LR{I_v\otimes e_\l}}:dP - \LR{\sum_{k=1}^c \LR{I_v\otimes\e_k}^TA\LR{I_v\otimes\e_k}}:dP \\ \grad{F}{P} &= \LR{\sum_{\l=1}^r \LR{I_v\otimes e_\l}^TB\LR{I_v\otimes e_\l}} - \LR{\sum_{k=1}^c \LR{I_v\otimes\e_k}^TA\LR{I_v\otimes\e_k}} \\ }$$ So that's the gradient expression if you decompose $I_c$ and $I_r$ rather than $M$.
In either case, a closed-form solution of the zero-gradient condition is impossible. However, since you have an exact formula for the gradient of $F$, it would be a simple exercise to use your favorite gradient-based solver (e.g. Conjugate Gradients, Barzilai-Borwein, etc) to numerically calculate the optimal value(s) of $P$.