A Stabilized Multilevel B-Spline-Based Fast Integral Method for the Solution of the Electric Field Integral Equation

Technical Summary: A Stabilized Multilevel B-Spline-Based Fast Integral Method for the EFIE

Problem Statement
The Electric Field Integral Equation (EFIE) is a standard tool for solving electromagnetic scattering and radiation problems, requiring only surface discretization. However, standard Boundary Element Method (BEM) discretizations result in $O(N^2)$ complexity for matrix assembly, storage, and matrix-vector products (MVP), where $N$ is the number of unknowns. Fast integral methods, such as the Multilevel Fast Multipole Method (MLFMM) and FFT-accelerated approaches like the Adaptive Integral Method (AIM), have been developed to reduce this complexity to $O(N)$ or near-$O(N)$.

A specific limitation exists in FFT-accelerated methods that rely on global Lagrange interpolation. These methods require equidistant interpolation points to leverage FFT structures. However, equidistant Lagrange interpolation suffers from Runge-type instabilities at high polynomial orders, leading to a breakdown in accuracy and preventing robust high-order compression. While B-splines offer a solution to Runge instabilities on equidistant grids, their application in multilevel fast integral methods is non-trivial. Unlike Lagrange polynomials, B-spline coefficients do not coincide with function values at interpolation points, necessitating the inversion of sampling matrices. These matrices can become severely ill-conditioned at high orders, compromising numerical stability. Furthermore, maintaining exact interlevel transfers (required for multilevel efficiency) while ensuring stability has been a significant challenge.

Methodology
The authors propose a stabilized multilevel B-spline-based fast integral method that combines FFT-compatible kernel interpolation with robust high-order interpolation. The core methodology involves three primary components:

1. B-Spline Kernel Interpolation:
   Instead of global Lagrange polynomials, the method employs B-splines defined on equidistant knot vectors. This choice mitigates Runge-type instabilities while maintaining compatibility with equidistant sampling required for FFT acceleration. The Green's function is factorized using tensor-product polynomials. For B-splines, the interpolation weights are not direct samples of the Green's function but are derived via a sampling matrix $B$. The tensor-product structure of the B-splines allows the inversion of the sampling matrix to be performed efficiently in $O(s^3)$ operations (where $s$ is the number of splines) rather than the naive $O(s^6)$, preserving the overall $O(s^3 \log s)$ complexity of the MVP.

2. Stabilization via Knot Removal:
   To address the ill-conditioning of the B-spline sampling matrices at high orders, the authors introduce a "knot-removal" stabilization strategy. Empirical analysis shows that removing specific knots near the boundaries of the equidistant knot vector significantly reduces the stability constant (the norm of the inverse sampling matrix) without destroying the equidistant structure required for the multilevel scheme.

   • Buffer Knots: To handle basis functions that protrude beyond their associated boxes (a common issue in BEM), the method extends the knot vector with "buffer knots." This ensures the support of all basis functions is contained within the interpolation domain.
   • Exact Interlevel Transfers: A critical requirement for multilevel methods is the ability to transfer moments between parent and child boxes exactly. The proposed stabilization strategy preserves the nested structure of the knot vectors. By using "knot insertion" (a standard spline refinement technique), the method constructs exact transfer matrices between levels. This ensures that moments computed on the stabilized (boundary-knot-removed) grid can be mapped exactly to the original equidistant grid used for kernel interpolation, maintaining machine-precision accuracy across levels.

3. Factorization Strategies:
   The paper introduces a new "star-based" factorization for the discretized EFIE, alongside existing direct and dyadic factorizations.

   • Star-Based Factorization: This approach exploits the representation of the scalar potential operator using a loop-star decomposition concept. It factorizes the scalar potential operator into a sparse star matrix and a single-layer operator. This strategy preserves the null space of the scalar potential operator up to machine precision, which is crucial for the effectiveness of low-frequency preconditioning techniques (e.g., Calderón preconditioning). It requires computing both RWG and piecewise-constant moments.
   • Direct and Dyadic Factorizations: These are reviewed as alternatives. The direct factorization interpolates the kernel directly, while the dyadic formulation uses a dyadic kernel. The star-based approach is highlighted for its superior stability properties regarding the null space.

Key Contributions

• Stabilized B-Spline Interpolation: The introduction of a knot-removal strategy that eliminates the Runge phenomenon and ill-conditioning associated with high-order equidistant B-spline interpolation, enabling robust high-accuracy compression.
• Exact Multilevel Transfers: The development of a mechanism to perform exact interlevel transfers using knot insertion, even when stabilization (knot removal) is applied. This resolves the conflict between the need for stable interpolation and the need for exact hierarchical transfers.
• Star-Based Factorization: The proposal of a new factorization strategy for the EFIE that preserves the null space of the scalar potential operator to machine precision, enhancing compatibility with low-frequency stabilization and preconditioning.
• Complexity Preservation: Demonstration that the proposed method maintains $O(N)$ complexity for MVP, assembly, and storage, leveraging the compact support of B-splines and the tensor-product structure of the sampling matrices.

Results
Numerical results were presented for both canonical (sphere) and realistic (aircraft) geometries across various electrical sizes (from $0.2\lambda$ to $4\lambda$).

• Stability and Accuracy: The B-spline-based method demonstrated stable high-order interpolation without the breakdown observed in Lagrange-based approaches. Error analysis showed convergence consistent with theoretical bounds as the polynomial degree $p$ and number of splines $s$ increased. The stabilization strategy successfully controlled errors even at high orders ($p=5, 6$), whereas Lagrange methods exhibited instability.
• Factorization Performance: The star-based factorization exhibited the most robust error behavior among the three tested strategies.
• Scalability: Performance analysis confirmed $O(N)$ scaling for storage, assembly time, and MVP time for both the sphere and aircraft models. For a sphere with $N \approx 3$ million unknowns, memory consumption ranged from 91 GB to 182 GB depending on $p$, and MVP times ranged from 7.6s to 54s. Similar linear scaling was observed for the aircraft model with up to $\approx 4$ million unknowns.

Significance
The paper claims that the proposed stabilized multilevel B-spline method provides a practical, high-accuracy alternative to Lagrange interpolation for fast integral equation solvers. By overcoming the Runge phenomenon and the ill-conditioning of sampling matrices, the method enables robust high-order compression that was previously unattainable with equidistant FFT-accelerated schemes. The preservation of the scalar potential's null space via the star-based factorization further enhances the method's utility for low-frequency problems. The authors note that preliminary results were presented at conferences and that the method can be combined with high-frequency techniques like MLFMM for broadband formulations, though the primary contribution is the stabilization of the low-to-mid frequency kernel interpolation framework.