Data-driven accelerated multigrid method for Riesz space fractional diffusion equations

Introduction

Fractional differential equations are powerful mathematical tools that can accurately describe certain complex phenomena, such as pricing options [1, 17], control problems [10, 23], image processing [4, 2]. Consider a two-dimensional Riesz space fractional diffusion equation as follows:

where $\Omega=(x_L,x_R)\times(y_L,y_R)$, $\partial\Omega$ denotes the boundary of $\Omega$, $f(x,y,t)$ is the source term, $d(x,y,t)\ge0$ is the coefficient function. $u_0(x,y)$ is the initial condition, $\gamma_1,\gamma_2\in(1,2)$. Specifically, Riesz fractional derivatives are defined as follows:

with $\sigma_{\gamma_l}=-\frac{1}{2\cos(\pi\gamma_l/2)}>0$ for $l=1,2$. Owing to the nonlocality of fractional differential operators, it is difficult to obtain the analytical solutions of fractional differential equations, and thus developing efficient numerical methods to acquire their numerical solutions has become a research focus in this field.

Discretizing equation (1) on a uniform grid using the shifted Grünwald difference scheme in space and the backward Euler scheme in time yields a sequence of Toeplitz-like linear systems. Due to the nonlocality of fractional differential operators, the coefficient matrices are generally dense, which renders the discrete linear systems computationally expensive for general purpose solvers. For solving such linear systems, [5] proposed a class of diagonal–Toeplitz splitting (DTS) iteration methods, implemented by alternately solving shifted diagonal and shifted Toeplitz subsystems. Building upon the DTS framework, [16] used a natural $\tau$ approximation of the symmetric Toeplitz component and proposed a sine-transform-based splitting preconditioner. More recently, [26] decomposed the coefficient matrix by extracting a scalar matrix, a diagonal compensation matrix, and a block $\tau$ matrix. The sum of these matrices was then employed to construct a class of $\tau$-splitting iteration methods.

For certain structured matrices, the multigrid method exhibits excellent solution efficiency[24, 15, 6]. Under suitable conditions, this method enjoys mesh-independent convergence. This advantageous property has drawn the attention of many researchers to investigate the application of multigrid methods to fractional partial differential equations. For instance, [21] employed the two-grid method to solve the dense Toeplitz linear system obtained after discretizing the one-dimensional constant-coefficient fractional diffusion equation, and utilized fast Toeplitz matrix–vector multiplication to reduce the computational overhead. [9] proved the convergence of the multigrid method for two-dimensional space-fractional diffusion equations with constant diffusion coefficients. Based on the generalized locally Toeplitz (GLT) sequences theory, they further extended this result to the case of uniformly bounded and positive variable diffusion coefficients. Furthermore, numerous studies adopt the GLT theory to analyze the spectral properties of coefficient matrices and thereby design multigrid methods[2, 25, 3]. In addition, [7, 8] extended the symbol-based multigrid analysis for scalar Toeplitz matrices to block-circulant and block-Toeplitz systems generated by matrix-valued symbols, and derived sufficient conditions for the two-grid and V-cycle methods.

The multigrid method consists of a smoother and a coarse-grid correction that complement each other. The coarse-grid correction is responsible for eliminating the low-frequency components of the error. For Toeplitz-like matrices, coarse-grid correction using bilinear interpolation as the transfer operator is effective and general [20]. Correspondingly, the smoother should rapidly attenuate the high-frequency components of the error that are difficult for the coarse grid to represent, and it is usually implemented by classical iterative methods. The weighted Jacobi method is widely adopted due to its simplicity and low computational cost per iteration. With an appropriate relaxation parameter, this method exhibits favorable smoothing ability for high-frequency components of the error[30, 32]. Therefore, the proper choice of the relaxation parameter is the key to constructing an efficient smoother.

Previously, research on the selection of relaxation parameters for the weighted Jacobi smoother has mainly focused on classical elliptic partial differential equations[12, 13]. [3]extended this study to dense Toeplitz-type linear systems arising from the discretization of tempered fractional diffusion equations. They used the spectral information provided by the generating function of the coefficient matrix to expand the range of parameters that guarantee the smoothing property of the Jacobi smoother, and gave estimation formulas for suitable relaxation parameters. However, while their estimation formulas yield feasible relaxation parameters, the resulting parameter may not always match the numerically optimal choice for finite-dimensional linear systems.

In recent years, the integration of machine learning with traditional numerical methods has demonstrated great potential. Machine learning techniques can identify nonlinear relationships between problem characteristics and numerical algorithms from data that are difficult to characterize explicitly, thereby providing new approaches to the design of complex numerical algorithms and the automatic selection of their parameters. Moreover, once trained, machine learning models can typically make rapid predictions for new problems at a low computational cost. In the field of multigrid methods, [11, 33, 31] have applied machine learning to the design and optimization of smoothers. In this work, we adopt Gaussian Process Regression (GPR) to learn and select the relaxation parameters of the weighted Jacobi smoother, so as to improve the convergence speed of the multigrid method.

GPR is a Bayesian nonparametric learning method that places a Gaussian process prior on the unknown mapping from problem features to the desired relaxation parameter. Given a small set of training samples, the correlations among the samples are characterized by a covariance kernel, and conditioning the prior on the observed data yields a posterior predictive distribution for a new problem instance. Its mean can be used as the predicted relaxation parameter, while its variance provides a quantitative measure of prediction uncertainty. It can capture nonlinear parameter dependence without prescribing a fixed functional form, perform well with relatively small datasets, and, once trained, deliver predictions at a low online cost.

The paper is organized as follows. In Section 2,we discretize the original equations and analyze the basic spectral properties of the coefficient matrix. In Section 3, with the aid of the Ruge–Stuben two grid framework, we design and analyze the convergence of the two grid method. In Section 4, we introduce the principles of GPR and its procedure for predicting parameters. In Section 5, we present numerical experiments to demonstrate the effectiveness and superiority of the multigrid method and GPR for parameter prediction. In Section 6, we summarize the conclusions and future work.

Discretization and matrix properties

In this section, we will introduce the specific discrete schemes and analyze the properties of the coefficient matrices obtained from the discretization, which will help us design the multigrid method. Below, we review the left- and right-sided Riemann–Liouville fractional derivatives:

where $\Gamma(\cdot)$ is the Gamma function. Let $N$, $M$ be positive integers. Take the step sizes $\Delta x=\frac{x_R-x_L}{N+1}$, $\Delta y=\frac{y_R-y_L}{N+1}$, $\Delta t=\frac{T}{M}$ to discretize the $x$-direction, $y$-direction and time, yielding

Define $u_{i,j}^{(m)}=u(x_i,y_j,t_m)$, $f_{i,j}^{(m)}=f(x_i,y_j,t_m)$, $d_{i,j}^{(m)}=d(x_i,y_j,t_m)$. Using the shifted Grünwald formula[19, 18], we can obtain the approximation of the fractional derivative at the point $(x_i,y_j,t_m)$:

and

where the coefficients $g_k^{(\gamma)}$ are given by the following lemma:

To facilitate the subsequent analysis, we introduce a matrix $T_N^\gamma\in\mathbb{R}^{N\times N}$ associated with $\gamma$:

According to Lemma 2.1, $T_N^\gamma$ is a matrix that is strictly diagonally dominant, symmetric, and has all positive diagonal entries. Therefore, $T_N^\gamma$ is a symmetric positive definite matrix. Let

We discretize (1.1) in the spatial direction using the shifted Grunwald formula and in the temporal direction using the backward Euler scheme to obtain the implicit finite-difference scheme:

Furthermore, we have

Let $\eta_1=\sigma_{\gamma_1}\frac{\Delta t}{\Delta x^{\gamma_1}}$, $\eta_2=\sigma_{\gamma_2}\frac{\Delta t}{\Delta y^{\gamma_2}}$. Consequently, we obtain

Dividing both sides of (2.2) by $\eta_1$, we have

Define

Then (2.2) is equivalent to

After fixing the partition number $N$, discretizing (1.1) yields a set of matrices $\{A_N^{(m)}\}$. When the coefficient $d(x,y,t)=0$, $D^{(m)}$ degenerates to the zero matrix, and in this case, $A_N^{(m)}$ is independent of the time level $m$ and is denoted by $A_N$. The matrix $A_N$ is symmetric positive definite and strictly diagonally dominant, and it is a Block-Toeplitz-Toeplitz-Block (BTTB) matrix.

BTTB matrices belong to the class of multilevel Toeplitz matrices, and their spectral information is closely related to the generating function of the corresponding BTTB matrix sequence. Let us now recall the definitions of Toeplitz matrix sequences and BTTB matrix sequences.

In order to obtain the generating function of $\{A_N\}_{N\in\mathbb{N}}$, we first give the following lemma.

Proof. It is noted that the generating function of $\{T_N^\gamma\}_{N\in\mathbb{N}}$ is

where the coefficients are given by

Furthermore, we obtain that

By virtue of the results in [9], we have

and

Owing to

When $0\leq\theta\leq\pi$,

When $-\pi<\theta\leq0$,

::: Combining the two gives (7). The Proof is completed.

With the above preparations, we obtain the generating function of the BTTB matrix sequence $\{A_N\}_{N\in\mathbb{N}}$.

Remark 2.1. If a bivariate function $f(\theta_1,\theta_2)$ satisfies $f(\theta_1,\theta_2)=f(\lvert\theta_1\rvert,\lvert\theta_2\rvert)$, then $f(\theta_1,\theta_2)$ is called an even function. It is easy to verify that $f_{\gamma_1,\gamma_2}(\theta_1,\theta_2)$ is a bivariate continuous real-valued even function.

For a real-valued continuous even function $f(\theta_1,\theta_2)$, the associated BTTB matrix $T_n^{(2)}(f)$ is real symmetric. To characterize the spectrum of the matrices generated by $f_{\gamma_1,\gamma_2}(\theta_1,\theta_2)$, we recall the following standard result, which bounds the eigenvalues in terms of the minimum and maximum of the generating function:

Lemma 2.3. [28] Let $f(\theta_1,\theta_2)\in C^1([-\pi,\pi]^2)$ be a real-valued even function, with its maximum and minimum denoted by $f_{\max}$ and $f_{\min}$, respectively. Then $T_n^{(2)}(f)$ is a real symmetric matrix, and its eigenvalues are bounded by the range of $f(\theta_1,\theta_2)$ in the sense that

Two Grid method

The Galerkin multigrid method belongs to the class of classical multigrid methods and is primarily applicable to symmetric positive definite linear systems $Au=b$. Its simplest form is the two grid method (TGM), which consists of a smoother and a coarse-grid correction. Multi-level schemes such as the V-cycle are realized by recursively cycling the TGM. The specific procedure of the TGM is shown in Table 1. In this paper, we choose the weighted Jacobi as the smoother, whose update formula is

1) Pre-smoothing: Apply $j_1$ iterations of $G$ to the system $Au=b$ to obtain $u^{(j_1)}$
2) Compute the residual: $r=b-Au^{(j_1)}$
3) Restriction: $b_H=I_h^H r$
4) Coarse grid equation solving: Solve the system $A_H\widehat{e}=b_H$ exactly to obtain $\widehat{e}$
5) Update solution $\widehat{u}=u^{(j_1)}+I_H^h\widehat{e}$
6) Post-smoothing: Starting from $\widehat{u}$, perform $j_2$ iterations using $G$

where $u^{(j)}$ is the iteration vector at the $j$-th step, and

The coarse grid correction is obtained by combining the restriction operator $I_h^H$ and the prolongation operator $I_H^h$, where $I_h^H$ is a rectangular matrix with full row rank. Specifically, we take

where

The prolongation operator is defined by

which defines a Galerkin coarse matrix $A_H = I_h^H A I_H^h$.

Take $j_1=0$, $j_2=1$. In this setting, the iteration matrix of the TGM is given as follows:

$E_{TG}$ is also referred to as the two grid error propagation operator. To characterize its convergence, we introduce the following inner product and the corresponding induced energy norm.

For the TGM, Ruge and Stüben established the following analytical framework, which reduces the overall convergence analysis of the method to the separate examination and characterization of the properties of the smoother and prolongation operators.

Ruge and Stüben also rigorously discussed the range of the relaxation parameter $\omega$ when the smoother $G$ is taken as the weighted Jacobi iteration method, in order to guarantee the smoothing property, namely the following lemma:

We first establish a TGM analysis for $A_N$, verifying that it satisfies both the smoothing property and the approximation property. Building upon this, by utilizing the perturbation property of positive diagonal matrices, we further derive the TGM analysis for $A_N^{(m)}$.

Proof. Since all diagonal entries of $A_N$ are $a_0$, we have

By Lemma 2.3, we obtain

Thus, from Lemma 3.2, a sufficient condition for $\omega$ is

The Proof is completed.

The approximation properties of prolongation operators were studied in [28]. In particular, sufficient conditions were established under which a BTTB matrix generated by a nonnegative even function ensures that the $I_H^h$ defined by (3.2) satisfies (3.8). More precisely, the following lemma holds.

Proof. Set

Define

Since $|\theta|\in[0,\pi]$ and $\gamma\in(1,2)$, we have

and

Then

Consequently, for every $(\theta_1,\theta_2)\ne(0,0)$,

Hence,

Therefore,

Owing to

and

Then

Thus, the quotient cannot attain its infimum near the origin. Since it is continuous away from the origin, its infimum is attained on a compact subset of $[-\pi,\pi]^2 \setminus\{(0,0)\}$. Therefore, the infimum may equivalently be written as a minimum.

The Proof is completed.

Several properties obtained by applying positive diagonal matrix perturbations are presented below. If $A$ and $B$ are symmetric positive definite matrices, the notation $A \le B$ means that $B-A$ is a nonnegative definite matrix.

The convergence of the TGM for $A_N$ has been analyzed previously. For the case $d(x,y,t)\geq0$, although $A_N^{(m)}$ may no longer be a BTTB matrix, the coefficient matrix at each time level can be expressed as $A_N^{(m)}=A_N+\frac{1}{\eta_1}D^{(m)}$, and it follows from Lemma 3.4 that

The above corollary indicates that regardless of how the conditions of equation (1) vary, as long as the relaxation parameter $\omega\in(0,1]$, the TGM is always convergent, which provides a conservative but safe interval for parameter selection. Furthermore, the surface of $\frac{2a_0}{\left\|f_{\gamma_1,\gamma_2}\right\|_\infty}$ corresponding to different values of $(\gamma_1,\gamma_2)$ when $N=64$ is shown in Figure 1. It can be observed that the surface representing the values of $\frac{2a_0}{\left\|f_{\gamma_1,\gamma_2}\right\|_\infty}$ lies strictly above the light blue plane, and $\frac{2a_0}{\left\|f_{\gamma_1,\gamma_2}\right\|_\infty}$ achieves an approximate maximum value of 1.0615 when $(\gamma_1,\gamma_2)=(1.445,1.445)$. In the following, we use an example to show that the choice of $\omega$ will significantly affect the iteration efficiency of the TGM.

and the source term $f$ is derived from the exact solution. Take $\gamma_1=1.1$, $\gamma_2=1.7$, $N=64$, $M=4$, and consider the linear system at the first time level. We apply the TGM to solve the linear system obtained from discretization, terminate the iteration when

and count the number of cycles. Specifically,

(1): The $d(x,y,t)=0$: We set $\omega=0.0001$ as the initial value and traverse to $\omega=1.1$ with a step size of 0.0001. Denoting one TGM iteration as one cycle, we record the cycle counts required to reach the stopping tolerance for different values of $\omega$, and plot the corresponding iteration curve, namely Figure 2a. A direct calculation gives

which is marked by a red line in Figure 2a.

As can be seen from Figure 2a, the choice of the relaxation parameter $\omega$ significantly affects the iteration efficiency of the TGM. For example, when $\omega=0.9335$, the required number of cycles is the smallest, only 25; whereas in other regions, the corresponding cycle counts increase markedly, even exceeding 100. In addition, all values of $\omega$ to the left of the red line ensure the convergence of the TGM, which verifies the effectiveness of Theorem 3.1.

Figure 2b shows the steep growth phenomenon near $\omega=1.0417$: when $\omega=1.04$, the number of cycles required to reach the convergence tolerance is 795; when $\omega=1.0417$, the required number of cycles increases sharply to 3043; and when $\omega=1.042$, the TGM still fails to converge after more than 5000 iterations. These phenomena indicate that the choice of $\omega$ has a great influence on both the convergence and the convergence speed of the TGM, so selecting an appropriate $\omega$ is necessary for accelerating convergence.

(2): $d(x,y,t)=(1+x)^7(1+y)^7$: Similar to the results obtained in (1), Figure 3a shows a sweep of the number of cycles required to reach convergence accuracy under different values of $\omega$. It can be seen that the number of iterations is minimized at $\omega=0.9231$, requiring only 14 cycles. In Figure 2b, the parameters to the right of $\omega=1.0417$ basically make the number of iterations required by TGM exceed 5000, whereas in Figure 3b, the parameters to the left of $\omega=1.0513$ all allow TGM to converge. It can be seen that the boundary parameter that ensures the smoothing property is shifted noticeably to the right, which verifies the validity of Corollary 3.1. In addition, when $\omega= 1.0512$, the required number of cycles is 1298, while when $\omega= 1.0513$, the iteration still fails to converge after more than 5000 steps. This steep growth phenomenon is far more pronounced than in the case of $d(x,y,t)=0$, reflecting that the parameter variation is more complicated when $d(x,y,t) \ge0$.

Data-driven acceleration strategy

As the discretization number $N$ increases, the discretization of the partial differential equation yields a family of matrices. When solving with parameterized iterative methods, there exists an optimal parameter for each problem size, which can be obtained via fine step grid search and corresponds to the minimal computational cost. As $N$ continues to grow, properties such as the matrix spectrum change, and these optimal parameters usually fluctuate. Such fluctuations generally exhibit nonlinear characteristics and cannot be characterized by rigorous mathematical expressions, determining them via exact eigenvalue computations incurs a high computational cost[29].

Gaussian process regression (GPR) provides an effective data-driven strategy for learning such parameters: it uses kernel functions to model the varying relationship between problem sizes and quasi-optimal parameters, and predicts near-optimal parameters at a low computational cost. Assume that we have a training set containing $d$ data:

where the pair $(N_i,\omega_i)$ denotes the $N_i$-th mesh partition number and the corresponding optimal relaxation parameter. The goal of Gaussian process regression is to establish a mapping between the input set $\mathbf{N}$ and the output set $\boldsymbol{\omega}$, represented as $\mathcal{F}(\mathbf{N}):\mathbb{R}\to\mathbb{R}$, and to forecast the potential output value $\tilde{\omega}_*=\mathcal{F}(N_*)$ at a new test point $N_*$. In practical linear regression scenarios, we model the observed value $\omega$ as the result of a latent function $\mathcal{F}(\mathbf{N})$ perturbed by additive noise $\xi$, i.e.,

Suppose that $\xi$ follows a Gaussian distribution with a zero mean and variance $\sigma^2$, that is $\xi\sim\mathcal{N}(0,\sigma^2)$. The common range for $\sigma$ is $[10^{-6},10^{-2}]$. For this investigation, we have set $\sigma=10^{-4}$. Without loss of generality, assume true function $\mathcal{F}(\mathbf{N})$ follows a zero mean Gaussian process:

where $\kappa(N_i,N_j)$ defines the covariance between the function values corresponding to the input points $N_i$ and $N_j$. The kernel function plays a pivotal role in GPR, as it is responsible for generating the covariance matrix that quantifies the similarity between any two input variables. The closer the distance between inputs, the higher the correlation between their corresponding output variables. Hence, it is essential to select or construct kernel functions that meet the specific demands of the problem at hand. In this paper, we consider the Gaussian kernel function, which is given by

where $\delta_{\mathcal{F}}$, $\delta_{\mathcal{L}}$ are hyperparameters,and the specific calculation of the hyperparameters is given in [14]. Define $K(\mathbf{N},\mathbf{N})=(\kappa(N_i,N_j))\in\mathbb{R}^{d\times d}$. Then, we can express the prior distribution of the observed values $\boldsymbol{\omega}$ as

The joint prior distribution of the observed value $\boldsymbol{\omega}$ and the predicted value $\omega_*$ is given by

Then the conditional expectation and variance of $\mathcal{F}(n)$ is given by

For the output $\omega_*$ in the test set, we use the mean of the Gaussian regression as its estimate, i.e. $\widehat{\omega}_*=\mathbb{E}[\mathcal{F}(N)|\boldsymbol{\omega}]$.

Numerical experiments

In Section 3, we discuss the convergence of the TGM. The TGM contains only one coarse level. When the size of the original matrix is too large, the order of the coarse matrix obtained by projection may still be high, and it remains difficult to solve the coarse-grid equation exactly. A natural approach to resolve this difficulty is to recursively apply the TGM on this coarse level to generate coarser levels successively, until the order of the coarsest level is sufficiently small. Except for the coarsest level, each of the remaining grid levels performs only one such recursive procedure, which yields the V-cycle of the multigrid method[24, 27, 32].

In this section, we present several numerical experiments to demonstrate the computational efficiency of the multigrid method and verify the acceleration effect of using GPR-predicted parameters on the multigrid solver. Specifically, we set the order of the matrix on the coarsest level of the multigrid method to be no more than 1000. Except for the coarsest level, one weighted Jacobi iteration is applied for pre-smoothing and post-smoothing on all remaining grid levels, and the same input relaxation parameter is used for all weighted Jacobi iterations. In the experiments, we consider three test cases for comparison: the multigrid method with $\omega= 1$ as the relaxation parameter, the multigrid method using GPR-predicted values as the relaxation parameter, and the unpreconditioned conjugate gradient method.

The training strategy of GPR is as follows: starting from $N = 5$ with a step size of 5 up to $N = 50$, we construct a training set consisting of 10 samples. For each problem size, the quasi-optimal parameter $\omega$ is obtained via grid search over the interval $[0.001, 1]$ with a step size of 0.001. For all algorithms, the zero vector is used as the initial iterate, with the stopping criterion

All experiments are carried out using MATLAB (version R2025b) on a laptop equipped with an Intel Core i5-9300H CPU (4 cores, 2.40 GHz) and 8 GB of memory.

Example 5.1. Consider the two-dimensional Riesz space fractional diffusion equation (1) defined on the domain $(0, 1) \times(0, 1) \times(0, 1)$. The exact solution is given by $u(x, y, t) = 2(t^3 + 0.01)xy(1 - x)(1 - y)$, and the initial condition is $u_0(x, y) = 0.02xy(1 - x)(1 - y)$. We take $d(x, y, t) = \ln(1 + x)(y + 1)^2$. The source term $f$ is derived from the exact solution.

It can be seen that the GPR parameter selection strategy is significantly superior to the Default strategy in most cases. For example, in Table 2 for $(\gamma_1,\gamma_2)=(1.4,1.4)$ and $N=512$, Default requires 18 iterations and 20.7812 s, whereas GPR requires only 14 iterations and 15.3594 s; in Table 2 for $(\gamma_1,\gamma_2)=(1.8,1.8)$ and $N=512$, Default requires 18 iterations and 20.4062 s, whereas GPR requires only 14 iterations and 15.9375 s. As $N$ increases, the iteration counts of GPR remain basically in the range of 13–31, exhibiting a favorable mesh-independent convergence trend; although the iteration counts of Default decrease with increasing $N$ in some test cases, they are still overall higher than those of GPR. The iteration counts and CPU times of CG grow fastest with $N$; at $N=512$, CG usually requires hundreds of iterations and hundreds of seconds of CPU time, which is far higher than those of the two multigrid methods. For a few parameter combinations, the CPU times of Default and GPR are already relatively close. However, over all test cases, the relaxation parameters predicted by GPR can effectively reduce both the iteration counts and the computing time.

$(\gamma_1,\gamma_2)$	Method	Index	$N=64$	$N=128$	$N=256$	$N=512$
$(1.4,1.4)$	Default	IT	27	21	19	18
		CPU	0.1710	1.0781	3.7969	20.7812
	GPR	IT	13	14	14	14
		CPU	0.0625	0.3596	2.1875	15.3594
	CG	IT	97	161	278	447
		CPU	0.2969	2.3125	32.7656	440.1875
$(1.5,1.6)$	Default	IT	30	24	20	17
		CPU	0.2656	1.1562	3.5781	21.0781
	GPR	IT	16	17	18	19
		CPU	0.1094	0.5156	2.8212	21.6406
	CG	IT	111	199	345	585
		CPU	0.1875	2.5940	42.6407	552.3438
$(1.7,1.9)$	Default	IT	72	60	44	33
		CPU	0.3906	2.0469	7.2656	39.1406
	GPR	IT	23	26	30	31
		CPU	0.0938	0.8125	4.3975	37.4217
	CG	IT	155	303	567	1072
		CPU	0.2812	3.7969	67.3438	1067.7187

Example 5.2. Consider the two-dimensional Riesz space fractional diffusion equation (1) defined on the domain $(0,1)\times(0,1)\times(0,1)$. The exact solution is given by $u(x,y,t)=e^{-2t}\sin(\pi x)\sin(\pi y)$, and the initial condition is $u_0(x,y)=\sin(\pi x)\sin(\pi y)$. We take $d(x,y,t)=10^{-6}$. The source term $f$ is derived from the exact solution.

$(\gamma_1,\gamma_2)$	Method	Index	$N=64$	$N=128$	$N=256$	$N=512$
$(1.2,1.2)$	Default	IT	39	29	22	21
		CPU	0.1875	1.2812	4.0469	23.8595
	GPR	IT	13	14	14	14
		CPU	0.0781	0.5312	2.8906	17.0781
	CG	IT	71	117	189	302
		CPU	0.2031	1.3906	22.8125	397.6406
$(1.5,1.5)$	Default	IT	30	24	17	16
		CPU	0.3281	1.0625	3.2188	18.4712
	GPR	IT	13	13	13	14
		CPU	0.1406	0.3125	2.4219	16.5000
	CG	IT	72	125	215	369
		CPU	0.1406	1.5469	24.6250	419.6719
$(1.8,1.8)$	Default	IT	59	49	34	18
		CPU	0.2970	1.7500	6.0938	20.4062
	GPR	IT	14	14	14	14
		CPU	0.0781	0.3906	2.7031	15.9375
	CG	IT	83	161	302	568
		CPU	0.2031	2.1406	36.0000	688.1875

It can be seen that the GPR parameter selection strategy is overall significantly superior to Default. For example, for $(\gamma_1,\gamma_2)=(1.2,1.2)$ and $N=512$, Default requires 21 iterations and 23.8595 s, whereas GPR requires only 14 iterations and 17.0781 s; for $(\gamma_1,\gamma_2)=(1.8,1.8)$ and $N=512$, Default requires 18 iterations and 20.4062 s, whereas GPR requires only 14 iterations and 15.9375 s. As $N$ increases, the iteration counts of GPR remain basically stable at 13–14, exhibiting a favorable mesh-independent convergence trend; although the iteration counts of Default also decrease with increasing $N$, they are still overall higher than those of GPR, and the gap is particularly evident at $N=64$ and $N=128$. The iteration counts and CPU times of CG grow fastest with $N$; at $N=512$, CG usually requires hundreds of iterations and hundreds of seconds of CPU time. For instance, for $(\gamma_1,\gamma_2)=(1.8,1.8)$, CG requires 568 iterations and 688.1875 s, which is far higher than those of the two multigrid methods. For a few parameter combinations (such as $(\gamma_1,\gamma_2)=(1.5,1.5)$ and $N=512$), the CPU times of Default and GPR are already relatively close, but GPR is still no worse than Default. Over all test cases, the relaxation parameters predicted by GPR can effectively reduce both the iteration counts and the computing time, further verifying the effectiveness and superiority of the data-driven parameter selection strategy for accelerating the multigrid method.

Conclusion

In this paper, the weighted Jacobi method is adopted as the smoother in a multigrid method for solving the two-dimensional Riesz space fractional diffusion equation. To address the problem of relaxation parameter selection, a GPR-based parameter selection strategy is proposed. Numerical experiments verify the correctness of the relevant theory and confirm the high computational efficiency of this data-driven algorithm. Moreover, the combination of neural networks and traditional numerical methods has become increasingly important in recent years. In fact, neural networks exhibit spectral bias, which makes them naturally complementary to multigrid methods. As a result, a growing number of researchers have combined the two and obtained many important results. Nevertheless, reliable mathematical explanations for these works remain lacking. In our subsequent work, we will provide corresponding analysis and design more efficient and robust algorithms.

Data availability: No data was used for the research described in the article.

Acknowledgements: The authors are thankful to the editors and the anonymous referees for their valuable comments to improve the paper.

Funding: This work was supported by the National Natural Science Foundation of China (12671450), the National Key Research and Development Program of China (2023YFB3001604), and the Graduate Research Innovation Project of Xiangtan University (XDCX2026Y276).