How can concentrated (profile) log marginal likelihood be used to optimize the mean and scale(outputscale) parameters in Gaussian Process Regression?

Question

The log marginal likelihood which is used in Gaussian Process Regression comes from a Multivariate Normal pdf Gaussian Processes for Machine Learning, p.19, eqn. 2.30, Surrogates, Chapter 5, eqn. 5.4

\begin{equation} \log \: p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2}\mathbf{\left(y - \mu\right)}^\top (\mathbf{K} + \sigma^2_n\mathbf{I})^{-1}\mathbf{\left(y - \mu\right)} \ - \frac{1}{2}\log |\mathbf{K} + \sigma^2_n\mathbf{I}| \ - \frac{n}{2} \log 2\pi \label{log_marginal_likelihood} \tag{1} \end{equation}

Both books linked above assume a zero mean GP (so $\mu$ is $\mathbf{0}$). Most GP implementations assume zero mean, as the training data (both inputs and labels) are standardized prior to training (and standardization is reversed when making predictions), but I have also seen constant mean in the documentation of libraries such as scikit-learn and gpytorch, where $\mu$ is also optimized by gradient descent as another model parameter.

However, this book A Taxonomy of Global Optimization Methods Based on Response Surfaces, page 357, eqns. 9-13 gives an analytic solution for $\mu$ by setting the derivative of the log likelihood function w.r.t. $\mu$ to 0 and solving for $\mu$:

\begin{equation} \hat \mu = \frac{\mathbf{1}^\top \: \left(\mathbf{R}\right)^{-1}\mathbf{y}} {\mathbf{1}^\top \: \left(\mathbf{R}\right)^{-1}\mathbf{1}} \label{predicted_mu} \tag{2} \end{equation}

Note that the likelihood in this book is given as (ignoring constant terms)

\begin{equation} -\frac{n}{2}\log{(\tau^2)} -\frac{1}{2}\log(|\mathbf{R}|) -\frac{\mathbf{\left(y - \mu\right)}^\top (\mathbf{R})^{-1}\mathbf{\left(y - \mu\right)}} {2\tau^2} \label{log_likelihood_357} \tag{3} \end{equation}

where $cov(Y)=\tau^2\mathbf{R}$. I am not sure why the outputscale is separated from the covariance matrix in Equation \ref{log_likelihood_357}, whereas (I assume) it is included in $\mathbf{K}$ in equation \ref{log_marginal_likelihood}.

$\hat \mu$ can then be plugged in to estimate the scale (outputscale) $\tau^2$ as follows

\begin{equation} \hat \tau^2 = \frac{\left(\mathbf{y}-\mathbf{1}\hat\mu \right)^\top \left(\mathbf{R}\right)^{-1} \left(\mathbf{y}-\mathbf{1}\hat\mu \right)}{n} \label{predicted_tau_squared} \tag{4} \end{equation}

where $n$ is the number of training samples. This is the same equation as in Surrogates (Equation 5.5), although $\mathbf{0}$ mean is assumed there, so $\hat\mu = 0$.

Finally, the estimated values of $\hat\mu$ and $\hat\tau^2$ are plugged in Equation \ref{log_likelihood_357} to give the concentrated (profile) log likelihood (ignoring constant terms). Accoding to the book, this is the function that is maximized to estimate model hyperparameters (lengthscales, etc.). I am not sure why the third term from Equation \ref{log_likelihood_357} is missing.

\begin{equation} -\frac{n}{2}\log(\hat\tau^2) -\frac{1}{2}\log(|\mathbf{R}|) \label{profile_likelihood} \tag{5} \end{equation}

My questions regarding the concentrated log likelihood are:

1. What are the benefits of estimating $\hat \mu$ and $\hat\tau^2$ this way, is it just to reduce the number of parameters that need to be optimized numerically?

2. In what order should parameters be optimized if using the profile likelihood \ref{profile_likelihood}? The book says to optimize lengthscales, noise, etc. first, and then plug the optimized values in Equations \ref{predicted_mu} and \ref{predicted_tau_squared}? This does not make a lot of sense to me as the profile likelihood has $\tau^2$ in it already.

3. Why is $\hat\mu$ in Equation \ref{predicted_mu} estimated this way? To me Equation \ref{predicted_mu} looks very much like generalized least squares, but I was not able to find other sources that use this equation or a derivation of it. Since the log-marginal likelihood comes from a MVN, then wouldn't $\hat \mu$ just be the Maximum Likelihood Estimate of the Multivariate Gaussian given as

\begin{equation} \bar y = \frac{1}{n}\sum_{i=1}^n y_i \tag{6} \label{mean_mvn} \end{equation}

as derived in another CrossValidated answer. Then the GP constant mean vector would just be $1\bar y$.

Answer 1

1. You can estimate the model parameters $\mu$ & $\tau^2$ using the MLE approach (along with other parameters of the kernel) or simply using the analytical GLS solutions.

If you choose MLE, while using gradient-based algorithms (like L-BFGS-B), you need to provide good initialization points for $\mu$ & $\tau^2$, in order to obtain robust estimates. The initialization is crucial since the non-convexity of the likelihood often causes early stopping and you may end up with very poor estimates for the parameters.

In the other approach (recommended in this article), you first initialize the lengthscale values from a grid and use the corresponding GLS solutions for the initialization of $\mu$ & $\tau^2$. Finally choose the one from the grid, that yields the best likelihood value.

2. see the last paragraph of the previous answer.
3. The derivation in CrossValidated answer. (they derive the MLE for the mean vector $\mu$ for a MVN distribution, given m observations) is totally different from the context here (for a constant mean GP we are interested in a scalar $\mu$). See the construction of the likelihood and it will be clear.