# Beta Coefficient for Gaussian Process with Non-Zero Mean

I'm reading through Rasmussen & Williams (2006)'s book on Gaussian Processes, specifically on section 2.7 on incorporating explicit basis function.

I am confused on how the analytical solution $$\beta$$ was derived in equation 2.41 for the conditional mean

$$\mu = H_\star^\top \beta + K_\star^\top K_y^{-1}(y - H^\top \beta)$$

where

$$\beta = (B^{-1} + HK_y^{-1}H^\top)^{-1}(HK_{y}^{-1}y + B^{-1}b)$$

I understand that in general, if $$y_1$$ and $$y_2$$ are jointly Gaussian, $$\begin{bmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \end{bmatrix} \sim \mathcal{N}\Big( \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{bmatrix} \Big)$$ then we have that the conditional mean is given by $$\Pr(\mathbf{y}_2 \big\vert \mathbf{y}_1, \mathbf{X}_1, \mathbf{X}_2) \sim \mathcal{N}(\mu_{2\vert 1}, \Sigma_{2\vert 1}) \implies \mu_{2\vert 1} = \mu_2 + \Sigma_{21}\Sigma_{11}^{-1}(y_1 - \mu_1)$$

This is similar to the one given by the Gaussian process posterior because we are assuming $$y_2 = g(X_2)$$ where $$\begin{gather} g(x) = h(x)^\top\beta + f(x) \\ \beta \sim \mathcal{N}(b, B) \\ f(x) \sim \mathcal{GP}(0, k(x, x')) \\ \end{gather}$$

I'm interested in how $$\beta$$ was analytically derived because I am trying to implement it in python and explicitly inverting $$K_y^{-1}$$ requires a lot of memory, while solving the linear system $$f(K_y)\beta = \text{something}$$ may be more memory efficient with scipy.

Edit 1: I realized, because of the prior distribution of $$\beta$$, this is close to the idea and closed form solution of Tikhonov_regularization.

Edit 2: Ah yes, it is the same. For a regularization regression problem $$\begin{gather} y = w^\top X + \epsilon \\ \epsilon \sim \mathcal{N}(0, \Pi) \\ w \sim \mathcal{N}(\mu, \Sigma) \\ \end{gather}$$ the optimal weight is given by $$w^\star = (\Sigma^{-1} + X^\top \Pi^{-1}X)^{-1}(\Sigma^{-1}\mu + X^\top \Pi^{-1}y)$$ This makes sense because in applying Gaussian Processes to a regression problem, where the residuals / errors are distributed jointly Gaussian, these two problems are equal (?). Still unknown whether it's possible to obtain $$\beta$$ without inverting the Kernel matrix.