NeurIPS 2026 (Accepted, according to the assigned list) Physics & Scientific Computing neural operators harmonic measure elliptic PDEs Walk-on-Spheres boundary integrals

Neural Harmonic Measure Operator¶

Conference: NeurIPS 2026 (Accepted, according to the assigned list)
arXiv: 2609.35752
Area: Physics
Keywords: neural operators, harmonic measure, elliptic PDEs, Walk-on-Spheres, boundary integrals

TL;DR¶

NHMO learns a geometry-only harmonic-measure density, integrates boundary data against its normalized kernel, and uses a learned lift for source contributions and approximation error, outperforming four compared baselines across all five 3D MCB-B Poisson categories while reducing repeated solves on one geometry to a cached matrix-vector product and a lift forward pass.

Background & Motivation¶

Repeatedly solving for different loads, heat sources, or boundary conditions on the same mechanical part is a natural application for neural operators. Traditional FEM/FDM requires volumetric discretization and numerical solution; Transolver, LNO, and UPT amortize solution costs but generally feed geometry, boundary values, and sources into one network to regress the solution field. This neither explicitly preserves a linear PDE's linear response to boundary data nor cleanly isolates a reusable geometry state. When coefficients leave the training range, the network can interpret a change in amplitude as a new pattern rather than a linear rescaling.

Potential theory offers a more specific object: the probability distribution of the first boundary exit of Brownian motion started at an interior query point, called harmonic measure. For a Dirichlet Laplace problem, the expected boundary value under this distribution gives the interior solution. The distribution depends on the query and geometry, not the boundary function. Compared with learning a volumetric Green's function, its density lives on “interior point × boundary point” and avoids differentiating a learned Green's function in the boundary-normal direction.

Poisson sources cannot be represented by a boundary average alone, so this is not a boundary network that solves every part of the PDE by itself. The classical balayage decomposition retains the geometric kernel channel while a separate field approximates the source contribution; the 2D realization also corrects boundary errors caused by inconsistent geometric representations. Core idea: separate the reusable geometric measure, which is linear in boundary data, from learned source and error corrections instead of treating the entire solution operator as an unstructured input-output map.

Method¶

Overall Architecture¶

Inputs comprise a representation of a bounded Lipschitz domain, boundary data, and a Poisson source; outputs are solution values at interior queries. The geometry boundary kernel encodes the shape and predicts a density over boundary points for each query. Normalized boundary integration supplies the boundary contribution, the source-residual lift supplies another term, and their sum forms the prediction. Kernel training needs only WoS exit samples, whereas lift training needs numerical reference fields: these distinct forms of supervision do not make the entire model FEM-free.

For fixed geometry and query/boundary samples, the effective kernel matrix used by normalized boundary integration can be built once and reused. The 2D lift reads boundary data, the source, and the kernel prediction to correct both source effects and boundary error. The 3D lift reads only geometry and sources, not boundary data or the kernel prediction. Dashed edges below denote training supervision, solid edges inference data flow; the kernel-prediction edge into the lift applies only in 2D.

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    G["Geometry + query"] --> K["Geometry boundary kernel"]
    W["WoS exit samples<br/>training only"] -.-> K
    K --> B["Normalized boundary integration<br/>build and cache matrix"]
    H["Boundary data"] --> B
    B -->|guide in 2D only| V["Source-residual lift"]
    F["Geometry + source<br/>also boundary data in 2D"] --> V
    R["Numerical reference fields<br/>training only"] -.-> V
    B --> O["Add to obtain solution prediction"]
    V --> O

Key Designs¶

1. Geometry boundary kernel: learn the exit distribution governing all boundary responses

Harmonic measure is a probability measure on the boundary, whereas the network predicts its density with respect to surface measure; these are not interchangeable objects. For a Laplace problem, the central representation is:

\[ u_h(p)=\int_{\partial\Omega}h(\zeta)\,d\omega_p(\zeta)=\mathbb{E}_p[h(B_\tau)]. \]

Brownian motion starts at an interior point, and its first exit determines the distribution. WoS approximates this process using uniform jumps on maximal inscribed spheres rather than many small Brownian steps. It samples the exit law without knowing the boundary function of the current problem, so training a kernel does not require regenerating labels for every boundary condition.

The network encodes boundary positions, outward normals, and Fourier features into a fixed number of slice tokens, then aggregates shape information with a Transformer. A cross-attention kernel head combines interior queries, boundary positions, normals, and displacements to predict scalar log-densities. The canonical 2D configuration uses width 256, 64 slice tokens, a 3-layer shape encoder, and a 2-layer kernel head; its kernel has approximately 4.8M parameters and the full model approximately 11.2M. The 3D kernel uses width 192 and approximately 2.77M parameters. Neither geometry encoding nor density prediction takes boundary values or sources as input: this is the prerequisite for structural reuse, not an incidental learned behavior.

2. Normalized boundary integration: turn density into a reusable convex combination

Positive density alone is insufficient: total mass must be one, or even constant boundary data will be incorrectly rescaled. NHMO normalizes density over the selected surface quadrature nodes and then integrates boundary values:

\[ K_{\theta}(p,\zeta_i;\Omega)=\frac{\widetilde K_{\theta}(p,\zeta_i;\Omega)}{\sum_j w_j\widetilde K_{\theta}(p,\zeta_j;\Omega)},\qquad u_h(p)=\sum_i w_iK_{\theta}(p,\zeta_i;\Omega)h(\zeta_i). \]

With nonnegative surface weights, each query output is a convex combination of sampled boundary values and therefore lies between their minimum and maximum. This does not establish harmonicity everywhere, exact boundary conditions, or satisfaction of the complete Poisson PDE. In particular, adding the lift removes the same convex-combination bound on the total prediction.

For a fixed sampling layout, the weighted kernel can be stored as an effective matrix. New boundary data requires only a matrix-vector product, and changing sources does not change the matrix either; changing queries, boundary discretization, or geometry requires rebuilding the corresponding matrix. The paper reports approximately 8 MB of cache per shape in its standard 2D setting. Preserving linear boundary response and enabling geometry caching are two consequences of the same factorization. NGF's geometry features also permit analogous caching in principle, so this is not an ability exclusive to NHMO.

3. Source-residual lift: avoid singular volume integration without turning a learned field into an analytic guarantee

For Poisson problems, the classical decomposition first uses a Newtonian potential for the source, then subtracts its boundary trace integrated against the same harmonic measure, yielding a zero-boundary source solution. NHMO does not directly evaluate this singular volume integral; a learned field approximates the contribution. Its implemented two-term model is:

\[ u_{\mathrm{pred}}(p)=u_h(p)+v_{\varphi}(p;\Omega,h,f). \]

The 2D lift is a U-Net receiving the interior mask, boundary data extended onto the grid, the source field, and the boundary-integral prediction. It outputs one residual channel, with the interior mask restricting the prediction region. The lift does more than fit sources: the learned boundary kernel and reference solver use different geometric discretizations, leaving an error that changes with boundary data and cannot be corrected by a network seeing only geometry and sources. Reading boundary data and the kernel prediction makes this error observable, at the cost of losing automatic linearity of the complete 2D output in boundary data.

The 3D lift uses cross-attention: Fourier-embedded queries attend to the frozen shape latent and source-sample tokens. It receives neither boundary data nor the boundary-integral prediction. Multiplying its output by the truncated negative interior SDF makes it vanish at the boundary. This more closely follows the classical source channel, but its source response remains a learned approximation, without guarantees of full PDE accuracy or extrapolation to arbitrary sources.

The paper also splits the 2D lift into an independent source lift and a boundary residual head. The former sees only the mask, SDF, and source; the latter sees only the mask, boundary values, and kernel prediction. This three-term variant explicitly distinguishes the corrections and improves accuracy, but its extra residual head remains a learned nonlinear network and its total capacity is larger. It does not prove that the complete operator is strictly linear.

A Worked Example¶

Suppose a kernel matrix has already been built for an MNIST digit domain with a hole, and a Laplace problem is solved for one set of boundary coefficients. Multiplying boundary values by the cached matrix gives the boundary contribution. Even with zero source, the 2D lift can output a boundary-error correction because the simplified polyline contour does not exactly match the numerical reference's raster mask.

Next, boundary coefficients move from the training range to the OOD evaluation range. The matrix stays unchanged and the kernel channel updates exactly linearly with the new boundary values, while the lift needs a new forward pass whose accuracy depends on generalization. Adding a source still reuses the matrix but requires evaluating the source-dependent lift again. “No retraining” does not mean “no new computation,” nor does it mean that every new boundary function and source has been experimentally validated.

Loss & Training¶

Training has two stages: first fit the kernel to WoS exits, then freeze it and train the lift with interior masked MSE between the composed prediction and reference solutions. Neither stage uses a PDE-residual loss. WoS stops at a boundary distance of 0.001 in normalized domain units and projects to the nearest boundary point; walks that do not terminate within 128 steps are masked. Finite termination thresholds, the step cap, and SDF discretization introduce approximations, so implemented supervision should not be described as perfectly unbiased.

The actual 2D setup uses 32 probes per shape, each with 10,000 precomputed walks, to form a Gaussian KDE with bandwidth 0.2% of domain width. Appendix D.1 specifies sampling 64 of 512 boundary nodes per step, normalizing both target and prediction, and optimizing KL plus 0.5 times L1 distance. It does not use mean-value, boundary-limit, or mass regularization. Kernel training runs for 60,000 steps and the standard lift for 10,000 steps.

In 3D, each step uses 8 probes and 4 fresh exits per probe. Valid-exit negative log-likelihood supervises the kernel alongside mean-value, near-boundary concentration, and total-mass regularizers. These respectively encourage agreement between interior spherical averages and center values, concentration near the corresponding boundary point as a query approaches it, and unit integral of the unnormalized density. The kernel trains for 30,000 steps; the lift uses FEM reference fields, with difficult categories reaching approximately 100,000 effective steps through warm-start rounds.

Main-text §4.5 jointly summarizes KDE/likelihood supervision and regularizers, whereas Appendix D.1 distinguishes the 2D and 3D implementations; the configuration above follows the appendix. The mean-value radius is also written as a positive multiple of the interior SDF despite the paper's negative-inside convention. The exact intended formula is not reconstructed here.

Key Experimental Results¶

Main Results¶

The 2D MNIST benchmark mixes Laplace and Poisson problems. Errors are relative L2 over all interior pixels, reported as percentages below; p95 is the 95th percentile of per-problem errors. Boundary coefficients train on [-1, +1] and are evaluated OOD on [+1, +2]. Headline results include only the poly3 and exp_mix boundary families and four source families. References use a five-point finite-difference solver at 256 × 256 and are downsampled to the models' 128 × 128 resolution.

Method	In-distribution median / mean (%)	In-distribution p95 (%)	OOD median / mean (%)	OOD p95 (%)
Transolver	4.4 / 5.3	10.7	36.4 / 36.2	55.2
LNO	5.2 / 6.4	Not reported	23.3 / 33.2	81.6
UPT	5.7 / 6.3	Not reported	26.5 / 26.6	Not reported
BENO	5.9 / 7.0	16.1	44.1 / 40.5	65.8
NGF (authors' 2D port)	2.0 / 3.9	16.9	3.9 / 4.2	6.0
NHMO, kernel only	5.4 / 7.7	18.2	5.6 / 6.8	14.2
NHMO, kernel + lift	2.0 / 2.1	3.3	2.5 / 2.6	4.2

NHMO has the lowest in-distribution mean and lighter error tails, with an OOD/in-distribution median ratio of 1.25. However, “the kernel alone beats every OOD baseline” is incorrect: it beats the four nonlinear end-to-end baselines but not NGF's OOD median, mean, or p95. NGF also combines geometry features with a read-out linear in the data and extrapolates substantially better than the end-to-end methods. This NGF row is the paper's aligned port of the official 3D implementation, not an officially released 2D experiment.

Each 3D MCB-B category contains 200 training shapes and 20 unseen test shapes, with 16 unseen boundary/source combinations per test shape, totaling 320 pairs. Errors are evaluated against FEM references at interior tetrahedral-mesh vertices. The following table reports relative L2 means without percentage units, unlike the preceding table.

Method	Nut	Gear	Motor	Fitting	Screws & Bolts
Transolver	0.320	0.281	0.407	0.180	0.221
LNO	0.372	0.466	0.528	0.259	0.239
UPT	0.516	0.507	0.765	0.392	0.358
NGF	0.275	0.243	0.338	0.160	0.189
NHMO	0.216	0.188	0.284	0.147	0.131

The four baseline rows come from NGF's table under the same evaluation protocol rather than all being retrained in this paper. NHMO is lower in every category. An additional 3D coefficient-OOD study uses only 40 problems per category: Poisson macro-means are 0.263 for NHMO and 0.615 for NGF, while Laplace-only values are 0.099 and 0.621. This is distinct from the official 320-pair test per category.

Runtime must be separated into stages. On one A100, 2D geometry construction takes approximately 6.6 s per shape, followed by 4.0 ms per cached problem; Nut takes 7.9 s + 7.5 ms and Motor 10.3 s + 7.9 ms. NGF's released Motor pipeline takes approximately 0.25 s per problem, or approximately 77 ms per forward pass with the mesh preloaded on the GPU; tetrahedral meshing separately takes approximately 48 s on the CPU. Grid-based 2D baselines have no corresponding kernel-build step, so 4.0 ms alone cannot establish an end-to-end speed advantage for one solve per new shape.

Ablation Study¶

The following results come from Appendix K.1, with mean/median errors in percentage units on the same mixed 2D test. The three-term variant uses an approximately 11.3M source lift and 6.4M residual head, exceeding the standard single lift's capacity.

Config	In-distribution mean / median (%)	OOD mean / median (%)	Note
Kernel only	7.7 / 5.4	6.8 / 5.6	No source or boundary-residual correction
Kernel + independent source lift	6.1 / 4.5	6.3 / 5.2	Lift cannot correct boundary-kernel errors without boundary data
Kernel + independent source lift + residual head	1.84 / 1.74	2.56 / 2.39	Explicitly separates source and boundary residual
Standard kernel + single lift	2.1 / 2.0	2.6 / 2.5	One network handles both corrections

Key Findings¶

On 205 Laplace test pairs, the standard lift's median correlation with kernel residuals is 0.97 and it removes 71% of the residual. Its utility is not restricted to nonzero sources. Appendix K.11 attributes the main error to mismatch between simplified polyline contours and raster masks.
Integrating the KDE targets themselves still gives 5.1% error, unchanged with 100 times more walks; WoS directly on mask geometry gives 0.07%. This favors fixing supervision geometry rather than simply increasing sampling.
In a separate 25-shape OOD quadrature scan, raising boundary nodes from 50 to 100 reduces median error from 3.12% to 2.48%; 200 and 400 nodes give 2.44% and 2.43%. These subset results should not be forced to match the headline 2.5% median.
The source contains conflicting 2D scale descriptions: Appendix G.1 states 991 / 50 / 50 shapes and 7500 / 408 / 397 problems, while K.5 and K.9 call a 5,000-shape corpus canonical. Their relationship cannot be confirmed from the cache alone; both descriptions are retained rather than reconciled by assumption.

Highlights & Insights¶

The reusable object is harmonic-measure density, not the final solution. This intermediate prediction can be checked independently using analytic harmonic functions instead of attributing all error to one black-box field.
Probability normalization both stabilizes constant-boundary response and bounds the boundary channel's amplitude. It is a local structural constraint, not a claim of strict physical correctness for the full model.
Lift ablations show that a purported “source network” can also correct discretization errors in a particular implementation. Ablations for other analytic/learned hybrid solvers should distinguish the ideal decomposition from the residual actually present in supervision.

Limitations & Future Work¶

The method currently targets Dirichlet elliptic problems; Neumann, Robin, spatially varying coefficients, and other PDEs require new measures or decompositions. Drift adaptation is a proof of concept with an additional adapter and residual head, not validation for arbitrary elliptic operators.
Linearity and boundary min–max bounds primarily belong to the kernel channel. The boundary-conditioned 2D lift can break total-output linearity; the 3D source channel does not read boundary values but does not guarantee strict source linearity or PDE satisfaction.
High-frequency trig1 and trig2 families were excluded from headline results because every learned method failed severely. The evidence concerns coefficient extrapolation within selected parameter families, not arbitrary boundary frequencies or function families.
Caching incurs per-shape precompute and storage costs and suits repeated queries on fixed geometry. Low-rank kernels could reduce construction and memory costs but introduce expressivity/error trade-offs.
Appendix J's per-problem complexity expression does not explicitly retain the general dense matrix-vector product's query-count × boundary-count term. Deployment should benchmark this operation separately rather than infer linear scaling directly from that expression.

vs NGF: NGF learns volumetric Green's functions/low-rank features and recovers solutions through volume integration and boundary read-out; NHMO directly learns boundary-measure density and amortizes the source field. Both geometry/data factorizations can help extrapolation. The distinction concerns learned objects, supervision, and source handling, not “structured versus wholly unstructured.”
vs Transolver, LNO, UPT, BENO: These end-to-end operators fuse problem data before field regression, whereas NHMO isolates reusable geometry. The 2D evidence supports this choice under the tested coefficient OOD shift, not a universal ranking on other tasks.
vs WoS / BEM: WoS re-estimates boundary expectations per query, while BEM typically discretizes boundary integral equations and solves a system. NHMO distills WoS exit distributions into reusable density. A boundary-linear residual corrector is a promising extension, but the current residual head does not impose that strict constraint.

Rating¶

Novelty: 4/5. Explicitly learning classical harmonic measure as a neural-operator object is a clear contribution.
Experimental Thoroughness: 4/5. Includes 2D, 3D, coefficient extrapolation, and extensive ablations, with high-frequency failures and scale discrepancies limiting conclusions.
Writing Quality: 4/5. Clearly bounds guarantees of the kernel and lift, but the training overview, scale descriptions, and some notation remain inconsistent.
Value: 4/5. Promising for repeated scientific-computing solves on fixed geometry; single-query new-geometry workloads require separate cost accounting.