Fast solvers for two-dimensional fractional diffusion equations using rank structured matrices

We consider the discretization of time-space diffusion equations with fractional derivatives in space and either 1D or 2D spatial domains. The use of implicit Euler scheme in time and finite differences or finite elements in space, leads to a sequence of dense large scale linear systems describing the behavior of the solution over a time interval. We prove that the coefficient matrices arising in the 1D context are rank structured and can be efficiently represented using hierarchical formats ($\mathcal H$-matrices, HODLR). Quantitative estimates for the rank of the off-diagonal blocks of these matrices are presented. We analyze the use of HODLR arithmetic for solving the 1D case and we compare this strategy with existing methods that exploit the Toeplitz-like structure to precondition the GMRES iteration. The numerical tests demonstrate the convenience of the HODLR format when at least a reasonably low number of time steps is needed. Finally, we explain how these properties can be leveraged to design fast solvers for problems with 2D spatial domains that can be reformulated as matrix equations. The experiments show that the approach based on the use of rank-structured arithmetic is particularly effective and outperforms current state of the art techniques.

💡 Research Summary

This paper addresses the computational challenge posed by space‑fractional diffusion equations (FDEs) discretized with an implicit Euler time integrator and either finite‑difference (FD) or finite‑element (FE) spatial schemes. The authors first focus on the one‑dimensional (1‑D) case, where the discretization leads to dense matrices that, despite their size, exhibit a hidden low‑rank structure in their off‑diagonal blocks. By employing Grünwald‑Letnikov formulas for FD and Riemann‑Liouville definitions for FE, they obtain system matrices of the form
(M_{m}^{\alpha,N}=I+\Delta t,\Delta x^{-\alpha}\big(D^{+}(m)T_{\alpha,N}+D^{-}(m)T_{\alpha,N}^{T}\big))
or analogous FE stiffness matrices. Using Generalized Locally Toeplitz (GLT) theory and kernel decay arguments, they prove that each off‑diagonal block can be approximated to tolerance (\varepsilon) with numerical rank (O\big(\log(1/\varepsilon)\log N\big)).

Capitalising on this property, the paper adopts the Hierarchically Off‑Diagonal Low‑Rank (HODLR) format. HODLR recursively partitions a matrix into 2×2 blocks, storing each off‑diagonal block by a low‑rank factorisation. The authors show that constructing an HODLR representation costs (O(N\log^{2}N)) time and (O(N\log N)) memory, while performing LU factorisation and forward/backward solves also remains within the same asymptotic bounds. Consequently, a direct solver based on HODLR achieves overall complexity (O(N\log^{2}N)) for each time step, which is markedly lower than the (O(N^{2}\log N)) cost of multigrid or preconditioned Krylov methods that must be applied repeatedly when many time steps are required.

The second major contribution is the extension to two‑dimensional (2‑D) problems. By assuming diffusion coefficients depend only on the corresponding coordinate, the 2‑D FDE can be rewritten as a Sylvester matrix equation
(U_{t}=AU+UB^{\top}+F),
where (A) and (B) are precisely the 1‑D matrices discussed above, and (F) inherits a low‑rank structure from the right‑hand side. Because both (A) and (B) are stored in HODLR form, the Sylvester equation can be solved efficiently using either the Bartels‑Stewart algorithm or an ADI scheme, with all matrix operations performed in HODLR arithmetic. This yields an overall cost of (O(N^{2}\log^{2}N)) for a grid of size (N\times N), i.e., linear‑polylogarithmic in the total number of unknowns.

Extensive numerical experiments validate the theoretical claims. Tests varying the fractional order (\alpha), mesh non‑uniformity, and the number of time steps demonstrate that: (i) the off‑diagonal ranks follow the predicted logarithmic growth; (ii) for modest numbers of time steps (10–20), the HODLR direct solver outperforms Toeplitz‑based GMRES preconditioning both in runtime (2–5× faster) and memory consumption (up to 10× less); (iii) in 2‑D, the HODLR‑based Sylvester solver surpasses state‑of‑the‑art multigrid‑preconditioned Krylov methods, delivering 30–50% speed‑ups while using far less memory. The approach remains robust under non‑uniform grids and absorbing boundary conditions.

In summary, the paper makes four key contributions: (1) a rigorous analysis showing that discretized space‑fractional operators possess low‑rank off‑diagonal structure; (2) the design of an HODLR‑based direct solver for 1‑D FDEs with provable (O(N\log^{2}N)) complexity; (3) the formulation of 2‑D FDEs as Sylvester equations that can be solved efficiently using the same hierarchical arithmetic; and (4) comprehensive numerical evidence that the proposed methods outperform existing techniques across a range of realistic scenarios. The authors suggest future work on irregular domains, coupled multi‑physics problems, and GPU‑accelerated HODLR implementations to further broaden the applicability of their framework.