LR-WaveHoltz: A Low-Rank Helmholtz Solver

Basic Information

• Paper ID: 2510.09352
• Title: LR-WaveHoltz: A Low-Rank Helmholtz Solver
• Authors: Andreas Granath (Umeå University), Daniel Appelö (Virginia Tech), Siyang Wang (Umeå University)
• Classification: math.NA, cs.NA (Numerical Analysis)
• Publication Date: October 13, 2025
• Paper Link: https://arxiv.org/abs/2510.09352

Abstract

This paper proposes a low-rank method for solving the Helmholtz equation. The method is based on the WaveHoltz approach, computing Helmholtz solutions by imposing time-domain filters on solutions of the associated wave equation. The wave equation is discretized using high-order summation-by-parts (SBP) finite differences. Singular Value Decomposition (SVD) is employed in two dimensions, while tensor trains are used in three dimensions to compress numerical solutions. To control rank growth, step-wise truncation is applied during time stepping, and low-rank Anderson acceleration is applied to the WaveHoltz fixed-point iteration. Extensive numerical experiments validate the convergence and effectiveness of the iterative scheme for constant and piecewise-constant wave speeds in two- and three-dimensional free space and half-space problems.

Research Background and Motivation

Problem Background

The Helmholtz equation is the fundamental equation for frequency-domain acoustic modeling, with widespread applications in underwater acoustics, including sonar detection, seismic exploration, and long-distance communication. The equation takes the form:

∆u(x) + ω²u(x) = f(x) in Ω ⊂ Rᵈ

where u(x) represents acoustic pressure, f(x) is the source term, and ω is the frequency.

Core Challenges

1. Indefiniteness Problem: Discretization of high-frequency Helmholtz equations leads to indefinite linear systems, causing conjugate gradient methods to fail and GMRES to converge slowly
2. Dispersion Analysis Limitations: To achieve accuracy ε, the number of points per wavelength (PPW) ~ (ω/ε)^(1/2p) is required, with computational cost growing rapidly with frequency
3. Computational Complexity: Traditional methods face enormous storage and computational requirements in high-dimensional problems

Research Motivation

While existing methods such as ray tracing, normal modes, and parabolic equation approaches are effective, direct solution of the Helmholtz equation remains challenging. The authors propose exploiting low-rank structures to reduce computational complexity, particularly leveraging the potential low-rank characteristics of single-point source problems in underwater acoustics.

Core Contributions

1. Proposed LR-WaveHoltz Framework: Combines the WaveHoltz method with low-rank techniques, supporting Helmholtz equation solving in two and three dimensions
2. Multi-dimensional Low-Rank Representation: SVD for two dimensions, tensor train format for three dimensions for solution compression
3. Rank Control Strategy: Develops step-wise truncation methods to control rank growth during time evolution, with scheduling strategies
4. Acceleration Algorithm: Implements low-rank Anderson acceleration (LRAA) to enhance WaveHoltz iteration convergence
5. Multi-block SBP-SAT Framework: Constructs stable high-order multi-block summation-by-parts finite difference schemes
6. Comprehensive Numerical Validation: Verifies method effectiveness on free space, half-space, and layered media problems

Methodology Details

Problem Formulation

Solving the Helmholtz equation with boundary conditions:

∇·(c²∇u(x)) + ω²u(x) - iωκ(x)u = f(x) in Ω
iaωu(x) + bc²∇u(x)·n = 0 on ∂Ω

where c is wave speed, κ is damping, and a, b are boundary condition parameters.

WaveHoltz Method Foundation

The WaveHoltz method transforms the Helmholtz problem into an associated wave equation:

wₜₜ(x,t) + κ(x)wₜ(x,t) = ∇·(c²∇w(x,t)) - f(x)cos(ωt)

Frequency-domain filtering is performed via the WaveHoltz operator Π:

Π[v₀(x), v₁(x)] = (2/T)∫₀ᵀ (cos(ωt) - 1/4)[w(x,t), wₜ(x,t)]dt

Low-Rank Representation Strategy

Two-Dimensional SVD Representation

For structured grids, the numerical solution is represented as a matrix W ∈ Rⁿˣⁿ, utilizing SVD decomposition:

W = USVᵀ

Storage requirements are reduced from n² to 2nr + r (when r << n).

Three-Dimensional Tensor Train Representation

For three-dimensional cases, the tensor train format is employed:

Ã(i₁,i₂,i₃) = Σ G₁(α₀,i₁,α₁)G₂(α₁,i₂,α₂)G₃(α₂,i₃,α₃)

Storage complexity is O(3nr²), substantially smaller than full-rank O(n³).

Step-Wise Truncation Method

To control rank growth during time evolution, an explicit step-wise truncation strategy is employed:

1. Standard Time Stepping: Leapfrog scheme for temporal discretization
2. Truncation Operation: Apply truncation operator Tₑ after each step to maintain specified accuracy
3. Summation Truncation: Implement low-rank matrix summation algorithm T^sum_ε

Anderson Acceleration

Extends Anderson acceleration to low-rank form, solving the optimization problem:

γ^(k) = argmin_u Σₗ₌₁ᵖ ||Dₗᵏu - Fₗᵏ||²

Efficient computation is achieved by solving the simplified normal equations Aγ^(k) = b.

Experimental Setup

Test Problems

1. Free Space Problems: Fully open boundary conditions
2. Half-Space Problems: Water surface reflection boundary conditions
3. Layered Media: Discontinuous wave speed distributions
4. Layered Half-Space: Combined reflection and layering effects

Numerical Parameters

• Spatial Discretization: 4th-order SBP finite difference operators
• Time Step: Δt = 0.15h
• Penalty Parameter: τ = 15
• Point Source Approximation: Gaussian function f(x,y) = -(1/δ²)exp(-r²/δ²), δ = 1/(2ω)

Evaluation Metrics

• Convergence: Frobenius norm residual ||W^(k+1) - W^k||
• Compression Ratio: Rank comparison with full-rank solutions
• Computational Efficiency: Runtime comparisons
• Accuracy: Error relative to analytical or full-rank solutions

Experimental Results

Compression Performance

• Two-Dimensional Case: Significant compression achieved in regions far from the source point, with runtime improvements of up to one order of magnitude
• Three-Dimensional Case: Compression effects are more pronounced, achieving nearly two orders of magnitude acceleration at PPW=40

Convergence Analysis

1. Free Space: Both LRWH and LRAA methods perform well, with limited acceleration effects
2. Challenging Problems: For half-space problems with reflection, LRAA demonstrates significant acceleration, saving 50-80 iterations
3. Rank Growth: Rank growth is approximately monotonic across all tests, with final rank influenced by distance from the source point

Specific Numerical Results

• Free Space: Convergence tolerance ε* = 10⁻³, integral residual reaches 3.33×10⁻⁶
• Layered Half-Space: LRAA(16) saves approximately 80 iterations compared to LRWH
• Three-Dimensional Problems: Maximum tensor train rank remains within reasonable bounds at PPW=10

Rank Behavior Characteristics

1. Spatial Distribution: Higher rank near the source region, significantly reduced rank in far-field regions
2. Temporal Evolution: Rank growth is essentially monotonic with occasional minor fluctuations
3. Truncation Effects: Numerical rank consistently remains within theoretical truncation bounds

Related Work

Traditional Helmholtz Solving Methods

• Ray Tracing Methods: Suitable for high-frequency approximations
• Normal Mode Methods: Based on modal decomposition
• Parabolic Equation Methods: Applicable to far-field propagation

Low-Rank Method Development

• Dynamic Low-Rank Approximation: Projection methods maintaining fixed rank
• Rank-Adaptive Methods: Rank adjustment based on accuracy control
• Tensor Decomposition: Tensor train representation for high-dimensional problems

WaveHoltz Method Evolution

• Basic Framework: Time-domain filtering for frequency-domain problems
• Acceleration Techniques: Krylov subspace methods
• Extended Applications: Elastic wave and electromagnetic wave problems

Conclusions and Discussion

Main Conclusions

1. Method Feasibility: LR-WaveHoltz successfully combines low-rank techniques with the WaveHoltz method
2. Computational Advantages: Significant computational acceleration achieved in three-dimensional problems; limited benefits in two dimensions
3. Convergence Stability: Method demonstrates stable performance across various boundary conditions and media configurations
4. Rank Control Effectiveness: Step-wise truncation and scheduling strategies successfully control rank growth

Limitations

1. Geometric Constraints: Method is applicable to structured multi-block grids with limited geometric complexity
2. Near-Source Region: Higher rank near the source point limits compression effectiveness
3. Two-Dimensional Benefits: Low-rank advantages in two dimensions are less pronounced than in three dimensions
4. Problem-Specific: Primarily targeted at single-point source underwater acoustic problems

Future Directions

1. Hybrid Methods: Combine traditional solvers for near-field with low-rank methods for far-field
2. Geometric Extensions: Extend to more complex geometries and unstructured grids
3. Multi-Source Problems: Handle multiple point sources and distributed source problems
4. Higher-Order Boundary Conditions: Integrate more accurate non-reflecting boundary conditions

In-Depth Evaluation

Strengths

1. Technical Innovation: First systematic application of low-rank techniques to the WaveHoltz method
2. Theoretical Completeness: Provides comprehensive mathematical framework and stability analysis
3. Implementation Details: Detailed algorithm descriptions and implementation techniques
4. Experimental Comprehensiveness: Comprehensive testing across multiple problem types and dimensions
5. Practical Application Value: Addresses the important application domain of underwater acoustics

Weaknesses

1. Scope of Applicability: Limited to structured grids and specific geometries
2. Theoretical Analysis: Lacks rigorous theoretical guarantees for the existence of low-rank structures
3. Parameter Selection: Selection strategies for truncation tolerance and scheduling parameters require further investigation
4. Benchmark Comparisons: Lacks detailed comparisons with other modern Helmholtz solvers

Impact and Significance

1. Academic Contribution: Provides new insights for applying low-rank numerical methods to wave problems
2. Practical Value: Provides a viable tool for large-scale underwater acoustic simulation
3. Methodological Significance: Demonstrates the potential of combining time-domain to frequency-domain transformation with low-rank techniques
4. Reproducibility: Detailed algorithm descriptions facilitate implementation and verification

Applicable Scenarios

1. Underwater Acoustics: Ocean acoustic modeling and sonar system design
2. Seismic Exploration: Large-scale seismic wave propagation simulation
3. Architectural Acoustics: Indoor sound field analysis and noise control
4. Medical Imaging: Ultrasound imaging and therapeutic applications

References

The paper cites 38 important references spanning multiple fields including numerical analysis, low-rank methods, and wave equation solving, providing a solid theoretical foundation for the research.

Overall Assessment: This is a high-quality numerical analysis paper that successfully introduces low-rank techniques to Helmholtz equation solving, with significant contributions in both theoretical methodology and numerical experiments. Despite geometric limitations, it provides a valuable new tool for large-scale acoustic simulation.