Mitigating Instability in High Residual Adaptive Sampling for PINNs via Langevin Dynamics
Authors: Minseok Jeong
Giup Seo
Euiseok Hwang
Venue: NeurIPS 2025 spotlight ()
Keywords: Langevin dynamics
Adaptive sampling
Physics-informed neural network
Learning stability
Primary area: machine_learning_for_sciences

Abstract: Recently, physics-informed neural networks (PINNs) have gained attention in the scientific community for their potential to solve partial differential equations (PDEs). However, they face challenges related to resource efficiency and slow convergence. Adaptive sampling methods, which prioritize collocation points with high residuals, improve both efficiency and accuracy. However, these methods often neglect points with medium or low residuals, which can affect stability as the complexity of the model increases. In this paper, we investigate this limitation and show that high residual-based approaches require stricter learning rate bounds to ensure stability. To address this, we propose a Langevin dynamics-based Adaptive Sampling (LAS) framework that is robust to various learning rates and model complexities. Our experiments demonstrate that the proposed method outperforms existing approaches in terms of relative $L^{2}$ error, and stability across a range of environments, including high-dimensional PDEs where Monte Carlo integration-based methods typically suffer from instability.

Mitigating Instability in High Residual Adaptive Sampling for PINNs via Langevin Dynamics
Minseok Jeong"; Giup Seo; Euiseok Hwang
Electrical Engineering and Computer Science Gwangju Institute of Science and Technology (GIST)
dsa950115@kaist. guseo 1 100@gm. gist.ac.ki kr, euiseokt@gist.ac
Abstract
Recently, physics-informed neural networks (PINNs) have gained attention in the scientific community for their potential to solve partial differential equations (PDEs). However, they face challenges related to resource efficiency and slow convergence. Adaptive sampling methods, which prioritize collocation points with high residuals, improve both efficiency and accuracy. However, these methods often neglect points with medium or low residuals, which can affect stability as the complexity the model increases. In this paper, investigate this limitation and show that high residual-based approaches require stricter learning rate bounds to ensure stability. To address this, we propose Langevin dynamicc-based Adaptive Sampling (LAS) framework that is robust to various learning rates and model complexities. Our experiments demonstrate that the proposed method outperforms existing approaches in terms of relative L error, and stability across range of envi- ronments, including digensimmesional PDEs where Monte Carlo integration- based methods typically suffer from instability. The implementation code is publicly available attp://ggm.t...ooo/ee//////a///a
Introduction
Partial differential equations (PDEs) describe a wide range of physical phenomena, including heat transfer [17 8). fluid flow 37 18 33), wave propagation [4 5), optics, and epidemiology [26] 36].
Accurate and efficient PDE solutions are vital across many industries. With recent advances in deep learning, physics-informed
promising approach for solving PDEs. PINNs train by minimizing errors from initial conditions (IC), boundary conditions (BC), and PDE residuals collocation points [31
informed neural networks (PINNs) have emerged
These error terms are treated as soft constraints, guiding the model to satisfy essential physical requirements.
49 12
23
This collocation-bait based learning method enhances the capability of PINNs by reducing the need for extensive experimental data collection across spatio-temporal ranges, demonstrating success in various industries promising alternative to traditional numerical methods such the finite PINN (hereafter referred to PINNs) encounter challenges in efficiently setting collocation points within the constraints of a limited sampling budget and in achieving fast convergence to accurate solutions. key challenge arises from the presence of small regions with abrupt changes, in contrast larger, smoother regions. This issue is particularly evident in stiff PDEs, which are often characterized by discontinuities, such sudden transitions jumps across the spatio-temporal domain.
difference method and the finite element method [50] 25 However, collocation-bate
"Equal contribution.
Work done while GIST. Now KAIST, ACSS. *Corresponding author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Adaptive sampling in PINNs mainly follows two strategies: residual distribution-badd based methods, which resample points proportionally to residual magnitudes, and high residual- based methods, which prioritize large-residual points while neglecting low-residual regions. The latter often yields strong empirical results but risks distorting the residual distribution. This raises a key question: should adaptive sampling focus exclusively on high residuals? Addressing this question requires a careful analysis of the trade-off's and risks inherent in different sampling schemes.
Langevin dynamics-based adap- tive sampling (LAS). Through a series of theoretical and empirical analyses, we demonstrate that LAS consistently achieves reliable relative L errors and stable convergence across diverse PDEs, architectures, and hyperparameters, outperforming existing methods in robustness and scalability.
In this work, we respond by proposing a principled alternative:
Background and Related Work
Physics-informed neural networks. The basic PINN framework [35] utilizes deep neural networks as function approximators fo to estimate the solution u of non-linear PDE:
4 + [u] = O, x xE C Rd u(x; 0) = h(x),
te [0,
(2.1)
x RR:
(2.2)
u(2, t) = g(z,)) t), E 0X CRA, tE 10.7,
(2.3)
where u(xr, i) denotes the hidden solution spatial and temporal coordinates 2r, [G is the non- linear differential operator, X' is the spatial domain, OX' is the boundary, and 7 is the time range. The spatio-temporal domain is f2 X [0, T], with collocation point x (z, f) Ω. The PDE
residuals (x) and loss function on collocation points P = {x Nude C Ω are calculated as:
θ Rx(x) J. fx()++ Nolooo(xx) 8
(2.4)
Npde
Spde (1xm): Ex) = | ~ Npde
Rx
(2.5)
where ((2) is the uniform distribution over Ω and Npde represents the number of sample points of PDE loss. Then, similar manner, the total loss function defined as:
[((xtotal); ); 0) Xpde ((x,); ): + Aiofic
(xx ((*h) 0) AbcLbb((x bc );
Hyperparameters Apde- Aic and Abe control the balance between the PDE, IC, and BC loss terms. Then, fo is trained to estimate appropriate solution for PDEs by minimizing the total loss L.
Adaptive sampling based on residual distribution. Classical PINNs generally adopt uniform collocation sampling. To improve efficiency, a residual-based adaptive sampling strategy was
[Ra(x.) proposed 31,, where each point x,, is drawn with probability p(x,.) = Σn. [X..]) Building on
this idea, the residual-baa -based adaptive distribution (RAD) 47 introduces hyperparameter c: p(x) x
[x,(x) i
[[]]1 + Here, & highlights high-residual regions, while enforces degree of uniformity. By tuning these parameters, RAD flexibly balances exploration and exploitation depending on the problem structure. More recently, Gaussian mixture distribution- -base adaptive sampling (GAS) has been proposed 20), where high-residual validation points serve the means of Gaussian mixture model (GMM), and the reciprocals of residual gradients determine the diagonal covariances.
Adaptive sampling focused on high residuals. Alongside residual distribution-based approaches, another major line of work explores sampling methods that do not explicitly approximate the underlying distribution. These methods can often be viewed special cases of RAD under extreme k, c settings, but we categorize them here based on whether they estimate a sampling distribution.
A representative approach is sesidual-based adaptive refinement (RAR) 27 which iteratively selects the top-M high-residual points until the mean residual falls below a tolerance. Although effective, RAR continually accumulates points and increases computational cost. To improve efficiency, the R3 method [T0] retains high-residual points, uniformly resamples for diversity, and discards low- residual points. Similarly, cuilur-infformed PINNs (FI-PINNs) adaptively sample from regions where

---

residuals exceed threshold, using strategies such as self-adaptive importance sampling or subset simulation [3]. Beyond these sampling strategies, recent work has questioned the appropriateness of the 1,2 loss for solving Hamilton-Jacobi-ellin
(HJB) equations [4] In response, adversarial training methods targeting the 1 norm have been introduced. we discuss later, this framework leverages partial gradient information and inherently emphasizes high-residual regions in sampling.
Theoretical analysis of the concentration effect. While many adaptive sampling methods have shown promising results, several theoretical aspects remain unclear. In particular, there is a lack of theoretical analysis regarding the concentration effect. Although numerous studies report success with adaptive sampling, there is limited analysis on the impact of focusing primarily on high residuals.
To bridge this analytical gap, investigate how learning stability is influenced by the extent to which high residuals are emphasized relative to model complexity. Our analysis reveals that an excessive focus on high residuals may lead to performance degradation in PINN training.
3 Analysis of the Learning Stability
3.1 The Effect of Sampling Concentration
In this work, we define sampling concentration as the spatial aggregation of collocation points within regions characterized by large residual magnitudes. Such strategies have been shown to improve accuracy, efficiency 22 24], stability 45), and physical consistency [2 46 in PINN training.
In this section, we analyze the effects of sampling concentration theoretically.
Setup. Consider the partial differential equation defined over the domain $ X x (0, 7]. Assume that we have N collocation points forming the sample population P [ C Ω. sampled from a uniform distribution U(2).
Assumption 3.1. For analytical simplicity, we assume that the residual error of the PDE at each collocation point X,, can be expressed as a linear combination of appropriate feature map o: Ω + R... Specifically, we represent the residual error follows:
vectors, given an
Rx(xn) 01 fo (x.) + NN
(3.1)
a o(xn;0)
(3.2)
D
an aa(0))b((xx:0))
(3.3)
We regard the sampling methodology weighting of each sample point depending on the residual
Rg (xx) and set & 2. Thus, we can represent the loss function ((P:0) >> wn (x,,)) Assume that we are solving for the solution based on the gradient descent (GD) algorithm. Then,
0++1 = -
Σ wh Rox (xm))pp)
(3.4)
where the weights assigned to each sample point for iteration l are determined as follows:
[ [x
(xm))²
wh exp
NJ.
(3.5)
g2
Additionally, w' is normalized to satisfy Exep w/, (x) = 1 and the parameter 3 >> 0 preceding the residual controls the concentration of sampling with respect to the residuals. Note that the parameters Øl = (Ø' ... Ø used to calculate the importance weights do not participate in the model parameter update process. Furthermore, in contexts where the meaning is clear, we will no longer explicitly indicate that & is parameterized by 0. i.e., denote o(x; 8) o(x).
For iteration 1, we focus on two extreme cases of interest: when 3 is too large (uniform sampling), most samples receive uniform weights, resulting in uniform sampling. Conversely, when ( is close to (high residual sampling), the effect dominated by the sample with the highest residual. To explore this in more depth, consider the following propositions.
3

---

Propasition 3.1 (Steepness of the uniform sampling). When the sampling concentration parameter B sufficiently large, the maximum eigenvalue of the Hessian of the loss function can be approximated as 2)max ( where Σ = (x)) (d(x)) o x) 'I and Xmax (2) is the maximum eigenvalue of Y.
Proposition 3.2 (Steepness of the high residual sampling). When the sampling concentration pa- rameter 3 is sufficiently small, the maximum eigenvalue of the Hessian of the loss function can be [Ro(x)) is the maximum residual point.
approximated as 21|Ø(x))) )ll?, where x arg maxxep
It is well known that to ensure the convergence of GD algorithms, the learning rate 7) must satisfy the following relationship with the largest eigenvalue Amax of the Hessian of the loss function: 17 X [6. Therefore, we consequently aim to examine the relationship of the largest eigenvalue in two extreme cases of 3. Before presenting the main result, we introduce two additional assumptions that formalize phenomena typically observed in neural network-based models as their complexity increases.
Detailed proof can be found in Appendix B.I B.2
Assumption 3.2. In high-dimensional feature space, [|x(()]l follows a heavy-tailed distribution. More specifically, P(Il@(x)II > () ~ s({) for large S, where ~ represents asymptotic equivalence,
ç"
g(C) satisfies Vt > 0, lime-ox g(t() and > 0 indicates the thickness of the tail.
g(C)
This assumption is substantiated by both empirical evidence and theoretical insights. The heavy-tailed nature of feature vectors has been documented in several studies [29 30 31, and theoretically. in high-dimensional settings with complex dependencies, classical assumptions of Gaussianity often break down, and heavy-tailed models provide more accurate fit to the observed distributional behavior 4 16 40. Next, we impose an assumption concerning the representative characteristics encoded in the norms of feature vectors.
Assumption 3.3. As the model complexity D increases, the feature vector with the maximal norm becomes increasingly representative, which results in the following: maxxep []x(x)]| = [|x((*)] where x" = arg maxxep 'Ro (x))² 2 is the maximum residual point.
This assumption can be seen as a concentration of measure phenomenon in figh-dimensional spaces [1 [3). We note that Assumptions 3.2 and 3.3 do not generally hold in isotropic settings. However, they become statistically valid under (i) the low- tew-tempertture limit 3 + O, and (ii) the high-dimensional regime D >> 1 with heavy-tailed feature norms. Specifically, writing (x(x)
a(0) Φ(x), we have [x) = o(x) Q (xx) where Q a())a(0)) is rank-one. In this context, we reason as follows to determine under what conditions <(x) QΦ(x) = o(x)* (xx) holds.
In the limit 3 0, the samples concentrate on regions with large (x),, thereby biasing ((x) toward alignment with (0)) and large norm. Moreover, for D > feature vectors are nearly orthogonal, so large [Ro (x)( ]2 occurs only when d(x) is both high-norm and well-aligned. Consequently, arg maxx (x)()) x arg maxx [[x{x)]] holds with high probability in this structured regime.
To support the validity of Assumptions 3.1 through 3.3, we provide consolidated empirical veri- fication in Appendix A demonstrating their general applicability across a broad range of settings. Finally, under the scaling assumption that the number of samples N grows linearly with model size D as N = cD, we derive our main theoretical result as follows.
Theorem 3.1. Given the heavy-tailed nature of [[x(x] and sufficiently large model complexity D. we have 2|| > 2Amax(E). This inequality establishes a tighter upper bound the learning rate for ensuring the convergence of the GD algorithm under the high residual sampling method.
The detailed proof can be found in the Appendix B.3
This indicates that the stability of the algorithm can vary significantly depending on the sampling strategy and model complexity. Specifically, in these two extreme cases (3 < 1 and | > 1), uniform sampling may struggle to find an appropriate solution due to the difficulty of the stiff PDE problems, while high residual sampling may fail due to instability in the learning process.
3.2 Characteristics of Adaptive Sampling Algorithms
Building upon the previous discussion, we now briefly review how existing adaptive sampling methods operate. In particular, we focus on how the trajectory of the sample population P evolves under these algorithms.

---

Berrtionn
berution
Ieeattin
Iteration
Eteration
Berntte
Figure 1: Schematic sampling diagram of the benchmark adaptive sampling methods: (a) R3 (first three figures), (b) _ (last three figures) where [Ra(x)* 0.3 x (x:0, [[1, (0.5, 11) - 0.7 x N (x; 5, [12. - 0.3], (-0.3,2]]]
RAD 47
The modeling of residual distribution is relatively straightforward and relies on
Monte Carlo integration (MCI) over the expectation p(x) o [R)(x))"
[])a(x)]* + c. In general,
larger value of l corresponds to high residual regime with smaller whereas larger value of e indicates a tendency toward uniform sampling regime with higher 8
R3 [0]] R3 employs strategy that consistently maintains high residuals, resulting in an excessive skew in the distribution of collocation points as iterations progress. Moreover, this approach may fail to effectively handle multi-moda landscapes in the long-term, which, as demonstrated in our previous theoretical analysis, results in a scenario where the sampling concentration parameter /3 becomes extremely small.
During the adversarial training, to estimate the inner maximal value supxx2 [*) _ iteratively utilizes gradient information sign (x)R((x)) allow- ing for some degree of access local modes. However, there is guarantee that the relative proportions between modes of different heights are preserved.
L 43
To facilitate an intuitive understanding of time evolving sampling methods (R3, L*),, we have illustrated the working mechanisms a schematic diagram shown in Figure [ For a detailed visualization of the sampling trajectories, we refer the reader to Appendix D
4 Proposed Approach: Langevin Adaptive Sampling (LAS)
Similar to other residual distribution-basee methodologies, our primary objective is to estimate the residual-based sampling distribution. However, unlike previous methods that directly model the distribution using residuals, we employ Langevin dynamics to model the target distribution. An intuitive visualization of our LAS framework depicted in Figure
Langevin Dynamics in PDE Domain ?
Gradient Descent Parameter Space ?
V RR (Xnew))2
N (0,
Xnew
Vaais Orremed
θ gold
V. RRoee (f
e Onew
X
V Ro new 1
Figure Bidirectional update: (Left) In the PDE domain, our LAS framework uses Langevin dynamics to adaptively update collocation points based on PDE residuals while keeping the PINN model fo fixed. (Right) In parameter space, the PINN model fo minimizes PDE residuals with the updated collocation points.
5

---

Algorithm Single LAS Sampling Iteration for Physics-Informed Neural Networks
1: Input: initial population P =Po with Npde collocation points.
2: Output: updated population P phi
3: for [= 0 to 1 do
4: for x,, pl do
5:
Calculate the residual gradient: (x(R)(xx,)1 = Vx 3 (x No
6: Sample white Gaussian noise: z/2 ~ NN 0.
7: Follow the Langevin dynamics: 1+1 t X, + (o + BVT2.
8: end for
9: Update collocation population: pl+ +- fx 1+1 Npés
In=l
10:
end for
4.1 Langevin Dynamics and Stationary Distriburion
The dynamics of the collocation points pl C S2 at the 1-th iteration in LAS are given as follows:
- = x,
+ T Vx Rg (xx) SV
(4.1)
2
where >> O is the Langevin step size, 2 N(z,, 0, I) represents the white Gaussian noise, and B the sampling concentration coefficient. Additionally, the residual exponent is set to 2. Unlike other methods that estimate the sampling distribution based residuals at every iteration, LAS dynamically updates the data points without requiring the estimation the sampling distribution. the Langevin dynamics are allowed to run for a sufficient number of iterations IL with sufficiently small step size т, we can theoretically derive the following result regarding the collocation points.
Theorem 4.1 (Stationary distribution). For fixed fe and concentration parameter B 0, sample population pl asymptotically follows lim:-+00 (x( (x)(x) exp [Ro(x))
Although the proof is well known 9.. to complete the formulation of the distribution under con- sideration, we include the full derivation in Appendix C.1 The detailed operational procedure is summarized in Algorithm [
Unlike R3 and RAD, which use MCI to resample collocation points independently at each step, LAS refines the sampling iteratively by reusing the current population. Similar to L exploits gradient information. However, instead of relying solely raw gradients, LAS injects noise into the signal design choice whose implications are outlined below.
Theorem 4.2 (LAS favors flat residual surfaces). Given a fixed residual landscape Rg and an initial set of randomly sampled collocation points Po, the LAS framework progressively refines the sampling towards flatter local muxima while avoiding less stable and sharp regions.
An intuition-based explanation is provided in Appendix C.2 as substitute for formal proof. This phenomenon reflects a desirable property of Langevin dynamics, which inherently favors flatter regions of the residual landscape. As result, when two local maxima exhibit similar residual values, the collocation points are more likely to concentrate near the flatter one. Consequently, models trained around such flat residual regions tend to exhibit more stable learning behavior.
4.2 Key Strengths (and Advantages) of the Proposed LAS Framework
Robustness in high-dimensional PDEs. In ligh-dimensional PDE settings, the MCI-based ex- pectation (x)) N >x' (x,))) l becomes unstable due to the curse dimensionality, which demands exponentially more samples for accurate estimation. In contrast, LAS leverages local gradient information, maintaining stability even under sigh-dimensional conditions.
Enhanced stability through noise injection. Stochastic gradient descent (SGD) is known to prefer flatter minima, which are associated with better generalization 22 19 48). Similarly, LAS injects noise into the gradient signal, resulting in significantly improved training stability compared to methods that rely solely on raw gradients (4]. reasonable to expect that sampling schemes biased toward flatter regions inherently promote more stable training dynamics.

---

LAS 10urs!
RAD (Wu et al., 2023)
R3 (Daw et al, 2023)
<" (Wang et 2922)
Random
80
4roo [8]
Steepeess
60
15000
40
10000
5000
20000
40000
600u0
80oei
100000
120000
140000
160000
180000
200000
20000
40000
60000
80tri
100000
120000
140000
160000
180000
200900
Number of terations
Number of Iterations
(a)
(b)
Figure 3: With fixed learning rate / 0.002 and 4 hidden layers, (a) The maximal eigenvalue of the Hessian (steepness) for the loss function, (b) the relative L error curve.
5
Experiments
This section presents an experimental evaluation of high-residual sampling under varying model complexities and learning rates, with the number of collocation points fixed. We compare the performance of our proposed LAS method against other adaptive sampling approaches, including RAD, R3, and ~ each evaluated under the default settings provided in their original papers, random sampling with resampling (Random-R), which essentially corresponds to uniform sampling.
Experimental setup. As the default settings, unless otherwise specified, the models utilized multilayer perceptron (MLP) with 128 nodes per layer and hidden layers, employing a hyperbolic langent activation function in each hidden layer. The Adam optimizer was utilized with the learning rate of 1) = 0.001 and a decay factor 0.9 applied every 5, 000 iterations. Training was conducted with 200, 000 iterations. and the number of collocation points was set to N Npde 1, 000. For the LAS configuration, the residual exponent = 2, the Langevin step size T = 0.002 for 1-2D PDEs and 0.01 for 4-8D PDEs, and the concentration parameter 3 0.2. This hyperparameter setting represents the empirically obtained optimum, with further discussion presented in the following subsection 5.5
5.1 Ablation Studies
First and foremost, we sought to verify how the analytical results regarding stability and model complexity, presented in Section operate and apply to the functioning of each algorithm. In this context, we performed the following key ablation studies based on the Allen-Cahn equation using different random seeds.
Steepness of the loss landscape. In our stability analysis, we hypothesized that sampling algorithms targeting extremely high residuals-sud as R3 and = induce sharper loss landscapes. To test this, we tracked the maximum eigenvalue of the Hessian of the loss throughout training, as shown in Figure (a). The results support our hypothesis, indicating that high-residual- -esiiuull-ocused sampling leads increased steepness in the loss surface. The proposed LAS method, like L leverages gradient information but achieves more stable training by reducing residual surface steepness, similar to Random-R. This stability likely stems from the injected noise term. Overall, LAS outperforms Random-R and other baselines in relative I,2 error while maintaining the lowest steepness.
Different number of hidden layers. We employed MLP architectures with hidden layers ranging from 4 to 10 across all sampling methods, maintaining a learning rate of 0,001 and utilizing step scheduler. As illustrated in Figure (a), )), it can be observed that as the number of layers increases, the overall performance improves; however, for most sampling strategies at 10 layers, except for LAS, the performance diverges. particular, for R3, it evident that fails converge more prominently compared to other algorithms. These results indicate that the sensitivity to model complexity varies depending on the sampling method. In particular, high residual methods are more susceptible to increasing model complexity. For more clarity, regarding the loss curves during the training of PINNs based on model complexity, analyses are presented in Appendix not only from the perspective of depth but also from the perspective of width expansion.
7

---

LAS (Ours)
R3 (Daw et al., 2023)
RAD (Wu el. al., 2023)
1 (Weng el., 2022)
Random-A
100
Relalive error [5] (Log scale)
100
Relalive (Log
8
10
0.001
0.002
0,003
0.004
layers
Loaming
(a)
(b)
Relative error (Log scale)
100
Relative error [%] (Log scale)
100
0.00005
0.0001
0,0005
0.00075
0,0022
0.0024
0.0026
0.0028
Leaming (c)
Leaming (d)
Figure 4: Relative 2,2 error (Log scale) for the Allen-Cahn equation: (a) varying layers with 0.001 (with scheduler), (b)-(d) fixed layers with different learning rates (no scheduler). Each boxplot is based on random seeds.
Varying learning rate without decaying. We evaluated MLPs with four hidden layers across learning rates ranging from 0.001 to 0.004 without applying decay. As shown in Figuree4--br the benchmark algorithms demonstrated performance degradation 77 0.002 compared whereas LAS showed improvement. At r 0.003, all methods exhibited reduced performance; however, LAS was able partially mitigate this degradation. At r = 0.004, none of the methods produced correct solutions. In particular, we visualized the performance for very low learning rates Figure (c) and highlighted the range between n 0.002 and 0.003, where all algorithms begin to exhibit instability in Figure From this, we observe that learning does not proceed properly at very low learning rates, and for layer 4, most algorithms become unstable at learning rate as low as approximately 0.0022.
0. 0.001,
5.2 Experiments on Representative 1-Dimensional PDFs
The proposed LAS framework is further evaluated on representative 1-dimensional PDEs derived from various benchmark problems tackled by several established algorithms, including RAD. R3, _ Random-R. In these evaluations, we also employ the default experimental settings outlined earlier. The specific configurations for the PDE parameters and the hyperparameters of the baseline algorithms are detailed in Appendix F
Experimental results. We evaluated each sampling method on five benchmark PDEs Burgers' Convection, Allen-Cahm Korteweg-DD -D Vries, Schrödinger-usite five random seeds. As summarized the ID case of Table Random-R tended to outperform other adaptive sampling strategies in high-complexity model settings. Meanwhile, LAS consistently achieved either the best second-best relative T,2 errors across all cases. In particular, LAA outperformed all methods on the -len-Cahn and Schrödinger equations, while remaining competitive the others.

---

Table I: Relative L error (mean + std) across PDEs for increasing model complexity. Bold indicates best, underline second-best.
Sumpling
(Durs)
Namber of layen
8
10
3
8
10
10
10
Burgen
0.00
0.02 0.00
0.27 2 0.14
0.01 0.00
0.00
0.01
0.06 0.06
Alke-FShm
0.77 0.20
0.62 0.27
2.54 2.50
11:09: 19:85
0.99 0.29
1.36.1 0.11
0.97 0.23
34-47 17.64
0.76 1.6.1.10 0.10
10,95
19.
10
KdV
2 + 1.70
1.99 8 0.50
1.04 8: 0.6
288 1:00
7.44 :: 1,83
7.97 9: 1.45
3 243
1:00 + K.77
5.70 1.45
4.80 1.45
Schridinges
000:00001
2.89 =299=000
0.22
0.39 0/00
Convection
0.34 0.12
0.20 0.02
2: 0.05
0.41 2: 0.10
0.25 :: 0.02
0 0.00
0.30 4 0.24
0.27 4: 0.06
73,87 5.07
54.17 22.53
Burgor'
0.85 0.01
Heat
0.26 4 0.00
1.27 4 0.00
11.88 6 0101
0.20 004
1/46 4 0.17
1.54 0.15
8.7 1-7---190 1.95
5:56 + 4.30
14:79 + 1.01
13007 - 6.71
40
DF-Heca
1,72
2.14::0:13
7.73
78.91 41.66
5.72 4 0.42
73.20
37.00
2.15 0.25
80.13 a 34.24
1.46 0.67
$3.78 a
Di--ieet
3,49 0.21
11 16.6 19.14
5.53 8.081
100.00 0.00
6.99 6 1/14
100000
4 0.000
4.9 0:18
100100 4 0.00
6.10 9M
100.00 0.00
5D
DF-Heea
6.92 a 0.31
34.69 a: 38.21
19.63 a. 1.78
100,00 =2 0.00
67.50 a. 59.83
100.00
a 8.000
65.00 8: 43.01
100.00 a 0.00
13.21 a 4.46
100.00 0.00
5.3 Experiments on High Dimensional PDFs
Beyond the one-dimensional setting, we applied the existing approaches and the proposed algorithm to higher-dimensional PDE problems. Benchmark problems include 2D Burgers', , 2D heat [, and dimension flexible heat (DF-heat) 4 to 8D equations [49] These PDEs have analytics solutions, the performance of all the sampling methods could be fairly observed and evaluated in high dimensional PDE cases. Especially, we modified the DF-heat PDE to pose more challenging problem setting by increasing frequency. The details for PDEs are also described in Appendix [F) We first report the performance of all sampling methods under the default settings from the original papers (Table [].
that Random-R and LAS exhibited competitive performance. RAD, R3, and L on the other hand, experienced notable performance degradation, primarily because the 2D PDEs represent smooth cases rather than high- frequency or stiff regimes. For the 4D, 6D and 8D cases, LAS significantly outperformed all other adaptive sampling strategies. In the 4D, Random-R and 0 produced plausible solutions, but their performance remained inferior to that of LAS. Meanwhile, MCI-based methods such as RAD and R3 suffered from instability. Notably, for the 6D and 8D cases, all the other sampling approaches failed to find appropriate solutions. Although convergence was observed under different balance terms, the resulting performance was inferior and highly sensitive, as discussed in the following subsection.
Experimental results. In the 2D Burgers' and heat equations,
5.4 Further Discussion Loss Balancing, Model Hyperparameters, and Applicability
To ensure fair comparison, we conducted an extensive search for the optimal loss balance terms (Xie, Xbe: Xpdc) and hyperparameter configurations of the baseline methods, as detailed in Appendix G Additional experiments demonstrate that RAD and R3 outperform Random-R when equipped with the optimal loss balance terms and hyperparameter settings identified in Appendix G However, both RAD and R3 exhibit sensitivity to the problem dimensionality and the choice of loss weighting. In contrast, LAS consistently attains either the best or second-best relative L error across all evaluated scenarios, including smooth, stiff, dineenmnsi,nal,, and high- nigh-dimensional cases. In terms of practical application, we further evaluate robustness across different model architectures (Appendix H), factor that may be of particular relevance to practitioners.
5.5
LAS Hyperparameter Tuning and Computational Complexity
The practical applicability of LAS depends on the number of Langevin iterations 4,, which directly affects computational complexity. We therefore compare its cost with existing sampling methods in Appendix[[] [[] and the results are summarized in Figure []] The computational time and memory usage are in the order: Random-R, R3, RAD, LAS, and _ Gradient-basee based methods such as LAS and _ generally incur higher costs; in particular, L is much slower because it requires many iterations (20 by default) with re-initialization each epoch, whereas LAS needs only 4 = 1 without re-initialization. Thus, increasing the number of iterations should be considered with caution.
Detailed tuning results in Appendix J show that LAS with 4. = 1 and without re-initialization achieves competitive performance, indicating that LAS can applied in practice with computational cost comparable to Random-R, R3, and RAD.

---

LAS
RAD
R3
Random-R
10a
4.7
104
Nom 100
Nats 10, 10,000
Nose 100
Now
Myt St,
Myse i 100,
Now 20,
New So
Myw 100,
(A) 1D
dhl 4D
100
162
10a
Mayy 100
Nose 0000
ong
000
Nase 000
Now w 100
Now 000
oto
Nga 008
New 100, 000
Figure 5: Comparison of computational time (seconds per 1,000 epochs) for each method across varying problem dimensions on an RTX 4090. Empty bars denote out-of-memory failures. Subplots: (a) ID, (b) 4D, (c) 6D, and (d) 8D.
Conclusion and Future Research Directions
This paper investigates the stability of training physics-infomee neural networks (PINNs) under adaptive sampling strategies, particularly as model complexity increases. Theoretical analysis shows that methods focusing excessively high residuals may undermine stability, especially in deeper networks or with larger learning rates. To address this, we introduce Langevin dynamics-based Adap- tive Sampling (LAS), which updates collocation points using sssidull-weighte Langevin dynamics. Empirical results demonstrate that LAS ensures stable convergence even in high-dimensional PDEs, where conventional methods often fail.
Our Langevin-bart sampling scheme builds on a theoretically grounded MCMC formulation that both simple and effective. While the current approach shows strong performance, future work could explore more efficient variants for high-dimensional problems and develop principled strategies for hyperparameter selection.
Acknowledgments
This research was supported in part by the National Research Foundation of Korea(NRF) grant funded by the Korea @overnment(MSTTTyR.... -B-D)XXXS-22222004455522 and in part by the National Research Foundation (NRF) Grants, Basic Research Laboratory, under Grant RS-2023-0021880 09288
10

---

References
[1] Chandrajit Bajaj. Luke McLennan, Timothy Andeen, and Avik Roy. Recipes for when physics fails: recovering robust learning of physies informed neural networks. Machine learning: science and technology, 4(1):015013, 2023.
[2] Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116(31):15344-1553) 2019.
[3] Melih Barsbey, Milad Sefidgaran, Murat A Erdogdu, Gael Richard, and Umut Simsekli. Heavy tails in sgd and compressibility of overparametrized neural networks. Advances neural information processing systems, 34:29364-29378. 2021.
[4] Jan Beirlant, Yuri Goegebeur, Johan Segers, and Jozef L Teugels. Statistics of extremes: theory and applications. John Wiley & Sons, 2006.
[5] Umair bin Waheed, Ehsan Haghighat, Tariq Alkhalifah, Chao Song, and Qi Hao. Pin- neik: Eikonal solution using physics-informed neural networks. Computers Geosciences, 155:104833, 2021.
[6] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press. 2004.
[7] Shenglin Cai. Zhenyu Wang. Florian Fuest, Younjung Jeon, Clark Gray, and George Em Karni- adakís. BC-PINN: an adaptive physics informed neural network based on biased multiobjective coevolutionary algorithm. Neural Computing and Applications, 33(4):1 29911-19881 2021.
[8] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis. Physics-informed neural networks for heat transfer problems. Journal of Heat Transfer, 143(6);06080,, 2021.
[9] Tzuu-Shuh Chiang, Chii-Ruey Hwang, and Shuenn Jyi Sheu. Diffusion for global optimization in R SIAM Journal Control and Optimizationn, 25(3):737-753, 1987.
[10] Arka Daw, Jie Bu, Sifan Wang, Paris Perdikaris, and Anuj Karpatne. Mitigating propaga- tion fallures in physics-informed neural networks using ettin-resample-releess (13) sampling. International Conference on Machine Learning, 202, 2023.
[11] Devdatt Dubhashi and Alessandro Panconesi. Concentration of measure for the analysis of randomized algorithms. Cambridge University Press, 2009.
[12] Wenhan Gao and Chunmei Wang. Active learning based sampling for digh-dimensional nonlinear partial differential equations. Journal of Computational Physics, 475:1 11848, 2023,
[13] Zhiwei Gao, Tao Tang, Liang Yan, and Tao Zhou. Failure-informed adaptive sampling for pinns, part i: combining with re-sampling and subset simulation. Communications on Applied Mathematics and Computation, 6(3):1720-174)) 2024.
[14] Zhiwei Gao, Liang Yan, and Tao Zhou Failure-informen adaptive sampling for pinns. SIAM Journal on Scientific Computing, 1971)AA 2023.
[15] Shivam Gupta, Jasper CH Lee, and Eric Price. High-dimensional location estimation via norm concentration for subgamma vectors. In International Conference on Machine Learning, pages 12132-12164. PMLR, 2023.
[16] Laurens Haan and Ana Ferreira. Extreme value theory: an introduction, volume 3. Springer, 2006.
[17] Ehsan Haghighat, Maziar Raissi, Adrian Moure, Hector Gomez, and Ruben Juanes. A physics- informed deep learning framework for inversion and surrogate modeling in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 379:113741, 2021.
11

---

[18] Ameya D Jagtap. Zhiping Mao, Nikolaus Adams, and George Em Karniadakis. Physics- informed neural networks for inverse problems in supersonic flows. Journal of Computational Physics, 466:1 1402, 2022.
[19] Stanislaw Jastrzebski, Zachary Kenton, Devansh Arpit, Nicolas Ballas, Asja Fischer, Yoshua Bengio, and Amos Storkey. Finding flatter minima with sgd. In Proceedings of the 6th International Conference on Learning Representations (ICLR) Workshop, 2018.
[20] Yuling Jiao, Di Li, Xiliang Lu, Jerry Zhijian Yang, and Cheng Yuan. A gaussian mixture ceeributioonbaase adaptive sampling method for thysics-informed neural networks. Engineering Applicaations of Artificial Intelligence, 135:108770, 2024.
[21] George Em Karniadakis, loannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-4400 2021.
[22] Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. On large-batte batch training for deep learning: Generalization gap and sharp minima. In Proceedings of the 5th International Conference on Learning Representations (ICLR), 2017.
[23] Gregory Kang Ruey Lau, Apivich Hemachandra, See- Bee-Kiong Ng, and Bryan Kian Hsiang Low. PINNACLE: PINN adaptive collocation and experimental points selection. In The Twelfth International Conference on Learning Representations, 2024.
[24] Wensheng Li, Chao Zhang, Chuncheng Wang, Hanting Guan, and Dacheng Tao. Revisiting pinns: Generative adversarial physics-informed neural networks and point-weighting method. arXiv preprint 22Xiv:2205.0875:5 2022.
[25] Zongyi Li, Nikola Kovachki, Kamyar Azizzadeneshelii Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations, arxiv. arXiv preprint arXiv:2010.08899 2020.
[26] Shuning Lin and Yong Chen. two-stage physics-informed neural network method based on conserved quantities and applications in localized wave solutions. Journal of Computational Physics, 457:1 1053, 2022.
[27] Lu Lu. Xuhui Meng. Zhiping Mao, and George Em Karniadakis. DeepXDE: A deep learning library for solving differential equations. SIAM review, 63(1):208-2288 2021.
[28] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. Enhancing pinns for solving pdes via adaptive collocation point movement and adaptive loss weighting. Nonlinear Dynamics, 105(4):3499-350, 2021.
[29] Michael Mahoney and Charles Martin. Traditional and heavy tailed self regularization in neural network models. In International Conference on Machine Learning, pages 4284-4293, PMLR, 2019.
[30] Charles Martin and Michael W Mahoney. Heavy-tailed universality predicts trends in test accuracies for very large pre-trained deep neural networks. In Proceedings of the 2020 SIAM International Conference on Data Mining, pages 505-513. SIAM, 2020.
[31] Mohammad Amin Nabian, Rini Jasmine Gladstone, and Hadi Meidani. Efficient training of physics-informed neural networks via importance sampling. Computer-Aided Civil and Infrastructure Engineering. (6tteer)rinn.333):9662662 36(8):962-977, 2021.
[32] Kimia Nadjahi, Alain Durmus, Pierre E Jacob, Roland Badeau, and Umut Simsekli. Fast approximation of the sliced-wasserstein distance using concentration of random projections. Advances in Neural Information Processing Systems,
2021.
[33] Luis Fernando Nazari, Eduardo Camponogara, and Laio Oriel Seman. Physics-informed neural networks for modeling water flows in a river channel. IEEE Transactions on Artificial Intelligence, 5(3): 1003:10011100 1015, 2022.
[34] Chris I. Pettit and D Keith Wilson. A physics-informed neural network for sound propagation in the atmospheric boundary layer. In Proceedings of Meetings on Acoustics, volume 42. AIP Publishing, 2020.
12

---

[35] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
[36] Alexander Rodríguez, Jiaming Cui, Naren Ramakrishnan, Bijaya Adhikari, and B Aditya Prakash. Einns: epidemiologically-inonroo -nnforrme neural networks. In Proceedings of the AAAI conference artificial intelligence, volume 37, pages 14453-14460, 2023.
[37] Rongye Shi, Zhaobin Mo, and Xuan Di. Physics-informed deep learning for traffic state estimation: A hybrid paradigm informed by second-order traffic models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 540-547, 2021.
[38] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich- Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions low dimensional domains. Advances neural information processing systems, 33:7537-7547, 2020.
[39] Kejun Tang, Xiaoliang Wan, and Chao Yang. DAS-PINNs: A deep adaptive sampling method for solving high-dimensional partial differential equations. Journal of Computational Physics, 476: 1 1168, 2023.
[40] Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012.
[41] Maryam Toloubidokhti, Yubo Ye, Ryan Missel, Xiajun Jiang, Nilesh Kumar, Ruby Shrestha, and Linwei Wang. DATS: Difficult-awaaee task sampler for meta-learning physics-informed neural networks. In The Twelfth International Conference on Learning Representations, 2024.
[42] Roman Vershynin. High-dimensional probability: An introduction with applications data science, volume 47. Cambridge university press, 2018.
[43] Chuwei Wang. Shanda Li, Di He, and Liwei Wang. Is L2 physics informed loss always suitable for training physics informed neural network? Advances in Neural Information Processing Systems, 31:827888822888 2022.
[44] Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An expert's guide to training physics-informed neural networks, 2023.
[45] Sifan Wang. Hanwen Wang, and Paris Perdikaris. Learning the solution operator of para- metric partial differential equations with physics-informed deeponets. Science advances, 7(40):eabi8605, 2021.
[46] Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022.
[47] Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and ceiiuual-bsees adaptive sampling for wwwssicsinnfomed neural networks. Computer Methods in Applied Mechanics and Engineering, (03:115671, 1, 2023.
[48] Zeke Xie, Issei Sato, and Masashi Sugiyama. A diffusion theory for deep learning dynamics: Stochastic gradient descent exponentially favors flat minima. arXiv preprini arXiv:2002.03495. 2020.
[49] Shaojie Zeng, Zong Zhang, and Qingsong Zou. Adaptive deep neural networks methods for high-dimensional partial differential equations. Journal of Computational Physics, 463:11 1232, 2022.
[50] Yinhao Zhu, Nicholas Zabaras, Phaedon-Stedie Stelios Koutsourelakis, and Paris Perdikaris. Physics- constrained deep learning for high-dimensional surrogate modeling and uncertainty quantifica- tion without labeled data. Journal of Computational Physics, 394:56-81, 81, 2019.
13

---

NeurIPS Paper Checklist
The checklist is designed encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No] or [NA]. .
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justificatin right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation.
While "[Yes] " is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] "provided proper justification is given (e.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used"). ). general, answering "[No] or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS Paper Checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only the provided macros for your answers.
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: The abstract clearly states the main motivation, the proposed Langevin dynamics-based based sampling method, and its theoretical and empirical benefits. These claims are fully supported and expanded upon in the main body of the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
14

---

The paper acknowledges that Langevin sampling can introduce computational complexity, especially in high-dingg nigh-dimensional settings. To address this, we applied practical acceleration techniques and verified their effectiveness (Appendix D, D), while also noting in the conclusion that future work may develop more advanced variants improve efficiency and relative 1.2 performance beyond the foundational Langevin dynamics used here.
Justification:
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g. if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: All theoretical results are stated with explicit assumptions and are accompanied by complete proofs in the appendix B and C. We also provide empirical validation for key assumptions through dedicated experiments in appendix A.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
15

---

Answer: [Yes]
Justification: We provide all the information to reproduce the main and appendix results in Algorithm Experiments, and Appendix G sections.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c)
If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We provide all the PDE settings and instructions for implementing the main results in appendix G section. We also uploaded anonymized zip files including training and evaluation codes.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.c ec/ mablic/guidest httpl///uuides//oeeuubbnisssssss
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https //nips http://cwc..eeeeeeeee/////aaaaaa
for more details.
16

---

The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All experimental settings-incluen learning rates, optimizers, hyperparame- ters, and architectural details- are provided prior to the experiments section, with additional PDE-specific configurations summarized in Appendix G.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
We report standard deviations for key results in the main tables and provide box plots based on five random seeds to capture performance variability. For certain ablation studies (e.g., Figure 3), only mean values are presented to improve visual clarity, as the variance was not deemed critical for interpreting the observed trends. Guidelines:
Justification:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
17

---

8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: Appendix specifies the hardware setup, including GPU model and memory, and reports additional experiments measuring computational efficiency across different settings.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental rums well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http:////w.rpppp.cppccp//////tta http://wwbc...t.ccc....eeeee.eee
Answer:
[Yes]
Justification:
We do not have any violations the Code of Ethics
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The work is foundational and does not involve deployment, but its improved PDE-solving robustness may benefit scientific modeling in physics engineering domains.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
18

---

If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (c.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: The paper does not involve the release of pretrained models, sensitive datasets, tools with high risk of misuse, and thus safeguards necessary.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All external assets, including baseline implementations and benchmark datasets, properly cited with source papers. Their usage complies with the respec- tive open-source licenses such as MIT or CC-BY.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, pawerssipppswiiiieea. thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provide structured templates for the details of implementing our algorithms. All authors consent to provide documentations.
19

---

Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve any human subjects or crowdsourcing in data collection or evaluation.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve human subjects, crowdsourcing, or any study that requires IRB approval.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
20

---

Justification: An LLM was only used for writing and editing support: it was not involved in any part of the method development or scientific contributions.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non non-standard components.
Please refer to our LLM policy hattps://neurips.o 20tps/////ooeeeeeeee/22///////// for what should or should be described.
21

---

Appendix
A Empirical Validation of Assumptions
By definition, the residual Rs involves transforming the neural network output fo through operator within certain function space, This implies that the features the function, resulting from the neural network output combined with an additional operator, not explicitly defined. Consequently, we extract the feature vector of the resulting function through local feature vector estimation based on the linearization of the residual function. Before delving into the main discussion, we first explain the logic behind how the feature vector & is inferred.
A.1 Local Approximation of the Feature Vector
Let (x(x)= a fx( + IN represent [x()))((x). To validate the assumption that a suitable 8t linearization exists, our goal derive a proper linear approximation of (g(f))((x at a specific point x E Ω = X' x [0, T given a specific function fe. In this process, we will utilize a Taylor expansion for (x). It is important note that since g(fo) represents the behavior in the function space, understanding how responds to small perturbations in So is crucial. This analysis employs the Fréchet derivative.
To summarize briefly, s(f++ae) ~ g(fe) + (f(f)) - fa), which implies that the result can be linearized around a baseline function fo where Do(fo) limae-o Ila(Se ||:0|||
- Since our
focus is on the linearization of (x), it is essential to ensure that fa+ is a function close to So within the function space. To achieve this, small noise perturbations A0 are added the neural network Sa- In conclusion, to approximate the value at a specific point x, we proceed as follows:
[0([++>a)(xx) xs (s(Sa) + D, (fo) fa)](x)
(A.1)
= (x( + [D, fo)](x)
(A.2)
(x)
(x) fo(x))
(A.3)
Here, if fo is assumed to be a well-trained PINN model and perturbation A6 sufficiently small, we readily infer the following for the first term:
[o( [x(fo)](x) ~ O, vX
Consequently, the linear approximation of the function [a(Jo+ao)] can be expressed using the Fréchet derivative. The aspect that conflicts with assumption is that, in this context, the Fréchet derivative can act as a function dependent on x. Therefore, we refer to this as a local approximation.
Approximation of Fréchet derivative. According to the problem formulation of PINN, 9 is an operator that takes the function f as input and generates new values through partial deriva- tives such as 0ef, fuf Ōtzf, and their combinations. Thus, we can assume g(f) G(f, 0.f. fzf, Ouf. O Here, G is a multivariate function that combines the derivative terms. Next, considering scenario where slight perturbation Af is applied to fe. the Fréchet derivative can be approximated as follows:
Dy(fe+so) (fetse- fo)
(A.4)
~ x(fo+so) g(fo)
(A.5)
= G + Afo. 0.5o + Adrfo, 0.So 10 AQ.Jo, -.lfoo..))GG(ffoofffoo..
(A.6)
DG Afo a(afo)
DG Ad. 0(0.fo)
8G
Ad +...
(A.7)
O(fo}
Σ a(axfo)
8G Adefo,
k (0, 2, tt, tz. ...
(A.8)
a(0)To To.
(A.9)
ac where a(0) = ( J... N...J ... and Φ = (Afo. Ad.fo. AO.Jo....)
Mathematical details. Here, we provide a systematic summary of the considerations underlying the validity of the employed estimation method.
22

---

1. Reliability of the Fréchet derivative Dylf]: The existence of the Fréchet derivative requires the following sufficient conditions:
G must be differentiable.
must be sufficiently differentiable with respect to z and t.
Both conditions are naturally satisfied in the context of our problem. This ensures that we can extract vector that locally approximates the actual feature vector & for each point x E S2, thereby facilitating robust estimation process.
2. Condition for the constancy Da[ is important to note that G' generally function of "...". ), and thus implicitly depends on However, when the variables are not entangled with each other, the partial derivatives can exhibit constant behavior. For instance:
8G is always 1. In our case, J(0,J)
aG aG Partial derivatives in the z-direction such a(H,f)* (02af] depend on Nu
If the output of the differential operator NV entangles the partial derivatives in the x-direction (i.e, N²af[ is non-linear), the assumption that a acts as a constant may weaken.
We assume that (0) is constant or locally linear, enabling unified lincarized form of the residual. This simplification ensures consistency across Assumption 3.1 and Propositions 3. 3-122.2, and un- derpins curvature analysis. We will revise the manuscript to explicitly include this structural assumption in Assumption 3.1 and ensure that its implications are consistently reflected in all related propositions.
To provide further intuition, we illustrate this within the context of two representative PDEs.
For the convection equation,
g(f) = 0,f + hdff
a()) = (0, 1, 0,0,.... is constant.
For the Burgers' equation,
0.01
a(0) = OnJo. 1, fa, 0, 0, 0, 0.01
s - 0 + f0.
Ozzs
is O-dependent.
Under the assumption that the neural network is sufficiently expressive and smooth with respect 0, we can locally linearize a()) around a fixed 0.
When a()) is affine, its Jacobian
|aa;
Jo(a(t))
Ja
a0,
is constant and its Hessian
Ho(a(O)) 0,
where the derivatives are defined elementwise.
A.2 Heave-Tailed Behavior of the Norm of Feature Vectors
Initially, we visualized the histogram of the norms of the extracted feature vectors across all feasible grid points in Figure i.e., the histogram of {[]((x)) x Ω X' x [0, T]} for models with 4, 6, and 10 layers, respectively.
From the provided histograms, it is evident that for each PDE, the distribution increasingly exhibits heavy-tail behavior as the layer depth grows. This tendency is particularly emphasized in the following two aspects:
Heavy-tail characteristics resembling Pareto distribution: As the layer depth increases, the distribution's 's tail becomes thicker, consistent with the heavy-tail properties of the Pareto distribution. Pareto distribution, the tail probability follows the form P(X > x)) ox x = decaying slowly and exhibiting high frequency of extreme values. This is reflected in the histograms, where deeper layers show data concentrated certain regions while displaying more frequent extreme values.
23

---

2. Increased concentration and frequency of extreme values: As the number of layers increases, the data become densely concentrated within specific ranges (represented on the y-axis as frequency), while significantly more frequent occurrences of large values (depicted on the r-axii as extreme values) are observed. This behavior suggests a progressive shift towards heavy-tail distributions.
In addition to the previously obtained histograms, we also calculated two statistical estimates- Pareto tail index and Hill estimator-baat
the samples to provide more quantitative representation. Pareto tail index. The Pareto tail index, denoted by a, quantifies the heaviness the tail of distribution. For a random variable X with heavy-tailed distribution, the tail probability follows power-law:
P(X 2) ~
oC,
where O represents the tail index. Therefore, smaller value of corresponds thicker tail, indicating a slower decay of the tail probability and a higher likelihood extreme events. Conversely, larger value of corresponds to thinner tail, where the tail probability decays more rapidly and extreme events are less likely.
Hictogram of Lager
Hictogram Layer
Hictegran of Layer
Hietegrams of Layer 10
250
Margene
205
1881
Hirtogran Lapar
Hirtogram Laper
Hattgraa Layer
Histegran Layer
Cosuection
203
342
186
Hirisgram Layer
Hirtogram Layer
Harlogran Layer
Hirlsgram Layer
302
MiveCaan
102
Hivtegram Loyer
Hirtogran
Mirtegrm Layer
Hirlegram Layer
Do Vitee.
543
300
363
12
Histegram
J52
Schoddigu
103
120
22
Figure 6: Each plot represents a (PDE, layer) pair, where the row corresponds to the type of PDE being solved (e.g., Burgers', Convection, Allen-Cahn, etc.), and the column indicates the model size by the number of hidden layers the PINN (e.g., layer 4, 6, 8, 10). The histograms show the distributions of the feature vector norms Φ(x)] for each pair.
24

---

Hill estimator. The Hill estimator is specifically designed to estimate the inverse of the tail index, <= 1. Given sample of independent and identically distributed observations X2......... N,)=
sorted in descending order as ((1) 2 X ... 2 (n)) the Hill estimator is defined as:
X
E
log X(k+1)
where & is the number of upper order statistics used for the estimation.
Figure shows box plots of estimates across random seeds for various PDEs. The Pareto tail index decreases with depth, indicating heavier tails and a higher likelihood of extreme events. Across all PDEs, the index drops consistently from layer 4 to 10 and remains well below the heavy-tail threshold of 2. Even with fewer layers and 1,000 collocation points, feature norms display clear heavy-tailed behavior. The Hill estimator, the inverse measure, increases with depth and exceeds the 0.5 threshold, further confirming the trend.
Poreto
3:8
matimattor
J00F' Lipe's 3pd
Upt LiW''s
Figure 7: Two statistical estimates of the norms of the feature vectors.
A.3 Emerging Disparities with Increasing Model Complexity
In the previous subsection, we conducted an empirical analysis of the distributional characteristics of feature vector norms. In this subsection, aim to validate the hypothesis of the norm's emerging disparities with increasing model complexity (Assumption 3.1 For clarity, this relation can be
expressed mathematically as (maxxep [[)) Amax (Σ) > (maxxep [[)(x))))* - [|ø((*))l
where x arg maXxen [x) Here, due the dominant scale of max ( |] we transformed
the values into logarithmic scale to investigate the relationship between Amax and (xx((") )l
90
Layin
Lapent Layer Lopett
LayiTe
Topen Loger La
Loger
Figure 8: Disparity comparison of Amax (2) and o(x*) with respect to max (x)
The Figure illustrates the behavior of a logarithmic metric for various PDEs as the layer count increases. The x-axis represents the number of layers, shown as 4, 6, 8, and 10, while the y-axis represents log-based value, denoted as Jog(max (e(())))= = [|x(x")]]
log[max (x) 1 log (20) which captures ratio involving maximum values and scaled terms. Across all PDEs, the y-axis value decreases monotonically as the log-metric across all PDEs suggests that as the layer count grows, the denominator in the ratio scales disproportionately compared to the numerator. This behavior indicates that the underlying system dynamics or representation becomes increasingly dominated by the factors represented in the denominator.
number layers increases. This consistent decline
25

---

B Learning Rate Upper Bound Varying B
For the sake of simplicity, we will consider the situation iteration Hence, in the forthcoming proof, we will omit the upper index related to iteration. i.c., denote w as wn. Under the assumption that the neural network fo is sufficiently expressive and smooth with respect to 8, we can locally linearize (f) around fixed B.
By the Assumption 3.1, when a()) is affine, its Jacobian J)(a(()) = its Hessian H& (o)) where the derivatives defined elementwise.
aa. 3e,
= J. is constant and
We consider the weighted loss
(( =
Wn a T 2n a = a(8), on == ((x,).
(B.1)
Then the Hessian is
N
D
Hg 2
Σ
wn
(n) (I 2n)
(aT Σ
Una Holai)
(B.2)
n=1
If a a(0) is linear, the second term vanishes, and we obtain
N
HoC =
Σ
W On
J.
(B.3)
B.1 Proof of Proposition 3.1
Proof. When the sampling concentration coefficient /3 is sufficiently large, the weights w. are approximately uniform (w,, R }) Thus, the Hessian matrix Hg of the loss function with respect to & can be approximated by:
N
Ho(C) = Ho
Σ
w.n (oo.|R((x)) = Ho
Ro (x(xx))2
(B.4)
N
==1
Since Rx(xx) a TT and Hg (laTo,)) fies:
o.)
2J (ENY
Wn 2n
J. the Hessian of L satis-
N
Ho(C) R
1 L Ho (la 1n N
2
Σ
J
(B.5)
N
This matrix approximates the product of the Jacobian and the sample covariance matrix of the feature vector oxx). Therefore, for sufficiently large N, the maximum eigenvalue of the Hessian can be approximated as 2Amax(2)
B.2 Proof of Proposition 3.2
Proof. When the sampling concentration coefficient @ is sufficiently small, for the sample x" with the largest residual (i.e., x* = arg maxxep. 'Ro (x) )) we can consider all other weights to zero except for x". Therefore, the Hessian matrix of the loss function be expressed as follows:
Ho( Hg
Σ wn Wull Ro(x))) Ho ((Ro(x"))
(B.6)
The Hessian of the loss function can expressed as 2J (xxx))(xx)] Therefore, in similar manner, if we aim to compute Amax here, we may consider the eigenvalue equation Av = Xe Assuming J. is independent that (xx) is an eigenvector of A, with the corresponding eigenvalue given by 21Ø((*)) )2 Since A is rank-1 I matrix, this eigenvalue is uniquely determined.
the dimension D. let A (x)))((x) and v = d(x). follows
26

---

B.3 Proof of Theorem 3.1
Proof. Assume that the feature norm (xx(x) follows heavy-tailed distribution with tail index a > and that the number of samples N satisfies N = D.
Step 1 (Low-temperature regime). When the loss is dominated by a few extreme samples, the empirical Hessian satisfies
HgC ~ Wmxx Ja @max T
max Ja:
Pmax
max ISnEN
Ilo(x,)ll
(B.7)
By extreme value theory, the maximum feature norm scales as
Pmax x N1/a D1/0
(B.8)
hence
max x D2/0
(B.9)
Step 2 (Uniform sampling regime), In contrast, under uniform sampling without dominance by extremes, have
HoL 2J EJ.: -=B(o(x)o(x))
(B.10)
where E is the population covariance matrix of features.
Step 3 (Scaling of covariance cigenvalues). For a > 2. Σ is finite, and in normalized isotropic setting Amax ( O(1). Considering instead the sample covariance
Σ (xx))(x))T
(B.11)
results from random matrix theory for 2 < a imply
Amax 4
max
N2/a-1 D2/a-1
(B.12)
Step 4 (Comparison). With N x D, obtain
m x 2
and
Amax (#N) D2/a-1
(B.13)
Therefore,
mmax
x D 00 as D
(B.14)
Amax
NNN
which shows that 2h Amax ( ) for large D.
Finally, combining this fact with Assumption 33.1 which states ||o(x*)] Qmaxs where x* maxxer (2) we conclude the proof.
27

---

C Stationary distribution and Stability of Langevin Dynamics
C.1 Proof of Theorem 4.1
Proof. Starting from the discretized Langevin dynamics:
++ x² + Vx|Ro(x))1 + 3VFd, It N
(C.1)
We consider the continuous-time Langevin dynamics as the step size - 0.
2 1 Vx (x),)) ds + 3dB.:
dx.
(C.2)
where B, denotes standard Brownian motion in Rd and B > 0 is the diffusion coefficient. We assume the drift x Ro (x))
(x) locally Lipschitz and satisfies linear growth condition, ensuring the
well-posedness of
Step 1: Fokker-Planck equation. Let (x( denote the probability density of xs. The Fokker- Planck equation corresponding (C.2) is
Op,(x)
32
Vx
(x)R(( px (x(x))
+ AxPo(x).
(C.3)
Ōs
2
where Ax denotes the Laplacian.
Step 2: Stationary equation. At stationarity, 0sps = 0, so
C.3
becomes
0= Vx (~Vx|R((x) Pae (x)
B
+
AxPoo (x),
(C.4)
2
where Poo is the stationary density.
Step Candidate Gibbs-type solution. We claim that
[x)
(x)
exp
(C.5)
Z
B2
where Z is the normalizing constant, solves (C.4.
Step 4: Verification. Differentiating CC.5
1
(x) = Vx (x)) Pro (x),
(C.6)
g
(x) = g2 Ax Ro(x))
+
[<) (x).
(C.7)
B
Substituting (C.6 and C.7
into C.4
gives
g2
1?]
Ax Ro
Pro
+ Rol
+ VxPeo]
2
g2
Rol
+
Vx
Vx
Pso
Ax 22
+ g2
V
2212
Poo +
Ax
Rol
1
R
121
R
Pro
2
2
0. Thus Poo satisfies
C.4 exactly.
Step 5: Conclusion. By uniqueness of the stationary solution under the given assumptions, the stationary density of (C2 is given by (C.5.
28

---

C.2 Proof of Theorem 4.2
Proof. Rather than presenting a formal proof, we illustrate the underlying intuition using a simple toy example.
Preference for flat high-residual regions under Langevin dynamics. For fixed neural network parameter 0, we denote the residual as Ro = Rxxx and x for notational simplicity.
We analyze the residual function:
R(x) = -aaxx-60 2 eexp (-6x (x- 8)2 (x==))
(C.8)
2x²
which has two local maxima: a broader maximum at Xitat = and narrower, sharper maximum at Xaharp = 0 due to the additional Gaussian term.
We consider Langevin sampling governed by the stationary distribution:
iT (x) exp (pR(x)),
(C.9)
where p > the inverse temperature. This distribution favors regions with higher residual values.
To compare the probability sampling near each local maximum, we approximate the local proba- bility mass by Laplace's method. Around each maximum x", we expand the residual as:
R(x) = R(x") + R" (x" (xxxx))
(C.10)
where (x") and the local curvature is defined as = -x"(xx)) > O.
The probability mass near local maximum x* is approximated by:
J
2T
exp(pR(x)) dx x exp(pR(x')) Dear
(C.11)
pH
Let:
Hilst - R" (x((a))) = 8af?
(C.12)
Haharp R" (Xsharp) = 8a02 + σ²
(C.13)
so that Hita < Haharp for small a.
Assuming the residual values at both peaks are equal, i.e., R(x(s)) = R R(Xsharp), the ratio of local probability masses becomes:
(near Xftat)
Haharp Haat
(C.14)
P(near Xsharp)
After normalization, the respective probabilities are:
1/VHaat
Haharp
(near X(lat)
(C.15)
1/v WHat + 1/ Haharr
Hiat +
Habarp
Vi
(near Xsharp)
(C.16)
WHat
/Hsharp
This analysis shows that, despite equal peak heights, Langevin dynamics assign higher probability to the broader maximum Xaat due its lower curvature and larger effective volume in the residual landscape.
For example, with parameters 1, 2, 10, and 0.1, the curvature values Haat 32 and Haharp = 1032, yielding a probability ratio of approximately 5.68. After normalization, we obtain:
P(near X(tat) 2 0.85, P(nnarr Xsharp) x 0.15.
This quantitatively confirms that Langevin dynamics significantly prefer broader residual peaks under equal height conditions.
29

---

Gradient sensitivity induced by sharp residual regions. Residual peaks with high curvature, such those near Xgharp: induce steep gradients with respect the model parameters due the amplified sensitivity of R to perturbations in x. Formally, the chain rule
d Ro(x)=
0R dx
d6
dx
d#
indicates that when 0R/00 0 is large as in highly localized peaks even small variations in x can result unstable updates in 6.
Moreover, a concentration of samples around sharp residual peaks effectively sharpens the empirical loss landscape, making the optimization process highly sensitive to initialization and learning rate. In contrast, broader regions such as those around Xaat yield smoother gradients and reduced curvature in the loss surface, enhancing the robustness of training against both gradient noise and model perturbations.
This observation highlights that the inherent sampling bias of Langevin dynamics toward flatter high-residual regions contributes not only to generalization but also to numerical stability during optimization.
30

---

D Sample Trajectories of LAS and Benchmark Algorithms
We visualized the collocation point trajectories during the training under various adaptive sampling algorithms. The experimental settings follow the default configurations specified in the main text, with 4 layers and the number of Langevin iterations heatmap using the plasma colormap, represents the solution errors, where dark purple indicates low values and bright yellow indicates high values. White points represent the collocation points aligned with solution errors to intuitively show how each algorithm works.
4 The background, shown
D.1 Random-R Sample Trajectory
The figure below represents the sample trajectory of Random-R, R, where different collocation points uniformly sampled each iteration.
Iteration 1000
Iteration 2000
Iteration 3000
Iteration 4000
Iteration 5000
Iteration 6000
Iteration 7000
Iteration 8000
lteration 9000
lteration 10000
Iteration 11000
Iberation 12000
Iteration 13000
Iteration 14000
Iteration 15000
Iteration 16000
Iteration 17000
Iteration 18000
Iteration 19000
Iteration 20000
Iteration 21000
Iberation 22000
Iteration 23000
Iteration 24000
Iteration 25000
Iteration 26000
lteration 27000
Iteration 28000
lteration 29000
Iteration 30000
Figure 9: Random-R sample trajectory.
31

---

D.2 RAD Sample Trajectory
This figure presents the RAD sample trajectory over multiple iterations, demonstrating a relatively stable pattern of sample distribution. To ensure fair comparison, we fixed the number of collocation points for the denser set to be pde- As iterations progress, the sample points concentrate around regions of high errors, with some diversity maintained throughout. However, despite the overall stability, the RAD sampling method exhibits a distribution that not significantly different from the Random-R approach. The clustering becomes more pronounced in certain areas, but the overall spread and distribution of samples remain similar, suggesting that RAD does not offer a distinct advantage over random-R R in terms of improving sampling diversity.
lteration 1000
Iteration 2000
lteration 3000
Iteration 4000
lteration 5000
Iteration 6000
Iteration 7000
Iteration 8000
Iteration 9000
Iteration 10000
1.3
Iteration 11000
Iteration 12000
Iteration 13000
Iteration 14000
Iteration 15000
14
Iteration 16000
Iteration 17000
Iteration 18000
Iteration 19000
Iteration 20000
Iteration 21000
Iteration 22000
Iteration 23000
Iteration 24000
Iteration 25000
Iteration 26000
Iteration 27000
Iteration 28000
Iteration 29000
Iteration 30000
Figure 10: RAD sample trajectory.
32

---

D.3 R3 Sample Trajectory
This figure illustrates the evolution of sample trajectories in the R3 algorithm, showing clear concentration of samples in regions with high errors as the process progresses. While early iterations exhibit some scattering. the sample points increasingly cluster around specific areas of the solution error, leading to lack of diversity in later stages. Furthermore, this imbalance indicates instability the sampling strategy, as it fails to maintain a continuous, balanced shift in the sample population across the entire domain. The discontinuous change in the sample population may result in instability from the perspective of the learning process.
Iteration 1000
Iteration 2000
lteration 3000
lteration 4000
Iee ation 5000
Iteration 6000
Iteration 7000
lteration 8000
Iteration 9000
Iteration 10000
Iteration 11000
Iteration 12000
Iteration 13000
Iteration 14000
Iteration 15000
10
Iteration 16000
Iteration 17000
Ieeration 18000
Iteration 19000
Iteration 20000
14
Iteration 21000
Iteration 22000
Iteration 23000
Iteration 24000
Iteration 25000
10
00
Iteration 26000
Iteration 27000
Iteration 28000
Iteration 29000
Iteration 30000
Figure 11: R33 sample trajectory.
33

---

D.4 L Sample Trajectory
This figure illustrates the evolution of sample trajectories in the _ algorithm. As the number of iterations increases, the samples become overly concentrated in regions with high errors, leading to lack of diversity across the domain, particularly in areas with lower errors. This imbalance goes against the goal of maintaining well-distributed sample set proportional to the solution errors. While some adaptation occurs, the excessive focus on extreme errors (small 3 case) results in a skewed distribution, highlighting the need for more balanced and diverse sampling to improve the algorithm's s performance.
lteration 1000
Iteration 2000
lteration 3000
lteration 4000
Iteration 5000
Iteration 6000
Iteration 7000
Iteration 8000
Iteration 9000
Iteration 10000
Iteration 11000
Iteration 12000
Iteration 13000
Iteration 14000
Iteration 15000
10
Iteration 16000
Iteration 17000
Iteration 18000
Iteration 19000
Iteration 20000
8
Iteration 21000
Iteration 22000
Iteration 23000
Iteration 24000
Iteration 25000
00
Iteration 26000
Ineration 27000
Iteration 28000
Iteration 29000
Iteration 30000
Figure 12: I sample trajectory.
34

---

D.5 LAS Sample Trajectory
This figure depicts the sample trajectory of the proposed LAS algorithm. As iterations progress, the sample points are proportionally distributed according to the solution errors, maintaining diversity across the domain. Unlike other methods, proposed LAS algorithm avoids over-concentration in regions of high errors, instead ensuring that sample points are scattered in balanced manner. Additionally, the distribution adapts in line with the error peaks, with an appropriate portion of samples allocated based on the peak heights. This indicates that the LAS algorithm successfully addresses the key objectives of both proportionality and diversity sample distribution, improving stability and overall performance.
Iteration 1000
Iteration 2000
lteration 3000
Iteration 4000
lteration 5000
Iteration 6000
Iteration 7000
Iteration 8000
Iteration 9000
Iteration 10000
1.0
10
0
Iteration 11000
Iteration 12000
Iteration 13000
Iteration 14000
Iteration 15000
Iteration 16000
Iteration 17000
Iteration 18000
Iteration 19000
Iteration 20000
00
Iteration 21000
Iteration 22000
Iteration 23000
Iteration 24000
Iteration 25000
Iteration 26000
Iteration 27000
Iteration 28000
Iteration 29000
Iteration 30000
Figure 13: LAS sample trajectory.
35

---

E Relationship Between Learning Rate and Model Complexity
E.1 Learning Curve Variation with Increasing Depth
Here, we aim to visualize and interpret the learning curves for the Allen-Cahn equation observed during the training process of the models, as reported in Table with detailed experimental settings provided in Section 5] varying only the depth of the neural networks.
LAS (Ours)
RAD (Wu at al., 2023)
R3 (Daw at al., 2023)
L (Wang et al., 2022)
Random- R
Rei ative L. error
100
50
20000
40000
60000
80000
100000
120000
140060
164080
189000
200000
Relatve errrr (%)
Number of iterations
100
50
20000
40000
60090
80580
100000
120000
140000
160000
180000
200000
Relative L2 error [%]]
Number of iterationn
50
26000
40000
60000
80000
100000
120000
140000
160000
130000
200000
Re ative error (86
Number of iterations
100
50
20000
40000
60000
80000
100000
120000
140000
160000
180000
200000
Number of iterations ons
Figure 14: From top to bottom in the figure, the learning curves correspond to layers 4, 6, 8, and 10 for the Allen-Cahn equation.
Our primary observation is that most algorithms exhibit slow learning progression until the learning rate reaches a specific value (which, of course, varies depending on the algorithm). This phenomenon appears correlate with the degree of high residual concentration in the residual landscape of each algorithm. Specifically, the relative 7,2 error in the learning curve requires more iterations to drop below certain threshold (denoted as 50 in the figure) as the number of layers increases.
From an algorithmic perspective, most methods achieve the 50-threshold the relative L2 error crossing before iteration 40,000 with a layer network. However, as the number of layers increases, particularly with 8 and 10 layers, the preeshoddcrossingg delay is especially pronounced in algorithms such as R3 and L which highly focused on regions of extreme high residuals. This observation suggests that these algorithms are more affected by the increased complexity and residual concentration in deeper networks.
iterations are significantly delayed. This
To verify whether this phenomenon depends on overall model complexity, in the following subsection, also conducted experiments focusing increasing the width rather than the depth of the model.
36

---

E.2
Learning Curve Variation with Increasing Width
Simple calculations show that a neural network with hidden layers and 128 nodes per layer has the same number of parameters as neural network with 4 hidden layers and width of 203. However, comparisons based solely on parameter count are inadequate, depth introduces issues such as gradient vanishing.
Therefore, instead of viewing width solely from the perspective of parameter count, conducted experiments by progressively doubling the width. The results showed that, similar to depth, increasing width also led a gradual breakdown in learning stability. Consistent with the rankings observed in depth experiments, among adaptive sampling techniques, LAS reached the relative L error threshold of 50 the fastest. Interestingly, Random-R demonstrated robustness in this setting, particularly with wide neural networks.
[Ours)
RAD al.,
2023)
La (Wang al., 2022)
Random-R
Rei Relative arror [%]
100
20000
40000
60000
80000
100000
120000
140000
160000
189000
200000
Number of iterations
Relatve Lz aroor
50
20060
40000
60090
80060
100000
120000
140000
160000
180000
200000
Relative error [%)
Number of iterations
100
50
20000
40000
600g0
86000
108000
120000
140000
160000
1800g0
206000
Relative L2 error [%)
Number of Iterationn
20060
40000
60000
80000
100000
120000
140000
160000
180000
200000
Number of iterations
Figure 15: From top bottom in the figure, the learning curves correspond to 128, 256, 512, and 1024 nodes per hidden layer, each with 4 layers, for the Allen-Cahn equation.
Through the experimental analyses deseribed above, argue that the proposed LAS demonstrates superior performance in terms learning stability, particularly for models with high complexity.
37

---

Supplementary Details on Experimental Setup
F1
Details of Partial Differential Equations
Burgers':
We set (Aie, Abc Apde) = (100, 1, 1) to solve the equation
Ou du 0.01 0.u = 0, a x11--111], te [0,1];
0I + u dr
(F.I)
Ox²
u( - 1, 1) = u(1,))
(F.2)
u(r, 0) = sin N..
(F.3)
Convection: We set (Aie, Abe: Apde) = (100, 100, 1) to solve the equation
du
du
+ dr
0,
x E [0,22t] 2T),
E [0, 1];
(F.4)
Øt
u00,
u(2T, t);
(F.5)
w(,, 0) sin J.
(F.6)
Allen-Cahn: set (dic; Ane Apa))
(100, 1, 1) to solve the equation
du 01 0.0001
82u
5(u - u²) - O, x 1E1-1,1], E [0,]]
(F.7)
Ox.
t) )) = u(1,t);
(F.8)
u, u (-1, uz (1,
(F.9)
u(z, 0) x² cos πT.
(F.10)
Korteweg-DD De Vries: We set (Dic, Apde Abc) = (100, 1, to solve the equation
Ou
Ou
03u
Ot + u dr
+ 0.0025
= 0, x 1-111,1] te [0, 1];
(F.11)
3z
u(=1,t) = u(1,t);
(F.12)
u(x, 0) cos N...
(F.13)
Schrödinger: We set (Aies Abe, Apul) = (100, 1, 1) to solve the equation
Oh 8hh
0.5
+ |1|21 0, - t oo, 2
(F.14)
Øt Or
h(-5,t) h(5,t)) t);
(F.15)
h=(-5,t) = he (5,
(F.16)
h(2, 2sech(x).
(F.17)
Heat 2D: We set (Aie, Abc. Apde) = (1, 1, to solve the equation:
Ou 8u
Fu
(I, εΩ [0,
(F.18)
0t dr
0yz²
u(z, 0) = h(x,)), (z,8) Q
(F.19)
u(x, y, 1) g(a, U, (), (I, E 02 [0, T].
(F.20)
h(2, v) = 3 3sin(TT) sin( + (37x) sin(ry),
(F.21)
g(2, y, 1) ( sin (T) 2 sin(3 sin( 10x²t -10m"
(F.22)
u(x, y, t) 3 sin(Ty)e 2x sin (3 sin(ry)e 10x
(F.23)
Ω 00, 1)2, 1-
(F.24)
38

---

Burgers' 2D: We set (Aie, Abe, Apde) = (1, to solve the equation:
8em + u02u + v0," 0.05
uu + Ogy"), (1, EQ [0.7],
(F:25)
0,0 + uθ + vā,u 0.05
(Oanv + Ogyu), (2,3, EΩ x [0,T],
(F.26)
u(z,y, 0) = h1 v),
(x, E Ω
(F:27)
y,x, 0) h2( y),
(x,y) v) E Q2
(F.28)
0.25
h1 (1 3) 0.75
(F.29)
1 exp
(-4r + 4y)
0..55 32
0.25
h2 3) 0.75 +
(F.30)
+ exp
0.05 4r + 4y)
u(s, y. 1) g1 y, 1):
(z, () 0n [0,
(F.31)
u(I,y, U. 1) = 92(x1 y, 1),
(r, () E an x [0,T],
(F.32)
0.25
gi (I,,,,) 1) 0.75
(F.33)
π
1 1 exp
4x + 4y - t)
0.05
32
0.25
02(2, v.t) 0.75 +
(F.34)
1 + exp
π
( (-4r + 4y - t)
0..0
32
0.25
u(,,,,, = 0.75
(F.35)
1 + exp
π
(-4x + 4y - t)
0..5 - 32
0.25
v(2, y, t) = 0.75 +
(F:36)
1 exp
4r + 4y - t)
0.05 32
. = =0,1)² ==1.
(F.37)
Dimension-flexible Heat: We set (Aic, Ane, Apde) = 1,,1, to solve the equation:
du Au = [(x,(), (r, EQ x [0{'|oTT]..
(F.38)
Ø1
u(x, f) = Q(I,t),
(x,t) E 02 x [0,
(F.39)
u(z, 0) = h(z), ZEΩ,
(F.40)
f(z,t) w² 1) cos( Σ x, ex))xxp(t)) t).
(F.41)
d
g(2, 1) cosl
w
Σ
xz (x))exp(()).
(F.42)
7
1=1
W
h(z) = cos(
xz),
(F.4)
7
i=1
u(z,t) t) cos(
x:)cxp(-tt)
(F.44)
i=1
Ω= (-1,1)d, T = 1,
(F:45)
where
(w 20, d = 4),
(w = 30, dl =6), and
(w 40, d 8).
(F.6)
39

---

F.2 Implementation Details of Baseline Algorithms
For all PDEs, we conducted experiments by fixing the algorithms' hyperparameters to the values specified in the original baseline code. Specifically, for R3, we set number of Maxiter = 1, for Annealed R-FIPINN, we set 6 0.05, Cp 0.01, a 0.5, b 1,p 0.1, for GAS-T, we set NG = 20, No 20, X 1, for RAD, we set = 1, and for 1 we fixed the number of gradient steps at 20 and step size 0.05.
Equation
Solution
Number of N&
Number of Nic
Number of Npde
Evaluation Point
ID Burgers'
Approximated Approximated
100
100
1,000
N₂ 256, N, = 100
ID Allen-Cahn
100
100
1,000
N₂ 512, N, 200
ID KdV
Approximated
100
100
1,000
N. 200, N. 100
ID Schrödinger
Approximated
100
100
1,000
N, 256.
ID Convection
Analytic
100
100
1,000
N₂ = 512. N. 200
2D Burgers'
Analytic
100
100
1,000
N₂ 100. N. 20
2D Heat
Analytic
100
100
1,000
N₂ 100. N. 20
4D DF-heat
Analytic
100
100
1,000
N 20, N. 20
6D DF-heat
Analytic
250
250
1,000
N 8, N 8
8D DF-heat
Analytic
450
450
1,000
N₂ 5, N 5
Table 2: Summary of PDE solutions and points
40

---

G Sensitivity of Adaptive Sampling Methods
G.1
Different Loss Balance Terms (Nics Abe Apde)
Sampling methods
LAS
Kandom-K
RAD
Amecnced RIFIPIN
GAS-T
100.1.1
0.77 0:20
2.51 1:30
0.59 0.29
0.97
0.22
0.76
0.07
200,1.1
0 : 0.00
1.0 0.38
1.15 0.18
10.62
1: 19.56
0.77
0.05
ID Allen-Cahn
10.1.1
1.02 = 0.20
1.07 + 0.48
1-25-022 0:27
76.81
32.72
86.77
7.53
10.15
1.55 0.80
2006 1 0.40
2.79 0.24
65.51 18.70
96.77
0.79
7:30
100..1
1.72
0:23
0:73
6.64
0:53
2:322 1.30
8.87
0.05
3.50 1:70
200.1.1
1.82
0.42
8-87 + 1.67
6.72 - 0.42
6.00 1.19
3.06 = 1.87
9:25
0.70
3.74 2.43
4D De-hool
10.1.]
1.79 1. 0.20
7.98 37.59
24.98 37.59
16.08 1. 3.25
2.46 0.67
6.75
7.60 3.50
10.15
2.85 0.34
26.15 37.84
72.53 32
16:37 + 0
42.71 16, 79
29.86
36 25
14.47 8:07
1.85
6:34
BXL..1
25501
5.75=
200.1.1
4.18 0.58
6:21 0.82
10.63 0.93
30:51
1:39 0.61
5.80 0.72
37.20
60 DF-hest
10,1.1
4.79
0.33
6.53 + 0.81
31.81
34.30
40:47
26.37 37.05
5.40 0.34
100.00 0.00
10.15
7.33 = 1.37
70.30 + 86.30
10.0000000000
02:37 37
90:21 + 19.57
100.00 0.00
100.00 + 0.00
L.I.1
5.20 0.92
53.15 18 8 16.82
100.00 0:00
61.37 30
100.00 0.00
62:57 15.83
100,00
84.71
100000
49:22
100,00
20011.1
0.31
17:63
83.50
13.21
15.01 2.60
100.00 0.90
SD DF-heot
10.1.1
7.97
0.50
100.
0000
100,00 = 0.00
100 0.00
100.00 0.00
100,00 0.00
100.00 + 0.00
10,15
11.17
= 0.65
100 100.00- 0.00
100.00 0.00
100 00 10000000 0.00
100.00 0.00
100.00 0:000
100.00 +: 0.00
1,1,3
9.08
0.50
100.00 0.00
100.00 0:00
100.00 - 0.00
100.00 0.00
100,00 0.00
100.00 0.00
Table 3: Relative 2 error [%] performance comparison of different balance terms across varying PDEs when Npde = 1000. Bold indicates best, underline second-best.
Sampling methods
LAS (ours)
Random-R
RAD
R3
L
Number of layers
8
8
8
8
ic/bc/pde
(100/1/1)
(1/1/1)
(100/1/1)
(100/1/1)
(100/1/1)
ID Allen-Cahn
Rel L2 error [%]
0.77 1 0.20
1.07 1: 0.48
0.99 +:: 0.29
0.97 + 0.23
0.76 + 0.07
ic/bc/pde
(100/1/1)
(1/1/1)
(200/1/1)
(200/1/1)
(10/1/1)
4D DF-heat
Rel error [%]]
1.72 0,23
7.73 + 1.85
5.72 + 0.42
6.50 + 1.19
2.46 + 0.67
ic/bc/pde
(100/1/1)
(200/1/1)
(200/1/1)
(100/1/1)
(200/1/1)
6D DF-heat
Rel error [%]]
3.49
0.20
5.21 1: 0.82
10.63 + 0.93
29.91 +: 35.23
4.39 + 0.61
ic/bc/pde
(200/1/1)
(200/1/1)
(200/1/1)
(200/1/1)
(200/1/1)
8D De-heat
Rel L2 error [%]
6.92
0.31
17.63 + 1.78
83.59 +: 32.80
100.00 0.00
13.21 +: 4.46
Table 4: Relative 1,2 error [%] comparison across different PDEs, reported under the best balancing setting. Bold indicates best, underline second-bess Npde = 1000.
G.2 Different Loss Balance Terms with Smaller Number of Layers
Sampling methods
LAS (ours)
Random-R
RAD
R3
1,00
Numher of layers
100.1.1
12.46 = 3.86
28.39 + 5.90
44.14 1114:60 8.39
18.17 + 2.38
22.53 + 1.30
200
13.18 + 2.56
23.35 + 1.62
44.55 + 4.93
24.56 + 3.20
22.96 2.17
4D DF-heat
10,1,1
34.88 + 3.38
26.43 3.97
63.67 2.02
38.20 + 5.71
48.23 + 23.94
10,1,5
59.14 9.39
70.91 = 24.78
100.00
+ 0.00
68.34
+ 6.02
100.00 + 0.0
1,1,1
12.48 82.51
50.33 6.90
100.00 0.00
96.76 1 4.72
96.17 1 5.31 5
Number of layers
5
5
5
100.1,1
1.68 0.44
5.42 1.59
4.04 1 1.25
4.27 0.63
1.40 1 0.33
200.1.1
1.67 1 0.33
5.41 - 1.00
3,31 1 1.36
3.82
0.59
1.51 0.13
4D DF-heat
10,1,1
1.88 + 0.45
6.95 = 0.94
4.74 1.58
5.16 + 0.57
2.74 + 1.68
1.,1.5
3.91 1.28
8.13 1.37
11.70 - 4.32
7.46 + 2.26
4.33 + 1.93
1.1.1
2.45
0.44
6.58 1.09
9.36 1.11
5.63 0.68
3.66 + 2.16
Table 5: Relative 2 error [%] performance comparison of different balance terms in smaller number of layers when Npde 1000. Bold indicates best, underline second-best.
41

---

G.3
Different Hyperparameters of MCI Approaches on the Best ic/bc/pde Balance Terms
DF-beat 4D, io/boipde (200/1/1)
Hyper-paramekers
k=0.5, c=]
k=1, c=1
k=2
1=3,
k=4. c=1
1=5, c=1
[%]
5.72 0.42
1.73
11.78
1.19
DF-heat ie/bo/pde (200/1/1)
parame-rraasterss
k=005, c=1
k=l,
k=2, c=l
k=3, c=l
k=4. c=1
k=5, c=]
RAD
Rel 1.2 error [%]
6.39 1.14
10.63 1. 0.93
29.52 8.07
46.47 1. 27.28
66.73 1: 27.86
55.58 1. 31.73
DF-beat 8D, io/bc/pde (200/1/1)
Hyper-paranne seperrparameters
k=0.5, cml
k=1, c=1
k=2, c-l
kw3, cml
Iem4, cml
k=5, c=1
L2 error [%]
67.50 39.83
83.59 32.80
100.00 + 0.00
77.71 29.27
87.37 + 25.24
90.40 19.17
DF-hcat io/bc/pde (2001/11)
paraae-arrettrrs
Max
is3
Max
i0
L2 error
6.50
2.15 0.25
2.20 + 0.17
0.24:20.22
2.49 + 0.45
2.37 0.11
Di-heat 6D. lo/bo/pde (1001/11
Sypeerparameters
Max jel
Max in3
Max_ie5
Max Max_i=10
Max Max_ie15
Max_ i=20
Rel L2 crror [%]
29.91 + 35.23
4,98 + 0.18
5.31 + 0.86
5,91 + 0.41
6.75 + 1.02
20.55 # 28.04
DF-heat 8D, lo/bo/pde (2001/1)
Hyper-parameters
Max_i=1
Max_i=3
Max_i=5
Max_i=10
Max, Ma1_
Max_i=20
error [%]
100.00 0.00
65.03
34.35
Table 6: Effect of hyper-parameters on relative L error [%] across MCI approaches when the best ic/bc/pde balance terms of DF-heat equations and Ndde 1000 are considered.
G.4 Different Number Collocation Points
Sampling methods
LAS
Random-R
RAD
Npdc 100
49.26
100.00 0.00
100.00 0.00
100.00
88.08
23.83
Npte 500
8.62 0.42
66.77 40.75
93.37 1 13.24
100.00
49.85 40.94
Npte 1, 0000
6.92 0.31
17.63 + 1.78
83.59 = 32.80
100.00 0.00
13.21 4.46
Npte 4,
5.94 0.13
15.26 + 3.26
19.49 4.17
70.37 37.49
7.61 1.61
Nyw: 10, 000
5.90
0.39
20.09 + 0.84
16.03 1.94
35.45 32.83
6.66 1.38
Nytc 20,000 000
5.27
0.53
17.77 + 3.18
11.96 1.43
35.78 32.21
6.64 1.59
Table Effect of Npde on relative 2 error [%] across different sampling methods when ic/bc/pde balance term of DF-heat 8D equation is (200/1/1), Bold indicates best, underline second-best.
42

---

H
Compatibility Issues with Different Neural Network Architectures
We aimed experimentally evaluate the compatibility the proposed adaptive sampling technique with architectures beyond MLPs, including self-attention and modified-Mi MMLss [4]] Detailed de- scriptions of each architecture are provided in Table Furthermore, to address the spectral bias
issue commonly found in MLPs, we performed additional experiments incorporating random Fourier
blocks (FB), as detailed in [38]. The FB hyperparameters were set with Fourier feature scale of and Fourier block dimension of 64.
Table Parameter configuration for different architectures
Parameter
MLP
Self-attention
Modified MLP
Activation
Tanh
Tanh
Tanh
Embedding dimension
128
128
128
Number of layers
4
4
4
Multi-head number
N/A
4
N/A
Fully connected dimension
N/A
256
N/A
Attention dropout
N/A
0.1
N/A
Additional encoders U, V
N/A
N/A
Yes
The results are summarized Table All algorithms were evaluated using the default settings provided in the benchmark algorithm papers, and the same applies to approach, detailed in the experimental settings described the main text. Analyzing the results, we observe that the proposed LAS demonstrates high compatibility in terms of relative L error. Even in the less favorable architectures, such as MLP and self-attention, LAS achieved the recond-best performance. Notably, in scenarios incorporating FB, the proposed LAS consistently exhibited superior compatibility across all cases.
Table 9: Relative I,2 error comparison for different architectures across sampling methods. Bold indicates best, underline second-best.
Architecture
LAS
RAD
R3
L 00
Random-R
MLP
2.15 0.12
2.65 = 0.47
1.84 + 0.08
2.76 + 0.26
3.56 = 0.20
MLP + FB
0.56 = 0.14
0.69 = 0.05
0.61 0.04
0.81 0.18
0.94 0.06
Modified MLP
0.43 = 0.10
0.51
0.66 = 0.07
0.58
0,56 = 0.07
Modified MLP + FB
0.11 = 0.04
0.22 0.09
0.24 + 0.03
0.34 0.01
0.21 + 0.04
Self-attention
2.13 0.13
0.11
2.25 = 0.06
1.98 0.03
2.91 = 0.69
Self-attention FB
1.29 0.09
1.53 + 0.03
1.52 + 0.15
1.33 + 0.12
1.35 0.12
43

---

I Experimental Comparison of the Computational Complexities
key limitation of the MCMC-based sampling approach lies in its high computational complexity, primarily due to the large number of iterations required to achieve convergence of the sampling distri- bution. To demonstrate the practical feasibility of this method, we present supporting experimental results.
To evaluate computational complexities, we measured the computational costs for training deep neural networks using each algorithm. Specifically, to validate the scalability of the algorithms, conducted experiments analyze their computational requirements in terms of the number of collocation points Nps and the dimensionality of the PDE.
Based on observations, the runtime of the sampling algorithms was independent of the specific PDE. Thus, we utilized equations that allowed for a straightforward extension from 1D to 8D in dimensionality. More specifically, experimented with different sizes of collocation points (100, 1,000, 10,000, 50,000, and 100,000) for both ID, 2D Burgers' equations, and 4D to 8D De-heat equations.
As part of the detailed experimental process, we calculated the elapsed time over 1,000 epochs. The measurement was repeated 10 times using 10 different random seeds, and the mean and standard deviation were computed. For additional clarity, the clapsed time was measured excluding auxiliary operations such saving the model or storing data, focusing solely on the computations required to run the algorithms.
Hardware specification. NVIDIA RTX 4090 GPU with 24GB of memory.
Changes with NYDE. As NPDE increases, the computational cost grows for all methods. However, the growth rate varies significantly between methods. Gradient- -bsed algorithms such as LAS and 1 show a partieularly sharp increase in computational cost as Nppr grows. This is due to the iteration-intensive nature of their sampling processes. For example, with NppE = 50, 000 in 2D, the computational cost of LAS (22- 10) reaches 63.98 seconds, whereas simpler methods like Fixed Random-R R remain below 35 seconds. For NPDE = 100, 000, LAS and L run out of memory in the 2D case, highlighting their scalability limitations for very large PDE sample sizes.
Changes with dimensionality. Extending from 1D to 8D consistently increases the computational cost for all methods. While simple methods like Fixed Random-R exhibit a relatively modest increase in cost when transitioning from ID to 8D, gradient- based methods such as LAS and L show disproportionately higher computational times. For example, in the 2D case with NepE 1, 000, LAS (4. 10) takes 45.74 seconds, compared to only 16.97 seconds in ID. At NPDE = 10, 000, LAS (4. 10) takes 59.32 seconds in 2D versus 20.43 seconds in ID.
However, despite the overall computational expense of LAS for higher l values, the case of 4 = 1 demonstrates significantly lower computational costs, making it relatively practical and scalable. For example, at NNDE = 50, 000, LAS (4,- = 1) takes 63.98 seconds in 2D, which is manageable compared to the prohibitive 324.49 seconds for 4 10. Similarly, for smaller NPDE, such 1,000, LAS (42 = 1) shows competitive runtimes (e.g.. 26.62 seconds in 2D). An additional advantage of using 4 1 is that it avoids the out of memory issues observed for higher values of 4r even in large-scale scenarios such as NPpA = 100, 000 2D. Furthermore, the relative L2 error reported in our paper uses efficiency with accuracy. Thus, LAS (4.- = 1) remains a promising candidate for solving PDEs efficiently scale.
the baseline, demonstrating its effectiveness in balancing computational
Finally, while gradient- based methods like LAS exhibit high computational costs due to the overhead of gradient computation, future work optimizing the gradient operations could significantly enhance the scalability and practicality of these methods.
44

---

Table 10: Computational time (in seconds) for each method across varying problem dimensions. 0OM denotes an out-of-memory failure during execution.
Dimension
Method
100
N 1,000
= 10,000
N 50,000
= 100,000
Fixed Random-R
2.58 + 0.03
2.54 + 0.03
3.07 + 0.05
11,63 + 0.56
24.93 = 0.02
2.66 - 0.10
2.61 0.01
3.14
11.66 0.05
24.91 0.04
R3
0.01
2.90
3,20
1.97
25.34
RAD
3.40
3.46
3.90
14.92 0.03
31.81 0.09
1D
LAS (4.. 1)
3.90 0.03
+.09
5.01
21.19 0.06
46.36 0.05
LAS (L. 5)
11.59 0.03
9.70
11.90
56.86
130.95 0.08
LAS (4. 10)
15.55 + 0.05
16.97
2.11
2 0.117
000
28.77 + 1.99
29.03 2.20
37.54 0.44 0.48
197.24 0.17
448.95 0.27
Fixed
6.07 + 0.15
6,45
7.83
34.83 + 0.04
80.17 + 0.05
Random-R
6.22 + 0.21
6.16
7.85
34.87 0.04
80.26 0.04
R3
6.15
6.27
0.00
35.19 0.03
80.61 = 0.06
RAD
7.02
8.00
1.00
44.30 0.05
99.67 0.10
2D
LAS (4 1)
0.26
0.34
1.05
63.98 0.08
151.33 0.11
LAS (4 5)
2.56
2.37
3.14
179.69 0.20
0OM
LAS (4 10)
44.68 0.52
45.74
59.32
324.49 0.35
0OM
7
82.75 + 0.88
84.21 1.68
0.21 115.29
612.22 + 0.53
0OM
Fixed
11.64 + 0.50
12.01 0.06
15.06 0.13
67.23 0.03
151.58 + 0.02
Random-R
11.83 + 0.26 0.26
12.06 0.05
15.03 + 0.06
67.33 + 0.07
151.60 + 0.00
R3
11.89 + 0.34
12.23 0.04
15.43 + 0.32
67.58 + 0.03
152.03 + 0.03
RAD
15.53 + 0.70 0.70
17.13 + 0.98
19.56 + 0.05
83.29 + 0.02
OOM
4D
LAS (l2- = =)
19.17 + 0.81
19.73 0.63
25.33 0.65
124.24 + 0.16
OOM
LAS (4, 5)
42.31 2.41
42.11 #
61.20 0.23
OOM
OOM
LAS (4. 10)
68.65 + 0.70
11.24-2.0
1 0.27
OOM
OOM
L 00
133.81 1.75
139.14 2.19
193.54 + 0.17
0ОM
0OM
Fixed
17.65 + 0.64
18.64 0.61
21.10 0.12
96.50 0.08
218.68 + 0.20
Random-R
16.33 2.07
16.44 1.66
21.36 + 0.07
96.63 + 0.03
218.70 + 0.09
R3
18.38 + 0.20
1.01
21.46
96.85 0.07
219.15 + 0.01
RAD
20.69 0.65
20.05
25.75
119.51 0.11
OOM
6D
LAS (4_ = 1)
25.97 1.51
1.98
34.54
179.94 + 0.32
OOM
LAS (4. 5)
57.66 1.18
2.76
87.11
OOM
0OM
LAS (4. 10)
102.97
3.22
103.50 5.22
153.10 0.25
0OM
0OM
7
183.39
7.13
6.34 195.98
275.06 0.19
00M
00M
Fixed
20.19 + 1.93
21.93 0.38
26.55 0.25
127.65
0.16
286.65 0.16
Random-R
22.08 + 0.13
21.84 0.20
26.95 0.18
127.68
0.14
286.61 + 0.01
R3
20.99 + 0.04
22.51 0.31
27.15 + 0.09
127.85
0.14
287.26 0.05
RAD
26.48 0.88
26.15 1.01
32.15
157.70
0.11
OOM
8D
LAS (4.- = 1)
0.99
33.75 +
44.89
1.89-0.0 0.15
OOM
OOM
LAS (4.- 5)
74.12
78.77 + 1.76
114.37 0.27
OOM
0OM
LAS (l. 10)
128.65 5.59
134.49 4.29
201.64 0.26
OOM
0OM
D
229.60 1.97
240.58 + 5.01
357.80 0.14
0OM
0OM
45

---

J Additional Results with Varying LAS Hyperparameters
J.1 Practical Implementation
In general, achieving successful Langevin sampling requires careful selection of hyperparameters (Langevin step size т, number of Langevin iterations I,, etc.). While running Langevin dynamics for many iterations with small step size might allow sampling to be proportional to the actual residual landscape, can significantly slow down the training speed of the PINN model in practical applications. Thus, we considered the following concepts when setting the hyperparameters, which are crucial for practical utilization of Langevin dynamics in PINN training.
Adjusting Langevin step size and iteration. To increase the computational efficiency, we adopted strategy of increasing the step size and reducing the number of Langevin iterations 4, We anticipated that the temporal variation of the PINN model, fe. would exhibit smooth behavior when utilizing adaptive sampling strategies with large network size. Consequently, even with fewer Langevin iterations, minimal changes in the loss landscape suggest that the sample trajectory would resemble that of fixed landscape.
Normalizing the gradient size. Since LAS leverages gradient information from the residual land- scape, the step size T needs to be set smaller for stiff PDEs, i.e., the step size is dependent the PDE. Additionally, empirical observations indicate that the gradient of the residual landscape exhibits substantial variations at the beginning of training. In contrast, towards the end of training, the residual landscape is characterized by relatively small gradients. This discrepancy restricts the movement of sample points, thereby making challenging secure reliable quality. To address these challenges, we normalized the magnitude of all residual gradients at each iteration relative to the largest residual gradient.
J.2 Variation of 3 and 44. with Fixed 0.002
1. Instability of performance for small 3 values as the layer increases: As the layer depth increases, small 8 values lead to unstable performance. For instance, in layer 4, small 3 0.001 results in a relatively stable error value of 1.18 + 0.23 at 4 1, whereas in layer 10, the error rises significantly 8.20 = 15.63. This pattern suggests that small values hinder performance stability in deeper layers.
2. Increased instability with higher I, values: Generally, the performance deteriorates as H. increases, particularly for small / values. For example, in layer 6 with 3 0.001, the relative error increases from 0.58 + 0.08 at I.. to 22.73 + 22.08 I, = 20. This indicates that excessive Langevin iterations could lead to performance instability, especially when the concentration parameter is low.
Overall, the impact of and It on both relative T,2 error and variance limited, indicating robustness of the method across range of parameters. This robustness simplifies Langevin hyperparameter tuning, making the approach more practical for real-world applications.
J.3 Variation of and 4, with Fixed B
1. Stability across 7 values: The relative I,² error exhibits minor fluctuations across different values of for each layer suggesting a negligible dependency т.
2. Limited impact of l,= While increasing 4, slightly reduces the variance of L error in some cases, the effect is not consistent across layers showing only marginal improvement.
As a result, the main takeaway is that small values of are prone to instability as the layer depth increases and 4, becomes large. Conversely, higher 3 values can ensure stable performance, even with varying 4, values. In shallow layers, however, lower 3 values can be beneficial, providing more precise error at lower l values.
46

---

Layer and concentration parameter 3
Langevin iteration It.
Layer
3
0.001
1.18 = 0.23
1.63 + 0.40
1.66 = 0.20
6.34 + 2.69
0.05
1.53
1.4+
1.27 0.21
1.25 0.16
Layer 4
0.1
1.08 0.37
0.930.00
1.23 0.23
1.06 0.13
0.3
2.50 0.28
2.30 =
2.54 0.35
2.51 + 0.36
0.4
2.80 0.60
3.14
3.14 0.49
2.90 0.44
0.001
0.58 0.08
10.89 19.65
10.10
17.80
22.73 22.08
0.05
0.62 0.11
0.64 0.12
0.65 0.12
0.67 0.15
Layer 6
0.1
0.74 0.12
0.58 + 0.08
0.64 0.07
0.58 - 0.03
0.3
1.05 0.27
1.15 0.21
0.11
1.33 0.13
0.4
1.47 = 0.10
1.49 1 0.08
1.48 = 0.22
1.58 + 0.13
0.001
0.65 0.14
+.25
10.66 19.50
22.14 24.51
0.05
0.88 = 0.34
0.46 0.05
0.59 = 0.10
1.44 + 0.81
Layer 8
0.1
0.59
0.58
0.35
0.3
1.06 0.16
1.25
1.08 0.24
0.92 + 0.11
0.4
1.12 = 0.25
1.09 1 0.19
1.14 = 0.17
1.12 1 0.27
0.001
8.20 15.63
23.14 + 20.74
40.14 19.78 17.21
41.79 + 20.50
0.05
0.46 0.17
17.33 + 21.03
21.02
30.08
24.07
Layer
0.1
0.63 0.20
0.73
1 0.23
0.70
0.22
10.23 18.99
0.3
0.94 0.18
0.98
+ 0.15
0.83 0.25
1.02 + 0.30
0.4
0.94
0.96
0.37
1.12
Layer and Langevin step size
Langevin iteration lT.
Layer
5
10
20
0.0001
1.98 + 0.35
2.12
+ 0.14
1.66 + 0.06
1.74
- 0.19
0.0005
1.81 0.13
1.74
1 0.40
1.77
1 0.34
1.74
0.35
Layer
0,001
2.35 - 0.37
1.85
+.25
0.25
1.98
1.68 = 0.24
0.005
1.96 1 0.15
2.08 + 0.188
1.22
1.71 0.17
0.01
1.93 - 0.28
1.92 0.09
1.48
1.79
0.09
0.0001
0.96
0.96
1.00
0.14
0.0005
0.92 0.13
0.69
1.14
1.17
Layer
0,001
0.85 =
0.18
0.00
0.96
0.005
0.85 0.07
0.22
0.97
0.78 0.37
0.01
0.90 : 0.04
0.90
L.11
0.93 0.06
0.0001
0.62 + 0.04
0.79
0.92 0.100
0.82 + 0.05
0.0005
0.81 = 0.02
0.85
0.78 0.07
0.75 = 0.06
Layer
0.001
0.91 - 0.13
0.82
0.71 0.01
0.77 + 0.10
0.005
0.85 = 0.07
0.83
0.64 1.00
0.63 - 0.08
0.01
0.64 1 0.14
0.77 0.17
0.88 1.88
0.81 = 0.09
0.0001
0.41 + 0.18
0.75
0.08
0.82 0.82 0.02
0.68 + 0.12
0.0005
0.82 + 0.13
0.64 1 0.04
0.68 0.07
0.82 + 0.06
Layer 10
0,001
0.91 - 0.21
0.74
0.87 0.225
0.59 - 0.05
0.005
0.64 + 0.04
0.53 + 0.32
0.67 + 0.05
0.64 + 0.12
0.01
0.56 0.29
0.65 1
0.71 0.14
0.61 0.09
Table 11: Relative 2 error for varying 3 and and 4 across different layers of the Allen-Cahn equation.
47


Reviews:
[{"review_id": "aXFDO8tgtZ", "reviewer": "Reviewer_3GbE", "summary": "This paper proposes to use Langevin dynamics to update collocation points when training PINNs. The authors analyzed the gradient behaviors for training with high residual points and training with random samples, and showed that there could be big difference between these two approaches regarding convergence and stability. The authors consequently proposed to update the collocation points locally via Langevin dynamics, which may favor flat regions and enjoy better training stability. The proposed method is tested on several 1D PDEs as well as high-dimension PDEs.", "questions": "- What is $a(\\theta)$ in equation 3.2?", "limitations": "- Limitation is discussed adequately.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength\n\n- The paper is well written and easy to follow, despite having certain theoretical depth.\n\n- The paper empirically benchmarked several collocation sampling methods, which are tested on both simple and high-dimensional PDEs.\n\nWeakness\n\n- Although the overall experimental results show favorable performance for the proposed method. I feel the gain in the performance seems not very strong. For example, as shown in Table 1, for some PDEs, random sampling still yields the best results. The results for high-dimensional PDEs are significantly better for the proposed method. However, the error seems to increase with the dimensionality, and the benchmarked existing methods do not really work and they may not be tailored for such problems.\n\n- Nevertheless, although I am not expert on tuning PINNs, I think the proposed Langevin dynamics sampling is an interesting idea and its motivation is well analyzed.", "was_revised": true, "final_justification": "The authors clarified some derivation and experimental evaluations. For high-dimensional problems, although the authors explained that they have tried reasonable hyperparameters, I believe getting improved results for baseline methods would further improve the presentation. Overall I will keep my rating as accept.\n\nUpdate:\n\nThe updated presentation about the high-dimensional results addressed my concern about baseline results. I remain supportive for acceptance.", "rebuttal": "**Dear Reviewer 3GbE,**\n\nWe sincerely appreciate your interest in our work and your positive evaluation.\n\nAs you correctly pointed out, we observed that **Random-R** (uniform random sampling at each iteration) generally performs well when using a sufficiently large number of collocation points (in our case, $N_{PDE}$=1000) for relatively smooth PDEs, such as Burgers’, Schrödinger, and KdV equations. However, for **stiff PDEs** such as the Allen–Cahn equation, all adaptive sampling methods consistently outperformed Random-R. Moreover, we found that existing benchmark adaptive sampling methods struggled to generalize across different PDE types—smooth, stiff, and high-dimensional (although we did not exhaustively search the entire hyperparameter space, we explored reasonable hyperparameter configurations as reported in Appendix H, pp. 33–34). In contrast, our proposed **LAS** method demonstrated consistently generalized performance across different PDEs by simultaneously prioritizing both high-residual and relatively flat regions, which was particularly effective in **high-dimensional cases**.\n\n---\n\n## Response to Q1: What is a($\\theta$) in Equation (3.2)?\n\nWe thank the reviewer for asking for clarification regarding the definition of $a(\\theta)$.\n\nIn our formulation, $a(\\theta)$ is the coefficient vector that arises from linearizing the residual function $R_\\theta(\\mathbf{x})$ with respect to a fixed basis of differential features, denoted $\\phi(\\mathbf{x})$. Specifically, when the residual is expressed as\n$$\nR_\\theta(\\mathbf{x}) = (\\partial_t f_{\\theta} + N_x[f_\\theta])(\\mathbf{x}) \\triangleq g(f_\\theta)(\\mathbf{x}) \\approx a(\\theta)^\\top \\phi(\\mathbf{x}),\n$$\nthe vector $a(\\theta)$ contains the weights applied to each feature in this expansion.\n\nThis structure is motivated by a first-order linearization of the operator $g$ in function space. If the residual operator admits a smooth form such as\n$$\ng(f) = G(f, \\partial_t f, \\partial_x f, \\partial_{tt} f, \\partial_{tx} f, \\dots),\n$$\nthen applying the Fréchet derivative yields:\n$$\ng(f_{\\theta + \\Delta\\theta}) \\approx g(f_\\theta) + D_g(f_\\theta)(f_{\\theta + \\Delta\\theta} - f_\\theta)\n$$\nUsing a Taylor expansion in terms of the partial derivatives of $f_\\theta$, the residual can be approximated as:\n$$\nR_\\theta(\\mathbf{x}) \\approx \\sum_k \\frac{\\partial G}{\\partial (\\partial_k f_\\theta)} \\Delta(\\partial_k f_\\theta) \\equiv a(\\theta)^\\top \\phi(\\mathbf{x})\n$$\nHence, $a(\\theta)$ is not an independently learned parameter, but rather a structured function of the network $f_\\theta$ and its derivatives with respect to both input and parameters.\n\nTo ground this in concrete examples:\n\n- **Convection equation**: $g(f) = \\partial_t f + \\mu \\partial_x f$ → $a(\\theta) = (0, 1, \\mu, 0, 0, \\dots)$ is constant\n- **Burgers’ equation**: $g(f) = \\partial_t f + f \\partial_x f - \\frac{0.01}{\\pi} \\partial_{xx} f$ → $a(\\theta) = (\\partial_x f_\\theta, 1, f_\\theta, 0, \\dots, \\frac{0.01}{\\pi})$ is input- and $\\theta$-dependent\n\nA more detailed description is provided in Appendix A.1. We acknowledge that the explanation of this coefficient vector was missing in the current manuscript, and we will include this clarification in the revised version.\n\nWe hope these responses address your concerns and that our planned revisions appropriately reflect your valuable suggestions. Please feel free to share any further comments or questions. Thank you again for your thoughtful comments."}, {"review_id": "UPJgr7NzUl", "reviewer": "Reviewer_DvNh", "summary": "In this paper, a new adaptive sampling method called LAS is proposed to enhance the accuracy of Physics-Informed Neural Networks (PINNs). LAS addresses the issue of neglecting points with medium or low residuals, which can compromise stability. The experiments demonstrate that the proposed method outperforms several existing approaches in terms of relative error and stability across a variety of environments, including high-dimensional partial differential equations (PDEs), where Monte Carlo integration-based methods typically struggle with instability.", "questions": "See weakness.", "limitations": "More numerical comparision is needed,  see weakness.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "engths: The paper is well-written, and the motivation is clearly articulated.\n\nWeaknesses: The numerical analysis is insufficient to demonstrate the value of this paper, particularly due to the lack of comparisons with closely related works such as \n\n[1] Failure-informed adaptive sampling for PINNs\n[2] Failure-informed adaptive sampling for PINNs, Part II: Combining with Re-sampling and Subset Simulation\n[3] A Gaussian Mixture Distribution-Based Adaptive Sampling Method for Physics-Informed Neural Networks", "was_revised": true, "final_justification": "Thanks the authors' reply. I will increase my score.", "rebuttal": "**Dear reviewer DvNh**\n\nWe sincerely appreciate your comments\n\n---\n\n## Response to Weakness: Lack of comparisons with closely related works\n\nAs you correctly pointed out, we overlooked some closely related studies [1–3] that are highly relevant to our work. We will cite these works in _Related Work_ and explicitly discuss how they relate to both existing methods and our proposed approach.\n\n**[1] Failure-informed PINN (FI-PINNs):**  \nIn FI-PINNs, “failure” refers to collocation points where the PDE residual error exceeds a predefined threshold ​. A truncated Gaussian model for **Self-adaptive Importance Sampling (SAIS)** is employed to preferentially sample collocation points from these high-residual (failure) regions. However, this method suffers from two key limitations:\n\n1. The number of collocation points increases over time as the failure regions are updated, which leads to additional computational cost.\n\n2. It struggles to handle **multi-modal failure regions** effectively.\n\n**[2] Annealed R-FIPINN (Extensions of FI-PINNs):**\n\nTo address these limitations, two important extensions were introduced:\n\n1. **Cosine-annealing-based resampling technique**, which maintains a fixed training size while gradually adapting to failure regions.\n\n2. **Subset Simulation (SS)** based on Markov Chain Monte Carlo (MCMC) using a modified Metropolis–Hastings (MMA) sampler, which replaces the truncated Gaussian model and better captures multi-modal failure regions.\n\n**[3] Gaussian mixture distribution-based adaptive sampling (GAS):**  \nGAS improves the representation of multi-modal residuals by selecting the top NG​ high-residual validation points as GMM means and using the reciprocal of the residual gradient as the diagonal of the covariance matrix. However, similar to FI-PINNs, GAS continuously increases the total number of collocation points over time—despite maintaining the ratio between collocation and boundary points—which leads to progressively higher computational costs.\n\nWe compared our **LAS** method against **Annealed R-FIPINN** and **GAS-T** under the same experimental conditions. To ensure fairness, we maintained the same $N_{BC}$​, $N_{IC}$​, $N_{PDE}$​, learning rate, and total epochs across all methods. Importantly, for GAS-T, we fixed the total number of collocation points by replacing existing collocation points at each update rather than aggregating them. **Bold** indicates the best, while ***Bold + Italic*** denote the second-best alternative.\n\n\n| PDE | Balance terms (IC/BC/PDE) | LAS | Random-R | RAD | R3 | $L^{\\infty}$ | Annealed R-FIPINN | GAS-T |\n|-----|----------------------------|-----|-----------|-----|-----|--------------|---------------------|--------|\n| **DF-Heat 4D** | (100,1,1) | **1.72 ± 0.23** | 9.73 ± 0.44 | 6.64 ± 2.72 | 9.53 ± 5.27 | ***2.92 ± 1.95*** | 8.87 ± 0.65 | 3.80 ± 1.70 |\n|  | (200,1,1) | **1.82 ± 0.42** | 8.87 ± 1.67 | 5.72 ± 0.42 | 6.50 ± 1.19 | 4.06 ± 1.87 | 9.25 ± 0.70 | ***3.74 ± 2.43*** |\n|  | (10,1,1) | **1.79 ± 0.20** | 7.98 ± 37.59 | 24.98 ± 37.59 | 16.08 ± 3.25 | ***2.46 ± 0.67*** | 6.75 ± 1.97 | 7.60 ± 3.30 |\n|  | (10,1,5) | **2.85 ± 0.34** | 25.18 ± 37.84 | 72.53 ± 32.73 | 15.57 ± 9.22 | 42.71 ± 46.79 | 29.86 ± 36.25 | ***14.47 ± 8.07*** |\n|  | (1,1,1) | **1.86 ± 0.20** | ***7.73 ± 1.85*** | 12.31 ± 1.24 | 13.56 ± 5.24 | 10.95 ± 14.41 | 26.36 ± 36.89 | 16.59 ± 11.01 |\n| **DF-Heat 6D** | (100,1,1) | **3.49 ± 0.20** | 6.14 ± 0.42 | 13.84 ± 2.08 | 29.91 ± 35.23 | ***4.46 ± 0.42*** | 5.75 ± 1.05 | 100.00 ± 0.00 |\n|  | (200,1,1) | **4.18 ± 0.58** | 5.21 ± 0.82 | 10.63 ± 0.93 | 30.51 ± 34.90 | ***4.39 ± 0.61*** | 5.80 ± 0.72 | 37.20 ± 36.19 |\n|  | (10,1,1) | **4.79 ± 0.33** | 5.53 ± 0.81 | 31.81 ± 34.30 | 49.47 ± 41.28 | 26.37 ± 37.05 | ***5.49 ± 0.34*** | 100.00 ± 0.00 |\n|  | (10,1,5) | **7.33 ± 1.37** | 70.30 ± 36.59 | 100.00 ± 0.00 | ***69.37 ± 37.51*** | 90.21 ± 19.57 | 100.00 ± 0.00 | 100.00 ± 0.00 |\n|  | (1,1,1) | **5.20 ± 0.92** | *53.18 ± 46.82* | 100.00 ± 0.00 | 61.37 ± 39.22 | 100.00 ± 0.00 | 62.57 ± 45.83 | 100.00 ± 0.00 |\n| **DF-Heat 8D** | (100,1,1) | **7.72 ± 0.59** | 52.41 ± 39.33 | 84.71 ± 30.56 | 100.00 ± 0.00 | 100.00 ± 0.00 | ***49.22 ± 41.61*** | 100.00 ± 0.00 |\n|  | (200,1,1) | **6.92 ± 0.31** | 17.63 ± 1.78 | 83.59 ± 32.80 | 100.00 ± 0.00 | ***13.21 ± 4.46*** | 15.01 ± 2.69 | 100.00 ± 0.00 |\n|  | (10,1,1) | **7.97 ± 0.30** | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 |\n|  | (10,1,5) | **11.17 ± 0.65** | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 |\n|  | (1,1,1) | **9.08 ± 0.50** | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 | 100.00 ± 0.00 |\n\n---\n\n| Method | Sampling Hyperparameters |\n|--------|----------------------------|\n| **LAS** | $\\tau$=0.01, $\\beta$=0.2, $l_L$=1 |\n| **Random-R** | NA |\n| **RAD** | k=1, c=1 |\n| **R3** | MaxIter=1 |\n| **$L^{\\infty}$** | $\\eta$=0.05, K=20 |\n| **Annealed R-FIPINN** | $\\epsilon_r$=0.05, $\\epsilon_p$=0.01, a=0.5, b=1, p=0.1 |\n| **GAS-T** | $N_G$=20, $N_a$=20, λ=1 |\n\n---\n\n**Training Hyperparameters**: \n\nEpochs=10,000, $N_{PDE}$=1,000, $N_{IC}$=100/250/450 (4D/6D/8D), $N_{BC}$=100/250/450 (4D/6D/8D), Learning Rate=1e-3, # of Hidden Layers=8, # of nodes = 128, Optimizer: Adam\n\n---\nWe hope these responses have addressed your concerns and that our planned revisions clearly reflect your valuable suggestions. Please don’t hesitate to share any further feedback or questions. \n\nThank you once again for your thoughtful and constructive comments."}, {"review_id": "tVTFlE91PI", "reviewer": "Reviewer_b6R2", "summary": "This paper introduces a Langevin dynamics-based Adaptive\n Sampling (LAS) for the solution of PDEs under the PINN framework. This work is part of a wider literature on adaptive sampling strategies for PINNs colloration point residuals, roughly divided into distribution based sampling and high-residual based sampling. The authors focus on the latter, proceeding by providing a systematical analysis of the theoretical properties of high residual sampling and their connection to learning instability. They identify sharp gradients as one the min source of instability and suggest the use of Langevin dynamics to update the residuals. This is a change from previous approaches based on sampling independent residuals from a target distribution at every iteration. The use of noise in the LAS step has a smoothing effect which the authors argue leads to flatter local optima and improved algorithm stability.", "questions": "As I mentioned, the paper looks good. However, in my opinion it critically lacks two things:\n1. A discussion on the computational requirements. This should involve both comments on the hardware used, and viability of the algorithm on other kind of hardware (e.g. smaller GPUs, CPUs). Moreover, it should mention the computational complexity (empirical or theoretical) and how the algorithm performs against the baselines. Part of this is present in the appendix but, in my opinion, it is very important to give visibility to this issues in the main body with a high-level discussion. Details can be provided in the appendix.\n2. The presentation of the LAS performance compared to the baseline (or its own performance with 8 and 10 layers) for a smaller number of layers (1, and 5). This provides a fairer comparison in terms of compute, and provides evidence to practitioners on whether this approach can be used in more resource costraiend environments.", "limitations": "Limitations were not mentioned. Crucially, the authors failed to comment on the computational cost-accuracy tradeoff.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n- The paper is very well written, and very clear. Morevoer, while the main body is dense with results, both theoretical and empirical analysis, this does not affect the flow of the exposition.\n- The authors have been methodic and thorough in their analysis: they develop a theoretical analysis under some assumptions, which are validated eitehr with references to prior work, or in the appendix. The main theorems are also assessed empirical to show additional evidence (e.g. Figure 3).\n- The experimental evaluation of their model is similarly precise. Eight different PDE systems are assessed which comprise increasing complexities and a wide range of behaviors. The performance is systematically compared with four competing approaches.\n- The experimental results prove the validity of the approach, point out the theoretical features of the algorithm described in the first half of the paper, and highlight the pitfalls of the competing algorithms.\n- Extensive ablation studies are reported, and algorithm runtimes are provided.\n- Finally, clear implementation details are provided, which allow other researchers to replicate and build on their work, and practictioners to use their methodology.\n\nOverall, I am satisfied with the quality of this paper.\n\n**Weaknesses**\n- The authors fail to provide a clear paragraph covering the limitations of their proposed approach.\n- While algorithm runtimes are reported in the appendix, these are not commented in the text. However, LAS is found to be more computationally expensive than the other approaches, except than L^\\infty. Generally, LAS is 1.2 to 6 times slower, using 1 to 10 number of layers. However, results in the main body are only reported for 8 and 10 layers, which are significantly more expensive and more prone to out-of-memory (OOM) issues for high problem dimensions, as reported in Table 8. Note that OOM is recorded using a RTX 4090 GPU with 24GB memory, so for practitioners with more limited hardware, they may encounter memory bottlenecks much sooner.", "was_revised": true, "final_justification": "For my considerations above, I find the paper merits publication. The authors engaged in a constructive rebuttal process with both me and the other reviews, and they address all the minor concerns that I pointed out above.", "rebuttal": "**Dear reviewer b6R2**\n\nWe sincerely appreciate insightful your comments and positive evaluation of our work\n\n---\n\n## Response to Weaknesses 1: Lack of a clear paragraph discussing the limitations of the proposed approach\n\nThank you for pointing this out. First, we would like to clarify that $l_{L}$​, which denotes the number of Langevin iterations used to update collocation points, was introduced in Algorithm 1 and Section 4.1. However, we acknowledge that this was not clearly mentioned in Appendix I, and we apologize for the confusion.\n\nSpecifically, we analyzed the performance of LAS across different hyperparameter settings ($l_{L}$​, $\\tau$, $\\beta$) in Appendix D (pp. 21–22), with respect to different collocation point configurations in Appendix H.3 (p. 34), and computational costs in Appendix I (p. 36). Nevertheless, as you rightly noted, the main text lacks a proper discussion on the limitations and guidance for practitioners—particularly how to balance computational cost and hyperparameter tuning across layers and collocation points.\n\nWe will revise the main text to include this discussion, focusing on the trade-off between performance and computational cost as the PDE dimensionality and number of collocation points increase. We will also provide guidance for practitioners on effectively applying LAS in resource-constrained environments.\n\n---\n\n## Response to Weakness 2: Algorithm runtimes are only reported in the appendix\n\nAs you correctly observed, LAS becomes significantly slower with higher values of $l_{L}$ ​(e.g., from 1 to 10), except for the $L^{\\infty}$. However, we also found that setting $l_{L} = 1$ is already sufficient to outperform competing methods in most cases, thereby keeping computational complexity comparable to others. All results in our main experiments were obtained with $l_{L} = 1$, and further performance comparisons for different $l_{L}$ values are provided in Appendix D (pp. 21–22). Nonetheless, as the number of collocation points and the problem dimensionality increase, LAS encounters memory bottlenecks earlier than methods like R3 and RAD. We will explicitly mention this limitation in the main text and recommend appropriate settings considering available memory. Of course, the number of layers affects the computation complexity, but difference of computational cost among sampling methods mainly comes from collocation points, PDE dimension, and repentance for sampling collocation points. However, we acknowledge that this was not clearly mentioned in Appendix I, and we apologize again for the confusion.\n\n---\n\n## Response to Question 1: Computational requirements\n\nWe agree with your suggestion and will include a visual representation of the computational complexity and memory requirements for different GPU models. This will be shown across varying collocation point configurations and PDE dimensionalities.\n\n---\n\n## Response to Question 2: Performance for smaller numbers of layers (1 and 5)\n\nThank you for this helpful suggestion. Performance in low-layer scenarios is particularly relevant in constrained environments.\n\nWhile fewer layers tend to experience less training instability, the comparative performance remains informative.\n\nIn Figure 4(a), LAS achieved the second-best performance in the 4-layer setting and the best performance in the 10-layer case. Additionally, we include comparisons for the 4D heat equation with different balancing terms and with 1 and 5 layers. **Bold** indicates the best, while ***Bold + Italic*** denote the second-best alternative.\n\n| PDE               | Balance terms (IC/BC/PDE) | LAS              | Random-R         | RAD              | R3               | $L^{\\infty}$      |\n|------------------|----------------------------|------------------|------------------|------------------|------------------|----------------------|\n| **DF-Heat 4D (1 layer)** | (100,1,1)              | **12.49 ± 3.86**  | 28.39 ± 5.90     | 44.14 ± 8.39     | ***18.17 ± 2.38***   | 22.53 ± 1.30         |\n|                  | (200,1,1)                  | **13.18 ± 2.56**  | ***23.35 ± 1.62***   | 44.55 ± 4.93     | 24.56 ± 3.20     | 22.96 ± 2.17         |\n|                  | (10,1,1)                   | ***34.88 ± 3.38***  | **26.43 ± 3.97**   | 63.67 ± 2.02     | 38.20 ± 5.71     | 48.23 ± 23.94        |\n|                  | (10,1,5)                   | **59.14 ± 9.39**  | ***70.91 ± 24.78***  | 100.00 ± 0.00    | 68.34 ± 6.02     | 100.00 ± 0.00        |\n|                  | (1,1,1)                    | ***82.51 ± 12.48*** | **50.33 ± 6.90**   | 100.00 ± 0.00    | 96.76 ± 4.72     | 96.17 ± 5.31         |\n| **DF-Heat 4D (5 layers)** | (100,1,1)           | ***1.68 ± 0.44***   | 5.42 ± 1.59      | 4.04 ± 1.25      | 4.27 ± 0.63      | **1.40 ± 0.33**        |\n|                  | (200,1,1)                  | ***1.67 ± 0.33***   | 5.41 ± 1.00      | 3.31 ± 1.36      | 3.82 ± 0.59      | **1.51 ± 0.13**        |\n|                  | (10,1,1)                   | **1.88 ± 0.45**   | 6.95 ± 0.94      | 4.74 ± 1.58      | 5.16 ± 0.57      | ***2.74 ± 1.68***        |\n|                  | (10,1,5)                   | **3.91 ± 1.28**   | 8.13 ± 1.37      | 11.70 ± 4.32     | 7.46 ± 2.26      | ***4.33 ± 1.93***        |\n|                  | (1,1,1)                    | **2.45 ± 0.44**   | 6.58 ± 1.09      | 9.36 ± 1.11      | 5.63 ± 0.68      | ***3.66 ± 2.16***        |\n\n---\n\n| Method | Sampling Hyperparameters |\n|--------|----------------------------|\n| **LAS** | $\\tau$=0.01, $\\beta$=0.2, $l_L$=1 |\n| **Random-R** | NA |\n| **RAD** | k=1, c=1 |\n| **R3** | MaxIter=1 |\n| **$L^{\\infty}$** | $\\eta$=0.05, K=20 |\n\n---\n\n**Training Hyperparameters**: \n\nEpochs=10,000, $N_{PDE}$=1,000, $N_{IC}$=100, $N_{BC}$=100, Learning Rate=1e-3, # of Hidden Layers=1, 5, # of nodes = 128, Optimizer: Adam\n\n---\n\nWe hope these responses address your concerns and that our planned revisions appropriately reflect your valuable suggestions. Please feel free to share any further comments or questions.\n\nThank you again for your thoughtful review."}, {"review_id": "y6DwpE0vP5", "reviewer": "Reviewer_Jy8K", "summary": "This paper introduces a mechanism for training PINNs using Langevin sampling to learn the distribution of residuals. The authors provide theoretical guarantees on the convergence properties of the method (i.e. distributional convergence and flatness). A comprehensive set of experiments were conducted, and LAS shows state-of-the-art results in six out of ten benchmarks, including the notably challenging high-dimensional PDEs.\n\nAlthough I think this work is interesting from an experimental perspective and well-informed in the existing literature, I find the theoretical aspects of this paper concerning. I am willing to increase my score if the authors provide improvements in the proofs in a rebuttal.", "questions": "* Can the authors clarify the inconsistency between propositions 3.1 and 3.2, and assumption 3.1? There seems to be a contradiction.\n* Every theorem/proposition, by convention, should list all the underlying assumptions explicitly instead of making the reader guess which assumptions hold.\n* I think the authors are missing a proof when bridging assumptions 3.2 and 3.3. Heavy-tailed distributions, by definition, prevent tight concentration, so these two cannot simultaneously hold. What are the conditions that the weights or $\\beta$ must satisfy for both assumptions to hold?\n* Theorem 3.1 also seems unclear. If $\\max\\|\\phi\\|\\sim D^{1/\\alpha}$, then $(\\max\\|\\phi\\|)^2\\sim D^{2/\\alpha}\\sim\\lambda_{max}(\\Sigma)$, so I don't see why you argue that $(\\max\\|\\phi\\|)^2 >> \\lambda_{max}(\\Sigma)$", "limitations": "A practical concern is that the benchmarking methods require fewer neural network evaluations compared to the proposed method, which involves solving an SDE with neural network evaluations at every time step. But I appreciate that the authors addressed this in the paper.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n* I believe this paper explores some interesting directions that improve the training of PINNs. The search for collocation points for PINNs is relevant, especially in high dimensions.\n* The assumptions (3.1, 3.2, 3.3) are well-formulated and fairly weak so that they don't compromise the generalizability of the theoretical results.\n\n**Weaknesses**\n\n* This paper, although theoretically correct, is inconsistent with continuous-time and discrete-time formulations of Langevin dynamics. For example, Equation (4.1) is the Euler-Maruyama discretization of the Langevin equation $dX_t=\\nabla_{x}|R_{\\theta}(x_n^l)|^2+\\beta dB_t$. I should expect to see both formulations in the paper.\n* Although Theorem 4.1 is correct, its proof is wrong. Moreover, in a continuous-time FP equation, I should not expect to see any discretization in time ($\\tau$) in the proof. In the rebuttal, I want to see the authors use the continuous-time SDE formulation in the proof. I recommend the following sources [1] [2].\n* There are serious inconsistencies in the structure of the residual. In the proof of Proposition 3.1, you write that $R_\\theta(x)=\\theta^T\\phi(x)$, which is different from Equation (3.3) in Assumption 3.1. The proof of Proposition 3.2 has the same issue.\n* Theorem 3.1 also seems wrong (see questions)\n\n[1] Jordan, R., D. Kinderlehrer, and F. Otto (1998). “The Variational Formulation of the Fokker–Planck Equation”. SIAM Journal on Mathematical Analysis, vol. 29, no. 1, pp. 1–17.\n\n[2] Ambrosio, L., N. Gigli, and G. Savaré (2008). “Gradient Flows in Metric Spaces and in the Space of Probability Measures”. 2nd. Springer.", "was_revised": true, "final_justification": "Based on the rebuttal, the authors proposed improvements in the mathematical foundation of their work and posted proofs that I verified are correct. The consistency in their assumptions was justified in the rebuttal, and I expect to see this change in the paper.", "rebuttal": "## Response to Q1: Discrete vs. Continuous-Time Langevin Dynamics\n\n**Summary:  \nWe acknowledge the importance of presenting both discrete- and continuous-time Langevin dynamics. We will explicitly include both formulations and their theoretical connection in the revised version.**\n\nWe thank the reviewer for pointing out the inconsistency between the discrete-time and continuous-time Langevin formulations. We agree with the reviewer that both discrete-time and continuous-time formulations should be presented for consistency. We will include the full derivation in the appendix and summarize the key steps in the main text for clarity.\n\nOur experiments implement the following discrete-time Langevin dynamics via the Euler–Maruyama scheme. For $l \\in \\mathbb{N}$:\n\n$$\\mathbf{x_{\\textit{l+1}}} = \\mathbf{x_{\\textit{l}}} + \\frac{\\tau}{2} \\nabla_\\mathbf{x} |R_\\theta(\\mathbf{x_{\\textit{l}}})|^2 + \\beta \\sqrt{\\tau} z_l,\\quad z_l \\sim \\mathcal{N}(0, I).$$\n\nWe clarify that $\\mathbf{x} = (t,x)$ (time, space) denotes the PDE domain, and we use $s$ to index continuous time for clarity.\n\nAs the step size $\\tau \\to 0$, under mild regularity conditions such as Lipschitz continuity of the gradient, the above discrete-time process converges in distribution to the following continuous-time stochastic differential equation (SDE):\n\n$$\nd\\mathbf{x_{\\textit{s}}} = \\frac{1}{2} \\nabla_{\\mathbf{x}} |R_\\theta(\\mathbf{x_{\\textit{s}}})|^2ds + \\beta dB_s.\n$$\n\n\n\nwhere $B_s$ denotes standard Brownian motion. This standard continuous-time limit justifies our use of the Fokker–Planck equation. We will nonetheless include both formulations for clarity and completeness.\n\nThe corresponding Fokker–Planck equation describing the time evolution of the probability density $p_s(x)$ is:\n\n$$\n\\frac{\\partial p_s(\\mathbf{x})}{\\partial s} = -\\nabla_\\mathbf{x} \\cdot \\left( \\frac{1}{2} \\nabla_\\mathbf{x} |R_\\theta(\\mathbf{x})|^2 p_s(\\mathbf{x}) \\right) + \\frac{\\beta^2}{2} \\Delta_\\mathbf{x} p_s(\\mathbf{x})\n$$\n\nwhose stationary distribution satisfies:\n\n$$\np_\\infty(\\mathbf{x}) \\propto \\exp \\left( -\\frac{|R_\\theta(\\mathbf{x})|^2}{\\beta^2} \\right)\n$$\n\nWe will provide the full derivation in the appendix and highlight its connection to gradient flows in the space of probability measures [1,2].\n\n---\n\n## Response to Q2: Consistency of Residual Structure Across Assumption 3.1 and Propositions\n\n**Summary:  \nWe clarify that the apparent inconsistency arises from an implicit linearization in the residual formulation. We will make this assumption explicit and ensure consistency across the assumptions and propositions.**\n\nWe recognize that the formulation of $R_\\theta(\\mathbf{x})$ in Assumption 3.1 and its use in Proposition 3.1, 3.2 appear inconsistent at first glance. We take this opportunity to clarify the derivation and underlying rationale.\n\nWe first clarify the structure of the residual function $R_\\theta(\\mathbf{x}) = a(\\theta)^\\top \\phi(\\mathbf{x}; \\theta)$ introduced in Assumption 3.1, and actively utilize it in our response.\n\nIn general, a residual function can be expressed as:\n\n$$\nR_\\theta(\\mathbf{x}) = (\\partial_t f_\\theta + N_x[f_\\theta])(\\mathbf{x})\n$$\n\nwhere $\\partial_t$ denotes the partial derivative with respect to the PDE’s time domain, and $N_x$ is a nonlinear operator over the spatial domain.\n\nWe define a functional operator $g$ acting on the function space, such that:\n\n$$\nR_\\theta(\\mathbf{x}) = (g(f_\\theta))(\\mathbf{x})\n$$\n\nThen, for a small perturbation $\\Delta\\theta$, we approximate:\n\n$$\ng(f_{\\theta+\\Delta\\theta}) \\approx g(f_\\theta) + D_g(f_\\theta)(f_{\\theta+\\Delta\\theta} - f_\\theta)\n$$\n\nwhere $D_g$ is the Fréchet derivative of the operator $g$. Suppose that $g$ admits the following expansion:\n\n$$\ng(f) = G(f, \\partial_t f, \\partial_x f, \\partial_{tt} f, \\partial_{tx} f, \\dots)\n$$\n\nfor some smooth function $G$. Then the perturbation becomes:\n\n$$\nD_g(f_\\theta)(f_{\\theta+\\Delta\\theta} - f_\\theta) \\approx g(f_{\\theta+\\Delta\\theta}) - g(f_\\theta) \\approx \\sum_k \\frac{\\partial G}{\\partial (\\partial_k f_\\theta)} \\Delta(\\partial_k f_\\theta) \\equiv a(\\theta)^\\top \\phi\n$$\n\nWe assume that $a(\\theta)$ is constant or locally linear, enabling a unified linearized form of the residual. This simplification ensures consistency across Assumption 3.1 and Propositions 3.1–3.2, and underpins our curvature analysis.\n\nTo provide further intuition, we illustrate this within the context of two representative PDEs:\n\n- **Convection equation**:  \n $g(f) = \\partial_t f + \\mu \\partial_x f$ → $a(\\theta) = (0, 1, \\mu, 0, 0, \\dots)$ is constant.\n\n- **Burgers’ equation**:  \n  $g(f) = \\partial_t f + f \\partial_x f - \\frac{0.01}{\\pi} \\partial_{xx} f$ → $a(\\theta) = (\\partial_x f_\\theta, 1, f_\\theta, 0,0,0,-\\frac{0.01}{\\pi},\\dots)$ is $\\theta$-dependent.\n\nUnder the assumption that the neural network $f_\\theta$ is sufficiently expressive and smooth with respect to $\\theta$, we can locally linearize $a(\\theta)$ around a fixed $\\theta$.\n\n\nWhen $a(\\theta)$ is affine, its Jacobian $J_a=\\left(\\frac{\\partial{a_i}}{\\partial\\theta_j}\\right)$ is constant and Hessian $H_\\theta(a) = 0$, where the derivatives are defined elementwise.\n\nWe consider the weighted loss:\n\n$$\n\\mathcal{L}(\\theta) = \\sum_{n=1}^N w_n |a(\\theta)^\\top \\phi_n|^2,\\quad \\phi_n := \\phi(\\mathbf{x}_n)\n$$\n\nThen the Hessian is:\n\n$$\nH_\\theta \\mathcal{L}(\\theta) = 2 \\sum_{n=1}^N w_n \\left[ (J_a^\\top \\phi_n)(J_a^\\top \\phi_n)^\\top + (a^\\top \\phi_n) \\sum_{i=1}^D \\phi_{n,i} H_\\theta(a_i) \\right]\n$$\n\nIf $a(\\theta)$ is linear, the second term vanishes, and we obtain:\n\n$$\nH_\\theta \\mathcal{L}(\\theta) = 2 J_a^\\top \\left( \\sum_{n=1}^N w_n \\phi_n \\phi_n^\\top \\right) J_a\n$$\n\nThis Hessian structure allows us to recover similar proposition results as in the paper, with the inclusion of the Jacobian term $J_a.$ We will revise the manuscript to explicitly include this structural assumption in Assumption 3.1 and ensure that its implications are consistently reflected in all related propositions.\n\n---\n\n## Response to Q3: Consistency of Assumptions 3.2 and 3.3 in the $\\beta \\to 0$ Regime\n\n**Summary:  \nWhile the inequality in Assumption 3.3 may not generally hold, we demonstrate that it becomes statistically valid in the low-temperature, high-dimensional regime. We will clarify this condition in the revised version.**\n\nWe thank the reviewer for pointing out a potential tension between Assumptions 3.2 (heavy-tailed feature norms) and 3.3 (residual maximizer closely follows feature norm maximizer), especially in relation to the inequality $\\max \\||\\phi(\\mathbf{x})\\||^2 \\gg \\lambda_{\\max}(\\Sigma_\\phi)$.\n\nWe agree with the reviewer that this inequality may not hold in general under isotropic feature settings. However, under a highly structured regime, our analysis ensures that both assumptions can be approximately satisfied:\n\n(i) the **low-temperature limit** $\\beta \\to 0$ of the residual-weighted Langevin sampler, and  \n(ii) the **high-dimensional setting** where heavy-tailed distributions emerge naturally from overparameterized neural networks.\n\nIn this regime, Assumption 3.3 becomes statistically valid and meaningful. Our residual function is given by:\n\n$$\nR_\\theta(x) = a(\\theta)^\\top \\phi(\\mathbf{x}) \\quad \\Rightarrow \\quad |R_\\theta(x)|^2 = \\phi(\\mathbf{x})^\\top a(\\theta) a(\\theta)^\\top \\phi(\\mathbf{x}) = \\phi(\\mathbf{x})^\\top Q \\phi(\\mathbf{x}),\n$$\n\nwhere $Q$ is a rank-one positive semi-definite matrix. For $\\phi(\\mathbf{x})$ aligned with $a(\\theta)$, this inner product approximates $\\||\\phi(\\mathbf{x})\\||^2$. Thus, the key question becomes:\n\n$$\n\\text{Under what conditions does } \\phi(\\mathbf{x})^\\top Q \\phi(\\mathbf{x}) \\approx \\||\\phi(\\mathbf{x})\\||^2 \\text{ hold?}\n$$\n\n1. **Directional distortion by residual-weighted sampling**:  \n   In the limit $\\beta \\to 0$, Langevin dynamics increasingly samples from regions where $R_\\theta(\\mathbf{x})$ is large. Since $R_\\theta(\\mathbf{x})$ is linear in $\\phi(\\mathbf{x})$, this implicitly biases the sampling toward feature vectors aligned with $a(\\theta)$ and with large norm.\n\n2. **High-dimensional ($D \\gg 1$) feature decorrelation**:  \n   In high-dimensional settings, $\\phi(\\mathbf{x})$ vectors are nearly orthogonal. Therefore, the only way for $|R_\\theta(\\mathbf{x})|^2$ to be large is for $\\phi(\\mathbf{x})$ to be both high-norm and well-aligned, which makes $\\arg\\max |R_\\theta(\\mathbf{x})|^2 \\approx \\arg\\max \\||\\phi(\\mathbf{x})\\||^2$ statistically plausible.\n\n---\n\nThis assumption is both conceptually sound and crucial for the curvature analysis. Specifically:\n\n- In the **low-temperature regime**, the empirical Hessian becomes:\n\n$$\nH_\\theta \\mathcal{L} \\approx 2 w_{\\max} J_a^\\top \\phi_{\\max} \\phi_{\\max}^\\top J_a, \\quad \\phi_{\\max} = \\max \\||\\phi(\\mathbf{x})\\||\n$$\n\n- In contrast, under **uniform sampling**:\n\n$$\nH_\\theta \\mathcal{L} \\approx 2 J_a^\\top \\Sigma_\\phi J_a, \\quad \\Sigma_\\phi = \\mathbb{E}_\\mathbf{x}[\\phi(\\mathbf{x})\\phi(\\mathbf{x})^\\top]\n$$\n\n- For feature vectors with heavy-tailed norms (Assumption 3.2), one can show:\n\n$$\n\\phi_{\\max}^2 \\sim D^{2/\\alpha}, \\quad \\lambda_{\\max}(\\Sigma_\\phi) \\sim \\frac{D^{2/\\alpha}}{\\log D}\n$$\n\nwhich yields:\n\n$$\n\\|\\| \\phi(\\mathbf{x_{\\max}}) \\|\\|^2  \\approx \\phi_{\\max}^2 \\gg \\lambda_{\\max}(\\Sigma_\\phi), \\quad \\mathbf{x_{\\max}} =\\arg\\max |R_{\\theta}(\\mathbf{x})|^2\n$$\n\n\nas $D\\to\\infty$. This is supported by extreme value theory, where the maximum grows faster than the mean energy [3–5].\n\n*Appendix A.3 (Fig. 7)* empirically shows increasing alignment between $\\arg\\max R_\\theta$ and $\\arg\\max \\|\\|\\phi\\|\\|$ as model size increases, supporting Assumption 3.3.\n\nWe will revise the manuscript to make this connection explicit and emphasize that the inequality is not an inconsistency but rather a natural consequence of our modeling and sampling design.\n\n\n---\n\n[1] Jordan et al., SIAM J. Math. Anal., 1998.  \n[2] Ambrosio et al., Springer, 2008.  \n[3] Vershynin, *High-Dimensional Probability*, Cambridge Univ. Press, 2018.  \n[4] Bobkov & Ledoux, Geom. Funct. Anal., 2000.  \n[5] Dasgupta, FOCS, 1999."}]

Author rebuttal:

Meta-review: The paper presents an adaptive collocation point sampling method for physics-informed neural networks (PINNs). The authors carefully analyze existing high-residual adaptive sampling approaches and point out that they can lead to stability issues. Their proposed Langevin dynamics–based method is shown to be robust across a range of learning rates and model complexities.

All reviewers agree that this is a solid, well-motivated, and clearly written paper. The numerous interactions during the rebuttal provided valuable additional results—both theoretical and numerical—and the reviewers acknowledge that their concerns have been satisfactorily addressed. I therefore recommend acceptance.

I strongly encourage the authors to integrate the rebuttal content, in particular the newly updated proofs and discussions, into the camera-ready version of the paper.
Decision: Accept (spotlight)