We thank the reviewer for the review and the positive feedback. Should the reviewer require any further clarification, we would be delighted to provide it.

We sincerely thank the reviewer for his time and thoughtful comments and appreciate the global positive feedback. First, we would like to make it clear that the main objective of our contribution is to provide a novel theoretical framework that allows the development of founded and effective models, which we will clarify in the abstract and introduction.

Questions

Q1: Our theoretical framework formulates BFNO layers as solutions to Bregman regularized optimization problems. When all weights are zero, they naturally reduce exactly to the identity mapping (Remark 3.6), creating an implicit regularization effect distinct from standard residual connections.

We created a novel architecture, ResFNO, precisely to perform an ablation study with ResNet-style connection, confirming this distinction. While ResFNO still degrades with increased depth, BFNO maintains or improves performance. This difference manifests in Figure 5, where BFNO exhibits a distinctive Laplace-like weight distribution sharply peaked at zero, contrasting with the Gaussian distributions in other models.

From an ODE perspective (see Add. insights), the difference becomes clearer: BFNO applies its linear operator outside the activation function $\frac{dz(t)}{dt} = K(t)\sigma(z(t))$, while residual networks place it inside $\frac{dv(t)}{dt} = \sigma(K(t)v(t))$, creating different gradient flow dynamics.

Overall, BNFO's improvement can be analyzed as combination of these two elements: an implicit regularization and a different gradient flow. While a rigorous layer-wise analysis establishing a causal link between weight distribution and improved generalization would be valuable, such analysis remains an open challenge in the neural operator literature.

Q2: The primary goal of our experiments was to validate our theoretical framework and demonstrate that the properties derived from our Bregman proximal viewpoint translate into measurable performance advantages. While focusing on the theoretical foundations, we ensured rigorous experimental comparisons to establish practical relevance (even by considering WNO models in Appendix D.2).

To ensure a fair comparison, we conducted thorough validation through multiple runs with four fixed seeds and cross-validation of learning rates for all models. To assess statistical significance, we performed non-parametric one-tailed Wilcoxon signed-rank tests. Compared to FNO, BFNO is superior on 5 datasets of Figure 6 with a p-value of 0.0625, with p=0.125 for the last one 2D NS $(10^{-3})$, which shows that improvements over FNO are consistent across diverse PDE types. Similarly, BFNO achieves the same level of superiority over ResFNO on 4 datasets (Burgers, Darcy, NS $10^{-4}$ and $10^{-8}$), and on 3 datasets over F-FNO (Burgers, NS $10^{-4}$ and $10^{-8}$).

Q3: From a theoretical perspective, expressing activation functions as proximity operators imposes constraints: Non-monotonic functions like GeLU are incompatible, as they cannot be represented as proxs (page 5). For Classical Neural Operators (Sec 3.2), activation functions must be monotonic, making ReLU suitable despite its non-smoothness, as it can be expressed as a proximity operator $prox_g$ of $g=\imath_{[0,+\infty[}$. In contrast, Bregman Neural Operators (Sec 3.3) impose strict monotonicity, theoretically excluding ReLU.

However, our implementation provides significant practical flexibility. As noted in Remark 3.7, our architecture doesn't require computing inverse activations when composing layers since they cancel out with previous activations (enabled by adding an activation before the first operator layer). Therefore, any activation function can be used in practice.

Results below, completing Appendix D.4, confirm this flexibility, showing BFNO works effectively with both ReLU and GeLU, performing comparably to SoftPlus, which serves as a common smooth approximation to ReLU in the literature when differentiability is needed.

As suggested by the results, extending our theoretical framework to non-invertible and non-monotonic activation functions is an interesting research direction.

Remarks

On computational costs: BFNO shares the same parameter count and memory footprint as FNO by design. Therefore, the training time per epoch is in the same order across all tasks. This will be clarified in the final version.

On universal approximation: We indeed acknowledge this limitation in our conclusion. The theorem serves as a first theoretical result for our framework and represents, to the best of our knowledge, one of the first universal approximation results for neural operators with Kovachki et al. (2023).

As for the other remarks, we tried to address them in the answers above while staying within the allowed limit.

We thank the reviewer for his positive feedback. We provide below our answers to the two questions mentioned in the review.

1. Normally, in the context of optimization, the Bregman divergence needs to be a strongly convex function so that it represents a notion of distance. However, here it seems it only requires it to be strictly convex. I am wondering why strictly convex functions suffices?

We thank the reviewer for the thoughtful remark. Indeed, strict convexity of the Legendre function $\psi$ is sufficient to define a Bregman divergence $D_\psi(x, y)$, as it ensures that $D_\psi(x, y) \geq 0$ and that $D_\psi(x, y) = 0$ if and only if $x = y$. These properties hold as long as $\psi$ is strictly convex and differentiable on the interior of its domain — which is the standard requirement in the general theory of Bregman divergences.

That being said, we would like to emphasize that in our setting — as shown in Table 1 — all Legendre functions considered are actually 1-strongly convex on their respective domains. This ensures that the associated Bregman divergences are lower bounded by a quadratic form, i.e.,

which provides metric-like behavior. We have clarified this point in the revised version.

2. I am wondering what is the intuition of choosing the neural operator to be a solution of an optimization problem? In the context of optimization, proximal operators are often used for non-smooth problems since it has an explicit smoothing phenomenon. Is this one of the considerations for choosing a Bregman proximal operator?

We appreciate the reviewer's insightful question.

While it is true that proximal operators are classically used to handle non-smooth terms, especially in composite minimization settings, they are not restricted to non-smooth functions. In fact, proximal maps of smooth (even strongly convex) functions are well-defined and have been studied in the context of regularization and smoothing.

In our setting, the motivation for casting the neural operator as the solution to an optimization problem is not solely driven by the smoothing property of proximal operators-although this can be beneficial-but more importantly by their connection to activation functions.

It follows that, from our perspective, each layer performs a structured update—balancing proximity to the previous state with the minimization of an energy functional. This interpretation introduces a variational inductive bias into the architecture, which can encode prior knowledge and improve interpretability and performance. Using Bregman divergences further allows for this trade-off to be geometry-aware.

We hope this novel viewpoint will be leveraged by practitioners to encode structured prior knowledge into the architecture of neural operators through principled, optimization-based layers.