A Fast Direct Solver for Boundary Integral Equations Using Quadrature By Expansion

TL;DR

The paper addresses the challenge of efficiently solving dense linear systems from boundary-integral equations discretized by Quadrature by Expansion (QBX). It introduces a hierarchical direct solver based on the Hierarchical Semi-Separable (HSS) framework with proxy-based skeletonization to compress far-field interactions, coupled to a QBX-aware discretization. A complete error model is developed, combining multipole expansions and Interpolative Decomposition (ID) errors, and used to automatically choose key parameters such as proxy count and radius. The approach yields 2D linear-time scaling and 3D near $O(N^{3/2})$ scaling with strong empirical validation across geometries and kernels, and is implemented in open-source software for broader use. Overall, the work provides a QBX-compatible, scalable direct solver for boundary-integral problems with a rigorous error framework and practical parameter-tuning guidance.

Abstract

We construct and analyze a hierarchical direct solver for linear systems arising from the discretization of boundary integral equations using the Quadrature by Expansion (QBX) method. Our scheme builds on the existing theory of Hierarchical Semi-Separable (HSS) matrix operators that contain low-rank off-diagonal submatrices. We use proxy-based approximations of the far-field interactions and the Interpolative Decomposition (ID) to construct compressed HSS operators that are used as fast direct solvers for the original system. We describe a number of modifications to the standard HSS framework that enable compatibility with the QBX family of discretization methods. We establish an error model for the direct solver that is based on a multipole expansion of the QBX-mediated proxy interactions and standard estimates for the ID. Based on these theoretical results, we develop an automatic approach for setting scheme parameters based on user-provided error tolerances. The resulting solver seamlessly generalizes across two- and tree-dimensional problems and achieves state-of-the-art asymptotic scaling. We conclude with numerical experiments that support the theoretical expectations for the error and computational cost of the direct solver.

Figures (13)

• Figure 1: Expansion centers $\boldsymbol{c}$ (full) and corresponding target points $\boldsymbol{x}$ (empty) with the corresponding expansion disk (gray).
• Figure 2: Target points $X_i$ (purple) inside the gray disk of radius $r_{\text{tgt}, i}$ centered at $\boldsymbol{c}_{\text{tgt}, i}$, near source point $Y^{(\text{near})}_i$ (cyan) inside the proxy ball of radius $r_{\text{pxy}, i}$, far source points $Y^{(\text{far})}_i$ (red) outside the proxy ball, and proxy points $P_i$ (black) on the circle.
• Figure 3: Multilevel compression and clustering. Black blocks are unchanged and low-rank, while white blocks are full rank and correspond to the diagonals $\boldsymbol{D}^{(\ell)}$ in \ref{['eq:schemes:rec_skeletonization']} (reproduced from Figure 2 in Martinsson2005).
• Figure 4: Geometry in (a) 2D and (b) 3D. Different colors (not unique) denote different clusters and a representative proxy ball (gray) with a center $\boldsymbol{c}_{\text{src}}$ (red dot) is shown on each geometry. The inner ball denotes the cluster bounding ball and the lighter outer circle denotes the proxy ball.
• Figure 5: Sorted cluster size, where the dashed line denotes the mean cluster size.

Theorems & Definitions (12)

• definition \@thmcounterdefinition: Block Separability Gillman2012
• theorem 1: Interpolative Decomposition Martinsson2005
• remark \@thmcounterremark
• remark \@thmcounterremark
• theorem 2
• corollary \@thmcountercorollary
• corollary \@thmcountercorollary
• remark \@thmcounterremark
• lemma \@thmcounterlemma
• proof