Learn Singularly Perturbed Solutions via Homotopy Dynamics

Thank you for your thoughtful review and constructive suggestions. Below, we address the concerns you raised.

• Why do we use PINN for this type of problem instead of well-established numerical methods?
  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.

• Regarding the homotopy loss $L_{H_\varepsilon}$ (line 181):
  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.

• What happens if we use only $L_{H_\varepsilon}$ and omit the residual loss $L_H$?
  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.

• What is the effect of $\alpha$?
  It balances the three loss terms, playing a similar role as $\lambda$.

• Algorithm 1, line 199: Why do we need $\Delta \epsilon_k$ in Strategy 2?
  It ensures proper scaling of $L_{H_\varepsilon}$ and guides the homotopy step correctly.

• For the numerical examples, which strategy is used?
  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.

---

Theoretical Answer:

• 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.

• Item 4: Based on Weyl’s inequality:

  $\lambda_{\min}(A+B) \ge \lambda_{\min}(A) + \lambda_{\min}(B)$ Therefore, $\lambda_{\min}(K_\varepsilon(\theta(t))) \ge \lambda_{\min}(K_\varepsilon(\theta(t)) - K_\varepsilon(\theta(0))) + \lambda_{\min}(K_\varepsilon(\theta(0)))$

  $\ge \lambda_{\min}(K_\varepsilon(\theta(0))) - \sigma_{\max}(K_\varepsilon(\theta(t)) - K_\varepsilon(\theta(0)))$

  $\ge \lambda_{\min}(K_\varepsilon(\theta(0))) - ||K_\varepsilon(\theta(t)) - K_\varepsilon(\theta(0))||_F$

  $\ge \frac{1}{2}\lambda_{\min}(K_\varepsilon(\theta(0)))$ We will add this in the revised proof.

• Item 5: A large learning rate won’t solve the slow convergence issue for small $\varepsilon$, and may lead to instability [3].

• Item 6: Thank you—we will update the notation.

• Item 7: For small $\varepsilon$, $K_\varepsilon$ may be large but remains finite. Using small $\Delta \varepsilon$ ensures stable error control.

• Items 8–9: The index should be lowercase $n$, and the correct reference is [4].

• 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.

---

Comments on the Weaknesses:
We 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.

Experimental Comparison:
Based 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 show that our homotopy-based method consistently achieves better accuracy.

---

References:

[1] Yang et al., NeurIPS, 2023.
[2] E et al., Constructive Approximation, 2022.
[3] Sutskever et al., ICML, 2013.
[4] Atkinson et al., Wiley, 2009.
[5] Krishnapriyan et al., NeurIPS, 2021.
[6] Wang et al., Comput. Methods Appl. Mech. Eng., 2021.

Thank you for your valuable suggestions on our paper.

First, 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.

Below, we provide our responses to the questions and concerns you have raised.

• 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.

  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.

  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.

  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.

• 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 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.

  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.

  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.

• 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.

Finally, we would like to thank you again for your insightful comments and questions.

Thank you for your careful reading and valuable suggestions.

First, 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.

Below, we address the questions and concerns you have raised.

• Comparision to Related Work

  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).

  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].

  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.

  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.

  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.

  Finally, we have added additional experimental comparisons, which further support the advantages of homotopy-based training and Results show that our homotopy-based method consistently achieves better accuracy.

• Theoretical Clarifications

  We appreciate the reviewer’s comments regarding theoretical assumptions and will incorporate further details in the revised version.

1. 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.

2. 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.

3. 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.

• Other Questions

• Based on Example 5.1, we added comparisons with other methods. Results in Table 3 and Figures 1–3 show that our homotopy-based method consistently achieves better accuracy.

• The training time and inference time for the numerical experiments can be found in Table 1 and 2.

• 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.

References

[1] Krishnapriyan et al., NeurIPS, 2021.
[2] Wang et al., Comput. Methods Appl. Mech. Eng., 2021.
[3] Zhang et al., J. Comput. Phys., 2024.
[4] Li et al., Comput. Methods Appl. Mech. Eng., 2023.

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.

In response to the concerns you raised, we provide the following answers:

• 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.

  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), demonstrating significant computational speedups.

  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.

• 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.

  For the second inequality (Eq.(20)), the Lidskii–Mirsky–Wielandt theorem provides the following full inequality:

  $$\lambda_{\text{min}}(D_\varepsilon^{T} D_\varepsilon)\lambda_{\text{min}}(S^{T} S) \le \lambda_{\text{min}}(K_\varepsilon) \le \lambda_{\text{max}}(D_\varepsilon^{T} D_\varepsilon)\lambda_{\text{min}}(S^{T} S)$$

  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.

• 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.

Finally, we thank the reviewer for the constructive feedback and helpful suggestions, which have greatly improved the clarity and completeness of our paper.

Reference:

• [1]Yang et al., NeurIPS, 2023.
• [2]E et al., Constructive Approximation, 2022.
• [3]Wight et al., Commun. Comput. Phys., 2021.
• [4]Zhang et al., J. Comput. Phys., 2024.
• [5]Li et al., Comput. Methods Appl. Mech. Eng., 2023.