[Paper Review] $\mathcal{H}$-matrix approximability of inverses of discretizations of the fractional Laplacian

Michael Karkulik, Jens Markus Melenk|arXiv (Cornell University)|Aug 13, 2018

Advanced Mathematical Modeling in Engineering1 references3 citations

TL;DR

This paper proves that the inverse of the stiffness matrix from a Galerkin discretization of the fractional Laplacian on a bounded domain can be approximated by hierarchical matrices ($\mathcal{H}$-matrices) at an exponential convergence rate in the block rank. The method relies on the Caffarelli-Silvestre extension and blockwise low-rank approximation via singular value decomposition, establishing theoretical and numerical evidence for efficient solution of fractional PDEs using $\mathcal{H}$-matrix arithmetic.

ABSTRACT

The integral version of the fractional Laplacian on a bounded domain is discretized by a Galerkin approximation based on piecewise linear functions on a quasi-uniform mesh. We show that the inverse of the associated stiffness matrix can be approximated by blockwise low-rank matrices at an exponential rate in the block rank.

Motivation & Objective

• To establish the theoretical approximability of the inverse of the stiffness matrix from a Galerkin discretization of the fractional Laplacian using $\mathcal{H}$-matrices.
• To extend the applicability of $\mathcal{H}$-matrix techniques—previously used for classical boundary element matrices—to the non-local, fully populated systems arising from fractional differential equations.
• To provide a foundation for efficient iterative solvers and preconditioners in the context of fractional PDEs through $\mathcal{H}$-matrix compression of the inverse.
• To validate the theoretical findings numerically across different domains and fractional orders, demonstrating exponential convergence in practice.

Proposed method

• Discretization of the integral form of the fractional Laplacian using piecewise linear finite elements on a quasi-uniform triangulation.
• Application of the Caffarelli-Silvestre extension to represent the fractional Laplacian as a Dirichlet-to-Neumann map of a degenerate elliptic PDE in $\mathbb{R}^{d+1}$.
• Construction of a cluster tree and block cluster tree using geometric clustering with admissibility parameter $\eta=2$.
• Computation of blockwise low-rank approximations of the inverse matrix via singular value decomposition on admissible blocks.
• Use of standard $\mathcal{H}$-matrix structure and admissibility criterion (2.11) to ensure data-sparse representation.
• Numerical experiments on square and L-shaped domains with varying fractional orders $s \in \{0.25, 0.5, 0.75\}$ and mesh sizes to validate the exponential convergence rate.

Experimental results

Research questions

• RQ1 Can the inverse of the stiffness matrix arising from a Galerkin discretization of the fractional Laplacian be approximated by $\mathcal{H}$-matrices with exponential convergence in block rank?
• RQ2 Does the Caffarelli-Silvestre extension enable the theoretical justification of $\mathcal{H}$-matrix approximability for the inverse of the fractional Laplacian?
• RQ3 How does the convergence rate of the $\mathcal{H}$-matrix approximation of the inverse compare to the theoretical bound predicted by the analysis?
• RQ4 Can the $\mathcal{H}$-matrix format be extended to $\mathcal{H}^2$-matrices or $\mathcal{H}$-LU factorizations for the fractional Laplacian?
• RQ5 What is the empirical performance of $\mathcal{H}$-matrix approximation on different domains and fractional orders?

Key findings

• The inverse of the stiffness matrix for the fractional Laplacian is $\mathcal{H}$-matrix approximable at an exponential rate in the block rank, with theoretical error bound $\|\mathbf{A}^{-1} - \mathbf{B}_{\mathcal{H}}^{r}\|_2 \lesssim e^{-br^{1/4}}$.
• Numerical experiments show a faster empirical convergence rate than predicted, with error decaying as $\|\mathbf{A}^{-1} - \mathbf{B}_{\mathcal{H}}^{r}\|_2 \sim e^{-10r^{1/3}}$.
• The method is robust across different domains, including square and L-shaped domains, and for varying fractional orders $s = 0.25, 0.5, 0.75$.
• The $\mathcal{H}$-matrix approximation of the inverse enables log-linear complexity arithmetic, supporting efficient solution of fractional PDEs via preconditioning or direct solvers.
• The analysis generalizes to higher-order finite elements (piecewise polynomials of fixed degree $p$) and extends to $\mathcal{H}^2$-matrices and $\mathcal{H}$-LU factorizations.
• The numerical results confirm that the $\mathcal{H}$-matrix format is computationally viable for the inverse of the fractional Laplacian, with error decreasing rapidly even at moderate block ranks.