pith. sign in

arxiv: 2504.13809 · v5 · submitted 2025-04-18 · 🧮 math.NA · cs.NA

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

Alexandru Fikl , Andreas Kl\"ockner This is my paper

Pith reviewed 2026-05-22 18:47 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords Boundary integral equationsQuadrature by expansionHierarchical semi-separable matricesDirect solverInterpolative decompositionProxy approximationsError analysisNumerical linear algebra
0
0 comments X

The pith

A hierarchical direct solver for QBX-discretized boundary integral equations achieves state-of-the-art asymptotic scaling in two and three dimensions.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper constructs a fast direct solver for linear systems coming from boundary integral equations discretized with Quadrature by Expansion. The method adapts Hierarchical Semi-Separable matrices by approximating far-field interactions through proxy points and compressing them with the Interpolative Decomposition. A multipole-expansion error model allows automatic selection of compression ranks to meet user error tolerances. If the analysis holds, the solver delivers reliable accuracy and favorable scaling without iterative methods or dimension-specific tuning. Readers in computational science would value a reliable direct method for these often ill-conditioned problems.

Core claim

The central claim is that modifications to the standard HSS framework, including proxy-based far-field approximations compatible with QBX and an error model from multipole expansions combined with ID estimates, produce a direct solver that generalizes across two- and three-dimensional problems while attaining state-of-the-art asymptotic scaling.

What carries the argument

Compressed HSS operators constructed via proxy approximations of far-field interactions and the Interpolative Decomposition, with error controlled by a multipole expansion of QBX-mediated proxies.

If this is right

  • The solver generalizes seamlessly to both two- and three-dimensional problems.
  • It achieves state-of-the-art asymptotic scaling in computational cost.
  • An automatic parameter selection procedure based on the error model meets user-specified tolerances.
  • Numerical experiments validate the error and cost predictions of the theory.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The proxy interaction model might be adapted to other local expansion quadrature methods for integral equations.
  • This could lead to hybrid solvers combining the direct approach with iterative refinement for very large systems.
  • The dimension-independent design suggests applicability to higher-dimensional problems if suitable proxies are defined.

Load-bearing premise

The error model based on a multipole expansion of the QBX-mediated proxy interactions and standard estimates for the ID accurately captures the solver error for the constructed HSS operators.

What would settle it

A numerical experiment on a Laplace boundary integral problem where the computed solution error is substantially larger than the target tolerance when parameters are set according to the model.

Figures

Figures reproduced from arXiv: 2504.13809 by Alexandru Fikl, Andreas Kl\"ockner.

Figure 1
Figure 1. Figure 1: Expansion centers c (full) and corresponding target points x (empty) with the corresponding expansion disk (gray). The expansions are constructed in a neighborhood of the boundary (see [16] for details) for every target point x ∈ Σ (see [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Target points Xi (purple) inside the gray disk of radius rtgt,i centered at ctgt,i, near source point Y (near) i (cyan) inside the proxy ball of radius rpxy,i, far source points Y (far) i (red) outside the proxy ball, and proxy points Pi (black) on the circle. in [12]. The basic idea behind this procedure can be performed on a block-by￾block basis as described below. Locally, the geometry we are interested… view at source ↗
Figure 3
Figure 3. Figure 3: Multilevel compression and clustering. Black blocks are unchanged and [PITH_FULL_IMAGE:figures/full_fig_p014_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Geometry in (a) 2D and (b) 3D. Different colors (not unique) denote [PITH_FULL_IMAGE:figures/full_fig_p024_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Sorted cluster size, where the dashed line denotes the mean cluster size. [PITH_FULL_IMAGE:figures/full_fig_p025_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Comparison of a direct solver using the additional weight matrix [PITH_FULL_IMAGE:figures/full_fig_p026_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Comparison of empirical error data (dashed) and the error model from [PITH_FULL_IMAGE:figures/full_fig_p027_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Proxy count estimation based on the model from Corollary 2 as a [PITH_FULL_IMAGE:figures/full_fig_p028_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Forward error for a Laplace (a,c) single-layer potential and a (b,d) [PITH_FULL_IMAGE:figures/full_fig_p029_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Solution error for a Laplace (a,c) single-layer potential and a (b,d) [PITH_FULL_IMAGE:figures/full_fig_p030_10.png] view at source ↗
Figure 11
Figure 11. Figure 11: Timing (in seconds) for a solving a single-layer potential using the [PITH_FULL_IMAGE:figures/full_fig_p031_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: Timing (in seconds) for a solving a double-layer potential using the [PITH_FULL_IMAGE:figures/full_fig_p032_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: Convergence of the direct solver with grid size [PITH_FULL_IMAGE:figures/full_fig_p033_13.png] view at source ↗
read the original 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.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

1 major / 2 minor

Summary. The paper constructs and analyzes a hierarchical direct solver for linear systems from boundary integral equations discretized via Quadrature by Expansion (QBX). It extends the Hierarchical Semi-Separable (HSS) framework with proxy-based far-field approximations and Interpolative Decomposition (ID) compression, introduces QBX-specific modifications, derives an error model from multipole expansions of QBX proxy interactions plus standard ID bounds, develops automatic parameter selection from user error tolerances, and reports numerical experiments supporting the expected error and cost with state-of-the-art asymptotic scaling that generalizes to 2D and 3D problems.

Significance. If the error model is shown to control the total solver error, the work provides a theoretically grounded fast direct solver for QBX-discretized BIEs with automatic tuning and dimension-independent scaling. This would be a useful contribution to the literature on fast direct methods for integral equations, particularly for applications requiring high accuracy in complex geometries.

major comments (1)
  1. [§4] §4: The error model is constructed from a multipole expansion of the QBX-mediated proxy interactions together with classical ID rank-revealing bounds. It does not explicitly propagate the QBX quadrature error or the truncation error of the local expansions through the HSS block compression. Without this propagation, it is unclear whether the user-specified tolerance guarantees the claimed residual bound once the HSS tree depth increases, especially for 3D geometries.
minor comments (2)
  1. [Abstract] Abstract: 'two- and tree-dimensional' is a typographical error and should read 'two- and three-dimensional'.
  2. The description of the automatic parameter selection procedure would benefit from an explicit statement of how the tolerance is mapped to the ID rank and QBX expansion order in the presence of the proxy approximation.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for their careful reading and constructive comments. We agree that the error analysis can be strengthened by making the propagation of QBX quadrature and truncation errors explicit through the HSS hierarchy, and we will revise the manuscript accordingly.

read point-by-point responses

  1. Referee: [§4] §4: The error model is constructed from a multipole expansion of the QBX-mediated proxy interactions together with classical ID rank-revealing bounds. It does not explicitly propagate the QBX quadrature error or the truncation error of the local expansions through the HSS block compression. Without this propagation, it is unclear whether the user-specified tolerance guarantees the claimed residual bound once the HSS tree depth increases, especially for 3D geometries.

    Authors: We agree that an explicit propagation analysis would make the error model more complete. In the revised manuscript we will add a subsection deriving how QBX quadrature errors and local expansion truncation errors accumulate across HSS levels. Using the additive structure of the ID compression bounds together with the multipole expansion of the proxy interactions, we will show that the total residual remains bounded by a constant multiple of the user tolerance, with the constant independent of tree depth. Separate arguments will be given for 2D and 3D, exploiting the fact that the proxy-point far-field approximation and ID rank-revealing properties hold uniformly in dimension. We will also update the numerical experiments to confirm that the observed residuals respect the predicted bound for deeper trees in 3D geometries. revision: yes

Circularity Check

0 steps flagged

No significant circularity; derivation self-contained against external HSS/ID theory

full rationale

The paper constructs an HSS direct solver by adapting proxy-based far-field approximations and the ID to the QBX discretization. It explicitly states that the error model rests on a multipole expansion of QBX-mediated proxy interactions together with standard ID rank-revealing bounds; these are presented as derived results rather than fitted quantities or self-referential definitions. No load-bearing step reduces by construction to an author-defined input, a self-citation chain, or an ansatz imported from prior work by the same authors. The claimed state-of-the-art scaling follows from the established asymptotic properties of HSS matrices once the QBX-specific modifications are shown to preserve the low-rank structure. The derivation therefore remains independent of the present paper's own fitted values or unverified internal assumptions.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The solver rests on existing HSS and ID theory plus a new error model derived from multipole expansions of QBX proxy interactions; no free parameters or invented entities are described in the abstract.

axioms (1)
  • domain assumption Standard estimates for the Interpolative Decomposition apply to the compressed far-field blocks arising from QBX proxy interactions.
    Invoked to establish the error model for the direct solver.

pith-pipeline@v0.9.0 · 5720 in / 1252 out tokens · 40660 ms · 2026-05-22T18:47:07.289634+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

27 extracted references · 27 canonical work pages

  1. [1]

    Journal of Scientific Computing57, 477–501 (2013)

    Ambikasaran, S., Darve, E.: Ano(nlogn) fast direct solver for partial hierarchically semi-separable matrices. Journal of Scientific Computing57, 477–501 (2013). doi: 10.1007/s10915-013-9714-z

    work page doi:10.1007/s10915-013-9714-z 2013

  2. [2]

    Springer Berlin Heidelberg (2012)

    Atkinson, K.E.: Spherical Harmonics and Approximations on the Unit Sphere: An In- troduction. Springer Berlin Heidelberg (2012)

    work page 2012

  3. [3]

    SIAM Journal on Scientific Computing 36, A427–A451 (2014)

    Barnett, A.H.: Evaluation of layer potentials close to the boundary for Laplace and Helmholtz problems on analytic planar domains. SIAM Journal on Scientific Computing 36, A427–A451 (2014). doi: 10.1137/120900253

    work page doi:10.1137/120900253 2014

  4. [4]

    SIAM Journal on Matrix Analysis and Applications28, 603–622 (2006)

    Chandrasekaran, S., Gu, M., Pals, T.: A fast ULV decomposition solver for hierarchically semiseparable representations. SIAM Journal on Matrix Analysis and Applications28, 603–622 (2006). doi: 10.1137/s0895479803436652

    work page doi:10.1137/s0895479803436652 2006

  5. [5]

    SIAM Journal on Scientific Computing26, 1389–1404 (2005)

    Cheng, H., Gimbutas, Z., Martinsson, P.G., Rokhlin, V.: On the compression of low rank matrices. SIAM Journal on Scientific Computing26, 1389–1404 (2005). doi: 10.1137/030602678 QBX Direct Solver 35

    work page doi:10.1137/030602678 2005

  6. [6]

    Frontiers of Mathematics in China7, 217–247 (2012)

    Gillman, A., Young, P.M., Martinsson, P.G.: A direct solver witho(n) complexity for integral equations on one-dimensional domains. Frontiers of Mathematics in China7, 217–247 (2012). doi: 10.1007/s11464-012-0188-3

    work page doi:10.1007/s11464-012-0188-3 2012

  7. [7]

    Acta Numerica18, 243–275 (2009)

    Greengard, L., Gueyffier, D., Martinsson, P.G., Rokhlin, V.: Fast direct solvers for integral equations in complex three-dimensional domains. Acta Numerica18, 243–275 (2009). doi: 10.1017/s0962492906410011

    work page doi:10.1017/s0962492906410011 2009

  8. [8]

    Greengard, V

    Greengard, L., Rokhlin, V.: A fast algorithm for particle simulations. Journal of Com- putational Physics73, 325–348 (1987). doi: 10.1016/0021-9991(87)90140-9

    work page doi:10.1016/0021-9991(87)90140-9 1987

  9. [9]

    Acta Numerica6, 229–269 (1997)

    Greengard, L., Rokhlin, V.: A new version of the fast multipole method for the Laplace equation in three dimensions. Acta Numerica6, 229–269 (1997). doi: 10.1017/s0962492900002725

    work page doi:10.1017/s0962492900002725 1997

  10. [10]

    Part I: Introduction toH-matrices

    Hackbusch, W.: A sparse matrix arithmetic based onH-matrices. Part I: Introduction toH-matrices. Computing62, 89–108 (1999). doi: 10.1007/s006070050015

    work page doi:10.1007/s006070050015 1999

  11. [11]

    Com- puting69, 1–35 (2002)

    Hackbusch, W., B¨ orm, S.: Data-sparse approximation by adaptiveH 2-matrices. Com- puting69, 1–35 (2002). doi: 10.1007/s00607-002-1450-4

    work page doi:10.1007/s00607-002-1450-4 2002

  12. [12]

    SIAM Journal on Scientific Computing34, A2507–a2532 (2012)

    Ho, K.L., Greengard, L.: A fast direct solver for structured linear systems by recursive skeletonization. SIAM Journal on Scientific Computing34, A2507–a2532 (2012). doi: 10.1137/120866683

    work page doi:10.1137/120866683 2012

  13. [13]

    SIAM Journal on Matrix Analysis and Applications35, 725– 748 (2014)

    Ho, K.L., Greengard, L.: A fast semidirect least squares algorithm for hierarchically block separable matrices. SIAM Journal on Matrix Analysis and Applications35, 725– 748 (2014). doi: 10.1137/120902677

    work page doi:10.1137/120902677 2014

  14. [14]

    Ho, K.L., Ying, L.: Hierarchical interpolative factorization for elliptic operators: Integral equations (2013)

    work page 2013

  15. [15]

    Communications on Pure and Applied Mathematics69, 1415–1451 (2016)

    Ho, K.L., Ying, L.: Hierarchical interpolative factorization for elliptic operators: Differ- ential equations. Communications on Pure and Applied Mathematics69, 1415–1451 (2016). doi: 10.1002/cpa.21582

    work page doi:10.1002/cpa.21582 2016

  16. [16]

    Journal of Computational Physics252, 332–349 (2013)

    Kloeckner, A., Barnett, A., Greengard, L., O’Neil, M.: Quadrature by expansion: A new method for the evaluation of layer potentials. Journal of Computational Physics252, 332–349 (2013). doi: 10.1016/j.jcp.2013.06.027

    work page doi:10.1016/j.jcp.2013.06.027 2013

  17. [17]

    URLhttps://github.com/inducer/ pytential

    Kloeckner, A., Fikl, A., et al.: pytential (2025). URLhttps://github.com/inducer/ pytential

    work page 2025

  18. [18]

    Springer Science & Business Media (2013)

    Kress, R.: Linear Integral Equations. Springer Science & Business Media (2013)

    work page 2013

  19. [19]

    Journal of Computational Physics205, 1–23 (2005)

    Martinsson, P.G., Rokhlin, V.: A fast direct solver for boundary integral equations in two dimensions. Journal of Computational Physics205, 1–23 (2005). doi: 10.1016/j.jcp.2004.10.033

    work page doi:10.1016/j.jcp.2004.10.033 2005

  20. [20]

    Princeton University Press (1971)

    Stein, E.M., Weiss, G.L., Mather, G.: Introduction to Fourier Analysis on Euclidean Spaces. Princeton University Press (1971)

    work page 1971

  21. [21]

    SIAM Journal on Scientific Computing36, A267–a288 (2014)

    Vioreanu, B., Rokhlin, V.: Spectra of multiplication operators as a numerical tool. SIAM Journal on Scientific Computing36, A267–a288 (2014). doi: 10.1137/110860082

    work page doi:10.1137/110860082 2014

  22. [22]

    Journal of Computational Physics388, 655–689 (2019)

    Wala, M., Kl¨ ockner, A.: A fast algorithm for quadrature by expansion in three dimensions. Journal of Computational Physics388, 655–689 (2019). doi: 10.1016/j.jcp.2019.03.024

    work page doi:10.1016/j.jcp.2019.03.024 2019

  23. [23]

    Wala, M., Kl¨ ockner, A.: On the approximation of local expansions of Laplace potentials by the fast multipole method (2020)

    work page 2020

  24. [24]

    Womersley, R.S.: Efficient spherical designs with good geometric properties. In: J. Dick, F.Y. Kuo, H. Wo´ zniakowski (eds.) Contemporary Computational Mathematics, pp. 1243–1285. Springer International Publishing (2018). doi: 10.1007/978-3-319-72456-0 57

    work page doi:10.1007/978-3-319-72456-0 2018

  25. [25]

    Applied and Computational Harmonic Analysis49, 316–327 (2020)

    Xing, X., Chow, E.: Error analysis of an accelerated interpolative decomposition for 3D Laplace problems. Applied and Computational Harmonic Analysis49, 316–327 (2020). doi: 10.1016/j.acha.2019.11.003

    work page doi:10.1016/j.acha.2019.11.003 2020

  26. [26]

    SIAM Journal on Matrix Analysis and Applications41, 221–243 (2020)

    Xing, X., Chow, E.: Interpolative decomposition via proxy points for kernel matri- ces. SIAM Journal on Matrix Analysis and Applications41, 221–243 (2020). doi: 10.1137/19m1258700

    work page doi:10.1137/19m1258700 2020

  27. [27]

    SIAM Journal on Matrix Analysis and Applications41, 1059–1085 (2020)

    Ye, X., Xia, J., Ying, L.: Analytical low-rank compression via proxy point selec- tion. SIAM Journal on Matrix Analysis and Applications41, 1059–1085 (2020). doi: 10.1137/19m1247838

    work page doi:10.1137/19m1247838 2020