Matérn Functions

Go back to the Contents page.

---

Press Show to reveal the code chunks.

---

# Set seed for reproducibility
set.seed(1982) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = FALSE,    
  # Disable warnings printed by R code chunks
  warning = FALSE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = FALSE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}

# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(rSPDE)
library(MetricGraph)
library(grateful)
library(Matrix)

library(ggplot2)
library(reshape2)
library(plotly)

The Matérn covariance function is given by

\[ \varrho(h; \kappa, \nu, \sigma)=\dfrac{\sigma^{2}}{2^{\nu-1} \Gamma(\nu)}(\kappa|h|)^\nu K_\nu(\kappa|h|),\quad \sigma^2 = \dfrac{\Gamma(\nu)}{\Gamma(\nu+1/2)(4\pi)^{1/2}\kappa^{2\nu}\tau^2} \tag{1} \label{materncovariance} \]

and it is coded below as matern.covariance().

matern.covariance <- function(h, kappa, nu, sigma) {
    if (nu == 1 / 2) {
        C <- sigma^2 * exp(-kappa * abs(h))
    } else {
        C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
            ((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
        C[h == 0] <- sigma^2
    }
    return(C)
}

The derivatives of the Matérn covariance function can be computed using the following recursive relation. For example, the first derivative is given by

\[ \varrho'(h; \kappa, \nu, \sigma) = \begin{cases} 0, & h = 0 \\ -\dfrac{\kappa^2}{2(\nu - 1)} h \varrho(h; \kappa, \nu - 1, \sigma), & h \neq 0 \end{cases} \tag{2} \label{matern_derivative} \]

and it can be obtained by setting deriv = 1 in the function matern.derivative() below.

matern.derivative <- function(h, kappa, nu, sigma, deriv = 1) 
{
    if(deriv == 0) {
        C = matern.covariance(h, kappa = kappa, nu = nu, sigma = sigma)
        return(C)
    } else if (deriv == 1) {
        C = h * matern.covariance(h, kappa = kappa, nu = nu - 1, sigma = sigma)
        C[h == 0] = 0
        return(-(kappa^2/(2 * (nu - 1))) * C)
    }
    else if (deriv == 2) {
        C = matern.covariance(h, kappa = kappa, nu = nu - 1, sigma = sigma) + 
            h * matern.derivative(h, kappa = kappa, nu = nu - 1, 
                                  sigma = sigma, deriv = 1)
        return(-(kappa^2/(2 * (nu - 1))) * C)
    }
    else {
        C = (deriv - 1) * matern.derivative(h, kappa = kappa, 
                                            nu = nu - 1, sigma = sigma, 
                                            deriv = deriv - 2) + 
            h * matern.derivative(h, kappa = kappa, nu = nu - 
                                      1, sigma = sigma, deriv = deriv - 1)
        return(-(kappa^2/(2 * (nu - 1))) * C)
    }
    
}

From (Bolin, Mehandiratta, and Simas 2025, Proposition 1), we have the following decomposition

\[ \varrho_m^\alpha(h)=\varrho_{m, 0}^\alpha(h)+\sum_{i=1}^m \varrho_{m, i}^\alpha(h) \tag{3} \label{matern_p_decomposition} \]

where

\[ \varrho_{m, 0}^\alpha(h)=k \sigma^2 \cdot \begin{cases}\dfrac{c_\alpha \sqrt{4 \pi}}{\kappa} 1_{[h=0]}, & 0<\alpha<1 \\ \varrho\left(h ;\kappa,\lfloor\alpha\rfloor-\frac{1}{2}, \sqrt{ \frac{c_\alpha}{c_{\lfloor\alpha\rfloor}}}\right), & \alpha \geq 1\end{cases} \tag{4} \label{matern_0_covariance} \]

and

\[ \varrho_{m, i}^\alpha(h)=r_i \sigma^2 \cdot \begin{cases}\varrho\left(h ;\kappa_i, \frac{1}{2}, \sqrt{ \frac{c_\alpha \sqrt{\pi}}{\sqrt{1-p_i}}}\right), & 0<\alpha<1 \\ \frac{1}{p_i^{\lfloor\alpha\rfloor}} \varrho\left(h ; \kappa_i, \frac{1}{2}, \sqrt{\frac{c_\alpha \sqrt{\pi}}{\sqrt{1-p_i}}}\right)-\sum_{j=1}^{\lfloor\alpha\rfloor} \frac{1}{p_i^{\lfloor\alpha\rfloor+1-j}} \varrho\left(h ; \kappa, j-\frac{1}{2}, \sqrt{\frac{c_\alpha}{c_j}}\right), & \alpha \geq 1\end{cases} \tag{5} \label{matern_p_covariance} \]

and \(c_a:=\Gamma(a) / \Gamma(a-1 / 2), \kappa_i=\kappa \sqrt{1-p_i}\), and \(\varrho\) is the Matérn covariance \(\eqref{materncovariance}\).

Function matern.p() implements

\[ \dfrac{\varrho_{m, i}^\alpha(h)}{r_i\sigma^2} \] as given below. Function matern.p.deriv() implements the derivatives of this function.

Here are the details of what matern.p() does.

• If p==0, then matern.p() returns \(\varrho(h = s-t; \kappa = \kappa, \nu = \alpha-1/2, \sigma = 1) =\dfrac{1}{2^{\nu-1} \Gamma(\nu)}(\kappa|h|)^\nu K_\nu(\kappa|h|)\). So if you want the full matern, just multiply by \(\sigma^2\). That is, \(\sigma^2*\)matern.p(s,t,kappa,p = 0, alpha) =\(\varrho(h = s-t; \kappa = \kappa, \nu = \alpha-1/2, \sigma = \sigma)\).

• If p!=0, then

  • If \(0<\alpha<1\), then matern.p() returns \(\varrho(h = s-t; \kappa = \kappa\sqrt{1-p_i}, \nu = 1/2, \sigma = \sqrt{\frac{c_\alpha \sqrt{\pi}}{\sqrt{1-p_i}}}) = \frac{c_\alpha \sqrt{\pi}}{\sqrt{1-p_i}} \cdot \varrho(h = s-t; \kappa = \kappa\sqrt{1-p_i}, \nu = 1/2, \sigma=1)\)

The bottom line is

• \(\varrho_{m, i}^\alpha(h)\) = \(r_i\sigma^2\)matern.p(s,t,kappa,p = p, alpha).

• \(\dfrac{d^k}{dh^k}\varrho_{m, i}^\alpha(h)\) = \(r_i\sigma^2\)matern.p.deriv(s,t,kappa,p = p, alpha, deriv = k).

• \(\mathbf{r}_i\left(t_1, t_2\right)=r_i\sigma^2\)matern.p.joint(s,t,kappa,p,alpha).

• If p==0, then we are in the \(\alpha\in\mathbb{N}\) case so that \(\lfloor\alpha\rfloor=\alpha\) and hence \(c_\alpha/c_{\lfloor\alpha\rfloor} = 1\) and hence \(\varrho_{m, 0}^\alpha(h)=k \sigma^2\varrho(h;\kappa, \alpha-1/2,1)=k \sigma^2\)matern.p(s,t,kappa,p=0,alpha)

matern.p <- function(s,t,kappa,p,alpha){
    h <- s-t
    if(p==0){
        return(matern.covariance(h, kappa = kappa, nu = alpha - 1/2, sigma = 1))
    } else {
        ca <- gamma(alpha)/gamma(alpha-0.5)
        fa <- floor(alpha)
        kp <- kappa*sqrt(1-p)
        out <- matern.covariance(h, kappa = kp, nu = 1/2, 
                                 sigma = sqrt(ca*sqrt(pi)/sqrt(1-p)))
        if(alpha < 1) {
            return(out)
        } else {
            
            for(j in 1:fa) {
                out <- out - matern.covariance(h, kappa = kappa, nu = j-1/2, 
                                               sigma = sqrt(ca*gamma(j-0.5)/gamma(j)))/p^(1 - j)
            }
            out <- out/p^fa
            return(out)    
        }
    }
}

matern.p.deriv <- function(s,t,kappa,p,alpha,deriv = 0){
    h <- s-t
    if(deriv ==0){
        return(matern.p(s,t,kappa,p,alpha))
    } else {
        if(p==0){
            return(matern.derivative(h, kappa = kappa, nu = alpha - 1/2, 
                                     sigma = 1, deriv = deriv))
        } else {
            ca <- gamma(alpha)/gamma(alpha-0.5)
            fa <- floor(alpha)
            kp <- kappa*sqrt(1-p)
            out <- matern.derivative(h, kappa = kp, nu =  1/2, 
                                     sigma = sqrt(ca*sqrt(pi)/sqrt(1-p)), deriv = deriv)
            if(alpha < 1) {
                return(out)
            } else {
                out <- out/p^fa
                for(j in 1:fa) {
                    out <- out - matern.derivative(h, kappa = kappa, nu = j-1/2, 
                                                   sigma = sqrt(ca*gamma(j-0.5)/gamma(j)), 
                                                   deriv = deriv)/p^(fa + 1 - j)
                }
            }
            return(out)
        }    
    }
    
}

Function matern.p.joint() computes the joint covariance matrix of the process and its derivatives up to order floor(alpha) for the shifted Matérn covariance. That is, it computes

\[ \mathbf{r}_i: \mathbb{R} \times \mathbb{R} \mapsto \mathbb{R}^{\lceil\alpha\rceil\times \lceil\alpha\rceil}, \quad \mathbf{r}_i\left(t_1, t_2\right)=\left[\frac{\partial^{i-1}}{\partial t_2^{i-1}} \frac{\partial^{j-1}}{\partial t_1^{j-1}} \dfrac{\varrho_{m, i}^\alpha\left(t_1-t_2\right)}{r_i\sigma^2}\right]_{i j \in\{1,2, \ldots, \lceil\alpha\rceil\}} \tag{2} \label{r_matrix} \] which is the multivariate covariance function of \((u_i(t), u_i'(t), \ldots, u_i^{(\lfloor\alpha\rfloor)}(t))\) and \(u_i\) is a shifted Whittle–Matérn bridge process with covariance function \(\varrho_{m, i}^\alpha(h)/(r_i\sigma^2)\).

#Joint covariance of process and derivative for shifted Matern
matern.p.joint <- function(s,t,kappa,p, alpha = 1){
    
    if(alpha%%1 == 0) {
        fa <- alpha
    } else {
        fa <- floor(alpha) + 1    
    }
    mat <- matrix(0, nrow = fa, ncol = fa)
    for(i in 1:fa) {
        for(j in i:fa) {
            if(i==j) {
                mat[i,i] <- ((-1)^(i-1))*matern.p.deriv(s,t,kappa,p,alpha, deriv = 2*(i-1))
            } else {
                tmp <- matern.p.deriv(s,t,kappa,p,alpha, deriv = i-1 + j - 1)
                mat[i,j] <- (-1)^(j-1)*tmp
                mat[j,i] <- (-1)^(i-1)*tmp
            }
        }
    }
    
    return(mat)
}

{r}
mat1 <- matern.p.joint(s = 0.6, t = 0.5, kappa = 10, p = 0.3, alpha = 2)
mat2 <- matern.p.joint(s = 0.6, t = 0.5, kappa = 10, p = 0.3, alpha = 3)
mat1
mat2

exp_precision <- function(loc, kappa, boundary = "free") {
    n <- length(loc)
    l <- diff(loc)
    
    # Precompute reusable terms
    a <- exp(-2 * kappa * l)
    b <- 1 - a
    
    # Initialize matrix Q efficiently
    Q <- Matrix(0, nrow = n, ncol = n)
    
    # Set diagonal elements
    diag(Q)[1:(n-1)] <- 1/2 + a/b
    diag(Q)[2:n] <- diag(Q)[2:n] + 1/2 + a/b
    
    # Set off-diagonal elements
    #  indices <- c(which(row(Q) == col(Q) + 1), which(row(Q) == col(Q) - 1))
    # off_diag_values <- -exp(-kappa * l) /b
    # Q[indices] <- off_diag_values
    
    
    row_indices <- seq_len(n - 1)
    col_indices <- row_indices + 1  # Offsets for off-diagonal elements
    
    # Compute the off-diagonal values
    off_diag_values <- -exp(-kappa * l)/b
    
    # Set the off-diagonal elements in matrix Q
    Q[cbind(row_indices, col_indices)] <- off_diag_values
    Q[cbind(col_indices, row_indices)] <- off_diag_values
    
    # Adjust boundary conditions if required
    if (boundary == "free") {
        Q[1, 1] <- Q[1, 1] + 0.5
        Q[n, n] <- Q[n, n] + 0.5
    }
    
    return(2 * kappa * Q[])
}

{r}
# Example usage of exp_precision function
loc <- seq(0, 1, length.out = 10)
kappa <- 5
Q_matrix <- exp_precision(loc, kappa, boundary = "free")
print(Q_matrix)

matern.p.precision <- function(loc,kappa,p,equally_spaced = FALSE, alpha = 1) {
    
    n <- length(loc)
    if(alpha==1) {
        Q <- exp_precision(loc,kappa)/(2*kappa)
        A <- Diagonal(n)
        
        return(list(Q=Q,A=A))
    } else {
        if(alpha%%1 == 0) {
            fa <- alpha
        } else {
            fa <- floor(alpha) + 1    
        }
        
        if(fa == 1) {
            N <- n  + n - 1 
        } else {
            N <- n*fa^2 + (n-1)*fa^2 - n*fa*(fa -1)/2    
        }
        
        ii <- numeric(N)
        jj <- numeric(N)
        val <- numeric(N)
        
        if(equally_spaced){
            Q.1 <- solve(rbind(cbind(matern.p.joint(loc[1],loc[1],kappa,p,alpha), 
                                     matern.p.joint(loc[1],loc[1+1],kappa,p,alpha)),
                               cbind(matern.p.joint(loc[1+1],loc[1],kappa,p,alpha), 
                                     matern.p.joint(loc[1+1],loc[1+1],kappa,p,alpha))))
            
            Q.i <- solve(rbind(cbind(matern.p.joint(loc[1],loc[1],kappa,p,alpha), 
                                     matern.p.joint(loc[1],loc[2],kappa,p,alpha),
                                     matern.p.joint(loc[1],loc[3],kappa,p,alpha)),
                               cbind(matern.p.joint(loc[2],loc[1],kappa,p,alpha),
                                     matern.p.joint(loc[2],loc[2],kappa,p,alpha),
                                     matern.p.joint(loc[2],loc[3],kappa,p,alpha)),
                               cbind(matern.p.joint(loc[3],loc[1],kappa,p,alpha),
                                     matern.p.joint(loc[3],loc[2],kappa,p,alpha),
                                     matern.p.joint(loc[3],loc[3],kappa,p,alpha))))[-c(1:fa),-c(1:fa)]
        }
        
        
        for(i in 1:max((n-1),1)){
            if(i==1){
                if(!equally_spaced){
                    Q.1 <- solve(rbind(cbind(matern.p.joint(loc[i],loc[i],kappa,p,alpha),
                                             matern.p.joint(loc[i],loc[i+1],kappa,p,alpha)),
                                       cbind(matern.p.joint(loc[i+1],loc[i],kappa,p,alpha),
                                             matern.p.joint(loc[i+1],loc[i+1],kappa,p,alpha))))
                    
                } 
                counter <- 1
                for(ki in 1:fa) {
                    for(kj in ki:(2*fa)) {
                        ii[counter] <- ki
                        jj[counter] <- kj
                        val[counter] <- Q.1[ki,kj]
                        counter <- counter + 1
                    }
                }
            } else {
                if(!equally_spaced){
                    Q.i <- solve(rbind(cbind(matern.p.joint(loc[i-1],loc[i-1],kappa,p,alpha),
                                             matern.p.joint(loc[i-1],loc[i],kappa,p,alpha),
                                             matern.p.joint(loc[i-1],loc[i+1],kappa,p,alpha)),
                                       cbind(matern.p.joint(loc[i],loc[i-1],kappa,p,alpha),
                                             matern.p.joint(loc[i],loc[i],kappa,p,alpha),
                                             matern.p.joint(loc[i],loc[i+1],kappa,p,alpha)),
                                       cbind(matern.p.joint(loc[i+1],loc[i-1],kappa,p,alpha),
                                             matern.p.joint(loc[i+1],loc[i],kappa,p,alpha),
                                             matern.p.joint(loc[i+1],loc[i+1],kappa,p,alpha))))[-c(1:fa),-c(1:fa)]
                } 
                # Q[(2*n-1):(2*n), (2*n-3):(2*n)] = Q.i[3:4,]
                for(ki in 1:fa){
                    for(kj in ki:(2*fa)){
                        ii[counter] <- fa*(i-1) + ki
                        jj[counter] <- fa*(i-1) + kj
                        val[counter] <-Q.i[ki,kj]
                        counter <- counter + 1
                    }
                }     
            }
        }
        if(n<=2){
            Q.i <- Q.1
        }
        for(ki in 1:fa){
            for(kj in ki:fa){
                ii[counter] <- fa*(n-1) + ki
                jj[counter] <- fa*(n-1) + kj
                val[counter] <-Q.i[ki+fa,kj+fa]
                counter <- counter + 1
            }
        }    
        Q <- Matrix::sparseMatrix(i   = ii,
                                  j    = jj,
                                  x    = val,
                                  dims = c(fa*n, fa*n), symmetric=TRUE)
        
        A <-  Matrix::sparseMatrix(i   = 1:n,
                                   j    = seq(from=1,to=n*fa,by=fa),
                                   x    = rep(1,n),
                                   dims = c(n, fa*n))
        return(list(Q=Q,A=A))    
    }
}

{r}
aux <-  matern.p.precision(loc = c(0,1,2), kappa = 10, p = 0.3, equally_spaced = TRUE, alpha = 2.2)
aux$Q
aux$A

{r}
# ------------ DATA ----------------------------------------------------------
kappa <- 25
nu <- 3/2
alpha <- nu + 1/2
sigma <- 1
p <- 0.3
x <- seq(0, 1, length.out = 1001)
h <- abs(x - 0.5)

y1 <- matern.covariance(h, kappa = kappa, nu = nu, sigma = sigma)
y2 <- matern.derivative(h, kappa = kappa, nu = nu, sigma = sigma, deriv = 1)
y3 <- matern.derivative(h, kappa = kappa, nu = nu, sigma = sigma, deriv = 2)

y4 <- matern.p(s = x, t = 0.5, kappa = kappa, p = 0, alpha = alpha)
y5 <- matern.p(s = x, t = 0.5, kappa = kappa, p = p, alpha = alpha)
y6 <- sapply(x, function(si) matern.p.joint(s = si, t = 0.5, kappa = kappa, p = p, alpha = alpha)[1,1])
y7 <- matern.p.deriv(s = x, t = 0.5, kappa = kappa, p = p, alpha = alpha, deriv = 1)
y8 <- sapply(x, function(si) matern.p.joint(s = si, t = 0.5, kappa = kappa, p = p, alpha = alpha)[2,2])
# ------------ SCALING -------------------------------------------------------

s1 <- max(abs(y1))
s2 <- max(abs(y2))
s3 <- max(abs(y3))
s4 <- max(abs(y4))
s5 <- max(abs(y5))
s6 <- max(abs(y6))
s7 <- max(abs(y7))
s8 <- max(abs(y8))

z1 <- y1 / s1
z2 <- y2 / s2
z3 <- y3 / s3
z4 <- y4 / s4
z5 <- y5 / s5
z6 <- y6 / s6
z7 <- y7 / s7
z8 <- y8 / s8

# ------------ PLOT ----------------------------------------------------------

plot_ly() %>%
  add_trace(
    x = ~x, y = ~0, z = ~z1,
    type = "scatter3d", mode = "lines",
    name = paste0("Matérn (1/", round(s1, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~1, z = ~z2,
    type = "scatter3d", mode = "lines",
    name = paste0("Matérn' (1/", round(s2, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~2, z = ~z3,
    type = "scatter3d", mode = "lines",
    name = paste0("Matérn'' (1/", round(s3, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~3, z = ~z4,
    type = "scatter3d", mode = "lines",
    name = paste0("matern.0 (1/", round(s4, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~4, z = ~z5,
    type = "scatter3d", mode = "lines",
    name = paste0("matern.p (1/", round(s5, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~5, z = ~z6,
    type = "scatter3d", mode = "lines",
    name = paste0("matern.p joint (1/", round(s6, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~6, z = ~z7,
    type = "scatter3d", mode = "lines",
    name = paste0("matern.p' (1/", round(s7, 3), ")")
  ) %>%
  add_trace(
    x = ~x, y = ~7, z = ~z8,
    type = "scatter3d", mode = "lines",
    name = paste0("matern.p joint' (1/", round(s8, 3), ")")
  ) %>%
  layout(
    scene = list(
      xaxis = list(title = "x"),
      yaxis = list(title = "Curve index"),
      zaxis = list(title = "Scaled value [-1, 1]")
    ),
    legend = list(title = list(text = "<b>Functions</b>"))
  )

Bolin, David, Vaibhav Mehandiratta, and Alexandre B. Simas. 2025. “Linear Cost and Exponentially Convergent Approximation of Gaussian Matérn Processes on Intervals.” https://arxiv.org/abs/2410.13000.