AC-PKAN: Attention-Enhanced and Chebyshev Polynomial-Based Physics-Informed Kolmogorov–Arnold Networks

We thank the reviewer for the thoughtful comments.

Response to W1

We mildly disagree with the assertion that "The experiments and the comparisons are not challenging." The experiments we selected are also employed in other studies [1–4]. You mentioned in [5] that for the 1D-wave case, a relative $L^2$ norm of $1.7 \times 10^{-3}$ was achieved, which is of the same order as that of AC-KAN. However, this result was obtained using additional external learning strategies. In our comparisons with different baselines, we did not utilize extra techniques, opting instead for the simplest possible designs to highlight the intrinsic performance of the models through experiment. Nevertheless, to strengthen our statement, we have added an experiment on Navier–Stokes flow around a cylinder (see 2D Navier–Stokes PDE Experiment). Moreover, we will discuss the paper [5-6] in our revised paper.

[1] Zhao, Z., Ding, X., Prakash, B. A. "PINNsformer: A transformer-based framework for physics-informed neural networks." arXiv preprint arXiv:2307.11833, 2023.

[2] Wu, H., Luo, H., Ma, Y., et al. "RoPINN: Region Optimized Physics-Informed Neural Networks." arXiv preprint arXiv:2405.14369, 2024.

[3] Wang, Y., Sun, J., Bai, J., et al. "Kolmogorov Arnold Informed Neural Network: A physics-informed deep learning framework for solving PDEs based on Kolmogorov Arnold Networks." arXiv preprint arXiv:2406.11045, 2024.

[4] Hao, Z., Yao, J., Su, C., et al. "Pinnacle: A comprehensive benchmark of physics-informed neural networks for solving PDEs." arXiv preprint arXiv:2306.08827, 2023.

[5] Wang, S., Yu, X., Perdikaris, P. "When and why PINNs fail to train: A neural tangent kernel perspective." Journal of Computational Physics, 2022, 449: 110768.

[6]Wang S, Sankaran S, Wang H, et al. An expert's guide to training physics-informed neural networks[J]. arXiv preprint arXiv:2308.08468, 2023.

Response to W2

We appreciate your constructive feedback and agree with your suggestion. We will revise our manuscript to streamline the method introduction, incorporate the newly conducted Navier–Stokes flow around a cylinder experiment into the main text, and provide a more in-depth discussion of theoretical insights and experimental results.

Thanks for your patient and timely reply!

Response to 1

These two NTK methods are consistent. In fact, to establish an effective comparison, our experimental setup in this section is identical to that of Table 2 in [1]. Furthermore, our results align with the experimental findings presented in Table 2 of [1]. For detailed information, please refer to our code (https://anonymous.4open.science/r/1d_wave_pinn_ntk-188A). However, it is a surprise to find that even slight modifications to the PDE parameters caused significant changes in the performance of PINN when using the NTK method. Specifically, the PDE's rRMSE between β = 4 and β = 3 differed by two orders of magnitude. The underlying reasons for this phenomenon warrant further investigation, as it suggests that the NTK method may still has certain limitations. Enhancing the robustness of the NTK method under varying PDE parameter conditions is a potential avenue for future research, although it is beyond the scope of our current study.

Response to 2

First, although our AC-PKAN’s improvement in 2D NS eq may appear limited, it is nevertheless a meaningful advancement. In our study, the primary aim is not to demonstrate that AC-PKAN thoroughly surpasses other PINN models, but rather to offer novel perspectives and alternatives through the application of KAN while achieving desirable accuracy in solving PDE systems. Furthermore, AC-PKAN represents a further—albeit still preliminary—exploratory endeavor building upon the relatively nascent development of KAN. While KAN still has performance aspects awaiting enhancement, it has been widely recognized as a promising new direction for PINNs [7].

To provide details on the two-dimensional Navier–Stokes equations, we set the parameters $\lambda_1 = 1$ and $\lambda_2 = 0.01$. Since the system lacks an explicit analytical solution, we utilize the simulated solution provided in [8]. We focus on the prototypical problem of incompressible flow past a circular cylinder—a scenario known to exhibit rich dynamic behavior and transitions across different regimes of the Reynolds number, defined as $\text{Re} = u_\infty D / \nu$. By assuming a dimensionless free-stream velocity $u_\infty = 1$, a cylinder diameter $D = 1$, and a kinematic viscosity $\nu = 0.01$, the system exhibits a periodic steady-state behavior characterized by an asymmetrical vortex shedding pattern in the cylinder wake, commonly known as the Kármán vortex street. For more comprehensive details about this problem, please refer to Section 4.1.1, "Example (Navier–Stokes equation)," in [8].

[7] Toscano J D, Oommen V, Varghese A J, et al. From PINNs to PIKANs: Recent Advances in Physics-Informed Machine Learning[J]. arXiv preprint arXiv:2410.13228, 2024.

[8] Raissi, Maziar, Paris Perdikaris, and George E. Karniadakis. "Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations." Journal of Computational physics 378 (2019): 686-707.

We hope this message finds you well.

As the deadline for final evaluations is tomorrow, we would greatly appreciate any insights or suggestions you might have to our previous response. We would also be grateful for any consideration towards improving our current score based on potential revisions.

Thank you very much for your time and understanding.

We thank the reviewer for the comments.

Response to W1

As shown in Ablation Study 2, our ablation experiments demonstrate that AC-PKAN leverages our proposed RGA more effectively than other frameworks.

Response to W2

FNO-related methods and PINNs belong to entirely different categories of research. FNO falls under another major direction of AI4PDEs—operator learning, which aims to learn a mapping operator rather than solving a specific partial differential equation.

FNO is a purely data-driven deep learning model without incorporation of physical constraints, whereas PINNs are strongly physics-constrained models with minimal reliance on data. Their learning objectives and implementation approaches are fundamentally different. As a result, it is standard practice[1-5] not to compare FNO-related methods with PINNs in PINN-related studies, nor to use PINNs as baselines in operator learning research.

[1] Zhao, Z., Ding, X., Prakash, B. A. "PINNsformer: A transformer-based framework for physics-informed neural networks." arXiv preprint arXiv:2307.11833, 2023.

[2] Wu, H., Luo, H., Ma, Y., et al. "RoPINN: Region Optimized Physics-Informed Neural Networks." arXiv preprint arXiv:2405.14369, 2024.

[3] Baez A, Zhang W, Ma Z, et al. Guaranteeing Conservation Laws with Projection in Physics-Informed Neural Networks[J]. arXiv preprint arXiv:2410.17445, 2024.

[4] Bonev B, Kurth T, Hundt C, et al. Spherical fourier neural operators: Learning stable dynamics on the sphere[C]//International conference on machine learning. PMLR, 2023: 2806-2823.

[5] Tran A, Mathews A, Xie L, et al. Factorized fourier neural operators[J]. arXiv preprint arXiv:2111.13802, 2021.

Response to Q1

No, it isn’t. Given that the dataset is extremely small—if it exists at all—it is unnecessary to train it using different batches.

PINNs are strongly physics-constrained models that embed PDE constraints directly into the loss function for optimization, requiring minimal or no data. This characteristic represents a fundamental distinction from operator learning. For instance, in our five PDE-solving experiments, no additional data was used except for the Navier–Stokes equation, which required approximately 24.08 MB of additional training data. In all other cases, no external data constraints were applied. Even when data is used, the dataset size is extremely small, and the entire dataset is utilized in every training iteration of AC-PKAN, eliminating the need to divide it into different batches.

Dear Reviewer i8XJ,

We appreciate your time for reviewing, and we really want to have a further discussion with you to see if our response solves the concerns. We have addressed all the thoughtful questions raised by the reviewer (eg, technical descriptions and contributions), and we hope that our work’s impact and results are better highlighted with our responses. It would be great if the reviewer can kindly check our responses and provide feedback with further questions/concerns (if any). We would be more than happy to address them.

Thank you!

Best wishes,

Authors

Thank you for your thoughtful reply!

Response 1

Our RGA mechanism is composed of Gradient-Related Attention (GRA) and Residual-Based Attention (RBA).

GRA is employed to alleviate the stiffness problem in the gradient flow dynamics of the PDE loss term in PINNs, as discussed in detail in [6]. Operations within the Cheby1KAN layer exacerbate this stiffness due to:

• Higher-Order Polynomials: Computing PDE residuals requires high-order derivatives, which necessitates higher-degree Chebyshev polynomial terms, leading to large coefficients and derivative magnitudes. This increases the Hessian's maximum eigenvalue, intensifying gradient flow stiffness. Moreover, higher-order polynomial terms are particularly prone to generating rapidly fluctuating numerical values, especially as the input $x$ approaches the boundaries of the interval. In such scenarios, the amplified contributions of these terms lead to instability in parameter updates along the affected directions.

• Nonlinear Operations: The nonlinear functions $\cos$ and $\arccos$ have gradients can vary significantly in certain regions. The gradient of $\cos(x)$ is $-\sin(x)$, which approaches zero as x→π/2 (gradient vanishing), while the gradient of $\arccos(x)$ is $-\frac{1}{\sqrt{1 - x^2}}$, tending to infinity as x→±1 (gradient explosion).

• Increased Nonlinearity with Higher Degrees: Higher-degree Chebyshev polynomials increase network nonlinearity, leading to more imbalanced loss gradients and exacerbated dynamic stiffness.

RBA, a lightweight dynamic weighting method based on pointwise residuals, is well-suited for computations within the Cheby1KAN layer. While Cheby1KAN captures complex physical features and excels in regions with rapid changes, its complexity can hinder convergence in certain areas. RBA enhances fitting capability by amplifying weights where residuals are high, dynamically balancing residuals across regions and improving training stability. This complements Cheby1KAN's strengths in handling strong nonlinearity and complex distributions, such as high-frequency and singular solutions. Furthermore, RBA's weighting mechanism forms a synchronously adjusted feedback loop with Cheby1KAN, mitigating numerical issues in optimization and enhancing global convergence efficiency.

A rigorous mathematical derivation of the connection between GRA and Cheby1KAN is a future research direction but beyond the scope of our current study.

[6] Wang S, Teng Y, Perdikaris P. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, 2021, 43(5): A3055–A3081.

Response 2

Firstly, in the four PDE experiments presented in our original manuscript, we did not utilize any external datasets; the data were self-generated. Specifically, we sampled points within the plane and computed the corresponding physical constraint PDE losses, as well as initial or boundary condition losses. For detailed information, please refer to lines 978–994 in the original manuscript(lines 1100-1134 in the revised pdf). We would like to further supplement this section as follows:

Training and Test Sets: The training and test sets each consist of two distinct groups containing $101 \times 101 = 10,201$ collocation points. These points are generated using the data generation method described earlier.

Secondly, for the 2D Navier–Stokes equation experiment added during the rebuttal phase, the dataset used are explained as follows:

Training and Test Sets: From the total dataset of 1,000,000 data points ($N \times T = 5,000 \times 200$), we randomly selected 2,500 samples for training, which include coordinate positions, time, and corresponding velocity and pressure components. The test set consists of all spatial data at the 100th time step.

We thank the reviewer for the thoughtful comments. Your suggestions are very important for improving our paper.

Response to W1

We present the actual running time data(hours: minutes: seconds) for all five PDE experiments in the paper, including the Navier–Stokes equation that we additionally conducted.

We trained the models until convergence but did not exceed 50,000 epochs. As shown, AC-PKAN still has certain advantages among the KAN model variants, although the running times of all the KAN variants are relatively long. This is primarily because the KAN model is relatively new and in its primitive stages; although it is theoretically novel, its engineering implementation remains coarse and lacks deeper refinements. Moreover, while traditional neural networks benefit from relatively mature optimizers such as Adam and L-BFGS, optimization schemes specifically tailored for KAN have not yet been explored. We believe the performance of AC-PKAN will be further enhanced as the overall optimization performance of KAN variants improves.

Response to W2

Please refer to Ablation Study 2

Response to W3

We reframe the statement of Proposition 1 as:

Let $\mathcal{N}$ be an AC-PKAN model with $L$ layers ($L \geq 2$) and infinite width. Then, the output $y = \mathcal{N}(x)$ has non-zero derivatives of any finite-order with respect to the input $x$.

Response to W4 & Q2

Please refer to 2D Navier–Stokes PDE Experiment. Our AC-PKAN model continues to demonstrate strong performance.

Response to W5

Please refer to Ablation Study 1. Rather than increasing the scaling factor, we consider it more reasonable to multiply the IC/BC loss term ,with its large coefficient $\lambda^{\text{GRA}}$, by a very small scaling factor. This approach ensures that the overall product of IC/BC loss terms’ coefficients approximates the magnitude of overall coefficient of the PDE residual loss term whose scaling factor is set to 1.

The motivation for the logarithmic transformation originates from the Bode plot in control engineering, which uses a logarithmic frequency axis while directly labeling actual frequency values. This technique not only compresses a wide frequency range but also linearizes the system's gain and phase characteristics on a logarithmic scale, thereby alleviating imbalances caused by significant differences in data scales.

Response to Q1

In fact, all our experiments use the AdamW optimizer and we clearly stated in the appendix. Step 7 in Algorithm 1 refers to a generic model parameter update step and does not mandate the use of SGD. The focus of our algorithm lies in loss weighting rather than parameter updates. Actually, it is common practice in papers to use gradient descent in pseudo-code as a simplified representation of the parameter update step without specifying a particular optimizer.

Response to Q3

It is true that an MLP with sinusoidal activations guarantees non-zero derivatives in its output. However, this advantage specifically addresses improvements to KAN variants rather than MLP-based models. The fitting capability of KAN models relies on polynomial functions with learnable parameters. To ensure non-zero derivatives in the output, the original B-spline KAN includes an additional nonlinear bias function $b(x)$. In contrast, some other KAN variants, such as Cheby1KAN, rely solely on polynomial bases, which inevitably result in zero derivatives when the order of differentiation exceeds the polynomial degree. Furthermore, Proposition 1 was primarily proposed to provide a theoretical guarantee for AC-PKAN's ability to solve any finite-order PDE, analogous to how the universal approximation theorem theoretically establishes the universal fitting capability of neural networks.

Response to Q4

Other PINN methods do not require similar linear projections because they are inherently MLP-based models. In AC-PKAN, the linear layer is designed to achieve a hybridization of KAN and MLP architectures. Its role as an initial projection was inspired by the Spatio-Temporal Mixer linear layer in the PINNsformer model, which enhances spatiotemporal aggregation.

Dear Reviewer mYon,

We appreciate your time for reviewing, and we really want to have a further discussion with you to see if our response solves the concerns. We have addressed all the thoughtful questions raised by the reviewer (eg, technical descriptions and contributions), and we hope that our work’s impact and results are better highlighted with our responses. It would be great if the reviewer can kindly check our responses and provide feedback with further questions/concerns (if any). We would be more than happy to address them.

Thank you!

Best wishes,

Authors

We thank the reviewer for the thoughtful comments and appreciation for our work.

Response to Q1

The wavelet activation function, introduced in [A], is defined as $\text{Waveact}(x) = \omega_1 \sin(x) + \omega_2 \cos(x)$, with learnable parameters $\omega_1$ and $\omega_2$. This function combines the benefits of the Fourier feature embedding method from [B] and the sine activation function from [C].

In [B], input coordinates $x$ are mapped to high-frequency signals as:

$$\gamma(x) = \begin{bmatrix} \cos(Bx) \\ \sin(Bx) \end{bmatrix},$$

where $B \in \mathbb{R}^{m \times d}$ is a matrix with elements sampled from a Gaussian distribution $N(0, \sigma^2)$, and $\sigma$ controls the frequency range.

In [C], the First-Layer Sine (FLS) method applies $\sin(x)$ as the activation function in the first layer, improving gradient distribution in early training.

The wavelet activation function, employed in the initial encoders $U$ and $V$, retains the advantages of FLS for proper gradient initialization. Intuitively, compared to the pure sine activation function, the wavelet activation introduces additional phase and magnitude information. The arbitrary and learnable combination of phase and magnitude offers greater representational capacity than the fixed phase and magnitude of the pure sine activation function. Moreover, through learnable parameters, it enables adaptive Fourier feature embedding, enhancing the network's ability to capture periodic features and reduce spectral bias.

[A] Zhao Z, Ding X, Prakash B A. Pinnsformer: A transformer-based framework for physics-informed neural networks. arXiv preprint arXiv:2307.11833, 2023.
[B] Wang S, Sankaran S, Wang H, et al. An expert's guide to training physics-informed neural networks. arXiv preprint arXiv:2308.08468, 2023.
[C] Wong J C, Ooi C C, Gupta A, et al. Learning in sinusoidal spaces with physics-informed neural networks. IEEE Transactions on Artificial Intelligence, 2022, 5(3): 985-1000.

Response to Q2

1. Performance of Chebyshev Type 1 Polynomials in High-Frequency Scenarios

Chebyshev polynomials of the first kind, $T_n(x) = \cos(n \arccos(x))$, focus their spectrum on high frequencies, with the frequency increasing linearly with polynomial order $n$. This makes them suitable for capturing high-frequency oscillations, as their high-frequency components decay slowly. The even distribution of extrema further aids in capturing rapid variations, which is beneficial for high-frequency features.

B-splines, being piecewise polynomials, have a rapidly decaying spectrum, limiting their ability to capture high-frequency features.

2. Global Support and Orthogonality of Chebyshev Polynomials vs. Local Support of B-Splines

Chebyshev polynomials have both global support and global orthogonality over the interval $[-1, 1]$. The value of a Chebyshev polynomial at any point depends on all points within the interval, making them highly effective at capturing global features and high-frequency components. They satisfy the orthogonality relation:

$$\int_{-1}^1 \frac{T_m(x) \\, T_n(x)}{\sqrt{1 - x^2}} \\, dx = 0,\quad \text{for } m \ne n.$$

This orthogonality allows Chebyshev polynomials to achieve minimax approximation, minimizing the maximum error over the interval.

In contrast, B-splines have local support; each basis function is nonzero only within a specific subinterval. This local nature limits their ability to capture global high-frequency features. Additionally, B-splines lack global orthogonality, reducing their efficiency in approximating functions.

3. Memory Efficiency of Chebyshev Polynomials Compared to B-Splines

B-spline-based KANs require substantial memory due to the storage of grids and coefficient matrices that scale cubically with grid size and spline order. They store grids of size $(in\\_features , grid\\_size + 2 \times spline\\_order + 1)$ and coefficient matrices of size $(out\\_features , in\\_features , grid\\_size + spline\\_order)$. Since B-splines are piecewise polynomials, each segment requires maintaining basis function values and performing high-order interpolation within its support interval. This involves generating polynomial bases, solving linear systems (e.g., using torch.linalg.lstsq), and executing recursive updates, resulting in high computational and storage demands.

In contrast, Chebyshev polynomials are globally defined and require only a coefficient matrix of size $(input\\_dim , output\\_dim , degree + 1)$, eliminating terms directly related to grid size and spline order. The memory complexity grows linearly with the degree. Chebyshev polynomials eliminate the need for grid storage and do not require solving linear systems, interpolation, or recursive updates of piecewise basis functions, which significantly reduce computational and storage requirements.