A multilevel stochastic-gradient neural solver for boundary integral equations

Abstract: We develop a multilevel stochastic-gradient neural solver for boundary integral equations of the second kind. The unknown density is represented by a multilayer perceptron, trained by minimizing the Nyström-discretized residual on a ladder of refining quadrature grids, each level warm-started from the parameters of the previous one. Each step requires only dense matrix-vector products on mini-batches of collocation rows and network passes, operations that map directly onto GPU hardware. The residual contraction is governed by the empirical neural tangent kernel (NTK), the discrete sample of a single continuum kernel. On a fixed grid, training stalls once the residual concentrates in modes the network contracts slowly, the plateau described by the frequency principle; a spectral analysis explains, and experiments confirm, how refining the quadrature resolves more of the continuum kernel's spectrum and returns these modes to the optimizer's reach. Spectral bias, elsewhere an obstruction to neural network solvers, thus serves as the smoother of a multigrid-type iteration, with quadrature refinement in place of coarse-grid correction. Under a uniform regularity bound on the network, the total work is a constant multiple of the work on the finest grid, and the uniform conditioning of the discrete second-kind operator leaves the NTK as the sole rate-determining spectrum while converting the training residual into an a posteriori error bound. Experiments on interior Dirichlet Laplace/Poisson problems and exterior Neumann Helmholtz problems, using both parametric and signed-distance surface representations, demonstrate the effectiveness and efficiency of the proposed method compared with GMRES at comparable tolerances.

• The paper introduces a multilevel stochastic gradient method using an MLP to solve second-kind boundary integral equations via residual minimization.
• It leverages warm-starting, grid refinement, and neural tangent kernel analysis to systematically expand the effective tangent space and control spectral bias.
• Empirical results demonstrate 10–100× speedups over GMRES on complex geometries, highlighting the method’s scalability and GPU efficiency.

Multilevel Stochastic-Gradient Neural Solver for Boundary Integral Equations

Problem Formulation and Theoretical Foundations

The paper "A multilevel stochastic-gradient neural solver for boundary integral equations" (2607.00560) introduces an SGD-based, multilevel solver for boundary integral equations (BIEs) of the second kind in two and three dimensions. The unknown boundary density is parametrized by a multilayer perceptron (MLP) and optimized via residual minimization over a sequence of increasingly refined quadrature grids. The residual minimization is performed on each grid, with the neural parameters warm-started from the preceding level.

The BIEs under study arise from elliptic boundary value problems via second-kind Fredholm equations,

$$\mathcal{A}\rho := \tfrac{1}{2}\rho + \mathcal{K}\rho = g,\qquad (\mathcal{K}\rho)(x) = \int_{\Gamma} k(x,y)\rho(y)dS(y),$$

where $\mathcal{K}$ is compact, $\Gamma$ is a $C^2$ boundary, and $g$ is prescribed data. Solutions are discretized using Nyström quadrature, producing well-conditioned dense linear systems.

Traditional solvers (e.g., GMRES) accelerate matrix-vector products via hierarchical compression but still incur $\mathcal{O}(N^2)$ cost without FMM. The proposed approach positions itself for computational benefit from GPU-optimized matrix–vector operations and stochastic optimization, circumventing the need for grid-aligned prolongation and restriction schemes.

Neural Parametrization and Loss Landscape

The density $\rho$ is represented as $\rho_\theta$, an MLP. The discrete residual loss at each level,

$$L(\theta) = \frac{1}{2}\|A_N\rho_\theta - g_N\|_W^2,$$

is minimized by SGD (Adam) over random mini-batches of quadrature points. The chain rule analysis exposes a gradient flow that is preconditioned in density space by the empirical neural tangent kernel (NTK),

$T_{\theta, N} = J_{\theta,N} J_{\theta,N}^*,$

where $\mathcal{K}$0 is the Jacobian of network output w.r.t. $\mathcal{K}$1.

The convergence of the network’s output is governed by the spectral interaction between $\mathcal{K}$2 (associated with the discretized BIE operator) and $\mathcal{K}$3 (reflecting the network’s effective tangent space). The spectrum of $\mathcal{K}$4 is uniformly bounded thanks to the compactness of $\mathcal{K}$5, so the bottleneck for per-mode residual decay is entirely due to the NTK—a point established formally and borne out numerically.

Spectral Bias and the Multilevel Algorithmic Structure

Gradient-based training aligns with the frequency principle: low-frequency residual components are learned first, with high-frequency errors contracting slowly due to spectral decay in the NTK. On a fixed fine grid, this results in training plateaux.

The paper makes two key claims:

1. Spectral bias becomes a tool rather than an obstacle in a multilevel training regime: Mini-batch optimization on each coarse grid damps the residual components that are well-represented by the current NTK. Upon grid refinement, additional Fourier modes are exposed and incorporated into the effective tangent space, as the NTK spectrum expands correspondingly.
2. Full multigrid analogy: The combination of stochastic-gradient smoothing (SGS) and grid refinement implements a multilevel “band-by-band” error clearance similar to standard geometric multigrid, except prolongation is implicit via the shared neural representation, and the grid schedule is either progressive or cyclic.

This is illustrated in detail through spectral analyses.

Figure 1: Distributions showing the reduction of epoch counts at each training level, comparing the progressive and cyclic grid schedules.

Computational Details and Schedule Design

The paper implements two grid refinement schedules:

• Progressive: Sequential (coarse-to-fine), each finer level only visited once.
• Cyclic (FMG analog): Grids are revisited in non-monotonic order, enabling re-refinement and variance reduction.

For each stage, the MLP is trained to tolerance proportional to the quadrature error. SGD (Adam) is employed with mini-batch gradients; the per-iteration arithmetic is $\mathcal{K}$6 ($\mathcal{K}$7 is batch size), executed efficiently on GPU. The algorithm exploits that, unlike grid-based approaches, the neural density is a function, allowing direct transfer (“warm starting”) between node sets of different levels.

Theoretical and Numerical Analysis of Residual Dynamics

The paper provides explicit modal decomposition of gradient flow and quantifies residual transfer between levels, relying on a uniform network smoothness assumption. This yields a total work estimate matching (up to a constant) the cost of training on the finest grid alone,

$\mathcal{K}$8

where $\mathcal{K}$9 is the number of quadrature points at the finest level. The paper documents that training to matched quadrature error at each stage ensures the newly exposed residual at each grid is dominated by components newly representable by $\Gamma$0, making the convergence work nearly optimal (band-by-band contraction).

Numerical Experiments

Spectral Expansion Tracking

Experiments on a 2D flower-shaped domain demonstrate the core spectral mechanism. As training progresses up the quadrature ladder, the leading band of the empirical NTK widens, and cross-level reinitialization ensures each band of the spectrum is addressed efficiently.

Figure 2: NTK spectrum at grid level 1, showing confinement of the effective tangent space.

Figure 3: NTK spectrum one iteration after grid refinement—showing effective tangent space expansion and the inter-level transfer of high-frequency components.

Comparison with GMRES and Scalability

The solver is applied to 3D interior and exterior problems. On the Stanford Bunny and Helmholtz multi-sphere benchmarks, the MLSG achieves target residuals and relative $\Gamma$1 errors on large domains (up to $\Gamma$2), with $\Gamma$3 wall-clock speedups over multi-core GMRES. At higher wavenumbers, MLSG maintains nearly flat scaling with respect to $\Gamma$4 while GMRES slows considerably.

Complex geometries and topology (single chart coiled tori, linked tori, genus-3 surfaces) are handled without special grid construction, reflecting the flexibility of the network-based approach.

Schedule Design and Variance Reduction

Experiments demonstrate that cyclic (FMG-type) schedules result in lower variance in convergence epochs and wall-clock time compared to the progressive schedule, due to more robust spectrum refilling and residual smoothing across the grid hierarchy.

Implications and Future Directions

The MLSG validates that the spectral bias phenomenon in neural optimization, typically a barrier for PDE solvers, is rendered beneficial via grid refinement. For smooth-activation networks and quasi-uniform quadrature, neural residual minimization achieves work complexity and accuracy commensurate with the best grid-based multigrid schemes, but with algorithmic simplicity on arbitrary geometries, due to the function-based neural representation.

Practically, the method enables GPU-accelerated (and thus highly parallel) scaling to BIE problems beyond the limits of classical algorithms, especially in high-frequency and complex-geometry regimes. The trained network can furthermore provide a high-quality initialization for conventional iterative solvers when higher numerical accuracy is required.

From a theoretical perspective, the work puts the interplay of spectral bias, neural tangent kernel structure, and residual transfer on a rigorous foundation in the context of second-kind integral equations.

Conclusion

This work establishes the MLSG as an efficient, parallel-friendly alternative for large-scale BIEs, leveraging NTK spectral bias through multilevel quadrature. Both the spectral analysis and the empirical results highlight the importance of adaptive tangent space refilling and warm-started refinement in overcoming the limitations of fixed-grid neural PDE solvers. The synthesis of neural optimization dynamics and classical multilevel discretization points toward future directions in scalable implicit solvers utilizing advanced neural architectures and adaptive schedule designs.