# Hyperparameter tuning in Gaussian Process Regression

I am trying to tune the hyperparameters of the gaussian process regression algorithm I've implemented. I simply want to maximize the log marginal likelihood given by the formula $$\log(\mathbf{y}|X,\mathbf{\theta})=-\frac{1}{2} \mathbf{y}^TK_y^{-1}\mathbf{y}-\frac{1}{2}\log(\det(K))-\frac{n}{2}\log(2\pi)$$ where $K$ is the covariance matrix with the elements $$K_{ij}=k(x_i,x_j)=b^{-1}\exp(-\frac{1}{2}(x_i-x_j)^TM(x_i-x_j))+a^{-1}\delta_{ij}$$ where $M=lI$ and $a,b$ and $l$ are hyperparameters.

partial derivative of the log marginal likelihood w.r.t parameters is given by the following $$\frac{\log(\mathbf{y}|X,\mathbf{\theta})}{d\theta}=\frac{1}{2}\mathrm{trace}(K^{-1}\frac{dK}{d\theta})+\frac{1}{2}(\mathbf{y}\frac{dK}{d\theta}K^{-1}\frac{dK}{d\theta}\mathbf{y})$$

As the entries of $K$ depend on the parameters, so do derivatives and inverse of $K$. This means, when a gradient-based optimizer is employed, evaluating the gradient at a given point (parameter value) will require recomputation of the covariance matrix. In my application, this is not feasible because computing the covariance matrix from scratch and computing its inverse in every iteration of gradient ascent is too expensive. My question is what my options are to find a fairly good combination of these three parameters? and I also don't know which parameter to optimize first and I would appreciate any pointers on this issue as well.

• may i recommend local gaussian processes? cran.r-project.org/web/packages/laGP/index.html Feb 21 at 23:10
• otherwise, if youre comfortable with math, what's hot in the field nowadays are some kinda krylov subspace approach to linear system solution / determinant calculation such as arxiv.org/pdf/1803.06058.pdf Feb 21 at 23:11

You are right that you need a new covariance matrix computation on every iteration of gradient ascent. So if the matrix computation is not feasible for your setting, then, I think, you cannot use gradient-based marginal likelihood optimization.

My suggestion is to use gradient-free methods for hyperparameter tuning, such as grid search, random search, or Bayesian optimization-based search. These methods are widely used for optimization hyperparameters of other machine learning algorithms e.g. SVMs.

I suggest the grid search for your first try. You basically form a table (grid) of possible hyperparameters, try every one, and look for the best validation performance (or best marginal likelihood).

Grid search would yield a suboptimal set of hyperparameters, and you have to specify grid by yourself.(tip: make grid in a log scale) but far less computation is needed. (and you don't need gradient!)

If you are not familiar with grid search, you can look up Wikipedia:Hyperparameter Optimization - Grid Search

If solving the linear problem $$K\pmb{\alpha} = \textbf{y}$$ is too expensive for you at each step of your optimisation, you could resort to approximation techniques such as the Nystr$$\ddot{o}$$m method (e.g., Williams and Seeger 2001), such that if $$K \in \mathbb{R}^{n \times n}$$, the Nystr$$\ddot{o}$$m method computes a rank $$m$$ approximation of $$K$$, where $$m < n$$. Standard inference in the GP setting has run time complexity $$\mathcal{O}(n^3)$$, while in the Nystr$$\ddot{o}$$m approach it is $$\mathcal{O}(nm^2)$$.

In a bit of detail:

First, select $$m$$ random samples from your inputs, $$X_{mm}$$, which are then referred to as the active set. Next, compute the relevant (noise free) submatrix $$K_{mm} = k(X_{mm},X_{mm})$$. The eigenvectors and eigenvalues of $$K_{mm}$$, denoted $$\textbf{u}$$ and $$\lambda$$ respectively, can then be extended to all $$n$$ data points by the following:

$$\tilde{\lambda} = \dfrac{n}{m}\lambda \\ \tilde{\textbf{u}} = \sqrt{\dfrac{m}{n}}\lambda^{-1}K_{nm}\textbf{u}$$

where $$K_{nm} = k(X,X_{mm})$$.

Notice now how the reduced-rank approximation of $$K$$ is given by:

$$\tilde{K} = (\tilde{\lambda}\tilde{\textbf{u}})\tilde{\textbf{u}}^\text{T}$$

where in the case $$m = n$$, $$\tilde{K} = K$$.

Now, computing the inverse of $$\tilde{K} + \pmb{V}$$, where $$\pmb{V} = a\pmb{I}$$, simply requires use of the Woodbury identity:

$$(\pmb{V} + (\tilde{\lambda}\tilde{\textbf{u}})\tilde{\textbf{u}}^\text{T})^{-1} = \pmb{V}^{-1} - \pmb{V}^{-1}\tilde{\textbf{u}}\big(diag(\tilde{\lambda}^{-1}) + \tilde{\textbf{u}}^\text{T}\pmb{V}^{-1}\tilde{\textbf{u}}\big)^{-1}\tilde{\textbf{u}}^\text{T}\pmb{V}^{-1}$$

which is now much cheaper, because $$\big(diag(\tilde{\lambda}^{-1}) + \tilde{\textbf{u}}^\text{T}\pmb{V}^{-1}\tilde{\textbf{u}}\big)^{-1}$$ is of size $$m \times m$$.

Finally, the evaluation of $$det(\tilde{K} + \pmb{I}a)$$ (required for computing the log marginal likelihood), can also be derived cheaply. First let us define $$\tilde{L} = \sqrt{\tilde{\lambda}}\tilde{\textbf{u}}$$, such that $$\tilde{L}\tilde{L}^\text{T} = \tilde{K}$$. Now, we can use Sylvester's determinant theorem, which states that

$$det(\pmb{I}_na + \tilde{L}\tilde{L}^\text{T}) = det(\pmb{I}_ma + \tilde{L}^\text{T}\tilde{L})$$

The LHS of the above definition has complexity $$\mathcal{O}(n^3)$$, while on the RHS it is $$\mathcal{O}(m^3)$$.

• I would add as a comment here that in practice it may be necessary to set any negative values in $\lambda$ to something small like 1e-12. Furthermore, for the computation of the determinant I would recommend something like numpy.linalg.slogdet if using Python, as this handles very small (or large) determinants better, and the log of the determinant is anyway required for the log marginal likelihood Dec 10, 2021 at 13:51