Efficient Krylov solvers for inverse source problem in 2D space-time fractional diffusion equation

Technical Summary: Efficient Krylov Solvers for Inverse Source Problems in 2D Space-Time Fractional Diffusion Equations

Problem Statement
This work addresses the inverse source problem associated with a two-dimensional space-time fractional diffusion equation (FDE) featuring variable coefficients. The governing equation involves a tempered fractional derivative in time ($T D^{\xi,\rho}_t$) and Caputo fractional derivatives in two spatial dimensions ($CD^{\eta_1}_x, CD^{\eta_2}_y$). The objective is to reconstruct the unknown source term $f(x, y)$ given initial conditions and noisy final-time data.

The problem is inherently ill-posed in the sense of Hadamard, where small perturbations in the input data (noise level $\epsilon$) can lead to significant errors in the solution. To mitigate this, the authors employ the quasi-boundary value method, a nonlocal regularization technique that transforms the ill-posed problem into a parameter-dependent well-posed regularized problem.

Methodology
The proposed approach combines numerical discretization, spectral analysis via Generalized Locally Toeplitz (GLT) theory, and the design of efficient preconditioners for iterative solvers.

1. Discretization:

   • The regularized problem is discretized using a finite difference scheme. The tempered fractional time derivative is approximated using the L1 formula, while the spatial Caputo derivatives are also discretized via L1 approximations.
   • This process yields a large-scale "all-at-once" linear system of the form $A_{\tilde{N}} v_\lambda = z_\lambda$, where $\tilde{N} = (m+1)\nu(n)$ represents the total degrees of freedom (combining time steps $m$ and spatial grid points $\nu(n) = n_1 n_2$).
   • The resulting coefficient matrix $A_{\tilde{N}}$ exhibits a complex multilevel block structure. Due to the variable coefficient $a(x, y)$, the matrix is not strictly Toeplitz but possesses a multilevel block Toeplitz-like structure perturbed by diagonal sampling matrices.

2. Spectral Analysis (GLT Theory):

   • The authors perform an asymptotic spectral analysis of the matrix sequences using Generalized Locally Toeplitz (GLT) theory. This framework allows for the characterization of the singular value and eigenvalue distributions of large matrix sequences arising from discretized differential operators.
   • The analysis identifies the GLT symbol of the coefficient matrix sequence. The symbol is derived as the sum of the symbol associated with the time-stepping block (related to the tempered derivative) and the spatial operator block (involving the variable coefficient and spatial fractional derivatives).
   • The analysis confirms that the matrix sequence is distributed according to this symbol, providing the theoretical basis for preconditioning.

3. Preconditioning Strategy:

   • Based on the GLT analysis, the authors propose a block triangular preconditioner ($P_{\tilde{N}}$). This preconditioner approximates the original matrix $A_{\tilde{N}}$ by retaining the block lower triangular structure but simplifying the diagonal blocks.
   • Theoretical results (Theorem 5.2) demonstrate that the preconditioned matrix sequence $\{P_{\tilde{N}}^{-1} A_{\tilde{N}}\}$ is distributed as the identity function (symbol = 1) in the GLT sense. Consequently, the eigenvalues of the preconditioned system are weakly clustered around 1, which is a desirable property for the convergence of Krylov subspace methods like GMRES.
   • The solution of linear systems involving the preconditioner's diagonal blocks is handled via standard techniques for lower triangular Toeplitz-like systems (e.g., using circulant approximations or integral averages of the variable coefficient).

Key Contributions

• Formulation: The paper formulates the inverse source problem for a 2D space-time FDE with variable coefficients and applies the quasi-boundary value method to ensure stability against noise.
• Structural Analysis: It provides a rigorous spectral analysis of the resulting "all-at-once" discretization matrices using GLT theory, explicitly deriving the symbol for the multilevel block structure.
• Preconditioner Design: A novel block triangular preconditioner is constructed that preserves the multilevel structure of the discretization. The authors prove that this preconditioner leads to eigenvalue clustering around unity, theoretically guaranteeing efficient convergence.
• Computational Efficiency: The proposed method utilizes the matrix structure to reduce the computational cost of matrix-vector products and preconditioner applications to $O(m^2 \nu(n))$ operations (with logarithmic factors for Toeplitz solves), making it scalable for large grids.

Numerical Results
Numerical experiments were conducted on the domain $\Omega = (0, 1)^2$ with various fractional orders ($\xi, \eta_1, \eta_2$) and grid sizes.

• Convergence: The preconditioned GMRES method significantly outperformed the unpreconditioned version. For instance, with grid sizes $n_1=n_2=m=25$, the iteration count dropped from 186 (unpreconditioned) to 39 (preconditioned) for specific fractional parameters, representing an ~80% reduction.
• Robustness: The number of iterations for the preconditioned system remained nearly constant as the temporal discretization ($m$) increased and grew only moderately as spatial discretization ($n_1, n_2$) increased. In contrast, the unpreconditioned system showed rapid growth in iteration counts with mesh refinement.
• Cost: While the application of the preconditioner adds overhead, the total CPU time was substantially reduced due to the drastic decrease in the number of GMRES iterations required for convergence.

Significance and Claims
The paper claims that the proposed framework offers an efficient and scalable approach for solving inverse source problems in fractional diffusion models. By integrating regularization, structure-preserving discretization, and GLT-based spectral analysis, the method successfully accelerates the GMRES solver.

The authors emphasize that the preconditioner effectively mitigates the deterioration of spectral properties typically associated with mesh refinement in fractional problems. They note that while the clustering is "weak" (as is typical for multilevel problems), the resulting convergence behavior is robust across different fractional orders and discretization levels. The work extends previous unidimensional studies to the more complex two-dimensional case, demonstrating that advanced matrix structure exploitation can overcome the computational challenges posed by non-local fractional operators and ill-posed inverse problems.

The paper modestly suggests that future work could extend these methods to more complex geometries, higher-order finite elements, isogeometric analysis, and higher spatial dimensions, noting that the theoretical framework presented is adaptable to these scenarios.