The existence of maximum likelihood estimator on $\mathbb{R}^n$ $\newcommand{\tr}{\operatorname{tr}}$For a matrix-variate Gaussian distribution, the negative log marginal likelihood is
\begin{equation}\label{matrixLikelihood}
  \mathcal{L} = \frac{nd}{2}\ln(2\pi) + \frac{d}{2}\ln \det(K') + \frac{n}{2}\ln \det(\Omega) + \frac{1}{2}\tr((K')^{-1}Y\Omega^{-1}Y^{\mathrm{T}}),
\end{equation}
where $n$ is the number of input points, d is the dimension of output $Y \in \mathbb{R}^{n \times d}$, $K' = K + \sigma^2_n I$, $K \in \mathbb{R}^{n \times n }$ is the column
covariance matrix of matrix-variate distribution with many kernel parameters $\theta_i$, for example, squared exponential (SE) kernel, $K = [k_{ij}], k_{ij} = \theta_1 \exp( - \frac{- (x_i - x_j)^2}{2 \theta_2^2})$,
$\Omega \in \mathbb{R}^{d \times d}$ is row covariance matrix. In addition, both $K$ and $\Omega$ are postive semi-definite, and thus we can write $\Omega = \Phi \Phi^{\mathrm{T}}$, where
$$
\Phi =
\left[
\begin{matrix}
 \phi_{11}      & 0      & \cdots & 0      \\
 \phi_{21}      & \phi_{22}      & \cdots & 0      \\
 \vdots & \vdots & \ddots & \vdots \\
 \phi_{d1}      & \phi_{d2}      & \cdots & \phi_{dd}      \\
\end{matrix}
\right].
$$
To guarantee the uniqueness of $\Phi$, the diagonal elements are restricted to be positive and denote $\varphi_{ii} = \ln(\phi_{ii})$ for $i = 1,2,\cdots,d$.
How to prove the existence of a maximum likelihood estimator (MLE) with respect to parameter $\sigma_n^2$, $\theta_i$, $\phi_{ij}$ and $\varphi_{ii}$.
Here our parameter space are $\mathbb{R}^n_{+}$ for $\sigma_n^2$, $\theta_i$, and $\varphi_{ii}$ and $\mathbb{R}^n$ for $\phi_{ij}$.
I was trying to use the theorem that the continuous function on a compact set must have minimum and maximum. But the problem is now our parameter space is not compact. Therefore, I tried to check the monotonicity of the likelihood function at infinity, so I obtain the derivative:
with respect to parameter $\sigma_n^2$, $\theta_i$, $\phi_{ij}$ and $\varphi_{ii}$ are as follows
\begin{align*}
  \frac{\partial \mathcal{L}}{\partial \sigma_n^2} &=  \frac{d}{2}\tr((K')^{-1}) - \frac{1}{2}\tr(\alpha_{K'}\Omega^{-1}\alpha_{K'}^{\mathrm{T}}),\\
  \frac{\partial \mathcal{L}}{\partial \theta_i} & = \frac{d}{2}\tr\left((K')^{-1}\frac{\partial K_{\theta}}{\partial \theta_i}\right) - \frac{1}{2}\tr\left(\alpha_{K'}\Omega^{-1}\alpha_{K'}^{\mathrm{T}}\frac{\partial K_{\theta}}{\partial \theta_i}\right), \\
   \frac{\partial \mathcal{L}}{\partial \phi_{ij}} & = \frac{n}{2}\tr[\Omega^{-1}(\mathbf{E}_{ij}\Phi^{\mathrm{T}} + \Phi \mathbf{E}_{ij})] - \frac{1}{2}\tr[\alpha_{\Omega}(K')^{-1}\alpha_{\Omega}^{\mathrm{T}}(\mathbf{E}_{ij}\Phi^{\mathrm{T}} + \Phi \mathbf{E}_{ij})], \\
  \frac{\partial \mathcal{L}}{\partial \varphi_{ii}} & = \frac{n}{2}\tr[\Omega^{-1}(\mathbf{J}_{ii}\Phi^{\mathrm{T}} + \Phi \mathbf{J}_{ii})] - \frac{1}{2}\tr[\alpha_{\Omega}(K')^{-1}\alpha_{\Omega}^{\mathrm{T}}(\mathbf{J}_{ii}\Phi^{\mathrm{T}} + \Phi \mathbf{J}_{ii})],
\end{align*}
where $\alpha_{K'} = (K')^{-1}Y$, $\alpha_{\Omega} = \Omega^{-1}Y^{\mathrm{T}}$, $\mathbf{E}_{ij}$ is the $d \times d$ elementary matrix having unity in the (i,j)-th element and zeros elsewhere, and $\mathbf{J}_{ii}$ is the same as $\mathbf{E}_{ij}$ but with the unity being replaced by $e^{\varphi_{ii}}$.
But I don't know how to evaluate these derivatives positive or negative. 
Similarly, is it possible to prove the MLE existence for matrix-variate student-t distribution
\begin{eqnarray*}\label{MultiLikelihoodT}
% \nonumber % Remove numbering (before each equation)
  \mathcal{L}  &=& \frac{1}{2}(\nu+ d+n -1) \ln \det(\mathbf{I}_n + (K')^{-1}Y\Omega^{-1}Y^{\mathrm{T}}) + \frac{d}{2}\ln \det(K') + \frac{n}{2}\ln \det(\Omega) \\
               & & {} + \ln\Gamma_n \left(\frac{1}{2}(\nu + n -1)\right) - \ln \Gamma_n \left(\frac{1}{2}(\nu + d + n -1)\right) + \frac{1}{2}dn\ln\pi  \\
               &=& \frac{1}{2}(\nu+ d+n -1) \ln \det(K' +Y\Omega^{-1}Y^{\mathrm{T}}) - \frac{\nu + n -1}{2}\ln \det(K')  \\
               & & {} + \ln\Gamma_n \left(\frac{1}{2}(\nu + n -1)\right) - \ln \Gamma_n \left(\frac{1}{2}(\nu + d + n -1)\right)+ \frac{n}{2}\ln \det(\Omega)+ \frac{1}{2}dn\ln\pi,
\end{eqnarray*}
where $\nu$ is the degree of freedom of student-t distribution. 
\begin{align*}
  \frac{\partial \mathcal{L}}{\partial \nu} &= \frac{1}{2}\ln \det(U) - \frac{1}{2}\ln \det(K')  + \frac{1}{2}\psi_n(\frac{1}{2}\tau) - \frac{1}{2}\psi_n\left(\frac{1}{2}(\tau + d)\right), \\
  \frac{\partial \mathcal{L}}{\partial \sigma^2_n}
   & = \frac{(\tau+d)}{2}\tr(U^{-1}) - \frac{\tau}{2}\tr((K')^{-1}), \\
  \frac{\partial \mathcal{L}}{\partial \theta_i} & = \frac{(\tau+d)}{2}\tr\left(U^{-1} \frac{\partial K_{\theta}}{\partial \theta_i}\right) - \frac{\tau}{2}\tr\left(\Sigma^{-1} \frac{\partial K_{\theta}}{\partial \theta_i}\right),\\
    \frac{\partial \mathcal{L}}{\partial \phi_{ij}} & = - \frac{(\tau +d)}{2}\tr[U^{-1}\alpha_{\Omega}^{\mathrm{T}}(\mathbf{E}_{ij}\Phi^{\mathrm{T}} + \Phi \mathbf{E}_{ij})\alpha_{\Omega}] + \frac{n}{2}\tr[\Omega^{-1}(\mathbf{E}_{ij}\Phi^{\mathrm{T}} + \Phi \mathbf{E}_{ij})], \\
  \frac{\partial \mathcal{L}}{\partial \varphi_{ii}} & = -\frac{(\tau +d)}{2}\tr[U^{-1}\alpha_{\Omega}^{\mathrm{T}}(\mathbf{J}_{ii}\Phi^{\mathrm{T}} + \Phi \mathbf{J}_{ii})\alpha_{\Omega}] + \frac{n}{2}\tr[\Omega^{-1}(\mathbf{J}_{ii}\Phi^{\mathrm{T}} + \Phi \mathbf{J}_{ii})],
\end{align*}
where $U = K' + Y\Omega^{-1}Y^{\mathrm{T}}$, $\tau = \nu + n -1$ and $\psi_n(\cdot)$ is the derivative of the function $\ln \Gamma_n(\cdot)$ with respect to $\nu$.
 A: My best advice to prove the existence of min/max for a function of matrix parameters is to find the gradient with respect to relevant parameters (as you did), set it to 0, solve for the relevant parameters to find the stationary point/s, and then verify that, at that/those stationary point/s, the Hessian matrix (calculated at that/those stationary points) is positive definite for a minimum and negative definite for a maximum. So since it seems you have already calculated the gradient, what is missing is just that you show that you can solve $\frac{\delta f}{\delta params}=0$ . Then find the Hessian calculated for the optimal values of the params found at previous step. Then show that it is pd/nd by looking at the sign of the determinant. In case the determinant expression still depends on some values of the params, show that the determinant is always >0 or <0 for any value of such parameters. Almost The same holds for a constrained optimization where you take the gradient of the Lagrangian and set it to 0. Then however you have to find the bordered hessian instead of the simple hessian. But logically the process is the same.
Notice that you are fine working with the marginal log likelihood to simplify the problem of gradient calculations. But then remember that you want to max/min the loglikelihood/-loglikelihood, thus the gradient of loglikelihood/-loglikelihood will be the sum of the gradients of the marginal loglikelihood/-loglikelihood if the parameters are the same for each marginal loglikelihood/-loglikelihood as in your example.
See also this question for more conceptual review, and this for the issue of uniqueness of the solution
