Technical Summary: LR-WaveHoltz: A Low-Rank Helmholtz Solver

Problem Statement

The paper addresses the computational challenge of solving the Helmholtz equation in two and three dimensions, specifically modeling acoustic pressure $u(x)$ generated by a source $f(x)$ at frequency $\omega$:
$$\Delta u(x) + \omega^2 u(x) = f(x)$$
Direct frequency-domain discretization of this equation is notoriously difficult for high frequencies due to the indefinite nature of the resulting linear system, which hinders the effectiveness of standard iterative solvers (e.g., Conjugate Gradient fails, GMRES offers limited acceleration). Furthermore, classical dispersion analysis dictates that high-order methods require a number of points per wavelength (PPW) scaling as $O((\omega/\epsilon)^{1/2p})$ to maintain accuracy, leading to massive system sizes.

While low-rank methods have shown success for coercive elliptic problems, the Helmholtz solution is not generally expected to be low-rank; indeed, separation rank can grow rapidly with frequency. However, the authors argue that for underwater acoustics problems involving a single point source, specific low-rank structures can be exploited, particularly in regions far from the source where the solution exhibits "parabolic separation."

Methodology

The proposed method, LR-WaveHoltz (LRWH), combines the WaveHoltz time-domain filtering technique with low-rank tensor representations and Anderson acceleration.

1. The WaveHoltz Framework

Instead of solving the Helmholtz equation directly, the method solves a related time-domain wave equation over one period $T = 2\pi/\omega$. The solution is filtered using a specific kernel to extract the $\omega$-mode. This transforms the problem into a fixed-point iteration:
$$\mathbf{v}^{(n+1)} = \Pi[\mathbf{v}^{(n)}]$$
Crucially, the discretized WaveHoltz operator $\Pi$ yields a positive-definite system, making it amenable to Krylov subspace acceleration, unlike the indefinite Helmholtz system.

2. Spatial Discretization (SBP-SAT)

The wave equation is discretized using high-order Summation-By-Parts (SBP) finite difference operators combined with the Simultaneous Approximation Term (SAT) method.

• Multiblock Structure: The domain is partitioned into structured blocks. This allows the point source to be localized in specific blocks while maintaining a structured grid elsewhere.
• Matrix/Tensor Form: In 2D, the solution in each block is represented as a matrix $U(t) \in \mathbb{R}^{n \times n}$. In 3D, it is represented as a 3D array.
• Boundary Conditions: Non-reflecting and Neumann boundary conditions are imposed weakly via SAT terms. Interface conditions between blocks with varying wave speeds are handled via projection or weak enforcement.

3. Low-Rank Representation and Step-Truncation

To manage the computational cost, the solution is compressed at every time step:

• 2D (SVD): Solutions are stored in dyadic SVD form ($U S V^T$).
• 3D (Tensor Trains): Solutions are stored in the Tensor Train (TT) format.
• Step-Truncation: During time stepping (using a leapfrog scheme), the rank naturally grows due to the addition of terms (Laplacian, forcing, previous steps). The authors employ an explicit step-truncation strategy where the sum of low-rank terms is computed and then truncated using SVD (in 2D) or TT-rounding (in 3D) to a specified tolerance $\epsilon$.
• Scheduling: The truncation tolerance is dynamically adjusted based on the residual of the fixed-point iteration ($\epsilon_k \propto \rho_k$) to prevent unnecessary rank inflation while ensuring convergence.

4. Acceleration

The fixed-point iteration is accelerated using Low-Rank Anderson Acceleration (LRAA).

• The method uses a linear combination of previous iterates to minimize the residual.
• The minimization problem is reformulated to operate on the low-rank factors (cores) rather than the full matrices/tensors, reducing the complexity of solving for the acceleration coefficients from $O(mn^2)$ to $O(m^2)$.

Key Contributions

1. Algorithmic Framework: The paper introduces a unified low-rank framework for the WaveHoltz method, integrating SBP-SAT spatial discretization with step-truncation time stepping in both 2D (SVD) and 3D (Tensor Trains).
2. Control of Rank Growth: The authors demonstrate that by using step-truncation combined with residual-based scheduling, the numerical rank remains bounded and grows monotonically toward the final solution rank, avoiding the "intermediate rank inflation" often seen in low-rank time evolution.
3. Anderson Acceleration in Low-Rank: The paper extends Anderson acceleration to the low-rank setting, showing how to efficiently compute the acceleration coefficients using only the low-rank factors.
4. 3D Tensor Train Implementation: The work provides a practical implementation of the WaveHoltz method in 3D using Tensor Trains, a significant step given the memory constraints of full-rank 3D Helmholtz solvers.

Numerical Results

The authors present extensive numerical experiments for free-space, half-space, and layered media problems in 2D and 3D.

• Compression and Efficiency:
  • In 3D, the low-rank method offers substantial speedups (up to two orders of magnitude) compared to full-rank solvers for problems with 40 points per wavelength (PPW). The method allows for the simulation of problems with $2.7 \times 10^7$ degrees of freedom per block, which would be prohibitive for full-rank solvers.
  • In 2D, the low-rank format is efficient for high PPW and regions far from the source, but the gains are less dramatic than in 3D.
• Convergence:
  • The WaveHoltz iteration converges to the Helmholtz solution. The number of iterations scales as $O(\omega^{0.96})$, consistent with theoretical predictions.
  • Acceleration: For open (free-space) problems, Anderson acceleration provides only marginal improvements in iteration count because the WaveHoltz iteration is already efficient. However, for problems with trapped waves (e.g., half-space with reflecting boundaries or layered waveguides), LRAA provides significant acceleration, reducing iteration counts by 50–80 iterations for memory parameters $M=2$ to $M=32$.
• Rank Behavior:
  • The numerical rank is highest near the point source and decreases with distance, aligning with the theoretical expectation of low-rank separability in the far-field.
  • The rank remains bounded by the truncation tolerance and does not exhibit uncontrolled growth during the iteration.

Significance and Claims

The paper positions LR-WaveHoltz as a proof of concept for a fast volume discretization method for underwater acoustics. The authors explicitly state that the method is best suited for domains with low geometric complexity (structured multiblock meshes).

• Modest Claims: The authors acknowledge that for 2D open problems, the WaveHoltz iteration is already highly efficient, and Anderson acceleration does not significantly improve performance. They note that the low-rank benefits in 2D are limited compared to 3D.
• Future Directions: The paper suggests that a more efficient implementation would involve hybridizing traditional solvers near the source (where rank is high) with low-rank solvers in the far-field. This corresponds to separating near- and far-field terms.
• Limitations: The current method relies on structured grids. Extending this to curvilinear grids (necessary for realistic bathymetry) is identified as a major challenge because coordinate transformations destroy the separable structure required for SVD/TT truncation. The authors propose that cross-approximation techniques might be the solution for this future extension.

In summary, the paper demonstrates that low-rank representations, when combined with the WaveHoltz time-domain approach and step-truncation, can effectively solve 3D Helmholtz problems with significant compression, particularly for problems involving trapped waves or large domains, while maintaining the accuracy of high-order finite difference methods.