[Paper Review] A fast solver for spectral element approximation applied to fractional differential equations using hierarchical matrix approximation

Xianjuan Li, Zhiping Mao+3 more|arXiv (Cornell University)|Aug 8, 2018

Fractional Differential Equations Solutions16 references3 citations

TL;DR

This paper presents a fast iterative solver for spectral element discretizations of two-sided fractional diffusion equations using hierarchical matrix (H-matrix) approximation. By approximating the dense stiffness matrix via H-matrix compression and employing a hierarchical LU decomposition, the method reduces solution cost to O(R²Nd log²N + R³Nd log N + Nd²), achieving over 11 orders of magnitude reduction in condition number compared to unpreconditioned systems, while preserving the high-order accuracy of the spectral element method.

ABSTRACT

We develop a fast solver for the spectral element method (SEM) applied to the two-sided fractional diffusion equation on uniform, geometric and graded meshes. By approximating the singular kernel with a degenerate kernel, we construct a hierarchical matrix (H-matrix) to represent the stiffness matrix of the SEM and provide error estimates verified numerically. We can solve efficiently the H-matrix approximation problem using a hierarchical LU decomposition method, which reduces the computational cost to $O(R^2 N_d \log^2N) +O(R^3 N_d \log N)$, where $R$ it is the rank of submatrices of the H-matrix approximation, $N_d$ is the total number of degrees of freedom and $N$ is the number of elements. However, we lose the high accuracy of the SEM. Thus, we solve the corresponding preconditioned system by using the H-matrix approximation problem as a preconditioner, recovering the high order accuracy of the SEM. The condition number of the preconditioned system is independent of the polynomial degree $P$ and grows with the number of elements, but at modest values of the rank $R$ is below order 10 in our experiments, which represents a reduction of more than 11 orders of magnitude from the unpreconditioned system; this reduction is higher in the two-sided fractional derivative compared to one-sided fractional derivative. The corresponding cost is $O(R^2 N_d \log^2 N)+O(R^3 N_d \log N)+O(N_d^2)$. Moreover, by using a structured mesh (uniform or geometric mesh), we can further reduce the computational cost to $O(R^2 N_d\log^2 N) +O(R^3 N_d \log N)+ O(P^2 N\log N)$ for the preconditioned system. We present several numerical tests to illustrate the proposed algorithm using $h$ and $p$ refinements.

Motivation & Objective

• To address the high computational cost and ill-conditioning of dense linear systems arising from spectral element discretization of two-sided fractional diffusion equations.
• To preserve the high-order accuracy of spectral element methods while enabling efficient solution of the resulting dense systems.
• To develop a fast, scalable solver using hierarchical matrix approximation and preconditioning for uniform, geometric, and graded meshes.
• To achieve optimal complexity O(R²Nd log²N + R³Nd log N + Nd²) with condition number independent of polynomial degree P.

Proposed method

• Approximate the singular kernel in the fractional integral operator using a degenerate kernel to enable hierarchical matrix (H-matrix) representation of the stiffness matrix.
• Construct the H-matrix via hierarchical LU decomposition, exploiting the low-rank structure of off-diagonal blocks to reduce storage and computational cost.
• Use the H-matrix approximation as a preconditioner for the original spectral element system to recover high-order accuracy while reducing condition number.
• Leverage structured meshes (uniform or geometric) to further reduce cost to O(P²N log N) for the preconditioned system.
• Utilize block matrix decompositions involving Toeplitz-like structures and scaled integrals to compute H-matrix blocks efficiently.
• Apply a Petrov-Galerkin framework with Jacobi poly-fractonomials as basis functions to ensure high-order convergence on graded meshes.

Experimental results

Research questions

• RQ1 Can hierarchical matrix approximation be effectively applied to the dense stiffness matrices arising from spectral element discretization of two-sided fractional PDEs?
• RQ2 What is the computational complexity and condition number reduction achievable using H-matrix preconditioning in this context?
• RQ3 How does the method preserve the high-order accuracy of spectral element methods while achieving fast solution times?
• RQ4 To what extent can structured meshes further reduce computational cost in the H-matrix preconditioned system?
• RQ5 What is the dependence of the condition number on the polynomial degree P and number of elements N?

Key findings

• The H-matrix approximation reduces the solution cost to O(R²Nd log²N + R³Nd log N + Nd²), with the dominant term being the preconditioned system solve.
• The condition number of the preconditioned system is independent of the polynomial degree P and grows slowly with the number of elements, remaining below order 10 for modest rank R in numerical experiments.
• The preconditioning reduces the condition number by more than 11 orders of magnitude compared to the unpreconditioned system, with greater reduction observed in two-sided fractional derivatives than one-sided.
• For structured meshes (uniform or geometric), the cost of solving the preconditioned system is further reduced to O(P²N log N), enabling efficient p-refinement.
• Numerical tests confirm the method's high accuracy, with spectral element convergence rates preserved when using the H-matrix preconditioner.
• The method maintains high accuracy even for problems with end-point singularities, outperforming linear FEM and PGS methods in L∞-error for smooth right-hand sides.