# Gaussian process Optimization

The Gaussian process model (GP) is written as $$y(x)=h(x)^{t}\boldsymbol\beta + f(x)+\epsilon(x)\qquad[1]$$ where $h(x)^{t}$ is a regression function such as $\left [1,x,x^{2} \right ]$ ; $\boldsymbol\beta$ is the set of regression coefficients and $\epsilon(x)$ is normal with zero mean and $\sigma^{2}$ variance and $f_{i}(x)$ is a Gaussian process with a covariance matrix $K_{\theta}$ $$f_{i}(x)\sim GP(0,K_{\theta}(x,x^{'}))\qquad[2]$$ such that $\theta$ is the set of hyper parameters corresponding to the covariance matrix.

To estimate $\theta$ and $\boldsymbol\beta$ we need to maximize the multivariate Gaussian likelihood function where the joint distribution of $y(x)$ is $$y(x)\sim GP(h(x)^{t}\boldsymbol\beta,K_{\theta}(x,x^{'})+\sigma^{2}I)\qquad[3]$$

Where the best linear estimator of $\boldsymbol\beta$ is $${\widehat{\boldsymbol\beta}} = (\mathbf{F}^{t}K_{\theta}^{-1}\mathbf{F})^{-1}\mathbf{F}^{t}K_{\theta}^{-1}\mathbf{y}\qquad[4]$$ where $\mathbf{F}=h(\mathbf{x})^{t}$

Currently to estimate $\theta$ and $\boldsymbol\beta$, I optimize the multivariate Gaussian likelihood in equation $[3]$. I am wondering how can the results in $[4]$ be utilized in equation $[3]$ so that some simplification might be possible.

• Are you doing this numerically? It is possible your problem falls into the class of separable nonlinear least squares, in which case you could optimize on an objective function that only searches over $\theta$, but inside this function, you can solve for $\beta$ and $\sigma$ linearly, given $\theta$. – GeoMatt22 Sep 7 '16 at 4:19
• @GeoMatt22 this is a non linear unconstrained optimziation problem, the estimate of $\beta$ is indeed a least squares estimator . Can you please provide some more details on how i can solve for $\beta$ given $\theta$ and how is that linear ? Thank you – Wis Sep 7 '16 at 4:31
• If I understand your description, by for a given $\theta$, the covariance matrix $K_{\theta}$ is then known. Since the data $y$ is fixed, and the trend basis-function matrix $h$ is fixed, your equation 4 should then be equivalent to a linear least squares problem. – GeoMatt22 Sep 7 '16 at 4:38
• @GeoMatt22 so you are proposing to interchangeably optimize between the two functions. Also do you advise any optimization technique suitable for such blocked optimization – Wis Sep 7 '16 at 4:46
• Say you use whatever optimizer, and you are passing it an objective function $F[\theta,\beta]$. Typically this would be a function handle/pointer, associated to some coded up version of the actual math expression. The idea here is then to define a new function $f[\theta]$, which you pass to your optimizer. Inside the code for $f$, the first thing it does is solve for $\beta$. Then it computes $F$, which it returns, i.e. $f[\theta]=F[\theta,\beta [\theta]]$. So your optimizer only searches over $\theta$. – GeoMatt22 Sep 7 '16 at 4:55

I think that there are two possible solutions to adopt. The firstone is an EM-kind approach where you start with optimization with respect to $\theta$, than get an exact value for $\beta$ and so on. Another option is to estimate coefficients $\beta$ using an OLS estimation and then estimate $\theta$ using obtained residuals.

Due to complexity of the presented problem I think that the final choice of your kriging model and how to optimize it should come from results of cross-validation.

Also I recommend you to test whether Universal kriging model (one with trend parameters estimated) is better than Simple kriging (we normalize data by subtracting the mean value and then work in the case when trend part is zero everywhere). In some cases this choice is not obvious, as kriging with trend estimation provides broader class of models, but requires more parameters to be estimated - which can be a problem if sample size is small.