Auxiliary 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 = 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)
}

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

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

1 Theoretical stuff

We consider the equation

\[ \left(\kappa^{2}-\Delta\right)^{\alpha / 2}(\tau u)=\mathcal{W} \quad \text { on } \Gamma. \tag{0} \label{spde_equation} \]

The Matérn covariance function is given by

\[ \varrho_M(h)=\frac{\tau^{-2}}{2^{\nu-1} \Gamma(\nu+n / 2)(4 \pi)^{n / 2} \kappa^{2 \nu}}(\kappa|h|)^\nu K_\nu(\kappa|h|), \tag{1} \label{matern_cov} \]

where \(n=1\) and the parameters \(\tau, \kappa>0\) and \(0<\nu \leq 1 / 2\) control the variance, practical correlation range, and the sample path regularity, respectively. Further, \(K_\nu(\cdot)\) is a modified Bessel function of the second kind and \(\Gamma(\cdot)\) denotes the gamma function.

From (Bolin, Simas, and Wallin 2023c, Proposition 5), we know that

\[ \mathbf{r}: \mathbb{R} \times \mathbb{R} \mapsto \mathbb{R}^{\alpha \times \alpha}, \quad \mathbf{r}\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}} \varrho_M\left(t_1-t_2\right)\right]_{i j \in\{1,2, \ldots, \alpha\}} \tag{2} \label{r_matrix} \] is the covariance function of \(\mathbf{x}(t)=\left[x(t), x'(t), \ldots, x^{(\alpha-1)}(t)\right]^{\top}\) if \(x\) is a centered Gaussian process on \(\mathbb{R}\) with Matérn covariance function with \(\nu = \alpha-1/2\) and \(\alpha \in \mathbb{N}\).

Let \(\mathbf{r}(\cdot, \cdot)\) be given by \(\eqref{r_matrix}\) with \(\alpha \in \mathbb{N}\). Then, from (Bolin, Simas, and Wallin 2023c, Proposition 6), for \(\ell>0\),

\[ \widetilde{\mathbf{r}}_{\ell}\left(t_1, t_2\right)=\mathbf{r}\left(t_1, t_2\right)+\left[\mathbf{r}\left(t_1, 0\right) \mathbf{r}\left(t_1, \ell\right)\right]\left[\begin{array}{cc} \mathbf{r}(0,0) & -\mathbf{r}(0, \ell) \\ -\mathbf{r}(\ell, 0) & \mathbf{r}(0,0) \end{array}\right]^{-1}\left[\begin{array}{l} \mathbf{r}\left(t_2, 0\right) \\ \mathbf{r}\left(t_2, \ell\right) \end{array}\right] \tag{3} \label{r_tilde_matrix} \]

is the multivariate covariance function of the multivariate process \(\widetilde{\mathbf{x}}(t)=\left[\widetilde{x}(t), \widetilde{x}^{\prime}(t), \ldots, \widetilde{x}^{(\alpha-1)}(t)\right]^{\top}\) on the interval \([0, \ell]\), where \(\widetilde{x}\) is the boundaryless Whittle–Matérn process on \([0, \ell]\).

Let

\[ \mathcal{K}_\alpha(x)=\left\{\omega \in \Omega: \forall v \in \mathcal{V} \text { and each pair } e, \widetilde{e} \in \mathcal{E}_v, x_e^{(2 k)}(v, \omega)=x_{\widetilde{e}}^{(2 k)}(v, \omega) \text {, and } \sum_{e \in \mathcal{E}_v} \partial_e x_e^{2 k+1}(v, \omega)=0, k=0, \ldots,\lceil\alpha-1 / 2\rceil-1\right\} \tag{4} \label{K_alpha} \]

and \(\widetilde{u} = \{\widetilde{u}_e:e\in\mathcal{E}\}\) be a family of independent boundaryless Whittle–Matérn process on the edges. Further, let

\[ \widetilde{\mathbf{u}}(s)=\left[\widetilde{u}(s), \widetilde{u}^{\prime}(s),\ldots, \widetilde{u}^{(\alpha-1)}(s)\right]^{\top}=\sum_{e \in \mathcal{E}} \mathbb{I}(s \in e) \widetilde{\mathbf{u}}_e(s) \in\mathbb{R}^{\alpha}. \tag{5} \label{u_vector} \] From (Bolin, Simas, and Wallin 2023c, Theorem 4), if we define

\[ \mathbf{u}(s) = \widetilde{\mathbf{u}}(s)| \mathcal{K}_\alpha(\widetilde{u}), \tag{6} \label{u_vector_constrained} \]

then \(u(\cdot)\), the first entry of \(\mathbf{u}(\cdot)\), is a solution to \(\eqref{spde_equation}\).

The precision matrix of \([\widetilde{\mathbf{u}}(0), \widetilde{\mathbf{u}}(\ell)]\) is

\[ \widetilde{\mathbf{Q}}_e=\mathbf{Q}_e-\frac{1}{2}\left[\begin{array}{cc} \mathbf{r}(0,0)^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{r}(0,0)^{-1} \end{array}\right]\in\mathbb{R}^{2 \alpha \times 2 \alpha}, \tag{7} \label{Q_tilde_e} \]

where \(\mathbf{Q}_e\) is the precision matrix of \([\mathbf{u}(0), \mathbf{u}(\ell)]\).

Let \[ \mathbf{U}=\left[\mathbf{u}(\underline{e}_1)^{\top}, \mathbf{u}(\overline{e}_1)^{\top}, \mathbf{u}(\underline{e}_2)^{\top}, \mathbf{u}(\overline{e}_2)^{\top}, \ldots, \mathbf{u}(\underline{e}_{|\mathcal{E}|})^{\top}, \mathbf{u}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2 \alpha|\mathcal{E}|} \tag{8} \label{U_vector} \]

and

\[ \widetilde{\mathbf{U}}=\left[\widetilde{\mathbf{u}}_1(\underline{e}_1)^{\top}, \widetilde{\mathbf{u}}_1(\overline{e}_1)^{\top}, \widetilde{\mathbf{u}}_2(\underline{e}_2)^{\top}, \widetilde{\mathbf{u}}_2(\overline{e}_2)^{\top}, \ldots, \widetilde{\mathbf{u}}_{|\mathcal{E}|}(\underline{e}_{|\mathcal{E}|})^{\top}, \widetilde{\mathbf{u}}_{|\mathcal{E}|}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2 \alpha|\mathcal{E}|} \tag{9} \label{U_vector_tilde} \] We then have that

\[ \widetilde{\mathbf{U}}\sim \mathrm{N}\left(\mathbf{0}, \widetilde{\mathbf{Q}}^{-1}\right), \quad \widetilde{\mathbf{Q}} = \operatorname{diag}\left(\{\mathbf{Q}_e\}_{e\in\mathcal{E}}\right), \tag{10} \label{U_vector_tilde_distribution} \]

where \(\mathbf{Q}_e\) is the precision matrix of \([\widetilde{\mathbf{u}}_e(\underline{e}), \widetilde{\mathbf{u}}_e(\overline{e})]\).

The Kirchhoff conditions \(\eqref{K_alpha}\) can be written as a system of linear equations

\[ \mathbf{K} \widetilde{\mathbf{U}}=\mathbf{0}, \tag{11} \label{K_alpha_matrix} \]

where \(\mathbf{K}\in\mathbb{R}^{k\times 2\alpha|\mathcal{E}|}\) is a suitable matrix and \(k\) is the number of linear constraints given by \(\eqref{K_alpha}\).

Finally,

\[ \mathbf{U}=\widetilde{\mathbf{U}}| \{\mathbf{K} \widetilde{\mathbf{U}}=\mathbf{0}\}\sim\mathrm{N}\left(\mathbf{0}, \widetilde{\mathbf{Q}}^{-1}-\widetilde{\mathbf{Q}}^{-1} \mathbf{K}^{\top}\left(\mathbf{K} \widetilde{\mathbf{Q}}^{-1} \mathbf{K}^{\top}\right)^{-1} \mathbf{K} \widetilde{\mathbf{Q}}^{-1}\right). \tag{12} \label{U_vector_constrained} \]

Let \(\mathbf{A}\) (which is not unique) be a \(|\mathcal{V}| \times 2\alpha|\mathcal{E}|\) matrix such that

\[ \mathbf{U}_v = [u(v_1),u(v_2), \dots, u(v_{|\mathcal{V}|})]^{\top} = \mathbf{A} \mathbf{U}\sim\mathrm{N}\left(\mathbf{0},\mathbf{A}\left(\widetilde{\mathbf{Q}}^{-1}-\widetilde{\mathbf{Q}}^{-1} \mathbf{K}^{\top}\left(\mathbf{K} \widetilde{\mathbf{Q}}^{-1} \mathbf{K}^{\top}\right)^{-1} \mathbf{K} \widetilde{\mathbf{Q}}^{-1}\right) \mathbf{A}^{\top}\right), \tag{13} \label{A_matrix} \]

---

1.1 Case \(\alpha=1\)

---

For \(\alpha=1, \mathbf{U}_v\) in \(\eqref{A_matrix}\) satisfies

\[ \mathbf{U}_v \sim \mathrm{~N}\left(\mathbf{0}, \mathbf{Q}^{-1}\right), \tag{14} \label{U_v_distribution_alpha1} \]

where

\[ \mathbf{Q}_{i j}=2 \kappa \tau^2 . \begin{cases}\dfrac{1}{2}\mathbb{I}(\text{deg}(v_i) = 1) + \displaystyle\sum_{e \in \mathcal{E}_{v_i}}\left\{\left(\frac{1}{2}+\frac{e^{-2 \kappa \ell_e}}{1-e^{-2 \kappa \ell_e}}\right) \mathbb{I}(\bar{e} \neq \underline{e})+\tanh \left(\kappa \frac{l_e}{2}\right) \mathbb{I}(\bar{e}=\underline{e})\right\} & \text { if } i=j \\ \displaystyle\sum_{e \in \mathcal{E}_{v_i} \cap \mathcal{E}_{v_j}}-\frac{e^{-\kappa \ell_e}}{1-e^{-2 \kappa \ell_e}} & \text { if } i \neq j\end{cases} \tag{15} \label{Q_matrix} \]

and we have added \(\dfrac{1}{2}\mathbb{I}(\text{deg}(v_i) = 1)\) to correct for stationarity at vertices with degree one (see (Bolin, Simas, and Wallin 2023c, sec. 7)).

---

1.2 Case \(\alpha>1\)

---

Let \(\alpha>0\) and \(\mathbf{T}\) be the change-of-basis matrix such that

\[ \widetilde{\mathbf{U}}^\star = \mathbf{T} \widetilde{\mathbf{U}}, \tag{16} \label{T_matrix} \]

and the \(k\) constraints of \(\mathbf{K}\) are given by the first \(k\) rows of \(\widetilde{\mathbf{U}}^\star\). That is,

\[ \widetilde{\mathbf{U}}^\star\sim \mathrm{N}\left(\mathbf{0}, \mathbf{T}\widetilde{\mathbf{Q}}^{-1}\mathbf{T}^\top\right) \tag{17} \label{U_star_distribution} \]

Let \(\widetilde{\mathbf{Q}}^*=\mathbf{T} \widetilde{\mathbf{Q}} \mathbf{T}^{\top}\). Further, let \(\mathbf{T}_{\mathcal{U}}\) denote the matrix obtained by removing the first \(k\) rows from \(\mathbf{T}\), and let \(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\) denote the matrix obtained by removing the first \(k\) rows and the first \(k\) columns of \(\widetilde{\mathbf{Q}}^*\). Then

\[ \mathbf{U} \sim \mathrm{N}\left(\mathbf{0}, \mathbf{T}_{\mathcal{U}}^{\top}\left(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\right)^{-1} \mathbf{T}_{\mathcal{U}}\right) \mathbb{I}(\mathbf{K} \mathbf{U}=\mathbf{0}) \tag{18} \label{U_distribution_final} \] and

\[ \mathbf{U}_v \sim \mathrm{~N}\left(\mathbf{0}, \mathbf{A} \mathbf{T}_{\mathcal{U}}^{\top}\left(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\right)^{-1} \mathbf{T}_{\mathcal{U}} \mathbf{A}^{\top}\right). \tag{19} \label{U_v_distribution_final} \]

The matrices \(\mathbf{A}, \mathbf{T}_{\mathcal{U}}\) and \(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\) in \(\eqref{U_v_distribution_final}\) are sparse. Thus, we can simulate \(\mathbf{U}_v\) efficiently by simulating

\[ \mathbf{v} \sim \mathrm{N}\left(\mathbf{0},\left(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\right)^{-1}\right) \tag{20} \label{v_distribution} \] through sparse Cholesky factorization of \(\widetilde{\mathbf{Q}}_{\mathcal{U} \mathcal{U}}^*\) and then computing the sparse matrix vector product \(\mathbf{U}_v=\mathbf{A T}_{\mathcal{U}}^{\top} \mathbf{v}\).

---

2 Likelihood evaluation

---

Assume we have observations \(\mathbf{y} = [y_1,\dots,y_n]\) at locations \(s_1,\dots,s_n\in\Gamma\) such that (1) \(y_i = u(s_i)\) or (2) \(y_i|u(\cdot)\sim \mathrm{N}(u(s_i), \sigma^2)\).

---

2.1 Case \(\alpha=1\)

---

We add the locations \(s_1,\dots,s_n\) as vertices and called the extended graph as \(\overline{\Gamma}\). Let \(\mathbf{v} = (v_1, \dots, v_m)\) be the indices of the original vertices and \(\mathbf{s} = (s_1, \dots, s_n)\) the indices of the added locations.

3 Fractional case

We consider the covariance operator \(\mathcal{C} = \tau^{-2}L^{-\alpha}\), where \(L = \kappa^2 - \Delta_\Gamma\) and \(\alpha > 0\) is not necessarily an integer.

We consider an approximation

\[ \mathcal{C}_{m,\alpha} = \tau^{-2}L^{-\lfloor \alpha \rfloor} p(L^{-1})q(L^{-1})^{-1} = \tau^{-2}L^{-\lfloor \alpha \rfloor} \left(\sum_{i=1}^m r_i (L-p_iI)^{-1} + kI\right) = \tau^{-2}\left(kL^{-\lfloor \alpha \rfloor} + \sum_{i=1}^m r_i L^{-\lfloor \alpha \rfloor}(L-p_iI)^{-1}\right), \tag{21} \label{C_fractional_approx} \] So \(u(s) = u_0(s) + \sum_{i=1}^m u_i(s)\), where \(u_0\) and \(u_i\) satisfy

\[ \dfrac{1}{\sqrt{k}}(\kappa^2 - \Delta_\Gamma)^{\lfloor \alpha \rfloor / 2} (\tau u_0) = \mathcal{W},\text{ on }\Gamma \tag{22} \label{u0_equation} \]

and

\[ \dfrac{1}{\sqrt{r_i}}((\kappa^2 - \Delta_\Gamma)^{\lfloor \alpha \rfloor}(\kappa^2 - p_i-\Delta_\Gamma))^{1/2} (\tau u_i) = \mathcal{W}, \text{ on }\Gamma, \quad i = 1, \dots, m. \tag{23} \label{ui_equation} \]

3.1 Case \(\alpha \in (0,1)\)

In this case, \(\eqref{u0_equation}\) reduces to

\[ \dfrac{1}{\sqrt{k}}\tau u_0 = \mathcal{W},\text{ on }\Gamma \tag{24} \label{u0_equation_alpha01} \]

and \(\eqref{ui_equation}\) reduces to

\[ \dfrac{1}{\sqrt{r_i}}(\kappa^2 - p_i - \Delta_\Gamma)^{1/2} (\tau u_i) = \mathcal{W}, \text{ on }\Gamma, \quad i = 1, \dots, m. \tag{25} \label{ui_equation_alpha01} \]

Equation \(\eqref{ui_equation_alpha01}\) is just the case \(\alpha=1\) in the original setting but now with shifted parameter \(\kappa^2 - p_i\) instead of \(\kappa^2\) and with a constant factor \(1/\sqrt{r_i}\).

3.2 Case \(\alpha \in (1,2)\)

In this case, \(\eqref{u0_equation}\) reduces to \[ \dfrac{1}{\sqrt{k}}(\kappa^2 - \Delta_\Gamma)^{1/2}(\tau u_0) = \mathcal{W},\text{ on }\Gamma \tag{26} \label{u0_equation_alpha12} \]

and \(\eqref{ui_equation}\) reduces to

\[ \dfrac{1}{\sqrt{r_i}}((\kappa^2 - \Delta_\Gamma)(\kappa^2 - p_i - \Delta_\Gamma))^{1/2} (\tau u_i) = \mathcal{W}, \text{ on }\Gamma,\quad i = 1, \dots, m. \tag{27} \label{ui_equation_alpha12} \] Equation \(\eqref{u0_equation_alpha12}\) is just the case \(\alpha=1\) in the original setting but now with a constant factor \(1/\sqrt{k}\).

3.3 Case \(\alpha \in (2,3)\)

In this case, \(\eqref{u0_equation}\) reduces to

\[ \dfrac{1}{\sqrt{k}}(\kappa^2 - \Delta_\Gamma)(\tau u_0) = \mathcal{W}, \tag{28} \label{u0_equation_alpha23} \] and \(\eqref{ui_equation}\) reduces to \[ \dfrac{1}{\sqrt{r_i}}((\kappa^2 - \Delta_\Gamma)^2(\kappa^2 - p_i - \Delta_\Gamma))^{1/2} (\tau u_i) = \mathcal{W}, \quad i = 1, \dots, m. \tag{29} \label{ui_equation_alpha23} \]

Equation \(\eqref{u0_equation_alpha23}\) is just the case \(\alpha=2\) in the original setting but now with a constant factor \(1/\sqrt{k}\).

3.4 Matrix version

Let

\[ \mathbf{u}(s) = [u(s), u'(s), \dots, u^{(\lfloor \alpha \rfloor)}(s)]^\top\in\mathbb{R}^{\lfloor \alpha \rfloor + 1} \tag{30} \label{u_vector_fractional} \]

and

\[ \widetilde{\mathbf{u}}(s)=\left[\widetilde{u}(s), \widetilde{u}^{\prime}(s),\ldots, \widetilde{u}^{(\lfloor \alpha \rfloor)}(s)\right]^{\top} \in\mathbb{R}^{\lfloor \alpha \rfloor + 1} \tag{31} \]

and

\[ \mathbf{u}(\mathbf{s}) = \widetilde{\mathbf{u}}(\mathbf{s})| \mathcal{K}_{\lfloor\alpha\rfloor}(\widetilde{u}), \tag{32} \]

and

\[ \mathbf{u}(\mathbf{s}) = [\mathbf{u}(s_1), \mathbf{u}(s_2), \dots, \mathbf{u}(s_n)]^\top \in\mathbb{R}^{n(\lfloor \alpha \rfloor + 1)} \]

and

\[ \widetilde{\mathbf{u}}(\mathbf{s}) = [\widetilde{\mathbf{u}}(s_1), \widetilde{\mathbf{u}}(s_2), \dots, \widetilde{\mathbf{u}}(s_n)]^\top \in\mathbb{R}^{n(\lfloor \alpha \rfloor + 1)} \] and

\[ \mathbf{s} = [s_1, s_2, \dots, s_n]^\top\in\mathbb{R}^n,\quad s_1,s_2,\dots,s_n\in\Gamma \]

and

\[ \mathbf{U}=\left[\mathbf{u}(\underline{e}_1)^{\top}, \mathbf{u}(\overline{e}_1)^{\top}, \mathbf{u}(\underline{e}_2)^{\top}, \mathbf{u}(\overline{e}_2)^{\top}, \ldots, \mathbf{u}(\underline{e}_{|\mathcal{E}|})^{\top}, \mathbf{u}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2(\lfloor \alpha \rfloor + 1)|\mathcal{E}|} \tag{33} \]

and

\[ \widetilde{\mathbf{U}}=\left[\widetilde{\mathbf{u}}_1(\underline{e}_1)^{\top}, \widetilde{\mathbf{u}}_1(\overline{e}_1)^{\top}, \widetilde{\mathbf{u}}_2(\underline{e}_2)^{\top}, \widetilde{\mathbf{u}}_2(\overline{e}_2)^{\top}, \ldots, \widetilde{\mathbf{u}}_{|\mathcal{E}|}(\underline{e}_{|\mathcal{E}|})^{\top}, \widetilde{\mathbf{u}}_{|\mathcal{E}|}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2(\lfloor \alpha \rfloor + 1)|\mathcal{E}|} \tag{34} \]

We then have that

\[ \mathbf{u}_0(s) = \widetilde{\mathbf{u}}_0(s)| \mathcal{K}_{\lfloor\alpha\rfloor-1}(\widetilde{u}_0), \tag{35} \]

and

\[ \mathbf{u}_i(s) = \widetilde{\mathbf{u}}_i(s)| \mathcal{K}_{\lfloor\alpha\rfloor}(\widetilde{u}_i),\quad i=1,\dots,m, \tag{35} \]

where

\[ \widetilde{\mathbf{u}}_0(s)=\left[\widetilde{u}_0(s), \widetilde{u}_0^{\prime}(s),\ldots, \widetilde{u}_0^{(\lfloor \alpha \rfloor-1)}(s)\right]^{\top} \in\mathbb{R}^{\lfloor \alpha \rfloor} \tag{36} \]

has covariance function

Covariance of the boundaryless process on an interval \([0, \ell]\)

\[ \widetilde{\mathbf{r}}\left(t_1, t_2\right)=\mathbf{r}\left(t_1, t_2\right)+\left[\mathbf{r}\left(t_1, 0\right) \mathbf{r}\left(t_1, \ell\right)\right]\left[\begin{array}{cc} \mathbf{r}(0,0) & -\mathbf{r}(0, \ell) \\ -\mathbf{r}(\ell, 0) & \mathbf{r}(0,0) \end{array}\right]^{-1}\left[\begin{array}{l} \mathbf{r}\left(t_2, 0\right) \\ \mathbf{r}\left(t_2, \ell\right) \end{array}\right] \tag{37} \]

Covariance of the corresponding unrestricted process on \(\mathbb{R}\)

\[ \mathbf{r}: \mathbb{R} \times \mathbb{R} \mapsto \mathbb{R}^{\lfloor\alpha \rfloor\times \lfloor\alpha\rfloor}, \quad \mathbf{r}\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}} \varrho_M\left(t_1-t_2\right)\right]_{i j \in\{1,2, \ldots, \lfloor\alpha\rfloor\}} \tag{38} \]

and

\[ \widetilde{\mathbf{u}}_i(s)=\left[\widetilde{u}_i(s), \widetilde{u}_i^{\prime}(s),\ldots, \widetilde{u}_i^{(\lfloor \alpha \rfloor)}(s)\right]^{\top} \in\mathbb{R}^{\lfloor \alpha \rfloor+1} \tag{39} \] has covariance function

Covariance of the shifted boundaryless process on an interval \([0, \ell]\)

\[ \widetilde{\mathbf{r}}_i\left(t_1, t_2\right)=\mathbf{r}_i\left(t_1, t_2\right)+\left[\mathbf{r}_i\left(t_1, 0\right) \mathbf{r}_i\left(t_1, \ell\right)\right]\left[\begin{array}{cc} \mathbf{r}_i(0,0) & -\mathbf{r}_i(0, \ell) \\ -\mathbf{r}_i(\ell, 0) & \mathbf{r}_i(0,0) \end{array}\right]^{-1}\left[\begin{array}{l} \mathbf{r}_i\left(t_2, 0\right) \\ \mathbf{r}_i\left(t_2, \ell\right) \end{array}\right] \tag{40} \]

Covariance of the corresponding unrestricted process on \(\mathbb{R}\)

\[ \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{41} \]

The precision matrix of \([\widetilde{\mathbf{u}}_{0,e}(0), \widetilde{\mathbf{u}}_{0,e}(\ell)]\) is

\[ \widetilde{\mathbf{Q}}_{0,e}=\mathbf{Q}_{0,e}-\frac{1}{2}\left[\begin{array}{cc} \mathbf{r}(0,0)^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{r}(0,0)^{-1} \end{array}\right]\in\mathbb{R}^{2 \lfloor\alpha\rfloor \times 2 \lfloor\alpha\rfloor}, \]

where \(\mathbf{Q}_{0,e}\) is the precision matrix of \([\mathbf{u}_{0,e}(0), \mathbf{u}_{0,e}(\ell)]\).

The precision matrix of \([\widetilde{\mathbf{u}}_{i,e}(0), \widetilde{\mathbf{u}}_{i,e}(\ell)]\) is

\[ \widetilde{\mathbf{Q}}_{i,e}=\mathbf{Q}_{i,e}-\frac{1}{2}\left[\begin{array}{cc} \mathbf{r}_i(0,0)^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{r}_i(0,0)^{-1} \end{array}\right]\in\mathbb{R}^{2 \lceil\alpha\rceil \times 2 \lceil\alpha\rceil}, \]

where \(\mathbf{Q}_{i,e}\) is the precision matrix of \([\mathbf{u}_{i,e}(0), \mathbf{u}_{i,e}(\ell)]\).

Define

\[ \widetilde{\mathbf{U}}_0=\left[\widetilde{\mathbf{u}}_{0,1}(\underline{e}_1)^{\top}, \widetilde{\mathbf{u}}_{0,1}(\overline{e}_1)^{\top}, \widetilde{\mathbf{u}}_{i,2}(\underline{e}_2)^{\top}, \widetilde{\mathbf{u}}_{0,2}(\overline{e}_2)^{\top}, \ldots, \widetilde{\mathbf{u}}_{0,|\mathcal{E}|}(\underline{e}_{|\mathcal{E}|})^{\top}, \widetilde{\mathbf{u}}_{0,|\mathcal{E}|}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2\lfloor \alpha \rfloor|\mathcal{E}|} \tag{42} \]

and

\[ \widetilde{\mathbf{U}}_i=\left[\widetilde{\mathbf{u}}_{i,1}(\underline{e}_1)^{\top}, \widetilde{\mathbf{u}}_{i,1}(\overline{e}_1)^{\top}, \widetilde{\mathbf{u}}_{i,2}(\underline{e}_2)^{\top}, \widetilde{\mathbf{u}}_{i,2}(\overline{e}_2)^{\top}, \ldots, \widetilde{\mathbf{u}}_{i,|\mathcal{E}|}(\underline{e}_{|\mathcal{E}|})^{\top}, \widetilde{\mathbf{u}}_{i,|\mathcal{E}|}(\overline{e}_{|\mathcal{E}|})^{\top}\right]^{\top}\in\mathbb{R}^{2(\lfloor \alpha \rfloor + 1)|\mathcal{E}|},\quad i=0,\dots,m, \tag{43} \]

Then

\[ \widetilde{\mathbf{U}}_0\sim \mathrm{N}\left(\mathbf{0}, \widetilde{\mathbf{Q}}_0^{-1}\right), \quad \widetilde{\mathbf{Q}}_0 = \operatorname{diag}\left(\{\widetilde{\mathbf{Q}}_{0,e}\}_{e\in\mathcal{E}}\right), \]

and

\[ \widetilde{\mathbf{U}}_i\sim \mathrm{N}\left(\mathbf{0}, \widetilde{\mathbf{Q}}_i^{-1}\right), \quad\widetilde{\mathbf{Q}}_i = \operatorname{diag}\left(\{\widetilde{\mathbf{Q}}_{i,e}\}_{e\in\mathcal{E}}\right),\quad i=1,\dots,m. \]

4 Auxiliary functions

4.1 Function gets.graph.tadpole()

Given a mesh size h, function gets.graph.tadpole() builds a tadpole graph and creates a mesh.

# Function to build a tadpole graph and create a mesh
gets.graph.tadpole <- function(flip_edge = FALSE){
  if(flip_edge) {
    edge1 <- rbind(c(0,0),c(1,0))[c(2,1),]
    } else {
    edge1 <- rbind(c(0,0),c(1,0))}
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges <- list(edge1, edge2)
  graph <- metric_graph$new(edges = edges, verbose = 0)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  #graph$build_mesh(h = h)
  return(graph)
}

# Eigenfunctions for the tadpole graph
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}


# Function to compute the true covariance matrix
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

Qalpha1 <- function(theta, graph, BC = 1, build = TRUE) {
  
  kappa <- theta[2]
  tau <- theta[1]
  i_ <- j_ <- x_ <- rep(0, dim(graph$V)[1]*4) # it has length 4*nV
  count <- 0
  for(i in 1:graph$nE){ # loop over edges
    l_e <- graph$edge_lengths[i]
    c1 <- exp(-kappa*l_e)
    c2 <- c1^2
    one_m_c2 = 1-c2
    c_1 = 0.5 + c2/one_m_c2
    c_2 = -c1/one_m_c2
    
    if (graph$E[i, 1] != graph$E[i, 2]) { # This is for non-circular edges and codes I(overline(e) != underline(e))
      
      i_[count + 1] <- graph$E[i, 1]
      j_[count + 1] <- graph$E[i, 1]
      x_[count + 1] <- c_1
      
      i_[count + 2] <- graph$E[i, 2]
      j_[count + 2] <- graph$E[i, 2]
      x_[count + 2] <- c_1
      
      
      i_[count + 3] <- graph$E[i, 1]
      j_[count + 3] <- graph$E[i, 2]
      x_[count + 3] <- c_2
      
      i_[count + 4] <- graph$E[i, 2]
      j_[count + 4] <- graph$E[i, 1]
      x_[count + 4] <- c_2
      count <- count + 4
    }else{ # This is for circular edges and codes I(overline(e) = underline(e))
      i_[count + 1] <- graph$E[i, 1]
      j_[count + 1] <- graph$E[i, 1]
      x_[count + 1] <- tanh(0.5 * kappa * l_e)
      count <- count + 1
    }

  }
  if(BC == 1){
    #does this work for circle?
    i.table <- table(i_[1:count])
    index = as.integer(names(which(i.table < 3)))
    i_ <- c(i_[1:count], index)
    j_ <- c(j_[1:count], index)
    x_ <- c(x_[1:count], rep(0.5, length(index))) # here is where we add the 0.5 for degree one vertices
    count <- count + length(index)
    # print(i_)
    # print(j_)
  }else if(BC==2){
    
    dV <- graph$get_vertices()$degree
    index <- 1:length(dV)
    i_ <- c(i_[1:count], index)
    j_ <- c(j_[1:count], index)
    x_ <- c(x_[1:count], -0.5*dV + 1)
    count <- count + length(index)
    
  }
  if(build){
    Q <- Matrix::sparseMatrix(i = i_[1:count],
                              j = j_[1:count],
                              x = (2 * kappa * tau^2) * x_[1:count], # This is the 2kappa*tau^2 factor
                              dims = c(graph$nV, graph$nV))
    
    
    return(Q)
  } else {
    return(list(i = i_[1:count],
                j = j_[1:count],
                x = (2 * kappa * tau^2) * x_[1:count],
                dims = c(graph$nV, graph$nV)))
  }
}

gives.indices <- function(graph, factor, constant){
  index.obs1 <- sapply(graph$PtV, 
                       function(i){
                         idx_temp <- i == graph$E[,1]
                         idx_temp <- which(idx_temp)
                         return(idx_temp[1])}
                       )
  index.obs1 <- (index.obs1 - 1) * factor + 1
  index.obs4 <- NULL
  na_obs1 <- is.na(index.obs1)
  if(any(na_obs1)){
    idx_na <- which(na_obs1)
    PtV_NA <- graph$PtV[idx_na]
    index.obs4 <- sapply(PtV_NA, 
                         function(i){
                           idx_temp <- i == graph$E[,2]
                           idx_temp <- which(idx_temp)
                           return(idx_temp[1])}
                         )
    index.obs1[na_obs1] <- (index.obs4 - 1 ) * factor + constant                                                                      
  }
  return(index.obs1)
}

conditioning <- function(graph, alpha = 1){
  i_  =  rep(0, 2 * graph$nE)
  j_  =  rep(0, 2 * graph$nE)
  x_  =  rep(0, 2 * graph$nE)

  count_constraint <- 0
  count <- 0
  for (v in 1:graph$nV) {
    edges_leaving_v  <- which(graph$E[, 1] %in% v) 
    edges_entering_v  <- which(graph$E[, 2] %in% v)
    n_leaving_edges <- length(edges_leaving_v)
    n_entering_edges <- length(edges_entering_v)
    n_e <- n_leaving_edges + n_entering_edges
    if (n_e > 1) { # the alternative is n_e = 1, which means v is a degree one vertex and so no conditioning is needed 
      if (n_entering_edges == 0) {
        edges <- cbind(edges_leaving_v, 1)
      } else if(n_leaving_edges == 0){
        edges <- cbind(edges_entering_v, 2)
      }else{
        edges <- rbind(cbind(edges_leaving_v, 1),
                       cbind(edges_entering_v, 2))
      }
      for (i in 2:n_e) {
        i_[count + 1:2] <- count_constraint + 1
        j_[count + 1:2] <- c(2 * (edges[i-1,1] - 1) + edges[i-1, 2],
                             2 * (edges[i,1]   - 1) + edges[i,   2])
        x_[count + 1:2] <- c(1,-1)
        count <- count + 2
        count_constraint <- count_constraint + 1
      }
    }
  }
  K <- Matrix::sparseMatrix(i = i_[1:count],
                            j = j_[1:count],
                            x = x_[1:count],
                            dims = c(count_constraint, 2*graph$nE))
                         
  CB <- MetricGraph:::c_basis2(K)
  return(CB)
}

Function matern.p.joint() 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^{k-1}}{\partial t_2^{k-1}} \frac{\partial^{j-1}}{\partial t_1^{j-1}} \varrho_{i}\left(t_1-t_2\right)\right]_{k j \in\{1,2, \ldots, \lceil\alpha\rceil\}} \]

where \(\beta = \lfloor\alpha\rfloor\) if \(i=0\) and \(\beta = \lceil\alpha\rceil\) if \(i=1,\dots,m\), and

\[ \varrho_{i}(h)=\begin{cases}\dfrac{\varrho_{m, i}^\alpha(h)}{k \sigma^2}, & i=0 \\ \dfrac{\varrho_{m, i}^\alpha(h)}{r_i \sigma^2}, & i=1,\dots,m\end{cases} \]

and

\[ \varrho_{m, i}^\alpha(h)= \begin{cases} k \sigma^2 \varrho\left(h ;\lfloor\alpha\rfloor-\frac{1}{2}, \kappa, \frac{c_\alpha}{c_{\lfloor\alpha\rfloor}}\right), & i=0 \\ r_i \sigma^2 \left[\dfrac{1}{p_i^{\lfloor\alpha\rfloor}} \varrho\left(h ; \frac{1}{2}, \kappa_i, \frac{c_\alpha \sqrt{\pi}}{\sqrt{1-p_i}}\right)-\displaystyle\sum_{j=1}^{\lfloor\alpha\rfloor} \dfrac{1}{p_i^{\lfloor\alpha\rfloor+1-j}} \varrho\left(h ; j-\frac{1}{2}, \kappa, \frac{c_\alpha}{c_j}\right)\right], & i=1,\dots,m\end{cases} \]

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

\[ \varrho(h; \nu, \kappa, \sigma^2)=\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} \]

---

5 Case \(i=0\)

---

\[ \widetilde{\mathbf{Q}}_{0,e}=\mathbf{Q}_{0,e}-\frac{1}{2}\left[\begin{array}{cc} \mathbf{r}(0,0)^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{r}(0,0)^{-1} \end{array}\right]\in\mathbb{R}^{2 \lfloor\alpha\rfloor \times 2 \lfloor\alpha\rfloor}, \]

This is computed in R as follows

• \(\mathbf{r}(0,0) =\) matern.p.joint(s = 0, t = 0, kappa = kappa, p = 0, alpha = floor(alpha)).

• \(\mathbf{Q}_{0,e} =\) matern.p.precision(loc = c(0, L_e[e]), kappa = kappa, p = 0, equally_spaced = TRUE, alpha = floor(alpha))$Q.

• Constant correction: \(\widetilde{\mathbf{Q}}_{0,e}=\dfrac{2\tau^2\kappa^{2\alpha}}{k\kappa}\widetilde{\mathbf{Q}}_{0,e}\)

---

6 Case \(i=1,\dots,m\)

---

\[ \widetilde{\mathbf{Q}}_{i,e}=\mathbf{Q}_{i,e}-\frac{1}{2}\left[\begin{array}{cc} \mathbf{r}_i(0,0)^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{r}_i(0,0)^{-1} \end{array}\right]\in\mathbb{R}^{2 \lceil\alpha\rceil \times 2 \lceil\alpha\rceil}, \]

This is computed in R as follows

• \(\mathbf{r}_i(0,0) =\) matern.p.joint(s = 0, t = 0, kappa = kappa, p = p_i, alpha = alpha).

• \(\mathbf{Q}_{i,e} =\) matern.p.precision(loc = c(0, L_e[e]), kappa = kappa, p = p_i, equally_spaced = TRUE, alpha = alpha)$Q.

• Constant correction: \(\widetilde{\mathbf{Q}}_{i,e}=\dfrac{2\tau^2\kappa^{2\alpha}c_\alpha\sqrt{\pi}}{r_i \kappa}\widetilde{\mathbf{Q}}_{i,e}\)

Then

\[ \widetilde{\mathbf{Q}}_0 = \operatorname{diag}\left(\{\widetilde{\mathbf{Q}}_{0,e}\}_{e\in\mathcal{E}}\right), \]

and

\[ \widetilde{\mathbf{Q}}_i = \operatorname{diag}\left(\{\widetilde{\mathbf{Q}}_{i,e}\}_{e\in\mathcal{E}}\right),\quad i=1,\dots,m. \]

# This is the correct version, it is corrected the constants
gets_cov_mat_rat_approx_alpha_1_to_2 <- function(graph, kappa, tau, alpha, m){
  
  if(alpha == 2){
    Q_unconstraint <- MetricGraph:::Qalpha2(theta = c(tau, kappa),
                                         graph = graph,
                                         BC = 3000,
                                         build = TRUE)
    graph$buildC(alpha = alpha, edge_constraint = TRUE) # should always be TRUE
    COND <- graph$CoB
    Tc <- COND$T[-c(1:length(COND$S)), ]

    Q_U <-  Tc %*% Q_unconstraint %*% t(Tc)

    index.obs_i <- gives.indices(graph = graph, factor = 4, constant = 3)
    A <- t(Tc)[index.obs_i, ] 

    Sigma <- A %*% solve(Q_U, t(A))
    return(Sigma)
  }

  # get rational approximation coefficients
  coeff <- rSPDE:::interp_rational_coefficients(
    order = m, 
    type_rational_approx = "chebfun", 
    type_interp = "spline", 
    alpha = alpha)
  
  r <- coeff$r
  p <- coeff$p
  k <- coeff$k
  
  # compute parameters
  fa <- floor(alpha)
  ca <- ceiling(alpha)
  
  nu <- alpha - 1/2
  sigma <- sqrt(gamma(nu) / (tau^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))
  c_alpha <- gamma(alpha)/gamma(alpha - 0.5)
  c_1 <- gamma(fa)/gamma(fa - 0.5)
  
  # get edge lengths
  L_e <- graph$edge_lengths
  
  # initialize Qtilde_i, a list containing block diagonal matrices with blocks Qtilde_{i,e} for each i
  Qtilde_i <- list() 
  for(i in 1:m){
    
    # compute r_(0,0)
    r00 <- matern.p.joint(
      s = 0, 
      t = 0, 
      kappa = kappa, 
      p = p[i], 
      alpha = alpha)
    
    # compute r_(0,0)^(-1)
    r00_inverse <- solve(r00, Diagonal(ca))
    
    # define zero block 
    zero_block <- matrix(0, ca, ca)
    
    # build correction term
    correction_term <- rbind(
      cbind(r00_inverse, zero_block),
      cbind(zero_block, r00_inverse))
    
    # initialize Qtilde_i[[i]], a list containing Qtilde_{i,e} for each edge e
    Qtilde_i[[i]] <- list()
    for(e in 1:length(L_e)){
      
      # compute Q_{i,e}
      Q_e <- matern.p.precision(
        loc = c(0, L_e[e]),
        kappa = kappa, 
        p = p[i],
        equally_spaced = FALSE, 
        alpha = alpha)$Q
      
      # store Qtilde_{i,e}
      Qtilde_i[[i]][[e]] <- Q_e - 0.5 * correction_term
    }
    # build block diagonal matrix Qtilde_i[[i]]
    Qtilde_i[[i]] <- bdiag(Qtilde_i[[i]])
  }
  
  # --------------------------------------------------
  # CASE i = 0
  # --------------------------------------------------
  
  # When I use MetricGraph:::Qalpha1, I am assuming that tau and sigma are related by tau^2 = gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)) where nu = 1/2
  
  inv_factor_0 <- 1/(k*c_alpha/c_1)
  NU <- fa - 0.5
  TAU <- sqrt(gamma(NU) / (sigma^2 * kappa^(2*NU) * (4*pi)^(1/2) * gamma(NU + 1/2)))
    
  Qtilde_0_star_UU <- MetricGraph:::Qalpha1(
    theta = c(TAU, kappa), 
    graph = graph, 
    BC = 3000, 
    build = TRUE) * inv_factor_0
  
  A_0 <- graph$.__enclos_env__$private$A()

  # --------------------------------------------------
  # CASE i = 1,...,m
  # --------------------------------------------------
  
  # build conditioning matrix
  graph$buildC(alpha = 2, edge_constraint = TRUE) # should always be TRUE
  COND_i <- graph$CoB
  Tc <- COND_i$T[-c(1:length(COND_i$S)), ]
  
  inv_factor_i <- 1/(r*sigma^2)

  Qtilde_i_star_UU <- purrr::map2(
    Qtilde_i, 
    inv_factor_i, 
    function(Q, x) Tc %*% Q %*% t(Tc) * x)

  index.obs_i <- gives.indices(graph = graph, factor = 4, constant = 3)
  A_i <- t(Tc)[index.obs_i, ] 
  
  # Build matrix A and Q_UU
  A <- cbind(A_0, do.call(cbind, rep(list(A_i), m)))
  Q_UU <- bdiag(Qtilde_0_star_UU, bdiag(Qtilde_i_star_UU))
  # Return Sigma
  Sigma <- A %*% solve(Q_UU, t(A)) 
  return(Sigma)
}

gets_cov_mat_rat_approx_alpha_0_to_1 <- function(graph, kappa, tau, alpha, m){
  
  if(alpha == 1){
    I <- Matrix::Diagonal(graph$nV)  
    Q_U <- MetricGraph:::Qalpha1(
      theta = c(tau, kappa),
      graph = graph,
      BC = 3000,
      build = TRUE)
    Sigma <- solve(Q_U, I)
    return(Sigma)
  }

  coeff <- rSPDE:::interp_rational_coefficients(
    order = m, 
    type_rational_approx = "chebfun", 
    type_interp = "spline", 
    alpha = alpha)
  
  r <- coeff$r
  p <- coeff$p
  k <- coeff$k
  
  nu <- alpha - 1/2
  sigma <- sqrt(gamma(nu) / (tau^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))
  
  fa <- floor(alpha)
  ca <- ceiling(alpha)
  
  c_alpha <- gamma(alpha)/gamma(alpha - 0.5)
  
  # --------------------------------------------------
  # CASE i = 0
  # --------------------------------------------------
  
  inv_factor_0 <- 1/(k*c_alpha*sqrt(4*pi)*sigma^2/kappa)
  I <- Matrix::Diagonal(graph$nV)     
  Qtilde_0_star_UU <- I * inv_factor_0
  
  # --------------------------------------------------
  # CASE i = 1,...,m
  # --------------------------------------------------
  
  inv_factor_i <- 1/(r*c_alpha*sqrt(pi)/sqrt(1 - p))
    
  NU <- ca - 0.5
  
  Qtilde_i_star_UU <- list()
  for(i in 1:m){
    
    KAPPA <- kappa*sqrt(1 - p[i])
    TAU <- sqrt(gamma(NU) / (sigma^2 * KAPPA^(2*NU) * (4*pi)^(1/2) * gamma(NU + 1/2)))
    
    Qtilde_i_star_UU[[i]] <- MetricGraph:::Qalpha1(
      theta = c(TAU, KAPPA), 
      graph = graph, 
      BC = 3000, 
      build = TRUE) * inv_factor_i[i]
  
  }

  Q_UU <- c(list(Qtilde_0_star_UU), Qtilde_i_star_UU)
  
  Sigma <- Reduce(`+`, lapply(Q_UU, function(Qi) solve(Qi, I)))
  return(Sigma)
}

rat_covariance <- function(graph, kappa, tau, alpha, m){
  if(alpha <= 0.5){
    stop("alpha = ", alpha, ", alpha should be larger than 0.5")
  }
  else if(alpha > 0.5 && alpha <= 1){
    return(gets_cov_mat_rat_approx_alpha_0_to_1(graph = graph,
                                                kappa = kappa,
                                                tau = tau,
                                                alpha = alpha,
                                                m = m))
  }
  else if(alpha > 1 && alpha <= 2){
    return(gets_cov_mat_rat_approx_alpha_1_to_2(graph = graph,
                                                kappa = kappa,
                                                tau = tau,
                                                alpha = alpha,
                                                m = m))
  }
}

{r}
# Using V's implementation
gets_cov_mat_rat_approx_alpha_1_to_2 <- function(graph, kappa, tau, alpha, m){

  # get rational approximation coefficients
  coeff <- rSPDE:::interp_rational_coefficients(
    order = m, 
    type_rational_approx = "chebfun", 
    type_interp = "spline", 
    alpha = alpha)
  
  r <- coeff$r
  p <- coeff$p
  k <- coeff$k
  
  # reparameterization
  nu <- alpha - 1/2
  sigma <- sqrt(gamma(nu) / (tau^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))
  c_alpha <- gamma(alpha)/gamma(alpha - 0.5)
  c_1 <- gamma(floor(alpha))/gamma(floor(alpha) - 0.5)
  
  # get edge lengths
  L_e <- graph$edge_lengths
  
  Qtilde_i <- list() 
  for(order in 0:m){
    if(order == 0){
      P <- order
      ALPHA <- floor(alpha)
      FACTOR <- (2*tau^2)/(k*kappa)
      r00_inverse <- solve(matern.p.joint(s = 0, t = 0, kappa = kappa, p = P, alpha = ALPHA))
      correction_term <- rbind(cbind(r00_inverse, matrix(0, floor(alpha), floor(alpha))),
                               cbind(matrix(0, floor(alpha), floor(alpha)), r00_inverse))
    } else {
      P <- p[order]
      ALPHA <- alpha
      FACTOR <- (2*c_alpha*sqrt(pi)*tau^2)/(r[order] * kappa)
      r00_inverse <- solve(matern.p.joint(s = 0, t = 0, kappa = kappa, p = P, alpha = ALPHA))
      correction_term <- rbind(cbind(r00_inverse, matrix(0, ceiling(alpha), ceiling(alpha))),
                               cbind(matrix(0, ceiling(alpha), ceiling(alpha)), r00_inverse))
    }
    Qtilde_i[[paste0("m=",order)]] <- list()
    for(e in 1:length(L_e)){
      Q_e <- matern.p.precision(loc = c(0, L_e[e]), 
                                kappa = kappa, 
                                p = P,
                                equally_spaced = TRUE, 
                                alpha = ALPHA)$Q
      Qtilde_i[[paste0("m=",order)]][[e]] <- (Q_e - 0.5 * correction_term)*FACTOR*kappa^(2*alpha)
    }
    Qtilde_i[[paste0("m=",order)]] <- bdiag(Qtilde_i[[paste0("m=",order)]])
  }
  
  Qtilde_0 <- Qtilde_i[[paste0("m=",0)]] # extract Qtilde_0
  Qtilde_i <- Qtilde_i[-1] # remove Qtilde_0
  

  #####################################
  ## CASE m = 0
  #####################################
  COND_0 <- conditioning(graph = graph, alpha = 1)
  index.obs_0 <- gives.indices(graph = graph, factor = 2, constant = 2)
  nc_0 <- 1:length(COND_0$S) # number of constraints
  T_0 <- COND_0$T # change of basis matrix
  W_0 <- Diagonal(2*floor(alpha)*graph$nE)[,-nc_0] # matrix to remove constraints
  Qtilde_0_star_UU <- t(W_0) %*% t(T_0) %*% Qtilde_0 %*% (T_0) %*% W_0 
  A0 <- T_0[index.obs_0, -nc_0] # observation matrix after conditioning
  
  Qtilde_0_star_UU0 <- MetricGraph:::Qalpha1(theta = c(tau, kappa), 
                                            graph = graph, 
                                            BC = 3, 
                                            build = TRUE)*kappa^(-1)*c_1/(2 * k * c_alpha * sigma^2 * tau^2)
  A00 <- graph$.__enclos_env__$private$A()
  S0 <- A0 %*% solve(Qtilde_0_star_UU, t(A0))
  S00 <- A00 %*% solve(Qtilde_0_star_UU0, t(A00))
  print(S00/S0)
  #####################################
  ## CASE m > 0
  #####################################
  graph$buildC(alpha = 2, edge_constraint = TRUE)
  COND_i <- graph$CoB
  index.obs_i <- gives.indices(graph = graph, factor = 4, constant = 3)
  n_const <- length(COND_i$S)
  ind.const <- c(1:n_const)
  Tc <- COND_i$T[-ind.const, ]
  Qtilde_i_star_UU <- lapply(Qtilde_i, function(Q) Tc %*% Q %*% t(Tc)) 
  Ai <- t(Tc)[index.obs_i, ] # observation matrix after conditioning
  
  #####################################
  ## Build matrix A and Q_UU
  #####################################
  A <- cbind(A0, do.call(cbind, rep(list(Ai), m)))
  Q_UU <- bdiag(Qtilde_0_star_UU, do.call(bdiag, Qtilde_i_star_UU))
  # Return Sigma
  Sigma <- A %*% solve(Q_UU, t(A)) 
  return(Sigma)
}

# before I changed the constants

gets_cov_mat_rat_approx_alpha_1_to_2_old <- function(graph, kappa, tau, alpha, m){

  # get rational approximation coefficients
  coeff <- rSPDE:::interp_rational_coefficients(
    order = m, 
    type_rational_approx = "chebfun", 
    type_interp = "spline", 
    alpha = alpha)
  
  r <- coeff$r
  p <- coeff$p
  k <- coeff$k
  
  # compute parameters
  fa <- floor(alpha)
  ca <- ceiling(alpha)
  
  nu <- alpha - 1/2
  sigma <- sqrt(gamma(nu) / (tau^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))
  c_alpha <- gamma(alpha)/gamma(alpha - 0.5)
  c_1 <- gamma(fa)/gamma(fa - 0.5)
  
  # get edge lengths
  L_e <- graph$edge_lengths
  
  # initialize Qtilde_i, a list containing block diagonal matrices with blocks Qtilde_{i,e} for each i
  Qtilde_i <- list() 
  for(i in 1:m){
    
    # compute r_(0,0)
    r00 <- matern.p.joint(
      s = 0, 
      t = 0, 
      kappa = kappa, 
      p = p[i], 
      alpha = alpha)
    
    # compute r_(0,0)^(-1)
    r00_inverse <- solve(r00, Diagonal(ca))
    
    # define zero block 
    zero_block <- matrix(0, ca, ca)
    
    # build correction term
    correction_term <- rbind(
      cbind(r00_inverse, zero_block),
      cbind(zero_block, r00_inverse))
    
    # initialize Qtilde_i[[i]], a list containing Qtilde_{i,e} for each edge e
    Qtilde_i[[i]] <- list()
    for(e in 1:length(L_e)){
      
      # compute Q_{i,e}
      Q_e <- matern.p.precision(
        loc = c(0, L_e[e]),
        kappa = kappa, 
        p = p[i],
        equally_spaced = FALSE, 
        alpha = alpha)$Q
      
      # store Qtilde_{i,e}
      Qtilde_i[[i]][[e]] <- Q_e - 0.5 * correction_term
    }
    # build block diagonal matrix Qtilde_i[[i]]
    Qtilde_i[[i]] <- bdiag(Qtilde_i[[i]])
  }
  
  # --------------------------------------------------
  # CASE i = 0
  # --------------------------------------------------
  
  factor_0 <- c_1/(2 * k * c_alpha * kappa * sigma^2 * tau^2)

  
  Qtilde_0_star_UU <- MetricGraph:::Qalpha1(
    theta = c(tau, kappa), 
    graph = graph, 
    BC = 3000, 
    build = TRUE) * factor_0
  
  A_0 <- graph$.__enclos_env__$private$A()

  # --------------------------------------------------
  # CASE i = 1,...,m
  # --------------------------------------------------
  
  # build conditioning matrix
  graph$buildC(alpha = 2, edge_constraint = TRUE) # should always be TRUE
  COND_i <- graph$CoB
  Tc <- COND_i$T[-c(1:length(COND_i$S)), ]
  
  factor_i <- (2 * kappa^(2 * alpha - 1) * c_alpha * sqrt(pi) * tau^2)/r
  
  Qtilde_i_star_UU <- purrr::map2(
    Qtilde_i, 
    factor_i, 
    function(Q, x) Tc %*% Q %*% t(Tc) * x)

  index.obs_i <- gives.indices(graph = graph, factor = 4, constant = 3)
  A_i <- t(Tc)[index.obs_i, ] 
  
  # Build matrix A and Q_UU
  A <- cbind(A_0, do.call(cbind, rep(list(A_i), m)))
  Q_UU <- bdiag(Qtilde_0_star_UU, bdiag(Qtilde_i_star_UU))
  # Return Sigma
  Sigma <- A %*% solve(Q_UU, t(A)) 
  return(Sigma)
}

7 References

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

We used R version 4.5.2 (R Core Team 2025) and the following R packages: ggmap v. 4.0.2 (Kahle and Wickham 2013), gridExtra v. 2.3 (Auguie 2017), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), knitr v. 1.50 (Xie 2014, 2015, 2025), latex2exp v. 0.9.8 (Meschiari 2026), Matrix v. 1.7.3 (Bates, Maechler, and Jagan 2025), MetricGraph v. 1.5.0.9000 (Bolin, Simas, and Wallin 2023a, 2023b, 2024, 2025; Bolin et al. 2024), patchwork v. 1.3.1 (Pedersen 2025), plotly v. 4.11.0 (Sievert 2020), qs v. 0.27.3 (Ching 2025), rdist v. 0.0.5 (Blaser 2020), renv v. 1.1.7 (Ushey and Wickham 2026), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.30 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2025), rSPDE v. 2.5.2.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.4.0 (Wickham, Pedersen, and Seidel 2025), sf v. 1.1.0 (E. Pebesma 2018; E. Pebesma and Bivand 2023), slackr v. 3.4.0 (Kaye et al. 2025), sp v. 2.2.1 (E. J. Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), tidyverse v. 2.0.0 (Wickham et al. 2019), viridis v. 0.6.5 (Garnier et al. 2024), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024), xtable v. 1.8.8 (Dahl et al. 2026).

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

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

Auguie, Baptiste. 2017. gridExtra: Miscellaneous Functions for “Grid” Graphics.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.package.Matrix.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. https://asdar-book.org/.

Blaser, Nello. 2020. rdist: Calculate Pairwise Distances. https://doi.org/10.32614/CRAN.package.rdist.

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. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023b. “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.

———. 2025. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” Bernoulli. https://doi.org/10.48550/arXiv.2304.03190.

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.

Bolin, David, Alexandre Simas, and Jonas Wallin. 2023c. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” https://arxiv.org/abs/2304.10372.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://github.com/rstudio/htmltools.

Ching, Travers. 2025. qs: Quick Serialization of R Objects. https://doi.org/10.32614/CRAN.package.qs.

Dahl, David B., David Scott, Charles Roosen, Arni Magnusson, and Jonathan Swinton. 2026. xtable: Export Tables to LaTeX or HTML. https://doi.org/10.32614/CRAN.package.xtable.

Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2024. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4679423.

Kahle, David, and Hadley Wickham. 2013. “ggmap: Spatial Visualization with Ggplot2.” The R Journal 5 (1): 144–61. https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf.

Kaye, Matt, Bob Rudis, Andrie de Vries, and Jonathan Sidi. 2025. slackr: Send Messages, Images, r Objects and Files to “Slack” Channels/Users. https://github.com/mrkaye97/slackr.

Meschiari, Stefano. 2026. Latex2exp: Use LaTeX Expressions in Plots. https://doi.org/10.32614/CRAN.package.latex2exp.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://doi.org/10.32614/CRAN.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 J., and Roger Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.

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. 2025. patchwork: The Composer of Plots. https://doi.org/10.32614/CRAN.package.patchwork.

R Core Team. 2025. 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.

Ushey, Kevin, and Hadley Wickham. 2026. renv: Project Environments. https://doi.org/10.32614/CRAN.package.renv.

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. 2025. scales: Scale Functions for Visualization. https://scales.r-lib.org.

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/.

———. 2025. 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.