Learn Singularly Perturbed Solutions via Homotopy Dynamics
Authors: Chuqi CHEN
Yahong Yang
Yang Xiang
Wenrui Hao
Venue: ICML 2025 poster (2025)
Keywords: scientific machine learning
parametric problems
homotopy dynamics
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: We introduce a homotopy dynamics-based method that enhances neural network training for sharp interface PDEs, achieving faster convergence and higher accuracy.

Abstract: Solving partial differential equations (PDEs) using neural networks has become a central focus in scientific machine learning. Training neural networks for singularly perturbed problems is particularly challenging due to certain parameters in the PDEs that introduce near-singularities in the loss function. In this study, we overcome this challenge by introducing a novel method based on homotopy dynamics to effectively manipulate these parameters.  From a theoretical perspective, we analyze the effects of these parameters on training difficulty in these singularly perturbed problems and establish the convergence of the proposed homotopy dynamics method. Experimentally, we demonstrate that our approach significantly accelerates convergence and improves the accuracy of these singularly perturbed problems. These findings present an efficient optimization strategy leveraging homotopy dynamics, offering a robust framework to extend the applicability of neural networks for solving singularly perturbed differential equations.

Learn Singularly Perturbed Solutions via Homotopy Dynamics
Chuqi Chen' Yahong Yang? 2
Yang Xiang 13 Wenrui Hao
Abstract
Yang & Xiang, 2022: Hao al., 2025). However, despite these advancements, numerically solving such fundamental physical equations remains a challenging task. Existing neural network betwork-based PDE solvers can be broadly divided into two categories. The first is solution approximation, which focuses on directly approximating PDE solutions us- ing methods such as PINNs (Raissi et al., 2019; Kariiadaii et al., 2021: Cuomo et al., 2022), the Deep Ritz Method (E & Yu, 2018), and random feature models (Chen et al., 2022: Dong & Wang, 2023; Sun al., 2024; Chen et al., 2024b), The second category, operator learning, aims to approximate the input-to-solution mapping, with representative methods including DeepONet (Lu et al., 2021a) and FNO (Li et al., 2021), as well as various extensions for broader operator classes (He et al., 2024; Lan et al., 2023; Li et al., 2023; Geng et al., 2024).
Solving partial differential equations (PDEs) us- ing neural networks has become a central focus in scientific machine learning. Training neural networks for singularly perturbed problems is par- ticularly challenging due to certain parameters in the PDEs that introduce near-singularittes in the loss function. In this study, we overcome this chal- lenge by introducing novel method based on ho- motopy dynamics to effectively manipulate these parameters. From theoretical perspective, we analyze the effects of these parameters training difficulty in these singularly perturbed problems and establish the convergence of the proposed ho- motopy dynamies method. Experimentally, we demonstrate that our approach significantly ac- celerates convergence and improves the accuracy of these singularly perturbed problems. These findings present an efficient optimization strategy leveraging homotopy dynamies, offering robust framework to extend the applicability of neural networks for solving singularly perturbed differ- ential equations.
However, the optimization challenges in solving PDEs sig- nificantly limit the applicability and development of neu- ral network- based methods. Studies have shown that the loss functions for solving PDEs are often difficult to mini- mize, even in simple scenarios (Krishnapriyan et al., 2021; Rathore et al., 2024; Xu et al., 2024; Chen et al., 2024b;a). For example, small diffusion coefficients in the Allen-Call Cahn equation (Allen & Cahn, 1975), small viscosity terms in the Burgers equation (Burgers, 1948), and large wave numbers in the Helmholtz equation (Hilbert, 1985). In these types of equations, the parameters sigmificantly influence the solution behavior. In the parameter sharpens the solution, often resulting in near- singular structures. In the Helmholtz equation, increasing oscillations. These
1. Introduction
The study of Partial Differential Equations (PDEs) serves as a cornerstone for numerous scientific and engineering disciplines. In recent years, leveraging neural network ar- chitectures to solve PDEs has gained significant attention, particularly in handling complex domains and incorporating empirical data. Theoretically, neural networks have the po- tential to overcome the curse of dimensionality when solving PDEs (Han al., 2018; Siegel & Xu, 2020; Lu et al., 2021b;
and Burgers equations, decreasing
the parameter induces effects complicate the loss landscape, making optimization challenging and often causing slow convergence, inaccurate solutions, or even divergence.
The root of this challenge lies in the highly complex energy landscape of the loss function near singularities. which sig- nificantly exacerbates optimization difficulties (Karniadakis et al., 2021; Xu et al., 2024). To address these challenges, two main strategies have been proposed. The first strat- egy is resampling, which involves introducing additional collocation points in regions with low regularity to better capture the solution's complexity (Wight & Zhao, 2020; Gao et al., 2024: Zhang et al.. 2025). However, resampling- based methods typically require a large number of sample
'Dearrtment of Mathematies, The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong SAR, China DDepartment of Mathematics. The Pennsylvania State University, PA. USA Algorithms of Machine Learning and Autonomous Driv- ing Research Lab, HKUST Shenzhen-Hong Kong Collaborative Innovation Research Institute, Futian, Shenzhen, China Correspondence to: Yahong Yang <yxy5498@psul edu.edu>.
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Homotopy Dynamics
Loss Landscape for problem: (u( f(u)
Algorithm
Phase
Direct training for large &o
min 11 (1o(((e)) (Mea)
f(u) l1
= 0.1
Homotopy Evolution
Phase
== 0.03
Homotopy dynamics path tracking
With Hard train
H(ug, E) = L.(u,) - f(u))
E 0.01
de -H
Stuck in local minima
u(E = 0) =
Wrong solution
Ground truth solution
Figure 1. Framework of homotopy dynamies for solving singularly perturhed problems.
points, leading to substantial memory consumption, and the sampling process becomes increasingly complicated in high- dimensional settings. The second strategy is the design of multiscale neural network architectures (Wang, 2020; Liu, 2020; Liu et al., 2024b; Hao et al., 2024; Wang et al., 2021; Huang et al., 2025). These approaches generally require certain priori knowledge of the solution properties, impose specific constraints on the network design, and are highly sensitive the selection of hyperparameters.
theoretical analysis. In contrast, our work is the first theo- retically demonstrate that in singularly perturbed problems, smaller values lead to greater training difficulty (Theorem 1). While both approaches share the curriculum-style moti- vation, our method differs in design and rigor: we construct a continuous homotopy path in parameter space with con- vergence guarantees and introduce a principled strategy for choosing the homotopy step size Ae (Theorem 2), which is absent in (Krishnapriyan al., 2021).
In this paper, introduce novel approach based on ho- motopy dynamics to gradually reshape the complex energy landscape with respect to a specific coefficient. Rather than directly computing solutions near singularities, we leverage homotopy dynamics to trace a solution path that approxi- mates them more effectively. More specifically, we inves- tigate the training challenges introduced by a parameter in the PDE residual term within the loss functions. As decreases, the problem becomes more significantly difficult to solve. To understand this effect, we provide a theoretical analysis of how & influences the convergence the training process. To address this issue, we propose a novel method called Homotopy Dynamics. The key idea is to first train the neural network PDEs with large e, where the problem is easier to learn and training is more efficient. Then, grad- ually and adaptively adjust the neural network according to the evolution of the homotopy dynamics, guiding toward its target value (as illustrated in Figure 1). Although the homotopy approach has been used to train neural networks (Chen & Hao, 2019; Yang al., 2025), this work the first apply homotopy dynamics sharp interface problems in PDEs through the parameter
Contributions. Our key contributions are summarized as follows:
We propose Homotopy Dynamics, novel method for solving singularly perturbed PDFs with neural networks, achieving improved training performance (Section 3). We provide theoretical analysis of how the PDE param- eter e affects training difficulty, and establish the conver- gence of our method (Section 4).
We validate the method on diverse problems, including the Allen-Cahn equation, high-dimensional Helmholtz equation, and operator learning for Burgers' equation (Section 5).
2. Problem Setup
We begin by introducing the singularly perturbed problems studied in this work, followed by the neural network-based solution approach and the training challenges that motivate our method.
2.1. Singularly perturbed Problems
The form of the singularly perturbed problem defined as follows:
A related idea appears in (Krishnapriyan et al., 2021) as Cur- riculum PINN Regularization, where PDE parameters are observed to influence PINN performance, though without
{ -=== f()), in Ω (1)
Bu - g(x), on 82,
2

---

Homotopy Dynamics
where Lg is differential operator defining the PDE with certain parameters, B is operator associated with the boundary and/or initial conditions, and 52 S R In the con- sidered PDEs, the parameter e governs the complexity of the solution, with smaller values generally leading to more challenging behaviors. For example, the Allen-Cahn equation (4), e represents the interfacial width parameter, where smaller results in sharper transition layers. In the Burgers equation (28), corresponds to the viscosity efficient, with small values leading to steep gradients or shock-like structures. In the Helmholtz equation (27), is inversely related to the wave number, and decreasing yields freuueccy oscillations.
training error
100
0.03
0.1
157
30"
13"
Itcrations
Figure 2. Training curves for different values of in solving the Allen- Cahn steady-state equation. decreases, the training emror increases, indicating that the training process becomes pro- gressively more difficult.
In all cases, as becomes small, the solution exhibits increased complexity-whetheee through sharp interfaces, steep gradients, or freeerryyyency structures poocturre-ppoingg sig- nificant challenges for neural network-based solvers. More details will be provided in the following section.
2.3. Challenges in Training Neural Networks
2.2. Neural Networks for Solving PDEs
In this paper, we consider the following singularly perturbed elliptic problem:
In this section, we focus solution approximation rather than operator learning for simplicity, specifically using neural network to approximate the PDE solution. In Section 5, we will demonstrate that our Homotopy Dynamics can also generalize to the operator learning case. The PDE prob- lem typically reformulated the following non-linear least-squares problem, aiming determine the parameters 0 of the neural network u(x; 0) (commonly multi-layer perceptron, MLP):
AeAu(z) f(u), x 2 u(x) = g()),
(3)
x 0n,
where is problem parameter that influences both the structure of the solution and the difficulty of training neural network solvers. As a representative example, we focus on the steady-state Allen Mahn equation in one spatial dimen- sion:
Thes
{ eul(x) = =n3 u, E [0, u (0)=-1, u(1) 11
min I(0)) BER*
(Cuu(ar); ( [u((a) 0)))
(4)
2nres
Tas
where the parameter e controls the width of the internal interface. As e decreases, the interface becomes increasingly sharp, resulting in solution with higher gradients near the transition region. The analytic steady-state solution of this problem is given by
Bu(zz); 0) g(xx))) (2)
+ X 21The
Liee
Here Liss is the PDE residual loss, Loce is the boundary loss and À is a constant used to balance these two terms. The sets {xr) Fie represent represent the interior sample points, and { represent the boundary sample points. We also introduce the lg relative error (L2RE) to evaluate the discrepancy between the neural network solution and the ground truth, defined as
0.5
u(x) tanh
/
where the interface centered at x = 0.5. As shown in Figure 3, the solution becomes sharper as e becomes smaller.
To show the challenges in the optimization problem defined in (2), we present the training curves for varying values of in Figure 2. As decreases, training errors increase. This is due to the significantly increased training difficulty and slower convergence for smaller s. In the subsequent sections, we analyze the underlying reasons for this phenomenon and introduce homotopy hyynaiic-basee approach to address the challenge.
I|we
u ||2
L2RE
112
where ug is the neural network solution and u is the ground truth.
3

---

Homotopy Dynamics
3. Homotopy Dynamies
with the initial condition 0(eg) Oo and t stands for Moore-Penrose inverse (Ben-Israel & Greville, 2006). Thus, to solve the singularly perturbed problem (3) where is small, we can first solve (3) with large using the loss function (2). Then, by following the homotopy dynamics path tracking (8), we can progressively obtain the solution for smaller values of e, ultimately solving the singularly perturbed problem.
To address training difficulties in neural networks for sin- gularly perturbed problems, we introduce a novel approach termed homotopy dynamics.
3.1. Homotopy Path Tracking
First, we introduce the homotopy function below:
In particular, path tracking in homotopy dynamics reduces to solving initial value problems numerically, with the start points serving as the initial conditions. For different neural network architectures, we propose two strategies, which are summarized in Algorithm 1.
Eu, = L.u fu)) = O,
(5)
where E is the parameter in the PDEs. Specifically, this formulation represents the PDE problem - f()). In this context, e is treated as path-tracking parameter. At e = Eo, we assume that the solutions to H Eo) = 0 are either known or can be easily approximated by neural networks. These solutions are referred to as the starting points, At E = et, the original system we aim to solve is recovered, which referred to as the target system. Therefore, solving the target system involves tracking the solutions of Hu,, E) 0 from = Eo, where the solutions are known, to = ==,, where the solutions are sought.
One is to solve the initial value problem by using the forward Euler method, as follows:
((sx) (ex-1) - DEu Vou(Ex-1)) He, (9)
where Ack = Ek - Ek-1- This approach is effective for small neural networks, as the pseudo-inverse is easy to com- pute.
The other approach is introduce the Homotopy Loss in the optimization, formulated as:
The process of path tracking between en and at is governed by solving the Davidenko differential equation:
L (0(E)) =L + ALre + aLH,, (10) ((E))
min
dH e (H(u(e); e) du(e) aH(u(e),e) de
0,
Ou dc
Oe
where Lu is defined in Eq. (12), and LH, is the loss function from Homotopy Dynamics, which is
(6) with the initial condition u(so) ug- Thus, path track- ing reduces to numerically solving an initial value problem, with the starting points acting as the initial conditions. Addi- tionally, the boundary condition in (3) should be taken into account when solving the initial value problem numerically.
LH. 27i
Thes Σ
Au
Aek + He((u(())))))
i=1
This approach is suitable for large neural networks, as it does not require the computation of the pseudo-inverse, and Au = U0(e2) 11((ex-))*
3.2. Incorporating Homotopy Dynamics into Neural Network Training
Algorithm 1 Homotopy Dynamics Path Tracking
To enhance the neural network training process, we incor- porate homotopy dynamics by gradually transitioning the network from an easier problem (with a larger én) to the original target problem (with e*). This approach helps mit- igate the challenges associated with training networks for problems involving small values of e, where solutions be- come increasingly sharp or oscillation and harder to com- pute. Specifically, we denote the neural network solution for (3) as u(x; (E)). The homotopy path tracking for training neural networks can then be refined as:
input tolerance list of parameters eo, , Phase I: Directly train NN for large su while (((()))) > T do min L(O(Eu))
end while
Phase II: Homotopy dynamics path tracking
for k- 1,...,n do
AEc = Ek Ek-1
Strategy Numerical solution via Forward Euler:
(EE)) = 0(s--1) A Vou(Ex-I)'H,,
Strategy 2. OOptimization using homotopy loss:
while LHom [Hoe))(E > 7 do min [mlm(((Es)) end while end for output 10((s))
(E( H1 -0
Vou-
de
where H, = aH du H.
ah de and Vou represents the Jaco-
bian with respect to the neural network parameters 0. Thus we can derive the homotopy dynamics system as:
dd(E) = -HH. He, EE [Eo,E"], de
If LH and AL are omitted in strategy 2, then strategy can be viewed as an alternative approach to solving the linear
(8)

---

Homotopy Dynamics
4. Theoretical analysis
Table/ 1. Training loss and L2 error (L2RE) for classical training vs. homotopy dynamics under different E. Homotopy dynamics
In this section, we provide theoretical support for homotopy dynamics. In the first part, we demonstrate that for certain PDEs with small parameters, direct training using PINN methods is highly challenging. This analysis is based on the neural tangent kernel (NTK) framework (Allen-Zhu et al., 2019). In the second part, we show that homotopy dynamics will converge to the solution with small parameter e, pro- vided that the dynamic step size is sufficiently small and the initial solution has been well learned by the neural network.
achieves consistently lower loss and error.
g 0.1
0.03
=
0.01
Loss
LIRE
Loss
L2RE
Loss
LIRE
Classical
5.000-6
171o-2
7.76e-4
1.11
7.21
8.17e-1
Homotopy
5,00366
177622
7,45-88
9.83-3
4,63e-8
8,08e-3
Inibill 0.1
<- 0.03
End 0.01
Ground truth
4.1. Challenges in Training Neural Network with Small Certain Parameters
022
-335
-329
Let us consider training neural networks without homotopy dynamics. The corresponding loss function can be expressed
-375
-10
(o)
(T)
(12)
Figure 3. Evolution of the Homotopy dynamics for steady state ID Allen-Cahn equation. The 1.2RE for e 0.01 is 8.08e
2n
where {< represents the training data used to optimize the neural network. Here, we assume that the parameter e in the PDE appears only in the interior terms and not the boundary conditions. Therefore, in this section, we omit the effect of boundary conditions, as the behavior at the boundary remains unchanged for any given E.
system given in Eq. (7). However, for certain PDEs or larger neural networks, directly solving Eq. (7) is unstable because the term H Vou contains many small singular values, caus- ing conventional methods (e.g., using SVD) to incur large errors. Therefore, we opt to solve an optimization problem in the traditional manner. Within this framework, strategy 2 follows the same dynamic process as strategy 1, yielding more stable training. Including Lu and ALa. ensures that the obtained solution satisfies the target PDE. Conversely, even if H( Huu, E) Const F O, the solutions still follow the same homotopy dynamics since they share the same LHt.-
Furthermore, to simplify the notation, we use n instead of Thes and denote æ². simply as æ; comparing with Eq. (2).
In the classical approach, such a loss function is optimized using gradient descent, stochastic gradient descent, or Adam. Considering the training process of gradient descent in its continuous form, it can be expressed as:
de
= V (0) Example: 1D Allen-Cahn steady-state equation. We demonstrate our proposed method on the one-dimensional Allen-Cahn steady-state equation by defining the following homotopy function:
dt
"
Σ H Huuo( E00 (): e))
n
H (uo(x), S
(13)
H (uu, s) = u((((()))) + nn ug = 0.
(II)
where in this section is the time of the gradient descent flow, oug is the functional variational corresponding to ug. and
Following the homotopy dynamics in Eq. (8), we set the initial value at e = 0.1 and gradually decrease it to the final value e,, = 0.01. The initial solution, O(Eo), is obtained using the standard training process by directly minimizing (2). The results and the evolution process are presented in Table 1 and Figure 3. These results show that when is large, the original training method achieves relatively small error, leading to an accurate solution. However, as s decreases, the error increases, which reduces the accuracy of the solution. In contrast, the homotopy dynamics-bris based approach maintains accuracy effectively as e decreases.
H ())( =
((o(()))))
= =1g- -D
(14)
and
4 = H [uo (x, 0).=)] E E
RXIXR De Rax
(15)
where D. represents the discrete form of the variation of PDEs in different scenarios. Furthermore,
s [Vouo(x1).... Vouo(
(16)
5

---

Homotopy Dynamics
1e10
Largest Eigenvalue of Dr
et al., 2023) discusses the positivity of the gradient ker- nel in PINNs for solving heat equations. Therefore, can reasonably assume that SS' is strictly positive ma- trix. In Appendix A.1, we present specific scenario where > O holds with high probability.
25-
C- 0.01; Aras 322.04
0.03: Apar 21267.00 03
20
0:1: Arwy 3752588.92
Ami ( SSS
15
This theorem demonstrates that the smallest eigenvalue of the kernel directly affects the training speed. Equation (20) shows that the upper bound of Amio (K&) can be influenced by Amax (D.D. In many PDE settings, the maximum eigenvalue Amax D D.D! tends to be small when is small. For example, in this paper, we consider the Allen-Cahn equation, given by
00
0.6
08
100
2 Au + f(u) = =0,
Figure Largest cigenvalue of Dz (21) for different e. A smaller results in smaller largest eigenvalue of (21), leading to slower convergence rate and increased difficulty training.
where f(u) = u² u. In this case, De corresponds to the discrete form of the operator - + f'(u), which can be written as
(21) According to (Morton & Mayers, 2005), the discrete Lapla- cian - Adis is strictly positive. Specifically, in the one- dimensional case, its largest eigenvalue is given by
D: - diag (l,
Therefore, we obtain
d0 (L) VoLn (0) dt
dt
H (() H evo()),)
2
2=2n² n cos
2 2n + 1
- DDSS D.
(17)
n
which is close 4=2n² as n is large enough.
Hence, the kernel of the gradient descent update is given by
Moreover, since f'(u(x;)) lies within the interval [-1, 2]. when e is large (i.e., close to 1), the largest eigenvalue of De becomes very large regardless of the sampling loca- tions {=.]"_, (see Figure for the case n 200. Other equations exhibit similar behavior. Please refer to the Ap- pendix B for detail.) Consequently, by Theorem 4.1 the upper bound on the smallest eigenvalue of K2 also be- comes spgse-specfiially.. it is of order n in this case due to Weyl's inequalities. As a result, the training speed can achieve a rate of exp(-Cn't) (see Eq. (19)), which is rapid and indicates that training is relatively easy. In contrast, the lower bound for the training speed is given by exp(-Ct//
K, : D.SS D
(18)
The following theorem provides an upper bound for the smallest eigenvalue of the kernel and its role in the gradient descent dynamics:
Theorem 4.1 (Effectivenesss of Training via the Eigenvalue 0 and D3 is non- singular, and let = 2 0 be constant. Then, we have Amin K2) > 0, and there exists T > 0 such that
of the Kernel), Suppose Amin (SS
L (f())) < LH((00)) exp (- Amina (k)nn
(19)
The fastest rate is attained in the special situation where there exists a nonzero vector that is eigenvector corre- sponding to the largest eigenvalue of D D and, simultane- ously, D Dax: is an eigenvector corresponding to the smallest eigenvalue S S. This scenario may occur under particu- lar configurations of S and D., which in turn depend on the underlying PDE and the distribution of the sampling points.
for all t [), Furthermore,
Xmin (SS )amin Ammin (D.D.)
I Amin (K2) S Amin (SST) (D.D..) (20)
Remark 4.2. For TS previous works such (Luo & Yang, 2020; Allen-Zhu et al., 2019; Arora et al., 2019; Cao & Gu, 2020; Yang et al., 2025; Du et al., 2019; Li et al., 2020) demonstrate that it becomes positive when the width of the neural network sufficiently large with ReLU acti- vation functions or smooth functions. Additionally, (Gao
On the other hand, when is small (close to 0), the largest eigenvalue of D is only of order I with respect to n. Conse- quently, the upper bound on the smallest eigenvalue of K; no longer scales as a constant with respect to n. In this case, the training speed is reduced to exp(-Ct/n (as indicated in Eq. (19)), which is significantly slower and suggests that
6

---

Homotopy Dynamics
E(S) = 1.0
E(S) = s 0.1
{()) = s 0.08
E(s) = 0.05
8(S) 0.05; 0.001
Figure 5. 2D Allen Cahn Equaiton. (Top) Evolution of the Homotopy Dynamcis. (Bottom) Plot for Cross-section of u(x, y) at 0.5 ie.. u(2, y 0.5). The reference solution Uhu (2)) represents the ground truth steady state solution. The L2RE is 8.78e 3. Number of residual points is Thes = 50 x 50.
training becomes difficult. Therefore, while larger E may yield relatively easy training in some instances, a smaller s will invariably lead to challenging training conditions.
Assume there exists constant P such that (Ex Ek+1) S P. camea == [|u(Eo) (5o))l <
and
n
sup (Ex (=====)))uuu(elllll < 0SkSn
4.2. Convergence of Homotopy Dynamics
In this section, we aim to demonstrate that homotopy dynam- ies is a reasonable approach for obtaining the solution when e is small. Recall that Strategy merely an alternative ap- proach for solving the linear system in our framework, with the underlying principles remaining the same. Therefore, our theoretical analysis is primarily based on Strategy 1. For denote u(s) the exact solution of H (u, E) 0 and (s) as its numerical approximation in the simulation. Suppose H(((e), e) = O, and assume that (u(e).)) is invertible. Then, the dynamical system (6) can du be rewritten as
then we have
[|u(En) - -(()))lr(())))
Plo--e) + T(e - 1. (24)
2P
simplicity of notation,
The proof of Theorem 4.3 inspired by (Atkinson al.. 2009).
Theorem 4.3 shows that if is small and the step size (Ek+1 - Ek) is sufficiently small at each step and satisfies (# < P. cumén i.e., the training step size should depend on the Lipschitz constant h(u, s), ensuring stable training, then u(Ex) U(E)) remains small. In other words, when
du de -
(HH(u()),e) (H(E))) h(u(E), s).
d=
(22)
Pg may be bounded. Nonetheless, do not require small constant. For small Ps might be large but remains finite. In this case, one must choose sufficiently small steps 6k+1 E to ensure that the training error stays controlled. The initial error ég can be very small since we use a neural network to approximate the solution PDEs for large s, where learning is effective.
Applying Euler's method to this dynamical system, we ob- tain
be exactly zero; it only needs be
U(ER+1) U(Ex) + (#x41 - Ex)h(U(Ee)
(23)
The following theorem shows that if u(Eo) - (Eo) is small and the step size (#x+1 5x) is sufficiently small at each step, then u(Ex) - (EX) remains small.
The total error eo consists of approximation, generalization, and training emrors. The training error can be effectively controlled when s is large (Theorem 4.1), while the approxi- mation and generalization errors remain small if the sample size sufficiently large and the neural network expressive enough. Theoretical justifications are provided in (Yang
Theorem 4.3 (Convergence of Homotopy Dynamics). Sup- pose h(e u) is a continuous operator for 0 < en S Eo and E (2),), and
[M[ e) h(u2, ( (2 & Pellur - 42ll

---

Homotopy Dynamics
Table 2. Amin ( for different initialization for E = 0.01 in Equation (4).
Initializaiinn
Kavier
Home
0.1
Hom e 0.05
Hom e = 0.03
Hom 0.02
(kn
7
2.
3.77e-5
137e-
1.3se2
Table 3. 2D Allen-Cahn Equation. Relative L error comparison across various training strategies. Time Seq, (Wight Zhao, 2020; Matley Ghosh, 2022)
Method
Original PINN Curriculum (Krishnapriyan et al., 2021)
Resampling
Homotopy
L2RE
9.56c-1
8.890-1
8.95p-2
8.25c-1
8.78c-3
et al., 2023; Yang & He, 2024) and further discussed in Appendix A.4.
ues to evolve toward 0. As shown in Figure 5 and Table 3, our method achieves superior solution accuracy compared to the other approaches.
For Strategy 2, we enforee H; 0 by retraining the net- work from scratch at each homotopy step rather than by integrating the path with Euler's method. Each iteration uses the previous solution as the initial guess for the next, producing progressively better starting points. To illustrate, we consider the ID Allen-Cahn equation (Eq. (4)) with 8 = 0.01. As shown in Table 2, smaller increments As increase Amin (K&), indicating improved conditioning and 0.02 yielded
5.2. High--Dimensiona Helmholtz Equation
One of the advantages solving differential equations using neural networks their potential to overcome the curse of dimensionality and tackle high-dimensional problems. In this example, demonstrate this capability by comparing the performance of the standard PINN approach with our proposed homotopy- based training method on the following high-dimensional Helmholtz equation:
faster convergence. In our experiments, the best initialization for 0.01, highlighting that select- ing a homotopy parameter close to the final target provides the most effective initial guess-this encapsulates the key idea behind Strategy 2.
==AAu u a u O in S
== g on a., (27)
where Ω (-1, 1jo. This problem admits the exact solution
5. Experiments
1 d -
2
u(x) = sin
We conduct several experiments across different problem settings to assess the efficiency of our proposed method. In the following experiments, we adopt Strategy 2 for all train- ing procedures. Detailed descriptions of the experimental settings are provided in Appendix B.
The corresponding homotopy is defined by
1 H(u,)) = = + d
5.1. 2D Allen Cahn Equation
Here, we start the homotopy training by so 1. The numerical results for dimension d = 20 are reported in Ta- ble 4, where we compare the relative 1,2 enrors obtained by classical PINN training and the proposed homotopy dy- namics for different values of e. As shown in the table, the homotopy- based approach achieves consistently lower er- rors, especially for small E, where classical training suffers from significant accuracy degradation.
First, we consider the following time-dependent problem:
Wi = e² Au - u(u² 1), (x,0) E [-1, 1]
u(x y, 0) = sin (Tx) sin(Try)
(25)
u(-1,y, y, = u11 y. u(x - 1,t) u(, 1, O.
We aim to find the steady-state solution for this equation with 0.05 and define the homotopy as:
H Hu, (u, 8)
- 8) (s(6) Au uuu² - 1)) + s(u -u),
Table 4. Comparison of relative L errors achieved by classical training and homotopy dynamics for different values in high- dimensional Helmholtz equation.
where: E [0, 1]. Specifically, when s the initial condition Ug is automatically satisfied, and when s 0,
it recovers the steady-state problem. The function e(s) is given by
Dimension d = 20
= 1/2
== 1/20
1750
Classical Training
1.23e
7.21e
9.98e
Homotopy Dynamics
5.86c-1
5.00e-1
5.8%
8, E [0.05, 1]. E(5) {
{
(26)
0.05, s s [0, 0.05].
5.3. Burgers Equation
Here, E(S) varies with s during the first half of the evolution. Once Es)) reaches 0.05, it remains fixed, and only 8 contin-
In this example, we adopt the operator learning framework to solve for the steady-state solution of the Burgers equation,
8

---

Homotopy Dynamics
Towards steady state solution
8(5) 0.5
0.1
((S) = 0.05
Homotopy Evolution
(c)
(d)
0.001
0.01
8(S)
Figure Burgers' Equation (Operator Learning): Steady-state solutions different initializations Mg under varying viscosity s: (a)
0.5, (b) 0.1, (c) 0.05. The results demonstrate that all final test solutions converge the correct steady-state solution. (d) Illustration of the evolution of test initialization up following homotopy dynamics. The number of residual points is fhyes 128.
given by:
Table 5. Homotopy Dynamics Results on operator lcarning for Burgers Equation: Homotopy Loss, Relative L Error, Shock Localization Accuracy
ut + 122
EUge = π (0n(T))cc((TT))) æ E [0, 1]
2
= 0.5
s = 0.1
= 0.03
Homotopy Loss
7,5507
3.40e-7
7.77e-7
u(x, 0) = (I(),
(28)
L2KE
300-3
7.00c44
2.520-2
MSE Distance
1.75--8
9.14c-8
1,20.1
u(0,t) = u 11 = 0,
with Dirichlet boundary conditions, where to E (0;(0,1); is the initial condition and R the viscosity coefficient. We aim learn the operator map- ping the initial condition to the steady-state solution, Gf (0, R); + ((0, 1): R). defined by + Uoo for any r > 0. As shown in Theorem 2.2 of (Kreiss & Kreiss, 1986) and Theorems 2.5 and 2.7 of (Hao Yang, 2019), for any & > O, the steady-state solution independent of the initial condition, with single shock occurring at x, 0.5. Here,
we analyze the training difficulties caused by certain pa- rameters in the PDEs. To overcome these challenges, propose training method based on homotopy dynamics to avoid training original and singularly perturbed problems directly and to improve the training performance of neural networks to solve such problems. Our theoretical analysis supports the convergence of the proposed homotopy dy- namies. Experimental results demonstrate that our method performs well across range of singularly perturbed prob- lems. In solution approximation tasks, it accurately captures the steady-state solutions of the Allen-Cahn equation and effectively handles high- nigh-dimensional Helmholtz equations with large wave numbers. Moreover, in the context of opera- tor learning, it achieves strong performance on the Burgers' equation. Both the theoretical analysis and experimental re- sults consistently validate the effectiveness of our proposed method.
DeepONet (Lu et al., 2021a) the network architecture. The homotopy definition, similar to Equa- tion (25), can be found in Appendix B.6. The results can be found in Figure and Table 5. Experimental results show that the homotopy dynamics strategy performs well in the operator learning setting as well. As shown Appendix 8, DeepONet trained via homotopy dynamics achieves com- parable accuracy but significantly faster inference than the finite difference method.
Looking ahead, it will be valuable to explore whether homo- topy dynamics can be applied to wider range of practical problems. We believe that homotopy dynamics offers a natural entry point for comparing neural-network methods with traditional techniques.
6. Conclusion
In this work, we explore the challenges of using neural net- works solve singularly perturbed problems. Specifically,

---

Homotopy Dynamics
Acknowledgements
Chen, C., Zhou, Q Yang, Y., Xiang, Y., and Luo, T. Quanti- fying training difficulty and accelerating convergence in neural network-based pde solvers. arXiv preprint arXiv:2410.06300, 2024b.
Y.Y. and W.H. was supported by National Institute of General Medical Sciences through grant BBSGM1166894 The work of Y.X. was supported by the Project of Hetao Shenzhen-HKUST Innovation Cooperation Zone HZQB- KCZYB-2020083.
Chen, J., Chi, X., Yang, Z.,, al. Bridging traditional and machine learning- based al algorithms for solving PDEs: the random feature method. J Mach Learn, 1:268-98, 2022.
We would like to acknowledge helpful comments from the anonymous reviewers and area chairs, which have improved this submission.
Chen, Q. and Hao, W. A homotopy training algorithm for fully connected neural networks. Proceedings of the Royal Society A, 475(2231):2009066)) 2019.
Impact Statement
Cuomo, S., Di Cola, S.S Giampaolo, F., Rozza, G., Raissi, M., and Piccialli, F. Scientific Machine Learning Through Physics-Informed Neural Networks: Where We Are and What's Next. J. Sci. Comput., 92(3), 2022.
This paper presents work whose goal is to advance the field of scientific machine learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
2022.
Dong, S. and Wang, Y. A method for computing inverse parametric PDE problems with random-weight neural net- works. Journal of Computational Physics, 489; 12263, 2023.
References
Allen, S. M. and Cahn, W. Coherent and incoherent equilibria in iron-rich iron-aluminum alloys. Acta Metal- lurgica, 2017 1975.
Du, S., Lee, J., Li, H., Wang, L.,, and Zhai, X. Gradient descent finds global minima of deep neural networks. In International conference on machine learning, pp. 1675- 1685. PMLR, 2019.
Allen-Zhu, Z., Li, Y., and Song, 7. A convergence theory for deep learning via vreepparameteriaation. In International conference on machine learning, pp- 242-252. PMLR, 2019.
E, W. and Yu, The Deep Ritz Method: A Deep Learning- Based Numerical Algorithm for Solving Variational Prob- lems. Communications in Mathematics and Statistics, 6 (1):1-12, 2018.
Arora, S., Du, S. S., Hu, W., Li, Z., Salakhutdinov, R. R., and Wang, R. On exact computation with an infinitely wide neural net. Advances in neural information process- ing systems, 32, 2019.
Evans, L. C. Partial differential equations, volume 19. American Mathematical Society, 2022.
Gao, Y. Gu, Y., and Ng. M. Gradient descent finds the global optima of two-layer physics-informed neural net- works. In International Conference Machine Learning, pp. 10676-10707 PMLR, 2023.
Atkinson, K., Han, W., and Stewart, D. E. Numerical solu- tion of ordinary differential equations, volume 81. John Wiley & Sons, 2009.
Ben-Israel, A. and Greville, T. N. Generalized inverses: theory and applications. Springer Science & Business Media, 2006.
Gao, Z., Tang, T., Yan, L., and Zhou, T. Failure-informed adaptive sampling for pinns, part ii: combining with re- sampling and subset simulation. Communications on Applied Mathematics and Computation, 1700117201 2024.
Burgers, J. M. A mathematical model illustrating the theory of turbulence. Advances in applied mechanics, 1:171- 199, 1948.
Geng, Y., Teng, Y., Wang. Z., and Ju, L. A deep learn- ing method for the dynamics of classic and conservative Allen-Cahn equations based on fully-discrete operators. Journal of Computational Physics, 496: 1 12589, 2024.
Cao, Y. and Gu, Q. Generalization error bounds of gra- dient descent for learning over-parameterized deep relu networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp. 3349-3356, 2020.
Glorot, X. and Bengio, Y. Understanding the difficulty training deep feedforward neural networks. In Pro- ceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010.
Chen, C., Yang, Y., Xiang, Y., and Hao, W. Auto- matic differentiation is essential in training neural net- works for solving differential equations. arXiv preprint arXiv:2405. /4099, 2024a.
Grisvard, P. Elliptic problems in nonsmooth domains. SIAM, 2011.
10

---

Homotopy Dynamics
Han, J., Jentzen, A., and E, W. Solving high-dimensional partial differential equations using deep learning. Pro- ceedings of the National Academy of Sciences, 1 115(34); 8510, 2018.
Li, C.-K. and Mathias, R. The Haiskii-Miiskk-Wieladdd theorem-additiv and multiplicative versions. Nu- merische Mathematik, 81:377-413, 1999.
Li, W., Bazant, M. Z., and Zhu, Phase-field deeponet: thysics-informer deep operator neural network for fast simulations of pattern formation governed by gradient flows of free-energy functionals. Computer Methods in Applied Mechanics and Engineering, 416: 16299, 2023.
Hao, W. and Yang, Y. Convergence of a homotopy finite clement method for computing steady states of burgers' equation. ESAIM: Mathematical Modelling and Numeri- cal Analysis, 53(5): 2019.
Hao, W., Li, R. P., Xi, Y., Xu, T., and Yang, Y. Multi- scale neural networks for approximating green's 's func- tions. arXiv preprint (tt::2410.184333 2024.
Li, Y., Luo, T., and Yip, N. Towards an understanding of residual networks using neural tangent hierarchy (nth). arXiv preprint
2020.
Li, Z.,, Kovachki, N. B., Azizzadenesheli, K., liu, B.,, Bhat- tacharya, K., Stuart, A., and Anandkumar, Fourier Neural Operator for Parametric Partial Differential Equa- tions. In International Conference Learning Represen- lations, 2021.
Hao, W., Liu, X., and Yang. Y. Newton informed neural operator for solving nonlinear partial differential equa- tions. In The Shirty-eighth Annual Conference on Neural Information Processing Systems, 2025.
He, J., Liu, X., and Xu, J. Mgno:
Efficient parameterization of linear operators via multigrid. In 12th International Conference on Learning Representations, ICIR 2024, 2024.
Liu, H., Yang, H., Chen, M., Zhao, T., and Liao, W. Deep nonparametric estimation of operators between infinite dimensional spaces. Journal of Machine Learning Re- search, 25(24):1 202,2
Hilbert, D. Methods of mathematical physics. CUP Archive, 1985.
Liu, X., Xu, B., Cao, S., and Zhang. L. Mitigating spectral bias for the multiscale operator learning. Journal of Computational Physics, 50ysics,506:11130 12944, 2024b. ISSN 0021- 9991. doi: hutp:://o...g////00666/....22222 URL.
Huang, J., You, R., and Zhou, T. Frequency-adapiive multi- scale deep neural networks. Computer Methods in Ap- plied Mechanics and Engineering, 437:117751, 2025.
Karniadakis, G. E., Kevrekidis, G., Lu, L., Perdikaris, P., Wang, S., and Yang, L. Physics-informed machine learning. Nature Reviews Physics, 366):422-440, 2021.
Liu, Z. Multi-scale deep neural network (mscalednn) for solving miisson-bottzaam equation in complex domains. Communications in Computational Physics, 28(5), 2020.
Kreiss, G. and Kreiss, H.-O. Convergence steady state of solutions of burgers' equation. Applied Numerical Mathematics, 2(3):161-179, 1986. ISSN 0168-9274.
Lu, L., Jin, P., Pang, G., Zhang. Z., and Karniadakis, G. Learning nonlinear operators via DeepONet based the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218 2229 2021a.
doi: URL sttps:://www.scceeeedircre com/ wweenc/articee/pi////ee. Special Issue in Honor of Milt Rose's Sixtieth Birthday.
Lu, Y., Lu, J., and Wang, M. A priori generalization analysis of the deep Ritz method for solving high dimensional elliptic partial differential equations. In Conference on learning theory, pp. 3 3196-3241. PMLR, 2021b.
Krishnapriyan, A., Gholami, A., Zhe, S., Kirby, R., and Mahoney, M. W. Characterizing possible failure modes in physics-informet neural networks. In Advances in Neural Information Processing Systems, 2021.
Luo, T. and Yang, H. Two-layer neural networks for partial differential equations: Optimization and generalization theory. arXiv preprint arXiv:2006. 15733, 2020.
Lan, Y., Li, Z., Sun, J., and Xiang, black-box pde solver: When deep learning meets operator splitting. Journal of Computational Physics, 491: 123443 2023.
Y. Dosnet as a non-
Mattey, R. and Ghosh, S. A novel sequential method to train physics informed neural networks for allen cahn and cahn hilliard equations. Computer Methods in Applied Mechanics and Engineering, 390:1 14474, 2022.
Lanthaler, S., Mishra, S., and Karniadakis, G. E. Error estimates for deeponets: A deep learning framework in infinite dimensions. Transactions of Mathematics and Its Applications, 6(1 xttac0011, 2022.
Morton, K. W. and Mayers, D. F. Numerical solution of par- tial differential equations: an introduction. Cambridge university press, 2005.

---

Homotopy Dynamics
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E. Z., DeVito, Z., Raison, M., Tejani. A., Chilamkurthy, S., Steiner. B., Fang, I.., Bai, J., and Chintala, S. Py PyTorch: An Imperative Style, High-Performanm Deep Learning Library. arXiv 2019.
Yang, Y. and He, J. Deeper wider: A perspective from optimal generalization error with sobolev loss. Forty-first International Conference on Machine Learning, 2024,
Yang, Y. and Xiang, Y. Approximation of functionals by neural network without curse of dimensionality. J Mach Learn, 1 (4):342-372, 2022.
preprint arXiv:1912.011
Yang, Y., Wu, Y., Yang, H., and Xiang, Y. Nearly opti- mal approximation rates for deep super relu networks on sobolev spaces, arXiv preprint arXiv:2310.10760 2023.
Raissi, M., Perdikaris, P., and Karniadakis, G. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physies, 378:686-707, 2019.
Yang, Y., Chen, Q., and Hao, W. Homotopy relaxation train- ing algorithms for infinite- two-layer relu neural networks. Journal of Scientific Computing, 102(2):40, 2025.
Rathore, P.. Lei, W. Frangella, Z., Lu, L., and Udell, M. Challenges in training pinns: A loss landscape perspec- tive. arXiv preprint arXiv:2402.01866 2024.
Zhang. J. Zhang. S.. Shen, J., and Lin, G. Energy-dissipative evolutionary deep operator neural networks. Journal of Computational Physics, 498:112638, 2024. ISSN 0021-9991, doi: http://bbioogg/0o000......2223.3 URL https://ww..sceencddeeee.ccccc octenne/artiice///////2229999119
Siegel, J. W. and Xu, J. Approximation rates for neural net- works with general activation functions. Neural Networks, 128:313-321, 2020.
Sun, J., Dong, S., and Wang. F. Local randomized neural networks with discontinuous Galerkin methods for partial differential equations. Journal of Computational and Applied Mathematics, 445:1 15830, 2024.
Zhang, Z., Li, J., and Liu, B. Annealed adaptive im- portance sampling method in pinns for solving high dimensional partial differential equations. Journal of Computational Physics, 521:113561, 2025. ISSN 0021- 9991. doi: https://biiogg///000bo....222222 URL attene//arrade///////22922200000
Vershynin, R. High-dimensional probability: An introduc- tion with applications in data science, volume 47, Cam- bridge university press, 2018.
Wang, Multi-scale deep neural network (mscalednn) methods for oscillatory stokes flows in complex domains. Communications in Computational Physics, 28(5):
2020.
Wang. S., Wang, H., and Perdikaris, P.
On the eigenvec- bias of fourier feature networks: From regression to solving multi-scale pdes with physies-informed neural networks. Computer Methods in Applied Mechanics and Engineering.388 384: 13938, 2021.
Wight, C. and Zhao, J. Solving Allen-Cahn and Cahn- Hilliard equations using the adaptive physics informed neural networks. arXiv preprint 27X1vv2007.04542, 2020.
Xu, X. variational analysis for the moving finite element method for gradient flows. arXiv preprint arXiv:2009.01 01.0009..3993 393, 2020.
Xu, Z.-Q. J., Zhang, Y., and Luo, T. Overview frequency principle/spectral bias deep learning. Communica- tions on Applied Mathematics and Computation, pp- 1- 38, 2024.
Yang, Y. DeepONet for solving PDEs: Generalization analy- sis in Sobolev training. arXiv preprint arXiv:2410.03344, 2024.
12

---

Homotopy Dynamics
A. Proofs of Theorems 4.1 and 4.3
Before we prove Theorem 4.1, note that our method can be readily generalized to deep neural networks, as the underlying theoretical techniques remain the same; we merely need combine our approach with the results in (Du al., 2019). In this paper, we employ two-layer neural network simplify the notation and enhance readability, focusing on explaining why training becomes challenging when is small.
A.1. Amin (SS") >
In this subsection, we consider two-layer neural network defined as follows:
<(x; 0)
On a (w at),
(29)
k=1
where the activation function satisfies the following assumption:
Assumption A.1. The function o()) is analytic and not polynomial. Moreover, there exists positive constant such that
[x(x)] c/
for x.
Note that both (2: = In(1 + exp(x)) and (<x) I+exp(-x) satisfy Assumption A.1.
We assume that the weights and biases are sampled as follows:
N (O, an 1),
(30)
where (O0 1) denotes the standard Gaussian distribution.
The kernels characterizing the training dynamics take the following form:
kldl (x, x') (w a))o( a²)
klwl (x, x²) aa d d'ww" ar)o (w T x. x'.
(31)
The Gram matrices, denoted as k! and [w|, corresponding to an infinite-width two-layer network with the activation function a, can be expressed follows:
K. k (z; xz). Klo) (K!4) )nxn
K = xlw)( 2z), k(w) (K(W)) )nxn-
(32)
Lemma A.2 ((Du al., 2019)), Suppose that Assumption A.I holds and for any i E [n], * J. H æz Then the matrices K(w) and k are strictly positive, iee.
A = min (Ammi (K(wO). Amain (K(ca) > 0.
(33)
is easy to check that
k(w]__ Kl(a
lim SS
(34)
m=+00
based on the law of large numbers. Furthermore, we can show that the accuracy decreases exponentially as the width the neural network increases.
Definition A.3 ((Vershynin, 2018)). A random variable X is sub-exponential if and only its sub-exponential norm is finite i.e.
| := inf > 0 Ex/e/x//> 2.]
(35)
Furthermore, the chi-square random variable X sub-exponential random variable and Cy.d == |X ||
13

---

Homotopy Dynamics
Lemma A.4. Suppose that N (0,Ta) N(O,1) 1) and given w J. E Q. Then we have
(f) == (w x.) (x æz)), then |X S cdC di where is the constant shown in Assumption A.I. (ii) if' a'd' (w) a.) of (w ayg) my a.j. then 11 ad
Proof. (i) X < / dz and
inf {s > 0 Ex (xp([XX/s) - 2}
inf {s> > E exp o (w x. ( x;)l s 2]
inf {s
E exp (ed 1//s)) 2
inf (s
Ez exp Z1/8) 2)
idl (8 > 0 Ez exp (Z//8) 2
ed [d|||x ||x
cdCsvad
(ii) X & cdja|? z edZ and || S cdC
Proposition A.5 (sub-expon (sub-exponential Bernstein's inequality (Vershynin, 2018)), Suppose that XI...., ......,,,, Xm are i..dd sub- exponential random variables with EX1 μu, then for any 2 0 we have
x²
(-Commin
P
X
1 28
< 2exp
| 12
|
where Co is an absolute constant.
Proposition A.6. Suppose that Assumption A.1 holds, and given 5 E (0,1), w N (0, Ia), a a N(0, 1) and the sample
set {o C. Qt with æ, 's drawn i.i.d. with uniformly distributed with and any i.j E [m], + j, i wg. If
m 16n 2f 4n
log then with probability at least over the choice of (0), we have
Coxa
Amin (SST)
3
(mmin(Kll)) (K|al) Amin Amin(( (0)
Proof. Recall that
S = Vouo Vouo(an))
(36)
and 0 contain two parts, and parts, therefore SS can be rewrite S.S. + Sw where
S. = Vauo(x)),.... Vauo(an)), S.. = (0 V Vo (x)))
(37)
For any E 0, we define
s = (S.S. K(oo
< $
2.51 = {0l (Sw K|W=| I
(38)
n
Setting E S endCa.a. by Proposition A.5 and Lemma A.4, we have
mCge²
(- mCpe²
P(( 2exp (- "n'd'
P(N
-2exp
(39)
Due to inclusion-exclusian Principle, we have
P ((0) S.ST. (0) - k|al|| Se se) ST. - k(wl|
so}) 2 Σ
((2)))) P(()(a)) 2n² 1,
(40)
14

---

Homotopy Dynamics
therefore, with probability least
mC
1 - 4n exp
over the choice of 0, we have
-Klu S, -k)
||
(T) - Klw] SE
(41)
Hence by taking e 1 and 8 = 4n² exp
mCDdi
where A1 Amin (K(o)), Amin ( K (
16n-
Amin (SST)
2 Amin (Sa5%) Amin (SW
Σ Amin Klal) + Amin (K(W) (ww
S -Klo
SwS - -k.u
F
3
mmin (klu)) Amina (Klwl))
(42)
Combining Lemma A.2 and Proposition A.6, we obtain that under the conditions stated in Proposition A.6, the following holds with high probability:
Ami (SS >
(43)
A.2. Proof of Theorem 4.1
We can analysis the smallest eigenvalue of the problems based on the following lemma:
Lemma A.7 (LLi & Mathias, 1999)). Let A be an n x n Hermitian matrix and let À = AT. Then we have
Amin (A)
Amin (T*)) S
S Amaa (T*T).
Amin (A)
Proof of Theorem 4.1. We first show that Amin (k)) > 0, which follows directly from Lemma A.7:
Amio (Kg)
Amin (SS Amin (D.D.) > 0.
Therefore, the beginning of gradient descent, the kernel of the gradient descent step is strictly positive. We then define T
T it (tt) & (((())))
(44)
where
(()) = { (0())(((tt) K, (0(0))] S "Amin (min(Ks))
We now analyze the evolution of the loss function:
(L(()(()) = Vo d0
dt
dt D.D.S D.
n
2
Amai (K. (0())))(00)
(45)
n
where we use the fact that 1. 1- = (n
Furthermore, for E [0, T), we have
1
K, (0(t)) K. (((00))) S
2 Amiil (K)).
15

---

Homotopy Dynamics
This implies
Amin (K (O())))) 2 Amin
(x(O(t)) - (o(0)))) + Amin K (0(0))
Amin
(K(()))))
min (k( K (0()))
Amin
(x((()))) [k(((( - (0(0))
(K=(0(0)).
Amin
Therefore, we obtain
dL (0(1))
-- Amin (0))
(46)
di
fort [0, .. Solving this differential inequality yields
(o(tt) - Lu (0(0))) (- Amin (k))
(47)
for all t E [0. Finally, for the inequality
Amiin (SS M ( (D.D.) S Amia (K₂) S (Sn )m (SS")) (D.D.)
(48)
it follows directly from Lemma A.7.
A.3. Proof of Theorem 4.3
Proof of Theorem 4.3. First, we have
u(Ex+1) u(Ex) + (Ek+1 Ex))u -xx)x²
= =((Ex) + (Ek+1 Ex)H(ER, u(Ex)) (Ex+1-
(49)
where & lies between 1x+1 and Ex and depends Therefore, we obtain
c(Ex+1) = E(Ex) + (EK+1 - (x)(E( u(Ex)) N(EE, U(Ex))) + 5(ENt1 - uu")
(50)
where e(Ex) = u((x)) - U(E)). Then, we have
[1((E++))l
=1le(E))l(aF()) + (E++1 - Ex)] u(ex)) - M(,, ((E))))HHA))))
1 E0 Sn
<|le(E=)]l=<(()) + (Et+1 - ex) Pr, ([(())) (2))
2
<|lc(Ex)]l + P. Eo - |((Ex) (2())
1 1 En 7.
(51)
2
Recalling that eo = []u(50) (Eo) (( we obtain
12
[(((n))l se
1 + P.
Eg En
+
1 A.P..
Eo En
71
En En
+
Sa
-
1 1+p.
+
2
P
7(e
- 1)
ene P(Eo-En)
(52)
2P
16

---

Homotopy Dynamics
where the last step follows from the inequality
(1 a)m S ema
for a > 0.
Corollary A.8 (Convergence of Homotopy Functions). Suppose the assumptions Theorem 4.3 hold, and H (en; u) Lipschitz continuous in (f2) i.e.,
|| H (u),En)) - (u En)ll S L| - 42l u2llm
Then, we have
[H((((n).
L ege P(eosaa) + T(e - 1)
<1
(53)
21
Proof. The proof follows directly from the result in Theorem 4.3.
A.4. Discussion ep
In Theorem 4.3 and Corollary A.8, one important assumption is that assume ep is small. Here, we discuss why this assumption is reasonable.
First, we use physics- neural networks (PINNs) to solve the following equations:
L. Su), in Ω Bu g(x), on 8n
(54)
where C. is a differential operator defining the PDE with certain parameters, B is operator associated with the boundary and/or initial conditions, and Ω C Rd.
The corresponding continuum loss function is given by:
A / (Bu(æ; 0) g(x))² dæ. (.) = 1 [ (C.u(x; - f(u))² do + 2 2
1
(55)
We assume this loss function satisfies a regularity condition: Assumption A.9, Let u, the exact solution of Eq. (54). Then, there exists a constant C' such that
||u 0) u. (=)) S (L.(O).
(56)
The above assumption holds in many cases. For example, based (Grisvard, 2011), when L is a linear elliptic operator with smooth coefficients, and f (u) reduces to (x) T22) and if Ω is polygonal domain (e.g., [0. 1]4), then, provided the boundary conditions are always satisfied, the assumption holds.
Therefore, we only need to ensure that (0.) )) is sufficiently small, where 0, denotes the learned parameters at convergence. Here, (O.) can be divided into three sources of error: approximation error, generalization error, and training error:
/ (C.u(w; 0) - (((x)))) dx + / (Bu(x; 0) - g(x)))² dx,
0, = arg min (0)
arg min 0
0
0 arg min (()) arg min
1
(C.u(xt; ful(xt 0)))
1
(uulzf; - s(xx))
(57)
2nd
2nh
where ar.,, 2/, are sampled points defined in Eq. (2).
17

---

Homotopy Dynamics
The error decomposition can then be expressed as:
ELo(0,) S Lo(0.) + - (LLO)) EL(O). (L(O2) + (..) EL((.))
I Lc( +EL(O.) + EL, + EL(O,) - EL(Oa).
(58)
approximation ervor
generalization error
training
where the last inequality is due EL(O2) (.) < 0 based on the definition of Os-
The approximation error describes how closely the neural network approximates the exact solution of the PDEs. If f Lipschitz continuous function, Le is Lipschitz continuous from (S2) + L' (9t, and B is Lipschitz continuous from (02) + (822, with u(;; 0), u. E W2,00 (St) and On E (2)), then we have
L.(0)- / (C.u(x;; f(uxx))) (Ceu, f(u,))² dæ + À f (Bu(æ; 0) g(x)))² -BBu, - g(x))) dæ -)) u.)ll (If((u(x;; 0)) f(u,)) f(m)))lmm)) + C2 ((((xx;)) u. u.)ll.)(00))
on
s C1 (1(x:
$ Os]lu(<; 0) u. wlllw + Callu(x; 0) u. (.llw..(())
s <O)u((;;)) - u. (l))
(59)
where the second inequality follows from the trace theorem (Evans, 2022). Therefore, we conclude that (0) can be bounded by [|u( 0)- u. (o)) which has been widely studied the context of shallow neural networks (Siegel & Xu, 2020) and deep neural networks (Yang al., 2023). These results show that if the number of neurons sufficiently large, the error in this part becomes small.
For the generalization error, it arises from the fact that we have only a finite number of data points. This error can be bounded using Rademacher complexity (Yang et al., 2023; Luo & Yang, 2020), which leads to a bound of Tir other words, this error term small when the number of sample points large.
+
+ n
In
For the training error, Theorem 4.1 shows that when is large in certain PDEs, the loss function can decay efficiently, reducing the training error to small value.
B. Details on Experiments
B.1. Overall Experiments Settings
Examples. We conduct experiments on function learning case: ID Allen-Cahn equation, 2D Allen-Cahn equation, high dimension Helmholz equation, high frequency function approximation and operator learning for Burgers' equation. These equations have been studied in previous works investigating difficulties in solving numerically; use the formulations in Xu (2020); Zhang al. (2024); Hao & Yang (2019) for our experiments.
Network Structure. We use multilayer perceptrons (MLPs) with tanh activations and three hidden layers with width 30. We initialize these networks with the Xavier normal initialization (Glorot & Bengio, 2010) and all biases equal to zero.
Training. We use Adam to train the neural network and we tune the learning rate by a grid search on {10-5, 5, 10-4 10-3 10 10-2). All iterations continue until the loss stabilizes and no longer decreases significantly.
Device. We develop our experiments in PyTorch 1.12.1 (Paszke al., 2019) with Python 3...12. Each experiment is run on single NVIDIA 3070Ti GPU using CUDA 11.8. As summarized in Table 6, report the training cost (per epoch and total epochs) for experiments in our paper.
B.2. 1D Allen-Cahn Equation
Number of residual points Thess 200 and number of boundary points The 2. In this example, we use forward Euler method to numerically solve the homotopy dynamics. And = 0.1 and E,, = 0.01, here we choose Aeg- 0.001. Here, we use strategy to train the neural network.
The results for using original training for this example Figure 7. As shown in the figure, the original training method results in large training error, leading to poor accuracy.
18

---

Homotopy Dynamics
1.00
Ground truth
0.59
4.25
5.0
-.75
-1.00
00
02
0.4
0.0
0.8
Figure 7. Solution for ID Allen-Cahn equation for origin training.
Table Training Cost for Different Examples. The table summarizes training time per epoch and total number of epochs for various
PDE problems.
Example
ID Clen-CCahn Equation
Example 5,1
Example 5.2
Example 5.3
Training Time per Epoch Total Epochs (Steps)
0.05s
0.09s
0.0ls
0.4s 2.0 x 106
1.0 x 10s
4.0 x 106
4.0 x 106
B.3. 2D Allen-Cahn Equation
Number of residual points f4es 50 50 and number of boundary points Tiee 198. For fair comparison, all methods were implemented using the same neural network architecture, specifically fully connected network with layer sizes [2, 30, 30, 1]. In this example, we optimize using the Homotopy Loss. We set
1.0 and O. initially choosing As = 0.1, and later refining it to At = 0.01. When s = 0.05, {ss) = 0.05 we fix E = 0.05 and gradually decrease to O.
The reference ground truth solution is obtained using the finite difference method with 1000 x 1000 grid points. The result shown in Figure 8.
1.00
046
0.5
0.9
025
080
-00
4225
4330
0.75
+1.00
3.00
-035
-0.50
-0.25
0.0
025
050
0.5
1.00
-036
Figure Reference Solution for 2D Allen-Cahn equation.
The result obtained using PINN is shown in the Figure 9. It is evident that the solution still deviates significantly from the ground truth solution.
The result obtained using curriculum regularity strategy (Krishnapriyan et al., 2021) is shown in the Figure 10 below, where AE 0.01.
The result obtained using seateenncet---qqeenee training strategy (Wight & Zhao, 2020; Mattey & Ghosh, 2022) is shown in the Figure 11 below, where At = 0.01.
The result obtained using resampling strategy is shown in the Figure 12. In all resampling strategies, additional sample points are eventually concentrated near the sharp interface region. In our comparative experiments, we start with uniform
19

---

Homotopy Dynamics
-1.00
-.75
-1.00
-0.25
Figure 9. Solution for 2D Allen CaennCaahn equation for origin training.
Solution at T= 5, 5- 1.0
Solution at T= 5, 0.5
1.00
1.00
0.75
0.75
0.50
0.50
000020
0.25
0.25
-0,00075
0.00
0.00
-0 0.00110
0,00225
-00005
3.50
0,001
0.50
-9.0315
-0 0.0005
(1)
-1.00
-00093
-1.00
0.55
Solution at 7= 5, 0.1
Solution at 7= 5, 0.05
0.05
6.75
0.005
050
-U00/5
029
0.55
0075
8.00
-0.0255
0.00
0150
-0.0300
-0.75
3.50
-0.0375
0.50
-0.07 0.0000
-0090
-3.00
-00325
-1.00
0.0450
(3)
Figure 10. Numerical results for the 2D Allen Cahn equation using the Curriculum training strategy.
grid of 50 x 50 sample points and augment it by adding 5000 points in the vicinity of the sharp interface.
20

---

Homotopy Dynamics
Solution at T 0.5, 0.05
Solution at T 1.5, 0.05
1.00
100
1.00
0.25
0.50
0.5
3.0
0.25
025
325
0.5
0.25
0.25
0.2
0.25
0.50
0.50
-1.00
-1.00
+1.00
-1.00
(1)
-0.0
-0.2
0.75
1.00
-1.00
0.00
0.2
Solution at T= 2.5, 0.05
Solution T= 3.5, 0.05
1.50
1.00
1.00
1.00
0.50
0.00
350
0.25
025
325
0.25
0.00
25
-0.25
-325
-0.30
4.50
0..35
-0.7
-1.00
-1.00
-1.00
-1.00
1.00
0.73
0.50
0.25
0.00
0.23
0.50
0.75
1.00
(4)
1.00
00.00
0.25
0.00
0.55
0.50
0.75
1.00
Solution at T 4.5,
0.05
Solution
0.05
1.00
1.00
1.00
1.00
0.25
0.5
0.5
325
0.00
3.00
3.00
-0.55
-0.25
0.0
-0.50
0.50
-0.50
-075
-3.75
(5)
0.00
-0.75
-0.50
-0.25
0.00
0.5
0.50
0.75
1.40
(6)
-1.0
-0.75
-0.00
-0.55
0.00
0.50
0.7%
1.00
Figure Numerical results for Allen-Cahn equation using the Sequence-to-sequeree training strategy.
B.4. High Dimension Helmholtz Equation
Number of residual points 10000 and number of boundary points 2000. Neural network architecture is fully connected network with layer sizes [2, 20,00.30 30, 30, 1]. In this example, we optimize using the Homotopy Loss. We set Eo 1.0, initially choosing As = 0.1, and later refining it to As = 0.01 until sx, 30
Largest eigenvalue D.. As shown in Figure 13, a smaller e results in a smaller largest eigenvalue of (60), leading to a slower convergence rate and increased difficulty in training.
21

---

Homotopy Dynamics
Sample Points
Solution T 0.05
1.00
100
012
0.25
-.25
0.25
0.00
-050
50
0.00
0.75
-0.75
-0.09
(1)
-1.00
0.75
-0.50
-0.25
0.00
0.5
0.50
0.75
1.00
-1.00
-0.15
-0.20
-9.25
0.90
0.23
0.50
0.75
1.00
Figure Numerical results for the Allen Cahn equation using the Resampling training strategy.
Largest Eigenvalue of Dr
1e10
E = 0.1; Amax 3.74e6
3 = 0.5; Amax = 1.73e9
2.5
1.0; Amx 2.56e10
1.5
0.5
Figure 13. Largest eigenvalue of D. (60) for different e, A smaller e results in smaller largest eigenvalue of (60), leading to slower convergence rate and increased difficulty training.
D = - + diag (diag(1.....)) 1)
(60)
B.5. High Frequency Function Approximation
We aim to approximate the following function: s5n x [0, 1]. The homotopy is defined as (u, - (in Tz), where E E 50 1.1. Number of residual points nles = 300. In this example, we optimize using the Homotopy 1
1 and en = 50 the list for { is [15 20 25 30 35 40 45 Loss. We set cn = 15 the homotopy dynamics approach can also mitigate the slow training issue caused by the Frequency Principle (F-Principle) when neural networks approximate high-frequency functions.
1 J.l. From this example, we observe that
As shown in Figure 14, due the F-principle (Xu et al., 2024), training particularly challenging when approximating high-frequenc functions like sin(50...). The loss decreases slowly, resulting poor approximation performance. However, training based on homotopy dynamics significantly reduces the loss, leading to better approximation of high- -reeffreyency functions. This demonstrates that homotopy dynamics-bast bbsed training can effectively facilitate convergence when approx- imating high-frequency data. Additionally, we compare the loss for approximating functions with different frequencies 1/E using both methods. The results, presented in Table 7, show that the homotopy dynamics training method consistently performs well for high- -yeeurrey functions.
22

---

Homotopy Dynamics
Table Comparison of the lowest loss achieved by the classical training and homotopy dynamies for different values of in approximating sin mz)
3.29e-1 8 s 1/15
8 8=
1/35
8 = 1/50
Classical Loss
4.91e-6
7.21e-2
Homotopy Loss L
1.73e-6
1.91e-6
2.82e-5
Lows Curves
171
Evelutio
Origin Evalution
Homwtopy
Iteratons
Oc 1/15
1=21/
1/35
--- 1/50
Homotopy Evolution Evolutioo
Figure 14, High-frequency function sin(50 approximation: Comparison of loss curves between original evolution and homotopy evolution. The comparison shows that homotopy evolution effectively reduces the loss, successfully approximating the high- high-frequency function, while the original evolution fails. The number of residual points Thes 300.
B.6. Operator Learning 1D Burgers' Equation
In this example, we apply homotopy dynamics to operator learning. The neural network architecture follows the DeepONet structure:
P
n
Go(u)())
a
Σ ct (z,) + o (wx y + bx)
(61)
1=11
Here, a (we, y by) represents the trunk net, which takes the coordinates D' as input, and a (E u represents the branch net, which takes the discretion function as input. Rigorous error bounds for DeepONet are established in (Lanthaler et al., 2022; Liu al., 2024a; Yang. 2024), so we omit them here. We can interpret the trunk net as the basis functions for solving PDEs. For this example, the input is ug and the output is uoa- We still train using the homotopy loss. It is important emphasize that, unlike conventional operator learning, which typically follows supervised learning strategy, our approach adopts unsupervised learning paradigm. This makes the training process significantly more challenging. The initial condition uo(2) is generated from Gaussian random field with Riesz kernel, denoted by GRF ~ N' (0, (-A + 49 and A and represent the Laplacian and the identity. We utilize a spatial resolution of 128 grids to represent both the input and output functions.
We want find the steady state solution for this equation and 0.05. The homotopy is:
H(u,,,e)=(1 (1-s) (5
uss)unx (in() cos( s(u -w)).
(62)
where s E 1]. In particular, when s = 1, the initial condition un automatically satisfies and when becomes the steady state problem. And e(s) can be to
23

---

Homotopy Dynamics
Largest Eigenvalue of D.
1010
2.5
0.05; Amx 9..2e5
ε = 0.1; Amax = 5.42e6
2.5
0.5; Amax 1.76e9
155
0.5
0.0
0.2
0.4
0.6
1.0
Figure 15. Largest cigenvalue of D. (64) for different e. A smaller e results in smaller largest eigenvalue of (64), leading to slower convergence and increased difficulty in training.
Table 8. Comparison of accuracy and efficiency between Finite Difference Method (FDM) and DeepONet (trained via Homotopy
Dynamics).
At
Finite Difference Method (FDM)
DeepONet (trained by Homotopy)
L2RE
MSE (x.)
Comp. Time (s)
Loss LH
L2RE
MSE (x:,)
Inference Time (s)
0.5
5 x 10-5
1.63e-12
7.35e-13
239.98
7.55e-7
1.50e-3
1.75e-8
0.2
0.1
1 x 10-5
5.83e-4
1.57e-5
1239.77
3.40e-7
7.00e-4
9.14e-8
0.2
0.05
5 x 10 10-6
1.01e-2
4.20c-3
2416.23
7.770-7
2.52c-2
1.20c-3
0.2
8, SE [0.05, 1]. {
E(s) = 0.05 E [0, 0.05].
(63)
Here, s()) varies with during the first half of the evolution. Once es)) reaches 0.05, it is fixed at e(s) 0.05, and only s continues to evolve toward 0.
Largest eigenvalue of De. As shown in Figure 15, smaller results in smaller largest eigenvalue of (64), leading to slower convergence rate and increased difficulty in training.
du(x1)
du(xx)
D - + diag u(xx))
+
u(xn)
+
(64)
dx dis
dr
dr
dæ
Compare with tradition method.
The table reports both inference time and accuracy metrics across varying While FDM achieves high accuracy, its computational cost increases significantly for small due to CFL constraints. Moreover, its accuracy deteriorates under small possibly due to resolution limitations. In contrast, our DeepONet model yields substantially faster inference with only moderate accuracy degradation, making it well-suited for many-query scenarios such as uncertainty quantification or real-time control.
24


Reviews:
[{"review_id": "9CErU53ah7", "reviewer": "Reviewer_gjqL", "summary": "This paper proposes Homotopy Dynamics to train neural network for solving sharp interface problems. For sharp interface problems, the parameter $\\epsilon$ in the PDE affects the singularity of the solution. As $\\epsilon \\to 0$, the PDE becomes increasingly singular and thus the solution is difficult to compute. The authors first train the\nneural network on PDEs with a large $\\epsilon$, and then adjust the neural network according to the evolution of the homotopy dynamics until $\\epsilon$ decreases towards the target value. To validate the performance of the proposed methods, Numerical experiments on Allen-Cahn equation, high frequency function approximation, and Burgers equation are performed to validate the performance of the proposed method.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful reading and valuable suggestions. Our work addresses the core challenge of training neural networks to solve PDEs with sharp interfaces, where small parameters introduce near-singularities that hinder optimization. **From a theoretical perspective, we analyze how such parameters affect training convergence.** To overcome these difficulties, we propose a **homotopy dynamics-based training strategy** with rigorous convergence analysis. Experimentally, we show the method performs effectively on the 2D Allen–Cahn equation, mitigates spectral bias, and generalizes well to the operator learning setting.\n\nIn response to the concerns you raised, we provide the following answers:\n\n- **Motivation for using neural networks.** Neural networks are widely used as PDE solvers due to their strong approximation power[1] and ability to mitigate the curse of dimensionality [2]. They also benefit from automatic differentiation for efficient derivative computation. However, training neural networks for complex PDEs remains difficult. Our work targets these optimization challenges with a homotopy-based solution.\n\n  For time-dependent PDEs like Allen–Cahn, computing the steady-state solution ($T \\to \\infty$) with traditional solvers requires very small $\\Delta t$ as $\\epsilon \\to 0$, making them expensive. Neural networks, once trained, allow efficient inference—especially on GPUs. In Example 5.3, we add a comparison between our homotopy-trained DeepONet and classical finite difference methods across 1,000 equations (see [Table 1](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view?usp=sharing)), demonstrating significant computational speedups.\n\n  Neural network methods for these PDEs have gained attention recently[3,4,5]. Many physical models (e.g., phase-field dynamics) are governed by Allen–Cahn-type equations. Neural networks provide a path to integrate real experimental data into more accurate physical models.\n- **Clarification on Theorem 4.1.** In Theorem 4.1, we present two inequalities. The first inequality (Eq. (19)) characterizes the worst-case scenario, attaining equality when all components vanish except the one associated with the smallest eigenvalue. In practice, however, the decay of the other components is typically faster, rendering Eq. (19) nearly tight. An eigenvalue decomposition reveals that components along other eigen-directions diminish more rapidly, so that the smallest-eigenvalue component dominates after a short transient period.\n\n  For the second inequality (Eq.(20)), the Lidskii–Mirsky–Wielandt theorem provides the following full inequality:\n\n  $$\n  \\lambda_{\\text{min}}(D_\\varepsilon^{T} D_\\varepsilon)\\lambda_{\\text{min}}(S^{T} S) \\le \\lambda_{\\text{min}}(K_\\varepsilon) \\le \n  \\lambda_{\\text{max}}(D_\\varepsilon^{T} D_\\varepsilon)\\lambda_{\\text{min}}(S^{T} S)\n  $$\n\n  Thus, when $\\varepsilon$ is large, the training speed ranges from $\\exp(-Ct/n)$ to $\\exp(-Ctn^3)$. For small $\\varepsilon$, it consistently decays as $\\exp(-Ct/n)$.  The upper bound of the above inequality is attained in specific cases where a nonzero vector $x$ exists such that it is the eigenvector of the largest eigenvalue of $D_\\varepsilon^\\top D_\\varepsilon$, and $D_\\varepsilon x$ is the eigenvector of the smallest eigenvalue of $S^\\top S$. This occurs under certain $S$ and $D_\\varepsilon$ structures, depending on the PDE and sampling distribution. Hence, when $\\varepsilon$ is large, training can be relatively easy in some cases; however, when $\\varepsilon$ is small, training becomes universally difficult.\n- **Clarification on Example 5.3.** The objective in this operator learning task is to map the initial condition to its corresponding steady-state solution. In Example 5.3, we select a setup where the steady state is identical across initial conditions. This simplifies verification of whether the operator correctly maps diverse inputs to the same target. However, our homotopy-based training is not restricted to such cases and can be applied to settings where steady states differ. This example serves as a proof-of-concept showing the method's effectiveness even in an unsupervised operator learning setting, contrasting with the typical supervised approaches.\n\nFinally, we thank the reviewer for the constructive feedback and helpful suggestions, which have greatly improved the clarity and completeness of our paper.\n\n**Reference:**\n- [1][Yang et al., NeurIPS, 2023.](https://arxiv.org/abs/2305.08466)\n- [2][E et al., Constructive Approximation, 2022.](https://arxiv.org/abs/1906.08039)\n- [3][Wight et al., Commun. Comput. Phys., 2021.](https://arxiv.org/abs/2007.04542)\n- [4][Zhang et al., J. Comput. Phys., 2024.](https://doi.org/10.1016/j.jcp.2023.112638)\n- [5][Li et al., Comput. Methods Appl. Mech. Eng., 2023.](https://doi.org/10.1016/j.jcp.2025.113843)", "claims_and_evidence": "Most of the claims made in the submission are supported by clear and convincing evidence.\n\nQuestion:\n\n- In Theorem 4.1 the upper bound of $\\lambda_{\\min}(K_{\\epsilon})$ is given. In line 301, 'Consequently, the training speed can reach $\\exp(−Cn^3t$) based on Eq. (19), which is fast and implies that training is\neasy.' This statement cannot be derived from Theorem 4.1 unless the bound is tight. More explanation is needed.", "methods_and_evaluation_criteria": "The numerical experiments make sense for validating the performance of the proposed method.", "theoretical_claims": "I did not carefully check the correctness of proofs for theoretical claims in Appendix A.", "experimental_designs_or_analyses": "The experimental design appears sound and aligns well with the theoretical framework.", "relation_to_broader_scientific_literature": "The paper overcomes the optimization challenges lies in the training of neural network for solving PDEs; as introduced in Section 1. However, the motvation of using neural network to solve these PDEs are not adequately discussed. To my knowledge, the numerical solution of the example problems are well-studied using method other than neural network, to name just a few literature: some references in the paper including [Kreiss & Kreiss, 1986], [Hao & Yang, 2019] and\n\nKim, Yongho, Gilnam Ryu, and Yongho Choi. \"Fast and accurate numerical solution of Allen–Cahn equation.\" Mathematical Problems in Engineering 2021, no. 1 (2021): 5263989.\n\nShen, Jie, and Xiaofeng Yang. \"Numerical approximations of allen-cahn and cahn-hilliard equations.\" Discrete Contin. Dyn. Syst 28, no. 4 (2010): 1669-1691.\n\nJiwari, Ram. \"A hybrid numerical scheme for the numerical solution of the Burgers’ equation.\" Computer Physics Communications 188 (2015): 59-67.", "essential_references_not_discussed": "The paper adequately discusses the key related works necessary for understanding the context and its contributions.", "other_strengths_and_weaknesses": "Strengths:\n\n- The proposed method improve the accuracy of the solution compared with other neural network-based methods.\n\nWeaknesses:\n\n- The motivation of using neural network to solve this problem is not convincing enough. \n\n- In line 40. 'Leveraging neural network architectures to solve PDEs,... particularly in handling complex domains and incorporating empirical data'. Both these aspects were not emphasized throughout the paper.", "other_comments_or_suggestions": "Some typos:\n\n- In line 94 right column, 'represent represent'\n\n- In line 182 left column, the indices of $u$ are not consistent.", "questions_for_authors": "- As the neural network is designed and trained for a specific PDE, I think comparison with traditional methods should be included. For example, the computation cost (training time for using neural network versus the computation time using numerical solvers to solve the nonlinear equation); the accuracy of the solution, etc... to highlight the advantage of the proposed method.\n\n- In example 5.3, why not directly train and solve for the steady state solution since it is independent of the initial condition.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jvYMZbQbLJ", "reviewer": "Reviewer_Vjtr", "summary": "The authors look at the physics-informed neural network (PINNs) setting of solving a PDE via minimizing the PDE residual. They look at cases where there are “sharp” interfaces (introducing near singularities). They propose a method based on homotopy dynamics, which involves starting with an easier to learn problem and then moving towards a harder to learn problem (where “easy” and “hard” are characterized by the parameters in the PDE). They show this on the 1D Burgers equation and 1D and 2D Allen-Cahn equations.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful reading and valuable suggestions.\n\nFirst, we would like to emphasize the main motivation and contribution of our work.  **From a theoretical perspective, we analyze how such parameters affect training convergence speed**. To overcome these difficulties, we propose a **homotopy dynamics-based training strategy** with rigorous convergence analysis.\n\nBelow, we address the questions and concerns you have raised.\n\n- **Comparision to Related Work**\n\n  We would like to emphasize that, to the best of our knowledge, **our work is the first to provide a theoretical justification that in sharp interface problems**, the parameter $\\epsilon$ directly determines training difficulty—the smaller the $\\epsilon$, the harder the optimization (see **Theorem 4.1**).\n\n  While our approach shares a high-level idea with curriculum-based methods (progressing from easy to hard tasks), it differs significantly in design. Unlike [1], which lacks a systematic mechanism, our homotopy dynamics defines a continuous path in the PDE parameter space with convergence guarantees. Moreover, we provide a dynamical update rule and a theoretically grounded strategy for choosing the homotopy step size $\\Delta\\epsilon$ (see **Theorem 4.3**), which is not addressed in [1].\n\n  Regarding Fourier feature methods [2], their success depends on prior knowledge and sensitive tuning of the feature scale $\\sigma$. Our focus is on sharp interface problems, which differ from general multiscale settings. Example 5.2 shows that our method can generalize beyond sharp interfaces, demonstrating its versatility.\n\n  As for resampling-based methods [3,4], they often require large sample sizes and careful tuning, making them computationally expensive. In contrast, our approach achieves competitive accuracy with fewer collocation points and lower computational cost.\n\n  We further highlight these strengths in Example 5.1 (2D Allen–Cahn, $\\epsilon = 0.05$). Unlike [1], which needs 50 time steps ($\\Delta t = 0.1$), our homotopy strategy reaches the steady state in only 10 steps ($\\Delta s = 0.1$) and uses just 2,500 collocation points. This demonstrates both the efficiency and effectiveness of our method.\n\n\n  Finally, we have added additional experimental comparisons, which further support the advantages of homotopy-based training and [Results](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view) show that our homotopy-based method consistently achieves better accuracy.\n\n - **Theoretical Clarifications**\n\n    We appreciate the reviewer’s comments regarding theoretical assumptions and will incorporate further details in the revised version.\n\n1. **Network Depth and Generality**: While we present the theory using two-layer networks for simplicity, our framework can be readily extended to deep architectures, building on standard results such as those in [2]. Our theoretical analysis focuses on how a small $\\varepsilon$ induces optimization difficulties, regardless of the network depth, and the results remain valid. We use shallow networks solely to simplify the notation and enhance readability, thereby helping readers grasp our key points.\n\n2. **Activation Functions**: Although we use ReLU in our theoretical analysis, the results hold for other smooth activations such as $\\tanh$ and $\\text{ReLU}^k$ (except for Lemma A.1, which has a known analog in [2]). In the revision, we will clarify that our results are **not restricted to ReLU**, and we will present a more general analysis accordingly.\n\n3. **Clarification of Eq. (31)**: We note that Eq. (31) is **not an assumption**, but rather defines the **continuous kernel limit** as $m \\to \\infty$ (based on Law of Large Numbers). We will provide further explanation in the appendix to clarify this point.\n\n- **Other Questions**\n\n- Based on Example 5.1, we added comparisons with other methods. Results in [Table 3 and Figures 1–3](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view) show that our homotopy-based method consistently achieves better accuracy.\n\n- The training time and inference time for the numerical experiments can be found in  [Table 1 and 2](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view).\n\n- Our setting includes not only single PDEs, but also unsupervised operator learning, which is harder than the commonly studied supervised setup. Prior works [3,4] have explored this, but our method achieves higher accuracy. We believe our homotopy strategy can be extended to even more complex systems in future work.\n\n\n**References**\n\n[1] [Krishnapriyan et al., NeurIPS, 2021.](https://arxiv.org/abs/2109.01050)  \n[2] [Wang et al., Comput. Methods Appl. Mech. Eng., 2021.](https://www.sciencedirect.com/science/article/pii/S0045782521002759) \\\n[3] [Zhang et al., J. Comput. Phys., 2024.](https://doi.org/10.1016/j.jcp.2023.112638) \\\n[4] [Li et al., Comput. Methods Appl. Mech. Eng., 2023.](https://doi.org/10.1016/j.jcp.2025.113843)", "claims_and_evidence": "The authors claim that this method makes it easier to get better error on PDE problems with near singular features. They show this on 1D Burgers, and 1D and 2D Allen-Cahn equations and that this type of training (easy to hard parameters) gets better error than directly using the PINNs approach on the “hard” problem right away.", "methods_and_evaluation_criteria": "The method is based on the parameters of the PDE. They start with cases when the parameter in these different PDEs is higher, and the solution is easier to compute here. Then they train the model by going from this high parameter to the low parameter (and the original target problem).\n\nHowever, these PDEs are still quite easy, and the PINNs literature has come a long way since. These systems were being studied years ago, with similar errors, and the field should be moving to harder problems at this point. Additionally, there are a number of other methods for improving PINNs training and prediction (including very similar ones) that are not compared to at all. Two examples:\n\nThis method looks very similar to the curriculum regularization approach in [1], where the authors started with an easy to learn PDE problem and then slowly trained the model to solve the harder problem. The authors try to demonstrate their approach on high frequency problems, but many approaches to address high frequency function approximation already exist, a simple one is to add Fourier features [2]. Additionally, if looking at the operator learning setting, how would it look to use something like Fourier Neural Operator [3] instead? These papers are many years old at this point, and the field has progressed a lot since: new approaches should be looking at much harder problems, comparing appropriately to prior approaches, and going beyond methods that have already been proposed. \n\nThe authors also mention resampling: there are also adaptive weighting methods that sample places where the PDE residual is high [4]. However, it seems like this approach is also not that computationally cheap since it requires training the network for longer by starting with the easier-to-learn parameter and then going to the harder one.\n\n[1] Krishnapriyan et al. Characterizing possible failure modes of physics-informed neural networks. NeurIPS (2021)\n[2] Wang, Wang, Perdikaris. On the eigenvector bias of Fourier feature networks. CMAME (2021)\n[3] Li et al. Fourier Neural Operator. ICLR (2021)\n[4] C. Wu, M. Zhu, Q. Tan, Y. Martha, L. Lu. A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. CMAME (2023).", "theoretical_claims": "The authors present a theoretical analysis around the effectiveness of training PINNs via this homotopy dynamics approach. There are various places where a lot more steps would be helpful. There are also a lot of assumptions that are being made in these claims, such as only analyzing 2 layer neural networks, and looking at the width of the NN with ReLU activation functions. Does analysis based on ReLU activation functions actually apply here given that the neural network needs to be continuously differentiable to get derivatives and train with the PDE residual?\n\nAlso, how do you make the assumption in Eqn 31 in appendix?", "experimental_designs_or_analyses": "The authors set up the three different PDE problems. They train a PINN with and without using the homotopy dynamics approach (directly solving the problem vs going from easy to hard parameters). See the above comment that there are now many methods to train PINNs better, many of which have addressed similar problems that the authors are looking at. At this point, it is needed to show proof-of-concept on much more difficult problems, such as those that many current PINNs methods struggle with.", "relation_to_broader_scientific_literature": "There is a lot of work on PINNs and better training methods for PINNs, as well as a wide range of work on using ML to solve PDEs. This work needs to be better contextualized against this broad landscape, and a lot of the progress that has been made in the field.", "essential_references_not_discussed": "Work that proposes a similar idea, as well as other works that attempt to deal with the same problems (such as Fourier features for high-frequency learning), are discussed above. There are many off-shoots and follow-ups of these works that are also relevant.", "other_strengths_and_weaknesses": "See above for comments. The primary comments are that this work proposes something very similar to past work in the PINNs literature, doesn’t compare or contextualize against a vast literature of PINNs work that has been done to address many of the problems described here, and the experiment problems shown are relatively easy (compared to how far the field has come since).", "other_comments_or_suggestions": "See above for comments.", "questions_for_authors": "- How does this compare to adding Fourier features? How about the many other approaches that have been done for PINNs, such as the very similar curriculum regularization or adaptive sampling / adaptive weighting of the loss function? \n\n- What is the computational cost of this method?\n\n- These systems have been well-studied and are easy. Can these methods show proof-of-concept on much harder systems?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ta9oxWNgZ7", "reviewer": "Reviewer_m8xv", "summary": "In this paper, the authors present a training based method based on homotopy dynamics for handling sharp interface problems. The authors provide a proof of the convergence of the Homotopy dynamics for stable training. The experiments results demonstrate that the proposed method can help capture the sharp interfaces as well as approximating high-frequency functions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable suggestions on our paper. \n\nFirst, We would like to emphasize the main focus and contribution of our work. This paper addresses the core challenge of training neural networks to solve PDEs, particularly those involving sharp interfaces, where specific parameters in the PDE induce near-singularities and hinder optimization. **From a theoretical perspective, we provide a novel analysis of how such parameters affect the convergence behavior during training.** To overcome these difficulties, we propose a **homotopy dynamics-based training strategy** and rigorously establish its convergence properties. On the experimental side, we demonstrate that our method not only performs effectively on the 2D Allen–Cahn equation, but also alleviates the spectral bias commonly seen in neural network training. Furthermore, we show that this homotopy-based approach generalizes well to the operator learning setting, highlighting its versatility and broad applicability.\n\nBelow, we provide our responses to the questions and concerns you have raised.\n\n- We appreciate the reviewer’s suggestion regarding the extension of our method to larger-scale 3D cases. While the current examples in the paper are primarily 1D or 2D, it is important to emphasize that Example 5.3 is set in the operator learning framework, which is inherently more complex and challenging than standard PDE regression tasks.\n\n  In particular, unlike most existing operator learning methods that are trained in a supervised manner (i.e., with access to input-output solution pairs), our setting adopts a fully unsupervised training strategy based on homotopy dynamics, where the model learns the solution operator solely from the PDE structure. This significantly increases the difficulty of the learning problem.\n\n  Despite this challenge, our method achieves competitive and accurate results, highlighting its potential applicability not only to standard PINNs but also to more complex, unsupervised operator learning tasks.\n\n  We thank the reviewer again for the helpful suggestion, and in the revised version of the paper, we will incorporate higher-dimensional (e.g., 3D) examples to further demonstrate the effectiveness of our approach.\n\n\n- Regarding **training time**, we provide additional details here. All experiments were conducted on a single RTX 3070 Ti GPU. The corresponding computation times for training each epoch are summarized in the [Table 2](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view) \n    The detailed settings of each numerical experiment, including all parameter choices, are provided in Appendix B. We would like to emphasize that although our training time may appear relatively longer and the training procedure more involved, it enables us to achieve significantly higher accuracy. This level of precision cannot be attained by other methods, regardless of how long they are trained.\n    \n  Regarding **inference time**, we take Example 5.3 as a representative case in the operator learning setting. Specifically, we compare the inference efficiency of our trained DeepONet model with that of the traditional finite difference method, by solving 1,000 instances of the PDE using both approaches. The results are summarized in the [Table 1](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view).\n\n  As with other neural network-based methods, increasing the number of network parameters generally leads to longer training times. Thank you for your helpful suggestion — we will include these details and clarifications in the revised version of the paper.\n\n- The question you raised is very interesting. We believe that our proposed homotopy dynamics-based approach can be extended to cases involving multiple parameters. At present, our initial idea is that in the multi-parameter setting, the homotopy dynamics may need to update the parameters sequentially or in a coordinated manner. We consider this a promising direction and plan to explore it further as part of our future work.\n\nFinally, we would like to thank you again for your insightful comments and questions.", "claims_and_evidence": "The authors claim that the proposed method can improve the training process for sharp interface problems, high frequency function approximation and operator learning. These claims are validated by three examples in the experiments part. The authors also provide a theoretical proof for the convergence of Homotopy Dynamics for stable training.", "methods_and_evaluation_criteria": "The proposed methods are evaluated on three examples: 2D Allen Cahn Equation, High Frequency Function Approximation and Burgers Equation, which looks sensible to me.", "theoretical_claims": "While I have thoroughly reviewed the methodology presented in the paper, I did not perform an exhaustive line-by-line verification of all mathematical derivations and proofs.", "experimental_designs_or_analyses": "The experiments demonstrate the effectiveness of the proposed method on sharp interface problem. Here the problems have more shaper interface with smaller parameter **epsilon**. The experiments show that the proposed homotopy loss keeps low while classical loss increases dramatically as the parameter value decreases.", "relation_to_broader_scientific_literature": "The proposed method provides a training strategy for learning based method such as PINN, Operator learning to learn sharp interface problems.  This paper also present some intuitions on the training difficulties caused by certain parameters in PDEs learning.", "essential_references_not_discussed": "n/a", "other_strengths_and_weaknesses": "The paper is well-written and easy to follow. The detailed background materials provided helps enhance the reader’s comprehension of the paper.", "other_comments_or_suggestions": "n/a", "questions_for_authors": "- The examples shown in this paper are restricted in small scale 1D or 2D scenarios. Just wondering how would the method performed when applied to larger 3D cases? \n- What are the training/inference time costs for the proposed method? How sensitive are the time costs to the number of trainable parameters?\n- In this paper, the authors investigate the scenarios where a single parameter $\\epsilon$ quantifies how singular the system is. Does the proof still holds if there are multiple parameters? How would the method extend to handle these multi parameter scenarios?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "9VBPTcldM4", "reviewer": "Reviewer_fKP4", "summary": "This paper introduces homotopy dynamics as a strategy to solve PDEs with sharp interfaces using PINNs. The key idea is to start training with a larger interface width parameter $\\epsilon$ (corresponding to a smoother solution), then gradually decrease $\\epsilon$ to the desired sharp-interface regime. This approach is particularly relevant for PDEs such as the Allen-Cahn equation, where $\\epsilon$ controls the interface sharpness.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review and constructive suggestions. Below, we address the concerns you raised.\n- **Why do we use PINN for this type of problem instead of well-established numerical methods?**  \n  Neural networks offer strong approximation capabilities [1] and help mitigate the curse of dimensionality [2], making them well-suited for solving PDEs. They have been widely applied in this context [3,4], particularly in operator learning, where they can significantly accelerate computation. Moreover, many AI for Science models are governed by similar equations, and neural networks enable seamless integration of experimental data into the modeling process.\n\n- **Regarding the homotopy loss $L_{H_\\varepsilon}$ (line 181):**  \n  Without $L_{H_\\varepsilon}$, simply varying $\\varepsilon$ leads to instability and sensitivity to $\\Delta \\varepsilon$. The homotopy loss enables a stable Euler-forward-like path consistent with the homotopy dynamics. We will include supporting experiments in the revised version.\n\n- **What happens if we use only $L_{H_\\varepsilon}$ and omit the residual loss $L_H$?**  \n  A solution satisfying $H(u, \\varepsilon) = \\text{Const} \\neq 0$ still minimizes $L_{H_\\varepsilon}$. Thus the residual loss is necessary to enforce the target PDE constraint.\n\n- **What is the effect of $\\alpha$?**  \n  It balances the three loss terms, playing a similar role as $\\lambda$.\n\n- **Algorithm 1, line 199: Why do we need $\\Delta \\epsilon_k$ in Strategy 2?**  \n  It ensures proper scaling of $L_{H_\\varepsilon}$ and guides the homotopy step correctly.\n\n- **For the numerical examples, which strategy is used?**  \n  Strategy 1 is used for 1D Allen–Cahn and high-frequency examples. Strategy 2 is used for 2D Allen–Cahn and Burgers, where solving the linear system is more difficult.\n---\n**Theoretical Answer:**\n\n- **Items 1–3:** We agree. $L$ should be $L_H$, $l$ as $l_\\varepsilon$, and $K_\\varepsilon$ refers to $K_\\varepsilon(\\theta(0))$. These will be corrected.\n\n- **Item 4:** Based on Weyl’s inequality:\n\n  $\\lambda_{\\min}(A+B) \\ge \\lambda_{\\min}(A) + \\lambda_{\\min}(B)$\n  Therefore,\n  $\\lambda_{\\min}(K_\\varepsilon(\\theta(t))) \n  \\ge \\lambda_{\\min}(K_\\varepsilon(\\theta(t)) - K_\\varepsilon(\\theta(0))) + \\lambda_{\\min}(K_\\varepsilon(\\theta(0)))$\n\n  $\\ge \\lambda_{\\min}(K_\\varepsilon(\\theta(0))) - \\sigma_{\\max}(K_\\varepsilon(\\theta(t)) - K_\\varepsilon(\\theta(0)))$\n\n  $\\ge \\lambda_{\\min}(K_\\varepsilon(\\theta(0))) - ||K_\\varepsilon(\\theta(t)) - K_\\varepsilon(\\theta(0))||_F$\n\n  $\\ge \\frac{1}{2}\\lambda_{\\min}(K_\\varepsilon(\\theta(0)))$\nWe will add this in the revised proof.\n\n- **Item 5:** A large learning rate won’t solve the slow convergence issue for small $\\varepsilon$, and may lead to instability [3].\n\n- **Item 6:** Thank you—we will update the notation.\n\n- **Item 7:** For small $\\varepsilon$, $K_\\varepsilon$ may be large but remains finite. Using small $\\Delta \\varepsilon$ ensures stable error control.\n\n- **Items 8–9:** The index should be lowercase $n$, and the correct reference is [4].\n\n- **Item 10:** Without $L_H$ and $\\lambda L_{bc}$, Strategy 2 becomes a way to solve Eq. (7). However, due to small singular values in $H_u \\nabla_\\theta u$, solving directly may be unstable. Thus, we adopt optimization instead. Strategy 2 shares the same dynamics as Strategy 1. Including $L_H$ ensures the solution satisfies the PDE. Even if $H(u, \\varepsilon) = \\text{Const} \\neq 0$, it still minimizes $L_{H_\\varepsilon}$, hence both terms are necessary.\n\n---\n\n**Comments on the Weaknesses:**  \nWe thank the reviewer for the insightful comments. Compared with [5], our method introduces a homotopy loss that prevents convergence to bifurcation solutions. This key idea is not present in [5]. Also, the difficulty of training when $\\varepsilon \\to 0$ is not theoretically addressed in [6]; we believe ours is the first rigorous analysis. Strategy 2 is an alternative solver within the same framework, so our theoretical results are based on Strategy 1.\n\n**Experimental Comparison:**  \nBased on Example 5.1, we added comparisons with other methods. In this example, the homotopy is defined as $H(u,s,\\epsilon(s))$. Results in [Table 3 and Figures 1–3](https://drive.google.com/file/d/1fvQdXIZpm0c3fci_85Sb4CW-ZIjVHBlx/view) show that our homotopy-based method consistently achieves better accuracy. \n\n---\n**References:**\n\n[1] [Yang et al., NeurIPS, 2023.](https://proceedings.neurips.cc/paper_files/paper/2023/hash/449a016a6ce6fba3fe50d05482abf836-Abstract-Conference.html)  \n[2] [E et al., Constructive Approximation, 2022.](https://arxiv.org/abs/1906.08039)  \n[3] [Sutskever et al., ICML, 2013.](https://proceedings.mlr.press/v28/sutskever13.html)  \n[4] [Atkinson et al., Wiley, 2009.](https://homepage.math.uiowa.edu/~atkinson/papers/NAODE_Book.pdf)  \n[5] [Krishnapriyan et al., NeurIPS, 2021.](https://arxiv.org/abs/2109.01050)  \n[6] [Wang et al., Comput. Methods Appl. Mech. Eng., 2021.](https://www.sciencedirect.com/science/article/pii/S0045782521002759)", "claims_and_evidence": "The claims are supported by the provided evidence, but the support is not particularly strong, see Experimental Designs and Analyses.", "methods_and_evaluation_criteria": "The proposed method and evaluation criteria are reasonable for the problem at hand. But there are alternative methods that are not explored or discussed:\n\n- The bigger questions is, why do we use PINN for this type of problems instead of well-established numerical methods?\n- Regarding the the homotopy loss $L_{H_\\epsilon}$  (line 181). What happens if we don't use $L_{H_\\epsilon}$? That is, we only change $\\epsilon$. This would be similar to the curriculum regularization (Krishnapriyan et al. 2021).\n- Since $L_{H_\\epsilon}$ is supposed to mimic the forward Euler step in strategy 1, what happens if we don't use the residual loss $L_H$ and only use $L_{H_\\epsilon}$?\n- What is the effect of $\\alpha$?", "theoretical_claims": "A few clarifications would improve clarity\n- Equation (13) - (18), L should be $L_H$ in equation (12)? That is, the loss under analysis is the residual loss without homotopy dynamics?\n- Shouldn't $l$ depends on $\\epsilon$? That is, $l = l_\\epsilon$?\n- In the proof Appendix A.2, line 720 and 736, What is $K_\\epsilon$? Is it $K_\\epsilon(\\theta(0))$?\n- How do we go from 736 to 740? Is it Weyl's inequality? Please clarify.\n- From theorem 4.1 and the discussion afterward, if the difficulty only lies in the speed of convergence, will a larger learning rate help for original training?\n- Theorem 4.3, what is $K_\\epsilon$? Is it related to the $K_\\epsilon$ in Theorem 4.1? If not, using a different notation might improve clarity.\n- As $\\epsilon\\rightarrow 0$, the solution becomes more singular. Does $K_\\epsilon$ remain bounded? If not, then K might not exist (line 305 column 2).\n- What is $N$ in Theorem 4.3?\n- Appendix A.3, it seems (45) - (48) are more like standard analysis on Euler's method? How is it related to Antonakopoulos\net al., 2022?\n- Theorem 4.3 is only related to strategy 1, where Euler's method is used to evolve $\\theta$, but not strategy 2. The author should explain the connection between the two strategies.", "experimental_designs_or_analyses": "The experimental setup is reasonable, but the baseline comparison is relatively weak.\n\n- Original training are known to struggle with PDEs that have highly oscillatory or near-singular solutions. Prior works have proposed strategies to address this issue:\n  - **Curriculum regularization**: Krishnapriyan et al. (2021) suggest first learning simpler problems before tackling harder ones, similar to the homotopy dynamics approach. They also propose learning early time steps first, which relates to the design of the homotopy in Examples 5.1 and 5.3 (where $s = 1$ corresponds to the initial condition and $s = 0$ to the full time problem).\n  - **Neural Tangent Kernel Perspective**: Wang et al. (2021) analyze PINN training and show that PINNs are biased toward smooth solutions, which could be relevant to Theorem 4.1. The proposed remedy involves using random Fourier features.\n\nThe homotopy dynamics approach appears to provide a more principled way to implement some of these intuitions. However, a discussion on the similarities and differences between these methods, along with a performance comparison, would strengthen the paper.\n\n1. Krishnapriyan, A.S., Gholami, A., Zhe, S., Kirby, R.M., Mahoney, M.W., 2021. Characterizing possible failure modes in physics-informed neural networks.\n2. Wang, S., Wang, H., Perdikaris, P., 2021. On the eigenvector bias of Fourier feature networks: From regression to solving multi-scale PDEs with physics-informed neural networks.\n\nIn addition, there is also a difference between the experiments and the proposed method. When describing the proposed method, only $\\epsilon$ changes. However, for the experiments, another parameter $s$ is introduced, which also modify the residual loss.", "relation_to_broader_scientific_literature": "The proposed method improve PINN for a specific type of PDE problem.", "essential_references_not_discussed": "See Experimental Designs Or Analyses.", "other_strengths_and_weaknesses": "Strengths\n- The homotopy dynamics method is a principled and effective way to address challenges in training PINNs for sharp-interface PDEs. \n- The paper is well-written and easy to follow, making it accessible to readers with different levels of familiarity with PINNs and interface PDE\n- Theorems provide insights into why homotopy dynamics improves training. \n- The approach has the potential to benefit researchers working on PINN, particularly in problems involving sharp interface.\n\nWeakness:\n- The baseline is relatively week, see Experimental Designs Or Analyses. \n- Alternative strategies are not fully explored, see Methods And Evaluation Criteria.", "other_comments_or_suggestions": "- Algorithm 1, line 199, why do wee need $\\Delta \\epsilon_k$ in the inner loop of strategy 2\n- For the numerical examples, which strategy is used?", "questions_for_authors": "See previous comments, Methods And Evaluation Criteria, Experiment Design and Analysis.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper proposes an interesting approach that incorporates $\epsilon$, representing how sharp the phase transition in phase-field PDEs, into the training algorithm of a neural network approximator of the solution to phase-field problems. By the average score this is a typical borderline paper. The authors responded in the rebuttal with lots of extra experiments, however, I think this is an indication that this paper is not ready, or at least not ready for the targeted audiences. Here are a summary of questions or points I found in the reviews (including the comments in the AC-reviewer discussion) that would be helpful for a revision.

- How users should choose between the two strategies given in Algorithm 1? Is there any comparison?
- What is the strategy for choosing or designing the homotopy schedule? or any comparison or ablation on this matter.
- Some reviewers think that the experiments featured in this paper are "too easy", while I personally strongly disagree with this view, as in good science the experiments are ought to be illustrative of the theory. With this being said of my personal opinion, I respect the peer review process, and I would argue that adding baselines from previous PINN papers studying Allen-Cahn or Cahn-Hillard, such as [1] and [2], would make the case stronger.

Aside from reviewers' comments, I read the paper in details myself as well, here are my main concerns:
- Figure 1 definitely has some eye-catching drawing that is quite illustrative, however, as indicated by the figure, calling the SGD directly would result in "stuck in local minima", there is not a detailed ablation on this matter comparing the new training with a baseline training. In Figure 7, the caption reads "all final test solutions converge to the correct steady-state solution", yet does the vanilla SGD applying on the residual loss converge to an "incorrect" solution? This question remains unanswered throughout the paper. In Figure 2 there is a comparison of the loss curves of different epsilons, however, does a vanilla training strategy converge to an "incorrect" solution? or the large error is attributed simply to the solution being increasingly unsmooth?
- There is not an inner loop in the training algorithm, the $u_{\epsilon(s)} \to u_{\epsilon}$ is handled automatically by the homotopy dynamics. This is already vastly superior to most tailored training for PINN I read. However, the authors did not highlight this advantage over previous PINN training literature, for example [3].

Overall I think the idea is very good, and I have never seen anyone done this before (despite some of the reviewers' criticism). I believe this paper will be a very good contribution if the comments above can be added or addressed in a revision.


Minor typos:
- Page 2: "Homotopy dynamic"->"homotopy dynamics", and "the" should be capitalized after that.
- Page 6: "dynamic system"-> "dynamical system".
- Page 7: "a L2 error" -> "an L2 error" 

---
[1]: C. Wight, J. Zhao. Solving Allen-Cahn and Cahn-Hilliard Equations using the Adaptive Physics Informed Neural Networks. 

[2]: R. Mattey, S. Ghosh. A novel sequential method to train physics informed neural networks for Allen Cahn and Cahn Hilliard equations. CMAME (2022)

[3]: S. Wang et al. Understanding and mitigating gradient pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing (2021)
Decision: Accept (poster)