PaperlayerView exact revision on arXiv

Source-linked AI summary

An O(N) Direct Solver for Integral Equations on the Plane

Eduardo Corona, Per-Gunnar Martinsson, Denis Zorin

arXiv:1303.5466v3math.NAcs.NA

TL;DR

Iterative integral-equation solvers can be challenged by convergence speed, while prior direct approaches for planar volume problems have higher asymptotic cost. This paper uses hierarchical structured-matrix compression to build a direct solver with O(N) complexity for two-dimensional systems with non-oscillatory or low-frequency kernels, reporting practical performance at large sizes and high accuracy.

  • Problem

    Iterative solvers can become costly when convergence is slow, and prior direct-solver adaptations for planar volume equations require O(N 3/2) work to build an inverse approximation.

    Source: Optimal complexity iterative solvers… · p. 2

  • Method

    The solver combines recursive skeletonization with two levels of HSS structure and directly constructs a compressed representation of the inverse.

    Source: §1.2 · p. 3, §6 · p. 28

  • Results

    For translation-invariant kernels such as Laplace, N = 1 × 10^7 at ε = 10−10 takes half an hour and 5 GB to invert, with 20-second solves.

    Source: §1.2 · p. 3

  • Takeaways & Limitations

    The solver achieves O(N) complexity for two-dimensional integral-equation systems with non-oscillatory or low-frequency kernels, with high performance reported at large sizes and high target accuracies.

    Source: §1.2 · p. 3, §6 · p. 28

    As frequency increases, the stated rank-growth assumptions eventually fail, so a different approach is needed for moderate- and high-frequency regimes.

    Source: §5.3 · p. 24

Abstract

from arXiv · show

An efficient direct solver for volume integral equations with O(N) complexity for a broad range of problems is presented. The solver relies on hierarchical compression of the discretized integral operator, and exploits that off-diagonal blocks of certain dense matrices have numerically low rank. Technically, the solver is inspired by previously developed direct solvers for integral equations based on "recursive skeletonization" and "Hierarchically Semi-Separable" (HSS) matrices, but it improves on the asymptotic complexity of existing solvers by incorporating an additional level of compression. The resulting solver has optimal O(N) complexity for all stages of the computation, as demonstrated by both theoretical analysis and numerical examples. The computational examples further display good practical performance in terms of both speed and memory usage. In particular, it is demonstrated that even problems involving 10^{7} unknowns can be solved to precision 10^{-10} using a simple Matlab implementation of the algorithm executed on a single core.

1. Introduction.

The paper develops a direct solver for two-dimensional volume integral equations to address the convergence dependence of iterative methods and the higher cost of prior direct approaches. It combines hierarchical compression with additional structure to achieve O(N) complexity for non-oscillatory and low-frequency kernels.

1.1. Previous Work.

1.2. Overview of new results.

  • The direct solver achieves optimal O(N) complexity for two-dimensional integral-equation systems with non-oscillatory kernels or kernels in low-frequency mode.The authors report observed O(N) time and storage for all computation stages, with practical accuracy up to 10−10−10−12.Source: §1.2 · p. 3
  • Around N ∼ 10^5, the solver outperforms the O(N 3/2) algorithm, with a one-order-of-magnitude gain in inverse compression time and storage.For translation-invariant kernels, N = 1 × 10^7 at ε = 10−10 takes half an hour and 5 GB to invert, with 20-second solves.Source: §1.2 · p. 3
  • Two hierarchical levels represent the matrix and selected blocks in HSS form, while direct inverse construction avoids matrix-matrix multiplication.This structured algebra accelerates the recursive-skeletonization approach, whose volume-equation adaptation would otherwise have O(N 3/2) complexity.Source: §1.2 · p. 3
  • For N = 10^7 and ε = 10−10, matrix compression takes 5 minutes and 1 GB of storage, and each matrix apply takes under 10 seconds.The reported compression performance is for the accelerated HSS binary-tree construction.Source: §1.2 · p. 3
  • With minor modifications, the method retains optimal scaling and similar performance gains for two-dimensional volume scattering with low-frequency oscillatory kernels, though costs are considerably higher.Source: §1.2 · p. 3

2. Background.

This section reviews a prior hierarchical direct-solver framework and the notation and algorithms needed to motivate the paper’s approach. It covers HSS compression, matrix-vector application, and inversion, including the limits of the earlier method for two-dimensional volume problems.

2. Background.

  • The earlier solver is linear-complexity for integral equations on planar curves, but compression and inversion rise to O(N^3/2) for 2D volume equations.The review frames that superlinear cost as the limitation the paper addresses.Source: §2 · p. 3

2.1. Notation and Preliminaries.

  • The interpolative decomposition selects matrix columns as a skeleton and represents the remaining columns through interpolation; the discretization is organized in boxes within a binary tree.The domain is assumed to lie in a rectangle with regular-grid samples; other input representations are resampled first.Source: §2.1 · p. 4

2.2. Constructing hierarchically semi-separable matrices.

  • HSS representation recursively reblocks and recompresses off-diagonal interactions, while equivalent densities let the algorithm find block bases from smaller proxy-point matrices instead of costly full block rows.Directly computing the interpolative decompositions of the large block rows costs O(N^2); the proxy construction uses the elliptic-PDE structure.Source: Hierarchical compression of A · p. 5, Using equivalent densities to… · p. 6

2.3. HSS matrix-vector multiplication.

  • HSS matrix-vector multiplication uses an upward pass to aggregate charges onto skeletons and a downward pass to propagate outside fields and add sibling and self-interactions.The authors describe the procedure as similar to the fast multipole method, but simpler because near-field interactions are not treated separately.Source: §2.3 · p. 7, Upward pass · p. 7

2.4. Computing the HSS form of A−1.

  • For a nonsingular HSS matrix, inversion proceeds recursively from fine to coarse to compute a compressed representation of the inverse; the reviewed algorithm costs O(N^3/2) for 2D volume equations.The paper uses this inversion algorithm as the starting point for deriving an O(N) method.Source: §2.4 · pp. 7–8

3. O(N) Inverse Compression Algorithm.

For 2D volume integral equations, the existing algorithm typically has O(N^3/2) complexity. The section develops a faster inversion by exposing compressible inverse blocks and avoiding non-low-rank matrix products.

  • O(N^3/2) is the typical complexity of the existing algorithm for 2D volume integral equations.Source: §3 · p. 8
  • The reformulation exposes essential blocks of A^-1 as compressible operators on 1D-like point sets, avoiding non-low-rank matrix-matrix products.Source: §3 · p. 8

3.1. Efficient skeleton construction in 2D.

The section replaces costly repeated interpolative decompositions in 2D with skeleton selection restricted to thin boundary layers. For the considered kernels, the required layer depth remains small even at high accuracy.

3.1. Efficient skeleton construction in 2D.

  • For 2D, applying interpolative decomposition at every hierarchy level yields at least O(N^3/2) overall complexity.The corresponding 1D inversion estimate is O(N).Source: §3.1 · p. 9, Structure of skeleton sets in 2D · p. 10
  • For considered kernels, skeletons can be selected from thin box-boundary layers; ε ∼10^-5 uses one layer and ε ∼10^-10 uses two.The layer depth increases with target accuracy but does not grow too large in the reported experiments.Source: Structure of skeleton sets in 2D · p. 10
  • The algorithm restricts candidate skeleton points to m boundary layers and uses all parent-box boundary-layer points instead of decomposing the merged child skeletons.The section motivates the boundary-layer observation through potential-theory representation results and numerical experiments.Source: Structure of skeleton sets in 2D · p. 10
  • Fixed-width skeleton sets yield compressible matrices analogous to boundary-integral operators, enabling linear-complexity matrix algebra.These matrices act like boundary-to-boundary operators on each box.Source: Structure of skeleton sets in 2D · p. 10

3.2. Overview of the modified inversion algorithm.

The modified HSS inversion reformulates per-box computations so matrix blocks can be compressed as low-rank or one-dimensional HSS operators. It avoids constructing the blocks to be compressed in dense form.

3.2. Overview of the modified inversion algorithm.

  • Dense storage and inversion of the per-box E_i and F_i blocks, or dense construction of interpolation matrices, cannot yield linear complexity.The reformulation partitions these computations into block operations that allow low-rank or one-dimensional HSS compression.Source: §3.2 · p. 11
  • Blocks are constructed directly in compressed form, using low-rank and one-dimensional HSS representations rather than explicitly forming the dense blocks.The reformulation organizes E_i and F_i^-1 computations as block operations.Source: §3.2 · p. 11, Building Ei · p. 12
  • For non-leaf boxes, F_i acts on the merged child skeletons and combines child-local operators with sibling interactions.The section expects its rank structure to resemble that of the interaction operator restricted to the merged skeleton set.Source: Building and inverting Fi · p. 11
  • The reformulated algorithm’s output includes E_i, and its remaining block constructions use HSS and low-rank operations.The section labels the compression used for each output block as low rank or one-dimensional HSS.Source: Building Ei · p. 11

3.3. Compressed two-dimensional HSS inversion.

The section organizes 2D HSS inversion around fast arithmetic and compressed per-box operators. It builds interpolation operators in low-rank form, compresses F_i^-1 blockwise, and constructs E_i in dense-block HSS form.

3.3.1. Fast Arithmetic.

  • The per-box subroutines target costs and storage proportional to O(n^1/2), or O(n^1/2 log^q(n)), and use randomized interpolative decomposition to obtain the needed compression complexity.The section also arranges structured matrix products between dense-block HSS matrices and low-rank matrices, costing O(qk).Source: Operator Notation · p. 12, §3.3.1 · p. 13

3.3. Compressed two-dimensional HSS inversion.

  • The section frames per-box computation as fast when operations and storage scale with the box block size O(n^1/2), up to a logarithmic factor.Its three main subroutines build low-rank interpolation operators, compressed inverse blocks, and E_i in dense-block HSS form.Source: Operator Notation · p. 12

3.3.2. Interpolation Operators in low-rank form.

  • Interpolation operators T_i are built in low-rank form using proxy charges and randomized ID; the proxy-to-residual interaction is low rank when most residual points are well separated from the proxy.The skeleton-to-proxy operator has dense-block HSS structure, while the interfacial points closest to the proxy are an exception to the separation argument.Source: Operator Notation · p. 12, §3.3.2 · pp. 13–14

3.3.3. Compressed F −1.

  • For F_i^-1, the skeleton and residual sets are ordered along the box boundary and sibling interface; the diagonal blocks use dense-block HSS and the off-diagonal blocks use low rank.The off-diagonal interactions are compressible because the corresponding sets are close only at a few points.Source: §3.3.3 · p. 14

3.3.4. Compressed E.

  • E_i retains dense-block HSS structure because it is a low-rank perturbation of the Schur-complement operator S_i^rs.The construction factors the last three low-rank terms, recompresses them, and converts the result to dense-block HSS form.Source: Output · p. 15

3.4. Inverse matrix-vector multiplication.

The algorithm applies the compressed inverse to a right-hand side using an upward pass followed by a downward pass. Fast matrix-vector products replace applications of the compressed operators.

3.4. Inverse matrix-vector multiplication.

  • Fast matrix-vector multiplication replaces applications of the compressed operators and blocks during the inverse apply.The fast products use L_i, R_i, E_i, and blocks of F_i^−1.Source: §3.4 · pp. 15–16, §11 · p. 16
  • The output is σ = A−1f, computed using the binary tree, the HSS-compressed inverse blocks, and the input field values.The inputs include the tree with skeleton indices and associated R_i and L_i, per-box F_i and E_i blocks, and f.Source: §3.4 · p. 16
  • The algorithm computes an upward pass in fine-to-coarse order, then a downward pass to produce A−1f on each leaf box.The downward pass computes φdn_i; the supplied description identifies the result on leaf B_i.Source: §3.4 · p. 16

4. Complexity Estimates.

The analysis estimates work and storage level by level on a binary tree, then sums the per-box bounds. Under stated matrix-structure assumptions, compressed-block algorithms attain linear complexity.

4.1. Complexity of algorithms on binary tries.

  • For per-box work growing slower than linearly, total work is linear; linear growth adds a log N factor, while superlinear growth makes top boxes dominate.For translation-invariant kernels, one set of matrices is stored per level, so the top-level single-box bound reflects overall complexity.Source: §4.1 · p. 17

4.2. Assumptions on matrix structure.

  • The analysis assumes skeleton sizes grow as O(n_ℓ^1/2), interactions can be localized, skeletons lie near box boundaries, and compressed blocks have logarithmic rank growth.The compressed-block algorithms additionally require the fourth assumption; the first three apply to dense- and compressed-block algorithms.Source: §1 · p. 17, §4.2 · p. 17

4.3. Estimates.

  • Dense-block algorithms require O(N^3/2) work to build the tree and compress the inverse; compressed-block algorithms use O(N) work and storage for NTI and TI cases.For TI compressed-block algorithms, inverse compression work and storage are sublinear: O(N^1/2 log^2 N) and O(N^1/2 log^2 N), respectively.Source: Storage · p. 18
  • Optimal complexity can extend to outer HSS ranks growing as O(n_ℓ^p) for p < 1, if the most expensive operations also remain sublinear.The stated sufficient condition also requires low-rank blocks to have rank q_ℓ ∼ n_ℓ^min{1/3,(p−1)/2p}.Source: Storage · p. 18

5. Numerical Results.

The numerical experiments assess inverse compression for translation-invariant and non-translation-invariant kernels using a Matlab implementation on a single processor.

5. Numerical Results.

  • The experiments test the inverse compression algorithm on both translation-invariant and non-translation-invariant kernels.The implementation ran in Matlab on an Intel Xeon X5650 2.67 GHz processor without parallelization.Source: §5 · p. 19
  • The general formulation covers integral equations with smooth coefficient functions and a kernel related to a free-space Green’s function.Source: General formulation · p. 19
  • The test problems use a Nyström discretization on a regular grid in Ω = [−1, 1]^2 to form the dense system Aσ = f.The stated setup uses points with spacing h and a kernel related to a free-space Green’s function.Source: General formulation · p. 19

5.1. High accuracy: performance and scaling.

The experiments test high-accuracy inverse compression and application for TI and NTI kernels. Compressed-block methods scale better in the reported large problems, while inverse application remains fast.

5.1. High accuracy: performance and scaling.

  • Most tests use target accuracy ε = 10^-10 and compare tree construction, inverse construction and compression, and inverse matrix-vector multiplication.The NTI setup uses the 2D Laplace free-space Green’s function; the tests also cover TI kernels and report analogous behavior for other non-oscillatory Green’s functions.Source: §5.1 · p. 19

5.1.1. Non-translation-invariant kernel.

  • At N = 10^7, NTI accelerated inversion gains one order of magnitude in speed and uses 2.5× less memory, or 150 GB.Dense-block inversion and tree-build times scale as O(N^3/2), versus O(N) for the accelerated algorithm; memory scales as O(N log N) versus O(N).Source: §5.1.1 · pp. 19–20, Figure 5.1 · p. 20

5.1.2. Translation-invariant kernel.

  • At N = 10^7, TI compressed-block inverse compression is about 11× faster and uses 30× less memory than the NTI case.It takes about 29 minutes and uses 5 GB; building the binary tree takes 4.5 minutes and 1.3 GB.Source: §5.1.2 · p. 21
  • At N = 10^7, the TI algorithm is about 20× faster than dense-block HSS and uses 8× less memory, or 5 GB.Both methods take about 10 seconds near N = 50,000; the TI algorithm’s memory scaling is sublinear.Source: §5.1.2 · p. 21, Table 5.2 · p. 21, Figure 5.2 · p. 22

5.1.3. Inverse matrix-vector multiplication.

  • Inverse application scales linearly and stays well under a minute for N near 10^7; both algorithms meet the target accuracy on the tested random right-hand sides.The compressed-block apply becomes marginally faster around N ∼ 3 × 10^6.Source: §5.1.3 · p. 21, Table 5.3 · p. 22

5.2. The effect of varying accuracy.

The section compares how lowering target accuracy affects the dense-block and compressed-block algorithms. The measured changes broadly match the predicted accuracy-dependent costs, with lower accuracy reducing compression time and storage.

5.2. The effect of varying accuracy.

  • The analysis predicts accuracy-dependent factors of 16× for compressed-block tree build and inverse compression, and 8× for inverse apply and storage.The estimates assume logarithmic growth in skeleton and dense-block sizes with 1/ε.Source: Compressed-block algorithm · p. 23
  • As N grows, HSS-D compression-time ratios approach 8 and memory-use ratios approach 4 between ε = 10^-10 and ε = 10^-5, matching the estimates.The compressed-block results vary more, but the estimates are generally conservative.Source: Compressed-block algorithm · p. 23
  • Inverse-apply time ratios between the two accuracies are generally 2–3 across the tested cases, consistent with conservative estimates.Source: Compressed-block algorithm · p. 23
  • At ε = 10^-5, both algorithms compress the inverse 8 times faster and store it using 4 times less memory than at ε = 10^-10.The inverse can also be applied two or three times faster at the lower accuracy.Source: Compressed-block algorithm · p. 23

5.3. Low Frequency Oscillatory Kernels.

The solver handles low-frequency Helmholtz kernels with modified skeleton sets and shows linear-scaling inverse compression, while increasing wavenumber raises costs and eventually violates the rank-growth assumptions. Translation invariance improves performance, and inverse applies remain fast, though accuracy can decline at higher tested wavenumbers.

5.3. Low Frequency Oscillatory Kernels.

  • The method handles low-frequency oscillatory kernels with minor modifications, but conditioning can cost several digits of accuracy; one iterative-refinement round can recover them.The required skeletons also gain layers and retain interior points as wavenumber increases.Source: §5.3 · pp. 23–24
  • At κ = 16, 32, accurate compression requires at least three skeleton layers plus extra interior points per wavelength; higher wavenumbers ultimately break the rank-growth assumptions.At κ = 4, 8, two boundary layers suffice, similar to the Laplace case.Source: §5.3.1 · p. 24, §5.3 · p. 24, Figure 5.3 · p. 24

5.3.1. Inverse compression results: TI and NTI cases.

  • HSS-C inverse compression for the non-translation-invariant kernel approaches O(N) scaling; at N ≈ 10^6 it uses 3–4 times less memory than HSS-D.At that size, compression takes about 1.3 hours for κ = 4, 8 and 7–10 hours for κ = 16, 32, with storage reaching approximately 60 GB.Source: §5.3.1 · p. 24, Figure 5.4 · p. 25
  • The compressed-block algorithm beats HSS-D at break-even sizes that increase with κ: about 10^5 for κ = 4, 8 and near 10^6 for κ = 16, 32 in the NTI case.For TI, the break-even point ranges from below 10^5 at κ = 4 to around 3 × 10^5 at κ = 32.Source: §5.3.1 · pp. 24–25
  • By N ≈ 10^6, translation-invariant HSS-C compression is 5–6× faster and uses 15× less memory than its dense-block counterpart.Translation invariance gives gains at all tested wavenumbers; the dense-block algorithm remains O(N^3/2), while HSS-C exhibits sublinear scaling.Source: §5.3.1 · pp. 24–25, Figure 5.5 · p. 25
  • For the translation-invariant case at N ≈ 10^6, HSS-C takes about 15 minutes and 1.5 GB for κ = 4, 8, rising to 1.2–2.6 hours and 4.6–5 GB for κ = 16, 32.The higher-wavenumber cases take substantially longer than the Laplace TI case, while their storage differences are similar to the non-translation-invariant comparison.Source: Figure 5.5 · p. 25, §5.3.1 · p. 26

5.3.2. Inverse matrix-vector multiplication.

  • The HSS-C inverse apply remains fast and optimally scaled, although κ = 16, 32 typically reduces accuracy from the 10^-10 target to about 10^-8.Increasing compression accuracy or applying one round of iterative refinement can address this loss.Source: §5.3.2 · p. 26

5.4. 2D scattering problem: Lippmann-Schwinger equation.

The section tests the direct solver on a 2D acoustic Lippmann–Schwinger scattering problem with a spatially varying wave-number deviation. The discretized operator and inverse show the expected fourth-order convergence, and a translation-invariant Helmholtz inverse serves as a size-independent-iteration preconditioner.

Background on scattering

  • The scattering model writes the outgoing field as an integral-equation solve driven by the incoming point-source field and the deviation b(x), which is zero outside Ω.The test uses a nonnegative b(x) and assumes the wave speed is constant outside Ω.Source: Background on scattering · pp. 26–27

Test problem

  • The test uses an O(h^4) corrected trapezoidal rule; because its correction is diagonal and compression targets off-diagonal blocks, it does not affect solver performance.Higher-order quadratures up to O(h^10) are also applicable.Source: Background on scattering · p. 27

Numerical results

  • The approximate matrix and its inverse converge at approximately O(h^4) = O(N^-2) under grid refinement until errors approach the target accuracy.Source: Background on scattering · p. 27
  • For fixed κ, BiCGstab iteration counts for the preconditioned system are moderate and independent of problem size N.The preconditioner is the inverse of the translation-invariant operator from the low-frequency Helmholtz section.Source: Background on scattering · p. 27

6. Conclusions and Future Work.

The paper concludes that its compressed direct solver achieves O(N) complexity and practical efficiency for non-oscillatory and moderately oscillatory volume integral equations in the plane. It identifies 3D surface problems, parallelization, and moderate- to high-frequency solvers as future directions.

6. Conclusions and Future Work.

  • The solver achieves O(N) complexity and high practical efficiency for non-oscillatory or moderately oscillatory planar volume integral equations, including high-accuracy large problems.Translation invariance can yield sublinear inverse-compression time and storage because only one set of structured matrices is needed per tree level.Source: §6 · p. 28
  • A second HSS-compression level accelerates the recursive-skeletonization solver from O(N^3/2) to O(N), with better compression than the dense-block method and faster runtime beyond N ≈ 10^5.The solver directly builds a compressed inverse and avoids matrix-matrix multiplication.Source: §6 · p. 28
  • The tree-build routine can be adapted into an O(N) matrix-compression algorithm competitive with the FMM and usable with iterative methods such as GMRES.Source: §6 · p. 28
  • The authors propose extending the approach to 3D surfaces, parallel implementation, and direct solvers for moderate- and high-frequency oscillatory problems.The 3D extension may face more complex skeleton structures and geometry-dependent implementation challenges.Source: §6 · pp. 28–29
Loading 1303.5466v3…