Skip to contents

Derivatives for estimating Gaussian process parameters

Collin Erickson

2026-05-11

Source: vignettes/derivatives.Rmd
derivatives.Rmd

Estimating the parameters is often done through maximum likelihood estimation. In this vignette, derivatives of parameters are calculated to give the parameter gradient. This can then be used in numerical optimization, such as with BFGS, and should be much faster than calculating without having the gradient.

Deviance

In this vignette, I follow the lead of the GPfit paper. Instead of maximizing the likelihood, we find parameter estimates by minimizing the deviance. For correlation parameters θ\theta, the deviance is shown below.

−2logLθ∝log|R|+nlog[(Y−1nμ̂(θ))TR−1(Y−1nμ̂(θ))]=𝒟 -2\log L_{\theta} \propto \log |R| + n\log[(Y-1_n \hat{\mu}(\theta))^T R^{-1} (Y-1_n \hat{\mu}(\theta))] = \mathcal{D}

For now I will assume that μ̂\hat{\mu} does not depend on θ\theta, and will replace Y−μ̂(θ)Y-\hat{\mu}(\theta) with ZZ, as shown below.

𝒟=log|R|+nlog[ZTR−1Z] \mathcal{D} = \log |R| + n\log[Z^T R^{-1} Z]

Thus the only dependence on the correlation parameters is through RR.

Via the GPML book section A.3.1, we need the following equations:

∂∂θR−1=−R−1∂R∂θR−1 \frac{\partial }{\partial \theta} R^{-1} = -R^{-1} \frac{\partial R}{\partial \theta} R^{-1}

∂∂θlog|R|=tr(R−1∂R∂θ) \frac{\partial }{\partial \theta} \log |R| = \text{tr}(R ^ {-1}\frac{\partial R}{\partial \theta} )

Now we can calculate the derivative of the deviance in terms of ∂R∂θ\frac{\partial R}{\partial \theta}.

∂𝒟∂θ=tr(R−1∂R∂θ)−nZTR−1ZZTR−1∂R∂θR−1Z \frac{\partial \mathcal{D}}{\partial \theta} = \text{tr}(R ^ {-1}\frac{\partial R}{\partial \theta} ) - \frac{n}{Z^T R^{-1} Z} Z^T R^{-1} \frac{\partial R}{\partial \theta} R^{-1} Z

Now we just need to find ∂R∂θ\frac{\partial R}{\partial \theta}, which depends on the specific correlation function R(θ)R(\theta).

Nugget

The correlation function is usually augmented by adding a nugget term δ\delta to the diagonal:

R=R*+δI R = R^* + \delta I

The nugget accounts for noise in the responses, smoothing the predicted response function. It also helps with numerical stability, which is a serious problem when there is a lot of data. Often a small value is used even the function is noiseless.

For the nugget,

∂R∂δ=I\frac{\partial R}{\partial \delta} = I

Thus, the derivative for the deviance is very simple. ∂𝒟∂δ=tr(R−1)−nZTR−1ZZTR−1R−1Z \frac{\partial \mathcal{D}}{\partial \delta} = \text{tr}(R ^ {-1} ) - \frac{n}{Z^T R^{-1} Z} Z^T R^{-1} R^{-1} Z

This equation gives the derivative of the deviance with respect to the nugget regardless of the correlation function used.

Gaussian correlation

The Gaussian correlation function has parameter vector θ=(θ1,...,θd)\theta = (\theta_1, ..., \theta_d).

Ignoring the nugget term, the ii,jj entry of the correlation matrix is

Rij(x,y)=exp[−∑k=1dθk(xik−xjk)2] R_{ij}(x, y) = \exp\bigg[-\sum_{k=1}^{d} \theta_k (x_{ik} - x_{jk})^2 \bigg]

∂∂θlRij(x,y)=−exp[−∑k=1dθk(xik−xjk)2](xil−xjl)2 \frac{\partial}{\partial \theta_l} R_{ij}(x, y) = -\exp\bigg[-\sum_{k=1}^{d} \theta_k (x_{ik} - x_{jk})^2 \bigg] (x_{il} - x_{jl})^2

This will give the matrix ∂R∂θl\frac{\partial R}{\partial \theta_l}, which can be used with the previous equations to calculate ∂𝒟∂θl\frac{\partial \mathcal{D}}{\partial \theta_l}.

Lifted brownian covariance

The lifted brownian covariance function is c(x,x′)=ψ(x)+ψ(x′)+ψ(x−x′)−1 c(x, x') = \psi(x) + \psi(x') + \psi(x-x') - 1 where ψ(h)=(1+||h||a2γ)β \psi(h) = (1 + ||h||_a^{2\gamma} ) ^{\beta}

and ||h||a2=𝐡T[a10…00a2…0⋮⋮⋱⋮00…ad]𝐡=∑i=1daihi2 ||h||_a^2 = \boldsymbol{h}^T \begin{bmatrix} a_1 & 0 & \dots & 0 \\ 0 & a_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & a_d \end{bmatrix} \boldsymbol{h} = \sum_{i=1}^d a_i h_i^2

This is different from the others because it is not a correlation function, which ranges from 0 to 1, but is the covariance function itself. Thus we will have to use first deviance before inserting σ̂2\hat{\sigma}^2.

∂∂βψ(h)=(1+||h||a2γ)βlog(1+||h||a2γ)=ψ(h)log(1+||h||a2γ) \frac{\partial}{\partial \beta} \psi(h) = (1 + ||h||_a^{2\gamma} ) ^{\beta} \log (1 + ||h||_a^{2\gamma} ) = \psi(h) \log (1 + ||h||_a^{2\gamma} )

∂∂γψ(h)=β(1+||h||a2γ)β−1||h||a2γlog(||h||a2) \frac{\partial}{\partial \gamma} \psi(h) = \beta (1 + ||h||_a^{2\gamma} ) ^{\beta-1} ||h||_a^{2\gamma} \log(||h||_a^{2}) ∂∂ai||h||a2=hi2 \frac{\partial}{\partial a_i} ||h||_a^{2} = h_i^2

∂∂ai||h||a2γ=γ||h||a2(γ−1)hi2 \frac{\partial}{\partial a_i} ||h||_a^{2\gamma} = \gamma ||h||_a^{2(\gamma-1)} h_i^2

∂∂aiψ(h)=β(1+||h||a2γ)β−1γ||h||a2(γ−1)hi2 \frac{\partial}{\partial a_i} \psi(h) = \beta (1 + ||h||_a^{2\gamma} ) ^{\beta-1} \gamma ||h||_a^{2(\gamma-1)} h_i^2

Likelihood

The likelihood for data from a Gaussian process follows the standard multivariate normal probability distribution function (pdf). L=(2π)−k/2|Σ|−1/2exp[−12(Y−μ)Σ−1(Y−μ)] L = (2 \pi)^{-k/2} |\Sigma|^{-1/2} \exp[\frac{-1}{2}(Y - \mu) \Sigma^{-1} (Y - \mu)] The log likelihood is generally easier to work with.

logL=−k2log(2π)+−12log|Σ|+−12(Y−μ)Σ−1(Y−μ) \log L = \frac{-k}{2} \log(2 \pi) + \frac{-1}{2}\log|\Sigma| + \frac{-1}{2}(Y - \mu) \Sigma^{-1} (Y - \mu) To simplify, we can multiply it by -2, and call this the deviance, denoted here as 𝒟\mathcal{D}.

𝒟=−2logL=klog(2π)+log|Σ|+(Y−μ)Σ−1(Y−μ) \mathcal{D} = -2\log L = k \log(2 \pi) + \log|\Sigma| + (Y - \mu) \Sigma^{-1} (Y - \mu) kk can be ignored since it usually constant while optimizing parameters. 𝒟=−2logL∝log|Σ|+(Y−μ)Σ−1(Y−μ) \mathcal{D} = -2\log L \propto \log|\Sigma| + (Y - \mu) \Sigma^{-1} (Y - \mu)

∂𝒟∂θ=tr(Σ−1∂Σ∂θ)−ZTΣ−1∂Σ∂θΣ−1Z \frac{\partial \mathcal{D}}{\partial \theta} = \text{tr}(\Sigma ^ {-1}\frac{\partial \Sigma}{\partial \theta} ) - Z^T \Sigma^{-1} \frac{\partial \Sigma}{\partial \theta} \Sigma^{-1} Z