Convergence rates for the circle graph

Go back to the Convergence Rates page.

Let us set some global options for all code chunks in this document.

# 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 = TRUE,       
  # 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)
}

Import libraries

library(rSPDE)
library(MetricGraph)
library(Matrix)

library(dplyr)
library(plotly)
library(scales)
library(patchwork)
library(tidyr)
library(ggplot2)
library(reshape2)
library(sf)

library(here)
library(rmarkdown)
library(knitr)
library(grateful) # Cite all loaded packages

library(latex2exp)
library(plotrix)

Circle graph

Define utility functions

# Matern covariance function
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(as.matrix(C))
}

# Folded.matern.covariance.1d
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma, L = 1, N = 10, boundary = c("neumann",
                                                                                               "dirichlet", "periodic")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
    stop("The possible boundary conditions are 'neumann',
    'dirichlet' or 'periodic'!")
  }
  addi = t(outer(x, x, "+"))
  diff = t(outer(x, x, "-"))
  s1 <- sapply(-N:N, function(j) { # s1 is a matrix of size length(h)x(2N+1)
    diff + 2 * j * L
  })
  s2 <- sapply(-N:N, function(j) {
    addi + 2 * j * L
  })
  if (boundary == "neumann") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) +
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else if (boundary == "dirichlet") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) -
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else {
    C <- rowSums(matern.covariance(h = s1,
                                   kappa = kappa, nu = nu, sigma = sigma))
  }
  return(matrix(C, nrow = length(x)))
}

# Function to get the true covariance matrix
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
  h = c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[,2])
  true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
  return(true_cov_mat)
}

Define the graph

r <- 1/(pi)
theta <- seq(from = -pi, to = pi, length.out = 100)
edge <- cbind(1+r+r*cos(theta), r*sin(theta))
edges <- list(edge)
graph <- metric_graph$new(edges = edges)
graph$plot() + 
  ggtitle("Circle graph") + 
  theme_minimal() + 
  theme(text = element_text(family = "Palatino")) +
  coord_fixed(ratio = 1)

Figure 1: Circle graph.

Define the parameters

# parameters
h.ok <- 2^-10
type <- "covariance"
type_rational_approximation = "brasil"
rho <- 0.5
#m = 4
sigma <- 1
N.folded <- 10
boundary <- "periodic" # Do not change this

# Mesh sizes
h_aux <- seq(5.5, 4.5, by = -1/4)
h_vector <- 2^-h_aux
h_label <- paste0("2^-", h_aux, "")
h_label_latex <- sprintf("$2^{-%f}$", h_aux)

# Alpha values
alpha_aux <- c(6, 7, 8, 9, 12)
alpha_vector <- alpha_aux/8
alpha_label <- paste0(alpha_aux, "/8")
theoretical_rate <- pmin(2*alpha_vector-1/2,2)

Build the graph with overkill mesh

graph.ok <- graph$clone()
# Build graph with overkill mesh
graph.ok$build_mesh(h = h.ok)
graph.ok$compute_fem()

graph.ok$plot(mesh = TRUE) +
  ggtitle("Circle graph with overkill mesh") +
  theme_minimal() +
  theme(text = element_text(family = "Palatino")) +
  coord_fixed(ratio = 1)

Figure 2: Circle graph with overkill mesh.

# Get the overkill mesh locations
loc.ok <- graph.ok$mesh$VtE # or graph.ok$get_mesh_locations()

Each A[[i]] below is the projection matrix corresponding to the graph with mesh size h_vector[i] onto the overkill mesh.

# Initialize the list of graphs and the list of projection matrices
graphs <- list()
A <- list()
for(i in 1:length(h_vector)){
  graphs[[i]] <- graph$clone()
  graphs[[i]]$build_mesh(h = h_vector[i])
  A[[i]] <- graphs[[i]]$fem_basis(loc.ok)
}
# Print the dimensions of the projection matrices
print(lapply(A, dim))

## [[1]]
## [1] 2048   91
## 
## [[2]]
## [1] 2048   77
## 
## [[3]]
## [1] 2048   64
## 
## [[4]]
## [1] 2048   54
## 
## [[5]]
## [1] 2048   46

Compute the covariance error

cov.error.ok.mesh <- matrix(NA, nrow = length(h_vector), ncol = length(alpha_vector))
cov.error.folded <- matrix(NA, nrow = length(h_vector), ncol = length(alpha_vector))

par(family = "Palatino")
layout(matrix(1:(length(alpha_vector)*length(h_vector)),
              nrow = length(alpha_vector),
              byrow = TRUE))

m_values <- c()
for (j in 1:length(alpha_vector)) {
  alpha <- alpha_vector[j]
  fract_alpha <- alpha - floor(alpha)
  nu <- alpha - 0.5
  kappa <- sqrt(8*nu)/rho
  tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1

  Sigma.folded <- gets_true_cov_mat(graph = graph.ok, # using folded
                               kappa = kappa,
                               nu = nu,
                               sigma = sigma,
                               N = N.folded,
                               boundary = boundary)

  for (i in 1:length(h_vector)) {
    h <- h_vector[i]
    m <- min(20, 5*ceiling((min(2*alpha - 1/2,2) + 1/2)^2*log(h)^2/(4*pi^2*fract_alpha)))
    m_values <- c(m_values, m)

    Sigma.ok.mesh <- matern.operators(alpha = alpha, # using overkill mesh
                                kappa = kappa,
                                tau = tau,
                                m = m,
                                graph = graph.ok,
                                type = type,
                                type_rational_approximation = type_rational_approximation)$covariance_mesh()

    Sigma <- matern.operators(alpha = alpha,  # this is the approximation
                             kappa = kappa,
                             tau = tau,
                             m = m,
                             graph = graphs[[i]],
                             type = type,
                             type_rational_approximation = type_rational_approximation)$covariance_mesh()

    Sigma.approx <- A[[i]]%*%Sigma%*%t(A[[i]])

    cov.error.ok.mesh[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.ok.mesh - Sigma.approx)^2%*%graph.ok$mesh$weights))
    cov.error.folded[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.folded - Sigma.approx)^2%*%graph.ok$mesh$weights))

    aux <- dim(Sigma.approx)[1]
    reordering <- c(1,3:aux,2)
    auxhalf <- reordering[length(reordering) %/% 2 + 1]

    plot(Sigma.folded[auxhalf, reordering], col = "black", type = "l", main = bquote(alpha == .(alpha_label[j]) ~ "," ~ h == .(h_label[i])), xlab = "", ylab = "", lwd = 2, lty = 1)
    lines(Sigma.ok.mesh[auxhalf, reordering], col = "red", lty = 2, lwd = 2)
    lines(Sigma.approx[auxhalf, reordering], col = "blue", lty = 3, lwd = 2)
    legend("topleft", lwd = 2, col = c("black", "red", "blue"), legend = c("true", "ok.mesh", "approx"), lty = c(1,2,3))
  }
}

Figure 3: Covariance error for the circle graph.

print(m_values)

##  [1] 10 10  5  5  5 10 10 10  5  5 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20

slope_circle_folded <- numeric(length(alpha_vector))
slope_circle_ok_mesh <- numeric(length(alpha_vector))

# folded
for (u in 1:length(alpha_vector)) {
  slope_circle_folded[u] <- coef(lm(log(cov.error.folded[,u]) ~ log(h_vector)))[2]
  }
# ok.mesh
for (u in 1:length(alpha_vector)) {
  slope_circle_ok_mesh[u] <- coef(lm(log(cov.error.ok.mesh[,u]) ~ log(h_vector)))[2]
  }

transposed_df <- data.frame(t(data.frame(alpha = alpha_vector, theoretical_rate = theoretical_rate, slope1 = slope_circle_folded, slope2 = slope_circle_ok_mesh)))
rownames(transposed_df) <- c("alpha", "Theoretical rates", "Folded", "Overkill mesh")
colnames(transposed_df) <- NULL
# Display the transposed data frame
transposed_df |> paged_table()

---

Press the Show button below to reveal the code.

---

loglog_line_equation <- function(x1, y1, slope) {
  b <- log10(y1 / (x1 ^ slope))

  function(x) {
    (x ^ slope) * (10 ^ b)
  }
}

guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(alpha_vector))
for (j in 1:length(alpha_vector)) {
  guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
  for(k in 1:length(h_vector)){
    point_x1 <- h_vector[k]
    point_y1 <- cov.error.ok.mesh[k, j]
    slope <- theoretical_rate[j]
    line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
    guiding_lines_aux[,k] <- line(h_vector)
  }
  guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}

# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))

# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
  geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
            aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})

df <- as.data.frame(cbind(h_vector, cov.error.ok.mesh))
colnames(df) <- c("h_vector", alpha_vector)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")

p <- ggplot() +
  geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  plot_lines +
  labs(title = "Covariance error using overkill mesh",
    x = "h",
       y = "Covariance Error",
       color = expression(alpha)) +
  scale_x_log10(breaks = h_vector, labels = round(h_vector,3)) +
  scale_y_log10() +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

---

Press the Show button below to reveal the code.

---

guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(alpha_vector))
for (j in 1:length(alpha_vector)) {
  guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
  for(k in 1:length(h_vector)){
    point_x1 <- h_vector[k]
    point_y1 <- cov.error.folded[k, j]
    slope <- theoretical_rate[j]
    line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
    guiding_lines_aux[,k] <- line(h_vector)
  }
  guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}

# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))

# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
  geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
            aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})

df <- as.data.frame(cbind(h_vector, cov.error.folded))
colnames(df) <- c("h_vector", alpha_vector)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")

q <- ggplot() +
  geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  plot_lines +
  labs(title = "Covariance error using folded",
    x = "h",
       y = "Covariance Error",
       color = expression(alpha)) +
  scale_x_log10(breaks = h_vector, labels = round(h_vector,3)) +
  scale_y_log10() +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

p + q

Figure 4: Covariance error for the circle graph, using overkill mesh and folded.

References

cite_packages(output = "paragraph", out.dir = ".")

We used R version 4.4.1 (R Core Team 2024) and the following R packages: here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), knitr v. 1.48 (Xie 2014, 2015, 2024), latex2exp v. 0.9.6 (Meschiari 2022), Matrix v. 1.6.5 (Bates, Maechler, and Jagan 2024), MetricGraph v. 1.4.0.9000 (Bolin, Simas, and Wallin 2023b, 2023a, 2023c, 2024; Bolin et al. 2024), patchwork v. 1.2.0 (Pedersen 2024), plotly v. 4.10.4 (Sievert 2020), plotrix v. 3.8.4 (J 2006), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.28 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.4.0.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), sf v. 1.0.19 (Pebesma 2018; Pebesma and Bivand 2023), tidyverse v. 2.0.0 (Wickham et al. 2019), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).

Aden-Buie, Garrick, and Matthew T. Warkentin. 2024. xaringanExtra: Extras and Extensions for “xaringan” Slides. https://CRAN.R-project.org/package=xaringanExtra.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2024. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” Journal of Computational and Graphical Statistics 29 (2): 274–85. https://doi.org/10.1080/10618600.2019.1665537.

Bolin, David, Mihály Kovács, Vivek Kumar, and Alexandre B. Simas. 2024. “Regularity and Numerical Approximation of Fractional Elliptic Differential Equations on Compact Metric Graphs.” Mathematics of Computation 93 (349): 2439–72. https://doi.org/10.1090/mcom/3929.

Bolin, David, and Alexandre B. Simas. 2023. rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations. https://CRAN.R-project.org/package=rSPDE.

Bolin, David, Alexandre B. Simas, and Jonas Wallin. 2023a. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” arXiv Preprint arXiv:2304.03190. https://doi.org/10.48550/arXiv.2304.03190.

———. 2023b. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023c. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” arXiv Preprint arXiv:2304.10372. https://doi.org/10.48550/arXiv.2304.10372.

———. 2024. “Gaussian Whittle-Matérn Fields on Metric Graphs.” Bernoulli 30 (2): 1611–39. https://doi.org/10.3150/23-BEJ1647.

Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2024. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics 33 (1): 64–74. https://doi.org/10.1080/10618600.2023.2231051.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://CRAN.R-project.org/package=htmltools.

J, Lemon. 2006. “Plotrix: A Package in the Red Light District of r.” R-News 6 (4): 8–12.

Meschiari, Stefano. 2022. Latex2exp: Use LaTeX Expressions in Plots. https://CRAN.R-project.org/package=latex2exp.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://CRAN.R-project.org/package=here.

Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.

Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.

Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.

R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2023. scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2024. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.