Controlling Statistical, Discretization, and Truncation Errors in Learning Fourier Linear Operators

We are glad that the reviewer enjoyed reading our article. We address the reviewer's concern below.

1. We acknowledge that, in practice, rates faster than $1/\sqrt{n}$, and even faster than $1/n$, are often observed. This discrepancy likely occurs from the careful selection of training samples in practice, which deviates from the i.i.d. sampling assumption commonly used in statistical models. In the discussion and future work section, we highlight that alternative learning frameworks, such as active learning, could provide a more appropriate theoretical foundation for analyzing operator learning and may offer insights into these faster convergence rates.

2. Based on our lower bound of $1/n$, we require $K \geq n^{1/2s}$ and $N \geq n^{1/2s}$ to ensure that the statistical error does not dominate the truncation and discretization errors. However, as the reviewer correctly points out, the recommendation of $K \geq n^{1/4s}$ may not be optimal. We plan to revisit and refine these recommendations as we address the gap in the upper and lower bounds of the statistical rate.

3. It is indeed possible for the input space and output space to have different smoothness levels, $s_1$ and $s_2$, with $s_1 < s_2$. However, our analysis still require $s_1 > d/2$, and the statistical and truncation errors will be dominated by the less smooth space. Specifically, the upper bounds will be $N^{-s_1}$ and $K^{-2s_1}$. Since the more smooth space is contained in the less smooth space, we chose both spaces to have the lower smoothness level in our analysis.

4. We believe that bounds of the form $N^{-(s-d/2)}$ are typically obtained when analyzing the discretization error of DFT-based interpolants under Sobolev norms. In our analysis, we quantify the error under the $L^2$-norm, resulting in a rate of $N^{-s}$ for $s > d/2$. If the reviewer could provide specific references for comparison, we would be happy to review them and offer a more detailed discussion of how our bounds compare with existing results.

We thank the reviewer for noting that our article is easy to read, rigorously written, and provides many details about related topics in functional data analysis and neural operators. We address the reviewers' concerns below.

1. We agree with the reviewer that our problem involves estimating the error of a linear functional model with an integral kernel that is diagonal in the Fourier basis (though the true model itself need not be truncated—only our estimator is). We chose to motivate this problem using FNO because our target audience is the operator learning community. We are happy to move some of the details to the Appendix in the final version if necessary.

2. We acknowledge that the approximation capabilities of FNOs and our linear model are not directly comparable, as FNOs incorporate nonlinear activation functions capable of processing higher frequencies. Our analysis, however, focuses on understanding specific error components, such as truncation and discretization, which are fundamental to both linear and nonlinear models. While FNOs may effectively handle high frequencies in practice, it is not yet clear whether their discretization error can be controlled for such functions. The only available discretization error analysis for the forward pass through FNO (see https://arxiv.org/pdf/2405.02221) is in the regime $s>d/2$, where functions are dominated by low frequencies. As a result, any approximation advantage of FNOs may be offset by non-vanishing discretization errors when dealing with non-smooth functions. Since the lower bound on the discretization error for our model also applies to FNO, our model can be a valuable framework for exploring FNO’s limitations in this regard.

3. We understand the reviewer's concern regarding the simplicity of our setup. However, this was partly intentional. More importantly, our contributions extend beyond standard results by addressing the novel problem of multiresolution generalization, where operators trained on lower resolutions generalize effectively to higher resolutions as observed in practice. In particular, we establish the generalization error of operators trained on lower-resolution grids ($N^d$) that are evaluated at full resolution ($N \to \infty$). Since this is a novel problem, we do not believe that results addressing it currently exist in the literature. Additionally, while classical results often depend on the effective dimension $K^d$, our analysis establishes a $1/\sqrt{n}$ rate that is independent of $K$. This distinction is significant because such $K^d$-dependence can lead to a statistical sample complexity that is worse than the Monte Carlo rate when the truncation error is made arbitrarily small by taking $K \to \infty$. While some existing works achieve a $1/\sqrt{n}$ rate under specific assumptions, such as decay in the operator's spectrum (e.g., bounded Hilbert-Schmidt norm), to the best of our knowledge, no results establish this rate without such assumptions. In our case, this rate is achieved by using the smoothness of the input functions. Since most FDA literature assumes only $L^2$-integrability of input functions, it is unlikely that these results are present in the existing FDA literature.

4. ``I would like the authors to extend the discussion of the differences. For example, why is it significant that data is not discretized?" This is an excellent question. We believe the key factor that conceptually distinguishes operator learning from standard function-to-function regression in FDA is the role of discretization. In FDA, the grid size is typically predetermined and not within the user’s control. Even in rare cases where measurements can be taken at specific intervals, it is often difficult to make reliable assumptions about the smoothness of the functions. In contrast, operator learning uses training data generated by PDE solvers, enabling the user to determine both the grid size and the smoothness properties of the training functions. This control provides a unique learning-theoretic framework where discretization error plays a critical role in understanding the trade-offs between grid size, smoothness, and generalization. Our work emphasizes the importance of quantifying discretization error, particularly in multiresolution settings. Importantly, this error is not merely the error of the DFT for functions in a specific function space. While the lower bound on the DFT error is likely a lower bound for the discretization error, the upper bound on the DFT error does not necessarily provide an upper bound on the discretization error of the trained operator.

We thank the reviewer for their comments.

1. We agree that our problem can be framed as regression in the Fourier domain. However, quantifying the discretization error requires analysis beyond that of standard DFT analysis. Specifically, our focus is on quantifying the generalization error of an operator trained on a grid of size $N^d$ but evaluated at full resolution ($N \to \infty$). This formalizes the concept of multiresolution generalization, a phenomenon commonly observed in practice. To address this, we first established a non-trivial result showing that the discretization error of the estimator is bounded by the uniform deviation between the empirical risk at resolution $N^d$ and the empirical risk at full resolution ($N \to \infty$) (see the beginning of Appendix D). While DFT error bounds are a component of our analysis, they do not fully determine the trained operator’s discretization error. Specifically, while the lower bound on the DFT error is likely a lower bound for the discretization error, the upper bound on the DFT error may not necessarily translate to an upper bound on the discretization error of the trained operator. To the best of our knowledge, quantifying multiresolution generalization error remains an open question, particularly for general nonlinear operators.

2. We agree with the reviewer that classical statistical bounds apply after truncation. However, these bounds depend on the effective dimension $K^d$, leading to a statistical risk of $\sqrt{\frac{K^d}{n}}$, which becomes vacuous as $K \to \infty$. Our key contribution is establishing a $\frac{1}{\sqrt{n}}$ bound independent of $K^d$, a significant improvement that allows us to analyze the problem for arbitrary $K$. Achieving this required careful Rademacher analysis, as standard techniques such as covering numbers and parameter counting fail to eliminate $K^d$-dependence.

3. The $K$-free bound is significant because $K$ in FNOs is analogous to the width in standard neural networks. While generalization bounds for FNOs currently depend on $K$, it is well-established that such bounds for standard MLP networks do not depend on width (see https://arxiv.org/pdf/1712.06541). Our work provides preliminary evidence that a similar $K$-free generalization bound is possible for FNOs, contributing to the theoretical understanding of this widely used model.

4. ``The proposed method cannot be directly extended to learning nonlinear operators because the main idea is a simple linear regression in the Fourier domain.” We agree that our method is limited to linear operators and cannot be directly extended to nonlinear cases. However, this simplification is intentional, as it allows us to analyze the core component of Fourier Neural Operators (FNOs) in depth.

5. We thank the reviewer for bringing the JMLR paper (https://jmlr.org/papers/volume25/22-0719/22-0719.pdf) and other references to our attention, which we will cite in the final version. While this JMLR paper achieves impressive results, its analysis does not directly apply to our context. Specifically, the generalization bound in Corollary 11 (where the encoder/decoder are trigonometric polynomials) is $\lesssim n^{-\frac{2}{2+d_{X}}} \log^2{n} + d_{X}^{-\frac{2s}{D}} + d_Y^{-\frac{2s}{D}}$. In particular, the bound depends on the encoding dimension $d_X$, resulting in a slower decay rate, $n^{-\frac{2}{2+d_X}}$, compared to our $n^{-1/2}$ bound. This bound is only tight when $d_X \leq 4$, but even in this regime, the term $d_X^{-2s/D}$ does not vanish, making the overall bound vacuous. For larger $d_X$, the decay rate deteriorates further and becomes arbitrarily slower than the Monte-Carlo rate.

The dependence on $d_X$ in their analysis arises because the authors work with non-parametric operator classes and rely on covering numbers. Such covering number analysis inherently ties the statistical rate to the encoding dimension and cannot establish a bound independent of $d_X$, even for linear operator classes.

In our setting, $d_X$ corresponds to $K^d$, and one of our primary contributions is decoupling $K^d$ from the statistical error through a careful Rademacher analysis. This allows us to establish a $n^{-1/2}$ bound independent of $K^d$, a distinction that is not addressed in the JMLR paper. While we appreciate the depth of their analysis, we believe the goals of these two works differ slightly.

7. We acknowledge that the immediate practical implications of our results are limited. However, our work addresses foundational questions in operator learning, such as multiresolution generalization and $K$-free bounds, which are critical for bridging gaps between theory and practice. By highlighting the conceptual similarities with the FDA and identifying unique learning-theoretic challenges, we also hope to generate interest from both the functional data analysis and the learning theory communities.

4. Yes, the practical implementation of the linear core of the FNO includes a pointwise transformation $Wv(\cdot)$. If $v$ is a $\mathbb{R}^p$-valued function, $W$ is generally a $p \times p$ matrix (or $q \times p$ if the channel dimension changes across layers). However, for mappings between scalar-valued functions, $W$ reduces to a scalar constant (denoted as $c$), and the transformation $Wv$ becomes a scaled identity operator, $c \mathbf{I}$. Since the identity operator is diagonalizable in the Fourier basis, we have $c\mathbf{I} = c\sum_{m} \varphi_m \otimes \varphi_m$. In light of Proposition 2, this operator shares a similar structure as the Fourier layer $\mathcal{F}^{-1}(\Lambda \mathcal{F}(v))$. Thus, all our analysis should easily extend with minimal adjustments for the scalar-valued case. However, as the reviewer notes, the situation becomes more involved for vector-valued functions, as the matrix $W$ mixes information across channels. Addressing these challenges for vector-valued cases is an important direction for future work.

5. We have included experiments demonstrating that our estimator achieves a vanishing error at a reasonable rate. These results are presented in the final section of the revision, located in the appendix, and highlighted in blue for emphasis.

6. We will add a notation table in the appendix to clarify the symbols and their meanings. Additionally, we will separate the three error terms into distinct subsections or bullet points to address each error clearly.

We thank the reviewer for noting that our work is an important step in the direction of certifiable and trustworthy operator learning. We address reviewers' concern below.

1. We agree with the reviewer that FNOs use lifting and projection operators (via MLPs), which we do not explicitly account for. However, as stated in Section 1.3, our goal is to “conceptually separate the paradigm of operator learning from its commonly used instantiation using neural network architectures.” Thus, we omitted lifting and projection operators from our model. Although nonlinear projection/lifting is often used in practice, Remark 2.2 of (https://arxiv.org/pdf/2107.07562) notes that FNO retains universal approximation properties even when projection/lifting operators act linearly pointwise. In our case, since both input and output functions are scalar-valued, these lifting and projection operators reduce to simple scaling by constants $r$ and $q$. Because there is no nonlinearity in our model, these constants can be absorbed into the model parameters $\lambda$, and our results remain valid. However, we acknowledge that our results may not hold under nonlinear projection or lifting operations.

2. The reviewer is correct that many practical operator learning problems involve functions in $L^2(\mathbb{R}^d, \mathbb{R})$. However, since current operator learning methods are grid-based, they are effectively constrained to handling functions on a bounded subset $\Omega \subset \mathbb{R}^d$, mostly rectangular. For bounded rectangular domains, a periodic extension of the function can be performed, and for sufficiently smooth functions, any potential discontinuities introduced by this extension are confined to the boundary of $\Omega$. Since we quantify error using the $L^2$-norm (rather than pointwise deviation or uniform norms) and the boundary is a measure-zero set, such discontinuities are unlikely to affect our theoretical results. However, Gibbs-type errors may still arise in practice when computations are performed on a discrete grid. Thus, for non-periodic functions, alternative bases such as orthonormal polynomials or wavelet bases may be more suitable than the Fourier basis, as they are better equipped to handle non-periodicity. Given that our operator of interest has the general form $\sum_{m} \lambda_m \, \varphi_m \otimes \varphi_{-m}$, we could replace $\varphi_m$ with a different orthonormal basis and extend our analysis accordingly.

3. ``I’d like to see at least preliminary results with input functions that have limited smoothness." The truncation error bound of $\lesssim K^{-2s}$ holds for any $s > 0$ and does not require $s > d/2$. As discussed at the beginning of Section 4.3, $s > 0$ is necessary, as we can establish a non-vanishing lower bound when $s = 0$. Regarding discretization error, the DFT-based method is generally ineffective when $s \leq \frac{d}{2}$, as the DFT cannot reliably approximate the true Fourier transform in this regime. Specifically, for any mode $m$, it is possible to construct a function $v$ such that the error $\left| \text{DFT}(v^N)(m) - (\mathcal{F}v)(m) \right|$ does not vanish as $N \to \infty$. Finally, for the statistical error, when $s \leq d/2$, we can establish a bound of type $\sqrt{\frac{K^d}{n}}$. Then, to achieve a desired accuracy $\varepsilon$, choosing $K$ such that $K^{-2s} \leq \varepsilon$ implies the sample size must satisfy $n \geq \varepsilon^{-2 - \frac{d}{2s}}$ to make statistical error $\leq \varepsilon$. However, this is worse than our sample complexity of $n \geq \varepsilon^{-2}$ and can become arbitrarily worse when $d/2 \gg s$. If a $\frac{1}{\sqrt{n}}$ rate for the statistical error is desired without any $K$ dependence, it can be achieved under a mild assumption on the spectrum of the operator. If we assume the spectrum sequence $\lambda$ lies in $\ell^2$ rather than $\ell^{\infty}$, the operator becomes Hilbert-Schmidt. Then, we can use results from (https://arxiv.org/abs/1901.10076) and (https://arxiv.org/abs/2309.06548) to show a $1/\sqrt{n}$ upper bound on the statistical error, even when $s = 0$. We are happy to include a discussion of these results in the paper if the reviewer would like.

We thank the reviewer for engaging with us and increasing the score. We address their questions below.

1. Our goal was to remove the complexities associated with neural networks while retaining the focus on the linear core of one of the most prominent architectures in the field. By focusing on a single basis rather than a more general one, we aimed to keep the discussion centered on the interesting learning-theoretic aspects of the problem, particularly the interplay between various error types. While the analysis of statistical and truncation errors would remain largely similar for other bases, the discretization error analysis would require a separate analysis, as a different basis would require data on the different grids. For example, the estimator based on Discrete Chebyshev Transform for Chebyshev polynomials requires access to functions over a Chebyshev grid instead of the uniform grid.

2. We understand the reviewer's concern that the boundary terms may pollute the interior for non-periodic functions when DFT-based techniques are used. So, to ensure rigor and clarity, we will refrain from making any formal claims about non-periodic functions in this work and instead leave this analysis for future work.

3. We will make sure to include these discussions in the final version.

4. Including the transformation by a constant operator $c \mathbf{I}$ introduces only a minor adjustment to our analysis. Specifically, we can express $c \mathbf{I}$ as $c \mathbf{I} = c \sum_{m \in \mathbb{Z}^d} \varphi_{m} \otimes \varphi_{m}.$ Thus, the operator of interest takes the form: $T =\sum_{m \in \mathbb{Z}^d} c \varphi_{m} \otimes \varphi_{m}+ \sum_{m \in \mathbb{Z}^d} \lambda_m \varphi_m \otimes \varphi_{-m}.$

Let us call the first term $T_1$ and the second $T_2$. For the upper bound, as our analysis relies on controlling deviations between two types of risks, the terms involving $T_1$ and $T_2$ can be separated using the triangle inequality. The same steps applied to $T_2$ in the paper can be repeated verbatim for $T_1$, yielding at most an additional factor of $c$ times the previously established error. While this modification is straightforward, implementing it would result in changes throughout the paper. Thus, we will include it as a remark in the final version of the paper.