The method of polarized traces for the 2D Helmholtz equation

1 Introduction

• The proposed hybrid uses legacy direct solvers locally on large subdomains and couples them with accurate transmission, addressing direct methods’ memory-related limits on very large problems.The introduction situates the approach among domain-decomposition and sweeping methods, tracing a lineage from AILU to sweeping preconditioners and coarser layerings.Source: §1 · pp. 1–2

1.3 Complexity scalings

• Online complexity scales as O(N/L), strictly improving on O(N) when 1 << L << N, with empirical exponents depending on frequency scaling.The experiments observe γ = 5/8 for ω ∼ n^1/2 and γ = 7/8 for ω ∼ n; the authors prefer γ = 3/4 based on theoretical arguments, while noting the scalings lack proof for heterogeneous media of limited smoothness.Source: §1 · p. 2, §1.3 · p. 7

1.1 Polarization

• The method defines interface polarization through incomplete Green’s formulas with local Green’s functions, approximating conditions that would be exact with the global Green’s function.The paper’s formulation reduces the discrete Helmholtz system to an interface integral system without Schur complements, while allowing polarization into one-way components.Source: §1 · p. 2, §1.1 · pp. 3–4

1.2 Algorithm

• The algorithm uses polarized traces as interface unknowns and a block-triangular preconditioner whose forward and backward substitutions sweep transmitted and locally reflected waves.Polarization makes the matrix blocks adaptively low-rank compressible, enabling sublinear-time iterations while keeping GMRES iteration counts essentially bounded with frequency and layer count.Source: §1.2 · p. 5

1.3 Complexity scalings

• The authors extend the algorithm to 3D in principle, estimating online runtime O(N^4/3/L^2 + L N^5/6), or potentially sublinear runtime with butterfly application of Green’s integrals.The 3D estimate is conditional: the butterfly application is described as theoretically capable of O(N^2/3) time, up to logarithmic factors.Source: §1.3 · p. 7

1.4 Additional related work

• Related work includes transmission-based domain decomposition, multi-trace interface formulations, approximate-transmission Schwarz methods, and source-transfer approaches.The introduction also notes direct methods based on multifrontal elimination and H-matrices [31] and impedance-to-impedance maps, while raising uncertainty about offline linear complexity.Source: §1 · p. 1, §1.4 · p. 8
• Prior sweeping work includes 3D parallelization, spectrally accurate elastic and Maxwell solvers, and a recursive 3D variant that reduces offline cost to linear complexity.A domain-decomposition sweeping variant coupled with additive Schwarz was also developed; an earlier version introduced polarized traces using Green’s representation and reported favorable scalings.Source: §1.4 · p. 8
• Multigrid-related alternatives include optimized two-level coarse-grid corrections, wave-ray error elimination, relaxed multigrid preconditioning, eikonal-based reformulation, and algebraic multilevel preconditioning.The wave-ray method uses a geometric-optics approximation of the Green’s function; the alternative formulation solves eikonal and advection-diffusion-reaction equations with Krylov and multigrid methods.Source: §1.4 · pp. 8–9
• Complex-shifted Laplace preconditioners are described as simple but suboptimal, with performance dependent on the shift; variants use multilevel shifts, deflation, or optimal grids.The introduction also notes studies of optimal shifts, further extensions, large 3D multigrid implementations, and reviews of Helmholtz iterative methods and high-frequency difficulties.Source: §1.4 · p. 9

1.5 Organization

• The paper proceeds through the discrete integral reformulation, polarization, preconditioner definition, adaptive PLR compression, complexity analysis, and numerical experiments varying L and ω.The polarization trick produces a triangular block structure while doubling the number of unknowns.Source: §1.5 · p. 9

2.1 Discretization

• The 2D Helmholtz equation is discretized on a regular grid with a five-point stencil and absorbing boundary conditions implemented through a PML.The PML modifies the operator to damp outgoing waves; its thickness and absorption strength are tuning parameters.Source: §2.1 · pp. 9–10

2.2 Domain Decomposition

• The extended domain is partitioned into layers, with additional PML padding at internal boundaries and local problems matching the global equations inside each physical layer.The local PML profile is chosen to make the local and global equations coincide in the layer interior.Source: §2.2 · pp. 11–12, Figure 3 · p. 11

2.3 Discrete Green’s Representation Formula

• Local Green’s representation formulas express layer solutions using interface traces and forcing, and their operators need only be stored at interfaces.The incomplete Green’s operators are not traditional Schur complements.Source: §2.3 · pp. 13–14

2.4 Discrete Integral Equation

• Evaluating local Green’s formulas at layer boundaries and enforcing continuity yields an algebraically equivalent discrete integral system for interface unknowns.The system’s solution equals the global solution restricted to the interfaces.Source: §2.4 · pp. 14–15
• The discrete integral system recovers exactly the global discretized Helmholtz solution restricted to subdomain interfaces, without using traditional Schur complements.The formulation is algebraically exact.Source: §2.3 · p. 14, §2.4 · p. 15
• Offline computation factors local problems, computes interface Green’s functions by back substitution, and assembles the integral system for later right-hand sides.The local Green’s-function solves are independent and can be performed in parallel.Source: §2.4 · p. 15
• Directly solving the interface system can be prohibitively expensive, while its observed condition number of O(h^-2) can lead to many iterative steps.The paper therefore focuses on preconditioning the system for propagating-wave regimes.Source: §2.4 · p. 16
• Although the boundary-integral reduction extends to other partitions, the preconditioner’s polarizing conditions require a layered partition.The paper uses a layered decomposition to enforce those conditions.Source: §2.4 · p. 17

3.1 Polarized Wavefields

• Polarized fields are not purely directional: they include locally reflected as well as transmitted waves, with polarization quality inherited from the local Green’s functions.A polarized wave can be extrapolated outside the domain using only its Dirichlet boundary data.Source: §3.1 · p. 17

3.2 Polarized Traces

• Polarized traces split each interface field into up- and down-going components, whose sum gives the original field.The formulation increases unknowns without increasing the equation count, so additional conditions are needed to close the system.Source: §3.2 · p. 19
• Three equivalent closure strategies impose annihilation conditions, combine extrapolator relations with annihilation, or use jump conditions to simplify expressions.The second yields a decoupled block-triangular preconditioner but requires dense diagonal-block inversions; the third avoids that requirement.Source: §3.2 · p. 19

3.3 Annihilation relations

• The annihilation formulation completes the polarized system by imposing annihilation conditions on outgoing traces, alongside continuity conditions.A solution of the completed system gives a solution of the original interface system when the up- and down-going traces are summed.Source: §3.3 · pp. 19–20

3.4 Extrapolation conditions

• The extrapolation formulation encodes polarization with extrapolator relations and permits a row permutation that restores a useful sparsity structure.Its preconditioner uses block Jacobi, but constructing the extrapolator requires inverting small dense blocks, which can be costly for large interfaces.Source: §3.4 · pp. 20–21, §3.5 · p. 22

3.5 Jump condition

• The equivalent jump-condition formulation is benchmarked and avoids the unnecessary dense linear algebra required by extrapolator construction.The jump and extrapolation systems have the same solutions; the numerical claims use the jump system preconditioned by direct inversion of Djump.Source: §3.5 · pp. 22, 24

4.1 Gauss-Seidel Iteration

• Inverting the triangular block factors with identity diagonals performs top-down and bottom-up sweeps that transmit waves using incomplete Green’s formulas.The paper calls this procedure Gauss-Seidel, while noting it could also be called multiplicative Schwarz.Source: §4.1 · p. 24

4.2 GMRES

• GMRES with the Gauss-Seidel sweep as preconditioner has an iteration count roughly constant in the number of subdomains, unlike the standalone iteration’s often increasing count.The standalone count often grows as a fractional power of the number of subdomains.Source: §4.2 · p. 25
• The preconditioned eigenvalues are 1 and 1 ± √µi, and smaller |µi| values cluster the latter around 1, favoring faster GMRES convergence.The clustering intuition tolerates a small number of large-|µi| outliers, whereas the spectral radius is not robust to them in heterogeneous media.Source: §4.2 · p. 26
• The operator (D↑)^-1 L (D↓)^-1 U models successive within-layer reflections and transmissions across the domain, with its eigenvalue magnitudes tied to heterogeneity and PML quality.The authors expect eigenvalues of the preconditioned matrix to cluster around 1 when the medium has no resonant cavities.Source: §4.2 · p. 26

5 Partitioned low-rank matrices

• The interface system’s uncompressed update costs O(n^2L), while fast applications of its blocks could reduce the solve cost toward O(nL).This is the target before accounting for the further conditions needed for sublinear overall online complexity.Source: §5 · p. 27

5.1 Compression

• PLR compression stores low-rank factors in a quadtree of maximally large square blocks with bounded ϵ-rank.The paper implements PLR matrix-vector multiplication through sparse factorization to reduce recursive-operation overhead and use optimized sparse routines.Source: §5.1 · pp. 27, 29, Figure 9 · p. 29
• Adaptive partitioning is suited to interface kernels because nearby blocks have diagonal singularities while remote blocks have higher ranks from wave-vector diversity.The method selects the smallest admissible quadtree partition, with each block’s ϵ-rank bounded by rmax.Source: §5.1 · p. 28

5.2 Compression scalings

• For high-frequency Green’s kernels, low ranks generally require block diameters dA ≤ C/√ω, with smaller admissible blocks when A and B are close.The estimate assumes smooth media and single-valued traveltimes; representative blocks then motivate the compression and multiplication cost analysis.Source: §5.2 · p. 30
• For N ∼ n^2, PLR matrix-vector multiplication has asymptotic complexity N^5/8 when ρ = 1/2, and approximately N^3/4 when ρ = 1.Source: §5.2 · p. 31
• The theoretical scaling assumes known phase functions, but finite-difference Green’s-function errors impair compressibility at large ω.Increasing maximum ranks as √n can improve observed costs, but the reported bracketed scalings are pre-asymptotic and the unbracketed values are the expected asymptotic figures.Source: §5.2 · p. 31

6.1 Computational cost

• Local LU factorizations cost O((N/L)^3/2) and can run independently across L nodes; local solves cost O(N log(N)/L) per solve.The factorization and local-solve stages use nested dissection and back-substitution, respectively.Source: §6.1 · p. 31
• Green’s-function extraction requires O(n) independent solves per layer, with a trade-off between parallelism and communication.The authors identify extraction as the current complexity bottleneck and note that randomized algorithms may reduce the solve count.Source: §6.1 · pp. 31–32
• Numerical Green’s-function errors hinder compression at large ω, motivating milder frequency scaling such as ω ∼ √n.The second-order finite-difference discretization’s faithfulness is given as the reason for considering this scaling.Source: §6.1 · p. 32, §6.2 · p. 32
• With ω ∼ √n, each GMRES iteration costs O(LN^5/8) when the iteration count remains bounded.The block matrix-vector cost is O(N^5/8) under this frequency scaling.Source: §6.1 · p. 32
• When L ∼ N^δ with 0 < δ < 3/8, online execution time excluding communication is sublinear; specifically, δ = 3/16 yields O(N^13/16) runtime plus sublinear communication cost.Source: §6.1 · p. 32

6.2 Communication cost

• The online communication cost is O(αL + βLN^1/2), and communication was negligible relative to floating-point operations in the reported 2D experiments.This estimate assumes sources are distributed and traces are sent to and from the master node.Source: §6.2 · p. 33

7 Numerical Experiments

• The method satisfies the desired iteration stability with respect to subdomain count and frequency, and shows empirical sublinear-complexity behavior in GMRES and online runtime.For rough geophysical models, the third desired property—weak dependence on model roughness—is also observed in cases of interest to geophysicists.Source: §7 · p. 33

7.1 Precomputation

• Independent slab solves and Green's-function columns make precomputation parallelizable; the resulting operators are compressed in sparse PLR form for use in block solves.The experiments used pivoted sparse LU, stored per-slab factors, and randomized SVD compression.Source: §7.1 · p. 34

7.2 Smooth Velocity Model

• For smoothed Marmousi2, GMRES iteration counts are almost independent of frequency and slab count, while runtimes scale sublinearly with volume unknowns.Tables 5 and 6 average timings over 200 randomly generated right-hand sides.Source: §7.2 · p. 34, Table 5 · p. 35, Table 6 · p. 35

7.3 Rough Velocity Model

• On rough geophysical models, performance deteriorates only marginally as frequency and subdomain count increase, and convergence grows slowly.The tested layering direction places model discontinuities orthogonal to interfaces; reversing that direction increases iteration counts, with marginal deterioration as frequency and subdomain count rise.Source: §7.3 · pp. 35–36
• Average solve runtimes scale sublinearly with volume unknowns on Marmousi2, while high-frequency tests show slightly higher, slowly growing iteration counts and slightly sublinear runtimes.With fixed L, later LU costs are lower-bounded by O(N log(N)/L), and GMRES complexity may deteriorate to O(N^3/4) at sufficiently high frequencies.Source: §7.3 · p. 36

7.4 Cavities

• In resonant cavities, increasing frequency spreads preconditioned-system eigenvalues away from 1, including eigenvalues near zero, and the number far from 1 grows with frequency.The paper states that sharp interfaces and resonant cavities increase complexity mainly through more GMRES iterations, and that optimal scaling or accuracy is not expected for resonant cavities.Source: §7.4 · pp. 39–40

Computations for Eq. 32 in 1D

• In 1D, summation by parts supplies the differentiation rule used to compute solution and Green-function traces and derivatives at interfaces.The example uses the three-point stencil approximation of the second derivative.Source: Computations for Eq. 32 in 1D · p. 41

Computations for Eq. 32 in 2D

• For the 2D derivation, symmetric and unsymmetric PML formulations yield the same Green's representation formula, and the discrete formula extends Eq. 11.The derivation uses the symmetric formulation for proofs and applies the resulting relation to the unsymmetric formulation as well.Source: Computations for Eq. 32 in 2D · pp. 42–43

B Triangular and block triangular matrices

• For the 2D block-ordered Helmholtz operator, each block corresponds to a fixed depth, with tridiagonal symmetric diagonal blocks and constant-identity couplings.This structure supports the appendix's block-matrix analysis.Source: B Triangular and block triangular… · p. 45
• The recurrence-defined matrices yield an inverse formula for the proper block Hessenberg operator H.The proposition's stated properness follows from invertibility of the coupling blocks.Source: B Triangular and block triangular… · p. 45
• The generalized Casorati determinant for the recurrence sequences is constant, a property used in deriving the inverse of H.The proof invokes the constant determinant property in the inverse argument.Source: B Triangular and block triangular… · pp. 46–47

Proof of Lemma 1

• The discrete Green's representation formula follows by applying the discrete inner product and summation by parts to the Helmholtz solution and Green's function.The argument also covers source indices outside the domain, where the delta support lies outside the summation region.Source: Proof of Lemma 1 · p. 48
• The identity remains valid when the global Green's function is replaced by local Green's functions.The proof states that the local replacement preserves the result.Source: Proof of Lemma 1 · p. 48

Proof of Lemma 2

• Applying the extrapolator to the fixed-interface expression and computing its components shows equivalence to Eq. 48.The calculation uses the extrapolator definition and Lemma 4, then evaluates the remaining term using the matrix form of Gℓ.Source: Proof of Lemma 2 · p. 48

• The nullity theorem [89] reduces the proof to showing that Gℓ+1(z1, z1) is invertible.Source: Proof of Lemma 3 · p. 49
• After reordering Gℓ+1, the block Gℓ+1(z1, z1) occupies the upper-left corner as a square diagonal block.Source: Proof of Lemma 3 · p. 49
• Its nullity equals that of the complementary block of the inverse, which is full rank; hence Gℓ+1(z1, z1) is invertible.Source: Proof of Lemma 3 · p. 49

• The proof begins from the local Green’s function definition in Eq. 22 and invokes Def. 1 and Eq. 129.Source: Proof of Lemma 4 · p. 49
• It uses the layer-to-layer Green’s function definition and the extrapolator, including Proposition 4 at j = 0, to derive intermediate relations.Source: Proof of Lemma 4 · p. 49
• For physical-domain zk, Ck = −I/h^2; this lets the proof simplify Eq. 169, with a corresponding derivation from Eq. 141.Source: Proof of Lemma 4 · pp. 49–50

Proof of Lemma 5

• The Lemma 5 proof proceeds by induction, using the extrapolator to shift the evaluation index at each stage.Source: Proof of Lemma 5 · p. 50

Proof of Lemma 6

• For Lemma 6, the proof treats j = 0, assumes u solves the system in Def. 8, and follows the computation of Lemma 2.Source: Proof of Lemma 6 · p. 50

Proof of Proposition 1

• The Proposition 1 proof takes sufficiency from Lemma 6 and establishes necessity using the invertibility of the extrapolator; the case i = nℓ is analogous.Source: Proof of Proposition 1 · p. 51

• In 1D, the proof reconstructs local solutions and checks the difference equation at subdomain endpoints using information from neighboring subdomains.Source: Proof of Theorem 1 · pp. 51–52
• The 2D argument requires particular care with the PML and applies the same proof method with the stated interior-index range.Source: Proof of Theorem 1 · p. 52
• The concatenated solution satisfies the finite-difference equation at interfaces and, by construction, the homogeneous Dirichlet boundary conditions; uniqueness makes it the discrete PDE solution.Source: Proof of Theorem 1 · p. 52